text:
module NAFEMS_R0031_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using Statistics: mean
using LinearAlgebra: norm, cross
function NAFEMS_R0031_1()
    println("""
    Laminated Strip Under Three-Point Bending
    """)

    # Determine the central transverse displacement in a simply-supported seven
    # layer symmetric strip with a central line load. A 0/90/0/90/0/90/0
    # material lay-up is specified with the center ply being four times as
    # thick as the others.
    # Reference: NAFEMS Report R0031, Test No.1, 17-Dec-1998.

    # Because of the symmetries of the geometry and load, only the
    # first-quadrant   (in XY) quarter of the plate is modeled.

    # The coordinate system is centered at point E (at the difference with
    # respect to the original benchmark definition).  The  load is applied
    # along a curve passing through point C. The simple support is applied
    # along the curve passing through point B.

    # We realize the simple supports along the lines  A, B and the line load at
    # point C  are illegal from the point of view of convergence.  No
    # convergence can be hoped for as the stress underneath the load and above
    # the simple supports  is infinite in the limit (these locations are stress
    # singularities).   However, for relatively coarse meshes the results away
    # from the singularities are still meaningful.

    # The target quantities are displacement at the bottom surface at point E,
    # the tensile axial stress at the same point,  and of the transverse shear
    # stress at point D  in between the bottommost two layers (See figure 1).

    t0 = time()
    # Orthotropic material parameters of the material of the layers
    E1s = 100.0*phun("GPa")
    E2s = E3s = 5.0*phun("GPa")
    nu12s = 0.4
    nu13s = 0.3
    nu23s = 0.3
    G12s = 3.0*phun("GPa")
    G13s = G23s = 2.0*phun("GPa")
    CTE1 = 3.0e-6
    CTE2 = 2.0e-5
    CTE3 = 2.0e-5

    AB = 30.0*phun("mm") # span between simple supports
    CD = 4.0*phun("mm") # distance between the point D and the center
    OH = 10.0*phun("mm") # overhang
    W = 10.0*phun("mm") # width of the strip

    # Here we define the layout and the thicknesses of the layers.
    angles = vec([0 90 0 90 0 90 0]);
    ts = vec([0.1  0.1  0.1  0.4  0.1  0.1  0.1])*phun("mm"); # layer thicknesses
    TH = sum(ts); # total thickness of the plate

    tolerance = 0.0001*TH

    # The line load is in the negative Z direction.
    q0 = 10*phun("N/mm"); #    line load

    # Reference deflection under the load is
    wEref = -1.06*phun("mm");

    # The reference tensile stress at the bottom of the lowest layer is
    sigma11Eref = 684*phun("MPa");

    # Because we model the first-quadrant quarter of the plate using
    # coordinate axes centered  at the point E  the shear at the point D is
    # positive instead of negative as in the benchmark where the coordinate
    # system is located at the outer corner of the strip.
    sigma13Dref=4.1*phun("MPa");

    Refinement = 10
    # We select 8 elements spanwise and 2 elements widthwise.  The overhang
    # of the plate is given one element.
    nL = Refinement * 4; nO = Refinement * 1; nW = Refinement * 1;

    # Each layer is modeled with a single element.
    nts= Refinement * ones(Int, length(angles));# number of elements per layer

    xs = unique(vcat(collect(linearspace(0.0,AB/2,nL+1)), [CD],
    collect(linearspace(AB/2,AB/2+OH,nO+1))))
    xs = xs[sortperm(xs)]
    ys = collect(linearspace(0.0,W/2,nW+1));

    fens,fes = H8layeredplatex(xs, ys, ts, nts)

    # This is the material  model
    MR = DeforModelRed3D
    material = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    CTE1, CTE2, CTE3)

    # The material coordinate system function is defined as:
    function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        rotmat3!(csmatout, angles[fe_label]/180.0*pi* [0.0; 0.0; 1.0]);
    end

    # The vvolume integrals are evaluated using this rule
    gr = GaussRule(3, 2)

    # We will create two regions, one for the layers with 0°  orientation,
    # and one for the layers with 90° orientation.
    rl1 = vcat(selectelem(fens, fes, label = 1), selectelem(fens, fes, label = 3),
    selectelem(fens, fes, label = 5), selectelem(fens, fes, label = 7))
    rl2 = vcat(selectelem(fens, fes, label = 2), selectelem(fens, fes, label = 4),
    selectelem(fens, fes, label = 6))
    region1 = FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegDomain(subset(fes, rl1), gr), CSys(3, 3, updatecs!), material))
    region2 = FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegDomain(subset(fes, rl2), gr), CSys(3, 3, updatecs!), material))

    # File =  "NAFEMS-R0031-1-plate-r1.vtk"
    # vtkexportmesh(File, region1["femm"].integdomain.fes.conn, fens.xyz, FinEtools.MeshExportModule.H8)
    # # @async run(`"paraview.exe" $File`)
    # File =  "NAFEMS-R0031-1-plate-r2.vtk"
    # vtkexportmesh(File, region2["femm"].integdomain.fes.conn, fens.xyz, FinEtools.MeshExportModule.H8)
    # @async run(`"paraview.exe" $File`)

    # The essential boundary conditions are applied on the symmetry planes.
    # First the plane X=0;...
    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    ex0 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
    # ... and then the plane Y=0.
    ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    ey0 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>ly0 )
    # The transverse displacement is fixed along the line  passing through
    # point B. The nodes are fixed in the box along this line in the Z
    # direction.
    lz0 = selectnode(fens, box=[AB/2 AB/2 -Inf Inf -Inf Inf], inflate=tolerance)
    ez0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lz0 )

    # The traction boundary condition is applied  along  the edge of the
    # mesh passing through point C at the top surface of the strip.   First
    # we extract the boundary of the hexahedral mesh.
    bfes = meshboundary(fes)
    # From  the entire boundary we select those quadrilaterals that lie on the plane
    # X = 0
    xl = selectelem(fens, bfes, box = [0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    # Now we extract the boundary  of these selected quadrilaterals
    bbfes = meshboundary(subset(bfes, xl))
    # …  And from these  we extract the ones at the top
    zl = selectelem(fens, bbfes, box = [0.0 0.0 -Inf Inf TH TH], inflate=tolerance)
    # Note that  we have to apply only half of the line load given that
    # were modeling  just one quarter of the geometry and we are splitting
    # the line load  with the symmetry plane X=0. Also note that the
    # quadrature rule is one-dimensional  since we are integrating along
    # a curve.
    Trac = FDataDict("traction_vector"=>vec([0.0; 0.0; -q0/2]), "femm"=>FEMMBase(IntegDomain(subset(bbfes, zl), GaussRule(1, 3))))

    modeldata = FDataDict("fens"=>fens, "regions"=>[region1, region2], "essential_bcs"=>[ex0, ey0, ez0], "traction_bcs"=> [Trac])
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    modeldata["postprocessing"] = FDataDict("file"=>"NAFEMS-R0031-1-plate")
    modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

    u = modeldata["u"]
    geom = modeldata["geom"]

    # The results of the displacement and stresses will be reported at
    # nodes located at the appropriate points.
    nE = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
    nC = selectnode(fens, box=[0.0 0.0 0.0 0.0 TH TH], inflate=tolerance)
    nD = selectnode(fens, box=[CD CD 0.0 0.0 ts[1] ts[1]], inflate=tolerance)
    n0z = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 TH], inflate=tolerance)
    ix = sortperm(geom.values[n0z, 3])
    # println("ix = $(ix)")

    cdis = mean(u.values[nE, 3])
    println("")
    println("Normalized Center deflection: $(cdis/wEref)")

    extrap = :extraptrend
    # # extrap = :extrapmean
    inspectormeth = :averaging
    # extrap = :default
    # inspectormeth = :invdistance

    modeldata["postprocessing"] = FDataDict("file"=>"NAFEMS-R0031-1-plate-sx",
    "quantity"=>:Cauchy, "component"=>1, "outputcsys"=>CSys(3),
    "nodevalmethod"=>inspectormeth, "reportat"=>extrap)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    s = modeldata["postprocessing"]["exported"][1]["field"]
    println("sx@E = $(s.values[nE]/phun("MPa")) [MPa]")

    modeldata["postprocessing"] = FDataDict("file"=>"NAFEMS-R0031-1-plate-sxz",
    "quantity"=>:Cauchy, "component"=>5, "outputcsys"=>CSys(3),
    "nodevalmethod"=>inspectormeth, "reportat"=>extrap)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    s = modeldata["postprocessing"]["exported"][1]["field"]
    println("sxz@D_1 = $(s.values[nD]/phun("MPa")) [MPa]")
    s = modeldata["postprocessing"]["exported"][2]["field"]
    println("sxz@D_2 = $(s.values[nD]/phun("MPa")) [MPa]")

    #
    # s = fieldfromintegpoints(region1["femm"], geom, u, :Cauchy, 1;
    #     outputcsys = CSys(3), nodevalmethod = inspectormeth, reportat = extrap)
    # println("sx@E = $(s.values[nE]/phun("MPa")) [MPa]")
    # sx_z = s.values[n0z]/phun("MPa")
    # println("sx(z)_1 = $(sx_z)")
    #
    # s = fieldfromintegpoints(region1["femm"], geom, u, :Cauchy, 5;
    #     outputcsys = CSys(3), nodevalmethod = inspectormeth, reportat = extrap)
    # println("sxz@D_1 = $(s.values[nD]/phun("MPa")) [MPa]")
    # sxz_z_1 = s.values[n0z]/phun("MPa")
    # println("sxz(z)_1 = $(sxz_z_1)")
    # s = fieldfromintegpoints(region2["femm"], geom, u, :Cauchy, 5;
    #     outputcsys = CSys(3), nodevalmethod = inspectormeth, reportat = extrap)
    # println("sxz@D_2 = $(s.values[nD]/phun("MPa")) [MPa]")
    # sxz_z_2 = s.values[n0z]/phun("MPa")
    # println("sxz(z)_2 = $(sxz_z_2)")

    # function _inspector(idat, elnum, conn, xe,  out,  xq)
    #     # xe = coordinates of the nodes of the element
    #     # xq = coordinate of the quadrature point
    #     println("@$(xq): $(out/1.0e6)")
    #     return idat
    # end
    #
    # felist = selectelem(fens, region1["femm"].integdomain.fes,
    #     box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance, allin = false)
    #
    # inspectintegpoints(region1["femm"], geom, u, felist,
    #     _inspector, 0, quantity=:Cauchy, outputcsys = CSys(3))
    #
    # femm = deepcopy(region1["femm"])
    # femm.integdomain.fes = subset(femm.integdomain.fes, felist)
    # associategeometry!(femm, geom)
    # s = fieldfromintegpoints(femm, geom, u, :Cauchy, 5;
    #     outputcsys = CSys(3), nodevalmethod = inspectormeth, reportat = extrap)
    # println("sxz@D_1 = $(s.values[nD]/phun("MPa")) [MPa]")

    # felist = selectelem(fens, region2["femm"].integdomain.fes,
    #     box=[0.0 0.0 0.0 0.0 0.0 TH], inflate=tolerance, allin = false)
    #
    # inspectintegpoints(region2["femm"], geom, u, felist,
    #     _inspector, 0, quantity=:Cauchy, outputcsys = CSys(3))

    println("Done")
    true

end # NAFEMS_R0031_1

