text:
export set_ωpml!, set_boundft!, set_Npml!, set_kbloch!  # basic setter functions
export create_e⁻ⁱᵏᴸ
export clear_objs!  # plural because clearing both electric and magnetic quantities
export create_paramops, create_curls, create_πcmps
export create_iF′ₗgen, create_iHₗgen, create_iEₗgen
export create_βF′ₜgen, create_βHₜgen, create_βEₜgen
export create_A
export complete_fields

# Do not export Model; quality it with the package name MaxwellFDFD, because I would have
# similar types in other packages such as MaxwellSALT and MaxwellGuide.
Base.@kwdef mutable struct Model{K,Kₑ,Kₘ,
                                 Kₑₜ,Kₘₜ,Kₑₗ,Kₘₗ,
                                 K²,
                                 K₊₂,AK₊₂<:AbsArrComplexF{K₊₂}}
    # Frequency
    ωpml::Number = 0.0  # can be complex

    # Grid
    grid::Grid{K}
    cmpₛ::SInt{K}  # Cartesian components of shape dimension
    cmpₑₜ::SInt{Kₑₜ}  # Cartesian components of transverse E-field
    cmpₘₜ::SInt{Kₘₜ}  # Cartesian components of transverse H-field
    cmpₑₗ::SInt{Kₑₗ}  # Cartesian components of longitudinal E-field
    cmpₘₗ::SInt{Kₘₗ}  # Cartesian components of longitudinal H-field

    # Boundary properties
    boundft::SVec{K,FieldType} = SVec(ntuple(k->EE, Val(K)))
    Npml::Tuple2{SInt{K}} = (SVec(ntuple(k->0, Val(K))), SVec(ntuple(k->0, Val(K))))
    kbloch::SFloat{K} = SVec(ntuple(k->0.0, Val(K)))  # [kx_bloch, ky_bloch, kz_bloch]

    # Material parameter arrays
    εₜarr::AK₊₂ = create_param_array(grid.N, ncmp=Kₑₜ)  # filled with zeros
    μₜarr::AK₊₂ = create_param_array(grid.N, ncmp=Kₘₜ)  # filled with zeros

    εₗarr::AK₊₂ = create_param_array(grid.N, ncmp=Kₑ-Kₑₜ)  # filled with zeros
    μₗarr::AK₊₂ = create_param_array(grid.N, ncmp=Kₘ-Kₘₜ)  # filled with zeros

    # Temporary storages for assignment and smoothing of material parameters
    oind2shp::Vector{Shape{K,K²}} = Shape{K,K²}[]
    oind2εind::Vector{ParamInd} = ParamInd[]
    oind2μind::Vector{ParamInd} = ParamInd[]
    εind2ε::Vector{S²ComplexF3} = S²ComplexF3[]
    μind2μ::Vector{S²ComplexF3} = S²ComplexF3[]

    # Boolean flags indicating if the transverse E- and H-field dimensions are orthogonal to
    # the shape dimensions.  Used for material parameter smoothing in calc_matparams!().
    iseₜ˔shp::Bool
    ishₜ˔shp::Bool

    # Indexing scheme for DOFs
    order_cmpfirst::Bool = true
end

function Base.size(mdl::Model{K,Kₑ,Kₘ,Kₑₜ,Kₘₜ}, ft::FieldType) where {K,Kₑ,Kₘ,Kₑₜ,Kₘₜ}
    Kf = ft==EE ? Kₑₜ : Kₘₜ
    sz_grid = size(mdl.grid)
    sz = mdl.order_cmpfirst ? (Kf, sz_grid...) : (sz_grid..., Kf)

    return sz
end

Base.length(mdl::Model, ft::FieldType) = prod(size(mdl, ft))

# Basic setters
set_ωpml!(mdl::Model, ωpml::Number) = (mdl.ωpml = ωpml; nothing)
set_boundft!(mdl::Model{K}, boundft::AbsVec{FieldType}) where {K} = (mdl.boundft = SVec{K}(boundft); nothing)
set_Npml!(mdl::Model{K}, Npml::Tuple2{AbsVecInteger}) where {K} = (mdl.Npml = SVec{K}.(Npml); nothing)
set_kbloch!(mdl::Model{K}, kbloch::AbsVecReal) where {K} = (mdl.kbloch = SVec{K}(kbloch); nothing)

create_e⁻ⁱᵏᴸ(mdl::Model) = exp.(-im .* mdl.kbloch .* mdl.grid.L)

# Main functions
function clear_objs!(mdl::Model)
    mdl.εₜarr .= 0
    mdl.μₜarr .= 0

    mdl.εₗarr .= 0
    mdl.μₗarr .= 0

    empty!(mdl.oind2shp)
    empty!(mdl.oind2εind)
    empty!(mdl.oind2μind)
    empty!(mdl.εind2ε)
    empty!(mdl.μind2μ)

    return nothing
end

MaxwellBase.add_obj!(mdl::Model{K}, matname::String, shapes::AbsVec{<:Shape{K}};
                     ε::MatParam=1.0, μ::MatParam=1.0) where {K} =
    add_obj!(mdl, matname, ε=ε, μ=μ, tuple(shapes...))

function MaxwellBase.add_obj!(mdl::Model{K}, matname::String, shapes::Shape{K}...;
                              ε::MatParam=1.0, μ::MatParam=1.0) where {K}
    mat = Material(matname, ε=ε, μ=μ)
    isbidirectional(EE, mat, mdl) || @error "ε = $ε should not couple transverse and longitudinal fields."
    isbidirectional(HH, mat, mdl) || @error "μ = $μ should not couple transverse and longitudinal fields."

    for shp = shapes  # shapes is tuple
        obj = Object(shp, mat)
        add_obj!(mdl.oind2shp, (mdl.oind2εind,mdl.oind2μind), (mdl.εind2ε,mdl.μind2μ), obj)
    end

    return nothing
end

# Returns true if the material parameter does not couple the transverse and longitudinal
# components of the fields; false otherwise.
function isbidirectional(ft::FieldType, mat::Material, mdl::Model)
    nft = Int(ft)
    cmpₜ = (mdl.cmpₑₜ, mdl.cmpₘₜ)[nft]
    cmpₗ = (mdl.cmpₑₗ, mdl.cmpₘₗ)[nft]

    param = mat.param[nft]

    return iszero(param[cmpₜ,cmpₗ]) && iszero(param[cmpₗ,cmpₜ])
end

function create_stretched_∆ls(mdl::Model)
    ωpml = mdl.ωpml
    grid = mdl.grid
    Npml = mdl.Npml
    s∆l = create_stretched_∆l(ωpml, grid, Npml)
    s∆l⁻¹ = invert_∆l(s∆l)

    boundft = mdl.boundft
    gₑ = ft2gt.(EE, boundft)
    gₘ = ft2gt.(HH, boundft)

    s∆lₑ = t_ind(s∆l, gₑ)  # stretched ∆l's centered at E-field plane locations
    s∆lₘ = t_ind(s∆l, gₘ)  # stretched ∆l's centered at H-field plane locations
    s∆lₑ⁻¹ = t_ind(s∆l⁻¹, gₑ)
    s∆lₘ⁻¹ = t_ind(s∆l⁻¹, gₘ)

    return s∆lₑ, s∆lₘ, s∆lₑ⁻¹, s∆lₘ⁻¹
end

function create_paramops(mdl::Model{K,Kₑ,Kₘ,Kₑₜ,Kₘₜ}) where {K,Kₑ,Kₘ,Kₑₜ,Kₘₜ}
    s∆lₑ, s∆lₘ, s∆lₑ⁻¹, s∆lₘ⁻¹ = create_stretched_∆ls(mdl)
    calc_matparams!(mdl)  # assignment and smoothing; implemented for each specialized alias of Model

    boundft = mdl.boundft
    isbloch = mdl.grid.isbloch
    e⁻ⁱᵏᴸ = create_e⁻ⁱᵏᴸ(mdl)

    isfwd_inₑ = boundft.!=EE
    isfwd_inₘ = boundft.!=HH

    Pεₜ = Kₑₜ==1 ? create_paramop(mdl.εₜarr; mdl.order_cmpfirst) :
                   create_paramop(mdl.εₜarr, isfwd_inₑ, s∆lₘ, s∆lₑ⁻¹, isbloch, e⁻ⁱᵏᴸ; mdl.order_cmpfirst)
    Pμₜ = Kₘₜ==1 ? create_paramop(mdl.μₜarr; mdl.order_cmpfirst) :
                   create_paramop(mdl.μₜarr, isfwd_inₘ, s∆lₑ, s∆lₘ⁻¹, isbloch, e⁻ⁱᵏᴸ; mdl.order_cmpfirst)

    # The following lines work even if the model does not deal with the corresponding
    # longitudinal components.
    Pεₗ = create_paramop(mdl.εₗarr; mdl.order_cmpfirst)  # 0×0 matrix for Kₑₗ = 0 (no longitudinal E-field)
    Pμₗ = create_paramop(mdl.μₗarr; mdl.order_cmpfirst)  # 0×0 matrix for Kₘₗ = 0 (no longitudinal H-field)

    return (Pεₜ,Pμₜ), (Pεₗ,Pμₗ)
end

function create_curls(mdl::Model)
    _, _, s∆lₑ⁻¹, s∆lₘ⁻¹ = create_stretched_∆ls(mdl)

    boundft = mdl.boundft

    cmpₛ = mdl.cmpₛ
    cmpₑₜ, cmpₘₜ = mdl.cmpₑₜ, mdl.cmpₘₜ
    cmpₑₗ, cmpₘₗ = mdl.cmpₑₗ, mdl.cmpₘₗ

    isbloch = mdl.grid.isbloch
    e⁻ⁱᵏᴸ = create_e⁻ⁱᵏᴸ(mdl)

    isfwdₑ = boundft.==EE
    isfwdₘ = boundft.==HH

    # The following lines work even if the model does not deal with the corresponding
    # longitudinal components.
    Cₑₜ = create_curl(isfwdₑ, s∆lₘ⁻¹, isbloch, e⁻ⁱᵏᴸ; cmp_shp=cmpₛ, cmp_out=cmpₘₗ, cmp_in=cmpₑₜ, mdl.order_cmpfirst)
    Cₘₜ = create_curl(isfwdₘ, s∆lₑ⁻¹, isbloch, e⁻ⁱᵏᴸ; cmp_shp=cmpₛ, cmp_out=cmpₑₗ, cmp_in=cmpₘₜ, mdl.order_cmpfirst)

    Cₑₗ = create_curl(isfwdₑ, s∆lₘ⁻¹, isbloch, e⁻ⁱᵏᴸ; cmp_shp=cmpₛ, cmp_out=cmpₘₜ, cmp_in=cmpₑₗ, mdl.order_cmpfirst)
    Cₘₗ = create_curl(isfwdₘ, s∆lₑ⁻¹, isbloch, e⁻ⁱᵏᴸ; cmp_shp=cmpₛ, cmp_out=cmpₑₜ, cmp_in=cmpₘₗ, mdl.order_cmpfirst)

    return (Cₑₜ,Cₘₜ), (Cₑₗ,Cₘₗ)
end

