/*
 * Créé le 10 juin 2004
 *
 */
package traitement.moteur;

import information.DomaineIncompatibleException;
import information.magnitude.Magnitude;
import information.tpt.TPT;
import information.type.Type;

import jama.EigenvalueDecomposition;
import jama.Matrix;

import java.io.IOException;
import java.util.Vector;

import principal.copiercoller.HomonymeIntrouvableException;
import principal.copiercoller.InterfaceHomonyme;

import traitement.TraitementEvent;
import traitement.TraitementIncompletException;
import traitement.TraitementListener;
import traitement.modele.Modele;
import traitement.noeud.InterfaceNoeud;
import traitement.noeud.fonctionnel.NoeudFonctionnelEndogene;
import traitement.noeud.fonctionnel.NoeudFonctionnelExogene;
import traitement.noeud.fonctionnel.VectorNoeudFonctionnel;
import fichier.ContenuFichierIncompatibleException;
import fonction.magnitude.MagnitudeBruit;
import fonction.magnitude.VectorFonctionFilsMagnitude;
import gui.InterfaceAffichable;

/**
 * adaptation du filtre ŕ différences divisées du second ordre (DD2) de Magnus Norgaard.
 * 
 * @author	Vincent Labatut
 * @version	1
 * @see		traitement.modele.Modele
  */
