static HResults computeH()

in include/Algo-Direct-Common.h [215:281]
55 lines of code
8 McCabe index (conditional complexity)

    static HResults<T> computeH(const T *px, uint32 nx)
    {
        myassert((nx > Gap), "Array X too small");
        myassert(((Gap == 1) || (Gap == 2)), "Only tested for these values of Gap");

        const T x0 = px[0];
        const T xN = px[nx-1];

        const T range = xN - x0;
        myassert((range < std::numeric_limits<T>::max()), "range too large");

        // check that D_i are strictly increasing and compute minimum value D_{i+Offset}-D_i
        T deltaDMin = range;
        for (uint32 i = Gap; i < nx; ++i) {
            T Dnew = px[i] - x0;
            T Dold = px[i - Gap] - x0;
            myassert((Dnew > Dold),
                "Problem unfeasible: D_i sequence not strictly increasing"
                << " X[" << 0 << "]=" << x0
                << " X[" << i - Gap << "]=" << px[i - Gap]
                << " X[" << i << "]=" << px[i]
                << "\n"
            );
            T deltaD = Dnew - Dold;
            if (deltaD < deltaDMin)
                deltaDMin = deltaD;
        }

        // initial guess for H
        const T H0 = T(1.0) / deltaDMin;
        T H = H0;

        T cst0 = fun_t::cst0(H, x0);
        fun_t::checkH(H, cst0, xN);

        // adjust H by trial and error until succeed
        size_t nInc = 0;
        bool modified = false;
        size_t npasses = 0;
        T step = growStep(H);
        uint32 seg_already_checked_from = nx;
        do {
            myassert((npasses++ < 2), "verification failed\n");
            // if there has been an increase, then check only up to that point
            uint32 last_seg_to_be_checked = seg_already_checked_from - 1;
            modified = false;
            uint32 inew = 0;
            for (uint32 i = Gap; i <= last_seg_to_be_checked; ++i) {
                uint32 iold = fun_t::f(H, cst0, px[i-Gap]);
                uint32 inew = fun_t::f(H, cst0, px[i]);
                while (inew == iold) {
                    seg_already_checked_from = i;
                    last_seg_to_be_checked = nx-1;  // everything needs to be checked
                    modified = true;
                    H = H + step;
                    step *= 2;
                    // recalculate all constants and indices
                    cst0 = fun_t::cst0(H, x0);
                    fun_t::checkH(H, cst0, xN);
                    iold = fun_t::f(H, cst0, px[i - Gap]);
                    inew = fun_t::f(H, cst0, px[i]);
                }
            }
        } while (SAFETY_MULTI_PASS && modified);

        return HResults<T>(H, (((double)H) / H0) - 1.0, nInc);
    }