# Übungsblatt 3 - Evolutionary Algorithms (Paul Helstab & Yannik Bretschneider)

## Loading Data

```{r setwd}
# Get the current working directory
current_directory <- getwd()

# Define the folder name to check
folder_name <- "w3"

# Construct the full path to the folder
folder_path <- file.path(current_directory, folder_name)

# Check if the folder exists (means most likely in the repo root folder)
if (dir.exists(folder_path)) {
    setwd(folder_path)
    cat("Changed working directory to:", getwd(), "\n")
} else {
    cat("Folder 'w3' does not exist in the current working directory.
Staying where we are\n")
}
```


```{r read_data}
file_path <- "./data1.txt"

print(paste0("read data from '", file_path, "'"))
data <- read.table(
    file_path,
    header = FALSE,
    sep = " ",
    # colClasses = c("numeric", "numeric", "integer", "integer")
)
colnames(data) <- c("x1", "x2", "x_bias", "label")


evaluation_set_idx <- c(90:109)
```

```r
# possible, but unnecessary here:
# grab first 10 as evaluation set
evaluation_set <- data[evaluation_set_idx, ]

# grab the rest as training set
training_set <- data[-evaluation_set_idx, ]
```

## Aufgabe 1: Perzeptron

```{r peceptron-impl}
# randomized order
random_order <- sample(nrow(data))
vecs <- data[random_order, 1:3]
classes <- data[random_order, 4]

weights <- c(0, 0, 0)

# matrix for logging weights
wlog <- matrix(weights, ncol = 3)

infer <- function(x, wghts) {
    return(sign(x %*% wghts))
}

perceptron_infer <- function(x) {
    return(infer(x, weights))
}

while (TRUE) {
    mispredictions <- 0
    for (i in seq_len(nrow(vecs))) {
        vec <- as.numeric(vecs[i, ])
        class <- classes[i]

        # print(vec %*% weights)

        error <- as.numeric(class - perceptron_infer(vec))
        if (error != 0) {
            mispredictions <- mispredictions + 1

            print("adjusting")
            weights <- weights + error * vec
            wlog <- rbind(wlog, weights)
        }
    }


    if (mispredictions == 0) {
        print("No more mispredictions")
        break
    } else {
        print(paste0("Misclassifications: ", mispredictions))
    }
}

print("NOTE: for more interesting plots, run this chunk multiple times until multiple \"adjusting\" can be seen in the output")
```

```{r apply-to-data}
perceptron_reinfer <- t(sapply(seq_len(nrow(data)), function(i) {
    return(c(as.numeric(data[i, 1:3]), perceptron_infer((as.numeric(data[i, 1:3])))))
}))

# perceptron_reinfer
```


### Plots

The plots show the datapoints on a 2D plot, with the colors showing the classes (`green := 1` and `red := -1`)  
The separation plane defined by the normal of the weight vector is shown in increasing opacities, until the final 
weights are reached.

```{r plot-weights}
perceptron_plot <- function() {
    plot(data$x1, data$x2,
        col = ifelse(data$label == 1, "green", "red"),
        pch = 16, xlab = "x1", ylab = "x2", xlim = c(-5, 5), ylim = c(-5, 5)
    )

    legend("topright",
        legend = c("Class 1", "Class -1"),
        col = c("green", "red"), pch = 16
    )

    # Plot the separation planes
    num_steps <- nrow(wlog)
    alpha_step <- 1 / num_steps

    for (i in 1:num_steps) {
        a <- wlog[i, 1]
        b <- wlog[i, 2]
        c <- wlog[i, 3]

        # Calculate the line coordinates
        x <- seq(-5, 5, length.out = 100)
        y <- -(a * x + c) / b

        # Determine the line color based on the iteration
        alpha <- i * alpha_step
        line_color <- rgb(0, 0, 1, alpha)

        # Plot the line
        lines(x, y, col = line_color)
    }

    # Plot the final separation plane
    a <- weights[1]
    b <- weights[2]
    c <- weights[3]

    x <- seq(-5, 5, length.out = 100)
    y <- -(a * x + c) / b

    lines(x, y, col = "blue", lwd = 2)
}

perceptron_plot()
```

