text:
using Printf

# ------------------------------------
# output results
# ------------------------------------
function output_result(t, x, np, outdir)

    t = string(t)
    while length(t) < 8
        t = "0"*t
    end

    fff = outdir * "/result_" *t* ".dat"
    open(fff,"w") do f
        write(f,"result:x, y, z\n")
        for i in 1:np
            for l in 1:3
                a = @sprintf("%8.8e", x[i,l])
                write(f, a*" ")
            end
            write(f, "\n")
        end
    end
    println("\nwrite "*fff)
end

# ------------------------------
#  make directory
# ------------------------------
function make_dir(outdir)
    k = 0
    try rm(outdir,recursive=true)
    catch
        mkdir(outdir)
        k = 1
    end

    if k == 0
        mkdir(outdir)
    end
end

# ------------------------------
#  flux
# ------------------------------
function set_flux(Qbase, flux, m, Ix, Iy, Iz, Fx, Fy, Fz, Mx, My, Mz)

    flux[1] = Fx / m
    flux[2] = Fy / m
    flux[3] = Fz / m
    flux[4] = Qbase[1]  # v_x
    flux[5] = Qbase[2]  # v_y
    flux[6] = Qbase[3]  # v_z
    flux[7] = ( Iy - Iz ) * Qbase[8] * Qbase[9] / Ix + Mx / Ix # (Iy - Iz) omega_y omega_z / Ix + Mx/Ix
    flux[8] = ( Iz - Ix ) * Qbase[9] * Qbase[7] / Iy + My / Iy # (Iz - Ix) omega_z omega_x / Iy + My/Iy
    flux[9] = ( Ix - Iy ) * Qbase[7] * Qbase[8] / Iz + Mz / Iz # (Ix - Iy) omega_x omega_y / Iz + Mz/Iz
    #flux[10] = 0.5*(Qbase[13]*Qbase[7] - Qbase[12]*Qbase[8] + Qbase[11]*Qbase[9])  # (kai*omega_x - zeta*omega_y + eta*omega_z) / 2
    #flux[11] = 0.5*(Qbase[12]*Qbase[7] + Qbase[13]*Qbase[8] - Qbase[10]*Qbase[9])  # (kai*omega_x - zeta*omega_y + eta*omega_z) / 2
    #flux[12] = 0.5*(-Qbase[11]*Qbase[7] + Qbase[10]*Qbase[8] + Qbase[13]*Qbase[9])  # (kai*omega_x - zeta*omega_y + eta*omega_z) / 2
    #flux[13] = 0.5*(-Qbase[10]*Qbase[7] - Qbase[11]*Qbase[8] - Qbase[12]*Qbase[9])  # (kai*omega_x - zeta*omega_y + eta*omega_z) / 2

    #q0 = 0.5*( Qbase[13]*Qbase[7] - Qbase[12]*Qbase[8] + Qbase[11]*Qbase[9])  # (kai*omega_x - zeta*omega_y + eta*omega_z) / 2
    #q1 = 0.5*( Qbase[12]*Qbase[7] + Qbase[13]*Qbase[8] - Qbase[10]*Qbase[9])  # (kai*omega_x - zeta*omega_y + eta*omega_z) / 2
    #q2 = 0.5*(-Qbase[11]*Qbase[7] + Qbase[10]*Qbase[8] + Qbase[13]*Qbase[9])  # (kai*omega_x - zeta*omega_y + eta*omega_z) / 2
    #q3 = 0.5*(-Qbase[10]*Qbase[7] - Qbase[11]*Qbase[8] - Qbase[12]*Qbase[9])  # (kai*omega_x - zeta*omega_y + eta*omega_z) / 2

    q0 = 0.5*( Qbase[13]*Qbase[7] + Qbase[12]*Qbase[8] - Qbase[11]*Qbase[9])  # (kai*omega_x - zeta*omega_y + eta*omega_z) / 2
    q1 = 0.5*(-Qbase[12]*Qbase[7] + Qbase[13]*Qbase[8] + Qbase[10]*Qbase[9])  # (kai*omega_x - zeta*omega_y + eta*omega_z) / 2
    q2 = 0.5*( Qbase[11]*Qbase[7] - Qbase[10]*Qbase[8] + Qbase[13]*Qbase[9])  # (kai*omega_x - zeta*omega_y + eta*omega_z) / 2
    q3 = 0.5*(-Qbase[10]*Qbase[7] - Qbase[11]*Qbase[8] - Qbase[12]*Qbase[9])  # (kai*omega_x - zeta*omega_y + eta*omega_z) / 2

    n = ( q0^2 + q1^2 + q2^2 + q3^2 )^0.5

    flux[10] = q0
    flux[11] = q1
    flux[12] = q2
    flux[13] = q3
    #=
    if n == 0.0
        flux[10] = 0.0
        flux[11] = 0.0
        flux[12] = 0.0
        flux[13] = 0.0
    else
        flux[10] = q0 / n
        flux[11] = q1 / n
        flux[12] = q2 / n
        flux[13] = q3 / n
    end
    =#
    return flux
end

# ------------------------------
#  rotation
# ------------------------------
function rotation(x,y,z, Qbase)
    A = zeros(3, 3)
    q0 = Qbase[10]
    q1 = Qbase[11]
    q2 = Qbase[12]
    q3 = Qbase[13]

    #=

    A[1,1] = q0^2 + q1^2 - q2^2 - q3^2
    A[1,2] = 2*(q1*q2 + q0*q3)
    A[1,3] = 2*(q1*q3 - q0*q2)
    A[2,1] = 2*(q1*q2 - q0*q3)
    A[2,2] = q0^2 - q1^2 + q2^2 - q3^2
    A[2,3] = 2*(q2*q3 + q0*q1)
    A[3,1] = 2*(q1*q3 + q0*q2)
    A[3,2] = 2*(q2*q3 - q0*q1)
    A[3,3] = q0^2 - q1^2 - q2^2 + q3^2
    =#

    q1 = Qbase[13]
    q2 = Qbase[10]
    q3 = Qbase[11]
    q4 = Qbase[12]

    A[1,1] = 1 - 2*q3^2 - 2*q4^2
    A[1,2] = 2*q2*q3 + 2*q1*q4
    A[1,3] = 2*q2*q4 - 2*q1*q3

    A[2,1] = 2*q2*q3 - 2*q1*q4
    A[2,2] = 1 - 2*q4^2 - 2*q2^2
    A[2,3] = 2*q3*q4 + 2*q1*q2

    A[3,1] = 2*q2*q4 + 2*q1*q3
    A[3,2] = 2*q3*q4 - 2*q1*q2
    A[3,3] = 1 - 2*q2^2 - 2*q3^2

    newx = A[1,1]*x + A[1,2]*y + A[1,3]*z
    newy = A[2,1]*x + A[2,2]*y + A[2,3]*z
    newz = A[3,1]*x + A[3,2]*y + A[3,3]*z
    return newx, newy, newz
end

# ------------------------------
#  main
# ------------------------------
function main()

    m = 1.0 # [kg]

    Fx = 0.0 # [N]
    Fy = 0.0 # [N]
    Fz = 0.0 # [N]
    Mx = 0.0 # [N/m]
    My = 0.0 # [N/m]
    Mz = 0.0e0 # [N/m]

    init_v_x = 0.0
    init_v_y = 0.0
    init_v_z = 0.0
    init_gc_x = 0.0
    init_gc_y = 0.0
    init_gc_z = 0.0
    init_omega_x = 1.5e-4
    init_omega_y = 0.0
    init_omega_z = 0.0

    dt = 1.0e-6 # [s]
    nt = 1.0e6
    outstep = 1.0e4

    outdir = "post_result"
    make_dir(outdir)
    outdir = "result"
    make_dir(outdir)

    # quontanion
    # q = (1, 0, 0, 0)
    # q = (ux sin(theta/2), uy sin(theta/2), uz sin(theta/2), cos(theta/2))
    # これにより初期角度と軸を決定できる
    # u = (ux, uy, uz) ：回転軸方向の単位ベクトル
    init_q_0 = 0.0
    init_q_1 = 0.0
    init_q_2 = 0.0
    init_q_3 = 1.0

    # ------------------------------
    #  make block
    # ------------------------------
    np = 8
    x = zeros(np, 3)

    for i in 1:np
        if (i-1) % 8 < 4
            x[i, 1] = 1.0
        else
            x[i, 1] = -1.0
        end

        if (i-1) % 4 < 2
            x[i, 2] = 1.0
        else
            x[i, 2] = -1.0
        end

        if (i-1) % 2 < 1
            x[i, 3] = 1.0
        else
            x[i, 3] = -1.0
        end
    end

    # ------------------------------
    #  primitive
    # ------------------------------
    # vecVelo, GC, vecOmega, くっぉーたにオン
    nprim = 13
    Qbase = zeros(nprim)
    Qbase[1] = init_v_x
    Qbase[2] = init_v_y
    Qbase[3] = init_v_z
    Qbase[4] = init_gc_x
    Qbase[5] = init_gc_y
    Qbase[6] = init_gc_z
    Qbase[7] = init_omega_x
    Qbase[8] = init_omega_y
    Qbase[9] = init_omega_z
    Qbase[10] = init_q_0
    Qbase[11] = init_q_1
    Qbase[12] = init_q_2
    Qbase[13] = init_q_3

    flux = zeros(nprim)
    # ------------------------------
    # innertia moment
    # ------------------------------
    # 長さ2より
    Ix = m * (2^2 + 2^2) / 12
    Iy = m * (2^2 + 2^2) / 12
    Iz = m * (2^2 + 2^2) / 12

    # ------------------------------
    # main loop
    # ------------------------------
    for t in 1:Int(nt)

        flux = set_flux(Qbase, flux, m, Ix, Iy, Iz, Fx, Fy, Fz, Mx, My, Mz)

        # Euler explicit
        for i in 1:nprim
            Qbase[i] = Qbase[i] + dt*flux[i]
        end

        # 規格化
        #=
        q0 = Qbase[10]
        q1 = Qbase[11]
        q2 = Qbase[12]
        q3 = Qbase[13]
        n = ( q0^2 + q1^2 + q2^2 + q3^2 )^0.5
        Qbase[10] = q0*n
        Qbase[11] = q1*n
        Qbase[12] = q2*n
        Qbase[13] = q3*n
        =#

        #println(Qbase)

        # move all point
        for i in 1:np
            x[i,1] += Qbase[1] * dt
            x[i,2] += Qbase[2] * dt
            x[i,3] += Qbase[3] * dt

            # rotation
            x[i,1], x[i,2], x[i,3] = rotation(x[i,1], x[i,2], x[i,3], Qbase)
        end

        #=
        println("throw(UndefVarError(:x))")
        println(Qbase)
        println(flux)
        println(x)
        =#

        if t % Int(outstep) == 0
            output_result(t, x, np, outdir)
        end

        #=

        if t == 100
            Qbase[7] = -1.5e-3
        elseif t == 200
            Qbase[7] = 1.5e-3
        else
            Qbase[7] = 0.0
        end
        =#

        #=

        q0 = Qbase[10]
        q1 = Qbase[11]
        q2 = Qbase[12]
        q3 = Qbase[13]

        theta = acos(q3)*2
        ux    = q0 / sin(theta/2)
        uy    = q1 / sin(theta/2)
        uz    = q2 / sin(theta/2)
        println(" ")
        println(theta)
        println(ux)
        println(uy)
        println(uz)
        #throw(UndefVarError(:x))
        =#

        if isequal(x[1,1], NaN) == true
            println("\n -------------------- \n ")
            println(" ")
            println(" diverge ")
            println(" ")
            println("\n -------------------- \n ")
            throw(UndefVarError(:x))
        end

        if abs(x[1,1]) > 2
            println("\n -------------------- \n ")
            println(" ")
            println(" use small dt ")
            println(" ")
            println("\n -------------------- \n ")
            throw(UndefVarError(:x))
        end

    end

end
# ------------------------------
# ------------------------------
main()


text:
"""
    chr2ind(s, i, offset)

Find a byte index for character in string `s`.
We start at **byte index** `i` with an `offset` number of characters.
"""
function chr2ind(s::AbstractString, i::Integer, offset::Integer)
    k = i
    if offset < 0
        for j in 1:abs(offset)
            k = prevind(s, k)
        end
    elseif offset > 0
        for j in 1:offset
            k = nextind(s, k)
        end
    end
    clamp(k, 1, lastindex(s))
end

ext: jl


text:
# Hidden Markov models

# Copyright (c) 2016 Jeffrey W. Miller
# This software is released under the MIT License.
# If you use this software in your research, please cite:
#
#   Jeffrey W. Miller (2016). Lecture Notes on Advanced Stochastic Modeling. Duke University, Durham, NC.

# HMM with discrete emission distributions, assumed to be on {1,...,K}.
module DiscreteHMM
	using Distributions

	# Emission distribution specification
	log_q(x,p) = log(p[x])  # log density
	qrnd(p) = rand(Categorical(p))  # random sample

	# Initialize an array of emission distribution parameters
	function phi_init(m,x)
		K = maximum(x)
	    return [rand(Dirichlet(K,1)) for s = 1:m]
	end

	# Maximum likelihood estimate of array of emission distribution
	# parameters using data x with weights gamma.
	function phi_max!(phi,x,gamma)
	    m = length(phi)
	    n = length(x)
	    for s = 1:m
	    	ph = phi[s]
	    	ph[:] = 0.0
	        for j = 1:n
	            ph[x[j]] += gamma[j,s]
	        end
	        ph[:] /= sum(ph)
	    end
	end

	# Core HMM code
	include("hmm.jl")
end # module

# HMM with univariate normal emission distributions.
module NormalHMM
	using Distributions

	# Emission distribution specification
	log_q(x,D) = logpdf(D,x)  # log density
	qrnd(D) = rand(D)  # random sample

	# Initialize an array of emission distribution parameters
	function phi_init(m,x)
		D_m = Normal(mean(x),std(x))
		D_v = InverseGamma(1.0,var(x))
	    return [Normal(rand(D_m),sqrt(rand(D_v))) for s = 1:m]
	end

	# Maximum likelihood estimate of array of emission distribution
	# parameters using data x with weights gamma.
	function phi_max!(phi,x,gamma)
	    m = length(phi)
	    n = length(x)
	    for s = 1:m
	         weights = gamma[:,s] / sum(gamma[:,s])
             mu = dot(weights,x)
             sigma = sqrt(dot(weights,(x-mu).*(x-mu)))
             phi[s] = Normal(mu,sigma)
	    end
	end

	# Core HMM code
	include("hmm.jl")
end # module

