Recent posts (max 20) - Browse or Archive for more

OT function decorator

Hi,

Here is an easier way to pass python functions to openturns

First define once this decorator for python functions (could be proposed in a future version):

def FunctionDecorator(n, p):
  class otfunc(OpenTURNSPythonFunction):
      
    def __init__(self, func):
      OpenTURNSPythonFunction.__init__(self, n, p) 
      self.__class__.__name__ = func.__name__
      self._exec = func 
        
  return otfunc

Then apply it to your function by adding this one line:

@FunctionDecorator(2, 1)
def SIMPLEFUNC(X):
    Y = [ X[0] + X[1] ]    
    return Y

Then use it through OpenTURNS as usual:

myFunc = NumericalMathFunction( SIMPLEFUNC )

This is open for discussion ...

J.

Python Function for Finding the best Continuous Distribution among all OT Distributions

Hy hereafter, you will find the code of a file which contains a python function which is very usefull for finding the best distribution from a numerical sample. Its use is :

> FindBestContinuousDistribution(NS) # where NS is a NumericalSample


This function returns the continuous distribution in OT which best fits the emprirical distribution of the NumericalSample NS. Some distributions may have been not tested due to the nature of the data (distributions which requires positive value for example) or due to the estimation of the parameters which could fail (see the warnings which are printed during execution)

ps: I have not find an easy way to share a file, that's why i have made a copy of the whole file ...

import openturns as ot
import numpy as np
import warnings




def ListAllFactoryName():
    """ Return the name of all the DistributionFactory in OT """
    allDist = dir(ot.dist)
    AllFactory = []
    for elt in allDist:
        if ('Factory' in elt) and ('_' not in elt) and ('Histogram' not in elt):
            AllFactory.append(elt)
    return(AllFactory)



def ListAllFactory():
    """ Return a dictionnary with DistributionFactory objects of OT with the distinction
    'AllContinuousFactory' and 'AllDiscreteFactory' """
    result = {}
    AllContinuousFactoryName = []
    AllContinuousFactory =[]
    AllDiscreteFactoryName = []
    AllDiscreteFactory =[]    
    AllFactoryName = ListAllFactoryName()    
    for factoryName in AllFactoryName:
        factoryName = 'ot.' + factoryName
        StrDist = factoryName.replace('Factory','()')
        Dist = eval(StrDist)
        if Dist.isContinuous():
            factoryName = factoryName.replace('Factory','Factory()')            
            AllContinuousFactoryName.append(factoryName)
            factory = eval(factoryName)
            AllContinuousFactory.append(factory)
        else :
            factoryName = factoryName.replace('Factory','Factory()')            
            AllDiscreteFactoryName.append(factoryName)
            factory = eval(factoryName)
            AllDiscreteFactory.append(factory)            
    result['AllContinuousFactoryName']=AllContinuousFactoryName
    result['AllContinuousFactory']=AllContinuousFactory
    result['AllDiscreteFactoryName']=AllDiscreteFactoryName
    result['AllDiscreteFactory']=AllDiscreteFactory
    return(result)


