text:
"""
GTH algorithm
"""

"""
gth!(Q::Matrix{Tv})
gth(Q::Matrix{Tv})

Compute the stationary vector for the CTMC with kernel Q with GTH algorithm.
Both of gth! and gth return a vector.
In gth!, the matrix Q is used as a workspace.
Note that Q does not have any absorbing state.
"""

function gth!(Q::Matrix{Tv})::Vector{Tv} where {Tv}
    @inbounds begin
        m, n = size(Q)
        @assert m == n
        for l = n:-1:2
            tmp::Tv = 0
            for u = 1:l-1
                tmp += Q[l,u]
            end
            for j = 1:l-1
                for i = 1:l-1
                    if i != j
                        Q[i,j] += Q[l,j] * Q[i,l] / tmp
                    end
                end
            end
            for i = 1:l-1
                Q[i,l] /= tmp
            end
            for i = 1:l-1
                Q[l,i] = 0
            end
            Q[l,l] = -1
        end
        x = Vector{Tv}(undef, n)
        x[1] = 1.0
        for l = 2:n
            x[l] = 0.0
            for i = 1:l-1
                x[l] += x[i] * Q[i,l]
            end
        end
        x /= sum(x)
    end
end

function gth(Q::Matrix{Tv}) where {Tv}
    gth!(copy(Q))
end

"""
gth!(Q::Matrix{Tv}, index::Vector{Ti})
gth(Q::Matrix{Tv}, index::Vector{Ti})

Compute the stationary vector for the CTMC with kernel Q with GTH algorithm.
Both of gth! and gth return a vector.
The CTMC states are permutated with an index vector.
In gth!, the matrix Q is used as a workspace.
Note that Q does not have any absorbing state.
"""

function gth!(Q::Matrix{Tv}, index::Vector{Ti})::Vector{Tv} where {Tv,Ti}
    m, n = size(Q)
    @assert m == n
    for l = n:-1:2
        tmp::Tv = 0
        for u = 1:l-1
            tmp += Q[index[l],index[u]]
        end
        for j = 1:l-1
            for i = 1:l-1
                if i != j
                    Q[index[i],index[j]] += Q[index[l],index[j]] * Q[index[i],index[l]] / tmp
                end
            end
        end
        for i = 1:l-1
            Q[index[i],index[l]] /= tmp
        end
        for i = 1:l-1
            Q[index[l],index[i]] = 0
        end
        Q[index[l],index[l]] = -1
    end
    x = Vector{Tv}(undef, n)
    x[index[1]] = 1.0
    for l = 2:n
        x[index[l]] = 0.0
        for i = 1:l-1
            x[index[l]] += x[index[i]] * Q[index[i],index[l]]
        end
    end
    x /= sum(x)
end

function gth(Q::Matrix{Tv}, index::Vector{Ti}) where {Tv, Ti}
    gth!(copy(Q), index)
end


text:
module VectorWrapper

using JavaCall

include("LeastsquaresWrapper/LeastsquaresWrapper.jl")

end

ext: jl


