Tuesday, May 19, 2015

Intermezzo su list comprehensions e generators

Ho iniziato a usare Python circa 3 anni fa, dopo aver passato anni a scrivere codice in C++, Java e Visual Basic. All'inizio scrivevo come ero stato sempre abituato, con cicli for e while, scrivendo funzioni e un po' di programmazione per oggetti. Python mi ha fatto conoscere tipologie dai dati che prima non usavo, come le liste, le tuple, i dizionari e i set. C'è voluto altro tempo per capire che Python è un linguaggio con una ricca scelta di strumenti di programmazione che, grazie alla necessità di velocizzare l'esecuzione di alcune tipologie di software, mi ha costretto a considerare scelte differenti da quelle che caratterizzavano lo stile precedente. Se inizi a lavorare con le list comprehensions e i generatori ci vuole comunque un po' di tempo per capire che le puoi usare dappertutto, al posto di cicli for e minimizzando lo storage intermedio di dati. Oggi faccio una pausa nell'elaborazione del software per il calcolo astronomico e parlo un po' di questi strumenti.

Gli esempi che si trovano nei manuali sono utili per afferrare l'idea, ma sono troppo semplici per essere esaustivi delle possibilità che i metodi funzionali offrono, per cui ho pensato di scrivere un programmino d'esempio che complica un po' le cose ma spero faccia capire meglio come si puo' procedere.

Se aprite una console Python e scrivete:

import this

viene fuori un Easter Egg (Uovo di Pasqua), che è il seguente (per curiosità è la PEP 20):

The Zen of Python, by Tim Peters

Beautiful is better than ugly.
Explicit is better than implicit.
Simple is better than complex.
Complex is better than complicated.
Flat is better than nested.
Sparse is better than dense.
Readability counts.
Special cases aren't special enough to break the rules.
Although practicality beats purity.
Errors should never pass silently.
Unless explicitly silenced.
In the face of ambiguity, refuse the temptation to guess.
There should be one-- and preferably only one --obvious way to do it.
Although that way may not be obvious at first unless you're Dutch.
Now is better than never.
Although never is often better than *right* now.
If the implementation is hard to explain, it's a bad idea.
If the implementation is easy to explain, it may be a good idea.
Namespaces are one honking great idea -- let's do more of those!

Nella comunità Python questa piccola poesia è nota e citata quanto i Monty Python (che hanno dato il nome al linguaggio), e trovate spesso piccole citazioni che vanno a braccetto con spam e eggs e altre divertenti citazioni dei comici inglesi. Mi interessa farvi notare l'enfasi che viene data alla leggibilità del codice, alla sua semplicità costruttiva, e alla ricerca dell'unica soluzione che risulta ovvia (potremmo dire pythonica).

Il piccolo esercizio che segue discende da un software che ho scritto per gestire la mia spesa domestica, visto che la mia banca online mi permette di esportare in formato Excel un batch di operazioni bancarie. Python ha, ovviamente, delle belle librerie di terze parti per leggere file nel formato xls di Excel, ma preferisco, per vecchia abitudine, usare i file csv, che hanno il pregio della leggibilità diretta come testi.

