// TODO: Time to finish the bells and whistles here
// TODO: Use sparse matrix representation for Jacobian, (this is lower priority since typical uses will not produce very large matrices)
// TODO: Add residual instrumentation to tests
// TODO: Try extracting sub matrix and solving (Certain constraints are known to eliminate variables (Radius, Fixation), if the matrix can be permuted so that a submatrix on the block diagonal is lower-triangular, then a subsystem of equations can be solved independently of the remaining equations. Moreover, if the jacobian can be permuted into a block lower triangular form, subsets of equations can be solved blockwise)
// TODO: Extract newton_update over Verticies, Curves, and Constraints into a private method
// TODO: Extract armijo_update into a private method
std::cerr<<"Variables = "<<NumVariables<<", Equations = "<<NumEquations<<". Number of variables is less than the number of equations. The sketch may be over-constrained."<<std::endl;
// TODO: Try extracting sub matrix and solving (Certain constraints are known to eliminate variables (Radius, Fixation), if the matrix can be permuted so that a submatrix on the block diagonal is lower-triangular, then a subsystem of equations can be solved independently of the remaining equations. Moreover, if the jacobian can be permuted into a block lower triangular form, subsets of equations can be solved blockwise)
// TODO: Extract newton_update over Verticies, Curves, and Constraints into a private method
// TODO: Extract armijo_update into a private method
// TODO: Factor out rank deficiency perturbation to be independent of QR/LU/LDL^T choice