text:
import torch
from torch.autograd import Variable
import numpy as np
from torchvision.utils import save_image
from training_function_generator_pytorch  import training_function_generator
from compute_gradient_penalty import compute_gradient_penalty
from sample_image import sample_image

#See how many of these can be removed
from data_iter import GenDataIter, DisDataIter
from target_lstm import TargetLSTM
from rollout import Rollout
import torch.nn as nn
from sequential_generate_sample import generate_samples
from sequential_eval_epoch import eval_epoch
from sequential_train_epoch import train_epoch

def train_GAN(conf_data):

    """Training Process for GAN.

    Parameters
    ----------
    conf_data: dict
        Dictionary containing all parameters and objects.

    Returns
    -------
    conf_data: dict
        Dictionary containing all parameters and objects.

    """
    seq = conf_data['GAN_model']['seq']
    if seq == 1:
        pre_epoch_num = conf_data['generator']['pre_epoch_num']
        GENERATED_NUM = 10000
        EVAL_FILE = 'eval.data'
        POSITIVE_FILE = 'real.data'
        NEGATIVE_FILE = 'gene.data'
    temp = 1 #TODO Determines how many times is the discriminator updated. Take this as a value input
    epochs = int(conf_data['GAN_model']['epochs'])
    if seq == 0:
        dataloader = conf_data['data_learn']
    mini_batch_size = int(conf_data['GAN_model']['mini_batch_size'])
    data_label = int(conf_data['GAN_model']['data_label'])
    cuda = conf_data['cuda']
    g_latent_dim = int(conf_data['generator']['latent_dim'])
    classes = int(conf_data['GAN_model']['classes'])

    w_loss = int(conf_data['GAN_model']['w_loss'])

    clip_value = float(conf_data['GAN_model']['clip_value'])
    n_critic = int(conf_data['GAN_model']['n_critic'])

    lambda_gp = int(conf_data['GAN_model']['lambda_gp'])

    log_file = open(conf_data['performance_log']+"/log.txt","w+")
    #Covert these to parameters of the config data
    Tensor = torch.cuda.FloatTensor if cuda else torch.FloatTensor
    conf_data['Tensor'] = Tensor
    LongTensor = torch.cuda.LongTensor if cuda else torch.LongTensor
    conf_data['LongTensor'] = LongTensor
    FloatTensor = torch.cuda.FloatTensor if cuda else torch.FloatTensor
    conf_data['FloatTensor'] = FloatTensor

    conf_data['epochs'] = epochs

    #print ("Just before training")
    if seq == 1: #TODO: Change back to 1
        target_lstm = TargetLSTM(conf_data['GAN_model']['vocab_size'], conf_data['generator']['embedding_dim'], conf_data['generator']['hidden_dim'], conf_data['cuda'])
        if cuda == True:
            target_lstm = target_lstm.cuda()
        conf_data['target_lstm'] = target_lstm
        gen_data_iter = GenDataIter('real.data', mini_batch_size)
        generator = conf_data['generator_model']
        discriminator = conf_data['discriminator_model']
        g_loss_func = conf_data['generator_loss']
        d_loss_func = conf_data['discriminator_loss']
        optimizer_D = conf_data['discriminator_optimizer']
        optimizer_G = conf_data['generator_optimizer']
        #print('Pretrain with MLE ...')
        for epoch in range(pre_epoch_num): #TODO: Change the range
            loss = train_epoch(generator, gen_data_iter, g_loss_func, optimizer_G,conf_data,'g')
            print('Epoch [%d] Model Loss: %f'% (epoch, loss))
            generate_samples(generator, mini_batch_size, GENERATED_NUM, EVAL_FILE,conf_data)
            eval_iter = GenDataIter(EVAL_FILE, mini_batch_size)
            loss = eval_epoch(target_lstm, eval_iter, g_loss_func,conf_data)
            print('Epoch [%d] True Loss: %f' % (epoch, loss))

        dis_criterion = d_loss_func
        dis_optimizer = optimizer_D
        #TODO: Understand why the below two code line were there ?
        # if conf_data['cuda']:
        #     dis_criterion = dis_criterion.cuda()

        #print('Pretrain Dsicriminator ...')
        dis_data_iter = DisDataIter(POSITIVE_FILE, NEGATIVE_FILE, mini_batch_size)
        for epoch in range(5): #TODO: change back 5
            generate_samples(generator, mini_batch_size, GENERATED_NUM, NEGATIVE_FILE,conf_data)
            dis_data_iter = DisDataIter(POSITIVE_FILE, NEGATIVE_FILE, mini_batch_size)
            for _ in range(3): #TODO: change back 3
                loss = train_epoch(discriminator, dis_data_iter, dis_criterion, dis_optimizer,conf_data,'d')
                print('Epoch [%d], loss: %f' % (epoch, loss))
        conf_data['generator_model'] = generator
        conf_data['discriminator_model'] = discriminator
        torch.save(conf_data['generator_model'].state_dict(),conf_data['save_model_path']+'/Seq/'+'pre_generator.pt')
        torch.save(conf_data['discriminator_model'].state_dict(),conf_data['save_model_path']+'/Seq/'+'pre_discriminator.pt')

        conf_data['rollout'] = Rollout(generator, 0.8)

    for epoch in range(epochs):
        conf_data['epoch'] = epoch
        if seq == 0:
            to_iter = dataloader
        elif seq == 1: #TODO: Change this back to 1
            to_iter = [1]

        for i, iterator in enumerate(to_iter):
            optimizer_D = conf_data['discriminator_optimizer']
            optimizer_G = conf_data['generator_optimizer']

            generator = conf_data['generator_model']
            discriminator = conf_data['discriminator_model']

            g_loss_func = conf_data['generator_loss']
            d_loss_func = conf_data['discriminator_loss']

            # if aux = 1:

            #print ("Reached here --------------> ")
            conf_data['iterator'] = i
            if seq == 0:

                if data_label == 1:
                    imgs,labels = iterator
                else:
                    imgs = iterator
                # Adversarial ground truths
                valid = Variable(Tensor(imgs.size(0), 1).fill_(1.0), requires_grad=False)
                conf_data['valid']=valid
                fake = Variable(Tensor(imgs.size(0), 1).fill_(0.0), requires_grad=False)
                conf_data['fake']=fake
                # Configure input
                real_imgs = Variable(imgs.type(Tensor))

                if data_label == 1:
                    labels = Variable(labels.type(LongTensor))
                # Sample noise as generator input
                z = Variable(Tensor(np.random.normal(0, 1, (imgs.shape[0], g_latent_dim))))
                if classes > 0:
                    gen_labels = Variable(LongTensor(np.random.randint(0, classes, imgs.shape[0])))
                    conf_data['gen_labels'] = gen_labels
            # elif seq == 1: #If yes seqGAN
            #     # samples = generator.sample(mini_batch_size,conf_data['generator']['sequece_length'])
            #     # zeros = torch.zeros((mini_batch_size,1)).type(LongTensor)
            #     # imgs = Variable(torch.cat([zeros,samples.data]),dim=1)[:,:-1].contiguous() #TODO: change imgs to inps all, to make more sense of the code
            #     # targets = Variable(sample.data).contiguous().view((-1,))
            #     # rewards = rollout.get_reward(sample,16,discriminator)
            #     # rewards = Variable(Tensor(rewards))
            #     # prob = generator.forward(inputs)
            #     # loss = gen_gan_loss(prob)
            #     pass
            #     #optimizer_G

            # ---------------------
            #  Train Discriminator
            # ---------------------
            optimizer_D.zero_grad()
            if seq == 1:#TODO change this back to 1
                dis_data_iter = DisDataIter(POSITIVE_FILE, NEGATIVE_FILE, mini_batch_size)
            for i in range(temp): # TODO: Make this a parameter -> for x updates --> I am read the stored models here as well. Should I reamove this ???
                optimizer_D = conf_data['discriminator_optimizer']
                optimizer_G = conf_data['generator_optimizer']

                generator = conf_data['generator_model']
                discriminator = conf_data['discriminator_model']

                g_loss_func = conf_data['generator_loss']
                d_loss_func = conf_data['discriminator_loss']
                if classes <= 0:
                    #print ("Reached here 2 --------------> ")
                    if seq == 0:
                        gen_imgs = generator(z)
                        # Measure discriminator's ability to classify real from generated samples
                        #Real images
                        real_validity = discriminator(real_imgs)
                        #Fake images
                        fake_validity = discriminator(gen_imgs.detach())
                    if seq == 1:
                        generate_samples(generator,mini_batch_size,GENERATED_NUM,NEGATIVE_FILE,conf_data)
                        dis_data_iter = DisDataIter(POSITIVE_FILE,NEGATIVE_FILE,mini_batch_size)
                        loss = train_epoch(discriminator,dis_data_iter,d_loss_func,optimizer_D,conf_data,'d')
                        conf_data['d_loss'] = loss
                        #exit()

                else:
                    if seq == 0:
                        gen_imgs = generator(z,gen_labels)
                        real_validity = discriminator(real_imgs, labels)
                        fake_validity = discriminator(gen_imgs.detach(), labels)

                if seq == 0:
                    conf_data['gen_imgs'] = gen_imgs
                if seq == 0:
                    if w_loss == 0:
                        real_loss = d_loss_func.loss(real_validity, valid)
                        fake_loss = d_loss_func.loss(fake_validity, fake)
                        d_loss = (real_loss + fake_loss) / 2
                    elif w_loss == 1:
                        d_loss = -d_loss_func.loss(real_validity,valid) + d_loss_func.loss(fake_validity,fake)
                        if lambda_gp > 0:
                            conf_data['real_data_sample'] = real_imgs.data
                            conf_data['fake_data_sample'] = gen_imgs.data
                            conf_data = compute_gradient_penalty(conf_data)
                            gradient_penalty = conf_data['gradient_penalty']
                            d_loss = d_loss+ lambda_gp * gradient_penalty
                    conf_data['d_loss'] = d_loss
                    d_loss.backward()
                    optimizer_D.step()

                if clip_value > 0:
                    # Clip weights of discriminator
                    for p in discriminator.parameters():
                        p.data.clamp_(-clip_value, clip_value)

            # -----------------
            #  Train Generator
            # -----------------
            conf_data['generator_model'] = generator
            conf_data['discriminator_model'] = discriminator

            #Next 4 lines were recently added maybe have to remove this.
            conf_data['optimizer_G'] = optimizer_G
            conf_data['optimizer_D'] = optimizer_D
            conf_data['generator_loss'] = g_loss_func
            conf_data['discriminator_loss'] = d_loss_func
            if seq == 0:
                conf_data['noise'] = z

            if n_critic <= 0:
                conf_data = training_function_generator(conf_data)
            elif n_critic > 0:
                # Train the generator every n_critic iterations
                if i % n_critic == 0:
                    conf_data = training_function_generator(conf_data)
            #exit()

           # print ("------------------ Here (train_GAN.py)")

            if seq == 0:
                batches_done = epoch * len(dataloader) + i
                if batches_done % int(conf_data['sample_interval']) == 0:
                    if classes <= 0:
                        # print ("Here")
                        # print (type(gen_imgs.data[:25]))
                        # print (gen_imgs.data[:25].shape)
                        save_image(gen_imgs.data[:25], conf_data['result_path']+'/%d.png' % batches_done, nrow=5, normalize=True)
                    elif classes > 0:
                        sample_image(10,batches_done,conf_data)
        if seq == 0:
            log_file.write("[Epoch %d/%d] [D loss: %f] [G loss: %f] \n" % (epoch, epochs,conf_data['d_loss'].item(), conf_data['g_loss'].item()))
        elif seq == 1:
            # print ("Done")
            log_file.write("[Epoch %d/%d] [D loss: %f] [G loss: %f] \n" % (epoch, epochs,conf_data['d_loss'], conf_data['g_loss']))
    conf_data['generator_model'] = generator
    conf_data['discriminator_model'] = discriminator
    conf_data['log_file'] = log_file
    return conf_data

max_stars_repo_name: Sinha-Raunak/gan-toolkit

