function DataChromatix(varargin)
%
% DataChromatix(data, timeunit, yunit, displaywin, shift, memorywin, ...
% colorwin, nbins, ncolors, lambda, fs, videoname, minsec, cmap)
%
% Given any time series, creates a video of its chromatized version
% based on the probability distribution of the previous values
% on a certain window (memory). The memory can start from the beginning
% of the recording or be specified by the user.
%
% inputs:
% data --> 1-column matrix containing time series to colorize, in any units
% timeunit --> string containing the label for the time axis
%              if you choose 1 for minsec, it overrides this one
% yunit --> string containing the label for the y axis (physical units of
%           data)
% displaywin --> window to display, in same units as data
% shift --> shift from one window to the next, in same units as data
% memorywin --> memory used to colorize time series, in same units as data.
%               If left empty,infinite memory is employed
%               default=memory starting from beginning of recording
% colorwin --> length of window that is colorized at each iteration,
%              in same units as data. Must be >shift
%              if left empty, coincides with shift
% nbins --> number of bins of the histogram
%           if left empty, it is estimated via the Freedman-Diaconis rule
% ncolors --> number of colors used in the colorization (default=64)
% lambda --> smoothing parameter for histogram (default=10)
% fs --> sampling frequency of the data (default=1)
% videoname --> string containing name of the video to be saved
%               (default='Myvideo')
%               input 0 if you don't want the video, only the final figure
% minsec --> flag, if =1 shows time in mm:ss (default=0), assuming that
%            the original time series is in seconds
% cmap --> colormap you want to use (default='jet')
%
% Examples: please note that the wfdb library must be installed
% and in your Matlab path for the examples to work
%
% % Example 1: Fixed-length memory window
% wfdb2mat('ctu-uhb-ctgdb/1170')
% load 1170m
% data=val(1,:)';
% data=data/100;
% data(data<30 | data>240)=[];
% lambda=10;
% fs=4;
% T=600.0;      % sec time interval to display (10min)
% nbins=100;
% ncolors=64;
% colorwin=100;
% Tshift=60.0; % sec time shift
% DataChromatix(data, 'Time (mm:ss)', 'FHR (bpm)',T, Tshift, 1000, colorwin, nbins, ncolors, 10, fs, 'video1', 1)
%
% % Example 2: Memory window starting from the beginning of recording
% wfdb2mat('ctu-uhb-ctgdb/1170')
% load 1170m
% data=val(1,:)';
% data=data/100;
% data(data<30 | data>240)=[];
% lambda=10;
% fs=4;
% T=600.0;      % sec time interval to display (10min)
% Tshift=60.0; % sec time shift
% DataChromatix(data, 'Time (s)', 'FHR (bpm)',T, Tshift, [], [], [], [], 10, fs, 0, 0)
%
% Authors: Anton Burykin and Sara Mariani, 2015
% Last Modified: January 7th, 2015
%
% please report bugs/questions at sara.mariani7@gmail.com
%
% when using this script, please reference: 
% A. Burykin*, S. Mariani*, T. Henriques, T. Silva, 
% W. Schnettler, M. D. Costa** and A. L. Goldberger**, 
% "Remembrance of time series past: simple chromatic 
% method for visualizing trends in biomedical signals,"  
% Physiol.  Meas., vol. 36, pp. N95, 2015.

inputs={'data','timeunit','yunit','displaywin','shift','memorywin', ...
    'colorwin','nbins','ncolors','lambda', 'fs', 'videoname','minsec','cmap'};

if nargin>14
    error('Too many input arguments')
end

if nargin<5
    error('Not enough input arguments')
end

for n=1:nargin
    eval([inputs{n} '=varargin{n};'])
end
for n=nargin+1:14
    eval([inputs{n} '=[];'])
end

if isempty(lambda), lambda=10; end
if isempty(videoname), videoname='MyVideo'; end
if isempty(minsec), minsec=0; end
if isempty(cmap), cmap='jet'; end
if isempty(fs), fs=1; end
if isempty(ncolors), ncolors=64; end
if isempty(colorwin), colorwin=shift; end
if colorwin<shift
    error('The colorization window must be longer than the shift')
end

scrsz = get(0,'ScreenSize');

if videoname~=0
    vidObj = VideoWriter(videoname);
    vidObj.FrameRate = 1;
    open(vidObj);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
window=displaywin*fs;     % time interval in data points
shift=shift*fs;
memorywin=memorywin*fs;
colorwin=colorwin*fs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t=(1:length(data))'/fs;
% set max and min range for screen based on data range
lowl=min(data)-std(data)/2;
hil=max(data)+std(data)/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Lhr=length(data);
Jmax=fix((Lhr-window)/shift);