public class MoteurDD2 implements InterfaceMoteur
{	
//	----------------------------------------	
//	Moteur
//	----------------------------------------
	/**
	* Liste des TraitementListeners écoutant cette Moteur.
	* 
	 */
	private Vector vectorTraitementListener;
	/**
	* Modele parent de ce Moteur.
	* 
	 */
	private Modele parent;
	/**
	 * indique si cet objet a été modifé depuis sa derničre sauvegarde.
	  */
	private boolean modifie;
//	----------------------------------------	
//	Listes de Noeuds
//	----------------------------------------
	/**
	* Liste des NoeudExogenes utilisée pendant le traitement.
	* 
	 */
	private VectorNoeudFonctionnel vectorNoeudExogene;
	/**
	* Liste des NoeudEndogenes non observés pour l'activation, utilisée pendant le traitement.
	* 
	 */
	private VectorNoeudFonctionnel vectorNoeudEndogeneActivation;
	/**
	* Liste des NoeudEndogenes non observés pour l'émission, utilisée pendant le traitement.
	* 
	 */
	private VectorNoeudFonctionnel vectorNoeudEndogeneEmission;
	/**
	* Liste des NoeudEndogenes observés pour l'activation, utilisée pendant le traitement.
	* 
	 */
	private VectorNoeudFonctionnel vectorNoeudObserveActivation;
	/**
	* Liste des NoeudEndogenes observés pour l'émission, utilisée pendant le traitement.
	* 
	 */
	private VectorNoeudFonctionnel vectorNoeudObserveEmission;
//	----------------------------------------	
//	Listes de bruits
//	----------------------------------------
	/**
	* Liste des bruits de traitement pour l'activation, utilisée pendant le traitement.
	* 
	 */
	private VectorFonctionFilsMagnitude vectorBruitTraitementActivation;
	/**
	* Liste des bruits de traitement pour l'émission, utilisée pendant le traitement.
	* 
	 */
	private VectorFonctionFilsMagnitude vectorBruitTraitementEmission;
	/**
	* Liste des bruits de mesure pour l'activation, utilisée pendant le traitement.
	* 
	 */
	private VectorFonctionFilsMagnitude vectorBruitMesureActivation;
	/**
	* Liste des bruits de mesure pour l'émission, utilisée pendant le traitement.
	* 
	 */
	private VectorFonctionFilsMagnitude vectorBruitMesureEmission;
	
//	----------------------------------------	
//	Dimensions
//	----------------------------------------
	/**
	* dimension de la variable d'entrée (ie nbre de NoeudExogenes) 
	* 
	 */
	private int nu; 
	/**
	* dimension de la variable d'état (ie nbre de NoeudEndogenes non-observés) 
	* 
	 */
	private int nxA; 
	/**
	* dimension de la variable d'état (ie nbre de NoeudEndogenes non-observés) 
	* 
	 */
	private int nxE; 
	/**
	* dimension de la variable de sortie ou observée (ie nbre de NoeudEndogenes observés) 
	* 
	 */
	private int nyA;
	/**
	* dimension de la variable de sortie ou observée (ie nbre de NoeudEndogenes observés) 
	* 
	 */
	private int nyE;
	/**
	* dimension de la matrice de covariance de bruit (ie nbre de sources de bruit dans les NoeudEndogenes non-observés) 
	* 
	 */
	private int nvA;
	/**
	* dimension de la matrice de covariance de bruit (ie nbre de sources de bruit dans les NoeudEndogenes non-observés)
	* 
	 */
	private int nvE; 
	/**
	* dimension de la matrice de covariance de bruit observée (ie nbre de sources de bruit dans les NoeudEndogenes observés)
	* 
	 */
	private int nwA;
	/**
	* dimension de la matrice de covariance de bruit observée (ie nbre de sources de bruit dans les NoeudEndogenes observés)
	* 
	 */
	private int nwE; // 

//	----------------------------------------	
//	Divers
//	----------------------------------------
	/**
	* Squared divided-difference step size
	* 
	 */
	private double h2;
	/**
	* Divided difference step-size
	* 
	 */
	private double h;
	/**
	* A convenient scaling factor
	* 
	 */
	private double scal1;
	/**
	* Another scaling factor
	* 
	 */
	private double scal2;
	
//	----------------------------------------	
//	Matrices
//	----------------------------------------
	/**
	* Matrice d'état a posteriori
	* 
	 */
	private Matrix xhat;						
	/**
	* Matrice d'état a priori
	* 
	 */
	private Matrix xbar;						
	/**
	* Matrice de covariance triangulaire
	* 
	 */
	private Matrix Sx;
	/**
	* Matrice de covariance triangulaire
	* 
	 */
	private Matrix Sxbar;
	/**
	* Matrice de covariance triangulaire
	* 
	 */
	private Matrix SxxSxv;
	/**
	* Matrice de covariance triangulaire
	* 
	 */
	private Matrix SyxSyw;
	/**
	* Matrice de bruits 
	* 
	 */
	private Matrix wmean;
	/**
	* Matrice de bruits 
	* 
	 */
	private Matrix vmean;
	/**
	* Matrice de bruits 
	* 
	 */
	private Matrix hSv;
	/**
	* Matrice de bruits 
	* 
	 */
	private Matrix hSw;
	/**
	* Matrice d'historique
	* 
	 */
	private Matrix Smat;
	/**
	* Matrice d'historique
	* 
	 */
	private Matrix sidx; 
	/**
	* indices
	* 
	 */
	private int nxnv2;
	/**
	* indices
	* 
	 */
	private int nxnw2;
	/**
	* indices
	* 
	 */
	private int idx_syw2;
	/**
	* indices
	* 
	 */
	private int idx_sxv2;
		
//	----------------------------------------	
//	Constructeurs
//	----------------------------------------
	/**
	 * Constructeur pour créer un Moteur
	 * de parent le Modele
	 * passé en paramčtre.
	 * 
	 * @param	m	parent de ce Moteur.
	 * 
	  */
	public MoteurDD2(Modele m)
	{	parent = m;
		vectorTraitementListener = new Vector();
		setModifie();
   	}

//	----------------------------------------	
//	Parent
//	----------------------------------------
	/*
	 *
	 */
	public Modele getParent()
	{	return parent;
	}
	/*
	 *
	 */
	public void setParent(Modele p)
	{	parent = p;
	}	
	
//	----------------------------------------	
//	Modele
//	----------------------------------------
	/*
	 *
	 */
	public Modele getModele()
	{	return parent;
	}
//	----------------------------------------	
//	InterfaceTraitement
//	----------------------------------------
	/*
	 *
	 */
	public void addTraitementListener(TraitementListener i)
	{	vectorTraitementListener.add(i);
	}
	/*
	 *
	 */
	public void removeTraitementListener(TraitementListener i)
	{	vectorTraitementListener.remove(i);
	}
	/**
	 * annonce la suppression de cette Distance.
	 * Chaque TraitementListener doit se supprimer lui-męme de la liste des listeners de cette Distance.  
	 * 
	 * @param 	e	le TraitementEvent qui est émis vers les TraitementListeners. 
	  */
	private void fireTraitementSuppr(TraitementEvent e)
	{	
	    while (vectorTraitementListener.size()!=0)
			try
			{	((TraitementListener) vectorTraitementListener.get(0)).traitementSuppr(e);
			}
			catch (TraitementIncompletException e1)
			{	e1.printStackTrace();
			}		
	}
	/**
	 * annonce un raffraichissement de cette Distance.
	 * 
	 * @param 	e	le TraitementEvent qui est émis vers les TraitementListeners. 
	  */
	private void fireTraitementRaffr(TraitementEvent e)
	{	
	    for (int i=0;i<vectorTraitementListener.size();i++)
			try
			{	((TraitementListener) vectorTraitementListener.get(i)).traitementRaffr(e);
			}
			catch (TraitementIncompletException e1)
			{	e1.printStackTrace();
			}		
	}
	/**
	 * annonce un raffraichissement de cette Distance.
	 * 
	 * @param 	e	le TraitementEvent qui est émis vers les TraitementListeners. 
	  */
	private void fireTraitementRempl(TraitementEvent e)
	{	for (int i=0;i<vectorTraitementListener.size();i++)
			try
			{	((TraitementListener) vectorTraitementListener.get(i)).traitementRempl(e);
			}
			catch (TraitementIncompletException e1)
			{	e1.printStackTrace();
			}		
	}
	/**
	 * annonce un raffraichissement de cette Distance.
	 * 
	 * @param 	e	le TraitementEvent qui est émis vers les TraitementListeners. 
	  */
	private void fireTraitementCache(TraitementEvent e)
	{	for (int i=0;i<vectorTraitementListener.size();i++)
			try
			{	((TraitementListener) vectorTraitementListener.get(i)).traitementCache(e);
			}
			catch (TraitementIncompletException e1)
			{	e1.printStackTrace();
			}		
	}

//----------------------------------------	
//	This
//----------------------------------------
	/*
	 *
	 */
	public void remove()
	{	fireTraitementSuppr(new TraitementEvent(this));
		//
		parent = null;
		//
		
	}
	
//	----------------------------------------	
//	Simulation
//	----------------------------------------
	/*
	 *
	 */
	public boolean estComplet()
	{	return true;
	}
	/*
	 *
	 */
	public void preparerSimulation()
	{	
//NOTE tester la présence d'au moins un noeudEndogene ? non-observable ?
//NOTE tester la présence de bruit ?
		VectorNoeudFonctionnel vectorNoeudFonctionnel = parent.getVectorNoeudFonctionnel();
		// déterminer la liste de NoeudFonctionnelEndogene et celle de NoeudFonctionnelExogene
		// et la liste des sources de bruit
		vectorNoeudExogene = new VectorNoeudFonctionnel();
		vectorNoeudEndogeneActivation = new VectorNoeudFonctionnel();
		vectorNoeudEndogeneEmission = new VectorNoeudFonctionnel();
		vectorNoeudObserveActivation = new VectorNoeudFonctionnel();
		vectorNoeudObserveEmission = new VectorNoeudFonctionnel();
		vectorBruitTraitementActivation = new VectorFonctionFilsMagnitude();
		vectorBruitTraitementEmission = new VectorFonctionFilsMagnitude();
		vectorBruitMesureActivation = new VectorFonctionFilsMagnitude();
		vectorBruitMesureEmission = new VectorFonctionFilsMagnitude();
		for (int i=0;i<vectorNoeudFonctionnel.size();i++)
		{	InterfaceNoeud n = vectorNoeudFonctionnel.getNoeudFonctionnel(i);
			if (n instanceof NoeudFonctionnelEndogene)
			{	if (((NoeudFonctionnelEndogene) n).estObservableActivation())
				{	vectorNoeudObserveActivation.addNoeudFonctionnel((NoeudFonctionnelEndogene)n); 
					VectorFonctionFilsMagnitude temp = ((NoeudFonctionnelEndogene)n).getSectionActivation().getFonctionMagnitude().getMagnitudeBruit();
					for (int j=0;j<temp.size();j++)
						vectorBruitMesureActivation.addInterfaceFonctionFilsMagnitude(temp.getInterfaceFonctionFilsMagnitude(j));
				}
				else
				{	vectorNoeudEndogeneActivation.addNoeudFonctionnel((NoeudFonctionnelEndogene)n); 
					VectorFonctionFilsMagnitude temp = ((NoeudFonctionnelEndogene)n).getSectionActivation().getFonctionMagnitude().getMagnitudeBruit();
					for (int j=0;j<temp.size();j++)
						vectorBruitTraitementActivation.addInterfaceFonctionFilsMagnitude(temp.getInterfaceFonctionFilsMagnitude(j));
				} 				 
				if (((NoeudFonctionnelEndogene) n).estObservableEmission())
				{	vectorNoeudObserveEmission.addNoeudFonctionnel((NoeudFonctionnelEndogene)n); 
					VectorFonctionFilsMagnitude temp = ((NoeudFonctionnelEndogene)n).getSectionEmission().getFonctionMagnitude().getMagnitudeBruit();
					for (int j=0;j<temp.size();j++)
						vectorBruitMesureEmission.addInterfaceFonctionFilsMagnitude(temp.getInterfaceFonctionFilsMagnitude(j));
				}
				else
				{	vectorNoeudEndogeneEmission.addNoeudFonctionnel((NoeudFonctionnelEndogene)n); 
					VectorFonctionFilsMagnitude temp = ((NoeudFonctionnelEndogene)n).getSectionEmission().getFonctionMagnitude().getMagnitudeBruit();
					for (int j=0;j<temp.size();j++)
						vectorBruitTraitementEmission.addInterfaceFonctionFilsMagnitude(temp.getInterfaceFonctionFilsMagnitude(j));
				}
			} 				 
			else if (n instanceof NoeudFonctionnelExogene)
				vectorNoeudExogene.addNoeudFonctionnel((NoeudFonctionnelExogene)n);
		}
		
		// parametres divers
		h2 = 3;                 			// Squared divided-difference step size
		h = Math.sqrt(h2);          		// Divided difference step-size
		scal1 = 0.5/h;             			// A convenient scaling factor
		scal2 = Math.sqrt((h2-1)/(4*h2*h2));// Another scaling factor		
		
		// variables de dimension
		nu = vectorNoeudExogene.size();
		nxA = vectorNoeudEndogeneActivation.size();
		nxE = vectorNoeudEndogeneEmission.size();
		nyA = vectorNoeudObserveActivation.size();
		nyE = vectorNoeudObserveEmission.size();
		nvA = vectorBruitTraitementActivation.size();
		nvE = vectorBruitTraitementEmission.size();
		nwA = vectorBruitMesureActivation.size();
		nwE = vectorBruitMesureEmission.size();
		
		// matrices initiales
		Matrix p0,x0;
		p0 = new Matrix(nxA+nxE,nxA+nxE,0);
		x0 = new Matrix(nxA+nxE,1,0);
		for (int i=0;i<nxA;i++)
			{	Magnitude mag = ((NoeudFonctionnelEndogene)vectorNoeudEndogeneActivation.getNoeudFonctionnel(i)).getMagnitudeActivation(0);
				p0.set(i,i,mag.getEcartType());
				x0.set(i,0,mag.getMoyenne());
			}
		for (int i=0;i<nxE;i++)
		{	Magnitude mag = ((NoeudFonctionnelEndogene)vectorNoeudEndogeneEmission.getNoeudFonctionnel(i)).getMagnitudeEmission(0);
			p0.set(i+nxA,i+nxA,mag.getEcartType());
			x0.set(i+nxA,0,mag.getMoyenne());
		}
		xbar = x0;		
		Matrix q;
		q = new Matrix(nvA+nvE,nvA+nvE,0);
		for (int i=0;i<nvA;i++)
			q.set(i,i,((MagnitudeBruit)vectorBruitTraitementActivation.getInterfaceFonctionFilsMagnitude(i)).getMagnitude().getEcartType());
		for (int i=0;i<nvE;i++)
			q.set(i+nvA,i+nvA,((MagnitudeBruit)vectorBruitTraitementEmission.getInterfaceFonctionFilsMagnitude(i)).getMagnitude().getEcartType());
		Matrix r;
		r = new Matrix(nwA+nwE,nwA+nwE,0);
		for (int i=0;i<nwA;i++)
			r.set(i,i,((MagnitudeBruit)vectorBruitMesureActivation.getInterfaceFonctionFilsMagnitude(i)).getMagnitude().getEcartType());
		for (int i=0;i<nwE;i++)
			r.set(i+nwA,i+nwA,((MagnitudeBruit)vectorBruitMesureEmission.getInterfaceFonctionFilsMagnitude(i)).getMagnitude().getEcartType());

		// Square root of initial state covariance
		EigenvalueDecomposition eigen;
		Matrix v,d;
		if (p0.isEmpty())
			Sxbar = new Matrix(0,0);
		else
		{	eigen = p0.eig();
			// v et d peuvent ętre complexes
			v = eigen.getV();
			d = diagonaliser(eigen.getD()); 
			// le fait de diagonaliser permet de ne garder que les parties réelles
			// car les éventuelles parties imaginaires ne sont pas sur la diagonale
			// (cf la doc du package des Matrix)
			Sxbar = triag(v.times(sqrtElmt(d)));
		}
		// Square root of process noise covariance
		if (q.isEmpty())
			Sxbar = new Matrix(0,0);
		else
		{	eigen = q.eig();
			v = eigen.getV();
			d = diagonaliser(eigen.getD());
			Matrix sv = v.times(sqrtElmt(d));
			hSv = sv.times(h);
		}
		// Square root of measurement noise cov.
		if (r.isEmpty())
			Sxbar = new Matrix(0,0);
		else
		{	eigen = r.eig();
			v = eigen.getV();
			d = diagonaliser(eigen.getD());
			Matrix sw = v.times(sqrtElmt(d));
			hSw = sw.times(h);
		}

		// Allocate compund matrix consisting of Sxx and Syv
		SxxSxv = new Matrix(nxA+nxE,2*(nxA+nxE+nvA+nvE));   
		// Allocate compund matrix consisting of Syx and Syw
		SyxSyw = new Matrix(nyA+nyE,2*(nxA+nxE+nwA+nwE));  

		// Matrix for storing state estimates
		// Matrix xhat_data = new Matrix(samples+1,nxA+nxE);
		// pas besoin : stocké directement dans les noeuds	
		// Matrix for storing cov. matrices
		Smat = new Matrix(1,(int)(0.5*(nxA+nxE)*(nxA+nxE+1)));
		
		// Index to elem. in Sx
		//[I,J] = find(triu(reshape(1:nx*nx,nx,nx))');
		Matrix temp = new Matrix(nxA+nxE,nxA+nxE);
		for (int i=0;i<(nxA+nxE);i++)
			for (int j=0;j<(nxA+nxE);j++)
				temp.set(i,j,j*(nxA+nxE)+i+1);
		temp = triu(temp);
		temp = temp.transpose();
		Matrix temp2[] = find(temp);
		Matrix mI = temp2[0];
		Matrix mJ = temp2[1];
		// sidx      = sub2ind([nx nx],J,I); 
		sidx = sub2ind(nxA+nxE,nxA+nxE,mJ,mI); 
		// Index into y-vector
		// int yidx  = 1;                  
		// Mean of process noise
		vmean = new Matrix(nvA+nvE,1);       
		// Mean of measurement noise
		wmean = new Matrix(nwA+nwE,1);       
		
		// ----- Initialize state+output equations and linearization -----
		nxnv2 = nxA+nxE+nvA+nvE;
		nxnw2 = nxA+nxE+nwA+nwE;

		// Index to location of Sxv2 and Syw2 in SxxSxv and SyxSyw matrices
		idx_syw2 = 2*(nxA+nxE)+nwA+nwE-1;
		idx_sxv2 = 2*(nxA+nxE)+nvA+nwE-1;		
		
		SxxSxv = acollerDroite(SxxSxv, new Matrix(nxA+nxE,nxnv2));
		SyxSyw = acollerDroite(SyxSyw, new Matrix(nyA+nyE,nxnw2));
	    fireTraitementRaffr(new TraitementEvent(this));		
	}	
	/*
	 *
	 */
	public void evaluer(int instant) throws IOException, ClassNotFoundException, ContenuFichierIncompatibleException
	{	// algo général :
		// 1) mŕj Types et Magnitude d'émission des NoeudExogenes
		// 2) mŕj Types d'activation des NoeudEndogenes
		// 3) mŕj simultanée des Magnitudes d'activation et d'émission des NoeudEndogenes
		// 4) mŕj Types d'émission des NoeudEndogenes
		// 5) mŕj TPT (apprentissage) des NoeudEndogenes


		// ------------------------------------------------------
		// 1) mŕj Types et Magnitude d'émission des NoeudExogenes
		// ------------------------------------------------------
		for (int i=0;i<nu;i++)
		{	NoeudFonctionnelExogene n = (NoeudFonctionnelExogene)vectorNoeudExogene.getNoeudFonctionnel(i);
			if (instant != 0)
			{	Vector resultat;
				resultat = n.lireEntreeFichier();
				// Magnitude
				Magnitude resM = (Magnitude) resultat.get(0);
				n.addMagnitudeEmission(resM);
				// Type
				Type resT = (Type) resultat.get(1);
				try
				{	n.addTypeEmission(resT);
				}
				catch (DomaineIncompatibleException e)
				{	e.printStackTrace();
				}
			}
		}	  


		// ------------------------------------------------------
		// 2) mŕj Types d'activation des NoeudEndogenes
		// ------------------------------------------------------
		// non-observables
		try
		{	for (int i=0;i<nxA;i++)
			{	NoeudFonctionnelEndogene n = (NoeudFonctionnelEndogene)vectorNoeudEndogeneActivation.getNoeudFonctionnel(i);
				// Type d'activation
				Type resT = n.getSectionActivation().getFonctionType().evaluer(instant);
				n.addTypeActivation(resT);
				// mise ŕ jour des variables persistantes du Noeud
				n.getSectionApprentissage().getFonctionTPT().raffraichirNoeud(instant);
			}	  
			// observables (aucune différence en ce qui concerne les Types)
			for (int i=0;i<nyA;i++)
			{	NoeudFonctionnelEndogene n = (NoeudFonctionnelEndogene)vectorNoeudObserveActivation.getNoeudFonctionnel(i);
				// Type d'activation
				Type resT = n.getSectionActivation().getFonctionType().evaluer(instant);
				n.addTypeActivation(resT);
				// mise ŕ jour des variables persistantes du Noeud
				n.getSectionApprentissage().getFonctionTPT().raffraichirNoeud(instant);
			}	  
		}		
		catch (DomaineIncompatibleException e)
		{	e.printStackTrace();
		}

		// ------------------------------------------------------
		// 3) mŕj simultanée des Magnitudes d'activation et d'émission des NoeudEndogenes
		// ------------------------------------------------------
		// d'aprčs le script Matlab de Magnus Norgaard
		Matrix fxbar, sxp, sxm, svp, svm, y0, ybar, syp, sym, swp, swm, Sy, K, y;
		int kx2;
	
		// --- Time update (a'priori update) of state and covariance ---
		if (instant!=1) 
		{	// etat(instant-1) = xhat, bruit = vmean
			fxbar = fevalX(instant);
			xbar = fxbar.times((h2-nxnv2)/h2);
			kx2 = nxA+nxE+nvA+nvE;
			for (int kx=0;kx<nxA+nxE;kx++)
			{	// etat(instant-1) = xhat+h*Sx(:,kx), bruit = vmean
				changerEtat(instant-1,xhat.plus(Sx.getMatrix(0,Sx.getRowDimension()-1,kx,kx).times(h))); 
				sxp = fevalX(instant);
				// etat(instant-1) = xhat-h*Sx(:,kx), bruit = vmean
				changerEtat(instant-1,xhat.minus(Sx.getMatrix(0,Sx.getRowDimension()-1,kx,kx).times(h)));
				sxm = fevalX(instant);
				// SxxSxv(:,kx) = scal1*(sxp-sxm);
				SxxSxv.setMatrix(0,SxxSxv.getRowDimension()-1,kx,kx,sxp.minus(sxm).times(scal1));
				// SxxSxv(:,kx2+kx) = scal2*(sxp+sxm-2*fxbar);
				SxxSxv.setMatrix(0,SxxSxv.getRowDimension()-1,kx2+kx-1,kx2+kx-1,sxp.plus(sxm).minus(fxbar.times(2)).times(scal2));
				// xbar =  xbar + (sxp+sxm)/(2*h2);
				xbar = xbar.plus(sxp.plus(sxm).times(1/(2*h2)));
			}
			changerEtat(instant-1,xhat); 
			for (int kv=0;kv<nvA+nvE;kv++)
			{	// etat(instant-1) = xhat, bruit = hSv(:,kv)
				changerBruitTraitement(hSv.getMatrix(0,hSv.getRowDimension()-1,kv,kv));
				svp = fevalX(instant);
				// etat(instant-1) = xhat, bruit = -hSv(:,kv)
				changerBruitTraitement(hSv.getMatrix(0,hSv.getRowDimension()-1,kv,kv).times(-1));
				svm = fevalX(instant);
				// SxxSxv(:,nx+kv)       = scal1*(svp-svm);
				SxxSxv.setMatrix(0,SxxSxv.getRowDimension()-1,nxA+nxE+kv-1,nxA+nxE+kv-1,svp.minus(svm).times(scal1));
				// SxxSxv(:,idx_sxv2+kv) = scal2*(svp+svm-2*fxbar);
				SxxSxv.setMatrix(0,SxxSxv.getRowDimension()-1,idx_sxv2+kv-1,idx_sxv2+kv-1,(svp.plus(svm).minus(fxbar.times(2))).times(scal2));
				// xbar = xbar + (svp+svm)/(2*h2);
				xbar = xbar.plus(svp.plus(svm).times(1/(2*h2)));
			}

			// Cholesky factor of a'priori estimation error covariance
			Sxbar = triag(SxxSxv);
		}
  
		// --- Measurement update (a posteriori update) ---
		// etat(instant) = xbar, bruit = wmean
		creerEtat(xbar);
		changerBruitMesure(wmean);
		y0 = fevalY(instant);
		if (estObservable(instant))
		{	ybar = y0.times((h2-nxnw2)/h2);
			kx2 = nxA+nxE+nwA+nwE;
			for (int kx=0;kx<nxA+nxE;kx++)
			{	// etat(instant) = xbar+h*Sxbar(:,kx), bruit = wmean
				changerEtat(instant,xbar.plus(Sxbar.getMatrix(0,Sxbar.getRowDimension()-1,kx,kx).times(h)));
				syp = fevalY(instant);
				// etat(instant) = xbar-h*Sxbar(:,kx), bruit = wmean
				changerEtat(instant,xbar.minus(Sxbar.getMatrix(0,Sxbar.getRowDimension()-1,kx,kx).times(h)));
				sym = fevalY(instant);
				// SyxSyw(:,kx) = scal1*(syp-sym);
				SyxSyw.setMatrix(0,SyxSyw.getRowDimension()-1,kx,kx,syp.minus(sym).times(scal1));
				// SyxSyw(:,kx2+kx) = scal2*(syp+sym-2*y0);
				SyxSyw.setMatrix(0,SyxSyw.getRowDimension()-1,kx2+kx-1,kx2+kx-1,syp.plus(sym).minus(y0.times(2)).times(scal2));
				// ybar = ybar + (syp+sym)/(2*h2);
				ybar = ybar.plus(syp.plus(sym).times(1/(2*h2)));    
			}
			changerEtat(instant,xbar); 
			for (int kw=0;kw<nwA+nwE;kw++)
			{	// etat(instant) = xbar, bruit = hSw(:,kw)
				changerBruitMesure(hSw.getMatrix(0,hSw.getRowDimension()-1,kw,kw));
				swp = fevalY(instant);
				// etat(instant) = xbar, bruit = -hSw(:,kw)
				changerBruitMesure(hSw.getMatrix(0,hSw.getRowDimension()-1,kw,kw).times(-1));
				swm = fevalY(instant);
				// SyxSyw(:,nx+kw) = scal1*(swp-swm);
				SyxSyw.setMatrix(0,SyxSyw.getRowDimension()-1,nxA+nxE+kw-1,nxA+nxE+kw-1,swp.minus(swm).times(scal1));
				// SyxSyw(:,idx_syw2+kw) = scal2*(swp+swm-2*y0);
				SyxSyw.setMatrix(0,SyxSyw.getRowDimension()-1,idx_syw2+kw,idx_syw2+kw,swp.plus(swm).minus(y0.times(2)).times(scal2));
				// ybar = ybar + (swp+swm)/(2*h2);          
				ybar = ybar.plus(swp.plus(swm).times(1/(2*h2)));          
			}
    
			// Cholesky factor of a'posteriori output estimation error covariance
			Sy = triag(SyxSyw);
			// K = (Sxbar*SyxSyw(:,1:nx)')/(Sy*Sy');
			K = Sxbar.times(SyxSyw.getMatrix(0,SyxSyw.getRowDimension()-1,0,nxA+nxE-1).transpose()).times(Sy.times(Sy.transpose()).inverse());
			// State estimate
			y = getObservations(instant);
			// xhat = xbar + K*(y(yidx,:)'-ybar);
	  		xhat = xbar.plus(K.times(y.minus(ybar)));
	  	
			// Cholesky factor of a'posteriori estimation error covariance
			// Sx = triag([Sxbar-K*SyxSyw(:,1:nx) K*SyxSyw(:,nx+1:end)]);
			Sx = triag(acollerDroite(Sxbar.minus(K.times(SyxSyw.getMatrix(0,SyxSyw.getRowDimension()-1,0,nxA+nxE-1))),K.times(SyxSyw.getMatrix(0,SyxSyw.getRowDimension()-1,nxA+nxE+1-1,SyxSyw.getColumnDimension()-1))));
			//yidx = yidx + 1; 

		}
		// No observations available at this sampling instant
		else
		{  	// Copy a priori state estimate
			xhat = xbar;                       
			// Copy a priori covariance factor
			Sx   = Sxbar;
			//
			y = y0;         
		}

		// --- Store results ---
		// xhat_data(k,:) = xhat';
		// extraire(vecteur) permet d'extraire un vecteur contenant les valeurs 
		// de la matrice aux coordonnées spécifiées ds le vecteur d'entrée 
		// Smat(k,:) = Sx(sidx)';
		// NOTE ici je considčre que les valeurs de l'instant 0 sont dans la matrice
//		Smat.setMatrix(instant,instant,0,Smat.getColumnDimension()-1,extraire(Sx,sidx.transpose()));
		Smat = acollerBas(Smat,extraire(Sx,sidx).transpose());
		// stocker les moyennes ds les historiques des Noeuds
		changerEtat(instant,xhat);
		changerObservation(instant, y);
			
			
		// ------------------------------------------------------
		// 4) mŕj Types d'émission des NoeudEndogenes
		// ------------------------------------------------------
		// non-observables
		try
		{	for (int i=0;i<nxE;i++)
			{	NoeudFonctionnelEndogene n = (NoeudFonctionnelEndogene)vectorNoeudEndogeneEmission.getNoeudFonctionnel(i);
				Type resT = n.getSectionEmission().getFonctionType().evaluer(instant);
				n.addTypeEmission(resT);
			}	  
			// observables (aucune différence en ce qui concerne les Types)
			for (int i=0;i<nyE;i++)
			{	NoeudFonctionnelEndogene n = (NoeudFonctionnelEndogene)vectorNoeudObserveEmission.getNoeudFonctionnel(i);
				Type resT = n.getSectionEmission().getFonctionType().evaluer(instant);
				n.addTypeEmission(resT);
			}	  
		}		
		catch (DomaineIncompatibleException e)
		{	e.printStackTrace();
		}
		
		
		// ------------------------------------------------------
		// 5) mŕj TPT (apprentissage) des NoeudEndogenes
		// ------------------------------------------------------
		// non-observables
		try
		{	for (int i=0;i<nxE;i++)
			{	NoeudFonctionnelEndogene n = (NoeudFonctionnelEndogene)vectorNoeudEndogeneEmission.getNoeudFonctionnel(i);
				TPT res = n.getSectionApprentissage().getFonctionTPT().evaluer(instant);
				n.addTPT(res);
			}	  
			// observables (aucune différence en ce qui concerne les Types)
			for (int i=0;i<nyE;i++)
			{	NoeudFonctionnelEndogene n = (NoeudFonctionnelEndogene)vectorNoeudObserveEmission.getNoeudFonctionnel(i);
				TPT res = n.getSectionApprentissage().getFonctionTPT().evaluer(instant);
				n.addTPT(res);
			}
		}		
		catch (DomaineIncompatibleException e)
		{	e.printStackTrace();
		}
	}
	/*
	 *
	 */
	public void interrompreSimulation()
	{	// on efface les variables qui ont servi lors de la simulation
		hSv = null;
		hSw = null;
		sidx = null;
		Smat = null;
		Sx = null;
		Sxbar = null;
		SxxSxv = null;
		SyxSyw = null;
		vectorBruitMesureActivation = null;
		vectorBruitMesureEmission = null;
		vectorBruitTraitementActivation = null;
		vectorBruitTraitementEmission = null;
		vectorNoeudEndogeneActivation = null;
		vectorNoeudEndogeneEmission = null;
		vectorNoeudExogene = null;
		vectorNoeudObserveActivation = null;
		vectorNoeudObserveEmission = null;
		vmean = null;
		wmean = null;
		xbar = null;
		xhat = null;
	}	
	/*
	 *
	 */
	public void terminerSimulation()
	{	// SMAT2VAR(Smat) returns a matrix where each column is the variance
		// of a state estimate.
		Matrix var = smat2var(Smat);
		for (int i=0;i<var.getRowDimension();i++)
		{	for (int j=0;j<nxA;j++)
				// NOTE : vérif' le i+1. ça dépend si l'instant 0 est stocké ou pas dans la matrice des instants. suffit de vérifier dans la matrice des états
				((NoeudFonctionnelEndogene) vectorNoeudEndogeneActivation.getNoeudFonctionnel(j)).getMagnitudeActivation(i+1).setEcartType(var.get(i,j));
			for (int j=0;j<nxE;j++)
				// NOTE : vérif' ausi le i+1
				((NoeudFonctionnelEndogene) vectorNoeudEndogeneEmission.getNoeudFonctionnel(j)).getMagnitudeEmission(i+1).setEcartType(var.get(i,j+nxA));
		}
		// on efface les variables qui ont servi lors de la simulation
		hSv = null;
		hSw = null;
		sidx = null;
		Smat = null;
		Sx = null;
		Sxbar = null;
		SxxSxv = null;
		SyxSyw = null;
		vectorBruitMesureActivation = null;
		vectorBruitMesureEmission = null;
		vectorBruitTraitementActivation = null;
		vectorBruitTraitementEmission = null;
		vectorNoeudEndogeneActivation = null;
		vectorNoeudEndogeneEmission = null;
		vectorNoeudExogene = null;
		vectorNoeudObserveActivation = null;
		vectorNoeudObserveEmission = null;
		vmean = null;
		wmean = null;
		xbar = null;
		xhat = null;
	}	
	
