Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Exact Hessians with SymEngine lambda backend #249

Merged
merged 36 commits into from Nov 4, 2019
Merged

Conversation

@richardotis
Copy link
Collaborator

richardotis commented Oct 28, 2019

Fixes gh-235.

richardotis added 23 commits Jun 11, 2019
… for performance reasons
…prevent garbage collection"

This reverts commit 0d3c9c0.
…conditions, for now
@richardotis richardotis added this to the 0.8.1 milestone Oct 28, 2019
@richardotis richardotis requested a review from bocklund Oct 28, 2019
@bocklund

This comment has been minimized.

Copy link
Collaborator

bocklund commented Oct 28, 2019

I have found some odd behavior when different sets of parameters are provided:

tdb = """$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
$ Date: 2018-12-18 11:03
$ Components: CU, MG, VA
$ Phases: CUMG2, FCC_A1, HCP_A3, LAVES_C15, LIQUID
$ Generated by brandon (pycalphad 0.7.1.post2)
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

ELEMENT CU BLANK 0 0 0 !
ELEMENT MG BLANK 0 0 0 !
ELEMENT VA BLANK 0 0 0 !


FUNCTION GFCCCU 1.0 GHSERCU#; 10000.0 N !
FUNCTION GFCCMG 0.01 -2.18283562784578E-6*T**3 - 6.21802726222479E-5*T**2 +
   3.2404446864*T*LN(EXP(95.82831*T**(-1)) - 1.0) +
   21.6535319868*T*LN(EXP(247.8675*T**(-1)) - 1.0) + 0.246565697987779*T -
   8158.16393259455; 298.15 Y -1.393669E-6*T**3 + 0.0004858*T**2 -
   26.1849782*T*LN(T) + 142.775547*T - 5767.34 + 78950*T**(-1);
   922.205302616508 Y -34.3088*T*LN(T) + 203.816215*T - 11530.1836392866 +
   1.038192E+28*T**(-9); 3000.0 N !
FUNCTION GHCPCU 0.01 -3.38438862938597E-7*T**3 - 0.00121182291077191*T**2 +
   8.7685671186*T*LN(EXP(155.1404*T**(-1)) - 1.0) +
   16.1968683846*T*LN(EXP(290.9421*T**(-1)) - 1.0) - 0.321147237334052*T -
   10441.4393392344; 298.15 Y 1.29223E-7*T**3 - 0.00265684*T**2 -
   24.112392*T*LN(T) + 130.685235*T - 7170.458 + 52478*T**(-1); 1357.77 Y
   -31.38*T*LN(T) + 184.003828*T - 12942.0252504739 + 3.64167E+29*T**(-9);
   3200.0 N !
FUNCTION GHCPMG 1.0 GHSERMG#; 10000.0 N !
FUNCTION GHSERCU 0.01 -0.0010514335*T**2 +
   8.7685671186*T*LN(EXP(155.1404*T**(-1)) - 1.0) +
   16.1968683846*T*LN(EXP(290.9421*T**(-1)) - 1.0) - 11038.0904080745;
   103.57591 Y -2.15621953171362E-6*T**3 + 0.000288560900942072*T**2 -
   0.13879113947248*T*LN(T) + 8.7685671186*T*LN(EXP(155.1404*T**(-1)) - 1.0)
   + 16.1968683846*T*LN(EXP(290.9421*T**(-1)) - 1.0) + 0.574637617323048*T -
   11042.8822142647; 210.33309 Y -0.002432585*T**2 + 0.4335558862135*T*LN(T)
   + 8.7685671186*T*LN(EXP(155.1404*T**(-1)) - 1.0) +
   16.1968683846*T*LN(EXP(290.9421*T**(-1)) - 1.0) - 2.20049706600083*T -
   11002.7543747764; 1357.77 Y -31.38*T*LN(T) + 183.555483717662*T -
   12730.2995781851 + 7.42232714807953E+28*T**(-9); 3200.0 N !
FUNCTION GHSERMG 0.01 0.0047195465*T**2 +
   3.2404446864*T*LN(EXP(95.82831*T**(-1)) - 1.0) +
   21.6535319868*T*LN(EXP(247.8675*T**(-1)) - 1.0) - 10652.1012810789;
   36.71926 Y -1.53643054262276E-5*T**3 + 0.00810454205399037*T**2 -
   0.124294531845816*T*LN(T) + 3.2404446864*T*LN(EXP(95.82831*T**(-1)) - 1.0)
   + 21.6535319868*T*LN(EXP(247.8675*T**(-1)) - 1.0) + 0.385723396310737*T -
   10653.6226154894; 143.18844 Y -0.0050954035*T**2 + 1.765785080115*T*LN(T)
   + 3.2404446864*T*LN(EXP(95.82831*T**(-1)) - 1.0) +
   21.6535319868*T*LN(EXP(247.8675*T**(-1)) - 1.0) - 8.0518972866125*T -
   10563.4100984519; 922.205302616508 Y -34.3088*T*LN(T) + 204.341485347575*T
   - 13775.4156328263 + 9.4687256586798E+27*T**(-9); 10000.0 N !
FUNCTION GHSERVA 1 0; 10000 N !
FUNCTION GLIQCU 0.01 -3.40056501515466E-7*T**3 - 0.00121066539331185*T**2 +
   8.7685671186*T*LN(EXP(155.1404*T**(-1)) - 1.0) +
   16.1968683846*T*LN(EXP(290.9421*T**(-1)) - 1.0) - 10.033338832193*T +
   2379.36422209194; 298.15 Y -5.8489E-21*T**7 + 1.29223E-7*T**3 -
   0.00265684*T**2 - 24.112392*T*LN(T) + 120.973331*T + 5650.32106235287 +
   52478*T**(-1); 1357.77 Y -31.38*T*LN(T) + 173.881484*T + 409.498458129716;
   3200.0 N !
FUNCTION GLIQMG 0.01 -2.2050100179942E-6*T**3 - 4.63131660076452E-5*T**2 +
   3.2404446864*T*LN(EXP(95.82831*T**(-1)) - 1.0) +
   21.6535319868*T*LN(EXP(247.8675*T**(-1)) - 1.0) - 7.6943066168765*T -
   2555.58510336379; 298.15 Y -8.0176E-20*T**7 - 1.393669E-6*T**3 +
   0.0004858*T**2 - 26.1849782*T*LN(T) + 134.838617*T - 165.096999676889 +
   78950*T**(-1); 922.205302616508 Y -34.3088*T*LN(T) + 195.324057*T -
   5439.86911093575; 10000.0 N !
FUNCTION VV0000 1 -32539.5; 10000 N !
FUNCTION VV0001 1 8236.3; 10000 N !
FUNCTION VV0002 1 -14675.0; 10000 N !
FUNCTION VV0003 1 -24441.2; 10000 N !
FUNCTION VV0004 1 20149.6; 10000 N !
FUNCTION VV0005 1 46500.0; 10000 N !
FUNCTION VV0006 1 -39591.3; 10000 N !
FUNCTION VV0007 1 104160.0; 10000 N !
FUNCTION VV0008 1 21000.0; 10000 N !
FUNCTION VV0009 1 17772.0; 10000 N !
FUNCTION VV0010 1 21240.0; 10000 N !
FUNCTION VV0011 1 14321.1; 10000 N !
FUNCTION VV0012 1 -4923.18; 10000 N !
FUNCTION VV0013 1 -1962.8; 10000 N !
FUNCTION VV0014 1 -31626.6; 10000 N !

TYPE_DEFINITION % SEQ * !
DEFINE_SYSTEM_DEFAULT ELEMENT 2 !
DEFAULT_COMMAND DEFINE_SYSTEM_ELEMENT VA !

PHASE CUMG2 %  2 1 2 !
CONSTITUENT CUMG2 :CU:MG: !

PHASE FCC_A1 %  1 1 !
CONSTITUENT FCC_A1 :CU,MG: !

PHASE HCP_A3 %  1 1 !
CONSTITUENT HCP_A3 :CU,MG: !

PHASE LAVES_C15 %  2 2 1 !
CONSTITUENT LAVES_C15 :CU,MG:CU,MG: !

PHASE LIQUID %  1 1 !
CONSTITUENT LIQUID :CU,MG: !



$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
$                                     CU                                     $
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

PARAMETER G(FCC_A1,CU;0) 1 GFCCCU#; 10000 N !
PARAMETER G(HCP_A3,CU;0) 1 GHCPCU#; 10000 N !
PARAMETER G(LAVES_C15,CU:CU;0) 1 3*GHSERCU# + VV0005#; 10000 N !
PARAMETER G(LIQUID,CU;0) 1 GLIQCU#; 10000 N !


$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
$                                     MG                                     $
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

PARAMETER G(FCC_A1,MG;0) 1 GFCCMG#; 10000 N !
PARAMETER G(HCP_A3,MG;0) 1 GHCPMG#; 10000 N !
PARAMETER G(LAVES_C15,MG:MG;0) 1 3*GHSERMG# + VV0008#; 10000 N !
PARAMETER G(LIQUID,MG;0) 1 GLIQMG#; 10000 N !


$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
$                                   CU-MG                                    $
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

PARAMETER G(CUMG2,CU:MG;0) 1 GHSERCU# + 2*GHSERMG# + VV0000#; 10000 N !
PARAMETER L(FCC_A1,CU,MG;0) 1 VV0002#; 10000 N !
PARAMETER L(FCC_A1,CU,MG;1) 1 VV0001#; 10000 N !
PARAMETER L(HCP_A3,CU,MG;0) 1 VV0004#; 10000 N !
PARAMETER L(HCP_A3,CU,MG;1) 1 VV0003#; 10000 N !
PARAMETER G(LAVES_C15,CU:MG;0) 1 2*GHSERCU# + GHSERMG# + VV0006#; 10000 N !
PARAMETER G(LAVES_C15,MG:CU;0) 1 GHSERCU# + 2*GHSERMG# + VV0007#; 10000 N !
PARAMETER L(LAVES_C15,CU:CU,MG;0) 1 VV0009#; 10000 N !
PARAMETER L(LAVES_C15,CU,MG:MG;0) 1 VV0010#; 10000 N !
PARAMETER L(LIQUID,CU,MG;0) 1 VV0014#; 10000 N !
PARAMETER L(LIQUID,CU,MG;1) 1 VV0013#; 10000 N !
PARAMETER L(LIQUID,CU,MG;2) 1 VV0012#; 10000 N !
PARAMETER L(LIQUID,CU,MG;3) 1 VV0011#; 10000 N !
"""