max_issues_repo_name: Sinha-Raunak/gan-toolkit

max_forks_repo_name: Sinha-Raunak/gan-toolkit


text:
### This file is a part of the Syncpy library.
### Copyright 2015, ISIR / Universite Pierre et Marie Curie (UPMC)
### Main contributor(s): Giovanna Varni, Marie Avril,
### syncpy@isir.upmc.fr
###
### This software is a computer program whose for investigating
### synchrony in a fast and exhaustive way.
###
### This software is governed by the CeCILL-B license under French law
### and abiding by the rules of distribution of free software.  You
### can use, modify and/ or redistribute the software under the terms
### of the CeCILL-B license as circulated by CEA, CNRS and INRIA at the
### following URL "http://www.cecill.info".

### As a counterpart to the access to the source code and rights to
### copy, modify and redistribute granted by the license, users are
### provided only with a limited warranty and the software's author,
### the holder of the economic rights, and the successive licensors
### have only limited liability.
###
### In this respect, the user's attention is drawn to the risks
### associated with loading, using, modifying and/or developing or
### reproducing the software by the user in light of its specific
### status of free software, that may mean that it is complicated to
### manipulate, and that also therefore means that it is reserved for
### developers and experienced professionals having in-depth computer
### knowledge. Users are therefore encouraged to load and test the
### software's suitability as regards their requirements in conditions
### enabling the security of their systems and/or data to be ensured
### and, more generally, to use and operate it in the same conditions
### as regards security.
###
### The fact that you are presently reading this means that you have
### had knowledge of the CeCILL-B license and that you accept its terms.

"""
.. moduleauthor:: Marie Avril
"""

import numpy as np # For math operation
import pandas as pd # For DataFrame
import matplotlib.pyplot as plt # Plotting package
from math import ceil

from Methods.utils import Standardize
from Method import Method, MethodArgList

class S_Estimator(Method):
    """
    S_Estimator computes the S-Estimator, Genuine and Random Synchronization index among multiple monovariate signals (organized as a list of pandas DataFrame).

    **Reference :**

    * Cui, D. and al., Estimation of genuine and random synchronization in multivariate neural series. Neural Networks 23 (2010) 698-704.
    * Andrzejak, R. and al., Bivariate surrogate techniques: Necessity, strengths, and caveats. Physical Review E 68, 066202 (2003).
    * Schreiber, T. and al., Surrogate time series. Physica D 142 (2000) 346-382.

    :param surr_nb_iter:
        Number of surrogate iterations. Default : 100
    :type surr_nb_iter: int

    :param plot:
        if True the plot of surrogates signals is returned. Default: False
    :type plot: bool
    """
    argsList = MethodArgList()
    argsList.append('surr_nb_iter', 100, int,
                    'Number of surrogate iterations.')
    argsList.append('plot', False, bool, 'plot of surrogates signals is returned')

    ''' Constructor '''
    def __init__(self, surr_nb_iter = 100, plot = False, **kwargs):
        super(S_Estimator, self).__init__(plot, **kwargs)
        #In the constructor we can check that params have corrects values and initialize stuff

        ' Raise error if parameters are not in the correct type '
        try :
            if not(isinstance(surr_nb_iter, int)) : raise TypeError("Requires surr_nb_iter to be an integer")
            if not(isinstance(plot, bool))        : raise TypeError("Requires plot to be a boolean")
        except TypeError, err_msg:
            raise TypeError(err_msg)
            return

        ' Raise error if parameters do not respect input rules '
        try :
            if surr_nb_iter <= 0 : raise ValueError("Requires surr_nb_iter to be a strictly positive scalar")
        except ValueError, err_msg:
            raise ValueError(err_msg)
            return

        self._surr_nb_iter = surr_nb_iter
        self._plot = plot

        return

    def plot_result(self, result):
        """
        Plot surrogates signals

        :param result:
            S-estimator result from compute() (only 'surrogate_signal' used)
        :type result: dict

        :returns: plt.figure
         -- figure plot

        """
        ' Raise error if parameters are not in the correct type '
        try :
            if not(isinstance(result, dict)) : raise TypeError("Requires result to be a dictionary")
        except TypeError, err_msg:
            raise TypeError(err_msg)
            return

        ' Raise error if not the good dictionary '
        try :
            if not 'surrogate_signal' in result : raise ValueError("Requires dictionary to be the output of compute() method")
        except ValueError, err_msg:
            raise ValueError(err_msg)
            return

        nrows = int(ceil(result['surrogate_signal'].shape[1]))
        figure, ax = plt.subplots(nrows, sharex=True)

        n = np.arange(0,result['surrogate_signal'].shape[0])
        idx = 0
        for col in result['surrogate_signal'].columns :
            ax[idx].grid(True) # Display a grid
            ax[idx].set_title('Surrogate signal for : ' + col + ' variable')
            ax[idx].plot(n, result['surrogate_signal'][col].values)
            idx += 1

        ax[idx-1].set_xlabel('Samples')

        return figure

    def getSynchronizationIndex(self, lambda_i):
        """
         Compute Synchronization Index (SI)

        :param lambda_i:
            normalized eigenvalues (depending on the type of SI)
        :type lambda_i: np.array

        :returns: float
            -- Synchronization index
        """

        tmp = [lambda_i[k] * np.log(lambda_i[k]) for k in range(len(lambda_i))]

        SI = 1 + sum(tmp) / np.log(len(lambda_i))

        return SI

    def AAFT_surrogates(self,Xi):
        """
         Computes amplitude adjusted Fourier Transform (AAFT) method to create a surrogate signal.
         Get starting point / initial conditions for R surrogates

        :param Xi:
            signal to surrogate
        :type Xi: np.array

        :returns: np.array
            -- surrogated signal
        """

        Xi_fft = np.fft.rfft(Xi)
        Xi_ampli = np.abs(Xi_fft) # get original amplitude

        #  Create sorted Gaussian reference series
        gaussian = np.random.randn(Xi.shape[0])
        gaussian.sort()

        #  Rescale data to Gaussian distribution
        ranks = Xi.argsort()
        rescaled_data = np.zeros(Xi.shape[0])

        for i in xrange(Xi.shape[0]):
            rescaled_data[i] = gaussian[ranks[i]]

        #  Phase randomize rescaled data

        #  Get shapes
        N = Xi.shape[0]
        len_phase = Xi_fft.shape

        #  Generate random phases uniformly distributed in the
        #  interval [0, 2*Pi]
        phases = np.random.uniform(low=0, high=2 * np.pi)

        #  Add random phases uniformly distributed in the interval [0, 2*Pi]
        Xi_fft *= np.exp(1j * phases)

        #  Calculate IFFT and take the real part, the remaining imaginary part
        #  is due to numerical errors.
        phase_randomized_data = \
            np.ascontiguousarray(np.real(np.fft.irfft(Xi_fft)))

        #  Rescale back to amplitude distribution of original data
        sorted_original = Xi.copy()
        sorted_original.sort()

        ranks = phase_randomized_data.argsort()

        for i in xrange(Xi.shape[0]):
            rescaled_data[i] = sorted_original[ranks[i]]

        return rescaled_data

    def refined_AAFT_surrogate(self, X):
        """
         Computes an Iteratively refined amplitude adjusted Fourier Transform (AAFT) method to create a surrogate signal.

        .. note:: Assume that original signal X is already standardized.

        :param X:
            signal to surrogate
        :type X: pd.DataFrame

        :returns: pd.DataFrame
            -- surrogated signal
        :returns: np.array
            -- average eigvalues of surrogate signal among all iterations
        """

        ''' 1- Initialization '''
        X_fft = {}
        X_ampli = {}
        X_sorted = {}
        X_surr = {}
        X_surr_eig = {}
        X_surr_tot_ampli = {}
        for col in X.columns :
            # get real amplitude
            Xi_fft = np.fft.rfft(X[col].values)
            Xi_ampli = np.abs(Xi_fft)
            X_fft.update({col : Xi_fft})
            X_ampli.update({col : Xi_ampli})

            # get sorted real values
            sorted_Xi = X.copy()[col].values
            sorted_Xi.sort()
            X_sorted.update({col : sorted_Xi})

            # shuffle original data for surrogate initialization
            Xi_surr = X.copy()[col].values
            np.random.shuffle(Xi_surr)
            #Xi_surr = self.AAFT_surrogates(Xi)
            X_surr.update({col : Xi_surr})

        # save initial surrogate eigen values
        df_X_surr = pd.DataFrame(X_surr)
        R = np.dot((df_X_surr.values).T, df_X_surr.values) / float(X.index.size) # Correlation matrix
        eig_values_surr, eig_vectors_surr = np.linalg.eig(R) # Eigenvalues decomposition
        eig_values_surr = np.sort(np.real(eig_values_surr)) # sort by increasing order
        X_surr_eig.update({0 : eig_values_surr})

        for p in range(self._surr_nb_iter) :
            X_surr_ampli = {}
            for col in X_surr.keys() :
                Xi_surr = X_surr[col]

                ''' 2 - Filtering '''
                Xi_surr_fft = np.fft.rfft(Xi_surr)
                Xi_surr_ampli = np.abs(Xi_surr_fft)
                Xi_surr_phase = Xi_surr_fft / Xi_surr_ampli

                # Save surrogate amplitude
                X_surr_ampli.update({col : Xi_surr_ampli})

                # inverse fft with real amplitude
                s = np.fft.irfft(X_ampli[col] * Xi_surr_phase)

                ''' 3 - Rescaling '''
                ranks = s.argsort().argsort() # get sorted indexes
                for j in xrange(Xi_surr.shape[0]):
                    Xi_surr[j] = X_sorted[col][ranks[j]] # replace smallest Xi_surr by smallets Xi etc...
                X_surr.update({col : Xi_surr})

            ''' 4 - Save surrogates eigen values '''
            df_X_surr = pd.DataFrame(X_surr)
            R = np.dot((df_X_surr.values).T, df_X_surr.values) / float(X.index.size) # Correlation matrix
            eig_values_surr, eig_vectors_surr = np.linalg.eig(R) # Eigenvalues decomposition
            eig_values_surr = np.sort(eig_values_surr) # sort by increasing order
            X_surr_eig.update({(p+1) : eig_values_surr})

            # Save surrogate amplitude at each iteration to see if there is a convergence to real amplitude
            X_surr_tot_ampli.update({p : X_surr_ampli})

        # Compute average surrogate eigen values
        X_surr_average_eig = np.mean(X_surr_eig.values(), axis=0)

        return df_X_surr, X_surr_average_eig

    def compute(self, signals):
        """
         Computes SSI for multiple monovariate signals (organized as a list).
         If input signals are multivariates, only the first column of the signal is considered

        :param signals:
            list of signals, one per person.
        :type signals: list[pd.DataFrame]

        :returns: dict
            -- Synchronization indexes : S-Estimator (SSI), Genuine Synchronization Index (GSI) and Random Synchronization Index (RSI)
        """

        ' Raise error if parameters are not in the correct type '
        try :
            for i in range(len(signals)) :
                if not(isinstance(signals[i], pd.DataFrame)): raise TypeError("Requires signal " + str(i+1) + " to be a pd.DataFrame.")
        except TypeError, err_msg:
            raise TypeError(err_msg)
            return

        ' Raise error if DataFrames have not the same size or same indexes '
        try :
            for i in range(0,len(signals)):
                if len(signals[0]) != len(signals[i]) : raise ValueError("All the signals must have the same size. Signal " + str(i+1) + " does not have the same size as first signal.")
                if signals[0].index.tolist() != signals[i].index.tolist() : raise ValueError("All the signals must have the same time indexes. Signal " + str(i+1) + " does not have the same time index as first signal.")
        except ValueError, err_msg:
            raise ValueError(err_msg)
            return

        'Formate signals in one DataFrame for computing'
        # If input signals are multivariates, only the first column is considered
        x = pd.DataFrame()
        for i in range(0,len(signals)):
            if x.empty :
                x = pd.DataFrame(signals[i].iloc[:,0], signals[i].index)
                x.columns = [signals[i].columns[0]]
            else :
                x[signals[i].columns[0]] = signals[i].iloc[:,0]

        ''' Ignore last value if len(x) is odd (avoiding trouble with fft<=>ifft)'''
        if (x.shape[0]%2 != 0):
            x = x.iloc[0:x.shape[0]-1,:]

        ''' Standardize '''
        X = Standardize.Standardize(x)

        ''' Correlation matrix  '''
        C = np.dot((X.values).T, X.values) / float(x.index.size) # on original signal

        ''' Eigenvalues decomposition '''
        eig_values, eig_vectors = np.linalg.eig(C) # on original signal
        eig_values = np.sort(eig_values) # sort by increasing order

        ''' Surrogate signal '''
        df_X_surr, X_surr_average_eig = self.refined_AAFT_surrogate(X)

        ''' Get Synchronization Indexes '''
        lambda_1 = []
        lambda_2 = []
        lambda_3 = []
        for i in range(X_surr_average_eig.size):
            lambda_1.append( eig_values[i] / sum(eig_values) )

            lambda_2.append( (eig_values[i]/X_surr_average_eig[i]) / sum(eig_values/X_surr_average_eig) )

            lambda_3.append( X_surr_average_eig[i] / sum(X_surr_average_eig) )

        SI = dict()
        SI['SSI'] = self.getSynchronizationIndex(lambda_1)
        SI['GSI'] = self.getSynchronizationIndex(lambda_2)
        SI['RSI'] = self.getSynchronizationIndex(lambda_3)
        SI['surrogate_signal'] = df_X_surr

        if self._plot:
            self.plot_result(SI)

        return SI

    @staticmethod
    def getArguments():
        return S_Estimator.argsList.getMethodArgs()

    @staticmethod
    def getArgumentsAsDictionary():
        return S_Estimator.argsList.getArgumentsAsDictionary()

