Scipy´s Multidimensional Optimization Progress without "callback-function" - an example of using Sympy, Numpy, Scipy and Pandas

Using Sympy , Numpy, Scipy and Pandas for 3D - Plotting and Optimization¶

In [20]:

from sympy import init_printing
init_printing() #for pretty latex Output
import sympy as sp #Symbolic Math Library
import scipy as sc #Scientific Python Library for Multidimensional Optimization
from scipy.optimize import minimize #Importing Minimize 
from sympy import * #Imports all special functions including cos, acos...
import numpy as np  #Imports numeric python
import re           #imports String manipulation and filtering functions
import pandas as pd #Importing for Dataframe - Excel Like Array - Printing

Plotting works only if there are 2 Independants (3D Plot)¶

In [3]:

x0, x1=sp.symbols('x0 x1') #define how many independent variables expression has!

In [4]:

F_string="Abs(((cos(x0))**4+(cos(x1))**4-2*((cos(x0))**2)*((cos(x1))**2))/(sqrt(x0**2+2*x1**2)))"
F_string

Out[4]:

'Abs(((cos(x0))**4+(cos(x1))**4-2*((cos(x0))**2)*((cos(x1))**2))/(sqrt(x0**2+2*x1**2)))'

In [5]:

F_sym=sp.sympify(F_string)

In [23]:

F_sym

Out[23]:

$$\left\lvert{\frac{1}{\sqrt{x_{0}^{2} + 2 x_{1}^{2}}} \left(\cos^{4}{\left (x_{0} \right )} - 2 \cos^{2}{\left (x_{0} \right )} \cos^{2}{\left (x_{1} \right )} + \cos^{4}{\left (x_{1} \right )}\right)}\right\rvert$$

This Part will work if 2 independants were defined previously¶

In [7]:

F_lam=sp.lambdify((x0,x1),F_sym,modules='numpy')#Matplotlib needs a lambda function
x0_t=np.linspace(0.5,10)# sets ranges for a Plot
x1_t=np.linspace(0.5,10)

X0,X1=np.meshgrid(x0_t,x1_t)
Z=F_lam(X0,X1)

In [8]:

from mpl_toolkits.mplot3d.axes3d import Axes3D

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X0 , X1, Z, cmap='autumn', cstride=2, rstride=2)

plt.draw()

Building an Optimization Function from the String Function¶

Works with more than 2 Independants¶

In [9]:

#F_opt_str=F_string.replace("x0","x[0]").replace("x1","x[1]")

#Finds indices and converts to a scipy.optimize expression 
indices = sorted(set(int(n) for n in re.findall(r'\bx(\d+)', F_string)),reverse=True)#reverse is important because otherwise x[10] becomes x[1]0
#Set takes unique objects (erases dublicates)
t = F_string[:]
#global space
space=len(indices)
for ind in indices:
    t = t.replace('x' + str(ind), 'x[' + str(ind) + ']')# if variable begins with "1" lower second to "ind-1"
F_opt_str=t

In [10]:

F_opt_str

Out[10]:

'Abs(((cos(x[0]))**4+(cos(x[1]))**4-2*((cos(x[0]))**2)*((cos(x[1]))**2))/(sqrt(x[0]**2+2*x[1]**2)))'

In [11]:

F_opt_lam=eval('lambda x:'+F_opt_str)#establish a lambda function of string_function

In [12]:

space

Out[12]:

$$2$$

In [14]:

#allows to See Iteration Progress: Method "COBYLA" does not allow callback functions 
class iter_progress:
    x_it=[]#allowing to access objects "outside functions"
    y_it=[]
    iter_log=[]
    def __init__(self):
        pass#iter_progress.x_iter
    def build_iter(self,space):
        import numpy as np
        iter_progress.y_it=np.zeros((1,1), dtype=float)
        iter_progress.x_it=np.zeros((1,space), dtype=float)
        iter_progress.iter_log=np.zeros((1,space+1), dtype=float)
    def obj(self,x):
        #x=np.array([1,2])
        out=F_opt_lam(x)
       # print(out)
        out=np.array([[out]])
        x_copy=x.copy()[None]
        #print("x_it=")
        iter_progress.x_it=np.concatenate((iter_progress.x_it, x_copy))
        #print(iter_progress.x_it)
        #print(" ")
        iter_progress.y_it=np.concatenate((iter_progress.y_it, out))
        #print(iter_progress.y_it)
        #print(" ")
        #print(out)
        return out
    def mine(self,x0=None):
        if x0==None:# if no starting value is defined "a zero array is defined" does not work always
            x0=np.zeros(((1,space)))
        print("x0=%s" % x0)
        #ip=iter_progress()
        #ip.build_iter
        res=minimize(ip.obj,x0=x0,method='SLSQP')
        iter_progress.iter_log=np.concatenate((iter_progress.x_it, iter_progress.y_it),axis=1)
        #print(iter_progress.x_it)
        print(" ")
        #print(iter_progress.y_it)
        print(iter_progress.iter_log)
        print(" ")
        print(res.x)

