NumericalOptimizationForFuncDesignerModels

from FuncDesigner import *
from openopt import LLSP
 
# create some variables
a, b, c = oovars('a', 'b', 'c')
# or just a, b, c = oovars(3)
 
# start point is unused by lapack_dgelss and lapack_dgelss
# but it is required to set dimensions of variables
# also, it is used for some converters e.g. r = p.solve('nlp:ralg')
startPoint = {a:0, b:0, c:0} 
 
# overdetermined system of 4 linear equations with 3 variables
# you can use "for" cycle, operations of sum() eg 
# startPoint = {a:[0,10], b:0, c:0} 
# f = [sum(a) + i * b + 4 * sqrt(i) * c - 2 * i**2 for i in xrange(40)]
f = [2*a+3*b-4*c+5, 2*a+13*b+15, a+4*b+4*c-25, 20*a+30*b-4*c+50]
 
# assign prob
p = LLSP(f, startPoint)
 
# solve
r = p.solve('lapack_dgelss')
 
# print result
print r.xf
 
# Expected output:
# {a: array([-0.3091145]), b: array([-0.86376906]), c: array([ 4.03827441])}

# Example of creating LP in FuncDesigner
# and solving it via OpenOpt
 
from FuncDesigner import *
from openopt import LP
 
# Define some oovars
x, y, z = oovars(3)
 
# Let's define some linear functions
f1 = 4*x+5*y + 3*z + 5
f2 = f1.sum() + 2*x + 4*y + 15
f3 = 5*f1 + 4*f2 + 20
 
# Define objective; sum(a) and a.sum() are same as well as for numpy arrays
obj = x.sum() + y - 50*z + sum(f3) + 2*f2.sum() + 4064.6
 
# Define some constraints via Python list or tuple or set of any length, probably via while/for cycles
constraints = [x+5*y<15, x[0]<4, f1<[25, 35], f1>-100, 2*f1+4*z<[80, 800], 5*f2+4*z<100, -5<x,  x<1, -20<y,  y<20, -4000<z, z<4]
 
# Start point - currently matters only size of variables, glpk, lpSolve and cvxopt_lp use start val = all-zeros
startPoint = {x:[8, 15], y:25, z:80} # however, using numpy.arrays is more recommended than Python lists
 
# Create prob
p = LP(obj, startPoint, constraints = constraints)
 
# Solve
r = p.solve('glpk') # glpk is name of solver involved, see OOF doc for more arguments
 
# Decode solution
print('Solution: x = %s   y = %f  z = %f' % (str(x(r)), y(r), z(r)))
# Solution: x = [-4.25 -4.25]   y = -20.000000  z = 4.000000

from FuncDesigner import *
 
# create some variables
a, b, c = oovars('a', 'b', 'c')
# or just a, b, c = oovars(3)
 
# Python list of 3 linear equations with 3 variables
f = [2*a+3*b-2*c+5, 2*a+13*b+15, a+4*b+2*c-45]
# alternatively, you could pass equations (the example below is valid for FD 0.18+):
# f = [2*a+3*b-2*c == -5, 2*a+15 == -13*b, a+45==4*b+2*c]
# assign SLE
linSys = sle(f)
 
r = linSys.solve()
A, B, C =  r(a, b, c)
maxRes = r.ff # max residual
# Expected result:
# {b: array([-5.]), a: array([ 25.]), c: array([ 20.])}

from FuncDesigner import *
n = 100
# create some variables
a, b, c = oovar(), oovar(size=n), oovar(size=2*n)
 
# let's construct some linear functions R^i -> R^j
# in Python range(m, k) is [m, m+1, ..., m+k-1]
f1 = a + sum(b*range(5, n+5)) # R^(n+1) -> R
f2 = a + 2*b + c.sum() # R^(2n+1) -> R^n
 
# you could use size of oovars
f3 = a + a.size + 2*c.size # R^(2n+1) -> R; a.size and c.size will be resolved into 1 and 2*n
 
f4 = c + range(4, 2*n+4) + f1 + 0.5*f2.sum() + 4*f3
 
# We can use "for" cycle:
for i in xrange(4):
    f4 = 0.5*f4 + a + f1 + 1
 
# Also, we could use matrix multiplication, eg f5 = dot(someMatrix, f4):
rng = 1.5 + cos(range(2*n)).reshape(-1, 1) # define 2n x 1 vector
R = dot(rng, rng.T) # create a matrix of shape 2n x 2n
f4 = dot(R, f4) # involve matrix multiplication
 
