Re-write the script:

1. replace curve_fit with leastsq, which supports the analytical Jacobian
2. specify a "class" (contains both residum and jacobian) for each criterion.
3. add the calculation of Jacobian

processing/misc/yieldSurface.py

@@ -7,7 +7,7 @@ from scipy.optimize import curve_fit
 from scipy.linalg import svd
 from optparse import OptionParser
 import damask
-from damask.util import curve_fit_bound
+from damask.util import leastsqBound
 
 scriptID   = string.replace('$Id$','\n','\\n')
 scriptName = scriptID.split()[1][:-3]
@@ -77,52 +77,55 @@ def get_weight(ndim):
 # isotropic yield surfaces
 # ---------------------------------------------------------------------------------------------
 
-def Tresca(sigmas, sigma0):
+class Tresca(object):
   '''
     residuum of Tresca yield criterion (eq. 2.26)
   '''
-  lambdas = principalStresses(sigmas)
+  def fun(self,sigma0, ydata, sigmas):
-  r = np.amax(np.array([abs(lambdas[2,:]-lambdas[1,:]),\
+    lambdas = principalStresses(sigmas)
-                        abs(lambdas[1,:]-lambdas[0,:]),\
+    r = np.amax(np.array([abs(lambdas[2,:]-lambdas[1,:]),\
-                        abs(lambdas[0,:]-lambdas[2,:])]),0) - sigma0
+                          abs(lambdas[1,:]-lambdas[0,:]),\
-  return r.ravel()
+                          abs(lambdas[0,:]-lambdas[2,:])]),0) - sigma0
+    return r.ravel()
+  def jac(self,sigma0, ydata, sigmas):
+    return np.ones(len(ydata)) * (-1.0)
 
-
+class vonMises(object):
-def vonMises(sigmas, sigma0):
   '''
     residuum of Huber-Mises-Hencky yield criterion (eq. 2.37)
   '''
+  def fun(self, sigma0, ydata, sigmas):
+    return HosfordBasis(sigma0, 1.0,1.0,1.0, 2.0, sigmas)
+  def jac(self, sigma0, ydata, sigmas):
+    return HosfordBasis(sigma0, 1.0,1.0,1.0, 2.0, sigmas, Jac=True, nParas=1)
 
-  return Hosford(sigmas, sigma0, 2.0)
+class Drucker(object):
-
-
-def Drucker(sigmas, sigma0, C_D):
   '''
     residuum of Drucker yield criterion (eq. 2.41, F = sigma0)
   '''
+  def fun(sigma0, C_D, ydata, sigmas):
+    return DruckerBasis(sigma0, C_D, 1.0, ydata, sigmas)
+  def jac(sigma0, C_D, ydata, sigmas):
+    pass
 
-  return generalDrucker(sigmas, sigma0, C_D, 1.0)
+class generalDrucker(object):
-
-
-def generalDrucker(sigmas, sigma0, C_D, p):
   '''
     residuum of general Drucker yield criterion (eq. 2.42, F = sigma0)
   '''
-  Is = stressInvariants(principalStresses(sigmas))
+  def fun(sigma0, C_D, ydata, sigmas):
-  r  = (Is[1,:]**(3.0*p)-C_D*Is[2,:]**(2.0*p)) - sigma0
+    return DruckerBasis(sigma0, C_D, p, ydata, sigmas)
-  return r.ravel()
+  def jac(sigma0, C_D, ydata, sigmas):
+    pass
 
-
+class Hosford(object):
-def Hosford(sigmas, sigma0, a):
   '''
     residuum of Hershey yield criterion (eq. 2.43, Y = sigma0)
   '''
-  lambdas = principalStresses(sigmas)
+  def fun(self, (sigma0, a), ydata, sigmas):
-  r   = ((abs(lambdas[2,:]-lambdas[1,:]))**a\
+    return HosfordBasis(sigma0, 1.0,1.0,1.0, a, sigmas)
-      +  (abs(lambdas[1,:]-lambdas[0,:]))**a\
+  def jac(self, (sigma0, a), ydata, sigmas):
-      +  (abs(lambdas[0,:]-lambdas[2,:]))**a) **(1.0/a)\
+    return HosfordBasis(sigma0, 1.0,1.0,1.0, a, sigmas, Jac=True, nParas=2)
-      -2.0**(1.0/a)*sigma0
+
-  return r.ravel()
 
 #more to do
 # KarafillisAndBoyce
