cvex conjugate

1960 days ago by JAVIER

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