text:
function nfkbOde!(delta, x, v, t)
    ############################################################
    # ODE file for the whole model.
    ############################################################
    # We ran the ode15s solver with the default options [an empty []]
    #  See the ode15s help file for a list of options you can use

    # Load Previous Concentrations
    # [a,b] = size(x)
    # delta = zeros(a,b)

    # nf-kb part
    IkBa           = x[1]
    IkBan          = x[2]
    IkBaNFkB       = x[3]
    IkBaNFkBn      = x[4]
    IkBat          = x[5]
    IkBb           = x[6]
    IkBbn          = x[7]
    IkBbNFkB       = x[8]
    IkBbNFkBn      = x[9]
    IkBbt          = x[10]
    IkBe           = x[11]
    IkBen          = x[12]
    IkBeNFkB       = x[13]
    IkBeNFkBn      = x[14]
    IkBet          = x[15]
    NFkB           = x[16]
    NFkBn          = x[17]

    # TLR4&IKK part
    LPSo           = x[18]; # input LPS
    LPS            = x[19];
    LPSen          = x[20];
    TLR4           = x[21];
    TLR4en         = x[22];
    TLR4LPS        = x[23]
    TLR4LPSen      = x[24];
    MyD88          = x[25]
    MyD88s         = x[26]
    TRIF           = x[27];
    TRIFs          = x[28];
    TRAF6          = x[29]
    TRAF6s         = x[30]
    IKKK_off       = x[31];
    IKKK           = x[32];
    IKK_off        = x[33];
    IKK            = x[34];
    IKK_i          = x[35];

    # TNF part
    TNFnas        = x[36]
    TNFmRNA       = x[37]
    TNFpro        = x[38]
    TNF           = x[39]

    # TNFR part
    tnfrm          = x[40]
    TNFR           = x[41]
    TNFRtnf        = x[42]
    C1             = x[43]
    C1_off         = x[44]
    C1tnf          = x[45]
    C1tnf_off      = x[46]
    TTR            = x[47]
    # A20
    a20            = x[48]
    a20t           = x[49]

    ## Set iteration changes [delta] to zero
    delta_IkBa        = 0
    delta_IkBan       = 0
    delta_IkBaNFkB    = 0
    delta_IkBaNFkBn   = 0
    delta_IkBat       = 0
    delta_IkBb        = 0
    delta_IkBbn       = 0
    delta_IkBbNFkB    = 0
    delta_IkBbNFkBn   = 0
    delta_IkBbt       = 0
    delta_IkBe        = 0
    delta_IkBen       = 0
    delta_IkBeNFkB    = 0
    delta_IkBeNFkBn   = 0
    delta_IkBet       = 0
    delta_NFkB        = 0
    delta_NFkBn       = 0
    # TLR4 module
    delta_LPS         = 0
    delta_TLR4        = 0
    delta_TLR4LPS     = 0
    delta_MyD88       = 0
    delta_MyD88s      = 0
    delta_TRIF        = 0
    delta_TRIFs       = 0
    delta_LPSen       = 0
    delta_TLR4LPSen   = 0
    # IKK module
    delta_IKKK        = 0
    delta_IKKK_off    = 0
    delta_IKK         = 0
    delta_IKK_off     = 0
    delta_IKK_i       = 0
    # TLR4 continue()
    delta_TRAF6       = 0
    delta_TRAF6s      = 0
    delta_TLR4en      = 0

    delta_LPSo           = 0

    # TNF component
    delta_TNFnas  = 0
    delta_TNFmRNA = 0
    delta_TNFpro  = 0
    delta_TNF  = 0

    # TNFR part
    delta_tnfrm       = 0
    delta_TNFR        = 0
    delta_TNFRtnf     = 0
    delta_C1          = 0
    delta_C1_off      = 0
    delta_C1tnf       = 0
    delta_C1tnf_off   = 0
    delta_TTR         = 0

    # A20 NEW
    delta_a20         = 0
    delta_a20t        = 0

    ## set IKK flux

    IKK_flux = IKK / (v.IKK_TOTAL); # ikk modification

    # TNF mRNA stabilization
    fa = (TRIFs + v.TP[10]) / (TRIFs + v.TP[11])
    fi = v.TP[12] / (TRIFs + v.TP[12])

    kdegm = v.TP[5] * fi
    kpr = v.TP[6] * fa
    ktl = v.TP[7] * fa
    ksec = v.TP[9] * fa

    # golgi plug effect
    kdeg_tnf = v.TP[8]
    ktl_a = v.NP[16]
    ktl_b = v.NP[17]
    ktl_e = v.NP[18]

    if (v.flag_noTnfFeedback) & (t > 0)
        golgi = t / (t + v.TP[14])
        ksec = ksec * (1 - golgi)
        ktnfr = v.AP[2] * (1 - golgi)
        ktlr4 = v.IP[9] * (1 - golgi)
        kshut = v.IP[2] * (1 - golgi);
        kshut_com = v.IP[14] * (1 - golgi)
    else
        ktnfr = v.AP[2]
        ktlr4 = v.IP[9]
        kshut = v.IP[2];
        kshut_com = v.IP[14]
    end

    ## Calculate Iteration Fluxes

    # -- NFkB Activation Module
    flux_rsu_a         = v.NP[1]
    flux_rsu_b         = v.NP[2]
    flux_rsu_e         = v.NP[3]
    flux_rsu_f         = v.TP[1] ; # tnf basal transcription [NEW]

    flux_rsr_an        = v.NP[4]   * (NFkBn^v.NP[7])
    flux_rsr_bn        = v.NP[5]   * (NFkBn^v.NP[8])
    flux_rsr_en        = v.NP[6]   * (NFkBn^v.NP[9])

    flux_rsr_fn        = v.TP[2]  * (NFkBn^v.TP[3]) / (NFkBn^v.TP[3] + v.TP[4]^v.TP[3]); # tnf transcription
    flux_rp_f         = kpr * TNFnas; # tnf RNA processing

    flux_rd_a          = v.NP[13]  * IkBat
    flux_rd_b          = v.NP[14]  * IkBbt
    flux_rd_e          = v.NP[15]  * IkBet
    flux_rd_f          = kdegm * TNFmRNA; # tnf mRNA degradation

    flux_ps_c_a        = ktl_a	* IkBat
    flux_ps_c_b        = ktl_b	* IkBbt
    flux_ps_c_e        = ktl_e	* IkBet
    flux_ps_c_f        = ktl	* TNFmRNA;  # tnf protein synthsis
    flux_pd_f          = kdeg_tnf    * TNFpro;  # tnf protein deg. (inside the cell)
    flux_sc_c_f        = ksec   * TNFpro;  # tnf protein secrection

    flux_in_a          = v.NP[19]	* IkBa
    flux_in_b          = v.NP[20]	* IkBb
    flux_in_e          = v.NP[21]	* IkBe
    flux_in_n          = v.NP[22]	* NFkB
    flux_ex_a          = v.NP[23]	* IkBan
    flux_ex_b          = v.NP[24]	* IkBbn
    flux_ex_e          = v.NP[25]	* IkBen
    flux_ex_n          = v.NP[26]	* NFkBn
    flux_in_2an        = v.NP[27]	* IkBaNFkB
    flux_in_2bn        = v.NP[28]	* IkBbNFkB
    flux_in_2en        = v.NP[29]	* IkBeNFkB
    flux_ex_2an        = v.NP[30]	* IkBaNFkBn
    flux_ex_2bn        = v.NP[31]	* IkBbNFkBn
    flux_ex_2en        = v.NP[32]	* IkBeNFkBn

    flux_pd_c_a        = v.NP[33]	* IkBa
    flux_pd_c_b        = v.NP[34]	* IkBb
    flux_pd_c_e        = v.NP[35]	* IkBe
    flux_pd_n_a        = v.NP[36]	* IkBan
    flux_pd_n_b        = v.NP[37]	* IkBbn
    flux_pd_n_e        = v.NP[38]	* IkBen
    flux_pd_c_2an      = v.NP[39]	* IkBaNFkB
    flux_pd_c_2bn      = v.NP[40]	* IkBbNFkB
    flux_pd_c_2en      = v.NP[41]	* IkBeNFkB
    flux_pd_n_2an      = v.NP[42]	* IkBaNFkBn
    flux_pd_n_2bn      = v.NP[43]	* IkBbNFkBn
    flux_pd_n_2en      = v.NP[44]	* IkBeNFkBn

    flux_a_c_an        = v.NP[45]  * IkBa	* NFkB
    flux_a_c_bn        = v.NP[46]  * IkBb	* NFkB
    flux_a_c_en        = v.NP[47]  * IkBe	* NFkB
    flux_a_n_an        = v.NP[48]	* IkBan	* NFkBn
    flux_a_n_bn        = v.NP[49]	* IkBbn	* NFkBn
    flux_a_n_en        = v.NP[50]  * IkBen	* NFkBn
    flux_d_c_an        = v.NP[51]	* IkBaNFkB
    flux_d_c_bn        = v.NP[52]	* IkBbNFkB
    flux_d_c_en        = v.NP[53]	* IkBeNFkB
    flux_d_n_an        = v.NP[54]	* IkBaNFkBn
    flux_d_n_bn        = v.NP[55]	* IkBbNFkBn
    flux_d_n_en        = v.NP[56]	* IkBeNFkBn

    # IKK Mediated IkB Degradation [free & bound]
    flux_ph_c_a        = v.NP[57]  * IkBa * IKK_flux
    flux_ph_c_b        = v.NP[58]	* IkBb * IKK_flux
    flux_ph_c_e        = v.NP[59]	* IkBe * IKK_flux

    flux_ph_c_an       = v.NP[60]  * IkBaNFkB * IKK_flux
    flux_ph_c_bn       = v.NP[61]	* IkBbNFkB * IKK_flux
    flux_ph_c_en       = v.NP[62]	* IkBeNFkB * IKK_flux

    # A20 Fluxes
    flux_rsu_a20       = v.NP[63]
    flux_rsr_a20       = v.NP[64]  * (NFkBn^v.NP[65])
    flux_rd_a20        = v.NP[66]  * a20t
    flux_ps_c_a20      = v.NP[67]  * NFkBn
    flux_pd_c_a20      = v.NP[68]  * a20

    if ( t .> v.NP[69])  # disables inducible A20 txn to match exp data
        flux_rsr_a20 = 0
    end

    # -- IKK Activation Module & TLR4 module
    # LPS binding to the surface()
    flux_b_LPS         = v.IP[1]   * LPSo

    # ligand binding
    flux_b_TLR4LPS     = v.IP[5]    * TLR4 * LPS
    flux_ub_TLR4LPS    = v.IP[6]    * TLR4LPS
    flux_b_TLR4LPSen   = v.IP[7]    * TLR4en * LPSen
    flux_ub_TLR4LPSen  = v.IP[8]    * TLR4LPSen

    # Activation modified on 0426
    flux_a_MyD88       = v.IP[18]    * (TLR4LPS)^v.IP[19] / ((TLR4LPS)^v.IP[19] + (v.IP[20])^v.IP[19])  * MyD88
    flux_i_MyD88       = v.IP[21]    * MyD88s
    flux_a_TRIF        = v.IP[22]    * TLR4LPSen * TRIF
    flux_i_TRIF        = v.IP[23]   * TRIFs

    # shuttling; modified on 0426
    flux_in_LPS        = kshut   * LPS
    flux_out_LPS       = v.IP[3]   * LPSen
    flux_in_TLR4       = v.IP[12]   * TLR4
    flux_out_TLR4      = v.IP[13]   * TLR4en
    flux_in_TLR4LPS    = kshut_com   * TLR4LPS
    flux_out_TLR4LPS   = v.IP[15]   * TLR4LPSen

    # generation & degradation
    flux_g_TLR4        = ktlr4
    flux_d_TLR4        = v.IP[10]   * TLR4
    flux_d_TLR4en      = v.IP[11]   * TLR4en
    flux_d_LPSen       = v.IP[4]   * LPSen
    flux_d_TLR4LPSen   = v.IP[17]   * TLR4LPSen
    # modified 0515
    flux_d_TLR4LPS     = v.IP[16]   * TLR4LPS
    # -----------------------IKK module
    # TRAF6
    flux_a_TRAF6_MyD88s = v.IP[24]   * TRAF6 * MyD88s
    flux_a_TRAF6_TRIFs = v.IP[25]   * TRAF6 * TRIFs
    flux_i_TRAF6       = v.IP[26]   * TRAF6s

    # IKKK

    flux_IKKK_on       = v.IP[28]   * IKKK_off
    flux_IKKK_on_TRAF6s = v.IP[27]   * IKKK_off * TRAF6s
    flux_IKKK_off      = v.IP[29]   * IKKK

    # IKK
    flux_IKK_on        = v.IP[31]   * IKK_off
    flux_IKK_on_IKKK   = v.IP[30]   * IKK_off * IKKK
    flux_IKK_off       = v.IP[32]   * IKK
    flux_IKK_off_i     = v.IP[33]   * IKK
    flux_IKKi          = v.IP[34]   * IKK_i

    ## TNFR part
    flux_pd_m_tnf      = v.AP[1]    * TNF ; # pd_m_tnf 45' half life of exogenous TNF
    flux_syn_tnfrm     = ktnfr;  # tnfrm synthesis [txn, tsl, localization]
    flux_pd_tnfrm      = v.AP[3]    * tnfrm;# tnfrm --> deg  -- 120' halflife
    flux_a_tnfrm       = v.AP[4]    * tnfrm;# 3tnfrm --> TNFR
    flux_d_TNFR        = v.AP[5]    * TNFR;# TNFR   --> 3tnfrm
    flux_i_TNFR        = v.AP[6]    * TNFR;# TNFR internalization -- 30' halflife
    flux_a_C1_off      = v.AP[7]    * TNFR * TTR;# TNFR + TTR --> C1_off
    flux_d_C1_off      = v.AP[8]    * C1_off;# C1_off --> TNFR + TTR
    flux_a_C1          = v.AP[9]    * C1_off;# C1_off --> C1
    flux_C1_off        = v.AP[10]   * C1;# C1     --> C1_off
    flux_C1_A20        = v.AP[11]   * C1 * a20;# C1     --> C1_off [A20 Mediated]
    flux_d_C1          = v.AP[12]   * C1; # UPDATE # C1     --> TNFR + TTR
    flux_i_C1_off      = v.AP[13]   * C1_off; # C1_off internalization
    flux_i_C1          = v.AP[14]   * C1; # C1 internalization
    flux_a_tnfrmtnf    = v.AP[15]   * tnfrm * TNF; # 3tnfrm + tnf --> TNFRtnf
    flux_a_TNFRtnf     = v.AP[16]   * TNFR * TNF; # TNFR + tnf --> TNFRtnf
    flux_d_TNFRtnf     = v.AP[17]   * TNFRtnf; # TNFRtnf   --> TNFR + tnf
    flux_i_TNFRtnf     = v.AP[18]   * TNFRtnf;  # TNFRtnf internalization -- 30' halflife

    flux_a_C1tnf_off   = v.AP[19]   * TNFRtnf * TTR; # TNFRtnf + TTR --> C1tnf_off
    flux_d_C1tnf_off   = v.AP[20]   * C1tnf_off; # C1tnf_off --> TNFRtnf + TTR
    flux_a_C1tnf       = v.AP[21]   * C1tnf_off; # C1tnf_off --> C1tnf
    flux_C1tnf_off     = v.AP[22]   * C1tnf; # C1tnf     --> C1tnf_off
    flux_C1tnf_A20     = v.AP[23]   * C1tnf * a20;  # C1tnf     --> C1tnf_off [A20 Mediated]
    flux_d_C1tnf       = v.AP[24]   * C1tnf; # updated  # C1tnf     --> TNFRtnf + TTR
    flux_i_C1tnf       = v.AP[25]   * C1tnf; # C1tnf internalization
    flux_i_C1tnf_off   = v.AP[26]   * C1tnf_off;  # C1tnf_off internalization
    flux_d_tnf_C1_off  = v.AP[27]   * C1tnf_off; # C1tnf_off --> C1_off + tnf
    flux_a_tnf_C1_off  = v.AP[28]   * C1_off * TNF; # C1_off + tnf --> C1tnf_off
    flux_d_tnf_C1      = v.AP[29]   * C1tnf; # C1tnf    --> C1 + tnf
    flux_a_tnf_C1      = v.AP[30]   * C1 * TNF; # C1 + tnf --> C1tnf
    flux_IKKK_on_C1    = v.AP[31]   * IKKK_off * C1; # IKKK_off --> IKKK [C1 mediated]?500
    flux_IKKK_on_C1tnf = v.AP[32]   * IKKK_off * C1tnf;# IKKK_off --> IKKK [C1tnf mediated]

    ## Add Fluxes to component concentrations; then Save

    # IkB Transcrv.IPtion
    delta_IkBat     = delta_IkBat + flux_rsu_a
    delta_IkBbt     = delta_IkBbt + flux_rsu_b
    delta_IkBet     = delta_IkBet + flux_rsu_e
    delta_IkBat     = delta_IkBat + flux_rsr_an
    delta_IkBbt     = delta_IkBbt + flux_rsr_bn
    delta_IkBet     = delta_IkBet + flux_rsr_en
    delta_TNFnas    = delta_TNFnas + flux_rsu_f; # NEW
    delta_TNFnas    = delta_TNFnas + flux_rsr_fn; # NEW
    delta_TNFnas    = delta_TNFnas - flux_rp_f; # NEW
    delta_TNFmRNA   = delta_TNFmRNA + flux_rp_f; # NEW

    # IkB Transcrv.IPt Degradation
    delta_IkBbt     = delta_IkBbt - flux_rd_b
    delta_IkBet     = delta_IkBet - flux_rd_e
    delta_IkBat     = delta_IkBat - flux_rd_a
    delta_TNFmRNA   = delta_TNFmRNA  - flux_rd_f;   # tnft degradation

    # IkB Translation
    delta_IkBat     = delta_IkBat - flux_ps_c_a
    delta_IkBat     = delta_IkBat + flux_ps_c_a
    delta_IkBa      = delta_IkBa  + flux_ps_c_a

    delta_IkBbt     = delta_IkBbt - flux_ps_c_b
    delta_IkBbt     = delta_IkBbt + flux_ps_c_b
    delta_IkBb      = delta_IkBb  + flux_ps_c_b

    delta_IkBet     = delta_IkBet - flux_ps_c_e
    delta_IkBet     = delta_IkBet + flux_ps_c_e
    delta_IkBe      = delta_IkBe  + flux_ps_c_e

    delta_TNFpro    = delta_TNFpro + flux_ps_c_f;

    # TNF secretion + protein degradation
    delta_TNFpro     = delta_TNFpro - flux_sc_c_f;
    delta_TNF        = delta_TNF    + flux_sc_c_f;
    delta_TNFpro     = delta_TNFpro - flux_pd_f;

    # IkB:NFkB Shuttling [Free & Bound]
    delta_IkBa      = delta_IkBa  - flux_in_a
    delta_IkBan     = delta_IkBan + flux_in_a

    delta_IkBb      = delta_IkBb  - flux_in_b
    delta_IkBbn     = delta_IkBbn + flux_in_b

    delta_IkBe      = delta_IkBe  - flux_in_e
    delta_IkBen     = delta_IkBen + flux_in_e

    delta_NFkB      = delta_NFkB  - flux_in_n
    delta_NFkBn     = delta_NFkBn + flux_in_n

    delta_IkBan     = delta_IkBan - flux_ex_a
    delta_IkBa      = delta_IkBa  + flux_ex_a

    delta_IkBbn     = delta_IkBbn - flux_ex_b
    delta_IkBb      = delta_IkBb  + flux_ex_b

    delta_IkBen     = delta_IkBen - flux_ex_e
    delta_IkBe      = delta_IkBe  + flux_ex_e

    delta_NFkBn     = delta_NFkBn - flux_ex_n
    delta_NFkB      = delta_NFkB  + flux_ex_n

    delta_IkBaNFkB  = delta_IkBaNFkB  - flux_in_2an
    delta_IkBaNFkBn = delta_IkBaNFkBn + flux_in_2an

    delta_IkBbNFkB  = delta_IkBbNFkB  - flux_in_2bn
    delta_IkBbNFkBn = delta_IkBbNFkBn + flux_in_2bn

    delta_IkBeNFkB  = delta_IkBeNFkB  - flux_in_2en
    delta_IkBeNFkBn = delta_IkBeNFkBn + flux_in_2en

    delta_IkBaNFkBn = delta_IkBaNFkBn - flux_ex_2an
    delta_IkBaNFkB  = delta_IkBaNFkB  + flux_ex_2an

    delta_IkBbNFkBn = delta_IkBbNFkBn - flux_ex_2bn
    delta_IkBbNFkB  = delta_IkBbNFkB  + flux_ex_2bn

    delta_IkBeNFkBn = delta_IkBeNFkBn - flux_ex_2en
    delta_IkBeNFkB  = delta_IkBeNFkB  + flux_ex_2en

    # IkB:NFkB Association [Cytoplasm & Nucleus]
    delta_IkBa      = delta_IkBa - flux_a_c_an
    delta_NFkB      = delta_NFkB - flux_a_c_an
    delta_IkBaNFkB  = delta_IkBaNFkB + flux_a_c_an

    delta_IkBb      = delta_IkBb - flux_a_c_bn
    delta_NFkB      = delta_NFkB - flux_a_c_bn
    delta_IkBbNFkB  = delta_IkBbNFkB + flux_a_c_bn

    delta_IkBe      = delta_IkBe - flux_a_c_en
    delta_NFkB      = delta_NFkB - flux_a_c_en
    delta_IkBeNFkB  = delta_IkBeNFkB + flux_a_c_en

    delta_IkBan     = delta_IkBan - flux_a_n_an
    delta_NFkBn     = delta_NFkBn - flux_a_n_an
    delta_IkBaNFkBn = delta_IkBaNFkBn + flux_a_n_an

    delta_IkBbn     = delta_IkBbn - flux_a_n_bn
    delta_NFkBn     = delta_NFkBn - flux_a_n_bn
    delta_IkBbNFkBn = delta_IkBbNFkBn + flux_a_n_bn

    delta_IkBen     = delta_IkBen - flux_a_n_en
    delta_NFkBn     = delta_NFkBn - flux_a_n_en
    delta_IkBeNFkBn = delta_IkBeNFkBn + flux_a_n_en

    # IkB:NFkB Dissociation [Cytoplasm & Nucleus]
    delta_IkBaNFkB  = delta_IkBaNFkB - flux_d_c_an
    delta_IkBa      = delta_IkBa + flux_d_c_an
    delta_NFkB      = delta_NFkB + flux_d_c_an

    delta_IkBbNFkB  = delta_IkBbNFkB - flux_d_c_bn
    delta_IkBb      = delta_IkBb + flux_d_c_bn
    delta_NFkB      = delta_NFkB + flux_d_c_bn

    delta_IkBeNFkB  = delta_IkBeNFkB - flux_d_c_en
    delta_IkBe      = delta_IkBe + flux_d_c_en
    delta_NFkB      = delta_NFkB + flux_d_c_en

    delta_IkBaNFkBn = delta_IkBaNFkBn - flux_d_n_an
    delta_IkBan     = delta_IkBan + flux_d_n_an
    delta_NFkBn     = delta_NFkBn + flux_d_n_an

    delta_IkBbNFkBn = delta_IkBbNFkBn - flux_d_n_bn
    delta_IkBbn     = delta_IkBbn + flux_d_n_bn
    delta_NFkBn     = delta_NFkBn + flux_d_n_bn

    delta_IkBeNFkBn = delta_IkBeNFkBn - flux_d_n_en
    delta_IkBen     = delta_IkBen + flux_d_n_en
    delta_NFkBn     = delta_NFkBn + flux_d_n_en

    # Free IkB Degradation [Cytoplasm & Nucleus]
    delta_IkBa      = delta_IkBa  - flux_pd_c_a
    delta_IkBb      = delta_IkBb  - flux_pd_c_b
    delta_IkBe      = delta_IkBe  - flux_pd_c_e
    delta_IkBan     = delta_IkBan - flux_pd_n_a
    delta_IkBbn     = delta_IkBbn - flux_pd_n_b
    delta_IkBen     = delta_IkBen - flux_pd_n_e

    # IkB:NFkB Degradation [Cytoplasm & Nucleus]
    delta_IkBaNFkB  = delta_IkBaNFkB - flux_pd_c_2an
    delta_NFkB      = delta_NFkB + flux_pd_c_2an

    delta_IkBbNFkB  = delta_IkBbNFkB - flux_pd_c_2bn
    delta_NFkB      = delta_NFkB + flux_pd_c_2bn

    delta_IkBeNFkB  = delta_IkBeNFkB - flux_pd_c_2en
    delta_NFkB      = delta_NFkB + flux_pd_c_2en

    delta_IkBaNFkBn = delta_IkBaNFkBn - flux_pd_n_2an
    delta_NFkBn     = delta_NFkBn + flux_pd_n_2an

    delta_IkBbNFkBn = delta_IkBbNFkBn - flux_pd_n_2bn
    delta_NFkBn     = delta_NFkBn + flux_pd_n_2bn

    delta_IkBeNFkBn = delta_IkBeNFkBn - flux_pd_n_2en
    delta_NFkBn     = delta_NFkBn + flux_pd_n_2en

    # IKK Mediated IkB Degradation
    delta_IkBa      = delta_IkBa - flux_ph_c_a
    delta_IkBb      = delta_IkBb - flux_ph_c_b
    delta_IkBe      = delta_IkBe - flux_ph_c_e

    # IKK Mediated IkB:NFkB Degradation
    delta_IkBaNFkB  = delta_IkBaNFkB - flux_ph_c_an
    delta_NFkB      = delta_NFkB + flux_ph_c_an

    delta_IkBbNFkB  = delta_IkBbNFkB - flux_ph_c_bn
    delta_NFkB      = delta_NFkB + flux_ph_c_bn

    delta_IkBeNFkB  = delta_IkBeNFkB - flux_ph_c_en
    delta_NFkB      = delta_NFkB + flux_ph_c_en

    # Upstream Pathway  [LPS] --> IKK activation

    # Ligand binding
    delta_LPSo      = delta_LPSo - flux_b_LPS
    delta_LPS       = delta_LPS  + flux_b_LPS

    delta_LPS       = delta_LPS       -   flux_b_TLR4LPS
    delta_TLR4       = delta_TLR4       -   flux_b_TLR4LPS
    delta_TLR4LPS   = delta_TLR4LPS   +   flux_b_TLR4LPS
    delta_LPS       = delta_LPS       +   flux_ub_TLR4LPS
    delta_TLR4       = delta_TLR4     +   flux_ub_TLR4LPS
    delta_TLR4LPS   = delta_TLR4LPS   -   flux_ub_TLR4LPS

    delta_LPSen     = delta_LPSen     -   flux_b_TLR4LPSen
    delta_TLR4en     = delta_TLR4en     -   flux_b_TLR4LPSen
    delta_TLR4LPSen = delta_TLR4LPSen +   flux_b_TLR4LPSen
    delta_LPSen     = delta_LPSen     +   flux_ub_TLR4LPSen
    delta_TLR4en     = delta_TLR4en   +   flux_ub_TLR4LPSen
    delta_TLR4LPSen = delta_TLR4LPSen -   flux_ub_TLR4LPSen

    # Activation
    delta_MyD88     = delta_MyD88     -   flux_a_MyD88
    delta_MyD88s    = delta_MyD88s    +   flux_a_MyD88
    delta_MyD88     = delta_MyD88     +   flux_i_MyD88
    delta_MyD88s    = delta_MyD88s    -   flux_i_MyD88

    delta_TRIF      = delta_TRIF      -   flux_a_TRIF
    delta_TRIFs     = delta_TRIFs     +   flux_a_TRIF
    delta_TRIF      = delta_TRIF      +   flux_i_TRIF
    delta_TRIFs     = delta_TRIFs     -   flux_i_TRIF

    # shuttling
    delta_LPS       = delta_LPS       -   flux_in_LPS
    delta_LPSen     = delta_LPSen     +   flux_in_LPS
    # delta_LPSo       = delta_LPSo     +   flux_out_LPS
    delta_LPS       = delta_LPS       +   flux_out_LPS
    delta_LPSen     = delta_LPSen     -   flux_out_LPS

    delta_TLR4      = delta_TLR4      -   flux_in_TLR4
    delta_TLR4en    = delta_TLR4en    +   flux_in_TLR4
    delta_TLR4      = delta_TLR4      +   flux_out_TLR4
    delta_TLR4en    = delta_TLR4en    -   flux_out_TLR4

    delta_TLR4LPS      = delta_TLR4LPS      -   flux_in_TLR4LPS
    delta_TLR4LPSen    = delta_TLR4LPSen    +   flux_in_TLR4LPS
    delta_TLR4LPS      = delta_TLR4LPS      +   flux_out_TLR4LPS
    delta_TLR4LPSen     = delta_TLR4LPSen   -   flux_out_TLR4LPS

    # generation & degradation
    delta_TLR4      = delta_TLR4 + flux_g_TLR4
    delta_TLR4      = delta_TLR4 - flux_d_TLR4
    delta_TLR4en     = delta_TLR4en - flux_d_TLR4en
    delta_LPSen      = delta_LPSen - flux_d_LPSen
    delta_TLR4LPSen      = delta_TLR4LPSen - flux_d_TLR4LPSen
    # -- modified 0515
    delta_TLR4LPS      = delta_TLR4LPS - flux_d_TLR4LPS

    # TRAF6
    delta_TRAF6       = delta_TRAF6  - flux_a_TRAF6_MyD88s
    delta_TRAF6s      = delta_TRAF6s + flux_a_TRAF6_MyD88s
    delta_TRAF6       = delta_TRAF6  - flux_a_TRAF6_TRIFs
    delta_TRAF6s      = delta_TRAF6s + flux_a_TRAF6_TRIFs
    delta_TRAF6       = delta_TRAF6  + flux_i_TRAF6
    delta_TRAF6s      = delta_TRAF6s - flux_i_TRAF6

    # IKKK Regulation
    delta_IKKK      = delta_IKKK     + flux_IKKK_on
    delta_IKKK_off  = delta_IKKK_off - flux_IKKK_on

    delta_IKKK      = delta_IKKK     + flux_IKKK_on_TRAF6s
    delta_IKKK_off  = delta_IKKK_off - flux_IKKK_on_TRAF6s

    delta_IKKK_off  = delta_IKKK_off + flux_IKKK_off
    delta_IKKK      = delta_IKKK     - flux_IKKK_off

    # IKK Regulation
    delta_IKK       = delta_IKK     + flux_IKK_on
    delta_IKK_off   = delta_IKK_off - flux_IKK_on

    delta_IKK       = delta_IKK     + flux_IKK_on_IKKK
    delta_IKK_off   = delta_IKK_off - flux_IKK_on_IKKK

    delta_IKK       = delta_IKK     - flux_IKK_off
    delta_IKK_off   = delta_IKK_off + flux_IKK_off

    delta_IKK       = delta_IKK     - flux_IKK_off_i
    delta_IKK_i     = delta_IKK_i   + flux_IKK_off_i

    delta_IKK_off   = delta_IKK_off + flux_IKKi
    delta_IKK_i     = delta_IKK_i   - flux_IKKi

    ## TNR module [NEW] + A20 [NEW]
    # A20 Expression Transcrviption
    delta_a20t      = delta_a20t + flux_rsu_a20
    delta_a20t      = delta_a20t + flux_rsr_a20
    delta_a20t      = delta_a20t - flux_rd_a20

    # Translation
    delta_a20       = delta_a20  + flux_ps_c_a20

    # Degradation
    delta_a20       = delta_a20 - flux_pd_c_a20

    # Upstream Pathway  [TNF] --> IKK activation

    # TNF-independent Activation

    # tnfrm trimerization 3tnfrm<--> TNFR
    delta_tnfrm     = delta_tnfrm - 3 * flux_a_tnfrm
    delta_TNFR      = delta_TNFR +   flux_a_tnfrm

    delta_tnfrm     = delta_tnfrm + 3 * flux_d_TNFR
    delta_TNFR      = delta_TNFR -   flux_d_TNFR

    # TNF Receptor Metabolism
    delta_TNFR     = delta_TNFR     - flux_i_TNFR
    delta_TNFRtnf  = delta_TNFRtnf  - flux_i_TNFRtnf

    delta_tnfrm    = delta_tnfrm    - flux_pd_tnfrm
    delta_tnfrm    = delta_tnfrm    + flux_syn_tnfrm

    # Complex I Generation  TTR + TNFR <--> C1_off; C1
    delta_C1_off    = delta_C1_off  + flux_a_C1_off
    delta_TTR       = delta_TTR     - flux_a_C1_off
    delta_TNFR      = delta_TNFR    - flux_a_C1_off

    delta_C1_off    = delta_C1_off  - flux_d_C1_off
    delta_TTR       = delta_TTR     + flux_d_C1_off
    delta_TNFR      = delta_TNFR    + flux_d_C1_off

    delta_C1        = delta_C1      - flux_d_C1
    delta_TTR       = delta_TTR     + flux_d_C1
    delta_TNFR      = delta_TNFR    + flux_d_C1

    # C1_off  <--> C1
    delta_C1_off    = delta_C1_off  - flux_a_C1
    delta_C1        = delta_C1      + flux_a_C1

    delta_C1_off    = delta_C1_off  + flux_C1_A20
    delta_C1        = delta_C1      - flux_C1_A20

    delta_C1_off    = delta_C1_off  + flux_C1_off
    delta_C1        = delta_C1      - flux_C1_off

    # C1 & C1_off internalization [treated as deg]
    delta_C1        = delta_C1      - flux_i_C1
    delta_C1_off    = delta_C1_off  - flux_i_C1_off

    # TNF-Dependent Activation

    # TNF degradation
    delta_TNF       = delta_TNF     - flux_pd_m_tnf

    # TNF:TNFR binding
    delta_TNFRtnf      = delta_TNFRtnf    +   flux_a_tnfrmtnf
    delta_tnfrm     = delta_tnfrm   - 3 * flux_a_tnfrmtnf
    delta_TNF       = delta_TNF     -   flux_a_tnfrmtnf

    delta_TNFR      = delta_TNFR    - flux_a_TNFRtnf
    delta_TNF       = delta_TNF     - flux_a_TNFRtnf
    delta_TNFRtnf   = delta_TNFRtnf + flux_a_TNFRtnf

    delta_TNFR      = delta_TNFR    + flux_d_TNFRtnf
    delta_TNF       = delta_TNF     + flux_d_TNFRtnf
    delta_TNFRtnf   = delta_TNFRtnf - flux_d_TNFRtnf

    # Complex I Generation  TTR + TNFRtnf <--> C1tnf_off; C1tnf
    delta_C1tnf_off = delta_C1tnf_off   + flux_a_C1tnf_off
    delta_TTR       = delta_TTR         - flux_a_C1tnf_off
    delta_TNFRtnf   = delta_TNFRtnf     - flux_a_C1tnf_off

    delta_C1tnf_off = delta_C1tnf_off   - flux_d_C1tnf_off
    delta_TTR       = delta_TTR         + flux_d_C1tnf_off
    delta_TNFRtnf   = delta_TNFRtnf     + flux_d_C1tnf_off

    delta_C1tnf     = delta_C1tnf       - flux_d_C1tnf
    delta_TTR       = delta_TTR         + flux_d_C1tnf
    delta_TNFRtnf   = delta_TNFRtnf     + flux_d_C1tnf

    # C1tnf_off <--> C1tnf
    delta_C1tnf_off = delta_C1tnf_off   - flux_a_C1tnf
    delta_C1tnf     = delta_C1tnf       + flux_a_C1tnf

    delta_C1tnf_off = delta_C1tnf_off   + flux_C1tnf_off
    delta_C1tnf     = delta_C1tnf       - flux_C1tnf_off

    delta_C1tnf_off = delta_C1tnf_off   + flux_C1tnf_A20
    delta_C1tnf     = delta_C1tnf       - flux_C1tnf_A20

    # C1tnf & C1tnf_off internalization [treated as deg]
    delta_C1tnf_off = delta_C1tnf_off   - flux_i_C1tnf_off
    delta_C1tnf     = delta_C1tnf       - flux_i_C1tnf

    # C1tnf_off;C1tnf <--> C1_off;C1 + TNF
    delta_C1tnf_off = delta_C1tnf_off   - flux_d_tnf_C1_off
    delta_C1_off    = delta_C1_off      + flux_d_tnf_C1_off
    delta_TNF       = delta_TNF         + flux_d_tnf_C1_off

    delta_C1tnf_off = delta_C1tnf_off   + flux_a_tnf_C1_off
    delta_C1_off    = delta_C1_off      - flux_a_tnf_C1_off
    delta_TNF       = delta_TNF         - flux_a_tnf_C1_off

    delta_C1tnf     = delta_C1tnf       - flux_d_tnf_C1
    delta_C1        = delta_C1          + flux_d_tnf_C1
    delta_TNF       = delta_TNF         + flux_d_tnf_C1

    delta_C1tnf     = delta_C1tnf       + flux_a_tnf_C1
    delta_C1        = delta_C1          - flux_a_tnf_C1
    delta_TNF       = delta_TNF         - flux_a_tnf_C1

    # IKKK Regulation

    delta_IKKK      = delta_IKKK     + flux_IKKK_on_C1
    delta_IKKK_off  = delta_IKKK_off - flux_IKKK_on_C1

    delta_IKKK      = delta_IKKK     + flux_IKKK_on_C1tnf
    delta_IKKK_off  = delta_IKKK_off - flux_IKKK_on_C1tnf

    ## ADDED
    # A20 mediated TRAF6 inactivation
    delta_TRAF6    = delta_TRAF6    + v.TP[15] * a20 * TRAF6s
    delta_TRAF6s    = delta_TRAF6s  - v.TP[15] * a20 * TRAF6s

    ## Save concentrations for next time step
    delta[1]      = delta_IkBa
    delta[2]      = delta_IkBan
    delta[3]      = delta_IkBaNFkB
    delta[4]      = delta_IkBaNFkBn
    delta[5]      = delta_IkBat
    delta[6]      = delta_IkBb
    delta[7]      = delta_IkBbn
    delta[8]      = delta_IkBbNFkB
    delta[9]      = delta_IkBbNFkBn
    delta[10]     = delta_IkBbt
    delta[11]     = delta_IkBe
    delta[12]     = delta_IkBen
    delta[13]     = delta_IkBeNFkB
    delta[14]     = delta_IkBeNFkBn
    delta[15]     = delta_IkBet

    delta[16]     = delta_NFkB
    delta[17]     = delta_NFkBn
    #
    delta[18]     = delta_LPSo
    delta[19]     = delta_LPS
    delta[20]     = delta_LPSen
    delta[21]     = delta_TLR4
    delta[22]     = delta_TLR4en
    delta[23]     = delta_TLR4LPS
    delta[24]     = delta_TLR4LPSen
    delta[25]     = delta_MyD88
    delta[26]     = delta_MyD88s
    delta[27]     = delta_TRIF
    delta[28]     = delta_TRIFs
    delta[29]     = delta_TRAF6
    delta[30]     = delta_TRAF6s
    delta[31]     = delta_IKKK_off
    delta[32]     = delta_IKKK
    delta[33]     = delta_IKK_off
    delta[34]     = delta_IKK
    delta[35]     = delta_IKK_i

    # TNF component
    delta[36]     = delta_TNFnas
    delta[37]     = delta_TNFmRNA
    delta[38]     = delta_TNFpro
    delta[39]     = delta_TNF

    # TNFR part
    delta[40]     = delta_tnfrm
    delta[41]     = delta_TNFR
    delta[42]     = delta_TNFRtnf
    delta[43]     = delta_C1
    delta[44]     = delta_C1_off
    delta[45]     = delta_C1tnf
    delta[46]     = delta_C1tnf_off
    delta[47]     = delta_TTR

    # A20
    delta[48]     = delta_a20
    delta[49]     = delta_a20t

    # # if (v.flag_noTnfFeedback) && (t>0)
    # #     v.TP[8]=kdeg_tnf
    # # end