@@ -131,89 +134,57 @@ def Hosford(sigmas, sigma0, a):
 # isotropic yield surfaces
 # ---------------------------------------------------------------------------------------------
 
-def Hill1948(sigmas, F,G,H,L,M,N):
+class Hill1948(object):
   '''
     residuum of Hill 1948 quadratic yield criterion (eq. 2.48)
   '''
-  r =     F*(sigmas[1]-sigmas[2])**2.0\
+  def fun(self, (F,G,H,L,M,N), ydata, sigmas):
-    +     G*(sigmas[2]-sigmas[0])**2.0\
+    r = F*(sigmas[1]-sigmas[2])**2.0 + G*(sigmas[2]-sigmas[0])**2.0 + H*(sigmas[0]-sigmas[1])**2.0\
-    +     H*(sigmas[0]-sigmas[1])**2.0\
+      + 2.0*L*sigmas[4]**2.0 + 2.0*M*sigmas[5]**2.0 + 2.0*N*sigmas[3]**2.0 - 1.0
-    + 2.0*L* sigmas[4]**2.0\
+    return r.ravel()/2.0
-    + 2.0*M* sigmas[5]**2.0\
+  def jac(self, (F,G,H,L,M,N), ydata, sigmas):
-    + 2.0*N* sigmas[3]**2.0\
+    pass
-    - 1.0
-  return r.ravel()/2.0
 
 #more to do
 # Hill 1979
 
 # Hill 1990,1993 need special stresses to fit
 
-def generalHosford(sigmas, sigma0, a):
+class generalHosford(object):
   '''
-    residuum of Hershey yield criterion (eq. 2.104, sigma = sigma0)
+    residuum of Hershey yield criterion (eq. 2.104, sigmas = sigma0)
   '''
-  lambdas = principalStresses(sigmas)
+  def fun(self, (sigma0, F, G, H, a), ydata, sigmas, nParas=5):
-  r = np.amax(np.array([F*(abs(lambdas[:,1]-lambdas[:,2]))**a,\
+    return HosfordBasis(sigma0, F, G, H, a, sigmas)
-                        G*(abs(lambdas[:,2]-lambdas[:,0]))**a,\
+  def jac(self, (sigma0, F, G, H, a), ydata, sigmas):
-                        H*(abs(lambdas[:,0]-lambdas[:,1]))**a]),1) - sigma0**a
+    return HosfordBasis(sigma0, F,G,H, a, sigmas, Jac=True, nParas=5)
-  return r.ravel()
 
-
+class Barlat1991iso(object):
-def Barlat1991(sigmas, sigma0, order, a, b, c, f, g, h):
-  '''
-    residuum of Barlat 1997 yield criterion
-  '''
-  cos = np.cos; pi = np.pi; abs = np.abs
-  A = a*(sigmas[1] - sigmas[2])
-  B = b*(sigmas[2] - sigmas[0])
-  C = c*(sigmas[0] - sigmas[1])
-  F = f*sigmas[4]
-  G = g*sigmas[5]
-  H = h*sigmas[3]
-
-  I2 =   (F*F + G*G  +  H*H)/3.0  + ((A-C)**2+(C-B)**2+(B-A)**2)/54.0
-  I3 =   (C-B)*(A-C) * (B-A)/54.0 + F*G*H - \
-       ( (C-B)*F*F   + (A-C)*G*G  + (B-A)*H*H )/6.0
-  theta = np.arccos(I3/I2**1.5)
-  Phi = np.sqrt(3.0*I2)* (
-        (abs(2.0*cos((2.0*theta + pi)/6.0)))**order +
-        (abs(2.0*cos((2.0*theta + pi*3.0)/6.0)))**order +
-        (abs(2.0*cos(( 2.0*theta + pi*5.0)/6.0)))**order
-                         )**(1.0/order)
-  r   = Phi/2.0**(1.0/order) - sigma0
-
-
-  return r.ravel()
-
-def Barlat1991iso(sigmas, sigma0, m):
   '''
     residuum of isotropic Barlat 1991 yield criterion (eq. 2.37)
   '''
