C++ Problem::functionValue方法代码示例

    EndCriteria::Type HybridSimulatedAnnealing<Sampler, Probability, Temperature, Reannealing>::minimize(Problem &P, const EndCriteria &endCriteria) {
        EndCriteria::Type ecType = EndCriteria::None;
        P.reset();
        reannealing_.setProblem(P);
        Array x = P.currentValue();
        Size n = x.size();
        Size k = 1;
        Size kStationary = 1;
        Size kReAnneal = 1;
        Size kReset = 1;
        Size maxK = endCriteria.maxIterations();
        Size maxKStationary = endCriteria.maxStationaryStateIterations();
        bool temperatureBreached = false;
        Array currentTemperature(n, startTemperature_);
        Array annealStep(n, 1.0);
        Array bestPoint(x);
        Array currentPoint(x);
        Array startingPoint(x);
        Array newPoint(x);
        Real bestValue = P.value(bestPoint);
        Real currentValue = bestValue;
        Real startingValue = bestValue; //to reset to starting point if desired
        while (k <= maxK && kStationary <= maxKStationary && !temperatureBreached)
        {
            //Draw a new sample point
            sampler_(newPoint, currentPoint, currentTemperature);

            //Evaluate new point
            Real newValue = P.value(newPoint);

            //Determine if new point is accepted
            if (probability_(currentValue, newValue, currentTemperature)) {
                if (optimizeScheme_ == EveryNewPoint) {
                    P.setCurrentValue(newPoint);
                    P.setFunctionValue(newValue);
                    localOptimizer_->minimize(P, endCriteria);
                    newPoint = P.currentValue();
                    newValue = P.functionValue();
                }
                currentPoint = newPoint;
                currentValue = newValue;
            }

            //Check if we have a new best point
            if (newValue < bestValue) {
                if (optimizeScheme_ == EveryBestPoint) {
                    P.setCurrentValue(newPoint);
                    P.setFunctionValue(newValue);
                    localOptimizer_->minimize(P, endCriteria);
                    newPoint = P.currentValue();
                    newValue = P.functionValue();
                }
                kStationary = 0;
                bestValue = newValue;
                bestPoint = newPoint;
            }

            //Increase steps
            k++;
            kStationary++;
            for (Size i = 0; i < annealStep.size(); i++)
                annealStep[i]++;

            //Reanneal if necessary
            if (kReAnneal == reAnnealSteps_) {
                kReAnneal = 0;
                reannealing_(annealStep, currentPoint, currentValue, currentTemperature);
            }
            kReAnneal++;

            //Reset if necessary
            if (kReset == resetSteps_) {
                kReset = 0;
                switch (resetScheme_) {
                case NoResetScheme:
                    break;
                case ResetToBestPoint:
                    currentPoint = startingPoint;
                    currentValue = startingValue;
                    break;
                case ResetToOrigin:
                    currentPoint = bestPoint;
                    currentValue = bestValue;
                    break;
                }
            }
            kReset++;

            //Update the current temperature according to current step
            temperature_(currentTemperature, currentTemperature, annealStep);

            //Check if temperature condition is breached
            for (Size i = 0; i < n; i++)
                temperatureBreached = temperatureBreached && currentTemperature[i] < endTemperature_;
        }
        //Change end criteria type if appropriate
        if (k > maxK)
            ecType = EndCriteria::MaxIterations;
        else if (kStationary > maxKStationary)
            ecType = EndCriteria::StationaryPoint;
//.........这里部分代码省略.........

    Real GoldsteinLineSearch::operator()(Problem& P, 
                                         EndCriteria::Type& ecType, 
                                         const EndCriteria& endCriteria, 
                                         const Real t_ini) 
    {
        Constraint& constraint = P.constraint();
        succeed_=true;
        bool maxIter = false;
        Real /*qtold,*/ t = t_ini; // see below, this is never used ?
        Size loopNumber = 0;

        Real q0 = P.functionValue();
        Real qp0 = P.gradientNormValue();
        
        Real tl = 0.0;
        Real tr = 0.0;

        qt_ = q0;
        qpt_ = (gradient_.empty()) ? qp0 : -DotProduct(gradient_,searchDirection_);

        // Initialize gradient
        gradient_ = Array(P.currentValue().size());
        // Compute new point
        xtd_ = P.currentValue();
        t = update(xtd_, searchDirection_, t, constraint);
        // Compute function value at the new point
        qt_ = P.value (xtd_);

        while ((qt_ - q0) < -beta_*t*qpt_ || (qt_ - q0) > -alpha_*t*qpt_) {
            if ((qt_ - q0) > -alpha_*t*qpt_)
                tr = t;
            else
                tl = t;
            ++loopNumber;
            
            // calculate the new step
            if (close_enough(tr, 0.0))
                t *= extrapolation_;
            else
                t = (tl + tr) / 2.0;

            // Store old value of the function
            // qtold = qt_; // this is never used ?
            // New point value
            xtd_ = P.currentValue();
            t = update(xtd_, searchDirection_, t, constraint);

            // Compute function value at the new point
            qt_ = P.value (xtd_);
            P.gradient (gradient_, xtd_);
            // and it squared norm
            maxIter = endCriteria.checkMaxIterations(loopNumber, ecType);

            if (maxIter)
                break;
        }

        if (maxIter)
            succeed_ = false;

        // Compute new gradient
        P.gradient(gradient_, xtd_);
        // and it squared norm
        qpt_ = DotProduct(gradient_, gradient_);

        // Return new step value
        return t;
    }

    EndCriteria::Type
    LineSearchBasedMethod::minimize(Problem& P,
                                    const EndCriteria& endCriteria) {
        // Initializations
        Real ftol = endCriteria.functionEpsilon();
        Size maxStationaryStateIterations_
            = endCriteria.maxStationaryStateIterations();
        EndCriteria::Type ecType = EndCriteria::None;   // reset end criteria
        P.reset();                                      // reset problem
        Array x_ = P.currentValue();              // store the starting point
        Size iterationNumber_ = 0;
        // dimension line search
        lineSearch_->searchDirection() = Array(x_.size());
        bool done = false;

        // function and squared norm of gradient values;
        Real fnew, fold, gold2;
        Real fdiff;
        // classical initial value for line-search step
        Real t = 1.0;
        // Set gradient g at the size of the optimization problem
        // search direction
        Size sz = lineSearch_->searchDirection().size();
        Array prevGradient(sz), d(sz), sddiff(sz), direction(sz);
        // Initialize cost function, gradient prevGradient and search direction
        P.setFunctionValue(P.valueAndGradient(prevGradient, x_));
        P.setGradientNormValue(DotProduct(prevGradient, prevGradient));
        lineSearch_->searchDirection() = -prevGradient;

        bool first_time = true;
        // Loop over iterations
        do {
            // Linesearch
            if (!first_time)
                prevGradient = lineSearch_->lastGradient();
            t = (*lineSearch_)(P, ecType, endCriteria, t);
            // don't throw: it can fail just because maxIterations exceeded
            //QL_REQUIRE(lineSearch_->succeed(), "line-search failed!");
            if (lineSearch_->succeed())
            {
                // Updates

                // New point
                x_ = lineSearch_->lastX();
                // New function value
                fold = P.functionValue();
                P.setFunctionValue(lineSearch_->lastFunctionValue());
                // New gradient and search direction vectors

                // orthogonalization coef
                gold2 = P.gradientNormValue();
                P.setGradientNormValue(lineSearch_->lastGradientNorm2());

                // conjugate gradient search direction
                direction = getUpdatedDirection(P, gold2, prevGradient);

                sddiff = direction - lineSearch_->searchDirection();
                lineSearch_->searchDirection() = direction;
                // Now compute accuracy and check end criteria
                // Numerical Recipes exit strategy on fx (see NR in C++, p.423)
                fnew = P.functionValue();
                fdiff = 2.0*std::fabs(fnew-fold) /
                        (std::fabs(fnew) + std::fabs(fold) + QL_EPSILON);
                if (fdiff < ftol ||
                    endCriteria.checkMaxIterations(iterationNumber_, ecType)) {
                    endCriteria.checkStationaryFunctionValue(0.0, 0.0,
                        maxStationaryStateIterations_, ecType);
                    endCriteria.checkMaxIterations(iterationNumber_, ecType);
                    return ecType;
                }
                P.setCurrentValue(x_);      // update problem current value
                ++iterationNumber_;         // Increase iteration number
                first_time = false;
            } else {
                done = true;
            }
        } while (!done);
        P.setCurrentValue(x_);
        return ecType;
    }