text:
import numpy as np
import sklearn.datasets as skdata
import sklearn.metrics as skmetrics
import sklearn.preprocessing as skpreprocess
from sklearn.linear_model import Ridge as RidgeRegression
from matplotlib import pyplot as plt

'''
Name: kanagaraj_kanimozhi

Summary:
Ridge regression is a linear regssion with regualrization term added to it

1. Plotted the graph for scikit-learn RidgeRegrssion model
    a. assigned the x values as log scale of alpha
    b. assigned the y values as score of MSe for MSE plot and scores of R-squared for r-squared plot
    c. assigned the x_limits from 0  to maximum of x values + 1
    d. assigned the y_limits from 0 10 40 for MSE and 0 to 1 for R-squared
2. Implemented Ridge regression model
    a. took the transpose of dataset to get the shape (d,N)
    b. intialized the class
    c. Trained hte model for evry alpha for a dgree of polynomial 2
    d. trained the model on training set using fit function
3. Implemented the fit function
    a. initialized the weights w* = (Z^T.Z + alpha(lambda)*I)^-1.Z^T.y,
        inorder to uniformaly apply lambda(alpha)[scalar] to all weights we multiply with identity matrix
        the inverse gives the solution of w*, if alpha is small the w* is big(more weights used) and if alpha is big w* is small(fewer weights used)
    b. Compute the loss until it statisfies both data fidelity and regularizrion condition
        loss = 1/N((Z.w-y)^T(Z.w-y) + alpha(lambda)w^T.w)
4. Calculate the mean_squared error and r-sqaured scores for training, validation and testing set
5. implented predict, score , mean squeared and r-swaured functions
6. Plotted the MSE and R-sqaure graph for own ridge regession model

Report your scores here. For example,

Results for scikit-learn RidgeRegression model with alpha=1.0
Training set mean squared error: 6.3724
Training set r-squared scores: 0.9267
Validation set mean squared error: 9.6293
Validation set r-squared scores: 0.8626
Testing set mean squared error: 19.2863
Testing set r-squared scores: 0.7531
Results for scikit-learn RidgeRegression model with alpha=10.0
Training set mean squared error: 6.9915
Training set r-squared scores: 0.9195
Validation set mean squared error: 10.5660
Validation set r-squared scores: 0.8493
Testing set mean squared error: 18.0993
Testing set r-squared scores: 0.7683
Results for scikit-learn RidgeRegression model with alpha=100.0
Training set mean squared error: 7.8843
Training set r-squared scores: 0.9093
Validation set mean squared error: 11.9197
Validation set r-squared scores: 0.8300
Testing set mean squared error: 18.5883
Testing set r-squared scores: 0.7620
Results for scikit-learn RidgeRegression model with alpha=1000.0
Training set mean squared error: 8.8610
Training set r-squared scores: 0.8980
Validation set mean squared error: 11.7491
Validation set r-squared scores: 0.8324
Testing set mean squared error: 15.2857
Testing set r-squared scores: 0.8043
Results for scikit-learn RidgeRegression model with alpha=10000.0
Training set mean squared error: 10.0741
Training set r-squared scores: 0.8841
Validation set mean squared error: 11.7167
Validation set r-squared scores: 0.8329
Testing set mean squared error: 13.5444
Testing set r-squared scores: 0.8266
Results for scikit-learn RidgeRegression model with alpha=100000.0
Training set mean squared error: 11.4729
Training set r-squared scores: 0.8680
Validation set mean squared error: 12.5270
Validation set r-squared scores: 0.8213
Testing set mean squared error: 10.8895
Testing set r-squared scores: 0.8606
Results for our RidgeRegression model with alpha=1.0
Training Loss: 6.664
Data Fidelity Loss: 6.413  Regularization Loss: 0.252
Training set mean squared error: 6.4127
Training set r-squared scores: 0.9262
Validation set mean squared error: 8.9723
Validation set r-squared scores: 0.8720
Testing set mean squared error: 18.4835
Testing set r-squared scores: 0.7633
Results for our RidgeRegression model with alpha=10.0
Training Loss: 7.415
Data Fidelity Loss: 7.026  Regularization Loss: 0.389
Training set mean squared error: 7.0258
Training set r-squared scores: 0.9191
Validation set mean squared error: 9.5386
Validation set r-squared scores: 0.8639
Testing set mean squared error: 16.1997
Testing set r-squared scores: 0.7926
Results for our RidgeRegression model with alpha=100.0
Training Loss: 8.347
Data Fidelity Loss: 7.930  Regularization Loss: 0.417
Training set mean squared error: 7.9301
Training set r-squared scores: 0.9087
Validation set mean squared error: 10.6471
Validation set r-squared scores: 0.8481
Testing set mean squared error: 16.3874
Testing set r-squared scores: 0.7902
Results for our RidgeRegression model with alpha=1000.0
Training Loss: 9.429
Data Fidelity Loss: 8.911  Regularization Loss: 0.517
Training set mean squared error: 8.9114
Training set r-squared scores: 0.8974
Validation set mean squared error: 11.2366
Validation set r-squared scores: 0.8397
Testing set mean squared error: 14.5313
Testing set r-squared scores: 0.8139
Results for our RidgeRegression model with alpha=10000.0
Training Loss: 10.707
Data Fidelity Loss: 10.042  Regularization Loss: 0.665
Training set mean squared error: 10.0420
Training set r-squared scores: 0.8844
Validation set mean squared error: 11.8909
Validation set r-squared scores: 0.8304
Testing set mean squared error: 13.8512
Testing set r-squared scores: 0.8226
Results for our RidgeRegression model with alpha=100000.0
Training Loss: 12.984
Data Fidelity Loss: 11.598  Regularization Loss: 1.385
Training set mean squared error: 11.5984
Training set r-squared scores: 0.8665
Validation set mean squared error: 13.1313
Validation set r-squared scores: 0.8127
Testing set mean squared error: 11.8234
Testing set r-squared scores: 0.8486
'''

'''
Implementation of ridge regression
'''
class RidgeRegressionClosedForm(object):

    def __init__(self,alpha):

        # Define private variables
        self.__weights = None
        self.__alpha = alpha

    def fit(self, z, y, alpha =0.0):
        '''
        Fits the model to x and y using closed form solution

        Args:
            z : numpy
                d x N feature vector
            y : numpy
                1 x N ground-truth label
            alpha : float
                weight (lambda) of regularization term
        '''

        # TODO: Implement the fit function
        #computing the weights
        # convert from (d,N) to (N,d) by taking its transpose
        Z = z.T
        #print(z.shape,Z.shape,y.shape)
        I = np.identity(Z.shape[1])
        N = np.size(y) #no of features
        #print(I.shape)
        #w* = (Z^T.Z + alpha * I)^-1.Z^T.y
        inverse = np.linalg.inv(np.matmul(Z.T,Z) + self.__alpha * I) # inverse of (Z^T.Z + alpha*I)
        second_term = np.matmul(Z.T,y)
        self.__weights = np.matmul(inverse, second_term)

        # TODO: Compute loss
        # loss = 1/N((Z.w-y)^T(Z.w-y) + alpha(lambda)w^T.w)
        #wehre , loss of data fidelity = 1/N(Z.w-y)^T(Z.w-y)
        #loss of regularization = 1/N alpha(lambda)W^T.w
        loss_data_fidelity = 1/N * (np.matmul(np.transpose(np.matmul(Z, self.__weights)- y),(np.matmul(Z, self.__weights)- y)))
        loss_regularization = 1/N * (np.matmul((self.__alpha * self.__weights.T), self.__weights))

        loss = loss_data_fidelity + loss_regularization

        print('Training Loss: {:.3f}'.format(loss))
        print('Data Fidelity Loss: {:.3f}  Regularization Loss: {:.3f}'.format(
            loss_data_fidelity, loss_regularization))

    def predict(self, z):
        '''
        Predicts the label for each feature vector x

        Args:
            z : numpy
                d x N feature vector

        Returns:
            numpy : d x 1 label vector
        '''

        # TODO: Implement the predict function

        predictions = np.matmul(self.__weights.T,z)

        return predictions

    def __score_r_squared(self, y_hat, y):
        '''
        Measures the r-squared score from groundtruth y

        Args:
            y_hat : numpy
                1 x N predictions
            y : numpy
                1 x N ground-truth label

        Returns:
            float : r-squared score
        '''

        # TODO: Implement the __score_r_squared function
        # u : sum(y_hat - y)^2
        sum_squared_errors = np.sum((y_hat - y) ** 2)

        #v total variance : sum(y - mean(y))^2
        sum_variance = np.sum((y - np.mean(y)) ** 2)

        return 1.0 - (sum_squared_errors / sum_variance)

    def __score_mean_squared_error(self, y_hat, y):
        '''
        Measures the mean squared error (distance) from groundtruth y

        Args:
            y_hat : numpy
                1 x N predictions
            y : numpy
                1 x N ground-truth label

        Returns:
            float : mean squared error (mse)
        '''

        # TODO: Implement the __score_mean_squared_error function

        return np.mean((y_hat - y) ** 2)

    def score(self, x, y, scoring_func='r_squared'):
        '''
        Predicts real values from x and measures the mean squared error (distance)
        or r-squared from groundtruth y

        Args:
            x : numpy
                d x N feature vector
            y : numpy
                1 x N ground-truth label
            scoring_func : str
                scoring function: r_squared, mean_squared_error

        Returns:
            float : mean squared error (mse)
        '''

        # TODO: Implement the score function
        predictions = self.predict(x)
        if scoring_func == 'r_squared':
            return self.__score_r_squared(predictions, y)
        elif scoring_func == 'mean_squared_error':
            return self.__score_mean_squared_error(predictions, y)
        else:
            raise ValueError('Encountered unsupported scoring_func: {}'.format(scoring_func))

'''
Utility functions to compute error and plot
'''
def score_mean_squared_error(model, x, y):
    '''
    Scores the model on mean squared error metric

    Args:
        model : object
            trained model, assumes predict function returns N x d predictions
        x : numpy
            N x d numpy array of features
        y : numpy
            N x 1 groundtruth vector
    Returns:
        float : mean squared error
    '''

    # Implement the score mean squared error function
    predictions = model.predict(x)
    mse = skmetrics.mean_squared_error(predictions, y)
    return mse

def plot_results(axis,
                 x_values,
                 y_values,
                 labels,
                 colors,
                 x_limits,
                 y_limits,
                 x_label,
                 y_label):
    '''
    Plots x and y values using line plot with labels and colors

    Args:
        axis :  pyplot.ax
            matplotlib subplot axis
        x_values : list[numpy]
            list of numpy array of x values
        y_values : list[numpy]
            list of numpy array of y values
        labels : str
            list of names for legend
        colors : str
            colors for each line
        x_limits : list[float]
            min and max values of x axis
        y_limits : list[float]
            min and max values of y axis
        x_label : list[float]
            name of x axis
        y_label : list[float]
            name of y axis
    '''

    # Iterate through x_values, y_values, labels, and colors and plot them
    # with associated legend
    for x, y, label, color in zip(x_values, y_values, labels, colors):
        axis.plot(x, y, marker='o', color=color, label=label)
        axis.legend(loc='best')

    # Set x and y limits
    axis.set_xlim(x_limits)
    axis.set_ylim(y_limits)

    # Set x and y labels
    axis.set_xlabel(x_label)
    axis.set_ylabel(y_label)