```{r plot-weights-pdf}
pdf("perceptron-plot.pdf", width = 8, height = 6)
perceptron_plot()
dev.off()
```



## Aufgabe 2: Perzeptron Evolutionary Algorithm (EA)


```{r peceptron-ea-impl}
# randomized order
ea_random_order <- sample(nrow(data))
ea_vecs <- data[ea_random_order, 1:3]
ea_classes <- data[ea_random_order, 4]

NUM_INDIVIDUALS <- 3
MUTATE_ALL_GENES <- FALSE
MUTATE_STEP <- function() {
    return(0.5)
}

# MUTATE_STEP <- function() {
#     return(runif(1, min = -1, max = 1)[1])
# }

weights <- matrix(0, ncol = 3, nrow = NUM_INDIVIDUALS)

# matrix for logging weights
ea_wlog <- matrix(0, ncol = 3)

fittest <- c(0, 0, 0)


mutate <- function(genes, step_size_fn, mutate_all) {
    if (mutate_all == TRUE) {
        for (i in seq_along(genes)) {
            step <- step_size_fn() * sample(c(1, -1), size = 1)[1] # mutation in step_size, -step_size or not at all
            genes[i] <- genes[i] + step
        }
    } else {
        idx_to_mutate <- sample(seq_along(genes), size = 1)[1]
        genes[idx_to_mutate] <- genes[idx_to_mutate] + step_size_fn() * sample(c(1, -1), size = 1)[1] # mutation in step_size or not AT all
    }

    return(genes)
}

infer <- function(x, wghts) {
    return(sign(x %*% wghts))
}

fitness <- function(wghts) {
    mispredictions <- 0


    for (i in seq_len(nrow(ea_vecs))) {
        vec <- as.numeric(ea_vecs[i, ])
        class <- ea_classes[i]

        # print(vec %*% weights)

        error <- as.numeric(class - infer(vec, wghts))
        if (error != 0) {
            mispredictions <- mispredictions + 1
        }
    }

    return(mispredictions)
}

while (TRUE) {
    i_fittest <- weights[1, ]
    i_fit <- fitness(weights[1, ])

    for (individual in 2:NUM_INDIVIDUALS) {
        this_fitness <- fitness(weights[individual, ])
        if (this_fitness < i_fit) {
            i_fittest <- weights[individual, ]
            i_fit <- this_fitness
        }
    }

    # add fittest to weights
    ea_wlog <- rbind(ea_wlog, i_fittest)

    if (i_fit == 0) {
        print("No more mispredictions")
        fittest <- i_fittest
        break
    } else {
        print(paste0("Misclassifications: ", i_fit))

        for (individual in 1:NUM_INDIVIDUALS) {
            weights[individual, ] <- mutate(i_fittest, step_size_fn = MUTATE_STEP, mutate_all = MUTATE_ALL_GENES)
        }
    }
}
```


```{r ea-apply-to-data}
ea_reinfer <- t(sapply(seq_len(nrow(data)), function(i) {
    return(c(as.numeric(data[i, 1:3]), infer((as.numeric(data[i, 1:3])), fittest)))
}))

# ea_reinfer
```

### Plots

The plots show the datapoints on a 2D plot, with the colors showing the classes (`green := 1` and `red := -1`)  
The separation plane defined by the normal of the weight vector is shown in increasing opacities, until the final 
weights are reached.

