% despiking algorithm profiling (manuscript: "Computational Efficiency for the Surface Renewal Method")
% Jason Kelley, NEWAg lab
% January 2017
% contact: kelleyja@oregonstate.edu

%% load data
%cd('D:\')
cd('C:\Users\kelleyja\Data')
[file, path] = uigetfile; cd(path);
EC_DATA = load(file,'data');
EC_DATA = EC_DATA.data;
EC_DATA = errNaN(EC_DATA,[3 10 0]);%SR data  %[1 5 0]);%CO data

%    RhoV = EC_DATA.H2O;
%  Tsonic = EC_DATA.Tsonic;

%% time series parameters and subsample the data
 TIME.freq = 10;                  	% frequency (Hz)
tpd = TIME.freq*60*60*24;
tph = TIME.freq*60*60;

sec = 60;
min = 60;
hrs = 8;
day = 1;
 
%{x
    % advance data to desired start timestamp here
            % frequency x secs x min x hrs x days
      tt =  1:TIME.freq * sec * min * hrs * day ;
    %tt =  1440*600+1:1440*600+TIME.freq *  60  * 60  * 4  *  1  ; % 1:10^6;
    DATA = truncStruct(EC_DATA,tt);
    rhov = EC_DATA.H2O(tt); Ts = EC_DATA.Tsonic(tt);
    timestamp = DATA.TIMESTAMP;
%}

%%%%%%%%%%%%%% end of user parameters %%%%%%%%%%%%%%%%

TIME.tsmp = 1/TIME.freq;          	% sampling time (s)
TIME.beg  = datevec(timestamp(1));
TIME.fin  = datevec(timestamp(end));
TIME.n    = size(timestamp,1);
TIME.ttot = TIME.n*TIME.tsmp/60;	% length of record (minutes)

% Calculate density of air and heat capacity (to convert units of fluxes)
Tc = [DATA.cell_tmpr];
TK = Tc + 273.15;
DATA.celP = [DATA.cell_press]; 
ea = (rhov./1000).*0.287.*TK./0.622;            % [kPa] partial pressure water
rhod = (DATA.celP-ea)./0.287./TK;             	% [kg/m^3] density of dry air
rho = rhov./1000+rhod;                          % [kg/m^3] density of moist air
q = rhov./1000./rho;                            % [kg/kg] specific humidity
esat = 0.6112.*(exp(17.27.*Tc./TK));
wind = [DATA.Ux, DATA.Uy, DATA.Uz];
% EC_DATA.RH  = ea./esat.*100;
% EC_DATA.RCP = rho_cp(Tc, celP, rhov);            	
clear Tc TK P rhov ea rhod rho esat

%clear EC_DATA

if not(exist('piter','var')); piter = 1; end

%% despiking raw time trace data - Vickers Mahrt method

for piter = 60:62;

profile clear
profile('on','-detail','builtin','-memory','-timer','cpu')
tic
dsVM = despikeVM(timestamp,Ts,200,4,10);
% catalog performance
 %cond = ['dsVM-lg',num2str(maxlag),'-av',num2str(TIME.tavSR),'-np',num2str(TIME.npSR)];
