import catmap
from catmap import ReactionModelWrapper
from catmap.model import ReactionModel
from catmap.functions import smooth_piecewise_linear
from catmap.functions import parse_constraint
import pylab as plt
import numpy as np

try:
    from scipy.optimize import fmin
    from scipy import integrate
except ImportError:
    fmin = None
    integrate = None

from itertools import combinations_with_replacement as combinations

class FirstOrderInteractions(ReactionModelWrapper):
    """Class for implementing 'first-order adsorbate interaction model. 
    Should be sub-classed by scaler."""

    def __init__(self,reaction_model=None):
        if reaction_model is None:
            reaction_model = ReactionModel()
        self._rxm = reaction_model
        defaults = dict(
                cross_interaction_mode='geometric_mean',
                transition_state_cross_interaction_mode='intermediate_state',
                default_self_interaction_parameter = 0,
                interaction_response_function = 'linear',
                interaction_response_parameters = {'slope':1,'cutoff':0.25,
                    'smoothing':0.05},
                interaction_fitting_mode=None,
                input_overrides_fit = True, #user inputs override fitted values
                interaction_strength = 1,
                #weight interaction parameters by this
                numerical_delta_theta = None
                )
        self._rxm.update(defaults)
        self._required = {'cross_interaction_mode':str,
                'transition_state_cross_interaction_mode':str,
                'interaction_fitting_mode':None
                }
        self._log_strings = {'interaction_fitting_success':
                "interaction parameter ${param}@${surface} = ${paramval} (avg_error=${error})"}

    def parameterize_interactions(self):
        self._parameterized = True
        self.get_interaction_transition_state_scaling_matrix()
        
        if self.interaction_fitting_mode is not None:
            if '+' in self.interaction_fitting_mode:
                self.interaction_fitting_mode = self.interaction_fitting_mode.split('+')
            try:
                self.fit()
            except ValueError as inst:
                #This is due to poor error handling elsewhere... if self.fit() raises
                # a ValueError it looks like a convergence issue. While this hack 
                # is not ideal, it at least isolates the problem to within self.fit()
                raise UserWarning('Parameter fitting failed with ValueError: '+inst.message)

    def get_energy_error(self, epsilon_ij, theta, Ediff, Eint, 
            parameter_name, surface_name):

        i_surf = self.surface_names.index(surface_name)
        all_ads = self.adsorbate_names + self.transition_state_names
        E0_list = [self.species_definitions[sp]['formation_energy'][i_surf] 
                    for sp in all_ads]

        if None in E0_list:
            no_Es = [a for a,E in zip(all_ads,E0_list) if E is None]
            raise UserWarning('No adsorption energy specified for: '+str(no_Es))

        Ediff,Ediff_ads = Ediff
        ads_a, ads_b = parameter_name.split('&')
        if Ediff_ads == ads_a:
            i = all_ads.index(ads_a)
            j = all_ads.index(ads_b)
        elif Ediff_ads == ads_b:
            i = all_ads.index(ads_b)
            j = all_ads.index(ads_a)
        else:
            raise UserWarning('Differential adsorbate not in parameter name')

        if self.interaction_cross_term_names:
            param_names = self.adsorbate_names + self.interaction_cross_term_names
        else:
            param_names = self.adsorbate_names
        
        info = self.thermodynamics.adsorbate_interactions.get_interaction_info()
        params = [info[pi][i_surf] for pi in param_names]
        int_strength = self.interaction_strength
        self.interaction_strength = 1.0
        eps_matrix = self.params_to_matrix(E0_list+params)
        self.interaction_strength = int_strength
        eps_matrix[i,j] = eps_matrix[j,i] = epsilon_ij
        eps_list = list(eps_matrix.ravel())
        
        if 'differential' in self.interaction_fitting_mode:
            include_derivative = True
        else:
            include_derivative = False
        include_integral = False
        if 'integral' in self.interaction_fitting_mode:
            include_integral = True
        if 'numerical_differential' in self.interaction_fitting_mode:
            include_integral = True
        
        Eint_m,Ediff_m,Ejac = self.interaction_function(theta,E0_list,eps_list,
                    self.interaction_response_function,include_derivative,
                                                          include_integral)

        if 'numerical_differential' in self.interaction_fitting_mode:
            if 'differential' in self.interaction_fitting_mode:
                raise UserWarning('interaction_fitting_mode cannot contain both '
                        'numerical_differential and differential.')

            if 'delta_theta' in self.species_definitions[all_ads[i]]:
                delta_theta = self.species_definitions[all_ads[i]]['delta_theta']

            elif self.numerical_delta_theta is None:
                raise UserWarning('The change in coverage used to numerically compute '
                        'differential binding energies (numerical_delta_theta) '
                        'must be specified.')
            else:
                delta_theta = self.numerical_delta_theta

            new_theta = [x for x in theta]
            new_theta[i] -= delta_theta
            Eint_i = self.interaction_function(new_theta,E0_list,eps_list,
                    self.interaction_response_function,False,True)[0]
            Ediff_model = (Eint_m - Eint_i)/(delta_theta)
            diff_err = (Ediff-Ediff_model)

        elif 'differential' in self.interaction_fitting_mode:
            Ediff_model = Ediff_m[i]
            diff_err = (Ediff-Ediff_model)
        else:
            diff_err = None

        if 'integral' in self.interaction_fitting_mode:
            if Eint_m is None:
                raise AttributeError('Interaction model has no function for ' +\
                        'computing integral interacting energies. Fit to differential '+\
                        'energies instead.')
            int_err = (Eint-Eint_m)
        else:
            int_err = None

        return diff_err, int_err

    def error_norm(self,diff_err, int_err):
        if diff_err is None:
            return abs(int_err)
        elif int_err is None:
            return abs(diff_err)
        elif self.integral_error_fitting_weight:
            w = self.integral_error_fitting_weight
            return abs(diff_err) + w*abs(int_err)
        else:
            return abs(diff_err) + abs(int_err)

    def fit_interaction_parameter(self,theta_list,E_diffs,E_ints,param_name,surf_name):
        def target_function(eps_ij, theta_list, E_diffs,E_ints,param_name,surf_name):
            error = 0
            for theta,Ed,Ei in zip(theta_list,E_diffs,E_ints):
                d_err, i_err = self.get_energy_error(eps_ij, theta,Ed,Ei,
                        param_name,surf_name)
                err_norm = self.error_norm(d_err, i_err)
                error += err_norm
            return error/len(theta_list)
        
        minimize = lambda x: target_function(x,theta_list,E_diffs,E_ints,
                param_name,surf_name)

        x0 = getattr(self,'cross_interaction_initial_guess',0)
        if x0 is None:
            x0 = 0
         
        verbose_fitting = getattr(self,'verbose_interaction_parameter_fitting',False)
        eps_ij = fmin(minimize,[x0],disp=verbose_fitting)[0]
        err = minimize(eps_ij)
        self.log('interaction_fitting_success',
                param = param_name,
                paramval = str(eps_ij),
                error = str(round(err,2)),
                surface = surf_name
                )
        return eps_ij

    def fit(self):

        def required_params_to_fitting_info(param, key, required_params, fitting_info):
            #helper function to transfer data from "required_params_dict" to "fitting_info"
            theta_list,Ediffs,Eints = zip(*required_params_dict[key])
            eps_ij = self.fit_interaction_parameter(theta_list,Ediffs,Eints,param,surf)
            fitting_info[param] = [required_params_dict[key],eps_ij]
            del required_params_dict[key]

        def sync_with_species_defs(fitting_info,surf_id):
            #helper function to dump information into species_definitions
            for key in [k for k in fitting_info.keys() if '&' in k]:
                eps_ij = fitting_info[key][-1]
                ads_a,ads_b = sorted(key.split('&'))
                if ads_a == ads_b:
                    old_param = self.species_definitions[ads_a].get('self_interaction_parameter',None)
                    if (old_param is None):
                        if old_param is None:
                            old_param = [None]*len(self.surface_names)
                        old_param[surf_id] = eps_ij    
                        self.species_definitions[ads_a]['self_interaction_parameter'] = old_param
                else:
                    cross_ints_a = self.species_definitions[ads_a].get('cross_interaction_parameters',None)
                    cross_ints_b = self.species_definitions[ads_b].get('cross_interaction_parameters',None)
                    if cross_ints_a == cross_ints_b == None:
                        cross_ints = {}
                        cross_key = ads_a
                    elif cross_ints_b is None:
                        cross_ints = cross_ints_a
                        cross_key = ads_a
                    elif cross_ints_a is None:
                        cross_ints = cross_ints_b
                        cross_key = ads_b
                    else:
                        if ads_a in cross_ints_b:
                            cross_ints = cross_ints_b
                            cross_key = ads_b
                        else: 
                            cross_ints = cross_ints_a
                            cross_key = ads_a

                    if cross_key == ads_a:
                        other_key = ads_b
                    else:
                        other_key = ads_a

                    none_vec = [None]*len(self.surface_names)
                    if cross_key == ads_a:
                        cross_ints_vec = cross_ints.get(ads_b,none_vec)
                    elif cross_key == ads_b:
                        cross_ints_vec = cross_ints.get(ads_a,none_vec)
                    cross_ints_vec[surf_id] = eps_ij
                    cross_ints[other_key] = cross_ints_vec
                   
                    self.species_definitions[cross_key]['cross_interaction_parameters'] = cross_ints

        all_ads = self.adsorbate_names + self.transition_state_names
        all_fitting_info = {}

        for surf_id,surf in enumerate(self.surface_names):
            #collect all coverage/energy pairs
            all_cvgs = []
            for key in self.species_definitions:
                info = self.species_definitions[key]
                if info.get('coverage_dependent_energy',[None]*(surf_id+1))[surf_id]:
                    for cvg_i in info['coverage_dependent_energy'][surf_id]:
                        cvg,Ediff,Eint= cvg_i
                        Ediff = [Ediff,key]
                        all_cvgs.append([cvg,Ediff,Eint])

            required_params_dict = {}
            fitting_info = {}
            #organize coverages by which parameters the correspond to
            for cvg_Ed_Eint_i in all_cvgs:
                cvg_i,Ed_i,Eint_i = cvg_Ed_Eint_i
                reqd = str(self.required_interaction_parameters(cvg_i))
                if reqd not in required_params_dict:
                    required_params_dict[reqd] = [cvg_Ed_Eint_i]
                else:
                    if cvg_Ed_Eint_i not in required_params_dict[reqd]:
                        required_params_dict[reqd] += [cvg_Ed_Eint_i]
            
            for key in list(required_params_dict.keys()):
                print(key)
                keylist = eval(key)
                cvgs = required_params_dict[key]
                if len(keylist) == 1:
                    ads = keylist[0].split('&')[0]
                    #pull out low-coverage limit calcs since they are already used
                    if len(cvgs) == 1:
                        fitting_info[ads] = [cvgs,cvgs[0][1][0]]
                        del required_params_dict[key]
                    else:
                        for cvg in cvgs:
                            theta,Ed,Ei = cvg
                            Ed,Ed_ads = Ed
                            Ef = self.species_definitions[Ed_ads]['formation_energy'][surf_id]
                            if Ed == Ef:
                                fitting_info[ads] = [[cvg],Ed]
                                cvgs.remove(cvg)

            self_keys = [k for k in required_params_dict if len(eval(k)) == 1]
            cross_keys_3 = [k for k in required_params_dict if len(eval(k)) == 3]
            cross_keys_6 = [k for k in required_params_dict if len(eval(k)) == 6]

            user_inputs = []
            #get user-input values
            for ads in all_ads:
                param = '&'.join([ads,ads])
                default_vec = [self.default_self_interaction_parameter]*len(self.surface_names)
                if param not in fitting_info:
                    val = self.species_definitions[ads].get('self_interaction_parameter',
                                                    None)
                    if val is not None:
                        user_inputs.append(param)
                        val = val[surf_id]
                    else:
                        val = default_vec[surf_id]
                    #ensure that all self-interactions have a value
                    fitting_info[param] = [[],val]

            #get user-input cross-interactions
            for ads in all_ads:
                cross_ints = self.species_definitions[ads].get('cross_interaction_parameters',None)
                if cross_ints:
                    for key in cross_ints:
                        param1 = '&'.join([ads,key])
                        param2 = '&'.join([key,ads])
                        val = cross_ints[key][surf_id]
                        if val is not None:
                            user_inputs += [param1,param2]

            if not self.input_overrides_fit:
                user_inputs = []

            #fit known cross-interactions
            for key in self_keys:
                param = eval(key)[0]
                if param not in user_inputs:
                    required_params_to_fitting_info(param,key,required_params_dict,fitting_info)
          
            sync_with_species_defs(fitting_info,surf_id) 

            for key in cross_keys_3:
                params = eval(key)
                params = [pi for pi in params if pi not in list(fitting_info.keys())+user_inputs]
                assert len(params) <= 1 #if not then something is very weird.
                if params:
                    param = params[0]
                    if param not in user_inputs:
                        required_params_to_fitting_info(param,key,required_params_dict,fitting_info)
                else:
                    print('Using user-input value for '+eval(key)[1])

            sync_with_species_defs(fitting_info,surf_id) 

            underdetermined = cross_keys_6
            new_underdetermined = []
            initialized = False
            while underdetermined != new_underdetermined: 
                #hacky way of back-substituting... this can be cleaned up.
                if not initialized:
                    new_underdetermined = [k for k in underdetermined]
                    initialized = True
                underdetermined = [k for k in new_underdetermined]
                for key in underdetermined:
                    params = eval(key)
                    params = [pi for pi in params if pi not in fitting_info.keys()+user_inputs]
                    if len(params) == 1: #then we can do fitting
                        param = params[0]
                        if param not in user_inputs:
                            required_params_to_fitting_info(param,key,required_params_dict,fitting_info)
                        new_underdetermined.remove(key)
                    elif len(params) == 0:
                        cvg_states = []
                        for cvg,Ed,Ei in required_params_dict[key]:
                            cvg_state = ','.join(['='.join([a,str(c)]) for a,c in zip(all_ads,cvg) if c])
                            cvg_state += ' (E_'+Ed[-1]+'='+str(Ed[0])+',E_int='+str(Ei)+')'
                            cvg_states.append(cvg_state)
                        print('Warning: Not fitting to coverage-depedent energies for :'+'\n'.join(cvg_states)+'. '
                               'Parameters are already fully determined.')
                        new_underdetermined.remove(key)
                    else:
                        pass

                sync_with_species_defs(fitting_info,surf_id) 
            
            for key in underdetermined:
                missing = [k for k in eval(key) if k not in fitting_info]
                cvg_states = []
                for cvg,Ed,Ei in required_params_dict[key]:
                    cvg_state = ','.join(['='.join([a,str(c)]) for a,c in zip(all_ads,cvg) if c])
                    cvg_state += ' (E_'+Ed[-1]+'='+str(Ed[0])+',E_int='+str(Ei)+')'
                    cvg_states.append(cvg_state)
                state_str = '\n'.join(cvg_states)
                missing_str = ' and/or '.join(missing)
                print('Warning: Not fitting to coverage-depedent energies for:'+state_str+'. '
                        'Parameters are under-determined; add calculations to determine ' + \
                        missing_str +' to make this data useful.')

            all_fitting_info[surf] = fitting_info

        self.get_interaction_info() 
        self.fitting_info = all_fitting_info #store as attribute so that it can be double-checked externally
        
    def required_interaction_parameters(self,cvg):
        nonzeros = [i for i,c in enumerate(cvg) if c]
        all_combos = []
        for combo in combinations(nonzeros,2):
            all_combos.append('&'.join([self.parameter_names[i] for i in combo]))
        return all_combos

    def get_interaction_info(self):
        interaction_dict = {}
        n_ads = len(self.adsorbate_names)
        cross_term_names = []
        for a in self.adsorbate_names:
            interaction_dict[a] = self.species_definitions[a].get('self_interaction_parameter',None)
            if not interaction_dict[a]:
                eps_ii0 = self.default_self_interaction_parameter
                if a not in getattr(self,'_self_interaction_warned',[]):
                    print('Warning: No self-interaction parameter specified for '+ \
                            a+'. Assuming default ('+str(eps_ii0)+')')
                    self._self_interaction_warned = getattr(self,
                            '_self_interaction_warned',[]) + [a]
                interaction_dict[a] = [eps_ii0]*len(self.surface_names)
            cross_params = self.species_definitions[a].get(
                    'cross_interaction_parameters',{})
            for cp in cross_params:
                if cp not in self.adsorbate_names+self.transition_state_names:
                    raise ValueError(
                            'Cross parameter name must be in adsorbate names. The '+\
                            'name ' + cp + ' is not in ' +\
                            str(self.adsorbate_names+self.transition_state_names))
                params = cross_params[cp]
                if len(params) != len(self.surface_names):
                    raise ValueError('Cross parameters must be specified for each surface. '+\
                            str(len(params)) + ' parameters were specified, but there are '+\
                            str(len(self.surface_names)) + 'surfaces for '+ a +','+cp+'.' )
                ads_a, ads_b = sorted([a,cp]) #sort to avoid duplicates 
                name = '&'.join([ads_a,ads_b])
                if name not in cross_term_names:
                    cross_term_names.append(name)
                    interaction_dict[name] = params

        if cross_term_names:
            cross_term_names = tuple(cross_term_names) 
            self.interaction_cross_term_names = cross_term_names
    
        return interaction_dict

    def get_interaction_scaling_matrix(self):
        interaction_dict = self.get_interaction_info()

        cross_names = self.interaction_cross_term_names
        if cross_names:
            param_names = self.adsorbate_names + cross_names
        else:
            param_names = self.adsorbate_names

        constraint_dict = {}
        if not self.interaction_scaling_constraint_dict:
            self.interaction_scaling_constraint_dict = self.scaling_constraint_dict.copy()
            for ads in self.interaction_scaling_constraint_dict:
                if '-' not in ads:
                    constr = self.interaction_scaling_constraint_dict[ads]
                    new_constr = []
                    #preserve only 0 constraints
                    for ci in constr:
                        if ci != 0:
                            new_constr.append(None)
                        else:
                            new_constr.append(0)
                    constraint_dict[ads] = new_constr

        for ads in self.interaction_scaling_constraint_dict:
            if '&' in ads:
                a,b = ads.split('&')
                a,b = sorted([a,b])
                new_ads = '&'.join([a,b])
                constraint_dict[new_ads] = self.interaction_scaling_constraint_dict[ads]
            else:
                constraint_dict[ads] = self.interaction_scaling_constraint_dict[ads]

        #get mins/maxs
        interaction_mins = []
        interaction_maxs = []
        if self.default_interaction_constraints is None:
            self.default_interaction_constraints = [None]*(len(self.descriptor_names)+1)
        for p in param_names:
            if p not in constraint_dict:
                constr = self.default_interaction_constraints
            else:
                constr = constraint_dict[p]
            minvs,maxvs = parse_constraint(constr,p)
            interaction_mins.append(minvs)
            interaction_maxs.append(maxvs)

        #get coefficient matrix
        C = catmap.functions.scaling_coefficient_matrix(
                interaction_dict, self.descriptor_dict, 
                self.surface_names, 
                param_names,
                interaction_mins,
                interaction_maxs)
        self.interaction_coefficient_matrix = C.T

        return self.interaction_coefficient_matrix

    def get_interaction_transition_state_scaling_matrix(self):
        #get TS scaling matrix equivalent
        if self.transition_state_cross_interaction_mode == 'transition_state_scaling':
            if self.transition_state_scaling_matrix is None:
                raise AttributeError('Transition state scaling can only be used '+\
                ' for interactions if the transition_state_scaling_matrix is defined.')
            TS_weight_matrix = []
            for params in list(self.transition_state_scaling_matrix):
                TS_weight_matrix.append(params[:-1])

        else:
            #take TS interaction parameters as some linear combination of initial/final state
            if self.transition_state_cross_interaction_mode == 'initial_state':
                weight = 0
            elif self.transition_state_cross_interaction_mode == 'final_state':
                weight = 1
            elif self.transition_state_cross_interaction_mode.startswith('intermediate_state'):
                if self.transition_state_cross_interaction_mode == 'intermediate_state':
                    weight = 0.5
                else:
                    mname = self.transition_state_cross_interaction_mode
                    junk,w = mname.rsplit('(',1)
                    w,junk = w.rsplit(')',1)
                    try:
                        weight = float(w)
                    except:
                        raise AttributeError('Intermediate state weight is not floatable:'+w)
            else:
                raise AttributeError('Undefined transition_state_cross_interaction_mode: ' +\
                        self.transition_state_cross_interaction_mode)
            TS_weight_matrix = self.get_TS_weight_matrix(weight)

        self.interaction_transition_state_scaling_matrix = TS_weight_matrix

    def get_interaction_matrix(self,descriptors):
        full_descriptors = list(descriptors) + [1.]
        param_vector = np.dot(self.coefficient_matrix,full_descriptors)
        eps_vector = self.params_to_matrix(param_vector)
        return eps_vector

    def params_to_matrix(self,param_vector):
        n_ads = len(self.adsorbate_names)
        n_TS = len(self.transition_state_names)
        all_names = self.adsorbate_names + self.transition_state_names
        n_tot = len(all_names)
        epsilon_matrix = np.zeros((n_tot,n_tot))
        self_interactions = param_vector[n_tot:n_tot+n_ads]
        cross_interactions = param_vector[n_tot+n_ads:]
        for ads,param in zip(self.adsorbate_names,self_interactions):
            max_val = self.species_definitions[ads].get('max_self_interaction',None)
            if max_val is not None:
                param = min(param,max_val)
                self_interactions[self.adsorbate_names.index(ads)] = param

        for i,e_ii in enumerate(self_interactions):
            epsilon_matrix[i,i] = e_ii

        for i, e_i in enumerate(self_interactions):
            for j, e_j in enumerate(self_interactions):
                if not epsilon_matrix[i,j]:
                    if self.cross_interaction_mode == 'geometric_mean':
                        epsilon_matrix[i,j] = np.sqrt(abs(e_i)*abs(e_j))
                    elif self.cross_interaction_mode == 'arithmetic_mean':
                        epsilon_matrix[i,j] = (e_i+e_j)/2.
                    elif self.cross_interaction_mode == 'neglect':
                        epsilon_matrix[i,j] = 0

        if self.non_interacting_site_pairs:
            for site_1, site_2 in self.non_interacting_site_pairs:
                idxs1 = [self.adsorbate_names.index(sp) for 
                        sp in self.adsorbate_names if
                        self.species_definitions[sp]['site']==site_1]
                idxs2 = [self.adsorbate_names.index(sp) for 
                        sp in self.adsorbate_names if
                        self.species_definitions[sp]['site']==site_2]
                for i in idxs1:
                    for j in idxs2:
                        epsilon_matrix[i,j] = epsilon_matrix[j,i] = 0
        
        for i,TS_params in enumerate(list(self.interaction_transition_state_scaling_matrix)):
            i += n_ads
            TS_params = list(TS_params)
            for j, epsilon_params in enumerate(list(epsilon_matrix[0:n_ads,0:n_ads])):
                e_TS = np.dot(TS_params,epsilon_params)
                epsilon_matrix[i,j] = epsilon_matrix[j,i] = e_TS

        if len(cross_interactions)>0:
            for name, param in zip(self.interaction_cross_term_names,cross_interactions):
                a,b = name.split('&')
                i = all_names.index(a)
                j = all_names.index(b)
                epsilon_matrix[i,j] = epsilon_matrix[j,i] = param

        epsilon_matrix *= self.interaction_strength
        self._interaction_matrix = epsilon_matrix.tolist()
        return epsilon_matrix

    def get_TS_weight_matrix(self,weight):
        """Helper function to get `weights' of how
        to distribute TS-cross interactions between IS/FS.
        Should not be called externally."""

        weight_matrix = []
        for TS in self.transition_state_names:
            weight_TS = [0]*len(self.adsorbate_names)
            IS = None
            FS = None
            for rxn in self.elementary_rxns:
                if TS in rxn[1]:
                    if IS or FS:
                        print('Warning: Ambiguous initial/final state for '+TS)
                    IS = rxn[0]
                    FS = rxn[-1]
            for ads in IS:
                if ads in self.adsorbate_names:
                    idx = self.adsorbate_names.index(ads)
                    weight_TS[idx] += (1-weight)
            for ads in FS:
                if ads in self.adsorbate_names:
                    idx = self.adsorbate_names.index(ads)
                    weight_TS[idx] += weight
            weight_matrix.append(weight_TS)
        return weight_matrix

    @staticmethod
    def smooth_piecewise_linear_response(*args,**kwargs):
        return smooth_piecewise_linear(*args,**kwargs)[:2]

    @staticmethod
    def piecewise_linear_response(*args,**kwargs):
        kwargs['smoothing'] = 0
        return smooth_piecewise_linear(*args,**kwargs)[:2]

    @staticmethod
    def linear_response(*args,**kwargs):
        kwargs['smoothing'] = 0
        kwargs['cutoff'] = 0
        return smooth_piecewise_linear(*args,**kwargs)[:2]

    def fit_old(self): ##REMOVE AFTER WORKING
        all_adsnames = self.adsorbate_names+self.transition_state_names
        E0_list = []

        for i, surf in enumerate(self.surface_names):
            fitting_info = {}
            for sp in self.species_definitions:
                cvg_dep = self.species_definitions[sp].get('coverage_dependent_energy',
                        [None]*len(self.surface_names))[i]
                if cvg_dep:
                    for theta_E in cvg_dep:
                        theta,Ed,Ei = theta_E
                        n_theta = sum(ti != 0 for ti in theta)
                        ads_idx = all_adsnames.index(sp)
                        if n_theta == 1:
                            if sp in fitting_info:
                                if theta_E not in fitting_info[sp]:
                                    fitting_info[sp].append(theta_E)
                            else:
                                fitting_info[sp] = [theta_E]
                        elif n_theta == 2:
                            coads_idx = [j for j,cvg in enumerate(theta) if cvg and j != ads_idx][0]
                            coads = all_adsnames[coads_idx]
                            ads_a, ads_b = sorted([sp,coads])
                            name = '&'.join([ads_a,ads_b])
                            if name in fitting_info:
                                if theta_E not in fitting_info[name]:
                                    fitting_info[name].append(theta_E)
                            else:
                                fitting_info[name] = [theta_E]
                        else:
                            print('Warning: Ignoring coverage dependent entry for '+sp+\
                                   '. Only 1-2 adsorbates can have non-zero coverages')
            
            for key in fitting_info.keys():
                #clean up fitting_info
                if '&' not in key and len(fitting_info[key]) == 1:
                    #move into cross-parameter info if possible
                    #This may not be totally general, but the goal
                    #is to ensure that cross-parameters are
                    #also fit to the low-coverage limit
                    cross_keys = [k for k in fitting_info.keys() if (
                        '&' in k and key in k.split('&'))]
                    for k in cross_keys:
                        fitting_info[k] = fitting_info[key] + fitting_info[k]

                    del fitting_info[key] #a single entry is the same one used for Ef
                    
                elif self.input_overrides_fit: 
                    #user inputs should over-ride fit, therefore delete any fitting
                    #data for parameters the user already input
                    if '&' not in key:
                        if self.species_definitions[key].get(
                                'self_interaction_parameter',
                                [None]*len(self.surface_names))[i] is not None:
                            del fitting_info[key]
                    else:
                        ads_a, ads_b = key.split('&')
                        #have to check for cross interaction specified either way
                        ab_crossint = self.species_definitions[ads_a].get(
                               'cross_interaction_parameters',{})
                        if ads_b in ab_crossint:
                            del fitting_info[key]
                        ba_crossint = self.species_definitions[ads_b].get(
                               'cross_interaction_parameters',{})
                        if ads_a in ba_crossint:
                            del fitting_info[key]

            self_params = [k for k in fitting_info if '&' not in k]
            cross_params = [k for k in fitting_info if '&' in k]
            E0_list = [self.species_definitions[sp]['formation_energy'][i]  \
                    for sp in all_adsnames]
            for key in self_params:
                thetas,Ediffs,Eints = zip(*fitting_info[key])
                eps_ii = self.fit_interaction_parameter(thetas,Ediffs,Eints,key,surf)
                eps = self.species_definitions[key].get('self_interaction_parameter',
                        [None]*len(self.surface_names))
                eps[i] = eps_ii
                self.species_definitions[key]['self_interaction_parameter'] = eps
            
            for key in cross_params:
                thetas,Ediffs,Eints = zip(*fitting_info[key])
                ads_a, ads_b = key.split('&')
                idx_a = all_adsnames.index(ads_a)
                idx_b = all_adsnames.index(ads_b)
                thetas_a = [ti[idx_a] for ti in thetas]
                thetas_b = [ti[idx_b] for ti in thetas]
                stdev_a = np.std(thetas_a)
                stdev_b = np.std(thetas_b)
                if 0 not in [stdev_a,stdev_b]:
                    if 'differential' in self.interaction_fitting_mode:
                        print('Warning: no constant coverage for fitting '+\
                                'cross-interaction parameter '+key+'. Using '+\
                                'adsorbate with minimum standard deviation '+\
                                'of coverages as differential adsorbate')

                if stdev_b < stdev_a:    
                    #switch adsorbate so that the one at constant coverage
                    #is the adsorbate used for differential adsorption energies.
                    ads_a, ads_b = ads_b, ads_a

                fit_param = '&'.join([ads_a,ads_b])
                eps_ij = self.fit_interaction_parameter(thetas,Ediffs,Eints,
                        fit_param,surf)
                eps_dict_b = self.species_definitions[ads_b].get(
                        'cross_interaction_parameters',
                        {})
                if ads_a not in eps_dict_b:
                    eps_dict_b[ads_a] = [None]*len(self.surface_names)
                eps_dict_b[ads_a][i] = eps_ij
                self.species_definitions[ads_b][
                        'cross_interaction_parameters'] = eps_dict_b
        
        #make sure that new parameters get incorporated into interaction matrix
        self.get_interaction_info()
        self._fitting_info = fitting_info