function NAFEMS_R0031_1_H20()
    println("""
    Laminated Strip Under Three-Point Bending
    """)

    # Determine the central transverse displacement in a simply-supported seven
    # layer symmetric strip with a central line load. A 0/90/0/90/0/90/0
    # material lay-up is specified with the center ply being four times as
    # thick as the others.
    # Reference: NAFEMS Report R0031, Test No.1, 17-Dec-1998.

    # Because of the symmetries of the geometry and load, only the
    # first-quadrant   (in XY) quarter of the plate is modeled.

    # The coordinate system is centered at point E (at the difference with
    # respect to the original benchmark definition).  The  load is applied
    # along a curve passing through point C. The simple support is applied
    # along the curve passing through point B.

    # We realize the simple supports along the lines  A, B and the line load at
    # point C  are illegal from the point of view of convergence.  No
    # convergence can be hoped for as the stress underneath the load and above
    # the simple supports  is infinite in the limit (these locations are stress
    # singularities).   However, for relatively coarse meshes the results away
    # from the singularities are still meaningful.

    # The target quantities are displacement at the bottom surface at point E,
    # the tensile axial stress at the same point,  and of the transverse shear
    # stress at point D  in between the bottommost two layers (See figure 1).

    t0 = time()
    # Orthotropic material parameters of the material of the layers
    E1s = 100.0*phun("GPa")
    E2s = E3s = 5.0*phun("GPa")
    nu12s = 0.4
    nu13s = 0.3
    nu23s = 0.3
    G12s = 3.0*phun("GPa")
    G13s = G23s = 2.0*phun("GPa")
    CTE1 = 3.0e-6
    CTE2 = 2.0e-5
    CTE3 = 2.0e-5

    AB = 30.0*phun("mm") # span between simple supports
    CD = 4.0*phun("mm") # distance between the point D and the center
    OH = 10.0*phun("mm") # overhang
    W = 10.0*phun("mm") # width of the strip

    # Here we define the layout and the thicknesses of the layers.
    angles = vec([0 90 0 90 0 90 0]);
    ts = vec([0.1  0.1  0.1  0.4  0.1  0.1  0.1])*phun("mm"); # layer thicknesses
    TH = sum(ts); # total thickness of the plate

    tolerance = 0.0001*TH

    # The line load is in the negative Z direction.
    q0 = 10*phun("N/mm"); #    line load

    # Reference deflection under the load is
    wEref = -1.06*phun("mm");

    # The reference tensile stress at the bottom of the lowest layer is
    sigma11Eref = 684*phun("MPa");

    # Because we model the first-quadrant quarter of the plate using
    # coordinate axes centered  at the point E  the shear at the point D is
    # positive instead of negative as in the benchmark where the coordinate
    # system is located at the outer corner of the strip.
    sigma13Dref=4.1*phun("MPa");

    Refinement = 8
    # We select 8 elements spanwise and 2 elements widthwise.  The overhang
    # of the plate is given one element.
    nL = Refinement * 4; nO = Refinement * 1; nW = Refinement * 1;

    # Each layer is modeled with a single element.
    nts= Refinement * ones(Int, length(angles));# number of elements per layer

    xs = unique(vcat(collect(linearspace(0.0,AB/2,nL+1)), [CD],
    collect(linearspace(AB/2,AB/2+OH,nO+1))))
    xs = xs[sortperm(xs)]
    ys = collect(linearspace(0.0,W/2,nW+1));

    fens,fes = H8layeredplatex(xs, ys, ts, nts)
    fens,fes = H8toH20(fens,fes)

    # This is the material  model
    MR = DeforModelRed3D
    material = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    CTE1, CTE2, CTE3)

    # The material coordinate system function is defined as:
    function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        rotmat3!(csmatout, angles[fe_label]/180.0*pi* [0.0; 0.0; 1.0]);
    end

    # The vvolume integrals are evaluated using this rule
    gr = GaussRule(3, 3)

    # We will create two regions, one for the layers with 0°  orientation,
    # and one for the layers with 90° orientation.
    rl1 = vcat(selectelem(fens, fes, label = 1), selectelem(fens, fes, label = 3),
    selectelem(fens, fes, label = 5), selectelem(fens, fes, label = 7))
    rl2 = vcat(selectelem(fens, fes, label = 2), selectelem(fens, fes, label = 4),
    selectelem(fens, fes, label = 6))
    region1 = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rl1), gr), CSys(3, 3, updatecs!), material))
    region2 = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rl2), gr), CSys(3, 3, updatecs!), material))

    # File =  "NAFEMS-R0031-1-plate-r1.vtk"
    # vtkexportmesh(File, region1["femm"].integdomain.fes.conn, fens.xyz, FinEtools.MeshExportModule.H8)
    # # @async run(`"paraview.exe" $File`)
    # File =  "NAFEMS-R0031-1-plate-r2.vtk"
    # vtkexportmesh(File, region2["femm"].integdomain.fes.conn, fens.xyz, FinEtools.MeshExportModule.H8)
    # @async run(`"paraview.exe" $File`)

    # The essential boundary conditions are applied on the symmetry planes.
    # First the plane X=0;...
    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    ex0 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
    # ... and then the plane Y=0.
    ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    ey0 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>ly0 )
    # The transverse displacement is fixed along the line  passing through
    # point B. The nodes are fixed in the box along this line in the Z
    # direction.
    lz0 = selectnode(fens, box=[AB/2 AB/2 -Inf Inf 0.0 0.0], inflate=tolerance)
    ez0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lz0 )

    # The traction boundary condition is applied  along  the edge of the
    # mesh passing through point C at the top surface of the strip.   First
    # we extract the boundary of the hexahedral mesh.
    bfes = meshboundary(fes)
    # From  the entire boundary we select those quadrilaterals that lie on the plane
    # X = 0
    xl = selectelem(fens, bfes, box = [0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    # Now we extract the boundary  of these selected quadrilaterals
    bbfes = meshboundary(subset(bfes, xl))
    # …  And from these  we extract the ones at the top
    zl = selectelem(fens, bbfes, box = [0.0 0.0 -Inf Inf TH TH], inflate=tolerance)
    # Note that  we have to apply only half of the line load given that
    # were modeling  just one quarter of the geometry and we are splitting
    # the line load  with the symmetry plane X=0. Also note that the
    # quadrature rule is one-dimensional  since we are integrating along
    # a curve.
    Trac = FDataDict("traction_vector"=>vec([0.0; 0.0; -q0/2]),    "femm"=>FEMMBase(IntegDomain(subset(bbfes, zl), GaussRule(1, 3))))

    modeldata = FDataDict("fens"=>fens, "regions"=>[region1, region2], "essential_bcs"=>[ex0, ey0, ez0], "traction_bcs"=> [Trac])
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    modeldata["postprocessing"] = FDataDict("file"=>"NAFEMS-R0031-1-plate")
    modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

    u = modeldata["u"]
    geom = modeldata["geom"]

    # The results of the displacement and stresses will be reported at
    # nodes located at the appropriate points.
    nE = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
    nC = selectnode(fens, box=[0.0 0.0 0.0 0.0 TH TH], inflate=tolerance)
    nD = selectnode(fens, box=[CD CD 0.0 0.0 ts[1] ts[1]], inflate=tolerance)

    cdis = mean(u.values[nE, 3])
    println("")
    println("Normalized Center deflection: $(cdis/wEref)")

    modeldata["postprocessing"] = FDataDict("file"=>"NAFEMS-R0031-1-plate-sx",
    "quantity"=>:Cauchy, "component"=>1, "outputcsys"=>CSys(3))
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    s = modeldata["postprocessing"]["exported"][1]["field"]
    println("sx@E = $(s.values[nE]/phun("MPa")) [MPa]")

    modeldata["postprocessing"] = FDataDict("file"=>"NAFEMS-R0031-1-plate-sxz",
    "quantity"=>:Cauchy, "component"=>5, "outputcsys"=>CSys(3))
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    s = modeldata["postprocessing"]["exported"][1]["field"]
    println("sxz@D = $(s.values[nD]/phun("MPa")) [MPa]")
    s = modeldata["postprocessing"]["exported"][2]["field"]
    println("sxz@D = $(s.values[nD]/phun("MPa")) [MPa]")

    println("Done")
    true

end # NAFEMS_R0031_1_H20

function NAFEMS_R0031_2_both()
    println("""
    R0031(2): Wrapped thick cylinder under pressure and thermal loading.
    """)

    t0 = time()
    pu = ustring -> phun(ustring; system_of_units = :SIMM)
    # Orthotropic material parameters of the outer cylinder
    E1s = 130.0*pu("GPa")
    E2s = 5.0*pu("GPa")
    E3s = E2s
    nu12s = nu13s = 0.25
    nu23s = 0.0
    G12s = 10.0*pu("GPa")
    G13s = G12s
    G23s = 5.0*pu("GPa")
    CTE1 = 3.0e-6
    CTE2 = 2.0e-5
    CTE3 = 2.0e-5
    # Isotropic material parameters of the inner cylinder
    E = 2.1e5*pu("MPa")
    nu = 0.3
    CTE = 2.0e-5

    L = 200.0*pu("mm") # length of the cylinder
    ri = 23.0*pu("mm") # inner radius of the cylinder
    ti = 2.0*pu("mm") # thickness of the inner cylinder
    te = 2.0*pu("mm") # thickness of the outer cylinder
    q0 = 200.0*pu("MPa") # inner pressure
    dT = 130*pu("K") # temperature rise

    tolerance = 0.0001*ti

    # Generate mesh
    nL = 16 # number of elements lengthwise
    nc = 8 # number of elements circumferentially
    xs = collect(linearspace(0.0, L/2, nL+1))
    ys = collect(linearspace(0.0, pi/2, nc+1))
    ts = [ti; te];# layer thicknesses
    nts= 6*ones(Int, length(ts));# number of elements per layer
    fens,fes = H8layeredplatex(xs, ys, ts, nts)
    fens,fes = H8toH20(fens,fes)
    bfes = meshboundary(fes)
    # inner surface  for the pressure loading
    intl = selectelem(fens, bfes; facing=true, direction = [0.0 0.0 -1.0])
    # Shape into a cylinder
    for i = 1:count(fens)
        z = fens.xyz[i,1]; a = fens.xyz[i,2]; t = fens.xyz[i,3];
        fens.xyz[i,:] = [(ri+t)*cos(pi/2-a) (ri+t)*sin(pi/2-a) z];
    end

    MR = DeforModelRed3D
    outermaterial = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    CTE1, CTE2, CTE3)
    innermaterial = MatDeforElastIso(MR,
    0.0, E, nu, CTE)

    function cylcs!(csmatout::FFltMat, XYZ::FFltMat)
        csmatout[:, 2] = [0.0 0.0 1.0]
        csmatout[:, 3] = XYZ
        csmatout[3, 3] = 0.0
        csmatout[:, 3] = csmatout[:, 3]/norm(csmatout[:, 3])
        csmatout[:, 1] = cross(csmatout[:, 2], csmatout[:, 3])
    end

    function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        cylcs!(csmatout, XYZ)
    end

    gr = GaussRule(3, 3)

    rli = selectelem(fens, fes, label=1)
    innerregion = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rli), gr), innermaterial))
    rle = selectelem(fens, fes, label=2)
    outerregion = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rle), gr), CSys(3, 3, updatecs!), outermaterial))

    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    lz0 = selectnode(fens, box=[-Inf Inf -Inf Inf 0.0 0.0], inflate=tolerance)

    ex0 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
    ey0 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>ly0 )
    ez0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lz0 )

    function getpr!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        csmatout = zeros(3, 3)
        cylcs!(csmatout, XYZ)
        copy!(forceout, q0*csmatout[:, 3])
    end

    Trac = FDataDict("traction_vector"=>getpr!, "femm"=>FEMMBase(IntegDomain(subset(bfes, intl), GaussRule(2, 3))))

    modeldata = FDataDict("fens"=>fens, "regions"=>[innerregion, outerregion], "essential_bcs"=>[ex0, ey0, ez0], "traction_bcs"=>[Trac], "temperature_change"=>FDataDict("temperature"=>dT))
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    u = modeldata["u"]
    geom = modeldata["geom"]
    # File =  "NAFEMS-R0031-2-plate.vtk"
    # vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H20;
    #     scalars = [("Layer", fes.label)], vectors = [("displacement", u.values)])
    # @async run(`"paraview.exe" $File`)

    modeldata["postprocessing"] = FDataDict("file"=>"NAFEMS-R0031-2", "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy, "component"=>6)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    modeldata["postprocessing"] = FDataDict("file"=>"NAFEMS-R0031-2-elem", "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy, "component"=>6)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)

    println("Done")
    true

