#################################################
### Konstrukte ('Flow Control')
#################################################

# Häufig sollen Ausdrücke (Berechnungen, Schätzungen, Simulationen, Plots...) nur unter bestimmten
# Bedingungen und/oder sehr oft ausgeführt werden.
# Manchmal ist eine Operation auch so komplex, dass man sie in wiederholten Schritten leichter 
# durchdenken und programmieren kann.

# Man braucht dann sog. Konstrukte; hier: 
# - Bedingte Ausführung (if, else)
# - Schleifen (for, while)

# Oder, als schnellere Alternative zur Schleife die Anwendung von Funktionen auf
# viele Elemente gleichzeitig mit ...apply()-Befehlen verschiedener Art
# Vectorize() modifiziert Funktionen, so dass sie auf mehrere Elemente 
# gleichzeitig angewendet werden kann



#################################################
### Bedingte Ausführung
#################################################


### if und else: Bedingte Ausführung
### falls "if" Bedingung wahr ist, wird Ausdruck ausgeführt,
### wenn nicht die Befehle nach "else"

sport   <- TRUE
rauchen <- TRUE


if (!rauchen){
  print("Gut so!")
}  else{
  print("Treib Sport")
}

if(sport & !rauchen) "Weiter so!"

if (sport){
  if (!rauchen) "Weiter so!"
}
# Auch hier wichtig: Mehrere Befehle als EIN Ausdruck mit "{}" 

y <- c(3,5,4,6,9,7,4,2,-0.1)

y>0

if(y>0) "Hierher!" # wertet nur den ersten Eintrag von y aus

any(y<=0) 

!all(y>0)


if (all(y<0)) # all() ist "TRUE", falls Argument="Vektor mit lauter TRUE"
{
    print("Alle y sind größer als 0")
    ly <- log(y)
    plot(ly, main="Logarithmus von y")
} 
else {
    print("Mindestens ein Element von y ist kleiner oder gleich 0!")
    plot(y, main="Werte von y")
}

z <- c(3,5,4,6,9,7,4,2,-0.1)

if(max(z)<=5){
  print("Alle Einträge sind kleiner als 5")
} else {
  print("Ein Eintrag ist größer als 5")
}

### Aufgabe E1
# Benutzen Sie die Funktion is.numeric() und is.character() um  zwei Objekte x und y
# auf einen numerischen und auf einen Buchstabenwert mit Hilfe einer selbstgeschriebenen Funktion
# if_test zu prüfen. 
# Falls sowohl x als auch y dem entsprechen, geben Sie "Super" aus, ansonsten sagen sie,
# ob bei x oder y oder bei beiden der Fehler liegt. 
# Testen Sie mit:
# 1. 
x<-5
y<-"char"
# 2.
x<-"char"
y<-5
# 3. 
x<-list(1,2)
y<-matrix("char")





#################################################
### Schleifen
#################################################

### for-Schleife
### Eine Variable nimmt bei den Wiederholungen die Werte eines
### vorgegebenen Vektors an. Dann wird abgebrochen.


vec <- c("Eins","Zwei","Drei")
# Syntax for - ( Variable in Vektor ) - Auszuführender Befehl

for (v in vec) print(v)

## Eine for Schleife:

for (i in 1:10) {
  print(i)
}

## bewirkt, dass folgender Code ausgeführt wird

i <- 1
print(i)
i <- 2 
print(i)
# .
# .
# .
i <- 10
print(i)

## d.h. nach dem Ausführen ist der Laufindex (hier i) mit dem letzten Wert belegt (hier 10):

if (exists("i")) rm("i") ## löscht i falls es existiert

for (i in 1:10) {
  print(i)
}

i


# Häufig: 
# Variable stellt einen numerischen Laufindex dar (for (i in 1:1000) ... )
N <- 10
x <- rep(NA,times=N)


for (i in 1:N) {
	x[i] <- sin(2*pi*i/N)
}

plot(x,type="l")



# Es können auch mehrere Schleifen ineinander geschrieben werden

X <- matrix(NA, nrow= 4, ncol= 3)
X

for (i in 1:4){ # i Laufindex der Zeilen
    for (j in 1:3) # j Laufindex der Spalten
    {
	X[i,j] <- i/j
    }
}
X