from pycalphad import Database, equilibrium, variables as v, __version__ as pycalphad_version
print(f'pycalphad version: {pycalphad_version}')
dbf= Database(tdb)
comps = ['CU', 'MG']

parameters={'VV0000': -33134.699474175846, 'VV0001': 7734.114029426941, 'VV0002': -13498.542175596054, 'VV0003': -26555.048975092268, 'VV0004': 20777.637577083482, 'VV0005': 41915.70425630003, 'VV0006': -34525.21964215504, 'VV0007': 95457.14639216446, 'VV0008': 21139.578967453144, 'VV0009': 19047.833726419598, 'VV0010': 20468.91829601273, 'VV0011': 19601.617855958328, 'VV0012': -4546.9325861738, 'VV0013': -1640.6354331231278, 'VV0014': -35682.950005357634}
print('1 run: all parameters')
# I am getting "Status: -2 b"Restoration phase failed, algorithm doesn't know how to proceed.""
eq = equilibrium(dbf, comps, ['HCP_A3'], {v.X('CU'): 0.0001052, v.P: 101325.0, v.T: 743.15, v.N: 1}, parameters=parameters, verbose=True)
print('\n\n2 run: some parameters')
# I am getting "Status: -1 b'Maximum number of iterations exceeded (can be specified by an option).'"
eq = equilibrium(dbf, comps, ['HCP_A3'], {v.X('CU'): 0.0001052, v.P: 101325.0, v.T: 743.15, v.N: 1}, parameters={'VV0003': -26555.048975092268, 'VV0004': 20777.637577083482}, verbose=True)
print('\n\n3 run: parameters pre-subbed')
dbf.symbols.update(parameters)
# I am getting "Status: 0 b'Algorithm terminated successfully at a locally optimal point, satisfying the convergence tolerances (can be specified by options).'"
eq = equilibrium(dbf, comps, ['HCP_A3'], {v.X('CU'): 0.0001052, v.P: 101325.0, v.T: 743.15, v.N: 1}, verbose=True)