Dovendo assemblare vari pezzi (il sito della mia banca non mi permette di estrarre tutto l'archivio storico in un colpo solo, per cui sono costretto a unire degli estratti parziali), mi sono posto il problema di come eliminare tutte le linee di intestazione e commento e ridurmi alle sole linee che contengono le operazioni bancarie.

Per avere un codice più essenziale ho usato estesamente le list comprehension e i generatori, come faccio sempre allego il codice poche righe alla volta e le commento,

from datetime import date, timedelta
import re, string
from random import randrange, shuffle

# helper functions:
# rand_string() generates a random shuffled string from all ascii chars
# europ_date converts a date in a string conforming to european format

def rand_shuffled_string():
    all_letters = list(string.ascii_letters)
    shuffle(all_letters)
    return "".join(all_letters)

def european_date(y,m,d, diff):
    return (date(y,m,d) + timedelta(diff)).strftime("%d/%m/%Y")

In breve: importo alcuni moduli essenziali: date e timedelta da datetime per la gestione della data, re per le regular expressions che uso più avanti, string per generare dei set di caratteri ascii dalle funzioni built-in e randrange e shuffle dal modulo random, per generare delle stringhe pseudocasuali.

Le due funzioni helper, banalmente, servono a creare una stringa pseudocasuale e una stringa per la data in formato europeo.

# To create a text file on hard disk we use a list generator
# to generate a text file where each row contains a date and a string

file_out = open("prova_gen.dat","w")
[file_out.write(
    "#".join(
        (european_date(2000,1,1,t),
         rand_shuffled_string(),
         '\n')
        )
    )
    if randrange(0,100) < 90
    else file_out.write("<----------------------- header\n")
    for t in xrange(0,3650,15)
]
file_out.close()

La list comprehension è formata da una funzione applicata ad una o più variabili che entrano in un ciclo e puo' comprendere un filtro di tipo if ... else. Nel caso specifico, applico la funzione helper european_date() alla variabile t che rappresenta lo scarto temporale tra la data di inizio e l'estremo del ciclo for. In pratica creo una data in formato europeo ogni quindici giorni a partire dal 1 gennaio 2000, la associo ad una stringa pseudocasuale e salvo la riga su disco se un numero causale tra 0 e 100 rimane compreso tra 0 e 90, quindi per il 90% dei casi, altrimenti scrivo sul file una pseudo riga di intestazione.

E' interessante notare che l'inclusione di funzioni esterne da me create, quindi presenti nel namespace, mi permette di creare list comprehension molto elaborate ma ancora molto leggibili.

# reload the file just written, filtering rows not begininng with a date
pattern = '\d{2}\/\d{2}\/\d{4}'
less_ = (i for i in open("prova_gen.dat") if re.match(pattern, i) != None)
all_set = (i.split("#") for i in less_)

# generates a list of tuples with date and event 
events = [(i[0], i[1]) for i in all_set]
# prints the first 10 tuples
print events[:10]

Nelle ultimo righe faccio il lavoro di rilettura del file e di filtraggio delle righe. Vorrei farvi notare come la funzione open(<file>) funziona come un generatore, mettendo a disposizione, uno alla volta, i contenuti del file in forma di righe singole. Viene ciclato tutto il file per estrarre le sole righe che contengono una data (a questo scopo uso l'espressione regolare del pattern, che corriponde ad una definita sequenza di caratteri). Il tutto, anziché come in precedenza, in una lista comprehension, in un generatore di espressioni, che quindi non genera una lista e, di fatto, non impegna memoria. Se provo a stampare la variabile less_ l'interprete Python mi risponde <generator object at 0x7f2f39985b90>, quindi in questa fase non ho ancora nulla di leggibile. Idem per la variabile all_set, che splitta ogni singola riga usando il separatore # che ho definito in fase di creazione del file

Infine la variabile event è generata attraverso una list comprehension, quindi è una lista, contenente tuple di due elementi (data e stringa pseudocasuale), che posso stampare o utilizzare altrimenti come insieme di dati. Si noti che, in tutto questo procedimento, è solo alla fine che viene generata una variabile che impegna memoria. Non ho dovuto infatti usare storage intermedi, pur avendo scomposto il procedimento in più fasi. In questo modo spero di avervi mostrato un po' di programmazione funzionale e di lazy evaluation (tipico dei generatori). Si dice lazy (pigro) perchè non pretende di buttare fuori l'output tutto insieme, ma fornisce un risultato alla volta quando richiesto da un'altra funzione o istruzione.

Il prossimo post ritorna sulle VSOP2013 e sull'elaborazione delle matrici di dati che abbiamo reso disponibili. Sarà un'occasione per riprendere il discorso sull'uso dei generatori. Come vedremo, sono anche una delle possibilità tecniche offerte per la velocizzazione del codice.

Monday, May 11, 2015

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.

Wednesday, May 6, 2015

VSOP2013 - profiling e scelta dei metodi - 1 parte: rilettura dei file di parametri orbitali

La convinzione diffusa che Python sia un ambiente di sviluppo caratterizzato da lentezza esecutiva, almeno in rapporto ai compilatori tradizionalmente ritenuti più efficienti, come C, C++ e Fortran, è fondata solo in parte.

Giocano sicuramente a sfavore alcune caratteristiche intrinseche del linguaggio, che è fondamentamentalmente un linguaggio di scripting, a tipizzazione dinamica, con alcune strutture dati, quali le liste, che impongono all'interprete grandi spostamenti in memoria, una volta superate certe dimensioni. Per esempio, aggiornare una lista inserendo elementi al suo interno, anzichè accodarli in append, puo' risultare molto pesante quando la lista è costituita da migliaia di termini, dovendosi procedere, ad ogni nuovo inserimento, allo spostamento globale in memoria di tutta la lista. Non è detto pero' che la libreria di base non fornisca metodi utili a evitare sia overhead che massiccio impegno di memoria. Per ogni necessità, possiamo confrontare tra loro più metodi, utilizzando allo scopo degli appositi strumenti di misura, che vanno sotto il nome di profiler.

Oggi continuo il lavoro iniziato nel post precedente, in cui abbiamo ottenuto un dizionario contenente riferimenti ai limiti di inizio e fine di ogni blocco dati all'interno dei file di parametri VSOP2013.

La creazione del dizionario e il suo storage, anche se lenti, sono stati un'operazione da fare una volta sola, della durata di circa un minuto, per cui non mi sono preoccupato particolarmente della velocità esecutiva.

Diverso è il problema della rilettura ed elaborazione dei dati, al fine di ricavare i parametri da sottoporre poi a calcolo per ottenere, alla fine, le variazioni temporali da applicare ai parametri orbitali per gli otto pianeti più il baricentro terra-luna. Non è un compito semplice, voglio arrivare alla fine dell'estrazione ad avere una o più matrici numeriche con cui effettuare i calcoli con i metodi dell'algebra matriciale e, in particolare, con la libreria numpy, le cui routine critiche sono scritte in C e Fortran.

Ho scomposto il lavoro di definizione in più tempi e in più alternative. Oggi provo a scrivere la routine di lettura usando diversi metodi, per confrontarli poi con le funzioni del modulo di misura integrato nella libreria di Python cProfile, che naturalmente dovro' importare nel mio codice, se voglio effettuare le misure dall'interno del codice stesso. Per ovviare alla velocità (comunque elevata) di esecuzione, utilizzo il modulo timeit per ognuna delle tre funzioni

Il listato allegato di seguito, che contiene l'intero programma, usa tre metodi diversi per leggere i file dati (ho scelto quello relativo al pianeta Mercurio, alla sola prima variabile, al primo blocco e a una soglia di arresto pari a 1e-20, praticamente l'intera serie). In questa fase non faccio trattamento delle stringhe per ricavarne i dati numerici, mi accontento di verificare eventuali differenze nella velocità di recupero delle stesse. In coda al listato si trovano una def main() che riepiloga le tre routine e, infine, i metodi di profiling e il timing con il metodo timeit con 100 ripetizioni per ogni funzione.

import cPickle
from itertools import islice
import cProfile, pstats, timeit

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

def timing01():
    chunk = []
    file_in = open('VSOP2013p1.dat')
    all_lines = file_in.readlines()
    start, end = limits[1][1][0][20]
    chunk = all_lines[start:end+1]
    file_in.close()
    return chunk
    
def timing02():
    chunk=[]
    file_in = open('VSOP2013p1.dat')
    start, end = limits[1][1][0][20]
    chunk = [line for line in islice(file_in, start, end+1)]
    file_in.close()    
    return chunk
    
def timing03():
    chunk = []
    file_in = open('VSOP2013p1.dat')
    start, end = limits[1][1][0][20]
    chunk = list(islice(file_in, start, end+1))
    file_in.close()
    return chunk
        
def main():
    result1 = timing01()
    result2 = timing02()
    result3 = timing03()
    return (result1, result2, result3)

if __name__ == '__main__':
    
    profile = cProfile.run("main()","prof.prof")
    p = pstats.Stats("prof.prof")
    p.sort_stats("cumulative").print_stats(15)

    result = main()
    print "timing 100 repliche timing01():", timeit.timeit("timing01()", setup = "from __main__ import timing01", number=100)
    print "timing 100 repliche timing02():", timeit.timeit("timing02()", setup = "from __main__ import timing02", number=100)
    print "timing 100 repliche timing03():", timeit.timeit("timing03()", setup = "from __main__ import timing03", number=100)
    print "-" * 80
    print "primo elemento della lista 1:", result[0][0]
    print "primo elemento della lista 2:", result[1][0]
    print "primo elemento della lista 3:", result[2][0]
    
    print "ultimo elemento della lista 1:", result[0][-1]
    print "ultimo elemento della lista 2:", result[1][-1]
    print "ultimo elemento della lista 3:", result[2][-1]
    

Nel mio PC i risultati sono i seguenti (possono differire in modo significativo su altre macchine):

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 ================================
>>> 
Wed May  6 17:23:08 2015    prof.prof

         13 function calls in 0.072 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.004    0.004    0.072    0.072 <string>:1(<module>)
        1    0.010    0.010    0.068    0.068 /home/ubuntu/Scrivania/vsop2013/solution/timing.py:32(main)
        1    0.001    0.001    0.047    0.047 /home/ubuntu/Scrivania/vsop2013/solution/timing.py:7(timing01)
        1    0.046    0.046    0.046    0.046 {method 'readlines' of 'file' objects}
        1    0.006    0.006    0.006    0.006 /home/ubuntu/Scrivania/vsop2013/solution/timing.py:16(timing02)
        1    0.006    0.006    0.006    0.006 /home/ubuntu/Scrivania/vsop2013/solution/timing.py:24(timing03)
        3    0.000    0.000    0.000    0.000 {method 'close' of 'file' objects}
        3    0.000    0.000    0.000    0.000 {open}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}