if __name__ == '__main__':

    boston_housing_data = skdata.load_boston()
    x = boston_housing_data.data
    y = boston_housing_data.target

    # 80 percent train, 10 percent validation, 10 percent test split
    train_idx = []
    val_idx = []
    test_idx = []
    for idx in range(x.shape[0]):
        if idx and idx % 10 == 9:
            val_idx.append(idx)
        elif idx and idx % 10 == 0:
            test_idx.append(idx)
        else:
            train_idx.append(idx)

    x_train, x_val, x_test = x[train_idx, :], x[val_idx, :], x[test_idx, :]
    y_train, y_val, y_test = y[train_idx], y[val_idx], y[test_idx]

    '''
    Trains and tests Ridge regression model from scikit-learn
    '''
    # Initialize polynomial expansion of degree 2
    poly_transform = skpreprocess.PolynomialFeatures(degree=2)

    # Compute the polynomial terms needed for the data
    poly_transform.fit(x_train)

    # Transform the data by nonlinear mapping
    x_poly_train = poly_transform.transform(x_train)
    x_poly_val = poly_transform.transform(x_val)
    x_poly_test = poly_transform.transform(x_test)

    # Initialize empty lists to store scores for MSE and R-squared
    scores_mse_ridge_scikit_train = []
    scores_r2_ridge_scikit_train = []
    scores_mse_ridge_scikit_val = []
    scores_r2_ridge_scikit_val = []
    scores_mse_ridge_scikit_test = []
    scores_r2_ridge_scikit_test = []

    alphas = [1.0, 10.0, 100.0, 1000.0, 10000.0, 100000.0]

    for alpha in alphas:

        # Initialize scikit-learn ridge regression model
        model_ridge_scikit = RidgeRegression(alpha=alpha)

        # Trains scikit-learn ridge regression model
        model_ridge_scikit.fit(x_poly_train, y_train)

        print('Results for scikit-learn RidgeRegression model with alpha={}'.format(alpha))

        # Test model on training set
        score_mse_ridge_scikit_train = score_mean_squared_error(model_ridge_scikit, x_poly_train, y_train)
        print('Training set mean squared error: {:.4f}'.format(score_mse_ridge_scikit_train))

        score_r2_ridge_scikit_train = model_ridge_scikit.score(x_poly_train, y_train)
        print('Training set r-squared scores: {:.4f}'.format(score_r2_ridge_scikit_train))

        # Save MSE and R-squared training scores
        scores_mse_ridge_scikit_train.append(score_mse_ridge_scikit_train)
        scores_r2_ridge_scikit_train.append(score_r2_ridge_scikit_train)

        # Test model on validation set
        score_mse_ridge_scikit_val = score_mean_squared_error(model_ridge_scikit, x_poly_val, y_val)
        print('Validation set mean squared error: {:.4f}'.format(score_mse_ridge_scikit_val))

        score_r2_ridge_scikit_val = model_ridge_scikit.score(x_poly_val, y_val)
        print('Validation set r-squared scores: {:.4f}'.format(score_r2_ridge_scikit_val))

        # Save MSE and R-squared validation scores
        scores_mse_ridge_scikit_val.append(score_mse_ridge_scikit_val)
        scores_r2_ridge_scikit_val.append(score_r2_ridge_scikit_val)

        # Test model on testing set
        score_mse_ridge_scikit_test = score_mean_squared_error(model_ridge_scikit, x_poly_test, y_test)
        print('Testing set mean squared error: {:.4f}'.format(score_mse_ridge_scikit_test))

        score_r2_ridge_scikit_test = model_ridge_scikit.score(x_poly_test, y_test)
        print('Testing set r-squared scores: {:.4f}'.format(score_r2_ridge_scikit_test))

        # Save MSE and R-squared testing scores
        scores_mse_ridge_scikit_test.append(score_mse_ridge_scikit_test)
        scores_r2_ridge_scikit_test.append(score_r2_ridge_scikit_test)

    # Convert each scores to NumPy arrays
    scores_mse_ridge_scikit_train = np.array(scores_mse_ridge_scikit_train)
    scores_mse_ridge_scikit_val = np.array(scores_mse_ridge_scikit_val)
    scores_mse_ridge_scikit_test = np.array(scores_mse_ridge_scikit_test)
    scores_r2_ridge_scikit_train = np.array(scores_r2_ridge_scikit_train)
    scores_r2_ridge_scikit_val = np.array(scores_r2_ridge_scikit_val)
    scores_r2_ridge_scikit_test = np.array(scores_r2_ridge_scikit_test)

    # Clip each set of MSE scores between 0 and 40
    scores_mse_ridge_scikit_train = np.clip(scores_mse_ridge_scikit_train, 0.0, 40.0)
    scores_mse_ridge_scikit_val = np.clip(scores_mse_ridge_scikit_val, 0.0, 40.0)
    scores_mse_ridge_scikit_test = np.clip(scores_mse_ridge_scikit_test, 0.0, 40.0)

    # Clip each set of R-squared scores between 0 and 1
    scores_r2_ridge_scikit_train = np.clip(scores_r2_ridge_scikit_train, 0.0, 1.0)
    scores_r2_ridge_scikit_val = np.clip(scores_r2_ridge_scikit_val, 0.0, 1.0)
    scores_r2_ridge_scikit_test = np.clip(scores_r2_ridge_scikit_test, 0.0, 1.0)

    # Create figure for training, validation and testing scores for different features
    n_experiments = scores_mse_ridge_scikit_train.shape[0]
    fig = plt.figure()

    labels = ['Training', 'Validation', 'Testing']
    colors = ['blue', 'red', 'green']

    # Create the first subplot of a 1 by 2 figure to plot MSE for training, validation, testing
    ax = fig.add_subplot(1, 2, 1)

    # TODO: Set x (alpha in log scale) and y values (MSE)
    x_values = [np.log(np.asarray(alphas) + 1.0)] * n_experiments
    y_values = [
        scores_mse_ridge_scikit_train,
        scores_mse_ridge_scikit_val,
        scores_mse_ridge_scikit_test
    ]

    # TODO: Plot MSE scores for training, validation, testing sets
    # Set x limits to 0 to max of x_values + 1 and y limits between 0 and 40
    # Set x label to 'alpha (log scale)' and y label to 'MSE',
    plot_results(
        axis=ax,
        x_values=x_values,
        y_values=y_values,
        labels=labels,
        colors=colors,
        x_limits=[0.0,np.max(x_values) +1],
        y_limits=[0.0, 40.0],
        x_label='alpha',
        y_label='MSE')

    # Create the second subplot of a 1 by 2 figure to plot R-squared for training, validation, testing
    ax = fig.add_subplot(1, 2, 2)

    # TODO: Set x (alpha in log scale) and y values (R-squared)
    x_values = [np.log(np.asarray(alphas) + 1.0)] * n_experiments
    y_values = [
        scores_r2_ridge_scikit_train,
        scores_r2_ridge_scikit_val,
        scores_r2_ridge_scikit_test
    ]

    # TODO: Plot R-squared scores for training, validation, testing sets
    # Set x limits to 0 to max of x_values + 1 and y limits between 0 and 1
    # Set x label to 'alpha (log scale)' and y label to 'R-squared',
    plot_results(
        axis=ax,
        x_values=x_values,
        y_values=y_values,
        labels=labels,
        colors=colors,
        x_limits=[0.0,np.max(x_values) +1],
        y_limits=[0.0, 1.0],
        x_label='alpha',
        y_label='R-squared')
    # TODO: Create super title 'Scikit-Learn Ridge Regression on Training, Validation and Testing Sets'
    plt.suptitle('Scikit-Learn Ridge Regression on Training, Validation and Testing Sets')

    '''
    Trains and tests our ridge regression model using different alphas
    '''

    # Initialize empty lists to store scores for MSE and R-squared
    scores_mse_ridge_ours_train = []
    scores_r2_ridge_ours_train = []
    scores_mse_ridge_ours_val = []
    scores_r2_ridge_ours_val = []
    scores_mse_ridge_ours_test = []
    scores_r2_ridge_ours_test = []

    # TODO: convert dataset (N x d) to correct shape (d x N)
    x_poly_train = np.transpose(x_poly_train, axes=(1, 0))
    x_poly_val = np.transpose(x_poly_val, axes=(1, 0))
    x_poly_test = np.transpose(x_poly_test, axes=(1, 0))

    # For each alpha, train a ridge regression model on degree 2 polynomial features
    for alpha in alphas:

        # TODO: Initialize our own ridge regression model
        model_our_ridge = RidgeRegressionClosedForm(alpha=alpha)
        print('Results for our RidgeRegression model with alpha={}'.format(alpha))

        # TODO: Train model on training set
        model_our_ridge.fit(x_poly_train, y_train)

        # TODO: Test model on training set using mean squared error and r-squared

        score_mse_ridge_ours_train = score_mean_squared_error(model_our_ridge, x_poly_train, y_train)
        print('Training set mean squared error: {:.4f}'.format(score_mse_ridge_ours_train))

        score_r2_ridge_ours_train = model_our_ridge.score(x_poly_train, y_train)
        print('Training set r-squared scores: {:.4f}'.format(score_r2_ridge_ours_train))

        # TODO: Save MSE and R-squared training scores
        scores_mse_ridge_ours_train.append(score_mse_ridge_ours_train)
        scores_r2_ridge_ours_train.append(score_r2_ridge_ours_train)

        # TODO: Test model on validation set using mean squared error and r-squared
        score_mse_ridge_ours_val = score_mean_squared_error(model_our_ridge, x_poly_val, y_val)
        print('Validation set mean squared error: {:.4f}'.format(score_mse_ridge_ours_val))

        score_r2_ridge_ours_val = model_our_ridge.score(x_poly_val, y_val)
        print('Validation set r-squared scores: {:.4f}'.format(score_r2_ridge_ours_val))

        # TODO: Save MSE and R-squared validation scores
        scores_mse_ridge_ours_val.append(score_mse_ridge_ours_val)
        scores_r2_ridge_ours_val.append(score_r2_ridge_ours_val)

        # TODO: Test model on testing set using mean squared error and r-squared
        score_mse_ridge_ours_test = score_mean_squared_error(model_our_ridge, x_poly_test, y_test)
        print('Testing set mean squared error: {:.4f}'.format(score_mse_ridge_ours_test))

        score_r2_ridge_ours_test = model_our_ridge.score(x_poly_test, y_test)
        print('Testing set r-squared scores: {:.4f}'.format(score_r2_ridge_ours_test))

        # TODO: Save MSE and R-squared testing scores
        scores_mse_ridge_ours_test.append(score_mse_ridge_ours_test)
        scores_r2_ridge_ours_test.append(score_r2_ridge_ours_test)

    # Convert each scores to NumPy arrays
    scores_mse_ridge_ours_train = np.array(scores_mse_ridge_ours_train)
    scores_mse_ridge_ours_val = np.array(scores_mse_ridge_ours_val)
    scores_mse_ridge_ours_test = np.array(scores_mse_ridge_ours_test)
    scores_r2_ridge_ours_train = np.array(scores_r2_ridge_ours_train)
    scores_r2_ridge_ours_val = np.array(scores_r2_ridge_ours_val)
    scores_r2_ridge_ours_test = np.array(scores_r2_ridge_ours_test)

    # TODO: Clip each set of MSE scores between 0 and 40
    scores_mse_ridge_ours_train = np.clip(scores_mse_ridge_ours_train, 0.0, 40.0)
    scores_mse_ridge_ours_val = np.clip(scores_mse_ridge_ours_val, 0.0, 40.0)
    scores_mse_ridge_ours_test = np.clip(scores_mse_ridge_ours_test, 0.0, 40.0)

    # TODO: Clip each set of R-squared scores between 0 and 1
    scores_r2_ridge_ours_train = np.clip(scores_r2_ridge_ours_train, 0.0, 1.0)
    scores_r2_ridge_ours_val = np.clip(scores_r2_ridge_ours_val, 0.0, 1.0)
    scores_r2_ridge_ours_test = np.clip(scores_r2_ridge_ours_test, 0.0, 1.0)

    # Create figure for training, validation and testing scores for different features
    n_experiments = scores_mse_ridge_ours_train.shape[0]
    fig = plt.figure()

    labels = ['Training', 'Validation', 'Testing']
    colors = ['blue', 'red', 'green']

    # Create the first subplot of a 1 by 2 figure to plot MSE for training, validation, testing
    ax = fig.add_subplot(1, 2, 1)

    # TODO: Set x (alpha in log scale) and y values (MSE)
    x_values = [np.log(np.asarray(alphas) + 1.0)] * n_experiments
    y_values = [
        scores_mse_ridge_ours_train,
        scores_mse_ridge_ours_val,
        scores_mse_ridge_ours_test
    ]

    # TODO: Plot MSE scores for training, validation, testing sets
    # Set x limits to 0 to max of x_values + 1 and y limits between 0 and 40
    # Set x label to 'alpha (log scale)' and y label to 'MSE',
    plot_results(
        axis=ax,
        x_values=x_values,
        y_values=y_values,
        labels=labels,
        colors=colors,
        x_limits=[0.0, np.max(x_values)+ 1],
        y_limits=[0.0, 40.0],
        x_label='alpha (log scale)',
        y_label='MSE'
    )

    # Create the second subplot of a 1 by 2 figure to plot R-squared for training, validation, testing
    ax = fig.add_subplot(1, 2, 2)

    # TODO: Set x (alpha in log scale) and y values (R-squared)
    x_values = [np.log(np.asarray(alphas) + 1.0)] * n_experiments
    y_values = [
        scores_r2_ridge_ours_train,
        scores_r2_ridge_ours_val,
        scores_r2_ridge_ours_test
    ]

    # TODO: Plot R-squared scores for training, validation, testing sets
    # Set x limits to 0 to max of x_values + 1 and y limits between 0 and 1
    # Set x label to 'alpha (log scale)' and y label to 'R-squared',
    plot_results(
        axis=ax,
        x_values=x_values,
        y_values=y_values,
        labels=labels,
        colors=colors,
        x_limits=[0.0, np.max(x_values)+ 1],
        y_limits=[0.0, 1.0],
        x_label='alpha (log scale)',
        y_label='R-squared'
    )

    # TODO: Create super title 'Our Ridge Regression on Training, Validation and Testing Sets'
    plt.suptitle('Our Ridge Regression on Training, Validation and Testing Sets')
    plt.show()