	//----------------------------------------	
	//	Fonctions d'évaluation
	//----------------------------------------
	/**
	 * réalise l'évaluation des états ŕ un instant donné.
	 *
	 * @param	instant	instant de l'évaluation ŕ réaliser. 
	 * @return	une matrice (vecteur) contenant le résultat de l'évaluation.
	  */
	private Matrix fevalX(int instant)
	{	Matrix resultat = new Matrix(nxA+nxE,1);
		for (int i=0;i<nxA;i++)
			{	NoeudFonctionnelEndogene n = (NoeudFonctionnelEndogene)vectorNoeudEndogeneActivation.getNoeudFonctionnel(i);
				Magnitude resM = n.getSectionActivation().getFonctionMagnitude().evaluer(instant);
				n.setMagnitudeActivation(instant,resM);
				resultat.set(i,0,resM.getMoyenne());
			}	  
			for (int i=0;i<nxE;i++)
			{	NoeudFonctionnelEndogene n = (NoeudFonctionnelEndogene)vectorNoeudEndogeneEmission.getNoeudFonctionnel(i);
				Magnitude resM = n.getSectionEmission().getFonctionMagnitude().evaluer(instant);
				n.setMagnitudeEmission(instant,resM);
				resultat.set(i+nxA,0,resM.getMoyenne());
			}
		return resultat;
	}
	/**
	 * réalise l'évaluation des observations ŕ un instant donné.
	 *
	 * @param	instant	instant de l'évaluation ŕ réaliser. 
	 * @return	une matrice (vecteur) contenant le résultat de l'évaluation.
	  */
	private Matrix fevalY(int instant)
	{	Matrix resultat = new Matrix(nyA+nyE,1);
		for (int i=0;i<nyA;i++)
			{	NoeudFonctionnelEndogene n = (NoeudFonctionnelEndogene)vectorNoeudObserveActivation.getNoeudFonctionnel(i);
				Magnitude resM = n.getSectionActivation().getFonctionMagnitude().evaluer(instant);
				resultat.set(i,0,resM.getMoyenne());
			}	  
			for (int i=0;i<nyE;i++)
			{	NoeudFonctionnelEndogene n = (NoeudFonctionnelEndogene)vectorNoeudObserveEmission.getNoeudFonctionnel(i);
				Magnitude resM = n.getSectionEmission().getFonctionMagnitude().evaluer(instant);
				resultat.set(i+nyA,0,resM.getMoyenne());
			}

		return resultat;
	}

