Add the following features to LLS

  sum of squares expression now support adding a constant scalar and
  construction with a scalar
  constraints lists can now be initialized as empty array []

docs/_build/html/_sources/julia_examples.txt

@@ -238,7 +238,7 @@ The following code builds and solves our control example:
 
   # Create a problem instance
   mu = 1
-  constraints = EqConstraint[]
+  constraints = []
 
   # Add constraints on our variables
   for i in 1 : T - 1
@@ -462,8 +462,8 @@ indicative of sports articles and which were most indicative of nonsports.
     print(words[indices[length(words) - i]])
   end
 
-The 5 strings with largest positive weights were "play", "peopl", "olymp", "nativ", and "fan".
-The 5 strings with largest negative weights were "get", "ml", "issu", "professor", and "student".
+Each run will yield different words, but it'll be clear which words
+come from sports articles.
 
 
 Time Series Analysis
@@ -518,7 +518,7 @@ The following code uses LLS to find and plot the model:
   # draw(PNG("melbourne.png", 16cm, 12cm), p)
 
   yearly = Variable(n)
-  eq_constraints = EqConstraint[]
+  eq_constraints = []
   for i in 365 + 1 : n
     eq_constraints += yearly[i] == yearly[i - 365]
   end

docs/_build/html/_sources/julia_tutorial.txt

@@ -175,11 +175,13 @@ to a list using the ``+`` operator.
   constraint_list = [A * x == randn(4, 1), 3 == x[1:2]]
   constraint_list += x[3] == 1.6
 
-An empty list of constraints can be created with
+An empty list of constraints can be created with ``[]``. You can add to an empty
+list with the same syntax.
 
 .. code-block:: none
 
-  empty_list = EqConstraint[]
+  new_list = []
+  new_list += x[2] == 1.2
 
 
 The ``solve!`` Function
@@ -216,6 +218,7 @@ conjunction with the following rules:
 
 #. Two sum of squares expressions can be added
 #. A sum of squares expression can be multiplied or divided by a postive, scalar constant.
+#. A nonnegative scalar constant may be added to a sum of squares expression.
 
 Note that sum of squares expression cannot be subtracted from each other,
 or multiplied or divided by a negative number. LLS will issue an error message if
@@ -240,6 +243,16 @@ assuming ``x`` has been populated with a value:
 
   println(evaluate(reg_least_squares))
 
+Often you'll find it useful to first initialize a sum of squares expression
+to ``0`` and then add on more sum of squares expressions in a for loop.
+
+.. code-block::none
+
+  error_term = 0
+  for i in 1:3
+    error_term += rand() * sum_squares(A[i, :] * x + b[i])
+  end
+
 The variance of the entries of an affine expression ``X`` can be expressed as
 ``sum_squares(mean(X) - X) / (m * n)``, where ``m`` and ``n`` are the number of rows
 and number of columns of ``X``, respectively. For convenience, the function ``var``

docs/_build/html/julia_examples.html

@@ -244,7 +244,7 @@ <h2>Control<a class="headerlink" href="#control" title="Permalink to this headli
 
 # Create a problem instance
 mu = 1
-constraints = EqConstraint[]
+constraints = []
 
 # Add constraints on our variables
 for i in 1 : T - 1
@@ -439,8 +439,8 @@ <h3>Binary Classification<a class="headerlink" href="#binary-classification" tit
 end
 </pre></div>
 </div>
-<p>The 5 strings with largest positive weights were &#8220;play&#8221;, &#8220;peopl&#8221;, &#8220;olymp&#8221;, &#8220;nativ&#8221;, and &#8220;fan&#8221;.
-The 5 strings with largest negative weights were &#8220;get&#8221;, &#8220;ml&#8221;, &#8220;issu&#8221;, &#8220;professor&#8221;, and &#8220;student&#8221;.</p>
+<p>Each run will yield different words, but it&#8217;ll be clear which words
+come from sports articles.</p>
 </div>
 </div>
 <div class="section" id="time-series-analysis">
@@ -483,7 +483,7 @@ <h2>Time Series Analysis<a class="headerlink" href="#time-series-analysis" title
 # draw(PNG(&quot;melbourne.png&quot;, 16cm, 12cm), p)
 
 yearly = Variable(n)
-eq_constraints = EqConstraint[]
+eq_constraints = []
 for i in 365 + 1 : n
   eq_constraints += yearly[i] == yearly[i - 365]
 end

docs/_build/html/julia_tutorial.html

@@ -196,8 +196,10 @@ <h2>Linear Equality Constraints<a class="headerlink" href="#linear-equality-cons
 constraint_list += x[3] == 1.6
 </pre></div>
 </div>
