Implementare CSPICE - Alcune costanti e funzioni comuni

Ora che abbiamo creato una libreria dinamica, è necessario creare l'interfaccia per il linguaggio Python. Il modulo ctypes della libreria standard è il metodo che prediligo, ma ce ne sono altri:

• SWIG, per esempio, utilizzabile con più linguaggi di programmazione
• le Python/C API, che consentono di scrivere codice python dentro un sorgente C (le trovo piuttosto scomode da usare)
• Cython, che è un compilatore statico per Python e per un linguaggio proprio derivato da Pyrex, molto più user friendly delle Python/C API

Io perferisco ctypes perchè diretto e relativamente facile da usare, come vi dimostrero' a breve.

Naturalmente qualcuno potrebbe obiettare che, avendo delle librerie SPICE scritte in C, tanto vale usare il C. Il problema è che C va bene per scrivere routine e libreria, ma è molto complesso da utilizzare in ambiente grafico e molto meno produttivo di Python. Inoltre il debug di C può essere molto più complicato.

Iniziamo dai collegamenti fondamentali.

Linking della libreria dinamica cspice

La volta scorsa abbiamo creato una shared library che abbiamo chiamato libcspice.so, equivalente in ambiente windows di una dll. Creiamo una directory testSPICE, che useremo per lo sviluppo del codice sorgente.

Oltre al file libcspice.so, copiamo nella directory uno o più file de.bsp che recuperiamo dal sito ftp della NASA, sono i cosiddetti DAF/SPK Kernel, che contengono i coefficienti dei polinomi interpolatori usati per calcolare le posizioni di oggetti nello spazio, in particolare i pianeti. Durante la compilazione della libreria abbiamo creato anche una serie di file eseguibili che troviamo nella directory cspice/exe. Uno di questi è brief, che possiamo usare per trovare le specifiche di ognuno di questi file. Per esempio, se abbiamo scaricato de421.bsp e vogliamo sapere cosa contiene, copiamo brief nella directory test e lo lanciamo in esecuzione:

./brief de421.bsp

L'output del programma sarà il seguente:

BRIEF -- Version 4.0.0, September 8, 2010 -- Toolkit Version N0066
 
 
Summary for: de421.bsp
 
Bodies: MERCURY BARYCENTER (1)  SATURN BARYCENTER (6)   MERCURY (199)
        VENUS BARYCENTER (2)    URANUS BARYCENTER (7)   VENUS (299)
        EARTH BARYCENTER (3)    NEPTUNE BARYCENTER (8)  MOON (301)
        MARS BARYCENTER (4)     PLUTO BARYCENTER (9)    EARTH (399)
        JUPITER BARYCENTER (5)  SUN (10)                MARS (499)
        Start of Interval (ET)              End of Interval (ET)
        -----------------------------       -----------------------------
        1899 JUL 29 00:00:00.000            2053 OCT 09 00:00:00.000

I corpi celesti di cui è possibile avere un'effemeride sono i pianeti o i loro baricentri, il sole, la terra e la luna. Tra parentesi sono indicati i codici NAIF corrispondenti. L'ultima riga contiene gli estremi temporali, in questo caso si va dal 29 luglio 1899 00:00:00 Ephemeris Time al 9 ottobre 2053 00:00:00. Se vogliamo fare ricerche storiche più lontane nel tempo dobbiamo caricare un file effemeride più grande e quindi necessariamente più voluminoso. Il sito della NASA contiene molte informazioni, tutorial e letture tecniche per capire come si usano, come si costruiscono, come si importano ed esportano i file SPK e come è possibile modificarli, per esempio per creane un estratto e salvarlo in un nuovo file. Le utilità della sottodirectory exe sono utili anche a questo scopo

La prima cosa che dobbiamo fare è creare un modulo che useremo nel seguito per creare i nostri sorgenti python.

Possiamo aprire il nostro editor Python (idle, Wingware, Geany, PyCharm sono i più diffusi in ambiente Ubuntu/Mint), e cominciare a buttare giù qualcosa.

#!/usr/bin/env python
import ctypes as ct
from tools import *

Il file tools.py è quello che ho creato due post addietro, contiene un po' di funzioni per lavorare con gradi e radianti, funzioni di calendario ecc. Se tornate indietro nel blog trovate tutto fatto. Ricordatevi di mette nella cartella un file vuoto chiamato __init__.py, che serve a considerare i file contenenti le funzioni di utilità come moduli.

Adesso generiamo il link alla libreria libcspice.so, che abbiamo opportunamente messo nella directory.

ct.shared_lib = ct.CDLL('./libcspice.so')

Da questo momento in poi ct.shared_lib è l'anello di congiunzione tra Python e C. A questa funzione collegheremo tutte le funzioni della libreria cspice, man mano che ci serviranno.

Inziamo da alcune costanti:

SpiceDouble pi_c ( void )
SpiceDouble dpr_c ( void )
SpiceDouble spd_c ( void )

la prima è semplicemente il valore di pi greco, la seconda è il coefficiente di conversione radianti->gradi, la terza è il numero di secondi trascorsi dall'epoca J2000 (1° gennaio 2000, tempo delle effemeridi) in un giorno. Per usare queste costanti dobbiamo creare l'equivalente in python.

_spd_c = ct.shared_lib.spd_c
_pi_c = ct.shared_lib.pi_c
_dpr_c = ct.shared_lib.dpr_c

Ho usato l'underscore iniziale perchè riservo il nome in piano alla funzione python che scrivero' successivamente.

Creiamo le funzioni proprie in Python:

def spd_c():
    _spd_c.restype = ct.c_double
    return _spd_c()

def pi_c():
    _pi_c.restype = ct.c_double
    return _pi_c()

def dpr_c():
    _dpr_c.restype = ct.c_double
    return _dpr_c()

Fin qui è tutto molto semplice: dalla definizione C delle costanti desumo che non ricevono parametri di input(per forza, sono costanti) e restituiscono invece un valore SpiceDouble, che corrisponde ad un double di python. Per cui ognuna di queste costanti è chiamata dal codice Python come funzionePython, che gestisce all'interno la funzioneC con underscore e la restituisce con un return _funzioneC.

