<jupyter_start><jupyter_text>RegressionLiz Strong, 4/1/2019, modified 4/24/2020<jupyter_code>from __future__ import absolute_import, division, print_function

import pathlib

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import seaborn as sns

import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers

print(tf.__version__)

from IPython.display import set_matplotlib_formats
set_matplotlib_formats('retina')
%matplotlib notebook<jupyter_output><empty_output><jupyter_text>Import data set using pandas  concatenate all of the various datasets to make the test and the training datasets<jupyter_code>def concatenate_dfs(df_old, df_new):
    old_length = len(df_old)
    df_old.index = range(old_length)
    
    new_length = len(df_new)
    df_new.index = range(old_length, old_length+new_length)
    df_old = df_old.append(df_new)
    return df_old<jupyter_output><empty_output><jupyter_text>test dataset<jupyter_code>tdf = pd.read_pickle(r"../oam_sim/random_testset1.pkl")
tdf2 = pd.read_pickle(r"../oam_sim/random_testset2.pkl")
test_dataset = concatenate_dfs(tdf, tdf2)

# Concatenate all of the imported dfs
test_dataset = concatenate_dfs(tdf, tdf2)<jupyter_output><empty_output><jupyter_text>training dataset<jupyter_code>df0 = pd.read_pickle(r"../oam_sim/train_set1.pkl")
df1 = pd.read_pickle(r"../oam_sim/train_set2.pkl")
df = concatenate_dfs(df, test_df)
train_dataset = df.copy()<jupyter_output><empty_output><jupyter_text>Ensure all angles theta are between $\frac{\pi}{2l}\le\theta\le\frac{3\pi}{2l}$<jupyter_code>def symmetric_theta_maker(theta, l=4, lower_wedge_angle=np.pi/(2*4), upper_wedge_angle=3*np.pi/(2*4)):
    """
    Make all the simulations have arcs which are located in the same wedge defined by symmetry. 
    do this by changing the angle (but keeping the radius) of the (R,theta) pair defining the center of the circle
    so that it's appropriately placed within the wedge and not outside of it
    
    Args:
        theta (float): angle between center of light and center of orbit in simulations
        l (float): OAM azimuthal mode number
        lower_wedge_angle (float): lower limit of transformed angles as defined by symmetry
        upper_wedge_angle (float): upper limit of transformed angles as defined by symmetry
        
    Returns:
        theta (float): angle beween center of light and center of orbit in simulations, adjusted to take a value between lower_wdge_angle and upper_wedge_angle
    """
    # if angle is negative, make it positive, theta=-theta
    if theta<0:
        theta *= -1

    # if angle is smaller than lower threshold, add the wedge angular spacing to bring it into wedge
    if theta <= lower_wedge_angle:
        theta += np.pi/l

    # if angle is larger than upper threshold, subtract wedge angular spacing till it's within wedge
    while theta > upper_wedge_angle:
        theta -= np.pi/l

    # cause an error if angle is smaller than lower 
    if theta < lower_wedge_angle:
      raise Exception("angle is smaller than lower permissible wedge angle")
    if theta > upper_wedge_angle:
      raise Exception("angle is larger than upper permissible wedge angle")
    
    return theta
train_dataset['theta'] = train_dataset['theta'].apply(symmetric_theta_maker)
test_dataset['theta'] = test_dataset['theta'].apply(symmetric_theta_maker)<jupyter_output><empty_output><jupyter_text>set the preferences to display all of the columns<jupyter_code>pd.set_option('display.max_columns', None)<jupyter_output><empty_output><jupyter_text>Use only the fits that have less than 20% max error<jupyter_code>dropvals = train_dataset[train_dataset['max_error_normalized']>.2].index
train_dataset = train_dataset.drop(dropvals)