%dirst = fullfile(path,'\despiketest\',cond);

%mkdir(dirst);
%profsave(profile('info'),dirst)
display('Vickers-Mahrt Method')
perftime = toc; toc

%profview
PERFSTATS = profile('info');
Prow = find(ismember({PERFSTATS.FunctionTable.FunctionName},'despikeVM'));

PERF(piter,1).n = TIME.n;
PERF(piter,1).t = perftime; 
PERF(piter,1).TotalTime = PERFSTATS.FunctionTable(Prow).TotalTime;
PERF(piter,1).TotalMemAllocated = PERFSTATS.FunctionTable(Prow).TotalMemAllocated;
PERF(piter,1).PeakMem = PERFSTATS.FunctionTable(Prow).PeakMem;

 %%%%%%%% Convolution method %%%%%%%

profile clear
profile('on','-detail','builtin','-memory','-timer','cpu')
tic

 dsCV = despike(Ts,200,4,3);
%  cond = ['dsCONV-lg',num2str(maxlag),'-av',num2str(TIME.tavSR),'-np',num2str(TIME.npSR)];
% dirst = fullfile(path,'\despiketest\',cond);

%mkdir(dirst);
%profsave(profile('info'),dirst)
display('Convolution Method')
perftime = toc; toc
%profview

PERFSTATS = profile('info');
Prow = find(ismember({PERFSTATS.FunctionTable.FunctionName},'despike'));

PERF(piter,2).n = TIME.n;
PERF(piter,2).t = perftime; 
PERF(piter,2).TotalTime = PERFSTATS.FunctionTable(Prow).TotalTime;
PERF(piter,2).TotalMemAllocated = PERFSTATS.FunctionTable(Prow).TotalMemAllocated;
PERF(piter,2).PeakMem = PERFSTATS.FunctionTable(Prow).PeakMem;


clear cond dirst perftime PERFSTATS Prow

end

%piter = piter+1;
%fprintf('Completed despiking for %0.2f hours, %i days data\nNext iteration %i\n',hrs+min/60,day-1,piter)

%% various plots

figure; semilogy([PERF(:,1).n],[PERF(:,1).TotalMemAllocated],'*'); hold('on'); semilogy([PERF(:,2).n],[PERF(:,2).TotalMemAllocated],'o'); grid
figure; semilogy([PERF(:,1).n],[PERF(:,1).PeakMem],'*'); hold('on'); semilogy([PERF(:,2).n],[PERF(:,2).PeakMem],'o'); grid
figure; semilogy([PERF(:,1).n],[PERF(:,1).TotalTime],'*'); hold('on'); semilogy([PERF(:,2).n],[PERF(:,2).TotalTime],'o'); grid
figure; semilogy([PERF(:,1).n],[PERF(:,1).t],'*'); hold('on'); semilogy([PERF(:,2).n],[PERF(:,2).t],'ro'); grid
figure; loglog([PERF(:,1).n],[PERF(:,1).t],'*'); hold('on'); loglog([PERF(:,2).n],[PERF(:,2).t],'ro'); grid
figure; loglog([PERF(:,1).n],[PERF(:,1).TotalTime],'*'); hold('on'); loglog([PERF(:,2).n],[PERF(:,2).TotalTime],'o'); grid
%%

figure; loglog([PERF(:,1).n]./tph,[PERF(:,1).TotalTime],'*'); hold('on'); loglog([PERF(:,2).n]./tph,[PERF(:,2).TotalTime],'o'); grid
%%
figure; loglog([PERF(:,1).n]./tph,[PERF(:,1).TotalTime]./[PERF(:,1).n],'*'); hold('on'); 
loglog([PERF(:,2).n]./tph,[PERF(:,2).TotalTime]./[PERF(:,2).n],'o'); grid
%%
figure; loglog([PERF(:,1).n]./tph,[PERF(:,1).TotalMemAllocated]./[PERF(:,1).n],'*'); hold('on'); 
loglog([PERF(:,2).n]./tph,[PERF(:,2).TotalMemAllocated]./[PERF(:,2).n],'o'); grid
%%
figure; loglog([PERF(:,1).n]./tph,[PERF(:,1).PeakMem]./[PERF(:,1).n],'*'); hold('on'); 
loglog([PERF(:,2).n]./tph,[PERF(:,2).PeakMem]./[PERF(:,2).n],'o'); grid
%%
figure; loglog([PERF(:,1).n]./tph,[PERF(:,1).t]./[PERF(:,1).n],'*'); hold('on'); 
loglog([PERF(:,2).n]./tph,[PERF(:,2).t]./[PERF(:,2).n],'o'); grid
%%
figure; loglog([PERF(:,1).n]./tph,[PERF(:,1).t]./[PERF(:,1).TotalTime],'*'); hold('on'); 
loglog([PERF(:,2).n]./tph,[PERF(:,2).t]./[PERF(:,2).TotalTime],'o'); grid

%% tick labels
ax = gca;
ax.XTick = [ 0.25 1 4 12 24 48];
ax.XTickLabel = { '15min'; '1h'; '4h'; '12h'; '24h'; '48h'};

%% sort and average by length of data
uts = unique([PERF(:,1).n]);
for u = 1:length(uts)
    uu = find([PERF(:,1).n]==uts(u));
    PERFave(u,1).n = uts(u);
    PERFave(u,1).t = nanmean([PERF(uu,1).t]);
    PERFave(u,1).tstd = nanstd([PERF(uu,1).t]);
    PERFave(u,1).ttot = nanmean([PERF(uu,1).TotalTime]);
    PERFave(u,1).ttstd = nanstd([PERF(uu,1).TotalTime]);
	PERFave(u,2).n = uts(u);
    PERFave(u,2).t = nanmean([PERF(uu,2).t]);
    PERFave(u,2).tstd = nanstd([PERF(uu,2).t]);
    PERFave(u,2).ttot = nanmean([PERF(uu,2).TotalTime]);
    PERFave(u,2).ttstd = nanstd([PERF(uu,2).TotalTime]);

end
%%
% figure; loglog([PERFave(:,1).n]./tph,[PERFave(:,1).ttot],'k^'); hold('on'); 
% loglog([PERFave(:,2).n]./tph,[PERFave(:,2).ttot],'rv'); grid

% figure; loglog([PERFave(:,1).n]./tph,[PERFave(:,1).t]./[PERFave(:,1).ttot],'k^'); hold('on'); 
% loglog([PERFave(:,2).n]./tph,[PERFave(:,2).t]./[PERFave(:,2).ttot],'rv'); grid

figure; loglog([PERFave(:,1).n]./tph,[PERFave(:,1).ttot]./[PERFave(:,1).n].*10^9,'k^'); hold('on'); 
loglog([PERFave(:,2).n]./tph,[PERFave(:,2).ttot]./[PERFave(:,2).n].*10^9,'rv'); grid

ax = gca;
ax.XTick = [ 0.25 1 4 12 24 48];
ax.XTickLabel = { '15min'; '1h'; '4h'; '12h'; '24h'; '48h'};