IMPEM_BEGIN_NAMESPACE

float 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 here
ProteomicsEMAlignmentAtomic::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_NAMESPACE

ResultAlign2D 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_NAMESPACE

TransformMover::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_NAMESPACE

Particle* 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);
}