function [barcodeGen_periodic,barcodeGen_notperiodic,barcodeGen_edgy,barcodeGen_all] = periodic_finder(barcodeGen,minreps,plots)
% Periodic finder
% This is a second version which uses the following procedure:
% - Autocorrelation function to find repetitive stuff
% - low coef of variation between peaks in the autocorr for repetitive vs non
% - Fast fourier transform to double check period
% - Rolling window corr to find repetitive stuff within repetitives
% - findchangepts for partial vs fullrepeats

%% Parameters

if nargin <2
  minreps=5;% was 3 before but many regular barcodes have short rep regions in within, this solves it
end
if nargin <3
  plots=0;% plots used for testing and validation
end

w=150; % w from always... 150px for comparisons, used only to classify between fully rep and partially rep ( if regular region <w then fully rep)
cov_thres=0.2; % 20% we use coef of variation instead of std, maybe could go lower
mincorr=0.4;% works between 0.4 and 0.5
minfreq=0.08; % might drop FFT since I will not be able to comfortably handle longer rep lengths
laglimit=110; %how many lags we take in ACF? %arbitrary, easier to interpret if we take less than this (not too much less!)
maxchanges=2;% how many segments will a molecule have? regular-repetitive-regular or rep-reg-rep at most?
minpeakacf=-2; % mean peak height in ACF, virtually -inf
acfpeaks=6; % how many peaks to take from ACF
edgeignore=6; %pixels to ignore from the sides
period_tolerance=0.1; % 10% range tolerance of period agreement between ACF and FFT
detrend_bool=1;

%% Alocating variables
barcodeGen_periodic = [];
barcodeGen_notperiodic = [];
barcodeGen_edgy = [];
barcodeGen_all = [];
rep_count = 0;
notrep_count = 0;
edgy_count = 0;

%% Actual loop 

warning('off', 'all');

for  k= 1:length(barcodeGen) 
% Step 0: Ignore edges
if length(barcodeGen)==1
x=barcodeGen.rawBarcode(edgeignore+1:end-edgeignore);
else
x=barcodeGen{k}.rawBarcode(edgeignore+1:end-edgeignore); 
end
% Step 1: Compute the autocorrelation
[acf, lags] = autocorr(x, 'NumLags', length(x)-1);

% The first many lags make sense, then stops making sense, so crop down

if length(x)<laglimit
laglimit=length(x);
end
acf=acf(1:laglimit);
lags=lags(1:laglimit);

% Sinusoidal curves in acf are periodic barcodes, but maxima can be hard to
% detect in some, so we detrend them, a more serious detrend could be
% better
if detrend_bool==1
acf=detrend(acf);
end


% Step 2: Detect peaks in the ACF
[peaks_acf, peak_lags] = findpeaks(acf, lags, 'MinPeakHeight', minpeakacf);

% I am only taking the first couple peaks, it kinda looses meaning
% afterwards
if length(peak_lags)>acfpeaks
    peak_lags=peak_lags(1:acfpeaks); % peak position
    peaks_acf=peaks_acf(1:acfpeaks); % peak magnitude (not really neded unless plotting)
end


% Step 3: Calculate the period (average distance between peaks)

if length(peak_lags) > 2 % If not even 3 peaks not perioidc
    period = mean(diff(peak_lags));
    % fprintf('Detected period by acf: %.2f\n', period);

   if std(diff(peak_lags))/ mean(diff(peak_lags))<cov_thres % peaks should be around equally spaced, if so, at least its edgy, if not, then not periodic
    periodic_flag=1;
    else
    periodic_flag=0;
    end    
else
    period=100; % just a place holder so there is never not a period
    % fprintf('No significant periodicity detected.\n');
    periodic_flag=0;
end

%  Period has been detected (or not) by autocorr
%  it should be around the same as FFT

% Step 4: Compute FFT
X = fft(x-mean(x));  % Detrending by removing the mean
n = length(X);
frequencies = (0:n-1)*(1/n);  % Frequency vector
amplitude = abs(X/n);

% Step 5: Identify dominant frequency
% Low Freqs will show peaks because of PSF, this while FFT approach could
% be removed in a future version to accept larger repeats, in place try to
% go for a "local" autocorr
halfamplitude=amplitude(1:floor(length(amplitude)/2)); % Freq domain is symmetric, so let's just take the first half
[peaks, locs] = findpeaks(halfamplitude,'MinPeakProminence',5,'Annotate','extents');