### Aufgabe E2 
#a Berechnen Sie das Matrixprodukt aus A und B
#
A <- matrix(1:12,4,3)
B <- matrix(1:9,3,3)
C <- A%*%B
#
# ohne die %*% Funktion. 
# Verwenden Sie eine oder mehrere Schleifen
# und etwa sum(A[1,]*B[,1]) für C[1,1]...


# b) Erstellen Sie nun aus obigen Schleifen
# eine Funktion "mat_mult", die bei Eingabe 
# zweier Matrizen das Produkt berechnet, aber vorher
# mit if testet, ob die Zeilenlänge der einen Matrix
# die Spaltenlänge der anderen Matrix ist.
# Testen Sie das mit den Matrizen von oben.


mat_mult(A,B)

ProdAB <- mat_mult(A,B)
any( ProdAB != A %*% B)




### Weitere Schleifen: repeat und while (anschauen)
?"while"
s<-1


while(s<10){
  #k <- s+1 
  print(s)
  s <- s+1
  
}
print(s)



## Beispiel "Optimierung" einer Funktion

## Idee: wir wollen numerisch das Maximum einer nach unten geöffneten Parabel finden.
##       Wir starten bei einem Wert, wo wir sicher wissen, dass die Parabel noch ansteigt,
##       bewegen und schrittweise nach rechts überprüfen dabei, ob der Wert der Parabel noch
##       ausreichend größer wird. Ist dies nicht der Fall, stoppen wir.

parabel <- function(x) -x^2
plot(parabel , -2,2)

## (Sehr naive) Implementierung:
xx <- -2          ## hier starten wir
yy <- parabel(xx) ## welchen Wert hat die Parabel hier?

dY <- 1           ## wir initialisieren einen Wert, der in dem while () -statement TRUE ergibt
while( dY > 0.0000001){
  xx <- xx + 0.01       ## neuen xx einen schritt weiter rechts
  yyNew <- parabel(xx)  ## welchen wert hat die Parabel an dem neuen xx?
  dY <- yyNew - yy      ## um wieviel ist der Wert der Parabel durch diesen Schritt angestiegen
                        ##    dieses dY wird vor dem nächsten Schleifendurchlauf benutzt um zu
                        ##    überprüfen ob die Bedingung dY >0.01 noch erfüllt ist.
  yy <- yyNew           ## mach das yy für den nächsten Schleifendurchlauf zum neuen Wert
}

xx  ##  noch nicht ganz beim Maximum!


?"repeat"


#################################################
### apply-Befehle
#################################################

#### lapply: Funktionsanwendung auf Listen und Vektoren
?lapply

# Aufgabe: Wende eine Funktion (hier: sum) auf alle Listenelemente an

Liste <- list(a=c(4,8,7),
  b=seq(0,100,5),
  c=c(TRUE,TRUE,FALSE,TRUE))

Liste

# In kurzen Listen möglich:

sum(Liste$a) # oder äquivalent: sum(Liste[[1]])
sum(Liste$b)
sum(Liste$c)

# Bei längeren Listen denkbar: Schleife

for(i in 1:3) print(sum(Liste[[i]]))

# Eleganter (und mit großen Listen schneller): lapply und sapply


#### sapply
?lapply
?sapply

?apply

llist <- lapply(Liste, sum)
llist
class(llist)


slist <- sapply(Liste, sum)
slist[2]
class(slist)

sapply(Liste,sum)    # Ergebnis wird möglichst in Vektor/Matrix umgewandelt (s steht für simplify)

# Beispiel mit selbst geschriebener Funktion
funnyFun <- function(x, m, std) sum(x)/( 2 * rnorm(1, mean=m, sd=std)) 

sapply(X=Liste, FUN=funnyFun, m=2, std=5)

# Auch weitere Argumente zur auszuführenden Funktion sind möglich
# (hier: type, main, lwd)
dev.off()
par(mfrow=c(length(Liste),1))
lapply(Liste, plot, main="Titel", type="l", lwd=2)




### apply: Funktionen auf einzelne Zeilen/Spalten von Matrizen anwenden 
### (auch arrays mit >2 dims!)


A <- matrix(rnorm(12),4)    # rnorm() generiert standardnormalverteilte ZV
A

