Repository: deepcharles/ruptures (language: python)

Issue:
Mac C library makes KernalCPD non-deterministic, but it is deterministic in Linux container
Repro case:

import ruptures as rpt
import numpy as np

new_list = [-0.0155, 0.0194, 0.0289, 0.0071, -0.0059, -0.0102, 0.0046, 0.0218, 0.0153, 0.0491, 0.016, 0.0365, 0.0388, 0.0516, 0.0222, 0.0019, -0.0418, 0.0, -0.0262, 0.0468, 0.0, 0.0311, 0.0341, -0.0, 0.0569, 0.0206, 0.0336, 0.0615]
trend_error = np.asarray(new_list)

results = set({})
for _ in range(10000):
    change_points = (
        rpt.KernelCPD(kernel="rbf", min_size=7)
        .fit(trend_error)
        .predict(pen=1.0)
    )
    results.add(len(change_points))

print(results)

When running this on a mac I'll get two different values in results : 2, 4
But when running this in a linux container it will consistently yield just 2

When running  rpt.Pelt(model="rbf", jump=1, min_size=7)  I repeatedly get a deterministic result (which matches the result in 2).
 



Discussion / hints:
This is indeed an issue. I am looking into it.
Hi @josh-boehm ,

Thanks for you interest in `ruptures` ! 

@deepcharles told me he can replicate the issue, so he will investigate on this. On my side, I could not replicate the issue. So I am just dumping here the setup infos of the two setups I tested in case it can help in the investigation if this is somehow related to an underlying architecture issue. 

I made your code sample runs on both my MacOS laptop and a linux server and they both returned the same : 
```python
>>> print(results)
{2}
```

Here are the setups :

1. MacOS Catalina