-  return Barlat1991(sigmas, sigma0, m, 1.0,1.0,1.0,1.0,1.0,1.0)
+  def fun(self, (sigma0, m), ydata, sigmas):
+    return Barlat1991Basis(sigma0, 1.0,1.0,1.0,1.0,1.0,1.0, m, sigmas)
+  def jac(self, (sigma0, m), ydata, sigmas):
+    pass
 
-def Barlat1991aniso(sigmas, sigma0, a,b,c,f,g,h, m):
+class Barlat1991aniso(object):
   '''
     residuum of anisotropic Barlat 1991 yield criterion (eq. 2.37)
   '''
-  return Barlat1991(sigmas, sigma0, m, a,b,c,f,g,h)
+  def fun(self, (sigma0, a,b,c,f,g,h, m), ydata, sigmas):
+    return Barlat1991Basis(sigma0, a,b,c,f,g,h, m, sigmas)
+  def jac(self, (sigma0, a,b,c,f,g,h, m), ydata, sigmas):
+    pass
 
-
+def Cazacu_Barlat3D(sigma0,a1,a2,a3,a4,a5,a6, b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11, c,
-def Barlat1994(sigmas, sigma0, a):
+                    ydata, sigmas):
-  '''
-    residuum of Hershey yield criterion (eq. 2.104, sigma_e = sigma0)
-  '''
-
-  return None
-
-def Cazacu_Barlat3D(sigmas, sigma0,
-                    a1,a2,a3,a4,a5,a6, b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11, c):
   '''
     residuum of the Cazacu–Barlat (CZ) yield criterion
   '''
   s11 = sigmas[0]; s22 = sigmas[1]; s33 = sigmas[2]
   s12 = sigmas[3]; s23 = sigmas[4]; s31 = sigmas[5]
-
+  
   J20 = ( a1*(s22-s33)**2 + a2*(s33-s11)**2 + a3*(s11-s22)**2 )/6.0 + \
           a4* s23**2      + a5* s31**2      + a6* s12**2
 
@@ -224,14 +195,13 @@ def Cazacu_Barlat3D(sigmas, sigma0,
           ( 2.0*b10*s33 - b5*s22  - (2*b10-b5)*s11 )*s12**2 +
           ( (b6+b7)*s11 - b6*s22  - b7*s33         )*s23**2
         )/3.0
-
   f0 = (J20**3 - c*J30**2)**(1.0/6.0)
   k2 = (sigma0/3.0) *18.0 **(1.0/6.0)
   r  = f0/k2 - 1.0
   return r.ravel()
 
-def Cazacu_Barlat2D(sigmas, sigma0,
+def Cazacu_Barlat2D(sigma0,a1,a2,a3,a6, b1,b2,b3,b4,b5,b10, c,
-                    a1,a2,a3,a6, b1,b2,b3,b4,b5,b10, c):
+                    ydata, sigmas):
   '''
     residuum of the Cazacu–Barlat (CZ) yield criterion for plain stress
   '''
