bool CMultivariatePrior::probabilityOfLessLikelySamples()

in lib/maths/common/CMultivariatePrior.cc [107:222]
100 lines of code
18 McCabe index (conditional complexity)

bool CMultivariatePrior::probabilityOfLessLikelySamples(maths_t::EProbabilityCalculation calculation,
                                                        const TDouble10Vec1Vec& samples,
                                                        const TDouble10VecWeightsAry1Vec& weights,
                                                        const TSize10Vec& coordinates,
                                                        TDouble10Vec2Vec& lowerBounds,
                                                        TDouble10Vec2Vec& upperBounds,
                                                        TTail10Vec& tail) const {
    if (coordinates.empty()) {
        lowerBounds.clear();
        upperBounds.clear();
        tail.clear();
        return true;
    }

    lowerBounds.assign(2, TDouble10Vec(coordinates.size(), 0.0));
    upperBounds.assign(2, TDouble10Vec(coordinates.size(), 1.0));
    tail.assign(coordinates.size(), maths_t::E_UndeterminedTail);

    if (samples.empty()) {
        LOG_ERROR(<< "Can't compute distribution for empty sample set");
        return false;
    }
    if (!this->check(samples, weights)) {
        return false;
    }
    if (this->isNonInformative()) {
        lowerBounds.assign(2, TDouble10Vec(coordinates.size(), 1.0));
        upperBounds.assign(2, TDouble10Vec(coordinates.size(), 1.0));
        return true;
    }

    using TDouble1Vec = core::CSmallVector<double, 1>;
    using TDoubleWeightsAry1Vec = maths_t::TDoubleWeightsAry1Vec;
    using TJointProbabilityOfLessLikelySamplesVec =
        core::CSmallVector<CJointProbabilityOfLessLikelySamples, 10>;

    static const TSize10Vec NO_MARGINS;
    static const TSizeDoublePr10Vec NO_CONDITIONS;

    TJointProbabilityOfLessLikelySamplesVec lowerBounds_[2]{
        TJointProbabilityOfLessLikelySamplesVec(coordinates.size()),
        TJointProbabilityOfLessLikelySamplesVec(coordinates.size())};
    TJointProbabilityOfLessLikelySamplesVec upperBounds_[2]{
        TJointProbabilityOfLessLikelySamplesVec(coordinates.size()),
        TJointProbabilityOfLessLikelySamplesVec(coordinates.size())};

    std::size_t d = this->dimension();
    TSize10Vec marginalize(d - 1);
    TSizeDoublePr10Vec condition(d - 1);
    TDouble1Vec sc(1);
    TDoubleWeightsAry1Vec wc(1);

    for (std::size_t i = 0; i < coordinates.size(); ++i) {
        std::size_t coordinate = coordinates[i];

        std::copy_if(boost::make_counting_iterator(std::size_t(0)),
                     boost::make_counting_iterator(d), marginalize.begin(),
                     [coordinate](std::size_t j) { return j != coordinate; });
        TUnivariatePriorPtr margin(this->univariate(marginalize, NO_CONDITIONS).first);
        if (!margin) {
            return false;
        }

        for (std::size_t j = 0; j < samples.size(); ++j) {
            for (std::size_t k = 0, l = 0; k < d; ++k) {
                if (k != coordinate) {
                    condition[l++] = std::make_pair(k, samples[j][k]);
                }
            }

            sc[0] = samples[j][coordinate];
            for (std::size_t k = 0; k < weights[j].size(); ++k) {
                wc[0][k] = weights[j][k][coordinate];
            }

            double lb[2], ub[2];
            maths_t::ETail tc[2];

            if (!margin->probabilityOfLessLikelySamples(calculation, sc, wc,
                                                        lb[0], ub[0], tc[0])) {
                LOG_ERROR(<< "Failed to compute probability for coordinate " << coordinate);
                return false;
            }
            LOG_TRACE(<< "lb(" << coordinate << ") = " << lb[0] << ", ub("
                      << coordinate << ") = " << ub[0]);

            TUnivariatePriorPtr conditional(
                this->univariate(NO_MARGINS, condition).first);
            if (!conditional->probabilityOfLessLikelySamples(calculation, sc, wc,
                                                             lb[1], ub[1], tc[1])) {
                LOG_ERROR(<< "Failed to compute probability for coordinate " << coordinate);
                return false;
            }
            LOG_TRACE(<< "lb(" << coordinate << "|.) = " << lb[1] << ", ub("
                      << coordinate << "|.) = " << ub[1]);

            lowerBounds_[0][i].add(lb[0]);
            upperBounds_[0][i].add(ub[0]);
            lowerBounds_[1][i].add(lb[1]);
            upperBounds_[1][i].add(ub[1]);
            tail[i] = static_cast<maths_t::ETail>(tail[i] | tc[1]);
        }
    }

    for (std::size_t i = 0; i < coordinates.size(); ++i) {
        if (!lowerBounds_[0][i].calculate(lowerBounds[0][i]) ||
            !upperBounds_[0][i].calculate(upperBounds[0][i]) ||
            !lowerBounds_[1][i].calculate(lowerBounds[1][i]) ||
            !upperBounds_[1][i].calculate(upperBounds[1][i])) {
            LOG_ERROR(<< "Failed to compute probability for coordinate " << coordinates[i]);
            return false;
        }
    }

    return true;
}