apply(A,1,var)             # Varianz der einzelnen Zeilen
var(A[1,])
apply(A,1,sum)
A
apply(A,2,var)             # Varianz der einzelnen Spalten




### tapply: Funktionen auf Kategorien einer diskreten Variable Anwenden


# Kurzes Beispiel:
himmelfarbe <- factor(c("blau","weiß","weiß","blau","blau-weiß"))
himmelfarbe
temperatur <- c(23,16,22,30,26)

tapply(temperatur, himmelfarbe, mean)

tapply(temperatur, himmelfarbe, plot, pch=16, cex=2)



### Aufgabe E3
# Erstellen Sie eine 10x10 Matrix, 
# in der als ij-tes Element (i*j)^(1/j) steht.
# Benutzen Sie hierfür zwei 
# ineinandergeschriebene Schleifen.
# Berechnen Sie dann die Spalten- und Zeilensummen 
# dieser Matrix mit einem ...apply-Befehl.





#################################################
### Paket parallel 
### *apply-Funktionen auf mehreren Rechnerkernen gleichzeitig
#################################################

# Hintergrund: Für rechenaufwendige Aufgaben ist es hilfreich, 
# Berechnungen auf verschiedene Prozessoren verteilen zu können
# (=parallel computing). Die Nutzung mehrerer Rechnerkerne auf
# EINEM PC schöpft die Ressourcen moderner Multicore-Prozessoren aus.
# Daneben kann man auch VERSCHIEDENE PC's oder Prozessoren von
# Hochleistungsclustern miteinander kommunizieren lassen. 
# (HPC-Cluster der Uni Regensburg: Athene, siehe RZ-Homepage)
# Ersteres ist mit R sehr einfach möglich (Paket parallel), 
# letzteres ist komplizierter.
# (z.B. Paket Rmpi für die Nutzung des "Message Processing Interface" unter R)

library(parallel)
?clusterApply

# Cluster erstellen, öffnet R-Prozesse (Knoten) im Hintergrund
cl <- makeCluster(2)
cl

# Analog zu lapply:
lapply(Liste,sum)
parLapply(cl, Liste, sum)

# Knoten muss wieder geschlossen werden.
stopCluster(cl)

# Vorsicht: In den Knoten sind Objekte des Workspace nicht verfügbar!
# Erzeuge a (nur auf localem R Process)
a <- 10

# ... lapply findets
lapply(1:2,function(x) a+x)

# ... parLapply nicht!
cl <- makeCluster(2)
parLapply(cl,1:2,function(x) a+x)

# -> benötigte Objekte müssen zuerst exportiert werden!
clusterExport(cl, "a")
parLapply(cl,1:2,function(x) a+x)

# Einen Befehl überall ausführen (automatisch Ergebnis anzeigen)
clusterEvalQ(cl,library(AER))
clusterEvalQ(cl, x <- rnorm(10))
clusterEvalQ(cl, ls())

stopCluster(cl)


#################################################
# Automatische Benennung von Objekten
#################################################
install.packages("AER")
library(AER)
data(CPS1985)
attach(CPS1985)
CPS1985
summary(CPS1985)


### assign() und get(): Automatisches benennen und verwenden von Objekten

?assign
x <- 1
assign("x",2)
x



char <- "ab"
assign(char , 20)
ab

get("x")
get(char)



## Beispiel: je mach Wert der Variable x soll entweder mean() oder sum() ausgeführt werden
myFunc1 <- mean
myFunc2 <- sum

x <- 1
get(paste("myFunc",x,sep=""))(1:10)
x <- 2
get(paste("myFunc",x,sep=""))(1:10)



names(CPS1985)
for(variable in names(CPS1985))
{
    assign(paste("summary_", variable, sep = ""), summary(get(variable)))
    print(paste("summary_", variable, sep = ""))
}


summary(wage)
wageNamen <- "wage" # Speicher Variable mit dem Namen der anderen Variablen
summary(get(wageNamen))
get(wageNamen)

summary_wage

# Was geschieht in der obigen schleife?
# Schleife läuft über folgende Variablen
names(CPS1985)

