diffusivity_fitter.py

import numpy as np

import itertools
import scipy.optimize as optimize
from sklearn.metrics import r2_score

import matplotlib.pyplot as plt
import matplotlib.cm as cm

from decimal import *

getcontext().prec = 256
plt.close('all')


def poly_matrix(x, y, order):
    """ generate Matrix use with lstsq """
    ncols = (order + 1)**2
    G = np.zeros((x.size, ncols))
    ij = itertools.product(range(order+1), range(order+1))
    for k, (i, j) in enumerate(ij):
        G[:, k] = x**i * y**j
    return G

def expfunc(inputs, a, b, c): 
    
     x, y = np.log(inputs[0][:]), np.log(inputs[1][:])
#    x, y = (inputs[0][:]), (inputs[1][:])
    
     return a*np.exp(b*x)*np.exp(c*y)

def linfitter(x,y,z, ordr):
    
    # make Matrix:
    G = poly_matrix(x, y, ordr)
    # Solve for np.dot(G, m) = z:
    m = np.linalg.lstsq(G,z,rcond=None)[0]
    
    out = np.reshape(np.dot(G,m),x.shape)

    return G, m, out

def linpred(xx, yy, G, m, ordr): 
    
    GG = poly_matrix(xx.ravel(), yy.ravel(), ordr)
    return np.reshape(np.dot(GG, m), xx.shape)

def newfunc(inputs, a, b, c):
    
    x, y = (inputs[0][:]), (inputs[1][:])
    return a*x**b*y**c


###############################################################################
#
###############################################################################
#
#xdat = np.array([18,18, 130.077,130.77,130.077,130.77, 191.976,191.976,191.976,
#               230.093,230.093,230.093, 304.212, #434.5,434.5,
#               543.525, 558.64,558.64,558.64,558.64, 747.953,747.953,747.953,
#               786.98, 822.94,822.94, 912.085,
#               1435.6 #,42700
#               ])
#ydat = np.array([15,52,24,38,72,104,10,18,54, 34,48,96, 90, #18,54,
#        10, 8,32,37,45,12,46,120,#17,54,190,
#        17,5,10,17, 40 # ,22
#        ])
#zdat = np.array([8.8e-10,4.33e-10, 5e-16,3e-16,3.5e-16,2.25e-16, 5e-17,
#                 3.8e-17,7e-18, 1.1e-17,7e-18,2.5e-18, 5e-18, # 1e-17,3e-18,
#                 3e-18, 5e-19,7.5e-20,1e-20,2.5e-21, 5e-20,1.5e-20,5e-20, 1e-20,
#                 4.76e-19,4.44e-19, 5e-21,1e-21 # ,1.64e-21
#        ])
    
    
xdat = np.array([130.77,#467,
                 787,434.5,434.5,
                 225.21,1435.63,#302.236,
                 747.953,747.953,
                 747.953,359.35,
                 558.64,558.64,558.64,
                 558.64,130,130,304,543,
                 230,230,230,18,
                 18])
ydat = np.array([104,#12,
                 12,18,54,
                 30,40,#60,
                 12,46,
                 120,120,
                 8,32,37,
                 45,72,24,105,10,
                 34,48,96,15,
                 52])  
zdat = np.array([95e-18,#1e-21,
                 0.75e-20,1.75e-17,4e-18,
                 1e-19,1e-21,#3e-21,
                 2e-20,10e-21,
                 1.5e-21,5e-21,
                 2e-19,0.45e-19,0.1e-19,
                 2.5e-21,4.5e-16,0.5e-16,3e-18,6e-20,
                 0.6e-17,0.5e-17,0.8e-17,8.8e-10,
                 4.33e-10])     
                 
                 
               
                
                 



zdat_log=np.log(zdat)
#zdat_log=np.log10(zdat)

###############################################################################
#
###############################################################################

A = np.asarray([list(a) for a in zip(xdat,ydat,zdat_log)])

xi = np.linspace(min(xdat)*0.25,max(xdat)*1.25,len(xdat))
yi = np.linspace(min(ydat)*0.25,max(ydat)*1.1,len(ydat))
xx, yy  = np.meshgrid(xi, yi)
 
#------------------------------------------------------------------------------

zs = []
rsqs = []
for idx, ordr in enumerate([1,2,3,4,5]): 
    G, m, out = linfitter(xdat,ydat,zdat_log,ordr)
    zs.append(linpred(xx,yy,G,m,ordr))
    rsqs.append(r2_score(zdat_log, out))
    print('R\u00b2 (linear, order={0}): {1:1.4f}'.format(ordr, rsqs[idx]))
    print('Linear Params (order={0}): '.format(ordr), m) 

expp, pcov = optimize.curve_fit(expfunc, (xdat, ydat), zdat_log,ftol=1e-15, xtol=1e-15, maxfev=10000)
zexp = expfunc((xdat,ydat), *expp)
rsq = r2_score(zdat_log, zexp)
zz_exp = expfunc((xx,yy),*expp)
print('\na*exp(b*log(x))*exp(c*log(y))')
print('R\u00b2 (exponential): {0:1.4f}'.format(rsq))
print('Exponential Params: A={0}, B={1}, C={2}'.format(*expp))



print('\nz = exp( a(x^b)*(y^c) )')
newp, pcov_new = optimize.curve_fit(newfunc, (xdat, ydat), zdat, ftol=1e-15, xtol=1e-15, maxfev=800000)

