text:
"""
    pyis(x, y)

True if `x` and `y` are the same Python object. Equivalent to `x is y` in Python.
"""
pyis(x, y) = @autopy x y getptr(x_) == getptr(y_)
export pyis

pyisnot(x, y) = !pyis(x, y)

"""
    pyrepr(x)

Equivalent to `repr(x)` in Python.
"""
pyrepr(x) = pynew(errcheck(@autopy x C.PyObject_Repr(getptr(x_))))
pyrepr(::Type{String}, x) = (s=pyrepr(x); ans=pystr_asstring(s); pydel!(s); ans)
export pyrepr

"""
    pyascii(x)

Equivalent to `ascii(x)` in Python.
"""
pyascii(x) = pynew(errcheck(@autopy x C.PyObject_ASCII(getptr(x_))))
pyascii(::Type{String}, x) = (s=pyascii(x); ans=pystr_asstring(s); pydel!(s); ans)
export pyascii

"""
    pyhasattr(x, k)

Equivalent to `hasattr(x, k)` in Python.
"""
pyhasattr(x, k) = errcheck(@autopy x k C.PyObject_HasAttr(getptr(x_), getptr(k_))) == 1
export pyhasattr

"""
    pygetattr(x, k, [d])

Equivalent to `getattr(x, k)` or `x.k` in Python.

If `d` is specified, it is returned if the attribute does not exist.
"""
pygetattr(x, k) = pynew(errcheck(@autopy x k C.PyObject_GetAttr(getptr(x_), getptr(k_))))
function pygetattr(x, k, d)
    ptr = @autopy x k C.PyObject_GetAttr(getptr(x_), getptr(k_))
    if iserrset(ptr)
        errclear()
        return d
    else
        return pynew(ptr)
    end
end
export pygetattr

"""
    pysetattr(x, k, v)

Equivalent to `setattr(x, k, v)` or `x.k = v` in Python.
"""
pysetattr(x, k, v) = (errcheck(@autopy x k v C.PyObject_SetAttr(getptr(x_), getptr(k_), getptr(v_))); nothing)
export pysetattr

"""
    pydelattr(x, k)

Equivalent to `delattr(x, k)` or `del x.k` in Python.
"""
pydelattr(x, k) = (errcheck(@autopy x k C.PyObject_SetAttr(getptr(x_), getptr(k_), C.PyNULL)); nothing)
export pydelattr

"""
    pyissubclass(s, t)

Test if `s` is a subclass of `t`. Equivalent to `issubclass(s, t)` in Python.
"""
pyissubclass(s, t) = errcheck(@autopy s t C.PyObject_IsSubclass(getptr(s_), getptr(t_))) == 1
export pyissubclass

"""
    pyisinstance(x, t)

Test if `x` is of type `t`. Equivalent to `isinstance(x, t)` in Python.
"""
pyisinstance(x, t) = errcheck(@autopy x t C.PyObject_IsInstance(getptr(x_), getptr(t_))) == 1
export pyisinstance

"""
    pyhash(x)

Equivalent to `hash(x)` in Python, converted to an `Integer`.
"""
pyhash(x) = errcheck(@autopy x C.PyObject_Hash(getptr(x_)))
export pyhash

"""
    pytruth(x)

The truthyness of `x`. Equivalent to `bool(x)` in Python, converted to a `Bool`.
"""
pytruth(x) = errcheck(@autopy x C.PyObject_IsTrue(getptr(x_))) == 1
export pytruth

"""
    pynot(x)

The falsyness of `x`. Equivalent to `not x` in Python, converted to a `Bool`.
"""
pynot(x) = errcheck(@autopy x C.PyObject_Not(getptr(x_))) == 1
export pynot

"""
    pylen(x)

The length of `x`. Equivalent to `len(x)` in Python, converted to an `Integer`.
"""
pylen(x) = errcheck(@autopy x C.PyObject_Length(getptr(x_)))
export pylen

"""
    pygetitem(x, k)

Equivalent to `getitem(x, k)` or `x[k]` in Python.
"""
pygetitem(x, k) = pynew(errcheck(@autopy x k C.PyObject_GetItem(getptr(x_), getptr(k_))))
export pygetitem

"""
    pysetitem(x, k, v)

Equivalent to `setitem(x, k, v)` or `x[k] = v` in Python.
"""
pysetitem(x, k, v) = (errcheck(@autopy x k v C.PyObject_SetItem(getptr(x_), getptr(k_), getptr(v_))); nothing)
export pysetitem

"""
    pydelitem(x, k)

Equivalent to `delitem(x, k)` or `del x[k]` in Python.
"""
pydelitem(x, k) = (errcheck(@autopy x k C.PyObject_DelItem(getptr(x_), getptr(k_))); nothing)
export pydelitem

"""
    pydir(x)

Equivalent to `dir(x)` in Python.
"""
pydir(x) = pynew(errcheck(@autopy x C.PyObject_Dir(getptr(x_))))
export pydir

pycallargs(f) = pynew(errcheck(@autopy f C.PyObject_CallObject(getptr(f_), C.PyNULL)))
pycallargs(f, args) = pynew(errcheck(@autopy f args C.PyObject_CallObject(getptr(f_), getptr(args_))))
pycallargs(f, args, kwargs) = pynew(errcheck(@autopy f args kwargs C.PyObject_Call(getptr(f_), getptr(args_), getptr(kwargs_))))

"""
    pycall(f, args...; kwargs...)

Call the Python object `f` with the given arguments.
"""
pycall(f, args...; kwargs...) =
    if !isempty(kwargs)
        args_ = pytuple_fromiter(args)
        kwargs_ = pystrdict_fromiter(kwargs)
        ans = pycallargs(f, args_, kwargs_)
        pydel!(args_)
        pydel!(kwargs_)
        ans
    elseif !isempty(args)
        args_ = pytuple_fromiter(args)
        ans = pycallargs(f, args_)
        pydel!(args_)
        ans
    else
        pycallargs(f)
    end
export pycall

"""
    pyeq(x, y)
    pyeq(Bool, x, y)

Equivalent to `x == y` in Python. The second form converts to `Bool`.
"""
pyeq(x, y) = pynew(errcheck(@autopy x y C.PyObject_RichCompare(getptr(x_), getptr(y_), C.Py_EQ)))
"""
    pyne(x, y)
    pyne(Bool, x, y)

Equivalent to `x != y` in Python. The second form converts to `Bool`.
"""
pyne(x, y) = pynew(errcheck(@autopy x y C.PyObject_RichCompare(getptr(x_), getptr(y_), C.Py_NE)))
"""
    pyle(x, y)
    pyle(Bool, x, y)

Equivalent to `x <= y` in Python. The second form converts to `Bool`.
"""
pyle(x, y) = pynew(errcheck(@autopy x y C.PyObject_RichCompare(getptr(x_), getptr(y_), C.Py_LE)))
"""
    pylt(x, y)
    pylt(Bool, x, y)

Equivalent to `x < y` in Python. The second form converts to `Bool`.
"""
pylt(x, y) = pynew(errcheck(@autopy x y C.PyObject_RichCompare(getptr(x_), getptr(y_), C.Py_LT)))
"""
    pyge(x, y)
    pyge(Bool, x, y)

Equivalent to `x >= y` in Python. The second form converts to `Bool`.
"""
pyge(x, y) = pynew(errcheck(@autopy x y C.PyObject_RichCompare(getptr(x_), getptr(y_), C.Py_GE)))
"""
    pygt(x, y)
    pygt(Bool, x, y)

Equivalent to `x > y` in Python. The second form converts to `Bool`.
"""
pygt(x, y) = pynew(errcheck(@autopy x y C.PyObject_RichCompare(getptr(x_), getptr(y_), C.Py_GT)))
pyeq(::Type{Bool}, x, y) = errcheck(@autopy x y C.PyObject_RichCompareBool(getptr(x_), getptr(y_), C.Py_EQ)) == 1
pyne(::Type{Bool}, x, y) = errcheck(@autopy x y C.PyObject_RichCompareBool(getptr(x_), getptr(y_), C.Py_NE)) == 1
pyle(::Type{Bool}, x, y) = errcheck(@autopy x y C.PyObject_RichCompareBool(getptr(x_), getptr(y_), C.Py_LE)) == 1
pylt(::Type{Bool}, x, y) = errcheck(@autopy x y C.PyObject_RichCompareBool(getptr(x_), getptr(y_), C.Py_LT)) == 1
pyge(::Type{Bool}, x, y) = errcheck(@autopy x y C.PyObject_RichCompareBool(getptr(x_), getptr(y_), C.Py_GE)) == 1
pygt(::Type{Bool}, x, y) = errcheck(@autopy x y C.PyObject_RichCompareBool(getptr(x_), getptr(y_), C.Py_GT)) == 1
export pyeq, pyne, pyle, pylt, pyge, pygt

pyconvert_rule_object(::Type{Py}, x) = pyconvert_return(Py(x))

"""
    pycontains(x, v)

Equivalent to `v in x` in Python.
"""
pycontains(x, v) = errcheck(@autopy x v C.PySequence_Contains(getptr(x_), getptr(v_))) == 1
export pycontains

"""
    pyin(v, x)

Equivalent to `v in x` in Python.
"""
pyin(v, x) = pycontains(x, v)
export pyin

pynotin(v, x) = !pyin(v, x)


text:
using Yao, LinearAlgebra, YaoBlocks

Ugate(nbit::Int, i::Int) = put(nbit, i=>chain(Rz(0), Rx(0), Rz(0)));
Ugate(i::Int) = put(i=>chain(Rz(0), Rx(0), Rz(0)));

function build_circuit()
    # circuit = chain(2, control(1, 2=>H), Ugate(1))
    circuit = chain(2, control(1, 2=>H), Ugate(1))
    return circuit
end

struct InitStateMatrix
    m1::Matrix{ComplexF64}
    m2::Matrix{ComplexF64}
end

function build_state(mats::InitStateMatrix)
    # return zero_state(2) |> chain(2, put(1=>chain(Rz(1), Rx(2))), put(2=>chain(Rz(1), Rx(2))))
    return zero_state(2) |> chain(2, put(2, (1)=>matblock(mats.m1)), put(2, (2)=>matblock(mats.m2)))
end

# @show probs(build_state())

function overlap(x::InitStateMatrix)
    return dot(apply!(build_state(x), chain(2, control(1, 2=>X), put(1=>H))) |> probs,[1, 1, 1, -1])
end

function cost(trainX, trainY, circuit, c)
    actual_out = [dot(i |> circuit |> probs, c) for i in map(build_state, trainX)]
    return sum((trainY-actual_out).^2)/length(trainY)
end

function gradient(circuit, trainX, trainY, delta)
    n = nqubits(circuit)
    grad = zeros(Float64, nparameters(circuit))
    params = parameters(circuit)

    count = 1
    for k in 1:nparameters(circuit)
        params[k] += delta
        dispatch!(circuit, params)
        cost_pos = cost(trainX, trainY, circuit, c)

        params[k] -= 2 * delta
        dispatch!(circuit, params)
        cost_neg = cost(trainX, trainY, circuit, c)

        params[k] += delta
        dispatch!(circuit, params)

        # @show (cost_pos - cost_neg)

        grad[count] = (cost_pos - cost_neg)/(2 * delta)
        count += 1
    end
    return grad
end

function generate_gates_set(n::Int64)
    rtn = []
    for i in 1:n
        for j in 1:n
            if i != j
                push!(rtn, control(n, i, j=>X))
            end
        end
    end

    for i in 1:n
        push!(rtn, Ugate(n, i))
    end

    return rtn
end

# 检测生成的circuit是否合法
function validate_circuit(circuit)
    last = circuit[1]
    for i in 2:length(circuit)
        now = circuit[i]
        if isa(last, ControlBlock) && last == now
            return false
        end
        if (!isa(last, ControlBlock)) && (!isa(now, ControlBlock)) && occupied_locs(last) == occupied_locs(now)
            return false
        end
        last = now
    end

    return true
end

function generate_init_circuit(n::Int64, d::Int64)
    gate_set = generate_gates_set(n)
    rtn = chain(n)
    while length(rtn) < d
        push!(rtn, rand(gate_set))
        if !validate_circuit(rtn)
            pop!(rtn)
        end
    end
    return rtn
end

function random_modify(circuit, gate_set)
    circuit2 = deepcopy(circuit)
    index = rand(1:length(circuit2))
    # @show index
    circuit2[index] = rand(gate_set)
    if validate_circuit(circuit2)
        return circuit2
    end
    return random_modify(circuit, gate_set)
end

using Flux.Optimise

c = [1, 1, 1, -1];
# circuit = build_circuit()
trainX = [InitStateMatrix(rand_unitary(2), rand_unitary(2)) for _ in 1:8]
trainY = map(overlap, trainX)
opt = Descent(3.0)

function train(circuit)
    history = Float64[]
    last_cost = now_cost = 0
    for i in 1:200
        last_cost = now_cost
        now_cost = cost(trainX, trainY, circuit, c)
        push!(history, now_cost)
        ps = parameters(circuit)
        # @show (i, ps, now_cost)
        if abs(now_cost - last_cost) < now_cost * .005 || now_cost < 1e-7
            # @show now_cost
            break
        end
        Optimise.update!(opt, ps, gradient(circuit, trainX, trainY, .001))
        popdispatch!(circuit, ps)
    end
    return (now_cost, history)
end

function find_algo(annealing_param, d::Int64)
    gate_set = generate_gates_set(2)
    circuit = generate_init_circuit(2, d)
    (best_cost, history) = train(circuit)
    count = 0
    circuit_backup = circuit
    while best_cost > 1e-7 && count < 1000
        count += 1;
        circuit_backup = deepcopy(circuit)
        circuit = random_modify(circuit, gate_set)
        (new_cost, history) = train(circuit)

        accept_prob = exp(-annealing_param * (new_cost - best_cost))
        if new_cost < best_cost
            @show (count, new_cost)
            best_cost = new_cost
            continue
        elseif rand() > accept_prob
            @show "fallback"
            circuit = circuit_backup
        end
        @show (count, accept_prob)
    end
    return circuit
end

@show find_algo(10.0, 2)
# (now_cost, history) = train(circuit);

# using Plots
#
# function plot_history(history)
#     pyplot() # Choose a backend
#     fig1 = plot(history; legend=nothing)
#     title!("training history")
#     xlabel!("steps"); ylabel!("cost")
# end

ext: jl


text:
typealias Stats Dict{AbstractString,Vector{UInt64}}

macro stats(name,expr)
    name = esc(name)
    expr = esc(expr)
    stats = esc(:stats)
    quote
        local gc_stats = Base.gc_num()
        local elapsedtime = Base.time_ns()
        $expr
        elapsedtime = Base.time_ns() - elapsedtime
        local diff = Base.GC_Diff(Base.gc_num(),gc_stats)
        if $name in keys($stats)
            $stats[$name] += [elapsedtime,diff.allocd,diff.total_time,Base.gc_alloc_count(diff)]
        else
            $stats[$name] = [elapsedtime,diff.allocd,diff.total_time,Base.gc_alloc_count(diff)]
        end
    end
