diff --git a/docs/development/usability-plan.md b/docs/development/usability-plan.md index e14da21..4cbbc62 100644 --- a/docs/development/usability-plan.md +++ b/docs/development/usability-plan.md @@ -527,12 +527,12 @@ Design: Tasks: -- [ ] `main/interpolation_size.py` with `summarise_data` and `suggest_nelements`. -- [ ] Setting, settings page and preference test. -- [ ] Model manager: one method for the interpolator arguments, used in all +- [x] `main/interpolation_size.py` with `summarise_data` and `suggest_nelements`. +- [x] Setting, settings page and preference test. +- [x] Model manager: one method for the interpolator arguments, used in all build paths. Fix the use of the class defaults. -- [ ] Feature panel: show the number, and an "Automatic" check box. -- [ ] Docs: say how the number is chosen, and what the user can change. +- [x] Feature panel: show the number, and an "Automatic" check box. +- [x] Docs: say how the number is chosen, and what the user can change. Acceptance: a model with 20 contact points and a model with 5 000 points get different numbers of elements, both inside the limits. A saved fixed number is @@ -636,13 +636,34 @@ user can select any combination that is in the table below. frame). This adds the axial surface constraint (`fold_orientation`). 2. **Fold axis.** The user selects the source of the fold axis: - constant: plunge and azimuth (as now); - - average intersection lineation of the folded feature and the axial - surface (`av_fold_axis`, as now); - - lineation data: a layer of fold axis or intersection lineation - measurements. The plugin fits the fold axis rotation angle to - coordinate 1 (the axis S-plot). + - lineations. The user selects one of two lineation sources: + 1. **Lineation layer**: a point layer (for example a shapefile) with + measured fold axes or intersection lineations. The user selects the + plunge field and the trend (plunge direction) field. + 2. **Calculated intersection lineations**: the plugin calculates a + lineation at each orientation point of the folded feature. The + lineation is the intersection of the folded foliation and the axial + foliation (coordinate 0 of the fold frame) at that point + (`FoldFrame.calculate_intersection_lineation`). This source needs an + axial surface. + + Then the user selects how the plugin uses the lineations: + - **average**: the fold axis is the mean of the lineations, and it is + constant in the model (`av_fold_axis` for the calculated lineations, + as now); + - **fit**: the plugin fits the fold axis rotation angle of the + lineations to coordinate 1 (the axis S-plot), so that the fold axis + can change in the model. This needs an axial surface. This adds the fold axis constraint (`fold_axis_w`). + + | Lineation source | Average | Fit (axis S-plot) | Needs an axial surface | + |---|---|---|---| + | Lineation layer | yes | yes | only for "fit" | + | Calculated intersection lineations | yes | yes | yes | + + The UI shows the lineations of the selected source in the 3D view and in + the axis S-plot, so that the user can compare the two sources. 3. **S-plot.** The user controls the rotation angle profiles: the profile type (Fourier series, trigonometric), the wavelength, and fixed values for the profile parameters. Without this control, the plugin fits the @@ -652,11 +673,11 @@ user can select any combination that is in the table below. | Axial surface | Fold axis | S-plot | Result | |---|---|---|---| | - | - | - | A standard foliation. No fold constraint. | -| x | - | - | Fold frame and DFI. Axial surface constraint on. The fold axis is the average intersection lineation, but the fold axis constraint is off (`fold_axis_w = None`). Automatic profiles. | -| - | x | - | No fold frame. The plugin adds the fold axis as tangent constraints on a regular grid in the bounding box (gradient . axis = 0). The interpolator of the feature does not change. Only a constant fold axis is possible. | +| x | - | - | Fold frame and DFI. Axial surface constraint on. The fold axis is the average of the calculated intersection lineations, but the fold axis constraint is off (`fold_axis_w = None`). Automatic profiles. | +| - | x | - | No fold frame. The plugin adds the fold axis as tangent constraints on a regular grid in the bounding box (gradient . axis = 0). The interpolator of the feature does not change. Only a constant fold axis is possible: plunge and azimuth, or the average of a lineation layer. | | x | x | - | Fold frame and DFI. Both constraints on. Automatic profiles. | | x | - | x | As "axial surface only", but with the limb profile of the user. | -| x | x | x | All constraints on. The limb profile of the user. With "lineation data", also the axis profile of the user. | +| x | x | x | All constraints on. The limb profile of the user. With lineations and "fit", also the axis profile of the user. | | - | - | x, or - x x | Not possible. An S-plot needs a fold frame coordinate. The S-plot check box is disabled until the user selects an axial surface. The tooltip says why. | Rules: @@ -719,8 +740,11 @@ with orientations: `s2`, `s1` and `s0`. stratigraphic group: - `fold_event`: the name of the fold event, or `None`; - `axial_surface`: on or off, and `fold_orientation` weight; - - `fold_axis`: off, `constant` (plunge, azimuth), `average`, or `data` - (layer dicts), and `fold_axis_w` weight; + - `fold_axis`: `source` (off, `constant` or `lineations`), the plunge and + azimuth for `constant`, the `fold_axis_w` weight, and for `lineations`: + - `lineation_source`: `layer` (layer dict with the plunge and trend + fields) or `intersection` (calculated); + - `use`: `average` or `fit`; - `splot`: off or on; for the limb and for the axis profile: the type, the wavelength (or "automatic"), and the fixed parameters; - `fold_normalisation`, `fold_norm`, `fold_regularisation`. @@ -741,7 +765,8 @@ with orientations: `s2`, `s1` and `s0`. The S-plot needs the fold frame before the folded feature is built. Thus: - "Calculate rotation angles" builds the fold frame (and the fold events that - fold it) if it is not current. Then it calculates the rotation angles of the + fold it) if it is not current. Then it calculates the lineations (for the + "calculated intersection lineations" source) and the rotation angles of the data of the feature. It does not build the folded feature. Run it in a background task with progress. - The S-plot panel shows the data points, the fitted curve, the S-variogram @@ -791,8 +816,15 @@ Do the parts in this order. Each part is a separate pull request. the three controls with their check boxes and "Advanced" weights. - [ ] Apply the rules of the combination table. Disable the S-plot control when there is no axial surface. -- [ ] Fold axis "lineation data": a layer picker with the trend and plunge - fields. +- [ ] Fold axis "lineations": a choice of the two lineation sources. For + "lineation layer", a point layer picker with the plunge and trend + fields. For "calculated intersection lineations", no input (disabled + without an axial surface). +- [ ] Fold axis "average" or "fit" for the lineations. Disable "fit" without + an axial surface. +- [ ] Convert the plunge and trend of the layer to vectors, and give them + with their points (N x 6) to the fold axis calculation + (`main/fold_spec.py`). - [ ] Fold axis without axial surface: make the tangent constraints on a grid (`main/fold_spec.py`). The grid step comes from the bounding box and the number of elements. @@ -827,6 +859,11 @@ Tests: F2 -> F1 -> S0; a cycle is an error; each row of the combination table gives the correct arguments (a control that is off gives `None`); the S-plot is refused without an axial surface; old specs with `folded_feature_name` load. +- Unit tests for the two lineation sources: the plunge and trend of a layer + give the correct vectors; the calculated intersection lineation of a known + folded foliation and a known axial foliation is their cross product; the + "average" of both sources gives the same axis for the same data; "fit" and + the calculated source are refused without an axial surface. - A QGIS test that builds the refolded fold from the test layers and compares the S0 scalar field on a coarse grid with the result of the LoopStructural calls of the example (same data, same arguments). @@ -837,6 +874,22 @@ depend on 7.1. 7.4 depends on 7.3 (the S-plot control). 7.5 comes last. Phase 7 depends on 6.1 (the workflow mode), 6.3 (`SectionStack`) and 6.4 (the number of elements of each feature, which is also used for the fold frames). +### Phase 8: Advanced 3D viewer + +A second 3D viewer for hard 3D problems. It is a standalone application +(Rust, Bevy) with a live link to QGIS. The PyVista viewer stays as the default +viewer. The new viewer shows geological data objects that are similar to +those of Geoscience ANALYST (`geoh5` types): points, curves, surfaces, +sections, block models, drillholes and orientations. + +This phase is large, and most of the work is in a separate repository. The +tasks, the protocol, the risks and the open questions are in the +[viewer plan](viewer-plan.md). + +Acceptance: from step 5, a user opens the advanced viewer. The viewer shows +the model and its input data, and updates after each build. A pick in the +viewer selects the feature in the dock and shows the point on the map. + ## Risks - **Large UI change.** Users of the current version must learn the new layout. @@ -897,3 +950,20 @@ elements of each feature, which is also used for the fold frames). Does a fault need a different rule from a foliation? 8. (6.1) Must a manual unconformity or fold in the "constraints" mode use the column? Recommendation: no. It uses only the features of that mode. +9. (7) How does the plugin give the lineations of a layer to LoopStructural? + `FoldFrame.calculate_fold_axis_rotation` has a `fold_axis` argument (an + N x 6 array of points and lineations), but + `FoldedFeatureBuilder.set_fold_axis` does not use it, and there is no build + argument for it. "Average" of a lineation layer is easy: the plugin + calculates the mean and gives it as `fold_axis`. For "fit", either the + plugin calls `calculate_fold_axis_rotation` itself and sets + `fold.fold_axis_rotation`, or LoopStructural gets a build argument for the + lineations. Recommendation: add the build argument to LoopStructural, and + use the plugin call until that version is released. +10. (7) Is "fold axis without axial surface" (tangent constraints on a grid) + good enough, or must the plugin make a simple fold frame from the fold + axis? Recommendation: tangent constraints first. Compare the two on the + single fold example (`load_noddy_single_fold`). +11. (7) Must a fault that cuts a fold frame also cut the folded features? + LoopStructural calculates the rotation angles in the restored space. + Recommendation: yes, use the same faults. Test it before 7.5. diff --git a/docs/usage/interface.md b/docs/usage/interface.md index 54611f7..a9709a3 100644 --- a/docs/usage/interface.md +++ b/docs/usage/interface.md @@ -115,4 +115,11 @@ Units that don't match the stratigraphic column will have null values, helping y Once the layers have been selected, stratigraphic column defined and the fault topology relationships set, the LoopStructural model can be initialised. Initialise model will create a LoopStructural model with all of the geological features in the model. For each feature in the model the number of interpolation elements (degrees of freedom), the weighting of the regularisation, contact points and orientation weight can be changed. + +#### Number of interpolation elements + +By default the number of elements is "Automatic". The plugin chooses it for each feature from the data of the feature: more constraints, more surfaces and a wider spread of orientations give more elements. The number is rounded to 1 000 and kept between 5 000 and 250 000. The panel of the feature shows the number that was used. + +- To keep a number for one feature, clear the "Automatic" check box in the panel of the feature and enter the number. The feature keeps it, and the project file saves it, until you select "Automatic" again. +- To use one fixed number for all features, clear "Automatic" in the plugin settings and enter the number. Features that you set by hand keep their own number. ![Model Parameters](../static/model-setup.png) diff --git a/loopstructural/gui/dlg_settings.py b/loopstructural/gui/dlg_settings.py index 439768a..d719475 100644 --- a/loopstructural/gui/dlg_settings.py +++ b/loopstructural/gui/dlg_settings.py @@ -81,6 +81,10 @@ def __init__(self, parent): if hasattr(self, "btn_open_debug_directory"): self.btn_open_debug_directory.pressed.connect(self._open_debug_directory) + self.n_elements_auto_check_box.toggled.connect( + lambda checked: self.n_elements_spin_box.setEnabled(not checked) + ) + # load previously saved settings self.load_settings() @@ -92,6 +96,7 @@ def apply(self): settings.debug_mode = self.opt_debug.isChecked() settings.separate_dock_widgets = self.opt_separate_dock_widgets.isChecked() settings.interpolator_nelements = self.n_elements_spin_box.value() + settings.interpolator_nelements_auto = self.n_elements_auto_check_box.isChecked() settings.interpolator_npw = self.npw_spin_box.value() settings.interpolator_cpw = self.cpw_spin_box.value() settings.interpolator_regularisation = self.regularisation_spin_box.value() @@ -121,6 +126,8 @@ def load_settings(self): self.lbl_version_saved_value.setText(settings.version) # self.interpolator_type_combo.setCurrentText(settings.interpolator_type) self.n_elements_spin_box.setValue(settings.interpolator_nelements) + self.n_elements_auto_check_box.setChecked(settings.interpolator_nelements_auto) + self.n_elements_spin_box.setEnabled(not settings.interpolator_nelements_auto) self.regularisation_spin_box.setValue(settings.interpolator_regularisation) self.cpw_spin_box.setValue(settings.interpolator_cpw) self.npw_spin_box.setValue(settings.interpolator_npw) diff --git a/loopstructural/gui/dlg_settings.ui b/loopstructural/gui/dlg_settings.ui index 477f844..6a7c50e 100644 --- a/loopstructural/gui/dlg_settings.ui +++ b/loopstructural/gui/dlg_settings.ui @@ -280,6 +280,16 @@ + + + + Choose the number of elements of each feature from its data. If this is off, every feature uses the fixed number. + + + Automatic + + + diff --git a/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py b/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py index 675eebc..300f59b 100644 --- a/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py +++ b/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py @@ -4,6 +4,7 @@ from qgis.gui import QgsCollapsibleGroupBox, QgsMapLayerComboBox from qgis.PyQt.QtCore import Qt from qgis.PyQt.QtWidgets import ( + QCheckBox, QComboBox, QDoubleSpinBox, QFormLayout, @@ -134,8 +135,16 @@ def __init__(self, parent=None, *, feature=None, model_manager=None, data_manage self.n_elements_spinbox.setRange(100, 1000000) self.n_elements_spinbox.setValue(self.getNelements(feature)) self.n_elements_spinbox.setPrefix("Number of Elements: ") + self.n_elements_auto_check = QCheckBox("Automatic") + self.n_elements_auto_check.setToolTip( + "Choose the number of elements from the data of the feature. " + "Clear the check box to keep the number that you enter." + ) + self.n_elements_auto_check.setChecked(not self._has_nelements_override()) + self.n_elements_spinbox.setEnabled(not self.n_elements_auto_check.isChecked()) self.n_elements_spinbox.valueChanged.connect(self.updateNelements) + self.n_elements_auto_check.toggled.connect(self._on_nelements_auto_toggled) table_group_box = QgsCollapsibleGroupBox('Data Layers') self.layer_table = LayerSelectionTable( @@ -160,6 +169,7 @@ def __init__(self, parent=None, *, feature=None, model_manager=None, data_manage form_layout = QFormLayout() form_layout.addRow(self.interpolator_type_label, self.interpolator_type_combo) form_layout.addRow("Number of Elements:", self.n_elements_spinbox) + form_layout.addRow("", self.n_elements_auto_check) form_layout.addRow('Regularisation', self.regularisation_spin_box) form_layout.addRow('Contact points weight', self.cpw_spin_box) form_layout.addRow('Orientation point weight', self.npw_spin_box) @@ -172,6 +182,7 @@ def __init__(self, parent=None, *, feature=None, model_manager=None, data_manage summary=lambda: ( f"{self.interpolator_type_combo.currentText()}, " f"{int(self.n_elements_spinbox.value())} elements" + f"{' (automatic)' if self.n_elements_auto_check.isChecked() else ''}" ), ) self.layout.add_section(self._build_preview_widget(), 'preview', 'Preview', collapsed=True) @@ -565,8 +576,38 @@ def _on_bounding_box_updated(self, bounding_box): except Exception: pass - def updateNelements(self, value): - """Update the number of elements in the feature's interpolator.""" + def _has_nelements_override(self): + """Return True if the user set the number of elements of this feature.""" + manager = self.model_manager + name = getattr(self.feature, 'name', None) + return manager is not None and name in getattr(manager, 'nelements_overrides', {}) + + def _on_nelements_auto_toggled(self, automatic): + """Use the automatic number again, or keep the number that is shown.""" + manager = self.model_manager + name = getattr(self.feature, 'name', None) + self.n_elements_spinbox.setEnabled(not automatic) + if manager is None or name is None: + return + if automatic: + manager.set_nelements_override(name, None) + used = manager.nelements_used.get(name) + if used is not None: + self.n_elements_spinbox.blockSignals(True) + self.n_elements_spinbox.setValue(used[0]) + self.n_elements_spinbox.blockSignals(False) + self.updateNelements(used[0], keep=False) + else: + manager.set_nelements_override(name, int(self.n_elements_spinbox.value())) + + def updateNelements(self, value, keep=True): + """Update the number of elements in the feature's interpolator. + + With `keep`, the feature keeps this number in the next builds, until + the user selects "Automatic". + """ + if keep and self.model_manager is not None and self.feature is not None: + self.model_manager.set_nelements_override(self.feature.name, int(value)) if self.feature: if issubclass(type(self.feature), StructuralFrame): for i in range(3): @@ -587,6 +628,11 @@ def updateNelements(self, value): def getNelements(self, feature): """Get the number of elements from the feature's interpolator.""" if feature: + used = getattr(self.model_manager, 'nelements_used', {}).get( + getattr(feature, 'name', None) + ) + if used is not None: + return used[0] if issubclass(type(feature), StructuralFrame): return feature[0].interpolator.n_elements elif feature.interpolator is not None: diff --git a/loopstructural/main/data_manager.py b/loopstructural/main/data_manager.py index 7e7fd7d..281bd18 100644 --- a/loopstructural/main/data_manager.py +++ b/loopstructural/main/data_manager.py @@ -1848,6 +1848,7 @@ def save_state(self, filepath): state['manual_foliations'] = self._model_manager.manual_foliations_to_dict() state['detached_features'] = self._model_manager.detached_to_dict() state['parametric_faults'] = self._model_manager.parametric_faults_to_dict() + state['nelements_overrides'] = self._model_manager.nelements_overrides_to_dict() with open(path, 'w') as f: json.dump(state, f, indent=2) @@ -1882,6 +1883,7 @@ def load_state(self, filepath): self._model_manager.manual_foliations_from_dict(state.get('manual_foliations', {})) self._model_manager.detached_from_dict(state.get('detached_features', {})) self._model_manager.parametric_faults_from_dict(state.get('parametric_faults', {})) + self._model_manager.nelements_overrides_from_dict(state.get('nelements_overrides', {})) # the data was just read from the layers self._changed_layer_ids.clear() self.refresh_layer_watchers() diff --git a/loopstructural/main/interpolation_size.py b/loopstructural/main/interpolation_size.py new file mode 100644 index 0000000..c01819b --- /dev/null +++ b/loopstructural/main/interpolation_size.py @@ -0,0 +1,132 @@ +"""Choose the number of interpolation elements from the data of a feature. + +The module does not import QGIS, so the tests run in the fast tests/unit/ job. +""" + +from dataclasses import dataclass + +import numpy as np + +MIN_NELEMENTS = 5_000 +MAX_NELEMENTS = 250_000 +ROUND_TO = 1_000 +ELEMENTS_PER_EQUATION = 25 +MAX_SURFACES = 10 +SURFACE_STEP = 0.1 + + +@dataclass +class DataSummary: + """The measures of the data of one feature. + + `spread` is 0 for parallel planes and 1 for orientations in all + directions. + """ + + n_value: int = 0 + n_orientation: int = 0 + n_surfaces: int = 0 + spread: float = 0.0 + + +def _normals_from_strike_dip(strike, dip) -> np.ndarray: + """Return the unit normals of planes (right-hand rule strike, dip in degrees).""" + strike = np.radians(np.asarray(strike, dtype=float)) + dip = np.radians(np.asarray(dip, dtype=float)) + dip_direction = strike + np.pi / 2 + return np.column_stack( + [ + -np.sin(dip) * np.sin(dip_direction), + -np.sin(dip) * np.cos(dip_direction), + np.cos(dip), + ] + ) + + +def orientation_spread(normals) -> float: + """Return the spread of orientations: 0 (parallel) to 1 (all directions). + + The sign of a normal does not change the result: the measure uses the + largest eigenvalue of the mean orientation tensor, which is 1 for parallel + planes and 1/3 for uniform directions. + """ + normals = np.asarray(normals, dtype=float).reshape(-1, 3) + normals = normals[np.all(np.isfinite(normals), axis=1)] + lengths = np.linalg.norm(normals, axis=1) + normals = normals[lengths > 0] / lengths[lengths > 0, None] + if len(normals) < 2: + return 0.0 + tensor = normals.T @ normals / len(normals) + largest = np.linalg.eigvalsh(tensor)[-1] + return float(np.clip((1.0 - largest) / (2.0 / 3.0), 0.0, 1.0)) + + +def _has(df, columns) -> np.ndarray: + """Return a mask of the rows where all `columns` exist and are not NaN.""" + mask = np.ones(len(df), dtype=bool) + for column in columns: + if column not in df.columns: + return np.zeros(len(df), dtype=bool) + mask &= df[column].notna().to_numpy() + return mask + + +def summarise_data(df) -> DataSummary: + """Return the `DataSummary` of a feature data frame (or None for no data). + + Value, interface and inequality rows are value constraints. Rows with a + normal (nx, ny, nz or gx, gy, gz), a strike and a dip, or a tangent (tx, + ty, tz) are orientation constraints. + """ + if df is None or len(df) == 0: + return DataSummary() + value = _has(df, ['val']) | _has(df, ['interface']) | _has(df, ['l', 'u']) + normal = _has(df, ['nx', 'ny', 'nz']) + gradient = _has(df, ['gx', 'gy', 'gz']) & ~normal + strike_dip = _has(df, ['strike', 'dip']) & ~normal & ~gradient + tangent = _has(df, ['tx', 'ty', 'tz']) + orientation = normal | gradient | strike_dip | tangent + + vectors = [] + if normal.any(): + vectors.append(df.loc[normal, ['nx', 'ny', 'nz']].to_numpy(dtype=float)) + if gradient.any(): + vectors.append(df.loc[gradient, ['gx', 'gy', 'gz']].to_numpy(dtype=float)) + if strike_dip.any(): + vectors.append( + _normals_from_strike_dip( + df.loc[strike_dip, 'strike'].to_numpy(), df.loc[strike_dip, 'dip'].to_numpy() + ) + ) + spread = orientation_spread(np.vstack(vectors)) if vectors else 0.0 + + n_surfaces = 0 + if _has(df, ['val']).any(): + n_surfaces = int(df.loc[_has(df, ['val']), 'val'].nunique()) + if 'interface' in df.columns and df['interface'].notna().any(): + n_surfaces = max(n_surfaces, int(df['interface'].nunique())) + return DataSummary( + n_value=int(value.sum()), + n_orientation=int(orientation.sum()), + n_surfaces=n_surfaces, + spread=spread, + ) + + +def suggest_nelements(summary: DataSummary) -> int: + """Return the number of elements for a feature with the data `summary`. + + elements = equations x 25 x surface factor x spread factor, where an + orientation is two equations, the surface factor is 1 + 0.1 for each + surface after the first (at most 10 surfaces) and the spread factor is + 1 + spread. The result is rounded to 1 000 and limited to 5 000 .. 250 000. + """ + equations = summary.n_value + 2 * summary.n_orientation + if equations <= 0: + return MIN_NELEMENTS + surfaces = min(max(summary.n_surfaces, 1), MAX_SURFACES) + surface_factor = 1 + SURFACE_STEP * (surfaces - 1) + spread_factor = 1 + min(max(summary.spread, 0.0), 1.0) + elements = equations * ELEMENTS_PER_EQUATION * surface_factor * spread_factor + elements = int(round(elements / ROUND_TO)) * ROUND_TO + return int(min(max(elements, MIN_NELEMENTS), MAX_NELEMENTS)) diff --git a/loopstructural/main/model_manager.py b/loopstructural/main/model_manager.py index a241295..cc59778 100644 --- a/loopstructural/main/model_manager.py +++ b/loopstructural/main/model_manager.py @@ -27,11 +27,12 @@ from LoopStructural.utils.observer import Observable from LoopStructural import GeologicalModel -from loopstructural.toolbelt.preferences import PlgSettingsStructure +from loopstructural.toolbelt.preferences import PlgOptionsManager from ..main import constraints, parametric_fault from ..main.data_types import FaultEntry, StratigraphyEntry from ..main.helpers import qgisAttributeIsNone +from ..main.interpolation_size import suggest_nelements, summarise_data from ..main.workflow_mode import WORKFLOW_MODE_CONSTRAINTS, WORKFLOW_MODE_MAP, WORKFLOW_MODES @@ -236,6 +237,13 @@ def __init__(self, debug_manager=None): self._cancel_requested = False # How the user builds the model; see `set_workflow_mode`. self.workflow_mode = WORKFLOW_MODE_MAP + # feature name -> number of elements that the user set for the + # feature. A feature that is not here uses the automatic number or the + # fixed number of the settings; see `interpolator_arguments`. + self.nelements_overrides: Dict[str, int] = {} + # feature name -> (number of elements, True if automatic) of the last + # build, for the feature panel. + self.nelements_used: Dict[str, tuple] = {} @property def uses_map_data(self) -> bool: @@ -309,6 +317,8 @@ def reset(self): self.detached = {} self.extra_constraints = {} self.parametric_faults = {} + self.nelements_overrides = {} + self.nelements_used = {} self.dem_function = lambda x, y: 0 self._topology_dirty = False self._data_dirty = False @@ -977,14 +987,53 @@ def _build_domain_fault_boundary(self, fault_name, groupname, flipped=False): self.model.data = data_for_fault self.model.create_and_add_domain_fault( fault_name, - nelements=PlgSettingsStructure.interpolator_nelements, - npw=PlgSettingsStructure.interpolator_npw, - cpw=PlgSettingsStructure.interpolator_cpw, - regularisation=PlgSettingsStructure.interpolator_regularisation, + **self.interpolator_arguments(fault_name, data_for_fault), ) self._strip_domain_fault_region_from_boundary_above(fault_name) return True + def interpolator_arguments(self, name, data) -> dict: + """Return `nelements`, `npw`, `cpw` and `regularisation` for the + feature `name` that is built from the data frame `data`. + + The number of elements is, in this order: the number that the user set + for the feature, the number from the data (if the saved setting + "automatic" is on), or the saved fixed number. All values come from the + saved settings, not from the class defaults. + """ + settings = PlgOptionsManager.get_plg_settings() + automatic = False + if self.nelements_overrides.get(name): + nelements = int(self.nelements_overrides[name]) + elif settings.interpolator_nelements_auto: + nelements = suggest_nelements(summarise_data(data)) + automatic = True + else: + nelements = int(settings.interpolator_nelements) + self.nelements_used[name] = (nelements, automatic) + return { + 'nelements': nelements, + 'npw': settings.interpolator_npw, + 'cpw': settings.interpolator_cpw, + 'regularisation': settings.interpolator_regularisation, + } + + def nelements_overrides_to_dict(self) -> dict: + """Return `nelements_overrides` in a form that `json.dump` can write.""" + return {name: int(n) for name, n in self.nelements_overrides.items()} + + def nelements_overrides_from_dict(self, overrides: dict): + """Replace `nelements_overrides` with the numbers from `nelements_overrides_to_dict`.""" + self.nelements_overrides = {name: int(n) for name, n in (overrides or {}).items()} + + def set_nelements_override(self, name: str, nelements: Optional[int]): + """Keep `nelements` for the feature `name` in the next builds, or use + the automatic (or fixed) number again if `nelements` is None.""" + if nelements is None: + self.nelements_overrides.pop(name, None) + else: + self.nelements_overrides[name] = int(nelements) + def _strip_domain_fault_region_from_boundary_above(self, fault_name): """Work around a LoopStructural core bug that crops the boundary above a domain fault. @@ -1064,10 +1113,7 @@ def update_foliation_features(self): data=data, force_constrained=True, **self._extra_kwargs(groupname), - nelements=PlgSettingsStructure.interpolator_nelements, - npw=PlgSettingsStructure.interpolator_npw, - cpw=PlgSettingsStructure.interpolator_cpw, - regularisation=PlgSettingsStructure.interpolator_regularisation, + **self.interpolator_arguments(groupname, data), ) fault_name, flipped = self._base_fault_boundary(group) if fault_name is None or not self._build_domain_fault_boundary( @@ -1413,10 +1459,7 @@ def update_fault_features(self): fault_pitch=pitch, data=data, faults=self._cutting_faults_for(fault_name), - nelements=PlgSettingsStructure.interpolator_nelements, - npw=PlgSettingsStructure.interpolator_npw, - cpw=PlgSettingsStructure.interpolator_cpw, - regularisation=PlgSettingsStructure.interpolator_regularisation, + **self.interpolator_arguments(fault_name, data), ) self._build_parametric_faults() self.apply_fault_abutting_relationships() @@ -1458,13 +1501,11 @@ def _build_parametric_faults(self): for name, spec in self.parametric_faults.items(): self._report_progress(f"Building fault '{name}'") try: + fault_data = parametric_fault.frame_data(spec) self.model.create_and_add_fault( name, - data=parametric_fault.frame_data(spec), - nelements=PlgSettingsStructure.interpolator_nelements, - npw=PlgSettingsStructure.interpolator_npw, - cpw=PlgSettingsStructure.interpolator_cpw, - regularisation=PlgSettingsStructure.interpolator_regularisation, + data=fault_data, + **self.interpolator_arguments(name, fault_data), **parametric_fault.fault_arguments(spec), ) except Exception as e: @@ -2105,7 +2146,9 @@ def _create_foliation( ): """Create the foliation feature in the model; see `add_foliation`.""" data, kwargs = self._foliation_data(name, data, sampler, use_z_coordinate) - foliation = self.model.create_and_add_foliation(name, data=data, **kwargs) + foliation = self.model.create_and_add_foliation( + name, data=data, **kwargs, **self.interpolator_arguments(name, data) + ) if not restrict_to_stratigraphic_domain: foliation.regions = [ r for r in foliation.regions if not isinstance(r, UnconformityFeature) diff --git a/loopstructural/toolbelt/preferences.py b/loopstructural/toolbelt/preferences.py index 083ce8b..158f453 100644 --- a/loopstructural/toolbelt/preferences.py +++ b/loopstructural/toolbelt/preferences.py @@ -32,6 +32,7 @@ class PlgSettingsStructure: version: str = __version__ interpolator_type: str = 'FDI' interpolator_nelements: int = 50000 + interpolator_nelements_auto: bool = True interpolator_regularisation: float = 1.0 interpolator_cpw: float = 1.0 interpolator_npw: float = 1.0 diff --git a/tests/unit/test_interpolation_size.py b/tests/unit/test_interpolation_size.py new file mode 100644 index 0000000..4eb9f8b --- /dev/null +++ b/tests/unit/test_interpolation_size.py @@ -0,0 +1,107 @@ +"""Pytest tests for the choice of the number of interpolation elements. + +The module does not import QGIS, so the tests run in the fast tests/unit/ job. +""" + +import numpy as np +import pandas as pd + +from loopstructural.main.interpolation_size import ( + MAX_NELEMENTS, + MIN_NELEMENTS, + DataSummary, + orientation_spread, + suggest_nelements, + summarise_data, +) + + +class TestSuggestNelements: + def test_no_data_gives_the_minimum(self): + assert suggest_nelements(DataSummary()) == MIN_NELEMENTS + + def test_more_data_gives_more_or_equal_elements(self): + previous = 0 + for n in (1, 10, 100, 1_000, 10_000, 100_000): + value = suggest_nelements(DataSummary(n_value=n, n_orientation=n, n_surfaces=3)) + assert value >= previous + previous = value + + def test_limits(self): + assert suggest_nelements(DataSummary(n_value=1)) == MIN_NELEMENTS + assert suggest_nelements(DataSummary(n_value=10**7)) == MAX_NELEMENTS + + def test_result_is_a_multiple_of_1000(self): + assert suggest_nelements(DataSummary(n_value=777, n_orientation=33)) % 1000 == 0 + + def test_spread_gives_more_elements(self): + parallel = suggest_nelements(DataSummary(n_value=200, n_orientation=200, spread=0)) + spread = suggest_nelements(DataSummary(n_value=200, n_orientation=200, spread=1)) + assert spread > parallel + + def test_surfaces_give_more_elements_and_stop_at_ten(self): + one = suggest_nelements(DataSummary(n_value=200, n_surfaces=1)) + five = suggest_nelements(DataSummary(n_value=200, n_surfaces=5)) + ten = suggest_nelements(DataSummary(n_value=200, n_surfaces=10)) + twenty = suggest_nelements(DataSummary(n_value=200, n_surfaces=20)) + assert one < five < ten + assert ten == twenty + + def test_small_and_large_data_differ(self): + small = summarise_data(_values(20)) + large = summarise_data(_values(5_000)) + assert suggest_nelements(small) < suggest_nelements(large) + + +def _values(n): + return pd.DataFrame({'X': range(n), 'val': [0.0] * n, 'feature_name': 'a'}) + + +class TestOrientationSpread: + def test_parallel_planes_have_no_spread(self): + assert orientation_spread(np.tile([0, 0, 1], (10, 1))) == 0 + + def test_all_directions_have_spread_near_one(self): + rng = np.random.default_rng(0) + normals = rng.normal(size=(5_000, 3)) + assert orientation_spread(normals) > 0.95 + + def test_sign_does_not_matter(self): + rng = np.random.default_rng(1) + normals = rng.normal(size=(50, 3)) + flipped = normals.copy() + flipped[::2] *= -1 + assert np.isclose(orientation_spread(normals), orientation_spread(flipped)) + + def test_one_orientation_has_no_spread(self): + assert orientation_spread([[1, 0, 0]]) == 0 + + +class TestSummariseData: + def test_none_and_empty(self): + assert summarise_data(None) == DataSummary() + assert summarise_data(pd.DataFrame()) == DataSummary() + + def test_counts_and_surfaces(self): + df = pd.DataFrame( + { + 'val': [0.0, 0.0, 10.0, np.nan], + 'nx': [np.nan, np.nan, np.nan, 0.0], + 'ny': [np.nan, np.nan, np.nan, 0.0], + 'nz': [np.nan, np.nan, np.nan, 1.0], + } + ) + summary = summarise_data(df) + assert summary.n_value == 3 + assert summary.n_orientation == 1 + assert summary.n_surfaces == 2 + + def test_strike_dip_spread(self): + parallel = pd.DataFrame({'strike': [10.0] * 5, 'dip': [30.0] * 5}) + mixed = pd.DataFrame({'strike': [0.0, 90.0, 180.0, 270.0], 'dip': [60.0] * 4}) + assert summarise_data(parallel).spread == 0 + assert summarise_data(mixed).spread > 0.3 + + def test_opposite_dip_directions_spread(self): + a = summarise_data(pd.DataFrame({'strike': [0.0, 180.0], 'dip': [30.0, 30.0]})) + assert a.spread > 0