# HMM with multivariate normal emission distributions,
# with diagonal covariance matrices.
module MVNdiagHMM
	using Distributions

	# Emission distribution specification
	log_q(x,D) = logpdf(D,x)  # log density
	qrnd(D) = rand(D)  # random sample

	# Initialize an array of emission distribution parameters
	function phi_init(m,x)
		d = length(x[1])  # dimension
		mu = mean(x)
		sigma = Float64[std(Float64[xi[j] for xi in x]) for j=1:d]
		randsig(d) = sqrt(rand(InverseGamma(1,1),d))
		phi = [MvNormal(randn(d).*sigma+mu, randsig(d).*sigma) for s=1:m]
	end

	# Maximum likelihood estimate of array of emission distribution
	# parameters using data x with weights gamma.
	function phi_max!(phi,x,gamma)
	    m = length(phi)
	    n = length(x)
		d = length(x[1])  # dimension
	    for s = 1:m
	         weights = gamma[:,s] / sum(gamma[:,s])
	         mu = zeros(d)
	         for j = 1:n
	         	mu += weights[j]*x[j]
	         end
	         sigma2 = zeros(d)
	         for j = 1:n
	         	sigma2 += weights[j]*(x[j]-mu).*(x[j]-mu)
	         end
             phi[s] = MvNormal(mu,sqrt(sigma2))
	    end
	end

	# Core HMM code
	include("hmm.jl")
end # module

max_stars_repo_name: jwmi/HMM

max_issues_repo_name: jwmi/HMM

max_forks_repo_name: jwmi/HMM


text:
facts("HTML Table") do
  for tblname in tblnames
    println(tblname)
    pb = PipeBuffer()
    rendercells(pb, tocellarray(eval(tblname)))
    @fact takebuf_string(pb) == readall("./html/examples/$(tblname).txt") => true
  end
end


text:
module TestNodeCommand

using ComoniconTypes
using ComoniconTargetExpr
using ComoniconTargetExpr: emit_expr, emit_body, emit_norm_body, emit_dash_body
using Test

function called()
    @test true
end

cmd = Entry(;
    version=v"1.2.0",
    root=NodeCommand(;
        name="node",
        subcmds=Dict(
            "cmd1" => LeafCommand(;
                fn=called,
                name="cmd1",
            ),
            "cmd2" => LeafCommand(;
                fn=called,
                name="cmd2",
            )
        )
    )
)

eval(emit_expr(cmd))

@testset "test node" begin
    @test command_main(["cmd3"]) == 1
    @test command_main(["cmd1", "foo"]) == 1
    @test command_main(["cmd1", "foo", "-h"]) == 0
    @test command_main(["cmd1", "foo", "-V"]) == 0
    @test command_main(String[]) == 1
end

end


text:
using GridInterpolations
using LocalFunctionApproximation
using LocalApproximationValueIteration
using POMDPs
using POMDPModels
using POMDPModelTools
using StaticArrays
using JLD2, FileIO
using Random
using Logging
using Distributions
using CMSSPs

# Read parameters from file
scale_file = ARGS[1]
simtime_file = ARGS[2]
cost_file = ARGS[3]
policy_name = ARGS[4]
poly_or_exp = ARGS[5]
beta = parse(Float64, ARGS[6])

# scale_file = "./paramsets/scale-1.toml"
# simtime_file = "./paramsets/simtime-1.toml"
# cost_file = "./paramsets/cost-1.toml"
# policy_name = "dreamr-cf-params1-betapt75"
# poly_or_exp = "poly"

rng = MersenneTwister(10)

params = parse_params(scale_file=scale_file, simtime_file=simtime_file, cost_file=cost_file)

cf_mdp = get_cf_mdp(MultiRotorUAVState, MultiRotorUAVAction, params, beta)

if poly_or_exp == "poly"

    xy_spacing = polyspace_symmetric(params.scale_params.XY_LIM, params.scale_params.XY_AXISVALS)
    xydot_spacing = polyspace_symmetric(params.scale_params.XYDOT_LIM, params.scale_params.XYDOT_AXISVALS)

elseif poly_or_exp == "exp"

    xy_spacing = log2space_symmetric(params.scale_params.XY_LIM, params.scale_params.XY_AXISVALS)
    xydot_spacing = log2space_symmetric(params.scale_params.XYDOT_LIM, params.scale_params.XYDOT_AXISVALS)

end

state_grid = RectangleGrid(xy_spacing, xy_spacing, xydot_spacing, xydot_spacing)
@show length(vertices(state_grid))
lfa = LocalGIFunctionApproximator(state_grid)

compute_terminalcost_localapprox!(cf_mdp, lfa)

cf_policy = finite_horizon_VI_localapprox!(cf_mdp, lfa, true, params.scale_params.MC_GENERATIVE_NUMSAMPLES, rng)

compute_min_value_per_horizon_localapprox!(cf_policy)

inhor_fn = string(policy_name, "-inhor.jld2")
outhor_fn = string(policy_name, "-outhor.jld2")
save_modal_horizon_policy_localapprox(cf_policy, inhor_fn,outhor_fn, cf_mdp)


text:
using Pkg
pkg"precompile"
using DoubleFloats

using MPI
using Test

import MPI: mpiexec_path

# Code coverage command line options; must correspond to src/julia.h
# and src/ui/repl.c
const JL_LOG_NONE = 0
const JL_LOG_USER = 1
const JL_LOG_ALL = 2
const coverage_opts =
    Dict{Int, String}(JL_LOG_NONE => "none",
                      JL_LOG_USER => "user",
                      JL_LOG_ALL => "all")

# Files to run with mpiexec -n 1
singlefiles = ["test_spawn.jl"]

function runtests()
    nprocs = clamp(Sys.CPU_THREADS, 2, 4)
    exename = joinpath(Sys.BINDIR, Base.julia_exename())
    testdir = dirname(@__FILE__)
    istest(f) = endswith(f, ".jl") && startswith(f, "test_")
    testfiles = sort(filter(istest, readdir(testdir)))
    mpiexec_args = split(get(ENV, "JULIA_MPIEXEC_ARGS", ""))

    if !Sys.iswindows() && occursin( "OpenRTE", read(`$mpiexec_path --version`, String))
        push!(mpiexec_args,"--oversubscribe")
    end

    nfail = 0
    printstyled("Running MPI.jl tests\n"; color=:white)

    for f in testfiles
        coverage_opt = coverage_opts[Base.JLOptions().code_coverage]
        if f ∈ singlefiles
            run(`$mpiexec_path $mpiexec_args -n 1 $exename --code-coverage=$coverage_opt $(joinpath(testdir, f))`)
        else
            run(`$mpiexec_path $mpiexec_args -n $nprocs $exename --code-coverage=$coverage_opt $(joinpath(testdir, f))`)
        end
        Base.with_output_color(:green,stdout) do io
            println(io,"\tSUCCESS: $f")
        end
    end
    return nfail
end

exit(runtests())


text:
module InfiniteArrays
using Base, Statistics, LinearAlgebra, FillArrays, Infinities, LazyArrays, DSP, ArrayLayouts

import Base: *, +, -, /, \, ==, isinf, isfinite, sign, signbit, angle, show, isless,
            fld, cld, div, min, max, minimum, maximum, mod,
            <, ≤, >, ≥, promote_rule, convert, unsafe_convert, copy,
            size, step, isempty, length, first, last, tail,
            getindex, setindex!, intersect, @_inline_meta,
            sort, sort!, issorted, sortperm, sum, in, broadcast,
            eltype, elsize, parent, parentindices, reinterpret,
            unaliascopy, dataids,
            real, imag, conj, transpose,
            exp, log, sqrt, cos, sin, tan, csc, sec, cot,
            cosh, sinh, tanh, csch, sech, coth, acos, asin, atan, acsc, asec, acot,
            acosh, asinh, atanh, acsch, asech, acoth, (:),
            AbstractMatrix, AbstractArray, checkindex, unsafe_length, unsafe_indices, OneTo,
            to_shape, _sub2ind, print_matrix, print_matrix_row, print_matrix_vdots,
            checkindex, Slice, IdentityUnitRange, @propagate_inbounds, @_propagate_inbounds_meta,
         	_in_range, _range, _rangestyle, Ordered,
         	ArithmeticWraps, ArithmeticUnknown, floatrange, reverse, unitrange_last,
         	AbstractArray, AbstractVector, Array, Vector, Matrix,
         	axes, (:), _sub2ind_recurse, broadcast, promote_eltypeof,
         	diff, cumsum, show_delim_array, show_circular, Int,
         	similar, _unsafe_getindex, string, zeros, ones, fill, permutedims,
         	cat_similar, vcat, one, zero,
		 	reshape, ReshapedIndex, ind2sub_rs, _unsafe_getindex_rs,
            searchsorted, searchsortedfirst, searchsortedlast, Ordering, lt, Fix2, findfirst,
            cat_indices, cat_size, cat_similar, __cat, _ind2sub_recurse, union, intersect

using Base.Broadcast
import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcasted,
                        @nexprs, @ncall, combine_eltypes, DefaultArrayStyle, instantiate, axistype

import LinearAlgebra: BlasInt, BlasFloat, norm, diag, diagm, ishermitian, issymmetric,
                             det, logdet, istriu, istril, adjoint, tr, AbstractTriangular,
                             norm2, norm1, normp, AdjOrTrans, HermOrSym

import Statistics: mean, median

import FillArrays: AbstractFill, getindex_value, fill_reshape, RectDiagonal, Fill, Ones, Zeros, Eye
import LazyArrays: LazyArrayStyle, AbstractBandedLayout, MemoryLayout, LazyLayout, UnknownLayout,
                    ZerosLayout, AbstractArrayApplyStyle, AbstractCachedVector, CachedArray, CachedVector, ApplyLayout, LazyMatrix,
                    reshapedlayout, sub_materialize, LayoutMatrix, LayoutVector, _padded_sub_materialize, PaddedLayout

import DSP: conv
import ArrayLayouts: RangeCumsum
import Infinities: ∞, Infinity, InfiniteCardinal

export ∞, ℵ₀, Hcat, Vcat, Zeros, Ones, Fill, Eye, BroadcastArray, cache

if VERSION ≥ v"1.6-"
   import Base: unitrange, oneto
end

include("infrange.jl")
include("infarrays.jl")
include("reshapedarray.jl")

if VERSION < v"1.6-"
   ##
   # Temporary hacks for base support
   ##
   Base.OneTo(::Infinity) = OneToInf()
   Base.OneTo(::InfiniteCardinal{0}) = OneToInf()
   function Base.OneTo(x::ComplexInfinity)
      iszero(x.angle) && return oneto(∞)
      throw(ArgumentError("Cannot create infinite range with negative length"))
   end
   function Base.OneTo(x::RealInfinity)
      signbit(x) || return oneto(∞)
      throw(ArgumentError("Cannot create infinite range with negative length"))
   end
   Base.OneTo{T}(::Infinity) where T<:Integer = OneToInf{T}()
   Base.UnitRange(start::Integer, ::Infinity) = InfUnitRange(start)
   Base.UnitRange{T}(start::Integer, ::Infinity) where T<:Real = InfUnitRange{T}(start)
   Base.OneTo{T}(::InfiniteCardinal{0}) where T<:Integer = OneToInf{T}()
   Base.UnitRange(start::Integer, ::InfiniteCardinal{0}) = InfUnitRange(start)
   Base.UnitRange{T}(start::Integer, ::InfiniteCardinal{0}) where T<:Real = InfUnitRange{T}(start)
   Base.OneTo(a::OneToInf) = a
   Base.OneTo{T}(::OneToInf) where T<:Integer = OneToInf{T}()

   Base.Int(::Infinity) = ℵ₀
   Base.Int(::InfiniteCardinal{0}) = ℵ₀

   unitrange(a, b) = UnitRange(a, b)
   oneto(n) = Base.OneTo(n)
end

##
# Fill FillArrays
##

length(::Ones{<:Any,1,Tuple{OneToInf{Int}}}) = ℵ₀
length(::Fill{<:Any,1,Tuple{OneToInf{Int}}}) = ℵ₀
length(::Zeros{<:Any,1,Tuple{OneToInf{Int}}}) = ℵ₀
length(::Ones{<:Any,2,Tuple{OneToInf{Int},OneToInf{Int}}}) = ℵ₀
length(::Ones{<:Any,2,<:Tuple{OneToInf{Int},<:Any}}) = ℵ₀
length(::Ones{<:Any,2,<:Tuple{<:Any,OneToInf{Int}}}) = ℵ₀
length(::Fill{<:Any,2,Tuple{OneToInf{Int},OneToInf{Int}}}) = ℵ₀
length(::Fill{<:Any,2,<:Tuple{OneToInf{Int},<:Any}}) = ℵ₀
length(::Fill{<:Any,2,<:Tuple{<:Any,OneToInf{Int}}}) = ℵ₀
length(::Zeros{<:Any,2,Tuple{OneToInf{Int},OneToInf{Int}}}) = ℵ₀
length(::Zeros{<:Any,2,<:Tuple{OneToInf{Int},<:Any}}) = ℵ₀
length(::Zeros{<:Any,2,<:Tuple{<:Any,OneToInf{Int}}}) = ℵ₀

for op in (:norm2, :norm1)
   @eval $op(a::Zeros{T,N,NTuple{N,OneToInf{Int}}}) where {T,N} = norm(getindex_value(a))
end

normp(a::Zeros{T,N,NTuple{N,OneToInf{Int}}}, p) where {T,N} = norm(getindex_value(a))

for N=1:3
   for op in (:norm2, :norm1)
      @eval function $op(a::AbstractFill{T,$N,NTuple{$N,OneToInf{Int}}}) where {T,N}
         z = norm(getindex_value(a))
         iszero(z) && return z
         typeof(z)(Inf)
      end
   end
   @eval function normp(a::AbstractFill{T,$N,NTuple{$N,OneToInf{Int}}}, p) where {T,N }
      z = norm(getindex_value(a))
      iszero(z) && return z
      typeof(z)(Inf)
   end
end

for Typ in (:Number, :AbstractVector)
   @eval begin
      vcat(a::$Typ, b::AbstractFill{<:Any,1,<:Tuple{OneToInf}}) = Vcat(a, b)
      vcat(a::$Typ, c::CachedVector{<:Any,<:Any,<:AbstractFill{<:Any,1,<:Tuple{OneToInf}}}) =
         CachedArray(vcat(a, view(c.data,1:c.datasize[1])), c.array)
   end
end

vcat(a::AbstractVector, b::AbstractVector, c::AbstractFill{<:Any,1,<:Tuple{OneToInf}}) = Vcat(vcat(a,b), c)
vcat(a::AbstractVector, b::AbstractVector, c::AbstractVector, d::AbstractFill{<:Any,1,<:Tuple{OneToInf}}) = Vcat(vcat(a,b,c), d)

vcat(a::AbstractMatrix, b::AbstractFill{<:Any,2,<:Tuple{OneToInf,OneTo}}) = Vcat(a, b)

cat_similar(A, ::Type{T}, shape::Tuple{PosInfinity}) where T = zeros(T,∞)
if VERSION < v"1.6-"
   cat_similar(A::AbstractArray, ::Type{T}, shape::Tuple{PosInfinity}) where T =
      Base.invoke(cat_similar, Tuple{AbstractArray, Any, Any}, A, T, shape)
else
   cat_similar(A::AbstractArray, ::Type{T}, shape::Tuple{PosInfinity}) where T =
      Base.invoke(cat_similar, Tuple{AbstractArray, Type{T}, Any}, A, T, shape)