	//----------------------------------------	
	//	Manipulation des États/Bruits
	//----------------------------------------
	/**
	 * insčre la Magnitude définissant un nouvel instant, pour chaque noeud.
	 *
	 * @param	m		matrice contenant les nouveaux états. 
	  */
	private void creerEtat(Matrix m)
	{	int index = 0;
		for (int i=0;i<nxA;i++)
			{	NoeudFonctionnelEndogene noeud = (NoeudFonctionnelEndogene)vectorNoeudEndogeneActivation.getNoeudFonctionnel(i);
		    	Magnitude magnitude = new Magnitude(noeud.getSectionActivation());
		    	magnitude.setMoyenne(m.get(index,0));
		    	noeud.addMagnitudeActivation(magnitude);
				index ++;	
			}
			for (int i=0;i<nxE;i++)
			{	NoeudFonctionnelEndogene noeud = (NoeudFonctionnelEndogene)vectorNoeudEndogeneEmission.getNoeudFonctionnel(i);
	    		Magnitude magnitude = new Magnitude(noeud.getSectionEmission());
	    		magnitude.setMoyenne(m.get(index,0));
	    		noeud.addMagnitudeEmission(magnitude);
				index ++;	
			}
	}
	/**
	 * échange les états de l'instant passé en paramčtre par
	 * ceux contenus dans la matrice.
	 *
	 * @param	instant	instant des états ŕ échanger. 
	 * @param	m		matrice contenant les nouveaux états. 
	  */
	private void changerEtat(int instant, Matrix m)
	{	int index = 0;
		for (int i=0;i<nxA;i++)
			{	((NoeudFonctionnelEndogene) vectorNoeudEndogeneActivation.getNoeudFonctionnel(i)).getMagnitudeActivation(instant).setMoyenne(m.get(index,0));
				index ++;	
			}
			for (int i=0;i<nxE;i++)
			{	((NoeudFonctionnelEndogene) vectorNoeudEndogeneEmission.getNoeudFonctionnel(i)).getMagnitudeEmission(instant).setMoyenne(m.get(index,0));
				index ++;	
			}
	}
	/**
	 * échange les observations de l'instant passé en paramčtre par
	 * celles contenues dans la matrice.
	 *
	 * @param	instant	instant des observations ŕ échanger. 
	 * @param	m		matrice contenant les nouvelles observations. 
	  */
	private void changerObservation(int instant, Matrix m)
	{	int index = 0;
		for (int i=0;i<nyA;i++)
			{	((NoeudFonctionnelEndogene) vectorNoeudObserveActivation.getNoeudFonctionnel(i)).getMagnitudeActivation(instant).setMoyenne(m.get(index,0));
				index ++;	
			}
			for (int i=0;i<nyE;i++)
			{	((NoeudFonctionnelEndogene) vectorNoeudObserveEmission.getNoeudFonctionnel(i)).getMagnitudeEmission(instant).setMoyenne(m.get(index,0));
				index ++;	
			}
	}
	/**
	 * échange les bruits de traitement par ceux contenus dans la matrice.
	 *
	 * @param	m	matrice contenant les nouveaux bruits de traitement. 
	  */
	private void changerBruitTraitement(Matrix m)
	{	int index = 0;
		for (int i=0;i<nvA;i++)
		{	((MagnitudeBruit) vectorBruitTraitementActivation.getInterfaceFonctionFilsMagnitude(i)).getMagnitude().setMoyenne(m.get(index,0));
			index ++;	
		}
		for (int i=0;i<nvE;i++)
		{	((MagnitudeBruit) vectorBruitTraitementEmission.getInterfaceFonctionFilsMagnitude(i)).getMagnitude().setMoyenne(m.get(index,0));
			index ++;	
		}
	}
	/**
	 * échange les bruits de mesure par ceux contenus dans la matrice.
	 *
	 * @param	m	matrice contenant les nouveaux bruits de mesure. 
	  */
	private void changerBruitMesure(Matrix m)
	{	int index = 0;
		for (int i=0;i<nwA;i++)
		{	((MagnitudeBruit) vectorBruitMesureActivation.getInterfaceFonctionFilsMagnitude(i)).getMagnitude().setMoyenne(m.get(index,0));
			index ++;	
		}
		for (int i=0;i<nwE;i++)
		{	((MagnitudeBruit) vectorBruitMesureEmission.getInterfaceFonctionFilsMagnitude(i)).getMagnitude().setMoyenne(m.get(index,0));
			index ++;	
		}
	}
	
