diff --git a/README.md b/README.md index 4af2a0f..21c6ef3 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ [![PyPI license](https://img.shields.io/pypi/l/attpc_spyral.svg)](https://pypi.python.org/pypi/attpc_spyral/) [![DOI](https://zenodo.org/badge/528950398.svg)](https://doi.org/10.5281/zenodo.14143006) -Spyral is an analysis library for data from the Active Target Time Projection Chamber (AT-TPC). Spyral provides a flexible analysis pipeline, transforming the raw trace data into physical observables over several tunable steps. The analysis pipeline is also extensible, supporting a diverse array of datasets. Sypral can process multiple data files in parallel, allowing for scalable performance over larger experiment datasets. +Spyral is an analysis library for data from the Active Target Time Projection Chamber (AT-TPC). Spyral provides a flexible analysis pipeline, transforming the raw trace data into physical observables over several tunable steps. The analysis pipeline is also extensible, supporting a diverse array of datasets. Spyral can process multiple data files in parallel, allowing for scalable performance over larger experiment datasets. ## Installation @@ -51,6 +51,7 @@ from spyral import ( DetectorParameters, ClusterParameters, OverlapJoinParameters, + TripclustParameters, SolverParameters, EstimateParameters, DEFAULT_MAP, @@ -63,7 +64,7 @@ workspace_path = Path("/some/workspace/path/") trace_path = Path("/some/trace/path/") run_min = 94 -run_max = 94 +run_max = 97 n_processes = 4 pad_params = PadParameters( @@ -104,16 +105,42 @@ det_params = DetectorParameters( cluster_params = ClusterParameters( min_cloud_size=50, - min_points=3, - min_size_scale_factor=0.05, - min_size_lower_cutoff=10, - cluster_selection_epsilon=10.0, - overlap_join=OverlapJoinParameters( - min_cluster_size_join=15.0, - circle_overlap_ratio=0.25, + hdbscan_parameters = None, + # hdbscan_parameters = HdbscanParameters( + # min_points=3, + # min_size_scale_factor=0.03, + # min_size_lower_cutoff=10, + # cluster_selection_epsilon=10.0), + # overlap_join=OverlapJoinParameters( + # min_cluster_size_join=15, + # circle_overlap_ratio=0.25, + # ), + # continuity_join=None, + continuity_join = ContinuityJoinParameters( + join_radius_fraction=0.4, + join_z_fraction=0.2), + overlap_join=None, + outlier_scale_factor=0.1, + direction_threshold=0.5, + # tripclust_parameters=None, + tripclust_parameters=TripclustParameters( + r=6, + rdnn=True, + k=12, + n=3, + a=0.03, + s=0.3, + sdnn=True, + t=0.0, + tauto=True, + dmax=0.0, + dmax_dnn=False, + ordered=True, + link=0, + m=50, + postprocess=False, + min_depth=25, ), - continuity_join=None, - outlier_scale_factor=0.05, ) estimate_params = EstimateParameters( @@ -172,7 +199,7 @@ The core of Spyral is the Pipeline. A Pipeline in a complete description of an a ### Parallel Processing -Spyral is capable of running multiple data files in parallel. This is acheived through the python `multiprocessing` library. In the `start_pipeline` function a parameter named `n_processors` indicates to Spyral the *maximum* number of processors which can be spawned. Spyral will then inspect the data load that was submitted in the configuration and attempt to balance the load across the processors as equally as possible. +Spyral is capable of running multiple data files in parallel. This is achieved through the python `multiprocessing` library. In the `start_pipeline` function a parameter named `n_processors` indicates to Spyral the *maximum* number of processors which can be spawned. Spyral will then inspect the data load that was submitted in the configuration and attempt to balance the load across the processors as equally as possible. Some notes about parallel processing: @@ -183,7 +210,7 @@ Some notes about parallel processing: ### Logs and Output -Spyral creates a set of logfiles when it is run (located in the log directory of the workspace). These logfiles can contain critical information describing the state of Spyral. In particular, if Spyral has a crash, the logfiles can be useful for determining what went wrong. A logfile is created for each process (including the parent process). The files are labeld by process number (or as parent in the case of the parent). +Spyral creates a set of logfiles when it is run (located in the log directory of the workspace). These logfiles can contain critical information describing the state of Spyral. In particular, if Spyral has a crash, the logfiles can be useful for determining what went wrong. A logfile is created for each process (including the parent process). The files are labeled by process number (or as parent in the case of the parent). ## Notebooks diff --git a/docs/user_guide/config/cluster.md b/docs/user_guide/config/cluster.md index f1fecc8..9f6f442 100644 --- a/docs/user_guide/config/cluster.md +++ b/docs/user_guide/config/cluster.md @@ -2,38 +2,48 @@ The Cluster parameters which control the clustering, joining, and outlier detection algorithms. -The default recommended settings when using the continuity join method are: +The default recommended settings for HDBSCAN when using the continuity join method are: ```python cluster_params = ClusterParameters( min_cloud_size=50, - min_points=5, - min_size_scale_factor=0.0, - min_size_lower_cutoff=5, - cluster_selection_epsilon=13.0, + # hdbscan_parameters = None + hdbscan_parameters = HdbscanParameters( + min_points=5, + min_size_scale_factor=0.0, + min_size_lower_cutoff=5, + cluster_selection_epsilon=13.0), overlap_join=None, continuity_join=ContinuityOverlapParamters( join_radius_fraction=0.3, join_z_fraction=0.2, ), + # continuity_join = None outlier_scale_factor=0.05, + direction_threshold = 0.5, + tripclust_parameters = None, ) ``` -The default recommended settings when using the circle overlap join method are: +The default recommended settings for HDBSCAN when using the circle overlap join method are: ```python cluster_params = ClusterParameters( min_cloud_size=50, - min_points=3, - min_size_scale_factor=0.05, - min_size_lower_cutoff=10, - cluster_selection_epsilon=10.0, + # hdbscan_parameters = None + hdbscan_parameters = HdbscanParameters( + min_points=3, + min_size_scale_factor=0.05, + min_size_lower_cutoff=10, + cluster_selection_epsilon=10.0), + # overlap_join=None, overlap_join=OverlapJoinParameters( min_cluster_size_join=15.0, circle_overlap_ratio=0.25, ), continuity_join=None, + direction_threshold = 0.5, + tripclust_parameters = None, outlier_scale_factor=0.05, ) ``` @@ -46,56 +56,148 @@ This is the minimum size a point cloud must be (in number of points) to be sent ## minimum_points -The minimum number of samples (points) in a neighborhood for a point to be a core point. This is a re-exposure of the `min_samples` parameter of -[scikit-learn's HDBSCAN](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.HDBSCAN.html#sklearn.cluster.HDBSCAN). See -their documentation for more details. Larger values will make the algorithm more likely to identify points as noise. See the original +The minimum number of samples (points) in a neighborhood for a point to be a core point. This is a re-exposure of the `min_samples` parameter of +[scikit-learn's HDBSCAN](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.HDBSCAN.html#sklearn.cluster.HDBSCAN). See +their documentation for more details. Larger values will make the algorithm more likely to identify points as noise. See the original [HDBSCAN docs](https://hdbscan.readthedocs.io/en/latest/parameter_selection.html#) for details on why this parameter is important and how it can impact the data. ## minimum_size_scale_factor -HDBSCAN requires a minimum size (the hyper parameter `min_cluster_size` in +HDBSCAN requires a minimum size (the hyper parameter `min_cluster_size` in [scikit-learn's HDBSCAN](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.HDBSCAN.html#sklearn.cluster.HDBSCAN)) in terms of samples for a group to -be considered a valid cluster. AT-TPC point clouds vary dramatically in size, from tens of points to thousands. To handle this wide scale, we use a scale factor -to determine the appropriate minimum size, where `min_cluster_size = minimum_size_scale_factor * n_cloud_points`. The default value was found through some testing, +be considered a valid cluster. AT-TPC point clouds vary dramatically in size, from tens of points to thousands. To handle this wide scale, we use a scale factor +to determine the appropriate minimum size, where `min_cluster_size = minimum_size_scale_factor * n_cloud_points`. The default value was found through some testing, and may need serious adjustment to produce best results. Note that the scale factor should be *small*. Update: in the case of continuity joining, scaling was not needed. ## minimum_size_lower_cutoff -As discussed in the above `minimum_size_scale_factor`, we need to scale the `min_cluster_size` parameter to the size of the point cloud. However, there must be -a lower limit (i.e. you can't have a minimum cluster size of 0). This parameter sets the lower limit; that is any `min_cluster_size` calculated using the scale factor -that is smaller than this cutoff is replaced with the cutoff value. As an example, if the cutoff is set to 10 and the calculated value is 50, the calculated value would +As discussed in the above `minimum_size_scale_factor`, we need to scale the `min_cluster_size` parameter to the size of the point cloud. However, there must be +a lower limit (i.e. you can't have a minimum cluster size of 0). This parameter sets the lower limit; that is any `min_cluster_size` calculated using the scale factor +that is smaller than this cutoff is replaced with the cutoff value. As an example, if the cutoff is set to 10 and the calculated value is 50, the calculated value would be used. However, if the calculated value is 5, the cutoff would be used instead. ## cluster_selection_epsilon -A re-exposure of the `cluster_selection_epsilon` paramter of -[scikit-learn's HDBSCAN](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.HDBSCAN.html#sklearn.cluster.HDBSCAN). This parameter will merge clusters that -are less than epsilon apart. Note that this epsilon must be on the scale of the scaled data (i.e. it is not in normal units). The impact of this parameter is large, and -small changes to this value can produce dramatically different results. Larger values will bias the clustering to assume the point cloud is onesingle cluster (or all noise), -while smaller values will cause the algorithm to revert to the default result of HDBSCAN. See the original +A re-exposure of the `cluster_selection_epsilon` parameter of +[scikit-learn's HDBSCAN](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.HDBSCAN.html#sklearn.cluster.HDBSCAN). This parameter will merge clusters that +are less than epsilon apart. Note that this epsilon must be on the scale of the scaled data (i.e. it is not in normal units). The impact of this parameter is large, and +small changes to this value can produce dramatically different results. Larger values will bias the clustering to assume the point cloud is one single cluster (or all noise), +while smaller values will cause the algorithm to revert to the default result of HDBSCAN. See the original [HDBSCAN docs](https://hdbscan.readthedocs.io/en/latest/parameter_selection.html#) for details on why this parameter is important and how it can impact the data. +The recommended parameters when using the TRIPCLUST (Dalitz) clustering algorithm are: +(see also the publication [Dalitz clustering](https://doi.org/10.1016/j.cpc.2018.09.010) for more information) + +```python +# tripclust_parameters=None, +tripclust_parameters=TripclustParameters( + r=6, + rdnn=True, + k=12, + n=3, + a=0.03, + s=0.3, + sdnn=True, + t=0.0, + tauto=True, + dmax=0.0, + dmax_dnn=False, + ordered=True, + link=0, + m=50, + postprocess=False, + min_depth=25, +``` +Breakdown of each of the parameters: + +## r + +The neighbor distance for smoothing. + +## rdnn + +When this boolean is set to True the value of r is calculated automatically using dnn. + +## k + +The number of tested neighbors for each triplet mid point. + +## n + +The maximum number of triplets retained for each mid point. + +## a + +The maximum value of angle between the two triplet branches, expressed as 1-cos(alpha). + +## s + +The distance scale factor in the metric of triplet distance. + +## sdnn + +When set to True, the value of s is calculated automatically using dnn. + +## t + +The threshold for the "distance" between triplets. + +## tauto + +When set to True, the value of t is set automatically. + +## dmax + +The maximum gap width. + +## dmax_dnn + +Use dnn for dmax. + +## ordered + +When True the point cloud is ordered in chronological order. + +## link + +Linkage method for hierarchical clustering of the triplets. + +## m + +The minimum number of triplets per cluster. + +## postprocess + +When set to True, the post process algorithm is attempted. + +## min_depth + +Value of the minimum depth used in the post processing algorithm. + + + ## outlier_scale_factor -We use [scikit-learn's LocalOutlierFactor](https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.LocalOutlierFactor.html) as a last round of noise elimination on -a cluster-by-cluster basis. This algorithim requires a number of neighbors to search over (the `n_neighbors` parameter). As with the `min_cluster_size` in HDBSCAN, we need to -scale this value off the size of the cluster. This factor multiplied by the size of the cluster gives the number of neighbors to search over -(`n_neighbors = outlier_scale_factor * cluster_size`). This value tends to have a "sweet spot" where it is most effective. If it is too large, every point has basically the -same outlier factor as you're including the entire cluster for every point. If it is too small the variance between neighbors can be too large and the results will be +We use [scikit-learn's LocalOutlierFactor](https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.LocalOutlierFactor.html) as a last round of noise elimination on +a cluster-by-cluster basis. This algorithm requires a number of neighbors to search over (the `n_neighbors` parameter). As with the `min_cluster_size` in HDBSCAN, we need to +scale this value off the size of the cluster. This factor multiplied by the size of the cluster gives the number of neighbors to search over +(`n_neighbors = outlier_scale_factor * cluster_size`). This value tends to have a "sweet spot" where it is most effective. If it is too large, every point has basically the +same outlier factor as you're including the entire cluster for every point. If it is too small the variance between neighbors can be too large and the results will be unpredictable. Note that if the value of `outlier_scale_factor * cluster_size` is less than 2, `n_neighbors` will be set to 2 as this is the minimum allowed value. +This algorithm is also used to clean clusters prior to the joining process, using the same parameter. ## Overlap Join Parameters ### min_cluster_size_join -The minimum size of a cluster for it to be considered in the joining step of the clustering. After HDBSCAN has made the initial clusters we attempt to combine any clusters which +The minimum size of a cluster for it to be considered in the joining step of the clustering. After HDBSCAN has made the initial clusters we attempt to combine any clusters which have overlapping circles in the 2-D projection (see `circle_overlap_ratio`). However, many times, small pockets of noise will be clustered and often sit within the larger trajectory. To avoid these being joined we require a cluster to have a minimum size. ### circle_overlap_ratio -The minimum amount of overlap between circles fit to two clusters for the clusters to be joined together into a single cluster. HDBSCAN often fractures trajectories into multiple +The minimum amount of overlap between circles fit to two clusters for the clusters to be joined together into a single cluster. HDBSCAN often fractures trajectories into multiple clusters as the point density changes due to the pad size, gaps, etc. These fragments are grouped together based on how much circles fit on their 2-D projections (X-Y) overlap. ## Continuity Join Parameters diff --git a/pyproject.toml b/pyproject.toml index 536387c..7c9aaa2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,13 +1,14 @@ [project] name = "attpc_spyral" -version = "1.0.0" +version = "1.1.1" description = "AT-TPC analysis pipeline" authors = [ {name = "gwm17", email = "gordonmccann215@gmail.com"}, {name = "turinath", email = "turi@frib.msu.edu"}, + {name = "DBazin", email = "bazin@frib.msu.edu"}, ] dependencies = [ - "spyral-utils>=2.0.0", + "spyral-utils>=2.1.0", "contourpy>=1.2.1", "h5py>=3.11.0", "lmfit>=1.3.0", diff --git a/src/spyral/__init__.py b/src/spyral/__init__.py index 49c136d..17bf30d 100644 --- a/src/spyral/__init__.py +++ b/src/spyral/__init__.py @@ -20,6 +20,7 @@ ClusterParameters, OverlapJoinParameters, ContinuityJoinParameters, + TripclustParameters, EstimateParameters, SolverParameters, DEFAULT_MAP, @@ -52,6 +53,7 @@ "ClusterParameters", "OverlapJoinParameters", "ContinuityJoinParameters", + "TripclustParameters", "EstimateParameters", "SolverParameters", "DEFAULT_MAP", diff --git a/src/spyral/core/cluster.py b/src/spyral/core/cluster.py index 4f5fcea..7603427 100644 --- a/src/spyral/core/cluster.py +++ b/src/spyral/core/cluster.py @@ -4,9 +4,27 @@ from dataclasses import dataclass from sklearn.neighbors import LocalOutlierFactor from scipy.interpolate import BSpline, make_smoothing_spline +from enum import Enum EMPTY_DATA = np.empty(0, dtype=float) +class Direction(Enum): + """Enum for the direction of a trajectory + + Attributes + ---------- + NONE: int + Invalid value (-1) + FORWARD: int + Trajectory traveling in the positive z-direction (0) + BACKWARD: int + Trajectory traveling in the negative z-direction (1) + """ + + NONE = -1 # type: ignore + FORWARD = 0 # type: ignore + BACKWARD = 1 # type: ignore + @dataclass class LabeledCloud: @@ -23,6 +41,7 @@ class LabeledCloud: """ label: int # default is noise label + direction: Direction # default is None point_cloud: PointCloud parent_indicies: np.ndarray @@ -45,6 +64,8 @@ class Cluster: The event number label: int The cluster label from the algorithm + direction: Direction + The direction of the cluster data: ndarray The point cloud data (trimmed down). Contains position, integrated charge x_spline: BSpline | None @@ -68,10 +89,12 @@ def __init__( self, event: int = -1, label: int = -1, + direction: Direction = Direction.NONE, data: np.ndarray = EMPTY_DATA, ): self.event = event self.label = label + self.direction = direction self.data = data self.x_spline: BSpline | None = None self.y_spline: BSpline | None = None @@ -175,6 +198,6 @@ def convert_labeled_to_cluster( data[:, :3] = cloud.point_cloud.data[:, :3] # position data[:, 3] = cloud.point_cloud.data[:, 4] # peak integral data[:, 4] = cloud.point_cloud.data[:, 7] # scale (big or small) - cluster = Cluster(cloud.point_cloud.event_number, cloud.label, data) + cluster = Cluster(cloud.point_cloud.event_number, cloud.label, cloud.direction, data) outliers = cluster.drop_outliers(params.outlier_scale_factor) return (cluster, outliers) diff --git a/src/spyral/core/clusterize.py b/src/spyral/core/clusterize.py index e705499..d83c34f 100644 --- a/src/spyral/core/clusterize.py +++ b/src/spyral/core/clusterize.py @@ -1,21 +1,69 @@ from .point_cloud import PointCloud -from .cluster import LabeledCloud, Cluster, convert_labeled_to_cluster +from .cluster import LabeledCloud, Cluster, convert_labeled_to_cluster, Direction from .config import ClusterParameters from ..geometry.circle import least_squares_circle import sklearn.cluster as skcluster import numpy as np +import tripclust +from enum import Enum +from scipy.interpolate import BSpline, make_smoothing_spline +from sklearn.neighbors import LocalOutlierFactor +from sklearn.ensemble import IsolationForest -NOISE_LABEL: int = -1 +NOISE_LABEL: int = -1 class ClusterError(Exception): pass +# Credit Daniela Ramirez-Chavez 2025 ###### +def get_direction(cluster: LabeledCloud, params: ClusterParameters): + """Determine the direction of a cluster based on the + differences in angle between adjacent points in z. + Smoothing is applied to iron out fluctuations and + a statistical determination is based on the fraction + of angles greater than a chosen threshold + + Parameters + ---------- + cluster: LabeledCloud + the cluster in consideration + params: ClusterParameters + contains the parameters related to clustering + """ + cluster_data = cluster.point_cloud.data.copy() + if len(cluster_data) > 5: + x_spline = make_smoothing_spline(cluster_data[:, 2], cluster_data[:, 0], lam=10.0) + y_spline = make_smoothing_spline(cluster_data[:, 2], cluster_data[:, 1], lam=10.0) + cluster_data[:, 0] = x_spline(cluster_data[:, 2]) # type: ignore + cluster_data[:, 1] = y_spline(cluster_data[:, 2]) # type: ignore + cx, cy, _, _ = least_squares_circle(cluster_data[:,0], cluster_data[:,1]) + + angle_right = np.arctan2(cluster_data[:,1] - cy, cluster_data[:,0] - cx) + + for i in range(0,len(angle_right)): + if angle_right[i]<0: + angle_right[i]= 2*np.pi + angle_right[i] + + da_right = np.diff(np.rad2deg(angle_right)) + + positive_fraction = sum(1 for num in da_right if num > 0) / len(da_right) + negative_fraction = 1 - positive_fraction + if positive_fraction > params.direction_threshold: + cluster.direction = Direction.FORWARD + elif negative_fraction > params.direction_threshold: + cluster.direction = Direction.BACKWARD + else: + cluster.direction = Direction.NONE + def join_clusters( clusters: list[LabeledCloud], params: ClusterParameters, labels: np.ndarray ) -> tuple[list[LabeledCloud], np.ndarray]: + # First determine direction for each cluster + for cluster in clusters: + get_direction(cluster, params) """Detect which joining algorithm to use and dispatch""" if params.continuity_join is not None: return join_clusters_continuity(clusters, params, labels) @@ -175,7 +223,7 @@ def join_clusters_overlap_step( smaller_area = min(area, comp_area) if ( area_overlap > params.overlap_join.circle_overlap_ratio * smaller_area # type: ignore - ) and (cidx not in groups_index[cluster.label]): + ) and (cidx not in groups_index[cluster.label]) and (cluster.direction == comp_cluster.direction): comp_indicies = groups_index.pop(comp_cluster.label) comp_labels = groups_label.pop(comp_cluster.label) for subs in comp_indicies: @@ -240,12 +288,12 @@ def join_clusters_continuity( after = 0 while before != after and len(jclusters) > 1: before = len(jclusters) - # Order clusters by start time, scaled by mean charge, pad size + # Order clusters by start time jclusters = sorted( jclusters, - key=lambda x: x.point_cloud.data[0, 2] - * np.mean(x.point_cloud.data[:, 3]) - * np.mean(x.point_cloud.data[:, 7]), + key=lambda x: x.point_cloud.data[0, 2], + # * np.mean(x.point_cloud.data[:, 3]) + # * np.mean(x.point_cloud.data[:, 7]), ) jclusters, labels = join_clusters_continuity_step(jclusters, params, labels) after = len(jclusters) @@ -298,59 +346,38 @@ def join_clusters_continuity_step( if cluster.label == NOISE_LABEL or used_clusters[idx]: continue for cidx, comp_cluster in enumerate(clusters): - if comp_cluster.label == NOISE_LABEL or cidx == idx or used_clusters[cidx]: + # if clusters have different directions or already used, try the next + if comp_cluster.label == NOISE_LABEL or cidx == idx or used_clusters[cidx] or comp_cluster.direction != cluster.direction : continue # downstream - avg_pos = np.median(cluster.point_cloud.data[-3:, :3], axis=0) # downstream - rhos = np.linalg.norm(cluster.point_cloud.data[:, :2], axis=1) + # determine center of cluster from circle fit and take positions in that circle ref + centered = cluster.point_cloud.data.copy() + cx, cy, radius, _ = least_squares_circle(cluster.point_cloud.data[:,0], cluster.point_cloud.data[:,1]) + centered[:, 0] -= cx + centered[:, 1] -= cy + avg_pos = np.median(centered[-3:, :3], axis=0) # downstream min_z = np.min(cluster.point_cloud.data[:, 2]) max_z = np.max(cluster.point_cloud.data[:, 2]) - min_rho = np.min(rhos) - max_rho = np.max(rhos) avg_rho = np.linalg.norm(avg_pos[:2]) - avg_phi = np.arctan2(avg_pos[1], avg_pos[0]) - if avg_phi < 0.0: - avg_phi += 2.0 * np.pi # upstream - avg_pos_comp = np.median( - comp_cluster.point_cloud.data[:3, :3], axis=0 - ) # upstream - rhos_comp = np.linalg.norm(comp_cluster.point_cloud.data[:, :2], axis=1) + # determine center of cluster from circle fit and take positions in that circle ref + centered = comp_cluster.point_cloud.data.copy() + cx, cy, comp_radius, _ = least_squares_circle(comp_cluster.point_cloud.data[:,0], comp_cluster.point_cloud.data[:,1]) + centered[:, 0] -= cx + centered[:, 1] -= cy + avg_pos_comp = np.median(centered[:3, :3], axis=0) # upstream min_z_comp = np.min(comp_cluster.point_cloud.data[:, 2]) max_z_comp = np.max(comp_cluster.point_cloud.data[:, 2]) - min_rho_comp = np.min(rhos_comp) - max_rho_comp = np.max(rhos_comp) avg_rho_comp = np.linalg.norm(avg_pos_comp[:2]) - avg_phi_comp = np.arctan2(avg_pos_comp[1], avg_pos_comp[0]) - if avg_phi_comp < 0.0: - avg_phi_comp += 2.0 * np.pi - # compare - is_near_beam_region = ( - np.fabs(avg_rho - 25.0) < 10.0 and np.fabs(avg_rho_comp - 25.0) < 10.0 - ) - z_thresh = ( - max(max_z, max_z_comp) - min(min_z, min_z_comp) - ) * params.continuity_join.join_z_fraction # type: ignore - rho_thresh = ( - max(max_rho, max_rho_comp) - min(min_rho, min_rho_comp) - ) * params.continuity_join.join_radius_fraction # type: ignore + # compare z and rho differences to thresholds + z_thresh = (max(max_z, max_z_comp) - min(min_z, min_z_comp)) * params.continuity_join.join_z_fraction # type: ignore + rho_thresh = (radius + comp_radius) / 2.0 * params.continuity_join.join_radius_fraction z_diff = np.fabs(avg_pos[2] - avg_pos_comp[2]) rho_diff = np.fabs(avg_rho - avg_rho_comp) - phi_diff = np.fabs(avg_phi - avg_phi_comp) - if phi_diff > np.pi: - phi_diff = 2.0 * np.pi - phi_diff - - if ( - (is_near_beam_region and z_diff < z_thresh) - or ( - rho_diff < rho_thresh - and z_diff < z_thresh - and phi_diff < np.pi * 0.25 - ) - ) and (avg_pos[2] < avg_pos_comp[2] or z_diff < 10.0): + if (rho_diff < rho_thresh and z_diff < z_thresh ): comp_indicies = groups_index.pop(comp_cluster.label) comp_labels = groups_label.pop(comp_cluster.label) for subs in comp_indicies: @@ -371,9 +398,10 @@ def join_clusters_continuity_step( continue new_cluster = LabeledCloud( - g, PointCloud(event_number, np.zeros((0, 8))), np.empty(0) + g, Direction.NONE, PointCloud(event_number, np.zeros((0, 8))), np.empty(0) ) for idx in groups_index[g]: + new_cluster.direction = clusters[idx].direction new_cluster.point_cloud.data = np.concatenate( (new_cluster.point_cloud.data, clusters[idx].point_cloud.data), axis=0 ) @@ -436,8 +464,107 @@ def cleanup_clusters( return (cleaned, labels) +def clean_clusters(clusters: list[LabeledCloud], params: ClusterParameters, labels: np.ndarray + ) -> tuple[list[LabeledCloud], np.ndarray]: + """Attempt to clean the clusters using the LocalOutlierFactor algorithm + + Parameters + ---------- + clusters: list[LabeledCloud] + clusters to clean + params: ClusterParameters + Configuration parameters controlling the clustering algorithms + labels: numpy.ndarray + The cluster label for each point in the original point cloud + + Returns + ------- + tuple[list[LabeledCloud], numpy.ndarray] + A two element tuple, the first the list of cleaned clusters, + the second being an updated list of labels for each point in + the point cloud + + """ + cleaned = [] + for cluster in clusters: + if cluster.label == NOISE_LABEL or len(cluster.point_cloud) < 2: + continue + cleaned_cluster, outliers = clean_cluster(cluster, params.outlier_scale_factor) + cleaned.append(cleaned_cluster) + if len(outliers) == 0: + continue + orig_indicies = (cluster.parent_indicies[outliers]).astype(dtype=np.int32) + labels[orig_indicies] = NOISE_LABEL + return (cleaned, labels) + +def clean_cluster(cluster: LabeledCloud, scale) -> tuple[np.ndarray, np.ndarray]: + """Attempt to clean an individual cluster using the LocalOutlierFactor algorithm + + Parameters + ---------- + cluster: LabeledCloud + cluster to clean + scale: Float + Scale factor applied to the algorithm + + Returns + ------- + tuple[numpy.ndarray, numpy.ndarray] + A two element tuple, the first the cleaned cluster, + the second being an array of the outliers + + """ + neighbors = int(scale * len(cluster.point_cloud.data)) # 0.05 default + if neighbors < 2: + neighbors = 2 + test_data = cluster.point_cloud.data[:, :3].copy() + # clf = IsolationForest(contamination=0.1, random_state=42) + # clf.fit(test_data) + # result = clf.predict(test_data) + neigh = LocalOutlierFactor(n_neighbors=neighbors) + result = neigh.fit_predict(test_data) + mask = result > 0 + cleaned_cluster = LabeledCloud(cluster.label, cluster.direction, cluster.point_cloud, cluster.parent_indicies) + cleaned_cluster.point_cloud.data = cleaned_cluster.point_cloud.data[mask] # label=-1 is an outlier + cleaned_cluster.parent_indicies = cleaned_cluster.parent_indicies[mask] + outliers = np.flatnonzero(~mask) + return (cleaned_cluster, outliers) + def form_clusters( pc: PointCloud, params: ClusterParameters +) -> tuple[list[LabeledCloud], np.ndarray]: + """Apply clustering algorithm based on the values of the parameters + + Parameters + ---------- + pc: PointCloud + The point cloud to be clustered + params: ClusterParameters + Configuration parameters controlling the clustering algorithms + + Returns + -------- + tuple[list[LabeledCloud], numpy.ndarray] + Two element tuple, the first being a ist of clusters found by the algorithm with labels + the second being an array of length of the point cloud with each element containing + that point's label. + """ + + if len(pc) < params.min_cloud_size: + return ([], np.empty(0)) + + """Detect which clustering algorithm to use and dispatch""" + if params.hdbscan_parameters is not None: + return hdbscan_clusters(pc, params) + elif params.tripclust_parameters is not None: + return tripclust_clusters(pc, params) + else: + raise ClusterError( + "No clustering parameters were given! Clustering cannot proceed." + ) + +def hdbscan_clusters( + pc: PointCloud, params: ClusterParameters ) -> tuple[list[LabeledCloud], np.ndarray]: """Apply the HDBSCAN clustering algorithm to a PointCloud @@ -460,17 +587,14 @@ def form_clusters( -------- tuple[list[LabeledCloud], numpy.ndarray] Two element tuple, the first being a ist of clusters found by the algorithm with labels - the second being an array of length of the point cloud with each element conatining + the second being an array of length of the point cloud with each element containing that point's label. """ - if len(pc) < params.min_cloud_size: - return ([], np.empty(0)) - n_points = len(pc) - min_size = int(params.min_size_scale_factor * n_points) - if min_size < params.min_size_lower_cutoff: - min_size = params.min_size_lower_cutoff + min_size = int(params.hdbscan_parameters.min_size_scale_factor * n_points) + if min_size < params.hdbscan_parameters.min_size_lower_cutoff: + min_size = params.hdbscan_parameters.min_size_lower_cutoff # Use spatial dimensions cluster_data = np.empty(shape=(len(pc), 3)) @@ -481,8 +605,8 @@ def form_clusters( clusterizer = skcluster.HDBSCAN( # type: ignore min_cluster_size=min_size, - min_samples=params.min_points, - cluster_selection_epsilon=params.cluster_selection_epsilon, + min_samples=params.hdbscan_parameters.min_points, + cluster_selection_epsilon=params.hdbscan_parameters.cluster_selection_epsilon, ) fitted_clusters = clusterizer.fit(cluster_data) @@ -495,8 +619,79 @@ def form_clusters( clusters.append( LabeledCloud( label, + Direction.NONE, PointCloud(pc.event_number, pc.data[mask]), np.flatnonzero(mask), ) ) return (clusters, fitted_clusters.labels_) + + +def tripclust_clusters( + pc: PointCloud, params: ClusterParameters +) -> tuple[list[LabeledCloud], np.ndarray]: + """Apply the tripclust clustering algorithm to a PointCloud + + Parameters + ---------- + pc: PointCloud + The point cloud to be clustered + params: ClusterParameters + Configuration parameters controlling the clustering algorithms + + Analyze a point cloud, and group the points into clusters which in principle should correspond to particle trajectories. + This function uses the tripclust (Dalitz) clustering algorithm, which groups points in triplets and does + a hierarchical clustering on them to find trajectories that are collinear. + For details see the publication from Dalitz et al. (Dalitz et al., Comp. Phys. Comm. 235 (2019), 159, + DOI: https://doi.org/10.1016/j.cpc.2018.09.010) + + Returns + -------- + tuple[list[LabeledCloud], numpy.ndarray] + Two element tuple, the first being a ist of clusters found by the algorithm with labels + the second being an array of length of the point cloud with each element conatining + that point's label. + """ + + if len(pc) < params.min_cloud_size: + return ([], np.empty(0)) + + # Create class to pass parameters and calls to C++ code + tc = tripclust.tripclust() + tc.set_r(params.tripclust_parameters.r) + tc.set_rdnn(params.tripclust_parameters.rdnn) + tc.set_k(params.tripclust_parameters.k) + tc.set_n(params.tripclust_parameters.n) + tc.set_a(params.tripclust_parameters.a) + tc.set_s(params.tripclust_parameters.s) + tc.set_sdnn(params.tripclust_parameters.sdnn) + tc.set_t(params.tripclust_parameters.t) + tc.set_tauto(params.tripclust_parameters.tauto) + tc.set_dmax(params.tripclust_parameters.dmax) + tc.set_dmax_dnn(params.tripclust_parameters.dmax_dnn) + tc.set_ordered(params.tripclust_parameters.ordered) + # tc.set_link(params.tripclust_parameters.link) + tc.set_m(params.tripclust_parameters.m) + tc.set_postprocess(params.tripclust_parameters.postprocess) + tc.set_min_depth(params.tripclust_parameters.min_depth) + + # Perform tripclust (Dalitz) clustering + tc.fill_pointcloud(pc.data) + tc.perform_clustering() + labels = tc.get_labels() + + # Below is based on a copy from form_clusters to get same output + # Select out data into clusters + clusters: list[LabeledCloud] = [] + ulabels = np.unique(labels) + for label in ulabels: + mask = labels == label + clusters.append( + LabeledCloud( + label, + Direction.NONE, + PointCloud(pc.event_number, pc.data[mask]), + np.flatnonzero(mask), + ) + ) + return (clusters, labels) diff --git a/src/spyral/core/config.py b/src/spyral/core/config.py index 8deec02..3b20508 100644 --- a/src/spyral/core/config.py +++ b/src/spyral/core/config.py @@ -191,13 +191,69 @@ class OverlapJoinParameters: @dataclass -class ClusterParameters: - """Parameters for clustering, cluster joining, and cluster cleaning +class TripclustParameters: + """Parameters for the tripclust (Dalitz) clustering algorithms + + Attributes + ---------- + r: float + maximum neighbour distance for smoothing (default 2) + rdnn: boolean + whether or not compute r with dnn (default true) + k: int + number of tested neighbours of triplet mid point (default 19) + n: int + max number of triplets to one mid point (default 2) + a: float + 1 - cos alpha, where alpha is the angle between the two triplet branches (default 0.03) + s: float + distance scale factor in metric (default 0.3) + sdnn: boolean + whether or not compute s with dnn (default true) + t: float + threshold for cdist in clustering (default 0.0) + tauto: boolean + whether or not auto generate t (default true) + dmax: float + maximum gap width (default 0.0) + dmax_dnn: boolean + whether or not use dnn for dmax (default false) + ordered: boolean + whether or not points are in chronological order (default false) + link: int + linkage method for clustering (default 0=SINGLE) + m: int + min number of triplets per cluster (default 5) + postprocess: boolean + whether or not post processing should be enabled (default false) + min_depth: int + minimum number of points making a branch in curve in post processing (default 25) + """ + + r: float + rdnn: bool + k: int + n: int + a: float + s: float + sdnn: bool + t: float + tauto: bool + dmax: float + dmax_dnn: bool + ordered: bool + link: int + m: int + postprocess: bool + min_depth: int + + +@dataclass +class HdbscanParameters: + """Parameters for clustering using the HDBSCAN algorithm Attributes ---------- - min_cloud_size: int - The minimum size for a point cloud to be clustered min_points: int min_samples parameter in scikit-learns' HDBSCAN algorithm min_size_scale_factor: int @@ -209,22 +265,43 @@ class ClusterParameters: cluster_selection_epsilon: float cluster_selection_epsilon parameter in scikit-learn's HDBSCAN algorithm. Clusters less than this distance apart are merged in the hierarchy - min_cluster_size_join: int - The minimum size of a cluster for it to be included in the joining algorithm - circle_overlap_ratio: float - minimum overlap ratio between two circles in the cluster joining algorithm - outlier_scale_factor: float - Factor which is multiplied by the number of points in a trajectory to set the number of neighbors parameter - for scikit-learns LocalOutlierFactor test """ - min_cloud_size: int min_points: int min_size_scale_factor: float min_size_lower_cutoff: int cluster_selection_epsilon: float + + +@dataclass +class ClusterParameters: + """Parameters for clustering, cluster joining, and cluster cleaning + + Attributes + ---------- + min_cloud_size: int + The minimum size for a point cloud to be clustered + hdbscan_parameters: HdbscanParameters + Parameters used when selecting the HDBSCAN algorithm + tripclust_parameters: TripclustParameters + Parameters used when selecting the Tripclust (Dalitz) algorithm + overlap_join: OverlapJoinParameters + Parameters used when selecting overlap joining + continuity_join: ContinuityJoinParameters + Parameters used when selecting continuitu joining + direction_threshold: float + Fraction threshold for the determination of the direction for each cluster + outlier_scale_factor: float + Factor which is multiplied by the number of points in a trajectory to set the number of neighbors parameter + for scikit-learns LocalOutlierFactor test + """ + + min_cloud_size: int + hdbscan_parameters: HdbscanParameters | None + tripclust_parameters: TripclustParameters | None overlap_join: OverlapJoinParameters | None continuity_join: ContinuityJoinParameters | None + direction_threshold: float outlier_scale_factor: float diff --git a/src/spyral/core/estimator.py b/src/spyral/core/estimator.py index a566ae0..899b40f 100644 --- a/src/spyral/core/estimator.py +++ b/src/spyral/core/estimator.py @@ -1,4 +1,4 @@ -from .cluster import Cluster +from .cluster import Cluster, Direction from .config import DetectorParameters, EstimateParameters from ..geometry.circle import generate_circle_points, least_squares_circle from .spy_log import spyral_warn @@ -6,28 +6,9 @@ import numpy as np import math from scipy.stats import linregress -from enum import Enum from dataclasses import dataclass -class Direction(Enum): - """Enum for the direction of a trajectory - - Attributes - ---------- - NONE: int - Invalid value (-1) - FORWARD: int - Trajectory traveling in the positive z-direction (0) - BACKWARD: int - Trajectory traveling in the negative z-direction (1) - """ - - NONE = -1 # type: ignore - FORWARD = 0 # type: ignore - BACKWARD = 1 # type: ignore - - @dataclass class EstimateResult: """Container for results of estimation algorithm with metadata @@ -167,7 +148,7 @@ def estimate_physics( return None # Run estimation where we attempt to guess the right direction - result, direction = estimate_physics_pass( + result = estimate_physics_pass( cluster_index, cluster, ic_amplitude, @@ -178,61 +159,9 @@ def estimate_physics( orig_event, detector_params, ) - - # If estimation was consistent or didn't meet valid criteria we're done - if result is not None or (result is None and direction == Direction.NONE): - return result - # If we made a bad guess, try the other direction - elif direction == Direction.FORWARD: - result, _ = estimate_physics_pass( - cluster_index, - cluster, - ic_amplitude, - ic_centroid, - ic_integral, - ic_multiplicity, - orig_run, - orig_event, - detector_params, - Direction.BACKWARD, - ) - else: - result, _ = estimate_physics_pass( - cluster_index, - cluster, - ic_amplitude, - ic_centroid, - ic_integral, - ic_multiplicity, - orig_run, - orig_event, - detector_params, - Direction.FORWARD, - ) return result -def choose_direction(cluster_data: np.ndarray) -> Direction: - """Choose a direction for the trajectory based on which end of the data is in-spiraling - - Parameters - ---------- - cluster_data: np.ndarray - T - - """ - rhos = np.linalg.norm(cluster_data[:, :2], axis=1) # cylindrical coordinates rho - direction = Direction.NONE - - # See if in-spiraling to the window or microgmegas, sets the direction and guess of z-vertex - if rhos[0] < rhos[-1]: - direction = Direction.FORWARD - else: - direction = Direction.BACKWARD - - return direction - - def estimate_physics_pass( cluster_index: int, cluster: Cluster, @@ -243,8 +172,7 @@ def estimate_physics_pass( orig_run: int, orig_event: int, detector_params: DetectorParameters, - chosen_direction: Direction = Direction.NONE, -) -> tuple[EstimateResult | None, Direction]: +) -> EstimateResult: """Estimate the physics parameters for a cluster which could represent a particle trajectory Estimation is an imprecise process (by definition), and as such this algorithm requires a lot of @@ -283,7 +211,9 @@ def estimate_physics_pass( """ - direction = chosen_direction + direction = cluster.direction # We already know the direction from the clustering phase + if direction == Direction.NONE: # We could not determine the direction of the track + return None vertex = np.array([0.0, 0.0, 0.0]) # reaction vertex center = np.array([0.0, 0.0, 0.0]) # spiral center # copy the data so we can modify it without worrying about side-effects @@ -291,8 +221,8 @@ def estimate_physics_pass( # If chosen direction is set to NONE, we want to have the algorithm # try to decide which direction the trajectory is going - if direction == Direction.NONE: - direction = choose_direction(cluster_data) + # if direction == Direction.NONE: + # direction = choose_direction(cluster_data) if direction == Direction.BACKWARD: cluster_data = np.flip(cluster_data, axis=0) @@ -327,23 +257,19 @@ def estimate_physics_pass( # Since we fit to rho_to_vertex, just find intercept point # Check to see if slope is zero, as this can lead to NaN's if fit.slope == 0.0: # type: ignore - return (None, Direction.NONE) + return None vertex[2] = -1.0 * fit.intercept / fit.slope # type: ignore center[2] = vertex[2] # Toss tracks whose verticies are not close to the origin in x,y if vertex_rho > detector_params.beam_region_radius: - return (None, Direction.NONE) + return None polar = math.atan(fit.slope) # type: ignore # We have a self consistency case here. Polar should match chosen Direction - if (polar > 0.0 and direction == Direction.BACKWARD) or ( - polar < 0.0 and direction == Direction.FORWARD - ): - return ( - None, - direction, - ) # Our direction guess was bad, we need to try again with the other direction + if (polar > 0.0 and direction == Direction.BACKWARD) or (polar < 0.0 and direction == Direction.FORWARD): + print("Direction={direction} and polar={polar} not consistent!") + return None # Our direction guess was bad, we need to try again with the other direction elif direction is Direction.BACKWARD: polar += math.pi @@ -372,7 +298,7 @@ def estimate_physics_pass( # arclength += np.linalg.norm(first_arc[idx + 1, :3] - first_arc[idx, :3]) charge_deposited += first_arc[idx + 1, 3] if charge_deposited == first_arc[0, 3]: - return (None, Direction.NONE) + return None # Use the splines to do a fine-grained line integral to calculate the distance points = np.empty((1000, 3)) @@ -384,8 +310,7 @@ def estimate_physics_pass( dEdx = charge_deposited / arclength # fill in our map - return ( - EstimateResult( + return EstimateResult( event=cluster.event, cluster_index=cluster_index, cluster_label=cluster.label, @@ -408,7 +333,5 @@ def estimate_physics_pass( sqrt_dEdx=np.sqrt(np.fabs(dEdx)), dE=charge_deposited, arclength=arclength, - direction=direction.value, - ), - direction, - ) + direction=direction, + ) diff --git a/src/spyral/phases/cluster_phase.py b/src/spyral/phases/cluster_phase.py index 846f4ec..1e084f6 100644 --- a/src/spyral/phases/cluster_phase.py +++ b/src/spyral/phases/cluster_phase.py @@ -3,7 +3,7 @@ from ..core.config import ClusterParameters, DetectorParameters from ..core.status_message import StatusMessage from ..core.point_cloud import PointCloud -from ..core.clusterize import form_clusters, join_clusters, cleanup_clusters +from ..core.clusterize import form_clusters, join_clusters, cleanup_clusters, tripclust_clusters, clean_clusters from ..core.spy_log import spyral_warn, spyral_error, spyral_info from ..core.run_stacks import form_run_string from .schema import POINTCLOUD_SCHEMA, CLUSTER_SCHEMA @@ -145,8 +145,19 @@ def run( # Here we don't need to use the labels array. # We just pass it along as needed. + # Perform clustering using algorithm of choice clusters, labels = form_clusters(cloud, self.cluster_params) - joined, labels = join_clusters(clusters, self.cluster_params, labels) + + # Clean our individual clusters before attempting the joining + clean, labels = clean_clusters(clusters, self.cluster_params, labels) + + # Run joining algorithm if either parameters present + if self.cluster_params.overlap_join is not None or self.cluster_params.continuity_join is not None: + joined, labels = join_clusters(clean, self.cluster_params, labels) + else: + joined = clean + + # Clean up noise points cleaned, _ = cleanup_clusters(joined, self.cluster_params, labels) # Each event can contain many clusters @@ -163,6 +174,7 @@ def run( for cidx, cluster in enumerate(cleaned): local_group = cluster_event_group.create_group(f"cluster_{cidx}") local_group.attrs["label"] = cluster.label + local_group.attrs["direction"] = cluster.direction.value local_group.create_dataset("cloud", data=cluster.data) spyral_info(__name__, f"Phase Cluster complete for run {payload.run_number}") diff --git a/src/spyral/phases/estimation_phase.py b/src/spyral/phases/estimation_phase.py index 2d8611d..8dd15d9 100644 --- a/src/spyral/phases/estimation_phase.py +++ b/src/spyral/phases/estimation_phase.py @@ -2,7 +2,7 @@ from ..core.schema import PhaseResult, ResultSchema from ..core.config import EstimateParameters, DetectorParameters from ..core.status_message import StatusMessage -from ..core.cluster import Cluster +from ..core.cluster import Cluster, Direction from ..core.estimator import estimate_physics from ..core.spy_log import spyral_warn, spyral_error, spyral_info from ..core.run_stacks import form_run_string @@ -148,6 +148,7 @@ def run( cluster = Cluster( idx, local_cluster.attrs["label"], # type: ignore + Direction(local_cluster.attrs["direction"]), # type: ignore local_cluster["cloud"][:].copy(), # type: ignore ) diff --git a/src/spyral/phases/interp_solver_phase.py b/src/spyral/phases/interp_solver_phase.py index d910959..f25c766 100644 --- a/src/spyral/phases/interp_solver_phase.py +++ b/src/spyral/phases/interp_solver_phase.py @@ -285,7 +285,7 @@ def run( # Select the particle group data, beam region of ic, convert to dictionary for row-wise operations estimates_gated = ( estimate_df.filter( - pl.struct([xaxis, yaxis]).map_batches(pid.cut.is_cols_inside) + pl.struct([xaxis, yaxis]).map_batches(pid.cut.is_cols_inside, return_dtype=bool) & (pl.col("ic_amplitude") > self.solver_params.ic_min_val) & (pl.col("ic_amplitude") < self.solver_params.ic_max_val) ) @@ -355,6 +355,7 @@ def run( cluster = Cluster( event, local_cluster.attrs["label"], # type: ignore + local_cluster.attrs["direction"], # type: ignore local_cluster["cloud"][:].copy(), # type: ignore ) orig_run = estimates_gated["orig_run"][row] diff --git a/src/spyral/trace/get_trace.py b/src/spyral/trace/get_trace.py index 6857c8b..0b39125 100644 --- a/src/spyral/trace/get_trace.py +++ b/src/spyral/trace/get_trace.py @@ -156,7 +156,7 @@ def find_peaks( ) for idx, p in enumerate(pks): peak = Peak() - peak.centroid = float(p) + rng.random() + peak.centroid = float(p) + rng.random()/1000.0 peak.amplitude = float(self.trace[p]) peak.positive_inflection = int(np.floor(props["left_ips"][idx])) peak.negative_inflection = int(np.ceil(props["right_ips"][idx])) diff --git a/tests/test_default_pipeline.py b/tests/test_default_pipeline.py index 7364173..d4551ac 100644 --- a/tests/test_default_pipeline.py +++ b/tests/test_default_pipeline.py @@ -57,16 +57,42 @@ ) cluster_params = ClusterParameters( min_cloud_size=50, - min_points=3, - min_size_scale_factor=0.05, - min_size_lower_cutoff=10, - cluster_selection_epsilon=10.0, + hdbscan_parameters = None, + # hdbscan_parameters = HdbscanParameters( + # min_points=3, + # min_size_scale_factor=0.03, + # min_size_lower_cutoff=10, + # cluster_selection_epsilon=10.0), + # overlap_join=OverlapJoinParameters( + # min_cluster_size_join=15, + # circle_overlap_ratio=0.25, + # ), + # continuity_join=None, + continuity_join = ContinuityJoinParameters( + join_radius_fraction=0.4, + join_z_fraction=0.2), overlap_join=None, - continuity_join=ContinuityJoinParameters( - join_radius_fraction=0.3, - join_z_fraction=0.2, + outlier_scale_factor=0.1, + direction_threshold=0.5, + # tripclust_parameters=None, + tripclust_parameters=TripclustParameters( + r=6, + rdnn=True, + k=12, + n=3, + a=0.03, + s=0.3, + sdnn=True, + t=0.0, + tauto=True, + dmax=0.0, + dmax_dnn=False, + ordered=True, + link=0, + m=50, + postprocess=False, + min_depth=25, ), - outlier_scale_factor=0.05, ) estimate_params = EstimateParameters( min_total_trajectory_points=30, smoothing_factor=100.0