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

    void DifferentialEvolution::fillInitialPopulation(
                                          std::vector<Candidate> & population,
                                          const Problem& p) const {

        // use initial values provided by the user
        population.front().values = p.currentValue();
        population.front().cost = p.costFunction().value(population.front().values);
        // rest of the initial population is random
        for (Size j = 1; j < population.size(); ++j) {
            for (Size i = 0; i < p.currentValue().size(); ++i) {
                Real l = lowerBound_[i], u = upperBound_[i];
                population[j].values[i] = l + (u-l)*rng_.nextReal();
            }
            population[j].cost = p.costFunction().value(population[j].values);
        }
    }

    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;
    }

    EndCriteria::Type LevenbergMarquardt::minimize(Problem& P,
                                                   const EndCriteria& endCriteria) {
        EndCriteria::Type ecType = EndCriteria::None;
        P.reset();
        Array x_ = P.currentValue();
        currentProblem_ = &P;
        initCostValues_ = P.costFunction().values(x_);
        int m = initCostValues_.size();
        int n = x_.size();
        boost::scoped_array<double> xx(new double[n]);
        std::copy(x_.begin(), x_.end(), xx.get());
        boost::scoped_array<double> fvec(new double[m]);
        boost::scoped_array<double> diag(new double[n]);
        int mode = 1;
        double factor = 1;
        int nprint = 0;
        int info = 0;
        int nfev =0;
        boost::scoped_array<double> fjac(new double[m*n]);
        int ldfjac = m;
        boost::scoped_array<int> ipvt(new int[n]);
        boost::scoped_array<double> qtf(new double[n]);
        boost::scoped_array<double> wa1(new double[n]);
        boost::scoped_array<double> wa2(new double[n]);
        boost::scoped_array<double> wa3(new double[n]);
        boost::scoped_array<double> wa4(new double[m]);
        // requirements; check here to get more detailed error messages.
        QL_REQUIRE(n > 0, "no variables given");
        QL_REQUIRE(m >= n,
                   "less functions (" << m <<
                   ") than available variables (" << n << ")");
        QL_REQUIRE(endCriteria.functionEpsilon() >= 0.0,
                   "negative f tolerance");
        QL_REQUIRE(xtol_ >= 0.0, "negative x tolerance");
        QL_REQUIRE(gtol_ >= 0.0, "negative g tolerance");
        QL_REQUIRE(endCriteria.maxIterations() > 0,
                   "null number of evaluations");

        // call lmdif to minimize the sum of the squares of m functions
        // in n variables by the Levenberg-Marquardt algorithm.
        MINPACK::LmdifCostFunction lmdifCostFunction = 
            boost::bind(&LevenbergMarquardt::fcn, this, _1, _2, _3, _4, _5);
        MINPACK::lmdif(m, n, xx.get(), fvec.get(),
                       static_cast<double>(endCriteria.functionEpsilon()),
                       static_cast<double>(xtol_),
                       static_cast<double>(gtol_),
                       static_cast<int>(endCriteria.maxIterations()),
                       static_cast<double>(epsfcn_),
                       diag.get(), mode, factor,
                       nprint, &info, &nfev, fjac.get(),
                       ldfjac, ipvt.get(), qtf.get(),
                       wa1.get(), wa2.get(), wa3.get(), wa4.get(),
                       lmdifCostFunction);
        info_ = info;
        // check requirements & endCriteria evaluation
        QL_REQUIRE(info != 0, "MINPACK: improper input parameters");
        //QL_REQUIRE(info != 6, "MINPACK: ftol is too small. no further "
        //                               "reduction in the sum of squares "
        //                               "is possible.");
        if (info != 6) ecType = QuantLib::EndCriteria::StationaryFunctionValue;
        //QL_REQUIRE(info != 5, "MINPACK: number of calls to fcn has "
        //                               "reached or exceeded maxfev.");
        endCriteria.checkMaxIterations(nfev, ecType);
        QL_REQUIRE(info != 7, "MINPACK: xtol is too small. no further "
                                       "improvement in the approximate "
                                       "solution x is possible.");
        QL_REQUIRE(info != 8, "MINPACK: gtol is too small. fvec is "
                                       "orthogonal to the columns of the "
                                       "jacobian to machine precision.");
        // set problem
        std::copy(xx.get(), xx.get()+n, x_.begin());
        P.setCurrentValue(x_);
        P.setFunctionValue(P.costFunction().value(x_));
        
        return ecType;
    }

    EndCriteria::Type DifferentialEvolution::minimize(Problem& P,
											const EndCriteria& endCriteria) {

		EndCriteria::Type ecType = EndCriteria::MaxIterations;
	    QL_REQUIRE(P.currentValue().size() == nParam_,
			"Number of parameters mismatch between problem and DE optimizer");
        P.reset();		
		init();

		Real bestCost = QL_MAX_REAL;
		Size bestPop  = 0;
		for (Size p = 0; p < nPop_; ++p) {
			Array tmp(currGen_[p].pop_);
			try {
				currGen_[p].cost_ = P.costFunction().value(tmp);
			} catch (Error&) {
				currGen_[p].cost_ = QL_MAX_REAL;
			}
			if (currGen_[p].cost_ < bestCost) {
				bestPop = p;
				bestCost = currGen_[p].cost_;
			}
		}

		Size lastChange = 0;
		Size lastParamChange = 0;
		for(Size i=0; i<endCriteria.maxIterations(); ++i) {

			Size newBestPop = bestPop;
			Real newBestCost = bestCost;

			for (Size p=0; p<nPop_; ++p) {
				// Find 3 different populations randomly
				Size r1;
				do {
					r1 = static_cast <Size> (uniformRng_.nextInt32() % nPop_);
				}		
				while(r1 == p || r1 == bestPop);

				Size r2;
				do {
					r2 = static_cast <Size> (uniformRng_.nextInt32() % nPop_);
				}
				while ( r2 == p || r2 == bestPop || r2 == r1);

				Size r3;
				do {
					r3 = static_cast <Size> (uniformRng_.nextInt32() % nPop_);
				} while ( r3 == p || r3 == bestPop || r3 == r1 || r3 == r2);

				for(Size j=0; j<nParam_; ++j) {
					nextGen_[p].pop_[j] = currGen_[p].pop_[j];
				}

				Size j = static_cast <Size> (uniformRng_.nextInt32() % nParam_);
				Size L = 0;
				do {
					const double tmp = 
						currGen_[      p].pop_[j] * a0_
					  + currGen_[     r1].pop_[j] * a1_
					  + currGen_[     r2].pop_[j] * a2_
					  + currGen_[     r3].pop_[j] * a3_
					  + currGen_[bestPop].pop_[j] * aBest_;

					nextGen_[p].pop_[j] =
						std::min(maxParams_[j], std::max(minParams_[j], tmp));

					j = (j+1)%nParam_;
					++L;
				} while ((uniformRng_.nextReal() < CR_) && (L < nParam_));

				// Evaluate the new population
				Array tmp(nextGen_[p].pop_);
				try {
					nextGen_[p].cost_ = P.costFunction().value(tmp);
                } catch (Error&) {
					nextGen_[p].cost_ = QL_MAX_REAL;
                }

				// Not better, discard it and keep the old one.
				if (nextGen_[p].cost_ >= currGen_[p].cost_) {
					nextGen_[p] = currGen_[p];
				}
				// Better, keep it.
				else {
					// New best?
					if (nextGen_[p].cost_ < newBestCost) {
						newBestPop = p;
						newBestCost = nextGen_[p].cost_;
					}
				}
			}

			if(std::abs(newBestCost-bestCost) > endCriteria.functionEpsilon()) {
				lastChange = i;
			}
			const Array absDiff = Abs(nextGen_[newBestPop].pop_-currGen_[bestPop].pop_);
			if(*std::max_element(absDiff.begin(), absDiff.end()) > endCriteria.rootEpsilon()) {
				lastParamChange = i;
			}
//.........这里部分代码省略.........