How to use the pycalphad.variables.T function in pycalphad

We just have a template database that has the phase defined. We then hot
    patch the Model object to have the GM from the fixed model we printed out
    and the data we printed out. The hot patch is needed because this is
    formation enthalpy data and the model needs to have the lower order terms
    in composition.

    One possible issue is that the new GM in the fixed model does not have any
    individual contributions, so it cannot be used to test excluded model
    contributions. The only excluded model contributions in this data are
    idmix, but the property we are testing is HM_FORM, so the feature transform
    of the idmix property should be zero.

    """
    # Hack the namespace to make the copy-pasted Gibbs energy function work
    from sympy import log, Piecewise
    T = v.T

    data = [{'components': ['AL', 'NI', 'VA'], 'phases': ['BCC_B2'], 'solver': {'mode': 'manual', 'sublattice_occupancies': [[1.0, [0.5, 0.5], 1.0], [1.0, [0.75, 0.25], 1.0]], 'sublattice_site_ratios': [0.5, 0.5, 1.0], 'sublattice_configurations': (('AL', ('NI', 'VA'), 'VA'), ('AL', ('NI', 'VA'), 'VA')), 'comment': 'BCC_B2 sublattice configuration (2SL)'}, 'conditions': {'P': 101325.0, 'T': np.array([300.])}, 'reference_state': 'SGTE91', 'output': 'HM_FORM', 'values': np.array([[[-40316.61077, -56361.58554]]]), 'reference': 'C. Jiang 2009 (constrained SQS)', 'excluded_model_contributions': ['idmix']}, {'components': ['AL', 'NI', 'VA'], 'phases': ['BCC_B2'], 'solver': {'mode': 'manual', 'sublattice_occupancies': [[1.0, [0.5, 0.5], 1.0], [1.0, [0.75, 0.25], 1.0]], 'sublattice_site_ratios': [0.5, 0.5, 1.0], 'sublattice_configurations': (('AL', ('NI', 'VA'), 'VA'), ('AL', ('NI', 'VA'), 'VA')), 'comment': 'BCC_B2 sublattice configuration (2SL)'}, 'conditions': {'P': 101325.0, 'T': np.array([300.])}, 'reference_state': 'SGTE91', 'output': 'HM_FORM', 'values': np.array([[[-41921.43363, -57769.49473]]]), 'reference': 'C. Jiang 2009 (relaxed SQS)', 'excluded_model_contributions': ['idmix']}]
    NEW_GM = 8.3145*T*(0.5*Piecewise((v.SiteFraction("BCC_B2", 0, "AL")*log(v.SiteFraction("BCC_B2", 0, "AL")), v.SiteFraction("BCC_B2", 0, "AL") > 1.0e-16), (0, True))/(0.5*v.SiteFraction("BCC_B2", 0, "AL") + 0.5*v.SiteFraction("BCC_B2", 0, "NI") + 0.5*v.SiteFraction("BCC_B2", 1, "AL") + 0.5*v.SiteFraction("BCC_B2", 1, "NI")) + 0.5*Piecewise((v.SiteFraction("BCC_B2", 0, "NI")*log(v.SiteFraction("BCC_B2", 0, "NI")), v.SiteFraction("BCC_B2", 0, "NI") > 1.0e-16), (0, True))/(0.5*v.SiteFraction("BCC_B2", 0, "AL") + 0.5*v.SiteFraction("BCC_B2", 0, "NI") + 0.5*v.SiteFraction("BCC_B2", 1, "AL") + 0.5*v.SiteFraction("BCC_B2", 1, "NI")) + 0.5*Piecewise((v.SiteFraction("BCC_B2", 0, "VA")*log(v.SiteFraction("BCC_B2", 0, "VA")), v.SiteFraction("BCC_B2", 0, "VA") > 1.0e-16), (0, True))/(0.5*v.SiteFraction("BCC_B2", 0, "AL") + 0.5*v.SiteFraction("BCC_B2", 0, "NI") + 0.5*v.SiteFraction("BCC_B2", 1, "AL") + 0.5*v.SiteFraction("BCC_B2", 1, "NI")) + 0.5*Piecewise((v.SiteFraction("BCC_B2", 1, "AL")*log(v.SiteFraction("BCC_B2", 1, "AL")), v.SiteFraction("BCC_B2", 1, "AL") > 1.0e-16), (0, True))/(0.5*v.SiteFraction("BCC_B2", 0, "AL") + 0.5*v.SiteFraction("BCC_B2", 0, "NI") + 0.5*v.SiteFraction("BCC_B2", 1, "AL") + 0.5*v.SiteFraction("BCC_B2", 1, "NI")) + 0.5*Piecewise((v.SiteFraction("BCC_B2", 1, "NI")*log(v.SiteFraction("BCC_B2", 1, "NI")), v.SiteFraction("BCC_B2", 1, "NI") > 1.0e-16), (0, True))/(0.5*v.SiteFraction("BCC_B2", 0, "AL") + 0.5*v.SiteFraction("BCC_B2", 0, "NI") + 0.5*v.SiteFraction("BCC_B2", 1, "AL") + 0.5*v.SiteFraction("BCC_B2", 1, "NI")) + 0.5*Piecewise((v.SiteFraction("BCC_B2", 1, "VA")*log(v.SiteFraction("BCC_B2", 1, "VA")), v.SiteFraction("BCC_B2", 1, "VA") > 1.0e-16), (0, True))/(0.5*v.SiteFraction("BCC_B2", 0, "AL") + 0.5*v.SiteFraction("BCC_B2", 0, "NI") + 0.5*v.SiteFraction("BCC_B2", 1, "AL") + 0.5*v.SiteFraction("BCC_B2", 1, "NI")) + Piecewise((v.SiteFraction("BCC_B2", 2, "VA")*log(v.SiteFraction("BCC_B2", 2, "VA")), v.SiteFraction("BCC_B2", 2, "VA") > 1.0e-16), (0, True))/(0.5*v.SiteFraction("BCC_B2", 0, "AL") + 0.5*v.SiteFraction("BCC_B2", 0, "NI") + 0.5*v.SiteFraction("BCC_B2", 1, "AL") + 0.5*v.SiteFraction("BCC_B2", 1, "NI"))) + (45262.9*v.SiteFraction("BCC_B2", 0, "AL")*v.SiteFraction("BCC_B2", 0, "NI")*v.SiteFraction("BCC_B2", 1, "AL")*v.SiteFraction("BCC_B2", 2, "VA") + 45262.9*v.SiteFraction("BCC_B2", 0, "AL")*v.SiteFraction("BCC_B2", 1, "AL")*v.SiteFraction("BCC_B2", 1, "NI")*v.SiteFraction("BCC_B2", 2, "VA"))/(0.5*v.SiteFraction("BCC_B2", 0, "AL") + 0.5*v.SiteFraction("BCC_B2", 0, "NI") + 0.5*v.SiteFraction("BCC_B2", 1, "AL") + 0.5*v.SiteFraction("BCC_B2", 1, "NI")) + (1.0*v.SiteFraction("BCC_B2", 0, "AL")*v.SiteFraction("BCC_B2", 1, "AL")*v.SiteFraction("BCC_B2", 2, "VA")*Piecewise((10083 - 4.813*T, (T >= 298.15) & (T < 2900.0)), (0, True)) + v.SiteFraction("BCC_B2", 0, "AL")*v.SiteFraction("BCC_B2", 1, "NI")*v.SiteFraction("BCC_B2", 2, "VA")*(9.52839e-8*T**3 + 0.00123463*T**2 + 0.000871898*T*log(T) + 1.31471*T - 64435.3 + 23095.2/T) + v.SiteFraction("BCC_B2", 0, "AL")*v.SiteFraction("BCC_B2", 1, "VA")*v.SiteFraction("BCC_B2", 2, "VA")*(10.0*T + 16432.5) + v.SiteFraction("BCC_B2", 0, "NI")*v.SiteFraction("BCC_B2", 1, "AL")*v.SiteFraction("BCC_B2", 2, "VA")*(9.52839e-8*T**3 + 0.00123463*T**2 + 0.000871898*T*log(T) + 1.31471*T - 64435.3 + 23095.2/T) + 1.0*v.SiteFraction("BCC_B2", 0, "NI")*v.SiteFraction("BCC_B2", 1, "NI")*v.SiteFraction("BCC_B2", 2, "VA")*Piecewise((8715.084 - 3.556*T, (T >= 298.15) & (T < 3000.0)), (0, True)) + 32790.6*v.SiteFraction("BCC_B2", 0, "NI")*v.SiteFraction("BCC_B2", 1, "VA")*v.SiteFraction("BCC_B2", 2, "VA") + v.SiteFraction("BCC_B2", 0, "VA")*v.SiteFraction("BCC_B2", 1, "AL")*v.SiteFraction("BCC_B2", 2, "VA")*(10.0*T + 16432.5) + 32790.6*v.SiteFraction("BCC_B2", 0, "VA")*v.SiteFraction("BCC_B2", 1, "NI")*v.SiteFraction("BCC_B2", 2, "VA"))/(0.5*v.SiteFraction("BCC_B2", 0, "AL") + 0.5*v.SiteFraction("BCC_B2", 0, "NI") + 0.5*v.SiteFraction("BCC_B2", 1, "AL") + 0.5*v.SiteFraction("BCC_B2", 1, "NI"))

    dbf = Database("""$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    $ Date: 2019-12-08 18:05
    $ Components: AL, NI, VA
    $ Phases: BCC_B2
    $ Generated by brandon (pycalphad 0.8.1.post1)
    $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

    ELEMENT AL FCC_A1 26.982 4577.3 28.322 !
    ELEMENT NI FCC_A1 58.69 4787.0 29.796 !
    ELEMENT VA VACUUM 0.0 0.0 0.0 !

    TYPE_DEFINITION % SEQ * !

