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

IMPEM_BEGIN_NAMESPACEfloat CoarseCC::calc_score(DensityMap *em_map, SampledDensityMap *model_map,                           float scalefac, bool recalc_rms, bool resample,                           FloatPair norm_factors) {  // resample the map for the particles provided  if (resample) {    model_map->resample();  }  if (recalc_rms) {    em_map->calcRMS();    // determine a threshold for calculating the CC    model_map->calcRMS();  // This function adequately computes the dmin value,                           // the safest value for the threshold  }  emreal voxel_data_threshold = model_map->get_header()->dmin - EPS;  // rmss have already been calculated  float escore = cross_correlation_coefficient(      em_map, model_map, voxel_data_threshold, false, norm_factors);  IMP_LOG_VERBOSE("CoarseCC::evaluate parameters:  threshold:"                  << voxel_data_threshold << std::endl);  IMP_LOG_VERBOSE("CoarseCC::evaluate: the score is:" << escore << std::endl);  escore = scalefac * (1. - escore);  return escore;}

SettingsData *read_settings(const char *filename) {  std::fstream in;  in.open(filename, std::fstream::in);  if (!in.good()) {    IMP_THROW("Problem opening file " << filename << " for reading "                                      << std::endl,              ValueException);  }  std::string line;  IMP_NEW(SettingsData, header, ());  getline(in, line);  // skip header line  int status = 0;  while (!in.eof()) {    getline(in, line);  // skip header line    std::vector<std::string> line_split;    boost::split(line_split, line, boost::is_any_of("|"));    if ((line_split.size() == 10) && (status == 0)) {  // protein  line      IMP_LOG_VERBOSE("parsing component line:" << line << std::endl);      header->add_component_header(parse_component_line(filename, line));    } else if (status == 0) {  // map header line      status = 1;    } else if (status == 1) {  // map line      IMP_LOG_VERBOSE("parsing EM line:" << line << std::endl);      header->set_assembly_header(parse_assembly_line(filename, line));      status = 2;    } else if (line.length() > 0) {  // don't warn about empty lines      IMP_WARN("the line was not parsed:" << line << "| with status:" << status                                          << std::endl);    }  }  in.close();  header->set_assembly_filename(filename);  header->set_data_path(".");  return header.release();}

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

void get_segmentation(em::DensityMap *dmap, double apix,                      double density_threshold, int num_means,                      const std::string pdb_filename,                      const std::string cmm_filename,                      const std::string seg_filename,                      const std::string txt_filename) {  IMP_LOG_VERBOSE("start setting trn_em" << std::endl);  IMP_NEW(DensityDataPoints, ddp, (dmap, density_threshold));  IMP_LOG_VERBOSE("initialize calculation of initial centers" << std::endl);  statistics::internal::VQClustering vq(ddp, num_means);  vq.run();  DataPointsAssignment assignment(ddp, &vq);  AnchorsData ad(assignment.get_centers(), *(assignment.get_edges()));  multifit::write_pdb(pdb_filename, assignment);  // also write cmm string into a file:  if (cmm_filename != "") {    write_cmm(cmm_filename, "anchor_graph", ad);  }  if (seg_filename != "") {    write_segments_as_mrc(dmap, assignment, apix, apix, density_threshold,                          seg_filename);  }  if (txt_filename != "") {    write_txt(txt_filename, ad);  }}

void FitRestraint::resample() const {  // TODO - first check that the bounding box of the particles  // matches that of the sampled ones.  // resample the map containing all non rigid body particles  // this map has all of the non rigid body particles.  if (not_part_of_rb_.size() > 0) {    none_rb_model_dens_map_->resample();    none_rb_model_dens_map_->calcRMS();    model_dens_map_->copy_map(none_rb_model_dens_map_);  } else {    model_dens_map_->reset_data(0.);  }  for (unsigned int rb_i = 0; rb_i < rbs_.size(); rb_i++) {    IMP_LOG_VERBOSE("Rb model dens map size:"                    << get_bounding_box(rb_model_dens_map_[rb_i], -1000.)                    << "/n Target size:"                    << get_bounding_box(target_dens_map_, -1000.) << "/n");    algebra::Transformation3D rb_t =        algebra::get_transformation_from_first_to_second(            rbs_orig_rf_[rb_i], rbs_[rb_i].get_reference_frame());    Pointer<DensityMap> transformed =        get_transformed(rb_model_dens_map_[rb_i], rb_t);    IMP_LOG_VERBOSE("transformed map size:"                    << get_bounding_box(transformed, -1000.) << std::endl);    model_dens_map_->add(transformed);    transformed->set_was_used(true);  }}

//TODO - do not use ProteinsAnchorsSamplingSpace.//you are not going to use the paths hereProteomicsEMAlignmentAtomic::ProteomicsEMAlignmentAtomic(                   const ProteinsAnchorsSamplingSpace &mapping_data,                   multifit::SettingsData *asmb_data,                   const AlignmentParams &align_param) :  base::Object("ProteomicsEMAlignmentAtomic%1%"),  mapping_data_(mapping_data),  params_(align_param),  order_key_(IntKey("order")),  asmb_data_(asmb_data){  fast_scoring_=false;  std::cout<<"start"<<std::endl;  std::cout<<"here0.2/n";  //initialize everything  mdl_=new Model();  IMP_LOG_VERBOSE("get proteomics data/n");  std::cout<<"get proteomics data/n";  prot_data_=mapping_data_.get_proteomics_data();  fit_state_key_ = IntKey("fit_state_key");  load_atomic_molecules();  std::cout<<"here1"<<std::endl;  IMP_LOG_VERBOSE("set NULL /n");  pst_=nullptr;  restraints_set_=false;states_set_=false;filters_set_=false;  ev_thr_=0.001;//TODO make a parameter  IMP_LOG_VERBOSE("end initialization/n");}

void KMCentersNodeSplit::get_neighbors(const Ints &cands, KMPointArray *sums,                                       KMPoint *sum_sqs, Ints *weights) {  if (cands.size() == 1) {    IMP_LOG_VERBOSE("KMCentersNodeSplit::get_neighbors the data points are"                    << " associated to center : " << cands[0] << std::endl);    // post points as neighbors    post_neighbor(sums, sum_sqs, weights, cands[0]);  }      // get cloest candidate to the box represented by the node      else {    Ints new_cands;    IMP_LOG_VERBOSE(        "KMCentersNodeSplit::get_neighbors compute close centers for node:/n");    IMP_LOG_WRITE(VERBOSE, show(IMP_STREAM));    compute_close_centers(cands, &new_cands);    for (unsigned int i = 0; i < new_cands.size(); i++) {      IMP_LOG_VERBOSE(new_cands[i] << "  | ");    }    IMP_LOG_VERBOSE("/nKMCentersNodeSplit::get_neighbors call left child with "                    << new_cands.size() << " candidates/n");    children_[0]->get_neighbors(new_cands, sums, sum_sqs, weights);    IMP_LOG_VERBOSE("KMCentersNodeSplit::get_neighbors call right child with "                    << new_cands.size() << " candidates/n");    children_[1]->get_neighbors(new_cands, sums, sum_sqs, weights);  }}

void get_projection(em2d::Image *img, const ParticlesTemp &ps,                    const RegistrationResult &reg,                    const ProjectingOptions &options, MasksManagerPtr masks,                    String name) {  IMP_LOG_VERBOSE("Generating projection in a em2d::Image" << std::endl);  if (masks == MasksManagerPtr()) {    masks = MasksManagerPtr(        new MasksManager(options.resolution, options.pixel_size));    masks->create_masks(ps);    IMP_LOG_VERBOSE("Masks generated from get_projection()" << std::endl);  }  algebra::Vector3D translation = options.pixel_size * reg.get_shift_3d();  algebra::Rotation3D R = reg.get_rotation();  do_project_particles(ps, img->get_data(), R, translation, options, masks);  if (options.normalize) em2d::do_normalize(img, true);  reg.set_in_image(img->get_header());  img->get_header().set_object_pixel_size(options.pixel_size);  if (options.save_images) {    if (name.empty()) {      IMP_THROW("get_projection: File name string is empty ", IOException);    }    if (options.srw == Pointer<ImageReaderWriter>()) {      IMP_THROW(          "The options class does not have an "          "ImageReaderWriter assigned. Create an ImageReaderWriter "          "and assigned to the srw member of ProjectingOptions.",          IOException);    }    img->write(name, options.srw);  }}

InteractionGraph get_interaction_graph(ScoringFunctionAdaptor rsi,                                       const ParticlesTemp &ps) {  if (ps.empty()) return InteractionGraph();  InteractionGraph ret(ps.size());  Restraints rs =      create_decomposition(rsi->create_restraints());  // Model *m= ps[0]->get_model();  boost::unordered_map<ModelObject *, int> map;  InteractionGraphVertexName pm = boost::get(boost::vertex_name, ret);  DependencyGraph dg = get_dependency_graph(ps[0]->get_model());  DependencyGraphVertexIndex index = IMP::get_vertex_index(dg);  /*IMP_IF_LOG(VERBOSE) {    IMP_LOG_VERBOSE( "dependency graph is /n");    IMP::internal::show_as_graphviz(dg, std::cout);    }*/  for (unsigned int i = 0; i < ps.size(); ++i) {    ParticlesTemp t = get_dependent_particles(        ps[i], ParticlesTemp(ps.begin(), ps.end()), dg, index);    for (unsigned int j = 0; j < t.size(); ++j) {      IMP_USAGE_CHECK(map.find(t[j]) == map.end(),                      "Currently particles which depend on more "                          << "than one particle "                          << "from the input set are not supported."                          << "  Particle /"" << t[j]->get_name()                          << "/" depends on /"" << ps[i]->get_name()                          << "/" and /""                          << ps[map.find(t[j])->second]->get_name() << "/"");      map[t[j]] = i;    }    IMP_IF_LOG(VERBOSE) {      IMP_LOG_VERBOSE("Particle /"" << ps[i]->get_name() << "/" controls ");      for (unsigned int i = 0; i < t.size(); ++i) {        IMP_LOG_VERBOSE("/"" << t[i]->get_name() << "/" ");      }      IMP_LOG_VERBOSE(std::endl);    }    pm[i] = ps[i];  }  IMP::Restraints all_rs = IMP::get_restraints(rs);  for (Restraints::const_iterator it = all_rs.begin();       it != all_rs.end(); ++it) {    ModelObjectsTemp pl = (*it)->get_inputs();    add_edges(ps, pl, map, *it, ret);  }  /* Make sure that composite score states (eg the normalizer for     rigid body rotations) don't induce interactions among unconnected     particles.*/  ScoreStatesTemp ss = get_required_score_states(rs);  for (ScoreStatesTemp::const_iterator it = ss.begin(); it != ss.end(); ++it) {    ModelObjectsTemps interactions = (*it)->get_interactions();    for (unsigned int i = 0; i < interactions.size(); ++i) {      add_edges(ps, interactions[i], map, *it, ret);    }  }  IMP_INTERNAL_CHECK(boost::num_vertices(ret) == ps.size(),                     "Wrong number of vertices " << boost::num_vertices(ret)                                                 << " vs " << ps.size());  return ret;}

/*  Calculate the line segment PaPb that is the shortest route between  two lines P1P2 and P3P4. Calculate also the values of mua and mub where  Pa = P1 + mua (P2 - P1)  Pb = P3 + mub (P4 - P3)  Return FALSE if no solution exists.*/Segment3D get_shortest_segment(const Segment3D &sa,                               const Segment3D &sb) {  const double eps= .0001;  algebra::Vector3D va= sa.get_point(1) - sa.get_point(0);  algebra::Vector3D vb= sb.get_point(1) - sb.get_point(0);  double ma= va*va;  double mb= vb*vb;  // if one of them is too short to have a well defined direction  // just look at an endpoint  if (ma < eps && mb < eps) {    return Segment3D(sa.get_point(0), sb.get_point(0));  } else if (ma < eps) {    return get_reversed(get_shortest_segment(sb, sa.get_point(0)));  } else if (mb < eps) {    return get_shortest_segment(sa, sb.get_point(0));  }  algebra::Vector3D vfirst = sa.get_point(0)- sb.get_point(0);  IMP_LOG_VERBOSE( vfirst << " | " << va << " | " << vb << std::endl);  double dab = vb*va;  double denom = ma * mb - dab * dab;  // they are parallel or anti-parallel  if (std::abs(denom) < eps) {    return get_shortest_segment_parallel(sa, sb);  }  double dfb = vfirst*vb;  double dfa = vfirst*va;  double numer = dfb * dab - dfa * mb;  double fa = numer / denom;  double fb = (dfb + dab * fa) / mb;  /*    denom = d2121 * d4343 - d4321 * d4321;    numer = d1343 * d4321 - d1321 * d4343;    *mua = numer / denom;    *mub = (d1343 + d4321 * (*mua)) / d4343;    pa->x = p1.x + *mua * p21.x;    pa->y = p1.y + *mua * p21.y;    pa->z = p1.z + *mua * p21.z;    pb->x = p3.x + *mub * p43.x;    pb->y = p3.y + *mub * p43.y;    pb->z = p3.z + *mub * p43.z;  */  algebra::Vector3D ra=get_clipped_point(sa, fa);  algebra::Vector3D rb=get_clipped_point(sb, fb);  IMP_LOG_VERBOSE( fa << " " << fb << std::endl);  return Segment3D(ra, rb);}

ResultAlign2D get_complete_alignment_no_preprocessing(    const cv::Mat &input, const cv::Mat &INPUT, const cv::Mat &POLAR1,    cv::Mat &m_to_align, const cv::Mat &POLAR2, bool apply) {  IMP_LOG_TERSE("starting complete 2D alignment with no preprocessing"                << std::endl);  cv::Mat aux1, aux2, aux3, aux4;  // auxiliary matrices  cv::Mat AUX1, AUX2, AUX3;        // ffts  algebra::Transformation2D transformation1, transformation2;  double angle1 = 0, angle2 = 0;  ResultAlign2D RA = get_rotational_alignment_no_preprocessing(POLAR1, POLAR2);  angle1 = RA.first.get_rotation().get_angle();  get_transformed(m_to_align, aux1, RA.first);  // rotate  get_fft_using_optimal_size(aux1, AUX1);  RA = get_translational_alignment_no_preprocessing(INPUT, AUX1);  algebra::Vector2D shift1 = RA.first.get_translation();  transformation1.set_rotation(angle1);  transformation1.set_translation(shift1);  get_transformed(m_to_align, aux2, transformation1);  // rotate  double ccc1 = get_cross_correlation_coefficient(input, aux2);  // Check the opposed angle  if (angle1 < PI) {    angle2 = angle1 + PI;  } else {    angle2 = angle1 - PI;  }  algebra::Rotation2D R2(angle2);  algebra::Transformation2D tr(R2);  get_transformed(m_to_align, aux3, tr);  // rotate  get_fft_using_optimal_size(aux3, AUX3);  RA = get_translational_alignment_no_preprocessing(INPUT, AUX3);  algebra::Vector2D shift2 = RA.first.get_translation();  transformation2.set_rotation(angle2);  transformation2.set_translation(shift2);  get_transformed(m_to_align, aux3, transformation2);  double ccc2 = get_cross_correlation_coefficient(input, aux3);  if (ccc2 > ccc1) {    if (apply) {      aux3.copyTo(m_to_align);    }    IMP_LOG_VERBOSE(" Align2D complete Transformation= "                    << transformation2 << " cross_correlation = " << ccc2                    << std::endl);    return ResultAlign2D(transformation2, ccc2);  } else {    if (apply) {      aux3.copyTo(m_to_align);    }    IMP_LOG_VERBOSE(" Align2D complete Transformation= "                    << transformation1 << " cross_correlation = " << ccc1                    << std::endl);    return ResultAlign2D(transformation1, ccc1);  }}

RigidBodyMover::RigidBodyMover(RigidBody d,                               Float max_translation, Float max_angle):  MonteCarloMover(d->get_model(), d->get_name()+" mover"){  IMP_LOG_VERBOSE("start RigidBodyMover constructor");  max_translation_=max_translation;  max_angle_ =max_angle;  pi_ = d.get_particle_index();  IMP_LOG_VERBOSE("finish mover construction" << std::endl);}

void Optimizer::set_is_optimizing_states(bool tf) const {  IMP_LOG_VERBOSE("Reseting OptimizerStates " << std::flush);  for (OptimizerStateConstIterator it = optimizer_states_begin();       it != optimizer_states_end(); ++it) {    IMP_CHECK_OBJECT(*it);    (*it)->set_is_optimizing(tf);    IMP_LOG_VERBOSE("." << std::flush);  }  IMP_LOG_VERBOSE("done." << std::endl);}

IMPEM2D_BEGIN_NAMESPACEResultAlign2D get_complete_alignment(const cv::Mat &input, cv::Mat &m_to_align,                                     bool apply) {  IMP_LOG_TERSE("starting complete 2D alignment " << std::endl);  cv::Mat autoc1, autoc2, aux1, aux2, aux3;  algebra::Transformation2D transformation1, transformation2;  ResultAlign2D RA;  get_autocorrelation2d(input, autoc1);  get_autocorrelation2d(m_to_align, autoc2);  RA = get_rotational_alignment(autoc1, autoc2, false);  double angle1 = RA.first.get_rotation().get_angle();  get_transformed(m_to_align, aux1, RA.first);  // rotate  RA = get_translational_alignment(input, aux1);  algebra::Vector2D shift1 = RA.first.get_translation();  transformation1.set_rotation(angle1);  transformation1.set_translation(shift1);  get_transformed(m_to_align, aux2, transformation1);  double ccc1 = get_cross_correlation_coefficient(input, aux2);  // Check for both angles that can be the solution  double angle2;  if (angle1 < PI) {    angle2 = angle1 + PI;  } else {    angle2 = angle1 - PI;  }  // rotate  algebra::Rotation2D R2(angle2);  algebra::Transformation2D tr(R2);  get_transformed(m_to_align, aux3, tr);  RA = get_translational_alignment(input, aux3);  algebra::Vector2D shift2 = RA.first.get_translation();  transformation2.set_rotation(angle2);  transformation2.set_translation(shift2);  get_transformed(m_to_align, aux3, transformation2);  double ccc2 = get_cross_correlation_coefficient(input, aux3);  if (ccc2 > ccc1) {    if (apply) {      aux3.copyTo(m_to_align);    }    IMP_LOG_VERBOSE(" Transformation= " << transformation2                                        << " cross_correlation = " << ccc2                                        << std::endl);    return em2d::ResultAlign2D(transformation2, ccc2);  } else {    if (apply) {      aux2.copyTo(m_to_align);    }    IMP_LOG_VERBOSE(" Transformation= " << transformation1                                        << " cross_correlation = " << ccc1                                        << std::endl);    return em2d::ResultAlign2D(transformation1, ccc1);  }}

IMPPMI_BEGIN_NAMESPACETransformMover::TransformMover(Model *m,                               Float max_translation, Float max_angle)    : MonteCarloMover(m, "Transform mover") {  IMP_LOG_VERBOSE("start TransformMover constructor");  max_translation_ = max_translation;  max_angle_ = max_angle;  constr_=0;  IMP_LOG_VERBOSE("finish mover construction" << std::endl);}

em::FittingSolutions pca_based_rigid_fitting(    ParticlesTemp ps, em::DensityMap *em_map, Float threshold, FloatKey,    algebra::PrincipalComponentAnalysis dens_pca_input) {  // find the pca of the density  algebra::PrincipalComponentAnalysis dens_pca;  if (dens_pca_input != algebra::PrincipalComponentAnalysis()) {    dens_pca = dens_pca_input;  } else {    algebra::Vector3Ds dens_vecs = em::density2vectors(em_map, threshold);    dens_pca = algebra::get_principal_components(dens_vecs);  }  // move the rigid body to the center of the map  core::XYZs ps_xyz = core::XYZs(ps);  algebra::Transformation3D move2center_trans = algebra::Transformation3D(      algebra::get_identity_rotation_3d(),      dens_pca.get_centroid() - core::get_centroid(core::XYZs(ps_xyz)));  for (unsigned int i = 0; i < ps_xyz.size(); i++) {    ps_xyz[i].set_coordinates(        move2center_trans.get_transformed(ps_xyz[i].get_coordinates()));  }  // find the pca of the protein  algebra::Vector3Ds ps_vecs;  for (core::XYZs::iterator it = ps_xyz.begin(); it != ps_xyz.end(); it++) {    ps_vecs.push_back(it->get_coordinates());  }  algebra::PrincipalComponentAnalysis ps_pca =      algebra::get_principal_components(ps_vecs);  IMP_IF_LOG(VERBOSE) {    IMP_LOG_VERBOSE("in pca_based_rigid_fitting, density PCA:" << std::endl);    IMP_LOG_WRITE(VERBOSE, dens_pca.show());    IMP_LOG_VERBOSE("particles PCA:" << std::endl);    IMP_LOG_WRITE(VERBOSE, ps_pca.show());  }  algebra::Transformation3Ds all_trans =      algebra::get_alignments_from_first_to_second(ps_pca, dens_pca);  em::FittingSolutions fs =      em::compute_fitting_scores(ps, em_map, all_trans, false);  fs.sort();  // compose the center translation to the results  em::FittingSolutions returned_fits;  for (int i = 0; i < fs.get_number_of_solutions(); i++) {    returned_fits.add_solution(        algebra::compose(fs.get_transformation(i), move2center_trans),        fs.get_score(i));  }  // move protein to the center of the map  algebra::Transformation3D move2center_inv = move2center_trans.get_inverse();  for (unsigned int i = 0; i < ps_xyz.size(); i++) {    ps_xyz[i].set_coordinates(        move2center_inv.get_transformed(ps_xyz[i].get_coordinates()));  }  return returned_fits;}

void Optimizer::update_states() const{  IMP_LOG_VERBOSE(          "Updating OptimizerStates " << std::flush);  for (OptimizerStateConstIterator it = optimizer_states_begin();       it != optimizer_states_end(); ++it) {    IMP_CHECK_OBJECT(*it);    (*it)->update();    IMP_LOG_VERBOSE( "." << std::flush);  }  IMP_LOG_VERBOSE( "done." << std::endl);}

IMPATOM_BEGIN_INTERNAL_NAMESPACEParticle* atom_particle(Model* m, const std::string& atom_type,                        Element element, bool is_hetatm,                        int atom_index, int residue_index,                        double x, double y, double z,                        double occupancy, double temp_factor) {  AtomType atom_name;  std::string string_name = atom_type;  // determine AtomType  if (is_hetatm) {    string_name = "HET:" + string_name;    if (!get_atom_type_exists(string_name)) {      atom_name = add_atom_type(string_name, element);    } else {      atom_name = AtomType(string_name);    }  } else {  // ATOM line    boost::trim(string_name);    if (string_name.empty()) {      string_name = "UNK";    }    if (!AtomType::get_key_exists(string_name)) {      IMP_LOG_VERBOSE("ATOM record type not found: /""                      << string_name << "/" in PDB file " << std::endl);      atom_name = add_atom_type(string_name, element);    } else {      atom_name = AtomType(string_name);    }  }  // new particle  Particle* p = new Particle(m);  p->add_attribute(get_pdb_index_key(), atom_index);  algebra::Vector3D v(x, y, z);  // atom decorator  Atom d = Atom::setup_particle(p, atom_name);  std::ostringstream oss;  oss << "Atom " + atom_name.get_string() << " of residue " << residue_index;  p->set_name(oss.str());  core::XYZ::setup_particle(p, v).set_coordinates_are_optimized(true);  d.set_input_index(atom_index);  d.set_occupancy(occupancy);  d.set_temperature_factor(temp_factor);  d.set_element(element);  // check if the element matches  Element e2 = get_element_for_atom_type(atom_name);  if (element != e2) {    IMP_LOG_VERBOSE(        "AtomType element and PDB line elements don't match. AtomType "        << e2 << " vs. determined from PDB " << element << std::endl);  }  return p;}

RigidBodyMover::RigidBodyMover(RigidBody d, Float max_translation,                               Float max_angle)    : MonteCarloMover(d->get_model(), d->get_name() + " mover") {  IMP_USAGE_CHECK(      d.get_coordinates_are_optimized(),      "Rigid body passed to RigidBodyMover"          << " must be set to be optimized. particle: " << d->get_name());  IMP_LOG_VERBOSE("start RigidBodyMover constructor");  max_translation_ = max_translation;  max_angle_ = max_angle;  pi_ = d.get_particle_index();  IMP_LOG_VERBOSE("finish mover construction" << std::endl);}

TransformMover::TransformMover(Model *m, IMP::ParticleIndexAdaptor p1i, IMP::ParticleIndexAdaptor p2i,                               Float max_translation, Float max_angle)    : MonteCarloMover(m, "Transform mover") {  IMP_LOG_VERBOSE("start TransformMover constructor");  //this constructor defines a rotation about an axis defined by two particles  p1i_ = p1i;  p2i_ = p2i;  max_translation_ = max_translation;  max_angle_ = max_angle;  constr_=2;  called_=0;  not_accepted_=0;  IMP_LOG_VERBOSE("finish mover construction" << std::endl);}

void RestraintCache::validate() const {#if IMP_HAS_CHECKS >= IMP_INTERNAL  IMP_OBJECT_LOG;  IMP_LOG_VERBOSE("Validating cache...." << std::endl);  for (Cache::ContentIterator it = cache_.contents_begin();       it != cache_.contents_end(); ++it) {    double score = it->value;    double new_score = cache_.get_generator()(it->key, cache_);    IMP_LOG_VERBOSE("Validating " << it->key << std::endl);    IMP_INTERNAL_CHECK_FLOAT_EQUAL(score, new_score,                                   "Cached and computed scores don't match "                                       << score << " vs " << new_score);  }#endif}

void Fine2DRegistrationRestraint::set_subject_image(em2d::Image *subject) {  // Read the registration parameters from the subject images  algebra::Vector3D euler = subject->get_header().get_euler_angles();  algebra::Rotation3D R =      algebra::get_rotation_from_fixed_zyz(euler[0], euler[1], euler[2]);  algebra::Vector3D origin = subject->get_header().get_origin();  algebra::Vector3D translation(origin[0] * params_.pixel_size,                                origin[1] * params_.pixel_size, 0.0);  subject_->set_data(subject->get_data());  // deep copy, avoids leaks  unsigned int rows = subject_->get_header().get_number_of_rows();  unsigned int cols = subject_->get_header().get_number_of_columns();  if (projection_->get_header().get_number_of_columns() != cols ||      projection_->get_header().get_number_of_rows() != rows) {    projection_->set_size(rows, cols);  }  PP_.set_rotation(R);  PP_.set_translation(translation);  double s = params_.pixel_size;  algebra::Vector3D min_values(-s * rows, -s * cols, 0.0);  algebra::Vector3D max_values(s * rows, s * cols, 0.0);  PP_.set_proper_ranges_for_keys(this->get_model(), min_values, max_values);  IMP_LOG_VERBOSE("Subject set for Fine2DRegistrationRestraint" << std::endl);}

double WeightedExcludedVolumeRestraint::unprotected_evaluate(    DerivativeAccumulator *accum) const{    bool calc_deriv = accum? true: false;    IMP_LOG_VERBOSE("before resample/n");    // //generate the transformed maps    // std::vector<DensityMap*> transformed_maps;    // for(int rb_i=0;rb_i<rbs_.size();rb_i++){    //   DensityMap *dm=create_density_map(    //       atom::get_bounding_box(atom::Hierarchy(rbs_[rb_i])),spacing);    //   get_transformed_into(    //       rbs_surface_maps_[rb_i],    //       rbs_[rb_i].get_transformation()*rbs_orig_trans_[rb_i],    //       dm,    //       false);    //   transformed_maps.push_back(dm);    // }    double score=0.;    // MRCReaderWriter mrw;    // for(int i=0;i<transformed_maps.size();i++){    //   std::stringstream ss;    //   ss<<"transformed_map_"<<i<<".mrc";    //   std::stringstream s1;    //   s1<<"transformed_pdb_"<<i<<".pdb";    //   atom::write_pdb(atom::Hierarchy(rbs_[i]),s1.str().c_str());    //   write_map(transformed_maps[i],ss.str().c_str(),mrw);    //   for(int j=i+1;j<transformed_maps.size();j++){    //     if (get_interiors_intersect(transformed_maps[i],    //                                 transformed_maps[j])){    //     score += CoarseCC::cross_correlation_coefficient(    //                              *transformed_maps[i],    //                              *transformed_maps[j],1.,false,true);    //     }    //   }    // }    em::SurfaceShellDensityMaps resampled_surfaces;    for(unsigned int i=0; i<rbs_.size(); i++) {        kernel::ParticlesTemp rb_ps=rb_refiner_->get_refined(rbs_[i]);        resampled_surfaces.push_back(new em::SurfaceShellDensityMap(rb_ps,1.));    }    for(unsigned int i=0; i<rbs_.size(); i++) {        for(unsigned int j=i+1; j<rbs_.size(); j++) {            if (get_interiors_intersect(resampled_surfaces[i],                                        resampled_surfaces[j])) {                score += em::CoarseCC::cross_correlation_coefficient(                             resampled_surfaces[i],                             resampled_surfaces[j],1.,true,FloatPair(0.,0.));            }        }    }    if (calc_deriv) {        IMP_WARN("WeightedExcludedVolumeRestraint currently"                 <<" does not support derivatives/n");    }    /*for(int i=resampled_surfaces.size()-1;i<0;i--){      delete(resampled_surfaces[i]);      }*/    return score;}

MonteCarloMoverResult MonteCarlo::do_move() {    ParticleIndexes ret;    double prob = 1.0;    for (MoverIterator it = movers_begin(); it != movers_end(); ++it) {        IMP_LOG_VERBOSE("Moving using " << (*it)->get_name() << std::endl);        IMP_CHECK_OBJECT(*it);        {            // IMP_LOG_CONTEXT("Mover " << (*it)->get_name());            MonteCarloMoverResult cur = (*it)->propose();            ret += cur.get_moved_particles();            prob *= cur.get_proposal_ratio();        }        IMP_LOG_VERBOSE("end/n");    }    return MonteCarloMoverResult(ret, prob);}

em2d::Images get_projections(const ParticlesTemp &ps,                             const RegistrationResults &registration_values,                             int rows, int cols,                             const ProjectingOptions &options, Strings names) {  IMP_LOG_VERBOSE("Generating projections from registration results"                  << std::endl);  if (options.save_images && (names.size() < registration_values.size())) {    IMP_THROW("get_projections: Insufficient number of image names provided",              IOException);  }  unsigned long n_projs = registration_values.size();  em2d::Images projections(n_projs);  // Precomputation of all the possible projection masks for the particles  MasksManagerPtr masks(      new MasksManager(options.resolution, options.pixel_size));  masks->create_masks(ps);  for (unsigned long i = 0; i < n_projs; ++i) {    IMP_NEW(em2d::Image, img, ());    img->set_size(rows, cols);    img->set_was_used(true);    String name = "";    if (options.save_images) name = names[i];    get_projection(img, ps, registration_values[i], options, masks, name);    projections[i] = img;  }  return projections;}

double Fine2DRegistrationRestraint::unprotected_evaluate(    DerivativeAccumulator *accum) const {  IMP_UNUSED(accum);  calls_++;  IMP_USAGE_CHECK(accum == nullptr,                  "Fine2DRegistrationRestraint: This restraint does not "                  "provide derivatives ");  // projection_ needs to be mutable, son this const function can change it.  // project_particles changes the matrix of projection_  ProjectingOptions options(params_.pixel_size, params_.resolution);  double score = 0;  try {    do_project_particles(ps_, projection_->get_data(), PP_.get_rotation(),                         PP_.get_translation(), options, masks_);    score = score_function_->get_score(subject_, projection_);  }  catch (cv::Exception &e) {    IMP_LOG(WARNING,            "Fine2DRegistration. Error computing the score: "            "Returning 1 (maximum score). Most probably because projecting "            "out of the image size."                << e.what() << std::endl);    score = 1.0;  }  IMP_LOG_VERBOSE("Fine2DRegistration. Score: " << score << std::endl);  return score;}