<commit_before>#! /usr/bin/env python
"""Creates source code for calculating analytical Jacobian matrix.
"""

# Python 2 compatibility
from __future__ import division
from __future__ import print_function

# Standard libraries
import sys
import math
from argparse import ArgumentParser
import os

# Local imports
import chem_utilities as chem
import mech_interpret as mech
import rate_subs as rate
import utils
import mech_auxiliary as aux
import CUDAParams
import cache_optimizer as cache
import shared_memory as shared


def calculate_shared_memory(rind, rxn, specs, reacs, rev_reacs, pdep_reacs):
    # need to figure out shared memory stuff
    variable_list = []
    usages = []
    fwd_usage = 2
    rev_usage = 0 if not rxn.rev else 2
    pres_mod_usage = 0 if not (rxn.pdep or rxn.thd_body) else (4 if rxn.thd_body else 2)
    reac_usages = [0 for i in range(len(rxn.reac))]
    prod_usages = [0 for i in range(len(rxn.prod))]
    # add variables
    variable_list.append(utils.get_array('cuda', 'fwd_rates', rind))
    if rxn.rev:
        variable_list.append(utils.get_array('cuda', 'rev_rates', rev_reacs.index(rind)))
    if rxn.pdep or rxn.thd_body:
        variable_list.append(utils.get_array('cuda', 'pres_mod', pdep_reacs.index(rind)))
    for sp in set(rxn.reac + rxn.prod + [x[0] for x in rxn.thd_body_eff]):
        variable_list.append(utils.get_array('cuda', 'conc', sp))

    alphaij_count = 0
    # calculate usages
    if rxn.thd_body or rxn.pdep:
        fwd_usage += 1
        if rxn.rev:
            rev_usage += 1
        for i, thd in enumerate(rxn.thd_body_eff):
            # check alphaij
            alphaij = thd[1]
            if alphaij is not None and alphaij != 1.0:
                alphaij_count += 1
        fwd_usage += alphaij_count
        if rxn.rev:
            rev_usage += alphaij_count

    for i, sp in enumerate(rxn.reac):
        nu = rxn.reac_nu[i]
        if nu - 1 > 0:
            reac_usages[i] += 1
            if rxn.thd_body:
                pres_mod_usage += 1
        for i2, sp2 in enumerate(rxn.reac):
            if sp == sp2:
                continue
            reac_usages[i2] += 1

    if rxn.rev:
        for i, sp in enumerate(rxn.prod):
            nu = rxn.prod_nu[i]
            if nu - 1 > 0:
                prod_usages[i] += 1
            for i2, sp2 in enumerate(rxn.prod):
                if sp == sp2:
                    continue
                prod_usages[i2] += 1

    usages.append(fwd_usage)
    if rxn.rev:
        usages.append(rev_usage)
    if rxn.pdep or rxn.thd_body:
        usages.append(pres_mod_usage)
    for sp in set(rxn.reac + rxn.prod + [x[0] for x in rxn.thd_body_eff]):
        u = 0
        if sp in rxn.reac:
            u += reac_usages[rxn.reac.index(sp)]
        if sp in rxn.prod:
            u += prod_usages[rxn.prod.index(sp)]
        if sp in rxn.thd_body_eff:
            u += 1
        usages.append(u)

    if CUDAParams.JacRateStrat == CUDAParams.JacRatesCacheStrat.Exclude:
        if rxn.rev:
            usages[0] = 0
            usages[1] = 0
        else:
            usages[0] = 0

    return variable_list, usages


def get_net_rate_string(lang, rxn, rind, rev_reacs, get_array):
    jline = ''
    if rxn.rev:
        jline += '(' + get_array(lang, 'fwd_rates', rind)
        jline += ' - ' + \
                 get_array(lang, 'rev_rates', rev_reacs.index(rind))
        jline += ')'
    else:
        jline += get_array(lang, 'fwd_rates', rind)
    return jline


