Bug in GPUFOHelasCallWriter format_coupling (wrong params2order for SUSY processes)

[valassi]

Hi @oliviermattelaer I think I hit a bug in your original implementation of GPUFOHelasCallWriter format_coupling (which I then overloaded in the cudacpp plugin, but keeping the same logic and hence the same bug).

The result of the bug is that params2order has weird values like the following, where two parameters have the same index,

{'mdl_Msu3': 3, 'mdl_Wsu3': 1, 'mdl_Msu6': 2, 'mdl_Wsu6': 3}

while it should instead be (I think)

{'mdl_Msu3': 0, 'mdl_Wsu3': 1, 'mdl_Msu6': 2, 'mdl_Wsu6': 3}

The problem is this code

    def format_coupling(self, call):
        """Format the coupling so any minus signs are put in front"""
        import re
        ###print(call) # FOR DEBUGGING
        model = self.get('model')
        newcoup = False
        if not hasattr(self, 'couplings2order'):
            self.couplings2order = {}
            self.params2order = {}
        if not hasattr(self, 'couporderdep'):
            self.couporderdep = {}
            self.couporderindep = {}
        for coup in re.findall(self.findcoupling, call):
            if coup == 'ZERO':
                ###call = call.replace('pars->ZERO', '0.')
                call = call.replace('m_pars->ZERO', '0.') # AV
                continue
            sign = ''
            if coup.startswith('-'):
                sign = '-'
                coup = coup[1:]
            try:
                param = model.get_parameter(coup)
            except KeyError:
                param = False
            if param:
                alias = self.params2order
                name = 'cIPD'
            elif model.is_running_coupling(coup):
                alias = self.couporderdep
                name = 'cIPC'
            else:
                alias = self.couporderindep
                name = 'cIPC'
            if coup not in alias:
                if alias == self.couporderindep:
                    if not len(alias):
                        alias[coup] = len(self.couporderdep)
                    else:
                        alias[coup] = alias[list(alias)[-1]]+1
                else:
                    alias[coup] = len(alias)
                if alias == self.couporderdep:
                    for k in self.couporderindep:
                        self.couporderindep[k] += 1
                newcoup = True
            if name == 'cIPD':
                call = call.replace('m_pars->%s%s' % (sign, coup),
                                    '%s%s[%s]' % (sign, name, alias[coup]))                        
            elif model.is_running_coupling(coup):
                ###call = call.replace('m_pars->%s%s' % (sign, coup),
                ###                    '%scxmake( cIPC[%s], cIPC[%s] )' %
                ###                    (sign, 2*alias[coup],2*alias[coup]+1))
                ###misc.sprint(name, alias[coup])
                # AV from cIPCs to COUP array (running alphas #373)
                # OM fix handling of 'unary minus' #628
                call = call.replace('CI_ACCESS', 'CD_ACCESS')
                call = call.replace('m_pars->%s%s' % (sign, coup),
                                    'COUPs[%s], %s' % (alias[coup], '1.0' if not sign else '-1.0')) 
            else:
                call = call.replace('CD_ACCESS', 'CI_ACCESS')
                call = call.replace('m_pars->%s%s' % (sign, coup),
                                    'COUPs[ndcoup + %s], %s' % (alias[coup]-len(self.couporderdep), '1.0' if not sign else '-1.0'))
            if newcoup:
                self.couplings2order = self.couporderdep | self.couporderindep
        return call

The bug specifically are the two checks alias == self.couporderindep and alias == self.couporderdep. While alias is essentially a pointer or a reference to an existing dictionary, in python this becomes a check on the contents of the two dictionaries. So if alias points to another dictionary, but both the pointed alias and the checked dictionary are empty, then the result is True while (I THINK!) it is meant to be False.

My fix is simply the following:

            if param:
                alias = self.params2order
                aliastxt = 'PARAM'
                name = 'cIPD'
            elif model.is_running_coupling(coup):
                alias = self.couporderdep
                aliastxt = 'COUPD'
                name = 'cIPC'
            else:
                alias = self.couporderindep
                aliastxt = 'COUPI'
                name = 'cIPC'
            if coup not in alias:
                ###if alias == self.couporderindep: # bug! this is incorrectly true when all all dictionaries are empty!
                if aliastxt == 'COUPI':
                    if not len(alias):
                        alias[coup] = len(self.couporderdep)
                    else:
                        alias[coup] = alias[list(alias)[-1]]+1
                else:
                    alias[coup] = len(alias)
                ###if alias == self.couporderdep: # bug! this is incorrectly true when all all dictionaries are empty!
                if aliastxt == 'COUPD':
                    for k in self.couporderindep:
                        self.couporderindep[k] += 1
                newcoup = True

which is ugly but effective.

This will be part of an upcoming PR for #820, extending PR #625 (or most likely I will put everything in that one). But I will also make this a small separate PR.

Or in any case I want to be sure and ask you, do I understand correctly the logic of the above and is my fix correct?

Thanks Andrea

cc @roiser

[valassi]

This is fixed in PR #822 and/or PR #824, closing for simplicity