text:
function loadstreaming(src::File)
    sfinfo = SF_INFO()
    fname = filename(src)
    # sf_open fills in sfinfo
    filePtr = sf_open(fname, SFM_READ, sfinfo)

    SndFileSource(fname, filePtr, sfinfo)
end

function loadstreaming(src::Stream, filelen=inferlen(src))
    sfinfo = SF_INFO()
    io = LengthIO(stream(src), filelen)
    # sf_open fills in sfinfo
    filePtr = sf_open(io, SFM_READ, sfinfo)

    SndFileSource(io, filePtr, sfinfo)
end

for T in (:File, :Stream), fmt in supported_formats
    @eval @inline load(src::$T{$fmt}, args...) = load_helper(src, args...)
end

# convert a `load` call into a `loadstreaming` call that properly
# cleans up the stream
function load_helper(src::Union{File, Stream}, args...)
    str = loadstreaming(src, args...)
    buf = try
        read(str)
    finally
        close(str)
    end

    buf
end

function savestreaming(src::Union{File{T}, Stream{T}}, nchannels, samplerate, elemtype) where T
    sfinfo = SF_INFO()

    sfinfo.samplerate = samplerate
    sfinfo.channels = nchannels
    sfinfo.format = formatcode(T)
    io = if src isa Stream
        LengthIO(stream(src), 0)
    else
        filename(src)
    end
    # TODO: should we auto-convert 32-bit integer samples to 24-bit?
    if T == format"FLAC" && elemtype != PCM16Sample
        error("LibSndFile.jl: FLAC only supports 16-bit integer samples")
    end
    # will probably need to figure out how this would interact with
    # ogg opus files
    if T == format"OGG"
        sfinfo.format |= SF_FORMAT_VORBIS
    else
        sfinfo.format |= subformatcode(elemtype)
    end

    filePtr = sf_open(io, SFM_WRITE, sfinfo)

    SndFileSink(io, filePtr, sfinfo)
end

for T in (:File, :Stream), fmt in supported_formats
    @eval @inline save(src::$T{$fmt}, args...; kwargs...) =
        save_helper(src, args...; kwargs...)
end

function save_helper(src, buf::SampleBuf)
    sfinfo = SF_INFO()
    stream = savestreaming(src, nchannels(buf), samplerate(buf), eltype(buf))

    try
        frameswritten = write(stream, buf)
        if frameswritten != nframes(buf)
            error("Only wrote $frameswritten frames, expected $(nframes(buf))")
        end
    finally
        # make sure we close the file even if something goes wrong
        close(stream)
    end

    nothing
end

function save_helper(src, buf::AbstractArray; samplerate=nothing)
    if samplerate === nothing
        throw(ArgumentError(string(
            "Can't save `Array` with unknown samplerate. ",
            "Specify the `samplerate` keyword argument or ",
            "use a SampleBuf")))
    end
    save_helper(src, SampleBuf(buf, samplerate))
end


text:
"""
测试能量截断的求导
"""

using Test
using MARY_fRG.Fermi
using MARY_fRG.FlowEquation

@testset "能量截断" begin
    model = common_square_lattice(0.)
    Γ4 = ECGamma4(model, 4.0, 8, 40)
    Γ4.V .+= 1.0
    println(Γ4.patches)
    bubb_pp, bubb_fs, bubb_nfs, bubb_ex, bubb_nex =
    all_bubble_ec_mt(Γ4, 1.0)
    dl = dl_ec_mt(Γ4, bubb_pp, bubb_fs, bubb_nfs, bubb_ex, bubb_nex)
    @test isapprox(dl[1, 1, 1, 1, 2, 2, 2], -0.1289112, atol=1e-6)
    @test isapprox(dl[1, 1, 1, 1, 2, 2, 3], -0.0304875, atol=1e-6)
    @test isapprox(dl[1, 1, 1, 1, 2, 2, 4], 0.0032081, atol=1e-6)
    @test isapprox(dl[1, 1, 1, 1, 2, 2, 5], 0.2625107, atol=1e-6)
    @test isapprox(dl[1, 1, 1, 1, 4, 2, 2], 0.0336956, atol=1e-6)
    @test isapprox(dl[1, 1, 1, 1, 4, 2, 3], 0.0304875, atol=1e-6)
    @test isapprox(dl[1, 1, 1, 1, 4, 2, 4], -0.098423, atol=1e-6)
    @test isapprox(dl[1, 1, 1, 1, 4, 2, 5], 0., atol=1e-6)
end

@testset

ext: jl


text:
# TODO: SInt128 promotions are type unstable!
#(::Type{T})(x::Float16) where {T<:Integer} = T(Float32(x))

#Bool(x::Real) = x==0 ? false : x==1 ? true : throw(InexactError(:Bool, Bool, x))

Base.reinterpret(::Type{<:SUnsigned}, x::SFloat16) = reinterpret(SUInt16, x)
Base.reinterpret(::Type{<:SSigned}, x::SFloat16) = reinterpret(SInt16, x)

Base.reinterpret(::Type{<:SInt16}, x::SFloat16{X}) where X = SInt16{Base.bitcast(Int16, X::Float16)}()
Base.reinterpret(::Type{<:SUInt16}, x::SFloat16{X}) where X = SUInt16{Base.bitcast(UInt16, X::Float16)}()

Base.reinterpret(::Type{<:SFloat32}, x::SUInt32{X}) where X = SFloat32{Base.bitcast(Float32, X::UInt32)}()
Base.reinterpret(::Type{<:SFloat32}, x::SInt32{X}) where X = SFloat32{Base.bitcast(Float32, X::Int32)}()

#=

function _sfloat32(exp::SUInt32ZeroType, sig::SUInt32, sign::SUInt32)
    __sfloat32(SOne, SUInt16(0x0200), sig, sign)
end
function __sfloat32(n_bit::SInt, bit::SUInt16, sig::SUInt32ZeroType, sign::SUInt32)
    if (bit & sig) == 0
        __sfloat32(n_bit, bit, sig, sign)
    else
        return (sign << SInt(31)) | (((-SInt(14) - n_bit + SInt(127)) << SInt(23)) % SUInt32) | (((sig & (~bit)) << n_bit) << (SInt(23) - SInt(10)))
    end
end

#=
ival = reinterpret(UInt16, Float16(3))
(ival & 0x7c00) >> 10

ival = reinterpret(SUInt16, SFloat16(1))
exp = (ival & SUInt16{0x7c00}()) >> SInt{10}()
sig = (ival & SUInt16{0x3ff}()) >> SZero
sign = (ival & SUInt16{0x8000}()) >> SInt{15}()

(ival & SUInt16(0x7c00)) >> SInt(10)
(ival & 0x7c00) >> 10
=#

_sfloat32(exp::SUInt32ZeroType,     sig::SUInt32ZeroType, sign::SUInt32) = (sign << SInt(31)) | exp | sig
_sfloat32(exp::SUInt32{0x0000001f}, sig::SUInt32ZeroType, sign::SUInt32ZeroType) = SUInt32(0x7f800000)
_sfloat32(exp::SUInt32{0x0000001f}, sig::SUInt32ZeroType, sign::SUInt32) = 0xff800000
_sfloat32(exp::SUInt32{0x0000001f}, sig::SUInt32,         sign::SUInt32) = SUInt32(0x7fc00000) | (sign<<SInt(31)) | (sig<<(SInt(13)))
_sfloat32(exp::SUInt32,             sig::SUInt32,         sign::SUInt32) = (sign << SInt(31)) | (((exp + SInt(112)) << SInt(23)) % SUInt32) | (sig << SInt(13))

function SFloat32(val::SFloat16)
    local ival::SUInt32 = reinterpret(SUInt16, val)
    return reinterpret(SFloat32, _sfloat32((ival & SUInt16{0x7c00}()) >> SInt{10}(),
                                           (ival & SUInt16{0x3ff}()) >> SZero,
                                           (ival & SUInt16{0x8000}()) >> SInt{15}()))
end
=#

Base.sign_mask(::Type{<:SFloat64})        = SUInt64(0x8000_0000_0000_0000)
Base.exponent_mask(::Type{<:SFloat64})    = SUInt64(0x7ff0_0000_0000_0000)
Base.exponent_one(::Type{<:SFloat64})     = SUInt64(0x3ff0_0000_0000_0000)
Base.exponent_half(::Type{<:SFloat64})    = SUInt64(0x3fe0_0000_0000_0000)
Base.significand_mask(::Type{<:SFloat64}) = SUInt64(0x000f_ffff_ffff_ffff)

Base.sign_mask(::Type{<:SFloat32})        = SUInt32(0x8000_0000)
Base.exponent_mask(::Type{<:SFloat32})    = SUInt32(0x7f80_0000)
Base.exponent_one(::Type{<:SFloat32})     = SUInt32(0x3f80_0000)
Base.exponent_half(::Type{<:SFloat32})    = SUInt32(0x3f00_0000)
Base.significand_mask(::Type{<:SFloat32}) = SUInt32(0x007f_ffff)

Base.sign_mask(::Type{<:SFloat16}) =        SUInt16(0x8000)
Base.exponent_mask(::Type{<:SFloat16}) =    SUInt16(0x7c00)
Base.exponent_one(::Type{<:SFloat16}) =     SUInt16(0x3c00)
Base.exponent_half(::Type{<:SFloat16}) =    SUInt16(0x3800)
Base.significand_mask(::Type{<:SFloat16}) = SUInt16(0x03ff)

# integer size of float
Base.uinttype(::Type{<:SFloat64}) = SUInt64
Base.uinttype(::Type{<:SFloat32}) = SUInt32
Base.uinttype(::Type{<:SFloat16}) = SUInt16

const IEEESFloat = Union{SFloat16,SFloat32,SFloat64}

significand_bits(::Type{T}) where T<:IEEESFloat = SInt(trailing_ones(Base.significand_mask(T)))
exponent_bits(::Type{T}) where T<:IEEESFloat = SInt(sizeof(eltype(T))*8 - significand_bits(SFloat64) - SOne)
exponent_bias(::Type{T}) where T<:IEEESFloat = SInt(Base.exponent_one(T) >> significand_bits(T))
exponent_max(::Type{T}) where T<:IEEESFloat = SInt(Base.exponent_mask(T) >> significand_bits(T)) - exponent_bias(T)
exponent_raw_max(::Type{T}) where T<:IEEESFloat = SInt(Base.exponent_mask(T) >> significand_bits(T))

#Base.:(-)(x::SFloat32{V1}, y::SFloat16{V2}) where {V1,V2} = SFloat32{(-)(V1::Float32, V2::Float16)::Float32}()

for (ST,BT) in SF2BF
    @eval begin
        Base.ldexp(::$ST{X}, n::SInteger) where X = $ST{ldexp(X::$BT, values(n))}()
    end
end

for (ST,BT) in SF2BF
    for (ST2,BT2) in SI2BI
        @eval begin
            Base.trunc(::Type{$BT2}, x::$ST{X}) where {T,X} = $ST2{trunc($BT2, X::$BT)}()
        end
    end
end

Base.unsigned(x::SFloat) = Base.unsigned(SInt(x))