text:
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

from matplotlib.patches import Ellipse
from sklearn.datasets import make_blobs
from scipy.stats import multivariate_normal

sns.set()

# Set random seed for reproducibility
np.random.seed(1000)

nb_samples = 250
nb_unlabeled = 200
threshold = 1e-4

# First Gaussian
m1 = np.array([-2.0, -2.5])
c1 = np.array([[1.0, 1.0],
               [1.0, 2.0]])
q1 = 0.5

# Second Gaussian
m2 = np.array([1.0, 3.0])
c2 = np.array([[2.0, -1.0],
               [-1.0, 3.5]])
q2 = 0.5

m1_old = np.zeros((2,))
c1_old = np.zeros((2, 2))
q1_old = 0

m2_old = np.zeros((2,))
c2_old = np.zeros((2, 2))
q2_old = 0

def total_norm():
    global m1, m1_old, m2, m2_old, c1, c1_old, c2, c2_old, q1, q1_old, q2, q2_old
    return np.linalg.norm(m1 - m1_old) + \
           np.linalg.norm(m2 - m2_old) + \
           np.linalg.norm(c1 - c1_old) + \
           np.linalg.norm(c2 - c2_old) + \
           np.linalg.norm(q1 - q1_old) + \
           np.linalg.norm(q2 - q2_old)

def show_dataset():
    w1, v1 = np.linalg.eigh(c1)
    w2, v2 = np.linalg.eigh(c2)

    nv1 = v1 / np.linalg.norm(v1)
    nv2 = v2 / np.linalg.norm(v2)

    a1 = np.arccos(np.dot(nv1[:, 1], [1.0, 0.0]) / np.linalg.norm(nv1[:, 1])) * 180.0 / np.pi
    a2 = np.arccos(np.dot(nv2[:, 1], [1.0, 0.0]) / np.linalg.norm(nv2[:, 1])) * 180.0 / np.pi

    fig, ax = plt.subplots(figsize=(20, 15))

    ax.scatter(X[Y == 0, 0], X[Y == 0, 1], s=120, marker='o', label='Class 1')
    ax.scatter(X[Y == 1, 0], X[Y == 1, 1], s=120, marker='d', label='Class 2')
    ax.scatter(X[Y == -1, 0], X[Y == -1, 1], s=120, marker='x', label='Unlabeled')

    g1 = Ellipse(xy=m1, width=w1[1] * 3, height=w1[0] * 3, fill=False, linestyle='dashed', angle=a1, color='black',
                 linewidth=1)
    g1_1 = Ellipse(xy=m1, width=w1[1] * 2, height=w1[0] * 2, fill=False, linestyle='dashed', angle=a1, color='black',
                   linewidth=2)
    g1_2 = Ellipse(xy=m1, width=w1[1] * 1.4, height=w1[0] * 1.4, fill=False, linestyle='dashed', angle=a1,
                   color='black', linewidth=3)

    g2 = Ellipse(xy=m2, width=w2[1] * 3, height=w2[0] * 3, fill=False, linestyle='dashed', angle=a2, color='black',
                 linewidth=1)
    g2_1 = Ellipse(xy=m2, width=w2[1] * 2, height=w2[0] * 2, fill=False, linestyle='dashed', angle=a2, color='black',
                   linewidth=2)
    g2_2 = Ellipse(xy=m2, width=w2[1] * 1.4, height=w2[0] * 1.4, fill=False, linestyle='dashed', angle=a2,
                   color='black', linewidth=3)

    ax.set_xlabel(r'$x_0$', fontsize=18)
    ax.set_ylabel(r'$x_1$', fontsize=18)

    ax.add_artist(g1)
    ax.add_artist(g1_1)
    ax.add_artist(g1_2)
    ax.add_artist(g2)
    ax.add_artist(g2_1)
    ax.add_artist(g2_2)

    ax.legend(fontsize=18)

    plt.show()

if __name__ == '__main__':
    # Generate dataset
    X, Y = make_blobs(n_samples=nb_samples, n_features=2, centers=2, cluster_std=1.25, random_state=100)

    unlabeled_idx = np.random.choice(np.arange(0, nb_samples, 1), replace=False, size=nb_unlabeled)
    Y[unlabeled_idx] = -1

    # Show the dataset with the initial Gaussians
    show_dataset()

    # Training process
    while total_norm() > threshold:
        m1_old = m1.copy()
        c1_old = c1.copy()
        q1_old = q1

        m2_old = m2.copy()
        c2_old = c2.copy()
        q2_old = q2

        Pij = np.zeros((nb_samples, 2))

        # E Step
        for i in range(nb_samples):
            if Y[i] == -1:
                p1 = multivariate_normal.pdf(X[i], m1, c1, allow_singular=True) * q1
                p2 = multivariate_normal.pdf(X[i], m2, c2, allow_singular=True) * q2
                Pij[i] = [p1, p2] / (p1 + p2)

            else:
                Pij[i, :] = [1.0, 0.0] if Y[i] == 0 else [0.0, 1.0]

        # M Step
        n = np.sum(Pij, axis=0)
        m = np.sum(np.dot(Pij.T, X), axis=0)

        m1 = np.dot(Pij[:, 0], X) / n[0]
        m2 = np.dot(Pij[:, 1], X) / n[1]

        q1 = n[0] / float(nb_samples)
        q2 = n[1] / float(nb_samples)

        c1 = np.zeros((2, 2))
        c2 = np.zeros((2, 2))

        for t in range(nb_samples):
            c1 += Pij[t, 0] * np.outer(X[t] - m1, X[t] - m1)
            c2 += Pij[t, 1] * np.outer(X[t] - m2, X[t] - m2)

        c1 /= n[0]
        c2 /= n[1]

    # Show the final Gaussians
    show_dataset()

    # Check some points
    print('The first 10 unlabeled samples:')
    print(np.round(X[Y == -1][0:10], 3))

    print('\nCorresponding Gaussian assigments:')
    print(np.round(Pij[Y == -1][0:10], 3))

ext: py

reason: import numpy,from scipy


text:
import os.path
import numpy as np

home = os.path.expanduser("~")
data_dir = os.path.join(home, "sdtw_data")
ucr_dir = os.path.join(data_dir, "UCR_TS_Archive_2015")

def _parse_ucr(filename):
    y = []
    X = []
    for line in open(filename):
        line = line.strip()
        arr = line.split(",")
        label = int(arr[0])
        feat = list(map(float, arr[1:]))
        feat = np.array(feat).reshape(-1, 1)
        y.append(label)
        X.append(feat)
    return X, np.array(y)

def list_ucr():
    return sorted(os.listdir(ucr_dir))

def load_ucr(name):
    folder = os.path.join(ucr_dir, name)
    tr = os.path.join(folder, "%s_TRAIN" % name)
    te = os.path.join(folder, "%s_TEST" % name)

    try:
        X_tr, y_tr = _parse_ucr(tr)
        X_te, y_te = _parse_ucr(te)
    except IOError:
        raise IOError("Please copy UCR_TS_Archive_2015/ to $HOME/sdtw_data/. "
                      "Download from www.cs.ucr.edu/~eamonn/time_series_data.")

    y_tr = np.array(y_tr)
    y_te = np.array(y_te)
    X_tr = np.array(X_tr)
    X_te = np.array(X_te)

    return X_tr, y_tr, X_te, y_te

max_stars_repo_name: sophiaas/alphacsc

max_issues_repo_name: sophiaas/alphacsc

max_forks_repo_name: sophiaas/alphacsc

reason: import numpy


text:
import numpy as np
import scipy.sparse
import smurff
import pytest

verbose = 0
seed = 1234

# 4 different types of side info
def no_side_info(U):
    return None

def sparse_side_info(U):
    return smurff.make_sparse(U, 0.5, seed=seed)

def binary_side_info(U):
    F = np.digitize(U, bins = [.0])
    F = scipy.sparse.coo_matrix(F)
    return F

def dense_side_info(U):
    return U

# 5 different noise configs
def noise_fixed5():
    return smurff.FixedNoise(5.0)

def noise_fixed10():
    return smurff.FixedNoise(10.0)

def noise_adaptive1():
    return smurff.AdaptiveNoise(1.0, 10)

def noise_adaptive10():
    return smurff.AdaptiveNoise(10.0, 100.0)

def noise_probit():
    return smurff.ProbitNoise(.0)

def train_test(density, nmodes, side_info):
    np.random.seed(seed)
    Us = [ np.random.randn(i*4, 1) for i in range(1,nmodes+1) ]
    subscripts = [ [i+1, 0] for i in range(nmodes) ]
    Y = np.einsum(*[j for i in zip(Us,subscripts) for j in i])
    if density < 1.:
        _, Y = smurff.make_train_test(Y, density, seed=seed)
    Ytrain, Ytest = smurff.make_train_test(Y, 0.5, seed=seed)
    return Ytrain, Ytest, side_info(Us[0])

