1 //                                               -*- C++ -*-
2 /**
3  *  @brief The MixedHistogramUserDefined distribution
4  *
5  *  Copyright 2005-2021 Airbus-EDF-IMACS-ONERA-Phimeca
6  *
7  *  This library is free software: you can redistribute it and/or modify
8  *  it under the terms of the GNU Lesser General Public License as published by
9  *  the Free Software Foundation, either version 3 of the License, or
10  *  (at your option) any later version.
11  *
12  *  This library is distributed in the hope that it will be useful,
13  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  *  GNU Lesser General Public License for more details.
16  *
17  *  You should have received a copy of the GNU Lesser General Public License
18  *  along with this library.  If not, see <http://www.gnu.org/licenses/>.
19  *
20  */
21 
22 #include <cmath>
23 #include "openturns/MixedHistogramUserDefined.hxx"
24 #include "openturns/ComposedDistribution.hxx"
25 #include "openturns/Dirac.hxx"
26 #include "openturns/Histogram.hxx"
27 #include "openturns/Uniform.hxx"
28 #include "openturns/DistFunc.hxx"
29 #include "openturns/OSS.hxx"
30 #include "openturns/Tuples.hxx"
31 #include "openturns/PersistentObjectFactory.hxx"
32 #include "openturns/ResourceMap.hxx"
33 #include "openturns/RandomGenerator.hxx"
34 
35 
36 BEGIN_NAMESPACE_OPENTURNS
37 
38 CLASSNAMEINIT(MixedHistogramUserDefined);
39 
40 static const Factory<MixedHistogramUserDefined> Factory_MixedHistogramUserDefined;
41 
42 /* Default constructor */
MixedHistogramUserDefined()43 MixedHistogramUserDefined::MixedHistogramUserDefined()
44   : DistributionImplementation()
45   , ticksCollection_(1, Point(1))
46   , kind_(1, DISCRETE)
47   , probabilityTable_(1, 1.0)
48   , discreteIndices_(1, 0)
49   , allIndices_(1, 1, Indices(1, 0))
50   , normalizedProbabilityTable_(1, 1.0)
51 {
52   setName("MixedHistogramUserDefined");
53   setDimension(1);
54   computeRange();
55 }
56 
57 /* Parameters constructor */
MixedHistogramUserDefined(const PointCollection & ticksCollection,const Indices & kind,const Point & probabilityTable)58 MixedHistogramUserDefined::MixedHistogramUserDefined(const PointCollection & ticksCollection,
59     const Indices & kind,
60     const Point & probabilityTable)
61   : DistributionImplementation()
62   , ticksCollection_(ticksCollection)
63   , kind_(kind)
64   , probabilityTable_(probabilityTable)
65 {
66   setName("MixedHistogramUserDefined");
67   const UnsignedInteger dimension = kind.getSize();
68   // Check the ticks
69   if (ticksCollection.getSize() != dimension) throw InvalidArgumentException(HERE) << "Error: expected a collection of ticks of size=" << dimension << ", got size=" << ticksCollection.getSize();
70   // Check the probability table
71   // kind[i] == 0 -> the ith marginal is discrete
72   // kind[i] == 1 -> the ith marginal is continuous
73   UnsignedInteger totalSize = 1;
74   for (UnsignedInteger i = 0; i < dimension; ++i)
75   {
76     if (kind[i] > 1)
77       throw InvalidArgumentException(HERE) << "Kind must be in [[0, 1]]";
78     if (!(ticksCollection[i].getSize() >= 1))
79       throw InvalidArgumentException(HERE) << "Empty ticks";
80     if (ticksCollection[i].getSize() == kind[i])
81       throw InvalidArgumentException(HERE) << "Need at least 2 ticks for continuous variable";
82     totalSize *= ticksCollection[i].getSize() - kind[i];
83   }
84   if (probabilityTable.getSize() != totalSize) throw InvalidArgumentException(HERE) << "Error: expected a probability table of size=" << totalSize << ", got size=" << probabilityTable.getSize();
85 
86   // cache some data
87   Indices discretization(dimension);
88   for (UnsignedInteger i = 0; i < dimension; ++i)
89     // Here, kind[i] == 0 <-> kind[i] is false <-> i is discrete
90     discretization[i] = ticksCollection_[i].getSize() - kind_[i];
91   allIndices_ = Tuples(discretization).generate();
92   for (UnsignedInteger j = 0; j < dimension; ++ j)
93     if (kind_[j] == DISCRETE)
94       discreteIndices_.add(j);
95     else
96       continuousIndices_.add(j);
97   Scalar weightSum = 0.0;
98   for (UnsignedInteger i = 0; i < probabilityTable_.getSize(); ++ i)
99     weightSum += probabilityTable_[i];
100   normalizedProbabilityTable_ = probabilityTable_ / weightSum;
101 
102   setDimension(dimension);
103   computeRange();
104 }
105 
106 /* Comparison operator */
operator ==(const MixedHistogramUserDefined & other) const107 Bool MixedHistogramUserDefined::operator ==(const MixedHistogramUserDefined & other) const
108 {
109   if (this == &other) return true;
110   return (ticksCollection_ == other.ticksCollection_) && (kind_ == other.kind_) && (probabilityTable_ == other.probabilityTable_);
111 }
112 
equals(const DistributionImplementation & other) const113 Bool MixedHistogramUserDefined::equals(const DistributionImplementation & other) const
114 {
115   const MixedHistogramUserDefined* p_other = dynamic_cast<const MixedHistogramUserDefined*>(&other);
116   return p_other && (*this == *p_other);
117 }
118 
119 /* String converter */
__repr__() const120 String MixedHistogramUserDefined::__repr__() const
121 {
122   OSS oss(true);
123   oss << "class=" << MixedHistogramUserDefined::GetClassName()
124       << " name=" << getName()
125       << " dimension=" << getDimension()
126       << " ticksCollection=" << ticksCollection_
127       << " kind=" << kind_
128       << " probabilityTable=" << probabilityTable_;
129   return oss;
130 }
131 
__str__(const String & offset) const132 String MixedHistogramUserDefined::__str__(const String & offset) const
133 {
134   OSS oss(false);
135   oss << offset << getClassName() << "(ticksCollection = " << ticksCollection_ << ", kind = " << kind_ << ", probabilityTable = " << probabilityTable_ << ")";
136   return oss;
137 }
138 
139 /* Virtual constructor */
clone() const140 MixedHistogramUserDefined * MixedHistogramUserDefined::clone() const
141 {
142   return new MixedHistogramUserDefined(*this);
143 }
144 
145 /* Compute the numerical range of the distribution given the parameters values */
computeRange()146 void MixedHistogramUserDefined::computeRange()
147 {
148   const UnsignedInteger dimension = getDimension();
149   Point lowerBound(dimension);
150   Point upperBound(dimension);
151   for (UnsignedInteger j = 0; j < dimension; ++j)
152   {
153     const Point ticks(ticksCollection_[j]);
154     lowerBound[j] = ticks[0];
155     upperBound[j] = ticks[0];
156     for (UnsignedInteger k = 1; k < ticks.getSize(); ++ k)
157     {
158       if (ticks[k] < lowerBound[j])
159         lowerBound[j] = ticks[k];
160       if (ticks[k] > upperBound[j])
161         upperBound[j] = ticks[k];
162     }
163   }
164   setRange(Interval(lowerBound, upperBound));
165 }
166 
167 
168 /* Get one realization of the distribution */
getRealization() const169 Point MixedHistogramUserDefined::getRealization() const
170 {
171   const UnsignedInteger dimension = getDimension();
172   if (!base_.getSize())
173     (void) DistFunc::rDiscrete(normalizedProbabilityTable_, base_, alias_);
174   const UnsignedInteger index = DistFunc::rDiscrete(base_, alias_);
175   Point realization(dimension);
176   for (UnsignedInteger j = 0; j < discreteIndices_.getSize(); ++j)
177   {
178     const UnsignedInteger jDiscrete = discreteIndices_[j];
179     const UnsignedInteger k = allIndices_(index, jDiscrete);
180     const Point ticks(ticksCollection_[jDiscrete]);
181     realization[jDiscrete] = ticks[k];
182   }
183   for (UnsignedInteger j = 0; j < continuousIndices_.getSize(); ++j)
184   {
185     const UnsignedInteger jContinuous = continuousIndices_[j];
186     const UnsignedInteger k = allIndices_(index, jContinuous);
187     const Point ticks(ticksCollection_[jContinuous]);
188     realization[jContinuous] = ticks[k] + (ticks[k + 1] - ticks[k]) * RandomGenerator::Generate();
189   }
190   return realization;
191 }
192 
193 /* Get a sample of the distribution */
getSample(const UnsignedInteger size) const194 Sample MixedHistogramUserDefined::getSample(const UnsignedInteger size) const
195 {
196   return DistributionImplementation::getSample(size);
197 }
198 
199 /* Get the PDF of the distribution */
computePDF(const Point & point) const200 Scalar MixedHistogramUserDefined::computePDF(const Point & point) const
201 {
202   const UnsignedInteger dimension = getDimension();
203   if (point.getDimension() != dimension)
204     throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=" << dimension << ", here dimension=" << point.getDimension();
205 
206   // build the list of discrete ticks indices, with early exit if no tick matches
207   Indices discreteTicksIndices(discreteIndices_.getSize());
208   for (UnsignedInteger j = 0; j < discreteIndices_.getSize(); ++ j)
209   {
210     const Point ticks(ticksCollection_[discreteIndices_[j]]);
211     const UnsignedInteger index = ticks.find(point[discreteIndices_[j]]);
212     if (index >= ticks.getSize())
213       return 0.0;
214     discreteTicksIndices[j] = index;
215   }
216 
217   // loop over the cpt
218   Scalar pdfValue = 0.0;
219   const UnsignedInteger totalSize = probabilityTable_.getSize();
220   for (UnsignedInteger i = 0; i < totalSize; ++i)
221   {
222     Bool skip = false;
223     // first, loop over the discrete components
224     for (UnsignedInteger j = 0; j < discreteIndices_.getSize(); ++j)
225     {
226       const UnsignedInteger k = allIndices_(i, discreteIndices_[j]);
227       if (discreteTicksIndices[j] != k)
228       {
229         skip = true;
230         break;
231       }
232     }
233 
234     // exclude non-matching discrete terms
235     if (skip)
236       continue;
237 
238     skip = false;
239 
240     // now compute the pdf over continuous components
241     Scalar pdfI = 1.0;
242     for (UnsignedInteger j = 0; j < continuousIndices_.getSize(); ++j)
243     {
244       const UnsignedInteger k = allIndices_(i, continuousIndices_[j]);
245       const Point ticks(ticksCollection_[continuousIndices_[j]]);
246       const Scalar x = point[continuousIndices_[j]];
247       if ((x <= ticks[k]) || (x > ticks[k + 1]))
248       {
249         skip = true;
250         break;
251       }
252       pdfI *= 1.0 / (ticks[k + 1] - ticks[k]);
253     }
254 
255     // exclude non-matching continuous terms
256     if (skip)
257       continue;
258 
259     pdfValue += normalizedProbabilityTable_[i] * pdfI;
260   }
261   return pdfValue;
262 }
263 
264 
265 /* Get the CDF of the distribution */
computeCDF(const Point & point) const266 Scalar MixedHistogramUserDefined::computeCDF(const Point & point) const
267 {
268   const UnsignedInteger dimension = getDimension();
269   if (point.getDimension() != dimension)
270     throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=" << dimension << ", here dimension=" << point.getDimension();
271 
272   // build the list of discrete ticks, with early exit if no tick matches
273   Indices discreteTicksIndices(discreteIndices_.getSize());
274   for (UnsignedInteger j = 0; j < discreteIndices_.getSize(); ++ j)
275   {
276     const Point ticks(ticksCollection_[discreteIndices_[j]]);
277     const Scalar x = point[discreteIndices_[j]];
278     UnsignedInteger index = ticks.getSize();
279     for (UnsignedInteger i = 0; i < ticks.getSize(); ++i)
280       if (ticks[i] <= x)
281         index = i;
282     if (index >= ticks.getSize())
283       return 0.0;
284     discreteTicksIndices[j] = index;
285   }
286 
287   // loop over the cpt
288   Scalar cdfValue = 0.0;
289   const UnsignedInteger totalSize = probabilityTable_.getSize();
290   for (UnsignedInteger i = 0; i < totalSize; ++i)
291   {
292     Bool skip = false;
293     // first, loop over the discrete components
294     for (UnsignedInteger j = 0; j < discreteIndices_.getSize(); ++j)
295     {
296       const UnsignedInteger k = allIndices_(i, discreteIndices_[j]);
297       if (k > discreteTicksIndices[j])
298       {
299         skip = true;
300         break;
301       }
302     }
303 
304     // exclude non-matching discrete terms
305     if (skip)
306       continue;
307 
308     skip = false;
309 
310     // now compute the cdf over continuous components
311     Scalar cdfI = 1.0;
312     for (UnsignedInteger j = 0; j < continuousIndices_.getSize(); ++j)
313     {
314       const UnsignedInteger k = allIndices_(i, continuousIndices_[j]);
315       const Point ticks(ticksCollection_[continuousIndices_[j]]);
316       const Scalar x = point[continuousIndices_[j]];
317       if (x <= ticks[k])
318       {
319         skip = true;
320         break;
321       }
322       else if (x < ticks[k + 1])
323         cdfI *= (x - ticks[k]) / (ticks[k + 1] - ticks[k]);
324       //else (x>ticks[k + 1]: nothing to do
325     }
326 
327     // exclude non-matching continuous terms
328     if (skip)
329       continue;
330 
331     cdfValue += normalizedProbabilityTable_[i] * cdfI;
332   }
333   return cdfValue;
334 }
335 
computeComplementaryCDF(const Point & point) const336 Scalar MixedHistogramUserDefined::computeComplementaryCDF(const Point & point) const
337 {
338   return DistributionImplementation::computeComplementaryCDF(point);
339 }
340 
341 /* Get the characteristic function of the distribution, i.e. phi(u) = E(exp(I*u*X)) */
computeCharacteristicFunction(const Scalar x) const342 Complex MixedHistogramUserDefined::computeCharacteristicFunction(const Scalar x) const
343 {
344   return DistributionImplementation::computeCharacteristicFunction(x);
345 }
346 
347 /* Get the quantile of the distribution */
computeQuantile(const Scalar prob,const Bool tail) const348 Point MixedHistogramUserDefined::computeQuantile(const Scalar prob,
349     const Bool tail) const
350 {
351   return DistributionImplementation::computeQuantile(prob, tail);
352 }
353 
354 /* Get the i-th marginal distribution */
getMarginal(const UnsignedInteger index) const355 Distribution MixedHistogramUserDefined::getMarginal(const UnsignedInteger index) const
356 {
357   const UnsignedInteger dimension = getDimension();
358   if (index >= dimension) throw InvalidArgumentException(HERE) << "The index of a marginal distribution must be in the range [0, dim-1]";
359   if (dimension == 1) return clone();
360 
361   // contract probability table
362   const Point ticks(ticksCollection_[index]);
363   const UnsignedInteger size = ticks.getSize();
364   Point marginalProbabilityTable((kind_[index] == DISCRETE) ? size : size - 1);
365   const UnsignedInteger totalSize = probabilityTable_.getSize();
366   for (UnsignedInteger i = 0; i < totalSize; ++i)
367   {
368     const UnsignedInteger k = allIndices_(i, index);
369     marginalProbabilityTable[k] += probabilityTable_[i];
370   }
371 
372   Distribution marginal;
373   if (kind_[index] == DISCRETE)
374   {
375     SampleImplementation support(size, 1);
376     support.setData(ticks);
377     marginal = UserDefined(support, marginalProbabilityTable);
378   }
379   else
380   {
381     marginal = Histogram(ticks, marginalProbabilityTable);
382   }
383   marginal.setDescription(Description(1, getDescription()[index]));
384   return marginal;
385 }
386 
387 /* Get the distribution of the marginal distribution corresponding to indices dimensions */
getMarginal(const Indices & indices) const388 Distribution MixedHistogramUserDefined::getMarginal(const Indices & indices) const
389 {
390   const UnsignedInteger dimension = getDimension();
391   if (!indices.check(dimension)) throw InvalidArgumentException(HERE) << "The indices of a marginal distribution must be in the range [0, dim-1] and must be different";
392 
393   Indices full(dimension);
394   full.fill();
395   if (indices == full) return clone();
396 
397   // contract probability table
398   Indices marginalKind;
399   PointCollection marginalTicksCollection;
400   UnsignedInteger marginalTotalSize = 1;
401   Description description(getDescription());
402   Description marginalDescription;
403   Indices discretization;
404   for (UnsignedInteger j = 0; j < indices.getSize(); ++j)
405   {
406     const UnsignedInteger index = indices[j];
407     marginalKind.add(kind_[index]);
408     marginalTicksCollection.add(ticksCollection_[index]);
409     const UnsignedInteger size = ticksCollection_[index].getSize();
410     discretization.add((kind_[index] == DISCRETE) ? size : size - 1);
411     marginalTotalSize *= discretization[j];
412     marginalDescription.add(description[index]);
413   }
414   IndicesCollection marginalAllIndices(Tuples(discretization).generate());
415   Point marginalProbabilityTable(marginalTotalSize);
416   const UnsignedInteger totalSize = probabilityTable_.getSize();
417 
418   // compute the base of discretization to quickly retrieve the global marginal index
419   Indices productDiscretization(indices.getSize(), 1);
420   for (UnsignedInteger j = 1; j < indices.getSize(); ++j)
421   {
422     productDiscretization[j] = productDiscretization[j - 1] * discretization[j - 1];
423   }
424 
425   for (UnsignedInteger i = 0; i < totalSize; ++i)
426   {
427     for (UnsignedInteger j = 0; j < indices.getSize(); ++j)
428     {
429       // find the global marginal index
430       UnsignedInteger marginalProbabilityTableIndex = 0;
431       for (UnsignedInteger k = 0; k < indices.getSize(); ++k)
432       {
433         marginalProbabilityTableIndex += allIndices_(i, indices[k]) * productDiscretization[k];
434       }
435       marginalProbabilityTable[marginalProbabilityTableIndex] += probabilityTable_[i];
436     }
437   }
438 
439   MixedHistogramUserDefined marginal(marginalTicksCollection, marginalKind, marginalProbabilityTable);
440   marginal.setDescription(marginalDescription);
441   return marginal;
442 } // getMarginal(Indices)
443 
444 /* Check if the distribution is continuous */
isContinuous() const445 Bool MixedHistogramUserDefined::isContinuous() const
446 {
447   const UnsignedInteger size = kind_.getSize();
448   for (UnsignedInteger i = 0; i < size; ++i)
449     if (kind_[i] == DISCRETE) return false;
450   return true;
451 }
452 
453 /* Check if the distribution is discrete */
isDiscrete() const454 Bool MixedHistogramUserDefined::isDiscrete() const
455 {
456   const UnsignedInteger size = kind_.getSize();
457   for (UnsignedInteger i = 0; i < size; ++i)
458     if (kind_[i] == CONTINUOUS) return false;
459   return true;
460 }
461 
462 /* Check if the distribution is integral */
isIntegral() const463 Bool MixedHistogramUserDefined::isIntegral() const
464 {
465   const Scalar epsilon = ResourceMap::GetAsScalar("DiscreteDistribution-SupportEpsilon");
466   const UnsignedInteger size = kind_.getSize();
467   for (UnsignedInteger i = 0; i < size; ++i)
468   {
469     if (kind_[i] == CONTINUOUS) return false;
470     const UnsignedInteger supportSize = ticksCollection_[i].getSize();
471     for (UnsignedInteger j = 0; j < supportSize; ++j)
472     {
473       const Scalar x = ticksCollection_[i][j];
474       if (std::abs(x - floor(x + 0.5)) > epsilon) return false;
475     }
476   } // i
477   return true;
478 }
479 
480 
481 /* Compute the mean of the distribution */
computeMean() const482 void MixedHistogramUserDefined::computeMean() const
483 {
484   const UnsignedInteger dimension = getDimension();
485   mean_ = Point(dimension);
486   const UnsignedInteger totalSize = probabilityTable_.getSize();
487   for (UnsignedInteger i = 0; i < totalSize; ++i)
488   {
489     Point meani(dimension);
490     for (UnsignedInteger j = 0; j < dimension; ++j)
491     {
492       const UnsignedInteger k = allIndices_(i, j);
493       const Point ticks(ticksCollection_[j]);
494       if (kind_[j] == DISCRETE)
495         meani[j] = ticks[k];
496       else
497         meani[j] = 0.5 * (ticks[k] + ticks[k + 1]);
498     }
499     mean_ += meani * normalizedProbabilityTable_[i];
500   }
501   isAlreadyComputedMean_ = true;
502 }
503 
504 /* Get the standard deviation of the distribution */
getStandardDeviation() const505 Point MixedHistogramUserDefined::getStandardDeviation() const
506 {
507   const UnsignedInteger dimension = getDimension();
508   Point standardDeviation(dimension);
509   for (UnsignedInteger i = 0; i < dimension; ++i)
510     standardDeviation[i] = getMarginal(i).getStandardDeviation()[0];
511   return standardDeviation;
512 }
513 
514 /* Get the skewness of the distribution */
getSkewness() const515 Point MixedHistogramUserDefined::getSkewness() const
516 {
517   const UnsignedInteger dimension = getDimension();
518   Point skewness(dimension);
519   for (UnsignedInteger i = 0; i < dimension; ++i)
520     skewness[i] = getMarginal(i).getSkewness()[0];
521   return skewness;
522 }
523 
524 /* Get the kurtosis of the distribution */
getKurtosis() const525 Point MixedHistogramUserDefined::getKurtosis() const
526 {
527   const UnsignedInteger dimension = getDimension();
528   Point kurtosis(dimension);
529   for (UnsignedInteger i = 0; i < dimension; ++i)
530     kurtosis[i] = getMarginal(i).getKurtosis()[0];
531   return kurtosis;
532 }
533 
534 /* Compute the covariance of the distribution */
computeCovariance() const535 void MixedHistogramUserDefined::computeCovariance() const
536 {
537   const UnsignedInteger dimension = getDimension();
538   covariance_ = CovarianceMatrix(dimension);
539   for (UnsignedInteger j = 0; j < dimension; ++j)
540     covariance_(j, j) = 0.0;
541   // First, compute E(X.X^t)
542   const UnsignedInteger totalSize = probabilityTable_.getSize();
543   for (UnsignedInteger i = 0; i < totalSize; ++i)
544   {
545     Point meanI(dimension);
546     Point varianceI(dimension);
547     for (UnsignedInteger j = 0; j < dimension; ++j)
548     {
549       const UnsignedInteger k = allIndices_(i, j);
550       const Point ticks(ticksCollection_[j]);
551       if (kind_[j] == DISCRETE)
552         meanI[j] = ticks[k];
553       else
554       {
555         meanI[j] = 0.5 * (ticks[k] + ticks[k + 1]);
556         const Scalar eta = ticks[k + 1] - ticks[k];
557         varianceI[j] = eta * eta / 12.0;
558       }
559     }
560     for(UnsignedInteger row = 0; row < dimension; ++row)
561       for(UnsignedInteger column = 0; column <= row; ++column)
562         covariance_(row, column) += normalizedProbabilityTable_[i] * (((row == column) ? varianceI[row] : 0.0) + meanI[row] * meanI[column]);
563   }
564   // Then, subtract E(X).E(X)^t
565   const Point mean(getMean());
566   for(UnsignedInteger row = 0; row < dimension; ++row)
567     for(UnsignedInteger column = 0; column <= row; ++column)
568       covariance_(row, column) -= mean[row] * mean[column];
569   isAlreadyComputedCovariance_ = true;
570 }
571 
572 /* Get the moments of the standardized distribution */
getStandardMoment(const UnsignedInteger n) const573 Point MixedHistogramUserDefined::getStandardMoment(const UnsignedInteger n) const
574 {
575   const UnsignedInteger dimension = getDimension();
576   Point standardMoment(dimension);
577   for (UnsignedInteger i = 0; i < dimension; ++i)
578     standardMoment[i] = getMarginal(i).getStandardMoment(n)[0];
579   return standardMoment;
580 }
581 
582 /* Get the standard representative in the parametric family, associated with the standard moments */
getStandardRepresentative() const583 Distribution MixedHistogramUserDefined::getStandardRepresentative() const
584 {
585   return clone();
586 }
587 
588 /* Ticks collection accessor */
setTicksCollection(const Collection<Point> & ticksCollection)589 void MixedHistogramUserDefined::setTicksCollection(const Collection<Point> & ticksCollection)
590 {
591   *this = MixedHistogramUserDefined(ticksCollection, kind_, probabilityTable_);
592   computeRange();
593 }
594 
getTicksCollection() const595 Collection<Point> MixedHistogramUserDefined::getTicksCollection() const
596 {
597   return ticksCollection_;
598 }
599 
600 /* Kind accessor */
setKind(const Indices & kind)601 void MixedHistogramUserDefined::setKind(const Indices & kind)
602 {
603   *this = MixedHistogramUserDefined(ticksCollection_, kind, probabilityTable_);
604   computeRange();
605 }
606 
getKind() const607 Indices MixedHistogramUserDefined::getKind() const
608 {
609   return kind_;
610 }
611 
612 /* Probability table accessor */
setProbabilityTable(const Point & probabilityTable)613 void MixedHistogramUserDefined::setProbabilityTable(const Point & probabilityTable)
614 {
615   *this = MixedHistogramUserDefined(ticksCollection_, kind_, probabilityTable);
616   computeRange();
617 }
618 
getProbabilityTable() const619 Point MixedHistogramUserDefined::getProbabilityTable() const
620 {
621   return probabilityTable_;
622 }
623 
624 /* Conversion as a Mixture */
asMixture() const625 Mixture MixedHistogramUserDefined::asMixture() const
626 {
627   const UnsignedInteger dimension = getDimension();
628   const UnsignedInteger totalSize = probabilityTable_.getSize();
629   Mixture mixture;
630   // Special case: dimension 1
631   if (dimension == 1)
632   {
633     const Point ticks(ticksCollection_[0]);
634     if (kind_[0] == DISCRETE)
635     {
636       const UnsignedInteger size = ticks.getSize();
637       SampleImplementation support(size, 1);
638       support.setData(ticks);
639       mixture = Mixture(Mixture::DistributionCollection(1, UserDefined(support, probabilityTable_)));
640     }
641     // Continuous
642     else
643       mixture = Mixture(Mixture::DistributionCollection(1, Histogram(ticks, probabilityTable_)));
644   } // dimension == 1
645   else
646   {
647     Bool allDiscrete = (discreteIndices_.getSize() == dimension);
648     // Multivariate discrete
649     if (allDiscrete)
650     {
651       Sample support(totalSize, dimension);
652       for (UnsignedInteger i = 0; i < totalSize; ++i)
653       {
654         for (UnsignedInteger j = 0; j < dimension; ++j)
655           support(i, j) = ticksCollection_[j][allIndices_(i, j)];
656       }
657       mixture = Mixture(Mixture::DistributionCollection(1, UserDefined(support, probabilityTable_)));
658     } // allDiscrete
659     else
660     {
661       Mixture::DistributionCollection atoms(totalSize);
662       for (UnsignedInteger i = 0; i < totalSize; ++i)
663       {
664         Mixture::DistributionCollection subAtoms(dimension);
665         for (UnsignedInteger j = 0; j < dimension; ++j)
666         {
667           const UnsignedInteger k = allIndices_(i, j);
668           const Point ticks = ticksCollection_[j];
669           if (kind_[j] == DISCRETE)
670             subAtoms[j] = Dirac(Point(1, ticks[k]));
671           // Continuous
672           else
673             subAtoms[j] = Uniform(ticks[k], ticks[k + 1]);
674         } // j
675         atoms[i] = ComposedDistribution(subAtoms);
676       } // i
677       mixture = Mixture(atoms, probabilityTable_);
678     } // At least one continuous
679   } // dimension > 1
680   mixture.setDescription(getDescription());
681   return mixture;
682 }
683 
save(Advocate & adv) const684 void MixedHistogramUserDefined::save(Advocate & adv) const
685 {
686   DistributionImplementation::save(adv);
687   adv.saveAttribute( "ticksCollection_", ticksCollection_ );
688   adv.saveAttribute( "kind_", kind_ );
689   adv.saveAttribute( "probabilityTable_", probabilityTable_ );
690 }
691 
692 /* Method load() reloads the object from the StorageManager */
load(Advocate & adv)693 void MixedHistogramUserDefined::load(Advocate & adv)
694 {
695   DistributionImplementation::load(adv);
696   adv.loadAttribute( "ticksCollection_", ticksCollection_ );
697   adv.loadAttribute( "kind_", kind_ );
698   adv.loadAttribute( "probabilityTable_", probabilityTable_ );
699   computeRange();
700 }
701 
702 /* Description accessor */
setDescription(const Description & description)703 void MixedHistogramUserDefined::setDescription(const Description & description)
704 {
705   DistributionImplementation::setDescription(description);
706 }
707 
708 END_NAMESPACE_OPENTURNS
709