Base.maxintfloat(x::T) where {T<:SFloat} = maxintfloat(T)

Base.widen(::Type{SFloat16}) = SFloat32
Base.widen(::Type{SFloat32}) = SFloat64

Base.AbstractFloat(x::SFloat) = x


text:
#=
    Vector valued constraints
=#
constrdict(model::LinQuadOptimizer, ::VLCI{MOI.Nonnegatives}) = cmap(model).nonnegatives
constrdict(model::LinQuadOptimizer, ::VLCI{MOI.Nonpositives}) = cmap(model).nonpositives
constrdict(model::LinQuadOptimizer, ::VLCI{MOI.Zeros})        = cmap(model).zeros

function MOI.add_constraint(model::LinQuadOptimizer, func::VecLin, set::S) where S <: VecLinSets
    __assert_supported_constraint__(model, VecLin, S)
    @assert MOI.dimension(set) == length(func.constants)
    nrows = get_number_linear_constraints(model)
    add_linear_constraint(model, func, backend_type(model, set))
    nrows2 = get_number_linear_constraints(model)
    model.last_constraint_reference += 1
    index = VLCI{S}(model.last_constraint_reference)
    dict = constrdict(model, index)
    dict[index] = collect(nrows+1:nrows2)
    for i in 1:MOI.dimension(set)
        push!(model.constraint_primal_solution, NaN)
        push!(model.constraint_dual_solution, NaN)
        push!(model.constraint_constant, func.constants[i])
    end
    return index
end

function add_linear_constraint(model::LinQuadOptimizer, func::VecLin, sense::Cchar)

    # sort terms by output index, then by variable index
    func = copy(func)
    MOIU.sort_and_compress!(func.terms, MOI.term_indices, t -> true, MOIU.unsafe_add)
    columns       = [get_column(model, term.scalar_term.variable_index) for term in func.terms]
    coefficients  = [term.scalar_term.coefficient for term in func.terms]

    # Compute the row pointers in the compressed sparse row matrix. The row
    # pointers are defined recursively:
    #  r[1] = 1
    #  r[i] = r[i - 1] + (number of nonzero elements in the (i - 1)th row)
    #
    # To compute this, we first count up the number of nonzero elements in each
    # row (i - 1), storing the result in r[i]. Then we perform a cumsum on r,
    # storing the result back in r.
    num_rows = length(func.constants)
    row_pointers = fill(0, num_rows+1)
    row_pointers[1] = 1
    for term in func.terms
        row = term.output_index
        if row < 1 || row > num_rows
            throw(ArgumentError("Output index $row out of range [1, $num_rows]"))
        end
        if row < num_rows
            row_pointers[row + 1] += 1
        end
    end
    cumsum!(row_pointers, row_pointers)
    row_pointers[end]=length(coefficients)+1
    A = CSRMatrix{Float64}(row_pointers, columns, coefficients)
    add_linear_constraints!(model, A, fill(sense, num_rows), -func.constants)
end

function MOI.modify(model::LinQuadOptimizer, index::VLCI{<: VecLinSets},
                     change::MOI.VectorConstantChange{Float64})
    rows = model[index]
    @assert length(change.new_constant) == length(rows)
    for (row, constant) in zip(rows, change.new_constant)
        change_rhs_coefficient!(model, row, -constant)
        model.constraint_constant[row] = constant
    end
end

function MOI.modify(model::LinQuadOptimizer, index::VLCI{<: VecLinSets},
                     change::MOI.MultirowChange{Float64})
    column = get_column(model, change.variable)
    for (row, coef) in change.new_coefficients
        change_matrix_coefficient!(model, row, column, coef)
    end
end

_constant(f::Linear) = f.constant
_constant(f::VecLin) = f.constants

"""
    _matching_sparsity_pattern(f1::T, f2::T) where {T <: Union{Linear, VecLin}}

Internal function, not intended for external use.

Determines whether functions `f1` and `f2` have the same sparsity pattern
w.r.t. their constraint rows and columns. Assumes both functions are already
in canonical form.
"""
function _matching_sparsity_pattern(f1::T, f2::T) where {T <: Union{Linear, VecLin}}
    axes(f1.terms) == axes(f2.terms) || return false
    axes(_constant(f1)) == axes(_constant(f2)) || return false
    for i in eachindex(f1.terms)
        MOI.term_indices(f1.terms[i]) == MOI.term_indices(f2.terms[i]) || return false
    end
    return true
end

"""
    _replace_with_matching_sparsity!(model::LinQuadOptimizer,
        previous::VecLin, replacement::VecLin, row)

Internal function, not intended for external use.

Change the linear constraint function at index `row` in `model` from
`previous` to `replacement`. This function assumes that `previous` and
`replacement` have exactly the same sparsity pattern w.r.t. which variables
they include and that both constraint functions are in canonical form (as
returned by `MOIU.canonical()`. Neither assumption is checked within the body
of this function.
"""
function _replace_with_matching_sparsity!(model::LinQuadOptimizer, previous::VecLin, replacement::VecLin, constraint_indices)
    rows = [constraint_indices[t.output_index] for t in previous.terms]
    cols = [model.variable_mapping[t.scalar_term.variable_index] for t in previous.terms]
    coefs = MOI.coefficient.(replacement.terms)
    change_matrix_coefficients!(model, rows, cols, coefs)
end

"""
    _replace_with_different_sparsity!(model::LinQuadOptimizer,
        previous::VecLin, replacement::VecLin, row)

Internal function, not intended for external use.

Change the linear constraint function at index `row` in `model` from
`previous` to `replacement`. This function assumes that `previous` and
`replacement` may have different sparsity patterns.

This function (and `_replace_with_matching_sparsity!` above) are necessary
because the LQOI API currently *only* allows linear constraint modification
through the `change_matrix_coefficient!` and `change_matrix_coefficients!`
functions. In order to fully replace a linear constraint, we have to zero out
the current matrix coefficients and then set the new matrix coefficients. When
the sparsity patterns match, the zeroing-out step can be skipped.
"""
function _replace_with_different_sparsity!(model::LinQuadOptimizer, previous::VecLin, replacement::VecLin, constraint_indices)
    # First, zero out the old constraint function terms
    rows = [constraint_indices[t.output_index] for t in previous.terms]
    cols = [model.variable_mapping[t.scalar_term.variable_index] for t in previous.terms]
    coefs = fill(0.0, length(previous.terms))
    change_matrix_coefficients!(model, rows, cols, coefs)

    # Next, set the new constraint function terms
    rows = [constraint_indices[t.output_index] for t in replacement.terms]
    cols = [model.variable_mapping[t.scalar_term.variable_index] for t in replacement.terms]
    coefs = MOI.coefficient.(replacement.terms)
    change_matrix_coefficients!(model, rows, cols, coefs)
end

MOI.supports(::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{VLCI{S}}) where {S <: VecLinSets} = true
function MOI.set(model::LinQuadOptimizer, attr::MOI.ConstraintFunction, c::VLCI{S}, replacement::VecLin) where {S <: VecLinSets}
    replacement = MOI.Utilities.canonical(replacement)
    previous = MOI.get(model, attr, c)
    MOI.Utilities.canonicalize!(previous)
    # If the previous and replacement constraint functions have exactly
    # the same sparsity pattern, then we can take a faster path by just
    # passing the replacement terms to the model. But if their sparsity
    # patterns differ, then we need to first zero out the previous terms
    # and then set the replacement terms.
    constraint_indices = model[c]
    if _matching_sparsity_pattern(previous, replacement)
        _replace_with_matching_sparsity!(model, previous, replacement, constraint_indices)
    else
        _replace_with_different_sparsity!(model, previous, replacement, constraint_indices)
    end
    for i in eachindex(constraint_indices)
        row = constraint_indices[i]
        model.constraint_constant[row] = replacement.constants[i]
        change_rhs_coefficient!(model, row, -replacement.constants[i])
    end
end

function MOI.delete(model::LinQuadOptimizer, index::VLCI{<:VecLinSets})
    __assert_valid__(model, index)
    delete_constraint_name(model, index)
    dict = constrdict(model, index)
    # we delete rows from largest to smallest here so that we don't have
    # to worry about updating references in a greater numbered row, only to
    # modify it later.
    for row in sort(dict[index], rev=true)
        delete_linear_constraints!(model, row, row)
        deleteat!(model.constraint_primal_solution, row)
        deleteat!(model.constraint_dual_solution, row)
        deleteat!(model.constraint_constant, row)
        # shift all the other references
        shift_references_after_delete_affine!(model, row)
    end
    delete!(dict, index)
end

#=
    Constraint set of Linear function
=#

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::VLCI{S}) where S <: VecLinSets
    constraint_indices = model[index]
    S(length(constraint_indices))
end

#=
    Constraint function of Linear function
=#

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintFunction, index::VLCI{<: VecLinSets})
    rows = model[index]
    constants = Float64[]
    terms = MOI.VectorAffineTerm{Float64}[]
    for (i, row) in enumerate(rows)
        rhs = get_rhs(model, row)
        push!(constants, -rhs)
        columns, coefficients = get_linear_constraint(model, row)
        for (column, coefficient) in zip(columns, coefficients)
            push!(terms, MOI.VectorAffineTerm{Float64}(
                    i,
                    MOI.ScalarAffineTerm{Float64}(
                        coefficient, model.variable_references[column]
                    )
                )
            )
        end
    end
    return MOI.VectorAffineFunction(terms, constants)
end


text:
using .Turing, Distributions, Test, Random
using DynamicPPL: VarInfo

Random.seed!(129)

