--------------------------
#  Kurzeinfuehrung in R  #
--------------------------

# In R gibt es im wesentlichen 4 verschiedene 
# Datentypen: Vektoren, Matrizen, Tabellen oder 
# data frames und Listen. Zahlen sind Vektoren 
# der Laenge 1:
x = 7.7
y = 8.8
x+y


###   VEKTOREN   ###

# Fuer das Anlegen von Vektoren gibt es im wesent-
# lichen 3 Befehle, c(), seq() und rep():
v1 = c(2.22,-4.56,29)    # "c" fuer "concatenate"
v1
v2 = seq(from=10,to=30,by=2)
v2
v3 = seq(from=0,to=2*pi,length=100)
v3
plot(sin(v3))            # R rechnet immer element-weise
plot(sin(v3),type="l")
points(cos(v3))          # addiert Punkte zu einem bestehendem plot
lines(cos(v3),col="red") # dasselbe wie points, nur in linienform
?plot                    # Hilfe-Seiten
v4 = rep(2,10)
v4
v5 = rep(c(77,33),4)
v5
# fuer Vektoren mit Schrittweite 1 (oder -1):
v6 = 1:10
v6
v7 = 10:30    # hat immer Schrittweite 1 (oder -1)
v7
v8 = 4:(-4)
v8
v8[1]
v8[2]
v8[7:9]
v8[7:11]

# viele Funktionen geben automatisch Vekoren zurueck:
x = runif(1000,min=-1,max=1)   # 1000 auf dem Intervall [-1,1] 
                               # gleichverteilte Zufallszahlen
x
plot(x)
plot(x,ylim=c(-2,2))
hist(x)
hist(x,xlim=c(-2,2))

# R ist sehr performant und hat viele Funktionalitaeten 
# fuer das Rechnen mit Zufallszahlen; das Kapitel 5 in 
# RSkript_UniGiessen bietet auf 3 Seiten eine sehr nuetz-
# liche und kompakte Uebersicht. Kann man insbesondere fuer
# Stochastik I + II und Oekonometrie gut brauchen.



###   MATRIZEN   ###

x = 1:20
x
mat1 = matrix(x,nrow=4,ncol=5)
mat1
mat2 = matrix(x,nrow=5,ncol=4)
mat2
mat3 = matrix(x,nrow=5,ncol=4,byrow=TRUE)
mat3
mat3[2,4]
mat3[4:5,2]
mat3[4:5,2:3]

# R rechnet typischerweise immer element-weise:
x
x^2
x+0.33
1/x
mat3
mat3^2   # das ist nicht Matrix-Multiplikation
mat3+0.77
1/mat3   # das ist nicht Matrix-Inverse

# Matrix-Multiplikation:
mat4 = mat1 %*% mat2
mat4
det(mat4)  # nicht invertierbar

v1=c(1,1,1)
v2=c(1,2,3)
v3=c(1,4,9)
mat5 = rbind(v1,v2,v3)   # "rowbind", analog mit cbind, "columnbind"
mat5
det(mat5)
mat5inv = solve(mat5)    # das ist Matrix-Inverse
mat5inv
# check:
mat5 %*% mat5inv
zapsmall(mat5 %*% mat5inv)


##  Eigenwerte und Eigenvektoren einer Matrix  ##

n = 5                                # we repeat this for n=1000
x = runif(n*n,min=-1,max=1)          # n^2 Zufallszahlen
x
mat = matrix(x,nrow=n,ncol=n)        # nxn Zufallsmatrix
mat
eigen(mat)
res = eigen(mat)
str(res)
names(res)

# das Resultat von eigen(mat) 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
n = 1000
# redo above (-> quite fast calculation) and then:
plot(res$values)

# Offensichtlich kann R mit komplexen Zahlen rechnen:
z = 0.5 + sqrt(3)/2 *1i
z
Re(z)
Im(z)
abs(z)
arg(z)
Arg(z)          # R ist case-sensitive
Arg(z)/pi       # 60 Grad

# Neben Vektoren, Matrizen und Listen gibt es noch den Datentyp 
# Tabelle oder data frame. Wenn man Daten aus einem externen file 
# einliest, werden sie typischerweise als data frame angelegt. 
# Das ist, wie eine Matrix, ein 2-dimensionales Daten-Array, deren 
# Elemente mit NameOfArray[i,j] angesprochen werden koennen. 
# Allerdings, waehrend bei einer Matrix alle Elemente denselben 
# Datentyp haben muessen (alles Zahlen oder alles Text-Variablen), 
# kann bei einem data frame der Datentyp von Spalte zu Spalte 
# variieren.