function create_πcmps(mdl::Model{K}) where {K}
    grid = mdl.grid
    cmpₑₜ, cmpₘₜ = mdl.cmpₑₜ, mdl.cmpₘₜ

    permute = reverse(SVec(ntuple(identity, Val(K))))
    scale₀ = SFloat(1,-1)

    πcmpₑ = create_πcmp(grid.N, permute, scale=scale₀[cmpₑₜ]; mdl.order_cmpfirst)
    πcmpₘ = create_πcmp(grid.N, permute, scale=scale₀[cmpₘₜ]; mdl.order_cmpfirst)

    return πcmpₑ, πcmpₘ
end

create_iHₗgen(ω, Ps::Tuple22{AbsMatNumber}, Cs::Tuple22{AbsMatNumber}) = create_iF′ₗgen(EE, ω, Ps, Cs)
create_iEₗgen(ω, Ps::Tuple22{AbsMatNumber}, Cs::Tuple22{AbsMatNumber}) = create_iF′ₗgen(HH, ω, Ps, Cs)

create_iF′ₗgen(ft, ω, Ps::Tuple22{AbsMatNumber}, Cs::Tuple22{AbsMatNumber}) =
    create_iF′ₗgen(ft, ω, Ps[2], Cs[1])

# Create the operator that generates im * F′ₗ from Fₜ, where F′ is the field complementary
# to F (e.g., F′ = H-field for F = E-field).
function create_iF′ₗgen(ft::FieldType,  # type of input field Fₜ, not output field F′ₗ
                        ω::Number,
                        Pₗs::Tuple2{AbsMatNumber},
                        Cₜs::Tuple2{AbsMatNumber})
    nft = Int(ft)
    nft′ = alter(nft)

    α = ft==EE ? -1/ω : 1/ω

    P′ₗ = Pₗs[nft′]
    Cₜ = Cₜs[nft]

    L = α .* (P′ₗ \ Cₜ)

    return L
end

create_βHₜgen(ω, Ps, Cs, πcmps) = create_βF′ₜgen(EE, ω, Ps, Cs, πcmps)
create_βEₜgen(ω, Ps, Cs, πcmps) = create_βF′ₜgen(HH, ω, Ps, Cs, πcmps)

# Create the operator that generates β * F′ₜ, where F′ is the field complementary to F (e.g.,
# F′ = H-field for F = E-field).
create_βF′ₜgen(ft, ω, Ps, Cs, πcmps) = sum(create_βF′ₜgen_comps(ft, ω, Ps, Cs, πcmps))

function create_βF′ₜgen_comps(ft::FieldType,  # type of input field Fₜ, not output field F′ₗ
                              ω::Number,
                              Ps::Tuple22{AbsMatNumber},
                              Cs::Tuple22{AbsMatNumber},
                              πcmps::Tuple2{AbsMatNumber})
    nft = Int(ft)
    nft′ = alter(nft)

    α = ft==EE ? ω : -ω

    Pₜs, Pₗs = Ps
    Cₜs, Cₗs = Cs

    Pₜ = Pₜs[nft]
    C′ₗ = Cₗs[nft′]
    πcmp = πcmps[nft]

    βF′ₜgen₁ = (α .* πcmp) * Pₜ

    if length(Pₗs[nft′]) > 0
        iF′ₗgen = create_iF′ₗgen(ft, ω, Pₗs, Cₜs)
        βF′ₜgen₂ = (πcmp * C′ₗ) * iF′ₗgen
    else
        βF′ₜgen₂ = spzeros(eltype(βF′ₜgen₁), size(βF′ₜgen₁)...)
    end

    return βF′ₜgen₁, βF′ₜgen₂
end

function create_A(ft::FieldType,  # type of input field Fₜ
                  ω::Number,
                  Ps::Tuple22{AbsMatNumber},
                  Cs::Tuple22{AbsMatNumber},
                  πcmps::Tuple2{AbsMatNumber})
    ft′ = alter(ft)

    # Old implementation.  Though more straightforward, this implementation does not
    # gurantee the vector calculus identity βFₜgen₂ * βF′ₜgen₂ = 0 due to rounding errors.
    #
    # βF′ₜgen = create_βF′ₜgen(ft, ω, Ps, Cs, πcmps)
    # βFₜgen = create_βF′ₜgen(ft′, ω, Ps, Cs, πcmps)
    # A = βFₜgen * βF′ₜgen

    βF′ₜgen₁, βF′ₜgen₂ = create_βF′ₜgen_comps(ft, ω, Ps, Cs, πcmps)
    βFₜgen₁, βFₜgen₂ = create_βF′ₜgen_comps(ft′, ω, Ps, Cs, πcmps)

    A = βFₜgen₁ * βF′ₜgen₁ + βFₜgen₁ * βF′ₜgen₂ + βFₜgen₂ * βF′ₜgen₁  # βFₜgen₂ * βF′ₜgen₂ = 0 always

    return A
end

# Implement the following functions for each specific Model.
function calc_matparams! end
function complete_fields end

include("full.jl")
include("te.jl")
include("tm.jl")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "GTK3"
version = v"3.24.11"

