diff --git a/ortools/com/google/ortools/sat/CpModel.java b/ortools/com/google/ortools/sat/CpModel.java
index ef68017da3..335c58b01a 100644
--- a/ortools/com/google/ortools/sat/CpModel.java
+++ b/ortools/com/google/ortools/sat/CpModel.java
@@ -62,6 +62,40 @@ public class CpModel {
     modelBuilder = CpModelProto.newBuilder();
   }
 
+  // Domains
+
+  /** Create a flattened list of intervals from a list of integer values.*/
+  public long[] domainFromValues(long[] values) {
+    return SatHelper.domainFromValues(values);
+  }
+
+  /** Create a flattened list of intervals from a list of integer values.*/
+  public long[] domainFromValues(int[] values) {
+    return SatHelper.domainFromValues(toLongArray(values));
+  }
+
+  /** Consolidate a flattened list of intervals. */
+  public long[] domainFromIntervals(long[] bounds) {
+    long[] starts = new long[bounds.length / 2];
+    long[] ends = new long[bounds.length / 2];
+    for (int i = 0; i < bounds.length / 2; ++i) {
+      starts[i] = bounds[2 * i];
+      ends[i] = bounds[2 * i + 1];
+    }
+    return SatHelper.domainFromStartsAndEnds(starts, ends);
+  }
+
+  /** Consolidate a flattened list of intervals. */
+  public long[] domainFromIntervals(int[] bounds) {
+    long[] starts = new long[bounds.length / 2];
+    long[] ends = new long[bounds.length / 2];
+    for (int i = 0; i < bounds.length / 2; ++i) {
+      starts[i] = bounds[2 * i];
+      ends[i] = bounds[2 * i + 1];
+    }
+    return SatHelper.domainFromStartsAndEnds(starts, ends);
+  }
+
   // Integer variables.
 
   /** Creates an integer variable with domain [lb, ub]. */
@@ -70,16 +104,25 @@ public class CpModel {
   }
 
   /**
-   * Creates an integer variable with an enumerated domain.
+   * Creates an integer variable with an array of values.
    *
-   * @param bounds a flattened list of disjoint intervals
+   * @param values a list of integer values
    * @param name the name of the variable
-   * @return a variable whose domain is the union of [bounds[2 * i] .. bounds[2 * i + 1]]
-   *     <p>To create a variable with enumerated domain [1, 2, 3, 5, 7, 8], pass in the array [1, 3,
-   *     5, 5, 7, 8].
+   * @return a variable whose domain is the set of values passed as argument.
    */
-  public IntVar newEnumeratedIntVar(long[] bounds, String name) {
-    return new IntVar(modelBuilder, bounds, name);
+  public IntVar newIntVarFromValues(long[] values, String name) {
+    return new IntVar(modelBuilder, domainFromValues(values), name);
+  }
+
+  /**
+   * Creates an integer variable with an array of values.
+   *
+   * @param values a list of integer values
+   * @param name the name of the variable
+   * @return a variable whose domain is the set of values passed as argument.
+   */
+  public IntVar newIntVarFromValues(int[] values, String name) {
+    return new IntVar(modelBuilder, domainFromValues(values), name);
   }
 
   /**
@@ -91,8 +134,21 @@ public class CpModel {
    *     <p>To create a variable with enumerated domain [1, 2, 3, 5, 7, 8], pass in the array [1, 3,
    *     5, 5, 7, 8].
    */
-  public IntVar newEnumeratedIntVar(int[] bounds, String name) {
-    return new IntVar(modelBuilder, toLongArray(bounds), name);
+  public IntVar newIntVarFromIntervals(long[] bounds, String name) {
+    return new IntVar(modelBuilder, domainFromIntervals(bounds), name);
+  }
+
+  /**
+   * Creates an integer variable with an enumerated domain.
+   *
+   * @param bounds a flattened list of disjoint intervals
+   * @param name the name of the variable
+   * @return a variable whose domain is the union of [bounds[2 * i] .. bounds[2 * i + 1]]
+   *     <p>To create a variable with enumerated domain [1, 2, 3, 5, 7, 8], pass in the array [1, 3,
+   *     5, 5, 7, 8].
+   */
+  public IntVar newIntVarFromIntervals(int[] bounds, String name) {
+    return new IntVar(modelBuilder, domainFromIntervals(bounds), name);
   }
 
   /** Creates a Boolean variable with the given name. */
diff --git a/ortools/sat/csharp/CpModel.cs b/ortools/sat/csharp/CpModel.cs
index a727023d6c..d2a45b6fff 100644
--- a/ortools/sat/csharp/CpModel.cs
+++ b/ortools/sat/csharp/CpModel.cs
@@ -15,6 +15,7 @@ namespace Google.OrTools.Sat
 {
 using System;
 using System.Collections.Generic;
+using System.Linq;
 
 /// <summary>
 ///  Wrapper class around the cp_model proto.
@@ -39,6 +40,22 @@ public class CpModel
     return -index - 1;
   }
 
+  // Domains
+
+  public long[] DomainFromValues(long[] values) {
+    return SatHelper.DomainFromValues(values);
+  }
+
+  public long[] DomainFromIntervals(long[] bounds) {
+    long[] starts = new long[bounds.Length / 2];
+    long[] ends = new long[bounds.Length / 2];
+    for (int i = 0; i < bounds.Length / 2; ++i) {
+      starts[i] = bounds[2 * i];
+      ends[i] = bounds[2 * i + 1];
+    }
+    return SatHelper.DomainFromStartsAndEnds(starts, ends);
+  }
+
   // Integer variables and constraints.
 
   public IntVar NewIntVar(long lb, long ub, string name)
@@ -47,9 +64,14 @@ public class CpModel
     return new IntVar(model_, bounds, name);
   }
 
-  public IntVar NewEnumeratedIntVar(IEnumerable<long> bounds, string name)
+  public IntVar NewIntVarFromValues(IEnumerable<long> values, string name)
   {
-    return new IntVar(model_, bounds, name);
+    return new IntVar(model_, DomainFromValues(values.ToArray()), name);
+  }
+
+  public IntVar NewIntVarFromIntervals(IEnumerable<long> bounds, string name)
+  {
+    return new IntVar(model_, DomainFromIntervals(bounds.ToArray()), name);
   }
 
   // Constants (named or not).
@@ -129,6 +151,24 @@ public class CpModel
     return ct;
   }
 
+  public Constraint AddLinearSumWithBounds(
+      IEnumerable<IntVar> variables, IEnumerable<long> bounds)
+  {
+    Constraint ct = new Constraint(model_);
+    LinearConstraintProto lin = new LinearConstraintProto();
+    foreach (IntVar var in variables)
+    {
+      lin.Vars.Add(var.Index);
+      lin.Coeffs.Add(1);
+    }
+    foreach (long b in bounds)
+    {
+      lin.Domain.Add(b);
+    }
+    ct.Proto.Linear = lin;
+    return ct;
+  }
+
   public Constraint AddLinearConstraintWithBounds(
       IEnumerable<Tuple<IntVar, long>> terms, IEnumerable<long> bounds)
   {
diff --git a/ortools/sat/csharp/sat.i b/ortools/sat/csharp/sat.i
index 766b6c76d6..87253166e6 100644
--- a/ortools/sat/csharp/sat.i
+++ b/ortools/sat/csharp/sat.i
@@ -18,9 +18,11 @@ using System.Collections;
 %}
 
 %include "stdint.i"
+%include "std_vector.i"
 
 %include "ortools/base/base.i"
 %include "ortools/util/csharp/proto.i"