Mi fermo qui con questo post, nel prossimo continuiamo l'implentazione e facciamo qualche prova di collaudo. A presto.

Don't reinvent the wheel - non reinventiamo la ruota, CSPICE toolkit è perfetto anche per un astrologo

Nei post passati ho posto spesso l'enfasi sull'autoproduzione di software. Se avete cominciato a programmare, come me, alla fine degli anni '80, ricorderete che i compilatori dell'epoca non erano sempre facilmente disponibili e, almeno in ambiente Windows, nemmeno gratuiti. Il primo software di astrologia che ho usato per il DOS di allora era copiato da un numero della rivista del CIDA, scritto in basic, che avevo dotato di una misera interfaccia grafica, completato da alcune routine scritte da me e compilato con Quick Basic. Ne ho ancora la copia, che gira sotto DosBox, e questa a fianco è la sua immagine in esecuzione.

Sono andato sempre un po' di corsa, per cui il programma non ha un titolo decente, l'avevo chiamato Prova70 in base al numero di release, e tale è rimasto.

Da metà anni 80 in avanti ho iniziato a usare visual basic versione 5 poi 6, che è rimasto un po' l'unica piattaforma fino all'avvento di .NET e Visual Basic 2005, 2008 e successivi. Oggi uso soprattutto linux Debian, Ubuntu e derivate, che mettono a disposizione una quantità di software di elevata qualità che fa buona concorrenza a sistemi operativi più diffusi presso gli utenti desktop (mi riferisco a Windows e Mac, ovviamente), linguaggio di programmazione python per quasi tutto (ma pasticcio volentieri con Tcl/Tk, l'ho scoperto da poco meno di un anno ma me ne sono innamorato subito), Android Studio per i sistemi mobile e Debian ARMhf per la raspberry pi come ambiente di sperimentazione (ve ne parlo tra poco, a proposito).

Non c'è maggiore soddisfazione che farsi le cose da sè ma, come scrivo nel titolo, non è che dobbiamo sempre reinventare la ruota, proprio oggi che internet ci dà accesso alle migliori librerie software disponibili. Per questo motivo oggi vi mostro come sono riuscito a installare il toolkit CSPICE, che è lo standard per il calcolo astronomico e non solo, visto che ci mette in grado, volendo, di calcolare il movimento della nostra astronave e farcela parcheggiare, con discreta approssimazione, nel giardino di casa.

CSPICE è la versione in linguaggio C del toolkit elaborato dalla NASA/JPL come "an ancillary information system that provides scientists and engineers the capability to include space geometry and event data into mission design, science observation planning, and science data analysis software."

E' molto esteso e complesso, il solo core contiene più di 1500 funzioni, che spaziano dalla gestione e conversione di data e ora, al calcolo di posizione e movimento di oggetti quali pianeti, asteroidi, satelliti e astronavi, anche nelle relazioni geometriche tra loro e altri oggetti vicini e lontani, alla conversione tra sistemi di coordinate celesti e altro ancora. Per un'esposizione dettagliata delle funzioni disponibili vi invito a leggere la pagina della documentazione di CSPICE.

A questo punto pero' devo fare un

DISCLAIMER - AVVERTENZA

Non so quanto l'astrologia sia considerata conforme alla dignità scientifica delle istituzioni che hanno creato il software. Alcune importanti effemeridi basate su kernel JPL sono utilizzate da software open source e commerciali (vedi Swiss Ephemeris, di cui vi ho parlato in precedenza), per cui mi sento relativamente sicuro di non abusare dei termini della licenza. In ogni caso non fornisco software basato su librerie NASA/JPL ma mi limito a dare qualche suggerimento per un libero e disinteressato utilizzo (volendo anche in campo astronomico, l'astrologo si accontenta di poco, qualcuno dei metodi che descrivo puo' essere benefico anche per astronomi dilettanti e professionisti). Se le mie procedure risultassero in qualche modo difettose, comunque, non me ne faccio carico oltre i limiti ragionevoli della sperimentazione. Ognuno di voi è libero di prendere quello che faccio, copiarlo, modificarlo, venderlo (ammesso che abbia un valore commerciale), sono anche disponibile a dare qualche suggerimento in privato, ma non assumo la responsabilità che deriva dall'uso di queste procedure, ripeto, sono anch'io uno sperimentatore e un dilettante.

Passiamo alle cose divertenti.

Premessa

Non so perchè la NASA/JPL mette a disposizione toolkit in C e FORTRAN per ambienti x86 e x64, per PC, MAC Solaris e SPARC e si dimentica dei processori ARM, se volete compilare il toolkit in una raspberry PI 2 dovete fare qualche piccolo ritocco alla procedura di compilazione. Per prendere il caso più generale (se siete fanatici come me e volete installare il toolkit in uno smartphone cinese da 100 euro, come ho fatto io, siete i benvenuti) dovete in primo luogo trovare il software e compilarlo nel vostro ambiente operativo. Facciamo qualche caso, guardando la tabella delle versioni disponibili per il linguaggio C:

Download e decompressione

• Mac/Intel, OSX, Apple C, 32bit
• Mac/Intel, OSX, Apple C, 64bit
• PC, CYGWIN, gCC, 32bit
• PC, CYGWIN, gCC, 64bit
• PC, Linux, gCC, 32bit
• PC, Linux, gCC, 64bit
• PC, Windows, MS Visual C, 32bit
• PC, Windows, MS Visual C, 64bit
• Sun/Intel, Solaris, Sun C, 32bit
• Sun/Intel, Solaris, Sun C, 64bit
• Sun/SPARC, Solaris, gCC, 32bit
• Sun/SPARC, Solaris, gCC, 64bit
• Sun/SPARC, Solaris, Sun C, 32bit
• Sun/SPARC, Solaris, Sun C, 64bit

Io uso una Linux Mint 18.2 64bit, basata su Ubuntu 16.04, per cui la versione PC; Linux, gCC, 64bit mi va benissimo. Cliccando il link si apre una nuova pagina, dove vediamo, in prima posizione nell'elenco, la voce cspice.tar.Z. Facciamo il download dove preferiamo, magari sul desktop.

Apriamo una finestra terminale nella directory di salvataggio e lanciamo il comando:

tar xvf cspice.tar.Z