znew = newfunc((xdat,ydat), *newp)
rsq_new = r2_score(zdat, znew)
print('R\u00b2 (new): {0:1.4f}'.format(rsq_new))
znew_list = np.asarray([list(a) for a in zip(zdat, znew)])

rsq_new_log = r2_score(np.log(zdat), np.log(znew))
print('R\u00b2 (new): {0:1.4f}'.format(rsq_new_log))

print('\nz = exp( a(x^b)*(y^c) )')
print('New Params: a={0}, b={1}, c={2}'.format(*newp))

zz = newfunc((xx,yy),*newp)

#------------------------------------------------------------------------------




#
#cmap = cm.RdBu
#
#fig, axes = plt.subplots(nrows=1, ncols=6,subplot_kw={'projection': '3d'})
#for idx, ax in enumerate(axes.flat): 
#    if idx < len(zs):
#        
#        im = ax.plot_surface(xx, yy, zs[idx], color=cmap(rsqs[idx]), alpha=0.2, vmin=0, vmax=1)
#        ax.set_title('R\u00b2 (lin, order={0}):{1:8.4f}'.format(idx+1, rsqs[idx]))
#    else: 
#        im = ax.plot_surface(xx,yy,zz,color=cmap(rsq), alpha=.2, vmin=0, vmax=1)
#        ax.set_title('R\u00b2 (exp):{1:8.4f}'.format(ordr, rsq))
#    
#    ax.set_zlim(-20, 0)
#    ax.scatter(xdat, ydat, zdat_log, c='k')
#    ax.set_zlim(-20, 0)
#
#
#m = cm.ScalarMappable(cmap=cm.RdBu)
#m.set_array(np.linspace(0,1,1000))
#plt.colorbar(m,alpha=1)
#
#
#fig.tight_layout()
#plt.show()

m = cm.ScalarMappable(cmap=cm.Blues)
m.set_array(np.linspace(0,1,1000))

fig = plt.figure(1,figsize=(9,6))
ax = plt.axes(projection='3d')
ax.set_xlabel('Drug Mw (Da)',fontsize=14)
ax.set_ylabel('Polymer Mw (kDa)',fontsize=14)
ax.set_zlabel('Diffusivity ($log_{e}(m^2/s)$)',fontsize=14)
title = '$R^2={0:1.4f}$ $D=exp(ax^by^c)$\na={1:1.4f}, b={2:1.4f}, c={3:1.4f}'.format(rsq_new_log, newp[0], newp[1], newp[2])
ax.set_title(title,fontsize=20)
ax.plot_surface(xx,yy,np.log(zz),cmap='Blues', edgecolor='none',alpha=0.5)
ax.scatter(xdat,ydat,zdat_log,color='k',s=30 )
ax.view_init(elev=25,azim=75)
plt.colorbar(m,alpha=1)
fig.tight_layout()



fig = plt.figure(2,figsize=(9,6))
ax = plt.axes(projection='3d')
ax.set_xlabel('Drug Mw (Da)',fontsize=14)
ax.set_ylabel('Polymer Mw (kDa)',fontsize=14)
ax.set_zlabel('Diffusivity ($log_{e}(m^2/s)$)',fontsize=14)
title = '$R^2={0:1.4f}'.format(rsq)+'$ log(D)=$ae^{bMw_{Poly}}10^{cMw_{Drug}}$'+'\na={1:1.4f}, b={2:1.4f}, c={3:1.4f}'.format(rsq, expp[0], expp[1], expp[2])

ax.set_title(title,fontsize=20)
ax.plot_surface(xx, yy, zz_exp,cmap='Blues', edgecolor='none',alpha=0.5)
ax.scatter(xdat,ydat,zdat_log,color='k',s=30 )
plt.colorbar(m,alpha=1)
ax.view_init(elev=25,azim=75)
plt.show()




def Diff(Mw_drug,Mw_poly):
    return np.exp(expp[0]*    Mw_drug**expp[1]  * Mw_poly**expp[2])
def Diff2(Mw_drug,Mw_poly):
    return newp[0]*(Mw_drug**newp[1])*(Mw_poly**newp[2])

print('MW \t \t Guess \t \t Actual')
print('6.5 kDa \t %1.2E \t %1.2E' %(Diff(331,6.5),40e-18))
print('21 kDa \t \t %1.2E \t %1.2E' %(Diff(331,21),4e-18))
print('34 kDa \t \t %1.2E \t %1.2E' %(Diff(331,34),3e-18))
print('49 kDa \t \t %1.2E \t %1.2E' %(Diff(331,49),2.5e-18))
print('67 kDa \t \t %1.2E \t %1.2E' %(Diff(331,67),1.5e-18))

print('MW \t \t Guess \t \t Actual')
print('6.5 kDa \t %1.2E \t %1.2E' %(Diff2(331,6.5),40e-18))
print('21 kDa \t \t %1.2E \t %1.2E' %(Diff2(331,21),4e-18))
print('34 kDa \t \t %1.2E \t %1.2E' %(Diff2(331,34),3e-18))
print('49 kDa \t \t %1.2E \t %1.2E' %(Diff2(331,49),2.5e-18))
print('67 kDa \t \t %1.2E \t %1.2E' %(Diff2(331,67),1.5e-18))