-------------------------
#   LinOpt-Beispiele:   #
-------------------------


# Wir loesen noch einmal die Aufgabe 2 von Uebungsblatt1, 
# und zwar
#
# a) mit Hilfe der graphischen Loesungsmethode, indem wir 
#    die entsprechenden Geraden und Hoehenlinien in R plotten
#    -> wir ueben etwas die Graphik-Befehle
#
# b) mit Hilfe des Simplex-Algorithmus, der bereits in dem 
#    package "linprog" implementiert ist.
#    -> wir muessen das package installieren (d.h. die Software 
#       von einem externen R-Server auf unseren Computer holen, 
#       das muss man nur einmal machen) und dann laden (d.h. 
#       der jeweiligen R-Session mitteilen, dass man diese 
#       zusaetzliche Software jetzt auch benutzen moechte, das 
#       muss man bei jeder R-Session neu machen).
#    Typischerweise hat jedes package eine Bedienungsanleitung,  
#    ein zusaetzliches pdf-file, was wir uns dann also ein bischen 
#    genauer anschauen muessen.

##  a) Graphische Loesungsmethode:  ##

# Wir plotten die Geraden
# y = -x/2 + 85
# y = -x + 150
# y = 60
# fuer x in [0,150] in einem plot-Fenster:

x = 0:150
x
y1 = -x/2 +85
y2 = -x + 150
y3 = 0*x + 60
y3                    # ok, ist Vektor, keine Zahl
plot(x,y1,type="l")
lines(x,y2)
lines(x,y3)

# jetzt brauchen wir noch die Hoehenlinien: 
# 300x + 500y = c  <=>  y = -3/5 * x + c/500
# kurze Erinnerung: das Max war 49000, wir waehlen deshalb 
# das c etwa aus 0,1000,2000,3000,...,58000,59000,60000 

c = seq(from=0,to=60000,by=1000)
c

# wir fuegen jetzt 61 Hoehenlinien hinzu mit einem loop:
for(cc in c)
{
  y = -3/5 * x + cc/500
  lines(x,y,col="red")
}


##  b) mit "linprog"-package:  ##

# Wir muessen zunaechst das package installieren:
# Menue-tab Packages->Install , es erscheint eine Liste mit 
# R-Servern, vermutlich kann man einen beliebigen waehlen, 
# waehlen wir einen aus Deutschland: Goettingen.
# Dann erscheint eine Liste mit saemtlichen zur Verfuegung 
# stehenden R-packages: wir waehlen das linprog package.

# Jetzt muessen wir das linprog package noch fuer unsere 
# aktuelle R-Session laden:

require("linprog")

# Jetzt kann man alle Befehle des packages benutzen. Schauen 
# wir uns die Bedienungsanleitung etwas genauer an:
# Ok, wir muessen die Vektoren b, c und die Matrix A 
# anlegen:

c = c(300,500)
c               # ok, die c()-Funktion funktioniert noch, haben 
                # wir durch das obige c nicht kaputt gemacht..
b = c(170,150,180)
b

v1 = c(1,2)
v2 = c(1,1)
v3 = c(0,3)
A = rbind(v1,v2,v3)
A

# jetzt maximieren wir mit der solveLP()-Funktion aus dem package:

solveLP(c,b,A,TRUE)	# ok, funktioniert!!



# machen wir noch einen kurzen performance-Test:
# 50 Variablen mit 500 Nebenbedingungen: wir nehmen 
# Zufallszahlen fuer b und A:

n = 50
m = 500
a = runif(n*m)	  # Zufallszahlen in [0,1]
A = matrix(a,ncol=n,nrow=m)

c = rep(1,n)      # die Summe x1 + x2 + ... + x50 soll maximiert 
                  # werden

b = runif(m,min=20,max=30)	# Zufallszahlen zwischen 20 und 30

# wir maximieren mit der solveLP()-Funktion:

solveLP(c,b,A,TRUE)	# ..das dauert ein paar Sekunden!
                        # Machen Sie das mal mit Bleistift und Papier...


# wie viele der 500 Ungleichungen sind "am Anschlag"?
# -> wir brauchen die Zahlen in der free-Spalte:
res = solveLP(c,b,A,TRUE)
names(res)
constraints = res$con
constraints
class(constraints)		# ist ein data frame
names(constraints)
free = constraints$free
free
test = ifelse(free<0.0000001,1,0)
test
sum(test)	# mal 34, mal 37, mal 40..