def FindBestContinuousDistribution(NS):
    """ Return the continuous distribution in OT which best fits the emprirical distribution of the NumericalSample NS.
    Some distributions may have been not tested due to the nature of the data (distributions which requires positive value for example)
    or due to the estimation of the parameters which could fail (see the warnings)."""
    print '\n Return the continuous distribution in OT which best fits the emprirical distribution of the NumericalSample NS.'
    print 'Some distributions may have been not tested due to the nature of the data (distributions which requires positive value for example)'
    print 'or due to the estimation of the parameters which could fail (see the warnings).'
    if not isinstance(NS,ot.statistics.NumericalSample):
        warnings.warn( '-- The argument is not NumericalSample')
        return
    if NS.getDimension()>1:
        warnings.warn( '-- The NumericalSample must be of dimension 1')
        return
    AllFactory = ListAllFactory()
    ContinuousFactory = AllFactory['AllContinuousFactory']
    ContinuousFactoryName = AllFactory['AllContinuousFactoryName']

    DistTestedOK =[]
    DistTestedNOK = []
    ListFactoryOK = []
    ListFactoryNOK = []
    TestResu = []
    TestPValue = []    
    i = 0
    
    print('\n-------------- NOT TESTED DISTRIBUTIONS ----------------\n')
    for factory in ContinuousFactory:
        Name = ContinuousFactoryName[i]
        Name = Name.replace('Factory()','')
        Name = Name.replace('ot.','')
        try:
            EstimatedDistribution = factory.build(NS)
            ListFactoryOK.append(factory)
            DistTestedOK.append(Name)
            Test = ot.FittingTest.Kolmogorov(NS,EstimatedDistribution)
            TestResu.append(Test.getBinaryQualityMeasure())
            TestPValue.append(Test.getPValue())
        except Exception as e:
            ListFactoryNOK.append(factory)
            DistTestedNOK.append(Name)
            print(Name + ' Distribution Not tested')
            print(e)
            print('----\n')
            #pass
        i = i+1    
    print(DistTestedNOK)
    print('-----------------------------------------------------------------------\n')
    
    print('\n---------------- TESTED DISTRIBUTIONS  -------------------------------')    
    print(DistTestedOK)
    print('--------------------------------------\n')

    print('\n------- ACCEPTED DISTRIBUTIONS WITH KS TEST (at level 5%) -----------\n')
    DistTestedOK = np.array(DistTestedOK)
    DistAccepted = DistTestedOK[np.array(TestResu)==True]
    TestPValue = np.array(TestPValue)
    DistAcceptedPValue = TestPValue[np.array(TestResu)==True]
    Tab = [DistAccepted,DistAcceptedPValue]
    print 'KS-pvalue' + '\t\t', 'Ditribution' 
    print('------------------------------------------')
    for i in range(len(DistAccepted)):
      print DistAcceptedPValue[i],'\t\t',DistAccepted[i]
    #print(np.array(Tab).transpose())
    print('\n--------------------------------------\n')
    
    
    collectionFactory = ot.DistributionFactoryCollection(ListFactoryOK)
    bestDistributionKolmogorov = ot.FittingTest.BestModelKolmogorov(NS,collectionFactory)
    bestDistributionBIC =  ot.FittingTest_BestModelBIC(NS,collectionFactory)
    print("\n--- Best Distribution with Kolmogorov Test ----------------------------")
    print(bestDistributionKolmogorov)
    print('-------------------------------------------------------------------------\n')    

    print("\n--- Best Distribution with BIC ----------------------------------------")
    print(bestDistributionBIC)
    print('-------------------------------------------------------------------------\n')    


    return bestDistributionBIC

SIMUREX OpenTURNS presentation & démonstration

Bonjour,

Voici les présentations faites sur le thème des incertitudes lors du workshop SIMUREX (Cargèse, avril 2012) dédié à la simulation des performances énergétiques du bâtiment:

  • Méthodologie Incertitudes
  • Présentation d'OpenTURNS
  • Présentation d'un cas d'étude

Le cas d'étude consiste à évaluer les besoins en chauffage d'une maison basse consommation pour un profil de confort donné tout en prenant en compte les incertitudes sur les caractéristiques de la maison et de son environnement.

L'étude est réalisée à l'aide de deux modèles différents de puissance instantanée. On évalue les caractéristiques probabilistes des besoins de chauffage dans différents régimes, mettant ainsi en oeuvre plusieurs approches disponibles dans OpenTURNS.

Le modèle de calcul est joint (exécutable Windows 32 bits), mais une approximation de haute fidélité est également disponible pour une utilisation sous Linux.

Tous les scripts et les supports sont attachés à ce post (cliquer sur le titre du post!)

L'étude s'appuie en partie sur une base de données du modèle de calcul qui peut être récupérée ici:  http://dl.free.fr/hq9a0jMTO

Le mot de passe est SIMUREX_OPENTURNS, le délai de garde est de 30j donc ne tardez pas!

Anne D.

OT Plotting capabilities using Matplotlib

Hi,

The work on a new module has begun to plot OpenTURNS graphs through matplotlib.

Those who want a preview of the capabilities can check it out here: otmpl

This development is not finalized yet, but we welcome any comments.
To test it, just download this single python script and run it, it will show the following examples.
Beware, it requires a decent version of matplotlib (>= 1.0) and numpy (=> 1.4 for contours)