def build_piecewise_matrix(sympy_obj, cur_exponents, low_temp, high_temp, output_matrix, symbol_matrix, param_symbols):
    if sympy_obj == v.T:
        sympy_obj = Mul(1.0, v.T)
    elif sympy_obj == v.P:
        sympy_obj = Mul(1.0, v.P)
    if isinstance(sympy_obj, (Float, Integer, Rational)) or \
            (isinstance(sympy_obj, (log, exp)) and isinstance(sympy_obj.args[0], (Float, Integer, Rational))):
        result = float(sympy_obj)
        if result != 0.0:
            output_matrix.append([low_temp, high_temp] + list(cur_exponents) + [result])
    elif isinstance(sympy_obj, Piecewise):
        piece_args = [i for i in sympy_obj.args if i.expr != S.Zero]
        intervals = [to_interval(i.cond) for i in piece_args]
        if (len(intervals) > 1) and Intersection(intervals) != EmptySet():
            raise ValueError('Overlapping intervals cannot be represented: {}'.format(intervals))
        if not isinstance(Union(intervals), Interval):
            raise ValueError('Piecewise intervals must be continuous')
        if not all([arg.cond.free_symbols == {v.T} for arg in piece_args]):

                      "HM_MIX": lambda GM: GM - v.T*sympy.diff(GM, v.T),
                      "HM": lambda GM: GM - v.T*sympy.diff(GM, v.T)}

