# 6/10/04
# per definire un'array, es 2x3x4:
z<-1:(2*3*4)     # riempiamo il vettore con numeri consecutivi
dim(z)<-c(2,3,4) # lo informiamo delle dimensioni
z                # guardiamolo
dim(z)<-c(6,4)   # le dimensioni possono essere cambiate
z                # riguardiamolo
dim(z)<-c(5,4)   # questo comando da' errore in quanto 5x5 !=  2x3x4

z <- 0           # attenzione: questo comando non azzera le
                 # componenti, ma ridefinisce z !
                 # [ricominciare da capo..o a ridefinislo]
z <- z-z         # cosi' invece si azzera sicuramente.
                 # Oppure, per accedere simultaneamente a tutte 
                 # le componenti:
z[] <- 0         # pone tutte le compionenti a zero

z <- rep(1,6*4)  # crea z con 6x4 '1'
                 # si potevano mettere anche dei NaN: "z <- rep(NaN,6*4)"
dim(z)<-c(6,4)   # fissa le dimensioni
dim(z)<-c(4,6)   # ovviamente la matrice 6x4 non e' la stessa di una 4x6

z[1,]            # prima riga
z[,1]            # prima colonna

# Index arrays (Sect. 5.3 "An introduction to R")
x<-array(1:20, c(4,5))        # il primo argomento e' la lista di valori
			      # il secondo le dimensioni
i<-array(c(1:3,3:1), c(3,2))  # idem per i
                              # risultato va visto...
> i
     [,1] [,2]
[1,]    1    3                # (1,3)
[2,]    2    2                # (2,2)
[3,]    3    1                # (3,1)
> z[i]
[1] 9 6 3                     # z[1,3], z[2,2], z[3,1]

## function array() [Sec. 5.4]
Z<-array(0,c(3,3))            # fa una matrice 3x3 di zeri 
                              # in array(), se la lista di valori e'
			      # inferiore alle dimensioni, essa viene
			      # 'ciclata', altrimenti viene troncata, 
                              # esempi.: 
Z<-array(1:3,c(3,3))
Z<-array(1:2,c(3,3))
Z<-array(1:4,c(3,3))

dim(Z)                        # fornisce un vettori con le dimensioni di Z
length(dim(Z))
Z[] <- 2                      # pone a 2 tutti gli elementi
Z <- Z + 1                    # aggiunge 1 a tutti gli elementi
                              # MA 
Z <- 2                        # !!!   resetta Z ****

##  5.5 : outer product of two arrays
a <- 1:2
b <- 3:5
c <- a %o% b      # l'outer product da' a_i*b_j
                  #  a[1]*b[1], a[1]*b[2], a[1]*b[3]
                  #  a[2]*b[1], a[2]*b[2], a[2]*b[3]
                  #  a[3]*b[1], a[3]*b[2], a[3]*b[3]
> c
     [,1] [,2] [,3]
[1,]    3    4    5
[2,]    6    8   10
 
# in generale:
c <- outer(a,b,"*")     # per il 'prodotto'
                        # che pero' puo' essere sostituito con altri
     			# operatori: nota la differenza fra
c <- outer(a,b,"-")     # e
c <- outer(b,a,"-")
                        # inoltre accetta, piu' in generale, una
			# funzione al posto dell'operatore:
f<-function(x,y) x^2+y
c <- outer(b,a,f)
> c
     [,1] [,2]
[1,]   10   11
[2,]   17   18
[3,]   26   27

# esempio sulla distribuzione dei 10000 determinanti di una matrice
# 2x2, i cui elementi possono assumere valori discreti fra 0 e 9
# [vedi pag. 23, 5.5]
e <-0:9                # valori possibili degli elementi
p <- outer(e,e)        # i 100 possibili prodotti (matrice 10x10)            
d <- outer(p, p, "-")  # le 10000 possibili differenze 
                       # dei prodotti, ovvero i determ (matrice 10x10x10x10)   
dim(p); length(p)      # per verificare
dim(d); length(d) 
fr<-table(d)           # fa una tabella di frequenze dei valori di d
                       # che razza di oggetto e'? -> si ricorda del 
                       # nome 'd' e dei valori, ma apparentemente 
                       # e' un oggetto lungo 163, di dimensione 163
                       #: ECCO: 
names(fr)              # da' le occorrenze di d, viste come nomi
fr                     # da' la frequenza di ciascuna occorrenza
as.numeric(names(fr))  # traduce in numeri le occorrenze (probabilmente
                       # table e' pensata per statistici, per lavorare
                       # anche su attributi non numerici:
# [ provare con 
x<-c("a","b","c","b","c","c","b","c")
table(x)
# ]               OK!
                       # quindi: 
n<-as.numeric(names(fr))
plot(n,fr,type="h")    
                       # ovvero, come sul documento:
plot(as.numeric(names(fr)),fr,type="h",xlab="Determinant",ylab="Frequency")

# in outer si puo' usare una funzione qualsiasi: es.
> ff<-function(x1,x2) x1 + 2*x2 + 3
> outer(1:3,1:3,"ff")
     [,1] [,2] [,3]
[1,]    6    8   10
[2,]    7    9   11
[3,]    8   10   12
# e quindi
> one <- function(x1,x2) 1
> outer(1:2,1:2,"one")
     [,1] [,2]
[1,]    1    1
[2,]    1    1

##################################################################
# 5.6 permutazione degli indici di un'array.
                     # data 
z<-1:(2*3*4)         # con
dim(z)<-c(2,3,4) 
                     # la funzione aperm() scambia gli indici
                     # le modo indicato dal secondo argomento
w<-aperm(z,c(1,2,3)) # : non fa niente: indici: 1->1, 2->2, 3->3
w<-aperm(z,c(3,2,1)) # : indici 1->3, 2->2, 3->1
w[3,3,2]==z[2,3,3]   # verifica (risultato TRUE)
w[1,3,2]==z[2,3,1]

v<-aperm(z,c(2,3,1)) # rotazione ciclica degli indici
v[2,3,1]==z[1,2,3]   # verifica
v[3,4,2]==z[2,3,4]

# nel caso di array bidimensionali, ovvero matrici:
dim(z)<-c(6,4)
aperm(z,c(2,1))      # corrisponde a fare la trasposta, operazione che
                     # si puo' anche fare con t()
aperm(z,c(2,1))==t(z)# verifica


# 5.7 matrici, ovvero array 2x2
A<-array(1:6,c(2,3)
B<-array(1,c(2,3))

# alternativa
A<-matrix(1:6,2,3)
B<-matrix(1,2,3)   # e quindi per creare una matrice vuota 
                   # (ovvero piena di NA):
matrix(,2,3)

# operazioni elemento ad elemento: 
A+B
A-B
A*B   # etc.

# prodotto di matrici: %*%
A%*%B     # da' errore in quanto non sono matrici quadrate

> B%*%t(A)
     [,1] [,2]
[1,]    9   12
[2,]    9   12

> t(A)%*%B
     [,1] [,2] [,3]
[1,]    3    3    3
[2,]    7    7    7
[3,]   11   11   11
                         # quest'ultima operazione si puo' effettuare
                         # con la function crossprod:
t(A)%*%B==crossprod(A,B) # verifica
C<-crossprod(A,B)        #
diag(C)                  # diagonale; funziona anche con matrici non
                         # quadrate:
diag(A)
diag(C) <- c(2,2,2)      # diag() non e' una banale funzione tipo
                         # funzioni fortran o C !!

# ATTENZIONE!
> k<-2
> diag(k)
     [,1] [,2]
[1,]    1    0
[2,]    0    1           # a cosa puo' essere utile? venuto per caso?

solve(C)                  # matrice inversa
solve(C)%*%C              # verifica      
solve(C)%*%C-C%*%solve(C) # altra verifica (gli zeri sono c.a O(1.e-16))

# il motivo del nome 'solve' e' che l'inversione e' usualmente
# un byproduct della soluzione di equazioni lineari

b<-c(3,-1,0)     # vettore delle costanti
solve(C,b)       # soluzione del sistema di equazioni lineari
                 # "vettore delle x"
C%*%solve(C,b)   # verifica

# matrici simmetriche
# data
S<-matrix(c(1,2,3,2,4,5,3,5,2),3,3)
ev<-eigen(S)     # risolve il problema degli autovalori, 
ev$val           # autovalori   (ev$val sta per ev$values)
ev$vec           # autovettori  (ev$vec sta per ev$vectors)

ev$vec[,1]       # autovettore nr 1, etc
                 # check:
S %*% ev$vec[,1] 
ev$values[1] * ev$vec[,1]
                 # attenzione: R distingue fra un'array 3x1 da un
		 # vettore di dimensione 3


# 5.7.4 singular value decomposition of a generic matrix:

s<-svd(C)        # ??? che ci facciamo


# costruzione di matrici mettendo insieme righe o colonne
# Sec. 5.8
> cbind(c(1,2,3),c(11,22,33),c(101,102,103))
     [,1] [,2] [,3]
[1,]    1   11  101
[2,]    2   22  102
[3,]    3   33  103
> rbind(c(1,2,3),c(11,22,33),c(101,102,103))
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]   11   22   33
[3,]  101  102  103
# in realta', queste funzioni non si aspettano di avere esattemente
# righe intere o colonne intere, ma concatenano tutto quello che viene
# passato

# per 'sciogliere' un'array in un vettore
as.vector(z)   # oppure, semplicemente
c(z)