Method run

src/modules/stats/robust_variance_coxph.cpp:146–229

// ----------------------------------------------------------------------

AnyType rb_coxph_step_transition::run(AnyType& args)
{
    RBCoxPHTransitionState<MutableArrayHandle<double> > state = args[0];
    if (args[1].isNull() || args[2].isNull()) { return args[0]; }
    MappedColumnVector x;
    try {
        // an exception is raised in the backend if input data contains nulls
        MappedColumnVector xx = args[1].getAs<MappedColumnVector>();
        // x is a const reference, we can only rebind to change its pointer
        x.rebind(xx.memoryHandle(), xx.size());
    } catch (const ArrayWithNullException &e) {
        // independent variable array contains NULL. We skip this row
        return args[0];
    }
    double y = args[2].getAs<double>();
    bool status;
    if (args[3].isNull()) {
        // by default we assume that the data is uncensored => status = TRUE
        status = true;
    } else {
        status = args[3].getAs<bool>();
    }

    MappedColumnVector H = args[6].getAs<MappedColumnVector>();
    double S = args[7].getAs<double>();

    // The following check was added with MADLIB-138.
    if (!dbal::eigen_integration::isfinite(x))
        throw std::domain_error("Design matrix is not finite.");

    if (x.size() > std::numeric_limits<uint16_t>::max())
        throw std::domain_error(
            "Number of independent variables cannot be larger than 65535.");

    if (state.numRows == 0) {
        state.initialize(*this, static_cast<uint16_t>(x.size()),
                         args[4].getAs<MappedColumnVector>().data(),
                         args[5].getAs<MappedMatrix>().data());
    }

    state.numRows++;

    /** In case of a tied time of death or in the first iteration:
        We must only perform the "pre compuations". When the tie is resolved
        we add up all the precomputations once in for all. This is
        an implementation of Breslow's method.
        The time of death for two records are considered "equal" if they
        differ by less than 1.0e-6.
        Also, in case status = 0, the observation must be censored so no
        computations are required
    */

    if (std::abs(y-state.y_previous) > 1.0e-6 && state.numRows != 1) {
        state.M = state.M - (state.tie12 + trans(state.tie12)) * state.A
            + (state.tie13 * trans(state.B) + state.B * trans(state.tie13))
            - (state.tie23 * trans(state.B) + state.B * trans(state.tie23)) * state.A
            + state.tie22 * state.A * state.A + state.tie33 * state.B * trans(state.B);
        state.tie12.setZero();