FitVarianceTrend.hpp

Go to the documentation of this file.

#ifndef SCRAN_FIT_TREND_VAR_H
#define SCRAN_FIT_TREND_VAR_H
 
#include "../utils/macros.hpp"
 
#include <algorithm>
#include <vector>
#include "WeightedLowess/WeightedLowess.hpp"
 
namespace scran {
 
class FitVarianceTrend {
public:
    struct Defaults {
        static constexpr double minimum_mean = 0.1;
 
        static constexpr bool filter = true;
 
        static constexpr bool transform = true;
 
        static constexpr double span = 0.3;
 
        static constexpr bool use_fixed_width = false;
 
        static constexpr double fixed_width = 1;
 
        static constexpr int minimum_window_count = 200;
    };
 
public:
    FitVarianceTrend& set_span(double s = Defaults::span) {
        span = s;
        return *this;
    }
 
    FitVarianceTrend& set_minimum_mean(double m = Defaults::minimum_mean) {
        min_mean = m;
        return *this;
    }
 
    FitVarianceTrend& set_filter(bool f = Defaults::filter) {
        filter = f;
        return *this;
    }
 
    FitVarianceTrend& set_transform(bool t = Defaults::transform) {
        transform = t;
        return *this;
    }
 
    FitVarianceTrend& set_use_fixed_width(bool u = Defaults::use_fixed_width) {
        use_fixed_width = u;
        return *this;
    }
 
    FitVarianceTrend& set_fixed_width(double f = Defaults::fixed_width) {
        fixed_width = f;
        return *this;
    }
 
    FitVarianceTrend& set_minimum_window_count(int c = Defaults::minimum_window_count) {
        minimum_window_count = c; 
        return *this;
    }
 
private:
    double span = Defaults::span;
    double min_mean = Defaults::minimum_mean;
    bool filter = Defaults::filter;
    bool transform = Defaults::transform;
 
    bool use_fixed_width = Defaults::use_fixed_width;
    bool fixed_width = Defaults::fixed_width;
    int minimum_window_count = Defaults::minimum_window_count;
 
    static double quad(double x) {
        return x*x*x*x;
    }
 
public:
    void run(size_t n, const double* mean, const double* variance, double* fitted, double* residuals) const {
        std::vector<double> xbuffer(n), ybuffer(n);
 
        size_t counter = 0;
        for (size_t i = 0; i < n; ++i) {
            if (!filter || mean[i] >= min_mean) {
                xbuffer[counter] = mean[i];
                if (transform) {
                    ybuffer[counter] = std::pow(variance[i], 0.25); // Using the same quarter-root transform that limma::voom uses.
                } else {
                    ybuffer[counter] = variance[i];
                }
                ++counter;
            }
        }
 
        if (counter < 2) {
            throw std::runtime_error("not enough observations above the minimum mean");
        }
 
        // Determining the left edge. This needs to be done before
        // run_in_place, which mutates the input xbuffer.
        size_t left_index = std::min_element(xbuffer.begin(), xbuffer.begin() + counter) - xbuffer.begin();
        double left_x = xbuffer[left_index];
 
        WeightedLowess::WeightedLowess<> smoother;
        if (use_fixed_width) {
            smoother.set_span(minimum_window_count);
            smoother.set_span_as_proportion(false);
            smoother.set_min_width(fixed_width);
        } else {
            smoother.set_span(span);
        }
 
        std::vector<double> fbuffer(counter), rbuffer(counter);
        smoother.run(counter, xbuffer.data(), ybuffer.data(), NULL, fbuffer.data(), rbuffer.data());
 
        // Identifying the left-most fitted value.
        double left_fitted = (transform ? quad(fbuffer[left_index]) : fbuffer[left_index]);
 
        counter = 0;
        for (size_t i = 0; i < n; ++i) {
            if (!filter || mean[i] >= min_mean) {
                fitted[i] = (transform ? quad(fbuffer[counter]) : fbuffer[counter]);
                ++counter;
            } else {
                fitted[i] = mean[i] / left_x * left_fitted; // draw a y = x line to the origin from the left of the fitted trend.
            }
            residuals[i] = variance[i] - fitted[i];
        }
        return;
    }
 
public:
    struct Results {
        Results() {}
 
        Results(size_t n) : fitted(n), residuals(n) {}
        std::vector<double> fitted;
 
        std::vector<double> residuals;
    };
 
    Results run(size_t n, const double* mean, const double* variance) const {
        Results output(n);
        run(n, mean, variance, output.fitted.data(), output.residuals.data());
        return output;
    }
};
 
}
 
#endif