```{r ea-plot-weights}
ea_plot <- function() {
    plot(data$x1, data$x2,
        col = ifelse(data$label == 1, "green", "red"),
        pch = 16, xlab = "x1", ylab = "x2", xlim = c(-5, 5), ylim = c(-5, 5)
    )

    legend("topright",
        legend = c("Class 1", "Class -1"),
        col = c("green", "red"), pch = 16
    )

    # Plot the separation planes
    num_steps <- nrow(ea_wlog)
    alpha_step <- 1 / num_steps

    for (i in 1:num_steps) {
        a <- ea_wlog[i, 1]
        b <- ea_wlog[i, 2]
        c <- ea_wlog[i, 3]

        # Calculate the line coordinates
        x <- seq(-5, 5, length.out = 100)
        y <- -(a * x + c) / b

        # Determine the line color based on the iteration
        alpha <- i * alpha_step
        line_color <- rgb(0, 0, 1, alpha)

        # Plot the line
        lines(x, y, col = line_color)
    }

    # Plot the final separation plane
    a <- fittest[1]
    b <- fittest[2]
    c <- fittest[3]

    x <- seq(-5, 5, length.out = 100)
    y <- -(a * x + c) / b

    lines(x, y, col = "blue", lwd = 2)
}

ea_plot()
```


```{r plot-weights-pdf}
pdf("ea-plot_static-delta_one.pdf", width = 8, height = 6)
ea_plot()
dev.off()
```

### Evaluation

Gaussian (Mutate All) with a high number of individuals generally acheives a correct result first-try. In our very constrained individuum pool tests with 3 individuals, this also seems to be the best choice for parameters for the EA. Its efficacy seems to increase even more with a larger gene pool. 


## Aufgabe 3: Exploration vs Exploitation

```{r setup}
# generate a random binary vector of length 10
generate_individual <- function() {
    sample(c(0, 1), 10, replace = TRUE)
}

# 1-point crossover
one_point_crossover <- function(parent1, parent2) {
    point <- sample(1:9, 1)
    child1 <- c(parent1[1:point], parent2[(point + 1):10])
    child2 <- c(parent2[1:point], parent1[(point + 1):10])
    return(list(child1, child2))
}

# 2-point crossover
two_point_crossover <- function(parent1, parent2) {
    points <- sort(sample(1:9, 2))
    child1 <- c(parent1[1:points[1]], parent2[(points[1] + 1):points[2]], parent1[(points[2] + 1):10])
    child2 <- c(parent2[1:points[1]], parent1[(points[1] + 1):points[2]], parent2[(points[2] + 1):10])
    return(list(child1, child2))
}

# uniform crossover
uniform_crossover <- function(parent1, parent2) {
    mask <- sample(c(0, 1), 10, replace = TRUE)
    child1 <- ifelse(mask == 1, parent1, parent2)
    child2 <- ifelse(mask == 1, parent2, parent1)
    return(list(child1, child2))
}

# Function to convert binary vector to decimal
binary_to_decimal <- function(binary_vector) {
    sum(binary_vector * 2^(rev(seq_along(binary_vector) - 1)))
}
```


```{r initialize_parents}
# init parents
parent1 <- generate_individual()
parent2 <- generate_individual()

print(parent1)
print(parent2)

# Set num gens
num_generations <- 10000

# for tracking unique solutions
unique_solutions <- list()
unique_solutions[[binary_to_decimal(parent1)]] <- TRUE
unique_solutions[[binary_to_decimal(parent2)]] <- TRUE
```


```{r run_genetic_algorithm}
for (generation in 1:num_generations) {
    # Choose crossover method
    crossover_method <- sample(c("one_point", "two_point", "uniform"), 1)

    if (crossover_method == "one_point") {
        children <- one_point_crossover(parent1, parent2)
    } else if (crossover_method == "two_point") {
        children <- two_point_crossover(parent1, parent2)
    } else {
        children <- uniform_crossover(parent1, parent2)
    }



    # Replace
    parent1 <- children[[1]]
    parent2 <- children[[2]]

    # Track sols
    unique_solutions[[binary_to_decimal(parent1)]] <- TRUE
    unique_solutions[[binary_to_decimal(parent2)]] <- TRUE
}

# Calculate number of unique solutions observed
num_unique_solutions <- length(unique_solutions)

# Calculate number of possible solutions covered
num_possible_solutions <- 2^10
num_solutions_covered <- length(unique_solutions)

# Print results
cat("Number of unique solutions observed:", num_unique_solutions, "\n")
cat("Number of possible solutions covered:", num_solutions_covered, "\n")
cat("Percentage of 10-bit numbers covered:", round(num_solutions_covered / num_possible_solutions * 100, 2), "%\n")
```