Skip to content

Commit

Permalink
Refactor FunctionalChaosValidation based on MetaModelValidation.
Browse files Browse the repository at this point in the history
  • Loading branch information
@mbaudin47
mbaudin47 committed Jul 1, 2023
1 parent a2938eb commit 00fcddc
Show file tree
Hide file tree
Showing 10 changed files with 181 additions and 201 deletions.
2 changes: 1 addition & 1 deletion lib/etc/openturns.conf.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -782,7 +782,7 @@
<!-- OT::LinearModelAnalysis parameters -->
<LinearModelAnalysis-Identifiers value_int="3" />

<!-- OT::LinearModelValidation parameters -->
<!-- OT::LinearModelValidation parameters -->
<LinearModelValidation-DefaultKFoldParameter value_int="10" />

<!-- OT::LinearModelStepwiseAlgorithm parameters -->
Expand Down
8 changes: 8 additions & 0 deletions lib/src/Uncertainty/Algorithm/MetaModel/FunctionalChaos/FunctionalChaosResult.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -270,6 +270,14 @@ Function FunctionalChaosResult::getComposedMetaModel() const
return composedMetaModel_;
}

/** Residuals accessor */
Sample FunctionalChaosResult::getSampleResiduals() const
{
const Sample predictionsSample(metaModel_(inputSample_));
const Sample residualsSample(outputSample_ - predictionsSample);
return residualsSample;
}

/* Method save() stores the object through the StorageManager */
void FunctionalChaosResult::save(Advocate & adv) const
{
Expand Down
227 changes: 91 additions & 136 deletions lib/src/Uncertainty/Algorithm/MetaModel/FunctionalChaos/FunctionalChaosValidation.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,14 +42,26 @@ FunctionalChaosValidation::FunctionalChaosValidation()
}

/* Default constructor */
FunctionalChaosValidation::FunctionalChaosValidation(const FunctionalChaosResult & functionalChaosResult)

FunctionalChaosValidation::FunctionalChaosValidation(const FunctionalChaosResult & functionalChaosResult,
const CrossValidationMethod method,
const UnsignedInteger & kParameter)
: PersistentObject()
, functionalChaosResult_(functionalChaosResult)
, isInitialized_(false)
{
// Nothing to do
if ((method != LEAVEONEOUT) and (method != KFOLD))
throw InvalidArgumentException(HERE) << "The method " << method << " is not available.";
cvMethod_ = method;
if (kParameter < 1)
throw InvalidArgumentException(HERE) << "Cannot set k parameter of K-Fold method to " << kParameter << " which is lower than 1";
const UnsignedInteger sampleSize = functionalChaosResult_.getInputSample().getSize();
if (kParameter > sampleSize)
throw InvalidArgumentException(HERE) << "Cannot set k parameter of K-Fold method to " << kParameter
<< " which is larger than the sample size =" << sampleSize;
kParameter_ = kParameter;
}


/* Virtual constructor */
FunctionalChaosValidation * FunctionalChaosValidation::clone() const
{
Expand All @@ -65,114 +77,64 @@ String FunctionalChaosValidation::__repr__() const
return oss;
}

/* Initialize */
void FunctionalChaosValidation::initialize()
{
if (!isInitialized_)
{
const Sample metamodelPredictions(computeMetamodelCrossValidationPredictions());
Sample outputSample(functionalChaosResult_.getOutputSample());
validation_ = MetaModelValidation(outputSample, metamodelPredictions); // True for LOO. And K-Fold ?
isInitialized_ = true;
}
}

/* Get result*/
FunctionalChaosResult FunctionalChaosValidation::getResult() const
{
return functionalChaosResult_;
}

/* Compute LOO mean squared error */
Point FunctionalChaosValidation::computeLeaveOneOutMeanSquaredError() const
/* Compute mean squared error */
Point FunctionalChaosValidation::computeMeanSquaredError()
{
const UnsignedInteger sampleSize = functionalChaosResult_.getInputSample().getSize();
const UnsignedInteger outputDimension = functionalChaosResult_.getOutputSample().getDimension();
Point leaveOneOutMSE(outputDimension);
Sample leaveOneOutResiduals(computeLeaveOneOutResiduals());
for (UnsignedInteger j = 0; j < outputDimension; ++j)
{
Sample squaredCorrectedResiduals(sampleSize, 1);
for (UnsignedInteger i = 0; i < sampleSize; ++i)
squaredCorrectedResiduals(i, 0) = std::pow(leaveOneOutResiduals(i, j), 2);
leaveOneOutMSE[j] = squaredCorrectedResiduals.computeMean()[0];
} // For output marginal indices
return leaveOneOutMSE;
initialize();
return validation_.computeMeanSquaredError();
}