end

function print_stats(stats::Stats)
    for k in sort(collect(keys(stats)), by=k->stats[k][1], rev=true)  # sort by elapsed time
        @printf "%-30s" k
        Base.time_print(stats[k]...)
    end
end

max_stars_repo_name: johncsnyder/Element


text:
using GLVisualize, GeometryTypes, Colors, GLAbstraction

if !isdefined(:runtests)
	window = glscreen()
	timesignal = bounce(linspace(0f0,1f0,360))
end

t = const_lift(*, timesignal, 20pi)
n = 50
const yrange = linspace(0.03, 0.3, n)
trange = linspace(0, 10pi, 200)

function gen_points(timesignal, y)
    x = sin(timesignal+(y*60*pi*y)+y)*y*5
    z = cos((timesignal+pi)+(y*60*pi*y)+y)*y*5
    Point3f0(x,y*60f0,z)
end
function gen_points(timesignal)
    Point3f0[gen_points(timesignal, y) for y in yrange]
end

positions = map(gen_points, t)
scale     = map(Vec3f0, linspace(0.05, 0.6, n))
primitive = centered(Sphere)
color     = map(RGB{Float32}, colormap("RdBu", n))
points 	  = visualize((primitive, positions), scale=scale, color=color)

_view(points, window)

if !isdefined(:runtests)
	renderloop(window)
end


text:
using Base.Test

using Rifraf

import Rifraf.random_seq
import Rifraf.sample_from_template
import Rifraf.sample_reference

srand(1)

@testset "sample" begin
    const seq_errors = ErrorModel(1.0, 1.0, 1.0, 0.0, 0.0)
    const ref_errors = ErrorModel(1.0, 0.0, 0.0, 1.0, 1.0)

    const params = Dict(:ref_error_rate => 0.05,
                        :ref_errors => ref_errors,
                        :error_rate => 0.01,
                        :alpha => 1.0,
                        :phred_scale => 5.0,
                        :actual_std => 3.0,
                        :reported_std => 1.0,
                        :seq_errors => seq_errors)

    template = random_seq(102)
    reference = sample_reference(template, 0.1, ref_errors)

    @testset "sample_from_template" begin
        template = random_seq(102)
        template_error_p = 0.05 * ones(102)
        seq, actual, reported = sample_from_template(template,
                                                     template_error_p,
                                                     params[:seq_errors],
                                                     params[:phred_scale],
                                                     params[:actual_std],
                                                     params[:reported_std])
    end

    @testset "test_sample()" begin
        nseqs = 3
        len = 900
        (ref, template, template_error_p,
         seqs, actual, phreds) = sample_sequences(nseqs, len; params...)
        @test mean(template_error_p) ≈ params[:error_rate] atol=0.1
    end

    @testset "test_sample_mixture" begin
        nseqs = (3, 3)
        len = 900
        n_diffs = 3
        (ref, template, template_error_p,
         seqs, actual, phreds) = sample_mixture(nseqs, len, n_diffs;
                                                params...)
        @test length(seqs) == sum(nseqs)
    end
end


text:
module libsdl2
	using CBinding
	using SDL2_jll

	c`$([
		Main.SYSROOT...,

		"-I$(SDL2_jll.artifact_dir)/include",
		"-I$(SDL2_jll.artifact_dir)/include/SDL2",
		"-L$(dirname(SDL2_jll.libsdl2_path))", "-lSDL2",
	])`

	const c"int8_t" = Int8
	const c"int16_t" = Int16
	const c"int32_t" = Int32
	const c"int64_t" = Int64
	const c"uint8_t" = UInt8
	const c"uint16_t" = UInt16
	const c"uint32_t" = UInt32
	const c"uint64_t" = UInt64
	const c"size_t" = Csize_t
	const c"wchar_t" = Cwchar_t

	# define them as Cvoid since they are usually used as opaque
	const c"FILE"    = Cvoid
	const c"va_list" = Cvoid

	c"""
		#include <SDL2/SDL.h>
	"""ji
end

module libsdl2_image
	using CBinding
	using ..libsdl2
	using SDL2_image_jll

	c`$([
		Main.SYSROOT...,

		"-I$(libsdl2.SDL2_jll.artifact_dir)/include",
		"-I$(libsdl2.SDL2_jll.artifact_dir)/include/SDL2",
		"-L$(dirname(libsdl2.SDL2_jll.libsdl2_path))", "-lSDL2",

		"-I$(SDL2_image_jll.artifact_dir)/include",
		"-I$(SDL2_image_jll.artifact_dir)/include/SDL2",
		"-L$(dirname(SDL2_image_jll.libsdl2_image_path))", "-lSDL2_image",
	])`

	const c"int8_t" = Int8
	const c"int16_t" = Int16
	const c"int32_t" = Int32
	const c"int64_t" = Int64
	const c"uint8_t" = UInt8
	const c"uint16_t" = UInt16
	const c"uint32_t" = UInt32
	const c"uint64_t" = UInt64
	const c"size_t" = Csize_t
	const c"wchar_t" = Cwchar_t

	# define them as Cvoid since they are usually used as opaque
	const c"FILE"    = Cvoid
	const c"va_list" = Cvoid

	c"""
		#include <SDL2/SDL.h>
	"""s
	c"""
		#include <SDL2/SDL_image.h>
	"""ji
end

module libsdl2_mixer
	using CBinding
	using ..libsdl2
	using SDL2_mixer_jll

	c`$([
		Main.SYSROOT...,

		"-I$(libsdl2.SDL2_jll.artifact_dir)/include",
		"-I$(libsdl2.SDL2_jll.artifact_dir)/include/SDL2",
		"-L$(dirname(libsdl2.SDL2_jll.libsdl2_path))", "-lSDL2",

		"-I$(SDL2_mixer_jll.artifact_dir)/include",
		"-I$(SDL2_mixer_jll.artifact_dir)/include/SDL2",
		"-L$(dirname(SDL2_mixer_jll.libsdl2_mixer_path))", "-lSDL2_mixer",
	])`

	const c"int8_t" = Int8
	const c"int16_t" = Int16
	const c"int32_t" = Int32
	const c"int64_t" = Int64
	const c"uint8_t" = UInt8
	const c"uint16_t" = UInt16
	const c"uint32_t" = UInt32
	const c"uint64_t" = UInt64
	const c"size_t" = Csize_t
	const c"wchar_t" = Cwchar_t

	# define them as Cvoid since they are usually used as opaque
	const c"FILE"    = Cvoid
	const c"va_list" = Cvoid

	c"""
		#include <SDL2/SDL.h>
	"""s
	c"""
		#include <SDL2/SDL_mixer.h>
	"""ji
end

module libsdl2_gfx
	using CBinding
	using ..libsdl2
	using SDL2_gfx_jll

	c`$([
		Main.SYSROOT...,

		"-I$(libsdl2.SDL2_jll.artifact_dir)/include",
		"-I$(libsdl2.SDL2_jll.artifact_dir)/include/SDL2",
		"-L$(dirname(libsdl2.SDL2_jll.libsdl2_path))", "-lSDL2",

		"-I$(SDL2_gfx_jll.artifact_dir)/include",
		"-I$(SDL2_gfx_jll.artifact_dir)/include/SDL2",
		"-L$(dirname(SDL2_gfx_jll.libsdl2_gfx_path))", "-lSDL2_gfx",
	])`

	const c"int8_t" = Int8
	const c"int16_t" = Int16
	const c"int32_t" = Int32
	const c"int64_t" = Int64
	const c"uint8_t" = UInt8
	const c"uint16_t" = UInt16
	const c"uint32_t" = UInt32
	const c"uint64_t" = UInt64
	const c"size_t" = Csize_t
	const c"wchar_t" = Cwchar_t

	# define them as Cvoid since they are usually used as opaque
	const c"FILE"    = Cvoid
	const c"va_list" = Cvoid

	c"""
		#include <SDL2/SDL.h>
	"""s
	c"""
		#include <SDL2/SDL2_framerate.h>
		#include <SDL2/SDL2_gfxPrimitives.h>
		#include <SDL2/SDL2_imageFilter.h>
		#include <SDL2/SDL2_rotozoom.h>
	"""ji
end

module libsdl2_ttf
	using CBinding
	using ..libsdl2
	using SDL2_ttf_jll

	c`$([
		Main.SYSROOT...,

		"-I$(libsdl2.SDL2_jll.artifact_dir)/include",
		"-I$(libsdl2.SDL2_jll.artifact_dir)/include/SDL2",
		"-L$(dirname(libsdl2.SDL2_jll.libsdl2_path))", "-lSDL2",

		"-I$(SDL2_ttf_jll.artifact_dir)/include",
		"-I$(SDL2_ttf_jll.artifact_dir)/include/SDL2",
		"-L$(dirname(SDL2_ttf_jll.libsdl2_ttf_path))", "-lSDL2_ttf",
	])`

	const c"int8_t" = Int8
	const c"int16_t" = Int16
	const c"int32_t" = Int32
	const c"int64_t" = Int64
	const c"uint8_t" = UInt8
	const c"uint16_t" = UInt16
	const c"uint32_t" = UInt32
	const c"uint64_t" = UInt64
	const c"size_t" = Csize_t
	const c"wchar_t" = Cwchar_t

	# define them as Cvoid since they are usually used as opaque
	const c"FILE"    = Cvoid
	const c"va_list" = Cvoid

	c"""
		#include <SDL2/SDL.h>
	"""s
	c"""
		#include <SDL2/SDL_ttf.h>
	"""ji
end


text:
using ComplexMixtures
using PDBTools

# Load pdb file of the system
system = readPDB("./simulation/equilibrated.pdb")

# The trajectory file can be downloaded from:
# https://drive.google.com/file/d/12TT5tblkFp1NtFOAQgjjGhmnYaXA8vQi/view?usp=sharing
trajectory_file = "../trajectories/traj_POPC.dcd"

# Select the atoms corresponding to glycerol and water
popc = select(system,"resname POPC")
water = select(system,"water")
ethanol = select(system,"resname ETOH")

# Set the complete membrane as the solute. We use nmols=1 here such
# that the membrane is considered a single solute in the calculation.
solute = Selection(popc,nmols=1)

# Compute water-POPC distribution and KB integral
solvent = Selection(water,natomspermol=3)
traj = Trajectory(trajectory_file,solute,solvent)

# We want to get reasonably converged KB integrals, which usually
# require large solute domains. Distribution functions converge
# rapidly (~10Angs or less), on the other side.
opt = ComplexMixtures.Options(dbulk=20.)

# Compute the MDDF
mddf_water_POPC = mddf(traj,opt)

# Save results for later use
save(mddf_water_POPC,"./results/mddf_water_POPC.json")

# Compute water-POPC distribution and KB integral
solvent = Selection(ethanol,natomspermol=9)
traj = Trajectory(trajectory_file,solute,solvent)
opt = ComplexMixtures.Options(dbulk=20.)
mddf_ethanol_POPC = mddf(traj,opt)

# Save results for later use
save(mddf_ethanol_POPC,"./results/mddf_ethanol_POPC.json")


text: export Spline, Bearing, bearing!, controlForX, yVal


text:
# computes the errors as outlined in the notes. Use these to generate
# moment conditions for estimation
function DSGEmoments(θ, data)
        # break out variables
        y = data[:,1]
        c = data[:,2]
        n = data[:,3]
        r = data[:,4]
        w = data[:,5]
        # break out params
        α = 0.33
        δ = 0.025
        β = θ[1]
        γ = θ[2]
        ρ_z = θ[3]
        σ_z = θ[4]
        ρ_η = θ[5]
        σ_η = θ[6]
        nss = θ[7]
        # recover ψ
        c1 = ((1.0/β + δ - 1.0)/α)^(1.0/(1.0-α))
        kss = nss/c1
        iss = δ*kss
        yss = kss^α * nss^(1-α)
        css = yss - iss
        ψ =  (css^(-γ)) * (1-α) * (kss^α) * (nss^(-α))
        # use MPL-MRS to get η, the preference shock
        η = log.(w) .- γ*log.(c) .- log.(ψ)
        u = (η[2:end]-η[1:end-1]*ρ_η)/σ_η
        e1 = η[1:end-1].*u
        e2 = u.^2.0 .- 1.0
        pref_shock = u
        # now the Euler eqn
        e3 = (100*((1.0 .+ r .- δ).*β.*(c.^(-γ)) .- lag(c,1).^(-γ)))[2:end]
        # recover K from MPK/MPL
        k = (α/(1.0-α))*n.*w./r
        # get z, the production shock, from the production function
        z = log.(y) - α*log.(k) - (1.0-α)*log.(n)
        u = (z[2:end]-z[1:end-1]*ρ_z)/σ_z
        e4 = z[1:end-1].*u
        e5 = u.^2.0 .- 1.0
        tech_shock = u
        # make moment conditions
        data = log.(data)[1:end-1,:]
        errors = [e1 e2 e3 e4 e5 pref_shock.*data tech_shock.*data data.*e3]
        return errors
end

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
# Our tests attempt translating strings to French, so set the LANGUAGE
# accordingly.
old_language = get(ENV, "LANGUAGE", nothing)
ENV["LANGUAGE"] = "fr"

using Gettext
using Base.Test
using Formatting

# Set up gettext
trdir = realpath(joinpath(Pkg.dir("Gettext"), "po"))
@test isfile(joinpath(trdir, "fr", "LC_MESSAGES", "sample.mo"))
bindtextdomain("sample", trdir)
textdomain("sample")

# Test basic macros
@test _"Hello, world!" == "Bonjour le monde !"
@test N_"Hello, world!" == "Hello, world!"

# Test ngettext
daystr(n) = format(ngettext("{1} day", "{1} days", n), n)
@test daystr(1) == "1 jour"
@test daystr(3) == "3 jours"

# Test dgettext and dngettext
@test dgettext("sample", "Hello, world!") == "Bonjour le monde !"
@test dngettext("sample", "{1} day", "{1} days", 1) == "{1} jour"

# Set the language back to normal.
if old_language != nothing
    ENV["LANGUAGE"] = old_language
else
    pop!(ENV, "LANGUAGE")
end