ext: py

max_stars_repo_name: syncpy/SyncPy

max_issues_repo_name: syncpy/SyncPy

max_forks_repo_name: syncpy/SyncPy

reason: import numpy


text:
import cv2
import numpy as np

# Build Face Detection Function to Clean Up Code
def faceDetect(net, frame):

    # DNN Model parameters
    in_width = 300
    in_height = 300
    mean = [104, 117, 123]
    conf_threshold = 0.7

    # Frame Parameters
    frame_height = frame.shape[0]
    frame_width = frame.shape[1]

    # Create a 4D blob from a frame.
    blob = cv2.dnn.blobFromImage(frame, 1.0, (in_width, in_height), mean, swapRB = False, crop = False)

    # Run the Face Detection Model provided by OpenCV
    net.setInput(blob)
    detections = net.forward()

    for i in range(detections.shape[2]):
        confidence = detections[0, 0, i, 2]
        if confidence > conf_threshold:
            x_left_bottom = int(detections[0, 0, i, 3] * frame_width)
            y_left_bottom = int(detections[0, 0, i, 4] * frame_height)
            x_right_top = int(detections[0, 0, i, 5] * frame_width)
            y_right_top = int(detections[0, 0, i, 6] * frame_height)

            cv2.rectangle(frame, (x_left_bottom, y_left_bottom), (x_right_top, y_right_top), (0, 255, 0))
            label = "Confidence: %.4f" % confidence
            label_size, base_line = cv2.getTextSize(label, cv2.FONT_HERSHEY_SIMPLEX, 0.5, 1)

            cv2.rectangle(frame, (x_left_bottom, y_left_bottom - label_size[1]),
                                (x_left_bottom + label_size[0], y_left_bottom + base_line),
                                (255, 255, 255), cv2.FILLED)
            cv2.putText(frame, label, (x_left_bottom, y_left_bottom),
                        cv2.FONT_HERSHEY_SIMPLEX, 0.5, (0, 0, 0))

            # t, _ = net.getPerfProfile()
            # label = 'Inference time: %.2f ms' % (t * 1000.0 / cv2.getTickFrequency())
            center = [(x_right_top - x_left_bottom)/2 + x_left_bottom, \
                      (y_right_top - y_left_bottom)/2 + y_left_bottom]
            center[0] = int(center[0])
            center[1] = int(center[1])
            cv2.circle(frame, (center[0], center[1]), 10, (255, 0, 255), 2)
            label1 = 'Face Center at: %s px across, ' % center[0]
            label2 = '%s px down' % center[1]
            label = label1 + label2
            cv2.putText(frame, label, (0, 15), cv2.FONT_HERSHEY_SIMPLEX, 0.5, (0, 255, 0))
            DETECT_FLAG = True

        else:
            center = [int(frame.shape[0]/2), int(frame.shape[1]/2)]
            DETECT_FLAG = False

    # Return new frame and face detection coordinates
    u_face = center[0]
    v_face = center[1]
    return frame, DETECT_FLAG, u_face, v_face

#-------------------------------------------------------------------------------------------------------------
# Low-Grade "Control" System

def centerFace(frame, u_face, v_face):

    # Frame Parameters
    frame_height = frame.shape[0]
    frame_width = frame.shape[1]
    u_center = frame_height / 2
    v_center = frame_width / 2

    # Determine Error From Center
    du = u_face - u_center
    dv = v_face - v_center

    # Determine Az Adjustment (1 degree adjustments for now)
    if (np.abs(dv) <= 20):
        az_adj = 0
    elif (dv > 0):
        az_adj = -1
    elif (dv < 0):
        az_adj = 1
    else:
        az_adj = 0

    # Determine El Adjustment (1 degree adjustments for now)
    if (np.abs(du) <= 20):
        el_adj = 0
    elif (du > 0):
        el_adj = -1
    elif (du < 0):
        el_adj = 1
    else:
        el_adj = 0

    # Return Function Outputs
    return az_adj, el_adj


text:
import cv2, torch, torch.nn as nn
import numpy as np
import kornia as k, os, glob
from data.augmentation import Augmentation_random_crop

class Dataset(torch.utils.data.Dataset):

    def __init__(self, opt, mode):

        self.mode = 'train' if mode == 'train' else 'eval'

        self.seq_length = opt.Data['sequence_length']
        self.do_aug = opt.Data['aug']
        self.data_path = opt.Data['data_path'] + opt.Data['texture'] + '/' + self.mode + '/'
        self.img_format = opt.Data['image_format']

        self.videos = []; self.num_frames = []
        video_list = os.listdir(self.data_path)
        for vid in video_list:
            length = len(glob.glob(self.data_path + vid + f'/*.{self.img_format}'))
            for _ in range(opt.Data[f'iter_{mode}']):
                self.videos.append(vid)
                self.num_frames.append(length)

        self.length = len(self.videos)
        if self.mode == 'train' and self.do_aug:
            self.aug = Augmentation_random_crop(opt.Data['img_size'], opt.Data.Augmentation)
        else:
            self.aug = torch.nn.Sequential(
                        k.Resize(size=(opt.Data['img_size'], opt.Data['img_size'])),
                        k.augmentation.Normalize(0.5, 0.5))

    def __len__(self):
        return self.length

    def load_img(self, video, frame):
        img = cv2.imread(self.data_path + video + f'/frame{frame}.{self.img_format}')
        return k.image_to_tensor(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))/255.0

    def __getitem__(self, idx):
        video  = self.videos[idx]
        frames = np.arange(0, self.num_frames[idx])

        ## Sample random starting point in the sequence
        start_rand = np.random.randint(0, len(frames) - self.seq_length + 1) if self.mode == 'train' else 0
        seq = torch.stack([self.load_img(video, frames[start_rand + i]) for i in range(self.seq_length)], dim=0)
        return {'seq': self.aug(seq)}


text:
"""
Meta-parametrized SRF layer with alpha and sigma linearly dependent on t.
"""
# Import general dependencies
import numpy as np
import torch
import torch.nn as nn
from srf.gaussian_basis_filters import gaussian_basis_filters_shared
import torch.nn.functional as F

class Srf_layer_alpha_sigma_lin(nn.Module):
    def __init__(self,
                inC,
                outC,
                init_k,
                init_order,
                init_scale,
                learn_sigma,
                use_cuda):
        super(Srf_layer_alpha_sigma_lin, self).__init__()

        self.init_k = init_k
        self.init_order = init_order
        self.init_scale = init_scale
        self.inC = inC
        self.outC = outC
        #---------------
        F = int((self.init_order + 1) * (self.init_order + 2) / 2)

        """ Create weight variables. """
        self.use_cuda = use_cuda
        self.device = torch.device("cuda" if use_cuda else "cpu")

        self.alphas = torch.nn.Parameter(torch.zeros([F, inC, outC], \
                            device=self.device), requires_grad=True)
        torch.nn.init.normal_(self.alphas, mean=0.0, std=0.1)

        self.alphas_bias = torch.nn.Parameter(torch.zeros([F, inC, outC], \
                            device=self.device), requires_grad=True)
        torch.nn.init.normal_(self.alphas_bias, mean=0.0, std=0.05)

        if learn_sigma:
            self.scales = torch.nn.Parameter(torch.zeros([1],\
                            device=self.device), requires_grad=True)
            torch.nn.init.normal_(self.scales, mean=-0.0, std=2./3.)

            self.scale_bias = torch.nn.Parameter(torch.zeros([1],\
                            device=self.device), requires_grad=True)
            torch.nn.init.normal_(self.scale_bias, mean=-0.0, std=1./10.)
        else:
            self.scales = torch.nn.Parameter(torch.tensor(np.full((1), \
                            self.init_scale), device=self.device,\
                            dtype=torch.float32), requires_grad=False)
        self.extra_reg = 0

    def forward(self, data, t):
        # Define sigma from the scale
        self.sigma = 2.0**(self.scales*t+self.scale_bias)
        filters, _ = gaussian_basis_filters_shared(
                                            order=self.init_order, \
                                            sigma=self.sigma, \
                                            k=self.init_k, \
                                            alphas=self.alphas*t+self.alphas_bias,\
                                            use_cuda=self.use_cuda)
        final_conv = F.conv2d(
                        input=data, # NCHW
                        weight=filters, # KCHW
                        bias=None,
                        stride=1,
                        padding=int(filters.shape[2]/2))

        return final_conv


text:
import numpy as np
import cv2
import glob
import pickle
import matplotlib.pyplot as plt

def main():
    # Load calibration images
    images = glob.glob('camera_cal/*.jpg')
    mtx, dist = calib_cam(images, num_checker_pts=(9, 6))

    # Save the camera calibration result for later use
    dist_pickle = {}
    dist_pickle["mtx"] = mtx
    dist_pickle["dist"] = dist
    pickle.dump(dist_pickle, open("calib_cam_dist_pickle.p", "wb" ) )

    # Show an example of undistortion
    show_test_undist(mtx, dist,
                     src_path='camera_cal/calibration3.jpg', dst_path='output_images/test_undist.jpg')
    return

def calib_cam(image_paths, num_checker_pts):
    """ Do camera calibration given object points and image points
    https://github.com/udacity/CarND-Camera-Calibration
    """

    # prepare object points, like (0,0,0), (1,0,0), (2,0,0) ....,(6,5,0)
    n_px = num_checker_pts[0]
    n_py = num_checker_pts[1]
    objp = np.zeros((n_px * n_py, 3), np.float32)
    objp[:, :2] = np.mgrid[0:n_px, 0:n_py].T.reshape(-1, 2)

    # Arrays to store object points and image points from all the images.
    objpoints = []  # 3d points in real world space
    imgpoints = []  # 2d points in image plane.

    # Step through the list and search for chessboard corners
    for idx, fname in enumerate(image_paths):
        img = cv2.imread(fname)
        gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

        # Find the chessboard corners
        ret, corners = cv2.findChessboardCorners(gray, (n_px, n_py), None)

        # If found, add object points, image points
        if ret:
            objpoints.append(objp)
            imgpoints.append(corners)

    img_size = (img.shape[1], img.shape[0])
    ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(objpoints, imgpoints, img_size, None, None)
    return mtx, dist

