Hi Regis, Thanks for this long and well detailed answer! The code you provided seems to work as expected.
However during my tests I noticed that the memory was not freed correctly. Once the class FunctionalChaosAlgorithm is called, there is a memory bump and even after calling del and gc.collect(), memory is still not freed (using memory_profiler for that). Might be a memory leak? Kind regards, Pamphile ROY Chercheur doctorant en Quantification d’Incertitudes CERFACS - Toulouse (31) - France +33 (0) 5 61 19 31 57 +33 (0) 7 86 43 24 22 > Le 7 oct. 2017 à 19:59, regis lebrun <[email protected]> a > écrit : > > Hi Pamphil, > > You were almost right: the AdaptiveStieltjesAlgorithm is very close to what > you are looking for, but not exactly what you need. It is the algorithmic > part of the factory of orthonormal polynomials, the class you have to use is > StandardDistributionPolynomialFactory, ie a factory (=able to build > something) and not an algorithm (=something able to compute something). You > have all the details here: > > http://openturns.github.io/openturns/master/user_manual/_generated/openturns.StandardDistributionPolynomialFactory.html > > I agree on the fact that the difference is quite subtle, as it can be seen by > comparing the API of the two classes. The distinction was made at a time were > several algorithms were competing for the task (GramSchmidtAlgorithm, > ChebychevAlgorithm) but in fact the AdaptiveStieltjesAlgorithm proved to be > much more accurate and reliable than the other algorithms, and now it is the > only orthonormalization algorithm available. > > Another subtle trick is the following. > > If you create a basis this way: > basis = ot.StandardDistributionPolynomialFactory(dist) > you will get the basis associated to the *standard representative* > distribution in the parametric family to which dist belongs. It means the > distribution with zero mean and unit variance, or with support equals to [-1, > 1], or dist itself if no affine transformation is able to reduce the number > of parameters of the distribution. > It is managed automatically within the FunctionalChaosAlgorithm, but can be > disturbing if you do things by hand. > > If you create a basis this way: > basis = > ot.StandardDistributionPolynomialFactory(ot.AdaptiveStieltjesAlgorithm(dist)) > then the distribution is preserved, and you get the orthonormal polynomials > corresponding to dist. Be aware of the fact that the algorithm may have hard > time to build the polynomials if your distribution is far away from its > standard representative, as it may involves the computation of recurrence > coefficients with a much wider range of variation. The benefit is that the > orthonormality measure is exactly your distribution, assuming that its copula > is the independent one, so you don't have to introduce a marginal > transformation between both measures. > > Some additional remarks: > + it looks like dist has dimension>1, as you extract its marginal > distributions later on. AdaptiveStieltjesAlgorithm and > StandardDistributionPolynomialFactory only work with 1D distributions (it is > not checked by the library, my shame). What you have to do is: > > basis = > ot.OrthogonalProductPolynomialFactory([ot.StandardDistributionPolynomialFactory(ot.AdaptiveStieltjesAlgorithm(dist.getMarginal(i))) > for i in range(dist.getDimension())]) > Quite a long line, I know... > It will build a multivariate polynomial basis orthonormal with respect to the > product distribution (ie with independent copula) sharing the same 1D > marginal distributions as dist. > > > After that, everything will work as expected and you will NOT have to build > the transformation (if you build it it will coincide with the identity > function). If you encounter performance issues (the polynomials of high > degrees take ages to be built as in http://trac.openturns.org/ticket/885, or > there is an overflow, or the numerical precision is bad) then use: > basis = > ot.OrthogonalProductPolynomialFactory([ot.StandardDistributionPolynomialFactory(dist.getMarginal(i)) > for i in range(dist.getDimension())]) > and build the transformation the way you do it. > > + if you use the FunctionalChaosAlgorithm class by providing an input sample > and an output sample, you also have to provide the weights of the input > sample EVEN IF the experiment given in the projection strategy would allow to > recompute them. It is because the fact that you provide the input sample > overwrite the weighted experiment of the projection stratey by a > FixedExperiment doe. > > I attached two complete examples: one using the exact marginal distributions > and the other using the standard representatives. > > Best regards > > Régis > > ________________________________ > De : roy <[email protected]> > À : regis lebrun <[email protected]> > Cc : users <[email protected]> > Envoyé le : Vendredi 6 octobre 2017 14h22 > Objet : Re: [ot-users] Sample transformation > > > > Hi Regis, > > Thank you for this detailed answer. > > - I am using the latest release from conda (OT 1.9, python 3.6.2, latest > numpy, etc.) , > - For the sample, I need it to generate externally the output (cost code that > cannot be integrated into OT as model), > - I have to convert ot.Sample into np.array because it is then used by other > functions to create the simulations, etc. > > If I understood correctly, I can create the projection strategy using this > snippet: > > basis = ot.AdaptiveStieltjesAlgorithm(dist) > measure = basis.getMeasure() > quad = ot.Indices(in_dim) > for i in range(in_dim): > quad[i] = degree + 1 > > comp_dist = ot.GaussProductExperiment(measure, quad) > proj_strategy = ot.IntegrationStrategy(comp_dist) > > inv_trans = > ot.Function(ot.MarginalTransformationEvaluation([measure.getMarginal(i) for i > in range(in_dim)], distributions)) > sample = np.array(inv_trans(comp_dist.generate())) > > > It seems to work. Except that the basis does not work with > ot.FixedStrategy(basis, dim_basis). I get a non implemented method error. > > After I get the sample and the corresponding output, what is the way to go? > Which arguments do I need to use for the > ot.FunctionalChaosAlgorithm? > > I am comparing the Q2 and on Ishigami and I was only able to get correct > results using: > > pc_algo = ot.FunctionalChaosAlgorithm(sample, output, dist, trunc_strategy) > > But for least square strategy I had to use this: > > pc_algo = ot.FunctionalChaosAlgorithm(sample, output) > > > Is it normal? > > > Pamphile ROY > Chercheur doctorant en Quantification d’Incertitudes > CERFACS - Toulouse (31) - France > +33 (0) 5 61 19 31 57 > +33 (0) 7 86 43 24 22 > > > > Le 5 oct. 2017 à 15:40, regis lebrun <[email protected]> a > écrit : >> >> Hi Pamphile, >> >> >> 1) The problem: >> The problem you get is due to the fact that in your version of OpenTURNS >> (1.7 I suppose), the GaussProductExperiment class has a different way to >> handle the input distribution than the other WeightedExperiment classes: it >> generates the quadrature rule of the *standard representatives* of the >> marginal distributions instead of the marginal distributions. It does not >> change the rate of convergence of the PCE algorithm and allows to use >> specific algorithms for distributions with known orthonormal polynomials. It >> is not explained in the documentation and if you ask the doe for its >> distribution it will give you the initial distribution instead of the >> standardized one. >> >> 2) The mathematical background: >> The generation of quadrature rules for arbitrary 1D distributions is a badly >> conditioned problem. Even if the quadrature rule is well-defined (existence >> of moments of any order, distribution characterized by these moments), the >> application that maps the recurrence coefficients of the orthogonal >> polynomials to their value can have a very large condition number. As a >> result, the adaptive integration used to compute the recurrence coefficients >> of order n, based on the values of the polynomials of degree n-1 and n-2, >> can lead to wrong values and all the process falls down. >> >> 3) The current state of the software: >> Since version 1.8, OpenTURNS no more generates the quadrature rule of the >> standard representatives, but the quadrature rule of the actual marginal >> distributions. The AdaptiveStieltjesAlgorithm class, introduced in release >> 1.8, is much more robust than the previous orthonormalization algorithms and >> is able to handle even stiff problems. There are still difficult situations >> (distributions with discontinuous PDF inside of the range, fixed in OT 1.9, >> or really badly conditioned distributions, hopefully fixed when ticket#885 >> will be solved) but most usual situations are under control even with >> marginal degrees of order 20. >> >> 4) The (probable) bug in your code and the way to solve it >> You must be aware of the fact that the distribution you put into your >> WeightedExperiment object will be superseded by the distribution >> corresponding to your OrthogonalBasisFactory inside of the >> FunctionalChaosAlgorithm. If you need to have the input sample before to run >> the functional chaos algorithm, then you have to build your transformation >> by hand. Assuming that you already defined your projection basis called >> 'myBasis', your marginal integration degrees 'myDegrees' and your marginal >> distributions 'myMarginals', you have to write (in OT 1.7): >> >> # Here the explicit cast into a NumericalMathFunction is to be able to >> evaluate the transformation over a sample >> myTransformation = >> ot.NumericalMathFunction(ot.MarginalTransformationEvaluation([myBasis.getDistribution().getMarginal(i) >> for i in range(dimension), myMarginals)) >> sample = >> myTransformation(ot.GaussProductExperiment(myBasis.getDistribution(), >> myDegrees).generate()) >> >> >> You should avoid to cast OT objects into np objects as much as possible, and >> if you cannot avoid these casts you should do them only in the sections >> where they are needed. They can be expansive for large objects, and if the >> sample you get from generate() is used only as an argument of a >> NumericalMathFunction, then it will be converted back into a NumericalSample! >> >> Best regards >> >> Régis >> ________________________________ >> De : roy <[email protected]> >> À : users <[email protected]> >> Envoyé le : Jeudi 5 octobre 2017 11h13 >> Objet : [ot-users] Sample transformation >> >> >> >> Hi, >> >> I am facing consistency concerns in the API regarding distributions and >> sampling. >> >> The initial goal was to get the sampling for Polynomial Chaos as I must not >> use the model variable. >> So for least square strategy I do something like this: >> >> proj_strategy = ot.LeastSquaresStrategy(montecarlo_design) >> sample = np.array(proj_strategy.getExperiment().generate()) >> >> sample is correct as the bounds of each feature lie in the corresponding >> ranges. >> >> But now if I want to use IntegrationStrategy: >> >> ot.IntegrationStrategy(ot.GaussProductExperiment(dists, list)) >> sample = np.array(proj_strategy.getExperiment().generate()) >> >> sample’s outputs lie between [-1, 1] which does not corresponds to the >> distribution I have initially. >> >> So I used the conversion class but it does not work well with >> GaussProductExperiment as it requires [0, 1] instead of [-1, 1]. >> >> Thus I use this hack: >> >> # Convert from [-1, 1] -> input distributions >> marg_inv_transf = ot.MarginalTransformationEvaluation(distributions, 1) >> sample = (proj_strategy.getExperiment().generate() + 1) / 2. >> >> >> Is it normal that the distribution classes are not returning in the same >> intervals? >> >> >> Thanks for your support! >> >> >> Pamphile ROY >> Chercheur doctorant en Quantification d’Incertitudes >> CERFACS - Toulouse (31) - France >> +33 (0) 5 61 19 31 57 >> +33 (0) 7 86 43 24 22 >> >> >> _______________________________________________ >> OpenTURNS users mailing list >> [email protected] >> http://openturns.org/mailman/listinfo/users > <example.py><example_standard.py>
_______________________________________________ OpenTURNS users mailing list [email protected] http://openturns.org/mailman/listinfo/users