text:
@testset "[nl_write] Temp file handling" begin
    # Turn on debug mode so files persist
    old_debug = AmplNLWriter.debug
    AmplNLWriter.setdebug(true)

    filename = "test"
    filepath = joinpath(AmplNLWriter.solverdata_dir, filename)
    probfile = "$filepath.nl"
    solfile = "$filepath.sol"
    AmplNLWriter.clean_solverdata()

    @testset "all temp files deleted successfully" begin
        @test length(readdir(AmplNLWriter.solverdata_dir)) == 1
    end

    MOI = AmplNLWriter.MathOptInterface
    solver = AmplNLWriter.Optimizer(Ipopt.amplexe, filename=filename)
    x = MOI.add_variable(solver)
    MOI.add_constraint(
        solver,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
        MOI.GreaterThan(0.0))
    MOI.set(
        solver,
        MOI.ObjectiveFunction{MOI.SingleVariable}(),
        MOI.SingleVariable(x))
    MOI.set(solver, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(solver)

    @testset "temp files present after solve in debug mode" begin
        @test length(readdir(AmplNLWriter.solverdata_dir)) == 3
    end
    @testset "temp files used custom name" begin
        @test isfile(probfile) == true
        @test isfile(solfile) == true
    end

    # Reset debug mode and clean up
    AmplNLWriter.setdebug(old_debug)
    AmplNLWriter.clean_solverdata()
end


text:
random = rand(1:100)

println("Guess a number between 1 and 100.")

while true
    data = readLine()
    if data == ""
        println("Goodbye")
        break
    end
    guess = parse(Int64, data)
    if guess < random
        println("Too low")
    elseif guess > random
        println("Too high")
    else
        println("That's right!")
        break
end
# Might work


text:
using StructuralCausalModels

ProjDir = @__DIR__
cd(ProjDir)

d = OrderedDict(
  :b => :a,
  [:c, :d] => :b,
  :e => [:c, :d],
  :e => :f
);

d = OrderedDict(
  :b => :a,
  [:c, :d] => :b,
  :e => [:c, :d, :f]
);

dag = DAG("fig_3_1", d)

fname = scm_path("..", "examples", "Shipley", "fig3.1", "fig3.1.dot")
to_graphviz(dag, fname)
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

ggm = to_ggm(d)
display(ggm)

bs = basis_set(dag)
display(bs)

adjustment_sets(dag, :c, :e) |> display


text:
using LinearAlgebra
using CUDA.CUSPARSE, SparseArrays
using CUDA

@testset "sparse" begin
    n, m = 4, 4
    I = [1,2,3] |> cu
    J = [2,3,4] |> cu
    V = Float32[1,2,3] |> cu

    dense = rand(3,3) |> cu

    # check defaults
    @test sparse(I, J, V) isa CuSparseMatrixCSC
    @test sparse(dense) isa CuSparseMatrixCSC

    for (fmt, T) in  [(:coo, CuSparseMatrixCOO),
                      (:csc, CuSparseMatrixCSC),
                      (:csr, CuSparseMatrixCSR),
                      (:bsr, CuSparseMatrixBSR)
                     ]
        if fmt != :bsr # bsr not supported
            x = sparse(I, J, V; fmt=fmt)
            @test x isa T{Float32}
            @test size(x) == (3, 4)

            x = sparse(I, J, V, m, n; fmt=fmt)
            @test x isa T{Float32}
            @test  size(x) == (4, 4)
        end

        if fmt != :coo # dense to COO not implemented
            x = sparse(dense; fmt=fmt)
            @test x isa T{Float32}
            @test collect(x) == collect(dense)
        end
    end
end

@testset "unsorted sparse (CUDA.jl#1407)" begin
    I = [1, 1, 2, 3, 3, 4, 5, 4, 6, 4, 5, 6, 6, 6]
    J = [4, 6, 4, 5, 6, 6, 6, 1, 1, 2, 3, 3, 4, 5]

    # ensure we cover both the CUSPARSE-based and native COO row sort
    for typ in (Float16, Float32)
        A = sparse(I, J, ones(typ, length(I)), 6, 6)
        Agpu = sparse(I |> cu, J |> cu, ones(typ, length(I)) |> cu, 6, 6)
        @test Array(Agpu) == A
    end
end

@testset "CuSparseMatrix(::Diagonal)" begin
    X = Diagonal(rand(10))
    dX = cu(X)
    dY = CuSparseMatrixCSC{Float64, Int32}(dX)
    dZ = CuSparseMatrixCSR{Float64, Int32}(dX)
    @test SparseMatrixCSC(dY) ≈ SparseMatrixCSC(dZ)
    @test SparseMatrixCSC(CuSparseMatrixCSC(X)) ≈ SparseMatrixCSC(CuSparseMatrixCSR(X))
end


text:
#+foo_doc
#
# You can redefine the foo function:
#
# #+BEGIN_SRC julia
# # julia code example
# foo(x::Float64) = 1
# #+END_SRC
#
foo(x::Int) = 2*x


text:
# put unsuccessful code here
# just in case

# from imps.jl
# imps uninjectivity issue. tried to regulize it with finite power order
function pseudo_ovlp(psi::TensorMap{ComplexSpace, 2, 1}, phi1::TensorMap{ComplexSpace, 2, 1}, phi2::TensorMap{ComplexSpace, 2, 1}, power_order::Integer)
    ln_ovlp1 = log(ovlp(psi, phi1))
    ln_ovlp2 = log(ovlp(psi, phi2))
    ln_ovlp_tot = ln_ovlp1 + log1p(exp(power_order*(ln_ovlp2-ln_ovlp1))) / power_order
    return ln_ovlp_tot

end

function variational_truncate(phi::TensorMap{ComplexSpace, 2, 1}, chi::Integer)
    _, phi = left_canonical(phi)

    _, psi = iTEBD_truncate(phi, chi)
    psi_arr = toarray(psi)

    function _f(_arr)
        _psi = arr_to_TensorMap(_arr)
        up = norm(ovlp(_psi, phi))^2
        dn = ovlp(_psi, _psi)
        return -(up/dn) |> real
    end
    function _g!(_grad, _arr)
        _grad .= gradient(_f, _arr)[1]
    end

    res_f = optimize(_f, _g!, psi_arr, LBFGS(), Optim.Options(show_trace=false, iterations=100))
    psi = Optim.minimizer(res_f) |> arr_to_TensorMap

    return psi
end

# from cmps.jl
# periodic uniform cMPS compression issue. tried to construct the projector into the cMPS tangent space. not working.
function C_matrix(D::TensorMap{ComplexSpace, 1, 1}, L::Real)
    Cmat = similar(D)
    w_vec = diag(D.data)
    N = length(w_vec)
    for ix in 1:N
        for iy in 1:N
            if w_vec[ix] ≈ w_vec[iy]
                Cmat[ix, iy] = exp(w_vec[ix])
            else
                Cmat[ix, iy] = (exp(w_vec[ix] * L) - exp(w_vec[iy] * L)) / (w_vec[ix] - w_vec[iy])
            end
        end
    end
    return Cmat
end

function delta_tensor(chi2::Integer)
    δ = TensorMap(zeros, ComplexF64, ℂ^chi2*ℂ^chi2, ℂ^chi2)
    for ix in 1:chi2
        δ[ix, ix, ix] = 1
    end
    return δ
end

function company_tensor(psi::cmps)
    A = convert_to_array(psi)
    chi = get_chi(psi)
    A[:, 1, :] = Matrix{ComplexF64}(I, chi, chi)
    return convert_to_tensormap(A, 2)
end

"""
    gram_matrix(psi::cmps, L::Real)

    calculate the gram matrix.
"""
function gram_matrix(psi::cmps, L::Real)
    chi = get_chi(psi)

    K = K_mat(psi, psi)
    K = permute(K, (2, 1), (4, 3))
    D, VR = eig(K)
    VL = inv(VR)' # left eigenvalues

    VR = permute(VR, (1,), (3, 2))
    VL = permute(VL, (1,), (3, 2))

    δ = delta_tensor(chi^2)
    C = C_matrix(D, L)
    A = company_tensor(psi)
    Iph = id(ℂ^(get_d(psi)+1))
    Iph[1,1] = 0
    expD = exp(D*L)

    @tensor gram[-1, -2, -3; -4, -5, -6] :=
        δ[3, 2, 1] * VL'[1, 7, -5] * A'[8, 7, -6] * VR[-4, 4, 8] * VR[9, 2, -3] * A[10, -2, 9] * VL'[5, -1, 10] * δ'[4, 5, 6] * C[6, 3] +
        VL'[2, -1, -5] * Iph[-2, -6] * VR[-4, 1, -3] * expD[1, 2]
    return gram
end

"""
    tangent_proj(phi::cmps, psi::cmps, L::Real)

    project cmps `phi` into the tangent space of cmps `psi`. `L` is the length of the cmps.
"""
function tangent_proj(phi::cmps, psi::cmps, L::Real, tol::Real)

    gram_psi = gram_matrix(psi, L)
    gram_psi_inv = pinv(gram_psi; rtol=tol)
    gram_psi_inv = permute(gram_psi_inv, (6,2,3), (4,5,1))

    K = transf_mat(psi, phi)
    KT = transf_mat_T(psi, phi)
    V0 = TensorMap(rand, ComplexF64, ℂ^get_chi(phi), ℂ^get_chi(psi))

    _, VR = eigsolve(K, V0, 1, :LR)
    VR = VR[1]
    _, VL = eigsolve(KT, V0, 1, :LR)
    VL = VL[1]

    phi_tMap = company_tensor(phi)

    @tensor B_proj[-1, -2; -3] := phi_tMap[1, 5, 2] * VL'[3, 1] * VR[2, 4] * gram_psi_inv[4, -1, -2, 3, 5, -3]

    return B_proj
end


text:
function assert_enoughvaluesforepochs(vname::String, v::Vector, epochs::Int)
   if length(v) < epochs
      error("Not enough `$vname` (vector of length $(length(v))) " .*
         "for training epochs ($epochs)")
   end
end

function assert_initoptimizers(optimizer::AbstractOptimizer,
      optimizers::Vector{<:AbstractOptimizer}, bm::BM,
      learningrates::Vector{Float64}, sdlearningrates::Vector{Float64},
      epochs::Int
      ) where {BM<:AbstractBM}

   if isempty(optimizers)
      if optimizer === NoOptimizer()
         # default optimization algorithm
         optimizers = map(
               i -> defaultoptimizer(bm, learningrates[i], sdlearningrates[i]),
               1:epochs)
      else
         optimizers = fill(initialized(optimizer, bm), epochs)
      end
   else
      assert_enoughvaluesforepochs("optimizers", optimizers, epochs)
      initalizedoptimizers = Vector(undef, epochs)
      for i in 1:epochs
         initalizedoptimizers[i] = initialized(optimizers[i], bm)
      end
      optimizers = converttomostspecifictype(initalizedoptimizers)
   end

   optimizers
end

function defaultoptimizer(rbm::R,
      learningrate::Float64, sdlearningrate::Float64) where {R<: AbstractRBM}

   loglikelihoodoptimizer(rbm; learningrate = learningrate,
         sdlearningrate = sdlearningrate)
end

function defaultoptimizer(dbm::MultimodalDBM, learningrate::Float64,
      sdlearningrate::Float64)

   StackedOptimizer(map(
         rbm -> defaultoptimizer(rbm, learningrate, sdlearningrate),
         dbm))
end

"""
    fitrbm(x; ...)
Fits an RBM model to the data set `x`, using Stochastic Gradient Descent (SGD)
with Contrastive Divergence (CD), and returns it.

# Optional keyword arguments (ordered by importance):
* `rbmtype`: the type of the RBM that is to be trained
   This must be a subtype of `AbstractRBM` and defaults to `BernoulliRBM`.
* `nhidden`: number of hidden units for the returned RBM
* `epochs`: number of training epochs
* `learningrate`/`learningrates`: The learning rate for the weights and biases
   can be specified as single value, used throughout all epochs, or as a vector
   of `learningrates` that contains a value for each epoch. Defaults to 0.005.
* `batchsize`: number of samples that are used for making one step in the
   stochastic gradient descent optimizer algorithm. Default is 1.
* `pcd`: indicating whether Persistent Contrastive Divergence (PCD) is to
   be used (true, default) or simple CD that initializes the Gibbs Chain with
   the training sample (false)
* `cdsteps`: number of Gibbs sampling steps for (persistent)
   contrastive divergence, defaults to 1
* `monitoring`: a function that is executed after each training epoch.
   It takes an RBM and the epoch as arguments.
   See also `monitored_fitrbm` for another way of monitoring.
* `categories`: only relevant if `rbmtype = Softmax0BernoulliRBM`.
   The number of categories as `Int`, if all variables have the same number
   of categories, or as `Vector{Int}` that contains the number of categories
   of the i'th categorical variable in the i'th entry.
* `upfactor`, `downfactor`: If this function is used for pretraining a part of
   a DBM, it is necessary to multiply the weights of the RBM with factors.
* `sdlearningrate`/`sdlearningrates`: learning rate(s) for the
   standard deviation if training a `GaussianBernoulliRBM` or
   `GaussianBernoulliRBM2`. Ignored for other types of RBMs.
   It usually must be much smaller than the learning rates for
   the weights. By default it is 0.0, which means that the standard deviation
   is not learned.
* `startrbm`: start training with the parameters of the given RBM.
   If this argument is specified, `nhidden` and `rbmtype` are ignored.
* `optimizer`/`optimizers`: an object of type `AbstractOptimizer` or a vector of
   them for each epoch. If specified, the optimization is performed as implemented
   by the given optimizer type. By default, the `LoglikelihoodOptimizer`
   with the `learningrate`/`learningrates` and `sdlearningrate`/`sdlearningrates`
   is used. For other types of optimizers, the learning rates must be specified
   in the `optimizer`. For more information on how to write your own optimizer,
   see `AbstractOptimizer`.
See also: `monitored_fitrbm` for a convenient monitoring of the training.
"""
function fitrbm(x::Matrix{Float64};
      nhidden::Int = size(x,2),
      epochs::Int = 10,
      upfactor::Float64 = 1.0,
      downfactor::Float64 = 1.0,
      learningrate::Float64 = 0.005,
      learningrates::Vector{Float64} = fill(learningrate, epochs),
      pcd::Bool = true,
      cdsteps::Int = 1,
      batchsize::Int = 1,
      rbmtype::DataType = BernoulliRBM,
      startrbm::AbstractRBM = NoRBM(),
      monitoring::Function = emptyfunc,

      # this argument is only relevant for Softmax0BernoulliRBMs:
      categories::Union{Int, Vector{Int}} = 0,

      # these arguments are only relevant for GaussianBernoulliRBMs:
      sdlearningrate::Float64 = 0.0,
      sdlearningrates::Vector{Float64} = fill(sdlearningrate, epochs),
      sdinitfactor::Float64 = 0.0,

      optimizer::AbstractOptimizer = NoOptimizer(),
      optimizers::Vector{<:AbstractOptimizer} = Vector{AbstractOptimizer}(),
      sampler::AbstractSampler = (cdsteps == 1 ? NoSampler() : GibbsSampler(cdsteps - 1)))

   if startrbm === NoRBM()
      if rbmtype == Softmax0BernoulliRBM
         if categories == 0
            error("Argument `categories` required for `Softmax0BernoulliRBM`s.")
         end
         rbm = initsoftmax0rbm(x, nhidden, categories)
      else
         rbm = initrbm(x, nhidden, rbmtype)
      end
   else
      rbm = deepcopy(startrbm)
      nhidden = nhiddennodes(startrbm)
   end

   assert_enoughvaluesforepochs("learningrates", learningrates, epochs)
   assert_enoughvaluesforepochs("sdlearningrates", sdlearningrates, epochs)

   if sdinitfactor > 0 &&
         (rbmtype == GaussianBernoulliRBM || rbmtype == GaussianBernoulliRBM2)
      rbm.sd .*= sdinitfactor
   end

   optimizers = assert_initoptimizers(optimizer, optimizers, rbm,
         learningrates, sdlearningrates, epochs)

   if pcd
      chainstate = rand(batchsize, nhidden)
   else
      chainstate = Matrix{Float64}(undef, 0, 0)
   end

   # allocate space for trainrbm!
   nvisible = size(x, 2)
   h = Matrix{Float64}(undef, batchsize, nhidden)
   hmodel = Matrix{Float64}(undef, batchsize, nhidden)
   vmodel = Matrix{Float64}(undef, batchsize, nvisible)

   for epoch = 1:epochs

      # Train RBM on data set
      trainrbm!(rbm, x, cdsteps = cdsteps, chainstate = chainstate,
            upfactor = upfactor, downfactor = downfactor,
            learningrate = learningrates[epoch],
            sdlearningrate = sdlearningrates[epoch],
            optimizer = optimizers[epoch],
            sampler = sampler,
            batchsize = batchsize,
            h = h, hmodel = hmodel, vmodel = vmodel)

      # Evaluation of learning after each training epoch
      monitoring(rbm, epoch)
   end

   rbm
end

"""
    initrbm(x, nhidden)
    initrbm(x, nhidden, rbmtype)
Creates a RBM with `nhidden` hidden units and initalizes its weights for
training on dataset `x`.
`rbmtype` can be a subtype of `AbstractRBM`, default is `BernoulliRBM`.
"""
function initrbm(x::Array{Float64,2}, nhidden::Int,
      rbmtype::DataType = BernoulliRBM)

   nsamples, nvisible = size(x)
   weights = randn(nvisible, nhidden)/sqrt(nvisible)
   hidbias = zeros(nhidden)

   if rbmtype == BernoulliRBM
      visbias = initvisiblebias(x)
      return BernoulliRBM(weights, visbias, hidbias)

   elseif rbmtype == GaussianBernoulliRBM
      visbias = vec(mean(x, dims = 1))
      sd = vec(std(x, dims = 1))
      #sd = fill(0.05, length(vec(std(x, 1))))
      GaussianBernoulliRBM(weights, visbias, hidbias, sd)

   elseif rbmtype == GaussianBernoulliRBM2
      visbias = vec(mean(x, dims = 1))
      #sd = fill(0.05, length(vec(std(x, 1))))
      sd = vec(std(x, dims = 1))
      weights .*= sd
      return GaussianBernoulliRBM2(weights, visbias, hidbias, sd)

   elseif rbmtype == BernoulliGaussianRBM
      visbias = initvisiblebias(x)
      return BernoulliGaussianRBM(weights, visbias, hidbias)

   elseif rbmtype == Binomial2BernoulliRBM
      visbias = initvisiblebias(x/2)
      return Binomial2BernoulliRBM(weights, visbias, hidbias)

   elseif rbmtype == Softmax0BernoulliRBM
      error("For initialization of `Softmax0BernoulliRBM`s, please use `initsoftmax0rbm`.")

   else
      error(string("Datatype for RBM is unsupported: ", rbmtype))
   end
end

function initsoftmax0rbm(x::Matrix{Float64}, nhidden::Int, ncategories::Int)
   initsoftmax0rbm(x, nhidden, fill(ncategories, div(size(x, 2), ncategories - 1)))
end

function initsoftmax0rbm(x::Matrix{Float64}, nhidden::Int, nscategories::Vector{Int})
   if size(x, 2) != sum(nscategories) - length(nscategories)
      error("Number of variables ($(size(x, 2))) is not consisten with the " .*
            "given categories ($(nscategories)).")
   end
   rbm = initrbm(x, nhidden)
   Softmax0BernoulliRBM(rbm.weights, rbm.visbias, rbm.hidbias, nscategories)
end

"""
    initvisiblebias(x)
Returns sensible initial values for the visible bias for training an RBM
on the data set `x`.
"""
function initvisiblebias(x::Array{Float64,2})
   nvisible = size(x,2)
   initbias = zeros(nvisible)
   for j=1:nvisible
      empprob = mean(x[:,j])
      if empprob > 0
         initbias[j] = log(empprob/(1-empprob))
      end
   end
   initbias
end

"""
    randombatchmasks(nsamples, batchsize)
Returns BitArray-Sets for the sample indices when training on a dataset with
`nsamples` samples using minibatches of size `batchsize`.
"""
function randombatchmasks(nsamples::Int, batchsize::Int)
   if batchsize > nsamples
      error("Batchsize ($batchsize) exceeds number of samples ($nsamples).")
   end

   nfullbatches, nremainingbatch = divrem(nsamples, batchsize)
   batchsizes = fill(batchsize, nfullbatches)
   if nremainingbatch > 0
      push!(batchsizes, nremainingbatch)
   end
   randomsampleindices = randperm(nsamples)
   map(batchrange -> randomsampleindices[batchrange],
         ranges(batchsizes))
end

function sdupdateterm(gbrbm::GaussianBernoulliRBM,
         v::Matrix{Float64}, h::Matrix{Float64})
   vec(mean((v .- gbrbm.visbias').^2 ./ (gbrbm.sd' .^ 3) -
         h * gbrbm.weights' .* (v ./ (gbrbm.sd' .^ 2)), dims = 1))