/* Compute K-Fold mean squared error */
Point FunctionalChaosValidation::computeKFoldMeanSquaredError(const UnsignedInteger & kParameter) const
/* Compute R2 score */
Point FunctionalChaosValidation::computeR2Score()
{
const UnsignedInteger sampleSize = functionalChaosResult_.getInputSample().getSize();
const UnsignedInteger outputDimension = functionalChaosResult_.getOutputSample().getDimension();
// Compute K-Fold error
Sample foldSquaredError(kParameter, outputDimension);
Sample residualsKFold(computeKFoldResiduals(kParameter));
Point kFoldMSE(outputDimension);
KFoldSplitter splitter(sampleSize, kParameter);
for (UnsignedInteger foldIndex = 0; foldIndex < kParameter; ++foldIndex)
{
Indices indicesTest;
const Indices indicesTrain(splitter.generate(indicesTest));
const UnsignedInteger foldDimension = indicesTest.getSize();
for (UnsignedInteger j = 0; j < outputDimension; ++j)
{
for (UnsignedInteger i = 0; i < foldDimension; ++i)
foldSquaredError(foldIndex, j) += std::pow(residualsKFold(indicesTest[i], j), 2);
foldSquaredError(foldIndex, j) /= foldDimension;
} // For output marginal indices
} // For fold indices
kFoldMSE = foldSquaredError.computeMean();
return kFoldMSE;
initialize();
return validation_.computeR2Score();
}

/* Compute LOO R2 score */
Point FunctionalChaosValidation::computeLeaveOneOutR2Score() const
/* Compute residuals */
Sample FunctionalChaosValidation::getResidualSample()
{
Point meanSquaredError(computeLeaveOneOutMeanSquaredError());
const Sample outputSample(functionalChaosResult_.getOutputSample());
Point r2Score(MetaModelValidation::ComputeR2Score(meanSquaredError, outputSample));
return r2Score;
initialize();
return validation_.getResidualSample();
}

/* Compute K-Fold R2 score */
Point FunctionalChaosValidation::computeKFoldR2Score(const UnsignedInteger & kParameter) const
/* Get residual distribution */
Distribution FunctionalChaosValidation::getResidualDistribution(const Bool smooth)
{
Point meanSquaredError(computeKFoldMeanSquaredError(kParameter));
const Sample outputSample(functionalChaosResult_.getOutputSample());
Point r2Score(MetaModelValidation::ComputeR2Score(meanSquaredError, outputSample));
return r2Score;
initialize();
return validation_.getResidualDistribution(smooth);
}

