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

    void FireflyAlgorithm::startState(Problem &P, const EndCriteria &endCriteria) {
        N_ = P.currentValue().size();
        x_.reserve(M_);
        xI_.reserve(M_);
        xRW_.reserve(M_);
        values_.reserve(M_);
        uX_ = P.constraint().upperBound(P.currentValue());
        lX_ = P.constraint().lowerBound(P.currentValue());
        Array bounds = uX_ - lX_;

        //Random initialization is done by Sobol sequence
        SobolRsg sobol(N_);

        //Prepare containers
        for (Size i = 0; i < M_; i++) {
            const SobolRsg::sample_type::value_type &sample = sobol.nextSequence().value;
            x_.push_back(Array(N_, 0.0));
            xI_.push_back(Array(N_, 0.0));
            xRW_.push_back(Array(N_, 0.0));
            Array& x = x_.back();
            for (Size j = 0; j < N_; j++) {
                //Assign X=lb+(ub-lb)*random
                x[j] = lX_[j] + bounds[j] * sample[j];
            }
            //Evaluate point
            values_.push_back(std::make_pair(P.value(x), i));
        }

        //init intensity & randomWalk
        intensity_->init(this);
        randomWalk_->init(this);
    }

    EndCriteria::Type DifferentialEvolution::minimize(Problem& p, const EndCriteria& endCriteria) {
        EndCriteria::Type ecType;

        upperBound_ = p.constraint().upperBound(p.currentValue());
        lowerBound_ = p.constraint().lowerBound(p.currentValue());
        currGenSizeWeights_ = Array(configuration().populationMembers,
                                    configuration().stepsizeWeight);
        currGenCrossover_ = Array(configuration().populationMembers,
                                  configuration().crossoverProbability);

        std::vector<Candidate> population(configuration().populationMembers,
                                          Candidate(p.currentValue().size()));
        fillInitialPopulation(population, p);

        std::partial_sort(population.begin(), population.begin() + 1, population.end(),
                          sort_by_cost());
        bestMemberEver_ = population.front();
        Real fxOld = population.front().cost;
        Size iteration = 0, stationaryPointIteration = 0;

        // main loop - calculate consecutive emerging populations
        while (!endCriteria.checkMaxIterations(iteration++, ecType)) {
            calculateNextGeneration(population, p.costFunction());
            std::partial_sort(population.begin(), population.begin() + 1, population.end(),
                              sort_by_cost());
            if (population.front().cost < bestMemberEver_.cost)
                bestMemberEver_ = population.front();
            Real fxNew = population.front().cost;
            if (endCriteria.checkStationaryFunctionValue(fxOld, fxNew, stationaryPointIteration,
                                                         ecType))
                break;
            fxOld = fxNew;
        };
        p.setCurrentValue(bestMemberEver_.values);
        p.setFunctionValue(bestMemberEver_.cost);
        return ecType;
    }

    Real Simplex::extrapolate(Problem& P,
                              Size iHighest,
                              Real &factor) const {

        Array pTry;
        do {
            Size dimensions = values_.size() - 1;
            Real factor1 = (1.0 - factor)/dimensions;
            Real factor2 = factor1 - factor;
            #if defined(QL_ARRAY_EXPRESSIONS)
            pTry = sum_*factor1 - vertices_[iHighest]*factor2;
            #else
            // composite expressions fail to compile with gcc 3.4 on windows
            pTry = sum_*factor1;
            pTry -= vertices_[iHighest]*factor2;
            #endif
            factor *= 0.5;
        } while (!P.constraint().test(pTry) && std::fabs(factor) > QL_EPSILON);
        if (std::fabs(factor) <= QL_EPSILON) {
            return values_[iHighest];
        }
        factor *= 2.0;
        Real vTry = P.value(pTry);
        if (vTry < values_[iHighest]) {
            values_[iHighest] = vTry;
            #if defined(QL_ARRAY_EXPRESSIONS)
            sum_ += pTry - vertices_[iHighest];
            #else
            sum_ += pTry;
            sum_ -= vertices_[iHighest];
            #endif
            vertices_[iHighest] = pTry;
        }
        return vTry;

    }

    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 Simplex::minimize(Problem& P,
                                        const EndCriteria& endCriteria) {
        // set up of the problem
        //Real ftol = endCriteria.functionEpsilon();    // end criteria on f(x) (see Numerical Recipes in C++, p.410)
        Real xtol = endCriteria.rootEpsilon();          // end criteria on x (see GSL v. 1.9, http://www.gnu.org/software/gsl/)
        Size maxStationaryStateIterations_
            = endCriteria.maxStationaryStateIterations();
        EndCriteria::Type ecType = EndCriteria::None;
        P.reset();
        Array x_ = P.currentValue();
        Integer iterationNumber_=0;

        // Initialize vertices of the simplex
        bool end = false;
        Size n = x_.size(), i;
        vertices_ = std::vector<Array>(n+1, x_);
        for (i=0; i<n; i++) {
            Array direction(n, 0.0);
            direction[i] = 1.0;
            P.constraint().update(vertices_[i+1], direction, lambda_);
        }
        // Initialize function values at the vertices of the simplex
        values_ = Array(n+1, 0.0);
        for (i=0; i<=n; i++)
            values_[i] = P.value(vertices_[i]);
        // Loop looking for minimum
        do {
            sum_ = Array(n, 0.0);
            Size i;
            for (i=0; i<=n; i++)
                sum_ += vertices_[i];
            // Determine the best (iLowest), worst (iHighest)
            // and 2nd worst (iNextHighest) vertices
            Size iLowest = 0;
            Size iHighest, iNextHighest;
            if (values_[0]<values_[1]) {
                iHighest = 1;
                iNextHighest = 0;
            } else {
                iHighest = 0;
                iNextHighest = 1;
            }
            for (i=1;i<=n; i++) {
                if (values_[i]>values_[iHighest]) {
                    iNextHighest = iHighest;
                    iHighest = i;
                } else {
                    if ((values_[i]>values_[iNextHighest]) && i!=iHighest)
                        iNextHighest = i;
                }
                if (values_[i]<values_[iLowest])
                    iLowest = i;
            }
            // Now compute accuracy, update iteration number and check end criteria
            //// Numerical Recipes exit strategy on fx (see NR in C++, p.410)
            //Real low = values_[iLowest];
            //Real high = values_[iHighest];
            //Real rtol = 2.0*std::fabs(high - low)/
            //    (std::fabs(high) + std::fabs(low) + QL_EPSILON);
            //++iterationNumber_;
            //if (rtol < ftol ||
            //    endCriteria.checkMaxIterations(iterationNumber_, ecType)) {
            // GSL exit strategy on x (see GSL v. 1.9, http://www.gnu.org/software/gsl
            Real simplexSize = computeSimplexSize(vertices_);
            ++iterationNumber_;
            if (simplexSize < xtol ||
                endCriteria.checkMaxIterations(iterationNumber_, ecType)) {
                endCriteria.checkStationaryPoint(0.0, 0.0,
                maxStationaryStateIterations_, ecType);   // PC this is probably not meant like this ? Use separate counter ?
                endCriteria.checkMaxIterations(iterationNumber_, ecType);
                x_ = vertices_[iLowest];
                Real low = values_[iLowest];
                P.setFunctionValue(low);
                P.setCurrentValue(x_);
                return ecType;
            }
            // If end criteria is not met, continue
            Real factor = -1.0;
            Real vTry = extrapolate(P, iHighest, factor);
            if ((vTry <= values_[iLowest]) && (factor == -1.0)) {
                factor = 2.0;
                extrapolate(P, iHighest, factor);
            } else if (std::fabs(factor) > QL_EPSILON) {
                if (vTry >= values_[iNextHighest]) {
                    Real vSave = values_[iHighest];
                    factor = 0.5;
                    vTry = extrapolate(P, iHighest, factor);
                    if (vTry >= vSave && std::fabs(factor) > QL_EPSILON) {
                        for (Size i=0; i<=n; i++) {
                            if (i!=iLowest) {
                                #if defined(QL_ARRAY_EXPRESSIONS)
                                vertices_[i] =
                                    0.5*(vertices_[i] + vertices_[iLowest]);
                                #else
                                vertices_[i] += vertices_[iLowest];
                                vertices_[i] *= 0.5;
                                #endif
                                values_[i] = P.value(vertices_[i]);
                            }
                        }
//.........这里部分代码省略.........

    EndCriteria::Type FireflyAlgorithm::minimize(Problem &P, const EndCriteria &endCriteria) {
        QL_REQUIRE(!P.constraint().empty(), "Firefly Algorithm is a constrained optimizer");
        EndCriteria::Type ecType = EndCriteria::None;
        P.reset();
        Size iteration = 0;
        Size iterationStat = 0;
        Size maxIteration = endCriteria.maxIterations();
        Size maxIStationary = endCriteria.maxStationaryStateIterations();
        
        startState(P, endCriteria);

        bool isFA = Mfa_ > 0 ? true : false;
        //Variables for DE
        Array z(N_, 0.0);
        Size indexBest, indexR1, indexR2;

        //Set best value & position
        Real bestValue = values_[0].first;
        Size bestPosition = 0;
        for (Size i = 1; i < M_; i++) {
            if (values_[i].first < bestValue) {
                bestPosition = i;
                bestValue = values_[i].first;
            }
        }
        Array bestX = x_[bestPosition];

        //Run optimization
        do {
            iteration++;
            iterationStat++;
            //Check if stopping criteria is met
            if (iteration > maxIteration || iterationStat > maxIStationary)
                break;

            //Divide into two subpopulations
            //First sort values
            std::sort(values_.begin(), values_.end());

            //Differential evolution
            if(Mfa_ < M_){
                Size indexBest = values_[0].second;
                Array& xBest = x_[indexBest];
                for (Size i = Mfa_; i < M_; i++) { 
                    if (!isFA) {
                        //Pure DE requires random index
                        indexBest = drawIndex_();
                        xBest = x_[indexBest];
                    }
                    do { 
                        indexR1 = drawIndex_(); 
                    } while(indexR1 == indexBest);
                    do { 
                        indexR2 = drawIndex_(); 
                    } while(indexR2 == indexBest || indexR2 == indexR1);
                    
                    Size index = values_[i].second;
                    Array& x   = x_[index];
                    Array& xR1 = x_[indexR1];
                    Array& xR2 = x_[indexR2];
                    for (Size j = 0; j < N_; j++) {
                        if (rng_.nextReal() <= crossover_) {
                            //Change x[j] according to crossover
                            z[j] = xBest[j] + mutation_*(xR1[j] - xR2[j]);
                        } else {
                            z[j] = x[j];
                        }
                    }
                    Real val = P.value(z);
                    if (val < values_[index].first) {
                        //Accept new point
                        x = z;
                        values_[index].first = val;
                        //mark best
                        if (val < bestValue) {
                            bestValue = val;
                            bestX = x;
                            iterationStat = 0;
                        }
                    }
                }
            }
                
            //Firefly algorithm
            if(isFA){
                //According to the intensity, determine best global position
                intensity_->findBrightest();

                //Prepare random walk
                randomWalk_->walk();

                //Loop over particles
                for (Size i = 0; i < Mfa_; i++) {
                    Size index = values_[i].second;
                    Array& x   = x_[index];
                    Array& xI  = xI_[index];
                    Array& xRW = xRW_[index];

                    //Loop over dimensions
                    for (Size j = 0; j < N_; j++) {
//.........这里部分代码省略.........