std::vector<double> MultivariateFNormalSufficient::get_norms() const
{
    if (!flag_norms_)
    {
        //Eigen::LLT<MatrixXd, Eigen::Upper> ldlt(get_ldlt());
        Eigen::LDLT<MatrixXd, Eigen::Upper> ldlt(get_ldlt());
        // determinant and derived constants
        //MatrixXd L(ldlt.matrixU());
        //std::cout << "L: " << L << std::endl << std::endl;
        //std::cout << "D: " << ldlt.vectorD() << std::endl << std::endl;
        //double logDetSigma=2*L.diagonal().array().log().sum() ;
        double logDetSigma = ldlt.vectorD().array().abs().log().sum();
        IMP_LOG_TERSE( "MVN:   det(Sigma) = " <<exp(logDetSigma) << std::endl);
        double norm=pow(2*IMP::PI*IMP::square(factor_), -double(N_*M_)/2.0)
                    * exp(-double(N_)/2.0*logDetSigma);
        double lnorm=double(N_*M_)/2 * log(2*IMP::PI*IMP::square(factor_))
            + double(N_)/2 * logDetSigma;
        IMP_LOG_TERSE( "MVN:   norm = " << norm << " lnorm = "
                << lnorm << std::endl);
        const_cast<MultivariateFNormalSufficient *>(this)
            ->set_norms(norm,lnorm);
    }
    std::vector<double> norms;
    norms.push_back(norm_);
    norms.push_back(lnorm_);
    return norms;
}

IMPEM_BEGIN_NAMESPACE

EnvelopePenetrationRestraint::EnvelopePenetrationRestraint(
   Particles ps,
   DensityMap *em_map,Float threshold
   ): Restraint("Envelope penetration restraint")
{
  IMP_LOG_TERSE("Load envelope penetration with the following input:"<<
          "number of particles:"<<ps.size()<<
           "\n");
  threshold_=threshold;
  target_dens_map_ = em_map;
  IMP_IF_CHECK(USAGE) {
    for (unsigned int i=0; i< ps.size(); ++i) {
      IMP_USAGE_CHECK(core::XYZR::particle_is_instance(ps[i]),
                      "Particle " << ps[i]->get_name()
                      << " is not XYZR"
                      << std::endl);
    }
  }
  add_particles(ps);
  ps_=ps;
  IMP_LOG_TERSE("after adding particles"<<std::endl);
  IMP_LOG_TERSE( "Finish initialization" << std::endl);
}

ComplementarityRestraint::GridObject *ComplementarityRestraint::get_grid_object(
    core::RigidBody rb, const kernel::ParticlesTemp &a, ObjectKey ok,
    double thickness, double value, double interior_thickness,
    double voxel) const {
  IMP_USAGE_CHECK(!a.empty(), "No particles passed for excluded volume");
  for (unsigned int i = 1; i < a.size(); ++i) {
    IMP_USAGE_CHECK(core::RigidMember(a[0]).get_rigid_body() ==
                        core::RigidMember(a[i]).get_rigid_body(),
                    "Not all particles are from the same rigid body.");
  }
  if (!rb->has_attribute(ok)) {
    IMP_LOG_TERSE("Creating grid for rigid body " << rb->get_name()
                                                  << std::endl);
    IMP::algebra::DenseGrid3D<float> grid =
        get_grid(a, thickness, value, interior_thickness, voxel);
    IMP_LOG_TERSE("Grid has size " << grid.get_number_of_voxels(0) << ", "
                                   << grid.get_number_of_voxels(1) << ", "
                                   << grid.get_number_of_voxels(2)
                                   << std::endl);
    base::Pointer<GridObject> n(new GridObject(
        GridPair(rb.get_reference_frame().get_transformation_to(), grid)));
    rb->add_cache_attribute(ok, n);
  }
  IMP_CHECK_OBJECT(rb->get_value(ok));
  IMP_INTERNAL_CHECK(dynamic_cast<GridObject *>(rb->get_value(ok)),
                     "The saved grid is not a grid.");
  return dynamic_cast<GridObject *>(rb->get_value(ok));
}