end

"""
    trainrbm!(rbm, x)
Trains the given `rbm` for one epoch using the data set `x`.
(See also function `fitrbm`.)

# Optional keyword arguments:
* `learningrate`, `cdsteps`, `sdlearningrate`, `upfactor`, `downfactor`,
   `optimizer`:
   See documentation of function `fitrbm`.
* `chainstate`: a matrix for holding the states of the RBM's hidden nodes. If
   it is specified, PCD is used.
"""
function trainrbm!(rbm::AbstractRBM, x::Array{Float64,2};
      chainstate::Matrix{Float64} = Matrix{Float64}(undef, 0, 0),
      upfactor::Float64 = 1.0,
      downfactor::Float64 = 1.0,
      learningrate::Float64 = 0.005,
      cdsteps::Int = 1,
      batchsize::Int = 1,
      sdlearningrate::Float64 = 0.0,
      optimizer::AbstractOptimizer = LoglikelihoodOptimizer(rbm;
            learningrate = learningrate, sdlearningrate = sdlearningrate),
      sampler::AbstractSampler = (cdsteps == 1 ? NoSampler() : GibbsSampler(cdsteps - 1)),

      # write-only arguments for reusing allocated space:
      v::Matrix{Float64} = Matrix{Float64}(undef, batchsize, length(rbm.visbias)),
      h::Matrix{Float64} = Matrix{Float64}(undef, batchsize, length(rbm.hidbias)),
      hmodel::Matrix{Float64} = Matrix{Float64}(undef, batchsize, length(rbm.hidbias)),
      vmodel::Matrix{Float64} = Matrix{Float64}(undef, batchsize, length(rbm.visbias)))

   nsamples = size(x, 1)

   # perform PCD if a chain state is provided as parameter
   pcd = !isempty(chainstate)

   batchmasks = randombatchmasks(nsamples, batchsize)
   nbatchmasks = length(batchmasks)

   normalbatchsize = true

   # go through all samples or thorugh all batches of samples
   for batchindex in eachindex(batchmasks)
      batchmask = batchmasks[batchindex]

      normalbatchsize = (batchindex < nbatchmasks || nsamples % batchsize == 0)

      if normalbatchsize
         v .= view(x, batchmask, :)
      else
         v = x[batchmask, :]
         thisbatchsize = nsamples % batchsize
         h = Matrix{Float64}(undef, thisbatchsize, nhiddennodes(rbm))
         if !pcd
            vmodel = Matrix{Float64}(undef, size(v))
            hmodel = Matrix{Float64}(undef, size(h))
         end # in case of pcd, vmodel and hmodel are not downsized
      end

      # Calculate potential induced by visible nodes, used for update
      hiddenpotential!(h, rbm, v, upfactor)

      # In case of CD, start Gibbs chain with the hidden state induced by the
      # sample. In case of PCD, start Gibbs chain with
      # previous state of the Gibbs chain.
      if pcd
         hmodel = chainstate # note: state of chain will be visible by the caller
      else
         copyto!(hmodel, h)
      end
      samplehiddenpotential!(hmodel, rbm)

      # Additional sampling steps, may be customized
      sample!(vmodel, hmodel, sampler, rbm, upfactor, downfactor)

      # Do not sample in last step to avoid unnecessary sampling noise
      visiblepotential!(vmodel, rbm, hmodel, downfactor)
      hiddenpotential!(hmodel, rbm, vmodel, upfactor)

      if !normalbatchsize && pcd
         # remove additional samples
         # that are unnecessary for computing the gradient
         vmodel = vmodel[1:thisbatchsize, :]
         hmodel = hmodel[1:thisbatchsize, :]
      end

      computegradient!(optimizer, v, vmodel, h, hmodel, rbm)
      updateparameters!(rbm, optimizer)
   end

   rbm
end


text:
module Universes

using Setfield

using ...PCG.Types
using ...PCG.Storages: get_node, apply_changes!, Storage, StorageIndex, StorageNode
using ...PCG.Topologies: Topology, TopologyIndex, coordinates
using ...PCG.Storages
using PCG.ArraysCaches

export Element, AreaElements, Universe, finish_recording!, AreaCache, construct_current_element, construct_new_element, enabled, storage_index, storage_size, complete_turn!

abstract type AbstractUniverse end

function storage_index end
function storage_size end

struct Element{U<:AbstractUniverse, TI<:TopologyIndex, SI<:StorageIndex, P<:AbstractProperties}
    enabled::Bool # TODO: is other way exists to emulate Union{Element, Nothing} for not isbits element?
    universe::U
    topology_index::TI
    storage_index::SI
    properties::P
end

function Types.isenabled(element::Element)
    return element.enabled
end

function Types.disable(element::E) where E
    # TODO: simplify?
    return E(false, element.universe, element.topology_index, element.storage_index, element.properties)
end

Base.convert(::Type{Bool}, element::Element) = isenabled(element)

const AreaElements{U, TI, SI, N} = Vector{Element{U, TI, SI, N}}

Base.convert(::Type{Bool}, elements::AreaElements) = any(convert(Bool, element) for element in elements)

const AreaCache{U, TI, SI, P} = ArraysCache{AreaElements{U, TI, SI, P}}

mutable struct Universe{S<:Storage, T<:Topology} <: AbstractUniverse
    storage::S
    topology::T

    turn::Turn
    recorders::Recorders

    cache::Union{AreaCache, Nothing}
end

Universe(storage::S, topology::T, recorders::Recorders) where {S, T} = Universe(storage,
                                                                              topology,
                                                                              0,
                                                                              recorders,
                                                                              nothing)

# TODO: lazy calculation of fields "storage_index" & "node"?
function Element(universe::Universe, i::TI) where {TI<:TopologyIndex}
    index = storage_index(universe.storage, universe.topology, i)
    node = get_node(universe.storage, index)
    return Element(true, universe, i, index, node.current)
end

# TODO refactor construct_*_element functions into one

# TODO: why constructor call by template variable does not see short constructor?
function Types.construct_current_element(::Type{Element{U, TI, SI, P}}, universe::U, i::TI)::Element{U, TI, SI, P} where {U<:AbstractUniverse, TI<:TopologyIndex, SI<:StorageIndex, P<:AbstractProperties}
    index = storage_index(universe.storage, universe.topology, i)
    node = get_node(universe.storage, index)
    return Element{U, TI, SI, P}(true, universe, i, index, node.current)
end

function Types.construct_new_element(::Type{Element{U, TI, SI, P}}, universe::U, i::TI)::Element{U, TI, SI, P} where {U<:AbstractUniverse, TI<:TopologyIndex, SI<:StorageIndex, P<:AbstractProperties}
    index = storage_index(universe.storage, universe.topology, i)
    node = get_node(universe.storage, index)
    return Element{U, TI, SI, P}(true, universe, i, index, node.new)
end

function Types.reserve_area!(element::Element{U, TI, SI, P}, size::Int64)::AreaElements{U, TI, SI, P} where {U, TI, SI, P}
    return reserve!(element.universe.cache::AreaCache{U, TI, SI, P}, size)
end

function recorders(universe::Universe)
    return universe.recorders
end

function complete_turn!(universe::Universe)

    apply_changes!(universe.storage, universe.turn)

    universe.turn += 1

    record_state!(universe)
end

function record_state!(universe::Universe)
    for recorder in recorders(universe)
        record_state!(universe, recorder)
    end
end

function finish_recording!(universe::Universe)
    for recorder in recorders(universe)
        finish_recording!(recorder)
    end
end

function (universe::Universe)(callback::Function; complete_turn::Bool=true, turns::Turn=1)
    for turn in 1:turns
        for i in coordinates(universe.topology)
            element = Element(universe, i)
            callback(element)
            release_all!(universe.cache)
        end

        if complete_turn
            complete_turn!(universe)
        end
    end
end

function ArraysCaches.release_all!(universe::Universe)
    release_all!(universe.cache)
end

# TODO: construct new node on base of already changed new node?
function Base.:<<(element::Element, new_values::NamedTuple)
    new_node = StorageNode(element.properties,
                           setproperties(element.properties, new_values),
                           element.universe.turn)

    set_node!(element.universe.storage, element.storage_index, new_node)
end

end

max_stars_repo_name: Tiendil/pcg-julia

max_issues_repo_name: Tiendil/pcg-julia

max_forks_repo_name: Tiendil/pcg-julia


text:
using LoopVectorization, Test

function multiassign!(y, x)
  @assert length(y) + 3 == length(x)
  @inbounds for i ∈ eachindex(y)
    x₁, ((x₂, x₃), (x₄, x₅)) = x[i], (sincos(x[i+1]), (x[i+2], x[i+3]))
    y[i] = x₁ * x₄ - x₂ * x₃
  end
  y
end
multiassign(x) = multiassign!(similar(x, length(x) - 3), x)
function multiassign_turbo!(y, x)
  @assert length(y) + 3 == length(x)
  @turbo for i ∈ eachindex(y)
    x₁, ((x₂, x₃), (x₄, x₅)) = x[i], (sincos(x[i+1]), (x[i+2], x[i+3]))
    y[i] = x₁ * x₄ - x₂ * x₃
  end
  y
end
multiassign_turbo(x) = multiassign_turbo!(similar(x, length(x) - 3), x)

multistorefunc(x) = exp(x), sincos(x)
function multistore!(x, y, z, a)
  @inbounds for i ∈ eachindex(x, y, z, a)
    x[i], (y[i], z[i]) = multistorefunc(a[i])
  end