-<p>An empty list of constraints can be created with</p>
-<div class="highlight-none"><div class="highlight"><pre>empty_list = EqConstraint[]
+<p>An empty list of constraints can be created with <tt class="docutils literal"><span class="pre">[]</span></tt>. You can add to an empty
+list with the same syntax.</p>
+<div class="highlight-none"><div class="highlight"><pre>new_list = []
+new_list += x[2] == 1.2
 </pre></div>
 </div>
 </div>
@@ -232,6 +234,7 @@ <h2>Sum of Squares Expressions<a class="headerlink" href="#sum-of-squares-expres
 <ol class="arabic simple">
 <li>Two sum of squares expressions can be added</li>
 <li>A sum of squares expression can be multiplied or divided by a postive, scalar constant.</li>
+<li>A nonnegative scalar constant may be added to a sum of squares expression.</li>
 </ol>
 <p>Note that sum of squares expression cannot be subtracted from each other,
 or multiplied or divided by a negative number. LLS will issue an error message if
@@ -252,6 +255,8 @@ <h2>Sum of Squares Expressions<a class="headerlink" href="#sum-of-squares-expres
 <div class="highlight-none"><div class="highlight"><pre>println(evaluate(reg_least_squares))
 </pre></div>
 </div>
+<p>Often you&#8217;ll find it useful to first initialize a sum of squares expression
+to <tt class="docutils literal"><span class="pre">0</span></tt> and then add on more sum of squares expressions in a for loop.</p>
 <p>The variance of the entries of an affine expression <tt class="docutils literal"><span class="pre">X</span></tt> can be expressed as
 <tt class="docutils literal"><span class="pre">sum_squares(mean(X)</span> <span class="pre">-</span> <span class="pre">X)</span> <span class="pre">/</span> <span class="pre">(m</span> <span class="pre">*</span> <span class="pre">n)</span></tt>, where <tt class="docutils literal"><span class="pre">m</span></tt> and <tt class="docutils literal"><span class="pre">n</span></tt> are the number of rows
 and number of columns of <tt class="docutils literal"><span class="pre">X</span></tt>, respectively. For convenience, the function <tt class="docutils literal"><span class="pre">var</span></tt>

docs/julia_examples.rst

@@ -238,7 +238,7 @@ The following code builds and solves our control example:
 
   # Create a problem instance
   mu = 1
-  constraints = EqConstraint[]
+  constraints = []
 
   # Add constraints on our variables
   for i in 1 : T - 1
@@ -462,8 +462,8 @@ indicative of sports articles and which were most indicative of nonsports.
     print(words[indices[length(words) - i]])
   end
 
-The 5 strings with largest positive weights were "play", "peopl", "olymp", "nativ", and "fan".
-The 5 strings with largest negative weights were "get", "ml", "issu", "professor", and "student".
+Each run will yield different words, but it'll be clear which words
+come from sports articles.
 
 
 Time Series Analysis
@@ -518,7 +518,7 @@ The following code uses LLS to find and plot the model:
   # draw(PNG("melbourne.png", 16cm, 12cm), p)
 
   yearly = Variable(n)
-  eq_constraints = EqConstraint[]
+  eq_constraints = []
   for i in 365 + 1 : n
     eq_constraints += yearly[i] == yearly[i - 365]
   end

docs/julia_tutorial.rst

@@ -175,11 +175,13 @@ to a list using the ``+`` operator.
   constraint_list = [A * x == randn(4, 1), 3 == x[1:2]]
   constraint_list += x[3] == 1.6
 
-An empty list of constraints can be created with
+An empty list of constraints can be created with ``[]``. You can add to an empty
+list with the same syntax.
 
 .. code-block:: none
 
-  empty_list = EqConstraint[]
+  new_list = []
+  new_list += x[2] == 1.2
 
 
 The ``solve!`` Function
@@ -216,6 +218,7 @@ conjunction with the following rules:
 
 #. Two sum of squares expressions can be added
 #. A sum of squares expression can be multiplied or divided by a postive, scalar constant.
+#. A nonnegative scalar constant may be added to a sum of squares expression.
 
 Note that sum of squares expression cannot be subtracted from each other,
 or multiplied or divided by a negative number. LLS will issue an error message if
@@ -240,6 +243,16 @@ assuming ``x`` has been populated with a value:
 
   println(evaluate(reg_least_squares))
 
+Often you'll find it useful to first initialize a sum of squares expression
+to ``0`` and then add on more sum of squares expressions in a for loop.
+
+.. code-block::none
+
+  error_term = 0
+  for i in 1:3
+    error_term += rand() * sum_squares(A[i, :] * x + b[i])
+  end
+
 The variance of the entries of an affine expression ``X`` can be expressed as
 ``sum_squares(mean(X) - X) / (m * n)``, where ``m`` and ``n`` are the number of rows
 and number of columns of ``X``, respectively. For convenience, the function ``var``