# Create Python list of linear equations 
f = [a+f4+5, 2*a+b*range(10, n+10)+15, a+4*b.sum()+2*c.sum()-45]
# alternatively, you could pass equations:
#f = [a+f4==-5, 2*a+b==-15, a==-4*b.sum()-2*c.sum()+45]
 
linSys = sle(f)
r = linSys.solve()
 
# get result
A, B, C =  r(a,b,c)
# A, B, C will be numpy arrays of sizes 1,  n,  2*n
# if failed to solve - A, B, C will contain numpy.nan(s)
 
print('A=%f B[0]=%f C[15]=%f' % (A, B[0], C[15]))
# for n=100 it prints A=-5.000000 B[0]=-0.500000 C[15]=-1883724.947909
 
# you could get residuals
# here it will be Python list of length 3 
# that is number of (vectorized) equations
# with elements of type numpy.array
maxResidual = r.ff
residuals = [F(r) for F in f]

from FuncDesigner import *
from time import time
t = time()
 
n = 1000
 
a, b, c = oovar(), oovar(size=n), oovar(size=2*n)
 
f1 = a + sum(b*range(5, n+5))
f2 = a + 2*b + c.sum() 
f3 = a + a.size + 2*c.size 
f4 = c + range(4, 2*n+4) + 4*f3
 
f = [a+f4+5, 2*a+b*range(10, n+10)+15, a+4*b.sum()+2*c.sum()-45]
# alternatively, you could pass equations:
#f = [a+f4==-5, 2*a+b==-15, a==-4*b.sum()-2*c.sum()+45]
linSys = sle(f)
 
r = linSys.solve() # i.e. using autoselect - solver numpy.linalg.solve for dense (for current numpy 1.4 it's LAPACK dgesv)
# and scipy.sparse.linalg.spsolve for sparse SLEs (for current scipy 0.8 it's LAPACK dgessv)
 
A, B, C =  r(a,b,c)
 
print('A=%f B[4]=%f B[first]=%f C[last]=%f' % (A, B[4], B[0], C[-1]))
 
maxResidual = r.ff
 
# Note - time may differ due to different matrices obtained from SLE rendering
# because Python 2.6 doesn't has ordered sets (they are present in Python 3.x)
# maybe I'll implement fixed rendering in future for 3.x, I don't want to deal 
# with quite difficult walkaround for 2.6 
print('time elapsed: %f' % (time()-t))

from FuncDesigner import *
 
# Create some variables
a, b = oovars('a', 'b')
c = oovar('c', size=2)
 
# Python list of linear equations
f = [2*a + 3*b - 2*(c+[1, 2]) * [3, 4] + 5, # R^nVars - > R^2
     2*a+13*b+15, # R^nVars - > R
     a+4*b+2*c.sum()-45# R^nVars - > R
     ]  
# Alternatively, you could pass equations:
# f = [2*a + 3*b - 2*(c+[1, 2]) * [3, 4]==-5, 
#            2*a+15==-13*b, a==-4*b-2*c.sum()+45]
 
# Assign SLE
linSys = sle(f)
 
# Rendered matrices for Ax=b are linSys.A (numpy.array of shape 4 x 4) and linSys.b (numpy.array of length 4)
A, B = linSys.A, linSys.b # B instead b for to prevent overriding the oovar b defined above
# However, for large sparse problems they can have other type than numpy.ndarray, 
# e.g. scipy.sparse lil_matrix or coo_matrix (provided you have scipy installed);
# elseware sometimes they can be mere out of your computer RAM.
# You could cast it to ndarray for example via myArr = myMatrix.toarray()
 
"""
Now let's solve it, for this example I use numpy.linalg.solve 
(that is default dense solver for FD SLEs)
but you can use any other linear solver 
e.g. cvxopt, scipy.lapack or scipy.sparse solvers, PyAMG etc.
Also you can code it (via scipy.sparse or somehow else) into 
one of common sparse matrices formats,
save it into a file and solve it via C, Fortran etc
(BTW f2py or Cython could be useful here)
Current drawback: for huge sparse SLEs you should have enough memory 
to handle dense A before casting to sparse will be performed.
Maybe I'll fix it in future, currently for my purposes I have no deal with these situations
"""
from numpy import linalg
xf = linalg.solve(A, B)
 
# Decode solution
sol = linSys.decode(xf)
print(c(sol)[0])
print(sol)
# Expected output:
# 10.2
# {b: array([-7.72]), c: array([ 10.2,   6.4]), a: array([ 42.68])}