dropvals = test_dataset[test_dataset['max_error_normalized']>.2].index
test_dataset = test_dataset.drop(dropvals)<jupyter_output><empty_output><jupyter_text>Use only counterclockwise==1<jupyter_code>dropvals = test_dataset[test_dataset['clockwise']==1].index
test_dataset = test_dataset.drop(dropvals)
dropvals = train_dataset[train_dataset['clockwise']==1].index
train_dataset = train_dataset.drop(dropvals)
print(len(test_dataset), len(train_dataset))<jupyter_output><empty_output><jupyter_text>make a reference dataset that we can use later to analyze the data<jupyter_code>reference_dataset = test_dataset.copy()<jupyter_output><empty_output><jupyter_text>Delete some data that doesn't go into model<jupyter_code># delete time step 
train_dataset.pop(32);
test_dataset.pop(32);
reference_dataset.pop(32)

# delete max_error & RMSE since we're trying to eliminate absolute amplitude measurements
train_dataset.pop('max_error');
train_dataset.pop('RMSE');
test_dataset.pop('max_error');
test_dataset.pop('RMSE');
reference_dataset.pop('max_error');
reference_dataset.pop('RMSE');

# delete column against which distance is measured 
train_dataset.pop(15);
test_dataset.pop(15);
reference_dataset.pop(15);<jupyter_output><empty_output><jupyter_text>Normalize data<jupyter_code># set up indices.
amp_start = 5 # This starting value is the position of the column called '7'. EG if there is only ang_vel in front of it, then this number should be 1.
amp_end = amp_start + 8
mu_start = amp_end
mu_end = mu_start + 7
sigma_start = mu_end
sigma_end = sigma_start + 8
vert_offset_position = sigma_end
def normalize_data(data, amp_start, amp_end, mu_start, mu_end, sigma_start, sigma_end, vert_offset_position, normalized, delta_t=1e-5, l=4):
    """
    Normalize the data which feeds into the model. 
    Amplitudes normalized by distance between vertical offset and largest amplitude, then scaled by 100.
    Mus normalized by sample rate and a factor of 1/20
    Sigmas normalized by sample rate and a factor of 1/10
    Thetas normalized by pi/(2*l)*30
    D scaled by a factor of 1/10
    R scaled by a factor of 1/10
    ang_vel scaled by a factof of 1/10
    offsetx scaled by a factof of 1/10
    offsety scaled by a factor of of 1/10
    r2 scaled by a factor of 10
    
    Args:
        data (dataframe): data as packaged from multi-Gaussian fitting module
        amp_start (int): index in dataframe where the list of the amplitudes begins
        amp_end (int): index in dataframe where the list of the amplitudes ends
        mu_start (int): index in dataframe where the list of the mu begins
        mu_end (int): index in dataframe where the list of the mu ends
        sigma_start  (int): index in dataframe where the list of the sigma begins
        sigma_end (int): index in dataframe where the list of the sigma ends
        vert_offset_position (int): index in dataframe corresponding to the vertical offset
        normalized (bool): indicates if the data has previously been normalized or not. If true, normalization not repeated.
        delta_t (float): time between samples. sampling rate=1/delta_t
        l (float): OAM azimuthal mode index
    
    Returns:
        normalized (bool): indicates that the data has already been normalized. returns True
        data (dataframe): dataframe with normalized data
    
    """
    data_copy = data.copy()

    if normalized == True:
        print("previously normalized, nothing done.")
        return normalized
    else:
        print("normalizing...")
        ### 1) Treat the amplitudes first
        amplitudes = data.iloc[:, amp_start:amp_end] 
        vert_offset = data.iloc[:, vert_offset_position]
        
        # find the actual height of the peaks, subtracting the vertical offset (baseline) from the max_amps     
        amp_subtr = [amplitudes[a]-vert_offset for a in amplitudes.columns]
        for i in range(8):
            data.iloc[:, amp_start+i] = amp_subtr[i]
            
        # change the vertical offset to 0 since above we've adjusted the amplitudes 
        # such that they're referenced to this vertical offset instead of 0.
        vert_offset -= vert_offset

        # normalize the amplitudes
        max_amplitudes = amplitudes.max(axis=1)
        normalized_amplitudes = amplitudes.div(max_amplitudes, axis=0)*100
        
        data.iloc[:,amp_start:amp_end] = normalized_amplitudes
        
        ### 2) Normalize the mus with the sample rate
        data.iloc[:,mu_start:mu_end] = data.iloc[:,mu_start:mu_end]/delta_t / 20
        
        ### 3) Normalize the sigmas with the sample rate
        data.iloc[:,sigma_start:sigma_end] = data.iloc[:,sigma_start:sigma_end]/delta_t / 10
        
        ### 4) Normalize theta
        data["theta"] = data["theta"] / (np.pi/(2*l))*30
        
        ### 5) Normalize D
        data["D"] = data["D"] / 10
        
        ### 6) Normalize R
        data["orbit_R"] = data["orbit_R"] / 10
        
        ### 7) Make angular velocity smaller
        data["ang_vel"] = data["ang_vel"]/10
        
        ### 8) Make offsets smaller
        data["offsetx"] = data["offsetx"]/10
        data["offsety"] = data["offsety"]/10

        ### 9) Make r2 bigger
        data["r2"] = data["r2"]*10
        
        
        ### 10) Change status of data set normalization from False to True so we don't normalize more than once.
        normalized = True
        
        
        
        
        # if the amplitudes were 0 to start, make them 0 again.
        data = data.where(data_copy != 0, 0.0)
        return normalized, data
