EMA

clasp/cli/clasp_cli_options.inl

@@ -1,5 +1,5 @@
 //
-// Copyright (c) 2013-2017 Benjamin Kaufmann
+// Copyright (c) 2013-present Benjamin Kaufmann
 //
 // This file is part of Clasp. See http://www.cs.uni-potsdam.de/clasp/
 //
@@ -374,6 +374,7 @@ OPTION(block_restarts  , ""   , ARG(arg("<arg>")), "Use glucose-style blocking r
        return (arg.off() || (arg>>n>>opt(R = 1.4)>>opt(x = 10000) && n && R >= 1.0 && R <= 5.0))\
          && SET(SELF.blockWindow, n) && SET(SELF.blockScale, (float)R) && SET_OR_FILL(SELF.blockFirst, x);},\
        GET_IF(SELF.blockWindow, SELF.blockWindow, SELF.blockScale, SELF.blockFirst))
+OPTION(block_smooth    , "!,@3" , ARG(flag()), "Use smooth ema strategy for blocking", STORE_FLAG(SELF.blockSmooth), GET(SELF.blockSmooth))
 OPTION(shuffle         , "!"  , ARG(arg("<n1>,<n2>")), "Shuffle problem after <n1>+(<n2>*i) restarts\n", FUN(arg) { uint32 n1 = 0; uint32 n2 = 0;\
        return (arg.off() || (arg>>n1>>opt(n2) && n1)) && SET_OR_FILL(SELF.shuffle, n1) && SET_OR_FILL(SELF.shuffleNext, n2);},\
        GET_IF(SELF.shuffle, SELF.shuffle, SELF.shuffleNext))

clasp/solver_strategies.h

@@ -1,5 +1,5 @@
 //
-// Copyright (c) 2006-2017 Benjamin Kaufmann
+// Copyright (c) 2006-present Benjamin Kaufmann
 //
 // This file is part of Clasp. See http://www.cs.uni-potsdam.de/clasp/
 //
@@ -300,6 +300,7 @@ struct RestartParams {
 	uint32 shuffleNext:14;   /**< Re-Shuffle program every shuffleNext restarts (0: disable). */
 	uint32 upRestart  : 2;   /**< How to update restart sequence after a model was found (one of SeqUpdate). */
 	uint32 cntLocal   : 1;   /**< Count conflicts globally or relative to current branch? */
+	uint32 blockSmooth: 1;   /**< Use simple (0) or smoothing (1) ema for blocking restarts. */
 };
 
 //! Reduce strategy used during solving.

clasp/solver_types.h

