// Copyright (C) 2009 - CEA - Jean-Marc Martinez
//
// This file must be used under the terms of the GNU Lesser General Public License license :
// http://www.gnu.org/copyleft/lesser.html

// Modèle mathématique à analyser
// fonctions noyaux de Dirichlet (voir bench JMM Mascot Num)
function y = MyFunction (x , alpha)
  y=1;
  for i=1:size(x,2)
     if (x(i) == 0) then
       a = 2 * i + 1;
     else
       a = sin((2*i+1) * %pi * x(i))/sin(%pi * x(i));
     end
     b = (a - 1) / sqrt(2 * i);
     y = y * (1 + alpha(i) * b);
  end
endfunction
// Dimension
nx = 3;
// Coefficients du modele (utilises par MyFunction)
for i=1:nx
  alpha(i) = 1. / i;
end
// Calcul de la moyenne de Y
muy = TODO
// Calcul de la variance de Y
vay = TODO
// Calcul des indices du premier ordre et des indices de sensibilité totale
sx = ones(nx,1);
TODO
sg = ones(nx,1);
TODO
// On regroupe les variables
srvu = setrandvar_new();
for i=1:nx
  rvu(i) = randvar_new("Uniforme",0,1);
  setrandvar_addrandvar(srvu, rvu(i));
end
// Réalisation d'un échantillon via la méthode SRS
np = 5000;
setrandvar_buildsample(srvu,"MonteCarlo",np);

// Plan d'expériences : matrice A
A = setrandvar_getsample(srvu);
// Plan d'expériences : matrice B
setrandvar_buildsample(srvu,"MonteCarlo",np);
B = setrandvar_getsample(srvu);
// Réalisation des plans A et B
for k=1:np
  ya(k) = MyFunction(A(k,:), alpha);
  yb(k) = MyFunction(B(k,:), alpha);
end
// Estimation des indices par la méthode de Sobol
mprintf("\nSensitivity analysis\n");
for i=1:nx
  C=B;
  C(:,i) = A(:,i);
  TODO
  mprintf("\nSensitivity index of variable X[%d]\n",i);
  mprintf("First index is : %12.4e - (expected %12.4e)\n", s1, sx(i));
  mprintf("Total index is : %12.4e - (expected %12.4e)\n", st, sg(i));
end