text:
module TestCancellableContainers

using ArgCheck: @check
using ProgressLogging: Progress
using Reagents.Internal: Waiting, @trace, offerid
using UUIDs: uuid4
using ..Utils: random_sleep, run_concurrently

include("../../../examples/cancellablecontainers.jl")

blocking_treiberstack(T) = Blocking(TreiberStack{T}())
blocking_msqueue(T) = Blocking(MSQueue{T}())

function test_cancellation_token_internal()
    items = blocking_msqueue(Int)
    token = CancellationToken()

    t = @task cancellable_take!(items, token)
    yield(t)

    # At this point, the `offer` of `cancellable_take!` is registered into the
    # channels for the item taker *and* the cancellation listener:
    @test items.send.dual.next.value.offer === token.receive.msgs.next.value.offer
    # And it is waiting:
    @test items.send.dual.next.value.offer.state[] isa Waiting
end

function test_repeat_cancellation(; parentid = uuid4(), kwargs...)
    id = uuid4()
    i = 0
    try
        @testset for spawn in [false, true]
            @debug Progress(
                id,
                i / 2;
                name = "`test_repeat_cancellation`: spawn=$spawn",
                parentid,
            )
            i += 1

            test_repeat_cancellation(spawn; parentid = id, kwargs...)
        end
    finally
        @debug Progress(id; done = true)
    end
end

function test_repeat_cancellation(spawn::Bool; parentid = uuid4(), kwargs...)
    id = uuid4()
    i = 0
    try
        @testset for f in [blocking_msqueue, blocking_treiberstack]
            @debug Progress(id, i / 2; name = "`test_repeat_cancellation`: $f", parentid)
            i += 1

            test_repeat_cancellation(f, spawn; parentid = id, kwargs...)
        end
    finally
        @debug Progress(id; done = true)
    end
end

function test_repeat_cancellation(
    constructor,
    spawn::Bool;
    parentid = uuid4(),
    nrepeat = 10,
    kwargs...,
)
    id = uuid4()
    try
        for i in 1:nrepeat
            @debug Progress(
                id,
                (i - 1) / nrepeat;
                name = "`test_repeat_cancellation`: trial $i",
                parentid,
            )
            check_repeat_cancellation(constructor, spawn; kwargs...)
        end
    finally
        @debug Progress(id; done = true)
    end
    @test true
end

function check_repeat_cancellation(constructor, spawn::Bool; nitems = 100, randomize = true)
    items = constructor(Int)
    token = CancellationToken()
    global ITEMS = items
    global TOKEN = token

    received = [Int[] for _ in 1:5]
    receivers = map(eachindex(received), received) do i, dest
        function receiver_task()
            @trace(label = :start_receiver, i, taskid = objectid(current_task()))
            while true
                @trace(label = :taking, i, taskid = objectid(current_task()))
                y = cancellable_take!(items, token)
                if y isa Cancelled
                    @trace(label = :received_cancelled, i, taskid = objectid(current_task()))
                    break
                end
                @trace(label = :took, i, value = y, taskid = objectid(current_task()))
                push!(dest, y)
                randomize && random_sleep()
            end
            @debug "`check_repeat_cancellation`: Reciever $i done"
        end
    end

    nsenders = 5
    nfinished = Threads.Atomic{Int}(0)
    senders = map(1:nsenders) do i
        function sender_task()
            @trace(label = :start_sender, i, taskid = objectid(current_task()))
            for j in 1:nitems
                put!(items, nsenders * (j - 1) + i)
                randomize && random_sleep()
            end
            if Threads.atomic_add!(nfinished, 1) + 1 == nsenders
                # Note: using `nfinished` instead of using a task that waits for
                # all senders, so that we can use a single `run_concurrently`
                # call (to make debugging easier).

                sleep(0.2)  # TODO: don't do this; support `close`
                # `sleep` here is just a "hint" but it's still bad...

                @debug "`check_repeat_cancellation`: All sent"
                for msg in token.receive.msgs
                    @trace(
                        label = :check_token,
                        i,
                        taskid = objectid(current_task()),
                        offerid = offerid(msg.offer)
                    )
                end
                @trace(label = :sending_cancelled, i, taskid = objectid(current_task()))
                cancel!(token)
                @trace(label = :sent_cancelled, i, taskid = objectid(current_task()))
            end
            @debug "`check_repeat_cancellation`: Sender $i done"
        end
    end
    allsent = 1:nsenders*nitems

    run_concurrently([senders; receivers]; spawn)
    @debug "`check_repeat_cancellation`: All sendres and receivers done"
    @check token.iscancelled[]

    allreceived = reduce(vcat, received)
    nleft = 0
    ncleaned = 0
    while true
        ref = Reagents.try(taking(items))
        ref === nothing && break  # empty
        ncleaned += 1
        y = something(ref)[]
        y === nothing && continue  # taken
        push!(allreceived, something(y))
        nleft += 1
    end
    @debug "`check_repeat_cancellation`: $nleft items not recieved; $ncleaned refs cleaned"

    # Not using `@test` to avoid overhead in while "fuzzing"
    @check length(allreceived) == length(allsent)
    @check sort!(allreceived) == allsent
end

end  # module


text:
using DReal

import DReal:  opensmt_init, opensmt_mk_context,
               opensmt_mk_real_var, opensmt_mk_eq, opensmt_assert, opensmt_check,
               opensmt_mk_sin, opensmt_mk_cos, opensmt_get_ub, opensmt_get_lb,
               opensmt_set_precision, opensmt_mk_ite, opensmt_mk_num, opensmt_mk_gt

using Compat
# A = Var(Float64, "A", 0.0, 1.0)
# B = Var(Float64, "B", 0.0, 1.0)
# # C = Var(Float64, "C", 0.0, 1.0)
# add!(B > ifelse(A>0.4, 0.3,0.4))
# add!(A>0.4)
# model(A)
# model(B)

#
ctx = opensmt_mk_context(@compat UInt32(1))
opensmt_set_precision(ctx, 0.00001)
A = opensmt_mk_real_var( ctx, "A" , 0.0, 1.0)
B = opensmt_mk_real_var( ctx, "B" , 0.0, 1.0)

ite = opensmt_mk_ite(ctx, opensmt_mk_gt(ctx, A, opensmt_mk_num(ctx, 0.4)),
                          opensmt_mk_num(ctx, 0.3),
                          opensmt_mk_num(ctx, 0.4))
eq = opensmt_mk_gt(ctx, B, ite)
opensmt_assert( ctx, eq)
res = opensmt_check( ctx );
println("res is ", res == -1 ? "unsat" : "sat")

ext: jl


text:
function linshrink(F::Union{UniformScaling, AbstractMatrix},
                   S::AbstractMatrix, λ::Real)
    return (1.0 - λ)*S + λ*F
end

function linshrink!(F::AbstractMatrix, S::AbstractMatrix, λ::Real)
    F .*= λ
    F .+= (1.0 - λ) * S
    return F
end


text:
#-------------------------------------------------------------------------------
# File: output.jl
# Description: This files contains all functions and strucutes that are used for
# outputs.
# Date: November 17, 2019
# Authors: Killian Fretaud, Rémi Garcia,
#-------------------------------------------------------------------------------

mutable struct Machine
    starting_times::Vector{Int}
    processing_times::Vector{Int}
end

mutable struct Robot
    moves::Vector{Int}
    set_up_times::Vector{Int}
    transport_times::Vector{Int}
end

mutable struct Solution
    C_max::Int
    order_of_jobs::Vector{Int}
    machines::Vector{Machine}
    robot::Robot
end

function solve_to_represent(instance::Instance, z, sol)
    m = instance.nb_machines
    n = instance.nb_jobs

    jobs_ordered = Vector{Int}(undef, n)
    machines = Vector{Machine}(undef, m)
    for j in 1:m
        times = zeros(Int, n)
        for i in 1:n
            times[i] = round(Int, value.(sol.sm[i, j]))
        end
        machines[j] = Machine(times, instance.processing_times[:, j])
        #machines[j] = Machine(value.(sol.sm[:, j]), instance.processing_times[:, j])
    end

    precedence = sum(value.(sol.y), dims=2)
    index = 1 ; next_is = n-1
    while index <= n
        job = 1
        while precedence[job] != next_is && job <= n
            job += 1
        end
        jobs_ordered[index] = job
        index += 1 ; next_is -= 1
    end

    moves = zeros(Int, instance.nb_robot)
    occurence = ones(Int, n) # sur quelle prochaine machine doit aller la tâche
    st = zeros(Int, instance.nb_robot)

    # For all move
    for pos in 1:instance.nb_robot
        st[pos] = round(Int, value.(sol.st[pos]))
        for job in 1:n
            # This is the job taken
            if occurence[job] <= (m + 1)
                if Bool(round(Int, value.(sol.x[job, occurence[job], pos])))
                    moves[pos] = job
                    occurence[job] += 1
                end
            end
        end
    end

    return Solution(round(Int, z), jobs_ordered, machines, Robot(moves, st, instance.travel_times))
end

function print_solution(solution::Solution)
    println("\t= Solution is =============================================")
    println("\tCmax = ", solution.C_max)
    println("\tMachines: processed job (starting time, ending time)")
    for j in 1:length(solution.machines)
        print("\t", j, ":")
        for i in solution.order_of_jobs
            start = solution.machines[j].starting_times[i]
            print("\tjob ", i, "(", start, ", ", start + solution.machines[j].processing_times[i], ")")
        end
        println()
    end
    occurences = ones(Int, length(solution.order_of_jobs))
    println("\tRobot: job moved (set up time + transport time)")
    print("\t\t")
    for k in 1:length(solution.robot.moves)
        job = solution.robot.moves[k]
        print("J", job, "(", solution.robot.set_up_times[k], "+", solution.robot.transport_times[occurences[job]], ") ")
        occurences[job] += 1
    end
    println()
    println("\t===========================================================")
end

max_stars_repo_name: k-fretaud/Flowshop-robotic-cell

max_issues_repo_name: k-fretaud/Flowshop-robotic-cell

max_forks_repo_name: k-fretaud/Flowshop-robotic-cell


text:
using Plots
using Random
Random.seed!(1)

# Differential Dynamic Programming
include_ddp()
include(joinpath(pwd(), "src/differential_dynamic_programming/multiple_models.jl"))

# Model
include_model("double_integrator")

function f(model::DoubleIntegratorContinuous, x, u, w)
    [x[2] + w[1];
	 (1.0 + w[3]) * u[1] + w[2]]
end

model = DoubleIntegratorContinuous{Midpoint, FixedTime}(2, 1, 3)
n = model.n
m = model.m

# Models
N = 1000#2 * n + 1

# Time
T = 101
h = 0.01
tf = h * (T - 1)
t = range(0, stop = tf, length = T)

# Reference position trajectory
# z = range(0.0, stop = 3.0 * 2.0 * π, length = T)
# p_ref = 1.0 * cos.(1.0 * z)
# plot(z, p_ref)
p_ref = [0.0 for t = 1:T]

# Initial conditions, controls, disturbances
# x1 = [p_ref[1]; 0.0]
x1 = [1.0; 0.0]
_ū = [0.1 * rand(model.m) for t = 1:T-1]
ū = [_ū for i = 1:N]

W = Distributions.MvNormal(zeros(model.d), Diagonal([0.0, 0.0, 0.1]))
wi = [vec(rand(W, 1)) for i = 1:N]
w = [[zeros(model.d) for t = 1:T-1], [[wi[i] for t = 1:T-1] for i = 1:N-1]...]
# w = [[vec(rand(W, 1)) for t = 1:T-1] for i = 1:N]
# w = [[zeros(model.d) for t = 1:T-1] for i = 1:N]

# Rollout
x̄ = [rollout(model, x1, ū[i], w[1], h, T) for i = 1:N]

# x_ref = [sum([x̄i[t] for x̄i in x̄]) ./ N for t = 1:T]
# u_ref = [sum([ūi[t] for ūi in ū]) ./ N for t = 1:T-1]

p_idx = (1:1)
plt = plot()
for i = 1:N
	plt = plot!(t, hcat(x̄[i]...)[p_idx, :]',
		label = "")
end
display(plt)

# Objective
Q = [t < T ? h * Diagonal([1.0e-1; 1.0e-1]) : Diagonal([1000.0; 1000.0]) for t = 1:T]
R = [h * Diagonal(1.0e-2 * ones(model.m)) for t = 1:T-1]