# verwende einen Durchlauf der Schleife, um diese zu prüfen: "wage"
variable <- names(CPS1985)[1]
variable
get(variable)
summary(get(variable))
paste("summary_", variable, sep = "")
assign(paste("summary_", variable, sep = ""), summary(get(variable)))

summary_wage
get("summary_wage")

summary(CPS1985)


#### Aufgabe E4
## a) Ersetzen Sie die Schleife über variable mit einem apply-Befehl

## b) Wie lässt sich für das vorliegende Problem die Verwendung von
##    assign und get vermeiden?




### Längeres Beispiel (in Ruhe anschauen...)
# Verschiedene Hypothesen zur Bundestagswahl

## Einlesen Datensatz Bundestagswahlergebnis
# http://www.bundeswahlleiter.de/de/bundestagswahlen/BTW_BUND_09/veroeffentlichungen/

setwd("G:/PMR")

BTW <- read.csv("data/bund09.csv", header=TRUE) # Einlesen: Komma getrennt, mit Spaltenüberschrift
summary(BTW)   # Kurzer Überblick
head(BTW)

# Wähle nur die Wahlkreisergebnisse (1 - 16 in Land-Variable)
BTW <- subset(BTW, Land!=99, select=c(2,3,4,7,12,16,20,24,28,32,36,40))  
names(BTW) <- c("Kreis","Land","Wahlberechtigte","Wähler","SPD","CDU","FDP","LINKE","GRÜNE","CSU","NPD","Piraten")
str(BTW)   # Zeigt die Struktur eines R-Objekts an

## Bundesland als 'factor' hinzufügen, Kreis als 'character'

BuLand <- c("Schleswig-Holstein","Hamburg","Niedersachsen","Bremen","NRW","Hessen","Rheinland-Pfalz","BW","Bayern","Saarland",
    "Berlin","Brandenburg","Mecklenburg-Vopo","Sachsen","Sachsen-Anhalt","Thüringen")

BTW$Land <- BuLand[BTW$Land]
class(BTW$Land)
summary(BTW$Land)                   # Vorsicht: BTW$Land ist character, kann nicht geplottet werden
plot(BTW$Land)

BTW$Land <- factor(BTW$Land)        # Umwandlung in factor, empfohlen für KATEGORIALE Daten!
plot(BTW$Land)                      # Auswirkung auf plot (s. letzte Stunde)...
summary(BTW$Land)                   # ...und auf summary-Befehl

BTW$Kreis <- as.character(BTW$Kreis)

attach(BTW)

## Graphische Analyse

# Wie groß sind die Wahlkreise
hist(Wahlberechtigte,probability=TRUE,breaks=30)
lines(density(Wahlberechtigte),col="darkred")

hist(SPD/Wähler,breaks=30)           # Wie sehr variiert der SPD- Anteil in den Kreisen

boxplot(I(LINKE/Wähler)~Land)   # Unterscheidet sich der Anteil der LINKEn unter den Ländern?

plot((BTW$CDU+BTW$CSU)/BTW$Wähler,BTW$Wahlberechtigte)  # Gibt es einen Zusammenhang zwischen CDU/CSU-Anteil und Wahlkreisgröße?

par(mfrow=c(4,4))                  # Gibt es einen solchen Zusammenhang in den Ländern?
for (l in BuLand)
{
    plot((CDU[Land==l]+CSU[Land==l])/Wähler[Land==l],Wahlberechtigte[Land==l],xlab="Stimmenanteil CDU/CSU", ylab="Wahlkreisgröße",main=l)
}

## Berechne Landesergebnisse

Land.Wähler <- tapply(Wähler,Land,sum)
Land.Berechtigte <- tapply(Wahlberechtigte,Land,sum)

Land.Stimmen <- matrix(NA,16,8)
colnames(Land.Stimmen) <- names(BTW)[-(1:4)]
rownames(Land.Stimmen) <- names(Land.Wähler)
Land.Anteile <- Land.Stimmen

Land.Stimmen <- as.data.frame(Land.Stimmen)
Land.Anteile <- as.data.frame(Land.Anteile)

for (i in 1:8)
{
	Land.Stimmen[,i] <- tapply(BTW[,i+4],Land,sum)
	Land.Anteile[,i] <- Land.Stimmen[,i]/Land.Berechtigte
}

