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_NAMESPACE

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

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

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

    }
  }
}

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

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

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

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

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

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

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

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

}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

IMPISD_BEGIN_NAMESPACE

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

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