from FuncDesigner import *
from openopt import NLP
 
a, b, c = oovars('a', 'b', 'c')
f = sum(a*[1, 2])**2+b**2+c**2
startPoint = {a:[100, 12], b:2, c:40} # however, you'd better use numpy arrays instead of Python lists
p = NLP(f, startPoint)
p.constraints = [(2*c+a-10)**2 < 1.5 + 0.1*b, (a-10)**2<1.5, a[0]>8.9, a+b > [ 7.97999836, 7.8552538 ], \
a < 9, (c-2)**2 < 1, b < -1.02, c > 1.01, (b + c * log10(a).sum() - 1) ** 2==0]
r = p.solve('ralg')
print(r.xf)
a_opt, b_opt, c_opt = a(r), b(r), c(r)
# For FuncDesigner v >= 0.17.dev you can use any of the following as well: 
# a_opt,  b_opt, c_opt = r(a, b, c)
# a_opt,  b_opt, c_opt = r(a), r(b), r(c)
# r('a'), r('b'), r('c') (provided you have assigned the names to oovars as above)
# r[a], r[b], r[c]
# r['a'], r['b'], r['c']
# a(r), b(r), c(r)
 
"""
Expected output:
...
objFunValue: 717.75631 (feasible, max constraint =  7.44605e-07)
{a: array([ 8.99999792,  8.87525277]), b: array([-1.01999971]), c: array([ 1.0613562])}
"""

# Example for NSP 
# |x[0]| + 1.2*|x[1]| + 1.44*|x[2]| + ... + 1.2^74*|x[74]| + |y-15| +|y+15| + y^2 +  -> min
# coded in FuncDesigner
# all 1st derivatives are obtained via Automatic Differentiation,
# not to be confused with finite-difference derivatives approximation
 
from FuncDesigner import *
from openopt import NSP
from numpy import cos, arange
 
x, y = oovars('x', 'y')
 
N = 75
koeffs = arange(1, N+1) ** 1.2
 
f = (abs(x) * koeffs).sum() + abs(y-15) + abs(y+15) + y**2
 
startPoint = {x: cos(1+arange(N)), y:80}
 
p = NSP(f, startPoint, maxIter = 1e5)
 
r = p.solve('ralg')
x_opt, y_opt = x(r), y(r)
print(max(abs(x_opt)), y_opt)
 
# Expected output:
# 1129  3.000e+01 
# istop:  4 (|| F[k] - F[k-1] || < ftol)
# Solver:   Time Elapsed = 9.27 	CPU Time Elapsed = 8.62
# objFunValue: 30.002138
# (9.6909792750048539e-06, array([  1.03949811e-05]))

# Solving system of equations:
# x**3 + y**3 - 9 = 0
# x - 0.5*y = 0
# cos(z) + x - 1.5 = 0
 
from FuncDesigner import *
from openopt import NLSP
 
x, y, z = oovars(3)
 
f = (x**3 + y**3 - 9, x - 0.5*y, cos(z) + x - 1.5)
 
startPoint = {x:8, y:15, z:80}
 
p = NLSP(f, startPoint)
 
# optional: we could set some constraints
p.constraints = [z<70,  z>50,  z + sin(x) < 60]
 
r = p.solve('nssolve') # nssolve is name of solver involved, see OOF doc for more arguments
print('Solution: x = %f   y = %f  z = %f' % (x(r), y(r), z(r)))
# Expected:
# Solution: x = 1.000001   y = 2.000000  z = 57.595865
 
"""                 More complex example                 """
# Solving non-linear system of size 152
from FuncDesigner import *
from openopt import NLSP
from numpy import arange, ones
N = 150
 
x = oovar(size=N)
y,z = oovars(2)
 
f = (x.sum()**2 + y - 100*N, 1e-3*sum(x**2) - 5*y - 100*N - (z-1)**2, (x-15)*arange(N)/(100*N) - 4*z + 5)
 
startPoint = {x: ones(N), y:15, z:80}
 
p = NLSP(f, startPoint, maxIter=1e7)
 
r = p.solve('scipy_fsolve') # nssolve is name of solver involved, see OOF doc for more arguments
print('x[5]:%f  y:%f  z:%f' % (x(r)[5],  y(r),  z(r)))
#x[5]:15.000000  y:-2107.773183  z:1.250000

# Example of creating MILP in FuncDesigner
# and solving it via OpenOpt
 
