Index: kernel/include/IMP/restraints/ConnectivityRestraint.h =================================================================== --- kernel/include/IMP/restraints/ConnectivityRestraint.h (revision 349) +++ kernel/include/IMP/restraints/ConnectivityRestraint.h (working copy) @@ -20,86 +20,43 @@ namespace IMP { +class PairScore; //! Connectivity restraint. -/** Restrict max distance between at least one pair of particles of any - two distinct types. +/** The restraint takes several sets (added using the add_set method) + and ensures that the sets of points are connected. More precisely, + the restraint scores based on a minimum spanning tree on the set + of particles in the union of the sets. Edges between two particles + from the same set have weight 0. Edges between two nodes from + different sets have weights given by the PairScore for that pair. - Calculate the distance restraints for the given particles. Use the - smallest restraints that will connect one particle of each type - together (i.e. a minimum spanning tree with nodes corresponding to - particle types and the edge weights corresponding to restraint - violation scores). + \ingroup restraint */ class IMPDLLEXPORT ConnectivityRestraint : public Restraint { public: - ConnectivityRestraint(Model* model, - std::vector& particle_indexes, - IntKey type, - BasicScoreFuncParams* score_func_params); + ConnectivityRestraint(PairScore* ps); - ConnectivityRestraint(Model* model, - std::vector& particle_indexes, - IntKey type, - FloatKey attr_name, - BasicScoreFuncParams* score_func_params); + //! Add a new set of particles. + void add_set(const Particles &ps); + //! clear all the sets and start over + void clear_sets(); + virtual ~ConnectivityRestraint(); IMP_RESTRAINT("0.5", "Daniel Russel") protected: - // switch to using rsr_scores to allow STL sorting - class RestraintScore - { - public: - RestraintScore() {} - ~RestraintScore() {} - void evaluate() { - score_ = rsr_->evaluate(false); - } - bool operator<(const RestraintScore& rs) const { - return score_ < rs.score_; - } - int part1_type_; - int part2_type_; - DistanceRestraint* rsr_; - Float score_; - }; + std::auto_ptr ps_; - //! restraints and their scores - std::list rsr_scores_; - - //! Set up initial values for the constructors. - /** \param[in] model Pointer to the model. - \param[in] particle_indexes Vector of particle indices. \todo Should use - Particle pointers instead. - \param[in] type The attribute used to determine if particles - are equivalent. + // the indices for the first particle in each set + /* + set_offset_[i] is the first index of set i and set_offset_[i+1] is + the first index not in i. It is always defined if i is valid. */ - void set_up(Model* model, std::vector& particle_indexes, - IntKey type); - - //! variables to determine the particle type - IntKey type_; - - //! number of particles in the restraint - int num_particles_; - - //! maximum type (type can be from 0 to max_type-1) - int max_type_; - //! number of particle types - int num_types_; - //! particle types - std::vector particle_type_; - - //! total number of restraints being tested - int num_restraints_; - - //! each unconnected tree has a non-zero id - std::vector tree_id_; + std::vector set_offsets_; }; } // namespace IMP Index: kernel/src/restraints/ConnectivityRestraint.cpp =================================================================== --- kernel/src/restraints/ConnectivityRestraint.cpp (revision 349) +++ kernel/src/restraints/ConnectivityRestraint.cpp (working copy) @@ -7,314 +7,135 @@ * Copyright 2007-8 Sali Lab. All rights reserved. * */ - -#include - #include "IMP/Model.h" #include "IMP/Particle.h" #include "IMP/log.h" #include "IMP/restraints/ConnectivityRestraint.h" -#include "../mystdexcept.h" +#include "IMP/PairScore.h" +#include +#include + namespace IMP { -//! Constructor using a given mean for particle-particle expected distance. -/** - \param[in] model Pointer to the model. - \param[in] particle_indexes Vector of particle indices. \todo Should use - Particle pointers instead. - \param[in] type The attribute used to determine if particles are equivalent. - \param[in] score_func_params Parameters for creating a score function. - */ -ConnectivityRestraint::ConnectivityRestraint(Model* model, - std::vector& particle_indexes, IntKey type, - BasicScoreFuncParams* score_func_params) -{ - std::list::iterator rs_iter; - set_up(model, particle_indexes, type); - - // set up the restraints - rs_iter = rsr_scores_.begin(); - for (int i = 0; i < num_particles_ - 1; i++) { - for (int j = i + 1; j < num_particles_; j++) { - if (particle_type_[j] != particle_type_[i]) { - if (rs_iter == rsr_scores_.end()) { - IMP_failure("Over ran the calculated number of restraints in " - "ConnectivityRestraint", - std::out_of_range("Over ran the calculated number of " - "restraints")); - } else { - rs_iter->part1_type_ = particle_type_[i]; - rs_iter->part2_type_ = particle_type_[j]; - - // create the restraint - UnaryFunction *sf = score_func_params->create_score_func(); - rs_iter->rsr_ = new DistanceRestraint(get_particle(i), - get_particle(j), - sf); - ++rs_iter; - } - } - } - } - - IMP_LOG(VERBOSE, - "Number of types: " << num_types_ << " max_type_: " << max_type_ << - " num_restraints_: " << num_restraints_ << " num_particles_: " << - num_particles_ << std::endl); +ConnectivityRestraint::ConnectivityRestraint(PairScore *ps): ps_(ps) +{ + clear_sets(); } -//! Constructor using a given attribute for particle-particle expected distance. -/** \param[in] model Pointer to the model. - \param[in] particle_indexes Vector of particle indices. \todo Should use - Particle pointers instead. - \param[in] type The attribute used to determine if particles are equivalent. - \param[in] attr_name Name to get radii to calculate the mean distance. - \param[in] score_func_params Parameters for creating a score function. - */ -ConnectivityRestraint::ConnectivityRestraint(Model* model, - std::vector& particle_indexes, IntKey type, - FloatKey attr_name, BasicScoreFuncParams* score_func_params) -{ - std::list::iterator rs_iter; - IMP_LOG(VERBOSE, "ConnectivityRestraint constructor"); - set_up(model, particle_indexes, type); - // set up the restraints - Float actual_mean; - rs_iter = rsr_scores_.begin(); - for (int i = 0; i < num_particles_ - 1; i++) { - for (int j = i + 1; j < num_particles_; j++) { - if (particle_type_[j] != particle_type_[i]) { - if (rs_iter == rsr_scores_.end()) { - IMP_failure("Over ran the calculated number of restraints in " - "ConnectivityRestraint", - std::out_of_range("Over ran the calculated number of " - "restraints")); - } else { - rs_iter->part1_type_ = particle_type_[i]; - rs_iter->part2_type_ = particle_type_[j]; - // Use those radii to calculate the expected distance - Float fi = get_particle(i)->get_value(attr_name); - Float fj = get_particle(j)->get_value(attr_name); - actual_mean = fi + fj; - - score_func_params->set_mean(actual_mean); - - // create the restraint - UnaryFunction *sf = score_func_params->create_score_func(); - rs_iter->rsr_ = new DistanceRestraint(get_particle(i), - get_particle(j), - sf); - ++rs_iter; - } - } - } - } - - IMP_LOG(VERBOSE, - "Number of types: " << num_types_ << " max_type_: " << max_type_ << - " num_restraints_: " << num_restraints_ << " num_particles_: " << - num_particles_ << std::endl); +ConnectivityRestraint::~ConnectivityRestraint() +{ } - -//! Set up initial values for the constructors. -/** \param[in] model Pointer to the model. - \param[in] particle_indexes Vector of particle indices. \todo Should use - Particle pointers instead. - \param[in] type The attribute used to determine if particles are equivalent. - */ -void ConnectivityRestraint::set_up(Model* model, - std::vector& particle_indexes, - IntKey type) -{ - Particle* p1; - - IMP_LOG(VERBOSE, "init ConnectivityRestraint"); - - //IMP_LOG(VERBOSE, "got model data " << model_data_); - - /** number of particles in the restraint */ - num_particles_ = particle_indexes.size(); - particle_type_.resize(num_particles_); - - IMP_LOG(VERBOSE, "set up particle types"); - type_ = type; - // set up the particles, their position indexes, and their type indexes - for (int i = 0; i < num_particles_; i++) { - p1 = model->get_particle(particle_indexes[i]); - add_particle(p1); - //type_.push_back(p1->get_attribute(type)); +struct Pair { + typedef Pair This; + Pair(){p_[0]=-1; p_[1]=-1;} + Pair(int a, int b){ + p_[0]= std::min(a,b); + p_[1]= std::max(a,b); } - - IMP_LOG(VERBOSE, "Size of particles: " << number_of_particles() - << " num_particles_: " - << num_particles_ << " particle_type_: " << particle_type_.size()); - IMP_LOG(VERBOSE, "Figure out number of types"); - // figure out how many types there are - int next_type; - max_type_ = 0; - num_types_ = 0; - for (int i = 0; i < num_particles_; i++) { - next_type = get_particle(i)->get_value(type_); - if (max_type_ < next_type) { - max_type_ = next_type; - } - - num_types_++; - particle_type_[i] = next_type; - // if this type already was used, decrement numbe of types - for (int j = 0; j < i; j++) { - if (particle_type_[j] == next_type) { - num_types_--; - break; - } - } + int operator[](unsigned int i) const { + return p_[i]; } - - // types can range from 0 to max_type_ - 1 - max_type_++; - IMP_LOG(VERBOSE, "Max types: " << max_type_); - - IMP_LOG(VERBOSE, "Figure out number of restraints"); - // figure out how many restraints there are - // Could use num_restraints = sum((S-si)si)) where is total number of - // particles and si is number of particles of type i (summed over all types) - num_restraints_ = num_particles_ * (num_particles_ - 1) / 2; - for (int i = 0; i < num_particles_ - 1; i++) { - for (int j = i + 1; j < num_particles_; j++) { - if (particle_type_[j] == particle_type_[i]) { - num_restraints_--; - } - } + bool is_default() const { + return p_[0] < 0 && p_[1] < 0; } + IMP_COMPARISONS_2(p_[0], p_[1]); + int p_[2]; +}; - IMP_LOG(VERBOSE, "Num restraints: " << num_restraints_); - // use number of restraints and number of types to set up some of the - // ... main arrays - tree_id_.resize(max_type_); - rsr_scores_.resize(num_restraints_); +void ConnectivityRestraint::add_set(const Particles &ps) { + IMP_check(!ps.empty(), "Cannot add empty set to ConnectivityRestraint", + InvalidStateException("Empty set")); + unsigned int sz= ps.size(); + Restraint::add_particles(ps); + set_offsets_.push_back(set_offsets_.back()+sz); } -//! Destructor -ConnectivityRestraint::~ConnectivityRestraint() -{ - std::list::iterator rs_iter; - - for (rs_iter = rsr_scores_.begin(); rs_iter != rsr_scores_.end(); ++rs_iter) { - delete(rs_iter->rsr_); - } +void ConnectivityRestraint::clear_sets() { + set_offsets_.clear(); + Restraint::clear_particles(); + set_offsets_.push_back(0); } - -//! Evaluate this restraint and return the score. -/** Calculate the distance restraints for the given particles. Use the smallest - restraints that will connect one particle of each type together (i.e. a - minimum spanning tree with nodes corresponding to particle types and the - edge weights corresponding to restraint violation scores). - - \param[in] accum If not NULL, use this object to accumulate partial first - derivatives. - \return score associated with this restraint for the given state of - the model. - */ Float ConnectivityRestraint::evaluate(DerivativeAccumulator *accum) { - IMP_LOG(VERBOSE, "evaluate ConnectivityRestraint"); - - std::list::iterator rs_iter; - - // calculate the scores for all of the restraints - IMP_LOG(VERBOSE, "calculate restraint scores"); - for (rs_iter = rsr_scores_.begin(); rs_iter != rsr_scores_.end(); ++rs_iter) { - rs_iter->evaluate(); - } - - // sort by the scores - rsr_scores_.sort(); - - // calculate the minmimum spanning tree - IMP_LOG(VERBOSE, "calc spanning tree"); - int max_tree_id = 0; - - // initially there are no trees - for (int i = 0; i < max_type_; i++) { - tree_id_[i] = 0; - } - - // start with the lowest score_ edges and work our way up - Float score = 0.0; - int num_edges = 0; - int type1, type2; - for (rs_iter = rsr_scores_.begin(); - (rs_iter != rsr_scores_.end()) && (num_edges < num_types_ - 1); - ++rs_iter) { - type1 = rs_iter->part1_type_; - type2 = rs_iter->part2_type_; - - IMP_LOG(VERBOSE, "Test for " << type1 << " " << type2); - // if neither particle is in a tree, create a new tree - if (!tree_id_[type1] && !tree_id_[type2]) { - max_tree_id++; - tree_id_[type1] = max_tree_id; - tree_id_[type2] = max_tree_id; - - num_edges++; - IMP_LOG(VERBOSE, "Evaluate for " << type1 << " " << type2); - score += rs_iter->rsr_->evaluate(accum); - } - - // if only one particle is already in a tree, add the other particle - // to that tree - else if (!tree_id_[type1] || !tree_id_[type2]) { - if (tree_id_[type1]) { - tree_id_[type2] = tree_id_[type1]; - } else { - tree_id_[type1] = tree_id_[type2]; + IMP_CHECK_OBJECT(ps_.get()); + unsigned int num_sets= set_offsets_.size()-1; + typedef std::vector > Edges; + typedef std::map EMap; + EMap edges; + for (unsigned int i=0; i< num_sets; ++i) { + for (unsigned int j=0; j< i; ++j) { + Pair cp(i,j); + for (unsigned int ii= set_offsets_[i]; ii != set_offsets_[i+1]; ++ii) { + for (unsigned int ij= set_offsets_[j]; ij != set_offsets_[j+1]; ++ij) { + Pair pp(ii, ij); + float d= ps_->evaluate(get_particle(ii), get_particle(ij), NULL); + edges[cp].push_back(std::make_pair(pp,d)); + } } - - num_edges++; - IMP_LOG(VERBOSE, "Evaluate for " << type1 << " " << type2); - score += rs_iter->rsr_->evaluate(accum); } + } - // both particles are already in trees - // if they are already in the same tree, do nothing - // otherwise, merge the two trees - else { - if (tree_id_[type1] != tree_id_[type2]) { - int old_tree_num = tree_id_[type1]; - int new_tree_num = tree_id_[type2]; - - for (int j = 0; j < max_type_; j++) { - if (tree_id_[j] == old_tree_num) { - tree_id_[j] = new_tree_num; - } + std::set in; + in.insert(0); + Edges picked_edges; + while (in.size() != num_sets) { + float best=std::numeric_limits::max(); + Pair best_edge, best_sets; + IMP_assert(!edges.empty(), "No more edges :-("); + for (EMap::iterator it = edges.begin(); it != edges.end(); ++it) { + Pair pit= it->first; + /*IMP_LOG(VERBOSE, "Looking at edges between " << pit[0] << " " + << pit[1] << std::endl);*/ + IMP_assert(in.find(pit[0]) == in.end() + || in.find(pit[1]) == in.end(), + "Edge only involves already visited nodes"); + if (in.find(pit[0]) == in.end() && in.find(pit[1]) == in.end()) { + // want to maintain a connected tree + continue; + } + for (unsigned int i=0; i< it->second.size(); ++i) { + if (it->second[i].second < best) { + best= it->second[i].second; + best_edge= it->second[i].first; + IMP_assert(best_edge != Pair(), "Bad pair"); + best_sets= it->first; } - - num_edges++; - IMP_LOG(VERBOSE, "Evaluate for " << type1 << " " << type2); - score += rs_iter->rsr_->evaluate(accum); } } + int newi=best_sets[0]; + if (in.find(newi) != in.end()) newi= best_sets[1]; + + IMP_LOG(VERBOSE, "Adding set " << newi << " with score " + << best << std::endl); + picked_edges.push_back(std::make_pair(best_edge, best)); + in.insert(newi); + IMP_assert(edges.find(best_sets) != edges.end(), + "Where did they go?"); + /*IMP_LOG(VERBOSE, "Removing edges " << best_sets[0] << " " + << best_sets[1] << std::endl);*/ + edges.erase(best_sets); } - return score; + float sum=0; + // could be more clever if accum is NULL + for (unsigned int i=0; i< picked_edges.size(); ++i) { + sum+= ps_->evaluate(get_particle(picked_edges[i].first[0]), + get_particle(picked_edges[i].first[1]), + accum); + } + return sum; } -/** - Show the current restraint. - - \param[in] out Stream to send restraint description to. - */ - void ConnectivityRestraint::show(std::ostream& out) const { if (get_is_active()) { @@ -325,7 +146,6 @@ out << "version: " << version() << " "; out << "last_modified_by: " << last_modified_by() << std::endl; - out << " num particles:" << num_particles_; } } // namespace IMP Index: kernel/test/connectivity/test_connectivity.py =================================================================== --- kernel/test/connectivity/test_connectivity.py (revision 349) +++ kernel/test/connectivity/test_connectivity.py (working copy) @@ -7,190 +7,74 @@ class ConnectivityTests(IMP.test.TestCase): """Class to test connectivity restraints""" - def setUp(self): - """Build test model and optimizer""" - self.imp_model = IMP.Model() - self.particles = [] - self.restraint_sets = [] - self.rsrs = [] + def min_distance(self, psa, psb): + md=100 + for pa in psa: + for pb in psb: + da= IMP.XYZDecorator.cast(pa) + db= IMP.XYZDecorator.cast(pb) + d= IMP.distance(da,db) + if (d < md): + md=d + return md - for p in range(12): - self.particles.append(IMP.utils.XYZParticle(self.imp_model, - 0., 0., 0.)) - p1 = self.particles[0] - - p1.add_attribute(IMP.FloatKey("radius"), 1.0, False) - p1.add_attribute(IMP.IntKey("protein"), 1) - p1.add_attribute(IMP.IntKey("id"), 1) - - p1 = self.particles[1] - p1.add_attribute(IMP.FloatKey("radius"), 1.0, False) - p1.add_attribute(IMP.IntKey("protein"), 1) - p1.add_attribute(IMP.IntKey("id"), 2) - - p1 = self.particles[2] - p1.add_attribute(IMP.FloatKey("radius"), 1.0, False) - p1.add_attribute(IMP.IntKey("protein"), 1) - p1.add_attribute(IMP.IntKey("id"), 3) - - p1 = self.particles[3] - p1.add_attribute(IMP.FloatKey("radius"), 1.5, False) - p1.add_attribute(IMP.IntKey("protein"), 2) - p1.add_attribute(IMP.IntKey("id"), 4) - - p1 = self.particles[4] - p1.add_attribute(IMP.FloatKey("radius"), 1.5, False) - p1.add_attribute(IMP.IntKey("protein"), 2) - p1.add_attribute(IMP.IntKey("id"), 5) - - p1 = self.particles[5] - p1.add_attribute(IMP.FloatKey("radius"), 1.5, False) - p1.add_attribute(IMP.IntKey("protein"), 2) - p1.add_attribute(IMP.IntKey("id"), 6) - - p1 = self.particles[6] - p1.add_attribute(IMP.FloatKey("radius"), 1.5, False) - p1.add_attribute(IMP.IntKey("protein"), 2) - p1.add_attribute(IMP.IntKey("id"), 7) - - p1 = self.particles[7] - p1.add_attribute(IMP.FloatKey("radius"), 2.0, False) - p1.add_attribute(IMP.IntKey("protein"), 3) - p1.add_attribute(IMP.IntKey("id"), 8) - - p1 = self.particles[8] - p1.add_attribute(IMP.FloatKey("radius"), 2.0, False) - p1.add_attribute(IMP.IntKey("protein"), 3) - p1.add_attribute(IMP.IntKey("id"), 9) - - p1 = self.particles[9] - p1.add_attribute(IMP.FloatKey("radius"), 2.0, False) - p1.add_attribute(IMP.IntKey("protein"), 3) - p1.add_attribute(IMP.IntKey("id"), 10) - - p1 = self.particles[10] - p1.add_attribute(IMP.FloatKey("radius"), 2.0, False) - p1.add_attribute(IMP.IntKey("protein"), 3) - p1.add_attribute(IMP.IntKey("id"), 11) - - p1 = self.particles[11] - p1.add_attribute(IMP.FloatKey("radius"), 2.0, False) - p1.add_attribute(IMP.IntKey("protein"), 3) - p1.add_attribute(IMP.IntKey("id"), 12) - - self.opt = IMP.ConjugateGradients() - self.opt.set_threshold(1e-4) - self.opt.set_model(self.imp_model) - def test_connectivity(self): """Test connectivity restraint. All particles in a single protein should be connected, and all proteins should be connected, either directly or indirectly through other proteins.""" - self.randomize_particles(self.particles, 50.0) + IMP.set_log_level(IMP.VERBOSE) + m = IMP.Model() - rs = IMP.RestraintSet("connect") - self.restraint_sets.append(rs) - self.imp_model.add_restraint(rs) + rk= IMP.FloatKey("radius") - # add connectivity restraints - - particle_indexes = IMP.Ints() - rsrs = [] - - score_func_params_ub = IMP.BasicScoreFuncParams("harmonic_upper_bound", - 0.0, 0.1) - - # set up exclusion volumes - IMP.utils.set_up_exclusion_volumes(self.imp_model, self.particles, - IMP.FloatKey("radius"), rsrs) - - # connect 3 proteins together - particle_indexes.clear() + p0=IMP.Particles() + p1=IMP.Particles() + p2=IMP.Particles() for i in range(12): - particle_indexes.push_back(i) - rsrs.append(IMP.ConnectivityRestraint(self.imp_model, particle_indexes, - IMP.IntKey("protein"), IMP.FloatKey("radius"), - score_func_params_ub)) + p= IMP.Particle() + m.add_particle(p) + d= IMP.XYZDecorator.create(p) + d.set_coordinates_are_optimized(True) + if (i %3 == 0): + p0.append(p) + p.add_attribute(rk, 1); + elif (i %3 == 1): + p1.append(p) + p.add_attribute(rk, 2); + else : + p2.append(p) + p.add_attribute(rk, 3); - # connect particles in protein1 together - particle_indexes.clear() - for i in range(3): - particle_indexes.push_back(i) - rsrs.append(IMP.ConnectivityRestraint(self.imp_model, particle_indexes, - IMP.IntKey("id"), IMP.FloatKey("radius"), - score_func_params_ub)) + o = IMP.ConjugateGradients() + o.set_threshold(1e-4) + o.set_model(m) + self.randomize_particles(m.get_particles(), 50.0) - # connect particles in protein2 together - particle_indexes.clear() - for i in range(3, 7): - particle_indexes.push_back(i) - rsrs.append(IMP.ConnectivityRestraint(self.imp_model, particle_indexes, - IMP.IntKey("id"), IMP.FloatKey("radius"), - score_func_params_ub)) + # add connectivity restraints - # connect particles in protein3 together - particle_indexes.clear() - for i in range(7, 12): - particle_indexes.push_back(i) - rsrs.append(IMP.ConnectivityRestraint(self.imp_model, particle_indexes, - IMP.IntKey("id"), IMP.FloatKey("radius"), - score_func_params_ub)) + ub = IMP.HarmonicUpperBound(0.0, 0.1) + ss= IMP.SphereDistancePairScore(ub) + r= IMP.ConnectivityRestraint(ss) + m.add_restraint(r) + r.add_set(p0) + r.add_set(p1) + r.add_set(p2) + o.optimize(1000) + d01= self.min_distance(p0, p1) + d02= self.min_distance(p0, p2) + d12= self.min_distance(p1, p2) + ok01= (d01 < 3.5) + ok02= (d02 < 4.5) + ok12= (d12 < 5.5) + print d01 + print d02 + print d12 + score= m.evaluate(False) + self.assert_(ok01 and ok02 or ok01 and ok12 or ok12 and ok02, + "Particles are too far") + self.assert_(score < 10, "Score too high"); - # add restraints - for i in range(len(rsrs)): - rs.add_restraint(rsrs[i]) - self.randomize_particles(self.particles, 50.0) - self.opt.optimize(55) - - # min distances - for i in range(len(self.particles)): - p = self.particles[i] - icoord = (p.x(), p.y(), p.z()) - irad = p.get_value(IMP.FloatKey("radius")) - for j in range(i+1,len(self.particles)): - p = self.particles[j] - jcoord = (p.x(), p.y(), p.z()) - jrad = p.get_value(IMP.FloatKey("radius")) - self.assert_(self.check_min_distance(icoord, jcoord, - irad + jrad - 0.05), - "min distance for any pair condition") - - # max distances - d12 = 10000000 - d13 = 10000000 - d23 = 10000000 - for i in range(len(self.particles)): - p = self.particles[i] - icoord = (p.x(), p.y(), p.z()) - irad = p.get_value(IMP.FloatKey("radius")) - t1 = p.get_value(IMP.IntKey("protein")) - for j in range(i+1,len(self.particles)): - p = self.particles[j] - jcoord = (p.x(), p.y(), p.z()) - jrad = p.get_value(IMP.FloatKey("radius")) - t2 = p.get_value(IMP.IntKey("protein")) - d = self.get_distance(icoord, jcoord) - irad - jrad - if t1 == 1 and t2 == 2: - if d < d12: - d12 = d - if t1 == 1 and t2 == 3: - if d < d13: - d13 = d - if t1 == 2 and t2 == 3: - if d < d23: - d23 = d - - sum = 0; - max_dist = 0.05 - if d12 < max_dist: - sum = sum + 1 - if d13 < max_dist: - sum = sum + 1 - if d23 < max_dist: - sum = sum + 1 - self.assert_(sum >= 2, "min spanning tree for three particle types") - if __name__ == '__main__': unittest.main()