'massfuncs': {},
        'massgradfuncs': {},
        'masshessfuncs': {},
        'callables': {},
        'grad_callables': {},
        'hess_callables': {},
        'internal_cons': {},
        'internal_jac': {},
        'internal_cons_hess': {},
        'mp_cons': {},
        'mp_jac': {},
    }

    state_variables = get_state_variables(models=models)
    state_variables |= additional_statevars
    if state_variables != {v.T, v.P, v.N}:
        warnings.warn("State variables in `build_callables` are not {{N, P, T}}, "
                      "but {}. Be sure you know what you are doing. "
                      "State variables can be added with the `additional_statevars` "
                      "argument.".format(state_variables))
    state_variables = sorted(state_variables, key=str)

    for name in phases:
        mod = models[name]
        site_fracs = mod.site_fractions
        try:
            out = getattr(mod, output)
        except AttributeError:
            raise AttributeError('Missing Model attribute {0} specified for {1}'
                                 .format(output, mod.__class__))

        # Build the callables of the output

for pth in pths:
        if scatter:
            ax.scatter(pth[0], pth[1], c=pth[2], edgecolor='None', s=3, zorder=2)
        else:
            ax.plot(pth[0], pth[1], c=pth[2], markersize=3, zorder=2)
        if tielines:
            ax.add_collection(tieline_coll)
    box = ax.get_position()
    ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
    ax.legend(handles=legend_handles, loc='center left', bbox_to_anchor=(1, 0.5))
    ax.tick_params(axis='both', which='major', labelsize=14)
    ax.grid(True)
    plot_title = '-'.join([component for component in sorted(zpf_boundary_sets.components) if component != 'VA'])
    ax.set_title(plot_title, fontsize=20)
    ax.set_xlabel(_axis_label(zpf_boundary_sets.indep_comp), labelpad=15, fontsize=20)
    ax.set_ylabel(_axis_label(v.T), fontsize=20)
    ax.set_xlim(0, 1)

conditions = dict()
        if 'T' in row:
            conditions[v.T] = row['T']
        if 'P' in row:
            conditions[v.P] = row['P']
        for comp in global_comps[:-1]:
            conditions[v.X(comp)] = row['X('+comp+')']
        statevars = dict((str(key), value) for key, value in conditions.items() if key in [v.T, v.P])
        eq = equilibrium(dbf, comps, row['Phase'], conditions,
                         model=mods, callables=iter_callables,
                         grad_callables=iter_grad_callables, verbose=False)
        # TODO: Support for miscibility gaps, i.e., FCC_A1#1 specification
        eq_values = eq['Y'].sel(vertex=0).values
        #print(eq_values)
        # TODO: All the needed 'Types' should be precalculated and looked up
        variables = sorted(mods[row['Phase']].energy.atoms(v.StateVariable).union({v.T, v.P}), key=str)
        output_callables = {row['Phase']: functools.partial(make_callable(getattr(mods[row['Phase']],
                                                                                  row['Type']),
                                                            itertools.chain(param_names, variables)),
                                                                            *parvalues)}
        calculated_value = calculate(dbf, comps, row['Phase'], output=row['Type'],
                                     model=mods, callables=output_callables,
                                     points=eq_values, **statevars)
        res[idx] = float(row['Value']) - float(calculated_value[row['Type']].values)
        #print('res', idx, res[idx])
    return res