def show_test_undist(mtx, dist, src_path, dst_path):
    img = cv2.imread(src_path)
    img_size = (img.shape[1], img.shape[0])
    dst = cv2.undistort(img, mtx, dist, None, mtx)
    cv2.imwrite(dst_path, dst)
    f, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 10))
    ax1.imshow(img)
    ax1.set_title('Original Image', fontsize=30)
    ax2.imshow(dst)
    ax2.set_title('Undistorted Image', fontsize=30)
    plt.show()
    return

if __name__ == '__main__':
    main()


text:
''' This module contains experiments to handle data '''

import os
import numpy as np
import scipy.ndimage
from PIL import Image
import scipy
import sys
from matplotlib import pyplot as plt

class MnistHandler(object):

    ''' Provides a convenient interface to manipulate MNIST data '''

    def __init__(self):

        # Download data if needed
        self.X_train, self.y_train, self.X_val, self.y_val, self.X_test, self.y_test = self.load_dataset()

        # Load Lena image to memory
        self.lena = Image.open('resources/lena.jpg')

    def load_dataset(self):
        # Credit for this function: https://github.com/Lasagne/Lasagne/blob/master/examples/mnist.py

        # We first define a download function, supporting both Python 2 and 3.
        if sys.version_info[0] == 2:
            from urllib import urlretrieve
        else:
            from urllib.request import urlretrieve

        def download(filename, source='http://yann.lecun.com/exdb/mnist/'):
            print("Downloading %s" % filename)
            urlretrieve(source + filename, filename)

        # We then define functions for loading MNIST images and labels.
        # For convenience, they also download the requested files if needed.
        import gzip

        def load_mnist_images(filename):
            if not os.path.exists(filename):
                download(filename)
            # Read the inputs in Yann LeCun's binary format.
            with gzip.open(filename, 'rb') as f:
                data = np.frombuffer(f.read(), np.uint8, offset=16)
            # The inputs are vectors now, we reshape them to monochrome 2D images,
            # following the shape convention: (examples, channels, rows, columns)
            data = data.reshape(-1, 1, 28, 28)
            # The inputs come as bytes, we convert them to float32 in range [0,1].
            # (Actually to range [0, 255/256], for compatibility to the version
            # provided at http://deeplearning.net/data/mnist/mnist.pkl.gz.)
            return data / np.float32(256)

        def load_mnist_labels(filename):
            if not os.path.exists(filename):
                download(filename)
            # Read the labels in Yann LeCun's binary format.
            with gzip.open(filename, 'rb') as f:
                data = np.frombuffer(f.read(), np.uint8, offset=8)
            # The labels are vectors of integers now, that's exactly what we want.
            return data

        # We can now download and read the training and test set images and labels.
        X_train = load_mnist_images('resources/train-images-idx3-ubyte.gz')
        y_train = load_mnist_labels('resources/train-labels-idx1-ubyte.gz')
        X_test = load_mnist_images('resources/t10k-images-idx3-ubyte.gz')
        y_test = load_mnist_labels('resources/t10k-labels-idx1-ubyte.gz')

        # We reserve the last 10000 training examples for validation.
        X_train, X_val = X_train[:-10000], X_train[-10000:]
        y_train, y_val = y_train[:-10000], y_train[-10000:]

        # We just return all the arrays in order, as expected in main().
        # (It doesn't matter how we do this as long as we can read them again.)
        return X_train, y_train, X_val, y_val, X_test, y_test

    def process_batch(self, batch, batch_size, image_size=28, color=False, rescale=True):

        # Resize from 28x28 to 64x64
        if image_size == 64:
            batch_resized = []
            for i in range(batch.shape[0]):
                # resize to 64x64 pixels
                batch_resized.append(scipy.ndimage.zoom(batch[i, :, :], 2.3, order=1))
            batch = np.stack(batch_resized)

        # Convert to RGB
        batch = batch.reshape((batch_size, 1, image_size, image_size))
        batch = np.concatenate([batch, batch, batch], axis=1)

        # Modify images if color distribution requested
        if color:

            # Binarize images
            batch[batch >= 0.5] = 1
            batch[batch < 0.5] = 0

            # For each image in the mini batch
            for i in range(batch_size):

                # Take a random crop of the Lena image (background)
                x_c = np.random.randint(0, self.lena.size[0] - image_size)
                y_c = np.random.randint(0, self.lena.size[1] - image_size)
                image = self.lena.crop((x_c, y_c, x_c + image_size, y_c + image_size))
                image = np.asarray(image).transpose((2, 0, 1)) / 255.0

                # Randomly alter the color distribution of the crop
                for j in range(3):
                    image[j, :, :] = (image[j, :, :] + np.random.uniform(0, 1)) / 2.0

                # Invert the color of pixels where there is a number
                image[batch[i, :, :, :] == 1] = 1 - image[batch[i, :, :, :] == 1]
                batch[i, :, :, :] = image

        # Rescale to range [-1, +1]
        if rescale:
            batch = batch * 2 - 1

        # Channel last
        batch = batch.transpose((0, 2, 3, 1))

        return batch

    def get_batch(self, subset, batch_size, image_size=28, color=False, rescale=True):

        # Select a subset
        if subset == 'train':
            X = self.X_train
            y = self.y_train
        elif subset == 'valid':
            X = self.X_val
            y = self.y_val
        elif subset == 'test':
            X = self.X_test
            y = self.y_test

        # Random choice of samples
        idx = np.random.choice(X.shape[0], batch_size)
        batch = X[idx, 0, :].reshape((batch_size, 28, 28))

        # Process batch
        batch = self.process_batch(batch, batch_size, image_size, color, rescale)

        # Image label
        labels = y[idx]

        return batch.astype('float32'), labels.astype('int32')

    def get_batch_by_labels(self, subset, labels, image_size=28, color=False, rescale=True):

        # Select a subset
        if subset == 'train':
            X = self.X_train
            y = self.y_train
        elif subset == 'valid':
            X = self.X_val
            y = self.y_val
        elif subset == 'test':
            X = self.X_test
            y = self.y_test

        # Find samples matching labels
        idxs = []
        for i, label in enumerate(labels):

            idx = np.where(y == label)[0]
            idx_sel = np.random.choice(idx, 1)[0]
            idxs.append(idx_sel)

        # Retrieve images
        batch = X[np.array(idxs), 0, :].reshape((len(labels), 28, 28))

        # Process batch
        batch = self.process_batch(batch, len(labels), image_size, color, rescale)

        return batch.astype('float32'), labels.astype('int32')

    def get_n_samples(self, subset):

        if subset == 'train':
            y_len = self.y_train.shape[0]
        elif subset == 'valid':
            y_len = self.y_val.shape[0]
        elif subset == 'test':
            y_len = self.y_test.shape[0]

        return y_len

class MnistGenerator(object):

    ''' Data generator providing MNIST data '''

    def __init__(self, batch_size, subset, image_size=28, color=False, rescale=True):

        # Set params
        self.batch_size = batch_size
        self.subset = subset
        self.image_size = image_size
        self.color = color
        self.rescale = rescale

        # Initialize MNIST dataset
        self.mnist_handler = MnistHandler()
        self.n_samples = self.mnist_handler.get_n_samples(subset)
        self.n_batches = self.n_samples // batch_size

    def __iter__(self):
        return self

    def __next__(self):
        return self.next()

    def __len__(self):
        return self.n_batches

    def next(self):

        # Get data
        x, y = self.mnist_handler.get_batch(self.subset, self.batch_size, self.image_size, self.color, self.rescale)

        # Convert y to one-hot
        y_h = np.eye(10)[y]

        return x, y_h

class SortedNumberGenerator(object):

    ''' Data generator providing lists of sorted numbers '''

    def __init__(self, batch_size, subset, terms, positive_samples=1, predict_terms=1, image_size=28, color=False, rescale=True):

        # Set params
        self.positive_samples = positive_samples
        self.predict_terms = predict_terms
        self.batch_size = batch_size
        self.subset = subset
        self.terms = terms
        self.image_size = image_size
        self.color = color
        self.rescale = rescale

        # Initialize MNIST dataset
        self.mnist_handler = MnistHandler()
        self.n_samples = self.mnist_handler.get_n_samples(subset) // terms
        self.n_batches = self.n_samples // batch_size

    def __iter__(self):
        return self

    def __next__(self):
        return self.next()

    def __len__(self):
        return self.n_batches

    def next(self):

        # Build sentences
        image_labels = np.zeros((self.batch_size, self.terms + self.predict_terms))
        sentence_labels = np.ones((self.batch_size, 1)).astype('int32')
        positive_samples_n = self.positive_samples
        for b in range(self.batch_size):

            # Set ordered predictions for positive samples
            seed = np.random.randint(0, 10)
            sentence = np.mod(np.arange(seed, seed + self.terms + self.predict_terms), 10)

            if positive_samples_n <= 0:

                # Set random predictions for negative samples
                # Each predicted term draws a number from a distribution that excludes itself
                numbers = np.arange(0, 10)
                predicted_terms = sentence[-self.predict_terms:]
                for i, p in enumerate(predicted_terms):
                    predicted_terms[i] = np.random.choice(numbers[numbers != p], 1)
                sentence[-self.predict_terms:] = np.mod(predicted_terms, 10)
                sentence_labels[b, :] = 0

            # Save sentence
            image_labels[b, :] = sentence

            positive_samples_n -= 1

        # Retrieve actual images
        images, _ = self.mnist_handler.get_batch_by_labels(self.subset, image_labels.flatten(), self.image_size, self.color, self.rescale)

        # Assemble batch
        images = images.reshape((self.batch_size, self.terms + self.predict_terms, images.shape[1], images.shape[2], images.shape[3]))
        x_images = images[:, :-self.predict_terms, ...]
        y_images = images[:, -self.predict_terms:, ...]

        # Randomize
        idxs = np.random.choice(sentence_labels.shape[0], sentence_labels.shape[0], replace=False)

        return [x_images[idxs, ...], y_images[idxs, ...]], sentence_labels[idxs, ...]

class SameNumberGenerator(object):

    ''' Data generator providing lists of similar numbers '''

    def __init__(self, batch_size, subset, terms, positive_samples=1, predict_terms=1, image_size=28, color=False, rescale=True):

        # Set params
        self.positive_samples = positive_samples
        self.predict_terms = predict_terms
        self.batch_size = batch_size
        self.subset = subset
        self.terms = terms
        self.image_size = image_size
        self.color = color
        self.rescale = rescale

        # Initialize MNIST dataset
        self.mnist_handler = MnistHandler()
        self.n_samples = self.mnist_handler.get_n_samples(subset) // terms
        self.n_batches = self.n_samples // batch_size

    def __iter__(self):
        return self

    def __next__(self):
        return self.next()

    def __len__(self):
        return self.n_batches

    def next(self):

        # Build sentences
        image_labels = np.zeros((self.batch_size, self.terms + self.predict_terms))
        sentence_labels = np.ones((self.batch_size, 1)).astype('int32')
        positive_samples_n = self.positive_samples
        for b in range(self.batch_size):

            # Set positive samples
            seed = np.random.randint(0, 10)
            sentence = seed * np.ones(self.terms + self.predict_terms)

            if positive_samples_n <= 0:

                # Set random predictions for negative samples
                sentence[-self.predict_terms:] = np.mod(sentence[-self.predict_terms:] + np.random.randint(1, 10, self.predict_terms), 10)
                sentence_labels[b, :] = 0

            # Save sentence
            image_labels[b, :] = sentence

            positive_samples_n -= 1

        # Retrieve actual images
        images, _ = self.mnist_handler.get_batch_by_labels(self.subset, image_labels.flatten(), self.image_size, self.color, self.rescale)

        # Assemble batch
        images = images.reshape((self.batch_size, self.terms + self.predict_terms, images.shape[1], images.shape[2], images.shape[3]))
        x_images = images[:, :-self.predict_terms, ...]
        y_images = images[:, -self.predict_terms:, ...]

        # Randomize
        idxs = np.random.choice(sentence_labels.shape[0], sentence_labels.shape[0], replace=False)

        return [x_images[idxs, ...], y_images[idxs, ...]], sentence_labels[idxs, ...]

