#################################################
### 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"))