text:
"""
Resample curve

    resamplecurve(x, N=100)
    :param x: array describing curve (n,T)
    :param N: Number of samples to re-sample curve, N usually is > T
"""
function resamplecurve(x, N=100)
    n, T = size(x);
    xn = zeros(n, N);

    delta = zeros(T);

    for r in 2:T
        delta[r] = norm(x[:,r]-x[:,r-1]);
    end

    cumdel = cumsum(delta)/sum(delta);
    newdel = linspace(0,1,N);

    for r in 1:n
        s = Spline1D(cumdel,vec(x[r,:]));
        xn[r,:] = evaluate(s,newdel);
    end

    return xn
end

"""
Calculate centroid of curve

    calculatecentroid(beta)
    :param beta: array describing curve (n,T)
"""
function calculatecentroid(beta)
    n, T = size(beta);
    betadot = zeros(n,T);
    for i = 1:n
        betadot[i,:] = gradient(vec(beta[i,:]), 1.0/(T-1));
    end
    normbetadot = zeros(T);
    integrand = zeros(n,T);
    for i = 1:T
        normbetadot[i] = norm(betadot[:,i]);
        integrand[:,i] = beta[:,i]*normbetadot[i];
    end
    scale = trapz(collect(linspace(0,1,T)), normbetadot);
    centroid = trapz(collect(linspace(0,1,T)), integrand, 2)./scale;

    return centroid
end

"""
Convert curve to square-root velocity function (srvf)

    curve_to_q(beta)
    :param beta: array describing curve (n,T)
"""
function curve_to_q(beta, wscale::Bool=false)
    n, T = size(beta);
    v = zeros(n,T);
    for i = 1:n
        v[i,:] = gradient(vec(beta[i,:]), 1.0/(T-1));
    end

    len = sum(sqrt(sum(v.*v)))/T;
    v = v./len;
    q = zeros(n,T);
    for i = 1:T
        L = sqrt(norm(v[:,i]));
        if (L > 0.0001)
            q[:,i] = v[:,i]./L;
        else
            q[:,i] = v[:,i]*0.0001;
        end
    end

    if (!wscale)
        q = q./sqrt(innerprod_q2(q, q));
    end

    return q
end

"""
Convert srvf to curve

    q_to_curve(q)
    :param q: array describing srvf (n,T)
"""
function q_to_curve(q)
    T = size(q,2);
    qnorm = zeros(T);
    for i = 1:T
        qnorm[i] = norm(q[:,i]);
    end

    integrand = zeros(2,T);
    integrand[1,:] = vec(q[1,:]).*qnorm;
    integrand[2,:] = vec(q[2,:]).*qnorm;
    beta = cumtrapz(collect(1.:T), integrand, 2)./T;

    return beta
end