end
function multistore_turbo!(x, y, z, a)
  @turbo for i ∈ eachindex(x, y, z, a)
    x[i], (y[i], z[i]) = multistorefunc(a[i])
  end
end

@testset "Multiple assignments" begin
  @show @__LINE__
  x = rand(111)
  @test multiassign(x) ≈ multiassign_turbo(x)
  a0 = similar(x)
  b0 = similar(x)
  c0 = similar(x)
  a1 = similar(x)
  b1 = similar(x)
  c1 = similar(x)
  multistore!(a0, b0, c0, x)
  multistore_turbo!(a1, b1, c1, x)
  @test a0 ≈ a1
  @test b0 ≈ b1
  @test c0 ≈ c1
end


text:
"""
```
pseudo_measurement(m::Model1002{T}, TTT::Matrix{T}, RRR::Matrix{T},
    CCC::Vector{T}; reg::Int = 1) where {T<:AbstractFloat}
```

Assign pseudo-measurement equation (a linear combination of states):

```
x_t = ZZ_pseudo*s_t + DD_pseudo
```
"""
function pseudo_measurement(m::Model1002{T},
                            TTT::Matrix{T},
                            RRR::Matrix{T},
                            CCC::Vector{T};
                            reg::Int = 1) where {T<:AbstractFloat}

    endo      = m.endogenous_states
    endo_addl = m.endogenous_states_augmented
    pseudo    = m.pseudo_observables

    # Initialize pseudo ZZ and DD matrices
    _n_states = n_states_augmented(m)
    _n_pseudo = n_pseudo_observables(m)

    ZZ_pseudo = zeros(_n_pseudo, _n_states)
    DD_pseudo = zeros(_n_pseudo)

    # Set parameters
    for para in m.parameters
        if !isempty(para.regimes)
            ModelConstructors.toggle_regime!(para, reg)
        end
    end

    # Handle integrated series
    no_integ_inds = inds_states_no_integ_series(m)

    if haskey(m.endogenous_states, :pgap_t)
        no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:pgap_t]])
    end
    if haskey(m.endogenous_states, :ygap_t)
        no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:ygap_t]])
    end

    if haskey(get_settings(m), :integrated_series) || haskey(m.endogenous_states, :pgap_t) ||
        haskey(m.endogenous_states, :ygap_t)
        TTT = @view TTT[no_integ_inds, no_integ_inds]
    end

    # Compute TTT^10, used for Expected10YearRateGap, Expected10YearRate, and Expected10YearNaturalRate
    TTT10 = (1/40)*((UniformScaling(1.) - TTT)\(TTT - TTT^41))

    if get_setting(m, :add_laborproductivity_measurement)
        # Construct pseudo-obs from integrated states first
        ZZ_pseudo[pseudo[:laborproductivity], endo[:y_t]] = 1.
        ZZ_pseudo[pseudo[:laborproductivity], endo[:L_t]] = -1.
        # ZZ_pseudo[pseudo[:laborproductivity], endo_addl[:cum_z_t]] = 1.
        DD_pseudo[pseudo[:laborproductivity]] = 100. * log(m[:ystar] / m[:Lstar])

        # Remove integrated states (e.g. states w/unit roots)
        # RRR and CCC aren't used, so we don't do anything with them
    end

    if get_setting(m, :add_nominalgdp_level)
        ZZ_pseudo[pseudo[:NominalGDPLevel], endo_addl[:cum_y_t]]     = 1.
        ZZ_pseudo[pseudo[:NominalGDPLevel], endo_addl[:cum_z_t]]     = 1.
        ZZ_pseudo[pseudo[:NominalGDPLevel], endo_addl[:cum_e_gdp_t]] = 1.
        ZZ_pseudo[pseudo[:NominalGDPLevel], endo_addl[:cum_π_t]]     = 1.
    end

    if get_setting(m, :add_nominalgdp_growth)
        ZZ_pseudo[pseudo[:NominalGDPGrowth], endo[:y_t]]           = 1.
        ZZ_pseudo[pseudo[:NominalGDPGrowth], endo_addl[:y_t1]]     = -1.
        ZZ_pseudo[pseudo[:NominalGDPGrowth], endo_addl[:z_t]]      = 1.
        ZZ_pseudo[pseudo[:NominalGDPGrowth], endo_addl[:e_gdp_t]]  = 1.
        ZZ_pseudo[pseudo[:NominalGDPGrowth], endo_addl[:e_gdp_t1]] = -m[:me_level]
        ZZ_pseudo[pseudo[:NominalGDPGrowth], endo[:π_t]]           = 1.
        DD_pseudo[pseudo[:NominalGDPGrowth]]                       = 100. * (exp(m[:z_star] - 1.) + (m[:π_star] - 1.))
    end

    if get_setting(m, :add_cumulative)
        ZZ_pseudo[pseudo[:AccumOutputGap], endo_addl[:cum_y_t]]   = 1.
        ZZ_pseudo[pseudo[:AccumOutputGap], endo_addl[:cum_y_f_t]] = -1.

        ZZ_pseudo[pseudo[:GDPLevel], endo_addl[:cum_y_t]]     = 1.
        ZZ_pseudo[pseudo[:GDPLevel], endo_addl[:cum_z_t]]     = 1.
        ZZ_pseudo[pseudo[:GDPLevel], endo_addl[:cum_e_gdp_t]] = 1.

        ZZ_pseudo[pseudo[:FlexibleGDPLevel], endo_addl[:cum_y_f_t]] = 1.
        ZZ_pseudo[pseudo[:FlexibleGDPLevel], endo_addl[:cum_z_t]]   = 1.

        ZZ_pseudo[pseudo[:ConsumptionLevel], endo_addl[:cum_c_t]]     = 1.
        ZZ_pseudo[pseudo[:ConsumptionLevel], endo_addl[:cum_z_t]]     = 1.

        ZZ_pseudo[pseudo[:FlexibleConsumptionLevel], endo_addl[:cum_c_f_t]] = 1.
        ZZ_pseudo[pseudo[:FlexibleConsumptionLevel], endo_addl[:cum_z_t]]   = 1.

        ZZ_pseudo[pseudo[:InvestmentLevel], endo_addl[:cum_i_t]]     = 1.
        ZZ_pseudo[pseudo[:InvestmentLevel], endo_addl[:cum_z_t]]     = 1.

        ZZ_pseudo[pseudo[:FlexibleInvestmentLevel], endo_addl[:cum_i_f_t]] = 1.
        ZZ_pseudo[pseudo[:FlexibleInvestmentLevel], endo_addl[:cum_z_t]]   = 1.
    end

    if get_setting(m, :add_flexible_price_growth)
        ZZ_pseudo[pseudo[:FlexibleGDPGrowth], endo[:y_f_t]]       = 1.
        ZZ_pseudo[pseudo[:FlexibleGDPGrowth], endo_addl[:y_f_t1]] = -1.
        ZZ_pseudo[pseudo[:FlexibleGDPGrowth], endo[:z_t]]         = 1.
        DD_pseudo[pseudo[:FlexibleGDPGrowth]]                     = 100. * (exp(m[:z_star]) - 1.)

        ZZ_pseudo[pseudo[:FlexibleInvestmentGrowth], endo[:i_f_t]]       = 1.
        ZZ_pseudo[pseudo[:FlexibleInvestmentGrowth], endo_addl[:i_f_t1]] = -1.
        ZZ_pseudo[pseudo[:FlexibleInvestmentGrowth], endo[:z_t]]         = 1.
        DD_pseudo[pseudo[:FlexibleInvestmentGrowth]]                     = 100. * (exp(m[:z_star]) - 1.)

        ZZ_pseudo[pseudo[:FlexibleConsumptionGrowth], endo[:c_f_t]]       = 1.
        ZZ_pseudo[pseudo[:FlexibleConsumptionGrowth], endo_addl[:c_f_t1]] = -1.
        ZZ_pseudo[pseudo[:FlexibleConsumptionGrowth], endo[:z_t]]         = 1.
        DD_pseudo[pseudo[:FlexibleConsumptionGrowth]]                     = 100. * (exp(m[:z_star]) - 1.)
    end

    ##########################################################
    ## PSEUDO-OBSERVABLE EQUATIONS
    ##########################################################

    ## Output
    ZZ_pseudo[pseudo[:y_t],endo[:y_t]] = 1.

    ## Flexible Output
    ZZ_pseudo[pseudo[:y_f_t],endo[:y_f_t]] = 1.

    ## Pseudo GDP Growth
    if haskey(m.settings, :add_pseudo_gdp)
        if get_setting(m, :add_pseudo_gdp) && subspec(m) in ["ss59", "ss60", "ss61"]
            ZZ_pseudo[pseudo[:PseudoGDP], endo[:y_t]]          = 1.0
            ZZ_pseudo[pseudo[:PseudoGDP], endo_addl[:y_t1]]     = -1.0
            ZZ_pseudo[pseudo[:PseudoGDP], endo[:z_t]]          = 1.0
            ZZ_pseudo[pseudo[:PseudoGDP], endo_addl[:e_gdp_t]]  = 1.0
            ZZ_pseudo[pseudo[:PseudoGDP], endo_addl[:e_gdp_t1]] = -m[:me_level]
            DD_pseudo[pseudo[:PseudoGDP]]                      = 100*(exp(m[:z_star])-1)
        end
    end

    ## Natural Rate
    ZZ_pseudo[pseudo[:NaturalRate],endo[:r_f_t]] = 1.
    DD_pseudo[pseudo[:NaturalRate]]              = 100.0*(m[:rstar]-1.0)

    ## π_t
    ZZ_pseudo[pseudo[:π_t],endo[:π_t]] = 1.
    DD_pseudo[pseudo[:π_t]]            = 100*(m[:π_star]-1);

    ## Output Gap
    ZZ_pseudo[pseudo[:OutputGap],endo[:y_t]] = 1;
    ZZ_pseudo[pseudo[:OutputGap],endo[:y_f_t]] = -1;

    ## Ex Ante Real Rate
    ZZ_pseudo[pseudo[:ExAnteRealRate],endo[:R_t]]  = 1;
    ZZ_pseudo[pseudo[:ExAnteRealRate],endo[:Eπ_t]] = -1;
    DD_pseudo[pseudo[:ExAnteRealRate]]             = m[:Rstarn] - 100*(m[:π_star]-1);

    ## Long Run Inflation
    ZZ_pseudo[pseudo[:LongRunInflation],endo[:π_star_t]] = 1.
    DD_pseudo[pseudo[:LongRunInflation]]                 = 100. *(m[:π_star]-1.)

    ## Marginal Cost
    ZZ_pseudo[pseudo[:MarginalCost],endo[:mc_t]] = 1.

    ## Wages
    ZZ_pseudo[pseudo[:Wages],endo[:w_t]] = 1.

    ## Flexible Wages
    ZZ_pseudo[pseudo[:FlexibleWages],endo[:w_f_t]] = 1.

    # ## b Wages
    # ZZ_pseudo[pseudo[:b_t],endo[:b_t]] = 1.

    ## Hours
    ZZ_pseudo[pseudo[:Hours],endo[:L_t]] = 1.

    ## Flexible Hours
    ZZ_pseudo[pseudo[:FlexibleHours],endo[:L_f_t]] = 1.

    ## z_t
    ZZ_pseudo[pseudo[:z_t], endo[:z_t]] = 1.

    ## Expected 10-Year Rate Gap
    ZZ_pseudo[pseudo[:Expected10YearRateGap], no_integ_inds] = TTT10[endo[:R_t], :] - TTT10[endo[:r_f_t], :] - TTT10[endo[:Eπ_t], :]

    ## Nominal FFR
    ZZ_pseudo[pseudo[:NominalFFR], endo[:R_t]] = 1.
    DD_pseudo[pseudo[:NominalFFR]] = m[:Rstarn]

    ## Expected 10-Year Interest Rate
    ZZ_pseudo[pseudo[:Expected10YearRate], no_integ_inds] = TTT10[endo[:R_t], :]
    DD_pseudo[pseudo[:Expected10YearRate]]    = m[:Rstarn]

    ## Expected 10-Year Natural Rate
    ZZ_pseudo[pseudo[:Expected10YearNaturalRate], no_integ_inds] = TTT10[endo[:r_f_t], :] + TTT10[endo[:Eπ_t], :]
    DD_pseudo[pseudo[:Expected10YearNaturalRate]]    = m[:Rstarn]

    ## Expected Nominal Natural Rate
    ZZ_pseudo[pseudo[:ExpectedNominalNaturalRate], endo[:r_f_t]] = 1.
    ZZ_pseudo[pseudo[:ExpectedNominalNaturalRate], endo[:Eπ_t]]  = 1.
    DD_pseudo[pseudo[:ExpectedNominalNaturalRate]]               = m[:Rstarn]

    ## Nominal Rate Gap
    ZZ_pseudo[pseudo[:NominalRateGap], endo[:R_t]]   = 1.
    ZZ_pseudo[pseudo[:NominalRateGap], endo[:r_f_t]] = -1.
    ZZ_pseudo[pseudo[:NominalRateGap], endo[:Eπ_t]]  = -1.

    ## Labor Productivity Growth
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo[:y_t]]           = 1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:y_t1]]     = -1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo[:z_t]]           = 1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:e_gdp_t]]  = 1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:e_gdp_t1]] = -m[:me_level]
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo[:L_t]]           = -1
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:L_t1]]     = 1.
    DD_pseudo[pseudo[:LaborProductivityGrowth]]                       = 100*(exp(m[:z_star]) - 1)

    ## u_t
    ZZ_pseudo[pseudo[:u_t], endo[:u_t]] = 1.

    ## Nominal Wage Growth
    if haskey(m.settings, :add_NominalWageGrowth)
        if get_setting(m, :add_NominalWageGrowth)
            ZZ_pseudo[pseudo[:NominalWageGrowth],endo[:w_t]] = 1.
            ZZ_pseudo[pseudo[:NominalWageGrowth],endo_addl[:w_t1]] = -1.
            ZZ_pseudo[pseudo[:NominalWageGrowth],endo[:z_t]] = 1.
            ZZ_pseudo[pseudo[:NominalWageGrowth],endo[:π_t]] = 1.
            DD_pseudo[pseudo[:NominalWageGrowth]]            = 100*(m[:π_star]-1) + 100*(exp(m[:z_star])-1)
        end
    end

    # ## i_f_t
    # ZZ_pseudo[pseudo[:i_f_t], endo[:i_f_t]] = 1.

    # ## R_t
    # ZZ_pseudo[pseudo[:R_t], endo[:R_t]] = 1.
    # DD_pseudo[pseudo[:R_t]] = 100.0*(m[:rstar]-1.0)

    # ## c_f_t
    # ZZ_pseudo[pseudo[:c_f_t], endo[:c_f_t]] = 1.

    # ## c_t
    # ZZ_pseudo[pseudo[:c_t], endo[:c_t]] = 1.

    # ## qk_f_t
    # ZZ_pseudo[pseudo[:qk_f_t], endo[:qk_f_t]] = 1.

    # ## k_f_t
    # ZZ_pseudo[pseudo[:k_f_t], endo[:k_f_t]] = 1.

    # ## r_f_t
    # ZZ_pseudo[pseudo[:r_f_t], endo[:r_f_t]] = 1.

    # ## kbar_f_t
    # ZZ_pseudo[pseudo[:kbar_f_t], endo[:kbar_f_t]] = 1.

    # ## u_f_t
    # ZZ_pseudo[pseudo[:u_f_t], endo[:u_f_t]] = 1.

    # ## rk_f_t
    # ZZ_pseudo[pseudo[:rk_f_t], endo[:rk_f_t]] = 1.

    # ## w_f_t
    # ZZ_pseudo[pseudo[:w_f_t], endo[:w_f_t]] = 1.

    # ## L_f_t
    # ZZ_pseudo[pseudo[:L_f_t], endo[:L_f_t]] = 1.

    # ## rktil_f_t
    # ZZ_pseudo[pseudo[:rktil_f_t], endo[:rktil_f_t]] = 1.

    # ## n_f_t
    # ZZ_pseudo[pseudo[:n_f_t], endo[:n_f_t]] = 1.

    ## labor share
    if haskey(m.settings, :add_laborshare_measurement)
        if get_setting(m, :add_laborshare_measurement)
            ZZ_pseudo[pseudo[:laborshare_t], endo[:w_t]] = 1.
            ZZ_pseudo[pseudo[:laborshare_t], endo[:L_t]] = 1.
            ZZ_pseudo[pseudo[:laborshare_t], endo[:y_t]] = -1.
            DD_pseudo[pseudo[:laborshare_t]] = 100. * log(m[:wstar] * m[:Lstar] / m[:ystar])
        end
    end

    ## Fundamental inflation related pseudo-obs
    if subspec(m) in ["ss13", "ss14", "ss15", "ss16", "ss17", "ss18", "ss19"]
        # Compute coefficient on Sinf
        betabar = exp((1-m[:σ_c] ) * m[:z_star]) * m[:β]
        κ = ((1 - m[:ζ_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))*
             (1 - m[:ζ_p]))/(m[:ζ_p]*((m[:Φ]- 1)*m[:ϵ_p] + 1))/
        (1 + m[:ι_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
        κcoef = κ * (1 + m[:ι_p] * betabar)

        ZZ_pseudo[pseudo[:Sinf_t], endo_addl[:Sinf_t]] = 1.
        ZZ_pseudo[pseudo[:Sinf_w_coef_t], endo_addl[:Sinf_t]] = κcoef
        DD_pseudo[pseudo[:ι_p]] = m[:ι_p]
        ZZ_pseudo[pseudo[:πtil_t], endo_addl[:πtil_t]] = 1.
        DD_pseudo[pseudo[:πtil_t]] = 100 * (m[:π_star] - 1)
        ZZ_pseudo[pseudo[:e_tfp_t], endo_addl[:e_tfp_t]] = 1.
        if subspec(m) in ["ss14", "ss15", "ss16", "ss18", "ss19"]
            ZZ_pseudo[pseudo[:e_tfp_t1], endo_addl[:e_tfp_t1]] = 1.
        end
    end

    ## Exogenous processes
    if subspec(m) == "ss12"
        to_add = [:g_t, :b_t, :μ_t, :z_t, :λ_f_t, :λ_w_t, :rm_t, :σ_ω_t, :μ_e_t,
                  :γ_t, :π_star_t]
        to_add_addl = [:e_lr_t, :e_tfp_t, :e_gdpdef_t, :e_corepce_t, :e_gdp_t, :e_gdi_t]
        for i in to_add
            ZZ_pseudo[pseudo[i], endo[i]] = 1.
        end
        for i in to_add_addl
            ZZ_pseudo[pseudo[i], endo_addl[i]] = 1.
        end
    end

    if haskey(m.settings, :add_covid_pseudoobs)
        if get_setting(m, :add_covid_pseudoobs) && subspec(m) in ["ss59", "ss60", "ss61"]
            ZZ_pseudo[pseudo[:ziid], endo[:ziid_t]] = 1.
            ZZ_pseudo[pseudo[:varphiiid], endo[:φ_t]] = 1.
            ZZ_pseudo[pseudo[:biidc], endo[:biidc_t]] = 1.
        end
    end

    if haskey(m.settings, :add_ztil)
        if get_setting(m, :add_ztil)
            ZZ_pseudo[pseudo[:ztil], endo[:ztil_t]] = 1.
        end
    end

    if haskey(m.settings, :add_zp)
        if get_setting(m, :add_zp)
            ZZ_pseudo[pseudo[:zp], endo[:zp_t]] = 1.
        end
    end

    if haskey(m.settings, :add_pgap)
        if get_setting(m, :add_pgap)
            ZZ_pseudo[pseudo[:pgap], endo[:pgap_t]] = 1.
        end
    end

    if haskey(m.settings, :add_ygap)
        if get_setting(m, :add_ygap)
            ZZ_pseudo[pseudo[:ygap], endo[:ygap_t]] = 1.
        end
    end

    if haskey(m.settings, :add_urhat)
        if get_setting(m, :add_urhat)
            ZZ_pseudo[pseudo[:urhat], endo[:L_t]] = -0.7366220966444673
            DD_pseudo[pseudo[:urhat]] = -32.257825316004364 + -0.7366220966444673 * m[:Lmean]
        end
    end

    if haskey(m.settings, :add_rw)
        if get_setting(m, :add_rw)
            ZZ_pseudo[pseudo[:rw], endo[:rw_t]]     = 1.
            ZZ_pseudo[pseudo[:Rref], endo[:Rref_t]] = 1.
            DD_pseudo[pseudo[:Rref]]                = m[:Rstarn]
        end
    end

    for para in m.parameters
        if !isempty(para.regimes)
            ModelConstructors.toggle_regime!(para, 1)
        end
    end

    return PseudoMeasurement(ZZ_pseudo, DD_pseudo)
end

function pseudo_measurement(m::Model1002{T},
                            TTTs::Vector{Matrix{T}},
                            RRRs::Vector{Matrix{T}},
                            CCCs::Vector{Vector{T}}) where {T<:AbstractFloat}

    endo      = m.endogenous_states
    endo_addl = m.endogenous_states_augmented
    pseudo    = m.pseudo_observables

    _n_states = n_states_augmented(m)
    _n_pseudo = n_pseudo_observables(m)

    n_reg = length(TTTs)

    ZZ_pseudos = Vector{Matrix{Float64}}(undef, n_reg)
    DD_pseudos = Vector{Vector{Float64}}(undef, n_reg)

    for reg in 1:n_reg
        # Initialize pseudo ZZ and DD matrices
        ZZ_pseudos[reg] = zeros(_n_pseudo, _n_states)
        DD_pseudos[reg] = zeros(_n_pseudo)

        for para in m.parameters
            if !isempty(para.regimes)
                ModelConstructors.toggle_regime!(para, reg)
            end
        end

        no_integ_inds = inds_states_no_integ_series(m)
        if haskey(m.endogenous_states, :pgap_t)
            no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:pgap_t]])
        end
        if haskey(m.endogenous_states, :ygap_t)
            no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:ygap_t]])
        end

        if haskey(get_settings(m), :integrated_series) || haskey(m.endogenous_states, :pgap_t) || haskey(m.endogenous_states, :ygap_t)
            TTT = @view TTTs[reg][no_integ_inds, no_integ_inds]
        else
            TTT = TTTs[reg]
        end

        # Compute TTT^10, used for Expected10YearRateGap, Expected10YearRate, and Expected10YearNaturalRate
        TTT10 = (1/40)*((UniformScaling(1.) - TTT)\(TTT - TTT^41))

        if get_setting(m, :add_laborproductivity_measurement)
            # Construct pseudo-obs from integrated states first
            ZZ_pseudos[reg][pseudo[:laborproductivity], endo[:y_t]] = 1.
            ZZ_pseudos[reg][pseudo[:laborproductivity], endo[:L_t]] = -1.
            # ZZ_pseudos[reg][pseudo[:laborproductivity], endo_addl[:cum_z_t]] = 1.
            DD_pseudos[reg][pseudo[:laborproductivity]] = 100. * log(m[:ystar] / m[:Lstar])

            # Remove integrated states (e.g. states w/unit roots)
            # RRR and CCC aren't used, so we don't do anything with them
        end

        if get_setting(m, :add_nominalgdp_level)
            ZZ_pseudos[reg][pseudo[:NominalGDPLevel], endo_addl[:cum_y_t]]     = 1.
            ZZ_pseudos[reg][pseudo[:NominalGDPLevel], endo_addl[:cum_z_t]]     = 1.
            ZZ_pseudos[reg][pseudo[:NominalGDPLevel], endo_addl[:cum_e_gdp_t]] = 1.
            ZZ_pseudos[reg][pseudo[:NominalGDPLevel], endo_addl[:cum_π_t]]     = 1.
        end

        if get_setting(m, :add_nominalgdp_growth)
            ZZ_pseudos[reg][pseudo[:NominalGDPGrowth], endo[:y_t]]           = 1.
            ZZ_pseudos[reg][pseudo[:NominalGDPGrowth], endo_addl[:y_t1]]     = -1.
            ZZ_pseudos[reg][pseudo[:NominalGDPGrowth], endo_addl[:z_t]]      = 1.
            ZZ_pseudos[reg][pseudo[:NominalGDPGrowth], endo_addl[:e_gdp_t]]  = 1.
            ZZ_pseudos[reg][pseudo[:NominalGDPGrowth], endo_addl[:e_gdp_t1]] = -m[:me_level]
            ZZ_pseudos[reg][pseudo[:NominalGDPGrowth], endo[:π_t]]           = 1.
            DD_pseudos[reg][pseudo[:NominalGDPGrowth]]                       = 100. * (exp(m[:z_star] - 1.) + (m[:π_star] - 1.))
        end

        if get_setting(m, :add_cumulative)
            ZZ_pseudos[reg][pseudo[:AccumOutputGap], endo_addl[:cum_y_t]]   = 1.
            ZZ_pseudos[reg][pseudo[:AccumOutputGap], endo_addl[:cum_y_f_t]] = -1.

            ZZ_pseudos[reg][pseudo[:GDPLevel], endo_addl[:cum_y_t]]     = 1.
            ZZ_pseudos[reg][pseudo[:GDPLevel], endo_addl[:cum_z_t]]     = 1.
            ZZ_pseudos[reg][pseudo[:GDPLevel], endo_addl[:cum_e_gdp_t]] = 1.

            ZZ_pseudos[reg][pseudo[:FlexibleGDPLevel], endo_addl[:cum_y_f_t]] = 1.
            ZZ_pseudos[reg][pseudo[:FlexibleGDPLevel], endo_addl[:cum_z_t]]   = 1.

            ZZ_pseudos[reg][pseudo[:ConsumptionLevel], endo_addl[:cum_c_t]]     = 1.
            ZZ_pseudos[reg][pseudo[:ConsumptionLevel], endo_addl[:cum_z_t]]     = 1.

            ZZ_pseudos[reg][pseudo[:FlexibleConsumptionLevel], endo_addl[:cum_c_f_t]] = 1.
            ZZ_pseudos[reg][pseudo[:FlexibleConsumptionLevel], endo_addl[:cum_z_t]]   = 1.

            ZZ_pseudos[reg][pseudo[:InvestmentLevel], endo_addl[:cum_i_t]]     = 1.
            ZZ_pseudos[reg][pseudo[:InvestmentLevel], endo_addl[:cum_z_t]]     = 1.

            ZZ_pseudos[reg][pseudo[:FlexibleInvestmentLevel], endo_addl[:cum_i_f_t]] = 1.
            ZZ_pseudos[reg][pseudo[:FlexibleInvestmentLevel], endo_addl[:cum_z_t]]   = 1.
        end

        if get_setting(m, :add_flexible_price_growth)
            ZZ_pseudos[reg][pseudo[:FlexibleGDPGrowth], endo[:y_f_t]]       = 1.
            ZZ_pseudos[reg][pseudo[:FlexibleGDPGrowth], endo_addl[:y_f_t1]] = -1.
            ZZ_pseudos[reg][pseudo[:FlexibleGDPGrowth], endo[:z_t]]         = 1.
            DD_pseudos[reg][pseudo[:FlexibleGDPGrowth]]                     = 100. * (exp(m[:z_star]) - 1.)

            ZZ_pseudos[reg][pseudo[:FlexibleInvestmentGrowth], endo[:i_f_t]]       = 1.
            ZZ_pseudos[reg][pseudo[:FlexibleInvestmentGrowth], endo_addl[:i_f_t1]] = -1.
            ZZ_pseudos[reg][pseudo[:FlexibleInvestmentGrowth], endo[:z_t]]         = 1.
            DD_pseudos[reg][pseudo[:FlexibleInvestmentGrowth]]                     = 100. * (exp(m[:z_star]) - 1.)

            ZZ_pseudos[reg][pseudo[:FlexibleConsumptionGrowth], endo[:c_f_t]]       = 1.
            ZZ_pseudos[reg][pseudo[:FlexibleConsumptionGrowth], endo_addl[:c_f_t1]] = -1.
            ZZ_pseudos[reg][pseudo[:FlexibleConsumptionGrowth], endo[:z_t]]         = 1.
            DD_pseudos[reg][pseudo[:FlexibleConsumptionGrowth]]                     = 100. * (exp(m[:z_star]) - 1.)
        end

        ##########################################################
        ## PSEUDO-OBSERVABLE EQUATIONS
        ##########################################################

        ## Output
        ZZ_pseudos[reg][pseudo[:y_t],endo[:y_t]] = 1.

        ## Flexible Output
        ZZ_pseudos[reg][pseudo[:y_f_t],endo[:y_f_t]] = 1.

        ## Pseudo GDP Growth
        if haskey(m.settings, :add_pseudo_gdp)
            if get_setting(m, :add_pseudo_gdp) && subspec(m) in ["ss59", "ss60", "ss61"]
                ZZ_pseudos[reg][pseudo[:PseudoGDP], endo[:y_t]]          = 1.0
                ZZ_pseudos[reg][pseudo[:PseudoGDP], endo_addl[:y_t1]]     = -1.0
                ZZ_pseudos[reg][pseudo[:PseudoGDP], endo[:z_t]]          = 1.0
                ZZ_pseudos[reg][pseudo[:PseudoGDP], endo_addl[:e_gdp_t]]  = 1.0
                ZZ_pseudos[reg][pseudo[:PseudoGDP], endo_addl[:e_gdp_t1]] = -m[:me_level]
                DD_pseudos[reg][pseudo[:PseudoGDP]]                      = 100*(exp(m[:z_star])-1)
            end
        end
        ## Natural Rate
        ZZ_pseudos[reg][pseudo[:NaturalRate],endo[:r_f_t]] = 1.
        DD_pseudos[reg][pseudo[:NaturalRate]]              = 100.0*(m[:rstar]-1.0)

        ## π_t
        ZZ_pseudos[reg][pseudo[:π_t],endo[:π_t]] = 1.
        DD_pseudos[reg][pseudo[:π_t]]            = 100*(m[:π_star]-1);

        ## Output Gap
        ZZ_pseudos[reg][pseudo[:OutputGap],endo[:y_t]] = 1;
        ZZ_pseudos[reg][pseudo[:OutputGap],endo[:y_f_t]] = -1;

        ## Ex Ante Real Rate
        ZZ_pseudos[reg][pseudo[:ExAnteRealRate],endo[:R_t]]  = 1;
        ZZ_pseudos[reg][pseudo[:ExAnteRealRate],endo[:Eπ_t]] = -1;
        DD_pseudos[reg][pseudo[:ExAnteRealRate]]             = m[:Rstarn] - 100*(m[:π_star]-1);

        ## Long Run Inflation
        ZZ_pseudos[reg][pseudo[:LongRunInflation],endo[:π_star_t]] = 1.
        DD_pseudos[reg][pseudo[:LongRunInflation]]                 = 100. *(m[:π_star]-1.)

        ## Marginal Cost
        ZZ_pseudos[reg][pseudo[:MarginalCost],endo[:mc_t]] = 1.

        ## Wages
        ZZ_pseudos[reg][pseudo[:Wages],endo[:w_t]] = 1.

        # ## b Wages
        # ZZ_pseudos[reg][pseudo[:b_t],endo[:b_t]] = 1.

        ## Flexible Wages
        ZZ_pseudos[reg][pseudo[:FlexibleWages],endo[:w_f_t]] = 1.

        ## Hours
        ZZ_pseudos[reg][pseudo[:Hours],endo[:L_t]] = 1.

        ## Flexible Hours
        ZZ_pseudos[reg][pseudo[:FlexibleHours],endo[:L_f_t]] = 1.

        ## z_t
        ZZ_pseudos[reg][pseudo[:z_t], endo[:z_t]] = 1.

        ## Expected 10-Year Rate Gap
        ZZ_pseudos[reg][pseudo[:Expected10YearRateGap], no_integ_inds] = TTT10[endo[:R_t], :] - TTT10[endo[:r_f_t], :] - TTT10[endo[:Eπ_t], :]

        ## Nominal FFR
        ZZ_pseudos[reg][pseudo[:NominalFFR], endo[:R_t]] = 1.
        DD_pseudos[reg][pseudo[:NominalFFR]] = m[:Rstarn]

        ## Expected 10-Year Interest Rate
        ZZ_pseudos[reg][pseudo[:Expected10YearRate], no_integ_inds] = TTT10[endo[:R_t], :]
        DD_pseudos[reg][pseudo[:Expected10YearRate]]    = m[:Rstarn]

        ## Expected 10-Year Natural Rate
        ZZ_pseudos[reg][pseudo[:Expected10YearNaturalRate], no_integ_inds] = TTT10[endo[:r_f_t], :] + TTT10[endo[:Eπ_t], :]
        DD_pseudos[reg][pseudo[:Expected10YearNaturalRate]]    = m[:Rstarn]

        ## Expected Nominal Natural Rate
        ZZ_pseudos[reg][pseudo[:ExpectedNominalNaturalRate], endo[:r_f_t]] = 1.
        ZZ_pseudos[reg][pseudo[:ExpectedNominalNaturalRate], endo[:Eπ_t]]  = 1.
        DD_pseudos[reg][pseudo[:ExpectedNominalNaturalRate]]               = m[:Rstarn]

        ## Nominal Rate Gap
        ZZ_pseudos[reg][pseudo[:NominalRateGap], endo[:R_t]]   = 1.
        ZZ_pseudos[reg][pseudo[:NominalRateGap], endo[:r_f_t]] = -1.
        ZZ_pseudos[reg][pseudo[:NominalRateGap], endo[:Eπ_t]]  = -1.

        ## Labor Productivity Growth
        ZZ_pseudos[reg][pseudo[:LaborProductivityGrowth], endo[:y_t]]           = 1.
        ZZ_pseudos[reg][pseudo[:LaborProductivityGrowth], endo_addl[:y_t1]]     = -1.
        ZZ_pseudos[reg][pseudo[:LaborProductivityGrowth], endo[:z_t]]           = 1.
        ZZ_pseudos[reg][pseudo[:LaborProductivityGrowth], endo_addl[:e_gdp_t]]  = 1.
        ZZ_pseudos[reg][pseudo[:LaborProductivityGrowth], endo_addl[:e_gdp_t1]] = -m[:me_level]
        ZZ_pseudos[reg][pseudo[:LaborProductivityGrowth], endo[:L_t]]           = -1
        ZZ_pseudos[reg][pseudo[:LaborProductivityGrowth], endo_addl[:L_t1]]     = 1.
        DD_pseudos[reg][pseudo[:LaborProductivityGrowth]]                       = 100*(exp(m[:z_star]) - 1)

        ## u_t
        ZZ_pseudos[reg][pseudo[:u_t], endo[:u_t]] = 1.

        ## Nominal Wage Growth
        if haskey(m.settings, :add_NominalWageGrowth)
            if get_setting(m, :add_NominalWageGrowth)
                ZZ_pseudos[reg][pseudo[:NominalWageGrowth],endo[:w_t]] = 1.
                ZZ_pseudos[reg][pseudo[:NominalWageGrowth],endo_addl[:w_t1]] = -1.
                ZZ_pseudos[reg][pseudo[:NominalWageGrowth],endo[:z_t]] = 1.
                ZZ_pseudos[reg][pseudo[:NominalWageGrowth],endo[:π_t]] = 1.
                DD_pseudos[reg][pseudo[:NominalWageGrowth]]            = 100*(m[:π_star]-1) + 100*(exp(m[:z_star])-1)
            end
        end

        # ## i_f_t
        # ZZ_pseudos[reg][pseudo[:i_f_t], endo[:i_f_t]] = 1.

        # ## R_t
        # ZZ_pseudos[reg][pseudo[:R_t], endo[:R_t]] = 1.
        # DD_pseudos[reg][pseudo[:R_t]] = 100.0*(m[:rstar]-1.0)

        # ## c_f_t
        # ZZ_pseudos[reg][pseudo[:c_f_t], endo[:c_f_t]] = 1.

        # ## c_t
        # ZZ_pseudos[reg][pseudo[:c_t], endo[:c_t]] = 1.

        # ## r_f_t
        # ZZ_pseudos[reg][pseudo[:r_f_t], endo[:r_f_t]] = 1.

        # ## qk_f_t
        # ZZ_pseudos[reg][pseudo[:qk_f_t], endo[:qk_f_t]] = 1.

        # ## k_f_t
        # ZZ_pseudos[reg][pseudo[:k_f_t], endo[:k_f_t]] = 1.

        # ## kbar_f_t
        # ZZ_pseudos[reg][pseudo[:kbar_f_t], endo[:kbar_f_t]] = 1.

        # ## u_f_t
        # ZZ_pseudos[reg][pseudo[:u_f_t], endo[:u_f_t]] = 1.

        # ## rk_f_t
        # ZZ_pseudos[reg][pseudo[:rk_f_t], endo[:rk_f_t]] = 1.

        # ## w_f_t
        # ZZ_pseudos[reg][pseudo[:w_f_t], endo[:w_f_t]] = 1.

        # ## L_f_t
        # ZZ_pseudos[reg][pseudo[:L_f_t], endo[:L_f_t]] = 1.

        # ## rktil_f_t
        # ZZ_pseudos[reg][pseudo[:rktil_f_t], endo[:rktil_f_t]] = 1.

        # ## n_f_t
        # ZZ_pseudos[reg][pseudo[:n_f_t], endo[:n_f_t]] = 1.

        ## labor share
        if haskey(m.settings, :add_laborshare_measurement)
            if get_setting(m, :add_laborshare_measurement)
                ZZ_pseudos[reg][pseudo[:laborshare_t], endo[:w_t]] = 1.
                ZZ_pseudos[reg][pseudo[:laborshare_t], endo[:L_t]] = 1.
                ZZ_pseudos[reg][pseudo[:laborshare_t], endo[:y_t]] = -1.
                DD_pseudos[reg][pseudo[:laborshare_t]] = 100. * log(m[:wstar] * m[:Lstar] / m[:ystar])
            end
        end

        ## Labor Productivity Growth, no measurement error
        if haskey(m.settings, :add_laborproductivitygrowth_nome_measurement)
            if get_setting(m, :add_laborproductivitygrowth_nome_measurement)
                ZZ_pseudos[reg][pseudo[:LaborProductivityGrowthNoME], endo[:y_t]]       = 1.
                ZZ_pseudos[reg][pseudo[:LaborProductivityGrowthNoME], endo_addl[:y_t1]] = -1.
                ZZ_pseudos[reg][pseudo[:LaborProductivityGrowthNoME], endo[:z_t]]       = 1.
                ZZ_pseudos[reg][pseudo[:LaborProductivityGrowthNoME], endo[:L_t]]       = -1
                ZZ_pseudos[reg][pseudo[:LaborProductivityGrowthNoME], endo_addl[:L_t1]] = 1.
                DD_pseudos[reg][pseudo[:LaborProductivityGrowthNoME]]                   = 100*(exp(m[:z_star]) - 1)
            end
        end

        if haskey(m.settings, :add_covid_pseudoobs)
            if get_setting(m, :add_covid_pseudoobs) && subspec(m) in ["ss59", "ss60", "ss61"]
                ZZ_pseudos[reg][pseudo[:ziid], endo[:ziid_t]] = 1.
                ZZ_pseudos[reg][pseudo[:varphiiid], endo[:φ_t]] = 1.
                ZZ_pseudos[reg][pseudo[:biidc], endo[:biidc_t]] = 1.
            end
        end
        if haskey(m.settings, :add_ztil)
            if get_setting(m, :add_ztil)
                ZZ_pseudos[reg][pseudo[:ztil], endo[:ztil_t]] = 1.
            end
        end
        if haskey(m.settings, :add_zp)
            if get_setting(m, :add_zp)
                ZZ_pseudos[reg][pseudo[:zp], endo[:zp_t]] = 1.
            end
        end

        if haskey(m.settings, :add_ztil)
            if get_setting(m, :add_ztil)
                ZZ_pseudos[reg][pseudo[:ztil], endo[:ztil_t]] = 1.
            end
        end
        if haskey(m.settings, :add_zp)
            if get_setting(m, :add_zp)
                ZZ_pseudos[reg][pseudo[:zp], endo[:zp_t]] = 1.
            end
        end

        ## Labor Productivity Growth, no measurement error
        if haskey(m.settings, :add_Epi_t_measurement)
            if get_setting(m, :add_Epi_t_measurement)
                ZZ_pseudos[reg][pseudo[:Epi_t], endo[:Eπ_t]]       = 1.
                DD_pseudos[reg][pseudo[:Epi_t]]                    = 100. * (m[:π_star] - 1.)
            end
        end
        if haskey(m.settings, :add_pgap)
            if get_setting(m, :add_pgap)
                ZZ_pseudos[reg][pseudo[:pgap], endo[:pgap_t]] = 1.
            end
        end
        if haskey(m.settings, :add_ygap)
            if get_setting(m, :add_ygap)
                ZZ_pseudos[reg][pseudo[:ygap], endo[:ygap_t]] = 1.
            end
        end

        if haskey(m.settings, :add_urhat)
            if get_setting(m, :add_urhat)
                ZZ_pseudos[reg][pseudo[:urhat], endo[:L_t]] = -0.7366220966444673
                DD_pseudos[reg][pseudo[:urhat]] = -32.257825316004364 + -0.7366220966444673 * m[:Lmean]
            end
        end

        if haskey(m.settings, :add_rw)
            if get_setting(m, :add_rw)
                ZZ_pseudos[reg][pseudo[:rw], endo[:rw_t]] = 1.
                ZZ_pseudos[reg][pseudo[:Rref], endo[:Rref_t]] = 1.
                DD_pseudos[reg][pseudo[:Rref]]                = m[:Rstarn]
            end
        end

        ## Fundameantal inflation related pseudo-obs
        if subspec(m) in ["ss13", "ss14", "ss15", "ss16", "ss17", "ss18", "ss19", "ss20"] #,
            # "ss21", "ss22", "ss23","ss24", "ss25", "ss26"]
            # Compute coefficient on Sinf
            betabar = exp((1-m[:σ_c] ) * m[:z_star]) * m[:β]
            κ = ((1 - m[:ζ_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))*
                 (1 - m[:ζ_p]))/(m[:ζ_p]*((m[:Φ]- 1)*m[:ϵ_p] + 1))/
            (1 + m[:ι_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))

            κcoef = κ * (1 + m[:ι_p] * betabar)

            ZZ_pseudos[reg][pseudo[:Sinf_t], endo_addl[:Sinf_t]] = 1.
            ZZ_pseudos[reg][pseudo[:Sinf_w_coef_t], endo_addl[:Sinf_t]] = κcoef
            DD_pseudos[reg][pseudo[:ι_p]] = m[:ι_p]
            ZZ_pseudos[reg][pseudo[:πtil_t], endo_addl[:πtil_t]] = 1.
            DD_pseudos[reg][pseudo[:πtil_t]] = 100 * (m[:π_star] - 1)
            ZZ_pseudos[reg][pseudo[:e_tfp_t], endo_addl[:e_tfp_t]] = 1.
            if subspec(m) in ["ss14", "ss15", "ss16", "ss18", "ss19"]
                ZZ_pseudos[reg][pseudo[:e_tfp_t1], endo_addl[:e_tfp_t1]] = 1.
            end
        end

        ## Exogenous processes
        if subspec(m) == "ss12"
            to_add = [:g_t, :b_t, :μ_t, :z_t, :λ_f_t, :λ_w_t, :rm_t, :σ_ω_t, :μ_e_t,
                      :γ_t, :π_star_t]
            to_add_addl = [:e_lr_t, :e_tfp_t, :e_gdpdef_t, :e_corepce_t, :e_gdp_t, :e_gdi_t]
            for i in to_add
                ZZ_pseudos[reg][pseudo[i], endo[i]] = 1.
            end
            for i in to_add_addl
                ZZ_pseudos[reg][pseudo[i], endo_addl[i]] = 1.
            end
        end
    end

    for para in m.parameters
        if !isempty(para.regimes)
            ModelConstructors.toggle_regime!(para, 1)
        end
    end

    return [PseudoMeasurement(ZZ_pseudos[i], DD_pseudos[i]) for i in 1:n_reg]