	//----------------------------------------	
	//	Observations
	//----------------------------------------
	/**
	* détermine si tous les Noeuds observables le sont bien ŕ 
	* l'instant passé en paramčtre.
	*
	* @param	instant					instant ŕ tester. 
	* @return	vrai si tous les Noeuds observables le sont ŕ cet instant, faux sinon.
	* @throws	IOException				si une erreur se produit lors de l'accčs au fichier.
	* @throws	ClassNotFoundException	si une erreur se produit lors de l'accčs au fichier.	
	 */
	private boolean estObservable(int instant) throws IOException, ClassNotFoundException
	{	boolean resultat = true;
		for (int i=0;i<nyA;i++)
				resultat = resultat && ((NoeudFonctionnelEndogene) vectorNoeudObserveActivation.getNoeudFonctionnel(i)).estObservableActivation(instant);
			for (int i=0;i<nyE;i++)
				resultat = resultat && ((NoeudFonctionnelEndogene) vectorNoeudObserveEmission.getNoeudFonctionnel(i)).estObservableEmission(instant);
		return resultat;
	} 
	/**
	* renvoie la matrice (vecteur) des valeurs observées
	* ŕ l'instant passé en paramčtre.
	*
	* @param	instant					pour lequel on veut les observations. 
	* @return	observations ŕ l'instant donné.
	* @throws	IOException				si une erreur se produit lors de l'accčs au fichier.
	* @throws	ClassNotFoundException	si une erreur se produit lors de l'accčs au fichier.	
	 */
	private Matrix getObservations(int instant) throws IOException, ClassNotFoundException
	{	Matrix resultat = new Matrix(nyA+nyE,1);
		for (int i=0;i<nyA;i++)
		{	NoeudFonctionnelEndogene noeud = ((NoeudFonctionnelEndogene)vectorNoeudObserveActivation.getNoeudFonctionnel(i));
		    Magnitude magnitude = noeud.getObservationActivationValeur(instant);
		    noeud.setMagnitudeActivation(instant, magnitude);
			resultat.set(i,0,magnitude.getMoyenne());
		}
		for (int i=0;i<nyE;i++)
		{	NoeudFonctionnelEndogene noeud = ((NoeudFonctionnelEndogene)vectorNoeudObserveEmission.getNoeudFonctionnel(i));
		    Magnitude magnitude = noeud.getObservationEmissionValeur(instant);
		    noeud.setMagnitudeEmission(instant, magnitude);
			resultat.set(i+nyA,0,magnitude.getMoyenne());
		}
		return resultat;	
	}