Diamo invio e aspettiamo che venga costruita una nuova directory contenente tutti i file del toolkit.

Compilazione statica

Nella directory cspice appena creata troviamo un file di installazione makeall.csh. Vi prego di notare l'estensione csh, si tratta di uno script per la shell csh, che è un po' diversa dalla bash che usiamo abitualmente in ambiente linux. Procediamo al suo download, la useremo forse solo per questa procedura, ma ne vale la pena:

sudo apt install csh

Ora possiamo lanciare la compilazione, con il comando

sudo csh ./makeall.csh

Ci vorrà un po', abbiate pazienza.

Al termine del processo, nella sottodirectory lib troviamo il file cspice.a, che è la libreria statica per il linker di C. Non è quella che ci serve per Python, dobbiamo compilare diversamente. Il lancio dello script ha però consentito di decomprimere ulteriormente i file, creandone di nuovi. Di questi abbiamo bisogno per il prossimo passaggio procedurale.

Variante per processori ARM

In una raspberry pi è installato, di solito, il sistema operativo linux Debian jessie sotto distribuzione Raspbian. Non è affatto diverso dalla Debian per x86 nel funzionamento, ma lo è in relazione all'architettura del sistema. Se lanciamo lo script di cui al paragrafo precedente su raspberry pi va subito in crash. Lo script dovrebbe consentire di rimuovere i riferimenti a tale tipo di architettura, ma dopo alcuni tentativi ho scelto una via più semplice e funzionale. Entrate da teminale nella directory cspice e lanciate il seguente comando:

grep -rl "m64" . --include "*.csh" | sudo xargs sed -i "s/\-m64//g"

oppure, se avete scaricato la versione 32 bit del toolkit

grep -rl "m32" . --include "*.csh" | sudo xargs sed -i "s/\-m32//g"

Il comando rimuoverà in modo ricorsivo tutti i riferimenti del compilatore alle architetture 64 o 32bit, per cui abiliterà il funzionamento dello script makeall.csh che abbiamo già visto in precedenza. In questo modo avremo effettuato una compilazione statica del toolkit anche in ambiente ARM, con grande gioia di tutti gli astrofili possessori di raspberry pi (o di smartphone cinesi).

Compilazione dinamica

Per usare python nel calcolo astronomico ci sono diverse interessanti possibilità:

• La libreria AstroPy è una community based library, complessa e articolata, decisamente orientata all'osservazione astronomica e al relativo imaging. Per un astrologo è eccessiva e, a mio parere, piuttosto difficile da capire e usare

• La libreria Skyfield di Brandon Rhodes, elegante e pitonica, utilizza i kernel (ne parliamo presto) della JPL e quindi assicura un'elevata precisione di calcolo. Fa seguito alla libreria pyephem dello stesso autore, prevalentemente basata sul wrapping di funzioni c, mentre l'attuale è basata su pure python e utilizza la libreria numerica numpy per l'elevata efficienza nel calcolo numerico.

• Utile segnalare anche SpiceyPy di Andrew Annex, che è un python wrapper per lo SPICE toolkit.

Ho parlato, in un precedente post, del modulo ctypes di Python (incluso nella standard library) e ho fatto vedere come si puo' compilare dinamicamente le Swiss Ephemeris e utilizzare ctypes per scrivere delle funzioni wrapper. Il metodo che illustrero', applicato al JPL toolkit, è proprio questo. Pur compilando l'intera libreria, avremo una shared library di basso peso (meno di cinque Mb nella mia raspberry, con tutte le 1500 e oltre funzioni del toolkit, ma potremo artigianalmente scrivere solo i wrapper che ci interessano per il calcolo di posizione, lasciando agli astronomi funzioni più elaborate e complesse.

Per fare la compilazione dinamica della libreria CSPICE entriamo, tramite terminale, nella directory src, quindi nella sottodirectory cspice, dove troviamo molti file sorgenti in c (estensione .c) e file header (estensione .h). Per creare una libreria dinamica non dobbiamo fare altro che usare il seguente comando:

gcc -Wall -fPIC -shared -o libcspice.so *.c -lm -ldl

Vedrete un sacco di warnings ma nulla che blocchi la compilazione, probabilmente ci sono piccole anomalie non corrette nei sorgenti. Anche stavolta, portate pazienza, ci vorrà un po' di tempo. In questo caso sia che lavoriate su PC sia che usiate una raspberry, il comando è identico.

Una volta creata la libreria libcspice.so, copiatela in una nuova directory a cui darete il nome che preferite, relativo al vostro progetto python.

Per iniziare, ci servono altri due file, che metterete nella stessa directory:

• uno dei file kernel planetari di JPL, suggerisco il de421.bsp, che è abbastanza leggero (16 mega) e che trovate a questo indirizzo.

• un kernel che contiene informazioni relative al deltat, lo scarto temporale fra tempo degli orologi atomici e il tempo terrestre, che è fissato periodicamente da una autorità (IERS), ve ne ho anticipato in un post precedente, attualmente vale 37 secondi (leapseconds). Lo trovate a quest'altro indirizzo.

Con questi strumenti nella directory abbiamo quanto ci serve per cominciare a sviluppare codice python usando il JPL toolkit. Il prossimo post (credo più d'uno) sarà proprio dedicato a questo. Al prossimo post.

Un po' di tools di base

La libreria standard di Python è fornitissima di funzioni pronte all'uso, ma per il nostro modesto lavoro di programmatori di software astrologico è bene che ci dotiamo di strumenti fatti apposta per lo scopo.

Per prima cosa apriremo un file tools.py con un editor a nostra scelta, in cui raccoglieremo tutte le funzioni che andremo via via creando in modo da poterle richiamare facilmente.

Per prime definiremo delle funzioni goniometriche di supporto alle operazioni sulle coordinate.

La libreria standard ha già queste funzioni, pero' si applicano ai gradi misurati in radianti (in cui l'angolo piatto vale 2 pigreco, per intenderci), mentre molte delle formule disponibili usano la misura in gradi.

