text:
module mcsys_0
using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    XYZ = [0.0 0.0 0.0]
    J0 =  [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    labl = 1
    CM = CompositeLayupModule
    cs = CM.cartesian_csys((1, 2, 3))
    updatecsmat!(cs, reshape(XYZ, 1, 3), J0, labl);
    @test norm(cs.csmat -
        [1.0 0.0 0.0;
        0.0 1.0 0.0;
        0.0 0.0 1.0]) < 1.0e-15
    cs = CM.cartesian_csys((2, -1, 3))
    updatecsmat!(cs, reshape(XYZ, 1, 3), J0, labl);
    @test norm(cs.csmat -
            [0.0 -1.0 0.0;
            1.0 0.0 0.0;
            0.0 0.0 1.0]) < 1.0e-15
    cs = CM.cartesian_csys((2, 1, -3))
    updatecsmat!(cs, reshape(XYZ, 1, 3), J0, labl);
    @test norm(cs.csmat -
            [0.0 1.0 0.0;
            1.0 0.0 0.0;
            0.0 0.0 -1.0]) < 1.0e-15
    true
end
end
using .mcsys_0
mcsys_0.test()

module mlayup1
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    E = 200e3*phun("MPa")
    nu = 0.3;
    rho= 8000*phun("KG/M^3");
    thickness = 0.005*phun("m");
    angle = 45.0
    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
    cm = CompositeLayupModule
    ply1 = cm.Ply("ply1", mater, thickness, angle)
    @test ply1.angle == 45
    true
end
end
using .mlayup1
mlayup1.test()

module mlayup2
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using Test

function test()
    E = 200e3*phun("MPa")
    nu = 0.3;
    rho= 8000*phun("KG/M^3");
    thickness = 0.005*phun("m");
    angle = 45.0
    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
    CM = CompositeLayupModule
    ply1 = CM.Ply("ply1", mater, thickness, 45.0)
    ply2 = CM.Ply("ply2", mater, thickness, -45.0)
    cl = CM.CompositeLayup("sample", [ply1, ply2], CM.cartesian_csys((1, 2, 3))  )
    @test cl.plies[2].angle == -45
    true
end
end
using .mlayup2
mlayup2.test()

module mlayup3
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    mtol = 1.0e-15
    for dangle in [5 34 68 73 535 0 90 45 -135. -87]
        angle = dangle/180*pi
        Tm = fill(0.0, 3, 3)
        CM.plane_stress_T_matrix!(Tm, angle)
        Tinvm = fill(0.0, 3, 3)
        CM.plane_stress_Tinv_matrix!(Tinvm, angle)
        @test norm(Tm*Tinvm - I) < mtol
    end
    true
end
end
using .mlayup3
mlayup3.test()

module mlayup4
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    mtol = 1.0e-15
    R = [1 0 0; 0 1 0; 0 0 2]
    Tm = fill(0.0, 3, 3)
    Tme = fill(0.0, 3, 3)
    for dangle in [5 34 68 73 535 0 90 45 -135. -87]
        angle = dangle/180*pi
        CM.plane_stress_T_matrix!(Tm, angle)
        CM.plane_stress_T_matrix_eng!(Tme, angle)
        # @show R * Tm / R, Tme
        @test norm(R * Tm / R - Tme) < mtol
        # Tinvm = fill(0.0, 3, 3)
        # CM.plane_stress_Tinv_matrix!(DeforModelRedModule.DeforModelRed2DStress, Tinvm, angle)
        # @test norm(Tm*Tinvm - I) < mtol
    end
    true
end
end
using .mlayup4
mlayup4.test()

module mlayup5
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    mtol = 1.0e-15
    R = [1 0 0; 0 1 0; 0 0 2]
    Tm = fill(0.0, 3, 3)
    Tme = fill(0.0, 3, 3)
    for dangle in [5 34 68 73 535 0 90 45 -135. -87]
        angle = dangle/180*pi
        CM.plane_stress_Tinv_matrix!(Tm, angle)
        CM.plane_stress_Tinv_matrix_eng!(Tme, angle)
        @test norm(R * Tm / R - Tme) < mtol
    end
    true
end
end
using .mlayup5
mlayup5.test()

module mlayup6
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    mtol = 1.0e-15
    R = [1 0 0; 0 1 0; 0 0 2]
    Tm = fill(0.0, 3, 3)
    Tinvm = fill(0.0, 3, 3)
    for dangle in [-8.]
        angle = dangle/180*pi
        CM.plane_stress_T_matrix_eng!(Tm, angle)
        CM.plane_stress_Tinv_matrix_eng!(Tinvm, angle)
        @test norm(Tm*Tinvm - I) < mtol
    end
    true
end
end
using .mlayup6
mlayup6.test()

module mlayup7
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    mtol = 1.0e-15
    R = [1 0 0; 0 1 0; 0 0 2]
    Tm = fill(0.0, 3, 3)
    Tinvm = fill(0.0, 3, 3)
    for dangle in [5 34 68 73 535 0 90 45 -135. -87]
        angle = dangle/180*pi
        CM.plane_stress_T_matrix_eng!(Tm, -angle)
        CM.plane_stress_Tinv_matrix_eng!(Tinvm, angle)
        @test norm(Tm - Tinvm) < mtol
    end
    true
end
end
using .mlayup7
mlayup7.test()

module mlayup8
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    mtol = 1.0e-15
    R = [1 0 0; 0 1 0; 0 0 2]
    Tm = fill(0.0, 3, 3)
    Tinvm = fill(0.0, 3, 3)
    for dangle in [-55] # see Barbero, Introduction to composite materials, Example 5.3
        angle = dangle/180*pi
        CM.plane_stress_T_matrix!(Tm, angle)
        @test norm(Tm - [0.32898992833716556 0.6710100716628343 -0.9396926207859083; 0.6710100716628343 0.32898992833716556 0.9396926207859083; 0.46984631039295416 -0.46984631039295416 -0.34202014332566877]   ) < mtol
    end
    true
end
end
using .mlayup8
mlayup8.test()

module mlayup9
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    mtol = 1.0e-15
    R = [1 0 0; 0 1 0; 0 0 2]
    Tme = fill(0.0, 3, 3)
    Tm = fill(0.0, 3, 3)
    Tinvm = fill(0.0, 3, 3)
    for dangle in [5 34 68 73 535 0 90 45 -135. -87]
        angle = dangle/180*pi
        CM.plane_stress_T_matrix_eng!(Tme, angle)
        CM.plane_stress_T_matrix!(Tm, angle)
        # Verify the identity above equation 5.54 in Barbero
        @test norm(inv(Tm') - Tme) < mtol
    end
    true
end
end
using .mlayup9
mlayup9.test()

module mlayup10
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    E = 200e3*phun("MPa")
    nu = 0.3;
    rho= 8000*phun("KG/M^3");
    thickness = 0.005*phun("m");

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
    CM = CompositeLayupModule
    A = fill(0.0, 3, 3)
    B = fill(0.0, 3, 3)
    C = fill(0.0, 3, 3)
    Atrue = E/(1-nu^2) * [1 nu 0; nu 1 0; 0 0 (1-nu)/2] * thickness

    for angle in [0.0 47.0 90.0 129.0 180.0]
        ply1 = CM.Ply("ply1", mater, thickness, angle)
        cl = CM.CompositeLayup("sample", [ply1, ], CM.cartesian_csys((1, 2, 3)))
        @test cl.plies[1].angle == angle
        A, B, C = CM.laminate_stiffnesses!(cl, A, B, C)
        @test norm(A - Atrue) < 1.0e-15 * norm(Atrue)
    end
    true
end
end
using .mlayup10
mlayup10.test()

module mlayup11
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    E = 200e3*phun("MPa")
    nu = 0.3;
    rho= 8000*phun("KG/M^3");
    thickness = 0.005*phun("m");
    angle = 45.0
    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
    CM = CompositeLayupModule
    # From Barbero's Introduction ... book Example 5.6
    ply1 = CM.Ply("ply1", mater, 1.0, -55.0)
    cl = CM.CompositeLayup("sample", [ply1, ], CM.cartesian_csys((1, 2, 3)))
    cl.plies[1]._Dps .= [20.874 3.260 0; 3.260 11.898 0; 0 0 3.789]
    A = fill(0.0, 3, 3)
    B = fill(0.0, 3, 3)
    C = fill(0.0, 3, 3)
    A, B, C = CM.laminate_stiffnesses!(cl, A, B, C)
    Atrue = [12.401509954148374 5.709503642606021 -1.217123826408347; 5.709503642606021 15.471482760639574 -3.000216655678807; -1.2171238264083473 -3.000216655678808 6.2385036426060205]
    @test norm(A - Atrue) < 1.0e-15 * norm(Atrue)
    true
end
end
using .mlayup11
mlayup11.test()

module mlayup12
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    # From Barbero's Introduction ... book Example 5.7
    E1 = 67192*phun("MPa")
    E2 = 12139*phun("MPa")
    G12 = 7638*phun("MPa")
    nu12 = 0.365;
    rho= 8000*phun("KG/M^3");
    thickness = 1.0*phun("m");
    mater = MatDeforElastOrtho(DeforModelRed3D, rho, E1, E2, E2, nu12, nu12, nu12, G12, G12, G12, 0.0, 0.0, 0.0)
    CM = CompositeLayupModule
    A = fill(0.0, 3, 3)
    B = fill(0.0, 3, 3)
    C = fill(0.0, 3, 3)

    for angle in [45]
        ply1 = CM.Ply("ply1", mater, thickness, angle)
        # @show ply1._Dps./phun("MPa")
        @test norm(ply1._Dps./phun("MPa") - [68849.10243290635 4540.006665496835 0.0;
                                                4540.006665496834 12438.374426018723 0.0;
                                                0.0 0.0 7637.999999999999]) < 1.0e-15 * E2
        cl = CM.CompositeLayup("sample", [ply1, ], CM.cartesian_csys((1, 2, 3)))
        @test cl.plies[1].angle == angle
        A, B, C = CM.laminate_stiffnesses!(cl, A, B, C)
        # @show A./phun("MPa")
        @test norm(A ./ phun("MPa") - [30229.872547479692 14953.872547479687 14102.68200172191;
                                       14953.872547479687 30229.872547479685 14102.682001721907;
                                       14102.682001721909 14102.682001721909 18051.865881982852]) < 1.0e-15 * E2
    end
    true
end
end
using .mlayup12
mlayup12.test()

module mlayup13
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    # From Barbero's Introduction ... book Example 5.7
    E1 = 67192*phun("MPa")
    E2 = 12139*phun("MPa")
    G12 = 7638*phun("MPa")
    nu12 = 0.365;
    rho= 8000*phun("KG/M^3");
    thickness = 1.0*phun("m");
    CM = CompositeLayupModule
    mater = CM.lamina_material(E1, E2, nu12, G12, G12, G12)
    A = fill(0.0, 3, 3)
    B = fill(0.0, 3, 3)
    C = fill(0.0, 3, 3)

    for angle in [45]
        ply1 = CM.Ply("ply1", mater, thickness, angle)
        # @show ply1._Dps./phun("MPa")
        @test norm(ply1._Dps./phun("MPa") - [68849.10243290635 4540.006665496835 0.0;
                                                4540.006665496834 12438.374426018723 0.0;
                                                0.0 0.0 7637.999999999999]) < 1.0e-15 * E2
        cl = CM.CompositeLayup("sample", [ply1, ], CM.cartesian_csys((1, 2, 3)))
        @test cl.plies[1].angle == angle
        A, B, C = CM.laminate_stiffnesses!(cl, A, B, C)
        # @show A./phun("MPa")
        @test norm(A ./ phun("MPa") - [30229.872547479692 14953.872547479687 14102.68200172191;
                                       14953.872547479687 30229.872547479685 14102.682001721907;
                                       14102.682001721909 14102.682001721909 18051.865881982852]) < 1.0e-15 * E2
    end
    true
end
end
using .mlayup13
mlayup13.test()

module mlayup14
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    # From Barbero's Introduction ... book Example 5.7
    E1 = 67192*phun("MPa")
    E2 = 12139*phun("MPa")
    G12 = 7638*phun("MPa")
    nu12 = 0.365;
    rho= 8000*phun("KG/M^3");
    thickness = 0.5*phun("m");
    CM = CompositeLayupModule
    mater = CM.lamina_material(E1, E2, nu12, G12, G12, G12)
    A = fill(0.0, 3, 3)
    B = fill(0.0, 3, 3)
    C = fill(0.0, 3, 3)
    # "balanced" fabric
    ply1 = CM.Ply("ply1", mater, thickness, 45)
    ply2 = CM.Ply("ply2", mater, thickness, -45)
    @test norm(ply1._Dps./phun("MPa") -
        [68849.10243290635 4540.006665496835 0.0;
        4540.006665496834 12438.374426018723 0.0;
        0.0 0.0 7637.999999999999]) < 1.0e-15 * E2
    @test norm(ply2._Dps./phun("MPa") -
        [68849.10243290635 4540.006665496835 0.0;
        4540.006665496834 12438.374426018723 0.0;
        0.0 0.0 7637.999999999999]) < 1.0e-15 * E2
    cl = CM.CompositeLayup("sample", [ply1, ply2], CM.cartesian_csys((1, 2, 3)))
    A, B, C = CM.laminate_stiffnesses!(cl, A, B, C)
        # @show A./phun("MPa")
    @test norm(A ./ phun("MPa") -
        [30229.872547479692 14953.872547479687 0.0;
        14953.872547479687 30229.872547479685 0.0;
        0.0 0.0 18051.865881982852] ) < 1.0e-15 * E2

    true
end
end
using .mlayup14
mlayup14.test()

module mlayup15
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    # From Barbero's Introduction ... book Example 5.7
    E1 = 67192*phun("MPa")
    E2 = 12139*phun("MPa")
    G12 = 7638*phun("MPa")
    nu12 = 0.365;
    rho= 8000*phun("KG/M^3");
    thickness = 0.001*phun("m");
    CM = CompositeLayupModule
    mater = CM.lamina_material(E1, E2, nu12, G12, G12, G12)
    A = fill(0.0, 3, 3)
    B = fill(0.0, 3, 3)
    C = fill(0.0, 3, 3)
    # "cross-ply" fabric: no coupling of extension and bending
    ply1 = CM.Ply("ply1", mater, thickness, 0)
    ply2 = CM.Ply("ply2", mater, thickness, 90)
    ply3 = CM.Ply("ply3", mater, thickness, 0)
    ply4 = CM.Ply("ply4", mater, thickness, 90)
    ply5 = CM.Ply("ply5", mater, thickness, 0)
    cl = CM.CompositeLayup("sample", [ply1, ply2, ply3, ply4, ply5], CM.cartesian_csys((1, 2, 3)))
    A, B, C = CM.laminate_stiffnesses!(cl, A, B, C)
        # @show A./phun("MPa")
    @test norm(B) < 1.0e-15 * E2

    true
end
end
using .mlayup15
mlayup15.test()

module mlayup16
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    # From Barbero's Introduction ... book Example 5.1
    E1 = 19981*phun("MPa")
    E2 = 11389*phun("MPa")
    G12 = 3789*phun("MPa")
    G13 = G12
    nu12 = 0.274;
    nu23 = 0.385
    G23 = E2 / 2 / (1 + nu23)
    rho= 8000*phun("KG/M^3");
    thickness = 0.001*phun("m");
    CM = CompositeLayupModule
    mater = CM.lamina_material(E1, E2, nu12, G12, G13, G23)
    A = fill(0.0, 3, 3)
    B = fill(0.0, 3, 3)
    C = fill(0.0, 3, 3)
    H = fill(0.0, 2, 2)
    # "cross-ply" fabric: no coupling of extension and bending
    ply1 = CM.Ply("ply1", mater, thickness, -55)
    cl = CM.CompositeLayup("sample", [ply1, ], CM.cartesian_csys((1, 2, 3)))
    A, B, C = CM.laminate_stiffnesses!(cl, A, B, C)
        # @show A./phun("GPa")/thickness
    @test norm(A./phun("GPa")/thickness - [
        12.401649598137476 5.70964888928013 -1.2171315230336193;
        5.70964888928013 15.47166067503043 -3.000261532585858;
        -1.2171315230336193 -3.0002615325858573 6.238554718598339]) < 1.0e-15 * E2
    H = CM.laminate_transverse_stiffness!(cl, H)
    # @show H./phun("GPa")/thickness
    @test norm(H./phun("GPa")/thickness - [
        3.3378632276560145 -0.1262916916205752;
        -0.12629169162057502 3.245930394485983]) < 1.0e-15 * E2
    true
end
end
using .mlayup16
mlayup16.test()

module mlayup17
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test()
    CM = CompositeLayupModule
    # From Barbero's Introduction ... book Example 5.1
    E1 = 19981*phun("MPa")
    E2 = 11389*phun("MPa")
    G12 = 3789*phun("MPa")
    G13 = G12
    nu12 = 0.274;
    nu23 = 0.385
    G23 = E2 / 2 / (1 + nu23)
    rho= 8000*phun("KG/
        M^3");
    lamina_thickness = 0.000635*phun("m");
    CM = CompositeLayupModule
    mater = CM.lamina_material(E1, E2, nu12, G12, G13, G23)
    A = fill(0.0, 3, 3)
    B = fill(0.0, 3, 3)
    C = fill(0.0, 3, 3)
    H = fill(0.0, 2, 2)
    # "angle-ply" fabric: coupling of extension and bending
    cl = CM.CompositeLayup("sample", [
        CM.Ply("ply1", mater, lamina_thickness, +55),
        CM.Ply("ply2", mater, lamina_thickness, -55)], CM.cartesian_csys((1, 2, 3)))
    A, B, C = CM.laminate_stiffnesses!(cl, A, B, C)
        # @show A./phun("GPa*mm")
        # @show B./phun("GPa*mm^2")
        # @show C./phun("GPa*mm^3")
    @test norm(A ./ phun("GPa*mm") - [15.750094989634595 7.251254089385766 0.0; 7.251254089385766 19.649009057288648 0.0; 0.0 0.0 7.922964492619891]) < 1.0e-6
    @test norm(B ./ phun("GPa*mm^2") - [0.0 0.0 -0.49077785837523136; 0.0 0.0 -1.2097804564769328; -0.4907778583752312 -1.2097804564769326 0.0]) < 1.0e-6
    @test norm(C ./ phun("GPa*mm^3") - [2.1169440173984704 0.9746289767308587 0.0; 0.9746289767308587 2.6409905590417386 0.0; 0.0 0.0 1.0649124525122187]) < 1.0e-6
    true
end
end
using .mlayup17
mlayup17.test()

module mlayup18
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using Test
function test(npairs = 1)
    # From Barbero's Introduction ... book Example 3.1
    # ASFD/9310
    E1 = 133860*phun("MPa")
    E2 = 7706*phun("MPa")
    G12 = 4306*phun("MPa")
    G13 = G12
    nu12 = 0.301;
    nu23 = 0.396
    G23 = 2760*phun("MPa")
    npairs = 1
    thickness = 10*phun("mm");
    des = Dict(1 => 1585193, 5 => 317039, 20 => 79260)

    CM = CompositeLayupModule

    mater = CM.lamina_material(E1, E2, nu12, G12, G13, G23)

    A = fill(0.0, 3, 3)
    B = fill(0.0, 3, 3)
    C = fill(0.0, 3, 3)

    plies = []
    for j in 1:npairs
        push!(plies, CM.Ply("ply1", mater, thickness/npairs/2, 0))
        push!(plies, CM.Ply("ply1", mater, thickness/npairs/2, 90))
    end
    cl = CM.CompositeLayup("sample", plies, CM.cartesian_csys((1, 2, 3)))

    A, B, C = CM.laminate_stiffnesses!(cl, A, B, C)
    # @show B
    v = des[npairs]
    Btrue = [-v 0 0; 0 v 0; 0 0 0]

    @test norm(B - Btrue) < 1.0e-6 * v

    true
end
end
using .mlayup18
mlayup18.test(1)
mlayup18.test(5)
mlayup18.test(20)


text:
using DashBootstrapComponents

button = dbc_button("Block button", color="primary", block=true);

ext: jl

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
import MetNets: av, va

# MetState interface
av(out::FBAOut) = out.v
va(out::FBAOut) = zeros(length(out.v))


text:
####################
# DerivativeResult #
####################

type DerivativeResult{V,D} <: ForwardDiffResult
    value::V
    derivative::D
end

DerivativeResult(x) = DerivativeResult(copy(x), copy(x))

value(result::DerivativeResult) = result.value
derivative(result::DerivativeResult) = result.derivative

###############
# API methods #
###############

derivative(f, x) = extract_derivative(f(Dual(x, one(x))))

function derivative!(out, f, x)
    result = f(Dual(x, one(x)))
    load_derivative_value!(out, result)
    load_derivative!(out, result)
    return out
end

#####################
# result extraction #
#####################

@inline extract_derivative(n::Real) = partials(n, 1)
@inline extract_derivative(arr) = load_derivative!(similar(arr, numtype(eltype(arr))), arr)

function load_derivative!(out, arr)
    for i in eachindex(out)
        out[i] = extract_derivative(arr[i])
    end
    return out
end

@inline function load_derivative!(out::DerivativeResult, arr)
    load_derivative!(out.derivative, arr)
    return out
end

@inline function load_derivative!(out::DerivativeResult, n::Real)
    out.derivative = extract_derivative(n)
    return out
end

@inline load_derivative_value!(out, result) = out

function load_derivative_value!(out::DerivativeResult, arr)
    val = out.value
    for i in eachindex(val)
        val[i] = value(arr[i])
    end
    return out
end

@inline function load_derivative_value!(out::DerivativeResult, n::Real)
    out.value = value(n)
    return out
end


text:
function softabs(hessian::RealMatrix, a::Float64=1000.0)
  λ, Q = eig(hessian)
  Q*diagm(λ./tanh.(a*λ))*Q'
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -6.103907844236108 <= q <= 8.378675476959124)

                     add_NL_constraint(m, :(sigmoid(-0.8297607126794757 + -0.5064858067335005*sigmoid(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2])) + -0.8619658571743893*sigmoid(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2])) + 0.8415580445654567*sigmoid(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))) + sigmoid(-0.0809505831099564 + 0.6185835236357913*sigmoid(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2])) + -0.5288114683363174*sigmoid(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2])) + -0.42070009030061817*sigmoid(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))) + sigmoid(0.16011867292674786 + -0.14898607305663658*sigmoid(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2])) + -0.8220978033005912*sigmoid(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2])) + -0.8518875556726564*sigmoid(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# agent-based economy under covid infection
using Agents,Random,Plots, Statistics

snapshot = false # set true if snapshot should be made at end of run
loadsnapshot = true # set true if initialization comes from file "snapshot.dat" else random initialization
snapname = "snapshot100kr1.dat"

#snapshot = true # set true if snapshot should be made at end of run
#loadsnapshot = false # set true if initialization comes from file "snapshot.dat" else random initialization

#Random.seed!(20042020)
#filename_prefix = ""

#parameters
k1 = 1 # dimension1 of space
k2 = 1 # dimension2 of space
nhh = 100000 # number of hh

T = 744  # time periods
virustime =14 # time where infection starts
vacctime = virustime + 547
enablebankruptciestime = 1000 #1000
endfcadj = 1250
resetfirmsavingstime = 500

datat = 7 # interval of data collection
bankruptcypossible = false #true
if loadsnapshot   #adjust fixed costs in both directions only if no snapshot is loaded
    updatefixedcosts = false
else
    updatefixedcosts = true
end

#policy data
alpha_lockdown = 1.0
alpha_open = 0.0
tpol = virustime + 21
policies = Dict{Int64,String}()
policies[tpol] = "policy_lockdown_2.jl"
policies[vacctime] = "policy_allout.jl"
detfrac = 0.15 # detection frequency
adaptivepolicythresholdoff = 5 * nhh / (detfrac*100000) # Threshold infected/week
adaptivepolicythresholdon = 100 * nhh / (detfrac*100000) # Threshold infected/week
poladjfrac = 0.6 # factor by which gap in policy parameter is closeed during phase in / phase out every week

tadaptivepolicystart = tpol + 21 # Enable adaptive policy response
tadaptivepolicyend = vacctime - 7 # Enable adaptive policy response
adaptivepolicygood = "policy_open_2.jl" # policy acitvated if number of infected above threshold
adaptivepolicybad = "policy_lockdown_2.jl" #  policy activated if number of infected below threshold

fracy = .75 # fraction of young households
avempl = 15 # average number of employees per firm
nf =Int(floor(fracy*nhh/avempl)) # number of firms
nsec = 3 # number of sectors: man / ser / food
pubfrac = 0.12 # employment share in public sector
fraccons = [0.21, 0.5, 0.29] # fraction consumption man / ser / food
essentialsec = 3
essentialadj = 0.01
baseprod = [97, 62, 48, 62] # average productivities in 1000€ man / ser / food /pub
#fracemp = [0.14, 0.49, 0.37] # initial fraction employment man / ser / food
fracemp = (fraccons ./ baseprod[1:3]) / sum(fraccons ./ baseprod[1:3])

#phome = 0.9 # prob of being in own region
ptravel = 0.0 # prob of being in other region, only for young
socialmaxyy = 5 # max no. people met outside firm young_young
socialmaxoo = 2 # max no. people met outside firm old_old
socialmaxoy = 4 # max no. people met outside firm old_old
socialmaxyo = 2 # max no. people met outside firm young_old
workmeet = [8,8,8,8] # max number of people met at work: man / ser / food /public
shopmeet = [10,28,10] # max number of people met at each shopping instant  man / ser / food
homeofficefrac = [.45,.3,.0 ,0.75] # data taken from Fadinger/Schymik: man / ser / food /public
phomeoffice = 0

pinf = 1.25*[1,1.5] # infection prob at meeting before/after mutation
muttime =9999 # time of mutation
mutnum = 5 # number of mutating
trec = 21 # recovery time
corlatent = 5 # latency time
corinf = 5 # infectious period
hcap = Int(0.5*round(0.0003*(nhh))) # no of intensive beds based on 30 per 100.000
mortl = detfrac*[0.0066/trec, 0.16/trec] # mortality young, old, normal cap, data from Germany
morth = detfrac*[0.018/trec, 0.5/trec] # mortality young, old, overcap, data from Italy
icufrac = detfrac*0.085 # fraction of infected needing intensive care unit

#firm markup
etmin = [0.18, 0.18, 0.18, 0] #profit markups after fixed costs are introduced
etmax = [0.25, 0.25, 0.25, 0] #profit markups after fixed costs are introduced

#et = 0.15
#inventory buffer for firms  # man / ser / food

buff = .25*[0.1 / 7., 0.03 / 7, 0.05 / 7, 0] # changes
#buff = .005*[0.1 / 7., 0.05 / 7, 0.05 / 7, 0]

#buff = [0.2 / 7., 0.05 / 7, 0.1 / 7, 0]

# weekly depreciation of inventory
de = [0.01, 1, 0.5, 0] # man / service / food
rhodem = 0.5 # demand expectation smoothing parameter
zerostockboost = 8 # factor to boost expansion if stock is zero
divrate = 0.7 # dividend rate if savings are not too high
savtar = [1,.5,.5,0] # target of firm savings relative to av. revenues during last 4 weeks

#fcvcratio = 0.1*[0.47, 0.3, 0.3, 0.3] #fixed cost to variable cost ratio
fcvcratio = 0.16*[0.47, 0.3, 0.3, 0.3]

#productivity distribution
#prodlb = [.9, 1. / 1.1, .9, 1.]    # man / ser / food
#produb = [1.1, 1.1, 1.1, 1.]    # man / ser / food
prodlb = [.9, 1.05 , .9, 1.]    # man / ser / food
produb = [1.1, 0.95, 1.1, 1.]    # man / ser / food

wage = baseprod .* (1 .- buff .* de) ./ ((1 .+ fcvcratio) .* (1 .+ etmin)) # weekly wage
pension = 0.5 * mean(wage)
unempbenefit = 0.6 .* wage[1:nsec]
shorttimewage = 0.7 .* wage[1:nsec]

#consumption budget rule parameters
kap = 0.1 / 4. #marginal  prop to consume from exzess saving
phi = 4*8.  # target savings / weekly income ratio (16.67 per month in Eurace)
rhoinc = 0.4 # income expectation smoothing parameter
pshop = [[1,1,1],[1,1,1]] # prob of weekly shopping in man/ser/food
numshops = [4,4,4] # numer of local suppliers consumer checks each shopping: man/ser/food
gac = 16 # intensity of choice parameters for consumers

# government
tauadj = 0.05 # speed of adjusment of tax rate
debtfrac = 1. / (52. * 10.) # fraction of public debt erased'/added in one week

#calculate agent numbers
nhhy = Int(floor(fracy*nhh)) # number of young hh
nhho = nhh - nhhy # number of old hh
nfsec = [0 for tt=1:nsec+1]
for tt=1:nsec
    nfsec[tt] = Int(floor((1-pubfrac)*fracemp[tt]*nf)) # number  of firms in sector tt
end
nfsec[nsec+1] = nf - sum(nfsec[1:nsec]) # number of public sector offices

#initialization
consumption = zeros(nsec, k1, k2)
inifrac = 0.0001 # fraction of pop infected when virus starts, Germany 9.3. 1/10000 taken into account undetected
ininumy = Int(ceil(fracy*inifrac*nhh))
ininumo = Int(ceil((1-fracy)*inifrac*nhh))
iniinf = [[[ininumy,ininumo]]]
trigger = true # initial infection can start
unemp = 0.1 # initial unemployment
incomehh = (1-unemp)*wage
savinghh = phi*incomehh
incomehho = pension
savinghho = phi*incomehho
savingf  = 100
consbud = [0.,0.,0.] #man / service / food
prevconsbud = [0.,0.,0.] #man / service / food
shoprepini = [true, true, true] #  possible repetition of shopping trip
daycounter = 6
tau = (fracy*unemp*mean(unempbenefit) + (1-fracy)*pension) / ((1-unemp)*fracy*mean(wage .* (1 .+ etmin))) # initial tax rate
pubacc = 0. # public account
pastrevini = 10000*[1,1,1,1] # intitializations of past revenue
virus = false # initially there is no virus
genhomeoffice = false # initially no home-office
bailoutprogram = false # initially no bailouts
end_bailout_after_weeks = 78
shorttimeprogram = false # initially no short term program
shorttimeperiod = 78
shorttimeprob = 0.9
demandexpini = [1700,1000,600,0]

# list for total consumption i regions
regcons = [] # saves regional consumption every day

cas = [[0,0] for i=1:k1, j=1:k2] # casualties every day

# save default values
pinf_o = pinf
socialmaxyy_o = socialmaxyy
socialmaxoy_o = socialmaxoy
socialmaxoo_o = socialmaxoo
socialmaxyyh = socialmaxyy
socialmaxoyh = socialmaxoy
socialmaxooh = socialmaxoo
ptravel_o = ptravel
pshop_o = pshop
tauadj_o = tauadj
phomeoffice_o = phomeoffice
workmeet_o = workmeet
shopmeet_o = shopmeet
workmeeth = workmeet
shopmeeth = shopmeet
#set initial targets equal to current values
pinft = pinf
socialmaxyyt = socialmaxyy
socialmaxoyt = socialmaxoy
socialmaxoot = socialmaxoo
pshopt = pshop
phomeofficet = phomeoffice
workmeett = workmeet
shopmeett = shopmeet
socialmaxyo_o = socialmaxyo
socialmaxyoh = socialmaxyo
socialmaxyot = socialmaxyo


text:
include("../tools/plotting.jl")
include("../settings.jl")

"""
Structure for storing summary statistics
"""
struct sts
    median::Float64
    intwidth::Float64
end

"""
Compute MSE for LOOCV
"""
function compute_loss(train::trainingData{A, B}, rangeθ::B, rangeλ::B, transform, quadtype = "Gaussian", priortype = "Uniform", graph = "true", num_pts=100, strt = 0, endpt=2.5, tol = 1e-5) where A<:Array{Float64, 2} where B<:Array{Float64, 1}
    X = train.X; z = train.z; s = train.s
    medians = similar(z)
    x = collect(range(strt, stop = endpt, length = num_pts)) #define mesh
    Xs = repeat(x, outer = [1, length(z)]) #preallocate space
    Ys = Array{Float64}(undef, num_pts, length(z))

    for i = 1:length(z)
        println("i: ", i)
        ind = [collect(1:i-1); collect(i+1:length(z))]
        train_cur = trainingData(s[ind, :], X[ind, :], z[ind])
        test_cur = testingData(s[i:i, :], X[i:i, :])
        (pdf, cdf, _) = getBtgDensity(train_cur, test_cur, rangeθ, rangeλ, transform, quadtype, priortype)
        med_cur = median(pdf, cdf) #first median-finding routine
        medians[i] = med_cur
        #println("cdf(medcur)-0.5", abs(cdf(med_cur)-0.5))
        #println("median: ", med_cur)
        #println("actual z: ", z[i])
        #if (abs(cdf(med_cur)-0.5)>tol)
        if graph == "true"#one expensive loop
            for j = 1:num_pts
                Ys[j, i] = pdf(x[j])
            end
        end
    end
    return (medians, Xs, Ys)

    #norm(medians .- z)^2/length(z)
end

"""
Computes (x, y)-values for plotting distributions arising from LOOCV (Leave one out cross validation).
If n is size of training set, then return values have following dimensions
    Xs: num_pts by n matrix
    Ys: num_pts by n matrix
"""

function cross_validate(train::trainingData{A, B}, rangeθ::B, rangeλ::B, transform, quadtype = "Gaussian", priortype = "Uniform", num_pts=100, strt = 0, endpt=1.5) where A<:Array{Float64, 2} where B<:Array{Float64, 1}
    X = train.X; z = train.z; s = train.s
    x = collect(range(strt, stop = endpt, length = num_pts)) #define mesh
    Xs = repeat(x, outer = [1, length(z)]) #preallocate space
    Ys = Array{Float64}(undef, num_pts, length(z))
    for i=1:length(z)
        #println(i)
        ind = [collect(1:i-1); collect(i+1:length(z))]
        #@time begin
        train_cur = trainingData(s[ind, :], X[ind, :], z[ind])
        test_cur = testingData(s[i:i, :], X[i:i, :])
        pdf, cdf = getBtgDensity(train_cur, test_cur, rangeθ, rangeλ, transform, quadtype, priortype)
        for j = 1:num_pts
            Ys[j, i] = pdf(x[j])
        end
    end
        return Xs, Ys
end

"""
Single point cross-validation. Currently does not use confidence intervals or median finding.
"""
function cross_validate_artifact(X, s, g, gprime, pθ, pλ, z, rangeθ, rangeλ, num_pts=200, strt = 0, endpt=20)
    ss = Array{sts}(undef, length(z))
    x = collect(range(strt, stop = endpt, length = num_pts)) #define mesh
    Xs = repeat(x, outer = [1, length(z)]) #preallocate space
    Ys = Array{Float64}(undef, num_pts, length(z))
    for i=1:length(z)
        println(i)
        ind = [collect(1:i-1); collect(i+1:length(z))]
        @time begin
        X0 = size(X, 2)==1 ? X[1:1, :] : X[i:i, :] # constant mean?
        pdf, cdf = model(X[ind, :], X0, s[ind, :], s[i:i, :], boxCox, boxCoxPrime, pθ, pλ, z[ind], range_theta, range_lambda) #X0 is [1]
        constant = cdf(30)
        pdfn = x -> pdf(x)/constant
        end
        #@time begin
        #med = bisection(x -> cdfn(x)-0.5, 1e-3, 50, 1e-3, 15) #compute median
        #end
        #@time begin
        #intwidth = confidence(pdfn, med) #compute width of credible interval corresponding to 95% confidence level
        #end
        #st = sts(med, intwidth); ss[i]=st
        for j = 1:num_pts
            Ys[j, i] = pdfn(x[j])
        end
    end
    return ss, Xs, Ys
end

"""
Delete-one-group cross validation.
"""
function cross_validate_groups(X, s, g, gprime, pθ, pλ, z, rangeθ, rangeλ, num_group = 5, num_pts=200, strt = 0, endpt=20)
    if num_group>length(z)
        error("number of splits must be less than or equal to number of data points")
    end
    x = collect(range(strt, stop = endpt, length = num_pts)) #define mesh
    Xs = repeat(x, outer = [1, length(z)]) #preallocate space to store
    Ys = Array{Float64}(undef, num_pts, length(z))
    num_per_group = Integer(ceil(length(z)/num_group))
    for i=1:num_group
        ind_train = [collect(1:(i-1)*num_per_group); collect(i*num_per_group+1:length(z))] #delete the ith group
        println(ind_train)
        for k=1:num_per_group
            cur = (i-1)*num_per_group + k
            println(cur)
            if cur <= length(z)
                X0 = size(X, 2)==1 ? X[1:1, :] : X[cur:cur, :] # set design matrix
                @time begin
                pdf, cdf = model(X[ind_train, :], X0, s[ind_train, :], s[cur:cur, :], boxCox, boxCoxPrime, pθ, pλ, z[ind_train], range_theta, range_lambda)
                constant = cdf(30) #compute unnormalized cdf
                cdfn = x -> cdf(x)/constant; pdfn = x -> pdf(x)/constant  #compute normalized pdf and cdf
                end
                for j = 1:num_pts
                    Ys[j, cur] = pdfn(x[j])
                end
            end
        end
    end
    return Xs, Ys
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
`Transition{T<:AbstractFloat}`

The transition equation of the state-space model takes the form

    `s_t = TTT*s_{t-1} + RRR*ϵ_t + CCC`

The `Transition` type stores the coefficient `Matrix{T}`s (`TTT`, `RRR`) and constant `Vector{T} CCC`.
"""
mutable struct Transition{T<:AbstractFloat}
    TTT::Matrix{T}
    RRR::Matrix{T}
    CCC::Vector{T}
end
function Transition(TTT::Matrix{T}, RRR::Matrix{T}) where T<:AbstractFloat
    CCC = zeros(eltype(TTT), size(TTT, 1))
    Transition{T}(TTT, RRR, CCC)
end
function Transition(TTT::Matrix{T}, RRR::Matrix{T}, CCC::Matrix{T}) where T<:AbstractFloat
    Transition{T}(TTT, RRR, collect(CCC))
end
function Base.getindex(eq::Transition, d::Symbol)
    if d in (:TTT, :RRR, :CCC)
        return getfield(eq, d)
    else
        throw(KeyError(d))
    end
end

"""
`Measurement{T<:AbstractFloat}`

The measurement equation of the state-space model takes the form

   `y_t = ZZ*s_t + DD + u_t`

where the error `u_t` is the measurement error, which is uncorrelated with the
shocks in the transition equation `ϵ_t`.

### Fields

If `Ns` is the number of states `s_t`, `Ny` is the number of
observables `y_t`, and `Ne` is the number of shocks `ϵ_t`:

- `ZZ`: the `Ny` x `Ns`  measurement matrix
- `DD`: the `Ny` x 1 constant vector
- `QQ`: the `Ne` x `Ne` covariance matrix for the shocks `ϵ_t`
- `EE`: the `Ny` x `Ny` covariance matrix for the measurement error `η_t`
"""
mutable struct Measurement{T<:AbstractFloat}
    ZZ::Matrix{T}
    DD::Vector{T}
    QQ::Matrix{T}
    EE::Matrix{T}
end

function Base.getindex(M::Measurement, d::Symbol)
    if d in (:ZZ, :DD, :QQ, :EE)
        return getfield(M, d)
    else
        throw(KeyError(d))
    end
end

function measurement(m::AbstractDSGEModel, trans::Transition; shocks::Bool=true)
    TTT = trans[:TTT]
    RRR = trans[:RRR]
    CCC = trans[:CCC]
    measurement(m, TTT, RRR, CCC; shocks=shocks)
end

"""
```
PseudoMeasurement{T<:AbstractFloat}
```

The pseudo-measurement equation of the state-space model takes the form

   `x_t = ZZ_pseudo*s_t + DD_pseudo`

### Fields

Let `Ns` be the number of states `s_t` and `Nx` be the number of
pseudo-observables `x_t`:

- `ZZ_pseudo`: the `Nx` x `Ns` pseudo-measurement matrix
- `DD_pseudo`: the `Nx` x 1 constant vector
"""
mutable struct PseudoMeasurement{T<:AbstractFloat}
    ZZ_pseudo::Matrix{T}
    DD_pseudo::Vector{T}
end

function Base.getindex(M::PseudoMeasurement, d::Symbol)
    if d in (:ZZ_pseudo, :DD_pseudo)
        return getfield(M, d)
    else
        throw(KeyError(d))
    end
end

"""
`System{T<:AbstractFloat}`

A mutable struct containing the transition and measurement equations for a
state-space model. The matrices may be directly indexed: `sys[:TTT]`
returns `sys.transition.TTT`, `sys[:ZZ]` returns `sys.measurement.ZZ`, etc.
"""
mutable struct System{T<:AbstractFloat}
    transition::Transition{T}
    measurement::Measurement{T}
    pseudo_measurement::PseudoMeasurement{T}
end

function System(transition::Transition{T}, measurement::Measurement{T}) where T<:AbstractFloat
    # Initialize empty pseudo-measurement equation
    _n_states = size(transition.TTT, 1)
    _n_pseudo = 0
    ZZ_pseudo = zeros(_n_pseudo, _n_states)
    DD_pseudo = zeros(_n_pseudo)
    pseudo_measurement = PseudoMeasurement(ZZ_pseudo, DD_pseudo)

    return System(transition, measurement, pseudo_measurement)
end

function Base.getindex(system::System, d::Symbol)
    if d in (:transition, :measurement, :pseudo_measurement)
        return getfield(system, d)
    elseif d in fieldnames(typeof(system.transition))
        return getfield(system.transition, d)
    elseif d in fieldnames(typeof(system.measurement))
        return getfield(system.measurement, d)
    elseif d in fieldnames(typeof(system.pseudo_measurement))
        return getfield(system.pseudo_measurement, d)
    else
        throw(KeyError(d))
    end
end

function Base.copy(system::System)
    trans = Transition(system[:TTT], system[:RRR], system[:CCC])
    meas  = Measurement(system[:ZZ], system[:DD], system[:QQ], system[:EE])
    pseudo_meas = PseudoMeasurement(system[:ZZ_pseudo], system[:DD_pseudo])
    return System(trans, meas, pseudo_meas)
end

"""
```
compute_system(m; apply_altpolicy = false)
```

Given the current model parameters, compute the state-space system
corresponding to model `m`. Returns a `System` object.
"""
function compute_system(m::AbstractDSGEModel{T}; apply_altpolicy = false,
                        verbose::Symbol = :high) where T<:AbstractFloat

    solution_method = get_setting(m, :solution_method)

    # Solve model
    if solution_method == :gensys

        TTT, RRR, CCC = solve(m; apply_altpolicy = apply_altpolicy, verbose = verbose)
        transition_equation = Transition(TTT, RRR, CCC)

        # Solve measurement equation
        measurement_equation = measurement(m, TTT, RRR, CCC)

    elseif solution_method == :klein
        # Unpacking the method from solve to hang on to TTT_jump
        if m.spec == "het_dsge"
            TTT_jump, TTT_state, eu = klein(m)
        else
            TTT_jump, TTT_state, eu = klein(m)
        end
        if eu==-1
            throw(KleinError())
        end

        TTT, RRR = klein_transition_matrices(m, TTT_state, TTT_jump)
        CCC = zeros(n_model_states(m))

        if m.spec == "real_bond_mkup"
            GDPeqn = construct_GDPeqn(m, TTT_jump)
            TTT, RRR, CCC = augment_states(m, TTT, TTT_jump, RRR, CCC, GDPeqn)
            # Measurement (needs the additional TTT_jump argument)
            measurement_equation = measurement(m, TTT, TTT_jump, RRR, CCC, GDPeqn)
        elseif m.spec == "het_dsge" || m.spec == "rep_dsge"
            TTT, RRR, CCC = augment_states(m, TTT, RRR, CCC)
            measurement_equation = measurement(m, TTT, RRR, CCC)
        else
            TTT, RRR, CCC        = augment_states(m, TTT, RRR, CCC)
            measurement_equation = measurement(m, TTT, RRR, CCC)
        end

        transition_equation = Transition(TTT, RRR, CCC)

    else
        throw("solution_method provided does not exist.")
    end

    type_tuple = (typeof(m), Matrix{T}, Matrix{T}, Vector{T})
    if hasmethod(pseudo_measurement, type_tuple)
        # Solve pseudo-measurement equation
        pseudo_measurement_equation = pseudo_measurement(m, TTT, RRR, CCC)
        return System(transition_equation, measurement_equation, pseudo_measurement_equation)
    else
        return System(transition_equation, measurement_equation)
    end
end

"""
```
compute_system(m; apply_altpolicy = false,
               regime_switching = false, n_regimes = 2,
               check_system = false, get_system = false,
               get_population_moments = false, use_intercept = false,
               verbose = :high)
compute_system(m, data; apply_altpolicy = false,
               regime_switching = false, n_regimes = 2,
               check_system = false, get_system = false,
               get_population_moments = false,
               verbose = :high)
```
Given the current model parameters, compute the DSGE-VAR or DSGE-VECM system
corresponding to model `m`. If a matrix `data` is also passed, then
the VAR is estimated on `data` using the DSGE `m` as a prior
with weight λ.

### Keyword Arguments
* `check_system::Bool`: see `?compute_system` that takes the input `m::AbstractDSGEModel`
    and `system::System`.
* `get_system::Bool`: see Outputs
* `get_population_moments::Bool`: see Outputs
* `use_intercept::Bool`: use an intercept term when computing the OLS estimate of the VAR system.

### Outputs
* If `get_system = true`:
    Returns the updated `system` whose measurement matrices `ZZ`, `DD`, and `QQ` correspond
    to the VAR or VECM specified by `m`. If `m` is an `AbstractDSGEVECMModel`,
    then the `system` and the vector implied by additional cointegrating relationships
    are returned as a 2-element tuple.
* If `get_population_moments = true`:
    Returns the limit cross product matrices that describe the DSGE implied
    population moments between the observables and their lags. If `data` is
    also passed as an input, then the sample population moments are also returned.
* Otherwise:
    Returns `β` and `Σ`, the coefficients and observables covariance matrix of the VAR or VECM.
    If `data` is passed in, then `β` and `Σ` are estimated from the data using `m`
    as a prior with weight λ. Otherwise, `β` and `Σ` comprise the VECM approximation
    of the DSGE `m`.
"""
function compute_system(m::AbstractDSGEVARModel{T}; apply_altpolicy::Bool = false,
                        regime_switching::Bool = false, regime::Int = 1, n_regimes::Int = 2,
                        check_system::Bool = false, get_system::Bool = false,
                        get_population_moments::Bool = false, use_intercept::Bool = false,
                        verbose::Symbol = :high) where {T<:Real}
    dsge = get_dsge(m)
    if regime_switching
        regime_system = compute_system(dsge; apply_altpolicy = apply_altpolicy,
                                       regime_switching = regime_switching, n_regimes = n_regimes,
                                       verbose = verbose)
        system = System(regime_system, regime)
    else
        system = compute_system(dsge; verbose = verbose)

    end
    # Update the system to reflect the desired observables using information
    # from the underlying DSGE object
    system = compute_system(dsge, system; observables = collect(keys(get_observables(m))),
                            shocks = collect(keys(get_shocks(m))), check_system = check_system)

    if get_system
        return system
    else
        EE, MM = measurement_error(m)

        return var_approx_state_space(system[:TTT], system[:RRR], system[:QQ],
                                      system[:DD], system[:ZZ], EE, MM, n_lags(m);
                                      get_population_moments = get_population_moments,
                                      use_intercept = use_intercept)
    end
end

function compute_system(m::AbstractDSGEVARModel{T}, data::Matrix{T};
                        apply_altpolicy::Bool = false,
                        regime_switching::Bool = false, regime::Int = 1, n_regimes::Int = 2,
                        check_system::Bool = false, get_system::Bool = false,
                        get_population_moments::Bool = false,
                        verbose::Symbol = :high) where {T<:Real}

    if get_λ(m) == Inf
        # Then we just want the VAR approximation of the DSGE
        return compute_system(m; apply_altpolicy = apply_altpolicy,
                              regime_switching = regime_switching, regime = regime, n_regimes = n_regimes,
                              check_system = check_system, get_system = get_system,
                              get_population_moments = get_population_moments, use_intercept = true,
                              verbose = verbose)
    else
        # Create a system using the method for DSGE-VARs and λ = ∞
        system = compute_system(m; apply_altpolicy = apply_altpolicy,
                                regime_switching = regime_switching, regime = regime,
                                n_regimes = n_regimes, check_system = check_system,
                                get_system = true, use_intercept = true,
                                verbose = verbose)

        if get_system
            return system
        else
            EE, MM = measurement_error(m)

            lags = n_lags(m)
            YYYY, XXYY, XXXX =
                compute_var_population_moments(data, lags; use_intercept = true)
            out = var_approx_state_space(system[:TTT], system[:RRR], system[:QQ],
                                         system[:DD], system[:ZZ], EE, MM, n_lags(m);
                                         get_population_moments = true,
                                         use_intercept = true)

            if get_population_moments
                return out..., YYYY, XXYY, XXXX
            else
                # Compute prior-weighted population moments
                λ = get_λ(m)
                YYYYC = YYYY + λ .* out[1]
                XXYYC = XXYY + λ .* out[2]
                XXXXC = XXXX + λ .* out[3]

                # Draw stationary VAR system
                n_periods = size(data, 2) - lags
                β, Σ =  draw_stationary_VAR(YYYYC, XXYYC, XXXXC,
                                            convert(Int, n_periods + λ * n_periods),
                                            size(data, 1), lags)

                return β, Σ
            end
        end
    end
end

# Same functions as above but for AbstractDSGEVECMModel types
function compute_system(m::AbstractDSGEVECMModel{T}; apply_altpolicy::Bool = false,
                        regime_switching::Bool = false, regime::Int = 1, n_regimes::Int = 2,
                        check_system::Bool = false, get_system::Bool = false,
                        get_population_moments::Bool = false, use_intercept::Bool = false,
                        verbose::Symbol = :high) where {T<:Real}
    dsge = get_dsge(m)
    if regime_switching
        regime_system = compute_system(dsge; apply_altpolicy = apply_altpolicy,
                                       regime_switching = regime_switching, n_regimes = n_regimes,
                                       verbose = verbose)
        system = System(regime_system, regime)
    else
        system = compute_system(dsge; verbose = verbose)

    end
    # Use wrapper compute_system for AbstractDSGEVECMModel types
    # (as opposed to compute_system(m::AbstractDSGEModel, system::System; ...))
    system, DD_coint_add = compute_system(m, system; observables = collect(keys(get_observables(m))),
                                          cointegrating = collect(keys(get_cointegrating(m))),
                                          cointegrating_add = collect(keys(get_cointegrating_add(m))),
                                          shocks = collect(keys(get_shocks(m))), check_system = check_system,
                                          get_DD_coint_add = true)

    if get_system
        return system, DD_coint_add
    else
        EE, MM = measurement_error(m)

        return vecm_approx_state_space(system[:TTT], system[:RRR], system[:QQ],
                                       system[:DD], system[:ZZ], EE, MM, n_observables(m),
                                       n_lags(m), n_cointegrating(m), n_cointegrating_add(m),
                                       DD_coint_add;
                                       get_population_moments = get_population_moments,
                                       use_intercept = use_intercept)
    end
end

function compute_system(m::AbstractDSGEVECMModel{T}, data::Matrix{T};
                        apply_altpolicy::Bool = false,
                        regime_switching::Bool = false, regime::Int = 1, n_regimes::Int = 2,
                        check_system::Bool = false, get_system::Bool = false,
                        get_population_moments::Bool = false,
                        verbose::Symbol = :high) where {T<:Real}

    if get_λ(m) == Inf
        # Then we just want the VECM approximation of the DSGE
        # with no additional cointegration
        return compute_system(m; apply_altpolicy = apply_altpolicy,
                              regime_switching = regime_switching, regime = regime, n_regimes = n_regimes,
                              check_system = check_system, get_system = get_system,
                              get_population_moments = get_population_moments, use_intercept = true,
                              verbose = verbose)
    else
        # Create a system using the method for DSGE-VECMs and λ = ∞
        system, DD_coint_add = compute_system(m; apply_altpolicy = apply_altpolicy,
                                regime_switching = regime_switching, regime = regime,
                                n_regimes = n_regimes, check_system = check_system,
                                get_system = true, use_intercept = true,
                                verbose = verbose)

        if get_system
            return system
        else
            EE, MM = measurement_error(m)

            lags = n_lags(m)
            YYYY, XXYY, XXXX =
                compute_var_population_moments(data, lags; use_intercept = true)
            out = vecm_approx_state_space(system[:TTT], system[:RRR], system[:QQ],
                                         system[:DD], system[:ZZ], EE, MM, size(data, 1),
                                          n_lags(m), n_cointegrating(m),
                                          n_cointegrating_add(m), DD_coint_add;
                                          get_population_moments = true,
                                          use_intercept = true)

            if get_population_moments
                return out..., YYYY, XXYY, XXXX
            else
                # Compute prior-weighted population moments
                λ = get_λ(m)
                YYYYC = YYYY + λ .* out[1]
                XXYYC = XXYY + λ .* out[2]
                XXXXC = XXXX + λ .* out[3]

                # Draw VECM system
                n_periods = size(data, 2) - lags
                β, Σ =  draw_VECM(YYYYC, XXYYC, XXXXC,
                                  convert(Int, n_periods + λ * n_periods),
                                  size(data, 1), lags, n_cointegrating(m))

                return β, Σ
            end
        end
    end
end

"""
```
compute_system(m::PoolModel{T})
```

Given the current model parameters, compute the state-space system
corresponding to the PoolModel model `m`.

Outputs

```
Φ: state transition function
Ψ: likelihood function, given weights on underlying models (the states) and predictive densities
F_ϵ: structural shock distribution
F_u: likelihood function measurement error distribution
F_λ: initial distribution of λ for state transition function
"""
function compute_system(m::PoolModel{T}; verbose::Symbol = :high,
                        regime_switching::Bool = false,
                        n_regimes::Int = 1) where T<:AbstractFloat
    Φ, F_ϵ, F_λ = transition(m)
    Ψ, F_u = measurement(m)
    return Φ, Ψ, F_ϵ, F_u, F_λ
end

"""
```
compute_system(m::AbstractDSGEModel, system::System;
        observables::Vector{Symbol} = collect(keys(m.observables)),
        pseudo_observables::Vector{Symbol} = collect(keys(m.pseudo_observables)),
        states::Vector{Symbol} = vcat(collect(keys(m.endogenou_states)),
                                      collect(keys(m.endogenous_states_augmented)))
        shocks::Vector{Symbol} = collect(keys(m.exogenous_shocks)),
        zero_DD = false, zero_DD_pseudo = false)
compute_system(m::AbstractDSGEVECMModel, system::System;
        observables::Vector{Symbol} = collect(keys(m.observables)),
        pseudo_observables::Vector{Symbol} = collect(keys(m.pseudo_observables)),
        cointegrating::Vector{Symbol} = collect(keys(m.cointegrating)),
        states::Vector{Symbol} = vcat(collect(keys(m.endogenou_states)),
                                      collect(keys(m.endogenous_states_augmented)))
        shocks::Vector{Symbol} = collect(keys(m.exogenous_shocks)),
        zero_DD = false, zero_DD_pseudo = false)
```
computes the corresponding transition and measurement
equations specified by the keywords (e.g. `states`, `pseudo_observables`)
using the existing `ZZ`, `ZZ_pseudo`, and `TTT` matrices in `system`.

Note that this function does not update the EE matrix, which is
set to all zeros. To incorporate measurement errors, the user
must specify the EE matrix after applying compute_system.

### Keywords
* `observables`: variables that should be
    entered into the new `ZZ` and `DD` matrices.
    The `observables` can be both Observables and PseudoObservables,
    but they must be an element of the system already.
* `pseudo_observables`: variables that should be
    entered into the new `ZZ_pseudo` and `DD_pseudo` matrices.
    The `observables` can be both Observables and PseudoObservables,
    but they must be an element of the system already.
* `cointegrating`: variables that should be
    entered into the new `ZZ` and `DD` matrices as cointegrating relationships.
    The `observables` can be both Observables and PseudoObservables,
    but they must be an element of the system already.
* `states`: variables that should be
    entered into the new `TTT` and `RRR` matrices as states.
    They must be existing states.
* `shocks`: variables that should be
    entered into the new `RRR` and `QQ` matrices as shocks.
    They must be existing exogenous shocks.
"""
function compute_system(m::AbstractDSGEModel{S}, system::System;
                        observables::Vector{Symbol} = collect(keys(m.observables)),
                        pseudo_observables::Vector{Symbol} =
                        collect(keys(m.pseudo_observables)),
                        states::Vector{Symbol} =
                        vcat(collect(keys(m.endogenous_states)),
                             collect(keys(m.endogenous_states_augmented))),
                        shocks::Vector{Symbol} = collect(keys(m.exogenous_shocks)),
                        zero_DD::Bool = false, zero_DD_pseudo::Bool = false,
                        check_system::Bool = true) where {S<:Real}

    # Set up indices
    oid  = m.observables # observables indices dictionary
    pid  = m.pseudo_observables # pseudo observables indices dictionary
    sid  = m.endogenous_states
    said = m.endogenous_states_augmented # states augmented

    # Find shocks to keep
    shock_inds = map(k -> m.exogenous_shocks[k], shocks)
    Qout = system[:QQ][shock_inds, shock_inds]

    # Find states to keep
    if !issubset(states, vcat(collect(keys(sid)), collect(keys(said))))
        false_states = setdiff(states, vcat(collect(keys(sid)), collect(keys(said))))
        error("The following states in keyword `states` do not exist in `system`: " *
              join(string.(false_states), ", "))
    elseif !isempty(setdiff(vcat(collect(keys(sid)), collect(keys(said))), states))
        which_states = Vector{Int}(undef, length(states))
        for i = 1:length(states)
            which_states[i] = haskey(sid, states[i]) ? sid[states[i]] : said[states[i]]
        end
        Tout = system[:TTT][which_states, which_states]
        Rout = system[:RRR][which_states, shock_inds]
        Cout = system[:CCC][which_states]
    else
        which_states = 1:n_states_augmented(m)
        Tout = copy(system[:TTT])
        Rout = system[:RRR][:, shock_inds]
        Cout = copy(system[:CCC])
    end

    # Compute new ZZ and DD matrices if different observables than current system
    if !isempty(symdiff(observables, collect(keys(oid))))
        Zout = zeros(S, length(observables), size(Tout, 1))
        Dout = zeros(S, length(observables))
        for (i, obs) in enumerate(observables)
            Zout[i, :], Dout[i] = if haskey(oid, obs)
                system[:ZZ][oid[obs], which_states], zero_DD ? zero(S) : system[:DD][oid[obs]]
            elseif haskey(pid, obs)
                system[:ZZ_pseudo][pid[obs], which_states], zero_DD ? zero(S) : system[:DD_pseudo][pid[obs]]
            else
                error("Observable/PseudoObservable $obs cannot be found in the DSGE model $m")
            end
        end
    else
        Zout = copy(system[:ZZ])[:, which_states]
        Dout = zero_DD ? zeros(S, size(Zout, 1)) : copy(system[:DD])
    end

    Eout = zeros(S, length(observables), length(observables)) # measurement errors are set to zero

    # Compute new ZZ_pseudo, DD_pseudo if different pseudo_observables than current system
    if !isempty(symdiff(pseudo_observables, collect(keys(pid))))
        Zpseudoout = zeros(S, length(pseudo_observables), size(Tout, 1))
        Dpseudoout = zeros(S, length(pseudo_observables))
        for (i, pseudoobs) in enumerate(pseudo_observables)
            Zpseudoout[i, :], Dpseudoout[i] = if haskey(oid, pseudoobs)
                system[:ZZ][oid[pseudoobs], which_states], zero_DD_pseudo ?
                    zero(S) : system[:DD][oid[pseudoobs]]
            elseif haskey(pid, pseudoobs)
                system[:ZZ_pseudo][pid[pseudoobs], which_states], zero_DD_pseudo ?
                    zero(S) : system[:DD_pseudo][pid[pseudoobs]]
            else
                error("Observable/PseudoObservable $pseudoobs cannot be found in the DSGE model $m")
            end
        end
    else
        Zpseudoout = copy(system[:ZZ_pseudo])[:, which_states]
        Dpseudoout = zero_DD_pseudo ? zeros(S, size(Zpseudoout, 1)) : copy(system[:DD_pseudo])
    end

    if check_system
        @assert size(Zout, 2) == size(Tout, 1) "Dimension 2 of new ZZ ($(size(Zout,2))) does not match dimension of states ($(size(Tout,1)))."
        @assert size(Qout, 1) == size(Rout, 2) "Dimension 2 of new RRR ($(size(Zout,2))) does not match dimension of shocks ($(size(Qout,1)))."
    end

    # Construct new system
    return System(Transition(Tout, Rout, Cout),
                  Measurement(Zout, Dout, Qout, Eout),
                  PseudoMeasurement(Zpseudoout, Dpseudoout))
end

function compute_system(m::AbstractDSGEVECMModel{S}, system::System;
                        observables::Vector{Symbol} = collect(keys(get_observables(get_dsge(m)))),
                        cointegrating::Vector{Symbol} = Vector{Symbol}(undef, 0),
                        cointegrating_add::Vector{Symbol} = Vector{Symbol}(undef, 0),
                        pseudo_observables::Vector{Symbol} =
                        collect(keys(get_dsge(m).pseudo_observables)),
                        states::Vector{Symbol} =
                        vcat(collect(keys(get_dsge(m).endogenous_states)),
                             collect(keys(get_dsge(m).endogenous_states_augmented))),
                        shocks::Vector{Symbol} = collect(keys(get_dsge(m).exogenous_shocks)),
                        zero_DD::Bool = false, zero_DD_pseudo::Bool = false,
                        get_DD_coint_add::Bool = false,
                        check_system::Bool = true) where {S<:Real}
    # Cointegrating relationships should exist as observables/pseudo_observables already
    # in the underlying DSGE. We assume cointegrating relationships come after normal observables.
    # Default behavior is to recreate the underlying DSGE's state space representation, however.
    sys = compute_system(get_dsge(m), system; observables = vcat(observables, cointegrating),
                         pseudo_observables = pseudo_observables,
                         states = states, shocks = shocks, zero_DD = zero_DD,
                         zero_DD_pseudo = zero_DD_pseudo, check_system = check_system)
    if get_DD_coint_add
        mtype = typeof(m)
        DD_coint_add = if hasmethod(compute_DD_coint_add, (mtype, Vector{Symbol}))
            compute_DD_coint_add(m, cointegrating_add)
        elseif hasmethod(compute_DD_coint_add, (mtype, ))
            compute_DD_coint_add(m)
        else
            compute_DD_coint_add(m, sys, cointegrating_add)
        end
        return sys, DD_coint_add
    else
        return sys
    end
end

"""
```
function compute_DD_coint_add(m::AbstractDSGEVECMModel{S}, system::System,
                              cointegrating_add::Vector{Symbol}) where {S <: Real}
```
computes `DD_coint_add` for a `DSGEVECM` model. This vector
holds additional cointegrating relationships that do not require
changes to the `ZZ` matrix.

### Note
We recommend overloading this function if there are
cointegrating relationships which a user does not want
to add to the underlying DSGE. The function `compute_system`
first checks for a method `compute_DD_coint_add` that takes
a Type tuple of `(model_type, Vector{Symbol})` and then `(model_type, )`
before calling this method.

This function is generally intended to be internal. As an example of
other such functions, `eqcond` must be user-defined but
is primarily used internally and not directly called by the user in a script.
"""
function compute_DD_coint_add(m::AbstractDSGEVECMModel{S}, system::System,
                              cointegrating_add::Vector{Symbol}) where {S <: Real}
    if !isempty(cointegrating_add)
        oid = get_observables(get_dsge(m))
        pid = get_pseudo_observables(get_dsge(m))
        Dout = zeros(S, length(cointegrating_add))
        for (i, obs) in enumerate(cointegrating_add)
            Dout[i] = if haskey(oid, obs)
                system[:DD][oid[obs]]
            elseif haskey(pid, obs)
                system[:DD_pseudo][pid[obs]]
            else
                error("Observable/PseudoObservable $obs cannot be found in the DSGE model $m")
            end
        end
        return Dout
    else
        @warn "No additional cointegrating relationships specified. Returning an empty vector."
        return Vector{S}(undef, 0)
    end
end

"""
```
compute_system_function(system::System{S}) where S<:AbstractFloat
```

### Inputs

- `system::System`: The output of compute_system(m), i.e. the matrix outputs from solving a given model, m.

### Outputs

- `Φ::Function`: transition equation
- `Ψ::Function`: measurement equation
- `F_ϵ::Distributions.MvNormal`: shock distribution
- `F_u::Distributions.MvNormal`: measurement error distribution
"""
function compute_system_function(system::System{S}) where S<:AbstractFloat
    # Unpack system
    TTT    = system[:TTT]
    RRR    = system[:RRR]
    CCC    = system[:CCC]
    QQ     = system[:QQ]
    ZZ     = system[:ZZ]
    DD     = system[:DD]
    EE     = system[:EE]

    # Define transition and measurement functions
    @inline Φ(s_t1::Vector{S}, ϵ_t::Vector{S}) = TTT*s_t1 + RRR*ϵ_t + CCC
    @inline Ψ(s_t::Vector{S}) = ZZ*s_t + DD

    # Define shock and measurement error distributions
    nshocks = size(QQ, 1)
    nobs    = size(EE, 1)
    F_ϵ = Distributions.MvNormal(zeros(nshocks), QQ)
    F_u = Distributions.MvNormal(zeros(nobs),    EE)

    return Φ, Ψ, F_ϵ, F_u
end

function zero_system_constants(system::System{S}) where S<:AbstractFloat
    system = copy(system)

    system.transition.CCC = zeros(size(system[:CCC]))
    system.measurement.DD = zeros(size(system[:DD]))
    system.pseudo_measurement.DD_pseudo = zeros(size(system[:DD_pseudo]))

    return system
end

"""
```
var_approx_state_space(TTT, RRR, QQQ, DD, ZZ, EE, MM, p; get_population_moments = false,
    use_intercept = false) where {S<:Real}
```
computes the VAR(p) approximation of the linear state space system

```
sₜ = TTT * sₜ₋₁ + RRR * ϵₜ,
yₜ = ZZ * sₜ + DD + uₜ,
```
where the disturbances are assumed to follow
```
ϵₜ ∼ 𝒩 (0, QQ),
uₜ = ηₜ + MM * ϵₜ,
ηₜ ∼ 𝒩 (0, EE).
```
The `MM` matrix implies
```
cov(ϵₜ, uₜ) = QQ * MM'.
```

### Outputs
If `get_population_moments = false`:
* `β`: VAR(p) coefficients
* `Σ`: innovations variance-covariance matrix for the VAR(p) representation
```
yₜ = Xₜβ + μₜ
```
where `Xₜ` appropriately stacks the constants and `p` lags of `yₜ`, and `μₜ ∼ 𝒩 (0, Σ)`.

If `get_population_moments = true`: we return the limit cross product matrices.
* `yyyyd`: 𝔼[y,y]
* `XXXXd`: 𝔼[y,X(lag rr)]
* `XXyyd`: 𝔼[X(lag rr),X(lag ll)]

Using these matrices, the VAR(p) representation is given by
```
β = XXXXd \\ XXyyd
Σ = yyyyd - XXyyd' * β
```

The keyword `use_intercept` specifies whether or not to use an
intercept term in the VAR approximation.
"""
function var_approx_state_space(TTT::AbstractMatrix{S}, RRR::AbstractMatrix{S},
                                QQ::AbstractMatrix{S}, DD::AbstractVector{S},
                                ZZ::AbstractMatrix{S}, EE::AbstractMatrix{S},
                                MM::AbstractMatrix{S}, p::Int;
                                get_population_moments::Bool = false,
                                use_intercept::Bool = false) where {S<:Real}

    n_obs = size(ZZ, 1)

    HH = EE + MM * QQ * MM'
    VV = QQ * MM'

    ## Compute p autocovariances

    ## Initialize Autocovariances
    GAMM0 = zeros(S, n_obs ^ 2, p + 1)
    GA0 =  solve_discrete_lyapunov(TTT, RRR * QQ * RRR')
    Gl   = ZZ * GA0 * ZZ' + ZZ * RRR * VV + (ZZ * RRR * VV)' + HH
    GAMM0[:, 1] = vec(Gl)
    TTl = copy(TTT)
    GA0ZZ = GA0 * ZZ'
    RRRVV = RRR * VV
    for l = 1:p
        Gl = ZZ * TTl * (GA0ZZ + RRRVV) # ZZ * (TTl * GA0Z) * ZZ' + ZZ * (TTl * RRR * VV)
        GAMM0[:, l+1] = vec(Gl)
        TTl = TTl * TTT
    end

    ## Create limit cross product matrices
    yyyyd = zeros(S, n_obs, n_obs)
    if use_intercept
        XXXXd = zeros(S, 1 + p * n_obs, 1 + p * n_obs)
        yyXXd = zeros(S, n_obs, 1 + p * n_obs)

        XXXXd[1, 1] = 1.
        XXXXd[1, 2:1 + p * n_obs] = repeat(DD', 1, p) # same as kron(ones(1, p), DD')
        XXXXd[2:1 + p * n_obs, 1] = repeat(DD, p)     # same as kron(ones(p), DD)
        yyXXd[:, 1] = DD
    else
        XXXXd = zeros(S, p * n_obs, p * n_obs)
        yyXXd = zeros(S, n_obs, p * n_obs)
    end
    DDDD = DD * DD'
    yyyyd = reshape(GAMM0[:, 1], n_obs, n_obs) + DDDD

    shift = use_intercept ? 1 : 0 # for constructing XXXXd, to select the right indices
    for rr = 1:p
        # 𝔼[yy,x(lag rr)]
        yyXXd[:, n_obs * (rr - 1) + 1 + shift:n_obs * rr + shift] =
            reshape(GAMM0[:, rr + 1], n_obs, n_obs) + DDDD

        # 𝔼[x(lag rr),x(lag ll)]
        for ll = rr:p
            yyyydrrll = reshape(GAMM0[:, ll - rr + 1], n_obs, n_obs) + DDDD
            XXXXd[n_obs * (rr - 1) + 1 + shift:n_obs * rr + shift,
                  n_obs * (ll - 1) + 1 + shift:n_obs * ll + shift] = yyyydrrll
            XXXXd[n_obs * (ll - 1) + 1 + shift:n_obs * ll + shift,
                  n_obs * (rr - 1) + 1 + shift:n_obs * rr + shift] = yyyydrrll'
        end
    end

    XXyyd = convert(Matrix{S}, yyXXd')

    if get_population_moments
        return yyyyd, XXyyd, XXXXd
    else
        β = \(XXXXd, XXyyd)
        Σ = yyyyd - XXyyd' * β
        return β, Σ
    end
end

"""
```
vecm_approx_state_space(TTT, RRR, QQQ, DD, ZZ, EE, MM, n_obs, p, n_coint,
    n_coint, n_coint_add, DD_coint_add; get_population_moments = false,
    use_intercept = false) where {S<:Real}
```
computes the VECM(p) approximation of the linear state space system

```
sₜ = TTT * sₜ₋₁ + RRR * ϵₜ,
yₜ = ZZ * sₜ + DD + uₜ,
```
where the disturbances are assumed to follow
```
ϵₜ ∼ 𝒩 (0, QQ),
uₜ = ηₜ + MM * ϵₜ,
ηₜ ∼ 𝒩 (0, EE).
```
The `MM` matrix implies
```
cov(ϵₜ, uₜ) = QQ * MM'.
```

### Outputs
If `get_population_moments = false`:
* `β`: VECM(p) coefficients. The first `n_coint + n_coint_add`
    coefficients for each observable comprise the error correction terms,
    while the following `1 + p * n_obs` terms are the VAR terms.
* `Σ`: innovations variance-covariance matrix for the VECM(p) representation
```
Δyₜ = eₜβₑ + Xₜβᵥ + μₜ
```
where `βₑ` are the coefficients for the error correction terms;
`eₜ` are the error correction terms specifying the cointegrating relationships;
`βᵥ` are the coefficients for the VAR terms; `Xₜ` appropriately stacks
the constants and `p` lags of `Δyₜ`; and `μₜ ∼ 𝒩 (0, Σ)`.

Note that the error correction terms satisfy the mapping
`eₜ' = C * yₜ₋₁`, where `C` is a matrix.

If `get_population_moments = true`: we return the limit cross product matrices.
* `yyyyd`: 𝔼[y,y]
* `XXXXd`: 𝔼[y,X(lag rr)]
* `XXyyd`: 𝔼[X(lag rr),X(lag ll)]

Note that in the rows of `XXyyd` and the rows and columns of `XXXXd`,
the cointegrating relationships are stacked above the constants and
lagged `Δyₜ`.

Using these matrices, the VAR(p) representation is given by
```
β = XXXXd \\ XXyyd
Σ = yyyyd - XXyyd' * β,
```
where `β` has dimensions `n_obs × (n_coint + n_coint_add + 1 + p * n_obs)`,
and `Σ` has dimensions `n_obs × n_obs`.

The keyword `use_intercept` specifies whether or not to use an
intercept term in the VECM approximation.
"""
function vecm_approx_state_space(TTT::AbstractMatrix{S}, RRR::AbstractMatrix{S},
                                 QQ::AbstractMatrix{S}, DD::AbstractVector{S},
                                 ZZ::AbstractMatrix{S}, EE::AbstractMatrix{S},
                                 MM::AbstractMatrix{S}, n_obs::Int, p::Int,
                                 n_coint::Int, n_coint_add::Int = 0,
                                 DD_coint_add::AbstractVector{S} = Vector{S}(undef, 0);
                                 get_population_moments::Bool = false,
                                 use_intercept::Bool = false,
                                 test_GA0::AbstractMatrix{S} =
                                 Matrix{S}(undef, 0, 0)) where {S <: Real}

    # n_obs is number of observables, n_coint is number of cointegrating relationships
    n_coint_all = n_coint + n_coint_add

    # Create variance-covariance matrices w/measurement error included
    HH = EE + MM * QQ * MM'
    VV = QQ * MM'

    ## Compute p autocovariances

    ## Initialize Autocovariances
    GAMM0 = zeros(S, (n_obs + n_coint) ^ 2, p + 1)
    GA0 = isempty(test_GA0) ? solve_discrete_lyapunov(TTT, RRR * QQ * RRR') : test_GA0 # Matlab code uses a different numerical procedure -> some difference in this matrix
    Gl   = ZZ * GA0 * ZZ' + ZZ * RRR * VV + (ZZ * RRR * VV)' + HH
    GAMM0[:, 1] = vec(Gl)
    TTl = copy(TTT)
    GA0ZZ = GA0 * ZZ'
    RRRVV = RRR * VV
    for l = 1:p
        Gl = ZZ * TTl * (GA0ZZ + RRRVV) # ZZ * (TTl * GA0) * ZZ' + ZZ * (TTl * RRR * VV)
        GAMM0[:, l+1] = vec(Gl)
        TTl = TTl * TTT
    end

    ## Create limit cross product matrices
    DDDD = DD * DD'
    yyyyd_coint0 = reshape(GAMM0[:, 1], n_obs + n_coint, n_obs + n_coint) + DDDD
    yyyyd_coint1 = reshape(GAMM0[:, 2], n_obs + n_coint, n_obs + n_coint) + DDDD
    yyyyd = yyyyd_coint0[1:n_obs, 1:n_obs]

    # n_coint_add are treated as the first set of variables in XX
    # n_coint    are treated as the second set of variables in XX
    # composition: n_coint_add - n_coint - constant - lags
    if use_intercept
        XXXXd = zeros(S, 1 + p * n_obs + n_coint_all, 1 + p * n_obs + n_coint_all)
        yyXXd = zeros(S, n_obs, 1 + p * n_obs + n_coint_all)

        # 𝔼[x(n_coint), x(n_coint)]
        XXXXd[n_coint_add + 1:n_coint_all, n_coint_add + 1:n_coint_all] = yyyyd_coint0[n_obs + 1:n_obs + n_coint, n_obs + 1:n_obs + n_coint]

        # 𝔼[x(const), x(const)]
        XXXXd[n_coint_all + 1, n_coint_all + 1] = 1.

        # 𝔼[x(n_coint), x(const)]
        XXXXd[n_coint_add + 1:n_coint_all, n_coint_all + 1] = DD[n_obs + 1:n_obs + n_coint]
        XXXXd[n_coint_all + 1, n_coint_add + 1:n_coint_all] = DD[n_obs + 1:n_obs + n_coint]

        # 𝔼[x(const), x(lags)]
        XXXXd[n_coint_all + 1, n_coint_all + 2:n_coint_all + 1 + p * n_obs] = repeat(DD[1:n_obs]', 1, p)
        XXXXd[n_coint_all + 2:n_coint_all + 1 + p * n_obs, n_coint_all + 1] =
            XXXXd[n_coint_all + 1, n_coint_all + 2:n_coint_all + 1 + p * n_obs]' # same as kron(ones(p), DD[1:n_obs]) but avoids calculation

        # 𝔼[yy, x(n_coint)]
        yyXXd[:, n_coint_add + 1:n_coint_all] = yyyyd_coint1[1:n_obs, n_obs + 1:n_obs + n_coint]

        # 𝔼[yy, x(const)]
        yyXXd[:, n_coint_all + 1] = DD[1:n_obs]
    else
        XXXXd = zeros(S, p * n_obs + n_coint_all, p * n_obs + n_coint_all)
        yyXXd = zeros(S, n_obs, p * n_obs + n_coint_all)

        # 𝔼[x(n_coint), x(n_coint)]
        XXXXd[n_coint_add + 1:n_coint_all, n_coint_add + 1:n_coint_all] =
            yyyyd_coint0[n_obs + 1:n_obs + n_coint, n_obs + 1:n_obs + n_coint]

        # 𝔼[yy, x(n_coint)]
        yyXXd[:, n_coint_add + 1:n_coint_all] = yyyyd_coint1[1:n_obs, n_obs + 1:n_obs + n_coint]
    end

    if n_coint_add > 0
        DD_coint_add_div2 = DD_coint_add ./ 2
        if use_intercept
            # 𝔼[yy, x(n_coint_add)]
            yyXXd[:, 1:n_coint_add] = DD[1:n_obs] * DD_coint_add_div2

            # 𝔼[x(n_coint_add), x(n_coint_add)]
            XXXXd[1:n_coint_add, 1:n_coint_add] = DD_coint_add * DD_coint_add' ./ 3

            # 𝔼[x(n_coint_add), x(n_coint)]
            XXXXd[1:n_coint_add, 1 + n_coint_add:n_coint_all] =
                DD_coint_add_div2 * DD[n_obs + 1: n_obs + n_coint]'
            XXXXd[1 + n_coint_add:n_coint_all, 1:n_coint_add] =
                XXXXd[1:n_coint_add, 1 + n_coint_add:n_coint_all]' # transpose of the previous line

            # 𝔼[x(n_coint_add), x(const)]
            XXXXd[1:n_coint_add, n_coint_all + 1] = DD_coint_add_div2
            XXXXd[n_coint_all + 1, 1:n_coint_add] = DD_coint_add_div2'
        else
            # 𝔼[yy, x(n_coint_add)]
            yyXXd[:, 1:n_coint_add] = DD[1:n_obs] * DD_coint_add_div2

            # 𝔼[x(n_coint_add), x(n_coint_add)]
            XXXXd[1:n_coint_add, 1:n_coint_add] = DD_coint_add * DD_coint_add' ./ 3

            # 𝔼[x(n_coint_add), x(n_coint)]
            XXXXd[1:n_coint_add, 1 + n_coint_add:n_coint_all] =
                DD_coint_add_div2 * DD[n_obs + 1: n_obs + n_coint]'
            XXXXd[1 + n_coint_add:n_coint_all, 1:n_coint_add] =
                XXXXd[1:n_coint_add, 1 + n_coint_add:n_coint_all]' # transpose of the previous line
        end
    end

    shift = use_intercept ? 1 : 0 # for constructing XXXXd, to select the right indices
    for rr = 1:p
        # 𝔼[yy, x(lag rr)]
        yyyyd_cointrr = reshape(GAMM0[:, rr + 1], n_obs + n_coint, n_obs + n_coint) + DDDD
        yyXXd[:, n_coint_all + 1 + n_obs * (rr - 1) + shift:n_coint_all + n_obs * rr + shift] =
            yyyyd_cointrr[1:n_obs, 1:n_obs]

        if n_coint_add > 0
            # 𝔼[x(n_coint_add), x(lag rr)]
            XXXXd[1:n_coint_add, n_coint_all + 1 + n_obs * (rr - 1) + shift:n_coint_all + n_obs * rr + shift] =
                DD_coint_add_div2 * DD[1:n_obs]'
            XXXXd[n_coint_all + 1 + n_obs * (rr - 1) + shift:n_coint_all + n_obs * rr + shift, 1:n_coint_add] =
                XXXXd[1:n_coint_add, n_coint_all + 1 + n_obs * (rr - 1) + shift:n_coint_all + n_obs * rr + shift]'
        end

        # 𝔼[x(n_coint), x(lag rr)]
        yyyyd_cointrr1 = reshape(GAMM0[:, rr], n_obs + n_coint, n_obs + n_coint) + DDDD
        XXXXd[n_coint_add + 1:n_coint_all, n_coint_all + 1 + n_obs * (rr - 1) + shift:n_coint_all + n_obs * rr + shift] =
            yyyyd_cointrr1[n_obs + 1:n_obs + n_coint, 1:n_obs]
        XXXXd[n_coint_all + 1 + n_obs * (rr - 1) + shift:n_coint_all + n_obs * rr + shift, n_coint_add + 1:n_coint_all] =
            yyyyd_cointrr1[n_obs + 1:n_obs + n_coint, 1:n_obs]'

        # 𝔼[x(lag rr), x(lag ll)]
        for ll = rr:p
            yyyyd_cointrrll = reshape(GAMM0[:, ll - rr + 1], n_obs + n_coint, n_obs + n_coint) + DDDD
            XXXXd[n_coint_all + 1 + n_obs * (rr - 1) + shift:n_coint_all + n_obs * rr + shift,
                  n_coint_all + 1 + n_obs * (ll - 1) + shift:n_coint_all + n_obs * ll + shift] = yyyyd_cointrrll[1:n_obs, 1:n_obs]
            XXXXd[n_coint_all + 1 + n_obs * (ll - 1) + shift:n_coint_all + n_obs * ll + shift,
                  n_coint_all + 1 + n_obs * (rr - 1) + shift:n_coint_all + n_obs * rr + shift] = yyyyd_cointrrll[1:n_obs, 1:n_obs]'
        end
    end

    XXyyd = convert(Matrix{S}, yyXXd')

    if get_population_moments
        return yyyyd, XXyyd, XXXXd
    else
        β = XXXXd \ XXyyd
        Σ = yyyyd - XXyyd' * β
        return β, Σ
    end
end


text:
type Process
  id::Int32
  participant::Bool
  leader_process_id::Int32
  next_process::Process
  Process(id) = (x = new(); x.id = id; x.participant = false; x.leader_process_id = 0; x.next_process = x)
end

function send_election_message(p::Process)
  handle_election_message(p.next_process, p.id)
end

function handle_election_message(p::Process, msg::Int32)
  if msg > p.id
    handle_election_message(p.next_process, msg)
  elseif msg < p.id && p.participant == false
    send_election_message(p)
  elseif msg == p.id
    p.leader_process_id = p.id
    send_elected_message(p)
  end
end

function send_elected_message(p::Process)
  p.participant = false
  handle_elected_message(p.next_process, p.id)
end

function handle_elected_message(p::Process, msg::Int32)
  if msg == p.id
    return
  end

  p.participant = false
  p.leader_process_id = msg
  handle_elected_message(p.next_process, msg)
end

function elect_leader(p::Process)
  send_election_message(p)
end

function create_ring(n)
  uids = shuffle(1:n)

  ring = Process[]
  for i = uids
    p = Process(i)
    push!(ring, p)
  end

  for i = 1:n-1
    ring[i].next_process = ring[i+1]
  end
  ring[end].next_process = ring[1]

  ring
end

# This is the implementation of a Leader Election algorithm in a ring
# proposed by Chang and Roberts (1979) (https://en.wikipedia.org/wiki/Chang_and_Roberts_algorithm)
# The final outcome is, that the process with the highest UID is elected as a leader.
# The processes do not require knowledge about the amount of processes in the ring

ring = create_ring(10)
elect_leader(ring[1])
for r = ring
  @printf("For %d the leader is: %d\n", r.id, r.leader_process_id)
end

max_stars_repo_name: CarbonDDR/al-go-rithms


text:
@sumProductRule(:node_type     => Categorical,
                :outbound_type => Message{Categorical},
                :inbound_types => (Nothing, Message{PointMass}),
                :name          => SPCategoricalOutNP)

@sumProductRule(:node_type     => Categorical,
                :outbound_type => Message{Dirichlet},
                :inbound_types => (Message{PointMass}, Nothing),
                :name          => SPCategoricalIn1PN)

@naiveVariationalRule(:node_type     => Categorical,
                      :outbound_type => Message{Categorical},
                      :inbound_types => (Nothing, ProbabilityDistribution),
                      :name          => VBCategoricalOut)

@naiveVariationalRule(:node_type     => Categorical,
                      :outbound_type => Message{Dirichlet},
                      :inbound_types => (ProbabilityDistribution, Nothing),
                      :name          => VBCategoricalIn1)