partial_filtered_pricing.h

来自「CGAL is a collaborative effort of severa」· C头文件 代码 · 共 567 行 · 第 1/2 页

            if ( mu >= et_0) {
                vout() << "entering variable defeated by exact check\n";
                min_k = -1; min_j = -1; nt_min_mu = nt_0;
            }
        }
    
        if ( min_k < 0) {
    
        // --------------------------------------------------------------------
        vout() << "no entering variable found so far, test remaining variables"
        // --------------------------------------------------------------------
               << std::endl;
    
            // loop over all remaining non-basic variables
            for ( k = s; k < (int)N.size(); ++k) {
    
                j = N[ k];
    
                // compute mu_j
                if ( is_phase_I) {      // phase I
                    if ( j < n) {          // original variable
                        nt_mu = std::inner_product(
                            lambda.begin(), lambda.end(),
                            solve.a_begin()[ j],
                            nt_d * solve.c_auxiliary_begin()[ j]);
                    } else {               // artificial variable
                        nt_mu = std::inner_product(
                            lambda.begin(), lambda.end(),
                            solve.a_artificial_begin()[ j-n],
                            nt_d * solve.c_auxiliary_begin()[ j]);
                    }
                } else {                // phase II
                    nt_mu = std::inner_product(
                        lambda.begin(), lambda.end(),
                        solve.a_begin()[ j],
                        nt_d * solve.c_begin()[ j]);
                    // is QP?
                    if ( ! CGAL::check_tag( Is_lp())) {
                        nt_mu += nt_2 * std::inner_product(
                            x_B.begin(), x_B.end(),
                            D_Bj_iterator( solve.basic_variables_index_begin(),
                                           Access_D_Bj( solve.d_begin()[ j])),
                            nt_0);
                    }
                }
    
                CGAL_optimisation_debug {
                    vout() << "nt_mu_" << j << ": " << nt_mu << std::endl;
                }
    
                // improving variable?
                if ( nt_mu < nt_0) {
                    std::swap( N[ k], N[ s]);
    
                    // new minimum?
                    if ( ( nt_mu < nt_min_mu) ||
                         ( ( min_j >= n) && ( j < n) &&
                           ( nt_mu == nt_min_mu))) {
                        min_k  = s;
                        min_j  = j;
                        nt_min_mu = nt_mu;
                    }
    
                    ++s;
                }
            }
    
            // exact check of entering variable
            if ( min_k >= 0) {
                j = N[ min_k];
                if ( is_phase_I) {      // phase I
                    if ( j < n) {          // original variable
                        mu = std::inner_product(
                            solve.lambda_numerator_begin(),
                            solve.lambda_numerator_end(),
                            solve.a_begin()[ j],
                            d * solve.c_auxiliary_begin()[ j]);
                    } else {               // artificial variable
                        mu = std::inner_product(
                            solve.lambda_numerator_begin(),
                            solve.lambda_numerator_end(),
                            solve.a_artificial_begin()[ j-n],
                            d * solve.c_auxiliary_begin()[ j]);
                    }
                } else {                // phase II
                    mu = std::inner_product(
                        solve.lambda_numerator_begin(),
                        solve.lambda_numerator_end(),
                        solve.a_begin()[ j],
                        d * solve.c_begin()[ j]);
                    // is QP?
                    if ( ! CGAL::check_tag( Is_lp())) {
                        mu += et_2 * std::inner_product(
                            solve.basic_variables_numerator_begin(),
                            solve.basic_variables_numerator_end(),
                            D_Bj_iterator( solve.basic_variables_index_begin(),
                                           Access_D_Bj( solve.d_begin()[ j])),
                            et_0);
                    }
                }
                if ( mu >= et_0) {
                    vout() << "entering variable defeated by exact check\n";
                    min_k = -1; nt_min_mu = nt_0;
                }
            }
        }
        if ( min_k < 0) {
    
            // ----------------------------------------------------------------
            vout()
            << "no entering variable found so far, revert to exact arithmetic"
            // ----------------------------------------------------------------
            << std::endl;
    
            // compute first error bound
            k = m+b+1;
            NT q = ldexp( 1.015625*k*(k+1), -53);
            NT max_1 = nt_d * row_max_c;
            NT max_2 = nt_d;
            NT z;
            for ( i = 0; i < m; ++i) {
                z = CGAL_NTS abs( lambda[ i]) * row_max_A[ i];
                if ( z > max_1) max_1 = z;
                z = CGAL_NTS abs( lambda[ i]);
                if ( z > max_2) max_2 = z;
            }
            if ( ! CGAL::check_tag( Is_lp())) {
                typename std::iterator_traits<D_iterator>::value_type  row_D;
                for ( i = 0; i < b; ++i) {
                    k = solve.basic_variables_index_begin()[ i];
                    row_D = solve.d_begin()[ k];
                    if ( ! row_valid[ k]) {
                        NT  max = nt_0;
                        for ( j = 0; j < n; ++j) {
                            z = CGAL_NTS abs( row_D[ j]);
                            if ( z >     max    )     max     = z;
                            if ( z > col_max[ j]) col_max[ j] = z;
                        }
                        row_max_D[ k] = max;
                    }
                    z = CGAL_NTS abs( x_B[ i]) * row_max_D[ k];
                    if ( z > max_1) max_1 = z;
                    z = CGAL_NTS abs( x_B[ i]);
                    if ( z > max_2) max_2 = z;
                }
            }
            NT bound_1 = max_1 * q, max_u_q = max_2 * q, bound_2;
            CGAL_optimisation_debug {
                vout() << "[ first bound: " << bound_1 << " ]" << std::endl;
            }
    
            // loop again over all non-basic variables to verify optimality
            k = 0;
            while ( k < (int)N.size() && min_k < 0) {
    
                j = N[ k];
    
                // compute mu_j (inexact)
                if ( is_phase_I) {      // phase I
                    if ( j < n) {          // original variable
                        nt_mu = std::inner_product(
                            lambda.begin(), lambda.end(),
                            solve.a_begin()[ j],
                            nt_d * solve.c_auxiliary_begin()[ j]);
                    } else {               // artificial variable
                        nt_mu = std::inner_product(
                            lambda.begin(), lambda.end(),
                            solve.a_artificial_begin()[ j-n],
                            nt_d * solve.c_auxiliary_begin()[ j]);
                    }
                } else {                // phase II
                    nt_mu = std::inner_product(
                        lambda.begin(), lambda.end(),
                        solve.a_begin()[ j],
                        nt_d * solve.c_begin()[ j]);
                    // is QP?
                    if ( ! CGAL::check_tag( Is_lp())) {
                        nt_mu += nt_2 * std::inner_product(
                            x_B.begin(), x_B.end(),
                            D_Bj_iterator( solve.basic_variables_index_begin(),
                                           Access_D_Bj( solve.d_begin()[ j])),
                            nt_0);
                    }
                }
    
                CGAL_optimisation_debug {
                    vout() << "nt_mu_" << j << ": " << nt_mu;
                }
    
                // check against first bound
                if ( nt_mu >= bound_1) {
                    CGAL_optimisation_debug {
                        vout() << " [ certified by first bound ]" << std::endl;
                    }
                } else {
                    // compute second bound
                    bound_2 = col_max[ j] * max_u_q;
                    if ( nt_mu >= bound_2) {
                        CGAL_optimisation_debug {
                            vout() << " [ certified by second bound: "
                                   << bound_2 << " ]" << std::endl;
                        }
                    } else {
    
                        // compute mu_j (exact)
                        if ( is_phase_I) {      // phase I
                            if ( j < n) {          // original variable
                                mu = std::inner_product(
                                    solve.lambda_numerator_begin(),
                                    solve.lambda_numerator_end(),
                                    solve.a_begin()[ j],
                                    d * solve.c_auxiliary_begin()[ j]);
                            } else {               // artificial variable
                                mu = std::inner_product(
                                    solve.lambda_numerator_begin(),
                                    solve.lambda_numerator_end(),
                                    solve.a_artificial_begin()[ j-n],
                                    d * solve.c_auxiliary_begin()[ j]);
                            }
                        } else {                // phase II
                            mu = std::inner_product(
                                solve.lambda_numerator_begin(),
                                solve.lambda_numerator_end(),
                                solve.a_begin()[ j],
                                d * solve.c_begin()[ j]);
                            // is QP?
                            if ( ! CGAL::check_tag( Is_lp())) {
                                mu += et_2 * std::inner_product(
                                    solve.basic_variables_numerator_begin(),
                                    solve.basic_variables_numerator_end(),
                            D_Bj_iterator( solve.basic_variables_index_begin(),
                                           Access_D_Bj( solve.d_begin()[ j])),
                                    et_0);
                            }
                        }
    
                        CGAL_optimisation_debug {
                            vout() << " [ exact computation needed: "
                                   << mu << " ]" << std::endl;
                        }
    
                        if ( mu < et_0) min_k = k;
                    }
                }
                ++k;
            }
        }
        vout() << std::endl;
    
        // return index of entering variable
        if ( min_k >= 0) {
            j = N[ min_k];
            --s;
            N[ min_k] = N[ s];
            N[ s] = N.back();
            N.pop_back();
            return j;
        }
        return -1;
    }
    
    void  leaving_basis( int i)
    {
        if ( s == (int)N.size()) {
            N.push_back( i);
        } else {
            N.push_back( N[ s]);
            N[ s] = i;
        }
        ++s;
    }
    
    
};
  

CGAL_END_NAMESPACE
                  

#endif // CGAL_PARTIAL_FILTERED_PRICING_H

// ===== EOF ==================================================================