dir = splitdir(splitdir(pathof(DynamicPPL))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

Random.seed!(129)

@testset "prob_macro" begin
    @testset "scalar" begin
        @model demo(x) = begin
            m ~ Normal()
            x ~ Normal(m, 1)
        end

        mval = 3
        xval = 2
        iters = 1000

        logprior = logpdf(Normal(), mval)
        loglike = logpdf(Normal(mval, 1), xval)
        logjoint = logprior + loglike

        @test logprob"m = mval | model = demo" == logprior
        @test logprob"m = mval | x = xval, model = demo" == logprior
        @test logprob"x = xval | m = mval, model = demo" == loglike
        @test logprob"x = xval, m = mval | model = demo" == logjoint

        varinfo = VarInfo(demo(xval))
        @test logprob"m = mval | model = demo, varinfo = varinfo" == logprior
        @test logprob"m = mval | x = xval, model = demo, varinfo = varinfo" == logprior
        @test logprob"x = xval | m = mval, model = demo, varinfo = varinfo" == loglike
        varinfo = VarInfo(demo(missing))
        @test logprob"x = xval, m = mval | model = demo, varinfo = varinfo" == logjoint

        chain = sample(demo(xval), IS(), iters; save_state = true)
        chain2 = Chains(chain.value, chain.logevidence, chain.name_map, NamedTuple())
        lps = logpdf.(Normal.(vec(chain["m"].value), 1), xval)
        @test logprob"x = xval | chain = chain" == lps
        @test logprob"x = xval | chain = chain2, model = demo" == lps
        varinfo = VarInfo(demo(xval))
        @test logprob"x = xval | chain = chain, varinfo = varinfo" == lps
        @test logprob"x = xval | chain = chain2, model = demo, varinfo = varinfo" == lps
    end

    @testset "vector" begin
        n = 5
        @model demo(x, n = n, ::Type{T} = Float64) where {T} = begin
            m = Vector{T}(undef, n)
            @. m ~ Normal()
            @. x ~ Normal.(m, 1)
        end
        mval = rand(n)
        xval = rand(n)
        iters = 1000

        logprior = sum(logpdf.(Normal(), mval))
        like(m, x) = sum(logpdf.(Normal.(m, 1), x))
        loglike = like(mval, xval)
        logjoint = logprior + loglike

        @test logprob"m = mval | model = demo" == logprior
        @test logprob"x = xval | m = mval, model = demo" == loglike
        @test logprob"x = xval, m = mval | model = demo" == logjoint

        varinfo = VarInfo(demo(xval))
        @test logprob"m = mval | model = demo, varinfo = varinfo" == logprior
        @test logprob"x = xval | m = mval, model = demo, varinfo = varinfo" == loglike
        # Currently, we cannot easily pre-allocate `VarInfo` for vector data

        chain = sample(demo(xval), HMC(0.5, 1), iters; save_state = true)
        chain2 = Chains(chain.value, chain.logevidence, chain.name_map, NamedTuple())
        lps = like.([[chain["m[$i]"].value[j] for i in 1:n] for j in 1:iters], Ref(xval))
        @test logprob"x = xval | chain = chain" == lps
        @test logprob"x = xval | chain = chain2, model = demo" == lps
        @test logprob"x = xval | chain = chain, varinfo = varinfo" == lps
        @test logprob"x = xval | chain = chain2, model = demo, varinfo = varinfo" == lps
    end
end


text:
"""
Smetswoutersorig{T} <: AbstractDSGEModel{T}

The `SmetsWoutersOrig` type implements the *original* model from \"Shocks and
Frictions in US Business Cycles: A Bayesian DSGE Approach\" (2007) by Frank
Smets and Rafael Wouters.

In this version of the model (`SmetsWoutersOrig`), variables are detrended using
the deterministic trend ℯ^{γt}. Contrast with the `SmetsWouters` version, which
detrends variables using Z_t ℯ^{(γ + α/(1-α) log Υ)t}.

### Fields

#### Parameters and Steady-States

* `parameters::Vector{AbstractParameter}`: Vector of all time-invariant model
  parameters.

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

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

#### Inputs to Measurement and Equilibrium Condition Equations

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

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

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

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

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

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

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

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

#### Model Specifications and Settings

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

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

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

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

#### Other Fields

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

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

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

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

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

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

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

description(m::SmetsWoutersOrig) = "Smets-Wouters Original Model"

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

Initializes indices for all of `m`'s states, shocks, and equilibrium conditions.

### Arguments

`m:: SmetsWoutersOrig`: a model object
"""
function init_model_indices!(m::SmetsWoutersOrig)
    # Endogenous states
    endogenous_states = [[
        :g_t, :b_t, :μ_t, :z_t, :λ_f_t, :λ_w_t, :rm_t,           # exogenous processes
        :y_t, :c_t, :i_t, :qk_t, :k_t, :kbar_t, :u_t, :rk_t,     # production and capital
        :mc_t, :π_t, :μ_ω_t, :w_t, :L_t, :R_t, :λ_f_t1, :λ_w_t1, # prices, wages, interest rate
        :Ec_t, :Eqk_t, :Ei_t, :Eπ_t, :EL_t, :Erk_t, :Ew_t,       # expectations
        :y_f_t, :c_f_t, :i_f_t, :qk_f_t, :k_f_t, :kbar_f_t, :u_f_t, :rk_f_t, # flex-price production/capital
        :w_f_t, :L_f_t, :r_f_t,                                  # flex-price prices, wages, interest rate
        :Ec_f_t, :Eqk_f_t, :Ei_f_t, :EL_f_t, :Erk_f_t];          # flex-price expectations
        [Symbol("rm_tl$i") for i = 1:n_anticipated_shocks(m)]]

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

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

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

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

    # Observables
    observables = keys(m.observable_mappings)

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

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

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

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

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

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

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

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

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

    # Initialize parameters
    init_parameters!(m)

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

    return m
end

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

Initializes the model's parameters, as well as empty values for the steady-state
parameters (in preparation for `steadystate!(m)` being called to initialize
those).
"""
function init_parameters!(m::SmetsWoutersOrig)
    m <= parameter(:S′′, 6.3325, (-15., 15.), (-15., 15.), ModelConstructors.Untransformed(), Distributions.Normal(4., 1.5), fixed=false,
                   description="S'': The second derivative of households' cost of adjusting investment.",
                   tex_label="S''")

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    m <= parameter(:Lmean, 1.2918, (-5., 5.), (-5., 5.), ModelConstructors.Untransformed(), Distributions.Normal(0.0, 2.0), fixed=false,
                   description="Lmean: Mean level of hours.",
                   tex_label="\\bar{L}")

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

    m <= parameter(:α, 0.24, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Distributions.Normal(0.30, 0.05), fixed=false,
                   description="α: Capital elasticity in the intermediate goods sector's Cobb-Douglas production function.",
                   tex_label="\\alpha")

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

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

    m <= parameter(:ρ_b, 0.2703, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.2), fixed=false,
                   description="ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label="\\rho_b")
#added 4 9's to upper bound for the following parameter
    m <= parameter(:ρ_g, 0.9930, (1e-5, 0.9999999), (1e-5, 0.9999999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.2), fixed=false,
                   description="ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label="\\rho_g")

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    # steady states
    m <= SteadyStateParameter(:Rstarn, NaN, tex_label="\\R_*_n")
    m <= SteadyStateParameter(:rkstar, NaN, tex_label="\\r^k_*")
    m <= SteadyStateParameter(:wstar, NaN, tex_label="\\w_*")
    m <= SteadyStateParameter(:i_k, NaN, tex_label="")
    m <= SteadyStateParameter(:l_k, NaN, tex_label="")
    m <= SteadyStateParameter(:k_y, NaN, tex_label="")
    m <= SteadyStateParameter(:i_y, NaN, tex_label="\\frac{i_*}{y_*}")
    m <= SteadyStateParameter(:c_y, NaN, tex_label="\\frac{c_*}{y_*}")
    m <= SteadyStateParameter(:u_y, NaN, tex_label="\\frac{r^k_* k_*}{y_*}")
    m <= SteadyStateParameter(:wl_c, NaN, tex_label="wl_c")
end

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

Calculates the model's steady-state values. `steadystate!(m)` must be called whenever the parameters of `m` are updated.
"""
function steadystate!(m::SmetsWoutersOrig)

    #trying to match with fortran
    #m[:γ] = m[:γ]/100 + 1
    #   m[:β] = 1/(m[:β]/100+1)
    #   m[:π_star] = m[:π_star]/100+1

    m[:Rstarn]   = 100 * (m[:γ]^m[:σ_c] * m[:π_star] / m[:β] - 1)
    m[:rkstar]   = m[:γ]^m[:σ_c] / m[:β] - (1 - m[:δ])
    m[:wstar]    = (m[:α]^m[:α] * (1-m[:α])^(1-m[:α]) / (m[:Φ] * m[:rkstar]^m[:α]))^(1/(1-m[:α]))
    m[:i_k]      = (1 - (1-m[:δ])/m[:γ]) * m[:γ]
    m[:l_k]      = (1 - m[:α])/m[:α] * m[:rkstar]/m[:wstar]
    m[:k_y]      = m[:Φ] * m[:l_k]^(m[:α] - 1)
    m[:i_y]      = (m[:γ] - 1 +m[:δ]) * m[:k_y]
    m[:c_y]      = 1 - m[:g_star] - m[:i_y]
    m[:u_y]      = m[:rkstar] * m[:k_y]
    m[:wl_c]     = 1 / m[:λ_w] * (1-m[:α]) / m[:α] * m[:rkstar] * m[:k_y] / m[:c_y]

    return m
end

function model_settings!(m::SmetsWoutersOrig)

    default_settings!(m)

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

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

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


text:
# Test struct for `rand_tangent` and `difference`.
struct Foo
    a::Float64
    b::Int
    c::Any
 end

# to avoid deprecation spam (and actually test deprecations) we will define a wrapper `rand_tangent` function for testing
rand_tangent(args...) = @test_deprecated FiniteDifferences.rand_tangent(args...)

@testset "rand_tangent" begin
    rng = MersenneTwister(123456)

    @testset "Primal: $(typeof(x)), Tangent: $T_tangent" for (x, T_tangent) in [

        # Things without sensible tangents.
        ("hi", NoTangent),
        ('a', NoTangent),
        (:a, NoTangent),
        (true, NoTangent),
        (4, NoTangent),
        (FiniteDifferences, NoTangent),  # Module object
        # Types (not instances of type)
        (Foo, NoTangent),
        (Union{Int, Foo}, NoTangent),
        (Union{Int, Foo}, NoTangent),
        (Vector, NoTangent),
        (Vector{Float64}, NoTangent),
        (Integer, NoTangent),
        (Type{<:Real}, NoTangent),

        # Numbers.
        (5.0, Float64),
        (5.0 + 0.4im, Complex{Float64}),
        (big(5.0), BigFloat),

        # StridedArrays.
        (fill(randn(Float32)), Array{Float32, 0}),
        (fill(randn(Float64)), Array{Float64, 0}),
        (randn(Float32, 3), Vector{Float32}),
        (randn(Complex{Float64}, 2), Vector{Complex{Float64}}),
        (randn(5, 4), Matrix{Float64}),
        (randn(Complex{Float32}, 5, 4), Matrix{Complex{Float32}}),
        ([randn(5, 4), 4.0], Vector{Any}),

        # Wrapper Arrays
        (randn(5, 4)', Adjoint{Float64, Matrix{Float64}}),
        (transpose(randn(5, 4)), Transpose{Float64, Matrix{Float64}}),

        # Tuples.
        ((4.0, ), Tangent{Tuple{Float64}}),
        ((5.0, randn(3)), Tangent{Tuple{Float64, Vector{Float64}}}),

        # NamedTuples.
        ((a=4.0, ), Tangent{NamedTuple{(:a,), Tuple{Float64}}}),
        ((a=5.0, b=1), Tangent{NamedTuple{(:a, :b), Tuple{Float64, Int}}}),

        # structs.
        (Foo(5.0, 4, rand(rng, 3)), Tangent{Foo}),
        (Foo(4.0, 3, Foo(5.0, 2, 4)), Tangent{Foo}),
        (sin, NoTangent),
        # all fields NoTangent implies NoTangent
        (Pair(:a, "b"), NoTangent),
        (1:10, NoTangent),
        (1:2:10, NoTangent),

        # LinearAlgebra types (also just structs).
        (
            UpperTriangular(randn(3, 3)),
            Tangent{UpperTriangular{Float64, Matrix{Float64}}},
        ),
        (
            Diagonal(randn(2)),
            Tangent{Diagonal{Float64, Vector{Float64}}},
        ),
        (
            SVector{2, Float64}(1.0, 2.0),
            Tangent{typeof(SVector{2, Float64}(1.0, 2.0))},
        ),
        (
            SMatrix{2, 2, ComplexF64}(1.0, 2.0, 3.0, 4.0),
            Tangent{typeof(SMatrix{2, 2, ComplexF64}(1.0, 2.0, 3.0, 4.0))},
        ),
        (
            Symmetric(randn(2, 2)),
            Tangent{Symmetric{Float64, Matrix{Float64}}},
        ),
        (
            Hermitian(randn(ComplexF64, 1, 1)),
            Tangent{Hermitian{ComplexF64, Matrix{ComplexF64}}},
        ),
    ]
        @test rand_tangent(rng, x) isa T_tangent
        @test rand_tangent(x) isa T_tangent
    end

    @testset "erroring cases" begin
        # Ensure struct fallback errors for non-struct types.
        @test_throws ArgumentError invoke(
            FiniteDifferences.rand_tangent, Tuple{AbstractRNG, Any}, rng, 5.0
        )
    end

    @testset "compsition of addition" begin
        x = Foo(1.5, 2, Foo(1.1, 3, [1.7, 1.4, 0.9]))
        @test x + rand_tangent(x) isa typeof(x)
        @test x + (rand_tangent(x) + rand_tangent(x)) isa typeof(x)
    end

    # Julia 1.6 changed to using Ryu printing algorithm and seems better at printing short
    VERSION > v"1.6" && @testset "niceness of printing" begin
        for i in 1:50
            @test length(string(rand_tangent(1.0))) <= 6
            @test length(string(rand_tangent(1.0 + 1.0im))) <= 12
            @test length(string(rand_tangent(1f0))) <= 12
            @test length(string(rand_tangent(big"1.0"))) <= 12
        end
    end
end


text:
@testset "BFARadPowerModel" begin

    @testset "CIGRE TNEP single-period" begin
        data = _FP.parse_file(normpath(@__DIR__,"../test/data/cigre_mv_eu/cigre_mv_eu_unit_test.m"))
        _FP.add_dimension!(data, :hour, 1)
        _FP.add_dimension!(data, :year, 1)
        data = _FP.make_multinetwork(data)
        result = _FP.flex_tnep(data, _FP.BFARadPowerModel, cbc)
        sol = result["solution"]["nw"]["1"]
        @test result["termination_status"] == _PM.OPTIMAL
        @test result["objective"] ≈ 4360.45 rtol=1e-3
        @test sol["branch"]["16"]["pf"] ≈ -sol["branch"]["16"]["pt"] rtol=1e-3 # Zero active power losses in OLTC branch
        @test sol["branch"]["16"]["qf"] ≈ -sol["branch"]["16"]["qt"] rtol=1e-3 # Zero reactive power losses in OLTC branch
        @test sol["branch"]["17"]["pf"] ≈ -sol["branch"]["17"]["pt"] rtol=1e-3 # Zero active power losses in frb branch
        @test sol["branch"]["17"]["qf"] ≈ -sol["branch"]["17"]["qt"] rtol=1e-3 # Zero reactive power losses in frb branch
        @test sol["branch"]["1"]["pf"] ≈ -sol["branch"]["1"]["pt"] rtol=1e-3 # Zero active power losses in regular branch
        @test sol["branch"]["1"]["qf"] ≈ -sol["branch"]["1"]["qt"] rtol=1e-3 # Zero reactive power losses in regular branch
        @test sol["ne_branch"]["1"]["built"] ≈ 0.0 atol=1e-1 # Unused OLTC ne_branch
        @test sol["ne_branch"]["2"]["built"] ≈ 0.0 atol=1e-1 # Unused frb ne_branch
        @test sol["ne_branch"]["3"]["built"] ≈ 0.0 atol=1e-1 # Unused regular ne_branch
        @test sol["ne_branch"]["1"]["pf"] ≈ 0.0 atol=1e-2 # Zero active power in unused OLTC ne_branch
        @test sol["ne_branch"]["1"]["qf"] ≈ 0.0 atol=1e-2 # Zero reactive power in unused OLTC ne_branch
        @test sol["ne_branch"]["2"]["pf"] ≈ 0.0 atol=1e-2 # Zero active power in unused frb ne_branch
        @test sol["ne_branch"]["2"]["qf"] ≈ 0.0 atol=1e-2 # Zero reactive power in unused frb ne_branch
        @test sol["ne_branch"]["3"]["pf"] ≈ 0.0 atol=1e-2 # Zero active power in unused regular ne_branch
        @test sol["ne_branch"]["3"]["qf"] ≈ 0.0 atol=1e-2 # Zero reactive power in unused regular ne_branch
        @test sum(g["pg"] for g in values(sol["gen"])) ≈ sum(l["pflex"] for l in values(sol["load"])) rtol=1e-3 # Zero overall active power losses
        @test sum(g["qg"] for g in values(sol["gen"])) ≈ sum(l["qflex"] for l in values(sol["load"])) rtol=1e-3 # Zero overall reactive power losses

        data["nw"]["1"]["load"]["1"]["pd"] += 10.0 # Bus 1. Changes reactive power demand too, via `pf_angle`.
        data["nw"]["1"]["load"]["12"]["pd"] += 4.0 # Bus 13. Changes reactive power demand too, via `pf_angle`.
        data["nw"]["1"]["branch"]["12"]["rate_a"] = data["nw"]["1"]["branch"]["12"]["rate_b"] = data["nw"]["1"]["branch"]["12"]["rate_c"] = 0.0
        result = _FP.flex_tnep(data, _FP.BFARadPowerModel, cbc)
        sol = result["solution"]["nw"]["1"]
        @test result["termination_status"] == _PM.OPTIMAL
        @test result["objective"] ≈ 5764.48 rtol=1e-3
        @test sol["ne_branch"]["1"]["built"] ≈ 1.0 atol=1e-1 # Replacement OLTC ne_branch
        @test sol["ne_branch"]["2"]["built"] ≈ 1.0 atol=1e-1 # frb ne_branch added in parallel
        @test sol["ne_branch"]["3"]["built"] ≈ 1.0 atol=1e-1 # Replacement regular ne_branch
        @test sol["ne_branch"]["4"]["built"] ≈ 0.0 atol=1e-1 # Unused ne_branch
        @test sol["ne_branch"]["1"]["pf"] ≈ -sol["ne_branch"]["1"]["pt"] rtol=1e-3 # Zero active power losses in OLTC ne_branch
        @test sol["ne_branch"]["1"]["qf"] ≈ -sol["ne_branch"]["1"]["qt"] rtol=1e-3 # Zero reactive power losses in OLTC ne_branch
        @test sol["ne_branch"]["2"]["pf"] ≈ -sol["ne_branch"]["2"]["pt"] rtol=1e-3 # Zero active power losses in frb ne_branch
        @test sol["ne_branch"]["2"]["qf"] ≈ -sol["ne_branch"]["2"]["qt"] rtol=1e-3 # Zero reactive power losses in frb ne_branch
        @test sol["ne_branch"]["3"]["pf"] ≈ -sol["ne_branch"]["3"]["pt"] rtol=1e-3 # Zero active power losses in regular ne_branch
        @test sol["ne_branch"]["3"]["qf"] ≈ -sol["ne_branch"]["3"]["qt"] rtol=1e-3 # Zero reactive power losses in regular ne_branch
        @test sol["branch"]["16"]["pf"] ≈ 0.0 atol=1e-2 # Zero active power in replaced OLTC branch
        @test sol["branch"]["16"]["qf"] ≈ 0.0 atol=1e-2 # Zero reactive power in replaced OLTC branch
        @test sol["branch"]["17"]["pf"] ≈ 0.0 atol=1e-2 # Zero active power in replaced frb branch
        @test sol["branch"]["17"]["qf"] ≈ 0.0 atol=1e-2 # Zero reactive power in replaced frb branch
        @test sol["branch"]["13"]["pf"] ≈ 0.0 atol=1e-2 # Zero active power in replaced regular branch
        @test sol["branch"]["13"]["qf"] ≈ 0.0 atol=1e-2 # Zero reactive power in replaced regular branch
        @test sol["branch"]["12"]["pf"] ≈ 0.0 atol=1e-2 # Zero active power in branch having zero rating
        @test sol["branch"]["12"]["qf"] ≈ 0.0 atol=1e-2 # Zero reactive power in branch having zero rating
    end

end;


text:
module ScatteredArrays

using VectorizationBase, SIMDPirates, Base.Cartesian, PaddedMatrices
import PaddedMatrices: type_length

export ScatteredArray, ScatteredVector, ScatteredMatrix, LinearStorage, ChunkedArray

const MultiDimIndex{N} = Union{Tuple{Vararg{<:Integer, N}}, NTuple{N,<:Integer}, CartesianIndex{N}}

include("utilities.jl")
include("types.jl")
include("scattered_struct_array.jl")
include("chunked_struct_array.jl")

end # module


text:
# PROX_MATRIX    The proximal operator of a matrix function.
#
#   Suppose F is a orthogonally invariant matrix function such that
#   F(X) = f(s(X)), where s is the singular value map and f is some
#   absolutely symmetric function. Then
#
#     X = prox_matrix(V,lambda,prox_f)
#
#   evaluates the proximal operator of F via the proximal operator
#   of f. Here, it must be possible to evaluate prox_f as prox_f(v,lambda).
#
#   For example,
#
#     prox_matrix(V,lambda,prox_l1)
#
#   evaluates the proximal operator of the nuclear norm at V
#   (i.e., the singular value thresholding operator).
function prox_matrix(v, lambda, prox_f)
    U, S, V] = svd(v, 'econ')
    x = U * diag(prox_f(diag(S), lambda)) * V'
	return x
end


text:
# Load modules
using Plots
using NCDatasets
using Statistics: mean

# Setup run-time enviromnent
ENV["GKSwstype"] = "100"

# Specify needed directories and filenames: these will be replaced automatically by the bash run script or can be modified manually
CLIMA_NETCDF = "../netcdf/";

# Get the current and previous  GCM netcdf file names in the CLIMA_NETCDF directory
fnames = filter(x -> occursin(".nc", x), readdir( CLIMA_NETCDF ) );
#print(fnames)

filename = "$CLIMA_NETCDF"fnames[end] # selects the last file on list

# extract data
ds = NCDataset(filename, "r");
print(ds)
x = ds["x"][:]
y = ds["y"][:]
z = ds["z"][:]
time = ds["time"][:]

ρqliq = ds["moisture.q_liq"][:] # e.g. select q_liq to visualise clouds

function get_min_max(var);
    vmax = maximum(filter(!isnan,var))
    vmin = minimum(filter(!isnan,var))
    return vmin,vmax
    none
end

# Animation
clims = ( get_min_max(ρqliq) )
diag_dt_days =  (time[2] - time[1]).value / (1000*60*60*24) # get simtime

anim = @animate for t_i in 1:length(time)
    one_plot = contourf( x, z, (ρqliq[:,y_i,:,t_i] )',  xlabel="x (m)", ylabel="z (m)", clims = clims)
    plot(one_plot, layout=(1,1), size=(400, 400) )
end
mp4(anim, "plot_y_slice_anim.gif", fps = 3) # NB: mp4 can be more flaky than gif


text:
module ComposeBackend
using GLAbstraction, GLFW, GLWindow, Measures, GLVisualize, GeometryTypes, Colors, Reactive, ModernGL
importall Compose

type GLVisualizePropertyState
    stroke::RGBA{Float32}
    fill::RGBA{Float32}
    stroke_dash::Array{Float32,1}
    stroke_linecap::Compose.LineCap
    stroke_linejoin::Compose.LineJoin
    visible::Bool
    linewidth::AbsoluteLength
    fontsize::AbsoluteLength
    font::AbstractString
    clip::Nullable{Compose.ClipPrimitive}
end
type GLVisualizePropertyFrame
    # Vector properties in this frame.
    vector_properties::Dict{Type, Compose.Property}

    # True if this property frame has scalar properties. Scalar properties are
    # emitted as a group (<g> tag) that must be closed when the frame is popped.
    has_scalar_properties::Bool

    function GLVisualizePropertyFrame()
        return new(Dict{Type, Compose.Property}(), false)
    end
end

type GLVisualizeBackend <: Compose.Backend
    screen
    stroke::RGBA{Float32}
    fill::RGBA{Float32}
    stroke_dash::Array{Float32,1}
    stroke_linecap::Compose.LineCap
    stroke_linejoin::Compose.LineJoin
    visible::Bool
    linewidth::AbsoluteLength
    fontsize::AbsoluteLength
    font::AbstractString
    clip::Nullable{Compose.ClipPrimitive}
    ppmm::Vec2f0

    # Keep track of property
    state_stack::Vector{GLVisualizePropertyState}
    property_stack::Vector{GLVisualizePropertyFrame}
    vector_properties::Dict{Type, Nullable{Compose.Property}}
    function GLVisualizeBackend(screen)
        img = new()
        img.screen = screen
        img.stroke = Compose.default_stroke_color == nothing ?
                        RGBA{Float32}(0, 0, 0, 0) : convert(RGBA{Float32},  Compose.default_stroke_color)
        img.fill   =  Compose.default_fill_color == nothing ?
                        RGBA{Float32}(0, 0, 0, 0) : convert(RGBA{Float32},  Compose.default_fill_color)
        img.stroke_dash = []
        img.stroke_linecap = Compose.LineCapButt()
        img.stroke_linejoin = Compose.LineJoinMiter()
        img.visible = true
        img.linewidth = Compose.default_line_width
        img.fontsize = Compose.default_font_size
        img.font = Compose.default_font_family
        img.clip = Nullable{Compose.ClipPrimitive}()

        img.state_stack = Array(GLVisualizePropertyState, 0)
        img.property_stack = Array(GLVisualizePropertyFrame, 0)
        img.vector_properties = Dict{Type, Nullable{Compose.Property}}()

        m = GLFW.GetPrimaryMonitor()
        props = GLWindow.MonitorProperties(m)
        img.ppmm = Vec2f0(props.dpi/25.4f0)

        img
    end

end
function Measures.width(img::GLVisualizeBackend)
    ((img.screen.inputs[:framebuffer_size].value[1]) / img.ppmm[1]) * mm
end

function Measures.height(img::GLVisualizeBackend)
    ((img.screen.inputs[:framebuffer_size].value[2]) / img.ppmm[2]) * mm
end

function Compose.root_box(img::GLVisualizeBackend)
    BoundingBox(Measures.width(img), Measures.height(img))
end
function absolute_native_units(img::GLVisualizeBackend, u::Float64)
    Float32(img.ppmm[1] * u)
end
function absolute_native_units{T}(img::GLVisualizeBackend, u::Length{:mm, T})
    Float32(img.ppmm[1] * u.value)
end
function absolute_native_units{T}(img::GLVisualizeBackend, u::Tuple{Length{:mm, T},Length{:mm, T}})
    Point2f0(u[1].value, (Measures.width(img)-u[2]).value).*Point2f0(img.ppmm)
end
relative_native_units{T}(img::GLVisualizeBackend, u::Tuple{Length{:mm, T},Length{:mm, T}}) = Point2f0(u[1].value, u[2].value).*Point2f0(img.ppmm)

export absolute_native_units

function Compose.push_property_frame(img::GLVisualizeBackend, properties::Vector{Compose.Property})
    if isempty(properties)
        return
    end

    frame = GLVisualizePropertyFrame()
    applied_properties = Set{Type}()
    scalar_properties = Array(Compose.Property, 0)
    for property in properties
        if Compose.isscalar(property) && !(typeof(property) in applied_properties)
            push!(scalar_properties, property)
            push!(applied_properties, typeof(property))
            frame.has_scalar_properties = true
        elseif !Compose.isscalar(property)
            frame.vector_properties[typeof(property)] = property
            img.vector_properties[typeof(property)] = property
        end
    end
    push!(img.property_stack, frame)
    if isempty(scalar_properties)
        return
    end

    save_property_state(img)
    for property in scalar_properties
        apply_property(img, property.primitives[1])
    end
end

function Compose.pop_property_frame(img::GLVisualizeBackend)
    @assert !isempty(img.property_stack)
    frame = pop!(img.property_stack)

    if frame.has_scalar_properties
        restore_property_state(img)
    end

    for (propertytype, property) in frame.vector_properties
        img.vector_properties[propertytype] = Nullable{Compose.Property}()
        for i in length(img.property_stack):-1:1
            if haskey(img.property_stack[i].vector_properties, propertytype)
                img.vector_properties[propertytype] =
                    img.property_stack[i].vector_properties[propertytype]
            end
        end
    end
end

# Return true if the vector properties need to be pushed and popped, rather
# than simply applied.
function vector_properties_require_push_pop(img::GLVisualizeBackend)
    for (propertytype, property) in img.vector_properties
        propertytype
        if in(propertytype, [Compose.Property{Compose.FontPrimitive},
                             Compose.Property{Compose.FontSizePrimitive},
                             Compose.Property{Compose.ClipPrimitive}])
            return true
        end
    end
    return false
end

function save_property_state(img::GLVisualizeBackend)
    push!(img.state_stack,
        GLVisualizePropertyState(
            img.stroke,
            img.fill,
            img.stroke_dash,
            img.stroke_linecap,
            img.stroke_linejoin,
            img.visible,
            img.linewidth,
            img.fontsize,
            img.font,
            img.clip))
end

function restore_property_state(img::GLVisualizeBackend)
    state = pop!(img.state_stack)
    img.stroke = state.stroke
    img.fill = state.fill
    img.stroke_dash = state.stroke_dash
    img.stroke_linecap = state.stroke_linecap
    img.stroke_linejoin = state.stroke_linejoin
    img.visible = state.visible
    img.linewidth = state.linewidth
    img.fontsize = state.fontsize
    img.font = state.font
    img.clip = state.clip
end

function push_vector_properties(img::GLVisualizeBackend, idx::Int)
    save_property_state(img)
    for (propertytype, property) in img.vector_properties
        if isnull(property)
            continue
        end
        primitives = get(property).primitives
        if idx > length(primitives)
            error("Vector form and vector property differ in length. Can't distribute.")
        end
        apply_property(img, primitives[idx])
    end
end

apply_property(img::GLVisualizeBackend, p) = nothing # ignore unsupported properties

function pop_vector_properties(img::GLVisualizeBackend)
    restore_property_state(img)
end

function apply_property(img::GLVisualizeBackend, p::Compose.StrokePrimitive)
    img.stroke = p.color
end

function apply_property(img::GLVisualizeBackend, p::Compose.FillPrimitive)
    img.fill = p.color
end

function apply_property(img::GLVisualizeBackend, p::Compose.FillOpacityPrimitive)
    img.fill = RGBA{Float32}(color(img.fill), p.value)
end

function apply_property(img::GLVisualizeBackend, p::Compose.StrokeOpacityPrimitive)
    img.stroke = RGBA{Float32}(color(img.stroke), p.value)
end

function apply_property(img::GLVisualizeBackend, p::Compose.StrokeDashPrimitive)
    img.stroke_dash = map(v -> absolute_native_units(img, v.value), p.value)
end

function apply_property(img::GLVisualizeBackend, p::Compose.StrokeLineCapPrimitive)
    img.stroke_linecap = p.value
end

function apply_property(img::GLVisualizeBackend, p::Compose.StrokeLineJoinPrimitive)
    img.stroke_linejoin = p.value
end

function apply_property(img::GLVisualizeBackend, p::Compose.VisiblePrimitive)
    img.visible = p.value
end

function apply_property(img::GLVisualizeBackend, property::Compose.LineWidthPrimitive)
    img.linewidth = property.value
end

function apply_property(img::GLVisualizeBackend, property::Compose.FontPrimitive)
    img.font = property.family
    #=
    font_desc = ccall((:pango_layout_get_font_description, Cairo._jl_libpango),
                      Ptr{Void}, (Ptr{Void},), img.ctx.layout)

    if font_desc == C_NULL
        size = absolute_native_units(img, default_font_size.value)
    else
        size = ccall((:pango_font_description_get_size, Cairo._jl_libpango),
                     Cint, (Ptr{Void},), font_desc)
    end

    Cairo.set_font_face(img.ctx,
        @sprintf("%s %0.2fpx", property.family, size / PANGO_SCALE))
    =#
end

function apply_property(img::GLVisualizeBackend, property::Compose.FontSizePrimitive)
    img.fontsize = property.value

    #=
    Cairo.set_font_face(img.ctx,
        @sprintf("%s %.2fpx",
            family, absolute_native_units(img, property.value.value)))
    =#
end

function apply_property(img::GLVisualizeBackend, property::Compose.ClipPrimitive)
    if isempty(property.points); return; end
end

Compose.isfinished(backend::GLVisualizeBackend) = false
Compose.finish(backend::GLVisualizeBackend) = nothing

Compose.iswithjs(img::GLVisualizeBackend) = false
Compose.iswithousjs(img::GLVisualizeBackend) = true

function Compose.draw{T <: Compose.CirclePrimitive}(img::GLVisualizeBackend, form::Compose.Form{T})
    r = form.primitives[1].radius
    radius = absolute_native_units(img, r)
    positions = Point2f0[absolute_native_units(img, elem.center) for elem in form.primitives]
    _view(visualize((Circle{Float32}(Point2f0(0), radius), positions),
    	color=img.fill,
    	stroke_color=img.stroke,
    	visible=img.visible
    ), img.screen, camera=:orthographic_pixel)
end
#=
function Compose.draw{T <: Compose.RectanglePrimitive}(img::GLVisualizeBackend, form::Compose.Form{T})
    positions = [absolute_native_units(img, elem.corner) for elem in form.primitives]
    wh = form.primitives[1].width, form.primitives[1].height
    scale = absolute_native_units(img, wh)
    _view(visualize(positions,
    	scale=scale,
    	shape=Cint(RECTANGLE),
    	model=translationmatrix(Vec3f0(scale/2, 0)),
    	style=Cint(FILLED),
    	color=img.fill,
    	stroke_color=img.stroke,
    	visible=img.visible
    ), img.screen)
end
=#
function Compose.draw{T}(img::GLVisualizeBackend, form::Compose.Form{T})

	if vector_properties_require_push_pop(img)
        for (idx, primitive) in enumerate(form.primitives)
            push_vector_properties(img, idx)
            draw(img, primitive)
            pop_vector_properties(img)
        end
    else
        for (idx, primitive) in enumerate(form.primitives)
            for (propertytype, property) in img.vector_properties
                if isnull(property)
                    continue
                end
                primitives = get(property).primitives
                if idx > length(primitives)
                    error("Vector form and vector property differ in length. Can't distribute.")
                end
                apply_property(img, primitives[idx])
            end
            Compose.draw(img, primitive)
        end
    end
end

function _add_line(points::Vector{Point2f0}, color, thickness)
    c = fill(RGBA{Float32}(color), length(points))
    robj      = visualize(points, :lines, color=c, thickness=thickness)
    _view(robj, camera=:orthographic_pixel)
end
function _add_line(points::LineSegment{Point2f0}, color, thickness)
    robj = visualize([points],
        color     = color,
        thickness = thickness
    )
    _view(robj, camera=:orthographic_pixel)
end
function add_line(points::Vector{Point2f0}, color, thickness)
    N = length(points)
    N <= 1 && return
    if N == 2
        _add_line(LineSegment{Point2f0}(points...), color, thickness)
    elseif N == 3
        error("Not implemented yet")
    else
        _add_line(points, color, thickness)
    end
end
function Compose.draw(img::GLVisualizeBackend, prim::Compose.LinePrimitive)
	N = length(prim.points)
    N <= 1 && return
    points = Point2f0[absolute_native_units(img, p) for p in prim.points]
    add_line(points, img.stroke, absolute_native_units(img, img.linewidth))
end

function Compose.draw(img::GLVisualizeBackend, form::Compose.FormBatch)
    println("FormBatch: ", form)
end

let positions    = 0,
    scales       = 0,
    fillcolors   = 0,
    strokecolors = 0,
    indices      = 0

function Compose.draw(img::GLVisualizeBackend, prim::Compose.RectanglePrimitive)
    wh = relative_native_units(img, (prim.width, prim.height))
    xy = absolute_native_units(img, prim.corner)-Point2f0(0, wh[2])
    stroke_width = img.stroke.alpha > 0 ? 5f0 : 0f0
    fillcolor = img.fill
    stroke_color = img.stroke

        positions = GLBuffer([xy])
        scales    = GLBuffer([wh])
        fillcolors = GLBuffer([fillcolor])
        strokecolors = GLBuffer([stroke_color])
        _view(visualize((RECTANGLE, positions),
        	color=fillcolors,
            scale=scales,
        	stroke_color=strokecolors,
            stroke_width=stroke_width,
        ), img.screen, camera=:orthographic_pixel)

end
end

function Compose.draw(img::GLVisualizeBackend, prim::Compose.PolygonPrimitive)
    warn("PolygonPrimitive NOT IMPLEMENTED YET")
    # NOT SUPPORTED YET
end

function Compose.draw(img::GLVisualizeBackend, prim::Compose.EllipsePrimitive)
	warn("EllipsePrimitive NOT IMPLEMENTED YET")
    # NOT SUPPORTED YET
end

function Compose.draw(img::GLVisualizeBackend, prim::Compose.CurvePrimitive)
	warn("CURVE NOT IMPLEMENTED YET")
    # NOT SUPPORTED YET
end

function Compose.draw(img::GLVisualizeBackend, prim::Compose.BitmapPrimitive)
    xyz = Vec3f0(prim.corner.x.abs, prim.corner.y.abs, 0)
    scale = Vec3f0(prim.width.abs, prim.height.abs, 1)
    _view(visualize(colorim(prim.data), model=translationmatrix(xyz)*scalematrix(scale)), img.screen, camera=:orthographic_pixel)
end

function gen_text(text, atlas, font, position, scale, offset)
    pos = GLVisualize.calc_position(text, position, scale, font, atlas)
    s = Vec2f0[GLVisualize.glyph_scale!(atlas, c, font) .* scale for c in text]
    off = GLVisualize.calc_offset(text, scale, font, atlas) .+ Point2f0(0, offset)
    uvwidth = Vec4f0[GLVisualize.glyph_uv_width!(atlas, c, font) for c in text]
    pos, s, off, uvwidth
end

function pango_parse_markup(text)
    c_stripped_text = Array(Ptr{UInt8}, 1)
    c_attr_list = Array(Ptr{Void}, 1)
    ret = ccall((:pango_parse_markup, Compose.libpango),
        Int32,
        (Ptr{UInt8}, Int32, UInt32, Ptr{Ptr{Void}},
        Ptr{Ptr{UInt8}}, Ptr{UInt32}, Ptr{Void}),
        Compat.String(text), -1, 0, c_attr_list, c_stripped_text,
        C_NULL, C_NULL
    )
    if ret == 0
        error("Could not parse pango markup.")
    end
    str = Compat.String(c_stripped_text[1])

    # TODO: do c_stripped_text and c_attr_list need to be freed?

    return convert(Vector{UInt8}, str), c_attr_list[]
end

function parse_pango(text::AbstractString, scale)
    text, c_attr_list = pango_parse_markup(text)
    GLV = GLVisualize

    atlas          = GLV.get_texture_atlas()
    font           = GLV.DEFAULT_FONT_FACE

    last_idx = 1
    last_offset = 0.0
    last_scale = Vec2f0(1.0)
    last_position = Point2f0(0)

    positions       = Point2f0[]
    scales          = Vec2f0[]
    offset          = Point2f0[]
    uv_offset_width = Vec4f0[]
    for (idx, attr) in Compose.unpack_pango_attr_list(c_attr_list)
        current_text = Compat.String(text[last_idx:idx])
        last_idx = idx+1

        pos, sa, oa, uvwidth = gen_text(current_text,
            atlas, font, last_position, last_scale.*scale, last_offset
        )
        append!(scales, sa)
        append!(positions, pos)
        append!(offset, oa)
        append!(uv_offset_width, uvwidth)
        last_position = GLVisualize.calc_position(
            last(positions), Point2f0(0), atlas,
            last(current_text), font, last_scale.*scale
        )

        last_offset = 0.0
        if !(attr.rise === nothing)
            last_offset = ((attr.rise / Compose.PANGO_SCALE))
        end
        last_scale = Vec2f0(1)
        if !(attr.scale === nothing)
            last_scale = Vec2f0(attr.scale)
        end
    end
    last_scale = Vec2f0(1)
    if last_idx <= length(text)
        current_text = Compat.String(text[last_idx:end])
        pos, sa, oa, uvwidth = gen_text(current_text,
            atlas, font, last_position, last_scale.*scale, last_offset
        )
        append!(scales, sa)
        append!(positions, pos)
        append!(offset, oa)
        append!(uv_offset_width, uvwidth)
    end
    positions, scales, offset, uv_offset_width
end

function Compose.draw(img::GLVisualizeBackend, prim::Compose.TextPrimitive)
	pos 	= absolute_native_units(img, prim.position)
	s1 		= absolute_native_units(img, img.fontsize)/20f0
	s 		= Vec2f0(s1)
    positions, scales, offset, uv_offset_width = parse_pango(prim.value, s)
    atlas = GLVisualize.get_texture_atlas();
	obj 	= visualize((DISTANCEFIELD, positions),
        distancefield=atlas.images,
        scale=scales,
        offset=offset,
        uv_offset_width=uv_offset_width,
        color=img.fill
    );
	bb 		= value(GLAbstraction.boundingbox(obj))
    w,h,_   = widths(bb)
	x,y,_ 	= minimum(bb)
    pos     -= Point2f0(x, y)
    transmat = eye(Mat{4,4,Float32})
	if prim.rot.theta != 0.0
        pivot = Vec3f0(absolute_native_units(img, prim.rot.offset), 0)
        rot = GLAbstraction.rotationmatrix_z(Float32(-prim.rot.theta))
        transmat *= translationmatrix(pivot)*rot*translationmatrix(-pivot)
    end
    if prim.halign != Compose.hleft || prim.valign != Compose.vbottom
        if prim.halign == Compose.hcenter
            pos = (pos[1] - (w/2f0), pos[2])
        elseif prim.halign == Compose.hright
            pos = (pos[1] - w, pos[2])
        end
        if prim.valign == Compose.vcenter
            pos = (pos[1], pos[2] - h/2f0)
        elseif prim.valign == Compose.vtop
            pos = (pos[1], pos[2] - h)
        end
    end
    transmat *= translationmatrix(Vec3f0(pos..., 0))
	GLAbstraction.transformation(obj, transmat)
	_view(obj, img.screen, camera=:orthographic_pixel)
end

end


text:
function trace_to_seq(trace::Trace, incl_disconnect=true)
    sequence = map(enumerate(trace[:strokes])) do (i, (dir, x))
        s = ["D", "U", "L", "R"][dir]
        if incl_disconnect && trace[:strokes => i => :disconnect]
            s = "#" * s
        end
        return s
    end
    return join(sequence)
end

function stroke_edit_dist(trace::Trace, sequence::AbstractString, incl_disconnect=true)
    if !incl_disconnect sequence = filter(c -> c != '#', sequence) end
    return evaluate(Levenshtein(), trace_to_seq(trace, incl_disconnect), sequence)
end

function stroke_edit_dist(state::Gen.ParticleFilterState, sequence::AbstractString,
                          incl_disconnect=true)
    dists = (stroke_edit_dist(tr, sequence, incl_disconnect) for tr in get_traces(state))
    return sum(get_norm_weights(state) .* dists)
end

function stroke_rmse(trace::Trace, targets::AbstractMatrix)
    points = reduce(hcat, last.(trace[:strokes]))
    T = min(size(points)[2], size(targets)[2])
    mse = sum((points[:,1:T] - targets[:,1:T]).^2) ./ T
    return mse ^ 0.5
end

function stroke_rmse(state::Gen.ParticleFilterState, targets::AbstractMatrix)
    rmses = (stroke_rmse(tr, targets) for tr in get_traces(state))
    return sum(get_norm_weights(state) .* rmses)
end

function log_ml(state::Gen.ParticleFilterState)
    log_lls = (project(tr, select(:obs)) for tr in get_traces(state))
    return sum(get_norm_weights(state) .* log_lls)
end


text:
include("utils.jl")

@testset "Oblique manifold" begin
    M = Oblique(3, 2)
    @testset "Oblique manifold Basics" begin
        @test Sphere(2)^3 === Oblique(3, 3)
        @test Sphere(2)^(3,) === PowerManifold(Sphere(2), 3)
        @test ^(Sphere(2), 2) === Oblique(3, 2)
        @test typeof(^(Sphere(2), 2)) == Oblique{3,2,ℝ,2}
        @test repr(M) == "Oblique(3,2; field = ℝ)"
        @test representation_size(M) == (3, 2)
        @test manifold_dimension(M) == 4
        p = [2, 0, 0]
        p2 = [p p]
        @test !is_point(M, p)
        @test_throws DomainError is_point(M, p, true)
        @test !is_point(M, p2)
        @test_throws CompositeManifoldError is_point(M, p2, true)
        @test !is_vector(M, p2, 0.0)
        @test_throws CompositeManifoldError{ComponentManifoldError{Int64,DomainError}} is_vector(
            M,
            p2,
            [0.0, 0.0, 0.0],
            true,
        )
        @test !is_vector(M, p2, [0.0, 0.0, 0.0])
        @test_throws DomainError is_vector(M, p, [0.0, 0.0, 0.0], true) # p wrong
        @test injectivity_radius(M) ≈ π
        x = [1.0 0.0 0.0; 1.0 0.0 0.0]'
        @test_throws DomainError is_vector(M, x, [0.0, 0.0, 0.0], true) # tangent wrong
        y = [1.0 0.0 0.0; 1/sqrt(2) 1/sqrt(2) 0.0]'
        z = [1/sqrt(2) 1/sqrt(2) 0.0; 1.0 0.0 0.0]'
        basis_types = (DefaultOrthonormalBasis(),)
        transports = [ParallelTransport(), PowerVectorTransport(ParallelTransport())]
        test_manifold(
            M,
            [x, y, z],
            test_forward_diff=false,
            test_reverse_diff=false,
            test_vector_spaces=true,
            test_project_tangent=false,
            test_musical_isomorphisms=true,
            test_default_vector_transport=true,
            vector_transport_methods=transports,
            basis_types_to_from=basis_types,
            exp_log_atol_multiplier=1,
        )
    end
end


text:
"""
    VectorSpaceType

Abstract type for tangent spaces, cotangent spaces, their tensor products,
exterior products, etc.

Every vector space `VS` is supposed to provide:
* a method of constructing vectors,
* basic operations: addition, subtraction, multiplication by a scalar
  and negation (unary minus),
* [`zero_vector!(VS, v, x)`](@ref) to construct zero vectors at point `x`,
* `similar(v, T)` for vector `v`,
* `copyto!(v, w)` for vectors `v` and `w`,
* `eltype(v)` for vector `v`,
* [`vector_space_dimension(::VectorBundleFibers{<:typeof(VS)}) where VS`](@ref).

Optionally:
* inner product via `inner` (used to provide Riemannian metric on vector
  bundles),
* [`flat!`](@ref) and [`sharp!`](@ref),
* `norm` (by default uses `inner`),
* [`project_vector!`](@ref) (for embedded vector spaces),
* [`representation_size`](@ref) (if support for [`ProductArray`](@ref) is desired),
* broadcasting for basic operations.
"""
abstract type VectorSpaceType end

struct TangentSpaceType <: VectorSpaceType end
struct CotangentSpaceType <: VectorSpaceType end

TCoTSpaceType = Union{TangentSpaceType, CotangentSpaceType}

const TangentSpace = TangentSpaceType()
const CotangentSpace = CotangentSpaceType()

"""
    TensorProductType(spaces::VectorSpaceType...)

Vector space type corresponding to the tensor product of given vector space
types.
"""
struct TensorProductType{TS<:Tuple} <: VectorSpaceType
    spaces::TS
end

TensorProductType(spaces::VectorSpaceType...) = TensorProductType{typeof(spaces)}(spaces)

"""
    VectorBundleFibers(VS::VectorSpaceType, M::Manifold)

Type representing a family of vector spaces (fibers) of a vector bundle over `M`
with vector spaces of type `VS`. In contrast with `VectorBundle`, operations
on `VectorBundleFibers` expect point-like and vector-like parts to be
passed separately instead of being bundled together. It can be thought of
as a representation of vector spaces from a vector bundle but without
storing the point at which a vector space is attached (which is specified
separately in various functions).
"""
struct VectorBundleFibers{TVS<:VectorSpaceType, TM<:Manifold}
    VS::TVS
    M::TM
end

base_manifold(B::VectorBundleFibers) = base_manifold(B.M)

"""
    VectorSpaceAtPoint(fiber::VectorBundleFibers, x)

A vector space (fiber type `fiber` of a vector bundle) at point `x` from
the manifold `fiber.M`.
"""
struct VectorSpaceAtPoint{TFiber<:VectorBundleFibers, TX}
    fiber::TFiber
    x::TX
end

base_manifold(B::VectorSpaceAtPoint) = base_manifold(B.fiber)

"""
    TangentSpaceAtPoint(M::Manifold, x)

Return an object of type [`VectorSpaceAtPoint`](@ref) representing tangent
space at `x`.
"""
function TangentSpaceAtPoint(M::Manifold, x)
    return VectorSpaceAtPoint(TangentBundleFibers(M), x)
end

"""
    CotangentSpaceAtPoint(M::Manifold, x)

Return an object of type [`VectorSpaceAtPoint`](@ref) representing cotangent
space at `x`.
"""
function CotangentSpaceAtPoint(M::Manifold, x)
    return VectorSpaceAtPoint(CotangentBundleFibers(M), x)
end

"""
    FVector(type::VectorSpaceType, data)

Decorator indicating that the vector `data` from a fiber of a vector bundle
is of type `type`.
"""
struct FVector{TType<:VectorSpaceType,TData}
    type::TType
    data::TData
end

size(x::FVector) = size(x.data)
Base.@propagate_inbounds getindex(x::FVector, i) = getindex(x.data, i)
Base.@propagate_inbounds setindex!(x::FVector, val, i) = setindex!(x.data, val, i)

(+)(v1::FVector, v2::FVector) = FVector(v1.type, v1.data + v2.data)
(-)(v1::FVector, v2::FVector) = FVector(v1.type, v1.data - v2.data)
(-)(v::FVector) = FVector(v.type, -v.data)
(*)(a::Number, v::FVector) = FVector(v.type, a*v.data)

eltype(::Type{FVector{TType,TData}}) where {TType<:VectorSpaceType,TData} = eltype(TData)
similar(x::FVector) = FVector(x.type, similar(x.data))
similar(x::FVector, ::Type{T}) where {T} = FVector(x.type, similar(x.data, T))

@doc doc"""
    flat!(M::Manifold, v::FVector, x, w::FVector)

Compute the flat isomorphism (one of the musical isomorphisms) of vector `w`
from the vector space of type `M` at point `x` from manifold `M.M`.

The function can be used for example to transform vectors
from the tangent bundle to vectors from the cotangent bundle
$\flat \colon TM \to T^{*}M$
"""
function flat!(M::Manifold, v::FVector, x, w::FVector)
    error("flat! not implemented for vector bundle fibers space " *
        "of type $(typeof(M)), vector of type $(typeof(v)), point of " *
        "type $(typeof(x)) and vector of type $(typeof(w)).")
end

function flat(M::Manifold, x, w::FVector)
    v = similar_result(M, flat, w, x)
    flat!(M, v, x, w)
    return v
end

function similar_result(M::Manifold, ::typeof(flat), w::FVector{TangentSpaceType}, x)
    return FVector(CotangentSpace, similar(w.data))
end

@doc doc"""
    sharp!(M::Manifold, v::FVector, x, w::FVector)

Compute the sharp isomorphism (one of the musical isomorphisms) of vector `w`
from the vector space of type `M` at point `x` from manifold `M.M`.

The function can be used for example to transform vectors
from the cotangent bundle to vectors from the tangent bundle
$\sharp \colon T^{*}M \to TM$
"""
function sharp!(M::Manifold, v::FVector, x, w::FVector)
    error("sharp! not implemented for vector bundle fibers space " *
        "of type $(typeof(M)), vector of type $(typeof(v)), point of " *
        "type $(typeof(x)) and vector of type $(typeof(w)).")
end

function sharp(M::Manifold, x, w::FVector)
    v = similar_result(M, sharp, w, x)
    sharp!(M, v, x, w)
    return v
end

function similar_result(M::Manifold, ::typeof(sharp), w::FVector{CotangentSpaceType}, x)
    return FVector(TangentSpace, similar(w.data))
end

"""
    similar_result_type(B::VectorBundleFibers, f, args::NTuple{N,Any}) where N

Returns type of element of the array that will represent the result of
function `f` for representing an operation with result in the vector space `VS`
for manifold `M` on given arguments (passed at a tuple).
"""
function similar_result_type(B::VectorBundleFibers, f, args::NTuple{N,Any}) where N
    T = typeof(reduce(+, one(eltype(eti)) for eti ∈ args))
    return T
end

"""
    similar_result(B::VectorBundleFibers, f, x...)

Allocates an array for the result of function `f` that is
an element of the vector space of type `B.VS` on manifold `B.M`
and arguments `x...` for implementing the non-modifying operation
using the modifying operation.
"""
function similar_result(B::VectorBundleFibers, f, x...)
    T = similar_result_type(B, f, x)
    return similar(x[1], T)
end

"""
    distance(B::VectorBundleFibers, x, v, w)

Distance between vectors `v` and `w` from the vector space at point `x`
from the manifold `M.M`, that is the base manifold of `M`.
"""
distance(B::VectorBundleFibers, x, v, w) = norm(B, x, v-w)

"""
    inner(B::VectorBundleFibers, x, v, w)

Inner product of vectors `v` and `w` from the vector space of type `B.VS`
at point `x` from manifold `B.M`.
"""
function inner(B::VectorBundleFibers, x, v, w)
    error("inner not defined for vector space family of type $(typeof(B)), " *
        "point of type $(typeof(x)) and " *
        "vectors of types $(typeof(v)) and $(typeof(w)).")
end

function inner(B::VectorBundleFibers{<:TangentSpaceType}, x, v, w)
    return inner(B.M, x, v, w)
end

function inner(B::VectorBundleFibers{<:CotangentSpaceType}, x, v, w)
    return inner(B.M,
                 x,
                 sharp(B.M, x, FVector(CotangentSpace, v)).data,
                 sharp(B.M, x, FVector(CotangentSpace, w)).data)
end

"""
    norm(B::VectorBundleFibers, x, v)

Norm of the vector `v` from the vector space of type `B.VS`
at point `x` from manifold `B.M`.
"""
norm(B::VectorBundleFibers, x, v) = sqrt(inner(B, x, v, v))

norm(B::VectorBundleFibers{<:TangentSpaceType}, x, v) = norm(B.M, x, v)

"""
    vector_space_dimension(B::VectorBundleFibers)

Dimension of the vector space of type `B`.
"""
function vector_space_dimension(B::VectorBundleFibers)
    error("vector_space_dimension not implemented for vector space family $(typeof(B)).")
end

vector_space_dimension(B::VectorBundleFibers{<:TCoTSpaceType}) = manifold_dimension(B.M)

function vector_space_dimension(B::VectorBundleFibers{<:TensorProductType})
    dim = 1
    for space in B.VS.spaces
        dim *= vector_space_dimension(VectorBundleFibers(space, B.M))
    end
    return dim
end

function representation_size(B::VectorBundleFibers{<:TCoTSpaceType})
    representation_size(B.M)
end

"""
    zero_vector!(B::VectorBundleFibers, v, x)

Save the zero vector from the vector space of type `B.VS` at point `x`
from manifold `B.M` to `v`.
"""
function zero_vector!(B::VectorBundleFibers, v, x)
    error("zero_vector! not implemented for vector space family of type $(typeof(B)).")
end

function zero_vector!(B::VectorBundleFibers{<:TangentSpaceType}, v, x)
    zero_tangent_vector!(B.M, v, x)
    return v
end

"""
    zero_vector(B::VectorBundleFibers, x)

Compute the zero vector from the vector space of type `B.VS` at point `x`
from manifold `B.M`.
"""
function zero_vector(B::VectorBundleFibers, x)
    v = similar_result(B, zero_vector, x)
    zero_vector!(B, v, x)
    return v
end

"""
    project_vector!(B::VectorBundleFibers, v, x, w)

Project vector `w` from the vector space of type `B.VS` at point `x`
and save the result to `v`.
"""
function project_vector!(B::VectorBundleFibers, v, x, w)
    error("project_vector! not implemented for vector space family of type $(typeof(B)), output vector of type $(typeof(v)) and input vector at point $(typeof(x)) with type of w $(typeof(w)).")
end

function project_vector!(B::VectorBundleFibers{<:TangentSpaceType}, v, x, w)
    project_tangent!(B.M, v, x, w)
    return v
end

"""
    VectorBundle(M::Manifold, type::VectorSpaceType)

Vector bundle on manifold `M` of type `type`.
"""
struct VectorBundle{TVS<:VectorSpaceType, TM<:Manifold} <: Manifold
    type::TVS
    M::TM
    VS::VectorBundleFibers{TVS, TM}
end

base_manifold(B::VectorBundle) = base_manifold(B.M)

function VectorBundle(VS::TVS, M::TM) where {TVS<:VectorSpaceType, TM<:Manifold}
    return VectorBundle{TVS, TM}(VS, M, VectorBundleFibers(VS, M))
end

function representation_size(B::VectorBundle)
    len_manifold = prod(representation_size(B.M))
    len_vs = prod(representation_size(B.VS))
    return (len_manifold + len_vs,)
end

TangentBundleFibers{M} = VectorBundleFibers{TangentSpaceType,M}
TangentBundleFibers(M::Manifold) = VectorBundleFibers(TangentSpace, M)

CotangentBundleFibers{M} = VectorBundleFibers{CotangentSpaceType,M}
CotangentBundleFibers(M::Manifold) = VectorBundleFibers(CotangentSpace, M)

TangentBundle{M} = VectorBundle{TangentSpaceType,M}
TangentBundle(M::Manifold) = VectorBundle(TangentSpace, M)

CotangentBundle{M} = VectorBundle{CotangentSpaceType,M}
CotangentBundle(M::Manifold) = VectorBundle(CotangentSpace, M)

"""
    bundle_projection(B::VectorBundle, x::ProductRepr)

Projection of point `x` from the bundle `M` to the base manifold.
Returns the point on the base manifold `B.M` at which the vector part
of `x` is attached.
"""
function bundle_projection(B::VectorBundle, x)
    return submanifold_component(x, 1)
end

function isapprox(B::VectorBundle, x, y; kwargs...)
    return isapprox(B.M, x.parts[1], y.parts[1]; kwargs...) &&
        isapprox(x.parts[2], y.parts[2]; kwargs...)
end

function isapprox(B::VectorBundle, x, v, w; kwargs...)
    return isapprox(B.M, v.parts[1], w.parts[1]; kwargs...) &&
        isapprox(B.M, x.parts[1], v.parts[2], w.parts[2]; kwargs...)
end

manifold_dimension(B::VectorBundle) = manifold_dimension(B.M) + vector_space_dimension(B.VS)

@doc doc"""
    inner(B::VectorBundle, x, v, w)

Inner product of tangent vectors `v` and `w` at point `x` from the
vector bundle `B` over manifold `B.VS` (denoted $M$).

Notation:
  * The point $x = (p_x, \xi_x)$ where $p_x \in M$ and $\xi_x$ belongs to the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.
  * The tangent vector $v = (\xi_{v,M}, \xi_{v,F}) \in T_{x}B$ where
    $\xi_{v,M}$ is a tangent vector from the tangent space $T_{p_x}M$ and
    $\xi_{v,F}$ is a tangent vector from the tangent space $T_{\xi_x}F$ (isomorphic to $F$).
    Similarly for the other tangent vector $w = (\xi_{w,M}, \xi_{w,F}) \in T_{x}B$.

The inner product is calculated as

$\langle v, w \rangle_{B} = \langle \xi_{v,M}, \xi_{w,M} \rangle_{M} + \langle \xi_{v,F}, \xi_{w,F} \rangle_{F}.$
"""
function inner(B::VectorBundle, x, v, w)
    return inner(B.M, x.parts[1], v.parts[1], w.parts[1]) +
           inner(B.VS, x.parts[2], v.parts[2], w.parts[2])
end

@doc doc"""
    distance(B::VectorBundle, x, y)

Distance between points $x$ and $y$ from the
vector bundle `B` over manifold `B.VS` (denoted $M$).

Notation:
  * The point $x = (p_x, \xi_x)$ where $p_x \in M$ and $\xi_x$ belongs to the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.
    Similarly, $y = (p_y, \xi_y)$.

The distance is calculated as

$d_B(x, y) = \sqrt{d_M(p_x, p_y)^2 + d_F(\xi_x, \xi_{y\to x})^2}$

where $d_M$ is the distance on manifold $M$, $d_F$ is the distance
between two vectors from the fiber $F$ and $\xi_{y\to x}$ is the result
of parallel transport of vector $\xi_y$ to point $p_x$. The default
behavior of [`vector_transport_to`](@ref) is used to compute the vector
transport.
"""
function distance(B::VectorBundle, x, y)
    dist_man = distance(B.M, x.parts[1], y.parts[1])
    vy_x = vector_transport_to(B.M, y.parts[1], y.parts[2], x.parts[1])
    dist_vec = distance(B.VS, x.parts[1], x.parts[2], vy_x)

    return sqrt(dist_man^2 + dist_vec^2)
end

@doc doc"""
    exp!(B::VectorBundle, y, x, v)

Exponential map of tangent vector $v$ at point $x$ from
vector bundle `B` over manifold `B.VS` (denoted $M$).
The result overwrites the point `y`.

Notation:
  * The point $x = (p_x, \xi_x)$ where $p_x \in M$ and $\xi_x$ belongs to the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.
  * The tangent vector $v = (\xi_{v,M}, \xi_{v,F}) \in T_{x}B$ where
    $\xi_{v,M}$ is a tangent vector from the tangent space $T_{p_x}M$ and
    $\xi_{v,F}$ is a tangent vector from the tangent space $T_{\xi_x}F$ (isomorphic to $F$).

The exponential map is calculated as

$\exp_{x}(v) = (\exp_{p_x}(\xi_{v,M}), \xi_{\exp})$

where $\xi_{\exp}$ is the result of vector transport of $\xi_x + \xi_{v,F}$
to the point $\exp_{p_x}(\xi_{v,M})$.
The sum $\xi_x + \xi_{v,F}$ corresponds to the exponential map in the vector space $F$.
"""
function exp!(B::VectorBundle, y, x, v)
    exp!(B.M, y.parts[1], x.parts[1], v.parts[1])
    vector_transport_to!(B.M, y.parts[2], x.parts[1], x.parts[2] + v.parts[2], y.parts[1])
    return y
end

@doc doc"""
    log!(B::VectorBundle, v, x, y)

Logarithmic map of the point $y$ at point $x$ from
vector bundle `B` over manifold `B.VS` (denoted $M$).
The result overwrites the tangent vector `v`.

Notation:
  * The point $x = (p_x, \xi_x)$ where $p_x \in M$ and $\xi_x$ belongs to the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.
    Similarly, $y = (p_y, \xi_y)$.

The logarithmic map is calculated as

$\log_{x}(y) = (\log_{p_x}(p_y), \xi_{\log} - \xi_x)$

where $\xi_{\log}$ is the result of vector transport of $\xi_y$
to the point $p_x$.
The difference $\xi_{\log} - \xi_x$ corresponds to the logarithmic map in the vector space $F$.
"""
function log!(B::VectorBundle, v, x, y)
    log!(B.M, v.parts[1], x.parts[1], y.parts[1])
    vector_transport_to!(B.M, v.parts[2], y.parts[1], y.parts[2], x.parts[1])
    copyto!(v.parts[2], v.parts[2] - x.parts[2])
    return v
end

@doc doc"""
    zero_tangent_vector!(B::VectorBundle, v, x)

Zero tangent vector at point $x$ from the vector bundle `B`
over manifold `B.VS` (denoted $M$). The zero vector belongs to the space $T_{x}B$
The result overwrites the tangent vector `v`.

Notation:
  * The point $x = (p_x, \xi_x)$ where $p_x \in M$ and $\xi_x$ belongs to the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.

The zero vector is calculated as

$\mathbf{0}_{x} = (\mathbf{0}_{p_x}, \mathbf{0}_F)$

where $\mathbf{0}_{p_x}$ is the zero tangent vector from $T_{p_x}M$ and
$\mathbf{0}_F$ is the zero element of the vector space $F$.
"""
function zero_tangent_vector!(B::VectorBundle, v, x)
    zero_tangent_vector!(B.M, v.parts[1], x.parts[1])
    zero_vector!(B.VS, v.parts[2], x.parts[2])
    return v
end

@doc doc"""
    project_point!(B::VectorBundle, x)

Project the point $x$ from the ambient space of the vector bundle `B`
over manifold `B.VS` (denoted $M$) to the vector bundle.

Notation:
  * The point $x = (p_x, \xi_x)$ where $p_x$ belongs to the ambient space of $M$
    and $\xi_x$ belongs to the ambient space of the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.

The projection is calculated by projecting the point $p_x$ to the manifold $M$
and then projecting the vector $\xi_x$ to the tangent space $T_{p_x}M$.
"""
function project_point!(B::VectorBundle, x)
    project_point!(B.M, x.parts[1])
    project_tangent!(B.M, x.parts[2], x.parts[1], x.parts[2])
    return x
end

@doc doc"""
    project_tangent!(B::VectorBundle, w, x, v)

Project the element $v$ of the ambient space of the tangent space $T_x B$
to the tangent space $T_x B$.
The results overwrites $w$.

Notation:
  * The point $x = (p_x, \xi_x)$ where $p_x \in M$ and $\xi_x$ belongs to the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.
  * The vector $x = (\xi_{v,M}, \xi_{v,F})$ where $p_x$ belongs to the ambient space of $T_{p_x}M$
    and $\xi_{v,F}$ belongs to the ambient space of the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.

The projection is calculated by projecting $\xi_{v,M}$ to tangent space $T_{p_x}M$
and then projecting the vector $\xi_{v,F}$ to the fiber $F$.
"""
function project_tangent!(B::VectorBundle, w, x, v)
    project_tangent!(B.M, w.parts[1], x.parts[1], v.parts[1])
    project_tangent!(B.M, w.parts[2], x.parts[1], v.parts[2])
    return w
end