r-barnes / faster-unmixer

A faster implementation of the sediment inverse/unmixing scheme proposed in Lipp et al (2021).
5 stars 1 forks source link

Add a way to solve for a 'continuous' upstream map, not discrete basins #24

Closed AlexLipp closed 1 year ago

AlexLipp commented 1 year ago

This adds new classes to geochem_inverse_optimize.py which allow a continuous version of the inverse problem to be solved for. Instead of returning an upstream map that is discretised into individual sub-basins it finds the best-fitting smoothest surface upstream. In this way, it allows the type of outputs from the 'original' unmixing algorithms to be generated efficiently using the novel convex/network approach. This would close issue #14.

The most significant changes are detailed below:

InverseNode

This new class defines a single node on an inversion grid. It has following attributes:

InverseGrid

This new class defines a regularly spaced, rectangular, grid of InverseNode objects. It is initiated by providing the number of rows (ny) and columns (nx) which the grid should have, i.e., the desired resolution. It also requires the map of sub-basins with their corresponding labels (i.e., the integer label of the basin) and the sample_network (this is used to match the area labels to the sample numbers). It stores basic information about the grid as attributes (ny, nx, locations of notes, area map).

In addition it also stores two dictionaries which allow easy access to the InverseNode objects which make up the grid.

  1. node_coord_dict maps the coordinates of a node to the InverseNode object. This dictionary makes it easy to calculate the 'roughness' of the grid.
  2. node_sample_dict maps a sample name to the list of InverseNode objects which lie in the sub-catchment corresponding to that sample name. This makes it easy to mix all the nodes in a catchment to compare it against observation sample points.

SampleNetworkContinuous

This new class is, by design, very similar to SampleNetwork in behaves in much the same way, with many of the same methods. It differs in that it now has as an InverseGrid attribute. Instead of optimising variables for each node in sample_network. it instead optimises the variables in each InverseNode within the InverseGrid attribute. The output is thus a continuous, smooth, grid of values

Demonstration

The new functionality added can be seen by running the new script, e.g,: python unmix_continuous_mwe.py which generates a 60x60 output grid. A crude exploration has shown the runtime to be ~quadratic w.r.t to the number of nodes in the grid.

The output from this new implementation is compared to the original 'discrete' below for the upstream concentration of Mg

Continuous

continuous

Discrete

discrete

AlexLipp commented 1 year ago

@r-barnes One thing I felt was that SampleNetworkContinuous could possibly be a sub-class of SampleNetwork inheriting/modifying some of the attributes and methods. However sub-classes are not something I really have much familiarity with, so I have left it for now. What do you think? It could make things simpler if we ever modify SampleNetwork in the future.

r-barnes commented 1 year ago

@AlexLipp : Is it possible that you could make use of the existing SampleNetwork class without further change by modeling the continuous grid as a sample network?

r-barnes commented 1 year ago

@AlexLipp : Also, would this code have any significant simplifications if we were doing the C++ stuff in Python instead?

AlexLipp commented 1 year ago

@r-barnes I've adapted the SampleNetwork class so that it now takes a continuous flag which indicates whether to solve 'continuously' or 'discretely'. This flag changes the init behaviour only generating an InverseGrid if continuous is True. Similarly, the behaviour of various methods is now slightly different dependent on the value of continuous.

I've also simplified the InverseGrid and InverseNode classes as you suggest.

AlexLipp commented 1 year ago

The minor changes suggested above have also been implemented.