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

std::string create_temporary_file_name(std::string prefix, std::string suffix) {
  char *env = getenv("IMP_BUILD_ROOT");
  std::string imp_tmp;
  if (env) {
    imp_tmp = internal::get_concatenated_path(env, "build/tmp");
  }
#if defined _MSC_VER
  std::string tpathstr;
  if (imp_tmp.empty()) {
    TCHAR tpath[MAX_PATH];
    DWORD dwRetVal = GetTempPath(MAX_PATH, tpath);
    if (dwRetVal > MAX_PATH || (dwRetVal == 0)) {
      IMP_THROW("Unable to find temporary path", IOException);
    }
    tpathstr = tpath;
  } else {
    tpathstr = imp_tmp;
  }
  char filename[MAX_PATH];
  if (GetTempFileName(tpathstr.c_str(), prefix.c_str(), 0, filename) == 0) {
    IMP_THROW("Unable to create temp file in " << tpathstr, IOException);
  }
  return std::string(filename) + suffix;
#else
  std::string pathprefix;
  if (imp_tmp.empty()) {
    pathprefix = "/tmp";
  } else {
    pathprefix = imp_tmp;
  }
  std::string templ =
      internal::get_concatenated_path(pathprefix, prefix + ".XXXXXX");
  boost::scoped_array<char> filename;
  filename.reset(new char[templ.size() + suffix.size() + 1]);
  std::copy(templ.begin(), templ.end(), filename.get());
#ifdef __APPLE__
  std::copy(suffix.begin(), suffix.end(), filename.get() + templ.size());
  filename[templ.size() + suffix.size()] = '\0';
  int fd = mkstemps(filename.get(), suffix.size());
  if (fd == -1) {
    IMP_THROW("Unable to create temporary file: " << filename.get(),
              IOException);
  }
  close(fd);
#else
  filename[templ.size()] = '\0';
  int fd = mkstemp(filename.get());
  if (fd == -1) {
    IMP_THROW("Unable to create temporary file: " << filename.get(),
              IOException);
  }
  close(fd);
  std::copy(suffix.begin(), suffix.end(), filename.get() + templ.size());
  filename[templ.size() + suffix.size()] = '\0';
#endif
  return std::string(filename.get());
#endif
}

bool MSConnectivityScore::check_assignment(NNGraph &G, unsigned int node_handle,
                                           Assignment const &assignment,
                                           EdgeSet &picked) const {
  MSConnectivityRestraint::ExperimentalTree::Node const *node =
      tree_.get_node(node_handle);
  MSConnectivityRestraint::ExperimentalTree::Node::Label const &lb =
      node->get_label();
  Vector<Tuples> new_tuples;
  Ints empty_vector;
  for (unsigned int i = 0; i < lb.size(); ++i) {
    int prot_count = lb[i].second;
    unsigned int id = lb[i].first;
    while (new_tuples.size() < id)
      new_tuples.push_back(Tuples(empty_vector, 0));
    if (prot_count > 0) {
      if (!assignment[id].empty()) {
        Ints const &configuration = assignment[id].get_tuple();
        if (prot_count > int(configuration.size())) {
          IMP_THROW("Experimental tree is inconsistent",
                    IMP::ValueException);
        }
        new_tuples.push_back(Tuples(configuration, prot_count));
      } else {
        IMP_THROW("Experimental tree is inconsistent",
                  IMP::ValueException);
      }
    } else
      new_tuples.push_back(Tuples(empty_vector, 0));
  }
  while (new_tuples.size() <
         restraint_.particle_matrix_.get_number_of_classes())
    new_tuples.push_back(Tuples(empty_vector, 0));
  Assignment new_assignment(new_tuples);
  if (new_assignment.empty()) return false;
  do {
    NNGraph ng = build_subgraph_from_assignment(G, new_assignment);
    if (is_connected(ng)) {
      EdgeSet n_picked;
      bool good = true;
      for (unsigned int i = 0; i < node->get_number_of_children(); ++i) {
        unsigned int child_handle = node->get_child(i);
        if (!check_assignment(ng, child_handle, new_assignment, n_picked)) {
          good = false;
          break;
        }
      }
      if (good) {
        add_edges_to_set(ng, n_picked);
        picked.insert(n_picked.begin(), n_picked.end());
        return true;
      }
    }
  } while (new_assignment.next());
  return false;
}