end
function Base.__cat(A, shape::NTuple{N,PosInfinity}, catdims, X...) where N
   offsets = zeros(Union{Int,InfiniteCardinal{0}}, N)
   inds = Vector{Union{UnitRange{Int},InfUnitRange{Int}}}(undef, N)
   concat = copyto!(zeros(Bool, N), catdims)
   for x in X
       for i = 1:N
           if concat[i]
               inds[i] = offsets[i] .+ cat_indices(x, i)
               offsets[i] += cat_size(x, i)
           else
               inds[i] = 1:shape[i]
           end
       end
       I::NTuple{N, Union{InfUnitRange{Int},UnitRange{Int}}} = (inds...,)
       if x isa AbstractArray
           copyto!(view(A, I...), x)
       else
           fill!(view(A, I...), x)
       end
   end
   return A
end

reshape(parent::AbstractArray, shp::Tuple{OneToInf, Vararg{Union{Integer,OneTo,OneToInf}}}) =
   reshape(parent, to_shape(shp))
reshape(parent::AbstractArray, shp::Tuple{Union{Integer,OneTo}, OneToInf, Vararg{Union{Integer,OneTo,OneToInf}}}) =
   reshape(parent, to_shape(shp))

# cat_similar(A, T, ::Tuple{PosInfinity}) = zeros(T, ∞)

axistype(::OneTo{T}, ::OneToInf{V}) where {T,V} = OneToInf{promote_type(T,V)}()
axistype(::OneToInf{V}, ::OneTo{T}) where {T,V} = OneToInf{promote_type(T,V)}()

# sort.jl
# returns the range of indices of v equal to x
# if v does not contain x, returns a 0-length range
# indicating the insertion point of x
function searchsorted(v::AbstractVector, x, ilo::Int, ::PosInfinity, o::Ordering)
    lo = ilo-1
    hi = ℵ₀
    @inbounds while lo < hi-1
        m = isinf(hi) ? lo + 1000 : (lo+hi)>>>1
        if lt(o, v[m], x)
            lo = m
        elseif lt(o, x, v[m])
            hi = m
        else
            a = searchsortedfirst(v, x, max(lo,ilo), m, o)
            b = searchsortedlast(v, x, m, hi, o)
            return a : b
        end
    end
    return (lo + 1) : (hi - 1)
end

# index of the first value of vector a that is greater than or equal to x;
# returns length(v)+1 if x is greater than all values in v.
function searchsortedfirst(v::AbstractVector, x, lo::Int, hi::PosInfinity, o::Ordering)
   u = 1
   lo = lo - u
   hi = ℵ₀
   @inbounds while lo < hi - u
      m = isinf(hi) ? lo + 1000 : (lo+hi)>>>1
      if lt(o, v[m], x)
         lo = m
      else
         hi = m
      end
   end
   return hi
end

# index of the last value of vector a that is less than or equal to x;
# returns 0 if x is less than all values of v.
function searchsortedlast(v::AbstractVector, x, lo::Int, hi::PosInfinity, o::Ordering)
   u = 1
   lo = lo - u
   hi = ℵ₀
   @inbounds while lo < hi - u
       m = isinf(hi) ? lo + 1000 : (lo+hi)>>>1
       if lt(o, x, v[m])
           hi = m
       else
           lo = m
       end
   end
   return lo
end

# special case for Vcat
@inline function LazyArrays.searchsortedlast_recursive(::PosInfinity, x, a, args...)
    n = sum(map(length,args))
    r = searchsortedlast(a, x)
    r > 0 && return n + r
    return LazyArrays.searchsortedlast_recursive(n, x, args...)
end

##
# lazy sub_materialize
##

const InfAxes = Union{AbstractInfUnitRange,Slice{<:AbstractInfUnitRange},IdentityUnitRange{<:AbstractInfUnitRange}}

sub_materialize(_, V, ::Tuple{InfAxes}) = V
sub_materialize(_, V, ::Tuple{InfAxes,InfAxes}) = V
sub_materialize(_, V, ::Tuple{<:Any,InfAxes}) = V
sub_materialize(_, V, ::Tuple{InfAxes,Any}) = V

sub_materialize(::ApplyLayout{typeof(vcat)}, V::AbstractVector, ::Tuple{InfAxes}) = ApplyArray(V)
sub_materialize(::ApplyLayout{typeof(vcat)}, V::AbstractMatrix, ::Tuple{InfAxes,InfAxes}) = ApplyArray(V)
sub_materialize(::ApplyLayout{typeof(vcat)}, V::AbstractMatrix, ::Tuple{<:Any,InfAxes}) = ApplyArray(V)
sub_materialize(::ApplyLayout{typeof(vcat)}, V::AbstractMatrix, ::Tuple{InfAxes,Any}) = ApplyArray(V)

sub_materialize(::PaddedLayout, v::AbstractVector{T}, ::Tuple{InfAxes}) where T =
    _padded_sub_materialize(v)

end # module


text:
@deprecate VarName(sym::Symbol) VarName{sym}()
@deprecate VarName(sym::Symbol, indexing::Tuple) VarName{sym}(indexing)


text:
struct JLSize{name, T<:Tuple{Symbol, Vararg{Symbol}}}
    labels::T
end
JLSize{name}(labels...) where {name} = JLSize{name, typeof(labels)}(labels)
Base.nameof(::JLSize{name}) where name = name
Base.length(x::JLSize) = length(x.labels)

"""
    FieldInfo

Metadata for fields of the various MuJoCo structs.
Sizes are declared in the same order as C/row-major order
and should be reversed in Julia. Trailing/leading singleton
dimensions are dropped.

# Fields
- mjsize: the MuJoCo C size declared
- dsize: the dynamic size for pre-allocated, variable length fields (e.g. `mjData.contact`
    is declared as `(nconmax, )` but should be viewed as `(ncon,)`.
- jlsize: The size for "jl" prefixed structs e.g. `jlData`. By default this is
    `mjsize`, but can be provide new sizes
    (e.g. `mjData.xmat` is `(9, nbody)`` while `jlData.xmat` is `(:xmat, nbody)`
    which allows for named indexing like `jlData.xmat[:xz, :foobody]`
"""
struct FieldInfo
    mjsize::Tuple{Vararg{Union{Int,Symbol}}}
    jlsize::Union{Tuple{Vararg{Union{Int,Symbol, JLSize}}}}
    dsize::Union{Tuple{Vararg{Union{Int,Symbol, JLSize}}},Nothing}
end
FieldInfo(mjsize; jlsize=mjsize, dsize=nothing) = FieldInfo(mjsize, jlsize, dsize)

const Quat = JLSize{:quat}(:w, :x, :y, :z)
const XYZ = JLSize{:xyz}(:x, :y, :z)
const RGBA = JLSize{:rgba}(:r, :b, :g, :a)
const Range = JLSize{:range}(:min, :max)
const FrameAxis = JLSize{:frameaxis}(:x, :y, :z)
const JLSizes = (Quat, XYZ, RGBA, Range, FrameAxis)

frame() = (FrameAxis, XYZ)

# source is the thing which holds the runtime size e.g. mjModel
function lookupsize(structname::Symbol, source, fieldname::Symbol)
    jl_size = reverse(STRUCT_INFO[structname][fieldname].mjsize) # reverse b/c Julia is column order
    Base.map(jl_size) do sz
        sz = sz isa Symbol ? getfield(source, sz) : sz
        Int(sz)
    end
end

function lookupsize(structname::Symbol, fieldname::Symbol)
    jl_size = reverse(STRUCT_INFO[structname][fieldname].mjsize) # reverse b/c Julia is column order
    Base.map(jl_size) do sz
        sz isa Symbol && error("Found Symbol, try lookupsize(structname, source, fieldname)")
        Int(sz)
    end
end

@generated function lookup_dsize(
    ::Val{structname},
    source,
    ::Val{fieldname},
) where {structname,fieldname}
    jl_size = reverse(STRUCT_INFO[structname][fieldname].dsize) # reverse b/c Julia is column order
    args = Base.map(jl_size) do sz
        sz = sz isa Symbol ? :(Int(getfield(source, $(QuoteNode(sz))))) : sz
    end
    N = length(args)
    :($(Expr(:tuple, args...))::NTuple{$N,Int})
end

