How to use the pycalphad.equilibrium function in pycalphad

error = 0
    if len(activity_datasets) == 0:
        return error

    for ds in activity_datasets:
        acr_component = ds['output'].split('_')[1]  # the component of interest
        # calculate the reference state equilibrium
        ref = ds['reference_state']
        # data_comps and data_phases ensures that we only do calculations on
        # the subsystem of the system defining the data.
        data_comps = ds['components']
        species = list(map(v.Species, data_comps))
        data_phases = filter_phases(dbf, species, candidate_phases=phases)
        ref_conditions = {_map_coord_to_variable(coord): val for coord, val in ref['conditions'].items()}
        ref_result = equilibrium(dbf, data_comps, ref['phases'], ref_conditions,
                                 model=phase_models, parameters=parameters,
                                 callables=callables)

        # calculate current chemical potentials
        # get the conditions
        conditions = {}
        # first make sure the conditions are paired
        # only get the compositions, P and T are special cased
        conds_list = [(cond, value) for cond, value in ds['conditions'].items() if cond not in ('P', 'T')]
        # ravel the conditions
        # we will ravel each composition individually, since they all must have the same shape
        for comp_name, comp_x in conds_list:
            P, T, X = ravel_conditions(ds['values'], ds['conditions']['P'], ds['conditions']['T'], comp_x)
            conditions[v.P] = P
            conditions[v.T] = T
            conditions[_map_coord_to_variable(comp_name)] = X

def time_equilibrium_al_ni(self):
        equilibrium(self.db_alfe, ['AL', 'FE', 'VA'], ['LIQUID', 'B2_BCC'], {v.X('AL'): 0.25, v.T: (300, 2000, 50), v.P: 101325}, pbar=False)

print("Entering region {}".format(start_pt))
        T_current = start_pt.temperature + delta
        x_current = start_pt.composition
        T_backtracks = 0; total_T_backtrack_factor = 1;
        converged = False
        while not converged:
            if (T_current < T_min) or (T_current > T_max):
                converged = True
                end_point = StartPoint(T_current, opposite_direction(curr_direction), compsets)
                start_points.add_end_point(end_point)
                if veryverbose:
                    print("Adding end point {}".format(end_point))
                continue
            curr_conds[v.T] = T_current
            curr_conds[x_cond] = x_current
            eq = equilibrium(dbf, comps, phases, curr_conds, **eq_kwargs)
            compsets = get_compsets(eq, indep_comp=indep_comp, indep_comp_index=comp_idx)
            if veryverbose:
                print("found compsets {} at T={}K X={:0.3f} eq_phases={}".format(compsets, T_current, x_current, eq.Phase.values.flatten()))
            if len(compsets) == 1:
                found_str = "Found single phase region {} at T={}K X={:0.3f}".format(compsets[0].phase_name, T_current, x_current)
                if T_backtracks < max_T_backtracks:
                    T_backtracks += 1
                    total_T_backtrack_factor *= T_backtrack_factor
                    T_backtrack = T_current - delta/total_T_backtrack_factor
                    if veryverbose:
                        print(found_str + " Backtracking in temperature from {}K to {}K ({}/{})".format(T_current, T_backtrack, T_backtracks, max_T_backtracks))
                    T_current = T_backtrack
                    continue
                elif not backtrack_raise:
                    converged = True
                    end_point = StartPoint(T_current, opposite_direction(curr_direction), compsets)

conds_list = [(cond, value) for cond, value in ds['conditions'].items() if cond not in ('P', 'T')]
        # ravel the conditions
        # we will ravel each composition individually, since they all must have the same shape
        for comp_name, comp_x in conds_list:
            P, T, X = ravel_conditions(ds['values'], ds['conditions']['P'], ds['conditions']['T'], comp_x)
            conditions[v.P] = P
            conditions[v.T] = T
            conditions[_map_coord_to_variable(comp_name)] = X
        # do the calculations
        # we cannot currently turn broadcasting off, so we have to do equilibrium one by one
        # invert the conditions dicts to make a list of condition dicts rather than a condition dict of lists
        # assume now that the ravelled conditions all have the same size
        conditions_list = [{c: conditions[c][i] for c in conditions.keys()} for i in range(len(conditions[v.T]))]
        current_chempots = []
        for conds in conditions_list:
            sample_eq_res = equilibrium(dbf, data_comps, data_phases, conds,
                                        model=phase_models, parameters=parameters,
                                        callables=callables)
            current_chempots.append(sample_eq_res.MU.sel(component=acr_component).values.flatten()[0])
        current_chempots = np.array(current_chempots)

        # calculate target chempots
        samples = np.array(ds['values']).flatten()
        target_chempots = target_chempots_from_activity(acr_component, samples, conditions[v.T], ref_result)
        # calculate the error
        weight = ds.get('weight', 1.0)
        pe = chempot_error(current_chempots, target_chempots, std_dev=std_dev/data_weight/weight)
        error += np.sum(pe)
        logging.debug('Activity error - data: {}, chemical potential difference: {}, probability: {}, reference: {}'.format(samples, current_chempots-target_chempots, pe, ds["reference"]))

    # TODO: write a test for this
    if np.any(np.isnan(np.array([error], dtype=np.float64))):  # must coerce sympy.core.numbers.Float to float64

