write_tdb function outputs incorrect tdb due to symengine

[wahab2604]

The write_tdb function produces an unreadable tdb due to how symengine evaluates EXP and LN.
symengine.sympify is used here and is where the evaluation takes places:
https://github.com/pycalphad/pycalphad/blob/ebcfbdb4dadfcce98a40db41c444a19842d8849e/pycalphad/io/tdb.py#L62

I have attached a two line script which opens a tdb then writes it out again, as you can see the input (which pycalphad reads and plot correctly ) is very different to the output tdb

According to Brandon the issue may be related to this:
symengine/symengine#1828

test.zip

[richardotis]

Does pycalphad produce the same phase diagram for both TDBs?

For example this input:

FUNCTION CAO0 1 (1 - EXP(-314.8272896/T)); 10000 N !
FUNCTION CAO1 1 (1 - EXP(-125.12812657/T)); 10000 N !
FUNCTION CAO2 1 (1 - EXP(-598.37301449/T)); 10000 N !
FUNCTION CAOL0 1 (1 - EXP(-629.27558/T)); 10000 N !
FUNCTION CAOL1 1 (1 - EXP(-234.0578/T)); 10000 N !

and this output:

FUNCTION CAO0 1.0 1.0-1.0 * 2.71828182845905**(-314.8272896 * T**(-1.0));
   10000.0 N !
FUNCTION CAO1 1.0 1.0-1.0 * 2.71828182845905**(-125.12812657 * T**(-1.0));
   10000.0 N !
FUNCTION CAO2 1.0 1.0-1.0 * 2.71828182845905**(-598.37301449 * T**(-1.0));
   10000.0 N !
FUNCTION CAOL0 1.0 1.0-1.0 * 2.71828182845905**(-629.27558 * T**(-1.0));
   10000.0 N !
FUNCTION CAOL1 1.0 1.0-1.0 * 2.71828182845905**(-234.0578 * T**(-1.0));
   10000.0 N !

These two should be identical for many significant figures, and I would expect pycalphad to produce identical phase diagrams for both. This is still a bug because most other Calphad implementations will not accept this output, but I want to understand if this is only a compatibility bug or if it's also a correctness bug in pycalphad.

[bocklund]

Related: symengine/symengine#1828 (comment), in particular:

We evaluate things like sin(2.0). Reason is that 2.0 is only accurate to 53 bits, so it doesn't make sense to evaluate sin(2.0) to more precision than the precision of 2.0. That's one reason to do it. The other is that it prevents explosion of the symbolic tree.
So, the policy that symengine uses is that if you send in floating point number to a function like (sin, exp, etc) you get back a floating point number.

[wahab2604]

The binplot in Pycalphad doesn’t produce a plot for the output tdb but does produce one for the input.

It’s specifically the EXP causing the issue, removing those functions will produce a TDB that will plot in pycalphad.

[richardotis]

We have logic in the TDB writer that is supposed to detect this case:

pycalphad/pycalphad/io/tdb.py

Lines 469 to 478 in ebcfbdb

	elif isinstance(expr, Pow):
	if expr.args[0] == E:
	# This is the exponential function
	terms = 'exp(' + self._stringify_expr(expr.args[1]) + ')'
	else:
	argument = self._stringify_expr(expr.args[0])
	if isinstance(expr.args[0], (Add, Mul)):
	argument = '( ' + argument + ' )'
	terms = argument + '**' + '(' + self._stringify_expr(expr.args[1]) + ')'
	return terms

We just need to work out why it isn't being triggered here.

[wahab2604]

I believe that part is working as intended, but from my small investigation its the reading of the TDB that might be causing the issue specifically the tdb_grammar() function:

pycalphad/pycalphad/io/tdb.py

Line 207 in ebcfbdb

func_expr.setParseAction(_make_piecewise_ast) + \

which then points to:

pycalphad/pycalphad/io/tdb.py

Lines 107 to 116 in ebcfbdb

	def _make_piecewise_ast(toks):
	"""
	Convenience function for converting tokens into a piecewise symengine object.
	"""
	cur_tok = 0
	expr_cond_pairs = []
	# Only one token: Not a piecewise function; just return the AST
	if len(toks) == 1:
	return _sympify_string(toks[0].strip(' ,'))

and then points to the sympify evaluation:

pycalphad/pycalphad/io/tdb.py

Line 62 in ebcfbdb

return sympify(expr_string).xreplace(v.supported_variables_in_databases).n()

The logic you referenced doesn't work in this case because the expr has already been evaluated upon creation of a Database object, so the 2.71** is written out according to the else:

pycalphad/pycalphad/io/tdb.py

Lines 473 to 478 in ebcfbdb

	else:
	argument = self._stringify_expr(expr.args[0])
	if isinstance(expr.args[0], (Add, Mul)):
	argument = '( ' + argument + ' )'
	terms = argument + '**' + '(' + self._stringify_expr(expr.args[1]) + ')'
	return terms

[richardotis]

It doesn't solve the general problem, but can we tweak the writing logic for Pow to treat a base number of ~2.71 as the exponential function, for writing purposes?