C++ Problem类代码示例

TEST_F(AnalysisDriverFixture,RuntimeBehavior_StopCustomAnalysis) {
  // Tests for stopping time < 20s.

  // RETRIEVE PROBLEM
  Problem problem = retrieveProblem("UserScriptContinuous",true,false);

  // DEFINE SEED
  Model model = model::exampleModel();
  openstudio::path p = toPath("./example.osm");
  model.save(p,true);
  FileReference seedModel(p);

  // CREATE ANALYSIS
  Analysis analysis("Stop Custom Analysis",
                    problem,
                    seedModel);
  // generate 100 random points
  boost::mt19937 mt;
  typedef boost::uniform_real<> dist_type;
  typedef boost::variate_generator<boost::mt19937&, dist_type > gen_type;
  InputVariableVector variables = problem.variables();
  ContinuousVariable cvar = variables[0].cast<ContinuousVariable>();
  gen_type generator0(mt,dist_type(cvar.minimum().get(),cvar.maximum().get()));
  cvar = variables[1].cast<ContinuousVariable>();
  gen_type generator1(mt,dist_type(cvar.minimum().get(),cvar.maximum().get()));
  cvar = variables[2].cast<ContinuousVariable>();
  gen_type generator2(mt,dist_type(cvar.minimum().get(),cvar.maximum().get()));

  for (int i = 0, n = 100; i < n; ++i) {
    std::vector<QVariant> values;
    double value = generator0();
    values.push_back(value);
    value = generator1();
    values.push_back(value);
    value = generator2();
    values.push_back(value);
    OptionalDataPoint dataPoint = problem.createDataPoint(values);
    ASSERT_TRUE(dataPoint);
    ASSERT_TRUE(analysis.addDataPoint(*dataPoint));
  }

  // RUN ANALYSIS
  ProjectDatabase database = getCleanDatabase("StopCustomAnalysis");
  AnalysisDriver analysisDriver(database);
  AnalysisRunOptions runOptions = standardRunOptions(analysisDriver.database().path().parent_path());
  runOptions.setQueueSize(2);
  StopWatcher watcher(analysisDriver);
  watcher.watch(analysis.uuid());
  CurrentAnalysis currentAnalysis = analysisDriver.run(analysis,runOptions);
  EXPECT_EQ(2,currentAnalysis.numQueuedJobs());
  EXPECT_EQ(0,currentAnalysis.numQueuedDakotaJobs());
  EXPECT_EQ(100,currentAnalysis.totalNumJobsInOSIteration());
  EXPECT_EQ(0,currentAnalysis.numCompletedJobsInOSIteration());
  analysisDriver.waitForFinished();
  EXPECT_FALSE(analysisDriver.isRunning());
  EXPECT_GE(watcher.nComplete(),watcher.stopNum());
  EXPECT_LE(watcher.stoppingTime(),openstudio::Time(0,0,0,20));

  // check conditions afterward
  RunManager runManager = analysisDriver.database().runManager();
  EXPECT_FALSE(runManager.workPending());
  BOOST_FOREACH(const Job& job,runManager.getJobs()) {
    EXPECT_FALSE(job.running());
    EXPECT_FALSE(job.treeRunning());
  }
  EXPECT_TRUE(currentAnalysis.numCompletedJobsInOSIteration() > 0);
  EXPECT_TRUE(currentAnalysis.analysis().dataPointsToQueue().size() > 0u);
  EXPECT_TRUE(currentAnalysis.analysis().dataPointsToQueue().size() < 100u);
  EXPECT_EQ(0u,analysisDriver.currentAnalyses().size());
}

  /*
   * BIG MACHINE
   */
  std::list<LabelPosition*>* Pal::labeller( int nbLayers, char **layersName , double *layersFactor, double scale, double bbox[4], PalStat **stats, bool displayAll )
  {
#ifdef _DEBUG_
    std::cout << "LABELLER (selection)" << std::endl;
#endif

    Problem *prob;

    SearchMethod old_searchMethod = searchMethod;

    if ( displayAll )
    {
      setSearch( POPMUSIC_TABU );
    }

#ifdef _VERBOSE_
    clock_t start = clock();
    double create_time;
    std::cout << std::endl << "bbox: " << bbox[0] << " " << bbox[1] << " " << bbox[2] << " " << bbox[3] << std::endl;
#endif

#ifdef _EXPORT_MAP_
    // TODO this is not secure
    std::ofstream svgmap( "pal-map.svg" );

    svgmap << "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>" << std::endl
    << "<svg" << std::endl
    << "xmlns:dc=\"http://purl.org/dc/elements/1.1/\"" << std::endl
    << "xmlns:cc=\"http://creativecommons.org/ns#\"" << std::endl
    << "xmlns:rdf=\"http://www.w3.org/1999/02/22-rdf-syntax-ns#\"" << std::endl
    << "xmlns:svg=\"http://www.w3.org/2000/svg\"" << std::endl
    << "xmlns=\"http://www.w3.org/2000/svg\"" << std::endl
    << "xmlns:sodipodi=\"http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd\"" << std::endl
    << "xmlns:inkscape=\"http://www.inkscape.org/namespaces/inkscape\"" << std::endl
    << "width=\"" << convert2pt( bbox[2] - bbox[0], scale, dpi )  << "\"" << std::endl
    << "height=\"" << convert2pt( bbox[3] - bbox[1], scale, dpi )  << "\">" << std::endl; // TODO xmax ymax
#endif

    QTime t;
    t.start();

    // First, extract the problem
    // TODO which is the minimum scale ? (> 0, >= 0, >= 1, >1 )
    if ( scale < 1 || ( prob = extract( nbLayers, layersName, layersFactor, bbox[0], bbox[1], bbox[2], bbox[3], scale,
#ifdef _EXPORT_MAP_
                                        & svgmap
#else
                                        NULL
#endif
                                      ) ) == NULL )
    {

#ifdef _VERBOSE_
      if ( scale < 1 )
        std::cout << "Scale is 1:" << scale << std::endl;
      else
        std::cout << "empty problem... finishing" << std::endl;
#endif

#ifdef _EXPORT_MAP_
      svgmap << "</svg>" << std::endl;
      svgmap.close();
#endif

      // nothing to be done => return an empty result set
      if ( stats )
        ( *stats ) = new PalStat();
      return new std::list<LabelPosition*>();
    }

    std::cout << "PAL EXTRACT: " << t.elapsed() / 1000.0 << " s" << std::endl;
    t.restart();

    // reduce number of candidates
    // (remove candidates which surely won't be used)
    prob->reduce();

#ifdef _VERBOSE_
    std::cerr << prob->nblp << "\t"
              << prob->nbOverlap;
#endif


    prob->displayAll = displayAll;

#ifdef _VERBOSE_
    create_time = double( clock() - start ) / double( CLOCKS_PER_SEC );

    std::cout << std::endl << "Problem : " << prob->nblp << " candidates for " << prob->nbft << " features makes " << prob->nbOverlap << " overlaps" << std::endl;
    std::cout << std::endl << "Times:"  << std::endl << "    to create problem:  " << create_time << std::endl;
#endif

    // search a solution
    if ( searchMethod == FALP )
      prob->init_sol_falp();
    else if ( searchMethod == CHAIN )
      prob->chain_search();
//.........这里部分代码省略.........

bool solve_translations_problem_l2_chordal
(
  const int* edges,
  const double* poses,
  const double* weights,
  int num_edges,
  double loss_width,
  double* X,
  double function_tolerance,
  double parameter_tolerance,
  int max_iterations
)
{
  // seed the random number generator
  std::srand( std::time( nullptr ) );

  // re index the edges to be a sequential set
  std::vector<int> reindexed_edges(edges, edges+2*num_edges);
  std::vector<int> reindexed_lookup;
  reindex_problem(&reindexed_edges[0], num_edges, reindexed_lookup);
  const int num_nodes = reindexed_lookup.size();

  // Init with a random guess solution
  std::vector<double> x(3*num_nodes);
  for (int i=0; i<3*num_nodes; ++i)
    x[i] = (double)rand() / RAND_MAX;

  // add the parameter blocks (a 3-vector for each node)
  Problem problem;
  for (int i=0; i<num_nodes; ++i)
    problem.AddParameterBlock(&x[3*i], 3);

  // set the residual function (chordal distance for each edge)
  for (int i=0; i<num_edges; ++i) {
    CostFunction* cost_function =
      new AutoDiffCostFunction<ChordFunctor, 3, 3, 3>(
      new ChordFunctor(poses+3*i, weights[i]));

    if (loss_width == 0.0) {
      // No robust loss function
      problem.AddResidualBlock(cost_function, nullptr, &x[3*reindexed_edges[2*i+0]], &x[3*reindexed_edges[2*i+1]]);
    } else {
      problem.AddResidualBlock(cost_function, new ceres::HuberLoss(loss_width), &x[3*reindexed_edges[2*i+0]], &x[3*reindexed_edges[2*i+1]]);
    }
  }

  // Fix first camera in {0,0,0}: fix the translation ambiguity
  x[0] = x[1] = x[2] = 0.0;
  problem.SetParameterBlockConstant(&x[0]);

  // solve
  Solver::Options options;
#ifdef OPENMVG_USE_OPENMP
  options.num_threads = omp_get_max_threads();
#if CERES_VERSION_MAJOR < 2
  options.num_linear_solver_threads = omp_get_max_threads();
#endif
#endif // OPENMVG_USE_OPENMP
  options.minimizer_progress_to_stdout = false;
  options.logging_type = ceres::SILENT;
  options.max_num_iterations = max_iterations;
  options.function_tolerance = function_tolerance;
  options.parameter_tolerance = parameter_tolerance;

  // Since the problem is sparse, use a sparse solver iff available
  if (ceres::IsSparseLinearAlgebraLibraryTypeAvailable(ceres::SUITE_SPARSE))
  {
    options.sparse_linear_algebra_library_type = ceres::SUITE_SPARSE;
    options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
  }
  else if (ceres::IsSparseLinearAlgebraLibraryTypeAvailable(ceres::CX_SPARSE))
  {
    options.sparse_linear_algebra_library_type = ceres::CX_SPARSE;
    options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
  }
  else if (ceres::IsSparseLinearAlgebraLibraryTypeAvailable(ceres::EIGEN_SPARSE))
  {
    options.sparse_linear_algebra_library_type = ceres::EIGEN_SPARSE;
    options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
  }
  else
  {
    options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY;
  }

  Solver::Summary summary;
  Solve(options, &problem, &summary);

  std::cout << summary.FullReport() << "\n";

  if (summary.IsSolutionUsable())
  {
    // undo the re indexing
    for (int i=0; i<num_nodes; ++i) {
      const int j = reindexed_lookup[i];
      X[3*j+0] = x[3*i+0];
      X[3*j+1] = x[3*i+1];
      X[3*j+2] = x[3*i+2];
    }
  }
//.........这里部分代码省略.........

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

    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;
			}
//.........这里部分代码省略.........

    void TDefense::problemAction()
{
        if (myProblems.size() > 0) {
            for (int i = 0; i < myProblems.size(); i++)
            {
                Problem *curProblem = myProblems[i];
                curProblem->setVX(0);
                curProblem->setVY(1);

                double x = curProblem->getX();
                double y = curProblem->getY();
                double vx = curProblem->getVX();
                double vy = curProblem->getVY();
                int height = curProblem->getImage().height();
                int width = curProblem->getImage().width();

                if ((y <= 0 && vy < 0) || (x <= 0 && vx < 0 || (gpHeight <= (y + height)) && y > 0 && vy > 0) || ((gpWidth == (x + width)) && vx > 0 && vx > 0))
                {
                    if (!curProblem->isDropped)
                    {
                        numEnemiesDropped++;
                        extensions--;
                        IQ--;
                        curProblem->isDropped = true;
                    }
                }

                curProblem->move();
                curProblem->setVX(vx);
                curProblem->setVY(vy);

            } //end outer for
        }
}

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

TEST_F(ProjectFixture,AnalysisRecord_AddAndRemoveDataPoints) {
  // create an analysis with data points, no results
  Analysis analysis("My Analysis",
                    Problem("My Problem",VariableVector(),runmanager::Workflow()),
                    FileReferenceType::OSM);
  Problem problem = analysis.problem();

  DiscretePerturbationVector perturbations;
  int pi = 1;
  std::stringstream ss;
  for (int i = 0; i < 3; ++i) {
    perturbations.push_back(NullPerturbation());
    for (int j = 0; j < 4; ++j) {
      ss << "measure" << pi << ".rb";
      perturbations.push_back(RubyPerturbation(toPath(ss.str()),
                                               FileReferenceType::OSM,
                                               FileReferenceType::OSM,true));
      ss.str("");
      ++pi;
    }
    ss << "Variable " << i+1;
    problem.push(DiscreteVariable(ss.str(),perturbations));
    perturbations.clear();
    ss.str("");
  }

  EXPECT_EQ(3u,analysis.problem().variables().size());
  ASSERT_FALSE(problem.variables().empty());
  EXPECT_EQ(5u,problem.variables()[0].cast<DiscreteVariable>().numPerturbations(true));

  problem.pushResponse(
        LinearFunction("Energy Use",
                       VariableVector(1u,OutputAttributeVariable("Energy Use","Total.Energy.Use"))));
  EXPECT_EQ(1u,problem.responses().size());

  std::vector<QVariant> values(3u,0);
  OptionalDataPoint dataPoint = problem.createDataPoint(values);
  ASSERT_TRUE(dataPoint);
  bool test = analysis.addDataPoint(*dataPoint);
  EXPECT_TRUE(test);
  values[1] = 3;
  dataPoint = problem.createDataPoint(values);
  ASSERT_TRUE(dataPoint);
  test = analysis.addDataPoint(*dataPoint);
  EXPECT_TRUE(test);
  values[0] = 1;
  values[1] = 2;
  values[2] = 4;
  dataPoint = problem.createDataPoint(values);
  ASSERT_TRUE(dataPoint);
  test = analysis.addDataPoint(*dataPoint);
  EXPECT_TRUE(test);

  EXPECT_TRUE(analysis.isDirty());
  EXPECT_TRUE(problem.isDirty());
  BOOST_FOREACH(const Variable& variable, problem.variables()) {
    EXPECT_TRUE(variable.isDirty());
    ASSERT_TRUE(variable.optionalCast<DiscreteVariable>());
    BOOST_FOREACH(const DiscretePerturbation& perturbation,variable.cast<DiscreteVariable>().perturbations(false)) {
      EXPECT_TRUE(perturbation.isDirty());
    }
  }
  BOOST_FOREACH(const Function& response, problem.responses()) {
    EXPECT_TRUE(response.isDirty());
  }
  BOOST_FOREACH(const DataPoint& dataPoint, analysis.dataPoints()) {
    EXPECT_TRUE(dataPoint.isDirty());
  }

  // save to database, make sure everything is there, make sure clean signal filters down
  ProjectDatabase database = getCleanDatabase("AnalysisRecord_AddAndRemoveDataPoints");
  {
    bool transactionStarted = database.startTransaction();
    EXPECT_TRUE(transactionStarted);
    AnalysisRecord analysisRecord(analysis,database);
    database.save();
    test = database.commitTransaction();
    EXPECT_TRUE(test);
    EXPECT_EQ(analysis.uuid(),analysisRecord.handle());
    EXPECT_EQ(analysis.versionUUID(),analysisRecord.uuidLast());
    ASSERT_NO_THROW(analysisRecord.problemRecord());
    ProblemRecord problemRecord = analysisRecord.problemRecord();
    InputVariableRecordVector variableRecords = problemRecord.inputVariableRecords();
    EXPECT_EQ(3u,variableRecords.size());
    ASSERT_FALSE(variableRecords.empty());
    ASSERT_TRUE(variableRecords[0].optionalCast<DiscreteVariableRecord>());
    DiscretePerturbationRecordVector perturbationRecords = variableRecords[0].cast<DiscreteVariableRecord>().discretePerturbationRecords(false);
    EXPECT_EQ(5u,perturbationRecords.size());
    ASSERT_TRUE(perturbationRecords.size() > 1);
    EXPECT_TRUE(perturbationRecords[0].optionalCast<NullPerturbationRecord>());
    EXPECT_TRUE(perturbationRecords[1].optionalCast<RubyPerturbationRecord>());
    EXPECT_EQ(3u,analysisRecord.dataPointRecords().size());
    EXPECT_TRUE(analysisRecord.completeDataPointRecords().empty());
  }

  analysis.clearDirtyFlag();
  EXPECT_FALSE(analysis.isDirty());
  EXPECT_FALSE(problem.isDirty());
  BOOST_FOREACH(const Variable& variable, problem.variables()) {
    EXPECT_FALSE(variable.isDirty());
//.........这里部分代码省略.........

    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++) {
//.........这里部分代码省略.........

  Result solve(const Problem& iProblem) const {
    // set up initial (empty) result
    Result result;
    result.mSuccess = false;
    result.mNumIterations = 0;

    // ensure that there are enough data points to proceed
    const int sampleSize = iProblem.getSampleSize();
    const int n = iProblem.getNumDataPoints();
    if (n < sampleSize) {
      return result;
    }

    const double epsilon = 1e-10;

    // best results are currently invalid
    double bestScore = -1;
    bool success = false;

    // start number of iterations as infinite, then reduce as we go
    double numIterationsNeeded = 1e10;
    int iterationCount = 0;
    int skippedSampleCount = 0;

    // for random sample index generation
    std::vector<int> allIndices(n);

    // iterate until adaptive number of iterations are exceeded
    while (iterationCount < numIterationsNeeded) {

      // determine random sample indices
      for (int i = 0; i < n; ++i) {
        allIndices[i] = i;
      }
      for (int i = 0; i < sampleSize; ++i) {
        int randIndex = std::rand() % n;
        std::swap(allIndices[i], allIndices[randIndex]);
      }
      std::vector<int> sampleIndices(allIndices.begin(),
                                     allIndices.begin() + sampleSize);

      // compute solution on minimal set
      typename Problem::Solution solution = iProblem.estimate(sampleIndices);

      // compute errors over all data points
      std::vector<double> errors2 = iProblem.computeSquaredErrors(solution);

      // check whether this is a valid sample
      // TODO: this should be done via a method in Problem class, but would
      // require changing all existing usages to include that method
      if (errors2.size() == 0) {
        ++skippedSampleCount;
        if (skippedSampleCount >= mMaximumIterations) break;
        continue;
      }
      skippedSampleCount = 0;

      // compute error threshold to be applied to each term
      double thresh = mMaximumError;
      if (thresh < 0) {
        std::sort(errors2.begin(), errors2.end());
        double median = (n % 2 == 0) ?
          (0.5*(errors2[n/2]+errors2[n/2+1])) : errors2[n/2];
        thresh = 1.4826*std::sqrt(median)*4.6851;
      }
      thresh *= thresh;

      // determine inliers
      std::vector<int> inliers;
      inliers.reserve(n);
      for (int i = 0; i < n; ++i) {
        if (errors2[i] <= thresh) {
          inliers.push_back(i);
        }
      }

      // if this is the best score, update solution and convergence criteria
      double score = inliers.size();
      if (score > bestScore) {
        bestScore = score;
        result.mInliers = inliers;
        result.mSolution = solution;
        success = true;
        double inlierProbability = double(inliers.size()) / n;
        double anyOutlierProbability = 1 - pow(inlierProbability,sampleSize);
        anyOutlierProbability = std::min(anyOutlierProbability, 1-epsilon);
        anyOutlierProbability = std::max(anyOutlierProbability, epsilon);
        numIterationsNeeded =
          log(1-mGoodSolutionProbability) / log(anyOutlierProbability);
      }

      // bump up iteration count and terminate if it exceeds hard max
      ++iterationCount;
      if (iterationCount > mMaximumIterations) {
        break;
      }
    }

    // finish off result params
    result.mSuccess = success;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            {
                dateiname = (char*) Malloc( (int) modelRoot.size() + 1 );
                dateiname = strcpy( dateiname, modelRoot.c_str() );
            }
        DisplayMsgLn(dateiname);
    }
    //WW  DisplayMsgLn("");
    //WW  DisplayMsgLn("");
    // ----------23.02.2009. WW-----------------

    // LB Check if file exists
    std::string tmpFilename = dateiname;
    tmpFilename.append(".pcs");
    if(!IsFileExisting(tmpFilename))
    {
        std::cout << " Error: Cannot find file " << dateiname << "\n";
        return 1;
    }

    // If no option is given, output files are placed in the same directory as the input files
    if (defaultOutputPath.empty()) defaultOutputPath = pathDirname(std::string(dateiname));

    FileName = dateiname;
    size_t indexChWin, indexChLinux;
    indexChWin = indexChLinux = 0;
    indexChWin = FileName.find_last_of('\\');
    indexChLinux = FileName.find_last_of('/');
    //
    if(indexChWin != std::string::npos)
        FilePath = FileName.substr(0,indexChWin) + "\\";
    else if(indexChLinux != std::string::npos)
        FilePath = FileName.substr(0,indexChLinux) + "/";
    // ---------------------------WW
    Problem* aproblem = new Problem(dateiname);
#ifdef USE_PETSC
    aproblem->setRankandSize(rank, r_size);
#endif
#if defined(USE_MPI) || defined(USE_MPI_PARPROC) || defined(USE_MPI_REGSOIL) || defined(USE_MPI_GEMS)  || defined(USE_MPI_KRC)
    aproblem->setRankandSize(myrank, mysize);

    if (myrank != MPI_UNDEFINED)
    {
#endif
        aproblem->Euler_TimeDiscretize();
        delete aproblem;
        aproblem = NULL;
#if defined(USE_MPI) || defined(USE_MPI_PARPROC) || defined(USE_MPI_REGSOIL) || defined(USE_MPI_GEMS)  || defined(USE_MPI_KRC)
    }

#ifdef OGS_FEM_IPQC
    //sending killing signals to ranks of group_IPQC, only when the group exists
    if (splitcomm_flag == true) {
        int signal = -1, rank_IPQC, mysize_IPQC = np - nb_ddc;
        for (int i=0; i< mysize_IPQC; i++) {
            rank_IPQC = mysize + i;
            MPI_Send(&signal, 1, MPI_INT, rank_IPQC, 0, MPI_COMM_WORLD);
        }
    }
#endif

#endif


    if(ClockTimeVec.size()>0)
        ClockTimeVec[0]->PrintTimes();  //CB time
    DestroyClockTime();

void TDefense::timerEvent(QTimerEvent *event)
{
    //if (isLaunched)
    //{
    std::cout << "\nWeek: " << week;

        launchModerator++;
        int round = week - 1;

        if (numEnemiesLaunched <= 10 + (round * 10))
        {
            if (launchModerator % 50 - (round * 5) == 0)
            {
                Problem *p = new Problem(gpWidth / 2, 0, problem);
                myProblems.push_back(p);
                numEnemiesLaunched++;
            }

            if (week > 1)
            {
                Tower *r = new Tower(500, 300, stapler);
                myTowers.push_back(r);
                if (launchModerator % 80 - (round * 10) == 0)
                {
                    Assignment *a = new Assignment(gpWidth / 2, 0, assignment1);
                    myAssignments.push_back(a);
                    numEnemiesLaunched++;
                } // end 'launchModerator' if
            } // end 'week' if

            if (week % 3 == 0)
            {
                if (launchModerator % 150 - (round * 15) == 0)
                {
                    Paper *r = new Paper(gpWidth / 2, 0, paper1);
                    myPapers.push_back(r);
                    numEnemiesLaunched++;
                } // end 'launchModerator' if
            } // 'week' if

            if (week == 5 || week == 9)
            {
                if (launchModerator % 300 - (round * 20) == 0)
                {
                    Midterm *m = new Midterm(gpWidth / 2, 0, midterm5);
                    myMidterms.push_back(m);
                    numEnemiesLaunched++;
                } // end 'launchModerator' if
            } // 'week' if
        } // end 'numEnemiesLaunched' if

        enemyAction();

        if (myProblems.size() > 0) {
            Problem *curProblem = myProblems.back();
            std::cout << "myProblems size: " << myProblems.size();
            std::cout << "\nproblem y-pos: " << curProblem->getY();
        }

        towerAction();
        if (studentShotFired)
            studentAction();

        if (extensions <= 0)
            gameOver();

        if ((numEnemiesShot + numEnemiesDropped) >= 10 + (round * 10))
        {
            clearEnemies();
            numEnemiesDropped = 0;
            numEnemiesLaunched = 1;
            numEnemiesShot = 0;
            isLaunched = false;
            IQ += 10;
            week++;
        }

        if (week == 10)
            semester++;
        if (semester == 2)
        {
            semester = 1;
            year++;
        }

        repaint();
    //} // end 'isLaunch' if
} // end method

void TDefense::paintEvent(QPaintEvent *event)
{
    QPainter painter(this);
    QColor* color = new QColor(255, 255, 255);
    painter.fillRect(*windowRect, *color);
    painter.drawRect(*windowRect);

    if (isOver) {

      QFont font("Courier", 15, QFont::DemiBold);
      QFontMetrics fm(font);
      int textWidth = fm.width("GAME OVER");

      painter.setFont(font);
      int h = height();
      int w = width();

      painter.translate(QPoint(w/2, h/2));
      painter.drawText(-textWidth/2, 0, "GAME OVER");
    } else if(isWon) {

      QFont font("Courier", 15, QFont::DemiBold);
      QFontMetrics fm(font);
      int textWidth = fm.width("VICTORY");

      painter.setFont(font);
      int h = height();
      int w = width();

      painter.translate(QPoint(w/2, h/2));
      painter.drawText(-textWidth/2, 0, "VICTORY");
    } else if (isPaused) {

        QFont font("Courier", 15, QFont::DemiBold);
        QFontMetrics fm(font);
        int textWidth = fm.width("PAUSED");

        painter.setFont(font);
        int h = height();
        int w = width();

        painter.translate(QPoint(w/2, h/2));
        painter.drawText(-textWidth/2, 0, "PAUSED");
    }
    else
    {

        QFont font("Brush Script MT", 15);
        painter.setFont(font);

        painter.translate(QPoint(30, 30));
        painter.drawText(0, 0, "WWeek: " + week);
        painter.translate(QPoint(13, 15));
        painter.drawText(0, 0, "SSemester: " + semester);
        painter.translate(QPoint(30, 15));
        painter.drawText(0, 0, "YYear: " + year);

        painter.translate(QPoint(-50, 300));
        painter.drawText(0, 0, "IQ: ");
        painter.translate(QPoint(10, 17));
        painter.drawText(0, 0, "Funds: " + funds);
        painter.translate(QPoint(10, 17));
        painter.drawText(0, 0, "EExtensionExtensions: " + extensions);


        painter.drawImage(myStudent->getRect()->center(), myStudent->getImage());
        /*for (int i = 0; i < myTowers.size(); i++)
        {
            Tower *curTower = myTowers[i];
            painter.drawImage(*curTower->getRect(), curTower->getImage());
        }*/

        //if (isLaunched)
        //{

            for (int i = 0; i < myProblems.size(); i++)
            {
                Problem *curProblem = myProblems[i];
                painter.drawImage(*curProblem->getRect(), curProblem->getImage());
            } // end for

            /*
            for (int i = 0; i < myAssignments.size(); i++)
            {
                Assignment *curAssignment = myAssignments[i];
                painter.drawImage(*curAssignment->getRect(), curAssignment->getImage());
            } // end for
            for (int i = 0; i < myPapers.size(); i++)
            {
                Paper *curPaper = myPapers[i];
                painter.drawImage(*curPaper->getRect(), curPaper->getImage());
            } // end for
            for (int i = 0; i < myMidterms.size(); i++)
            {
                Midterm *curMidterm = myMidterms[i];
                painter.drawImage(*curMidterm->getRect(), curMidterm->getImage());
            } // end for
            */

            if (studentShotFired)
//.........这里部分代码省略.........

// --------------------------------------------------------------------------------------------------------------------
// --------------------------------------------------------------------------------------------------------------------
// --------------------------------------------------------------------------------------------------------------------
ProblemFactory::ProblemFactory(bool _skipInit)
{
    Problem* newProb = nullptr;
    // Null
    newProb = new NullProblem();
    if (_skipInit || newProb->Init()) {
        if (!_skipInit) {
            newProb->Shutdown();
        }
        mProblems.push_back(newProb);
        mSolutions[mProblems.back()->GetName()].push_back(new NullSolution());
    } else {
        newProb->Shutdown();
        SafeDelete(newProb);
        console::error("Unable to create the Null Problem--exiting.");
    }

    // DynamicStreaming
    newProb = new DynamicStreamingProblem();
    if (_skipInit || newProb->Init()) {
        if (!_skipInit) {
            newProb->Shutdown();
        }
        mProblems.push_back(newProb);
        mSolutions[mProblems.back()->GetName()].push_back(new DynamicStreamingGLBufferSubData());
        mSolutions[mProblems.back()->GetName()].push_back(new DynamicStreamingGLMapUnsynchronized());
        mSolutions[mProblems.back()->GetName()].push_back(new DynamicStreamingGLMapPersistent());
        #if WITH_D3D11
            mSolutions[mProblems.back()->GetName()].push_back(new DynamicStreamingD3D11MapNoOverwrite());
            mSolutions[mProblems.back()->GetName()].push_back(new DynamicStreamingD3D11UpdateSubresource());
        #endif
    } else {
        newProb->Shutdown();
        SafeDelete(newProb);
    }

    // UntexturedObjects
    newProb = new UntexturedObjectsProblem();
    if (_skipInit || newProb->Init()) {
        if (!_skipInit) {
            newProb->Shutdown();
        }
        mProblems.push_back(newProb);
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLUniform());
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLDrawLoop());
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLMultiDraw(true));
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLMultiDraw(false));
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLMultiDrawBuffer(true));
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLMultiDrawBuffer(false));
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLBindless());
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLBindlessIndirect());
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLBufferRange());
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLBufferStorage(true));
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLBufferStorage(false));
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLDynamicBuffer());
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLMapUnsynchronized());
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLMapPersistent());
        mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsGLTexCoord());
        #if WITH_D3D11
            mSolutions[mProblems.back()->GetName()].push_back(new UntexturedObjectsD3D11Naive());
        #endif
    } else {
        newProb->Shutdown();
        SafeDelete(newProb);
    }

    // Textured Quads
    newProb = new TexturedQuadsProblem();
    if (_skipInit || newProb->Init()) {
        if (!_skipInit) {
            newProb->Shutdown();
        }
        mProblems.push_back(newProb);
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLBindless());
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLBindlessMultiDraw());
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLNaive());
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLNaiveUniform());
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLNoTex());
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLNoTexUniform());
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLSparseBindlessTextureArray());
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLSparseBindlessTextureArrayMultiDraw(true));
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLSparseBindlessTextureArrayMultiDraw(false));
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLTextureArray());
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLTextureArrayUniform());
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLTextureArrayMultiDraw(true));
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLTextureArrayMultiDraw(false));
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLTextureArrayMultiDrawBuffer(true));
        mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsGLTextureArrayMultiDrawBuffer(false));
        #if WITH_D3D11
            mSolutions[mProblems.back()->GetName()].push_back(new TexturedQuadsD3D11Naive());
        #endif

    } else {
        newProb->Shutdown();
        SafeDelete(newProb);
    }
}

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