def write_dr_dy(file, lang, rev_reacs, rxn, rind, pind, nspec, get_array):
    # write the T_Pr and T_Fi terms if needed
    if rxn.pdep:
        jline = utils.line_start + 'pres_mod_temp = '
        jline += get_array(lang, 'pres_mod', pind) + ' * ('
        # dPr/dYj contribution
        if rxn.low:
            # unimolecular/recombination
            jline += '(1.0 / (1.0 + Pr))'
        elif rxn.high:
            # chem-activated bimolecular
            jline += '(-Pr / (1.0 + Pr))'
        if rxn.troe:
            jline += (' - Pr * log(Fcent) * 2.0 * A * (B * '
                      '{:.16}'.format(1.0 / math.log(10.0)) +
                      ' + A * '
                      '{:.16}) / '.format(0.14 / math.log(10.0)) +
                      '(B * B * B * (1.0 + A * A / (B * B)) '
                      '* (1.0 + A * A / (B * B)))'
                      )
        elif rxn.sri:
            jline += ('- Pr * X * X * '
                      '{:.16} * '.format(2.0 / math.log(10.0)) +
                      'log10(Pr) * '
                      'log({:.4} * '.format(rxn.sri_par[0]) +
                      'exp({:.4} / T) + '.format(-rxn.sri_par[1]) +
                      'exp(T / {:.4}))'.format(-rxn.sri_par[2])
                      )

        jline += ') * ' + get_net_rate_string(lang, rxn, rind, rev_reacs, get_array)
        file.write(jline + utils.line_end[lang])

    file.write('  j_temp = -mw_avg * rho_inv * (')
    jline = ''
    # next, contribution from dR/dYj
    # namely the T_dy independent term
    if rxn.pdep or rxn.thd_body:
        jline += get_array(lang, 'pres_mod', pind)
        jline += ' * ('

    reac_nu = 0
    prod_nu = 0
    if rxn.thd_body and not rxn.pdep:
        reac_nu = 1
        if rxn.rev:
            prod_nu = 1

    # get reac and prod nu sums
    reac_nu += sum(rxn.reac_nu)

    if rxn.rev:
        prod_nu += sum(rxn.prod_nu)

    if reac_nu != 0:
        if reac_nu == -1:
            jline += '-'
        elif reac_nu != 1:
            jline += '{} * '.format(float(reac_nu))
        jline += '' + get_array(lang, 'fwd_rates', rind)

    if prod_nu != 0:
        if prod_nu == 1:
            jline += ' - '
        elif prod_nu == -1:
            jline += ' + '
        else:
            jline += ' - {} * '.format(float(prod_nu))
        jline += '' + get_array(lang, 'rev_rates', rev_reacs.index(rind))

    # find alphaij_hat
    alphaij_hat = None
    counter = {}
    counter[1.0] = 0
    if rxn.thd_body_eff:
        for spec, efficiency in rxn.thd_body_eff:
            if not efficiency in counter:
                counter[efficiency] = 0
            counter[efficiency] += 1
        counter[1.0] += (nspec - len(rxn.thd_body_eff))
        alphaij_hat = max(counter.keys(), key=lambda x: counter[x])

    # now handle third body / pdep parts if needed
    if rxn.thd_body and not rxn.pdep:<commit_msg>unnecessary<commit_after>#! /usr/bin/env python
"""Creates source code for calculating analytical Jacobian matrix.
"""

# Python 2 compatibility
from __future__ import division
from __future__ import print_function

# Standard libraries
import sys
import math
from argparse import ArgumentParser
import os

# Local imports
import chem_utilities as chem
import mech_interpret as mech
import rate_subs as rate
import utils
import mech_auxiliary as aux
import CUDAParams
import cache_optimizer as cache
import shared_memory as shared


def calculate_shared_memory(rind, rxn, specs, reacs, rev_reacs, pdep_reacs):
    # need to figure out shared memory stuff
    variable_list = []
    usages = []
    fwd_usage = 2
    rev_usage = 0 if not rxn.rev else 2
    pres_mod_usage = 0 if not (rxn.pdep or rxn.thd_body) else (4 if rxn.thd_body else 2)
    reac_usages = [0 for i in range(len(rxn.reac))]
    prod_usages = [0 for i in range(len(rxn.prod))]
    # add variables
    variable_list.append(utils.get_array('cuda', 'fwd_rates', rind))
    if rxn.rev:
        variable_list.append(utils.get_array('cuda', 'rev_rates', rev_reacs.index(rind)))
    if rxn.pdep or rxn.thd_body:
        variable_list.append(utils.get_array('cuda', 'pres_mod', pdep_reacs.index(rind)))
    for sp in set(rxn.reac + rxn.prod + [x[0] for x in rxn.thd_body_eff]):
        variable_list.append(utils.get_array('cuda', 'conc', sp))

    alphaij_count = 0
    # calculate usages
    if rxn.thd_body or rxn.pdep:
        fwd_usage += 1
        if rxn.rev:
            rev_usage += 1
        for i, thd in enumerate(rxn.thd_body_eff):
            # check alphaij
            alphaij = thd[1]
            if alphaij is not None and alphaij != 1.0:
                alphaij_count += 1
        fwd_usage += alphaij_count
        if rxn.rev:
            rev_usage += alphaij_count

    for i, sp in enumerate(rxn.reac):
        nu = rxn.reac_nu[i]
        if nu - 1 > 0:
            reac_usages[i] += 1
            if rxn.thd_body:
                pres_mod_usage += 1
        for i2, sp2 in enumerate(rxn.reac):
            if sp == sp2:
                continue
            reac_usages[i2] += 1

    if rxn.rev:
        for i, sp in enumerate(rxn.prod):
            nu = rxn.prod_nu[i]
            if nu - 1 > 0:
                prod_usages[i] += 1
            for i2, sp2 in enumerate(rxn.prod):
                if sp == sp2:
                    continue
                prod_usages[i2] += 1

    usages.append(fwd_usage)
    if rxn.rev:
        usages.append(rev_usage)
    if rxn.pdep or rxn.thd_body:
        usages.append(pres_mod_usage)
    for sp in set(rxn.reac + rxn.prod + [x[0] for x in rxn.thd_body_eff]):
        u = 0
        if sp in rxn.reac:
            u += reac_usages[rxn.reac.index(sp)]
        if sp in rxn.prod:
            u += prod_usages[rxn.prod.index(sp)]
        if sp in rxn.thd_body_eff:
            u += 1
        usages.append(u)

    if CUDAParams.JacRateStrat == CUDAParams.JacRatesCacheStrat.Exclude:
        if rxn.rev:
            usages[0] = 0
            usages[1] = 0
        else:
            usages[0] = 0

    return variable_list, usages