Per prime le funzioni seno, coseno, tangente e arctangente2 (è la funzione inversa della tangente, ma a differenza della funzione arctangente normale effettua alcuni controlli sul segno degli operandi, risolvendo il dubbio sul corretto posizionamento dell'angolo nei quadranti, come si vede nellla tabella seguente):

Python 2.7.13 (default, Jan 19 2017, 14:48:08) 
[GCC 6.3.0 20170118] on linux2
Type "copyright", "credits" or "license()" for more information.
>>> 
======== RESTART: /home/ubuntu/Scrivania/proveStroBlog/provaatan2.py ========
arctan(x) for x varying from 0° to 360° step 30°

x              atan(sin(x)/cos(x))     atan2(sin(x),cos(x))
 30.000               30.000               30.000
 60.000               60.000               60.000
 90.000               90.000               90.000
120.000              -60.000              120.000
150.000              -30.000              150.000
180.000               -0.000              180.000
210.000               30.000             -150.000
240.000               60.000             -120.000
270.000               90.000              -90.000
300.000              -60.000              -60.000
330.000              -30.000              -30.000
360.000               -0.000               -0.000
>>> 

Ecco le prime funzioni goniometriche in gradi e una funzione per riportare un angolo all'angolo giro fondamentale che toglie il segno meno se presente:

import math

def sin(x):
        return math.sin(math.radians(x))

def cos(x):
 return math.cos(math.radians(x))

def tan(x):
        return sin(x)/cos(x)

def atan(x):
        return math.degrees(math.atan(x))

def atan2(y , x):
 return math.degrees(math.atan2(y, x))

def reduce360(x):
 return x % 360.0

Ora ci servono due funzioni per convertire un angolo dalla notazione sessagesimale alla decimale e viceversa. Decido che il segno che compare davanti alle ore darà il segno a tutta la variable decimale, mentre, al contrario, nella conversione da decimale a sessagesimale solo le ore riceveranno un eventuale segno negativo. La funzione, oltre che alle ore, puo' applicarsi immodificata anche agli angoli, sempre in notazione sessagesimale.

def dms2ddd(hour, minute, second):
    """ from sexagesimal to decimal
        the sign of hour variable is automatically applied to minutes and seconds
    """
    if hour < 0:
        sign = -1
        hour *= sign
    else:
        sign = 1
    return (hour+minute/60.0+second/3600.0)*sign

def ddd2dms(dec_hour):
    """ from decimal to sexagesimal representation of hours and angles.
        the sign of dec_hour variable is applied only to hours variable
        see the dms2ddd function for comparison
    """
    if dec_hour < 0:
        sign = -1
        dec_hour *= sign
    else:
        sign = 1
    total_seconds = int(dec_hour * 3600.0+.5)
    seconds = total_seconds % 60 
    total_minutes = int((total_seconds - seconds)/60.0)
    minutes = total_minutes % 60 
    hours = int((total_minutes - minutes)/60.0)
    return (hours * sign, minutes, seconds)

Ci serve inoltre una funzione che converta una data dal formato giorno, mese, anno, ora minuto e secondo, scomodissimo per il calcolo, in quella comoda notazione che corrisponde alla data juliana.

def cal2jul(year, month, day, hour=0, minute=0, second=0):
    """ converts calendar date to julian date
        this routine and the following are built following Duffet Smith /Zwart instructions
        as given in Peter Duffett-Smith-Zwart Practical Astronomy with your Calculator or Spreadsheet
        Fourth Edition, Cambridge University Press, Fourth Ed. 2011
        For an easier use of the function, hours minutes and seconds are defaulted to 0, so it's
        not necessary to give them as parameters when the hour is 00:00:00 
    """
    month2 = month
    year2 = year
    if month2 <= 2:
        year2 -= 1
        month2 += 12
    else:
        pass
    if (year*10000 + month*100 + day) >= 15821015:
        a = math.trunc(year2/100.0)
        b = 2 - a + math.trunc(a/4.0)
    else:
        a = 0
        b = 0
    if year < 0:
        c = math.trunc((365.25 * year2)-0.75)
    else:
        c = math.trunc(365.25 * year2)
    d = math.trunc(30.6001 *(month2 + 1))
    return b + c + d + day + hour / 24.0 + minute / 1440.0 + second / 86400.0 + 1720994.5

Infine ecco due funzioni che servono per il calcolo del tempo siderale di Greenwich e del tempo siderale locale, necessari ad allineare il tempo locale con un riferimento alle cosiddette "stelle fisse" (che tanto fisse non sono, ma per noi astrologi il riferimento statico è l'eclittica e la sua divisione in segni zodiacali, che, per inciso, non c'entrano niente con le costellazioni).

def ut2gst(year, month, day, hour, minute, second):
    """ Sidereal time is a time-keeping system astronomers use to keep track of the direction to point
        their telescopes to view a given star in the night sky.
        Briefly, sidereal time is a "time scale that is based on the Earth's rate of rotation measured
        relative to the fixed stars." (source Wikipedia)
        This routine converts Universal Time to Sidereal Time for Greenwich (Greenwich Sidereal Time)
    """
    jd = cal2jul(year, month, day)
    S = jd - 2451545.0
    T = S/36525.0
    T0 = (6.697374558 + 2400.051336 * T+ 0.000025862 *T*T) % 24
    UT = dms2ddd(hour, minute, second)*1.002737909
    GST = (UT + T0) % 24
    return GST

e

def gst2lst( gst, long_degree, long_minute, long_second=0):
    """ Corrects GST for a different location on the Earth
    """
    lg = dms2ddd(long_degree, long_minute, long_second)/15.
    lst = (gst + lg) % 24.
    return lst

Nel prossimo post vedremo qual è il significato di queste due ultime funzioni.

Trovare le coordinate e la timezone di un luogo

Prima di entrare nel merito delle tecniche di calcolo della domificazione, dobbiamo risolvere un problema: come trovo, in modo programmatico, le coordinate di un luogo geografico? Posso usare un datebase di località, naturalmente, vi ho già mostrato come ricavare una semplice tabella con milioni di località e trasformarla in un database sqlite, qualche post addietro, abbiamo anche costruito un widget apposito. Il problema è che i database sono, in genere, molto pesanti da portarsi dietro, si arriva facilmente a qualche decina di mega e non è detto che siamo sempre in condizioni di installarli nel memoria di un dispositivo, se abbiamo qualche limitazione di risorse.

Considerato che il web è un enorme deposito di informazioni, comprese quelle di natura geografica, possiamo pensare, volendo, a utilizzare un servizio REST come quello che mette a disposizione Google: Google Maps Geocoding API, che è gratuito per usi limitati e richiede l'apertura di un account, una API key, e il pagamento oltre una certa soglia di utilizzo, oppure quello di geonames.org che è il sito che mette a disposizione quell'enorme database che ho usato per costruire il mio widget e sicuramente molti altri.

Oggi voglio invece provare a usare un metodo più ruvido e gratuito, basato su wikipedia e su alcune librerie python. Ve ne parlo succintamente, un trafiletto per ogni libreria.

Wikipedia

Per ogni luogo nel mondo, per quanto piccolo, Wikipedia fornisce una pagina in cui sono riportate le coordinate geografiche del luogo. Spulciando nel sorgente della pagina, non è difficile individuarle ed estrarle usando le tecniche delle espressioni regolari, presenti nella Python Standard Library. Preferisco, fra tutte, l'edizione inglese, cosa di cui dobbiamo tener conto quando cerchiamo un nome di città. In quel grande contenitore di librerie python che è Pypi ho trovato un'API per wikipedia pronta per l'uso, che dovremo quindi installare preliminarmente :

# Se non l'avete già fatto:
sudo pip install python-pip <Invio>
# quindi usate pip per installare wikipedia
pip install wikipedia <Invio>

Per la necessità o meno di usare i diritti di amministratore o per Windows o Mac vi prego di controllare la documentazione relativa all'uso di pip. La libreria python wikipedia è soggetta a licenza MIT.

Requests: HTTP for Humans

Se avete mai usato urllib e urlib2 apprezzerete questo piccolo gioiello di Kenneth Reitz, che consente di accedere al contenuto di una pagina web, anche in forma di xml o json dove previsto, in un soffio.

Per installare questa libreria si procede come prima:

pip install requests

La licenza d'uso è Apache 2.0

timezonefinder

Se Wikipedia ci aiuta moltissimo a trovare le coordinate geografiche, usa una rappresentazione della timezone del luogo basata sullo scarto orario da UTC, che diventa complicata da gestire. Con questa libreria posso facilmente ottenere, in base alle coordinate geografiche, la timezone del luogo espressa come "Europe/Rome" e simili, lasciando ad altre librerie il compito di gestire l'ora estiva

Installazione:

pip install timezonefinder

La libreria è di J. Michelfeit, ed è pubblicata su Pypi con licenza MIT

MyLittleTownFinder

Veniamo dunque al codice con cui gestire tutte le complesse operazioni di ricerca ed estrazione delle coordinate geografiche.

Wikipedia è abbastanza regolare nel catalogare i nomi di città, per cui di fronte a toponimi unici, di solito basta comporre l'url così: https://en.wikipedia.org/wiki/nomedelposto e siamo nella pagina giusta. Non sempre è così semplice, per esempio se un nome di città è riportato più volte, per esempio Paris è la capitale della Francia, ma anche una ridente località del Texas (che ha dato il titolo ad un bellissimo film di Wim Wenders).

In questo e simili casi, dobbiamo escogitare un trucco per individuare la pagina di interesse tra le molte possibili. La prima funzione che vi presento, debitamente commentata, fa proprio questo.

import sys
import re
import requests
import wikipedia as wk
from timezonefinder import TimezoneFinder as tfz


def find_a_place_wikipage():

    # first input the name of a place, in English, if possible 

    town = raw_input("city, town, or village: ")

    # the wikipedia search function gets a list of all the available pages

    urls = wk.search(town)

    # the list is printed to the screen, one element at a time,
    # each element preceded by its ordinal numbers. The user can choose an
    # element or type 99 to exit 

    for x,y in enumerate(urls):
        print x,y
    print "99 to quit"
    n = 9999
    while n not in range(len(urls)):
        n = int (raw_input("select a link : "))
        if (n==99):
            sys.exit()
    # the choosen page is returned to the caller
    return urls[n]

Questa funzione cerca di inviduare la pagina corrispondente alla città che cerchiamo. Nel caso non trovi una voce che ci sembra quella corretta, perchè per esempio cerchiamo Venice ma non quella italiana, magari quella in Florida o in Illinois, possiamo aggiungere queste specificazioni geografiche alla stringa di ricerca, il motore di ricerca punterà con maggiore facilità alla voce giusta.

La seconda funzione è quella che gestisce il recupero dei dati dalla pagina che abbiamo individuato. E' abbastanza semplice e commentata, per cui ve la lascio esaminare senza ulteriori discussioni.

def get_coordinates(location):

    # location is given by the caller
    base_url = 'https://en.wikipedia.org/wiki/'

    # a complete url is formed
    url = base_url + location

    # a page request is sent, if HTTP request is successful
    # then the procedure goes on, else is stopped
    page = requests.get(url)
    if page.status_code == requests.codes.ok:
        # the html text is retrieved
        b = page.text
        pattern = '(wgCoordinates":{"lat":)([+-]?\d*\.?\d*)(,"lon":)([+-]?\d*\.?\d*)'
        m = re.search(pattern, b) 
        # a regular expression is compared to the text
        # if successful, the coordinates are extracted and converted to float
        if m:
            latitudine = float(m.group(2))
            longitudine = float(m.group(4))
            # timezonefinder is used to precisely locate the related timezone
            tf = tfz()
            timezone = tf.timezone_at(lat=latitudine, lng=longitudine)
            # all three results are returned back to the caller as a tuple
            return (latitudine, longitudine, timezone)
        else:
            # unsuccessful operation  
            print"not found, sorry"    
            return None
    else:
        # error in HTTP request 
        print "requests error"
        return None

Per finire qualche riga di codice per eseguire in serie le due funzioni.

if __name__ == '__main__':
    location = find_a_place_wikipage()
    result = get_coordinates(location)
    if result:
        lat, lon, timezone = result
        print "latitude {0} longitude {1} timezone {2} ".format(lat, lon, timezone) 

Bene, facciamo qualche prova.

Python 2.7.13 (default, Jan 19 2017, 14:48:08) 
[GCC 6.3.0 20170118] on linux2
Type "copyright", "credits" or "license()" for more information.
>>> 
========= RESTART: /home/ubuntu/Scrivania/proveStroBlog/cityurls.py =========
città: Roma
0 Roma
1 Representation oligonucleotide microarray analysis
2 A.S. Roma
3 Roma (mythology)
4 Ruska Roma
5 Roma, Queensland
6 Servitka Roma
7 Atletico Roma F.C.
8 Colonia Roma
9 Town of Roma
99 to quit
select a link : 0
not found, sorry
>>> 
========= RESTART: /home/ubuntu/Scrivania/proveStroBlog/cityurls.py =========
città: Rome
0 Rome
1 Rome, Georgia
2 Founding of Rome
3 History of Rome
4 Rome (TV series)
5 Rome Prize
6 Third Rome
7 Diocese of Rome
8 Sack of Rome
9 Province of Rome
99 to quit
select a link : 0
latitude 41.9 longitude 12.5 timezone Europe/Rome 
>>> 

Come potete vedere, ho inserito il nome di città nella dicitura italiana, ma sulla wikipedia in lingua inglese non trova niente. Una volta inserito il nome inglese, tutto fila liscio.

Un'altra prova:

Python 2.7.13 (default, Jan 19 2017, 14:48:08) 
[GCC 6.3.0 20170118] on linux2
Type "copyright", "credits" or "license()" for more information.
>>> 
========= RESTART: /home/ubuntu/Scrivania/proveStroBlog/cityurls.py =========
città: moscow
0 Moscow
1 MoSCoW method
2 Moscow, Idaho
3 Moscow Metro
4 FC Moscow
5 Administrative divisions of Moscow
6 Moscow Conservatory
7 Mayor of Moscow
8 Grand Duchy of Moscow
9 Moscow Time
99 to quit
select a link : 0
latitude 55.75 longitude 37.6166666667 timezone Europe/Moscow 
>>> 

Notate che c'è una Moscow anche nell'Idaho, USA. Se selezionassi quel link avrei questo risultato:

Python 2.7.13 (default, Jan 19 2017, 14:48:08) 
[GCC 6.3.0 20170118] on linux2
Type "copyright", "credits" or "license()" for more information.
>>> 
========= RESTART: /home/ubuntu/Scrivania/proveStroBlog/cityurls.py =========
città: moscow
0 Moscow
1 MoSCoW method
2 Moscow, Idaho
3 Moscow Metro
4 FC Moscow
5 Administrative divisions of Moscow
6 Moscow Conservatory
7 Mayor of Moscow
8 Moscow Time
9 Grand Duchy of Moscow
99 to quit
select a link : 2
latitude 46.73 longitude -117.0 timezone America/Los_Angeles 
>>> 

Facciamo direttamente la ricerca di Paris, Texas per vedere come se la cava.

Python 2.7.13 (default, Jan 19 2017, 14:48:08) 
[GCC 6.3.0 20170118] on linux2
Type "copyright", "credits" or "license()" for more information.
>>> 
========= RESTART: /home/ubuntu/Scrivania/proveStroBlog/cityurls.py =========
città: paris texas
0 Paris, Texas
1 Paris, Texas (film)
2 Paris, Texas (band)
3 Lamar County, Texas
4 Paris Independent School District
5 Paris, Texas (disambiguation)
6 The Paris News
7 Paris (disambiguation)
8 Eiffel Tower (Paris, Texas)
9 First Church of Christ, Scientist (Paris, Texas)
99 to quit
select a link : 0
latitude 33.662508 longitude -95.547692 timezone America/Chicago 
>>> 

Mi sembra che funzioni benino. A questo punto penso che abbiamo un po' di materiale per andare oltre e iniziare a parlare della domificazione. A fra un po'.

Il qui e ora in astrologia

Ogni volta che osserviamo l'ora nel nostro orologio, per qualunque necessità, astrologia compresa, dobbiamo essere consapevoli che l'ora non viene generata nè misurata dal nostro orologio. Tutto quello che l'orologio sa fare, con metodo meccanico o usando un oscillatore al quarzo, è generare eventi periodici.

Quando parliamo di ora civile facciamo riferimento ad accordi e convenzioni che discendono da fenomeni geografici e astronomici, ma anche da scelte amministrative e politiche.

Se osserviamo un orologio in Italia, l'ora che leggiamo è unica in tutto l'ambito territoriale nazionale italiano, ed è adottata anche da un certo numero di altri Paesi che, da Nord a Sud, dalla Norvegia in Nord Europa alla Namibia in Africa, condividono, all'incirca, la nostra posizione a est rispetto ad un meridiano di riferimento, come diciamo comunemente, hanno il nostro stesso fuso orario. I Peasi con estensione territoriale maggiore, come gli Stati Uniti, la Russia e la Cina, non possono che adottare più fusi orari e quindi ore diverse nei diversi territori che li compongono. Per poter viaggiare da un Paese all'altro, nello stesso fuso orario, dobbiamo inoltre tener conto delle regole associate all'ora estiva, che possono cambiare da Paese a Paese.

L'appartenenza ad un fuso orario piuttosto che all'altro genera strani effetti. Spagna e Norvegia adottano il fuso orario dell'Europa Centrale, con un disallineamento rispetto alla posizione del Sole, in ritardo e rispettivamente in anticipo, che ha ripercussioni sulla percezione della distribuzione della luce del giorno in questi paesi. Con un'escursione da 8° Ovest a 23° Est, anche se l'orologio indica la stessa ora, ci sono due ore di differenza nella percezione, per esempio, del sorgere e tramonto del Sole, anche alla stessa latitudine.

Per usare l'ora in Astrologia è pertanto necessario riportarla alla sua fondamentale natura di tempo locale, allineato con la reale posizione del Sole per un dato luogo. Quindi l'ora dell'alba, del mezzogiorno, del tramonto e della mezzanotte astrologici sono ricalcolati rispetto all'ora convenzionale per riportare la rappresentazione dell'ora alla sua fondamentale realtà geografica e astronomica. Per fare questo, tienamo conto della longitudine e della latitudine geografica. Vediamo in breve come fare, visto che siamo programmatori, partendo dall'ora segnata dal nostro PC.

Ora del computer

La scheda madre del vostro computer è dotata di circuiti elettronici oscillanti, che riproducono con buona precisione lo scorrere del tempo e, di conseguenza, della rappresentazione dell'ora nel computer. Si imposta nel BIOS, il sistema operativo di base del computer, e rimane stabile salvo che non ci sia la sincronizzazione con un server remoto, o server NTP. In tal caso l'orologio del PC, anche a livello di BIOS, viene sincronizzato con il server remoto.

Ora civile

Il vostro PC riconosce l'ora standard del Paese in cui vi trovate, il cosiddetto locale che è fondamentalmente riconducibile alla timezone della vostra località geografica. Se installate una qualsiasi distribuzione Linux, durante l'installazione vi viene chiesta la conferma della timezone individuata automaticamente (di solito per l'Italia compare Europe/Rome ma qualche volta mi è successo di leggere Europe/Vatican City. L'elenco completo delle timezone è contenuto nel PC, compilato da quello che, in onore del suo primo curatore, viene chiamato Olson database. Se ne volete un esempio e avete, come me, una Debian o una Ubuntu, aprite un terminale e provate a digitare tzselect . Se provo a seguire le istruzioni, viene fuori questo dialogo:

ubuntu@ubuntu:~$ tzselect
Please identify a location so that time zone rules can be set correctly.
Please select a continent, ocean, "coord", or "TZ".
 1) Africa
 2) Americas
 3) Antarctica
 4) Asia
 5) Atlantic Ocean
 6) Australia
 7) Europe
 8) Indian Ocean
 9) Pacific Ocean
