<commit_before>#!/usr/bin/python
# -*- coding: utf-8 -*-

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals

from math import exp, log

import numpy as np
from scipy.special import gammaln
import sys

def generate_basis_sets(n=1001, nbf=500, verbose=True):
    """ 
    Generate the basis set for the BASEX method. 

    This function was adapted from the a matlab script BASIS2.m,
    with some optimizations.

    Parameters:
    -----------
      n : integer : problem size ?
      nbf: integer: number of basis functions ?

    Returns:
    --------
      M, Mc : np.matrix
    """
    if n % 2 == 0:
        raise ValueError('The n parameter must be odd (more or less sure about it).')


    Rm = n//2 + 1

    I = np.arange(1,n+1)

    R2 = (I - Rm)**2
    # R = I - Rm
    M = np.zeros((n, nbf))
    Mc = np.zeros((n, nbf))

    M[:,0] = 2*np.exp(-R2)
    Mc[:,0] = np.exp(-R2)

    gammaln_0o5 = gammaln(0.5) 

    if verbose:
        print('Generating BASEX basis sets for n = {}, nbf = {} (this can take a while):\n'.format(n, nbf))
        sys.stdout.write('0')
        sys.stdout.flush()


    for k in range(1, nbf):
        k2 = k*k # so we don't recalculate it all the time
        log_k2 = log(k2) 
        angn = exp(
                    k2-2*k2*log(k) +
                    #np.log(np.arange(0.5, k2 + 0.5)).sum() # original version
                    gammaln(k2 + 0.5) - gammaln_0o5  # optimized version
                    )
        M[Rm-1, k] =  2*angn

        for l in range(1, n-Rm+1):
            l2 = l*l
            log_l2 = log(l2)
            # this was factorised from the argument of the exp in the
            # 3rd inner loop (allow to reuse the computed value 
            aux_factor = k2 - l2 - k2*log_k2  + gammaln(k2+1) - gammaln_0o5

            val = exp(k2 - l2 + 2*k2*log((1.0*l)/k))
            Mc[l-1+Rm, k] = val
            Mc[Rm-l-1, k] = val

            aux = val + angn*Mc[l+Rm-1, 0]
            #
            # A. version with a loop that we can completly vectorize 
            #
            #for p in range(max(1, l2 - 100), min(k2 - 1,  l2 + 100)+1):

            #    aux += exp(
            #            #k2 - l2 - k2*log_k2  # original 
            #            aux_factor            # factor containting elements computed outside of the loop
            #            + p*log_l2
            #            # version 1 : 
            #            #
            #            # np.log(np.arange(p+1, k**2+1)).sum() + \
            #            # np.log(np.arange(0.5, k**2 - p)).sum() - \
            #            # np.log(np.arange(1, k**2 - p + 1)).sum() 
            #            #
            #            # version 2 : (optimized) 
            #            #
            #            # we use here the fact that 
            #            # np.log(np.arange(p, k)).sum() == gammaln(k) - gammaln(p)
            #            # and put as much elements as possible out of this loop
            #            - gammaln(p+1)            # + gammaln(k2+1)   # factorized in  aux_factor
            #            + gammaln(k2 - p + 0.5)     # - gammaln_0o5   # factorized in  aux_factor  
            #            - gammaln(k2 - p + 1)       # since gammaln(1) == 0
            #            )
            # 
            # B. Vectorized version

            p = np.arange(max(1, l2 - 100), min(k2 - 1,  l2 + 100)+1)

            # see the non vecotorized version above for the origin of this expression
            aux += np.exp(aux_factor + p*log_l2
                      - gammaln(p+1) + gammaln(k2 - p + 0.5) 
                      - gammaln(k2 - p + 1)).sum()

            # End of vectorized version

            aux *= 2

            M[l+Rm-1, k] = aux
            M[Rm-l-1, k] = aux

        if verbose and k % 50 == 0:
            sys.stdout.write('...{}'.format(k))
            sys.stdout.flush()

    if verbose:
        print("...{}".format(k+1))


    return M.view(np.matrix), Mc.view(np.matrix)
<commit_msg>Cleaned up the code of `BASEX/basis.py`<commit_after>#!/usr/bin/python
# -*- coding: utf-8 -*-

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals

from math import exp, log

import numpy as np
from scipy.special import gammaln
import sys

def generate_basis_sets(n=1001, nbf=500, verbose=True):
    """ 
    Generate the basis set for the BASEX method. 

    This function was adapted from the a matlab script BASIS2.m,
    with some optimizations.

    Parameters:
    -----------
      n : integer : problem size ?
      nbf: integer: number of basis functions ?

    Returns:
    --------
      M, Mc : np.matrix
    """
    if n % 2 == 0:
        raise ValueError('The n parameter must be odd (more or less sure about it).')


    Rm = n//2 + 1

    I = np.arange(1,n+1)

    R2 = (I - Rm)**2
    # R = I - Rm
    M = np.zeros((n, nbf))
    Mc = np.zeros((n, nbf))

    M[:,0] = 2*np.exp(-R2)
    Mc[:,0] = np.exp(-R2)

    gammaln_0o5 = gammaln(0.5) 

    if verbose:
        print('Generating BASEX basis sets for n = {}, nbf = {}:\n'.format(n, nbf))
        sys.stdout.write('0')
        sys.stdout.flush()


    for k in range(1, nbf):
        k2 = k*k # so we don't recalculate it all the time
        log_k2 = log(k2) 
        angn = exp(
                    k2-2*k2*log(k) +
                    #np.log(np.arange(0.5, k2 + 0.5)).sum() # original version
                    gammaln(k2 + 0.5) - gammaln_0o5  # optimized version
                    )
        M[Rm-1, k] =  2*angn

        for l in range(1, n-Rm+1):
            l2 = l*l
            log_l2 = log(l2)

            val = exp(k2 - l2 + 2*k2*log((1.0*l)/k))
            Mc[l-1+Rm, k] = val
            Mc[Rm-l-1, k] = val

            aux = val + angn*Mc[l+Rm-1, 0]

            p = np.arange(max(1, l2 - 100), min(k2 - 1,  l2 + 100)+1)

            # We use here the fact that for p, k real and positive
            #
            #  np.log(np.arange(p, k)).sum() == gammaln(k) - gammaln(p) 
            #
            # where gammaln is scipy.misc.gammaln (i.e. the log of the Gamma function)
            #
            # The following line corresponds to the vectorized third
            # loop of the original BASIS2.m matlab file.

            aux += np.exp(k2 - l2 - k2*log_k2 + p*log_l2
                      + gammaln(k2+1) - gammaln(p+1) 
                      + gammaln(k2 - p + 0.5) - gammaln_0o5
                      - gammaln(k2 - p + 1)).sum()

            # End of vectorized third loop

            aux *= 2

            M[l+Rm-1, k] = aux
            M[Rm-l-1, k] = aux

        if verbose and k % 50 == 0:
            sys.stdout.write('...{}'.format(k))
            sys.stdout.flush()

    if verbose:
        print("...{}".format(k+1))


    return M.view(np.matrix), Mc.view(np.matrix)

{'fraction_non_alphanumeric': 0.13388096267473698, 'fraction_numerical': 0.03343421242253927, 'mean_word_length': 1.3997233748271092, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7006540', 'n_tokens_mistral': 2576, 'n_tokens_neox': 2404, 'n_words': 640}