Tuesday, March 25, 2014

Il modulo ctypes e la libreria libswe.so

Siccome la mia class SVG2Tk è ancora un po' indietro di cottura, apro una piccola parentesi sull'accesso diretto alle librerie condivise scritte in linguaggio C.

La libreria pyswisseph, che ho usato nei primi post, è un lavoro splendido e ci permette di effettuare chiamate alla libreria con istruzioni in Python nativo. Se pero' volessimo usare la libreria originaria, installata dal software center di Ubuntu, dovremmo confrontarci con un problema: non è scritta in Python ma in C e compilata verosimilmente con gcc, il compilatore GNU. La stessa pyswisseph è un insieme di funzioni wrapper (involucro) che fanno da ponte fra i due linguaggi. Anche noi, però, possiamo farlo direttamente dal nostro codice.

Tra le librerie integrate nell'installazione di Python c'è un modulo realizzato appositamente per chiamare una libreria in C, a patto che si conoscano le specifiche della libreria, i nomi delle funzioni e i parametri in input e in output.

Fortunatamente la libreria swisseph è molto ben documentata e con un po' di sforzo possiamo effettuare delle chiamate come se utilizzassimo una libreria pitonica. Vediamo come, prendendo a spunto un post del curatore delle pyswisseph Stanislas Marquis link.

Per prima cosa apriamo IDLE e scriviamo l'import di ctypes e di datetime. A questo punto creiamo una variabile shared_lib che costituirà il link alla libreria, attraverso il metodo CDLL di ctypes (o, in alternativa, il metodo cdll.LoadLibrary, sovrapponibile al precedente). L'installazione di swisseph da ubuntu software center mette la libreria libswe.so in una specifica posizione sul disco. Qualora volessimo spostarla nella stessa directory del sorgente Python sarà necessario cambiare il path della libreria stessa. Aggiungiamo inoltre due dizionari per i nomi dei pianeti e dei segni zodiacali.


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
# -*- coding: utf-8 -*-
import datetime
import ctypes as ct 

shared_lib = ct.CDLL('/usr/lib/x86_64-linux-gnu/libswe.so')

pianeti = {0:'Sole',
           1:'Luna',
           2:'Mercurio',
           3:'Venere',
           4:'Marte',
           5:'Giove',
           6:'Saturno',
           7:'Urano',
           8:'Nettuno',
           9:'Plutone',
           10:'Nodo lunare medio',
           11:'Nodo lunare vero'}

segni = {0:'Ariete',
         1:'Toro',
         2:'Gemelli',
         3:'Cancro',
         4:'Leone',
         5:'Vergine',
         6:'Bilancia',
         7:'Scorpione',
         8:'Sagittario',
         9:'Capricorno',
         10:'Acquario',
         11:'Pesci'}

Della libreria swisseph, in questo programma, uso solo due funzioni, una per il calcolo della data giuliana, l'altra per il calcolo della posizione dei pianeti. Vediamo la prima.


1
2
3
4
5
6
7
_julday=shared_lib.swe_julday
_julday.argtypes=[ct.c_int, ct.c_int, ct.c_int, ct.c_double, ct.c_int]
_julday.restype=ct.c_double

def julday(year, month, day, hour, gregflag=1):
    return _julday(ct.c_int(year), ct.c_int(month), ct.c_int(day),
                   ct.c_double(hour), ct.c_int(gregflag))

La prima parte è molto semplice. La libreria ctypes consente di accedere alla funzione nascosta nella shared library swe_julday definendo prima gli argtypes, cioè le tipologie di variabili da usare come parametri di input, poi il restype, cioè il tipo di variabile di ritorno dalla funzione (ricordiamo che il linguaggio C non consente la tipizzazione dinamica, cioè l'assegnazione arbitraria in runtime di un tipo all'atto della creazione della variabile. Essendo un linguaggio compilato, richiede, ovviamente, la tipizzazione di tipo statico).

La chiamata da codice Python risulta molto naturale, in questo caso i parametri passati alla funzione julday vengono trasferiti ai tipi definiti per l'interfaccia al linguaggio C. La funzione restituisce una variabile double, in Python corrisponde a una float e fornisce il valore della data giuliana calcolata da data e ora.

Un po' più complessa risulta, invece, la seconda funzione, che calcola le posizioni planetarie. Di seguito un estratto della guida Programming Interface to the Swiss Ephemeris:

2.1. The call parameters
swe_calc_ut() was introduced with Swisseph version 1.60 and makes planetary calculations a bit simpler. For the steps required, see the chapter  The programming steps to get a planet’s position.
swe_calc_ut() and swe_calc() work exactly the same way except that swe_calc() requires Ephemeris Time( more accurate: Dynamical Time ) as a parameter whereas swe_calc_ut() expects Universal Time. For common astrological calculations, you will only need swe_calc_ut() and will not have to think anymore about the conversion between Universal Time and Ephemeris Time.
swe_calc_ut() and swe_calc() compute positions of planets, asteroids, lunar nodes and apogees. They are defined as follows:
 
int swe_calc_ut ( double tjd_ut, int ipl, int iflag, double* xx, char* serr),
where
tjd_ut     =Julian day, Universal Time
ipl       =body number
iflag    =a 32 bit integer containing bit flags that indicate what kind of computation is wanted
xx       =array of 6 doubles for longitude, latitude, distance, speed in long., speed in lat., and speed in dist.
serr[256] =character string to return error messages in case of error.
 