I get the output:

pycalphad version: 0.8+27.g19bf29a6
1 run: all parameters
<snip>
Status: -2 b"Restoration phase failed, algorithm doesn't know how to proceed."


2 run: some parameters
<snip>
Status: -1 b'Maximum number of iterations exceeded (can be specified by an option).'


3 run: parameters pre-subbed
<snip>
Status: 0 b'Algorithm terminated successfully at a locally optimal point, satisfying the convergence tolerances (can be specified by options).'

This is interesting because the HCP_A3 only has VV0003 and VV0004 parameters in it. VV0003 and VV0004 are the same across all versions, but something about having other parameters (or no parameters) passed is messing with the solution.

@bocklund

This comment has been minimized.

Copy link
Collaborator

bocklund commented Oct 28, 2019

Also, looks like the test failures on Travis's PR build are caused by NumPy Python 3 syntax for matrix multiplication (A @ B).

Based on our release cadence, I expect that this will be the last release where our upstream dependencies are officially supporting Python 2.

@bocklund

This comment has been minimized.

Copy link
Collaborator

bocklund commented Oct 28, 2019

I compared all the outputs of phase_record for a dof array corresponding to the starting point of the equilibrium calculations above. Prx1 and Prx2 have identical entries for VV0003 and VV0004, but different entries for other parameters. Since, VV0003 and VV0004 are the only parameters that show up in the Gibbs energy of HCP_A3, the phase records outputs should all be equal. This script confirms that they are, so something else must be causing the discrepancy in the different equilibrium calls (above). I'm thinking it could be either

  1. My phase records are not built consistently with how equilibrium is building them (I believe they are, but my logic could use a second look)
  2. Something happening in how Problem is putting together the outputs from the phase records with the two composition set starting point is off

