function varargout=corrint(varargin)
%
% [y1,y2,y3]=corrint(x,embeddedDim,timeLag,timeStep,distanceThreshold,neighboorSize,estimationMode,findScaling)
%
% Correlation integral analysis of a time series. Based on:
%
% [1] Kaplan, Daniel, and Leon Glass. Understanding nonlinear dynamics. Vol. 19. Springer, 1995.
% [2] Kantz, Holger, and Thomas Schreiber. Nonlinear time series analysis. Cambridge university press, 2004.
%
% Required input parameter:
% x
% Nx1 matrix (doubles) of time series to be analyzed.
%
% Optional Parameters are:
%
% embeddedDim
% 1x1 Integer specifying the embedded dimension size to use (default
% =2).
%
% timeLag
% 1x1 Integer specifying the minimum time lag distance (in samples) of the point to
% be estimated. Default is 2. If timeLag=-1 the timeLag is estimated
% from the first zero-crossing point of the autocorrelation of x.
%
% timeStep
% 1x1 Integer specifying time lag distance (in samples) within
% each point used in the embeddedDimm vector. For example, if embeddedDim
% is 3 and timeStep =2, then the embedded dimension vector will consists of
% 3 samples separated by 2 samples each, covering a window of size of 7 samples.
%
% distanceThreshold
% 1x1 double specifying the distance threshold between embedded
% points. The points who's distance is less than distanceThreshold are considered in
% the same neighborhood and used for either prediction, recurrence, or the
% estimation of the embedded dimension.
%
% neighboorSize
% 1x1 Integer specifying the number of neighbors to be used for
% prediction and smoothing (see 'estimationMode' parameter).
%
% estimationMode
% String specifying what analysis type to be done in the time series.
% Options are:
% 'recurrence' -Calculates recurrence data to be used
% in for recurrence plots (default).
% 'dimension' -Generates statistics for the estimation of the correlation dimension
% of the time series and it's scaling
% regions.
% 'prediction' -Predicts second half of the time
% series using the first half as a model
% and neighboorSize nearest points.
% 'smooth' -Predicts all point of the times
% series using all other points as a
% model and neighboorSize nearest
% points.
%
% findScaling
% 1x1 Boolean flag to be passed when using 'dimension' mode. If set to
% true, the scaling region will be searched automatically, using
% r1=std(x)/4 and r2 -> C(r1)/C(r2) ~ 5. Default value is false.
%
% The output returned by CORRINT is dependent on the 'estimationMode'
% parameter, so that the description of the output below is broken down into the
% different possible options for the 'estimationMode' parameter.
%
% Output Parameters - 'recurrence' mode
% y1
% Lx1 Vector of integers for state i.
% y2
% Lx1 Vector of integers for state state j that is a neighbor of state i (first column).
%
% Output Parameters - 'dimension' mode
% y1
% Lx1 Vector of doubles of log(distanceThreshold).
% y2
% Lx1 Vector of doubles for log(neighborhood size) given the distanceThreshold used in column 1.
% y3
% 1x1 double. Optional, estimated slope of y1 and y2
%
% Output Parameters - 'prediction' mode
% y1
% Lx1 Vector of doubles of estimated second half of the time series.
% y2
% Lx1 Vector of doubles for original second half of the time series.
% y3
% 1x1 double. Optional, variance of the prediction error divided by variance of the second half of the time series.
%
%
% Output Parameters - 'smooth' mode
% y1
% Lx1 Vector of doubles of smoothed the time series.
% y2
% Lx1 Vector of doubles for original time series.
% y3
% 1x1 double. Optional, variance of the prediction error divided by variance of the time series.
%
%
% %%% Beging Example %%%
%
% N=500; %Number of points for each process
% model_names={'linearModel','nonlinearModel'};
%
% %Linear Auto Regressive model with measurement noise
% linearModel=zeros(N,1);
% x=77;
% linearModel(1)=x;
% for n=2:N
% x=4 + 0.95*x;
% linearModel(n)= x + randn(1)*2;
% end
%
% %Non-linear model of dimension ~ 3.9
% nonlinearModel=zeros(N,1);
% x=0.2;y=0.2;z=0.2;v=0.2;model_five(1)=x;
% for n=2:N
% m=0.4 - 6/(1+ x^2 + y^2);
% xold=x;yold=y;zold=z;vold=v;
% x= 1 + 0.7*(xold*cos(m)-yold*sin(m)) + 0.2*zold;
% y=0.7*(xold*sin(m) + yold*cos(m));
% z=1.4 + 0.3*vold - zold^2;
% v=zold;
% nonlinearModel(n)= x + 0.3*z + randn(1)*0.05;
% end
%
% %Plot time series
% figure(1)
% for i=1:2
% subplot(2,1,i)
% eval(['plot(' model_names{i} ');legend(''' model_names{i} ''')'])
% title('Time Plot');xlabel('time')
% end
%
% %Plot cross correlation
% figure(2)
% for i=1:2
% subplot(2,1,i)
% eval(['x=' model_names{i} ';'])
% R=xcorr(x-mean(x),'coeff');
% plot(R(round(N):end))
% eval(['legend(''' model_names{i} ''')'])
% title('Autocorelation'); xlabel('lag')
% end
%
% %Plot Phase Plots
% figure(3)
% for i=1:2
% subplot(2,1,i)
% eval(['x=' model_names{i} ';'])
% scatter(x(1:end-1),x(2:end))
% eval(['legend(''' model_names{i} ''')'])
% title('Phase Plot');xlabel('x(t)');ylabel('x(t+1)')
% end
%
% %Plot prediction errors vs surrogate
% timeLag=1;
% timeStep=1;
% distanceThreshold=[];
% embeddedDim=4;
% estimationMode='smooth';
% figure(4)
% K=[1:20 25 30 50 70 100];
% D=length(K);
% surrN=10;
% for i=1:2
% eval(['x=' model_names{i} ';'])
% err=zeros(D,1)+NaN;
% surr_data=zeros(D,surrN);
% SURR=surrogate(x,surrN);
% for d=1:D;
% neighboorSize=K(d);
% [y1,y2,y3]=corrint(x,embeddedDim,timeLag,timeStep,distanceThreshold,neighboorSize,estimationMode);
% err(d)=y3;
% for s=1:surrN
% [y1,y2,y3]=corrint(SURR(:,s),embeddedDim,timeLag,timeStep,distanceThreshold,neighboorSize,estimationMode);
% surr_data(d,s)=y3;
% end
% end
% subplot(2,1,i)
% plot(K,err,'o-');hold on
% errorbar(K,mean(surr_data,2),var(surr_data,[],2)./sqrt(10),'r')
% eval(['legend(''' model_names{i} ''',''surrogate'')'])
% xlabel('Embedded Dimension')
% ylabel('err/var')
%
% end
%
% %%% End Example %%%
%
% Written by Ikaro Silva, 20134
% Last Modified: November 23, 2014
% Version 1.0
%
% Since 0.9.8
%
%
% See also SURROGATE, DFA, MSENTROPY
%endOfHelp
persistent javaWfdbExec config
if(isempty(javaWfdbExec))
[javaWfdbExec,config]=getWfdbClass('corrint');
end
%Set default pararamter values
inputs={'x','embeddedDim','timeLag','timeStep','distanceThreshold','neighboorSize','estimationMode','findScaling'};
outputs={'y1','y2','y3'};
embeddedDim=[];
timeLag=[];
timeStep=[];
distanceThreshold=[];
neighboorSize=[];
estimationMode='recurrence';
findScaling=0;
wfdb_argument={};
y1=[];
y2=[];
y3=[];
for n=1:nargin
if(~isempty(varargin{n}))
eval([inputs{n} '=varargin{n};'])
end
end
if(~isempty(embeddedDim))
wfdb_argument{end+1}='-d';
wfdb_argument{end+1}=num2str(embeddedDim);
end
if(~isempty(timeLag))
wfdb_argument{end+1}='-t';
wfdb_argument{end+1}=num2str(timeLag);
end
if(~isempty(timeStep))
wfdb_argument{end+1}='-s';
wfdb_argument{end+1}=num2str(timeStep);
end
if(~isempty(distanceThreshold))
wfdb_argument{end+1}='-r';
wfdb_argument{end+1}=num2str(distanceThreshold);
end
if(~isempty(neighboorSize))
wfdb_argument{end+1}='-n';
wfdb_argument{end+1}=num2str(neighboorSize);
end
switch estimationMode
case 'recurrence'
wfdb_argument{end+1}='-p';
case 'dimension'
wfdb_argument{end+1}='-v';
y3=' ';
if(findScaling)
wfdb_argument{end+1}='-a';
end
case 'prediction'
wfdb_argument{end+1}='-P';
case 'smooth'
wfdb_argument{end+1}='-S';
y3=' ';
otherwise
error(['Unkown estimation mode: ' estimationMode])
end
javaWfdbExec.setArguments(wfdb_argument);
if(config.inOctave)
x=cellstr(num2str(x));
x=java2mat(javaWfdbExec.execWithStandardInput(x));
Nx=x.size;
out=cell(Nx,1);
for n=1:Nx
out{n}=x.get(n-1);
end
else
out=cell(javaWfdbExec.execWithStandardInput(x).toArray);
end
M=length(out);
if(~isempty(y3))
y3=out{end};
out(end)=[];
M=M-1;
if(strcmp(estimationMode,'smooth'))
tmp=y3;
sep=regexp(tmp,'\s');
y3=str2num(tmp(sep(end):end));
end
end
if(~isempty(strfind(out{1},'Possibly')))
warning(out{1})
out(1)=[];
M=M-1;
end
y1=zeros(M,1)+NaN;
y2=zeros(M,1)+NaN;
for m=1:M
str=out{m};
sep=regexp(str,'\s');
y1(m)=str2num(str(1:sep));
y2(m)=str2num(str(sep(1):end));
end
for n=1:nargout
eval(['varargout{n}=y' num2str(n) ';'])
end