@pytest.mark.parametrize('density', [1.0, 0.5])
@pytest.mark.parametrize('nmodes', [2, 3, 4])
@pytest.mark.parametrize('side_info', [no_side_info, sparse_side_info, binary_side_info, dense_side_info])
@pytest.mark.parametrize('noise_model', [noise_probit, noise_fixed5, noise_fixed10, noise_adaptive10, noise_adaptive1])
def test_noise_model(density, nmodes, side_info, noise_model):
    Ytrain, Ytest, si = train_test(density, nmodes, side_info)
    nm = noise_model()

    priors = ['normal'] * nmodes
    if si is not None:
        priors[0] = 'macau'

    trainSession = smurff.TrainSession(priors = priors, num_latent=8, burnin=20, nsamples=20, threshold=.0, seed=seed, verbose=verbose)

    trainSession.addTrainAndTest(Ytrain, Ytest, nm)
    if not si is None:
        trainSession.addSideInfo(0, si, smurff.SampledNoise(1.), direct=True)

    trainSession.init()
    while trainSession.step():
        pass

    predictions = trainSession.getTestPredictions()
    assert Ytest.nnz == len(predictions)
    if isinstance(nm, smurff.ProbitNoise):
        assert trainSession.getStatus().auc_avg <= 1.
        assert trainSession.getStatus().auc_avg >= 0.
    else:
        assert trainSession.getRmseAvg() < 10.
    return predictions

max_stars_repo_name: ExaScience/smurff

max_issues_repo_name: ExaScience/smurff

max_forks_repo_name: ExaScience/smurff

reason: import numpy,import scipy


text:
"""
python3 test.py
"""

######################
# Import main modules
######################

import unittest
import argparse
import logging
import time

import librosa
import sounddevice as sd

import numpy as np
import tensorflow as tf
import io
import csv

from utils import *
from model_interface import YAMNetInterface

######################
# Globals
######################

SAMPLERATE_HZ = 16000 # requirement of YAMNet

######################
# Unit test
######################

class MyTestCase(unittest.TestCase):
    def test_yamnet_sample(self):
        ''' Test a single prediction '''
        print("Single prediction test")

        # Get parse args
        args = get_parser()

        # Declare model interface
        model_interface = YAMNetInterface(args)

        # Get sample with open_data_source
        y, _ = librosa.load("sample_data/street_music_sample.wav",SAMPLERATE_HZ)

        # Make a prediction
        yh_k, yh_conf = model_interface.predict(y)

        for i in range(args.TOP_K):
            print("rank."+str(i+1)+".class", yh_k[i])
            print("rank."+str(i+1)+".prob", yh_conf[i])

        # Check outputs
        self.assertEqual(yh_k[0],'Music')
        self.assertEqual(yh_k[1],'Hip hop music')
        self.assertEqual(yh_k[2],'Electronic music')

    def test_local_mode_a(self):
        ''' Test for mode a done with local mic '''
        print("Local live stream test of mode a")

        # Get parse args
        args = get_parser()

        # Declare model interface
        model_interface = YAMNetInterface(args)

        def predict_sound(indata,frames,time,status):
            indata = np.squeeze(indata,axis=-1)
            yh_k, yh_conf = model_interface.predict(indata)
            print(yh_k[0], yh_conf[0])
        start_time = time.time()
        with sd.InputStream(samplerate=SAMPLERATE_HZ,
                            dtype='float32',
                            blocksize=SAMPLERATE_HZ*2,
                            callback=predict_sound):
            sd.sleep(10000)
        print("--- %s seconds ---" % (time.time() - start_time))

    def test_local_mode_b(self):
        ''' Test for mode b done with a local mic '''
        print("Local live stream mode b")

        # Get parse args
        args = get_parser()
        args.DURATION_S = 5
        args.MODE = 'b'
        args.WATCH_SOUNDS = ['Music','Ice cream truck, ice cream van']
        args.TOP_K = 30

        # Declare model interface
        model_interface = YAMNetInterface(args)

        def predict_sound(indata,frames,time,status):
            indata = np.squeeze(indata,axis=-1)
            yh_k, yh_conf = model_interface.predict(indata)
            matched_sounds = list(set(args.WATCH_SOUNDS).intersection(set(yh_k)))
            print(matched_sounds)

        start_time = time.time()
        with sd.InputStream(samplerate=SAMPLERATE_HZ,
                            dtype='float32',
                            blocksize=SAMPLERATE_HZ*2,
                            callback=predict_sound):
            sd.sleep(30000)
        print("--- %s seconds ---" % (time.time() - start_time))

######################

if __name__ == "__main__":
    unittest.main()

max_stars_repo_name: gemblerz/yamnet-plugin

max_issues_repo_name: gemblerz/yamnet-plugin

max_forks_repo_name: gemblerz/yamnet-plugin


text:
# -*- coding: utf-8 -*-
"""
Created on Tue May  5 12:44:03 2020

@author: mahjaf
"""

# -*- coding: utf-8 -*-
"""
Created on Wed Apr 22 10:07:37 2020

@author: mahjaf

Using this code, one can directly feed in EDF data and select channels of interest
to perform classification.

Please Note: we recommend to use "EDF_to_h5.py" to firstly convert EDF into
a lighter data format to accelerate computations; however, this code is only
meant for those that want to skip this conversion and directly choose EDF data
as input.
"""
#%% Reading EDF section
#####===================== Importiung libraries =========================#####

import mne
import numpy as np
from   numpy import loadtxt
import h5py
import time
import os
from ssccoorriinngg import ssccoorriinngg
import matplotlib.pyplot as plt
import pickle
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import RandomizedSearchCV
from sklearn.metrics import confusion_matrix, make_scorer, accuracy_score, precision_score, recall_score, f1_score

#####===================== Defining project folder========================#####

project_folder = "C:/Users/mahda/Documents/ssccoorriinngg/required_info/"
#project_folder = "P:"

#####===================== Reading EDF data files=========================#####

pat_labels = loadtxt(project_folder + "patient_labels.txt", delimiter="\t", skiprows = 1)

#####============= Distinguishing patients from control group=============#####

gp = loadtxt(project_folder + "grouping.txt", delimiter="\t", skiprows = 1, dtype = 'str')
subj_c = [] # Control
subj_p = [] # Patients

for indx, c in enumerate(gp):
    if c[1] == 'C':
        subj_c.append(int(c[0]))
    elif c[1] == 'CC':
        pass
    else:
        subj_p.append(int(c[0]))

# Initialization
subjects_dic     = {}
hyp_dic          = {}
metrics_per_fold = {}

#####============= create an object of ssccoorriinngg class ==============#####

Object = ssccoorriinngg(filename='', channel='', fs = 200, T = 30)
"""

tic_tot = time.time()

# Igonre unnecessary warnings
np.seterr(divide='ignore', invalid='ignore')

#####============= Iterate through each subject to find data =============#####

for idx, c_subj in enumerate(subj_c):
    print (f'Analyzing Subject Number: {c_subj}')
    ## Read in data
    file     = project_folder + "/3013080.02/ml_project/test_data/LK_" + str(int(c_subj)) + "_1.EDF"
    tic      = time.time()
    data     = mne.io.read_raw_edf(file)

    raw_data = data.get_data()
    print('Time to read EDF: {}'.format(time.time()-tic))

#####=================Retrieving information from data====================#####

    DataInfo          = data.info
    AvailableChannels = DataInfo['ch_names']
    fs                = int(DataInfo['sfreq'])

#####==================Choosing channels of interest======================#####

    # 1. The channels that need to be referenced
    Mastoids          = ['Fp2'] # Reference electrodes
    RequiredChannels  = ['Fp1'] # main electrodes
    # 2. Channels that don't need to be referenced: --> Deactive

    Idx               = []
    Idx_Mastoids      = []

#####================= Find index of required channels ===================#####

    for indx, c in enumerate(AvailableChannels):
        if c in RequiredChannels:
            Idx.append(indx)
        elif c in Mastoids:
            Idx_Mastoids.append(indx)

#####===== Sampling rate is 200hz; thus 1 epoch(30s) is 6000 samples =====#####

    T = 30 #secs
    len_epoch   = fs * T
    start_epoch = 0
    n_channels  =  len(AvailableChannels)

#####============ Cut tail; use modulo to find full epochs ===============#####

    raw_data = raw_data[:, 0:raw_data.shape[1] - raw_data.shape[1]%len_epoch]

#####========== Reshape data [n_channel, len_epoch, n_epochs] ============#####
    data_epoched = np.reshape(raw_data,
                              (n_channels, len_epoch,
                               int(raw_data.shape[1]/len_epoch)), order='F' )

#####===================== Reading hypnogram data ========================#####

    hyp = loadtxt(project_folder +"/3013065.04/Depressed_Loreta/hypnograms/LK_" +
                str(int(c_subj)) + ".txt", delimiter="\t")

    ### Create sepereate data subfiles based on hypnogram (N1, N2, N3, NREM, REM)
    tic      = time.time()
#####================= Concatenation of selected channels ================#####
    # Calculate referenced channels:
    data_epoched_selected = data_epoched[Idx] - data_epoched[Idx_Mastoids]
    # show picked channels for analysis
    picked_channels = []
    for jj,kk in enumerate(Idx):
        picked_channels = np.append(picked_channels, AvailableChannels[kk])

    print(f'subject LK {c_subj} --> detected channels: {picked_channels}')
    # Add non-referenced channels:  --> Deactive
    #data_epoched_selected_ = np.concatenate([data_epoched_selected, data_epoched_nonref], 0)
    print('Time to split sleep stages per epoch: {}'.format(time.time()-tic))

    #%% Analysis section
#####================= remove bad chanbnels and arousals =================#####

    # assign the proper data and labels
    x_tmp_init = data_epoched_selected
    y_tmp_init = hyp
    # Remove bad signals (hyp == -1)
    x_tmp, y_tmp =  Object.remove_bad_signals(hypno_labels = y_tmp_init,
                                              input_feats = x_tmp_init)

#####============= remove stages contaminated with arousal ===============#####

    x_tmp, y_tmp = Object.remove_arousals(hypno_labels = y_tmp, input_feats = x_tmp)
    # Create binary labels array
    yy = Object.binary_labels_creator_categories(y_tmp)

    # Initialize feature array:
    Feat_all_channels = np.empty((np.shape(x_tmp)[-1],0))

#####================== Extract the relevant features ====================#####

    for k in np.arange(np.shape(data_epoched_selected)[0]):

        feat_temp         = Object.FeatureExtraction_per_subject(Input_data = x_tmp[k,:,:])
        Feat_all_channels = np.column_stack((Feat_all_channels,feat_temp))

    toc = time.time()
    print(f'Features of subject {c_subj} were successfully extracted in: {toc-tic} secs')

    # Defining dictionary to save features PER SUBJECT
    subjects_dic["subject{}".format(c_subj)] = Feat_all_channels

    # Defining dictionary to save hypnogram PER SUBJECT
    hyp_dic["hyp{}".format(c_subj)] = yy

    #####============================= Save per case =========================#####
    path     = project_folder +"/3013080.02/ml_project/scripts/1D_TimeSeries/features/percase"
    filename = "/LK__"+ str(int(c_subj))
    Object.save_dictionary(path, filename, hyp_dic, subjects_dic)

#####=============== Removing variables for next iteration ===============#####
    del x_tmp, y_tmp, feat_temp, yy
    toc = time.time()
    print(f'Features and hypno of subject {c_subj} were successfully added to dictionary')

    print('Feature extraction of subject {c_subj} has been finished.')

print('Total feature extraction of subjects took {tic_tot - time.time()} secs.')
#%% Save created features and labels

#####====================== Save extracted features ======================#####

path     = project_folder +"/3013080.02/ml_project/scripts/1D_TimeSeries/features/"
filename = 'sleep_scoring_NoArousal_Fp1-Fp2_full_'
Object.save_dictionary(path, filename, hyp_dic, subjects_dic)

#%% Load featureset and labels
"""
path                  =  project_folder
filename              = "sleep_scoring_NoArousal_Fp1-Fp2_080520"
subjects_dic, hyp_dic = Object.load_dictionary(path, filename)