void ProteomicsEMAlignmentAtomic::load_atomic_molecules(){
  IMP_LOG_TERSE("load atomic molecules \n");
    IMP_NEW(atom::ATOMPDBSelector,sel,());
  //  IMP_NEW(atom::CAlphaPDBSelector,sel,());
  for(int i=0;i<prot_data_->get_number_of_proteins();i++) {
    IMP_LOG_TERSE("going to load molecule "<<
            asmb_data_->get_component_header(i)->get_filename()<<"\n");
    //            prot_data_.get_protein_filename(i)<<"|\n";
    atom::Hierarchy mh =
      atom::read_pdb(asmb_data_->get_component_header(i)->get_filename(),
                     mdl_,sel);
    mh->set_name(asmb_data_->get_component_header(i)->get_name());
    mh->set_was_used(true);
    mhs_.push_back(mh);
    std::cout<<"create pdb"<<std::endl;
    std::cout<<"are subunits rigid?"
             <<params_.get_fragments_params().subunit_rigid_<<std::endl;
    if (params_.get_fragments_params().subunit_rigid_) {
      std::cout<<"create rigid body"<<std::endl;
      rbs_.push_back(atom::create_rigid_body(mh));
      rbs_[rbs_.size()-1]->set_name(mh->get_name());
      rbs_[rbs_.size()-1]->add_attribute(fit_state_key_,-1);
      rbs_[rbs_.size()-1]->add_attribute(order_key_,i);

    }
  }
}