end # NAFEMS_R0031_2_both

function NAFEMS_R0031_2_pressure()
    println("""
    R0031(2): Wrapped thick cylinder under pressure and thermal loading.
    """)

    t0 = time()
    # Orthotropic material parameters of the outer cylinder
    E1s = 130.0*phun("GPa")
    E2s = 5.0*phun("GPa")
    E3s = E2s
    nu12s = nu13s = 0.25
    nu23s = 0.0
    G12s = 10.0*phun("GPa")
    G13s = G12s
    G23s = 5.0*phun("GPa")
    CTE1 = 3.0e-6
    CTE2 = 2.0e-5
    CTE3 = 2.0e-5
    # Isotropic material parameters of the inner cylinder
    E = 2.1e5*phun("MPa")
    nu = 0.3
    CTE = 2.0e-5

    L = 200.0*phun("mm") # length of the cylinder
    ri = 23.0*phun("mm") # inner radius of the cylinder
    ti = 2.0*phun("mm") # thickness of the inner cylinder
    te = 2.0*phun("mm") # thickness of the outer cylinder
    q0 = 200.0*phun("MPa") # inner pressure
    dT = 0*phun("K") # NO temperature rise

    tolerance = 0.0001*ti

    # Generate mesh
    nL = 18 # number of elements lengthwise
    nc = 18 # number of elements circumferentially
    xs = collect(linearspace(0.0, L/2, nL+1))
    ys = collect(linearspace(0.0, pi/2, nc+1))
    ts = [ti; te];# layer thicknesses
    nts= 6*ones(Int, length(ts));# number of elements per layer
    fens,fes = H8layeredplatex(xs, ys, ts, nts)
    fens,fes = H8toH20(fens,fes)
    bfes = meshboundary(fes)
    # inner surface  for the pressure loading
    intl = selectelem(fens, bfes; facing=true, direction = [0.0 0.0 -1.0])
    # Shape into a cylinder
    for i = 1:count(fens)
        z = fens.xyz[i,1]; a = fens.xyz[i,2]; t = fens.xyz[i,3];
        fens.xyz[i,:] = [(ri+t)*cos(pi/2-a) (ri+t)*sin(pi/2-a) z];
    end

    MR = DeforModelRed3D
    outermaterial = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    CTE1, CTE2, CTE3)
    innermaterial = MatDeforElastIso(MR,
    0.0, E, nu, CTE)

    function cylcs!(csmatout::FFltMat, XYZ::FFltMat)
        csmatout[:, 2] = [0.0 0.0 1.0]
        radial = XYZ; radial[3] = 0.0
        csmatout[:, 3] = radial/norm(radial)
        csmatout[:, 1] = cross(csmatout[:, 2], csmatout[:, 3])
    end

    function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        cylcs!(csmatout, XYZ)
    end

    gr = GaussRule(3, 3)

    rli = selectelem(fens, fes, label=1)
    innerregion = FDataDict("femm"=>FEMMDeforLinear(MR,
    IntegDomain(subset(fes, rli), gr), innermaterial))
    rle = selectelem(fens, fes, label=2)
    outerregion = FDataDict("femm"=>FEMMDeforLinear(MR,
    IntegDomain(subset(fes, rle), gr), CSys(3, 3, updatecs!), outermaterial))

    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    lz0 = selectnode(fens, box=[-Inf Inf -Inf Inf 0.0 0.0], inflate=tolerance)

    ex0 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
    ey0 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>ly0 )
    ez0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lz0 )

    function getpr!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        csmatout = zeros(3, 3)
        cylcs!(csmatout, XYZ)
        copy!(forceout, q0*csmatout[:, 3])
    end

    Trac = FDataDict("traction_vector"=>getpr!, "femm"=>FEMMBase(IntegDomain(subset(bfes, intl), GaussRule(2, 3))))

    modeldata = FDataDict("fens"=>fens,
    "regions"=>[innerregion, outerregion],
    "essential_bcs"=>[ex0, ey0, ez0],
    "traction_bcs"=>[Trac],
    "temperature_change"=>FDataDict("temperature"=>dT)
    )
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    u = modeldata["u"]
    geom = modeldata["geom"]
    # lcenter = selectnode(fens, box=[a/2 a/2  b/2 b/2 -Inf Inf], inflate=tolerance)
    # cdis = abs(mean(u.values[lcenter, 3]))
    # println("")
    # println("Normalized Center deflection: $(cdis/wc_analytical)")

    File =  "NAFEMS-R0031-2-plate.vtk"
    vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H20; scalars = [("Layer", fes.label)], vectors = [("displacement", u.values)])
    @async run(`"paraview.exe" $File`)

    println("Done")
    true

end # NAFEMS_R0031_2_pressure

function NAFEMS_R0031_3()
    println("""
    NAFEMS publication R0031/3 Composite plate test.
    Simply supported on all four edges.  Uniform transverse  loading.
    The modeled part is one quarter of the full plate here.
    """)

    # This is a test recommended by the National Agency for Finite Element Methods
    # and Standards (U.K.): Test R0031/3 from NAFEMS publication R0031, “Composites
    # Benchmarks,” February 1995.
    t0 = time()
    # Skin (face) material parameters
    E1s = 1.0e7*phun("psi")
    E2s = 0.4e7*phun("psi")
    E3s = 0.4e7*phun("psi")
    nu12s = 0.3
    nu13s = 0.3
    nu23s = 0.3
    G12s = 0.1875e7*phun("psi")
    G13s = 0.1875e7*phun("psi")
    G23s = 0.1875e7*phun("psi")
    # Core material parameters
    E1c = 10.0*phun("psi")
    E2c = 10.0*phun("psi")
    E3c = 10e4.*phun("psi")
    nu12c = 0.
    nu13c = 0.
    nu23c = 0.
    G12c = 10.0*phun("psi")
    G13c = 3.0e4*phun("psi")
    G23c = 1.2e4*phun("psi")
    L = 10.0*phun("in") # side of the square plate
    nL = 8 # number of elements along the side of the plate
    tolerance = 0.0001*phun("in")
    xs = collect(linearspace(0.0, L/2, nL+1))
    ys = collect(linearspace(0.0, L/2, nL+1))
    ts = [0.028; 0.75; 0.028]*phun("in")
    nts = [2; 3;  2; ] # number of elements through the thickness
    tmag = 100*phun("psi")

    # Generate mesh
    fens,fes = H8layeredplatex(xs, ys, ts, nts)
    fens,fes = H8toH20(fens,fes)

    MR = DeforModelRed3D
    skinmaterial = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    0.0, 0.0, 0.0)
    corematerial = MatDeforElastOrtho(MR,
    0.0, E1c, E2c, E3c,
    nu12c, nu13c, nu23c,
    G12c, G13c, G23c,
    0.0, 0.0, 0.0)

    gr = GaussRule(3, 3)

    rl1 = selectelem(fens, fes, label=1)
    skinbot = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rl1), gr), skinmaterial))

    rl3 = selectelem(fens, fes, label=3)
    skintop = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rl3), gr), skinmaterial))

    rl2 = selectelem(fens, fes, label=2)
    core = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rl2), gr), corematerial))

    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    lxL2 = selectnode(fens, box=[L/2 L/2 -Inf Inf -Inf Inf], inflate=tolerance)
    ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    lyL2 = selectnode(fens, box=[-Inf Inf L/2 L/2 -Inf Inf], inflate=tolerance)

    ex0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
    exL2 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lxL2 )
    ey0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>ly0 )
    eyL2 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lyL2 )

    bfes = meshboundary(fes)
    ttopl = selectelem(fens, bfes; facing=true, direction = [0.0 0.0 1.0])
    Trac = FDataDict("traction_vector"=>[0.0; 0.0; -tmag], "femm"=>FEMMBase(IntegDomain(subset(bfes, ttopl), GaussRule(2, 3))))

    modeldata = FDataDict("fens"=>fens, "regions"=>[skinbot, core, skintop], "essential_bcs"=>[ex0, exL2, ey0, eyL2], "traction_bcs"=> [Trac])
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    u = modeldata["u"]
    geom = modeldata["geom"]
    lcenter = selectnode(fens, box=[L/2 L/2  L/2 L/2 -Inf Inf], inflate=tolerance)
    cdis = mean(u.values[lcenter, 3])/phun("in")
    println("Center node displacements $(cdis) [in]; NAFEMS-R0031-3 lists –0.123	[in]")
    println("")

    File =  "NAFEMS-R0031-3-plate.vtk"
    vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H20; scalars = [("Layer", fes.label)], vectors = [("displacement", u.values)])
    @async run(`"paraview.exe" $File`)

    println("Done")
    true

end # NAFEMS_R0031_3

function allrun()
    println("#####################################################")
    println("# NAFEMS_R0031_1 ")
    NAFEMS_R0031_1()
    println("#####################################################")
    println("# NAFEMS_R0031_1_H20 ")
    NAFEMS_R0031_1_H20()
    println("#####################################################")
    println("# NAFEMS_R0031_2_pressure ")
    NAFEMS_R0031_2_pressure()
    println("#####################################################")
    println("# NAFEMS_R0031_2_both ")
    NAFEMS_R0031_2_both()
    println("#####################################################")
    println("# NAFEMS_R0031_3 ")
    NAFEMS_R0031_3()
    return true
end # function allrun

end # module NAFEMS_R0031_examples


text:
# Copyright 2019 Tobias Frilling
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

@doc """
    mutable struct Settings

A struct that encapsulates commonly used settings for Umatrix.

If no explicit settings are given, most functions implicitly use
`Umatrix.defaultSettings`. There are a couple of ways to change settings:

```julia
# change default settings for all function calls
Umatrix.defaultSettings.some_setting = 20
# create new settings with the kw constructor
s1 = Umatrix.Settings(setting1 = 23, setting2 = 42)
some_function(..., settings = s1)
# copy existing settings and change it
s2 = copy(s1)
s2.setting3 = 99
some_function(..., settings = s2)
```

# Fields and defaults
```julia
latticeSize           :: Tuple{Int64,Int64}     = (50, 80)
toroid                :: Bool                   = true
epochs                :: Int64                  = 24
learningRate          :: Tuple{Float64,Float64} = (0.5, 0.1)
learningRateCooling   :: Symbol                 = :linear
radius                :: Tuple{Int64,Int64}     = (24, 1)
radiusCooling         :: Symbol                 = :linear
initMethod            :: Symbol                 = :uniform_min_max
neighbourhoodKernel   :: Symbol                 = :gauss
distance              :: Distances.PreMetric    = Distances.Euclidean()
shiftToHighestDensity :: Bool                   = true
```
""" Settings
Base.@kwdef mutable struct Settings
    latticeSize::Tuple{Int, Int}         = (50,80)
    toroid::Bool                         = true
    epochs::Int                          = 24
    learningRate::Tuple{Float64,Float64} = (0.5,0.1)
    learningRateCooling::Symbol          = :linear
    radius::Tuple{Int,Int}               = (24, 1)
    radiusCooling::Symbol                = :linear
    initMethod::Symbol                   = :uniform_min_max
    neighbourhoodKernel::Symbol          = :gauss
    distance::PreMetric                  = Euclidean()
    shiftToHighestDensity::Bool          = true