#%% ================================Training part==============================

# Training perentage
train_size = .75
n_train = round(train_size * len(subj_c))

#######=== Randomly shuffle subjects to choose train and test splits ===#######

subj_c = np.random.RandomState(seed=42).permutation(subj_c)

#######=============== Initialize train and test arrays ================#######

X_train = np.empty((0, np.shape(subjects_dic['subject14'])[1]))
X_test  = np.empty((0, np.shape(subjects_dic['subject14'])[1]))
y_train = np.empty((0, np.shape(hyp_dic['hyp14'])[1]))
y_test  = np.empty((0, np.shape(hyp_dic['hyp14'])[1]))

########======= Picking the train subjetcs and concatenate them =======########
tic = time.time()

for c_subj in subj_c[0:n_train]:

    # train hypnogram
    str_train_hyp  = 'hyp' + str(c_subj)

    # train featureset
    str_train_feat = 'subject' + str(c_subj)

    # create template arrays for featurs and label
    tmp_x          =  subjects_dic[str_train_feat]
    tmp_y          =  hyp_dic[str_train_hyp]

    # Concatenate features and labels
    X_train = np.row_stack((X_train, tmp_x))
    y_train = np.row_stack((y_train, tmp_y))
    del tmp_x, tmp_y
print('Training set was successfully created in : {} secs'.format(time.time()-tic))

#%% ================================Test part==============================%%#

########======== Picking the test subjetcs and concatenate them =======########
tic = time.time()

for c_subj in subj_c[n_train:]:

    # test hypnogram
    str_test_hyp  = 'hyp' + str(c_subj)

    # test featureset
    str_test_feat = 'subject' + str(c_subj)

    # create template arrays for featurs and  label
    tmp_x         =  subjects_dic[str_test_feat]
    tmp_y         =  hyp_dic[str_test_hyp]

    # Concatenate features and labels
    X_test = np.row_stack((X_test, tmp_x))
    y_test = np.row_stack((y_test, tmp_y))
    del tmp_x, tmp_y
print('Test set was successfully created in : {} secs'.format(time.time()-tic))

print(f'Raw train and test data were created.')

########================== Replace any probable NaN ===================########

X_train = Object.replace_NaN_with_mean(X_train)
X_test  = Object.replace_NaN_with_mean(X_test)

########================== Replace any probable inf ===================########

X_train = Object.replace_inf_with_mean(X_train)
X_test  = Object.replace_inf_with_mean(X_test)

########==================== Z-score of features ======================########

X_train, X_test = Object.Standardadize_features(X_train, X_test)

########========== select features only on first iteration ============########

ranks, Feat_selected, selected_feats_ind = Object.FeatSelect_Boruta(X_train,
                                                    y_train[:,0], max_depth = 7)

########=================== Apply selected features ===================########

X_train = X_train[:, selected_feats_ind]
X_test  = X_test[:, selected_feats_ind]

########======= Add time dependence to the data classification ========########

td = 6 # Time dependence: # epochs of memory
X_train_td = Object.add_time_dependence_to_features(X_train, n_time_dependence=td)
X_test_td  = Object.add_time_dependence_to_features(X_test,  n_time_dependence=td)

########======== Temporary truncate first and last "td" epochs ========########

X_train_td = X_train_td[td:len(X_train_td)-td,:]
X_test_td  = X_test_td[td:len(X_test_td)-td,:]
y_train_td    = y_train[td:len(y_train)-td,:]
y_test_td     = y_test[td:len(y_test)-td,:]

########====================== Feature Selection ======================########

y_train_td = Object.binary_to_single_column_label(y_train_td)

########============== Define classifier of interest ==================########

'''y_pred = Object.RandomForest_Modelling(X_train, y_train, X_test, y_test, n_estimators = 10)
y_pred = Object.Extra_randomized_trees(X_train, y_train, X_test,y_test, n_estimators= 1000, max_depth = None, min_samples_split =2,
                               max_features="sqrt")

#y_pred = Object.ADAboost_Modelling(X_train, y_train,X_test, y_test, n_estimators = 1000)
y_pred = Object.gradient_boosting_classifier(X_train, y_train,X_test, y_test,
                                     n_estimators = 500, learning_rate= 1.0, max_depth=1)
'''

y_pred = Object.XGB_Modelling(X_train, y_train,X_test, y_test, n_estimators = 250,
                      max_depth=3, learning_rate=.1)

########===== Metrics to assess the model performance on test data ====########

Acc, Recall, prec, f1_sc,cm = Object.multi_label_confusion_matrix(y_test, y_pred)

########====================== Plot hypnogram =========================########

Object.plot_confusion_matrix(y_test, y_pred, target_names = ['Wake','N1','N2','SWS','REM'],
                          title='Confusion matrix of ssccoorriinngg algorithm',
                          cmap='jet',
                          normalize=True)

########=============== Concatenate metrics to store ==================########

all_metrics = [Acc, Recall, prec, f1_sc]

########========================== Hypnogram ==========================########

hyp_test = Object.binary_to_single_column_label(y_test)
Object.plot_hyp(hyp = hyp_test, mark_REM = 'active')

########================== Comparative hypnogram ======================########

hyp_pred = Object.binary_to_single_column_label(y_pred)
Object.plot_comparative_hyp(hyp_true = hyp_test, hyp_pred = hyp_pred, mark_REM = 'active')

########========================== Save figure =======================#########
Object.save_figure(saving_format = '.png',
                   directory = '/project/3013080.02/Mahdad/Github/ssccoorriinngg/Plots/v0.2/Fp1-Fp2/',
                   saving_name = 'test_subject_all' + str(c_subj), dpi = 900,
                   full_screen = False)
"""
#########======================= Randomized search ===================#########
# Define sccoring metrics
scoring = {'accuracy' : make_scorer(accuracy_score),
           'precision' : make_scorer(precision_score),
           'recall' : make_scorer(recall_score),
           'f1_score' : make_scorer(f1_score)}
# Apply randomized search
Object.RandomSearchRF(X = X_train, y = y_train, scoring = scoring, estimator = RandomForestClassifier(),
                        n_estimators = [int(x) for x in np.arange(10, 1000, 50)],
                        max_features = ['log2', 'sqrt'],
                        max_depth = [int(x) for x in np.arange(10, 100, 30)],
                        min_samples_split = [2, 5, 10],
                        min_samples_leaf = [1, 2, 4],
                        bootstrap = [True, False],
                        n_iter = 100, cv = 10)
#%% Save results
path = '/project/3013080.02/Mahdad/Github/ssccoorriinngg/Plots/v0.1/Fp1-Fp2/'
fname = 'LOOCV_results_Fp1-Fp2'
with open(path+fname+'.pickle',"wb") as f:
            pickle.dump(metrics_per_fold, f)

#%% load results

path     = 'P:/3013080.02/Mahdad/Github/ssccoorriinngg/Plots/v0.1/'
filename = 'LOOCV_results'
with open(path + filename + '.pickle', "rb") as f:
   Outputs = pickle.load(f)
#%% Show final average results
Object.Mean_leaveOneOut(Outputs)
"""

max_stars_repo_name: MahdadJafarzadeh/ssccoorriinngg

max_issues_repo_name: MahdadJafarzadeh/ssccoorriinngg

max_forks_repo_name: MahdadJafarzadeh/ssccoorriinngg


text:
"""Plot to test line styles"""
import matplotlib.pyplot as plt
import numpy as np
import mpld3
from mpld3 import plugins

def create_plot():
    fig, ax = plt.subplots()
    points = ax.plot(range(10), 'o', ms=20)
    plugins.connect(fig, plugins.PointLabelTooltip(points[0],
                                                   location="top left"))
    return fig

def test_tooltips_basic():
    fig = create_plot()
    html = mpld3.fig_to_html(fig)
    plt.close(fig)

if __name__ == "__main__":
    mpld3.show(create_plot())


text:
import argparse
from PIL import Image, ImageDraw, ImageFilter
import random
from scipy.spatial import distance
import csv
import os
import numpy as np
import sys
import glob

parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument('--data_case', help='either bright-field (bf) or phase-contrast (pc)', type=str, default="pc")
parser.add_argument('--ratio', help='multiplier of generated data', type=float, default=1.0)
args = parser.parse_args()

data_type = args.data_case #bf for bright-field  # pc for phase-contrast

fill_color = [180,180,180] if data_type=="bf" else [50,50,50] #default: 180 for bf, 50 for pc
fill_fluctuation = 30 if data_type=="bf" else 10 #default: 30 for bf, 10 for pc
edge_color = [80,80,80] if data_type=="bf" else [160,160,160] #default: bf: 80 default, 40 high_contrast, pc: 160 for 0.1 blur, 190 for 0.5 blur, 200 for 0.75 blur, 210 for 1 blur, +20 high contrast
edge_fluctuation = 50
SIZE = 128 #and height
num_circles_min = 1
num_circles_max = 31 #26 #exclusive, set to one higher than needed
r_min = 3
r_max = 8 if data_type=="bf" else 10 #default: 8 for bf, 10 for pc
crop_x_left = 26 if data_type=="bf" else 12 #default: 26 for bf, 12 for pc
crop_x_right = 20 if data_type=="bf" else 11 #default: 20 for bf, 11 for pc
crop_y_top = 12 if data_type=="bf" else 11 #default: 12 for bf, 11 for pc
crop_y_bottom = 6 if data_type=="bf" else 10 #default: 6 for bf, 10 for pc
allowedOverlap = 0 if data_type=="bf" else 1 #defautl: 0 for bf, 1 for pc
deformed_min = -2
deformed_max = 2

train_images_count = np.round(2469*args.ratio) if data_type=="bf" else np.round(3152*args.ratio) #2469 for bf train, 514 for test  # 3182 for pc train, 794 for test
test_images_count = 1000 #no constraint to match nat test data amount, but if wanted: 514 if data_type=="bf" else 792
foldername = "../data/128p_"+data_type+"_dyn"
bluriness = 0.5 #default: 0.5 for both
proportion_mod = 5

os.makedirs(foldername, exist_ok=True)
x = 0
y = 0
r = 0
created_circles = []

def get_white_noise_image(size):
	noise_map = Image.new("RGBA", (size, size), 255)
	random_grid = map(lambda x: (int(random.random() * 256),int(random.random() * 256),int(random.random() * 256)), [0] * size * size)
	noise_map.putdata(list(random_grid))
	noise_map = noise_map.convert("L") #conversion to luminosity noise
	noise_map = noise_map.convert("RGBA")
	return noise_map

