(step-detect) feat: add auto-detect for step counters

Author: root-hbx

StarPerf.py

@@ -12,11 +12,11 @@ def main():
     print("Starting StarPerf...")
 
 
-    print("Starting XML Constellations Testing...")
-    # test the core module functionality of various XML constellations
-    import samples.XML_constellation.XML_constellation_test_cases
-    samples.XML_constellation.XML_constellation_test_cases.XML_constellation_test_cases()
-    print("END.")
+    # print("Starting XML Constellations Testing...")
+    # # test the core module functionality of various XML constellations
+    # import samples.XML_constellation.XML_constellation_test_cases
+    # samples.XML_constellation.XML_constellation_test_cases.XML_constellation_test_cases()
+    # print("END.")
 
     print("Starting TLE Constellations Testing")
     # test the core module functionality of various TLE constellations

src/constellation_generation/by_TLE/satellite_to_orbit_mapping.py

@@ -20,29 +20,191 @@
 import src.TLE_constellation.constellation_entity.orbit as ORBIT
 import matplotlib.pyplot as plt
 import numpy as np
+import ruptures as rpt
+from scipy.ndimage import gaussian_filter1d
+from scipy.signal import find_peaks
+from sklearn.cluster import DBSCAN
 
 # Function to automatically detect the number of steps/orbits in RAAN data
 def detect_orbit_number(raans):
     """
     Automatically detect the number of orbits based on RAAN distribution steps
+    Using multiple advanced methods for maximum accuracy
     """
     if len(raans) < 10:
         return 1
 
-    # Calculate differences to find jumps/steps
-    diffs = np.diff(raans)
+    raans = np.array(raans)
+    unique_raans = len(np.unique(raans))
 