end


text:
using Celeste.DeterministicVIImagePSF

using Celeste.DeterministicVI.CanonicalParams
using Celeste.SensitiveFloats.zero_sensitive_float_array
using Celeste.SensitiveFloats.SensitiveFloat
using Celeste.SensitiveFloats.clear!

using ForwardDiff
using Base.Test

# Overload some base functions for testing with ForwardDiff
import Base.floor
function floor{N}(x::ForwardDiff.Dual{N,Float64})
    return floor(x.value)
end

function floor{N}(x::ForwardDiff.Dual{N,ForwardDiff.Dual{N,Float64}})
    return floor(x.value)
end

import Base.sinc
function sinc{N}(x::ForwardDiff.Dual{N,Float64})
    # Note that this won't work with x = 0, but it's ok for testing.
    return sin(x * pi) / (x * pi)
end

function sinc{N}(x::ForwardDiff.Dual{N,ForwardDiff.Dual{N,Float64}})
    # Note that this won't work with x = 0.
    return sin(x * pi) / (x * pi)
end

function test_convolve_sensitive_float_matrix()
    # Use the FSMSensitiveFloatMatrices because it initializes all the
    # sizes for us automatically.
    fsms = DeterministicVIImagePSF.FSMSensitiveFloatMatrices();
    psf_image = zeros(3, 3);
    psf_image[2, 2] = 0.5;
    psf_image[2, 1] = psf_image[1, 2] = 0.25;
    DeterministicVIImagePSF.initialize_fsm_sf_matrices_band!(
        fsms, 1, 1, 1,
        1, 1, 3, 3, psf_image)

    sf = SensitiveFloat{Float64}(length(GalaxyPosParams), 1, true, true)
    sf.v[] = 3;
    sf.d[:, 1] = rand(size(sf.d, 1))
    h = rand(size(sf.h))
    sf.h[:] = h * h';
    fsms.fs1m_image_padded[2, 2] = deepcopy(sf);
    DeterministicVIImagePSF.convolve_sensitive_float_matrix!(
        fsms.fs1m_image_padded, fsms.psf_fft, fsms.fs1m_conv_padded);
    h_indices = (1:3) + fsms.pad_pix_h
    w_indices = (1:3) + fsms.pad_pix_w
    conv_image =
        Float64[ sf.v[] for sf in fsms.fs1m_conv_padded ][h_indices, w_indices];
    @test_approx_eq(sf.v[] * psf_image, conv_image)

    for ind in 1:size(sf.d, 1)
        conv_image =
            Float64[ sf.d[ind] for sf in fsms.fs1m_conv_padded ][h_indices, w_indices];
        @test_approx_eq(sf.d[ind] * psf_image, conv_image)
    end

    for ind1 in 1:size(sf.h, 1), ind2 in 1:size(sf.h, 2)
        conv_image =
            Float64[ sf.h[ind1, ind2] for sf in fsms.fs1m_conv_padded ][h_indices, w_indices];
        @test_approx_eq(sf.h[ind1, ind2] * psf_image, conv_image)
    end