def generate_images(numImages, isTrain):
	generated_images = 0

	filename = ""
	if isTrain:
		filename = '/train_gen.csv'
	else:
		filename = '/test_gen.csv'

	with open(foldername + filename, mode='w', newline='') as mycsv:
		writer = csv.writer(mycsv, delimiter=',')

		props = list(range(num_circles_min, num_circles_max))
		props.reverse()
		proportions = [j+proportion_mod for j in props]
		total = sum(proportions)

		images_folder = ""
		if isTrain:
			images_folder = "/train_gen/"
		else:
			images_folder = "/test_gen/"
		os.makedirs(foldername+images_folder, exist_ok=True)
		for f in glob.glob(foldername+images_folder+"*"):
			os.remove(f)

		while(generated_images < numImages):
			c = 0
			image = Image.open("../data_generator/empty128.png")
			image = image.convert("RGBA")
			draw = ImageDraw.Draw(image)
			i = 0
			 #num_circles = random.randrange(num_circles_min, num_circles_max)

			sel = random.randrange(1, total)
			z = 0
			while sel > 0:
				sel -= proportions[z]
				z += 1
			num_circles = z

			 #num_circles = random.randrange(1,466) #1+...+40=(1+40)*20=820 or 1+...+31=(1+31)*15=465? +1 for numerical reasons
			 #c = 1
			 #while num_circles > num_circles_max-c:
			 #	if num_circles_max-c < 0:
			 #		print("error in num_circles calculation")
			 #	num_circles -= (num_circles_max-c)
			 #	c += 1
			 #num_circles = c

			while(i < num_circles):
				r = random.randrange(r_min, r_max)
				x = random.randrange(r + crop_x_left, SIZE - r - crop_x_right, 2)
				y = random.randrange(r + crop_y_top, SIZE - r - crop_y_bottom, 2)
				accectable_position = True
				for elem in created_circles:
					dist = distance.euclidean([x,y],[elem[0],elem[1]])
					if(dist < r+elem[2]-allowedOverlap):
						accectable_position = False
						 #print("invalid position")
						break
				if(accectable_position):
					deform_x = random.randrange(deformed_min, deformed_max)
					deform_y = random.randrange(deformed_min, deformed_max)
					if r <=3:
						deform_x /= 2
						deform_y /= 2
					side_selector = random.randrange(0,1)
					if side_selector == 0:
						left = x-r+deform_x
						right = x+r
					else:
						left = x-r
						right = x+r+deform_x
					side_selector = random.randrange(0,1)
					if side_selector == 0:
						top = y-r+deform_y
						bottom = y+r
					else:
						top = y-r
						bottom = y+r+deform_y
					b = random.randrange(0, edge_fluctuation)
					e_color = np.copy(edge_color)
					e_color[0] += b
					e_color[1] += b
					e_color[2] += b
					draw.ellipse((left, top, right, bottom), fill=tuple(e_color))
					b = random.randrange(0, fill_fluctuation)
					cell_color = np.copy(fill_color)
					cell_color[0] += b
					cell_color[1] += b
					cell_color[2] += b
					draw.ellipse((left + 1, top + 1, right - 1, bottom - 1), fill=tuple(cell_color))
					i += 1
					created_circles.append([x,y,r])

			file = foldername + images_folder + str(SIZE) + "_" + str(generated_images) + ".png"

			blurred = image.filter(ImageFilter.GaussianBlur(radius=1))
			mean = Image.open("../data_generator/mean" + str(SIZE) + "_" + data_type + ".png")
			mean = mean.convert("RGBA")

			 #SHOULD NOT BE REQUIRED DUE TO NOISE OPTIONS IN dataloader.py
			 #noise = get_white_noise_image(SIZE)
			 #blend = Image.blend(image, noise, 0.1) #hellfeld
			 #blend = Image.blend(image, noise, 0.0) #pc

			 #blend = Image.blend(image, blurred, bluriness) #hellfeld
			blend = Image.blend(image, blurred, bluriness) #pc
			blend = blend.convert("RGBA")

			result = Image.alpha_composite(mean, blend)

			result.save(file)
			writer.writerow([file, num_circles])
			generated_images += 1
			created_circles.clear()

			c += 1

generate_images(train_images_count, True)
generate_images(test_images_count, False)


text:
#!/usr/bin/env python

import numpy as np
from scipy import special

from ..routines import median, mahalanobis, gamln, psi

from nose.tools import assert_true
from numpy.testing import assert_almost_equal, assert_equal, TestCase

class TestAll(TestCase):

    def test_median(self):
        x = np.random.rand(100)
        assert_almost_equal(median(x), np.median(x))

    def test_median2(self):
        x = np.random.rand(101)
        assert_equal(median(x), np.median(x))

    def test_median3(self):
        x = np.random.rand(10, 30, 11)
        assert_almost_equal(np.squeeze(median(x,axis=1)), np.median(x,axis=1))

    def test_mahalanobis(self):
        x = np.random.rand(100) / 100
        A = np.random.rand(100, 100) / 100
        A = np.dot(A.transpose(), A) + np.eye(100)
        mah = np.dot(x, np.dot(np.linalg.inv(A), x))
        assert_almost_equal(mah, mahalanobis(x, A), decimal=1)

    def test_mahalanobis2(self):
        x = np.random.rand(100,3,4)
        Aa = np.zeros([100,100,3,4])
        for i in range(3):
            for j in range(4):
                A = np.random.rand(100,100)
                A = np.dot(A.T, A)
                Aa[:,:,i,j] = A
        i = np.random.randint(3)
        j = np.random.randint(4)
        mah = np.dot(x[:,i,j], np.dot(np.linalg.inv(Aa[:,:,i,j]), x[:,i,j]))
        f_mah = (mahalanobis(x, Aa))[i,j]
        assert_true(np.allclose(mah, f_mah))

    def test_gamln(self):
        for x in (0.01+100*np.random.random(50)):
            scipy_gamln = special.gammaln(x)
            my_gamln = gamln(x)
            assert_almost_equal(scipy_gamln, my_gamln)

    def test_psi(self):
        for x in (0.01+100*np.random.random(50)):
            scipy_psi = special.psi(x)
            my_psi = psi(x)
            assert_almost_equal(scipy_psi, my_psi)

if __name__ == "__main__":
    import nose
    nose.run(argv=['', __file__])

max_stars_repo_name: yarikoptic/nipy

max_issues_repo_name: fabianp/nipy

max_forks_repo_name: fabianp/nipy

reason: import numpy,from numpy,from scipy


text:
# utility functions for som processing
import numpy as np
import math
from math import sqrt

def parseBMFile(bmfile):
    bmfp = open(bmfile)
    nrow, ncol = bmfp.readline()[1:-1].split()
    nrow = int(nrow)
    ncol = int(ncol)
    nsamples = int(bmfp.readline()[1:-1])
    bmus = [None]*nsamples
    for line in bmfp:
        data = line[:-1].split()
        bmus[int(data[0])] = [int(data[1]),int(data[2])]
    bmfp.close()
    return nrow, ncol, bmus

def mapBMUs(bmfile):
	bmfp = open(bmfile)
	nrow, ncol = bmfp.readline()[1:-1].split()
	nrow = int(nrow)
	ncol = int(ncol)
	nsamples = int(bmfp.readline()[1:-1])
	bmuMap = [[[] for x in range(ncol)] for y in range(nrow)]
	for line in bmfp:
		data = line[:-1].split()
		bmuMap[int(data[1]),int(data[2])].append(int(data[0]))
	return nrow, ncol, bmuMap

def parseLRNHeader(lrnfile):
    lrnfp = open(lrnfile)
    samplenames = lrnfp.readline()[1:-1].split()
    lrnfp.readline()
    lrnfp.readline()
    lmask = lrnfp.readline()[1:-1].split()
    featnames = lrnfp.readline()[1:-1].split()
    featnames = [featnames[i] for i in range(len(featnames)) if lmask[i] == '1']
    lrnfp.close()
    return samplenames, featnames

def parseWTSFile(wtsfile):
    wtsfp = open(wtsfile)
    nrow,ncol = wtsfp.readline()[1:-1].split()
    numfeatures = int(wtsfp.readline()[1:-1])
    nrow,ncol = int(nrow),int(ncol)

    map = [[np.zeros(numfeatures) for y in range(ncol)] for x in range(nrow)]
    for k in range(nrow):
        for l in range(ncol):
            map[k][l] = np.asfarray(wtsfp.readline()[:-1].split(), float)
    wtsfp.close()
    return nrow,ncol,numfeatures,map

"""
def calculateQuantizationError(lrnfile, wtsfile, bmfile):
    bnrow, bncol, bmus = parseBMFile(bmfile)
    nrow,ncol,nfeatures,map = parseWTSFile(wtsfile)
    lrnfp = open(lrnfile)
    samplenames = lrnfp.readline()[1:-1].split()
    lrnfp.readline()
    lrnfp.readline()
    lmask = lrnfp.readline()[1:-1].split()
    featnames = lrnfp.readline()[1:-1].split()
    featnames = [featnames[i] for i in range(len(featnames)) if lmask[i] == '1']
    error = 0.
    for line in lrnfp:
        data = line[:-1].split()
        index = int(data[0])
        sample = np.asfarray(data[1:], float)
        error += np.linalg.norm(map[bmu[index][0]][bmu[index][1]] - sample)
        print(error)
        if index % 100 == 0:
            print(index)
    return error
"""
def getbmu1and2(map, nrow, ncol, sample):
    bmu1 = [0,0]
    bmu2 = [0,0]
    kmin, lmin, kmin2, lmin2 = 0,0,0,0
    distmin, distmin2 = 1.e6, 1.e6

    for k in range(nrow):
        for l in range(ncol):
            dist = np.linalg.norm(map[k][l] - sample)
            if dist < distmin:
                distmin2 = distmin
                kmin2 = kmin
                lmin2 = lmin
                distmin = dist
                kmin = k
                lmin = l
    return distmin, [kmin, lmin], distmin2, [kmin2, lmin2]

def areNeighbors(x1, y1, x2, y2, nrow, ncol, topo):
    if topo == "rectangular":
        if (abs(x1 - x2) + abs(y1-y2)) > 1:
            return 0 # not neighbors
        else:
            return 1
    elif topo == "rectangular-toroid":
        if (min(abs(x1-x2), nrow - abs(x1-x2)) + min(abs(y1-y2), ncol - abs(y1-y2))) > 1:
                return 0
        else:
            return 1

    elif topo == "hexagonal":
        xdist = abs(x1-x2)
        ydist = abs(y1-y2)
        ymin = min(y1,y2)
        if (ydist & 1) :
            if (ymin & 1):
                xdist -= 0.5
            else:
                xdist += 0.5
        dist = sqrt( xdist*xdist + ydist * ydist * 0.75)
        if dist > 1.:
            return 0
        else:
            return 1

    elif topo == "hexagonal-toroid":
        xdist = min(abs(x1 - x2), abs(x1-x2+nrow))
        ydist = min(abs(y1-y2), abs(y1-y2+ncol))
        ymin = min(y1, y2)
        if (ydist & 1) :
            if (ymin & 1):
                xdist -= .05
            else:
                xdist += 0.5
        dist = sqrt(xdist*xdist + ydist*ydist*0.75)
        if dist > 1.:
            return 0
        else:
            return 1