end

max_stars_repo_name: jiweiqi/NFkB

max_issues_repo_name: jiweiqi/NFkB

max_forks_repo_name: jiweiqi/NFkB


text:
"""
    write_dot_file(g, file; attributes, path, colors)

Export graph `g` to DOT-Format and store it in file `file`.

#### Arguments
- `g::AbstractSimpleWeightedGraph`: a graph representation to export
- `filename::AbstractString`: the filename to store (i.e. "graph.dot")
- (optional) `attributes::AttributeDict`: Dot-Attributes like node color stored in a dictionary
- (optional) `path = []`: Int-Array of nodes. Nodes and their edges are drawn in red color (i.e. shortest path)
- (optional) `colors = zeros(Int, nv(mat))`: Components-colors vector representated by a color number each node.
"""
function write_dot_file(graph::AbstractSimpleWeightedGraph, filename::AbstractString;
    attributes=GraphvizAttributes(graph), path=[], colors=zeros(Int, nv(graph)))

    if !isempty(path)
        color_path!(attributes, path, graph)
    end

    if (!is_all_zero(colors))
        color_nodes!(attributes, colors)
    end

    open(filename, "w") do f
        dot(graph, f, attributes)
    end
end

# internal function to get the dot representation of a graph as a string.
function string_dot(graph::AbstractSimpleWeightedGraph, attributes=GraphvizAttributes(g), path=[], colors=zeros(Int, nv(g)))

    if !isempty(path)
        color_path!(attributes, path, graph)
    end

    if (!is_all_zero(colors))
        color_nodes!(attributes, colors)
    end

    str = IOBuffer()
    dot(graph, str, attributes)
    String(take!(str)) #takebuf_string(str)