# set this to false to normalize. once the data has been normalized, this will take the value True so so it can't be normalized again.
normalized_dataset = False
normalized_testset = False
normalized_referenceset = False
normalized_dataset, train_dataset = normalize_data(train_dataset, amp_start, amp_end, mu_start, mu_end, sigma_start, sigma_end, vert_offset_position, normalized_dataset)
normalized_testset, test_dataset = normalize_data(test_dataset, amp_start, amp_end, mu_start, mu_end, sigma_start, sigma_end, vert_offset_position, normalized_testset)
normalized_refset, reference_dataset = normalize_data(reference_dataset, amp_start, amp_end, mu_start, mu_end, sigma_start, sigma_end, vert_offset_position, normalized_referenceset)
train_stats = train_dataset.describe()
train_stats = train_stats.transpose()
train_stats<jupyter_output><empty_output><jupyter_text>Visualize parameter space<jupyter_code>plt.figure(figsize=[4,4])
for i in train_dataset.orbit_R.unique():
    ds = train_dataset.loc[train_dataset['orbit_R']==i]
    unique_ang_vels = ds.ang_vel.unique()
    plt.plot(i*np.ones(len(unique_ang_vels)), unique_ang_vels, 'k.')
plt.ylabel('angular velocity, $\Omega$ [1/s]')
plt.xlabel('Orbit Radius, R [pixels]')
plt.title('R-$\Omega$ parameter space')

"""
for i in [1050, 1100, 1150]:
    unique_ang_vels = dataset.ang_vel.unique()
    plt.plot(i*np.ones(len(unique_vorticities)), unique_ang_vels, 'rx')
"""<jupyter_output><empty_output><jupyter_text>Get rid of identifying data<jupyter_code># get rid of column with orbit radius--for now--we'll use this to sort later.
#dataset.pop('orbit_R')
train_dataset.pop('offsetx')
train_dataset.pop('offsety')
train_dataset.pop('clockwise')
train_dataset.pop('counterclockwise')
train_dataset.pop('avg_SNR')
#train_dataset.pop("D")
#train_dataset.pop('theta')
train_dataset.pop(31) # vert offset
train_dataset.tail()


# get rid of column with orbit radius--for now--we'll use this to sort later.
#dataset.pop('orbit_R')
test_dataset.pop('offsetx')
test_dataset.pop('offsety')
test_dataset.pop('clockwise')
test_dataset.pop('counterclockwise')
test_dataset.pop('avg_SNR')
#test_dataset.pop("D")
#test_dataset.pop('theta')
test_dataset.pop(31) # vert offset
test_dataset.tail()
# delete data with R<60
dropvals = test_dataset[test_dataset['orbit_R']>60].index
test_dataset = test_dataset.drop(dropvals)
dropvals = train_dataset[train_dataset['orbit_R']>60].index
train_dataset = train_dataset.drop(dropvals)
dropvals = ref_dataset[ref_dataset['orbit_R']>60].index
ref_dataset = ref_dataset.drop(dropvals)<jupyter_output><empty_output><jupyter_text>count the number of unique angular velocities<jupyter_code>print("number unique angular velocities in training dataset", train_dataset['ang_vel'].nunique())