/* Compute leave-one-out residuals */
Sample FunctionalChaosValidation::computeLeaveOneOutResiduals() const
/* Draw */
GridLayout FunctionalChaosValidation::drawValidation()
{
const Sample inputSample(functionalChaosResult_.getInputSample());
const UnsignedInteger sampleSize = inputSample.getSize();
const FunctionCollection reducedBasis(functionalChaosResult_.getReducedBasis());
const UnsignedInteger reducedBasisSize = reducedBasis.getSize();
if (reducedBasisSize > sampleSize)
throw InvalidArgumentException(HERE) << "Error: the sample size is: " << sampleSize <<
" which is lower than the basis size: " << reducedBasisSize;
const Function transformation(functionalChaosResult_.getTransformation());
const Sample standardSample(transformation(inputSample));
DesignProxy designProxy(standardSample, reducedBasis);
Indices allIndices(reducedBasisSize);
allIndices.fill();
const MatrixImplementation designMatrix(designProxy.computeDesign(allIndices));
// The method name is set to the default one, given by ResourceMap
const String methodName(ResourceMap::GetAsString("LeastSquaresExpansion-DecompositionMethod"));
LeastSquaresMethod leastSquaresMethod(LeastSquaresMethod::Build(methodName, designProxy, allIndices));
leastSquaresMethod.update(Indices(0), allIndices, Indices(0));
const Point diagonalH = leastSquaresMethod.getHDiag();
const Sample outputSample(functionalChaosResult_.getOutputSample());
const Function metamodel(functionalChaosResult_.getMetaModel());
const Sample outputMetamodel(metamodel(inputSample));
const UnsignedInteger outputDimension = outputSample.getDimension();
Sample residualsLOO(sampleSize, outputDimension);
for (UnsignedInteger j = 0; j < outputDimension; ++j)
{
const Sample marginalOutputSample(outputSample.getMarginal(j));
const Point marginalOutputPoint(marginalOutputSample.asPoint());
const Sample outputMetamodelMarginal(outputMetamodel.getMarginal(j));
const Point residuals(marginalOutputPoint - outputMetamodelMarginal.asPoint());
for (UnsignedInteger i = 0; i < sampleSize; ++i)
residualsLOO(i, j) = residuals[i] / (1.0 - diagonalH[i]);
} // For output marginal indices
return residualsLOO;
initialize();
return validation_.drawValidation();
}

/* Compute K-Fold residuals */
Sample FunctionalChaosValidation::computeKFoldResiduals(const UnsignedInteger & kParameter) const