obj = StageCosts([QuadraticCost(Q[t], nothing,
	t < T ? R[t] : nothing, nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
	    R = obj.cost[t].R
        return (x - [p_ref[t]; 0.0])' * Q * (x - [p_ref[t]; 0.0]) + u' * R * u
    elseif t == T
		Q = obj.cost[T].Q
        return (x - [p_ref[T]; 0.0])' * Q * (x - [p_ref[T]; 0.0])
    else
        return 0.0
    end
end

function ddp_solve!(prob::ProblemData;
    max_iter = 10,
    grad_tol = 1.0e-5,
    verbose = true)

	println()
    verbose && println("Differential Dynamic Programming")

	# data
	p_data = prob.p_data
	m_data = prob.m_data
	s_data = prob.s_data

    # compute objective
    s_data.obj = objective(m_data, mode = :nominal)

	N = length(prob.m_data)

    for i = 1:max_iter
        # derivatives
        derivatives!(m_data)

        # backward pass
        backward_pass!(p_data, m_data)

		# # set nominal trajectories
		# x_ref = [sum([m.x̄[t] for m in prob.m_data]) ./ N for t = 1:T]
		# u_ref = [sum([m.ū[t] for m in prob.m_data]) ./ N for t = 1:T-1]
		#
		# # disturbances
		# w = [[rand(W, 1) for t = 1:T-1] for i = 1:N]
		#
		# for i = 1:N
		# 	# m_data[i].x̄ .= deepcopy(x_ref)
		# 	# m_data[i].ū .= deepcopy(u_ref)
		# 	m_data[i].w .= deepcopy(w[i])
		# end

        # forward pass
        forward_pass!(p_data, m_data, s_data)

		if i > 1
			for i = 2:N
				m_data[i].x̄ .= deepcopy(m_data[1].x̄)
				m_data[i].ū .= deepcopy(m_data[1].ū)
				# m_data[i].w .= deepcopy(w[i])
			end
		end
		# # set nominal trajectories
		# x_ref = [sum([m.x̄[t] for m in prob.m_data]) ./ N for t = 1:T]
		# u_ref = [sum([m.ū[t] for m in prob.m_data]) ./ N for t = 1:T-1]

        # check convergence
        grad_norm = norm(s_data.gradient, Inf)
        verbose && println("     iter: $i
             cost: $(s_data.obj)
			 grad_norm: $(grad_norm)")
		grad_norm < grad_tol && break
        !s_data.status && break
    end
end

models_data = models(model, obj, deepcopy(x̄), deepcopy(ū), w, h, T)

prob = problem_data(models_data)

# Solve
# for j = 1:10
# 	wi = [vec(rand(W, 1)) for i = 1:N]
# 	w = [[zeros(model.d) for t = 1:T-1], [[wi[i] for t = 1:T-1] for i = 1:N-1]...]
# 	# w = [[vec(rand(W, 1)) for t = 1:T-1] for i = 1:N]
# 	# w = [[zeros(model.d) for t = 1:T-1] for i = 1:N]
# 	for i = 1:N
# 		prob.m_data[i].w .= w[i]
# 		prob.m_data[i].x̄ .= rollout(model,
# 			prob.m_data[i].x̄[1], prob.m_data[1].ū,
# 			prob.m_data[1].w, h, T)
# 	end
# 	# Rollout
#
# 	@time ddp_solve!(prob,
# 	    max_iter = 1000, verbose = true,
# 		grad_tol = 1.0e-8)
# end
@time ddp_solve!(prob,
	max_iter = 1000, verbose = true,
	grad_tol = 1.0e-8)

x = [m.x for m in prob.m_data]
u = [m.u for m in prob.m_data]

x̄ = [m.x̄ for m in prob.m_data]
ū = [m.ū for m in prob.m_data]

x_ref = deepcopy(prob.m_data[1].x̄)#[sum([m.x̄[t] for m in prob.m_data]) ./ N for t = 1:T]
u_ref = deepcopy(prob.m_data[1].ū)#[sum([m.ū[t] for m in prob.m_data]) ./ N for t = 1:T-1]

# Visualize
idx = (1:2)
_plt = plot(hcat([[p_ref[t]; 0.0][idx] for t = 1:T]...)',
    width = 2.0, color = :black, label = "")
for i = 1:N
	_plt = plot!(hcat(x[i]...)[idx, :]', color = :magenta, label = "")
end
_plt = plot!(hcat(x_ref...)[idx, :]', color = :orange, width = 5.0, label = "")
display(_plt)
# plot(hcat(u..., u[end])', linetype = :steppost)

# Simulate policy
include_dpo()
include(joinpath(pwd(), "examples/direct_policy_optimization/simulate.jl"))

# Model
model_sim = model
x1_sim = copy(x1)
T_sim = 10 * (T-1) + 1

# Disturbance distributions

# Time
t_sim = range(0, stop = tf, length = T_sim)
dt_sim = tf / (T_sim - 1)

# Policy
policy = linear_feedback(model.n, model.m)
# K, P = tvlqr(model, x̄, ū, h, [Q for t = 1:T], [R for t = 1:T-1])
K = [-K for K in prob.p_data.K]

# Simulate
N_sim = 100
x_sim = []
u_sim = []
J_sim = []
Random.seed!(1)
W_sim = W#Distributions.MvNormal(zeros(model.d), Diagonal([0.0, 0.0, 1.0e-1]))

for i = 1:N_sim
	println("sim: $i")
	wi_sim = vec(rand(W_sim, 1))
	w_sim = hcat([wi_sim for t = 1:T_sim]...)

	x_ddp, u_ddp, J_ddp, Jx_ddp, Ju_ddp = simulate(
		model_sim,
		linear_feedback(model.n, model.m), K,
	    x_ref, u_ref,
		Q, R,
		T_sim, h,
		x1_sim,
		w_sim)

	push!(x_sim, x_ddp)
	push!(u_sim, u_ddp)
	push!(J_sim, J_ddp)
end

# Visualize
idx = (1:2)
_plt = plot(t, hcat(x_ref...)[idx, :]',
    width = 2.0, color = :black, label = "")

for (i, xi) in enumerate(x_sim)
	_plt = plot!(t_sim, hcat(xi...)[idx, :]',
    	width = 1.0, color = :magenta, label = i == 1 ? "sim" : "")
end
_plt = plot!(t, hcat(x_ref...)[idx, :]',
    width = 2.0, color = :black, label = "ref")
_plt = plot!(; xlabel = "time (s)",
	ylabel = "position",
	legend = :topleft,
	title = "N_mc = $N, N_sim = $(N_sim), J_avg = $(round(mean(J_sim), digits = 3))")
display(_plt)

_plt = plot(t, hcat(u_ref..., u_ref[end])[1:1, :]',
    width = 2.0, color = :black, label = "")
for (i, ui) in enumerate(u_sim)
	_plt = plot!(t_sim, hcat(ui..., ui[end])[1:1, :]',
    	width = 1.0, color = :magenta, label = i == 1 ? "sim" : "")
end
_plt = plot!(t, hcat(u_ref..., u_ref[end])[1:1, :]',
    width = 2.0, color = :black, label = "ref", linetype = :steppost)
_plt = plot!(; xlabel = "time (s)",
	ylabel = "control",
	title = "N_mc = $N, N_sim = $(N_sim), J_avg = $(round(mean(J_sim), digits = 3))",
	linetype = :steppost)
display(_plt)
@show mean(J_sim)
round(mean(J_sim), digits = 3)


text:
# ellipse_im.jl
# Copyright 2019-03-05, Jeff Fessler, University of Michigan

#using MIRT: MIRT_image_geom, jim
using Plots

"""
`phantom = ellipse_im(ig, params;
rot=0, oversample=1, hu_scale=1, replace=false, return_params=false)`

Generate ellipse phantom image from parameters:

`[x_center y_center x_radius y_radius angle_degrees amplitude]`

in
* `ig`		from `image_geom()`
* `params`	[ne 6]	ellipse parameters

option
* `rot`			rotate ellipses by this amount [degrees]
* `oversample`	oversampling factor, for grayscale boundaries
* `hu_scale`	use 1000 to scale shepp-logan to HU
* `replace`		replace ellipse values if true, else add
* `how`			:fast is the only option
* `return_params` if true return (phantom, params)

out
* `phantom`		[nx ny]	image
* `params`		[ne 6] parameters (only if return_params=true)

note: `op ellipse` in aspire with `nsub=3` is `oversample=4 = 2^(3-1)` here

"""
function ellipse_im(ig::MIRT_image_geom,
	params::AbstractArray{<:Real,2};
	rot::Real=0,
	oversample::Integer=1,
	hu_scale::Real=1,
	replace::Bool=false,
	how::Symbol=:fast, # todo
	return_params::Bool=false)

	params[:,6] .*= hu_scale

#if how == :fast && oversample == 1
#	warning("ignoring :fast option for oversample=1")
#	how = :slow
#end

	if oversample > 1
		ig = ig.over(oversample)
	end
	args = (ig.nx, ig.ny, params, ig.dx, ig.dy, ig.offset_x, ig.offset_y,
	rot, oversample, replace)

	if how == :fast
		phantom = ellipse_im_fast(args...)
#	elseif how == :slow
#		phantom = ellipse_im_slow(args...)
	else
		error("bad how")
	end

	if oversample > 1
		phantom = downsample2(phantom, oversample)
	end

	if return_params
		return (phantom, params)
	end
	return phantom
end

#%
#% ellipse_im_slow()
#% brute force fine grid - can use lots of memory
#%
#function phantom = ellipse_im_slow(nx, ny, params, dx, dy, ...
	#offset_x, offset_y, rot, over, replace)
#
#if size(params,2) ~= 6
	#fail('bad ellipse parameter vector size')
#end

#% optional rotation of ellipse parameters
#if rot ~= 0
	#th = deg2rad(rot);
	#cx = params[:,1);
	#cy = params[:,2);
	#params[:,1) = cx * cos(th) + cy * sin(th);
	#params[:,2) = -cx * sin(th) + cy * cos(th);
	#params[:,5) = params[:,5) + rot;
	#clear cx cy th
#end
#
#wx = (nx*over-1)/2 + offset_x * over;
#wy = (ny*over-1)/2 + offset_y * over;
#xx = ((0:nx*over-1) - wx) / over * dx;
#yy = ((0:ny*over-1) - wy) / over * dy;
#[xx yy] = ndgrid(xx, yy); % fine grid, equally spaced
#
#phantom = zeros(nx*over, ny*over, 'single'); % fine array
#
#ticker reset
#ne = nrow(params);
#for ie = 1:ne
	#ticker(mfilename, ie, ne)
#
	#ell = params[ie, :);
	#cx = ell(1);	rx = ell(3);
	#cy = ell(2);	ry = ell(4);
	#theta = deg2rad(ell(5));
	#[xr yr] = rot2(xx-cx, yy-cy, theta);
	#tmp = (xr / rx).^2 + (yr / ry).^2 <= 1;
#
	#if replace
		#phantom(tmp > 0) = ell(6);
	#else
		#phantom = phantom + ell(6) * tmp;
	#end
#end
#
#phantom = downsample2(phantom, over);

"""
`phantom = ellipse_im_fast()`
"""
function ellipse_im_fast(nx, ny, params_in, dx, dy,
	offset_x, offset_y, rot, over, replace)

	params = copy(params_in)
	if size(params,2) != 6
		error("bad ellipse parameter vector size")
	end

	# optional rotation
	if rot != 0
		th = rot * 180/pi
		cx = params[:,1]
		cy = params[:,2]
		params[:,1] = cx * cos(th) + cy * sin(th)
		params[:,2] = -cx * sin(th) + cy * cos(th)
		params[:,5] .+= rot
	end

	phantom = zeros(Float32, nx, ny)

	wx = (nx-1)/2 + offset_x
	wy = (ny-1)/2 + offset_y
	x1 = ((0:nx-1) .- wx) * dx
	y1 = ((0:ny-1) .- wy) * dy
	(xx, yy) = ndgrid(x1, y1)

	hx = abs(dx) / 2
	hy = abs(dy) / 2

	for ie = 1:size(params,1)

		ell = params[ie, :]
		cx = ell[1];	rx = ell[3]
		cy = ell[2];	ry = ell[4]
		theta = ell[5] * pi/180

		xs = xx .- cx # shift per ellipse center
		ys = yy .- cy

		# coordinates of "outer" corner of each pixel, relative to ellipse center
		xo = xs + sign.(xs) * hx
		yo = ys + sign.(ys) * hy

		# voxels that are entirely inside the ellipse:
		(xr, yr) = rot2(xo, yo, theta)
		is_inside = (xr / rx).^2 + (yr / ry).^2 .<= 1

		if replace
			phantom(is_inside) .= ell[6]
		else
			phantom += ell[6] * Float32.(is_inside)
		end

	end # ie loop

	return phantom
end # ellipse_im_fast()

"""
`phantom = ellipse_im(nx, dx, params; kwarg...)`

square image of size `nx` by `nx`,
specifying pixel size `dx` and ellipse `params`
"""
function ellipse_im(nx::Integer, dx::Real, params; kwarg...)
	ig = image_geom(nx=nx, dx=1)
	return ellipse_im(ig, params; kwarg...)
end

"""
`phantom = ellipse_im(nx::Integer, params; kwarg...)`

square image of size `nx` by `nx` with
pixel size `dx=1` and ellipse `params`
"""
function ellipse_im(nx::Integer, params; kwarg...)
	return ellipse_im(nx, 1., params; kwarg...)
end

"""
`phantom = ellipse_im(nx::Integer; ny::Integer=nx, dx::Real=1, kwarg...)`

image of size `nx` by `ny` (default `nx`) with specified `dx` (default 1),
defaults to `:shepplogan_emis`
"""
function ellipse_im(nx::Integer; ny::Integer=nx, dx::Real=1, kwarg...)
	ig = image_geom(nx=nx, ny=ny, dx=dx)
	return ellipse_im(ig, :shepplogan_emis; kwarg...)
end

"""
`phantom = ellipse_im(nx::Integer, ny::Integer; kwarg...)`

`:shepplogan_emis` of size `nx` by `ny`
"""
function ellipse_im(nx::Integer, ny::Integer; kwarg...)
	return ellipse_im(nx, ny=ny, dx=1.; kwarg...)
end

"""
`phantom = ellipse_im(ig, code, kwarg...)`

`code = :shepplogan | :shepplogan_emis | :shepplogan_brainweb`
"""
function ellipse_im(ig::MIRT_image_geom, params::Symbol; kwarg...)
	fov = ig.fov
	if params == :shepplogan
		params = shepp_logan_parameters(fov, fov, case=:kak)
	elseif params == :shepplogan_emis
		params = shepp_logan_parameters(fov, fov, case=:emis)
	elseif params == :shepplogan_brainweb
		params = shepp_logan_parameters(fov, fov, case=:brainweb)
	else
		error("bad phantom symbol")
	end
	return ellipse_im(ig, params; kwarg...)
end

"""
`phantom = ellipse_im(ig; kwarg...)`

`:shepplogan` (default) for given image geometry `ig`
"""
function ellipse_im(ig::MIRT_image_geom; kwarg...)
	return ellipse_im(ig, :shepplogan; kwarg...)
end

"""
`(xx,yy) = ndgrid(x,y)`
"""
function ndgrid(x::AbstractVector{<:Number},
				y::AbstractVector{<:Number})
	return (repeat(x, 1, length(y)), repeat(y', length(x), 1))
end

"""
`(xr,yr) = rot2(x, y, theta)`
2D rotation
"""
function rot2(x, y, theta)
	xr = cos(theta) * x + sin(theta) * y
	yr = -sin(theta) * x + cos(theta) * y
	return (xr, yr)
end

"""
`params = shepp_logan_parameters(xfov, yfov)`

parameters from Kak and Slaney text, p. 255

the first four columns are unitless "fractions of field of view"
"""
function shepp_logan_parameters(xfov::Real, yfov::Real; case::Symbol=:kak)
	params = [
	0		0		0.92	0.69	90	2;
	0		-0.0184	0.874	0.6624	90	-0.98;
	0.22	0		0.31	0.11	72	-0.02;
	-0.22	0		0.41	0.16	108	-0.02;
	0		0.35	0.25	0.21	90	0.01;
	0		0.1		0.046	0.046	0	0.01;
	0		-0.1	0.046	0.046	0	0.01;
	-0.08	-0.605	0.046	0.023	0	0.01;
	0		-0.605	0.023	0.023	0	0.01;
	0.06	-0.605	0.046	0.023	90	0.01];

	params[:,[1,3]] .*= xfov/2
	params[:,[2,4]] .*= yfov/2

	if case == :emis
		params[:,6] = [1, 1, -2, 2, 3, 4, 5, 6, 1, 1]
	elseif case == :brainweb
		params[:,6] = [1, 0, 2, 3, 4, 5, 6, 7, 8, 9] # brainweb uses index 1-10
	elseif case != :kak
		error("bad phantom case")
	end

	return params
end

#%
#% ellipse_im_profile()
#%
#function ellipse_im_profile
#ig = image_geom('nx', 2^9, 'ny', 2^9+2', 'fov', 250);
#profile on
#x0 = ellipse_im(ig, [], 'oversample', 3, 'type', 'fast');
#profile off
#profile report

"""
`ellipse_im()`

show docstring(s)
"""
function ellipse_im()
	@doc ellipse_im
end

"""
`ellipse_im_show()`
"""
function ellipse_im_show()
	ig = image_geom(nx=2^8, ny=2^8+2, fov=250)
#	ig.offset_y = 75.6 / ig.dy

	over = 2^2
	x0 = ellipse_im(ig, oversample=over)
	p1 = jim(ig.x(), ig.y(), x0, title="Shepp Logan", clim=[0.9,1.1])

	x1 = ellipse_im(ig, :shepplogan_emis, oversample=over)
	p4 = jim(x1, title="Shepp Logan Emission")

	x3 = ellipse_im(ig, :shepplogan_brainweb)
	p5 = jim(x3, title="Shepp Logan Brainweb")

	plot(p1,p4,p5)
end

# compare to aspire
function ellipse_im_aspire()
	nx = 2^6;
	ig = image_geom(nx=nx, ny=nx+2, fov=2^7)
	ell = [10, 20, 30, 40, 50, 1]

	ir_test_dir = "/tmp/" # todo
	file = ir_test_dir * "t.fld"
	over = 2^2
	com = sprintf("echo y | op -chat 99 ellipse %s %d %d  %g %g %g %g %g %g %d",
		file, ig.nx, ig.ny, ell ./ [ig.dx, ig.dx, ig.dx, ig.dx, 1, 1],
		log2(over)+1)
	os_run(com) # todo: not done
	asp = fld_read(file)
	jim(asp, title="aspire")

	area_asp = sum(asp[:]) * abs(ig.dx * ig.dy)

	jul = ellipse_im(ig, ell, oversample=over) # 'type', types{ii})
	area_jul = sum(jul[:]) * abs(ig.dx * ig.dy)

	jim(jul, title="jul")
	jim((jul-asp)*over^2, title="difference: (jul-asp)*over^2")

	@show maximum(abs.(jul - mat)) / ell[6] * over^2
	@assert isapprox(jul, mat)
#	max_percent_diff(mat, asp)

	area_real = pi * ell[3] * ell[4] * ell[6]
	@show area_real, area_asp, area_jul

	true
end

# ellipse_im_test()
function ellipse_im_test()
	fov = 100
	shepp_logan_parameters(fov, fov)
	shepp_logan_parameters(fov, fov, case=:kak)
	shepp_logan_parameters(fov, fov, case=:emis)
	shepp_logan_parameters(fov, fov, case=:brainweb)

	# test various ways of calling
	ellipse_im(20)
	ellipse_im(30, :shepplogan_emis, oversample=2)

	ig = image_geom(nx=80, dx=1)
	fov = ig.fov
	params = shepp_logan_parameters(fov, fov, case=:brainweb)
	ellipse_im(ig, params)
	ellipse_im(ig, params, oversample=2)

	if false
		x1 = ellipse_im(100, :shepplogan_emis)
		x2 = ellipse_im(100, :shepplogan_emis, oversample=2)
		plot(jim(x1), jim(x2), jim(x2-x1))
	end

	ellipse_im_show()

#	if has_aspire
#		ellipse_im_aspire() # todo
#	end
	true
end

"""
`ellipse_im(:test)`

`ellipse_im(:show)`

run tests
"""
function ellipse_im(test::Symbol)
	if test == :show
		return ellipse_im_show()
	end
	@assert test == :test
	ellipse_im()
	ellipse_im(:show)
	ellipse_im_test()
end


text:
function test_styles(io::IO = stdout)
    for style in (:bold,
                  :faint,
                  :italics,
                  :underline,
                  :blink,
                  :negative,
                  :conceal,
                  :strikethrough)
        print(io, Crayon(;style => true), "Printed with $style = true", Crayon(reset = true))
        style == :conceal && print(io, "  <- This is concealed = true")
        println(io)
    end
end

function test_system_colors(io::IO = stdout)
    for col in map(first,sort(collect(COLORS), by=last))
        print(io, Crayon(foreground = col), lpad("$col", 15, ' '), " ", Crayon(reset = true))
        print(io, Crayon(background = col), col, Crayon(reset = true), "\n")
    end
end

test_256_colors(codes::Bool = true) = test_256_colors(stdout, codes)
test_256_colors(io::IO) = test_256_colors(io, true)
function test_256_colors(io::IO, codes::Bool)
    println(io, "System colors (0..15):")

    for c in 0:15
        str = codes ? string(lpad(c, 3, '0'), " ") : "██"
        print(io, Crayon(foreground = c), str, Crayon(reset = true))
        (c + 1) % 8 == 0 && println(io)
    end
    print(io, "\n\n")

    println(io, "Color cube, 6×6×6 (16..231):")
    for c in 16:231
        str = codes ? string(lpad(c, 3, '0'), " ") : "██"
        print(io, Crayon(foreground = c), str, Crayon(reset = true))
        (c - 16) %  6 ==  5 && println(io)
        (c - 16) % 36 == 35 && println(io)
    end

    println(io, "Grayscale ramp (232..255):")
    for c in 232:255
        str = codes ? string(lpad(c, 3, '0'), " ") : "██"
        print(io, Crayon(foreground = c), str, Crayon(reset = true))
        (c - 232) %  6 == 5 && println(io)
    end
end

test_24bit_colors(codes::Bool = true) = test_24bit_colors(stdout, codes)
test_24bit_colors(io::IO) = test_24bit_colors(io, true)
function test_24bit_colors(io::IO, codes::Bool)
    steps = 0:30:255
    for r in steps
        for g in steps
            for b in steps
                str = codes ? string(lpad(r, 3, '0'), "|", lpad(g, 3, '0'), "|", lpad(b, 3, '0'), " ") : "██"
                print(io, Crayon(; foreground = (r, g, b)), str, Crayon(reset = true))
            end
            println(io)
        end
        println(io)
    end
end


text:
using CPIDataBase
using CPIDataBase.TestHelpers
using InflationFunctions
using Test

@testset "InflationFunctions.jl" begin
    # Write your tests here.
end

# Pruebas sobre medidas de inflación. Se debe probar instanciar los tipos y que
# definan sus métodos sobre los objetos de CPIDataBase
@testset "InflationSimpleMean" begin
    # Instanciar un tipo
    simplefn = InflationSimpleMean()
    @test simplefn isa InflationSimpleMean

    # Probar que esté definido el método para obtener su nombre
    @test measure_name(simplefn) isa String
    @test measure_tag(simplefn) isa String

    # Probar con bases del IPC con variaciones intermensuales iguales a cero.
    # Estas pruebas ayudan a verificar que la función de inflación se pueda
    # llamar sobre los tipos correctos

    zero_base = getzerobase()
    m_traj_infl = simplefn(zero_base)
    # Probamos que el resumen intermensual sea igual a cero
    @test all(m_traj_infl .≈ 0)

    # Obtenemos un UniformCountryStructure con dos bases y todas las variaciones
    # intermensuales iguales a cero
    zero_cst = getzerocountryst()
    traj_infl = simplefn(zero_cst)

    # Probamos que la trayectoria de inflación sea más larga que el resumen
    # intermensual de una sola base
    @test length(traj_infl) > length(m_traj_infl)

    # Probamos que la trayectoria de inflación sea igual a cero
    @test all(traj_infl .≈ 0)
end

## Pruebas para Inflación de Exclusión Fija de Gastos Básicos
@testset "InflationFixedExclusionCPI" begin
    # Creación de vectores de exclusión de prueba
    exc00 = [10, 100, 200, 218]
    exc10 = [20, 120, 220, 279]
    # Instanciar un tipo
    simplefn = InflationFixedExclusionCPI(exc00, exc10)
    @test simplefn isa InflationFixedExclusionCPI

    # Probar que esté definido el método para obtener su nombre
    @test measure_name(simplefn) isa String
    @test measure_tag(simplefn) isa String

    # Probar con bases del IPC con variaciones intermensuales iguales a cero.
    # Estas pruebas ayudan a verificar que la función de inflación se pueda
    # llamar sobre los tipos correctos

    zero_base = getzerobase()

    m_traj_infl = simplefn(zero_base,1)
    # Probamos que el resumen intermensual sea igual a cero
    @test all(isapprox.(m_traj_infl, 0; atol = 0.0001))

    # Obtenemos un UniformCountryStructure con dos bases y todas las variaciones
    # intermensuales iguales a cero
    zero_cst = getzerocountryst()
    traj_infl = simplefn(zero_cst)

    # Probamos que la trayectoria de inflación sea más larga que el resumen
    # intermensual de una sola base
    @test length(traj_infl) > length(m_traj_infl)

    @test all(isapprox.(traj_infl, 0; atol = 0.0001))

end

# Función de inflación por percentiles equiponderados
# prueba con percentil 70
@testset "InflationPercentileEq" begin
    # Instanciar un tipo
    percEqfn = InflationPercentileEq(70)
    @test percEqfn isa InflationPercentileEq

    # Probar que esté definido el método para obtener su nombre
    @test measure_name(percEqfn) isa String
    @test measure_tag(percEqfn) isa String

    # Probar con bases del IPC con variaciones intermensuales iguales a cero.
    # Estas pruebas ayudan a verificar que la función de inflación se pueda
    # llamar sobre los tipos correctos

    zero_base = getzerobase()
    m_traj_infl = percEqfn(zero_base)
    # Probamos que el resumen intermensual sea igual a cero
    @test all(m_traj_infl .≈ 0)

    # Obtenemos un UniformCountryStructure con dos bases y todas las variaciones
    # intermensuales iguales a cero
    zero_cst = getzerocountryst()
    traj_infl = percEqfn(zero_cst)

    # Probamos que la trayectoria de inflación sea más larga que el resumen
    # intermensual de una sola base
    @test length(traj_infl) > length(m_traj_infl)

    # Probamos que la trayectoria de inflación sea igual a cero
    @test all(traj_infl .≈ 0)
end

# Función de inflación por percentiles ponderados
# prueba con percentil 70
@testset "InflationPercentileWeighted" begin
    # Instanciar un tipo
    percfn = InflationPercentileWeighted(70)
    @test percfn isa InflationPercentileWeighted

    # Probar que esté definido el método para obtener su nombre
    @test measure_name(percfn) isa String
    @test measure_tag(percfn) isa String

    # Probar con bases del IPC con variaciones intermensuales iguales a cero.
    # Estas pruebas ayudan a verificar que la función de inflación se pueda
    # llamar sobre los tipos correctos

    zero_base = getzerobase()
    m_traj_infl = percfn(zero_base)
    # Probamos que el resumen intermensual sea igual a cero
    @test all(m_traj_infl .≈ 0)

    # Obtenemos un UniformCountryStructure con dos bases y todas las variaciones
    # intermensuales iguales a cero
    zero_cst = getzerocountryst()
    traj_infl = percfn(zero_cst)

    # Probamos que la trayectoria de inflación sea más larga que el resumen
    # intermensual de una sola base
    @test length(traj_infl) > length(m_traj_infl)

    # Probamos que la trayectoria de inflación sea igual a cero
    @test all(traj_infl .≈ 0)
end

@testset "InflationDynamicExclusion" begin
    # Instanciar un tipo
    dynExfn = InflationDynamicExclusion(2,2)
    @test dynExfn isa InflationDynamicExclusion

    # Probar que esté definido el método para obtener su nombre
    @test measure_name(dynExfn) isa String
    @test measure_tag(dynExfn) isa String

    # Probar con bases del IPC con variaciones intermensuales iguales a cero.
    # Estas pruebas ayudan a verificar que la función de inflación se pueda
    # llamar sobre los tipos correctos

    zero_base = getzerobase()
    m_traj_infl = dynExfn(zero_base)
    # Probamos que el resumen intermensual sea igual a cero
    @test all(m_traj_infl .≈ 0)

    # Obtenemos un UniformCountryStructure con dos bases y todas las variaciones
    # intermensuales iguales a cero
    zero_cst = getzerocountryst()
    traj_infl = dynExfn(zero_cst)

    # Probamos que la trayectoria de inflación sea más larga que el resumen
    # intermensual de una sola base
    @test length(traj_infl) > length(m_traj_infl)

    # Probamos que la trayectoria de inflación sea igual a cero
    @test all(traj_infl .≈ 0)
end


text:
export Tile, render

import Base: convert, writemime

@doc """
A `Tile` is the basic currency in Escher.
Most of the functions in the Escher API take `Tile`s
among other things as arguments, and return a `Tile` as the result.

Tiles are immutable: once created there is no way to mutate them.
""" ->
abstract Tile

render{T <: Tile}(x::T, state) =
    error("$T cannot be rendered.")

immutable Leaf <: Tile
    element::Elem
end
render(x::Elem, state) = x
render(x::Leaf, state) = x.element

convert(::Type{Tile}, x::AbstractString) = Leaf(Elem(:span, x))
convert(::Type{Tile}, x::Char) = Leaf(Elem(:span, string(x)))
convert{ns, tag}(::Type{Tile}, x::Elem{ns, tag}) = Leaf(x)

bestmime(val) = begin
  for mime in ("text/html", "image/svg+xml", "image/png", "text/plain")
    mimewritable(mime, val) && return MIME(symbol(mime))
  end
  error("Cannot render $val.")
end

render_fallback(m::MIME"text/plain", x) = Elem(:div, stringmime(m, x))
render_fallback(m::MIME"text/html", x) = Elem(:div, innerHTML=stringmime(m, x))
render_fallback(m::MIME"image/svg+xml", x) = Elem(:div, innerHTML=stringmime(m, x))
render_fallback(m::MIME"image/png", x) =
    Elem(:img, src="data:image/png;base64," * stringmime(m, x))

render(x::AbstractFloat, state) = render((@sprintf "%0.3f" x), state)
render(x::Symbol, state) = render(string(x), state)
render(x::AbstractString, state) = Elem(:span, x)

render{T}(x::T, state) =
    try
        render(convert(Tile, x), state)
    catch err
        if !(isa(err, MethodError) &&
             err.f === convert &&
             err.args[1] === Tile)
             rethrow(err)
        end
        render_fallback(bestmime(x), x)
    end

@doc """
`Empty` is handy tile that is well... Empty.
use `empty` constant exported by Escher in your code.
""" ->
immutable Empty <: Tile
end

@doc """
An Empty element.
""" ->
const empty = Empty()

render(t::Empty, state) = Elem(:div)

#  writemime(io::IO, m::MIME"text/html", x::Tile) =
#      writemime(io, m, Elem(:div, Escher.render(x, Dict()), className="escherRoot"))
#
#  writemime{T <: Tile}(io::IO, m::MIME"text/html", x::Signal{T}) =
#      writemime(io, m, map(Escher.render, Patchwork.Elem, x))

# Note a TileList is NOT a Tile
immutable TileList
    tiles::AbstractArray
end

convert(::Type{TileList}, xs::AbstractArray) =
    TileList(xs)
convert(::Type{TileList}, xs::Tuple) =
    TileList([x for x in xs])
convert(::Type{TileList}, x::TileList) = x
convert(::Type{TileList}, x) =
    TileList([x])

render(t::TileList, state) =
    map(x -> render(x, state), t.tiles)

render(t::TileList, wrap, state) =
    Elem(wrap, Elem[render(x, state) for x in t.tiles])

"""
SignalWrap lets you pretend that a signal of tiles is also a tile
"""
immutable SignalWrap <: Tile
    signal::Signal
end

convert(::Type{Tile}, x::Signal) = SignalWrap(x)

render(tile::SignalWrap, state) = begin
    id = "signal-" * makeid(tile.signal)
    state["embedded_signals"][id] = tile.signal
    Elem("signal-container", signalId=id)
end


text:
#
# SurfaceNets in Julia from Meshing.jl
#
# Ported from the Javascript implementation by Mikola Lysenko (C) 2012
# https://github.com/mikolalysenko/mikolalysenko.github.com/blob/master/Isosurface/js/surfacenets.js
# MIT License
#
# Based on: S.F. Gibson, "Constrained Elastic Surface Nets". (1998) MERL Tech Report.
#

# Precompute edge table, like Paul Bourke does.
# This saves a bit of time when computing the centroid of each boundary cell
#

const cube_edges = (0, 1, 0, 2, 0, 4, 1, 3, 1, 5, 2, 3, 2, 6, 3, 7, 4, 5, 4, 6, 5, 7, 6, 7)
const sn_edge_table = [0, 7, 25, 30, 98, 101, 123, 124, 168, 175, 177, 182, 202, 205, 211,
                       212, 772, 771, 797, 794, 870, 865, 895, 888, 940, 939, 949, 946, 974,
                       969, 983, 976, 1296, 1303, 1289, 1294, 1394, 1397, 1387, 1388, 1464,
                       1471, 1441, 1446, 1498, 1501, 1475, 1476, 1556, 1555, 1549, 1546,
                       1654, 1649, 1647, 1640, 1724, 1723, 1701, 1698, 1758, 1753, 1735,
                       1728, 2624, 2631, 2649, 2654, 2594, 2597, 2619, 2620, 2792, 2799,
                       2801, 2806, 2698, 2701, 2707, 2708, 2372, 2371, 2397, 2394, 2342,
                       2337, 2367, 2360, 2540, 2539, 2549, 2546, 2446, 2441, 2455, 2448,
                       3920, 3927, 3913, 3918, 3890, 3893, 3883, 3884, 4088, 4095, 4065,
                       4070, 3994, 3997, 3971, 3972, 3156, 3155, 3149, 3146, 3126, 3121,
                       3119, 3112, 3324, 3323, 3301, 3298, 3230, 3225, 3207, 3200, 3200,
                       3207, 3225, 3230, 3298, 3301, 3323, 3324, 3112, 3119, 3121, 3126,
                       3146, 3149, 3155, 3156, 3972, 3971, 3997, 3994, 4070, 4065, 4095,
                       4088, 3884, 3883, 3893, 3890, 3918, 3913, 3927, 3920, 2448, 2455,
                       2441, 2446, 2546, 2549, 2539, 2540, 2360, 2367, 2337, 2342, 2394,
                       2397, 2371, 2372, 2708, 2707, 2701, 2698, 2806, 2801, 2799, 2792,
                       2620, 2619, 2597, 2594, 2654, 2649, 2631, 2624, 1728, 1735, 1753,
                       1758, 1698, 1701, 1723, 1724, 1640, 1647, 1649, 1654, 1546, 1549,
                       1555, 1556, 1476, 1475, 1501, 1498, 1446, 1441, 1471, 1464, 1388,
                       1387, 1397, 1394, 1294, 1289, 1303, 1296, 976, 983, 969, 974, 946,
                       949, 939, 940, 888, 895, 865, 870, 794, 797, 771, 772, 212, 211, 205,
                       202, 182, 177, 175, 168, 124, 123, 101, 98, 30, 25, 7, 0]

"""
Generate a mesh using naive surface nets.
This takes the center of mass of the voxel as the vertex for each cube.
"""
function surface_nets(data::Vector{T}, dims, eps, scale, origin) where {T}
    vertices = Point{3,T}[]
    faces = Face{4,Int}[]

    sizehint!(vertices, ceil(Int, maximum(dims)^2 / 2))
    sizehint!(faces, ceil(Int, maximum(dims)^2 / 2))

    n = 0
    x = [0, 0, 0]
    R = Array{Int}([1, (dims[1] + 1), (dims[1] + 1) * (dims[2] + 1)])
    buf_no = 1

    buffer = fill(zero(Int), R[3] * 2)

    v = Vector{T}([0.0, 0.0, 0.0])

    #March over the voxel grid
    x[3] = 0
    @inbounds while x[3] < dims[3] - 1

        # m is the pointer into the buffer we are going to use.
        # This is slightly obtuse because javascript does not have good support for packed data structures, so we must use typed arrays :(
        # The contents of the buffer will be the indices of the vertices on the previous x/y slice of the volume
        m = 1 + (dims[1] + 1) * (1 + buf_no * (dims[2] + 1))

        x[2] = 0
        @inbounds while x[2] < dims[2] - 1
            x[1] = 0
            @inbounds while x[1] < dims[1] - 1

                # Read in 8 field values around this vertex and store them in an array
                # Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later
                mask = 0x00
                @inbounds grid = (data[n+1], data[n+2], data[n+dims[1]+1],
                                  data[n+dims[1]+2],
                                  data[n+dims[1]*2+1+dims[1]*(dims[2]-2)],
                                  data[n+dims[1]*2+2+dims[1]*(dims[2]-2)],
                                  data[n+dims[1]*3+1+dims[1]*(dims[2]-2)],
                                  data[n+dims[1]*3+2+dims[1]*(dims[2]-2)])

                signbit(grid[1]) && (mask |= 0x01)
                signbit(grid[2]) && (mask |= 0x02)
                signbit(grid[3]) && (mask |= 0x04)
                signbit(grid[4]) && (mask |= 0x08)
                signbit(grid[5]) && (mask |= 0x10)
                signbit(grid[6]) && (mask |= 0x20)
                signbit(grid[7]) && (mask |= 0x40)
                signbit(grid[8]) && (mask |= 0x80)

                # Check for early termination if cell does not intersect boundary
                if mask == 0x00 || mask == 0xff
                    x[1] += 1
                    n += 1
                    m += 1
                    continue
                end

                #Sum up edge intersections
                edge_mask = sn_edge_table[mask+1]
                v[1] = 0.0
                v[2] = 0.0
                v[3] = 0.0
                e_count = 0

                #For every edge of the cube...
                @inbounds for i in 0:11

                    #Use edge mask to check if it is crossed
                    if (edge_mask & (1 << i)) == 0
                        continue
                    end

                    #If it did, increment number of edge crossings
                    e_count += 1

                    #Now find the point of intersection
                    e0 = cube_edges[(i<<1)+1]       #Unpack vertices
                    e1 = cube_edges[(i<<1)+2]
                    g0 = grid[e0+1]                 #Unpack grid values
                    g1 = grid[e1+1]
                    t  = g0 - g1                 #Compute point of intersection
                    if abs(t) > eps
                        t = g0 / t
                    else
                        continue
                    end

                    #Interpolate vertices and add up intersections (this can be done without multiplying)
                    k = 1
                    for j in 1:3
                        a = e0 & k
                        b = e1 & k
                        (a != 0) && (v[j] += 1.0)
                        if a != b
                            v[j] += (a != 0 ? -t : t)
                        end
                        k <<= 1
                    end
                end # edge check

                #Now we just average the edge intersections and add them to coordinate
                s = 1.0 / e_count
                for i in 1:3
                    @inbounds v[i] = (x[i] + s * v[i]) * scale[i] + origin[i]
                end

                #Add vertex to buffer, store pointer to vertex index in buffer
                buffer[m+1] = length(vertices)
                push!(vertices, Point{3,T}(v[1], v[2], v[3]))

                #Now we need to add faces together, to do this we just loop over 3 basis components
                for i in 0:2
                    #The first three entries of the edge_mask count the crossings along the edge
                    if (edge_mask & (1 << i)) == 0
                        continue
                    end

                    # i = axes we are point along.  iu, iv = orthogonal axes
                    iu = (i + 1) % 3
                    iv = (i + 2) % 3

                    #If we are on a boundary, skip it
                    if (x[iu+1] == 0 || x[iv+1] == 0)
                        continue
                    end

                    #Otherwise, look up adjacent edges in buffer
                    du = R[iu+1]
                    dv = R[iv+1]

                    #Remember to flip orientation depending on the sign of the corner.
                    if (mask & 0x01) != 0x00
                        push!(faces,
                              Face{4,Int}(buffer[m+1] + 1, buffer[m-du+1] + 1,
                                          buffer[m-du-dv+1] + 1, buffer[m-dv+1] + 1))
                    else
                        push!(faces,
                              Face{4,Int}(buffer[m+1] + 1, buffer[m-dv+1] + 1,
                                          buffer[m-du-dv+1] + 1, buffer[m-du+1] + 1))
                    end
                end
                x[1] += 1
                n += 1
                m += 1
            end
            x[2] += 1
            n += 1
            m += 2
        end
        x[3] += 1
        n += dims[1]
        buf_no = xor(buf_no, 1)
        R[3] = -R[3]
    end
    #All done!  Return the result

    return vertices, faces # faces are quads, indexed to vertices
end


text:
"""
    findminimum(v, n)

FINDMINIMUM finds the first n smallest elements' indices.

# Input Arguments
- `v::Array{Float64}`: the candidate values for selection.
- `n::Int`: number of smallest elements for consideration.

# Output Argument
- `idx::Array{Int}`: n smallest elements' indices.

"""
function findminimum(v, n)
    idx = sortperm(v)[1:n]
    return idx
end

"""
    spike(i, N)

SPIKE gives the N-dim spike vector with i-th element equals 1.

# Input Arguments
- `i::Int`: index for one.
- `N::Int`: dimension of the target spike vector.

# Output Argument
- `v::Array{Float64}`: the N-dim spike vector with i-th element equals 1.

"""
function spike(i, N)
    v = zeros(N)
    v[i] = 1
    return v
end

"""
    characteristic(list, N)

CHARACTERISTIC gives the characteristic function in n-dim vector space with
    values of index in list equal to 1.

# Input Arguments
- `list::Array{Int}`: list of indices.
- `N::Int`: dimension of the target vector.

# Output Argument
- `v::Array{Float64}`: the n-dim characteristic vector with values of index in list equal to 1.

"""
function characteristic(list, N)
    v = zeros(N)
    v[list] .= 1.0
    return v
end

function χ(list, N)
    v = zeros(N)
    v[list] .= 1.0
    return v
end

"""
    heat_sol(f0,𝚽,Σ,t)

HEAT\\_SOL gives the solution of heat partial differential equation with initial condition u(⋅, 0) = f0

# Input Arguments
- `f0::Array{Float64}`: initial condition vector.
- `𝚽::Matrix{Float64}`: graph Laplacian eigenvectors, served as graph Fourier transform matrix
- `Σ::Array{Int}`: diagonal matrix of eigenvalues.
- `t::Float`: time elapse.

# Output Argument
- `u::Array{Float64}`: the solution vector at time t

"""
function heat_sol(f0, 𝚽, Σ, t)
    u = 𝚽 * (exp.(-t .* Σ) .* 𝚽' * f0)
    return u
end

"""
    freq_band_matrix(ls,n)

FREQ\\_BAND\\_MATRIX provides characteristic diagonal matrix, which is useful for spectral graph filters design.

# Input Arguments
- `ls::Array{Int}`: list of indices.
- `n::Int`: dimension of the target vector.

# Output Argument
- `D::Array{Float64}`: the zero/one diagonal matrix.

"""
function freq_band_matrix(ls, n)
    f = characteristic(list, n)
    return Diagonal(f)
end

"""
    scatter_gplot(X; marker = nothing, ms = 4, plotOrder = :normal, c = :viridis, subplot = 1)

SCATTER\\_GPLOT generates a scatter plot figure, which is for quick viewing of a graph signal.
SCATTER\\_GPLOT!(X; ...) adds a plot to `current` one.

# Input Arguments
- `X::Matrix{Float64}`: points locations, can be 2-dim or 3-dim.
- `marker::Array{Float64}`: default is `nothing`. Present different colors given
    different signal value at each node.
- `ms::Array{Float64}`: default is `4`. Present different node sizes given
    different signal value at each node.
- `shape::Symbol`: default is `:none`. Shape of the markers.
- `mswidth::Number`: default is `0`. Width of the marker stroke border.
- `msalpha::Number`: default is `nothing`. The opacity for the marker stroke.
- `plotOrder::Symbol`: default is `:normal`. Optional choices `:s2l` or `:l2s`, i.e.,
    plots from the smallest value of `marker` to the largest value or the other way around,
    or `:propabs`, ms is automatically set to be proportional to the absolute value
    of the marker values.
- `c::Symbol`: default is `:viridis`. Colors of the markers.
- `subgplot::Int`: default is `1`. The subplot index.

"""
function scatter_gplot(X; marker = nothing, ms = 4, shape = :none, mswidth = 0,
                          msalpha = nothing, plotOrder = :normal, c = :viridis, subplot = 1)
    N, dim = size(X)
    if !isnothing(marker)
        if size(marker) == (N,) || size(marker) == (N, 1)
            marker = marker[:]  # reshape N x 1 matrix to a vector of length N
        else
            error("marker only accepts a vector of length $(N) or a matrix of size $(N) x 1.")
        end
        if plotOrder == :normal
            idx = 1:N
        elseif plotOrder == :s2l
            idx = sortperm(marker)
        elseif plotOrder == :l2s
            idx = sortperm(marker, rev = true)
        elseif plotOrder == :propabs
            idx = sortperm(abs.(marker))
        else
            error("plotOrder only supports for :normal, :s2l, :l2s, or :propabs.")
        end
        X = X[idx, :]
        marker = marker[idx]
        if length(ms) > 1
            ms = ms[idx]
        end
        # if plotOrder is :propabs, set ms to be a vector prop to
        # abs.(marker) automatically.
        # The maximum marker size is set to be maximum(ms).
        if plotOrder == :propabs
            maxmarkersize = maximum(ms)
            markerabsmax = maximum(abs.(marker))
            ms = maxmarkersize .* log.(1 .+ abs.(marker)) ./ (log(1 + markerabsmax) + 10^3 * eps())
        end
    end
    if dim == 2
        scatter(X[:, 1], X[:, 2], marker_z = marker, ms = ms, shape = shape, c = c,
                legend = false, mswidth = mswidth, msalpha = msalpha, cbar = true,
                aspect_ratio = 1, grid = false, subplot = subplot)
    elseif dim == 3
        scatter(X[:, 1], X[:, 2], X[:, 3], marker_z = marker, ms = ms, shape = shape, c = c,
                legend = false, mswidth = mswidth, msalpha = msalpha, cbar = true,
                aspect_ratio = 1, grid = false, subplot = subplot)
    else
        error("Dimension Error: scatter_gplot only supports for 2-dim or 3-dim scatter plots.")
    end
end

function scatter_gplot!(X; marker = nothing, ms = 4, shape = :none, mswidth = 0,
                          msalpha = nothing, plotOrder = :normal, c = :viridis, subplot = 1)
    N, dim = size(X)
    if !isnothing(marker)
        if size(marker) == (N,) || size(marker) == (N, 1)
            marker = marker[:]  # reshape N x 1 matrix to a vector of length N
        else
            error("marker only accepts a vector of length $(N) or a matrix of size $(N) x 1.")
        end
        if plotOrder == :normal
            idx = 1:N
        elseif plotOrder == :s2l
            idx = sortperm(marker)
        elseif plotOrder == :l2s
            idx = sortperm(marker, rev = true)
        elseif plotOrder == :propabs
            idx = sortperm(abs.(marker))
        else
            error("plotOrder only supports for :normal, :s2l, :l2s, or :propabs.")
        end
        X = X[idx, :]
        marker = marker[idx]
        if length(ms) > 1
            ms = ms[idx]
        end
        # if plotOrder is :propabs, set ms to be a vector prop to
        # abs.(marker) automatically
        # The maximum marker size is set to be maximum(ms).
        if plotOrder == :propabs
            maxmarkersize = maximum(ms)
            markerabsmax = maximum(abs.(marker))
            ms = maxmarkersize .* log.(1 .+ abs.(marker)) ./ (log(1 + markerabsmax) + 10^3 * eps())
        end
    end
    if dim == 2
        scatter!(X[:, 1], X[:, 2], marker_z = marker, ms = ms, shape = shape, c = c,
                legend = false, mswidth = mswidth, msalpha = msalpha, cbar = true,
                aspect_ratio = 1, grid = false, subplot = subplot)
    elseif dim == 3
        scatter!(X[:, 1], X[:, 2], X[:, 3], marker_z = marker, ms = ms, shape = shape, c = c,
                legend = false, mswidth = mswidth, msalpha = msalpha, cbar = true,
                aspect_ratio = 1, grid = false, subplot = subplot)
    else
        error("Dimension Error: scatter_gplot! only supports for 2-dim or 3-dim scatter plots.")
    end
end

"""
    cat_plot(X; marker = nothing, ms = 4)

CAT\\_PLOT generates a scatter plot figure for cat example, which is for quick
viewing of a graph signal within a specific range (i.e., xlims, ylims, zlims).
CAT\\_PLOT!(X; ...) adds a plot to `current` one.

# Input Arguments
- `X::Matrix{Float64}`: 3-dim points.
- `marker::Array{Float64}`: default is nothing. Present different colors given
    different signal value at each node.
- `ms::Array{Float64}`: default is 4. Present different node sizes given
    different signal value at each node.

"""
function cat_plot(X; marker = nothing, ms = 4)
    scatter(X[:,1],X[:,2],X[:,3], marker_z = marker, ms = ms, c = :viridis,
        legend = false, cbar = true, aspect_ratio = 1, xlims = [-100, 100],
        ylims = [-100, 100], zlims = [-100, 100])
end

function cat_plot!(X; marker = nothing, ms = 4)
    scatter!(X[:,1],X[:,2],X[:,3], marker_z = marker, ms = ms, c = :viridis,
        legend = false, cbar = true, aspect_ratio = 1, xlims = [-100, 100],
        ylims = [-100, 100], zlims = [-100, 100])
end

"""
    approx_error_plot(DVEC::Vector{Vector{Float64}}; frac::Float64 = 0.50)

draw relative approx. error w.r.t. fraction of coefficients retained (FCR)

# Input Arguments
- `DVEC::Vector{Vector{Float64}}`: a list of expansion coefficients.
- `fraction::Float64`: default is 0.5.

"""
function approx_error_plot(DVEC::Vector{Vector{Float64}}; frac::Float64 = 0.50)
    plot(xaxis = "Fraction of Coefficients Retained",
            yaxis = "Relative Approximation Error")
    T = ["Laplacian", "HGLET", "Haar", "Walsh", "GHWT_c2f", "GHWT_f2c", "eGHWT",
            "PC-NGWP", "VM-NGWP"]
    L = [(:dashdot, :red), (:solid, :brown), (:dashdot, :orange),
            (:dashdot, :pink), (:solid, :gray), (:solid,:green),
            (:solid,:blue), (:solid,:purple), (:solid,:black)]
    for i = 1:length(DVEC)
        dvec = DVEC[i]
        N = length(dvec)
        dvec_norm = norm(dvec,2)
        dvec_sort = sort(dvec.^2) # the smallest first
        er = max.(sqrt.(reverse(cumsum(dvec_sort)))/dvec_norm, 1e-12)
        p = Int64(floor(frac*N)) + 1 # upper limit
        plot!(frac*(0:(p-1))/(p-1), er[1:p], yaxis=:log, xlims = (0.,frac),
                label = T[i], line = L[i], linewidth = 2, grid = false)
    end
end

"""
    getall_expansioncoeffs(G_Sig::GraphSig, GP_star::GraphPart, VM_NGWP::Array{Float64,3}, PC_NGWP::Array{Float64,3}, 𝚽::Matrix{Float64})

get all expansion coefficients of `f` via all methods in NGWP.jl and MTSG.jl

# Input Arguments
- `G_Sig::GraphSig`: GraphSig of the primal graph
- `GP_star::GraphPart`: GraphPart of the dual graph
- `VM_NGWP::Array{Float64,3}`: varimax NGWP.
- `PC_NGWP::Array{Float64,3}`: pair-clustering NGWP.
- `𝚽::Matrix{Float64}`: the matrix of graph Laplacian eigenvectors.

"""
function getall_expansioncoeffs(G_Sig::GraphSig, GP_star::GraphPart, VM_NGWP::Array{Float64,3}, PC_NGWP::Array{Float64,3}, 𝚽::Matrix{Float64})
    ############# VM_NGWP
    dmatrix_VM = ngwp_analysis(G_Sig, VM_NGWP)
    dvec_vm_ngwp, BS_vm_ngwp = ngwp_bestbasis(dmatrix_VM, GP_star)
    ############# PC_NGWP
    dmatrix_PC = ngwp_analysis(G_Sig, PC_NGWP)
    dvec_pc_ngwp, BS_pc_ngwp = ngwp_bestbasis(dmatrix_PC, GP_star)
    ############# Laplacian
    dvec_Laplacian = 𝚽' * G_Sig.f
    ############# plain Laplacian eigenvectors coefficients
    GP = partition_tree_fiedler(G_Sig)
    dmatrixH = HGLET_Analysis_All(G_Sig, GP)[1]
    dvec_hglet, BS_hglet, _ = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixH, cfspec = 1)
    dmatrix = ghwt_analysis!(G_Sig, GP = GP)
    ############# Haar
    BS_haar = bs_haar(GP)
    dvec_haar = dmatrix2dvec(dmatrix, GP, BS_haar)
    ############# Walsh
    BS_walsh = bs_walsh(GP)
    dvec_walsh = dmatrix2dvec(dmatrix, GP, BS_walsh)
    ############# GHWT_c2f
    dvec_c2f, BS_c2f = ghwt_c2f_bestbasis(dmatrix, GP)
    ############# GHWT_f2c
    dvec_f2c, BS_f2c = ghwt_f2c_bestbasis(dmatrix, GP)
    ############# eGHWT
    dvec_eghwt, BS_eghwt = ghwt_tf_bestbasis(dmatrix, GP)
    DVEC = [dvec_Laplacian[:], dvec_hglet[:], dvec_haar[:], dvec_walsh[:],
            dvec_c2f[:], dvec_f2c[:], dvec_eghwt[:], dvec_pc_ngwp[:],
            dvec_vm_ngwp[:]]
    return DVEC
end

using Clustering
"""
    spectral_clustering(𝚽, M)

SPECTRAL_CLUSTERING return M graph clusters, i.e., {Vₖ| k = 1,2,...,M}.

# Input Argument
- `𝚽::Matrix{Float64}`: the matrix of graph Laplacian eigenvectors.
- `M::Int`: the number of graph clusters.

# Output Argument
- `clusters::Vector{Vector{Int}}`: graph cluster indices.

"""
function spectral_clustering(𝚽, M)
    if M < 2
        return [1:size(𝚽, 1)]
    end
    cluster_indices = assignments(kmeans(𝚽[:, 2:M]', M))
    clusters = Vector{Vector{Int}}[]
    for k in 1:M
        push!(clusters, findall(cluster_indices .== k)[:])
    end
    return clusters
end

"""
    transform2D(X; s = 1, t = [0,0])

TRANSFORM2D dilate each point of `X` by scale s and translate by 2D vector t.
"""
function transform2D(X; s = 1, t = [0,0])
    X1 = X .* s
    X2 = zeros(size(X))
    for i in 1:size(X,1)
        X2[i,1] = X1[i,1] + t[1]
        X2[i,2] = X1[i,2] + t[2]
    end
    return X2
end

"""
    NN_rendering(X, Img_Mat)

NN\\_RENDERING generates a rendering signal at each point of `X` from the image
`Img_Mat` by nearest neighbor method.
"""
function NN_rendering(X, Img_Mat)
    N = size(X,1)
    f = zeros(N)
    for i in 1:N
        nn_x, nn_y = Int(round(X[i, 2])), Int(round(X[i, 1]))
        if nn_x < 1 || nn_x > size(Img_Mat, 2) || nn_y < 1 || nn_y > size(Img_Mat, 1)
            print("Error: pixel out of boundary!")
            return
        end
        f[i] = Img_Mat[nn_x, nn_y]
    end
    return f
end

"""
    Bilinear_rendering(X, Img_Mat)

NN\\_RENDERING generates a rendering signal at each point of `X` from the image
`Img_Mat` by bilinear interpolation method.
"""
function Bilinear_rendering(X, Img_Mat)
    N = size(X,1)
    f = zeros(N)
    for i in 1:N
        x1, x2, y1, y2 = Int(floor(X[i, 2])), Int(floor(X[i, 2])) + 1, Int(floor(X[i, 1])), Int(floor(X[i, 1])) + 1
        x, y = X[i,2], X[i,1]
        F = [Img_Mat[x1,y1] Img_Mat[x1,y2]
            Img_Mat[x2,y1] Img_Mat[x2,y2]]
        prod_res = 1/((x2 - x1) * (y2 - y1)) * [x2-x x-x1] * F * [y2-y y-y1]'
        f[i] = prod_res[1,1]
    end
    return f
end

"""
    dct1d(k, N)

DCT1D returns k-th 1D DCT basis vector in Rᴺ.

# Input Arguments
- `k::Int64`: ord of DCT basis vector. k = 1,2,...,N.
- `N::Int64`: vector dimension.

# Output Argument
- `φ::Array{Float64}`: k-th 1D DCT basis vector in Rᴺ. (k is 1-indexed)
"""
function dct1d(k, N)
    φ = [cos(π * (k - 1) * (l + 0.5) / N) for l = 0:(N - 1)]
    return φ ./ norm(φ, 2)
end

"""
    dct2d_basis(N1, N2)

DCT2D\\_BASIS returns 2D DCT basis vectors in [0,1] x [0,1] with N1-1 and N2-1 subintervals respectively.

# Input Arguments
- `N1::Int64`: number of nodes in x-axis.
- `N2::Int64`: number of nodes in y-axis.

# Output Argument
- `𝚽::Matrix{Float64}`: 2D DCT basis vectors.
"""
function dct2d_basis(N1, N2)
    N = N1 * N2
    𝚽 = zeros(N, N)
    ind = 1
    for i in 1:N1, j in 1:N2
        φ₁, φ₂ = dct1d(i, N1), dct1d(j, N2)
        φ = reshape(φ₁ * φ₂', N)
        𝚽[:, ind] = φ
        ind += 1
    end
    return 𝚽
end

"""
    alternating_numbers(n)

ALTERNATING\\_NUMBERS e.g., n = 5, returns [1,5,2,4,3]; n = 6, returns [1,6,2,5,3,4]

# Input Arguments
- `n::Int64`: number of nodes in x-axis.

# Output Argument
- `arr::Array{Int64}`: result array.
"""
function alternating_numbers(n)
    mid = Int(ceil(n / 2))
    arr1 = 1:mid
    arr2 = n:-1:(mid + 1)
    arr = Array{Int64}(zeros(n))
    p1, p2 = 1, 1
    for i = 1:n
        if i % 2 == 1
            arr[i] = arr1[p1]
            p1 += 1
        else
            arr[i] = arr2[p2]
            p2 += 1
        end
    end
    return arr
end

"""
    compute_SNR(f, g)

COMPUTE\\_SNR, g = f + ϵ, SNR = 20 * log10(norm(f)/norm(g-f)).

# Input Arguments
- `f::Array{Float64}`: original signal.
- `g::Array{Float64}`: noisy signal.

# Output Argument
- `SNR::Float64`: SNR value.
"""
function compute_SNR(f, g)
    SNR = 20 * log10(norm(f) / norm(g - f))
    return SNR
end

"""
    sort_wavelets(A; onlyByLoc = false)

sort A's column wavelet vectors based on their focused nodes' indices. flip
signs via the cross correlation.

# Input Arguments
- `A::Matrix{Float64}`: whose column vectors are wavelets.

# Output Argument
- `A::Matrix{Float64}`: a matrix with sorted and sign flipped column.
"""
function sort_wavelets(A; onlyByLoc = false)
    ord = findmax(abs.(A), dims = 1)[2][:]
    idx = sortperm([j[1] for j in ord])
    A = A[:, idx]
    if onlyByLoc
        return A
    end

    N = size(A, 1)
    sgn = ones(size(A, 2))
    mid = Int(round(size(A, 2) / 2))
    A[:, mid] .*= (maximum(A[:, mid]) > -minimum(A[:, mid])) * 2 - 1
    for i in 1:size(A, 2)
        cor_res = crosscor(A[:, mid], A[:, i], -Int(ceil(N / 2)):Int(floor(N / 2)))
        if maximum(cor_res) < -minimum(cor_res)
            sgn[i] = -1
        end
    end
    A = A * Diagonal(sgn)
    return A
end

"""
    standardize_eigenvectors!(𝚽::Matrix{Float64})

standardize the signs of the eigenvectors such that 1) the term with the largest
magnitude is positive; 2) if the max == -min, then make sure the first non-zero
entry of the eigenvector is positive.

# Input Arguments
- `𝚽::Matrix{Float64}`: matrix of graph Laplacian eigenvectors

"""
function standardize_eigenvectors!(𝚽::Matrix{Float64})
    N, nev = size(𝚽)
    tol = 10^3 * eps()
    for l in 1:nev
        𝛟max = maximum(𝚽[:, l])
        𝛟min = minimum(𝚽[:, l])
        if abs(𝛟max + 𝛟min) < tol
            row = 1
            standardized = false
            while !standardized
                if 𝚽[row, l] > tol
                    standardized = true
                elseif 𝚽[row, l] < -tol
                    𝚽[:, l] = -𝚽[:, l]
                else
                    row += 1
                end
            end
        elseif 𝛟max < -𝛟min
            𝚽[:, l] *= -1
        end
    end
end

"""
getall_expansioncoeffs2(G_Sig::GraphSig, GP_star::GraphPart,
                                 GP_star_Lsym::GraphPart,
                                 VM_NGWP::Array{Float64,3},
                                 PC_NGWP::Array{Float64,3},
                                 LP_NGWP::Array{Float64,3},
                                 VM_NGWP_Lsym::Array{Float64,3},
                                 PC_NGWP_Lsym::Array{Float64,3},
                                 LP_NGWP_Lsym::Array{Float64,3},
                                 𝚽::Matrix{Float64},
                                 𝚽sym::Matrix{Float64})

get all expansion coefficients of `f` in dissertation via all methods in NGWP.jl
and MTSG.jl

# Input Arguments
- `G_Sig::GraphSig`: GraphSig of the primal graph
- `GP_star::GraphPart`: GraphPart of the dual graph
- `VM_NGWP::Array{Float64,3}`: varimax NGWP.
- `PC_NGWP::Array{Float64,3}`: pair-clustering NGWP.
- `LP_NGWP::Array{Float64,3}`: lapped NGWP.
- `VM_NGWP_Lsym::Array{Float64,3}`: Lsym version of varimax NGWP.
- `PC_NGWP_Lsym::Array{Float64,3}`: Lsym version of pair-clustering NGWP.
- `LP_NGWP_Lsym::Array{Float64,3}`: Lsym version of lapped NGWP.
- `𝚽::Matrix{Float64}`: the matrix of `L` eigenvectors.
- `𝚽sym::Matrix{Float64}`: the matrix of `Lsym` eigenvectors.

"""
function getall_expansioncoeffs2(G_Sig::GraphSig, GP_star::GraphPart,
                                 GP_star_Lsym::GraphPart,
                                 VM_NGWP::Array{Float64,3},
                                 PC_NGWP::Array{Float64,3},
                                 LP_NGWP::Array{Float64,3},
                                 VM_NGWP_Lsym::Array{Float64,3},
                                 PC_NGWP_Lsym::Array{Float64,3},
                                 LP_NGWP_Lsym::Array{Float64,3},
                                 𝚽::Matrix{Float64},
                                 𝚽sym::Matrix{Float64})
    ############# VM_NGWP
    dmatrix_VM = ngwp_analysis(G_Sig, VM_NGWP)
    dvec_vm_ngwp, BS_vm_ngwp = ngwp_bestbasis(dmatrix_VM, GP_star)
    ############# PC_NGWP
    dmatrix_PC = ngwp_analysis(G_Sig, PC_NGWP)
    dvec_pc_ngwp, BS_pc_ngwp = ngwp_bestbasis(dmatrix_PC, GP_star)
    ############# LP_NGWP
    dmatrix_LP = ngwp_analysis(G_Sig, LP_NGWP)
    dvec_lp_ngwp, BS_lp_ngwp = ngwp_bestbasis(dmatrix_LP, GP_star)
    ############# VM_NGWP_Lsym
    dmatrix_VM_Lsym = ngwp_analysis(G_Sig, VM_NGWP_Lsym)
    dvec_vm_ngwp_Lsym, BS_vm_ngwp_Lsym = ngwp_bestbasis(dmatrix_VM_Lsym, GP_star_Lsym)
    ############# PC_NGWP_Lsym
    dmatrix_PC_Lsym = ngwp_analysis(G_Sig, PC_NGWP_Lsym)
    dvec_pc_ngwp_Lsym, BS_pc_ngwp_Lsym = ngwp_bestbasis(dmatrix_PC_Lsym, GP_star_Lsym)
    ############# LP_NGWP
    dmatrix_LP_Lsym = ngwp_analysis(G_Sig, LP_NGWP_Lsym)
    dvec_lp_ngwp_Lsym, BS_lp_ngwp_Lsym = ngwp_bestbasis(dmatrix_LP_Lsym, GP_star_Lsym)
    ############# Laplacian L
    dvec_Laplacian = 𝚽' * G_Sig.f
    ############# Laplacian Lsym
    dvec_Laplacian_sym = 𝚽sym' * G_Sig.f
    ############# HGLET
    GP = partition_tree_fiedler(G_Sig)
    dmatrixH, _, dmatrixHsym = HGLET_Analysis_All(G_Sig, GP)
    dvec_hglet, BS_hglet, trans_hglet = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixH, dmatrixHsym = dmatrixHsym, cfspec = 1)
    ############# LP-HGLET
    dmatrixsH, dmatrixsHsym = LPHGLET_Analysis_All(G_Sig, GP; ϵ = 0.3)
    dvec_lphglet, BS_lphglet, trans_lphglet = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixsH, dmatrixHsym = dmatrixsHsym, cfspec = 1)
    ############# GHWT dictionaries
    dmatrix = ghwt_analysis!(G_Sig, GP = GP)
    ############# Haar
    BS_haar = bs_haar(GP)
    dvec_haar = dmatrix2dvec(dmatrix, GP, BS_haar)
    ############# Walsh
    BS_walsh = bs_walsh(GP)
    dvec_walsh = dmatrix2dvec(dmatrix, GP, BS_walsh)
    ############# GHWT_c2f
    dvec_c2f, BS_c2f = ghwt_c2f_bestbasis(dmatrix, GP)
    ############# GHWT_f2c
    dvec_f2c, BS_f2c = ghwt_f2c_bestbasis(dmatrix, GP)
    ############# eGHWT
    dvec_eghwt, BS_eghwt = ghwt_tf_bestbasis(dmatrix, GP)
    DVEC = [dvec_Laplacian[:], dvec_Laplacian_sym[:], dvec_hglet[:], dvec_lphglet[:],
            dvec_haar[:], dvec_walsh[:], dvec_c2f[:], dvec_f2c[:], dvec_eghwt[:],
            dvec_pc_ngwp[:], dvec_pc_ngwp_Lsym[:],
            dvec_vm_ngwp[:], dvec_vm_ngwp_Lsym[:],
            dvec_lp_ngwp[:], dvec_lp_ngwp_Lsym[:]]
    return DVEC
end

"""
    approx_error_plot2(DVEC::Vector{Vector{Float64}}; frac::Float64 = 0.50)

draw relative approx. error w.r.t. fraction of coefficients retained (FCR)
in dissertation

# Input Arguments
- `DVEC::Vector{Vector{Float64}}`: a list of expansion coefficients.
- `fraction::Float64`: default is 0.5.

"""
function approx_error_plot2(DVEC::Vector{Vector{Float64}}; frac::Float64 = 0.50)
    plot(xaxis = "Fraction of Coefficients Retained",
            yaxis = "Relative Approximation Error", size = (600, 500))
    T = ["eigenbasis-L", "eigenbasis-Lsym", "HGLET", "LP-HGLET", "Haar", "Walsh",
         "GHWT_c2f", "GHWT_f2c", "eGHWT", "PC-NGWP", "PC-NGWP-Lsym",
         "VM-NGWP", "VM-NGWP-Lsym", "LP-NGWP", "LP-NGWP-Lsym"]
    L = [(:dot, :red), (:dot, :magenta), (:solid, :teal), (:dashdot, :teal),
         (:dashdot, :orange), (:dashdot, :pink), (:solid, :gray),
         (:solid, :green), (:solid, :blue), (:solid, :purple), (:dash, :purple),
         (:solid, :black), (:dash, :black), (:solid, :orange), (:dash, :orange)]
    for i = 1:length(DVEC)
        if i in [5, 6, 7, 8]
            continue
        end
        dvec = DVEC[i]
        N = length(dvec)
        dvec_norm = norm(dvec,2)
        dvec_sort = sort(dvec.^2) # the smallest first
        er = max.(sqrt.(reverse(cumsum(dvec_sort)))/dvec_norm, 1e-12)
        p = Int64(floor(frac*N)) + 1 # upper limit
        plot!(frac*(0:(p-1))/(p-1), er[1:p], yaxis=:log, xlims = (0.,frac),
                label = T[i], line = L[i], linewidth = 2, grid = false)
    end
end

"""
    stem(f; x = 1:length(f[:]), c = :black, ms = 5, shape = :circle, lw = 1)
Display a length N discrete sequence data in a stem plot.

# Input Arguments
- `f`: function values could be a vector or an N x 1 matrix.
- `x`: default is `1:N`. x values.
- `c::Symbol`: default is `:black`. Color of the markers.
- `ms::Number`: default is `5`. Size of markers.
- `shape::Symbol`: default is `:circle`. Shape of the markers.
- `lw::Number`: default is `1`. The line width of the stems.
- `subgplot::Int`: default is `1`. The subplot index.
"""
function stem(f; x = 1:length(f[:]), c = :black, ms = 5, shape = :circle, lw = 1, subplot = 1)
    if (ndims(f) == 2 && size(f, 2) > 1) || ndims(f) > 2
        error("input only accepts a vector or a matrix of size (N, 1).")
        return
    end
    f = f[:]
    N = length(f)
    plot(x, f, seriestype = :stem, color = c, lw = lw, subplot = subplot)
    scatter!(x, f, marker = shape, color = c, ms = ms, mswidth = 0, subplot = subplot)
    plot!(x, zeros(N), color = :black, lw = lw, subplot = subplot)
    plot!(grid = false, legend = false, ratio = 1, frame = :none, subplot = subplot)
end

function stem!(f; x = 1:length(f[:]), c = :black, ms = 5, shape = :circle, lw = 1, subplot = 1)
    if (ndims(f) == 2 && size(f, 2) > 1) || ndims(f) > 2
        error("input only accepts a vector or a matrix of size (N, 1).")
        return
    end
    f = f[:]
    N = length(f)
    plot!(x, f, seriestype = :stem, color = c, lw = lw, subplot = subplot)
    scatter!(x, f, marker = shape, color = c, ms = ms, mswidth = 0, subplot = subplot)
    plot!(x, zeros(N), color = :black, lw = lw, subplot = subplot)
    plot!(grid = false, legend = false, ratio = 1, frame = :none, subplot = subplot)
end


text:
function find_α(params, α_min, α_max, dα)
    α = Float64[]
    n = Int[]

    α₀ = [0.0]
    j = 0
    while !isempty(α₀) && α₀[1] < α_max
        println("Computing zeros of Fₙ(α) for n = $j")
        α₀ = find_zeros(α -> α_func(params, α, j), α_min, α_max)
        append!(α, α₀)
        append!(n, fill(j, size(α₀)))
        j = j + 1
    end

    inds = sortperm(α)
    α = α[inds]
    n = n[inds]

    # Number of layers
    L = length(params.D)

    inds_keep = dα / 100 .< α .< α_max
    α = α[inds_keep]
    n = n[inds_keep]

    if params.W[L] < 1e-12 && (length(params.W) == L || params.W[end] < 1e-12)
        # Add ground eigenmode explicitly
        pushfirst!(α, 0)
        pushfirst!(n, 0)
    end

    α, n
end


text:
@testset "Retirement tests" begin

mpSim = ManpowerSimulation( "sim" )
setSimulationLength!( mpSim, 300 )

attribute = Attribute( "attribute A" )
addInitialAttributeValue!( attribute, "value A", 1 )
addInitialAttributeValue!( attribute, "value B", 1 )
addSimulationAttribute!( mpSim, attribute )

node = BaseNode( "node A" )
addNodeRequirement!( node, "attribute A", "value A" )
addSimulationBaseNode!( mpSim, node )

@testset "Retirement at tenure" begin
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 12 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentFixed!( recruitment, 10 )
    addSimulationRecruitment!( mpSim, recruitment )

    retirement = Retirement()
    setRetirementCareerLength!( retirement, 120 )
    setRetirementSchedule!( retirement, 12 )
    setSimulationRetirement!( mpSim, retirement )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig = false )
    report = nodeFluxReport( mpSim, 12, :out, "active" )["active"]
    @test all( report[:, Symbol( "active => external" )] .==
        vcat( zeros( 10 ), fill( 10, 16 ) ) ) &&
        all( report[:, Symbol( "retirement: node A => external" )] .==
        vcat( zeros( 10 ), fill( 10, 16 ) ) )
end  # @testset "Retirement at tenure"

@testset "Retirement schedule" begin
    retirement = Retirement()
    setRetirementCareerLength!( retirement, 120 )
    setRetirementSchedule!( retirement, 24 )
    setSimulationRetirement!( mpSim, retirement )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig = false )
    report = nodeFluxReport( mpSim, 12, :out, "active" )["active"]
    @test all( report[:, Symbol( "active => external" )] .==
        vcat( zeros( 10 ), 10, 0, repeat( [20, 0], 7 ) ) )

    retirement = Retirement()
    setRetirementCareerLength!( retirement, 120 )
    setRetirementSchedule!( retirement, 24, 12 )
    setSimulationRetirement!( mpSim, retirement )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig = false )
    report = nodeFluxReport( mpSim, 12, :out, "active" )["active"]
    @test all( report[:, Symbol( "active => external" )] .==
        vcat( zeros( 10 ), repeat( [0, 20], 8 ) ) )
end  # @testset "Retirement schedule"

@testset "Retirement at age" begin
    clearSimulationRecruitment!( mpSim )
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 12 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentFixed!( recruitment, 10 )
    setRecruitmentAgeFixed!( recruitment, 60 )
    addSimulationRecruitment!( mpSim, recruitment )

    retirement = Retirement()
    setRetirementAge!( retirement, 180 )
    setRetirementSchedule!( retirement, 12 )
    setSimulationRetirement!( mpSim, retirement )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig = false )
    report = nodeFluxReport( mpSim, 12, :out, "active" )["active"]
    @test all( report[:, Symbol( "active => external" )] .==
        vcat( zeros( 10 ), fill( 10, 16 ) ) )
end  # @testset "Retirement at age"

@testset "Retirement at age/tenure" begin
    node = BaseNode( "node B" )
    addNodeRequirement!( node, "attribute A", "value B" )
    addSimulationBaseNode!( mpSim, node )

    clearSimulationRecruitment!( mpSim )
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 12 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentFixed!( recruitment, 10 )
    addSimulationRecruitment!( mpSim, recruitment )

    recruitment = Recruitment( "recruitment B" )
    setRecruitmentSchedule!( recruitment, 12 )
    setRecruitmentTarget!( recruitment, "node B" )
    setRecruitmentFixed!( recruitment, 10 )
    setRecruitmentAgeFixed!( recruitment, 60 )
    addSimulationRecruitment!( mpSim, recruitment )

    retirement = Retirement()
    setRetirementAge!( retirement, 144 )
    setRetirementCareerLength!( retirement, 120 )
    setRetirementSchedule!( retirement, 12 )
    setSimulationRetirement!( mpSim, retirement )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig = false )
    report = nodeFluxReport( mpSim, 12, :out, "active" )["active"]
    @test all( report[:, Symbol( "retirement: node A => external" )] .==
        vcat( fill( 0, 10 ), fill( 10, 16 ) ) ) &&
        all( report[:, Symbol( "retirement: node B => external" )] .==
        vcat( fill( 0, 7 ), fill( 10, 19 ) ) )

    retirement = Retirement()
    setRetirementAge!( retirement, 144 )
    setRetirementCareerLength!( retirement, 120 )
    setRetirementIsEither!( retirement, false )
    setRetirementSchedule!( retirement, 12 )
    setSimulationRetirement!( mpSim, retirement )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig = false )
    report = nodeFluxReport( mpSim, 12, :out, "active" )["active"]
    @test all( report[:, Symbol( "retirement: node A => external" )] .==
        vcat( fill( 0, 12 ), fill( 10, 14 ) ) ) &&
        all( report[:, Symbol( "retirement: node B => external" )] .==
        vcat( fill( 0, 10 ), fill( 10, 16 ) ) )
end  # @testset "Retirement at age/tenure"

end  # @testset "Retirement tests"


text:
using Parameters
using Statistics
using Distributions
using Plots
using StatsPlots
using CSV
using DataFrames

# Standard library
using Random

#using .BotMarket

function save_positions(investors, path, assets=nothing)
    if isnothing(assets)
        assets = [
            :generic_currency,
            :generic_stock,
        ]
    end

    df = DataFrame(
        name = [inv.name for inv in investors],
    )
    for asset in assets
        df[:, string(asset) * "_position"] = [get(inv.positions, eval(asset), 0) for inv in investors]
        df[:, string(asset) * "_reserve"] = [get(inv.reserved, eval(asset), 0) for inv in investors]
    end
    CSV.write(path, df) # "results/pos_start.csv"
    df
end

function save_trades(trades, trd_day_of_trade)
    df = DataFrame(
        timestamp = [trade.timestamp for trade in trades],
        trading_day =  trd_day_of_trade,
        price = [trade.price for trade in trades],
        quantity = [trade.quantity for trade in trades],
        seller = [trade.seller.name for trade in trades],
        buyer = [trade.buyer.name for trade in trades],
    )
    CSV.write("results/trades.csv", df)
    df
end

function save_trades(trades, trd_day_of_trade; path)
    df_main = DataFrame()
    for market in keys(trades)

        df = DataFrame(
            market = market.name,
            from = market.currency.name,
            to = market.traded_asset.name,
            timestamp = [trade.timestamp for trade in trades[market]],
            price = [trade.price for trade in trades[market]],
            quantity = [trade.quantity for trade in trades[market]],
            seller = [trade.seller.name for trade in trades[market]],
            buyer = [trade.buyer.name for trade in trades[market]],
        )
        df_main = vcat(df_main, df)
    end
    sort!(df_main, [:timestamp]);
    begin
        df_main[!, :trading_day] .= trd_day_of_trade
    end
    CSV.write(path, df_main)
    df_main
end

function get_market_books(markets::Array{T, 1} where {T<:AbstractMarket})
    df_books = DataFrame()
    for market in markets
        for (side, book) in (("bid", market.bid_limit_orders), ("ask", market.ask_limit_orders))
            if isempty(book)
                continue
            end
            df = DataFrame(
                [(side, order.dealer.name, order.from.name, order.to.name, order.price, order.quantity, market.name) for order in book],
            )
            rename!(df, [:side, :trader, :from, :to, :price, :quantity, :market])
            df_books = vcat(df_books, df)
        end
    end
    return df_books
end

function run_price_discovery(investors::Array{T, 1} where {T<:AbstractInvestor}, markets::Array{T, 1} where {T<:AbstractMarket})

    # Placing
    for investor in investors, market in markets

        place!(investor, market)
    end

    # Clearing
    trades = Dict{AbstractMarket, Array{Trade, 1}}(market => [] for market in markets)
    for market in markets
        trds = clear!(market)
        if ~isempty(trds)
            append!(trades[market], trds)
        end
    end
    return trades
end

function run_day(investors::Array{T, 1} where {T<:AbstractInvestor}, markets::Array{T, 1} where {T<:AbstractMarket})
    trades = Dict{AbstractMarket, Array{Trade, 1}}(market => [] for market in markets)
    for investor in investors, market in markets[shuffle(1:end)]

        place!(investor, market)
        trds = clear!(market)
        if ~isempty(trds)
            append!(trades[market], trds)
        end
    end

    return trades
end

function run_simulation(
        investors::Array{T, 1} where {T<:AbstractInvestor},
        markets::Array{T, 1} where {T<:AbstractMarket},
        assets::Array{T, 1} where {T<:AbstractAsset},
        trading_days::Int64, speak_ratio::Float64,
        name::String
    )
    save_positions(investors, "results/$name/pos_start.csv")
    trades = Dict{AbstractMarket, Array{Trade, 1}}(market => [] for market in markets)
    df_books = DataFrame()

    day_index = []
    for day in 1:trading_days
        # Trading day
        #println("Trading day ", day)

        session_trades = run_day(sample(investors, Int64(round(length(investors) * speak_ratio)), replace = false), markets)
        pay_cashflows!(investors, assets)
        update!(world)

        for market in keys(session_trades)
            trades[market] = vcat(trades[market], session_trades[market])
            end_price = length(session_trades[market]) > 0 ? session_trades[market][end].price : "<no trades>"
            println("Price of the day $day: $(end_price) ($(market.name))")
        end
        n_trades = sum([length(val) for val in values(session_trades)])
        day_index = vcat(day_index, repeat([day], n_trades))

        # Order book and its content
        df_book = get_market_books(markets)
        if ~isempty(df_book)
            df_book[!, :day] .= day
            df_books = vcat(df_books, df_book)
        end
    end

    save_positions(investors, "results/$name/pos_end.csv")

    save_trades(trades, day_index, path="results/$name/trades.csv")
    CSV.write("results/$name/books.csv", df_books)
end


text:
include("definitions.jl")

#Learn variational parameters for theta from one document, leaving topics fixed.
#Non-sparse version, for reference, not used.
function learnThetaForDoc(VBparams::collapsed_VB_params_type, documents::Array{Document,1}, docInd::Int, burnInPerDoc::Int, tau2::Float64, kappa2::Float64, meanchangeThreshold::Float64)
    oldDTC = VBparams.documentTopicCounts[docInd,:];
    docLength = length(documents[docInd].wordVector);
    probs::Array{Float64,1} = Array(Float64, length(VBparams.alpha))

    eta = VBparams.eta
    alpha = VBparams.alpha
    topicCounts = VBparams.topicCounts
    etaSum = VBparams.etaSum
    wordTopicCounts = VBparams.wordTopicCounts
    documentTopicCounts = VBparams.documentTopicCounts
    K = length(alpha)

    for burn=1:burnInPerDoc
        for wordInd = 1:docLength
            word::Int = documents[docInd].wordVector[wordInd];
            stepSize2::Float64 = (wordInd + (burn-1) .* docLength + tau2)^-kappa2;

            for topic = 1:K
                probs[topic] = (wordTopicCounts[word, topic] + eta[word]) .*
                        (documentTopicCounts[docInd,topic] + alpha[topic]) ./
                        (topicCounts[topic] + etaSum);
            end

            probs /= sum(probs);
            #update params for theta
            for topic = 1:K
                documentTopicCounts[docInd,topic] = (1 - stepSize2) * documentTopicCounts[docInd,topic] + stepSize2 .* docLength .* probs[topic];
            end
        end
    end
    VBparams
end

#Learn variational parameters for theta and phi from one document
#Non-sparse version, for reference, not used.
function learnTopicsForDoc(VBparams::collapsed_VB_params_type, documents::Array{Document,1}, docInd, tau2, kappa2, wordTopicCounts_hat, topicCounts_hat)
    docLength = length(documents[docInd].wordVector);
    probs::Array{Float64,1} = Array(Float64, length(VBparams.alpha))

    eta = VBparams.eta
    alpha = VBparams.alpha
    topicCounts = VBparams.topicCounts
    etaSum = VBparams.etaSum
    wordTopicCounts = VBparams.wordTopicCounts
    documentTopicCounts = VBparams.documentTopicCounts
    K = length(alpha)

    for wordInd = 1:docLength
        word::Int = documents[docInd].wordVector[wordInd];
        stepSize2::Float64 = (wordInd + tau2)^-kappa2;

        for topic = 1:K
            probs[topic] = (wordTopicCounts[word, topic] + eta[word]) .*
                    (documentTopicCounts[docInd,topic] + alpha[topic]) ./
                    (topicCounts[topic] + etaSum);
        end

        probs /= sum(probs);

        #update sufficient statistics
        for topic = 1:K
            VBparams.documentTopicCounts[docInd,topic] = (1 - stepSize2) * documentTopicCounts[docInd,topic] + stepSize2 .* docLength .* probs[topic];
            wordTopicCounts_hat[word, topic] = wordTopicCounts_hat[word, topic] + probs[topic];
            topicCounts_hat[topic] += probs[topic];
        end
    end
    VBparams, wordTopicCounts_hat, topicCounts_hat
end

#Learn variational parameters for theta, leaving topics fixed.
#sparse, with clumping
function learnThetaForDoc_Sparse(VBparams, documents, docInd, burnInPerDoc, tau2, kappa2, meanchangeThreshold)
    oldDTC = VBparams.documentTopicCounts[docInd,:];
    docLength = documents[docInd].docLength;
    numDistinctTokens = length(documents[docInd].tokenIndex);
    stepSizeCounter = 1;
    stepSize2 = (stepSizeCounter + tau2)^-kappa2;

    probs::Array{Float64,1} = Array{Float64}(length(VBparams.alpha))

    eta = VBparams.eta
    alpha = VBparams.alpha
    topicCounts = VBparams.topicCounts
    etaSum = VBparams.etaSum
    wordTopicCounts = VBparams.wordTopicCounts
    documentTopicCounts = VBparams.documentTopicCounts
    K = length(alpha)

    for burn=1:burnInPerDoc
        for tokenInd = 1:numDistinctTokens
            word = documents[docInd].tokenIndex[tokenInd];

            for topic = 1:K
                probs[topic] = (wordTopicCounts[word, topic] + eta[word]) .*
                        (documentTopicCounts[docInd,topic] + alpha[topic]) ./
                        (topicCounts[topic] + etaSum);
            end

            sum_probs = sum(probs)

            #update params for theta
            tokenCountInDoc = documents[docInd].tokenCount[tokenInd];

            left_step = (1 - stepSize2)^tokenCountInDoc
            right_step = 1-(1-stepSize2)^tokenCountInDoc
            for topic = 1:K
                probs[topic] /= sum_probs;
                documentTopicCounts[docInd,topic] = left_step .* documentTopicCounts[docInd,topic] +
                                                    docLength .* probs[topic] .* right_step;
            end

            stepSizeCounter = stepSizeCounter + tokenCountInDoc;
            stepSize2 = (stepSizeCounter + (burn-1) .* docLength + tau2)^-kappa2;
        end
        meanchange = mean(abs.(VBparams.documentTopicCounts[docInd,:] - oldDTC));
        if meanchange < meanchangeThreshold
            #fprintf("converged after #d iterations\n", burn);
            break;
        end
        oldDTC = VBparams.documentTopicCounts[docInd,:];
    end
    VBparams
end

#Learn variational parameters for theta and phi
#sparse, with clumping
function updateBoth_SparseFull(VBparams, documents, docInd, tau2, kappa2, wordTopicCounts_hat, topicCounts_hat)
    docLength = documents[docInd].docLength;
    numDistinctTokens = length(documents[docInd].tokenIndex);
    stepSizeCounter = 1;
    stepSize2 = (stepSizeCounter + tau2)^-kappa2;

    probs::Array{Float64,1} = Array{Float64}(length(VBparams.alpha))

    eta = VBparams.eta
    alpha = VBparams.alpha
    topicCounts = VBparams.topicCounts
    etaSum = VBparams.etaSum
    wordTopicCounts = VBparams.wordTopicCounts
    documentTopicCounts = VBparams.documentTopicCounts
    K = length(alpha)

    for tokenInd = 1:numDistinctTokens
        word = documents[docInd].tokenIndex[tokenInd];

        #CVB0 estimate of gamma for current token
        for topic = 1:K
            probs[topic] = (wordTopicCounts[word, topic] + eta[word]) .*
                    (documentTopicCounts[docInd,topic] + alpha[topic]) ./
                    (topicCounts[topic] + etaSum);
        end

        sum_probs = sum(probs);

        tokenCountInDoc = documents[docInd].tokenCount[tokenInd];

        left_step = (1 - stepSize2)^tokenCountInDoc
        right_step = 1-(1-stepSize2)^tokenCountInDoc

        for topic = 1:K
            probs[topic] /= sum_probs
            documentTopicCounts[docInd,topic] = left_step .* documentTopicCounts[docInd,topic] + docLength .* probs[topic] .* right_step;
            wordTopicCounts_hat[word, topic] = wordTopicCounts_hat[word, topic] + probs[topic] * tokenCountInDoc;
            topicCounts_hat[topic] += probs[topic]  * tokenCountInDoc;
        end

        stepSizeCounter = stepSizeCounter + tokenCountInDoc;
        stepSize2 = (stepSizeCounter + tau2)^-kappa2;
    end
    VBparams, wordTopicCounts_hat, topicCounts_hat
end


text:
using Plots
using FFMPEG

module VisDyn

include("colorMap.jl")
export cMap, cMap_def, cNormalize, qcMap_def
export plotTrace!, plotScatter!, plotTrajectory!
include("plotGradient.jl")
export plotVF!, plotNullCline!, plotGradient!, plotValongCycle!, circularCycle

# Write your package code here.
end


text:
using NiLang, NiLang.AD

include("symlib.jl")
NiLang.AD.isvar(sym::Basic) = true
NiLang.AD.GVar(sym::Basic) = GVar(sym, zero(sym))

# a patch for symbolic IROT
@i @inline function NiLang.IROT(a!::GVar{<:Basic}, b!::GVar{<:Basic}, θ::GVar{<:Basic})
    IROT(a!.x, b!.x, θ.x)
    NEG(θ.x)
    θ.x -= Basic(π)/2
    ROT(a!.g, b!.g, θ.x)
    θ.g += a!.x * a!.g
    θ.g += b!.x * b!.g
    θ.x += Basic(π)/2
    NEG(θ.x)
    ROT(a!.g, b!.g, Basic(π)/2)
end

NiLang.INC(x::Basic) = x + one(x)
NiLang.DEC(x::Basic) = x - one(x)
@inline function NiLang.ROT(i::Basic, j::Basic, θ::Basic)
    a, b = rot(i, j, θ)
    a, b, θ
end
@inline function NiLang.IROT(i::Basic, j::Basic, θ::Basic)
    i, j, _ = ROT(i, j, -θ)
    i, j, θ
end
Base.sincos(x::Basic) = (sin(x), cos(x))

function printall()
    syms = [Basic(:a), Basic(:b), Basic(:c)]

    for (subop, nargs) in [(identity, 2), (*, 3), (/, 3), (^, 3), (exp, 2), (log, 2), (sin, 2), (cos, 2)]
        for opm in [PlusEq, MinusEq]
            op = opm(subop)
            @show op
            printone(op, syms, nargs)
        end
    end
    for (op, nargs) in [(-, 1), (ROT, 3), (IROT, 3)]
        printone(op, syms, nargs)
    end
    # abs, conj
end

@i function jf1(op, x)
    op(x[1])
end

@i function jf2(op, x)
    op(x[1], x[2])
end

@i function jf3(op, x)
    op(x[1], x[2], x[3])
end

"""print the jacobian of one operator"""
function printone(op, syms, n)
    if n==1
        jac = jacobian_repeat(jf1, op, syms[1:1]; iin=2, iout=2)
    elseif n==2
        jac = jacobian_repeat(jf2, op, syms[1:2]; iin=2, iout=2)
    elseif n==3
        jac = jacobian_repeat(jf3, op, syms[1:3]; iin=2, iout=2)
    end
    println("------ $op ------")
    pretty_print_matrix(jac)
end

printall()