Self contained MWE:

tdb = """$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
$ Date: 2018-12-18 11:03
$ Components: CU, MG, VA
$ Phases: CUMG2, FCC_A1, HCP_A3, LAVES_C15, LIQUID
$ Generated by brandon (pycalphad 0.7.1.post2)
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

ELEMENT CU BLANK 0 0 0 !
ELEMENT MG BLANK 0 0 0 !
ELEMENT VA BLANK 0 0 0 !


FUNCTION GFCCCU 1.0 GHSERCU#; 10000.0 N !
FUNCTION GFCCMG 0.01 -2.18283562784578E-6*T**3 - 6.21802726222479E-5*T**2 +
   3.2404446864*T*LN(EXP(95.82831*T**(-1)) - 1.0) +
   21.6535319868*T*LN(EXP(247.8675*T**(-1)) - 1.0) + 0.246565697987779*T -
   8158.16393259455; 298.15 Y -1.393669E-6*T**3 + 0.0004858*T**2 -
   26.1849782*T*LN(T) + 142.775547*T - 5767.34 + 78950*T**(-1);
   922.205302616508 Y -34.3088*T*LN(T) + 203.816215*T - 11530.1836392866 +
   1.038192E+28*T**(-9); 3000.0 N !
FUNCTION GHCPCU 0.01 -3.38438862938597E-7*T**3 - 0.00121182291077191*T**2 +
   8.7685671186*T*LN(EXP(155.1404*T**(-1)) - 1.0) +
   16.1968683846*T*LN(EXP(290.9421*T**(-1)) - 1.0) - 0.321147237334052*T -
   10441.4393392344; 298.15 Y 1.29223E-7*T**3 - 0.00265684*T**2 -
   24.112392*T*LN(T) + 130.685235*T - 7170.458 + 52478*T**(-1); 1357.77 Y
   -31.38*T*LN(T) + 184.003828*T - 12942.0252504739 + 3.64167E+29*T**(-9);
   3200.0 N !
FUNCTION GHCPMG 1.0 GHSERMG#; 10000.0 N !
FUNCTION GHSERCU 0.01 -0.0010514335*T**2 +
   8.7685671186*T*LN(EXP(155.1404*T**(-1)) - 1.0) +
   16.1968683846*T*LN(EXP(290.9421*T**(-1)) - 1.0) - 11038.0904080745;
   103.57591 Y -2.15621953171362E-6*T**3 + 0.000288560900942072*T**2 -
   0.13879113947248*T*LN(T) + 8.7685671186*T*LN(EXP(155.1404*T**(-1)) - 1.0)
   + 16.1968683846*T*LN(EXP(290.9421*T**(-1)) - 1.0) + 0.574637617323048*T -
   11042.8822142647; 210.33309 Y -0.002432585*T**2 + 0.4335558862135*T*LN(T)
   + 8.7685671186*T*LN(EXP(155.1404*T**(-1)) - 1.0) +
   16.1968683846*T*LN(EXP(290.9421*T**(-1)) - 1.0) - 2.20049706600083*T -
   11002.7543747764; 1357.77 Y -31.38*T*LN(T) + 183.555483717662*T -
   12730.2995781851 + 7.42232714807953E+28*T**(-9); 3200.0 N !
FUNCTION GHSERMG 0.01 0.0047195465*T**2 +
   3.2404446864*T*LN(EXP(95.82831*T**(-1)) - 1.0) +
   21.6535319868*T*LN(EXP(247.8675*T**(-1)) - 1.0) - 10652.1012810789;
   36.71926 Y -1.53643054262276E-5*T**3 + 0.00810454205399037*T**2 -
   0.124294531845816*T*LN(T) + 3.2404446864*T*LN(EXP(95.82831*T**(-1)) - 1.0)
   + 21.6535319868*T*LN(EXP(247.8675*T**(-1)) - 1.0) + 0.385723396310737*T -
   10653.6226154894; 143.18844 Y -0.0050954035*T**2 + 1.765785080115*T*LN(T)
   + 3.2404446864*T*LN(EXP(95.82831*T**(-1)) - 1.0) +
   21.6535319868*T*LN(EXP(247.8675*T**(-1)) - 1.0) - 8.0518972866125*T -
   10563.4100984519; 922.205302616508 Y -34.3088*T*LN(T) + 204.341485347575*T
   - 13775.4156328263 + 9.4687256586798E+27*T**(-9); 10000.0 N !
FUNCTION GHSERVA 1 0; 10000 N !
FUNCTION GLIQCU 0.01 -3.40056501515466E-7*T**3 - 0.00121066539331185*T**2 +
   8.7685671186*T*LN(EXP(155.1404*T**(-1)) - 1.0) +
   16.1968683846*T*LN(EXP(290.9421*T**(-1)) - 1.0) - 10.033338832193*T +
   2379.36422209194; 298.15 Y -5.8489E-21*T**7 + 1.29223E-7*T**3 -
   0.00265684*T**2 - 24.112392*T*LN(T) + 120.973331*T + 5650.32106235287 +
   52478*T**(-1); 1357.77 Y -31.38*T*LN(T) + 173.881484*T + 409.498458129716;
   3200.0 N !
FUNCTION GLIQMG 0.01 -2.2050100179942E-6*T**3 - 4.63131660076452E-5*T**2 +
   3.2404446864*T*LN(EXP(95.82831*T**(-1)) - 1.0) +
   21.6535319868*T*LN(EXP(247.8675*T**(-1)) - 1.0) - 7.6943066168765*T -
   2555.58510336379; 298.15 Y -8.0176E-20*T**7 - 1.393669E-6*T**3 +
   0.0004858*T**2 - 26.1849782*T*LN(T) + 134.838617*T - 165.096999676889 +
   78950*T**(-1); 922.205302616508 Y -34.3088*T*LN(T) + 195.324057*T -
   5439.86911093575; 10000.0 N !
FUNCTION VV0000 1 -32539.5; 10000 N !
FUNCTION VV0001 1 8236.3; 10000 N !
FUNCTION VV0002 1 -14675.0; 10000 N !
FUNCTION VV0003 1 -24441.2; 10000 N !
FUNCTION VV0004 1 20149.6; 10000 N !
FUNCTION VV0005 1 46500.0; 10000 N !
FUNCTION VV0006 1 -39591.3; 10000 N !
FUNCTION VV0007 1 104160.0; 10000 N !
FUNCTION VV0008 1 21000.0; 10000 N !
FUNCTION VV0009 1 17772.0; 10000 N !
FUNCTION VV0010 1 21240.0; 10000 N !
FUNCTION VV0011 1 14321.1; 10000 N !
FUNCTION VV0012 1 -4923.18; 10000 N !
FUNCTION VV0013 1 -1962.8; 10000 N !
FUNCTION VV0014 1 -31626.6; 10000 N !

TYPE_DEFINITION % SEQ * !
DEFINE_SYSTEM_DEFAULT ELEMENT 2 !
DEFAULT_COMMAND DEFINE_SYSTEM_ELEMENT VA !

PHASE CUMG2 %  2 1 2 !
CONSTITUENT CUMG2 :CU:MG: !

PHASE FCC_A1 %  1 1 !
CONSTITUENT FCC_A1 :CU,MG: !

PHASE HCP_A3 %  1 1 !
CONSTITUENT HCP_A3 :CU,MG: !

PHASE LAVES_C15 %  2 2 1 !
CONSTITUENT LAVES_C15 :CU,MG:CU,MG: !

PHASE LIQUID %  1 1 !
CONSTITUENT LIQUID :CU,MG: !



$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
$                                     CU                                     $
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

PARAMETER G(FCC_A1,CU;0) 1 GFCCCU#; 10000 N !
PARAMETER G(HCP_A3,CU;0) 1 GHCPCU#; 10000 N !
PARAMETER G(LAVES_C15,CU:CU;0) 1 3*GHSERCU# + VV0005#; 10000 N !
PARAMETER G(LIQUID,CU;0) 1 GLIQCU#; 10000 N !


$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
$                                     MG                                     $
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

PARAMETER G(FCC_A1,MG;0) 1 GFCCMG#; 10000 N !
PARAMETER G(HCP_A3,MG;0) 1 GHCPMG#; 10000 N !
PARAMETER G(LAVES_C15,MG:MG;0) 1 3*GHSERMG# + VV0008#; 10000 N !
PARAMETER G(LIQUID,MG;0) 1 GLIQMG#; 10000 N !


$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
$                                   CU-MG                                    $
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

PARAMETER G(CUMG2,CU:MG;0) 1 GHSERCU# + 2*GHSERMG# + VV0000#; 10000 N !
PARAMETER L(FCC_A1,CU,MG;0) 1 VV0002#; 10000 N !
PARAMETER L(FCC_A1,CU,MG;1) 1 VV0001#; 10000 N !
PARAMETER L(HCP_A3,CU,MG;0) 1 VV0004#; 10000 N !
PARAMETER L(HCP_A3,CU,MG;1) 1 VV0003#; 10000 N !
PARAMETER G(LAVES_C15,CU:MG;0) 1 2*GHSERCU# + GHSERMG# + VV0006#; 10000 N !
PARAMETER G(LAVES_C15,MG:CU;0) 1 GHSERCU# + 2*GHSERMG# + VV0007#; 10000 N !
PARAMETER L(LAVES_C15,CU:CU,MG;0) 1 VV0009#; 10000 N !
PARAMETER L(LAVES_C15,CU,MG:MG;0) 1 VV0010#; 10000 N !
PARAMETER L(LIQUID,CU,MG;0) 1 VV0014#; 10000 N !
PARAMETER L(LIQUID,CU,MG;1) 1 VV0013#; 10000 N !
PARAMETER L(LIQUID,CU,MG;2) 1 VV0012#; 10000 N !
PARAMETER L(LIQUID,CU,MG;3) 1 VV0011#; 10000 N !
"""


