convergence_tools.cpp

Go to the documentation of this file.

/*
 *  SPDX-FileCopyrightText: Copyright 2021, Siavash Ameli <sameli@berkeley.edu>
 *  SPDX-License-Identifier: BSD-3-Clause
 *  SPDX-FileType: SOURCE
 *
 *  This program is free software: you can redistribute it and/or modify it
 *  under the terms of the license found in the LICENSE.txt file in the root
 *  directory of this source tree.
 */
 
 
// =======
// Headers
// =======
 
#include "./convergence_tools.h"
#include <cmath>  // std::sqrt, std::abs, INFINITY, NAN, isnan
#include <algorithm>  // std::max
#include "./special_functions.h"  // erf_inv
#include "../_c_arithmetics/c_arithmetics.h"  // c_arithmetics
 
 
// =================
// check convergence
// =================
 
 
template <typename DataType>
FlagType ConvergenceTools<DataType>::check_convergence(
        DataType** samples,
        const IndexType min_num_samples,
        const IndexType num_inquiries,
        const IndexType* processed_samples_indices,
        const IndexType num_processed_samples,
        const DataType confidence_level,
        const DataType error_atol,
        const DataType error_rtol,
        DataType* error,
        IndexType* num_samples_used,
        FlagType* converged)
{
    FlagType all_converged;
    IndexType j;
 
    // If number of processed samples are not enough, set to not converged yet.
    // This is essential since in the first few iterations, the standard
    // deviation of the cumulative averages are still too small.
    if (num_processed_samples < min_num_samples)
    {
        // Skip computing error. Fill outputs with trivial initial values
        for (j=0; j < num_inquiries; j++)
        {
            error[j] = INFINITY;
            converged[j] = 0;
            num_samples_used[j] = num_processed_samples;
        }
        all_converged = 0;
        return all_converged;
    }
 
    IndexType i;
    DataType summand;
    DataType mean;
    DataType mean_abs;
    DataType std;
    DataType mean_discrepancy;
    DataType data;
 
    // Quantile of normal distribution (usually known as the "z" coefficient)
    DataType standard_z_score = std::sqrt(2) * \
        static_cast<DataType>(erf_inv(static_cast<double>(confidence_level)));
 
    // For each column of samples, compute error of all processed rows
    for (j=0; j < num_inquiries; ++j)
    {
        // Do not check convergence if j-th column already converged
        if (converged[j] == 0)
        {
            // mean of j-th column using all processed rows of j-th column
            summand = 0.0;
            for (i=0; i < num_processed_samples; ++i)
            {
                summand += samples[processed_samples_indices[i]][j];
            }
            mean = summand / num_processed_samples;
 
            // mean of absolute values of j-th column using all processed rows
            // of j-th column
            summand = 0.0;
            for (i=0; i < num_processed_samples; ++i)
            {
                summand += std::abs(samples[processed_samples_indices[i]][j]);
            }
            mean_abs = summand / num_processed_samples;
 
            // If mean of absolute values is zero, do not scale data with it as
            // it causes divide by zero. In this case, all sample data are
            // zero, and there is no need to scale them.
            DataType zero = 0.0;
            if (c_arithmetics::is_equal(mean_abs, zero))
            {
                mean_abs = 1.0;
            }
 
            // std of j-th column using all processed rows of j-th column
            if (num_processed_samples > 1)
            {
                summand = 0.0;
                for (i=0; i < num_processed_samples; ++i)
                {
                    data = samples[processed_samples_indices[i]][j];
                    mean_discrepancy = data - mean;
 
                    // Normalize to the mean of absolute values to avoid
                    // underflow and overflow, but later, re-scale it back. The
                    // underflow or overflow is caused by taking the square two
                    // lines below.
                    mean_discrepancy /= mean_abs;
                    
                    summand += mean_discrepancy * mean_discrepancy;
                }
                std = std::sqrt(summand / (num_processed_samples - 1.0));
 
                // Re-scale back by mean of absolute values
                std *= mean_abs;
            }
            else
            {
                std = INFINITY;
            }
 
            // Compute error based of std and confidence level
            error[j] = standard_z_score * std / \
                std::sqrt(num_processed_samples);
 
            // Check error with atol and rtol to find if j-th column converged
            if (error[j] < std::max(error_atol, error_rtol*mean))
            {
                converged[j] = 1;
            }
 
            // Update how many samples used so far to average j-th column
            num_samples_used[j] = num_processed_samples;
        }
    }
 
    // Check convergence is reached for all columns (all inquiries)
    all_converged = 1;
    for (j=0; j < num_inquiries; ++j)
    {
        if (converged[j] == 0)
        {
            // The j-th column not converged.
            all_converged = 0;
            break;
        }
    }
 
    return all_converged;
}
 
 
// =================
// average estimates
// =================
 
 
template <typename DataType>
void ConvergenceTools<DataType>::average_estimates(
        const DataType confidence_level,
        const DataType outlier_significance_level,
        const IndexType num_inquiries,
        const IndexType max_num_samples,
        const IndexType* num_samples_used,
        const IndexType* processed_samples_indices,
        DataType** samples,
        IndexType* num_outliers,
        DataType* trace,
        DataType* error)
{
    IndexType i;
    IndexType j;
    DataType summand;
    DataType mean;
    DataType mean_abs;
    DataType std;
    DataType mean_discrepancy;
    DataType outlier_half_interval;
 
    // Flag which samples are outliers
    FlagType* outlier_indices = new FlagType[max_num_samples];
 
    // Quantile of normal distribution (usually known as the "z" coefficient)
    DataType error_z_score = std::sqrt(2) * erf_inv(confidence_level);
 
    // Confidence level of outlier is the complement of significance level
    DataType outlier_confidence_level = 1.0 - outlier_significance_level;
 
    // Quantile of normal distribution area where is not considered as outlier
    DataType outlier_z_score = std::sqrt(2.0) * \
        erf_inv(outlier_confidence_level);
 
    for (j=0; j < num_inquiries; ++j)
    {
        // Initialize outlier indices for each column of samples
        for (i=0; i < max_num_samples; ++i)
        {
            outlier_indices[i] = 0;
        }
        num_outliers[j] = 0;
 
        // Compute mean of the j-th column
        summand = 0.0;
        for (i=0; i < num_samples_used[j]; ++i)
        {
            summand += samples[processed_samples_indices[i]][j];
        }
        mean = summand / num_samples_used[j];
 
        // Compute mean of the absolute values of j-th column
        summand = 0.0;
        for (i=0; i < num_samples_used[j]; ++i)
        {
            summand += std::abs(samples[processed_samples_indices[i]][j]);
        }
        mean_abs = summand / num_samples_used[j];
 
        // If mean of absolute values is zero, do not scale data with it as it
        // causes divide by zero. In this case, all sample data are zero, and
        // there is no need to scale them.
        DataType zero = 0.0;
        if (c_arithmetics::is_equal(mean_abs, zero))
        {
            mean_abs = 1.0;
        }
 
        // Compute std of the j-th column
        if (num_samples_used[j] > 1)
        {
            summand = 0.0;
            for (i=0; i < num_samples_used[j]; ++i)
            {
                mean_discrepancy = \
                    samples[processed_samples_indices[i]][j] - mean;
 
                // Normalize to the mean of absolute values to avoid underflow
                // and overflow, but later, re-scale it back. The underflow or
                // overflow is caused by taking the square two lines below.
                mean_discrepancy /= mean_abs;
 
                summand += mean_discrepancy * mean_discrepancy;
            }
            std = std::sqrt(summand / (num_samples_used[j] - 1.0));
 
            // Re-scale back by mean of absolute values
            std *= mean_abs;
        }
        else
        {
            std = INFINITY;
        }
 
        // Outlier half interval
        outlier_half_interval = outlier_z_score * std;
 
        // Find outliers by the difference of each element from mean
        for (i=0; i < num_samples_used[j]; ++i)
        {
            mean_discrepancy = samples[processed_samples_indices[i]][j] - mean;
            if (std::abs(mean_discrepancy) > outlier_half_interval)
            {
                // Outlier detected
                outlier_indices[i] = 1;
                num_outliers[j] += 1;
            }
        }
 
        // Re-evaluate mean but leave out outliers
        summand = 0.0;
        for (i=0; i < num_samples_used[j]; ++i)
        {
            if (outlier_indices[i] == 0)
            {
                summand += samples[processed_samples_indices[i]][j];
            }
        }
        mean = summand / (num_samples_used[j] - num_outliers[j]);
 
        // Re-evaluate std but leave out outliers
        if (num_samples_used[j] > 1 + num_outliers[j])
        {
            summand = 0.0;
            for (i=0; i < num_samples_used[j]; ++i)
            {
                if (outlier_indices[i] == 0)
                {
                    mean_discrepancy = \
                        samples[processed_samples_indices[i]][j] - mean;
 
                // Normalize to the mean of absolute values to avoid underflow
                // and overflow, but later, re-scale it back. The underflow or
                // overflow is caused by taking the square two lines below.
                mean_discrepancy /= mean_abs;
 
                    summand += mean_discrepancy * mean_discrepancy;
                }
            }
            std = std::sqrt(
                summand/(num_samples_used[j] - num_outliers[j] - 1.0));
 
            // Re-scale back by mean of absolute values
            std *= mean_abs;
        }
        else
        {
            std = INFINITY;
        }
 
        // trace and its error
        trace[j] = mean;
        error[j] = error_z_score * std / \
            std::sqrt(num_samples_used[j] - num_outliers[j]);
    }
 
    delete[] outlier_indices;
}
 
 
// ===============================
// Explicit template instantiation
// ===============================
 
template class ConvergenceTools<float>;
template class ConvergenceTools<double>;
template class ConvergenceTools<long double>;