L'economia nell'uso delle risorse è importante in qualsiasi ambiente di sviluppo, ma diventa fondamentale quando si usa Python. La sua lentezza operativa puo', nella maggior parte dei casi, essere attributa alla scorretta scelta di tools, alla scrittura di codice non preottimizzato, magari perchè si è di vecchia scuola come me, abituato ad usare il Visual Basic o il C++ o Java, linguaggi compilati e naturalmente più performanti. Python ha pero' il vantaggio di essere un linguaggio di altissimo livello, consentendo di scrivere programmi non solo più puliti e in modo più veloce, ma anche di costruire in pochissime righe strutture che, in linguaggi più classici, richiederebbero molto più codice.
Se torniamo alle prove del post precedente, vediamo che ho associato la ricerca di migliori performance alla scelta di librerie e tools differenti. Ho sostituito mpmath con gmpy2, passando da una libreria in pure Python ad una in C compilato. Ora pero' devo iniziare a intervenire sul codice, per gestire meglio l'accesso ai file e le routine di calcolo.
Se guardate i file di parametri che sono presenti nel sito ftp dell'IMCCE potete verificare che si tratta di nove file, denominati VSOP2013pX.dat, dove la X sta per il numero corrispondente al pianeta, 1 per Mercurio, 2 per Venere, 3 per il baricentro Terra-Luna eccetera. Ogni file è strutturato in questo modo:
- un header di una riga che contiene alcuni numeri interi che rappresentano il pianeta, la variabile (di 6) che corrisponde ad un parametro orbitale, la potenza di a cui va elevato il tempo in millenni dalla data JD2000 di riferimento, il numero di righe successive da leggere e una serie di informazioni testuali che corrispondono agli integer iniziali.
- una serie di righe testuali che contengono i termini da utilizzare per il calcolo in relazione al pianeta scelto, alla variabile e alla potenza da applicare al tempo, contenuti, come si è detto, nell'header.
Ho detto in un post precedente che per evitare di passare le ore ad aspettare la fine del calcolo, si puo' predefinire il livello di precisione, riducendo il numero di termini da utilizzare. Ogni riga consente di stimare tale livello utilizzando gli ultimi quattro elementi della riga, che sono i parametri S e C e relativi esponenti in base 10 per la notazione scientifica. Da notare che questa verifica, indicata dagli autori dell'IMCCE, non è dipendente dal tempo, per cui è possibile effettuarla una volta per tutte, visto che fare in tempo reale calcolo e verifica di precisione è sicuramente dispendioso. In questo post cerco di scrivere un piccolo programma che fa il calcolo una volta per tutte e salva in un file esterno i risultati.
Iniziamo con l'intestazione e due parole di spiegazione.
import gmpy2 as gmp gmp.get_context().precision=200 from itertools import islice import cPickle
In queste prime righe importo alcuni moduli fondamentali:
- gmpy2, libreria di alta precisione molto veloce,l'abbiamo già discussa nel post precedente. La precisione è fissata arbitrariamente in 200 bits, contro i 53 standard di Python
- itertools è un modulo della libreria standard. Contiene funzioni per gestire gli iterators, cioè oggetti che rappresentano uno stream di dati. Sono utili per scorrere sequenzialmente un dato che puo' essere una lista, una tupla, un dictionary o anche, come vediamo in seguito, un file su disco. Delle funzioni disponibili mi interessa, in questo contesto, solo la funzione islice, che permette di generare un sottoinsieme da un iterator.
- Infine cPickle è un modulo che consente di serializzare (cioè trasformare in stream di byte) un oggetto di Python e anche deserializzarlo, cioè fare l'operazione inversa. Il valore di questo modulo è inestimabile, perchè se l'oggetto ha una struttura complessa e non risolvibile immediatamente in una sequenza, Python lo fa per noi, consentendo di salvare su file qualcosa che possiamo rileggere in altro contesto e ritrovare identico a come è stato costruito.
La prima cosa che voglio fare è trovare i puntatori agli header nei file di parametri. Leggiamo insieme il codice seguente, parte iniziale di una funzione (def):
def find_limits(): # create an empty dictionary limits={} for planet in xrange(1,10): # read date files, one per planet file_in = open('VSOP2013p'+str(planet)+'.dat') # create a list with the header position and values for variables planet, var, power elenco = [(s[0], s[1].split()) for s in enumerate(file_in) if "VSOP" in s[1]] file_in.close()
Uso un dizionario strutturato (limits) per conservare i riferimenti ai blocchi di ogni file che sono riferiti al pianeta, al parametro orbitale (variabile), alla potenza del tempo e alla soglia di precisione desiderata.
Inizio con un dizionario vuoto {}, che estendero' progressivamente di livello. Preliminarmente, pero', devo individuare gli headers nei file originari dell'IMCCE.
Le righe successive mostrano come aprire i file di parametri e come ricavarne gli header. La variabile elenco è pensata per elencare le sole righe del file che contengono la sequenza di testo 'VSOP2013'. Per fare questo uso una particolare struttura di Python chiamata list comprehension che rientra nella programmazione di tipo funzionale. In pratica la variabile file_in, riferita al file aperto all'inizio, viene scandita con la funzione enumerate, che restituisce insieme la posizione in lista e l'elemento della lista (in questo caso le righe del file) e, usando come filtro il termine 'VSOP2013' si accumulano nella lista elenco tutte le righe che soddisfano la condizione. Rispetto al ciclo for usuale, strutturato in più righe, la lista comprehension risulta più compatta come spazio e molto più performante in termini di velocità esecutiva.
# extend the dictionary with a subdictionary for each planet limits[planet]={} for var in xrange(1,7): # extend the dictionary with a subdictionary for each variable limits[planet][var] = {} for power in xrange(0,21): # extend the dictionary with a subdictionary for each power limits[planet][var][power] = {} for threshold in xrange(0,21): # extend the dictionary with a two values tuple (min and max) # for each sensibility threshold # and for each item in list assign (0,0) to the tuple limits[planet][var][power][threshold]=(0,0)
Nelle righe successive costruiamo il nostro dizionario, estendonolo progressivamente di livello. Usando un ciclo for nidificato creo tanti sottodizionari quante sono le modalità per ogni pianeta, ogni variabile, ogni potenza e ogni soglia di precisione. Al termine inizializzo a (0,0) il dettaglio più fine del dizionario.
A questo punto dispongo di due strumenti fondamentali, il dizionario contenitore e l'elenco degli header (per ogni file di parametri il dizionario non viene rigenerato, bensì ampliato. Andiamo oltre con il codice.
# create a tuple with range limits for each planet, variable, power for i in xrange(len(elenco)): start = int(elenco[i][0])+1 planet = int(elenco[i][1][1]) var = int(elenco[i][1][2]) power = int(elenco[i][1][3]) end = int(elenco[i][1][4])+start # load a file chunk for each voice in list and find extreme value for each threshold file_in = open('VSOP2013p'+str(planet)+'.dat') segmento = islice(file_in, start, end) for j, k in enumerate(segmento): cursor = j line = k.split() Sb = float(line[-4]) Se = int(line[-3]) S = Sb * 10 ** Se Cb =float(line[-2]) Ce = int(line[-1]) C = Cb * 10 ** Ce ro = gmp.sqrt(S*S+C*C) threshold = -int(gmp.floor(gmp.log10(ro))) limits[planet][var][power][threshold]=(start, start+cursor) file_in.close()
Inizio una scansione sequenziale della lista in <elenco>, da cui, riga per riga, ottengo il punto di inizio e fine nonchè le variabili e i parametri di riferimento nel file. Una volta ottenuti i riferimenti di inizio e fine, li uso con la funzione islice di itertools per ottenere un iteratore che è un sottoinsieme del file. Da ogni riga ricavo i parametri S e C e ottengo il livello di precisione corrispondente (dato dalla radice quadrata della somma dei quadrati di S e C). Trasformo il valore ottenuto nel logaritmo in base 10, cambio il segno e trasformo il risultato in integer. Questo è il livello di precisione che mi serve per indicizzare il mio dizionario. In corrispondenza della chiave pianeta, variabile, potenza e soglia (indicata dall'integer ottenuto dalla procedura precedente) indico i nuovi estremi di posizione nel file. A questo punto il lavoro è quasi concluso.
Ultime righe di codice:
for planet in xrange(1,10): for var in xrange(1,7): for power in xrange(0,21): for threshold in xrange(1,21): prec_ = limits[planet][var][power][threshold-1] subs_ = limits[planet][var][power][threshold] if subs_ == (0,0) and prec_ !=0: limits[planet][var][power][threshold] = limits[planet][var][power][threshold-1] return limits
Un'ultima procedura prima di restituire il dizionario al termine della funzione consiste nel riempire i buchi: se non è stato assegnato un valore alla tupla nello specifico dettaglio, la tupla rimane (0,0), il che mi permetterà, nell'uso successivo, di considerare inutile il dettaglio per la precisione richiesta. Per le soglie di precisione superiori non toccate dalla procedura di riempimento, ho scelto di assumere il valore di precisione inferiore se diverso da (0,0).
Vediamo ora la parte di salvataggio e caricamento del dizionario.
if __name__ == '__main__': limits = find_limits() cPickle.dump(limits, open('VSOP2013_ranges.pickle','wb')) # reload limits from disk load_limits = cPickle.load(open('VSOP2013_ranges.pickle','rb'))
La funzione dump di cPickle permette di serializzare il nostro dizionario e salvarlo su disco. La funzione load effettua il caricamento e la deserializzazione. In questo modo, dopo circa un minuto di lavoro (nella mia macchina, almeno), avremo un file illeggibile in modo diretto che potremo rileggere in programmi successivi per ricreare lo stesso dizionario, con tempi di elaborazione molto inferiori.
Salvo in github il programma sotto il nome di VSOP2013_ranges.py, nonchè il file serializzato che contiene il dizionario.Nel prossimo post vedremo come utilizzare in pratica il dizionario e misureremo le performance del nuovo codice. Al prossimo post.
