Merge branch 'teh/constrained' of git://github.com/JuliaOpt/Optim.jl …

…into teh/constrained

src/cg.jl

@@ -148,7 +148,9 @@ function cg{T}(df::Union(DifferentiableFunction,
     y = similar(x)
 
     # Store f(x) in f_x
-    f_x_previous, f_x = NaN, df.fg!(x, gr)
+    f_x = df.fg!(x, gr)
+    @assert typeof(f_x) == T
+    f_x_previous = nan(T)
     f_calls, g_calls = f_calls + 1, g_calls + 1
     copy!(gr_previous, gr)
 
@@ -206,6 +208,8 @@ function cg{T}(df::Union(DifferentiableFunction,
 
         # Refresh the line search cache
         clear!(lsr)
+        @assert typeof(f_x) == T
+        @assert typeof(dphi0) == T
         push!(lsr, zero(T), f_x, dphi0)
 
         alphamax = interior ? toedge(x, s, constraints) : inf(T)

src/constraints.jl

@@ -15,8 +15,8 @@ immutable ConstraintsBox{T,N} <: AbstractConstraints
     end
 end
 ConstraintsBox{T,N}(l::AbstractArray{T,N}, u::AbstractArray{T,N}) = ConstraintsBox{T,N}(l, u)
-ConstraintsBox{T,N}(l::AbstractArray{T,N}, u::Nothing) = ConstraintsBox{T,N}(l, infs(T, size(l)))
-ConstraintsBox{T,N}(l::Nothing, u::AbstractArray{T,N}) = ConstraintsBox{T,N}(-infs(T,size(u)), u)
+ConstraintsBox{T,N}(l::AbstractArray{T,N}, u::Nothing) = ConstraintsBox{T,N}(l, fill(inf(T), size(l)))
+ConstraintsBox{T,N}(l::Nothing, u::AbstractArray{T,N}) = ConstraintsBox{T,N}(fill(-inf(T),size(u)), u)
 
 immutable ConstraintsL{T,M<:AbstractMatrix,N} <: AbstractConstraints
     A::M

src/interior.jl

@@ -24,7 +24,7 @@ function barrier_f{T}(x::Array{T}, bounds::ConstraintsBox)
             phi += log(ui-xi)
         end
     end
-    -phi
+    -phi::T
 end
 
 function barrier_g!{T}(x::Array{T}, g::Array{T}, bounds::ConstraintsBox)
@@ -58,7 +58,7 @@ function barrier_fg!{T}(x::Array{T}, g::Array{T}, bounds::ConstraintsBox)
             g[i] += 1/(ui-xi)
         end
     end
-    -phi
+    -phi::T
 end
 
 function barrier_h!{T}(x::Vector{T}, H::Matrix{T}, bounds::ConstraintsBox)
@@ -116,7 +116,7 @@ function barrier_fg!{T}(x::Array{T}, g::Array{T}, constraints::ConstraintsL)
     for i = 1:lenx
         g[i] += y3[i]
     end
-    -phi
+    -phi::T
 end
 
 barrier_f{T}(x::Array{T}, constraints::ConstraintsL) = barrier_fg!(x, T[], constraints)
@@ -322,7 +322,7 @@ function interior{T}(objective::Union(DifferentiableFunction, TwiceDifferentiabl
         if gc == 0
             t = one(T)  # FIXME: bounds constraints, starting at exact center. This is probably not the right guess.
         else
-            t = gc/(0.1*go)
+            t = gc/(convert(T,0.1)*go)
         end
     end
     if method == :newton

src/linesearch/hz_linesearch.jl

@@ -120,21 +120,22 @@ function alphatry{T}(alpha::T,
         lsr.nfailures += 1
         iterfinite += 1
         if iterfinite >= iterfinitemax
-            error("Failed to achieve finite test value")
+            return zero(T), true, f_calls, g_calls
+#             error("Failed to achieve finite test value; alphatest = ", alphatest)
         end
     end
-    a = (phitest - alphatest * dphi0 - phi0) / alphatest^2 # quadratic fit
+    a = ((phitest-phi0)/alphatest - dphi0)/alphatest  # quadratic fit
     if detailed_trace & ALPHAGUESS > 0
         println("quadfit: alphatest = ", alphatest,
                 ", phi0 = ", phi0,
                 ", phitest = ", phitest,
                 ", quadcoef = ", a)
     end
     mayterminate = false
-    if a > 0 && phitest <= phi0
+    if isfinite(a) && a > 0 && phitest <= phi0
         alpha = -dphi0 / 2 / a # if convex, choose minimum of quadratic
         if alpha == 0
-            error("alpha is zero. dphi0 = ", dphi0, ", a = ", a)
+            error("alpha is zero. dphi0 = ", dphi0, ", phi0 = ", phi0, ", phitest = ", phitest, ", alphatest = ", alphatest, ", a = ", a)
         end
         if alpha <= alphamax
             mayterminate = true
@@ -180,6 +181,9 @@ function hz_linesearch!{T}(df::Union(DifferentiableFunction,
                            psi3::Real = convert(T,0.1),
                            iterfinitemax::Integer = iceil(-log2(eps(T))),
                            detailed_trace::Integer = 0)
+    if detailed_trace & LINESEARCH > 0
+        println("New linesearch")
+    end
     f_calls = 0
     g_calls = 0
 
@@ -190,7 +194,6 @@ function hz_linesearch!{T}(df::Union(DifferentiableFunction,
     philim = phi0 + epsilon * abs(phi0)
     @assert c > 0
     @assert isfinite(c) && c <= alphamax
-    # phic = phi(gphi, c)
     phic, dphic, f_up, g_up = linefunc!(df, x, s, c, xtmp, g, true, constraints)
     f_calls += f_up
     g_calls += g_up
@@ -200,7 +203,6 @@ function hz_linesearch!{T}(df::Union(DifferentiableFunction,
         lsr.nfailures += 1
         iterfinite += 1
         c *= psi3
-        # phic = phi(gphi, c) # TODO: Fix
         phic, dphic, f_up, g_up = linefunc!(df, x, s, c, xtmp, g, true, constraints)
         f_calls += f_up
         g_calls += g_up
@@ -263,17 +265,20 @@ function hz_linesearch!{T}(df::Union(DifferentiableFunction,
             cold = c
             c *= rho
             if c > alphamax
+                c = (alphamax + cold)/2
                 if detailed_trace & BRACKET > 0
-                    println("bracket: exceeding alphamax, bisecting")
+                    println("bracket: exceeding alphamax, bisecting: alphamax = ", alphamax,
+                            ", cold = ", cold, ", new c = ", c)
+                end
+                if c == cold || nextfloat(c) >= alphamax
+                    return cold, f_calls, g_calls
                 end
-                c = (alphamax + cold)/2
             end
-            # phic = phi(gphi, c) # TODO: Replace
             phic, dphic, f_up, g_up = linefunc!(df, x, s, c, xtmp, g, true, constraints)
             f_calls += f_up
             g_calls += g_up
             iterfinite = 1
-            while !isfinite(phic) && c > cold && iterfinite < iterfinitemax
+            while !isfinite(phic) && c > nextfloat(cold) && iterfinite < iterfinitemax
                 alphamax = c
                 lsr.nfailures += 1
                 iterfinite += 1
@@ -285,7 +290,9 @@ function hz_linesearch!{T}(df::Union(DifferentiableFunction,
                 f_calls += f_up
                 g_calls += g_up
             end
-            if (dphic < 0 && c == alphamax) || !isfinite(phic)
+            if !isfinite(phic)
+                return cold, f_calls, g_calls
+            elseif dphic < 0 && c == alphamax
                 # We're on the edge of the allowed region, and the
                 # value is still decreasing. This can be due to
                 # roundoff error in barrier penalties, a barrier
@@ -299,23 +306,7 @@ function hz_linesearch!{T}(df::Union(DifferentiableFunction,
                             ", phic = ", phic,
                             ", dphic = ", dphic)
                 end
-                ic = length(lsr)
-                while !isfinite(phic)
-                    ic -= 1
-                    c = lsr.alpha[ic]
-                    phic = lsr.value[ic]
-                    if isfinite(phic)
-                        println("Using c = ", c, ", phic = ", phic)
-                    end
-                    # Re-evaluate at current position. This is important if
-                    # reportfunc makes use of cached storage, and that cache
-                    # has been corrupted by NaN/Inf
-                    # phic = phi(gphi, c) # TODO: Replace
-                    phic, dphic, f_up, g_up = linefunc!(df, x, s, c, xtmp, g, true, constraints)
-                    f_calls += f_up
-                    g_calls += g_up
-                end
-                return c, f_calls, g_calls # phic
+                return c, f_calls, g_calls
             end
             push!(lsr, c, phic, dphic)
         end
@@ -332,7 +323,6 @@ function hz_linesearch!{T}(df::Union(DifferentiableFunction,
                     ", b = ", b,
                     ", phi(a) = ", lsr.value[ia],
                     ", phi(b) = ", lsr.value[ib])
-            println("phi(a) == phi(b) ? ", lsr.value[ia] == lsr.value[ib])
         end
         if b - a <= eps(b)
             return a, f_calls, g_calls # lsr.value[ia]
@@ -378,6 +368,11 @@ function hz_linesearch!{T}(df::Union(DifferentiableFunction,
         end
         iter += 1
     end
+    @show ia, ib
+    @show lsr.alpha[ia], lsr.alpha[ib]
+    @show lsr.value[ia], lsr.value[ib]
+    @show lsr.slope[ia], lsr.slope[ib]
+    @show alphamax, iter
     error("Linesearch failed to converge")
 end
 
@@ -428,8 +423,8 @@ function secant2!{T}(df::Union(DifferentiableFunction,
     b = lsr.alpha[ib]
     dphia = lsr.slope[ia]
     dphib = lsr.slope[ib]
-    @assert dphia < 0
-    @assert dphib >= 0
+    @assert dphia < 0 && isfinite(dphia)
+    @assert dphib >= 0 && isfinite(dphib)
     c = secant(a, b, dphia, dphib)
     if detailed_trace & SECANT2 > 0
         println("secant2: a = ", a, ", b = ", b, ", c = ", c)
@@ -439,7 +434,6 @@ function secant2!{T}(df::Union(DifferentiableFunction,
     phic, dphic, f_up, g_up = linefunc!(df, x, s, c, xtmp, g, true, constraints)
     f_calls += f_up
     g_calls += g_up
-    @assert isfinite(phic)
     push!(lsr, c, phic, dphic)
     ic = length(lsr)
     if satisfies_wolfe(c, phic, dphic, phi0, dphi0, philim, delta, sigma)
@@ -569,7 +563,12 @@ function bisect!{T}(df::Union(DifferentiableFunction,
     # Debugging (HZ, conditions shown following U3)
     @assert lsr.slope[ia] < 0
     @assert lsr.value[ia] <= philim
-    @assert lsr.slope[ib] < 0       # otherwise we wouldn't be here
+    # if !(lsr.slope[ib] < 0)
+    #     @show a, b
+    #     @show lsr.value[ia], lsr.value[ib]
+    #     @show lsr.slope[ia], lsr.slope[ib]
+    # end
+    # @assert lsr.slope[ib] < 0       # otherwise we wouldn't be here
     @assert lsr.value[ib] > philim
     @assert b > a
     while b - a > eps(b)
@@ -580,7 +579,6 @@ function bisect!{T}(df::Union(DifferentiableFunction,
         phid, gphi, f_up, g_up = linefunc!(df, x, s, d, xtmp, g, true, constraints)
         f_calls += f_up
         g_calls += g_up
-        @assert isfinite(phid)
         push!(lsr, d, phid, gphi)
         id = length(lsr)
         if gphi >= 0