text:
"""
Compose to Image Reparameterizations

    apply_gam_to_gam(gamnew::Array, gam::Array)
    :param gamnew: 3-D Array describing new gamma
    :param gam:: 3-D Array describing current gamma
"""
function apply_gam_to_gam(gamnew::Array, gam::Array)
    # gam \circ gam0
    m, n, D = size(gam);
    md = 8;
    mt = md*m;
    nt = md*n;
    U = LinRange(0,1,m);
    V = LinRange(0,1,n);
    Ut = LinRange(0,1,mt);
    Vt = LinRange(0,1,nt);
    gam_tmp = zeros(mt, nt, D);
    gam_new_tmp = zeros(mt, nt, D);
    for i = 1:D
        spl = Spline2D(U,V,gam[:,:,i],kx=1,ky=1);
        gam_tmp[:,:,i] = evalgrid(spl,Ut,Vt);
        spl = Spline2D(U,V,gamnew[:,:,i],kx=1,ky=1);
        gam_new_tmp[:,:,i] = evalgrid(spl,Ut,Vt);
    end

    gam_cum_tmp = zeros(mt,nt,D);
    for i = 1:D
        spl = Spline2D(Ut,Vt,gam_new_tmp[:,:,i],kx=1,ky=1);
        for j = 1:mt
            gam_cum_tmp[j,:,i] = spl(vec(gam_tmp[j,:,2]), vec(gam_tmp[j,:,1]));
        end
    end

    gam_cum = zeros(m,n,D);
    for i = 1:D
        spl = Spline2D(Ut,Vt,gam_cum_tmp[:,:,i],kx=1,ky=1);
        gam_cum[:,:,i] = evalgrid(spl,U,V);
    end

    return gam_cum
end

"""
Warp img by gam

    apply_gam_to_imag(img::Array, gam::Array)
    :param img: 2-D or 3-D Array
    :param gam: 3-D Array
"""
function apply_gam_to_imag(img::Array, gam::Array)
    if ndims(img) == 3
        m,n,d = size(img);
        img_new = zeros(m,n,d);
    elseif ndims(img) == 2
        m,n = size(img);
        d = 1;
        img_new = zeros(m,n);
    end

    U = LinRange(0,1,m);
    V = LinRange(0,1,n);

    if d == 1
        spl = Spline2D(U,V,img);
        for j = 1:m
            img_new[j,:] = spl(vec(gam[j,:,2]),vec(gam[j,:,1]));
        end
    else
        for i = 1:d
            spl = Spline2D(U,V,img[:,:,i]);
            for j = 1:m
                img_new[j,:,i] = spl(vec(gam[j,:,2]),vec(gam[j,:,1]));
            end
        end
    end

    return img_new
end

"""
Compose Image Reparameterization with Identity

    apply_gam_gamid(gamid::Array,gaminc::Array)
    :param gamid: 3-D Array
    :param gaminc: 3-D Array
"""
function apply_gam_gamid(gamid::Array,gaminc::Array)
    m,n,d = size(gamid);
    U = LinRange(0,1,m);
    V = LinRange(0,1,n);

    gam_cum = zeros(m,n,d);
    for j = 1:d
        spl = Spline2D(U,V,gamid[:,:,j]);
        for i = 1:m
            gam_cum[i,:,j] = spl(vec(gaminc[i,:,2]),vec(gaminc[i,:,1]));
        end
    end

    return gam_cum
end

"""
Find 2-D gradient

    compgrad2D(f::Array)
    :param f: Array

    :return dfdu: Array
    :return dfdv: Array
"""
function compgrad2D(f::Array)
    dims = ndims(f);
    if dims < 3
        n,t = size(f);
        d = 1;
        dfdu = zeros(n,t);
        dfdv = zeros(n,t);
    else
        n,t,d = size(f);
        dfdu = zeros(n,t,d);
        dfdv = zeros(n,t,d);
    end

    @cxx findgrad2D(dfdu, dfdv, f, n, t, d)

    return dfdu, dfdv
end

"""
Compute Gram Schmidt of basis using c library gradient

    gram_schmidt_c(b)
    :param b: basis
"""
function gram_schmidt_c(b)
    m,n,tmp,N = size(b);
    ds = 1./((m-1)*(n-1));

    cnt = 1;
    G = zeros(m,n,tmp,N);
    G[:,:,:,cnt] = b[:,:,:,cnt];

    dvx, dvy = compgrad2D(G[:,:,:,cnt]);
    l = vec(dvx)'*vec(dvx)*ds + vec(dvy)'*vec(dvy)*ds;

    for i = 2:N
        G[:,:,:,i] = b[:,:,:,i];
        dv1x, dv1y = compgrad2D(G[:,:,:,i]);
        for j = 1:i-1
            dv2x, dv2y = compgrad2D(G[:,:,:,j]);
            t = vec(dv1x)'*vec(dv2x)*ds + vec(dv1y)'*vec(dv2y)*ds;
            G[:,:,:,i] = G[:,:,:,i]-t.*G[:,:,:,j];
        end

        v = G[:,:,:,i];
        l = vec(v)'*vec(v)*ds;
        l = l[1];

        if l>0
            cnt += 1;
            G[:,:,:,cnt] = G[:,:,:,cnt]./sqrt(l);
        end
    end

    return G
end

"""
Find jacobian of image F

    jacob_imag(F::Array)
    :param F: Array
"""
function jacob_imag(F::Array)
    m,n,d = size(F);

    dfdu, dfdv = compgrad2D(F);

    mult_factor = zeros(m,n);
    if (d==2)
        mult_factor = dfdu[:,:,1]*dfdv[:,:,2] - dfdu[:,:,2]*dfdv[:,:,1];
        mult_factor = abs(mult_factor);
    elseif (d==3)
        mult_factor = (dfdu[:,:,2]*dfdv[:,:,3] - dfdu[:,:,3]*dfdv[:,:,2]).^2
            + (dfdu[:,:,1]*dfdv[:,:,3] - dfdu[:,:,3]*dfdv[:,:,1]).^2
            + (dfdu[:,:,1]*dfdv[:,:,2] - dfdu[:,:,2]*dfdv[:,:,1]).^2;
        mult_factor = sqrt(mult_factor);
    end

    return mult_factor
end

"""
Form identity image warping

    makediffeoid(nrow::Integer,ncol::Integer)
    :param nrow: number of rows
    :param ncol: number of columns
"""
function makediffeoid(nrow::Integer,ncol::Integer)
    D = 2;
    gamid = zeros(nrow,ncol,D);
    U, V = meshgrid(LinRange(0,1,ncol),LinRange(0,1,nrow));

    gamid[:,:,1] = U;
    gamid[:,:,2] = V;

    return gamid
end

"""
Convert image to q-map

    imag_to_q(F::Array)
    :param F: Array of Image
"""
function imag_to_q(F::Array)
    dims = ndims(F);
    if dims < 3
        error("Data dimension is wrong!")
    end
    d = size(F,3);
    if d < 2
        error("Data dimension is wrong!")
    end
    q = copy(F);

    sqrtmultfact = sqrt(jacob_imag(F));
    for i = 1:d
        q[:,:,i] = sqrtmultfact.*F[:,:,i];
    end

    return q
end

"""
Compute distance between to q-maps

    comp_dist(q1::Array, q2::Array)
    :param q1: q-map1
    :param q2: q-map2
"""
function comp_dist(q1::Array, q2::Array)
    m = size(q1,1);
    n = size(q1,2);
    ds = 1./((m-1)*(n-1));

    tmp = q1-q2;

    H = sum(sqrt(vec(tmp).*vec(tmp))*ds);

    return H
end

"""
Compute energy of gradient

    comp_energy(q1::Array,q2::Array)
    :param q1: q-map1
    :param q2: q-map2
"""
function comp_energy(q1::Array,q2::Array)
    m,n,d = size(q1);
    ds = 1./(m-1)/(n-1);

    tmp = q1-q2;
    H = vec(tmp)'*vec(tmp)*ds;
    H = H[1];

    return H
end

"""
Update warping function

    update_gam(qt, qm, b)
"""
function update_gam(qt, qm, b)
    v = qt - qm;
    w = findphistar(qt,b);

    gamupdate, innp = findupdategam(v,w,b);

    return gamupdate
end

"""
Find phi*

    findphistar(q,b)
"""
function findphistar(q, b)

    n,t,d = size(q);
    K = size(b,4);

    w = zeros(n,t,d,K);

    @cpp ccall((:findphistar, libfdaqmap), Void, (Ptr{Float64}, Ptr{Float64},
               Ptr{Float64}, Int32, Int32, Int32, Int32), w, q, b, n, t, d, K);
    return w
end

"""
Find updated warping function

    findupdategam(v, w, b)
"""
function findupdategam(v, w, b)
    m,n,D,K = size(b);
    ds = 1./((m-1)*(n-1));

    innp = zeros(K);

    gamupdate = zeros(m,n,D);

    for k = 1:K
        vt = w[:,:,:,k];
        tmp = vec(v)'*vec(vt)*ds;
        innp[k] = tmp[1];

        gamupdate = gamupdate + innp[k].*b[:,:,:,k];
    end

    return gamupdate, innp
end

"""
Check if warping function is a diffeo

    check_crossing(F)
"""
function check_crossing(f)
    n,t,D = size(f);

    if D!=2
        error("Third dimension of first argument expected to be 2")
    end

    diffeo = @cpp ccall((:check_crossing, libfdaqmap), Int32, (Ptr{Float64},
                        Int32, Int32, Int32), f, n, t, D);

    if diffeo == 0
        is_diffeo = false;
    elseif diffeo == 1
        is_diffeo = true;
    else
        error("check_crossing error")
    end

    return is_diffeo
end

"""
Compute Basis Tid

    formbasisTid(M,m,n;basis_type="tso")
    :param M: number of basis elements
    :param m:
    :param n:
    :param basis_type: type of basis (options are "t","s","o","i")
"""
function formbasisTid(M,m,n;basis_type="t")
    U, V = meshgrid(LinRange(0,1,n), LinRange(0,1,m));

    idx = 1;

    if basis_type=="t"
        b = zeros(m,n,2,2*M);
        for s = 1:M
            c = sqrt(2)*pi*s;
            sPI2 = 2*pi*s;

            b[:,:,1,idx] = zeros(m,n);
            b[:,:,2,idx] = (cos(sPI2*V)-1)/c;

            b[:,:,1,idx+1] = (cos(sPI2*U)-1)/c;
            b[:,:,2,idx+1] = zeros(m,n);

            idx += 2;
        end
    end

    if basis_type == "s"
        b = zeros(m,n,2,2*M);
        for s = 1:M
            c = sqrt(2)*pi*s;
            sPI2 = 2*pi*2;

            b[:,:,1,idx] = zeros(m,n);
            b[:,:,2,idx] = sin(sPI2*U)/c;

            b[:,:,1,idx+1] = sin(sPI2*U)/c;
            b[:,:,2,idx+1] = zero(m,n);

            idx += 2;
        end
    end

    if basis_type == "i"
        b = zeros(m,n,2,M*M*8);
        for s1 = 1:M
            s1PI2 = 2*pi*s1;
            for s2 = 1:M
                s2PI2 = 2*pi*s2;
                c = pi * sqrt(s1^2+3*s2^2);
                b[:,:,1,idx] = (cos(s1PI2*U)-1).*(cos(s2PI2*V))/c;
                b[:,:,2,idx] = zeros(m,n);
                b[:,:,1,idx+2] = ((cos(s1PI2*U)-1).*sin(s2PI2*V))/c;
                b[:,:,2,idx+2] = zeros(m,n);
                c = pi*sqrt(s1^2+s2^2);
                b[:,:,1,idx+4] = sin(s1PI2*U).*(cos(s2PI2*V))/c;
                b[:,:,2,idx+4] = zeros(m,n);
                b[:,:,1,idx+6] = (sin(s1PI2*U).*sin(s2PI2*V))/c;
                b[:,:,2,idx+6] = zeros(m,n);

                c = pi*sqrt(s1^2+3*s2^2);
                b[:,:,1,idx+1] = zeros(m,n);
                b[:,:,2,idx+1] = (cos(s1PI2*V)-1).*(cos(s2PI2*U))/c;
                b[:,:,1,idx+3] = zeros(m,n);
                b[:,:,2,idx+3] = ((cos(s1PI2*V)-1).*sin(s2PI2*U))/c;
                c = pi*sqrt(s1^2+s2^2);
                b[:,:,1,idx+5] = zeros(m,n);
                b[:,:,2,idx+5] = sin(s1PI2*V).*(cos(s2PI2*U))/c;
                b[:,:,1,idx+7] = zeros(m,n);
                b[:,:,2,idx+7] = (sin(s1PI2*V).*sin(s2PI2*U))/c;
                idx += 8;
            end
        end
    end

    if basis_type == "o"
        b = zeros(m,n,2,M*4);
        for s = 1:M
            c = sqrt((4*pi^2*s^2+9)/6);
            sPI2 = 2*pi*s;
            b[:,:,1,idx] = (cos(sPI2*U)-1)*V/c;
            b[:,:,2,idx] = zeros(m,n);
            b[:,:,2,idx+1] = (cos(sPI2*V)-1).*U/c;
            b[:,:,1,idx+1] = zeros(m,n);
            b[:,:,1,idx+2] = sin(sPI2*U).*V/c;
            b[:,:,2,idx+2] = zeros(m,n);
            b[:,:,2,idx+3] = sin(sPI2*V).*U/c;
            b[:,:,1,idx+3] = zeros(m,n);
            idx += 4;
        end
    end

    return b, U, V