@@ -242,13 +212,12 @@ def Cazacu_Barlat2D(sigmas, sigma0,
   J30 = ( (b1     + b2    )*s11**3  + (b3    +b4    )*s22**3 )/27.0- \
         ( (b1*s11 + b4*s22)*s11*s22  )/9.0 + \
         (  b5*s22 + (2*b10-b5)*s11 )*s12**2/3.0
-
   f0 = (J20**3 - c*J30**2)**(1.0/6.0)
   k2 = (sigma0/3.0) *18.0 **(1.0/6.0)
   r  = f0/k2 - 1.0
   return r.ravel()
 
-def BBC2003(sigmas, sigma0, a,b,c, d,e,f,g, k):
+def BBC2003(sigma0, a,b,c, d,e,f,g, k, ydata, sigmas):
   '''
     residuum of the BBC2003 yield criterion for plain stress
   '''
@@ -263,6 +232,84 @@ def BBC2003(sigmas, sigma0, a,b,c, d,e,f,g, k):
   r = sBar/sigma0 - 1.0
   return r.ravel()
 
+def DruckerBasis(sigma0, C_D, p, ydata, sigmas):
+  s11 = sigmas[0]; s22 = sigmas[1]; s33 = sigmas[2]
+  s12 = sigmas[3]; s23 = sigmas[4]; s31 = sigmas[5]
+  I1  = s11 + s22 + s33
+  I2  = s11*s22 + s22*s33 + s11*s33 - s12**2 - s23**2 - s31**2
+  I3  = s11*s22*s33 + 2.0*s12*s23*s31 - s12**2*s33 - s23**2*s11 - s31**2*s22
+  J2  = I1**2/3.0 - I2
+  J3  = I1**3/13.5 - I1*I2/3.0 + I3
+  r  = (J2**(3.0*p) - C_D*J3**(2.0*p))*27/(sigma0**6.0) - 1.0
+  return r.ravel()
+
+def HosfordBasis(sigma0, F,G,H, a, sigmas, Jac=False, nParas=1):
+  '''
+    residuum of Hershey yield criterion (eq. 2.43, Y = sigma0)
+  '''
+  lambdas = principalStresses(sigmas)
+  diff23  = abs(lambdas[1,:] - lambdas[2,:])
+  diff31  = abs(lambdas[2,:] - lambdas[0,:])
+  diff12  = abs(lambdas[0,:] - lambdas[1,:])
+  base    = F*diff23**a + G*diff31**a + H*diff12**a; expo = 1.0/a
+  left    = base**expo;    right = 2.0**expo*sigma0
+  
+  if not Jac:
+    if nParas == 1: return (left - right).ravel()
+    else:           return (left/right - 1.0).ravel()
+  else:
+    if nParas > 1:
+      ln   = lambda x : np.log(x + 1.0e-32)
+      dbda = F*ln(diff23)*diff23**a + G*ln(diff31)*diff31**a + H*ln(diff12)*diff12**a
+      deda = -expo*expo; dldb = expo*left/base; drda = sigma0*(2.0**expo)*ln(2.0)*deda
+      ones = np.ones(np.shape(sigmas)[1]);      jac  = []
+    
+    if   nParas == 1:     # von Mises
+      return ones*(-2.0**0.5)
+    elif nParas == 2:     # isotropic Hosford
+      j1 = ones*(-2.0**expo)                                              # d[]/dsigma0
+      j2 = dldb*dbda + left*ln(base)*deda - drda                          # d[]/da
+      for a,b in zip(j1, j2): jac.append([a,b])
+      return np.array(jac)
+    elif nParas == 5:           # anisotropic Hosford
+      j1 = -left/right/sigma0   #ones*(-2.0**expo)                                               # d[]/dsigma0
+      j2 = dldb*diff23**a/right; j3 = dldb*diff31**a/right; j4 = dldb*diff12**a/right
+      j5 =(dldb*dbda + left*ln(base)*deda)/right + left*(-right**(-2))*drda    # d[]/da
+      for a,b,c,d,e in zip(j1, j2,j3,j4,j5): jac.append([a,b,c,d,e])
+      return np.array(jac)
+      
+def Barlat1991Basis(sigma0, a, b, c, f, g, h, order, sigmas):
+  '''
+    residuum of Barlat 1997 yield criterion
+  '''
+  cos = np.cos; pi = np.pi; abs = np.abs
+  A = a*(sigmas[1] - sigmas[2])
+  B = b*(sigmas[2] - sigmas[0])
+  C = c*(sigmas[0] - sigmas[1])
+  F = f* sigmas[4]
+  G = g* sigmas[5]
+  H = h* sigmas[3]
+
+  I2 =   (F*F + G*G  +  H*H)/3.0  + ((A-C)**2+(C-B)**2+(B-A)**2)/54.0
+  I3 =   (C-B)*(A-C) * (B-A)/54.0 + F*G*H - \
+       ( (C-B)*F*F   + (A-C)*G*G  + (B-A)*H*H )/6.0
+  theta = np.arccos(I3/I2**1.5)
+  Phi = np.sqrt(3.0*I2)* (
+        (abs(2.0*cos((2.0*theta + pi)/6.0)))**order +
+        (abs(2.0*cos((2.0*theta + pi*3.0)/6.0)))**order +
+        (abs(2.0*cos(( 2.0*theta + pi*5.0)/6.0)))**order
+                         )**(1.0/order)
+#  r   = Phi/2.0**(1.0/order) - sigma0
+  r   = Phi/2.0**(1.0/order)/sigma0 - 1.0
+  # Phi = (3.0*I2)**(order/2.0) * (
+        # (abs(2.0*cos((2.0*theta + pi)/6.0)))    **order +
+        # (abs(2.0*cos((2.0*theta + pi*3.0)/6.0)))**order +
+        # (abs(2.0*cos((2.0*theta + pi*5.0)/6.0)))**order
+                                # )
+  # r   = (Phi - 2.0*sigma0**order)**(1.0/order)
+
+  return r.ravel()
+
 fittingCriteria = {
   'tresca'         :{'func' : Tresca,
                      'num'  : 1,'err':np.inf,
@@ -278,18 +325,25 @@ fittingCriteria = {
                      'text' : '\nCoefficient of Huber-Mises-Hencky criterion:\nsigma0: ',
                      'error': 'The standard deviation error is: '
                     },
-  'hosford'        :{'func' : Hosford,
+  'hosfordiso'     :{'func' : Hosford,
                      'num'  : 2,'err':np.inf,
-                     'name' : 'Gerenal Hosford',
+                     'name' : 'Gerenal isotropic Hosford',
-                     'paras': 'Initial yield stress:',
+                     'paras': 'Initial yield stress, a:',
                      'text' : '\nCoefficients of Hosford criterion:\nsigma0, a: ',
                      'error': 'The standard deviation errors are: '
                     },
+  'hosfordaniso'   :{'func' : generalHosford,
+                     'num'  : 5,'err':np.inf,
+                     'name' : 'Gerenal isotropic Hosford',
+                     'paras': 'Initial yield stress, F, G, H, a:',
+                     'text' : '\nCoefficients of Hosford criterion:\nsigma0, F, G, H, a: ',
+                     'error': 'The standard deviation errors are: '
+                    },
   'hill1948'       :{'func' : Hill1948,
                      'num'  : 6,'err':np.inf,
                      'name' : 'Hill1948',
-                     'paras': 'Normalized [F, G, H, L, M, N]',
+                     'paras': 'Normalized [F, G, H, L, M, N]:',
-                     'text' : '\nCoefficients of Hill1948 criterion:\n[F, G, H, L, M, N]:',
+                     'text' : '\nCoefficients of Hill1948 criterion:\n[F, G, H, L, M, N]:'+' '*16,
                      'error': 'The standard deviation errors are: '
                     },
   'drucker'        :{'func' : Drucker,
@@ -299,6 +353,13 @@ fittingCriteria = {
                      'text' : '\nCoefficients of Drucker criterion:\nsigma0, C_D: ',
                      'error': 'The standard deviation errors are: '
                     },
+  'gdrucker'       :{'func' : generalDrucker,
+                     'num'  : 3,'err':np.inf,
+                     'name' : 'General Drucker',
+                     'paras': 'Initial yield stress, C_D, p:',
+                     'text' : '\nCoefficients of Drucker criterion:\nsigma0, C_D, p: ',
+                     'error': 'The standard deviation errors are: '
+                    },
   'barlat1991iso'  :{'func' : Barlat1991iso,
                      'num'  : 2,'err':np.inf,
                      'name' : 'Barlat1991iso',
@@ -309,20 +370,20 @@ fittingCriteria = {
   'barlat1991aniso':{'func' : Barlat1991aniso,
                      'num'  : 8,'err':np.inf,
                      'name' : 'Barlat1991aniso',
-                     'paras': 'Initial yield stress, m, a, b, c, f, g, h:',
+                     'paras': 'Initial yield stress, a, b, c, f, g, h, m:',
                      'text' : '\nCoefficients of anisotropic Barlat 1991 criterion:\nsigma0, a, b, c, f, g, h, m:\n',
                      'error': 'The standard deviation errors are: '
                     },
   'bbc2003'        :{'func' : BBC2003,
                      'num'  : 9,'err':np.inf,
-                     'name' : 'Barlat1991aniso',
+                     'name' : 'Banabic-Balan-Comsa 2003',
                      'paras': 'Initial yield stress, a, b, c, d, e, f, g, k:',
                      'text' : '\nCoefficients of anisotropic Barlat 1991 criterion:\nsigma0, a, b, c, d, e, f, g, k:\n',
                      'error': 'The standard deviation errors are: '
                     },
   'Cazacu_Barlat2D':{'func' : Cazacu_Barlat2D,
                      'num'  : 12,'err':np.inf,
-                     'name' : 'Barlat1991aniso',
+                     'name' : 'Cazacu Barlat for plain stress',
                      'paras': 'Initial yield stress, a1,a2,a3,a6; b1,b2,b3,b4,b5,b10; c:',
                      'text' : '\nCoefficients of Cazacu Barlat yield criterion for plane stress: \
                                \n Y, a1,a2,a3,a6; b1,b2,b3,b4,b5,b10; c:\n',
@@ -330,7 +391,7 @@ fittingCriteria = {
                     },
   'Cazacu_Barlat3D':{'func' : Cazacu_Barlat3D,
                      'num'  : 19,'err':np.inf,
-                     'name' : 'Barlat1991aniso',
+                     'name' : 'Cazacu Barlat',
                      'paras': 'Initial yield stress, a1,a2,a3,a4,a5,a6; b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11; c:',
                      'text' : '\nCoefficients of Cazacu Barlat yield criterion for plane stress: \
                                \n Y, a1,a2,a3,a4,a5,a6; b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11; c\n',
@@ -393,7 +454,7 @@ class Criterion(object):
   def __init__(self,name='worst'):
     self.name = name
     self.results = fittingCriteria
-    
+
     if self.name.lower() not in map(str.lower, self.results.keys()):
       raise Exception('no suitable fitting criterion selected')
     else:
@@ -401,9 +462,9 @@ class Criterion(object):
 
   def fit(self,stress):
     global fitResults
-    
+
     nameCriterion = self.name.lower()
-    funResidum    = fittingCriteria[nameCriterion]['func']
+    criteriaClass = fittingCriteria[nameCriterion]['func'];  criteria = criteriaClass()
     numParas      = fittingCriteria[nameCriterion]['num']
     textParas     = fittingCriteria[nameCriterion]['text'] + formatOutput(numParas)
     textError     = fittingCriteria[nameCriterion]['error']+ formatOutput(numParas,'%-14.8f')+'\n'
@@ -413,14 +474,21 @@ class Criterion(object):
     if fitResults == [] : initialguess = guess0
     else                : initialguess = np.array(fitResults[-1])
     weight = get_weight(np.shape(stress)[1])
+    ydata  = np.zeros(np.shape(stress)[1])
     try:
-      popt, pcov = \
+      popt, pcov, infodict, errmsg, ierr = \
-        curve_fit_bound(funResidum, stress, np.zeros(np.shape(stress)[1]),
+         leastsqBound (criteria.fun,  initialguess,        args=(ydata,stress),
-                        initialguess, weight, bounds)
+                       bounds=bounds, full_output=True)
+      if ierr not in [1, 2, 3, 4]: raise RuntimeError("Optimal parameters not found: " + errmsg)
+      if (len(ydata) > len(initialguess)) and pcov is not None:
+        s_sq = (criteria.fun(popt, *(ydata,stress))**2).sum()/(len(ydata)-len(initialguess))
+        pcov = pcov * s_sq
       perr = np.sqrt(np.diag(pcov))
       fitResults.append(popt.tolist())
+
       print (textParas%array2tuple(popt))
       print (textError%array2tuple(perr))
+      print('Number of function calls =', infodict['nfev'])
     except Exception as detail:
       print detail
       pass