end

function test_sinc()
    function sinc_with_derivatives_fd{T <: Number}(x::Vector{T})
        v, d, h = DeterministicVIImagePSF.sinc_with_derivatives(x[1])
        return v
    end

    x = 0.7
    fd_v = sinc_with_derivatives_fd(Float64[ x ])
    fd_d = ForwardDiff.gradient(sinc_with_derivatives_fd, Float64[ x ])[1]
    fd_h = ForwardDiff.hessian(sinc_with_derivatives_fd, Float64[ x ])[1, 1]

    v, d, h = DeterministicVIImagePSF.sinc_with_derivatives(x)

    @test_approx_eq sinc(x) v
    @test_approx_eq fd_v v
    @test_approx_eq fd_d d
    @test_approx_eq fd_h h
end

function test_lanczos_kernel()
    kernel_width = 2.0
    function lanczos_kernel_fd{NumType <: Number}(x_vec::Vector{NumType})
        v, d, h = DeterministicVIImagePSF.lanczos_kernel_with_derivatives_nocheck(
            x_vec[1], kernel_width)
        return v
    end

    for x in Float64[-2.2, -1.2, -0.2, 0, 0.2, 1.2, 2.2]
        fd_v = lanczos_kernel_fd([ x ])
        fd_d = ForwardDiff.gradient(lanczos_kernel_fd, Float64[ x ])[1]
        fd_h = ForwardDiff.hessian(lanczos_kernel_fd, Float64[ x ])[1, 1]

        v, d, h = DeterministicVIImagePSF.lanczos_kernel_with_derivatives_nocheck(x, kernel_width)

        @test_approx_eq fd_v v
        @test_approx_eq fd_d d
        @test_approx_eq fd_h h
    end