end

const defaultSettings = Settings()

# macro to unpack a settings value into local scope
macro unpack(exp)
    keys = isa(exp, Symbol) ? [exp] : exp.args
    assigments = [:( $key = getproperty(settings, $(Expr(:quote, key))) ) for key in keys]
    esc(Expr(:block, assigments...))
end

Base.copy(s::Settings) = deepcopy(s)

function Base.show(io::IO, ::MIME"text/plain", settings::Settings)
    print(io, "Umatrix settings:")
    for name in fieldnames(Settings)
        print(io, "\n$(name) = ", getproperty(settings, name))
    end
end

ext: jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Clang_assert_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Clang_assert")
JLLWrappers.@generate_main_file("Clang_assert", UUID("72d370af-b0e4-5787-9b36-2e8c6f54fb93"))
end  # module Clang_assert_jll


text:
mutable struct Coloring
    hue::Luxor.RGB
    alpha::Float64
end


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()
##

using Plots
using MCPhylo
using CSV
using DataFrames
using LinearAlgebra

pyplot()
##

function getEqulibirium(rates)
      q = zeros(4, 4)
      indices = [
            (2, 1),
            (3, 1),
            (4, 1),
            (1, 2),
            (3, 2),
            (4, 2),
            (1, 3),
            (2, 3),
            (4, 3),
            (1, 4),
            (2, 4),
            (3, 4),
      ]
      q[[CartesianIndex(x...) for x in indices]] = rates
      dia = sum(q, dims = 2)
      q[diagind(4, 4)] = -dia
      eq = svd(q).U[:, end]
      eq /= sum(eq)
end

##

fpairs = CSV.read("../data/fpairs.txt", DataFrame, header = false)[:, 1]

nrc = findall(occursin.("NRc", fpairs))
##

theme(:solarized_light)
#
upscale = 2
fntsm = Plots.font("sans-serif", pointsize = round(10.0 * upscale))
fntlg = Plots.font("sans-serif", pointsize = round(15.0 * upscale))
default(titlefont = fntlg, guidefont = fntlg, tickfont = fntsm, legendfont = fntsm)
default(size = (1000 * upscale, 1000 * upscale)) #Plot canvas size
#sl_palette = PlotThemes._themes[:solarized_light].defaults[:palette]

##

sh = 0.05 # shift of arrows
rml = 1 # rate multiplier
brd = 0.3 # border around circles
dm = 30 * upscale # diameter

arPars = [2,1]

x = [0, 0, 1, 1]
y = [0, 1, 0, 1]

##
ctmc_ = []

for charNum = 1:28
      charName = fpairs[charNum]
      sim = read(
            "modelFitting/output/universal_$(lpad(charNum, 2, "0")).jls",
            ModelChains,
      )
      samples = vcat([sim[:, :s_rates, :].value[:, :, k] for k in sim.chains]...)
      if charNum ∈ nrc
            samples = samples[:, [4, 6, 5, 1, 3, 2, 11, 10, 12, 8, 7, 9]]
      end
      rates = exp.(mapslices(median, samples, dims = 1))[1, :]
      eqs = getEqulibirium(rates)
      features = string.(split(charName, "-"))
      xNames = [features[1][1:1], features[1][2:end]]
      yNames = [features[2][1:1], features[2][2:end]]
      p = plot(
            title = charName,
            legend = false,
            xlim = [-.2, 1.2],
            ylim = [-.2, 1.2],
            aspectratio = :equal,
            xticks = ([0, 1], [xNames[2] * xNames[1], xNames[1] * xNames[2]]),
            yticks = ([0, 1], [yNames[2] * yNames[1], yNames[1] * yNames[2]]),
            size=(4*500, 500*7),
      )
      for i = 1:4
            scatter!(x[i:i], y[i:i], markersize = dm * sqrt.(eqs[i:i]))
      end
      s, d = 2, 1
      plot!(
            [x[s], x[d]] .- sh,
            [y[s] - brd, y[d] + brd],
            color = :black,
            arrow = arrow(arPars...),
            linewidth = rml * rates[1],
      )
      s, d = 3, 1
      plot!(
            [x[s] - brd, x[d] + brd],
            [y[s], y[d]] .- sh,
            color = :black,
            arrow = arrow(arPars...),
            linewidth = rml * rates[2],
      )
      s, d = 4, 1
      plot!(
            [x[s] - brd / sqrt(2), x[d] + brd / sqrt(2)] .- sh / sqrt(2),
            [y[s] - brd / sqrt(2), y[d] + brd / sqrt(2)] .- -sh / sqrt(2),
            color = :black,
            arrow = arrow(arPars...),
            linewidth = rml * rates[3],
      )
      s, d = 1, 2
      plot!(
            [x[s], x[d]] .+ sh,
            [y[s] + brd, y[d] - brd],
            color = :black,
            arrow = arrow(arPars...),
            linewidth = rml * rates[4],
      )
      s, d = 3, 2
      plot!(
            [x[s] - brd / sqrt(2), x[d] + brd / sqrt(2)] .- sh / sqrt(2),
            [y[s] + brd / sqrt(2), y[d] - brd / sqrt(2)] .+ -sh / sqrt(2),
            color = :black,
            arrow = arrow(arPars...),
            linewidth = rml * rates[5],
      )
      s, d = 4, 2
      plot!(
            [x[s] - brd, x[d] + brd],
            [y[s], y[d]] .- sh,
            color = :black,
            arrow = arrow(arPars...),
            linewidth = rml * rates[6],
      )
      s, d = 1, 3
      plot!(
            [x[s] + brd, x[d] - brd],
            [y[s], y[d]] .+ sh,
            color = :black,
            arrow = arrow(arPars...),
            linewidth = rml * rates[7],
      )
      s, d = 2, 3
      plot!(
            [x[s] + brd / sqrt(2), x[d] - brd / sqrt(2)] .+ sh / sqrt(2),
            [y[s] - brd / sqrt(2), y[d] + brd / sqrt(2)] .- -sh / sqrt(2),
            color = :black,
            arrow = arrow(arPars...),
            linewidth = rml * rates[8],
      )
      s, d = 4, 3
      plot!(
            [x[s], x[d]] .- sh,
            [y[s] - brd, y[d] + brd],
            color = :black,
            arrow = arrow(arPars...),
            linewidth = rml * rates[9],
      )
      s, d = 1, 4
      plot!(
            [x[s] + brd / sqrt(2), x[d] - brd / sqrt(2)] .+ sh / sqrt(2),
            [y[s] + brd / sqrt(2), y[d] - brd / sqrt(2)] .+ -sh / sqrt(2),
            color = :black,
            arrow = arrow(arPars...),
            linewidth = rml * rates[10],
      )
      s, d = 2, 4
      plot!(
            [x[s] + brd, x[d] - brd],
            [y[s], y[d]] .+ sh,
            color = :black,
            arrow = arrow(arPars...),
            linewidth = rml * rates[11],
      )
      s, d = 3, 4
      plot!(
            [x[s], x[d]] .+ sh,
            [y[s] + brd, y[d] - brd],
            color = :black,
            arrow = arrow(arPars...),
            linewidth = rml * rates[12],
      )
      push!(ctmc_,p)
end
##
plot(ctmc_...,layout=grid(7,4, heights=repeat([1/7],28)))

savefig("../img/ctmc.pdf")

max_stars_repo_name: gerhardJaeger/phylogeneticTypology

max_issues_repo_name: gerhardJaeger/phylogeneticTypology

max_forks_repo_name: gerhardJaeger/phylogeneticTypology


text:
using JSON
using Profile
using Random
using Distributed
using StatProfilerHTML
include("src/Flips.jl")
using .Flips

function main()
    graph_data = JSON.parsefile("horizontal.json")
    graph = IndexedGraph(graph_data, "population");
    plan = Plan(graph, "district")
    t = MersenneTwister(0);
    n = 8
    unique_plans = 1
    for i in 1:10000
        println("--- ", i, " ---")
        flip, _, reasons = recom_until_step(graph, plan, 355, 365, t)
        update!(plan, graph, flip)
        unique_plans += 1
        println(reasons)
        if unique_plans % 100 == 0
            println(plan.assignment)
        end
    end
end

main()


text:
#TODO: the inverse transform logic need be prepared
struct pipeline
    transforms
end

Pipeline(Fn) = pipeline(Tuple([Fn]))

function Pipeline( FnStack... )
    pipe = Array{Any,1}( undef, length( FnStack ))
    for (i, fn) in enumerate( FnStack )
        pipe[ i ] = fn
    end
    return pipeline(Tuple( pipe ))
end

function (P::pipeline)(X)
	foldl((X, p) -> p(X), P.transforms, init = X)
end

struct Dataset
    data::AbstractVector
    tfms::NTuple
end

#TODO: Code can be better at checking the
#no. of pipelines
function Dataset(data; tfms=nothing)
    if length(tfms) == 2
        Dataset(data, (Pipeline(tfms[1]...), Pipeline(tfms[2]...)))
    else
        Dataset(data, (Pipeline(tfms[1]...)))
    end
end

nobs(dataset::Dataset) = length(dataset.data)

#TODO:Replace indexing with first and last
function getobs(dataset::Dataset, idx::Int64)
    data = dataset.data[idx]

    if dataset.tfms[1]==nothing || dataset.tfms[2]==nothing
       (dataset.tfms[1](data))
    else
       (dataset.tfms[1](data), dataset.tfms[2](data))
    end
end

function (dataset::Dataset)(idx::Int64)
	getobs(dataset, idx)
end


text:
export LiftedHRepresentation, LiftedVRepresentation

# H-Represenation

# No copy since I do not modify anything and a copy is done when building a polyhedron
mutable struct LiftedHRepresentation{N, T, MT<:AbstractMatrix{T}} <: MixedHRep{N, T}
    # Ax >= 0, it is [b -A] * [z; x] where z = 1
    A::MT
    linset::BitSet

    function LiftedHRepresentation{N, T, MT}(A::MT, linset::BitSet=BitSet()) where {N, T, MT}
        if !isempty(linset) && last(linset) > size(A, 1)
            error("The elements of linset should be between 1 and the number of rows of A")
        end
        if size(A, 2) != N+1
            error("dimension does not match")
        end
        new{N, T, MT}(A, linset)
    end
end

similar_type(::Type{LiftedHRepresentation{M, S, MT}}, ::FullDim{N}, ::Type{T}) where {M, S, N, T, MT} = LiftedHRepresentation{N, T, similar_type(MT, T)}
hvectortype(p::Type{LiftedHRepresentation{N, T, MT}}) where {N, T, MT} = vectortype(MT)

LiftedHRepresentation{N, T}(A::AbstractMatrix{T}, linset::BitSet=BitSet()) where {N, T} = LiftedHRepresentation{N, T, typeof(A)}(A, linset)
LiftedHRepresentation{N, T}(A::AbstractMatrix, linset::BitSet=BitSet()) where {N, T} = LiftedHRepresentation{N, T}(AbstractMatrix{T}(A), linset)
LiftedHRepresentation(A::AbstractMatrix{T}, linset::BitSet=BitSet()) where T = LiftedHRepresentation{size(A, 2) - 1, T}(A, linset)

LiftedHRepresentation(h::HRepresentation{N, T}) where {N, T} = LiftedHRepresentation{N, T}(h)
LiftedHRepresentation{N, T}(h::HRepresentation{N}) where {N, T} = LiftedHRepresentation{N, T, hmatrixtype(typeof(h), T)}(h)