def write_parameter(param_to_write):
        constituents = ':'.join([','.join(sorted([i.upper() for i in subl]))
                         for subl in param_to_write.constituent_array])
        # TODO: Handle references
        paramx = param_to_write.parameter
        if not isinstance(paramx, Piecewise):
            # Non-piecewise parameters need to be wrapped to print correctly
            # Otherwise TC's TDB parser will fail
            paramx = Piecewise((paramx, And(v.T >= 1, v.T < 10000)))
        exprx = TCPrinter().doprint(paramx).upper()
        if ';' not in exprx:
            exprx += '; N'
        if param_to_write.diffusing_species is not None:
            ds = "&" + param_to_write.diffusing_species
        else:
            ds = ""
        return "PARAMETER {}({}{},{};{}) {} !\n".format(param_to_write.parameter_type.upper(),
                                                        param_to_write.phase_name.upper(),
                                                        ds,
                                                        constituents,
                                                        param_to_write.parameter_order,
                                                        exprx)
    if groupby == 'subsystem':

    enthalpy = HM = property(lambda self: self.GM - v.T*self.GM.diff(v.T))
    heat_capacity = CPM = property(lambda self: -v.T*self.GM.diff(v.T, v.T))

String name of the property, e.g. 'HM_MIX'
    candidate_models : list
        List of SymPy parameters that can be fit for this property.
    desired_data : dict
        Full dataset dictionary containing values, conditions, etc.

    Returns
    -------
    numpy.ndarray
        An MxN matrix of M samples (from desired data) and N features.

    """
    transformed_features = sympy.Matrix([feature_transforms[prop](i) for i in candidate_models])
    all_samples = get_samples(desired_data)
    feature_matrix = np.empty((len(all_samples), len(transformed_features)), dtype=np.float)
    feature_matrix[:, :] = [transformed_features.subs({v.T: temp, 'YS': ys, 'V_I': v_i, 'V_J': v_j, 'V_K': v_k}).evalf() for temp, (ys, (v_i, v_j, v_k)) in all_samples]
    return feature_matrix

# To work around differences in sublattice models, relax the internal dof
            global_comps = sorted(set(data_group[0].json['components']) - set(['VA']))
            compositions = \
                _map_internal_dof(input_database,
                                  sorted(data_group[0].json['components']),
                                  data_group[0].json['phases'][0],
                                  np.atleast_2d(
                                      data_group[0].json['solver']['sublattice_configuration']).astype(np.float))
            # Tiny perturbation to work around a bug in lower_convex_hull (gh-28)
            compare_conds = {v.X(comp): np.add(compositions[:, idx], 1e-4).flatten().tolist()
                             for idx, comp in enumerate(global_comps[:-1])}
            compare_conds.update({v.__dict__[key]: value for key, value in conds.items()})
            # We only want to relax the internal dof at the lowest temperature
            # This will help us capture the most related sublattice config since solver mode=manual
            # probably means this is first-principles data
            compare_conds[v.T] = 300.
            eqres = equilibrium(dbf, data_group[0].json['components'],
                                str(data_group[0].json['phases'][0]), compare_conds, verbose=False)
            internal_dof = sum(map(len, dbf.phases[data_group[0].json['phases'][0]].constituents))
            largest_phase_fraction = eqres['NP'].values.argmax()
            eqpoints = eqres['Y'].values[..., largest_phase_fraction, :internal_dof]
            result = calculate(dbf, data_group[0].json['components'],
                               str(data_group[0].json['phases'][0]), output=y,
                               points=eqpoints, **conds)
            # Don't show CALPHAD results below 300 K because they're meaningless right now
            result = result.sel(T=slice(300, None))
            figure.gca().plot(result[x].values.flatten(), result[y].values.flatten(), label=label)
        label_mapping = dict(x=x, y=y)
        label_mapping['CPM'] = 'Molar Heat Capacity (J/mol-atom-K)'
        label_mapping['SM'] = 'Molar Entropy (J/mol-atom-K)'
        label_mapping['HM'] = 'Molar Enthalpy (J/mol-atom)'
        label_mapping['T'] = 'Temperature (K)'