end

function test_bspline_kernel()
    function bspline_kernel_fd{NumType <: Number}(x_vec::Vector{NumType})
        v, d, h = DeterministicVIImagePSF.bspline_kernel_with_derivatives(x_vec[1])
        return v
    end

    for x in Float64[-2.2, -1.2, -0.2, 0, 0.2, 1.2, 2.2]
        fd_v = bspline_kernel_fd([ x ])
        fd_d = ForwardDiff.gradient(bspline_kernel_fd, Float64[ x ])[1]
        fd_h = ForwardDiff.hessian(bspline_kernel_fd, Float64[ x ])[1, 1]

        v, d, h = DeterministicVIImagePSF.bspline_kernel_with_derivatives(x)

        @test_approx_eq fd_v v
        @test_approx_eq fd_d d
        @test_approx_eq fd_h h
    end
end

function test_cubic_kernel()
    kernel_param = -0.75

    function cubic_kernel_fd{NumType <: Number}(x_vec::Vector{NumType})
        v, d, h = DeterministicVIImagePSF.cubic_kernel_with_derivatives(
            x_vec[1], kernel_param)
        return v
    end

    for x in Float64[-2.2, -1.2, -0.2, 0, 0.2, 1.2, 2.2]
        fd_v = cubic_kernel_fd([ x ])
        fd_d = ForwardDiff.gradient(cubic_kernel_fd, Float64[ x ])[1]
        fd_h = ForwardDiff.hessian(cubic_kernel_fd, Float64[ x ])[1, 1]

        v, d, h = DeterministicVIImagePSF.cubic_kernel_with_derivatives(x, kernel_param)

        @test_approx_eq fd_v v
        @test_approx_eq fd_d d
        @test_approx_eq fd_h h
    end
end

function test_interpolate()
    psf_image = zeros(Float64, 5, 5);
    psf_image[3, 3] = 0.5
    psf_image[2, 3] = psf_image[3, 2] = psf_image[4, 3] = psf_image[3, 4] = 0.125

    local wcs_jacobian = Float64[0.9 0.2; 0.1 0.8]
    local world_loc = Float64[5.1, 5.2]
    local kernel_width = 2

    image_size = (11, 11)

    function test_kernel(kernel_fun)
        function interpolate_loc{T <: Number}(
            world_loc::Vector{T}, calculate_gradient::Bool)
            local image = zero_sensitive_float_array(T, length(StarPosParams), 1,
                                                                image_size...)
            local pixel_loc = Celeste.Model.linear_world_to_pix(
                wcs_jacobian, Float64[0., 0.], Float64[1.0, 0.5], world_loc)
            DeterministicVIImagePSF.interpolate!(
                kernel_fun, kernel_width, image, psf_image, pixel_loc, wcs_jacobian,
                calculate_gradient, calculate_gradient)
            return image
        end

        image = interpolate_loc(world_loc, true)
        for test_pix in prod(image_size)
            function interpolate_loc_fd{T <: Number}(world_loc::Vector{T})
                local image = interpolate_loc(world_loc, false)
                return image[test_pix].v[]
            end

            fd_v = interpolate_loc_fd(world_loc)
            fd_d = ForwardDiff.gradient(interpolate_loc_fd, world_loc)
            fd_h = ForwardDiff.hessian(interpolate_loc_fd, world_loc)

            @test_approx_eq image[test_pix].v[] fd_v
            @test_approx_eq image[test_pix].d fd_d
            @test_approx_eq image[test_pix].h fd_h
        end
    end

    println("Testing cubic kernel")
    test_kernel(x ->
        DeterministicVIImagePSF.cubic_kernel_with_derivatives(x, -0.75))

    println("Testing bspline kernel")
    test_kernel(DeterministicVIImagePSF.bspline_kernel_with_derivatives)

    println("Testing lanczos kernel")
    test_kernel(x ->
        DeterministicVIImagePSF.lanczos_kernel_with_derivatives(
            x, Float64(kernel_width)))
end

test_sinc()
test_lanczos_kernel()
test_bspline_kernel()
test_cubic_kernel()
test_convolve_sensitive_float_matrix()
test_interpolate()


text:
# # Ordinary least squares
#
# ## Packages

using CairoMakie
using CalibrationErrors
using CalibrationErrorsDistributions
using CalibrationTests
using Distributions
using StatsBase

using Random

using CairoMakie.AbstractPlotting.ColorSchemes: Dark2_8

## set random seed
Random.seed!(1234)

## create path before saving
function wsavefig(file, fig=current_figure())
    mkpath(dirname(file))
    return save(file, fig)
end

# ## Regression problem
#
# We consider a regression problem with scalar feature $X$ and scalar
# target $Y$ with input-dependent Gaussian noise that is inspired by a
# problem by [Gustafsson, Danelljan, and Schön](https://arxiv.org/abs/1906.01620).
# Feature $X$ is distributed uniformly at random in $[-1, 1]$,
# and target $Y$ is distributed according to
# ```math
# Y \sim \sin(\pi X) + | 1 + X | \epsilon,
# ```
# where $\epsilon \sim \mathcal{N}(0, 0.15^2)$.

# We start by generating a data set consisting of 100 i.i.d. pairs of
# feature $X$ and target $Y$:

xs = rand(Uniform(-1, 1), 100)
ys = rand.(Normal.(sinpi.(xs), 0.15 .* abs.(1 .+ xs)))

# ## Ordinary least squares regression
#
# We perform ordinary least squares regression for this nonlinear heteroscedastic
# regression problem, and train a model $P$ with homoscedastic variance. The fitted
# parameters of the model are

bs = hcat(ones(length(xs)), xs) \ ys

# and the standard deviation of the model is given by

stddev = std(bs[1] .+ bs[2] .* xs .- ys)

# The following plot visualizes the training data set and model $P$, together
# with the function $f(x) = \mathbb{E}[Y | X = x] = \sin(\pi x)$.

fig = Figure(; resolution=(960, 450))

## plot the data generating distribution
ax1 = Axis(fig[1, 1]; title="ℙ(Y|X)", xlabel="X", ylabel="Y")
heatmap!(
    -1:0.01:1,
    -2:0.01:2,
    (x, y) -> pdf(Normal(sinpi(x), 0.15 * abs(1 + x)), y);
    colorrange=(0, 1),
)
scatter!(xs, ys; color=Dark2_8[2])
tightlimits!(ax1)

## plot the predictions of the model
ax2 = Axis(fig[1, 2]; title="P(Y|X)", xlabel="X", ylabel="Y")
heatmap!(
    -1:0.01:1,
    -2:0.01:2,
    let offset = bs[1], slope = bs[2], stddev = stddev
        (x, y) -> pdf(Normal(offset + slope * x, stddev), y)
    end;
    colorrange=(0, 1),
)
scatter!(xs, ys; color=Dark2_8[2])
tightlimits!(ax2)

## link axes and hide y labels and ticks of the second plot
linkaxes!(ax1, ax2)
hideydecorations!(ax2; grid=false)

## add a colorbar
Colorbar(fig[1, 3]; label="density", width=30)

## adjust space
colgap!(fig.layout, 50)

wsavefig("figures/ols/heatmap.pdf") #jl
#!jl wsavefig("figures/ols/heatmap.svg");
#!jl # ![](figures/ols/heatmap.svg)

# ## Validation
#
# We evaluate calibration of the model with a validation data set of $n = 50$ i.i.d. pairs
# of samples $(X_1, Y_1), \ldots, (X_n, Y_n)$ of $(X, Y)$.

valxs = rand(Uniform(-1, 1), 50)
valys = rand.(Normal.(sinpi.(valxs), 0.15 .* abs.(1 .+ valxs)))

# For these validation data points we compute the predicted distributions $P(Y | X = X_i)$.

valps = Normal.(bs[1] .+ bs[2] .* valxs, stddev)

# ## Quantile calibration
#
# We evaluate the predicted cumulative probability $\tau_i = P(Y \leq Y_i | X = X_i)$ for
# each validation data point.

τs = cdf.(valps, valys)

# The following plot visualizes the empirical cumulative distribution function of the
# predicted quantiles.

fig = Figure(; resolution=(600, 450))

ax = Axis(
    fig[1, 1];
    xlabel="quantile level",
    ylabel="cumulative probability",
    xticks=0:0.25:1,
    yticks=0:0.25:1,
    autolimitaspect=1,
    rightspinevisible=false,
    topspinevisible=false,
    xgridvisible=false,
    ygridvisible=false,
)

## plot the ideal
lines!([0, 1], [0, 1]; label="ideal", linewidth=2, color=Dark2_8[1])

## plot the empirical cdf
sort!(τs)
ecdf_xs = vcat(0, repeat(τs; inner=2), 1)
ecdf_ys = repeat(range(0, 1; length=length(τs) + 1); inner=2)
lines!(ecdf_xs, ecdf_ys; label="data", linewidth=2, color=Dark2_8[2])

## add legend
Legend(fig[1, 2], ax; valign=:top, framevisible=false)

## set limits and aspect ratio
colsize!(fig.layout, 1, Aspect(1, 1))
tightlimits!(ax)

wsavefig("figures/ols/quantiles.pdf") #jl
#!jl wsavefig("figures/ols/quantiles.svg");
#!jl # ![](figures/ols/quantiles.svg)

# ## Calibration test
#
# We compute a $p$-value estimate of the null hypothesis that model $P$ is calibrated using
# an estimation of the quantile of the asymptotic distribution of
# $n \widehat{\mathrm{SKCE}}_{k,n}$ with 100000 bootstrap samples on the validation data
# set. Kernel $k$ is chosen as the tensor product kernel
# ```math
# \begin{aligned}
# k\big((p, y), (p', y')\big) &= \exp{\big(- W_2(p, p')\big)} \exp{\big(-(y - y')^2/2\big)} \\
# &= \exp{\big(-\sqrt{(m_p - m_{p'})^2 + (\sigma_p - \sigma_{p'})^2}\big)} \exp{\big( - (y - y')^2/2\big)},
# \end{aligned}
# ```
# where $W_2$ is the 2-Wasserstein distance and $m_p, m_{p'}$ and $\sigma_p, \sigma_{p'}$
# denote the mean and the standard deviation of the normal distributions $p$ and $p'$.

## define kernel
kernel = WassersteinExponentialKernel() ⊗ SqExponentialKernel()

## compute p-value estimate using bootstrapping
pvalue(AsymptoticSKCETest(kernel, valps, valys); bootstrap_iters=100_000)

# We obtain $p < 0.05$ in our experiment, and hence the calibration test rejects $H_0$ at
# the significance level $\alpha = 0.05$.