VSOP2013 - profiling e scelta dei metodi - 2 parte: conversione di blocchi di stringhe in array

Ora che abbiamo una procedura per rileggere velocemente i file di parametri VSOP2013, dobbiamo individuare un metodo veloce per convertire le stringhe nei singoli elementi numerici costitutivi. Dalla guida dell' IMCCE e dalla lettura del file in Fortran che abbiamo tradotto in Python qualche post addietro, sappiamo che ogni riga di testo contiene 17 coefficienti numerici, che vanno moltiplicati per i corrispondenti elementi del vettore lambda, sommati insieme per formare l'argomento phi, il cui seno e coseno vengono moltiplicati per i coefficienti S e C (presenti in fondo a ogni riga), risommati insieme e moltiplicati per una potenza del tempo, espresso in millenni, trascorso dalla J2000.

Il primo problema quindi è tradurre una stringa contenente i coefficienti in singoli elementi numerici. Il file è facilmente leggibile in Fortran, molto più macchinoso in Python. Il linguaggio Python, sia nella libreria di base, sia in Numpy, ha delle funzioni di lettura da file che consentono di ricostruire una tabella dati e convertirla in array, pero' il file dell'IMCCE presenta un difetto: non è a campi delimitati, perchè i primi coefficienti interi di ogni riga sono appiccicati insieme. La lettura va quindi fatta mediante un procedimento a campi a lunghezza fissa, che in Numpy si chiama genfromtxt, che ho provato, funziona, ma mi limita nella possibilità di settare la precisione (in termini di byte utilizzati) dei termini float. Ho pertanto risolto con un approccio differente.

Dal post precedente abbiamo ottenuto un generatore tramite la funzione islice di itertools. Vorrei evitare di creare storage temporanei sotto forma di stringa, visto che alla fine dell'intero procedimento, per ogni pianeta e variabile mi aspetto un risultato singolo di tipo float (possibilmente ad alta precisione), che va sommato per tutte le potenze di 10 del tempo fino ad ottenere il parametro orbitale di interesse.

Ho pensato di procedere in questo modo:

1. creo una funzione helper convert() che ha lo scopo di spezzettare la riga nei singoli elementi (17 coefficienti + S e C (base ed esponente) e di convertirli in float
2. applico questa funzione al generatore che legge il blocco di righe nel file e ottengo, mediante una list comprehension, l'input per la funzione numpy asarray che crea una matrice di float ad alta precisione formata da tante righe quante sono le righe del blocco dati nel file e da 21 colonne.
3. restituisco la matrice alla funzione chiamante

Allego di seguito un po' di codice:

Inizializzazione

import cPickle
import numpy as np
from itertools import islice
from fastnumbers import fast_int, fast_float
import struct
import cProfile, pstats 

limits = cPickle.load(open('VSOP2013_ranges.pickle','rb'))

fmt = struct.Struct("""6s 3s 3s 3s 3s x 3s 3s 3s 3s 3s x 4s 4s 4s 4s x
                       6s x 3s 3s 3s 20s x 3s 20s x 3s x""")

Parte di questi import è già stato visto nei post precedenti. Fastnumbers è una libreria, presente nel repository PyPI, che consente la creazione veloce di integer e float a partire dalle stringhe, sostituendo in modo più efficiente le funzioni di casting int() e float(). Si installa facilmente con pip (sudo pip install fastnumbers).

Di cPickle abbiamo già parlato a sufficienza.

Ho definito una stringa di formato fmt da utilizzare con il metodo struct per spezzettare la lista, con il metodo indicato di seguito.

Funzione helper convert()

def convert(i):
    
    a = fmt.unpack(i)
    return [fast_float(a[j]) for j in xrange(1,22)]

Il modulo di libreria struct serve a gestire la conversione tra valori Python e strutture C rappresentate come stringhe. E' molto raffinato, qui io lo uso in modo elementare per compilare una sequenza di delimitatori di stringa che mi servono per spezzare la riga in singoli elementi. Si traduce facilmente così: nella stringa fmt s sta per carattere e il numero che precede la s indica il numero di caratteri da estrarre, mentre la x rappresenta uno spazio.

Una volta fatto l'unpack della stringa passata alla funzione convert, ho una lista di 17 valori float distinti e 2 basi con esponenti separati, che restituisco alla funzione chiamante. Scarto il numeratore di riga, che non è più necessario. La precisione dei float è quella base di Python (64 bit), sufficiente per la fase di lettura.

La nuova funzione slice_file()

def slice_file(planet, var, power, precision):
    
    file_in = open("VSOP2013p{}.dat".format(planet))
    start, end = limits[planet][var][power][precision]
    data = np.asarray([convert(i) for i in islice(file_in, start, end+1)], np.longdouble)
    file_in.close()
    coeffs = data[:,:17]
    S = 10**data[:,18:19]*data[:,17:18]
    C = 10**data[:,20:21]*data[:,19:20]
    return (coeffs, S, C)

la funzione acquisisce come parametri il codice numerico del pianeta, la variabile corrispondente al parametro orbitale, la potenza di 10 e la precisione richiesta. Apertura e chiusura del file sono poco intensivi, per cui non li considero ai fini del timing. la variabile data corrisponde ad un array di numpy contenente i valori elementari delle variabili. Una volta creato l'array lo divido in una matrice di n per 17 elementi (coefficienti), un vettore colonna S e un vettore colonna C di n elementi, ottenuti con funzioni veloci di numpy, e li restituisco alla funzione chiamante, che in questo caso è il profiler.

Profiling

Per finire eseguo il timing della funzione islice con il metodo cProfile.run che abbiamo già usato in un post precedente:

if __name__ == '__main__':
    profile = cProfile.run("slice_file(1,1,0,12)","prof.prof")
    p = pstats.Stats("prof.prof")
    p.sort_stats("cumulative").print_stats(10)

Ed ecco i risultati con il mio computer:

Python 2.7.9 (default, Apr  2 2015, 15:33:21) 
[GCC 4.9.2] on linux2
Type "copyright", "credits" or "license()" for more information.
>>> ================================ RESTART ================================
>>> 
Mon May 11 16:47:37 2015    prof.prof

         36279 function calls in 0.023 seconds

   Ordered by: cumulative time
   List reduced from 11 to 10 due to restriction <10>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.023    0.023 <string>:1(<module>)
        1    0.004    0.004    0.023    0.023 /home/ubuntu/Scrivania/vsop2013/solution/VSOP2013_load_range7.py:18(slice_file)
     1577    0.009    0.000    0.017    0.000 /home/ubuntu/Scrivania/vsop2013/solution/VSOP2013_load_range7.py:13(convert)
    33117    0.007    0.000    0.007    0.000 {fastnumbers.fast_float}
        1    0.000    0.000    0.003    0.003 /usr/lib/python2.7/dist-packages/numpy/core/numeric.py:392(asarray)
        1    0.003    0.003    0.003    0.003 {numpy.core.multiarray.array}
     1577    0.001    0.000    0.001    0.000 {method 'unpack' of 'Struct' objects}
        1    0.000    0.000    0.000    0.000 {method 'close' of 'file' objects}
        1    0.000    0.000    0.000    0.000 {open}
        1    0.000    0.000    0.000    0.000 {method 'format' of 'str' objects}


>>> 

La routine, per una precisione richiesta di 1e-12, appare abbastanza veloce, considerato che genera le matrici in circa 2 centesimi di secondo. Per avere una reale stima dei tempi esecutivi, è necessario ora inserirla nel contesto più generale del programma di calcolo dei parametri per ogni pianeta. Di questo cominciamo ad occuparci nel prossimo post.