void IncrementalScoringFunction::create_scoring_functions() {
  IMP_OBJECT_LOG;
  IMP_LOG_TERSE("Creating scoring functions" << std::endl);
  if (flattened_restraints_.empty()) return;

  boost::unordered_map<Restraint *, int> mp;
  IMP_LOG_TERSE("All restraints are " << flattened_restraints_ << std::endl);
  for (unsigned int i = 0; i < flattened_restraints_.size(); ++i) {
    mp[flattened_restraints_[i]] = i;
  }

  Restraints drs;
  for (unsigned int i = 0; i < nbl_.size(); ++i) {
    // This ensures that the score states needed for the non-bonded terms
    drs.push_back(nbl_[i]->get_dummy_restraint());
  }

  Vector<RestraintsTemp> crs;
  IMP_FOREACH(ParticleIndex pi, all_) {
    RestraintsTemp cr = get_dependent_restraints(get_model(), pi);
    /* Remove any duplicates in cr (could happen with rigid bodies) */
    std::sort(cr.begin(), cr.end());
    cr.erase(std::unique(cr.begin(), cr.end()), cr.end());
    crs.push_back(cr);
  }

void ProjectionFinder::set_projections(const em2d::Images &projections) {
  IMP_LOG_TERSE("ProjectionFinder: Setting projections" << std::endl);

  if(projections.size()==0) {
    IMP_THROW("Passing empty set of projections",ValueException);
  }

  if(polar_params_.get_is_setup() == false) {
    polar_params_.setup(projections[0]->get_data().rows,
                             projections[0]->get_data().cols);
    polar_params_.set_estimated_number_of_angles(
                    projections[0]->get_header().get_number_of_columns());
    polar_params_.create_maps_for_resampling();
  }

  projections_.resize(projections.size());
  unsigned int n_projections = projections_.size();
  PROJECTIONS_POLAR_AUTOC_.clear();
  PROJECTIONS_POLAR_AUTOC_.resize(n_projections);
  projections_cog_.resize(n_projections);
  boost::timer preprocessing_timer;
  for (unsigned int i = 0; i < n_projections; ++i) {
    projections_[i] = projections[i]; // does not copy
    std::ostringstream oss;
    oss << "Projection" << i;
    projections_[i]->set_name(oss.str());
    do_preprocess_projection(i);
  }
  preprocessing_time_ = preprocessing_timer.elapsed();
  IMP_LOG_TERSE("ProjectionFinder: Projections set: "
          << projections_.size() << std::endl);

}

double Em2DRestraint::unprotected_evaluate(DerivativeAccumulator *accum) const {
  IMP_UNUSED(accum);
  IMP_USAGE_CHECK(!accum, "No derivatives provided");
  IMP_NEW(Model, model, ());
  model = get_model();
  // Project the model
  RegistrationResults regs =
      get_evenly_distributed_registration_results(params_.n_projections);
  unsigned int rows = em_images_[0]->get_header().get_number_of_rows();
  unsigned int cols = em_images_[0]->get_header().get_number_of_columns();

  ProjectingOptions options(params_.pixel_size, params_.resolution);
  Images projections = get_projections(particles_container_->get_particles(),
                                       regs, rows, cols, options);
  finder_->set_projections(projections);

  if (only_coarse_registration_) {
    IMP_LOG_TERSE("Em2DRestraint::unprotected::evaluate: Coarse registration"
                  << std::endl);
    finder_->get_coarse_registration();
  } else {
    if (fast_optimization_mode_) {
      IMP_LOG_TERSE("Em2DRestraint::unprotected::evaluate: Fast Mode "
                    << number_of_optimized_projections_
                    << " projections optimized" << std::endl);

      finder_->set_fast_mode(number_of_optimized_projections_);
    }
    IMP_LOG_TERSE("Em2DRestraint::unprotected::evaluate: Complete registration"
                  << std::endl);
    finder_->get_complete_registration();
  }
  return finder_->get_global_score();
}

void ProjectionFinder::set_subjects(const em2d::Images &subjects) {
  IMP_LOG_TERSE("ProjectionFinder: Setting subject images" << std::endl);
  if (subjects.size() == 0) {
    IMP_THROW("Passing empty set of subjects", ValueException);
  }
  if (polar_params_.get_is_setup() == false) {
    polar_params_.setup(subjects[0]->get_data().rows,
                        subjects[0]->get_data().cols);
    polar_params_.set_estimated_number_of_angles(
        subjects[0]->get_header().get_number_of_columns());
    polar_params_.create_maps_for_resampling();
  }
  boost::timer preprocessing_timer;
  subjects_.resize(subjects.size());
  unsigned int n_subjects = subjects_.size();
  registration_results_.clear();
  registration_results_.resize(n_subjects);
  SUBJECTS_.clear();
  SUBJECTS_.resize(n_subjects);
  SUBJECTS_POLAR_AUTOC_.clear();
  SUBJECTS_POLAR_AUTOC_.resize(n_subjects);
  subjects_cog_.resize(n_subjects);
  for (unsigned int i = 0; i < n_subjects; ++i) {
    subjects_[i] = subjects[i];  // doesn't not deep copy
    std::ostringstream oss;
    oss << "Image subject " << i;
    subjects_[i]->set_name(oss.str());
    subjects_[i]->set_was_used(true);
    do_preprocess_subject(i);
  }
  preprocessing_time_ = preprocessing_timer.elapsed();
  IMP_LOG_TERSE("ProjectionFinder: Subject images set" << std::endl);
}

void RestraintCache::add_restraint_internal(
    kernel::Restraint *r, unsigned int index, kernel::RestraintSet *parent,
    double parent_max, Subset parent_subset, const DepMap &dependencies) {
  IMP_OBJECT_LOG;
  IMP_LOG_TERSE("Processing " << Showable(r) << " with " << parent_max
                              << std::endl);
  r->set_was_used(true);
  // fix using PST
  Subset cur_subset = get_subset(r, dependencies);
  double cur_max = r->get_maximum_score();
  if (parent) {
    cur_max = std::min(parent_max / r->get_weight(), cur_max);
  }

  if (cur_max < std::numeric_limits<double>::max()) {
    IMP_LOG_TERSE("Adding restraint " << Showable(r) << " with max " << cur_max
                                      << " and subset " << cur_subset
                                      << std::endl);
    known_restraints_[r] = cur_subset;
    restraint_index_[r] = index;
  }
  add_restraint_set_child_internal(r, cur_subset, parent, parent_max,
                                   parent_subset);
  kernel::RestraintSet *rs = dynamic_cast<kernel::RestraintSet *>(r);
  if (rs) {
    add_restraint_set_internal(rs, index, cur_subset, cur_max, dependencies);
  } else {
    if (cur_max < std::numeric_limits<double>::max()) {
      cache_.access_generator().add_restraint(r, cur_subset, cur_max);
    }
  }
}

void GaussianProcessInterpolationScoreState::do_before_evaluate() {
  IMP_LOG_TERSE("GPISS: do_before_evaluate()" << std::endl);
  GaussianProcessInterpolation *gpi_;
  gpi_ = gpir_->gpi_;
  MultivariateFNormalSufficient *mvn_;
  mvn_ = gpir_->mvn_;
  //
  gpi_->update_flags_covariance();  // O(1)
  gpi_->update_flags_mean();        // O(1)
  if (!(gpi_->flag_m_gpir_))        // gpi says that our m is not up to date
  {
    mvn_->set_FM(gpi_->get_m());  // O(M_)
    gpi_->flag_m_gpir_ = true;
    IMP_LOG_TERSE(" updated mean");
  }
  if (!(gpi_->flag_Omega_gpir_)) {
    mvn_->set_Sigma(gpi_->get_Omega());  // O(M^2)
    gpi_->flag_Omega_gpir_ = true;
    IMP_LOG_TERSE(" updated covariance");
  }
  IMP_LOG_TERSE(std::endl);
  /*ParticlesTemp tmp(gpir_->get_input_particles());
  std::cout << "values: ";
  for (unsigned i=0; i<tmp.size(); i++)
  {
      std::cout << i << " = " << Scale(tmp[i]).get_nuisance() << " ";
  }
  std::cout <<std::endl;*/
}

double
get_diffusion_coefficient(const algebra::Vector3Ds &displacements,
                          double dt) {
  algebra::Vector3D Ds;
  for (unsigned int i=0; i< 3; ++i) {
    Ds[i]= get_diffusion_coefficient(displacements.begin(),
                                     displacements.end(), i,  dt);
  }
  IMP_LOG_TERSE( "Diffusion coefficients are " << Ds << std::endl);
  int len=displacements.size()/2;
  algebra::Vector3D Ds0;
  for (unsigned int i=0; i< 3; ++i) {
    Ds0[i]= get_diffusion_coefficient(displacements.begin(),
                                      displacements.begin()+len, i, dt);
  }
  algebra::Vector3D Ds1;
  for (unsigned int i=0; i< 3; ++i) {
    Ds1[i]= get_diffusion_coefficient(displacements.begin()+len,
                                      displacements.end(), i, dt);
  }
  IMP_LOG_TERSE( "Partial coefficients are " << Ds0 << " and "
          << Ds1 << std::endl);
  return std::accumulate(Ds1.coordinates_begin(),
                         Ds1.coordinates_end(), 0.0)/3.0;
}

/** \param[in] acc If true (not nullptr), partial first derivatives should be
                          calculated.
    \return score associated with this restraint for the given state of
            the model.
*/
double Restraint::unprotected_evaluate(DerivativeAccumulator* acc) const {

  IMP_LOG_TERSE("SAXS Restraint::evaluate score\n");

  IMP_NEW(Profile, model_profile, ());
  handler_->compute_profile(model_profile);
  double score = profile_fitter_->compute_score(model_profile);
  bool calc_deriv = acc ? true : false;
  if (!calc_deriv) return score;

  IMP_LOG_TERSE("SAXS Restraint::compute derivatives\n");

  // do we need to resample the curve since it's just been created??
  // yes, since it does not correspond to the experimental one
  IMP_NEW(Profile, resampled_profile, ());
  profile_fitter_->resample(model_profile, resampled_profile);

  Vector<double> effect_size;  // Gaussian model-specific derivative weights
  double offset = 0.0;
  double c = profile_fitter_->compute_scale_factor(model_profile);
  derivative_calculator_->compute_gaussian_effect_size(model_profile, c, offset,
                                                       effect_size);

  handler_->compute_derivatives(derivative_calculator_, model_profile,
                                effect_size, acc);

  IMP_LOG_TERSE("SAXS Restraint::done derivatives, score " << score
                                                                   << "\n");
  return score;
}

MatrixXd MultivariateFNormalSufficient::get_PW() const
{
    if (!flag_PW_)
    {
        ////PW
        timer_.start(SOLVE);
        MatrixXd PW(M_,M_);
        if (N_==1)
        {
            IMP_LOG_TERSE( "MVN:   W=0 => PW=0" << std::endl);
            PW.setZero();
        } else {
            IMP_LOG_TERSE( "MVN:   solving for PW" << std::endl);
            if (use_cg_)
            {
                if (first_PW_)
                {
                    PW = compute_PW_direct();
                    (*const_cast<bool *>(&first_PW_))=false;
                } else {
                    PW = compute_PW_cg();
                }
            } else {
                PW = compute_PW_direct();
            }
        }
        const_cast<MultivariateFNormalSufficient *>(this)->set_PW(PW);
        timer_.stop(SOLVE);
    }
    return PW_;
}

IMPISD_BEGIN_NAMESPACE

GaussianProcessInterpolationRestraintSparse::
    GaussianProcessInterpolationRestraintSparse(
        GaussianProcessInterpolationSparse *gpi) : gpi_(gpi)
{
    //number of observation points
    IMP_LOG_TERSE( "GPIR: init" << std::endl);
    M_ = gpi_->n_obs_.size();
    //number of repetitions
    int N=gpi_->n_obs_[0];
    //check if the number of repetitions is the same for every point
    IMP_IF_CHECK(USAGE) {
        for (unsigned i=1; i<M_; i++)
            IMP_USAGE_CHECK(N == gpi_->n_obs_[i],
                "must have the same number of repetitions for each point!");
    }
    // build multivariate normal with
    // mean : prior mean
    // covariance : prior covariance
    // observed at : the original observations
    IMP_LOG_TERSE( "GPIR: multivariate normal()" << std::endl);
    //args are: sample mean, jacobian, true mean,
    // nobs, sample variance, true variance
    c_ = gpi_->get_cholmod_common();
    mvn_ = new MultivariateFNormalSufficientSparse(gpi_->get_I(), 1.0,
            gpi_->get_m(), N, gpi_->get_S(), gpi_->get_W(), c_);
    IMP_LOG_TERSE( "GPIR: done init" << std::endl);
}

void GaussianProcessInterpolation::compute_OmiIm() {
  IMP_Eigen::VectorXd I(get_I());
  IMP_Eigen::VectorXd m(get_m());
  IMP_Eigen::MatrixXd Omi(get_Omi());
  IMP_LOG_TERSE("OmiIm ");
  OmiIm_ = ldlt_.solve(I - m);
  IMP_LOG_TERSE(std::endl);
}

void GaussianProcessInterpolation::compute_I(Floats mean) {
  I_ = IMP_Eigen::VectorXd(M_);
  IMP_LOG_TERSE("I: ");
  for (unsigned i = 0; i < M_; i++) {
    I_(i) = mean[i];
    IMP_LOG_TERSE(I_(i) << " ");
  }
  IMP_LOG_TERSE(std::endl);
}

void read_pdb(Model *model, const std::string file,
              std::vector<std::string>& pdb_file_names,
              std::vector<IMP::Particles>& particles_vec,
              bool residue_level, bool heavy_atoms_only, int multi_model_pdb,
              bool explicit_water) {
  IMP::atom::Hierarchies mhds;
  IMP::atom::PDBSelector* selector;
  if (residue_level)  // read CA only
    selector = new IMP::atom::CAlphaPDBSelector();
  else if (heavy_atoms_only) {  // read without hydrogens
    if (explicit_water)
      selector = new IMP::atom::NonHydrogenPDBSelector();
    else
      selector = new IMP::atom::NonWaterNonHydrogenPDBSelector();
  } else { // read with hydrogens
    if (explicit_water)
      selector = new IMP::atom::NonAlternativePDBSelector();
    else
      selector = new IMP::atom::NonWaterPDBSelector();
  }

  if (multi_model_pdb == 2) {
    mhds = read_multimodel_pdb(file, model, selector, true);
  } else {
    if (multi_model_pdb == 3) {
      IMP::atom::Hierarchy mhd =
          IMP::atom::read_pdb(file, model, selector, false, true);
      mhds.push_back(mhd);
    } else {
      IMP::atom::Hierarchy mhd =
          IMP::atom::read_pdb(file, model, selector, true, true);
      mhds.push_back(mhd);
    }
  }

  for (unsigned int h_index = 0; h_index < mhds.size(); h_index++) {
    IMP::ParticlesTemp ps =
        get_by_type(mhds[h_index], IMP::atom::ATOM_TYPE);
    if (ps.size() > 0) {  // pdb file
      std::string pdb_id = file;
      if (mhds.size() > 1) {
        pdb_id = trim_extension(file) + "_m" +
                 std::string(boost::lexical_cast<std::string>(h_index + 1)) +
                 ".pdb";
      }
      pdb_file_names.push_back(pdb_id);
      particles_vec.push_back(IMP::get_as<IMP::Particles>(ps));
      if (mhds.size() > 1) {
        IMP_LOG_TERSE(ps.size() << " atoms were read from PDB file " << file
                      << " MODEL " << h_index + 1 << std::endl);
      } else {
        IMP_LOG_TERSE(ps.size() << " atoms were read from PDB file " << file
                      << std::endl);
      }
    }
  }
}

double ComplementarityRestraint::unprotected_evaluate_if_good(
    DerivativeAccumulator *accum, double max) const {
  IMP_OBJECT_LOG;
  IMP_USAGE_CHECK_VARIABLE(accum);
  IMP_USAGE_CHECK(!accum,
                  "ComplementarityRestraint does not support derivatives.");
  double vol = cube(voxel_size_);

  internal::ComplementarityParameters params;
  params.maximum_separation = maximum_separation_;
  params.maximum_penetration_score =
      std::min(maximum_penetration_score_ / vol, max);
  // std::cout<<"max penet score:"<<params.maximum_penetration_score<<"(" <<
  // maximum_penetration_score_<<","<<vol<<","<<max<<")"<<std::endl;
  base::Pointer<GridObject> ga =
      get_grid_object(rba_, a_, ok_, complementarity_thickness_,
                      complementarity_value_, interior_thickness_, voxel_size_);
  base::Pointer<GridObject> gb =
      get_grid_object(rbb_, b_, ok_, complementarity_thickness_,
                      complementarity_value_, interior_thickness_, voxel_size_);
  algebra::Transformation3D tra =
      ga->get_data().first *
      rba_.get_reference_frame().get_transformation_from();
  algebra::Transformation3D trb =
      rbb_.get_reference_frame().get_transformation_to() / gb->get_data().first;
  // transform a by tra and b by trb
  // same as transforming b by na/oa Ma= oa/ nai nb/ob p
  algebra::Transformation3D tr = tra * trb;
  IMP_LOG_TERSE("Transformation is " << tr << " between "
                                     << rba_.get_reference_frame() << " and "
                                     << rbb_.get_reference_frame()
                                     << std::endl);

  IMP::multifit::internal::FitScore ps =
      IMP::multifit::internal::get_fit_scores(
          ga->get_data().second, gb->get_data().second, tr, params);
  IMP_LOG_TERSE("Scores are " << ps.penetration_score << ", "
                              << ps.complementarity_score << " and "
                              << ps.boundary_score << std::endl);
  /*  std::cout<<"Scores are " << ps.penetration_score << ", "
      << ps.complementarity_score << " and "
      << ps.boundary_score
      << std::endl;*/
  if (!score_acceptable(ps)) {
    //    std::cout<<"scores are not acceptable"<<std::endl;
    return std::numeric_limits<double>::max();
  }
  double score = penetration_coef_ * ps.penetration_score +
                 complementarity_coef_ * ps.complementarity_score +
                 boundary_coef_ * ps.boundary_score;
  return score * vol;
}

void GaussianProcessInterpolation::compute_S(Floats std) {
  // if you modify this routine so that
  // S is not diagonal check the GPIR to make sure it still needs
  // to call set_W_nonzero of MVN.
  IMP_Eigen::VectorXd v(M_);
  IMP_LOG_TERSE("S: ");
  for (unsigned i = 0; i < M_; i++) {
    v(i) = IMP::square(std[i]);
    IMP_LOG_TERSE(v(i) << " ");
  }
  S_ = v.asDiagonal();
  IMP_LOG_TERSE(std::endl);
}

double MonteCarlo::do_optimize(unsigned int max_steps) {
    IMP_OBJECT_LOG;
    IMP_CHECK_OBJECT(this);
    if (get_number_of_movers() == 0) {
        IMP_THROW("Running MonteCarlo without providing any"
                  << " movers isn't very useful.",
                  ValueException);
    }

    ParticleIndexes movable = get_movable_particles();

    // provide a way of feeding in this value
    last_energy_ = do_evaluate(movable);
    if (return_best_) {
        best_ = new Configuration(get_model());
        best_energy_ = last_energy_;
    }
    reset_statistics();
    update_states();

    IMP_LOG_TERSE("MC Initial energy is " << last_energy_ << std::endl);

    for (unsigned int i = 0; i < max_steps; ++i) {
        if (get_stop_on_good_score() &&
                get_scoring_function()->get_had_good_score()) {
            break;
        }
        do_step();
        if (best_energy_ < get_score_threshold()) break;
    }

    IMP_LOG_TERSE("MC Final energy is " << last_energy_ << std::endl);
    if (return_best_) {
        // std::cout << "Final score is " << get_model()->evaluate(false)
        //<< std::endl;
        best_->swap_configuration();
        IMP_LOG_TERSE("MC Returning energy " << best_energy_ << std::endl);
        IMP_IF_CHECK(USAGE) {
            IMP_LOG_TERSE("MC Got " << do_evaluate(get_movable_particles())
                          << std::endl);
            /*IMP_INTERNAL_CHECK((e >= std::numeric_limits<double>::max()
                                && best_energy_ >= std::numeric_limits<double>::max())
                               || std::abs(best_energy_ - e)
                               < .01+.1* std::abs(best_energy_ +e),
                               "Energies do not match "
                               << best_energy_ << " vs " << e << std::endl);*/
        }

        return do_evaluate(movable);
    } else {
        return last_energy_;