from FuncDesigner import *
from openopt import MILP
 
# Define some oovars
x, y, z = oovars(3)
 
# Let's define some linear functions
f1 = 4*x+5*y + 3*z + 5
f2 = f1.sum() + 2*x + 4*y + 15
f3 = 5*f1 + 4*f2 + 20
 
# Define objective; sum(a) and a.sum() are same as well as for numpy arrays
obj = x.sum() + y - 50*z + sum(f3) + 2*f2.sum() + 4064.6
 
# Start point - currently matters only size of variables
startPoint = {x:[8, 15], y:25, z:80} # however, using numpy.arrays is more recommended than Python lists
 
# Create prob
p = MILP(obj, startPoint, intVars = [y, z])
 
# Define some constraints
p.constraints = [x+5*y<15, x[0]<-5, f1<[25, 35], f1>-100, 2*f1+4*z<[80, 800], 5*f2+4*z<100, [-5.5, -4.5]<x,  x<1, -17<y,  y<20, -4000<z, z<4]
 
# Solve
r = p.solve('glpk') # glpk is name of the solver involved, see OOF doc for more arguments
 
# Decode solution
print('Solution: x = %s   y = %f  z = %f' % (str(x(r)), y(r), z(r)))
# Solution: x = [-5.  -4.5]   y = -17.000000  z = 0.000000

"""
GLP (GLobal Problem from OpenOpt set) example for FuncDesigner model:
searching for global minimum of the func 
(x-1.5)**2 + sin(0.8 * y ** 2 + 15)**4 + cos(0.8 * z ** 2 + 15)**4 + (t-7.5)**4
subjected to some constraints
See http://openopt.org/GLP for more info and examples.
"""
from openopt import GLP
from FuncDesigner import *
 
x, y, z, t = oovars(4)
 
# define objective
f = (x-1.5)**2 + sin(0.8 * y ** 2 + 15)**4 + cos(0.8 * z ** 2 + 15)**4 + (t-7.5)**4
 
# define some constraints
constraints = [x<1, x>-1, y<1, y>-1, z<1, z>-1, t<1, t>-1, x+2*y>-1.5,  sinh(x)+cosh(z)+sinh(t) <2.0]
 
# add some more constraints via Python "for" cycle
M = 10
for i in xrange(M):
    func = i*x+(M-i)*y+sinh(z)+cosh(t)
    constraints.append(func < i+1)
 
# define start point. You can use variables with length > 1 as well
startPoint = {x:0, y:0, z:0, t:0}
 
# assign prob
p = GLP(f, startPoint, constraints=constraints,  maxIter = 1e3,  maxFunEvals = 1e5,  maxTime = 5,  maxCPUTime = 5)
 
#optional: graphic output
#p.plot = 1 or p.solve(..., plot=1) or p = GLP(..., plot=1)
 
# solve
r = p.solve('de', plot=1) # try other solvers: galileo, pswarm
 
x_opt, y_opt, z_opt, t_opt = r(x, y, z, t)
# or 
# x_opt, y_opt, z_opt, t_opt = x(optPoint), y(optPoint), z(optPoint), t(optPoint)

from FuncDesigner import *
from openopt import NLP
 
a, b, c = oovars('a', 'b', 'c')
startPoint = {a:[100, 12], b:-0.9, c:40} # however, you'd better use numpy arrays instead of Python lists
 
objective = ((a*[1, 2])**2).sum() + b**2 + c**2
 
# Optional: set come constraints
constraints = [a>[1, 1], # it means a[0]>=1, a[1]>=1
               (a.sum()+b)**2<10, 
               (a+b+c)**2 < 15] # it means (a[i]+b+c)**2 < 15 for  i in {0,1}
 
fixedVars = b # or fixedVars = [b] / c / [b,c] / [a,c] etc
p = NLP(objective, startPoint, fixedVars = fixedVars, constraints = constraints)
 
r = p.solve('ralg')
print r.xf
 
# Alternatively, you can set optVars instead of fixedVars: 
optVars = [a, c] # or optVars = [a] / c / [b,c] / [a,c] etc
p = NLP(objective, startPoint, optVars = optVars, constraints = constraints)
r = p.solve('ralg')
print r.xf
 
"""
Expected output:
...
objFunValue: 5.809995 (feasible, max constraint =  7.62729e-07)
{b: -0.90000000000000002, a: array([ 1.00000032,  0.99999924]), c: array([ 0.00068568])}
"""