diff --git a/src/IntervalConstraintProgramming.jl b/src/IntervalConstraintProgramming.jl
index 463f989..1736530 100644
--- a/src/IntervalConstraintProgramming.jl
+++ b/src/IntervalConstraintProgramming.jl
@@ -10,6 +10,8 @@ using FixedSizeArrays: setindex
 import Base:
     show, ∩, ∪, !, ⊆, setdiff
 
+import ValidatedNumerics: sqr
+
 export
     @contractor,
     Separator, separator, @separator, @constraint,
@@ -17,8 +19,7 @@ export
     SubPaving, Paving,
     pave, refine!,
     Vol,
-    show_code,
-    plus_rev, mul_rev 
+    show_code
 
 
 include("reverse_mode.jl")
@@ -31,5 +32,9 @@ include("setinversion.jl")
 include("volume.jl")
 include("functions.jl")
 
+import Base.∪
+import ValidatedNumerics: IntervalBox
+∪{N,T}(X::IntervalBox{N,T}, Y::IntervalBox{N,T}) =
+    IntervalBox( [(x ∪ y) for (x,y) in zip(X, Y)]... )
 
 end # module
diff --git a/src/ast.jl b/src/ast.jl
index 67e0ba8..c4f8088 100644
--- a/src/ast.jl
+++ b/src/ast.jl
@@ -1,4 +1,3 @@
-
 const symbol_numbers = Dict{Symbol, Int}()
 
 doc"""Return a new, unique symbol like _z3_"""
@@ -11,10 +10,11 @@ function make_symbol(s::Symbol = :z)  # default is :z
 end
 
 
-function make_symbols(n::Integer)
-    [make_symbol() for i in 1:n]
+function make_symbols(n::Integer, s::Symbol = :z)
+    [make_symbol(s) for i in 1:n]
 end
 
+# The following function is not used
 doc"""Check if a symbol like `:a` has been uniqued to `:_a_1_`"""
 function isuniqued(s::Symbol)
     ss = string(s)
@@ -23,6 +23,8 @@ end
 
 # Types for representing a flattened AST:
 
+# Combine Assignment and FunctionAssignment ?
+
 immutable Assignment
     lhs
     op
@@ -35,88 +37,98 @@ immutable FunctionAssignment
     args
 end
 
-immutable GeneratedFunction
-    input_arguments::Vector{Symbol}
-    output_arguments::Vector{Symbol}
-    code::Expr
-end
-
 # Close to single assignment form
-type FlattenedAST
+type FlatAST
+    top  # topmost variable(s)
     input_variables::Set{Symbol}
     variables::Vector{Symbol}  # cleaned version
     intermediate::Vector{Symbol}  # generated vars
     code # ::Vector{Assignment}
 end
 
-import Base.show
-function show(io::IO, flatAST::FlattenedAST)
+
+function Base.show(io::IO, flatAST::FlatAST)
+    println(io, "top: ", flatAST.top)
     println(io, "input vars: ", flatAST.input_variables)
     println(io, "intermediate vars: ", flatAST.intermediate)
-    println(io, "code: ")
-    println(io, flatAST.code)
+    println(io, "code: ", flatAST.code)
 end
 
-FlattenedAST() = FlattenedAST(Set{Symbol}(), [], [], [])
+FlatAST() = FlatAST([], Set{Symbol}(), [], [], [])
 
-export FlattenedAST
+export FlatAST
 
 ##
 
+set_top!(flatAST::FlatAST, vars) = flatAST.top = vars  # also returns vars
+
+add_variable!(flatAST::FlatAST, var) = push!(flatAST.input_variables, var)
+
+add_intermediate!(flatAST::FlatAST, var::Symbol) = push!(flatAST.intermediate, var)
+add_intermediate!(flatAST::FlatAST, vars::Vector{Symbol}) = append!(flatAST.intermediate, vars)
+
+add_code!(flatAST::FlatAST, code) = push!(flatAST.code, code)
+
 export flatten
 
 
 function flatten(ex)
-    flatAST = FlattenedAST()
-    top_var = flatten!(flatAST, ex)
+    flatAST = FlatAST()
+    top = flatten!(flatAST, ex)
 