10) coord - I want to use geographical coordinates.
11) TZ - I want to specify the time zone using the Posix TZ format.
#? 4
Please select a country whose clocks agree with yours.
 1) Afghanistan    18) Israel      35) Palestine
 2) Armenia    19) Japan      36) Philippines
 3) Azerbaijan    20) Jordan      37) Qatar
 4) Bahrain    21) Kazakhstan     38) Russia
 5) Bangladesh    22) Korea (North)     39) Saudi Arabia
 6) Bhutan    23) Korea (South)     40) Singapore
 7) Brunei    24) Kuwait      41) Sri Lanka
 8) Cambodia    25) Kyrgyzstan     42) Syria
 9) China    26) Laos      43) Taiwan
10) Cyprus    27) Lebanon      44) Tajikistan
11) East Timor    28) Macau      45) Thailand
12) Georgia    29) Malaysia      46) Turkmenistan
13) Hong Kong    30) Mongolia      47) United Arab Emirates
14) India    31) Myanmar (Burma)     48) Uzbekistan
15) Indonesia    32) Nepal      49) Vietnam
16) Iran    33) Oman      50) Yemen
17) Iraq    34) Pakistan
#? 9
Please select one of the following time zone regions.
1) Beijing Time
2) Xinjiang Time
#? 1