def plot_sequences(x, y, labels=None, output_path=None):

    ''' Draws a plot where sequences of numbers can be studied conveniently '''

    images = np.concatenate([x, y], axis=1)
    n_batches = images.shape[0]
    n_terms = images.shape[1]
    counter = 1
    for n_b in range(n_batches):
        for n_t in range(n_terms):
            plt.subplot(n_batches, n_terms, counter)
            plt.imshow(images[n_b, n_t, :, :, :])
            plt.axis('off')
            counter += 1
        if labels is not None:
            plt.title(labels[n_b, 0])

    if output_path is not None:
        plt.savefig(output_path, dpi=600)
    else:
        plt.show()

if __name__ == "__main__":

    # Test SortedNumberGenerator
    ag = SortedNumberGenerator(batch_size=8, subset='train', terms=4, positive_samples=4, predict_terms=4, image_size=64, color=True, rescale=False)
    for (x, y), labels in ag:
        plot_sequences(x, y, labels, output_path=r'resources/batch_sample_sorted.png')
        break

reason: import numpy,import scipy


text:
import numpy as np
from collections import namedtuple
from ssm.utils.model_tools import integ_psf_on_pixel
import hickle as hkl
import pickle

class PixelResponseModel:
    def __init__(self, response, xi, yi, pixel_size):
        from scipy.interpolate import RectBivariateSpline

        if np.max(xi) != -np.min(xi) or np.max(yi) != -np.min(yi):
            raise ValueError("The response needs to be defined in a square")

        self.response = response
        self.xi = xi
        self.yi = yi
        self.pixel_size = pixel_size
        self.response_spl = RectBivariateSpline(yi[:, 0], xi[0, :], response)

        self.model_size = np.max(xi)
        self.psf = None

    @classmethod
    def from_file(cls, filename):
        loaded = pickle.load(open(filename, "rb"))
        pix_m = loaded["pix_m"]
        cons = cls(pix_m.response, pix_m.xi, pix_m.yi, pix_m.pixel_size)
        cons.psf = loaded["psf"]
        return cons

    @classmethod
    def from_psf(cls, psf, pixel_size, nbins=40):
        from tqdm.auto import tqdm

        model_size = pixel_size / 2 + psf.radialsize
        xi, yi = np.meshgrid(
            np.linspace(-model_size, model_size, nbins),
            np.linspace(-model_size, model_size, nbins),
        )
        d = np.diff(xi.flatten())[0]
        path = list(zip(xi.flatten(), yi.flatten()))
        response = np.empty(len(path))
        for i, p in tqdm(enumerate(path), total=len(path)):
            response[i] = integ_psf_on_pixel(psf, p, pixel_size, epsabs=0.1)[0]
        response /= np.max(response)
        cons = cls(response.reshape(xi.shape), xi, yi, pixel_size)
        cons.psf = psf
        return cons

    def save_model(self, filename):
        pickle.dump({"psf": self.psf, "pix_m": self}, open(filename, "wb"))

    def __call__(self, x, y, grid=False):
        return self.response_spl(x, y, grid=grid)

max_stars_repo_name: sflis/SSM-analysis

max_issues_repo_name: sflis/SSM-analysis

max_forks_repo_name: sflis/SSM-analysis

reason: import numpy,from scipy


text:
import numpy as np  # to deal with array and matrix calculation easily

from sklearn.datasets.samples_generator import make_blobs  # generate data

import matplotlib.pyplot as plt  # to plot
import matplotlib.animation as animation  # to animate the plot
import matplotlib.patches as mpatches  # to add legends easily

# region variables

# K-means variables
X = None  # features matrix
N_SAMPLES = 500  # number of sample to create artificially
K = 5  # value of k for the k-means algorithm
centroids = None  # will be the centroids array
last_centroids = None  # will keep in memory the previous array of centroids

# elbow method variables
ELBOW_METHOD = False  # to choose to launch and see the elbow method for the current dataset
K_MAX = 20  # max value of k to test with the elbow method ( will test k from 1 to K_MAX)

# plot variables
INTERVAL = 500  # speed of the animation frame (ms)
fig = plt.figure(1, figsize=(8, 4.5))  # to plot
ax = fig.add_subplot(1, 1, 1)  # the axes
ani = None  # to animate the plot
started = False  # indicates if the program is started
first_frame = -1  # indicates the first frame number to know the current iteration value
# choose the color respectively to the centroids (add color to have more than K = 10)
colors = {0: '#FF00FF', 1: '#999999', 2: '#2F9599', 3: 'red', 4: 'blue',
          5: 'yellow', 6: 'pink', 7: 'orange', 8: 'brown', 9: 'purple'}
# endregion

# region load dataset
def load_dataset():
    """
    Make an artificial dataset with random centered values
    :return: the dataset with point of coordinates (x, y) in [0, 30]x[0, 30]
    """
    X, _ = make_blobs(n_samples=N_SAMPLES, centers=25, n_features=2, cluster_std=2)
    return X + 15  # + 15 to have positive values (not needed of course)
# endregion

# region K-means
def euclidean_distance(X, centroid):
    """
    Calculate the euclidean distance between the samples and a given centroid
    :param X: matrix of sample coordinates
    :param centroid: given centroid from which the distance is calculated
    :return: array of each distance between the samples a the centroid
    """
    # Calculation of the distance of each coordinate's dimension to the centroid coordinate value associated
    coord_difference_to_square = (X - np.full(X.shape[1], centroid, dtype=np.ndarray))**2
    # Return the euclidean distance value
    return coord_difference_to_square.sum(axis=1, keepdims=True)**0.5

def nearest_centroid(X, centroids):
    """
    Indicates the associated centroid (index of the centroid in the array) for each sample of X
    :param X: matrix of sample coordinates
    :param centroids: array of centroids
    :return: index array of the associated centroid for each sample
    """
    # Create an array
    # (used to be stacked with arrays of euclidean distances from each sample to each centroid)
    dist_to_centroids = np.full((X.shape[0], 1), float('inf'))
    # Calculate for each centroid the euclidean distance from each sample
    for i, centroid in enumerate(centroids):
        # Stack the distance arrays to form a matrix
        dist_to_centroids = np.hstack((dist_to_centroids, euclidean_distance(X, centroid)))
    # Return the index of the nearest centroid for each sample
    return np.argmin(dist_to_centroids, axis=1) - 1 # - 1 as the first index if the unused fulled-'inf' array

def adjust_centroid(X, nearest_centroid_array, current_centroids):
    """
    Calculate the mean coordinates of the associated samples for each centroid
    :param X: matrix of sample coordinates
    :param nearest_centroid_array: array of index of the nearest centroid for each sample
    :param current_centroids: current centroid coordinates
    :return: the updated coordinates for each centroid as an array
    """
    # Calculation of the updated coordinate for each centroid
    for i, centroid in enumerate(np.unique(nearest_centroid_array)):  # for each centroid
        # pick up the indexes of nearest samples
        nearest_samples_index = np.where(nearest_centroid_array == centroid)
        # pick up their coordinates
        nearest_samples_coord = X[:][nearest_samples_index]
        # Calculate the mean coordinate for new coordinate of centroid
        new_coord = np.mean(nearest_samples_coord, axis=0)
        # change the centroid
        current_centroids[i] = new_coord
    return current_centroids

def calculate_cost(X, current_centroids):
    """
    Calculation the cost with the current centroids
    :param X: matrix of sample coordinates
    :param current_centroids: current centroids
    :return: cost of the k_means step
    """
    # Work as the the nearest_centroid() method
    # but return the sum of min values of euclidean distances for each samples to the centroids
    # (not the array of indexes of nearest centroid for each sample )

    # Create an array
    # (used to be stacked with arrays of euclidean distances from each sample to each centroid)
    dist_to_centroids = np.full((X.shape[0], 1), float('inf'))
    for i, centroid in enumerate(current_centroids):
        # Stack the distance arrays to form a matrix
        dist_to_centroids = np.hstack((dist_to_centroids, euclidean_distance(X, centroid)))
    # Return the sum of min values of euclidean distances for each samples to the centroids
    return sum(np.amin(dist_to_centroids, axis=1))

def k_means(frame_number):
    """
    K-means algorithm (main function)
    :param frame_number: frame number of the matplotlib.animation
    """
    global centroids, first_frame, last_centroids
    # Wait for the program to start
    if started:
        # To know if it is the first time the k_means() function is called
        if first_frame == -1:  # keep the first frame value
            first_frame = frame_number  # indicates the first frame value
            centroids = np.random.rand(K, 2) * np.array([0, 30])  # create random centroids on the x2-axis
            last_centroids = centroids[:] / 2  # to make sure to have different values

        # If an update is necessary (or while not converging)
        if not np.array_equal(centroids, last_centroids):
            # display the current centroid and associated sampels
            display(frame_number - first_frame)
            # keep the current centroid before the adjustments to know if there is a convergence
            last_centroids = np.array(centroids)
            # pick up the nearest centroid indexes of each samples
            nearest = nearest_centroid(X, centroids)
            # adjust the current centroids
            centroids = adjust_centroid(X, nearest, centroids)
# endregion

# region Elbow method
def elbow_method(k_max=K_MAX):
    """
    Elbow method : calculate the cost at the convergence of k_means algo for a range of k values
    :param k_max: indicate the range of k values to test (from 1 to k_max)
    """
    # Initialize the cost array
    cost_array = np.full(k_max, 0)

    # calculation the cost at the convergence of k_means algo for each k value
    for k, cost in enumerate(range(1, k_max + 1)):
        k_centroids =  np.random.rand(k + 1, 2) * np.array([0, 30])  # create random centroids on the x2-axis
        # Initialize an array to keep the last centroid values to know when there is a convergence
        last_centroids = np.array([])
        # Update the centroid while not convergence
        while not np.array_equal(k_centroids, last_centroids):
            # keep the current centroid before the adjustments to know if there is a convergence
            last_centroids = np.array(k_centroids)
            # pick up the nearest centroid indexes of each samples
            nearest = nearest_centroid(X, k_centroids)
            # adjust the current centroids
            k_centroids = adjust_centroid(X, nearest, k_centroids)
            # calculation of the current cost
            cost = calculate_cost(X, k_centroids)
        # keep the cost at the convergence
        cost_array[k] = cost
    # display the elbow method curve
    display_elbow(cost_array)
# endregion

# region display
def display(iteration=0):
    """
    Display the k_means algorithm
    :param iteration: number of current iteration/step
    """
    global centroids

    # clear the plot
    ax.clear()

    # set the labels to describe the plot
    ax.set_xlabel("x1")  # first feature
    ax.set_ylabel("x2")  # second feature
    ax.set_xlim(0, 30)  # fix the x1 axis
    ax.set_ylim(0, 30)  # fix the x2 axis

    # wait for the program to start
    if started:
        # add a big legend to describe the state of the k_means algorithm
        label = 'K-means :\n'
        label += 'K : {0}\n'.format(K)
        label += 'Iteration : {0}\n'.format(iteration)
        label += 'Cost : {0}\n'.format(round(calculate_cost(X, centroids), 3))

        # add the created legend
        plt.legend(handles=[mpatches.Patch(label=label)])

        # pick up the nearest centroid to have the respective color of each sample
        nearest = nearest_centroid(X, centroids)

        # plot each sample
        for i, sample in enumerate(X):
            ax.scatter(sample[0], sample[1], c=colors[nearest[i]])

        # plot each centroid
        for i, centroid in enumerate(centroids):
            ax.scatter(centroid[0], centroid[1], c=colors[i], edgecolors='black')

    else:  # if the program is not started
        # plot each sample as black points (no centroid yet)
        for i, sample in enumerate(X):
            ax.scatter(sample[0], sample[1], c='black')

def display_elbow(cost_array):
    """
    Display the elbow method algorithm
    :param cost_array: array of cost for each k
    """
    fig_elbow = plt.figure(2)  # to plot the elbow method
    ax_elbow = fig_elbow.add_subplot(1, 1, 1)  # the axes of the elbow method graph

    # set the labels to describe the plot
    ax_elbow.set_title('Elbow method for k to 1 to' + str(K_MAX))
    ax_elbow.set_xlabel('k')
    ax_elbow.set_ylabel('cost')
    ax_elbow.set_ylim(0, np.max(cost_array)*1.1)  # fix the cost axis

    # plot the elbow curve
    ax_elbow.plot(range(1, len(cost_array) + 1), cost_array)  # the axes