def calc_output(phase_rec, verbose=False):
    dof = np.array([1.0, 101325.0, 743.15, 0.79759519, 0.20240481])
    nspecies = len(phase_rec.components)
    wrt_size = len(phase_rec.state_variables) + phase_rec.phase_dof
    int_cons_size = phase_rec.num_internal_cons
    mp_cons_size = phase_rec.num_multiphase_cons

    # Placeholder arrays
    obj1 = np.zeros(1)
    grad1 = np.zeros(wrt_size)
    hess1 = np.zeros((wrt_size, wrt_size))
    mass_obj1 = np.zeros(1)
    mass_grad1 = np.zeros(wrt_size)
    mass_hess1 = np.zeros((wrt_size, wrt_size))
    internal_obj1 = np.zeros(int_cons_size)
    internal_grad1 = np.zeros((wrt_size, int_cons_size))
    internal_hess1 = np.zeros((wrt_size, wrt_size, int_cons_size))
    multiphase_obj1 = np.zeros(mp_cons_size)
    multiphase_grad1 = np.zeros((wrt_size + 1, mp_cons_size)) # + 1 because NP variable
    # Calculate
    ## Objectives
    if verbose:
      print("Objectives...", end=' ')
    phase_rec.obj(obj1, np.atleast_2d(dof))
    phase_rec.grad(grad1, dof)
    phase_rec.hess(hess1, dof)
    ## Mass
    if verbose:
      print("Mass...", end=' ')
    mass_obj1 =  np.array([(phase_rec.mass_obj(mass_obj1, np.atleast_2d(dof), i), copy.deepcopy(mass_obj1))[1] for i in range(nspecies)])
    mass_grad1 = np.array([(phase_rec.mass_grad(mass_grad1, dof, i), copy.copy(mass_grad1))[1] for i in range(nspecies)])
    mass_hess1 = np.array([(phase_rec.mass_hess(mass_hess1, dof, i), copy.copy(mass_hess1))[1] for i in range(nspecies)])
    ## Internal constraints
    if verbose:
      print("Internal constraints...", end=' ')
    phase_rec.internal_constraints(internal_obj1, dof)
    phase_rec.internal_jacobian(internal_grad1, dof)
    phase_rec.internal_cons_hessian(internal_hess1, dof)
    ## Multiphase constraints
    mpdof = np.concatenate([dof, [1.0]], axis=0)
    if verbose:
      print("Multiphase constraints...")
      print(mpdof)
    phase_rec.multiphase_constraints(multiphase_obj1, mpdof)
    phase_rec.multiphase_jacobian(multiphase_grad1, mpdof)
    # Results
    if verbose:
      print(f"objective {obj1}")
      print(f"gradient {grad1}")
      print(f"hessian {hess1}")
      print(f"mass objective {mass_obj1}")
      print(f"mass gradient {mass_grad1}")
      print(f"mass hessian {mass_hess1}")
      print(f"internal objective {internal_obj1}")
      print(f"internal gradient {internal_grad1}")
      print(f"internal hessian {internal_hess1}")
      print(f"multiphase objective {multiphase_obj1}")
      print(f"multiphase gradient {multiphase_grad1}")
    return obj1, grad1, hess1, mass_obj1, mass_grad1, mass_hess1, internal_obj1, internal_grad1, internal_hess1, multiphase_obj1, multiphase_grad1