Otherwise you can copy it somewhere at PYTHONPATH reach (e.g. /usr/local/lib/python2.6/dist-packages) where it can be used as a module as in:

import openturns as ot
from otmpl import View

# create a graph
graph = ot.Normal().drawPDF()

# create a View, accepted keywords are figure.add_subplot+axes.plot/axes.pie/axes.bar/axes.contour+axes.clabel/axes.step ones depending on the kind of graph
view = View(graph, color='blue')

# save to a file, accepted keywords are figure.savefig ones
view.save('curve.png', dpi=100)

# show the graph on the screenn, accepted keywords are pyplot.show ones
view.show(block=False)


Here are a few samples:
/raw-attachment/blog/Matplotlib%20module/curve1.png /raw-attachment/blog/Matplotlib%20module/curve2.png
/raw-attachment/blog/Matplotlib%20module/curve3.png /raw-attachment/blog/Matplotlib%20module/curve4.png
/raw-attachment/blog/Matplotlib%20module/curve5.png /raw-attachment/blog/Matplotlib%20module/curve6.png
/raw-attachment/blog/Matplotlib%20module/curve7.png /raw-attachment/blog/Matplotlib%20module/curve8.png
/raw-attachment/blog/Matplotlib%20module/curve9.png /raw-attachment/blog/Matplotlib%20module/curve10.png

J.

Fourth OpenTURNS Users Day slides

OpenTURNS Users Day slides are available here http://trac.openturns.org/extra/share/ju_ot_4.tgz .

Module Open TURNS " Integral Compound Poisson Distribution"

Thanks to Anne Dutfoy, a new Open TURNS module is now available. It aims to create a compound Poisson distribution. It is the sum of a Poisson-distributed number of independent identically-distributed random variables.

I invite you to read the documentation shared with the sources of the module to understand the theory and how to use it.

Sources are viewable here: http://trac.openturns.org/browser/openturns-modules/integralcompoundpoisson/trunk

Svn repository is here:  https://svn.openturns.org/openturns-modules/integralcompoundpoisson/trunk

Amélioration du chaos - Postprocessing - « Pretty print » des résultats

Ce travail porte sur une présentation lisible des résultats principaux d’un calcul par chaos polynomial : affichage des moments statistiques de la réponse, coefficients du développement en fonction des indices des polynômes, indices de Sobol’.

La fonction printchaos.py (cf. attachments) permet de réaliser un affichage clair des différents résultats du chaos. Elle contient les fonctions :

  • print_moments : dont l'argument est un objet FunctionalChaosResult. En sortie, la fonction affiche la valeur des moments.
import printchaos as pc
...
polynomialChaosAlgorithm.run()
polynomialChaosResult = polynomialChaosAlgorithm.getResult()
pc.print_moments(polynomialChaosResult)
  • print_coefficients : dont l'argument est un objet FunctionalChaosResult. En option, il est possible de spécifier l'objet ComposedDistribution qui a servi à réaliser le calcul par chaos. En sortie, la fonction affiche le type de polynôme relatif à chaque variable si l'objet ComposedDistribution a été spécifié :
    • He : polynôme d'Hermite
    • Le : polynôme de Legendre
    • Ja : polynôme de Jacobi
    • La : polynôme de Laguerre

La fonction affiche ensuite l'ensemble des degrés des polynômes unidimensionnels de chaque variable formant la base de développement du chaos et la valeur du coefficient associée.

import printchaos as pc
...
polynomialChaosAlgorithm.run()
polynomialChaosResult = polynomialChaosAlgorithm.getResult()
pc.print_coefficients(polynomialChaosResult, PM = ProbaModel)
  • print_plot_sensitivities : dont l'argument est un objet FunctionalChaosResult. En option, il est possible de spécifier l'objet ComposedDistribution qui a servi à réaliser le calcul par chaos si les variables aléatoires ont été nommées (méthode setName). La fonction affiche la liste de tous les indices de Sobol’ groupés par ordre en précisant à quelles variables ils sont associés. Enfin la fonction trace un diagramme circulaire des indices de Sobol’ totaux normalisés.
import printchaos as pc
...
polynomialChaosAlgorithm.run()
polynomialChaosResult = polynomialChaosAlgorithm.getResult()
pc.print_plot_sensitivities(polynomialChaosResult, PM = ProbaModel)

