function [ Binary,S_n,Q_n] = INADILIC_LCH(Trajectory ,T,Method )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This Matlab function implements the Local Convex Hull method
% by Y. Lanoiselee and D. S. Grebenkov (C).
%
% The function is a part of the non-commercial research project
% INADILIC (see www.inadilic.fr). This function can be freely 
% redistributed and used for non-commercial research purposes, 
% under the condition of keeping credits to Y. Lanoiselee and 
% D. S. Grebenkov and citing the paper, which introduces and 
% describes the method:
%
% Y. Lanoiselee and D. S. Grebenkov, Unravelling intermittent 
% features in single particle trajectories by a local convex 
% hull method, Phys. Rev. E 96, 022144 (2017).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function classifies points of a trajectory in any dimension 
% into two phases by using either the diameter, or the volume of 
% its local convex hull (the volume is realized only for dimension <=3).
%
% INPUT:  
%         If no input => Example mode : Classification of a Brownian motion
%         alterning between two Diffusion coeffients D1=1/2 and D2=2 with
%         phase duration=100. Are plotteds the Score function Sn in 'black', The
%         binary vector of recognition in 'red' and the true binary vector
%         of states in 'blue'. One state is when the binary vectors==0 and
%         the other test is when is does not equal zero.
%         Trajectory - matrix [N x d] containing N successive positions 
%	  of the analyzed trajectory in d dimensions
%         T - the half-size of the window used for the test (T = 10 by default)
%         Method - 'diameter' or 'volume' of the Local Convex Hull ('diameter' by default)
% OUTPUT: Binary - binary vector containing decision from belonging
%         to one phase or one other. The first 2T and the last 2T 
% 	  points of the trajectory remain untreated (the value 2 is set)
%         Sn - Time serie Sn, of averaged contibution of point n to n-T to n+T of Qn
%         Qn - Time serie of functional evaluated along the trajectory 
%	  (diameter or volume of the Local Convex Hull).
%

%% Initialisation
if nargin==2
    Method='diameter';
elseif nargin==1;
    Method='diameter';
    T=10;
elseif nargin==0 
%  Initialisation of Example
N=1000;
sigma=zeros(N,1);
sigma(sin((1:N)./100*pi)>0,1)=1;
sigma(sin((1:N)./100*pi)<=0,1)=2;
Trajectory=cumsum(repmat(sigma,1,2).*randn(N,2),1);
T=10;
Method='diameter';
end;

N=size(Trajectory,1);
Q_n=zeros(N,1);
if strcmp(char(Method),'diameter')==1
    List_couple=nchoosek(1:(2*T+1),2);
end;

for n=T+1:N-T
    if strcmp(char(Method),'diameter')==1
        %% if 'diameter'
        sample=Trajectory(n-T:(n+T),:);
        Q_n(n,1)  = max((sum((sample(List_couple(:,1),:)-sample(List_couple(:,2),:)  ).^2,2)).^0.5);
    elseif strcmp(char(Method),'volume')==1
        %% if 'volume'
        [K,Volume]=  convhull( Trajectory(n-T:(n+T),:) );
        Q_n(n,1)=Volume;
    else
        disp('Error, possible options : ''diameter'', ''volume'' ')
        break;
    end;
end;

%% Calculation of the score fucntion S_n(\tau)
S_n=zeros(N,1);
for i=T+1:N-T
    S_n(i-T:i+T,1)=S_n(i-T:i+T,1)+Q_n(i,1)./(2*T+1);%Time(1,pos_T);%mean(Historique_Hull(i,:));
end;
%% Criterion : active if S_n>quantile(S_n), passive anywhere else
Binary=zeros(N,1);
Mean=mean(S_n(2*T+1:N-2*T,1));
Binary(S_n>Mean)=1;
Binary(S_n<=Mean)=0;
Binary(1:2*T,1)=2;
Binary(end-2*T:end,1)=2;
%% Example : Plot
if nargin==0 
plot((sigma(2:end,1)-min(sigma(2:end,1)))./(max(sigma)-min(sigma(2:end,1))).*max(S_n),'r');
hold on
plot((2*T+1:N-2*T-1),Binary((2*T+1:end-2*T-1),1)*max(S_n),'b')
plot(S_n,'black')
legend('True State','detected state','S_n')
xlabel('t')
ylabel('S_n or State')
end;