hr_disp=[]';
t_disp=[]';
c_disp=[]';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:Jmax
    ind1=fix(1+shift*(j-1));
    ind2=fix(shift*(j-1)+window);
    if isempty(memorywin) % infinite memory
        hr_mem=data(1:ind2); % this is my memory window
    else
        % if the memory is finite
        if (ind2<=memorywin)
            ind_s=1; % like infinite memory
        else
            ind_s=ind2-memorywin;
        end
        hr_mem=data(ind_s:ind2);   % window you use for histogram
    end
    % histogram bin number estimation via Freedman-Diaconis rule
    if isempty(nbins)
        binw=2*iqr(hr_mem)/length(hr_mem)^(1/3);
        nbins=(max(hr_mem)-min(hr_mem))/binw;
    end
    H1=hist(hr_mem,nbins);
    H1=smooth1D(H1',lambda);
    H1=H1/max(H1);
    
    Nmax=length(H1);
    z=zeros(size(H1));
    Zmin=min(hr_mem);
    Zmax=max(hr_mem);
    
    for n=1:Nmax
        z(n)=((n-1)/(Nmax-1))*(Zmax-Zmin)+Zmin;
    end
    
    color1=zeros(size(hr_mem));
    Kmax=length(hr_mem);
    
    for k=1:Kmax
        for n=1:Nmax-1
            if (hr_mem(k)>=z(n)) && (hr_mem(k)<z(n+1))
                color1(k)=H1(n);
            end
        end
    end
    
    if j==1
        hr_disp=[hr_disp' data(ind1:ind2)']';
        t_disp=[t_disp' t(ind1:ind2)']';
        c_disp=[c_disp' color1(ind1:ind2)']';
    else
        hr_disp=[hr_disp(1:end-(colorwin-shift+1))' data(ind2-colorwin:ind2)']';
        t_disp=[t_disp(1:end-(colorwin-shift+1))' t(ind2-colorwin:ind2)']';
        if isempty(memorywin)
            c_disp=[c_disp(1:end-(colorwin-shift+1))' color1(ind2-colorwin:ind2)']';
        else
            c_disp=[c_disp(1:end-(colorwin-shift+1))' color1(end-colorwin:end)']';
        end
    end
    
    if videoname~=0
        fig1=figure('Position',...
            [0.05*scrsz(3) 0.2*scrsz(4) 0.89*scrsz(3) 0.60*scrsz(4)],...
            'Color',[0.7 0.7 0.7]);
        scatter(t_disp(ind1:ind2),hr_disp(ind1:ind2),20,...
            c_disp(ind1:ind2),'filled')
        eval(['colormap(flipud(' cmap '(ncolors)));']);
        colorbar('location','EastOutside');
        caxis([0 1])
        xlim([t_disp(ind1) t_disp(ind2)]);
        ylim([lowl hil]);
        ylabel(yunit,'fontsize',24)
        if minsec
            tstart=ceil(t_disp(ind1)/60)*60;
            tend=floor(t_disp(ind2)/60)*60;
            set(gca,'Color',[0.7 0.7 0.7],'FontSize',24,'XTick', [tstart:60:tend],'xticklabel',seconds2time([tstart:60:tend]))
            xlabel('Time (mm:ss)','fontsize',24)
        else
            set(gca,'Color',[0.7 0.7 0.7],'FontSize',24);
            xlabel(timeunit,'fontsize',24)
        end
        grid on
        box on
        
        Fr1=getframe(fig1);
        writeVideo(vidObj,Fr1);
        
        clear Fr1 Im1 Im1_map im_file fr_file;
        close(fig1);
    end
end
if videoname~=0
    close(vidObj);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display whole signal
fig1=figure('Position',...
    [0.05*scrsz(3) 0.05*scrsz(4) 0.89*scrsz(3) 0.89*scrsz(4)],...
    'Color',[1 1 1]);
ah1=subplot(211);
plot(t_disp,hr_disp,'.k','linewidth',4)
ylabel(yunit,'fontsize',24)
xlim([t_disp(1) t_disp(end)])
ylim([lowl hil]);
set(gca,'fontsize',24)
grid on
box on

ah2=subplot(212);
scatter(t_disp,hr_disp,20,c_disp,'filled')
eval(['colormap(flipud(' cmap '(ncolors)));']);
xlim([t_disp(1) t_disp(end)])
ylim([lowl hil]);
set(gca,'fontsize',24)
grid on
box on
ylabel(yunit,'fontsize',24)

xlabel(timeunit);

h=colorbar('location','EastOutside','ytick',[0:0.25:1],'fontsize',24);

%# find current position [x,y,width,height]
pos2 = get(ah2,'Position');
pos1 = get(ah1,'Position');
%# set width of second axes equal to first
pos2(3) = pos1(3);
set(ah2,'Position',pos2)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [array]=seconds2time(seconds)
for i=1:length(seconds)
    mm=floor(seconds(i)/60);
    ss=seconds(i)-mm*60;
    string=[num2str(mm) ':' num2str(ss)];
    if numel(num2str(mm))==1
        string=['0' string];
    end
    if numel(num2str(ss))==1
        string=[string(1:3) '0' string(4)];
    end
    array{i}=string;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Z = smooth1D(Y,lambda)
[m,n] = size(Y);
E = eye(m);
D1 = diff(E,1);
D2 = diff(D1,1);
P = lambda.^2 .* D2'*D2 + 2.*lambda .* D1'*D1;
Z = (E + P) \ Y;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%