The following information has been given:

 China
 Beijing Time

Therefore TZ='Asia/Shanghai' will be used.
Local time is now: Sun Jun 11 17:03:24 CST 2017.
Universal Time is now: Sun Jun 11 09:03:24 UTC 2017.
Is the above information OK?
1) Yes
2) No
#? 

Come si vede dal dialogo, io posso benissimo fingere di trovarmi a Beijiing in Cina, il sistema operativo non discute. Quello che mi interessa di più, pero', non è l'ora cinese ma quella che viene riportata come Universal Time, che, guarda caso, oggi che siamo in piena ora estiva, è esattamente due ore prima di quella del mio ufficio. L'Universal Time è quello che ci serve, in prima istanza, per calcolare l'ora dell'oroscopo.

UTC o Universal Coordinated Time

UTC non è un vero acronimo, ma una sintesi fra definizioni in lingua inglese e francese. Rappresenta l'ora di riferimento e viene ottenuto correggendo il TAI (Tempo Atomico Universale), ottenuto dalla misura delle oscillazioni dell'atomo di Cesio 133 in diversi orologi atomici in giro per il mondo, mediante una correzione periodica.

Da metà Ottocento fino al 1 gennaio 1972 si è utilizzato il meridiano di Greenwich (Regno Unito) come l'origine dei meridiani (meridiano 0) e l'ora corrispondente, misurata come tempo solare medio, veniva chiamata Greenwich Medium Time (GMT). In pratica doveva rappresentare il momento della massima elevazione solare sull'orizzonte, che necessitava tuttavia di una correzione fino a 16 minuti (equazione del tempo). L'ora così definita veniva contata a partire dalla mezzanotte (è stata rinominata per accordo comune come UT1 o tempo universale). L'UTC che abbiamo citato prima è molto vicino al valore di UT1 ma necessita di una correzione prevedibile ma non calcolabile esattamente che si chiama DUT1 e viene dichiarata dall'IERS (International Earth Rotation and Reference Systems Service) in un bollettino periodico che si trova anche sul web.