const STRUCT_INFO = Dict{Symbol,Dict{Symbol,<:FieldInfo}}(
    :mjModel => Dict(
        :tex_rgb => FieldInfo((:ntexdata,)),
        :jnt_bodyid => FieldInfo((:njnt,)),
        :tendon_group => FieldInfo((:ntendon,)),
        :mesh_faceadr => FieldInfo((:nmesh,)),
        :name_excludeadr => FieldInfo((:nexclude,)),
        :tuple_objid => FieldInfo((:ntupledata,)),
        :actuator_forcerange => FieldInfo((:nu, 2), jlsize=(:nu, Range)),
        :sensor_noise => FieldInfo((:nsensor,)),
        :tendon_range => FieldInfo((:ntendon, 2), jlsize=(:ntendon, Range)),
        :body_ipos => FieldInfo((:nbody, 3), jlsize=(:nbody, XYZ)),
        :body_dofadr => FieldInfo((:nbody,)),
        :geom_solmix => FieldInfo((:ngeom,)),
        :skin_vert => FieldInfo((:nskinvert, 3), jlsize=(:nskinvert, XYZ)),
        :wrap_type => FieldInfo((:nwrap,)),
        :name_lightadr => FieldInfo((:nlight,)),
        :skin_bonenum => FieldInfo((:nskin,)),
        :sensor_type => FieldInfo((:nsensor,)),
        :tuple_objtype => FieldInfo((:ntupledata,)),
        :hfield_nrow => FieldInfo((:nhfield,)),
        :body_dofnum => FieldInfo((:nbody,)),
        :key_qpos => FieldInfo((:nkey, :nq)),
        :names => FieldInfo((:nnames,)),
        :tendon_margin => FieldInfo((:ntendon,)),
        :skin_bonebodyid => FieldInfo((:nskinbone,)),
        :light_ambient => FieldInfo((:nlight, 3)),
        :tuple_adr => FieldInfo((:ntuple,)),
        :geom_priority => FieldInfo((:ngeom,)),
        :dof_solimp => FieldInfo((:nv, 5)),
        :geom_size => FieldInfo((:ngeom, 3)),
        :mesh_texcoordadr => FieldInfo((:nmesh,)),
        :eq_active => FieldInfo((:neq,)),
        :actuator_forcelimited => FieldInfo((:nu,)),
        :body_mass => FieldInfo((:nbody,)),
        :skin_bonebindpos => FieldInfo((:nskinbone, 3), jlsize=(:nskinbone, XYZ)),
        :hfield_ncol => FieldInfo((:nhfield,)),
        :cam_fovy => FieldInfo((:ncam,)),
        :name_geomadr => FieldInfo((:ngeom,)),
        :light_castshadow => FieldInfo((:nlight,)),
        :tendon_invweight0 => FieldInfo((:ntendon,)),
        :jnt_limited => FieldInfo((:njnt,)),
        :eq_type => FieldInfo((:neq,)),
        :body_mocapid => FieldInfo((:nbody,)),
        :body_inertia => FieldInfo((:nbody, 3)),
        :skin_faceadr => FieldInfo((:nskin,)),
        :name_textadr => FieldInfo((:ntext,)),
        :skin_texcoord => FieldInfo((:nskintexvert, 2)),
        :skin_face => FieldInfo((:nskinface, 3)),
        :mesh_facenum => FieldInfo((:nmesh,)),
        :sensor_objid => FieldInfo((:nsensor,)),
        :body_invweight0 => FieldInfo((:nbody, 2)),
        :cam_bodyid => FieldInfo((:ncam,)),
        :light_pos => FieldInfo((:nlight, 3), jlsize=(:nlight, XYZ)),
        :light_attenuation => FieldInfo((:nlight, 3)),
        :site_pos => FieldInfo((:nsite, 3), jlsize=(:nsite, XYZ)),
        :dof_jntid => FieldInfo((:nv,)),
        :mesh_texcoord => FieldInfo((:nmeshtexvert, 2)),
        :jnt_solimp => FieldInfo((:njnt, 5)),
        :name_skinadr => FieldInfo((:nskin,)),
        :mesh_vertadr => FieldInfo((:nmesh,)),
        :geom_matid => FieldInfo((:ngeom,)),
        :dof_bodyid => FieldInfo((:nv,)),
        :tex_width => FieldInfo((:ntex,)),
        :tendon_stiffness => FieldInfo((:ntendon,)),
        :body_jntadr => FieldInfo((:nbody,)),
        :tex_type => FieldInfo((:ntex,)),
        :jnt_dofadr => FieldInfo((:njnt,)),
        :jnt_range => FieldInfo((:njnt, 2)),
        :tendon_solref_lim => FieldInfo((:ntendon, 2)),
        :pair_solimp => FieldInfo((:npair, 5)),
        :cam_poscom0 => FieldInfo((:ncam, 3), jlsize=(:ncam, XYZ)),
        :cam_mat0 => FieldInfo((:ncam, 9)),
        :light_cutoff => FieldInfo((:nlight,)),
        :key_qvel => FieldInfo((:nkey, :nv)),
        :hfield_adr => FieldInfo((:nhfield,)),
        :site_bodyid => FieldInfo((:nsite,)),
        :skin_bonevertadr => FieldInfo((:nskinbone,)),
        :body_weldid => FieldInfo((:nbody,)),
        :geom_dataid => FieldInfo((:ngeom,)),
        :jnt_pos => FieldInfo((:njnt, 3), jlsize=(:njnt, XYZ)),
        :name_pairadr => FieldInfo((:npair,)),
        :site_type => FieldInfo((:nsite,)),
        :name_eqadr => FieldInfo((:neq,)),
        :tendon_length0 => FieldInfo((:ntendon,)),
        :dof_parentid => FieldInfo((:nv,)),
        :name_texadr => FieldInfo((:ntex,)),
        :light_dir => FieldInfo((:nlight, 3)),
        :dof_solref => FieldInfo((:nv, 2)),
        :skin_bonevertnum => FieldInfo((:nskinbone,)),
        :text_size => FieldInfo((:ntext,)),
        :mat_rgba => FieldInfo((:nmat, 4), jlsize=(:nmat, RGBA)),
        :body_sameframe => FieldInfo((:nbody,)),
        :exclude_signature => FieldInfo((:nexclude,)),
        :mesh_vert => FieldInfo((:nmeshvert, 3), jlsize=(:nmeshvert, XYZ)),
        :body_geomadr => FieldInfo((:nbody,)),
        :dof_M0 => FieldInfo((:nv,)),
        :mat_reflectance => FieldInfo((:nmat,)),
        :geom_quat => FieldInfo((:ngeom, 4)),
        :jnt_qposadr => FieldInfo((:njnt,)),
        :actuator_dynprm => FieldInfo((:nu, 10)),
        :pair_gap => FieldInfo((:npair,)),
        :actuator_gainprm => FieldInfo((:nu, 10)),
        :skin_bonevertid => FieldInfo((:nskinbonevert,)),
        :name_bodyadr => FieldInfo((:nbody,)),
        :cam_pos0 => FieldInfo((:ncam, 3), jlsize=(:ncam, XYZ)),
        :light_bodyid => FieldInfo((:nlight,)),
        :mesh_normal => FieldInfo((:nmeshvert, 3)),
        :name_hfieldadr => FieldInfo((:nhfield,)),
        :mesh_vertnum => FieldInfo((:nmesh,)),
        :mesh_graphadr => FieldInfo((:nmesh,)),
        :mesh_face => FieldInfo((:nmeshface, 3)),
        :skin_facenum => FieldInfo((:nskin,)),
        :tendon_num => FieldInfo((:ntendon,)),
        :hfield_data => FieldInfo((:nhfielddata,)),
        :text_adr => FieldInfo((:ntext,)),
        :cam_user => FieldInfo((:ncam, :nuser_cam)),
        :tendon_rgba => FieldInfo((:ntendon, 4)),
        :geom_pos => FieldInfo((:ngeom, 3)),
        :text_data => FieldInfo((:ntextdata,)),
        :geom_bodyid => FieldInfo((:ngeom,)),
        :skin_bonevertweight => FieldInfo((:nskinbonevert,)),
        :body_quat => FieldInfo((:nbody, 4)),
        :tendon_lengthspring => FieldInfo((:ntendon,)),
        :light_mode => FieldInfo((:nlight,)),
        :skin_bonebindquat => FieldInfo((:nskinbone, 4)),
        :name_tendonadr => FieldInfo((:ntendon,)),
        :light_diffuse => FieldInfo((:nlight, 3)),
        :tuple_size => FieldInfo((:ntuple,)),
        :skin_rgba => FieldInfo((:nskin, 4)),
        :geom_group => FieldInfo((:ngeom,)),
        :key_act => FieldInfo((:nkey, :na)),
        :geom_margin => FieldInfo((:ngeom,)),
        :body_iquat => FieldInfo((:nbody, 4)),
        :dof_invweight0 => FieldInfo((:nv,)),
        :cam_quat => FieldInfo((:ncam, 4)),
        :skin_vertadr => FieldInfo((:nskin,)),
        :actuator_gaintype => FieldInfo((:nu,)),
        :geom_solimp => FieldInfo((:ngeom, 5)),
        :dof_frictionloss => FieldInfo((:nv,)),
        :body_geomnum => FieldInfo((:nbody,)),
        :sensor_needstage => FieldInfo((:nsensor,)),
        :actuator_group => FieldInfo((:nu,)),
        :sensor_adr => FieldInfo((:nsensor,)),
        :eq_data => FieldInfo((:neq, 7)),
        :geom_user => FieldInfo((:ngeom, :nuser_geom)),
        :body_user => FieldInfo((:nbody, :nuser_body)),
        :skin_boneadr => FieldInfo((:nskin,)),
        :jnt_solref => FieldInfo((:njnt, 2)),
        :body_pos => FieldInfo((:nbody, 3)),
        :sensor_objtype => FieldInfo((:nsensor,)),
        :dof_Madr => FieldInfo((:nv,)),
        :site_quat => FieldInfo((:nsite, 4)),
        :sensor_datatype => FieldInfo((:nsensor,)),
        :dof_damping => FieldInfo((:nv,)),
        :site_user => FieldInfo((:nsite, :nuser_site)),
        :actuator_dyntype => FieldInfo((:nu,)),
        :pair_geom1 => FieldInfo((:npair,)),
        :light_specular => FieldInfo((:nlight, 3)),
        :jnt_axis => FieldInfo((:njnt, 3)),
        :mat_texid => FieldInfo((:nmat,)),
        :name_numericadr => FieldInfo((:nnumeric,)),
        :cam_ipd => FieldInfo((:ncam,)),
        :site_size => FieldInfo((:nsite, 3)),
        :geom_sameframe => FieldInfo((:ngeom,)),
        :tendon_limited => FieldInfo((:ntendon,)),
        :tendon_solref_fri => FieldInfo((:ntendon, 2)),
        :eq_obj2id => FieldInfo((:neq,)),
        :site_sameframe => FieldInfo((:nsite,)),
        :tex_adr => FieldInfo((:ntex,)),
        :tendon_damping => FieldInfo((:ntendon,)),
        :name_matadr => FieldInfo((:nmat,)),
        :geom_gap => FieldInfo((:ngeom,)),
        :sensor_dim => FieldInfo((:nsensor,)),
        :dof_armature => FieldInfo((:nv,)),
        :geom_conaffinity => FieldInfo((:ngeom,)),
        :name_actuatoradr => FieldInfo((:nu,)),
        :cam_targetbodyid => FieldInfo((:ncam,)),
        :geom_type => FieldInfo((:ngeom,)),
        :actuator_biastype => FieldInfo((:nu,)),
        :key_time => FieldInfo((:nkey,)),
        :body_subtreemass => FieldInfo((:nbody,)),
        :skin_inflate => FieldInfo((:nskin,)),
        :mat_shininess => FieldInfo((:nmat,)),
        :site_rgba => FieldInfo((:nsite, 4)),
        :jnt_user => FieldInfo((:njnt, :nuser_jnt)),
        :tendon_solimp_lim => FieldInfo((:ntendon, 5)),
        :cam_mode => FieldInfo((:ncam,)),
        :jnt_margin => FieldInfo((:njnt,)),
        :pair_signature => FieldInfo((:npair,)),
        :name_meshadr => FieldInfo((:nmesh,)),
        :wrap_prm => FieldInfo((:nwrap,)),
        :numeric_data => FieldInfo((:nnumericdata,)),
        :sensor_cutoff => FieldInfo((:nsensor,)),
        :site_matid => FieldInfo((:nsite,)),
        :hfield_size => FieldInfo((:nhfield, 4)),
        :sensor_user => FieldInfo((:nsensor, :nuser_sensor)),
        :actuator_trnid => FieldInfo((:nu, 2)),
        :name_camadr => FieldInfo((:ncam,)),
        :geom_solref => FieldInfo((:ngeom, 2)),
        :qpos_spring => FieldInfo((:nq,)),
        :cam_pos => FieldInfo((:ncam, 3)),
        :name_jntadr => FieldInfo((:njnt,)),
        :actuator_acc0 => FieldInfo((:nu,)),
        :jnt_stiffness => FieldInfo((:njnt,)),
        :body_jntnum => FieldInfo((:nbody,)),
        :actuator_length0 => FieldInfo((:nu,)),
        :pair_geom2 => FieldInfo((:npair,)),
        :name_siteadr => FieldInfo((:nsite,)),
        :actuator_trntype => FieldInfo((:nu,)),
        :dof_simplenum => FieldInfo((:nv,)),
        :pair_friction => FieldInfo((:npair, 5)),
        :name_keyadr => FieldInfo((:nkey,)),
        :tendon_user => FieldInfo((:ntendon, :nuser_tendon)),
        :actuator_ctrlrange => FieldInfo((:nu, 2)),
        :body_simple => FieldInfo((:nbody,)),
        :actuator_biasprm => FieldInfo((:nu, 10)),
        :tendon_frictionloss => FieldInfo((:ntendon,)),
        :actuator_cranklength => FieldInfo((:nu,)),
        :light_active => FieldInfo((:nlight,)),
        :mat_specular => FieldInfo((:nmat,)),
        :numeric_size => FieldInfo((:nnumeric,)),
        :tendon_adr => FieldInfo((:ntendon,)),
        :actuator_ctrllimited => FieldInfo((:nu,)),
        :geom_rgba => FieldInfo((:ngeom, 4)),
        :light_directional => FieldInfo((:nlight,)),
        :body_parentid => FieldInfo((:nbody,)),
        :light_dir0 => FieldInfo((:nlight, 3)),
        :tex_height => FieldInfo((:ntex,)),
        :pair_dim => FieldInfo((:npair,)),
        :mat_emission => FieldInfo((:nmat,)),
        :eq_solref => FieldInfo((:neq, 2)),
        :mesh_graph => FieldInfo((:nmeshgraph,)),
        :jnt_group => FieldInfo((:njnt,)),
        :light_poscom0 => FieldInfo((:nlight, 3)),
        :skin_texcoordadr => FieldInfo((:nskin,)),
        :mat_texrepeat => FieldInfo((:nmat, 2)),
        :jnt_type => FieldInfo((:njnt,)),
        :mat_texuniform => FieldInfo((:nmat,)),
        :geom_rbound => FieldInfo((:ngeom,)),
        :qpos0 => FieldInfo((:nq,)),
        :light_targetbodyid => FieldInfo((:nlight,)),
        :actuator_lengthrange => FieldInfo((:nu, 2)),
        :geom_contype => FieldInfo((:ngeom,)),
        :site_group => FieldInfo((:nsite,)),
        :pair_solref => FieldInfo((:npair, 2)),
        :light_pos0 => FieldInfo((:nlight, 3)),
        :tuple_objprm => FieldInfo((:ntupledata,)),
        :tendon_matid => FieldInfo((:ntendon,)),
        :actuator_user => FieldInfo((:nu, :nuser_actuator)),
        :actuator_gear => FieldInfo((:nu, 6)),
        :geom_condim => FieldInfo((:ngeom,)),
        :eq_obj1id => FieldInfo((:neq,)),
        :numeric_adr => FieldInfo((:nnumeric,)),
        :skin_matid => FieldInfo((:nskin,)),
        :tendon_width => FieldInfo((:ntendon,)),
        :tendon_solimp_fri => FieldInfo((:ntendon, 5)),
        :eq_solimp => FieldInfo((:neq, 5)),
        :name_sensoradr => FieldInfo((:nsensor,)),
        :body_rootid => FieldInfo((:nbody,)),
        :geom_friction => FieldInfo((:ngeom, 3)),
        :skin_vertnum => FieldInfo((:nskin,)),
        :name_tupleadr => FieldInfo((:ntuple,)),
        :pair_margin => FieldInfo((:npair,)),
        :wrap_objid => FieldInfo((:nwrap,)),
        :light_exponent => FieldInfo((:nlight,)),
    ),
    :mjData => Dict(
        :qLD => FieldInfo((:nM,)),
        :efc_aref => FieldInfo((:njmax,)),
        :efc_type => FieldInfo((:njmax,)),
        :userdata => FieldInfo((:nuserdata,)),
        :efc_force => FieldInfo((:njmax,)),
        :qfrc_unc => FieldInfo((:nv,)),
        :ten_wrapnum => FieldInfo((:ntendon,)),
        :cinert => FieldInfo((:nbody, 10)),
        :efc_JT_rowadr => FieldInfo((:nv,)),
        :ximat => FieldInfo((:nbody, 9)),
        :actuator_moment => FieldInfo((:nu, :nv)),
        :efc_AR_rownnz => FieldInfo((:njmax,)),
        :efc_J_rowsuper => FieldInfo((:njmax,)),
        :efc_J_rowadr => FieldInfo((:njmax,)),
        :efc_id => FieldInfo((:njmax,)),
        :qacc_unc => FieldInfo((:nv,)),
        :efc_margin => FieldInfo((:njmax,)),
        :mocap_quat => FieldInfo((:nmocap, 4), jlsize=(:nmocap, Quat)),
        :ten_length => FieldInfo((:ntendon,)),
        :xipos => FieldInfo((:nbody, 3)),
        :actuator_velocity => FieldInfo((:nu,)),
        :light_xpos => FieldInfo((:nlight, 3)),
        :xmat => FieldInfo((:nbody, 9), jlsize=(:nbody, frame()...)),
        :act => FieldInfo((:na,)),
        :efc_D => FieldInfo((:njmax,)),
        :qfrc_applied => FieldInfo((:nv,)),
        :ten_wrapadr => FieldInfo((:ntendon,)),
        :efc_frictionloss => FieldInfo((:njmax,)),
        :efc_JT_rownnz => FieldInfo((:nv,)),
        :cdof_dot => FieldInfo((:nv, 6)),
        :ten_J_rowadr => FieldInfo((:ntendon,)),
        :cdof => FieldInfo((:nv, 6)),
        :efc_JT => FieldInfo((:nv, :njmax)),
        :efc_b => FieldInfo((:njmax,)),
        :qfrc_bias => FieldInfo((:nv,)),
        :qacc => FieldInfo((:nv,)),
        :qfrc_passive => FieldInfo((:nv,)),
        :efc_AR => FieldInfo((:njmax, :njmax)),
        :geom_xpos => FieldInfo((:ngeom, 3)),
        :mocap_pos => FieldInfo((:nmocap, 3)),
        :qacc_warmstart => FieldInfo((:nv,)),
        :efc_KBIP => FieldInfo((:njmax, 4)),
        :cam_xpos => FieldInfo((:ncam, 3)),
        :qvel => FieldInfo((:nv,)),
        :ctrl => FieldInfo((:nu,)),
        :subtree_com => FieldInfo((:nbody, 3)),
        :efc_pos => FieldInfo((:njmax,)),
        :site_xmat => FieldInfo((:nsite, 9)),
        :wrap_xpos => FieldInfo((:nwrap, 6)),
        :qM => FieldInfo((:nM,)),
        :cam_xmat => FieldInfo((:ncam, 9)),
        :geom_xmat => FieldInfo((:ngeom, 9)),
        :xaxis => FieldInfo((:njnt, 3)),
        :efc_diagApprox => FieldInfo((:njmax,)),
        :crb => FieldInfo((:nbody, 10)),
        :efc_JT_colind => FieldInfo((:nv, :njmax)),
        :efc_J => FieldInfo((:njmax, :nv)),
        :qfrc_constraint => FieldInfo((:nv,)),
        :contact => FieldInfo((:nconmax,), dsize=(:ncon,)),
        :efc_state => FieldInfo((:njmax,)),
        :wrap_obj => FieldInfo((:nwrap,)),
        :qfrc_inverse => FieldInfo((:nv,)),
        :cfrc_int => FieldInfo((:nbody, 6)),
        :sensordata => FieldInfo((:nsensordata,)),
        :qpos => FieldInfo((:nq,)),
        :qfrc_actuator => FieldInfo((:nv,)),
        :cacc => FieldInfo((:nbody, 6)),
        :light_xdir => FieldInfo((:nlight, 3), jlsize=(:nlight, XYZ)),
        :act_dot => FieldInfo((:na,)),
        :efc_AR_rowadr => FieldInfo((:njmax,)),
        :cfrc_ext => FieldInfo((:nbody, 6)),
        :efc_J_colind => FieldInfo((:njmax, :nv)),
        :efc_AR_colind => FieldInfo((:njmax, :njmax)),
        :actuator_force => FieldInfo((:nu,)),
        :site_xpos => FieldInfo((:nsite, 3), jlsize=(:nsite, XYZ)),
        :subtree_angmom => FieldInfo((:nbody, 3)),
        :qLDiagSqrtInv => FieldInfo((:nv,)),
        :xquat => FieldInfo((:nbody, 4), jlsize=(:nbody, Quat)),
        :ten_J_colind => FieldInfo((:ntendon, :nv)),
        :efc_JT_rowsuper => FieldInfo((:nv,)),
        :xpos => FieldInfo((:nbody, 3), jlsize=(:nbody, XYZ)),
        :cvel => FieldInfo((:nbody, 6)),
        :qLDiagInv => FieldInfo((:nv,)),
        :efc_J_rownnz => FieldInfo((:njmax,)),
        :xanchor => FieldInfo((:njnt, 3), jlsize=(:njnt, XYZ)),
        :ten_velocity => FieldInfo((:ntendon,)),
        :ten_J_rownnz => FieldInfo((:ntendon,)),
        :xfrc_applied => FieldInfo((:nbody, 6)),
        :ten_J => FieldInfo((:ntendon, :nv)),
        :efc_vel => FieldInfo((:njmax,)),
        :subtree_linvel => FieldInfo((:nbody, 3)),
        :actuator_length => FieldInfo((:nu,)),
        :efc_R => FieldInfo((:njmax,)),
    ),
)


