Build tool to help generate symbolic expressions for unary plugins

[bocklund]

Unary plugins (see the unary skeleton repo) need to be written with the energy functions as valid SymPy expressions using pycalphad StateVariables. Many users have unary Gibbs energies and pure element lattice stabilities implemented as functions in TDB files.

A tool should exist to convert these TDB files valid Python/SymPy expressions so that it's easier to create reference state plugins for ESPEI. Here's some example code based on an old notebook that I wrote

from pycalphad import Database
dbf = Database('Mg-piecewise-example.tdb')
from sympy import Symbol
from sympy.utilities import topological_sort
from pycalphad import variables as v

def extract_symbol_value(symbol):
    """Return the value of a Piecewise if it only has one expression
    """
    # A Piecewise with a length of 2 has the expression over the range and then 0 otherwise.
    # We want to return the expression in this case, otherwise the expression is a "real" Piecewise.
    if len(symbol.args) == 2:
        return symbol.args[0].expr
    return symbol

def sorted_symbols(symbols):
    """Takes a dictionary of {Symbol: Piecewise} and returns the topological sort replacement
    """
    # Topological sorting requires the verticies and edges of the a graph of symbol dependencies.
    verticies = list(symbols.keys())
    edges = []
    for vertex in verticies:
        free_symbols = symbols[vertex].free_symbols
        # Filter out the StateVariables.
        # This will ensure P, T variables are not included in the graph because they can be variables in TDBs.
        free_symbols = [symbol for symbol in free_symbols if not isinstance(symbol, v.StateVariable)]
        # construct the edges (dependencies) of the symbols if there are any
        if len(free_symbols):
            edges.extend([(str(vertex), str(symbol)) for symbol in free_symbols])
    return topological_sort((verticies, edges))

def get_reduced_expression(dbf, symbol):
    """Return an expression for the passed symbol (as a string) reduced to only a function of state variables
    """
    # To reduce an expression, we have to pass a list of substitions to the expression in the order
    # we want to substitute them to get to having only StateVaribles as degrees of freedom.
    # Topological sorting orders the replacements in the correct order.
    topologically_sorted_symbols = sorted_symbols(dbf.symbols)
    # Then we construct the values that we will replace each symbol with, deconstructing the
    # symbols that are Piecewise, if necessary.
    replacements = [(sym, extract_symbol_value(dbf.symbols[sym])) for sym in topologically_sorted_symbols]
    return dbf.symbols[symbol].subs(replacements)

This can be run via

print(get_reduced_expression(dbf, 'GHSERMg'))

Which for, example, may print:

Piecewise((0.0047195465*T**2 + 3.2404446864*T*log(exp(95.82831/T) - 1.0) + 21.6535319868*T*log(exp(247.8675/T) - 1.0) - 10652.1012810789, And(T < 36.71926, T >= 0.01)), (-1.53643054262276e-5*T**3 + 0.00810454205399037*T**2 - 0.124294531845816*T*log(T) + 3.2404446864*T*log(exp(95.82831/T) - 1.0) + 21.6535319868*T*log(exp(247.8675/T) - 1.0) + 0.385723396310737*T - 10653.6226154894, And(T < 143.18844, T >= 36.71926)), (-0.0050954035*T**2 + 1.765785080115*T*log(T) + 3.2404446864*T*log(exp(95.82831/T) - 1.0) + 21.6535319868*T*log(exp(247.8675/T) - 1.0) - 8.0518972866125*T - 10563.4100984519, And(T < 922.205302616508, T >= 143.18844)), (-34.3088*T*log(T) + 204.341485347575*T - 13775.4156328263 + 9.4687256586798e+27/T**9, And(T < 3000.0, T >= 922.205302616508)), (0, True))

This output could be directly pasted into a module starting with:

from sympy import *
from pycalphad import *

# <paste output here>