end

"""
Generate basis elements on T_id(\gamma)

    run_basis(Ft, M=10, basis_type="t", is_orthon=true)
    :param Ft: Image Array
    :param M: number of basis elements
    :param basis_type: type of basis (option "t", "s", "i", "o")
    :param is_orthon: make basis orthonormal

    :return b: basis elements
    :retrun gamid: identity diffeo
"""
function run_basis(Ft, M=10, basis_type="t", is_orthon=true)
    m = size(Ft,1);
    n = size(Ft,2);

    gamid = makediffeoid(m,n);

    b, U, V = formbasisTid(M, m, n, basis_type=basis_type);
    if is_orthon
        b = gram_schmidt_c(b);
    end

    return b, gamid
end


text:
# Autogenerated wrapper script for poly2tri_jll for aarch64-linux-musl-cxx11
export libpoly2tri

JLLWrappers.@generate_wrapper_header("poly2tri")
JLLWrappers.@declare_library_product(libpoly2tri, "libpoly2tri.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libpoly2tri,
        "lib/libpoly2tri.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
function plot_network(case, spring_constant=1e-3)
    data = layout_graph_vega(case, spring_constant)
    remove_information!(data)
    df = form_df(data)
    p = @vlplot(
        width=500,
        height=500,
        config={view={stroke=nothing}},
        x={axis=nothing},
        y={axis=nothing},
        color={
            :ComponentType,
            scale={
                domain=[
                    "branch","bus","gen","connector",
                ],
                range=[
                    :green,
                    :blue,
                    :red,
                    :gray
                ]
            }
        },
    ) +
    @vlplot(
        mark ={
            :rule,
            tooltip=("content" => "data"),
            opacity =  1.0
        },
        data=df["branch"],
        x = :xcoord_1,
        x2 = :xcoord_2,
        y = :ycoord_1,
        y2 = :ycoord_2,
        size={value=5},
    ) +
    @vlplot(
        mark ={
            :rule,
            "tooltip" =("content" => "data"),
            opacity =  1.0
        },
        data=df["connector"],
        x = :xcoord_1,
        x2 = :xcoord_2,
        y = :ycoord_1,
        y2 = :ycoord_2,
        size={value=3},
        strokeDash={value=[4,4]}
    ) +
    @vlplot(
        data = df["bus"],
        mark ={
            :circle,
            "tooltip" =("content" => "data"),
            opacity =  1.0
        },
        x={:xcoord_1,},
        y={:ycoord_1,},
        size={value=1e2},
    )+
    @vlplot(
        data = df["gen"],
        mark ={
            :circle,
            "tooltip" =("content" => "data"),
            opacity =  1.0
        },
        x={:xcoord_1,},
        y={:ycoord_1,},
        size={value=5e1},
    )
    return p
end

function plot_power_flow(case;spring_constant=1e-3)

    for (gen_id,gen) in case["gen"]
        if !haskey(gen,"pg")
            @warn "Generator $(gen_id) does not have key `pg`"
        end
    end
    for (branch_id,branch) in case["branch"]
        if !haskey(branch,"pt")
            @warn "Branch $(branch_id) does not have key `pt`"
        end
        if !haskey(branch,"pf")
            @warn "Branch $(branch_id) does not have key `pf`"
        end
    end
    for (bus_id,bus) in case["bus"]
        if !haskey(bus,"vm")
            @warn "Bus $(bus_id) does not have key `vm`"
        end
    end

    data = layout_graph_vega(case,spring_constant)
    remove_information!(data)
    df = form_df(data)
    p = @vlplot(
        width=500,
        height=500,
        config={view={stroke=nothing}},
        x={axis=nothing},
        y={axis=nothing},
        color={"PercentRated:q",scale={"range"= ["black", "black", "red"]},title="Percent of Rated Capacity"}
    ) +
    @vlplot(
        mark ={
            :rule,
            "tooltip" =("content" => "data"),
            opacity =  1.0
        },
        data=df["branch"],
        transform=[
            {
                calculate="abs(datum.pf) /datum.rate_a * 100",
                as="PercentRated"
            }
        ],
        x = :xcoord_1,
        x2 = :xcoord_2,
        y = :ycoord_1,
        y2 = :ycoord_2,
        size={value=5},
    ) +
    @vlplot(
        mark ={
            :rule,
            "tooltip" =("content" => "data"),
            opacity =  1.0
        },
        data=df["connector"],
        x = :xcoord_1,
        x2 = :xcoord_2,
        y = :ycoord_1,
        y2 = :ycoord_2,
        size={value=3},
        strokeDash={value=[4,4]},
        color={value="gray"}
    ) +
    @vlplot(
        data = df["bus"],
        transform=[
            {
                calculate="abs(datum.vm) /( datum.vmax+datum.vmin) * 100",
                as="PercentRated"
            }
        ],
        mark ={
            :point,
            "tooltip" =("content" => "data"),
            opacity =  1.0,
            filled=true
        },
        x={:xcoord_1,},
        y={:ycoord_1,},
        size={value=2e2},
        shape="ComponentType"
    )+
    @vlplot(
        data = df["gen"],
        transform=[
            {
                calculate="abs(datum.pg) /datum.pmax * 100",
                as="PercentRated"
            }
        ],
        mark ={
            :point,
            "tooltip" =("content" => "data"),
            opacity =  1.0,
            filled=true
        },
        x={:xcoord_1,},
        y={:ycoord_1,},
        size={value=5e1},
        shape="ComponentType"
    )
    return p
end

function remove_information!(data)
    invalid_keys = Dict("branch"  => ["mu_angmin", "mu_angmax", "mu_sf", "shift", "rate_b", "rate_c", "g_to", "g_fr", "mu_st", "source_id", "f_bus", "br_status", "t_bus",  "qf", "angmin", "angmax", "qt", "transformer", "tap"],#["b_fr","b_to", "xcoord_1", "xcoord_2", "ycoord_1", "ycoord_2", "pf", "src","dst","rate_a","br_r","br_x","index"],
                        "bus"     => ["mu_vmax", "lam_q", "mu_vmin", "source_id", "area","lam_p","zone", "bus_i"],#["xcoord_1", "ycoord_1", "bus_type", "name", "vmax",  "vmin", "index", "va", "vm", "base_kv"],
                        "gen"     => ["gen_status","vg","gen_bus","cost","ncost", "qc1max","qc2max", "ramp_agc", "qc1min", "qc2min", "pc1", "ramp_q", "mu_qmax", "ramp_30", "mu_qmin","model", "shutdown", "startup","ramp_10","source_id", "mu_pmax", "pc2", "mu_pmin","apf",],#["xcoord_1", "ycoord_1",  "pg", "qg",  "pmax",   "mbase", "index", "cost", "qmax",  "qmin", "pmin", ]
    )
    for comp_type in ["bus","branch","gen"]
        for (id, comp) in data[comp_type]
            for key in keys(comp)
                if (key in invalid_keys[comp_type])
                    delete!(comp,key)
                end
            end
        end
    end
end

function plot_power_flow_geo(case; spring_constant = 1e-3)
    for (gen_id,gen) in case["gen"]
        if !haskey(gen,"pg")
            @warn "Generator $(gen_id) does not have key `pg`"
        end
    end
    for (branch_id,branch) in case["branch"]
        if !haskey(branch,"pt")
            @warn "Branch $(branch_id) does not have key `pt`"
        end
        if !haskey(branch,"pf")
            @warn "Branch $(branch_id) does not have key `pf`"
        end
    end
    for (bus_id,bus) in case["bus"]
        if !haskey(bus,"vm")
            @warn "Bus $(bus_id) does not have key `vm`"
        end
    end

    data = layout_graph_vega(case,spring_constant)
    remove_information!(data)
    df = form_df(data)
    p = @vlplot(
        config={view={stroke=nothing}},
        x={axis=nothing},
        y={axis=nothing},
        projection={type=:albersUsa},
        color={"PercentRated:q",scale={"range"= ["black", "black", "red"]},title="Percent of Rated Capacity"}
    ) +
    @vlplot(
        mark ={
            :rule,
            "tooltip" =("content" => "data"),
            opacity =  1.0
        },
        data=df["branch"],
        transform=[
            {
                calculate="abs(datum.pf) /datum.rate_a * 100",
                as="PercentRated"
            }
        ],
        longitude="ycoord_1:q",
        latitude="xcoord_1:q",
        longitude2="ycoord_2:q",
        latitude2="xcoord_2:q",
        size={value=3},
    ) +
    @vlplot(
        mark ={
            :rule,
            "tooltip" =("content" => "data"),
            opacity =  1.0
        },
        data=df["connector"],
        longitude="ycoord_1:q",
        latitude="xcoord_1:q",
        longitude2="ycoord_2:q",
        latitude2="xcoord_2:q",
        size={value=2},
        strokeDash={value=[4,4]},
        color={value="gray"}
    ) +
    @vlplot(
        data = df["bus"],
        transform=[
            {
                calculate="abs(datum.vm) /( datum.vmax+datum.vmin) * 100",
                as="PercentRated"
            }
        ],
        mark ={
            :point,
            tooltip =("content" => "data"),
            opacity =  1.0,
            filled=true
        },
        longitude="ycoord_1:q",
        latitude="xcoord_1:q",
        size={value=1e2},
        shape="ComponentType"
    ) +
    @vlplot(
        data = df["gen"],
        transform=[
            {
                calculate="abs(datum.pg) /datum.pmax * 100",
                as="PercentRated"
            }
        ],
        mark ={
            :point,
            "tooltip" =("content" => "data"),
            opacity =  1.0,
            filled=true
        },
        longitude="ycoord_1:q",
        latitude="xcoord_1:q",
        size={value=2e1},
        shape="ComponentType"
    )
    return p
end

function plot_WI(;size=(800,600))
    # download from  https://github.com/deldersveld/topojson/blob/master/countries/us-states/WI-55-wisconsin-counties.json
    file_path = "$(joinpath(dirname(pathof(PowerVega)), ".."))/test/networks/WI-55-wisconsin-counties.json"
    WI = VegaDatasets.VegaJSONDataset(JSON.parsefile(file_path), Path(file_path))

    (width,height) = size
    p =  @vlplot(width=width, height=height) +
     @vlplot(
         mark={
             :geoshape,
             fill=:lightgray,
             stroke=:white
         },
         data={
             values= WI,
             format={
                 type=:topojson,
                 feature=:cb_2015_wisconsin_county_20m
             }
         },
         projection={type=:albersUsa},
     )
     return p
 end


text:
######################################################################
# infotheoretic.jl
# has things like mutual information and stuff
######################################################################

# sum over all possible jammer locations
function p_obs(m::SearchDomain, xv::Float64, yv::Float64, o::Obs)
	prob = 0.0
	for xj = 1:m.num_cells
		for yj = 1:m.num_cells
			prob += m.b[xj,yj] * O(m, xv, yv, (xj-1,yj-1), o)
		end
	end
	return prob
end

# computes mutual information for a specific vehicle location
function mutual_information(m::SearchDomain, xv::Float64, yv::Float64)
	H_o = 0.0
	H_o_t = 0.0
	for o = 0:35
		po = p_obs(m, xv, yv, o)
		if po > 0.0
			H_o -= po * log(po)
		end

		# sum over theta
		for xj = 1:m.num_cells
			for yj = 1:m.num_cells
				pot = O(m, xv, yv, (xj-1,yj-1), o)
				if pot > 0.0
					H_o_t -= pot * m.b[xj,yj] * log(pot)
				end
			end
		end
	end
	return H_o - H_o_t
end

# computes mutual information for all locations
function mutual_information(m::SearchDomain)
	mut_info = zeros(m.num_cells, m.num_cells)
	for xv = 1:m.num_cells
		for yv = 1:m.num_cells
			mut_info[xv,yv] = mutual_information(m, xv-1.0, yv-1.0)
		end
	end
	return mut_info
end


text:
# DeclarationName
DeclarationName() = DeclarationName(clang_DeclarationName_create())

function DeclarationName(x::IdentifierInfo)
    @check_ptrs x
    return DeclarationName(clang_DeclarationName_createFromIdentifierInfo(x))
end

dump(x::DeclarationName) = clang_DeclarationName_dump(x)

isEmpty(x::DeclarationName) = clang_DeclarationName_isEmpty(x)

function getAsString(x::DeclarationName)
    str = clang_DeclarationName_getAsString(x)
    s = unsafe_string(str)
    clang_DeclarationName_disposeString(str)
    return s
end


text:
const CFV = CompositeFieldVector

length(::CFV{T,L,N}) where {T,L,N} = N
size(::CFV{T,L,N}) where {T,L,N} = (N,)
vec(v::CFV) = getindex(v, :)
copy(v::CFV) = typeof(v)([getfield(v, f) for f in fieldnames(v)]...)
endof(v::CFV{T,L,N}) where {T,L,N} = N

getindex(v::CFV{T,L}, I::AbstractArray) where {T,L} = [lookupget(v, T, L, i) for i in I]
getindex(v::CFV{T,L,N}, I::Colon) where {T,L,N} = [lookupget(v, T, L, i) for i in 1:N]
getindex(v::CFV{T,L}, i::Int) where {T,L} = lookupget(v, T, L, i)

setindex!(v::CFV{T,L}, x, i::Int) where {T,L} = lookupset!(v, T, L, i, x)

# Build functions programmatically to handle up to 200 fields
for num_args = 1:200
    types = join([Symbol.(string.("_T",i)) for i = 1:num_args], ",")

    ex = parse("lookupget(v::V, t::Type{T}, f::Type{Tuple{$types}}, i) where {T, V,$types} =
                    compositeget(v, t, ($types)[i])::T")
    eval(ex)

    ex = parse("lookupset!(v::V, t::Type{T}, f::Type{Tuple{$types}}, i, x) where {T, V,$types} =
                    compositeset!(v, t, ($types)[i], x)::T")
    eval(ex)
end

get_string(t::Tuple{T,Vararg}) where T = "getf($(get_string(Base.tail(t))), $(t[1]))"
get_string(t::Tuple{}) = "v"

# Maximum tree depth is 10
for num_args = 1:20
    types = tuple([Symbol.(string.("_T",i)) for i = 1:num_args]...)
    types_str = join(types, ",")
    revtypes = reverse(types)

    getstr = get_string(revtypes)
    ex = parse("compositeget(v::V, t, f::Type{Tuple{$types_str}}) where {V,$types_str} = $getstr")
    eval(ex)

    getstr = "setf!($(get_string(Base.tail(revtypes))), $(revtypes[1]), x)"
    ex = parse("compositeset!(v::V, t, f::Type{Tuple{$types_str}}, x) where {V,$types_str} = $getstr")
    eval(ex)
end

getf(v::V, f) where V <: Tuple = v[f]
getf(v, f) = strip(getfield(v, f))
strip(f) = f
setf!(v::V, f, x) where V = setftyped!(v, f, fieldtype(V, f), x)
setftyped!(v, f, t, x) = setfield!(v, f, x)

@require Unitful begin
    strip(f::T) where T <: Quantity = f.val
    setftyped!(v, f, t::Type{T}, x) where T <: Quantity = begin
        setfield!(v, f, x * oneunit(t))
        x
    end
end


text:
# This file is a part of BayesianBlocks.jl, licensed under the MIT License (MIT).

using StatsBase, ProgressMeter

"""
    partition(data, logfitness=:cash, logprior=:p0, progress=false, ...)

Return an array of optimal change points for a set of one-dimensional data
(observations or an histogram of them). This is the implementation of the
bayesian blocks algorithm, as outlined in [^1].

# Arguments
- `data`: numeric array or a `StatsBase.Histogram`
- `logfitness`: log of the block fitness function to be used, choose between
  [:cash]
- `logprior`: log of the prior distribution on the number of blocks to be used,
   choose between [:gamma, :p0]
- `gamma`, `p0`...: set the parameter value for the specified prior distribution
- `progress`: display a progress bar for long computations

# Example
```julia
using Distributions, StatsBase, Plots, LinearAlgebra

data = vcat(rand(Normal(0),1000),rand(Cauchy(5),1000))
data = data[(data .> -5) .& (data .< 10)]

h = fit(Histogram, data, -5:0.1:10)

# choose to use all data or an histogram of it!
hb = fit(Histogram, data, BayesianBlocks.partition(data))
hb = fit(Histogram, data, BayesianBlocks.partition(h))

plot(data, normalized=true, st=:stephist, nbins=1000)
plot!(normalize(hb), st=:step, w = 3)
```

## Performance tips
You can convert your data container to a less precise representation to improve
the performance a bit, e.g.
```julia
x::Array{Float32} = [1.1, π, (√5-1)/2]
```

[^1]: Scargle, J et al. (2012) [https://doi.org/10.1088/0004-637X/764/2/167]
"""
function partition(x::Vector{<:Real}; kwargs...)
    # take care of repeated data
    x_sorted = sort(x)
    x_unique = [x_sorted[1]]
    x_weight::Vector{Int32} = [1]
    for i in 2:length(x_sorted)
        if x_sorted[i] == x_sorted[i-1]
            x_weight[end] += 1
        else
            push!(x_unique, x_sorted[i])
            push!(x_weight, 1)
        end
    end

    partition(x_unique, x_weight; kwargs...)
end

function partition(h::Histogram{<:Integer,1}; kwargs...)
    v = collect(h.edges[1])
    x_unique = [0.5 * (v[i] + v[i+1]) for i in 1:length(v) - 1]
    x_weight::Vector{Int32} = h.weights

    # delete empty bins
    deleteat!(x_unique, findall(iszero, x_weight))
    deleteat!(x_weight, findall(iszero, x_weight))

    partition(x_unique, x_weight; kwargs...)
end

function partition(x_unique::Vector{<:Real}, x_weight::Vector{Int32};
                   logfitness::Symbol=:cash, logprior::Symbol=:p0,
                   gamma=0.01, p0=0.01, progress::Bool=false)

    # final number of data points
    N = length(x_unique)

    # pre-defined (log)fitness functions
    logf_dict = Dict(
        :cash => (N_k, T_k) -> N_k * log(N_k/T_k)
    )

    # pre-defined (log)prior distributions on Nblocks
    logp_dict = Dict(
        # simple γ^Nblocks prior
        :gamma => (γ=gamma) -> log(γ),
        # Note that there was a mistake in this equation in the original Scargle
        # paper (the "log" was missing). The following corrected form is taken
        # from https://arxiv.org/abs/1304.2818
        :p0    => (Np=N, p=p0) -> log(73.53 * p * Np^(-0.478)) - 4
    )

    # check input
    !haskey(logf_dict, logfitness) && error("$logfitness function not defined!")
    !haskey(logp_dict, logprior)   && error("$logprior function not defined!")

    # save prior value for later computation
    ncp_prior = logp_dict[logprior]()

    # array of (all possible) block edges
    edges = vcat(x_unique[1],
                 0.5f0(x_unique[1:end-1]+x_unique[2:end]),
                 x_unique[end])

    # see Sec. 2.6 in [^1]
    best = Vector{Float32}()
    last = Vector{Int32}()

    # display progress bar for long computations
    # total number of steps: ∑n(n-1) = N(N^2-1)/3
    if progress && Sys.WORD_SIZE == 32
        @warn "Progress bar not supported on 32-bit Julia, disabling"
        progress = false
    elseif progress
        p = Progress(UInt(N*(N^2-1)/3), 2); m = 0
    end

    for k in 1:N
        # define nice alias to mimic the notation used in [^1]
        F(r) = logf_dict[logfitness](cumsum(x_weight[r:k])[end], edges[k+1]-edges[r]) + ncp_prior

        # compute all possible configurations (Eq. (8) in [^1])
        A = [F(r) + (r == 1 ? 0 : best[r-1]) for r in 1:k]

        # save best configuration
        push!(last, argmax(A))
        push!(best, maximum(A))

        progress && update!(p, m += k*(k-1))
    end

    # extract changepoints by iteratively peeling off the last block
    cp = Vector{Int32}()
    i = N+1
    while i != 0
        push!(cp, i)
        i = (i == 1 ? 0 : last[i-1])
    end

    return [edges[j] for j in cp[end:-1:1]]
end


text:
export rand_eg, empty_eg, chain_eg, ring_eg, K_eg, disconnected_cliques_eg
export petersen_graph, c60_graph, get_petersen_graph, get_c60_graph
export inverse_eg, smallworld_eg, random_regular_eg

"""
    K_eg(nv::Int, [nw]) -> EliminateGraph

fully connected graph, if `nw` provided, compute `K(nv,nw)` graph.
"""
K_eg(nv::Int) = rand_eg(nv, 1.0)
K_eg(m::Int, n::Int) = EliminateGraph(m+n,vcat([vi=>vj for vi=1:m, vj=m+1:m+n]...))

"""
Inverse graph, arXiv: 1807.03739
"""
function inverse_eg(n::Int)
    EliminateGraph(n, vcat([x=>mod1(x+1, n) for x=1:n], [x=>mod(-mod_inverse(x-1, n), n)+1 for x=2:n]))
end

"""

"""
function smallworld_eg(n::Int; ρ::Real)
    longlinks = Pair{Int,Int}[]
    for x=1:n, y=1:n
        if abs(mod(x-y, n)) > 1 && rand() < ρ
            push!(longlinks, x=>y)
        end
    end
    EliminateGraph(n, vcat([x=>mod(x, n)+1 for x=1:n], longlinks))
end

function random_regular_eg(nv::Int, nreg::Int)
    eg = empty_eg(nv)
    g = random_regular_graph(nv, nreg)
    for i in 1:nv
        for j=g.fadjlist[i]
            unsafe_connect!(eg,i,j)
        end
    end
    return eg
end

"""
    mod_inverse(x::Int, N::Int) -> Int
Return `y` that `(x*y)%N == 1`, notice the `(x*y)%N` operations in Z* forms a group and this is the definition of inverse.
"""
function mod_inverse(x::Int, N::Int)
    for i=1:N
        (x*i)%N == 1 && return i
    end
    throw(ArgumentError("Can not find the inverse, $x is probably not in Z*($N)!"))
end

"""
A graph consists of disconnected cliques.
"""
function disconnected_cliques_eg(nvs::Int...)
    N = sum(nvs)
    tbl = zeros(Bool, N, N)
    offset = 0
    for nv in nvs
        for j=offset+1:offset+nv
            for i=offset+1:offset+nv
                i!=j && (@inbounds tbl[i,j] = true)
            end
        end
        offset += nv
    end
    return EliminateGraph(tbl)
end

"""
    empty_eg(nv::Int) -> EliminateGraph

fully disconnected graph
"""
empty_eg(nv::Int) = EliminateGraph(zeros(Bool,nv,nv))

"""
    rand_eg(nv::Int, density::Real) -> EliminateGraph

Generate a random `EliminateGraph`.
"""
function rand_eg(nv::Int, density::Real)
    tbl = rand(nv, nv) .< density
    copyltu!(tbl)
    for i=1:nv
        tbl[i,i] = false
    end
    EliminateGraph(tbl)
end

"""chain graph"""
chain_eg(nv::Int) = EliminateGraph(nv,[vi=>vi+1 for vi=1:nv-1])

"""ring graph"""
ring_eg(nv::Int) = EliminateGraph(nv,[vi=>mod1(vi+1,nv) for vi=1:nv])

# Special Graphs
get_petersen_graph() = EliminateGraph(10,
                    [1=>2, 2=>3, 3=>4, 4=>5, 5=>1,
                    1=>6, 2=>7, 3=>8, 4=>9, 5=>10,
                    6=>8, 7=>9, 8=>10, 9=>1, 10=>2])
const petersen_graph = get_petersen_graph()

get_c60_graph() = EliminateGraph(60, [1=>10, 1=>41, 1=>59, 2=>12, 2=>42, 2=>60, 3=>6, 3=>
        43, 3=>57, 4=>8, 4=>44, 4=>58, 5=>13, 5=>56, 5=>
        57, 6=>10, 6=>31, 7=>14, 7=>56, 7=>58, 8=>12, 8=>
        32, 9=>23, 9=>53, 9=>59, 10=>15, 11=>24, 11=>53, 11=>
        60, 12=>16, 13=>14, 13=>25, 14=>26, 15=>27, 15=>
        49, 16=>28, 16=>50, 17=>18, 17=>19, 17=>54, 18=>
        20, 18=>55, 19=>23, 19=>41, 20=>24, 20=>42, 21=>
        31, 21=>33, 21=>57, 22=>32, 22=>34, 22=>58, 23=>
        24, 25=>35, 25=>43, 26=>36, 26=>44, 27=>51, 27=>
        59, 28=>52, 28=>60, 29=>33, 29=>34, 29=>56, 30=>
        51, 30=>52, 30=>53, 31=>47, 32=>48, 33=>45, 34=>
        46, 35=>36, 35=>37, 36=>38, 37=>39, 37=>49, 38=>
        40, 38=>50, 39=>40, 39=>51, 40=>52, 41=>47, 42=>
        48, 43=>49, 44=>50, 45=>46, 45=>54, 46=>55, 47=>
        54, 48=>55])
const c60_graph = get_c60_graph()


text:
function UnivariateDerivative()
    return UnivariateDerivative(())
end

function get_derivative(obj::UnivariateDerivative, arg0::jint)
    return jcall(obj, "getDerivative", jdouble, (jint,), arg0)
end

function get_free_parameters(obj::UnivariateDerivative)
    return jcall(obj, "getFreeParameters", jint, ())
end

function get_partial_derivative(obj::UnivariateDerivative, arg0::Vector{jint})
    return jcall(obj, "getPartialDerivative", jdouble, (Vector{jint},), arg0)
end

function to_derivative_structure(obj::UnivariateDerivative)
    return jcall(obj, "toDerivativeStructure", DerivativeStructure, ())
end


text:
# export downloader

# ====================================================
# =================    Downloader    =================
# ====================================================
# export downloader

# download only audio
function downloader(audio, video::Nothing, path, metadata, opts)
    dest = string(path, '.', audio.ext)
    dl_file(audio, dest, metadata, opts)
end

# download audio and video with muxing
function downloader(audio, video::Video, path, metadata, opts)
    a_temp = string(path, "___audiotemp.", audio.temp_ext)
    dest = string(path, '.', video.ext)
    dl_file(audio, a_temp, opts)
    usetemp = get(opts, :use_tempfile, true)
    if usetemp
        v_temp = string(path, "___videotemp.", video.temp_ext)
        complete = dl_file(video, v_temp, opts)
        if complete
            mux(a_temp, v_temp, dest, metadata, opts)
        end
        rm(v_temp)
    else
        complete = dl_file(video, dest, a_temp, metadata, opts)
        complete || rm(dest)
    end
    rm(a_temp)
    return nothing
end

# ====================================================
# ==================    Payload    ===================
# ====================================================
# Audio temporary file
function dl_file(audio::Audio, dest, opts)
    ffpath, chunk_size = dl_options(opts)
    io = `$ffpath -i pipe:0 -c copy $dest`
    get_payload(audio, io, chunk_size)
end

# Audio file
function dl_file(audio::Audio, dest, metadata, opts)
    ffpath, chunk_size = dl_options(opts)
    io = `$ffpath -i pipe:0 -c copy $metadata $dest`
    get_payload(audio, io, chunk_size)
end

# Video temporary file
function dl_file(video::Video, dest, opts)
    ffpath, chunk_size = dl_options(opts)
    get_payload(video, dest, chunk_size)
end

# Video file (muxing video stream and audio temporary file)
function dl_file(video::Video, dest, temp, metadata, opts)
    ffpath, chunk_size = dl_options(opts)
    qs = get(opts, :queue_size, 10)
    io = `$ffpath -thread_queue_size $qs -i pipe:0 -i $temp -c copy $metadata $dest`
    get_payload(video, io, chunk_size)
end

#
function get_payload(av, io_target, chunk_size)
    open(io_target, "w") do io
        chunk = chunk_size * 1024 * 1024
        start = 0
        stop = chunk - 1
        if iszero(av.length)
            h = try
                resp = HTTP.head(av.url)
                Dict(resp.headers)
            catch e
                printstyled("Target file not found. Status:$(e.status)\n", color=:red)
                return false
            end
            len = h["Content-Length"]
        else
            len = av.length
        end
        while stop ≤ len
            chunked_payload(av.url, start, stop, io)
            start = stop + 1
            stop += chunk
        end
        chunked_payload(av.url, start, len, io)
        return true
    end
end

# youtube doesn't send data using Transfer-Encoding:chunked
# for accelerating download speed, YTDL.jl use RangeRequest instead
function chunked_payload(url, start, stop, io)
    header = ["Range"=>"bytes=$(start)-$(stop)"]
    HTTP.open("GET", url, header) do stream
        write(io, stream)
    end
end

# set default value to download option
function dl_options(opts)
    ff = get(opts, :ffmpeg_path, "ffmpeg")
    size = get(opts, :chunk_size, 10)
    ff, size
end

# temporary file post processsing
function mux(a, v, dest, metadata, opts)
    ffpath = get(opts, :ffmpeg_path, "ffmpeg")
    run(`$ffpath -i $v -i $a -c copy $metadata $dest`)
end


text:
module ModuleDraw
using MultipleViewGeometry.ModuleTypes, MultipleViewGeometry.ModuleOperators, MultipleViewGeometry.ModuleMathAliases
using StaticArrays
using Makie
using GeometryTypes
#using Plots #, Plotly

#export draw!, EpipolarLineGraphic, LineSegment3D, PlaneSegment3D, Camera3D
#export WorldCoordinateSystem3D
#include("draw.jl")
export draw!
include("visualize.jl")
end


text:
mutable struct InfiniteMatrixProductState{T<:RLorCX}
    d  :: Int
    Γ  :: Array{T, 3}
    Λ  :: Vector{Float64}
end

function InfiniteMatrixProductState(M::Array{T, 3}) where{T<:RLorCX}
    diml, d, dimr = size(M)
    @assert (diml == dimr)
    fact = svd(reshape(M, diml, d*dimr), full=false)
    Vt = reshape(fact.Vt, diml, d, dimr)

    @tensor Γ[l,d,r] := Vt[l,d,m] * fact.U[m,r]

    InfiniteMatrixProductState{T}(d, (Γ,), fact.S/norm(fact.S))
end

"""
    poshermitsqrt(A)

return the square root of the positvie hermitian matrix `A` using
eigenvalue decomposition.

"""
function poshermitsqrt(A::Matrix{T}) where{T<:Number}
    !isapprox(A, A') && error("Matrix is not hermitian! ", norm(A-A'))
    S, U = eigen(Hermitian(A))
    if !all(S .> 0)
        !all(S .< 0) && error("Matrix is not positive!", S)
        S = -S
    end

    U * Diagonal(sqrt.(S)) * U'
end

function _transfermatvec(v::Matrix{T}, uc::Array{T, 4}, direction::Symbol) where{T<:Number}
    if direction == :LEFT
        @tensor w[u,d] := (v[u',d'] * uc[u',o1,o2,u]) * conj(uc)[d',o1,o2,d]
        w
    elseif direction == :RIGHT
        @tensor w[u,d] := (v[u',d'] * uc[u,o1,o2,u']) * conj(uc)[d,o1,o2,d']
        w
    else
        error("unknown direction!")
    end
end

function InfiniteMatrixProductState(A::Array{T, 3},
                                    Λ::Vector{Float64},
                                    B::Array{T, 3},
                                    Λex::Vector{Float64}) where {T<:RLorCX}
    n = size(A, 1)
    d = size(A, 2)
    @assert size(A) == (n, d, n)
    @assert length(Λ) == n
    @assert size(B) == (n, d, n)
    @assert length(Λex) == n

    # @tensor LB[l,o,r] := Diagonal(Λ)[l,m] * B[m,o,r]
    # Q, Λr = qr(reshape(LB, n*d, n))
    # Q = Matrix(Q)
    @tensor Luc[l,o1,o2,r] := (A[l,o1,ml] * Diagonal(Λ)[ml,mm]) *
        B[mm,o2,mr] * Diagonal(1 ./Λex)[mr,r]

    es, vs, info = eigsolve(v->_transfermatvec(v, Luc, :LEFT),
                            rand(n,n), 1, :LR)
    if !(es[1] ≈ 1.0)
       println("The largest left eigenvalue of transfer matrix is not unit : ", es[1])
    end
    v = vs[1]
    if !(v ≈ conj(v))
        println("The left eigenvector is not real!", norm(v-conj(v)))
        v = real.(v)
    end
    X = poshermitsqrt(real.(v))

    # @tensor AL[l,o,r] := A[l,o,m] * Diagonal(Λ)[m,r]
    # Q, Λl = qr(transpose(reshape(AL, n, d*n)))
    # Q = transpose(Matrix(Q))
    # Λl = transpose(ΛL)
    @tensor Ruc[l,o1,o2,r] := (Diagonal(1 ./Λex)[l,ml] * A[ml,o1,mm]) *
        Diagonal(Λ)[mm, mr] * B[mr,o2,r]

    es, vs, info = eigsolve(v->_transfermatvec(v, Ruc, :RIGHT),
                            rand(n,n), 1, :LR)
    if !(es[1] ≈ 1.0)
        println("The largest right eigenvalue of transfer matrix is not unit : ", es[1])
    end
    v = vs[1]
    if !(v ≈ conj(v))
        println("The right eigenvector is not real!", norm(v-conj(v)))
    end
    Y = poshermitsqrt(real.(v))

    ###NOTE:
    ## Left normalzied canonical unitcell is (XA) Λ (BY) inv(Y Λex X)
    ## Right noralized canonical unitcell is inv(Y Λex X) (XA) Λ (BY)
    ## so naming A <- XA and B <-BY and Λex <- Y Λex X works!

    #println(X - conj(X))
    @tensor XA[l,o,r] := X[l,m] * A[m,o,r]
    Q, R = qr(reshape(XA, n*d, n))
    Γ1 = reshape(Matrix(Q), n,d,n)

    @tensor BY[l,o,r] := (R[l,ml] * Diagonal(Λ)[ml,mm]) * B[mm,o,mr] * Y[mr,r]
    Q, R = qr(reshape(BY, n*d, n))
    Γ2 = reshape(Matrix(Q), n,d,n)

    @show norm(R - Matrix(I, n, n))

    L = Y*Diagonal(Λex)*X
    norm(L - Diagonal(diag(L)))
    InfiniteMatrixProductState{T}(d, Γ1, inv(Y*Diagonal(Λex)*X))
end

function ketstate2imps(vstate::Vector{T}, lx::Int, d::Int;
                       maxdim::Int=d^div(lx, 2)) where{T<:RLorCX}
    @assert lx % 2 == 1
    @assert length(vstate) == d^lx

    l = div(lx, 2)
    fact = svd(reshape(vstate, d^(l+1), d^l), full=false)

    S, n, ratio = truncate(fact.S, maxdim=maxdim)
    Vt = fact.Vt[1:n, :]

    ##NOTE: Using inversion symmetry the left and right l sites have
    ##the same vector, so we only need to put an identity in the
    ##middle, conjugate and then contract to get rid of all the left
    ##and right tensors except the one in the middle (which we are
    ##looking for)

    psi_b = reshape(Vt, n, 1, d*ones(Int, l)...)
    psi_a = permutedims(psi_b, [1,2,(l+2:-1:3)...])

    ten_b = reshape(conj(psi_b), n, 1, d^l)
    ten_a = reshape(conj(psi_a), n, 1, d^l)
    I_ext = reshape(Matrix(1.0I, d, d), 1, d, 1, d)
    vstate_aib = reshape(vstate, d^l, 2, d^l)

    @tensor M[a,d,b] := ((vstate_aib[a',d', b'] * ten_a[a, il, a']) *
                         ten_b[b,ir, b']) * I_ext[il, d, ir, d']

    InfiniteMatrixProductState(M)
end

function measure_bond(imps::InfiniteMatrixProductState{T},
                      bond::Array{T, 4}) where{T<:Number}
    d = imps.d
    @assert size(bond) == (d,d,d,d)

    Γ = imps.Γs#[1]
    @tensor sbb[l,d1,d2,r] :=
        Diagonal(imps.Λ)[l, ml] * Γ[ml, d1, mr] * Γ[mr, d2, r]
    @tensor v = scalar(sbb[l,d1',d2',r] * bond[d1,d2,d1',d2'] * conj(sbb)[l,d1,d2,r])
    real(v)
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
mutable struct Section <: Element
    contents::Vector{Element}
    planning::Union{Planning, Nothing}
    properties::Union{PropertyDrawer, Nothing}
end

Section(contents::Vector=Element[]) = Section(Vector{Element}(contents), nothing, nothing)

mutable struct Heading <: Element
    level::Integer
    keyword::Union{SubString{String}, Nothing}
    priority::Union{SubString{String}, Nothing}
    title::Vector{Object}
    tags::Vector{SubString{String}}
    section::Union{Section, Nothing}
end

Heading(level::Integer, title::Vector{Object}, section::Union{Section, Nothing}=nothing;
        keyword::Union{String, Nothing}=nothing,
        priority::Union{String, Nothing}=nothing,
        tags::Vector{String}=String[]) =
            Heading(level, if !isnothing(keyword) SubString(keyword) end,
                    if !isnothing(priority) SubString(priority) end,
                    title, SubString.(tags), section)

Heading(level::Integer, title::String, section::Union{Section, Nothing}=nothing;
        keyword::Union{String, Nothing}=nothing,
        priority::Union{String, Nothing}=nothing,
        tags::Vector{String}=String[]) =
            Heading(level, Object[TextPlain(title)], section; keyword, priority, tags)


text:
include("gridworld.jl")
include("chain.jl")
include("cartpole.jl")
include("pinball/pinball.jl")
include("bicycle.jl")
include("mountaincar.jl")
include("acrobot.jl")
include("lightbox.jl")

export GridWorld, Chain
export CartPole, CartPoleParams
export PinBall, PinBallConfig
export Bicycle, BicycleParams
export MountainCar, MountainCarParams
export Acrobot, AcrobotParams
export LightBox

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
package ln

struct Scene
    Shapes::Array{Shape}
    Tree::Tree
end

type Scene struct {
	Shapes []Shape
	Tree   *Tree
}

func (s *Scene) Compile() {
	for _, shape := range s.Shapes {
		shape.Compile()
	}
	if s.Tree == nil {
		s.Tree = NewTree(s.Shapes)
	}
}

func (s *Scene) Add(shape Shape) {
	s.Shapes = append(s.Shapes, shape)
}

func (s *Scene) Intersect(r Ray) Hit {
	return s.Tree.Intersect(r)
}

func (s *Scene) Visible(eye, point Vector) bool {
	v := eye.Sub(point)
	r := Ray{point, v.Normalize()}
	hit := s.Intersect(r)
	return hit.T >= v.Length()
}

func (s *Scene) Paths() Paths {
	var result Paths
	for _, shape := range s.Shapes {
		result = append(result, shape.Paths()...)
	}
	return result
}

func (s *Scene) Render(eye, center, up Vector, width, height, fovy, near, far, step float64) Paths {
	aspect := width / height
	matrix := LookAt(eye, center, up)
	matrix = matrix.Perspective(fovy, aspect, near, far)
	return s.RenderWithMatrix(matrix, eye, width, height, step)
}

func (s *Scene) RenderWithMatrix(matrix Matrix, eye Vector, width, height, step float64) Paths {
	s.Compile()
	paths := s.Paths()
	if step > 0 {
		paths = paths.Chop(step)
	}
	paths = paths.Filter(&ClipFilter{matrix, eye, s})
	if step > 0 {
		paths = paths.Simplify(1e-6)
	}
	matrix = Translate(Vector{1, 1, 0}).Scale(Vector{width / 2, height / 2, 0})
	paths = paths.Transform(matrix)
	return paths
}


text:
#
# Propagator for the NBody problem.
#

"""
    nbody_tic

Currently not exported. Used for n-body numerical integration.
"""
function NBodyTic!(∂u, u, p, t=0)
    for i = 1:length(u.body)
        ∂u.body[i].r = u.body[i].v
        ∂u.body[i].v = zero.(∂u.body[i].v)
        for j = 1:length(u.body)
            if i ≠ j
                ∂u.body[i].v += ((p.m[i] * p.m[j]) / norm(u.body[j].r .- u.body[i].r)^3 * (u.body[j].r .- u.body[i].r))
            end
        end
        ∂u.body[i].v *= (p.G / p.m[i])
    end
    return nothing
end

"""
Uses OrdinaryDiffEq solvers to propagate `sys` Δt into the future.
All keyword arguments are passed directly to OrdinaryDiffEq solvers.

References:
* [1] https://diffeq.sciml.ai/v4.0/tutorials/ode_example.html
* [2] https://github.com/SciML/DifferentialEquations.jl/issues/393#issuecomment-658210231
* [3] https://discourse.julialang.org/t/smart-kwargs-dispatch/14571/15
"""
function propagate(sys::NBodySystem{N,T}, Δt::Unitful.Quantity; kwargs...) where N where T

    # Integration options
    defaults = (;  reltol=1e-14, abstol=1e-14)
    options = merge(defaults, kwargs)

    # Initial conditions
    u₀ = ComponentArray(body=(map(b -> ComponentArray((r=ustrip.(u"m", b.r), v=ustrip.(u"m/s", b.v))), sys.bodies)))
    ts = T.(ustrip.(u"s", (zero(Δt), Δt)))
    p  = ComponentArray((G=6.6743e-11, m=map(b->ustrip(u"kg",b.m), sys.bodies)))

    # Integrate!
    sols = solve(ODEProblem(NBodyTic!, u₀, ts, p); options...)

    # Unpack and return
    bodies = map(x->NBodySystem(
                        map(i->Body(u"m" * x.body[i].r, u"m/s" * x.body[i].v, sys.body[i].m),
                        1:length(sys.body))),
                sols.u)

    return Trajectory(
        bodies,
        u"s" * sols.t,
        sols.retcode
    )

end


text:
using Plots
using SpikingNeuralNetworks
SNN.@load_units

ear = SNN.Rate(;N = 2)
ear.r = zeros(ear.N) # exp(- (0:ear.N-1).^2 / 10^2)
ppc = SNN.Rate(;N = 200)
ear_ppc = SNN.PINningSynapse(ear, ear; σ = 1.5, p = 1.0)
ppc_ppc = SNN.PINningSynapse(ppc, ppc; σ = 1.5, p = 1.0)
P, C = [ppc], [ear_ppc, ppc_ppc]

SNN.monitor(ppc_ppc, [(:g, [1])])

A = 1.3 / 1.5; fr = 1 / 60ms
f(t) = (A /1.0) * sin(1π * fr * t) +  (A / 2.0) * sin(2π * fr * t) +
(A / 6.0) * sin(3π * fr * t) +   (A / 3.0) * sin(4π * fr * t)

ts = 0:0.1ms:1440ms
for (i, t) in enumerate(ts)
    ppc_ppc.f .= [f(t); ppc_ppc.g[2:end]]
    SNN.train!(P, C, 0.1ms, t)
end

SNN.vecplot(ppc_ppc, :g); plot!(f.(ts))


text:
include(joinpath("..","..","src","MaxwellSolver","FDTD1D.jl"))
import .FDTD1D

function test_fdtd_1d_propagation(Δx, Δt, cfl, f₀, f₁, exp_norm_dev, args...; plot_needed=false, kwargs...)
    if plot_needed
        plt = plot(f₀.ey, label = "initial")
        plot!(f₁.ey, label = "expected")
    end

    shape(t,x) = 0.0

    pulse_shape = (
        y = shape,
        z = shape
    )

    f = deepcopy(f₀)

    advance_fields! = FDTD1D.make_advance_fields(f, cfl, pulse_shape, Δt, Δx, 0, FDTD1D.PML(0,1.0,Δx,Δt))

    j = (
        y = zeros(length(f.ey)),
        z = zeros(length(f.ey))
    )

    for t in 1:10
        advance_fields!(t*Δt, j)
    end

    if plot_needed
        plot!(f.ey, label = "calculated")
    end

    s = sum(@. (f.ey[10:end-1] - f₁.ey[10:end-1])^2)
    println("FDTD1D propagation $args $(values(values(kwargs))): $(s)")
    @test s ≈ 0 atol=exp_norm_dev

    if plot_needed
        display(plt)
    end
end

function test_fdtd_1d_pml(Δx, Δt, cfl, f₀, f₁, exp_norm_dev, args...; plot_needed=false, kwargs...)
    if plot_needed
        plt = plot(f₁.ey, label = "expected")
    end

    shape(t,x) = 0.0

    pulse_shape = (
        y = shape,
        z = shape
    )

    f = deepcopy(f₀)

    advance_fields! = FDTD1D.make_advance_fields(f, cfl, pulse_shape, Δt, Δx, 0, FDTD1D.PML(10,1e3,Δx,Δt))

    j = (
        y = zeros(length(f.ey)),
        z = zeros(length(f.ey))
    )

    for t in 1:200
        advance_fields!(0, j)
    end

    if plot_needed
        plot!(f.ey, label = "calculated")
    end

    s = sum(@. (f.ey - f₁.ey)^2)
    println("FDTD1D PML $args $(values(values(kwargs))): $(s)")
    @test s ≈ 0 atol=exp_norm_dev

    if plot_needed
        display(plt)
    end
end

function test_fdtd_1d_generation(Δx, Δt, cfl, f₀, f₁, pulse_shape, exp_norm_dev, args...; plot_needed=false, kwargs...)
    if plot_needed
        plt = plot(f₀.ey, label = "initial")
        plot!(f₁.ey, label = "expected")
    end

    f = deepcopy(f₀)

    advance_fields! = FDTD1D.make_advance_fields(f, cfl, pulse_shape, Δt, Δx, 0, FDTD1D.PML(0,1.0,Δx,Δt))

    j = (
        y = zeros(length(f.ey)),
        z = zeros(length(f.ey))
    )

    for t in 1:100
        advance_fields!(t*Δt, j)
    end

    if plot_needed
        plot!(f.ey, label = "calculated")
    end

    s = sum(@. (f.ey - f₁.ey)^2)
    println("FDTD1D generation $args $(values(values(kwargs))): $(s)")
    @test s ≈ 0 atol=exp_norm_dev

    if plot_needed
        display(plt)
    end
end

function test_fdtd_1d_current(Δx, Δt, cfl, f₀, f₁, j, exp_norm_dev, args...; plot_needed=false, kwargs...)
    if plot_needed
        plt = plot(f₀.ey, label = "initial")
        plot!(f₁.ey, label = "expected")
    end

    shape(t,x) = 0.0

    pulse_shape = (
        y = shape,
        z = shape
    )

    f = deepcopy(f₀)

    advance_fields! = FDTD1D.make_advance_fields(f, cfl, pulse_shape, Δt, Δx, 0, FDTD1D.PML(0,1.0,Δx,Δt))

    for t in 1:100
        advance_fields!(t*Δt, j(t*Δt))
    end

    if plot_needed
        plot!(f.ey, label = "calculated")
    end

    s = sum(@. (f.ey - f₁.ey)^2)
    println("FDTD1D current $args $(values(values(kwargs))): $(s)")
    # @test s ≈ 0 atol=exp_norm_dev

    if plot_needed
        display(plt)
    end
end

@testset "Test 1D FDTD solvers" begin
    Δx = 0.01
    Δt = 0.8*Δx

    f₀ = FDTD1D.YeeMesh1D{Float64}(110)
    f₀.ey[2:102] = [sin(2π*i*Δx) for i = 0:100]
    f₀.hz[2:101] = [sin(2π*((i+0.5)*Δx-0.5*Δt)) for i = 0:99]

    f₁ = FDTD1D.YeeMesh1D{Float64}(110)
    f₁.ey[10:110] = [sin(2π*(i*Δx)) for i = 0:100]
    f₁.hz[10:109] = [sin(2π*((i+0.5)*Δx-0.5*Δt)) for i = 0:99]

    test_fdtd_1d_propagation(Δx, Δt, Δt/Δx, f₀, f₁, 1e-3)

    f₁ = FDTD1D.YeeMesh1D{Float64}(110)

    test_fdtd_1d_pml(Δx, Δt, Δt/Δx, f₀, f₁, 1e-5)

    f₀ = FDTD1D.YeeMesh1D{Float64}(100)

    f₁ = FDTD1D.YeeMesh1D{Float64}(100)
    f₁.ey[3:81] = [-sin(2π*(i*Δx-100*Δt)) for i = 1:79]
    f₁.hz[3:81] = [-sin(2π*((i+0.5)*Δx-100.5*Δt)) for i = 1:79]

    pulse_shape = (
        y = (t,x) -> sin(2π*(x-t)),
        z = (t,x) -> 0.0
    )

    test_fdtd_1d_generation(Δx, Δt, Δt/Δx, f₀, f₁, pulse_shape, 0.02)

    f₀ = FDTD1D.YeeMesh1D{Float64}(200)

    f₁ = FDTD1D.YeeMesh1D{Float64}(200)
    f₁.ey[100:178] = [-sin(2π*(i*Δx-100*Δt)) for i = 0:78]
    f₁.ey[22:100] = [-sin(2π*(i*Δx-100*Δt)) for i = 78:-1:0]
    f₁.hz[3:81] = [-sin(2π*((i+0.5)*Δx-100.5*Δt)) for i = 1:79]

    function j(t)
        jy = zeros(length(f₀.ey))
        jy[end÷2] = π/2*sin(2π*t)
        return (
            y = jy,
            z = zeros(length(f₀.ey))
        )
    end

    test_fdtd_1d_current(Δx, Δt, Δt/Δx, f₀, f₁, j, 0.1)
end


text:
## Specialized solver for sparse symmetric integer linear systems with exact rational solution

include("Modulos.jl")
using .Modulos
import Base.GMP: MPZ

using Base: OneTo
import LinearAlgebra: BlasInt, checknonsingular, LU, tril!, triu!, ipiv2perm,
                      lu, lu!, ldiv!, Factorization, issuccess
using SparseArrays
import SparseArrays: getcolptr

using BigRationals

function compat_lu(::Val{Tf}, B, maxvec, info) where Tf
    @static if VERSION < v"1.8-"
        LU{Tf,SparseMatrixCSC{Tf,Int}}(B, Vector{BlasInt}(1:maxvec), convert(BlasInt, info))
    else
        LU{Tf,SparseMatrixCSC{Tf,Int},OneTo{Int}}(Tf.(B), 1:maxvec, info)
    end
end

function compat_lu_convert(::Val{Tf}, B, maxvec, info) where Tf
    @static if VERSION < v"1.8-"
        LU{Tf,SparseMatrixCSC{Tf,Int}}(Tf.(B), Vector{BlasInt}(1:maxvec), convert(BlasInt, info))
    else
        LU{Tf,SparseMatrixCSC{Tf,Int},OneTo{Int}}(Tf.(B), 1:maxvec, info)
    end
end

function rational_lu!(B::SparseMatrixCSC, col_offset::Vector{Int}, check=true)
    Tf = eltype(B)
    m, n = size(B)
    minmn = min(m, n)
    info = 0
    @inbounds begin
        for k in 1:minmn
            ipiv = getcolptr(B)[k] + col_offset[k]
            piv = nonzeros(B)[ipiv]
            if iszero(piv)
                check && checknonsingular(k-1, Val(false)) # TODO update with Pivot
                return compat_lu(Val(Tf), B, minmn, k-1)
            end
            Bkkinv = inv(piv)
            @simd for i in ipiv+1:getcolptr(B)[k+1]-1
                nonzeros(B)[i] *= Bkkinv
            end
            for j in k+1:n
                r1 = getcolptr(B)[j]
                r2 = getcolptr(B)[j+1]-1
                r = searchsortedfirst(rowvals(B), k, r1, r2, Base.Forward)
                ((r > r2) || (rowvals(B)[r] != k)) && continue
                Bkj = nonzeros(B)[r]
                for i in ipiv+1:getcolptr(B)[k+1]-1
                    Bik = nonzeros(B)[i]
                    l = i - ipiv

                    while rowvals(B)[l+r] < rowvals(B)[i]
                        r += 1
                    end
                    nonzeros(B)[l+r] -= Bik * Bkj
                end
            end
        end
    end
    check && checknonsingular(info, Val(false))
    return compat_lu(Val(Tf), B, minmn, info)
end

# function lu!(B::SparseMatrixCSC{<:Rational}, ::Val{Pivot} = Val(false);
#                            col_offset, check::Bool = true) where Pivot
function rational_lu!(B::SparseMatrixCSC{BigRational}, col_offset::Vector{Int}, check::Bool=true)
    Tf = Rational{BigInt}
    m, n = size(B)
    minmn = min(m, n)
    info = 0
    Bkkinv = BigRational()
    tmp = BigRational()
    @inbounds begin
        for k in 1:minmn
            ipiv = getcolptr(B)[k] + col_offset[k]
            piv = nonzeros(B)[ipiv]
            if iszero(piv)
                check && checknonsingular(k-1, Val(false)) # TODO update with Pivot
                return compat_lu_convert(Val(Tf), B, minmn, k-1)
            end
            BigRationals.MPQ.inv!(Bkkinv, piv)
            @simd for i in ipiv+1:getcolptr(B)[k+1]-1
                BigRationals.MPQ.mul!(nonzeros(B)[i], Bkkinv)
            end
            for j in k+1:n
                r1 = getcolptr(B)[j]
                r2 = getcolptr(B)[j+1]-1
                r = searchsortedfirst(rowvals(B), k, r1, r2, Base.Forward)
                ((r > r2) || (rowvals(B)[r] != k)) && continue
                Bkj = nonzeros(B)[r]
                for i in ipiv+1:getcolptr(B)[k+1]-1
                    Bik = nonzeros(B)[i]
                    l = i - ipiv

                    while rowvals(B)[l+r] < rowvals(B)[i]
                        r += 1
                    end
                    # Base.GMP.MPZ.mul!(tmp, Bik, Bkj)
                    # Base.GMP.MPZ.sub!(nonzeros(B)[l+r], tmp)
                    BigRationals.MPQ.mul!(tmp, Bik, Bkj)
                    BigRationals.MPQ.sub!(nonzeros(B)[l+r], tmp)
                end
            end
        end
    end
    check && checknonsingular(info, Val(false))
    return compat_lu_convert(Tf, B, minmn, info)
end

# function lu(A::SparseMatrixCSC{<:Rational}, pivot::Union{Val{false}, Val{true}} = Val(false); check::Bool = true)
function rational_lu(A::SparseMatrixCSC, check::Bool=true, ::Type{Ti}=BigRational) where {Ti}
    Tf = Ti == BigRational ? Rational{BigInt} : Ti

    Base.require_one_based_indexing(A)
    _I, _J, _V = findnz(A)
    I, J, V = issorted(_J) ? (_I, _J, _V) : begin
        _indices = sortperm(_J)
        @inbounds (_I[_indices], _J[_indices], _V[_indices])
    end
    # @inbounds if !issorted(_J)
    #     indices = sortperm(J)
    #     I = I[indices]; J = J[indices]; V = V[indices]
    # end
    isempty(J) && return compat_lu_convert(Val(Tf), A, 0, 0)
    m, n = size(A)
    minmn = min(m, n)
    if J[1] != 1 || I[1] != 1
        check && checknonsingular(1, Val(false)) # TODO update with Pivot
        # return LU{eltype(A), typeof(A)}(A, collect(1:minmn), convert(BlasInt, 1))
        return compat_lu_convert(Val(Tf), A, minmn, 1)
    end

    col_offset = zeros(Int, minmn) # for each col, index of the pivot element
    idx_cols = [[I[i] for i in getcolptr(A)[col+1]-1:-1:getcolptr(A)[col]] for col in 1:minmn]
    # For each column, indices of the non-zeros elements
    in_idx_colscol = falses(n)
    for col in 2:minmn
        sort!(idx_cols[col-1]; rev=true)
        # All idx_cols[x] are sorted by decreasing order for x < col
        # @show idx_cols[col]
        idx_colscol = idx_cols[col]
        in_idx_colscol[idx_colscol] .= true
        for row_j in idx_colscol
            row_j >= col && continue
            col_offset[col] += 1
            # @show idx_cols[row_j]
            idx_colsj = idx_cols[row_j]
            sizcol = length(idx_colscol)
            for row_i in idx_colsj
                if row_i ≤ row_j
                    break # Because the row_i are sorted in decreasing order
                end
                if !in_idx_colscol[row_i]
                    push!(idx_colscol, row_i)
                    in_idx_colscol[row_i] = true
                end
            end
            countadd = length(idx_colscol) - sizcol
            if countadd > 0
                siz = length(I)
                resize!(I, siz + countadd)
                resize!(J, siz + countadd)
                resize!(V, siz + countadd)
                for i in 1:countadd
                    row_i = idx_colscol[sizcol+i]
                    _idx = siz + i
                    J[_idx] = col
                    I[_idx] = row_i
                    V[_idx] = 0
                end
            end
        end
        in_idx_colscol[idx_colscol] .= false
    end
    B = sparse(I, J, Ti.(V)) # TODO update with Pivot
    # lu!(B, col_offset, check)
    rational_lu!(B, col_offset, check)
end

function forward_substitution!(L::SparseMatrixCSC, b)
    _, n = size(L)
    _, m = size(b)
    @inbounds for col in 1:n
        k = getcolptr(L)[col]
        rowvals(L)[k] == col || checknonsingular(col, Val(true))
        invnzLk = inv(nonzeros(L)[k])
        x = invnzLk .* b[col,:]
        b[col,:] .= x
        #=@simd=# for i in (k+1):getcolptr(L)[col+1]-1
            nzLi = nonzeros(L)[i]
            rvLi = rowvals(L)[i]
            @simd for j in 1:m
                b[rvLi,j] -= nzLi*x[j]
            end
        end
    end
    nothing
end

function backward_substitution!(U::SparseMatrixCSC, b)
    _, n = size(U)
    _, m = size(b)
    @inbounds for col in n:-1:1
        k = getcolptr(U)[col+1]-1
        rowvals(U)[k] == col || checknonsingular(col, Val(true))
        invnzUk = inv(nonzeros(U)[k])
        x = invnzUk .* b[col,:]
        b[col,:] .= x
        #=@simd=# for i in getcolptr(U)[col]:(k-1)
            nzUi = nonzeros(U)[i]
            rvUi = rowvals(U)[i]
            @simd for j in 1:m
                b[rvUi,j] -= nzUi*x[j]
            end
        end
    end
    nothing
end

function linsolve!(F::LU, B::Matrix)
    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
    BB = similar(B, TFB, size(B))
    copyto!(BB, B)
    m, n = size(F)
    minmn = min(m,n)
    L = tril!(getfield(F, :factors)[1:m, 1:minmn])
    for i = 1:minmn; L[i,i] = 1; end
    forward_substitution!(L, BB)
    x = triu!(getfield(F, :factors)[1:minmn, 1:n])
    backward_substitution!(x, BB)
    return BB
end

function rational_solve(::Val{N}, A::SparseMatrixCSC{Int,Int}, Y::Matrix{Int}) where N
    B = rational_lu(A, false)
    if !issuccess(B)
        error("Singular exception while equilibrating. Is the graph connected?")
    end
    Z = linsolve!(B, Rational{BigInt}.(Y))
    return hcat(zeros(Rational{Int128}, N), adjoint(Rational{Int128}.(Z)))
    # Rational{Int64} is not enough for tep for instance.
end

function copyuntil(j, oldZ, ::Type{T}) where T
    Z = similar(oldZ, T)
    for i in eachindex(Z)
        i == j && return Z
        Z[i] = oldZ[i]
    end
    error("Invalid failure of _inner_dixon_p!, please report this error")
    return Z # Does not matter
end

function _inner_dixon_p!(indices::Vector{Int}, Z::Matrix{Rational{T}}, h::BigInt,
                         x̄::Matrix{BigInt}, sqh::BigInt, tmp::BigInt) where T
    while !isempty(indices)
        j = pop!(indices)
        ua = MPZ.set(h)
        ub = deepcopy(@inbounds x̄[j])
        va = Int128(0)
        vb = Int128(1)
        k = 0
        while ub >= sqh
            k += 1
            # cpua = deepcopy(ua)
            # cpub = deepcopy(ub)
            MPZ.tdiv_qr!(tmp, ua, ua, ub)
            ua, ub = ub, ua
            # @toggleassert tmp == cpua ÷ cpub
            # @toggleassert ua == cpub
            # @toggleassert ub == cpua - tmp * cpub
            # cpuc = deepcopy(va)
            if typemin(Clong) < vb < typemax(Clong)
                MPZ.mul_si!(tmp, vb % Clong)
            else
                tmp *= vb
            end
            flag = signbit(va)
            va = abs(va)
            if va < typemax(Culong)
                if flag
                    MPZ.sub_ui!(tmp, va)
                else
                    MPZ.add_ui!(tmp, va)
                end
                va, vb = vb, Int128(tmp)
            else
                va, vb = vb, va + tmp
            end
            #= or replace all of the above since if typeomin(Clong) < ... by
            MPZ.mul!(tmp, vb)
            MPZ.add!(tmp, va)
            va, vb, tmp = vb, tmp, va
            =#
            # @toggleassert vb == cpuc + tmp * va
        end

        uv::Tuple{T,T} = if T === BigInt
            Base.divgcd(ub, vb)
        else
            ud, vd = Base.divgcd(ub, vb)
            m = typemin(T)
            M = typemax(T)
            if !(m < ud < M && m < vd < M)
                push!(indices, j)
                return false
            end
            (ud % T, vd % T)
        end

        @inbounds Z[j] = (-2*isodd(k)+1) * Base.checked_den(uv[1], uv[2])
        # @show Z[j]
        # @toggleassert mod((-1)^isodd(k) * ub, h) == mod(vb * x̄[j], h)
    end
    return true
end

function dixon_p(::Val{N}, A::SparseMatrixCSC{Int,Int}, C::LU{Modulo{p,Int32}}, Y::Matrix{Int}) where {N,p}
    λs = [norm(x) for x in eachcol(A)]
    append!(λs, [norm(x) for x in eachcol(Y)])
    partialsort!(λs, N)
    for _ in 1:N
        popfirst!(λs)
    end
    δ::BigFloat = prod(BigFloat, λs; init=one(BigFloat))
    # @show δ
    # @show p
    m = ceil(Int, 2*log(δ / (MathConstants.φ - 1))/log(p))
    @toggleassert m ≥ 1
    # @show m
    B = copy(Y)
    x̄ = BigInt.(linsolve!(C, B))
    X = copy(x̄)
    @toggleassert A * Modulo{p,Int32}.(X) == B
    h = one(BigInt) # = p^i
    tmp = BigInt()
    for i in 1:m-1
        MPZ.mul_si!(h, p)
        B .= (B .- A*Integer.(X)) .÷ p
        X .= Integer.(linsolve!(C, B))
        @toggleassert A * Modulo{p,Int32}.(X) == B
        # x̄ .+= h .* X
        @inbounds for j in eachindex(x̄)
            MPZ.mul!(tmp, X[j], h)
            MPZ.add!(x̄[j], tmp)
        end
    end
    MPZ.mul_si!(h, p) # h = p^m
    @toggleassert mod.(A * x̄, h) == mod.(Y, h)
    sqh = MPZ.sqrt(h) # h = p^{m/2}
    typeofZ = Union{Matrix{Rational{Int64}},Matrix{Rational{Int128}},Matrix{Rational{BigInt}}}
    Z::typeofZ = similar(Y, Rational{Int64})
    indices = collect(reverse(eachindex(Z)))
    success = _inner_dixon_p!(indices, Z, h, x̄, sqh, tmp)
    if !success
        Z = copyuntil(first(indices), Z, Rational{Int128})
        success = _inner_dixon_p!(indices, Z, h, x̄, sqh, tmp)
        if !success
            Z = copyuntil(first(indices), Z, Rational{BigInt})
            success = _inner_dixon_p!(indices, Z, h, x̄, sqh, tmp)
            @toggleassert success
        end
    end

    @toggleassert eltype(Y).(A * big.(Z)) == Y
    return Z::typeofZ
    # return hcat(zeros(Rational{Int128}, N), Rational{Int128}.(x̄)')
    # Rational{Int64} is not enough for tep for instance.
end

@static if VERSION < v"1.8-"
    const typeofB = Union{
        LU{Modulo{2147483647,Int32},SparseMatrixCSC{Modulo{2147483647,Int32},Int}},
        LU{Modulo{2147483629,Int32},SparseMatrixCSC{Modulo{2147483629,Int32},Int}},
        LU{Modulo{2147483587,Int32},SparseMatrixCSC{Modulo{2147483587,Int32},Int}}
    }
else
    const typeofB = Union{
        LU{Modulo{2147483647,Int32},SparseMatrixCSC{Modulo{2147483647,Int32},Int},OneTo{Int}},
        LU{Modulo{2147483629,Int32},SparseMatrixCSC{Modulo{2147483629,Int32},Int},OneTo{Int}},
        LU{Modulo{2147483587,Int32},SparseMatrixCSC{Modulo{2147483587,Int32},Int},OneTo{Int}}
    }
end

const typeofZ = Union{Matrix{Rational{Int64}},Matrix{Rational{Int128}},Matrix{Rational{BigInt}}}

"""
    dixon_solve(::Val{N}, A, Y) where N

Specialized solver for the linear system `A*X = Y` where `A` is a sparse symmetric
integer `n×n` matrix and `Y` is a dense integer `n×N` matrix, using Dixon's method.

Return `X` as a matrix of `Rational{Int128}`.
"""
function dixon_solve(::Val{N}, A::SparseMatrixCSC{Int,Int}, Y::Matrix{Int}) where N
    # @show time_ns()
    B::typeofB = rational_lu(A, false, Modulo{2147483647,Int32})
    if issuccess(B)
        Z = dixon_p(Val(N), A, B, Y)'
    else
        B = rational_lu(A, false, Modulo{2147483629,Int32})
        if issuccess(B)
            Z = dixon_p(Val(N), A, B, Y)'
        else
            B = rational_lu(A, false, Modulo{2147483587,Int32})
            if issuccess(B)
                Z = dixon_p(Val(N), A, B, Y)'
            else # The probability of this being required is *extremely* low
                return rational_solve(Val(N), A, Y)
            end
        end
    end
    return hcat(zeros(eltype(Z), N), Z)::typeofZ
end


text:
#
# 1. Types supporting parameterized Timestep and Clock objects
#

abstract type AbstractTimestep end

struct FixedTimestep{FIRST, STEP, LAST} <: AbstractTimestep
    t::Int
end

struct VariableTimestep{TIMES} <: AbstractTimestep
    t::Int
    current::Int

    function VariableTimestep{TIMES}(t::Int = 1) where {TIMES}
        # The special case below handles when functions like next_step step beyond
        # the end of the TIMES array.  The assumption is that the length of this
        # last timestep, starting at TIMES[end], is 1.
        current::Int = t > length(TIMES) ? TIMES[end] + 1 : TIMES[t]

        return new(t, current)
    end
end

mutable struct Clock{T <: AbstractTimestep}
	ts::T

	function Clock{T}(FIRST::Int, STEP::Int, LAST::Int) where T
		return new(FixedTimestep{FIRST, STEP, LAST}(1))
    end

    function Clock{T}(TIMES::NTuple{N, Int} where N) where T
        return new(VariableTimestep{TIMES}())
    end
end

mutable struct TimestepArray{T_TS <: AbstractTimestep, T, N}
	data::Array{T, N}

    function TimestepArray{T_TS, T, N}(d::Array{T, N}) where {T_TS, T, N}
		return new(d)
	end

    function TimestepArray{T_TS, T, N}(lengths::Int...) where {T_TS, T, N}
		return new(Array{T, N}(undef, lengths...))
	end
end

# Since these are the most common cases, we define methods (in time.jl)
# specific to these type aliases, avoiding some of the inefficiencies
# associated with an arbitrary number of dimensions.
const TimestepMatrix{T_TS, T} = TimestepArray{T_TS, T, 2}
const TimestepVector{T_TS, T} = TimestepArray{T_TS, T, 1}

#
# 2. Dimensions
#

abstract type AbstractDimension end

const DimensionKeyTypes   = Union{AbstractString, Symbol, Int, Float64}
const DimensionRangeTypes = Union{UnitRange{Int}, StepRange{Int, Int}}

struct Dimension{T <: DimensionKeyTypes} <: AbstractDimension
    dict::OrderedDict{T, Int}

    function Dimension(keys::Vector{T}) where {T <: DimensionKeyTypes}
        dict = OrderedDict(collect(zip(keys, 1:length(keys))))
        return new{T}(dict)
    end

    function Dimension(rng::T) where {T <: DimensionRangeTypes}
        return Dimension(collect(rng))
    end

    Dimension(i::Int) = Dimension(1:i)

    # Support Dimension(:foo, :bar, :baz)
    function Dimension(keys::T...) where {T <: DimensionKeyTypes}
        vector = [key for key in keys]
        return Dimension(vector)
    end
end

#
# Simple optimization for ranges since indices are computable.
# Unclear whether this is really any better than simply using
# a dict for all cases. Might scrap this in the end.
#
mutable struct RangeDimension{T <: DimensionRangeTypes} <: AbstractDimension
    range::T
 end

#
# 3. Types supporting Parameters and their connections
#
abstract type ModelParameter end

# TBD: rename ScalarParameter, ArrayParameter, and AbstractParameter?

mutable struct ScalarModelParameter{T} <: ModelParameter
    value::T

    function ScalarModelParameter{T}(value::T) where T
        new(value)
    end

    function ScalarModelParameter{T1}(value::T2) where {T1, T2}
        try
            new(T1(value))
        catch err
            error("Failed to convert $value::$T2 to $T1")
        end
    end
end

mutable struct ArrayModelParameter{T} <: ModelParameter
    values::T
    dimensions::Vector{Symbol} # if empty, we don't have the dimensions' name information

    function ArrayModelParameter{T}(values::T, dims::Vector{Symbol}) where T
        new(values, dims)
    end
end

ScalarModelParameter(value) = ScalarModelParameter{typeof(value)}(value)

Base.convert(::Type{ScalarModelParameter{T}}, value::Number) where {T} = ScalarModelParameter{T}(T(value))

Base.convert(::Type{T}, s::ScalarModelParameter{T}) where {T} = T(s.value)

ArrayModelParameter(value, dims::Vector{Symbol}) = ArrayModelParameter{typeof(value)}(value, dims)

abstract type AbstractConnection end

struct InternalParameterConnection <: AbstractConnection
    src_comp_name::Symbol
    src_var_name::Symbol
    dst_comp_name::Symbol
    dst_par_name::Symbol
    ignoreunits::Bool
    backup::Union{Symbol, Nothing} # a Symbol identifying the external param providing backup data, or nothing
    offset::Int

    function InternalParameterConnection(src_comp::Symbol, src_var::Symbol, dst_comp::Symbol, dst_par::Symbol,
                                         ignoreunits::Bool, backup::Union{Symbol, Nothing}=nothing; offset::Int=0)
        self = new(src_comp, src_var, dst_comp, dst_par, ignoreunits, backup, offset)
        return self
    end
end

struct ExternalParameterConnection  <: AbstractConnection
    comp_name::Symbol
    param_name::Symbol      # name of the parameter in the component
    external_param::Symbol  # name of the parameter stored in md.external_params
end

#
# 4. Types supporting structural definition of models and their components
#

# To identify components, we create a variable with the name of the component
# whose value is an instance of this type, e.g.
# const global adder = ComponentId(module_name, comp_name)
struct ComponentId
    module_name::Symbol
    comp_name::Symbol
end

# Indicates that the object has a `name` attribute
abstract type NamedDef end

# Supertype for vars and params
# abstract type DatumDef <: NamedDef end

# The same structure is used for variables and parameters
mutable struct DatumDef <: NamedDef
    name::Symbol
    datatype::DataType
    dimensions::Vector{Symbol}
    description::String
    unit::String
    datum_type::Symbol          # :parameter or :variable
    default::Any                # used only for Parameters

    function DatumDef(name::Symbol, datatype::DataType, dimensions::Vector{Symbol},
                      description::String, unit::String, datum_type::Symbol,
                      default::Any=nothing)
        self = new()
        self.name = name
        self.datatype = datatype
        self.dimensions = dimensions
        self.description = description
        self.unit = unit
        self.datum_type = datum_type
        self.default = default
        return self
    end

end

mutable struct ComponentDef  <: NamedDef
    name::Symbol
    comp_id::ComponentId
    variables::OrderedDict{Symbol, DatumDef}
    parameters::OrderedDict{Symbol, DatumDef}
    dimensions::OrderedDict{Symbol, Union{Nothing, Dimension}}
    first::Union{Nothing, Int}
    last::Union{Nothing, Int}

    # ComponentDefs are created "empty"; elements are subsequently added
    # to them via addvariable, add_dimension!, etc.
    function ComponentDef(comp_id::ComponentId)
        self = new()
        self.name = comp_id.comp_name
        self.comp_id = comp_id
        self.variables  = OrderedDict{Symbol, DatumDef}()
        self.parameters = OrderedDict{Symbol, DatumDef}()
        self.dimensions = OrderedDict{Symbol, Union{Nothing, Dimension}}()
        self.first = self.last = nothing
        return self
    end
end

# Declarative definition of a model, used to create a ModelInstance
mutable struct ModelDef
    module_name::Symbol     # the module in which this model was defined

    # Components keyed by symbolic name, allowing a given component
    # to occur multiple times within a model.
    comp_defs::OrderedDict{Symbol, ComponentDef}

    dimensions::Dict{Symbol, Dimension}

    number_type::DataType

    internal_param_conns::Vector{InternalParameterConnection}
    external_param_conns::Vector{ExternalParameterConnection}

    # Names of external params that the ConnectorComps will use as their :input2 parameters.
    backups::Vector{Symbol}

    external_params::Dict{Symbol, ModelParameter}

    sorted_comps::Union{Nothing, Vector{Symbol}}

    is_uniform::Bool

    function ModelDef(number_type=Float64)
        self = new()
        self.module_name = nameof(@__MODULE__)                  # TBD: fix this; should by module model is defined in
        self.comp_defs = OrderedDict{Symbol, ComponentDef}()
        self.dimensions = Dict{Symbol, Dimension}()
        self.number_type = number_type
        self.internal_param_conns = Vector{InternalParameterConnection}()
        self.external_param_conns = Vector{ExternalParameterConnection}()
        self.external_params = Dict{Symbol, ModelParameter}()
        self.backups = Vector{Symbol}()
        self.sorted_comps = nothing
        self.is_uniform = true
        return self
    end
end

#
# 5. Types supporting instantiated models and their components
#

# Supertype for variables and parameters in component instances
abstract type ComponentInstanceData end

struct ComponentInstanceParameters{NT <: NamedTuple} <: ComponentInstanceData
    nt::NT

    function ComponentInstanceParameters{NT}(nt::NT) where {NT <: NamedTuple}
        return new{NT}(nt)
    end
end

function ComponentInstanceParameters(names, types, values)
    NT = NamedTuple{names, types}
    ComponentInstanceParameters{NT}(NT(values))
end

function ComponentInstanceParameters{NT}(values::T) where {NT <: NamedTuple, T <: AbstractArray}
    ComponentInstanceParameters{NT}(NT(values))
end

struct ComponentInstanceVariables{NT <: NamedTuple} <: ComponentInstanceData
    nt::NT

    function ComponentInstanceVariables{NT}(nt::NT) where {NT <: NamedTuple}
        return new{NT}(nt)
    end
end

function ComponentInstanceVariables{NT}(values::T) where {NT <: NamedTuple, T <: AbstractArray}
    ComponentInstanceVariables{NT}(NT(values))
end

function ComponentInstanceVariables(names, types, values)
    NT = NamedTuple{names, types}
    ComponentInstanceVariables{NT}(NT(values))
end

# A container class that wraps the dimension dictionary when passed to run_timestep()
# and init(), so we can safely implement Base.getproperty(), allowing `d.regions` etc.
struct DimDict
    dict::Dict{Symbol, Vector{Int}}
end

# Special case support for Dicts so we can use dot notation on dimension.
# The run_timestep() and init() funcs pass a DimDict of dimensions by name
# as the "d" parameter.
@inline function Base.getproperty(dimdict::DimDict, property::Symbol)
    return getfield(dimdict, :dict)[property]
end

# TBD: try with out where clause, i.e., just obj::ComponentInstanceData
nt(obj::T)  where {T <: ComponentInstanceData} = getfield(obj, :nt)
Base.names(obj::T)  where {T <: ComponentInstanceData} = keys(nt(obj))
Base.values(obj::T) where {T <: ComponentInstanceData} = values(nt(obj))
types(obj::T) where {T <: ComponentInstanceData} = typeof(nt(obj)).parameters[2].parameters

# An instance of this type is passed to the run_timestep function of a
# component, typically as the `p` argument. The main role of this type
# is to provide the convenient `p.nameofparameter` syntax.
mutable struct ComponentInstance{TV <: ComponentInstanceVariables, TP <: ComponentInstanceParameters}
    comp_name::Symbol
    comp_id::ComponentId
    variables::TV
    parameters::TP
    dim_dict::Dict{Symbol, Vector{Int}}

    first::Int
    last::Int

    init::Union{Nothing, Function}
    run_timestep::Union{Nothing, Function}

    function ComponentInstance{TV, TP}(comp_def::ComponentDef,
                               vars::TV, pars::TP,
                               first::Int, last::Int,
                               name::Symbol=name(comp_def)) where {TV <: ComponentInstanceVariables,
                                                                   TP <: ComponentInstanceParameters}
        self = new{TV, TP}()

        self.comp_id = comp_id = comp_def.comp_id
        self.comp_name = name
        self.dim_dict = Dict{Symbol, Vector{Int}}()    # set in "build" stage

        self.variables = vars
        self.parameters = pars
        self.first = first
        self.last = last

        comp_name   = comp_id.comp_name
        module_name = comp_id.module_name
        comp_module = Base.eval(Main, module_name)

        # TBD: use FunctionWrapper here?
        function get_func(name)
            func_name = Symbol("$(name)_$(comp_name)")
            try
                Base.eval(comp_module, func_name)
            catch err
                # No need to warn about this...
                # @warn "Failed to evaluate function name $func_name in module $comp_module"
                nothing
            end
        end

        self.init = get_func("init")
        self.run_timestep = get_func("run_timestep")

        return self
    end
end

# This type holds the values of a built model and can actually be run.
mutable struct ModelInstance
    md::ModelDef

    # Ordered list of components (including hidden ConnectorComps)
    components::OrderedDict{Symbol, ComponentInstance}

    firsts::Vector{Int}        # in order corresponding with components
    lasts::Vector{Int}

    function ModelInstance(md::ModelDef)
        self = new()
        self.md = md
        self.components = OrderedDict{Symbol, ComponentInstance}()
        self.firsts = Vector{Int}()
        self.lasts = Vector{Int}()
        return self
    end
end

#
# 6. User-facing Model types providing a simplified API to model definitions and instances.
#

"""
    Model

A user-facing API containing a `ModelInstance` (`mi`) and a `ModelDef` (`md`).
This `Model` can be created with the optional keyword argument `number_type` indicating
the default type of number used for the `ModelDef`.  If not specified the `Model` assumes
a `number_type` of `Float64`.
"""
mutable struct Model
    md::ModelDef
    mi::Union{Nothing, ModelInstance}

    function Model(number_type::DataType=Float64)
        return new(ModelDef(number_type), nothing)
    end

    # Create a copy of a model, e.g., to create marginal models
    function Model(m::Model)
        return new(copy(m.md), nothing)
    end
end

"""
    MarginalModel

A Mimi `Model` whose results are obtained by subtracting results of one `base` Model
from those of another `marginal` Model` that has a difference of `delta`.
"""
struct MarginalModel
    base::Model
    marginal::Model
    delta::Float64

    function MarginalModel(base::Model, delta::Float64=1.0)
        return new(base, Model(base), delta)
    end
end

function Base.getindex(mm::MarginalModel, comp_name::Symbol, name::Symbol)
    return (mm.marginal[comp_name, name] .- mm.base[comp_name, name]) ./ mm.delta
end

#
# 7. Reference types provide more convenient syntax for interrogating Components
#

"""
    ComponentReference

A container for a component, for interacting with it within a model.
"""
struct ComponentReference
    model::Model
    comp_name::Symbol
end

"""
    VariableReference

A container for a variable within a component, to improve connect_param! aesthetics,
by supporting subscripting notation via getindex & setindex .
"""
struct VariableReference
    model::Model
    comp_name::Symbol
    var_name::Symbol
end


text:
function atmosphere(alt)

    R0 = 2.377e-3  # Sea level density  (slug/ft³)

    # alt (ft)
    tfac = 1.0 - 0.703e-5 * alt
    # temperature (Rankine)
    alt >= 35000.0 ? T = 390.0 : T = 519.0 * tfac
    # density (slug / ft^3)
    ρ = R0 * tfac^4.14
    # sound velocity (ft/s)
    a = sqrt(1.4 * 1716.3 * T)
    # pressure (psf)
    p = 1715.0 * ρ * T
    return T, ρ, a, p
end

function atmosphere_f16(alt)
    T, ρ, a, p = atmosphere(alt * M2FT)
    return T * RANK2KEL, ρ * SLUGFT32KGM3, a * FT2M, p * PSF2PA
end


text:
mutable struct Body
    vertices

    position
    velocity

    angle  # orientation
    angular_velocity

    force
    torque

    restitution
    friction

    mass
    inertia
end

"""
    translate!(body::Body, translation::Vector)

Translate a body by a given vector without imparting velocity.
"""
function translate!(body::Body, translation::Vector)
    body.position += translation
end

"""
    rotate!(body::Body, rotation)

Rotate a body by a given angle without imparting angular velocity.
"""
function rotate!(body::Body, rotation)
    body.angle += rotation
end

"""
    update!(body::Body, dt)

Update the position, velocity, angle, and angular velocity of the body.
"""
function update!(body::Body, dt)

    # update velocity
    body.velocity.x += (body.force.x / body.mass) * dt
    body.velocity.y += (body.force.y / body.mass) * dt

    # TODO: add friction/drag

    # update position
    body.position.x += body.velocity.x
    body.position.y += body.velocity.y

    # update angular velocity
    body.angular_velocity += (body.torque / body.inertia) * dt

    # TODO: add friction/drag

    # update angle
    body.angle += body.angular_velocity

    # update vertices
    translate!(body.vertices, body.velocity)
    rotate!(body.vertices, body.angular_velocity)
end

"""
    apply_force!(body::Body, force, position)

Apply a force to a body at a given position.
"""
function apply_force!(b::Body, force, position)
    body.force += force
    offset = position - body.position
    body.torque += offset × force  # "cross product" = "determinant"
end