% Test code for runlength.m
%
% 2005 Jan 10
% David S. Ciochetto
%
% This m-file was used by commenting/un-commneting lines of code and activating them with the 'f9' key.
% Comments on the resutls left commented here should follow the important lines of code.
%
% Q: Does the result of generating random numbers evenly distributed reproduce the theoretical PDF?

clear; close all;

% Pick a number of points n, signal length to generate
n = 1e5;

% Generate a vector length n of random data
G = rand(n,1);
% Make the data mean zero (rand produces n uniformly distributed random numbers from 0 to 1)
G = G - 0.5; %mean(G); % Check the mean value compared to 0.5.
% Try Gaussian Distributed Noise
%G = randn(n,1); G = G - 0.5; %mean(G);
% Does not apply to Gaussian distributed noise, only uniformly distributed noise.

% Call runlength and get the result in res
res = runlength(G);

% Visualize the results
% Make an x-axis vector
N=[1:max(res)];
% Put the data into a histogram with integer bins (based on integer values set by N)
H=hist(res,N);
%bar(N,H/n); %y = H./n; % Wrong way to normalize
% Normalize the results to the total area, make it a probability from 0 - 1.
y = H/sum(H);
bar(N,y);
set(gca,'FontSize',12);
title(['Test of runlength.m with n = ' num2str(n) ' points. Yellow is theoretical; P(n) = 2^{-n}.'])
xlabel('n'); ylabel('Probability of run of length n, P(n)  [0-1]')
hold on; plot(N,2.^(-N),'y.-'); hold off;
% How does it look plotted on semi-log paper?
pause
semilogy(N,y,'kx-',N,2.^(-N),'y.-')
set(gca,'FontSize',12);
title(['Test of runlength.m with n = ' num2str(n) ' points. Yellow is theoretical; P(n) = 2^{-n}.'])
xlabel('n'); ylabel('Probability of run of length n, P(n)  [0-1]')

clear n G res N y H




% Q: What about a linear density profile with random noise added, will it return the predicted PDF?
close all

% Set up a realistic stratification
g = 9.81;			% Gravitational acceleration [m s^-2]
rho_r=1026;		% Reference Density [kg m^-3]
drho=2;				% Density difference over depth difference dz [kg m^-3]
dz=-10;				% Depth differnece for drho [m]

% CHECK: Calculate N^2 for the setting above. Linear profile so a constant N^2
N2 = -1*g/rho_r*drho/dz;

% Generate a linear profile: Profile Settings 
rho_1 = 1023;						% Surface density [kg m^-3]
z = [-100:0];						% Independent depth axis [m] (Increase this to add more points to the analsis
rho=rho_1 +drho/dz.*z;	% Density profile (dependent variable) (This is the sorted profile)

% Plot it to see
plot(rho,z,'k:')
set(gca,'FontSize',12);
% Add random noise that is mean zero.
% This is the fake signal for a linearly stratified fluid with random noise.
amp = 2;	% Allow one to adjust the amplitude of the noise.
% Q: Does adjusting the amplitude of the random noise increase the occurance of points falling outside of
% the threshold value for small sample sizes (large n mainly)?
r = rho+amp*(rand(1,length(rho))-0.5);

% Show the signal superimposed on the non-noisy linear stratification.
hold on; plot(r,z,'r-'); hold off
title('Simulated density signal for a linear stratification with uniformly distributed random noise (solid)')
ylabel('z [m]'); xlabel('\rho [kg m^{-3}] (non-noise profile is the dotted line)')
xL = get(gca,'XLim'); yL = get(gca,'YLim');
x = min(xL) + diff(xL)/10; y = min(yL) + diff(yL)/10;
text(x,y,['N = ' tsigfig(sqrt(N2),4) '[rad s^{-1}],  T = ' tsigfig(2*pi/sqrt(N2),4) '[s]'],'FontSize',12);

pause
% Check the runlength for the simulated signal.
close
% Calculate the runlength of the signal minus the trend (re-ordered profile)
res = runlength(r-rho);
% Plot a histogram of the resulting PDF
N=[1:max(res)]; H=hist(res,N); y = H/sum(H); bar(N,y);
set(gca,'FontSize',12);
title('Runlength analysis of fluctuating density signal for linear stratification with uniform noise')
xlabel('n'); ylabel('Probability of run of length n, P(n)  [0-1]')
% Add the theoretical distribution in yellow
hold on; plot(N,2.^(-N),'y.-'); hold off;
pause

% Make a semilog plot of this PDF including the threshold value.
semilogy(N,y,'kx-',N,2.^(-N),'y.-',N,2*2.^(-N),'r:')
set(gca,'FontSize',12);
title('Runlength analysis of fluctuating density signal for linear stratification with uniform noise')
xlabel('n'); ylabel('Probability of run of length n, P(n)  [0-1]')
%xL = get(gca,'XLim'); yL = get(gca,'YLim');
%x = min(xL) + diff(xL)*0.05; y = min(yL) + diff(yL)*0.002;
%text(x,y,'Solid Line with dots is theoretical; P(n) = 2^{-n}, dotted line is GK96 threshold value','FontSize',12);
xlabel('n (Solid Line with dots is P(n) = 2^{-n}, dotted line is GK96 threshold value.)')

% May have troubles with this method when the sample sizes are smaller. There is definitely increased
% uncertainties if they could be calculated.