text:
using BIOMUtils: BIOM
using HDF5
using Test
using SparseArrays

@testset "io" begin

end

@testset "parse" begin
    @testset "func" begin
        @test 1 + 1 == 2
    end
end

@testset "process" begin
    @testset "func" begin
        @test 1 + 1 == 2
    end
end

@testset "types" begin
    @testset "BIOM" begin
        data = [
            1 1 0 0
            0 1 1 0
            0 0 1 1
        ]
        sampleIds = ["S1", "S2", "S3"]
        otuIds = ["otu1", "otu2", "otu3", "otu4"]
        taxonomy = [
            "k__Bacteria"            "k__Bacteria"           "k__Bacteria"           "k__Bacteria"
            "p__Proteobacteria"      "p__Proteobacteria"     "p__Proteobacteria"     "p__Proteobacteria"
            "c__Alphaproteobacteria" "c__Betaproteobacteria" "c__Betaproteobacteria" "c__Gammaproteobacteria"
            "o__MND"                 "o__MND1"               "o__MND2"               "o__"
            "f__mnd"                 "f__mnd"                "f__"                   "f__"
            "g__mnd"                 "g__"                   "g__"                   "g__"
            "s__"                    "s__"                   "s__"                   "s__"
        ]
        #= Test Types =#
        @test BIOM(data) isa BIOM
        @test BIOM(data).data isa SparseMatrixCSC
        @test BIOM(data).sample isa NamedTuple
        @test BIOM(data).observation isa NamedTuple
        @test BIOM(data).sample.ids isa Array
        @test BIOM(data).sample.metadata isa Dict
        @test BIOM(data).sample.groupmetadata isa Dict
        @test BIOM(data).observation.ids isa Array
        @test BIOM(data).observation.metadata isa Dict
        @test BIOM(data).observation.groupmetadata isa Dict

        @test BIOM(data).sample == BIOM(sparse(data)).sample

    end
end

@testset "utils" begin
    @testset "func" begin
        @test 1 + 1 == 2
    end
end


text:
for ch in [1, 3, 16, 64]
    x = rand(Float32, 64, 64, ch, 64)
    model = Conv((3,3), ch=>ch)
    println("CPU ch=$ch")
    run_benchmark(model, x, cuda=false)
    println("CUDA ch=$ch")
    run_benchmark(model, x, cuda=true)
end


text:
@test mysignif(SpecificV(0.255_837_018E2, 650.0), TestDigits) ≈ 0.002
@test mysignif(SpecificV(0.255_837_018E2u"MPa", 650.0u"K"), TestDigits) ≈ 0.002u"m^3/kg"

@test mysignif(SpecificV(0.222_930_643E2, 650.0), TestDigits) ≈ 0.005
@test mysignif(SpecificV(0.222_930_643E2u"MPa", 650.0u"K"), TestDigits) ≈ 0.005u"m^3/kg"

@test mysignif(SpecificV(0.783_095_639E2, 750.0), TestDigits) ≈ 0.002
@test mysignif(SpecificV(0.783_095_639E2u"MPa", 750.0u"K"), TestDigits) ≈ 0.002u"m^3/kg"

@test SpecificH(0.255_837_018E2, 650.0) ≈ 0.186_343_019E4
@test SpecificH(0.255_837_018E2u"MPa", 650.0u"K") ≈ 0.186_343_019E4u"kJ/kg"

@test SpecificH(0.222_930_643E2, 650.0) ≈ 0.237_512_401E4
@test SpecificH(0.222_930_643E2u"MPa", 650.0u"K") ≈ 0.237_512_401E4u"kJ/kg"

@test SpecificH(0.783_095_639E2, 750.0) ≈ 0.225_868_845E4
@test SpecificH(0.783_095_639E2u"MPa", 750.0u"K") ≈ 0.225_868_845E4u"kJ/kg"

@test SpecificU(0.255_837_018E2, 650.0) ≈ 0.181_226_279E4
@test SpecificU(0.255_837_018E2u"MPa", 650.0u"K") ≈ 0.181_226_279E4u"kJ/kg"

@test SpecificU(0.222_930_643E2, 650.0) ≈ 0.226_365_868E4
@test SpecificU(0.222_930_643E2u"MPa", 650.0u"K") ≈ 0.226_365_868E4u"kJ/kg"

@test SpecificU(0.783_095_639E2, 750.0) ≈ 0.210_206_932E4
@test SpecificU(0.783_095_639E2u"MPa", 750.0u"K") ≈ 0.210_206_932E4u"kJ/kg"

@test SpecificS(0.255_837_018E2, 650.0) ≈ 0.405_427_273E1
@test SpecificS(0.255_837_018E2u"MPa", 650.0u"K") ≈ 0.405_427_273E1u"kJ/kg/K"

@test SpecificS(0.222_930_643E2, 650.0) ≈ 0.485_438_792E1
@test SpecificS(0.222_930_643E2u"MPa", 650.0u"K") ≈ 0.485_438_792E1u"kJ/kg/K"

@test SpecificS(0.783_095_639E2, 750.0) ≈ 0.446_971_906E1
@test SpecificS(0.783_095_639E2u"MPa", 750.0u"K") ≈ 0.446_971_906E1u"kJ/kg/K"

@test SpecificCP(0.255_837_018E2, 650.0) ≈ 0.138_935_717E2
@test SpecificCP(0.255_837_018E2u"MPa", 650.0u"K") ≈ 0.138_935_717E2u"kJ/kg/K"

@test SpecificCP(0.222_930_643E2, 650.0) ≈ 0.446_579_373E2
@test SpecificCP(0.222_930_643E2u"MPa", 650.0u"K") ≈ 0.446_579_373E2u"kJ/kg/K"

@test SpecificCP(0.783_095_639E2, 750.0) ≈ 0.634_165_359E1
@test SpecificCP(0.783_095_639E2u"MPa", 750.0u"K") ≈ 0.634_165_359E1u"kJ/kg/K"

@test SpeedOfSound(0.255_837_018E2, 650.0) ≈ 0.502_005_554E3
@test SpeedOfSound(0.255_837_018E2u"MPa", 650.0u"K") ≈ 0.502_005_554E3u"m/s"

@test SpeedOfSound(0.222_930_643E2, 650.0) ≈ 0.383_444_594E3
@test SpeedOfSound(0.222_930_643E2u"MPa", 650.0u"K") ≈ 0.383_444_594E3u"m/s"

@test SpeedOfSound(0.783_095_639E2, 750.0) ≈ 0.760_696_041E3
@test SpeedOfSound(0.783_095_639E2u"MPa", 750.0u"K") ≈ 0.760_696_041E3u"m/s"


text:
#=
mri_objects.jl
2019-11 Connor Martin
based on mri_objects.m
Copyright 2007-6-28, Jeff Fessler, University of Michigan
=#

export mri_objects

using Plots
using FFTW: fft, fftshift, ifftshift
#using MIRT: max_percent_diff, jinc, rect
#using MIRT: jim, image_geom, prompt, interp1

RealArray = AbstractArray{<:Real}
PairList = Vector{<:Tuple} # Vector{Tuple{Symbol,<:Any}}

# struct for wrapper definition
struct MIRT_mri_object
    pairs::PairList
    image::Function
    kspace::Function
end

# the two primary functions
# nvp is a vector of tuples (type::Symbol, parameter vector)
# for 2d, just make the third coord 0

function mri_objects_image(nvp::PairList, x::RealArray, y::RealArray, z ;
        dx::Real=0, dy::Real=0, dz::Real=0)
    out = zeros(size(x)...)
    for pair in nvp
        out += imageFuns[pair[1]](pair[2], x, y, z ; dx=dx, dy=dy, dz=dz)
    end
    return out
end

function mri_objects_kspace(nvp::PairList, u::RealArray, v::RealArray, w)
    out = zeros(size(u)...)
    for pair in nvp
        out += kspaceFuns[pair[1]](pair[2], u, v, w)
    end
    return out
end

# 2D
function mri_objects_image(data, x::RealArray, y::RealArray ; kwargs...)
    return mri_objects_image(data, x, y, 0 ; kwargs...)
end

function mri_objects_kspace(data, u::RealArray, v::RealArray)
    return mri_objects_kspace(data, u, v, 0)
end

# The big chunk of functions for each type (image & kspace)

function mri_objects_trap(z, dz, len)
    dz = abs(dz)
    (dz == 0 || len < dz) && return rect.(z / len)
    interp1([-len/2-dz/2, -len/2+dz/2, len/2-dz/2, len/2+dz/2], [0, 1, 1, 0], z)
end

function mri_objects_image_gauss2(data::RealArray, x::RealArray, y::RealArray,
        z ; kwargs...)
    size(data,2) != 5 && throw(DomainError(data,"gauss2 requires 5 parameters"))
    z = zeros(size(data,1))
    zw = z .+ Inf
    inputParams = [data[:,1:2] z data[:,3:4] zw data[:,5]]
    return mri_objects_image_gauss3(inputParams,x,y,0)
end

function mri_objects_image_gauss3(data::RealArray, x::RealArray, y::RealArray,
        z ; kwargs...)
    size(data,2) != 7 && throw(DomainError(data,"gauss3 requires 7 parameters"))
    out = zeros(size(x))
    for i=1:size(data,1)
        par = data[i,:]
        xc = par[1]
        yc = par[2]
        zc = par[3]
        xw = par[4] ./ sqrt(log(256)) .* sqrt(2*pi)
        yw = par[5] ./ sqrt(log(256)) .* sqrt(2*pi)
        zw = par[6] ./ sqrt(log(256)) .* sqrt(2*pi)
        param1 = exp.(-pi .* ((x .- xc) ./ xw).^2)
        param2 = exp.(-pi .* ((y .- yc) ./ yw).^2)
        param3 = exp.(-pi .* ((z .- zc) ./ zw).^2)
        out = out .+ par[7] .* param1 .* param2 .* param3
    end
    return out
end

function mri_objects_image_circ2(data::RealArray, x::RealArray, y::RealArray,
        zr ; kwargs...)
    size(data,2) != 4 && throw(DomainError(data,"circ2 requires 4 parameters"))
    z = zeros(size(data,1))
    inputParams = [data[:,1:2] z data[:,3] (z .+ 1) data[:,4]]
    return mri_objects_image_cyl3(inputParams, x,y,zr ; kwargs...)
end

function mri_objects_image_cyl3(data::RealArray, x::RealArray, y::RealArray,
        z ; dx::Real=0, dy::Real=0, dz::Real=0)
    size(data,2) != 6 && throw(DomainError(data,"cyl3 requires 6 parameters"))
    out = zeros(size(x))
    for i = 1:size(data,1)
        par = data[i,:]
        xc = par[1]
        yc = par[2]
        zc = par[3]
        rad = par[4]
        len = par[5]
        out += par[6] * ((x .- xc).^2 + (y .- yc).^2 .< rad^2) .*
            mri_objects_trap(z .- zc, dz, len)
    end
    return out
end

function mri_objects_image_rect2(data::RealArray, x::RealArray, y::RealArray,
        zr ; dx::Real=0, dy::Real=0, dz::Real=0)
    size(data,2) != 5 && throw(DomainError(data,"rect2 requires 5 parameters"))
    z = zeros(size(data,1))
    inputParams = [data[:,1:2] z data[:,3:4] (z .+ 1) data[:,5]]
    return mri_objects_image_rect3(inputParams, x,y,zr ; dx=dx, dy=dy, dz=dz)
end

function mri_objects_image_rect3(data::RealArray, x::RealArray, y::RealArray,
        z ; dx::Real=0, dy::Real=0, dz::Real=0)
    out = zeros(size(x))
    for i=1:size(data,1)
        par = data[i,:]
        xc = par[1]
        yc = par[2]
        zc = par[3]
        xw = par[4]
        yw = par[5]
        zw = par[6]
        out = out .+ par[7] *
            mri_objects_trap(x .- xc, dx, xw) .*
            mri_objects_trap(y .- yc, dy, yw) .*
            mri_objects_trap(z .- zc, dz, zw)
    end
    return out
end

function mri_objects_image_dirac2(data::RealArray, x::RealArray, y::RealArray,
        z ; kwargs...)
    size(data,2) != 3 && throw(DomainError(data,"dirac2 requires 3 parameters"))
    #@show data
    inputParams = [data[:,1:2] zeros(size(data,1),1) data[:,3]]
    return mri_objects_image_dirac3(inputParams, x,y,z ; kwargs...)
end

function mri_objects_image_dirac3(data::RealArray, x::RealArray, y::RealArray,
        z ; kwargs...)
    size(data,2) != 4 && throw(DomainError(data,"dirac3 requires 4 parameters"))
    out = zeros(size(x))
    for i=1:size(data,1)
        out += data[i,4] .* ((x .== data[i,1]) .* (y .== data[i,2]) .* (z .== data[i,3]))
    end
#TODO: what does out(out != 0) = inf, warn do? Ask fess. todo: test
    replace!(x -> x==0 ? 0 : Inf, out) #no temp array
    return out
end

function mri_objects_kspace_gauss2(data::RealArray, u::RealArray, v::RealArray, w)
    size(data,2) != 5 && throw(DomainError(data,"gauss2 requires 5 parameters"))
    z = zeros(size(data,1))
    zw = (1 .+ z) .* sqrt(log(256)) ./ sqrt(2*pi)
    inputParams = [data[:,1:2] z data[:,3:4] zw data[:,5]]
    return mri_objects_kspace_gauss3(inputParams,u,v,w)
end

function mri_objects_kspace_gauss3(data::RealArray, u::RealArray, v::RealArray, w)
    size(data,2) != 7 && throw(DomainError(data,"gauss3 requires 7 parameters"))
    out = 0
    for ii=1:size(data,1)
        par = data[ii,:]
        xc = par[1]
        yc = par[2]
        zc = par[3]
        xw = par[4] / sqrt(log(256)) * sqrt(2*pi)
        yw = par[5] / sqrt(log(256)) * sqrt(2*pi)
        zw = par[6] / sqrt(log(256)) * sqrt(2*pi)
        out = out .+ par[7] * xw * yw * zw *
            exp.(-pi * (u .* xw) .^2) .*
            exp.(-pi * (v .* yw) .^2) .*
            exp.(-pi * (w .* zw) .^2) .*
            exp.(-2im * pi * (u .* xc .+ v .* yc .+ w .* zc))
    end
    return out
end