In pratica, indipendentemente dalla complessità delle operazioni di rettifica e riallineamento dell'ora rispetto ai fenomeni terrestri e celesti, in astrologia usiamo, per le finalità di calcolo della domificazione, l'UTC che, come abbiamo visto, si ottiene facilmente dal nostro computer o dalle librerie dei principali linguaggi, compreso Python che è il linguaggio che ho assunto come riferimento.

Ora locale

Arriviamo quindi finalmente a calcolare un'ora del giorno che sia congruente con la nostra longitudine terrestre. L'operazione, semplicissima, consiste nell'ottenere la longitudine del luogo per cui redigiamo la carta natale, convertirla in ore e sommarla all'UTC se siamo a longitudine EST, sottrarla se siamo a OVEST.

Facciamo un esempio:

Ho sbagliato percorso e mi ritrovo a girare per le strade di Mumbay. Il mio orologio portato dall'Italia segna le 8.15, ma si schiatta del caldo e credo, data la posizione del Sole, che siamo più vicini a mezzogiorno di quanto il mio orologio possa pensare. So che Mumbay si trova a 72.825833333333 di longitudine est e voglio calcolare che ora sarebbe se non ci fosse un'ora standard e dovessi basarmi su un'ora calcolata per il luogo.

Presto fatto: in Italia è estate e ci sono due ore di differenza rispetto all'UTC, quindi posso calcolare un UTC=6.15. La longitudine di Mumbay, calcolata in ore, è 72.82 ecc / 15, considerando 15 gradi corrispondenti a un'ora di orologio. Bene, questo significa che devo aggiungere 4.8550 ore (4h 51' 18") all'UTC per avere l'ora locale di Mumbay.

Quindi:

8.15 ora estiva italiana - 2 ore 6.15 UTC + 4 ore 51 m 11.06 ora locale di Mumbay

In realtà Mumbay usa l'UTC+5:30 che mi darebbe un'ora civile di 11:45, estesa a tutto il fuso orario di riferimento.

Se fossi a Kolkata (Calcutta) l'orologio indicherebbe sempre le 11:45, ma l'ora locale sarebbe un po' diversa, considerata la longitudine (88°22' cioè +5 ore 53'), quindi sarebbero le 6.15+5.53 = 12.09. La differenza di circa 15 gradi in longitudine si esprime, come prevedibile, in circa un'ora di differenza nell'ora locale, a parità di ora civile

.

Spero che questa breve digressione sull'ora e il luogo sia stata di qualche utilità. Nel prossimo post parleremo di domificazione, cioè di rappresentazione di una rappresentazione della suddivisione del cerchio zodiacale alternativa allo zodiaco, e fondamentale per la centratura dell'oroscopo rispetto alla reale posizione del Sole sull'orizzonte locale. A presto.

E' tempo di ricominciare

Due anni di pausa per un blog sono tanti, e non sono passati per mancanza di interesse, quanto per la necessità di esplorare altri aspetti della programmazione e dei software disponibili per il calcolo astronomico che possono risultare utili anche ad un astrologo dilettante come me. Non sono un astrologo professionista, per cui mi limito a trattare di aspetti tecnici del calcolo, rinviando eventuali discussioni di merito sull'astrologia in sè come arte esplorativa del mondo ad un futuro non ancora precisato.

Per fare un po' il punto, prima di riprendere a trattare i temi del software astrologico, desidero riprendere da dove ero arrivato e, possibilmente, andare oltre.

Negli ultimi anni è esploso il mercato consumer degli smartphone e dei single board computer, di cui le schede Raspberry PI sono l'esempio più interessante, per il gran numero di elementi venduti. Rispetto ai classici computer desktop, questi dispositivi presentano dei vantaggi (gli smartphone e i tablet sono portabili ovunque, le SBC sono dispositivi di bassissimo consumo energetico, tanto da poter rimanere accese 24 ore senza gravare sulla bolletta elettrica) e dei potenziali inconvenienti (il linguaggio Python, che ho usato fino all'ultimo post, è scarsamente implementato sugli smartphone, in cui non gira nativamente e soprattutto non usa le librerie grafiche proprie dei sistemi operativi android e IOS, dovendosi quindi, per un buon risultato finale e una buona usabilità, usare i linguaggi propri di questi sistemi operativi (Java e Objective C), la Raspberry PI è invece basata su processore ARM, quindi la compilazione di codice deve considerare la particolarità della piattaforma). In margine alle mie sperimentazioni in Python su ambiente desktop, avro' ogni tanto occasione di parlare di come riesco a implementare funzioni e programmi anche in questi differenti ambienti operativi.

