diff --git a/MAP.m b/MAP.m
new file mode 100644
index 0000000..ef5e9f4
--- /dev/null
+++ b/MAP.m
@@ -0,0 +1,34 @@
+function [klas, Z] = MAP(PostProbs)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% function [klas, Z] = MAP(PostProbs)
+%   
+% calculate a partition by applying the Maximum A Posteriori Bayes
+% allocation rule
+%
+%
+% Inputs : 
+%   PostProbs, a matrix of dimensions [n x K] of the posterior
+%  probabilities of a given sample of n observations arizing from K groups
+%
+% Outputs:
+%   klas: a vector of n class labels (z_1, ...z_n) where z_i =k \in {1,...K}
+%       klas(i) = arg   max (PostProbs(i,k)) , for all i=1,...,n
+%                     1<=k<=K
+%               = arg   max  p(zi=k|xi;theta)
+%                     1<=k<=K
+%               = arg   max  p(zi=k;theta)p(xi|zi=k;theta)/sum{l=1}^{K}p(zi=l;theta) p(xi|zi=l;theta)
+%                     1<=k<=K
+%
+%
+%       Z : Hard partition data matrix [nxK] with binary elements Zik such
+%       that z_ik =1 iff z_i = k
+%
+%%%%%%%%%%%%%%%%%%%%%%%%% Faicel Chamroukhi %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+[n, K] = size(PostProbs);
+% maximum a posteriori rule
+[~,klas] = max(PostProbs,[],2);
+partition_MAP = (klas*ones(1,K))==(ones(n,1)*[1:K]);
+% double   
+Z = double(partition_MAP);
+
diff --git a/forwards_backwards.m b/forwards_backwards.m
new file mode 100644
index 0000000..e8c7144
--- /dev/null
+++ b/forwards_backwards.m
@@ -0,0 +1,72 @@
+function [tau_tk, xi_tkl, alpha_tk, beta_tk, loglik, xi_summed] = forwards_backwards(prior, transmat, f_tk)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% function [tau_tk, xi_tlk, alpha_tk, beta_tk, loglik, xi_summed] = forwards_backwards(prior, transmat, f_tk)
+% forwards_backwards : calculates the E-step of the EM algorithm for an HMM
+% (Gaussian HMM)
+
+% Inputs :
+%
+%         prior(k) = Pr(z_1 = k)
+%         transmat(\ell,k) = Pr(z_t=k | z_{t-1} = \ell)
+%         f_tk(t,k) = Pr(y_t | z_y=k;\theta) %gaussian
+%
+% Outputs:
+%
+%        tau_tk(t,k) = Pr(z_t=k | X): post probs (smoothing probs)
+%        xi_tk\elll(t,k,\ell)  = Pr(z_t=k, z_{t-1}=\ell | Y) t =2,..,n
+%        with Y = (y_1,...,y_n);
+%        alpha_tk(k,t): [Kxn], forwards probs: Pr(y1...yt, zt=k)
+%        beta_tk(k,t): [Kxn], backwards probs: Pr(yt+1...yn|zt=k)
+%
+%
+%
+% Faicel Chamroukhi
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if nargin < 6, filter_only = 0; end
+
+T = size(f_tk, 2);
+K = length(prior);
+
+if size(prior,2)~=1
+    prior = prior';
+end
+scale = ones(1,T);%pour que loglik = sum(log(scale)) part de zero
+
+prior = prior(:);
+tau_tk = zeros(K,T);
+xi_tkl = zeros(K,K,T-1);
+xi_summed = zeros(K,K);
+alpha_tk = zeros(K,T);
+beta_tk = zeros(K,T);
+
+%% forwards: calculate the alpha_tk's
+t = 1;
+alpha_tk(:,1) = prior.* f_tk(:,t);
+[alpha_tk(:,t), scale(t)] = normalize(alpha_tk(:,t));
+
+for t=2:T
+    [alpha_tk(:,t), scale(t)] = normalize((transmat'*alpha_tk(:,t-1)) .* f_tk(:,t));
+    %     filtered_prob (:,:,t-1)= normalize((transmat'*alpha(:,t-1)).*fit(:,t));%
+    % NB: filtered_prob(t,k,l) => filter_prb = squeeze(sum(filters_prob, 2));
+end
+% %loglik (technique du scaling)
+loglik = sum(log(scale));
+
+if filter_only
+    beta_tk = [];
+    tau_tk = alpha_tk;
+    return;
+end
+%% backwards: calculate the beta_tk's, the tau_tk's (and the xi_tkl's)
+%t=T;
+
+beta_tk(:,T) = ones(1,K);
+tau_tk(:,T) = normalize(alpha_tk(:,T) .* beta_tk(:,T));
+
+for t=T-1:-1:1
+    beta_tk(:,t) = normalize(transmat * (beta_tk(:,t+1) .* f_tk(:,t+1)));
+    tau_tk(:,t) = normalize(alpha_tk(:,t) .* beta_tk(:,t));
+    xi_tkl(:,:,t) = normalize((transmat .* (alpha_tk(:,t) * (beta_tk(:,t+1) .* f_tk(:,t+1))')));
+    xi_summed = xi_summed + sum(xi_tkl(:,:,t),3);
+end
+
diff --git a/init_MixFHMM.m b/init_MixFHMM.m
new file mode 100644
index 0000000..0926f21
--- /dev/null
+++ b/init_MixFHMM.m
@@ -0,0 +1,197 @@
+function mixFHMM = init_MixFHMM(Y, K, R, ...
+    variance_type, order_constraint, init_kmeans, try_algo)
+%
+%
+%
+%
+%
+%
+%%%%%%%%%%%%%%%%%%%%%% FC %%%%%%%%%%%%%%
+
+[n, m]=size(Y);
+
+% % 1. Initialization of cluster weights
+mixFHMM.param.w_k=1/K*ones(K,1);
+% Initialization of the model parameters for each cluster
+if init_kmeans
+    max_iter_kmeans = 400;
+    n_tries_kmeans = 20;
+    verbose_kmeans = 0;
+    res_kmeans = myKmeans(Y, K, n_tries_kmeans, max_iter_kmeans,verbose_kmeans);
+    for k=1:K
+        Yk = Y(res_kmeans.klas==k ,:); %if kmeans
+        
+        mixFHMM_init =  init_gauss_hmm(Yk, R, order_constraint, variance_type, try_algo);
+        
+        % 2. Initialisation de \pi_k
+        mixFHMM.param.pi_k(:,k) = mixFHMM_init.initial_prob;%[1;zeros(R-1,1)];
+        
+        % 3. Initialisation de la matrice des transitions
+        mixFHMM.param.A_k(:,:,k)  =  mixFHMM_init.trans_mat;
+        if order_constraint
+            mixFHMM.stats.mask = mixFHMM_init.mask;
+        end
+        
+        % 4. Initialisation des moyennes
+        mixFHMM.param.mu_kr(:,k) = mixFHMM_init.mur;
+        if strcmp(variance_type,'common')
+            mixFHMM.param.sigma_k(k) = mixFHMM_init.sigma;
+        else
+            mixFHMM.param.sigma_kr(:,k) = mixFHMM_init.sigma2r;
+        end
+    end
+else
+    ind = randperm(n);
+    for k=1:K
+        if k<K
+            Yk = Y(ind((k-1)*round(n/K) +1 : k*round(n/K)),:);
+        else
+            Yk = Y(ind((k-1)*round(n/K) +1 : end),:);
+        end
+        mixFHMM_init =  init_gauss_hmm(Yk, R, order_constraint, variance_type, try_algo);
+        
+        % 2. Initialisation de \pi_k
+        mixFHMM.param.pi_k(:,k) = mixFHMM_init.initial_prob;%[1;zeros(R-1,1)];
+        
+        % 3. Initialisation de la matrice des transitions
+        mixFHMM.param.A_k(:,:,k)  =  mixFHMM_init.trans_mat;
+        if order_constraint
+            mixFHMM.mask.mask = mixFHMM_init.mask;
+        end
+        % 4. Initialisation des moyennes
+        mixFHMM.param.mu_kr(:,k) = mixFHMM_init.mur;
+        if strcmp(variance_type,'common')
+            mixFHMM.param.sigma_k(k) = mixFHMM_init.sigma;
+        else
+            mixFHMM.param.sigma_kr(:,k) = mixFHMM_init.sigma2r;
+        end
+    end
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function param =  init_gauss_hmm(Y, R, order_constraint, variance_type, try_EM)
+% init_gauss_hmm  estime les paramètres initiaux d'un hmm où la loi conditionnelle des observations est une gaussienne
+%
+% Entrees :
+%
+%        Y(i,:,nsignal) = x(i) : observation à l'instant i du signal
+%        (séquence) nsignal (notez que pour la partie parametrisation des
+%        signaux les observations sont monodimentionnelles)
+%        R : nbre d'états (classes) cachés
+%
+% Sorties :
+%
+%         model : parametres initiaux du modele. structure
+%         contenant les champs: para: structrure with the fields:
+%         * le HMM initial
+%         1. initial_prob (k) = Pr(Z(1) = k) avec k=1,...,K. loi initiale de z.
+%         2. trans_mat(\ell,k) = Pr(z(i)=k | z(i-1)=\ell) : matrice des transitions
+%         *
+%         3.1. mur : moyenne de l'état k
+%         3.2 sigma2r(k) = variance de x(i) sachant z(i)=k; sigma2r(j) =
+%         sigma^2_r.
+%         mu(:,k) = Esperance de x(i) sachant z(i) = k ;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+if order_constraint
+    % % Tnitialisation en tenant compte de la contrainte:
+    
+    % Initialisation de la matrice des transitions
+    mask = eye(R);%mask d'ordre 1
+    for r=1:R-1
+        ind = find(mask(r,:) ~= 0);
+        mask(r,ind+1) = 1;
+    end
+    % Initialisation de la loi initiale de la variable cachee
+    param.initial_prob = [1;zeros(R-1,1)];
+    param.trans_mat = normalize(mask,2);%
+    param.mask = mask;
+else
+    % Initialisation de la loi initiale de la variable cachee
+    param.initial_prob = 1/R*ones(R,1);
+    param.trans_mat = mk_stochastic(rand(R,R));
+end
+
+%  Initialisation des moyennes et des variances.
+param_gauss = init_gauss_param_hmm(Y, R, variance_type, try_EM);
+
+param.mur = param_gauss.mur;
+if strcmp(variance_type,'common')
+    param.sigma = param_gauss.sigma;
+else
+    param.sigma2r = param_gauss.sigma2r;
+end
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function param = init_gauss_param_hmm(Y, R, variance_type, try_EM)
+% init_regression_model estime les parametres de la loi conditionnelle
+% des observations : une gaussienne d'un hmm homogène d'ordre 1
+%
+% Entrees :
+%
+%        Y : [nxm]
+%        nsignal (notez que pour la partie parametrisation des signaux les
+%        observations sont monodimentionnelles)
+%        R : nbre d'états (classes) cachés
+% Sorties :
+%
+%
+%         para : parametres initiaux de la loi cond de chaque état
+%         2. sigma2r(r) = variance de y(t) sachant z(t)=r; sigmar(j) =
+%         sigma^2_r.
+%         3. mu(:,r) : E[y(t)|z(t) =r] ;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+[n, m] = size(Y);
+
+if strcmp(variance_type,'common'),  s=0; end
+
+if (try_EM) ==1
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %decoupage de l'echantillon (signal) en K segments
+    zi = round(m/R)-1;
+    for r=1:R
+        i = (r-1)*zi+1;
+        j = r*zi;
+        
+        Yij = Y(:,i:j);
+        Yij = reshape(Yij',[],1);
+        param.mur(r) = mean(Yij);
+        if strcmp(variance_type,'common')
+            s=s+ sum((Yij-param.mur(r)).^2);
+            param.sigma = s/(n*m);
+        else
+            m_r = j-i+1 ;
+            param.sigma2r(r) = sum((Yij-param.mur(r)).^2)/(n*m_r);
+        end
+    end
+else % initialisation aléatoire
+    Lmin= 2;%round(m/(K+1));%nbr pts min dans un segments
+    tr_init = zeros(1,R+1);
+    tr_init(1) = 0;
+    R_1=R;
+    for r = 2:R
+        R_1 = R_1-1;
+        temp = tr_init(r-1)+Lmin:m-R_1*Lmin;
+        ind = randperm(length(temp));
+        tr_init(r)= temp(ind(1));
+    end
+    tr_init(R+1) = m;
+    for r=1:R
+        i = tr_init(r)+1;
+        j = tr_init(r+1);
+        Yij = Y(:,i:j);
+        Yij = reshape(Yij',[],1);
+        
+        param.mur(r) = mean(Yij);
+        
+        if strcmp(variance_type,'common')
+            s=s+ sum((Yij-param.mur(r)).^2);
+            param.sigma = s/(n*m);
+        else
+            m_r = j-i+1 ;
+            param.sigma2r(r) = sum((Yij-param.mur(r)).^2)/(n*m_r);
+        end
+    end
+end
diff --git a/learn_MixFHMM.m b/learn_MixFHMM.m
new file mode 100644
index 0000000..3e4c9ac
--- /dev/null
+++ b/learn_MixFHMM.m
@@ -0,0 +1,302 @@
+function mixFHMM =  learn_MixFHMM(data, K, R, ...
+    variance_type, order_constraint, total_EM_tries, max_iter_EM, init_kmeans, threshold, verbose)
+%
+% The EM algorithm for parameter estimation of the mixture of Hidden Markov
+% Models for clustering and segmentation of time series with regime changes
+%
+% Inputs:
+% data: a set of n time series with m observations (dim: [n x m]
+% K: number of clusters
+% R: number of regimes (states)
+% options
+%
+%
+%
+% faicel chamroukhi (septembre 2009)
+%
+%% Please cite the following references for this code
+%
+% @InProceedings{Chamroukhi-IJCNN-2011,
+%   author = {F. Chamroukhi and A. Sam\'e  and P. Aknin and G. Govaert},
+%   title = {Model-based clustering with Hidden Markov Model regression for time series with regime changes},
+%   Booktitle = {Proceedings of the International Joint Conference on Neural Networks (IJCNN), IEEE},
+%   Pages = {2814--2821},
+%   Adress = {San Jose, California, USA},
+%   year = {2011},
+%   month = {Jul-Aug},
+%   url = {https://chamroukhi.com/papers/Chamroukhi-ijcnn-2011.pdf}
+% }
+%
+% @PhdThesis{Chamroukhi_PhD_2010,
+% author = {Chamroukhi, F.},
+% title = {Hidden process regression for curve modeling, classification and tracking},
+% school = {Universit\'e de Technologie de Compi\`egne},
+% month = {13 december},
+% year = {2010},
+% type = {Ph.D. Thesis},
+% url ={https://chamroukhi.com/papers/FChamroukhi-Thesis.pdf}
+% }
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+warning off
+[n, m] = size(data);%n  nbre de signaux (individus); m: nbre de points pour chaque signal
+%
+Y=reshape(data',[],1);
+
+%
+try_EM = 0;
+best_loglik = -inf;
+cputime_total = [];
+
+while try_EM < total_EM_tries
+    try_EM = try_EM +1;
+    fprintf('EM try n° %d\n',try_EM);
+    time = cputime;
+    %%%%%%%%%%%%%%%%%%%
+    %  Initialization %
+    %%%%%%%%%%%%%%%%%%%
+    mixFHMM = init_MixFHMM(data, K, R,...
+        variance_type, order_constraint, init_kmeans, try_EM);
+    
+    %     Psi = zeros(nu,1);% vecteur parametre
+    iter = 0;
+    converge = 0;
+    loglik = 0;
+    prev_loglik=-inf;
+    stored_loglik=[];
+    
+    % main algorithm
+    % % EM %%%%
+    while ~converge && (iter< max_iter_EM)
+        
+        %
+        exp_num_trans_ck  = zeros(R,R,n);
+        exp_num_trans_from_l_ck = zeros(R,n);
+        %
+        exp_num_trans = zeros(R,R,n,K);
+        exp_num_trans_from_l = zeros(R,n,K);
+        %
+        w_k_fyi = zeros(n,K);
+        log_w_k_fyi = zeros(n,K);
+        
+        %%%%%%%%%%
+        % E-Step %
+        %%%%%%%%%%
+        gamma_ikjr = zeros(n*m,R,K);
+        for k=1:K
+            % run a hmm for each sequence
+            log_fkr_yij =zeros(R,m);
+            %
+            Li = zeros(n,1);% to store the loglik for each example
+            %
+            mu_kr = mixFHMM.param.mu_kr(:,k);
+            
+            for i=1:n
+                y_i = data(i,:);
+                for r = 1:R
+                    mukr = mu_kr(r);
+                    %sk = sigma_kr(k);
+                    if strcmp(variance_type,'common')
+                        sigma_kr = mixFHMM.param.sigma_k(k);
+                        sk = sigma_kr;
+                    else
+                        sigma_kr = mixFHMM.param.sigma_kr(:,k);
+                        sk = sigma_kr(r);
+                    end
+                    z=((y_i-mukr*ones(1,m)).^2)/sk;
+                    log_fkr_yij(r,:) = -0.5*ones(1,m).*(log(2*pi)+log(sk)) - 0.5*z;% pdf cond à  c_i = g et z_i = k de yij
+                    fkr_yij(r,:) = normpdf(y_i, mukr*ones(1,m), sqrt(sk));
+                end
+                %                         log_fkr_yij  = min(log_fkr_yij,log(realmax));
+                %                         log_fkr_yij = max(log_fkr_yij ,log(realmin));
+                %                         fkr_yij =  exp(log_fkr_yij);
+                
+                % calcul de p(y) : forwards backwards
+                
+                [gamma_ik, xi_ik, fwd_ik, backw_ik, loglik_i] = forwards_backwards(mixFHMM.param.pi_k(:,k), mixFHMM.param.A_k(:,:,k), fkr_yij);
+                %
+                
+                
+                Li(i) = loglik_i; % loglik of the ith curve
+                %
+                gamma_ikjr((i-1)*m+1:i*m,:,k) = gamma_ik';%[n*m K G]
+                % xi_ikjrl(:,:,(i-1)*(m-1)+1:i*(m-1),g) =  xi_ik;%[KxK n*m G]
+                %
+                exp_num_trans_ck(:,:,i) = sum(xi_ik,3); % [K K n]
+                exp_num_trans_from_l_ck(:,i) = gamma_ik(:,1);%[K x n]
+                %
+            end
+            
+            exp_num_trans_from_l(:,:,k) = exp_num_trans_from_l_ck;%[K n G]
+            exp_num_trans(:,:,:,k) = exp_num_trans_ck;%[K K n G]
+            
+            % for the MAP partition:  the numerator of the cluster post
+            % probabilities
+            num_log_post_prob(:,k) = log(mixFHMM.param.w_k(k)) + Li;
+            
+            % for computing the global loglik
+            w_k_fyi(:,k) = mixFHMM.param.w_k(k)*exp(Li);%[nx1]
+            
+            log_w_k_fyi(:,k) = log(mixFHMM.param.w_k(k)) + Li;
+        end
+        
+        log_w_k_fyi = min(log_w_k_fyi,log(realmax));
+        log_w_k_fyi = max(log_w_k_fyi,log(realmin));
+        
+        tau_ik = exp(log_w_k_fyi)./(sum(exp(log_w_k_fyi),2)*ones(1,K));
+        
+        % % log-likelihood
+        loglik = sum(log(sum(exp(log_w_k_fyi),2)));
+        
+        
+        %%%%%%%%%%
+        % M-Step %
+        %%%%%%%%%%
+        
+        % Maximization of Q1 w.r.t w_k
+        mixFHMM.param.w_k = sum(tau_ik,1)'/n;
+        for k=1:K
+            
+            if strcmp(variance_type,'common'), s=0; end
+            
+            weights_cluster_k = tau_ik(:,k);
+            % Maximization of Q2 w.r.t \pi^g
+            exp_num_trans_k_from_l =   (ones(R,1)*weights_cluster_k').*exp_num_trans_from_l(:,:,k);%[K x n]
+            mixFHMM.param.pi_k(:,k) = (1/sum(tau_ik(:,k)))*sum(exp_num_trans_k_from_l,2);% sum over i
+            % Maximization of Q3 w.r.t A^g
+            for r=1:R
+                if n==1
+                    exp_num_trans_k(r,:,:) = (ones(R,1)*weights_cluster_k)'.*squeeze(exp_num_trans(r,:,:,k));
+                else
+                    %exp_num_trans_k(k,:,:,g)
+                    exp_num_trans_k(r,:,:) = (ones(R,1)*weights_cluster_k').*squeeze(exp_num_trans(r,:,:,k));
+                end
+            end
+            if n==1
+                temp = exp_num_trans_k;
+            else
+                temp = sum(exp_num_trans_k,3);%sum over i
+            end
+            mixFHMM.param.A_k(:,:,k) = mk_stochastic(temp);
+            % if HMM with order constraints
+            if order_constraint
+                mixFHMM.param.A_k(:,:,k) = mk_stochastic(mixFHMM.stats.mask.*mixFHMM.param.A_k(:,:,k));
+            end
+            
+            % Maximisation de Q4 par rapport aux muk et sigmak
+            % each sequence i (m observations) is first weighted by the cluster weights
+            weights_cluster_k =  repmat((tau_ik(:,k))',m,1);
+            weights_cluster_k = weights_cluster_k(:);
+            
+            % secondly, the m observations of each sequance are weighted by the
+            % wights of each segment k (post prob of the segments for each
+            % cluster g)
+            gamma_ijk = gamma_ikjr(:,:,k);% [n*m K]
+            
+            nm_kr=sum(gamma_ijk,1);% cardinal nbr of the segments k,k=1,...,K within each cluster g, at iteration q
+            
+            sigma_kr = zeros(R,1);
+            for r=1:R
+                nmkr = nm_kr(r);%cardinal nbr of segment k for the cluster g
+                % % Maximization w.r.t muk
+                weights_seg_k = gamma_ijk(:,r);
+                
+                mu_kr(r) = (1/sum(weights_cluster_k.*weights_seg_k))*sum((weights_cluster_k.*weights_seg_k).*Y);
+                % % Maximization w.r.t sigmak :
+                z = sqrt(weights_cluster_k.*weights_seg_k).*(Y-ones(n*m,1)*mu_kr(r));
+                if strcmp(variance_type,'common')
+                    s = s + z'*z;
+                    ngm = sum(sum((weights_cluster_k*ones(1,R)).*gamma_ijk));
+                    sigma_k = s/ngm;
+                else
+                    ngmk = sum(weights_cluster_k.*weights_seg_k);
+                    sigma_kr(r)=  z'*z/(ngmk);
+                end
+            end
+            mixFHMM.param.mu_kr(:,k) = mu_kr;
+            if strcmp(variance_type,'common')
+                mixFHMM.param.sigma_k(k) = sigma_k;
+            else
+                mixFHMM.param.sigma_kr(:,k) = sigma_kr;
+            end
+        end
+        iter=iter+1;
+        
+        if prev_loglik-loglik > 1e-3, fprintf(1, '!!!!! EM log-lik is decreasing from %6.4f to %6.4f!\n', prev_loglik, loglik);end
+        if verbose
+            fprintf(1,'EM : Iteration : %d log-likelihood : %f \n',  iter,loglik);
+        end
+        
+        converge =  abs((loglik-prev_loglik)/prev_loglik) <= threshold;
+        prev_loglik = loglik;
+        stored_loglik = [stored_loglik loglik];
+        
+    end % end of EM  loop
+    cputime_total = [cputime_total cputime-time];
+    
+    mixFHMM.param = mixFHMM.param;
+    if strcmp(variance_type,'common')
+        mixFHMM.stats.Psi = [mixFHMM.param.w_k(:); mixFHMM.param.A_k(:); mixFHMM.param.pi_k(:); mixFHMM.param.mu_kr(:); mixFHMM.param.sigma_k(:)];
+    else
+        mixFHMM.stats.Psi = [mixFHMM.param.w_k(:); mixFHMM.param.A_k(:); mixFHMM.param.pi_k(:); mixFHMM.param.mu_kr(:); mixFHMM.param.sigma_kr(:)];
+    end
+    mixFHMM.stats.tau_ik = tau_ik;
+    mixFHMM.stats.gamma_ikjr = gamma_ikjr;
+    mixFHMM.stats.loglik = loglik;
+    mixFHMM.stats.stored_loglik = stored_loglik;
+    mixFHMM.stats.log_w_k_fyi = log_w_k_fyi;
+    
+    if mixFHMM.stats.loglik > best_loglik
+        best_loglik = mixFHMM.stats.loglik;
+        best_mixFHMM.stats = mixFHMM.stats;
+    end
+    
+    if try_EM>=1,  fprintf('log-lik at convergence: %f \n', mixFHMM.stats.loglik); end
+    
+end
+mixFHMM.stats.loglik = best_loglik;
+%
+if try_EM>1,  fprintf('log-lik max: %f \n', mixFHMM.stats.loglik); end
+
+mixFHMM.stats = best_mixFHMM.stats;
+% Finding the curve partition by using the MAP rule
+[klas, Cik] = MAP(mixFHMM.stats.tau_ik);% MAP partition of the n sequences
+mixFHMM.stats.klas = klas;
+
+% cas ou on prend la moyenne des gamma ijkr
+smoothed = zeros(m,K);
+mean_curves = zeros(m,R,K);
+%mean_gamma_ijk = zeros(m,R,K);
+for k=1:K
+    weighted_segments = sum(gamma_ikjr(:,:,k).*(Y*ones(1,R)),2);
+    %weighted_segments = sum(gamma_ikjr(:,:,g).*(ones(n*m,1)*mixFHMM.param.mu_kr(:,k)'),2);
+    
+    %
+    weighted_segments = reshape(weighted_segments,m,n);
+    weighted_clusters = (ones(m,1)*mixFHMM.stats.tau_ik(:,k)').* weighted_segments;
+    smoothed(:,k) = (1/sum(mixFHMM.stats.tau_ik(:,k)))*sum(weighted_clusters,2);
+end
+
+mixFHMM.stats.smoothed = smoothed;
+mixFHMM.stats.mean_curves = mean_curves;
+%mixFHMM.stats.mean_gamma_ijk = mean_gamma_ijk;
+
+mixFHMM.stats.cputime = mean(cputime_total);
+
+% % Segmentation of each cluster using the MAP rule
+% for k=1:K
+%     [segments_k Zjk] = MAP(mixFHMM.stats.mean_gamma_ijk(:,:,k));%MAP segmentation of each cluster of sequences
+%     mixFHMM.stats.segments(:,k) = segments_k;
+% end
+
+nu = length(mixFHMM.stats.Psi);
+% BIC AIC et ICL*
+mixFHMM.stats.BIC = mixFHMM.stats.loglik - (nu*log(n)/2);%n*m/2!
+mixFHMM.stats.AIC = mixFHMM.stats.loglik - nu;
+% ICL*
+% Compute the comp-log-lik
+cik_log_w_k_fyi = (Cik).*(mixFHMM.stats.log_w_k_fyi);
+comp_loglik = sum(sum(cik_log_w_k_fyi,2));
+mixFHMM.stats.ICL1 = comp_loglik - nu*log(n)/2;%n*m/2!
+
+
diff --git a/main_MixFHMM.m b/main_MixFHMM.m
new file mode 100644
index 0000000..d5f34ea
--- /dev/null
+++ b/main_MixFHMM.m
@@ -0,0 +1,77 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Clustering and segmentation of time series (including with regime
+% changes) by mixture of gaussian Hidden Markov Models (MixFHMMs) and the EM algorithm; i.e functional data
+% clustering and segmentation
+%
+%
+%
+%
+% by Faicel Chamroukhi, 2009
+%
+%% Please cite the following references for this code
+% 
+% @InProceedings{Chamroukhi-IJCNN-2011,
+%   author = {F. Chamroukhi and A. Sam\'e  and P. Aknin and G. Govaert},
+%   title = {Model-based clustering with Hidden Markov Model regression for time series with regime changes},
+%   Booktitle = {Proceedings of the International Joint Conference on Neural Networks (IJCNN), IEEE},
+%   Pages = {2814--2821},
+%   Adress = {San Jose, California, USA},
+%   year = {2011},
+%   month = {Jul-Aug},
+%   url = {https://chamroukhi.com/papers/Chamroukhi-ijcnn-2011.pdf}
+% }
+% 
+% @PhdThesis{Chamroukhi_PhD_2010,
+% author = {Chamroukhi, F.},
+% title = {Hidden process regression for curve modeling, classification and tracking},
+% school = {Universit\'e de Technologie de Compi\`egne},
+% month = {13 december},
+% year = {2010},
+% type = {Ph.D. Thesis},
+% url ={https://chamroukhi.com/papers/FChamroukhi-Thesis.pdf}
+% }
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+clear 
+close all;
+clc;
+
+
+%% Model specification
+K = 3;% number of clusters
+R = 3;% number of regimes (HMM states)
+
+% % options
+variance_type = 'common';
+%variance_type = 'free'; 
+ordered_states = 1;
+total_EM_tries = 1;
+max_iter_EM = 1000;
+init_kmeans = 1;
+threshold = 1e-6;
+verbose = 1; 
+
+
+%% toy time series with regime changes
+load simulated_data.mat
+Y; % 
+% [n, m] = size(Y);
+% n: nbr of time sries
+% m: number of observations
+
+
+%%
+mixFHMM =  learn_MixFHMM(Y, K , R, ...
+    variance_type, ordered_states, total_EM_tries, max_iter_EM, init_kmeans, threshold, verbose);        
+ 
+%%
+show_mixFHMM_results(Y, mixFHMM)
+
+
+ 
+
+
+
+
+
+
diff --git a/mk_stochastic.m b/mk_stochastic.m
new file mode 100644
index 0000000..59f3728
--- /dev/null
+++ b/mk_stochastic.m
@@ -0,0 +1,27 @@
+function T = mk_stochastic(T)
+% MK_STOCHASTIC Ensure the argument is a stochastic matrix, i.e., the sum over the last dimension is 1.
+% T = mk_stochastic(T)
+%
+% If T is a vector, it will sum to 1.
+% If T is a matrix, each row will sum to 1.
+% If T is a 3D array, then sum_k T(i,j,k) = 1 for all i,j.
+
+% Set zeros to 1 before dividing
+% This is valid since S(j) = 0 iff T(i,j) = 0 for all j
+
+if (ndims(T)==2) & (size(T,1)==1 | size(T,2)==1) % isvector
+  T = normalise(T);
+elseif ndims(T)==2 % matrix
+  S = sum(T,2); 
+  S = S + (S==0);
+  norm = repmat(S, 1, size(T,2));
+  T = T ./ norm;
+else % multi-dimensional array
+  ns = size(T);
+  T = reshape(T, prod(ns(1:end-1)), ns(end));
+  S = sum(T,2);
+  S = S + (S==0);
+  norm = repmat(S, 1, ns(end));
+  T = T ./ norm;
+  T = reshape(T, ns);
+end
diff --git a/myKmeans.m b/myKmeans.m
new file mode 100644
index 0000000..4d9806c
--- /dev/null
+++ b/myKmeans.m
@@ -0,0 +1,100 @@
+function solution = myKmeans(X, K, nbr_runs, nbr_iter_max, verbose)
+%   function res = myKmeans(X, K, nbr_runs, nbr_iter_max, verbose)
+%
+%   Algorithme des K-means
+%
+%
+%
+% Faicel CHAMROUKHI Septembre 2008 (mise a jour)
+%
+%
+%
+% distance euclidienne
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if nargin<5; verbose=0;end
+if nargin<4; nbr_iter_max = 300;end
+if nargin<3; nbr_runs = 20;end
+
+
+[n, p] = size(X);
+% if one class
+global_mean = mean(X,1);
+
+if K==1
+    dmin = sum((X-ones(n, 1)*global_mean).^2,2);
+    solution.muk = global_mean;
+    klas = ones(n,1);
+    solution.klas = klas;
+    solution.err = sum(dmin);
+    solution.Zik = ones(n,1);
+    return;
+end
+
+
+nbr_run = 0;
+best_solution.err = inf;
+while (nbr_run<nbr_runs)
+    nbr_run=nbr_run + 1;
+    if (nbr_runs>1 && verbose); fprintf('Kmeans run n° : %d  \n',nbr_run);end
+    
+    iter = 0;
+    converged = 0;
+    previous_err = -inf;
+    Zik = zeros(n,K);% partition
+    %% 1. Initialization of the centres
+    rnd_indx = randperm(n);
+    centres = X(rnd_indx(1:K),:);
+    while (iter<nbr_iter_max && ~converged)
+        iter = iter+1;
+        old_centres = centres;
+        
+        % The Euclidean distances
+        eucld_dist = zeros(n, K);
+        for k = 1:K
+            muk = centres(k,:);
+            eucld_dist(:,k) = sum((X-ones(n,1)*muk).^2,2);
+        end
+        %% classification step
+        
+        [dmin, klas] = min(eucld_dist,[],2);
+        
+        Zik = (klas*ones(1,K))==(ones(n,1)*[1:K]);
+        %% relocation step
+        for k=1:K
+            ind_ck = find(klas==k);
+            %if empty classes
+            if isempty(ind_ck)
+                centres(k,:)= old_centres(k,:);
+            else
+                % update the centres
+                centres(k,:) = mean(X(ind_ck,:),1);
+            end
+        end
+        
+        % test of convergence
+        current_err = sum(sum(Zik.*eucld_dist,2)); % the distorsion measure
+
+        if (abs(current_err-previous_err))/previous_err <1e-6
+            converged = 1;
+        end
+        previous_err = current_err;
+        if verbose
+            fprintf('Kmeans : Iteration  %d  Objective: %6f  \n', iter, current_err);
+        end
+        solution.stored_err(iter) = current_err;
+    end% one run
+    %%
+    solution.muk = centres;
+    solution.Zik = Zik;
+    solution.klas = klas;
+    solution.err = current_err;
+    %
+    if current_err < best_solution.err
+        best_solution = solution;
+    end
+end %en of the Kmeans runs
+
+solution = best_solution;
+
+
diff --git a/normalize.m b/normalize.m
new file mode 100644
index 0000000..ffac04b
--- /dev/null
+++ b/normalize.m
@@ -0,0 +1,34 @@
+function [M, z] = normalize(A, dim)
+% NORMALISE Make the entries of a (multidimensional) array sum to 1
+% [M, c] = normalise(A)
+% c is the normalizing constant
+%
+% [M, c] = normalise(A, dim)
+% If dim is specified, we normalise the specified dimension only,
+% otherwise we normalise the whole array.
+
+if nargin < 2
+  z = sum(A(:));
+  % Set any zeros to one before dividing
+  % This is valid, since c=0 => all i. A(i)=0 => the answer should be 0/1=0
+  s = z + (z==0);
+  M = A / s;
+elseif dim==1 % normalize each column
+  z = sum(A);
+  s = z + (z==0);
+  %M = A ./ (d'*ones(1,size(A,1)))';
+  M = A ./ repmat(s, size(A,1), 1);
+else
+  % Keith Battocchi - v. slow because of repmat
+  z=sum(A,dim);
+  s = z + (z==0);
+  L=size(A,dim);
+  d=length(size(A));
+  v=ones(d,1);
+  v(dim)=L;
+  %c=repmat(s,v);
+  c=repmat(s,v');
+  M=A./c;
+end
+
+
diff --git a/show_mixFHMM_results.m b/show_mixFHMM_results.m
new file mode 100644
index 0000000..a2e93bf
--- /dev/null
+++ b/show_mixFHMM_results.m
@@ -0,0 +1,50 @@
+function show_mixFHMM_results(data,mixFHMM)
+%
+%
+%
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+set(0,'defaultaxesfontsize',14);
+
+
+% symboles = {'b--o','b:s','b->','b.-x', '+', 'o', '*' ,'p', 'h' ,'o', '<', '^' ,'v'};
+% couleurs = {'b--o', 'r:s','g->','m-s','black-s'};
+% %
+[n, m]=size(data);
+
+t = 0:m-1;
+scrsz = get(0,'ScreenSize');
+figure('Position',[10 scrsz(4)/2 550 scrsz(4)/2.15]);
+plot(t,data')
+xlabel('t')
+ylabel('y(t)')
+title('original time series')
+
+K = length(mixFHMM.param.w_k);
+
+colors = {'r','g','b','k','m','y','c','r','g','b','k','m','y','c'};
+figure('Position',[scrsz(4) scrsz(4)/2 550 scrsz(4)/2.15]);
+%clf(1)
+
+for k=1:K
+    plot(t,data(mixFHMM.stats.klas==k,:)','color',colors{k})
+    % hold on, plot(t,solution.smoothed(:,k),'color',colors{K+k},'linewidth', 2.5);
+    hold on
+end
+xlabel('t')
+ylabel('y(t)')
+title('Clustered time series')
+
+K = size(mixFHMM.stats.tau_ik,2);
+for k=1:K
+    %subplot(G,1,g),plot(solution.tau_ijgk(:,:,g))
+    %title(['\tau_{ ijgk} , g = ',int2str(g)])
+    %     hold on, subplot(G,1,g),
+    figure,plot(t,data(mixFHMM.stats.klas==k,:)','color',colors{k})
+    hold on, plot(t,mixFHMM.stats.smoothed(:,k),'color',colors{K+k},'linewidth',2.5);
+    xlabel('t')
+ylabel('y(t)')
+title(['Cluster ',int2str(k)])
+end
+
+
diff --git a/simulated_data.mat b/simulated_data.mat
new file mode 100644
index 0000000..6cc8017
Binary files /dev/null and b/simulated_data.mat differ