From e2fb4bda912ddc0480ce415c7c42c12f8f992865 Mon Sep 17 00:00:00 2001 From: Tim Atherton <54725421+softmattertheory@users.noreply.github.com> Date: Thu, 11 Apr 2024 10:34:56 -0400 Subject: [PATCH 01/16] Lineintegral_hessian --- src/geometry/functional.c | 5 ++++- .../lineintegral/lineintegral_hessian.morpho | 13 +++++++++++++ .../scalarpotential_hessian.morpho | 15 +++++++++++++++ 3 files changed, 32 insertions(+), 1 deletion(-) create mode 100644 test/functionals/lineintegral/lineintegral_hessian.morpho create mode 100644 test/functionals/scalarpotential/scalarpotential_hessian.morpho diff --git a/src/geometry/functional.c b/src/geometry/functional.c index 74bd30ab..9eddda18 100644 --- a/src/geometry/functional.c +++ b/src/geometry/functional.c @@ -4243,6 +4243,8 @@ FUNCTIONAL_METHOD(LineIntegral, total, MESH_GRADE_LINE, integralref, integral_pr FUNCTIONAL_METHOD(LineIntegral, gradient, MESH_GRADE_LINE, integralref, integral_prepareref, functional_mapnumericalgradient, lineintegral_integrand, NULL, GRADSQ_ARGS, SYMMETRY_NONE); +FUNCTIONAL_METHOD(LineIntegral, hessian, MESH_GRADE_LINE, integralref, integral_prepareref, functional_mapnumericalhessian, lineintegral_integrand, NULL, GRADSQ_ARGS, SYMMETRY_NONE) + /** Initialize a LineIntegral object */ value LineIntegral_init(vm *v, int nargs, value *args) { objectinstance *self = MORPHO_GETINSTANCE(MORPHO_SELF(args)); @@ -4321,7 +4323,8 @@ MORPHO_METHOD(MORPHO_INITIALIZER_METHOD, LineIntegral_init, BUILTIN_FLAGSEMPTY), MORPHO_METHOD(FUNCTIONAL_INTEGRAND_METHOD, LineIntegral_integrand, BUILTIN_FLAGSEMPTY), MORPHO_METHOD(FUNCTIONAL_TOTAL_METHOD, LineIntegral_total, BUILTIN_FLAGSEMPTY), MORPHO_METHOD(FUNCTIONAL_GRADIENT_METHOD, LineIntegral_gradient, BUILTIN_FLAGSEMPTY), -MORPHO_METHOD(FUNCTIONAL_FIELDGRADIENT_METHOD, LineIntegral_fieldgradient, BUILTIN_FLAGSEMPTY) +MORPHO_METHOD(FUNCTIONAL_FIELDGRADIENT_METHOD, LineIntegral_fieldgradient, BUILTIN_FLAGSEMPTY), +MORPHO_METHOD(FUNCTIONAL_HESSIAN_METHOD, LineIntegral_hessian, BUILTIN_FLAGSEMPTY) MORPHO_ENDCLASS /* ---------------------------------------------- diff --git a/test/functionals/lineintegral/lineintegral_hessian.morpho b/test/functionals/lineintegral/lineintegral_hessian.morpho new file mode 100644 index 00000000..e07b4a5e --- /dev/null +++ b/test/functionals/lineintegral/lineintegral_hessian.morpho @@ -0,0 +1,13 @@ + +import meshtools +import "../numericalderivatives.morpho" + +var m = LineMesh(fn (t) [t,0,0], 0..1:0.5) + +// A line integral with only spatial dependence +var lc = LineIntegral(fn (x) (x[0]*(1-x[0]))^2) + +var h = lc.hessian(m) +var h2 = numericalhessian(lc, m) + +print (Matrix(h) - h2).norm() < 1e-5 // expect: true diff --git a/test/functionals/scalarpotential/scalarpotential_hessian.morpho b/test/functionals/scalarpotential/scalarpotential_hessian.morpho new file mode 100644 index 00000000..c411bb43 --- /dev/null +++ b/test/functionals/scalarpotential/scalarpotential_hessian.morpho @@ -0,0 +1,15 @@ + +import "../numericalderivatives.morpho" + +var m = Mesh("tetrahedron.mesh") + +fn phi(x,y,z) { + return x^2+y^2+z^2 +} + +var a = ScalarPotential(phi) + +var h = a.hessian() +var h2 = numericalhessian(a, m) + +//print (Matrix(h) - h2).norm() < 1e-5 // expect: true From dc4f2ccc1a03a2ba256132ded19c9e1e8dd949a5 Mon Sep 17 00:00:00 2001 From: Tim Atherton <54725421+softmattertheory@users.noreply.github.com> Date: Thu, 11 Apr 2024 12:43:57 -0400 Subject: [PATCH 02/16] Modernize scalar potential and add hessian --- src/geometry/functional.c | 57 ++++++------------- .../scalarpotential_hessian.morpho | 6 +- 2 files changed, 20 insertions(+), 43 deletions(-) diff --git a/src/geometry/functional.c b/src/geometry/functional.c index 9eddda18..8717bf3b 100644 --- a/src/geometry/functional.c +++ b/src/geometry/functional.c @@ -1935,10 +1935,20 @@ MORPHO_ENDCLASS static value scalarpotential_functionproperty; static value scalarpotential_gradfunctionproperty; +typedef struct { + value fn; +} scalarpotentialref; + +bool scalarpotential_prepareref(objectinstance *self, objectmesh *mesh, grade g, objectselection *sel, scalarpotentialref *ref) { + ref->fn=MORPHO_NIL; + return (objectinstance_getpropertyinterned(self, scalarpotential_functionproperty, &ref->fn) && + MORPHO_ISCALLABLE(ref->fn)); +} + /** Evaluate the scalar potential */ bool scalarpotential_integrand(vm *v, objectmesh *mesh, elementid id, int nv, int *vid, void *ref, double *out) { double *x; - value fn = *(value *) ref; + value fn = ((scalarpotentialref *) ref)->fn; value args[mesh->dim]; value ret; @@ -1975,7 +1985,6 @@ bool scalarpotential_gradient(vm *v, objectmesh *mesh, elementid id, int nv, int return false; } - /** Initialize a scalar potential */ value ScalarPotential_init(vm *v, int nargs, value *args) { objectinstance_setproperty(MORPHO_GETINSTANCE(MORPHO_SELF(args)), functional_gradeproperty, MORPHO_INTEGER(MESH_GRADE_VERTEX)); @@ -1996,25 +2005,11 @@ value ScalarPotential_init(vm *v, int nargs, value *args) { return MORPHO_NIL; } -/** Integrand function */ -value ScalarPotential_integrand(vm *v, int nargs, value *args) { - functional_mapinfo info; - value out=MORPHO_NIL; +FUNCTIONAL_METHOD(ScalarPotential, integrand, MESH_GRADE_VERTEX, scalarpotentialref, scalarpotential_prepareref, functional_mapintegrand, scalarpotential_integrand, NULL, SCALARPOTENTIAL_FNCLLBL, SYMMETRY_NONE) - if (functional_validateargs(v, nargs, args, &info)) { - value fn; - if (objectinstance_getpropertyinterned(MORPHO_GETINSTANCE(MORPHO_SELF(args)), scalarpotential_functionproperty, &fn)) { - info.g = MESH_GRADE_VERTEX; - info.integrand = scalarpotential_integrand; - info.ref = &fn; - if (MORPHO_ISCALLABLE(fn)) { - functional_mapintegrand(v, &info, &out); - } else morpho_runtimeerror(v, SCALARPOTENTIAL_FNCLLBL); - } else morpho_runtimeerror(v, VM_OBJECTLACKSPROPERTY, SCALARPOTENTIAL_FUNCTION_PROPERTY); - } - if (!MORPHO_ISNIL(out)) morpho_bindobjects(v, 1, &out); - return out; -} +FUNCTIONAL_METHOD(ScalarPotential, total, MESH_GRADE_VERTEX, scalarpotentialref, scalarpotential_prepareref, functional_sumintegrand, scalarpotential_integrand, NULL, SCALARPOTENTIAL_FNCLLBL, SYMMETRY_NONE) + +FUNCTIONAL_METHOD(ScalarPotential, hessian, MESH_GRADE_VERTEX, scalarpotentialref, scalarpotential_prepareref, functional_mapnumericalhessian, scalarpotential_integrand, NULL, SCALARPOTENTIAL_FNCLLBL, SYMMETRY_NONE) /** Evaluate a gradient */ value ScalarPotential_gradient(vm *v, int nargs, value *args) { @@ -2051,30 +2046,12 @@ value ScalarPotential_gradient(vm *v, int nargs, value *args) { return out; } -/** Total function */ -value ScalarPotential_total(vm *v, int nargs, value *args) { - functional_mapinfo info; - value out=MORPHO_NIL; - - if (functional_validateargs(v, nargs, args, &info)) { - value fn; - if (objectinstance_getpropertyinterned(MORPHO_GETINSTANCE(MORPHO_SELF(args)), scalarpotential_functionproperty, &fn)) { - info.g = MESH_GRADE_VERTEX; - info.integrand = scalarpotential_integrand; - info.ref = &fn; - if (MORPHO_ISCALLABLE(fn)) { - functional_sumintegrand(v, &info, &out); - } else morpho_runtimeerror(v, SCALARPOTENTIAL_FNCLLBL); - } else morpho_runtimeerror(v, VM_OBJECTLACKSPROPERTY, SCALARPOTENTIAL_FUNCTION_PROPERTY); - } - return out; -} - MORPHO_BEGINCLASS(ScalarPotential) MORPHO_METHOD(MORPHO_INITIALIZER_METHOD, ScalarPotential_init, BUILTIN_FLAGSEMPTY), MORPHO_METHOD(FUNCTIONAL_INTEGRAND_METHOD, ScalarPotential_integrand, BUILTIN_FLAGSEMPTY), MORPHO_METHOD(FUNCTIONAL_GRADIENT_METHOD, ScalarPotential_gradient, BUILTIN_FLAGSEMPTY), -MORPHO_METHOD(FUNCTIONAL_TOTAL_METHOD, ScalarPotential_total, BUILTIN_FLAGSEMPTY) +MORPHO_METHOD(FUNCTIONAL_TOTAL_METHOD, ScalarPotential_total, BUILTIN_FLAGSEMPTY), +MORPHO_METHOD(FUNCTIONAL_HESSIAN_METHOD, ScalarPotential_hessian, BUILTIN_FLAGSEMPTY) MORPHO_ENDCLASS /* ---------------------------------------------- diff --git a/test/functionals/scalarpotential/scalarpotential_hessian.morpho b/test/functionals/scalarpotential/scalarpotential_hessian.morpho index c411bb43..8899dc62 100644 --- a/test/functionals/scalarpotential/scalarpotential_hessian.morpho +++ b/test/functionals/scalarpotential/scalarpotential_hessian.morpho @@ -4,12 +4,12 @@ import "../numericalderivatives.morpho" var m = Mesh("tetrahedron.mesh") fn phi(x,y,z) { - return x^2+y^2+z^2 + return x^2+y^2+0.5*z^2 + x*y } var a = ScalarPotential(phi) -var h = a.hessian() +var h = a.hessian(m) var h2 = numericalhessian(a, m) -//print (Matrix(h) - h2).norm() < 1e-5 // expect: true +print (Matrix(h) - h2).norm() < 1e-6 // expect: true From 572174c6b7a5b48b4ff79dc7f60932d3c5966407 Mon Sep 17 00:00:00 2001 From: Tim Atherton <54725421+softmattertheory@users.noreply.github.com> Date: Thu, 11 Apr 2024 17:10:01 -0400 Subject: [PATCH 03/16] Parallelized hessian calculations --- src/geometry/functional.c | 227 ++++++++++++++---- .../linecurvaturesq/hessian.morpho | 21 ++ 2 files changed, 205 insertions(+), 43 deletions(-) create mode 100644 test/functionals/linecurvaturesq/hessian.morpho diff --git a/src/geometry/functional.c b/src/geometry/functional.c index 8717bf3b..bfdc0634 100644 --- a/src/geometry/functional.c +++ b/src/geometry/functional.c @@ -28,9 +28,13 @@ value functional_fieldproperty; * Utility functions * ********************************************************************** */ +double fddelta1, // = pow(MORPHO_EPS, 1.0/3.0), + fddelta2; // = pow(MORPHO_EPS, 1.0/4.0); + // Estimates the correct stepsize for cell centered finite differences -double functional_fdstepsize(double x) { - double h = pow(MORPHO_EPS, 1.0/3.0); +double functional_fdstepsize(double x, int order) { + double h = fddelta1; + if (order==2) h = fddelta2; // h should be multiplied by an estimate of the lengthscale over which f changes, // | f / f''' | ^ (1/3) @@ -473,7 +477,7 @@ bool functional_numericalgradient(vm *v, objectmesh *mesh, elementid i, int nv, matrix_getelement(mesh->vert, k, vid[j], &x0); - eps=functional_fdstepsize(x0); + eps=functional_fdstepsize(x0, 1); matrix_setelement(mesh->vert, k, vid[j], x0+eps); if (!(*integrand) (v, mesh, i, nv, vid, ref, &fp)) return false; @@ -487,36 +491,6 @@ bool functional_numericalgradient(vm *v, objectmesh *mesh, elementid i, int nv, return true; } -/* Calculates a numerical gradient for a remote vertex */ -/*static bool functional_numericalremotegradientold(vm *v, functional_mapinfo *info, objectsparse *conn, elementid remoteid, elementid i, int nv, int *vid, objectmatrix *frc) { - objectmesh *mesh = info->mesh; - double f0,fp,fm,x0,eps=1e-10; // Should use sqrt(machineeps)*(1+|x|) here - - int *rvid=(info->g==0 ? &remoteid : NULL), - rnv=(info->g==0 ? 1 : 0); // The vertex indices - - if (conn) sparseccs_getrowindices(&conn->ccs, remoteid, &rnv, &rvid); - - // Loop over vertices in element - for (unsigned int j=0; jdim; k++) { - matrix_getelement(frc, k, vid[j], &f0); - - matrix_getelement(mesh->vert, k, vid[j], &x0); - matrix_setelement(mesh->vert, k, vid[j], x0+eps); - if (!(*info->integrand) (v, mesh, remoteid, rnv, rvid, info->ref, &fp)) return false; - matrix_setelement(mesh->vert, k, vid[j], x0-eps); - if (!(*info->integrand) (v, mesh, remoteid, rnv, rvid, info->ref, &fm)) return false; - matrix_setelement(mesh->vert, k, vid[j], x0); - - matrix_setelement(frc, k, vid[j], f0+(fp-fm)/(2*eps)); - } - } - - return true; -}*/ - static bool functional_numericalremotegradient(vm *v, functional_mapinfo *info, objectsparse *conn, elementid remoteid, elementid i, int nv, int *vid, objectmatrix *frc) { objectmesh *mesh = info->mesh; double f0,fp,fm,x0,eps=1e-6; @@ -526,7 +500,7 @@ static bool functional_numericalremotegradient(vm *v, functional_mapinfo *info, matrix_getelement(frc, k, remoteid, &f0); matrix_getelement(mesh->vert, k, remoteid, &x0); - eps=functional_fdstepsize(x0); + eps=functional_fdstepsize(x0, 1); matrix_setelement(mesh->vert, k, remoteid, x0+eps); if (!(*info->integrand) (v, mesh, i, nv, vid, info->ref, &fp)) return false; @@ -768,7 +742,7 @@ bool functional_mapnumericalfieldgradientX(vm *v, functional_mapinfo *info, valu int k=field->offset[g]+id*field->psize*field->dof[g]+j; double fld=field->data.elements[k]; - eps=functional_fdstepsize(fld); + eps=functional_fdstepsize(fld, 1); field->data.elements[k]+=eps; @@ -801,7 +775,7 @@ bool functional_mapnumericalfieldgradientX(vm *v, functional_mapinfo *info, valu * @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_mapnumericalhessian(vm *v, functional_mapinfo *info, value *out) { +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"); @@ -980,8 +954,11 @@ bool functional_mapfn_elements(void *arg) { /** Dispatches tasks to threadpool */ bool functional_parallelmap(int ntasks, functional_task *tasks) { + int nthreads = morpho_threadnumber(); + if (!nthreads) nthreads=1; + if (!functional_poolinitialized) { - functional_poolinitialized=threadpool_init(&functional_pool, morpho_threadnumber()); + functional_poolinitialized=threadpool_init(&functional_pool, nthreads); if (!functional_poolinitialized) return false; } @@ -1274,7 +1251,7 @@ bool functional_numericalgrad(vm *v, objectmesh *mesh, elementid eid, elementid matrix_getelement(mesh->vert, k, i, &x0); - eps=functional_fdstepsize(x0); + eps=functional_fdstepsize(x0, 1); matrix_setelement(mesh->vert, k, i, x0+eps); if (!(*integrand) (v, mesh, eid, nv, vid, ref, &fp)) return false; matrix_setelement(mesh->vert, k, i, x0-eps); @@ -1436,14 +1413,12 @@ bool functional_numericalfieldgrad(vm *v, objectmesh *mesh, elementid eid, objec int k=field->offset[g]+i*field->psize*field->dof[g]+j; double f0=field->data.elements[k]; - eps=functional_fdstepsize(f0); + eps=functional_fdstepsize(f0, 1); field->data.elements[k]+=eps; - if (!(*integrand) (v, mesh, eid, nv, vid, ref, &fr)) return false; field->data.elements[k]=f0-eps; - if (!(*integrand) (v, mesh, eid, nv, vid, ref, &fl)) return false; field->data.elements[k]=f0; @@ -1542,6 +1517,167 @@ bool functional_mapnumericalfieldgradient(vm *v, functional_mapinfo *info, value return success; } +/* ---------------------------- + * Map numerical hessians + * ---------------------------- */ + +/** Adds a value to an element of a sparse matrix */ +bool functional_sparseaccumulate(objectsparse *A, int i, int j, double val) { + double f0 = 0.0; + value h0; + if (sparsedok_get(&A->dok, i, j, &h0)) { + if (!morpho_valuetofloat(h0, &f0)) return false; + } + + sparsedok_insert(&A->dok, i, j, MORPHO_FLOAT(f0+val)); +} + +/** Computes the contribution to the hessian of element eid with respect to vertices i and j */ +bool functional_numericalhess(vm *v, objectmesh *mesh, elementid eid, elementid i, elementid j, int nv, int *vid, functional_integrand *integrand, void *ref, objectsparse *hess) { + double x0,y0,epsx=1e-4,epsy=1e-4; + + for (unsigned int k=0; kdim; k++) { // Loop over coordinates in vertex i + matrix_getelement(mesh->vert, k, i, &x0); + epsx=functional_fdstepsize(x0, 2); + + if (i==j) { // Use a special formula for diagonal elements + double fc, fr, fl; + if (!(*integrand) (v, mesh, eid, nv, vid, ref, &fc)) return false; + + matrix_setelement(mesh->vert, k, i, x0+epsx); + if (!(*integrand) (v, mesh, eid, nv, vid, ref, &fr)) return false; + + matrix_setelement(mesh->vert, k, i, x0-epsx); + if (!(*integrand) (v, mesh, eid, nv, vid, ref, &fl)) return false; + + matrix_setelement(mesh->vert, k, i, x0); // Restore vertex to original position + + functional_sparseaccumulate(hess, i*mesh->dim+k, i*mesh->dim+k, (fr + fl - 2*fc)/(epsx*epsx)); + } + + // Loop over coordinates in vertex j + for (unsigned int l=(i==j? k+1 : k); // Detect whether we're in an off diagonal block + ldim; l++) { + double fll,frr,flr,frl; + + matrix_getelement(mesh->vert, l, j, &y0); + epsy=functional_fdstepsize(y0, 2); + + matrix_setelement(mesh->vert, k, i, x0+epsx); + matrix_setelement(mesh->vert, l, j, y0+epsy); + if (!(*integrand) (v, mesh, eid, nv, vid, ref, &frr)) return false; + + matrix_setelement(mesh->vert, l, j, y0-epsy); + if (!(*integrand) (v, mesh, eid, nv, vid, ref, &frl)) return false; + + matrix_setelement(mesh->vert, k, i, x0-epsx); + if (!(*integrand) (v, mesh, eid, nv, vid, ref, &fll)) return false; + + matrix_setelement(mesh->vert, l, j, y0+epsy); + if (!(*integrand) (v, mesh, eid, nv, vid, ref, &flr)) return false; + + matrix_setelement(mesh->vert, k, i, x0); // Restore vertices to original position + matrix_setelement(mesh->vert, l, j, y0); + + functional_sparseaccumulate(hess, i*mesh->dim+k, j*mesh->dim+l, (frr + fll - flr - frl)/(4*epsx*epsy)); + functional_sparseaccumulate(hess, j*mesh->dim+l, i*mesh->dim+k, (frr + fll - flr - frl)/(4*epsx*epsy)); + } + } + + return true; +} + +/** Computes the gradient of element id with respect to its constituent vertices and any dependencies */ +bool functional_numericalhessianmapfn(vm *v, objectmesh *mesh, elementid id, int nv, int *vid, void *ref, void *out) { + bool success=true; + functional_mapinfo *info=(functional_mapinfo *) ref; + + for (int i=0; iintegrand, info->ref, out)) return false; + } + } + + // Now handle dependencies + if (info->dependencies) { + varray_elementid dependencies; + varray_elementidinit(&dependencies); + + // Get list of vertices this element depends on + if ((info->dependencies) (info, id, &dependencies)) { + for (int j=0; jintegrand, info->ref, out)) success=false; + } + } + + varray_elementidclear(&dependencies); + } + + return success; +} + +/** Compute the hessian numerically */ +bool functional_mapnumericalhessian(vm *v, functional_mapinfo *info, value *out) { + int success=false; + int ntask=morpho_threadnumber(); + if (ntask==0) ntask = 1; + functional_task task[ntask]; + + varray_elementid imageids; + varray_elementidinit(&imageids); + + objectsparse *new[ntask]; // Create an output matrix for each thread + for (int i=0; imesh->dim*mesh_nvertices(info->mesh); + + // Create one output matrix per thread + new[i]=object_newsparse(&N, &N); + if (!new[i]) { morpho_runtimeerror(v, ERROR_ALLOCATIONFAILED); goto functional_maphessian_cleanup; } + + // Clone the vertex matrix for each thread + meshclones[i]=*info->mesh; + meshclones[i].vert=object_clonematrix(info->mesh->vert); + task[i].mesh=&meshclones[i]; + + task[i].ref=(void *) info; // Use this to pass the info structure + task[i].mapfn=functional_numericalhessianmapfn; + task[i].result=(void *) new[i]; + } + + functional_parallelmap(ntask, task); + + /* Then add up all the matrices */ + sparse_checkformat(new[0], SPARSE_CCS, true, true); + for (int i=1; isym==SYMMETRY_ADD) functional_symmetrysumforces(info->mesh, new[0]); + + // ...and return the result + *out = MORPHO_OBJECT(new[0]); + +functional_maphessian_cleanup: + // Free the temporary copies of the vertex matrices + for (int i=0; i Date: Thu, 11 Apr 2024 20:27:30 -0400 Subject: [PATCH 04/16] Correct bounds for hessian --- src/geometry/functional.c | 6 ++++-- test/functionals/length/length_hessian.morpho | 4 ++++ test/functionals/lineintegral/lineintegral_hessian.morpho | 2 +- 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/geometry/functional.c b/src/geometry/functional.c index bfdc0634..381e16ea 100644 --- a/src/geometry/functional.c +++ b/src/geometry/functional.c @@ -1530,6 +1530,7 @@ bool functional_sparseaccumulate(objectsparse *A, int i, int j, double val) { } sparsedok_insert(&A->dok, i, j, MORPHO_FLOAT(f0+val)); + return true; } /** Computes the contribution to the hessian of element eid with respect to vertices i and j */ @@ -1556,8 +1557,9 @@ bool functional_numericalhess(vm *v, objectmesh *mesh, elementid eid, elementid } // Loop over coordinates in vertex j - for (unsigned int l=(i==j? k+1 : k); // Detect whether we're in an off diagonal block + for (unsigned int l=0; //(i==j? k+1 : k); // Detect whether we're in an off diagonal block ldim; l++) { + if (i==j && k==l) continue; double fll,frr,flr,frl; matrix_getelement(mesh->vert, l, j, &y0); @@ -1580,7 +1582,7 @@ bool functional_numericalhess(vm *v, objectmesh *mesh, elementid eid, elementid matrix_setelement(mesh->vert, l, j, y0); functional_sparseaccumulate(hess, i*mesh->dim+k, j*mesh->dim+l, (frr + fll - flr - frl)/(4*epsx*epsy)); - functional_sparseaccumulate(hess, j*mesh->dim+l, i*mesh->dim+k, (frr + fll - flr - frl)/(4*epsx*epsy)); + //functional_sparseaccumulate(hess, j*mesh->dim+l, i*mesh->dim+k, (frr + fll - flr - frl)/(4*epsx*epsy)); } } diff --git a/test/functionals/length/length_hessian.morpho b/test/functionals/length/length_hessian.morpho index c78a89bf..b9fbe887 100644 --- a/test/functionals/length/length_hessian.morpho +++ b/test/functionals/length/length_hessian.morpho @@ -11,5 +11,9 @@ print abs(a.total(m) - 4*sqrt(2)) < 1e-10 print (a.gradient(m)-numericalgradient(a, m)).norm()<1e-6 // expect: true +//print Matrix(a.hessian(m)) +//print "" +//print numericalhessian(a, m, eps=1e-4) + print (Matrix(a.hessian(m))-numericalhessian(a, m, eps=1e-4)).norm() < 1e-2 // expect: true diff --git a/test/functionals/lineintegral/lineintegral_hessian.morpho b/test/functionals/lineintegral/lineintegral_hessian.morpho index e07b4a5e..407ec6e4 100644 --- a/test/functionals/lineintegral/lineintegral_hessian.morpho +++ b/test/functionals/lineintegral/lineintegral_hessian.morpho @@ -2,7 +2,7 @@ import meshtools import "../numericalderivatives.morpho" -var m = LineMesh(fn (t) [t,0,0], 0..1:0.5) +var m = LineMesh(fn (t) [t,0,0], 0..1:1) // A line integral with only spatial dependence var lc = LineIntegral(fn (x) (x[0]*(1-x[0]))^2) From 3e3e12b38114539bcfd768e1735e60ae99b82dd7 Mon Sep 17 00:00:00 2001 From: Tim Atherton <54725421+softmattertheory@users.noreply.github.com> Date: Fri, 12 Apr 2024 09:20:25 -0400 Subject: [PATCH 05/16] Correct LineCurvatureSq hessian --- src/geometry/functional.c | 12 ++++++++--- .../linecurvaturesq/gradient.morpho | 21 +++---------------- .../linecurvaturesq/hessian.morpho | 11 ++++++---- 3 files changed, 19 insertions(+), 25 deletions(-) diff --git a/src/geometry/functional.c b/src/geometry/functional.c index 381e16ea..8eab3909 100644 --- a/src/geometry/functional.c +++ b/src/geometry/functional.c @@ -1607,9 +1607,11 @@ bool functional_numericalhessianmapfn(vm *v, objectmesh *mesh, elementid id, int // Get list of vertices this element depends on if ((info->dependencies) (info, id, &dependencies)) { - for (int j=0; jintegrand, info->ref, out)) success=false; + for (int i=0; iintegrand, info->ref, out)) success=false; + } } } @@ -2638,6 +2640,8 @@ bool equielement_dependencies(functional_mapinfo *info, elementid id, varray_ele varray_elementid nbrs; varray_elementidinit(&nbrs); + // varray_elementidwrite(out, id); // EquiElement is a vertex element, and hence depends on itself + if (mesh_findneighbors(mesh, MESH_GRADE_VERTEX, id, eref->grade, &nbrs)>0) { for (unsigned int i=0; i0) { for (unsigned int i=0; i Date: Fri, 12 Apr 2024 10:30:05 -0400 Subject: [PATCH 06/16] Update functional.c --- src/geometry/functional.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/geometry/functional.c b/src/geometry/functional.c index 8eab3909..2209628f 100644 --- a/src/geometry/functional.c +++ b/src/geometry/functional.c @@ -1594,6 +1594,8 @@ bool functional_numericalhessianmapfn(vm *v, objectmesh *mesh, elementid id, int bool success=true; functional_mapinfo *info=(functional_mapinfo *) ref; + // TODO: Exploit symmetry of hessian to reduce work + for (int i=0; iintegrand, info->ref, out)) return false; From 7f7f9a57503e6d80cadebc97b8cac2bf5cfdc9b4 Mon Sep 17 00:00:00 2001 From: Tim Atherton <54725421+softmattertheory@users.noreply.github.com> Date: Fri, 12 Apr 2024 10:44:27 -0400 Subject: [PATCH 07/16] LineTorsionSq hessian --- src/geometry/functional.c | 4 +++- .../functionals/linetorsionsq/gradient.morpho | 23 ++++--------------- test/functionals/linetorsionsq/hessian.morpho | 18 +++++++++++++++ 3 files changed, 25 insertions(+), 20 deletions(-) create mode 100644 test/functionals/linetorsionsq/hessian.morpho diff --git a/src/geometry/functional.c b/src/geometry/functional.c index 2209628f..b782b16f 100644 --- a/src/geometry/functional.c +++ b/src/geometry/functional.c @@ -3009,12 +3009,14 @@ FUNCTIONAL_INIT(LineTorsionSq, MESH_GRADE_LINE) FUNCTIONAL_METHOD(LineTorsionSq, integrand, MESH_GRADE_LINE, curvatureref, curvature_prepareref, functional_mapintegrand, linetorsionsq_integrand, NULL, FUNCTIONAL_ARGS, SYMMETRY_NONE) FUNCTIONAL_METHOD(LineTorsionSq, total, MESH_GRADE_LINE, curvatureref, curvature_prepareref, functional_sumintegrand, linetorsionsq_integrand, NULL, FUNCTIONAL_ARGS, SYMMETRY_NONE) FUNCTIONAL_METHOD(LineTorsionSq, gradient, MESH_GRADE_LINE, curvatureref, curvature_prepareref, functional_mapnumericalgradient, linetorsionsq_integrand, linetorsionsq_dependencies, FUNCTIONAL_ARGS, SYMMETRY_ADD) +FUNCTIONAL_METHOD(LineTorsionSq, hessian, MESH_GRADE_LINE, curvatureref, curvature_prepareref, functional_mapnumericalhessian, linetorsionsq_integrand, linetorsionsq_dependencies, FUNCTIONAL_ARGS, SYMMETRY_ADD) MORPHO_BEGINCLASS(LineTorsionSq) MORPHO_METHOD(MORPHO_INITIALIZER_METHOD, LineTorsionSq_init, BUILTIN_FLAGSEMPTY), MORPHO_METHOD(FUNCTIONAL_INTEGRAND_METHOD, LineTorsionSq_integrand, BUILTIN_FLAGSEMPTY), MORPHO_METHOD(FUNCTIONAL_GRADIENT_METHOD, LineTorsionSq_gradient, BUILTIN_FLAGSEMPTY), -MORPHO_METHOD(FUNCTIONAL_TOTAL_METHOD, LineTorsionSq_total, BUILTIN_FLAGSEMPTY) +MORPHO_METHOD(FUNCTIONAL_TOTAL_METHOD, LineTorsionSq_total, BUILTIN_FLAGSEMPTY), +MORPHO_METHOD(FUNCTIONAL_HESSIAN_METHOD, LineTorsionSq_hessian, BUILTIN_FLAGSEMPTY) MORPHO_ENDCLASS /* ---------------------------------------------- diff --git a/test/functionals/linetorsionsq/gradient.morpho b/test/functionals/linetorsionsq/gradient.morpho index 043506da..bac61b99 100644 --- a/test/functionals/linetorsionsq/gradient.morpho +++ b/test/functionals/linetorsionsq/gradient.morpho @@ -2,6 +2,8 @@ import constants import meshtools +import "../numericalderivatives.morpho" + var np=20 var R=1 @@ -11,23 +13,6 @@ var m = LineMesh(fn (t) [R*cos(t), R*sin(t), t], 0..2*Pi:2*Pi/np, closed=false) var lc = LineTorsionSq() var grad = lc.gradient(m) +var ngrad = numericalgradient(lc, m) -var dim = grad.dimensions() -var ngrad = Matrix(dim[0], dim[1]) - -// Manually calculate the gradient -var vert = m.vertexmatrix() -var eps = 1e-8 - -for (i in 0...dim[0]) { - for (j in 0...dim[1]) { - var v = vert[i, j] - vert[i, j] = v + eps - var fp = lc.total(m) - vert[i, j] = v - eps - var fm = lc.total(m) - ngrad[i,j] = (fp-fm)/(2*eps) - } -} - -print (grad-ngrad).norm()/grad.count() < 1e-4 // expect: true +print (grad-ngrad).norm()/grad.count() < 1e-6 // expect: true diff --git a/test/functionals/linetorsionsq/hessian.morpho b/test/functionals/linetorsionsq/hessian.morpho new file mode 100644 index 00000000..c420afd4 --- /dev/null +++ b/test/functionals/linetorsionsq/hessian.morpho @@ -0,0 +1,18 @@ +// Line torsion Sq +import constants +import meshtools + +import "../numericalderivatives.morpho" + +var np=4 +var R=1 + +var m = LineMesh(fn (t) [R*cos(t), R*sin(t), t], 0..2*Pi:2*Pi/np, closed=false) + +// Create the manifold +var lc = LineTorsionSq() + +var h = lc.hessian(m) +var h2 = numericalhessian(lc, m) + +print (Matrix(h) - h2).norm()/h2.count() < 1e-5 // expect: true \ No newline at end of file From 92be96f14296476769f097c60c78808f350e73da Mon Sep 17 00:00:00 2001 From: Tim Atherton <54725421+softmattertheory@users.noreply.github.com> Date: Fri, 12 Apr 2024 10:46:36 -0400 Subject: [PATCH 08/16] Remove old hessian code --- src/geometry/functional.c | 129 -------------------------------------- 1 file changed, 129 deletions(-) diff --git a/src/geometry/functional.c b/src/geometry/functional.c index b782b16f..972f44ca 100644 --- a/src/geometry/functional.c +++ b/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; jdim; 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; mdim; 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; nvert, 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; iccs, 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; jdependencies) varray_elementidclear(&dependencies); - if (!ret) object_free((object *) hess); - - return ret; -} - /* ********************************************************************** * Multithreaded map functions * ********************************************************************** */ From 9b05a3cc00f2177147a6f421bd21deb9ea03e5a6 Mon Sep 17 00:00:00 2001 From: Tim Atherton <54725421+softmattertheory@users.noreply.github.com> Date: Tue, 21 May 2024 10:16:57 -0400 Subject: [PATCH 09/16] Correct sparse matrix addition --- src/geometry/functional.c | 31 ++++++++++++++++--- src/geometry/functional.h | 3 ++ src/linalg/sparse.h | 5 +++ test/functionals/length/length_hessian.morpho | 4 --- 4 files changed, 34 insertions(+), 9 deletions(-) diff --git a/src/geometry/functional.c b/src/geometry/functional.c index 972f44ca..e5d6d354 100644 --- a/src/geometry/functional.c +++ b/src/geometry/functional.c @@ -1505,10 +1505,13 @@ bool functional_mapnumericalhessian(vm *v, functional_mapinfo *info, value *out) varray_elementidinit(&imageids); objectsparse *new[ntask]; // Create an output matrix for each thread - for (int i=0; imesh; meshclones[i].vert=object_clonematrix(info->mesh->vert); + if (!meshclones[i].vert) { morpho_runtimeerror(v, ERROR_ALLOCATIONFAILED); goto functional_maphessian_cleanup; } task[i].mesh=&meshclones[i]; task[i].ref=(void *) info; // Use this to pass the info structure @@ -1530,15 +1534,32 @@ bool functional_mapnumericalhessian(vm *v, functional_mapinfo *info, value *out) functional_parallelmap(ntask, task); - /* Then add up all the matrices */ - sparse_checkformat(new[0], SPARSE_CCS, true, true); - for (int i=1; iccs); + new[0]->ccs = out.ccs; + } else { + morpho_runtimeerror(v, SPARSE_OPFAILEDERR); + goto functional_maphessian_cleanup; + } + } success=true; // Use symmetry actions //if (info->sym==SYMMETRY_ADD) functional_symmetrysumforces(info->mesh, new[0]); + sparsedok_clear(&new[0]->dok); // Remove dok info + // ...and return the result *out = MORPHO_OBJECT(new[0]); diff --git a/src/geometry/functional.h b/src/geometry/functional.h index 1cd26bf0..d0389f04 100644 --- a/src/geometry/functional.h +++ b/src/geometry/functional.h @@ -83,6 +83,9 @@ #define FUNC_INTEGRAND_MESH "FnctlIntMsh" #define FUNC_INTEGRAND_MESH_MSG "Method 'integrand' requires a mesh as the argument." +#define FUNC_INTEGRAND_MESH "FnctlIntMsh" +#define FUNC_INTEGRAND_MESH_MSG "Method 'integrand' requires a mesh as the argument." + #define FUNC_ELNTFND "FnctlELNtFnd" #define FUNC_ELNTFND_MSG "Mesh does not provide elements of grade %u." diff --git a/src/linalg/sparse.h b/src/linalg/sparse.h index 984b5d37..36431590 100644 --- a/src/linalg/sparse.h +++ b/src/linalg/sparse.h @@ -76,6 +76,9 @@ typedef struct { /** Gets the object as a sparse matrix */ #define MORPHO_GETSPARSE(val) ((objectsparse *) MORPHO_GETOBJECT(val)) +/** @brief Use to create static sparse matrices on the C stack. Note that the entries should be initialized */ +#define MORPHO_STATICSPARSE() { .obj.type=OBJECT_SPARSE, .obj.status=OBJECT_ISUNMANAGED, .obj.next=NULL } + objectsparse *object_newsparse(int *nrows, int *ncols); objectsparse *sparse_sparsefromarray(objectarray *array); @@ -130,6 +133,7 @@ bool sparsedok_copy(sparsedok *src, sparsedok *dest); bool sparsedok_copyat(sparsedok *src, sparsedok *dest, int row0, int col0); bool sparsedok_copymatrixat(objectmatrix *src, sparsedok *dest, int row0, int col0); bool sparsedok_copytomatrix(sparsedok *src, objectmatrix *dest, int row0, int col0); +void sparsedok_print(vm *v, sparsedok *dok); /* *************************************** * Compressed Column Storage Format @@ -149,6 +153,7 @@ bool sparseccs_doktoccs(sparsedok *in, sparseccs *out, bool copyvals); bool sparseccs_copy(sparseccs *src, sparseccs *dest); bool sparseccs_copytodok(sparseccs *src, sparsedok *dest, int row0, int col0); bool sparseccs_copytomatrix(sparseccs *src, objectmatrix *dest, int row0, int col0); +void sparseccs_print(vm *v, sparseccs *ccs); typedef enum { SPARSE_DOK, SPARSE_CCS } objectsparseformat; diff --git a/test/functionals/length/length_hessian.morpho b/test/functionals/length/length_hessian.morpho index b9fbe887..c78a89bf 100644 --- a/test/functionals/length/length_hessian.morpho +++ b/test/functionals/length/length_hessian.morpho @@ -11,9 +11,5 @@ print abs(a.total(m) - 4*sqrt(2)) < 1e-10 print (a.gradient(m)-numericalgradient(a, m)).norm()<1e-6 // expect: true -//print Matrix(a.hessian(m)) -//print "" -//print numericalhessian(a, m, eps=1e-4) - print (Matrix(a.hessian(m))-numericalhessian(a, m, eps=1e-4)).norm() < 1e-2 // expect: true From 4287bc3e7dde49c07d512d1f84065bf0e15cc9d7 Mon Sep 17 00:00:00 2001 From: Chaitanya Joshi Date: Tue, 21 May 2024 16:13:30 -0400 Subject: [PATCH 10/16] Add automated testing + CI with two threads --- .../workflows/buildandtestmultithreaded.yml | 42 +++++++++++++++++++ .gitignore | 2 +- test/test.py | 18 ++++++-- 3 files changed, 58 insertions(+), 4 deletions(-) create mode 100644 .github/workflows/buildandtestmultithreaded.yml diff --git a/.github/workflows/buildandtestmultithreaded.yml b/.github/workflows/buildandtestmultithreaded.yml new file mode 100644 index 00000000..f03e3691 --- /dev/null +++ b/.github/workflows/buildandtestmultithreaded.yml @@ -0,0 +1,42 @@ +name: Build and Test + +on: + push: + branches: [ "main" ] + pull_request: + branches: [ "main" ] + +jobs: + build: + + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v3 + - name: configure + run: | + sudo apt update + sudo apt install libsuitesparse-dev + sudo apt install liblapacke-dev + sudo apt install libunistring-dev + python -m pip install --upgrade pip + python -m pip install regex colored + - name: make + run: | + mkdir build + cd build + cmake -DCMAKE_BUILD_TYPE=Release .. + sudo make install + sudo mkdir /usr/local/lib/morpho + - name: getcli + run: | + git clone https://github.com/Morpho-lang/morpho-cli.git + cd morpho-cli + mkdir build + cd build + cmake .. + sudo make install + - name: test + run: | + cd test + python3 test.py -c -m diff --git a/.gitignore b/.gitignore index 8301e574..e09f8565 100644 --- a/.gitignore +++ b/.gitignore @@ -58,7 +58,7 @@ dkms.conf .entitlements .vscode/settings.json -test/FailedTests.txt +test/FailedTests*.txt *.png *.out manual/src/manual.lyx~ diff --git a/test/test.py b/test/test.py index 049bfe9d..81e5d994 100755 --- a/test/test.py +++ b/test/test.py @@ -183,11 +183,23 @@ def test(file,testLog,CI): # look for a command line arguement that says # this is being run for continous integration CI = False -if (len(sys.argv) > 1): - CI = sys.argv[1] == '-c' +# Also look for a command line argument that says this is being run with multiple threads +MP = False +if (len(sys.argv) == 2): + CI = sys.argv[1] == '-c' # if the argument is -c, then we are running in CI mode + MP = sys.argv[1] == '-m' # if the argument is -m, then we are running in multi-thread mode +elif (len(sys.argv) == 3): + CI = sys.argv[1] == '-c' or sys.argv[2] == '-c' + MP = sys.argv[1] == '-m' or sys.argv[2] == '-m' + +failedTestsFileName = "FailedTests.txt" +if MP: + failedTestsFileName = "FailedTestsMultiThreaded.txt" + command = "morpho6 -w2" + print("Running tests with 2 threads") files=glob.glob('**/**.'+ext, recursive=True) -with open("FailedTests.txt",'w') as testLog: +with open(failedTestsFileName,'w') as testLog: for f in files: # print(f) From 178e2e7636893dec43446480523fdd3db392854b Mon Sep 17 00:00:00 2001 From: Tim Atherton <54725421+softmattertheory@users.noreply.github.com> Date: Tue, 21 May 2024 21:14:26 -0400 Subject: [PATCH 11/16] Fix gradvector test --- src/geometry/functional.c | 9 +++++++++ test/functionals/areaintegral/gradvector.morpho | 15 +++++---------- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/geometry/functional.c b/src/geometry/functional.c index e5d6d354..e54d951e 100644 --- a/src/geometry/functional.c +++ b/src/geometry/functional.c @@ -1494,6 +1494,12 @@ bool functional_numericalhessianmapfn(vm *v, objectmesh *mesh, elementid id, int return success; } +static int _sparsecmp(const void *a, const void *b) { + objectsparse *aa = *(objectsparse **) a; + objectsparse *bb = *(objectsparse **) b; + return bb->dok.dict.count - aa->dok.dict.count; +} + /** Compute the hessian numerically */ bool functional_mapnumericalhessian(vm *v, functional_mapinfo *info, value *out) { int success=false; @@ -1534,6 +1540,8 @@ bool functional_mapnumericalhessian(vm *v, functional_mapinfo *info, value *out) functional_parallelmap(ntask, task); + qsort(new, ntask, sizeof(objectsparse *), _sparsecmp); + if (!sparse_checkformat(new[0], SPARSE_CCS, true, true)) { morpho_runtimeerror(v, SPARSE_OPFAILEDERR); goto functional_maphessian_cleanup; @@ -1541,6 +1549,7 @@ bool functional_mapnumericalhessian(vm *v, functional_mapinfo *info, value *out) /* Then add up all the matrices */ for (int i=1; idok.dict.count) continue; objectsparse out = MORPHO_STATICSPARSE(); sparsedok_init(&out.dok); sparseccs_init(&out.ccs); diff --git a/test/functionals/areaintegral/gradvector.morpho b/test/functionals/areaintegral/gradvector.morpho index 3b407ce5..6f0b70c4 100644 --- a/test/functionals/areaintegral/gradvector.morpho +++ b/test/functionals/areaintegral/gradvector.morpho @@ -10,22 +10,17 @@ mb.addface([0,1,2]) var m = mb.build() var f = Field(m, fn (x,y) Matrix([x,2*y])) -var r + +var r = [ Matrix([[1],[0]]), Matrix([[0],[2]]) ] +var out = true fn integrand(x, n) { var g = grad(n) - r = [] - for (u in g) r.append(u.clone()) + for (gg, k in g) if ((gg-r[k]).norm()>1e-4) out = false return 0 } print AreaIntegral(integrand, f).total(m) // expect: 0 -print r[0] -// expect: [ 1 ] -// expect: [ 0 ] - -print r[1] -// expect: [ 0 ] -// expect: [ 2 ] +print out // expect: true From 74ffd286dfc3acebde954a9356c22ad96ea05a0b Mon Sep 17 00:00:00 2001 From: Tim Atherton <54725421+softmattertheory@users.noreply.github.com> Date: Tue, 21 May 2024 21:49:55 -0400 Subject: [PATCH 12/16] Ensure meshes have CCS matrices --- src/geometry/functional.c | 11 ++++++++++- test/test.py | 2 +- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/src/geometry/functional.c b/src/geometry/functional.c index e54d951e..9807253b 100644 --- a/src/geometry/functional.c +++ b/src/geometry/functional.c @@ -881,6 +881,15 @@ int functional_preparetasks(vm *v, functional_mapinfo *info, int ntask, function int bins[ntask+1]; functional_binbounds(cmax, ntask, bins); + /* Ensure all mesh topology matrices have CCS */ + int maxgrade=mesh_maxgrade(info->mesh); + for (int i=0; i<=maxgrade; i++) { + for (int j=0; j<=maxgrade; j++) { + objectsparse *s = mesh_getconnectivityelement(info->mesh, i, j); + if (s) sparse_checkformat(s, SPARSE_CCS, true, false); + } + } + /* Find any image elements so they can be skipped */ functional_symmetryimagelist(info->mesh, info->g, true, imageids); if (info->field) field_addpool(info->field); @@ -1549,7 +1558,7 @@ bool functional_mapnumericalhessian(vm *v, functional_mapinfo *info, value *out) /* Then add up all the matrices */ for (int i=1; idok.dict.count) continue; + if (!new[i]->dok.dict.count) continue; objectsparse out = MORPHO_STATICSPARSE(); sparsedok_init(&out.dok); sparseccs_init(&out.ccs); diff --git a/test/test.py b/test/test.py index 049bfe9d..9bbeeb58 100755 --- a/test/test.py +++ b/test/test.py @@ -16,7 +16,7 @@ from colored import stylize # define what command to use to invoke the interpreter -command = 'morpho6' +command = 'morpho6 -w4' # define the file extension to test ext = 'morpho' From 057dcd98fdcb8de3dec26929211211dc60aeb755 Mon Sep 17 00:00:00 2001 From: Chaitanya Joshi Date: Wed, 22 May 2024 08:37:54 -0400 Subject: [PATCH 13/16] Simplify sys argv and change to 4 threads --- test/test.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/test/test.py b/test/test.py index 81e5d994..eb22fb6b 100755 --- a/test/test.py +++ b/test/test.py @@ -185,18 +185,17 @@ def test(file,testLog,CI): CI = False # Also look for a command line argument that says this is being run with multiple threads MP = False -if (len(sys.argv) == 2): - CI = sys.argv[1] == '-c' # if the argument is -c, then we are running in CI mode - MP = sys.argv[1] == '-m' # if the argument is -m, then we are running in multi-thread mode -elif (len(sys.argv) == 3): - CI = sys.argv[1] == '-c' or sys.argv[2] == '-c' - MP = sys.argv[1] == '-m' or sys.argv[2] == '-m' +for arg in sys.argv: + if arg == '-c': # if the argument is -c, then we are running in CI mode + CI = True + if arg == '-m': # if the argument is -m, then we are running in multi-thread mode + MP = True failedTestsFileName = "FailedTests.txt" if MP: failedTestsFileName = "FailedTestsMultiThreaded.txt" - command = "morpho6 -w2" - print("Running tests with 2 threads") + command += " -w4" + print("Running tests with 4 threads") files=glob.glob('**/**.'+ext, recursive=True) with open(failedTestsFileName,'w') as testLog: From 6e9c82bd2fa27fe4d38e8f3e528cd9b7690184a8 Mon Sep 17 00:00:00 2001 From: Chaitanya Joshi Date: Wed, 22 May 2024 08:42:23 -0400 Subject: [PATCH 14/16] Add dev branch for CI --- .github/workflows/buildandtestmultithreaded.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/buildandtestmultithreaded.yml b/.github/workflows/buildandtestmultithreaded.yml index f03e3691..22979cde 100644 --- a/.github/workflows/buildandtestmultithreaded.yml +++ b/.github/workflows/buildandtestmultithreaded.yml @@ -2,9 +2,9 @@ name: Build and Test on: push: - branches: [ "main" ] + branches: [ "main", "dev" ] pull_request: - branches: [ "main" ] + branches: [ "main", "dev" ] jobs: build: From c227091be3640769d94966cedd46033642bb78d6 Mon Sep 17 00:00:00 2001 From: Chaitanya Joshi Date: Wed, 22 May 2024 08:49:41 -0400 Subject: [PATCH 15/16] Update test name --- .github/workflows/buildandtestmultithreaded.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/buildandtestmultithreaded.yml b/.github/workflows/buildandtestmultithreaded.yml index 22979cde..11815d11 100644 --- a/.github/workflows/buildandtestmultithreaded.yml +++ b/.github/workflows/buildandtestmultithreaded.yml @@ -1,4 +1,4 @@ -name: Build and Test +name: TestMultiThreaded on: push: From 0af760d8fb3926a1bc00a47d890171b47b962c37 Mon Sep 17 00:00:00 2001 From: Chaitanya Joshi Date: Wed, 22 May 2024 08:51:10 -0400 Subject: [PATCH 16/16] Change MP to MT to match "MultiThreaded" --- test/test.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test.py b/test/test.py index eb22fb6b..12d047af 100755 --- a/test/test.py +++ b/test/test.py @@ -184,15 +184,15 @@ def test(file,testLog,CI): # this is being run for continous integration CI = False # Also look for a command line argument that says this is being run with multiple threads -MP = False +MT = False for arg in sys.argv: if arg == '-c': # if the argument is -c, then we are running in CI mode CI = True if arg == '-m': # if the argument is -m, then we are running in multi-thread mode - MP = True + MT = True failedTestsFileName = "FailedTests.txt" -if MP: +if MT: failedTestsFileName = "FailedTestsMultiThreaded.txt" command += " -w4" print("Running tests with 4 threads")