text:
@testset "Transfer Entropy" begin
    @test_throws MethodError TransferEntropy(0.5, 1.0; k=2)
    @test_throws TypeError TransferEntropy(2, 2; k=2.0)
    @test_throws ArgumentError TransferEntropy(2, 2; k=0)
    for b in -2:1
        @test_throws ArgumentError TransferEntropy(b, 2; k=2)
        @test_throws ArgumentError TransferEntropy(2, b; k=2)
    end

    @test_throws ArgumentError TransferEntropy(Int[], Int[]; k=2)
    @test_throws ArgumentError TransferEntropy([1,2,2], Int[]; k=2)
    @test_throws ArgumentError TransferEntropy(Int[], [1,2,3]; k=2)
    @test_throws ArgumentError TransferEntropy([0,1,2], [1,2,2]; k=2)
    @test_throws ArgumentError TransferEntropy([1,1,2], [1,0,2]; k=2)

    let te = TransferEntropy([1,2,1,2,1,2], [1,1,2,2,1,1]; k=2)
        clear!(te)
        @test te.states == zeros(Int, 16)
        @test te.histories == zeros(Int, 4)
        @test te.sources == zeros(Int, 8)
        @test te.predicates == zeros(Int, 8)
        @test te.N == zero(Int)
    end

    @test_throws ArgumentError observe!(TransferEntropy(2,2; k=2), [1,1], [1,2,1,1,2])
    @test_throws ArgumentError observe!(TransferEntropy(2,2; k=2), [1,2], [2,2])

    let as = [2,2,2,1,1], bs = [2,2,1,1,2]
        @test transferentropy(as, as; k=2) ≈ 0.0 atol=1e-6
        @test transferentropy(bs, bs; k=2) ≈ 0.0 atol=1e-6
        @test transferentropy(as, bs; k=2) ≈ 0.0 atol=1e-6
        @test transferentropy(bs, as; k=2) ≈ 2/3 atol=1e-6
    end

    let as = [1,1,2,2,2,1,1,1,1,2]
        bs = [2,2,1,1,1,1,1,1,2,2]
        @test transferentropy(as, as; k=2) ≈ 0.0 atol=1e-6
        @test transferentropy(bs, bs; k=2) ≈ 0.0 atol=1e-6
        @test transferentropy(as, bs; k=2) ≈ 0.106844 atol=1e-6
        @test transferentropy(bs, as; k=2) ≈ 1/2 atol=1e-6
    end

    let as = [1,2,1,2,1,1,2,2,1,1]
        bs = [1,1,2,1,2,2,2,1,2,2]
        @test transferentropy(as, as; k=2) ≈ 0.0 atol=1e-6
        @test transferentropy(bs, bs; k=2) ≈ 0.0 atol=1e-6
        @test transferentropy(as, bs; k=2) ≈ 1/4 atol=1e-6
        @test transferentropy(bs, as; k=2) ≈ 0.344361 atol=1e-6
    end

    let as = [1,1,1,2,2], bs = [1,1,2,2,1]
        @test transferentropy(as, as; k=2) ≈ 0.0 atol=1e-6
        @test transferentropy(bs, bs; k=2) ≈ 0.0 atol=1e-6
        @test transferentropy(as, bs; k=2) ≈ 0.0 atol=1e-6
        @test transferentropy(bs, as; k=2) ≈ 2/3 atol=1e-6
    end

    let as = [2,2,1,1,1,2,2,2,2,1]
        bs = [1,1,2,2,2,2,2,2,1,1]
        @test transferentropy(as, as; k=2) ≈ 0.0 atol=1e-6
        @test transferentropy(bs, bs; k=2) ≈ 0.0 atol=1e-6
        @test transferentropy(as, bs; k=2) ≈ 0.106844 atol=1e-6
        @test transferentropy(bs, as; k=2) ≈ 1/2 atol=1e-6
    end

    let as = [2,1,2,1,2,2,1,1,2,2]
        bs = [2,2,1,2,1,1,1,2,1,1]
        @test transferentropy(as, as; k=2) ≈ 0.0 atol=1e-6
        @test transferentropy(bs, bs; k=2) ≈ 0.0 atol=1e-6
        @test transferentropy(as, bs; k=2) ≈ 1/4 atol=1e-6
        @test transferentropy(bs, as; k=2) ≈ 0.344361 atol=1e-6
    end

    let as = [2,1,2,1,2,2,1,1,2,2]
        bs = [2,2,1,2,1,1,1,2,1,1]
        te = TransferEntropy(2, 2; k=2)
        @test transferentropy!(te, as, bs) ≈ 1/4 atol=1e-6
    end

    let as = [4,4,4,3,2,1,1,1,2]
        bs = [3,3,4,4,1,2,2,1,3]
        @test transferentropy(as, bs; k=2) ≈ 0.000000 atol=1e-6
        @test transferentropy(bs, as; k=2) ≈ 0.571429 atol=1e-6
    end

    let as = [4,2,1,1,3,2,1,3,3]
        bs = [3,3,2,2,2,2,3,1,4]
        @test transferentropy(as, bs; k=2) ≈ 0.393555 atol=1e-6
        @test transferentropy(bs, as; k=2) ≈ 0.571429 atol=1e-6
    end

    let as = [2 2 2 1 2 1 1 1 2;
              2 2 2 2 1 1 1 1 1]
        bs = [1 1 2 2 1 2 2 1 1;
              2 2 2 2 1 1 1 1 2]
        @test transferentropy(as, bs; k=2) ≈ 0.000000 atol=1e-6
        @test transferentropy(bs, as; k=2) ≈ 0.571429 atol=1e-6
    end

    let as = [2 2 1 1 1 2 1 1 1;
              2 1 1 1 2 1 1 2 2]
        bs = [1 1 2 2 2 2 1 1 2;
              2 2 1 1 1 1 2 1 2]
        @test transferentropy(as, bs; k=2) ≈ 0.393555 atol=1e-6
        @test transferentropy(bs, as; k=2) ≈ 0.571429 atol=1e-6
    end

    let as = zeros(Int, 2, 9, 2)
        as[:,:,1] = [2 2 2 1 2 1 1 1 2;
                     2 2 2 2 1 1 1 1 1]
        as[:,:,2] = [2 2 1 1 1 2 1 1 1;
                     2 1 1 1 2 1 1 2 2]

        bs = zeros(Int, 2, 9, 2)
        bs[:,:,1] = [1 1 2 2 1 2 2 1 1;
                     2 2 2 2 1 1 1 1 2]
        bs[:,:,2] = [1 1 2 2 2 2 1 1 2;
                     2 2 1 1 1 1 2 1 2]

        @test transferentropy(as, bs; k=2) ≈ 0.428571 atol=1e-6
        @test transferentropy(bs, as; k=2) ≈ 0.679270 atol=1e-6
    end

    let as = [2 2 1 1 1 2 2 2 2 1;
              2 1 2 1 2 2 1 1 2 2]
        bs = [1 1 2 2 2 2 2 2 1 1;
              2 2 1 2 1 1 1 2 1 1]
        te = TransferEntropy(2, 2; k=2)
        @test transferentropy!(te, as[1,:], bs[1,:]) ≈ 0.106844 atol=1e-6
        @test transferentropy!(te, as[2,:], bs[2,:]) ≈ 0.047181 atol=1e-6
    end
end


text:
import GeometryBasics
import Observables
import AbstractPlotting

export Sprite
export getsurface
export buildsurface
export updatesurface

"""
    Represents a shape in a scene.
"""
abstract type Sprite end

"""
    getsurface(observables, segments)

Get the array of points of a 2-surface with the given tuple of `observables` and the number
of `segments`.
"""
function getsurface(observable::Tuple{Observables.Observable{Array{Float64,2}},
                                      Observables.Observable{Array{Float64,2}},
                                      Observables.Observable{Array{Float64,2}}},
                    segments::Int)
    y₁, y₂, y₃ = map(x -> Observables.to_value(x), observable)
    value = Array{ℝ³,2}(undef, segments, segments)
    for i in 1:segments
        for j in 1:segments
            index = (j - 1) * segments + i
            value[i, j] = ℝ³(y₁[index], y₂[index], y₃[index])
        end
    end
    value
end

"""
    getsurface(observables, segments1, segments2)

Get the array of points of a 2-surface with the given tuple of `observables`, and the number
of `segments1` and `segments2`.
"""
function getsurface(observable::Tuple{Observables.Observable{Array{Float64,2}},
                                      Observables.Observable{Array{Float64,2}},
                                      Observables.Observable{Array{Float64,2}}},
                    segments1::Int,
                    segments2::Int)
    y₁, y₂, y₃ = map(x -> Observables.to_value(x), observable)
    value = Array{ℝ³,2}(undef, segments1, segments2)
    for i in 1:segments1
        for j in 1:segments2
            index = (j - 1) * segments1 + i
            value[i, j] = ℝ³(y₁[index], y₂[index], y₃[index])
        end
    end
    value
end

"""
    buildsurface(scene, value, color)

Build a surface with the given `scene`, `value` and `color`. An optional argument is
`transparency`.
"""
function buildsurface(scene::AbstractPlotting.Scene,
                      value::Array{ℝ³,2},
                      color::Array{AbstractPlotting.ColorTypes.RGBA{Float32},2};
                      transparency::Bool = false)
    x = Observables.Observable(map(x -> vec(x)[1] , value))
    y = Observables.Observable(map(x -> vec(x)[2] , value))
    z = Observables.Observable(map(x -> vec(x)[3] , value))
    AbstractPlotting.surface!(scene, x, y, z, color = color, transparency = transparency)
    x, y, z
end

"""
    buildsurface(scene, value, color)

Build a surface with the given `scene`, `value` and `color`. An optional argument is
`transparency`.
"""
function buildsurface(scene::AbstractPlotting.Scene,
                      value::Array{ℝ³,2},
                      color::Observables.Observable{Array{AbstractPlotting.ColorTypes.RGBA{
                                                          Float32},2}};
                      transparency::Bool = false)
    x = Observables.Observable(map(x -> vec(x)[1] , value))
    y = Observables.Observable(map(x -> vec(x)[2] , value))
    z = Observables.Observable(map(x -> vec(x)[3] , value))
    AbstractPlotting.surface!(scene, x, y, z, color = color, transparency = transparency)
    x, y, z
end

"""
    buildsurface(scene, value, color)

Build a surface with the given `scene`, `value` and `color`. An optional argument is
`transparency`.
"""
function buildsurface(scene::AbstractPlotting.Scene,
                      value::Array{ℝ³,2},
                      color::Any;
                      transparency::Bool = false)
    x = Observables.Observable(map(x -> vec(x)[1] , value))
    y = Observables.Observable(map(x -> vec(x)[2] , value))
    z = Observables.Observable(map(x -> vec(x)[3] , value))
    AbstractPlotting.surface!(scene, x, y, z, color = color, transparency = transparency)
    x, y, z
end

"""
    updatesurface(value, observable)

Update a surface with the given `value` and `observable`.
"""
function updatesurface(value::Array{ℝ³,2},
                       observable::Tuple{Observables.Observable{Array{Float64,2}},
                                         Observables.Observable{Array{Float64,2}},
                                         Observables.Observable{Array{Float64,2}}})
    x, y, z = observable
    x[] = map(x -> vec(x)[1], value)
    y[] = map(x -> vec(x)[2], value)
    z[] = map(x -> vec(x)[3], value)
end

ext: jl


text:
##
using Random, MAT, Statistics, Printf

##
function forward(m,x)
    d=m[1].+m[2].*x + m[3]*x.^2
    return d
end

function log_likelihood(d,d_obs,d_std)
    logL=-0.5*sum(((d_obs.-d).^2)./(d_std.^2))
end

function sample_prior()
    m1 = randn(1).*sqrt(10)
    m2 = (rand(1).-0.5).*20
    m3 = (rand(1).-0.5).*0.2
    m=[m1[1],m2[1],m3[1]]
    return m
end

## LOAD SOME DATA
#vars = matread("linefit_nd05.mat")
vars = matread("linefit_nd35.mat")
x_obs = vars["x_obs"][:]
d_obs = vars["d_obs"][:]
d_std = vars["d_std"][:]
nd=length(d_obs)
nm=3

## Rejection sampler
N_ite=2000000
N_acc=0
logL_max = -250

m_post=zeros((nm,N_ite))
#println("Start rejection sampling")
t_start=time()
for i = 1:N_ite
    m = sample_prior()
    d = forward(m,x_obs)
    logL=log_likelihood(d,d_obs,d_std)

    if logL_max<logL
        global logL_max = logL*.99
        #println("logL_max = ",logL_max, ", i=",i)
        global N_acc = 0
    end

    Pacc = exp(logL-logL_max)
    if rand(1)[1]<Pacc
        N_acc = N_acc+1
        m_post[:,N_acc]=m
     end
end
t_stop=time()
#println("End rejection sampling")
t_elapsed  = t_stop-t_start

@printf("%34s%6s: t=%6.2fs, N_ite=%8d, %8d iterations/s\n", "JULIA", VERSION, t_elapsed, N_ite, ceil(N_ite/t_elapsed))

## Posterior statistics
showPostStat=0
if showPostStat==1
    m_post = m_post[:,1:N_acc]
    m_post_mean = mean(m_post, dims=2)
    m_post_std = std(m_post, dims=2)
    @printf("m1_mean=%4.2f, m2_mean=%4.2f, m3_mean=%4.3f\n", m_post_mean[1], m_post_mean[2], m_post_mean[3])
    @printf("m1_std=%4.2f, m2_std=%4.2f, m3_std=%4.3f\n", m_post_std[1], m_post_std[2], m_post_std[3])
end


text:
using LinearAlgebra, PyPlot, Random, FFTW, Distributions, CircularArrays, Printf, Statistics

close(:all)

ci, co, N, ps, b = 16, 16, 256, 64, 50
plot_I = false

eb = randn(Float32, co*N, ps)

chunk(arr, n) = [arr[i:min(i + n - 1, end)] for i in 1:n:length(arr)]
simil(x, y) = 100*dot(x[:], y[:])/(norm(x)*norm(y))
error(x, y) = norm(x - y)/(norm(x)+norm(y))