# Collection of sources required to build GTK
sources = [
    ArchiveSource("http://ftp.gnome.org/pub/gnome/sources/gtk+/$(version.major).$(version.minor)/gtk+-$(version).tar.xz",
                  "dba7658d0a2e1bfad8260f5210ca02988f233d1d86edacb95eceed7eca982895"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gtk+-*/

# We need to run some commands with a native Glib
apk add glib-dev gtk+3.0

if [[ "${target}" == *-linux-* ]]; then
    # We need to run `wayland-scanner` on the build system only when Wayland is
    # available, i.e. only on Linux
    apk add wayland-dev
fi

# This is awful, I know
ln -sf /usr/bin/glib-compile-resources ${prefix}/bin/glib-compile-resources
ln -sf /usr/bin/glib-compile-schemas ${prefix}/bin/glib-compile-schemas
ln -sf /usr/bin/gdk-pixbuf-pixdata ${prefix}/bin/gdk-pixbuf-pixdata
# Remove wayland-scanner when present, so that we can call the native one
rm -f ${prefix}/bin/wayland-scanner

atomic_patch -p1 $WORKSPACE/srcdir/patches/gdkwindow-quartz_c.patch
atomic_patch -p1 $WORKSPACE/srcdir/patches/meson_build.patch

FLAGS=()
if [[ "${target}" == *-apple-* ]]; then
    FLAGS+=(-Dx11_backend=false -Dwayland_backend=false)
elif [[ "${target}" == *-freebsd* ]]; then
    FLAGS+=(-Dwayland_backend=false)
fi

mkdir build-gtk && cd build-gtk
meson .. \
    -Dintrospection=false \
    -Ddemos=false \
    -Dexamples=false \
    -Dtests=false \
    "${FLAGS[@]}" \
    --cross-file="${MESON_TARGET_TOOLCHAIN}"
ninja -j${nproc}
ninja install

# Remove temporary links
rm ${prefix}/bin/gdk-pixbuf-pixdata ${prefix}/bin/glib-compile-{resources,schemas}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = Product[
    LibraryProduct("libgailutil-3", :libgailutil3),
    LibraryProduct("libgdk-3", :libgdk3),
    LibraryProduct("libgtk-3", :libgtk3),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"),
    Dependency("Glib_jll"),
    Dependency("Cairo_jll"),
    Dependency("Pango_jll"),
    Dependency("FriBidi_jll"),
    Dependency("FreeType2_jll"),
    Dependency("gdk_pixbuf_jll"),
    Dependency("Libepoxy_jll"),
    Dependency("ATK_jll"),
    Dependency("HarfBuzz_jll"),
    Dependency("xkbcommon_jll"),
    Dependency("iso_codes_jll"),
    Dependency("Wayland_jll"),
    Dependency("Xorg_libXrandr_jll"),
    Dependency("Xorg_libX11_jll"),
    Dependency("Xorg_libXrender_jll"),
    Dependency("Xorg_libXi_jll"),
    Dependency("Xorg_libXext_jll"),
    Dependency("Xorg_libXcursor_jll"),
    Dependency("Xorg_libXdamage_jll"),
    Dependency("Xorg_libXfixes_jll"),
    Dependency("Xorg_libXcomposite_jll"),
    Dependency("Xorg_libXinerama_jll"),
    Dependency("Fontconfig_jll"),
    Dependency("at_spi2_atk_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

ext: jl


text:
import CIM
test_files_folder = joinpath(@__DIR__, "CIM_files", "CGMES_v2.4.15_MiniGridTestConfiguration_BaseCase_v3")
@show test_files_folder
EQ_file = CIM.PyCIM.import_file(joinpath(test_files_folder, "MiniGridTestConfiguration_BC_EQ_v3.0.0.xml"))
#SSH_file = CIM.PyCIM.import_file(joinpath(test_files_folder, "MiniGridTestConfiguration_BC_SSH_v3.0.0.xml"))
@test length(EQ_file) == 233

EQ_file = CIM.CIMpy.import_file(joinpath(test_files_folder, "MiniGridTestConfiguration_BC_EQ_v3.0.0.xml"))
@test length(EQ_file["topology"]) == 238

merged_file = CIM.CIMpy.import_file(joinpath.(test_files_folder, readdir(test_files_folder)))
@test length(merged_file["topology"]) == 436


text:
using PSCOPF

using Test
using Dates

@testset "test_check_timesteps" begin

    @testset "definition_example" begin
        TS = [DateTime("2015-01-01T11:00:00"),
            DateTime("2015-01-01T11:15:00"),
            DateTime("2015-01-01T11:30:00"),
            DateTime("2015-01-01T11:45:00")]

        @test PSCOPF.check_target_timepoints(TS)
    end

    @testset "not_sorted" begin
        TS = [DateTime("2015-01-01T11:15:00"),
            DateTime("2015-01-01T11:00:00"),
            DateTime("2015-01-01T11:30:00"),
            DateTime("2015-01-01T11:45:00")]

        @test !PSCOPF.check_target_timepoints(TS)
    end

    @testset "duplicate_value" begin
        TS = [DateTime("2015-01-01T11:00:00"),
            DateTime("2015-01-01T11:00:00"),
            DateTime("2015-01-01T11:30:00"),
            DateTime("2015-01-01T11:45:00")]

        @test !PSCOPF.check_target_timepoints(TS)
    end

    #=
    # FIXME Do we assume it is 1 hour with 4 equally spaced Timesteps ?
    # not for now, to allow simple test cases
    # i.e.:
    #       length(TS) == 4
    #       TS[i] - TS[i-1] == 15 mins
    @testset "duplicate_value" begin
        TS = [DateTime("2015-01-01T11:00:00"),
            DateTime("2015-01-01T11:00:00"),
            DateTime("2015-01-01T11:30:00"),
            DateTime("2015-01-01T11:45:00")]

        @test !PSCOPF.check_target_timepoints(TS)
    end
    =#

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
export tbaa_make_child, tbaa_addrspace

function tbaa_make_child(name::String; constant::Bool=false, ctx::LLVM.Context)
    tbaa_root = MDNode([MDString("custom_tbaa"; ctx)]; ctx)
    tbaa_struct_type =
        MDNode([MDString("custom_tbaa_$name"; ctx),
                tbaa_root,
                ConstantInt(0; ctx)]; ctx)
    tbaa_access_tag =
        MDNode([tbaa_struct_type,
                tbaa_struct_type,
                ConstantInt(0; ctx),
                ConstantInt(constant ? 1 : 0; ctx)]; ctx)

    return tbaa_access_tag
end

tbaa_addrspace(as; ctx::LLVM.Context) = tbaa_make_child("addrspace($(as))"; ctx)


text:
"""
    SumBasis(b1, b2...)

Similar to [`CompositeBasis`](@ref) but for the [`directsum`](@ref) (⊕)
"""
struct SumBasis{S,B} <: Basis
    shape::S
    bases::B
end
SumBasis(bases) = SumBasis(Int[length(b) for b in bases], bases)
SumBasis(shape, bases::Vector) = (tmp = (bases...,); SumBasis(shape, tmp))
SumBasis(bases::Vector) = SumBasis((bases...,))
SumBasis(bases::Basis...) = SumBasis((bases...,))

==(b1::T, b2::T) where T<:SumBasis = equal_shape(b1.shape, b2.shape)
==(b1::SumBasis, b2::SumBasis) = false
length(b::SumBasis) = sum(b.shape)

"""
    directsum(b1::Basis, b2::Basis)

Construct the [`SumBasis`](@ref) out of two sub-bases.
"""
directsum(b1::Basis, b2::Basis) = SumBasis(Int[length(b1); length(b2)], Basis[b1, b2])
directsum(b::Basis) = b
directsum(b::Basis...) = reduce(directsum, b)
function directsum(b1::SumBasis, b2::Basis)
    shape = [b1.shape;length(b2)]
    bases = [b1.bases...;b2]
    return SumBasis(shape, (bases...,))
end
function directsum(b1::Basis, b2::SumBasis)
    shape = [length(b1);b2.shape]
    bases = [b1;b2.bases...]
    return SumBasis(shape, (bases...,))
end
function directsum(b1::SumBasis, b2::SumBasis)
    shape = [b1.shape;b2.shape]
    bases = [b1.bases...;b2.bases...]
    return SumBasis(shape, (bases...,))
end

const ⊕ = directsum

"""
    directsum(x::Ket, y::Ket)

Construct a spinor via the [`directsum`](@ref) of two [`Ket`](@ref)s.
The result is a [`Ket`](@ref) with data given by `[x.data;y.data]` and its
basis given by the corresponding [`SumBasis`](@ref).
**NOTE**: The resulting state is not normalized!
"""
directsum(x::Ket, y::Ket) = Ket(directsum(x.basis, y.basis), [x.data; y.data])
directsum(x::StateVector...) = reduce(directsum, x)

"""
    getblock(x::Ket{<:SumBasis}, i)

For a [`Ket`](@ref) defined on a [`SumBasis`](@ref), get the state as it is defined
on the ith sub-basis.
"""
function getblock(x::Ket{B}, i) where B<:SumBasis
    b_i = x.basis.bases[i]
    inds = cumsum([0;length.(x.basis.bases[1:i])...])
    return Ket(b_i, x.data[inds[i]+1:inds[i+1]])
end

"""
    setblock!(x::Ket{<:SumBasis}, val::Ket, i)

Set the data of `x` on the ith sub-basis equal to the data of `val`.
"""
function setblock!(x::Ket{B}, val::Ket, i) where B<:SumBasis
    check_samebases(x.basis.bases[i], val)
    inds = cumsum([0;length.(x.basis.bases[1:i])...])
    x.data[inds[i]+1:inds[i+1]].data .= val.data
    return x
end

"""
    directsum(x::DataOperator, y::DataOperator)

Compute the direct sum of two operators. The result is an operator on the
corresponding [`SumBasis`](@ref).
"""
function directsum(a::DataOperator, b::DataOperator)
    dType = promote_type(eltype(a),eltype(b))
    data = zeros(dType, size(a,1)+size(b,1), size(a,2)+size(b,2))
    data[1:size(a,1),1:size(a,2)] = a.data
    data[size(a,1)+1:end, size(a,2)+1:end] = b.data
    return Operator(directsum(a.basis_l, b.basis_l), directsum(a.basis_r, b.basis_r), data)
end
function directsum(a::SparseOpType, b::SparseOpType)
    dType = promote_type(eltype(a),eltype(b))
    data = spzeros(dType, size(a,1)+size(b,1), size(a,2)+size(b,2))
    data[1:size(a,1),1:size(a,2)] = a.data
    data[size(a,1)+1:end, size(a,2)+1:end] = b.data
    return Operator(directsum(a.basis_l, b.basis_l), directsum(a.basis_r, b.basis_r), data)
end
directsum(a::AbstractOperator...) = reduce(directsum, a)

"""
    setblock!(op::DataOperator{<:SumBasis,<:SumBasis}, val::DataOperator, i, j)

Set the data of `op` corresponding to the block `(i,j)` equal to the data of `val`.
"""
function setblock!(op::DataOperator{<:SumBasis,<:SumBasis}, val::DataOperator, i, j)
    (bases_l,bases_r) = op.basis_l.bases, op.basis_r.bases
    check_samebases(bases_l[i], val.basis_l)
    check_samebases(bases_r[j], val.basis_r)
    inds_i = cumsum([0;length.(bases_l[1:i])...])
    inds_j = cumsum([0;length.(bases_r[1:j])...])
    op.data[inds_i[i]+1:inds_i[i+1],inds_j[j]+1:inds_j[j+1]] = val.data
    return op
end

"""
    getblock(op::Operator{<:SumBasis,<:SumBasis}, i, j)

Get the sub-basis operator corresponding to the block `(i,j)` of `op`.
"""
function getblock(op::DataOperator{BL,BR}, i, j) where {BL<:SumBasis,BR<:SumBasis}
    (bases_l,bases_r) = op.basis_l.bases, op.basis_r.bases
    inds_i = cumsum([0;length.(bases_l[1:i])...])
    inds_j = cumsum([0;length.(bases_r[1:j])...])
    data = op.data[inds_i[i]+1:inds_i[i+1],inds_j[j]+1:inds_j[j+1]]
    return Operator(bases_l[i],bases_r[j], data)
end

"""
    LazyDirectSum <: AbstractOperator

Lazy implementation of `directsum`
"""
mutable struct LazyDirectSum{BL,BR,T} <: AbstractOperator{BL,BR}
    basis_l::BL
    basis_r::BR
    operators::T
end

# Methods
LazyDirectSum(op1::AbstractOperator, op2::AbstractOperator) = LazyDirectSum(directsum(op1.basis_l,op2.basis_l),directsum(op1.basis_r,op2.basis_r),(op1,op2))
LazyDirectSum(op1::LazyDirectSum, op2::AbstractOperator) = LazyDirectSum(directsum(op1.basis_l,op2.basis_l),directsum(op1.basis_r,op2.basis_r),(op1.operators...,op2))
LazyDirectSum(op1::AbstractOperator, op2::LazyDirectSum) = LazyDirectSum(directsum(op1.basis_l,op2.basis_l),directsum(op1.basis_r,op2.basis_r),(op1,op2.operators...))
LazyDirectSum(op1::LazyDirectSum, op2::LazyDirectSum) = LazyDirectSum(directsum(op1.basis_l,op2.basis_l),directsum(op1.basis_r,op2.basis_r),(op1.operators...,op2.operators...))
LazyDirectSum(op::AbstractOperator...) = reduce(LazyDirectSum, op)

# Algebra
dense(x::LazyDirectSum) = directsum(dense.(x.operators)...)
*(op1::LazyDirectSum{B1,B2},op2::LazyDirectSum{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis}= LazyDirectSum(op1.basis_l,op2.basis_r,Tuple(op1.operators[i]*op2.operators[i] for i=1:length(op1.operators)))
-(op::LazyDirectSum) = LazyDirectSum([-op.operators[1];op.operators[2:end]]...)
+(op1::LazyDirectSum{B1,B2},op2::LazyDirectSum{B1,B2}) where {B1<:Basis,B2<:Basis} = LazyDirectSum((op1.operators .+ op2.operators)...)
dagger(op::LazyDirectSum) = LazyDirectSum(op.basis_r, op.basis_l, dagger.(op.operators))
Base.eltype(x::LazyDirectSum) = promote_type(eltype.(x.operators)...)

directsum(op1::AbstractOperator,op2::LazyDirectSum) = LazyDirectSum(op1,op2)
directsum(op1::LazyDirectSum,op2::AbstractOperator) = LazyDirectSum(op1,op2)
directsum(op1::LazyDirectSum,op2::LazyDirectSum) = LazyDirectSum(op1,op2)

# Use lazy sum for FFTOperator
transform(b1::SumBasis,b2::SumBasis; kwargs...) = LazyDirectSum([transform(b1.bases[i],b2.bases[i];kwargs...) for i=1:length(b1.bases)]...)
directsum(op1::FFTOperator, op2::FFTOperator) = LazyDirectSum(op1,op2)

# Fast in-place multiplication
function mul!(result::Ket{B1},M::LazyDirectSum{B1,B2},b::Ket{B2},alpha_,beta_) where {B1,B2}
    alpha = convert(ComplexF64, alpha_)
    beta = convert(ComplexF64, beta_)
    bases_r = b.basis.bases
    bases_l = M.basis_l.bases
    index = cumsum([0;length.(bases_r)...])
    for i=1:length(index)-1
        tmpket = Ket(bases_r[i], b.data[index[i]+1:index[i+1]])
        tmpresult = Ket(bases_l[i], result.data[index[i]+1:index[i+1]])
        mul!(tmpresult,M.operators[i],tmpket,alpha,beta)
        result.data[index[i]+1:index[i+1]] = tmpresult.data
    end
    return result
end
function mul!(result::Bra{B2},b::Bra{B1},M::LazyDirectSum{B1,B2},alpha_,beta_) where {B1,B2}
    alpha = convert(ComplexF64, alpha_)
    beta = convert(ComplexF64, beta_)
    bases_l = b.basis.bases
    bases_r = M.basis_r.bases
    index = cumsum([0;length.(bases_r)...])
    for i=1:length(index)-1
        tmpbra = Bra(bases_l[i], b.data[index[i]+1:index[i+1]])
        tmpresult = Bra(bases_r[i], result.data[index[i]+1:index[i+1]])
        mul!(tmpresult,tmpbra,M.operators[i],alpha,beta)
        result.data[index[i]+1:index[i+1]] = tmpresult.data
    end
    return result
end


text:
export TrackEvent, MetaEvent, MIDIEvent, SysexEvent

"""
    TrackEvent <: Any
Abstract supertype for all MIDI events.

All track events begin with a variable length time value
(see [`readvariablelength`](@ref)) and have a field named `dT` which contains it.
This number notes after how many ticks since the last event does the current
even takes place.

`MIDIEvent`s then resume with a MIDI channel message
defined in `constants.jl`. They're followed by 1 or 2 bytes, depending on the
channel message (see `MIDI.EVENTTYPETOLENGTH`).
If no valid channel message is identified,
the previous seen channel message is used. After that the MIDI command is encoded.

`MetaEvent`s and `SysexEvent`s both resume with a specific byte (see `constants.jl`).
"""
abstract type TrackEvent end

#####################################################################################
# MetaEvent
#####################################################################################
"""
    MetaEvent <: TrackEvent
See [`TrackEvent`](@ref).
"""
mutable struct MetaEvent <: TrackEvent
    dT::Int
    metatype::UInt8
    data::Array{UInt8,1}
end

@inline ismetaevent(b::UInt8) = b == 0xFF

function readmetaevent(dT::Int, f::IO)
    # Meta events are 0xFF - type (1 byte) - variable length data length - data bytes
    skip(f, 1) # Skip the 0xff that starts the event
    metatype = read(f, UInt8)
    datalength = readvariablelength(f)
    data = read!(f, Array{UInt8}(undef, datalength))

    MetaEvent(dT, metatype, data)
end

function writeevent(f::IO, event::MetaEvent)
    if event.dT < 0
        error("Negative deltas are not allowed. Please reorder your events.")
    end
    writevariablelength(f, event.dT)
    write(f, META)
    write(f, event.metatype)
    writevariablelength(f, convert(Int, length(event.data)))
    write(f, event.data)
end

#####################################################################################
# MidiEvent
#####################################################################################
"""
    MIDIEvent <: TrackEvent
See [`TrackEvent`](@ref).
"""
mutable struct MIDIEvent <: TrackEvent
    dT::Int
    status::UInt8
    data::Array{UInt8,1}
end

function isstatusbyte(b::UInt8)
    (b & 0b10000000) == 0b10000000
end

function isdatabyte(b::UInt8)
    !isstatusbyte(b)
end

function isMIDIevent(b::UInt8)
    !ismetaevent(b) && !issysexevent(b)
end

function channelnumber(m::MIDIEvent)
    0x0F & m.status
end

function readMIDIevent(dT::Int, f::IO, laststatus::UInt8)
    statusbyte = read(f, UInt8)
    highnybble = statusbyte & 0b11110000

    toread = 0
    if haskey(EVENTTYPETOLENGTH, highnybble)
        laststatus = statusbyte
        toread = EVENTTYPETOLENGTH[highnybble]
    else # Running status is in use
        toread = EVENTTYPETOLENGTH[laststatus & 0b11110000]
        statusbyte = laststatus
        skip(f, -1)
    end

    data = read!(f, Array{UInt8}(undef, toread))

    MIDIEvent(dT, statusbyte, data)
end

function writeevent(f::IO, event::MIDIEvent, writestatus::Bool)
    event.dT < 0 && error(
    "Negative deltas are not allowed. Please reorder your events.")

    writevariablelength(f, event.dT)

    if writestatus
        write(f, event.status)
    end

    write(f, event.data)
end

function writeevent(f::IO, event::MIDIEvent)
    writeevent(f, event, true)
end

#####################################################################################
# SysexEvent
#####################################################################################
"""
    SysexEvent <: TrackEvent
See [`TrackEvent`](@ref).
"""
mutable struct SysexEvent <: TrackEvent
    dT::Int
    data::Array{UInt8,1}
end

issysexevent(b::UInt8) = b == 0xF0

function readsysexevent(dT::Int, f::IO)
    data = UInt8[]
    read(f, UInt8) # Eat the SYSEX that's on top of f
    datalength = readvariablelength(f)
    b = read(f, UInt8)
    while isdatabyte(b)
        push!(data, b)
        b = read(f, UInt8)
    end

    if b != 0xF7
        error("Invalid sysex event, did not end with 0xF7")
    end
    # The last byte of sysex event is F7. We leave that out of the data, and write it back when we write the event.

    if length(data) + 1 != datalength
        error("Invalid sysex event. Expected $(datalength) bytes, received $(length(data) + 1)")
    end
    SysexEvent(dT, data)
end

function writeevent(f::IO, event::SysexEvent)
    if event.dT < 0
        error("Negative deltas are not allowed. Please reorder your events.")
    end

    writevariablelength(f, event.dT)
    write(f, SYSEX)
    writevariablelength(f, length(event.data) + 1) # +1 for the ending F7
    write(f, event.data)
    write(f, 0xF7)
end


text:
"""
    SeisSWCGP(in,operator,param,iter,lp,imax)

Robust Stagewise weak Orthognal Matching Pursuit algorithm with the coefficient updated by the conjugate gradient pursuit

# Arguments
- `dn`: Noisy data
- `operator`: Any operator passed the dot test
- `param`: parameters for the operator
- `iter`: Max iterations for the matching pursuit
- `a`: parameter alpha used to calcualte the threshold value
- `lp`: lp-norm space inner product
- `imax`: Max interation for the lp-norm inner product

"""
function SeisSWCGP(dn,operator,param,iter,a,lp,imax)
    m=operator(dn,param,-1,lp,imax);
    r=dn
    n=1
    ind_r=[];
    ind_c=[];
    v=[];
    d=[];
    dr=zeros(size(dn));
    for i=1:iter
        g=zeros(size(m))
        gn=operator(r,param,-1,lp,imax)
        thres=a*maximum(abs.(vec(gn)));
        ind=findall(abs.(gn).>= thres)
        ind_ri=getindex.(ind, 1)
        ind_ci=getindex.(ind, 2)
        push!(ind_r,ind_ri)
        push!(ind_c,ind_ci)
        for j=length(ind_r)
            g[ind_r[j],ind_c[j]]=gn[ind_r[j],ind_c[j]]
        end
        if n==1
            d=g;
            v=operator(d,param,1,lp,imax);
        else
            w=operator(g,param,1,lp,imax)
            u=-LP_norm(vec(v),vec(w),lp,imax,1)
            d=g+u*d
            v=w+u*v
        end
        beta=LP_norm(vec(v),vec(r),lp,imax,1)
        dr=dr+beta*v
        r=r-beta*v
        n=n+1
    end
    return dr
end

max_stars_repo_name: Ji-seismic/SeisTest

max_issues_repo_name: Ji-seismic/SeisTest

max_forks_repo_name: Ji-seismic/SeisTest


text:
##
using Pkg, Revise, BenchmarkTools
Pkg.activate(".")
##
import Images
import NPZ: npzread
using GLMakie
using ZStackAnalysis
##
file_content = npzread("./data/Fed_X63_Z2_SIM.npz")
#file_content = npzread("./SIM_Jing.npz")

zstack = file_content["data"] |> Preprocessing.convert_to_float32 |> Preprocessing.equalize_zstack
scale = file_content["scale"]
# TODO: create a struct and/or use NamedArrays / AxisArrays to store channel names & relevant scale
# Also make all the downstream functions dispatch on this struct

##
size(zstack)
##
xy_region = (1201:1500, 201:500)
actin = interpolate_zstack(zstack[1, :, xy_region...], scale)
desmin = interpolate_zstack(zstack[2, :, xy_region...], scale)

##
volume(actin, algorithm = :mip, colormap = :Greens_9, transparancy = true)
volume(desmin, algorithm = :mip, colormap = :Reds_9, transparancy = true)
current_figure()
##
z_ind = 50
actin_I, actin_θ = LineFilterTransform(actin[z_ind, :, :], 10, 7, 20)
desmin_I, desmin_θ = LineFilterTransform(desmin[z_ind, :, :], 10, 7, 20)
actin_oft = OrientationFilterTransform(actin_I, actin_θ, 10, 7, 20)
desmin_oft = OrientationFilterTransform(desmin_I, desmin_θ, 10, 7, 20)
##
binarize(x) = x > zero(x) ? one(x) : zero(x)
oft_mask =

##
using FileIO
jing = load("./Image 6_Out crop.tif")
jing_reshaped = reshape(jing, 2560, 2560, 4, 21)
size(jing_reshaped)
##

##
size(zstack)
##

##

image(jing_reshaped[:, :, 3, 20])
##
jing_overlay = Float32.(mean(jing_reshaped, dims = 4)[1000:1500, 1000:1500, :, 1])
image(Images.adjust_histogram(jing_overlay[:, :,3], Images.Equalization(nbins = 256, minval = 0.0, maxval = 0.5)))
##
typeof(jing_overlay)

##
image(@view zstack[2, 1, :, :])
##
vol = zstack[1, :, 1201:1500, 201:500]
image(vol[1, :, :])
##
axes(zstack, 1)

##
xy_region = (1201:1500, 201:500)
actin = interpolate_zstack(zstack[1, :, xy_region...], scale)
desmin = interpolate_zstack(zstack[2, :, xy_region...], scale)
##
jing3D = interpolate_zstack(Float32.(jing_reshaped[1000:1500, 1000:1500, :, :]), [3.13e-8, 3.13e-8, 1.1e-7])
size(jing3D)
##

##

##
jing_I, jing_θ = LineFilterTransform(jing_overlay[:, :, 3], 10, 15, 20)
##
jing_oft = OrientationFilterTransform(jing_I, jing_θ, 10, 15, 20)
##
angle_to_color(angle) = Images.colorsigned(Images.RGB(1.0, 0.0, 0.0), Images.RGB(0.5, 0.0, 0.5), Images.RGB(0.0, 0.0, 1.0))(Images.scalesigned(π/2)(angle))
##
fig = Figure(resolution = (1500, 1000))
ax1 = fig[1, 1] = GLMakie.Axis(fig, title = "Raw Image")
img1 = image!(ax1, jing_overlay[:, :, 1])
ax2 = fig[2, 1] = GLMakie.Axis(fig, title = "Line Filter Transform (intensity)")
img2 = image!(ax2, Images.scaleminmax(0, maximum(jing_I)).(jing_I))
ax3 = fig[2, 2] = GLMakie.Axis(fig, title = "Line Filter Transform (orientation)")

img3 = image!(ax3, angle_to_color.(jing_θ))
ax4 = fig[1, 2] = GLMakie.Axis(fig, title = "Orientation Filter Transform")
img4 = image!(ax4, Images.scaleminmax(0, maximum(jing_oft)).(jing_oft))
fig
##

##
fig = Figure(resolution = (1000, 700))
ax1 = fig[1, 1] = GLMakie.Axis(fig, title = "Actin (Equalized)")
img1 = image!(ax1, actin[z_ind, :, :])
ax2 = fig[2, 1] = GLMakie.Axis(fig, title = "Actin (LFT intensity)")
img2 = image!(ax2, Images.scaleminmax(0, maximum(actin_I)).(actin_I))
ax3 = fig[2, 2] = GLMakie.Axis(fig, title = "Actin (LFT orientation)")
img3 = image!(ax3, angle_to_color.(actin_θ))
ax4 = fig[1, 2] = GLMakie.Axis(fig, title = "Actin (OFT)")
img4 = image!(ax4, Images.scaleminmax(0, maximum(actin_oft)).(actin_oft))
fig
##
fig = Figure(resolution = (1000, 700))
ax1 = fig[1, 1] = GLMakie.Axis(fig, title = "Desmin (Equalized)")
img1 = image!(ax1, desmin[z_ind, :, :])
ax2 = fig[2, 1] = GLMakie.Axis(fig, title = "Desmin (LFT intensity)")
img2 = image!(ax2, Images.scaleminmax(0, maximum(desmin_I)).(desmin_I))
ax3 = fig[2, 2] = GLMakie.Axis(fig, title = "Desmin (LFT orientation)")
img3 = image!(ax3, angle_to_color.(desmin_θ))
ax4 = fig[1, 2] = GLMakie.Axis(fig, title = "Desmin (OFT)")
img4 = image!(ax4, Images.scaleminmax(0, maximum(desmin_oft)).(desmin_oft))
fig
##
fig = Figure(resolution = (800, 700))
ax1 = fig[1, 1] = GLMakie.Axis(fig, title = "Overlay (Equalized)")
img1 = image!(ax1, Images.colorview(Images.RGB, desmin, actin, Images.zeroarray)[z_ind, :, :])
ax2 = fig[2, 1] = GLMakie.Axis(fig, title = "Actin (LFT intensity)")
img2 = image!(ax2, Images.colorview(Images.RGB, Images.scaleminmax(0, maximum(desmin_I)).(desmin_I), Images.scaleminmax(0, maximum(actin_I)).(actin_I), Images.zeroarray))
fig
##
import ImageBinarization: binarize, Otsu
lft_bin = binarize(desmin_I, Otsu())
image(lft_bin)
##
size(actin)
##
analysis_region = (1:100, 1:200, 1:200)
fig = Figure(resolution = (1500, 1000))
ax1 = fig[1, 1] = Axis3(fig, title = "Actin 3D stack")
ax2 = fig[1, 2] = Axis3(fig, title = "Desmin 3D stack")
volume!(ax1, actin[analysis_region...], algorithm = :mip, colormap = :Greens_9, transparancy = true)
volume!(ax2, desmin[analysis_region...], algorithm = :mip, colormap = :Reds_9, transparancy = true)
fig
##
actin3dI, actin3dθ = LineFilterTransform(actin[analysis_region...], 10, 5, 10)
##
desmin3dI, desmin3dθ = LineFilterTransform(desmin[analysis_region...], 10, 5, 10)
##
fig = Figure(resolution = (1000, 700))
ax1 = fig[1, 1] = Axis3(fig, title = "Actin LFT intensity")
ax2 = fig[1, 2] = Axis3(fig, title = "Desmin LFT intensity")
volume!(ax1, Images.scaleminmax(0, maximum(actin3dI)).(actin3dI), algorithm = :mip, colormap = :Greens_9, transparancy = true)
volume!(ax2, Images.scaleminmax(0, maximum(desmin3dI)).(desmin3dI), algorithm = :mip, colormap = :Reds_9, transparancy = true)
fig
##
desmin3doft = OrientationFilterTransform(desmin3dI, desmin3dθ, 10, 5, 10)
##
actin3doft = OrientationFilterTransform(actin3dI, actin3dθ, 10, 5, 10)
##
##
fig = Figure(resolution = (1000, 700))
ax1 = fig[1, 1] = Axis3(fig, title = "Actin OFT intensity")
ax2 = fig[1, 2] = Axis3(fig, title = "Desmin OFT intensity")
volume!(ax1, Images.scaleminmax(0.0, maximum(actin3doft)).(actin3doft), algorithm = :mip, colormap = :Greens_9, transparancy = true)
volume!(ax2, Images.scaleminmax(0.0, maximum(desmin3doft)).(desmin3doft), algorithm = :mip, colormap = :Reds_9, transparancy = true)
fig
##
fig = Figure(resolution = (800, 700))
ax1 = fig[1, 1] = GLMakie.Axis(fig, title = "Overlay (Equalized)")
img1 = image!(ax1, Images.colorview(Images.RGB, desmin_adj, actin_adj, Images.zeroarray)[z_ind, :, :])
ax2 = fig[2, 1] = GLMakie.Axis(fig, title = "Actin (LFT intensity)")
img2 = image!(ax2, Images.colorview(Images.RGB, desmin3dI[z_ind, :, :], actin3dI[z_ind, :, :], Images.zeroarray))
fig


text:
"""
`init_subspec!(m::HetDSGE)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::HetDSGE)
    if subspec(m) == "ss0"
        return
    # All but g_shock
    elseif subspec(m) == "ss1"
        return ss1!(m)
    # All but b_shock
    elseif subspec(m) == "ss2"
        return ss2!(m)
    # All but μ shock
    elseif subspec(m) == "ss3"
        return ss3!(m)
    # All but z shock
    elseif subspec(m) == "ss4"
        return ss4!(m)
    # All but λ_f shock
    elseif subspec(m) == "ss5"
        return ss5!(m)
    # All but λ_w shock
    elseif subspec(m) == "ss6"
        return ss6!(m)
    # All but rm_shock
    elseif subspec(m) == "ss7"
        return ss7!(m)
    # All but steady state parameters
    elseif subspec(m) == "ss8"
        return ss8!(m)
    else
        error("This subspec should be a 0")
    end
end

"""
```
ss1!(m::HetDSGE)
```

Initializes subspec 1 of `HetDSGE`.
This shuts down all shocks except the government spending process (g shock).
"""
function ss1!(m::HetDSGE)
    m <= parameter(:spp, 4., (-15., 15.), (-15., 15.), Untransformed(),
                   Normal(4., 1.5), fixed = true,
                   description = "S'': The second derivative of households' cost of adjusting investment.",
                   tex_label = "S''")
    m <= parameter(:ρR , 0.75, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   BetaAlt(0.75, 0.10), fixed = true,
                   description = "ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho_R")
    m <= parameter(:ψπ , 1.5, (1e-5, 10.), (1e-5, 10.0), ModelConstructors.Exponential(),
                   Normal(1.5, 0.25), fixed = true,
                   description = "ψ1: Weight on inflation gap in monetary policy rule.",
                   tex_label = "\\psi_1")

    m <= parameter(:ρ_G, 0.5, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_g")

    m <= parameter(:ρ_B, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label = "\\rho_B")
    m <= parameter(:ρ_μ, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label = "\\rho_{\\mu}")
    m <= parameter(:ρ_z, 0.0,(1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_z: AR(1) coefficient in the technology process.",
                   tex_label = "\\rho_z")
    m <= parameter(:ρ_lamf, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_f}")
    m <= parameter(:ρ_lamw, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_w}")
    m <= parameter(:ρ_mon, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_rm: AR(1) coefficient in the monetary policy shock process.",
                   tex_label = "\\rho_{r^m}")

    m <= parameter(:σ_b, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label = "\\sigma_{b}")
    m <= parameter(:σ_μ, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label = "\\sigma_{\\mu}")
    m <= parameter(:σ_z, 0.0,(1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_z: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label = "\\sigma_z")
    m <= parameter(:σ_λ_f, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   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.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_λ_w",
                   tex_label = "\\sigma_{\\lambda_w}")
    m <= parameter(:σ_rm, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_r_m: The standard deviation of the monetary policy shock.", tex_label = "\\sigma_{r^m}")

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

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

    m <= parameter(:κ_p, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3),
                 fixed = true, description = "κ_p : The slope of the Price Phillips curve",
                 tex_label = "\\kappa")

    m <= parameter(:κ_w, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3),
                   fixed = true, description = "κ_w: The slope of the Wage Phillips curve",
                   tex_label = "\\kappa")

    m <= parameter(:ψy , 0.5, (-0.5, 0.5), (-0.5, 0.5), Untransformed(),
                   Normal(0.12, 0.05), fixed = true,
                   description = "ψy: Weight on output growth in monetary policy rule")

    m <= parameter(:ϕh, 2., (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(),
                   Normal(2, 0.75), fixed = true, description = "inverse frisch elasticity",
                   tex_label = "\\phi_h")
end

"""
```
ss2!(m::HetDSGE)
```

Initializes subspec 2 of `HetDSGE`.
This shuts down all shocks except the preference shifter process (b shock).
"""
function ss2!(m::HetDSGE)
    m <= parameter(:spp, 4., (-15., 15.), (-15., 15.), Untransformed(),
                   Normal(4., 1.5), fixed = true,
                   description = "S'': The second derivative of households' cost of adjusting investment.",
                   tex_label = "S''")
    m <= parameter(:ρR , 0.75, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   BetaAlt(0.75, 0.10), fixed = true,
                   description = "ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho_R")
    m <= parameter(:ψπ , 1.5, (1e-5, 10.), (1e-5, 10.0), ModelConstructors.Exponential(),
                   Normal(1.5, 0.25), fixed = true,
                   description = "ψ1: Weight on inflation gap in monetary policy rule.",
                   tex_label = "\\psi_1")

    m <= parameter(:ρ_G, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_g")

    m <= parameter(:ρ_B, 0.5, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_b: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_b")

    m <= parameter(:ρ_μ, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label = "\\rho_{\\mu}")
    m <= parameter(:ρ_z, 0.0,(1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_z: AR(1) coefficient in the technology process.",
                   tex_label = "\\rho_z")
    m <= parameter(:ρ_lamf, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_f}")
    m <= parameter(:ρ_lamw, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_w}")
    m <= parameter(:ρ_mon, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_rm: AR(1) coefficient in the monetary policy shock process.",
                   tex_label = "\\rho_{r^m}")

    m <= parameter(:σ_g, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_g: The standard deviation of the government spending process.",
                   tex_label = "\\sigma_{g}")
    m <= parameter(:σ_μ, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label = "\\sigma_{\\mu}")
    m <= parameter(:σ_z, 0.0,(1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_z: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label = "\\sigma_z")
    m <= parameter(:σ_λ_f, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   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.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_λ_w",
                   tex_label = "\\sigma_{\\lambda_w}")
    m <= parameter(:σ_rm, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_r_m: The standard deviation of the monetary policy shock.", tex_label = "\\sigma_{r^m}")

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

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

    m <= parameter(:κ_p, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3),
                 fixed = true, description = "κ_p : The slope of the Price Phillips curve",
                 tex_label = "\\kappa")

    m <= parameter(:κ_w, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3),
                   fixed = true, description = "κ_w: The slope of the Wage Phillips curve",
                   tex_label = "\\kappa")

    m <= parameter(:ψy , 0.5, (-0.5, 0.5), (-0.5, 0.5), Untransformed(),
                   Normal(0.12, 0.05), fixed = true,
                   description = "ψy: Weight on output growth in monetary policy rule")

    m <= parameter(:ϕh, 2., (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(),
                   Normal(2, 0.75), fixed = true, description = "inverse frisch elasticity",
                   tex_label = "\\phi_h")
end

"""
```
ss3!(m::HetDSGE)
```

Initializes subspec 3 of `HetDSGE`.
This shuts down all shocks except the capital adjustment cost process (μ shock).
"""
function ss3!(m::HetDSGE)
    m <= parameter(:spp, 4., (-15., 15.), (-15., 15.), Untransformed(),
                   Normal(4., 1.5), fixed = true,
                   description = "S'': The second derivative of households' cost of adjusting investment.",
                   tex_label = "S''")
    m <= parameter(:ρR , 0.75, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   BetaAlt(0.75, 0.10), fixed = true,
                   description = "ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho_R")
    m <= parameter(:ψπ , 1.5, (1e-5, 10.), (1e-5, 10.0), ModelConstructors.Exponential(),
                   Normal(1.5, 0.25), fixed = true,
                   description = "ψ1: Weight on inflation gap in monetary policy rule.",
                   tex_label = "\\psi_1")

    m <= parameter(:ρ_G, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_g")
    m <= parameter(:ρ_B, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label = "\\rho_B")

    m <= parameter(:ρ_μ, 0.5, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_μ: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_\\mu")

    m <= parameter(:ρ_z, 0.0,(1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_z: AR(1) coefficient in the technology process.",
                   tex_label = "\\rho_z")
    m <= parameter(:ρ_lamf, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_f}")
    m <= parameter(:ρ_lamw, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_w}")
    m <= parameter(:ρ_mon, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_rm: AR(1) coefficient in the monetary policy shock process.",
                   tex_label = "\\rho_{r^m}")

    m <= parameter(:σ_g, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_g: The standard deviation of the government spending process.",
                   tex_label = "\\sigma_{g}")
    m <= parameter(:σ_b, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label = "\\sigma_{b}")
    m <= parameter(:σ_z, 0.0,(1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_z: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label = "\\sigma_z")
    m <= parameter(:σ_λ_f, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   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.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_λ_w",
                   tex_label = "\\sigma_{\\lambda_w}")
    m <= parameter(:σ_rm, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_r_m: The standard deviation of the monetary policy shock.", tex_label = "\\sigma_{r^m}")

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

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

    m <= parameter(:κ_p, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3),
                 fixed = true, description = "κ_p : The slope of the Price Phillips curve",
                 tex_label = "\\kappa")

    m <= parameter(:κ_w, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3),
                   fixed = true, description = "κ_w: The slope of the Wage Phillips curve",
                   tex_label = "\\kappa")

    m <= parameter(:ψy , 0.5, (-0.5, 0.5), (-0.5, 0.5), Untransformed(),
                   Normal(0.12, 0.05), fixed = true,
                   description = "ψy: Weight on output growth in monetary policy rule")

    m <= parameter(:ϕh, 2., (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(),
                   Normal(2, 0.75), fixed = true, description = "inverse frisch elasticity",
                   tex_label = "\\phi_h")
end

"""
```
ss4(m::HetDSGE)
```

Initializes subspec 4 of `HetDSGE`.
This shuts down all shocks except the technology process (z shock).
"""
function ss4!(m::HetDSGE)
    m <= parameter(:spp, 4., (-15., 15.), (-15., 15.), Untransformed(),
                   Normal(4., 1.5), fixed = true,
                   description = "S'': The second derivative of households' cost of adjusting investment.",
                   tex_label = "S''")
    m <= parameter(:ρR , 0.75, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   BetaAlt(0.75, 0.10), fixed = true,
                   description = "ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho_R")
    m <= parameter(:ψπ , 1.5, (1e-5, 10.), (1e-5, 10.0), ModelConstructors.Exponential(),
                   Normal(1.5, 0.25), fixed = true,
                   description = "ψ1: Weight on inflation gap in monetary policy rule.",
                   tex_label = "\\psi_1")

    m <= parameter(:ρ_G, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_g")
    m <= parameter(:ρ_B, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label = "\\rho_B")
    m <= parameter(:ρ_μ, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label = "\\rho_{\\mu}")

    m <= parameter(:ρ_z, 0.5, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_z: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_z")

    m <= parameter(:ρ_lamf, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_f}")
    m <= parameter(:ρ_lamw, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_w}")
    m <= parameter(:ρ_mon, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_rm: AR(1) coefficient in the monetary policy shock process.",
                   tex_label = "\\rho_{r^m}")

    m <= parameter(:σ_g, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_g: The standard deviation of the government spending process.",
                   tex_label = "\\sigma_{g}")
    m <= parameter(:σ_b, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label = "\\sigma_{b}")
    m <= parameter(:σ_μ, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label = "\\sigma_{\\mu}")
    m <= parameter(:σ_λ_f, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   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.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_λ_w",
                   tex_label = "\\sigma_{\\lambda_w}")
    m <= parameter(:σ_rm, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_r_m: The standard deviation of the monetary policy shock.", tex_label = "\\sigma_{r^m}")

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

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

    m <= parameter(:κ_p, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3),
                 fixed = true, description = "κ_p : The slope of the Price Phillips curve",
                 tex_label = "\\kappa")

    m <= parameter(:κ_w, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3),
                   fixed = true, description = "κ_w: The slope of the Wage Phillips curve",
                   tex_label = "\\kappa")

    m <= parameter(:ψy , 0.5, (-0.5, 0.5), (-0.5, 0.5), Untransformed(),
                   Normal(0.12, 0.05), fixed = true,
                   description = "ψy: Weight on output growth in monetary policy rule")

    m <= parameter(:ϕh, 2., (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(),
                   Normal(2, 0.75), fixed = true, description = "inverse frisch elasticity",
                   tex_label = "\\phi_h")

end

"""
```
ss5!(m::HetDSGE)
```

Initializes subspec 5 of `HetDSGE`.
This shuts down all shocks except the price mark-up shock process (λ_f shock).
"""
function ss5!(m::HetDSGE)
    m <= parameter(:spp, 4., (-15., 15.), (-15., 15.), Untransformed(),
                   Normal(4., 1.5), fixed = true,
                   description = "S'': The second derivative of households' cost of adjusting investment.",
                   tex_label = "S''")
    m <= parameter(:ρR , 0.75, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   BetaAlt(0.75, 0.10), fixed = true,
                   description = "ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho_R")
    m <= parameter(:ψπ , 1.5, (1e-5, 10.), (1e-5, 10.0), ModelConstructors.Exponential(),
                   Normal(1.5, 0.25), fixed = true,
                   description = "ψ1: Weight on inflation gap in monetary policy rule.",
                   tex_label = "\\psi_1")

    m <= parameter(:ρ_G, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_g")
    m <= parameter(:ρ_B, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label = "\\rho_B")
    m <= parameter(:ρ_μ, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label = "\\rho_{\\mu}")
    m <= parameter(:ρ_z, 0.0,(1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_z: AR(1) coefficient in the technology process.",
                   tex_label = "\\rho_z")

    m <= parameter(:ρ_lamf, 0.5, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_lamf: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_{\\lambda_f}")

    m <= parameter(:ρ_lamw, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_w}")
    m <= parameter(:ρ_mon, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_rm: AR(1) coefficient in the monetary policy shock process.",
                   tex_label = "\\rho_{r^m}")

    m <= parameter(:σ_g, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_g: The standard deviation of the government spending process.",
                   tex_label = "\\sigma_{g}")
    m <= parameter(:σ_b, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label = "\\sigma_{b}")
    m <= parameter(:σ_μ, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label = "\\sigma_{\\mu}")
    m <= parameter(:σ_z, 0.0,(1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_z: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label = "\\sigma_z")
    m <= parameter(:σ_λ_w, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_λ_w",
                   tex_label = "\\sigma_{\\lambda_w}")
    m <= parameter(:σ_rm, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_r_m: The standard deviation of the monetary policy shock.", tex_label = "\\sigma_{r^m}")

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

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

    m <= parameter(:κ_p, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3),
                 fixed = true, description = "κ_p : The slope of the Price Phillips curve",
                 tex_label = "\\kappa")

    m <= parameter(:κ_w, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3),
                   fixed = true, description = "κ_w: The slope of the Wage Phillips curve",
                   tex_label = "\\kappa")

    m <= parameter(:ψy , 0.5, (-0.5, 0.5), (-0.5, 0.5), Untransformed(),
                   Normal(0.12, 0.05), fixed = true,
                   description = "ψy: Weight on output growth in monetary policy rule")

    m <= parameter(:ϕh, 2., (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(),
                   Normal(2, 0.75), fixed = true, description = "inverse frisch elasticity",
                   tex_label = "\\phi_h")

end

"""
```
ss6!(m::HetDSGE)
```

Initializes subspec 6 of `HetDSGE`.
This shuts down all shocks except the wage mark-up shock process (λ_w shock).
"""
function ss6!(m::HetDSGE)
    m <= parameter(:spp, 4., (-15., 15.), (-15., 15.), Untransformed(),
                   Normal(4., 1.5), fixed = true,
                   description = "S'': The second derivative of households' cost of adjusting investment.",
                   tex_label = "S''")
    m <= parameter(:ρR , 0.75, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   BetaAlt(0.75, 0.10), fixed = true,
                   description = "ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho_R")
    m <= parameter(:ψπ , 1.5, (1e-5, 10.), (1e-5, 10.0), ModelConstructors.Exponential(),
                   Normal(1.5, 0.25), fixed = true,
                   description = "ψ1: Weight on inflation gap in monetary policy rule.",
                   tex_label = "\\psi_1")

    m <= parameter(:ρ_G, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_g")
    m <= parameter(:ρ_B, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label = "\\rho_B")
    m <= parameter(:ρ_μ, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label = "\\rho_{\\mu}")
    m <= parameter(:ρ_z, 0.0,(1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_z: AR(1) coefficient in the technology process.",
                   tex_label = "\\rho_z")
    m <= parameter(:ρ_lamf, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_f}")
    m <= parameter(:ρ_lamw, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_λ_w: AR(1) coefficient in the price mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_w}")

    m <= parameter(:ρ_mon, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_rm: AR(1) coefficient in the monetary policy shock process.",
                   tex_label = "\\rho_{r^m}")

    m <= parameter(:σ_g, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_g: The standard deviation of the government spending process.",
                   tex_label = "\\sigma_{g}")
    m <= parameter(:σ_b, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label = "\\sigma_{b}")
    m <= parameter(:σ_μ, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label = "\\sigma_{\\mu}")
    m <= parameter(:σ_z, 0.0,(1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_z: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label = "\\sigma_z")
    m <= parameter(:σ_λ_f, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   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(:σ_rm, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_r_m: The standard deviation of the monetary policy shock.", tex_label = "\\sigma_{r^m}")

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

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

    m <= parameter(:κ_p, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3),
                 fixed = true, description = "κ_p : The slope of the Price Phillips curve",
                 tex_label = "\\kappa")

    m <= parameter(:κ_w, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3),
                   fixed = true, description = "κ_w: The slope of the Wage Phillips curve",
                   tex_label = "\\kappa")

    m <= parameter(:ψy , 0.5, (-0.5, 0.5), (-0.5, 0.5), Untransformed(),
                   Normal(0.12, 0.05), fixed = true,
                   description = "ψy: Weight on output growth in monetary policy rule")

    m <= parameter(:ϕh, 2., (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(),
                   Normal(2, 0.75), fixed = true, description = "inverse frisch elasticity",
                   tex_label = "\\phi_h")

end

"""
```
ss7!(m::HetDSGE)
```

Initializes subspec 7 of `HetDSGE`.
This shuts down all shocks except the monetary policy shock process (rm shock).
"""
function ss7!(m::HetDSGE)
    m <= parameter(:spp, 4., (-15., 15.), (-15., 15.), Untransformed(),
                   Normal(4., 1.5), fixed = true,
                   description = "S'': The second derivative of households' cost of adjusting investment.",
                   tex_label = "S''")
    m <= parameter(:ρR , 0.75, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   BetaAlt(0.75, 0.10), fixed = true,
                   description = "ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho_R")
    m <= parameter(:ψπ , 1.5, (1e-5, 10.), (1e-5, 10.0), ModelConstructors.Exponential(),
                   Normal(1.5, 0.25), fixed = true,
                   description = "ψ1: Weight on inflation gap in monetary policy rule.",
                   tex_label = "\\psi_1")

    m <= parameter(:ρ_G, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_g")
    m <= parameter(:ρ_B, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label = "\\rho_B")
    m <= parameter(:ρ_μ, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label = "\\rho_{\\mu}")
    m <= parameter(:ρ_z, 0.0,(1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_z: AR(1) coefficient in the technology process.",
                   tex_label = "\\rho_z")
    m <= parameter(:ρ_lamf, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_f}")
    m <= parameter(:ρ_lamw, 0.0, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_w}")

    m <= parameter(:ρ_mon, 0.5, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_rm: AR(1) coefficient in the monetary policy shock process.",
                   tex_label = "\\rho_{r^m}")

    m <= parameter(:σ_g, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_g: The standard deviation of the government spending process.",
                   tex_label = "\\sigma_{g}")
    m <= parameter(:σ_b, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label = "\\sigma_{b}")
    m <= parameter(:σ_μ, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label = "\\sigma_{\\mu}")
    m <= parameter(:σ_z, 0.0,(1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_z: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label = "\\sigma_z")
    m <= parameter(:σ_λ_f, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   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.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description = "σ_λ_w",
                   tex_label = "\\sigma_{\\lambda_w}")

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

    m <= parameter(:Lmean, -45.9364, (-1000., 1000.), (-1e3, 1e3), Untransformed(),
                   Normal(-45., 5.), fixed=true,
                   description="Lmean: Mean level of hours.",
                   tex_label="\\bar{L}")
    m <= parameter(:κ_p, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3),
                 fixed = true, description = "κ_p : The slope of the Price Phillips curve",
                 tex_label = "\\kappa")

    m <= parameter(:κ_w, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3),
                   fixed = true, description = "κ_w: The slope of the Wage Phillips curve",
                   tex_label = "\\kappa")

    m <= parameter(:ψy , 0.5, (-0.5, 0.5), (-0.5, 0.5), Untransformed(),
                   Normal(0.12, 0.05), fixed = true,
                   description = "ψy: Weight on output growth in monetary policy rule")

    m <= parameter(:ϕh, 2., (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(),
                   Normal(2, 0.75), fixed = true, description = "inverse frisch elasticity",
                   tex_label = "\\phi_h")
end

"""
```
ss8!(m::HetDSGE)
```

Initializes subspec 8 of `HetDSGE`.
Estimates all non-steadystate parameters.
Right now it's empty becaus we do this by default
"""
function ss8!(m::HetDSGE)

end


text:
# power analysis functions for mean for normal distribution with known variance

"""
    powerNormTest(d::Real = 0, n::Real = 0, alpha::Float64 = 0.05, alternative = "two.sided")

Compute power of a sample of `n` observations to test the effect size `d`
at type I error = `alpha` (default: 0.05). The effect size d = μ - μ₀. The option
`alternative` is `two.sided` for H₀: μ = μ₀ vs H₁: μ ≠ μ₀, `less` for H₁: μ < μ₀, and `greater`
for H₁: μ > μ₀.
"""
function powerNormTest(;
    d::Real = 0.0,
    n::Real = 0,
    alpha::Float64 = 0.05,
    alternative::String = "two")

    check_args(d=d,n=n,alpha=alpha)

    if alternative == "less"
        tside = 1
    elseif alternative in ("two","two.sided","two-sided")
        tside = 2
        d = abs(d)
    elseif alternative == "greater"
        tside = 3
    end

    if tside == 1
        return cdf(Normal(),quantile(Normal(),alpha) - d*sqrt(n))
    elseif tside == 2
        return ccdf(Normal(),cquantile(Normal(),alpha/2) - d*sqrt(n)) + cdf(Normal(),quantile(Normal(),alpha/2) - d*sqrt(n))
    elseif tside == 3
        return ccdf(Normal(),cquantile(Normal(),alpha) - d*sqrt(n))
    end
end

"""
    samplesizeNormTest(d::Real = 0, alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute the sample size that can test the effect size `d` with power > `power` (default: 0.8)
at type I error = `alpha` (default: 0.05). The effect size d = μ - μ₀. The option
`alternative` is `two.sided` for H₀: μ = μ₀ vs H₁: μ ≠ μ₀, `less` for H₁: μ < μ₀, and `greater`
for H₁: μ > μ₀.
"""
function samplesizeNormTest(;
    d::Real = 0.0,
    alpha::Real = 0.05,
    power::Real = 0.0,
    alternative::String = "two")

    check_args(d=d,alpha=alpha,power=power)

    return ceil(Int64,fzero(x->powerNormTest(n = x, d = d, alpha = alpha, alternative = alternative) - power, 2.0, 10.0^7))
end

"""
    effectsizeNormTest(alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute the effect size that a sample size `n` can test with power > `power` (default: 0.8)
at type I error = `alpha` (default: 0.05). The option
`alternative` is `two.sided` for H₀: μ = μ₀ vs H₁: μ ≠ μ₀, `less` for H₁: μ < μ₀, and `greater`
for H₁: μ > μ₀.
"""
function effectsizeNormTest(;
    n::Real = 0,
    alpha::Real = 0.05,
    power::Real = 0.0,
    alternative::String = "two"
    )

    check_args(n=n,alpha=alpha,power=power)

    return fzero(x -> powerNormTest(n = n, d = x, alpha = alpha, alternative = alternative) - power,.001,100)
end

"""
    alphaNormTest(d::Real = 0, n::Real = 0, power::Float64 = 0.8, alternative = "two.sided")

Compute the probability of type I error that a sample size `n` has in testing
the effect size `d` with power > `power` (default: 0.8) at type I error = `alpha` (default: 0.05).
The effect size d = μ - μ₀. The option `alternative` is `two.sided` for H₀: μ = μ₀ vs H₁: μ ≠ μ₀,
`less` for H₁: μ < μ₀, and `greater` for H₁: μ > μ₀.
"""
function alphaNormTest(;
    n::Real = 0,
    d::Real = 0.0,
    power::Real = 0.0,
    alternative::String = "two"
    )

    check_args(d=d,n=n,power=power)

    return fzero(x->powerNormTest(n = n, d = d, alpha = x, alternative = alternative) - power, 1e-10, 1 - 1e-10)
end

"""
    pwr.NormTest(d::Real = 0, n::Real = 0, alpha::Float64 = 0.0, power::Float64 = 0.0, alternative = "two.sided")

Compute one of the test parameters such as sample size (`n`),
effect size (`d`), type I error (`alpha`), and power (`power`).
The parameter to be estimated must be set to zero (default).
The effect size d = μ - μ₀. The option `alternative` is `two.sided` for H₀: μ = μ₀ vs H₁: μ ≠ μ₀,
`less` for H₁: μ < μ₀, and `greater` for H₁: μ > μ₀.
"""
function NormTest(;
    n::Real = 0,
    d::Real = 0.0,
    alpha::Real = 0.05,
    power::Real = 0.0,
    alternative::String = "two")

    if sum([x == 0 for x in (n,d,alpha,power)]) != 1
        error("exactly one of n, d, power, and alpha must be zero")
    end

    if power == 0.0
        power = powerNormTest(n = n, d = d, alpha = alpha, alternative = alternative)
    elseif alpha == 0.0
        alpha = alphaNormTest(n = n, d = d, power = power, alternative = alternative)
    elseif d == 0.0
        d = effectsizeNormTest(n = n, alpha = alpha, power = power, alternative = alternative)
    elseif n == 0
        n = samplesizeNormTest(d = d, alpha = alpha, power = power, alternative = alternative)
    end

    alt = Dict("two" => "two-sided", "two.sided" => "two-sided", "less" => "less", "greater" => "greater")

    return htest(
        string("Mean power calculation for normal distribution with known variance"),
        OrderedDict(
            "n" => n,
            "d" => d,
            "alpha" => alpha,
            "power" => power,
            "alternative" => alt[alternative]
            )
        )
end

max_stars_repo_name: mwsohn/pwr

max_issues_repo_name: mwsohn/pwr

max_forks_repo_name: mwsohn/pwr


text:
using Images, IJulia, Plots
gr() # FASTER!

grad = ColorGradient(:grays)
# grad = ColorGradient(RGBA{Float64}[RGBA{Float64}(0.05,0.05,0.05,1.0), RGBA{Float64}(0.95,0.95,0.95,1.0)], [0.0, 1.0])
heatmap(res["last"],
    c = grad, # seems to be the same as :grays, tbh.
    legend = nothing,
    xticks = nothing,
    yticks = nothing,
    aspect_ratio = 1.0)

function test_anim4()
    anim = @animate for i in 1:150
        plot(1:i, 1:i, color_palette = :grays)
    end

    gif(anim, "line_anim.gif", fps = 15)
    gif(anim, "line_anim.mp4", fps = 15)
    mp4(anim, "line_anim_.mp4", fps = 15)
    gif(anim, "line_anim.mov", fps = 15)
end

"""
Julia animate version
"""
function test_anim(curGrid)
    # define correct backend for plotting
    pyplot()
    # create animation object that can be saved as mp4 or displayed as gif
    anim = @animate for i=1:20
        #plt[:figure](); plt[:axis]("off"); plt[:imshow](curGrid, cmap=cm.Greys_r, vmin = 0.0, vmax = 1.0)
        Plots.plot(rand(5,5),linewidth=2,title="My Plot")
    end

    gif(anim, "test_anim.rand.gif", fps = 3)

    anim = @animate for i=1:20
        #plt[:figure](); plt[:axis]("off"); plt[:imshow](curGrid, cmap=cm.Greys_r, vmin = 0.0, vmax = 1.0)
        Plots.plot(curGrid, linewidth=2, title="My Plot")
    end

    gif(anim, "test_anim.curGrid.gif", fps = 3)
    return(anim);
end

function test_anim2(curGrid)
    function init()
    end

    function animate(i)
        im[:set_data](curGrid)
        plt[:draw]()
    end

    # create base figure
    fig = figure("MyFigure", figsize=(512, 512))

    myanim = anim.FuncAnimation(im, animate, init_func = init, frames = 100, interval=20)
    return(myanim)
end

# Does not seem to work for the real deal
function test_anim3(curGrid; nframes = 20, save_type = "mp4")
    function showanim(filename)
        base64_video = base64encode(open(filename))
        display("text/html", """<video controls src="data:video/x-m4v;base64,$base64_video">""")
    end

    # set up writer for movie / image files
    if save_type == "mp4"
        f_ending = ".mp4"
        #writer = anim.FFMpegWriter(fps=15, extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])
        writer = anim.FFMpegWriter(fps=15)
    elseif save_type == "gif"
        f_ending = ".gif"
        #writer = "imagemagick"
        writer = anim.ImageMagickWriter()
    else
        throw(error("Wrong save type!"))
    end

    fig = figure(figsize=(4,4))
    axis("off")

    function make_frame(i)
        #curGrid = min.(1, curGrid .* i)
        imshow(curGrid, cmap=cm.Greys_r, vmin = 0.0, vmax = 1.0)
    end

    withfig(fig) do
        myanim = anim.FuncAnimation(fig, make_frame, frames=nframes, interval=20)
        myanim[:save]("test3" * f_ending, dpi = 100, writer=writer)
    end
    close()
end

max_stars_repo_name: Xentrics/SmoothLifeJulia

max_issues_repo_name: Xentrics/SmoothLifeJulia

max_forks_repo_name: Xentrics/SmoothLifeJulia


text:
using Pkg

unregistered = [
   ("https://github.com/sdwfrost/Gillespie.jl","master"),
   ("https://github.com/augustinas1/MomentClosure.jl","main")
]

for u in unregistered
  Pkg.add(PackageSpec(url=u[1], rev=u[2]))
end


text:
module Test_meutil

using LinearAlgebra
using MatrixEquations
using Test

n = 500
ar, ur, = schur(rand(n,n))
cr = rand(n,n)
Qr = cr'*cr

@testset "Matrix Equations Utilities" begin
x = copy(Qr); @time utqu!(x,ur);
@time y = ur'*Qr*ur;
@test norm(x-y) < sqrt(eps(1.))

x = copy(Qr); @time utqu!(x,ur');
@time y = ur*Qr*ur';
@test norm(x-y) < sqrt(eps(1.))

@time x = utqu(Qr,ur);
@time y = ur'*Qr*ur;
@test norm(x-y) < sqrt(eps(1.))

@time x = utqu(Qr,ur');
@time y = ur*Qr*ur';
@test norm(x-y) < sqrt(eps(1.))

@time x = vec2triu(triu2vec(Qr, rowwise = false, her = false),rowwise = false, her = false)
@test x == triu(Qr)

@time x = vec2triu(triu2vec(Qr, rowwise = false, her = true),rowwise = false, her = true)
@test x == Qr

@time x = vec2triu(triu2vec(Qr, rowwise = true, her = false),rowwise = true, her = false)
@test x == triu(Qr)

@time x = vec2triu(triu2vec(Qr, rowwise = true, her = true),rowwise = true, her = true)
@test x == Qr

@test isschur(rand(2,3)) == false &&
      isschur(rand(1,1)) == true &&
      isschur(rand(3,3)) == false &&
      isschur(rand(3,3),rand(2,2)) == false

 u = [1 2;3 4]
 q = [1 3; 3 1]

try
    utqu(q,u)
    utqu(q//2,u)
    utqu(q//2,u//2)
    utqu!(q,u)
    utqu!(q//2,u)
    utqu(q//2,u//2)
   @test true
catch
    @test false
end

try
    utqu!(q,u//2)
    @test false
catch
    @test true
end

try
    utqu!(rand(3,3),rand(3,3))
    @test false
catch
    @test true
end

try
    utqu!(Symmetric(rand(3,3)),rand(2,2))
    @test false
catch
    @test true
end

try
    utqu(rand(3,3),rand(3,3))
    @test false
catch
    @test true
end

try
    utqu(Symmetric(rand(3,3)),rand(2,2))
    @test false
catch
    @test true
end

try
    utqu(Symmetric(rand(3,3)),rand(2,2)')
    @test false
catch
    @test true
end

try
    qrupdate!(triu(rand(4,4)), rand(2,3))
    @test false
catch
    @test true
end

try
    rqupdate!(triu(rand(4,4)), rand(2,3))
    @test false
catch
    @test true
end

try
    triu2vec(rand(2,2), her = true)
    @test false
catch
    @test true
end

end

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module dSFMT

import Base: copy, copy!, ==

export DSFMT_state, dsfmt_get_min_array_size, dsfmt_get_idstring,
       dsfmt_init_gen_rand, dsfmt_init_by_array, dsfmt_gv_init_by_array,
       dsfmt_fill_array_close_open!, dsfmt_fill_array_close1_open2!,
       win32_SystemFunction036!

"Mersenne Exponent"
const MEXP = 19937
"DSFMT internal state array size of N 128-bit integers."
const N    = floor(Int, ((MEXP - 128) / 104 + 1))
"""Julia DSFMT state representation size counted in 32-bit integers.

Size: (DSFMT state array of Int128 + 1)*4 + Int32 index + Int32 padding
"""
const JN32 = (N+1)*4+1+1
"Jump polynomial for 10^20 steps for dSFMT with exponent 19937"
const JPOLY1e21  = "[encoded data, 5046 characters]"

mutable struct DSFMT_state
    val::Vector{Int32}

    DSFMT_state(val::Vector{Int32} = zeros(Int32, JN32)) =
        new(length(val) == JN32 ? val : throw(DomainError()))
end

copy!(dst::DSFMT_state, src::DSFMT_state) = (copy!(dst.val, src.val); dst)
copy(src::DSFMT_state) = DSFMT_state(copy(src.val))

==(s1::DSFMT_state, s2::DSFMT_state) = s1.val == s2.val

function dsfmt_get_idstring()
    idstring = ccall((:dsfmt_get_idstring,:libdSFMT),
                     Ptr{UInt8},
                     ())
    return unsafe_string(idstring)
end

function dsfmt_get_min_array_size()
    min_array_size = ccall((:dsfmt_get_min_array_size,:libdSFMT),
                           Int32,
                           ())
end

const dsfmt_min_array_size = dsfmt_get_min_array_size()

function dsfmt_init_gen_rand(s::DSFMT_state, seed::UInt32)
    ccall((:dsfmt_init_gen_rand,:libdSFMT),
          Void,
          (Ptr{Void}, UInt32,),
          s.val, seed)
end

function dsfmt_init_by_array(s::DSFMT_state, seed::Vector{UInt32})
    ccall((:dsfmt_init_by_array,:libdSFMT),
          Void,
          (Ptr{Void}, Ptr{UInt32}, Int32),
          s.val, seed, length(seed))
end

function dsfmt_gv_init_by_array(seed::Vector{UInt32})
    ccall((:dsfmt_gv_init_by_array,:libdSFMT),
          Void,
          (Ptr{UInt32}, Int32),
          seed, length(seed))
end

function dsfmt_fill_array_close1_open2!(s::DSFMT_state, A::Ptr{Float64}, n::Int)
    @assert Csize_t(A) % 16 == 0 # the underlying C array must be 16-byte aligned
    @assert dsfmt_min_array_size <= n && iseven(n)
    ccall((:dsfmt_fill_array_close1_open2,:libdSFMT),
          Void,
          (Ptr{Void}, Ptr{Float64}, Int),
          s.val, A, n)
end

function dsfmt_fill_array_close_open!(s::DSFMT_state, A::Ptr{Float64}, n::Int)
    @assert Csize_t(A) % 16 == 0 # the underlying C array must be 16-byte aligned
    @assert dsfmt_min_array_size <= n && iseven(n)
    ccall((:dsfmt_fill_array_close_open,:libdSFMT),
          Void,
          (Ptr{Void}, Ptr{Float64}, Int),
          s.val, A, n)
end

# dSFMT jump
function dsfmt_jump(s::DSFMT_state, jp::AbstractString)
    val = s.val
    nval = length(val)
    index = val[nval - 1]
    work = zeros(UInt64, JN32 >> 1)
    dsfmt = Vector{UInt64}(nval >> 1)
    ccall(:memcpy, Ptr{Void}, (Ptr{UInt64}, Ptr{Int32}, Csize_t),
          dsfmt, val, (nval - 1) * sizeof(Int32))
    dsfmt[end] = UInt64(N*2)

    for c in jp
        bits = parse(UInt8,c,16)
        for j in 1:4
            (bits & 0x01) != 0x00 && dsfmt_jump_add!(work, dsfmt)
            bits = bits >> 0x01
            dsfmt_jump_next_state!(dsfmt)
        end
    end

    work[end] = index
    return DSFMT_state(reinterpret(Int32, work))
end

function dsfmt_jump_add!(dest::Vector{UInt64}, src::Vector{UInt64})
    dp = dest[end] >> 1
    sp = src[end] >> 1
    diff = ((sp - dp + N) % N)
    i = 1
    while i <= N-diff
        j = i*2-1
        p = j + diff*2
        dest[j]   ⊻= src[p]
        dest[j+1] ⊻= src[p+1]
        i += 1
    end
    while i <= N
        j = i*2-1
        p = j + (diff - N)*2
        dest[j]   ⊻= src[p]
        dest[j+1] ⊻= src[p+1]
        i += 1
    end
    dest[N*2+1] ⊻= src[N*2+1]
    dest[N*2+2] ⊻= src[N*2+2]
    return dest
end

function dsfmt_jump_next_state!(mts::Vector{UInt64})
    POS1 = 117
    SL1  = 19
    SR   = 12
    MSK1 = 0x000ffafffffffb3f
    MSK2 = 0x000ffdfffc90fffd

    idx = (mts[end] >> 1) % N

    a = idx*2+1
    b = ((idx + POS1) % N)*2+1
    u = N*2+1

    t0 = mts[a]
    t1 = mts[a+1]
    L0 = mts[u]
    L1 = mts[u+1]
    mts[u]   = xor(t0 << SL1, L1 >> 32, L1 << 32, mts[b])
    mts[u+1] = xor(t1 << SL1, L0 >> 32, L0 << 32, mts[b+1])
    mts[a]   = xor(mts[u]   >> SR, mts[u]   & MSK1, t0)
    mts[a+1] = xor(mts[u+1] >> SR, mts[u+1] & MSK2, t1)

    mts[end] = (mts[end] + 2) % (N*2)
    return mts
end

## Windows entropy

if is_windows()
    function win32_SystemFunction036!(a::Array)
        ccall((:SystemFunction036, :Advapi32), stdcall, UInt8, (Ptr{Void}, UInt32), a, sizeof(a))
    end
end

end # module

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
# ---
# title: 686. Repeated String Match
# id: problem686
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/repeated-string-match/description/>
# hidden: true
# ---
#
# Given two strings `a` and `b`, return the minimum number of times you should
# repeat string `a` so that string `b` is a substring of it. If it is impossible
# for `b` to be a substring of `a` after repeating it, return `-1`.
#
# **Notice:**  string `"abc"` repeated 0 times is `""`,  repeated 1 time is
# `"abc"` and repeated 2 times is `"abcabc"`.
#
#
#
# **Example 1:**
#
#
#
#     Input: a = "abcd", b = "cdabcdab"
#     Output: 3
#     Explanation: We return 3 because by repeating a three times "ab **cdabcdab** cd", b is a substring of it.
#
#
# **Example 2:**
#
#
#
#     Input: a = "a", b = "aa"
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: a = "a", b = "a"
#     Output: 1
#
#
# **Example 4:**
#
#
#
#     Input: a = "abc", b = "wxyz"
#     Output: -1
#
#
#
#
# **Constraints:**
#
#   * `1 <= a.length <= 104`
#   * `1 <= b.length <= 104`
#   * `a` and `b` consist of lower-case English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
__precompile__()

module ResultTypes

using Nullables

export Result, ErrorResult, unwrap, unwrap_error, iserror

struct Result{T, E<:Exception}
    result::Nullable{T}
    error::Nullable{E}
end

Result(x::T) where {T} = Result{T, ErrorException}(Nullable{T}(x), Nullable{ErrorException}())
Result(x::T, ::Type{E}) where {T, E} = Result{T, E}(Nullable{T}(x), Nullable{E}())

function ErrorResult(::Type{T}, e::E) where {T, E<:Exception}
    Result{T, E}(Nullable{T}(), Nullable{E}(e))
end

function ErrorResult(::Type{T}, e::AbstractString="") where T
    Result{T, ErrorException}(Nullable{T}(), Nullable{ErrorException}(ErrorException(e)))
end

function unwrap(r::Result{T, E})::T where {T, E}
    if !isnull(r.result)
        return get(r.result)
    elseif !isnull(r.error)
        throw(get(r.error))
    else
        error("Empty Result{$T, $E} type")
    end
end

function unwrap_error(r::Result{T, E})::E where {T, E}
    if !isnull(r.error)
        return get(r.error)
    else
        error("$r is not an ErrorResult")
    end
end

function Base.convert(::Type{T}, r::Result{S, E})::T where {T, S, E}
    unwrap(r)
end

function Base.convert(::Type{Result{S, E}}, x::T) where {T, S, E}
    Result{S, E}(Nullable{S}(convert(S, x)), Nullable{E}())
end

function Base.convert(::Type{Result{T, E}}, e::E) where {T, E}
    Result{T, E}(Nullable{T}(), Nullable{E}(e))
end

function Base.show(io::IO, r::Result{T, E}) where {T, E}
    if iserror(r)
        print(io, "ErrorResult(", T, ", ", unwrap_error(r), ")")
    else
        print(io, "Result(", unwrap(r), ")")
    end
end

iserror(r::Result) = !isnull(r.error)

end # module


text:
using SpecialFunctions
gamma(1/2) - sqrt(pi)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module FastMarching

include("libmsfm.jl")
# include("msfm2d.jl")
include("ndgrid.jl")
include("pointmin.jl")
include("rk4.jl")
include("s1.jl")
include("shortestpath.jl")

end # module