def calculateErrors(lrnfile, wtsfile, topo):
    # since need to recalculate bmus anyway, don't use bmfile for this
    nrow,ncol,nfeatures,map = parseWTSFile(wtsfile)
    lrnfp = open(lrnfile)
    samplenames = lrnfp.readline()[1:-1].split()
    lrnfp.readline()
    lrnfp.readline()
    lmask = np.asfarray(lrnfp.readline()[1:-1].split(), int)
    featnames = np.array(lrnfp.readline()[1:-1].split())
    featnames = featnames[lmask==1]

    # account for masking
    quantizationError = 0.
    topologicalError = 0
    index = 0
    for line in lrnfp:
        sample = np.asfarray(line[:-1].split(), float)
        sample = sample[lmask==1]
        distmin, bmu1, distmin2, bmu2 = getbmu1and2(map, nrow, ncol, sample)
        quantizationError += distmin
        topologicalError += (1 - areNeighbors(bmu1[0], bmu1[1], bmu2[0], bmu2[1], nrow, ncol, topo))
        index += 1
        if index % 500 == 0:
            print("Processed "+str(index)+" out of "+str(len(samplenames))+" sample errors.")
    topologicalError /= (len(samplenames))
    quantizationError /= (len(samplenames))

    return quantizationError, topologicalError


text:
"""
#;+
#; NAME:
#; fits
#;    Version 1.0
#;
#; PURPOSE:
#;    Module for FITS utilities
#;   10-Dec-2014 by JXP
#;-
#;------------------------------------------------------------------------------
"""
from __future__ import print_function, absolute_import, division, unicode_literals

import numpy as np
import gzip
import subprocess, glob

from astropy import units as u
from astropy.io import ascii, fits
from astropy.table import QTable, Table, Column

from xastropy.xutils import xdebug as xdb

# def bintab_to_table(fits_fil,exten=1, silent=False):
# def table_to_fits(table, outfil, compress=False, comment=None):

#
def bintab_to_table(fits_fil,exten=1, silent=False):
    ''' Read a binary FITS table into an astropy QTable

    Parameters
    ---------
    fits_fil: sting
      File
    exten: int (1)
      Extension of the FITS file
    '''
    # Look for file
    fil = glob.glob(fits_fil+'*')
    if len(fil) == 1:
        if not silent:
            print('x.fits.bintab_to_table: Reading {:s}'.format(fil[0]))
        infil = fil[0]
    elif len(fil) == 0:
        raise IOError('x.fits.bintab_to_table: No file {:s}'.format(fits_fil))
    else:
        # Taking specified one
        infil = fits_fil

    return QTable( fits.open(infil)[exten].data )

#
def table_to_fits(itable, outfil, compress=False, comment=None):
    ''' Write an astropy Table as a FITS binary table
    ---
    Parameters
    ---------
    table: astropy.Table
    outfil: string
    compress: bool (False)
       gzip compress?
    comment: string
       Comment to insert into the header
    '''

    # Comments
    table = Table(itable)
    if not comment is None:
        table.meta['comments'] = [comment]
    table.write(outfil, format='fits', overwrite=True)

    # Compress?
    if compress:
        subprocess.call(["gzip", "-f", outfil])

max_stars_repo_name: mneeleman/xastropy

max_issues_repo_name: mneeleman/xastropy

max_forks_repo_name: mneeleman/xastropy

reason: import numpy,from astropy


text:
"""
feeder. polls exchanges and publishes
"""

from datetime import datetime

import archon.broker as broker
import archon.exchange.exchanges as exc
import archon.exchange.bitmex.bitmex as mex
import archon.exchange.deribit.Wrapper as deri
import archon.broker.facade as facade
import archon.model.models as models
from archon.util.custom_logger import setup_logger, remove_loggers

from archon.util import *

import pandas as pd
import numpy

import time
import json
import atexit
import logging
import threading
import time
import redis
from .topics import *

#default_symbols = {exc.BITMEX:mex.instrument_btc_perp}

class Feeder(threading.Thread):

    def __init__(self, abroker):
        threading.Thread.__init__(self)
        setup_logger(logger_name="Feeder", log_file='feeder.log')
        self.log = logging.getLogger("Feeder")
        self.abroker = abroker
        self.log.info("init feeder")
        self.redisclient = abroker.get_redis()

        #TODO in config
        mex_sym = mex.instrument_btc_perp
        self.mex_sym = mex_sym

    def pub_set(self, topic, data):
        """ publish and set """
        d = {"topic":topic,"data":data}
        jdata = json.dumps(d)
        self.redisclient.publish(topic, jdata)
        t = rep + topic[4:]
        self.redisclient.set(t, jdata)

        db = self.abroker.get_db()
        tt = topic[4:]
        #store
        db[tt].insert_one(d)

    def publish_bitmex(self):
        self.log.info("publish_bitmex")
        wait = 0.1 # bitmex rate limit 300 per 300 seconds
        pos = self.abroker.afacade.position(exc.BITMEX)
        pos = {"position": pos}
        time.sleep(wait)
        self.log.info("pos %s"%pos)
        #print ("position ",pos)
        #if pos == []: pos = {}

        self.pub_set(SUB_TOPIC_POS_BITMEX, pos)

        oo = self.abroker.afacade.openorders(exc.BITMEX, mex.instrument_btc_perp)
        time.sleep(wait)
        self.pub_set(SUB_TOPIC_ORDERS_BITMEX, oo)

        book = self.abroker.afacade.orderbook(self.mex_sym, exc.BITMEX)
        time.sleep(wait)
        self.pub_set(SUB_TOPIC_MARKET_BOOK_BITMEX, book)

    def publish_facil(self):
        pass

    def publish_deribit(self):
        pos = self.abroker.afacade.position(exc.DERIBIT)
        if pos == None or pos == []: pos = {}
        self.pub_set(SUB_TOPIC_POS_DERIBIT, pos)

        oo = self.abroker.afacade.openorders(exc.DERIBIT, deri.instrument_btc_perp)
        if oo == None or oo == []: oo = {}
        self.pub_set(SUB_TOPIC_ORDERS_DERIBIT, oo)
        data = self.abroker.afacade.orderbook(deri.instrument_btc_perp, exc.DERIBIT)
        self.pub_set(SUB_TOPIC_MARKET_BOOK_DERIBIT, data)

    def run(self):
        self.log.info("start feeder")
        while True:
            self.log.info("feeder loop")
            #TODO openorders
            #TODO position

            self.publish_bitmex()
            #self.publish_deribit()

            time.sleep(0.5)

max_stars_repo_name: tebzzz/archon

max_issues_repo_name: tebzzz/archon

max_forks_repo_name: tebzzz/archon


text:
#  Copyright (C) 2014-2021 Syntrogi Inc dba Intheon. All rights reserved.

from collections import deque, namedtuple
import json
import numpy as np
from qtpy import QtGui
import pyqtgraph as pg
from stream_viewer.renderers.data.base import RendererDataTimeSeries
from stream_viewer.renderers.display.pyqtgraph import PGRenderer

MarkerMap = namedtuple('MarkerMap', ['source_id', 'timestamp', 'item'])