print("number unique angular velocities in test dataset", test_dataset['ang_vel'].nunique())<jupyter_output><empty_output><jupyter_text>clean the data of any unknown values<jupyter_code>test_dataset.isna().sum();
train_dataset.isna().sum();
# drop the rows with na
train_dataset = train_dataset.dropna()
test_dataset = test_dataset.dropna()
train_stats = train_dataset.describe()
#train_stats.pop(1)
#train_stats.pop('ang_vel')
#train_stats.pop('orbit_R')

train_stats = train_stats.transpose()
train_stats<jupyter_output><empty_output><jupyter_text>The model<jupyter_code>def build_model():
    model = keras.Sequential([
        layers.Dense(64, activation=tf.nn.relu, input_shape=[len(train_dataset.keys())]),
        layers.Dense(64, activation=tf.nn.relu),
        layers.Dense(len(train_labels.keys()))
    ])
    optimizer = tf.keras.optimizers.RMSprop(0.001)
    model.compile(loss='mean_squared_error',
                 optimizer=optimizer,
                 metrics=['mean_absolute_error','mean_squared_error','mean_absolute_percentage_error','acc'])
    return model
# Display training progress by printing a single dot for each completed epoch
class PrintDot(keras.callbacks.Callback):
    def on_epoch_end(self, epoch, logs):
        if epoch % 100 == 0: print('')
        print('.',end='')
def plot_history(history):
    hist = pd.DataFrame(history.history)
    hist['epoch'] = history.epoch

    plt.figure()
    plt.xlabel('Epoch')
    plt.ylabel('Mean Abs Error [1/s]')
    plt.plot(hist['epoch'], hist['val_mean_absolute_error'],
           label = 'Val Error')
    plt.plot(hist['epoch'], hist['mean_absolute_error'],
           label='Train Error')

    #plt.ylim([0,5])
    plt.legend()

    plt.figure()
    plt.xlabel('Epoch')
    plt.ylabel('Mean Square Error [$1/s^2$]')

    plt.plot(hist['epoch'], hist['val_mean_squared_error'],
           label = 'Val Error')
    plt.plot(hist['epoch'], hist['mean_squared_error'],
           label='Train Error')
    #plt.ylim([0,20])
    plt.show()
    
    plt.figure()
    plt.plot(hist['epoch'],hist['val_mean_absolute_percentage_error'],label='val mape')
    plt.plot(hist['epoch'],hist['mean_absolute_percentage_error'],label='training mape')

    plt.xlabel('Epoch')
    plt.ylabel('mean absolute percentage error')
    plt.show()


labels_to_sep = ["ang_vel", "orbit_R"]

all_labels_to_eliminate = ["ang_vel", "orbit_R", "D", "theta"]



# separate the features from the labels
train_labels = train_dataset[labels_to_sep]
test_labels = test_dataset[labels_to_sep]


for l in all_labels_to_eliminate:
    train_dataset.pop(l)
    test_dataset.pop(l)

# dont apply more normalization for now
normed_train_data = train_dataset
normed_test_data = test_dataset

# build the model
model = build_model()

# inspect the model
example_batch = normed_train_data[:10]
example_result = model.predict(example_batch)

# train the model

# patience parameter is the number of epochs to check for improvement
early_stop = keras.callbacks.EarlyStopping(monitor='val_loss', patience=30)
EPOCHS = 400
history = model.fit(normed_train_data, train_labels, epochs=EPOCHS,
                    validation_split = 0.2, verbose=0, callbacks=[early_stop, PrintDot()])

plot_history(history)
model.save('model')<jupyter_output><empty_output><jupyter_text>Now we'll see how the model generalizes using the test set<jupyter_code>model = keras.models.load_model('model')
# mae: mean averaged error, mse: mean squared error
loss, mae, mse, mpe, acc = model.evaluate(normed_test_data, test_labels, verbose=0)
print("Testing set Mean Abs Error: {:5.2f} 1/s".format(mae))
print("\nTesting set mean square error: {:5.2f} 1/s^2".format(mse))
print("\nTesting set mean absolute percentage error: {:5.2f}".format(mpe))
# use model to make predictions
test_predictions = model.predict(normed_test_data)
predicted_ang_vels = [test_predictions[i][0] for i in range(len(test_predictions))]
ref_dataset = reference_dataset.copy()