and
int swe_calc(double tjd_et, int ipl, int iflag, double *xx, char *serr),
same but
tjd_et     =     Julian day, Ephemeris time,  where tjd_et = tjd_ut + swe_deltat(tjd_ut)

Nella dichiarazione della funzione swe_calc compaiono due puntatori (in Python non esistono, o meglio le variabili Python sono dei puntatori, ma non dilunghiamoci), cioè gli indirizzi di memoria nei quali possiamo trovare le variabili xx (che contiene un array di 6 doubles, in Python diremmo una tupla di 6 float) e serr (per l'eventuale messaggio di errore).

I tipi di variabili ctypes possono essere inizializzati nel codice Python con modalità differenti.

Per i tipi più semplici, ctypes.c_int(42) ctypes.c_float(3.1415), ad esempio, il valore tra parentesi viene conservato in una locazione di memoria da cui possono essere recuperati con l'espressione .value. Per esempio:

>>> a = ctypes.c_int(42)
>>> print a
c_int(42)
>>> print a.value
42
>>> b = ctypes.c_float(3.1415)
>>> print b
c_float(3.1414999961853027)
>>> print b.value
3.14149999619

Vediamo ora come si puo' effettuare la traduzione in Python della funzione swe_calc di swisseph utilizzando ctypes.


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
_calc = shared_lib.swe_calc
_calc.argtypes = [ct.c_double, ct.c_int, ct.c_int, ct.POINTER(ct.c_double*6), ct.c_char_p]
_calc.restype = ct.c_int

def calc(julday, planet, flag = 0):

    xx = (ct.c_double * 6)()
    err = ct.c_char_p('')
    if _calc(ct.c_double(julday), ct.c_int(planet), ct.c_int(flag),
             ct.byref(xx), err) == 0:
        return [x for x in xx]
    return err.value

Come potete vedere, il tipo POINTER viene usato per tradurre il puntatore ad una variabile (preceduta da *) di C. Per la variabile serr si usa il tipo predefinito ctypes.c_char_p che è esso stesso un puntatore. La funzione accetta come parametri un double che riprende il risultato del calcolo della data giuliana, una variabile c_int per il numero di posizione del pianeta nella serie e infine il tipo c_int(flag), inizializzato a 0 perchè "If no bits are set, i.e. if iflag == 0, swe_calc() computes what common astrological ephemerides (as available in book shops) supply". Se il valore di ritorno di _calc è zero, cioè non viene resituito un errore, l'array di 6 double va a popolare una list di 6 float contenente vari elementi relativi al pianeta considerato. Di questi ci interessa solo l'elemento con indice 0, riferito alla longitudine in gradi. La lista viene generata con una list comprehension [x for x in xx].

L'ultima parte del programma contiene istruzioni per l'uso delle funzioni C che abbiamo creato in precedenza.


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
def d2hm(x):
    lon30 = x % 30
    segno = int(x) / 30
    grado = int(lon30)
    minuto = int((lon30-grado)*60)
    return (grado, segno, minuto)


adesso = datetime.datetime.now()
in_questo_momento = julday(adesso.year,
       adesso.month,
       adesso.day,
       adesso.hour + adesso.minute/60.0 + adesso.second/3600.0)

for i in range(0,12):
    longitudine = calc(in_questo_momento,i)[0]
    pos = d2hm(longitudine)
    print "%18s %2d° %10s %2d'" % (pianeti[i], pos[0], segni[pos[1]], pos[2])

Lanciando il programma si ottiene il seguente risultato:

              Sole  4°     Ariete 49'
              Luna 24° Capricorno 33'
          Mercurio  9°      Pesci 42'
            Venere 18°   Acquario 18'
             Marte 23°   Bilancia 56'
             Giove 11°     Cancro  1'
           Saturno 22°  Scorpione 53'
             Urano 12°     Ariete  1'
           Nettuno  6°      Pesci  8'
           Plutone 13° Capricorno 28'
 Nodo lunare medio 29°   Bilancia 50'
  Nodo lunare vero 28°   Bilancia 36'

Riferito ad oggi 25 marzo alle 13.37 tempo di Greenwich

Spero che questo post sia stato utile per capire come interfacciarsi ad una shared library compilata in linguaggio C. Ulteriori informazioni sull'uso del modulo ctypes sono rintracciabili alla pagina http://docs.python.org/2/library/ctypes.html#

Se volete creare voi stessi la libreria libswe.so, scaricate il file a questo link. Decomprimetelo e lanciate l'istruzione make, che crea i file oggetto. Quindi create la libreria condivisa, nella stessa directory, con l'istruzione gcc -g -O9 -Wall -fPIC -lm -shared -o libswe.so swedate.o swehouse.o swejpl.o swemmoon.o swemplan.o swepcalc.o sweph.o swepdate.o swephlib.o swecl.o swehel.o. La libreria libswe.so comparirà nella directory. Copiatela dove avete collocato il sorgente .py di cui a questo post. Modificate il riferimento alla libreria nel vostro codice in shared_lib = ct.CDLL('libswe.so') e lanciate in esecuzione il codice python, funzionerà perfettamente.

No comments:

Post a Comment

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...