diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5cfdd35 --- /dev/null +++ b/.gitignore @@ -0,0 +1,33 @@ +torch/data/* +torch/runs/* +torch/checkpoints/* +*/runs/* +*/data/* +*/checkpoints/* +.directory + +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] + +# Distribution / packaging +bin/ +build/ +lib/ +lib64/ +var/ + +# Logging: +*.log +*.out + +models/ +weights/ +rl_models/ +rf_models/ +runs/ +ga_runs/ +process_control_model/params/ + +preprocessing/.ipynb_checkpoints/ +preprocessing/checkpoints/ diff --git a/README.md b/README.md new file mode 100644 index 0000000..308f8a9 --- /dev/null +++ b/README.md @@ -0,0 +1,39 @@ +# adaptive-spreading + +Code for the paper "Learning Controllers for Adaptive Spreading ofCarbon Fiber Tows". +Download data (preprocessed) and a pretrained process model here: +https://figshare.com/s/1a3e9b1ac16362b46cf9 + +## Preprocessing + +Scripts to preprocess the data are located in the subdirectory ./preprocessing. +The major preprocessing steps are: ++ Removing NANs ++ Applying Savitzky-Golay-Filter ++ Building the average of consecutive measurements ++ Shifting the tow to the middle of each measurement + +## Tow Prediction / Process Model + +### Feedforward Neural Networks + +### Random Forests + +## Process Control Model + +### Reward/Fitness Function + +The reward function consists of three parts: ++ Target height ++ Target width ++ Bar movement + +`cost = -(k_h * abs(target_height - mean_current_height) + k_w * abs(target_width - current_width) + k_m * total_bar_movement)` + +The three criteria can be scaled individually. + +### Algorithms + + ++ Genetic Algorithm (start_neuroevolution.py) + + Implements a rather simplistic GA with mutation. No crossover capabilities are provided. diff --git a/bar_setups.txt b/bar_setups.txt new file mode 100644 index 0000000..762142a --- /dev/null +++ b/bar_setups.txt @@ -0,0 +1,20 @@ +BAR_CONFIG = { + 'Setup0': [], + 'Setup1': [17.2, 16.1, 17.4, 16.1, 17.2], + 'Setup2': [17.2, 16.1, 22.4, 16.1, 17.2], + 'Setup3': [17.2, 11.2, 22.4, 11.1, 17.2], + 'Setup4': [22.2, 11.2, 22.4, 11.1, 22.4], + 'Setup5': [17.1, 11.2, 27.5, 11.1, 17.3], + 'Setup6': [17.2, 16.1, 17.2, 16.1, 17.3], + 'Setup7': [17.2, 16.1, 27.3, 16.1, 17.3], + 'Setup8': [17.2, 16.1, 32.2, 16.1, 17.3], + 'Setup9': [17.2, 16.1, 37.2, 16.1, 17.3], + 'Setup10': [22.1, 16.1, 37.2, 16.1, 22.1], + 'Setup11': [22.1, 16.1, 32.2, 16.1, 22.1], + 'Setup12': [27.3, 16.1, 17.6, 16.1, 27.3], + 'Setup13': [27.3, 11.0, 12.5, 11.0, 27.3], + 'Setup14': [27.3, 11.0, 31.2, 16.2, 22.2], + 'Setup15': [27.3, 19.5, 31.2, 7.8, 22.2], + 'Setup16': [27.3, 11.0, 12.5, 11.0, 27.3], + 'Setup17': [22.1, 16.1, 37.2, 16.1, 22.1], + 'Setup18': [17.2, 16.1, 32.2, 16.1, 17.3]} \ No newline at end of file diff --git a/plotter/__init__.py b/plotter/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/plotter/input_output_plotter.py b/plotter/input_output_plotter.py new file mode 100644 index 0000000..612d7c4 --- /dev/null +++ b/plotter/input_output_plotter.py @@ -0,0 +1,20 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt + + +def plot_profile(x): + plt.ylim((0., 2.)) + plt.plot(x, linewidth=1.4) + plt.axis('off') + plt.show() + + +if __name__ == '__main__': + input_file = '/media/julia/Data/datasets/lufpro/real/simple_preprocess3/Razieh05.09_1Tow_Setup8_6,0m_min_3.csv' + df = pd.read_csv(input_file, header=None) + x_data = df.loc[:, 5:804].to_numpy() + y_data = df.loc[:, 805:].to_numpy() + del df + plot_profile(x_data[30]) + plot_profile(y_data[30]) diff --git a/plotter/preprocessing_plotter.py b/plotter/preprocessing_plotter.py new file mode 100644 index 0000000..bf1e78e --- /dev/null +++ b/plotter/preprocessing_plotter.py @@ -0,0 +1,32 @@ +from pathlib import Path +import pandas as pd +import numpy as np +from preprocessing.simple_preprocessing import subproc +from process_model.model_comparison import plot_samples + +if __name__ == '__main__': + filename = Path('/media/julia/Data/datasets/lufpro/real/Razieh11.09_1Tow_Setup10_8,0m_min_3.csv') + #preprocessed_data = subproc(filename) + sensor_dim = 800 + params_dim = 17 + + data = pd.read_csv(filename) + data.columns = range(data.shape[1]) + num = data._get_numeric_data() + num[num < -100] = np.nan + xmean_val = np.nanmean(num.loc[:, params_dim:sensor_dim + params_dim - 1]) + ymean_val = np.nanmean(num.loc[:, sensor_dim + params_dim:]) + x = num.values + params_dim = x.shape[-1] - 2 * sensor_dim + + y = x[:, (sensor_dim + params_dim):] #- ymean_val + x = x[:, params_dim:(sensor_dim + params_dim)] #- xmean_val + + # x[x < xmean_val-0.5] = np.nan + # x[x > xmean_val+0.5] = np.nan + # y[y < ymean_val-0.5] = np.nan + # y[y > ymean_val+0.5] = np.nan + select_data = np.random.choice(x.shape[0], 8, replace=False) + # preprocessed_data[select_data, 5:805] + plot_samples([x[select_data], y[select_data]], legend=['raw', 'preprocessed'], + title='Pre-processing', plot_edges=False) diff --git a/plotter/profile_plotter.py b/plotter/profile_plotter.py new file mode 100644 index 0000000..777585c --- /dev/null +++ b/plotter/profile_plotter.py @@ -0,0 +1,31 @@ +import numpy as np +import torch +import matplotlib.pyplot as plt + + +def plot_samples(model_name, sensor_dim, data, meta=None): + num_samples = 10 + ixs = np.random.randint(len(data[0]), size=num_samples) + preds = data[0].view(-1, sensor_dim) + targets = data[1].view(-1, sensor_dim) + print(ixs) + plot_i = 1 + for i in ixs: + if meta: + plt.subplot(num_samples, 1, plot_i, + title=f'{meta[0][i]}, line {int(meta[1][i].item())}') + plt.subplots_adjust(hspace=.5) + else: + plt.subplot(num_samples, 1, plot_i) + plt.plot(targets[i], linewidth=0.4) + plt.plot(preds[i], linewidth=0.4) + plot_i += 1 + plt.legend(['target', 'prediction']) + plt.suptitle(model_name) + plt.show() + + +def plot_samples_np_wrapper(preds, targets, sensor_dim, name='Random forest', meta=None): + if meta is not None: + meta = meta[0], torch.tensor(meta[1]) + plot_samples(name, sensor_dim, (torch.tensor(preds), torch.tensor(targets)), meta=meta) diff --git a/plotter/rl_plotter.py b/plotter/rl_plotter.py new file mode 100644 index 0000000..216c39e --- /dev/null +++ b/plotter/rl_plotter.py @@ -0,0 +1,36 @@ +from _tkinter import TclError +import matplotlib.pyplot as plt + +from preprocessing.tape_detection import get_tape_edges +from process_control_model.rl_utils import create_target + + +def render(tape_in, tape_out, bar_positions, target_width, target_height, sensor_dim=800, setup_dim=5): + target = create_target(target_width, target_height, dim=1, sensor_dim=sensor_dim).squeeze() + fig, axs = plt.subplots(3) + fig.set_size_inches(16, 14) + plt.show(block=False) + raising_edges, falling_edges = get_tape_edges(tape_out) + for i, (start, row, e1, e2, action) in enumerate(zip(tape_in, tape_out, raising_edges, falling_edges, + bar_positions)): + try: + axs[0].clear() + axs[0].plot(start, '-b', label='start values') + axs[0].legend() + axs[0].set_ylim([-0.3, 0.6]) + axs[1].clear() + axs[1].plot(target, color='red', label='target values') + axs[1].plot(row, '-bo', markerfacecolor='r', + markevery=[e1, e2], label='end values') + axs[1].legend() + axs[1].set_ylim([-0.3, 0.6]) + axs[2].clear() + axs[2].plot(action, '-ko', markerfacecolor='r', markevery=[i for i in range(setup_dim)], + label='bar position') + axs[2].legend() + axs[2].set_ylim([10, 40]) + fig.canvas.draw() + plt.pause(0.001) + except (KeyboardInterrupt, TclError): + break + plt.close(fig) diff --git a/preprocessing/__init__.py b/preprocessing/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/preprocessing/datasetloader.py b/preprocessing/datasetloader.py new file mode 100644 index 0000000..0f3b565 --- /dev/null +++ b/preprocessing/datasetloader.py @@ -0,0 +1,62 @@ +import logging +from pathlib import Path +import torch +from torch.utils.data.dataloader import DataLoader + +from preprocessing.datasets import DefaultDataset, ExtendedMetaDataset +from preprocessing.split_datasets import get_test_set_respecting_files,\ + split_data_respecting_files_including_meta, get_test_set_including_meta +import utils.paths as dirs + + +LOGGER = logging.getLogger(__name__) + + +def load_data_and_meta(device, batch_size, *args): + train_data, test_data, val_data = split_data_respecting_files_including_meta(*args) + + def get_dataloader(data_x, data_y, data_filenames, data_index, shuffle=False, batchsize=1): + return DataLoader(ExtendedMetaDataset(torch.from_numpy(data_x).float().to(device), + torch.from_numpy(data_y).float().to(device), + data_filenames, data_index), + batch_size=batchsize, shuffle=shuffle) + train_loader = get_dataloader(*train_data, shuffle=True, batchsize=batch_size) + validation_loader = get_dataloader(*val_data) + test_loader = get_dataloader(*test_data) + return train_loader, validation_loader, test_loader + + +def load_test_data_and_meta(device, datapath, sensor_dim, setup_dim, random_seed=42): + test_x, test_y, filenames, indices = get_test_set_including_meta(datapath, sensor_dim, + setup_dim, random_seed) + return DataLoader(ExtendedMetaDataset(torch.from_numpy(test_x).float().to(device), + torch.from_numpy(test_y).float().to(device), + filenames, indices)) + + +def load_data(device, split_function, batch_size, *args): + train_x, test_x, val_x, train_y, test_y, val_y = split_function(*args) + LOGGER.debug(f'size training set: {train_x.shape}, test set: {test_y.shape}') + train_loader = DataLoader(DefaultDataset(torch.from_numpy(train_x).float().to(device), + torch.from_numpy(train_y).float().to(device)), + batch_size=batch_size, shuffle=True) + validation_loader = DataLoader(DefaultDataset(torch.from_numpy(val_x).float().to(device), + torch.from_numpy(val_y).float().to(device))) + test_loader = DataLoader(DefaultDataset(torch.from_numpy(test_x).float().to(device), + torch.from_numpy(test_y).float().to(device))) + return train_loader, validation_loader, test_loader + + +def load_test_data(device, datapath, sensor_dim, setup_dim, random_seed=42): + test_x, test_y = get_test_set_respecting_files( + datapath, sensor_dim, setup_dim, random_seed) + return DataLoader(DefaultDataset(torch.from_numpy(test_x).float().to(device), + torch.from_numpy(test_y).float().to(device))) + + +if __name__ == '__main__': + DEVICE = torch.device("cpu") + SPLIT_ARGS = [Path(dirs.DATA), 800, 5, False, 1111] + _, _, TEST_LOADER = load_data_and_meta(DEVICE, 20, *SPLIT_ARGS) + _, _, FILENAME, INDEX = TEST_LOADER.__iter__().__next__() + print(FILENAME[0], INDEX[0].item()) diff --git a/preprocessing/datasets.py b/preprocessing/datasets.py new file mode 100644 index 0000000..e7961a6 --- /dev/null +++ b/preprocessing/datasets.py @@ -0,0 +1,32 @@ +from torch.utils.data import Dataset + + +class DefaultDataset(Dataset): + def __init__(self, input_profiles, target_profiles): + super(DefaultDataset, self).__init__() + self.input_profiles = input_profiles + self.target_profiles = target_profiles + self.datasetsize = len(self.input_profiles) + + def __getitem__(self, index): + return self.input_profiles[index], self.target_profiles[index] + + def __len__(self): + return self.datasetsize + + +class ExtendedMetaDataset(Dataset): + def __init__(self, input_profiles, target_profiles, filenames, lines): + super(ExtendedMetaDataset, self).__init__() + self.input_profiles = input_profiles + self.target_profiles = target_profiles + self.filenames = filenames + self.lines = lines + self.dataset_size = len(self.input_profiles) + + def __getitem__(self, index): + return self.input_profiles[index], self.target_profiles[index], \ + self.filenames[index], self.lines[index] + + def __len__(self): + return self.dataset_size diff --git a/preprocessing/preprocessing_helper.py b/preprocessing/preprocessing_helper.py new file mode 100644 index 0000000..9f1be5a --- /dev/null +++ b/preprocessing/preprocessing_helper.py @@ -0,0 +1,125 @@ +import math +import re +from datetime import datetime +import numpy as np +import torch + + +MAX_ACTION = 60. + +ACTION_RANGES = torch.tensor([[17., 27.], [11., 19.], [12., 37.], [11., 16.], [17., 27.]]) + +# bar positions according to setup files (see bar_setups.txt oder Setup.pdf) + +BAR_CONFIG = { + 'Setup1': [17.2, 16.1, 17.4, 16.1, 17.2], + 'Setup2': [17.2, 16.1, 22.4, 16.1, 17.2], + 'Setup3': [17.2, 11.2, 22.4, 11.1, 17.2], + 'Setup4': [22.2, 11.2, 22.4, 11.1, 22.4], + 'Setup5': [17.1, 11.2, 27.5, 11.1, 17.3], + 'Setup6': [17.2, 16.1, 17.2, 16.1, 17.3], + 'Setup7': [17.2, 16.1, 27.3, 16.1, 17.3], + 'Setup8': [17.2, 16.1, 32.2, 16.1, 17.3], + 'Setup9': [17.2, 16.1, 37.2, 16.1, 17.3], + 'Setup10': [22.1, 16.1, 37.2, 16.1, 22.1], + 'Setup11': [22.1, 16.1, 32.2, 16.1, 22.1], + 'Setup12': [27.3, 16.1, 17.6, 16.1, 27.3], + 'Setup13': [27.3, 11.0, 12.5, 11.0, 27.3], + 'Setup14': [27.3, 11.0, 31.2, 16.2, 22.2], + 'Setup15': [27.3, 19.5, 31.2, 7.8, 22.2], + 'Setup16': [27.3, 11.0, 12.5, 11.0, 27.3], + 'Setup17': [22.1, 16.1, 37.2, 16.1, 22.1], + 'Setup18': [17.2, 16.1, 32.2, 16.1, 17.3]} + + +FILENAME_DELIMITER = '_' +FILENAME_SETUP_POS = 0 +FILENAME_VELO_POS = 1 + +DATE_FORMAT = '%m/%d/%Y %I:%M:%S %p' + + +BAD_FILES = ['Setup10_2,7,0m_min_1.csv', + 'Setup10_2,7,0m_min_2.csv', + 'Setup10_2,7,0m_min_3.csv', + 'Setup7_2,7m_min_1.csv'] + +BAR_GAP = 8.0 +BASE = 5.0 # measure tape for 5mm +FREQUENCY = 500 + + +def decode_filename(name): + comp = name.split(FILENAME_DELIMITER) + setup = comp[FILENAME_SETUP_POS] + velo = comp[FILENAME_VELO_POS] + velo = float(velo[:-1].replace(',', '.')) + return setup, velo + + +def calc_offset_index(offset, time_entries): + end = time_entries[0] + offset + for i, time in enumerate(time_entries): + if time >= end: + return i + return len(time_entries) - 1 + + +def calc_offset_time(setup, velocity): + bar_pos = BAR_CONFIG[setup] + # meter/minute to cm/second + velocity *= 100.0 + velocity /= 60.0 + c = 0 + for i in range(len(bar_pos) - 1): + # pythagoras + a = BAR_GAP + b = abs(bar_pos[i] - bar_pos[i + 1]) + c += math.sqrt(a ** 2 + b ** 2) + return c / velocity + + +def normalize_time(timestamp): + time_normalized = [] + first_day = timestamp[0][0] + first_day = datetime.strptime(first_day, DATE_FORMAT) + for time in timestamp: + current_day = time[0] + current_day = datetime.strptime(current_day, DATE_FORMAT) + difference = (current_day - first_day).total_seconds() + time_normalized.append(difference) + return time_normalized + + +def add_ms(time): + time_precise = [] + frequencies = check_frequencies(time) + base_ms = 1 / FREQUENCY + current_entry = 0 + for i in range(frequencies[0]): + precise = time[current_entry] + \ + ((i + (FREQUENCY - frequencies[0])) * base_ms) + time_precise.append(precise) + current_entry += 1 + for i in frequencies[1:]: + for j in range(i): + precise = time[current_entry] + (j * base_ms) + time_precise.append(precise) + current_entry += 1 + return time_precise + + +def check_frequencies(time): + frequency_count = np.bincount(time) + return frequency_count + + +def calc_windowsize(velocity): + velo_precise = (velocity / 60) * 1000 + p = BASE / velo_precise + window_size = int(p * FREQUENCY) + return window_size + + +def get_bar_positions(setup): + return np.array(BAR_CONFIG[setup]) diff --git a/preprocessing/simple_preprocessing.py b/preprocessing/simple_preprocessing.py new file mode 100644 index 0000000..a199c6f --- /dev/null +++ b/preprocessing/simple_preprocessing.py @@ -0,0 +1,229 @@ +import multiprocessing +from pathlib import Path +import numpy as np +import pandas as pd +from scipy import signal, interpolate +from preprocessing.preprocessing_helper import decode_filename, BAD_FILES, get_bar_positions, \ + calc_offset_time, normalize_time, add_ms, calc_offset_index, calc_windowsize +import utils.paths as dirs + + +def calculate_bar(x, support_points=20): + ixs = list(range(support_points)) + \ + list(range(x.shape[-1] - support_points, x.shape[-1])) + return interpolate.interp1d(ixs, x[ixs]) + + +def get_indices_over_threshold(x, threshold=None, merge_distance=None): + if np.isnan(x).all(): + return np.array([np.nan, np.nan]) + if not threshold: + threshold = np.nanmean(x) + if not merge_distance: + merge_distance = int(np.nanstd(x) * 165) + if merge_distance > 200: + merge_distance = 29 + # print(threshold, merge_distance) + # index of all values greater threshold + idx = x > threshold + if not np.any(idx): + return np.array([0, 0]) + # array of start and end points of areas containing values greater threshold + idx = np.argwhere(idx != np.roll(idx, -1)) + # merge near areas + even_idx = idx[0::2][np.where( + np.abs(idx[0::2] - np.roll(idx, 1)[0::2]) > merge_distance)] + odd_idx = idx[1::2][np.where( + np.abs(np.roll(idx, -1)[1::2] - idx[1::2]) > merge_distance)] + idx = np.array(list(zip(even_idx, odd_idx))).flatten() + # remove very small areas (peaks) + even_idx = idx[1::2][np.where(idx[1::2] - np.roll(idx, 1)[1::2] > 50)] + odd_idx = idx[0::2][np.where(np.roll(idx, -1)[0::2] - idx[0::2] > 50)] + if not len(even_idx) or not len(odd_idx): + return np.array([0, 0]) + # return start and end point of the tape hopefully + return np.array([odd_idx[0], even_idx[-1]]) + + +def preprocess_series(data): + sensor_dim = 800 + params_dim = 17 + data.columns = range(data.shape[1]) + timestamp = add_ms(normalize_time(data.loc[:, 0].to_numpy()[:, np.newaxis])) + num = data._get_numeric_data() + num[num < -100] = np.nan + xmean_val = np.nanmean(num.loc[:, params_dim:sensor_dim + params_dim - 1]) + ymean_val = np.nanmean(num.loc[:, sensor_dim + params_dim:]) + x = num.values + params_dim = x.shape[-1] - 2 * sensor_dim + print(x.shape) + + y = x[:, (sensor_dim + params_dim):] - ymean_val + x = x[:, params_dim:(sensor_dim + params_dim)] - xmean_val + + x[x < -0.3] = np.nan + x[x > 0.3] = np.nan + y[y < -0.2] = np.nan + y[y > 0.2] = np.nan + + x = df_fillna(pd.DataFrame(x)) + y = df_fillna(pd.DataFrame(y)) + x, y = smooth_profile(x, y) + + x_sav = signal.savgol_filter(x, 15, 3) + y_sav = signal.savgol_filter(y, 15, 3) + return timestamp, x_sav, y_sav + + +def replace_peaks(profile): + peaks = signal.find_peaks(profile, threshold=0.01)[0] + widths = signal.peak_widths(profile, peaks, rel_height=0.4) + widths = zip(widths[2].astype(int), widths[3].astype(int) + 1) + widths = [index for peak_start, peak_end in widths + for index in range(peak_start, peak_end) if peak_end - peak_start < 150] + profile[widths] = np.nan + return profile + + +def smooth_profile(x, y): + for row in x: + replace_peaks(row) + for row in y: + replace_peaks(row) + x = df_fillna(pd.DataFrame(x)) + y = df_fillna(pd.DataFrame(y)) + return x, y + + +def df_fillna(data): + data = data.interpolate(limit_direction='both', method='linear', axis=1) + return data.values + + +def correct_tape_idx(tape_idx, normalized, width_bounds, height_thresholds=(0.045, 0.075), + min_distances=(29, 9)): + tape_width = tape_idx[:, -1] - tape_idx[:, 0] + tape_width_in_range = np.where(tape_width < width_bounds[0]) + if tape_width_in_range[0].size: + tape_idx[tape_width_in_range] = np.apply_along_axis(get_indices_over_threshold, axis=1, + arr=normalized[tape_width_in_range], + threshold=height_thresholds[0], + merge_distance=min_distances[0]) + tape_width = tape_idx[:, -1] - tape_idx[:, 0] + tape_width_in_range = np.where(tape_width > width_bounds[1]) + if tape_width_in_range[0].size: + tape_idx[tape_width_in_range] = np.apply_along_axis(get_indices_over_threshold, axis=1, + arr=normalized[tape_width_in_range], + threshold=height_thresholds[1], + merge_distance=min_distances[1]) + + +def subproc(f): + sensor_dim = 800 + x_values = np.arange(sensor_dim) + print(f.name) + setup, velocity = decode_filename(f.name) + # remove major flaws + time, x, y = preprocess_series(pd.read_csv(f)) + # rotate data so that the underlying bar is horizontal + x_bars = np.array([calculate_bar(sample)(x_values) for sample in x]) + y_bars = np.array([calculate_bar(sample)(x_values) for sample in y]) + normalized_x = x - x_bars + normalized_y = y - y_bars + print('normalized shapes', normalized_x.shape, normalized_y.shape) + # remove time offset of x and y data + offset_time = calc_offset_time(setup, velocity) + offset_index = calc_offset_index(offset_time, time) + normalized_x = normalized_x[:-offset_index] + normalized_y = normalized_y[offset_index:] + print('after offset', normalized_x.shape, normalized_y.shape) + + # average data + normalized_x, normalized_y = average_data(normalized_x, normalized_y, velocity) + normalized_x, normalized_y = remove_nans(normalized_x, normalized_y) + print('after averaging and removing nans', normalized_x.shape, normalized_y.shape) + # get start and end point of each tape measurement + x_tape_idx = np.apply_along_axis( + get_indices_over_threshold, axis=1, arr=normalized_x) + y_tape_idx = np.apply_along_axis( + get_indices_over_threshold, axis=1, arr=normalized_y) + # try to fix start and end points where tape width is out of bounds + correct_tape_idx(x_tape_idx, normalized_x, (180, 400)) + correct_tape_idx(y_tape_idx, normalized_y, (270, 450)) + + # center data + center_data(x_tape_idx, y_tape_idx, normalized_x, normalized_y) + bar_positions = np.repeat(np.expand_dims( + get_bar_positions(setup), axis=0), normalized_x.shape[0], axis=0) + data = np.concatenate((bar_positions, normalized_x, normalized_y), axis=1) + # write to file + # pd.DataFrame(data).to_csv(output_dir / f.name, index=False, + # header=False, float_format='%.3f') + return data + + +def remove_nans(normalized_x, normalized_y): + x_nans = ~np.isnan(normalized_x).all(axis=1) + normalized_x = normalized_x[x_nans] + normalized_y = normalized_y[x_nans] + y_nans = ~np.isnan(normalized_y).all(axis=1) + normalized_x = normalized_x[y_nans] + normalized_y = normalized_y[y_nans] + return normalized_x, normalized_y + + +def average_data(normalized_x, normalized_y, velocity): + win_size = calc_windowsize(velocity) + + def grouped_avg(a, n=2.): + result = np.cumsum(a, 0)[n - 1::n] / n + result[1:] = result[1:] - result[:-1] + return result + + normalized_x = grouped_avg(normalized_x, win_size) + normalized_y = grouped_avg(normalized_y, win_size) + return normalized_x, normalized_y + + +def center_data(x_tape_idx, y_tape_idx, x_data, y_data): + tape_width = x_tape_idx[:, -1] - x_tape_idx[:, 0] + centered_x_idx = np.empty_like(x_tape_idx) + centered_x_idx[:, 0] = x_data.shape[-1] // 2 - (tape_width // 2) + centered_x_idx[:, 1] = centered_x_idx[:, 0] + tape_width + tape_width = y_tape_idx[:, -1] - y_tape_idx[:, 0] + centered_y_idx = np.empty_like(y_tape_idx) + centered_y_idx[:, 0] = y_data.shape[-1] // 2 - (tape_width // 2) + centered_y_idx[:, 1] = centered_y_idx[:, 0] + tape_width + x_data[((centered_x_idx[:, 0][:, np.newaxis] - 1 < np.arange(x_data.shape[-1])) & + (centered_x_idx[:, -1][:, np.newaxis] + 1 > np.arange(x_data.shape[-1])))] = x_data[ + ((x_tape_idx[:, 0][:, np.newaxis] - 1 < np.arange(x_data.shape[-1])) & + (x_tape_idx[:, -1][:, np.newaxis] + 1 > np.arange(x_data.shape[-1])))] + + y_data[((centered_y_idx[:, 0][:, np.newaxis] - 1 < np.arange(y_data.shape[-1])) & + (centered_y_idx[:, -1][:, np.newaxis] + 1 > np.arange(y_data.shape[-1])))] = y_data[ + ((y_tape_idx[:, 0][:, np.newaxis] - 1 < np.arange(y_data.shape[-1])) & + (y_tape_idx[:, -1][:, np.newaxis] + 1 > np.arange(y_data.shape[-1])))] + + x_data[~((centered_x_idx[:, 0][:, np.newaxis] - 1 < np.arange(x_data.shape[-1])) & + (centered_x_idx[:, -1][:, np.newaxis] + 1 > np.arange(x_data.shape[-1])))] = 0 + y_data[~((centered_y_idx[:, 0][:, np.newaxis] - 1 < np.arange(y_data.shape[-1])) & + (centered_y_idx[:, -1][:, np.newaxis] + 1 > np.arange(y_data.shape[-1])))] = 0 + + +if __name__ == '__main__': + data_path = Path(dirs.DATA) + output_dir = data_path / 'preprocessed' + if not output_dir.exists(): + output_dir.mkdir() + files = [f for f in data_path.iterdir() if f.suffix == + '.csv' and f.name not in BAD_FILES] + print(len(files)) + files = [f for f in files if not (output_dir / f.name).exists()] + print(len(files)) + subproc(Path('/media/julia/Data/datasets/lufpro/real/Razieh27.08_1Tow_Setup5_6,0m_min_3.csv')) + # preprocess one file after another (single core usage) + # for f in files: + # subproc(f) + # or use multiprocessing + # pool = multiprocessing.Pool(6) + # pool.map(subproc, files) diff --git a/preprocessing/split_datasets.py b/preprocessing/split_datasets.py new file mode 100644 index 0000000..7eb5509 --- /dev/null +++ b/preprocessing/split_datasets.py @@ -0,0 +1,173 @@ +from multiprocessing import Pool +from pathlib import Path +import pandas as pd +from sklearn.model_selection import train_test_split +from sklearn.utils import shuffle +import numpy as np +from preprocessing.preprocessing_helper import MAX_ACTION +import utils.paths as dirs + + +NUM_PROCESSES = 6 + + +def scale_setup_values(data, setup_dim): + data[:, :setup_dim] /= MAX_ACTION + + +def shuffled_x_y_split(data, sensor_dim, setup_dim, profile_only): + if profile_only: + train_x, test_x, train_y, test_y = train_test_split( + data[:, setup_dim: sensor_dim + setup_dim], + data[:, sensor_dim + setup_dim:], test_size=0.2, random_state=42) + train_x, val_x, train_y, val_y = train_test_split( + train_x, train_y, test_size=0.2, random_state=42) + else: + scale_setup_values(data, setup_dim) + train_x, test_x, train_y, test_y = train_test_split( + data[:, : sensor_dim + setup_dim], data[:, sensor_dim + setup_dim:], + test_size=0.2, random_state=42) + train_x, val_x, train_y, val_y = train_test_split( + train_x, train_y, test_size=0.2, random_state=42) + return train_x, test_x, val_x, train_y, test_y, val_y + + +def x_y_split(data, sensor_dim, setup_dim, profile_only): + if profile_only: + data_x, data_y = data[:, setup_dim: sensor_dim + setup_dim],\ + data[:, sensor_dim + setup_dim:] + else: + scale_setup_values(data, setup_dim) + data_x, data_y = data[:, : sensor_dim + setup_dim], data[:, sensor_dim + setup_dim:] + return data_x, data_y + + +def get_meta_data(filenames, data, sensor_dim, setup_dim, profile_only): + index = [np.arange(len(df)) for df in data] + filenames = [f.name for i, f in enumerate(filenames) for _ in range(len(index[i]))] + index = np.concatenate(index) + data_x, data_y = x_y_split(np.concatenate(data), sensor_dim, setup_dim, profile_only) + return data_x, data_y, filenames, index + + +def ordered_split_respecting_files(training_files, test_files, sensor_dim, setup_dim, profile_only): + training_files, validation_files = train_test_split(training_files, test_size=0.2, + random_state=42) + training_data = mp_read_file_lists(training_files) + testing_data = mp_read_file_lists(test_files) + validation_data = mp_read_file_lists(validation_files) + train_x, train_y = x_y_split(np.concatenate(training_data), sensor_dim, setup_dim, profile_only) + test_x, test_y = x_y_split(np.concatenate(testing_data), sensor_dim, setup_dim, profile_only) + val_x, val_y = x_y_split(np.concatenate(validation_data), sensor_dim, setup_dim, profile_only) + return train_x, test_x, val_x, train_y, test_y, val_y + + +def ordered_split_and_meta(training_files, test_files, sensor_dim, setup_dim, profile_only): + training_files, validation_files = train_test_split(training_files, test_size=0.2, + random_state=42) + training_data = mp_read_file_lists(training_files) + testing_data = mp_read_file_lists(test_files) + validation_data = mp_read_file_lists(validation_files) + return get_meta_data(training_files, training_data, sensor_dim, setup_dim, profile_only),\ + get_meta_data(test_files, testing_data, sensor_dim, setup_dim, profile_only),\ + get_meta_data(validation_files, validation_data, sensor_dim, setup_dim, profile_only) + + +def mp_read_file_lists(files): + pool = Pool(NUM_PROCESSES) + data = pool.map(_mp_read_csv, files) + pool.close() + pool.join() + return data + + +def _mp_read_csv(filename): + return pd.read_csv(filename, header=None).to_numpy() + + +def read_csvs_from_list(filenames): + data = None + for filename in filenames: + temp_df = pd.read_csv(filename, header=None) + data = temp_df if data is None else data.append(temp_df) + return data + + +def load_one_file_and_split(filepath, sensor_dim, setup_dim, profile_only=False): + data = pd.read_csv(filepath, header=None).to_numpy() + return shuffled_x_y_split(data, sensor_dim, setup_dim, profile_only) + + +def shuffle_all_data_and_split(data_path, sensor_dim, setup_dim, profile_only=False): + files = [f for f in data_path.iterdir() if f.suffix == '.csv'] + data = read_csvs_from_list(files).to_numpy() + return shuffled_x_y_split(data, sensor_dim, setup_dim, profile_only) + + +def split_data_excluding_multiple_setups(data_path, sensor_dim, setup_dim, + profile_only=False, setups=(8, 18)): + substrings = [f'Setup{setup}' for setup in setups] + training_files = shuffle([f for f in data_path.iterdir() + if f.suffix == '.csv' + and not any(s in f.name for s in substrings)]) + test_files = [f for f in data_path.iterdir() if f.suffix == '.csv' + and any(s in f.name for s in substrings)] + return ordered_split_respecting_files(training_files, test_files, sensor_dim, + setup_dim, profile_only) + + +def split_data_excluding_one_setup(data_path, sensor_dim, setup_dim, profile_only=False, setup=8): + substring = f'Setup{setup}_' + training_files = shuffle([f for f in data_path.iterdir() + if f.suffix == '.csv' + and substring not in f.name]) + test_files = [f for f in data_path.iterdir() if f.suffix == '.csv' + and substring in f.name] + return ordered_split_respecting_files(training_files, test_files, sensor_dim, + setup_dim, profile_only) + + +def split_data_respecting_files_including_meta(data_path, sensor_dim, setup_dim, + profile_only=False, random_seed=42): + files = shuffle([f for f in data_path.iterdir() if f.suffix == '.csv'], + random_state=random_seed) + training_files, test_files = train_test_split(files, test_size=0.2, + random_state=random_seed) + return ordered_split_and_meta(training_files, test_files, sensor_dim, + setup_dim, profile_only) + + +def get_test_set_including_meta(data_path, sensor_dim, setup_dim, random_seed=42): + files = shuffle([f for f in data_path.iterdir() if f.suffix == '.csv'], + random_state=random_seed) + _, test_files = train_test_split(files, test_size=0.2, random_state=random_seed) + return get_meta_data(test_files, mp_read_file_lists(test_files), sensor_dim, + setup_dim, profile_only=False) + + +def split_data_respecting_files(data_path, sensor_dim, setup_dim, + profile_only=False, random_seed=42): + files = shuffle([f for f in data_path.iterdir() if f.suffix == '.csv'], + random_state=random_seed) + training_files, test_files = train_test_split(files, test_size=0.2, + random_state=random_seed) + return ordered_split_respecting_files(training_files, test_files, sensor_dim, + setup_dim, profile_only) + + +def get_test_set_respecting_files(data_path, sensor_dim, setup_dim, random_seed=42): + files = shuffle([f for f in data_path.iterdir() if f.suffix == '.csv'], + random_state=random_seed) + _, test_files = train_test_split(files, test_size=0.2, random_state=random_seed) + test_data = np.concatenate(mp_read_file_lists(test_files)) + test_x, test_y = x_y_split(test_data, sensor_dim, setup_dim, profile_only=False) + return test_x, test_y + + +if __name__ == '__main__': + FUNC = split_data_respecting_files + TRAIN_X, TEST_X, VAL_X, TRAIN_Y, TEST_Y, VAL_Y = FUNC( + Path(dirs.DATA), + 800, 5) + print(FUNC.__name__, TRAIN_X.shape, TRAIN_Y.shape, TEST_X.shape, + TEST_Y.shape, VAL_X.shape, VAL_Y.shape) diff --git a/preprocessing/subset_generation.py b/preprocessing/subset_generation.py new file mode 100644 index 0000000..3d21c42 --- /dev/null +++ b/preprocessing/subset_generation.py @@ -0,0 +1,66 @@ +from multiprocessing import Pool +from pathlib import Path +import numpy as np +import pandas as pd +import utils.paths as dirs + + +SOURCE_DIR = Path(dirs.DATA) +RL_START_FILE = Path('../process_control_model/start_values.csv') +NON_TOUCHING_DIR = Path(dirs.DATA) + +SETUP_DIM = 5 +SENSOR_DIM = 800 + +PROCESS_COUNT = 6 + + +def _get_random_start_values_from_file(filename, num_lines=2): + df = pd.read_csv(filename, header=None) + random_idx = [np.random.randint(low=0, high=df.shape[0]) for _ in range(num_lines)] + return df.loc[random_idx, SETUP_DIM:(SETUP_DIM + SENSOR_DIM - 1)].to_numpy() + + +def _get_random_lines_from_file(filename, num_lines=2): + df = pd.read_csv(filename, header=None) + random_idx = [np.random.randint(low=0, high=df.shape[0]) for _ in range(num_lines)] + return df.loc[random_idx, :].to_numpy() + + +def _mp_get_subset(func, lines_per_file=2): + files = [(file, lines_per_file) for file in SOURCE_DIR.iterdir() if file.suffix == '.csv'] + pool = Pool(PROCESS_COUNT) + data = pool.starmap(func, files) + pool.close() + pool.join() + return data + + +def generate_rl_start_values(): + data = np.concatenate(_mp_get_subset(_get_random_start_values_from_file, lines_per_file=15)) + print(data.shape) + pd.DataFrame(data).to_csv(RL_START_FILE, header=False, index=False) + + +def generate_non_touching_subset(): + data = np.concatenate(_mp_get_subset(_get_random_lines_from_file, lines_per_file=25)) + print(data.shape) + bars_below = (0, 2, 4) + bars_above = (1, 3) + # set output values = input values + data[:, SETUP_DIM + SENSOR_DIM:] = data[:, SETUP_DIM: SENSOR_DIM + SETUP_DIM] + # random bar setups without touching the tape + data[:, bars_below] = np.random.randint(0, 11, size=data[:, bars_below].shape) + data[:, bars_above] = np.random.randint(18, 45, size=data[:, bars_above].shape) + np.random.shuffle(data) + num_output_files = 12 + split_idx = [i * len(data) // num_output_files for i in range(num_output_files)] + for i in range(len(split_idx) - 1): + pd.DataFrame(data[split_idx[i]: split_idx[i+1]]).to_csv(NON_TOUCHING_DIR / f'Setup0_{i}.csv', + header=False, index=False) + pd.DataFrame(data[split_idx[-1]:]).to_csv(NON_TOUCHING_DIR / f'Setup0_{num_output_files - 1}.csv', + header=False, index=False) + + +if __name__ == '__main__': + generate_rl_start_values() diff --git a/preprocessing/tape_detection.py b/preprocessing/tape_detection.py new file mode 100644 index 0000000..97d615d --- /dev/null +++ b/preprocessing/tape_detection.py @@ -0,0 +1,17 @@ +import torch + + +def get_tape_edges(profiles, axis=1): + nonzeros = (profiles > torch.unsqueeze(torch.mean(profiles, dim=1), -1)) + tape_counter = nonzeros.cumsum(axis) + tape_counter[~nonzeros] = 0 + temp_tape_profile, temp_tape_start_idx = (tape_counter == 1).max(axis) + temp_tape_start_idx[temp_tape_profile == 0] = 0 + temp_tape_profile, temp_tape_end_idx = tape_counter.max(axis) + temp_tape_end_idx[temp_tape_profile == 0] = 0 + return temp_tape_start_idx, temp_tape_end_idx + + +def get_tape_width(profiles, axis=1): + temp_tape_start_idx, temp_tape_end_idx = get_tape_edges(profiles, axis) + return temp_tape_end_idx - temp_tape_start_idx diff --git a/process_control_model/__init__.py b/process_control_model/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/process_control_model/fixed_setups.py b/process_control_model/fixed_setups.py new file mode 100644 index 0000000..4ae0396 --- /dev/null +++ b/process_control_model/fixed_setups.py @@ -0,0 +1,82 @@ +import sys +from collections import OrderedDict +from pathlib import Path +import logging +import torch +import numpy as np +from preprocessing.split_datasets import get_test_set_respecting_files +from preprocessing.preprocessing_helper import BAR_CONFIG +from utils.tape_width_evaluation import evaluate_tape_width +from utils import logging_utils +from utils.logging_utils import set_multi_file_logging_parameters +from process_control_model.model_adapter import init_world_model +from process_control_model.rl_utils import MAX_ACTION, ACTION_RANGES, create_target +from utils.paths import NN_BACKEND, RF_BACKEND +from process_control_model.rl_parameters import TapeTargets +import process_control_model.rl_parameters as rn +import utils.paths as dirs + + +def test_fixed_setup(setup='Setup9'): + if setup is 'MAX': + actions = ACTION_RANGES[:, -1].float().to(device) + elif setup is 'MIN': + actions = ACTION_RANGES[:, 0].float().to(device) + else: + actions = torch.tensor(BAR_CONFIG[setup]).float().to(device) + test_x[:, :RL_PARAM_DICT[rn.setup_dim]] = actions / MAX_ACTION + resulting_tapes = world_model.predict(test_x.to(device)) + target = create_target(RL_PARAM_DICT[rn.target_width], RL_PARAM_DICT[rn.target_height], + dim=test_x.shape[0]) + evaluate_tape_width(target, resulting_tapes.cpu(), logger=LOGGER) + #choice = torch.randint(0, test_x.shape[0], (200,)) + #render(test_x[choice, RL_PARAM_DICT[rn.setup_dim]:].cpu(), resulting_tapes[choice].cpu(), + # actions[choice].cpu(), RL_PARAM_DICT[rn.target_width], RL_PARAM_DICT[rn.target_height]) + + +if __name__ == '__main__': + RL_PARAM_DICT = OrderedDict() + # default values + RL_PARAM_DICT[rn.random_seed] = 42 + RL_PARAM_DICT[rn.data_seed] = 1111 + RL_PARAM_DICT[rn.batch_size] = 64 + RL_PARAM_DICT[rn.sensor_dim] = 800 + RL_PARAM_DICT[rn.setup_dim] = 5 + # define target + if len(sys.argv) > 1: + tape_targets = TapeTargets(int(sys.argv[1]), 0.15) + else: + tape_targets = TapeTargets(240, 0.15) + RL_PARAM_DICT[rn.target_width] = tape_targets.target_width + RL_PARAM_DICT[rn.target_height] = tape_targets.target_height + # world model + RL_PARAM_DICT[rn.world_model_type] = 'nn' + RL_PARAM_DICT[rn.world_model_path] = NN_BACKEND if RL_PARAM_DICT[rn.world_model_type] is 'nn'\ + else RF_BACKEND + + torch.manual_seed(RL_PARAM_DICT[rn.random_seed]) + np.random.seed(RL_PARAM_DICT[rn.random_seed]) + device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu") + print('load world model') + world_model = init_world_model(RL_PARAM_DICT[rn.world_model_type], + RL_PARAM_DICT[rn.world_model_path], device) + print('load data') + test_x, _ = get_test_set_respecting_files(Path(dirs.DATA), + RL_PARAM_DICT[rn.sensor_dim], RL_PARAM_DICT[rn.setup_dim], random_seed=1111) + test_x = torch.from_numpy(test_x).float().to(device) + setups = range(1, 17) + for i in setups: + setup = f'Setup{i}' + RL_PARAM_DICT[rn.rl_model_path] = f'fixed_actions={setup}_on_' \ + f'{RL_PARAM_DICT[rn.world_model_type]}_target_w=' \ + f'{RL_PARAM_DICT[rn.target_width]}_' \ + f'target_h={RL_PARAM_DICT[rn.target_height]}' \ + f'_new_eval_simpleprep3' + RL_PARAM_DICT[rn.logging_path] = dirs.RL_LOGS + set_multi_file_logging_parameters(setup, RL_PARAM_DICT[rn.rl_model_path]) + LOGGER = logging.getLogger(setup) + print(RL_PARAM_DICT[rn.rl_model_path]) + LOGGER.info(f'Create {RL_PARAM_DICT[rn.rl_model_path]}') + rn.save_parameters(RL_PARAM_DICT, RL_PARAM_DICT[rn.rl_model_path]) + + test_fixed_setup(setup=setup) diff --git a/process_control_model/model_adapter.py b/process_control_model/model_adapter.py new file mode 100644 index 0000000..6c705c9 --- /dev/null +++ b/process_control_model/model_adapter.py @@ -0,0 +1,33 @@ +from joblib import load +import torch + +from process_control_model.rl_utils import load_pretrained_ff_model + + +def init_world_model(identifier, path, device): + if identifier == 'nn': + adapter = NNAdapter(path, device) + else: + adapter = RFAdapter(path, device) + return adapter + + +class RFAdapter: + def __init__(self, path, device): + self.model = load(path) + self.model.verbose = 0 + self.device = device + + def predict(self, data): + resulting_tapes = self.model.predict(data.detach().cpu().numpy()) + return torch.from_numpy(resulting_tapes).float().to(self.device) + + +class NNAdapter: + def __init__(self, path, device): + self.device = device + self._model = torch.load(path).to(device) + + def predict(self, data): + with torch.no_grad(): + return self._model(data) diff --git a/process_control_model/rl_parameters.py b/process_control_model/rl_parameters.py new file mode 100644 index 0000000..37b2ea7 --- /dev/null +++ b/process_control_model/rl_parameters.py @@ -0,0 +1,69 @@ +import json +from collections import namedtuple +from pathlib import Path +import utils.paths as dirs + +random_seed = 'random_seed' +data_seed = 'data_seed' +sensor_dim = 'sensor_dim' +setup_dim = 'setup_dim' +target_width = 'target_width' +target_height = 'target_height' +k_h = 'k_h' +k_w = 'k_w' +k_m = 'k_m' +world_model_type = 'world_model_type' +world_model_path = 'world_model_path' +hidden_dims = 'hidden_dims' +learning_rate = 'learning_rate' +batch_size = 'batch_size' +rl_model_path = 'rl_model_path' +logging_path = 'logging_path' +critic_architecture = 'critic_architecture' +actor_architecture = 'actor_architecture' +population_size = 'population_size' +ga_noise = 'ga_noise' +reproducing_parents = 'reproducing_parents' + + +def create_rl_model_path(prefix, param_dict): + return '_'.join([prefix] + [str(i) for i in param_dict[hidden_dims]] + + [param_dict[world_model_type], f'lr={param_dict[learning_rate]}', + f'bs={param_dict[batch_size]}', + 'rescaled_actions', f'target_w={param_dict[target_width]}', + f'target_h={param_dict[target_height]}', f'k_h={param_dict[k_h]}', + f'k_w={param_dict[k_w]}', + f'k_m={param_dict[k_m]}']) + + +def create_ga_model_path(prefix, param_dict): + return '_'.join([prefix] + [str(i) for i in param_dict[hidden_dims]] + + [param_dict[world_model_type], + f'population_size={param_dict[population_size]}', + f'noise={param_dict[ga_noise]}', + f'rep_parents={param_dict[reproducing_parents]}', + f'target_w={param_dict[target_width]}', + f'target_h={param_dict[target_height]}', + f'k_h={param_dict[k_h]}', f'k_w={param_dict[k_w]}', + f'k_m={param_dict[k_m]}']) + + +RewardFactors = namedtuple('RewardFactors', ['height_factor', 'width_factor', 'movement_factor']) +TapeTargets = namedtuple('TapeTargets', ['target_width', 'target_height']) + + +PARAMETERS_DIR = Path(dirs.RL_PARAMS) + + +def save_parameters(param_dict, filename): + if not PARAMETERS_DIR.exists(): + PARAMETERS_DIR.mkdir() + with (PARAMETERS_DIR / (filename + '.json')).open(mode='w') as output_file: + json.dump(param_dict, output_file) + + +def load_parameters(filename): + rl_param_dict = json.load((PARAMETERS_DIR / Path(filename)).open()) + tape_targets = TapeTargets(rl_param_dict[target_width], rl_param_dict[target_height]) + reward_factors = RewardFactors(rl_param_dict[k_h], rl_param_dict[k_w], rl_param_dict[k_m]) + return rl_param_dict, reward_factors, tape_targets diff --git a/process_control_model/rl_utils.py b/process_control_model/rl_utils.py new file mode 100644 index 0000000..ee38861 --- /dev/null +++ b/process_control_model/rl_utils.py @@ -0,0 +1,79 @@ +import numpy as np +import torch +from torch.utils.data import Dataset +from process_model.feedforward import Feedforward +from preprocessing.tape_detection import get_tape_edges +from preprocessing.preprocessing_helper import MAX_ACTION, ACTION_RANGES + + +def create_target(target_width, target_height, dim=1, sensor_dim=800): + target = torch.zeros((dim, sensor_dim)) + begin_target_tape = sensor_dim // 2 - target_width // 2 + target[:, begin_target_tape: begin_target_tape + target_width] = target_height + return target + + +def load_pretrained_ff_model(path, device, hidden_dims=(300, 300, 300)): + ff_model = Feedforward(805, hidden_dims, 800) + ff_model.load_state_dict(torch.load(path, map_location=device)) + ff_model.to(device) + for param in ff_model.parameters(): + param.requires_grad = False + ff_model.eval() + return ff_model + + +class RLDataSet(Dataset): + def __init__(self, observations): + self._observations = observations + + def __len__(self): + return self._observations.shape[0] + + def __getitem__(self, item): + return self._observations[item] + + +def clear_setup(data): + data[:, :5] = np.array([17., 11., 22., 11., 17.]) / MAX_ACTION + + +def _get_height_difference(rows, edges, target_height, e_min=0, e_max=799, min_len=180): + edges[0][edges[0] < e_min] = e_min + edges[1][edges[0] == edges[1]] += min_len + edges[1][edges[1] > e_max] = e_max + device = 'cuda:0' if rows.is_cuda else 'cpu' + avg_heights = torch.empty(rows.shape[0]).to(device) + for i, row in enumerate(rows): + avg_heights[i] = torch.mean(row[edges[0][i]: edges[1][i]]) + return -torch.abs(avg_heights - target_height) + + +def _get_width_difference(edges, target_width, e_min=0, e_max=799, max_width=800): + edges[0][edges[0] < e_min] = e_min + edges[1][edges[1] > e_max] = e_max + return -(abs((edges[1] - edges[0]) - target_width).float() / max_width) + + +def _get_bar_movement(actions, previous_actions, max_movement=250): + total_movement = torch.sum(torch.abs(previous_actions - actions), dim=1) + return -total_movement / max_movement + + +def get_reward(resulting_tapes, actions, previous_actions, reward_factors, tape_targets): + edges = get_tape_edges(resulting_tapes) + y_height = _get_height_difference(resulting_tapes, edges, tape_targets.target_height) + y_width = _get_width_difference(edges, tape_targets.target_width) + y_bar_movement = _get_bar_movement(actions, previous_actions) + y = reward_factors.height_factor * y_height +\ + reward_factors.width_factor * y_width +\ + reward_factors.movement_factor * y_bar_movement + if torch.isnan(y).any(): + print('NAN') + return y + + +def transform_actions_prob(actions): + device = 'cuda:0' if actions.is_cuda else 'cpu' + actions *= (ACTION_RANGES[:, 1] - ACTION_RANGES[:, 0]).to(device) + actions += ACTION_RANGES[:, 0].to(device) diff --git a/process_control_model/start_neuroevolution.py b/process_control_model/start_neuroevolution.py new file mode 100644 index 0000000..2d89218 --- /dev/null +++ b/process_control_model/start_neuroevolution.py @@ -0,0 +1,254 @@ +import sys +from pathlib import Path +import copy +from collections import OrderedDict +import logging +import torch +import torch.nn as nn +from torch.utils.tensorboard import SummaryWriter +import numpy as np + +from preprocessing.split_datasets import split_data_respecting_files, get_test_set_respecting_files +from utils.tape_width_evaluation import evaluate_tape_width +import utils.logging_utils as logging_config +from utils.nn_utils import count_trainable_parameters +from process_control_model.model_adapter import init_world_model +from process_control_model.rl_utils import MAX_ACTION, create_target, clear_setup, \ + transform_actions_prob, get_reward +from process_control_model.rl_parameters import RewardFactors, TapeTargets +import process_control_model.rl_parameters as rp +from utils.paths import NN_BACKEND, RF_BACKEND +from plotter.rl_plotter import render +import utils.paths as dirs + + +class Net(nn.Module): + def __init__(self, state_dim, action_dim, hidden_dims): + super(Net, self).__init__() + act = nn.ReLU() + layers = OrderedDict() + layers['layer_0'] = nn.Linear(state_dim, hidden_dims[0]) + layers['act_0'] = act + for i, h_d in enumerate(hidden_dims[1:]): + layers[f'layer_{i + 1}'] = nn.Linear(hidden_dims[i], h_d) + layers[f'act_{i + 1}'] = act + layers[f'layer_{len(hidden_dims)}'] = nn.Linear(hidden_dims[-1], action_dim) + layers[f'act_{len(hidden_dims)}'] = nn.Sigmoid() + self.base = nn.Sequential(layers) + + def forward(self, x): + return self.base(x) + + +def create_new_population(): + networks = [Net(RL_PARAM_DICT[rp.sensor_dim] + RL_PARAM_DICT[rp.setup_dim], + RL_PARAM_DICT[rp.setup_dim], + RL_PARAM_DICT[rp.hidden_dims]).to(device) + for _ in range(RL_PARAM_DICT[rp.population_size])] + print('validate') + select_data = np.random.choice(train_x.shape[0], 512, replace=False) + population = [(determine_fitness(net, train_x[select_data]), net) + for net in networks] + return population + + +def rank_genotypes(population): + population.sort(key=lambda p: p[0], reverse=True) + + +def reduce_mutation_rate(factor=0.6, min_rate=0.01): + RL_PARAM_DICT[rp.ga_noise] *= factor + if RL_PARAM_DICT[rp.ga_noise] < min_rate: + RL_PARAM_DICT[rp.ga_noise] = min_rate + + +def save_genotype(genotype): + if not dirs.RL_WEIGHTS.exists(): + dirs.RL_WEIGHTS.mkdir() + torch.save(genotype.state_dict(), + dirs.RL_WEIGHTS / (RL_PARAM_DICT[rp.rl_model_path] + '.pth')) + if not dirs.RL_MODELS.exists(): + dirs.RL_MODELS.mkdir() + torch.save(genotype, dirs.RL_MODELS / (RL_PARAM_DICT[rp.rl_model_path] + '.pth')) + + +def set_up_first_generation(): + population = create_new_population() + rank_genotypes(population) + best_fitness = determine_fitness(population[0][1], val_x, 0).item() + logger.info(str(population[0][1])) + logger.info(f'Trainable Parameters: {count_trainable_parameters(population[0][1])}') + return population, best_fitness + + +def start_evolution(): + population, best_fitness = set_up_first_generation() + patience = 10 + stagnating_epochs = 0 + for i in range(1000): + best_parent = population[0] + select_parents = np.random.choice(RL_PARAM_DICT[rp.reproducing_parents], len(population)) + mutants = [mutate(population[j][1]) for j in select_parents] + select_data = np.random.choice(train_x.shape[0], 512, replace=False) + population = [(determine_fitness(net, train_x[select_data]), + net) for net in mutants] + population.append(best_parent) + rank_genotypes(population) + print(f'generation {i + 1}') + print(f'best on train: {population[0][0].item()}') + logger.info(f'generation {i + 1}') + logger.info(f'best on train: {population[0][0].item()}') + writer.add_scalar('train reward', population[0][0].item(), i) + population_fitness = determine_fitness(population[0][1], val_x, i, True).item() + print(f'best on validation: {population_fitness}') + logger.info(f'best on validation: {population_fitness}') + writer.add_scalar('validation reward', population_fitness, i) + if population_fitness > best_fitness: + save_genotype(population[0][1]) + print('save states') + stagnating_epochs = 0 + best_fitness = population_fitness + else: + stagnating_epochs += 1 + if stagnating_epochs > patience: + print(f'early stopping kicked in after {patience} epochs without improvement') + logger.info(f'early stopping kicked in after {patience} epochs without improvement') + test_fitness = determine_fitness(population[0][1], test_x, -1, True, False).item() + print(f'best on test: {test_fitness}') + logger.info(f'best on test: {test_fitness}') + break + if i % 5 == 0: + reduce_mutation_rate() + logger.info('reduced mutation rate') + print('reduced mutation rate') + + +def determine_fitness(net, data, generation=1, evaluate_resulting_tapes=False, plot_results=False): + with torch.no_grad(): + actions = net(data) + scaled_actions = actions.clone().detach() + transform_actions_prob(scaled_actions) + normed_actions = scaled_actions.clone() / MAX_ACTION + world_model_input = data.clone().detach() + world_model_input[:, :RL_PARAM_DICT[rp.setup_dim]] = normed_actions + resulting_tapes = world_model.predict(world_model_input) + rewards = get_reward(resulting_tapes.detach(), scaled_actions, + data[:, :RL_PARAM_DICT[rp.setup_dim]].detach() * MAX_ACTION, + REWARD_FACTORS, TAPE_TARGETS) + if evaluate_resulting_tapes: + target = create_target(RL_PARAM_DICT[rp.target_width], + RL_PARAM_DICT[rp.target_height], + dim=resulting_tapes.shape[0]) + avg_tape_width_dev = evaluate_tape_width(target, resulting_tapes.cpu()) + if generation >= 0: + writer.add_scalar('validation width ratio', avg_tape_width_dev, generation) + if plot_results: + choice = torch.randint(0, data.shape[0], (200,)) + render(data[choice, RL_PARAM_DICT[rp.setup_dim]:].cpu(), + resulting_tapes[choice].cpu(), + scaled_actions[choice].cpu(), + RL_PARAM_DICT[rp.target_width], + RL_PARAM_DICT[rp.target_height]) + + return torch.mean(rewards) + + +def mutate(parent_net): + net = copy.deepcopy(parent_net) + for p in net.parameters(): + p.data += RL_PARAM_DICT[rp.ga_noise] * torch.randn_like(p.data) + return net + + +def test_trained_model(): + if (dirs.RL_WEIGHTS / (RL_PARAM_DICT[rp.rl_model_path] + '.pth')).exists(): + net = torch.load(dirs.RL_MODELS / (RL_PARAM_DICT[rp.rl_model_path] + '.pth')) + else: + net = Net(RL_PARAM_DICT[rp.sensor_dim] + RL_PARAM_DICT[rp.setup_dim], RL_PARAM_DICT[rp.setup_dim], + RL_PARAM_DICT[rp.hidden_dims]) + net.load_state_dict(torch.load(dirs.RL_WEIGHTS / (RL_PARAM_DICT[rp.rl_model_path] + '.pth'))) + net.to(device) + test_fitness = determine_fitness(net, test_x, -1, True, True).item() + print(f'reward on test set: {test_fitness}') + logger.info(f'reward on test set: {test_fitness}') + + +def set_parameters(): + rl_param_dict = OrderedDict() + # default values + rl_param_dict[rp.random_seed] = 42 + rl_param_dict[rp.data_seed] = 1111 + rl_param_dict[rp.sensor_dim] = 800 + rl_param_dict[rp.setup_dim] = 5 + # define target + tape_targets = TapeTargets(240, 0.15) + rl_param_dict[rp.target_width] = tape_targets.target_width + rl_param_dict[rp.target_height] = tape_targets.target_height + # specify reward function + reward_factors = RewardFactors(1.0, 4.0, 0.5) + rl_param_dict[rp.k_h] = reward_factors.height_factor + rl_param_dict[rp.k_w] = reward_factors.width_factor + rl_param_dict[rp.k_m] = reward_factors.movement_factor + # world model parameters + rl_param_dict[rp.world_model_type] = 'nn' + rl_param_dict[rp.world_model_path] = NN_BACKEND if rl_param_dict[rp.world_model_type]\ + == 'nn' else RF_BACKEND + # ga parameters + rl_param_dict[rp.hidden_dims] = [128, 16] + # population size > reproducing parents + rl_param_dict[rp.population_size] = 300 + rl_param_dict[rp.reproducing_parents] = 15 + rl_param_dict[rp.ga_noise] = 0.2 + rl_param_dict[rp.rl_model_path] = rp.create_ga_model_path( + 'ga_wo_env_first_step_noise_reduction_new_eval_simpleprep3_small_backend', rl_param_dict) + rl_param_dict[rp.logging_path] = dirs.RL_LOGS + print(rl_param_dict[rp.rl_model_path]) + rp.save_parameters(rl_param_dict, rl_param_dict[rp.rl_model_path]) + return rl_param_dict, reward_factors, tape_targets + + +if __name__ == '__main__': + TRAINING = True + if len(sys.argv) > 1: + RL_PARAM_DICT, REWARD_FACTORS, TAPE_TARGETS = rp.load_parameters(sys.argv[1]) + else: + RL_PARAM_DICT, REWARD_FACTORS, TAPE_TARGETS = set_parameters() + torch.manual_seed(RL_PARAM_DICT[rp.random_seed]) + np.random.seed(RL_PARAM_DICT[rp.random_seed]) + device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu") + print('load world model') + world_model = init_world_model(RL_PARAM_DICT[rp.world_model_type], + RL_PARAM_DICT[rp.world_model_path], device) + print('load data') + if TRAINING: + train_x, test_x, val_x = split_data_respecting_files(Path(dirs.DATA), + RL_PARAM_DICT[rp.sensor_dim], RL_PARAM_DICT[rp.setup_dim], + random_seed=RL_PARAM_DICT[rp.data_seed])[:3] + clear_setup(train_x) + clear_setup(test_x) + clear_setup(val_x) + train_x = torch.from_numpy(train_x).float().to(device) + test_x = torch.from_numpy(test_x).float().to(device) + val_x = torch.from_numpy(val_x).float().to(device) + else: + test_x = get_test_set_respecting_files(Path(dirs.DATA), + RL_PARAM_DICT[rp.sensor_dim], RL_PARAM_DICT[rp.setup_dim], + random_seed=RL_PARAM_DICT[rp.data_seed])[0] + clear_setup(test_x) + test_x = torch.from_numpy(test_x).float().to(device) + + logging_config.set_rl_logging_parameters(RL_PARAM_DICT[rp.rl_model_path]) + logger = logging.getLogger(__name__) + + logger.info(f'target_width={RL_PARAM_DICT[rp.target_width]}') + logger.info(f'target_height={RL_PARAM_DICT[rp.target_height]}') + + logger.info(f'reward factors: height_diff * {RL_PARAM_DICT[rp.k_h]} + ' + f'width_diff * {RL_PARAM_DICT[rp.k_w]}' + f' + movement * {RL_PARAM_DICT[rp.k_m]}') + logger.info(f'Create {RL_PARAM_DICT[rp.rl_model_path]} on {device}') + if TRAINING: + writer = SummaryWriter(Path(dirs.RL_TENSORBOARD, RL_PARAM_DICT[rp.rl_model_path])) + start_evolution() + else: + test_trained_model() diff --git a/process_model/__init__.py b/process_model/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/process_model/feedforward.py b/process_model/feedforward.py new file mode 100644 index 0000000..b248afa --- /dev/null +++ b/process_model/feedforward.py @@ -0,0 +1,21 @@ +import torch.nn as nn + + +class Feedforward(nn.Module): + def __init__(self, input_size, hidden_dimensions, output_size): + super(Feedforward, self).__init__() + assert len(hidden_dimensions) > 0 + self.input_layer = nn.Linear(input_size, hidden_dimensions[0]) + self.input_activation = nn.ReLU() + self.hidden_layers = nn.ModuleList([nn.Linear(hidden_dimensions[i], h_size) + for i, h_size in enumerate(hidden_dimensions[1:])]) + self.hidden_activations = nn.ModuleList([nn.ReLU() for _ in hidden_dimensions[1:]]) + self.output_layer = nn.Linear(hidden_dimensions[-1], output_size) + + def forward(self, x): + x = self.input_activation(self.input_layer(x)) + for i, layer in enumerate(self.hidden_layers): + x = self.hidden_activations[i](layer(x)) + x = self.output_layer(x) + return x + diff --git a/process_model/ff_trainer.py b/process_model/ff_trainer.py new file mode 100644 index 0000000..03988bb --- /dev/null +++ b/process_model/ff_trainer.py @@ -0,0 +1,26 @@ +from pathlib import Path +import torch +from torch.optim import lr_scheduler +from process_model.feedforward import Feedforward +from process_model.trainer import Trainer + + +class FFTrainer(Trainer): + def __init__(self, sensor_dim=800, param_dim=5, + hidden_dimensions=(300, 300, 300), learning_rate=0.001, profile_only_input=False, + model_path=Path('model.pth')): + # network + if profile_only_input: + model = Feedforward(sensor_dim, hidden_dimensions, sensor_dim) + else: + model = Feedforward(sensor_dim + param_dim, + hidden_dimensions, sensor_dim) + optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate) + scheduler = lr_scheduler.ReduceLROnPlateau( + optimizer, + factor=0.5, + patience=5, + verbose=True, + threshold=1e-5, + min_lr=1e-6) + super(FFTrainer, self).__init__(model, sensor_dim, model_path, optimizer, scheduler) diff --git a/process_model/model_comparison.py b/process_model/model_comparison.py new file mode 100644 index 0000000..381b57e --- /dev/null +++ b/process_model/model_comparison.py @@ -0,0 +1,103 @@ +from pathlib import Path +import torch +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import rc + +from preprocessing.preprocessing_helper import ACTION_RANGES +from preprocessing.tape_detection import get_tape_edges +from preprocessing.split_datasets import get_test_set_respecting_files +from utils.logging_utils import get_best_results, get_results_from_logs +import utils.paths as dirs +from process_control_model.model_adapter import NNAdapter, RFAdapter +from process_control_model.rl_utils import transform_actions_prob + + +def _get_best_models(): + nn_best_mean, nn_best_mean_std = get_best_results(get_results_from_logs( + ['relu', 'fixed_seed'], path=dirs.NN_LOGS)) + nn_model = NNAdapter(dirs.NN_MODELS / (nn_best_mean_std.filename + '.pth'), 'cpu') + rf_best_mean, rf_best_mean_std = get_best_results(get_results_from_logs( + ['random_forest'], path=dirs.RF_LOGS)) + rf_model = RFAdapter(dirs.RF_MODELS / (rf_best_mean_std.filename + '.joblib'), 'cpu') + return nn_model, rf_model + + +def _get_test_tensors(num_samples=10): + test_x, test_y = get_test_set_respecting_files(Path(dirs.DATA), 800, 5, 1111) + test_x = torch.from_numpy(test_x).float() + test_y = torch.from_numpy(test_y).float() + select_data = np.random.choice(test_x.shape[0], num_samples, replace=False) + return test_x[select_data], test_y[select_data] + + +def plot_samples(samples_list, legend, title, plot_edges=False): + if plot_edges: + edges = [get_tape_edges(samples_sublist) for samples_sublist in samples_list] + num_samples = samples_list[0].shape[0] + + font_size = 22 + # rc('text', usetex=True) + # rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']}) + plt.rc('font', size=font_size) # controls default text sizes + plt.rc('axes', titlesize=font_size) # fontsize of the axes title + plt.rc('axes', labelsize=font_size) # fontsize of the x and y labels + plt.rc('xtick', labelsize=font_size-4) # fontsize of the tick labels + plt.rc('ytick', labelsize=font_size-4) # fontsize of the tick labels + plt.rc('legend', fontsize=font_size) # legend fontsize + plt.rc('figure', titlesize=font_size) # fontsize of the figure title + fig, axs = plt.subplots(num_samples, sharex=True, sharey=True) + for i, samples_sublist in enumerate(samples_list): + plot_i = 1 + for j in range(num_samples): + # plt.subplot(num_samples, 1, plot_i) + if plot_edges: + axs[j].plot(samples_sublist[j], '-o', linewidth=1.4, + markevery=[edges[i][0][j], edges[i][1][j]]) + # plt.plot(samples_sublist[j], '-o', linewidth=0.8, + # markevery=[edges[i][0][j], edges[i][1][j]]) + else: + # plt.plot(samples_sublist[j], linewidth=0.8) + axs[j].plot(samples_sublist[j], linewidth=1.4) + plot_i += 1 + fig.text(0.5, 0.045, 'Sensor pixels', ha='center') + fig.text(0.08, 0.5, 'Tow height', va='center', rotation='vertical') + plt.legend(legend, loc='right') + ttl = fig.suptitle(title) + ttl.set_position([.5, 0.92]) + plt.show() + + +def compare_random_actions(num_samples=10): + nn_model, rf_model = _get_best_models() + test_x, test_y = _get_test_tensors(num_samples=num_samples) + # random actions + # actions = torch.rand_like(test_x[:, :5]) + # transform_actions_prob(actions) + # or max actions to see the differences + actions = torch.empty_like(test_x[:, :5]) + actions[:, 0::2] = ACTION_RANGES[0::2, 1] + actions[:, 1::2] = ACTION_RANGES[1::2, 0] + # or min actions to see the differences + # actions = torch.empty_like(test_x[:, :5]) + # actions[:, 0::2] = ACTION_RANGES[0::2, 0] + # actions[:, 1::2] = ACTION_RANGES[1::2, 1] + test_x[:, :5] = actions / 60. + nn_preds = nn_model.predict(test_x) + rf_preds = rf_model.predict(test_x) + plot_samples([nn_preds, rf_preds], legend=['ffnn_predictions', 'rf_predictions'], + title='Process models - maximum spreading', plot_edges=True) + + +def compare_on_real_data(num_samples=10): + nn_model, rf_model = _get_best_models() + test_x, test_y = _get_test_tensors(num_samples=num_samples) + nn_preds = nn_model.predict(test_x) + rf_preds = rf_model.predict(test_x) + plot_samples([nn_preds, rf_preds, test_y], legend=['ffnn_predictions', 'rf_predictions', 'target'], + title='Process models', plot_edges=False) + + +if __name__ == '__main__': + # compare_on_real_data() + compare_random_actions() diff --git a/process_model/random_forest_regression.py b/process_model/random_forest_regression.py new file mode 100644 index 0000000..dfe8beb --- /dev/null +++ b/process_model/random_forest_regression.py @@ -0,0 +1,49 @@ +import logging +from pathlib import Path +from sklearn.ensemble import RandomForestRegressor, ExtraTreesRegressor +from sklearn.metrics import mean_squared_error as mse +from joblib import dump, load +from utils.tape_width_evaluation import evaluate_tape_width_np_wrapper +import utils.paths as dirs + + +LOGGER = logging.getLogger(__name__) + + +def load_model(name): + return load(dirs.RF_MODELS / (name + '.joblib')) + + +def save_model(reg, name): + if not dirs.RF_MODELS.exists(): + dirs.RF_MODELS.mkdir() + dump(reg, dirs.RF_MODELS / (name + '.joblib')) + + +def train_rf(train_x, train_y, estimators=10, min_samples_split=2, min_samples_leaf=1, + max_depth=None, random_state=1111, max_features='auto', jobs=8): + reg = RandomForestRegressor(n_estimators=estimators, min_samples_split=min_samples_split, + min_samples_leaf=min_samples_leaf, max_depth=max_depth, + n_jobs=jobs, verbose=20, random_state=random_state, + max_features=max_features, bootstrap=True) + LOGGER.info(str(reg)) + reg.fit(train_x, train_y) + return reg + + +def train_extra_trees(train_x, train_y, estimators=10, min_samples_split=2, min_samples_leaf=1, + max_depth=None, random_state=1111, max_features='auto', jobs=8): + reg = ExtraTreesRegressor(n_estimators=estimators, min_samples_split=min_samples_split, + min_samples_leaf=min_samples_leaf, max_depth=max_depth, n_jobs=jobs, + verbose=20, random_state=random_state, max_features=max_features, + bootstrap=True) + LOGGER.info(str(reg)) + reg.fit(train_x, train_y) + return reg + + +def validate_rf(reg, val_x, val_y, plot_stats=False): + preds = reg.predict(val_x) + tape_width_dev = evaluate_tape_width_np_wrapper(val_y, preds, plot_stats) + LOGGER.info(f'tape width dev: {tape_width_dev:.4f}') + return mse(val_y, preds), tape_width_dev, preds diff --git a/process_model/start_ff.py b/process_model/start_ff.py new file mode 100644 index 0000000..37e2e4c --- /dev/null +++ b/process_model/start_ff.py @@ -0,0 +1,70 @@ +import sys +import logging +from pathlib import Path +import torch +import numpy as np +from preprocessing.split_datasets import split_data_respecting_files +from process_model.ff_trainer import FFTrainer +from preprocessing.datasetloader import load_data_and_meta, load_test_data_and_meta +from plotter.profile_plotter import plot_samples +from utils.nn_utils import count_trainable_parameters +from utils.logging_utils import set_logging_parameters +import utils.paths as dirs + +# execute in adaptive-spreading directory: python3 -m process_model.start_ff + + +def _start_ff(): + set_logging_parameters(MODEL_NAME) + logger = logging.getLogger(__name__) + logger.info(f'Create {MODEL_NAME} on {DEVICE}') + trainer = FFTrainer(sensor_dim=SENSOR_DIM, param_dim=PARAMS_DIM, hidden_dimensions=HIDDEN_DIMS, + learning_rate=LEARNING_RATE, + profile_only_input=PROFILE_ONLY_INPUT, model_path=MODEL_PATH) + logger.info(str(trainer.model)) + trainer.model.to(DEVICE) + logger.info(f'Trainable Parameters: {count_trainable_parameters(trainer.model)}') + logger.info(f'Splitting strategy: {SPLIT_FUNCTION.__name__} {SPLIT_PARAM}') + args = [DATA_PATH, SENSOR_DIM, PARAMS_DIM, PROFILE_ONLY_INPUT] + if SPLIT_PARAM: + args.append(SPLIT_PARAM) + if IS_TRAINING: + train_loader, validation_loader, test_loader = load_data_and_meta( + DEVICE, BATCH_SIZE, *args) + trainer.train_network(train_loader, validation_loader) + else: + test_loader = load_test_data_and_meta(DEVICE, DATA_PATH, SENSOR_DIM, PARAMS_DIM, + SPLIT_PARAM) + _, all_targets, all_preds, all_filenames, all_indices = trainer.test_network(test_loader) + # plot_samples(MODEL_NAME, SENSOR_DIM, (all_preds, all_targets), + # (all_filenames, all_indices)) + + +if __name__ == '__main__': + IS_TRAINING = False + PROFILE_ONLY_INPUT = False + if len(sys.argv) > 1: + HIDDEN_DIMS = [int(i) for i in sys.argv[1].strip('[]').split(',')] + else: + HIDDEN_DIMS = [100, 100, 100] + DATA_PATH = Path(dirs.DATA) + SPLIT_FUNCTION = split_data_respecting_files + SPLIT_PARAM = 1111 + SENSOR_DIM = 800 + PARAMS_DIM = 5 + LEARNING_RATE = 5e-04 + BATCH_SIZE = 64 + torch.manual_seed(42) + np.random.seed(SPLIT_PARAM) + MODEL_NAME = '_'.join([str(i) for i in HIDDEN_DIMS] + ['relu', f'lr={LEARNING_RATE}', + f'batch_size={BATCH_SIZE}', + SPLIT_FUNCTION.__name__, + 'new_eval_fixed_seed']) + MODEL_NAME = '100_100_100_relu_lr=0.0001_batch_size=64_split_data_respecting_files_new_eval' + if SPLIT_PARAM: + MODEL_NAME += f'_param={SPLIT_PARAM}' + if PROFILE_ONLY_INPUT: + MODEL_NAME += '_profile_only' + DEVICE = torch.device("cuda:0" if torch.cuda.is_available() else "cpu") + MODEL_PATH = dirs.NN_WEIGHTS / f'{MODEL_NAME}.pth' + _start_ff() diff --git a/process_model/start_rf.py b/process_model/start_rf.py new file mode 100644 index 0000000..e7cf037 --- /dev/null +++ b/process_model/start_rf.py @@ -0,0 +1,56 @@ +import logging +from pathlib import Path +from preprocessing.split_datasets import split_data_respecting_files_including_meta, \ + get_test_set_including_meta +from process_model.random_forest_regression import train_rf, validate_rf, save_model, load_model +from plotter.profile_plotter import plot_samples_np_wrapper +from utils.logging_utils import set_logging_parameters +import utils.paths as dirs + + +def _start_rf(): + set_logging_parameters(MODEL_NAME) + logger = logging.getLogger(__name__) + logger.info(MODEL_NAME) + logger.info('Load data...') + if IS_TRAINING: + train_data, test_data, val_data = split_data_respecting_files_including_meta( + DATA_PATH, SENSOR_DIM, PARAMS_DIM, PROFILE_ONLY_INPUT, RANDOM_SEED) + logger.info( + f'Loaded! Trainset size: {train_data[0].shape}, Test set size: {test_data[0].shape}') + logger.info('Build random forest...') + reg = train_rf(train_data[0], train_data[1], random_state=RANDOM_SEED, + estimators=ESTIMATORS, min_samples_split=MIN_SAMPLES_SPLIT, + max_features=MAX_FEAT, max_depth=MAX_DEPTH, jobs=10) + save_model(reg, MODEL_NAME) + mse, _, _ = validate_rf(reg, train_data[0], train_data[1]) + logger.info(f'Train loss: {mse}') + mse, _, preds = validate_rf(reg, val_data[0], val_data[1]) + logger.info(f'Val loss: {mse}') + else: + test_data = get_test_set_including_meta( + DATA_PATH, SENSOR_DIM, PARAMS_DIM, RANDOM_SEED) + logger.info(f'Loaded! Test set size: {test_data[0].shape}') + logger.info('Load pretrained model...') + reg = load_model(MODEL_NAME) + mse, _, preds = validate_rf(reg, test_data[0], test_data[1], plot_stats=False) + logger.info(f'Test loss: {mse}') + plot_samples_np_wrapper(preds, test_data[1], SENSOR_DIM, MODEL_NAME, + meta=(test_data[2], test_data[3])) + + +if __name__ == '__main__': + IS_TRAINING = False + DATA_PATH = Path(dirs.DATA) + SENSOR_DIM = 800 + PARAMS_DIM = 5 + PROFILE_ONLY_INPUT = False + RANDOM_SEED = 1111 + ESTIMATORS = 50 + MAX_DEPTH = 12 + MIN_SAMPLES_SPLIT = 150 + MAX_FEAT = 300 + MODEL_NAME = f'est_{ESTIMATORS}_maxdepth_{MAX_DEPTH}_min_split_{MIN_SAMPLES_SPLIT}_' \ + f'max_feats_{MAX_FEAT}_random_forest_regressor_new_eval' + MODEL_NAME = 'est_50_maxdepth_15_min_split_150_max_feats_200_random_forest_regressor_new_eval' + _start_rf() diff --git a/process_model/trainer.py b/process_model/trainer.py new file mode 100644 index 0000000..a936b51 --- /dev/null +++ b/process_model/trainer.py @@ -0,0 +1,117 @@ +import logging +from pathlib import Path +import torch +import torch.nn as nn +from torch.utils.tensorboard import SummaryWriter +from utils.tape_width_evaluation import evaluate_tape_width +import utils.paths as dirs + + +class Trainer: + def __init__(self, model, output_dim, model_path, optimizer, lr_scheduler=None): + self.output_dim = output_dim + self.model = model + self.model_path = model_path + self.optimizer = optimizer + self.scheduler = lr_scheduler + self.patience = 10 + self.num_stagnating_epochs = 0 + self.lowest_error = 10000 + # tensorboard + self.writer = SummaryWriter(Path(dirs.NN_TENSORBOARD, self.model_path.name)) + self.logger = logging.getLogger(__name__) + + def save_model_weights(self): + if not self.model_path.parents[0].exists(): + self.model_path.parents[0].mkdir() + torch.save(self.model.state_dict(), self.model_path) + if not dirs.NN_MODELS.exists(): + dirs.NN_MODELS.mkdir() + torch.save(self.model, dirs.NN_MODELS / self.model_path.name) + + def early_stopper(self, val_loss): + """ + Returns flag to stop training when validation loss did not decrease over the last + epochs (patience). Saves the network parameters after the best epoch (lowest + validation loss) + :param val_loss: Loss on validation set during current epoch + :return: True to stop, False to continue training + """ + if val_loss < self.lowest_error: + self.lowest_error = val_loss + self.save_model_weights() + self.logger.info('save states') + self.num_stagnating_epochs = 0 + else: + self.num_stagnating_epochs += 1 + if self.num_stagnating_epochs >= self.patience: + return True + return False + + def train_network(self, train_loader, val_loader): + criterion = nn.MSELoss() + loss_history, val_loss_history = [], [] + for i in range(2000): + losses = 0 + self.model.train() + for sample, target, _, _ in train_loader: + self.optimizer.zero_grad() + output = self.model(sample) + loss = criterion(output, target) + loss.backward() + self.optimizer.step() + losses += loss.item() + losses /= len(train_loader) + self.writer.add_scalar("training_loss", losses, i) + self.logger.info(f'Training loss after epoch {i + 1}: {losses}') + loss_history.append(losses) + val_loss, _, _, _, _ = self.test_network(val_loader, True, i) + val_loss_history.append(val_loss) + self.writer.add_scalar("validation_loss", val_loss, i) + + if self.early_stopper(val_loss): + break + if self.scheduler: + self.scheduler.step(val_loss) + + return loss_history, val_loss_history + + def test_network(self, loader, validating=False, epoch=0): + """ + Evaluates network using the given data. + :param loader: Dataloader containing validation or testing data + :param validating: Indicates if the function is called for validating or testing. + :param epoch: epoch counter when validating + :return: Accumulated loss on the given data, target and prediction tensor + """ + if not validating: + # for testing: restore parameters with best performance on + # validation set + self.model.load_state_dict(torch.load(self.model_path)) + self.model.eval() + loader_size = len(loader) + + all_preds = torch.empty( + (loader_size, loader.batch_size, self.output_dim)) + all_targets = torch.empty_like(all_preds) + i = 0 + # evaluate network without gradient computation + with torch.no_grad(): + criterion = nn.MSELoss() + for sample, target, _, _ in loader: + pred = self.model(sample) + all_preds[i] = pred + all_targets[i] = target + i += 1 + + loss = criterion(all_preds, all_targets).item() + tape_width_dev = evaluate_tape_width(all_targets, all_preds, plot_stats=False) + self.logger.info(f'tape width dev: {tape_width_dev:.4f}') + _, _, all_filenames, all_indices = loader.dataset[:] + if validating: + self.logger.info(f'Val Loss: {loss}') + self.writer.add_scalar("avg_width_dev", tape_width_dev, epoch) + else: + self.logger.info(f'Test Loss: {loss}') + + return loss, all_targets, all_preds, all_filenames, all_indices diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..0889e86 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,8 @@ +matplotlib==3.0.2 +vispy==0.6.3 +pandas==0.25.3 +joblib==0.13.2 +scipy==1.3.1 +numpy==1.17.4 +torch==1.3.1 +scikit_learn==0.21.3 diff --git a/utils/__init__.py b/utils/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/utils/data_finetuning.py b/utils/data_finetuning.py new file mode 100644 index 0000000..c7cab28 --- /dev/null +++ b/utils/data_finetuning.py @@ -0,0 +1,66 @@ +import shutil +from pathlib import Path +import pandas as pd +import matplotlib.pyplot as plt +import utils.paths as dirs + +SOURCE_PATH = Path(dirs.DATA) +TARGET_PATH = Path(dirs.DATA, 'shady') +SENSOR_DIM = 800 +PARAMS_DIM = 5 + + +def _create_target_dir(): + if not TARGET_PATH.exists(): + TARGET_PATH.mkdir() + + +def plot_line(filename, line): + if not (SOURCE_PATH / filename).exists(): + print("file doesn't exist") + return + data = pd.read_csv(SOURCE_PATH / filename, header=None).to_numpy() + print('file size:', len(data)) + plt.subplot(2, 1, 1) + plt.plot(data[line, PARAMS_DIM:SENSOR_DIM+PARAMS_DIM]) + plt.subplot(2, 1, 2) + plt.plot(data[line, SENSOR_DIM+PARAMS_DIM:]) + plt.show() + + +def remove_lines(filename, lines): + if not (SOURCE_PATH / filename).exists(): + print("file doesn't exist") + return + data = pd.read_csv(SOURCE_PATH / filename, header=None) + len_before = len(data) + dropped_lines = data.loc[lines, :] + print(len(dropped_lines)) + idx = [i for i in range(len(data)) if i not in lines] + data = data.loc[idx, :] + print(f'dropped {len_before - len(data)} lines') + data.to_csv(SOURCE_PATH / filename, index=False, header=False) + _create_target_dir() + if (TARGET_PATH / filename).exists(): + print('Appended lines to', TARGET_PATH / filename) + part_two = pd.read_csv(TARGET_PATH / filename, header=None).append(dropped_lines) + part_two.to_csv(TARGET_PATH / filename, header=False, index=False) + else: + print('Moved lines to', TARGET_PATH / filename) + dropped_lines.to_csv(TARGET_PATH / filename, header=False, index=False) + + +def move_file(filename): + if not (SOURCE_PATH / filename).exists(): + print("file doesn't exist") + return + _create_target_dir() + if (TARGET_PATH / filename).exists(): + print('Append file to', TARGET_PATH / filename) + part_one = pd.read_csv(SOURCE_PATH / filename, header=None) + part_two = pd.read_csv(TARGET_PATH / filename, header=None).append(part_one) + part_two.to_csv(TARGET_PATH / filename, header=False, index=False) + (SOURCE_PATH / filename).unlink() + else: + print('Moved file to', TARGET_PATH / filename) + shutil.move(SOURCE_PATH / filename, TARGET_PATH / filename) diff --git a/utils/logging_utils.py b/utils/logging_utils.py new file mode 100644 index 0000000..8191222 --- /dev/null +++ b/utils/logging_utils.py @@ -0,0 +1,85 @@ +import logging +import datetime +from collections import namedtuple +from pathlib import Path +import utils.paths as dirs + + +def _get_logger_path(filename, path=None): + if not path: + path = Path(dirs.NN_LOGS) + path = Path(path) + if not path.exists(): + path.mkdir() + if not filename: + filename = 'output_' + str(datetime.datetime.now().strftime('%Y-%m-%d')) + return path / (filename + '.log') + + +def set_logging_parameters(filename, path=None): + logging.basicConfig(level=logging.INFO, + filename=_get_logger_path(filename, path), + format='%(levelname)s:%(name)s: %(message)s') + + +def set_rl_logging_parameters(filename): + set_logging_parameters(filename, dirs.RL_LOGS) + + +def set_multi_file_logging_parameters(logger_name, filename): + log_setup = logging.getLogger(logger_name) + file_handler = logging.FileHandler(_get_logger_path(filename, dirs.RL_LOGS), mode='a') + stream_handler = logging.StreamHandler() + log_setup.setLevel(logging.INFO) + log_setup.addHandler(file_handler) + log_setup.addHandler(stream_handler) + + +def log_to_file(msg, logfile): + log = logging.getLogger(logfile) + log.info(msg) + + +def get_results_from_logs(substrings, path=dirs.RF_LOGS): + files = [f for f in Path(path).iterdir() if f.suffix == '.log' and + all(s in f.name for s in substrings)] + BestValueFiles = namedtuple('BestValueFiles', ['mean', 'std', 'lower_mean_std', + 'upper_mean_std', 'filename']) + results = [] + for f in files: + std = 1.5 + for line in reversed(list(f.open())): + if 'ratio geo_mean' in line: + mean = float(line.split()[-1]) + results.append(BestValueFiles(mean, std, mean / std, mean * std, f.stem)) + break + if 'ratio geo_std' in line: + std = float(line.split()[-1]) + return results + + +def get_best_results(results): + results.sort(key=lambda p: p[0]) + best_mean = results[0] + results.sort(key=lambda p: p[3]) + return best_mean, results[0] + + +def _compare_best_models(shared_ids, unique_ids, path=dirs.RL_LOGS): + print('find results with', shared_ids) + for identifier in unique_ids: + results = get_results_from_logs(shared_ids + identifier, path) + if len(results) == 0: + print('no reults for', identifier) + else: + best_mean, best_mean_std = get_best_results(results) + print('best mean for', identifier, best_mean) + print('best mean*std for', identifier, best_mean_std) + + +if __name__ == '__main__': + # best world models (ff nn and rf) + # _compare_best_models([], [['relu'], ['random_forest']], path=dirs.NN_LOGS) + # best action model (fixed actions, neuroevolution) + _compare_best_models(['new_eval', 'nn', 'target_w=340'], [['fixed_actions'], + ['ga_wo_env']]) diff --git a/utils/nn_utils.py b/utils/nn_utils.py new file mode 100644 index 0000000..c525673 --- /dev/null +++ b/utils/nn_utils.py @@ -0,0 +1,4 @@ + + +def count_trainable_parameters(model): + return sum(p.numel() for p in model.parameters() if p.requires_grad) diff --git a/utils/paths.py b/utils/paths.py new file mode 100644 index 0000000..45044ea --- /dev/null +++ b/utils/paths.py @@ -0,0 +1,21 @@ +from pathlib import Path + + +DATA = '/media/julia/Data/datasets/lufpro/real/adaptive_spreading_preprocessed_data' +NN_WEIGHTS = Path('weights/') +NN_MODELS = Path('models/') +NN_LOGS = 'logs/' +NN_TENSORBOARD = 'runs/' +RF_MODELS = Path('rf_models/') +RF_LOGS = 'logs/' +RL_WEIGHTS = Path('rl_weights/') +RL_MODELS = Path('rl_models/') +RL_LOGS = 'process_control_model/logs/' +RL_PARAMS = 'process_control_model/params/' +RL_TENSORBOARD = 'ga_runs/' +# NN_BACKEND = 'models/300_300_300_relu_lr=0.0001' \ +# '_batch_size=64_split_data_respecting_files_simpleprep3_param=1111.pth' +NN_BACKEND = 'models/100_100_100_relu_lr=0.0001' \ + '_batch_size=64_split_data_respecting_files_new_eval_param=1111.pth' +RF_BACKEND = 'rf_models/est_50_maxdepth_12_min_split_150_max_feats_200_' \ + 'simpleprep3_random_forest_regressor.joblib' diff --git a/utils/tape_width_evaluation.py b/utils/tape_width_evaluation.py new file mode 100644 index 0000000..352b730 --- /dev/null +++ b/utils/tape_width_evaluation.py @@ -0,0 +1,85 @@ +import logging +import torch +import numpy as np +import matplotlib.pyplot as plt +from sklearn.metrics import r2_score +from preprocessing.tape_detection import get_tape_width + + +def geo_mean(stats): + # use log to avoid overflow + # drop zero values (to avoid zero output) TODO other solution? + log_stats = torch.log(stats[stats > 0]) + return torch.exp(torch.sum(log_stats) / log_stats.shape[0]) + + +def geo_standard_dev(stats): + return torch.exp(torch.sqrt( + torch.sum(torch.log(stats[stats > 0] / geo_mean(stats[stats > 0])) ** 2) + / stats[stats > 0].shape[0])) + + +def get_fraction_above_threshold(widths, threshold): + return (torch.sum(widths > threshold).to(dtype=torch.float) / widths.shape[0]).item() + + +def get_fraction_in_area(widths, lower_bounds, step=0.05, title=''): + widths = widths.detach().numpy() + bins = np.empty(len(lower_bounds)) + lower_bounds -= step / 2 + for i, lower_bound in enumerate(lower_bounds): + bins[i] = np.sum((widths >= lower_bound) & (widths < lower_bound + step)) / widths.shape[0] + lower_bounds += step / 2 + hist, bin_edges = np.histogram(widths, bins=12, range=(lower_bounds[0], lower_bounds[-1])) + hist = hist / widths.shape[0] + plt.subplot(1, 2, 1) + plt.bar(["{:.3f}".format(l_b) for l_b in lower_bounds], bins) + plt.subplot(1, 2, 2) + plt.hist(bin_edges[:-1], bin_edges, weights=hist) + plt.suptitle(title) + plt.show() + + +def evaluate_tape_width(targets, preds, plot_stats=False, logger=None): + if not logger: + logger = logging.getLogger(__name__) + targets_widths = get_tape_width( + targets.view(targets.shape[0], targets.shape[-1])) + preds_widths = get_tape_width(preds.view(preds.shape[0], preds.shape[-1])) + logger.info(f'r2 score: {r2_score(targets_widths, preds_widths)}') + print(f'r2 score: {r2_score(targets_widths, preds_widths)}') + dev = torch.abs(targets_widths - preds_widths).to(dtype=torch.float) / targets_widths + width_std, width_mean = torch.std_mean(dev) + logger.info('stats for (target_w - pred_w)/target_w') + logger.info(f'std, mean {width_std.item()} {width_mean.item()}') + logger.info(f'geo_mean dev (> 0) {geo_mean(dev).item()}') + logger.info(f'geo_std dev (> 0) {geo_standard_dev(dev).item()}') + ratio = 1. + torch.abs(1 - preds_widths.to(dtype=torch.float) / + targets_widths.to(dtype=torch.float)) + logger.info('stats for pred_w/target_w') + logger.info(f'ratio mean {torch.mean(ratio).item()}') + logger.info(f'ratio geo_mean {geo_mean(ratio).item()}') + logger.info(f'ratio geo_std {geo_standard_dev(ratio).item()}') + print('stats for pred_w/target_w') + print(f'ratio geo_mean {geo_mean(ratio).item()}') + print(f'ratio geo_std {geo_standard_dev(ratio).item()}') + print(f'average target width={targets_widths.float().mean().item()} ' + f'prediction width={preds_widths.float().mean().item()}') + print(f'std target width={targets_widths.float().std().item()} ' + f'prediction width={preds_widths.float().std().item()}') + logger.info(f'average target width={targets_widths.float().mean().item()} ' + f'prediction width={preds_widths.float().mean().item()}') + logger.info(f'std target width={targets_widths.float().std().item()} ' + f'prediction width={preds_widths.float().std().item()}') + if plot_stats: + thresholds = np.arange(0., 0.3, 0.025) + get_fraction_in_area(dev, thresholds, step=0.025, + title='Distribution (pred_w - target_w) / target_w') + thresholds = np.arange(0.7, 1.3, 0.025) + get_fraction_in_area(ratio, thresholds, step=0.025, title='Distribution pred_w / target_w') + + return geo_mean(ratio).item() + + +def evaluate_tape_width_np_wrapper(targets, preds, plot_stats=False): + return evaluate_tape_width(torch.tensor(targets), torch.tensor(preds), plot_stats)