end

function save_dot_as(graph::AbstractSimpleWeightedGraph, filename::AbstractString;
    attributes=GraphvizAttributes(graph), path=[], colors=zeros(Int, nv(graph)))

    ## ToDo

end

max_stars_repo_name: tragisch/PlotGraphviz

max_issues_repo_name: tragisch/PlotGraphviz


text:
############################################################################################
"""
$(TYPEDEF)

Utility struct to help assign boundaries to parameter.

# Fields
$(TYPEDFIELDS)
"""
struct CorrelationMatrix <: AbstractConstraint
    function CorrelationMatrix()
        return new()
    end
end

############################################################################################
function construct_transform(constraint::CorrelationMatrix, val)
    transform = Bijectors.CorrBijector()
    return transform, Bijectors.inverse(transform)
end

############################################################################################
function _check(
    _rng::Random.AbstractRNG,
    constraint::CorrelationMatrix,
    val::Matrix{R},
) where {R<:Real}
    ArgCheck.@argcheck all(LinearAlgebra.diag(val) .== 1.0)
    return true
end

############################################################################################
#= !NOTES:
    Unconstrained will always be 0 everywhere except upper diagonal elements. All other entries do not matter for constrain/unconstrain.
    Constrained will always have unit variance.
=#
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::UnflattenStrict,
    constraint::C,
    x::Matrix{R},
) where {
    T<:AbstractFloat,
    F<:FlattenTypes,
    R<:Real,
    C<:Union{CorrelationMatrix, Distributions.LKJ, Bijectors.CorrBijector}
}
    #!NOTE: CorrBijector seems to unconstrain to a Upper Diagonal Matrix
    idx_upper = tag(x, true, false)
    len = length(x)
    len_unflat = sum(idx_upper)
    #!NOTE: Buffer should be of type R, not T, as we want same type back afterwards
    function flatten_CorrMatrix(x::AbstractMatrix{R}) where {R<:Real}
        ArgCheck.@argcheck length(x) == len
        return Vector{T}(flatten_Symmetric(x, idx_upper))
    end
    buffer_unflat = ones(R, size(x))
    function CorrMatrix_from_vec(v::Union{<:Real,AbstractVector{<:Real}})
        ArgCheck.@argcheck length(v) == len_unflat
        return Symmetric_from_flatten!(buffer_unflat, v, idx_upper)
    end
    return flatten_CorrMatrix, CorrMatrix_from_vec
end
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::UnflattenFlexible,
    constraint::C,
    x::Matrix{R},
) where {
    T<:AbstractFloat,
    F<:FlattenTypes,
    R<:Real,
    C<:Union{CorrelationMatrix, Distributions.LKJ, Bijectors.CorrBijector}
}
    #!NOTE: CorrBijector seems to unconstrain to a Upper Diagonal Matrix
    idx_upper = tag(x, true, false)
    len = length(x)
    len_unflat = sum(idx_upper)
    function flatten_CorrMatrix_AD(x::AbstractMatrix{R}) where {R<:Real}
        ArgCheck.@argcheck length(x) == len
        return Vector{R}(flatten_Symmetric(x, idx_upper))
    end
    function CorrMatrix_from_vec_AD(v::Union{<:Real,AbstractVector{<:Real}})
        ArgCheck.@argcheck length(v) == len_unflat
        return Symmetric_from_flatten!(ones(eltype(v), size(x)), v, idx_upper)
    end
    return flatten_CorrMatrix_AD, CorrMatrix_from_vec_AD
end

############################################################################################
#=
!NOTE: For this constraint, we use a Bijector as Transformer, so we do not need to add a new functors
1. MyTransformer <: AbstractTransformer, MyInverseTransformer <: AbstractTransformer
2. define a function construct_transform(MyConstraint, val) -> MyTransformer, MyInverseTransformer
3. overload unconstrain and log_abs_det_jac on MyTransformer.
4. overload constrain on MyInverseTransformer.
=#

############################################################################################
#Export
export
    CorrelationMatrix,
    construct_flatten,
    construct_transform,
    constrain,
    unconstrain,
    log_abs_det_jac


text:
# function any_method_exists(func::Symbol, types, mod::Module=Base)::Bool
#   subtypes_any = subtypes(Any)
#   for concrete_types in product((t === Any ? subtypes_any : [t] for t in types)...)
#     try
#       if eval(mod, :(method_exists($func, $concrete_types)))
#         return true  # short cycle as soon as possible
#       end
#
#     catch exc
#       if isa(exc, TypeError)
#         # This happens for the Vararg type, hence Varargs are currently regarded as not existent which is defintely worth some improvements
#         continue
#       else
#         rethrow(exc)
#       end
#     end
#   end
#   false
# end


text:
cube, angles = BetaPictoris[:cube, :pa]

@testset "Interface" begin
    @test TPCA(5) == TPCA(ncomps=5)
    @test isnothing(TPCA().ncomps)
end

@testset "Decomposition" begin

    # get sizes correct for ncomps
    for N in [1, 3, 5]
        A, w = ADI.fit(TPCA(N), cube)
        @test size(A) == (N, 101 * 101)
        @test size(w) == (size(cube, 1), N)
    end

    # default is to use whole cube
    S = reconstruct(TPCA(), cube)
    @test size(S) == size(cube)
end

# @testset "ADI Trivial" begin
#     cube = ones(10, 100, 100)
#     angles = zeros(10)

#     reduced_5 = TPCA(5)(cube, angles)
#     reduced_10 = TPCA(10)(cube, angles)

#     @test size(reduced_5) == size(reduced_10) == (100, 100)
#     @test all(x -> isapprox(x, 0, atol = 1e-9), reduced_5)
#     @test all(x -> isapprox(x, 0, atol = 1e-9), reduced_10)
# end

# @testset "RDI Trivial" begin
#     data = randn(rng, 30, 101, 101)
#     angles = sort!(90rand(rng, 30)) |> normalize_par_angles

#     S = subtract(TPCA(2), data; ref=zeros(10, 101, 101))
#     @test S ≈ data rtol=2e-2
# end