Sobol.png, janv. 2011: /extra/oldshare/AmeliorationChaos/Sobol.png

  • print_plot_regression : dont les arguments sont l’échantillon de sortie (NumericalSample) obtenu par à la main (OT 0.13.2) ou par la méthode getProjectionStrategy().getOutputSample() (OT version supérieure) et l’objet FunctionalChaosResult associé au plan d’expérience. La fonction trace un « scatter plot » des résultats du plan d’expérience avec en abscisse les valeurs du vrai modèle et en ordonnées les valeurs du méta-modèle. Dans le terminal, les coefficients de détermination R2 et Q2 sont donnés.
import printchaos as pc
...
polynomialChaosAlgorithm.run()
polynomialChaosResult = polynomialChaosAlgorithm.getResult()
out = polynomialChaosAlgorithm.getProjectionStrategy().getOutputSample()
pc.print_plot_regression(out, polynomialChaosResult)

Model_vs_MetaModel.png, janv. 2011: /extra/oldshare/AmeliorationChaos/Model_vs_MetaModel.png

Amélioration du chaos - Postprocessing - Obtention de la distribution de la réponse d’un chaos polynomial

Il serait utile de pouvoir obtenir, pour la réponse du chaos, un objet Distribution, soit à partir d'un FunctionalChaosResult, soit d'un FunctionalChaosRandomVector. L'objectif est de pouvoir utiliser la distribution de la sortie comme une entrée pour un autre chaos, ou calcul en général. Cette réponse n’est possible qu'avec un lissage à noyaux (KernelSmoothing).

La fonction kernelchaos.py (cf. attachments) permet de créer la distribution de la variable de sortie développée sur le chaos. Elle contient la fonction getOutputDistribution dont l’argument est un objet FunctionalChaosResult. En sortie la fonction retourne une liste de distribution :

  • la distribution lissée par un noyau gaussien,
  • La meilleure distribution usuelle d’après le test de Kolmogorov si le test est validé.
import kernelchaos as kc
...
polynomialChaosAlgorithm.run()
polynomialChaosResult = polynomialChaosAlgorithm.getResult()
ouputDistributions = kc.getOutputDistribution(polynomialChaosResult)

En outre, elle crée 3 graphiques :

  • Output_distribution_kernel_smoothing.png : présentant un échantillon de la variable de sortie sous la forme d’un histogramme ainsi que le lissage à noyau de l’échantillon

Output_distribution_kernel_smoothing.png|C|Output_distribution_kernel_smoothing.png, janv. 2011: /extra/oldshare/AmeliorationChaos/Output_distribution_kernel_smoothing.png

  • Output_distribution_Kolmogorov_test.png : présentant un échantillon de la variable de sortie sous la forme d’un histogramme ainsi que la distribution usuelle validant au mieux le test de Kolmogorov, dans le cas où le test est validé pour un « factory ».

Output_distribution_Kolmogorov_test.png|C|Output_distribution_Kolmogorov_test.png, janv. 2011: /extra/oldshare/AmeliorationChaos/Output_distribution_Kolmogorov_test.png

  • Output_distribution_KK.png : présentant la comparaison entre la distribution lissée par noyau et la distribution usuelle résultat du test de Kolmogorov si ce dernier est validé pour un « factory ».

Output_distribution_KK.png, janv. 2011 /extra/oldshare/AmeliorationChaos/Output_distribution_KK.png

Présentations de la journée utilisateur #3

Bonjour,

Voici les supports des présentations exposées lors de la troisième journée utilisateur:

  • présentation principale
  • présentation des réseaux bayésiens
  • présentation du module mixmod

Le fichier PDF de la présentation des réseaux bayésiens a été coupé en 6 parties. Faire cat Expose_LIP6_ot2010_BN.pdf.* > Expose_LIP6_ot2010_BN.pdf sous linux, ou utiliser hjsplit sous windows pour reconstruire le fichier

Bonne lecture!

Anne

edit par François : voici le lien pour une version light de l'exposé sur les réseaux bayésiens :

openturns-agrum: a new module for OpenTURNS, that allows to create a distribution from a Bayesian Network using the aGrUM library