function LiftedHRepresentation{N, T, MT}(hyperplanes::ElemIt{<:HyperPlane{N, T}}, halfspaces::ElemIt{<:HalfSpace{N, T}}) where {N, T, MT}
    nhyperplane = length(hyperplanes)
    nhrep = nhyperplane + length(halfspaces)
    A = emptymatrix(MT, nhrep, N+1)
    linset = BitSet(1:nhyperplane)
    for (i, h) in enumerate(hyperplanes)
        A[i,2:end] = -h.a
        A[i,1] = h.β
    end
    for (i, h) in enumerate(halfspaces)
        A[nhyperplane+i,2:end] = -h.a
        A[nhyperplane+i,1] = h.β
    end
    LiftedHRepresentation{N, T}(A, linset)
end

Base.copy(ine::LiftedHRepresentation{N,T}) where {N,T} = LiftedHRepresentation{N,T}(copy(ine.A), copy(ine.linset))

Base.isvalid(hrep::LiftedHRepresentation{N, T}, idx::HIndex{N, T}) where {N, T} = 0 < idx.value <= size(hrep.A, 1) && (idx.value in hrep.linset) == islin(idx)
Base.done(idxs::HIndices{N, T, <:LiftedHRepresentation{N, T}}, idx::HIndex{N, T}) where {N, T} = idx.value > size(idxs.rep.A, 1)
Base.get(hrep::LiftedHRepresentation{N, T}, idx::HIndex{N, T}) where {N, T} = valuetype(idx)(-hrep.A[idx.value,2:end], hrep.A[idx.value,1])

# V-Representation

mutable struct LiftedVRepresentation{N, T, MT<:AbstractMatrix{T}} <: MixedVRep{N, T}
    R::MT # each row is a vertex if the first element is 1 and a ray otherwise
    linset::BitSet

    function LiftedVRepresentation{N, T, MT}(R::MT, linset::BitSet=BitSet([])) where {N, T, MT}
        if length(R) > 0 && size(R, 2) != N+1
            error("dimension does not match")
        end
        if !isempty(linset) && last(linset) > size(R, 1)
            error("The elements of linset should be between 1 and the number of rows of R")
        end
        new{N, T, MT}(R, linset)
    end
end

similar_type(::Type{LiftedVRepresentation{M, S, MT}}, ::FullDim{N}, ::Type{T}) where {M, S, N, T, MT} = LiftedVRepresentation{N, T, similar_type(MT, T)}
vvectortype(p::Type{LiftedVRepresentation{N, T, MT}}) where {N, T, MT} = vectortype(MT)

LiftedVRepresentation{N, T}(R::AbstractMatrix{T}, linset::BitSet=BitSet()) where {N, T} = LiftedVRepresentation{N, T, typeof(R)}(R, linset)
LiftedVRepresentation{N, T}(R::AbstractMatrix, linset::BitSet=BitSet()) where {N, T} = LiftedVRepresentation{N, T}(AbstractMatrix{T}(R), linset)
LiftedVRepresentation(R::AbstractMatrix{T}, linset::BitSet=BitSet()) where T = LiftedVRepresentation{size(R, 2) - 1, T}(R, linset)

LiftedVRepresentation(v::VRepresentation{N, T}) where {N, T} = LiftedVRepresentation{N, T}(v)
LiftedVRepresentation{N, T}(v::VRepresentation{N}) where {N, T} = LiftedVRepresentation{N, T, vmatrixtype(typeof(v), T)}(v)

function LiftedVRepresentation{N, T, MT}(vits::VIt{N, T}...) where {N, T, MT}
    points, lines, rays = fillvits(FullDim{N}(), vits...)
    npoint = length(points)
    nline = length(lines)
    nray = length(rays)
    nvrep = npoint + nline + nray
    R = emptymatrix(MT, nvrep, N+1)
    linset = BitSet()
    function _fill(offset, z, ps)
        for (i, p) in enumerate(ps)
            R[offset + i,2:end] = coord(p)
            R[offset + i,1] = z
            if islin(p)
                push!(linset, offset + i)
            end
        end
    end
    _fill(0, one(T), points)
    _fill(npoint, zero(T), lines)
    _fill(npoint+nline, zero(T), rays)
    LiftedVRepresentation{N, T}(R, linset)
end

Base.copy(ext::LiftedVRepresentation{N,T}) where {N,T} = LiftedVRepresentation{N,T}(copy(ext.R), copy(ext.linset))

nvreps(ext::LiftedVRepresentation) = size(ext.R, 1)

function isrowpoint(ext::LiftedVRepresentation{N,T}, i) where {N,T}
    ispoint = ext.R[i,1]
    @assert ispoint == zero(T) || ispoint == one(T)
    ispoint == one(T)
end

function Base.isvalid(vrep::LiftedVRepresentation{N, T}, idx::VIndex{N, T}) where {N, T}
    isp = isrowpoint(vrep, idx.value)
    isl = (idx.value in vrep.linset)
    @assert !isp || !isl # if isp && isl, it is a symmetric point but it is not allowed to mix symmetric points and points
    0 < idx.value <= nvreps(vrep) && isp == ispoint(idx) && isl == islin(idx)
end
Base.done(idxs::VIndices{N, T, <:LiftedVRepresentation{N, T}}, idx::VIndex{N, T}) where {N, T} = idx.value > size(idxs.rep.R, 1)
Base.get(vrep::LiftedVRepresentation{N, T}, idx::VIndex{N, T}) where {N, T} = valuetype(idx)(vrep.R[idx.value,2:end])


text:
using RecordedArrays: TimeSeries, PhasePortrait

@testset "Scalar" begin
    c = DiscreteClock(10)
    s = recorded(DynamicEntry, c, fill(0))
    @test size(s) == ()
    for t in c
        s[] = t
        @test s[] == t
    end
    e1 = getentries(s)
    @test e1 == RecordedArrays.getrecord(s)[]
    @test getts(e1) == 0:10
    @test getvs(e1) == 0:10
    @test ArrayInterface.parent_type(getrecord(s)) == DynamicEntry{Int,Int}
end

@testset "Vector" begin
    @testset "Static" begin
        c = ContinuousClock(10)
        v = recorded(StaticEntry, c, [1])
        increase!(c, 1) # t = 1
        resize!(v, 2)[2] = 2
        increase!(c, 1) # t = 2
        @test push!(v, 3) == 1:3
        increase!(c, 1) # t = 3
        @test resize!(v, (Bool[0, 0, 1],)) == [3]
        increase!(c, 1) # t = 4
        resize!(v, (3,))
        v[2:3] = 4:5 # this a setindex!
        increase!(c, 1) # t = 5
        @test deleteat!(v, 2) == [3, 5]
        increase!(c, 1) # t = 6
        @test insert!(v, 2, 6) == [3, 6, 5]
        increase!(c, 1) # t = 7
        @test append!(v, [7]) == [3, 6, 5, 7]
        increase!(c, 1) # t = 8
        @test append!(v, 8) == [3, 6, 5, 7, 8]
        increase!(c, 1) # t = 9
        @test resize!(v, 1, Bool[1, 0, 0, 0, 1]) == [3, 8]
        increase!(c, 1) # t = 10
        @test resize!(v, 1) == [3]
        es = getentries(v)
        @test getts(es[1]) == [0, 3]
        @test getts(es[2]) == [1, 3]
        @test getts(es[3]) == [2]
        @test getts(es[4]) == [4, 5]
        @test getts(es[5]) == [4, 9]
        @test getts(es[6]) == [6, 9]
        @test getts(es[7]) == [7, 9]
        @test getts(es[8]) == [8, 10]
        @test gettime(es, 0) == [1, 0, 0, 0, 0, 0, 0, 0]
        @test gettime(es, 1) == [1, 2, 0, 0, 0, 0, 0, 0]
        @test gettime(es, 2) == [1, 2, 3, 0, 0, 0, 0, 0]
        @test gettime(es, 3) == [1, 2, 3, 0, 0, 0, 0, 0]
        @test gettime(es, 4) == [0, 0, 3, 4, 5, 0, 0, 0]
        @test gettime(es, 5) == [0, 0, 3, 4, 5, 0, 0, 0]
        @test gettime(es, 6) == [0, 0, 3, 0, 5, 6, 0, 0]
        @test gettime(es, 7) == [0, 0, 3, 0, 5, 6, 7, 0]
        @test gettime(es, 8) == [0, 0, 3, 0, 5, 6, 7, 8]
        @test gettime(es, 9) == [0, 0, 3, 0, 5, 6, 7, 8]
        @test gettime(es, 10) == [0, 0, 3, 0, 0, 0, 0, 8]
        @test gettime(es, 11) == [0, 0, 3, 0, 0, 0, 0, 0]
        @test gettime(es, -0.5) == [0, 0, 0, 0, 0, 0, 0, 0]
        @test gettime(es, 1.5) == [1, 2, 0, 0, 0, 0, 0, 0]
        @test gettime(es, 3.5) == [0, 0, 3, 0, 0, 0, 0, 0]
        @test gettime(es, 9.5) == [0, 0, 3, 0, 0, 0, 0, 8]
        ts = -0.5:0.5:10.5
        m = gettime(es, ts)
        for (i, t) in enumerate(ts)
            @test m[i, :] ==
                  gettime(es, t) ==
                  map(&, gettime(es, round(t, RoundUp)), gettime(es, round(t, RoundDown)))
        end
    end
    @testset "Dynamic" begin
        c = ContinuousClock(10)
        v = recorded(DynamicEntry, c, [1])
        increase!(c, 1) # t = 1
        resize!(v, 2)[:] = 2:3
        @test v == 2:3
        increase!(c, 1) # t = 2
        v[2] = 4
        @test v == [2, 4]
        increase!(c, 1) # t = 3
        @test resize!(v, (1,)) == [2]
        increase!(c, 1) # t = 4
        v[1] = 5
        increase!(c, 1) # t = 5
        resize!(v, (3,))[2:3] = 6:7
        @test v == 5:7
        es = getentries(v)
        @test getts(es[1]) == [0, 1, 4]
        @test getts(es[2]) == [1, 2]
        @test getts(es[3]) == [5]
        @test getts(es[4]) == [5]
        @test getvs(es[1]) == [1, 2, 5]
        @test getvs(es[2]) == [3, 4]
        @test getvs(es[3]) == [6]
        @test getvs(es[4]) == [7]

        @test gettime(es, 0) == [1, 0, 0, 0]
        @test gettime(es, 1) == [2, 3, 0, 0]
        @test gettime(es, 2) == [2, 4, 0, 0]
        @test gettime(es, 3) == [2, 4, 0, 0]
        @test gettime(es, 4) == [5, 4, 0, 0]
        @test gettime(es, 5) == [5, 4, 6, 7]
        @test gettime(es, -0.5) == [0, 0, 0, 0]
        @test gettime(es, 1.5) == [2, 3, 0, 0]
        @test gettime(es, 2.5) == [2, 4, 0, 0]
        @test gettime(es, 3.5) == [2, 4, 0, 0]
        @test gettime(es, 0:2:4) == [1 0 0 0; 2 4 0 0; 5 4 0 0]
        @test gettime(es, -1:2:5) == [0 0 0 0; 2 3 0 0; 2 4 0 0; 5 4 6 7]
        ts = 0:0.5:5.5
        a = gettime(es, ts)
        for (i, t) in enumerate(ts)
            @test a[i, :] == gettime(es, t) == gettime(es, round(t, RoundDown))
        end
    end
end

