#--------------------------------------#
#                                      #
#    week1b: Rechnen mit Matrizen;     #
#      Schleifen und Funktionen        #
#                                      #
#--------------------------------------#


Die wichtigsten Matrix-Operationen sind:

 - %*%        fuer Matrix-Multiplikation
 - solve()    fuer Matrix-Inversion
 - t()        Transponieren
 - det()      Determinante
 - eigen()    Eigenwerte und Eigenvektoren
   
Schauen wir uns das an:
   
   
###  Start R-Session  ###

# Generelle Regel fuer das Rechnen 
# mit Vektoren und Matrizen in R:
#
# R versucht immer, ELEMENTWEISE zu rechnen:

x = 1:5
x
1/x 
x - 7
x^2
exp(x)

# fuer Matrizen gilt dasselbe:
A = matrix(1:10,2,5)
A
1/A
A - 7
A^2
A*A                         # das ist keine Matrix-Multiplikation
exp(A)                      # das ist kein Matrix-Exponential, sondern 
                            # elementweise: exp(a11), exp(a12), exp(a13),...

# Damit stellt sich die Frage: Was ist die Syntax 
# fuer Matrix-Multiplikation und das Matrix-Inverse?

AT = t(A)
AT                          # A transponiert

# Matrix-Multiplikation: %*%

M1 = A %*% AT
M1
M2 = AT %*% A
M2

det(M1)                     # Determinante
det(M2)                     # ziemlich nullig..

# Matrix-Inverse: solve()

M1inv = solve(M1)
M1inv
M1inv %*% M1                # Einheitsmatrix
M1 %*% M1inv 
1/M1                        # NICHT das Inverse
M1 * 1/M1
M1 %*% (1/M1) 
M1 %*% 1/M1                 # Fehlermeldung

det(M2)                     # approx 0, nicht invertierbar
solve(M2)                   # Fehlermeldung




# EIGENWERTE und EIGENVEKTOREN

# Wenn Sie etwas groessere Programme in R schreiben wollen, 
# wollen Sie natuerlich nicht Zeile fuer Zeile in das Consolen- 
# Fenster eingeben, wo ja, wenn Sie return gedrueckt haben, die 
# Zeile, der Code, dann verschwunden ist, Sie sehen das dann ja  
# nicht mehr. 
# Also Sie brauchen ein anderes Fenster, ein Programm-Editor, 
# wo Sie erstmal Ihren ganzen Code reinschreiben koennen und 
# dann ganz am Ende wollen Sie das dann executen. 
#
# Dazu, wenn Sie wie ich in Windows mit dem RGui arbeiten, 
# muessen Sie auf 
#                    File->New script
#  
# gehen, und dann oeffnet sich ein neues Fenster. Dort koennen 
# Sie Ihren R-Code reinschreiben, und wenn Sie den Code dann 
# ausfuehren wollen, selektieren Sie die entsprechenden Code-
# Zeilen mit der Maus (also so dass das dann blau hervorgehoben 
# oder blau unterlegt ist) und dann muessen Sie 
#
#                      die Taste F5 
#
# druecken. Dann springen genau die selektierten Code-Zeilen in 
# das R-Console-Fenster und werden dann also ausgefuehrt. 
#
# Oeffnen Sie jetzt ein neues Skript-Fenster, geben Sie den 
# folgenden Code ein und selektieren und executen Sie dann 
# mit der Taste F5:

n = 5                                    # we repeat this for n=1000
x = rnorm(n*n)                           # n^2 Zufallszahlen
x
M = matrix(x,nrow=n,ncol=n)              # n x n Zufallsmatrix
M

# Eigenwerte und Eigenvektoren:
eigen(M)

# wir schreiben das Resultat in eine Variable:
res = eigen(M)

# Info zur Datenstruktur:
mode(res)
str(res)
class(res)

names(res)

# das Resultat von eigen(M) ist eine Liste mit 2 Elementen, 
# res[[1]]=res$values ist vom Typ vector und res[[2]]=res$vectors 
# ist vom Typ Matrix.
# Listen sind im wesentlichen dazu da, um Objekte mit unterschied-
# lichen Datentypen in einem Objekt zusammenfassen zu koennen. 
# So ist der Rueckgabewert von R-Funktionen haeufig vom Typ Liste, 
# damit verschieden Arten von Information in einem Objekt zurueck-
# gegeben werden koennen. 

res[[1]]
res$values

res[[2]]
res$vectors

# folgendes scheint auch zu gehen:
res[1]          # anstatt res[[1]]
res[2]

# Machen wir das jetzt nochmal fuer eine 1000 x 1000 Matrix 
# mit 1 Mio Zufallszahlen:

n = 1000
x = rnorm(n*n)
M = matrix(x,nrow=n,ncol=n)
res = eigen(M)

# das geht noch recht fix, messen wir mal die Zeit:
t0 = Sys.time()
res = eigen(M)
Sys.time()-t0                 # etwa 5 Sekunden auf meinem altem Geraet, 
                              # 0.5 Sekunden mit Intel-MathKernelLibrary

plot(res$values)

# Vermutung: die Eigenwerte sind gleichverteilt in 
# einem Kreis mit Radius sqrt(n).
# Checken wir das: Offensichtlich kann R mit komplexen 
# Zahlen rechnen: Wurzel aus -1 = i wird in R durch die 
# Zeichenkette 1i dargestellt:

z = 1i
z^2       # ok, das passt..

# Jetzt machen wir einen Kreis in der komplexen Ebene 
# mit Radius Wurzel(n):

phi = seq( from=0 , to=2*pi , length=500 )
z = sqrt(n) * exp( 1i * phi )

plot(res$values)
lines( z , col="red" )        # fuegt Punkte zu einem bestehen Plot 
                              # hinzu, im Linien-Format

# probieren wir noch groessere n:
n = 2000
x = rnorm(n*n)
M = matrix(x,nrow=n,ncol=n)
t0 = Sys.time()
res = eigen(M)
Sys.time()-t0                 # etwa 32 Sekunden auf meinem alten Geraet,
                              # 2.5 Sekunden mit Intel-MathKernelLibrary

# wir wollen nur die Eigenwerte wissen, aber nicht die 
# Eigenvektoren:
t0 = Sys.time()
res = eigen(M, only.values=TRUE )
Sys.time()-t0                 # etwa 15 Sekunden auf meinem alten Geraet


# und nochmal das Bild:
plot(res$values)
phi = seq( from=0 , to=2*pi , length=500 )
z = sqrt(n) * exp( 1i * phi )
lines( z , col="red" )            # sehr huebsch..



# was gibt die eigen()-Funktion im Fall 
# von nicht-diagonalisierbaren Matrizen?

test = cbind(c(1,1),c(0,1))
test                              # nicht diagonalisierbar..
eigen(test)



#-----------------#
#    Schleifen    #
#-----------------# 

# Wir wollen eine n x n - Matrix mit Eintraegen a_ij = i*j 
# anlegen: Geht am schnellsten mit dem outer-Befehl:

n = 6 
A = outer(1:n,1:n)
A

# oder aequivalent
calctime = Sys.time()            # Rechenzeit mit angeben..
A = outer(1:n,1:n,FUN="*")
Sys.time() - calctime
A

# Einfachste oder elementarste Methode: Mit Doppel-Loop:
B = matrix(0,n,n)          # Matrix mit Nullen initialisieren
B

calctime = Sys.time()
for(i in 1:n)
{
  for(j in 1:n)
  {
    B[i,j] = i*j
  }
}
Sys.time() - calctime
B                          # geht auch..

# Vergleichen wir die Rechenzeiten fuer grosse n:
n = 10000                          # 100 Mio Matrix-Elemente..
# A: etwa 0.1 sec
# B: etwa 7 sec, loops dauern etwas..




#--------------------------#
#    Benutzerdefinierte    #
#        Funktionen        #
#--------------------------#

MeinQuadrat = function( x )
{
  result = x^2
  return(result)
}

x = 7
MeinQuadrat(7)

y = 5:15
y
MeinQuadrat(y)

M = diag(y)
M
M = M+2
M
MeinQuadrat(M)


TestVermutung = function( n )
{
  x = rnorm(n^2)
  M = matrix(x,n,n)
  evals = eigen(M,only.values=TRUE)[[1]]
  phi = seq(from=0,to=2*pi,length=500)
  z = sqrt(n)*exp(1i*phi)
  plot(evals)
  lines(z,col="red")            # hat kein return, plottet nur was..
}

TestVermutung(1000)
TestVermutung(2000)


f         # Error: object 'f' not found

# ok, das f ist noch nicht verbraucht, koennen wir etwa als eigenen 
# Funktionsnamen verwenden:

f = function(x,y)
{
  res = 1/(1+x^2+3*y^2)
  return(res)
}
f(0,0)
f(0,1)
f(0)   #error-message..

f = function(x,y=1)      # y hat den default-Wert 1
{
  res = 1/(1+x^2+3*y^2)
  return(res)
}

f(0,1)
f(0)        # geht jetzt..

f(3,8) 
f(y=8,x=3)  # geht auch, falls man nicht die genaue Reihenfolge weiss..
f(x=3,y=8)  # dasselbe
f(8,3)      # was anderes