from pycalphad import Database, equilibrium, variables as v, __version__ as pycalphad_version
from pycalphad.core.utils import instantiate_models
from pycalphad.codegen.callables import build_phase_records, unpack_components
from pycalphad import Model
import numpy as np
import copy

print(f'pycalphad version: {pycalphad_version}')
dbf= Database(tdb)
comps = ['CU', 'MG']
phases = ['HCP_A3']
species = unpack_components(dbf, comps)

# Lots of parameters
parameters = {'VV0000': -33134.699474175846, 'VV0001': 7734.114029426941, 'VV0002': -13498.542175596054, 'VV0003': -26555.048975092268, 'VV0004': 20777.637577083482, 'VV0005': 41915.70425630003, 'VV0006': -34525.21964215504, 'VV0007': 95457.14639216446, 'VV0008': 21139.578967453144, 'VV0009': 19047.833726419598, 'VV0010': 20468.91829601273, 'VV0011': 19601.617855958328, 'VV0012': -4546.9325861738, 'VV0013': -1640.6354331231278, 'VV0014': -35682.950005357634}
models = instantiate_models(dbf, comps, phases, parameters=parameters)
prxs = build_phase_records(dbf, species, ['HCP_A3'], {v.X('CU'): 0.0001052, v.P: 101325.0, v.T: 743.15, v.N: 1}, models, parameters=parameters, build_gradients=True, build_hessians=True)
prx1 = prxs['HCP_A3']
outp1 = calc_output(prx1)