-    # Find significant jumps (larger than median + 2*std)
-    diff_threshold = np.median(diffs) + 2 * np.std(diffs)
-    significant_jumps = np.where(diffs > diff_threshold)[0]
+    # If all values are the same, return 1
+    if unique_raans <= 1:
+        return 1
+    
+    # Method 1: Enhanced Ruptures with optimal parameters
+    def ruptures_detection():
+        results = []
+        
+        # Try different models and parameters
+        models = ["rbf", "l1", "l2", "normal"]
+        
+        for model in models:
+            try:
+                # PELT with adaptive penalty
+                algo = rpt.Pelt(model=model, min_size=max(3, len(raans) // 100))
+                algo.fit(raans.reshape(-1, 1))
+                
+                # Try different penalty strategies
+                penalties = [
+                    np.log(len(raans)) * 2,  # BIC-like
+                    len(raans) ** 0.5,       # Square root penalty
+                    len(raans) * 0.01,       # Linear penalty
+                ]
+                
+                for pen in penalties:
+                    try:
+                        changepoints = algo.predict(pen=pen)
+                        if len(changepoints) > 1:  # Exclude trivial case
+                            results.append(len(changepoints))
+                    except Exception:
+                        continue
+                        
+            except Exception:
+                continue
+        
+        return results
+    
+    # Method 2: Histogram-based plateau detection with smoothing
+    def histogram_detection():
+        # Create adaptive number of bins
+        n_bins = max(20, min(len(raans) // 10, 100))
+        hist, bin_edges = np.histogram(raans, bins=n_bins)
+        
+        # Apply Gaussian smoothing
+        smoothed = gaussian_filter1d(hist.astype(float), sigma=1.0)
+        
+        # Find local maxima (peaks)
+        peaks, _ = find_peaks(smoothed, 
+                             height=np.max(smoothed) * 0.05,  # 5% of max height
+                             distance=max(1, len(smoothed) // 50))
+        
+        return len(peaks) if len(peaks) > 0 else 1
+    
+    # Method 3: Gap analysis with statistical significance
+    def gap_analysis():
+        if len(raans) < 20:
+            return 1
+            
+        # Calculate gaps
+        sorted_raans = np.sort(raans)
+        gaps = np.diff(sorted_raans)
+        
+        # Use statistical tests to find significant gaps
+        # Method: Find gaps that are outliers using z-score
+        if len(gaps) > 0:
+            z_scores = np.abs((gaps - np.mean(gaps)) / np.std(gaps))
+            significant_gaps = np.sum(z_scores > 2.0)  # 2 standard deviations
+            return max(1, significant_gaps + 1)
+        
+        return 1
+    
+    # Method 4: Density-based clustering approach
+    def density_clustering():
+        try:            
+            # Reshape for clustering
+            raans_reshaped = raans.reshape(-1, 1)
+            
+            # Use DBSCAN to find dense regions
+            eps = np.std(raans) * 0.1  # Adaptive epsilon
+            db = DBSCAN(eps=eps, min_samples=max(2, len(raans) // 100))
+            labels = db.fit_predict(raans_reshaped)
+            
+            # Count number of clusters (excluding noise: label -1)
+            n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
+            return max(1, n_clusters)
+        except Exception:
+            return 1
+    
+    # Method 5: Jenks optimization with adaptive classes
+    def jenks_optimization():
+        best_gvf = 0
+        best_k = 1
+        
+        max_k = min(len(np.unique(raans)), len(raans) // 5)
+        
+        for k in range(1, max_k + 1):
+            try:
+                breaks = jenkspy.jenks_breaks(raans, n_classes=k)
+                gvf = jenkspy.goodness_of_variance_fit(raans, breaks)
+                
+                # Look for elbow point
+                if k == 1:
+                    best_gvf = gvf
+                    best_k = k
+                elif gvf - best_gvf > 0.03:  # Significant improvement
+                    best_gvf = gvf
+                    best_k = k
+                elif gvf > 0.85 and (gvf - best_gvf) < 0.01:  # Diminishing returns
+                    break
+            except Exception:
+                continue
+                
+        return best_k
+    
+    # Run all methods
+    all_results = []
+    
+    try:
+        # Ruptures (multiple results)
+        ruptures_results = ruptures_detection()
+        all_results.extend(ruptures_results)
+        
+        # Histogram
+        hist_result = histogram_detection()
+        all_results.extend([hist_result] * 2)  # Give it some weight
+        
+        # Gap analysis
+        gap_result = gap_analysis()
+        all_results.append(gap_result)
+        
+        # Density clustering
+        density_result = density_clustering()
+        all_results.append(density_result)
+        
+        # Jenks optimization
+        jenks_result = jenks_optimization()
+        all_results.extend([jenks_result] * 2)  # Give it more weight
+        
+    except Exception as e:
+        print(f"Error in advanced detection: {e}")
+    
+    # Robust combination of results
+    if all_results:
+        # Remove extreme outliers using IQR
+        all_results = np.array(all_results)
+        q75, q25 = np.percentile(all_results, [75, 25])
+        iqr = q75 - q25
+        lower_bound = q25 - 1.5 * iqr
+        upper_bound = q75 + 1.5 * iqr
+        
+        filtered_results = all_results[(all_results >= lower_bound) & (all_results <= upper_bound)]
+        
+        if len(filtered_results) > 0:
+            # Use median for robustness
+            final_result = int(np.median(filtered_results))
+            return max(1, final_result)
 
-    # Number of orbits = number of significant jumps + 1
-    orbits_number = len(significant_jumps) + 1
+    # Final fallback: simple but reliable method
+    sorted_raans = np.sort(raans)
+    gaps = np.diff(sorted_raans)
 
-    # Ensure reasonable bounds (1-20 orbits)
-    orbits_number = max(1, min(orbits_number, 20))
+    if len(gaps) > 0:
+        # Use 95th percentile as threshold
+        threshold = np.percentile(gaps, 95)
+        significant_gaps = np.sum(gaps > threshold)
+        return max(1, significant_gaps + 1)
 
-    return orbits_number
+    return 1
 
 # Parameter :
 # shells : a collection of shell objects that have established corresponding relationships
@@ -57,9 +219,9 @@ def satellite_to_orbit_mapping(shells):
             raans.append(sat.tle_json["RA_OF_ASC_NODE"])
         raans = sorted(raans)
 
-        plt.plot(raans)
-        plt.ylabel('RAANS')
-        plt.show()
+        # plt.plot(raans)
+        # plt.ylabel('RAANS')
+        # plt.show()
 
         # Automatically detect the number of orbits
         orbits_number = detect_orbit_number(raans)