timing 100 repliche timing01(): 5.62393903732
timing 100 repliche timing02(): 0.71301984787
timing 100 repliche timing03(): 0.654221773148
--------------------------------------------------------------------------------
primo elemento della lista 1:     1   0  0  0  0   0  0  0  0  0    0   0   0   0      0   0  0  0  0.0000000000000000 +00  0.3870983098840000 +00

primo elemento della lista 2:     1   0  0  0  0   0  0  0  0  0    0   0   0   0      0   0  0  0  0.0000000000000000 +00  0.3870983098840000 +00

primo elemento della lista 3:     1   0  0  0  0   0  0  0  0  0    0   0   0   0      0   0  0  0  0.0000000000000000 +00  0.3870983098840000 +00

ultimo elemento della lista 1: 32240   3  0  0  0   0  0  0  0 -7    0   0   0   0      0   0  0  0 -0.7232436513837180 -16 -0.4617790227683706 -15

ultimo elemento della lista 2: 32240   3  0  0  0   0  0  0  0 -7    0   0   0   0      0   0  0  0 -0.7232436513837180 -16 -0.4617790227683706 -15

ultimo elemento della lista 3: 32240   3  0  0  0   0  0  0  0 -7    0   0   0   0      0   0  0  0 -0.7232436513837180 -16 -0.4617790227683706 -15