Hi, in collaboration with Pierre-Henri WUILLEMIN (LIP6, aGrUM team member) and Anne DUTFOY (EDF R&D, OpenTURNS team member), I have written a new OpenTURNS module that allows to create 1D OpenTURNS distributions as conditional marginal distributions of a Bayesian network built using aGrUM ( http://agrum.lip6.fr/doku.php). It is very useful when one has conditional information about an uncertain quantity but not statistical data directly related to this quantity.

This module is compatible with OpenTURNS 0.13.2 and aGrUM revision 382.

Here an example of a python script to use this module:

from openturns import * from pyAgrum import * from otagrum import *

  1. Create an empty BN

myBN = BayesNet()

  1. Create light variable

light = LabelizedVar("Light", "Light", 0) light.addLabel("Dim") light.addLabel("Bright") lightSize = len(light)

  1. Create moisture vriable

moisture = LabelizedVar("Moisture", "Moisture", 0) moisture.addLabel("Dry") moisture.addLabel("Wet") moistureSize = len(moisture)

  1. Create height variable

height = DiscretizedVar("Height", "Height") data = NumericalPoint(range(10, 100, 10)) heightWhenDimAndDry = Uniform(0.0, 40.0) heightWhenDimAndWet = Uniform(10.0, 60.0) heightWhenBrightAndDry = Uniform(0.0, 20.0) heightWhenBrightAndWet = Uniform(40.0, 100.0) data = BayesNetAgrum.AdaptGrid(DistributionCollection(heightWhenDimAndDry, heightWhenDimA...), data) heightSize = data.getSize() - 1 for i in list(data):

height.addTick(i)

  1. Add variables to the net

indexLight = myBN.add(light) indexMoisture = myBN.add(moisture) indexHeight = myBN.add(height)

  1. Create arcs

myBN.insertArc(indexLight, indexMoisture) myBN.insertArc(indexLight, indexHeight) myBN.insertArc(indexMoisture, indexHeight)

  1. Create conditional probability tables

myBN.cpt(indexLight):= 0.5, 0.5 myBN.cpt(indexMoisture){'Light' : 'Dim'} = 0.2, 0.8 myBN.cpt(indexMoisture){'Light' : 'Bright'} = 0.6, 0.4 myBN.cpt(indexHeight){'Light': 'Dim', 'Moisture': 'Dry'} = BayesNetAgrum().Discretize(heightWhenDimAndDry, data) myBN.cpt(indexHeight){'Light': 'Bright', 'Moisture': 'Dr... = BayesNetAgrum().Discretize(heightWhenBrightAndDry, data) myBN.cpt(indexHeight){'Light': 'Dim', 'Moisture': 'Wet'} = BayesNetAgrum().Discretize(heightWhenDimAndWet, data) myBN.cpt(indexHeight){'Light': 'Bright', 'Moisture': 'We... = BayesNetAgrum().Discretize(heightWhenBrightAndWet, data)

  1. Create a BayesNetAgrum object

otbn = BayesNetAgrum(myBN)

  1. Export to BIF file

otbn.exportToBIFFile("test.bif")

  1. Perform inference

for evidenceValue in "Bright", "Dim":

   # First, set evidence
   otbn.setEvidence("Light", evidenceValue)
   # Second, get a marginal distribution of interrest
   heightDistribution = otbn.getMarginal("Height")
   heightDistribution.setName("Height_Light_" + evidenceValue)
   # Use the distribution the way you want
   print "heightDistribution=", heightDistribution
   heightDistribution.drawPDF(-10.0, 110.0, 512).draw("PDF_" + heightDistribution.getName())
   heightDistribution.drawCDF(-10.0, 110.0, 512).draw("CDF_" + heightDistribution.getName())

New version of Eficas for Open TURNS 0.13.2

Here are tarball of Eficas for use with Open TURNS 0.13.2 or later:

This tool will help you in writing python script to run uncertainty studies with Open TURNS 0.13.2 or later.

You can install it this way:

  • decompress your tarball:

either tar -xzvf EficasV2.0.1.1_OpenTURNS-0.13.2-source.tar.gz

or tar -xjvf EficasV2.0.1.1_OpenTURNS-0.13.2-source.tar.bz2

  • change to newly created directory EficasV2.0.1.1_OpenTURNS-0.13.2

cd EficasV2.0.1.1_OpenTURNS-0.13.2

  • configure the installation (we make the assumption that you plan to install Eficas in <EficasDIR>, to use Python 2.4 and that OpenTURNS is installed in <OpenTURNSDIR>:

cmake . -DCMAKE_INSTALL_PREFIX=<EficasDIR> -DWITH_OPENTURNS=ON -DOPENTURNS_DIR= <OpenTURNSDIR> -DPYTHON_EXECUTABLE=/usr/bin/python2.4

  • build and install:

make make install

  • run Eficas to make wrapper description files:

python2.4 <EficasDIR>/Openturns_Study/qtEficas_openturns_study.py

  • run Eficas to make probabilistic studies:

python2.4 <EficasDIR>/Openturns_Wrapper/qtEficas_openturns_wrapper.py

See the README file for more installation instructions (if you don't have cmake)

Ivan DUTKA-MALEN

Building a metamodel with Polynomial Chaos Expansion using a database

Hello all,

Recently, a colleague of mine asked me if it was possible to build a Polynomial Chaos Expansion metamodel in OpenTURNS with only a database of inputs versus outputs. After a little bit of scratching on the top of my head, I reached Regis and after a small chat with him, I obtained the following script. I have not tested it with such decoration (name of the files, etc.) so maybe it will not work straightforwardly but the core of the problem is treated here.

Enjoy!

# -*- coding: iso-8859-15 -*-
from openturns import *
from math import *
 
#This script allows to build the metamodel by Polynomial Chaos Expansion using a database as input. 
#Y = f(X)
#For the moment, it is mandatory to go through an intermediate variable Z
#The hypothesis is that each variable varies in an interval [a,b]
#Z = phi(X), such that each Zj is Uniform(0.,1.). phi is affine
#Y = g(Z)
#We do the PCE on g and finally f = metamodel o phi
 
#The input file in InputFile.txt
 
 
 
#this function takes a input in the database Z and returns the associated output
def norm2(X,Y):
  norm = 0.
  for i in range(Y.getDimension()):
    norm+=(X[i]-Y[i])**2
  return sqrt(norm)
 
 
#this function takes a input in the database Z and returns the associated output
class statModel(OpenTURNSPythonFunction):
    def __init__(self):
        OpenTURNSPythonFunction.__init__(self, globals()["dim"], 1)
    def f(self,X):
      for i in range(globals()["siz"]):
	if norm2(X,globals()["Z"][i])<1.e-6:
	  v=i
	  break
      try:
	Y = globals()["spOutputTR"][v][0]
      except:
	Y = 0.
      return [Y]
 
 
 
#this function is the affine funtion that transforms X into Z
class affineTransform(OpenTURNSPythonFunction):
	def __init__(self):
	  OpenTURNSPythonFunction.__init__(self, globals()["dim"], globals()["dim"])
	  #change input dimension
	def f(self,X):
	  interval = globals()["interv"]
	  Z = []
	  for j in range(len(X)):
	    Z.append(X[j] - interval[j][0])/(interval[j][1]-interval[j][0])
	  return Z
 
 
 
 
# Get training input from InputFile.txt
inVar = NumericalSample.ImportFromCSVFile("InputFile.txt")
inVar.setName("Input Variable")
 
 
#get the dimension and size of the sample
dim = inVar.getDimension()
siz = inVar.getSize()
 
# Define the intervals where the input variables are defined
#list of min values for each variable, to be specified by user
a1 = 0.
a2 = 0.
a3 = 0.
a4 = 0.
Min = [a1,a2,a3,a4]
#list of max values for each variable, to be specified by user
b1 = 1.
b2 = 2.
b3 = 3.
b4 = 4.
Max = [b1,b2,b3,b4]
 
interv = []
for i in range(len(Min)):
  interv.append([Min[i],Max[i]])
 
# Create Z the transformed input vector
Z = NumericalSample(siz,dim)
for i in range(siz):
  for j in range(dim):
    Z[i][j] = (inVar[i][j] - interv[j][0])/(interv[j][1]-interv[j][0])
 
# Get training output from OutputFile.txt
outVar = NumericalSample.ImportFromCSVFile("OutputFile.txt")
outVar.setName("Output Variable")
 
# Create Polynomial Chaos Expansion of outVar = g(Z)
polyColl = PolynomialFamilyCollection(dim)
 
legendreFamily = LegendreFactory()
for i in range(dim):
  polyColl[i] = OrthogonalUniVariatePolynomialFamily(legendreFamily)
 
multivariateBasis = OrthogonalProductPolynomialFactory(polyColl,EnumerateFunction(dim))
maximumConsideredTerms = 500
mostSignificant = 10
significanceFactor = 1.e-6
truncatureBasisStrategy = CleaningStrategy(OrthogonalBasis(multivariateBasis),maximumConsideredTerms, mostSignificant, significanceFactor, True)
evaluationCoeffStrategy = LeastSquaresStrategy(FixedExperiment(Z))
model = NumericalMathFunction(statModel())
polynomialChaosAlgorithm = FunctionalChaosAlgorithm(model, Distribution(dist), AdaptiveStrategy(truncatureBasisStrategy), ProjectionStrategy(evaluationCoeffStrategy))
 
polynomialChaosAlgorithm.run()
polynomialChaosResult = polynomialChaosAlgorithm.getResult()
 
coefficients = polynomialChaosResult.getCoefficients()
metaModel = polynomialChaosResult.getMetaModel()
 
# Get Response Surface : composition between Y = g(Z) and Z = aff(X) => our model is (g o aff)
proj = NumericalMathFunction(affineTransform())
modelfull = NumericalMathFunction(metaModel,proj)
 
# Get estimation input from InputCheck.txt
inpCheck = NumericalSample.ImportFromCSVFile("InputCheck.txt")
inpCheck.setName("Input variable for checking")
 
 
# Get estimation output
outCheck = modelfull(inpCheck)
outCheck.exportToCSVFile("OutputCheck.csv"," ")
 
fOutputCheck = open("OutputCheck.txt", 'w')
 
for i in range(outCheck.getSize()):
    fOutputCheck.write("%f\n" % (outCheck[i])[0])
 
fOutputCheck.close()

ot-mixmod: a new OT-module to create a Mixture Factory composed of Normal Distributions

Hi, I have written a small OpenTURNS module to create a Mixture Factory composed of Normal Distributions from a Numerical Sample. The sources of the module are available in the openturns-module repository. However, you can download here the version 0.1.0.

This module uses MixMod with the default parameters. This module requires version 0.13.2 of OpenTURNS and MixMod executable.

Here an example of a python script to use this module:

from openturns import * from otmixmod import * draw = True dim = 2 size = 2000 coll = DistributionCollection(0) R = CorrelationMatrix(dim) for i in range(dim-1):

R[i,i+1] = 0.1

coll.add(Distribution(Normal(NumericalPoint(dim, -2.0), NumericalPoint(dim, 1.2), R))) for i in range(dim-1):

R[i,i+1] = -0.2

coll.add(Distribution(Normal(NumericalPoint(dim, 2.0), NumericalPoint(dim, 0.8), R))) mean = NumericalPoint(dim, 5.0) mean[0] = -5.0 coll.add(Distribution(Normal(mean, NumericalPoint(dim, 1.4), CorrelationMatrix(dim)))) mixture = Mixture(coll) sample = mixture.getNumericalSample(size) factory = MixtureFactory(3) estimatedDistribution = factory.buildImplementation(sample) print "Estimated distribution=", estimatedDistribution if draw:

g = estimatedDistribution.drawPDF() c = Cloud(sample) c.setColor("red") c.setPointStyle("bullet") ctmp = g.getDrawable(0) g.setDrawable(Drawable(c), 0) g.addDrawable(ctmp) g.draw("testMixMod")

Automatic quantification of a marginal based on a 1D sample

Hello all,

I propose you the following python script. It allows to simply quantify a 1D variable from a sample. First, you create the family of parametric distribution on which you want to test your sample. If the best distribution does not fulfill the quality criteria, you go to a non parametric quantification by Kernel Smoothing. The output of the function is the distribution and the value of the test.

# -*- coding: iso-8859-15 -*-
from openturns import *
 
def quantification(sample1D):
  # Creation of a factory collection
  factoryCollection = FactoryCollection()
  factoryCollection.add(DistributionFactory(UniformFactory()))
  factoryCollection.add(DistributionFactory(NormalFactory()))
  factoryCollection.add(DistributionFactory(BetaFactory()))
  #factoryCollection.add(DistributionFactory(WeibullFactory()))
  #factoryCollection.add(DistributionFactory(LogNormalFactory()))
  #factoryCollection.add(DistributionFactory(BetaFactory()))
  #factoryCollection.add(DistributionFactory(GammaFactory()))
  #factoryCollection.add(DistributionFactory(ExponentialFactory()))
  #factoryCollection.add(DistributionFactory(GumbelFactory()))
  #factoryCollection.add(DistributionFactory(WeibullFactory()))
  #factoryCollection.add(DistributionFactory(TruncatedNormalFactory()))
  #factoryCollection.add(DistributionFactory(StudentFactory()))
 
  try:
    BestDistributionMarg = Distribution(FittingTest().BestModelKolmogorov(sample1D, factoryCollection))
    print BestDistributionMarg
  except:
    print "impossible to do the parametric quantification"
## Study of the result by comparing the empirical law of the sample and the distribution obtained by parametric
## if the test is ok, we keep the parametric law, if not, go to non parametric quantification
  try:
    TestMarg = FittingTest_Kolmogorov(SampleInputMarg,BestDistributionMarg).getBinaryQualityMeasure()
    print "Parametric quantification"
  except:
    TestMarg = False
## test not ok
  if TestMarg == False :
    kernel=KernelSmoothing()
    BestDistributionMarg = Distribution(kernel.buildImplementation(sample1D))
    print "Quantification by kernel smoothing"
  ## optional visual tests
  QQPlot = VisualTest().DrawQQplot(sample1D,BestDistributionMarg,1000)
  QQPlot.draw("QQPlotMargin", 640, 480)
  QQPlot.getBitmap()
 
  return BestDistributionMarg, TestMarg

Annonce workshop OPUS le 25 novembre 2009 à Suresnes

Dans le cadre du projet ANR OPUS, un workshop intitulé "Méthodes spectrales et Chaos polynomial" se tiendra le 25 novembre 2009 chez EADS à Suresnes. Seront traités au cours de cette journée aussi bien des aspects théoriques (le matin) que des problématiques d'implémentations logicielles (l'après-midi).

A partir de 9h30, Olivier Lemaitre (LIMSI) et Fabio Nobile (Politecnico di Milano) feront deux présentations sur la théorie. L'après midi, Géraud Blatman et Thierry Crestaux, deux doctorants, présenteront leurs travaux de thèse en rapport avec le sujet. Enfin, sur les problèmes d'implémentations informatiques, Michael Bodin (Scilab) et Jean-Marc Martinez (CEA) présenteront leur toolbox Scilab NISP ; Marc Berveiller (EDF) et Régis Lebrun (EADS) présenteront ensuite les développements sur le chaos fonctionnel dans OpenTURNS.

Un plan d'accès à EADS Suresnes est disponible. Pour des raisons de sécurité inhérente à l'hôte de ce workshop, il est obligatoire de se pré-inscrire auprès de Jayant Sen Gupta (jayant.sengupta_at_eads_dot_net) en précisant le nom, le prénom et la société. De plus, le jour du workshop, il sera demandé de présenter une pièce d'identité (carte d'identité ou passeport).

Eficas is a Python/QT4 GUI to ease writing uncertainty studies and wrapper description files

Here is a tarball of Eficas

This tool will help you in writing python script to run uncertainty studies with Open TURNS 0.13.2 or later. See the README file for installation instructions.

Ivan

How to use LHSExperiment ?

I want to use the method LHSExperiment from OpenTURNS version 0.13.1.

I use LHSExperiment(distribution, size) under this form but I have an error (wrong number of argument) and I don't understand my mistake.

How to get back a model from a saved study?

I want to obtain a model that was saved in a Study Object and I do not know how to do it. The use case guide indicates how to do it for a NumericalPoint but not for a model. I do not achieve to derive it for a NumericalMathFunction. Thanks for your help

Fabien

Présentation "Chaos Polynômial" de la deuxième journée utilisateur

Courte introduction au chaos polynômial et à sa mise en oeuvre dans Open TURNS

Planches de la journée utilisateur #2 à EADS Innovation Works