<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) else (4 if rxn.thd 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(rxn)))
    if rxn.pdep or rxn.thd:
        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]):
        sp_real = next(spec for spec in specs if spec.name == sp)
        variable_list.append(utils.get_array('cuda', 'conc', specs.index(sp_real)))

    alphaij_count = 0
    #calculate usages
    if rxn.thd or rxn.pdep:
        fwd_usage +=1
        if rxn.rev:
            rev_usage += 1
        for i, thd in enumerate(rxn.thd_body):
            #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 sp_name in rxn.reac:
        nu = rxn.reac_nu[rxn.reac.index(sp_name)]
        if nu - 1 > 0:
            reac_usages[rxn.reac.index(sp_name)] += 1
            if rxn.thd:
                pres_mod_usage += 1
        for sp_name_2 in rxn.reac:
            if sp_name == sp_name_2:
                continue
            reac_usages[rxn.reac.index(sp_name_2)] += 1

    if rxn.rev:
        for sp_name in rxn.prod:
            nu = rxn.prod_nu[rxn.prod.index(sp_name)]
            if nu - 1 > 0:
                prod_usages[rxn.prod.index(sp_name)] += 1
            for sp_name_2 in rxn.prod:
                if sp_name == sp_name_2:
                    continue
                prod_usages[rxn.prod.index(sp_name_2)] += 1

    usages.append(fwd_usage)
    if rxn.rev:
        usages.append(rev_usage)
    if rxn.pdep or rxn.thd:
        usages.append(pres_mod_usage)
    for sp in set(rxn.reac + rxn.prod + [x[0] for x in rxn.thd_body]):
        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:
            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(rxn))
        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 = '  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 += (' - log(Fcent) * 2.0 * A * (B * '
                      '{:.6}'.format(1.0 / math.log(10.0)) +
                      ' + A * '
                      '{:.6}) / '.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 += ('- X * X * '
                      '{:.6} * '.format(2.0 / math.log(10.0)) +
                      'log10(Pr) * '
                      'log({:.4} * '.format(rxn.sri[0]) +
                      'exp(-{:4} / T) + '.format(rxn.sri[1]) +
                      'exp(-T / {:.4}))'.format(rxn.sri[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:
        jline += get_array(lang, 'pres_mod', pind)
        jline += ' * ('

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

    #get reac and prod nu sums
    for sp_name in rxn.reac:
        nu = rxn.reac_nu[rxn.reac.index(sp_name)]
        if nu == 0:
            continue
        reac_nu += nu

    if rxn.rev:
        for sp_name in rxn.prod:
            nu = rxn.prod_nu[rxn.prod.index(sp_name)]
            if nu == 0:
                continue
            prod_nu += 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(rxn))
  
    #find alphaij_hat
    alphaij_hat = None
    counter = {}
    counter[1.0] = 0
    if rxn.thd_body:
        for spec, efficiency in rxn.thd_body:
            if not efficiency in counter:
                counter[efficiency] = 0
            counter[efficiency] += 1
        counter[1.0] += (nspec - sum(counter.values()))
        alphaij_hat = max(counter.keys(), key=lambda x:counter[x])

    #now handle third body / pdep parts if needed
    if rxn.thd and not rxn.pdep:
        jline += '))'
        if alphaij_hat is not None:
            if alphaij_hat == 1.0:
                jline += ' + '
            elif alphaij_hat == -1.0:
                jline += ' - '
            else:
                jline += ' + {} * '.format(alphaij_hat) 
            jline += get_net_rate_string(lang, rxn, rind, rev_reacs, get_array)
    elif rxn.pdep:
        jline += ') + pres_mod_temp'
        jline += ')'
        if alphaij_hat is not None:
            if alphaij_hat == 1.0:
                jline += ' + '
            elif alphaij_hat == -1.0:
                jline += ' - '
            else:
                jline += ' + {} * '.format(alphaij_hat) 
            jline += '(pres_mod_temp / conc_temp)'
    else:
        jline += ')'
    file.write(jline + utils.line_end[lang])

    if rxn.pdep and not rxn.thd and not rxn.pdep_sp:
        file.write('  pres_mod_temp /= conc_temp' + utils.line_end[lang])
    return alphaij_hat

def write_rates(file, lang, rxn):
    file.write('  kf = ' + rate.rxn_rate_const(rxn.A, rxn.b, rxn.E) +
                utils.line_end[lang])
    if rxn.rev and not rxn.rev_par:
        file.write('  kr = kf / Kc' + utils.line_end[lang])
    elif rxn.rev_par:
        file.write('  kr = ' +
        rate.rxn_rate_const(rxn.rev_par[0],
                                         rxn.rev_par[1],
                                         rxn.rev_par[2]
                                         ) +
        utils.line_end[lang])
                                      
def write_dr_dy_species(lang, specs, rxn, pind, j_sp, sp_j, alphaij_hat, rind, rev_reacs, get_array):
    jline = 'j_temp'
    if rxn.pdep and not rxn.pdep_sp:
        alphaij = next((thd[1] for thd in rxn.thd_body
                        if thd[0] == sp_j.name), None)
        if alphaij is None:
            alphaij = 1.0
        if alphaij != alphaij_hat:
            diff = alphaij - alphaij_hat
            if diff != 1:
                if diff == -1:
                    jline += ' - pres_mod_temp'
                else:
                    jline += ' + {} * pres_mod_temp'.format(diff)
            else:
                jline += ' + pres_mod_temp'
    elif rxn.pdep and rxn.pdep_sp and rxn.pdep_sp == sp_j.name:
        jline += ' + pres_mod_temp / (rho * {})'.format(get_array(lang, 'y', j_sp))
    elif rxn.thd and not rxn.pdep:
        alphaij = next((thd[1] for thd in rxn.thd_body
                        if thd[0] == sp_j.name), None)
        if alphaij is None:
            alphaij = 1.0
        if alphaij != alphaij_hat:
            diff = alphaij - alphaij_hat
            if diff != 1:
                if diff == -1:
                    jline += ' - '
                else:
                    jline += ' + {} * '.format(diff)
            else:
                jline += ' + '
            jline += get_net_rate_string(lang, rxn, rind, rev_reacs, get_array)

    if (rxn.pdep or rxn.thd) and (sp_j.name in rxn.reac or (rxn.rev and sp_j.name in rxn.prod)):
        jline += ' + ' + get_array(lang, 'pres_mod', pind)
        jline += ' * '

    if sp_j.name in rxn.reac:
        if not rxn.pdep and not rxn.thd:
            jline += ' + '
        nu = rxn.reac_nu[rxn.reac.index(sp_j.name)]
        if nu != 1:
            if nu == -1:
                jline += '-'
            else:
                jline += '{} * '.format(float(nu))
        jline += 'kf'
        nu_temp = rxn.reac_nu[rxn.reac.index(sp_j.name)]
        if (nu_temp - 1) > 0:
            if isinstance(nu_temp - 1, float):
                jline += ' * pow(' + \
                    get_array(lang, 'conc', j_sp)
                jline += ', {})'.format(nu_temp - 1)
            else:
                # integer, so just use multiplication
                for i in range(nu_temp - 1):
                    jline += ' * ' + \
                        get_array(lang, 'conc', j_sp)

        # loop through remaining reactants
        for sp_reac in rxn.reac:
            if sp_reac == sp_j.name:
                continue

            nu_temp = rxn.reac_nu[rxn.reac.index(sp_reac)]
            isp = next(i for i in range(len(specs))
                       if specs[i].name == sp_reac)
            if isinstance(nu_temp, float):
                jline += ' * pow(conc' + \
                    get_array(lang, isp)
                jline += ', ' + str(nu_temp) + ')'
            else:
                # integer, so just use multiplication
                for i in range(nu_temp):
                    jline += ' * ' + \
                        get_array(lang, 'conc', isp)

    if rxn.rev and sp_j.name in rxn.prod:
        if not rxn.pdep and not rxn.thd and not sp_j.name in rxn.reac:
            jline += ' + '
        elif sp_j.name in rxn.reac:
            jline += ' + '

        nu = rxn.prod_nu[rxn.prod.index(sp_j.name)]
        if nu != -1:
            if nu == 1:
                jline += '-'
            else:
                jline += '{} * '.format(float(-1 * nu))


        jline += 'kr'
        temp_nu = rxn.prod_nu[rxn.prod.index(sp_j.name)]
        if (temp_nu - 1) > 0:
            if isinstance(temp_nu - 1, float):
                jline += ' * pow(' + \
                    get_array(lang, 'conc', j_sp)
                jline += ', {})'.format(temp_nu - 1)
            else:
                # integer, so just use multiplication
                for i in range(temp_nu - 1):
                    jline += ' * ' + \
                        get_array(lang, 'conc', j_sp)

        # loop through remaining products
        for sp_prod in rxn.prod:
            if sp_prod == sp_j.name:
                continue

            temp_nu = rxn.prod_nu[rxn.prod.index(sp_prod)]
            isp = next(i for i in range(len(specs))
                       if specs[i].name == sp_prod)
            if isinstance(temp_nu, float):
                jline += ' * pow(' + \
                    get_array(lang, 'conc', isp)
                jline += ', {})'.format(temp_nu)
            else:
                # integer, so just use multiplication
                jline += ''.join([' * ' + get_array(lang, 'conc', isp) for i in range(temp_nu)])
    
    return jline

def write_kc(file, lang, specs, rxn):
    sum_nu = 0
    coeffs = {}
    for sp_name in set(rxn.reac + rxn.prod):
        spec = next((spec for spec in specs
            if spec.name == sp_name), None)
        nu = get_nu(spec, rxn)

        if nu == 0:
            continue

        sum_nu += nu

        lo_array = [utils.round_sig(nu, 3)] + [utils.round_sig(x, 9) for x in [
                    spec.lo[6], spec.lo[0], spec.lo[0] - 1.0, spec.lo[1] / 2.0,
                    spec.lo[2] / 6.0, spec.lo[3] / 12.0, spec.lo[4] / 20.0,
                    spec.lo[5]]
                ]
        lo_array = [x * lo_array[0] for x in [lo_array[1] - lo_array[2]] + lo_array[3:]]

        hi_array = [utils.round_sig(nu, 3)] + [utils.round_sig(x, 9) for x in [
                    spec.hi[6], spec.hi[0], spec.hi[0] - 1.0, spec.hi[1] / 2.0,
                    spec.hi[2] / 6.0, spec.hi[3] / 12.0, spec.hi[4] / 20.0,
                    spec.hi[5]]
                ]
        hi_array = [x * hi_array[0] for x in [hi_array[1] - hi_array[2]] + hi_array[3:]]                     
        if not spec.Trange[1] in coeffs:
            coeffs[spec.Trange[1]] = lo_array, hi_array
        else:
            coeffs[spec.Trange[1]] = [lo_array[i] + coeffs[spec.Trange[1]][0][i] for i in range(len(lo_array))], \
                                    [hi_array[i] + coeffs[spec.Trange[1]][1][i] for i in range(len(hi_array))]

    isFirst = True
    for T_mid in coeffs:
        # need temperature conditional for equilibrium constants
        line = '  if (T <= {:})'.format(T_mid)
        if lang in ['c', 'cuda']:
            line += ' {\n'
        elif lang == 'fortran':
            line += ' then\n'
        elif lang == 'matlab':
            line += '\n'
        file.write(line)

        lo_array, hi_array = coeffs[T_mid]
        
        if isFirst:
            line = '    Kc = '
        else:
            if lang in ['cuda', 'c']:
                line = '    Kc += '
            else:
                line = '    Kc = Kc + '
        line += ('({:.8e} + '.format(lo_array[0]) + 
                 '{:.8e} * '.format(lo_array[1]) + 
                 'logT + T * ('
                 '{:.8e} + T * ('.format(lo_array[2]) + 
                 '{:.8e} + T * ('.format(lo_array[3]) + 
                 '{:.8e} + '.format(lo_array[4]) + 
                 '{:.8e} * T))) - '.format(lo_array[5]) + 
                 '{:.8e} / T)'.format(lo_array[6]) + 
                 utils.line_end[lang]
                 )
        file.write(line)
        
        if lang in ['c', 'cuda']:
            file.write('  } else {\n')
        elif lang in ['fortran', 'matlab']:
            file.write('  else\n')
        
        if isFirst:
            line = '    Kc = '
        else:
            if lang in ['cuda', 'c']:
                line = '    Kc += '
            else:
                line = '    Kc = Kc + '
        line += ('({:.8e} + '.format(hi_array[0]) + 
                 '{:.8e} * '.format(hi_array[1]) + 
                 'logT + T * ('
                 '{:.8e} + T * ('.format(hi_array[2]) + 
                 '{:.8e} + T * ('.format(hi_array[3]) + 
                 '{:.8e} + '.format(hi_array[4]) + 
                 '{:.8e} * T))) - '.format(hi_array[5]) + 
                 '{:.8e} / T)'.format(hi_array[6]) + 
                 utils.line_end[lang]
                 )
        file.write(line)
        
        if lang in ['c', 'cuda']:
            file.write('  }\n\n')
        elif lang == 'fortran':
            file.write('  end if\n\n')
        elif lang == 'matlab':
            file.write('  end\n\n')
        isFirst = False

    line = '  Kc = '
    if sum_nu != 0:
        num = (chem.PA / chem.RU) ** sum_nu
        line += '{:.8e} * '.format(num)
    line += 'exp(Kc)' + utils.line_end[lang]
    file.write(line)

def get_nu(sp_k, rxn):
    if sp_k.name in rxn.prod and sp_k.name in rxn.reac:
        nu = (rxn.prod_nu[rxn.prod.index(sp_k.name)] -
              rxn.reac_nu[rxn.reac.index(sp_k.name)])
        # check if net production zero
        if nu == 0:
            return 0
    elif sp_k.name in rxn.prod:
        nu = rxn.prod_nu[rxn.prod.index(sp_k.name)]
    elif sp_k.name in rxn.reac:
        nu = -rxn.reac_nu[rxn.reac.index(sp_k.name)]
    else:
        # doesn't participate in reaction
        return 0
    return nu

def get_infs(rxn):
    if rxn.low:
        # unimolecular/recombination fall-off
        beta_0minf = rxn.low[1] - rxn.b
        E_0minf = rxn.low[2] - rxn.E
        k0kinf = rate.rxn_rate_const(rxn.low[0] / rxn.A,
                                     beta_0minf, E_0minf
                                     )
    elif rxn.high:
        # chem-activated bimolecular rxn
        beta_0minf = rxn.b - rxn.high[1]
        E_0minf = rxn.E - rxn.high[2]
        k0kinf = rate.rxn_rate_const(rxn.A / rxn.high[0],
                                     beta_0minf, E_0minf
                                     )
    return beta_0minf, E_0minf, k0kinf

def write_dt_comment(file, lang, rind):
    if lang in ['c', 'cuda']:
        line = '  //'
    elif lang == 'fortran':
        line = '  !'
    elif lang == 'matlab':
        line = '  %'
    line += ('partial of rxn ' + str(rind) + ' wrt T' +
             utils.line_end[lang]
             )
    file.write(line)

def write_dy_comment(file, lang, rind):
    if lang in ['c', 'cuda']:
        line = '  //'
    elif lang == 'fortran':
        line = '  !'
    elif lang == 'matlab':
        line = '  %'
    line += ('partial of rxn ' + str(rind) + ' wrt species' +
             utils.line_end[lang]
             )
    file.write(line)

def write_dy_y_finish_comment(file, lang):
    if lang in ['c', 'cuda']:
        line = '  //'
    elif lang == 'fortran':
        line = '  !'
    elif lang == 'matlab':
        line = '  %'
    line += 'Finish dYk / Yj\'s\n'
    file.write(line)

def write_dt_t_comment(file, lang, sp):
    if lang in ['c', 'cuda']:
        line = '  //'
    elif lang == 'fortran':
        line = '  !'
    elif lang == 'matlab':
        line = '  %'
    line += 'partial of dT wrt T for spec {}'.format(sp.name) + '\n'
    file.write(line)

def write_dt_y_comment(file, lang, sp):
    if lang in ['c', 'cuda']:
        line = '  //'
    elif lang == 'fortran':
        line = '  !'
    elif lang == 'matlab':
        line = '  %'
    line += 'partial of T wrt Y_{}'.format(sp.name) + '\n'
    file.write(line)

def write_rxn_params_dt(file, rxn, rev=False):
    jline = ''
    if rev:
        if (abs(rxn.rev_par[1]) > 1.0e-90
            and abs(rxn.rev_par[2]) > 1.0e-90):
            jline += ('{:.8e} + '.format(rxn.rev_par[1]) +
                      '({:.8e} / T)'.format(rxn.rev_par[2])
                      )
        elif abs(rxn.rev_par[1]) > 1.0e-90:
            jline += '{:.8e}'.format(rxn.rev_par[1])
        elif abs(rxn.rev_par[2]) > 1.0e-90:
            jline += '({:.8e} / T)'.format(rxn.rev_par[2])
        jline += '{}1.0 - '.format(' + ' if (abs(rxn.rev_par[1]) > 1.0e-90) or (
            abs(rxn.rev_par[2]) > 1.0e-90) else '')
    else:
        if (abs(rxn.b) > 1.0e-90) and (abs(rxn.E) > 1.0e-90):
            jline += '{:.8e} + ({:.8e} / T)'.format(rxn.b, rxn.E)
        elif abs(rxn.b) > 1.0e-90:
            jline += '{:.8e}'.format(rxn.b)
        elif abs(rxn.E) > 1.0e-90:
            jline += '({:.8e} / T)'.format(rxn.E)
        jline += '{}1.0 - '.format(' + ' if (abs(rxn.b)
                                             > 1.0e-90) or (abs(rxn.E) > 1.0e-90) else '')
    file.write(jline)

def write_db_dt_def(file, lang, specs, reacs, rev_reacs, dBdT_flag):
    if lang == 'cuda':
        if len(rev_reacs):
            file.write('  double dBdT[{}]'.format(len(specs)) + utils.line_end[lang])
        template = 'dBdT[{}]'
    else:
        template = 'dBdT_{}'
    for rxn in rev_reacs:
        # only reactions with no reverse Arrhenius parameters
        if rxn.rev_par:
            continue

        # all participating species
        for rxn_sp in (rxn.reac + rxn.prod):
            sp_ind = next((specs.index(s) for s in specs
                           if s.name == rxn_sp), None)

            # skip if already printed
            if dBdT_flag[sp_ind]:
                continue

            dBdT_flag[sp_ind] = True

            if lang in ['c', 'cuda']:
                dBdT = template.format(sp_ind)
            elif lang in ['fortran', 'matlab']:
                dBdT = template.format(sp_ind + 1)
            # declare dBdT
            if lang != 'cuda':
                file.write('  double ' + dBdT + utils.line_end[lang])

            # dB/dT evaluation (with temperature conditional)
            line = '  if (T <= {:})'.format(specs[sp_ind].Trange[1])
            if lang in ['c', 'cuda']:
                line += ' {\n'
            elif lang == 'fortran':
                line += ' then\n'
            elif lang == 'matlab':
                line += '\n'
            file.write(line)

            line = ('    ' + dBdT +
                    ' = ({:.8e}'.format(specs[sp_ind].lo[0] - 1.0) +
                    ' + {:.8e} / T) / T'.format(specs[sp_ind].lo[5]) +
                    ' + {:.8e} + T'.format(specs[sp_ind].lo[1] / 2.0) +
                    ' * ({:.8e}'.format(specs[sp_ind].lo[2] / 3.0) +
                    ' + T * ({:.8e}'.format(specs[sp_ind].lo[3] / 4.0) +
                    ' + {:.8e} * T))'.format(specs[sp_ind].lo[4] / 5.0) +
                    utils.line_end[lang]
                    )
            file.write(line)

            if lang in ['c', 'cuda']:
                file.write('  } else {\n')
            elif lang in ['fortran', 'matlab']:
                file.write('  else\n')

            line = ('    ' + dBdT +
                    ' = ({:.8e}'.format(specs[sp_ind].hi[0] - 1.0) +
                    ' + {:.8e} / T) / T'.format(specs[sp_ind].hi[5]) +
                    ' + {:.8e} + T'.format(specs[sp_ind].hi[1] / 2.0) +
                    ' * ({:.8e}'.format(specs[sp_ind].hi[2] / 3.0) +
                    ' + T * ({:.8e}'.format(specs[sp_ind].hi[3] / 4.0) +
                    ' + {:.8e} * T))'.format(specs[sp_ind].hi[4] / 5.0) +
                    utils.line_end[lang]
                    )
            file.write(line)

            if lang in ['c', 'cuda']:
                file.write('  }\n\n')
            elif lang == 'fortran':
                file.write('  end if\n\n')
            elif lang == 'matlab':
                file.write('  end\n\n')

def write_db_dt(file, lang, specs, rxn):
    if lang == 'cuda':
        template = 'dBdT[{}]'
    else:
        template = 'dBdT_{}'
    jline = ''
    notfirst = False
    # contribution from dBdT terms from
    # all participating species
    for sp in rxn.prod:
        if sp in rxn.reac:
            nu = (rxn.prod_nu[rxn.prod.index(sp)] -
                  rxn.reac_nu[rxn.reac.index(sp)]
                  )
        else:
            nu = rxn.prod_nu[rxn.prod.index(sp)]

        if (nu == 0):
            continue

        sp_ind = next((specs.index(s) for s in specs
                       if s.name == sp), None)

        if lang in ['c', 'cuda']:
            dBdT = template.format(sp_ind)
        elif lang in ['fortran', 'matlab']:
            dBdT = template.format(sp_ind + 1)

        if not notfirst:
            # first entry
            if nu == 1:
                jline += dBdT
            elif nu == -1:
                jline += '-' + dBdT
            else:
                jline += '{} * '.format(float(nu)) + dBdT
        else:
            # not first entry
            if nu == 1:
                jline += ' + '
            elif nu == -1:
                jline += ' - '
            else:
                if (nu > 0):
                    jline += ' + {}'.format(float(nu))
                else:
                    jline += ' - {}'.format(float(abs(nu)))
                jline += ' * '
            jline += dBdT
        notfirst = True

    for sp in rxn.reac:
        # skip species also in products, already counted
        if sp in rxn.prod:
            continue

        nu = -rxn.reac_nu[rxn.reac.index(sp)]

        sp_ind = next((specs.index(s) for s in specs
                       if s.name == sp), None)

        if lang in ['c', 'cuda']:
            dBdT = template.format(sp_ind)
        elif lang in ['fortran', 'matlab']:
            dBdT = template.format(sp_ind + 1)

        # not first entry
        if nu == 1:
            jline += ' + '
        elif nu == -1:
            jline += ' - '
        else:
            if (nu > 0):
                jline += ' + {}'.format(float(nu))
            else:
                jline += ' - {}'.format(float(abs(nu)))
            jline += ' * '
        jline += dBdT

    jline += '))'

    file.write(jline)

def write_pr(file, lang, specs, reacs, pdep_reacs, rxn, get_array, last_conc_temp=None):
    # print lines for necessary pressure-dependent variables
    line = '  {} = '.format('conc_temp' if rxn.thd_body else 'Pr')
    conc_temp_log = None
    if rxn.thd_body:
        #take care of the conc_temp collapsing
        conc_temp_log = []
    if rxn.pdep_sp:
        line += get_array(lang, 'conc', specs.index(rxn.pdep_sp))
    else:
        line += '(m'

        for thd_sp in rxn.thd_body:
            isp = specs.index(next((s for s in specs
                                    if s.name == thd_sp[0]), None))
            if thd_sp[1] > 1.0:
                line += ' + {} * '.format(thd_sp[1] - 1.0)
            elif thd_sp[1] < 1.0:
                line += ' - {} * '.format(1.0 - thd_sp[1])
            if thd_sp[1] != 1.0:
                line += get_array(lang, 'conc', isp)
                if conc_temp_log is not None:
                    conc_temp_log.append((thd_sp[0], thd_sp[1] - 1.0))
        line += ')'

        if last_conc_temp is not None:
            #need to update based on the last
            new_conc_temp = []
            for species, alpha in conc_temp_log:
                match = next((sp for sp in last_conc_temp if sp[0] == species), None)
                if match is not None:
                    coeff = alpha - match[1]
                else:
                    coeff = alpha
                if coeff != 0.0:
                    new_conc_temp.append((species, coeff))
            for species, alpha in last_conc_temp:
                match = next((sp for sp in conc_temp_log if sp[0] == species), None)
                if match is None:
                    new_conc_temp.append((species, -alpha))

            if len(new_conc_temp):
                #remake the line with the updated numbers
                line = '  conc_temp += ('

                for i, thd_sp in enumerate(new_conc_temp):
                    isp = specs.index(next((s for s in specs
                                            if s.name == thd_sp[0]), None))
                    if i > 0:
                        line += ' {}{} * '.format('- ' if thd_sp[1] < 0 else '+ ', abs(thd_sp[1]))
                    else:
                        line += '{} * '.format(thd_sp[1])
                    line += get_array(lang, 'conc', isp)
                line += ')'
            else:
                line = ''

    if rxn.thd_body and len(line):
        file.write(line + utils.line_end[lang])

    if rxn.pdep:
        if rxn.thd_body:
            line = '  Pr = conc_temp'
        beta_0minf, E_0minf, k0kinf = get_infs(rxn)
        # finish writing P_ri
        line += (' * (' + k0kinf + ')' +
                 utils.line_end[lang]
                 )
        file.write(line)

    return conc_temp_log


def write_troe(file, lang, rxn):
    line = ('  Fcent = '
            '{:.8e} * '.format(1.0 - rxn.troe_par[0]) +
            'exp(-T / {:.8e})'.format(rxn.troe_par[1]) +
            ' + {:.8e} * exp(T / '.format(rxn.troe_par[0]) +
            '{:.8e})'.format(rxn.troe_par[2])
            )
    if len(rxn.troe_par) == 4:
        line += ' + exp(-{:.8e} / T)'.format(rxn.troe_par[3])
    line += utils.line_end[lang]
    file.write(line)

    line = ('  A = log10(Pr) - 0.67 * log10(Fcent) - 0.4' +
            utils.line_end[lang]
            )
    file.write(line)

    line = ('  B = 0.806 - 1.1762 * log10(Fcent) - '
            '0.14 * log10(Pr)' +
            utils.line_end[lang]
            )
    file.write(line)

    line = ('  lnF_AB = 2.0 * log(Fcent) * '
            'A / (B * B * B * (1.0 + A * A / (B * B)) * '
            '(1.0 + A * A / (B * B)))' +
            utils.line_end[lang]
            )
    file.write(line)

def write_sri(file, lang):
    line = ('  X = 1.0 / (1.0 + log10(Pr) * log10(Pr))' +
                            utils.line_end[lang]
                            )
    file.write(line)

def write_pdep_dt(file, lang, rxn, rev_reacs, rind, pind, get_array):
    beta_0minf, E_0minf, k0kinf = get_infs(rxn)
    jline = '  j_temp = (' + get_array(lang, 'pres_mod', pind)
    jline += ' * ((' + ('-Pr' if rxn.high else '') #high -> chem-activated bimolecular rxn

    # dPr/dT
    jline += ('({:.4e} + ('.format(beta_0minf) +
              '{:.8e} / T) - 1.0) / '.format(E_0minf) +
              '(T * (1.0 + Pr)))'
              )

    if rxn.sri:
        jline += write_sri_dt(lang, rxn, beta_0minf, E_0minf, k0kinf)
    elif rxn.troe:
        jline += write_troe_dt(lang, rxn, beta_0minf, E_0minf, k0kinf)

    jline += ') * '

    if rxn.rev:
        # forward and reverse reaction rates
        jline += '(' + get_array(lang, 'fwd_rates', rind)
        jline += ' - ' + \
            get_array(lang, 'rev_rates', rev_reacs.index(rxn))
        jline += ')'
    else:
        # forward reaction rate only
        jline += '' + get_array(lang, 'fwd_rates', rind)

    jline += ' + (' + get_array(lang, 'pres_mod', pind)

    file.write(jline)

def write_sri_dt(lang, rxn, beta_0minf, E_0minf, k0kinf):
    jline = (' + X * ((('
          '{:.8} / '.format(rxn.sri[0] * rxn.sri[1]) +
          '(T * T)) * exp(-'
          '{:.8} / T) - '.format(rxn.sri[1]) +
          '{:.8e} * '.format(1.0 / rxn.sri[2]) +
          'exp(-T / {:.8})) / '.format(rxn.sri[2]) +
          '({:.8} * '.format(rxn.sri[0]) +
          'exp(-{:.8} / T) + '.format(rxn.sri[1]) +
          'exp(-T / {:.8})) - '.format(rxn.sri[2]) +
          'X * {:.8} * '.format(2.0 / math.log(10.0)) +
          'log10(Pr) * ('
          '{:.8e} + ('.format(beta_0minf) +
          '{:.8e} / T) - 1.0) * '.format(E_0minf) +
          'log({:8} * exp('.format(rxn.sri[0]) +
          '-{:.8} / T) + '.format(rxn.sri[1]) +
          'exp(-T / '
          '{:8})) / T)'.format(rxn.sri[2])
          )

    if len(rxn.sri) == 5:
        jline += ' + ({:.8} / T)'.format(rxn.sri[4])

    return jline

def write_troe_dt(lang, rxn, beta_0minf, E_0minf, k0kinf):<commit_msg>switch to 'get' methods instead of write, in order to enable plog/cheb<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) else (4 if rxn.thd 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(rxn)))
    if rxn.pdep or rxn.thd:
        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]):
        sp_real = next(spec for spec in specs if spec.name == sp)
        variable_list.append(utils.get_array('cuda', 'conc', specs.index(sp_real)))

    alphaij_count = 0
    #calculate usages
    if rxn.thd or rxn.pdep:
        fwd_usage +=1
        if rxn.rev:
            rev_usage += 1
        for i, thd in enumerate(rxn.thd_body):
            #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 sp_name in rxn.reac:
        nu = rxn.reac_nu[rxn.reac.index(sp_name)]
        if nu - 1 > 0:
            reac_usages[rxn.reac.index(sp_name)] += 1
            if rxn.thd:
                pres_mod_usage += 1
        for sp_name_2 in rxn.reac:
            if sp_name == sp_name_2:
                continue
            reac_usages[rxn.reac.index(sp_name_2)] += 1

    if rxn.rev:
        for sp_name in rxn.prod:
            nu = rxn.prod_nu[rxn.prod.index(sp_name)]
            if nu - 1 > 0:
                prod_usages[rxn.prod.index(sp_name)] += 1
            for sp_name_2 in rxn.prod:
                if sp_name == sp_name_2:
                    continue
                prod_usages[rxn.prod.index(sp_name_2)] += 1

    usages.append(fwd_usage)
    if rxn.rev:
        usages.append(rev_usage)
    if rxn.pdep or rxn.thd:
        usages.append(pres_mod_usage)
    for sp in set(rxn.reac + rxn.prod + [x[0] for x in rxn.thd_body]):
        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:
            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(rxn))
        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 = '  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 += (' - log(Fcent) * 2.0 * A * (B * '
                      '{:.6}'.format(1.0 / math.log(10.0)) +
                      ' + A * '
                      '{:.6}) / '.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 += ('- X * X * '
                      '{:.6} * '.format(2.0 / math.log(10.0)) +
                      'log10(Pr) * '
                      'log({:.4} * '.format(rxn.sri[0]) +
                      'exp(-{:4} / T) + '.format(rxn.sri[1]) +
                      'exp(-T / {:.4}))'.format(rxn.sri[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:
        jline += get_array(lang, 'pres_mod', pind)
        jline += ' * ('

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

    #get reac and prod nu sums
    for sp_name in rxn.reac:
        nu = rxn.reac_nu[rxn.reac.index(sp_name)]
        if nu == 0:
            continue
        reac_nu += nu

    if rxn.rev:
        for sp_name in rxn.prod:
            nu = rxn.prod_nu[rxn.prod.index(sp_name)]
            if nu == 0:
                continue
            prod_nu += 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(rxn))
  
    #find alphaij_hat
    alphaij_hat = None
    counter = {}
    counter[1.0] = 0
    if rxn.thd_body:
        for spec, efficiency in rxn.thd_body:
            if not efficiency in counter:
                counter[efficiency] = 0
            counter[efficiency] += 1
        counter[1.0] += (nspec - sum(counter.values()))
        alphaij_hat = max(counter.keys(), key=lambda x:counter[x])

    #now handle third body / pdep parts if needed
    if rxn.thd and not rxn.pdep:
        jline += '))'
        if alphaij_hat is not None:
            if alphaij_hat == 1.0:
                jline += ' + '
            elif alphaij_hat == -1.0:
                jline += ' - '
            else:
                jline += ' + {} * '.format(alphaij_hat) 
            jline += get_net_rate_string(lang, rxn, rind, rev_reacs, get_array)
    elif rxn.pdep:
        jline += ') + pres_mod_temp'
        jline += ')'
        if alphaij_hat is not None:
            if alphaij_hat == 1.0:
                jline += ' + '
            elif alphaij_hat == -1.0:
                jline += ' - '
            else:
                jline += ' + {} * '.format(alphaij_hat) 
            jline += '(pres_mod_temp / conc_temp)'
    else:
        jline += ')'
    file.write(jline + utils.line_end[lang])

    if rxn.pdep and not rxn.thd and not rxn.pdep_sp:
        file.write('  pres_mod_temp /= conc_temp' + utils.line_end[lang])
    return alphaij_hat

def write_rates(file, lang, rxn):
    file.write('  kf = ' + rate.rxn_rate_const(rxn.A, rxn.b, rxn.E) +
                utils.line_end[lang])
    if rxn.rev and not rxn.rev_par:
        file.write('  kr = kf / Kc' + utils.line_end[lang])
    elif rxn.rev_par:
        file.write('  kr = ' +
        rate.rxn_rate_const(rxn.rev_par[0],
                                         rxn.rev_par[1],
                                         rxn.rev_par[2]
                                         ) +
        utils.line_end[lang])
                                      
def write_dr_dy_species(lang, specs, rxn, pind, j_sp, sp_j, alphaij_hat, rind, rev_reacs, get_array):
    jline = 'j_temp'
    if rxn.pdep and not rxn.pdep_sp:
        alphaij = next((thd[1] for thd in rxn.thd_body
                        if thd[0] == sp_j.name), None)
        if alphaij is None:
            alphaij = 1.0
        if alphaij != alphaij_hat:
            diff = alphaij - alphaij_hat
            if diff != 1:
                if diff == -1:
                    jline += ' - pres_mod_temp'
                else:
                    jline += ' + {} * pres_mod_temp'.format(diff)
            else:
                jline += ' + pres_mod_temp'
    elif rxn.pdep and rxn.pdep_sp and rxn.pdep_sp == sp_j.name:
        jline += ' + pres_mod_temp / (rho * {})'.format(get_array(lang, 'y', j_sp))
    elif rxn.thd and not rxn.pdep:
        alphaij = next((thd[1] for thd in rxn.thd_body
                        if thd[0] == sp_j.name), None)
        if alphaij is None:
            alphaij = 1.0
        if alphaij != alphaij_hat:
            diff = alphaij - alphaij_hat
            if diff != 1:
                if diff == -1:
                    jline += ' - '
                else:
                    jline += ' + {} * '.format(diff)
            else:
                jline += ' + '
            jline += get_net_rate_string(lang, rxn, rind, rev_reacs, get_array)

    if (rxn.pdep or rxn.thd) and (sp_j.name in rxn.reac or (rxn.rev and sp_j.name in rxn.prod)):
        jline += ' + ' + get_array(lang, 'pres_mod', pind)
        jline += ' * '

    if sp_j.name in rxn.reac:
        if not rxn.pdep and not rxn.thd:
            jline += ' + '
        nu = rxn.reac_nu[rxn.reac.index(sp_j.name)]
        if nu != 1:
            if nu == -1:
                jline += '-'
            else:
                jline += '{} * '.format(float(nu))
        jline += 'kf'
        nu_temp = rxn.reac_nu[rxn.reac.index(sp_j.name)]
        if (nu_temp - 1) > 0:
            if isinstance(nu_temp - 1, float):
                jline += ' * pow(' + \
                    get_array(lang, 'conc', j_sp)
                jline += ', {})'.format(nu_temp - 1)
            else:
                # integer, so just use multiplication
                for i in range(nu_temp - 1):
                    jline += ' * ' + \
                        get_array(lang, 'conc', j_sp)

        # loop through remaining reactants
        for sp_reac in rxn.reac:
            if sp_reac == sp_j.name:
                continue

            nu_temp = rxn.reac_nu[rxn.reac.index(sp_reac)]
            isp = next(i for i in range(len(specs))
                       if specs[i].name == sp_reac)
            if isinstance(nu_temp, float):
                jline += ' * pow(conc' + \
                    get_array(lang, isp)
                jline += ', ' + str(nu_temp) + ')'
            else:
                # integer, so just use multiplication
                for i in range(nu_temp):
                    jline += ' * ' + \
                        get_array(lang, 'conc', isp)

    if rxn.rev and sp_j.name in rxn.prod:
        if not rxn.pdep and not rxn.thd and not sp_j.name in rxn.reac:
            jline += ' + '
        elif sp_j.name in rxn.reac:
            jline += ' + '

        nu = rxn.prod_nu[rxn.prod.index(sp_j.name)]
        if nu != -1:
            if nu == 1:
                jline += '-'
            else:
                jline += '{} * '.format(float(-1 * nu))


        jline += 'kr'
        temp_nu = rxn.prod_nu[rxn.prod.index(sp_j.name)]
        if (temp_nu - 1) > 0:
            if isinstance(temp_nu - 1, float):
                jline += ' * pow(' + \
                    get_array(lang, 'conc', j_sp)
                jline += ', {})'.format(temp_nu - 1)
            else:
                # integer, so just use multiplication
                for i in range(temp_nu - 1):
                    jline += ' * ' + \
                        get_array(lang, 'conc', j_sp)

        # loop through remaining products
        for sp_prod in rxn.prod:
            if sp_prod == sp_j.name:
                continue

            temp_nu = rxn.prod_nu[rxn.prod.index(sp_prod)]
            isp = next(i for i in range(len(specs))
                       if specs[i].name == sp_prod)
            if isinstance(temp_nu, float):
                jline += ' * pow(' + \
                    get_array(lang, 'conc', isp)
                jline += ', {})'.format(temp_nu)
            else:
                # integer, so just use multiplication
                jline += ''.join([' * ' + get_array(lang, 'conc', isp) for i in range(temp_nu)])
    
    return jline

def write_kc(file, lang, specs, rxn):
    sum_nu = 0
    coeffs = {}
    for sp_name in set(rxn.reac + rxn.prod):
        spec = next((spec for spec in specs
            if spec.name == sp_name), None)
        nu = get_nu(spec, rxn)

        if nu == 0:
            continue

        sum_nu += nu

        lo_array = [utils.round_sig(nu, 3)] + [utils.round_sig(x, 9) for x in [
                    spec.lo[6], spec.lo[0], spec.lo[0] - 1.0, spec.lo[1] / 2.0,
                    spec.lo[2] / 6.0, spec.lo[3] / 12.0, spec.lo[4] / 20.0,
                    spec.lo[5]]
                ]
        lo_array = [x * lo_array[0] for x in [lo_array[1] - lo_array[2]] + lo_array[3:]]

        hi_array = [utils.round_sig(nu, 3)] + [utils.round_sig(x, 9) for x in [
                    spec.hi[6], spec.hi[0], spec.hi[0] - 1.0, spec.hi[1] / 2.0,
                    spec.hi[2] / 6.0, spec.hi[3] / 12.0, spec.hi[4] / 20.0,
                    spec.hi[5]]
                ]
        hi_array = [x * hi_array[0] for x in [hi_array[1] - hi_array[2]] + hi_array[3:]]                     
        if not spec.Trange[1] in coeffs:
            coeffs[spec.Trange[1]] = lo_array, hi_array
        else:
            coeffs[spec.Trange[1]] = [lo_array[i] + coeffs[spec.Trange[1]][0][i] for i in range(len(lo_array))], \
                                    [hi_array[i] + coeffs[spec.Trange[1]][1][i] for i in range(len(hi_array))]

    isFirst = True
    for T_mid in coeffs:
        # need temperature conditional for equilibrium constants
        line = '  if (T <= {:})'.format(T_mid)
        if lang in ['c', 'cuda']:
            line += ' {\n'
        elif lang == 'fortran':
            line += ' then\n'
        elif lang == 'matlab':
            line += '\n'
        file.write(line)

        lo_array, hi_array = coeffs[T_mid]
        
        if isFirst:
            line = '    Kc = '
        else:
            if lang in ['cuda', 'c']:
                line = '    Kc += '
            else:
                line = '    Kc = Kc + '
        line += ('({:.8e} + '.format(lo_array[0]) + 
                 '{:.8e} * '.format(lo_array[1]) + 
                 'logT + T * ('
                 '{:.8e} + T * ('.format(lo_array[2]) + 
                 '{:.8e} + T * ('.format(lo_array[3]) + 
                 '{:.8e} + '.format(lo_array[4]) + 
                 '{:.8e} * T))) - '.format(lo_array[5]) + 
                 '{:.8e} / T)'.format(lo_array[6]) + 
                 utils.line_end[lang]
                 )
        file.write(line)
        
        if lang in ['c', 'cuda']:
            file.write('  } else {\n')
        elif lang in ['fortran', 'matlab']:
            file.write('  else\n')
        
        if isFirst:
            line = '    Kc = '
        else:
            if lang in ['cuda', 'c']:
                line = '    Kc += '
            else:
                line = '    Kc = Kc + '
        line += ('({:.8e} + '.format(hi_array[0]) + 
                 '{:.8e} * '.format(hi_array[1]) + 
                 'logT + T * ('
                 '{:.8e} + T * ('.format(hi_array[2]) + 
                 '{:.8e} + T * ('.format(hi_array[3]) + 
                 '{:.8e} + '.format(hi_array[4]) + 
                 '{:.8e} * T))) - '.format(hi_array[5]) + 
                 '{:.8e} / T)'.format(hi_array[6]) + 
                 utils.line_end[lang]
                 )
        file.write(line)
        
        if lang in ['c', 'cuda']:
            file.write('  }\n\n')
        elif lang == 'fortran':
            file.write('  end if\n\n')
        elif lang == 'matlab':
            file.write('  end\n\n')
        isFirst = False

    line = '  Kc = '
    if sum_nu != 0:
        num = (chem.PA / chem.RU) ** sum_nu
        line += '{:.8e} * '.format(num)
    line += 'exp(Kc)' + utils.line_end[lang]
    file.write(line)

def get_nu(sp_k, rxn):
    if sp_k.name in rxn.prod and sp_k.name in rxn.reac:
        nu = (rxn.prod_nu[rxn.prod.index(sp_k.name)] -
              rxn.reac_nu[rxn.reac.index(sp_k.name)])
        # check if net production zero
        if nu == 0:
            return 0
    elif sp_k.name in rxn.prod:
        nu = rxn.prod_nu[rxn.prod.index(sp_k.name)]
    elif sp_k.name in rxn.reac:
        nu = -rxn.reac_nu[rxn.reac.index(sp_k.name)]
    else:
        # doesn't participate in reaction
        return 0
    return nu

def get_infs(rxn):
    if rxn.low:
        # unimolecular/recombination fall-off
        beta_0minf = rxn.low[1] - rxn.b
        E_0minf = rxn.low[2] - rxn.E
        k0kinf = rate.rxn_rate_const(rxn.low[0] / rxn.A,
                                     beta_0minf, E_0minf
                                     )
    elif rxn.high:
        # chem-activated bimolecular rxn
        beta_0minf = rxn.b - rxn.high[1]
        E_0minf = rxn.E - rxn.high[2]
        k0kinf = rate.rxn_rate_const(rxn.A / rxn.high[0],
                                     beta_0minf, E_0minf
                                     )
    return beta_0minf, E_0minf, k0kinf

def write_dt_comment(file, lang, rind):
    if lang in ['c', 'cuda']:
        line = '  //'
    elif lang == 'fortran':
        line = '  !'
    elif lang == 'matlab':
        line = '  %'
    line += ('partial of rxn ' + str(rind) + ' wrt T' +
             utils.line_end[lang]
             )
    file.write(line)

def write_dy_comment(file, lang, rind):
    if lang in ['c', 'cuda']:
        line = '  //'
    elif lang == 'fortran':
        line = '  !'
    elif lang == 'matlab':
        line = '  %'
    line += ('partial of rxn ' + str(rind) + ' wrt species' +
             utils.line_end[lang]
             )
    file.write(line)

def write_dy_y_finish_comment(file, lang):
    if lang in ['c', 'cuda']:
        line = '  //'
    elif lang == 'fortran':
        line = '  !'
    elif lang == 'matlab':
        line = '  %'
    line += 'Finish dYk / Yj\'s\n'
    file.write(line)

def write_dt_t_comment(file, lang, sp):
    if lang in ['c', 'cuda']:
        line = '  //'
    elif lang == 'fortran':
        line = '  !'
    elif lang == 'matlab':
        line = '  %'
    line += 'partial of dT wrt T for spec {}'.format(sp.name) + '\n'
    file.write(line)

def write_dt_y_comment(file, lang, sp):
    if lang in ['c', 'cuda']:
        line = '  //'
    elif lang == 'fortran':
        line = '  !'
    elif lang == 'matlab':
        line = '  %'
    line += 'partial of T wrt Y_{}'.format(sp.name) + '\n'
    file.write(line)

def get_rxn_params_dt(rxn, rev=False):
    jline = ''
    if rev:
        if (abs(rxn.rev_par[1]) > 1.0e-90
            and abs(rxn.rev_par[2]) > 1.0e-90):
            jline += ('{:.8e} + '.format(rxn.rev_par[1]) +
                      '({:.8e} / T)'.format(rxn.rev_par[2])
                      )
        elif abs(rxn.rev_par[1]) > 1.0e-90:
            jline += '{:.8e}'.format(rxn.rev_par[1])
        elif abs(rxn.rev_par[2]) > 1.0e-90:
            jline += '({:.8e} / T)'.format(rxn.rev_par[2])
        jline += '{}1.0 - '.format(' + ' if (abs(rxn.rev_par[1]) > 1.0e-90) or (
            abs(rxn.rev_par[2]) > 1.0e-90) else '')
    else:
        if (abs(rxn.b) > 1.0e-90) and (abs(rxn.E) > 1.0e-90):
            jline += '{:.8e} + ({:.8e} / T)'.format(rxn.b, rxn.E)
        elif abs(rxn.b) > 1.0e-90:
            jline += '{:.8e}'.format(rxn.b)
        elif abs(rxn.E) > 1.0e-90:
            jline += '({:.8e} / T)'.format(rxn.E)
        jline += '{}1.0 - '.format(' + ' if (abs(rxn.b)
                                             > 1.0e-90) or (abs(rxn.E) > 1.0e-90) else '')
    return jline

def write_db_dt_def(file, lang, specs, reacs, rev_reacs, dBdT_flag):
    if lang == 'cuda':
        if len(rev_reacs):
            file.write('  double dBdT[{}]'.format(len(specs)) + utils.line_end[lang])
        template = 'dBdT[{}]'
    else:
        template = 'dBdT_{}'
    for rxn in rev_reacs:
        # only reactions with no reverse Arrhenius parameters
        if rxn.rev_par:
            continue

        # all participating species
        for rxn_sp in (rxn.reac + rxn.prod):
            sp_ind = next((specs.index(s) for s in specs
                           if s.name == rxn_sp), None)

            # skip if already printed
            if dBdT_flag[sp_ind]:
                continue

            dBdT_flag[sp_ind] = True

            if lang in ['c', 'cuda']:
                dBdT = template.format(sp_ind)
            elif lang in ['fortran', 'matlab']:
                dBdT = template.format(sp_ind + 1)
            # declare dBdT
            if lang != 'cuda':
                file.write('  double ' + dBdT + utils.line_end[lang])

            # dB/dT evaluation (with temperature conditional)
            line = '  if (T <= {:})'.format(specs[sp_ind].Trange[1])
            if lang in ['c', 'cuda']:
                line += ' {\n'
            elif lang == 'fortran':
                line += ' then\n'
            elif lang == 'matlab':
                line += '\n'
            file.write(line)

            line = ('    ' + dBdT +
                    ' = ({:.8e}'.format(specs[sp_ind].lo[0] - 1.0) +
                    ' + {:.8e} / T) / T'.format(specs[sp_ind].lo[5]) +
                    ' + {:.8e} + T'.format(specs[sp_ind].lo[1] / 2.0) +
                    ' * ({:.8e}'.format(specs[sp_ind].lo[2] / 3.0) +
                    ' + T * ({:.8e}'.format(specs[sp_ind].lo[3] / 4.0) +
                    ' + {:.8e} * T))'.format(specs[sp_ind].lo[4] / 5.0) +
                    utils.line_end[lang]
                    )
            file.write(line)

            if lang in ['c', 'cuda']:
                file.write('  } else {\n')
            elif lang in ['fortran', 'matlab']:
                file.write('  else\n')

            line = ('    ' + dBdT +
                    ' = ({:.8e}'.format(specs[sp_ind].hi[0] - 1.0) +
                    ' + {:.8e} / T) / T'.format(specs[sp_ind].hi[5]) +
                    ' + {:.8e} + T'.format(specs[sp_ind].hi[1] / 2.0) +
                    ' * ({:.8e}'.format(specs[sp_ind].hi[2] / 3.0) +
                    ' + T * ({:.8e}'.format(specs[sp_ind].hi[3] / 4.0) +
                    ' + {:.8e} * T))'.format(specs[sp_ind].hi[4] / 5.0) +
                    utils.line_end[lang]
                    )
            file.write(line)

            if lang in ['c', 'cuda']:
                file.write('  }\n\n')
            elif lang == 'fortran':
                file.write('  end if\n\n')
            elif lang == 'matlab':
                file.write('  end\n\n')

def get_db_dt(lang, specs, rxn):
    if lang == 'cuda':
        template = 'dBdT[{}]'
    else:
        template = 'dBdT_{}'
    jline = ''
    notfirst = False
    # contribution from dBdT terms from
    # all participating species
    for sp in rxn.prod:
        if sp in rxn.reac:
            nu = (rxn.prod_nu[rxn.prod.index(sp)] -
                  rxn.reac_nu[rxn.reac.index(sp)]
                  )
        else:
            nu = rxn.prod_nu[rxn.prod.index(sp)]

        if (nu == 0):
            continue

        sp_ind = next((specs.index(s) for s in specs
                       if s.name == sp), None)

        if lang in ['c', 'cuda']:
            dBdT = template.format(sp_ind)
        elif lang in ['fortran', 'matlab']:
            dBdT = template.format(sp_ind + 1)

        if not notfirst:
            # first entry
            if nu == 1:
                jline += dBdT
            elif nu == -1:
                jline += '-' + dBdT
            else:
                jline += '{} * '.format(float(nu)) + dBdT
        else:
            # not first entry
            if nu == 1:
                jline += ' + '
            elif nu == -1:
                jline += ' - '
            else:
                if (nu > 0):
                    jline += ' + {}'.format(float(nu))
                else:
                    jline += ' - {}'.format(float(abs(nu)))
                jline += ' * '
            jline += dBdT
        notfirst = True

    for sp in rxn.reac:
        # skip species also in products, already counted
        if sp in rxn.prod:
            continue

        nu = -rxn.reac_nu[rxn.reac.index(sp)]

        sp_ind = next((specs.index(s) for s in specs
                       if s.name == sp), None)

        if lang in ['c', 'cuda']:
            dBdT = template.format(sp_ind)
        elif lang in ['fortran', 'matlab']:
            dBdT = template.format(sp_ind + 1)

        # not first entry
        if nu == 1:
            jline += ' + '
        elif nu == -1:
            jline += ' - '
        else:
            if (nu > 0):
                jline += ' + {}'.format(float(nu))
            else:
                jline += ' - {}'.format(float(abs(nu)))
            jline += ' * '
        jline += dBdT

    jline += '))'

    return jline

def write_pr(file, lang, specs, reacs, pdep_reacs, rxn, get_array, last_conc_temp=None):
    # print lines for necessary pressure-dependent variables
    line = '  {} = '.format('conc_temp' if rxn.thd_body else 'Pr')
    conc_temp_log = None
    if rxn.thd_body:
        #take care of the conc_temp collapsing
        conc_temp_log = []
    if rxn.pdep_sp:
        line += get_array(lang, 'conc', specs.index(rxn.pdep_sp))
    else:
        line += '(m'

        for thd_sp in rxn.thd_body:
            isp = specs.index(next((s for s in specs
                                    if s.name == thd_sp[0]), None))
            if thd_sp[1] > 1.0:
                line += ' + {} * '.format(thd_sp[1] - 1.0)
            elif thd_sp[1] < 1.0:
                line += ' - {} * '.format(1.0 - thd_sp[1])
            if thd_sp[1] != 1.0:
                line += get_array(lang, 'conc', isp)
                if conc_temp_log is not None:
                    conc_temp_log.append((thd_sp[0], thd_sp[1] - 1.0))
        line += ')'

        if last_conc_temp is not None:
            #need to update based on the last
            new_conc_temp = []
            for species, alpha in conc_temp_log:
                match = next((sp for sp in last_conc_temp if sp[0] == species), None)
                if match is not None:
                    coeff = alpha - match[1]
                else:
                    coeff = alpha
                if coeff != 0.0:
                    new_conc_temp.append((species, coeff))
            for species, alpha in last_conc_temp:
                match = next((sp for sp in conc_temp_log if sp[0] == species), None)
                if match is None:
                    new_conc_temp.append((species, -alpha))

            if len(new_conc_temp):
                #remake the line with the updated numbers
                line = '  conc_temp += ('

                for i, thd_sp in enumerate(new_conc_temp):
                    isp = specs.index(next((s for s in specs
                                            if s.name == thd_sp[0]), None))
                    if i > 0:
                        line += ' {}{} * '.format('- ' if thd_sp[1] < 0 else '+ ', abs(thd_sp[1]))
                    else:
                        line += '{} * '.format(thd_sp[1])
                    line += get_array(lang, 'conc', isp)
                line += ')'
            else:
                line = ''

    if rxn.thd_body and len(line):
        file.write(line + utils.line_end[lang])

    if rxn.pdep:
        if rxn.thd_body:
            line = '  Pr = conc_temp'
        beta_0minf, E_0minf, k0kinf = get_infs(rxn)
        # finish writing P_ri
        line += (' * (' + k0kinf + ')' +
                 utils.line_end[lang]
                 )
        file.write(line)

    return conc_temp_log


def write_troe(file, lang, rxn):
    line = ('  Fcent = '
            '{:.8e} * '.format(1.0 - rxn.troe_par[0]) +
            'exp(-T / {:.8e})'.format(rxn.troe_par[1]) +
            ' + {:.8e} * exp(T / '.format(rxn.troe_par[0]) +
            '{:.8e})'.format(rxn.troe_par[2])
            )
    if len(rxn.troe_par) == 4:
        line += ' + exp(-{:.8e} / T)'.format(rxn.troe_par[3])
    line += utils.line_end[lang]
    file.write(line)

    line = ('  A = log10(Pr) - 0.67 * log10(Fcent) - 0.4' +
            utils.line_end[lang]
            )
    file.write(line)

    line = ('  B = 0.806 - 1.1762 * log10(Fcent) - '
            '0.14 * log10(Pr)' +
            utils.line_end[lang]
            )
    file.write(line)

    line = ('  lnF_AB = 2.0 * log(Fcent) * '
            'A / (B * B * B * (1.0 + A * A / (B * B)) * '
            '(1.0 + A * A / (B * B)))' +
            utils.line_end[lang]
            )
    file.write(line)

def write_sri(file, lang):
    line = ('  X = 1.0 / (1.0 + log10(Pr) * log10(Pr))' +
                            utils.line_end[lang]
                            )
    file.write(line)

def get_pdep_dt(lang, rxn, rev_reacs, rind, pind, get_array):
    beta_0minf, E_0minf, k0kinf = get_infs(rxn)
    jline = '  j_temp = (' + get_array(lang, 'pres_mod', pind)
    jline += ' * ((' + ('-Pr' if rxn.high else '') #high -> chem-activated bimolecular rxn

    # dPr/dT
    jline += ('({:.4e} + ('.format(beta_0minf) +
              '{:.8e} / T) - 1.0) / '.format(E_0minf) +
              '(T * (1.0 + Pr)))'
              )

    if rxn.sri:
        jline += write_sri_dt(lang, rxn, beta_0minf, E_0minf, k0kinf)
    elif rxn.troe:
        jline += write_troe_dt(lang, rxn, beta_0minf, E_0minf, k0kinf)

    jline += ') * '

    if rxn.rev:
        # forward and reverse reaction rates
        jline += '(' + get_array(lang, 'fwd_rates', rind)
        jline += ' - ' + \
            get_array(lang, 'rev_rates', rev_reacs.index(rxn))
        jline += ')'
    else:
        # forward reaction rate only
        jline += '' + get_array(lang, 'fwd_rates', rind)

    jline += ' + (' + get_array(lang, 'pres_mod', pind)

    return jline

def write_sri_dt(lang, rxn, beta_0minf, E_0minf, k0kinf):
    jline = (' + X * ((('
          '{:.8} / '.format(rxn.sri[0] * rxn.sri[1]) +
          '(T * T)) * exp(-'
          '{:.8} / T) - '.format(rxn.sri[1]) +
          '{:.8e} * '.format(1.0 / rxn.sri[2]) +
          'exp(-T / {:.8})) / '.format(rxn.sri[2]) +
          '({:.8} * '.format(rxn.sri[0]) +
          'exp(-{:.8} / T) + '.format(rxn.sri[1]) +
          'exp(-T / {:.8})) - '.format(rxn.sri[2]) +
          'X * {:.8} * '.format(2.0 / math.log(10.0)) +
          'log10(Pr) * ('
          '{:.8e} + ('.format(beta_0minf) +
          '{:.8e} / T) - 1.0) * '.format(E_0minf) +
          'log({:8} * exp('.format(rxn.sri[0]) +
          '-{:.8} / T) + '.format(rxn.sri[1]) +
          'exp(-T / '
          '{:8})) / T)'.format(rxn.sri[2])
          )

    if len(rxn.sri) == 5:
        jline += ' + ({:.8} / T)'.format(rxn.sri[4])

    return jline

def write_troe_dt(lang, rxn, beta_0minf, E_0minf, k0kinf):
{'fraction_non_alphanumeric': 0.151638174379143, 'fraction_numerical': 0.016714770528719593, 'mean_word_length': 1.3987751065124772, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 34, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5729225', 'n_tokens_mistral': 23313, 'n_tokens_neox': 22232, 'n_words': 4770}