	//----------------------------------------	
	//	Manipulation des Matrix
	//	Méthodes adaptées de Norgaard
	//----------------------------------------
//  NOTE possibilité de bug du ŕ l'apparition de valeurs complexes, cf appel de cette méthode	
	/**
	* méthode adaptée de celle de Magnus Norgaard.
	* <p>
	*     S=triag(A) uses a Householder transformation of the rectangular
	*     matrix A to produce a square and upper triangular matrix S with
	*     the property S*S' = A*A'.
	*<p>
	*     Literature:
	*     The algorithm is described in Grewal & Andrews: "Kalman Filtering - 
	*     Theory and Practice" Prentice-Hall, 1993, pp. 234. However, there 
	*     are certain cases in which their implementation fails. These have been
	*     handled as described in G.W. Stewart: "Matrix Algorithms, 
	*     Vol. 1: Basic Decompositions", SIAM 1988.
	*<p>
	* Written by: Magnus Norgaard, IMM/IAU, Technical University of Denmark
	* LastEditDate: Apr. 15, 2000
	* 
	* 
	*
	* @param	m	la matrice rectangulaire que l'on veut triangulariser. 
	* @return	la matrice triangulaire. 
	 */
	private Matrix triag(Matrix m)
	{	Matrix a = m.copy();
	    // Rows and columns of A
	    // [n,rn]=size(A);
		int n = a.getRowDimension();
		int rn = a.getColumnDimension(); 
		// r = rn-n;
		int r = rn-n;
		double ny;
		// v = zeros(n,1);
		Matrix v = new Matrix(n,1);
		// u = zeros(1,rn);
		Matrix u = new Matrix(1,rn);
		// S = zeros(n,n);
		Matrix S = new Matrix(n,n);

		// for k=n:-1:1,
		for (int k=(n-1);k>=0;k--)
		{	// u(1:r+k) = A(k,1:r+k);
			u.setMatrix(0,0,0,r+k,a.getMatrix(k,k,0,r+k));
			// ny = norm(u(1:r+k));
		   	ny = u.getMatrix(0,0,0,r+k).normF();
		    // if ny==0,
		   	if (ny==0)
		        // u(r+k) = sqrt(2);
			 	u.set(0,r+k,Math.sqrt(2));
		    // else
		   	else
			{	// u(1:r+k) = u(1:r+k)/ny;
				u.setMatrix(0,0,0,r+k,u.getMatrix(0,0,0,r+k).times(1/ny));
			    // if u(r+k)>=0,  
				if (u.get(0,r+k)>=0)  
				{	// u(r+k) = u(r+k) + 1;
				    u.set(0,r+k,u.get(0,r+k)+1);
			        // ny = -ny;
					ny = -ny;
				}
			    // else
			  	else
			        // u(r+k) = u(r+k) - 1;
					u.set(0,r+k,u.get(0,r+k)-1);
				//u(1:r+k) = u(1:r+k)/sqrt(abs(u(r+k)));
				u.setMatrix(0,0,0,r+k,u.getMatrix(0,0,0,r+k).times(1/Math.sqrt(Math.abs(u.get(0,r+k)))));
			}
			//v(1:k-1,1) = A(1:k-1,1:r+k)*u(1,1:r+k)';
			v.setMatrix(0,k-1,0,0,a.getMatrix(0,k-1,0,r+k).times(u.getMatrix(0,0,0,r+k).transpose()));
			//A(1:k-1,1:r+k) = A(1:k-1,1:r+k) - v(1:k-1,1)*u(1,1:r+k);
			a.setMatrix(0,k-1,0,r+k,a.getMatrix(0,k-1,0,r+k).minus(v.getMatrix(0,k-1,0,0).times(u.getMatrix(0,0,0,r+k))));
		    //S(1:k-1,k) = A(1:k-1,r+k);
		    S.setMatrix(0,k-1,k,k,a.getMatrix(0,k-1,r+k,r+k));
		   	//S(k,k)=ny;
		   	S.set(k,k,ny);
		}
		return S;
	}
	/**
	* équivalent ŕ la fonction smat2cov de Norgaard qui permet de passer 
	* d'une matrice triangularisée aux variances.
	* <p>
	* SMAT2VAR   Calculate variance estimate for each state.
	*    SMAT2VAR(Smat) returns a matrix where each column is the variance
	*    of a state estimate. 'Smat' is a matrix where each row contains elements 
	*    of (the upper triangular part of) the Cholesky factor of a covariance 
	*    matrix.
	*
	* @param	sMat	matrice triangulaire. 
	* @return	matrice de covariance. 
	 */
	private Matrix smat2var(Matrix sMat)
	{	// Number of states
		// states = round(-0.5 + sqrt(0.25 + 2*size(sMat,2)));
		int states = (int) Math.round(-0.5 + Math.sqrt(0.25 + 2*sMat.getColumnDimension())); 
		// Index to first element in each row
		// idx = cumsum ([1 states:-1:2]);
		Matrix idx = new Matrix(1,states);
		idx.set(0,0,1);
		for (int i=1;i<states;i++)
		    idx.set(0,i,states+1-i);
		for (int i=1;i<states;i++)
			idx.set(0,i,idx.get(0,i-1)+idx.get(0,i));
		// Index to last element in each row
		//idx2 = idx+[states-1:-1:0];
		Matrix temp = new Matrix(1,states);
		for (int i=0;i<states;i++)
			temp.set(0,i,states-1-i);
		Matrix idx2 = idx.plus(temp);
		// Allocate space for variance matrix
		// vars = zeros(size(sMat,1),states);
		Matrix vars = new Matrix(sMat.getRowDimension(),states);   
		// for k=1:states,
		//   vars(:,k) = sum(Smat(:,idx(k):idx2(k)).*Smat(:,idx(k):idx2(k)),2);
		// end
		for (int k=0;k<states;k++)
		{	//   Smat(:,idx(k):idx2(k))
		    Matrix temp2 = sMat.getMatrix(0,sMat.getRowDimension()-1,(int)idx.get(0,k)-1,(int)idx2.get(0,k)-1);
			//   Smat(:,idx(k):idx2(k)).*Smat(:,idx(k):idx2(k))
			temp2 = temp2.arrayTimes(sMat.getMatrix(0,sMat.getRowDimension()-1,(int)idx.get(0,k)-1,(int)idx2.get(0,k)-1));
			//   sum(Smat(:,idx(k):idx2(k)).*Smat(:,idx(k):idx2(k)),2);
			Matrix temp3 = new Matrix(sMat.getRowDimension(),1);
			for (int i=0;i<temp2.getRowDimension();i++)
				for (int j=0;j<temp2.getColumnDimension();j++)
					temp3.set(i,0,temp3.get(i,0)+temp2.get(i,j));
			//   vars(:,k) = sum(Smat(:,idx(k):idx2(k)).*Smat(:,idx(k):idx2(k)),2);
	   		vars.setMatrix(0,sMat.getRowDimension()-1,k,k,temp3);
		}
				
		return vars;
	}
	//----------------------------------------	
	//	Manipulation des Matrix
	//	Méthodes adaptées de Matlab
	//----------------------------------------
	/**
	* équivalent de la fonction 'find' de Matlab.
	* <p>
	* FIND   Find indices of nonzero elements.
	* [I,J] = FIND(X) returns the row and column indices of
	* the nonzero entries in the matrix X.  This is often used
	* with sparse matrices.
	* <p>
	* rque : version simplifiée par rapport ŕ celle de matlab, puisque
	* je ne prends qu'un seul paramčtre : la matrice ŕ traiter. De plus,
	* il y a aussi une restriction au niveau de ce qui est renvoyé, puisque
	* il s'agit seulement de coordonnées classique (pas les coordonnées linéaires
	* ou les valeurs qui sont non-nulles).
	*
	* @param	a	la matrice ŕ tester. 
	* @return	renvoie un tableau de 2 matrices (vecteurs) contenant les
	* 			coordonnées des éléments non nuls de la matrice passée en paramčtre. 
	 */
	private Matrix[] find(Matrix a)
	{	Matrix resultat[] = new Matrix[2];
		Vector tempI = new Vector();
		Vector tempJ = new Vector();
		for (int i=0;i<a.getRowDimension();i++)
			for (int j=0;j<a.getColumnDimension();j++)
				if (a.get(i,j) != 0)
				{	tempI.add(new Integer(i)); 
					tempJ.add(new Integer(j));
				} 
		resultat[0] = new Matrix(tempI.size(),1);
		resultat[1] = new Matrix(tempI.size(),1);
		for (int i=0;i<tempI.size();i++)
		{	resultat[0].set(i,0,((Integer)tempI.get(i)).intValue());
			resultat[1].set(i,0,((Integer)tempJ.get(i)).intValue());
		}
		return resultat;
	}
	/**
	* équivalent de la fonction 'sub2ind' de Matlab.
	* <p>
	* SUB2IND Linear index from multiple subscripts.
	* SUB2IND is used to determine the equivalent single index
	* corresponding to a given set of subscript values.
	* IND = SUB2IND(SIZ,I,J) returns the linear index equivalent to the
	* row and column subscripts in the arrays I and J for a matrix of
	* size SIZ.
	* <p>
	* rque : j'ai implémenté le cas particulier utilisé ici (matrice ŕ 2 dimensions), pas le cas général
	* ŕ la différence de la version matlab, ici on commence la numérotation ŕ 0. Par exemple, pour une matrice
	* 2x4, ça donne :
	* [ 0 2 4 6
	*   1 3 5 7 ]
	* et non pas (en matlab) :
	* [ 1 3 5 7
	*   2 4 6 8 ]
	*<p>
	* exemple :
	* [1],[1] -> [3]
	* [0,1],[1,1] -> [2,3]
	* [0,1],[1,1] -> [2,3]
	* [0,1 ,[0,0     [0,1
	*  0,1]  1,1] ->  2,3]
	* etc.
	* 
	* @param	siz1	nombre de lignes de la matrice ŕ numéroter. 
	* @param	siz2	nombre de colonnes de la matrice ŕ numéroter.
	* @param	mI		matrice des lignes des éléments dont on veut les coordonnées. 
	* @param	mJ		matrice des colonnes des éléments dont on veut les coordonnées. 
	* @return	la matrice de coordonnées linéaires des éléments dans la matrice de dimensions spécifiées.
	* 
	 */
	private Matrix sub2ind(int siz1, int siz2, Matrix mI, Matrix mJ)
	{	Matrix resultat = new Matrix(mI.getRowDimension(),mI.getColumnDimension());
			
		for (int i=0;i<mI.getRowDimension();i++)
		    for (int j=0;j<mI.getColumnDimension();j++)
		    resultat.set(i,j,mJ.get(i,j)*siz1+mI.get(i,j));
		return resultat;
	}