@testset "Matrix" begin
    @testset "Static" begin
        c = ContinuousClock(10)
        m = recorded(StaticEntry, c, fill(1, 1, 1))
        increase!(c, 1) # t = 1
        resize!(m, (2, 2))[2:4] = 2:4
        increase!(c, 1) # t = 2
        resize!(m, 2, 3)[5:6] = 5:6
        @test vec(m) == 1:6
        increase!(c, 1) # t = 3
        @test resize!(m, (:, 2:3)) == reshape(3:6, 2, 2)
        increase!(c, 1) # t = 4
        @test resize!(m, 1, 1) == [3 5]
        increase!(c, 1) # t = 5
        resize!(m, 1, 3)[not(1), :] = reshape(7:10, 2, 2)
        increase!(c, 1) # t = 6
        resize!(m, 2, 3)[:, 3] = reshape(11:13, 1, 3)
        @test m == [
            3 5 11
            7 9 12
            8 10 13
        ]
        increase!(c, 1) # t = 7
        @test resize!(m, (not(2), not(2))) |> vec == [3, 8, 11, 13]
        es = getentries(m)
        ks = keys(es)
        vs = values(es)
        for (k, v) in zip(ks, vs)
            @test es[k] == v
        end
        @test getts(es[1, 1]) == [0, 3] # 1
        @test getts(es[2, 1]) == [1, 3] # 2
        @test getts(es[1, 2]) == [1]    # 3
        @test getts(es[2, 2]) == [1, 4] # 4
        @test getts(es[1, 3]) == [2, 7] # 5
        @test getts(es[2, 3]) == [2, 4] # 6
        @test getts(es[3, 2]) == [5, 7] # 7
        @test getts(es[4, 2]) == [5]    # 8
        @test getts(es[3, 3]) == [5, 7] # 9
        @test getts(es[4, 3]) == [5, 7] # 10
        @test getts(es[1, 4]) == [6]    # 11
        @test getts(es[3, 4]) == [6, 7] # 12
        @test getts(es[4, 4]) == [6]    # 13
        @test gettime(es, 0) == [1 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 0]
        @test gettime(es, 1) == [1 3 0 0; 2 4 0 0; 0 0 0 0; 0 0 0 0]
        @test gettime(es, 2) == [1 3 5 0; 2 4 6 0; 0 0 0 0; 0 0 0 0]
        @test gettime(es, 3) == [1 3 5 0; 2 4 6 0; 0 0 0 0; 0 0 0 0]
        @test gettime(es, 4) == [0 3 5 0; 0 4 6 0; 0 0 0 0; 0 0 0 0]
        @test gettime(es, 5) == [0 3 5 0; 0 0 0 0; 0 7 9 0; 0 8 10 0]
        @test gettime(es, 6) == [0 3 5 11; 0 0 0 0; 0 7 9 12; 0 8 10 13]
        @test gettime(es, 7) == [0 3 5 11; 0 0 0 0; 0 7 9 12; 0 8 10 13]
        @test gettime(es, 8) == [0 3 0 11; 0 0 0 0; 0 0 0 0; 0 8 0 13]
        @test gettime(es, -0.5) == [0 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 0]
        @test gettime(es, 0.5) == gettime(es, 0)
        @test gettime(es, 1.5) == gettime(es, 1)
        @test gettime(es, 3.5) == gettime(es, 4)
        @test gettime(es, 4.5) == [0 3 5 0; 0 0 0 0; 0 0 0 0; 0 0 0 0]
        @test gettime(es, 6.5) == gettime(es, 6)
        ts = 0:0.5:7.5
        a = gettime(es, ts)
        for (i, t) in enumerate(ts)
            @test a[i, :, :] ==
                  gettime(es, t) ==
                  map(&, gettime(es, round(t, RoundUp)), gettime(es, round(t, RoundDown)))
        end
    end
    @testset "DynamicEntry" begin
        c = ContinuousClock(10)
        m = recorded(DynamicEntry, c, fill(1, 1, 1))
        increase!(c, 1) # t = 1
        resize!(m, (2, 2))[2:4] = 2:4
        increase!(c, 1) # t = 2
        m[1] = 5
        increase!(c, 1) # t = 3
        resize!(m, 2, 3)[5:6] = 6:7
        increase!(c, 1) # t = 4
        m[2, 3] = 8
        increase!(c, 1) # t = 5
        resize!(m, (:, not(2))) |> vec == [5, 2, 6, 8]
        increase!(c, 1) # t = 6
        m[2] = 9
        increase!(c, 1) # t = 7
        m[3] = 10
        @test m[:] == [5, 9, 10, 8]
        increase!(c, 1) # t = 8
        resize!(m, (3, :))[3, :] = 11:12
        increase!(c, 1) # t = 9
        m[3, 1] = 13
        @test m == [5 10; 9 8; 13 12]
        es = getentries(m)
        @test getts(es[1, 1]) == [0, 2]
        @test getts(es[2, 1]) == [1, 6]
        @test getts(es[1, 2]) == [1]
        @test getts(es[2, 2]) == [1]
        @test getts(es[1, 3]) == [3, 7]
        @test getts(es[2, 3]) == [3, 4]
        @test getts(es[3, 1]) == [8, 9]
        @test getts(es[3, 3]) == [8]

        @test getvs(es[1, 1]) == [1, 5]
        @test getvs(es[2, 1]) == [2, 9]
        @test getvs(es[1, 2]) == [3]
        @test getvs(es[2, 2]) == [4]
        @test getvs(es[1, 3]) == [6, 10]
        @test getvs(es[2, 3]) == [7, 8]
        @test getvs(es[3, 1]) == [11, 13]
        @test getvs(es[3, 3]) == [12]

        @test gettime(es, 0) == [1 0 0; 0 0 0; 0 0 0]
        @test gettime(es, 1) == [1 3 0; 2 4 0; 0 0 0]
        @test gettime(es, 2) == [5 3 0; 2 4 0; 0 0 0]
        @test gettime(es, 3) == [5 3 6; 2 4 7; 0 0 0]
        @test gettime(es, 4) == [5 3 6; 2 4 8; 0 0 0]
        @test gettime(es, 5) == [5 3 6; 2 4 8; 0 0 0]
        @test gettime(es, 6) == [5 3 6; 9 4 8; 0 0 0]
        @test gettime(es, 7) == [5 3 10; 9 4 8; 0 0 0]
        @test gettime(es, 8) == [5 3 10; 9 4 8; 11 0 12]
        @test gettime(es, 9) == [5 3 10; 9 4 8; 13 0 12]
        ts = 0:0.5:9.5
        a = gettime(es, ts)
        for (i, t) in enumerate(ts)
            @test a[i, :, :] == gettime(es, t) == gettime(es, round(t, RoundDown))
        end
    end
end

@testset "MISC" begin
    @testset "Recipes" begin
        v = recorded(DynamicEntry, ContinuousClock(10), [1, 2])
        test_recipe(TimeSeries((v,)), [([0], [1]), ([0], [2])])
        test_recipe(TimeSeries((getentries(v),)), [([0], [1]), ([0], [2])])
        test_recipe(TimeSeries((getentries(v)...,)), [([0], [1]), ([0], [2])])
        test_recipe(PhasePortrait((v,)), [([1], [2])])
        test_recipe(PhasePortrait((getentries(v),)), [([1], [2])])
        test_recipe(PhasePortrait((getentries(v)...,)), [([1], [2])])
        m = recorded(DynamicEntry, ContinuousClock(10), [1 2])
        # the order may be different, thus only test TimeSeries
        test_recipe(TimeSeries((m,)), [([0], [1]), ([0], [2])])
        test_recipe(TimeSeries((getentries(m),)), [([0], [1]), ([0], [2])])
        test_recipe(TimeSeries((getentries(m)...,)), [([0], [1]), ([0], [2])])
    end
    @testset "Show" begin
        v = recorded(DynamicEntry, ContinuousClock(10), [1, 2])
        test_show(v, "2-element recorded(::$(Vector{Int})):\n 1\n 2")
        r = getrecord(v)
        test_show(r, "2-element $(typeof(r))")
        e = getentries(v)[1]
        test_show(e, "$(typeof(e)) with timestamps:\n 0")
    end
    @testset "sizehint!" begin
        # run one time for compile
        v = recorded(DynamicEntry, ContinuousClock(10), [1, 2])
        sizehint!(v, 3)
        sizehint!(getrecord(v), 3)
        @test begin
            v = recorded(DynamicEntry, ContinuousClock(10), [1, 2])
            sizehint!(v, 3)
            @allocated resize!(v, 3)
        end < begin
            v = recorded(DynamicEntry, ContinuousClock(10), [1, 2])
            @allocated resize!(v, 3)
        end
        @test begin
            v = recorded(DynamicEntry, ContinuousClock(10), [1, 2])
            sizehint!(getrecord(v), 3)
            @allocated resize!(v, 3)
        end < begin
            v = recorded(DynamicEntry, ContinuousClock(10), [1, 2])
            @allocated resize!(v, 3)
        end
    end
end


text:
"""
SmetsWouters{T} <: AbstractRepModel{T}

The `SmetsWouters` type defines the structure of the New York Fed DSGE
model.

### Fields

#### Parameters and Steady-States
* `parameters::Vector{AbstractParameter}`: Vector of all time-invariant model
  parameters.

* `steady_state::Vector{AbstractParameter}`: Model steady-state values, computed
  as a function of elements of `parameters`.

* `keys::OrderedDict{Symbol,Int}`: Maps human-readable names for all model
  parameters and steady-states to their indices in `parameters` and
  `steady_state`.

#### Inputs to Measurement and Equilibrium Condition Equations

The following fields are dictionaries that map human-readable names to row and
column indices in the matrix representations of of the measurement equation and
equilibrium conditions.

* `endogenous_states::OrderedDict{Symbol,Int}`: Maps each state to a column in
  the measurement and equilibrium condition matrices.

* `exogenous_shocks::OrderedDict{Symbol,Int}`: Maps each shock to a column in
  the measurement and equilibrium condition matrices.

* `expected_shocks::OrderedDict{Symbol,Int}`: Maps each expected shock to a
  column in the measurement and equilibrium condition matrices.

* `equilibrium_conditions::OrderedDict{Symbol,Int}`: Maps each equlibrium
  condition to a row in the model's equilibrium condition matrices.

* `endogenous_states_augmented::OrderedDict{Symbol,Int}`: Maps lagged states to
  their columns in the measurement and equilibrium condition equations. These
  are added after `gensys` solves the model.

* `observables::OrderedDict{Symbol,Int}`: Maps each observable to a row in the
  model's measurement equation matrices.

* `pseudo_observables::OrderedDict{Symbol,Int}`: Maps each pseudo-observable to
  a row in the model's pseudo-measurement equation matrices.

#### Model Specifications and Settings

* `spec::String`: The model specification identifier, \"smets_wouters\",
cached here for filepath computation.

* `subspec::String`: The model subspecification number,
indicating that some parameters from the original model spec (\"ss0\")
are initialized differently. Cached here for filepath computation.

* `settings::Dict{Symbol,Setting}`: Settings/flags that affect
computation without changing the economic or mathematical setup of
the model.

* `test_settings::Dict{Symbol,Setting}`: Settings/flags for testing mode

#### Other Fields

* `rng::MersenneTwister`: Random number generator. Can be is seeded to ensure
  reproducibility in algorithms that involve randomness (such as
  Metropolis-Hastings).

* `testing::Bool`: Indicates whether the model is in testing mode. If `true`,
  settings from `m.test_settings` are used in place of those in `m.settings`.

* `observable_mappings::OrderedDict{Symbol,Observable}`: A dictionary that
  stores data sources, series mnemonics, and transformations to/from model
  units. DSGE.jl will fetch data from the Federal Reserve Bank of St. Louis's
  FRED database; all other data must be downloaded by the user. See `load_data`
  and `Observable` for further details.

* `pseudo_observable_mappings::OrderedDict{Symbol,PseudoObservable}`: A
  dictionary that stores names and transformations to/from model units. See
  `PseudoObservable` for further details.
"""
mutable struct SmetsWouters{T} <: AbstractRepModel{T}
    parameters::ParameterVector{T}                         # vector of all time-invariant model parameters
    steady_state::ParameterVector{T}                       # model steady-state values
    keys::OrderedDict{Symbol,Int}                          # human-readable names for all the model
                                                           # parameters and steady-states

    endogenous_states::OrderedDict{Symbol,Int}             # these fields used to create matrices in the
    exogenous_shocks::OrderedDict{Symbol,Int}              # measurement and equilibrium condition equations.
    expected_shocks::OrderedDict{Symbol,Int}               #
    equilibrium_conditions::OrderedDict{Symbol,Int}        #
    endogenous_states_augmented::OrderedDict{Symbol,Int}   #
    observables::OrderedDict{Symbol,Int}                   #
    pseudo_observables::OrderedDict{Symbol,Int}            #

    spec::String                                           # Model specification number (eg "m990")
    subspec::String                                        # Model subspecification (eg "ss0")
    settings::Dict{Symbol,Setting}                         # Settings/flags for computation
    test_settings::Dict{Symbol,Setting}                    # Settings/flags for testing mode
    rng::MersenneTwister                                   # Random number generator
    testing::Bool                                          # Whether we are in testing mode or not

    observable_mappings::OrderedDict{Symbol, Observable}
    pseudo_observable_mappings::OrderedDict{Symbol, PseudoObservable}