class LinePG(RendererDataTimeSeries, PGRenderer):

    plot_modes = ["Sweep", "Scroll"]
    gui_kwargs = dict(RendererDataTimeSeries.gui_kwargs, **PGRenderer.gui_kwargs,
                      offset_channels=bool, reset_colormap=bool,
                      line_width=float, antialias=bool, ylabel_as_title=bool,
                      ylabel_width=int)

    def __init__(self,
                 # Override inherited
                 auto_scale: str = 'none',
                 show_chan_labels: bool = True,
                 color_set: str = 'viridis',
                 # New
                 offset_channels: bool = True,
                 reset_colormap: bool = False,
                 line_width: float = 2.0,
                 antialias: bool = True,
                 ylabel_as_title: bool = False,
                 ylabel_width: int = None,
                 ylabel: str = None,
                 **kwargs):
        """
        Multi-channel timeseries visualization using pyqtgraph widgets. Channels originating from the same data source
         will be grouped together in a single axes.
        This is a slower but more flexible alternative to `LineVis`.

        Args:
            offset_channels: Set True to force channels into separate y-ranges within the same plot.
                Otherwise the channels are free to overlap.
            reset_colormap: Set True to reset the colormap iterator between data sources. i.e. channel i will have the
                same color for all data sources.
            line_width: Width of all lines.
            antialias: Set True to enable anti-aliasing. I haven't noticed much effect.
            ylabel_as_title: The ylabel is removed and instead a title is created with typical ylabel contents.
            ylabel_width: Force a minimum amount of real-estate to be allocated to the ylabel. This overcomes a problem
                with the required y-width not being calculated until after all the labels and ticks have been
                created.
            ylabel: The ylabel for the plot. If unspecified, the name of the LSL stream is automatically used.
            **kwargs:
        """
        self._offset_channels = offset_channels
        self._reset_colormap = reset_colormap
        self._line_width = line_width
        self._antialias = antialias
        self._ylabel_as_title = ylabel_as_title
        self._ylabel_width = ylabel_width
        self._ylabel = ylabel
        self._requested_auto_scale = auto_scale.lower()  # Actual auto-scale is different depending on n streams.
        self._widget = pg.GraphicsLayoutWidget()
        self._do_yaxis_sync = False
        self._src_last_marker_time = []
        self.marker_texts_pool = deque()
        self._marker_info = deque()
        self._t_expired = -np.inf  # Anything items (e.g. marker strings) older than this can be removed safely.
        super().__init__(show_chan_labels=show_chan_labels, color_set=color_set, **kwargs)
        self.reset_renderer()

    def reset_renderer(self, reset_channel_labels=True):
        # Clear existing elements
        self._widget.clear()
        self._widget.setBackground(self.parse_color_str(self.bg_color))
        self._src_last_marker_time = [-np.inf for _ in range(len(self._data_sources))]

        if len(self.chan_states) == 0:
            return

        chans_per_src = self.chan_states.loc[self.chan_states['vis']].groupby('src')['name'].nunique().values
        n_chans_colormap = np.max(chans_per_src) if self.reset_colormap else np.sum(chans_per_src)
        color_map = self.get_colormap(self.color_set, n_chans_colormap)

        labelStyle = {'color': '#FFF', 'font-size': str(self.font_size) + 'pt'}

        if self._requested_auto_scale == 'all' and (not self.offset_channels or np.all(chans_per_src <= 1)):
            # When doing a global auto scale, we don't actually need it if we are not offsetting channels or each
            #  source has at most 1 channel, because then we can rely on pg's auto scaling.
            self._auto_scale = 'none'
        else:
            self._auto_scale = self._requested_auto_scale

        row_offset = -1
        ch_offset_color = -1
        last_row = 0
        for src_ix, src in enumerate(self._data_sources):
            ch_states = self.chan_states[self.chan_states['src'] == src.identifier]
            n_vis_src = ch_states['vis'].sum()
            if n_vis_src == 0:
                continue

            offset_chans = self.offset_channels and n_vis_src > 1

            buff = self._buffers[src_ix]
            n_samples = buff._data.shape[-1]  # int(np.ceil(stats['srate'] * self.duration))
            t_vec = np.arange(n_samples, dtype=float)
            if src.data_stats['srate'] > 0:
                t_vec /= src.data_stats['srate']

            row_offset += 1
            pw = self._widget.addPlot(row=row_offset, col=0, antialias=self._antialias)
            last_row = row_offset

            if self.show_chan_labels and not offset_chans:
                legend_bg = QtGui.QColor(self.bg_color)
                legend_bg.setAlphaF(0.5)
                legend = pg.LegendItem(offset=(0, 1), brush=legend_bg)
                legend.setParentItem(pw)
            else:
                legend = None

            pw.showGrid(x=True, y=True, alpha=0.3)
            font = QtGui.QFont()
            font.setPointSize(self.font_size - 2.0)
            pw.setXRange(0, self.duration)
            pw.getAxis("bottom").setTickFont(font)
            pw.getAxis("bottom").setStyle(showValues=self.ylabel_as_title)
            yax = pw.getAxis('left')
            yax.setTickFont(font)
            stream_ylabel = json.loads(src.identifier)['name']
            if 'unit' in ch_states and ch_states['unit'].nunique() == 1:
                # I don't use the `units=` kwarg here because it prepends a magnitude prefix (u, m, k, M),
                # which we don't always want.
                stream_ylabel = stream_ylabel + ' (%s)' % ch_states['unit'].iloc[0]
            pw.setLabel('top' if self.ylabel_as_title else 'left', self._ylabel or stream_ylabel, **labelStyle)
            if self.ylabel_as_title:
                pw.getAxis("top").setStyle(showValues=False)

            # Set the y-range
            if (n_vis_src <= 1 and self._requested_auto_scale != 'none') \
                    or (self._requested_auto_scale == 'all' and not offset_chans):
                # We rely on pg auto range
                pw.enableAutoRange(axis='y')
            elif self._auto_scale == 'none' and not offset_chans:
                # Set the y range manually to a constant.
                pw.setYRange(self.lower_limit, self.upper_limit)
            else:
                # using -0.5 to 0.5, per channel
                major_ticks = []
                minor_ticks = []
                if offset_chans:
                    pw.setYRange(-0.5, n_vis_src - 0.5)
                    chan_ticks = list(zip(range(n_vis_src), ch_states['name']))
                    if self._auto_scale == 'none':
                        data_ticks = [(-0.5, str(self.lower_limit)), (0.5, self.upper_limit)]
                        data_ticks += [(_ + 0.5, '') for _ in range(1, n_vis_src)]
                    else:
                        data_ticks = []
                    if self.show_chan_labels:
                        major_ticks = chan_ticks
                        minor_ticks = data_ticks
                    else:
                        major_ticks = data_ticks
                else:
                    # Channels are overlapping within (-0.5, 0.5), but the data are auto-scaled. ticks don't make sense.
                    pw.setYRange(-0.5, 0.5)
                yax.setTicks([major_ticks, minor_ticks])

            ch_offset_row = -1
            for ch_ix, ch_state in ch_states.iterrows():
                if ch_state['vis']:
                    ch_offset_row += 1
                    ch_offset_color = ch_offset_row if self.reset_colormap else (ch_offset_color + 1)
                    pen = pg.mkPen(color_map[ch_offset_color], width=self.line_width)
                    curve = pg.PlotCurveItem(t_vec, buff._data[0], connect='finite', pen=pen, name=ch_state['name'])
                    curve.setPos(0, ch_offset_row if offset_chans else 0)
                    pw.addItem(curve)
                    # pdi = pg.PlotDataItem(t_vec, buff._data[0],
                    #                       antialias=self.antialias,
                    #                       pen=pen, name=ch_state['name'])
                    # pw.addItem(pdi)
                    if self.show_chan_labels and not offset_chans:
                        legend.addItem(curve, name=ch_state['name'])

        # Unhide bottom x-axis. Link all axes to bottom axis.
        bottom_pw = self._widget.getItem(last_row, 0)
        bottom_pw.setLabel('bottom', 'Time', units='s', **labelStyle)
        bottom_pw.getAxis("bottom").setStyle(showValues=True)

        # Link for zooming
        for row_ix in range(last_row):
            pw = self._widget.getItem(row_ix, 0)
            pw.setXLink(bottom_pw)

        self._widget.ci.setSpacing(10. if self.ylabel_as_title else 0.)
        # self._widget.setContentsMargins(0., 0., 0., 0.)

        self._do_yaxis_sync = True

    def sync_y_axes(self):
        # Get max width
        max_width = self._ylabel_width or 0.
        for src_ix in range(len(self._data_sources)):
            pw = self._widget.getItem(src_ix, 0)
            if pw is None:
                break
            yax = pw.getAxis('left')
            max_width = max(max_width, yax.minimumWidth())

        # Apply widths
        for src_ix in range(len(self._data_sources)):
            pw = self._widget.getItem(src_ix, 0)
            if pw is None:
                break
            pw.getAxis('left').setWidth(max_width)

        self._do_yaxis_sync = False

    def update_visualization(self, data: np.ndarray, timestamps: np.ndarray) -> None:
        if not any([np.any(_) for _ in timestamps[0]]):
            return

        for src_ix in range(len(data)):
            pw = self._widget.getItem(src_ix, 0)
            if pw is None:
                # Can happen if reset_renderer is slow.
                return
            dat, mrk = data[src_ix]
            ts, mrk_ts = timestamps[src_ix]

            if dat.size:
                if self._do_yaxis_sync:
                    self.sync_y_axes()

                # Time series data or a continuous-ified version of irregular data

                per_chan_range = (-0.5, 0.5)  # channel range if data are auto-scaled.

                offset_chans = self.offset_channels and dat.shape[0] > 1
                if self.auto_scale != 'none':
                    # dat auto-scaled between (0, 1). Each channel is allocated (-0.5, 0.5)
                    dat -= 0.5
                elif offset_chans:
                    # Data are not auto-scaled, but they need to be scaled here because of the channel offsets.
                    coef = (per_chan_range[1] - per_chan_range[0]) / (self.upper_limit - self.lower_limit)
                    dat = dat - self.lower_limit
                    np.multiply(dat, coef, out=dat)
                    np.add(dat, per_chan_range[0], out=dat)

                for ch_ix, _d in enumerate(dat):
                    curve = pw.curves[ch_ix + (1 if self._antialias else 0)]
                    curve.setData(ts % self.duration, _d)

                    # curve.setData(ts % self.duration, _d)

            if not self._buffers[src_ix]._tvec.size:
                continue
            lead_t = self._buffers[src_ix]._tvec[self._buffers[src_ix]._write_idx]
            self._t_expired = max(lead_t - self._duration, self._t_expired)

            # Remove any markers that are more than 1 sweep older than any time series.
            if isinstance(pw.items[-1], pg.TextItem):
                while (len(self._marker_info) > 0) and (self._marker_info[0].timestamp < self._t_expired):
                    pop_info = self._marker_info.popleft()
                    pw.removeItem(pop_info[2])
                    self.marker_texts_pool.append(pop_info[2])
                    # And update the pw.curves with zeros at pop_info[0]
                    for pw_ix in range(1 if self._antialias else 0, len(pw.curves)):
                        old_x = pw.curves[pw_ix].xData
                        old_y = pw.curves[pw_ix].yData
                        rem_ix = np.searchsorted(old_x, pop_info[1] % self._duration)
                        if old_y[rem_ix] == np.min(old_y) and old_y[rem_ix + 1] == np.nanmax(old_y):
                            rem_ix = rem_ix + 1
                        old_y[rem_ix] = np.nan
                        pw.curves[pw_ix].setData(old_x, old_y)

            if mrk.size:
                b_new = mrk_ts > self._src_last_marker_time[src_ix]
                # y_offset = pw.viewRange()[1][0]
                for _t, _m in zip(mrk_ts[b_new], mrk[b_new]):
                    if len(self.marker_texts_pool) > 0:
                        text = self.marker_texts_pool.popleft()
                        text.setText(_m)
                    else:
                        text = pg.TextItem(text=_m, angle=90)
                        font = QtGui.QFont()
                        font.setPointSize(self.font_size + 2.0)
                        text.setFont(font)
                    text.setPos(_t % self.duration, -1)  # y_offset)
                    pw.addItem(text)
                    self._marker_info.append(MarkerMap(src_ix, _t, text))

                if np.any(b_new):
                    self._src_last_marker_time[src_ix] = mrk_ts[b_new][-1]

    @property
    def line_width(self):
        return self._line_width

    @line_width.setter
    def line_width(self, value):
        self._line_width = value
        self.reset_renderer(reset_channel_labels=False)

    @property
    def offset_channels(self):
        return self._offset_channels

    @offset_channels.setter
    def offset_channels(self, value):
        self._offset_channels = value
        self.reset_renderer(reset_channel_labels=True)

    @property
    def reset_colormap(self):
        return self._reset_colormap

    @reset_colormap.setter
    def reset_colormap(self, value):
        self._reset_colormap = value
        self.reset_renderer(reset_channel_labels=False)

    @property
    def antialias(self):
        return self._antialias

    @antialias.setter
    def antialias(self, value):
        self._antialias = value
        self.reset_renderer(reset_channel_labels=False)

    @property
    def ylabel_as_title(self):
        return self._ylabel_as_title

    @ylabel_as_title.setter
    def ylabel_as_title(self, value):
        self._ylabel_as_title = value
        self.reset_renderer(reset_channel_labels=True)

    @property
    def ylabel_width(self):
        return self._ylabel_width

    @ylabel_width.setter
    def ylabel_width(self, value):
        self._ylabel_width = value
        self.reset_renderer(reset_channel_labels=True)

    @RendererDataTimeSeries.auto_scale.setter
    def auto_scale(self, value):
        self._requested_auto_scale = value.lower()
        self.reset_renderer(reset_channel_labels=False)


text:
import argparse

import cv2

import numpy as np

import os

import subprocess

from tqdm import tqdm

if __name__ == '__main__':
    # Argument parsing
    parser = argparse.ArgumentParser(description='SURF feature extraction script')

    parser.add_argument(
        '--image_path', type=str, required=True,
        help='path to the images'
    )

    parser.add_argument(
        '--max_edge', type=int, default=1600,
        help='maximum image size at octave 0'
    )

    parser.add_argument(
        '--output_extension', type=str, default='.surf',
        help='extension for the output'
    )

    args = parser.parse_args()

    print(args)

    # Define SURF.
    surf = cv2.xfeatures2d.SURF_create()
    surf.setExtended(True)
    surf.setHessianThreshold(500)

    # Extract features.
    for image_name in os.listdir(args.image_path):
        image = cv2.imread(os.path.join(args.image_path, image_name))
        if image is None:
            continue
        image_size = image.shape

        # Resize such that max edge is 1600.
        image = image[:, :, ::-1]  # BGR -> RGB
        if len(image.shape) == 2:
            image = image[:, :, np.newaxis].repeat(3, axis=-1)
        downscaling_factor = max(1, max(image.shape) / args.max_edge)
        image = cv2.resize(image, None, fx=(1 / downscaling_factor), fy=(1 / downscaling_factor), interpolation=cv2.INTER_AREA)
        gray_image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

        # Extract SURF.
        opencv_keypoints, opencv_descriptors = surf.detectAndCompute(gray_image, None)

        # Convert to numpy.
        keypoints = []
        scores = []
        for keypoint in opencv_keypoints:
            keypoints.append([keypoint.pt[0], keypoint.pt[1], keypoint.size, keypoint.angle])
            scores.append(keypoint.response)
        keypoints = np.array(keypoints)
        if keypoints.shape[0] == 0:
            keypoints = np.zeros([0, 4])
        keypoints[:, : 2] *= downscaling_factor
        scores = np.array(scores)
        descriptors = opencv_descriptors
        if keypoints.shape[0] == 0:
            descriptors = np.zeros([0, 128])

        print('[%s] %dx%d, downscaling factor %.4f; %d keypoints' % (
            image_name, image_size[0], image_size[1], downscaling_factor, keypoints.shape[0]
        ))

        with open(os.path.join(args.image_path, image_name) + args.output_extension, 'wb') as output_file:
            np.savez(
                output_file,
                keypoints=keypoints,
                scores=scores,
                descriptors=descriptors
            )