Ritorno un attimo sul tema delle librerie di alta precisione. Dopo aver tentato di implementare le VSOP2013 in ambiente Python e aver discusso delle difficoltà di usare Python per avere una efficienza elevata quando i termini da sommare sono migliaia, ho definitivamente rinunciato a quest'idea. In realtà ho visto ben poche implementazioni di queste librerie, mentre ho constatato che gran parte degli sviluppi del software vanno in altre direzioni. In particolare al posto dell'approccio analitico classico, basato sulle equazioni del moto dei corpi celesti, con cui ho elaborato alcuni software nei post più vecchi di questo blog, si è progressivamente affermata una metodologia basata sull'uso di effemeridi precompilate, la cui interpolazione con metodi numerici produce soluzioni della posizione dei pianeti e dei satelliti che non sono altrimenti risolvibili con le metodologie tradizionali. Non sono da considerare una novità, in quanto le prime effemeridi, elaborate per supportare le prime missioni spaziali, risalgono alla fine degli anni '60.

In questa pagina di Wikipedia è riassunta la storia delle JPL ephemeris. La versione de436t che è la più recente alla data attuale, copre un tempo che va dal 31 dicembre 1549 ald 25 gennaio 2650 e contiene le posizioni e le velocità del Sole, della Luna, dei pianeti da Mercurio a Plutone o dei loro baricentri. Versioni più estese nel tempo o che considerano più fattori astronomici sono presenti nel sito FTP della NASA all'indirizzo ftp://ssd.jpl.nasa.gov/pub/eph/planets/bsp/ da cui sono liberamente scaricabili. Un inconveniente serio dell'uso di queste effemeridi, rispetto all'approccio analitico classico, è la dimensione, per quelle che possono interessare un astrologo si va da un minimo di 15 mega fino ai 3.4 giga delle de431t. Di cio' dovremo tenere conto nel momento in cui decidiamo di usarle in dispositivi leggeri come smartphone e Single Board.

Oltre ai kernel elaborati dalla Caltech in California per la NASA, esistono altre autorevoli effemeridi prodotte da istituti accademici in Europa, cito quelle del progetto INPOP dell' Institut de Mecanique Celeste et des Calcul des Ephemerides (IMCCE)), che già conosciamo per le VSOP87 e 2013, e le EPM dell' Institute of Applied Astronomy of the Russian Academy of Sciences. A parte lievissime differenze, tutti questi sistemi di integrazione numerica sono sviluppati in modo coerente e presentati in formati compatibili fra loro, rendendo possibile l'utilizzo degli stessi software. Ne discuteremo ampiamente nei prossimi post.

Per riprendere in modo serio la discussione relativamente all'uso di metodologie astronomiche per la pratica astrologica, dobbiamo fermarci a considerare l'elemento tempo, che fin qui abbiamo dato per scontato. Il prossimo post sarà quindi dedicato al tema della misura del tempo in astronomia e quindi in astrologia, utilizzando fonti molto autorevoli. Alla prossima.

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.

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.