From: Anna Dehof <anne@bioinf.uni-sb.de>
Date: Mon, 26 Mar 2012 19:11:39 +0200
Subject: Added MAX_PENALTY option to bond order assignment

---
 include/BALL/DATATYPE/GRAPH/treeWidth.h            |   28 +--
 .../STRUCTURE/BONDORDERS/FPTBondOrderStrategy.h    |   16 +-
 include/BALL/STRUCTURE/assignBondOrderProcessor.h  |   30 ++-
 source/APPLICATIONS/UTILITIES/assign_bond_orders.C |   34 +++-
 source/DATATYPE/GRAPH/treeWidth.C                  |    2 +-
 source/STRUCTURE/BONDORDERS/FPTBondOrderStrategy.C |  199 +++++++++-----------
 source/STRUCTURE/assignBondOrderProcessor.C        |   44 +++--
 7 files changed, 186 insertions(+), 167 deletions(-)

diff --git a/include/BALL/DATATYPE/GRAPH/treeWidth.h b/include/BALL/DATATYPE/GRAPH/treeWidth.h
index f7687ef..53ecffc 100644
--- a/include/BALL/DATATYPE/GRAPH/treeWidth.h
+++ b/include/BALL/DATATYPE/GRAPH/treeWidth.h
@@ -100,14 +100,14 @@ namespace BALL
 
 			typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS,
 																		  boost::property<boost::vertex_bag_content_t, std::set<OriginalVertexType>,
-																			boost::property<boost::vertex_bag_special_t, OriginalVertexType, 
-																			boost::property<boost::vertex_bag_type_t, int> > >, 
+																			boost::property<boost::vertex_bag_special_t, OriginalVertexType,
+																			boost::property<boost::vertex_bag_type_t, int> > >,
 																	  boost::no_property> TreeDecompositionGraph;
 
 			typedef typename boost::graph_traits<TreeDecompositionGraph>::vertex_descriptor TreeDecompositionBag;
 
 			typedef boost::iterator_property_map<typename std::vector<TreeDecompositionBag>::iterator,
-							                             typename boost::property_map<TreeDecompositionGraph, boost::vertex_index_t>::type> 
+							                             typename boost::property_map<TreeDecompositionGraph, boost::vertex_index_t>::type>
 																					 TreeDecompositionParentMap;
 			typedef boost::graph_as_tree<TreeDecompositionGraph, TreeDecompositionParentMap> TreeDecomposition;
 
@@ -131,7 +131,7 @@ namespace BALL
 			class ComponentFilter_
 			{
 				public:
-					ComponentFilter_(ComponentMap cm, Position i) 
+					ComponentFilter_(ComponentMap cm, Position i)
 						: cm_(cm),
 							component_(i)
 					{ }
@@ -144,7 +144,7 @@ namespace BALL
 
 				protected:
 					ComponentMap  cm_;
-					Position      component_;	
+					Position      component_;
 			};
 
 			/** PropertyWriter for graphviz output.
@@ -167,7 +167,7 @@ namespace BALL
 			// TODO: would UndirectedGraph suffice here?
 			MolecularGraph const* input_;
 
-			std::vector<boost::shared_ptr<EditableGraph> > components_; 
+			std::vector<boost::shared_ptr<EditableGraph> > components_;
 
 			std::vector<boost::shared_ptr<TreeDecomposition> >      nice_tree_decompositions_;
 			std::vector<boost::shared_ptr<TreeDecompositionGraph> > nice_tree_decomposition_graphs_;
@@ -201,7 +201,7 @@ namespace BALL
 		 * @tparam Reducer the reducer which removes a vertex from the graph to reduce it's size
 		 */
 		template<class Criterion, class Reducer>
