自学教程：C++ IMP_LOG_TERSE函数代码示例

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_NAMESPACEEnvelopePenetrationRestraint::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;*/}

doubleget_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_NAMESPACEGaussianProcessInterpolationRestraintSparse::    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_;

void MovedSingletonContainer::do_score_state_before_evaluate() {  IMP_OBJECT_LOG;  IMP_CHECK_OBJECT(this);  IMP_CHECK_OBJECT(pc_);  if (pc_version_ != pc_->get_contents_hash()) {    pc_version_ = pc_->get_contents_hash();    IMP_LOG_TERSE("First call" << std::endl);    initialize();  } else {    ParticleIndexes mved = do_get_moved();    IMP_LOG_TERSE("Setting moved list: " << Showable(mved) << std::endl);    swap(mved);  }  IMP_IF_CHECK(USAGE_AND_INTERNAL) { validate(); }}

display::Geometries get_rigid_body_derivative_geometries(    Model *m, ParticleIndex pi) {    RigidBody d(m, pi);    display::Geometries ret;    Particles ms = get_as<Particles>(d.get_rigid_members());    algebra::Transformation3D otr =        d.get_reference_frame().get_transformation_to();    algebra::VectorD<4> rderiv = d.get_rotational_derivatives();    algebra::Vector3D tderiv = d.get_derivatives();    algebra::VectorD<4> rot = otr.get_rotation().get_quaternion();    IMP_LOG_TERSE("Old rotation was " << rot << std::endl);    Float scale = .1;    algebra::VectorD<4> dv = rderiv;    if (dv.get_squared_magnitude() > 0.00001) {        dv = scale * dv.get_unit_vector();    }    rot += dv;    rot = rot.get_unit_vector();    algebra::Rotation3D r(rot[0], rot[1], rot[2], rot[3]);    IMP_LOG_TERSE("Derivative was " << tderiv << std::endl);    IMP_LOG_TERSE("New rotation is " << rot << std::endl);    FloatRange xr =        d.get_particle()->get_model()->get_range(core::XYZ::get_xyz_keys()[0]);    Float wid = xr.second - xr.first;    algebra::Vector3D stderiv = scale * tderiv * wid;    algebra::Transformation3D ntr(        algebra::Rotation3D(rot[0], rot[1], rot[2], rot[3]),        stderiv + otr.get_translation());    for (unsigned int i = 0; i < ms.size(); ++i) {        core::RigidMember dm(ms[i]);        display::SegmentGeometry *tr = new display::SegmentGeometry(            algebra::Segment3D(dm.get_coordinates(), dm.get_coordinates() + tderiv),            /*xyzcolor_*/            display::Color(1, 0, 0));        ret.push_back(tr);        algebra::Vector3D ic =            r.get_rotated(dm.get_internal_coordinates()) + d.get_coordinates();        display::SegmentGeometry *rtr = new display::SegmentGeometry(            algebra::Segment3D(dm.get_coordinates(), ic), display::Color(0, 1, 0));        ret.push_back(rtr);        display::SegmentGeometry *nrtr = new display::SegmentGeometry(            algebra::Segment3D(dm.get_coordinates(),                               ntr.get_transformed(dm.get_internal_coordinates())),            display::Color(0, 0, 1));        ret.push_back(nrtr);    }    return ret;}

void Fine2DRegistrationRestraint::setup(    ParticlesTemp &ps, const ProjectingParameters &params,    Model *scoring_model,    //                       ScoreFunctionPtr score_function,    ScoreFunction *score_function, MasksManagerPtr masks) {  IMP_LOG_TERSE("Initializing Fine2DRegistrationRestraint" << std::endl);  ps_ = ps;  params_ = params;  // Generate all the projection masks for the structure  if (masks == MasksManagerPtr()) {    // Create the masks    masks_ =        MasksManagerPtr(new MasksManager(params.resolution, params.pixel_size));    masks_->create_masks(ps);    IMP_LOG_VERBOSE("Created " << masks_->get_number_of_masks()                               << " masks withing Fine2DRegistrationRestraint "                               << std::endl);  } else {    masks_ = masks;    IMP_LOG_VERBOSE("masks given to Fine2DRegistrationRestraint " << std::endl);  }  // Create a particle for the projection parameters to be optimized  subj_params_particle_ = new Particle(scoring_model);  PP_ = ProjectionParameters::setup_particle(subj_params_particle_);  PP_.set_parameters_optimized(true);  // Add an score state to the model  IMP_NEW(ProjectionParametersScoreState, pp_score_state,          (subj_params_particle_));  scoring_model->add_score_state(pp_score_state);  score_function_ = score_function;}

ResultAlign2D get_rotational_alignment(const cv::Mat &input,                          cv::Mat &m_to_align,bool apply) {  IMP_LOG_TERSE(          "starting 2D rotational alignment" << std::endl);  IMP_USAGE_CHECK((input.rows==m_to_align.rows) &&                  (input.cols==m_to_align.cols),                  "em2d::align_rotational: Matrices have different size.");  cv::Mat polar1,polar2,corr;  // Build maps for resampling  PolarResamplingParameters polar_params(input.rows,input.cols);  polar_params.set_estimated_number_of_angles(std::min(input.rows,input.cols));  polar_params.create_maps_for_resampling();  do_resample_polar(input,polar1,polar_params);  // subject image  do_resample_polar(m_to_align,polar2,polar_params); // projection image  ResultAlign2D RA= get_translational_alignment(polar1,polar2);  algebra::Vector2D shift=RA.first.get_translation();  // column shift[0] is the optimal angle to bring m_to_align to input  double angle =shift[0]*polar_params.get_angle_step();  RA.first.set_rotation(angle);  RA.first.set_translation(algebra::Vector2D(0.0,0.0));  // Apply the rotation if requested  if(apply) {    cv::Mat result;    get_transformed(m_to_align,result,RA.first);    result.copyTo(m_to_align);  }  IMP_LOG_VERBOSE("Rotational alingment: Transformation= "          << RA.first << " cross_correlation = " << RA.second << std::endl);  return RA;}