void
KinematicForest::add_edge(Joint* joint)
{
  joint->set_owner_kf( this );
  IMP::core::RigidBody parent_rb = joint->get_parent_node();
  IMP::core::RigidBody child_rb = joint->get_child_node();
  KinematicNode parent_kn, child_kn;

  // decorate parent and store here
  kernel::Particle* parent_p = parent_rb.get_particle();
  if(!KinematicNode::get_is_setup( parent_p ) ) {
    parent_kn = KinematicNode::setup_particle( parent_p, this );
    nodes_.insert( parent_kn );
    roots_.insert( parent_kn );
  } else {
    parent_kn = KinematicNode( parent_p );
    if( parent_kn.get_owner() != this ) {
      IMP_THROW( "the parent rigid body " << parent_rb
                 << " in the joint " << joint
                 << " was already stored in a different kinematic forest -"
                 << " this IMP version does not support such switching",
                 IMP::base::ValueException );
    }
  }

  // decorate child and store here
  kernel::Particle* child_p = child_rb.get_particle();
  if(!KinematicNode::get_is_setup( child_p ) ) {
    child_kn = KinematicNode::setup_particle( child_p, this, joint );
    nodes_.insert( child_kn );
  } else {
    child_kn = KinematicNode( child_p );
    if( child_kn.get_owner() != this ){
      IMP_THROW( "the child rigid body " << child_rb
                 << " in the joint " << joint
                 << " was already stored in a different kinematic forest -"
                 << " this IMP version does not support such switching",
                 IMP::base::ValueException );
    }

    if( roots_.find( child_kn) != roots_.end() ) {
      roots_.erase( child_kn ); // will no longer be a root
    } else {
      IMP_THROW( "IMP currently does not support switching of "
                 << " parents in a kinematic tree",
                 IMP::base::ValueException );
    }
  }

  // store joint
  parent_kn.add_out_joint( joint );
  child_kn.set_in_joint( joint );
  joints_.push_back( joint);
}

  void MultivariateFNormalSufficient::set_W(const MatrixXd& W)
{
    if (W.rows() != W_.rows() || W.cols() != W_.cols() || W != W_)
    {
        if (W.cols() != W.rows()) {
            IMP_THROW("need a square matrix!", ModelException);
            }
        if (W.rows() != M_) {
            IMP_THROW("size mismatch for W!" , ModelException);
            }
        W_ = W;
        IMP_LOG(TERSE, "MVN:   set W to new matrix"<< std::endl);
        flag_PW_ = false;
    }
    flag_W_ = true;
}

 MatrixXd MultivariateFNormalSufficient::evaluate_second_derivative_FM_FM()
   const
 {
     if (N_!=1) IMP_THROW("not implemented when N>1", ModelException);
     MatrixXd ret(get_P()/IMP::square(factor_));
     return ret;
 }

IMPINTEGRATIVEDOCKING_BEGIN_INTERNAL_NAMESPACE

void ResidueContent::read_content_file(const std::string& file_name) {
  std::ifstream in_file(file_name.c_str());
  if (!in_file) {
    IMP_THROW("Can't open file " << file_name, IMP::IOException);
  }

  std::string line;
  while (!in_file.eof()) {
    getline(in_file, line);
    boost::trim(line);  // remove all spaces
    // skip comments
    if (line[0] == '#' || line[0] == '\0') continue;
    std::vector<std::string> split_results;
    boost::split(split_results, line, boost::is_any_of("\t "),
                 boost::token_compress_on);
    if (split_results.size() != 2) continue;
    int counter = atoi(split_results[1].c_str());
    IMP::atom::ResidueType residue_type =
        IMP::atom::ResidueType(split_results[0]);
    residue_content_[residue_type] = counter;
  }
  in_file.close();
}

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

}

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

ResidueType get_residue_type(char c) {
  if (rp_map.find(c) == rp_map.end()) {
    IMP_THROW("Residue name not found " << c, ValueException);
  } else {
    return rp_map.find(c)->second;
  }
}

double get_volume_from_residue_type(ResidueType rt) {
  typedef std::pair<ResidueType, double> RP;
  static const RP radii[]={RP(ResidueType("ALA"), 2.516),
                     RP(ResidueType("ARG"), 3.244),
                     RP(ResidueType("ASN"), 2.887),
                     RP(ResidueType("ASP"), 2.866),
                     RP(ResidueType("CYS"), 2.710),
                     RP(ResidueType("GLN"), 3.008),
                     RP(ResidueType("GLU"), 2.997),
                     RP(ResidueType("GLY"), 2.273),
                     RP(ResidueType("HIS"), 3.051),
                     RP(ResidueType("ILE"), 3.047),
                     RP(ResidueType("LEU"), 3.052),
                     RP(ResidueType("LYS"), 3.047),
                     RP(ResidueType("MET"), 3.068),
                     RP(ResidueType("PHE"), 3.259),
                     RP(ResidueType("PRO"), 2.780),
                     RP(ResidueType("SER"), 2.609),
                     RP(ResidueType("THR"), 2.799),
                     RP(ResidueType("TRP"), 3.456),
                     RP(ResidueType("TYR"), 3.318),
                     RP(ResidueType("VAL"), 2.888)};
  static const IMP::base::map<ResidueType, double> radii_map(radii,
                                                                 radii
                                                                 +sizeof(radii)
                                                                 /sizeof(RP));
  if (radii_map.find(rt) == radii_map.end()) {
    IMP_THROW( "Can't approximate volume of non-standard residue "
               << rt, ValueException);
  }
  double r= radii_map.find(rt)->second;
  IMP_INTERNAL_CHECK(r>0, "Read garbage r for "<< rt);
  return algebra::get_volume(algebra::Sphere3D(algebra::get_zero_vector_d<3>(),
                                               r));
}

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 GaussianProcessInterpolationRestraint::get_logdet_hessian() const {
  // compute ldlt
  IMP_Eigen::LDLT<IMP_Eigen::MatrixXd, IMP_Eigen::Upper> ldlt(get_hessian());
  if (!ldlt.isPositive())
    IMP_THROW("Hessian matrix is not positive definite!", ModelException);
  return ldlt.vectorD().array().abs().log().sum();
}