def key_pressed(event):
    """
    To start to run the programme by enter key
    :param event: key_press_event
    """
    if event.key == 'enter':
        global started
        started = not started
# endregion

if __name__ == '__main__':
    X = load_dataset()  # load artificial dataset

    if ELBOW_METHOD:  # launch the elbow method if it is chosen
        elbow_method()

    # connect to the key press event to start the program
    fig.canvas.mpl_connect('key_press_event', key_pressed)
    # to animate the plot and launch the gradient descent update
    ani = animation.FuncAnimation(fig, k_means, interval=INTERVAL)
    display()  # first plot
    plt.show()  # show the plot

max_stars_repo_name: AxelThevenot/K-Means

max_issues_repo_name: AxelThevenot/K-Means

max_forks_repo_name: AxelThevenot/K-Means


text:
#!/usr/bin/python3
"""Curses-based application to aid calibrating joint friction.

Rotates the joints of the OneJoint robot with constant velocity using a simple
PI velocity controller.
Prints averaged torque/current and velocity measurements from which the
friction can be estimated.

See `--help` for options.
"""

import argparse
import curses
import os
import time
import numpy as np
from ament_index_python.packages import get_package_share_directory

from robot_interfaces import one_joint
import robot_fingers

N_JOINTS = 1

CURRENT_TO_TORQUE_FACTOR = 0.02 * 9

class AverageBuffer:
    def __init__(self, size):
        self.size = size
        self.buffer = []

    def append(self, entry):
        self.buffer.append(entry)
        if len(self.buffer) > self.size:
            self.buffer = self.buffer[1:]

    def mean(self):
        return np.array(self.buffer).mean(axis=0)

def run_application(stdscr, robot, velocity_radps, buffer_size):

    stdscr.nodelay(True)

    enabled = False

    desired_velocity_radps = velocity_radps
    step_duration_s = 0.001

    velocity_buffer = AverageBuffer(buffer_size)
    measured_torque_buffer = AverageBuffer(buffer_size)
    applied_torque_buffer = AverageBuffer(buffer_size)

    # PI controller parameters
    kp = 0.1
    ki = 0.2
    max_error_integal = 10
    error_integral = np.zeros(N_JOINTS)
    desired_velocity_radps = np.ones(N_JOINTS) * desired_velocity_radps

    # start by sending a zero torque command
    action = one_joint.Action()
    t = robot.append_desired_action(action)
    obs = robot.get_observation(t)

    last_update = 0
    while True:
        t = robot.append_desired_action(action)
        obs = robot.get_observation(t)
        applied_action = robot.get_applied_action(t)

        if enabled:
            # run simple PI velocity controller
            velocity_error = desired_velocity_radps - obs.velocity
            error_integral += velocity_error * step_duration_s
            error_integral = np.clip(
                error_integral, -max_error_integal, +max_error_integal
            )
            torque = kp * velocity_error + ki * error_integral

            action = one_joint.Action(torque=torque)
        else:
            action = one_joint.Action()

        velocity_buffer.append(obs.velocity)
        measured_torque_buffer.append(obs.torque)
        applied_torque_buffer.append(applied_action.torque)

        now = time.time()
        if now - last_update > 0.1:
            last_update = now

            # generate the output
            stdscr.clear()
            line = 0
            # help line at the bottom
            stdscr.addstr(
                curses.LINES - 1,
                0,
                " s: start/stop motors | q: quit ",
                curses.A_STANDOUT,
            )

            stdscr.addstr(
                line,
                0,
                "Values are averaged over {:.3f} seconds".format(
                    buffer_size * step_duration_s
                ),
            )

            line += 2
            stdscr.addstr(line, 0, "Status:", curses.A_BOLD)
            stdscr.addstr(
                line, len("Status:") + 1, "RUNNING" if enabled else "STOPPED"
            )

            line += 2
            stdscr.addstr(line, 0, "Commanded Torque:", curses.A_BOLD)
            for i, p in enumerate(applied_torque_buffer.mean()):
                line += 1
                stdscr.addstr(
                    line,
                    4,
                    "Joint {}: {: .3f} Nm / {: .3f} A".format(
                        i, p, p / CURRENT_TO_TORQUE_FACTOR
                    ),
                )

            line += 2
            stdscr.addstr(line, 0, "Measured Torque:", curses.A_BOLD)
            for i, p in enumerate(measured_torque_buffer.mean()):
                line += 1
                stdscr.addstr(
                    line,
                    4,
                    "Joint {}: {: .3f} Nm / {: .3f} A".format(
                        i, p, p / CURRENT_TO_TORQUE_FACTOR
                    ),
                )

            line += 2
            stdscr.addstr(
                line, 0, "Velocity (actual / desired):", curses.A_BOLD
            )
            for i, (p, d) in enumerate(
                zip(velocity_buffer.mean(), desired_velocity_radps)
            ):
                line += 1
                stdscr.addstr(
                    line, 4, "Joint {}: {:.3f} / {:3f}".format(i, p, d)
                )

            stdscr.refresh()

        c = stdscr.getch()
        if c == ord("q"):
            return
        if c == ord("s"):
            enabled = not enabled
        else:
            pass

def main():
    argparser = argparse.ArgumentParser(description=__doc__)
    argparser.add_argument(
        "--canport",
        type=str,
        required=True,
        help="The CAN port to which the joint is connected (e.g. 'can0').",
    )
    argparser.add_argument(
        "--velocity",
        type=float,
        default=np.pi / 2,
        help="Velocity at which joints are moving [rad/s].",
    )
    argparser.add_argument(
        "--average-window",
        type=int,
        default=1000,
        help="Number of values used for averaging.",
    )
    args = argparser.parse_args()

    # load the default config file
    config_file_path = os.path.join(
        get_package_share_directory("robot_fingers"),
        "config",
        "onejoint_friction_calibration.yml",
    )

    config = robot_fingers.OneJointConfig.load_config(config_file_path)
    config.can_ports = [args.canport]
    robot_data = one_joint.SingleProcessData()
    robot_backend = robot_fingers.create_one_joint_backend(robot_data, config)
    robot_frontend = one_joint.Frontend(robot_data)

    robot_backend.initialize()

    curses.wrapper(
        lambda stdscr: run_application(
            stdscr, robot_frontend, args.velocity, args.average_window
        )
    )

if __name__ == "__main__":
    main()


text:
import numpy as np
from layers.activation_layer import *
from layers.dense_layer import *
from layers.gradient_check import *

def affine_relu_forward(X, W, b):
    y, affine_cache = affine_forward(X, W, b)
    outputs, relu_cache = relu_forward(y)
    cache = (affine_cache, relu_cache)
    return outputs, cache

def affine_relu_backward(dout, cache):
    affine_cache, relu_cache = cache
    dy = relu_backward(dout, relu_cache)
    dx, dw, db = affine_backward(dy, affine_cache)
    return dx, dw, db

if __name__ == '__main__':
    np.random.seed(231)
    x = np.random.randn(2, 3, 4)
    w = np.random.randn(12, 10)
    b = np.random.randn(10)
    dout = np.random.randn(2, 10)

    out, cache = affine_relu_forward(x, w, b)
    dx, dw, db = affine_relu_backward(dout, cache)

    dx_num = eval_numerical_gradient_array(lambda x: affine_relu_forward(x, w, b)[0], x, dout)
    dw_num = eval_numerical_gradient_array(lambda w: affine_relu_forward(x, w, b)[0], w, dout)
    db_num = eval_numerical_gradient_array(lambda b: affine_relu_forward(x, w, b)[0], b, dout)

    print('Testing affine_relu_forward:')
    print('dx error: ', rel_error(dx_num, dx))
    print('dw error: ', rel_error(dw_num, dw))
    print('db error: ', rel_error(db_num, db))


text:
import unittest
from unittest import TestCase

from e2cnn.nn import *
from e2cnn.gspaces import *

import torch

import numpy as np

import random