function optimum_reparam(beta1::Array{Float64,2}, beta2::Array{Float64,2},
                         lam::Float64=0.0; method::String="DP", w=0.01,
                         rotated::Bool=true, isclosed::Bool=false)
    n1, M = size(beta2);
    timet = collect(linspace(0,1,M));
    skipm = 4;
    auto = 2;
    tau = 0;
    if (method == "DP")
        # Optimze over SO(n) x Gamma
        q1 = curve_to_q(beta1)

        # Optimzie over SO(n)
        beta2, R, tau = find_rotation_seed_coord(beta1, beta2);
        q2 = curve_to_q(beta2);

        # Optimzie over Gamma
        q1i = vec(reshape(q1, M*n1, 1));
        q2i = vec(reshape(q2, M*n1, 1));
        G = zeros(M);
        T = zeros(M);
        sizei = Cdouble[0];
        ccall((:DynamicProgrammingQ2, libfdasrsf), Void,
            (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int32,
            Int32, Int32, Ptr{Float64},Ptr{Float64}, Int32, Int32,
            Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Float64), q2i, timet,
            q1i, timet, n1, M, M, timet, timet, M, M, G, T, sizei, lam)
        G = G[1:round(Integer,sizei[1])];
        T = T[1:round(Integer,sizei[1])];
        yi = interpolate((T,), G, Gridded(Linear()))
        gam = yi[timet];

    elseif (method == "DP2")
        c1 = vec(reshape(beta1', M*n1, 1));
        c2 = vec(reshape(beta2', M*n1, 1));
        opt = zeros(M+n1*n1+1);
        swap = false;
        fopts = zeros(5);
        comtime = zeros(5);
        @cpp ccall((:optimum_reparam, libgropt), Void,
                   (Ptr{Float64}, Ptr{Float64}, Int32, Int32, Float64, Bool,
                    Bool, Bool, Int32, Int32, Ptr{Float64}, Bool,
                    Ptr{Float64}, Ptr{Float64}), c1, c2, M, n1, 0.0, true,
                    rotated, isclosed, skipm, auto, opt, swap, fopts,
                    comtime)

        gam = opt[1:end-5];
        R = reshape(opt[end-4:end-1],2,2);

        if swap
            gam = invertGamma(gam);
            R = R';
        end

    else
        c1 = vec(reshape(beta1', M*n1, 1));
        c2 = vec(reshape(beta2', M*n1, 1));
        opt = zeros(M+n1*n1+1);
        swap = false;
        fopts = zeros(5);
        comtime = zeros(5);
        @cpp ccall((:optimum_reparam, libgropt), Void,
                   (Ptr{Float64}, Ptr{Float64}, Int32, Int32, Float64, Bool,
                    Ptr{Float64}, Bool, Bool, Bool, Int32, Int32, Ptr{Float64},
                    Ptr{Float64}), c1, c2, M, n1, w, false, rotated, isclosed,
                    skipm, auto, opt, swap, fopts, comtime)

        if fopts[1] == 1000
            @cpp ccall((:optimum_reparam, libgropt), Void,
                       (Ptr{Float64}, Ptr{Float64}, Int32, Int32, Float64, Bool,
                        Ptr{Float64}, Bool, Bool, Bool, Int32, Int32,
                        Ptr{Float64}, Ptr{Float64}), c1, c2, M, n1, 0.0, true,
                        rotated, isclosed, skipm, auto, opt, swap, fopts,
                        comtime)
        end

        gam = opt[1:end-5];
        R = reshape(opt[end-4:end-1],2,2);

        if swap
            gam = invertGamma(gam);
            R = R';
        end
    end

    gam = (gam-gam[1]) ./ (gam[end] - gam[1]);

    return gam, R, tau
end

"""
Find best rotation between two srvfs

    find_best_roation(q1,q2)
    :param q1: array (n, T)
    :param q2: array (n, T)
"""
function find_best_rotation(q1, q2)
    epsilon = eps(Float64);
    n, T = size(q1);
    A = q1*q2';
    U, S, V = svd(A);
    if (abs(det(U) * det(V) - 1) < 10 * epsilon)
        S = eye(n)
    else
        S = eye(n)
        S[:,end] = -S[:,end];
    end
    R = U*S*V';
    q2new = R*q2;

    return q2new, R
end

"""
Calculate variance along curve

    calculate_variance(beta)
    :param beta: array (n, T)
"""
function calculate_variance(beta)
    n, T = size(beta);
    betadot = zeros(n, T);
    for i = 1:n
        betadot[i,:] = gradient(vec(beta[i,:]), 1.0/(T-1));
    end
    normbetadot = zeros(T);
    centroid = calculatecentroid(beta);
    integrand = zeros(n, n, T);
    t = collect(linspace(0,1,T));
    for i = 1:T
        normbetadot[i] = norm(betadot[:,i]);
        a1 = beta[:,i] - centroid;
        integrand[:, :, i] = a1 * a1' * normbetadot[i];
    end
    l = trapz(t, normbetadot);
    variance = trapz(t, integrand, 3);
    variance /= l;

    return variance
end

"""
Calculate psi

    psi(x,a,q)
    :param x: array (n, T) of curve
    :param a: vector (n) of centroid
    :param q: array (n, T) of srvf
"""
function psi(x, a, q)
    T = size(q,2);
    covmat = calculate_variance(x+repmat(a,1,T));
    psi1 = covmat[1,1] - covmat[2,2];
    psi2 = covmat[1,2];
    psi3 = x[1,end];
    psi4 = x[2,end];

    return psi1, psi2, psi3, psi4
end

"""
Find basis normal to srvf

    find_basis_normal(q)
    :param q: array (n, T)
"""
function find_basis_normal(q)
    n, T = size(q);

    f1 = zeros(n,T);
    f2 = zeros(n,T);
    for i = 1:T
        f1[:,i] = q[1,i]*q[:,i]/norm(q[:,i])+[norm(q[:,i]);0];
        f2[:,i] = q[2,i]*q[:,i]/norm(q[:,i])+[0;norm(q[:,i])];
    end
    h3 = copy(f1);
    h4 = copy(f2);
    integrandb3 = zeros(T);
    integrandb4 = zeros(T);
    for i = 1:T
        integrandb3[i] = (q[:,i]'*h3[:,i])[1];
        integrandb4[i] = (q[:,i]'*h4[:,i])[1];
    end
    b3 = h3 - q*trapz(collect(linspace(0,1,T)),integrandb3);
    b4 = h4 - q*trapz(collect(linspace(0,1,T)),integrandb4);

    basis = Array(Any,2);
    basis[1] = b3;
    basis[2] = b4;

    return basis
end

"""
Calculate Jacobian of basis

    calc_j(basis)
    :param basis: array of (2) from find_basis_normal
"""
function calc_j(basis)
    b1 = basis[1];
    b2 = basis[2];
    T = size(b1,2);

    integrand11 = zeros(T);
    integrand12 = zeros(T);
    integrand22 = zeros(T);

    for i = 1:T
        integrand11[i] = (b1[:,i]'*b1[:,i])[1];
        integrand12[i] = (b1[:,i]'*b2[:,i])[1];
        integrand22[i] = (b2[:,i]'*b2[:,i])[1];
    end

    j = zeros(2,2);
    j[1,1] = trapz(collect(linspace(0,1,T)), integrand11);
    j[1,2] = trapz(collect(linspace(0,1,T)), integrand12);
    j[2,2] = trapz(collect(linspace(0,1,T)), integrand22);
    j[2,1] = j[1,2];

    return j
end

"""
Shift curve f by tau

    shift_f(f, tau)
    :param f: array (n,T)
    :param tau: scalar
"""
function shift_f(f, tau)
    n, T = size(f);
    fn = zeros(n,T);
    fn[:,1:T-1] = circshift(f[:,1:T-1],[0 tau]);
    fn[:,T] = fn[:,1];

    return fn
end

"""
Finds best rotation and seed value between two curves beta1, beta2

    find_rotation_seed_coord(beta1,beta2)
    :param beta1: array (n,T)
    :param beta2: array (n,T)

    :return beta2new: rotated curve
    :return O_hat1: rotation matrix
    :return tau: seed value
"""
function find_rotation_seed_coord(beta1,beta2)
    n,T = size(beta1);
    q1 = curve_to_q(beta1);
    Ltwo = zeros(T);
    Rlist = zeros(n, n, T);
    for ctr = 1:T
        beta2n = shift_f(beta2, ctr);
        beta2new, R = find_best_rotation(beta1, beta2n);
        q2new = curve_to_q(beta2new);
        Ltwo[ctr] = innerprod_q2(q1-q2new, q1-q2new);
        Rlist[:,:,ctr] = R;
    end

    tau = indmin(Ltwo);
    O_hat = Rlist[:,:,tau];
    beta2new = shift_f(beta2, tau);
    beta2new = O_hat * beta2new;

    return beta2new, O_hat, tau
end

"""
Find best rotation and seed value between two srvfs, q1 and q2

    find_rotation_and_seed_q(q1,q2)
    :param q1: array (n,T)
    :param q2: array (n,T)

    :return q2new: rotated srvf
    :return O_hat: rotation matrix
    :return tau: seed value
"""
function find_rotation_and_seed_q(q1,q2)
    n, T = size(q1);
    Ltwo = zeros(T);
    Rlist = zeros(n, n, T);
    for ctr = 1:T
        q2n = shift_f(q2, ctr);
        q2new, R = find_best_rotation(q1, q2n);
        Ltwo[ctr] = innerprod_q2(q1-q2new, q1-q2new);
        Rlist[:, :, ctr] = R;
    end

    tau = indmin(Ltwo);
    O_hat = Rlist[:,:,tau];
    q2new = shift_f(q2,tau);
    q2new = O_hat * q2new;

    return q2new, O_hat, tau
end

"""
Warp srvf by gamma

    group_action_by_gamma(q, gamma)
    :param q: array (n,T)
    :param gamma: vector (T)
"""
function group_action_by_gamma(q, gamma)
    n, T = size(q);
    gammadot = gradient(gamma, 1.0/T);
    qn = zeros(n, T);
    timet = collect(linspace(0, 1, T));

    for j = 1:n
        s = Spline1D(timet, vec(q[j, :]));
        qn[j, :] = evaluate(s, gamma) .* sqrt(gammadot);
    end

    qn /= sqrt(innerprod_q2(qn,qn));

    return qn
end

"""
Warp curve f by gamma

    group_action_by_gamma_coord(f, gamma)
    :param f: array (n,T)
    :param gamma: vector (T)
"""
function group_action_by_gamma_coord(f, gamma)
    n, T = size(f);
    fn = zeros(n,T);
    timet = collect(linspace(0,1,T));

    for j = 1:n
        s = Spline1D(timet, vec(f[j, :]));
        fn[j,:] = evaluate(s, gamma);
    end

    return fn
end

"""
Project curve onto normal tangent space

    project_curve(q)
    :param q: array (n,T)
"""
function project_curve(q)
    T = size(q,2);
    tol = 1e-5;
    maxit = 200;
    itr = 0;
    delta = 0.5;
    x = q_to_curve(q);
    a = -1 * calculatecentroid(x);

    psi1, psi2, psi3, psi4 = psi(x,a,q);

    r = [psi3,psi4];
    rnorm = zeros(maxit+1);
    rnorm[1] = norm(r);

    while itr <= maxit
        basis = find_basis_normal(q);

        # calculate Jacobian
        J = calc_j(basis);

        # Newton-Raphson step to update q
        y = J\(-r);
        dq = delta * (y[1]*basis[1] + y[2]*basis[2]);
        normdq = sqrt(innerprod_q2(dq,dq));
        q = cos(normdq) * q + sin(normdq) * dq/normdq;
        q /= sqrt(innerprod_q2(q,q));

        # update x and a from the new q
        beta_new = q_to_curve(q);
        x = copy(beta_new);
        a = -1 * calculatecentroid(x);
        beta_new = x + repmat(a,1,T);

        # calculate the new value of psi
        psi1, psi2, psi3, psi4 = psi(x,a,q);
        r = [psi3,psi4];
        rnorm[itr+1] = norm(r);

        if norm(r) < tol
            break
        end

        itr += 1;

    end

    rnorm = rnorm[1:itr-1];

    return q
end

"""
Pre-processes curve by centering and projecting to tangent space

    pre_proc_curve(beta, T=100)
    :param beta: array (n,T)
    :param T: number of re-sample points

    :return betanew: re-sampled and centered curve
    :return qnew: projected srvf
    :return A: rotation matrix
"""
function pre_proc_curve(beta, T=100)
    beta = resamplecurve(beta, T);
    q = curve_to_q(beta);
    qnew = project_curve(q);
    x = q_to_curve(qnew);
    a = -1 * calculatecentroid(x);
    betanew = x + repmat(a,1,T);
    A = eye(2)

    return betanew, qnew, A
end

"""
Calculate shooting vector and distance between two curves beta1 and beta2

    inverse_exp_coord(beta1, beta2)
    :param beta1: array (n,T)
    :param beta2: array (n,T)
    :param method: optimization method to find warping, default is
                   Dynamic Programming ("DP"). Other options are
                   Coordinate Descent ("DP2"), Riemannian BFGS
                   ("RBFGS")
    :param wscale: keep scale (false)

    :return v: shooting vector
    :return dist: shape distance
"""
function inverse_exp_coord(beta1, beta2; method::AbstractString="DP",
                           wscale=false)
    T = size(beta1,2);
    centroid1 = calculatecentroid(beta1);
    beta1 -= repmat(centroid1,1,T);
    centroid2 = calculatecentroid(beta2);
    beta2 -= repmat(centroid2,1,T);

    q1 = curve_to_q(beta1);

    # Iteratively optimize over SO(n) x Gamma using old DP
    gam, R, tau = optimum_reparam(beta1, beta2, method=method);
    beta2 = R * shift_f(beta2, tau);
    gamI = invert_gamma(gam);
    beta2 = group_action_by_gamma_coord(beta2, gamI);
    beta2, R, tau = find_rotation_seed_coord(beta1, beta2);
    q2n = curve_to_q(beta2);

    # Compute geodesic distance
    q1dotq2 = innerprod_q2(q1, q2n);
    if wscale
        dist = sqrt(innerprod_q2(q1-q2n, q1-q2n))
    else
        dist = acos(q1dotq2);
    end

    # Compute shooting vector
    if q1dotq2 > 1
        q1dotq2 = 1.;
    end

    u = q2n - q1dotq2 * q1;
    normu = sqrt(innerprod_q2(u, u));

    if normu > 1e-4
        v = u * acos(q1dotq2)/normu;
        if wscale
            v = q2n - q1
        end
    else
        v = zeros(2, T);
    end

    return v, dist
end

"""
Calculate shooting vector between two srvfs q1 and q2

    inverse_exp(q1, q2, beta2)
    :param q1: array (n,T)
    :param q2: array (n,T)
    :param beta2: array (n,T)

    :return v: shooting vector
"""
function inverse_exp(q1, q2, beta2)
    T = size(q1,2);
    centroid1 = calculatecentroid(beta2);
    beta2 -= repmat(centroid1,1,T);

    # Optimize over SO(n) x Gamma
    beta1 = q_to_curve(q1);
    gam, R, tau = optimum_reparam(beta1, beta2);
    gamI = invert_gamma(gam);
    beta2 = R * shift_f(beta2, tau);

    # Applying optimal re-parameterization to the second curve
    beta2 = group_action_by_gamma_coord(beta2, gamI);
    q2 = curve_to_q(beta2);

    # Optimize over SO(n)
    q2, R2, tau = find_rotation_and_seed_q(q1,q2);

    # Compute geodesic distance
    q1dotq2 = innerprod_q2(q1, q2);
    dist = acos(q1dotq2);

    # Compute shooting vector
    if q1dotq2 > 1
        q1dotq2 = 1.;
    end

    u = q2 - q1dotq2 * q1;
    normu = sqrt(innerprod_q2(u, u));

    if normu > 1e-4
        v = u * acos(q1dotq2) / normu;
    else
        v = zeros(2, T);
    end

    return v
end

"""
Gram-Schmidt Orthogonalization of basis

    gram_schmidt(basis)
    :param basis: basis out of find_basis_normal
"""
function gram_schmidt(basis)
    b1 = basis[1];
    b2 = basis[2];

    basis1 = b1 / sqrt(innerprod_q2(b1, b1));
    b2 = b2 - innerprod_q2(basis1,b2)*basis1;
    basis2 = b2 / sqrt(innerprod_q2(b2, b2));

    basis_o = Array(Any,2);
    basis_o[1] = basis1;
    basis_o[2] = basis2;

    return basis_o
end

"""
Project srvf onto tangent space

    project_tangent(w, q, basis)
    :param w: shooting vector
    :param q: array (n,T)
    :param basis: basis out of find_basis_normal
"""
function project_tangent(w, q, basis)
    w = w - innerprod_q2(w, q) * q;
    bo = gram_schmidt(basis);

    wproj = w - innerprod_q2(w, bo[1])*bo[1] - innerprod_q2(w,bo[2])*bo[2];

    return wproj
end

"""
Scale curve to unit length

    scale_curve(beta)
    :param beta: array (n,T)
"""
function scale_curve(beta)
    n, T = size(beta);
    normbetadot = zeros(T);
    betadot = zeros(n,T);
    for i = 1:n
        betadot[i,:] = gradient(vec(beta[i,:]), 1.0/T);
    end
    for i = 1:T
        normbetadot[i] = norm(betadot[:,i]);
    end
    scale = trapz(collect(linspace(0,1,T)), normbetadot);
    beta_scaled = beta / scale;

    return beta_scaled, scale
end

"""
Parallel translate srvf  along tangent space

    parallel_translate(w, q1, q2, basis, mode='O')
    :param w: shooting vector
    :param q1: array (n,T)
    :param q2: array (n,T)
    :param basis: basis out of find_basis_normal
    :param mode: Open ('O') or Closed ('C') curves
"""
function parallel_translate(w, q1, q2, basis, mode='O')
    wtilde = w - 2*innerprod_q2(w,q2) / innerprod_q2(q1+q2, q1+q2) * (q1+q2);
    l = sqrt(innerprod_q2(wtilde, wtilde));

    if mode == 'C'
        wbar = project_tangent(wtilde, q2, basis);
        normwbar = sqrt(innerprod_q2(wbar, wbar));
        if normwbar > 1e-4
            wbar = wbar * l / normwbar;
        end
    else
        wbar = wtilde;
    end

    return wbar
end

"""
Calculate zero crossing for optimal regression between curves
"""
function curve_zero_crossing(Y, beta, bt, y_max, y_min, gmax, gmin)
    T = size(beta,2);
    betanu = q_to_curve(bt);
    max_itr = 100;
    a = zeros(max_itr);
    a[1] = 1.;
    f = zeros(max_itr);
    f[1] = y_max - Y;
    f[2] = y_min - Y;
    mrp = f[1];
    mrn = f[2];
    mrp_ind = 1;  # most recent positive index
    mrn_ind = 2;  # most recent negative index

    for ii = 3:max_itr
        x1 = a[mrp_ind];
        x2 = a[mrn_ind];
        y1 = mrp;
        y2 = mrn;
        a[ii] = (x1 * y1 - x2 * y2) / float(y2-y1);

        beta1, O_hat, tau = find_rotation_seed_coord(betanu, beta);
        gamma = a[ii] * gmax + (1-a[ii]) * gmin;

        beta1 = group_action_by_gamma_coord(beta1, gamma);
        beta1, O_hat1, tau = find_rotation_seed_coord(betanu, beta1);
        q1 = curve_to_q(beta1);
        f[ii] = innerprod_q2(q1, bt) - Y;

        if abs(f[ii]) < 1e-5
            break
        elseif f[ii] > 0
            mrp = f[ii];
            mrp_ind = ii;
        else
            mrn = f[ii];
            mrn_ind = ii;
        end
    end

   gamma = a[ii] * gmax + (1-a[ii]) * gmin;

   beta1, O_hat, tau = find_rotation_seed_coord(betanu, beta1);
   beta1 = group_action_by_gamma_coord(beta1, gamma);
   beta1, O_hat1, tau = find_rotation_seed_coord(betanu, beta1);
   O_hat = O_hat*O_hat1;

  return gamma, O_hat, tau
end

"""
Find rotation matrix for angle theta

    rot_mat(theta)
"""
function rot_mat(theta)
    O = [cos(theta) -1*sin(theta); sin(theta) cos(theta)];

    return O
end

"""
Find inner product between two srvfs, q1 and q2

    innerprod_q2(q1, q2)
    :param q1: array (n,T)
    :param q2: array (n,T)
"""
function innerprod_q2(q1, q2)
    T = size(q1,2);
    val = sum(sum(q1.*q2))/T;

    return val
end

"""
Calculate elastic shape distance between two curves beta1 and beta2

    calc_shape_dist(beta1, beta2)
    :param beta1: array (n,T)
    :param beta2: array (n,T)
    :param method: optimization method to find warping, default is
                   Dynamic Programming ("DP"). Other options are
                   Coordinate Descent ("DP2"), Riemannian BFGS
                   ("RBFGS")
    :param wscale: with scale (false)
"""
function calc_shape_dist(beta1::Array{Float64,2}, beta2::Array{Float64,2};
                         method::AbstractString="DP", wscale=false)

    T1 = size(beta1,2)
    T2 = size(beta2,2)

    if (T1 != T2)
        if (T1>T2)
            beta2 = resamplecurve(beta2, T1)
        else
            beta1 = resamplecurve(beta1, T2)
        end
    end

    v, d = inverse_exp_coord(beta1, beta2, wscale=wscale)

    return d
end

"""
Pairwise align two curves

    curve_pair_align(beta1::Array{Float64,2}, beta2::Array{Float64,2})
    :param beta1: array (n,T)
    :param beta2: array (n,T)

    Returns
    :return beta2n: aligned curve 2 to 1
    :return q2n: aligned srvf 2 to 1
    :return gam: warping function
    :return q1: srvf of curve 1
"""
function curve_pair_align(beta1::Array{Float64,2}, beta2::Array{Float64,2})

    T = size(beta1,2);
    centroid1 = calculatecentroid(beta1);
    beta1 -= repmat(centroid1,1,T);
    centroid2 = calculatecentroid(beta2);
    beta2 -= repmat(centroid2,1,T);

    q1 = curve_to_q(beta1);

    # Iteratively optimize over SO(n) x Gamma using old DP
    gam, R, tau = optimum_reparam(beta1, beta2);
    beta2n = R * shift_f(beta2, tau);
    gamI = invert_gamma(gam);
    beta2n = group_action_by_gamma_coord(beta2n, gamI);
    beta2n, R, tau = find_rotation_seed_coord(beta1, beta2n);
    q2n = curve_to_q(beta2n);

    return beta2n, q2n, gamI, q1
end

"""
Form geodesic between two curves

    curve_geodesic(beta1::Array{Float64,2}, beta2::Array{Float64,2}, k::Integer=5)
    :param beta1: array (n,T)
    :param beta2: array (n,T)
    :param k: number of curves along geodesic
    :param wscale: with scale (false)

    Returns
    :return geod: curves along geodesic (n,T,k)
    :return geod_q: srvf's along geodesic
"""
function curve_geodesic(beta1::Array{Float64,2}, beta2::Array{Float64,2},
                        k::Integer=5; wscale=false)

    n, T = size(beta1);
    beta1 = resamplecurve(beta1, T);
    beta2 = resamplecurve(beta2, T);
    centroid1 = calculatecentroid(beta1);
    beta1 -= repmat(centroid1,1,T);
    centroid2 = calculatecentroid(beta2);
    beta2 -= repmat(centroid2,1,T);

    q1 = curve_to_q(beta1);

    # Iteratively optimize over SO(n) x Gamma using old DP
    gam, R, tau = optimum_reparam(beta1, beta2);
    beta2n = R * shift_f(beta2, tau);
    gamI = invert_gamma(gam);
    beta2n = group_action_by_gamma_coord(beta2n, gamI);
    beta2n, R, tau = find_rotation_seed_coord(beta1, beta2n);
    q2n = curve_to_q(beta2n);

    # form geodesic between the registered curves
    if wscale
        dist = sqrt(innerprod_q2(q1-q2n,q1-q2n))
    else
        dist = acos(innerprod_q2(q1,q2n));
    end
    geod = zeros(n,T,k);
    geod_q = zeros(n,T,k);

    for tau = 1:k
        s = dist*(tau-1)/(k-1);
        geod_q[:,:,tau] = (sin(dist-s)*q1+sin(s)*q2n)/sin(dist);
        geod[:,:,tau] = q_to_curve(geod_q[:,:,tau]);
    end

    return geod, geod_q
end


text:
#=
## Network Dynamics without NetworkDynamics.jl
=#
#src In this example we want to explore the same problem as in [Integration with `NetworkDynamics.jl`](@ref).
#src But this time without `NetworkDynamics.jl` ...
#=
In this example we model a Kuramoto system on a complex network.
=#
using LightGraphs
using BlockSystems
using ModelingToolkit
using OrdinaryDiffEq
using Plots

#=
Dynamical systems on complex networks are best modelled on directed graphs, because in
general an interaction of node i and j may be asymmetric, e.g. when node i influences
node j but not vice versa.

Symmetries of the system such as (anti-)symmetric coupling functions may justify other
approaches, but those are of limited scope, for a detailed discussion see
[our paper](https://arxiv.org/abs/2012.12696), Section II.D.

In the following each edge is represented by a function:

```
     e₁₂ = f(1,2)
     .--->---.
   (1)       (2)
     ˙---<---˙
     e₂₁ = f(2,1)
```

Each edge function sees the values of the connected nodes, its source and its destination.
Nodes are modeled as functions as well and see all their *incoming* edges.

The goal is to generate IOBlocks for edges and vertices based on a given graph.

While defining the edge block is straightforward, we have to make the vertex blocks a bit
special: since MTK does not support vector inputs yet we need a special IOBlock which
depends on the number of connected edges.
=#
function gen_edge_block(name)
    @parameters t src(t) dst(t) K
    @variables out(t)
    IOBlock([out ~ K*sin(src-dst)], [src, dst], [out], name=Symbol(name))
end

## append subscript, `sybscript(:foo, 2) => :foo₂`
subscript(s, i) = Symbol(s, Char(0x02080 + i))

function gen_vertex_block(n_edges, name)
    @parameters t ω
    @variables ϕ(t)
    D = Differential(t)
    ## the way array variables work changed. This is a hack to retrieve the old behvaviour of
    ## this closely mimics the old @parameters edge[1:n_edges](t)
    edge = Num[]
    for i in 1:n_edges
        symname = subscript(:edge, i)
        append!(edge, @parameters $symname(t))
    end

    IOBlock([D(ϕ) ~ ω + (+)(edge...)],
            [edge...],
            [ϕ],
            name=Symbol(name))
end
nothing #hide

#=
For simplicity our graph will be a simple ring network specified with LightGraphs.jl

=#
N = 8
g = SimpleDiGraph(watts_strogatz(N,2,0)) # ring network
nothing #hide

# First we generate a list of all edge-blocks because they don't depend on the vertices.

edgelist = [(i=i, src=e.src, dst=e.dst, block=gen_edge_block("e_$(e.src)_$(e.dst)"))
            for (i, e) in enumerate(edges(g))]
edge_blocks = [e.block for e in edgelist]
nothing #hide

#=
Now we can generate vertex blocks based on their number of incoming edges.
We will also create the connections
```
e_i_k.out => node.edge₁
e_j_k.out => node.edge₂
...
```
for all edges that go from vertices i or j to vertex k.
=#
vert_blocks = IOBlock[]
connections = Pair[]

for i in vertices(g)
    ## collect the incoming edges for each node
    edges = filter(e -> e.dst == i, edgelist)

    node = gen_vertex_block(length(edges), "v$i")
    push!(vert_blocks, node)

    ## each node has the open inputs edge₁, edge₂, ...
    ## we need to connect the outputs of the edge-blocks to the
    ## inputs of the node like edge_j_to_1.out => node.edge₁
    for (i, edge) in enumerate(edges)
        node_input_i = getproperty(node, subscript(:edge, i))
        push!(connections, edge.block.out => node_input_i)
    end
end

# Once the vertices are generated we can plug the edges' src and dst to the output of the #
# corresponding vertex, in this case its the oscillators angle ϕ
for edge in edgelist
    push!(connections, vert_blocks[edge.src].ϕ => edge.block.src)
    push!(connections, vert_blocks[edge.dst].ϕ => edge.block.dst)
end

#=
We want the `connect_system` to get rid of the algebraic states for the edges. Therefore
we have to provide a list of `outputs` which only contains the outputs of the vertices. By
doing so the edge outputs will become internal `istates` of the `IOSystem` and upon
connection may be reduced.
=#
outputs = [block.ϕ for block in vert_blocks]

network = IOSystem(connections, vcat(vert_blocks, edge_blocks), outputs=outputs)

networkblock = connect_system(network, verbose=false)
nothing #hide

# As the output shows the system has be reduced to just N equations.
# Well now we can generate the functions...
gen = generate_io_function(networkblock,
                           f_states=[v.ϕ for v in vert_blocks],
                           f_params=vcat([v.ω for v in vert_blocks],
                                         [e.K for e in edge_blocks]),
                           warn=false);
nothing #hide

# ... enclose the `f_ip` to get rid of the empty `inputs` field...
odefun(du, u, p, t) = gen.f_ip(du, u, (), p, t)
nothing #hide

# ... set the starting conditions ...
ω = collect(1:N)./N
ω .-= sum(ω)/N
K = [3.0 for i in edge_blocks]
p = (ω..., K...)

x0 = collect(1:N)./N
x0 .-= sum(x0)./N
nothing #hide

#src using BenchmarkTools
#src @btime $odefun($x0, $x0, $p, 0.0)

# ... and solve the system!
tspan = (0., 10.)
prob = ODEProblem(odefun, x0, tspan, p)
sol = solve(prob, Tsit5())
plot(sol, ylabel="ϕ")

ext: jl


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Behavior: suggestmissingsteps, ParseOptions

if length(ARGS) !== 2 && length(ARGS) !== 3
    println("Usage: julia suggeststeps.jl <feature file> <steps root path> [--experimental]")
    exit(1)
end

featurefile = ARGS[1]
stepsrootpath = ARGS[2]

const use_experimental = length(ARGS) == 3 && ARGS[3] == "--experimental"

parseoptions = ParseOptions(allow_any_step_order=true, use_experimental=use_experimental)

suggestmissingsteps(featurefile, stepsrootpath, parseoptions=parseoptions)


text:
function locS(w::Real, b::Real, c::Real, d::Real)
    S = (w^2 + c^2)^(-d) + b^2
  return S
end

function locS(b::Real, c::Real, d::Real)
  S = NaN*zeros(10000); n_w = length(S)
  w_grid = collect(linspace(0.001, 15*π, n_w))
  i = 0
  for w in w_grid
    i += 1
    S[i] = (w^2 + c^2)^(-d) + b^2
  end
  return S
end

function locSs(b::Real, c::Real, d::Real)
  S = locS(b, c, d)
  n = length(S)
  w_grid = collect(linspace(0.001, 15π, n))
  tmp, i1 = findmin(abs(w_grid-6*π))
  Ss = S./S[i1]
  x = w_grid/π
  return Ss, x
end
#
function plot_locS(b::Real, c::Real, d::Real)
    title("Local to Zero Spectrum: S(ω;b,c,d)∝(c² + d²)^(-1) + b²")
    Ss, x = locSs(b, c, d)
    _plt = plot(x, log(Ss), label = "($(round(b, 2)), $(round(c, 2)), $(round(d, 2)))")
    legend()
    return _plt
end
#
function plot_locS(m::String)
  if m == "local level"
    println("$m model")
    b_grid = [1/50 1/10 1/2 10e16]
    for b in b_grid
      Ss, x = locSs(b, 0, 1)
      plot(x, log(Ss), label = "b = $b")
      legend()
    end

  end
  if m == "local to unity"
    println("$m model")
    c_grid = [0 3 10 30]
    for c in c_grid
      Ss, x = locSs(0, c, 1)
      plot(x, log(Ss), label = "c = $c")
    end
    legend()
  end
  if m == "fractional"
    println("$m model")
    d_grid = [-1/3 0 1/3 1 4/3]
    for d in d_grid
      Ss, x = locSs(0, 0, d)
      plot(x, log(Ss), label = "d = $d")
    end
    legend()
  end
  if m == "bcd"
    b_grid = [1/50 1/10 1/2 10e16]
    c_grid = [0 3 10 30]
    d_grid = [-1/3 0 1/3 1 4/3]
    for b in b_grid
      for c in c_grid
        for d in d_grid
          println("$b, $c, $d")
          Ss, x = locSs(b, c, d)
          plot(x, log(Ss))
        end
      end
    end
  end
end

#end #module


text:
export ◺, commutation, CopyMatrix, kron_axpy!, kr_axpy!, kr, mul!, vech,
    Ct_At_kron_A_KC, Ct_At_kron_A_KC!,
    Ct_A_kron_B_C, Ct_A_kron_B_C!,
    Ct_A_kr_B, Ct_A_kr_B!
import LinearAlgebra: mul!

"""
    ◺(n::Integer)

Triangular number `n * (n + 1) / 2`.
"""
@inline ◺(n::Integer) = (n * (n + 1)) >> 1

"""
    kron_axpy!(A, X, Y)

Overwrite `Y` with `A ⊗ X + Y`. Same as `Y += kron(A, X)` but
more memory efficient.
"""
function kron_axpy!(
    A::AbstractVecOrMat{T},
    X::AbstractVecOrMat{T},
    Y::AbstractVecOrMat{T}
) where T <: Real
    m, n = size(A, 1), size(A, 2)
    p, q = size(X, 1), size(X, 2)
    @assert size(Y, 1) == m * p
    @assert size(Y, 2) == n * q
    @inbounds for j in 1:n
        coffset = (j - 1) * q
        for i in 1:m
            a = A[i, j]
            roffset = (i - 1) * p
            for l in 1:q
                r = roffset + 1
                c = coffset + l
                for k in 1:p
                    Y[r, c] += a * X[k, l]
                    r += 1
                end
            end
        end
    end
    Y
end

"""
    kr_axpy!(A, X, Y)

Overwrite `Y` with `A ⊙ X + Y`, where `⊙` stands for the Khatri-Rao (columnwise
Kronecker) product. `A` and `X` need to have same number of columns.
"""
function kr_axpy!(
    A::AbstractVecOrMat{T},
    X::AbstractVecOrMat{T},
    Y::AbstractVecOrMat{T}
) where T <: Real
    @assert size(A, 2) == size(X, 2) == size(Y, 2)
    m, n, p = size(A, 1), size(A, 2), size(X, 1)
    @inbounds for j in 1:n
        r = 1
        for i in 1:m
            aij = A[i, j]
            for k in 1:p
                Y[r, j] += aij * X[k, j]
                r += 1
            end
        end
    end
    Y
end

kr(A::AbstractVecOrMat{T}, X::AbstractVecOrMat{T}) where T <: Real =
    kr_axpy!(A, X, zeros(T, size(A, 1) * size(X, 1), size(A, 2)))

struct CopyMatrix <: AbstractMatrix{Int}
    n::Int
end

Base.size(C::CopyMatrix) = (abs2(C.n), (C.n * (C.n + 1)) >> 1)

Base.IndexStyle(::Type{<:CopyMatrix}) = IndexCartesian()

function Base.getindex(C::CopyMatrix, i::Int, j::Int)
    r, c = CartesianIndices((1:C.n, 1:C.n))[i].I
    if r ≥ c && j == (c - 1) * C.n - ((c - 2) * (c - 1)) >> 1 + r - c + 1
        return 1
    else
        return 0
    end
end

"""
    mul!(result, A, C::CopyMatrix)

Right-multiplying a matrix `A` by a copying matrix is equivalent to keeping
the columns of `A` corresponding to the lower triangular indices.
"""
function LinearAlgebra.mul!(
    result :: AbstractVecOrMat,
    A      :: AbstractVecOrMat,
    C      :: CopyMatrix
)
    n = isqrt(size(A, 2))
    m = size(A, 1)
    @assert size(result, 1) == m
    @assert size(result, 2) == (n * (n + 1)) >> 1
    ac, rc = 0, 0
    @inbounds for j in 1:n, i in 1:n
        ac += 1
        i < j && continue
        rc += 1
        for k in 1:m
            result[k, rc] = A[k, ac]
        end
    end
    result
end

"""
    mul!(result, Ct::Transpose{Int, CopyMatrix}, A)

Left-multiplying a matrix `A` by transpose of a copying matrix is equivalent to
keeping the rows of `A` corresponding to the lower triangular indices.
"""
LinearAlgebra.mul!(
    result :: AbstractVecOrMat,
    Ct     :: Transpose{Int, CopyMatrix},
    A      :: AbstractVecOrMat
) = mul!(transpose(result), transpose(A), Ct.parent)

"""
    vech!(v::AbstractVector, A::AbstractVecOrMat)

Overwrite vector `v` by the entries from lower triangular part of `A`.
"""
function vech!(v::AbstractVector, A::AbstractVecOrMat)
    m, n = size(A, 1), size(A, 2)
    idx = 1
    @inbounds for j in 1:n, i in j:m
        v[idx] = A[i, j]
        idx += 1
    end
    v
end

function commutation(m::Integer, n::Integer)
    K = zeros(Int, m * n, m * n)
    colK = 1
    @inbounds for j in 1:n, i in 1:m
        rowK          = n * (i - 1) + j
        K[rowK, colK] = 1
        colK         += 1
    end
    K
end

"""
    vech(A::AbstractVecOrMat) -> AbstractVector

Return the entries from lower triangular part of `A` as a vector.
"""
function vech(A::AbstractVecOrMat)
    m, n = size(A, 1), size(A, 2)
    vech!(similar(A, n * m - (n * (n - 1)) >> 1), A)
end

"""
    Ct_At_kron_A_KC!(H, A, B)

Overwrite `H` by `H + C'(A'⊗A)KC`, where `K` is the commutation matrix and
`C` is the copying matrix.
"""
function Ct_At_kron_A_KC!(H::AbstractMatrix, A::AbstractMatrix)
    q = size(A, 1)
    @assert size(A, 2) == q
    @assert size(H, 1) == size(H, 2) == (q * (q + 1)) >> 1
    j = 1
    @inbounds for w in 1:q, s in w:q
        i = 1
        for r in 1:q, v in r:q
            H[i, j] += A[s, r] * A[v, w]
            i += 1
        end
        j += 1
    end
    H
end

function Ct_At_kron_A_KC(A)
    n◺ = ◺(size(A, 1))
    H = zeros(eltype(A), n◺, n◺)
    Ct_At_kron_A_KC!(H, A)
end

"""
    Ct_A_kron_B_C!(H, A, B)

Overwrite `H` by `H + C'(A⊗B)C`, where `C` is the copying matrix.
"""
function Ct_A_kron_B_C!(
    H::AbstractMatrix,
    A::AbstractMatrix,
    B::AbstractMatrix,
    )
    q = size(A, 1)
    @assert size(A, 2) == size(B, 1) == size(B, 2) == q
    j = 1
    @inbounds for s in 1:q, w in s:q
        i = 1
        for r in 1:q, v in r:q
            H[i, j] += A[r, s] * B[v, w]
            i += 1
        end
        j += 1
    end
    H
end

function Ct_A_kron_B_C(A, B)
    n◺ = ◺(size(A, 1))
    H = zeros(eltype(A), n◺, n◺)
    Ct_A_kron_B_C!(H, A, B)
end

"""
    Ct_A_kr_B!(H, A, B)

Overwrite `H` by `H + C'(A⊙B)`, where `C` is the copying matrix and `⊙` is the
Khatri-Rao (column-wise Kronecker) product.
"""
function Ct_A_kr_B!(H, A, B)
    @assert size(A) == size(B)
    (q, n) = size(A)
    @inbounds for c in 1:n
        r = 1
        for ia in 1:q
            a = A[ia, c]
            for ib in ia:q
                H[r, c] += a * B[ib, c]
                r += 1
            end
        end
    end
    H
end

function Ct_A_kr_B(A, B)
    @assert size(A) == size(B)
    (q, n) = size(A)
    H = zeros(eltype(A), ◺(q), n)
    Ct_A_kr_B!(H, A, B)
end


text:
using Logging
using Reexport

import JLD2
@reexport using ArgParse

const HASH_KEY="_HASH"
const SAVE_NAME_KEY="_SAVE"
const GIT_INFO_KEY="_GIT_INFO"

make_save_name(hashed, git_info; head="RP") = "$(head)_$(git_info)_0x$(string(hashed,base=16))"

get_save_dir(parsed::Dict) = parsed[keytype(parsed)(SAVE_NAME_KEY)]
get_hash(parsed::Dict) = parsed[keytype(parsed)(HASH_KEY)]
get_git_info(parsed::Dict) = parsed[keytype(parsed)(GIT_INFO_KEY)]

"""
    create_info!
"""
function create_info!(parsed_args::Dict,
                      save_dir::String;
                      filter_keys::Array{String,1} = Array{String,1}(),
                      use_git_info = false,
                      custom_folder_name = "RP",
                      HASHER=hash,
                      replace=true,
                      settings_file="settings.jld2")

    KEY_TYPE = keytype(parsed_args)

    unused_keys = KEY_TYPE.(filter_keys)
    hash_args = filter(k->(!(k[1] in unused_keys)), parsed_args)
    used_keys=keys(hash_args)

    hash_key = KEY_TYPE(HASH_KEY)
    save_name_key = KEY_TYPE(SAVE_NAME_KEY)
    git_info_key = KEY_TYPE(GIT_INFO_KEY)

    hashed = HASHER(hash_args)
    parsed_args[hash_key] = hashed

    git_info = use_git_info ? git_head() : "0"
    parsed_args[git_info_key] = git_info

    save_name = joinpath(save_dir, make_save_name(hashed, git_info; head=custom_folder_name))
    parsed_args[save_name_key] = save_name

    save_settings_path = save_name

    if !isdir(save_settings_path)
        mkpath(save_settings_path)
    else
        if replace
            @warn "Hash Conflict in Reproduce create_info! Overwriting data."
        else
            @info "Told not to replace. Exiting Experiment."
            throw("Hash Conflict.")
        end
    end

    # settings_dict = Dict("parsed_args"=>parsed_args, "used_keys"=>used_keys)
    save_settings_file = joinpath(save_settings_path, settings_file)
    JLD2.@save save_settings_file parsed_args used_keys
end

function create_info(arg_list::Vector{String},
                     settings::ArgParseSettings,
                     save_dir::AbstractString;
                     filter_keys::Array{String,1} = Array{String,1}(),
                     use_git_info = false,
                     custom_folder_name = "RP",
                     HASHER=hash,
                     replace=true,
                     settings_file="settings.jld2", kwargs...)
    parsed = parse_args(arg_list, settings; kwargs...)
    create_info!(parsed, save_dir;
                 filter_keys=filter_keys,
                 use_git_info=use_git_info,
                 HASHER=HASHER,
                 custom_folder_name=custom_folder_name,
                 replace=replace,
                 settings_file=settings_file)
    return parsed
end

function default_save_str(parsed, use_keys; save_dir="RP")
    strs = collect(zip(
        [string(key) for key in use_keys],
        [string(parsed[key]) for key in use_keys]))
    dir = joinpath(homedir, strs...)
    return dir
end

function create_custom_info!(parsed_args::Dict,
                             save_dir::String;
                             use_keys::Array{String,1} = Array{String,1}(),
                             make_save_str=default_save_str,
                             replace=true,
                             settings_file="settings.jld2")

    KEY_TYPE = keytype(parsed_args)

    save_name_key = KEY_TYPE(SAVE_NAME_KEY)
    save_path = make_save_str(parsed_args, KEY_TYPE.(use_keys); save_dir=save_dir)
    parsed_args[save_name_key] = save_path

    if !isdir(save_path)
        mkpath(save_path)
    else
        if replace
            @warn "Dir Conflict in Reproduce parse_args! Overwriting data."
        else
            @info "Told not to replace. Exiting Experiment."
            exit(0)
        end
    end

    # settings_dict = Dict("parsed_args"=>parsed_args, "used_keys"=>used_keys)
    save_settings_file = joinpath(save_path, settings_file)
    JLD2.@save save_settings_file parsed_args used_keys
end

function create_custom_info(arg_list::Vector{String},
                            settings::ArgParseSettings,
                            save_dir::AbstractString;
                            use_keys::Array{String,1} = Array{String,1}(),
                            make_save_str=default_save_str,
                            replace=true,
                            settings_file="settings.jld2", kwargs...)
    parsed = parse_args(arg_list, settings; kwargs...)
    create_custom_info!(parsed, save_dir;
                        use_keys=use_keys,
                        make_save_str=make_save_str,
                        replace=replace,
                        settings_file=settings_file)
    return parsed
end


text:
################################################################################
# MathProgBase solver interface - GMMEstimator
################################################################################

function MathProgBase.initialize(d::GMMEstimator, rf::Vector{Symbol})
    for feat in rf
        if !(feat in [:Grad, :Jac, :Hess])
            error("Unsupported feature $feat")
        end
    end
end

MathProgBase.features_available(e::GMMEstimator) = [:Grad, :Jac, :Hess]

function MathProgBase.eval_f(e::GMMEstimator, theta)
    idx = e.ist.n[1]
    gn = vec(sum(e.mf.s(theta), dims=1))
    LinearAlgebra.dot(gn, e.W[idx]*gn)
end

MathProgBase.eval_g(e::GMMEstimator{M, V, T, S}, g, theta) where {M, V, T<:Unconstrained, S} = nothing

function MathProgBase.eval_g(e::GMMEstimator{M, V, T, S}, g, theta) where {M, V, T<:Constrained, S}
    g[:] .= e.c.h(theta)
end

function MathProgBase.eval_grad_f(e::GMMEstimator{M, V, T, S}, grad_f, θ) where {M<:FADMomFun, V, T, S}
    idx = e.ist.n[1]
    sn(θ) = vec(sum(e.mf.s(θ), dims=1))
    Δsn = ForwardDiff.jacobian(sn, θ)::Matrix{Float64}  # m x k
    Wg  = e.W[idx]*sn(θ)
    mul!(grad_f, Δsn', Wg)
end

function MathProgBase.eval_grad_f(e::GMMEstimator{M, V, T, S}, grad_f, θ) where {M <: AnaMomFun, V, T, S}
    ##grad_f[:] = 2*d.Dsn(theta)'*(d.W*d.sn(theta))
    idx = e.ist.n[1]
    sn = vec(sum(e.mf.s(θ), 1))
    gemm!('T', 'N', 2.0,
          e.mf.Dsn(θ),
          e.W[idx]*sn, 0.0, grad_f)
end

MathProgBase.jac_structure(e::GMMEstimator{M, V, T, S}) where {M, V, T<:Unconstrained, S} = Int[],Int[]
MathProgBase.eval_jac_g(e::GMMEstimator{M, V, T, S}, J, x) where {M, V, T<:Unconstrained, S} = nothing

function MathProgBase.jac_structure(e::GMMEstimator{M, V, T, S}) where {M, V, T<:Constrained, S}
    nc = e.c.nc            ## Number of constraints
    n, k, m = size(e)
    ## The jacobian is a nc x k
    rows = Array(Int64, e.gele)
    cols = Array(Int64, e.gele)
    for j = 1:nc, r = 1:k
        @inbounds rows[r+(j-1)*k] = j
        @inbounds cols[r+(j-1)*k] = r
    end
    rows, cols
end

function MathProgBase.eval_jac_g(e::GMMEstimator{M, V, T, S}, J, θ) where {M, V, T<:Constrained, S}
    h(θ) = e.c.h(θ)
    J[:] = vec((ForwardDiff.jacobian(h, θ)::Matrix{Float64})')
end

MathProgBase.hesslag_structure(d::GMMEstimator) = Int[], Int[]
MathProgBase.eval_hesslag(d::GMMEstimator, H, x, σ, μ) = nothing


text:
module Tsys

export model_flux
export apparent_flux
export apparent_temperature
export tcal_onoff
export tsys_onoff

"""
Boltzmann constant (Joules/Kelvin)
"""
Kb = 1.380469e-23

"""
Janskys per Joule
"""
JY_PER_J = 1e26

"""
Temperature of the *cosmic microwave background* (CMB) in Kelvin
"""
Tcmb = 2.728

"""
`Tsys.ν1` is a `Ref{Float64}` that holds the value used as the default
calibrator flux model frequency, in Hz, for `ν1` in `flux`.  The initial value
is `1e6` (i.e. 1 MHz), but the user may change this to match the frequency
convention of the flux models they will be using.  For example, the flux models
presented in the 2017 paper by Perley and Butler use frequency in GHz, so the
user could set `Tsys.ν1[] = 1e9` to make it more convenient to work with those
flux models.  Of course, one can always pass an explicit value for `ν1` when
calling `model_flux`, which allows for comparing different flux models with
different frequency conventions.
"""
ν1 = Ref{Float64}(1e6)

"""
    model_flux(model, hz=ν1[], ν1=ν1[])

Compute the spectral flux density at frequency `hz` (in Hz) given flux model
coefficients `model` based on frequency units of `ν1` Hz.  `hz` and `ν1` default
to `Tsys.ν1[]` if not given.
"""
function model_flux(model, hz=ν1[], ν1=ν1[])
    10^evalpoly(log10(hz/ν1), model)
end

"""
    apparent_flux(k, diameter; η=1.0, τ=0.0, airmass=1.0)

Return the apparent spectral flux density in Janskys corresponding to black body
temperature `k` in Janskys and antenna diameter `diameter` in meters with
aperture efficiency of `η`, which defaults to 1.0.  Atmospheric opacity `τ`
defaults to 0.0 (totally transparent) and `airmass` defaults to 1.0 (airmass at
zenith).
"""
function apparent_flux(k, diameter; η=1.0, τ=0.0, airmass=1.0)
	a = η * π * diameter^2 / 4
    (2Kb * JY_PER_J * k) / (a * exp(-τ*airmass))
end

"""
    apparent_temperature(s, diameter; η=1.0, τ=0.0, airmass=1.0)

Return the apparent black body temperature in Kelvin corresponding to spectral
flux density `s` in Janskys and antenna diameter `diameter` in meters with
aperture efficiency of `η`, which defaults to 1.0.  Atmospheric opacity `τ`
defaults to 0.0 (totally transparent) and `airmass` defaults to 1.0 (airmass at
zenith).
"""
function apparent_temperature(s, diameter; η=1.0, τ=0.0, airmass=1.0)
	a = η * π * diameter^2 / 4
	s * a * exp(-τ*airmass) / (2Kb * JY_PER_J)
end

"""
    tcal_onoff(pon, poff, tsys; tsky=Tcmb, clip=true)

Calculates calibrator apparent temperature `Tcal`.  The `pon` and `poff`
parameters are the power measurements when pointing ON and OFF the calibrator
using a radio telescope with system temperature `tsys` (in Kelvin).  `tsky` is
the apparent temperature in Kelvin of the "blank sky" at the OFF position.  If
`clip` is true, any computed value less than 0.0 will be clamped to 0.0.

The calculation is derived by observing that:

    pon  = G(Tsky + Tsys + Tcal)
    poff = G(Tsky + Tsys)

Then:

     pon - poff        Tcal
    ------------ = -------------
         poff       Tsky + Tsys

Leads to:

            pon - poff
    Tcal = ------------ * (Tsky + Tsys)
               poff
"""
function tcal_onoff(pon, poff, tsys; tsky=Tcmb, clip=true)
    lim = clip ? 0.0 : -Inf
    y = (pon - poff) / poff
    clamp(y * (tsky + tsys), lim, Inf)
end

"""
    tsys_onoff(pon, poff, tcal; tsky=Tcmb, clip=true)

Calculates system temperature `Tsys`.  The `pon` and `poff` parameters are the
power measurements when pointing ON and OFF the calibrator with apparent
temperature `tcal` (in Kelvin).  `tsky` is the apparent temperature in Kelvin of
the "blank sky" at the OFF position.  If `clip` is true, any computed value less
than 0.0 will be clamped to 0.0.

The calculation is derived by observing that:

    pon  = G(Tsky + Tsys + Tcal)
    poff = G(Tsky + Tsys)

Then:

     pon - poff        Tcal
    ------------ = -------------
         poff       Tsky + Tsys

Leads to:

                     poff
    Tsys = Tcal * ------------ - Tsky
                   pon - poff
"""
function tsys_onoff(pon, poff, tcal; tsky=Tcmb, clip=true)
    lim = clip ? 0.0 : -Inf
    y = (pon - poff) / poff
    clamp(tcal / y - tsky, lim, Inf)
end

"""
    tsys_onoff(pon, poff, scal, diameter; kwargs...)

Returns Tsys by computing Tcal from `scal` and passing it to `tsys_onoff`.
Keyword arguments `kwargs` are the same as for `apparent_temperature` and the
three positional argument method of `tsys_onoff`.
"""
function tsys_onoff(pon, poff, scal, diameter;
              η=1.0, τ=0.0, airmass=1.0, tsky=Tcmb, clip=true)
    tcal = apparent_temperature(scal, diameter; η=η, τ=τ, airmass=airmass)
    tsys_onoff(pon, poff, tcal; tsky=tsky, clip=clip)
end

end # module Tsys


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@with_kw mutable struct GradientThresholdParameters{T <: AbstractFloat}
    τ::T = 1e-6
end

"""
    ObjectiveThreshold

Functor object for using a zero gradient threshold as termination criterion in a quasigradient algorithm. Create by supplying a [`AtGradientThreshold`](@ref) object through `terminate` to `QuasiGradient.Optimizer` or by setting the [`Termination`](@ref) attribute.

...
# Parameters
- `τ::AbstractFloat = 1e-6`: Numerical tolerance for zero gradient
...
"""
struct GradientThreshold{T <: AbstractFloat} <: AbstractTerminationCriterion
    parameters::GradientThresholdParameters{T}

    function GradientThreshold(::Type{T}; kw...) where T <: AbstractFloat
        return new{T}(GradientThresholdParameters{T}(; kw...))
    end
end

function Progress(termination::GradientThreshold, str::AbstractString)
    @unpack τ = termination.parameters
    return ProgressThresh(τ, 0.0, str)
end

function progress_value(::GradientThreshold, k::Integer, f::AbstractFloat, ∇f_norm::AbstractFloat)
    return ∇f_norm
end

function terminate(termination::GradientThreshold, k::Integer, f::Float64, x::AbstractVector, ∇f::AbstractVector)
    return norm(∇f) <= termination.τ
end

# API
# ------------------------------------------------------------
"""
    AtGradientThreshold

Factory object for [`GradientThreshold`](@ref). Pass to `terminate` in `Quasigradient.Optimizer` or set the [`Termination`](@ref) attribute. See ?GradientThreshold for parameter descriptions.

"""
mutable struct AtGradientThreshold <: AbstractTermination
    parameters::GradientThresholdParameters{Float64}
end
AtGradientThreshold(τ::AbstractFloat) = AtGradientThreshold(GradientThresholdParameters(; τ = Float64(τ)))
AtGradientThreshold(; kw...) = AtGradientThreshold(GradientThresholdParameters(; kw...))

function (criteria::AtGradientThreshold)(::Type{T}) where T <: AbstractFloat
    return GradientThreshold(T; type2dict(criteria.parameters)...)
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: compute_optimizer
using AWS.Compat
using AWS.UUIDs

"""
    DescribeRecommendationExportJobs()

Describes recommendation export jobs created in the last seven days. Use the ExportAutoScalingGroupRecommendations or ExportEC2InstanceRecommendations actions to request an export of your recommendations. Then use the DescribeRecommendationExportJobs action to view your export jobs.

# Optional Parameters
- `filters`: An array of objects that describe a filter to return a more specific list of export jobs.
- `jobIds`: The identification numbers of the export jobs to return. An export job ID is returned when you create an export using the ExportAutoScalingGroupRecommendations or ExportEC2InstanceRecommendations actions. All export jobs created in the last seven days are returned if this parameter is omitted.
- `maxResults`: The maximum number of export jobs to return with a single request. To retrieve the remaining results, make another request with the returned NextToken value.
- `nextToken`: The token to advance to the next page of export jobs.
"""
describe_recommendation_export_jobs(; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("DescribeRecommendationExportJobs"; aws_config=aws_config)
describe_recommendation_export_jobs(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("DescribeRecommendationExportJobs", args; aws_config=aws_config)

"""
    ExportAutoScalingGroupRecommendations()

Exports optimization recommendations for Auto Scaling groups. Recommendations are exported in a comma-separated values (.csv) file, and its metadata in a JavaScript Object Notation (.json) file, to an existing Amazon Simple Storage Service (Amazon S3) bucket that you specify. For more information, see Exporting Recommendations in the Compute Optimizer User Guide. You can have only one Auto Scaling group export job in progress per AWS Region.

# Required Parameters
- `s3DestinationConfig`: An object to specify the destination Amazon Simple Storage Service (Amazon S3) bucket name and key prefix for the export job. You must create the destination Amazon S3 bucket for your recommendations export before you create the export job. Compute Optimizer does not create the S3 bucket for you. After you create the S3 bucket, ensure that it has the required permission policy to allow Compute Optimizer to write the export file to it. If you plan to specify an object prefix when you create the export job, you must include the object prefix in the policy that you add to the S3 bucket. For more information, see Amazon S3 Bucket Policy for Compute Optimizer in the Compute Optimizer user guide.

# Optional Parameters
- `accountIds`: The IDs of the AWS accounts for which to export Auto Scaling group recommendations. If your account is the master account of an organization, use this parameter to specify the member accounts for which you want to export recommendations. This parameter cannot be specified together with the include member accounts parameter. The parameters are mutually exclusive. Recommendations for member accounts are not included in the export if this parameter, or the include member accounts parameter, is omitted. You can specify multiple account IDs per request.
- `fieldsToExport`: The recommendations data to include in the export file. For more information about the fields that can be exported, see Exported files in the Compute Optimizer User Guide.
- `fileFormat`: The format of the export file. The only export file format currently supported is Csv.
- `filters`: An array of objects that describe a filter to export a more specific set of Auto Scaling group recommendations.
- `includeMemberAccounts`: Indicates whether to include recommendations for resources in all member accounts of the organization if your account is the master account of an organization. The member accounts must also be opted in to Compute Optimizer. Recommendations for member accounts of the organization are not included in the export file if this parameter is omitted. This parameter cannot be specified together with the account IDs parameter. The parameters are mutually exclusive. Recommendations for member accounts are not included in the export if this parameter, or the account IDs parameter, is omitted.
"""
export_auto_scaling_group_recommendations(s3DestinationConfig; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("ExportAutoScalingGroupRecommendations", Dict{String, Any}("s3DestinationConfig"=>s3DestinationConfig); aws_config=aws_config)
export_auto_scaling_group_recommendations(s3DestinationConfig, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("ExportAutoScalingGroupRecommendations", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("s3DestinationConfig"=>s3DestinationConfig), args)); aws_config=aws_config)

"""
    ExportEC2InstanceRecommendations()

Exports optimization recommendations for Amazon EC2 instances. Recommendations are exported in a comma-separated values (.csv) file, and its metadata in a JavaScript Object Notation (.json) file, to an existing Amazon Simple Storage Service (Amazon S3) bucket that you specify. For more information, see Exporting Recommendations in the Compute Optimizer User Guide. You can have only one Amazon EC2 instance export job in progress per AWS Region.

# Required Parameters
- `s3DestinationConfig`: An object to specify the destination Amazon Simple Storage Service (Amazon S3) bucket name and key prefix for the export job. You must create the destination Amazon S3 bucket for your recommendations export before you create the export job. Compute Optimizer does not create the S3 bucket for you. After you create the S3 bucket, ensure that it has the required permission policy to allow Compute Optimizer to write the export file to it. If you plan to specify an object prefix when you create the export job, you must include the object prefix in the policy that you add to the S3 bucket. For more information, see Amazon S3 Bucket Policy for Compute Optimizer in the Compute Optimizer user guide.

# Optional Parameters
- `accountIds`: The IDs of the AWS accounts for which to export instance recommendations. If your account is the master account of an organization, use this parameter to specify the member accounts for which you want to export recommendations. This parameter cannot be specified together with the include member accounts parameter. The parameters are mutually exclusive. Recommendations for member accounts are not included in the export if this parameter, or the include member accounts parameter, is omitted. You can specify multiple account IDs per request.
- `fieldsToExport`: The recommendations data to include in the export file. For more information about the fields that can be exported, see Exported files in the Compute Optimizer User Guide.
- `fileFormat`: The format of the export file. The only export file format currently supported is Csv.
- `filters`: An array of objects that describe a filter to export a more specific set of instance recommendations.
- `includeMemberAccounts`: Indicates whether to include recommendations for resources in all member accounts of the organization if your account is the master account of an organization. The member accounts must also be opted in to Compute Optimizer. Recommendations for member accounts of the organization are not included in the export file if this parameter is omitted. Recommendations for member accounts are not included in the export if this parameter, or the account IDs parameter, is omitted.
"""
export_ec2_instance_recommendations(s3DestinationConfig; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("ExportEC2InstanceRecommendations", Dict{String, Any}("s3DestinationConfig"=>s3DestinationConfig); aws_config=aws_config)
export_ec2_instance_recommendations(s3DestinationConfig, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("ExportEC2InstanceRecommendations", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("s3DestinationConfig"=>s3DestinationConfig), args)); aws_config=aws_config)

"""
    GetAutoScalingGroupRecommendations()

Returns Auto Scaling group recommendations. AWS Compute Optimizer generates recommendations for Amazon EC2 Auto Scaling groups that meet a specific set of requirements. For more information, see the Supported resources and requirements in the AWS Compute Optimizer User Guide.

# Optional Parameters
- `accountIds`: The IDs of the AWS accounts for which to return Auto Scaling group recommendations. If your account is the master account of an organization, use this parameter to specify the member accounts for which you want to return Auto Scaling group recommendations. Only one account ID can be specified per request.
- `autoScalingGroupArns`: The Amazon Resource Name (ARN) of the Auto Scaling groups for which to return recommendations.
- `filters`: An array of objects that describe a filter that returns a more specific list of Auto Scaling group recommendations.
- `maxResults`: The maximum number of Auto Scaling group recommendations to return with a single request. To retrieve the remaining results, make another request with the returned NextToken value.
- `nextToken`: The token to advance to the next page of Auto Scaling group recommendations.
"""
get_auto_scaling_group_recommendations(; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("GetAutoScalingGroupRecommendations"; aws_config=aws_config)
get_auto_scaling_group_recommendations(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("GetAutoScalingGroupRecommendations", args; aws_config=aws_config)

"""
    GetEC2InstanceRecommendations()

Returns Amazon EC2 instance recommendations. AWS Compute Optimizer generates recommendations for Amazon Elastic Compute Cloud (Amazon EC2) instances that meet a specific set of requirements. For more information, see the Supported resources and requirements in the AWS Compute Optimizer User Guide.

# Optional Parameters
- `accountIds`: The IDs of the AWS accounts for which to return instance recommendations. If your account is the master account of an organization, use this parameter to specify the member accounts for which you want to return instance recommendations. Only one account ID can be specified per request.
- `filters`: An array of objects that describe a filter that returns a more specific list of instance recommendations.
- `instanceArns`: The Amazon Resource Name (ARN) of the instances for which to return recommendations.
- `maxResults`: The maximum number of instance recommendations to return with a single request. To retrieve the remaining results, make another request with the returned NextToken value.
- `nextToken`: The token to advance to the next page of instance recommendations.
"""
get_ec2_instance_recommendations(; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("GetEC2InstanceRecommendations"; aws_config=aws_config)
get_ec2_instance_recommendations(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("GetEC2InstanceRecommendations", args; aws_config=aws_config)

"""
    GetEC2RecommendationProjectedMetrics()

Returns the projected utilization metrics of Amazon EC2 instance recommendations.  The Cpu and Memory metrics are the only projected utilization metrics returned when you run this action. Additionally, the Memory metric is returned only for resources that have the unified CloudWatch agent installed on them. For more information, see Enabling Memory Utilization with the CloudWatch Agent.

# Required Parameters
- `endTime`: The time stamp of the last projected metrics data point to return.
- `instanceArn`: The Amazon Resource Name (ARN) of the instances for which to return recommendation projected metrics.
- `period`: The granularity, in seconds, of the projected metrics data points.
- `startTime`: The time stamp of the first projected metrics data point to return.
- `stat`: The statistic of the projected metrics.

"""
get_ec2_recommendation_projected_metrics(endTime, instanceArn, period, startTime, stat; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("GetEC2RecommendationProjectedMetrics", Dict{String, Any}("endTime"=>endTime, "instanceArn"=>instanceArn, "period"=>period, "startTime"=>startTime, "stat"=>stat); aws_config=aws_config)
get_ec2_recommendation_projected_metrics(endTime, instanceArn, period, startTime, stat, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("GetEC2RecommendationProjectedMetrics", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("endTime"=>endTime, "instanceArn"=>instanceArn, "period"=>period, "startTime"=>startTime, "stat"=>stat), args)); aws_config=aws_config)

"""
    GetEnrollmentStatus()

Returns the enrollment (opt in) status of an account to the AWS Compute Optimizer service. If the account is the master account of an organization, this action also confirms the enrollment status of member accounts within the organization.

"""
get_enrollment_status(; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("GetEnrollmentStatus"; aws_config=aws_config)
get_enrollment_status(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("GetEnrollmentStatus", args; aws_config=aws_config)

"""
    GetRecommendationSummaries()

Returns the optimization findings for an account. For example, it returns the number of Amazon EC2 instances in an account that are under-provisioned, over-provisioned, or optimized. It also returns the number of Auto Scaling groups in an account that are not optimized, or optimized.

# Optional Parameters
- `accountIds`: The IDs of the AWS accounts for which to return recommendation summaries. If your account is the master account of an organization, use this parameter to specify the member accounts for which you want to return recommendation summaries. Only one account ID can be specified per request.
- `maxResults`: The maximum number of recommendation summaries to return with a single request. To retrieve the remaining results, make another request with the returned NextToken value.
- `nextToken`: The token to advance to the next page of recommendation summaries.
"""
get_recommendation_summaries(; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("GetRecommendationSummaries"; aws_config=aws_config)
get_recommendation_summaries(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("GetRecommendationSummaries", args; aws_config=aws_config)

"""
    UpdateEnrollmentStatus()

Updates the enrollment (opt in) status of an account to the AWS Compute Optimizer service. If the account is a master account of an organization, this action can also be used to enroll member accounts within the organization.

# Required Parameters
- `status`: The new enrollment status of the account. Accepted options are Active or Inactive. You will get an error if Pending or Failed are specified.

# Optional Parameters
- `includeMemberAccounts`: Indicates whether to enroll member accounts of the organization if the your account is the master account of an organization.
"""
update_enrollment_status(status; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("UpdateEnrollmentStatus", Dict{String, Any}("status"=>status); aws_config=aws_config)
update_enrollment_status(status, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = compute_optimizer("UpdateEnrollmentStatus", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("status"=>status), args)); aws_config=aws_config)


text:
"""
    MetaGraph{T<:Integer,Label,Graph,VertexData,EdgeData,GraphData,WeightFunction,U<:Real} <: AbstractGraph{T}

A graph type with custom vertex labels containing vertex-, edge- and graph-level metadata.

Vertex labels have type `Label`, while vertex (resp. edge, resp. graph) metadata has type `VertexData` (resp. `EdgeData`, resp. `GraphData`).
It is recommended not to set `Label` to an integer type, so as to avoid confusion between vertex labels and vertex codes (which have type `T<:Integer`).

# Fields
- `g::Graph`: underlying, data-less graph with vertex indices of type `T`
- `vertex_labels::Dict{T,Label}`: dictionary mapping vertex codes to vertex labels
- `vertex_properties::Dict{Label,Tuple{T,VertexData}}`: dictionary mapping vertex labels to vertex codes & data
- `edge_data::Dict{Tuple{Label,Label},EdgeData}`: dictionary mapping edge labels such as `(label_u, label_v)` to edge metadata
- `graph_data::GraphData`: graph metadata
- `weight_function::WeightFunction`: function defining edge weight from edge metadata
- `default_weight::U`: default weight for the edges
"""
struct MetaGraph{
    T<:Integer,Label,Graph,VertexData,EdgeData,GraphData,WeightFunction,U<:Real
} <: AbstractGraph{T}
    graph::Graph
    vertex_labels::Dict{T,Label}
    vertex_properties::Dict{Label,Tuple{T,VertexData}}
    edge_data::Dict{Tuple{Label,Label},EdgeData}
    graph_data::GraphData
    weight_function::WeightFunction
    default_weight::U
end

"""
    MetaGraph(
        g;
        Label = Symbol,
        VertexData = Nothing,
        EdgeData = Nothing,
        graph_data = nothing,
        weight_function = edge_data -> 1.0,
        default_weight = 1.0
    )

Construct an empty `MetaGraph` with the given metadata types and weights.
"""
function MetaGraph(
    graph::AbstractGraph{T};
    Label=Symbol,
    VertexData=Nothing,
    EdgeData=Nothing,
    graph_data=nothing,
    weight_function=edata -> 1.0,
    default_weight=1.0,
) where {T}
    if Label <: Integer
        @warn "Constructing a MetaGraph with integer labels is not advised."
    elseif nv(graph) > 0
        @warn "Constructing a MetaGraph with a nonempty underlying graph is not advised."
    end
    return MetaGraph(
        graph,
        Dict{T,Label}(),
        Dict{Label,Tuple{T,VertexData}}(),
        Dict{Tuple{Label,Label},EdgeData}(),
        graph_data,
        weight_function,
        default_weight,
    )
end

"""
    arrange(g, label_1, label_2)

Sort two vertex labels in a default order (useful to uniquely express undirected edges).
"""
function arrange end

function Base.zero(
    g::MetaGraph{T,Label,Graph,VertexData,EdgeData}
) where {T,Label,Graph,VertexData,EdgeData}
    return MetaGraph(
        Graph();
        Label=Label,
        VertexData=VertexData,
        EdgeData=EdgeData,
        graph_data=g.graph_data,
        weight_function=g.weight_function,
        default_weight=g.default_weight,
    )
end

Base.:(==)(x::MetaGraph, y::MetaGraph) = x.graph == y.graph

Base.copy(g::MetaGraph) = deepcopy(g)

function Base.show(
    io::IO, g::MetaGraph{<:Any,Label,<:Any,VertexData,EdgeData}
) where {Label,VertexData,EdgeData}
    return print(
        io,
        "Meta graph based on a $(g.graph) with vertex labels of type $Label, vertex metadata of type $VertexData, edge metadata of type $EdgeData, graph metadata given by $(repr(g.graph_data)), and default weight $(g.default_weight)",
    )
end