## Berechne Spalte "Sonstige"

Sonstige <- Land.Wähler - apply(Land.Stimmen,1,sum)
Land.Stimmen$Sonstige <- Sonstige

Sonstige <- Sonstige/Land.Berechtigte
Land.Anteile$Sonstige <- Sonstige

## Berechne Spalte "Nichtwähler"

Nichtwähler <- Land.Berechtigte - Land.Wähler
Land.Stimmen$Nichtwähler <- Nichtwähler

Nichtwähler <- 1 - apply(Land.Anteile,1,sum)
Land.Anteile$Nichtwähler <- Nichtwähler

Land.Anteile


## Bundesergebnis

D.Berechtigte <- sum(Land.Berechtigte)
D.Wähler <- sum(Land.Wähler)

D.Stimmen <- apply(Land.Stimmen,2,sum)

(D.Anteile.1 <- D.Stimmen/D.Wähler)
(D.Anteile.2 <- D.Stimmen/D.Berechtigte)

## Grafik
dev.off()
pie(D.Anteile.1[c(2,6,3,1,4,5,7,8,9,10)],
	col=c("black","blue","yellow","red","darkred","green","brown","orange","pink","white"),
	main="Ergebnis mit Nichtwählern")

## Gesamtergebnis (über 5%)

pie(D.Anteile.2[c(2,6,3,1,4,5)],col=c("black","blue","yellow","red",
	"darkred","green"),main="Bundestagswahlergebnis")



#### Aufgabe E5
## Sie sollen grafisch die Hypothese der Einkommenskonvergenz zwischen den deutschen Kreisen und kreisfreien
## Städten untersuchen. Besorgen Sie von
## https://www.regionalstatistik.de
## das verfügbare Einkommen 2000 und möglichst aktuell. 
## Sie können die Tabellen als Text kopieren, in eine Datei speichern
## und in R importieren.
## Erstellen Sie aussagekräftige Grafiken, etwa die Veränderung des Einkommens auf der y-Achse vs. das Anfangsniveau
## auf der x-Achse.
## Führen Sie die Untersuchung auch auf Teilpopulationen durch.

## enthält einkommen von 2000 und 2014
gdp_nuts3 <- read.csv(file = "http://pc1011601230.ur.de/gdp_kreise.csv" ,
                      header = TRUE,
                      sep =";" , dec = "." ,
                      encoding = "latin1",
                      na.strings = "-" , as.is = TRUE)

View(gdp_nuts3)
names(gdp_nuts3)[c(1,2,6)] <- c("year", "kkz" , "gdppc")

View(gdp_nuts3)
gdp_nuts3$kkz <- as.numeric(gdp_nuts3$kkz)

selection <- subset( gdp_nuts3 ,as.numeric(kkz) > 1000 ,
                    select = c("year", "kkz" , "gdppc"))
selection$year <- as.numeric(selection$year)
selection$gdppc <- as.numeric(selection$gdppc) 
selection <- na.omit(selection)

View(selection)

gdp_nuts3_wd <- reshape( selection , direction = "wide" ,
                        v.names = "gdppc" , idvar= "kkz" , timevar = "year")

gdp_nuts3_wd$growth <- (gdp_nuts3_wd$gdppc.2014 / gdp_nuts3_wd$gdppc.2000)^(1/14)

with(gdp_nuts3_wd ,  plot( growth ~ gdppc.2000))
summary(est <- lm(growth ~ gdppc.2000 , data = gdp_nuts3_wd))
abline(est$coef)

## definiere neue Variable für neue Bundesländer:
gdp_nuts3_wd$ostwest <- as.numeric(gdp_nuts3_wd$kkz >= 11000)

with(gdp_nuts3_wd ,  plot( growth ~ gdppc.2000, col = as.factor(ostwest), pch = 16))
summary(est <- lm(growth ~ (gdppc.2000+ ostwest)^2 , data = gdp_nuts3_wd))
abline(est$coef[1:2] , col = "black" , lwd = 2)
abline(est$coef[1:2] + est$coef[3:4] , col = "red" ,lwd = 2)
legend( x = "topright" , legend = c("Ost" , "West") , pch = 16 , col = c("red" , "black"))