examples/control/control.jl

@@ -26,7 +26,7 @@ force = Variable(2, T - 1)
 
 # Create a problem instance
 mu = 1
-constraints = EqConstraint[]
+constraints = []
 
 # Add constraints on our variables
 for i in 1 : T - 1

examples/time_series/time_series.jl

@@ -10,7 +10,7 @@ p = plot(
 # draw(PNG("melbourne.png", 16cm, 12cm), p)
 
 yearly = Variable(n)
-eq_constraints = EqConstraint[]
+eq_constraints = []
 for i in 365 + 1 : n
   eq_constraints += yearly[i] == yearly[i - 365]
 end

src/atoms/add_subtract.jl

@@ -81,10 +81,25 @@ function +(x::SumSquaresExpr, y::SumSquaresExpr)
   multipliers = copy(x.multipliers)
   append!(affines, y.affines)
   append!(multipliers, y.multipliers)
-  this = SumSquaresExpr(:+, affines, multipliers)
+  this = SumSquaresExpr(:+, affines, multipliers, x.scalar + y.scalar)
   return this
 end
 
+function +(x::Number, y::SumSquaresExpr)
+  try
+    x = convert(Float64, x)
+  catch
+    error("Only real scalars can be added to Sum Squares expressions")
+  end
+  if x < 0
+    error("Only nonnegative scalars can be added to Sum Squares expressions")
+  end
+  this = SumSquaresExpr(:+, y.affines, y.multipliers, y.scalar + x)
+  return this
+end
+
++(x::SumSquaresExpr, y::Number) = +(y, x)
+
 
 # Binary Subtraction
 -(x::AffineExpr, y::AffineExpr) = +(x, -y)

src/atoms/multiply_divide.jl

@@ -118,7 +118,7 @@ function *(x::Number, y::SumSquaresExpr)
   if x < 0
     error("Sum Squares expressions can only be multiplied by nonegative scalars")
   end
-  return SumSquaresExpr(:*, y.affines, x * y.multipliers)
+  return SumSquaresExpr(:*, y.affines, x * y.multipliers, x * y.scalar)
 end
 
 *(x::SumSquaresExpr, y::Number) = *(y, x)
@@ -132,5 +132,5 @@ function /(x::SumSquaresExpr, y::Number)
   if y <= 0
     error("Sum Squares expressions can only be divided by positive scalars")
   end
-  return SumSquaresExpr(:/, x.affines, x.multipliers / y)
+  return SumSquaresExpr(:/, x.affines, x.multipliers / y, x.scalar / y)
 end

src/solve.jl

@@ -2,7 +2,7 @@ export minimize!, solve!
 
 # Eventually, this can take in arguments to specify solving method.
 # For now, backslash is the only supported method.
-function solve!(problem_type::Symbol, objective::SumSquaresExpr, constraints::Array{EqConstraint})
+function solve!(problem_type::Symbol, objective::SumSquaresExpr, constraints::Array{EqConstraint, 1})
   p = Problem(problem_type, objective, constraints)
   # return backslash_solve!(p)
   backslash_solve!(p)
@@ -16,11 +16,11 @@ function solve!(problem_type::Symbol, objective::SumSquaresExpr, constraints::Ar
 end
 
 
-minimize!(objective::SumSquaresExpr, constraints::Array{EqConstraint}) = solve!(:minimize, objective, constraints)
+minimize!(objective::SumSquaresExpr, constraints::Array{EqConstraint, 1}) = solve!(:minimize, objective, constraints)
 minimize!(objective::SumSquaresExpr, constraint::EqConstraint) = minimize!(objective, [constraint])
 minimize!(objective::SumSquaresExpr, constraints::EqConstraint...) = minimize!(objective, [constraints...])
 minimize!(objective::SumSquaresExpr) = minimize!(objective, EqConstraint[])
 
-solve!(constraints::Array{EqConstraint}) = solve!(:solve, SumSquaresExpr(:empty, AffineExpr[], Float64[]), constraints)
+solve!(constraints::Array{EqConstraint, 1}) = solve!(:solve, sum_squares(), constraints)
 solve!(constraint::EqConstraint) = solve!([constraint])
 solve!(constraints::EqConstraint...) = solve!([constraints...])

src/solvers/backslash_solver.jl

@@ -126,10 +126,9 @@ function backslash_solve!(p::Problem)
     p.status = "KKT singular"
   else
     p.status = "solved"
+    p.optval = p.objective.scalar
     if size(A, 1) > 0
-      p.optval = sum((A*solution[1:num_vars] + b).^2)
-    else
-      p.optval = 0.0
+      p.optval += sum((A*solution[1:num_vars] + b).^2)
     end
     populate_vars!(unique_vars_map, vars_to_ranges_map, solution)
   end

src/types/constraints.jl

@@ -12,6 +12,9 @@ type EqConstraint
     if lhs.head == :Constant && rhs.head == :Constant
       error ("Equality constraints between two constants are not allowed")
     end
+    if lhs.size != rhs.size && lhs.size != (1, 1) && rhs.size != (1, 1)
+      error ("LHS and RHS of constraints muust have the same size, or one needs to be scalar")
+    end
     this = new(:(==), lhs, rhs, lhs - rhs)
     return this
   end
@@ -21,5 +24,7 @@ end
 ==(lhs::Numeric, rhs::AffineExpr) = ==(convert(Constant, lhs), rhs)
 ==(lhs::AffineExpr, rhs::Numeric) = ==(lhs, convert(Constant, rhs))
 
-+(constraints::Array{EqConstraint}, new_constraints::Array{EqConstraint}) = append!(constraints, new_constraints)
-+(constraints::Array{EqConstraint}, new_constraint::EqConstraint) = push!(constraints, new_constraint)
++(constraints::Array{EqConstraint, 1}, new_constraints::Array{EqConstraint, 1}) = append!(constraints, new_constraints)
++(constraints::Array{EqConstraint, 1}, new_constraint::EqConstraint) = push!(constraints, new_constraint)
++(constraints::Array{None, 1}, new_constraints::Array{EqConstraint, 1}) = new_constraints
++(constraints::Array{None, 1}, new_constraint::EqConstraint) = [new_constraint]

src/types/expressions.jl

@@ -106,34 +106,51 @@ type SumSquaresExpr <: AbstractExpr
   value::ValueOrNothing
   affines::Array{AffineExpr}
   multipliers::Array{Float64}
+  scalar::Float64
   uid::Uint64
   evaluate::Function
 
-  function SumSquaresExpr(head::Symbol, affines::Array{AffineExpr}, multipliers::Array{Float64})
-    this = new(head, nothing, affines, multipliers)
+  function SumSquaresExpr(head::Symbol, affines::Array{AffineExpr}, multipliers::Array{Float64}, scalar::Float64)
+    this = new(head, nothing, affines, multipliers, scalar)
     this.uid = object_id(this)
     this.evaluate = ()->begin
       value = 0
       for i in 1:length(this.affines)
         value += multipliers[i] * sum((this.affines[i].evaluate()).^2)
       end
-      return value
+      return value + this.scalar
     end
     return this
   end
 end
 
+function sum_squares()
+  this = SumSquaresExpr(:sum_squares, AffineExpr[], Float64[], 0.0)
+end
+
+function sum_squares(num::Number)
+  try
+    num = convert(Float64, num)
+  catch
+    error("Sum Squares expressions can only be constructed from real numbers")
+  end
+  if num < 0
+    error("Sum Squares expressions can only be constructed from nonnegative scalars")
+  end
+  this = SumSquaresExpr(:sum_squares, AffineExpr[], Float64[], num^2)
+end
+
 function sum_squares(affine::AffineExpr)
   affines = AffineExpr[]
   push!(affines, affine)
-  this = SumSquaresExpr(:sum_squares, affines, ones(length(affines)))
+  this = SumSquaresExpr(:sum_squares, affines, ones(length(affines)), 0.0)
 end
 
 function evaluate(x::SumSquaresExpr)
   return (x.evaluate())[1]
 end
 
-endof(x::AffineOrConstant) = x.size[1]*x.size[2]
+endof(x::AffineOrConstant) = x.size[1] * x.size[2]
 function size(x::AffineOrConstant, dim::Integer)
   if dim < 1
     error("dimension out of range")

src/types/problems.jl

@@ -6,11 +6,11 @@ SumSquaresExprOrNothing = Union(SumSquaresExpr, Nothing)
 type Problem
   head::Symbol
   objective::SumSquaresExpr
-  constraints::Array{EqConstraint}
+  constraints::Array{EqConstraint, 1}
   status::ASCIIString
   optval::Float64OrNothing
 
-  function Problem(head::Symbol, objective::SumSquaresExprOrNothing, constraints::Array{EqConstraint})
+  function Problem(head::Symbol, objective::SumSquaresExprOrNothing, constraints::Array{EqConstraint, 1})
     this = new(head, objective, constraints, "not yet solved", nothing)
     return this
   end

test/unit_tests/constraints.jl

@@ -6,4 +6,12 @@ x + y == 3
 rand(4, 3) == x - 4*y
 sum(y, 1) == rand(1, 3)
 
+# lists of constraints
+list = []
+list += x + y == ones(4, 3)
+
+list2 = EqConstraint[]
+list2 += y == 3
+list2 += list
+
 info("All equality constraints tests passed")