# Only parameters in this phase
parameters={'VV0003': -26555.048975092268, 'VV0004': 20777.637577083482}
models = instantiate_models(dbf, comps, phases, parameters=parameters)
prxs = build_phase_records(dbf, species, ['HCP_A3'], {v.X('CU'): 0.0001052, v.P: 101325.0, v.T: 743.15, v.N: 1}, models, parameters=parameters, build_gradients=True, build_hessians=True)
prx2 = prxs['HCP_A3']
outp2 = calc_output(prx2)

qty_names = ['Obj', 'Jac', 'Hess', 'Mass Obj', 'Mass Jac', 'Mass Hess', 'Internal Cons Obj', 'Internal Cons Jac', 'Internal Cons Hess', 'Multiphase Cons Obj', 'Multiphase Cons Jac']
print("Comparing prx1 and prx2 (expected to match)")
num_matches = 0
for qty, o1, o2 in zip(qty_names,  outp1, outp2):
    if not np.all(o1 == o2):
        num_matches += 1
        print(f"{qty} not matched ", o1, o2)
else:
    if num_matches == 0:
      print("All outputs match")

print()
# Sanity check, compare prx2 and prx3 where 3 has different parameters that should affect the values
parameters={'VV0003': 1e5, 'VV0004': 1e5}
prxs = build_phase_records(dbf, species, ['HCP_A3'], {v.X('CU'): 0.0001052, v.P: 101325.0, v.T: 743.15, v.N: 1}, {'HCP_A3': Model(dbf, comps, 'HCP_A3', parameters=sorted(parameters.keys()))}, parameters=parameters, build_gradients=True, build_hessians=True)
prx3 = prxs['HCP_A3']
outp3 = calc_output(prx3)

print("Comparing prx2 and prx3 (expected to not match)")
num_matches = 0
for qty, o1, o2 in zip(qty_names,  outp2, outp3):
    if not np.all(o1 == o2):
        num_matches += 1
        print(f"{qty} not matched ", o1, o2)
else:
    if num_matches == 0:
      print("All outputs match")