NNGraph MSConnectivityScore::find_threshold() const {
  double max_dist = 1.1 * restraint_.particle_matrix_.max_distance(),
         min_dist = restraint_.particle_matrix_.min_distance();
  NNGraph g = create_nn_graph(min_dist);
  {
    std::set<std::pair<unsigned int, unsigned int> > picked;
    if (perform_search(g, picked)) {
      return pick_graph(picked);
    }
    g = create_nn_graph(max_dist);
    if (!perform_search(g, picked)) {
      IMP_THROW("Cannot build a nearest neighbor graph",
                IMP::ValueException);
    }
  }
  EdgeSet picked;
  while (max_dist - min_dist > eps_) {
    double dist = 0.5 * (max_dist + min_dist);
    g = create_nn_graph(dist);
    EdgeSet tmp_picked;
    if (perform_search(g, tmp_picked)) {
      picked.swap(tmp_picked);
      max_dist = dist;
    } else
      min_dist = dist;
  }
  return pick_graph(picked);
}

double GaussianProcessInterpolationRestraintSparse::unprotected_evaluate(
                         DerivativeAccumulator *accum) const
{
    //check if the functions have changed
    const_cast<GaussianProcessInterpolationRestraintSparse*>(this)->
        update_mean_and_covariance();

    double ene = mvn_->evaluate();

    if (accum)
    {
        cholmod_dense *dmv = mvn_->evaluate_derivative_FM();
        if (dmv->xtype != CHOLMOD_REAL)
            IMP_THROW("matrix type is not real, update code here first",
                    ModelException);
        double *dmvx=(double*) dmv->x;
        //derivatives for mean particles
        for (size_t i=0; i<M_; i++)
        {
            DerivativeAccumulator a(*accum, dmvx[i]);
            gpi_->mean_function_->add_to_derivatives(gpi_->x_[i], a);
        }
        cholmod_free_dense(&dmv, c_);

        //derivatives for covariance particles
        cholmod_sparse *tmp(mvn_->evaluate_derivative_Sigma());
        cholmod_triplet *dmvS = cholmod_sparse_to_triplet(tmp, c_);
        cholmod_free_sparse(&tmp, c_);
        if ((dmvS->itype != CHOLMOD_INT) && (dmvS->xtype != CHOLMOD_REAL))
            IMP_THROW("matrix type is not real or coefficients are not int!",
                    ModelException);
        int *dmvi=(int*) dmvS->i;
        int *dmvj=(int*) dmvS->j;
        dmvx=(double*) dmvS->x;
        for (size_t p=0; p<dmvS->nzmax; ++p)
        {
            int i=dmvi[p];
            int j=dmvj[p];
            double val=dmvx[p];
            DerivativeAccumulator a(*accum, val);
            gpi_->covariance_function_->add_to_derivatives(
                    gpi_->x_[i],gpi_->x_[j], a);
        }
        cholmod_free_triplet(&dmvS,c_);
    }
    return ene;
}

void MSConnectivityRestraint::ExperimentalTree::connect(unsigned int parent,
                                                        unsigned int child) {
  if (finalized_)
    IMP_THROW("Cannot add new edges to finalized tree",
              IMP::ValueException);
  nodes_[parent].children_.push_back(child);
  nodes_[child].parents_.push_back(parent);
}

 void MultivariateFNormalSufficientSparse::set_FX(const MatrixXd& FX,
         double cutoff)
 {
   if (FX.rows() != FX_.rows() || FX.cols() != FX_.cols() || FX != FX_){
       if (FX.rows() != N_) {
           IMP_THROW("size mismatch for FX in the number of repetitions: got "
                   << FX.rows() << " instead of "<<N_, ModelException);
           }
       if (FX.cols() != M_) {
           IMP_THROW("size mismatch for FX in the number of variables: got "
                   <<FX.cols() << " instead of "<<M_, ModelException);
           }
       FX_=FX;
       IMP_LOG(TERSE, "MVNsparse:   set FX to new matrix"<< std::endl);
       compute_sufficient_statistics(cutoff);
   }
 }

double HybridMonteCarlo::do_evaluate(const kernel::ParticleIndexes &) const
{
    if (get_use_incremental_scoring_function())
        IMP_THROW("Incremental scoring not supported", ModelException);
    double ekin = md_->get_kinetic_energy();
    double epot = get_scoring_function()->evaluate(false);
    return ekin + epot;
}

VectorXd MultivariateFNormalSufficient::get_Sigma_eigenvalues() const
{
    Eigen::SelfAdjointEigenSolver<MatrixXd> eigensolver(get_Sigma(),
            Eigen::EigenvaluesOnly);
    if (eigensolver.info() != Eigen::Success)
            IMP_THROW("Could not determine eigenvalues!", ModelException);
    return eigensolver.eigenvalues();
}