>>> 

E' abbastanza evidente il pessimo timing della prima funzione: 8 volte più lenta e inefficiente delle successive. Inoltre, la necessità di creare una lista per procedere al suo slicing, oltre a comportare un esagerato impegno di memoria, si rivela impegnativa anche sotto il profilo dei tempi di esecuzione, giustificando la quasi totalità del tempo impiegato dalla funzione. Le due funzioni successive sono decisamente più performanti, a dimostrare che una scelta accurata dei metodi puo' essere determinante, anche quando usiamo funzioni di libreria.

Ho previsto la stampa del primo e dell'ultimo elemento del blocco ottenuto da ogni funzione per mostrare che i risultati prodotti dalle funzioni sono identici. Nel post successivo utilizzeremo il metodo islice per l'estrazione delle righe del file dei parametri e inizieremo, sempre confrontando più metodi, a convertire le righe di testo in variabili numeriche.

Sunday, April 26, 2015

VSOP 2013 - Eliminare i task ripetitivi

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:

  1. 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
  2. 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.
  3. 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.

Friday, April 17, 2015

VSOP2013 - mpmath, gmpy2 e aumento della velocità di esecuzione

Per questo post e i successivi ho pensato di trattare un po' il tema della velocità di esecuzione del codice. Finora non ne ho sentito ll bisogno, in quanto, con routine statiche e limitata dipendenza da file di parametri su disco, anche nel caso delle VSOP87, i tempi di risposta rimanevano adeguati. Con le VSOP2013 i tempi sono significativamente più lunghi, avendo file di parametri con centinaia di migliaia di righe e la necessità di una lettura sequenziale degli stessi con verifica del punto di arresto ad una precisione prefissata. Il precedente post delineava un approccio puramente imitativo del file FORTRAN allegato dal gruppo dell'Observatoire de Paris ai file di parametri presenti sul sito ftp dell'IMCCE. Come scrivevo nello stesso post non pretendo di imitare con il linguaggio Python (interpretato, a tipizzazione dinamica) la stessa velocità esecutiva di FORTRAN (compilato e a tipizzazione statica), ma vorrei tentare comunque di ridurre il tempo di esecuzione, per consentire una fluidità operativa che attualmente non è disponibile.