% Not all frequencies are really needed, there are some peaks below 0.05, a
% the upper limit feels useless after the half amplitude... should I
% discard that?
valid_freq_range = locs > round(minfreq*length(amplitude)) &  locs <round(0.5*length(amplitude));
peaks_subset=peaks(valid_freq_range);
locs_subset=locs(valid_freq_range);
[~,maxloc]=max(peaks_subset);
max_idx = locs_subset(maxloc);
dominant_frequency = frequencies(max_idx );  

% Step 6: Calculate period from dominant frequency
if isempty(dominant_frequency)
dominant_frequency=0.05; % this default is just so the rest of calculations is done in 20px windows
end

period_fft = 1 / dominant_frequency;
%fprintf('Dominant Frequency: %.2f\n', dominant_frequency);
%fprintf('Estimated Period: %.2f pixels\n', period_fft);

% Here we do, if period by fft is not close to period autocorr by 10%ish,
% so if periods do not agree, not periodic!
if abs(period_fft-period)>period_tolerance*period
periodic_flag=0;
end

% If periodic so far then further analysis
if periodic_flag==1

% Step 7: Identify detected period
% Rolling window approach to calculate correlation
% each pixel gets a correlation value of its window to the next one
% pixels from the later half get compared to the window before 
% This makes the plot symmetric just in the mid 2*windowsize pixels

window_size = round(period_fft);  % Window size equal to one period
correlation_values = zeros(length(x),1);

 for i = 1:length(x)
     if i < length(x)/2
     segment = x(i:i + window_size - 1);
     correlation_values(i) = corr(segment', x(i+window_size:i + 2*window_size-1)');
     else 
     segment = x(i - window_size + 1:i);
     correlation_values(i) = corr(segment', x(i - 2*window_size+1:i-window_size)');
     end
 end

% Find the start and end indices of the periodic segments
% findchangepts finds changes in the mean corr value
% corr values are high for periodic stuff low for non periodic

[cp]=findchangepts(correlation_values,"MaxNumChanges",maxchanges);
% This slices our plots in several sections

a=[1;cp;length(correlation_values)];

% Step 8: Identify and index rep_region
% Take segments of at least X repetitions big (3-5)
% If region's meancorr > 0.5 or  "min mean" then its a rep region
% Cycle through all the inner segments of each barcode
rep_region=[];
meansegments={};
for seg = 1:(length(a)-1)
if length(a(seg):a(seg+1))>=window_size*minreps
    meansegments{seg}=mean(correlation_values(a(seg):a(seg+1)));  
    if meansegments{seg}>mincorr
    rep_t_region=a(seg):a(seg+1);
    rep_region=[rep_region,rep_t_region];
    end
else
    meansegments{seg}=mean(correlation_values(a(seg):a(seg+1)));
end

end

% Step 9: Finally classify as rep, edgy or normal 
% If there is at least 1 rep region, its at least an edgy barcode
if isempty(rep_region)==0 
% If the non rep region is small (smaller than w) then is full rep
    if length(x)-length(rep_region) < w
    periodic_flag=2;
    else
    periodic_flag=1;
    end
else
periodic_flag=0;    
end

% Highlight the periodic part on the original signal
% plotting purposes only
%periodic_part = x(rep_region);

% Step 10: Save rep_region as repmask
rep_region=rep_region+edgeignore;
repmask=zeros(1,length(x)+edgeignore*2);
repmask(rep_region)=1;

else % STUFF TO DO FOR NON REPETITIVES
rep_region=[];
rep_region=rep_region+edgeignore;
correlation_values=[];
repmask=zeros(1,length(x)+edgeignore*2);
repmask(rep_region)=1;
meansegments={};
end

