In [1]:
from pycalphad import Database, Model
from pycalphad.core.utils import NumPyPrinter
from pycalphad.core.autograd_utils import build_functions
import pycalphad.variables as v
from sympy import lambdify
import numpy as np
import algopy
from algopy import UTPM
import itertools
mod = Model(Database('2016-04-01-AlNi.tdb'), ['AL', 'NI', 'VA'], 'FCC_L12')
func = lambdify(tuple([v.T] + mod.site_fractions), mod.ast, dummify=True,
modules=[algopy])
print(tuple([v.T] + mod.site_fractions))
inp_arr = np.tile([[300,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,1], [600, 0.4, 0.6,0.4,0.6,0.4,0.6,0.4,0.6,1]], (400,1))
In [2]:
def extract_hessian(N, y):
H = np.zeros((y.data.shape[1], N,N), dtype=y.data.dtype)
for n in range(N):
for m in range(n):
a = sum(range(n+1))
b = sum(range(m+1))
k = sum(range(n+2)) - m - 1
#print 'k,a,b=', k,a,b
if n!=m:
tmp = (y.data[2, :, k] - y.data[2, :, a] - y.data[2, :, b])
H[:, m,n]= H[:, n,m]= tmp
a = sum(range(n+1))
H[:, n,n] = 2*y.data[2, :, a]
return H
In [3]:
%%timeit
# generate directions
N = len(mod.variables)
M = (N*(N+1))/2
L = (N*(N-1))/2
S = np.zeros((N,M))
s = 0
i = 0
for n in range(1,N+1):
S[-n:, s:s+n] = np.eye(n)
S[-n, s:s+n] = np.ones(n)
s += n
i += 1
#print(S)
S = S[::-1].T
print(S)
x = algopy.UTPM(np.zeros((3, inp_arr.shape[0]) + S.shape))
x.data[0, :, :, :] = inp_arr[..., None, :]
x.data[1, :, :] = S
y = func(*[x[..., i] for i in range(len(mod.variables))])
#obj_algopy = y.data[0, :, 0]
#grad_algopy = y.data[1, :, np.cumsum(np.arange(N), dtype=np.int)].T
#hess_algopy = extract_hessian(N, y)
#print(obj_algopy)
#print(grad_algopy)
#print(hess_algopy)
In [4]:
obj, grad, hess = build_functions(mod.ast, mod.variables)
In [5]:
%%timeit
o = obj(*inp_arr.T)
g = grad(*inp_arr.T)
h = hess(*inp_arr.T)
#print(o)
#print(g)
#print(h)
In [6]:
import numpy.testing
numpy.testing.assert_allclose(obj_algopy, o)
numpy.testing.assert_allclose(grad_algopy, g)
numpy.testing.assert_allclose(hess_algopy, h)
print('equivalent')
In [ ]:
In [ ]: