In [2]:

    

%load_ext Cython


    

    




/home/rotis/anaconda/envs/calphadpy3/lib/python3.5/site-packages/Cython/Distutils/old_build_ext.py:30: UserWarning: Cython.Distutils.old_build_ext does not properly handle dependencies and is deprecated.
  "Cython.Distutils.old_build_ext does not properly handle dependencies "

In [68]:

    

%%cython
from pycalphad.core.rksum import RedlichKisterSum
from pycalphad import Database
cimport numpy as np
import numpy as np
from tinydb import where
from sympy import Symbol

cdef np.ndarray[ndim=1, dtype=np.float64_t] _eval_rk_matrix_gradient(double[:,:] coef_mat, double[:,:] symbol_mat,
                                   double[:] eval_row, double[:] parameters):
    cdef np.ndarray[ndim=1, dtype=np.float64_t] out = np.zeros(eval_row.shape[0]-2)
    cdef double result = 0
    cdef double prod_result
    cdef int row_idx1 = 0
    cdef int row_idx2 = 0
    cdef int col_idx = 0
    cdef int dof_idx
    # eval_row order: P,T,ln(P),ln(T),y...
    # dof order: P,T,y...
    # coef_mat order: low_temp,high_temp,P,T,ln(P),ln(T),y...,constant_term,parameter_value
    for dof_idx in range(eval_row.shape[0]-2):
        if coef_mat.shape[1] > 0:
            for row_idx1 in range(coef_mat.shape[0]):
                if (eval_row[1] >= coef_mat[row_idx1, 0]) and (eval_row[1] < coef_mat[row_idx1, 1]):
                    if dof_idx < 2:
                        # special handling for state variables since they also can have a ln term
                        if coef_mat[row_idx1, 2+dof_idx] != 0:
                            prod_result = coef_mat[row_idx1, coef_mat.shape[1]-2] * coef_mat[row_idx1, coef_mat.shape[1]-1]
                            for col_idx in range(coef_mat.shape[1]-4):
                                if col_idx == dof_idx:
                                    prod_result *= (eval_row[col_idx]**(coef_mat[row_idx1, 2+col_idx]-1) * eval_row[2+col_idx]**(coef_mat[row_idx1, 4+col_idx]-1) * (coef_mat[row_idx1, 4+col_idx]+coef_mat[row_idx1, 2+col_idx]*(eval_row[2+col_idx])))
                                else:
                                    prod_result *= (eval_row[col_idx] ** coef_mat[row_idx1, 2+col_idx])
                            out[dof_idx] += prod_result
                    else:
                        if coef_mat[row_idx1, 4+dof_idx] != 0:
                            prod_result = coef_mat[row_idx1, coef_mat.shape[1]-2] * coef_mat[row_idx1, coef_mat.shape[1]-1]
                            for col_idx in range(coef_mat.shape[1]-4):
                                if col_idx == 2+dof_idx:
                                    prod_result *= (coef_mat[row_idx1, 4+dof_idx] * eval_row[col_idx] ** (coef_mat[row_idx1, 2+col_idx] - 1))
                                else:
                                    prod_result *= (eval_row[col_idx] ** coef_mat[row_idx1, 2+col_idx])
                            out[dof_idx] += prod_result
        if symbol_mat.shape[1] > 0:
            for row_idx2 in range(symbol_mat.shape[0]):
                if (eval_row[1] >= symbol_mat[row_idx2, 0]) and (eval_row[1] < symbol_mat[row_idx2, 1]):
                    if dof_idx < 2:
                        # special handling for state variables since they also can have a ln term
                        if coef_mat[row_idx2, 2+dof_idx] != 0:
                            prod_result = coef_mat[row_idx2, coef_mat.shape[1]-2] * parameters[<int>symbol_mat[row_idx2, symbol_mat.shape[1]-1]]
                            for col_idx in range(coef_mat.shape[1]-4):
                                if col_idx == dof_idx:
                                    prod_result *= (eval_row[col_idx]**(coef_mat[row_idx1, 2+col_idx]-1) * eval_row[2+col_idx]**(coef_mat[row_idx1, 4+col_idx]-1) * (coef_mat[row_idx1, 4+col_idx]+coef_mat[row_idx1, 2+col_idx]*(eval_row[2+col_idx])))
                                else:
                                    prod_result *= (eval_row[col_idx] ** coef_mat[row_idx1, 2+col_idx])
                            out[dof_idx] += prod_result
                    else:
                        if coef_mat[row_idx1, 4+dof_idx] != 0:
                            prod_result = coef_mat[row_idx2, coef_mat.shape[1]-2] * parameters[<int>symbol_mat[row_idx2, symbol_mat.shape[1]-1]]
                            for col_idx in range(coef_mat.shape[1]-4):
                                if col_idx == 2+dof_idx:
                                    prod_result *= (coef_mat[row_idx2, 4+dof_idx] * eval_row[col_idx] ** (coef_mat[row_idx2, 2+col_idx] - 1))
                                else:
                                    prod_result *= (eval_row[col_idx] ** coef_mat[row_idx2, 2+col_idx])
                            out[dof_idx] += prod_result
    return out

dbf = Database('alcocr-sandbox.tdb')
param_query = (
    (where('phase_name') == 'LIQUID') & \
    (where('parameter_type') == "G")
)
parameters = {}
all_symbols = dbf.symbols.copy()
# Convert string symbol names to sympy Symbol objects
# This makes xreplace work with the symbols dict
all_symbols = dict([(Symbol(s), val) for s, val in all_symbols.items()])
for param in parameters.keys():
    all_symbols.pop(param, None)
rks = RedlichKisterSum(['AL','CO','CR','VA'], dbf.phases['LIQUID'], dbf.search, param_query, list(parameters.keys()), all_symbols)
exx = _eval_rk_matrix_gradient(rks.output_matrix[2:4,:], np.array([[]]),
                        np.array([1013235, 300, np.log(101325), np.log(300), 0.3, 0.3, 0.4]), np.array([]))
print(exx)


    

    




[  0.00000000e+00  -3.55256010e+00  -1.15287101e+04   0.00000000e+00
   0.00000000e+00]

In [69]:

    

rks.output_matrix[2:4,:]


    

    
Out[69]:





array([[  2.98150000e+02,   9.33600000e+02,   0.00000000e+00,
          1.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          1.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          1.00000000e+00,  -1.18418670e+01],
       [  2.98150000e+02,   7.00000000e+02,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          1.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          1.00000000e+00,  -7.97615000e+03]])

In [70]:

    

-1.18418670e01 * 300 + -7.97615000e03


    

    
Out[70]:





-11528.7101

In [ ]: