Tre numeri a caso

Nel post precedente ho enunciato quello che sembra un problemino abbastanza semplice. Talmente semplice che non è che uno ci stia a pensare più di tanto, perché la soluzione sembra ultrabanale. Ecco il problema:

Estrarre tre numeri reali positivi a caso in modo che la loro somma sia pari a 1. Fare in modo che le triplette di numeri siano distribuite uniformemente, ovvero ognuna delle possibili triplette deve avere la stessa densità di probabilità di essere estratta a caso.

Sembra talmente semplice che di solito uno applica un metodo molto intuitivo e semplice senza stare a pensarci più di tanto:

Si estraggono a caso tre numeri reali positivi tra 0 e 1, siano x1, x2 e x3, con densità di probabilità uniforme (ovvero ogni possibile reale tra 0 e 1 ha la stessa probabilità di uscire). Quindi si sommano, sia s=x1+x2+x3. I tre numeri cercati sono U1=x1/s, U2=x2/s, U3=x3/s.

Semplice ed intuitivo. Peccato che sia sbagliato: facendo così, certe triplette vengono fuori più frequentemente di altre. Come è possibile? E’ quello che si sono chiesti i miei colleghi e co-autori Enrico e Giorgio nel paper citato in fondo a questo post.

Il problema si presenta spesso nelle nostre ricerche: ci occupiamo di algoritmi per la schedulazione dei processi in sistemi in tempo reale, e spesso dobbiamo misurare e comparare le loro performance. Per fare queste comparazioni si generano a caso dei processi e si cerca di applicare un test di schedulabilità per ogni algoritmo. Spesso, il parametro più importante di un processo per noi è il carico di lavoro che impone sulla CPU. Il carico si esprime come una percentuale, quindi un numero reale tra 0 e 1. Inoltre, una CPU non puà essere caricata per più del 100%, per cui ecco che salta fuori il problema: dobbiamo generare a caso 3 (o più) numeri, uno per ogni processore, in modo che la loro somma sia un valore predeterminato e fisso non superiore a uno. Per semplicità, assumiamo qui che tale valore deve essere esattamente pari a 1.

Enrico stava appunto cercando di rifare degli esperimenti comparativi fra i due algoritmi più famosi: Fixed Priority e Earliest Deadine First, e guardando i lavori in letteratura si accorse che tutti utilizzavano il metodo di cui sopra. Al mio amico venne il sospetto che tale metodo non fosse “fair”, cioè tendesse a favorire uno dei due algoritmi contendendi. Per cui si mise a studiare il problema. Per prima cosa, generò un buon numero di triplette con il metodo della normalizzazione. Tali numeri si possono graficare in 3D, e vanno tutti a finire su un triangolo che ha i sui vertici nei punti (1,0,0), (0,1,0), (0,0,1). Ecco il triangolo che venne fuori:

Come vedete con questo metodo i punti tendono a finire più spesso verso il centro del triangolo. Come mai? C’è una spiegazione geometrica semplice: Guardate la seguente sezione in 2D (cioè, cosa succederebbe se i numeri da estrarre fossero solo 2):

I punti x1 e x2 generati a caso tra 0 e 1 stanno nel quadrato. Dopo la normalizzazione U1 e U2 vengono riportati sulla retta retta obliqua che va da (0,1) a (1,0). Quindi, tutti i punti gialli andranno sul segmento in alto; tutti i punti nell’area rossa andranno sul segmento nel mezzo. Poiché l’area rossa è molto maggiore dell’area gialla, ne consegue che molti più punti finiranno sul segmento centrale rispetto a quelli che finiscono sui segmenti laterali. Torna? Beh, ora immaginatelo in 3D.

Qualcun’altro ha proposto questo metodo:

genero solo due numeri: U1 come distribuzione uniforme tra (0,1); U2 come distribuzione uniforme tra (0, 1-x1). Infine, U3 = 1-U1-U2.

Purtroppo non va meglio, ecco cosa esce fuori:

Stanno più frequentemente da una parte! Questo perché, dopo aver generato il primo, il secondo non ha distribuzione propriamente uniforme tra 0 e 1…

Un metodo corretto è quello del rejection sampling, proposto da Alb su friendfeed, e da hronir su okpanico. Si generano punti su un rettangolo nel piano 2D che ha base e altezza pari alla base del triangolo; se il punto finisce fuori dal triangolo, si butta via e se ne genera un altro. Se invece finsce dentro il triangolo, si proietta sui tre assi per ottenere i tre numeri originari. In questo modo siamo sicur di generare uniformemente. Come riprova, ecco la figura ottenuta:

Bella vero?

C’è un solo problema: se invece di generare solo 3 numeri ne vogliamo generare N, allora il numero di punti che cascano fuori dall’iper-triangolo in N-1 dimensioni sono molti di pù di quelli che cascano dentro. Talmente tanti di più che il metodo diventa improponibile. Quindi, al crescere di N tale metodo non è utilizzabile.

E questo era proprio il problema di Enrico e Giorgio: come generare N numeri a caso tali che la somma faccia un numero fisso minore di 1? Ci sono altri due metodi possibili. Il primo è il seguente:

Generiamo N-1 numeri in maniera uniforme compresi tra (0,1), chiamiamoli x(1), x(2), … x(N-1). Poi li ordiniamo dal più piccolo al più grande, e aggiungiamo x(0) = 0 in testa, e x(N)=1 in fondo. A questo punto, i numeri richiesti saranno:
U(1)=x(1)-x(0)=x(1),
U(2)=x(2)-x(1),

…

U(i) = x(i) – x(i-1),
…

U(N)=x(N)-x(N-1)=1-x(N-1).

Per capire come funziona, pensate di avere un segmento lungo 1, e di “gettare” N-1 punti a caso su di esso. La lunghezza dei segmenti delimitati tra due punti consecutivi sono i numeri che cerchiamo.

Per N=3 si ottiene un grafico 3D proprio uguale a quello mostrato qui su: provare per credere! Inoltre il metodo ha complessità O(N log(N)), quindi va abbastanza bene. Perché funziona? Non è facilissimo fare la dimostrazione. Bisogna dimostrare che la probabilità che un segmento abbia lunghezza minore di un certo L è proporzionale a L. Facendo i conti con le probabilità nel caso N=3 se ne dovrebbe uscire facilmente. L’altra sera l’abbiamo fatto con un mio studente e ce ne siamo venuti via con un paio di paginette… (se qualcuno ha tirato fuori un ragionamento migliore ce lo faccia sapere!)

L’ultimo metodo lo prendo dall’articolo, Enrico e Giorgio l’hanno chiamato UUniFast.

La distribuzione di probabilità della somma di N variabili casuali è la convoluzione delle distribuzioni delle variabili sommate. La pdf uniforme è una retta di inclinazione 1: . La convoluzione di N variabili uniformi è un poliomio di grado N: . Si può quindi sfruttare questo fatto come segue:

• Per prima cosa si genera un numero tra 0 e 1 uniformemente, e se ne fa la radice (N-1)-esima. Chiamiamo questo numero s(N-1).
• Quindi si genera un secondo mumero tra 0 e s(N-1), e se ne fa la radice (N-2)-esima. Chiamiamo questo numero s(N-2).
• E così via fino a s(1).

I numeri richiesti sono quindi:

U(1)=S(1),

U(2)=S(2)-S(1),
…
U(i)=S(i)-S(i-1),

…
U(N)=1-S(N-1).

Ok, panico, se qui non mi avete seguito è perché magari non vi ricordate bene la teoria delle probabilità! Comunque, tutti i dettagli sono nel paper qui sotto.

A proposito: Enrico e Giorgio dimostrarono che il metodo “normalizzato” era a favore dell’algoritmo EDF, e quindi rifacendo le simulazioni con il nuovo metodo di generazione ristabilirono la giustizia tra gli algoritmi! (vince sempre EDF, ma con un distacco minore…)

---

Bini, E., & Buttazzo, G. (2005). Measuring the Performance of Schedulability Tests Real-Time Systems, 30 (1-2), 129-154 DOI: 10.1007/s11241-005-0507-9

Problema: numeri a caso

Tempo di quiz! Ecco un problemino per chi vuole cimentarsi nel week-end (la soluzione la darò lunedì sera).

Voglio un metodo per generare 3 numeri (non negativi) a caso in maniera che la loro somma sia 1. Attenzione: voglio una distribuzione uniforme delle triplette: non ci deve essere una tripletta più probabile di un altra. Sapete fornirmi un metodo? E sapreste generalizzarlo a N numeri?

Mettete pure le vostre soluzioni nei commenti!

Errare (non) è (solo) umano

La calcolatrice elettronica pisana al museo.

I precedenti post del collega Juhan e di Zar mi hanno ricordato i bei vecchi tempi di quando ero uno studente pischello a Ingegneria a Pisa (nei lontanissimi anni ’90) e tra le visite a piazza dei Miracoli in cerca di belle turiste straniere (sempre andate male) e rovinose partite di calcetto (la nostra squadra si chiamava “Inseminators”, ma non finì il campionato tanto bene), ero costretto a studiare il corso di “Calcolo Numerico”. Lo teneva il prof. Lombardi, uno di quelli che  partecipò alla costruzione della mitica CEP (Calcolatrice Elettronica Pisana, ma guardate anche qui), che potete oggi ammirare in un museo sempre in quel di Pisa. I tempi del prof. Lombardi dovevano essere migliori dei miei, perché l’Italia era all’avanguardia nella ricerca sull’elettronica e l’informatica, l’Olivetti di Adriano Olivetti era una gran bella azienda, e gli italiani erano pieni di belle speranze. Come avrete capito il prof. Lombardi non era giovanissimo, e si perdeva spesso (almeno una volta a lezione) nei ricordi del bel tempo che fu. Ma non divaghiamo troppo.

Il corso di Calcolo Numerico, tenuto ancora oggi ed obbligatorio in un certo numero di corsi di laurea a Ingegneria, insegna metodi numerici per il calcolo: per esempio a come risolvere numericamente equazioni, sistemi di equazioni, fare integrali, ecc. Uno dei primi argomenti è il trattamento dell’errore di rappresentazione nei numeri.  Lungi da me rifare qui il corso: vorrei solo farvi notare alcune cosine, che magari vi faranno capire come non sia poi così facile far fare i calcoli a un computer. Dato che non mi ricordo più tanto bene, magari dirò qualche cappellata e allora mi correggerete (spero).

Come tutti saprete, non tutti i numeri reali sono rappresentabili con esattezza. In particolare, ogni numero all’interno del nostro amato PC è rappresentato con un numero finito di bit. La rappresentazione standard per numeri con la virgola è la rappresentazione in virgola mobile. Praticamente ogni numero viene rappresentato come una mantissa e un esponente. Ad esempio, se devo rappresentare il numero 0.00002 il ‘puter lo rappresenta come . Un certo numero di cifre sono riservate alla mantissa, il resto sono riservate all’esponente. Per fare gli esempi semplici, supporrò che 5 cifre decimali siano dedicate alla mantissa, e 2 cifre decimali all’esponente. Si può quindi arrivare in valore assoluto da a .

Quindi, per seguire l’esempio classico, se devo rappresentare posso farlo solo con un certo errore, nel caso in questione l’errore sarà . Qual’è l’errore massimo di rappresentazione che posso avere su un qualunque numero? E’ chiaro che se il numero che devo rappresentare non è troppo grande (cioè maggiore di ) né troppo piccolo (cioè minore di ), l’errore si ha solo sulla mantissa. Cioè, se indico con la rappresentazione del numero , si ha l’errore assoluto

.

Adesso scrivo x e r come mantissa ed esponente, utilizzando lo stesso esponente p:

e ,

dove ho indicato con la mantissa esatta e con la sua rappresentazione. Poi, invece di calcolare l’errore assoluto, calcolo l’errore relativo, cioè in percentuale sul numero:

.

Supponendo che facciate i calcoli facendo un’arrotondamento (e non un troncamento), la differenza dipenderà da quante cifre utilizzate per rappresentare la mantissa; nel caso ad esempio di 5 cifre, tale errore al massimo sarà pari a , e nel caso peggiore ( più grande possibile) avremo appunto . Quindi, possiamo scrivere che:

,

ovvero l’errore massimo relativo che posso fare nel rappresentare il numero è appunto pari a .

Fin qui è stata matematica noiosa e con un sacco di simboli (e probabilmente a questo punto il 99% dei lettori avrà smesso di leggere, e saremo rimasti in pochi, al limite sarò rimasto da solo. Va beh, meglio pochi ma buoni!). Era roba preparatoria, adesso però viene il bello: ci resta da vedere cosa succede quando facciamo dei calcoli. Analizzerò due operazioni: somma e moltiplicazione.

Supponiamo di dover sommare e . Le loro rappresentazioni sono e . Naturalmente, quando facciamo i calcoli sul PC, invece di sommare , sommiamo . Qual’è il massimo errore relativo sulla somma? Vediamo di fare i calcoli per bene, e consideriamo prima il caso in cui x e y abbiano lo stesso segno (per esempio positivo):

Ovvero, l’errore relativo sulla somma è lo stesso che l’errore sui termini originali. Bene, siamo confortati, possiamo sommare numeri senza paura di perdere precisione! O no? Aspettate, ci resta da vedere l’altro caso. Supponiamo che due termini x e y hanno segno diverso. Per semplificare i calcoli  supponiamo il primo sia positivo e il secondo negativo, e inoltre che il primo sia maggiore del secondo. Ecco come calcolare l’errore:

Purtroppo è successa una cosa non molto simpatica. Se x e y hanno segno diverso e sono simili fra di loro, l’errore relativo può aumentare e di molto. Infatti, al denominatore avremo una quantità piccola mentre al numeratore una quantità grande (ricordatevi che x è positivo e y è negativo!). Ad esempio: se prendiamo i numeri 5.02 e -5.01, ognuno dei quali con un errore assoluto massimo di 0.00001, la somma sarà 0.1, con un errore massimo di 0.00002. Però 0.1 è adesso circa 500 volte più grande di 5.01, e quindi l’errore è relativamente più grande sul risultato di quanto non fosse sui numeri originali!

Insomma, non si può neanche sommare in pace! Ma se ci pensate è giusto che sia così. Se sottraete due rappresentazioni uguali, come si fa a sapere se la differenza originale fa zero oppure no? In percentuale potreste anche fare il 100% o il 1000% di errore! Va bene, basta saperlo, e tenerne conto opportunamente quando andate a fare i calcoli.

Vediamo adesso la moltiplicazione e speriamo ci vada meglio. L’errore si esprime come:

Adesso proviamo a fare un po’ di magie algebriche:

Osserviamo il primo termine della somma a numeratore: è il prodotto degli errori assoluti fra i due numeri. Se ad esempio gli errori sono entrambi , il loro prodotto è molto più piccolo: . Insomma, fidatevi, lo possiamo trascurare. Quindi:

.

Insomma, l’errore relativo risultante è la somma dei due errori relativi. Qundi, ogni volta che moltiplichiamo due numeri con lo stesso errore relativo, l’errore risultante raddoppia!

Immaginate di avere un numero con un errore relativo di (ovvero circa , e immaginate di elevarlo alla ventesima potenza (come nell’equazione di Zar). L’errore relativo sarà circa 20 volte superiore ovvero , ovvero, ovvero circa . Ci siamo mangiati un po’ di precisione, ma neanche poi tanta.

Voglio mettere un attimo a fuoco queste cose: non si scappa da questi problemi. Essi hanno a che fare con due fatti non facilmente eliminabili:

1) la precisione nella rappresentazione dei numeri reali all’interno dei calcolatori

2) proprietà matematiche indipendenti dall’esistenza o meno dei calcolatori.

In parole povere: gli ingegneri devono rassegnarsi: se hanno da fare un calcolo complicato devono analizzare l’algoritmo cercando di eliminare il più possibile le operazioni inutili, operando semplificazioni a priori, usando formule invece di calcoli iterativi ogni volta che è possibile. Insomma, non possiamo dare niente per scontato, neanche avendo a disposizione una bella e potente macchina di calcolo.

(siete arrivati fino a qui? complimenti per la pazienza!)

L’alternativa ortodossa a calc — bc

no! non sto parlando di questo

Credevate che il post su calc della scorsa settimana fosse l’ultimo della serie? Allora credevate male perché c’è bc che, a differenza di calc, viene installato automaticamente con Linux (e con Unix, naturalmente).
L’autore di bc è Mark Hopkins (almeno secondo questa pagina), il man di Ubuntu da invece Philip A. Nelson.

bc è un vero e proprio linguaggio di programmazione. Il manuale completo al solito è disponibile con

man bc

o online, All’indirizzo http://man.cx/bc dove potete trovarne un’altra versione in inglese, spagnolo, finlandese, ungherese, giapponese, coreano, polacco e turco –no l’italiano non c’è– diversa da quella installata su Ubuntu. Inoltre le versioni nelle varie lingue sono diverse (OK, non ne conosco nessuna di quelle elencate, ho semplicemente contato le opzioni in synopsis).

L’uso più elementare di bc è un pochino quirkoso

cioè bisogna passargli l’operazione da eseguire attraverso echo, usando la pipe. E ancora

bc liscio determina la precisione (le cifre dopo il separatore decimale) in base agli operandi che riceve; è possibile modificare questo comportamento settando la variabile scale. Notare che in questo caso dobbiamo racchiudere la stringa da echoare tra doppi apici.

È naturalmente possibile salvare il risultato in una variabile, così

o, seguendo la sintassi antica

che pare quella ancora preferita da qualcuno. A proposito: state seguendo il corso di BASH degli amici Bit3Lux e Lightuono, qui o qui vero?

Naturalmente si vede subito che questo modo d’uso è scomodo se si vogliono eseguire più calcoli in sequenza, eventualmente salvando i risultati in variabili. Ecco un esempio

Ma si può fare di meglio: con l’opzione -l (o l’equivalente lunga –mathlib) scale viene automaticamente impostato a 20 e sono disponibili le funzioni standard, come a – arcotangente dell’angolo, espresso in radianti.

Copiando da man.cx ecco la definizione di una funzione, mye, valore approssimato dell’esponenziale e suo esempio d’uso. È fatta per scopi puramente didattici in quanto la funzione e – esponenziale è predefinita quando si specifica l’opzione -l. Ecco il file test_e.bc

scale = 20
define mye(x) {
	auto a, b, c, i, s
	a = 1
	b = 1
	s = 1
	for (i = 1; 1 == 1; i++) {
		a = a*x
		b = b*i
		c = a/b
		if (c == 0) {
			return(s)
		}
		s = s+c
	}
}

/* calcolo l'esponenziale dei
   primi dieci interi
*/

for (i = 0; i <= 10; i++) {
   mye(i)
}

quit

eseguendola otteniamo

E con questo credo di aver illustrato cosa si può fare con bc.
Storicamente –ma non in Linux– bc si appoggia su dc.
Già resta da parlare di dc, prossimamente, forse 8)

Grafici con Scilab

Scilab è ottimo. Per tante cose. Io l’ho scoperto per i grafici, lo sto utilizzando e in questo post vorrei riassumere cosa ho scoperto.
Intanto è ottimo, lo uso e lo consiglio a tutti.

Pero non è facilissimo, è un prodotto grosso (proteiforme), si evolve (giustamente) e qualche istruzione diventa obsoleta (com’è giusto che sia).
Non sempre i manuali sono aggiornati e per quel che ho trovato manca un “Getting Started” o un “How To” per le singole branche in cui può essere suddiviso. In particolare per il plotting ho dovuto scorrere bazillioni di files .html e .pdf e in nessuno ho trovato tutto quello che cercavo.

Tra i tanti considero fondamentali:

1. Scilab: Graphics di Anuradha Amrutkar dell’Indian Institute Of Technology, Bombay, 15th April, 2010, disponibile qui;
2. Graphics: Getting started, qui e naturalmente qui cercando Graphics Library (e c’è tutto il resto del manuale online, da perdersi).

Allora ecco cosa ho ottenuto. Il mio sogno sarebbe stato quello di costruire uno script che, lanciato da terminale, generasse automaticamente il grafico richiesto; non è esattamente quello che sono riuscito a fare, difatti lanciando l’interprete di Scilab ottengo:
e, ovviamente, si apre il terminale di Scilab. Trascinando uno script, qualunque script, su quest’ultimo mi compare questo messaggio d’errore:
Ci sono due cose: 1) un messaggio solo parzialmente visibile, che fa riferimento a ecanvas"* che non so cosa sia, e 2) la variabile mostra non definita. Strano che c’entri una variabile con il nome italiano, non citata da nessuna parte.
Ho provato a inserire istruzioni relative a driver(), xinit() e usecanvas() senza risultato. Per mostra basta settarla a %t o %f (non so qualse sia quella giusta).

In genere la prima esecuzione dello script abortisce, completamente o produce solo parte del grafico, basta ridroppare lo script una seconda volta e tutto funziona.

Questi sono alcuni esempi, presi dalla documentazione citata.

mostra = %T;
t=linspace(0,%pi,20);
a=gca();a.data_bounds=[t(1) -1.8;t($) 1.8];

plot2d(t,[cos(t'),cos(2*t'),cos(3*t')],[-5,2 3]);
e=gce();
e1=e.children(1);e1.thickness=2;e1.polyline_style=4;e1.arrow_size_factor = 1/2;
e2=e.children(2);e2.line_style=4;e2.thickness=3;
e3=e.children(3);e3.line_mode='on';e3.mark_background=5;

hl=legend(['cos(t)';'cos(2*t)';'cos(3*t)']);

xs2png(0, "e4.png")

mostra = %F;
t = linspace(-%pi, %pi,30);
function z = mysurface(x,y);
    z = x * sin(x) ^ 2 * cos(y);
endfunction
contour(t,t,mysurface,10)
xs2png(0, "aa2.png")

mostra = %F;
x = [0 : 0.1 : 2 * %pi]';
plot2d(x, [sin(x) sin(2*x) sin(3*x)], ..
       [1,2,3], leg = "L1@L2@L3", nax = [2,10,2,10],rect = [0,-2,2*%pi,2]);
xs2png(0, "scp4.png")

mostra = %f;
t=-%pi:0.1:%pi; m=sin(t)'*cos(t);
clf()
xset("colormap",jetcolormap(64))
colorbar(-1,1)
Sgrayplot(t,t,m, strf="041")
xs2png(0, "scp6.png")

mostra = %f;
n = 30;
nt = 100;
x = linspace(0,2*%pi,n);
y = linspace(0,%pi,n/2);
z = sin(x')*sin(y);
t = linspace(0,4*%pi,nt);
show_window(); clf()
f=gcf();
f.color_map=jetcolormap(64);
f.pixmap='on';
colorbar(-1,1)
Sgrayplot(x,y,cos(t(1))*z, strf="042", zminmax=[-1,1])
c=gce(),e=c.children
xtitle("Kaa''s eyes")
for i = 1:nt
  e.data(:,3)=matrix(cos(t(i))*z,-1,1);
  show_pixmap()
end
f.pixmap='off';
xs2png(0, "scp7.png")

(quest’ultima è dinamica, provate…).

Bon. Per adesso Scilab lo uso così.
Ho trascurato altre cose, per esempio lo sapete che è uscita la nuova versione di newLISP, la 10.3? Lo sapevo che non interessa nessuno
Ho trascurato anche Python, ma restate sintonizzati, prossimamente…
E poi tante altre cose nuove 8) e nuovi collaboratori 8)

Scilab, le funzioni e tutte le altre meraviglie

In questo post, continuazione del precedente,si coninua l’esame del pdf progscilab per un illustrazione veloce delle notevoli capacità del package. È utile per arrivare a una mia conclusione, giù in fondo in “Conclusione”

In Scilab le funzioni sono estremamente flessibili. Intanto si distinguono in macro, scritte nel linguaggio di Scilab e primitive, scritte in qualche altro linguaggio (per esempio C, C++, Fortran) e compilate.
Nella maggior parte dei casi non c’è differenza tra i due tipi e per questo si usa il termine generico function. Ovviamente la cosa cambia nel debug.
Le funzioni sono variabili: non male per un progetto partito avendo come riferimento il Fortran
D’altronde oggi è una cosa abbastanza comune, a me capita di farla anche senza volerlo: in Python basta dimenticare le parentesi. Ma è di una comodità unica, conosciuta come callback.
A p.57 e ssg viene spiegato molto bene, con esempi, come creare funzioni con numero di argomenti variabili, sia in input che in output, cose tecniche: se servono leggetevele e sapetele
A p.81 si parla dello scope delle variabili. Diversi casi sono possibili, attenzione! La soluzione più semplice (il livello dello stack) di p.82 può creare problemi, raccontati nelle pagine successive. A me casi del genere sono capitati più volte, provo una discreta Schadenfreude venire a conoscenza che succede anche a altri, tanto spesso che vale la pena di riportarli in un manuale 8)
A p.93 e ssg viene descritto come sia possibile generare dinamicamente le funzioni come stringhe cosa che conferisce a Scilab un livello di flessibilità altrimenti irraggiungibile (ma a me comincia a frullarmi in t’a capa il Lisp).
Le tecniche per rendere l’esecuzione veloce sono descritte a partire da p.99, Performances.
Da copiare l’idea della coppia di funzioni gemelle tic() e toc(): la prima fa partire il cronometro, la seconda lo blocca. In uno dei prossimi post le faccio in Python, newLISP e Fortran, forse, se nel frattempo non me ne sono dimenticato o sono stato attratto da qualcos’altro.
Bello l’esempio d’uso del profiling, da p.102 (profilaggio? la banda del nano usa “briffare” addirittura!).
A p.110 inizia la vettorizzazione, da leggere. Poi vengono i trucchi per ottimizzare il codice, p.115. Regola aurea: usare le funzioni di Scilab e evitare come la peste i cicli for.
Poi basta. Il discorso diventa troppo tecnico e mi sono addormentato

Il pdf 419922 An Introduction to Scilab (un’altra) di Satish Annigeri
(prof. per ingegneri civili) è conciso (55 pagine) e pratico: a p.35 ci dice come gestire file su disco, nella pagina successiva già si occupa di plottaggio di curve polinomiali e poi passa a leggere file di Excel.
A p.47 passa a un caso pratico da ingegneri civili, calcolo di una struttura reticolare piana arrivando a determinare le sollecitazioni nelle aste. Notevole.

Il pdf intro dello Scilab group a p.87 e ssg ci illustra come gestire i grafici.
L’esempio di p.81 mi da errore e allora ho provato a dare le istruzioni direttamente dentro al Terminale si Scilab. Comodissimo il comando xset() che apre la finestra grafica e permette di modificare il grafico. Ecco il mio risultato:

Proseguendo con gli esempi non tutti vengono, ci sono errori nei comandi, alcuni di difficile individuazione, i font sono diversi da quelli del manuale (Windoze strike again!), …
Ci sono però numerosi file di esempi, richiamabili direttamente anche all’interno dell’editor (alle volte sono davvero gnugnu).

Il capitolo 6 a p.101 ci introduce a una delle caratteristiche più avanzate di Scilab: Interfacing C or Fortran programs with Scilab. Chiaramente questo sarà usato solo da utenti skillati per task non banali ma sembra fattibile, semplice e ci sono esempi. No non ci provo, sono troppo vecchio per quest cose. Però: avete nai provato a farlo con M$ DOS (sto parlando quando ancora non si usava Windoze): frustrante. Molto meglio la filosofia *NIX/Linux.

Per la parte matematica sono naturalmente previsti collegamenti con Matlab e Maple.

Conclusione

Il manuale Programming in Scilab di Michaël Baudin, contenuto nel file progscilab-so-v0.6.pdf è fondamentale.
Riassumerlo è impossibile 8)
Ciò porta a una conclusione: questi post su Scilab non possono sostituire il manuale. Per diversi motivi:

• il manuale è ottimo;
• Scilab è un tool molto grande e completo, il suo linguaggio di programmazione anche;
• compito del blog non è fare un corso bensì segnalare le chicche che trovo, anzi che mi piacciono;
• quando possibile adoro proporre script che per questioni oggettive debbono essere piccoli, se piacciono sarà il lettore a completarli.

Prossimamente

Il modulo Plot2d mi attira. Ma ci sono anche altre cose nel mondo della programmazione. È da parecchio che non nomino il newLISP, Pyhton idem. Poi ci sarebbero altri linguaggi, non quelli predominanti C++ e Java, su quelli c’è già troppa roba, altri meno noti e pubblicizati ma intriganti. E poi ci sarebbe da considerare anche gli utenti di Windows: alcune cosa che con Linux si danno per scontate lì mancano, o co$tano.
Prossimamente, se non mi sto beccando l’influenza, sempre su questo canale

Scilab – polinomi e matrici

Dopo l’intro e aver visto come usare Scilab per fare un grafico oggi racconterò di come usarlo per fare calcoli. Scilab adora contare e allora via con i calcoli.

Intanto sto ancora imparando ma è disponibile tutto quello che serve, anzi molto di più. Parecchi documenti si trovano a Scilab Wiki. In particolare Tutorials and documentations dove c’è, tra un mucchio di altre cose, Programming in Scilab di Michaël Baudin.

A p.22 introduce i polinomi. Seguendolo definiamo un polinomio in funzione delle sue radici 1 e 2: p(x) = (x-1)(x-2) con

p = poly([1,2], "x")

Si può definire il polinomio in funzione dei suoi coefficienti (in ordine crescente): q(x) = 1 + 2x sarà

q = poly([1,2], "x", "coeff")

A questo punto i polinomi si possono sommare, sottrarre, moltiplicare:

-->m = p * q
 m  =

              2    3
    2 + x - 5x + 2x   

-->

Per valutare il valore del polinomio in un dato punto si usa horner. Ecco, niente panico!, un po’ di sano terrorismo matematico, qui su wikipedia, ma con Scilab niente panico!

:-->p = poly([1,2], "x")
p  =

          2
2 - 3x + x

-->horner(p, [0 1 2 3])
ans  =

2.    0.    0.    2.

che sono, ovviamente, i valori per x = 0, 1, 2, 3. Quindi dimenticate la pagina di wiki: con Scilab i calcoli (noiosi e (dai) astrusi) li fa lui.

roots ci da le radici, verifica:

-->p = poly([1,2,3,4,5], "x")
 p  =

                     2     3     4   5
  - 120 + 274x - 225x + 85x - 15x + x   

-->roots(p)
 ans  =

    1.
    2.
    3.
    4.
    5.    

OK, viso che funziona

-->p = poly([1,2,3,4,5], "x", "coeff")
 p  =

               2    3    4
    1 + 2x + 3x + 4x + 5x   

-->roots(p)
 ans  =

    0.1378323 + 0.6781544i
    0.1378323 - 0.6781544i
  - 0.5378323 + 0.3582847i
  - 0.5378323 - 0.3582847i  

numeri brutticomplessi, ma ecco:

-->t = poly([24, -50, 35, -10, 1], "x", "coeff")
 t  =

                  2     3   4
    24 - 50x + 35x - 10x + x   

-->roots(t)
 ans  =

    1.
    2.
    3.
    4.  

Scilab è sempre più simpatico
Uno (cioè a dire il vero io, fino a un istante fa) potrebbe pensare che Scilab tratti solo matrici bidimensionali: falso, c’è l’ipermatrice (hypermatrix)

This feature is familiar to Fortran developers, where the multi-dimensional array is one of the basic data structures.

Quindi la volta scorsa ho detto una cosa non vera: chiedo venia ma mi verrebbe voglia di dare la colpa al Fortran

A p.51 di progscilab si inizia a parlare delle funzioni. Potrebbe essere l’argomento di un prossimo post, forse…

Introduzione a Scilab

Scilab è una meraviglia. Molto di più di quello che avevo immaginato, ma c’è sempre un punto di partenza e che io sono n00b (niubbo) di prima categoria non è una novità. Per cui questo post va preso come un primo assaggio, assolutamente parziale e non esaustivo, diciamo un assaggino.

L’installazione di Scilab è davvero all’antica: scompattate il file tar o zip o quello che è (naturalmente è multipiattaforma: Linux, Windows e MacOS) dove vi pare e piace e avrete la subdirectory bin con lo script scilab per aprire il terminale.

È luogo comune trito e ritrito che Ankh-Morpork non è stata costruita in un giorno, ogni cosa richiede tempo. E bisogna ambientarsi. E con calma, arrabbiarsi non serva anzi è controproducente.

Allora io ho scompattato il tar il ~/apps; ho creato nella directory ~/bin il symlink per scilab trovato in bin di Scilab  ed eccomi pronto.
A dire il vero un geek lo installerebbe in /usr ma Ubuntu ti obbliga ogni sei mesi a reinstallare tutto; in questo caso basta backuppare la directory ~/apps 8)

Da quando K&R hanno inondato il mondo (la parte di esso che ci interessa) con “The C Programming Language” si deve partire con “Hello World!“. Non sono solo: il pdf introscilab.pdf scaricato dalla pagina della documentazione mi guida passo-passo. Quindi apro il terminale nella dir dove provo i miei pasticci e lancio il comando

scilab

che apre il Teminale di Scilab nel quale posso sbizzarrirmi come da figura

Visto che tutto è OK posso scrivere quit per uscire. E per oggi basta. No! fermi lì! scherzavo.

Intanto per non bloccare il terminale conviene dare il comando scilab & cosa di cui mi ricordo sempre un secondo dopo aver premuto Invio.

Dal Terminale di Scilab lancio Applicazioni/SciNotes che mi apre l’editor dentro al quale posso scrivere il mio script hello.sce, come in figura

A questo punto posso, dopo aver salvato lo script posso trascinare (draggare) il file da Nautilus alla finestra del Terminale che lo esegue, fin troppo bene: l’istruzione finale quit chiude il terminale stesso! Allora i casi sono due: 1) commentare l’istruzione quit inserrendo // all’inizio della riga in modo che la sua esecuzione nel Teminale (di Scilab) non lo chiuda; oppure 2) lanciare da terminale (quello della shell) il comando scilab -nwni -nb -f hello.sce
Chiaramente preferisco la soluzione 2. Le opzioni sono elencate in fig. 11, p. 20 di introscilab.

il mio desktop è sempre molto affollato

OK, funziona. I file .sce si possono benissimo editare con gedit, vi, Emacs anche se l’editor di Scilab colora la sintassi, ha l’help, etc.

Una delle caratteristiche peculiari di Scilab è il tipo di dati che predilige. Qui le alternative operate dai vari linguaggi sono innumeri: il Forth ha solo gli interi, il Basic preferisce i reali, il Lisp parla di atomi e liste, nel C tutto quello che non è intero è un puntatore (sì, lo so semplifico, dai! non ho mai capito a fondo la cosa ma non mi cazziate che tanto cancellerò i commenti). Scilab usa le matrici, p.31. Inoltre, anche se non detto esplicitamente sono matrici bidimensionali.
Ah se l’avessi avuto 25 anni fa! Quanti cicli mi sarei risparmiato! Inoltre è veloce a trattarle:

a Scilab script which performs the same operations with loops is
typically from 10 to 100 times slower.

Una cosa divertente: la matrice Identità come la chiamiamo? I può generare confusione, sembra una variabile o, peggio ancora, un indice. Pas de panique: la chiamiamo eye, ingénieux! vraiment.

Una cosa che può creare confusione è l’interpretazione come interi di numeri float (p.45): l’interpretazione avviene sempre per troncamento e non per arrotondamento. Da ricordarsene.

L’istruzione for (p.52) usa una sintassi sua quando si precisa il passo; d’altronde capita la stessa cosa con le sequenze. E qui anche Python è particolare e a volte me ne dimentico. Ma il più balzano di tutti era il Fortran di una volta (prima del 1978) che pretendeva di fare comunque sempre un giro anche quando non doveva e si doveva mettere un if di guardia con relativo goto: be tempi, sigh!

Le funzioni (p.55 e ssg) hanno anche loro una sintassi particolare, ma va bene, basta saperlo. A p.59 racconta come creare le proprie librerie e ci dice che inoltre ci sono i moduli che qui non affronta.

A p.66 inizia il capitolo Plotting. La decisione di esaminare Scilab da parte mia è nata da un post di Gianluigi Filippelli, qui. . Poi ho scoperto che Scilab è ben altro. E dice il dr. prof

Scilab viene usato un bel po’ nel mio laboratorio [...] anche come generatore di codice C per sistemini embedded. Ecco ad esempio un’applicazione:
http://www.evidence.eu.com/content/view/361/432/

Ma per fare grafici è ottimo, poi si possono esportare nel formato voluto (p.74).

Ho provato a giochicchiarci ottenendo questo. Ho dovuto usare il manuale per plot.

// esempio di grafico fatto quasi senza guardare il manuale

xt = linspace(0, 10, 11)
c1 = [0, 10, 8, 6, 4, 2, 4, 6, 8, 5, 8]
plot(xt, c1, "r-o")

xr = linspace(5, 10, 6)
c2 = [3, 2, 1, 5, 4, 6]
plot(xr, c2, "b-p")

m = [4, 3, 2, 4, 2, 2, 6, 8, 6, 5, 2]
plot(xt, m, "k:*")

ps = [4, 4, 3, 2, 5, 6, 4, 3, 5, 7, 5]
plot(xt, ps, "m--d")

xtitle("Ultimi indici rilevati", "tempo", "valore indice")
legend("c1", "c2", "m", "ps")
xs2png(0, "indici.png")

Quant’è grosso il manuale? 4001 pagine! Panico, paura, terrore

Aussi en français

Ho provato a usare livret_lycees_2010.pdf

L’objectif de ce livret est de vous guider pas à pas dans la découverte des différentes fonctionnalités de base du logiciel Scilab dans le cadre d’une utilisation en classes de lycées. Cette présentation se limite volontairement à l’essentiel pour permettre une prise en main facilitée de Scilab.

ho (re)imparato diverse parole nuove (crochet, suite), altre le sapevo già (logiciel, taper) ma poi mi sono perso con afficher (== disp) e histogramme (== histogram). Questi francesi! hanno tradotto anche i comandi; cosa vecchia peraltro, è una loro fissa

Conclusione

Scilab dell’INRIA è molto di più, ci sono anche altre cose, ad esempio Xcos. Il sito web è bello, chiaro e esaustivo: un modello. Credo che ne parlerò ancora in futuro.