%% Preliminaries
N = 1024; 
x = 32*pi*(1:N)'/N; 
u = cos(x) + 0.1*cos(x/16).*(1+2*sin(x/16));
dt = 1/4;
tmax = 1500;

%% Solve KS eqn w/ adapted kursiv.m from https://doi.org/10.1137/S1064827502410633
[U,t] = kursiv(x,u,dt,tmax);

%% Find stripe width using xcorr (xcorr2 is needlessly expensive)
foo = nan(1,numel(t));
for j = floor(numel(t)/2):numel(t)
    % Determine typical grid offset between minima and maxima for each time
    xc = xcorr(U(:,j),-U(:,j));
    [~,lag] = max(xc);
    foo(j) = abs(lag-numel(x));
end
stripeGridWidth = round(mean(foo(~isnan(foo))));

%% Plot solution of KS equation
figure; 
pcolor(t,x,U); 
shading flat;

%% Attempted stripe definition 1
% Stripes are defined by the bottom quartile and the top quartile after a
% suitable circular shift: besides fragments (which bwareafilt removes
% along with some regions on the [unidentified] boundary]), this leaves
% some short time gaps on occasion
below = U<quantile(U(:),.25);
aboveShifted = circshift(U,stripeGridWidth,1)>quantile(U(:),.75);
stripes = below|aboveShifted;
% Filter small areas while avoiding artificial boundaries
stripesA = bwareafilt(stripes,[1e3,Inf],4);
shift = round(numel(x)/2);
stripesB = circshift(bwareafilt(circshift(stripes,shift,1),[1e3,Inf],4),-shift,1);
stripes = stripesA|stripesB;
%
figure;
pcolor(t,x,1*stripes);
shading flat;
colormap(flipud(gray));

%% Attempted stripe definition 2
foo = U-circshift(U,-stripeGridWidth,1);
stripes2 = foo>-median(foo);
% Filter small areas while avoiding artificial boundaries
stripes2A = bwareafilt(stripes2,[1e3,Inf],4);
shift = round(numel(x)/2);
stripes2B = circshift(bwareafilt(circshift(stripes2,shift,1),[1e3,Inf],4),-shift,1);
stripes2 = stripes2A|stripes2B;
%
figure;
pcolor(t,x,1*stripes2);
shading flat;
colormap(flipud(gray));


% % Plot stripes
% figure;
% d1con = diff(U(2:end,:),1,1)<0; % constraint on first derivative
% isStripe = 1*islocalmax(diff(abs(U),2,1),1).*d1con;
% pcolor(t,x(2:end-1),isStripe); 
% shading flat;
% colorbar;
% colormap(flipud(gray));
% 
% % Plot number of stripes
% figure; 
% plot(sum(isStripe,1),'k');
% 
% figure;
% pcolor(t,x(2:end-1),isStripe); 
% shading flat;
% colorbar;
% colormap(flipud(gray));
% hold on;
% contour(t,x,U,0*[1,1],'g');
% 
% 
% figure;
% pcolor(t,x(2:end-1),isStripe); 
% shading flat;
% colorbar;
% colormap(flipud(gray));
% hold on;
% contour(t,x(2:end-1),U(2:end-1,:).*d1con,0*[1,1],'g');
% 
% 