function mri_objects_kspace_circ2(data::RealArray, u::RealArray, v::RealArray, w)
    size(data,2) != 4 && throw(DomainError(data,"circ2 requires 4 parameters"))
    z = zeros(size(data,1))
    inputParams = [data[:,1:2] z data[:,3] (z .+ 1) data[:,4]]
    return mri_objects_kspace_cyl3(inputParams,u,v,w)
end

function mri_objects_kspace_cyl3(data::RealArray, u::RealArray, v::RealArray, w)
    size(data,2) != 6 && throw(DomainError(data,"cyl3 requires 6 parameters"))
    out = zeros(size(u))
    for i = 1:size(data,1)
        par = data[i,:]
        xc = par[1]
        yc = par[2]
        zc = par[3]
        rad = par[4]
        len = par[5]
        out = out .+ par[6] * rad^2 * 4 .*
            jinc.(2 * sqrt.(u.^2 .+ v.^2) .* rad) .* (len .* sinc.(w .* len)) .*
            exp.(-2im * pi * (u .* xc .+ v .* yc .+ w .* zc))
    end
    return out
end

function mri_objects_kspace_rect2(data::RealArray, u::RealArray, v::RealArray, w)
    size(data,2) != 5 && throw(DomainError(data,"rect2 requires 5 parameters"))
    z = zeros(size(data,1),1)
    inputParams = [data[:,1:2] z data[:,3:4] (z .+ 1) data[:,5]]
    return mri_objects_kspace_rect3(inputParams,u,v,0)
end

function mri_objects_kspace_rect3(data::RealArray, u::RealArray, v::RealArray, w)
    size(data,2) != 7 && throw(DomainError(data,"rect3 requires 7 parameters"))
    out = 0
    for ii=1:size(data,1)
        par = data[ii,:]
        xc = par[1]
        yc = par[2]
        zc = par[3]
        xw = par[4]
        yw = par[5]
        zw = par[6]
        out = out .+ par[7] * xw * yw * zw *
            sinc.(u .* xw) .* sinc.(v .* yw) .* sinc.(w .* zw) .*
            exp.(-2im * pi * (u .* xc + v .* yc .+ w .* zc))
    end
    return out
end

function mri_objects_kspace_dirac2(data::RealArray, u::RealArray, v::RealArray, w)
    size(data,2) != 3 && throw(DomainError(data,"dirac2 requires 3 parameters"))
    inputParams = [data[:,1:2] zeros(size(data,1),1) data[:,3]]
    return mri_objects_kspace_dirac3(inputParams,u,v,w)
end

function mri_objects_kspace_dirac3(data::RealArray, u::RealArray, v::RealArray, w)
    size(data,2) != 4 && throw(DomainError(data,"dirac3 requires 4 parameters"))
    out = zeros(size(u))
    for i=1:size(data,1)
        out += data[i,4] * exp.(-2im * pi * (u .* data[i,1] .+ v .* data[i,2] .+ w .* data[i,3]))
    end
    return out
end

function mri_objects_case1( ; unit::Symbol = :mm)
    rp = [ # rect2
     0 0.0    200 200  1;
    -50 -50  40  40  1;
     50 -50  20  20  1;
     0 50    50  50  1;
    ]
    gp = [ # gauss2
     -70 0.0  1 1  1;
     -60 0  2 2  1;
     -50 0  3 3  1;
     -40 0  4 4  1;
     -20 0  5 5  1;
       0 0  6 6  1;
      20 0  7 7  1;
      50 0  8 8  1;
    ]
    rp
    if unit === :cm
        rp[:,1:4] ./= 10
        gp[:,1:4] ./= 10
    end
    return [(:rect2,rp), (:gauss2,gp)] # PairList
end

# image function dictionary
imageFuns = Dict([
    (:dirac2, mri_objects_image_dirac2),
    (:dirac3, mri_objects_image_dirac3),
    (:rect2, mri_objects_image_rect2),
    (:rect3, mri_objects_image_rect3),
    (:gauss2, mri_objects_image_gauss2),
    (:gauss3, mri_objects_image_gauss3),
    (:circ2, mri_objects_image_circ2),
    (:cyl3, mri_objects_image_cyl3),
])

# kspace function dictionary
kspaceFuns = Dict([
    (:dirac2, mri_objects_kspace_dirac2),
    (:dirac3, mri_objects_kspace_dirac3),
    (:rect2, mri_objects_kspace_rect2),
    (:rect3, mri_objects_kspace_rect3),
    (:gauss2, mri_objects_kspace_gauss2),
    (:gauss3, mri_objects_kspace_gauss3),
    (:circ2, mri_objects_kspace_circ2),
    (:cyl3, mri_objects_kspace_cyl3),
])

"""
    st = mri_objects([(:type1, params1), (:type2, params2), ...])
Generate object that describes image-domain objects and Fourier domain spectra
of simple structures such as rectangles, disks, superpositions thereof.
These functions are useful for simple "idealized" MRI simulations
where the data is modeled as analytical Fourier samples,
i.e., no field inhomogeneity and no relaxation effects.

in
- `[(type, params), ...]` e.g. `[(:rect2, params), (:gauss3, params), ...]`

type, params:
- `:dirac2 [N 3] [xcent ycent value]`
- `:dirac3 [N 4] [xcent ycent zcent value]`
- `:rect2 [N 5] [xcent ycent xwidth ywidth value]`
- `:rect3 [N 7] [xcent ycent zcent xwidth ywidth zwidth value]`
- `:gauss2 [N 5] [xcent ycent xwidth ywidth value]`
- `:gauss3 [N 7] [xcent ycent zcent xwidth ywidth zwidth value]`
- `:circ2 [N 5] [xcent ycent rad value]`
- `:cyl3 [N 6] [xcent ycent zcent xrad zwidth value]`

All types must be 2D or 3D, not mixed.

out
- `st` struct
* `st.image(x,y)`  returns 2D image-domain (sampled) picture
* `st.image(x,y,z)` same but 3D
* `st.kspace(u,v)` returns 2D Fourier-space samples
* `st.kspace(u,v,w)` same but 3D
"""
function mri_objects(nvp::PairList)
    type2 = (:dirac2, :rect2, :gauss2, :circ2)
    type3 = (:dirac3, :rect3, :gauss3, :cyl3)
    is3 = nvp[1][1] ∈ type3
    if is3
        for pair in nvp
            (pair[1] ∉ type3) && throw("mix of 2d and 3d types")
        end

        return MIRT_mri_object(nvp,
              (x,y,z ; kwargs...) -> mri_objects_image(nvp, x,y,z ; kwargs...),
              (u,v,w ; kwargs...) -> mri_objects_kspace(nvp, u,v,w ; kwargs...))
    else
        for pair in nvp
            (pair[1] ∉ type2) && throw("mix of 2d and 3d types")
        end
        return MIRT_mri_object(nvp,
            (x,y ; kwargs...) -> mri_objects_image(nvp, x,y ; kwargs...),
            (u,v ; kwargs...) -> mri_objects_kspace(nvp, u,v ; kwargs...))
    end
end

"""
    mri_objects((type,params)) for a single tuple
Here the `params` can be a single row or column vector.
"""
mri_objects(pair::Tuple) = mri_objects([(pair[1],reshape(pair[2],1,:)),])

# special cases

function mri_objects_case4(fov::AbstractVector{<:Real} ; unit::Symbol=:mm)
    length(fov) != 3 && throw("fov must have length 3")
    cp = [0 0 0 0.4*fov[1] fov[3] 1] # cyl3
    rp = Array{Float16}([ # rect3
      -50 -50  0  40 40 40  1;
       50 -50 40  20 20 50  1;
       0  50 -40  30 30 60  1;
    ])
    rp[:,1:3] = rp[:,1:3]/256 * fov[1]
    rp[:,4:6] = rp[:,4:6]/256 * fov[1]
    gp = Array{Float32}([ # gauss3
       -70 0 0  1 1 1  1;
       -60 0 0  2 2 2  1;
       -50 0 0  3 3 3  1;
       -40 0 0  4 4 4  1;
       -20 0 0  5 5 5  1;
        00 0 0  6 6 6  1;
        20 0 0  7 7 7  1;
        50 0 0  8 8 8  1;
    ])
    gp[:,1:3] = gp[:,1:3]/256 * fov[1]
    gp[:,4:6] = gp[:,4:6]/256 * fov[1]
    if (unit === :cm)
        rp[:,1:6] /= 10
        gp[:,1:6] /= 10
    end
    return [(:cyl3, cp), (:rect3, rp), (:gauss3, gp)] #PairList
end

