import time
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
|
|
def my_log_exp_plus_one(z):
x = np.atleast_1d(z)
y = x*0
x_plus = x >= 0
x_minus= x < 0
y[x_minus] = np.log(1+np.exp(x[x_minus]))
y[x_plus] = x[x_plus] + np.log(1+np.exp(-x[x_plus]))
return y
|
|
def phi(x):
return np.abs(x)
#return np.exp(x)
#return np.sqrt(1+x*x)
#return x**2
#return my_log_exp_plus_one(x)
|
|
def diff_fun(xi,x):
return xi*x - phi(x)
|
|
def inf_recursive_search(xi,level,L):
if level > 8:
return (np.inf,0)
else:
N = 100
x = np.linspace(-L,L,N)
y = diff_fun(xi,x)
arg_max = np.argmax(y)
y_max = y[arg_max]
#print arg_max, y_max, level, L
#print y
if y_max > 0 and (arg_max == 0 or arg_max == N-1):
return inf_recursive_search(xi,level+1,2*L)
else:
return (y[arg_max],x[arg_max])
|
|
def fine_recursive_search(xi,center,level,L):
old_y_max = diff_fun(xi,center)
if level > 10:
return old_y_max
else:
tol = 1e-6
N = 100
x = np.linspace(center-L,center+L,N)
y = diff_fun(xi,x)
arg_max = np.argmax(y)
y_max = y[arg_max]
new_center = x[arg_max]
#print arg_max, y_max, level, center
#print 'Diff:',np.abs(old_y_max - y_max)
if y_max - old_y_max > tol:
return fine_recursive_search(xi,new_center,level+1,L/2)
else:
return max(y_max,old_y_max)
|
|
def phi_star(xi):
L = 10
value,center = inf_recursive_search(xi,0,L)
if np.isinf(value):
return value
return fine_recursive_search(xi,center,0,L)
|
|
N = 500
A = -3
B = 3
x_data = np.linspace(A,B,N)
y_data = x_data*0
start = time.time()
for k,x in enumerate(x_data):
y_data[k] = phi_star(x)
end = time.time()
print 'Computation time:',end-start,'seconds'
#print y_data
fig = plt.figure()
axes = fig.add_subplot(111)
#---------------------------------------------
min_val = y_data.min()
if np.isinf(min_val):
mx = 1
else:
inf_data = np.isinf(y_data)
finite_data = np.logical_not(inf_data)
mx = max(y_data[finite_data]) + 1
#---------------------------------------------
axes.plot(x_data[finite_data],y_data[finite_data],'bo',\
markeredgecolor='none')
axes.plot(x_data[inf_data],x_data[inf_data]*0+mx,'ro',\
markeredgecolor='none',markersize=8)
#---------------------------------------------Test_Exp
#x_data = np.linspace(0.1,3,20)
#y_data = x_data * np.log(x_data) - x_data
#axes.plot(x_data,y_data,'y+',linewidth=2)
#---------------------------------------------
x_data = np.linspace(A,B,N)
y_data = phi(x_data)
axes.plot(x_data,y_data,'g+',linewidth=2)
#---------------------------------------------
print 'Green curve ---> original function'
print 'Blue curve ---> convex conjugate'
print 'Red curve ---> +Infinity'
plt.axis('equal')
fig.savefig('test.png')
plt.close()
Computation time: 0.555828809738 seconds Green curve ---> original function Blue curve ---> convex conjugate Red curve ---> +Infinity Computation time: 0.555828809738 seconds Green curve ---> original function Blue curve ---> convex conjugate Red curve ---> +Infinity 
|
|
|