Improvements and bugfixes

- Improved the runtime of the linear extension subroutine for small signals

Author: mstorath

cssd.m

@@ -263,7 +263,7 @@
                 end
                 %%% END PELT PRUNING
 
-            %%% BEGIN FPVVI PRUNING
+            %%% BEGIN FPVI PRUNING
             otherwise
                 for rb=3:N
 

demos/Ex_Stock_CV.m

@@ -73,11 +73,28 @@
 
 %% run again with the PELT pruning for runtime comparison
 rng(123)
-
-% starting the optimization with these fine-tuned initial values turned out to 
-% result in a reasonbly good CV score for this data
-p_init = 0.4702;  
-gamma_init = 0.0069;
 tic
 output_cv = cssd_cv(x, y, [], [], [], [p_init; gamma_init], 'verbose', 1, 'pruning', 'PELT');
-runtime_pelt = toc
+runtime_pelt = toc
+
+%% plot result
+xx = linspace(min(x), max(x),100000);
+fig = figure(1);
+clf;
+set(gcf, 'name', 'Stock', 'Color', 'white', 'units','normalized', 'position', [0,0,0.3,0.3]);
+plot([],[])
+hold on
+output.pcw_fun.plot(xx, '-', 'Linewidth', 2, 'Color', '#77AC30', 'DisplayName', 'CSSD')
+plot(x, y, '.k', 'Markersize', 0.5, 'DisplayName', 'Data')
+legend('Location', 'NorthWest')
+xlim([min(xx), max(xx)])
+for i = 1:numel(discont)
+    plot([discont(i), discont(i)], ylim, '--', 'Color', '#999999', 'DisplayName','')
+end
+xticks(discont);
+xtickangle(60)
+dates_ext_str = datestr(dates_ext);
+xticklabels(dates_ext_str(round(xticks),:));
+hold off
+leg = legend;
+legend(leg.String{1:3})

demos/ruptures_cssd.py

@@ -5,6 +5,8 @@
 import time
 import cProfile
 import pstats
+import math
+
 
 class CSSD_baseline(BaseCost):
 
@@ -64,8 +66,11 @@ def detect_changepoints(x,y, p, gamma, delta):
     cost.num_elems = 0
 
     # invoke PELT method
-    algo = rpt.Pelt(custom_cost=cost, min_size=1, jump=1).fit(y)
-    result = algo.predict(pen=gamma)
+    if math.isinf(gamma):
+        result = []
+    else:
+        algo = rpt.Pelt(custom_cost=cost, min_size=1, jump=1).fit(y)
+        result = algo.predict(pen=gamma)
     elapsed_time = time.time() - start_time
 
     return result, elapsed_time, cost.num_calls, cost.num_elems

subroutines/linext_pp.m

@@ -3,22 +3,28 @@
 %to the boundaries [l, r]
 
 assert( (l <= pp.breaks(1)) && (pp.breaks(end) <= r), 'New boundaries must be larger than the old ones for extension.')
+pp = embed_pptocubic(pp);
 
-% extends the boundaries by one in each direction and continues the spline
-% linearly (fnxtr is Matlab built-in)
-pp = fnxtr(pp,2);
+pp_deriv = pp;
+pp_deriv.coefs = pp_deriv.coefs(:, 1:3) .* [3,2,1];
+pp_deriv.order = 3;
+first = pp.breaks(1);
+last  = pp.breaks(end);
 
-% adjust the right bound 
-pp.breaks(end) = r;
+new_breaks = [l, pp.breaks, r];
+base_last = ppval(pp, last);
 
-% when changing the first break, the base point of the polynomial is changed
-% which requires a corrected     
-for i = 1:pp.dim % loop over all dimensions
-    a = pp.coefs(i,end-1);
-    pp.coefs(i,end) = pp.coefs(i,end) - a * (pp.breaks(1) - l);
-    % finally adjust the endpoints
-end
-pp.breaks(1) = l;
+slope_last = ppval(pp_deriv, last);
+
+base_first = ppval(pp, first);
+slope_first = ppval(pp_deriv, first);
+base_l = base_first + slope_first .* (l - first);
+
+new_coefs  = [zeros(pp.dim, 2), slope_first, base_l; pp.coefs; zeros(pp.dim, 2), slope_last, base_last];
+
+pp.breaks = new_breaks;
+pp.coefs = new_coefs;
+pp.pieces = pp.pieces + 2;
 
 end