def get_net_rate_string(lang, rxn, rind, rev_reacs, get_array):
    jline = ''
    if rxn.rev:
        jline += '(' + get_array(lang, 'fwd_rates', rind)
        jline += ' - ' + \
                 get_array(lang, 'rev_rates', rev_reacs.index(rind))
        jline += ')'
    else:
        jline += get_array(lang, 'fwd_rates', rind)
    return jline


def write_dr_dy(file, lang, rev_reacs, rxn, rind, pind, nspec, get_array):
    # write the T_Pr and T_Fi terms if needed
    if rxn.pdep:
        jline = utils.line_start + 'pres_mod_temp = '
        jline += get_array(lang, 'pres_mod', pind) + ' * ('
        # dPr/dYj contribution
        if rxn.low:
            # unimolecular/recombination
            jline += '(1.0 / (1.0 + Pr))'
        elif rxn.high:
            # chem-activated bimolecular
            jline += '(-Pr / (1.0 + Pr))'
        if rxn.troe:
            jline += (' - Pr * log(Fcent) * 2.0 * A * (B * '
                      '{:.16}'.format(1.0 / math.log(10.0)) +
                      ' + A * '
                      '{:.16}) / '.format(0.14 / math.log(10.0)) +
                      '(B * B * B * (1.0 + A * A / (B * B)) '
                      '* (1.0 + A * A / (B * B)))'
                      )
        elif rxn.sri:
            jline += ('- Pr * X * X * '
                      '{:.16} * '.format(2.0 / math.log(10.0)) +
                      'log10(Pr) * '
                      'log({:.4} * '.format(rxn.sri_par[0]) +
                      'exp({:.4} / T) + '.format(-rxn.sri_par[1]) +
                      'exp(T / {:.4}))'.format(-rxn.sri_par[2])
                      )

        jline += ') * ' + get_net_rate_string(lang, rxn, rind, rev_reacs, get_array)
        file.write(jline + utils.line_end[lang])

    file.write('  j_temp = -mw_avg * rho_inv * (')
    jline = ''
    # next, contribution from dR/dYj
    # namely the T_dy independent term
    if rxn.pdep or rxn.thd_body:
        jline += get_array(lang, 'pres_mod', pind)
        jline += ' * ('

    reac_nu = 0
    prod_nu = 0
    if rxn.thd_body and not rxn.pdep:
        reac_nu = 1
        if rxn.rev:
            prod_nu = 1

    # get reac and prod nu sums
    reac_nu += sum(rxn.reac_nu)

    if rxn.rev:
        prod_nu += sum(rxn.prod_nu)

    if reac_nu != 0:
        if reac_nu == -1:
            jline += '-'
        elif reac_nu != 1:
            jline += '{} * '.format(float(reac_nu))
        jline += '' + get_array(lang, 'fwd_rates', rind)

    if prod_nu != 0:
        if prod_nu == 1:
            jline += ' - '
        elif prod_nu == -1:
            jline += ' + '
        else:
            jline += ' - {} * '.format(float(prod_nu))
        jline += '' + get_array(lang, 'rev_rates', rev_reacs.index(rind))

    # find alphaij_hat
    alphaij_hat = None
    counter = {}
    counter[1.0] = 0
    if rxn.thd_body_eff:
        for spec, efficiency in rxn.thd_body_eff:
            if not efficiency in counter:
                counter[efficiency] = 0
            counter[efficiency] += 1
        counter[1.0] += (nspec - len(rxn.thd_body_eff))
        alphaij_hat = max(counter.keys(), key=lambda x: counter[x])

    # now handle third body / pdep parts if needed
    if rxn.thd_body:
{'fraction_non_alphanumeric': 0.12540294048274236, 'fraction_numerical': 0.016038996776476137, 'mean_word_length': 1.7598177478845736, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2311822', 'n_tokens_mistral': 4755, 'n_tokens_neox': 4530, 'n_words': 1076}