qoi = ref_dataset.iloc[:]['ang_vel'].to_numpy()
x_values, y_values = qoi, predicted_ang_vels
def analyze_means_plot_maker(x_values, y_values, num_bins=30, w=15, x_label='x', y_label=' y'):
    """
    Generate figure to summarize results statistically. Bin results into num_bins and within each, plot a mean and standard deviation.
    
    Args:
        x_values (array):
        y_values (array):
        num_bins (int): the number of segments within which we average.
        w (float): the width of the violin plots
        x_label (string):
        y_label (string): labels for the horizontal and vertical axes of the figure, respectively
        
    Returns:
        -
    """

    min_x, max_x = np.min(x_values), np.max(x_values)+.1 # add this .1 on so that points with this last value are included in the averaging below.
    bin_locations = np.linspace(min_x, max_x, num_bins+1, endpoint=True)

    #find all y values whose positions correspond to the x values of the bin locations
    args_initial = [np.argwhere((x_values>=bl1) & (x_values<bl2)).flatten() for bl1, bl2 in zip(bin_locations[:-1], bin_locations[1:])]
    bin_centers_initial = [np.mean([bl1, bl2]) for bl1, bl2 in zip(bin_locations[:-1], bin_locations[1:])]

    # make sure that if any of these bins made above don't have any elements in them, they get deleted.
    args = [a for a in args_initial if a.size!=0] #get rid of any empty ones
    bin_centers = [b for a,b in zip(args_initial, bin_centers_initial) if a.size!=0]

    predictions = [np.asarray(y_values)[a] for a in args]
    positions = np.array(bin_centers)
    std_devs = [np.std(p) for p in predictions]
    means = [np.mean(p) for p in predictions]
    

    # Create plot
    plt.figure(899, figsize=(10,7))
    plt.plot([0, 1200], [0,  1200], 'k--')
    plt.plot(positions, means, 'k.')
    plt.errorbar(positions, means, yerr=std_devs, ls='none', capsize=3)
    
    plt.xlabel(x_label, fontsize=14)
    plt.ylabel(y_label, fontsize=14)
    

    plt.xticks(fontsize=14)
    plt.yticks(fontsize=14)
    plt.xlim(0, 1200)   

   
    # Create plot
    plt.figure(900, figsize=(10,5))
    plt.plot([0, max_x*1.1], [0, 0], 'k--')
    plt.plot(positions, means-positions, '.')
     
    plt.xlabel(x_label, fontsize=14)
    plt.ylabel(y_label+' means', fontsize=14)
    plt.xlim(0, 1200)
    plt.xticks(fontsize=14)
    plt.yticks(fontsize=14)
    plt.ylim([-200,200])
    
    
    # Create plot
    plt.figure(901, figsize=(10,5))
    plt.plot(positions, np.asarray(std_devs)/np.asarray(positions)*100, '.')
    plt.xlabel(x_label, fontsize=14)
    plt.ylabel('standard deviations as a percentage of position', fontsize=14)
    plt.ylim(0, 100)#np.max(std_devs)*1.5)
    plt.xlim(0, 1200)
    plt.xticks(fontsize=14)
    plt.yticks(fontsize=14)

    
    print("error%:", np.asarray(std_devs)/np.asarray(positions)*100, "average error normalized:", np.mean(np.asarray(std_devs)/np.asarray(positions)*100))
    return
analyze_means_plot_maker(np.asarray(x_values)*10, 10*np.asarray(y_values), num_bins=25, x_label='actual angular velocity', y_label='predicted-actual angular velocity')
plt.figure(899)
plt.ylim(0,1200)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11083986911130733, 'fraction_numerical': 0.017451825689502935, 'mean_word_length': 4.270736381056666, 'pattern_counts': {'":': 0, '<': 74, '<?xml version=': 0, '>': 77, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '421175', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6566, 'n_tokens_neox': 6173, 'n_words': 1619, 'path': 'slvel/regression.ipynb', 'repo_name': 'Liz-Strong/slvel', 'size': '19253'}