Per prima cosa vorrei quindi presentarvi i tempi di esecuzione dell'eseguibile ottenuto dalla compilazione in FORTRAN e quindi quelli che raggiungo attualmente con Python. Ho fatto il timing su una macchina equipaggiata con processore AMD A8-Series APU A8-6600K quadcore con 8 giga di RAM :

VSOP2013.f compilato: time ./VSOP2013 real 0m18.422s user 0m18.260s sys 0m0.072s

Il file produce un output su disco (VSOP.out)

python e mpmath

timing del file mp_struct.py, traduzione in linguaggio Python con l'uso della libreria di alta precisione mpmath, eseguito con l'interprete Python2.7.9. La lettura sequenziale dei file prevede una precisione d'arresto calcolata dai parametri S e C mediante la formula

mp_struct.py interpretato con Python 2.7.9 (precisione 1e-12): real 118m56.059s user 118m40.096s sys 0m6.192s

Si tratta chiaramente di una tempistica inaccettabile per calcolare le posizioni di 9 pianeti in 11 date. Accettando una precisione inferiore (1e-7) le cose migliorano un po':

mp_struct.py interpretato con Python 2.7.9 (precisione 1e-7): real 2m2.971s user 2m2.004s sys 0m0.796s

Siamo ancora lontani da un risultato soddisfacente. Ho provato a utilizzare pypy, che è un compilatore jit (just in time) ad alta efficienza. Non permette di includere la libreria mpmath dai moduli aggiunti di Python, ma si puo' fare un'installazione dedicata partendo dai sorgenti, decomprimendoli e lanciando all'interno della directory così ottenuta il comando:

sudo pypy setup.py install

La libreria viene ricreata nell'ambiente di pypy, dove, almeno relativamente al mio codice, funziona perfettamente.

mp_struct.py interpretato con PyPy 2.5.0 (precisione 1e-7): real 0m43.148s user 0m41.128s sys 0m1.024s

Il miglioramento è evidente. L'ottimizzazione creata dal compilatore pypy è decisamente migliorativa portando il tempo esecutivo al 35% del timing iniziale. Rimane ancora un'opzione, prima di passare alla reingegnerizzazione del codice: cambiare libreria di alta precisione. Senza ristrutturare pesantemente il codice, cambiando solo i riferimenti alle classi, ho importato la libreria gmpy2, che è un wrapper python di librerie scritte in C, pertanto molto veloci nell'esecuzione:
gmpy2 is a C-coded Python extension module that supports multiple-precision arithmetic. gmpy2 is the successor to the original gmpy module. The gmpy module only supported the GMP multiple-precision library. gmpy2 adds support for the MPFR (correctly rounded real floating-point arithmetic) and MPC (correctly rounded complex floating-point arithmetic) libraries. gmpy2 also updates the API and naming conventions to be more consistent and support the additional functionality.

Per l'installazione in python si puo' utilizzare pip, previa installazione di alcune dipendenze. Rimando alla pagina dedicata nella documentazione di gmpy2: installazione di gmpy2 su linux

Allego il nuovo codice in forma integrale. Trovate questo, come il precedente scritto con mpmath, nel mio repository git sotto astrometry/vsop2013

import gmpy2 as gmp

import struct

gmp.get_context().precision=200

def cal2jul(year, month, day, hour=0, minute =0, second=0):
    month2 = month
    year2 = year
    if month2 <= 2:
        year2 -= 1
        month2 += 12
    if (year*10000 + month*100 + day) > 15821015:
        a = int(year2/100)
        b = 2 - a + int(a/4)
    else:
        a = 0
        b = 0
    if year < 0:
        c = int((365.25 * year2)-0.75)
    else:
        c = int(365.25 * year2)
    d = int(30.6001 *(month2 + 1))
    return b + c + d + day + hour / 24.0 + minute / 1440.0 + second / 86400.0 + 1720994.5        