with output:

pycalphad version: 0.8+27.g19bf29a6
Comparing prx1 and prx2 (expected to match)
All outputs match

Comparing prx2 and prx3 (expected to not match)
Obj not matched  [-32270.97962691] [-7321.40648521]
Jac not matched  [   0.            0.          -61.03871068 3915.09258244 6459.41400703] [ 0.00000000e+00  0.00000000e+00 -6.10387107e+01  3.06771974e+04
  8.43448730e+04]
Hess not matched  [[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00 -3.77598461e-02  1.08113490e+01
  -1.52456842e+00]
 [ 0.00000000e+00  0.00000000e+00  1.08113490e+01 -1.08329862e+04
  -2.12074884e+04]
 [ 0.00000000e+00  0.00000000e+00 -1.52456842e+00 -2.12074884e+04
   5.99690694e+04]] [[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00 -3.77598461e-02  1.08113490e+01
  -1.52456842e+00]
 [ 0.00000000e+00  0.00000000e+00  1.08113490e+01 -1.31264946e+04
   1.04016006e+05]
 [ 0.00000000e+00  0.00000000e+00 -1.52456842e+00  1.04016006e+05
  -2.97681245e+05]]
richardotis added 4 commits Oct 29, 2019
…rs. This change corrects a case where dof and parameter arrays could be misaligned during concatenation, when the caller specified an oversized dof array. See gh-249 for details.
@richardotis

This comment has been minimized.

Copy link
Collaborator Author

richardotis commented Oct 29, 2019

The bug discussed above should be fixed now. To understand the issue, a little background:
The PhaseRecord object is designed to insulate callers (notably Problem) from having to know about model parameters (VV0001, VV0002, etc.). The caller specifies the degrees of freedom (dof) array with state variables and site fractions, which then get passed to the underlying callables for the objective, gradient, Hessian, constraints, etc.

Now, it's important to note that the underlying callables require the input array (dof+parameters) to be contiguous in memory. The PhaseRecord object has two helper functions, alloc_dof_with_parameters and alloc_dof_with_parameters_vectorized, which its primary functions call to guarantee that the input array is contiguous in memory. Since the dof array is already contiguous by construction, these helpers only allocate new memory if the parameters array is nonzero in size, which is why this bug only shows up in that case.

The other important fact to know is that, in some parts of the Problem class, oversized dof arrays are allocated and reused, for performance reasons. At this point you're probably seeing the issue: PhaseRecord is concatenating an oversized dof array with the parameters array, causing the input array to be contiguous, but misaligned, with padding zeros crowding out the parameters. The solution is simple: slice the dof array to be the known correct size prior to concatenation.

In terms of impact, by luck it turns out we were using correctly-sized dof arrays in the objective, gradient, as well as all the chemical potential calculation code, so, if the system converged, the answers were correct. (The Hessian and certain constraint code use oversized dof arrays.)

Copy link
Collaborator

bocklund left a comment

Exciting performance and accuracy improvements here. Great work!

pycalphad/codegen/sympydiff_utils.py Show resolved Hide resolved
pycalphad/codegen/sympydiff_utils.py Show resolved Hide resolved
pycalphad/core/constraints.py Outdated Show resolved Hide resolved
pycalphad/core/constraints.py Outdated Show resolved Hide resolved
pycalphad/core/phase_rec.pyx Outdated Show resolved Hide resolved
pycalphad/core/solver.py Outdated Show resolved Hide resolved
pycalphad/model.py Outdated Show resolved Hide resolved
pycalphad/model.py Outdated Show resolved Hide resolved
pycalphad/model.py Outdated Show resolved Hide resolved
setup.py Show resolved Hide resolved
@richardotis richardotis requested a review from bocklund Nov 3, 2019
@bocklund

This comment has been minimized.

Copy link
Collaborator

bocklund commented Nov 4, 2019

I think the macOS failures are due to symengine 0.4 being installed on the latest commits (for some reason?). We probably just should update the deps to 'symengine>=0.5' and 'python-symengine>=0.5' everywhere

Copy link
Collaborator

bocklund left a comment

Approved pending passing tests

@richardotis richardotis merged commit 05ff0cf into develop Nov 4, 2019
5 checks passed
5 checks passed
continuous-integration/appveyor/branch AppVeyor build succeeded
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
continuous-integration/travis-ci/push The Travis CI build passed
Details
coverage/coveralls Coverage increased (+0.006%) to 86.05%
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
2 participants
You can’t perform that action at this time.