How to use the pycalphad.calculate function in pycalphad

def run_test():
    dbf = Database()
    dbf.elements = frozenset(['A'])
    dbf.add_phase('TEST', {}, [1])
    dbf.add_phase_constituents('TEST', [['A']])
    # add THETA parameters here
    dbf.add_parameter('THETA', 'TEST', [['A']], 0, 334.)
    conds = {v.T: np.arange(1.,800.,1), v.P: 101325}
    res = calculate(dbf, ['A'], 'TEST', T=conds[v.T], P=conds[v.P],
                    model=EinsteinModel, output='testprop')
    #res_TE = calculate(dbf, ['A'], 'TEST', T=conds[v.T], P=conds[v.P],
    #                model=EinsteinModel, output='einstein_temperature')
    import matplotlib.pyplot as plt
    plt.scatter(res['T'], res['testprop'])
    plt.xlabel('Temperature (K)')
    plt.ylabel('Molar Heat Capacity (J/mol-K)')
    plt.savefig('einstein.png')
    print(dbf.to_string(fmt='tdb'))

result = LightDataset({output: (prop_dims, np.full(prop_shape, np.nan))}, coords=coord_dict)
    # For each phase select all conditions where that phase exists
    # Perform the appropriate calculation and then write the result back
    for phase in active_phases:
        dof = sum([len(x) for x in dbf.phases[phase].constituents])
        current_phase_indices = (data.Phase == phase)
        if ~np.any(current_phase_indices):
            continue
        points = data.Y[np.nonzero(current_phase_indices)][..., :dof]
        statevar_indices = np.nonzero(current_phase_indices)[:len(indep_vals)]
        statevars = {key: np.take(np.asarray(vals), idx)
                     for key, vals, idx in zip(indep_vars, indep_vals, statevar_indices)}
        statevars.update(kwargs)
        if statevars.get('mode', None) is None:
            statevars['mode'] = 'numpy'
        calcres = calculate(dbf, comps, [phase], output=output, points=points, broadcast=False,
                            callables=callables, parameters=parameters, **statevars)
        result[output][np.nonzero(current_phase_indices)] = calcres[output].values
    if not per_phase:
        out = np.nansum(result[output] * data['NP'], axis=-1)
        dv_output = result.data_vars[output]
        result.remove(output)
        # remove the vertex coordinate because we summed over it
        result.add_variable(output, dv_output[0][:-1], out)
    else:
        dv_phase = data.data_vars['Phase']
        dv_np = data.data_vars['NP']
        result.add_variable('Phase', dv_phase[0], dv_phase[1])
        result.add_variable('NP', dv_np[0], dv_np[1])
    return result

ax.plot(temperatures, response_data,
                           label='This work', color='k', **extra_kwargs)
            ax.set_xlabel(plot_mapping.get(x, x))
            ax.set_ylabel(plot_mapping.get(y, y))
        else:
            bar_labels.append('This work')
            bar_data.append(response_data[0])
    elif len(endpoints) == 2:
        # Binary interaction parameter
        first_endpoint = _translate_endmember_to_array(endpoints[0], mod.ast.atoms(v.SiteFraction))
        second_endpoint = _translate_endmember_to_array(endpoints[1], mod.ast.atoms(v.SiteFraction))
        point_matrix = np.linspace(0, 1, num=100)[None].T * second_endpoint + \
            (1 - np.linspace(0, 1, num=100))[None].T * first_endpoint
        # TODO: Real temperature support
        point_matrix = point_matrix[None, None]
        predicted_quantities = calculate(dbf, comps, [phase_name], output=yattr,
                                         T=300, P=101325, points=point_matrix, model=mod, mode='numpy')
        response_data = predicted_quantities[yattr].values.flatten()
        if not bar_chart:
            extra_kwargs = {}
            if len(response_data) < 10:
                extra_kwargs['markersize'] = 20
                extra_kwargs['marker'] = '.'
                extra_kwargs['linestyle'] = 'none'
                extra_kwargs['clip_on'] = False
            ax.plot(np.linspace(0, 1, num=100), response_data, label='This work', color='k', **extra_kwargs)
            ax.set_xlim((0, 1))
            ax.set_xlabel(str(':'.join(endpoints[0])) + ' to ' + str(':'.join(endpoints[1])))
            ax.set_ylabel(plot_mapping.get(y, y))
        else:
            bar_labels.append('This work')
            bar_data.append(response_data[0])

# Remove derivatives wrt chemical potentials
                cast_grad = cast_grad[..., properties.MU.shape[-1]:]
                grad = grad - cast_grad
                plane_args = itertools.chain([properties.MU.values[..., i][..., None] for i in range(properties.MU.shape[-1])],
                                             [points[..., i] for i in range(points.shape[-1])])
                cast_hess = np.array(plane_hess(*plane_args), dtype=np.float)
                # Remove derivatives wrt chemical potentials
                cast_hess = cast_hess[..., properties.MU.shape[-1]:, properties.MU.shape[-1]:]
                cast_hess = -cast_hess + hess
                hess = cast_hess.astype(np.float, copy=False)
                try:
                    e_matrix = np.linalg.inv(hess)
                except np.linalg.LinAlgError:
                    print(hess)
                    raise
                current = calculate(dbf, comps, name, output='GM',
                                    model=models, callables=callable_dict,
                                    fake_points=False,
                                    points=points.reshape(points.shape[:len(indep_vals)] + (-1, points.shape[-1])),
                                    **grid_opts)
                current_plane = np.multiply(current.X.values.reshape(points.shape[:-1] + (len(components),)),
                                            properties.MU.values[..., np.newaxis, :]).sum(axis=-1)
                current_df = current.GM.values.reshape(points.shape[:-1]) - current_plane
                #print('Inv hess check: ', np.isnan(e_matrix).any())
                #print('grad check: ', np.isnan(grad).any())
                dy_unconstrained = -np.einsum('...ij,...j->...i', e_matrix, grad)
                #print('dy_unconstrained check: ', np.isnan(dy_unconstrained).any())
                proj_matrix = np.dot(e_matrix, constraint_jac.T)
                inv_matrix = np.rollaxis(np.dot(constraint_jac, proj_matrix), 0, -1)
                inv_term = np.linalg.inv(inv_matrix)
                #print('inv_term check: ', np.isnan(inv_term).any())
                first_term = np.einsum('...ij,...jk->...ik', proj_matrix, inv_term)