Version 10.15.7
Processor : 2.7 GHz Dual-Core Intel Core i5
Memory : 8 GB 1867 MHz DDR3
Python : 3.6.8
ruptures : uptodate with `main` branch of the [repo](https://github.com/deepcharles/ruptures/commit/3d7ec0eeebd9583a3bed204a199501c1b0fcb953)
numpy : 1.19.2
scipy : 1.5.2

2. Centos

CPU infos :
processor       : 0
vendor_id       : GenuineIntel
cpu family      : 6
model           : 60
model name      : Intel Core Processor (Haswell, no TSX)
stepping        : 1
microcode       : 0x1
cpu MHz         : 2399.996
cache size      : 16384 KB
physical id     : 0
siblings        : 1
core id         : 0
cpu cores       : 1
apicid          : 0
initial apicid  : 0
fpu             : yes
fpu_exception   : yes
cpuid level     : 13
wp              : yes
flags           : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 syscall nx rdtscp lm constant_tsc rep_good nopl xtopology cpuid tsc_known_freq pni pclmulqdq vmx ssse3 fma cx16 pcid sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand hypervisor lahf_lm abm cpuid_fault invpcid_single pti tpr_shadow vnmi flexpriority ept vpid ept_ad fsgsbase bmi1 avx2 smep bmi2 erms invpcid xsaveopt arat
bugs            : cpu_meltdown spectre_v1 spectre_v2 spec_store_bypass l1tf mds swapgs itlb_multihit srbds
bogomips        : 4799.99
clflush size    : 64
cache_alignment : 64
address sizes   : 40 bits physical, 48 bits virtual
power management:
 
Memory Info : 
MemTotal:        1811320 kB
MemFree:          175728 kB
MemAvailable:    1250400 kB
Buffers:             592 kB
Cached:          1256232 kB
SwapCached:            0 kB
Active:           626256 kB
Inactive:         723016 kB
Active(anon):     141852 kB
Inactive(anon):    59420 kB
Active(file):     484404 kB
Inactive(file):   663596 kB
Unevictable:           0 kB
Mlocked:               0 kB
SwapTotal:             0 kB
SwapFree:              0 kB
Dirty:                32 kB
Writeback:             0 kB
AnonPages:         90632 kB
Mapped:            82324 kB
Shmem:            108824 kB
KReclaimable:     113644 kB
Slab:             215664 kB
SReclaimable:     113644 kB
SUnreclaim:       102020 kB
KernelStack:        6256 kB
PageTables:        11136 kB
NFS_Unstable:          0 kB
Bounce:                0 kB
WritebackTmp:          0 kB
CommitLimit:      905660 kB
Committed_AS:    2725092 kB
VmallocTotal:   34359738367 kB
VmallocUsed:           0 kB
VmallocChunk:          0 kB
Percpu:             1544 kB
HardwareCorrupted:     0 kB
AnonHugePages:     20480 kB
ShmemHugePages:        0 kB
ShmemPmdMapped:        0 kB
HugePages_Total:       0
HugePages_Free:        0
HugePages_Rsvd:        0
HugePages_Surp:        0
Hugepagesize:       2048 kB
Hugetlb:               0 kB
DirectMap4k:      190320 kB
DirectMap2M:     1857536 kB

Python : 3.9.4
ruptures : 1.1.3
numpy : 1.20.2
scipy : 1.6.2
@oboulant this may be related to Mac OS Catalina vs Big Sur. @josh-boehm and I are both running Big Sur. Maybe @deepcharles is as well?

This could be related to the updated C libraries in Big Sur. 
Indeed I am using Big Sur. I will look into it.
FYI, my MacOS laptop is now on Big Sur. 

I now manage to get : 

* `{2}`
* `{2, 4}`
* `{1, 2}`

Surprisingly, I cannot get `{1, 2, 4}` even if I increase a lot the number of iteration of the for loop. 

Fix patch:
diff --git a/src/ruptures/detection/_detection/ekcpd_pelt_computation.c b/src/ruptures/detection/_detection/ekcpd_pelt_computation.c
index 2a70dd62..6ca5cda1 100644
--- a/src/ruptures/detection/_detection/ekcpd_pelt_computation.c
+++ b/src/ruptures/detection/_detection/ekcpd_pelt_computation.c
@@ -63,7 +63,7 @@ void ekcpd_pelt_compute(double *signal, int n_samples, int n_dims, double beta,
             c_r += kernel_value_by_name(&(signal[s * n_dims]), &(signal[(t - 1) * n_dims]), n_dims, kernelDescObj);
             S[s] += 2 * c_r - diag_element;
         }
-        c_cost = D[t] - D[s] - S[0] / t;
+        c_cost = D[t] - D[0] - S[0] / t;
         M_V[t] = c_cost + beta;
     }
 


Test patch:
diff --git a/tests/test_detection.py b/tests/test_detection.py
index 77d91f56..79560915 100644
--- a/tests/test_detection.py
+++ b/tests/test_detection.py
@@ -316,6 +316,25 @@ def test_kernelcpd_small_signal(signal_bkps_5D_n10, kernel):
     )
 
 
+@pytest.mark.parametrize("kernel", ["linear", "rbf", "cosine"])
+def test_kernelcpd_small_signal_same_result(signal_bkps_5D_n10, kernel):
+    signal, _ = signal_bkps_5D_n10
+    algo = KernelCPD(kernel=kernel)
+    list_of_segmentations = list()
+    n_iter = 100
+    for _ in range(n_iter):
+        bkps = algo.fit(signal=signal).predict(pen=1.0)
+        list_of_segmentations.append(bkps)
+
+    # test if all segmentations are equal
+    first_bkps = list_of_segmentations[0]
+    all_elements_are_equal = all(
+        first_bkps == other_bkps for other_bkps in list_of_segmentations[1:]
+    )
+    err_msg = "KernelCPD returns different segmentations on the same signal."
+    assert all_elements_are_equal, err_msg
+
+
 @pytest.mark.parametrize(
     "algo, model",
     product(