class VSOP2013():
    
    def __init__(self, t, planet, precision=1e-7):
        # calculate millennia from J2000
        self.JD = t
        self.t = gmp.div((t - cal2jul(2000,1,1,12)), 365250.0)
        # predefine powers of self.t
        self.power = []; self.power.append(gmp.mpfr(1.0)); self.power.append(self.t)
        for i in xrange(2,21):
            t = self.power[-1]
            self.power.append(gmp.mul(self.t,t))
        # choose planet file in a dict
        self.planet = planet
        self.planets = {'Mercury':'VSOP2013p1.dat',
                        'Venus'  :'VSOP2013p2.dat',
                        'EMB'    :'VSOP2013p3.dat',
                        'Mars'   :'VSOP2013p4.dat',
                        'Jupiter':'VSOP2013p5.dat',
                        'Saturn' :'VSOP2013p6.dat',
                        'Uranus' :'VSOP2013p7.dat',
                        'Neptune':'VSOP2013p8.dat',
                        'Pluto'  :'VSOP2013p9.dat'}
        # VSOP2013 routines precision
        self.precision = precision
        
        # lambda coefficients
        # l(1,13) : linear part of the mean longitudes of the planets (radian).
        # l(14): argument derived from TOP2013 and used for Pluto (radian).
        # l(15,17) : linear part of Delaunay lunar arguments D, F, l (radian).

        self.l = (
            (gmp.mpfr(4.402608631669), gmp.mpfr(26087.90314068555)),
            (gmp.mpfr(3.176134461576), gmp.mpfr(10213.28554743445)),
            (gmp.mpfr(1.753470369433), gmp.mpfr( 6283.075850353215)),
            (gmp.mpfr(6.203500014141), gmp.mpfr( 3340.612434145457)),
            (gmp.mpfr(4.091360003050), gmp.mpfr( 1731.170452721855)),
            (gmp.mpfr(1.713740719173), gmp.mpfr( 1704.450855027201)), 
            (gmp.mpfr(5.598641292287), gmp.mpfr( 1428.948917844273)),
            (gmp.mpfr(2.805136360408), gmp.mpfr( 1364.756513629990)),
            (gmp.mpfr(2.326989734620), gmp.mpfr( 1361.923207632842)),
            (gmp.mpfr(0.599546107035), gmp.mpfr(  529.6909615623250)),
            (gmp.mpfr(0.874018510107), gmp.mpfr(  213.2990861084880)),
            (gmp.mpfr(5.481225395663), gmp.mpfr(   74.78165903077800)),
            (gmp.mpfr(5.311897933164), gmp.mpfr(   38.13297222612500)),
            (gmp.mpfr(0.000000000000), gmp.mpfr(    0.3595362285049309)),
            (gmp.mpfr(5.198466400630), gmp.mpfr(77713.7714481804)),
            (gmp.mpfr(1.627905136020), gmp.mpfr(84334.6615717837)),
            (gmp.mpfr(2.355555638750), gmp.mpfr(83286.9142477147))
            )

        # planetary frequencies in longitude
        self.freqpla = {'Mercury'   :   gmp.mpfr(0.2608790314068555e5),  
                        'Venus'     :   gmp.mpfr(0.1021328554743445e5),
                        'EMB'       :   gmp.mpfr(0.6283075850353215e4),
                        'Mars'      :   gmp.mpfr(0.3340612434145457e4), 
                        'Jupiter'   :   gmp.mpfr(0.5296909615623250e3),  
                        'Saturn'    :   gmp.mpfr(0.2132990861084880e3),
                        'Uranus'    :   gmp.mpfr(0.7478165903077800e2),  
                        'Neptune'   :   gmp.mpfr(0.3813297222612500e2), 
                        'Pluto'     :   gmp.mpfr(0.2533566020437000e2)}  

        
            
        # target variables
        self.ax = gmp.mpfr(0.0)            # major semiaxis
        self.ml = gmp.mpfr(0.0)            # mean longitude
        self.kp = gmp.mpfr(0.0)            # e*cos(perielium longitude)
        self.hp = gmp.mpfr(0.0)            # e*sin(perielium longitude)
        self.qa = gmp.mpfr(0.0)            # sin(inclination/2)*cos(ascending node longitude)
        self.pa = gmp.mpfr(0.0)            # sin(inclination/2)*cos(ascending node longitude)
        
        self.tg_var = {'A':self.ax, 'L':self.ml, 'K':self.kp,
                       'H':self.hp, 'Q':self.qa, 'P':self.pa } 

        # eps  = (23.d0+26.d0/60.d0+21.41136d0/3600.d0)*dgrad
        self.eps = gmp.mpfr((23.0+26.0/60.0+21.411360/3600.0)*gmp.const_pi()/180.0)
        self.phi  = gmp.mpfr(-0.051880 * gmp.const_pi() / 180.0 / 3600.0)
        self.ceps = gmp.cos(self.eps)
        self.seps = gmp.sin(self.eps)
        self.cphi = gmp.cos(self.phi)
        self.sphi = gmp.sin(self.phi)
        
        # rotation of ecliptic -> equatorial rect coords
        self.rot = [[self.cphi, -self.sphi*self.ceps,  self.sphi*self.seps],
                    [self.sphi,  self.cphi*self.ceps, -self.cphi*self.seps],
                    [0.0,        self.seps,            self.ceps          ]]
        
        self.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""")

        self.gmp_ = {
       'Mercury' : gmp.mpfr(4.9125474514508118699e-11),
       'Venus'   : gmp.mpfr(7.2434524861627027000e-10),
       'EMB'     : gmp.mpfr(8.9970116036316091182e-10),
       'Mars'    : gmp.mpfr(9.5495351057792580598e-11),
       'Jupiter' : gmp.mpfr(2.8253458420837780000e-07),
       'Saturn'  : gmp.mpfr(8.4597151856806587398e-08),
       'Uranus'  : gmp.mpfr(1.2920249167819693900e-08),
       'Neptune' : gmp.mpfr(1.5243589007842762800e-08),
       'Pluto'   : gmp.mpfr(2.1886997654259696800e-12)
       }

        self.gmsol = gmp.mpfr(2.9591220836841438269e-04)
        self.rgm = gmp.sqrt(self.gmp_[self.planet]+self.gmsol)
        # run calculus routine
        self.calc()
        

    def __str__(self):
        
        vsop_out = "{:3.13} {:3.13} {:3.13} {:3.13} {:3.13} {:3.13}\n".format(
                                          self.tg_var['A'],
                                          self.tg_var['L'],
                                          self.tg_var['K'],
                                          self.tg_var['H'],
                                          self.tg_var['Q'],
                                          self.tg_var['P'])
        vsop_out += "{:3.13} {:3.13} {:3.13} {:3.13} {:3.13} {:3.13}\n".format(
                                          self.ecl[0],
                                          self.ecl[1],
                                          self.ecl[2],
                                          self.ecl[3],
                                          self.ecl[4],
                                          self.ecl[5])
        vsop_out += "{:3.13} {:3.13} {:3.13} {:3.13} {:3.13} {:3.13}\n".format(
                                          self.equat[0],
                                          self.equat[1],
                                          self.equat[2],
                                          self.equat[3],
                                          self.equat[4],
                                          self.equat[5])
                        

        return vsop_out



    def calc(self):
        with open(self.planets[self.planet]) as file_in:
            terms = []
            b = '*'
            while b != '':
                b = file_in.readline()
                if b != '':
                    if b[:5] == ' VSOP':
                        header = b.split()
                        #print header[3], header[7], header[8], self.t**int(header[3])
                        no_terms = int(header[4])
                        for i in xrange(no_terms):
                            #6x,4i3,1x,5i3,1x,4i4,1x,i6,1x,3i3,2a24
                            terms = file_in.readline()
                            a = self.fmt.unpack(terms)
                            S = gmp.mul(gmp.mpfr(a[18]),gmp.exp10(int(a[19])))
                            C = gmp.mul(gmp.mpfr(a[20]),gmp.exp10(int(a[21])))
                            if gmp.sqrt(S*S+C*C) < self.precision:
                                break
                            aa = 0.0; bb = 0.0;
                            for j in xrange(1,18):
                                aa += gmp.mul(gmp.mpfr(a[j]), self.l[j-1][0])
                                bb += gmp.mul(gmp.mpfr(a[j]), self.l[j-1][1])
                            arg = aa + bb * self.t
                            power = int(header[3])
                            comp = self.power[power] * (S * gmp.sin(arg) + C * gmp.cos(arg))
                            if header[7] == 'L' and power == 1 and int(a[0]) == 1:
                                pass
                            else:
                                self.tg_var[header[7]] += comp
        self.tg_var['L'] = self.tg_var['L'] + self.t * self.freqpla[self.planet]
        self.tg_var['L'] = self.tg_var['L'] % (2 * gmp.const_pi())
        if self.tg_var['L'] < 0:
            self.tg_var['L'] += 2*gmp.const_pi()
        print "Julian date {}".format(self.JD)
        file_in.close()
        ##print self.tg_var
    
        ####    def ELLXYZ(self):

        xa = self.tg_var['A']
        xl = self.tg_var['L']
        xk = self.tg_var['K']
        xh = self.tg_var['H']
        xq = self.tg_var['Q']
        xp = self.tg_var['P']
        
        # Computation
  
        xfi = gmp.sqrt(1.0 -xk * xk - xh * xh)
        xki = gmp.sqrt(1.0 -xq * xq - xp * xp)
        u = 1.0/(1.0 + xfi)
        z = complex(xk, xh)
        ex = abs(z)
        ex2 = ex  * ex
        ex3 = ex2 * ex
        z1 = z.conjugate()
        #
        gl = xl % (2*gmp.const_pi())
        gm = gl - gmp.atan2(xh, xk)
        e  = gl + (ex - 0.1250 * ex3) * gmp.sin(gm)
        e += 0.50 * ex2 * gmp.sin(2.0 * gm)
        e += 0.3750 * ex3 * gmp.sin(3.0 * gm)
        #
        while True:
            z2 = complex(0.0, e)
            zteta = gmp.exp(z2)
            z3 = z1 * zteta
            dl = gl - e + z3.imag
            rsa = 1.0 - z3.real
            e = e + dl / rsa
            if abs(dl) < 1e-15:
                break
        #
        z1  =  u * z * z3.imag
        z2  =  gmp.mpc(z1.imag, -z1.real)
        zto = (-z+zteta+z2)/rsa
        xcw = zto.real
        xsw = zto.imag
        xm  = xp * xcw - xq * xsw
        xr  = xa * rsa
        #
        self.ecl = []; self.equ = {}
        self.ecl.append(xr * (xcw -2.0 *xp * xm))
        self.ecl.append(xr * (xsw +2.0 *xq * xm))
        self.ecl.append(-2.0 * xr * xki * xm)
        #
        xms = xa *(xh + xsw) / xfi
        xmc = xa *(xk + xcw) / xfi
        xn  = self.rgm / xa ** (1.50)
        #
        self.ecl.append( xn *((2.0 * xp * xp - 1.0) * xms + 2.0 * xp * xq * xmc))
        self.ecl.append( xn *((1.0 -2.0 * xq * xq) * xmc -2.0 * xp * xq * xms))
        self.ecl.append( 2.0 * xn * xki * (xp * xms + xq * xmc))
        
        # Equatorial rectangular coordinates and velocity

        #         
        #
        # --- Computation ------------------------------------------------------
        #
        self.equat = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]    
        for i in xrange(3):
           for j in xrange(3):
               self.equat[i]   = self.equat[i]   + self.rot[i][j] * self.ecl[j]
               self.equat[i+3] = self.equat[i+3] + self.rot[i][j] * self.ecl[j+3]



if __name__ == '__main__':
    for planet in ('Mercury', 'Venus', 'EMB', 'Mars', 'Jupiter',
          'Saturn', 'Uranus', 'Neptune', 'Pluto'):
        print "PLANETARY EPHEMERIS VSOP2013  "+ planet + "(TDB)\n"+""" 
  1/ Elliptic   Elements:    a (au), lambda (radian), k, h, q, p        - Dynamical Frame J2000
  2/ Ecliptic   Heliocentric Coordinates:  X,Y,Z (au)  X',Y',Z' (au/d)  - Dynamical Frame J2000
  3/ Equatorial Heliocentric Coordinates:  X,Y,Z (au)  X',Y',Z' (au/d)  - ICRS Frame J2000
"""
        init_date = cal2jul(1890,6,26,12)
        set_date = init_date

        while set_date < init_date + 41000:
            v = VSOP2013(set_date, planet)
            print v
            set_date += 4000
          
        

gmp_struct.py interpretato con Python 2.7.9 (precisione 1e-7): real 0m27.228s user 0m26.184s sys 0m1.000s

Il risultato è incoraggiante: 67% rispetto alla performance di pypy sul file con mpmath. Rinuncio a tentare di compilare gmpy2 sotto pypy, non mi aspetto in questo momento grandi miglioramenti, visto che la gmpy2 è scritta in C. Penso che sia venuto il momento di provare a fare un po' di profiling ed eventualmente riscrivere un po' il codice. Alla prossima.

How to create a virtual linux machine with qemu under Debian or Ubuntu with near native graphics performance

It's been a long time since my latest post. I know, I'm lazy. But every now and then I like to publish something that other people c...