Remove old hessian code

src/geometry/functional.c

@@ -527,71 +527,6 @@ double functional_sumlist(double *list, unsigned int nel) {
     return sum;
 }
 
-/* Calculates a numerical hessian */
-static bool functional_numericalhessian(vm *v, objectmesh *mesh, elementid i, int nv, int *vid, functional_integrand *integrand, void *ref, objectsparse *hess) {
-    double eps=1e-4; // ~ (eps)^(1/4)
-    value f0;
-
-    // Finite difference rules from Abramowitz and Stegun 1972, p. 884
-    double d2xy[] = { 1.0, eps, eps, // Data for second derivative formula
-                     1.0,-eps,-eps,
-                    -1.0,-eps, eps,
-                    -1.0, eps,-eps};
-
-    double d2xx[] = { -1.0, 0, 2*eps, // Data for second derivative formula
-                     -1.0, 0, -2*eps,
-                    +16.0, 0, +eps,
-                    +16.0, 0, -eps,
-                    -30.0, 0,    0 };
-
-    double *d2,scale=1.0;
-    int neval, nevalxx=5, nevalxy=4;
-
-    // Loop over vertices in element
-    for (unsigned int j=0; j<nv; j++) {
-        for (unsigned int k=0; k<nv; k++) {
-
-            // Loop over coordinates
-            for (unsigned int l=0; l<mesh->dim; l++) {
-                double x0,y0; // x and y are the two variables currently being differentiated wrt; x=Xj[l], y=Xk[m]
-                for (unsigned int m=0; m<mesh->dim; m++) {
-                    double ff=0;
-                    matrix_getelement(mesh->vert, l, vid[j], &x0); // Get the initial values
-                    matrix_getelement(mesh->vert, m, vid[k], &y0);
-
-                    if (sparsedok_get(&hess->dok, vid[j]*mesh->dim+l, vid[k]*mesh->dim+m, &f0)) ff=MORPHO_GETFLOATVALUE(f0);
-
-                    if ((j==k) && (l==m)) { // Diagonal element
-                        d2=d2xx; neval=nevalxx;
-                        scale=1.0/(12.0*eps*eps);
-                    } else {
-                        d2=d2xy; neval=nevalxy;
-                        scale=1.0/(4.0*eps*eps);
-                    }
-
-                    double f[neval];
-                    for (int n=0; n<neval; n++) {
-                        matrix_setelement(mesh->vert, l, vid[j], x0+d2[3*n+1]);
-                        matrix_setelement(mesh->vert, m, vid[k], y0+d2[3*n+2]); // For j==l and l==m, this just overwrites the x value
-                        if (!(*integrand) (v, mesh, i, nv, vid, ref, &f[n])) return false;
-                        f[n]*=d2[3*n+0];
-                    }
-
-                    double d2f = functional_sumlist(f, neval);
-
-                    matrix_setelement(mesh->vert, l, vid[j], x0); // Reset element
-                    matrix_setelement(mesh->vert, m, vid[k], y0);
-
-                    f0=MORPHO_FLOAT(ff+d2f*scale);
-
-                    sparsedok_insert(&hess->dok, vid[j]*mesh->dim+l, vid[k]*mesh->dim+m, f0);
-                }
-            }
-        }
-    }
-    return true;
-}
-
 /* *************************
  * Map functions
  * ************************* */
@@ -770,70 +705,6 @@ bool functional_mapnumericalfieldgradientX(vm *v, functional_mapinfo *info, valu
     return ret;
 }
 
-/** Map numerical hessian over the elements
- * @param[in] v - virtual machine in use
- * @param[in] info - map info
- * @param[out] out - a matrix of integrand values
- * @returns true on success, false otherwise. Error reporting through VM. */
-bool functional_mapnumericalhessianX(vm *v, functional_mapinfo *info, value *out) {
-    objectmesh *mesh = info->mesh;
-    objectselection *sel = info->sel;
-    if (sel) UNREACHABLE("Selections not implemented in hessian");
-    grade g = info->g;
-    functional_integrand *integrand = info->integrand;
-    void *ref = info->ref;
-
-    objectsparse *conn=NULL;
-    objectsparse *hess=NULL;
-
-    bool ret=false;
-    int n=0;
-
-    varray_elementid dependencies;
-    if (info->dependencies) varray_elementidinit(&dependencies);
-
-    /* How many elements? */
-    if (!functional_countelements(v, mesh, g, &n, &conn)) return false;
-
-    /* Create the output matrix */
-    if (n>0) {
-        int ndof = mesh->vert->nrows*mesh->vert->ncols;
-        hess=object_newsparse(&ndof, &ndof);
-        if (!hess)  { morpho_runtimeerror(v, ERROR_ALLOCATIONFAILED); return false; }
-    }
-
-    int vertexid; // Use this if looping over grade 0
-    int *vid=(g==0 ? &vertexid : NULL),
-        nv=(g==0 ? 1 : 0); // The vertex indice
-
-    for (elementid i=0; i<n; i++) {
-        if (conn) sparseccs_getrowindices(&conn->ccs, i, &nv, &vid);
-        else vertexid=i;
-
-        if (vid && nv>0) {
-            if (!functional_numericalhessian(v, mesh, i, nv, vid, integrand, ref, hess)) goto functional_mapnumericalhessian_cleanup;
-
-            if (info->dependencies && // Loop over dependencies if there are any
-                (info->dependencies) (info, i, &dependencies)) {
-                for (int j=0; j<dependencies.count; j++) {
-                    if (functional_containsvertex(nv, vid, dependencies.data[j])) continue;
-                    //if (!functional_numericalremotegradient(v, info, s, dependencies.data[j], i, nv, vid, frc)) goto functional_mapnumericalhessian_cleanup;
-                }
-                dependencies.count=0;
-            }
-        }
-    }
-
-    *out = MORPHO_OBJECT(hess);
-    ret=true;
-
-functional_mapnumericalhessian_cleanup:
-    if (info->dependencies) varray_elementidclear(&dependencies);
-    if (!ret) object_free((object *) hess);
-
-    return ret;
-}
-
 /* **********************************************************************
  * Multithreaded map functions
  * ********************************************************************** */