The NISP Module
In this page, we present the Non Intrusive Spectral Projection (NISP) module available in Scilab since the 29th of September 2009.
Contents
Introduction
This toolbox allows to approximate a given model, which is associated with input random variables.
The module provides the following components :
- “nisp” provides function to configure the global behaviour of the toolbox. This allows to startup and shutdown the library, configure and quiery the verbose level or initialize the seed of the random number generator.
- “randvar” is the class which allows to manage a random variable. Various types of random variables are available, including uniform, normal, exponential, etc...
- “setrandvar”, is the class which allows to manage a set of random variables. Several methods are available to build a sampling from a set of random variables. We can use, for example, a Monte-Carlo sampling or a Sobol low discrepancy sequence. This feature allows to use the class as Design of Experiment tool (DOE).
- “polychaos” is the class which allows to manage a polynomial chaos expansion. The coefficients of the expansion are computed based on given numerical experiments which creates the association between the inputs and the outputs. Once computed, the expansion can be used as a regular function. The mean, standard deviation or quantile can also be directly retrieved.
The current toolbox provides an object-oriented approach of the C++ NISP library.
This toolbox has been created in the context of the OPUS project :
within the workpackage 2.1.1 "Construction de méta-modèles".
This project has received funding by Agence Nationale de la recherche :
http://www.agence-nationale-recherche.fr/
Features
Main Features:
- randvar:
- Manage various types of random variables
- uniform, normal, exponential, log-normal
- setrandvar:
- Manage various sampling methods for sets of random variables
- Monte-Carlo, Sobol Quasi-Random, Latin Hypercube Sampling, LHS Max Min sampling, and various samplings based on Smolyak Cubature points.
- polychaos:
- Manage polynomial chaos expansion and get specific outputs
- mean, variance, sensitivity indices, quantiles, Wilks quantiles, correlation, etc...
- Generate a stand-alone C source code which computes the output of the polynomial chaos expansion.
Configuration Functions:
- nisp_destroyall : Destroy all current objects.
- nisp_getpath : Returns the path to the current module.
- nisp_initseed : Sets the seed of the uniform random number generator.
- nisp_printall : Prints all current objects.
- nisp_shutdown : Shuts down the NISP toolbox.
- nisp_startup : Starts up the NISP toolbox.
- nisp_verboselevelget : Returns the current verbose level.
- nisp_verboselevelset : Sets the current verbose level.
Sensitivity Analysis Functions:
- nisp_bruteforcesa : Compute sensitivity indices by brute force.
- nisp_sobolsa : Compute sensitivity indices by Sobol, Ishigami, Homma.
Support functions:
- nisp_buildlhs : Creates a LHS design
- nisp_corrcoef : Returns the linear correlation coefficient of x and y.
- nisp_cov : Returns the empirical covariance matrix of x and y.
- nisp_erfcinv : Computes the inverse erfc function.
- nisp_expcdf : Computes the Exponential CDF.
- nisp_expinv : Computes the Exponential quantile.
- nisp_exppdf : Computes the Exponential PDF.
- nisp_lognormalcdf : Computes the Lognormal CDF.
- nisp_lognormalinv : Computes the Lognormal quantile.
- nisp_lognormalpdf : Computes the Lognormal PDF.
Test functions:
- nisp_ishigami : Returns the Ishigami function.
- nisp_ishigamisa : Exact sensitivity analysis for the Ishigami function
- nisp_product : Returns the value of the Product function
- nisp_productsa : Exact sensitivity analysis for the Product function
- nisp_sum : Returns the value of the Product function
- nisp_sumsa : Returns the sensitivity indices of the Sum function
How to get the module
Forge
The Scilab Forge hosts the sources of the module:
http://forge.scilab.org/index.php/p/nisp/
How to install the module
The module is available under ATOMS :
http://atoms.scilab.org/toolboxes/NISP
To install it, type the following statement in the Scilab console (Scilab version >= 5.1):
atomsInstall("NISP")
Documentation
See in the help provided in the help/en_US directory of the toolbox for more information about its use. Use cases are presented in the demos directory.
This module has been presented at the "42èmes Journées de Statistique", du 24 au 28 mai 2010 - Marseille, France. The title was "Polynômes de chaos sous Scilab via la librairie NISP", Michaël Baudin, Jean-Marc Martinez, 2010 and the talk has been given the 28th of May 2010.
Slides (25 pages) : NispScilab-slides-Baudin-Martinez-2010.pdf
Paper (6 pages) : NispScilab-paper-Baudin-Martinez-2010.pdf
The following User's Manual presents the toolbox :
Manual, v0.2, 2011 (101 pages) : nisp_toolbox_manual_v0.2.pdf
Manual, v0.1, 2009 (51 pages) : nisp_toolbox_manual.pdf
The NISP module for Scilab has been presented at the OPUS Final workshop, "Computer experiments and uncertainty analysis" Friday, October 21, 2011 Henri Poincaré Institute, France. The title was "Polynômes du chaos et analyse de sensibilité : zoom sur les outils disponibles dans Scilab", Michaël Baudin, Jean-Marc Martinez, 2011
Slides (26 pages) : slides-workshopOPUS-October2011.pdf
Teaching
The NISP module is used to teach sensitivity analysis to Master's students at INSTN. The master is "Master de Modélisation et Simulation" (M2S) :
http://www.maisondelasimulation.fr/m2s/
This session is within the course "Calcul hautes performances", Informatique scientifique approfondie (I1), Traitement des incertitudes (I1c).
Session 2009-2010
TP du 30 Novembre 2009 (PDF)
Examen du 17/12/2009 - Sujet : (PDF)
Examen du 17/12/2009 - Script Exercice 1 : examen_2009_sujet_ex1.sce
Examen du 17/12/2009 - Script Exercice 2 : examen_2009_sujet_ex2.sce
Examen du 17/12/2009 - Correction Exercice 1 : examen_2009_correction_ex1.sce
Examen du 17/12/2009 - Correction Exercice 2 : examen_2009_correction_ex2.sce
Session 2010-2011
Examen du 14/02/2011 - Sujet et Correction : (PDF)
Examen du 14/02/2011 - Script Exercice 1 : examen_2011_sujet_ex1.sce
Examen du 14/02/2011 - Script Exercice 2 : examen_2011_sujet_ex2.sce
Examen du 14/02/2011 - Correction Exercice 1 : examen_2011_correction_ex1.sce
Examen du 14/02/2011 - Correction Exercice 2 : examen_2011_correction_ex2.sce
Session 2011-2012
Examples
In this section, we present two simple examples of use of this module.
The product function
In the following example, we compute the polynomial chaos expansion of a very simple function, taking two input variables and returning one output variable. The output result is computed by the Exemple function as the product of the two inputs. The first variable is associated with a Normal law, while the second is associated with a Uniform law.
// Copyright (C) 2009 - CEA - Jean-Marc Martinez
// Copyright (C) 2008-2009 - INRIA - Michael Baudin
//
// This file must be used under the terms of the GNU Lesser General Public License license :
// http://www.gnu.org/copyleft/lesser.html
//
// Test the chaos polynomial decomposition on the product x1*x2 function.
// Two random variables:
// * a normal random variable with mu=1.0 and sigma=0.5,
// * a uniform random variable in [1,2.5].
// Create a Quadrature sampling and compute the coefficients of the polynomial
// by integration.
// Use a degree 2 polynomial.
function y = Exemple (x)
// Returns the output y of the product x1 * x2.
// Parameters
// x: a np-by-nx matrix of doubles, where np is the number of experiments, and nx=2.
// y: a np-by-1 matrix of doubles
y(:,1) = x(:,1).* x(:,2)
endfunction
// 1. Two stochastic (normalized) variables
vx1 = randvar_new("Normale");
vx2 = randvar_new("Uniforme");
// 2. A collection of stoch. variables.
srvx = setrandvar_new();
setrandvar_addrandvar ( srvx,vx1);
setrandvar_addrandvar ( srvx,vx2);
// 3. Two uncertain parameters
vu1 = randvar_new("Normale",1.0,0.5);
vu2 = randvar_new("Uniforme",1.0,2.5);
// 4. A collection of uncertain parameters
srvu = setrandvar_new();
setrandvar_addrandvar ( srvu,vu1);
setrandvar_addrandvar ( srvu,vu2);
// 5. Create the Design Of Experiments
degre = 2;
setrandvar_buildsample(srvx,"Quadrature",degre);
setrandvar_buildsample( srvu , srvx );
// 6. Create the polynomial
ny = 1;
pc = polychaos_new ( srvx , ny );
np = setrandvar_getsize(srvx);
polychaos_setsizetarget(pc,np);
// 7. Perform the DOE
inputdata = setrandvar_getsample(srvu);
outputdata = Exemple(inputdata);
polychaos_settarget(pc,outputdata);
// 8. Compute the coefficients of the polynomial expansion
polychaos_setdegree(pc,degre);
polychaos_computeexp(pc,srvx,"Integration");
// 9. Get the sensitivity indices
average = polychaos_getmean(pc);
var = polychaos_getvariance(pc);
mprintf("Mean = %f\n",average);
mprintf("Variance = %f\n",var);
S = polychaos_getindexfirst(pc);
mprintf("S1 = %f\n",S(1));
mprintf("S2 = %f\n",S(2));
ST = polychaos_getindextotal(pc);
mprintf("ST1 = %f\n",ST(1));
mprintf("ST2 = %f\n",ST(2));
// Clean-up
polychaos_destroy(pc);
randvar_destroy(vu1);
randvar_destroy(vu2);
randvar_destroy(vx1);
randvar_destroy(vx2);
setrandvar_destroy(srvu);
setrandvar_destroy(srvx);The previous script produces the following output.
Mean = 1.750000 Variance = 1.000000 S1 = 0.765625 S2 = 0.187500 ST1 = 0.812500 ST2 = 0.234375
The Ishigami function
In the following example, we compute the polynomial chaos expansion of the Ishigami function, which has 3 inputs and 1 output. The three input variables are uniform random variables in the interval [-Pi,Pi].
// Copyright (C) 2009 - CEA - Jean-Marc Martinez
// Copyright (C) 2009-2011 - INRIA - Michael Baudin
//
// This file must be used under the terms of the GNU Lesser General Public License license :
// http://www.gnu.org/copyleft/lesser.html
//
// Test the chaos polynomial decomposition on the Ishigami function.
// Three uniform variables in [-pi,pi].
// Create a Petras sampling and compute the coefficients of the polynomial
// by integration.
//
// 3 random variable uniform in [-pi,pi]
a=7.;
b=0.1;
exact = nisp_ishigamisa ( a , b );
// 1. Create a group of stochastic variables
nx = 3;
srvx = setrandvar_new( nx );
// 2. Create a group of uncertain variables
rvu1 = randvar_new("Uniforme",-%pi,%pi);
rvu2 = randvar_new("Uniforme",-%pi,%pi);
rvu3 = randvar_new("Uniforme",-%pi,%pi);
srvu = setrandvar_new();
setrandvar_addrandvar ( srvu, rvu1);
setrandvar_addrandvar ( srvu, rvu2);
setrandvar_addrandvar ( srvu, rvu3);
// 3. Create a Petras sampling
degre = 9;
setrandvar_buildsample(srvx,"Petras",degre);
setrandvar_buildsample( srvu , srvx );
// 4. Create the chaos polynomial
pc = polychaos_new ( srvx , 1 );
polychaos_setdegree(pc,degre);
// 5. Perform the experiments
np = setrandvar_getsize(srvu);
mprintf("Number of experiments = %d\n",np);
polychaos_setsizetarget(pc,np);
inputdata = setrandvar_getsample(srvu);
outputdata = nisp_ishigami(inputdata,a,b);
polychaos_settarget(pc,outputdata);
// 6. Compute the coefficients by integration
polychaos_computeexp(pc,srvx,"Integration");
// 7. Get the sensitivity indices
average = polychaos_getmean(pc);
var = polychaos_getvariance(pc);
mprintf("Mean = %f (expected=%f)\n",average,exact.expectation);
mprintf("Variance = %f (expected=%f)\n",var,exact.var);
S = polychaos_getindexfirst (pc);
mprintf("S1 = %f (expected=%f)\n",S(1),exact.S1);
mprintf("S2 = %f (expected=%f)\n",S(2),exact.S2);
mprintf("S3 = %f (expected=%f)\n",S(3),exact.S3);
ST = polychaos_getindextotal (pc);
mprintf("ST1 = %f (expected=%f)\n",ST(1),exact.ST1);
mprintf("ST2 = %f (expected=%f)\n",ST(2),exact.ST2);
mprintf("ST3 = %f (expected=%f)\n",ST(3),exact.ST3);
// Consider the group {X1,X2}
groupe = [1 3];
polychaos_setgroupempty ( pc );
polychaos_setgroupaddvar ( pc , groupe(1) );
polychaos_setgroupaddvar ( pc , groupe(2) );
mprintf("Group {X1,X3}\n");
ST13=exact.S1+exact.S3+exact.S13;
mprintf(" ST13 =%f (expected=%f)\n",polychaos_getgroupind(pc),ST13);
mprintf(" S13 =%f (expected=%f)\n",polychaos_getgroupinter(pc),exact.S13);
// Compute all sensitivity indices.
// The line [1 0 1] is associated with the group {X1,X3}.
// The sum of all sensitivity indices is equal to 1.
mprintf("Anova decomposition of the normalized variance.\n");
polychaos_getanova(pc);
// Compute the histogram of the output of the chaos polynomial.
polychaos_buildsample(pc,"Lhs",10000,0);
sample_output = polychaos_getsample(pc);
my_handle = scf(10001);
histplot(50,sample_output)
xtitle("Ishigami - Histogram");
// Plot the sensitivity indices.
for i=1:nx
indexfirst(i)=polychaos_getindexfirst(pc,i);
indextotal(i)=polychaos_getindextotal(pc,i);
end
my_handle = scf(10002);
bar(indextotal,0.2,'blue');
bar(indexfirst,0.15,'yellow');
legend(["Total" "First order"],pos=1);
xtitle("Ishigami - Sensitivity index");
//
// Clean-up
//
randvar_destroy ( rvu1 );
randvar_destroy ( rvu2 );
randvar_destroy ( rvu3 );
setrandvar_destroy ( srvu );
polychaos_destroy ( pc );
setrandvar_destroy ( srvx );The previous script produces the following output.
Number of experiments = 751
Mean = 3.500000 (expected=3.500000)
Variance = 13.842473 (expected=13.844588)
S1 = 0.313953 (expected=0.313905)
S2 = 0.442325 (expected=0.442411)
S3 = 0.000000 (expected=0.000000)
ST1 = 0.557675 (expected=0.557589)
ST2 = 0.442326 (expected=0.442411)
ST3 = 0.243721 (expected=0.243684)
Group {X1,X3}
ST13 =0.557674 (expected=0.557589)
S13 =0.243721 (expected=0.243684)
Anova decomposition of the normalized variance.
1 0 0 : 0.313953
0 1 0 : 0.442325
1 1 0 : 1.55229e-009
0 0 1 : 8.08718e-031
1 0 1 : 0.243721
0 1 1 : 8.43733e-031
1 1 1 : 1.6007e-007The previous script produces the following plots.
NISP v2.2
The main changes for this release are the following.
Lots of bugs have been fixed including the management of the internal indices, the use of the meta-model and the computeoutput / getoutput functions, the "SmolyakGauss" sampling.
- The help pages have been improved. All the functions are now correctly sorted, the output arguments are clearly documented, the size of the input and output arguments are clarified, the examples are now correct, all the support functions are documented. The help pages are now generated with the automation tool help_from_sci, which allows a fast update of the help pages.
- The User's Manual is improved. The the theory for the first order nonlinear sensitivity indices are added, the theory for product example is inserted.
- Some functions were removed for simplification: the polychaos_computeoutput, getoutput, setinput functions are replaced by polychaos_compute, the propagate input function is removed.
- New functions were created: the support functions nisp_getpath, nisp_destroyall, nisp_print have been created. Probability functions, cumulated probability functions and inverse cumulated probability functions which are used in the module and are not provided by Scilab were added: nisp_exppdf, nisp_expcdf, nisp_expinv, nisp_lognormalpdf, nisp_lognormalcdf, nisp_lognormalinv.
- The module now depends on the apifun and assert modules, which allow to improve the robustness of the functions and simplifies the unit tests.
Here is the detailed changelog for the v2.2:
Fixed getgroupinter: wrong call to GetGroupIndice in the gateway. 12 Nov 2009
- Un-commented the delete operators. Fixed unit tests accordingly (destroy in reverse order of the new)
- Created nisp_destroyall and nisp_print. 8 Dec 2009
- Fixed bug in setinput, computeoutput, getoutput
- Removed propagate input (backward compatibility lost).
- Updated getsample, set sample to use straightforward index management.
- Moved non regression test for getgroupinter into ticket 209.
- Fixed bug in the .xml formatting of polychaos.xml.
- Renamed NISP_library to src.
- Updated building system.
- Added non-regression test for ticket #210 : The computeoutput / getoutput functions do not work.
- Updated formatting of the LaTeX sources of the manual (added licence headers).
- Moved nisp_getpath to a regular macro.
- Restructured the manual.
- Created help pages for nisp configuration functions from pseudo macros.
- Moved directory for help pages of nisp functions.
- Added unit tests for nisp_* functions.
- Fixed description of the randvar_size, setrandvar_size and polychaos_size functions.
- Fixed ticket #222 : the polychaos_getindexfirst function did not work.
- Added non regression test for ticket #222.
- Updated licence headers.
Fixed ticket #211 : the "SmolyakGauss" sampling did not work.
- Fixed ticket #230 : The example in the help of polychaos was uncomplete.
- Fixed ticket #216 : The functions in the polychaos help page were not sorted.
Added examples for LhsMinMax design.
- Added nisp_buildlhs for a classical Lhs design.
- Added erfcinv, lognormalpdf, lognormalcdf, lognormalinv.
- Added exppdf, expcdf, expinv.
- Added comparison between naive method, Homma-Saltelli and Chaos Polynomials.
- Added unit test for nisp_erfcinv.
- Updated nisp_erfcinv for robustness.
- Added missing help pages for support functions.
- Added basic tests for nisp_exppdf, nisp_expcdf, nisp_expinv.
- Added basic tests for nisp_lognormalpdf, nisp_lognormalcdf, nisp_lognormalinv.
- Fixed ticket #212: the polychaos_readtarget function was undefined.
Fixed ticket #235 : a part of nisp_util.cpp moved in SetRandomVariable.
- Removed nisp_index.cpp: moved into nisp_util.cpp
- Removed nisp_random.cpp: moved into nisp_util.cpp
- Moved nisp_gc.cpp into nisp_pc.cpp.
- Removed nisp_gc.cpp.
- Moved smolyak building into thirdparties.
- Fixed ticket #229: the polychaos_setanova function was useless.
- Fixed ticket #217: the polychaos_computeoutput, getoutput, setinput functions are replaced by polychaos_compute. (Backward compatibility lost)
- Forgot to commit nisp_formqua.cpp and nisp_formqua.h
- Updated the manual: filled the theory for the first order nonlinear sensitivity indices.
- Updated the manual: filled the theory for product example.
- Forgot to commit sci_polychaos_compute.cpp
- Created separate directory for test functions.
- Created polychaos/ sub-directory for tests.
- Moved polychaos tests into separate directory.
- Disabled (pragma) the optimization for the Sobol sequence: failed to compile on some Windows.
- Updated unit tests: using functions from the assert module.
NISP v2.3
The main changes for this release are the following.
- Performances of the examples and demos are improved, based on vectorization. This makes us benefit from Scilab's multi-core capabilities. We can perform 100 000 simulations in less than a second.
- The help page are clearer and better organized.
- New functions have been created: nisp_buildlhs, nisp_corrcoef, nisp_cov.
- The User's Manual has been entirely updated. A complete chapter is now devoted to Sensitivity Analysis, with an introduction to the theory and many examples.
- New test functions have been created: nisp_ishigami, nisp_ishigamisa, nisp_product, nisp_productsa, nisp_sum, nisp_sumsa. We provide both the test function and the exact expectation, variance and first and total sensitivity indices. This simplifies the testing of the module.
- New methods for sensitivity analysis were included: a brute-force algorithm and the Sobol method for sensitivity indices. These functions are vectorized, which improves the performance.
- The module now depends on the Specfun module, which simplifies the Sobol method for sensitivity analysis.
The demonstrations were updated and new demos were created: the module now provides 25 demonstrations, which present the use of chaos polynomials, SRC indices, brute force and the Sobol method on 3 functions including the sum function, the product function, the Ishigami function. We included a Scilab port of the AxialStressedBeam demo of Open Turns.
Here is the detailed changelog for the v2.3:
- Updated manual: filled linear model.
- Fixed ticket #252: The nisp_exppdf unit test did not pass on Linux.
- Fixed ticket #254: The ticket_209 test does not pass on Linux.
- Fixed ticket #256 : The polychaos_getquantwilks test in polychaos failed.
- Fixed ticket #255: The polychaos_getquantile made the polychaos unit test fail.
- Fixed ticket #253: The polychaos1 unit test did not pass on Linux (getgroupind).
- Update manual: filled first part of ishigami analysis.
- Update manual: computed the variance of the ishigami model.
- Update manual: moved light png figures into pdf.
- Added comments in the demos.
- Vectorized all scripts.
- Added the ishigami demo with Sobol method.
- Update manual: added ishigami analysis with Sobol method.
- Created nisp_cov and nisp_corrcoef.
- Updated manual: added a section on the LHS samplings.
- Added demos on the LHS designs.
- Simplified the demos for the setrandvar samplings.
- Updated manual: added a separate note on performance.
- Updated manual: filled decomposition of the variance.
- Updated manual: filled sensitivity analysis of Ishigami.
- Updated demos: computed exact sensitivity indices of Ishigami.
- Added nisp_sobolsa, the Sobol method for sensitivity analysis.
- Updated manual: filled numerical results for Ishigami.
- Updated manual: partial fill for Sobol method of SA.
- Updated manual: filled particular case for product.
- Added 3 test functions and their exact SI: Sum, Product, Ishigami.
- Updated the tests, the demos according to the new functions.
- Improved robustness of the test functions.
- Created a bruteforce Sensitivity Indices function.
- Disabled (pragma) the optimization for the Sobol sequence: failed to compile on some Windows.
- Split the tests for polychaos.
- Fixed ticket #237: The unit test dataset for erfcinv was wrong.
- Fixed ticket #249: The parameter b of a normal variable in randvar was unclear.
- Improved robustness and flexibility of lognormal and exponential PDF, CDF and Inv.
- Fixed ticket #228: The output arguments of many methods of polychaos were undocumented.
- Created separated unit test directory for support functions.
- Created unit tests for nisp_buildlhs, nisp_corrcoef, nisp_cov.
- Fixed bug in nisp_corrcoef.
- Organized demos.
- Updated the Axial Stressed Beam demo.
Authors
- Copyright (C) 2008 - 2009 - INRIA - Michaël Baudin
- Copyright (C) 2009 - 2011 - DIGITEO - Michaël Baudin
- Copyright (C) 2009 - CEA-Saclay, DEN, DM2S, SFME - Jean-Marc Martinez
- Copyright (C) 2012 - Michael Baudin