@@ -77,7 +77,7 @@ struct DynamicLimit {
 	enum Type   { lbd_limit, level_limit };
 	enum QStrat { keep_never = 0, keep_restart = 1, keep_block = 2, keep_always = 3 };
 	//! Creates new limit with bounded queue of the given window size.
-	static DynamicLimit* create(uint32 window);
+	static DynamicLimit* create(uint32 window, bool smooth, bool gloEma);
 	//! Destroys this object and its bounded queue.
 	void   destroy();
 	//! Resets moving average and adjust limit.
@@ -99,10 +99,8 @@ struct DynamicLimit {
 	uint32 restart(uint32 maxLBD, float k);
 	//! Returns the number of updates since last restart.
 	uint32 runLen()  const { return num_;  }
-	//! Returns the maximal size of the bounded queue.
-	uint32 window()  const { return cap_; }
 	//! Returns whether it is time to restart.
-	bool   reached() const { return runLen() >= window() && (sma(adjust.type) * adjust.rk) > global.avg(adjust.type); }
+	bool   reached() const { return runLen() >= emaLbd_.initLim() && (ema(adjust.type) * adjust.rk) > globalAverage(); }
 	struct {
 		//! Returns the global lbd or conflict level average.
 		double avg(Type t) const { return ratio(sum[t], samples); }
@@ -118,20 +116,31 @@ struct DynamicLimit {
 		float  rk;      //!< LBD/CFL dynamic limit factor (typically < 1.0).
 		Type   type;    //!< Dynamic limit based on lbd or conflict level.
 	} adjust; //!< Data for dynamically adjusting margin ratio (rk).
+
+	double globalAverage() const {
+		return globalAverage(adjust.type);
+	}
+
+	double globalAverage(Type t) const {
+		if (!gloEma_)
+			return global.avg(adjust.type);
+		else
+			return (adjust.type == lbd_limit ? gloLbd_ : gloCfl_).get();
+	}
 private:
-	DynamicLimit(uint32 size);
+	DynamicLimit(uint32 size, bool smooth, bool globalEma);
 	DynamicLimit(const DynamicLimit&);
 	DynamicLimit& operator=(const DynamicLimit&);
-	double sma(Type t)  const { return sum_[t] / double(cap_); }
 	void resetRun(bool keepQ);
-	uint64 sum_[2];
-	uint32 cap_;
-	uint32 pos_;
+
+	double ema(Type t)  const { return (&emaLbd_ + t)->get(); }
+	EmaQ   emaLbd_;
+	EmaQ   emaCfl_;
+	EmaQ   gloLbd_;
+	EmaQ   gloCfl_;
 	uint32 num_;
 	QStrat qStrat_;
-POTASSCO_WARNING_BEGIN_RELAXED
-	uint32 buffer_[0];
-POTASSCO_WARNING_END_RELAXED
+	bool   gloEma_;
 };
 
 //! Type for implementing Glucose-style blocking of restarts.
@@ -140,21 +149,17 @@ POTASSCO_WARNING_END_RELAXED
  * \see A. Biere, A. Froehlich "Evaluating CDCL Restart Schemes"
  */
 struct BlockLimit {
-	explicit BlockLimit(uint32 windowSize, double R = 1.4);
+	explicit BlockLimit(uint32 windowSize, double R = 1.4, bool smooth = false);
 	bool push(uint32 nAssign) {
-		ema = n >= span
-			? exponentialMovingAverage(ema, nAssign, alpha)
-			: cumulativeMovingAverage(ema, nAssign, n);
+		emaQ.push(nAssign);
 		return ++n >= next;
 	}
 	//! Returns the exponential moving average scaled by r.
-	double scaled() const { return ema * r; }
-	double ema;   //!< Current exponential moving average.
-	double alpha; //!< Smoothing factor: 2/(span+1).
+	double scaled() const { return emaQ.get() * r; }
+	EmaQ   emaQ;  //!<
 	uint64 next;  //!< Enable once n >= next.
 	uint64 inc;   //!< Block restart for next inc conflicts.
 	uint64 n;     //!< Number of data points seen so far.
-	uint32 span;  //!< Minimum observation window.
 	float  r;     //!< Scale factor for ema.
 };
 ///////////////////////////////////////////////////////////////////////////////
@@ -183,8 +188,9 @@ struct CoreStats {
 	STAT(uint64 choices;    DOXY(number of choices)     , "choices"           , VALUE(choices)    , LHS.choices    += RHS.choices   )\
 	STAT(uint64 conflicts;  DOXY(number of conflicts)   , "conflicts"         , VALUE(conflicts)  , LHS.conflicts  += RHS.conflicts )\
 	STAT(uint64 analyzed;   DOXY(number of conflicts analyzed), "conflicts_analyzed", VALUE(analyzed)   , LHS.analyzed   += RHS.analyzed  )\
-	STAT(uint64 restarts;   DOXY(number of restarts)    , "restarts"          , VALUE(restarts)   , LHS.restarts   += RHS.restarts  )\
-	STAT(uint64 lastRestart;DOXY(length of last restart), "restarts_last"     , VALUE(lastRestart), LHS.lastRestart = std::max(LHS.lastRestart, RHS.lastRestart))
+	STAT(uint64 restarts;   DOXY(number of restarts)        , "restarts"        , VALUE(restarts)   , LHS.restarts   += RHS.restarts  )\
+	STAT(uint64 lastRestart;DOXY(length of last restart)    , "restarts_last"   , VALUE(lastRestart), LHS.lastRestart = std::max(LHS.lastRestart, RHS.lastRestart)) \
+	STAT(uint64 blRestarts; DOXY(number of blocked restarts), "restarts_blocked", VALUE(blRestarts), LHS.blRestarts = std::max(LHS.blRestarts, RHS.blRestarts))
 
 	CoreStats() { reset(); }
 	void reset();
@@ -300,7 +306,7 @@ struct SolverStats : public CoreStats {
 	~SolverStats();
 	bool enableExtended();
 	bool enable(const SolverStats& o) { return !o.extra || enableExtended(); }
-	void enableLimit(uint32 size);
+	void enableLimit(uint32 size, bool smooth, bool gloEma);
 	void reset();
 	void accu(const SolverStats& o);
 	void accu(const SolverStats& o, bool enableRhs);

clasp/util/misc_types.h

@@ -1,5 +1,5 @@
 //
-// Copyright (c) 2006-2017 Benjamin Kaufmann
+// Copyright (c) 2006-present Benjamin Kaufmann
 //
 // This file is part of Clasp. See http://www.cs.uni-potsdam.de/clasp/
 //
@@ -141,16 +141,67 @@ class RNG {
 /*!
  * Computes ema = currentEma + ((double(sample) - currentEma)*alpha);
  */
-template <class T>
-inline double exponentialMovingAverage(double currentEma, T sample, double alpha) {
-	return (static_cast<double>(sample) * alpha) + (currentEma * (1.0 - alpha));
+inline double exponentialMovingAverage(double currentEma, double sample, double alpha) {
+	return currentEma + (alpha * (sample - currentEma));
 }
 //! Updates the given moving average with the given sample.
 template <class T>
 inline double cumulativeMovingAverage(double currentAvg, T sample, uint64 numSeen) {
 	return (static_cast<double>(sample) + (currentAvg * numSeen)) / static_cast<double>(numSeen + 1);
 }
 
+class EmaQ {
+public:
+	EmaQ(uint32 window, bool initSmooth) : ema_(0.0), alpha_(), initN_(0), initLim_(window), initCma_(!initSmooth) {
+		assert(window > 0);
+		if (initCma_) {
+			alpha_ = 2.0 / static_cast<double>(window + 1);
+		}
+		else {
+			uint32 x = log2(window);
+			alpha_   = 1.0 / static_cast<double>(1u << x);
+		}
+	}
+
+	bool push(double val) {
+		if (initN_ == UINT32_MAX) {
+			ema_ = exponentialMovingAverage(ema_, val, alpha_);
+			return true;
+		}
+		else {
+			if (initCma_) {
+				ema_ = cumulativeMovingAverage(ema_, val, initN_);
+			}
+			else if (initN_ < 32) {
+				double a = std::max(alpha_, 1.0 / static_cast<double>(uint32(1) << initN_));
+				ema_ = exponentialMovingAverage(ema_, val, a);
+			}
+			else {
+				ema_ = exponentialMovingAverage(ema_, val, alpha_);
+			}
+			if (++initN_ == initLim_)
+				initN_ = UINT32_MAX;
+			return initN_ == UINT32_MAX;
+		}
+	}
+
+	double get()     const { return ema_; }
+	bool   valid()   const { return initN_ == UINT32_MAX; }
+	uint32 initLim() const { return initLim_; }
+
+	void clear() {
+		ema_   = 0;
+		initN_ = 0;
+	}
+
+private:
+	double ema_;        //!< Current exponential moving average.
+	double alpha_;      //!< Smoothing factor.
+	uint32 initN_;      //!< Number of data points (capped at initLim).
+	uint32 initLim_:31; //!< Number of data points needed for initialization.
+	uint32 initCma_: 1; //!< Init scheme (adjusted smoothing = 0, CMA = 1).
+};
+
 //! An unary operator function that calls p->destroy().
 struct DestroyObject {
 	template <class T> void operator()(T* p) const { if (p) p->destroy(); }

src/clasp_options.cpp

@@ -270,11 +270,33 @@ static int xconvert(const char* x, ScheduleStrategy& out, const char** errPos, i
 		std::pair<double, uint32> arg(0, 0);
 		if (*next != ',' || !xconvert(next+1, arg, &next, e)) { return 0; }
 		if (strncasecmp(x, "d", 1) == 0 && arg.first > 0.0) {
-			SET_OR_FILL(arg.second, arg.second << 2);
-			uint32 keep = 0;
-			if (*next == ',' && (!xconvert(next+1, keep, &next, e) || keep > 3)) { return 0; }
-			arg.second &= ~3u;
-			arg.second |= keep;
+			SET_OR_FILL(arg.second, arg.second << 4);
+			arg.second &= ~15u;
+			uint32 extra = 0;
+			if (*next == ',') {
+				if (!*++next)  { return 0; };
+				while (*next && *next != ',') {
+					uint32 temp = extra;
+					switch (tolower(static_cast<unsigned char>(*next))) {
+						case 'r':
+							if (store_set_bit(temp, 0) == extra) { return 0; }
+							break;
+						case 'b':
+							if (store_set_bit(temp, 1) == extra) { return 0; }
+							break;
+						case 's':
+							if (store_set_bit(temp, 2) == extra) { return 0; }
+							break;
+						case 'e':
+							if (store_set_bit(temp, 3) == extra) { return 0; }
+							break;
+						default: return 0;
+					}
+					extra = temp;
+					++next;
+				}
+			}
+			arg.second |= extra;
 			out = ScheduleStrategy(ScheduleStrategy::User, base, arg.first, arg.second);
 		}
 		else if (strncasecmp(x, "d", 1) != 0 && arg.first >= 1.0){ out = ScheduleStrategy::geom(base, arg.first, arg.second); }
@@ -305,7 +327,14 @@ static std::string& xconvert(std::string& out, const ScheduleStrategy& sched) {
 		case ScheduleStrategy::User:
 			out[t] = 'd';
 			xconvert(out.append(1, ','), (double)sched.grow);
-			return xconvert(out.append(1, ','), std::make_pair(sched.len >> 2, sched.len & 3u));
+			xconvert(out.append(1, ','), sched.len >> 4);
+			if ((sched.len & 15u) != 0) {
+				out.append(1, ',');
+				if (test_bit(sched.len, 0)) out.append(1, 'r');
+				if (test_bit(sched.len, 1)) out.append(1, 'b');
+				if (test_bit(sched.len, 2)) out.append(1, 's');
+				if (test_bit(sched.len, 3)) out.append(1, 'e');
+			}
 		default: POTASSCO_ASSERT(false, "xconvert(ScheduleStrategy): unknown type");
 	}
 }

src/clasp_output.cpp

@@ -1,5 +1,5 @@
 //
-// Copyright (c) 2009-2017 Benjamin Kaufmann
+// Copyright (c) 2009-present Benjamin Kaufmann
 //
 // This file is part of Clasp. See http://www.cs.uni-potsdam.de/clasp/
 //
@@ -1188,7 +1188,7 @@ void TextOutput::printStats(const Clasp::SolverStats& st) const {
 	printf(" (Analyzed: %" PRIu64")\n", st.backjumps());
 	printKeyValue("Restarts", "%-8" PRIu64"", st.restarts);
 	if (st.restarts) {
-		printf(" (Average: %.2f Last: %" PRIu64")", st.avgRestart(), st.lastRestart);
+		printf(" (Average: %.2f Last: %" PRIu64" Blocked: %" PRIu64")", st.avgRestart(), st.lastRestart, st.blRestarts);
 	}
 	printf("\n");
 	if (!st.extra) return;

src/solve_algorithms.cpp

@@ -116,11 +116,12 @@ BasicSolve::State::State(Solver& s, const SolveParams& p) {
 		dbRedInit = 0;
 	}
 	if (p.restart.dynamic()) {
-		s.stats.enableLimit(p.restart.sched.base);
-		s.stats.limit->reset();
+		bool smooth = test_bit(p.restart.sched.len, 2);
+		bool ema    = test_bit(p.restart.sched.len, 3);
+		s.stats.enableLimit(p.restart.sched.base, smooth, ema);
 	}
 	if (p.restart.blockScale > 0.0f && p.restart.blockWindow > 0) {
-		rsBlock.reset(new BlockLimit(p.restart.blockWindow, p.restart.blockScale));
+		rsBlock.reset(new BlockLimit(p.restart.blockWindow, p.restart.blockScale, p.restart.blockSmooth != 0));
 		rsBlock->inc  = std::max(p.restart.sched.base, uint32(50));
 		rsBlock->next = std::max(p.restart.blockWindow, p.restart.blockFirst);
 	}
@@ -205,7 +206,7 @@ ValueRep BasicSolve::State::solve(Solver& s, const SolveParams& p, SolveLimits*
 			}
 			if (sLimit.restart.dynamic) {
 				n = sLimit.restart.dynamic->runLen();
-				uint32 limLbd = rs.len >> 2;
+				uint32 limLbd = rs.len >> 4;
 				sLimit.restart.conflicts = sLimit.restart.dynamic->restart(limLbd ? limLbd : UINT32_MAX, (float)rs.grow);
 			}
 			else {

src/solver.cpp

@@ -1,5 +1,5 @@
 //
-// Copyright (c) 2006-2017 Benjamin Kaufmann
+// Copyright (c) 2006-present Benjamin Kaufmann
 //
 // This file is part of Clasp. See http://www.cs.uni-potsdam.de/clasp/
 //
@@ -1833,6 +1833,7 @@ ValueRep Solver::search(SearchLimits& limit, double rf) {
 						if (limit.restart.dynamic) { limit.restart.dynamic->resetBlock(); }
 						else                       { limit.restart.conflicts += block->inc; }
 						block->next = block->n + block->inc;
+						++stats.blRestarts;
 					}
 				} while (resolveConflict() && !propagate() && (++n, true));
 				limit.used += n;