text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mc01md(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    DPMAX = 20
    P = Array{Float64,1}(undef, DPMAX+1)
    Q = Array{Float64,1}(undef, DPMAX+1)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    DP = parse(BlasInt, vs[1])
    ALPHA = parse(Float64, replace(vs[2],'D'=>'E'))
    K = parse(BlasInt, vs[3])
if ( DP<=-1 || DP>DPMAX )
else

    vs = String[]
    _isz = DP+1
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    P[1:_isz] .= parsex.(Float64, vs)

# interp call 1

    INFO = SLICOT.mc01md!(DP, ALPHA, K, P, Q)
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
end # if
end # if
    close(f)
end # run_X()


text:
using PyCall
# @pyimport stormpy

struct ModelCheckingResult{S}
    mdp::MDP{S}
    labels::Dict{S, Vector{String}}
    property::String
    result::PyObject
end

struct Scheduler{S, A} <: Policy
    mdp::MDP{S, A}
    _scheduler::PyObject
    scheduler::Vector{A}
    action_map::Vector{A}
end

function Scheduler(mdp::MDP{S, A}, result::ModelCheckingResult) where {S,A}
    @assert result.result[:has_scheduler]
    py_scheduler = result.result[:scheduler]
    return Scheduler(mdp, py_scheduler)
end

function Scheduler(mdp::MDP{S, A}, py_scheduler::PyObject) where {S,A}
    action_map = ordered_actions(mdp)
    scheduler = Vector{A}(undef, n_states(mdp))
    for i=1:n_states(mdp)
        choice = py_scheduler[:get_choice](i-1)
        ai = choice[:get_deterministic_choice]() + 1
        scheduler[i] = action_map[ai]
    end
    return Scheduler(mdp, py_scheduler, scheduler, action_map)
end

function POMDPs.action(policy::Scheduler{S, A}, s::S) where {S,A}
    si = stateindex(policy.mdp, s)
    return policy.scheduler[si]
end

"""
Write the MDP model and its associated labels to a file and parse them using storm. It then runs the reachability analysis by computing for each state,
the probability of satisfying the desired property. To learn more about how to specify properties, look at the storm documentation.
This function returns a `ModelCheckingResult` object, from which the probability of satisfaction can be extracted as well as the scheduler maximizing the probability
of success.
**Arguments:**
- `mdp::MDP`, the MDP models defined according to the POMDPs.jl interface
- `labeling::Dict{S, Vector{String}}` , a dictionnary mapping the states of the MDP to their labels
- `property::String`, a LTL property to verify. (follow storm property definition)
- `transition_file_name::String = "mdp.tra"` the file name to read from or write the MDP model
- `labels_file_name::String = "mdp.lab"` the file name to read from or write the state labels
- `overwrite::Bool = false` if set to true, it will remove the existing transition and labels file and write new ones.

`model_checking{S}(mdp::MDP{S}, labeling::Dict{S, Vector{String}}, property::String; transition_file_name::String = "mdp.tra", labels_file_name::String = "mdp.lab", overwrite::Bool = false)`
"""
function model_checking(mdp::MDP{S}, labeling::Dict{S, Vector{String}}, property::String;
                        transition_file_name::String = "mdp.tra",
                        labels_file_name::String = "mdp.lab",
                        overwrite::Bool = false) where S
    model = parse_mdp_model(mdp, labeling, transition_file_name, labels_file_name, overwrite)
    properties = stormpy.parse_properties(property)
    result = stormpy.model_checking(model, properties[1],  only_initialstates=false, extract_scheduler=true)
    @assert result[:result_for_all_states]
    return ModelCheckingResult(mdp, labeling, property, result)
end

"""
    parse an MDP model with Storm, if the specified transition and labels file do not exist it will write the files first
"""
function parse_mdp_model(mdp::MDP{S}, labeling::Dict{S, Vector{String}},
                        transition_file_name::String = "mdp.tra",
                        labels_file_name::String = "mdp.lab",
                        overwrite::Bool = false) where S
        if overwrite
            warn("overwriting potential existing files!")
            try
                run(`rm $transition_file_name`)
                run(`rm $labels_file_name`)
            catch
            end
        end
        if !isfile(transition_file_name)
            write_mdp_transition(mdp, transition_file_name)
        end
        if !isfile(labels_file_name)
            write_mdp_labels(mdp, labeling, labels_file_name)
        end
        model = parse_mdp_model(transition_file_name, labels_file_name)
        return model
end

function parse_mdp_model(transition_file_name::String, labels_file_name::String)
    return stormpy.build_sparse_model_from_explicit(transition_file_name, labels_file_name)
end

struct StormPolicy{M<:MDP} <: Policy
    mdp::M
    risk_vec::Vector{Float64}
    risk_mat::Array{Float64, 2}
end

function StormPolicy(mdp::MDP, result::ModelCheckingResult)
    risk_vec = get_proba(mdp, result)
    risk_mat = get_state_action_proba(mdp, risk_vec)
    return StormPolicy(mdp, risk_vec, risk_mat)
end

function value_vector(policy::StormPolicy, s)
    si = stateindex(policy.mdp, s)
    return policy.risk_mat[si, :]
end

function value(policy::StormPolicy, s)
    si = stateindex(policy.mdp, s)
    return policy.risk_vec[si]
end

function action(policy::StormPolicy, s)
    acts = actions(mdp)
    vals = value_vector(policy, s)
    return acts[argmax(vals)]
end

"""
extract a vector P of size |S| where P(s) is the probability of satisfying a property
when starting in state s
`get_proba(mdp::MDP, result::ModelCheckingResult)`
"""
function get_proba(mdp::MDP, result::ModelCheckingResult)
    P = zeros(n_states(mdp))
    for (i, val) in enumerate(result.result[:get_values]())
        P[i] = val
    end
    return P
end

"""
Returns a matrix of dimension |S|x|A| where each element is the probability of satisfying an LTL formula for a given state action pair.
This algorithm is mathematically sound only for basic LTL property like "!a U b" !!! (should technically be over the product MDP)
Arguments:
- `mdp::MDP` the MDP model
- `result::ModelCheckingResult` result from model checking
 `get_state_action_proba(mdp::MDP, P::Vector{Float64})`
"""
function get_state_action(mdp::MDP, result::ModelCheckingResult)
    P = get_proba(mdp, result)
    P_sa = get_state_action_proba(mdp, P)
    return P_sa
end
function get_state_action_proba(mdp::MDP, P::Vector{Float64})
    P_sa = zeros(n_states(mdp), n_actions(mdp))
    states = ordered_states(mdp)
    actions = ordered_actions(mdp)
    for (si, s) in enumerate(states)
        P[si] == 0. ? continue : nothing
        for (ai, a) in enumerate(actions)
            dist = transition(mdp, s, a)
            for (sp, p) in  weighted_iterator(dist)
                p == 0.0 ? continue : nothing # skip if zero prob
                spi = stateindex(mdp, sp)
                P_sa[si, ai] += p * P[spi]
            end
        end
    end
    return P_sa
end


text:
import Ipopt
import JuMP
import PowerModels
import SymbolicAD

function power_model(case::String)
    pm = PowerModels.instantiate_model(
        joinpath(dirname(@__DIR__), "data", case),
        PowerModels.ACPPowerModel,
        PowerModels.build_opf,
    )
    return pm.model
end

import ProfileView

model = power_model("pglib_opf_case118_ieee.m")
JuMP.set_optimizer(model, Ipopt.Optimizer)
JuMP.set_optimizer_attribute(model, "print_timing_statistics", "yes")
JuMP.set_optimize_hook(model, SymbolicAD.optimize_hook)
JuMP.optimize!(model)

ProfileView.@profview JuMP.optimize!(model)


text:
@testset "PointCloud Model tests" begin
    x_test = rand(Float32, 3, 64, 2)
    _dims = size(x_test, 2)

    @testset "PointNet Model tests" begin
        MODEL = PointNet
        @info "Testing $(MODEL)..."
        for num_classes in [10, 40]
            @testset "PointNet num_classes: $(num_classes)" begin
                model = MODEL(num_classes)

                y_test = model(x_test)

                # Testing Forward Propagation
                @test y_test isa AbstractArray
                @test size(y_test) == (num_classes, size(x_test, 3))

                # Testing Backward Propagation
                # @test gradtest(x_test -> model(x_test), x_test)
            end
        end
    end

    @testset "DGCNN Model tests" begin
        MODEL = DGCNN
        @info "Testing $(MODEL)..."
        for num_classes in [10, 40]
            @testset "DGCNN num_classes: $(num_classes)" begin
                model = MODEL(num_classes, 10, _dims)

                y_test = model(x_test)

                # Testing Forward Propagation
                @test y_test isa AbstractArray
                @test size(y_test) == (num_classes, size(x_test, 3))

                # Testing Backward Propagation
                # @test gradtest(x_test -> model(x_test), x_test)
            end
        end
    end
end


text:
using Oanda
using Base.Test
@testset "pairs" begin
    #
    #
end
@testset "pricing" begin
    @test typeof(price("USD_JPY")) == HttpCommon.Response
    @test typeof(price_stream("USD_JPY",true)) == Requests.ResponseStream{MbedTLS.SSLContext}
    @test typeof(price(["USD_JPY"])) == HttpCommon.Response
    @test typeof(price_stream(["USD_JPY"],true)) == Requests.ResponseStream{MbedTLS.SSLContext}
end

max_stars_repo_name: cadwallader/Oanda

max_issues_repo_name: cadwallader/Oanda

max_forks_repo_name: cadwallader/Oanda


text:
import GaussQuadrature.legendre
using FractionalTimeDG
using ArgCheck
using OffsetArrays

#include("FEM1D.jl")

const L = 2.0
const max_t = 2.0
const C0 = 1.0
const Cf = 2.0
α = 0.6

u0(x) = C0 * x * (L-x)
f(x,t) = Cf * t * exp(-t)

function refsoln(x, t, α, N)
    return Bromich_integral(t, z -> uhat(x, z, α), N)
end

function uhat(x::T, z::Complex{T}, α::T) where T <: AbstractFloat
    ω = z^(α/2)
    ρ1(x) = ( ω * x * (L-x) - 2/ω ) * cosh(ω*x) + (2x-L) * sinh(ω*x) + 2 / ω
    ρ2(x) = cosh(ω*x) - 1
    return ( (C0/ω) * ( ρ1(x) * sinh(ω*(L-x)) + ρ1(L-x) * sinh(ω*x) )
           + Cf/(z+1)^2 * ( ρ2(x) * sinh(ω*(L-x)) + ρ2(L-x) * sinh(ω*x) )
          ) / ( z * sinh(ω*L) )
end

"""
    Bromich_integral(t, F, N)

Evaluates `1/(2πi)` times the integral of `exp(zt)F(z) dz` along a
Hankel contour, assuming `F(z)` is analytic in the cut plane
`|arg(z)|<π` and `F(conj(z)) = conj(F(z))`.
"""
function Bromich_integral(t::Float64, F::Function, N::Integer)
    μ = 4.492075 * N / t
    ϕ = 1.172104
    h = 1.081792 / N
    z0, w0 = integration_contour(0.0, μ, ϕ)
    Σ = real( w0 * exp(z0*t) * F(z0) ) / 2
    for n = 1:N
        zn, wn = integration_contour(n*h, μ, ϕ)
        Σ += real( wn * exp(zn*t) * F(zn) )
    end
    return h * Σ
end

function integration_contour(u, μ, ϕ)
    z = μ * ( 1 + sin(Complex(-ϕ, u)) )
    w = (μ/π) * cos(Complex(-ϕ, u))
    return z, w
end

function L2error(U::Array{Float64}, u::Function, t::OffsetArray{Float64},
                 τ::AbstractVector{Float64},
                 grid::SpatialGrid{Float64}, Mx::Integer)
    ξ, w = legendre(Mx)
    xvals, pcwise_t, pcwise_U = evaluate_pcwise_poly(U, t, grid, τ, ξ)
    pcwise_err = evaluate_pcwise_error(pcwise_U, u, pcwise_t, xvals)

    x = grid.x
    pts_per_time_interval = length(τ)
    Nt = length(t) - 1
    Nx = length(x) - 1
    err = Array{Float64}(undef, pts_per_time_interval, Nt)
    for nt = 1:Nt
        for i = 1:pts_per_time_interval
            m = 0
            s = 0.0
            for nx = 1:Nx
                hx = x[nx] - x[nx-1]
                for j = 1:Mx
                    m += 1
                    s += (hx/2) * w[j] * pcwise_err[m,i,nt]^2
                end
            end
            err[i,nt] = sqrt(s)
        end
    end
    return pcwise_t, err
end

function evaluate_jumps(U::Array{Float64}, U0::Vector{Float64},
                        t::OffsetVector{Float64},
                        grid::SpatialGrid{Float64}, Mx::Integer)
    x = grid.x
    Nx = length(x) - 1
    ξ, w = legendre(Mx)
    τ = [-1.0, 1.0]
    xvals, pcwise_t, pcwise_U = evaluate_pcwise_poly(U, t, grid, τ, ξ)
    U0vals = similar(xvals)
    evaluate_fem1d_soln!(U0vals, U0, grid, ξ)
    Nt = length(t) - 1
    jump = OffsetArray{Float64}(undef, 0:Nt-1)
    m = 0
    s = 0.0
    for nx = 1:Nx
        hx = x[nx] - x[nx-1]
        for j = 1:Mx
            m += 1
            s += (hx/2) * w[j] * ( pcwise_U[m,1,1] - U0vals[m] )^2
        end
    end
    jump[0] = sqrt(s)
    for nt = 1:Nt-1
        m = 0
        s = 0.0
        for nx = 1:Nx
            hx = x[nx] - x[nx-1]
            for j = 1:Mx
                m += 1
                s += (hx/2) * w[j] * ( pcwise_U[m,1,nt+1] - pcwise_U[m,2,nt] )^2
            end
        end
        jump[nt] = sqrt(s)
    end
    return jump
end

function reconstruction(U::Array{Float64}, U0::Vector{Float64},
                        store::Store{Float64})
    NgDoF, rt, Nt = size(U)
    Uhat = Array{Float64}(undef, NgDoF, rt+1, Nt)
    jumpU = OffsetArray{Float64}(undef, 1:NgDoF, 0:Nt-1)
    rmax = store.rmax
    pow = OffsetArray{Float64}(undef, 0:rmax)
    pow[0:2:end] .=  1.0
    pow[1:2:end] .= -1.0
    @argcheck rt+1 ≤ rmax
    U_left  = zeros(NgDoF)
    U_right = zeros(NgDoF)
    for j = 1:rt
        U_left .+= pow[j-1] * U[:,j,1]
        U_right .+= U[:,j,1]
    end
    jumpU[:,0] .= U_left - U0
    Uhat[:,1:rt,1] .= U[:,1:rt,1]
    Uhat[:,rt,1] .+= pow[rt] * jumpU[:,0] / 2
    Uhat[:,rt+1,1] .= - pow[rt] * jumpU[:,0] / 2
    for n = 2:Nt
        fill!(U_left, 0.0)
        for j = 1:rt
            U_left .+= pow[j-1] * U[:,j,n]
        end
        jumpU[:,n-1] .= U_left - U_right
        Uhat[:,1:rt,n] .= U[:,1:rt,n]
        Uhat[:,rt,n]  .+=  pow[rt] * jumpU[:,n-1] / 2
        Uhat[:,rt+1,n] .= -pow[rt] * jumpU[:,n-1] / 2
        fill!(U_right, 0.0)
        for j = 1:rt
            U_right .+= U[:,j,n]
        end
    end
    return Uhat, jumpU
end


text:
module LightGraphicalModels
using LightGraphs, NamedTupleTools, Random, Statistics, LinearAlgebra
export CausalModel
export rand

# design choice: kappa(x) or kappa(x...)

"Causal model with `dag` representing the model and `scm` holds the name of the variables and the SCM"
struct CausalModel{U, T, S, R} <: AbstractGraph{T}
    dag::SimpleDiGraph{T}  # Representing the DAG of the causal model
    kernel::S  # Name and function (or distribution) of each variable
    topo::R
end
CausalModel(::Type{U}, dag::SimpleDiGraph{T}, kernel::S) where {U,T,S} = CausalModel{U,T,S,Vector{T}}(dag, kernel, topological_sort_by_dfs(dag))

struct FactorModel{T, S<:AbstractVector}
    graph::SimpleGraph{T}
    laws::S
end

kernel_call(kernel, X, v, J) = kernel[v](X[J]...)

import Random.rand
function rand(rng::Random.AbstractRNG, m::CausalModel{T}) where {T}
    topo = m.topo
    dag = m.dag
    kernel = m.kernel
    X = zeros(T, nv(dag))
    for v in topo
        X[v] = rand(rng, kernel_call(kernel, X, v, inneighbors(dag, v)))
    end
    X
end

end


text:
@testset "Yen" begin
    g4 = path_digraph(5)
    d1 = float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])
    d2 = sparse(float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]))

    for g in testdigraphs(g4)
        x = @inferred(yen_k_shortest_paths(g, 5, 5))
        @test length(x.dists) == length(x.paths) ==  1
        @test x.dists[1] ==  0
        @test x.paths[1] ==  [5]

        y = @inferred(yen_k_shortest_paths(g, 2, 5, d1, 1))
        z = @inferred(yen_k_shortest_paths(g, 2, 5, d2, 1))

        @test length(y.dists) == length(z.dists) == 1
        @test length(y.paths) == length(z.paths) == 1
        @test y.dists == z.dists == [33.0]
        @test z.paths == y.paths == [[2, 3, 4, 5]]

        # Remove edge to test empty paths
        rem_edge!(g, 1, 2)
        x = @inferred(yen_k_shortest_paths(g, 1, 5))
        @test length(x.dists) == length(x.paths) ==  0
    end

    gx = path_graph(5)
    add_edge!(gx, 2, 4)
    d = ones(Int, 5, 5)
    d[2, 3] = 100

    for g in testgraphs(gx)
        z = @inferred(yen_k_shortest_paths(g, 1, 5, d, 1))
        @test length(z.dists) == 1
        @test z.dists == [3]
        @test z.paths[1] == [1, 2, 4, 5]
    end

    G = LightGraphs.Graph()
    add_vertices!(G, 4)
    add_edge!(G, 2, 1)
    add_edge!(G, 2, 3)
    add_edge!(G, 1, 4)
    add_edge!(G, 3, 4)
    add_edge!(G, 2, 2)
    w = [0. 3. 0. 1.;
        3. 0. 2. 0.;
        0. 2. 0. 3.;
        1. 0. 3. 0.]

    for g in testgraphs(G)
        ds = @inferred(yen_k_shortest_paths(g, 2, 4, w))
        @test ds.paths == [[2, 1, 4]]
        ds = @inferred(yen_k_shortest_paths(g, 2, 1, w))
        @test ds.paths == [[2, 1]]
        ds = @inferred(yen_k_shortest_paths(g, 2, 3, w))
        @test ds.paths == [[2, 3]]

        # Test with multiple paths
        ds = @inferred(yen_k_shortest_paths(g, 2, 4, w, 2))
        @test ds.paths == [[2, 1, 4], [2, 3, 4]]

        # here a selflink at source is introduced; it should not change the shortest paths
        w[2, 2] = 10.0
        ds = @inferred(yen_k_shortest_paths(g, 2, 4, w))
        @test ds.paths == [[2, 1, 4]]
        ds = @inferred(yen_k_shortest_paths(g, 2, 4, w, 3))
        @test ds.paths == [[2, 1, 4], [2, 3, 4]]

        # Test with new short link
        add_edge!(g, 2, 4)
        w[2, 4] = 1.
        w[4, 2] = 1.
        ds = @inferred(yen_k_shortest_paths(g, 2, 4, w, 2))
        @test ds.paths == [[2, 4], [2, 1, 4]]
        ds = @inferred(yen_k_shortest_paths(g, 2, 4, w, 3))
        @test ds.paths == [[2, 4], [2, 1, 4], [2, 3, 4]]
    end

    # Testing with the example in https://en.wikipedia.org/wiki/Yen%27s_algorithm
    G = LightGraphs.DiGraph()
    add_vertices!(G, 6)
    add_edge!(G, 1, 2)
    add_edge!(G, 1, 3)
    add_edge!(G, 2, 4)
    add_edge!(G, 3, 2)
    add_edge!(G, 3, 4)
    add_edge!(G, 3, 5)
    add_edge!(G, 4, 5)
    add_edge!(G, 4, 6)
    add_edge!(G, 5, 6)
    w = [0. 3. 2. 0. 0. 0.;
        0. 0. 0. 4. 0. 0.;
        0. 1. 0. 2. 3. 0.;
        0. 0. 0. 0. 2. 1.;
        0. 0. 0. 0. 0. 2.;
        0. 0. 0. 0. 0. 0.;]

    for g in testdigraphs(G)
        ds = @inferred(yen_k_shortest_paths(g, 1, 6, w, 3))
        @test ds.dists == [5., 7., 8.]
        @test ds.paths[1] == [1, 3, 4, 6]
        @test ds.paths[2] == [1, 3, 5, 6]
    end

    # Test all the paths
    add_edge!(G, 1, 4)
    w[1, 4] = 3
    w[4, 1] = 3
    for g in testdigraphs(G)
        ds = @inferred(yen_k_shortest_paths(G, 1, 6, w, 100))
        @test ds.dists == [4.0, 5.0, 7.0, 7.0, 8.0, 8.0, 8.0, 11.0, 11.0]

        ds = @inferred(yen_k_shortest_paths(G, 1, 6, w, 100, maxdist=7))
        @test ds.dists == [4.0, 5.0, 7.0, 7.0]
    end