end

description(m::SmetsWouters) = "Smets-Wouters Model"

"""
`init_model_indices!(m::SmetsWouters)`

Arguments:
`m:: SmetsWouters`: a model object

Description:
Initializes indices for all of `m`'s states, shocks, and equilibrium conditions.
"""
function init_model_indices!(m::SmetsWouters)
    # Endogenous states
    endogenous_states = [[
        :y_t, :c_t, :i_t, :qk_t, :k_t, :kbar_t, :u_t, :rk_t, :mc_t,
        :π_t, :μ_ω_t, :w_t, :L_t, :R_t, :g_t, :b_t, :μ_t, :z_t,
        :λ_f_t, :λ_f_t1, :λ_w_t, :λ_w_t1, :rm_t, :Ec_t, :Eqk_t, :Ei_t,
        :Eπ_t, :EL_t, :Erk_t, :Ew_t, :y_f_t, :c_f_t,
        :i_f_t, :qk_f_t, :k_f_t, :kbar_f_t, :u_f_t, :rk_f_t, :w_f_t,
        :L_f_t, :r_f_t, :Ec_f_t, :Eqk_f_t, :Ei_f_t, :EL_f_t, :Erk_f_t, :ztil_t];
        [Symbol("rm_tl$i") for i = 1:n_anticipated_shocks(m)]]

    # Exogenous shocks
    exogenous_shocks = [[
        :g_sh, :b_sh, :μ_sh, :z_sh, :λ_f_sh, :λ_w_sh, :rm_sh];
        [Symbol("rm_shl$i") for i = 1:n_anticipated_shocks(m)]]

    # Expectations shocks
    expected_shocks = [
        :Ec_sh, :Eqk_sh, :Ei_sh, :Eπ_sh, :EL_sh, :Erk_sh, :Ew_sh, :Ec_f_sh,
        :Eqk_f_sh, :Ei_f_sh, :EL_f_sh, :Erk_f_sh]

    # Equilibrium conditions
    equilibrium_conditions = [[
        :eq_euler, :eq_inv, :eq_capval, :eq_output, :eq_caputl, :eq_capsrv, :eq_capev,
        :eq_mkupp, :eq_phlps, :eq_caprnt, :eq_msub, :eq_wage, :eq_mp, :eq_res, :eq_g, :eq_b, :eq_μ, :eq_z,
        :eq_λ_f, :eq_λ_w, :eq_rm, :eq_λ_f1, :eq_λ_w1, :eq_Ec,
        :eq_Eqk, :eq_Ei, :eq_Eπ, :eq_EL, :eq_Erk, :eq_Ew, :eq_euler_f, :eq_inv_f,
        :eq_capval_f, :eq_output_f, :eq_caputl_f, :eq_capsrv_f, :eq_capev_f, :eq_mkupp_f, :eq_caprnt_f, :eq_msub_f,
        :eq_res_f, :eq_Ec_f, :eq_Eqk_f, :eq_Ei_f, :eq_EL_f, :eq_Erk_f, :eq_ztil];
        [Symbol("eq_rml$i") for i=1:n_anticipated_shocks(m)]]

    # Additional states added after solving model
    # Lagged states and observables measurement error
    endogenous_states_augmented = [
        :y_t1, :c_t1, :i_t1, :w_t1, :π_t1, :L_t1, :Et_π_t]

    # Observables
    observables = keys(m.observable_mappings)

    # Pseudo-observables
    pseudo_observables = keys(m.pseudo_observable_mappings)

    for (i,k) in enumerate(endogenous_states);           m.endogenous_states[k]           = i end
    for (i,k) in enumerate(exogenous_shocks);            m.exogenous_shocks[k]            = i end
    for (i,k) in enumerate(expected_shocks);             m.expected_shocks[k]             = i end
    for (i,k) in enumerate(equilibrium_conditions);      m.equilibrium_conditions[k]      = i end
    for (i,k) in enumerate(endogenous_states);           m.endogenous_states[k]           = i end
    for (i,k) in enumerate(endogenous_states_augmented); m.endogenous_states_augmented[k] = i+length(endogenous_states) end
    for (i,k) in enumerate(observables);                 m.observables[k]                 = i end
    for (i,k) in enumerate(pseudo_observables);          m.pseudo_observables[k]          = i end
end

function SmetsWouters(subspec::String="ss0";
                      custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}(),
                      testing = false)

    # Model-specific specifications
    spec               = split(basename(@__FILE__),'.')[1]
    subspec            = subspec
    settings           = Dict{Symbol,Setting}()
    test_settings      = Dict{Symbol,Setting}()
    rng                = MersenneTwister(0)

    # initialize empty model
    m = SmetsWouters{Float64}(
            # model parameters and steady state values
            Vector{AbstractParameter{Float64}}(), Vector{Float64}(), OrderedDict{Symbol,Int}(),

            # model indices
            OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(),

            spec,
            subspec,
            settings,
            test_settings,
            rng,
            testing,
            OrderedDict{Symbol,Observable}(),
            OrderedDict{Symbol,PseudoObservable}())

    # Set settings
    model_settings!(m)
    default_test_settings!(m)
    for custom_setting in values(custom_settings)
        m <= custom_setting
    end

    # Set observable and pseudo-observable transformations
    init_observable_mappings!(m)

    # Initialize parameters
    init_parameters!(m)

    init_model_indices!(m)
    init_subspec!(m)
    steadystate!(m)

    return m
end