class TestNonLinearitiesRotations(TestCase):

    def test_cyclic_norm_relu(self):
        N = 8
        g = Rot2dOnR2(N)

        r = FieldType(g, list(g.representations.values()) * 4)

        nnl = NormNonLinearity(r, function='n_relu')

        nnl.check_equivariance()

    def test_cyclic_norm_sigmoid(self):
        N = 8
        g = Rot2dOnR2(N)

        r = FieldType(g, list(g.representations.values()) * 4)

        nnl = NormNonLinearity(r, function='n_sigmoid')

        nnl.check_equivariance()

    def test_cyclic_pointwise_relu(self):
        N = 8
        g = Rot2dOnR2(N)

        reprs = [r for r in g.representations.values() if 'pointwise' in r.supported_nonlinearities]

        r = FieldType(g, reprs)

        # nnl = PointwiseNonLinearity(r, function='p_relu')
        nnl = ReLU(r)

        nnl.check_equivariance()

    def test_cyclic_pointwise_sigmoid(self):
        N = 8
        g = Rot2dOnR2(N)

        reprs = [r for r in g.representations.values() if 'pointwise' in r.supported_nonlinearities]

        r = FieldType(g, reprs)

        nnl = PointwiseNonLinearity(r, function='p_sigmoid')

        nnl.check_equivariance()

    def test_cyclic_gated_one_input_shuffled_gated(self):
        N = 8
        g = Rot2dOnR2(N)

        reprs = [r for r in g.representations.values() if 'gated' in r.supported_nonlinearities]
        ngates = len(reprs)
        reprs += [g.trivial_repr] * ngates

        gates = ['gated'] * ngates + ['gate'] * ngates

        r = FieldType(g, reprs)

        nnl = GatedNonLinearity1(r, gates=gates)
        nnl.check_equivariance()

    def test_cyclic_gated_one_input_sorted_gated(self):
        N = 8
        g = Rot2dOnR2(N)

        reprs = [r for r in g.representations.values() if 'gated' in r.supported_nonlinearities]

        r = FieldType(g, reprs).sorted()

        ngates = len(r)

        reprs = [g.trivial_repr] * ngates

        gates = ['gated'] * ngates + ['gate'] * ngates

        r = r + FieldType(g, reprs)

        nnl = GatedNonLinearity1(r, gates=gates)
        nnl.check_equivariance()

    def test_cyclic_gated_one_input_all_shuffled(self):
        N = 8
        g = Rot2dOnR2(N)

        reprs = [r for r in g.representations.values() if 'gated' in r.supported_nonlinearities]

        ngates = len(reprs)

        reprs += [g.trivial_repr] * ngates

        gates = ['gated'] * ngates + ['gate'] * ngates

        t = list(zip(reprs, gates))

        random.shuffle(t)

        reprs, gates = zip(*t)

        r = FieldType(g, reprs)

        nnl = GatedNonLinearity1(r, gates=gates)
        nnl.check_equivariance()

    def test_cyclic_gated_two_inputs_shuffled_gated(self):
        N = 8
        g = Rot2dOnR2(N)

        gated = [r for r in g.representations.values() if 'gated' in r.supported_nonlinearities]
        ngates = len(gated)
        gates = [g.trivial_repr] * ngates

        gates = FieldType(g, gates)
        gated = FieldType(g, gated)

        nnl = GatedNonLinearity2((gates, gated))
        nnl.check_equivariance()

    def test_cyclic_gated_two_inputs_sorted_gated(self):
        N = 8
        g = Rot2dOnR2(N)

        gated = [r for r in g.representations.values() if 'gated' in r.supported_nonlinearities]

        gated = FieldType(g, gated).sorted()

        ngates = len(gated)

        gates = [g.trivial_repr] * ngates
        gates = FieldType(g, gates)

        nnl = GatedNonLinearity2((gates, gated))
        nnl.check_equivariance()

    def test_cyclic_concat_relu(self):
        N = 8
        g = Rot2dOnR2(N)

        reprs = [r for r in g.representations.values() if 'concatenated' in r.supported_nonlinearities]

        for rep in reprs:
            print(rep.name)
            r = FieldType(g, [rep])
            nnl = ConcatenatedNonLinearity(r, function='c_relu')
            nnl.check_equivariance()

    def test_cyclic_vectorfield(self):
        N = 8
        g = Rot2dOnR2(N)

        reprs = [r for r in g.representations.values() if 'vectorfield' in r.supported_nonlinearities] * 8

        r = FieldType(g, reprs)
        nnl = VectorFieldNonLinearity(r)
        nnl.check_equivariance(atol=2e-6)

    def test_cyclic_induced_norm_relu(self):

        N = 15
        g = Rot2dOnR2(N)

        sg_id = 5
        sg, _, _ = g.fibergroup.subgroup(sg_id)

        r = FieldType(g, [g.induced_repr(sg_id, sg.irrep(k)) for k in range(1, int(sg.order() // 2))] * 4).sorted()
        nnl = InducedNormNonLinearity(r, function='n_relu')
        nnl.check_equivariance()

    def test_so2_norm_relu(self):

        g = Rot2dOnR2(-1, 10)

        r = FieldType(g, list(g.representations.values()) * 4)

        nnl = NormNonLinearity(r, function='n_relu')

        nnl.check_equivariance()

    def test_so2_norm_sigmoid(self):
        g = Rot2dOnR2(-1, 10)

        r = FieldType(g, list(g.representations.values()) * 4)

        nnl = NormNonLinearity(r, function='n_sigmoid')

        nnl.check_equivariance()

    def test_so2_pointwise_relu(self):
        g = Rot2dOnR2(-1, 10)

        reprs = [r for r in g.representations.values() if 'pointwise' in r.supported_nonlinearities]

        r = FieldType(g, reprs)

        nnl = PointwiseNonLinearity(r, function='p_relu')

        nnl.check_equivariance()

    def test_so2_pointwise_sigmoid(self):
        g = Rot2dOnR2(-1, 10)

        reprs = [r for r in g.representations.values() if 'pointwise' in r.supported_nonlinearities]

        r = FieldType(g, reprs)

        nnl = PointwiseNonLinearity(r, function='p_sigmoid')

        nnl.check_equivariance()

    def test_so2_gated_one_input_shuffled_gated(self):
        g = Rot2dOnR2(-1, 10)

        reprs = [r for r in g.representations.values() if 'gated' in r.supported_nonlinearities]
        ngates = len(reprs)
        reprs += [g.trivial_repr] * ngates

        gates = ['gated'] * ngates + ['gate'] * ngates

        r = FieldType(g, reprs)

        nnl = GatedNonLinearity1(r, gates=gates)
        nnl.check_equivariance()

    def test_so2_gated_one_input_sorted_gated(self):
        g = Rot2dOnR2(-1, 10)

        reprs = [r for r in g.representations.values() if 'gated' in r.supported_nonlinearities]

        r = FieldType(g, reprs).sorted()

        ngates = len(r)

        reprs = [g.trivial_repr] * ngates

        gates = ['gated'] * ngates + ['gate'] * ngates

        r = r + FieldType(g, reprs)

        nnl = GatedNonLinearity1(r, gates=gates)
        nnl.check_equivariance()

    def test_so2_gated_one_input_all_shuffled(self):
        g = Rot2dOnR2(-1, 10)

        reprs = [r for r in g.representations.values() if 'gated' in r.supported_nonlinearities]

        ngates = len(reprs)

        reprs += [g.trivial_repr] * ngates

        gates = ['gated'] * ngates + ['gate'] * ngates

        t = list(zip(reprs, gates))

        random.shuffle(t)

        reprs, gates = zip(*t)

        r = FieldType(g, reprs)

        nnl = GatedNonLinearity1(r, gates=gates)
        nnl.check_equivariance()

    def test_so2_gated_two_inputs_shuffled_gated(self):
        g = Rot2dOnR2(-1, 10)

        gated = [r for r in g.representations.values() if 'gated' in r.supported_nonlinearities]
        ngates = len(gated)
        gates = [g.trivial_repr] * ngates

        gates = FieldType(g, gates)
        gated = FieldType(g, gated)

        nnl = GatedNonLinearity2((gates, gated))
        nnl.check_equivariance()

    def test_so2_gated_two_inputs_sorted_gated(self):
        g = Rot2dOnR2(-1, 10)

        gated = [r for r in g.representations.values() if 'gated' in r.supported_nonlinearities]

        gated = FieldType(g, gated).sorted()

        ngates = len(gated)

        gates = [g.trivial_repr] * ngates
        gates = FieldType(g, gates)

        nnl = GatedNonLinearity2((gates, gated))
        nnl.check_equivariance()

    def test_cyclic_gated1_error(self):
        N = 8
        g = Rot2dOnR2(N)

        for r in g.representations.values():
            if 'gated' not in r.supported_nonlinearities:
                r1 = FieldType(g, [r, g.trivial_repr])
                gates = ['gated', 'gate']
                self.assertRaises(AssertionError, GatedNonLinearity1, r1, gates=gates)

        for r in g.representations.values():
            if 'gate' not in r.supported_nonlinearities:
                r1 = FieldType(g, [g.trivial_repr, r])
                gates = ['gated', 'gate']
                self.assertRaises(AssertionError, GatedNonLinearity1, r1, gates=gates)

    def test_cyclic_gated2_error(self):
        N = 8
        g = Rot2dOnR2(N)

        for r in g.representations.values():
            if 'gated' not in r.supported_nonlinearities:
                gated = FieldType(g, [r])
                gates = FieldType(g, [g.trivial_repr])
                self.assertRaises(AssertionError, GatedNonLinearity2, (gates, gated))

        for r in g.representations.values():
            if 'gate' not in r.supported_nonlinearities:
                gated = FieldType(g, [g.trivial_repr])
                gates = FieldType(g, [r])
                self.assertRaises(AssertionError, GatedNonLinearity2, (gates, gated))

    def test_cyclic_norm_error(self):
        N = 8
        g = Rot2dOnR2(N)

        for r in g.representations.values():

            if 'norm' not in r.supported_nonlinearities:
                r1 = FieldType(g, [r])
                self.assertRaises(AssertionError, NormNonLinearity, r1)

    def test_cyclic_pointwise_error(self):
        N = 8
        g = Rot2dOnR2(N)

        for r in g.representations.values():

            if 'pointwise' not in r.supported_nonlinearities:
                r1 = FieldType(g, [r])
                self.assertRaises(AssertionError, PointwiseNonLinearity, r1)

    def test_cyclic_concat_error(self):
        N = 8
        g = Rot2dOnR2(N)

        for r in g.representations.values():

            if 'concatenated' not in r.supported_nonlinearities:
                r1 = FieldType(g, [r])
                self.assertRaises(AssertionError, ConcatenatedNonLinearity, r1)

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


text:
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import xarray as xr

import keras.backend as K

from sklearn.preprocessing import StandardScaler

from ninolearn.IO.read_processed import data_reader
from ninolearn.plot.evaluation  import plot_correlation
from ninolearn.plot.prediction import plot_prediction
from ninolearn.learn.evaluation.skillMeasures import rmse
from ninolearn.learn.models.dem import DEM
from ninolearn.utils import print_header, include_time_lag, basin_means
from ninolearn.pathes import modeldir

K.clear_session()

#%% =============================================================================
# read data
# =============================================================================
reader = data_reader(startdate='1960-01', enddate='2017-12')

# NINO3.4 Index
nino34 = reader.read_csv('oni')
#%%
# Other indeces
iod = reader.read_csv('iod')
wwv = reader.read_csv('wwv_proxy')

# seasonal cycle
sc = np.cos(np.arange(len(nino34))/12*2*np.pi)

# SSH network metrics
network_ssh = reader.read_statistic('network_metrics', variable='zos',
                           dataset='ORAS4', processed="anom")

c2_ssh = network_ssh['fraction_clusters_size_2']
H_ssh = network_ssh['corrected_hamming_distance']

#wind stress
taux = reader.read_netcdf('taux', dataset='NCEP', processed='anom')
taux_WP_mean, taux_CP_mean, taux_EP_mean = basin_means(taux)

ssh = reader.read_netcdf('zos', dataset='ORAS4', processed='anom')

# decadel variation of leading eof
pca_dec = reader.read_statistic('pca', variable='dec_sst', dataset='ERSSTv5', processed='anom')['pca1']

#%% =============================================================================
#  process data
# =============================================================================
time_lag = 12
lead_time = 6
shift = 3

# TODO warm pool edge, how to better include temperature
feature_unscaled = np.stack((nino34, iod, wwv,  sc,#pca_dec, sc,
                             taux_WP_mean, #taux_CP_mean, taux_EP_mean,
                             H_ssh, c2_ssh,
                             pca_dec
                             ), axis=1)

#feature_unscaled = np.concatenate((feature_unscaled, sst_equator),
#                                 axis=1)

scaler = StandardScaler()
Xorg = scaler.fit_transform(feature_unscaled)

Xorg = np.nan_to_num(Xorg)

X = Xorg[:-lead_time-shift,:]
futureX = Xorg[-lead_time-shift-time_lag:,:]

X = include_time_lag(X, max_lag=time_lag)
futureX =  include_time_lag(futureX, max_lag=time_lag)

yorg = nino34.values
y = yorg[lead_time + time_lag + shift:]

timey = nino34.index[lead_time + time_lag + shift:]
futuretime = pd.date_range(start='2019-01-01',
                                        end=pd.to_datetime('2019-01-01')+pd.tseries.offsets.MonthEnd(lead_time+shift),
                                        freq='MS')

test_indeces = (timey>='2012-01-01') & (timey<='2017-12-01')

train_indeces = np.invert(test_indeces)

trainX, trainy, traintimey = X[train_indeces,:], y[train_indeces], timey[train_indeces]
testX, testy, testtimey = X[test_indeces,:], y[test_indeces], timey[test_indeces]

#%% =============================================================================
# Deep ensemble
# =============================================================================
model = DEM()

model.set_parameters(layers=1, neurons=16, dropout=[0.0,0.5],
                     noise_in=[0.0,0.5], noise_mu=0.0, noise_sigma=0.0, noise_alpha=0.,
                     l1_hidden=[0, 0.2], l2_hidden=[0.,0.2],
                     l1_mu=0., l2_mu=0.,
                     l1_sigma=0, l2_sigma=0.,
                     l1_alpha=0., l2_alpha=0.,
                     lr=[0.0001, 0.01], batch_size=100, epochs=500, n_segments=5, n_members_segment=1,
                     patience=30, verbose=0, pdf='normal')

model.fit_RandomizedSearch(trainX, trainy, n_iter=100)

#%%
pred_mean, pred_std = model.predict(testX)

loss = model.evaluate(testy, pred_mean, pred_std)
print_header(f"Loss: {loss}")

if model.pdf == "normal":
    ens_dir=f'ensemble_lead{lead_time}'
elif model.pdf == "skewed":
    ens_dir=f'ensemble_skew_lead{lead_time}'
else:
    ens_dir=f'simple_ensemble_lead{lead_time}'

model.save(location=modeldir, dir_name=ens_dir)

#%% =============================================================================
# Plots
# =============================================================================
plt.close("all")

# Scores during trianing
plt.subplots()

for h in model.history:
    plt.plot(h.history[f'val_{model.loss_name}'], c='k')
    plt.plot(h.history[f'{model.loss_name}'], c='r')

plt.plot(np.nan, c='k', label='test')
plt.plot(np.nan, c='r', label='train')
plt.legend()

#%% just for testing the loading function delete and load the model
plt.close("all")
del model
model = DEM()
model.load(location=modeldir, dir_name=ens_dir)

pred_mean, pred_std = model.predict(testX)
predtrain_mean, predtrain_std = model.predict(trainX)
predfuture_mean, predfuture_std =  model.predict(futureX)

# Predictions
plt.subplots(figsize=(15,3.5))
plt.axhspan(-0.5,
            -6,
            facecolor='blue',
            alpha=0.1,zorder=0)

plt.axhspan(0.5,
            6,
            facecolor='red',
            alpha=0.1,zorder=0)

plt.xlim(timey[0],futuretime[-1])
plt.ylim(-3,3)

# test
plot_prediction(testtimey, pred_mean, std=pred_std, facecolor='royalblue', line_color='navy')

# train
#predtrain_mean[traintimey=='2002-02-01'] = np.nan
plot_prediction(traintimey, predtrain_mean, std=predtrain_std, facecolor='lime', line_color='g')

# future
plot_prediction(futuretime, predfuture_mean, std=predfuture_std, facecolor='orange', line_color='darkorange')

# observation
plt.plot(timey, y, "k")

pred_rmse = round(rmse(testy, pred_mean),2)
plt.title(f"Lead time: {lead_time} month, RMSE (of mean): {pred_rmse}")
plt.grid()
plt.xlabel('Year')
plt.ylabel('NINO3.4 [K]')

# Seaonality of Standard deviations
plt.subplots()

xr_nino34 = xr.DataArray(nino34)
std_data = xr_nino34.groupby('time.month').std(dim='time')

pd_pred_std = pd.Series(data=pred_std, index = testtimey)
xr_pred_std = xr.DataArray(pd_pred_std)
std_pred = xr_pred_std.groupby('time.month').mean(dim='time')

std_data.plot()
std_pred.plot()

# plot explained variance
plot_correlation(testy, pred_mean, testtimey - pd.tseries.offsets.MonthBegin(1), title="")

# Error distribution
plt.subplots()
plt.title("Error distribution")
error = pred_mean - testy

plt.hist(error, bins=16)

#%%
in_or_out = np.zeros((len(pred_mean)))
in_or_out[(testy>pred_mean-pred_std) & (testy<pred_mean+pred_std)] = 1
in_frac_1std = np.sum(in_or_out)/len(testy)

in_or_out = np.zeros((len(pred_mean)))
in_or_out[(testy>pred_mean-2*pred_std) & (testy<pred_mean+2*pred_std)] = 1
in_frac_2std = np.sum(in_or_out)/len(testy)

#%%
summed_weights = np.abs(model.ensemble[0].get_weights()[0]).sum(axis=1)
for j in range(1, len(model.ensemble)):
    summed_weights+= np.abs(model.ensemble[j].get_weights()[0]).sum(axis=1)

selected_var = np.argwhere(summed_weights>0.3*np.max(summed_weights))

n_var = feature_unscaled.shape[1]
selected_var_lead = selected_var//n_var

for i in range(feature_unscaled.shape[1]):
    index = (i==selected_var%n_var)
    lead = selected_var_lead[index]
    print(f'var nr{i} at:')
    print(lead)

max_stars_repo_name: pjpetersik/ninolearn

max_issues_repo_name: pjpetersik/ninolearn

max_forks_repo_name: pjpetersik/ninolearn


text:
import json
import numpy as np
import pickle
import random
import tensorflow as tf
import tensorflow_datasets as tfds

from tensorflow import keras
from tensorflow.keras.models import Sequential, Model
from tqdm import tqdm

from PIL import Image

# model_path = "/home/ec2-user/gan_submission_1/mnist/mnist_v2/model_keras.h5"
# file_path = '/home/ec2-user/gan_submission_1/mnist/mnist_v2/'
# tfds_name = 'mnist'
# im_size = 28
# n_channels = 1
# n_samples = 10000  # 10000
# n_classes = 10

# model_path = "/home/ec2-user/gan_submission_1/svhn/svhn_v2/model_keras.h5"
# file_path = '/home/ec2-user/gan_submission_1/svhn/svhn_v2/'
# tfds_name = 'svhn_cropped'
# im_size = 32
# n_channels = 3
# n_samples = 26032
# n_classes = 10

model_path = "/home/ec2-user/gan_submission_1/cifar10/cifar10_v2/model_keras.h5"
file_path = '/home/ec2-user/gan_submission_1/cifar10/cifar10_v2/'
tfds_name = 'cifar10'
im_size = 32
n_channels = 3
n_samples = 10000  # 10000
n_classes = 10

baseline = 'closest'  # 'gan', 'closest', 'ranch'

def get_evaluation_curve(model, xs, ys_label, ys_target, saliency):
    # data_struct = {key: np.copy(xs) for key in ['zero', 'max', 'mid', 'flip', 'rand', 'rand33', 'rand55']}
    data_struct = {key: np.copy(xs) for key in ['flip']}

    saliency = np.mean(saliency, axis=-1)
    instructions = [
        sorted([(saliency[i, x, y], x, y) for x in range(im_size) for y in range(im_size)], reverse=True)
        for i in range(n_samples)]

    curve_prob_map = {k: [] for k in data_struct.keys()}
    mask33 = [[[True] * im_size for _ in range(im_size)] for _ in range(n_samples)]
    mask55 = [[[True] * im_size for _ in range(im_size)] for _ in range(n_samples)]

    for idx in tqdm(range(im_size * im_size)):
        for i in range(n_samples):
            _, x, y = instructions[i][idx]
            # data_struct['zero'][i, x, y, :] = 0.0
            # data_struct['max'][i, x, y, :] = 1.0
            # data_struct['mid'][i, x, y, :] = 0.5
            data_struct['flip'][i, x, y, :] = 1.0 - data_struct['flip'][i, x, y, :]
            # data_struct['rand'][i, x, y, :] = np.random.rand(n_channels)
            #
            # for dx in range(-1, 2):
            #     for dy in range(-1, 2):
            #         if 0 <= x + dx < 28 and 0 <= y + dy < 28 and mask33[i][x + dx][y + dy]:
            #             mask33[i][x + dx][y + dy] = False
            #             data_struct['rand33'][i, x + dx, y + dy, :] = np.random.rand(n_channels)
            #
            # for dx in range(-2, 3):
            #     for dy in range(-2, 3):
            #
            #         if 0 <= x + dx < 28 and 0 <= y + dy < 28 and mask55[i][x + dx][y + dy]:
            #             mask55[i][x + dx][y + dy] = False
            #             data_struct['rand55'][i, x + dx, y + dy, :] = np.random.rand(n_channels)

        for eval_type, xs_this in data_struct.items():
            ys_pred = model.predict(xs_this)

            orig_score = np.mean([ys_pred[i, l] for i, l in enumerate(ys_label)])
            target_score = np.mean([ys_pred[i, l] for i, l in enumerate(ys_target)])

            curve_prob_map[eval_type].append(orig_score - target_score)

    output = {}
    for eval_type, curve in curve_prob_map.items():
        output[eval_type] = curve
    return output

if __name__ == '__main__':
    # Load model
    model = keras.models.load_model(model_path)

    # Load data
    ds = tfds.load(tfds_name, split='test')
    dslist = list(tfds.as_numpy(ds.take(n_samples)))

    ys_target = np.random.RandomState(seed=222).randint(n_classes - 1, size=n_samples)
    xs, ys_label = [], []
    for ind, sample in enumerate(dslist):
        xs.append(sample['image'])
        ys_label.append(sample['label'])
        if ys_target[ind] >= sample['label']:
            ys_target[ind] += 1

    xs = np.array(xs)
    xs = xs / 255.5
    ys_label = np.array(ys_label)

    assert xs.shape == (n_samples, im_size, im_size, n_channels)
    assert ys_label.shape == (n_samples,)
    assert ys_label.shape == (n_samples,)

    # Load explanations
    if baseline == 'gan':
        explain_types = ["expected_intgrad", "intgrad", "intgrad_base", "deeplift", "deeplift_base", "occlusion", "occlusion_base"]
    else:
        explain_types = ["intgrad_base", "deeplift_base", "occlusion_base"]

    explain_map = {et: np.load('%s%s_%s_%s.npy' % (file_path, et, baseline, n_samples)) for et in explain_types}

    # Evaluation
    all_curves = {}
    for et, saliency in explain_map.items():
        all_curves[et] = get_evaluation_curve(model, xs, ys_label, ys_target, saliency)

        with open('%sevaluation_%s_%d_flip' % (file_path, baseline, n_samples), "wb") as fp:
            pickle.dump(all_curves, fp)

# with open('%sevaluation_%d.json' % (file_path, n_samples), 'r') as fp:
#     data = json.load(fp)

max_stars_repo_name: pinjutien/DeepExplain

max_issues_repo_name: pinjutien/DeepExplain

max_forks_repo_name: pinjutien/DeepExplain


text:
from collections import defaultdict
import numpy as np
import tensorflow as tf
import tensorflow_probability as tfp

def kde_tf(query, dataset, kernel_tf, bandwidth=1):
    @tf.function(experimental_relax_shapes=True)
    def model_fn(arg):
        return kernel_tf(x=arg, y=query, bandwidth=bandwidth)
    return tf.math.reduce_mean(tf.vectorized_map(model_fn, dataset))

def laplacian_kernel(x, y, bandwidth=1):
    "Laplacian kernel"
    return np.exp(-1. * np.linalg.norm((x-y), ord=1) / bandwidth)

def laplacian_kernel_tf(x, y, bandwidth=1):
    return tf.math.exp(-1. * tf.norm((x - y), ord=1) / bandwidth)

class LaplaceLSH_tf:
    def __init__(self, dimension, bandwidth):
        poisson_param = dimension * 1. / bandwidth
        self.rep = np.random.poisson(poisson_param)
        self.axes = tf.random.uniform([self.rep], minval=0, maxval=dimension, dtype=tf.int32)
        self.thresholds = tf.random.uniform([self.rep], minval=0, maxval=1)
    def hash(self, point):
        return tf.gather(point, self.axes) < self.thresholds

class BinningLaplaceLSH_tf:
    def __init__(self, dimension, bandwidth):
        self.dimension = dimension
        self.delta = tf.random.gamma((1, dimension), 2, 1 / bandwidth)[0]
        weight = tf.random.uniform((1, dimension), 0, 1)[0]
        self.shift = tf.math.multiply(self.delta, weight)
        self.max_len = tf.math.floor(tf.nn.relu(tf.math.divide(1 - self.shift, self.delta)))

    def hash(self, point):
        pt_len = tf.nn.relu(tf.math.ceil((point - self.shift) / self.delta))
        return  tf.cast(tf.math.minimum(pt_len, self.max_len), tf.int32)

class FastLaplacianKDE_tf:
    def __init__(self, dataset, bandwidth, L):
        shape = tf.shape(dataset).numpy()
        n_points = shape[0]
        dimension = shape[1]
        sizes = np.random.binomial(n_points, L*1./n_points, L)
        self.hashed_points = []
        self.lshs = []
        self.L = L
        self.bandwidth = bandwidth
        self.dataset = dataset

        for size in sizes:
            hash_tf = BinningLaplaceLSH_tf(dimension, 2 * bandwidth)
            sample = np.random.choice(n_points, size, replace=False)
            result = {}
            for row_num in sample:
                row = dataset[row_num, :]
                row_hash = hash_tf.hash(row)
                hash_tuple_hash = hash(tuple(row_hash.numpy()))
                if hash_tuple_hash not in result:
                    result[hash_tuple_hash] = []
                result[hash_tuple_hash].append(row_num)
            self.lshs.append(hash_tf)
            self.hashed_points.append(result)

    def laplacian_kernel(self, x, y, bandwidth=1):
        return np.exp(-1. * np.linalg.norm((x-y), ord=1) / bandwidth)

    def kde(self, query):
        estimators = []
        for j in range(self.L):
            query_hash_tf = self.lshs[j].hash(query)
            query_hash = hash(tuple(query_hash_tf.numpy()))

            if query_hash not in self.hashed_points[j]:
                estimators.append(0)
            else:
                bin_size = len(self.hashed_points[j][query_hash])
                point_ind = self.hashed_points[j][query_hash][np.random.randint(bin_size)]
                point = self.dataset[point_ind, :]
                estimators.append(self.laplacian_kernel(query, point, 2*self.bandwidth) * bin_size * 1. / self.L)
        return np.mean(estimators)