KL divergence seems to work in principle, not quite right yet.

prob_laplace.m

@@ -418,7 +418,31 @@ plot(ls, sum(p_o_i, 2));
 %% Save time points for comparison with python
 csvwrite('lb_in_out.csv', ls)
 csvwrite('prob_in_out.csv', p);
+%%
+KL = zeros(1, 20);
+for i = 1:20
+i_t = 1+i;
+j_t = 2+i;
+integrand = @(p, q) abs(p.*log(p./q));
+integrate = @(x, y) sum(diff(x)'.*(y(2:end)+y(1:end-1))/2);
+n_i = @(p1, p2, p3, p4) abs(integrate(ls, p1) + integrate(ls, p2)+...
+                            integrate(ls, p3) + integrate(ls, p4));
+N_i = n_i(p_i_i(:, i_t), p_o_o(:, i_t), p_i_o(:, i_t), p_o_i(:, i_t));
+N_j = n_i(p_i_i(:, j_t), p_o_o(:, j_t), p_i_o(:, j_t), p_o_i(:, j_t));
 
+KL(i) = integrate(ls, integrand(p_i_i(:, i_t)/N_i, p_i_i(:, j_t)/N_j))+...
+     integrate(ls, integrand(p_o_o(:, i_t)/N_i, p_o_o(:, j_t)/N_j))+...
+     integrate(ls, integrand(p_i_o(:, i_t)/N_i, p_o_i(:, j_t)/N_j))+...
+     integrate(ls, integrand(p_o_i(:, i_t)/N_i, p_i_o(:, j_t)/N_j));
+end
+plot(t(2:21), abs(KL));
+ylim([0, 1])
+%%
+figure; hold on;
+plot(ls, p_i_i(:, i_t)/N_i);
+plot(ls, p_o_o(:, i_t)/N_i);
+plot(ls, p_i_o(:, i_t)/N_i);
+plot(ls, p_o_i(:, i_t)/N_i);
 %% %%%%%%%%%%%%%%%%%%% INTEGRATION FUNCTION DEFINITIONS %%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%