In [15]:

ip=iter_progress()
ip.build_iter(space)
xt = np.array([[4,6]])
ip.mine(xt)#Define an optional starting x0 if you want - x0 = np.array([[1,2,...length space !]])

x0=[[4 6]]
 
[[0.0 0.0 0.0]
 [4.0 6.0 0.0260856758941630]
 [4.0 6.0 0.0260856758941630]
 [4.000000014901161 6.0 0.0260856774313287]
 [4.0 6.000000014901161 0.0260856766844131]
 [3.896842558402568 5.9469672127161175 0.0140613728881540]
 [3.896842558402568 5.9469672127161175 0.0140613728881540]
 [3.8968425733037293 5.9469672127161175 0.0140613740373019]
 [3.896842558402568 5.946967227617279 0.0140613735821784]
 [3.548518480797617 5.7371142278847085 0.00144343287586224]
 [3.548518480797617 5.7371142278847085 0.00144343287586224]
 [3.548518495698778 5.7371142278847085 0.00144343259829832]
 [3.548518480797617 5.73711424278587 0.00144343253498852]
 [3.626154616515002 5.792628075807261 2.74906978159724e-6]
 [3.626154616515002 5.792628075807261 2.74906978159724e-6]
 [3.626154631416163 5.792628075807261 2.74905616988759e-6]
 [3.626154616515002 5.792628090708422 2.74905605476187e-6]
 [3.629920154620364 5.795149474914002 6.51408045966049e-9]
 [3.629920154620364 5.795149474914002 6.51408045966049e-9]
 [3.629920169521525 5.795149474914002 6.51474623195347e-9]
 [3.629920154620364 5.795149489815163 6.51474597461352e-9]]
 
[ 3.62992015  5.79514947]

Simple Optimization without Using the Class¶

In [16]:

#Simple Optimization without Using the Class and Progress
#x0=np.array([3,4])
#res=minimize(F_opt_lam,x0=x0,method='SLSQP')
#res.x
#res

"Pandas" DataFrame for a "pretty" table-view¶

In [17]:

var_str=[]# strings of Columns for pandas
indices_help=indices[::1]#reverse indices
for ind in indices_help:
    var_str.append('x'+str(ind))
var_str.append('f(x_i,n)')

In [18]:

var_str

Out[18]:

['x1', 'x0', 'f(x_i,n)']

In [21]:

df=pd.DataFrame(iter_progress.iter_log, columns=var_str) #Pretty Version of Optimization
df

Out[21]:

	x1	x0	f(x_i,n)
0	0	0	0
1	4	6	0.0260856758941630
2	4	6	0.0260856758941630
3	4	6	0.0260856774313287
4	4	6	0.0260856766844131
5	3.896843	5.946967	0.0140613728881540
6	3.896843	5.946967	0.0140613728881540
7	3.896843	5.946967	0.0140613740373019
8	3.896843	5.946967	0.0140613735821784
9	3.548518	5.737114	0.00144343287586224
10	3.548518	5.737114	0.00144343287586224
11	3.548518	5.737114	0.00144343259829832
12	3.548518	5.737114	0.00144343253498852
13	3.626155	5.792628	2.74906978159724e-6
14	3.626155	5.792628	2.74906978159724e-6
15	3.626155	5.792628	2.74905616988759e-6
16	3.626155	5.792628	2.74905605476187e-6
17	3.62992	5.795149	6.51408045966049e-9
18	3.62992	5.795149	6.51408045966049e-9
19	3.62992	5.795149	6.51474623195347e-9
20	3.62992	5.795149	6.51474597461352e-9

21 rows × 3 columns

In [22]:

from mpl_toolkits.mplot3d.axes3d import Axes3D

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X0 , X1, Z, cmap='autumn', cstride=2, rstride=2)
ax.plot(iter_progress.x_it[1:,0],iter_progress.x_it[1:,1], zs=iter_progress.y_it[1:,0], zdir='z', label='path')
ax.view_init(20,50)
plt.axis([3.1,3.9,4.5,5.8])
plt.draw()