end


text:
"""
    continuousphase(signal, events) -> Vector

Calculate the continuous phase of `signal` by finding the angle of the Hilbert transform of the
amplitude centered signal, `signal`. The amplitude centering uses the average extrema between
pairs of `events`.

This method of continuous phase calculation is as recommended by Lamb and Stöckl (2014) [1]

[1] P. F. Lamb and M. Stöckl, “On the use of continuous relative phase: Review of current
    approaches and outline for a new standard,” Clinical Biomechanics, vol. 29, no. 5,
    pp. 484–493, May 2014, [doi](https://doi.org/10.1016%2Fj.clinbiomech.2014.03.008).
"""
function continuousphase(signal, events::AbstractVector{Int}; centerfun=_center)
    signalcent = centerfun(signal, events)
    Hsignal = hilbert(signalcent)

    return angle.(Hsignal)
end

function _center(signal, events)
    mi, ma = avgextrema(signal, events)
    signal .- Ref(mi + (ma - mi)/2)
end

function crpensemble(sigA, sigB, events; centerfun=_center)
    θA = continuousphase(sigA, events; centerfun)
    θB = continuousphase(sigB, events; centerfun)

    crp = unwrap(θA; range=2pi) - unwrap(θB; range=2pi)

    ens_avg, ens_std = limitcycle(crp, events, mean=circmean, std=circstd)

    return rad2deg.(ens_avg), rad2deg.(ens_std)
end

"""
    demean(x)

Subtract the mean of `x` from `x`
"""
function demean(x)
    return x .- mean(x)
end

"""
    demean!(x)

Mutate `x` by subtracting its mean
"""
function demean!(x)
    x .-= mean(x)
end

function _linreg(x::AbstractVector{T}, y::AbstractVector{U}) where {T,U}
    V = promote_type(T, U)
    _linreg(convert(Vector{V}, collect(x)), convert(Vector{V}, collect(y)))
end

# From @dmbates in https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/28
function _linreg(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:AbstractFloat}
    (N = length(x)) == length(y) || throw(DimensionMismatch())
    ldiv!(cholesky!(Symmetric([T(N) sum(x); zero(T) sum(abs2, x)], :U)), [sum(y), dot(x, y)])
end

"""
    detrend(y)

Remove the linear trend from `y`
"""
function detrend(y)
    β = linreg(1:length(y),y)
    @. return y - (β[1] + β[2]*y)
end

"""
    detrend!(y)

Mutate `y` by removing its linear trend
"""
function detrend!(y)
    a, b = linreg(1:length(y),y)
    @. y -= β[1] + β[2]*y
end

struct ForwardBackwardPad; end

"""
    centraldiff(x, order=1; dt=1, padding=NaN)

Take a numeric finite central difference along the first dimension of the array `x`, where
`dt` is the period between samples in `x`, and `padding` can be `nothing` (returning an
array which is 2 elements shorter than `x`), NaN, or `ForwardBackwardPad()` (which uses a
forward or backward finite difference for the first and last elements, respectively, of the
result).
"""
function centraldiff(x::AbstractVector{T}, order::Integer=1; dt=1, padding=NaN) where T
    if padding === nothing
        x′ = similar(x, T, length(x)-2)
        _centraldiff!(x′, x, Val(order), dt)
    else
        x′ = similar(x)
        _centraldiff!(x′, x, Val(order), dt, padding)
    end

    return x′
end

function centraldiff(
    x::Matrix{T}, order::Integer=1; dt=1, padding=NaN
) where T
    sz = size(x)
    dims=2
    if padding === nothing
        sz = ntuple(i -> i == 1 ? sz[1] - 2 : sz[i], length(sz))
        x′ = similar(x, sz)
        foreach(eachslice(x′; dims), eachslice(x; dims)) do vx′, vx
            _centraldiff!(vx′, vx, Val(order), dt)
        end
    else
        x′ = similar(x, sz)
        foreach(eachslice(x′; dims), eachslice(x; dims)) do vx′, vx
            _centraldiff!(vx′, vx, Val(order), dt, padding)
        end
    end

    return x′
end

@inline function _centraldiff!(x′::AbstractVector, x::AbstractVector, order, dt, padding)
    _centraldiff!(@view(x′[begin+1:end-1]), x, order, dt)
    x′[begin] = padding
    x′[end] = padding

    return nothing
end

@inline function _centraldiff!(x′::AbstractVector, x::AbstractVector, order::Val{1}, dt, padding::ForwardBackwardPad)
    _centraldiff!(@view(x′[begin+1:end-1]), x, order, dt)
    x′[begin] = (x[begin+1] - x[begin]) / dt
    x′[end] = (x[end] - x[end-1]) / dt
end

@inline function _centraldiff!(x′::AbstractVector, x::AbstractVector, order::Val{2}, dt, padding::ForwardBackwardPad)
    _centraldiff!(@view(x′[begin+1:end-1]), x, order, dt)
    x′[begin] = (x[begin+2] - 2*x[begin+1] + x[begin]) / dt^2
    x′[end] = (x[end] - 2*x[end-1] + x[end-2]) / dt^2

    return nothing
end

@inline function _centraldiff!(x′::AbstractVector, x::AbstractVector, order::Val{1}, dt)
    start = firstindex(x)
    last = lastindex(x) - 2
    @boundscheck (start, last) == (firstindex(x′), lastindex(x′)) || throw(BoundsError())

    @inbounds @simd ivdep for i in start:last
        x′[i] = (x[i+2] - x[i]) / 2dt
    end

    return nothing
end

@inline function _centraldiff!(x′::AbstractVector, x::AbstractVector, order::Val{2}, dt)
    start = firstindex(x)
    last = lastindex(x) - 2
    @boundscheck (start, last) == (firstindex(x′), lastindex(x′)) || throw(BoundsError())

    @inbounds @simd ivdep for i in start:last
        x′[i] = (x[i+2] - 2x[i+1] + x[i]) / dt^2
    end

    return nothing
end


text: import Compat: String, walkdir


text:
abstract type AbstractProductMeasure <: AbstractMeasure end

struct ProductMeasure{K,I} <: AbstractProductMeasure
    κ::K
    xs::I
end

marginals(μ::ProductMeasure) = mappedarray(μ.κ, μ.xs)

productmeasure(κ::AbstractKernel, xs::AbstractArray) = ProductMeasure(κ, xs)

function logdensity_def(μ::AbstractProductMeasure, x::AbstractArray)
    mapreduce(+, μ.xs, x) do (j,x)
        logdensity_def(μ.κ(j), x)
    end
end


text:
#
# - (2017.07.22) Copy of sim1.jl
#
# - IJCNN:
#
#   V.0.3 (27.10.2016) No ROC yet
#
# ROC curve for artificial signal.
# Simulation - generate changes and detect them.
# Obtain ROC curve.

include("./lib/metadata.jl")
include("./lib/pccf.jl")
include("./1naive_update.jl")
include("./2bayes_detector_replica.jl")
include("./lib/SignalGenerator.jl")
include("./lib/roc.jl")

using SignalGenerator, BayesChpOrig, Roc, Metadata, Pccf
--todo

"""
In this simulation we variate the lambda parameter
of the Bayesian detector and plot ROC curve.
"""
function run_sim1(outputfolder)
    println("... sim1.jl.run_sim1 > save results into ", outputfolder)
    s = Dict()
    #--- START SIMULATION SETTINGS
    max_delay        = 10  # 30
    n_sim_for_each_h = 200 # 10
    h_from           = 50  # 40
    h_to             = 300 # 300
    h_by             = 15  # 25  # 15
    #--- end SIMULATION SETTINGS

    avg_tp_rates_for_h = []
    avg_fa_rates_for_h = []
    avg_ppv_for_h      = []

    pccf_avg_tp_rates_for_h = []
    pccf_avg_fa_rates_for_h = []
    pccf_avg_ppv_for_h      = []

    mu0        = 100.0
    sd0        = 10.0
    n_chps     = 10
    chpStep    = 1.3
    threshold  = 0.3
    max_d      = 500 # 400

    n_pccf = round(Int64, mu0 * (n_chps + 5)) # Length of the prediction using PCCF (so that to not recalc on eachj iteration)
    probs  = scale01(pccf(n_pccf, mu0, sd0))

    ##?? srand(1235) - we need a new data on every iteration

    # Variate lambda value

    for h in h_from : h_by : h_to
        if mod(h, 10)==0
            @printf("%0.2f %s \n", h/h_to * 100, "%")
        end

        l_common = h

        tp_rates_for_i = []
        fa_rates_for_i = []
        ppv_for_i      = []

        pccf_tp_rates_for_i = []
        pccf_fa_rates_for_i = []
        pccf_ppv_for_i      = []

        for i_signals = 1 : n_sim_for_each_h
            # # genSig(n, mu, sigma, changeStep, sigStd)
            # helped : variate changeStep
            # ret = genSig(5, 100.0, 2.0, 1.5, 0.75, 1)
            # ret = genSig(10, 100.0, 10.0, 1.2, 0.75, 1) # seeems interesting
            ret          = genSig(n_chps, mu0, sd0, chpStep, 0.75)   # OK
            sig          = ret["sig"]
            changes      = ret["changes"]
            n_sig        = length(sig)

            ADDNOISE = true

            if ADDNOISE
                #: addNoise!(sig, changes, changeStep; width = 5)
                noise_ids = addNoise!(sig, changes, chpStep, width = 10)
                # -- stopped here (06 nov 2016)
                # inc mu ? to add more noise
                #noise_ids = addNoiseRec!(sig, changes, chpStep, width = 7, min_len = 50)
            end

            if n_sig > n_pccf
                # Update PCCF if its length is less than signal's length
                n_pccf = n_sig
                probs  = scale01(pccf(n_pccf, mu0, sd0))
            end

            maxes, mat   = bayes_detector(sig, lambda = h)
            detections   = extract_changes(maxes)

            # --- check max_d parameter value here
            pccf_maxes, pccf_mat = pccf_bayes_detector(sig, probs[1:n_sig], l_common + 100, l_common, threshold, h_common = l_common, max_d = max_d)
            pccf_detections      = extract_changes(pccf_maxes)

            frac = 0.1/1.0 # importance factor of TN's

            tpr          = tpRate(detections, changes, max_delay)
            fpr          = faRateScaled(detections, changes, max_delay, n_sig, tn_price_scale = frac)
            ppv_val      = ppv(detections, changes, max_delay)

            pccf_tpr     = tpRate(pccf_detections, changes, max_delay)
            pccf_fpr     = faRateScaled(pccf_detections, changes, max_delay, n_sig, tn_price_scale = frac)
            pccf_ppv_val = ppv(pccf_detections, changes, max_delay)

            push!(ppv_for_i     , ppv_val)
            push!(tp_rates_for_i, tpr)
            push!(fa_rates_for_i, fpr)

            push!(pccf_ppv_for_i     , pccf_ppv_val)
            push!(pccf_tp_rates_for_i, pccf_tpr)
            push!(pccf_fa_rates_for_i, pccf_fpr)

        end
        push!(avg_tp_rates_for_h, mean(tp_rates_for_i))
        push!(avg_fa_rates_for_h, mean(fa_rates_for_i))
        push!(avg_ppv_for_h     , mean(ppv_for_i))

        push!(pccf_avg_tp_rates_for_h, mean(pccf_tp_rates_for_i))
        push!(pccf_avg_fa_rates_for_h, mean(pccf_fa_rates_for_i))
        push!(pccf_avg_ppv_for_h     , mean(pccf_ppv_for_i))
    end
    ## OLD
    # sim1_tp_rates.txt, sim1_fa_rates.txt,...
    writedlm(outputfolder * "tp_rates.txt", avg_tp_rates_for_h)
    writedlm(outputfolder * "fa_rates.txt", avg_fa_rates_for_h)
    writedlm(outputfolder * "ppv.txt"     , avg_ppv_for_h)

    writedlm(outputfolder * "pccf_tp_rates.txt", pccf_avg_tp_rates_for_h)
    writedlm(outputfolder * "pccf_fa_rates.txt", pccf_avg_fa_rates_for_h)
    writedlm(outputfolder * "pccf_ppv.txt"     , pccf_avg_ppv_for_h)