+%include "ortools/util/csharp/vector.i"
 
 %{
 #include "ortools/sat/cp_model.pb.h"
@@ -48,6 +50,9 @@ PROTO_INPUT(operations_research::sat::CpSolverResponse,
 PROTO2_RETURN(operations_research::sat::CpSolverResponse,
               Google.OrTools.Sat.CpSolverResponse);
 
+%template(SatInt64Vector) std::vector<int64>;
+VECTOR_AS_CSHARP_ARRAY(int64, int64, long, SatInt64Vector);
+
 %ignoreall
 
 // SatParameters are proto2, thus not compatible with C# Protobufs.
@@ -62,6 +67,8 @@ PROTO2_RETURN(operations_research::sat::CpSolverResponse,
 %unignore operations_research::sat::SatHelper::ModelStats;
 %unignore operations_research::sat::SatHelper::SolverResponseStats;
 %unignore operations_research::sat::SatHelper::ValidateModel;
+%unignore operations_research::sat::SatHelper::DomainFromValues;
+%unignore operations_research::sat::SatHelper::DomainFromStartsAndEnds;
 
 %feature("director") operations_research::sat::SolutionCallback;
 %unignore operations_research::sat::SolutionCallback;
diff --git a/ortools/sat/doc/integer_arithmetic.md b/ortools/sat/doc/integer_arithmetic.md
index 47f66080c8..790fadf5e1 100644
--- a/ortools/sat/doc/integer_arithmetic.md
+++ b/ortools/sat/doc/integer_arithmetic.md
@@ -22,16 +22,19 @@ In Java, Python, and C#:
     value as in [5, 5].
 -   To create a variable with a single value domain, use the `NewConstant()` API
     (or `newConstant()` in Java).
--   To represent an enumerated list of values, for example {-5, -4, -3, 1, 3, 4,
-    6, 6}, you need to rewrite it as a list of intervals [-5, -3] U [1] U [3,
-    6], then flatten the list into a single list of integers. This gives `[-5,
-    -3, 1, 1, 3, 6]` in python, or `new long[] {-5, -3, 1, 1, 3, 6}` in Java or
-    C#.
 -   To create a variable with an enumerated domain, use the
-    `NewEnumeratedIntVar()` API as in:
-    -   Python: `model.NewEnumeratedIntVar([-5, -3, 1, 1, 3, 6], 'x')`
-    -   Java: `model.newEnumeratedIntVar(new long[] {-5, -3, 1, 1, 3, 6}, "x")`
-    -   C#: `model.NewEnumeratedIntVar(new long[] {-5, -3, 1, 1, 3, 6}, "x")`
+    `NewIntVarFromValues()` API as in:
+    -   Python: `model.NewIntVarFromValues([-5, -4, -3, 1, 3, 4], 'x')`
+    -   Java: `model.newIntVarFromValues(new long[] {-5, -4, -3, 1, 3, 4}, "x")`
+    -   C#: `model.NewIntVarFromValues(new long[] {-5, -4, -3, 1, 3, 4}, "x")`
+-   You can also create an enumerated domain by passing intervals instead of a
+    list of values by using the `NewIntVarFromIntervals` method. Note that the
+    Java and C# versions use a single one dimension integer array for
+    convenience. The following code samples build integer variables with domain
+    {-5, -4, -3, 1, 2, 3, .., 99, 100}:
+    -   Python: `model.NewIntVarFromIntervals([[-5, -3], [1, 100]], 'x')`
+    -   Java: `model.newIntVarFromIntervals(new long[] {-5, -3, 1, 100}, "x")`
+    -   C#: `model.NewIntVarFromIntervals(new long[] {-5, -3, 1, 100}, "x")`
 -   To exclude a single value, use int64min and int64max values as in [int64min,
     4, 6, int64max]:
     -   Python: `cp_model.INT_MIN` and `cp_model.INT_MAX`
@@ -52,23 +55,38 @@ In C++, domains use the Domain class.
     1}, {3, 6}})).WithName("x")`.
 -   To exclude a single value, use `Domain(5).Complement()`.
 
+
 ## Linear constraints
 
 In **C++** and **Java**, the model supports linear constraints as in:
 
     x <= y + 3 (also ==, !=, <, >=, >).
 
-as well as domain constraints as in:
+These are available through specific methods of the cp_model like
+`cp_model.AddEquality(x, 3)` in C++, `cp_model.addGreaterThan(x, 10)` in java.
+
+**Python** and **C\#** CP-SAT APIs support general linear arithmetic operators
+(+, *, -, ==, >=, >, <, <=, !=). You need to use the Add method of the cp_model
+as in `cp_model.Add(x != 3)`.
+
+In all language, you can use more advanced domains in linear constraints of the
+form:
 
     sum(ai * xi) in domain
+    sum(xi) in domain
 
-where domain uses the same encoding as integer variables. These are available
-through specific methods of the cp_model like `cp_model.AddEquality(x, 3)` in
-C++, `cp_model.addGreaterThan(x, 10)` in java.
+In Python, Java, and C#, the model class contains two methods to build these
+domains, similar from the IntVarFromValues and IntVarFromIntervals.
+
+-   Python: `model.DomainFromValues([-5, -4, -3, 1, 3, 4])`
+-   Java: `model.domainFromValues(new long[] {-5, -4, -3, 1, 3, 4})`
+-   C#: `model.DomainFromValues(new long[] {-5, -4, -3, 1, 3, 4})`
+-   Python: `model.DomainrFromIntervals([[-5, -3], [1, 100]])`
+-   Java: `model.domainFromIntervals(new long[] {-5, -3, 1, 100})`
+-   C#: `model.DomainFromIntervals(new long[] {-5, -3, 1, 100})`
+
+In C++, you can reuse the Domain class.
 
-**Python** and **C\#** CP-SAT APIs support general linear arithmetic (+, *, -,
-==, >=, >, <, <=, !=). You need to use the Add method of the cp_model as in
-`cp_model.Add(x != 3)`.
 
 ## Rabbits and Pheasants examples
 
@@ -665,13 +683,16 @@ def step_function_sample_sat():
 
   # expr == 0 on [5, 6] U [8, 10]
   b0 = model.NewBoolVar('b0')
-  model.AddLinearConstraintWithBounds([(x, 1)], [5, 6, 8, 10]).OnlyEnforceIf(b0)
+  model.AddSumConstraintWithBounds(
+      [x], model.DomainFromIntervals([(5, 6), (8, 10)])).OnlyEnforceIf(b0)
   model.Add(expr == 0).OnlyEnforceIf(b0)
 
   # expr == 2 on [0, 1] U [3, 4] U [11, 20]
   b2 = model.NewBoolVar('b2')
-  model.AddLinearConstraintWithBounds([(x, 1)],
-                                      [0, 1, 3, 4, 11, 20]).OnlyEnforceIf(b2)
+  model.AddSumConstraintWithBounds([x],
+                                   model.DomainFromIntervals([
+                                       (0, 1), (3, 4), (11, 20)
+                                   ])).OnlyEnforceIf(b2)
   model.Add(expr == 2).OnlyEnforceIf(b2)
 
   # expr == 3 when x == 7
diff --git a/ortools/sat/doc/scheduling.md b/ortools/sat/doc/scheduling.md
index fa5eae812c..493c6a657b 100644
--- a/ortools/sat/doc/scheduling.md
+++ b/ortools/sat/doc/scheduling.md
@@ -1187,11 +1187,14 @@ public class RankingSampleSat
 
 ## Intervals spanning over breaks in the calendar
 
-Sometimes, a task can be interrupted by a break (overnight, lunch break). In that context, although the processing time of the task is the same, the duration can vary.
+Sometimes, a task can be interrupted by a break (overnight, lunch break). In
+that context, although the processing time of the task is the same, the duration
+can vary.
 
-To implement this feature, we will have the duration of the task be a function of the start of the task. This is implemented using channeling constraints.
+To implement this feature, we will have the duration of the task be a function
+of the start of the task. This is implemented using channeling constraints.
 
-The following code outputs:
+The following code displays:
 
     start=8 duration=3 across=0
     start=9 duration=3 across=0
@@ -1238,22 +1241,22 @@ def SchedulingWithCalendarSampleSat():
   # The data is the following:
   #   Work starts at 8h, ends at 18h, with a lunch break between 13h and 14h.
   #   We need to schedule a task that needs 3 hours of processing time.
-  #   Total duration can be 3 or 4 (if it span across the lunch break.
+  #   Total duration can be 3 or 4 (if it spans the lunch break).
   #
-  # Because the duration is at least 3, work cannot start after 15h.
+  # Because the duration is at least 3 hours, work cannot start after 15h.
   # Because of the break, work cannot start at 13h.
 
-  start = model.NewEnumeratedIntVar([8, 12, 14, 15], 'start')
+  start = model.NewIntVarFromIntervals([(8, 12), (14, 15)], 'start')
   duration = model.NewIntVar(3, 4, 'duration')
   end = model.NewIntVar(8, 18, 'end')
   unused_interval = model.NewIntervalVar(start, duration, end, 'interval')
 
   # We have 2 states (spanning across lunch or not)
   across = model.NewBoolVar('across')
-  model.AddLinearConstraintWithBounds(
-      [(start, 1)], [8, 10, 14, 15]).OnlyEnforceIf(across.Not())
-  model.AddLinearConstraintWithBounds([(start, 1)],
-                                      [11, 12]).OnlyEnforceIf(across)
+  non_spanning_hours = model.DomainFromValues([8, 9, 10, 14, 15])
+  model.AddSumConstraintWithBounds([start], non_spanning_hours).OnlyEnforceIf(
+      across.Not())
+  model.AddSumConstraint([start], 11, 12).OnlyEnforceIf(across)
   model.Add(duration == 3).OnlyEnforceIf(across.Not())
   model.Add(duration == 4).OnlyEnforceIf(across)
 
@@ -1267,7 +1270,7 @@ def SchedulingWithCalendarSampleSat():
   # Force the solver to follow the decision strategy exactly.
   solver.parameters.search_branching = cp_model.FIXED_SEARCH
 
-  # Search and print out all solutions.
+  # Search and print all solutions.
   solution_printer = VarArraySolutionPrinter([start, duration, across])
   solver.SearchForAllSolutions(model, solution_printer)
 
@@ -1282,4 +1285,3 @@ SchedulingWithCalendarSampleSat()
 ## Convex hull of a set of intervals
 
 ## Reservoir constraint
-
diff --git a/ortools/sat/java/sat.i b/ortools/sat/java/sat.i
index 42db3e97d2..5d82a6903d 100644
--- a/ortools/sat/java/sat.i
+++ b/ortools/sat/java/sat.i
@@ -16,6 +16,7 @@
 %include "ortools/base/base.i"
 
 %include "ortools/util/java/proto.i"
+%include "ortools/util/java/vector.i"
 
 %{
 #include "ortools/sat/cp_model.pb.h"
@@ -56,6 +57,8 @@ PROTO2_RETURN(operations_research::sat::CpSolverResponse,
 %rename (modelStats) operations_research::sat::SatHelper::ModelStats;
 %rename (solverResponseStats) operations_research::sat::SatHelper::SolverResponseStats;
 %rename (validateModel) operations_research::sat::SatHelper::ValidateModel;
+%rename (domainFromValues) operations_research::sat::SatHelper::DomainFromValues;
+%rename (domainFromStartsAndEnds) operations_research::sat::SatHelper::DomainFromStartsAndEnds;
 
 // We use directors for the solution callback.
 %feature("director") operations_research::sat::SolutionCallback;
diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc
index 24f78cb1cb..b9b0489c8e 100644
--- a/ortools/sat/linear_programming_constraint.cc
+++ b/ortools/sat/linear_programming_constraint.cc
@@ -28,6 +28,7 @@
 #include "ortools/glop/preprocessor.h"
 #include "ortools/glop/status.h"
 #include "ortools/graph/strongly_connected_components.h"
+#include "ortools/lp_data/lp_types.h"
 #include "ortools/util/saturated_arithmetic.h"
 
 namespace operations_research {
@@ -339,6 +340,11 @@ bool LinearProgrammingConstraint::SolveLp() {
     simplex_.ClearStateForNextSolve();
     return false;
   }
+  UpdateDegeneracyData();
+  if (average_degeneracy_.CurrentAverage() >= 1000.0) {
+    VLOG(1) << "High average degeneracy: "
+            << average_degeneracy_.CurrentAverage();
+  }
 
   if (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL) {
     lp_solution_is_set_ = true;
@@ -651,9 +657,11 @@ bool LinearProgrammingConstraint::Propagate() {
   // Put an iteration limit on the work we do in the simplex for this call. Note
   // that because we are "incremental", even if we don't solve it this time we
   // will make progress towards a solve in the lower node of the tree search.
-  //
-  // TODO(user): Put more at the root, and less afterwards?
-  parameters.set_max_number_of_iterations(500);
+  if (trail_->CurrentDecisionLevel() == 0) {
+    parameters.set_max_number_of_iterations(2000);
+  } else {
+    parameters.set_max_number_of_iterations(500);
+  }
   if (sat_parameters_.use_exact_lp_reason()) {
     parameters.set_change_status_to_imprecise(false);
     parameters.set_primal_feasibility_tolerance(1e-7);
@@ -1324,6 +1332,21 @@ void LinearProgrammingConstraint::FillReducedCostsReason() {
   integer_trail_->RemoveLevelZeroBounds(&integer_reason_);
 }
 
+void LinearProgrammingConstraint::UpdateDegeneracyData() {
+  const int num_vars = integer_variables_.size();
+  int num_non_basic_with_zero_rc = 0;
+  for (int i = 0; i < num_vars; i++) {
+    const double rc = simplex_.GetReducedCost(glop::ColIndex(i));
+    if (rc != 0.0) continue;
+    if (simplex_.GetVariableStatus(glop::ColIndex(i)) ==
+        glop::VariableStatus::BASIC) {
+      continue;
+    }
+    num_non_basic_with_zero_rc++;
+  }
+  average_degeneracy_.AddData(num_non_basic_with_zero_rc);
+}
+
 void LinearProgrammingConstraint::FillDualRayReason() {
   integer_reason_.clear();
   const int num_vars = integer_variables_.size();
diff --git a/ortools/sat/linear_programming_constraint.h b/ortools/sat/linear_programming_constraint.h
index a9badfc2ee..9dd1b36ec6 100644
--- a/ortools/sat/linear_programming_constraint.h
+++ b/ortools/sat/linear_programming_constraint.h
@@ -29,6 +29,7 @@
 #include "ortools/sat/linear_constraint.h"
 #include "ortools/sat/linear_constraint_manager.h"
 #include "ortools/sat/model.h"
+#include "ortools/sat/util.h"
 #include "ortools/util/rev.h"
 #include "ortools/util/time_limit.h"
 
@@ -147,6 +148,11 @@ class LinearProgrammingConstraint : public PropagatorInterface,
   // Tie-breaking is done using the variable natural order.
   std::function<LiteralIndex()> LPReducedCostAverageBranching();
 
+  // Average number of nonbasic variables with zero reduced costs.
+  double average_degeneracy() const {
+    return average_degeneracy_.CurrentAverage();
+  }
+
  private:
   // Reinitialize the LP from a potentially new set of constraints.
   // This fills all data structure and properly rescale the underlying LP.
@@ -197,6 +203,10 @@ class LinearProgrammingConstraint : public PropagatorInterface,
   // computations, true otherwise.
   bool FillExactDualRayReason();
 
+  // Computes number of non basic variables with zero reduced costs and updates
+  // 'average_degeneracy_'.
+  void UpdateDegeneracyData();
+
   // From a set of row multipliers (at LP scale), scale them back to the CP
   // world and then make them integer (eventually multiplying them by a new
   // scaling factor returned in *scaling).
@@ -377,6 +387,9 @@ class LinearProgrammingConstraint : public PropagatorInterface,
   std::vector<double> sum_cost_down_;
   std::vector<int> num_cost_up_;
   std::vector<int> num_cost_down_;
+
+  // Defined as average number of nonbasic variables with zero reduced costs.
+  IncrementalAverage average_degeneracy_;
 };
 
 // A class that stores which LP propagator is associated to each variable.
diff --git a/ortools/sat/pseudo_costs.cc b/ortools/sat/pseudo_costs.cc
index 7f70b9ee0c..022b7fffa8 100644
--- a/ortools/sat/pseudo_costs.cc
+++ b/ortools/sat/pseudo_costs.cc
@@ -19,6 +19,7 @@
 #include "ortools/sat/integer.h"
 #include "ortools/sat/sat_decision.h"
 #include "ortools/sat/sat_parameters.pb.h"
+#include "ortools/sat/util.h"
 
 namespace operations_research {
 namespace sat {
@@ -27,15 +28,14 @@ PseudoCosts::PseudoCosts(Model* model)
     : integer_trail_(*model->GetOrCreate<IntegerTrail>()),
       parameters_(*model->GetOrCreate<SatParameters>()) {
   const int num_vars = integer_trail_.NumIntegerVariables().value();
-  pseudo_costs_.resize(num_vars, 0.0);
-  num_recordings_.resize(num_vars, 0);
+  pseudo_costs_.resize(num_vars);
 }
 
 void PseudoCosts::InitializeCosts(double initial_value) {
   if (pseudo_costs_initialized_) return;
   VLOG(1) << "Initializing pseudo costs";
   for (int i = 0; i < pseudo_costs_.size(); ++i) {
-    pseudo_costs_[IntegerVariable(i)] = initial_value;
+    pseudo_costs_[IntegerVariable(i)].Reset(initial_value);
   }
   pseudo_costs_initialized_ = true;
 }
@@ -44,12 +44,10 @@ void PseudoCosts::UpdateCostForVar(IntegerVariable var, double new_cost) {
   if (var >= pseudo_costs_.size()) {
     // Create space for new variable and its negation.
     const int new_size = std::max(var, NegationOf(var)).value() + 1;
-    pseudo_costs_.resize(new_size, initial_cost_);
-    num_recordings_.resize(new_size, 0);
+    pseudo_costs_.resize(new_size, IncrementalAverage(initial_cost_));
   }
   CHECK_LT(var, pseudo_costs_.size());
-  num_recordings_[var]++;
-  pseudo_costs_[var] += (new_cost - pseudo_costs_[var]) / num_recordings_[var];
+  pseudo_costs_[var].AddData(new_cost);
 }
 
 void PseudoCosts::UpdateCost(
@@ -96,15 +94,15 @@ IntegerVariable PseudoCosts::GetBestDecisionVar() {
     const IntegerValue lb = integer_trail_.LowerBound(positive_var);
     const IntegerValue ub = integer_trail_.UpperBound(positive_var);
     if (lb >= ub) continue;
-    if (num_recordings_[positive_var] + num_recordings_[negative_var] <
+    if (GetRecordings(positive_var) + GetRecordings(negative_var) <
         parameters_.pseudo_cost_reliability_threshold()) {
       continue;
     }
 
     // TODO(user): Experiment with different ways to merge the costs.
     const double current_merged_cost =
-        std::min(pseudo_costs_[positive_var], epsilon) *
-        std::min(pseudo_costs_[negative_var], epsilon);
+        std::min(GetCost(positive_var), epsilon) *
+        std::min(GetCost(negative_var), epsilon);
 
     if (current_merged_cost > best_cost) {
       chosen_var = positive_var;
@@ -114,7 +112,7 @@ IntegerVariable PseudoCosts::GetBestDecisionVar() {
 
   // Pick the direction with best pseudo cost.
   if (chosen_var != kNoIntegerVariable &&
-      pseudo_costs_[chosen_var] < pseudo_costs_[NegationOf(chosen_var)]) {
+      GetCost(chosen_var) < GetCost(NegationOf(chosen_var))) {
     chosen_var = NegationOf(chosen_var);
   }
   return chosen_var;
diff --git a/ortools/sat/pseudo_costs.h b/ortools/sat/pseudo_costs.h
index 1e21808308..868b011b61 100644
--- a/ortools/sat/pseudo_costs.h
+++ b/ortools/sat/pseudo_costs.h
@@ -17,6 +17,7 @@
 #include <vector>
 
 #include "ortools/sat/integer.h"
+#include "ortools/sat/util.h"
 
 namespace operations_research {
 namespace sat {
@@ -43,14 +44,14 @@ class PseudoCosts {
   // Returns the pseudo cost of given variable. Currently used for testing only.
   double GetCost(IntegerVariable var) const {
     CHECK_LT(var, pseudo_costs_.size());
-    return pseudo_costs_[var];
+    return pseudo_costs_[var].CurrentAverage();
   }
 
   // Returns the number of recordings of given variable. Currently used for
   // testing only.
   int GetRecordings(IntegerVariable var) const {
-    CHECK_LT(var, num_recordings_.size());
-    return num_recordings_[var];
+    CHECK_LT(var, pseudo_costs_.size());
+    return pseudo_costs_[var].NumRecords();
   }
 
  private:
@@ -70,8 +71,7 @@ class PseudoCosts {
 
   double initial_cost_ = 0.0;
 
-  gtl::ITIVector<IntegerVariable, double> pseudo_costs_;
-  gtl::ITIVector<IntegerVariable, int> num_recordings_;
+  gtl::ITIVector<IntegerVariable, IncrementalAverage> pseudo_costs_;
 };
 
 // Returns extracted information to update pseudo costs from the given
diff --git a/ortools/sat/python/cp_model.py b/ortools/sat/python/cp_model.py
index ba1b940461..c54612624e 100644
--- a/ortools/sat/python/cp_model.py
+++ b/ortools/sat/python/cp_model.py
@@ -10,11 +10,7 @@
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
-"""Propose a natural language on top of cp_model_pb2 python proto.
-
-This file implements a easy-to-use API on top of the cp_model_pb2 protobuf
-defined in ../ .
-"""
+"""Methods for building and solving CP-SAT models."""
 
 from __future__ import absolute_import
 from __future__ import division
@@ -66,6 +62,17 @@ AUTOMATIC_SEARCH = sat_parameters_pb2.SatParameters.AUTOMATIC_SEARCH
 FIXED_SEARCH = sat_parameters_pb2.SatParameters.FIXED_SEARCH
 PORTFOLIO_SEARCH = sat_parameters_pb2.SatParameters.PORTFOLIO_SEARCH
 LP_SEARCH = sat_parameters_pb2.SatParameters.LP_SEARCH
+"""The following sections describe methods for building and solving
+CP-SAT models, and related tasks:
+
+* [Create model](#ortools.sat.python.cp_model.CpModel): Methods for creating
+models, including variables and constraints.
+* [Solve](#ortools.sat.python.cp_model.CpSolver): Methods for solving
+a model and evaluating solutions.
+* [Solution callback](#ortools.sat.python.cp_model.CpSolverSolutionCallback):
+Create a callback that is invoked every time the solver finds a new solution.
+* [Solution printer](#ortools.sat.python.cp_model.ObjectiveSolutionPrinter):
+Print objective values and elapsed time for intermediate solutions."""
 
 
 def DisplayBounds(bounds):
@@ -82,7 +89,7 @@ def DisplayBounds(bounds):
 
 
 def ShortName(model, i):
-    """Returns a short name of an integer variable, or its negation."""
+    """Return a short name of an integer variable, or its negation."""
     if i < 0:
         return 'Not(%s)' % ShortName(model, -i - 1)
     v = model.variables[i]
@@ -339,11 +346,11 @@ class IntVar(LinearExpression):
         self.__negation = None
 
     def Index(self):
-        """Returns the index of the variable in the model."""
+        """Return the index of the variable in the model."""
         return self.__index
 
     def Proto(self):
-        """Returns the variable protobuf."""
+        """Return the variable protobuf."""
         return self.__var
 
     def __str__(self):
@@ -356,7 +363,7 @@ class IntVar(LinearExpression):
         return self.__var.name
 
     def Not(self):
-        """Returns the negation of a Boolean variable.
+        """Return the negation of a Boolean variable.
 
     This method implements the logical negation of a Boolean variable.
     It is only valid of the variable has a Boolean domain (0 or 1).
@@ -449,7 +456,7 @@ class Constraint(object):
         self.__constraint = constraints.add()
 
     def OnlyEnforceIf(self, boolvar):
-        """Adds an enforcement literal to the constraint.
+        """Add an enforcement literal to the constraint.
 
     Args:
         boolvar: A boolean literal or a list of boolean literals.
@@ -482,11 +489,11 @@ class Constraint(object):
         return self
 
     def Index(self):
-        """Returns the index of the constraint in the model."""
+        """Return the index of the constraint in the model."""
         return self.__index
 
     def Proto(self):
-        """Returns the constraint protobuf."""
+        """Return the constraint protobuf."""
         return self.__constraint
 
 
@@ -522,11 +529,11 @@ class IntervalVar(object):
             self.__ct.name = name
 
     def Index(self):
-        """Returns the index of the interval constraint in the model."""
+        """Return the index of the interval constraint in the model."""
         return self.__index
 
     def Proto(self):
-        """Returns the interval protobuf."""
+        """Return the interval protobuf."""
         return self.__ct.interval
 
     def __str__(self):
@@ -551,11 +558,11 @@ class IntervalVar(object):
 
 
 class CpModel(object):
-    """Wrapper class around the cp_model proto.
+    """Methods for building a CP model.
 
-  This class provides two types of methods:
-    - NewXXX to create integer, boolean, or interval variables.
-    - AddXXX to create new constraints and add them to the model.
+  Methods beginning with:
+  * ```New``` create integer, boolean, or interval variables.
+  * ```Add``` create new constraints and add them to the model.
   """
 
     def __init__(self):
@@ -566,10 +573,31 @@ class CpModel(object):
     # Domains
 
     def DomainFromValues(self, values):
-        """Builds a suitable domain from a list of values."""
+        """Build a suitable domain from a list of values.
+
+    Args:
+        values: a list of integer values.
+
+    Returns:
+        A flattened list of intervals representing the set of values.
+
+    This method build a flattened domain suitable to be used in
+    AddSumConstraintWithBounds and AddLinearConstraintWithBounds.
+    """
         return pywrapsat.SatHelper.DomainFromValues(values)
 
     def DomainFromIntervals(self, intervals):
+        """Build a suitable domain from a list of tuples (start, end).
+
+    Args:
+        intervals: a list of intervals (start, end) inclusive.
+
+    Returns:
+        A flattened list of intervals representing the union of the intervals.
+
+    This method build a flattened domain suitable to be used in
+    AddSumConstraintWithBounds and AddLinearConstraintWithBounds.
+    """
         starts = []
         ends = []
         for interval in intervals:
@@ -580,23 +608,27 @@ class CpModel(object):
     # Integer variable.
 
     def NewIntVar(self, lb, ub, name):
-        """Creates an integer variable with domain [lb, ub]."""
+        """Create an integer variable with domain [lb, ub]."""
         return IntVar(self.__model, [lb, ub], name)
 
-    def NewEnumeratedIntVar(self, bounds, name):
-        """Creates an integer variable with an enumerated domain.
+    def NewIntVarFromValues(self, values, name):
+        """Create an integer variable with domain {values}."""
+        return IntVar(self.__model, self.DomainFromValues(values), name)
+
+    def NewIntVarFromIntervals(self, intervals, name):
+        """Create an integer variable from a list of intervals.
 
     Args:
-        bounds: A flattened list of disjoint intervals.
+        intervals: A list of intervals (start, end) inclusive.
         name: The name of the variable.
 
     Returns:
-        a variable whose domain is union[bounds[2*i]..bounds[2*i + 1]].
+        a variable whose domain is the union of the intervals.
 
     To create a variable with domain [1, 2, 3, 5, 7, 8], pass in the
-    array [1, 3, 5, 5, 7, 8].
+    array [(1, 3), (5, 5), (7, 8)].
     """
-        return IntVar(self.__model, bounds, name)
+        return IntVar(self.__model, self.DomainFromIntervals(intervals), name)
 
     def NewBoolVar(self, name):
         """Creates a 0-1 variable with the given name."""
@@ -605,7 +637,7 @@ class CpModel(object):
     # Integer constraints.
 
     def AddLinearConstraint(self, terms, lb, ub):
-        """Adds the constraints lb <= sum(terms) <= ub, where term = (var, coef)."""
+        """Add the constraints lb <= sum(terms) <= ub, where term = (var, coef)."""
         ct = Constraint(self.__model.constraints)
         model_ct = self.__model.constraints[ct.Index()]
         for t in terms:
@@ -618,7 +650,7 @@ class CpModel(object):
         return ct
 
     def AddSumConstraint(self, variables, lb, ub):
-        """Adds the constraints lb <= sum(variables) <= ub."""
+        """Add the constraints lb <= sum(variables) <= ub."""
         ct = Constraint(self.__model.constraints)
         model_ct = self.__model.constraints[ct.Index()]
         for v in variables:
@@ -628,7 +660,7 @@ class CpModel(object):
         return ct
 
     def AddLinearConstraintWithBounds(self, terms, bounds):
-        """Adds the constraints sum(terms) in bounds, where term = (var, coef)."""
+        """Add the constraints sum(terms) in bounds, where term = (var, coef)."""
         ct = Constraint(self.__model.constraints)
         model_ct = self.__model.constraints[ct.Index()]
         for t in terms:
@@ -640,8 +672,20 @@ class CpModel(object):
         model_ct.linear.domain.extend(bounds)
         return ct
 
+    def AddSumConstraintWithBounds(self, variables, bounds):
+        """Add the constraints sum(variables) in bounds."""
+        ct = Constraint(self.__model.constraints)
+        model_ct = self.__model.constraints[ct.Index()]
+        for var in variables:
+            if not isinstance(var, IntVar):
+                raise TypeError('Wrong argument' + str(var))
+            model_ct.linear.vars.append(var.Index())
+            model_ct.linear.coeffs.append(1)
+        model_ct.linear.domain.extend(bounds)
+        return ct
+
     def Add(self, ct):
-        """Adds a LinearInequality to the model."""
+        """Add a LinearInequality to the model."""
         if isinstance(ct, LinearInequality):
             coeffs_map, constant = ct.Expression().GetVarValueMap()
             bounds = [cp_model_helper.CapSub(x, constant) for x in ct.Bounds()]
@@ -655,7 +699,7 @@ class CpModel(object):
             raise TypeError('Not supported: CpModel.Add(' + str(ct) + ')')
 
     def AddAllDifferent(self, variables):
-        """Adds AllDifferent(variables).
+        """Add AllDifferent(variables).
 
     This constraint forces all variables to have different values.
 
@@ -672,7 +716,7 @@ class CpModel(object):
         return ct
 
     def AddElement(self, index, variables, target):
-        """Adds the element constraint: variables[index] == target."""
+        """Add the element constraint: variables[index] == target."""
 
         if not variables:
             raise ValueError('AddElement expects a non empty variables array')
@@ -686,9 +730,9 @@ class CpModel(object):
         return ct
 
     def AddCircuit(self, arcs):
-        """Adds Circuit(arcs).
+        """Add Circuit(arcs).
 
-    Adds a circuit constraint from a sparse list of arcs that encode the graph.
+    Add a circuit constraint from a sparse list of arcs that encode the graph.
 
     A circuit is a unique Hamiltonian path in a subgraph of the total
     graph. In case a node 'i' is not in the path, then there must be a
@@ -721,7 +765,7 @@ class CpModel(object):
         return ct
 
     def AddAllowedAssignments(self, variables, tuples_list):
-        """Adds AllowedAssignments(variables, tuples_list).
+        """Add AllowedAssignments(variables, tuples_list).
 
     An AllowedAssignments constraint is a constraint on an array of variables
     that forces, when all variables are fixed to a single value, that the
@@ -761,7 +805,7 @@ class CpModel(object):
         return ct
 
     def AddForbiddenAssignments(self, variables, tuples_list):
-        """Adds AddForbiddenAssignments(variables, [tuples_list]).
+        """Add AddForbiddenAssignments(variables, [tuples_list]).
 
     A ForbiddenAssignments constraint is a constraint on an array of variables
     where the list of impossible combinations is provided in the tuples list.
@@ -793,7 +837,7 @@ class CpModel(object):
 
     def AddAutomaton(self, transition_variables, starting_state, final_states,
                      transition_triples):
-        """Adds an automaton constraint.
+        """Add an automaton constraint.
 
     An automaton constraint takes a list of variables (of size n), an initial
     state, a set of final states, and a set of transitions. A transition is a
@@ -865,7 +909,7 @@ class CpModel(object):
         return ct
 
     def AddInverse(self, variables, inverse_variables):
-        """Adds Inverse(variables, inverse_variables).
+        """Add Inverse(variables, inverse_variables).
 
     An inverse constraint enforces that if 'variables[i]' is assigned a value
     'j', then inverse_variables[j] is assigned a value 'i'. And vice versa.
@@ -898,7 +942,7 @@ class CpModel(object):
         return ct
 
     def AddReservoirConstraint(self, times, demands, min_level, max_level):
-        """Adds Reservoir(times, demands, min_level, max_level).
+        """Add Reservoir(times, demands, min_level, max_level).
 
     Maintains a reservoir level within bounds. The water level starts at 0, and
     at any time >= 0, it must be between min_level and max_level. Furthermore,
@@ -942,7 +986,7 @@ class CpModel(object):
 
     def AddReservoirConstraintWithActive(self, times, demands, actives,
                                          min_level, max_level):
-        """Adds Reservoir(times, demands, actives, min_level, max_level).
+        """Add Reservoir(times, demands, actives, min_level, max_level).
 
     Maintain a reservoir level within bounds. The water level starts at 0, and
     at
@@ -992,7 +1036,7 @@ class CpModel(object):
         return ct
 
     def AddMapDomain(self, var, bool_var_array, offset=0):
-        """Adds var == i + offset <=> bool_var_array[i] == true for all i."""
+        """Add var == i + offset <=> bool_var_array[i] == true for all i."""
 
         for i, bool_var in enumerate(bool_var_array):
             b_index = bool_var.Index()
@@ -1013,7 +1057,7 @@ class CpModel(object):
                 model_ct.linear.domain.extend([offset + i + 1, INT_MAX])
 
     def AddImplication(self, a, b):
-        """Adds a => b."""
+        """Add a => b."""
         ct = Constraint(self.__model.constraints)
         model_ct = self.__model.constraints[ct.Index()]
         model_ct.bool_or.literals.append(self.GetOrMakeBooleanIndex(b))
@@ -1021,7 +1065,7 @@ class CpModel(object):
         return ct
 
     def AddBoolOr(self, literals):
-        """Adds Or(literals) == true."""
+        """Add Or(literals) == true."""
         ct = Constraint(self.__model.constraints)
         model_ct = self.__model.constraints[ct.Index()]
         model_ct.bool_or.literals.extend(
@@ -1029,7 +1073,7 @@ class CpModel(object):
         return ct
 
     def AddBoolAnd(self, literals):
-        """Adds And(literals) == true."""
+        """Add And(literals) == true."""
         ct = Constraint(self.__model.constraints)
         model_ct = self.__model.constraints[ct.Index()]
         model_ct.bool_and.literals.extend(
@@ -1037,7 +1081,7 @@ class CpModel(object):
         return ct
 
     def AddBoolXOr(self, literals):
-        """Adds XOr(literals) == true."""
+        """Add XOr(literals) == true."""
         ct = Constraint(self.__model.constraints)
         model_ct = self.__model.constraints[ct.Index()]
         model_ct.bool_xor.literals.extend(
@@ -1045,7 +1089,7 @@ class CpModel(object):
         return ct
 
     def AddMinEquality(self, target, variables):
-        """Adds target == Min(variables)."""
+        """Add target == Min(variables)."""
         ct = Constraint(self.__model.constraints)
         model_ct = self.__model.constraints[ct.Index()]
         model_ct.int_min.vars.extend(
@@ -1054,7 +1098,7 @@ class CpModel(object):
         return ct
 
     def AddMaxEquality(self, target, args):
-        """Adds target == Max(variables)."""
+        """Add target == Max(variables)."""
         ct = Constraint(self.__model.constraints)
         model_ct = self.__model.constraints[ct.Index()]
         model_ct.int_max.vars.extend([self.GetOrMakeIndex(x) for x in args])
@@ -1062,7 +1106,7 @@ class CpModel(object):
         return ct
 
     def AddDivisionEquality(self, target, num, denom):
-        """Adds target == num // denom."""
+        """Add target == num // denom."""
         ct = Constraint(self.__model.constraints)
         model_ct = self.__model.constraints[ct.Index()]
         model_ct.int_div.vars.extend(
@@ -1072,7 +1116,7 @@ class CpModel(object):
         return ct
 
     def AddAbsEquality(self, target, var):
-        """Adds target == Abs(var)."""
+        """Add target == Abs(var)."""
         ct = Constraint(self.__model.constraints)
         model_ct = self.__model.constraints[ct.Index()]
         index = self.GetOrMakeIndex(var)
@@ -1081,7 +1125,7 @@ class CpModel(object):
         return ct
 
     def AddModuloEquality(self, target, var, mod):
-        """Adds target = var % mod."""
+        """Add target = var % mod."""
         ct = Constraint(self.__model.constraints)
         model_ct = self.__model.constraints[ct.Index()]
         model_ct.int_mod.vars.extend(
@@ -1091,7 +1135,7 @@ class CpModel(object):
         return ct
 
     def AddProdEquality(self, target, args):
-        """Adds target == PROD(args)."""
+        """Add target == PROD(args)."""
         ct = Constraint(self.__model.constraints)
         model_ct = self.__model.constraints[ct.Index()]
         model_ct.int_prod.vars.extend([self.GetOrMakeIndex(x) for x in args])
@@ -1158,7 +1202,7 @@ class CpModel(object):
                            is_present_index, name)
 
     def AddNoOverlap(self, interval_vars):
-        """Adds NoOverlap(interval_vars).
+        """Add NoOverlap(interval_vars).
 
     A NoOverlap constraint ensures that all present intervals do not overlap
     in time.
@@ -1176,7 +1220,7 @@ class CpModel(object):
         return ct
 
     def AddNoOverlap2D(self, x_intervals, y_intervals):
-        """Adds NoOverlap2D(x_intervals, y_intervals).
+        """Add NoOverlap2D(x_intervals, y_intervals).
 
     A NoOverlap2D constraint ensures that all present rectangles do not overlap
     on a plan. Each rectangle is aligned with the X and Y axis, and is defined
@@ -1198,7 +1242,7 @@ class CpModel(object):
         return ct
 
     def AddCumulative(self, intervals, demands, capacity):
-        """Adds Cumulative(intervals, demands, capacity).
+        """Add Cumulative(intervals, demands, capacity).
 
     This constraint enforces that:
       for all t:
@@ -1231,14 +1275,14 @@ class CpModel(object):
         return str(self.__model)
 
     def Proto(self):
-        """Returns the underling CpModelProto."""
+        """Return the underling CpModelProto."""
         return self.__model
 
     def Negated(self, index):
         return -index - 1
 
     def GetOrMakeIndex(self, arg):
-        """Returns the index of a variables, its negation, or a number."""
+        """Return the index of a variables, its negation, or a number."""
         if isinstance(arg, IntVar):
             return arg.Index()
         elif (isinstance(arg, _ProductCst) and
@@ -1252,7 +1296,7 @@ class CpModel(object):
                 'NotSupported: model.GetOrMakeIndex(' + str(arg) + ')')
 
     def GetOrMakeBooleanIndex(self, arg):
-        """Returns an index from a boolean expression."""
+        """Return an index from a boolean expression."""
         if isinstance(arg, IntVar):
             self.AssertIsBooleanVariable(arg)
             return arg.Index()
@@ -1332,7 +1376,7 @@ class CpModel(object):
         return self.__model.HasField('objective')
 
     def AddDecisionStrategy(self, variables, var_strategy, domain_strategy):
-        """Adds a search strategy to the model.
+        """Add a search strategy to the model.
 
     Args:
       variables: a list of variables this strategy will assign.
@@ -1350,11 +1394,11 @@ class CpModel(object):
         strategy.domain_reduction_strategy = domain_strategy
 
     def ModelStats(self):
-        """Returns some statistics on the model as a string."""
+        """Return some statistics on the model as a string."""
         return pywrapsat.SatHelper.ModelStats(self.__model)
 
     def Validate(self):
-        """Returns a string explaining the issue is the model is not valid."""
+        """Return a string explaining the issue is the model is not valid."""
         return pywrapsat.SatHelper.ValidateModel(self.__model)
 
     def AssertIsBooleanVariable(self, x):
@@ -1422,13 +1466,13 @@ class CpSolver(object):
         self.parameters = sat_parameters_pb2.SatParameters()
 
     def Solve(self, model):
-        """Solves the given model and returns the solve status."""
+        """Solve the given model and Return the solve status."""
         self.__solution = pywrapsat.SatHelper.SolveWithParameters(
             model.Proto(), self.parameters)
         return self.__solution.status
 
     def SolveWithSolutionCallback(self, model, callback):
-        """Solves a problem and pass each solution found to the callback."""
+        """Solve a problem and pass each solution found to the callback."""
         self.__solution = (
             pywrapsat.SatHelper.SolveWithParametersAndSolutionCallback(
                 model.Proto(), self.parameters, callback))
@@ -1463,55 +1507,55 @@ class CpSolver(object):
         return self.__solution.status
 
     def Value(self, expression):
-        """Returns the value of an linear expression after solve."""
+        """Return the value of an linear expression after solve."""
         if not self.__solution:
             raise RuntimeError('Solve() has not be called.')
         return EvaluateLinearExpression(expression, self.__solution)
 
     def BooleanValue(self, literal):
-        """Returns the boolean value of a literal after solve."""
+        """Return the boolean value of a literal after solve."""
         if not self.__solution:
             raise RuntimeError('Solve() has not be called.')
         return EvaluateBooleanExpression(literal, self.__solution)
 
     def ObjectiveValue(self):
-        """Returns the value of objective after solve."""
+        """Return the value of objective after solve."""
         return self.__solution.objective_value
 
     def BestObjectiveBound(self):
-        """Returns the best lower (upper) bound found when min(max)imizing."""
+        """Return the best lower (upper) bound found when min(max)imizing."""
         return self.__solution.best_objective_bound
 
     def StatusName(self, status):
-        """Returns the name of the status returned by Solve()."""
+        """Return the name of the status returned by Solve()."""
         return cp_model_pb2.CpSolverStatus.Name(status)
 
     def NumBooleans(self):
-        """Returns the number of boolean variables managed by the SAT solver."""
+        """Return the number of boolean variables managed by the SAT solver."""
         return self.__solution.num_booleans
 
     def NumConflicts(self):
-        """Returns the number of conflicts since the creation of the solver."""
+        """Return the number of conflicts since the creation of the solver."""
         return self.__solution.num_conflicts
 
     def NumBranches(self):
-        """Returns the number of search branches explored by the solver."""
+        """Return the number of search branches explored by the solver."""
         return self.__solution.num_branches
 
     def WallTime(self):
-        """Returns the wall time in seconds since the creation of the solver."""
+        """Return the wall time in seconds since the creation of the solver."""
         return self.__solution.wall_time
 
     def UserTime(self):
-        """Returns the user time in seconds since the creation of the solver."""
+        """Return the user time in seconds since the creation of the solver."""
         return self.__solution.user_time
 
     def ResponseStats(self):
-        """Returns some statistics on the solution found as a string."""
+        """Return some statistics on the solution found as a string."""
         return pywrapsat.SatHelper.SolverResponseStats(self.__solution)
 
     def ResponseProto(self):
-        """Returns the response object."""
+        """Return the response object."""
         return self.__solution
 
 
@@ -1531,7 +1575,7 @@ class CpSolverSolutionCallback(pywrapsat.SolutionCallback):
         self.on_solution_callback()
 
     def BooleanValue(self, lit):
-        """Returns the boolean value of a boolean literal.
+        """Return the boolean value of a boolean literal.
 
     Args:
         lit: A boolean variable or its negation.
diff --git a/ortools/sat/python/sat.i b/ortools/sat/python/sat.i
index be7ac6cb08..48ae6b3161 100644
--- a/ortools/sat/python/sat.i
+++ b/ortools/sat/python/sat.i
@@ -20,7 +20,7 @@
 
 // std::function utilities.
 %include "ortools/util/python/functions.i"
-%include "ortools/util/python/vector.i"
+%import "ortools/util/python/vector.i"
 
 %{
 #include "ortools/sat/cp_model.pb.h"
@@ -93,3 +93,4 @@ PY_PROTO_TYPEMAP(ortools.sat.sat_parameters_pb2,
 %include "ortools/sat/swig_helper.h"
 
 %unignoreall
+
diff --git a/ortools/sat/samples/StepFunctionSampleSat.cs b/ortools/sat/samples/StepFunctionSampleSat.cs
index a2c5f09064..5dfeab6562 100644
--- a/ortools/sat/samples/StepFunctionSampleSat.cs
+++ b/ortools/sat/samples/StepFunctionSampleSat.cs
@@ -58,18 +58,17 @@ public class StepFunctionSampleSat
 
     // expr == 0 on [5, 6] U [8, 10]
     ILiteral b0 = model.NewBoolVar("b0");
-    model.AddLinearConstraintWithBounds(
+    model.AddLinearSumWithBounds(
         new IntVar[] {x},
-        new long[] {1},
-        new long[] {5, 6, 8, 10}).OnlyEnforceIf(b0);
+        model.DomainFromIntervals(new long[] {5, 6, 8, 10})).OnlyEnforceIf(b0);
     model.Add(expr == 0).OnlyEnforceIf(b0);
 
     // expr == 2 on [0, 1] U [3, 4] U [11, 20]
     ILiteral b2 = model.NewBoolVar("b2");
-    model.AddLinearConstraintWithBounds(
+    model.AddLinearSumWithBounds(
         new IntVar[] {x},
-        new long[] {1},
-        new long[] {0, 1, 3, 4, 11, 20}).OnlyEnforceIf(b2);
+        model.DomainFromIntervals(
+            new long[] {0, 1, 3, 4, 11, 20})).OnlyEnforceIf(b2);
     model.Add(expr == 2).OnlyEnforceIf(b2);
 
     // expr == 3 when x == 7
diff --git a/ortools/sat/samples/binpacking_problem_sat.py b/ortools/sat/samples/binpacking_problem_sat.py
index 9e98984de1..9dd0101296 100644
--- a/ortools/sat/samples/binpacking_problem_sat.py
+++ b/ortools/sat/samples/binpacking_problem_sat.py
@@ -12,8 +12,6 @@
 # limitations under the License.
 """Solves a binpacking problem using the CP-SAT solver."""
 
-from __future__ import absolute_import
-from __future__ import division
 from __future__ import print_function
 
 from ortools.sat.python import cp_model
diff --git a/ortools/sat/samples/minimal_jobshop_sat.py b/ortools/sat/samples/minimal_jobshop_sat.py
index 337da83dee..ecf3287d17 100644
--- a/ortools/sat/samples/minimal_jobshop_sat.py
+++ b/ortools/sat/samples/minimal_jobshop_sat.py
@@ -13,7 +13,6 @@
 """Minimal jobshop example."""
 
 # [START program]
-from __future__ import absolute_import
 from __future__ import print_function
 
 import collections
diff --git a/ortools/sat/samples/nurses_sat.py b/ortools/sat/samples/nurses_sat.py
index a916c72943..6949522683 100644
--- a/ortools/sat/samples/nurses_sat.py
+++ b/ortools/sat/samples/nurses_sat.py
@@ -13,9 +13,8 @@
 """Example of a simple nurse scheduling problem."""
 
 # [START program]
-from __future__ import division
-from __future__ import print_function
 # [START import]
+from __future__ import print_function
 from ortools.sat.python import cp_model
 
 # [END import]
diff --git a/ortools/sat/samples/ranking_sample_sat.py b/ortools/sat/samples/ranking_sample_sat.py
index 4ff56d5559..0a922e6b3d 100644
--- a/ortools/sat/samples/ranking_sample_sat.py
+++ b/ortools/sat/samples/ranking_sample_sat.py
@@ -98,7 +98,7 @@ def RankingSampleSat():
         start = model.NewIntVar(0, horizon, 'start_%i' % t)
         duration = t + 1
         end = model.NewIntVar(0, horizon, 'end_%i' % t)
-        if t < num_tasks / 2:
+        if t < num_tasks // 2:
             interval = model.NewIntervalVar(start, duration, end,
                                             'interval_%i' % t)
             presence = True
diff --git a/ortools/sat/samples/schedule_requests_sat.py b/ortools/sat/samples/schedule_requests_sat.py
index e0e7a7e44f..41314b5aae 100644
--- a/ortools/sat/samples/schedule_requests_sat.py
+++ b/ortools/sat/samples/schedule_requests_sat.py
@@ -13,9 +13,8 @@
 """Nurse scheduling problem with shift requests."""
 
 # [START program]
-from __future__ import division
-from __future__ import print_function
 # [START import]
+from __future__ import print_function
 from ortools.sat.python import cp_model
 
 # [END import]
diff --git a/ortools/sat/samples/scheduling_with_calendar_sample_sat.py b/ortools/sat/samples/scheduling_with_calendar_sample_sat.py
index 66088cd057..37a73a96c8 100644
--- a/ortools/sat/samples/scheduling_with_calendar_sample_sat.py
+++ b/ortools/sat/samples/scheduling_with_calendar_sample_sat.py
@@ -44,23 +44,22 @@ def SchedulingWithCalendarSampleSat():
     # The data is the following:
     #   Work starts at 8h, ends at 18h, with a lunch break between 13h and 14h.
     #   We need to schedule a task that needs 3 hours of processing time.
-    #   Total duration can be 3 or 4 (if it span across the lunch break.
+    #   Total duration can be 3 or 4 (if it spans the lunch break).
     #
-    # Because the duration is at least 3, work cannot start after 15h.
+    # Because the duration is at least 3 hours, work cannot start after 15h.
     # Because of the break, work cannot start at 13h.
 
-    start = model.NewEnumeratedIntVar([8, 12, 14, 15], 'start')
+    start = model.NewIntVarFromIntervals([(8, 12), (14, 15)], 'start')
     duration = model.NewIntVar(3, 4, 'duration')
     end = model.NewIntVar(8, 18, 'end')
     unused_interval = model.NewIntervalVar(start, duration, end, 'interval')
 
     # We have 2 states (spanning across lunch or not)
     across = model.NewBoolVar('across')
-    model.AddLinearConstraintWithBounds([(start, 1)],
-                                        [8, 10, 14, 15]).OnlyEnforceIf(
-                                            across.Not())
-    model.AddLinearConstraintWithBounds([(start, 1)],
-                                        [11, 12]).OnlyEnforceIf(across)
+    non_spanning_hours = model.DomainFromValues([8, 9, 10, 14, 15])
+    model.AddSumConstraintWithBounds([start], non_spanning_hours).OnlyEnforceIf(
+        across.Not())
+    model.AddSumConstraint([start], 11, 12).OnlyEnforceIf(across)
     model.Add(duration == 3).OnlyEnforceIf(across.Not())
     model.Add(duration == 4).OnlyEnforceIf(across)
 
@@ -74,7 +73,7 @@ def SchedulingWithCalendarSampleSat():
     # Force the solver to follow the decision strategy exactly.
     solver.parameters.search_branching = cp_model.FIXED_SEARCH
 
-    # Search and print out all solutions.
+    # Search and print all solutions.
     solution_printer = VarArraySolutionPrinter([start, duration, across])
     solver.SearchForAllSolutions(model, solution_printer)
 
diff --git a/ortools/sat/samples/step_function_sample_sat.py b/ortools/sat/samples/step_function_sample_sat.py
index e5dfdfa95f..c3cd93fcba 100644
--- a/ortools/sat/samples/step_function_sample_sat.py
+++ b/ortools/sat/samples/step_function_sample_sat.py
@@ -59,14 +59,17 @@ def step_function_sample_sat():
 
     # expr == 0 on [5, 6] U [8, 10]
     b0 = model.NewBoolVar('b0')
-    model.AddLinearConstraintWithBounds([(x, 1)],
-                                        [5, 6, 8, 10]).OnlyEnforceIf(b0)
+    model.AddSumConstraintWithBounds([x],
+                                     model.DomainFromIntervals(
+                                         [(5, 6), (8, 10)])).OnlyEnforceIf(b0)
     model.Add(expr == 0).OnlyEnforceIf(b0)
 
     # expr == 2 on [0, 1] U [3, 4] U [11, 20]
     b2 = model.NewBoolVar('b2')
-    model.AddLinearConstraintWithBounds([(x, 1)],
-                                        [0, 1, 3, 4, 11, 20]).OnlyEnforceIf(b2)
+    model.AddSumConstraintWithBounds([x],
+                                     model.DomainFromIntervals(
+                                         [(0, 1), (3, 4),
+                                          (11, 20)])).OnlyEnforceIf(b2)
     model.Add(expr == 2).OnlyEnforceIf(b2)
 
     # expr == 3 when x == 7
diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto
index 09f50fd5c4..8cb6bfeeca 100644
--- a/ortools/sat/sat_parameters.proto
+++ b/ortools/sat/sat_parameters.proto
@@ -369,7 +369,8 @@ message SatParameters {
 
   // The maximum "deterministic" time limit to spend in probing. A value of
   // zero will disable the probing.
-  optional double presolve_probing_deterministic_time_limit = 57 [default = 30.0];
+  optional double presolve_probing_deterministic_time_limit = 57
+      [default = 30.0];
 
   // Whether we use an heuristic to detect some basic case of blocked clause
   // in the SAT presolve.
diff --git a/ortools/sat/swig_helper.h b/ortools/sat/swig_helper.h
index efa572aedd..95318e721b 100644
--- a/ortools/sat/swig_helper.h
+++ b/ortools/sat/swig_helper.h
@@ -210,8 +210,7 @@ class SatHelper {
   // Builds a flattened list of intervals from two parallel vectors of starts
   // and ends.
   static std::vector<int64> DomainFromStartsAndEnds(
-      const std::vector<int64>& starts,
-      const std::vector<int64>& ends) {
+      const std::vector<int64>& starts, const std::vector<int64>& ends) {
     std::vector<ClosedInterval> input;
     for (int i = 0; i < starts.size(); ++i) {
       ClosedInterval tmp;
diff --git a/ortools/sat/util.cc b/ortools/sat/util.cc
index 2d313bb424..32852ed3ab 100644
--- a/ortools/sat/util.cc
+++ b/ortools/sat/util.cc
@@ -55,5 +55,15 @@ int MoveOneUnprocessedLiteralLast(const std::set<LiteralIndex>& processed,
   return target_prefix_size;
 }
 
+void IncrementalAverage::Reset(double reset_value) {
+  num_records_ = 0;
+  average_ = reset_value;
+}
+
+void IncrementalAverage::AddData(double new_record) {
+  num_records_++;
+  average_ += (new_record - average_) / num_records_;
+}
+
 }  // namespace sat
 }  // namespace operations_research
diff --git a/ortools/sat/util.h b/ortools/sat/util.h
index 59b5a1ef28..d516f693f9 100644
--- a/ortools/sat/util.h
+++ b/ortools/sat/util.h
@@ -97,6 +97,27 @@ inline void RandomizeDecisionHeuristic(URBG* random,
                                                                       : 0.0);
 }
 
+// Manages incremental averages.
+class IncrementalAverage {
+ public:
+  // Initializes the average with 'initial_average' and number of records to 0.
+  explicit IncrementalAverage(double initial_average)
+      : average_(initial_average) {}
+  IncrementalAverage() {}
+
+  // Sets the number of records to 0 and average to 'reset_value'.
+  void Reset(double reset_value);
+
+  double CurrentAverage() const { return average_; }
+  int64 NumRecords() const { return num_records_; }
+
+  void AddData(double new_record);
+
+ private:
+  double average_ = 0.0;
+  int64 num_records_ = 0;
+};
+
 }  // namespace sat
 }  // namespace operations_research