function draw_e(ps::Integer, co::Integer, N::Integer)
    n = ps ÷ co
    if n < 8
        n = 8
        inds = chunk(randperm(n*co) .% ps .+ 1, n)
        overlap = 1 #n / (ps ÷ co)
    else
        inds = chunk(randperm(ps), n)
        overlap = 1
    end
    e = zeros(Float32, N*co, ps)
    for i=1:co
        iloc = inds[i]
        ni = length(iloc)
        e[(i-1)*N+1:i*N, iloc] .= Float32(sqrt(ps/(overlap*ni)))*randn(Float32, N, ni)
    end

    return e
end

eortho = draw_e(ps, co, N)

if plot_I
    plot_e_scale = ps/50
    figure()
    subplot(221)
    imshow(eb, vmin=-.1, vmax=.1, cmap="seismic", aspect="auto")
    title("Z, Random")
    subplot(222)
    imshow(eortho, vmin=-.1, vmax=.1, cmap="seismic", aspect="auto")
    title("Z, Random Orthogonalized")
    subplot(223)
    imshow(eb*eb', vmin=-plot_e_scale, vmax=plot_e_scale, cmap="seismic", aspect="auto")
    title(L"$ZZ^T$ Random")
    subplot(224)
    imshow(eortho*eortho', vmin=-plot_e_scale, vmax=plot_e_scale, cmap="seismic", aspect="auto")
    title(L"$ZZ^T$ Random Orthogonalized")
    tight_layout()
    savefig("./figures/zortho.png", bbox_inches="tight")
end

####### PROBE 3D Tensore with it

# Simple X Y^T struct to avoid blowing up memory
struct outer_LR
    L
    R
end

Base.:*(A::outer_LR, x) = A.L*(A.R'*x)
Base.:*(x, A::outer_LR) = (x*A.L)*A.R'
block(A::outer_LR, i, j) = view(A.L, (i-1)*N+1:i*N, :) * view(A.R, (j-1)*N+1:j*N, :)'

ch_scales = rand([.001f0, .01f0, .1f0, 1f0, 100f0, 100f0, 1000f0], ci) .* rand([-1f0, 1f0], ci)
a = vcat([ch_scales[i] * randn(Float32, N, b) for i=1:ci]...)

A = outer_LR(a, max.(0f0, 10 .* randn(Float32, N*co, b)))

function probe_trace(ci::Integer, co::Integer, N::Integer, A::outer_LR, e::Matrix{T}) where T
    # A * z
    xyz = A * e
    trl = zeros(T, ci, co)
    for cii = 1:ci
        for coi = 1:co
            eloc = view(e, (coi - 1)*N+1:coi*N, :)'
            xyxloc = view(xyz, (cii - 1)*N+1:cii*N, :)
            trl[cii, coi] = LinearAlgebra.tr(eloc * xyxloc)/ps
        end
    end
    return trl
end
true_tr = zeros(Float32, ci, co)
[true_tr[i, j] = LinearAlgebra.tr(block(A, i, j)) for i=1:ci for j=1:co]

# Single sample
trr = probe_trace(ci, co, N, A, eb)
tro = probe_trace(ci, co, N, A, eortho)

function plot_tr(trr, tro, true_tr)
    iplot = rand(1:ci, 4)
    figure()
    for (ip, i)=enumerate(iplot)
        subplot(2, 2, ip)
        title("Input channel $(i)")
        plot(trr[i, :], "--*r", label="random", linewidth=1, markersize=5)
        plot(tro[i, :], "--ob", label="ortho random", linewidth=1, markersize=5, markerfacecolor="none")
        plot(true_tr[i, :], "--xk", label="true", linewidth=1, markersize=5)
    end
    legend()
end

plot_tr(trr, tro, true_tr)

# Mean
np = 5000

convergence = zeros(Float32, np, 2)
similarity = zeros(Float32, np, 2)
convergence[1, 1] = error(true_tr, trr)
convergence[1, 2] = error(true_tr, tro)

similarity[1, 1] = simil(true_tr, trr)
similarity[1, 2] = simil(true_tr, tro)

i=1

@printf("%d/%d : Er = %2.2e, Sr = %2.2f, Eo = %2.2e, So = %2.2f \n",i, np,
        convergence[i, 1], similarity[i, 1], convergence[i, 2], similarity[i, 2])

for i=2:np
    local eeb = randn(Float32, co*N, ps)
    local eeortho = draw_e(ps, co, N)

    global trr = ((i-1)*trr + probe_trace(ci, co, N, A, eeb))/i
    global tro = ((i-1)*tro + probe_trace(ci, co, N, A, eeortho))/i
    convergence[i, 1] = error(true_tr, trr)
    convergence[i, 2] = error(true_tr, tro)
    similarity[i, 1] = simil(true_tr, trr)
    similarity[i, 2] = simil(true_tr, tro)
    @printf("%d/%d : Er = %2.2e, Sr = %2.2f, Eo = %2.2e, So = %2.2f \n",i, np,
            convergence[i, 1], similarity[i, 1], convergence[i, 2], similarity[i, 2])
end

ni(x) = x / norm(x, Inf)

plot_tr(trr, tro, true_tr)

figure()
plot((true_tr[:] .- trr[:])./true_tr[:], label="Relative error random")
plot((true_tr[:] .- tro[:])./true_tr[:], label=" Relative error ortho random")
legend()

figure()
subplot(211)
loglog(convergence[1:np, 1], label="Error random")
loglog(convergence[1:np, 2], label="Error random ortho")
legend()
subplot(212)
plot(similarity[1:np, 1], label="Similarity random")
plot(similarity[1:np, 2], label="Similarity random ortho")
yticks([i*10 for i=0:2:10], labels=["$(i*10)%" for i=0:2:10])
legend()

# independent
# r x ci x batch
# all
# r x batch

max_stars_repo_name: sheevy/XConv

max_issues_repo_name: sheevy/XConv

max_forks_repo_name: sheevy/XConv


text:
module Downloads

using ArgTools

include("Curl/Curl.jl")
using .Curl

## Base download API ##

export download

struct Downloader
    multi::Multi
    Downloader() = new(Multi())
end

const DEFAULT_DOWNLOADER = Ref{Union{Downloader, Nothing}}(nothing)

function default_downloader()::Downloader
    DEFAULT_DOWNLOADER[] isa Downloader && return DEFAULT_DOWNLOADER[]
    DEFAULT_DOWNLOADER[] = Downloader()
end

const Headers = Union{AbstractVector, AbstractDict}

"""
    download(url, [ output = tempfile() ]; [ headers ]) -> output

        url     :: AbstractString
        output  :: Union{AbstractString, AbstractCmd, IO}
        headers :: Union{AbstractVector, AbstractDict}

Download a file from the given url, saving it to `output` or if not specified,
a temporary path. The `output` can also be an `IO` handle, in which case the
body of the response is streamed to that handle and the handle is returned. If
`output` is a command, the command is run and output is sent to it on stdin.

If the `headers` keyword argument is provided, it must be a vector or dictionary
whose elements are all pairs of strings. These pairs are passed as headers when
downloading URLs with protocols that supports them, such as HTTP/S.
"""
function download(
    url::AbstractString,
    output::Union{ArgWrite, Nothing} = nothing;
    headers::Headers = Pair{String,String}[],
    downloader::Downloader = default_downloader(),
)
    yield() # prevents deadlocks, shouldn't be necessary
    arg_write(output) do io
        easy = Easy()
        set_url(easy, url)
        for hdr in headers
            hdr isa Pair{<:AbstractString, <:Union{AbstractString, Nothing}} ||
                throw(ArgumentError("invalid header: $(repr(hdr))"))
            add_header(easy, hdr)
        end
        add_handle(downloader.multi, easy)
        for buf in easy.buffers
            write(io, buf)
        end
        remove_handle(downloader.multi, easy)
        status = get_response_code(easy)
        status == 200 && return
        if easy.code == Curl.CURLE_OK
            message = get_response_headers(easy)[1]
        else
            message = GC.@preserve easy unsafe_string(pointer(easy.errbuf))
        end
        error(message)
    end
end

## experimental request API ##

export request, Request, Response

struct Request
    io::IO
    url::String
    headers::Vector{Pair{String,String}}
end

struct Response
    url::String
    status::Int
    response::String
    headers::Vector{Pair{String,String}}
end

function request(req::Request, multi = Multi(), progress = p -> nothing)
    yield() # prevents deadlocks, shouldn't be necessary
    easy = Easy()
    set_url(easy, req.url)
    for hdr in req.headers
        add_header(easy, hdr)
    end
    enable_progress(easy, true)
    add_handle(multi, easy)
    @sync begin
        @async for buf in easy.buffers
            write(req.io, buf)
        end
        @async for prog in easy.progress
            progress(prog)
        end
    end
    remove_handle(multi, easy)
    url = get_effective_url(easy)
    status = get_response_code(easy)
    response, headers = get_response_headers(easy)
    return Response(url, status, response, headers)
end

end # module


text:
# TODO: Remove when `MatrixReshaped` is removed
using Distributions, Test, Random, LinearAlgebra

rng = MersenneTwister(123456)

@testset "matrixreshaped.jl" begin

function test_matrixreshaped(rng, d1, sizes)
    x1 = rand(rng, d1)
    d1s = [@test_deprecated(MatrixReshaped(d1, s...)) for s in sizes]

@testset "MatrixReshaped $(nameof(typeof(d1))) tests" begin
    @testset "MatrixReshaped constructor" begin
        for d in d1s
            @test d isa MatrixReshaped
        end
    end
    @testset "MatrixReshaped constructor errors" begin
        @test_deprecated(@test_throws ArgumentError MatrixReshaped(d1, length(d1), 2))
        @test_deprecated(@test_throws ArgumentError MatrixReshaped(d1, length(d1)))
        @test_deprecated(@test_throws ArgumentError MatrixReshaped(d1, -length(d1), -1))
    end
    @testset "MatrixReshaped size" begin
        for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
            @test size(d) == s
        end
    end
    @testset "MatrixReshaped length" begin
        for d in d1s
            @test length(d) == length(d1)
        end
    end
    @testset "MatrixReshaped rank" begin
        for (d, s) in zip(d1s, sizes)
            @test rank(d) == minimum(s)
        end
    end
    @testset "MatrixReshaped insupport" begin
        for (i, d) in enumerate(d1s[1:end-1])
            for (j, s) in enumerate(sizes[1:end-1])
                @test (i == j) ⊻ !insupport(d, reshape(x1, s))
            end
        end
    end
    @testset "MatrixReshaped mean" begin
        for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
            @test mean(d) == reshape(mean(d1), s)
        end
    end
    @testset "MatrixReshaped mode" begin
        for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
            @test mode(d) == reshape(mode(d1), s)
        end
    end
    @testset "MatrixReshaped covariance" begin
        for (d, (n, p)) in zip(d1s[1:end-1], sizes[1:end-1])
            @test cov(d) == cov(d1)
            @test cov(d, Val(false)) == reshape(cov(d1), n, p, n, p)
        end
    end
    @testset "MatrixReshaped variance" begin
        for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
            @test var(d) == reshape(var(d1), s)
        end
    end
    @testset "MatrixReshaped params" begin
        for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
            @test params(d) == (d1, s)
        end
    end
    @testset "MatrixReshaped partype" begin
        for d in d1s
            @test partype(d) === partype(d1)
        end
    end
    @testset "MatrixReshaped eltype" begin
        for d in d1s
            @test eltype(d) === eltype(d1)
        end
    end
    @testset "MatrixReshaped logpdf" begin
        for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
            x = reshape(x1, s)
            @test logpdf(d, x) == logpdf(d1, x1)
        end
    end
    @testset "MatrixReshaped rand" begin
        for d in d1s
            x = rand(rng, d)
            @test insupport(d, x)
            @test insupport(d1, vec(x))
            @test logpdf(d, x) == logpdf(d1, vec(x))
        end
    end
    @testset "MatrixReshaped vec" begin
        for d in d1s
            @test vec(d) === d1
        end
    end
end
end

    # MvNormal
    σ = rand(rng, 16, 16)
    μ = rand(rng, 16)
    d1 = MvNormal(μ, σ * σ')
    sizes = [(4, 4), (8, 2), (2, 8), (1, 16), (16, 1), (4,)]
    test_matrixreshaped(rng, d1, sizes)

    # Dirichlet
    α = rand(rng, 36) .+ 1 # mode is only defined if all alpha > 1
    d1 = Dirichlet(α)
    sizes = [(6, 6), (4, 9), (9, 4), (3, 12), (12, 3), (1, 36), (36, 1), (6,)]
    test_matrixreshaped(rng, d1, sizes)
end


text:
function idmapping(nout,nin)
    nin == nout || @warn "Identity mapping not possible with nin != nout! Got nin=$nin, nout=$nout."
    return idmapping_nowarn(nout,nin)
end
idmapping_nowarn(nout, nin) = Matrix{Float32}(I, nout,nin)

function idmapping(h,w,nin,nout)
    nin == nout || @warn "Identity mapping not possible with nin != nout! Got nin=$nin, nout=$nout."
    w % 2 == 1 && h % 2 == 1 || @warn "Identity mapping requires odd kernel sizes! Got w=$w, h=$h."
    return idmapping_nowarn(h,w, nin,nout)
end
function idmapping_nowarn(h,w,nin,nout)
    center_w = h ÷ 2 + 1
    center_h = w ÷ 2 + 1
    weights = zeros(Float32,h,w,nin,nout)
    for i in 1:min(nin,nout)
        weights[center_w, center_h, i, i] = 1
    end
    return weights
end


text:
function build_primitive()
    # Build a simple metadata.
    metadata = Dict("mappable" => true)
    component = Component("primitive"; metadata = metadata)
    # Instantiate two output and two input ports.
    # Make the metadata for the two input ports slightly different to test routing
    add_port(component, "in", Input)
    add_port(component, "out", Output)

    check(component)
    return component
end

function build_tile()
    # Build a simple metadata.
    metadata = Dict("mappable" => false)
    component = Component("tile"; metadata = metadata)

    # Add general primitive
    add_child(component, build_primitive(), "u0")

    # Instantiate two output and two input ports.
    # Make the metadata for the two input ports slightly different to test routing
    dirs = ("north", "east", "south", "west", "up", "down")
    for dir in dirs
        add_port(component, "$(dir)_in", Input)
        add_port(component, "$(dir)_out", Output)
    end

    # instantiate a mux
    add_child(component, build_mux(7, 7), "mux")

    ## Links

    # Component to mux
    add_link(component, "u0.out", "mux.in[0]")
    add_link(component, "mux.out[0]", "u0.in")

    # Mux to tile i/o
    for (i, d) in enumerate(dirs)
        add_link(component, "mux.out[$i]", "$(d)_out")
        add_link(component, "$(d)_in", "mux.in[$i]")
    end

    check(component)
    return component
end

function build_arch()
    a = TopLevel{3}("arch_3d")
    tile = build_tile()

    # Instantiate tiles
    for (i, j, k) in Iterators.product(1:3, 1:3, 1:3)
        add_child(a, tile, CartesianIndex(i, j, k))
    end

    offsets = [
        Offset((-1, 0, 0), "north_out", "south_in"),
        Offset((1, 0, 0), "south_out", "north_in"),
        Offset((0, 1, 0), "east_out", "west_in"),
        Offset((0, -1, 0), "west_out", "east_in"),
        Offset((0, 0, 1), "up_out", "down_in"),
        Offset((0, 0, -1), "down_out", "up_in"),
    ]

    connection_rule(a, offsets)
    check(a)
    return a
end


text:
# Official spec : http://www.libpng.org/pub/png/spec/1.2/png-1.2.pdf
module PNGImage

using FileIO: @format_str, Stream, save
using ImageCore: Colorant, Gray, GrayA, RGB, RGBA, Normed

export PngImage, png_color_T

# workaround for FileIO pre 1.6
include("FileIO_workaround.jl")

"""
    PngImage

A wrapper around the binary encoding of a PNG image, holding it's attributes
"""
struct PngImage
    attr::Dict
    data::Vector{UInt8}
end

function PngImage(buffer::Base.IOBuffer)
    !check_valid_png(buffer.data) && return false
    attr  = read_info(buffer.data)
    return PngImage(attr, buffer.data)
end

"""
    check_valid_png(data)

Checks if the header of the provided buffer is a valid png header.
Also checks if the first chunk is valid according to the PNG standard

Refs:
[1] http://www.libpng.org/pub/png/spec/1.2/png-1.2.pdf
[2] https://en.wikipedia.org/wiki/Portable_Network_Graphics#File_header
"""
function check_valid_png(data)
    #Header https://en.wikipedia.org/wiki/Portable_Network_Graphics#File_header
    header_value = UInt8[0x89, 0x50, 0x4e, 0x47, 0x0d, 0x0a, 0x1a, 0x0a]

    # Avoid accessing memory that is not there
    length(data) < 12 && return false
    # Check the header
    data[1:8] != header_value && return false

    # First chunk: check length
    cnk_start = 8;
    cnk_len = reinterpret(Int32, data[cnk_start.+(1:4)]) |> first |> ntoh
    String(data[cnk_start.+(5:8)]) != "IHDR" && return false
    cnk_len <13 && return false

    return true
end

"""
    read_info(data)

Reads the content of the first chunk of a PNG, containing information such as
height, width and color type.
"""
function read_info(data)
    cnk_start = 8;
    cnk_data = cnk_start + 8;

    # Read chunk data
    width, height = ntoh.(reinterpret(Int32, data[cnk_data.+(1:8)]))
    bit_depth     = ntoh.(reinterpret(Int8,  data[cnk_data+9]))
    col_type      = ntoh.(reinterpret(Int8,  data[cnk_data+10]))
    compression   = ntoh.(reinterpret(Int8,  data[cnk_data+11]))
    filter        = ntoh.(reinterpret(Int8,  data[cnk_data+12]))
    interlace     = ntoh.(reinterpret(Int8,  data[cnk_data+12]))

    attr = Dict(
        :width=>width,
        :height=>height,
        :bitdepth=>bit_depth,
        :colorspace=>col_type,
        :compression=>compression,
        :filter=>filter,
        :interlace=>interlace
    )
    return attr
end

function Base.convert(::Type{PngImage}, img::AbstractArray{<:Colorant})
    pb = PipeBuffer()
    save(_format_stream(format"PNG", pb), img)
    return PngImage(pb)
end

# Utility function, for seeing those objects in Juno/iJulia.
Base.show(io::IO, ::MIME"image/png", x::PngImage) =
    write(io, x.data)

# Things for deserializing
const PNG_COLOR_MASK_PALETTE = 1
const PNG_COLOR_MASK_COLOR = 2
const PNG_COLOR_MASK_ALPHA = 4

const PNG_COLOR_TYPE_GRAY = 0
const PNG_COLOR_TYPE_PALETTE = PNG_COLOR_MASK_COLOR | PNG_COLOR_MASK_PALETTE
const PNG_COLOR_TYPE_RGB = PNG_COLOR_MASK_COLOR
const PNG_COLOR_TYPE_RGB_ALPHA = PNG_COLOR_MASK_COLOR | PNG_COLOR_MASK_ALPHA
const PNG_COLOR_TYPE_GRAY_ALPHA = PNG_COLOR_MASK_ALPHA
const PNG_COLOR_TYPE_RGBA = PNG_COLOR_TYPE_RGB_ALPHA
const PNG_COLOR_TYPE_GA = PNG_COLOR_TYPE_GRAY_ALPHA

function png_color_T(attributes)
    color_type = attributes[:colorspace]
    bit_depth  = attributes[:bitdepth]

    if color_type == PNG_COLOR_TYPE_GRAY
        colors_type = Gray{bit_depth > 8 ? Normed{UInt16,bit_depth} : Normed{UInt8,bit_depth}}
    elseif color_type == PNG_COLOR_TYPE_PALETTE
        colors_type = RGB{bit_depth == 16 ? N0f16 : N0f8}
    elseif color_type == PNG_COLOR_TYPE_RGB
        colors_type = RGB{bit_depth == 16 ? N0f16 : N0f8}
    elseif color_type == PNG_COLOR_TYPE_RGB_ALPHA
        colors_type = RGBA{bit_depth == 16 ? N0f16 : N0f8}
    elseif color_type == PNG_COLOR_TYPE_GRAY_ALPHA
        colors_type = GrayA{bit_depth > 8 ? Normed{UInt16,bit_depth} : Normed{UInt8,bit_depth}}
    else
        throw(error("Unknown color type: $color_type"))
    end
    return colors_type
end

end


text:
# # Leaf Level Photosynthesis
#
# Here, we will go one step further and look at the entire leaf-level response, including the impact of stomatal conductance as well.

# ## Stomatal Conductance
# Before, we focussed mainly on the demand-driven constraints through Rubisco and RuBP regeneration. Stomatal conductance is highly important as it drives the suppy-side of photosynthesis and interacts with the energy balance and leaf temperature as latent heat fluxes are a major factor in surface cooling.
#
# Before, we have derived net rates of photosynthesis $A_n$, which have to be matched with the supply side through CO$_2$ diffusion:
#
# $$A_n^{diff} = g_{leaf,C}(C_a-C_c)$$
#
# which can be separated into diffusion from the air to the leaf surface with a boundary layer conductance $g_{b,C}$, diffusion from the surface to the interstitial air-space with stomatal conductance $g_{s,C}$ and diffusion from the  interstitial air-space to the chloroplast with mesophyll conductance $g_{m,C}$ (the letter C stands for CO2 here, as the diffusion constants vary with species, e.g. H$_2$O stomatal conductance is a factor 1.6 higher than those for CO$_2$):
#
# $$A_n^{diff} = g_{b,C}(C_a-C_s) = g_{s,C}(C_s-C_i) = g_{m,C}(C_i-C_c) = g_{leaf,C}(C_a-C_c)$$
#
# $$g_{leaf,C} = \left(\frac{1}{g_{b,C}+g_{s,C}+g_{m,C}}\right)^{-1}$$
#
# The importance of stomatal conductance and its interplay with photosynthesis here is that the supply and demand rates determine $C_c$. A reduction in $C_c$ reduces the demand-limited rates while it increase the diffusion rates (at a given $g$). Most models run internal so-called A-$C_c$ iterations to ensure both rates are balanced and in steady-state ($\partial C_c/\partial t=0$). We implement this option but also opt to run stomatal conductance prognostically, as stomata don't open and close instantanously but have a response time of around 15 minutes.
#
# ### Empirical Stomatal Conductance models
#
# Currently, we can choose between the widely used Ball-Berry model[^1] or the Medlyn approach[^2].
#
# Ball and Berry derived the following empirical formulation based on leaf level data (RH is relative humidity):
#
# $$g_{s,w} = g_0 + g_1  \frac{A_n \cdot RH}{C_s}$$
#
# Medlyn derived the following equations based on stomatal optimization theory but in a similar form as the BB model:
#
# $$g_{s,w} = g_0 + \left(1 + \frac{g_1}{\sqrt{VPD}}\right) \frac{A_n}{C_s}$$
#
# Both models are proportional to $A_n$ and inversely proportional to $C_s$, with the main difference in the dependence on either relative humidity in the Ball-Berry model vs. vapor pressure deficit (VPD) in Medlyn et al. In both cases, $g_0$ is the residual conductivity even if stomata are fully closed and $g_1$ is related to the marginal water cost of plant carbon gain. Importantly, $g_1$ can't be inter-changed between the formulations, underlining again that parameters have to be optimized with respect to the model that is eventually being used.
#
# ### Stomatal Optimization Theories
#
# The empirical formulations only hold in well-watered soil and our main goal is to implement stomatal optimization models to take the entire soil-plant-atmosphere continuum into account[^3]. Here, we will just use the empirical models and steady-state photosynthesis to show the underlying principles.
#
# [^1]: Ball, J.T., Woodrow, I.E. and Berry, J.A., 1987. A model predicting stomatal conductance and its contribution to the control of photosynthesis under different environmental conditions. In Progress in photosynthesis research (pp. 221-224). Springer, Dordrecht.
#
# [^2]: Medlyn, B.E., Duursma, R.A., Eamus, D., Ellsworth, D.S., Prentice, I.C., Barton, C.V., Crous, K.Y., De Angelis, P., Freeman, M. and Wingate, L., 2011. Reconciling the optimal and empirical approaches to modelling stomatal conductance. Global Change Biology, 17(6), pp.2134-2144.
#
# [^3]: Wang, Y., Sperry, J.S., Anderegg, W.R., Venturas, M.D. and Trugman, A.T., 2020. A theoretical and empirical assessment of stomatal optimization modeling. New Phytologist.
#

## Add usual tools we use:
##using Revise
using BenchmarkTools
using PyPlot
#----------------------------------------------------------------------------

## load Photosynthesis module:
using Land.Photosynthesis
#----------------------------------------------------------------------------

## Specify Field Type
const FT = Float64
#----------------------------------------------------------------------------

## Create a standard leaf with defualt parameters
leaf3 = Leaf{FT}(APAR=1200, Vcmax25=90, Jmax25=90*1.9, Vpmax25=100, Rd25=1);
leaf4 = Leaf{FT}(APAR=1200, Vcmax25=90, Jmax25=90*1.9, Vpmax25=100, Rd25=1);
## Create a standard meteo structure:
envir  = AirLayer{FT}();
#----------------------------------------------------------------------------

## Setting some standard values (dynamic-state=false forces A-Cc iterations)
##leaf.dynamic_state = false
#----------------------------------------------------------------------------

# ---
# ## Modular stucture:
# The most important step is to define which submodules to use. There might be different implementations for Fluorescence, Photosynthesis (C3,C4,etc), respiration, stomatal conductance (e.g. Ball-Berry, Medlyn), the T-dependence of J$_{max}$, V$_{c,max}$ and Michaelis Menten constants as well as leaf boundary layer resistance (setting it to 0 here to mimic well vented laboratory leaf level measurements) and colimitation.
#
# Below, we create two different model setups with a C3 and C4 photosynthesis pathway. Note that we chose the Ball Berry stomatal conductance model with different $g_1$ for C3 and C4 (C4 typically much lower).

## use this as the boundary layer resistance (1/gb)
## deprecated as there are g_bw and g_bc terms in Leaf struct
## ra = FT(0.5)

## C3 Photosynthesis
# Use pre-defined parameter sets
modC3 = C3CLM(FT)
#modC3.Sto = Photosynthesis.OSMWang()
modC3.Sto = Photosynthesis.ESMBallBerry{FT}(g1 = 16)

## C4 Photosynthesis
modC4 = C4CLM(FT)
#modC4.Sto = Photosynthesis.OSMWang()
modC4.Sto = Photosynthesis.ESMBallBerry{FT}(g1 = 8)
#----------------------------------------------------------------------------

# ## Light response curves
# Now we can compute light response curves for C3 and C4 photosynthesis

##Again, this looks tedious and repetitive but is the easiest way to do this right now:

## Variable I want to save:
Ag_C3 = Float32[]; Ag_C4 = Float32[]
An_C3 = Float32[]; An_C4 = Float32[]
Aj_C3 = Float32[]; Aj_C4 = Float32[]
Ap_C3 = Float32[]; Ap_C4 = Float32[]
Ac_C3 = Float32[]; Ac_C4 = Float32[]
Cc_C3 = Float32[]; Cc_C4 = Float32[]
gs_C3 = Float32[]; gs_C4 = Float32[]

APAR = collect(FT, 0:10:1700)
for _APAR in APAR
    leaf3.APAR = _APAR;
    leaf4.APAR = _APAR;
    ## Run C3 photosynthesis (general photosynthesis model):
    leaf_photo_from_envir!(modC3, leaf3, envir, modC3.Sto);
    ## Save leaf variables:
    push!(An_C3, leaf3.An); push!(Ag_C3, leaf3.Ag);
    push!(Aj_C3, leaf3.Aj); push!(Ap_C3, leaf3.Ap);
    push!(Ac_C3, leaf3.Ac); push!(Cc_C3, leaf3.p_i);
    push!(gs_C3, leaf3.g_sw);

    ## Run C4 photosynthesis:
    leaf_photo_from_envir!(modC4, leaf4, envir, modC4.Sto);
    ## Save leaf variables:
    push!(An_C4, leaf4.An); push!(Ag_C4, leaf4.Ag);
    push!(Aj_C4, leaf4.Aj); push!(Ap_C4, leaf4.Ap);
    push!(Ac_C4, leaf4.Ac); push!(Cc_C4, leaf4.p_i);
    push!(gs_C4, leaf4.g_sw);
end
#----------------------------------------------------------------------------

## Testing some times, how long does this take (as we need to run it globally, it has to be efficient)?
## Slow for now, because of unnecessary allocations
## Will improve when the structs are cleaned up
#@btime leaf_photo_from_envir!(modC3, leaf3, envir, modC3.Sto);
#@btime leaf_photo_from_envir!(modC4, leaf4, envir, modC4.Sto);
#----------------------------------------------------------------------------

# ## C3 Light Response Curve

##plot(APAR, An,  label="An")
figure()
plot(APAR, Ag_C3,color=:black,lw=2, alpha=0.7, label="Ag C3")
plot(APAR, Ac_C3, ls="--", lw=2, label="Ac C3")
plot(APAR, Aj_C3, ls="--", lw=2, label="Aj C3")
plot(APAR, Ap_C3, ls="--", lw=2, label="Ap C3" )
xlabel("APAR [μmol/m2/s]")
ylabel("Aᵢ [μmol/m2/s]")
title("C3 photosynthesis light response")
legend()
gcf()
#----------------------------------------------------------------------------

# ## C4 Light Response Curve

##plot(APAR, An,  label="An")
figure()
plot(APAR, Ag_C4,color=:black, lw=2, alpha=0.7, label="Ag C4")
plot(APAR, Ac_C4, ls="--", lw=2, label="Ac C4")
plot(APAR, Aj_C4, ls="--", lw=2, label="Aj C4")
plot(APAR, Ap_C4, ls="--", lw=2, label="Ap C4" )
xlabel("APAR [μmol/m2/s]")
ylabel("Aᵢ [μmol/m2/s]")
title("C4 photosynthesis light response")
legend()
gcf()
#----------------------------------------------------------------------------

# ### Dependence of $C_c$ on light
# Increasing light should result in overall lower $C_c$ concentrations. C4 plants can work with much lower $C_c$ as for C4, the concentration in the bundle sheath cells count, not the mesophyll.

figure()
plot(APAR, Cc_C4/envir.p_a,color=:black ,lw=2, alpha=0.7, label="Cc/Ca C4")
plot(APAR, Cc_C3/envir.p_a,color=:orange,lw=2, alpha=0.7, label="Cc/Ca C3")
xlabel("APAR [μmol/m²/s]")
ylabel("Cc/Ca [-]")
legend()
gcf()
#----------------------------------------------------------------------------

# ### Dependence of stomatal coductance g$_s$

#This part has been broken by Yuije, please fix it
figure()
plot(APAR, gs_C4,color=:black ,lw=2, alpha=0.7, label="gs C4")
plot(APAR, gs_C3,color=:orange,lw=2, alpha=0.7, label="gs C3")
xlabel("APAR [μmol/m²/s]")
ylabel("gs")
legend()
gcf()
#----------------------------------------------------------------------------

# ## More complex examples:

##Again, this looks tedious and repetitive but is the easiest way to do this right now:

## Now I want to vary T, APAR and CO2:
## Start CO2 from 100 ppm to make sure it is higher than Γ*
APAR = [100.0, 250.0, 500.0, 1000.0, 1500.0]
CO2  = collect(100:20:800)
T    = collect(274:2:310)

n1 = length(APAR);
n2 = length(CO2);
n3 = length(T);

Ag_C3 = zeros(n1,n2,n3); Ag_C4 = zeros(n1,n2,n3)
An_C3 = zeros(n1,n2,n3); An_C4 = zeros(n1,n2,n3)
Aj_C3 = zeros(n1,n2,n3); Aj_C4 = zeros(n1,n2,n3)
Ap_C3 = zeros(n1,n2,n3); Ap_C4 = zeros(n1,n2,n3)
Ac_C3 = zeros(n1,n2,n3); Ac_C4 = zeros(n1,n2,n3)
Cc_C3 = zeros(n1,n2,n3); Cc_C4 = zeros(n1,n2,n3)
gs_C3 = zeros(n1,n2,n3); gs_C4 = zeros(n1,n2,n3)

#----------------------------------------------------------------------------

## Run this over all potential 3D dimensions:

## I really like the compact form of nested loops in Julia!
for iA in eachindex(APAR), iC in eachindex(CO2), iT in eachindex(T)
    #println(iA, "/", iC, "/", iT)
    envir.p_a  = CO2[iC]/10;
    leaf3.T    = T[iT];
    leaf3.APAR = APAR[iA];
    leaf4.T    = T[iT];
    leaf4.APAR = APAR[iA];

    ## Run C3 photosynthesis:
    leaf_photo_from_envir!(modC3, leaf3, envir, modC3.Sto);

    ## Save leaf variables:
    An_C3[iA,iC,iT]=leaf3.An;
    Ag_C3[iA,iC,iT]=leaf3.Ag;
    Aj_C3[iA,iC,iT]=leaf3.Aj;
    Ap_C3[iA,iC,iT]=leaf3.Ap;
    Ac_C3[iA,iC,iT]=leaf3.Ac;
    Cc_C3[iA,iC,iT]=leaf3.p_i;
    gs_C3[iA,iC,iT]=leaf3.g_sw;

    ## Run C4 photosynthesis:
    leaf_photo_from_envir!(modC4, leaf4, envir, modC4.Sto);

    ## Save leaf variables:
    An_C4[iA,iC,iT]=leaf4.An;
    Ag_C4[iA,iC,iT]=leaf4.Ag;
    Aj_C4[iA,iC,iT]=leaf4.Aj;
    Ap_C4[iA,iC,iT]=leaf4.Ap;
    Ac_C4[iA,iC,iT]=leaf4.Ac;
    Cc_C4[iA,iC,iT]=leaf4.p_i;
    gs_C4[iA,iC,iT]=leaf4.g_sw;
end
#----------------------------------------------------------------------------

# #### Plotting 2D contours of A$_n$ vs. CO$_2$ and Temperature

##Let's take one slice in APAR space:
i = 4

## and plot:
figure()
contourf(T.-273.15, CO2, An_C3[i,:,:])
xlabel("T [°C]")
ylabel("Ambient CO₂ [ppm]")
title("C3 , An [μmol/m²/s] at APAR=$(APAR[i])")
colorbar()
gcf()
#----------------------------------------------------------------------------

## Same for C4 plants, why is it so different??

figure()
contourf(T.-273.15, CO2, Cc_C4[i,:,:])
xlabel("T [°C]")
ylabel("Ambient CO₂ [ppm]")
title("C4 , Cc [Pa] at APAR=$(APAR[i])")
colorbar()
gcf()
#----------------------------------------------------------------------------

## Same for C4 plants, why is it so different??

figure()
contourf(T.-273.15, CO2, An_C4[i,:,:])
xlabel("T [°C]")
ylabel("Ambient CO₂ [ppm]")
title("C4 , An [μmol/m²/s] at APAR=$(APAR[i])")
colorbar()
gcf()
#----------------------------------------------------------------------------

# ## Ambient CO$_2$ response curves and limiting rates for C3 and C4

iA = 4; iT=12
figure()
plot( CO2, Ag_C3[iA,:,iT],color=:black,lw=2, alpha=0.7, label="Ag C3")
plot(CO2, Ac_C3[iA,:,iT], ls="--", lw=2, label="Ac C3")
plot(CO2, Aj_C3[iA,:,iT], ls="--", lw=2, label="Aj C3")
plot(CO2, Ap_C3[iA,:,iT], ls="--", lw=2, label="Ap C3" )
xlabel("CO₂ [ppm]")
ylabel("Aᵢ [μmol/m2/s]")
title("Ambient C3 CO₂ response, T=$(T[iT]-273), APAR=$(APAR[iA])")
legend()
gcf()

#----------------------------------------------------------------------------

iA = 4; iT=12
figure()
plot(CO2, Ag_C4[iA,:,iT],color=:black,lw=2, alpha=0.7, label="Ag C4")
plot(CO2, Ac_C4[iA,:,iT], ls="--", lw=2, label="Ac C4")
plot(CO2, Aj_C4[iA,:,iT], ls="--", lw=2, label="Aj C4")
plot(CO2, Ap_C4[iA,:,iT], ls="--", lw=2, label="Ap C4")
xlabel("CO₂ [ppm]")
ylabel("Aᵢ [μmol/m2/s]")
title("Ambient C4 CO₂ response, T=$(T[iT]-273), APAR=$(APAR[iA])")
legend()
gcf()

#----------------------------------------------------------------------------

# ---
# ## Summary:
# Play around a bit with different parameters. You can also try to change humidity at the leaf level or the boundary layer resistance. There are ## many things that can be tested, and we often like to go to the extremes (T, CO$_2$, ra) to check for bugs, which most often have to do with the ## A-C$_c$ iterations.
#
# Note of caution: we really plotted CO$_2$ response curves against ambient CO$_2$. At low C$_c$, leaves can be a source of CO$_2$ due to the CO$_2$ compensation point (i.e. A$_g$ can be negative!). Also, at CO$_2$=0, stomatal conductance using the empirical formulation goes to $\infty$, which can cause unwanted behavior.
#
# Test things out, learn by changing parameters and try to break it!

#----------------------------------------------------------------------------

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
# LibFTD2XX.jl - High Level Module
#
# By Reuben Hill 2019, Gowerlabs Ltd, reuben@gowerlabs.co.uk
#
# Copyright (c) Gowerlabs Ltd.
#
# This module contains methods and functions for interacting with D2XX devices.
# It calls functions from the submodule `Wrapper` which in turn call Functions
# from the FT D2XX library. See D2XX Programmer's Guide (FT_000071) for more
# information on that library.

module LibFTD2XX

# Enums
export FTOpenBy, OPEN_BY_SERIAL_NUMBER, OPEN_BY_DESCRIPTION
export FTWordLength, BITS_8, BITS_7
export FTStopBits, STOP_BITS_1, STOP_BITS_2
export FTParity, PARITY_NONE, PARITY_ODD, PARITY_EVEN, PARITY_MARK, PARITY_SPACE
export FTBitMode, MODE_RESET, MODE_ASYNC_BITBANG, MODE_MPSSE, MODE_SYNC_BITBANG,
       MODE_MCU_EMULATION, MODE_FAST_OPTO, MODE_CBUS_BITBANG, MODE_SCS_FIFO

# Types and Constructors
export FT_HANDLE # Exported by .Wrapper
export D2XXException
export D2XXDevice, D2XXDevices

# Port communication functions
export baudrate, datacharacteristics, timeouts, status, driverversion

# Library Functions
export libversion

# D2XXDevice Accessor Functions
export deviceidx, deviceflags, devicetype, deviceid, locationid, serialnumber,
       description, fthandle

include("util.jl")
include("wrapper.jl")

using .Wrapper

"""
    @enum(
      FTOpenBy,
      OPEN_BY_SERIAL_NUMBER = FT_OPEN_BY_SERIAL_NUMBER,
      OPEN_BY_DESCRIPTION = FT_OPEN_BY_DESCRIPTION)

For use with [`open`](@ref).
"""
@enum(
  FTOpenBy,
  OPEN_BY_SERIAL_NUMBER = FT_OPEN_BY_SERIAL_NUMBER,
  OPEN_BY_DESCRIPTION = FT_OPEN_BY_DESCRIPTION)

"""
    @enum(
      FTWordLength,
      BITS_8 = FT_BITS_8,
      BITS_7 = FT_BITS_7)

For use with [`datacharacteristics`](@ref).
"""
@enum(
  FTWordLength,
  BITS_8 = FT_BITS_8,
  BITS_7 = FT_BITS_7)

"""
    @enum(
      FTStopBits,
      STOP_BITS_1 = FT_STOP_BITS_1,
      STOP_BITS_2 = FT_STOP_BITS_2)

For use with [`datacharacteristics`](@ref).
"""
@enum(
  FTStopBits,
  STOP_BITS_1 = FT_STOP_BITS_1,
  STOP_BITS_2 = FT_STOP_BITS_2)

"""
    @enum(
      FTParity,
      PARITY_NONE = FT_PARITY_NONE,
      PARITY_ODD  = FT_PARITY_ODD,
      PARITY_EVEN = FT_PARITY_EVEN,
      PARITY_MARK = FT_PARITY_MARK,
      PARITY_SPACE = FT_PARITY_SPACE)

For use with [`datacharacteristics`](@ref).
"""
@enum(
  FTParity,
  PARITY_NONE = FT_PARITY_NONE,
  PARITY_ODD  = FT_PARITY_ODD,
  PARITY_EVEN = FT_PARITY_EVEN,
  PARITY_MARK = FT_PARITY_MARK,
  PARITY_SPACE = FT_PARITY_SPACE)

"""
    @enum(
      FTBitMode,
      MODE_RESET = FT_MODE_RESET,
      MODE_ASYNC_BITBANG = FT_MODE_ASYNC_BITBANG,
      MODE_MPSSE = FT_MODE_MPSSE,
      MODE_SYNC_BITBANG = FT_MODE_SYNC_BITBANG,
      MODE_MCU_EMULATION = FT_MODE_MCU_EMULATION,
      MODE_FAST_OPTO = FT_MODE_FAST_OPTO,
      MODE_CBUS_BITBANG = FT_MODE_CBUS_BITBANG,
      MODE_SCS_FIFO = FT_MODE_SCS_FIFO)

For use with [`bitmode`](@ref).
"""
@enum(
  FTBitMode,
  MODE_RESET = FT_MODE_RESET,
  MODE_ASYNC_BITBANG = FT_MODE_ASYNC_BITBANG,
  MODE_MPSSE = FT_MODE_MPSSE,
  MODE_SYNC_BITBANG = FT_MODE_SYNC_BITBANG,
  MODE_MCU_EMULATION = FT_MODE_MCU_EMULATION,
  MODE_FAST_OPTO = FT_MODE_FAST_OPTO,
  MODE_CBUS_BITBANG = FT_MODE_CBUS_BITBANG,
  MODE_SCS_FIFO = FT_MODE_SCS_FIFO)

"""
    D2XXException <: Exception

LibFTD2XX High-Level Library Error Type.
"""
struct D2XXException <: Exception
  str::String
end

"""
    struct D2XXDevice <: IO

Device identifier for a D2XX device.

See also: [`D2XXDevices`](@ref), [`deviceidx`](@ref), [`deviceflags`](@ref),
[`devicetype`](@ref), [`deviceid`](@ref), [`locationid`](@ref),
[`serialnumber`](@ref), [`description`](@ref), [`fthandle`](@ref).
"""
struct D2XXDevice <: IO
  idx::Int
  flags::Int
  typ::Int
  id::Int
  locid::Int
  serialnumber::String
  description::String
  fthandle::Ref{FT_HANDLE}
end

# D2XXDevice Constructors
#
"""
    D2XXDevice(deviceidx::Integer)

Construct a `D2XXDevice` without opening it. D2XX hardware must pre present to
work.
"""
D2XXDevice(deviceidx::Integer) = D2XXDevice(getdeviceinfodetail(deviceidx)...)

"""
    D2XXDevices()

Returns an array of available D2XX devices of type `D2XXDevice`.

Their state is not modified - if they are already open/closed they remain
open/closed.

See also: [`D2XXDevice`](@ref), [`open`](@ref)
"""
function D2XXDevices()
  numdevs = createdeviceinfolist() # NB numdevs is DWORD = Cuint
  devices = D2XXDevice[]
  for devidx = 0:(Int(numdevs)-1) # Int conversion to avoid devidx = 0x00000000:0xffffffff on Win32
    push!(devices, D2XXDevice(devidx))
  end
  devices
end

# Port communication functions
#
"""
    isopen(d::D2XXDevice) -> Bool

See also: [`D2XXDevice`](@ref)
"""
Base.isopen(d::D2XXDevice) = isopen(fthandle(d))

"""
    isopen(handle::FT_HANDLE) -> Bool

See also: [`FT_HANDLE`](@ref)
"""
function Base.isopen(handle::FT_HANDLE)
  open = true
  if Wrapper._ptr(handle) == C_NULL
    open = false
  else
    try
      FT_GetModemStatus(handle)
    catch ex
      if ex == FT_INVALID_HANDLE
        open = false
      else
        rethrow(ex)
      end
    end
  end
  open
end

"""
    Base.open(d::D2XXDevice)

Open a [`D2XXDevice`](@ref) for reading and writing using [`FT_OpenEx`](@ref).
Cannot be used to open the same device twice.

See also: [`isopen`](@ref), [`close`](@ref)
"""
function Base.open(d::D2XXDevice)
  isopen(d) && throw(D2XXException("Device already open."))
  fthandle(d, FT_Open(deviceidx(d)))
  return
end

"""
    open(str::AbstractString, openby::FTOpenBy) -> FT_HANDLE

Create an open [`FT_HANDLE`](@ref) for reading and writing using
[`FT_OpenEx`](@ref). Cannot be used to open the same device twice.

# Arguments
 - `str::AbstractString` : Device identifier. Type depends on `openby`
 - `openby::FTOpenBy` : Indicator of device identifier `str` type.

See also: [`isopen`](@ref), [`close`](@ref)
"""
Base.open(str::AbstractString, openby::FTOpenBy) =  FT_OpenEx(str, DWORD(openby))

"""
    close(d::D2XXDevice)

Close a [`D2XXDevice`](@ref) using [`FT_Close`](@ref). Does not perform a
[`flush`](@ref) first.
"""
Base.close(d::D2XXDevice) = close(fthandle(d))

"""
    close(handle::FT_HANDLE)

Close an [`FT_HANDLE`](@ref) using [`FT_Close`](@ref). Does not perform a
[`flush`](@ref) first.
"""
function Base.close(handle::FT_HANDLE)
  if isopen(handle)
    FT_Close(handle)
  end
  return
end

"""
    bytesavailable(d::D2XXDevice)

See also: [`D2XXDevice`](@ref), [`isopen`](@ref), [`open`](@ref),
[`readavailable`](@ref), [`read`](@ref)
"""
Base.bytesavailable(d::D2XXDevice) = bytesavailable(fthandle(d))

"""
    bytesavailable(handle::FT_HANDLE)

See also: [`FT_HANDLE`](@ref), [`isopen`](@ref), [`open`](@ref),
[`readavailable`](@ref), [`read`](@ref)
"""
function Base.bytesavailable(handle::FT_HANDLE)
  isopen(handle) || throw(D2XXException("Device must be open to check bytes available."))
  FT_GetQueueStatus(handle)
end

"""
    eof(d::D2XXDevice) -> Bool

Indicates if any bytes are available to be read from an open
[`D2XXDevice`](@ref). Non-blocking.

See also: [`isopen`](@ref), [`open`](@ref),
[`readavailable`](@ref), [`read`](@ref)
"""
Base.eof(d::D2XXDevice) = eof(fthandle(d))

"""
    eof(d::D2XXDevice) -> Bool

Indicates if any bytes are available to be read from an open
[`FT_HANDLE`](@ref). Non-blocking.

See also: [`isopen`](@ref), [`open`](@ref),
[`readavailable`](@ref), [`read`](@ref)
"""
Base.eof(handle::FT_HANDLE) = (bytesavailable(handle) == 0)

"""
    readbytes!(d::D2XXDevice, b::AbstractVector{UInt8}, nb=length(b))

See description for
[`readbytes!(stream::IO, b::AbstractVector{UInt8}, nb=length(b))`](@ref).
`d` must be open. Uses [`FTRead`](@ref).

See also: [`D2XXDevice`](@ref).
"""
Base.readbytes!(d::D2XXDevice, b::AbstractVector{UInt8}, nb=length(b)) =
readbytes!(fthandle(d), b, nb)

"""
    readbytes!(handle::FT_HANDLE, b::AbstractVector{UInt8}, nb=length(b))

See description for
[`readbytes!(stream::IO, b::AbstractVector{UInt8}, nb=length(b))`](@ref).
`handle` must be open. Uses [`FTRead`](@ref).

See also: [`FT_HANDLE`](@ref)
"""
function Base.readbytes!(handle::FT_HANDLE, b::AbstractVector{UInt8}, nb=length(b))
  isopen(handle) || throw(D2XXException("Device must be open to read."))
  if length(b) < nb
    resize!(b, nb)
  end
  nbrx = FT_Read(handle, b, nb)
end

"""
    readavailable(d::D2XXDevice)

Read all available data from an open [`D2XXDevice`](@ref). Does not block if
nothing is available.
"""
Base.readavailable(d::D2XXDevice) = readavailable(fthandle(d))

"""
    readavailable(handle::FT_HANDLE)

Read all available data from an open [`FT_HANDLE`](@ref). Does not block if
nothing is available.

See also: [`readbytes`](@ref) [`isopen`](@ref), [`open`](@ref)
"""
function Base.readavailable(handle::FT_HANDLE)
  b = Vector{UInt8}(undef, bytesavailable(handle))
  readbytes!(handle, b)
  b
end

"""
    write(d::D2XXDevice, buffer::Vector{UInt8})

Write `buffer` to an open [`D2XXDevice`](@ref) using [`FT_Write`](@ref).
"""
Base.write(d::D2XXDevice, buffer::Vector{UInt8}) = write(fthandle(d), buffer)

"""
    write(handle::FT_HANDLE, buffer::Vector{UInt8})

Write `buffer` to an open [`FT_HANDLE`](@ref) using [`FT_Write`](@ref).

See also: [`isopen`](@ref), [`open`](@ref)
"""
function Base.write(handle::FT_HANDLE, buffer::Vector{UInt8})
  isopen(handle) || throw(D2XXException("Device must be open to write."))
  FT_Write(handle, buffer, length(buffer))
end

"""
    baudrate(d::D2XXDevice, baud)

Set the baudrate of an open [`D2XXDevice`](@ref) using [`FT_SetBaudRate`](@ref).
"""
baudrate(d::D2XXDevice, baud) = baudrate(fthandle(d), baud)

"""
    baudrate(handle::FT_HANDLE, baud)

Set the baudrate of an open [`FT_HANDLE`](@ref) to `baud` using
[`FT_SetBaudRate`](@ref).

See also: [`isopen`](@ref), [`open`](@ref)
"""
function baudrate(handle::FT_HANDLE, baud)
  0 < baud || throw(DomainError("0 <= baud"))
  isopen(handle) || throw(D2XXException("Device must be open to set baudrate."))
  FT_SetBaudRate(handle, baud)
end

"""
    datacharacteristics(d::D2XXDevice;
                        wordlength::FTWordLength = BITS_8,
                        stopbits::FTStopBits = STOP_BITS_1,
                        parity::FTParity = PARITY_NONE)

Set the transmission and reception data characteristics for an open
[`D2XXDevice`](@ref) using [`FT_SetDataCharacteristics`](@ref).

# Arguments
 - `wordlength::FTWordLength` : Either BITS_7 or BITS_8
 - `stopbits::FTStopBits` : Either STOP_BITS_1 or STOP_BITS_2
 - `parity::FTParity` : PARITY_NONE, PARITY_ODD, PARITY_EVEN, PARITY_MARK, or
   PARITY_SPACE
"""
datacharacteristics(d::D2XXDevice;
                    wordlength::FTWordLength = BITS_8,
                    stopbits::FTStopBits = STOP_BITS_1,
                    parity::FTParity = PARITY_NONE) =
datacharacteristics(fthandle(d),
                    wordlength=wordlength,
                    stopbits=stopbits,
                    parity=parity)

"""
    datacharacteristics(handle::FT_HANDLE;
                        wordlength::FTWordLength = BITS_8,
                        stopbits::FTStopBits = STOP_BITS_1,
                        parity::FTParity = PARITY_NONE)

Set the transmission and reception data characteristics for an open
[`FT_HANDLE`](@ref) using [`FT_SetDataCharacteristics`](@ref).

See also: [`isopen`](@ref), [`open`](@ref)
"""
function datacharacteristics(handle::FT_HANDLE;
                             wordlength::FTWordLength = BITS_8,
                             stopbits::FTStopBits = STOP_BITS_1,
                             parity::FTParity = PARITY_NONE)
isopen(handle) || throw(D2XXException("Device must be open to set data characteristics."))
FT_SetDataCharacteristics(handle, DWORD(wordlength), DWORD(stopbits), DWORD(parity))
end

"""
    timeouts(d::D2XXDevice, timeout_rd, timeout_wr)

Set the timeouts of an open [`D2XXDevice`](@ref) using [`FT_SetTimeouts`](@ref).
"""
timeouts(d::D2XXDevice, timeout_rd, timeout_wr) =
  timeouts(fthandle(d) , timeout_rd, timeout_wr)

"""
    timeouts(handle::FT_HANDLE, timeout_rd, timeout_wr)

Set the timeouts of an open [`FT_HANDLE`](@ref) using [`FT_SetTimeouts`](@ref).

Behaviour is undefined for `timeout_rd` = 0 and `timeout_wr` = 0.
`timeout_rd` = 0 appears to cause [`read`](@ref) and [`readavailable`](@ref) to
block.

See also: [`isopen`](@ref), [`open`](@ref)
"""
function timeouts(handle::FT_HANDLE, timeout_rd, timeout_wr)
  0 <= timeout_rd || throw(DomainError("0 <= timeout_rd"))
  0 <= timeout_wr || throw(DomainError("0 <= timeout_wr"))
  isopen(handle) || throw(D2XXException("Device must be open to set timeouts."))
  FT_SetTimeouts(handle, timeout_rd, timeout_wr)
end

"""
    reset(d::D2XXDevice)

Reset an [`D2XXDevice`](@ref) using [`FT_ResetDevice`](@ref).
"""
resetdevice(d::D2XXDevice) = resetdevice(fthandle(d))

"""
    reset(handle::FT_HANDLE)

Reset an [`FT_HANDLE`](@ref) using [`FT_ResetDevice`](@ref).
"""
function resetdevice(handle::FT_HANDLE)
  if isopen(handle)
    FT_ResetDevice(handle)
  end
  return
end

"""
    usbparameters(d::D2XXDevice, transfersize_in, transfersize_out)

Set the USB input and output transfer buffer sizes in bytes, of an open [`D2XXDevice`](@ref).

Transfer sizes must be set to a multiple of 64 bytes between 64 bytes and 64k bytes.
`transfersize_in` and `transfersize_out` will be automatically rounded to the nearest
valid value. When usbparameters is called, the change comes into effect immediately
and any data that was held in the driver at the time of the change is lost.

Note that, at present, only transfersize_in is supported.
"""
usbparameters(d::D2XXDevice, transfersize_in, transfersize_out) =
  usbparameters(fthandle(d), transfersize_in, transfersize_out)

"""
    usbparameters(handle::FT_HANDLE, transfersize_in, transfersize_out)

Set the USB input and output transfer buffer sizes in bytes, of an open [`FT_HANDLE`](@ref).

Transfer sizes must be set to a multiple of 64 bytes between 64 bytes and 64k bytes.
`transfersize_in` and `transfersize_out` will be automatically rounded to the nearest
valid value. When `usbparameters` is called, the change comes into effect immediately
and any data that was held in the driver at the time of the change is lost.

Note that, at present, only transfersize_in is supported.
"""
function usbparameters(handle::FT_HANDLE, transfersize_in, transfersize_out)
  isopen(handle) || throw(D2XXException("Device must be open to set USB parameters."))
  transfersize_in = round(transfersize_in/64)*64
  transfersize_out = round(transfersize_out/64)*64
  64 <= transfersize_in  || throw(DomainError("transfersize_in < 64"))
  64 <= transfersize_out || throw(DomainError("transfersize_out < 64"))
  65536 >= transfersize_in  || throw(DomainError("transfersize_in > 65536"))
  65536 >= transfersize_out || throw(DomainError("transfersize_out > 65536"))
  FT_SetUSBParameters(handle, transfersize_in, transfersize_out)
end

"""
    characters(d::D2XXDevice, transfersize_in, transfersize_out)

Set the event- and error characters of an open [`D2XXDevice`](@ref) using [`FT_SetChars`](@ref).

This function allows for inserting specified characters in the data stream to
represent events firing or errors occurring.
"""
characters(d::D2XXDevice, event_char, event_enable, error_char, error_enable) =
  characters(fthandle(d), event_char, event_enable, error_char, error_enable)

"""
    characters(handle::FT_HANDLE, transfersize_in, transfersize_out)

Set the event- and error characters of an open [`FT_HANDLE`](@ref) using [`FT_SetChars`](@ref).

This function allows for inserting specified characters in the data stream to
represent events firing or errors occurring.
"""
function characters(handle::FT_HANDLE, event_char, event_enable, error_char, error_enable)
  isopen(handle) || throw(D2XXException("Device must be open to set event- and error characters."))
  FT_SetChars(handle, UInt8(event_char), UInt8(event_enable), UInt8(error_char), UInt8(error_enable))
end

"""
    latencytimer(d::D2XXDevice, timer_val)

Set the latency timer in milliseconds of an open [`D2XXDevice`](@ref) using [`FT_SetLatencyTimer`](@ref).

In the FT8U232AM and FT8U245AM devices, the receive buffer timeout that is used to flush
remaining data from the receive buffer was fixed at 16 ms. In all other FTDI devices, this
timeout is programmable and can be set in 1 ms intervals between 2ms and 255 ms.
This allows the device to be better optimized for protocols requiring faster response
times for short data packets.
"""
latencytimer(d::D2XXDevice, timer_val) = latencytimer(fthandle(d), timer_val)

"""
    latencytimer(handle::FT_HANDLE, timer_val)

Set the latency timer in milliseconds of an open [`FT_HANDLE`](@ref) using [`FT_SetLatencyTimer`](@ref).

In the FT8U232AM and FT8U245AM devices, the receive buffer timeout that is used to flush
remaining data from the receive buffer was fixed at 16 ms. In all other FTDI devices, this
timeout is programmable and can be set in 1 ms intervals between 2ms and 255 ms.
This allows the device to be better optimized for protocols requiring faster response
times for short data packets.
"""
function latencytimer(handle::FT_HANDLE, timer_val)
  isopen(handle) || throw(D2XXException("Device must be open to set latency timer."))
  isinteger(timer_val) || throw(DomainError("latency timer value not an integer"))
  2 <= timer_val  || throw(DomainError("latency timer < 2"))
  255 >= timer_val  || throw(DomainError("latency timer > 255"))
  FT_SetLatencyTimer(handle, UInt8(timer_val))
end

"""
    latencytimer(d::D2XXDevice)

Return the latency timer in milliseconds of an open [`D2XXDevice`](@ref) using [`FT_GetLatencyTimer`](@ref).
"""
latencytimer(d::D2XXDevice) = latencytimer(fthandle(d))

"""
    latencytimer(handle::FT_HANDLE)

Return the latency timer in milliseconds of an open [`FT_HANDLE`](@ref) using [`FT_GetLatencyTimer`](@ref).
"""
function latencytimer(handle::FT_HANDLE)
  isopen(handle) || throw(D2XXException("Device must be open to get latency timer."))
  return FT_GetLatencyTimer(handle)
end

"""
    status(d::D2XXDevice) ->
      mflaglist::Dict{String, Bool}, lflaglist::Dict{String, Bool}

Return `Bool` dictionaries of flags for  the modem status (`mflaglist`) and
line status (`lflaglist`) for an open [`D2XXDevice`](@ref) using
[`FT_GetModemStatus`](@ref).
"""
status(d::D2XXDevice) = status(fthandle(d))

"""
    status(d::D2XXDevice) ->
      mflaglist::Dict{String, Bool}, lflaglist::Dict{String, Bool}

Return `Bool` dictionaries of flags for  the modem status (`mflaglist`) and
line status (`lflaglist`) for an open [`FT_HANDLE`](@ref) using
[`FT_GetModemStatus`](@ref).

See also: [`isopen`](@ref), [`open`](@ref)
"""
function status(handle::FT_HANDLE)
  isopen(handle) || throw(D2XXException("Device must be open to check status."))
  flags = FT_GetModemStatus(handle)
  modemstatus = flags & 0xFF
  linestatus = (flags >> 8) & 0xFF
  mflaglist = Dict{String, Bool}()
  lflaglist = Dict{String, Bool}()
  mflaglist["CTS"]  = (modemstatus & 0x10) == 0x10
  mflaglist["DSR"]  = (modemstatus & 0x20) == 0x20
  mflaglist["RI"]   = (modemstatus & 0x40) == 0x40
  mflaglist["DCD"]  = (modemstatus & 0x80) == 0x89
  # Below is only non-zero for windows
  lflaglist["OE"]   = (linestatus  & 0x02) == 0x02
  lflaglist["PE"]   = (linestatus  & 0x04) == 0x04
  lflaglist["FE"]   = (linestatus  & 0x08) == 0x08
  lflaglist["BI"]   = (linestatus  & 0x10) == 0x10
  mflaglist, lflaglist
end

if Sys.iswindows()

  """
      driverversion(d::D2XXDevice)

  Get the driver version for an open [`D2XXDevice`](@ref) using
  [`FT_GetDriverVersion`](@ref). Windows only.
  """
  driverversion(d::D2XXDevice) = driverversion(fthandle(d))

  """
      driverversion(handle::FT_HANDLE)

  Get the driver version for an open [`FT_HANDLE`](@ref) using
  [`FT_GetDriverVersion`](@ref). Windows only.

  See also: [`isopen`](@ref), [`open`](@ref)
  """
  function driverversion(handle::FT_HANDLE)
    isopen(handle) || throw(D2XXException("Device must be open to check driver version"))
    version = FT_GetDriverVersion(handle)
    @assert (version >> 24) & 0xFF == 0x00 # 4th byte should be 0 according to docs
    Util.versionnumber(version)
  end

end # Sys.iswindows()

"""
    bitmode(d::D2XXDevice, direction, mode::FTBitMode)

Set the initial pin direction and bit [`FTBitMode`](@ref) for an open [`D2XXDevice`](@ref).
"""
bitmode(d::D2XXDevice, direction, mode::FTBitMode) = bitmode(fthandle(d), direction, mode)

"""
    bitmode(handle::FT_HANDLE, direction, mode::FTBitMode)

Set the initial pin direction and bit [`FTBitMode`](@ref) for an open [`D2XXDevice`](@ref).
"""
function bitmode(handle::FT_HANDLE, direction, mode::FTBitMode)
  isopen(handle) || throw(D2XXException("Device must be open to set bit mode."))
  FT_SetBitMode(handle, UInt8(direction), UInt8(mode))
  return
end

"""
    flush(d::D2XXDevice)

Clear the transmit and receive buffers for an open [`D2XXDevice`](@ref).
"""
Base.flush(d::D2XXDevice) = flush(fthandle(d))

"""
    flush(handle::FT_HANDLE)

Clear the transmit and receive buffers for an open [`FT_HANDLE`](@ref).

See also: [`isopen`](@ref), [`open`](@ref), [`bytesavailable`](@ref)
"""
function Base.flush(handle::FT_HANDLE)
  isopen(handle) || throw(D2XXException("Device must be open to flush."))
  FT_StopInTask(handle)
  FT_Purge(handle, FT_PURGE_TX|FT_PURGE_RX)
  readavailable(handle)
  FT_RestartInTask(handle)
end

# Other Functions
#
function createdeviceinfolist()
  numdevs = FT_CreateDeviceInfoList()
end

function getdeviceinfodetail(deviceidx)
  0 <= deviceidx || throw(DomainError("0 <= deviceidx"))
  deviceidx < createdeviceinfolist() || throw(D2XXException("Device index $deviceidx not in range."))
  idx, flags, typ, id, locid, serialnumber, description, fthandle = FT_GetDeviceInfoDetail(deviceidx)
end

# Library Functions
#

if Sys.iswindows()

  """
      libversion()

  Get a version number from a call to [`FT_GetLibraryVersion`](@ref). Windows
  only.
  """
  function libversion()
    version = FT_GetLibraryVersion()
    @assert (version >> 24) & 0xFF == 0x00 # 4th byte should be 0 according to docs
    Util.versionnumber(version)
  end

end # Sys.iswindows()

# D2XXDevice Accessor Functions
#
"""
    deviceidx(d::D2XXDevice)

Get D2XXDevice index.

See also: [`D2XXDevice`](@ref)
"""
deviceidx(d::D2XXDevice) = d.idx

"""
    deviceflags(d::D2XXDevice)

Get the D2XXDevice flags list.

See also: [`D2XXDevice`](@ref)
"""
deviceflags(d::D2XXDevice) = d.flags

"""
    devicetype(d::D2XXDevice)

Get the D2XXDevice device type.

See also: [`D2XXDevice`](@ref)
"""
devicetype(d::D2XXDevice) = d.typ

"""
  deviceid(d::D2XXDevice)

Get the D2XXDevice device id.

See also: [`D2XXDevice`](@ref)
"""
deviceid(d::D2XXDevice) = d.id

"""
    locationid(d::D2XXDevice)

Get the D2XXDevice location id. This is zero for windows devices.

See also: [`D2XXDevice`](@ref)
"""
locationid(d::D2XXDevice) = d.locid

"""
    serialnumber(d::D2XXDevice)

Get the D2XXDevice device serial number.

See also: [`D2XXDevice`](@ref)
"""
serialnumber(d::D2XXDevice) = d.serialnumber

"""
    description(d::D2XXDevice)

Get the D2XXDevice device description.

See also: [`D2XXDevice`](@ref)
"""
description(d::D2XXDevice) = d.description

"""
    fthandle(d::D2XXDevice)

Get the D2XXDevice device D2XX handle of type ::FT_HANDLE`.

See also: [`D2XXDevice`](@ref)
"""
fthandle(d::D2XXDevice) = d.fthandle[]

"""
    fthandle(d::D2XXDevice, fthandle::FT_HANDLE)

Set the D2XXDevice device D2XX handle of type ::FT_HANDLE`.

See also: [`D2XXDevice`](@ref)
"""
fthandle(d::D2XXDevice, fthandle::FT_HANDLE) = (d.fthandle[] = fthandle)

end # module LibFTD2XX


text:
# New factor interface, something perhaps like this

export calcFactorResidualTemporary

getFactorOperationalMemoryType(dfg::SolverParams) = CommonConvWrapper
# difficult type piracy case needing both types NoSolverParams and CommonConvWrapper.
getFactorOperationalMemoryType(dfg::NoSolverParams) = CommonConvWrapper

function _getDimensionsPartial(ccw::CommonConvWrapper)
  # @warn "_getDimensionsPartial not ready for use yet"
  ccw.partialDims
end
_getDimensionsPartial(data::GenericFunctionNodeData) = _getCCW(data) |> _getDimensionsPartial
_getDimensionsPartial(fct::DFGFactor) = _getDimensionsPartial(_getCCW(fct))
_getDimensionsPartial(fg::AbstractDFG, lbl::Symbol) = _getDimensionsPartial(getFactor(fg, lbl))

"""
    $SIGNATURES
Get `.factormetadata` for each CPT in CCW for a specific factor in `fg`.
"""
_getFMdThread(ccw::CommonConvWrapper,
              thrid::Int=Threads.threadid()) = ccw.cpt[thrid].factormetadata
#
_getFMdThread(fc::Union{GenericFunctionNodeData,DFGFactor},
              thrid::Int=Threads.threadid()) = _getFMdThread(_getCCW(fc), thrid)
#
_getFMdThread(dfg::AbstractDFG,
              lbl::Symbol,
              thrid::Int=Threads.threadid()) = _getFMdThread(_getCCW(dfg, lbl), thrid)
#

# Helper function to construct CF from a CCW
function CalcFactor(ccwl::CommonConvWrapper;
                    factor = ccwl.usrfnc!,
                    metadata = _getFMdThread(ccwl),
                    _sampleIdx = 0,
                    _measCount = length(ccwl.measurement),
                    _legacyMeas = ccwl.measurement,
                    _legacyParams = ccwl.params,
                    _allowThreads = true,
                    cache = ccwl.dummyCache )
  #
  # FIXME using ccwl.dummyCache is not thread-safe
  CalcFactor( factor,
              metadata,
              _sampleIdx,
              _measCount,
              _legacyMeas,
              _legacyParams,
              _allowThreads,
              cache )
end

"""
    $SIGNATURES

Sample the factor stochastic model `N::Int` times and store the samples in the preallocated `ccw.measurement` container.

DevNotes
- Use in place operations where possible and remember `measurement` is a `::Tuple`.
- TODO only works on `.threadid()==1` at present, see #1094
- Also see, JuliaRobotics/RoME.jl#465
"""
sampleFactor(cf::CalcFactor{<:AbstractFactor}, N::Int=1  ) = [getSample(cf) for _=1:N]

function Base.show(io::IO, x::CalcFactor)
  println(io, )
  printstyled(io, " CalcFactor:\n", color=:blue)
  println(io, "  .factor: ", typeof(x.factor))
end

Base.show(io::IO, ::MIME"text/plain", x::CalcFactor) = show(io, x)

"""
    $SIGNATURES

Function to calculate measurement dimension from factor sampling.

Notes
- Will not work in all situations, but good enough so far.
  - # TODO standardize via domain or manifold definition...??
"""
function calcZDim(cf::CalcFactor{T}) where {T <: AbstractFactor}
  #
  try
    M = getManifold(T)
    return manifold_dimension(M)
  catch
    try
      M = getManifold(cf.factor)
      return manifold_dimension(M)
    catch
      @warn "no method getManifold(::$(string(T))), calcZDim will attempt legacy length(sample) method instead"
    end
  end

  # NOTE try to make sure we get matrix back (not a vector)
  smpls = sampleFactor(cf, 2)[1]
  return length(smpls)
end

calcZDim(ccw::CommonConvWrapper) = calcZDim(CalcFactor(ccw))

calcZDim(cf::CalcFactor{<:GenericMarginal}) = 0

calcZDim(cf::CalcFactor{<:ManifoldPrior}) = manifold_dimension(cf.factor.M)

"""
    $SIGNATURES

Helper function for evaluating factor residual functions, by adding necessary `CalcFactor` wrapper.

Notes
- Factor must already be in a factor graph to work
- Will not yet properly support all multihypo nuances, more a function for testing
- Useful for debugging a factor.

Example
```julia
fg = generateGraph_Kaess()

residual = calcFactorResidual(fg, :x1x2f1, [1.0], [0.0], [0.0])
```

Related

[`calcFactorResidualTemporary`](@ref), [`_evalFactorTemporary!`](@ref), [`evalFactor`](@ref), [`approxConv`](@ref)
"""
calcFactorResidual(dfgfct::DFGFactor, args...; ccw::CommonConvWrapper=IIF._getCCW(dfgfct)) = CalcFactor(ccw)(args...)
calcFactorResidual(dfg::AbstractDFG, fctsym::Symbol, args...) = calcFactorResidual(getFactor(dfg, fctsym), args...)

"""
    $SIGNATURES

Evaluate the residual function for a single sample.

Notes
- Binary factors only at this stage, and `multihypo` does not have to be considered in this test
- Assumes calculation is for a single particle, so `meas::Tuple{Z,other}` is only a single particles value.

Example
```julia
residual = calcFactorResidualTemporary(Pose2Pose2(...), (RoME.Pose2,RoME.Pose2), (z_i,), (x1, x2))
```

Related

[`calcFactorResidual`](@ref), [`CalcResidual`](@ref), [`_evalFactorTemporary!`](@ref), [`approxConv`](@ref), [`_buildGraphByFactorAndTypes!`](@ref)
"""
function calcFactorResidualTemporary( fct::AbstractRelative,
                                      varTypes::Tuple,
                                      measurement,
                                      pts::Tuple;
                                      tfg::AbstractDFG = initfg(),
                                      _blockRecursion::Bool=false,
                                      doTime::Bool=false  )
  #

  # build a new temporary graph
  _, _dfgfct = _buildGraphByFactorAndTypes!(fct, varTypes, pts, dfg=tfg, _blockRecursion=_blockRecursion)

  # get a fresh measurement if needed
  _measurement = if measurement != [] #length(measurement) != 0
    measurement
  else
    # now use the CommonConvWrapper object in `_dfgfct`
    cfo = CalcFactor(_getCCW(_dfgfct))
    sampleFactor(cfo, 1)[1]
  end

  # assume a single sample point is being run
  res = if doTime
    @time res = calcFactorResidual(_dfgfct, _measurement, pts...)
    res
  else
    calcFactorResidual(_dfgfct, _measurement, pts...)
  end
  return res
end

## =============================================================================================
## FactorOperationalMemory helper constructors
## =============================================================================================

function ConvPerThread( X::AbstractVector{P},
                        zDim::Int,
                        factormetadata::FactorMetadata;
                        particleidx::Int=1,
                        activehypo= 1:length(params),
                        perturb=zeros(zDim),
                        res=zeros(zDim),
                        thrid_ = 0  ) where P
  #
  return ConvPerThread{typeof(res), typeof(factormetadata), Any}( thrid_,
                        particleidx,
                        factormetadata,
                        Int[activehypo;],
                        perturb,
                        X,
                        res )
end

function CommonConvWrapper( usrfnc::T,
                            X::AbstractVector{P},
                            zDim::Int,
                            varValsLink::NamedTuple,
                            factormetadata::FactorMetadata;
                            partial::Bool=false,
                            hypotheses::H=nothing,
                            certainhypo=nothing,
                            activehypo= 1:length(varValsLink),
                            nullhypo::Real=0,
                            varidx::Int=1,
                            measurement::AbstractVector=Vector(Vector{Float64}(),),
                            particleidx::Int=1,
                            xDim::Int=size(X,1),
                            partialDims::AbstractVector{<:Integer}=1:length(X),
                            perturb=zeros(zDim),
                            res::AbstractVector{<:Real}=zeros(zDim),
                            threadmodel::Type{<:_AbstractThreadModel}=MultiThreaded,
                            inflation::Real=3.0,
                            vartypes=typeof.(getVariableType.(factormetadata.fullvariables)),
                            gradients=nothing,
                            userCache::CT = nothing ) where {T<:AbstractFactor,P,H,CT}
  #

  return  CommonConvWrapper(usrfnc,
                            xDim,
                            zDim,
                            partial,
                            hypotheses,
                            certainhypo,
                            Float64(nullhypo),
                            varValsLink,
                            varidx,
                            measurement,
                            threadmodel,
                            (i->ConvPerThread(X, zDim,factormetadata, particleidx=particleidx,
                                              activehypo=activehypo,
                                              perturb=perturb, res=res )).(1:Threads.nthreads()),
                            inflation,
                            partialDims,
                            DataType[vartypes...],
                            gradients,
                            userCache )
end

function _resizePointsVector!(vecP::AbstractVector{P}, mkd::ManifoldKernelDensity, N::Int) where P
  #
  pN = length(vecP)
  resize!(vecP, N)
  for j in pN:N
    smp = AMP.sample(mkd, 1)[1]
    # @show j, smp, typeof(smp), typeof(vecP[j])
    vecP[j] = smp[1]
  end

  vecP
end

"""
    $(SIGNATURES)

Prepare the particle arrays `ARR` to be used for approximate convolution.
This function ensures that ARR has te same dimensions among all the parameters.
Function returns with ARR[sfidx] pointing at newly allocated deepcopy of the
existing values in getVal(Xi[.label==solvefor]).

Notes
- Return values `sfidx` is the element in ARR where `Xi.label==solvefor` and
- `maxlen` is length of all (possibly resampled) `ARR` contained particles.
- `Xi` is order sensitive.
- for initialization, solveFor = Nothing.
- `P = getPointType(<:InferenceVariable)`

DevNotes
- FIXME ARR internally should become a NamedTuple
"""
function _prepParamVec( Xi::Vector{<:DFGVariable},
                        solvefor::Union{Nothing, Symbol},
                        N::Int=0;
                        solveKey::Symbol=:default  ) where P
  #
  # FIXME ON FIRE, refactor to new NamedTuple instead
  varParamsAll = Vector{Vector{Any}}()

  LEN = Int[]
  maxlen = N # FIXME see #105
  count = 0
  sfidx = 0

  for xi in Xi
    vecP = getVal(xi, solveKey=solveKey)
    push!(varParamsAll, vecP)
    LEN = length.(varParamsAll)
    maxlen = maximum([N; LEN])
    count += 1
    if xi.label == solvefor
      sfidx = count #xi.index
    end
  end

  # resample variables with too few kernels (manifolds points)
  SAMP = LEN .< maxlen
  for i in 1:count
    if SAMP[i]
      Pr = getBelief(Xi[i], solveKey)
      _resizePointsVector!(varParamsAll[i], Pr, maxlen)
    end
  end

  # TODO --rather define reusable memory for the proposal
  # we are generating a proposal distribution, not direct replacement for existing memory and hence the deepcopy.
  if sfidx > 0
    varParamsAll[sfidx] = deepcopy(varParamsAll[sfidx])
  end

  # get solvefor manifolds
  # FIXME deprecate use of (:null,)
  mani = length(Xi)==0 || sfidx==0 ? (:null,) : getManifold(Xi[sfidx])

  varTypes = typeof.(getVariableType.(Xi)) # previous need to force unstable, ::Vector{DataType}

  tup = tuple(varParamsAll...)
  nms = tuple(getLabel.(Xi)...)
  ntp = NamedTuple{nms,typeof(tup)}(tup)

  # FIXME, forcing maxlen to N results in errors (see test/testVariousNSolveSize.jl) see #105
  # maxlen = N == 0 ? maxlen : N
  return ntp, maxlen, sfidx, mani, varTypes
end

"""
    $SIGNATURES
Internal method to set which dimensions should be used as the decision variables for later numerical optimization.
"""
function _setCCWDecisionDimsConv!(ccwl::Union{CommonConvWrapper{F},
                                              CommonConvWrapper{Mixture{N_,F,S,T}}} ) where {N_,F<:Union{AbstractManifoldMinimize, AbstractRelativeMinimize, AbstractRelativeRoots, AbstractPrior},S,T}
  #

  # NOTE should only be done in the constructor
  ccwl.partialDims = if ccwl.partial
    Int[ccwl.usrfnc!.partial...]
  else
    Int[1:ccwl.xDim...]
  end

  nothing
end

function _updateCPTs!(ccwl, sfidx)
  for thrid in 1:Threads.nthreads()
    cpt_ = ccwl.cpt[thrid]
    cpt_.X = ccwl.params[sfidx]
    # used in ccw functor for AbstractRelativeMinimize
    # TODO JT - Confirm it should be updated here. Testing in _prepCCW
    resize!(cpt_.res, ccwl.zDim)
    fill!(cpt_.res, 0.0)
  end
  nothing
end

function attemptGradientPrep( varTypes, usrfnc, varParamsAll, multihypo, meas_single, _blockRecursion )
  # prepare new cached gradient lambdas (attempt)
  try
    # https://github.com/JuliaRobotics/IncrementalInference.jl/blob/db7ff84225cc848c325e57b5fb9d0d85cb6c79b8/src/DispatchPackedConversions.jl#L46
    # also https://github.com/JuliaRobotics/DistributedFactorGraphs.jl/issues/590#issuecomment-891450762
    # FIXME, suppressing nested gradient propagation on GenericMarginals for the time being, see #1010
    if (!_blockRecursion) && usrfnc isa AbstractRelative && !(usrfnc isa GenericMarginal)
      # take first value from each measurement-tuple-element
      measurement_ = meas_single
      # compensate if no info available during deserialization
      # take the first value from each variable param
      pts_ = map(x->x[1], varParamsAll)
      # FIXME, only using first meas and params values at this time...
      # NOTE, must block recurions here, since FGC uses this function to calculate numerical gradients on a temp fg.
      # assume for now fractional-var in multihypo have same varType
      hypoidxs = _selectHypoVariables(pts_, multihypo)
      gradients = FactorGradientsCached!(usrfnc, tuple(varTypes[hypoidxs]...), measurement_, tuple(pts_[hypoidxs]...), _blockRecursion=true);

      return gradients
    end
  catch e
    @warn "Unable to create measurements and gradients for $usrfnc during prep of CCW, falling back on no-partial information assumption.  Enable ENV[\"JULIA_DEBUG\"] = \"IncrementalInference\" for @debug printing to see the error."
    # rethrow(e)
    @debug(e)
  end
  return nothing
end

"""
    $SIGNATURES

Notes
- Can be called with `length(Xi)==0`
"""
function _prepCCW(Xi::Vector{<:DFGVariable},
                  usrfnc::T;
                  multihypo::Union{Nothing, <:Distributions.Categorical}=nothing,
                  nullhypo::Real=0.0,
                  certainhypo = multihypo !== nothing ? collect(1:length(multihypo.p))[multihypo.p .== 0.0] : collect(1:length(Xi)),
                  inflation::Real=0.0,
                  solveKey::Symbol=:default,
                  threadmodel=MultiThreaded,
                  _blockRecursion::Bool=false,
                  userCache::CT = nothing  ) where {T <: AbstractFactor, CT}
  #
  length(Xi) !== 0 ? nothing : @debug("cannot prep ccw.param list with length(Xi)==0, see DFG #590")

  # TODO check no Anys, see #1321
  _varValsQuick, maxlen, sfidx, mani, varTypes = _prepParamVec( Xi, nothing, 0; solveKey ) # Nothing for init.

  # standard factor metadata
  sflbl = 0 == length(Xi) ? :null : getLabel(Xi[end])
  lbs = getLabel.(Xi)
  fmd = FactorMetadata(Xi, lbs, _varValsQuick, sflbl, nothing)

  # create a temporary CalcFactor object for extracting the first sample
  # TODO, deprecate this:  guess measurement points type
  # MeasType = Vector{Float64} # FIXME use `usrfnc` to get this information instead
  _cf = CalcFactor( usrfnc, fmd, 0, 1, nothing, _varValsQuick, false, userCache)

  # get a measurement sample
  meas_single = sampleFactor(_cf, 1)[1]

  #TODO preallocate measurement?
  measurement = Vector{typeof(meas_single)}()

  # partialDims are sensitive to both which solvefor variable index and whether the factor is partial
  partial = hasfield(T, :partial)
  partialDims = if partial
    Int[usrfnc.partial...]
  else
    Int[]
  end

  # as per struct CommonConvWrapper
  gradients = attemptGradientPrep( varTypes, usrfnc, _varValsQuick, multihypo, meas_single, _blockRecursion )

  # variable Types
  pttypes = getVariableType.(Xi) .|> getPointType
  PointType = 0 < length(pttypes) ? pttypes[1] : Vector{Float64}

  return CommonConvWrapper(
          usrfnc,
          PointType[],
          calcZDim(_cf),
          _varValsQuick,
          fmd;
          partial,
          measurement,
          hypotheses = multihypo,
          certainhypo,
          nullhypo,
          threadmodel,
          inflation,
          partialDims,
          vartypes = varTypes,
          gradients,
          userCache
        )
end

"""
    $(SIGNATURES)

Prepare a common functor computation object `prepareCommonConvWrapper{T}` containing
the user factor functor along with additional variables and information using during
approximate convolution computations.

DevNotes
- TODO consolidate with others, see https://github.com/JuliaRobotics/IncrementalInference.jl/projects/6
"""
function _updateCCW!( F_::Type{<:AbstractRelative},
                      ccwl::CommonConvWrapper{F},
                      Xi::AbstractVector{<:DFGVariable},
                      solvefor::Symbol,
                      N::Integer;
                      needFreshMeasurements::Bool=true,
                      solveKey::Symbol=:default  ) where {F <: AbstractFactor}
  #
  length(Xi) !== 0 ? nothing : @debug("cannot prep ccw.param list with length(Xi)==0, see DFG #590")

  # FIXME, order of fmd ccwl cf are a little weird and should be revised.
  # FIXME maxlen should parrot N (barring multi-/nullhypo issues)
  _varValsQuick, maxlen, sfidx, mani, varTypes = _prepParamVec( Xi, solvefor, N; solveKey)

  # NOTE should be selecting for the correct multihypothesis mode
  ccwl.params = _varValsQuick
  # some better consolidate is needed
  ccwl.vartypes = varTypes
  # FIXME ON FIRE, what happens if this is a partial dimension factor?  See #1246
  ccwl.xDim = getDimension(getVariableType(Xi[sfidx]))
  # TODO maybe refactor new type higher up?

  # setup the partial or complete decision variable dimensions for this ccwl object
  # NOTE perhaps deconv has changed the decision variable list, so placed here during consolidation phase
  # TODO, should this not be part of `prepareCommonConvWrapper` -- only here do we look for .partial
  _setCCWDecisionDimsConv!(ccwl)

  # get factor metadata -- TODO, populate, also see #784
  # TODO consolidate with ccwl??
  # FIXME do not divert Mixture for sampling
  cf = CalcFactor(ccwl; _allowThreads=true)

  # cache the measurement dimension
  @assert ccwl.zDim == calcZDim( cf ) "refactoring in progress, cannot drop assignment ccwl.zDim:$(ccwl.zDim) = calcZDim( cf ):$(calcZDim( cf ))"
  # ccwl.zDim = calcZDim( cf )  # CalcFactor(ccwl) )

  # set the 'solvefor' variable index -- i.e. which connected variable of the factor is being computed in this convolution.
  ccwl.varidx = sfidx

  # option to disable fresh samples
  if needFreshMeasurements
    # TODO refactor
    ccwl.measurement = sampleFactor(cf, maxlen)
  end

  # set each CPT
  _updateCPTs!(ccwl, sfidx)

  # calculate new gradients perhaps
  # J = ccwl.gradients(measurement..., pts...)

  return sfidx, maxlen, mani
end

function _updateCCW!( ccwl::Union{CommonConvWrapper{F},
                                                CommonConvWrapper{Mixture{N_,F,S,T}}},
                      Xi::AbstractVector{<:DFGVariable},
                      solvefor::Symbol,
                      N::Integer;
                      kw...  ) where {N_,F<:AbstractRelative,S,T}
  #
  _updateCCW!(F, ccwl, Xi, solvefor, N; kw...)
end

#


text:
struct Vertex{T<:Number, U<:Complex}
  position::AbstractVector{T}
  value::U
end
value(v::Vertex) = v.value
position(v::Vertex) = v.position
import Base.isequal
function Base.isequal(a::Vertex, b::Vertex)
  values_equal = value(a) == value(b) || (isnan(value(a)) && isnan(value(b)))
  positions_equal = all(position(a) .== position(b))
  return values_equal && positions_equal
end
import Base.hash
Base.hash(v::Vertex) = hash(v, hash(:Vertex))
Base.hash(v::Vertex, h::UInt64) = hash(hash.(v.position), hash(v.value, h))


text:
# TODO: https://github.com/JuliaLang/julia/issues/29420
const DEFAULT_EXTRA_DEFINITIONS = Dict{Symbol,AbstractJuliaType}(
    :ssize_t => JuliaCssize_t(),
    :size_t => JuliaCsize_t(),
    :ptrdiff_t => JuliaCptrdiff_t(),
    :uint64_t => JuliaCuint64_t(),
    :uint32_t => JuliaCuint32_t(),
    :uint16_t => JuliaCuint16_t(),
    :uint8_t => JuliaCuint8_t(),
    :int64_t => JuliaCint64_t(),
    :int32_t => JuliaCint32_t(),
    :int16_t => JuliaCint16_t(),
    :int8_t => JuliaCint8_t(),
    :uintptr_t => JuliaCuintptr_t(),
    :tm => JuliaCtm(),
    :wchar_t => JuliaCwchar_t(),
    :va_list => JuliaUnsupported(),
    :FILE => JuliaCFILE(),
    :__uint64_t => JuliaCuint64_t(),
    :__uint32_t => JuliaCuint32_t(),
    :__uint16_t => JuliaCuint16_t(),
    :__uint8_t => JuliaCuint8_t(),
    :__int64_t => JuliaCint64_t(),
    :__int32_t => JuliaCint32_t(),
    :__int16_t => JuliaCint16_t(),
    :__int8_t => JuliaCint8_t(),
)

const EXTRA_DEFINITIONS = Dict{Symbol,AbstractJuliaType}()

add_definition(x::Dict{Symbol,<:AbstractJuliaType}) = merge!(EXTRA_DEFINITIONS, x)
add_definition(x::Pair{Symbol,<:AbstractJuliaType}) = push!(EXTRA_DEFINITIONS, x)

get_definition() = EXTRA_DEFINITIONS

function reset_definition()
    empty!(EXTRA_DEFINITIONS)
    merge!(EXTRA_DEFINITIONS, DEFAULT_EXTRA_DEFINITIONS)
end

macro add_def(c_symbol::Symbol, type::Symbol=:AbstractJuliaSIT, typename::Symbol=c_symbol, translated::Symbol=c_symbol)
    arg1 = Expr(:(::), typename)
    arg2 = Expr(:kw, :options, :(Dict()))
    sig = Expr(:call, Expr(:., :Generators, QuoteNode(:translate)), arg1, arg2)
    func = Expr(:(=), sig, Expr(:block, QuoteNode(:($translated))))
    ret = quote
        struct $typename <: $type end
        add_definition($(QuoteNode(c_symbol)) => $(esc(typename))())
        $(esc(func))
    end
    return ret
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig <: SwaggerModel
    burst::Any # spec type: Union{ Nothing, Int64 } # spec name: burst
    qps::Any # spec type: Union{ Nothing, Int64 } # spec name: qps

    function IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig(;burst=nothing, qps=nothing)
        o = new()
        validate_property(IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig, Symbol("burst"), burst)
        setfield!(o, Symbol("burst"), burst)
        validate_property(IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig, Symbol("qps"), qps)
        setfield!(o, Symbol("qps"), qps)
        o
    end
end # type IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig

const _property_map_IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig = Dict{Symbol,Symbol}(Symbol("burst")=>Symbol("burst"), Symbol("qps")=>Symbol("qps"))
const _property_types_IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig = Dict{Symbol,String}(Symbol("burst")=>"Int64", Symbol("qps")=>"Int64")
Base.propertynames(::Type{ IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig }) = collect(keys(_property_map_IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig))
Swagger.property_type(::Type{ IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig[name]))}
Swagger.field_name(::Type{ IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig }, property_name::Symbol) =  _property_map_IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig[property_name]

function check_required(o::IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig)
    true
end

function validate_property(::Type{ IoK8sApiAuditregistrationV1alpha1WebhookThrottleConfig }, name::Symbol, val)
end


text:
include("common.jl")

using .common
using Plots
using LinearAlgebra
using ProgressMeter
gr()
theme(:solarized_light)

dt = 0.00001
m1 = 0.6
m2 = 1 - m1

function rk4!(p::Particle, a, h)
    k1_v = a
    k1_r = p.v

    k2_v = a .* (k1_r*h/2)
    k2_r = p.v .* k1_v*h/2

    k3_v = a .* (k2_r*h/2)
    k3_r = p.v .* k2_v*h/2

    k4_v = a .* (k3_r*h)
    k4_r = p.v .* k3_v*h

    p.r += h/6 * (k1_r + 2*k2_r + 2*k3_r + k4_r)
    p.v += h/6 * (k1_v + 2*k2_v + 2*k3_v + k4_v)
end

p1 = Particle([1., 0., 0.], [0, -1., 0], 1)

# initialize a 3D plot with 1 empty series
err_plt = plot(1, title= "Error")
plt = plot(1,  ylim=(-1.2, 1.2), xlim=(-1.2,1.2),
                title = "PJP <3", size = (800, 800),
                m = (:cross, 3, stroke(0)))

function step!(p::Particle)
    len_r2 = dot(p.r, p.r)
    len_r3 = len_r2 * sqrt(len_r2)
    a = - (p.r / len_r3);

    rk4!(p, a, dt)
end

n = 45
prog = Progress(n,1)
# build an animated gif by pushing new points to the plot, saving every 10th frame
for i=1:n
    for j=1:100000
        step!(p1)
    end
    push!(err_plt, dot(p1.r, p1.r) - 1)

    push!(plt, p1.r[1], p1.r[2])
    next!(prog)
end

plot(plt, err_plt, layout= grid(2, 1, heights=[0.7,0.3]), size = (500, 700))


text:
"Returns true if the first three elements of `names` are `(:u, :v, :w)`."
has_velocities(names) = :u == names[1] && :v == names[2] && :w == names[3]

# Tuples of length 0-2 cannot contain velocity fields
has_velocities(::Tuple{}) = false
has_velocities(::Tuple{X}) where X = false
has_velocities(::Tuple{X, Y}) where {X, Y} = false

tracernames(::Nothing) = ()
tracernames(name::Symbol) = tuple(name)
tracernames(names::NTuple{N, Symbol}) where N = has_velocities(names) ? names[4:end] : names
tracernames(::NamedTuple{names}) where names = tracernames(names)


text:
function retrieve_length(spec)
    chain = length_chain(spec, spec.len)
    @match length(chain) begin
        1 => vk_call(first(chain))
        GuardBy(>(1)) =>
            chain_getproperty(:($(wrap_identifier(first(chain))).vks), getproperty.(chain[2:end], :name))
    end
end

chain_getproperty(ex, props) = foldl((x, y) -> :($x.$y), props; init = ex)

function initialize_ptr(param::SpecFuncParam, count_ptr::SpecFuncParam)
    @assert is_data(param) && is_size(count_ptr)
    :($(param.name) = Libc.malloc($(count_ptr.name)[]))
end

function initialize_ptr(param::SpecFuncParam)
    rhs = @match param begin
        GuardBy(is_arr) => :(Vector{$(ptr_type(param.type))}(undef, $(retrieve_length(param))))
        _ => @match param.type begin
            :(Ptr{Cvoid}) => :(Ref{Ptr{Cvoid}}())
            _ => :(Ref{$(ptr_type(param.type))}())
        end
    end
    :($(param.name) = $rhs)
end

function initialize_array(param::SpecFuncParam, count_ptr::SpecFuncParam)
    :($(param.name) = Vector{$(ptr_type(param.type))}(undef, $(count_ptr.name)[]))
end


text:
include("../imports.jl")

struct ConstGrouper <: FastAI.ParamGrouper
    g
end
FastAI.group(cg::ConstGrouper, m) = Dict(cg.g => m)

@testset ExtendedTestSet "`finetune!`" begin
    learner = testlearner(Recorder())
    @test_nowarn finetune!(learner, 1; grouper = ConstGrouper(2))

end


text:
using PyPlot
using Aerospace

# Motor Initialization
wn = 20.0;
zeta = 0.707;
rpm_pos_lim = 20.0;
rpm_neg_lim = 0.0;
rpm_rate_pos_lim = 15.0;
rpm_rate_neg_lim = -15.0;

motor1 = Initialize_Act(dt, wn, zeta,
                        rpm_pos_lim,
                        rpm_neg_lim,
                        rpm_rate_pos_lim,
                        rpm_rate_neg_lim
                        );

# Simulation Control
time = 0.0;
dt = 0.01;
tfinal = 5.0;

## Number of Sample points
nn = (tfinal - time)/dt + 1;
NN = Int64(nn);

TT = zeros(NN);
UU = zeros(NN);
X1 = zeros(NN);
X2 = zeros(NN);

for i = 1:NN
    TT[i] = time;
    X1[i] = motor1.position;
    X2[i] = motor1.rate;

    ucmd = 15.0;

    motor1 = Update_Act(motor1, ucmd);
    time = time + dt;

end

figure()
plot(TT, X1, TT, X2)
grid()
xlabel("Time (s)")
ylabel("Actuator Output")
legend(["Position", "Rate"])


text:
struct ConvergenceStatus
    has_converged::Bool
    n_iterations::Int
    norm_residuals::Float64
    n_linear_solves::Int
end

mutable struct Spmat{VTI<:AbstractVector, VTF<:AbstractVector}
    colptr::VTI
    rowval::VTI
    nzval::VTF

    # create 2 Spmats from complex matrix
    function Spmat{VTI, VTF}(mat::SparseMatrixCSC{Complex{Float64}, Int}) where {VTI, VTF}
        matreal = new(VTI(mat.colptr), VTI(mat.rowval), VTF(real.(mat.nzval)))
        matimag = new(VTI(mat.colptr), VTI(mat.rowval), VTF(imag.(mat.nzval)))
        return matreal, matimag
    end
    # copy constructor
    function Spmat{VTI, VTF}(mat) where {VTI, VTF}
        return new(VTI(mat.colptr), VTI(mat.rowval), VTF(mat.nzval))
    end
end

# norm
xnorm(x::AbstractVector) = norm(x, 2)
xnorm(x::CuVector) = CUBLAS.nrm2(x)

# Array initialization
xzeros(S, n) = fill!(S(undef, n), zero(eltype(S)))
xones(S, n) = fill!(S(undef, n), one(eltype(S)))

# projection operator
function project!(w::VT, u::VT, u♭::VT, u♯::VT) where VT<:AbstractArray
    w .= max.(min.(u, u♯), u♭)
end

# Utils function to solve transposed linear system  A' x = y
# Source code taken from:
# https://github.com/JuliaGPU/CUDA.jl/blob/master/lib/cusolver/wrappers.jl#L78L111
function csclsvqr!(A::CuSparseMatrixCSC{Float64},
                    b::CuVector{Float64},
                    x::CuVector{Float64},
                    tol::Float64,
                    reorder::Cint,
                    inda::Char)
    n = size(A,1)
    desca = CUSPARSE.CuMatrixDescriptor(
        CUSPARSE.CUSPARSE_MATRIX_TYPE_GENERAL,
        CUSPARSE.CUSPARSE_FILL_MODE_LOWER,
        CUSPARSE.CUSPARSE_DIAG_TYPE_NON_UNIT, inda)
    singularity = Ref{Cint}(1)
    CUSOLVER.cusolverSpDcsrlsvqr(CUSOLVER.sparse_handle(), n, A.nnz, desca, A.nzVal, A.colPtr, A.rowVal, b, tol, reorder, x, singularity)

    if singularity[] != -1
        throw(SingularException(singularity[]))
    end

    x
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_index

"""
    ''_index(;kwargs...)

An index component.
DateTimePicker is a datetime input component.
The inputs can be initialized with the `defaultValue` property and the
input date, on ISO format, is specified with the `value` property.
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `value` (String; optional): The input's current value
- `defaultValue` (String; optional): The input's default value. Accepts values on the ISO format.
- `renderTimezone` (Bool; optional): Defines if the timezone input should be rendered. Defaults to true.
- `style` (Dict | String; optional): The container's style.
- `dateInputStyle` (Dict | String; optional): The date input's style.
- `timezoneInputStyle` (Dict; optional): The timezone input's container style.
"""
function ''_index(; kwargs...)
        available_props = Symbol[:id, :value, :defaultValue, :renderTimezone, :style, :dateInputStyle, :timezoneInputStyle]
        wild_props = Symbol[]
        return Component("''_index", "index", "dash_cool_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: raulrozza/dash-cool-components

max_issues_repo_name: raulrozza/dash-cool-components

max_forks_repo_name: raulrozza/dash-cool-components


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module PilotResponse

export
    AbstractPilotResponse,

    initialize,
    update,

    updatePilotResponse,

    SimplePilotResponse,
    SimplePRResolutionAdvisory,
    SimplePRCommand,

    StochasticLinearPR,
    StochasticLinearPRCommand,
    StochasticLinearPRRA,

    LLDetPR,
    LLDetPRAircraft,
    LLDetPRCommand,
    LLDetPRRA,
    LLDetPROutput

import CommonInterfaces: initialize, update

using AbstractPilotResponseImpl

using SimplePilotResponseImpl

using StochasticLinearPRImpl

using LLDetPRImpl

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES