Open bpaniagua opened 1 month ago
https://github.com/slicersalt/py-spharm-pdm/
Two attempts at optimization with branches torch
and torch-optim
- neither work very well for one reason or another.
Currently the two core challenges:
scipy.optimize
uses dense matrices for storing Jacobian and Hessian information, so its constrained optimizers use a prohibitive amount of memory.
torch
attempts to use pytorch backward autograd to compute these instead. Runtime improves but memory usage still prohibitive.torch-optim
attempts to use torch.optim
library, however it doesn't support constrained optimizers. I attempted a constraint penalty term, but it doesn't perform well: it optimizes for the constraint or the metric, not both. There are too many degrees of freedom and the constraints are too sparse.core.initial_parameterization
that causes self-intersecting initial parameterization; since the initial conditions are wrong the optimizer fails to produce any meaningful result.The problem with the python optimizers is we need sparse and constrained optimization. I'm curious to investigate the package pyOptSparse; it claims to support this. It doesn't provide an autograd interface but I can attempt with the same technique as in the torch
branch. I'm not super confident since it's not a terribly popular library but I think it's worth looking at before diving into a custom implementation (that's the current situation with EqualAreaParametricMeshNewtonIterator
!!)
The first issue to resolve is the discrepency in initial parameterization. My current Python implementation is based on the paper, not the code. I've isolated the relevant sections of EqualAreaParametricMeshNewtonIterator
and am comparing them now. I will follow up here with details, or schedule a breakout with you and Jared if the bug is too subtle.
Note for posterity: the relevant methods in EqualAreaParametricMeshNewtonIterator
are:
start_values
invokes general workflow.generate_matrix
populates laplacian matrix for latitude problem.set_lati_rhs
populates rhs vector for latitude problem.Ax = b
; extract latitudes from x
.modify_matrix
updates the matrix for the longitude problem.set_longi_rhs
populates rhs vector for longitude problem.Ax = b
; extract longitudes from x
.I believe I located and resolved the bug: https://github.com/slicersalt/py-spharm-pdm/commit/e5a47f5a007fe81d4d313a386da8169d2273eb0d Now working on more detailed tests and automated validation to be sure, but I wanted to report this progress.
While investigating this I realize that the legacy TSurfaceNet
and new vtkSurfaceNets3D
output have different vertex ordering; then the spherical parameterization and resulting harmonics will differ. Therefore I do not think it will be possible to maintain correspondence between models generated with this rework and the legacy implementation.
It may be possible to resolve it in certain cases with ellipsoid alignment but I'm not confident. Assuming not, I plan to use model to model distance as the metric to validate.
https://github.com/Kitware/SlicerSALT/blob/master/CMakeLists.txt#L399
@allemangD could you please link here the SPHARM-PDM fork that includes your WIP on the python implementation of SPHARM-PDM?
Plan is to include the rework version as an "experimental" alternative for SPHARM-PDM.