-------
    A phase diagram as a figure.

    Examples
    --------
    None yet.
    """
    eq_kwargs = eq_kwargs if eq_kwargs is not None else dict()
    indep_comp = [key for key, value in conds.items() if isinstance(key, v.Composition) and len(np.atleast_1d(value)) > 1]
    indep_pot = [key for key, value in conds.items() if (type(key) is v.StateVariable) and len(np.atleast_1d(value)) > 1]
    if (len(indep_comp) != 1) or (len(indep_pot) != 1):
        raise ValueError('binplot() requires exactly one composition and one potential coordinate')
    indep_comp = indep_comp[0]
    indep_pot = indep_pot[0]

    full_eq = equilibrium(dbf, comps, phases, conds, **eq_kwargs)
    return eqplot(full_eq, x=indep_comp, y=indep_pot, **plot_kwargs)

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)'
        label_mapping['P'] = 'Pressure (Pa)'

for sp in initial_start_points:
                start_points.add_start_point(sp)

    # find a starting point
    starting_T = 0.9*(T_max - T_min)+T_min
    time_start = time.time()
    max_startpoint_discrepancy = np.max([tol_zero_one, tol_same, dx])
    while len(start_points.remaining_start_points) < 1:
        curr_conds[v.T] = starting_T
        hull = convex_hull(dbf, comps, phases, curr_conds, **eq_kwargs)
        cs = find_two_phase_region_compsets(hull, starting_T, indep_comp, comp_idx, discrepancy_tol=max_startpoint_discrepancy)
        if len(cs) == 2:
            # verify that these show up in the equilibrium calculation
            specific_conds = deepcopy(curr_conds)
            specific_conds[x_cond] = BinaryCompSet.mean_composition(cs)
            eq_cs = get_compsets(equilibrium(dbf, comps, phases, specific_conds, **eq_kwargs), indep_comp=indep_comp, indep_comp_index=comp_idx)
            if len(eq_cs) == 2:
                # add a direction of dT > 0 and dT < 0
                zpf_boundaries.add_compsets(eq_cs)
                # shift starting_T so they start at the same place.
                start_points.add_start_point(StartPoint(starting_T - dT, pos, eq_cs))
                start_points.add_start_point(StartPoint(starting_T + dT, neg, eq_cs))

        if starting_T - dT > T_min:
            starting_T -= dT
        else:
            raise StartingPointError("Unable to find an initial starting point.")
    if verbose:
        print("Found start points {} in {:0.2f}s".format(start_points, time.time()-time_start))

    # Main loop
    while len(start_points.remaining_start_points) > 0:

Keyword arguments to eqplot().

    Returns
    -------
    A phase diagram as a figure.

    Examples
    --------
    None yet.
    """
    eq_kwargs = eq_kwargs if eq_kwargs is not None else dict()
    indep_comps = [key for key, value in conds.items() if isinstance(key, v.Composition) and len(np.atleast_1d(value)) > 1]
    indep_pots = [key for key, value in conds.items() if (type(key) is v.StateVariable) and len(np.atleast_1d(value)) > 1]
    if (len(indep_comps) != 2) or (len(indep_pots) != 0):
        raise ValueError('ternplot() requires exactly two composition coordinates')
    full_eq = equilibrium(dbf, comps, phases, conds, **eq_kwargs)
    # TODO: handle x and y as strings with #87
    x = x if x in indep_comps else indep_comps[0]
    y = y if y in indep_comps else indep_comps[1]
    return eqplot(full_eq, x=x, y=y, **plot_kwargs)

--------

    >>> from pycalphad import Database, variables as v  # doctest: +SKIP
    >>> from pycalphad.plot.eqplot import eqplot  # doctest: +SKIP
    >>> from espei.datasets import load_datasets, recursive_glob  # doctest: +SKIP
    >>> datasets = load_datasets(recursive_glob('.', '*.json'))  # doctest: +SKIP
    >>> dbf = Database('my_databases.tdb')  # doctest: +SKIP
    >>> my_phases = list(dbf.phases.keys())  # doctest: +SKIP
    >>> multiplot(dbf, ['CU', 'MG', 'VA'], my_phases, {v.P: 101325, v.T: 1000, v.X('MG'): (0, 1, 0.01)}, datasets)  # doctest: +SKIP

    """
    eq_kwargs = eq_kwargs or dict()
    plot_kwargs = plot_kwargs or dict()
    data_kwargs = data_kwargs or dict()

    eq_result = equilibrium(dbf, comps, phases, conds, **eq_kwargs)
    ax = eqplot(eq_result, **plot_kwargs)
    ax = eqdataplot(eq_result, datasets, ax=ax, plot_kwargs=data_kwargs)
    return ax

def time_equilibrium_al_fe(self):
        equilibrium(self.db_alni, ['AL', 'NI', 'VA'], ['LIQUID', 'FCC_L12'], {v.X('AL'): 0.10, v.T: (300, 2500, 20), v.P: 101325}, pbar=False)