	//----------------------------------------	
	//	Manipulation des Matrix
	//	Méthodes basiques
	//----------------------------------------
	/**
	* renvoie la partie réelle de la racine carrée de chaque élément
	* de la matrice passée en paramčtre.
	*
	* @param	a	matrice dont on veut les racines carrées. 
	* @return	la matrice contenant les racines carrées. 
	 */
	private Matrix sqrtElmt(Matrix a)
	{	Matrix resultat = new Matrix(a.getRowDimension(), a.getColumnDimension());
			
		for(int i=0;i<a.getRowDimension();i++)
			for(int j=0;j<a.getColumnDimension();j++)
			{	double valeur = a.get(i,j); 
				if (valeur>=0)
					resultat.set(i,j,Math.sqrt(valeur));
				else
					// si le nbre est négatif, la partie réelle
					// de la racine carrée (complexe) est toujours 0
					resultat.set(i,j,0);
			}
		return resultat;
	}	
	// 
	/**
	* renvoie la partie triangulaire supérieure d'une matrice
	* sous la forme d'une nouvelle matrice
	*
	* @param	a	la matrice originale. 
	* @return	la partie triangulaire supérieure de la matrice 'a'. 
	 */
	private Matrix triu(Matrix a)
	{	Matrix resultat = a.copy();
		for (int i=0;i<a.getRowDimension();i++)
			for (int j=0;j<a.getColumnDimension();j++)
				if (i>j)
					resultat.set(i,j,0);
		return resultat;
	}
	/**
	* accolle la seconde matrice ŕ gauche de la premiere.
	* Leurs dimensions doivent ętre compatibles
	*
	* @param	a	permičre matrice ŕ accoller. 
	* @param	b	seconde matrice ŕ accoller. 
	* @return	matrice résultant de l'accollage. 
	 */
	private Matrix acollerDroite(Matrix a, Matrix b)
	{	Matrix resultat = new Matrix(a.getRowDimension(), a.getColumnDimension()+b.getColumnDimension());
		for (int i=0;i<a.getRowDimension();i++)
			for (int j=0;j<a.getColumnDimension();j++)
				resultat.set(i,j,a.get(i,j));
		for (int i=0;i<b.getRowDimension();i++)
			for (int j=0;j<b.getColumnDimension();j++)
				resultat.set(i,a.getColumnDimension()+j,b.get(i,j));
		return resultat;
	}
	/**
	* accolle la seconde matrice en bas de la premiere.
	* Leurs dimensions doivent ętre compatibles
	*
	* @param	a	permičre matrice ŕ accoller. 
	* @param	b	seconde matrice ŕ accoller. 
	* @return	matrice résultant de l'accollage. 
	 */
	private Matrix acollerBas(Matrix a, Matrix b)
	{	Matrix resultat = new Matrix(a.getRowDimension()+b.getRowDimension(), a.getColumnDimension());
		for (int i=0;i<a.getRowDimension();i++)
			for (int j=0;j<a.getColumnDimension();j++)
				resultat.set(i,j,a.get(i,j));
		for (int i=0;i<b.getRowDimension();i++)
			for (int j=0;j<b.getColumnDimension();j++)
				resultat.set(a.getRowDimension()+i,j,b.get(i,j));
		return resultat;
	}
	/**
	* extrait un vecteur colonne contenant les valeurs de la matrice
	* 'm' aux coordonnées linéaires spécifiées dans le vecteur colonne 'coord'.
	* exemple : 
	* 	m = [ 1 2 3 4
	*         5 6 7 8 ]
	*   coord = [1
	*            2
	*            6 ]
	* 	resultat = [5
	*               2
	*               4 ]
	* La matrice m est numérotée comme dans la fonction sub2ind, c'est ŕ dire :
	* 	m = [ 0 2 4 6
	*         1 3 5 7 ]
	*
	* @param	m		matrice dont on veut extraire des valeurs. 
	* @param	coord	matrice (vecteur) contenant les coordonnées linéaires dans 'm' des valeurs ŕ extraire. 
	* @return	matrice (vecteur) contenant les valeurs extraites. 
	 */
	private Matrix extraire(Matrix m, Matrix coord)
	{	Matrix resultat = new Matrix(coord.getRowDimension(),1);
		
		for (int i=0;i<coord.getRowDimension();i++)
		{	int index = (int) coord.get(i,0);
			resultat.set(i,0,m.get(index%m.getRowDimension(),index/m.getRowDimension()));
		}
		
		return resultat;
	}	
	/**
	*	renvoie une matrice ne contenant que la diagonale de la matrice passée en paramčtre.
	* @param	m		matrice dont qu'on veut diagonaliser. 
	* @return	matrice contenant uniquement la diagonale de m, et des zéros ailleurs. 
	 */
	public Matrix diagonaliser(Matrix m)
	{	Matrix resultat = new Matrix(m.getRowDimension(),m.getColumnDimension(),0);
	    for (int i=0; i<m.getRowDimension(); i++)
	    	resultat.set(i,i,m.get(i,i));
	    return resultat;
	}

	
//	----------------------------------------	
//	  Copier/Coller
//	----------------------------------------
	/*
	 *
	 */
	public InterfaceMoteur cloner()
	{	MoteurDD2 resultat = new MoteurDD2(parent);
		resultat.setNonModifie();
	    return resultat;
	}
	/*
	 *
	 */
	public InterfaceHomonyme getHomonyme(Modele modele) throws HomonymeIntrouvableException
	{	InterfaceHomonyme resultat;
		if (modele == getModele())
		    resultat = this;
		else
		{	resultat = modele.getMoteur();
		}
		return resultat;
	}

//	----------------------------------------	
//	InterfaceAffichable
//	----------------------------------------
	/*
	 *
	 */
	public void remplacer(InterfaceAffichable ia)
	{	fireTraitementRempl(new TraitementEvent(this,(MoteurDD2)ia));
	}
	/*
	 *
	 */
	public void cacher()
	{	fireTraitementCache(new TraitementEvent(this));
	}
	
//	----------------------------------------	
//	Modification
//	----------------------------------------
  /*
   *
   */
  public void setModifie()
  {	modifie = true;
  }
  /*
   *
   */
  public boolean estModifie()
  {	return modifie;  
  }
  /*
   *
   */
  public void setNonModifie()
  {	modifie = false;
    fireTraitementRaffr(new TraitementEvent(this));
  }
}