end
run_sim1(outdir * "sim1/")

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
# This file was generated, do not modify it. # hide
using MLJ
import DataFrames: DataFrame
using PrettyPrinting
using StableRNGs
MLJ.color_off() # hide

rng = StableRNG(512)
Xraw = rand(rng, 300, 3)
y = exp.(Xraw[:,1] - Xraw[:,2] - 2Xraw[:,3] + 0.1*rand(rng, 300))
X = DataFrame(Xraw, :auto)

train, test = partition(eachindex(y), 0.7);

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
@testset "989.add-add_to_array_form-of-integer.jl" begin
    @test add_to_array_form([1, 2, 0, 0], 34) == [1, 2, 3, 4]
    @test add_to_array_form([2, 7, 4], 181) == [4, 5, 5]
    @test add_to_array_form([2, 1, 5], 806) == [1, 0, 2, 1]
    @test add_to_array_form([9, 9, 9, 9, 9, 9, 9, 9, 9, 9], 1) == [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
end


text:
using IRTools
using IRTools: IR, slots!, argnames!, varargs!, xcall, SSAValue
using ..Zygote: Argument, argmap, ignored, isexpr

function record!(ir::IR)
  ir = map(ex -> argmap(a -> Argument(a.n+2), ex), ir)
  grads = Dict()
  grad(x::Union{SSAValue,Argument}) = grads[x]
  grad(x) = nothing
  ks = keys(ir)
  for n = 1:length(ir.args)
    grads[Argument(n+2)] = pushfirst!(ir, xcall(Base, :getindex, Argument(2)))
  end
  for k in ks
    ex = ir[k].expr
    if isexpr(ex, :call) && !ignored(ir, ex)
      Δ = insert!(ir, k, xcall(Base, :tuple, grad.(ex.args)...))
      ir[k] = xcall(Forward, :_forward, Δ, ex.args...)
    end
  end
  return ir
end

@generated function _forward(Δ::Tuple, f, args...)
  T = Tuple{f,args...}
  meta = IRTools.meta(T)
  forw = record(IR(meta))
  pushfirst!(forw.args, Any, Any)
  argnames!(meta, Symbol("#self#"), :Δ, :f, :args)
  forw = varargs!(meta, forw, 3)
  forw = slots!(forw)
  return IRTools.update!(meta, forw)
end


text:
# This file implements the loss functions
# For now it only includes the generator loss
using Knet
using Statistics

function HingeAdversarialGenerator(fake)
    return mean(-1.*fake)
end

function AdversarialGenerator(fake)
    return mean(softplus(-1.*fake))
end

function softplus(x;beta=1)
    1/beta * log.(1 + exp.(beta.*x))
end


text:
#=
state_space.jl

    State space model abstract types and general fallback routines

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2022/02/08
=#

"""
    nparams(model)

Determine the number of hyper parameters.

#### Arguments
  - `model::StateSpaceModel`: state space model

#### Returns
  - `n_params::Integer` : number of hyper parameters
"""
function nparams end

"""
	get_params!(ψ, model)

Retrieve hyper parameters from the state space model and store them in `ψ`.

#### Arguments
  - `model::StateSpaceModel`: state space model

#### Returns
  - `ψ::AbstractVector`     : hyper parameters
"""
function get_params! end

"""
    store_params!(model, ψ)

Store hyper parameters `ψ` from the state space model in `model`.

#### Arguments
  - `ψ::AbstractVector`     : hyper parameters

#### Returns
  - `model::StateSpaceModel`: state space model
"""
function store_params! end

"""
    get_system!(sys, model, method)

Retrieve system matrices from the state space model and store them in `sys`.

#### Arguments
  - `model::StateSpaceModel`: state space model
  - `method::Symbol`		: filtering method

#### Returns
  - `sys::StateSpaceSystem` : state space system matrices
"""
function get_system! end

"""
    create_system(model, method)

Create system matrices from the state space model.

#### Arguments
  - `model::StateSpaceModel`: state space model
  - `method::Symbol`		: filtering method

#### Returns
  - `sys::StateSpaceSystem` : state space system matrices
"""
function create_system! end

"""
    init!(model, init, method)

Initialize the state space model hyper parameters as defined by `model` and the
initial conditions of the state space model in `sys`.

#### Arguments
  - `init::NamedTuple`	: initial hyper parameters
  - `method::Symbol`	: filtering method

#### Returns
  - `model::StateSpaceModel`: state space model
  - `sys::StateSpaceSystem` : state space system matrices
"""
function init! end

"""
    init_model!(model, init)

Initialize the state space model hyper parameters.

#### Arguments
  - `init::NamedTuple`	: initial hyper parameters

#### Returns
  - `model::StateSpaceModel`: state space model
"""
function init_model! end

"""
    fix_system!(sys, model, method)

Fix state space system components.

#### Arguments
  - `model::StateSpaceModel`: state space model
  - `method::Symbol`		: filtering method

#### Returns
  - `sys::StateSpaceSystem` : state space system matrices
"""
function fix_system! end

"""
    init_system!(sys, model)

Initialize state space system.

#### Arguments
  - `model::StateSpaceModel`: state space model

#### Returns
  - `sys::StateSpaceSystem` : state space system matrices
"""
function init_system! end

"""
	loglik(filter, sys, model, method)

Compute the log-likelihood for a linear Gaussian State Space model, given by
`model`, with Kalman filter output `filter`, based on filtering method
(`method`).

#### Arguments
  - `filter::KalmanFilter`	: Kalman filter output
  - `sys::StateSpaceSystem` : state space system matrices
  - `model::StateSpaceModel`: state space model
  - `method::Symbol`		: filtering method

#### Returns
  - `ll::Real`	: log-likelihood
"""
function loglik end


text:
# Metagraphs files are simply JLD2 files.

"""
    struct MGFormat <: AbstractGraphFormat end

You can save `MetaGraph`s in a `MGFormat`, currently based on `JLD2`.

```jldoctest
julia> using MetaGraphsNext

julia> using LightGraphs: Edge, Graph,  loadgraph, savegraph

julia> example = MetaGraph(Graph());

julia> mktemp() do file, io
            savegraph(file, example)
            loadgraph(file, "something", MGFormat()) == example
        end
true
```
"""
struct MGFormat <: AbstractGraphFormat end

"""
    struct DOTFormat <: AbstractGraphFormat end

If all metadata types support `pairs` or are `nothing`, you can save `MetaGraph`s in `DOTFormat`.

```jldoctest DotFormat
julia> using MetaGraphsNext

julia> using LightGraphs

julia> simple = MetaGraph(Graph());

julia> simple[:a] = nothing; simple[:b] = nothing; simple[:a, :b] = nothing;

julia> mktemp() do file, io
            savegraph(file, simple, DOTFormat())
            print(read(file, String))
        end
graph T {
    a
    b
    a -- b
}

julia> complicated = MetaGraph(DiGraph(),
            VertexMeta = Dict{Symbol, Int},
            EdgeMeta = Dict{Symbol, Int},
            gprops = (tagged = true,)
        );

julia> complicated[:a] = Dict(:code_1 => 1, :code_2 => 2);

julia> complicated[:b] = Dict(:code => 2);

julia> complicated[:a, :b] = Dict(:code => 12);

julia> mktemp() do file, io
            savegraph(file, complicated, DOTFormat())
            print(read(file, String))
        end
digraph G {
    tagged = true
    a [code_1 = 1, code_2 = 2]
    b [code = 2]
    a -> b [code = 12]
}
```
"""
struct DOTFormat <: AbstractGraphFormat end

function loadmg(fn::AbstractString)
    @load fn g
    return g
end

function savemg(fn::AbstractString, g::MetaGraph)
    @save fn g
    return 1
end

loadgraph(fn::AbstractString, gname::String, ::MGFormat) = loadmg(fn)
savegraph(fn::AbstractString, g::MetaGraph) = savemg(fn, g)

function show_meta_list(io::IO, meta)
    if meta !== nothing && length(meta) > 0
        next = false
        write(io, " [")
        for (key, value) in meta
            if next
                write(io, ", ")
            else
                next = true
            end
            write(io, key)
            write(io, " = ")
            show(io, value)
        end
        write(io, "]")
    end
end

function savedot(io::IO, g::MetaGraph)
    gprops = g.gprops
    metaindex = g.metaindex

    if is_directed(g)
        write(io, "digraph G {\n")
        dash = "->"
    else
        write(io, "graph T {\n")
        dash = "--"
    end

    if gprops !== nothing
        for (key, value) in pairs(gprops)
            write(io, "    ")
            write(io, key)
            write(io, " = ")
            show(io, value)
            write(io, '\n')
        end
    end

    for label in keys(g.vprops)
        write(io, "    ")
        write(io, label)
        show_meta_list(io, g[label])
        write(io, '\n')
    end

    for (label_1, label_2) in keys(g.eprops)
        write(io, "    ")
        write(io, label_1)
        write(io, ' ')
        write(io, dash)
        write(io, ' ')
        write(io, label_2)
        show_meta_list(io, g.eprops[arrange(g, label_1, label_2)])
        write(io, "\n")
    end
    write(io, "}\n")
end

function savegraph(fn::AbstractString, g::MetaGraph, ::DOTFormat)
    open(fn, "w") do fp
        savedot(fp, g)
    end
end


text:
using Base.Test
using WoodburyMatrices
using Compat:view

srand(123)
n = 5

for elty in (Float32, Float64, Complex64, Complex128, Int), AMat in (diagm,)

    elty = Float64

    a = rand(n); B = rand(n,2); D = rand(2,2); v = rand(n)

    if elty == Int
        v = rand(1:100, n)
        a = rand(1:100, n)
        B = rand(1:100, 2, n)
        D = rand(1:100, 2, 2)
    else
        v = convert(Vector{elty}, v)
        a = convert(Vector{elty}, a)
        B = convert(Matrix{elty}, B)
        D = convert(Matrix{elty}, D)
    end

    ε = eps(abs2(float(one(elty))))
    A = AMat(a)

    # Woodbury
    for W in (SymWoodbury(A, B, D), SymWoodbury(A, B[:,1][:], 2.))

        F = full(W)

        @test_approx_eq (2*W)*v 2*(W*v)
        @test_approx_eq W'*v W*v
        @test_approx_eq (W'W)*v full(W)*(full(W)*v)
        @test_approx_eq (W*W)*v full(W)*(full(W)*v)
        @test_approx_eq (W*W')*v full(W)*(full(W)*v)
        @test_approx_eq W[1:3,1:3]*v[1:3] full(W)[1:3,1:3]*v[1:3]
        @test_approx_eq sparse(W) full(W)
        @test W === W'
        @test_approx_eq W*eye(n) full(W)
        @test_approx_eq W'*eye(n) full(W)

        Z = randn(n,n)
        @test_approx_eq full(W*Z) full(W)*Z

        R = rand(n,n)

        for v = (rand(n, 1), view(rand(n,1), 1:n), view(rand(n,2),1:n,1:2))
            @test_approx_eq (2*W)*v 2*(W*v)
            @test_approx_eq (W*2)*v 2*(W*v)
            @test_approx_eq (W'W)*v full(W)*(full(W)*v)
            @test_approx_eq (W*W)*v full(W)*(full(W)*v)
            @test_approx_eq (W*W')*v full(W)*(full(W)*v)
            @test_approx_eq W[1:3,1:3]*v[1:3] full(W)[1:3,1:3]*v[1:3]
            @test_approx_eq full(WoodburyMatrices.conjm(W, R)) R*full(W)*R'
            @test_approx_eq full((copy(W)'W)*v) full(W)*(full(W)*v)
            @test_approx_eq full(W + A) full(W)+full(A)
            @test_approx_eq full(A + W) full(W)+full(A)
        end

        v = rand(n,1)

        W2 = convert(Woodbury, W)
        @test_approx_eq full(W2) full(W)

        if elty != Int
            @test_approx_eq inv(W)*v inv(full(W))*v
            @test_approx_eq W\v inv(full(W))*v
            @test_approx_eq liftFactor(W)(v) inv(W)*v
            @test_approx_eq WoodburyMatrices.partialInv(W)[1] inv(W).B
            @test_approx_eq WoodburyMatrices.partialInv(W)[2] inv(W).D
            @test_approx_eq det(W) det(full(W))
        end

    end

end

for elty in (Float32, Float64, Complex64, Complex128, Int)

    elty = Float64

    a1 = rand(n); B1 = rand(n,2); D1 = rand(2,2); v = rand(n)
    a2 = rand(n); B2 = rand(n,2); D2 = rand(2,2);

    if elty == Int
        v = rand(1:100, n)

        a1 = rand(1:100, n)
        B1 = rand(1:100, 2, n)
        D1 = rand(1:100, 2, 2)

        a2 = rand(1:100, n)
        B2 = rand(1:100, 2, n)
        D2 = rand(1:100, 2, 2)
    else
        v = convert(Vector{elty}, v)

        a1 = convert(Vector{elty}, a1)
        B1 = convert(Matrix{elty}, B1)
        D1 = convert(Matrix{elty}, D1)

        a2 = convert(Vector{elty}, a2)
        B2 = convert(Matrix{elty}, B2)
        D2 = convert(Matrix{elty}, D2)
    end

    ε = eps(abs2(float(one(elty))))

    # Woodbury
    A1 = diagm(a1)
    A2 = diagm(a2)

    W1 = SymWoodbury(A1, B1, D1)
    W2 = SymWoodbury(A2, B2, D2)

    W1r = SymWoodbury(A1, B1[:,1][:], 2.)
    W2r = SymWoodbury(A2, B2[:,1][:], 3.)

    for (W1, W2) = ((W1,W2), (W1r, W2), (W1, W2r), (W1r,W2r))
        @test_approx_eq (W1 + W2)*v (full(W1) + full(W2))*v
        @test_approx_eq (full(W1) + W2)*v (full(W1) + full(W2))*v
        @test_approx_eq (W1 + 2*diagm(a1))*v (full(W1) + full(2*diagm(a1)))*v
        @test_throws MethodError W1*W2
    end

end

# Sparse U and D

A = diagm(rand(n))
B = sprandn(n,2,1.)
D = sprandn(2,2,1.)
W = SymWoodbury(A, B, D)
v = randn(n)
V = randn(n,1)

@test size(W) == (n,n)
@test size(W,1) == n
@test size(W,2) == n

@test_approx_eq inv(W)*v inv(full(W))*v
@test_approx_eq (2*W)*v 2*(W*v)
@test_approx_eq (W'W)*v full(W)*(full(W)*v)
@test_approx_eq (W*W)*v full(W)*(full(W)*v)
@test_approx_eq (W*W')*v full(W)*(full(W)*v)
@test_approx_eq liftFactor(W)(v) inv(W)*v

@test_approx_eq inv(W)*V inv(full(W))*V
@test_approx_eq (2*W)*V 2*(W*V)
@test_approx_eq (W'W)*V full(W)*(full(W)*V)
@test_approx_eq (W*W)*V full(W)*(full(W)*V)
@test_approx_eq (W*W')*V full(W)*(full(W)*V)

# Mismatched sizes
@test_throws DimensionMismatch SymWoodbury(rand(5,5),rand(5,2),rand(2,3))
@test_throws DimensionMismatch SymWoodbury(rand(5,5),rand(5,2),rand(3,3))
@test_throws DimensionMismatch SymWoodbury(rand(5,5),rand(3),1.)


text:
using PyCall

py"""
import numpy as np

def carregar_npz(arquivo):
    ks = np.load(arquivo)
    return [ks[ks.files[i]] for i in range(len(ks.files))]

def carregar_npy(arquivo):
    return np.load(arquivo)
"""

# using NPZ
#
# function carregar_npz(arquivo)
#     arq = npzread(arquivo)
#     arq_def = [v[2] for v in arq]
# end
#
# function carregar_npy(arquivo)
#     npzread(arquivo)
# end

ctes = Dict(1 => "kelems",
            2 => "gls_elementos",
            3 => "gls_estrutura",
            4 => "apoios",
            5 => "forcas",
            6 => "RCM",
            7 => "volumes_elementos",
            8 => "nos")

"""
Corrige os números dos graus de liberdade da notação do Python para a notação
de Julia. Soma 1 aos graus de liberdade não impedidos.
"""
function corrigir_graus_lib(x)
    if x != -1
        x + 1
    else
        x
    end
end

"""
Retorna um dicionário com as leituras dos arquivos necessários para a resolução do sistema
de equações. Os arquivos já devem ter sido extraídos para que a leitura seja efetiva.
"""
function ler_arquivos_entrada(arq_kelems, arq_gls_elementos, arq_gls_estrutura,
    arq_apoios, arq_forcas, arq_rcm)
    # , arq_volumes_elementos, arq_nos

    dados = Dict(
        ctes[1] => py"carregar_npz"(arq_kelems),
        ctes[2] => py"carregar_npz"(arq_gls_elementos),
        ctes[3] => py"carregar_npy"(arq_gls_estrutura),
        ctes[4] => py"carregar_npy"(arq_apoios),
        ctes[5] => py"carregar_npy"(arq_forcas),
        ctes[6] => py"carregar_npy"(arq_rcm)
    )

    # dados = Dict(
    #     ctes[1] => carregar_npz(arq_kelems),
    #     ctes[2] => carregar_npz(arq_gls_elementos),
    #     ctes[3] => carregar_npy(arq_gls_estrutura),
    #     ctes[4] => carregar_npy(arq_apoios),
    #     ctes[5] => carregar_npy(arq_forcas),
    #     ctes[6] => carregar_npy(arq_rcm),
    #     ctes[7] => carregar_npy(arq_volumes_elementos),
    #     ctes[8] => carregar_npy(arq_nos)
    # )

    # A todos os graus de liberdade deve ser somado 1 para a compatibilização com Julia
    for (i, gl) in enumerate(dados[ctes[2]])
        dados[ctes[2]][i] = map(corrigir_graus_lib, gl)
    end

    # gls_estrutura
    dados[ctes[3]] = map(corrigir_graus_lib, dados[ctes[3]])
    # apoios
    dados[ctes[4]] = map(corrigir_graus_lib, dados[ctes[4]])
    # RCM
    dados[ctes[6]] = map(corrigir_graus_lib, dados[ctes[6]])

    # Return
    dados
end

max_stars_repo_name: lorranfoliveira/otm

max_issues_repo_name: lorranfoliveira/otm

max_forks_repo_name: lorranfoliveira/otm


text:
#
module Hymod

using DataFrames, Dates, Distributions, TOML

include("Utils.jl")
include("Calibration.jl")

    function get_random_params()
        param_bounds = Dict{Symbol,Dict}(
            :cmax => Dict{Symbol,Float64}(:lower => 1.0, :upper => 100),
            :bexp => Dict{Symbol,Float64}(:lower => 0.0, :upper => 2.0),
            :alpha => Dict{Symbol,Float64}(:lower => 0.2, :upper => 0.99),
            :ks => Dict{Symbol,Float64}(:lower => 0.01, :upper => 0.5),
            :kq => Dict{Symbol,Float64}(:lower => 0.5, :upper => 1.2)
        )

        out = Dict()

        for k in keys(param_bounds)
            minV = param_bounds[k][:lower]
            maxV = param_bounds[k][:upper]
            p = rand(Uniform(minV,maxV), 1)
            out[k] = p[1]
        end

        return out
    end

    function _power(x,y)
        x=abs(x) # Needed to capture invalid overflow with netgative values
        return x^y
    end

    function _excess(x_loss,cmax,bexp,prc,pet)
        xn_prev = x_loss
        ct_prev = cmax * (1 - _power((1 - ((bexp + 1) * (xn_prev) / cmax)), (1 / (bexp + 1))))
        # Calculate Effective rainfall 1
        ER1 = max((prc - cmax + ct_prev), 0.0)
        prc = prc- ER1
        dummy = min(((ct_prev + prc) / cmax), 1)
        xn = (cmax / (bexp + 1)) * (1 - _power((1 - dummy), (bexp + 1)))

        # Calculate Effective rainfall 2
        ER2 = max(prc - (xn - xn_prev), 0)

        # Alternative approach
        evap = (1 - (((cmax / (bexp + 1)) - xn) / (cmax / (bexp + 1)))) * pet  # actual ET is linearly related to the soil moisture state
        xn = max(xn - evap, 0)  # update state

        return ER1,ER2,xn
    end

    function linearReservoir(x_slow,inflow,Rs)
        # Linear reservoir routing
        x_slow = (1 - Rs) * x_slow + (1 - Rs) * inflow
        outflow = (Rs / (1 - Rs)) * x_slow
        return x_slow,outflow
    end

    function hargreaves(forcings; tminCol=:tmin,tmaxCol=:tmax,dtCol=:datetime)

        dts = forcings[!,dtCol]
        tmin = forcings[!,tminCol]
        tmax = forcings[!,tmaxCol]
        len = length(tmax)
        Gsc = 367
        lhov = 2.257

        doy = map(dayofyear, dts)

        tavg = map(mean,zip(tmin,tmax))

        eto = zeros(len)

        for (i,t) in enumerate(doy)
            b = 2 * pi * (Int16(t)/365)
            Rav = 1.00011 + 0.034221*cos(b) + 0.00128*sin(b) + 0.000719*cos(2*b) + 0.000077*sin(2*b)
            Ho = ((Gsc * Rav) * 86400)/1e6

            eto[i] = (0.0023 * Ho * (tmax[i]-tmin[i])^0.5 * (tavg[i]+17.8))
        end

        return eto
    end

    function simulate(forcings; precipCol=:precip,petCol=:pet,initFlow=true,cmax=1.0,bexp=0.0,alpha=0.2,ks=0.01,kq=0.5)
        p = forcings[!,precipCol]
        e = forcings[!,petCol]

        lt_to_m = 0.001

        # HYMOD PROGRAM IS SIMPLE RAINFALL RUNOFF MODEL
        x_loss = 0.0
        # Initialize slow tank state
        # value of 0 init flow works ok if calibration data starts with low discharge
        x_slow = initFlow ? (2.3503 / (ks * 22.5)) : 0
        # Initialize state(s) of quick tank(s)
        x_quick = zeros(Float64, 3)
        outflow = zeros(Float64, length(p))
        output = zeros(Float64, length(p))
        # START PROGRAMMING LOOP WITH DETERMINING RAINFALL - RUNOFF AMOUNTS

        for t in 1:length(p)
            Pval = p[t]
            PETval = e[t]
            # Compute excess precipitation and evaporation
            ER1, ER2, x_loss = _excess(x_loss, cmax, bexp, Pval, PETval)
            # Calculate total effective rainfall
            ET = ER1 + ER2
            #  Now partition ER between quick and slow flow reservoirs
            UQ = alpha * ET
            US = (1 - alpha) * ET
            # Route slow flow component with single linear reservoir
            x_slow, QS = linearReservoir(x_slow, US, ks)
            # Route quick flow component with linear reservoirs
            inflow = UQ

            for i = 1:length(x_quick)
                # Linear reservoir
                x_quick[i], outflow = linearReservoir(x_quick[i], inflow, kq)
                inflow = outflow
            end

            # Compute total flow for timestep
            output[t] = ((QS + outflow)/lt_to_m)
        end

        output[output.<0] .= 0

        return output
    end

    function calibrate(forcing, obs, paramSpace, samples; precipCol=:precip, petCol=:pet, saveResults=false)
        # implementation of a monte carlo sampling

        params = Calibration.monteCarlo(paramSpace,samples=samples)
        keyList = collect(keys(params))
        nIter = samples
        losses = zeros(nIter)

        print("Running $nIter iterations...")
        for i in 1:nIter
            vals = [params[k][i] for k in keyList]
            pars = Dict(keyList .=> vals)
            q = simulate(forcing, precipCol=:precip, petCol=:pet; pars...)
            losses[i] = Utils.nse(q[365:end],obs[365:end])
        end

        losses = convert(Vector{Union{Missing,Float64}}, losses)
        losses[map(isnan,losses)] .= missing

        finalLoss, idx = findmax(skipmissing(losses))
        vals = [params[k][idx] for k in keyList]
        finalPars = Dict(keyList .=> vals)
        finalQ = simulate(forcing, precipCol=:precip, petCol=:pet; finalPars...)

        if saveResults
            t = now()
            finalPars[:loss] = finalLoss
            tomlPars = Dict(String(k)=>v for (k,v) in finalPars)
            print(tomlPars)
            open("hymod_calibration_results_$t.toml","w") do io
                TOML.print(io, tomlPars)
            end
        end

        return finalQ, finalPars, finalLoss

    end

# end module
end

max_stars_repo_name: jsta/hymod

max_issues_repo_name: jsta/hymod

max_forks_repo_name: jsta/hymod


text:
using Documenter, Zalgo, Luxor

makedocs(
    modules = [Zalgo],
    sitename = "Zalgo",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = ["assets/zalgo-docs.css"]),
    pages = Any[
        "Introduction"  => "index.md",
        "Index"         => "functionindex.md"
    ]
)

deploydocs(
    repo = "github.com/cormullion/Zalgo.jl.git",
    target = "build"
)


text:
import Base.rand

export Singleton

"""
    Singleton{N, VN<:AbstractVector{N}} <: AbstractSingleton{N}

Type that represents a singleton, that is, a set with a unique element.

### Fields

- `element` -- the only element of the set
"""
struct Singleton{N, VN<:AbstractVector{N}} <: AbstractSingleton{N}
    element::VN
end

isoperationtype(::Type{<:Singleton}) = false
isconvextype(::Type{<:Singleton}) = true

# --- AbstractSingleton interface functions ---

"""
    element(S::Singleton)

Return the element of a singleton.

### Input

- `S` -- singleton

### Output

The element of the singleton.
"""
function element(S::Singleton)
    return S.element
end

"""
    element(S::Singleton, i::Int)

Return the i-th entry of the element of a singleton.

### Input

- `S` -- singleton
- `i` -- dimension

### Output

The i-th entry of the element of the singleton.
"""
function element(S::Singleton, i::Int)
    return S.element[i]
end

# --- LazySet interface functions ---

"""
    rand(::Type{Singleton}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)

Create a random singleton.

### Input

- `Singleton` -- type for dispatch
- `N`         -- (optional, default: `Float64`) numeric type
- `dim`       -- (optional, default: 2) dimension
- `rng`       -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`      -- (optional, default: `nothing`) seed for reseeding

### Output

A random singleton.

### Algorithm

The element is a normally distributed vector with entries of mean 0 and standard
deviation 1.
"""
function rand(::Type{Singleton};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing)
    rng = reseed(rng, seed)
    element = randn(rng, N, dim)
    return Singleton(element)
end

"""
    translate(S::Singleton, v::AbstractVector)

Translate (i.e., shift) a singleton by a given vector.

### Input

- `S` -- singleton
- `v` -- translation vector

### Output

A translated singleton.

### Algorithm

We add the vector to the point in the singleton.
"""
function translate(S::Singleton, v::AbstractVector)
    @assert length(v) == dim(S) "cannot translate a $(dim(S))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    return Singleton(element(S) + v)
end

"""
    rectify(S::Singleton)

Concrete rectification of a singleton.

### Input

- `S` -- singleton

### Output

The `Singleton` that corresponds to the rectification of `S`.
"""
function rectify(S::Singleton)
    return Singleton(rectify(element(S)))
end

"""
    project(S::Singleton, block::AbstractVector{Int}; kwargs...)

Concrete projection of a singleton.

### Input

- `S`     -- singleton
- `block` -- block structure, a vector with the dimensions of interest

### Output

A set representing the projection of the singleton `S` on the dimensions
specified by `block`.
"""
function project(S::Singleton, block::AbstractVector{Int}; kwargs...)
    return Singleton(element(S)[block])
end

function permute(S::Singleton, p::AbstractVector{Int})
    e = element(S)[p]
    return Singleton(e)
end