-    return top_var, flatAST
+    return top, flatAST
 end
 
 
 doc"""`flatten!` recursively converts a Julia expression into a "flat" (one-dimensional)
-structure, stored in a FlattenedAST object. This is close to SSA (single-assignment form,
- https://en.wikipedia.org/wiki/Static_single_assignment_form).
-
- Variables that are found are considered `input_variables`.
- Generated variables introduced at intermediate nodes are stored in
- `intermediate`.
- The function returns the variable that is
-at the top of the current piece of the tree."""
-# process numbers
-function flatten!(flatAST::FlattenedAST, ex)
+structure, stored in a FlatAST object. This is close to SSA (single-assignment form,
+https://en.wikipedia.org/wiki/Static_single_assignment_form).
+
+Variables that are found are considered `input_variables`.
+Generated variables introduced at intermediate nodes are stored in
+`intermediate`.
+Returns the variable at the top of the current piece of the tree."""
+
+# TODO: Parameters
+
+# numbers:
+function flatten!(flatAST::FlatAST, ex)
     return ex  # nothing to do to the AST; return the number
 end
 
-function flatten!(flatAST::FlattenedAST, ex::Symbol)  # symbols are leaves
-    push!(flatAST.input_variables, ex)  # add the discovered symbol as an input variable
+# symbols:
+function flatten!(flatAST::FlatAST, ex::Symbol)  # symbols are leaves
+    add_variable!(flatAST, ex)  # add the discovered symbol as an input variable
     return ex
 end
 
 
-function flatten!(flatAST::FlattenedAST, ex::Expr)
+function flatten!(flatAST::FlatAST, ex::Expr)
+    local top
 
     if ex.head == :$    # constants written as $a
-        process_constant!(flatAST, ex)
+        top = process_constant!(flatAST, ex)
 
     elseif ex.head == :call  # function calls
-        process_call!(flatAST, ex)
+        top = process_call!(flatAST, ex)
 
     elseif ex.head == :(=)  # assignments
-        process_assignment!(flatAST, ex)
+        top = process_assignment!(flatAST, ex)
 
     elseif ex.head == :block
-        process_block!(flatAST, ex)
+        top = process_block!(flatAST, ex)
 
     elseif ex.head == :tuple
-        process_tuple!(flatAST, ex)
+        top = process_tuple!(flatAST, ex)
 
     elseif ex.head == :return
-        process_return!(flatAST, ex)
+        top = process_return!(flatAST, ex)
     end
+
+    set_top!(flatAST, top)
 end
 
-function process_constant!(flatAST::FlattenedAST, ex)
+function process_constant!(flatAST::FlatAST, ex)
     return esc(ex.args[1])  # interpolate the value of the external constant
 end
 
@@ -125,7 +137,7 @@ end
 They are processed in order.
 """
 
-function process_block!(flatAST::FlattenedAST, ex)
+function process_block!(flatAST::FlatAST, ex)
     local top
 
     for arg in ex.args
@@ -136,9 +148,9 @@ function process_block!(flatAST::FlattenedAST, ex)
     return top  # last variable assigned
 end
 
-# function process_iterated_function!(flatAST::FlattenedAST, ex)
+# function process_iterated_function!(flatAST::FlatAST, ex)
 
-function process_tuple!(flatAST::FlattenedAST, ex)
+function process_tuple!(flatAST::FlatAST, ex)
     # println("Entering process_tuple")
     # @show flatAST
     # @show ex
@@ -147,8 +159,6 @@ function process_tuple!(flatAST::FlattenedAST, ex)
     top_args = []  # the arguments returned for each element of the tuple
     for arg in ex.args
         top = flatten!(flatAST, arg)
-        # @show flatAST
-
         push!(top_args, top)
     end
 
@@ -160,7 +170,7 @@ end
 The name a is currently retained.
 TODO: It should later be made unique.
 """
-function process_assignment!(flatAST::FlattenedAST, ex)
+function process_assignment!(flatAST::FlatAST, ex)
     # println("process_assignment!:")
     #  @show ex
     #  @show ex.args[1], ex.args[2]
@@ -185,10 +195,10 @@ function process_assignment!(flatAST::FlattenedAST, ex)
         vars = [var]
     end
 
-    append!(flatAST.intermediate, vars)
+    add_intermediate!(flatAST, vars)
 
     top_level_code = Assignment(vars, :(), top)  # empty operation
-    push!(flatAST.code, top_level_code)
+    add_code!(flatAST, top_level_code)
 
     # j@show flatAST
 
@@ -198,7 +208,7 @@ end
 
 """Processes something of the form `(f↑4)(x)` (write as `\\uparrow<TAB>`)
 by rewriting it to the equivalent set of iterated functions"""
-function process_iterated_function!(flatAST::FlattenedAST, ex)
+function process_iterated_function!(flatAST::FlatAST, ex)
     total_function_call = ex.args[1]
     args = ex.args[2:end]
 
@@ -217,7 +227,7 @@ function process_iterated_function!(flatAST::FlattenedAST, ex)
 
     # @show new_expr
 
-    flatten!(flatAST, new_expr)
+    flatten!(flatAST, new_expr)  # replace the current expression with the new one
 end
 
 """A call is something like +(x, y).
@@ -225,7 +235,7 @@ A new variable is introduced for the result; its name can be specified
     using the new_var optional argument. If none is given, then a new, generated
     name is used.
 """
-function process_call!(flatAST::FlattenedAST, ex, new_var=nothing)
+function process_call!(flatAST::FlatAST, ex, new_var=nothing)
 
     #println("Entering process_call!")
     #@show ex
@@ -245,17 +255,14 @@ function process_call!(flatAST::FlattenedAST, ex, new_var=nothing)
 
     end
 
-    # rewrite +(a,b,c) as +(a,+(b,c)):
+    # rewrite +(a,b,c) as +(a,+(b,c)) by recursive splitting
     # TODO: Use @match here!
 
-    # if op in (:+, :*) && length(ex.args) > 3
-    #     return insert_variables( :( ($op)($(ex.args[2]), ($op)($(ex.args[3:end]...) )) ))
-    # end
-
-    # new_code = quote end
-    # current_args = []  # the arguments in the current expression that will be added
-    # all_vars = Set{Symbol}()  # all variables contained in the sub-expressions
-    # generated_variables = Symbol[]
+    if op in (:+, :*) && length(ex.args) > 3
+        return flatten!(flatAST,
+            :( ($op)($(ex.args[2]), ($op)($(ex.args[3:end]...) )) )
+            )
+    end
 
     top_args = []
     for arg in ex.args[2:end]
@@ -281,7 +288,7 @@ function process_call!(flatAST::FlattenedAST, ex, new_var=nothing)
             new_var = make_symbol()
         end
 
-        push!(flatAST.intermediate, new_var)
+        add_intermediate!(flatAST, new_var)
 
         top_level_code = Assignment(new_var, op, top_args)
 
@@ -291,7 +298,7 @@ function process_call!(flatAST::FlattenedAST, ex, new_var=nothing)
             # make enough new variables for all the returned arguments:
             new_vars = make_symbols(length(registered_functions[op].generated))
 
-            append!(flatAST.intermediate, new_vars)
+            add_intermediate!(flatAST, new_vars)
 
             top_level_code = FunctionAssignment(new_vars, op, top_args)
 
@@ -305,7 +312,7 @@ function process_call!(flatAST::FlattenedAST, ex, new_var=nothing)
     end
 
 
-    push!(flatAST.code, top_level_code)
+    add_code!(flatAST, top_level_code)
 
     return new_var
 
diff --git a/src/code_generation.jl b/src/code_generation.jl
index 73fbbb0..9ee0a2d 100644
--- a/src/code_generation.jl
+++ b/src/code_generation.jl
@@ -37,19 +37,17 @@ function emit_forward_code(a::FunctionAssignment)
 end
 
 
-function emit_forward_code(code) #code::Vector{Assignment})
+function emit_forward_code(code)  # code::Vector{Assignment})
     new_code = quote end
     new_code.args = vcat([emit_forward_code(line) for line in code])
     return new_code
 end
 
 
-
 function emit_backward_code(a::Assignment)
 
     args = isa(a.args, Vector) ? a.args : [a.args]
 
-
     return_args = [a.lhs, args...]
     rev_op = rev_ops[a.op]  # find reverse operation
 
@@ -95,32 +93,41 @@ end
 
 
 
-function forward_pass(flatAST::FlattenedAST)
+function forward_backward(flatAST::FlatAST)
 
     # @show flatAST.input_variables
     # @show flatAST.intermediate
 
-    input_variables = sort(collect(flatAST.input_variables))
-    input_variables = setdiff(input_variables, flatAST.intermediate)  # remove local variables
-    flatAST.variables = input_variables  # why??
+    input = sort(collect(flatAST.input_variables))
 
-    generated_code = emit_forward_code(flatAST.code)
-    #make_function(input_variables, flatAST.intermediate, generated_code)
-    return GeneratedFunction(input_variables, flatAST.intermediate, generated_code)
-end
+    if isa(flatAST.top, Symbol)
+        output = [flatAST.top]
+    else
+        output = flatAST.top
+    end
+
+    #@show input
+    #@show flatAST.intermediate
+
+    input = setdiff(input, flatAST.intermediate)  # remove local variables
+    intermediate = setdiff(flatAST.intermediate, output)
+
+    flatAST.variables = input
 
-function backward_pass(flatAST::FlattenedAST)
+    code = emit_forward_code(flatAST.code)
+    forward = make_function(input, [output; intermediate], code)
 
-    generated_code = emit_backward_code(flatAST.code)
-    # make_function([flatAST.variables; flatAST.intermediate],
-    #                 flatAST.variables,
-    #                 generated_code)
-    # # reverse input_variables and intermediate?
+    code = emit_backward_code(flatAST.code)
+    backward = make_function([input; output; intermediate],
+                                input, code)
 
-    all_variables = [flatAST.variables; flatAST.intermediate]
-    return GeneratedFunction(all_variables,
-                            flatAST.variables,
-                            generated_code)
+# @show input
+# @show output
+# @show intermediate
+
+    # return GeneratedFunction(input, output, intermediate, code)
+
+    return (forward, backward)
 end
 
 
@@ -133,16 +140,10 @@ function make_function(input_args, output_args, code)
     input = make_tuple(input_args)  # make a tuple of the variables
     output = make_tuple(output_args)  # make a tuple of the variables
 
-    return quote
-                $input -> begin
-                            $code
-                            return $output
-                        end
-                end
-
-
-
+    quote
+        $input -> begin
+                    $code
+                    return $output
+                  end
+        end
 end
-
-make_function(f::GeneratedFunction) =
-        make_function(f.input_arguments, f.output_arguments, f.code)
diff --git a/src/contractor.jl b/src/contractor.jl
index 862321d..89b2b75 100644
--- a/src/contractor.jl
+++ b/src/contractor.jl
@@ -1,67 +1,20 @@
 
 
-type Contractor{F<:Function}
-    variables::Vector{Symbol}
-    constraint_expression::Expr
-    code::Expr
-    contractor::F  # function
+immutable Contractor{F1<:Function, F2<:Function}
+    variables::Vector{Symbol}  # input variables
+    num_outputs::Int
+    forward::F1
+    backward::F2
+    forward_code::Expr
+    backward_code::Expr
 end
 
 
-doc"""`parse_comparison` parses comparisons like `x >= 10`
-into the corresponding interval, expressed as `x ∈ [10,∞]`
-
-Returns the expression and the constraint interval
-
-TODO: Allow something like [3,4]' for the complement of [3,4]
-"""
-
-function parse_comparison(ex)
-    expr, limits =
-    @match ex begin
-       ((a_ <= b_) | (a_ < b_) | (a_ ≤ b_))   => (a, (-∞, b))
-       ((a_ >= b_) | (a_ > b_) | (a_ ≥ b_))   => (a, (b, ∞))
-
-       ((a_ == b_) | (a_ = b_))   => (a, (b, b))
-
-       ((a_ <= b_ <= c_)
-        | (a_ < b_ < c_)
-        | (a_ <= b_ < c)
-        | (a_ < b_ <= c))         => (b, (a, c))
-
-       ((a_ >= b_ >= c_)
-       | (a_ > b_ > c_)
-       | (a_ >= b_ > c_)
-       | (a_ > b_ >= c))          => (b, (c, a))
-
-       ((a_ ∈ [b_, c_])
-       | (a_ in [b_, c_])
-       | (a_ ∈ b_ .. c_)
-       | (a_ in b_ .. c_))        => (a, (b, c))
-
-       _                          => (ex, (-∞, ∞))
-
-   end
-
-   a, b = limits
-
-   return (expr, a..b)   # expr ∈ [a,b]
-
-end
-
-
-# new call syntax to define a "functor" (object that behaves like a function)
-@compat (C::Contractor)(x...) = C.contractor(x...)
-
-@compat (C::Contractor)(X::IntervalBox) = C.contractor(X...)
-#show_code(c::Contractor) = c.code
-
-
-function Base.show(io::IO, C::Contractor)
-    println(io, "Contractor:")
-    println(io, "  - variables: $(C.variables)")
-    print(io, "  - constraint: $(C.constraint_expression)")
-end
+# function Base.show(io::IO, C::Contractor)
+#     println(io, "Contractor:")
+#     println(io, "  - variables: $(C.variables)")
+#     print(io, "  - constraint: $(C.constraint_expression)")
+# end
 
 doc"""Usage:
 ```
@@ -77,15 +30,25 @@ TODO: Hygiene for global variables, or pass in parameters
 """
 
 macro contractor(ex)
-    # println("@contractor; ex=$ex")
-
     make_contractor(ex)
 end
 
 
+@compat function (C::Contractor{F1,F2}){F1,F2}(A, X) # X::IntervalBox)
+    z = IntervalBox( C.forward(IntervalBox(X...)...)... )
+    #z = [1:C.num_outputs] = tuple(IntervalBox(z[1:C.num_outputs]...) ∩ A
 
+    # @show z
+    constrained = IntervalBox(z[1:C.num_outputs]...) ∩ IntervalBox(A...)
+    #@show constrained
+    #@show z[(C.num_outputs)+1:end]
+    return IntervalBox( C.backward( X...,
+                                    constrained...,
+                                    z[(C.num_outputs)+1:end]...
+                                  )...
+                       )
+end
 
-#function Contractor(ex::Expr)
 function make_contractor(ex::Expr)
     # println("Entering Contractor(ex) with ex=$ex")
     expr, constraint_interval = parse_comparison(ex)
@@ -98,82 +61,15 @@ function make_contractor(ex::Expr)
 
     # @show top, linear_AST
 
-    forward = forward_pass(linear_AST)
-    backward = backward_pass(linear_AST)
-
-
-
-    # TODO: What about interval box constraints?
-    input_arguments = forward.input_arguments
-    augmented_input_arguments = [:_A_; forward.input_arguments]
-
-    # @show input_arguments
-    # @show augmented_input_arguments
-
-    # add constraint interval as first argument
-
-    input_variables = make_tuple(input_arguments)
-    augmented_input_variables = make_tuple(augmented_input_arguments)
-
-    forward_output = make_tuple(forward.output_arguments)
-
-    backward_output = make_tuple(backward.output_arguments)
-
-    # @show forward
-    # @show backward
-    #
-    # @show input_variables
-    # @show forward_output
-    # @show backward_output
-
-    if isa(top, Symbol)
-        # nothing
-    elseif length(top) == 1  # single variable
-        top = top[]
-
-    else
-        # TODO: implement what happens for multiple variables in the constraint
-        # using an IntervalBox and intersection of IntervalBoxes
-    end
-
-    top_args = make_tuple(top)
-
-    local intersect_code
-
-    if isa(top_args, Symbol)
-        intersect_code = :($top_args = $top_args ∩ _A_)  # check type stability
-    else
-        intersect_code = :($top_args = IntervalBox($top_args) ∩ _A_)  # check type stability
-    end
-
-
-    code =
-        #esc(quote
-        quote
-            $(augmented_input_variables) -> begin
-                forward = $(make_function(forward))
-                backward = $(make_function(backward))
-
-                $(forward_output) = forward($(forward.input_arguments...))
-
-                $intersect_code
-
-                $(backward_output) = backward($(backward.input_arguments...))
-
-                return $(input_variables)
-
-            end
+    forward, backward  = forward_backward(linear_AST)
 
-        end
+    num_outputs = isa(linear_AST.top, Symbol) ? 1 : length(linear_AST.top)
 
-    #  @show forward
-    #  @show backward
-    # #
-    #  @show code
+    :(Contractor($linear_AST.variables,
+                    $num_outputs,
+                    $forward,
+                    $backward,
+                    $(Meta.quot(forward)),
+                    $(Meta.quot(backward))))
 
-    return :(Contractor($(augmented_input_arguments),
-                        $(Meta.quot(expr)),
-                        $(Meta.quot(code)),
-                        $(code)
-                        ))
 end
diff --git a/src/functions.jl b/src/functions.jl
index d54ffe1..b028b82 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -1,21 +1,4 @@
-#=
 
-Want to process `@constraint f(f(x)) ∈ [0.3, 0.4]`
-where `f(x) = 4x * (1-x)`
-
-Given code `f(f(x))`, we need `f_forward` and `f_backward`.
-
-Each "copy" of `f` uses the same actual forward and back functions,
-`f.forward` and `f.backward`.
-
-```
-@function f(x) = 4x * (1-x)
-```
-should generate these forward and backward functions, and register the function
-`f`.
-
-"""
-=#
 
 doc"""
 A `ConstraintFunction` contains the created forward and backward
@@ -53,24 +36,10 @@ Example: `@function f(x, y) = x^2 + y^2`
 
 @eval macro ($(:function))(ex)   # workaround to define macro @function
 
-    # (f, args, code) = @match ex begin
-    #     ( f_(args__) = code_ ) => (f, args, code)
-    # end
-
     (f, args, code) = match_function(ex)
 
-    # @show f, args, code
-
     return_arguments, flatAST = flatten(code)
 
-    # @show return_arguments
-
-    # @show root, all_vars, generated, code2
-
-
-    # @show flatAST
-    # @show return_arguments
-
     # make into an array:
     if !(isa(return_arguments, Array))
         return_arguments = [return_arguments]
@@ -84,34 +53,23 @@ Example: `@function f(x, y) = x^2 + y^2`
     flatAST.intermediate = [return_arguments; flatAST.intermediate]
 
 
-    forward_code = forward_pass(flatAST) #root, all_vars, generated, code2)
-    backward_code = backward_pass(flatAST) #root, all_vars, generated, code2)
-
-    # @show forward_code, backward_code
+    registered_functions[f] = FunctionArguments(flatAST.variables, flatAST.intermediate, return_arguments)
 
-    # @show make_function(forward_code)
-    # @show make_function(backward_code)
+    forward, backward = forward_backward(flatAST)
 
-    registered_functions[f] = FunctionArguments(flatAST.variables, flatAST.intermediate, return_arguments)
 
-    forward_function = make_function(forward_code)
-    backward_function = make_function(backward_code)
 
     return quote
-        #$(esc(Meta.quot(f))) = ConstraintFunction($(all_vars), $(generated), $(forward_code), $(backward_code))
-        #$(esc(f)) =
+
         $(esc(f)) =
             ConstraintFunction($(flatAST.variables),
                                 $(flatAST.intermediate),
-                                $(forward_function),
-                                $(backward_function),
-                                $(Meta.quot(forward_function)),
-                                $(Meta.quot(backward_function))
+                                $(forward),
+                                $(backward),
+                                $(Meta.quot(forward)),
+                                $(Meta.quot(backward))
                                 )
 
-
-        #registered_functions[$(Meta.quot(f))] =  ConstraintFunction($(all_vars), $(generated), $(forward_code), $(backward_code))
-        #$(Meta.quot(f)) =  ConstraintFunction($(all_vars), $(generated), $(forward_code), $(backward_code))
     end
 
 end
@@ -120,10 +78,9 @@ end
 function match_function(ex)
 
     try
-        (f, args, body) =
-            @match ex begin
-             ( (f_(args__) = body_) |
-              (function f_(args__) body_ end) ) => (f, args, body)
+        @capture ex begin
+            (   (f_(args__) = body_)
+                  | (function f_(args__) body_ end) )
            end
 
          return (f, args, rmlines(body))  # rmlines is from MacroTools package
diff --git a/src/paving.jl b/src/paving.jl
index 50189a0..34970c3 100644
--- a/src/paving.jl
+++ b/src/paving.jl
@@ -1,7 +1,7 @@
 typealias SubPaving{N,T} Vector{IntervalBox{N,T}}
 
 type Paving{N,T}
-    separator::Separator
+    separator::Separator   # parametrize!
     inner::SubPaving{N,T}
     boundary::SubPaving{N,T}
     ϵ::Float64
diff --git a/src/reverse_mode.jl b/src/reverse_mode.jl
index 4321c97..ba13f33 100644
--- a/src/reverse_mode.jl
+++ b/src/reverse_mode.jl
@@ -1,18 +1,22 @@
 export plus_rev, minus_rev, mul_rev,
     power_rev, sqrt_rev, sqr_rev # export due to quoting issue
 
+# export sqr
+# sqr(x) = x^2
+
 const rev_ops = Dict(
                     :+     => :plus_rev,
+                    :-     => :minus_rev,
                     :*     => :mul_rev,
                     :^     => :power_rev,
-                    :-     => :minus_rev,
                     :sqrt  => :sqrt_rev,
+                    :sqr   => :sqr_rev,
                     :()    => :()
                     )
 
 
 function plus_rev(a::Interval, b::Interval, c::Interval)  # a = b + c
-    # a = a ∩ (b + c)
+    # a = a ∩ (b + c)  # add this line for plus contractor (as opposed to reverse function)
     b = b ∩ (a - c)
     c = c ∩ (a - b)
 
@@ -94,7 +98,7 @@ function sqr_rev(c, x)   # c = x^2;  refine x
     x1 = sqrt(c) ∩ x
     x2 = -(sqrt(c)) ∩ x
 
-    return hull(x1, x2)
+    return (c, hull(x1, x2))
 end
 
 sqr_rev(c) = sqr_rev(c, -∞..∞)
diff --git a/src/separator.jl b/src/separator.jl
index 7e12cd8..c16d1da 100644
--- a/src/separator.jl
+++ b/src/separator.jl
@@ -5,21 +5,19 @@ doc"""
 ConstraintSeparator is a separator that represents a constraint defined directly
 using `@constraint`.
 """
-# CHANGE TO IMMUTABLE AND PARAMETRIZE THE FUNCTION FOR EFFICIENCY
-type ConstraintSeparator{C, II} <: Separator
+immutable ConstraintSeparator{C, II} <: Separator
     variables::Vector{Symbol}
-    #separator::Function
     constraint::II  # Interval or IntervalBox
     contractor::C
     expression::Expr
 end
 
-ConstraintSeparator(constraint, contractor, expression) = ConstraintSeparator(contractor.variables[2:end], constraint, contractor, expression)
+ConstraintSeparator(constraint, contractor, expression) = ConstraintSeparator(contractor.variables, constraint, contractor, expression)
 
 doc"""CombinationSeparator is a separator that is a combination (union, intersection,
 or complement) of other separators.
 """
-type CombinationSeparator{F} <: Separator
+immutable CombinationSeparator{F} <: Separator
     variables::Vector{Symbol}
     separator::F
     expression::Expr
@@ -29,25 +27,25 @@ end
     C = S.contractor
     a, b = S.constraint.lo, S.constraint.hi
 
-    inner = C(a..b, X...)  # closure over the function C
+    inner = C(a..b, X)
 
     local outer
 
     if a == -∞
-        outer = C(b..∞, X...)
+        outer = C(b..∞, X)
 
     elseif b == ∞
-        outer = C(-∞..a, X...)
+        outer = C(-∞..a, X)
 
     else
 
-        outer1 = C(-∞..a, X...)
-        outer2 = C(b..∞, X...)
+        outer1 = C(-∞..a, X)
+        outer2 = C(b..∞, X)
 
-        outer = [ hull(x1, x2) for (x1,x2) in zip(outer1, outer2) ]
+        outer = outer1 ∪ outer2
     end
 
-    return (inner, (outer...))
+    return (inner, outer)
 end
 
 
@@ -135,32 +133,25 @@ function make_constraint(expr, constraint)
     code
 end
 
-macro constraint(ex::Expr)  # alternative name for constraint -- remove?
-    # @show ex
-    expr, constraint = parse_comparison(ex)
+doc"""Create a separator from a given constraint expression, written as
+standard Julia code.
+
+e.g. `C = @constraint x^2 + y^2 <= 1`
 
+The variables (`x` and `y`, in this case) are automatically inferred.
+External constants can be used as e.g. `$a`:
+
+```
+a = 3
+C = @constraint x^2 + y^2 <= $a
+```
+"""
+macro constraint(ex::Expr)
+    expr, constraint = parse_comparison(ex)
     make_constraint(expr, constraint)
 end
 
 
-# doc"""Create a separator from a given constraint expression, written as
-# standard Julia code.
-#
-# e.g. `C = @constraint x^2 + y^2 <= 1`
-#
-# The variables (`x` and `y`, in this case) are automatically inferred.
-# External constants can be used as e.g. `$a`:
-#
-# ```
-# a = 3
-# C = @constraint x^2 + y^2 <= $a
-# ```
-# """
-# macro constraint(ex::Expr)
-#     ex = Meta.quot(ex)
-#     :(ConstraintSeparator($ex))
-# end
-
 function show(io::IO, S::Separator)
     println(io, "Separator:")
     print(io, "- variables: ")
@@ -170,15 +161,9 @@ function show(io::IO, S::Separator)
     println(io, S.expression)
 end
 
-# show_code(S::ConstraintSeparator) = show_code(S.contractor)
 
-#@compat (S::ConstraintSeparator)(X) = S.separator(X)
 @compat (S::CombinationSeparator)(X) = S.separator(X)
 
-#@compat (S::ConstraintSeparator)(X) = S.separator(X)
-
-# show_code(S::Separator) = show_code(S.contractor)
-
 
 doc"Unify the variables of two separators"
 function unify_variables(vars1, vars2)
@@ -252,8 +237,8 @@ function ∩(S1::Separator, S2::Separator)
 
         # Treat as if had X[i] in the other directions, except if empty
 
-        inner = tuple( [x ∩ y for (x,y) in zip(inner1, inner2) ]... )
-        outer = tuple( [x ∪ y for (x,y) in zip(outer1, outer2) ]... )
+        inner = IntervalBox( [x ∩ y for (x,y) in zip(inner1, inner2) ]... )
+        outer = IntervalBox( [x ∪ y for (x,y) in zip(outer1, outer2) ]... )
 
         return (inner, outer)
 
@@ -299,8 +284,8 @@ function ∪(S1::Separator, S2::Separator)
         end
 
 
-        inner = tuple( [x ∪ y for (x,y) in zip(inner1, inner2) ]... )
-        outer = tuple( [x ∩ y for (x,y) in zip(outer1, outer2) ]... )
+        inner = IntervalBox( [x ∪ y for (x,y) in zip(inner1, inner2) ]... )
+        outer = IntervalBox( [x ∩ y for (x,y) in zip(outer1, outer2) ]... )
 
         return (inner, outer)
     end
diff --git a/src/setinversion.jl b/src/setinversion.jl
index 81e340a..a896cd1 100644
--- a/src/setinversion.jl
+++ b/src/setinversion.jl
@@ -37,17 +37,15 @@ function pave{N,T}(S::Separator, working::Vector{IntervalBox{N,T}}, ϵ)
         inner, outer = S(X)   # here inner and outer are reversed compared to Jaulin
         # S(X) returns the pair (contractor with respect to the inside of the constraing, contractor with respect to outside)
 
-        inner2 = IntervalBox(inner)
-        outer2 = IntervalBox(outer)
-
-        inside_list = setdiff(X, outer2)
+        #@show X, outer
+        inside_list = setdiff(X, outer)
 
         if length(inside_list) > 0
             append!(inner_list, inside_list)
         end
 
 
-        boundary = inner2 ∩ outer2
+        boundary = inner ∩ outer
 
         if isempty(boundary)
             continue
@@ -77,7 +75,7 @@ function pave{N,T}(S::Separator, X::IntervalBox{N,T}, ϵ = 1e-2)
 
     inner_list, boundary_list = pave(S, [X], ϵ)
 
-    return Paving{N,T}(S, inner_list, boundary_list, ϵ)
+    return Paving(S, inner_list, boundary_list, ϵ)
 
 end
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 09f35c0..11e6e2f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -22,7 +22,7 @@ end
 
     C = @contractor x^2 + y^2
 
-    @test C(-∞..1, x, y) == (-1..1, -1..1)
+    @test C(-∞..1, X) == IntervalBox(-1..1, -1..1)
 end
 
 @testset "Separators" begin
@@ -33,13 +33,13 @@ end
     @test typeof(S) <: IntervalConstraintProgramming.ConstraintSeparator
 
     inner, outer = S(X)
-    @test inner == (-1..1, -1..1)
-    @test outer == (II, II)
+    @test inner == IntervalBox(-1..1, -1..1)
+    @test outer == IntervalBox(II, II)
 
     X = IntervalBox(-∞..∞, -∞..∞)
     inner, outer = S(X)
-    @test inner == (-1..1, -1..1)
-    @test outer == (-∞..∞, -∞..∞)
+    @test inner == IntervalBox(-1..1, -1..1)
+    @test outer == IntervalBox(-∞..∞, -∞..∞)
 end
 
 @testset "pave" begin
@@ -60,7 +60,7 @@ end
 
     S3 = @constraint x^2 + y^2 <= 1
     X = IntervalBox(-∞..∞, -∞..∞)
-    paving = pave(S3, X, 1)
+    paving = pave(S3, X, 1.0)
 
     @test paving.inner == [IntervalBox(Interval(0.0, 0.5), Interval(0.0, 0.8660254037844386)),
                     IntervalBox(Interval(0.0, 0.5), Interval(-0.8660254037844386, 0.0)),
@@ -110,7 +110,7 @@ end
     A = 0.5..1
     x = 0..1
 
-    @test C1(A, x) == 0.125..0.25   # x such that 4x ∈ A=[0.5, 1]
+    @test C1(A, x) == IntervalBox(0.125..0.25)   # x such that 4x ∈ A=[0.5, 1]
 
 
     C2 = @constraint f(x) ∈ [0.5, 0.6]
@@ -137,12 +137,12 @@ end
     A = 0.5..1
     x = 0..1
 
-    @test C(A, x) == sqrt(A / 4)
+    @test C(A, x) == IntervalBox(sqrt(A / 4))
 
     @function g2(x) = ( a = f(f(x)); a^2 )
     C2 = @contractor g2(x)
 
-    @test C2(A, x) == sqrt(A / 16)
+    @test C2(A, x) == IntervalBox(sqrt(A / 16))
 
 end
 
@@ -156,23 +156,5 @@ end
     @function g3(x) = ( a = (f↑2)(x); a^2 )
     C3 = @contractor g3(x)
 
-    @test C3(A, x) == sqrt(A / 16)
-end
-
-@testset "Multidimensional functions" begin
-
-    @function g4(x, y) = (2x, 2y)
-
-    A = IntervalBox(0.5..1, 0.5..1)
-    x = y = 0..1
-
-    C4 = @contractor g4(x, y)
-
-    @test IntervalBox(C4(A, x, y)) == A / 2
-
-
-    C5 = @contractor (g4↑2)(x, y)
-
-    @test IntervalBox(C5(A, x, y)) == A / 4
-
+    @test C3(A, x) == IntervalBox(sqrt(A / 16))
 end