-		class GeneralLowerBoundAlgorithm 
+		class GeneralLowerBoundAlgorithm
 			: public UnaryFunctor<UndirectedGraph, Size>
 		{
 			public:
@@ -284,7 +284,7 @@ namespace BALL
 		 * @throw BALL::GRAPH::UnconnectedGraphException if called on unconnected graphs
 		 */
 		template<class Criterion>
-		class GreedyX 
+		class GreedyX
 			: public UnaryFunctor<UndirectedGraph, typename std::pair<
 				std::vector<boost::graph_traits<typename UndirectedGraph::vertex_descriptor> >, Size> >
 		{
@@ -296,7 +296,7 @@ namespace BALL
 		 * A criterium for GreedyFillIn which search for a vertex with
 		 * minimum number of additional edges after eliminating
 		 */
-		struct FillInHeuristic 
+		struct FillInHeuristic
 		{
 			VertexType& operator() (UndirectedGraph& graph);
 
@@ -323,8 +323,8 @@ namespace BALL
 				 */
 				enum SIMPLICIAL_TYPE
 				{
-					NOT_SIMPLICIAL, 
-					ALMOST_SIMPLICIAL, 
+					NOT_SIMPLICIAL,
+					ALMOST_SIMPLICIAL,
 					IS_SIMPLICIAL
 				};
 
@@ -484,7 +484,7 @@ namespace BALL
 				TreeDecompositionBag buildJoin_(TreeDecompositionBag node, TreeDecompositionBag left,
 				                                  TreeDecompositionBag right, bool do_forget);
 
-				TreeDecompositionBag buildSingle_(TreeDecompositionBag node, int node_type, 
+				TreeDecompositionBag buildSingle_(TreeDecompositionBag node, int node_type,
 				                                    TreeDecompositionBag child);
 
 				TreeDecompositionBag buildLinkage_(TreeDecompositionBag node, TreeDecompositionBag child);
@@ -492,8 +492,8 @@ namespace BALL
 				TreeDecompositionBag linkWithIntroduceNodes_(TreeDecompositionContent parent_set, TreeDecompositionBag child);
 				TreeDecompositionBag linkWithForgetNodes_   (TreeDecompositionContent parent_set, TreeDecompositionBag child);
 
-				TreeDecompositionBag branch_(TreeDecompositionBag node, int node_type, 
-				                               typename std::vector<TreeDecompositionBag>::iterator begin, 
+				TreeDecompositionBag branch_(TreeDecompositionBag node, int node_type,
+				                               typename std::vector<TreeDecompositionBag>::iterator begin,
 																			 typename std::vector<TreeDecompositionBag>::iterator end);
 
 				boost::shared_ptr<TreeDecomposition> tree_;
diff --git a/include/BALL/STRUCTURE/BONDORDERS/FPTBondOrderStrategy.h b/include/BALL/STRUCTURE/BONDORDERS/FPTBondOrderStrategy.h
index 55cef31..ba129b0 100644
--- a/include/BALL/STRUCTURE/BONDORDERS/FPTBondOrderStrategy.h
+++ b/include/BALL/STRUCTURE/BONDORDERS/FPTBondOrderStrategy.h
@@ -1117,6 +1117,11 @@ namespace BALL
 					Size max_heap_size_;
 
 					/**
+					 * The number of solutions produced so far.
+					 */
+					Size num_computed_solutions_;
+
+					/**
 					 * current upperbound. This algorithm will just iterate solutions which are better than this upperbound;
 					 */
 					Penalty upper_bound_;
@@ -1248,14 +1253,6 @@ namespace BALL
 							                    Size solution_number, Penalty upper_bound = infinite_penalty);
 
 					/**
-					 * Construct a DPBackTrackingCombiner with the given FPTBondOrder and the number of solutions
-					 * @param bondAssignments vector with the bond assignments. Call #compute before constructing
-					 * @param solutionNumber the maximum number of solutions you want to backtrack
-					 */
-					DPBackTrackingCombiner_(std::vector<FPTBondOrderAssignment_>& bond_assignments, 
-					                        Size solution_number, Penalty upper_bound = infinite_penalty);
-
-					/**
 					 * Copy constructor
 					 */
 					DPBackTrackingCombiner_(DPBackTrackingCombiner_ const& copy);
@@ -1407,9 +1404,6 @@ namespace BALL
 			 */
 			std::vector<std::vector<int> > const* atom_to_block_;
 
-			/// upper bound on the penalty values for which we compute solutions
-			Penalty upper_bound_;
-
 			/**
 			 * A shared pointer to the computing data, so that you can copy this instance without copying
 			 * all the computing data
diff --git a/include/BALL/STRUCTURE/assignBondOrderProcessor.h b/include/BALL/STRUCTURE/assignBondOrderProcessor.h
index f8065d6..0d14af8 100644
--- a/include/BALL/STRUCTURE/assignBondOrderProcessor.h
+++ b/include/BALL/STRUCTURE/assignBondOrderProcessor.h
@@ -183,17 +183,35 @@ namespace BALL
 
 				/** the maximal number of solutions to compute
 				 *
-				 *  If set to zero all optimal solutions will be computed.
+				 *  If set to zero all optimal or all up to \link MAX_PENALTY MAX_PENALTY \endlink 
+				 *  solutions will be computed.
 				 *
-				 *  @see Option::COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS
+				 *  @see Option::COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS 
+				 *  @see Option::MAX_PENALTY
 				 */
 				static const char* MAX_NUMBER_OF_SOLUTIONS;
 
-				/** compute also non-optimal solutions but not more than 
-				 *  \link MAX_NUMBER_OF_SOLUTIONS MAX_NUMBER_OF_SOLUTIONS \endlink solutions.
+				/** the maximal penalty score allowed
+				 *  
+				 *  This option respects option \link MAX_NUMBER_OF_SOLUTIONS MAX_NUMBER_OF_SOLUTIONS \endlink
+				 *  if specified.
 				 *
+				 *  If set to -1 this option will be ignored. 
+				 *
+				 *  @see Option::COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS
 				 *  @see Option::MAX_NUMBER_OF_SOLUTIONS
 				 */
+				static const char* MAX_PENALTY;
+
+				/** compute also sub-optimal solutions but not more than 
+				 *  \link MAX_NUMBER_OF_SOLUTIONS MAX_NUMBER_OF_SOLUTIONS \endlink solutions.
+				 *  
+				 *  Alternatively \link Option::MAX_PENALTY  Option::MAX_PENALTY  \endlink
+				 *  allows to specify a maximal penalty. 
+				 *
+				 *  @see Option::MAX_NUMBER_OF_SOLUTIONS 
+				 *  @see Option::MAX_PENALTY
+				 */
 				static const char* COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS;
 
 				/** weighting of bond length penalties wrt valence penalties.
@@ -227,6 +245,7 @@ namespace BALL
 				static const String INIFile;
 				static const int MAX_BOND_ORDER;
 				static const int MAX_NUMBER_OF_SOLUTIONS;
+				static const int MAX_PENALTY;
 				static const bool COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS;
 				static const float BOND_LENGTH_WEIGHTING;
 				static const bool APPLY_FIRST_SOLUTION;
@@ -716,6 +735,9 @@ namespace BALL
 			// the max number of solutions to compute 
 			int max_number_of_solutions_;
 
+			// the max penalty score 
+			int max_penalty_;
+
 			// flag to indicate, whether also non-optimal solutions should be computed 
 			bool compute_also_non_optimal_solutions_;
 
diff --git a/source/APPLICATIONS/UTILITIES/assign_bond_orders.C b/source/APPLICATIONS/UTILITIES/assign_bond_orders.C
index cefef9b..7053baa 100644
--- a/source/APPLICATIONS/UTILITIES/assign_bond_orders.C
+++ b/source/APPLICATIONS/UTILITIES/assign_bond_orders.C
@@ -16,10 +16,10 @@ using namespace BALL;
 
 int main(int argc, char** argv)
 {
-	if ((argc < 3) || (argc > 4))
+	if (argc < 3)
 	//	if (argc != 3)
 	{
-		Log << "Usage:" << argv[0] << " <Molecule infile> <Molecule outfile> [num_sol]" << endl;
+		Log << "Usage:" << argv[0] << " <Molecule infile> <Molecule outfile> [--num-sol=n] [--strategy=(fpt|a_star|ilp)] [--max-penalty=m]" << endl;
 		return 1;
 	}
 
@@ -49,24 +49,40 @@ int main(int argc, char** argv)
 	//    set the options:
 	//
 	// the solution strategy (A*, ILP or FPT)
-	abop.options.set(AssignBondOrderProcessor::Option::ALGORITHM, AssignBondOrderProcessor::Algorithm::A_STAR);
+	abop.options.set(AssignBondOrderProcessor::Option::ALGORITHM, AssignBondOrderProcessor::Algorithm::FPT);
 
 	// specify the inifile with the atomic valence penalties
 	abop.options.set(AssignBondOrderProcessor::Option::INIFile, AssignBondOrderProcessor::Default::INIFile);
 
 	// options for considering bond length as well 
-	abop.options.setReal(AssignBondOrderProcessor::Option::BOND_LENGTH_WEIGHTING, 0);
-	abop.options.setReal(AssignBondOrderProcessor::Option::USE_FINE_PENALTY, true);
+//	abop.options.setReal(AssignBondOrderProcessor::Option::BOND_LENGTH_WEIGHTING, 0);
+//	abop.options.setReal(AssignBondOrderProcessor::Option::USE_FINE_PENALTY, true);
 
 	// the combination of the following two options causes the computation of all optimal solutions
 	abop.options.setInteger(AssignBondOrderProcessor::Option::MAX_NUMBER_OF_SOLUTIONS, 0);
 	abop.options.setBool(AssignBondOrderProcessor::Option::COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS, false);
 
-	if (argc == 4)
+	if (argc > 3)
 	{
-		Log << "  with compute fixed number (" << argv[3] << ") solutions" << endl;
-		abop.options.setInteger(AssignBondOrderProcessor::Option::MAX_NUMBER_OF_SOLUTIONS, String(argv[3]).toInt());
-		abop.options.setBool(AssignBondOrderProcessor::Option::COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS, true);
+		for (Index i=3; i<argc; ++i)
+		{
+			if (String(argv[i]).hasPrefix("--strategy"))
+			{
+				abop.options.set(AssignBondOrderProcessor::Option::ALGORITHM, String(argv[i]).getField(1, "="));
+			}
+			else if (String(argv[i]).hasPrefix("--num-sol"))
+			{
+				Log << "  with compute fixed number (" << argv[i] << ") solutions" << endl;
+				abop.options.setInteger(AssignBondOrderProcessor::Option::MAX_NUMBER_OF_SOLUTIONS, String(argv[i]).getField(1, "=").toInt());
+				abop.options.setBool(AssignBondOrderProcessor::Option::COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS, true);
+			}
+			else if (String(argv[i]).hasPrefix("--max-penalty"))
+			{
+				abop.options.setInteger(AssignBondOrderProcessor::Option::MAX_PENALTY, String(argv[i]).getField(1, "=").toInt());
+				abop.options.setBool(AssignBondOrderProcessor::Option::COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS, true);
+				Log << "  with maximum penalty (" <<  abop.options.getInteger(AssignBondOrderProcessor::Option::MAX_PENALTY) << ")" << endl;
+			}
+		}
 	}
 
 	// define input and output properties
diff --git a/source/DATATYPE/GRAPH/treeWidth.C b/source/DATATYPE/GRAPH/treeWidth.C
index 7485305..cd0a879 100644
--- a/source/DATATYPE/GRAPH/treeWidth.C
+++ b/source/DATATYPE/GRAPH/treeWidth.C
@@ -12,7 +12,7 @@ namespace BALL
 		// find all vertices in the current bag
 		TreeDecompositionContent content = boost::get(boost::vertex_bag_content, *td_, v);
 
-		for (typename TreeDecompositionContent::const_iterator tdc_it = content.begin(); tdc_it != content.end(); ++tdc_it)
+		for (TreeDecompositionContent::const_iterator tdc_it = content.begin(); tdc_it != content.end(); ++tdc_it)
 		{
 			TreeWidth<MolecularGraph>::OriginalVertexType ov = *tdc_it;
 			Atom const* atom = boost::get(boost::vertex_atom_ptr, *original_graph_, ov);
diff --git a/source/STRUCTURE/BONDORDERS/FPTBondOrderStrategy.C b/source/STRUCTURE/BONDORDERS/FPTBondOrderStrategy.C
index 42040e5..107823e 100644
--- a/source/STRUCTURE/BONDORDERS/FPTBondOrderStrategy.C
+++ b/source/STRUCTURE/BONDORDERS/FPTBondOrderStrategy.C
@@ -12,10 +12,6 @@ namespace BALL
 	const FPTBondOrderStrategy::Penalty FPTBondOrderStrategy::infinite_penalty = 1e5;
 	const FPTBondOrderStrategy::Valence FPTBondOrderStrategy::max_valence      = 8;
 
-	String  FPTBondOrderStrategy::Option::UPPER_PENALTY_BOUND  = "upper_penalty_bound";
-
-	FPTBondOrderStrategy::Penalty FPTBondOrderStrategy::Default::UPPER_PENALTY_BOUND = infinite_penalty;
-
 	FPTBondOrderStrategy::FPTBondOrderStrategy(AssignBondOrderProcessor* parent)
 		: BondOrderAssignmentStrategy(parent)
 	{
@@ -33,15 +29,11 @@ namespace BALL
 
 	bool FPTBondOrderStrategy::readOptions(const Options& options)
 	{
-		upper_bound_ = options.get(Option::UPPER_PENALTY_BOUND).toFloat();
-
 		return true;
 	}
 
 	void FPTBondOrderStrategy::setDefaultOptions()
 	{
-		abop->options.setDefault(FPTBondOrderStrategy::Option::UPPER_PENALTY_BOUND,
-		                         FPTBondOrderStrategy::Default::UPPER_PENALTY_BOUND);
 	}
 
 	void FPTBondOrderStrategy::init()
@@ -64,16 +56,22 @@ namespace BALL
 		std::vector<boost::shared_ptr<TreeDecomposition> > & ntds = computing_data_->tw->getNiceTreeDecompositions();
 		std::vector<FPTBondOrderAssignment_*>& bond_assignments   = computing_data_->bond_assignments;
 
+		int max_penalty = abop->max_penalty_;
+		if (max_penalty == -1)
+			max_penalty = infinite_penalty;
+		else
+			max_penalty += 1;
+
 		bond_assignments.reserve(ntds.size());
 		for (Position i = 0; i < ntds.size(); ++i)
 		{
-			bond_assignments.push_back(new FPTBondOrderAssignment_(*this, ntds[i], upper_bound_));
-		  bond_assignments[i]->compute();
+			bond_assignments.push_back(new FPTBondOrderAssignment_(*this, ntds[i], max_penalty));
+		  Penalty result = bond_assignments[i]->compute();
 		}
 
 		// initialize backtracking
 		combiner_ = boost::shared_ptr<DPBackTrackingCombiner_>(
-		                   new DPBackTrackingCombiner_(bond_assignments, abop->max_number_of_solutions_, upper_bound_));
+		                   new DPBackTrackingCombiner_(bond_assignments, abop->max_number_of_solutions_, max_penalty));
 	}
 
 	boost::shared_ptr<BondOrderAssignment> FPTBondOrderStrategy::computeNextSolution()
@@ -144,7 +142,7 @@ namespace BALL
 			int penalty = (*penalties_)[start_index + position];
 
 			return (penalty < 0) ? infinite_penalty : static_cast<float>(penalty);
-		} 
+		}
 		else
 		{
 			return infinite_penalty;
@@ -155,20 +153,20 @@ namespace BALL
 	//*           DPConfig_                                                                               *
 	//*****************************************************************************************************
 
-	FPTBondOrderStrategy::DPConfig_::DPConfig_() 
-		: consumed_valences(0), 
+	FPTBondOrderStrategy::DPConfig_::DPConfig_()
+		: consumed_valences(0),
 		  bond_assignments(0)
 	{
 	}
 
-	FPTBondOrderStrategy::DPConfig_::DPConfig_(Size atoms, Size bonds) 
-		: consumed_valences(atoms, 0), 
+	FPTBondOrderStrategy::DPConfig_::DPConfig_(Size atoms, Size bonds)
+		: consumed_valences(atoms, 0),
 		  bond_assignments(bonds, 0)
 	{
 	}
 
-	FPTBondOrderStrategy::DPConfig_::DPConfig_(std::vector<Valence> const& v, std::vector<BondOrder> const& bo) 
-		: consumed_valences(v), 
+	FPTBondOrderStrategy::DPConfig_::DPConfig_(std::vector<Valence> const& v, std::vector<BondOrder> const& bo)
+		: consumed_valences(v),
 		  bond_assignments(bo)
 	{
 	}
@@ -258,12 +256,12 @@ namespace BALL
 	//*           DPTable_                                                                                *
 	//*****************************************************************************************************
 
-	FPTBondOrderStrategy::DPTable_::DPTable_() 
+	FPTBondOrderStrategy::DPTable_::DPTable_()
 	  : table()
 	{
 	}
 
-	FPTBondOrderStrategy::DPTable_::DPTable_(DPTable_ const& dptab) 
+	FPTBondOrderStrategy::DPTable_::DPTable_(DPTable_ const& dptab)
 		: table(dptab.table)
 	{
 	}
@@ -345,7 +343,7 @@ namespace BALL
 		if (insertion.second)
 		{
 			return true;
-		} 
+		}
 		else
 		{
 			if (insertion.first->second > penalty)
@@ -353,7 +351,7 @@ namespace BALL
 				insertion.first->second = penalty;
 
 				return true;
-			} 
+			}
 			else
 			{
 				return false;
@@ -365,14 +363,14 @@ namespace BALL
 	//*           AdditionalBagProperties_                                                                *
 	//*****************************************************************************************************
 
-	FPTBondOrderStrategy::AdditionalBagProperties_::AdditionalBagProperties_() 
-		: bonds(), 
-		  table(new DPTable_()) 
+	FPTBondOrderStrategy::AdditionalBagProperties_::AdditionalBagProperties_()
+		: bonds(),
+		  table(new DPTable_())
 	{
 	}
 
-	FPTBondOrderStrategy::AdditionalBagProperties_::AdditionalBagProperties_(AdditionalBagProperties_ const& copy) 
-		: bonds(copy.bonds), 
+	FPTBondOrderStrategy::AdditionalBagProperties_::AdditionalBagProperties_(AdditionalBagProperties_ const& copy)
+		: bonds(copy.bonds),
 			table(new DPTable_(*copy.table))
 	{
 	}
@@ -384,7 +382,7 @@ namespace BALL
 
 	FPTBondOrderStrategy::AdditionalBagProperties_& FPTBondOrderStrategy::AdditionalBagProperties_::operator=(AdditionalBagProperties_ const& copy)
 	{
-		if (&copy != this) 
+		if (&copy != this)
 		{
 			bonds = copy.bonds;
 			table = new DPTable_(*copy.table);
@@ -447,7 +445,7 @@ namespace BALL
 				{
 					throw Exception::NullPointer(__FILE__, __LINE__);
 				}
-			case TreeWidth<MolecularGraph>::LEAF_BAG: 
+			case TreeWidth<MolecularGraph>::LEAF_BAG:
 				break;
 			case TreeWidth<MolecularGraph>::END_BAG:
 				if (boost::num_vertices(*molecule_) == 0 && (bag == root(*ntd_)))
@@ -455,7 +453,7 @@ namespace BALL
 					// empty molecule -> empty table
 					bag_properties.table->insert(DPConfig_(), 0);
 					return bag_properties.table;
-				} 
+				}
 				else
 				{
 					// else nice tree decomposition is damaged
@@ -603,7 +601,7 @@ namespace BALL
 		return bonds;
 	}
 
-	void FPTBondOrderStrategy::FPTBondOrderAssignment_::computeIntroduceBag(TreeDecompositionBag& bag, DPTable_& child, 
+	void FPTBondOrderStrategy::FPTBondOrderAssignment_::computeIntroduceBag(TreeDecompositionBag& bag, DPTable_& child,
 	                                                                        AdditionalBagProperties_& property)
 	{
 		typedef TreeWidth<MolecularGraph>::TreeDecompositionContent TreeDecompositionContent;
@@ -614,7 +612,7 @@ namespace BALL
 
 		TreeDecompositionContent vertices = boost::get(boost::vertex_bag_content, *ntd_, bag);
 		Size num_valences = vertices.size();
-		
+
 		std::vector<MolecularGraphTraits::EdgeType>& bonds = property.bonds;
 		Size num_bonds = bonds.size();
 
@@ -654,7 +652,7 @@ namespace BALL
 				{
 					// remember the indizes of the new introduced bonds
 					indices.push_back(vindex);
-				} 
+				}
 				else
 				{
 					conf.bond_assignments[vindex] = child_entry.first.get().bond_assignments[cindex++];
@@ -662,7 +660,7 @@ namespace BALL
 			}
 
 			// if there are any introduced bonds we have to fill them with each possible value
-			if (indices.size() > 0) 
+			if (indices.size() > 0)
 			{
 				std::vector<int> bond_values(indices.size());
 
@@ -685,7 +683,7 @@ namespace BALL
 						++current_index;
 					}
 				}
-			} 
+			}
 			else
 			{
 				table.insert(conf, child_entry.second);
@@ -768,7 +766,7 @@ namespace BALL
 		return parent_->getPenaltyFor_(forgotten_vertex, forgotten_valence) + child_row.second;
 	}
 
-	void FPTBondOrderStrategy::FPTBondOrderAssignment_::computeForgetBag(TreeDecompositionBag& bag, DPTable_& child, 
+	void FPTBondOrderStrategy::FPTBondOrderAssignment_::computeForgetBag(TreeDecompositionBag& bag, DPTable_& child,
 																					                             AdditionalBagProperties_& property)
 	{
 		typedef TreeWidth<MolecularGraph>::TreeDecompositionContent TreeDecompositionContent;
@@ -795,7 +793,7 @@ namespace BALL
 		}
 	}
 
-	void FPTBondOrderStrategy::FPTBondOrderAssignment_::computeRootBag(TreeDecompositionBag& bag, 
+	void FPTBondOrderStrategy::FPTBondOrderAssignment_::computeRootBag(TreeDecompositionBag& bag,
 	                                                                   DPTable_& child, AdditionalBagProperties_& bag_properties)
 	{
 		DPConfig_ empty(0, 0);
@@ -816,7 +814,8 @@ namespace BALL
 	}
 
 	void FPTBondOrderStrategy::FPTBondOrderAssignment_::computeJoinBag(TreeDecompositionBag& bag,
-                                  			DPTable_& left_child, DPTable_& right_child, AdditionalBagProperties_& property)
+			                                                               DPTable_& left_child, DPTable_& right_child, 
+																																		 AdditionalBagProperties_& property)
 	{
 		typedef TreeWidth<MolecularGraph>::TreeDecompositionContent TreeDecompositionContent;
 
@@ -1268,22 +1267,23 @@ namespace BALL
 	//*           DPBackTracking_                                                                        *
 	//*****************************************************************************************************
 
-	FPTBondOrderStrategy::DPBackTracking_::DPBackTracking_(FPTBondOrderAssignment_& bond_assignment, Size max_num_solutions, 
-                                                         std::vector<MolecularGraphTraits::EdgeType> const& bonds, Penalty upper_bound) 
-		: bond_assignment_(&bond_assignment), 
-			current_state_(NULL), 
-			queue_(), 
-			max_num_solutions_(max_num_solutions), 
+	FPTBondOrderStrategy::DPBackTracking_::DPBackTracking_(FPTBondOrderAssignment_& bond_assignment, Size max_num_solutions,
+                                                         std::vector<MolecularGraphTraits::EdgeType> const& bonds, Penalty upper_bound)
+		: bond_assignment_(&bond_assignment),
+			current_state_(NULL),
+			queue_(),
+			max_num_solutions_(max_num_solutions),
 			bonds_(&bonds),
-			bags_(boost::shared_ptr<std::vector<TreeDecompositionBag> >(new std::vector<TreeDecompositionBag>)), 
-			max_heap_size_(max_num_solutions), 
+			bags_(boost::shared_ptr<std::vector<TreeDecompositionBag> >(new std::vector<TreeDecompositionBag>)),
+			max_heap_size_(1e6),
+			num_computed_solutions_(0),
 			upper_bound_(upper_bound)
 	{
 		// order bags in preorder
 		if (bond_assignment.ntd_ == NULL)
 		{
 			throw Exception::NullPointer(__FILE__, __LINE__);
-		} 
+		}
 		else
 		{
 			bags_->reserve(boost::num_vertices(bond_assignment.ntd_->_g));
@@ -1303,14 +1303,15 @@ namespace BALL
 		}
 	}
 
-	FPTBondOrderStrategy::DPBackTracking_::DPBackTracking_(DPBackTracking_ const& copy) 
-	 : bond_assignment_(copy.bond_assignment_), 
-		 current_state_(NULL), 
+	FPTBondOrderStrategy::DPBackTracking_::DPBackTracking_(DPBackTracking_ const& copy)
+	 : bond_assignment_(copy.bond_assignment_),
+		 current_state_(NULL),
 		 queue_(),
-		 max_num_solutions_(copy.max_num_solutions_), 
-		 bonds_(copy.bonds_), 
-		 bags_(copy.bags_), 
+		 max_num_solutions_(copy.max_num_solutions_),
+		 bonds_(copy.bonds_),
+		 bags_(copy.bags_),
 		 max_heap_size_(copy.max_heap_size_),
+		 num_computed_solutions_(copy.num_computed_solutions_),
 		 upper_bound_(copy.upper_bound_)
 	{
 		if (copy.current_state_ != NULL)
@@ -1345,6 +1346,7 @@ namespace BALL
 			bonds_ = copy.bonds_;
 			bags_ = copy.bags_;
 			max_heap_size_ = copy.max_heap_size_;
+			num_computed_solutions_ = copy.num_computed_solutions_;
 			upper_bound_ = copy.upper_bound_;
 		}
 
@@ -1365,6 +1367,7 @@ namespace BALL
 
 		std::vector<BackTrackingState_*> copy(queue_.begin(), queue_.end());
 		queue_.clear();
+		num_computed_solutions_ = 0;
 
 		for (std::vector<BackTrackingState_*>::iterator iter = copy.begin(); iter != copy.end(); ++iter)
 		{
@@ -1372,7 +1375,7 @@ namespace BALL
 		}
 	}
 
-	bool FPTBondOrderStrategy::DPBackTracking_::StateComparator_::operator() (BackTrackingState_ const * left, 
+	bool FPTBondOrderStrategy::DPBackTracking_::StateComparator_::operator() (BackTrackingState_ const * left,
 	                                                                          BackTrackingState_ const * right) const
 	{
 		return *left < *right;
@@ -1397,7 +1400,7 @@ namespace BALL
 		if (order < properties.size())
 		{
 			return properties[order];
-		} 
+		}
 		else
 		{
 			throw Exception::IndexOverflow(__FILE__, __LINE__, static_cast<Index>(order), properties.size());
@@ -1426,12 +1429,12 @@ namespace BALL
 
 	bool FPTBondOrderStrategy::DPBackTracking_::hasMoreSolutions() const
 	{
-		return (!queue_.empty() && (max_heap_size_ > 0));
+		return (!queue_.empty() && (!max_num_solutions_ || (num_computed_solutions_ <= max_num_solutions_)));
 	}
 
 	void FPTBondOrderStrategy::DPBackTracking_::nextSolution()
 	{
-		if (queue_.empty() || max_heap_size_ == 0)
+		if (queue_.empty() || max_heap_size_ == 0 || ((max_num_solutions_ > 0) && (num_computed_solutions_ > max_num_solutions_)))
 		{
 			throw Exception::OutOfRange(__FILE__, __LINE__);
 		}
@@ -1444,7 +1447,7 @@ namespace BALL
 		std::multiset<BackTrackingState_*, StateComparator_>::iterator first = queue_.begin();
 		current_state_ = *first;
 		queue_.erase(first);
-		--max_heap_size_;
+		++num_computed_solutions_;
 
 		TreeDecomposition& tree = *bond_assignment_->ntd_;
 
@@ -1469,7 +1472,7 @@ namespace BALL
 				{
 					TreeDecomposition::children_iterator c_2 = c_i;
 					++c_2;
-					visitJoin(*current_state_, bag, getTable(boost::get(boost::vertex_index, tree, *c_i)), 
+					visitJoin(*current_state_, bag, getTable(boost::get(boost::vertex_index, tree, *c_i)),
 																				getTable(boost::get(boost::vertex_index, tree, *c_2)));
 					break;
 				}
@@ -1769,7 +1772,7 @@ namespace BALL
 		if (d < bonds_->size())
 		{
 			return d;
-		} 
+		}
 		else
 		{
 			throw Exception::IndexOverflow(__FILE__, __LINE__, static_cast<Index>(d), bonds_->size());
@@ -1795,7 +1798,7 @@ namespace BALL
 
 	bool FPTBondOrderStrategy::DPBackTracking_::isSolutionNeeded(Penalty penalty)
 	{
-		if (max_heap_size_ == 0) {return false;}
+		if (max_heap_size_ == 0 || ((max_num_solutions_ > 0) && (num_computed_solutions_ > max_num_solutions_))) {return false;}
 		if (queue_.size() >= max_heap_size_) {return Maths::isLess(penalty, upper_bound_);}
 		return Maths::isLessOrEqual(penalty, upper_bound_);
 	}
@@ -1846,7 +1849,7 @@ namespace BALL
 		if (queue_.empty() || max_heap_size_ == 0)
 		{
 			return FPTBondOrderStrategy::infinite_penalty;
-		} 
+		}
 		else
 		{
 			return (*queue_.begin())->assignment.penalty;
@@ -1857,25 +1860,25 @@ namespace BALL
 	//*           DPBackTrackingCombiner_                                                                 *
 	//*****************************************************************************************************
 
-	FPTBondOrderStrategy::DPBackTrackingCombiner_::DPBackTrackingCombiner_(DPBackTrackingCombiner_ const& copy) 
-		: backtrackers_(copy.deepCopyOfBacktrackers_()), 
+	FPTBondOrderStrategy::DPBackTrackingCombiner_::DPBackTrackingCombiner_(DPBackTrackingCombiner_ const& copy)
+		: backtrackers_(copy.deepCopyOfBacktrackers_()),
 		  priority_queue_(copy.priority_queue_),
-			component_solutions_(copy.component_solutions_), 
+			component_solutions_(copy.component_solutions_),
 			assignment_(copy.assignment_),
-			solution_number_(copy.solution_number_), 
-			optimum_(copy.optimum_), 
+			solution_number_(copy.solution_number_),
+			optimum_(copy.optimum_),
 			upper_bound_(copy.upper_bound_)
 	{
 	}
 
-	FPTBondOrderStrategy::DPBackTrackingCombiner_::DPBackTrackingCombiner_(std::vector<FPTBondOrderAssignment_*>& bond_assignments, 
-		                                                                     Size solution_number, Penalty upper_bound) 
-		: backtrackers_(), 
-		  priority_queue_(), 
+	FPTBondOrderStrategy::DPBackTrackingCombiner_::DPBackTrackingCombiner_(std::vector<FPTBondOrderAssignment_*>& bond_assignments,
+		                                                                     Size solution_number, Penalty upper_bound)
+		: backtrackers_(),
+		  priority_queue_(),
 			component_solutions_(bond_assignments.size()),
-			assignment_(), 
-			solution_number_(solution_number), 
-			optimum_(FPTBondOrderStrategy::infinite_penalty), 
+			assignment_(),
+			solution_number_(solution_number),
+			optimum_(FPTBondOrderStrategy::infinite_penalty),
 			upper_bound_(upper_bound)
 	{
 		if (bond_assignments.empty())
@@ -1897,7 +1900,7 @@ namespace BALL
 		FPTBondOrderStrategy::EdgeComparator_ ec(&graph);
 		std::sort(sorted_edges.begin(), sorted_edges.end(), ec);
 
-		for (std::vector<FPTBondOrderAssignment_*>::const_iterator iter  = bond_assignments.begin(); 
+		for (std::vector<FPTBondOrderAssignment_*>::const_iterator iter  = bond_assignments.begin();
 		                                                           iter != bond_assignments.end();   ++iter)
 		{
 			backtrackers_.push_back(new DPBackTracking_(**iter, solution_number_, sorted_edges, upper_bound_));
@@ -1906,39 +1909,6 @@ namespace BALL
 		initialize_();
 	}
 
-	FPTBondOrderStrategy::DPBackTrackingCombiner_::DPBackTrackingCombiner_(std::vector<FPTBondOrderAssignment_>& bond_assignments, 
-		                                             Size solution_number, Penalty upper_bound) 
-		: backtrackers_(), 
-		  priority_queue_(), 
-			component_solutions_(bond_assignments.size()),
-			assignment_(), 
-			solution_number_(solution_number), 
-			optimum_(FPTBondOrderStrategy::infinite_penalty), 
-			upper_bound_(upper_bound)
-	{
-		backtrackers_.reserve(bond_assignments.size());
-
-		MolecularGraph& graph = *bond_assignments[0].molecule_;
-		sorted_edges.reserve(boost::num_edges(graph));
-
-		BGL_FORALL_EDGES(edge_it, graph, MolecularGraph)
-		{
-			sorted_edges.push_back(edge_it);
-		}
-
-		// sort bonds - the second vertex could be in false order
-		FPTBondOrderStrategy::EdgeComparator_ ec(&graph);
-		std::sort(sorted_edges.begin(), sorted_edges.end(), ec);
-
-		for (std::vector<FPTBondOrderAssignment_>::iterator iter  = bond_assignments.begin(); 
-		                                                    iter != bond_assignments.end();   ++iter)
-		{
-			backtrackers_.push_back(new DPBackTracking_(*iter, solution_number_, sorted_edges, upper_bound_));
-		}
-
-		initialize_();
-	}
-
 	FPTBondOrderStrategy::DPBackTrackingCombiner_::~DPBackTrackingCombiner_()
 	{
 		clear();
@@ -1946,14 +1916,14 @@ namespace BALL
 
 	void FPTBondOrderStrategy::DPBackTrackingCombiner_::clear()
 	{
-		for (std::vector<DPBackTracking_*>::iterator iter  = backtrackers_.begin(); 
+		for (std::vector<DPBackTracking_*>::iterator iter  = backtrackers_.begin();
 		                                             iter != backtrackers_.end();   ++iter)
 		{
 			delete *iter;
 		}
 	}
 
-	FPTBondOrderStrategy::DPBackTrackingCombiner_& 
+	FPTBondOrderStrategy::DPBackTrackingCombiner_&
 	FPTBondOrderStrategy::DPBackTrackingCombiner_::operator = (DPBackTrackingCombiner_ const& copy)
 	{
 		if (this != &copy)
@@ -1972,7 +1942,7 @@ namespace BALL
 		return *this;
 	}
 
-	std::pair<Size, FPTBondOrderStrategy::Penalty> 
+	std::pair<Size, FPTBondOrderStrategy::Penalty>
 	FPTBondOrderStrategy::DPBackTrackingCombiner_::getNextMinimumBackTracker_() const
 	{
 		if (backtrackers_.size() == 1)
@@ -2034,7 +2004,7 @@ namespace BALL
 				assignment_ = ass;
 
 				return;
-			} 
+			}
 			else
 			{
 				throw Exception::OutOfRange(__FILE__, __LINE__);
@@ -2061,7 +2031,7 @@ namespace BALL
 			assignment_ = priority_queue_.top();
 
 			priority_queue_.pop();
-		} 
+		}
 		else
 		{
 			throw Exception::OutOfRange(__FILE__, __LINE__);
@@ -2151,7 +2121,7 @@ namespace BALL
 		return assignment_;
 	}
 
-	std::vector<FPTBondOrderStrategy::DPBackTracking_*> 
+	std::vector<FPTBondOrderStrategy::DPBackTracking_*>
 	FPTBondOrderStrategy::DPBackTrackingCombiner_::deepCopyOfBacktrackers_() const
 	{
 		std::vector<DPBackTracking_*> ts;
@@ -2166,5 +2136,4 @@ namespace BALL
 		return ts;
 	}
 
-
 }
diff --git a/source/STRUCTURE/assignBondOrderProcessor.C b/source/STRUCTURE/assignBondOrderProcessor.C
index f8acc4f..9863a54 100644
--- a/source/STRUCTURE/assignBondOrderProcessor.C
+++ b/source/STRUCTURE/assignBondOrderProcessor.C
@@ -108,6 +108,9 @@ namespace BALL
 	const char* AssignBondOrderProcessor::Option::MAX_NUMBER_OF_SOLUTIONS = "max_number_of_solutions";
 	const int  AssignBondOrderProcessor::Default::MAX_NUMBER_OF_SOLUTIONS = 10;
 
+	const char* AssignBondOrderProcessor::Option::MAX_PENALTY = "max_penalty_score";
+	const int  AssignBondOrderProcessor::Default::MAX_PENALTY = -1;
+
 	const char* AssignBondOrderProcessor::Option::COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS = "compute_also_non_optimal_solutions";
 	const bool  AssignBondOrderProcessor::Default::COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS = false;
 
@@ -148,6 +151,7 @@ namespace BALL
 			max_bond_order_(),
 			alpha_(),
 			max_number_of_solutions_(),
+			max_penalty_(-1),
 			compute_also_non_optimal_solutions_(),
 			add_missing_hydrogens_(),
 			compute_also_connectivity_(),
@@ -237,13 +241,14 @@ namespace BALL
 			ret &= s_it->second->readOptions(options);
 		}
 
-		max_bond_order_ = options.getInteger(Option::MAX_BOND_ORDER);
-		alpha_ = options.getReal(Option::BOND_LENGTH_WEIGHTING);
-		max_number_of_solutions_ = options.getInteger(Option::MAX_NUMBER_OF_SOLUTIONS);
+		max_bond_order_            = options.getInteger(Option::MAX_BOND_ORDER);
+		alpha_                     = options.getReal(Option::BOND_LENGTH_WEIGHTING);
+		max_number_of_solutions_   = options.getInteger(Option::MAX_NUMBER_OF_SOLUTIONS);
+		max_penalty_               = options.getInteger(Option::MAX_PENALTY);
 		compute_also_non_optimal_solutions_ = options.getBool(Option::COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS);
-		add_missing_hydrogens_ = options.getBool(Option::ADD_HYDROGENS);
+		add_missing_hydrogens_     = options.getBool(Option::ADD_HYDROGENS);
 		compute_also_connectivity_ = options.getBool(Option::COMPUTE_ALSO_CONNECTIVITY);
-		use_fine_penalty_ = options.getBool(Option::USE_FINE_PENALTY);
+		use_fine_penalty_          = options.getBool(Option::USE_FINE_PENALTY);
 
 		if (max_bond_order_ <= 0)
 		{
@@ -257,6 +262,12 @@ namespace BALL
 				          << " : Error in options! Please check the option Option::MAX_NUMBER_OF_SOLUTIONS."  << endl;
 			ret = false;
 		}
+		if (max_penalty_ < -1)
+		{
+			Log.error() << __FILE__ << " " << __LINE__
+				          << " : Error in options! Please check the option Option::MAX_PENALTY."  << endl;
+			ret = false;
+		}
 		if ((alpha_ < 0) || ((alpha_ > 1)))
 		{
 			Log.error() << __FILE__ << " " << __LINE__
@@ -280,14 +291,13 @@ namespace BALL
 						|| (options.getBool(Option::OVERWRITE_SELECTED_BONDS)  == true)
 						|| (options.getBool(Option::OVERWRITE_SINGLE_BOND_ORDERS) == false)
 						|| (options.getBool(Option::OVERWRITE_DOUBLE_BOND_ORDERS) == false)
-						|| (options.getBool(Option::OVERWRITE_TRIPLE_BOND_ORDERS) == false)
-						|| (options.getInteger(Option::MAX_NUMBER_OF_SOLUTIONS) == 0))
+						|| (options.getBool(Option::OVERWRITE_TRIPLE_BOND_ORDERS) == false) )
 			 )
 		{
 			Log.error() << __FILE__ << " " << __LINE__
 				          << " : Error in options! FPT cannot be used with these options." << endl
-									<< " Switch to solution strategy ASTAR by setting Option::ALGORITHM to Algorithm::ASTAR." << endl;
-
+									<< " Consider switch to solution strategy ASTAR by setting Option::ALGORITHM to Algorithm::ASTAR." << endl
+									<< " Abort." << endl;
 			ret = false;
 		}
 
@@ -328,6 +338,7 @@ cout << " \t Penalty file " << options[Option::Option::INIFile] << endl;
 cout << " \t alpha: " << options[Option::BOND_LENGTH_WEIGHTING] << endl;
 cout << " \t max bond order: " << options[Option::MAX_BOND_ORDER] << endl;
 cout << " \t max number of solutions " << options[Option::MAX_NUMBER_OF_SOLUTIONS] << endl;
+cout << " \t max penalty " << options[Option::MAX_PENALTY] << endl;
 cout << " \t compute also non-optimal solutions: " << options.getBool(Option::COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS) << endl;
 cout << " \t compute also connectivity: " << options.getBool(Option::COMPUTE_ALSO_CONNECTIVITY) << endl;
 cout << " \t connectivity cutoff: " << options[Option::CONNECTIVITY_CUTOFF] << endl;
@@ -354,7 +365,7 @@ cout << endl;
 			// check the options
 
 			// What kind of composite do we have?
-			// Do we have a molecule? (Nothingelse is allowed)
+			// Do we have a molecule? (Nothing else is allowed)
 			if (RTTI::isKindOf<Molecule>(ac))
 			{
 				// Store the AtomContainer
@@ -517,7 +528,7 @@ cout << "preassignPenaltyClasses_:" << preassignPenaltyClasses_() << " precomput
 #ifdef BALL_HAS_LPSOLVE
 						strategy = strategies_["ILP"].get();
 #else
-						Log.error() << "Error: BALL was configured without lpsolve support! Try A_STAR instead!" <<
+						Log.error() << "Error: BALL was configured without lpsolve support! Try A_STAR or FPT instead!" <<
 												__FILE__ << " " << __LINE__<< std::endl;
 
 						return Processor::ABORT;
@@ -552,16 +563,20 @@ cout << "preassignPenaltyClasses_:" << preassignPenaltyClasses_() << " precomput
 						solutions_.push_back(*solution);
 
 						// Do we have to find more solutions?
-						bool found_another = true;
-						bool last_sol_is_optimal = true;
+						bool found_another        = true;
+						bool last_sol_is_optimal  = true;
+						bool consider_max_penalty = (max_penalty_ > -1);
+						double last_penalty       = getTotalPenalty(0);
 
 						while (    found_another
 						        && ((getNumberOfComputedSolutions() < max_number_of_solutions_) || (!max_number_of_solutions_))
 						        && (last_sol_is_optimal || (compute_also_non_optimal_solutions_))
+										&& ((last_penalty <= max_penalty_) || (!consider_max_penalty))
 						      )
 						{
 							found_another = computeNextSolution(options.getBool(Option::APPLY_FIRST_SOLUTION));
 							last_sol_is_optimal &= (fabs(getTotalPenalty(0) - getTotalPenalty(solutions_.size()-1)) < 1.e-4);
+							last_penalty  = getTotalPenalty(solutions_.size()-1);
 						}
 #if defined DEBUG_TIMER		
 					timer_.stop();
@@ -1206,6 +1221,9 @@ cout << " ~~~~~~~~ added hydrogen dump ~~~~~~~~~~~~~~~~" << endl;
 		options.setDefaultInteger(AssignBondOrderProcessor::Option::MAX_NUMBER_OF_SOLUTIONS,
 		                          AssignBondOrderProcessor::Default::MAX_NUMBER_OF_SOLUTIONS);
 
+		options.setDefaultInteger(AssignBondOrderProcessor::Option::MAX_PENALTY,
+		                          AssignBondOrderProcessor::Default::MAX_PENALTY);
+
 		options.setDefaultBool(AssignBondOrderProcessor::Option::COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS,
 		                       AssignBondOrderProcessor::Default::COMPUTE_ALSO_NON_OPTIMAL_SOLUTIONS);