# 2D test
#NOTE: Dirac will not plot very well. Just make sure it prints a gray image, and that's fine. FFT will also be screwed up.
function mri_objects_test2()
    ig = image_geom(nx = 2^7, ny = 2^7 + 2, dx = 4, offsets = :dsp)

    shift = [0.1, 0.2] .* ig.fovs
    sizes = [.15, .1] .* ig.fovs
    tests = ( # 3 separate tests:
                (:circ2, [shift' sizes[1] 2]),
                (:gauss2, [shift' sizes' 2]),
                (:rect2, [shift' sizes' 2]),
                (:dirac2,[shift' 2])
            )
    ntest = length(tests)
    pl = Matrix{Plots.Plot}(undef, ntest,3)

    for (i,test) = enumerate(tests)
        st = mri_objects(test)
        i2 = st.image(ig.xg, ig.yg, dx=ig.dx, dy=ig.dy)
        s2 = ifftshift(fft(fftshift(i2))) * abs(ig.dx * ig.dy)
        f2 = st.kspace(ig.fg...)
        @show max_percent_diff(f2,s2)
    #   @show round.(maximum.([i2, abs.(f2), abs.(s2)]), digits=2)
        pl[i,1] = jim(i2, title = "$(test[1])")
        pl[i,2] = jim(abs.(s2), title = "|fft|")
        pl[i,3] = jim(abs.(f2), title = "|kspace|")
    end
    plot(pl...)
    return plot(pl...)
end

function mri_objects_test_case1()
    ig = image_geom(nx = 2^7, ny = 2^7 + 2, dx = 4, offsets = :dsp)
    xt = mri_objects(:case1,unit = :cm).image(ig.xg, ig.yg)
    jim(xt, title = "case1cm")
end

function mri_objects_test_case4()
    ig = image_geom(nx = 2^5, ny = 2^5 + 2,nz = 2^5, dx = 1,dz = .5, offsets = :dsp)
    xt = mri_objects(:case4,unit = :cm).image(ig.xg, ig.yg, ig.zg)
    jim(xt, title = "case4cm")
end

# 3D test
#NOTE: Dirac will not plot very well. Just make sure it prints a gray image, and that's fine. FFT will also be screwed up.
function mri_objects_test3()
    ig = image_geom(nx = 2^7, nz = 2^5, offsets = :dsp, dx = 4, dz = 3)
    shift = [0.1, 0.2, 0.3] .* ig.fovs
    sizes = [0.15, 0.1, 0.2] .* ig.fovs
    tests = ((:cyl3, [shift' sizes[1] sizes[3] 2]),
             (:gauss3, [shift' sizes' 2]),
             (:rect3, [shift' sizes' 2]),
             (:dirac3,[shift' 2]))
    ntest = length(tests)
    pl = Matrix{Plots.Plot}(undef, ntest,3)
    for (i,test) in enumerate(tests)
        st = mri_objects(test)
        i3 = st.image(ig.xg,ig.yg,ig.zg, dx=ig.dx, dy=ig.dy, dz=ig.dz)
        s3 = ifftshift(fft(fftshift(i3)))
        s3 *= abs(ig.dx * ig.dy * ig.dz)
        fg = ig.fg
        f3 = st.kspace(fg...)
        @show max_percent_diff(f3,s3)
        pl[i,1] = jim(i3, title = "$(test[1])")
        pl[i,2] = jim(abs.(s3), title = "fft")
        pl[i,3] = jim(abs.(f3), title = "kspace")
        plot(pl[i,:]...)
        prompt()
    end
    return plot(pl...)
end

"""
    mri_objects(key::Symbol ; fov::Real=22, unit::Symbol=:mm)
Special cases

In:
- `fov` Store the following parameter in fov.
- `unit` specify units; default :mm
     Currently, only `:cm` (centimeters) supported as an alternate unit

key choices:
- `:test` run test suite (no other params)
- `:case1` predefined 2D test case
- `:case4` predefined 3D test case
"""
function mri_objects(key::Symbol ; fov::Real=22, unit::Symbol=:mm)
    key != :test && key != :case1 && key != :case4 && throw("bad key $key")
    if key == :test
        return mri_objects_test()
    elseif key == :case1
        return mri_objects(mri_objects_case1( ; unit=unit))
    elseif key == :case4
        return mri_objects(mri_objects_case4([fov,fov,fov] ; unit=unit))
    end
end

# plot the trapezoid
function mri_objects_trap_test()
    len = 6
    dz = 2
    z = LinRange(-8,8,201)
    trap = mri_objects_trap(z, dz, len)
    trap0 = mri_objects_trap(z, 0, len)
    plot(z, trap, label="trap")
    plot!(z, trap0, label="rect")
    plot!(xtick = [0 -len/2 len/2-dz/2 len/2 len/2+dz/2])
end

function mri_objects_test()
    @test mri_objects_trap_test() isa Plots.Plot
    prompt()

    @test mri_objects_test2() isa Plots.Plot
    prompt()

	@test mri_objects_test_case1() isa Plots.Plot
    prompt()

    @test mri_objects_test3() isa Plots.Plot

    @test mri_objects_test_case4() isa Plots.Plot
    prompt()

    true
end


text:
if Sys.isapple()
    using Homebrew
    Homebrew.brew(`cask install inkscape`)
    brew_path = joinpath(Homebrew.prefix(), "bin")
    ENV["PATH"] = string(ENV["PATH"], ":", brew_path)
end


text:
module PWS

using Distributions
using DiffEqBase: AbstractJumpProblem
using DiffEqJump: DiffEqBase
include("trajectories/trajectory.jl")
include("EventIter.jl")
include("MetropolisSampler.jl")
include("marginal_strategies/strategies.jl")
include("GaussianSystem.jl")
include("trajectories/distribution.jl")
include("EmpiricalDistribution.jl")
include("model_setup/DrivenJumpProblem.jl")
include("JumpNetwork.jl")
include("write_hdf5.jl")

include("example_systems.jl")

include("ParallelRun.jl")

using .DrivenJumpProblems

export JumpSystem,
    GaussianSystem, marginal_entropy, conditional_entropy,
    TIEstimate, AnnealingEstimate, DirectMCEstimate, SMCEstimate,
    generate_configuration, ComplexSystem,
    mutual_information

end # module


text:
mutable struct ParameterTypes{N}
	nn::Array{Int,1}
	c::NTuple{N,Number}
	d::NTuple{N,Number}
	C::NTuple{N,Number}
	noisestd::NTuple{N,Number}
	vr::NTuple{N,Number}
	cap::NTuple{N,Number}
	a::NTuple{N,Number}
	b::NTuple{N,Number}
	k::NTuple{N,Number}
	vt::NTuple{N,Number}
	dvoltage::NTuple{N,Number}
	synt::NTuple{N,Number}
	quantile::NTuple{N,Number}
	t₋::NTuple{N,Number}
	A₋::NTuple{N,Number}
	tconst::NTuple{N,Number}
	rev::NTuple{N,Number}
	Φ::NTuple{N,Number}
	τ::NTuple{N,Number}
	miniw::NTuple{N,Number}
	σ::NTuple{N,Number}
	init_weight::NTuple{N,Number}
	syn_density::NTuple{N,Number}
	weight_target::NTuple{N,Number}
	normalise_to::NTuple{N,Number}
	winners::NTuple{N,Number}
	da_on::Number
	δt::Number
end

struct MatrixTypes
	activation::NeuronLayers
	rec::NeuronLayers
	spiked::NeuronLayers
	spt::NeuronLayers
	I::NeuronLayers
	ACh::NeuronLayers
	input::NeuronLayers
	synapses::BrainTypes
	weights::BrainTypes
	γ::SynapseLayers
end

ContainerTypes = Union{ParameterTypes, MatrixTypes}

struct NeuroTransmitter
	da::Array
	BA::Array
end

function get_parameters()

	nn = [100, 100, 2]
	c = (-65, -65, -65)
	d = (8, 8, 8)
	C = (100, 4, 100)
	noisestd = (0.05, 0.05, 0.05)
	vr = (-60, -85, -60)
	cap = (100, 4, 100)
	a = (0.3, 0.01, 0.3)#(0.3, 0.01, 0.3)
	b = (-0.2, -0.3, -0.2)
	k = (2, 0.035, 2)#=k = (2, 0.015, 2)=#
	vt = (-40, -25, -40)
	dvoltage = (-60, -85, -60)
	synt = (3, 8, 8)
	quantile = (0.93, 8, 8)
	t₋ = (15, 15, 15)
	A₋ = (-1, -1, -1)
	tconst = (40, 40, 40)
	rev = (0, 0, 0)
	Φ = (0.5, 0.5, 0.5)
	τ = (20, 20, 20)
	miniw = (0.1, 0.1, 0.1)
	σ = (0.05, 0.05, 0.05)
	init_weight = (10, 20, 20)
	syn_density = (0.1, 1, 0.1)
	num_synapses = Int.(round.(nn .* syn_density))
	num_synapses = (4, nn[2], 1)
	weight_target = (1, 5 * nn[2], 1)
	normalise_to = (20000, 500, 200)
	winners = (10, 100, 10)
	da_on = 0.2#0.5#.0009
	δt = 1

	parameters = ParameterTypes(
		nn,
		c,
		d,
		C,
		noisestd,
		vr,
		cap,
		a,
		b,
		k,
		vt,
		dvoltage,
		synt,
		quantile,
		t₋,
		A₋,
		tconst,
		rev,
		Φ,
		τ,
		miniw,
		σ,
		init_weight,
		num_synapses,
		weight_target,
		normalise_to,
		winners,
		da_on,
		δt
		)

	return parameters
end

function initialise_matrices(nn, p)

	activation = NeuronLayers([NeuronLayer(-60.0, i) for i in nn])
	rec = NeuronLayers([NeuronLayer(0.0, i) for i in nn])
	spiked = NeuronLayers([NeuronLayer(false, i) for i in nn])
	spt = NeuronLayers([NeuronLayer(Int(-1), i) for i in nn])
	I = NeuronLayers([NeuronLayer(0.0, i) for i in nn])
	ACh = NeuronLayers([NeuronLayer(0.0, i) for i in nn])
	input = NeuronLayers([NeuronLayer(0.0, i) for i in nn])

	weights = create_synapses(SynapseLayers, nn, p.syn_density, p.init_weight)
	synapses = clone_synapses(create_MBON_synapses(nn, p.syn_density, p.init_weight))

	γ = fill_synapses(SynapseLayers, nn, 0.0)

	return activation, rec, spiked, spt, I, ACh, input, synapses, weights, γ
end

function initialise_matrices(nn, p, weights::SynapseLayers, synapses::BrainTypes)

	activation = NeuronLayers([NeuronLayer(-60.0, i) for i in nn])
	rec = NeuronLayers([NeuronLayer(0.0, i) for i in nn])
	spiked = NeuronLayers([NeuronLayer(false, i) for i in nn])
	spt = NeuronLayers([NeuronLayer(Int(-1), i) for i in nn])
	I = NeuronLayers([NeuronLayer(250.0, i) for i in nn])
	ACh = NeuronLayers([NeuronLayer(0.0, i) for i in nn])
	input = NeuronLayers([NeuronLayer(0.0, i) for i in nn])

	#synapses = synapses#clone_synapses(weights)
	γ = fill_synapses(SynapseLayers, nn, 0.0)

	return activation, rec, spiked, spt, I, ACh, input, synapses, weights, γ
end

function get_parameters(f::typeof(update_activation!), p::ParameterTypes, l::Int=1)

	return (p.vt[l], p.vr[l], p.C[l], p.a[l], p.b[l], p.c[l], p.d[l], p.k[l], p.σ[l], p.δt)
end

function get_parameters(f::typeof(update_weights!), p::ParameterTypes, l::Int=1)

	return (p.tconst[l], p.A₋[l], p.t₋[l], p.miniw[l], p.δt)
end

function get_parameters(f::typeof(update_ACh!), p::ParameterTypes, l::Int=1)

	return (p.Φ[l], p.synt[l], p.δt)
end

function get_matrices(f::typeof(update_activation!), m::MatrixTypes, l::Int=1)

	return (m.activation.layers[l], m.spiked.layers[l], m.spt.layers[l], m.rec.layers[l], m.input.layers[l])
end

function get_matrices(f::typeof(update_weights!), m::MatrixTypes, l::Int=1)

	return (m.weights.layers[l], m.γ.layers[l], m.synapses.layers[l], m.spt.layers[l], m.spt.layers[l+1])
end

function get_matrices(f::typeof(update_ACh!), m::MatrixTypes, l::Int=1)

	return (m.ACh.layers[l], m.spt.layers[l])
end

function get_matrices(f::typeof(calc_input!), m::MatrixTypes, l::Int=1)

	return (m.input.layers[l+1], m.weights.layers[l], m.ACh.layers[l], m.activation.layers[l+1])
end

Base.getproperty(cont::ContainerTypes, str::String) = getfield(cont, Symbol(str)) # allows to get field of struct in MATLAB-like way

Base.setproperty!(cont::ContainerTypes, str::String, new_value) = setfield!(cont, Symbol(str), new_value)

function string(p::ParameterTypes) # not so great b/c it returns an array not a string...

	fields = fieldnames(ParameterTypes)
	out = Array{String,1}(undef,(length(fields)))

	for (i, field) in enumerate(fields)

		out[i] = string(field) * ": " * string(getfield(p, field))

	end
	return out
end


text:
@testset "Test How to Stop II" begin
  ###############################################################################
  #
  # The Stopping structure eases the implementation of algorithms and the
  # stopping criterion.
  # We illustrate here the features of Stopping when the algorithm is used a
  # subStopping.
  #
  ###############################################################################
  #using Test, Stopping

  #Assume we want to solve "pb" starting from "x0" and solving at each step
  #of the algorithm the subproblem "subpb".
  # We can use this additional info to improve the stopping criterion.
  x0 = ones(2)
  pb = nothing
  subpb = nothing
  subsubpb = nothing

  ###############################################################################
  #Initialize a Stopping for the main pb
  main_stop = GenericStopping(pb, x0)
  #We can then, initialize another stopping to the subproblem, and providing
  #the main_stop as a keyword argument:
  sub_stop =
    GenericStopping(subpb, x0, main_stp = main_stop, tol_check = (atol, rtol, opt0) -> atol)
  #Note that by default main_stp is void
  @test main_stop.main_stp == VoidStopping()

  #The only difference appears in the event of a call to stop!, which now also
  #check the time and resources of the main_pb.
  OK = start!(sub_stop)
  @test OK == false #no reason to stop just yet.

  #Assume time is exhausted for the main_stop
  main_stop.meta.start_time = 0.0 #force a timing failure in the main problem
  stop!(sub_stop)

  @test status(sub_stop, list = true) == [:ResourcesOfMainProblemExhausted]
  @test sub_stop.meta.tired == false
  @test sub_stop.meta.main_pb == true

  #The same applies if there is now a third subproblem
  reinit!(main_stop)
  reinit!(sub_stop)
  subsub_stop =
    GenericStopping(subsubpb, x0, main_stp = sub_stop, tol_check = (atol, rtol, opt0) -> atol)
  main_stop.meta.start_time = 0.0 #force a timing failure in the main problem
  stop!(subsub_stop)

  @test status(subsub_stop, list = true) == [:ResourcesOfMainProblemExhausted]
  @test subsub_stop.meta.tired == false
  @test subsub_stop.meta.main_pb == true
  @test status(sub_stop, list = true) == [:ResourcesOfMainProblemExhausted]
  @test sub_stop.meta.tired == false
  @test sub_stop.meta.main_pb == true
  @test status(main_stop, list = true) == [:TimeLimit]
  @test main_stop.meta.tired == true
  @test main_stop.meta.main_pb == false
end


text:
using SchroedingerChess
using Test

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


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

name = "HMMER"
version = v"3.3.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://eddylab.org/software/hmmer/hmmer-$(version).tar.gz",
                  "92fee9b5efe37a5276352d3502775e7c46e9f7a0ee45a331eacb2a0cac713c69")
]

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

export CPPFLAGS="-I${includedir}"
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --with-gsl --enable-threads --enable-sse
make -j${nproc}
make install

install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
#
# Notes:
# - configure demands either sse (i686/x86_64) or vmx (powerpc) vector
#   instructions
# - windows compile fails, can't find syslog.h
# - compile fails on powerpc64le:
#   "HMMER3 Altivec/VMX only supports bigendian platforms: e.g. ppc64 not ppc64le"
platforms = supported_platforms(; exclude = p -> Sys.iswindows(p) || proc_family(p) != "intel")

# The products that we will ensure are always built
products = [
    ExecutableProduct("alimask", :alimask),
    ExecutableProduct("hmmalign", :hmmalign),
    ExecutableProduct("hmmbuild", :hmmbuild),
    ExecutableProduct("hmmconvert", :hmmconvert),
    ExecutableProduct("hmmemit", :hmmemit),
    ExecutableProduct("hmmfetch", :hmmfetch),
    ExecutableProduct("hmmlogo", :hmmlogo),
    ExecutableProduct("hmmpgmd", :hmmpgmd),
    ExecutableProduct("hmmpgmd_shard", :hmmpgmd_shard),
    ExecutableProduct("hmmpress", :hmmpress),
    ExecutableProduct("hmmscan", :hmmscan),
    ExecutableProduct("hmmsearch", :hmmsearch),
    ExecutableProduct("hmmsim", :hmmsim),
    ExecutableProduct("hmmstat", :hmmstat),
    ExecutableProduct("jackhmmer", :jackhmmer),
    ExecutableProduct("makehmmerdb", :makehmmerdb),
    ExecutableProduct("nhmmer", :nhmmer),
    ExecutableProduct("nhmmscan", :nhmmscan),
    ExecutableProduct("phmmer", :phmmer),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GSL_jll", uuid="1b77fbbe-d8ee-58f0-85f9-836ddc23a7a4"))
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function distances(x,y)
    xij = @. x - x'
    yij = @. y - y'
    return xij,yij
end

function periodic_distances!(x,y,Lx,Ly)
    make_periodic!(x,Lx)
    make_periodic!(y,Ly)
end

function make_periodic!(x,L)
    @. x[x>L/2]  -= L
    @. x[x<-L/2] += L
end

"""
    rij = sparse_distances(x,y)

Sparse form of unique lower triangular distance matrix for vectors `x`, `y` specifying points on the plane.
"""
function sparse_distances(x,y)
    xij = @. x - x'
    yij = @. y - y'
    rij = hypot.(xij,yij) |> tril |> sparse
end

"""
    tree = spanning_tree(x,y)

Minimal spanning tree for points specified by vectors `x`, `y`.
"""
function spanning_tree(x,y)
    rij = sparse_distances(x,y)
    # graph representation
    g = SimpleWeightedGraph(length(rij))
    for edge in eachindex(rij)
        n1,n2 = edge[1],edge[2]
        add_edge!(g,n1,n2,rij[n1,n2])
    end
    # get the spanning tree
    tree = prim_mst(g)
end

"""
    vortex_marker(c)

Define vortex marker using sign.
"""
function vortex_marker(c)
    if c==1.0
        return (4,0.6,:green, stroke(1.3, 1.0, :white))
    elseif c==-1.0
        return (4,0.6,:blue, stroke(1.3, 1.0, :white))
    else
    end
end

"""
    plot_vortices!(p,vdata)

Plot vortices defined in `vdata` on figure `p`.
"""
function plot_vortices!(p,vdata)
    for j in 1:size(vdata)[1]
        vx,vy = vdata[j,1],vdata[j,2]
        scatter!(p1,[vx],[vy],marker=vortex_marker(vdata[j,3]),label=false)
    end
    return p1
end

"""
    plot_cluster_vortices!(p,cluster::Cluster)

Plot cluster in figure `p` using `vortex_marker`.
"""
function plot_cluster_vortices!(p,cluster::Cluster)
    vdata = vortex_array(cluster.vortices)
    for j in 1:size(vdata)[1]
        vx,vy = vdata[j,1],vdata[j,2]
        scatter!(p,[vx],[vy],marker=vortex_marker(vdata[j,3]),label=false)
    end
    return p
end

"""
    plot_cluster_tree!(p,cluster::Cluster)

Plot minimal spanning tree on figure `p` for `cluster`.
"""
function plot_cluster_tree!(p,cluster::Cluster)
    vortices = cluster.vortices
    tree = cluster.tree
    va = vortex_array(vortices)
    xi = @view va[:,1]
    yi = @view va[:,2]
    ci = @view va[:,3]

    for i in 1:length(tree)
        edge = tree[i]
        src,dst = edge.src,edge.dst
        x1,y1 = xi[src],yi[src]
        x2,y2 = xi[dst],yi[dst]
        csign = vortices[1].qv
        col = csign == 1.0 ? :green : :blue
        plot!(p,[x1,x2],[y1,y2],label=false,c=col,alpha=0.4)
    end
    return p
end

"""
    plot_custer!(p,c)

Plot the vortices and minimal spanning tree.
"""
function plot_cluster!(p,c)
    plot_cluster_vortices!(p,c)
    plot_cluster_tree!(p,c)
    return p
end


text:
using JSON

# replace get_odata for testing
function CBSOData4.get_odata(url)
    if url == "http://beta-odata4.cbs.nl/CBS/Datasets"
        content = JSON.parse(String(open("meta_10.json") do io read(io) end))
        return content
    elseif url == "http://beta-odata4.cbs.nl/CBS"
        content = JSON.parse(String(open("meta_1.json") do io read(io) end))
        return content
    elseif url == "http://beta-odata4.cbs.nl/CBS/80072ned"
        content = JSON.parse(String(open("meta_2.json") do io read(io) end))
        return content
    elseif url == "http://beta-odata4.cbs.nl/CBS/80072ned/Observations?\$count=true&\$top=1"
        content = JSON.parse(String(open("meta_3.json") do io read(io) end))
        return content
    elseif url == "http://beta-odata4.cbs.nl/CBS/80072ned/Observations?"
        content = JSON.parse(String(open("meta_4.json") do io read(io) end))
        return content
    elseif url == "http://beta-odata4.cbs.nl/CBS/80072ned/MeasureCodes"
        content = JSON.parse(String(open("meta_5.json") do io read(io) end))
        return content
    elseif url == "http://beta-odata4.cbs.nl/CBS/80072ned/Dimensions"
        content = JSON.parse(String(open("meta_6.json") do io read(io) end))
        return content
    elseif url == "http://beta-odata4.cbs.nl/CBS/80072ned/Properties"
        content = JSON.parse(String(open("meta_7.json") do io read(io) end))
        return content
    elseif url == "http://beta-odata4.cbs.nl/CBS/80072ned/PeriodenGroups"
        content = JSON.parse(String(open("meta_8.json") do io read(io) end))
        return content
    elseif url == "http://beta-odata4.cbs.nl/CBS/80072ned/PeriodenCodes"
        content = JSON.parse(String(open("meta_9.json") do io read(io) end))
        return content
    elseif url == "http://beta-odata4.cbs.nl/CBS/80072ned/BedrijfskenmerkenSBI2008Groups"
        content = JSON.parse(String(open("meta_10.json") do io read(io) end))
        return content
    elseif url == "http://beta-odata4.cbs.nl/CBS/80072ned/BedrijfskenmerkenSBI2008Codes"
        content = JSON.parse(String(open("meta_11.json") do io read(io) end))
        return content
    else
        println("Fout:" * url)
        return Dict("value"=>nothing)
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ c4a4bfe6-79fc-11eb-389d-fb0dda01299f
begin
	using PlutoUI
	using NativeSVG
end

# ╔═╡ 6e3f9f5e-79fc-11eb-22ec-47b9055afb05
md"""# Case Study: Word Play

This chapter presents a second case study, which involves solving word puzzles by searching for words that have certain properties. For example, we’ll find the longest palindromes in English and search for words whose letters appear in alphabetical order. And I will present another program development plan: reduction to a previously solved problem."""

# ╔═╡ dc744022-7a03-11eb-015f-47cde5608878
md"""## Reading Word Lists

For the exercises in this chapter we need a list of English words. There are lots of word lists available on the web, but the one most suitable for our purpose is one of the word lists collected and contributed to the public domain by Grady Ward as part of the [Moby lexicon](http://www.gutenberg.org/ebooks/3201) project. It is a list of 113,809 official crosswords; that is, words that are considered valid in crossword puzzles and other word games. In the Moby collection, the filename is *113809of.fic*; you can download a copy, with the simpler name *words.txt*, from these lectures’ GitHub repository.

This file is in plain text, so you can open it with a text editor, but you can also read it from Julia. The built-in function open takes a name of the file as a parameter and returns a file stream you can use to read the file:

```julia
julia> fin = open("words.txt")
IOStream(<file words.txt>)
```

`fin` is a file stream used for input. When it is no longer needed, it has to be closed with `close(fin)`.

Julia provides several function for reading, including `readline`, which reads characters from the file until it gets to a *NEWLINE* and returns the result as a string:

```julia
julia> readline(fin)
"aa"
```

The first word in this particular list is “aa”, which is a kind of lava.

The file stream keeps track of where it is in the file, so if you call `readlin`e again, you get the next word:

```julia
julia> readline(fin)
"aah"
```

The next word is “aah”, which is a perfectly legitimate word, so stop looking at me like that.

You can also use a file as part of a `for` loop. This program reads *words.txt* and prints each word, one per line:

```julia
for line in eachline("words.txt")
	println(line)
end
```
"""

# ╔═╡ 1b1cb464-7a05-11eb-1069-a5d56ef55821
md"""## Exercises

#### Exercise 9-1

Write a program that reads *words.txt* and prints only the words with more than 20 characters (not counting whitespace).

#### Exercise 9-2

In 1939 Ernest Vincent Wright published a 50,000-word novel called Gadsby (Wetzel Publishing) that does not contain the letter *e*. Since *e* is the most common letter in English, that’s not easy to do.

In fact, it is difficult to construct a solitary thought without using that most common symbol. It is slow going at first, but with caution and hours of training you can gradually gain facility.

All right, I’ll stop now.

Write a function called `hasno_e` that returns true if the given word doesn’t have the
letter *e* in it.

Modify your program from the previous exercise to print only the words that have no *e* and compute the percentage of the words in the list that have no *e*.

#### Exercise 9-3

Write a function named `avoids` that takes a word and a string of forbidden letters, and that returns true if the word doesn’t use any of the forbidden letters.

Modify your program to prompt the user to enter a string of forbidden letters and then print the number of words that don’t contain any of them. Can you find a combination of five forbidden letters that excludes the smallest number of words?

#### Exercise 9-4

Write a function named `usesonly` that takes a word and a string of letters, and that returns true if the word contains only letters in the list. Can you make a sentence using only the letters *acefhlo*? Other than *"Hoe alfalfa"*?

#### Exercise 9-5

Write a function named `usesall` that takes a word and a string of required letters, and that returns true if the word uses all the required letters at least once. How many words are there that use all the vowels *aeiou*? How about *aeiouy*?

#### Exercise 9-6

Write a function called `isabecedarian` that returns `tru`e if the letters in a word appear in alphabetical order (double letters are okay). How many abecedarian words are there?
"""

# ╔═╡ 3ec4f68e-7a06-11eb-0942-0399471b3818
md"""## Search

All of the exercises in the previous section have something in common; they can be solved with the search pattern. The simplest example is:

```julia
function hasno_e(word)
	for letter in word
		if letter == 'e'
			return false
		end
	end
	true
end
```

The `for` loop traverses the characters in `word`. If we find the letter *e*, we can immediately return `false`; otherwise, we have to go to the next letter. If we exit the loop normally, that means we didn’t find an *e*, so we return `true`.

You could write this function more concisely using the `∉` (**`\notin TAB`**) operator, but I started with this version because it demonstrates the logic of the search pattern.

`avoids` is a more general version of `hasno_e`, but it has the same structure:

```julia
function avoids(word, forbidden)
	for letter in word
		if letter ∈ forbidden
			return false
		end
	end
	true
end
```

We can return `false` as soon as we find a forbidden letter; if we get to the end of the loop, we return `true`.

`usesonly` is similar except that the sense of the condition is reversed:

```julia
function usesonly(word, available)
	for letter in word
		if letter ∉ available
			return false
		end
	end
	true
end
```

Instead of an array of forbidden letters, we have an array of available letters. If we find a letter in word that is not in available, we can return `false`.

`usesall` is similar except that we reverse the role of the word and the string of letters:

```julia
function usesall(word, required)
	for letter in required
		if letter ∉ word
			return false
		end
	end
	true
end
```

Instead of traversing the letters in `word`, the loop traverses the required letters. If any of the required letters do not appear in `word`, we can return `false`.

If you were really thinking like a computer scientist, you would have recognized that `usesall` was an instance of a previously solved problem, and you would have written:

```julia
usesall(word, required) =  usesonly(required, word)
```

This is an example of a program development plan called *reduction to a previously solved problem*, which means that you recognize the problem you are working on as an instance of a solved problem and apply an existing solution.
"""

# ╔═╡ 07b18184-7a07-11eb-165a-ff024559c0a5
md"""## Looping with Indices

I wrote the functions in the previous section with for loops because I only needed the characters in the strings; I didn’t have to do anything with the indices.

For `isabecedarian` we have to compare adjacent letters, which is a little tricky with a `for` loop:

```julia
function isabecedarian(word)
	i = firstindex(word)
	previous = word[i]
	j = nextind(word, i)
	for c in word[j:end]
		if c < previous
			return false
		end
		previous = c
	end
	true
end
```

An alternative is to use recursion:

```julia
function isabecedarian(word)
	if length(word) <= 1
		return true
	end
	i = firstindex(word)
	j = nextind(word, i)
	if word[i] > word[j]
		return false
	end
	isabecedarian(word[j:end])
end
```

Another option is to use a `while` loop:

```julia
function isabecedarian(word)
	i = firstindex(word)
	j = nextind(word, 1)
	while j <= sizeof(word)
		if word[j] < word[i]
			return false
		end
		i = j
		j = nextind(word, i)
	end
	true
end
```

The loop starts at `i=1` and `j=nextind(word, 1)` and ends when `j>sizeof(word)`. Each time through the loop, it compares the ith character (which you can think of as the current character) to the `j`th character (which you can think of as the next).

If the next character is less than (alphabetically before) the current one, then we have discovered a break in the abecedarian trend, and we return `false`.

If we get to the end of the loop without finding a fault, then the word passes the test. To convince yourself that the loop ends correctly, consider an example like `"flossy"`.

Here is a version of `ispalindrome` that uses two indices; one starts at the beginning and goes up, and the other starts at the end and goes down:

```julia
function ispalindrome(word)
	i = firstindex(word)
	j = lastindex(word)
	while i < j
		if word[i] != word[j]
			return false
		end
		i = nextind(word, i)
		j = prevind(word, j)
	end
	true
end
```

Or we could reduce to a previously solved problem and write function

```julia
ispalindrome(word) = isreverse(word, word)
```

using `isreverse` from “Debugging”.
"""

# ╔═╡ d373aa0e-7a07-11eb-14d0-c7e1a531b225
md"""## Debugging

Testing programs is hard. The functions in this chapter are relatively easy to test because you can check the results by hand. Even so, it is somewhere between difficult and impossible to choose a set of words that tests for all possible errors.

Taking `hasno_e` as an example, there are two obvious cases to check: words that have an *e* should return `false`, and words that don’t should return `true`. You should have no trouble coming up with one of each.

Within each case, there are some less obvious subcases. Among the words that have an *e*, you should test words with an *e* at the beginning, the end, and somewhere in the middle. You should test long words, short words, and very short words, like the empty string. The empty string is an example of a *special case*, which is one of the nonobvious cases where errors often lurk.

In addition to the test cases you generate, you can also test your program with a word list like *words.txt*. By scanning the output, you might be able to catch errors, but be careful: you might catch one kind of error (words that should not be included, but are) and not another (words that should be included, but aren’t).

In general, testing can help you find bugs, but it is not easy to generate a good set of test cases, and even if you do, you can’t be sure your program is correct. According to a legendary computer scientist:

> Program testing can be used to show the presence of bugs, but never to show their absence!
> —Edsger W. Dijkstra
"""

# ╔═╡ 71ae6498-7a08-11eb-025e-03eb7a5c2d8b
md"""## Glossary

*file stream*:
A value that represents an open file.

*reduction to a previously solved problem*:
A way of solving a problem by expressing it as an instance of a previously solved problem.

*special case*:
A test case that is atypical or nonobvious (and less likely to be handled correctly).
"""

# ╔═╡ 9e542172-7a08-11eb-2d03-75c8a0ff571a
md"""## Exercise

#### Exercise 9-7

This question is based on a Puzzler that was broadcast on the radio program [Car Talk](https://www.cartalk.com/):

> Give me a word with three consecutive double letters. I’ll give you a couple of words that almost qualify, but don’t. For example, the word committee, c-o-m-m-i-t-t-e-e. It would be great except for the i that sneaks in there. Or Mississippi—M-i-s-s-i-s-s-i-p- p-i. If you could take out those i’s it would work. But there is a word that has three consecutive pairs of letters and to the best of my knowledge this may be the only word. Of course there are probably 500 more but I can only think of one. What is the word?

Write a program to find it.
"""

# ╔═╡ 3121945a-7a09-11eb-3274-553912ceed46
md"""#### Exercise 9-8

Here’s another Car Talk Puzzler:

> I was driving on the highway the other day recently and I happened to notice my odometer. Like most odometers nowadays, it shows six digits, in whole miles only—no tenths of a mile. So, if my car had 300,000 miles, for example, I’d see 3-0-0-0-0-0. ...
>
> Now, what I saw that day was very interesting. I noticed that the last 4 digits were pal‐ indromic, that is they read the same forwards as backwards. For example, “5-4-4-5” is a palindrome. So my odometer could have read 3-1-5-4-4-5 ... .
>
> One mile later, the last 5 numbers were palindromic. For example, it could have read 3-6-5-4-5-6.
>
> One mile after that, the middle 4 out of 6 numbers were palindromic. ... And you ready for this? One mile later, all 6 were palindromic! ...
>
> The question is, what did [I] see on the odometer when [I] first looked?

Write a Julia program that tests all the six-digit numbers and prints any numbers that satisfy these requirements.
"""

# ╔═╡ 6560b396-7a09-11eb-3a40-cf56131c10f6
md"""#### Exercise 9-9

Here’s a third Car Talk Puzzler that you can solve with a search:

> Recently I had a visit with my mom and we realized that the two digits that make up my age when reversed result in her age. For example, if she’s 73, I’m 37. We wondered how often this has happened over the years but we got sidetracked with other topics and we never came up with an answer.
>
> When I got home I figured out that the digits of our ages have been reversible six times so far. I also figured out that if we’re lucky it would happen again in a few years, and if we’re really lucky it would happen one more time after that. In other words, it would have happened 8 times over all. So the question is, how old am I now?

Write a Julia program that searches for solutions to this Puzzler.

!!! tip
    You might find the function `lpad` useful.
"""

# ╔═╡ Cell order:
# ╟─c4a4bfe6-79fc-11eb-389d-fb0dda01299f
# ╟─6e3f9f5e-79fc-11eb-22ec-47b9055afb05
# ╟─dc744022-7a03-11eb-015f-47cde5608878
# ╟─1b1cb464-7a05-11eb-1069-a5d56ef55821
# ╟─3ec4f68e-7a06-11eb-0942-0399471b3818
# ╟─07b18184-7a07-11eb-165a-ff024559c0a5
# ╟─d373aa0e-7a07-11eb-14d0-c7e1a531b225
# ╟─71ae6498-7a08-11eb-025e-03eb7a5c2d8b
# ╟─9e542172-7a08-11eb-2d03-75c8a0ff571a
# ╟─3121945a-7a09-11eb-3274-553912ceed46
# ╟─6560b396-7a09-11eb-3a40-cf56131c10f6


text:
using OCToolbox

using LinearAlgebra
using QuantumInformation

i2 = Matrix{Complex{Float64}}(I, 2, 2)
H_drift = 0*sz
H_ctrl = [π*sx, π * sy]

# we define our initial and final states
Ψ = [1+0.0im, 0.0]
ρ = [0+0.0im, 1.0]

# we set up a functional
function fn(controls)
    controls = complex.(real.(controls))
    #U = stack_props(controls, H_drift)
    U = pw_full_evolution(H_drift, H_ctrl, controls, Δt, i2)
    C2(ρ, U * Ψ)
end

K = 2
N = 10
T = 2
Δt = T/N

control_guess = rand(K, N).*0.001
fn(control_guess)

o = GRAPE(fn, control_guess, K, N)

fn(o.minimizer)