# 'calculate' accepts conditions through its keyword arguments
    grid_opts = calc_opts.copy()
    grid_opts.update({key: value for key, value in str_conds.items() if key in indep_vars})
    if 'pdens' not in grid_opts:
        grid_opts['pdens'] = 100

    coord_dict = str_conds.copy()
    coord_dict['vertex'] = np.arange(len(components))
    grid_shape = np.meshgrid(*coord_dict.values(),
                             indexing='ij', sparse=False)[0].shape
    coord_dict['component'] = components
    if verbose:
        print('Computing initial grid', end=' ')
    # TODO: vectorize this entire calculation over the conditions
    # TODO: Every condition-set should have its own grid
    grid = calculate(dbf, comps, active_phases, output='GM',
                     model=models, callables=callable_dict, fake_points=True, **grid_opts)

    if verbose:
        print('[{0} points, {1}]'.format(len(grid.points), sizeof_fmt(grid.nbytes)), end='\n')

    properties = Dataset({'NP': (list(str_conds.keys()) + ['vertex'],
                                 np.empty(grid_shape)),
                          'GM': (list(str_conds.keys()),
                                 np.empty(grid_shape[:-1])),
                          'MU': (list(str_conds.keys()) + ['component'],
                                 np.empty(grid_shape)),
                          'points': (list(str_conds.keys()) + ['vertex'],
                                     np.empty(grid_shape, dtype=np.int))
                          },
                          coords=coord_dict,
                          attrs={'hull_iterations': 1, 'solve_iterations': 0,

#print('Biggest step:', biggest_step)
                    #print('points', points)
                    #print('grad of points', grad)
                    #print('new_direction', new_direction)
                    #print('alpha', alpha)
                    #print('new_points', new_points)
                    pass
                points = new_points
                newton_iteration += 1
            new_points = points.reshape(points.shape[:len(indep_vals)] + (-1, points.shape[-1]))
            new_points = np.concatenate((current_site_fractions[..., :dof], new_points), axis=-2)
            points_dict[name] = new_points

        if verbose:
            print('Rebuilding grid', end=' ')
        grid = calculate(dbf, comps, active_phases, output='GM',
                         model=models, callables=callable_dict,
                         fake_points=True, points=points_dict, **grid_opts)
        if verbose:
            print('[{0} points, {1}]'.format(len(grid.points), sizeof_fmt(grid.nbytes)), end='\n')
        properties.attrs['hull_iterations'] += 1
    # Make sure all the verbose output appears to the user
    if verbose:
        print('Refining equilibrium')
        sys.stdout.flush()
    # One last call to ensure 'properties' and 'grid' are consistent with one another
    lower_convex_hull(grid, properties, verbose=verbose)
    indexer = []
    for idx, vals in enumerate(indep_vals):
        indexer.append(np.arange(len(vals), dtype=np.int)[idx * (np.newaxis,) + np.index_exp[:] + \
                                                          (len(conds.keys())-idx+1) * (np.newaxis,)])
    indexer.append(properties['points'].values[..., np.newaxis])

def compute_values(*args):
        prefill_callables = {key: functools.partial(*itertools.chain([func], args[:len(params)]))
                             for key, func in callables.items()}
        result = calculate(dbf, data['components'], data['phases'], output=data['output'],
                           points=np.atleast_2d(data['solver']['sublattice_configuration']).astype(np.float),
                           callables=prefill_callables, model=fit_models, **extra_conds)
        return result

mod_latticeonly.models = {key: value for key, value in mod_latticeonly.models.items()
                                          if key == 'ref'}
                temps = data['conditions'].get('T', 300)
                pressures = data['conditions'].get('P', 101325)
                points = build_sitefractions(phase_name, data['solver']['sublattice_configurations'],
                                             data['solver']['sublattice_occupancies'])
                for point_idx in range(len(points)):
                    missing_variables = mod_latticeonly.ast.atoms(v.SiteFraction) - set(points[point_idx].keys())
                    # Set unoccupied values to zero
                    points[point_idx].update({key: 0 for key in missing_variables})
                    # Change entry to a sorted array of site fractions
                    points[point_idx] = list(OrderedDict(sorted(points[point_idx].items(), key=str)).values())
                points = np.array(points, dtype=np.float)
                # TODO: Real temperature support
                points = points[None, None]
                stability = calculate(dbf, comps, [phase_name], output=data['output'][:-5],
                                      T=temps, P=pressures, points=points,
                                      model=mod_latticeonly, mode='numpy')
                response_data -= stability[data['output'][:-5]].values.squeeze()

        response_data += np.array(data['values'], dtype=np.float)
        response_data = response_data.flatten()
        if not bar_chart:
            extra_kwargs = {}
            extra_kwargs['markersize'] = 8
            extra_kwargs['linestyle'] = 'none'
            extra_kwargs['clip_on'] = False
            ref = data.get('reference', '')
            mark = symbol_map[ref]['markers']
            ax.plot(indep_var_data, response_data,
                    label=symbol_map[ref]['formatted'],
                    marker=mark['marker'],