/** Compute cross-validation metamodel predictions */
Sample FunctionalChaosValidation::computeMetamodelCrossValidationPredictions() const
{
const Sample residualsSample(functionalChaosResult_.getSampleResiduals());
const Sample inputSample(functionalChaosResult_.getInputSample());
const UnsignedInteger sampleSize = inputSample.getSize();
const FunctionCollection reducedBasis(functionalChaosResult_.getReducedBasis());
Expand All @@ -188,59 +150,52 @@ Sample FunctionalChaosValidation::computeKFoldResiduals(const UnsignedInteger &
DesignProxy designProxy(standardSample, reducedBasis);
Indices allIndices(reducedBasisSize);
allIndices.fill();
const MatrixImplementation designMatrix(designProxy.computeDesign(allIndices));
// The method name is set to the default one, given by ResourceMap
const String methodName(ResourceMap::GetAsString("LeastSquaresExpansion-DecompositionMethod"));
LeastSquaresMethod leastSquaresMethod(LeastSquaresMethod::Build(methodName, designProxy, allIndices));
leastSquaresMethod.update(Indices(0), allIndices, Indices(0));
const SymmetricMatrix projectionMatrix(leastSquaresMethod.getH());
// Compute K-Fold error
Sample residualsKFold(sampleSize, outputDimension);
Point kFoldMSE(outputDimension);
KFoldSplitter splitter(sampleSize, kParameter);
for (UnsignedInteger foldIndex = 0; foldIndex < kParameter; ++foldIndex)
Sample cvPredictions(sampleSize, outputDimension);
if (cvMethod_ == LEAVEONEOUT)
{
Indices indicesTest;
const Indices indicesTrain(splitter.generate(indicesTest));
const Sample inputTestKFold(inputSample.select(indicesTest));
const Sample outpoutTestKFold(outputSample.select(indicesTest));
const Sample predictionKFold(metamodel(inputTestKFold));
const UnsignedInteger foldDimension = indicesTest.getSize();
SymmetricMatrix projectionKFoldMatrix(foldDimension);
for (UnsignedInteger i1 = 0; i1 < foldDimension; ++i1)
for (UnsignedInteger i2 = 0; i2 < 1 + i1; ++i2)
projectionKFoldMatrix(i1, i2) = projectionMatrix(indicesTest[i1], indicesTest[i2]);
const IdentityMatrix identityMatrix(foldDimension);
const SymmetricMatrix reducedMatrix(identityMatrix - projectionKFoldMatrix);
Matrix multipleRightHandSide(foldDimension, outputDimension);
for (UnsignedInteger j = 0; j < outputDimension; ++j)
for (UnsignedInteger i = 0; i < foldDimension; ++i)
multipleRightHandSide(i, j) = outpoutTestKFold(i, j) - predictionKFold(i, j);
const Matrix correctedKFoldResidualsMatrix(reducedMatrix.solveLinearSystem(multipleRightHandSide));
const Point diagonalH = leastSquaresMethod.getHDiag();
for (UnsignedInteger j = 0; j < outputDimension; ++j)
for (UnsignedInteger i = 0; i < foldDimension; ++i)
residualsKFold(indicesTest[i], j) = correctedKFoldResidualsMatrix(i, j);
} // For fold indices
return residualsKFold;
}


/* Draw model vs metamodel validation graph */
GridLayout FunctionalChaosValidation::drawLeaveOneOutValidation() const
{
// Compute leave-one-out predictions
const UnsignedInteger sampleSize = functionalChaosResult_.getInputSample().getSize();
const Sample outputSample(functionalChaosResult_.getOutputSample());
const UnsignedInteger outputDimension = outputSample.getDimension();
Sample leaveOneOutResiduals(computeLeaveOneOutResiduals());
Sample leaveOneOutPredictions(sampleSize, outputDimension);
for (UnsignedInteger j = 0; j < outputDimension; ++j)
for (UnsignedInteger i = 0; i < sampleSize; ++i)
leaveOneOutPredictions(i, j) = outputSample(i, j) - leaveOneOutResiduals(i, j);

// Draw graph
MetaModelValidation validation(outputSample, leaveOneOutPredictions);
GridLayout grid(validation.drawValidation());
return grid;
for (UnsignedInteger i = 0; i < sampleSize; ++i)
cvPredictions(i, j) = outputSample(i, j) - residualsSample(i, j) / (1.0 - diagonalH[i]);
}
else if (cvMethod_ == KFOLD)
{
const SymmetricMatrix projectionMatrix(leastSquaresMethod.getH());
// Compute K-Fold error
Sample residualsKFold(sampleSize, outputDimension);
KFoldSplitter splitter(sampleSize, kParameter_);
for (UnsignedInteger foldIndex = 0; foldIndex < kParameter_; ++foldIndex)
{
Indices indicesTest;
const Indices indicesTrain(splitter.generate(indicesTest));
const UnsignedInteger foldSize = indicesTest.getSize();
// Take into account for different fold sizes
const Scalar foldCorrection = std::sqrt(float(sampleSize) / (kParameter_ * foldSize));
SymmetricMatrix projectionKFoldMatrix(foldSize);
for (UnsignedInteger i1 = 0; i1 < foldSize; ++i1)
for (UnsignedInteger i2 = 0; i2 < 1 + i1; ++i2)
projectionKFoldMatrix(i1, i2) = projectionMatrix(indicesTest[i1], indicesTest[i2]);
const IdentityMatrix identityMatrix(foldSize);
const SymmetricMatrix reducedMatrix(identityMatrix - projectionKFoldMatrix);
const Sample residualsSampleKFoldTest(residualsSample.select(indicesTest));
Matrix multipleRightHandSide(foldSize, outputDimension);
for (UnsignedInteger j = 0; j < outputDimension; ++j)
for (UnsignedInteger i = 0; i < foldSize; ++i)
multipleRightHandSide(i, j) = residualsSampleKFoldTest(i, j);
const Matrix correctedKFoldResidualsMatrix(reducedMatrix.solveLinearSystem(multipleRightHandSide));
for (UnsignedInteger j = 0; j < outputDimension; ++j)
for (UnsignedInteger i = 0; i < foldSize; ++i)
cvPredictions(indicesTest[i], j) = outputSample(indicesTest[i], j) - foldCorrection * correctedKFoldResidualsMatrix(i, j);
} // For fold indices
}
else
throw InvalidArgumentException(HERE) << "The method " << cvMethod_ << " is not available.";
return cvPredictions;
}

/* Method save() stores the object through the StorageManager */
Expand Down
3 changes: 3 additions & 0 deletions lib/src/Uncertainty/Algorithm/MetaModel/FunctionalChaos/openturns/FunctionalChaosResult.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,6 +100,9 @@ public:
/** Composed meta model accessor */
virtual Function getComposedMetaModel() const;

/** Residuals accessor */
virtual Sample getSampleResiduals() const;

/** Method save() stores the object through the StorageManager */
void save(Advocate & adv) const override;

Expand Down
Loading

0 comments on commit 00fcddc

Please sign in to comment.