1import numpy as np
2import time
3
4import benchmark1
5import use_adolc
6import use_scientific
7import use_uncertainties
8import use_algopy
9import use_numdifftools
10import use_funcdesigner
11import use_theano
12
13
14method = {'algopy_reverse_utps':0, 'algopy_reverse_utpm':1, 'pyadolc':2, 'scientific':3, 'uncertainties':4, 'numdifftools':5, 'algopy_forward_utpm':6, 'python':7, 'algopy_forward_utps':8, 'funcdesigner':9, 'theano':10}
15
16
17# FUNCTION COMPUTATION
18# --------------------
19
20function_N_list = [1,2,4,8,16,32,64,128,256]
21# function_N_list = [2]
22
23
24results_function_list = []
25for N in function_N_list:
26    print('N=',N)
27    results_function = np.zeros((11,3))
28
29    # pure python
30    f = benchmark1.F(N)
31    t = time.time(); pass ; preproc_time = time.time() - t
32    t = time.time();  ref_f =  f(3*np.ones(N));  run_time = time.time() - t
33    results_function[method['python']] = run_time, abs(ref_f - ref_f), preproc_time
34    print('ref_f=',ref_f)
35
36    # pyadolc
37    f = benchmark1.F(N)
38    t = time.time();  a = use_adolc.EVAL(f, np.ones(N), test='f'); preproc_time = time.time() - t
39    t = time.time();  f =  a.function(3*np.ones(N));  run_time = time.time() - t
40    results_function[method['pyadolc']] = run_time, abs(ref_f - f), preproc_time
41
42    # algopy utps version
43    f = benchmark1.F(N)
44    t = time.time();  a = use_algopy.EVAL(f, np.ones(N), test='f'); preproc_time = time.time() - t
45    t = time.time();  f =  a.function(3*np.ones(N));  run_time = time.time() - t
46    results_function[method['algopy_forward_utps']] = run_time, abs(ref_f - f), preproc_time
47
48    # algopy utpm version
49    f = benchmark1.G(N)
50    t = time.time();  a = use_algopy.EVAL2(f, np.ones(N), test='f'); preproc_time = time.time() - t
51    t = time.time();  f =  a.function(3*np.ones(N));  run_time = time.time() - t
52    results_function[method['algopy_forward_utpm']] = run_time, abs(ref_f - f), preproc_time
53
54    # funcdesigner
55    f = benchmark1.F(N)
56    t = time.time();  a = use_funcdesigner.EVAL(f, np.ones(N), test='f'); preproc_time = time.time() - t
57    t = time.time();  f =  a.function(3*np.ones(N));  run_time = time.time() - t
58    results_function[method['funcdesigner']] = run_time, abs(ref_f - f), preproc_time
59
60    # theano
61    f = benchmark1.F(N)
62    t = time.time();  a = use_theano.EVAL(f, np.ones(N), test='f'); preproc_time = time.time() - t
63    t = time.time();  f =  a.function(3*np.ones(N));  run_time = time.time() - t
64    results_function[method['theano']] = run_time, np.abs(ref_f - f), preproc_time
65    results_function_list.append(results_function)
66    # print f
67
68results_functions = np.array(results_function_list)
69print('results_functions=\n',results_functions)
70
71# GRADIENT COMPUTATION
72# --------------------
73
74gradient_N_list = [2,4,8,16,32,64,96]
75# gradient_N_list = [20]
76
77results_gradient_list = []
78for N in gradient_N_list:
79    print('N=',N)
80    results_gradient = np.zeros((11,3))
81
82    # pyadolc
83    f = benchmark1.F(N)
84    t = time.time();  a = use_adolc.EVAL(f, np.ones(N), test='g'); preproc_time = time.time() - t
85    t = time.time();  ref_g =  a.gradient(3*np.ones(N));  run_time = time.time() - t
86    results_gradient[method['pyadolc']] = run_time, 0, preproc_time
87
88    # scientifc
89    f = benchmark1.F(N)
90    t = time.time();  a = use_scientific.EVAL(f, np.ones(N), test='g'); preproc_time = time.time() - t
91    t = time.time();  g =  a.gradient(3*np.ones(N));  run_time = time.time() - t
92    results_gradient[method['scientific']] = run_time,  np.linalg.norm(g - ref_g)/np.linalg.norm(ref_g), preproc_time
93
94    # algopy, UTPS variant
95    f = benchmark1.F(N)
96    t = time.time();  a = use_algopy.EVAL(f, np.ones(N), test='g'); preproc_time = time.time() - t
97    t = time.time();  g =  a.gradient(3*np.ones(N));  run_time = time.time() - t
98    results_gradient[method['algopy_reverse_utps']] = run_time,  np.linalg.norm(g - ref_g)/np.linalg.norm(ref_g), preproc_time
99
100    # algopy, UTPM variant
101    f = benchmark1.G(N)
102    t = time.time();  a = use_algopy.EVAL2(f, np.ones(N), test='g'); preproc_time = time.time() - t
103    t = time.time();  g =  a.gradient(3*np.ones(N));  run_time = time.time() - t
104    results_gradient[method['algopy_reverse_utpm']] = run_time,  np.linalg.norm(g - ref_g)/np.linalg.norm(ref_g), preproc_time
105
106    # algopy, forward UTPM variant
107    f = benchmark1.G(N)
108    t = time.time();  a = use_algopy.EVAL2(f, np.ones(N), test='fg'); preproc_time = time.time() - t
109    t = time.time();  g = a.forwardgradient(3*np.ones(N));  run_time = time.time() - t
110    results_gradient[method['algopy_forward_utpm']] = run_time,  np.linalg.norm(g - ref_g)/np.linalg.norm(ref_g), preproc_time
111
112    # numdifftools
113    f = benchmark1.F(N)
114    t = time.time();  a = use_numdifftools.EVAL(f, np.ones(N), test='g'); preproc_time = time.time() - t
115    t = time.time();  g =  a.gradient(3*np.ones(N));  run_time = time.time() - t
116    results_gradient[method['numdifftools']] = run_time,  np.linalg.norm(g - ref_g)/np.linalg.norm(ref_g), preproc_time
117
118    # uncertainties
119    f = benchmark1.F(N)
120    t = time.time();  a = use_uncertainties.EVAL(f, np.ones(N), test='g'); preproc_time = time.time() - t
121    t = time.time();  g =  a.gradient(3*np.ones(N));  run_time = time.time() - t
122    results_gradient[method['uncertainties']] = run_time,  np.linalg.norm(g - ref_g)/np.linalg.norm(ref_g), preproc_time
123
124    # funcdesigner
125    f = benchmark1.F(N)
126    t = time.time();  a = use_funcdesigner.EVAL(f, np.ones(N), test='g'); preproc_time = time.time() - t
127    t = time.time();  g =  a.gradient(3*np.ones(N));  run_time = time.time() - t
128    results_gradient[method['funcdesigner']] = run_time,  np.linalg.norm(g - ref_g)/np.linalg.norm(ref_g), preproc_time
129
130    # theano
131    f = benchmark1.F(N)
132    t = time.time();  a = use_theano.EVAL(f, np.ones(N), test='g'); preproc_time = time.time() - t
133    t = time.time();  g =  a.gradient(3*np.ones(N));  run_time = time.time() - t
134    results_gradient[method['theano']] = run_time,  np.linalg.norm(g - ref_g)/np.linalg.norm(ref_g), preproc_time
135
136    results_gradient_list.append(results_gradient)
137
138results_gradients = np.array(results_gradient_list)
139print('results_gradients=\n',results_gradients)
140
141# HESSIAN COMPUTATION
142# -------------------
143print('starting hessian computation ')
144results_hessian_list = []
145hessian_N_list = [1,2,4,8,16,32,64]
146# hessian_N_list = [2]
147
148for N in hessian_N_list:
149    print('N=',N)
150    results_hessian = np.zeros((11,3))
151
152    # pyadolc
153    f = benchmark1.F(N)
154    t = time.time();  a = use_adolc.EVAL(f, np.ones(N), test='h'); preproc_time = time.time() - t
155    t = time.time();  ref_H =  a.hessian(3*np.ones(N));  run_time = time.time() - t
156    results_hessian[method['pyadolc']] = run_time, 0, preproc_time
157
158    # numdifftools
159    f = benchmark1.F(N)
160    t = time.time();  a = use_numdifftools.EVAL(f, np.ones(N), test='h'); preproc_time = time.time() - t
161    t = time.time();  H =  a.hessian(3*np.ones(N));  run_time = time.time() - t
162    results_hessian[method['numdifftools']] = run_time, np.linalg.norm( (H-ref_H).ravel())/ np.linalg.norm( (ref_H).ravel()), preproc_time
163
164    # algopy forward utpm variant
165    f = benchmark1.G(N)
166    t = time.time();  a = use_algopy.EVAL2(f, np.ones(N), test='fh'); preproc_time = time.time() - t
167    t = time.time();  H = a.forwardhessian(3*np.ones(N));  run_time = time.time() - t
168    results_hessian[method['algopy_forward_utpm']] = run_time, np.linalg.norm( (H-ref_H).ravel())/ np.linalg.norm( (ref_H).ravel()), preproc_time
169
170    # algopy forward/reverse utpm variant
171    f = benchmark1.G(N)
172    t = time.time();  a = use_algopy.EVAL2(f, np.ones(N), test='h'); preproc_time = time.time() - t
173    t = time.time();  H = a.hessian(3*np.ones(N));  run_time = time.time() - t
174    results_hessian[method['algopy_reverse_utpm']] = run_time, np.linalg.norm( (H-ref_H).ravel())/ np.linalg.norm( (ref_H).ravel()), preproc_time
175
176    # theano
177    f = benchmark1.F(N)
178    t = time.time();  a = use_theano.EVAL(f, np.ones(N), test='h'); preproc_time = time.time() - t
179    t = time.time();  H = a.hessian(3*np.ones(N));  run_time = time.time() - t
180    results_hessian[method['theano']] = run_time, np.linalg.norm( (H-ref_H).ravel())/ np.linalg.norm( (ref_H).ravel()), preproc_time
181
182    results_hessian_list.append(results_hessian)
183
184
185results_hessians = np.array(results_hessian_list)
186
187print(hessian_N_list)
188print('results_hessians=\n',results_hessians)
189
190
191# PLOT RESULTS
192
193print(results_gradients.shape)
194
195import matplotlib.pyplot as pyplot
196import prettyplotting
197
198# plot function run times
199pyplot.figure()
200pyplot.title('Function run times')
201pyplot.loglog(function_N_list, results_functions[:,method['pyadolc'],0], '-ko', markerfacecolor='None', label = 'pyadolc')
202pyplot.loglog(function_N_list, results_functions[:,method['theano'],0], '-.k1', markerfacecolor='None', label = 'theano')
203
204pyplot.loglog(function_N_list, results_functions[:,method['python'],0], '-ks', markerfacecolor='None', label = 'python')
205pyplot.loglog(function_N_list, results_functions[:,method['algopy_forward_utps'],0], '-.k<', markerfacecolor='None', label = 'algopy scalar')
206pyplot.loglog(function_N_list, results_functions[:,method['algopy_forward_utpm'],0], '-.k>', markerfacecolor='None', label = 'algopy matrix')
207pyplot.loglog(function_N_list, results_functions[:,method['funcdesigner'],0], '-.ks', markerfacecolor='None', label = 'funcdesigner')
208pyplot.ylabel('time $t$ [seconds]')
209pyplot.xlabel('problem size $N$')
210pyplot.grid()
211leg = pyplot.legend(loc=2)
212frame= leg.get_frame()
213frame.set_alpha(0.4)
214pyplot.savefig('function_runtimes.png')
215
216# plot gradient run times
217pyplot.figure()
218pyplot.title('Gradient run times')
219pyplot.loglog(gradient_N_list, results_gradients[:,method['pyadolc'],0], '-ko', markerfacecolor='None', label = 'pyadolc')
220pyplot.loglog(gradient_N_list, results_gradients[:,method['algopy_reverse_utps'],0], '-.k<', markerfacecolor='None', label = 'algopy reverse utps')
221pyplot.loglog(gradient_N_list, results_gradients[:,method['algopy_reverse_utpm'],0], '-.k>', markerfacecolor='None', label = 'algopy reverse utpm')
222pyplot.loglog(gradient_N_list, results_gradients[:,method['algopy_forward_utpm'],0], '-.kv', markerfacecolor='None', label = 'algopy forward utpm')
223pyplot.loglog(gradient_N_list, results_gradients[:,method['uncertainties'],0], '--kh', markerfacecolor='None', label = 'uncertainties')
224pyplot.loglog(gradient_N_list, results_gradients[:,method['numdifftools'],0], '--ks', markerfacecolor='None', label = 'numdifftools')
225pyplot.loglog(gradient_N_list, results_gradients[:,method['funcdesigner'],0], '--kd', markerfacecolor='None', label = 'funcdesigner')
226pyplot.loglog(gradient_N_list, results_gradients[:,method['theano'],0], '-.k1', markerfacecolor='None', label = 'theano')
227pyplot.ylabel('time $t$ [seconds]')
228pyplot.xlabel('problem size $N$')
229pyplot.grid()
230leg = pyplot.legend(loc=2)
231frame= leg.get_frame()
232frame.set_alpha(0.4)
233pyplot.savefig('gradient_runtimes.png')
234
235# plot hessian run times
236pyplot.figure()
237pyplot.title('Hessian run times')
238pyplot.loglog(hessian_N_list, results_hessians[:,method['pyadolc'],0], '-ko', markerfacecolor='None', label = 'pyadolc (fo/rev)')
239pyplot.loglog(hessian_N_list, results_hessians[:,method['algopy_forward_utpm'],0], '-.k>', markerfacecolor='None', label = 'algopy (fo)')
240pyplot.loglog(hessian_N_list, results_hessians[:,method['algopy_reverse_utpm'],0], '-.k<', markerfacecolor='None', label = 'algopy (fo/rev)')
241pyplot.loglog(hessian_N_list, results_hessians[:,method['numdifftools'],0], '--ks', markerfacecolor='None', label = 'numdifftools')
242pyplot.loglog(hessian_N_list, results_hessians[:,method['theano'],0], '-.k1', markerfacecolor='None', label = 'theano')
243pyplot.ylabel('time $t$ [seconds]')
244pyplot.xlabel('problem size $N$')
245pyplot.grid()
246leg = pyplot.legend(loc=2)
247frame= leg.get_frame()
248frame.set_alpha(0.4)
249pyplot.savefig('hessian_runtimes.png')
250
251
252# plot gradient preprocessing times
253pyplot.figure()
254pyplot.title('Gradient preprocessing times')
255pyplot.loglog(gradient_N_list, results_gradients[:,method['pyadolc'],2], '-ko', markerfacecolor='None', label = 'pyadolc')
256pyplot.loglog(gradient_N_list, results_gradients[:,method['algopy_reverse_utps'],2], '-.k<', markerfacecolor='None', label = 'algopy reverse utps')
257pyplot.loglog(gradient_N_list, results_gradients[:,method['algopy_reverse_utpm'],2], '-.k>', markerfacecolor='None', label = 'algopy reverse utpm')
258pyplot.loglog(gradient_N_list, results_gradients[:,method['algopy_forward_utpm'],2], '-.kv', markerfacecolor='None', label = 'algopy forward utpm')
259pyplot.loglog(gradient_N_list, results_gradients[:,method['uncertainties'],2], '--kh', markerfacecolor='None', label = 'uncertainties')
260pyplot.loglog(gradient_N_list, results_gradients[:,method['numdifftools'],2], '--ks', markerfacecolor='None', label = 'numdifftools')
261pyplot.loglog(gradient_N_list, results_gradients[:,method['funcdesigner'],2], '--kd', markerfacecolor='None', label = 'funcdesigner')
262pyplot.loglog(gradient_N_list, results_gradients[:,method['theano'],2], '-.k1', markerfacecolor='None', label = 'theano')
263pyplot.ylabel('time $t$ [seconds]')
264pyplot.xlabel('problem size $N$')
265pyplot.grid()
266leg = pyplot.legend(loc=2)
267frame= leg.get_frame()
268frame.set_alpha(0.4)
269pyplot.savefig('gradient_preprocessingtimes.png')
270
271# plot hessian preprocessing times
272pyplot.figure()
273pyplot.title('Hessian preprocessing times')
274pyplot.loglog(hessian_N_list, results_hessians[:,method['pyadolc'],2], '-ko', markerfacecolor='None', label = 'pyadolc (fo/rev)')
275pyplot.loglog(hessian_N_list, results_hessians[:,method['algopy_forward_utpm'],2], '-.k>', markerfacecolor='None', label = 'algopy (fo)')
276pyplot.loglog(hessian_N_list, results_hessians[:,method['algopy_reverse_utpm'],2], '-.k<', markerfacecolor='None', label = 'algopy (fo/rev)')
277pyplot.loglog(hessian_N_list, results_hessians[:,method['numdifftools'],2], '--ks', markerfacecolor='None', label = 'numdifftools')
278pyplot.loglog(hessian_N_list, results_hessians[:,method['theano'],2], '-.k1', markerfacecolor='None', label = 'theano')
279pyplot.ylabel('time $t$ [seconds]')
280pyplot.xlabel('problem size $N$')
281pyplot.grid()
282leg = pyplot.legend(loc=2)
283frame= leg.get_frame()
284frame.set_alpha(0.4)
285pyplot.savefig('hessian_preprocessingtimes.png')
286
287# plot gradient errors
288pyplot.figure()
289pyplot.title('Gradient Correctness')
290pyplot.loglog(gradient_N_list, results_gradients[:,method['numdifftools'],1], '--ks', markerfacecolor='None', label = 'numdifftools')
291pyplot.loglog(gradient_N_list, results_gradients[:,method['pyadolc'],1], '-ko', markerfacecolor='None', label = 'pyadolc')
292pyplot.loglog(gradient_N_list, results_gradients[:,method['algopy_reverse_utps'],1], '-.k<', markerfacecolor='None', label = 'algopy reverse utps')
293pyplot.loglog(gradient_N_list, results_gradients[:,method['algopy_reverse_utpm'],1], '-.k>', markerfacecolor='None', label = 'algopy reverse utpm')
294pyplot.loglog(gradient_N_list, results_gradients[:,method['algopy_forward_utpm'],1], '-.kv', markerfacecolor='None', label = 'algopy forward utpm')
295pyplot.loglog(gradient_N_list, results_gradients[:,method['uncertainties'],1], '--kh', markerfacecolor='None', label = 'uncertainties')
296pyplot.loglog(gradient_N_list, results_gradients[:,method['funcdesigner'],1], '--kd', markerfacecolor='None', label = 'funcdesigner')
297pyplot.loglog(gradient_N_list, results_gradients[:,method['theano'],1], '-.k1', markerfacecolor='None', label = 'theano')
298pyplot.ylabel(r'relative error $\|g_{ref} - g\|/\|g_{ref}\}$')
299pyplot.xlabel('problem size $N$')
300pyplot.grid()
301leg = pyplot.legend(loc=2)
302frame= leg.get_frame()
303frame.set_alpha(0.4)
304pyplot.savefig('gradient_errors.png')
305
306# plot hessian errors
307pyplot.figure()
308pyplot.title('Hessian Correctness')
309pyplot.loglog(hessian_N_list, results_hessians[:,method['numdifftools'],1], '--ks', markerfacecolor='None', label = 'numdifftools')
310pyplot.loglog(hessian_N_list, results_hessians[:,method['pyadolc'],1], '-ko', markerfacecolor='None', label = 'pyadolc (fo/rev)')
311pyplot.loglog(hessian_N_list, results_hessians[:,method['algopy_forward_utpm'],1], '-.k>', markerfacecolor='None', label = 'algopy (fo)')
312pyplot.loglog(hessian_N_list, results_hessians[:,method['algopy_reverse_utpm'],1], '-.k<', markerfacecolor='None', label = 'algopy (fo/rev)')
313pyplot.loglog(hessian_N_list, results_hessians[:,method['theano'],1], '-.k1', markerfacecolor='None', label = 'theano')
314pyplot.ylabel(r'relative error $\|H_{ref} - H\|/\|H_{ref}\|$')
315pyplot.xlabel('problem size $N$')
316pyplot.grid()
317leg = pyplot.legend(loc=2)
318frame= leg.get_frame()
319frame.set_alpha(0.4)
320pyplot.savefig('hessian_preprocessingtimes.png')
321
322
323
324# pyplot.show()
325
326
327