% Step 8: Save back to structure
   if length(barcodeGen)==1
   barcodeGen_all{1,k} = barcodeGen;
   else
   barcodeGen_all{1,k} = barcodeGen{1,k};
   end
   barcodeGen_all{1,k}.repmask = repmask;

   if periodic_flag == 2
    rep_count = rep_count + 1 ;
    if length(barcodeGen)==1
    barcodeGen_periodic{1,rep_count} = barcodeGen;
    else
    barcodeGen_periodic{1,rep_count} = barcodeGen{1,k};
    end
    barcodeGen_periodic{1,rep_count}.repmask = repmask;
    barcodeGen_periodic{1,rep_count}.rep_check = 2;
    barcodeGen_all{1,k}.rep_check = 2;  
   elseif periodic_flag == 1
    edgy_count= edgy_count + 1 ;
    if length(barcodeGen)==1
    barcodeGen_edgy{1,edgy_count} = barcodeGen;
    else
    barcodeGen_edgy{1,edgy_count} = barcodeGen{1,k};
    end    
    
    barcodeGen_edgy{1,edgy_count}.repmask = repmask;
    barcodeGen_edgy{1,edgy_count}.rep_check = 1; 
    barcodeGen_all{1,k}.rep_check = 1;  
   
   else 
    notrep_count= notrep_count + 1 ;
    if length(barcodeGen)==1
    barcodeGen_notperiodic{1,notrep_count} = barcodeGen;
    else
    barcodeGen_notperiodic{1,notrep_count} = barcodeGen{1,k};
    end   

    barcodeGen_notperiodic{1,notrep_count}.repmask = repmask;
    barcodeGen_notperiodic{1,notrep_count}.rep_check = 0;  
    barcodeGen_all{1,k}.rep_check = 0;  
   end 

%% Test plots %%
% ALL PLOTS
% 
if plots
figure;
t = tiledlayout(4,2,'TileSpacing','compact','Padding','compact');

% ---- Plot 1 ----
ax1 = nexttile(1);
plot(lags, acf); hold on;
plot(peak_lags, peaks_acf, 'ro');
%title('Autocorrelation and frequency estimation');
xlabel('Lag (px)');
ylabel('Autocorrelation');

lgd1 = legend(ax1, {'Autocorrelation function','Peaks'});
lgd1.Layout.Tile = 2;
lgd1.Location = 'southeastoutside';



% ---- Plot 2 ----
% fft and most likely freq
ax2 = nexttile(3);
if length(correlation_values)>2
    plot(frequencies, amplitude); hold on;
    plot(frequencies(max_idx), amplitude(max_idx),'ro');
end

%title('Fast fourier transform');

lgd2 = legend(ax2, {'Frequency spectrum','Repetitive frequency'});
lgd2.Layout.Tile = 4;
lgd2.Location = 'southeastoutside';
xlabel('Frequency (1/px)');
ylabel('Amplitude');


% ---- Plot 3 ----
% Let's show here how we decide which segment is repetitive
% rolling window autocorrelation
ax3 = nexttile(5);
plot(1:length(correlation_values), correlation_values,'-b'); 
hold on;

if length(correlation_values)>2
    yline(mincorr,'k-','LineWidth', 1);% Plot horizontal line for threshold
    for seg = 1:length(cp)+1
     if seg == 1
        line([a(seg) a(seg+1)],[meansegments{seg} meansegments{seg}],'Color','green','LineStyle','--', 'LineWidth', 1,'DisplayName','Segment mean'); % Plot horizontal lines for means of segments
     else
        line([a(seg) a(seg+1)],[meansegments{seg} meansegments{seg}],'Color','green','LineStyle','--', 'LineWidth', 1, 'HandleVisibility','off'); % Plot horizontal lines for means of segments        
     end
    end
    for i = 1:length(cp)
        xline(cp(i), '--', 'LineWidth', 1,'Color','red'); % Plot vertical lines at change points
    end
end

%title('Change Points in rolling window autocorrelation');
xlabel('Position (px)');
ylabel('Correlation Coefficient');
ylim(ax3,[-1,1]);
lgd3 = legend(ax3, {'Rolling window correlation','Correlation threshold','Segment mean','Change points'});
lgd3.Layout.Tile = 6;
lgd3.Location = 'southeastoutside';


% ---- Plot 4 ----
% Show the whole barcode with detected repeat last
ax4 = nexttile(7);
plot(1:length(x), x); hold on;
plot(rep_region-edgeignore, x(rep_region-edgeignore), 'r');
%title(sprintf('Barcode %.2f with detected repetitive region',k));
xlabel('Position (px)');
ylabel('Intensity (a.u.)');

lgd4 = legend(ax4, {'Non-repetitive','Repetitive region'});
lgd4.Layout.Tile = 8;
lgd4.Location = 'southeastoutside';
end % end of plots

end


fprintf('Total barcodes: %.2f\n', length(barcodeGen));
fprintf('Non-repetitive barcodes: %.2f\n', notrep_count);
fprintf('Edgy-repetitive barcodes: %.2f\n', edgy_count);
fprintf('Full repetitive barcodes: %.2f\n', rep_count);



warning('on', 'all');
end