// 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 : fonction produit
function y = MyFunction (x , alpha)
  y=1;
  for i=1:size(x,2)
    // %eps est la precison machine
    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
d = 3;

// Coefficients du modele (utilises par MyFunction)
for i=1:d
  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 totaux
sx = ones(d,1);
TODO
sg = ones(d,1);
TODO
// Par defaut, Nisp cree un groupe de d variables aleatoires uniforme [0,1]
srvx = setrandvar_new(d);
// Specification d'un plan : quadrature tensorisée
// formule exacte pour un polynôme de degré <= niveau
niveau = 10;
setrandvar_buildsample( srvx, "Quadrature", niveau);
// Polynome de chaos
noutput = 1;
pc = polychaos_new ( srvx , noutput );
// Réalisation du plan d'expériences numériques
np = setrandvar_getsize(srvx);
polychaos_setsizetarget(pc,np);
nx = polychaos_getdiminput(pc);
ny = polychaos_getdimoutput(pc);
inputdata = zeros(nx);
outputdata = zeros(ny);
for k=1:np
  inputdata = setrandvar_getsample(srvx,k);
  outputdata = MyFunction(inputdata, alpha);
  polychaos_settarget(pc,k,outputdata);
end
// Calcul des coefficients par intégration numérique
// rappel degre <= niveau de la formule d'integration
polychaos_setdegree(pc, degre);
polychaos_computeexp(pc,srvx,"Integration");
// Edition des indices de sensibilité du premier ordre
// et des indices de sensibilite totale et rappel des valeurs exactes
TODO