"""
```
init_parameters!(m::SmetsWouters)
```

Initializes the model's parameters, as well as empty values for the steady-state
parameters (in preparation for `steadystate!(m)` being called to initialize
those).
"""
function init_parameters!(m::SmetsWouters)
    m <= parameter(:α, 0.24, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Normal(0.30, 0.05), fixed=false,
                   description="α: Capital elasticity in the intermediate goods sector's Cobb-Douglas production function.",
                   tex_label="\\alpha")

    m <= parameter(:ζ_p, 0.7813, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1), fixed=false,
                   description="ζ_p: The Calvo parameter. In every period, intermediate goods producers optimize prices with probability (1-ζ_p). With probability ζ_p, prices are adjusted according to a weighted average of the previous period's inflation (π_t1) and steady-state inflation (π_star).",
                   tex_label="\\zeta_p")

    m <= parameter(:ι_p, 0.3291, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.15), fixed=false,
                   description="ι_p: The weight attributed to last period's inflation in price indexation. (1-ι_p) is the weight attributed to steady-state inflation.",
                   tex_label="\\iota_p")

    m <= parameter(:δ, 0.025, fixed=true,
                   description="δ: The capital depreciation rate.",
                   tex_label="\\delta" )

    m <= parameter(:Upsilon, 1.000, (0., 10.), (1e-5, 0.), ModelConstructors.Exponential(), GammaAlt(1., 0.5), fixed=true,
                   description="Υ: The trend evolution of the price of investment goods relative to consumption goods. Set equal to 1.",
                   tex_label="\\Upsilon")

    m <= parameter(:Φ, 1.4672, (1., 10.), (1.00, 10.00), ModelConstructors.Exponential(), Normal(1.25, 0.12), fixed=false,
                   description="Φ: Fixed costs.",
                   tex_label="\\Phi")

    m <= parameter(:S′′, 6.3325, (-15., 15.), (-15., 15.), ModelConstructors.Untransformed(), Normal(4., 1.5), fixed=false,
                   description="S'': The second derivative of households' cost of adjusting investment.",
                   tex_label="S''")

    m <= parameter(:h, 0.7205, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.7, 0.1), fixed=false,
                   description="h: Consumption habit persistence.",
                   tex_label="h")

    m <= parameter(:ppsi, 0.2648, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.15), fixed=false,
                   description="ppsi: Utilization costs.",
                   tex_label="\\psi")

    m <= parameter(:ν_l, 2.8401, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), Normal(2, 0.75), fixed=false,
                   description="ν_l: The coefficient of relative risk aversion on the labor term of households' utility function.", tex_label="\\nu_l")

    m <= parameter(:ζ_w, 0.7937, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1), fixed=false,
                   description="ζ_w: (1-ζ_w) is the probability with which households can freely choose wages in each period. With probability ζ_w, wages increase at a geometrically weighted average of the steady state rate of wage increases and last period's productivity times last period's inflation.",
                   tex_label="\\zeta_w")

    m <= parameter(:ι_w, 0.4425, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.15), fixed=false,
                   description="ι_w: The weight attributed to last period's wage in wage indexation. (1-ι_w) is the weight attributed to steady-state wages.",
                   tex_label="\\iota_w")

    m <= parameter(:λ_w, 1.5000, fixed=true,
                   description="λ_w: The wage markup, which affects the elasticity of substitution between differentiated labor services.",
                   tex_label="\\lambda_w")

    m <= parameter(:β, 0.7420, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.25, 0.1), fixed=false, scaling = x -> 1/(1 + x/100),
                   description="β: Discount rate.",
                   tex_label="\\beta ")

    m <= parameter(:ψ1, 1.7985, (1e-5, 10.), (1e-5, 10.00), ModelConstructors.Exponential(), Normal(1.5, 0.25), fixed=false,
                   description="ψ₁: Weight on inflation gap in monetary policy rule.",
                   tex_label="\\psi_1")

    m <= parameter(:ψ2, 0.0893, (-0.5, 0.5), (-0.5, 0.5), ModelConstructors.Untransformed(), Normal(0.12, 0.05), fixed=false,
                   description="ψ₂: Weight on output gap in monetary policy rule.",
                   tex_label="\\psi_2")

    m <= parameter(:ψ3, 0.2239, (-0.5, 0.5), (-0.5, 0.5), ModelConstructors.Untransformed(), Normal(0.12, 0.05), fixed=false,
                   description="ψ₃: Weight on rate of change of output gap in the monetary policy rule.",
                   tex_label="\\psi_3")

    m <= parameter(:π_star, 0.7000, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.62, 0.1), fixed=false, scaling = x -> 1 + x/100,
                   description="π_star: The steady-state rate of inflation.",
                   tex_label="\\pi_*")

    m <= parameter(:σ_c, 1.2312, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), Normal(1.5, 0.37), fixed=false,
                   tex_label="\\sigma_{c}")

    m <= parameter(:ρ, .8258, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.75, 0.10), fixed=false,
                   description="ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho_R")

    m <= parameter(:ϵ_p, 10.000, fixed=true,
                   description="ϵ_p: Curvature parameter in the Kimball aggregator for prices.",
                   tex_label="\\epsilon_{p}")

    m <= parameter(:ϵ_w, 10.000, fixed=true,
                   description="ϵ_w: Curvature parameter in the Kimball aggregator for wages.",
                   tex_label="\\epsilon_{w}")

    # exogenous processes - level
    m <= parameter(:γ, 0.3982, (-5.0, 5.0), (-5., 5.), ModelConstructors.Untransformed(), Normal(0.4, 0.1), fixed=false, scaling = x -> x/100,
                   description="γ: The log of the steady-state growth rate of technology.",
                   tex_label="\\gamma")

    m <= parameter(:Lmean, -45., (-1000., 1000.), (-1e3, 1e3), ModelConstructors.Untransformed(), Normal(-45, 5), fixed=false,
                   description="Lmean: Mean level of hours.",
                   tex_label="\\bar{L}")

    m <= parameter(:g_star, 0.1800, fixed=true,
                   description="g_star: 1 - (c_star + i_star)/y_star.",
                   tex_label="g_*")

    # exogenous processes - autocorrelation
    m <= parameter(:ρ_g, 0.9930, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.2), fixed=false,
                   description="ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label="\\rho_g")

    m <= parameter(:ρ_b, 0.2703, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.2), fixed=false,
                   description="ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label="\\rho_b")

    m <= parameter(:ρ_μ, 0.5724, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.2), fixed=false,
                   description="ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label="\\rho_{\\mu}")

    m <= parameter(:ρ_z, 0.9676, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.2), fixed=false,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")

    m <= parameter(:ρ_λ_f, 0.8692, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.2), fixed=false,
                   description="ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label="\\rho_{\\lambda_f}")

    m <= parameter(:ρ_λ_w, 0.9546, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.2), fixed=false,
                   description="ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.", # CHECK THIS
                   tex_label="\\rho_{\\lambda_w}")

    # monetary policy shock - see eqcond
    m <= parameter(:ρ_rm, 0.3000, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.2), fixed=false,
                   description="ρ_rm: AR(1) coefficient in the monetary policy shock process.", # CHECK THIS
                   tex_label="\\rho_{r^m}")

    # exogenous processes - standard deviation
    m <= parameter(:σ_g, 0.6090, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2., 0.10), fixed=false,
                   description="σ_g: The standard deviation of the government spending process.",
                   tex_label="\\sigma_{g}")

    m <= parameter(:σ_b, 0.1818, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2., 0.10), fixed=false,
                   description="σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label="\\sigma_{b}")

    m <= parameter(:σ_μ, 0.4601, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2., 0.10), fixed=false,
                   description="σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label="\\sigma_{\\mu}")

    m <= parameter(:σ_z, 0.4618, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2., 0.10), fixed=false,
                   description="σ_z: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label="\\sigma_{z}")

    m <= parameter(:σ_λ_f, 0.1455, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2., 0.10), fixed=false,
                   description="σ_λ_f: The mean of the process that generates the price elasticity of the composite good. Specifically, the elasticity is (1+λ_{f,t})/(λ_{f_t}).",
                   tex_label="\\sigma_{\\lambda_f}")

    m <= parameter(:σ_λ_w, 0.2089, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2., 0.10), fixed=false,
                   tex_label="\\sigma_{\\lambda_w}")

    m <= parameter(:σ_rm, 0.2397, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2., 0.10), fixed=false,
                   description="σ_r_m: The standard deviation of the monetary policy shock.",
                   tex_label="\\sigma_{r^m}")

    m <= parameter(:η_gz, 0.0500, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.50, 0.20), fixed=false,
                   description="η_gz: Correlate g and z shocks.",
                   tex_label="\\eta_{gz}")

    m <= parameter(:η_λ_f, 0.7652, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.50, 0.20), fixed=false,
                   description="η_λ_f: Moving average component in the price markup shock.",
                   tex_label="\\eta_{\\lambda_f}")

    m <= parameter(:η_λ_w, 0.8936, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.50, 0.20), fixed=false,
                   description="η_λ_w: Moving average component in the wage markup shock.",
                   tex_label="\\eta_{\\lambda_w}")

    m <= parameter(:e_y, 0.0, fixed = true, description = "e_y: Measurement error on GDP", tex_label = "e_y")
    m <= parameter(:e_L, 0.0, fixed = true, description = "e_L: Measurement error on hours worked", tex_label = "e_L")
    m <= parameter(:e_w, 0.0, fixed = true, description = "e_w: Measurement error on wages", tex_label = "e_w")
    m <= parameter(:e_π, 0.0, fixed = true, description = "e_π: Measurement error on GDP deflator", tex_label = "e_π")
    m <= parameter(:e_R, 0.0, fixed = true, description = "e_R: Measurement error on nominal rate of interest", tex_label = "e_R")
    m <= parameter(:e_c, 0.0, fixed = true, description = "e_c: Measurement error on consumption", tex_label = "e_c")
    m <= parameter(:e_i, 0.0, fixed = true, description = "e_i: Measurement error on investment", tex_label = "e_i")

    # steady states
    m <= SteadyStateParameter(:zstar, NaN, description="Steady-state growth rate of productivity", tex_label="\\z_*")
    m <= SteadyStateParameter(:rstar, NaN, tex_label="\\r_*")
    m <= SteadyStateParameter(:Rstarn, NaN, tex_label="\\R_*_n")
    m <= SteadyStateParameter(:rkstar, NaN, tex_label="\\r^k_")
    m <= SteadyStateParameter(:wstar, NaN, tex_label="\\w_*")
    m <= SteadyStateParameter(:Lstar, NaN, tex_label="\\L_*")
    m <= SteadyStateParameter(:kstar, NaN, description="Effective capital that households rent to firms in the steady state.", tex_label="\\k_*")
    m <= SteadyStateParameter(:kbarstar, NaN, description="Total capital owned by households in the steady state.", tex_label="\\bar{k}_*")
    m <= SteadyStateParameter(:istar, NaN, description="Detrended steady-state investment", tex_label="\\i_*")
    m <= SteadyStateParameter(:ystar, NaN, tex_label="\\y_*")
    m <= SteadyStateParameter(:cstar, NaN, tex_label="\\c_*")
    m <= SteadyStateParameter(:wl_c, NaN, tex_label="\\wl_c")
end

"""
```
steadystate!(m::SmetsWouters)
```

Calculates the model's steady-state values. `steadystate!(m)` must be called whenever the parameters of `m` are updated.
"""
function steadystate!(m::SmetsWouters)
    m[:zstar]    = log(1+m[:γ]) + m[:α]/(1-m[:α])*log(m[:Upsilon])
    m[:rstar]    = exp(m[:σ_c]*m[:zstar]) / m[:β]
    m[:Rstarn]   = 100*(m[:rstar]*m[:π_star] - 1)
    m[:rkstar]   = m[:rstar]*m[:Upsilon] - (1-m[:δ])
    m[:wstar]    = (m[:α]^m[:α] * (1-m[:α])^(1-m[:α]) * m[:rkstar]^(-m[:α]) / m[:Φ])^(1/(1-m[:α]))
    m[:Lstar]    = 1.
    m[:kstar]    = (m[:α]/(1-m[:α])) * m[:wstar] * m[:Lstar] / m[:rkstar]
    m[:kbarstar] = m[:kstar] * (1+m[:γ]) * m[:Upsilon]^(1 / (1-m[:α]))
    m[:istar]    = m[:kbarstar] * (1-((1-m[:δ])/((1+m[:γ]) * m[:Upsilon]^(1/(1-m[:α])))))
    m[:ystar]    = (m[:kstar]^m[:α]) * (m[:Lstar]^(1-m[:α])) / m[:Φ]
    m[:cstar]    = (1-m[:g_star])*m[:ystar] - m[:istar]
    m[:wl_c]     = (m[:wstar]*m[:Lstar])/(m[:cstar]*m[:λ_w])

    return m
end

function model_settings!(m::SmetsWouters)

    default_settings!(m)

    # Anticipated shocks
    m <= Setting(:n_mon_anticipated_shocks, 0)
    m <= Setting(:n_mon_anticipated_shocks_padding, 20)

    # Estimation
    m <= Setting(:reoptimize, true)
    m <= Setting(:recalculate_hessian, true)

    # Data
    m <= Setting(:data_vintage, "150827")
    m <= Setting(:data_id, 1, "Dataset identifier")
    m <= Setting(:cond_full_names, [:obs_gdp, :obs_nominalrate])
    m <= Setting(:cond_semi_names, [:obs_nominalrate])
end

"""
```
parameter_groupings(m::SmetsWouters)
```

Returns an `OrderedDict{String, Vector{Parameter}}` mapping descriptions of
parameter groupings (e.g. \"Policy Parameters\") to vectors of
`Parameter`s. This dictionary is passed in as a keyword argument to
`prior_table`.
"""
function parameter_groupings(m::SmetsWouters)
    policy     = [:ψ1, :ψ2, :ψ3, :ρ, :ρ_rm, :σ_rm]
    sticky     = [:ζ_p, :ι_p, :ϵ_p, :ζ_w, :ι_w, :ϵ_w]
    other_endo = [:γ, :α, :β, :σ_c, :h, :ν_l, :δ, :Φ, :S′′, :ppsi,
                  :Lmean, :λ_w, :π_star, :g_star]
    processes  = [:ρ_g, :ρ_b, :ρ_μ, :ρ_z, :ρ_λ_f, :ρ_λ_w, :η_λ_f, :η_λ_w,
                  :σ_g, :σ_b, :σ_μ, :σ_z, :σ_λ_f, :σ_λ_w, :η_gz]

    all_keys     = Vector[policy, sticky, other_endo, processes]
    all_params   = map(keys -> [m[θ]::Parameter for θ in keys], all_keys)
    descriptions = ["Policy", "Nominal Rigidities",
                    "Other Endogenous Propagation and Steady State",
                    "Exogenous Process"]

    groupings = OrderedDict{String, Vector{Parameter}}(zip(descriptions, all_params))

    # Ensure no parameters missing
    incl_params = vcat(collect(values(groupings))...)
    excl_params = [m[θ] for θ in [:Upsilon, :e_y, :e_L, :e_w, :e_π, :e_R, :e_c, :e_i]]
    @assert isempty(setdiff(m.parameters, vcat(incl_params, excl_params)))
    @assert isempty(setdiff(vcat(incl_params, excl_params), m.parameters))

    return groupings
end

"""
```
shock_groupings(m::SmetsWouters)
```

Returns a `Vector{ShockGroup}`, which must be passed in to
`plot_shock_decomposition`. See `?ShockGroup` for details.
"""
function shock_groupings(m::SmetsWouters)
    gov = ShockGroup("g", [:g_sh], RGB(0.70, 0.13, 0.13)) # firebrick
    bet = ShockGroup("b", [:b_sh], RGB(0.3, 0.3, 1.0))
    tfp = ShockGroup("z", [:z_sh], RGB(1.0, 0.55, 0.0)) # darkorange
    pmu = ShockGroup("p-mkp", [:λ_f_sh], RGB(0.60, 0.80, 0.20)) # yellowgreen
    wmu = ShockGroup("w-mkp", [:λ_w_sh], RGB(0.0, 0.5, 0.5)) # teal
    pol = ShockGroup("pol", vcat([:rm_sh], [Symbol("rm_shl$i") for i = 1:n_anticipated_shocks(m)]),
                     RGB(1.0, 0.84, 0.0)) # gold
    mei = ShockGroup("mu", [:μ_sh], :cyan)
    det = ShockGroup("dt", [:dettrend], :gray40)

    return [gov, bet, tfp, pmu, wmu, pol, mei, det]
end