text:
import numpy as np

from tbats.tbats import Components, Case, Context

class TestTBATSCase(object):

    def create_case(self, **components):
        return Case(Components(**components), Context())

    def test_fit_only_alpha(self):
        alpha = 0.7
        np.random.seed(333)
        T = 200

        l = 1
        y = [0] * T
        for t in range(0, T):
            d = np.random.normal()
            y[t] = l + d
            l = l + alpha * d

        bats_case = self.create_case(use_trend=False)
        model = bats_case.fit(y)
        # assert estimated values are close to actual values
        assert np.isclose(alpha, model.params.alpha, atol=0.1)

    def test_fit_with_trend(self):
        alpha = 0.8
        beta = 0.4
        np.random.seed(123)
        T = 400

        b = 0
        l = 1
        y = [0] * T
        for t in range(0, T):
            d = np.random.normal()
            y[t] = l + b + d
            l = l + b + alpha * d
            b = b + beta * d

        bats_case = self.create_case(use_trend=True)
        model = bats_case.fit(y)
        # assert estimated values are close to actual values
        assert np.isclose(alpha, model.params.alpha, atol=0.1)
        assert np.isclose(beta, model.params.beta, atol=0.2)

    def test_fit_trend_and_trigonometric_series(self):
        period_length = 6.5
        T = 300

        l = 0.2
        b = 1.
        s = 1.
        s_star = 0.
        y = [0] * T
        lam = 2 * np.pi / period_length
        for t in range(0, T):
            y[t] = l + s  # no noise
            l = l + b
            s_prev = s
            s = s_prev * np.cos(lam) + s_star * np.sin(lam)
            s_star = - s_prev * np.sin(lam) + s_star * np.cos(lam)

        case = self.create_case(use_trend=True, use_damped_trend=False,
                                seasonal_periods=[period_length], seasonal_harmonics=[1])
        fitted_model = case.fit(y)

        assert np.allclose(y, fitted_model.y_hat, atol=0.01)

max_stars_repo_name: series-temporais/tbats


text:
import numpy as np
import matplotlib.pyplot as plt
import sys

## Plot results
#  It plots the results in a unique figure
def plot_res(t,x,y,u,sp):
    plt.subplot(3,1,1)
    plt.plot(t.T,x.T)
    plt.xlabel('Time')
    plt.ylabel('States x')

    plt.subplot(3,1,2)
    plt.plot(t.T,y.T)
    plt.plot(t.T,sp.T,'k:')
    plt.xlabel('Time')
    plt.ylabel('Outputs y')

    plt.subplot(3,1,3)
    plt.plot(t.T,u.T)
    plt.xlabel('Time')
    plt.ylabel('Control u')
    umin = np.min(u)
    umax = np.max(u)
    axes = plt.gca()
    axes.set_ylim([umin*1.2,umax*1.2])

    plt.show()

def progress(val,end_val, bar_length=50):
    percent = float(val) / end_val
    hashes = '#' * int(round(percent * bar_length))
    spaces = ' ' * (bar_length - len(hashes))
    sys.stdout.write("\r    [{0}] {1}%".format(hashes + spaces, int(round(percent * 100))))
    sys.stdout.flush()

ext: py

reason: import numpy


text:
from keras.layers import Dense, LSTM, Embedding, Dropout, Activation
from keras import backend
from keras.models import *
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from utils import getActivationValue,layerName, hard_sigmoid
from keract import get_activations_single_layer
import os.path

def rmse(y_true, y_pred):
	return backend.sqrt(backend.mean(backend.square(y_pred - y_true), axis=-1))

class tsClass:
    def __init__(self):
        self.data = None
        self.X_orig = None
        self.X = None
        self.Y = None
        self.X_train = None
        self.y_train = None
        self.X_test = None
        self.y_test = None
        self.model = None
        self.unique_chars = None
        self.char_to_int = None
        self.int_to_char = None
        self.pad = 0
        self.numAdv = 0
        self.numSamples = 0
        self.perturbations = []

    def load_data(self, filename, seq_len, normalise_window):

        f = open(filename, 'r').read()
        data = f.split('\n')

        sequence_length = seq_len + 1
        result = []
        for index in range(len(data) - sequence_length):
            result.append(data[index: index + sequence_length])

        if normalise_window:
            result = self.normalise_windows(result)

        result = np.array(result)

        row = round(0.9 * result.shape[0])
        train = result[:int(row), :]
        np.random.shuffle(train)
        x_train = train[:, :-1]
        y_train = train[:, -1]
        x_test = result[int(row):, :-1]
        y_test = result[int(row):, -1]

        x_train = np.reshape(x_train, (x_train.shape[0], x_train.shape[1], 1))
        x_test = np.reshape(x_test, (x_test.shape[0], x_test.shape[1], 1))

        self.X_train = x_train
        self.y_train = y_train
        self.X_test = x_test
        self.y_test = y_test

        print(self.X_train.shape)
        print(self.y_train.shape)

        return [x_train, y_train, x_test, y_test]

    def normalise_windows(self, window_data):
        normalised_data = []
        for window in window_data:
            normalised_window = [((float(p) / float(window[0])) - 1) for p in window]
            normalised_data.append(normalised_window)
        return normalised_data

    def create_sequence(self, data):
        return self.normalise_windows(data)

    def load_model(self):
        filename = 'models/Ts.h5'

        if os.path.isfile(filename):
            self.model = load_model(filename,custom_objects={'rmse': rmse})
            self.model.compile(loss='mean_squared_error', optimizer='rmsprop',metrics=[rmse])
            self.model.summary()
        return self.model

    def layerName(self, layer):
        layerNames = [layer.name for layer in self.model.layers]
        return layerNames[layer]

    def train_model(self):
        self.model = Sequential()
        self.model.add(LSTM(input_shape=(None,1), units=50, return_sequences=True))
        self.model.add(Dropout(0.2))
        self.model.add(LSTM( 100, return_sequences=False))
        self.model.add(Dropout(0.2))
        self.model.add(Dense( units=1, activation="linear"))
        self.model.compile(loss="mse", optimizer="rmsprop", metrics=[rmse])

        print(self.model.summary())
        self.model.fit(self.X_train,self.y_train, validation_data=(self.X_test,self.y_test), epochs=300, batch_size=32)
        self.model.save('Ts.h5')

    def displayInfo(self,test):
        test = test[np.newaxis, :]
        smiles = self.vect_smile(test)
        output_value = np.squeeze(self.model.predict(test))
        print("current SMILES: ",smiles[0])
        print("current prediction: ", output_value)
        return output_value

    def updateSample(self,pred1,pred2,m,o):
        if abs(pred1-pred2) >= 1 and o == True:
            self.numAdv += 1
            self.perturbations.append(m)
        self.numSamples += 1
        self.displaySuccessRate()

    def displaySamples(self):
        print("%s samples are considered" % (self.numSamples))

    def displaySuccessRate(self):
        print("%s samples, within which there are %s adversarial examples" % (self.numSamples, self.numAdv))
        print("the rate of adversarial examples is %.2f\n" % (self.numAdv / self.numSamples))

    def displayPerturbations(self):
        if self.numAdv > 0:
            print("the average perturbation of the adversarial examples is %s" % (sum(self.perturbations) / self.numAdv))
            print("the smallest perturbation of the adversarial examples is %s" % (min(self.perturbations)))

    # calculate the lstm hidden state and cell state manually
    def cal_hidden_state(self, test):
        print('cal_hidden_state')
        print(test.shape)
        acx = get_activations_single_layer(self.model, np.array([test]), self.layerName(0))
        units = int(int(self.model.layers[1].trainable_weights[0].shape[1]) / 4)
        # print("No units: ", units)
        # lstm_layer = model.layers[1]
        W = self.model.layers[1].get_weights()[0]
        U = self.model.layers[1].get_weights()[1]
        b = self.model.layers[1].get_weights()[2]

        W_i = W[:, :units]
        W_f = W[:, units: units * 2]
        W_c = W[:, units * 2: units * 3]
        W_o = W[:, units * 3:]

        U_i = U[:, :units]
        U_f = U[:, units: units * 2]
        U_c = U[:, units * 2: units * 3]
        U_o = U[:, units * 3:]

        b_i = b[:units]
        b_f = b[units: units * 2]
        b_c = b[units * 2: units * 3]
        b_o = b[units * 3:]

        # calculate the hidden state value
        h_t = np.zeros((self.pad, units))
        c_t = np.zeros((self.pad, units))
        f_t = np.zeros((self.pad, units))
        h_t0 = np.zeros((1, units))
        c_t0 = np.zeros((1, units))

        for i in range(0, self.pad):
            f_gate = hard_sigmoid(np.dot(acx[i, :], W_f) + np.dot(h_t0, U_f) + b_f)
            i_gate = hard_sigmoid(np.dot(acx[i, :], W_i) + np.dot(h_t0, U_i) + b_i)
            o_gate = hard_sigmoid(np.dot(acx[i, :], W_o) + np.dot(h_t0, U_o) + b_o)
            new_C = np.tanh(np.dot(acx[i, :], W_c) + np.dot(h_t0, U_c) + b_c)
            c_t0 = f_gate * c_t0 + i_gate * new_C
            h_t0 = o_gate * np.tanh(c_t0)
            c_t[i, :] = c_t0
            h_t[i, :] = h_t0
            f_t[i, :] = f_gate

        return h_t, c_t, f_t


text:
import cv2
import numpy as np
import scipy.signal as sig

def blockshaped(arr, nrows, ncols):
    """
    Return an array of shape (n, nrows, ncols) where
    n * nrows * ncols = arr.size

    If arr is a 2D array, the returned array should look like n subblocks with
    each subblock preserving the "physical" layout of arr.
    """
    h, w = arr.shape
    assert h % nrows == 0, "{} rows is not evenly divisble by {}".format(h, nrows)
    assert w % ncols == 0, "{} cols is not evenly divisble by {}".format(w, ncols)
    return (arr.reshape(h//nrows, nrows, -1, ncols)
               .swapaxes(1,2)
               .reshape(-1, nrows, ncols))

def divide_in_tiles(frame, tile_size=10):
    return blockshaped(frame, tile_size, tile_size)

def resize(frame, dim):
    if len(frame.shape) >= 3:
        raise ValueError('frame must be 2 dimensional (grayscale), but current shape is: ', frame.shape)
    resized = cv2.resize(frame, dim)

    return resized

def gradients(frames):
    total_frames = frames.shape[2]

    merged_gradients = np.array([])

    framesDx = []
    framesDy = []
    framesDt = []

    convDt = sig.convolve(frames, np.array([[[1,0,-1]]]), mode='same')

    for frame_num in range(0, total_frames):
        frame_convDx = sig.convolve(frames[:,:,frame_num], np.array([[1, 0, -1]]), mode='same')
        frame_convDy = sig.convolve(frames[:,:,frame_num], np.array([[1], [0], [-1]]), mode='same')
        frame_convDt = convDt[:,:,frame_num]
        framesDx.append(frame_convDx)
        framesDy.append(frame_convDy)
        framesDt.append(frame_convDt)

        frame_gradient = np.stack((frame_convDx,frame_convDy,frame_convDt), axis=2).flatten()
        merged_gradients = np.concatenate((merged_gradients, frame_gradient), axis=None)

    return framesDx,framesDy,framesDt,merged_gradients

def spatiotemporal_cubes(frames, dim):
    tiled_frames = []
    for frame in frames:
        converted_frame = resize(frame, dim)
        tiled_frames.append(divide_in_tiles(converted_frame))

    return np.stack(tuple(tiled_frames), axis=3)

def video_sequence_cubes_and_gradients(frames,dim,cube_depth,tile_size):
    num_frames = len(frames)
    total_sequences = int(num_frames/cube_depth)

    total_tiles = int((dim[0]/tile_size) * (dim[1]/tile_size))
    size_gradients_flattened = tile_size*tile_size*cube_depth*3

    all_cubes = np.zeros((total_tiles,tile_size,tile_size,cube_depth))
    all_gradients = np.zeros((total_tiles, size_gradients_flattened))

    cubes = spatiotemporal_cubes(frames, dim)
    all_cubes = cubes

    cube_gradient = np.zeros((cubes.shape[0], size_gradients_flattened))
    for cube_number in range(0,cubes.shape[0]):
        cube_gradient[cube_number] = gradients(cubes[cube_number])[3]
    all_gradients = cube_gradient

    return all_cubes, all_gradients

def generate_features(frame_sequence,cube_depth,tile_size):
    #opencv is y,x
    #https://stackoverflow.com/questions/21248245/opencv-image-resize-flips-dimensions/22094421
    dimensions = [(20,20), (160,120), (40,30)]

    features = []

    for dim in dimensions:
        features.append(video_sequence_cubes_and_gradients(frame_sequence,dim,cube_depth,tile_size)[1])

    return np.concatenate(features, axis=0)

def process_frames(frames,cube_depth=5,tile_size=10):
    num_frames = len(frames)
    total_sequences = int(num_frames/cube_depth)

    sequence_features = []

    for i in range(0,total_sequences):
        frames_slice = frames[i*cube_depth:i*cube_depth + cube_depth]

        sequence_features.append(generate_features(frames_slice,cube_depth,tile_size))

    result = np.array(sequence_features) #shape (total_sequences, cubes, gradients) Lu et al applies PCA to reduce the gradients dimension
    return result.reshape(result.shape[0], result.shape[1]*result.shape[2]) #shape (total_sequences, features)

max_stars_repo_name: vinicius-pirees/phaino

max_issues_repo_name: vinicius-pirees/phaino

max_forks_repo_name: vinicius-pirees/phaino

reason: import numpy,import scipy


text:
from check_defn import *
import numpy as np
import random

def train(word_list, query_list, ensemble, iterations,l_rate,matrix,lamb):
    #counter is used to go to next ensemble
    counter = 0
    total_data = len(word_list)
    for iteration in range(iterations):
        #To keep the cost
        cost = 0
        #To store the entries for back propagation
        back = []
        #For keeping the accuracy
        accuracy = 0
        for q in range(ensemble):
            query_no = (q+counter)%total_data
            select_len = len(query_list[query_no])
            select = random.randint(1,100)%select_len
            query = query_list[query_no][select]
            #Gives a query for forward propagation
            (cost_minimum, found_word) = forward(query, word_list,matrix,lamb)
            cost+=cost_minimum
            #Backpropagation
            back.append((query, found_word, word_list[query_no]))
            if found_word == word_list[query_no]:
                accuracy+=1
        cost = cost/ensemble
        accuracy = accuracy*100/ensemble
        print('Iteration:',iteration,'\tAccuracy:',accuracy,'%', '\tCost:',cost)
        (correction,cost) = backward(back, matrix,lamb)
        print('\tSecondary Cost:',cost)
        matrix = matrix - l_rate*correction/(((accuracy+1)**2)*ensemble)
        counter+=ensemble
    return matrix

#Training dataset is wikipedia.txt
fp = open('wikipedia.txt','r')
#To store the words
word_list = []
#To store the queries
query_list = []
#Reads the data in the file
data = fp.readline()
while(data!=''):
    data = data.split(':')
    word_list.append(data[0])
    query_list.append(data[1].split())
    data = fp.readline()

#Initializes a matrix of all entries one of proper length
matrix = np.ones(len(create_vector(cost('a','a'))))
#Learning Rate : learning_rate
learning_rate = 1
#Ensemble or Batch size : ensemble_size
ensemble_size = len(word_list)//16
#Iterations : iterations
iterations = 32
#Lambda : lamb_da
lamb_da = 0

print('Learning rate:',learning_rate, 'Ensemble Size:', ensemble_size, 'Iterations:',iterations,'Lambda:',lamb_da)
matrix = train(word_list,query_list, ensemble_size, iterations, learning_rate, matrix, lamb_da)
np.save('matrix.npy', matrix)
np.savetxt('matrix.txt', matrix)

max_stars_repo_name: bridgesign/SpellCheck

max_issues_repo_name: bridgesign/SpellCheck

max_forks_repo_name: bridgesign/SpellCheck


text:
#!/usr/bin/python
# Dual polarisation beamforming: Track target and possibly calibrator for beamforming.

# The *with* keyword is standard in Python 2.6, but has to be explicitly
# imported in Python 2.5
from __future__ import with_statement

import time
import StringIO
import logging

import numpy as np

from katcorelib import (standard_script_options, verify_and_connect,
                        collect_targets, start_session, user_logger, ant_array)
import fbf_katcp_wrapper as fbf

class BeamformerReceiver(fbf.FBFClient):
    """KATCP client to beamformer receiver, with added metadata."""
    def __init__(self, name, server, rx_port, pol, meta_port, data_port, data_drive):
        user_logger.info('Connecting to server %r for beam %r...' % (server, name))
        logger = logging.getLogger('katcp')
        super(BeamformerReceiver, self).__init__(host=server, port=rx_port,
                                                 timeout=60, logger=logger)
        while not self.is_connected():
            user_logger.info('Waiting for TCP link to receiver server...')
            time.sleep(1)
        user_logger.info('Connected to server %r for beam %r' % (server, name))
        self.name = name
        self.pol = pol
        self.meta_port = meta_port
        self.data_port = data_port
        self.data_drive = data_drive
        self.obs_meta = {}

    def __repr__(self):
        return "<BeamformerReceiver %r -> %r at 0x%x>" % (self.name, self.pol, id(self))

    @property
    def inputs(self):
        return self.obs_meta.get('ants', [])

# Server where beamformer receivers are run
server = 'kat-dc2.karoo'
# beams = {'bf0': {'pol':'h', 'meta_port':'7152', 'data_port':'7150', 'rx_port':1235, 'data_drive':'/data1'},
#          'bf1': {'pol':'v', 'meta_port':'7153', 'data_port':'7151', 'rx_port':1236, 'data_drive':'/data2'}}
beams = [BeamformerReceiver('bf0', server, rx_port=1235, pol='h', meta_port=7152,
                            data_port=7150, data_drive='/data1'),
         BeamformerReceiver('bf1', server, rx_port=1236, pol='v', meta_port=7153,
                            data_port=7151, data_drive='/data2')]

def bf_inputs(cbf, bf):
    """Input labels associated with specified beamformer (*all* inputs)."""
    reply = cbf.req.dbe_label_input()
    return [] if not reply.succeeded else \
           [m.arguments[0] for m in reply.messages[1:] if m.arguments[3] == bf]

def select_ant(cbf, input, bf='bf0'):
    """Only use one antenna in specified beamformer."""
    # Iterate over *all* inputs going into the given beam
    for inp in bf_inputs(cbf, bf):
        status = 'beamformer input ' + inp + ':'
        weight = '1' if inp == input else '0'
        status += ' kept' if inp == input else ' zeroed'
        cbf.req.dbe_k7_beam_weights(bf, inp, *(1024 * [weight]))
        user_logger.info(status)

def get_weights(cbf):
    """Retrieve the latest gain corrections and their corresponding update times."""
    weights, times = {}, {}
    for sensor_name in cbf.sensor:
        if sensor_name.endswith('_gain_correction_per_channel'):
            sensor = cbf.sensor[sensor_name]
            input_name = sensor_name.split('_')[1]
            reading = sensor.get_reading()
            weights[input_name] = reading.value
            times[input_name] = reading.timestamp
    return weights, times

def phase_up(cbf, weights, inputs=None, bf='bf0', style='flatten'):
    """Phase up a group of antennas using latest gain corrections.

    The *style* parameter determines how the complex gain corrections obtained
    on the latest calibrator source will be turned into beamformer weights:

      - 'norm': Apply the complex gain corrections unchanged as weights,
        thereby normalising both amplitude and phase per channel.
      - 'phase': Only apply the phase correction, leaving the weight amplitudes
        equal to 1. This has the advantage of not boosting weaker inputs that
        will increase the noise level, but it also does not flatten the band.
      - 'flatten': Apply both amplitude and phase corrections, but preserve
        mean gain of each input. This flattens the band while also not boosting
        noise levels on weaker inputs.
      - 'scramble': Apply random phase corrections, just for the heck of it.

    Parameters
    ----------
    cbf : client object
        Object providing access to CBF (typically via proxy)
    weights : string->string mapping
        Gain corrections per input as returned by appropriate sensor
    inputs : None or sequence of strings, optional
        Names of inputs in use in given beamformer (default=all)
    bf : string, optional
        Name of beamformer instrument (one per polarisation)
    style : {'flatten', 'norm', 'phase', 'scramble'}, optional
        Processing done to gain corrections to turn them into weights

    """
    # Iterate over *all* inputs going into the given beam
    all_inputs = bf_inputs(cbf, bf)
    num_inputs = len(all_inputs)
    for inp in all_inputs:
        status = 'beamformer input ' + inp + ':'
        if (inputs is None or inp in inputs) and inp in weights and weights[inp]:
            # Extract array of complex weights from string representation
            weights_str = weights[inp]
            f = StringIO.StringIO(weights_str)
            orig_weights = np.loadtxt(f, dtype=np.complex, delimiter=' ')
            amp_weights = np.abs(orig_weights)
            phase_weights = orig_weights / amp_weights
            if style == 'norm':
                new_weights = orig_weights
                status += ' normed'
            elif style == 'phase':
                new_weights = phase_weights
                status += ' phased'
            elif style == 'flatten':
                # Get the average gain in the KAT-7 passband
                avg_amp = np.median(amp_weights[256:768])
                new_weights = orig_weights / avg_amp
                status += ' flattened'
            elif style == 'scramble':
                new_weights = np.exp(2j * np.pi * np.random.rand(1024))
                status += ' scrambled'
            else:
                raise ValueError('Unknown phasing-up style %r' % (style,))
            # Normalise weights by number of inputs to avoid overflow
            new_weights /= num_inputs
            # Reconstruct string representation of weights from array
            weights_str = ' '.join([('%+5.3f%+5.3fj' % (w.real, w.imag))
                                    for w in new_weights])
        else:
            # Zero the inputs that are not in use in the beamformer
            weights_str = ' '.join(1024 * ['0'])
            status += ' zeroed'
        cbf.req.dbe_k7_beam_weights(bf, inp, weights_str)
        user_logger.info(status)

def report_compact_traceback(tb):
    """Produce a compact traceback report."""
    print '--------------------------------------------------------'
    print 'Session interrupted while doing (most recent call last):'
    print '--------------------------------------------------------'
    while tb:
        f = tb.tb_frame
        print '%s %s(), line %d' % (f.f_code.co_filename, f.f_code.co_name, f.f_lineno)
        tb = tb.tb_next
    print '--------------------------------------------------------'

class BeamformerSession(object):
    """Context manager that ensures that beamformer is switched off."""
    def __init__(self, cbf, beams):
        self.cbf = cbf
        self.beams = beams

    def __enter__(self):
        """Enter the data capturing session."""
        return self

    def __exit__(self, exc_type, exc_value, traceback):
        """Exit the data capturing session, closing all streams."""
        if exc_value is not None:
            exc_msg = str(exc_value)
            msg = "Session interrupted by exception (%s%s)" % \
                  (exc_value.__class__.__name__,
                   (": '%s'" % (exc_msg,)) if exc_msg else '')
            if exc_type is KeyboardInterrupt:
                user_logger.warning(msg)
            else:
                user_logger.error(msg, exc_info=True)
        self.capture_stop()
        # Suppress KeyboardInterrupt so as not to scare the lay user,
        # but allow other exceptions that occurred in the body of with-statement
        if exc_type is KeyboardInterrupt:
            report_compact_traceback(traceback)
            return True
        else:
            return False

    def instrument_start(self, instrument):
        """Start given CBF instrument."""
        self.cbf.req.dbe_capture_start(instrument)
        user_logger.info('waiting 10s for stream %r to start' % (instrument,))
        time.sleep(10)

    def instrument_stop(self, instrument):
        """Stop given CBF instrument."""
        self.cbf.req.dbe_capture_stop(instrument)
        user_logger.info('waiting 10s for stream %r to stop' % (instrument,))
        time.sleep(10)

    def capture_start(self):
        """Enter the data capturing session, starting capture."""
        user_logger.info('Starting correlator (used for signal displays)')
        # Starting streams will issue metadata for capture
        # Allow long 10sec intervals to allow enough time to initiate data capture and to capture metadata
        # Else there will be collisions between the 2 beams
        for beam in self.beams:
            # Initialise receiver and setup server for data capture
            user_logger.info('Initialising receiver and stream for beam %r' %
                             (beam.name,))
            if not beam.rx_init(beam.data_drive, beam.obs_meta['half_band'],
                                beam.obs_meta['transpose']):
                raise RuntimeError('Could not initialise %r receiver' %
                                   (beam.name,))
            # Start metadata receiver before starting data transmit
            beam.rx_meta_init(beam.meta_port) # port
            self.instrument_start(beam.name)
            user_logger.info('beamformer metadata')
            beam.rx_meta(beam.obs_meta) # additional obs related info
            user_logger.info('waiting 10s to write metadata for beam %r' %
                             (beam.name,))
            time.sleep(10)
            # Start transmitting data
            user_logger.info('beamformer data for beam %r' % (beam.name,))
            beam.rx_beam(pol=beam.pol, port=beam.data_port)
            time.sleep(1)

    def capture_stop(self):
        """Exit the data capturing session, stopping the capture."""
        # End all receivers
        for beam in self.beams:
            user_logger.info('Stopping receiver and stream for beam %r' %
                             (beam.name,))
            beam.rx_stop()
            time.sleep(5)
            self.instrument_stop(beam.name)
            user_logger.info(beam.rx_close())
        user_logger.info('Stopping correlator (used for signal displays)')

# Set up standard script options
usage = "%prog [options] <'target'>"
description = "Perform a beamforming run on a specified target, optionally " \
              "visiting a gain calibrator beforehand to set beamformer weights."
parser = standard_script_options(usage, description)
# Add experiment-specific options
parser.add_option('-a', '--ants', default='all',
                  help="Antennas to include in beamformer (default='%default')")
parser.add_option('-b', '--buffercap', action='store_true',default=False,
              help="Use real-time dspsr pipeline (default='%default')")
parser.add_option('-t', '--target-duration', type='float', default=20,
                  help='Minimum duration to track the beamforming target, '
                       'in seconds (default=%default)')
parser.add_option('--half-band', action='store_true', default=False,
                  help='Use only inner 50% of output band')
parser.add_option('--reset', action="store_true", default=False,
                  help='Reset the gains to 160.')
# Set default value for any option (both standard and experiment-specific options)
parser.set_defaults(description='Beamformer observation', nd_params='off',
                    dump_rate=1.0, mode='bc16n400M1k')
# Parse the command line
opts, args = parser.parse_args()

# Check options and arguments and connect to KAT proxies and devices
if len(args) == 0:
    raise ValueError("Please specify the target")
with verify_and_connect(opts) as kat:
    cbf = kat.dbe7
    ants = kat.ants
    if opts.buffercap:  # set passband w.r.t. SPEAD rx
        bw, cfreq = [200000000, 100000000]
    else:
        bw, cfreq = [400000000, 200000000]

    for beam in ['bf0','bf1']:
        cbf.req.dbe_k7_beam_passband(beam, bw, cfreq)

    # We are only interested in the first target
    user_logger.info('Looking up main beamformer target...')
    target = collect_targets(kat, args[:1]).targets[0]
    # Ensure that the target is up
    target_elevation = np.degrees(target.azel()[1])
    if target_elevation < opts.horizon:
        raise ValueError("The desired target to be observed is below the horizon")

    # Start correlator capture session
    with start_session(kat, **vars(opts)) as corr_session:
        # Force delay tracking to be on
        opts.no_delays = False
        corr_session.standard_setup(**vars(opts))
        corr_session.capture_start()

        # Dictionary to hold observation metadata to send over to beamformer receiver
        for beam in beams:
            beam.obs_meta.update(vars(opts))
            beam.obs_meta['ants'] = [(ant.name + beam.pol) for ant in ants]
            beam.obs_meta['target'] = target.description
            if cal_target and len(ants) >= 4:
                beam.obs_meta['cal_target'] = cal_target.description

        if len(ants) > 1:
            user_logger.info('Setting beamformer weight to 1 for %d antennas' % (len(ants),))
            inputs = reduce(lambda inp, beam: inp + beam.inputs, beams, [])
            # set the beamformer weights to 1 as the phaseing is done in the f-engine
            weights = {}
            bf_weights_str = ' '.join(1024 * ['1'])
            for inp in inputs:
                weights[inp] = bf_weights_str
            for beam in beams:
                phase_up(cbf, weights, inputs=beam.inputs, bf=beam.name, style='norm')
                time.sleep(1)
        else:
            # The single-dish case does not need beamforming
            user_logger.info('Set beamformer weights to select single dish')
            for beam in beams:
                select_ant(cbf, input=beam.inputs[0], bf=beam.name)
                time.sleep(1)

        # Start beamformer session
        with BeamformerSession(cbf, beams) as bf_session:
            # Get onto beamformer target
            corr_session.label('track')
            corr_session.track(target, duration=0)
            # Only start capturing with beamformer once we are on target
            bf_session.capture_start()
            corr_session.track(target, duration=opts.target_duration)

max_stars_repo_name: ska-sa/katsdpscripts

max_issues_repo_name: ska-sa/katsdpscripts

max_forks_repo_name: ska-sa/katsdpscripts


text:
"""
Created on November 24, 2020

@author: Robin Yohannan
"""

import gym
import assistive_gym
import pybullet as p
import numpy as np
from PIL import Image
import utils
from pid import PID
from time import sleep
import time
import cv2
import dlib
import math

end_effector_id = 9
# States
FIND_POSE = 0
MASK_ON = 1
MASK_OFF = 2
DONE = 3

detector = dlib.get_frontal_face_detector()
predictor = dlib.shape_predictor('./shape_predictor_68_face_landmarks.dat')

def Start_video(env):
    detector = dlib.get_frontal_face_detector()
    predictor = dlib.shape_predictor('./shape_predictor_68_face_landmarks.dat')

    while True:
        frame = env.get_camera_frame()
        gray = cv2.cvtColor(frame, cv2.COLOR_RGBA2GRAY)

        ## Facial detection
        faces = detector(gray)
        # print(len(faces)) # =1

        if len(faces) == 0:
            print('No face found')
            continue

        for face in faces:
            x1 = face.left()
            y1 = face.top()
            x2 = face.right()
            y2 = face.bottom()
            cv2.rectangle(frame, (x1, y1), (x2, y2), (0, 255, 0), 1)  # 1 pixel thickness

            landmarks = predictor(gray, face)

            '''
            ### Find the ears upper and lower region
            #for n in range(0, 68):
            for n in [0, 2, 14, 16]:      ### Ears = [0,2,14,16]
                x = landmarks.part(n).x
                y = landmarks.part(n).y
                cv2.circle(frame, (x,y), 3, (255, 0, 0), -1) #radius=3, fill the circle = -1
            '''

        ### Drawing a line from point 28 to 8.
        x28 = landmarks.part(28).x
        y28 = landmarks.part(28).y
        x8 = landmarks.part(8).x
        y8 = landmarks.part(8).y

        ### Calculate angle between the 2 points
        angle = math.atan2((y8 - y28), (x8 - x28)) * 180 / 3.141592653;
        angle = angle + math.ceil(-angle / 360) * 360

        allx = []
        ally = []
        if angle > 85 and angle < 95:

            x = landmarks.part(66).x
            y = landmarks.part(66).y

            allx.append(x)
            ally.append(y)

            cv2.putText(frame, "Perfect, stay still!, mouth location = ({},{})".format(x, y),
                        (0, 25), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0, 255, 0), 2)

            cv2.circle(frame, (x, y), 4, (255, 0, 0), -1)  # radius=3, fill the circle = -1
            # return x, y

        else:
            ### Facial Center points forehead, lower lip center and chin. Points 28, 57 and 8.
            for n in [28, 57, 8]:
                x = landmarks.part(n).x
                y = landmarks.part(n).y
                cv2.circle(frame, (x, y), 4, (255, 0, 0), -1)  # radius=3, fill the circle = -1

            cv2.line(frame, (x28, y28), (x8, y8), (0, 0, 0), 2)

            ### Tell the user to adjust the face orientation
            # cv2.putText(img,'Hello World!', bottomLeftCornerOfText, font, fontScale,fontColor,lineType)
            cv2.putText(frame, 'Adjust face to 90 Deg, current angle = {}'.format(angle),
                        (0, 25), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0, 255, 0), 2)

        cv2.imshow("Window", frame)

        # This breaks on 'q' key
        if cv2.waitKey(1) & 0xFF == ord('q'):
            break

    cv2.destroyAllWindows()

    return math.floor(np.mean(allx)), math.floor(np.mean(ally))

def assisted_ctrl_loop(p, env, target_pos, target_orient, pid_ctrl):
    env.render()
    # IK to get new joint positions (angles) for the robot
    target_joint_positions = p.calculateInverseKinematics(env.robot, end_effector_id, target_pos, target_orient)
    target_joint_positions = target_joint_positions[:7]
    # Get the joint positions (angles) of the robot arm
    joint_positions, joint_velocities, joint_torques = env.get_motor_joint_states(env.robot)
    joint_positions = np.array(joint_positions)[:7]
    # Set joint action to be the error between current and target joint positions
    err = target_joint_positions - joint_positions
    u = pid_ctrl.update(err)
    observation, reward, done, info = env.step(u)

def raw_ctrl_loop(p, env, target_pos, target_orient, reduce_force=True):
    p.stepSimulation()
    joint_poses = p.calculateInverseKinematics(env.robot, end_effector_id, target_pos, target_orient)
    target_joint_positions = joint_poses[:7]

    joint_positions, joint_velocities, joint_torques = env.get_motor_joint_states(env.robot)
    joint_positions = np.array(joint_positions)[:7]

    if reduce_force:
        distance = np.linalg.norm(target_joint_positions - joint_positions)
        force_applied = distance * 0.25
    else:
        force_applied = 1.0

    p.setJointMotorControlArray(bodyIndex=env.robot,
                                jointIndices=list(range(1, 8)),
                                controlMode=p.POSITION_CONTROL,
                                targetPositions=joint_poses,
                                targetVelocities=np.zeros(len(joint_poses)),
                                forces=force_applied * np.ones(len(joint_poses)),
                                positionGains=0.03 * np.ones(len(joint_poses)),
                                velocityGains=np.ones(len(joint_poses)))
    sleep(0.05)

def run():
    env = gym.make("MaskPlacingJaco-v0")

    # Set head orientation
    env.set_head_orient(0.0, 0.0, -30.0)

    # Set robot base position (relative to wheelchair)
    env.set_robot_base([-0.35, -0.3, 0.3])  # x, y, z in meters

    # Set camera position
    # env.set_camera_position([0.25, -0.5, 0.75], [0.0, 0.0, 1.25])
    env.set_camera_position([0.0, -0.5, 1.25], [0.0, 0.0, 1.25])

    env.render()
    observation = env.reset()
    utils.print_joints(env, p)
    # p.addUserDebugLine(lineFromXYZ=[0.0, 0.0, 0.0], lineToXYZ=[0.0, 0.0, 2.0], lineColorRGB=[255,0,0])
    # p.addUserDebugLine(lineFromXYZ=[0.25, -0.5, 0.0], lineToXYZ=[0.25, -0.5, 0.75], lineColorRGB=[0,0,255])

    q_key = ord('q')
    r_key = ord('r')
    f_key = ord('f')
    b_key = ord('b')

    # FSM state variable
    curr_state = FIND_POSE

    target_pos = env.target_pos
    target_orient = env.target_orient

    start_pos, start_orient, _, _, _, _ = p.getLinkState(env.robot, linkIndex=end_effector_id)

    pid_ctrl = PID(kp=3.0, ki=0.0, kd=0.0, q_dim=7)

    # Main loop
    while True:
        if curr_state == FIND_POSE:
            # frame = env.get_camera_frame()
            # img = Image.fromarray(frame, 'RGBA')
            # img.save('head.png')
            # # x, y = Start_video(env)
            # # print(x, y)
            # curr_state = MASK_ON
            pass

        elif curr_state == MASK_ON:
            assisted_ctrl_loop(p, env, target_pos, target_orient, pid_ctrl)
            # raw_ctrl_loop(p, env, target_pos, target_orient, reduce_force=True)

        elif curr_state == MASK_OFF:
            # assisted_ctrl_loop(p, env, start_pos, start_orient, pid_ctrl)
            raw_ctrl_loop(p, env, start_pos, start_orient, reduce_force=False)

        else:
            pass

        # State transition logic
        keys = p.getKeyboardEvents()
        if q_key in keys and keys[q_key] & p.KEY_WAS_TRIGGERED:  # Quit
            break

        elif f_key in keys and keys[f_key] & p.KEY_WAS_TRIGGERED:  # Forward
            curr_state = MASK_ON

        elif b_key in keys and keys[b_key] & p.KEY_WAS_TRIGGERED:  # Backward
            curr_state = MASK_OFF

        elif r_key in keys and keys[r_key] & p.KEY_WAS_TRIGGERED:  # Reset Env
            env.set_head_orient(0.0, 30.0, 0.0)
            env.reset()
            target_pos = env.target_pos
            target_orient = env.target_orient
            curr_state = FIND_POSE

    p.disconnect()

if __name__ == "__main__":
    run()

max_stars_repo_name: xinchaosong/mask-placement-apparatus

max_issues_repo_name: xinchaosong/mask-placement-apparatus

max_forks_repo_name: xinchaosong/mask-placement-apparatus


text:
import nengo
import numpy as np
import pandas as pd
import scipy.linalg
from scipy.special import legendre
import pickle

# This code is completely taken from Terry Steward:
class DiscreteDelay(nengo.synapses.Synapse):
    def __init__(self, delay, size_in=1):
        self.delay = delay
        super().__init__(default_size_in=size_in, default_size_out=size_in)

    def make_state(self, shape_in, shape_out, dt, dtype=None, y0=None):
        return {}

    def make_step(self, shape_in, shape_out, dt, rng, state=None):
        steps = int(self.delay/dt)
        if steps == 0:
            def step_delay(t, x):
                return x
            return step_delay
        assert steps > 0

        state = np.zeros((steps, shape_in[0]))
        state_index = np.array([0])

        def step_delay(t, x, state=state, state_index=state_index):
            result = state[state_index]
            state[state_index] = x
            state_index[:] = (state_index + 1) % state.shape[0]
            return result

        return step_delay

# This code is completely taken from Terry Steward:
class LDN(nengo.Process):
    def __init__(self, theta, q, size_in=1):
        self.q = q  # number of internal state dimensions per input
        self.theta = theta  # size of time window (in seconds)
        self.size_in = size_in  # number of inputs

        # Do Aaron's math to generate the matrices A and B so that
        #  dx/dt = Ax + Bu will convert u into a legendre representation over a window theta
        #  https://github.com/arvoelke/nengolib/blob/master/nengolib/synapses/analog.py#L536
        A = np.zeros((q, q))
        B = np.zeros((q, 1))
        for i in range(q):
            B[i] = (-1.) ** i * (2 * i + 1)
            for j in range(q):
                A[i, j] = (2 * i + 1) * (-1 if i < j else (-1.) ** (i - j + 1))
        self.A = A / theta
        self.B = B / theta

        super().__init__(default_size_in=size_in, default_size_out=q * size_in)

    def make_step(self, shape_in, shape_out, dt, rng, state=None):
        state = np.zeros((self.q, self.size_in))

        # Handle the fact that we're discretizing the time step
        #  https://en.wikipedia.org/wiki/Discretization#Discretization_of_linear_state_space_models
        Ad = scipy.linalg.expm(self.A * dt)
        Bd = np.dot(np.dot(np.linalg.inv(self.A), (Ad - np.eye(self.q))), self.B)

        # this code will be called every timestep
        def step_legendre(t, x, state=state):
            state[:] = np.dot(Ad, state) + np.dot(Bd, x[None, :])
            return state.T.flatten()

        return step_legendre

    def get_weights_for_delays(self, r):
        # compute the weights needed to extract the value at time r
        # from the network (r=0 is right now, r=1 is theta seconds ago)
        r = np.asarray(r)
        m = np.asarray([legendre(i)(2 * r - 1) for i in range(self.q)])
        return m.reshape(self.q, -1).T

class DataParser:
    """This node will parse data from a pandas dataframe object as input.
    """
    def __init__(self, data_df, sample_freq, vars=[]):
        self.data_df = data_df
        self.sample_freq = sample_freq
        self.vars = vars

    def parse_data(self, t):
        r = [self.data_df[x].iloc[int(t * self.sample_freq)] for x in self.vars]
        return r

    def update_data(self, data_df):
        self.data_df = data_df

class SwitchNode:
    """This node will return the first half of its input for t <= t_switch and
    the second half of its input for t > t_switch.
    By default, the output will be clipped between -1 and 1.
    """
    def __init__(self, t_switch=4.0, stim_size=4, clip=True):
        self.t_switch = t_switch
        self.stim_size = stim_size
        self.clip=clip

    def step(self, t, x):
        if t <= self.t_switch:
            state = x[:self.stim_size]
        else:
            state = x[self.stim_size:]

        if self.clip:
            state = np.clip(state, -1, 1)

        return state

class ErrorSwitchNode:
    """This node will compute the error between a prediction and ground truth.
    The output of the node will be 0 for t < t_init and t > t_switch.
    """
    def __init__(self, t_init=0.1, t_switch=4.0, stim_size=4, error="normal"):
        self.t_init = t_init
        self.t_switch = t_switch
        self.stim_size = stim_size
        self.error = error
        if error.lower() == "normal":
            self.error_func = self.normal
        elif error.lower() == "mse":
            self.error_func = self.mse
        else:
            raise ValueError('This error type is not supported. Please use "mse" or "normal".')

    def mse(self, a, b):
        return np.sign(a - b) * (a - b) ** 2

    def normal(self, a, b):
        return a - b

    def step(self, t, x):
        if self.t_init <= t <= self.t_switch:
            return self.error_func(x[:self.stim_size], x[self.stim_size:])
        return np.zeros(self.stim_size)

class PredictiveModelAutoregressiveLMU:
    def __init__(self, seed=42, neurons_per_dim=100, sample_freq=50,
                 lmu_theta=0.1, lmu_q=1, radius=1.5, dt=0.001,
                 t_delay=0.02, learning_rate=5e-5, action_vars=["Q"],
                 state_vars=["angle_sin", "angle_cos", "angleD", "position", "positionD"],
                 action_df=None, state_df=None, weights=None, scales={},
                 predict_delta=True, error="normal", t_switch=4.0, t_init=0.1,
                 model_name="AutoregressiveLMU", *args, **kwargs):

        self.seed = seed
        self.neurons_per_dim = neurons_per_dim
        self.sample_freq = sample_freq
        self.t_delay = t_delay
        self.learning_rate = learning_rate
        self.action_vars = action_vars
        self.action_dim = len(action_vars)
        self.state_vars = state_vars
        self.state_dim = len(state_vars)
        self.radius = radius
        self.dt = dt
        self.lmu_q = lmu_q
        self.lmu_theta = lmu_theta
        self.weights = weights
        self.scales = scales
        self.predict_delta = predict_delta
        self.t_switch = t_switch
        self.t_init = t_init
        self.error = error
        self.model_name = model_name

        if action_df is None:
            self.action_df = pd.DataFrame(
                np.zeros((1, len(action_vars) + 1)),
                columns=["time"] + action_vars,
            )

        if state_df is None:
            self.state_df = pd.DataFrame(
                np.zeros((1, len(state_vars) + 1)),
                columns=["time"] + state_vars,
            )

        self.action_parser = DataParser(
            data_df=self.action_df,
            sample_freq=self.sample_freq,
            vars=self.action_vars
        )

        self.state_parser = DataParser(
            data_df=self.state_df,
            sample_freq=self.sample_freq,
            vars=self.state_vars
        )

        # this function streams the state signal from file to node
        def action_stim_func(t):
            return self.action_parser.parse_data(t)

        # this function streams the state signal from file to node
        def state_stim_func(t):
            return self.state_parser.parse_data(t)

        self.action_stim_func = action_stim_func
        self.state_stim_func = state_stim_func

        self.model, self.recordings = self.make_model()
        self.sim = nengo.Simulator(self.model, dt=self.dt, progress_bar=False)

    def set_inputs(self, action_df, state_df):

        for v in self.action_vars + ["time"]:
            assert v in action_df.columns
        for v in self.state_vars:
            assert v in state_df.columns

        self.action_df = action_df
        self.state_df = state_df

        self.action_parser.update_data(self.action_df)
        self.state_parser.update_data(self.state_df)

    def reset_sim(self):

        self.sim.reset(seed=self.seed)

    def set_weights(self, weights):

        weights = np.array(weights)
        assert weights.shape == (
            self.state_dim,
            self.neurons_per_dim * (self.state_dim + self.action_dim) * (1 + self.lmu_q)
        )

        self.weights = weights
        self.sim.signals[self.sim.model.sig[self.learned_connection]["weights"]] = self.weights

    def get_weights(self):

        weights = self.sim.signals[self.sim.model.sig[self.learned_connection]["weights"]]
        return np.array(weights)

    def process_files(self):

        t_max = self.action_df["time"].max()  # number of seconds to run
        self.sim.run(t_max)

        return self.recordings

    def get_state_dict(self):

        state_dict = {
            "seed": self.seed,
            "neurons_per_dim": self.neurons_per_dim,
            "sample_freq": self.sample_freq,
            "t_delay": self.t_delay,
            "learning_rate": self.learning_rate,
            "action_vars": self.action_vars,
            "state_vars": self.state_vars,
            "radius": self.radius,
            "dt": self.dt,
            "lmu_q": self.lmu_q,
            "lmu_theta": self.lmu_theta,
            "weights": self.get_weights(),
            "scales": self.scales,
            "predict_delta": self.predict_delta,
            "error": self.error,
            "t_switch": self.t_switch,
            "t_init": self.t_init,
            "model_name": self.model_name
        }

        return state_dict

    def save_state_dict(self, path="model_state.pkl"):

        state_dict = self.get_state_dict()
        with open(path, "wb") as f:
            pickle.dump(state_dict, f)

    def make_model(self):

        if self.weights is None:
            self.weights = np.zeros((
                self.state_dim,
                self.neurons_per_dim * (self.state_dim + self.action_dim) * (1 + self.lmu_q))
            )

        model = nengo.Network()
        with model:
            # set the default synapse to None (normal default is 0.005s)
            model.config[nengo.Connection].synapse = None

            # initialize input nodes
            true_action = nengo.Node(self.action_stim_func)
            true_state = nengo.Node(self.state_stim_func)

            # create a node that first outputs the true state and then switches to predicted state
            switch = SwitchNode(t_switch=self.t_switch, stim_size=self.state_dim)
            believed_state = nengo.Node(switch.step, size_in=self.state_dim*2, size_out=self.state_dim)
            nengo.Connection(true_state, believed_state[:self.state_dim])

            # make a node for the predicted future state
            predicted_future_state = nengo.Node(None, size_in=self.state_dim)

            # make LMU unit
            ldn = nengo.Node(LDN(theta=self.lmu_theta, q=self.lmu_q, size_in=self.state_dim + self.action_dim))
            nengo.Connection(true_action, ldn[:self.action_dim])
            nengo.Connection(believed_state, ldn[self.action_dim:])

            # make the hidden layer
            ens = nengo.Ensemble(
                n_neurons=self.neurons_per_dim * (self.state_dim + self.action_dim)*(1+self.lmu_q),
                dimensions=(self.state_dim + self.action_dim)*(1+self.lmu_q),
                neuron_type=nengo.LIFRate(),
                seed=self.seed
            )
            nengo.Connection(true_action, ens[:self.action_dim])
            nengo.Connection(believed_state, ens[self.action_dim:self.action_dim+self.state_dim])
            nengo.Connection(ldn, ens[self.action_dim+self.state_dim:])

            # if wanted, have the model predict the difference from the last state instead
            if self.predict_delta:
                print("predicting only the difference between the current and next state")
                nengo.Connection(believed_state, predicted_future_state)

            # make the output weights we can learn
            self.learned_connection = nengo.Connection(
                ens.neurons,
                predicted_future_state,
                transform=self.weights,  # change this if you have pre-recorded weights to use
                seed=self.seed,
                learning_rule_type=nengo.PES(
                    learning_rate=self.learning_rate,
                    pre_synapse=DiscreteDelay(self.t_delay)  # delay the activity value when updating weights
                )
            )

            # this is what the network predicted the current state to be in the past
            predicted_current_state = nengo.Node(None, size_in=self.state_dim)
            nengo.Connection(predicted_future_state, predicted_current_state, synapse=DiscreteDelay(self.t_delay))
            nengo.Connection(predicted_current_state, believed_state[self.state_dim:])

            # compute the error by subtracting the current measurement from a delayed version of the prediction
            error_node = ErrorSwitchNode(
                t_init=self.t_init,
                t_switch=self.t_switch,
                stim_size=self.state_dim,
                error=self.error
            )
            prediction_error = nengo.Node(error_node.step, size_in=self.state_dim*2, size_out=self.state_dim)
            nengo.Connection(predicted_current_state, prediction_error[:self.state_dim])
            nengo.Connection(true_state, prediction_error[self.state_dim:])

            # apply the error to the learning rule
            nengo.Connection(prediction_error, self.learned_connection.learning_rule)

            recordings = {
                "states" : nengo.Probe(true_state),
                "actions": nengo.Probe(true_action),
                "delay": self.t_delay,
                "predicted_current_states": nengo.Probe(predicted_current_state),
                "predicted_future_states": nengo.Probe(predicted_future_state),
                "prediction_errors": nengo.Probe(prediction_error),
            }

        return model, recordings

reason: import numpy,import scipy,from scipy


text:
def sizeEffect(listManual, listAuto):
    import numpy as np
    import math as m
    listM = [item for sublist in listManual for item in sublist]
    listA = [item for sublist in listAuto for item in sublist]
    diff = np.array(listA) - np.array(listM)
    meanD = np.mean(diff)
    stdD = np.std(diff)
    #Cohen size effect
    size_effect = meanD/stdD
    cohen_d = (np.mean(listA)-np.mean(listM)) / (m.sqrt((np.std(listA) ** 2 + np.std(listM) ** 2) / 2))
    return size_effect, cohen_d

def ttests(path_to_dict, name_dict='\\TOTresults'):
    import numpy as np
    import fnmatch
    from dictmanager import load_obj

    #NOTES TO UNDERSTAND NOTATIONS
    #d = list of distances from reference point
    #fl = list of FL
    #PAs = list of PA sup
    #PAi = list of PA inf
    #mt = list of MT
    #_s : simple images
    #_p : panoramic images
    #_m : manual
    #_a : automated
    #_filtered: matched fascicles only

    participants=['01_Kevin', '02_rafaelopes', '03_charlesbarrand', '04_guilhem',\
        '05_leandre', '06_thomasmartine', '10_victor',\
        '11_youssouf', '12_sufyan', '16_julien', '34_nicolas']

    '************************************************************************'
    '*****************************INITIALIZATION*****************************'

    d_s_m = [[] for par in range(len(participants))]
    mt_s_m = [[] for par in range(len(participants))]
    d_s_m_filtered = [[] for par in range(len(participants))]
    fl_s_m_filtered = [[] for par in range(len(participants))]
    PAs_s_m_filtered = [[] for par in range(len(participants))]
    PAi_s_m_filtered = [[] for par in range(len(participants))]

    d_s_a = [[] for par in range(len(participants))]
    mt_s_a = [[] for par in range(len(participants))]
    d_s_a_filtered = [[] for par in range(len(participants))]
    fl_s_a_filtered = [[] for par in range(len(participants))]
    PAs_s_a_filtered = [[] for par in range(len(participants))]
    PAi_s_a_filtered = [[] for par in range(len(participants))]

    d_p_m = [[] for par in range(len(participants))]
    mt_p_m = [[] for par in range(len(participants))]
    d_p_m_filtered = [[] for par in range(len(participants))]
    fl_p_m_filtered = [[] for par in range(len(participants))]
    PAs_p_m_filtered = [[] for par in range(len(participants))]
    PAi_p_m_filtered = [[] for par in range(len(participants))]

    d_p_a = [[] for par in range(len(participants))]
    mt_p_a = [[] for par in range(len(participants))]
    d_p_a_filtered = [[] for par in range(len(participants))]
    fl_p_a_filtered = [[] for par in range(len(participants))]
    PAs_p_a_filtered = [[] for par in range(len(participants))]
    PAi_p_a_filtered = [[] for par in range(len(participants))]

    #stats on the number of fascicles detected
    nb_fasc_tot_s = 0
    nb_fasc_in_s = 0
    nb_fasc_filt_s = 0
    nb_images_s = 0
    nb_fasc_tot_p = 0
    nb_fasc_in_p = 0
    nb_fasc_filt_p = 0
    nb_images_p = 0

    '************************************************************************'
    '*****************************DATA RETRIEVAL*****************************'

    dictio = load_obj(name_dict, path_to_dict)
    l2 = ['fasc*', 'fsc_*']

    for par in range(len(participants)):

        participant = participants[par]
        fam_folders = [str(d) for d in dictio[participant].keys()]

        s_manuFasc = []
        s_autoFasc = []
        p_manuFasc = []
        p_autoFasc = []

        for fam in fam_folders:

            ###################################################################
            # simple images
            dictioS = dictio[participant][fam]['BF']['simple']
            images = [str(im) for im in dictioS.keys()]
            for i in images:
                # if par == 9 and fam =='fam_2' and i=='img_2':
                #     print(par, fam, i)
                # else:
                nb_images_s = nb_images_s + 1

                ###############################################################
                # SIMPLE - manual
                dictioM = dictioS[i]['architecture manual']
                fascicles = [str(fa) for fa in dictioM if any(fnmatch.fnmatch(fa, p) for p in l2)]
                for f in fascicles:
                    dictioF = dictioM[f]
                    idf = fam + '/' + i + '/' + f
                    if len(dictioF.keys())>1:
                        s_manuFasc.append((idf, dictioF['dist from (0,0) of RGB image, in mm']))
                        d_s_m[par].append(dictioF['dist from (0,0) of RGB image, in mm'])
                ###############################################################
                # SIMPLE - automatic
                if ('architecture auto' in dictioS[i]):
                    dictioA = dictioS[i]['architecture auto']
                    midRow = np.mean(dictioA['crop']['lines'])
                    midCol = np.mean(dictioA['crop']['columns'])
                    if dictioA and ('MT' in dictioA):
                        fascicles = [fa for fa in dictioA if any(fnmatch.fnmatch(fa, p) for p in l2)]
                        nb_fasc_tot_s = nb_fasc_tot_s + len(fascicles)
                        for f in fascicles:
                            dictioF = dictioA[f]
                            idf = fam + '/' + i + '/' + f
                            if len(dictioF.keys())>1:
                                #keep the fascicles that are in the lower half of the image,
                                #to compare with manual data - often taken in that region
                                PAi = dictioF['PAinf']['intersection with apo']
                                PAs = dictioF['PAsup']['intersection with apo']
                                fasc_row = (PAs[0]-PAi[0])/(PAs[1]-PAi[1])*(midCol-PAs[1])+PAs[0]

                                if fasc_row <= midRow:
                                    s_autoFasc.append((idf, dictioF['dist from (0,0) of RGB image, in mm']))
                                    d_s_a[par].append(dictioF['dist from (0,0) of RGB image, in mm'])
                                    nb_fasc_in_s = nb_fasc_in_s + 1

                        if ('MT for labelled points' in dictioM['MT']):
                            for ind0 in range(len(dictioM['MT']['MT for labelled points'])):
                                elem = dictioM['MT']['MT for labelled points'][ind0]
                                if elem != 'error':
                                    mt_s_m[par].append(elem) #MT in mm

                            for ind0 in range(len(dictioA['MT']['MT for labelled points'])):
                                elem = dictioA['MT']['MT for labelled points'][ind0]
                                if elem != 'error':
                                    mt_s_a[par].append(elem)

            ###################################################################
            # panoramic images
            dictioP = dictio[participant][fam]['BF']['panoramic']
            images = [str(im) for im in dictioP.keys()]
            for i in images:
                nb_images_p = nb_images_p + 1

                ###############################################################
                # PANORAMIC - manual
                dictioM = dictioP[i]['architecture manual']
                fascicles = [fa for fa in dictioM if any(fnmatch.fnmatch(fa, p) for p in l2)]
                for f in fascicles:
                    dictioF = dictioM[f]
                    idf = fam + '/' + i + '/' + f
                    if len(dictioF.keys())>1:
                        p_manuFasc.append((idf, dictioF['dist from insertion in mm']))
                        d_p_m[par].append(dictioF['dist from insertion in mm'])

                ###############################################################
                # PANORAMIC - automatic
                if ('architecture auto' in dictioP[i]):
                    dictioA = dictioP[i]['architecture auto']
                    if dictioA and ('MT' in dictioA):
                        fascicles = [fa for fa in dictioA if any(fnmatch.fnmatch(fa, p) for p in l2)]
                        nb_fasc_tot_p = nb_fasc_tot_p + len(fascicles)
                        for f in fascicles:
                            dictioF = dictioA[f]
                            idf = fam + '/' + i + '/' + f
                            #only keep fascicles that are entirely within the cropped image,
                            #to compare with manually identified fascicles
                            if len(dictioF.keys())>1 and dictioF['FL']['in/out of the image'] == 'in image':
                                nb_fasc_in_p = nb_fasc_in_p + 1
                                p_autoFasc.append((idf, dictioF['dist from insertion in mm']))
                                d_p_a[par].append(dictioF['dist from insertion in mm'])

                        if ('MT for labelled points' in dictioM['MT']):
                            for ind0 in range(len(dictioM['MT']['MT for labelled points'])):
                                elem = dictioM['MT']['MT for labelled points'][ind0]
                                if elem != 'error':
                                    mt_p_m[par].append(elem) #MT in mm

                            for ind0 in range(len(dictioA['MT']['MT for labelled points'])):
                                elem = dictioA['MT']['MT for labelled points'][ind0]
                                if elem != 'error':
                                    mt_p_a[par].append(elem)

        '************************************************************************'
        '********************MATCHING AUTO & MANUAL FASCICLES*******************'

        listePair_manuF_s = []
        for n in range(len(s_manuFasc)):
            mf = s_manuFasc[n]
            subtr = [(tup,abs(tup[1]- mf[1])) for tup in s_autoFasc]
            subtr.sort(key=lambda x:x[1])
            closest = subtr[0]
            listePair_manuF_s.append((mf[0], closest[0][0], closest[1])) #tuple = ( ID manu fasc, ID auto fasc, distance entre les deux)
        listePair_manuF_s.sort(key=lambda x:x[1])
        uniqueMatching = []
        counterL = 0
        while counterL < len(listePair_manuF_s):
            currentAutoFasc = listePair_manuF_s[counterL][1]
            correspondingAutoFasc = [(listePair_manuF_s[counterL][0], listePair_manuF_s[counterL][2])]
            rank = counterL + 1
            while rank<len(listePair_manuF_s) and listePair_manuF_s[rank][1] == currentAutoFasc:
                correspondingAutoFasc.append((listePair_manuF_s[rank][0],listePair_manuF_s[rank][2]))
                rank = rank + 1
            correspondingAutoFasc.sort(key=lambda x:x[1])
            uniqueMatching.append((correspondingAutoFasc[0][0], currentAutoFasc, correspondingAutoFasc[0][1]))
            counterL = rank
        for element in uniqueMatching:
            pathA = element[1].split('/')
            pathM = element[0].split('/')
            nb_fasc_filt_s = nb_fasc_filt_s + 1
            d_s_m_filtered[par].append(dictio[participant][pathM[0]]['BF']['simple'][pathM[1]]['architecture manual'][pathM[2]]['dist from (0,0) of RGB image, in mm'])
            fl_s_m_filtered[par].append(dictio[participant][pathM[0]]['BF']['simple'][pathM[1]]['architecture manual'][pathM[2]]['FL']['length in mm'])
            PAs_s_m_filtered[par].append(dictio[participant][pathM[0]]['BF']['simple'][pathM[1]]['architecture manual'][pathM[2]]['PAsup']['value in degree'])
            PAi_s_m_filtered[par].append(dictio[participant][pathM[0]]['BF']['simple'][pathM[1]]['architecture manual'][pathM[2]]['PAinf']['value in degree'])
            d_s_a_filtered[par].append(dictio[participant][pathA[0]]['BF']['simple'][pathA[1]]['architecture auto'][pathA[2]]['dist from (0,0) of RGB image, in mm'])
            fl_s_a_filtered[par].append(dictio[participant][pathA[0]]['BF']['simple'][pathA[1]]['architecture auto'][pathA[2]]['FL']['length in mm'])
            PAs_s_a_filtered[par].append(dictio[participant][pathA[0]]['BF']['simple'][pathA[1]]['architecture auto'][pathA[2]]['PAsup']['value in degree'])
            PAi_s_a_filtered[par].append(dictio[participant][pathA[0]]['BF']['simple'][pathA[1]]['architecture auto'][pathA[2]]['PAinf']['value in degree'])

        listePair_manuF_p = []
        for n in range(len(p_manuFasc)):
            mf = p_manuFasc[n]
            subtr = [(tup,abs(tup[1]- mf[1])) for tup in p_autoFasc]
            subtr.sort(key=lambda x:x[1])
            closest = subtr[0]
            listePair_manuF_p.append((mf[0], closest[0][0], closest[1])) #tuple = ( ID manu fasc, ID auto fasc, distance entre les deux)
        listePair_manuF_p.sort(key=lambda x:x[1])
        uniqueMatching = []
        counterL = 0
        while counterL < len(listePair_manuF_p):
            currentAutoFasc = listePair_manuF_p[counterL][1]
            correspondingAutoFasc = [(listePair_manuF_p[counterL][0], listePair_manuF_p[counterL][2])]
            rank = counterL + 1
            while rank<len(listePair_manuF_p) and listePair_manuF_p[rank][1] == currentAutoFasc:
                correspondingAutoFasc.append((listePair_manuF_p[rank][0],listePair_manuF_p[rank][2]))
                rank = rank + 1
            correspondingAutoFasc.sort(key=lambda x:x[1])
            uniqueMatching.append((correspondingAutoFasc[0][0], currentAutoFasc, correspondingAutoFasc[0][1]))
            counterL = rank
        for element in uniqueMatching:
            pathA = element[1].split('/')
            pathM = element[0].split('/')
            nb_fasc_filt_p = nb_fasc_filt_p + 1
            d_p_m_filtered[par].append(dictio[participant][pathM[0]]['BF']['panoramic'][pathM[1]]['architecture manual'][pathM[2]]['dist from insertion in mm'])
            fl_p_m_filtered[par].append(dictio[participant][pathM[0]]['BF']['panoramic'][pathM[1]]['architecture manual'][pathM[2]]['FL']['length in mm'])
            PAs_p_m_filtered[par].append(dictio[participant][pathM[0]]['BF']['panoramic'][pathM[1]]['architecture manual'][pathM[2]]['PAsup']['value in degree'])
            PAi_p_m_filtered[par].append(dictio[participant][pathM[0]]['BF']['panoramic'][pathM[1]]['architecture manual'][pathM[2]]['PAinf']['value in degree'])
            d_p_a_filtered[par].append(dictio[participant][pathA[0]]['BF']['panoramic'][pathA[1]]['architecture auto'][pathA[2]]['dist from insertion in mm'])
            fl_p_a_filtered[par].append(dictio[participant][pathA[0]]['BF']['panoramic'][pathA[1]]['architecture auto'][pathA[2]]['FL']['length in mm'])
            PAs_p_a_filtered[par].append(dictio[participant][pathA[0]]['BF']['panoramic'][pathA[1]]['architecture auto'][pathA[2]]['PAsup']['value in degree'])
            PAi_p_a_filtered[par].append(dictio[participant][pathA[0]]['BF']['panoramic'][pathA[1]]['architecture auto'][pathA[2]]['PAinf']['value in degree'])

    #t_tests
    print('paired samples t-tests resuts: ')
    from scipy.stats.mstats import ttest_rel
    #NOTES: we cannot user '..._filtered' arrays directly because of their structure
    #we need to flatten them to 1-D lists
    t,p=ttest_rel([item for sublist in PAs_s_m_filtered for item in sublist],[item for sublist in PAs_s_a_filtered for item in sublist],axis=None)
    print('PAS s',p)
    t2,p2=ttest_rel([item for sublist in PAs_p_m_filtered for item in sublist],[item for sublist in PAs_p_a_filtered for item in sublist],axis=None)
    print('PAS p',p2)
    t3,p3=ttest_rel([item for sublist in PAi_s_m_filtered for item in sublist],[item for sublist in PAi_s_a_filtered for item in sublist],axis=None)
    print('PAI s',p3)
    t4,p4=ttest_rel([item for sublist in PAi_p_m_filtered for item in sublist],[item for sublist in PAi_p_a_filtered for item in sublist],axis=None)
    print('PAI p',p4)
    t5,p5=ttest_rel([item for sublist in fl_s_m_filtered for item in sublist],[item for sublist in fl_s_a_filtered for item in sublist],axis=None)
    print('FL s',p5)
    t6,p6=ttest_rel([item for sublist in fl_p_m_filtered for item in sublist],[item for sublist in fl_p_a_filtered for item in sublist],axis=None)
    print('FL p',p6)
    t7,p7=ttest_rel([item for sublist in mt_s_m for item in sublist],[item for sublist in mt_s_a for item in sublist],axis=None)
    print('mt s',p7)
    t7_2,p7_2=ttest_rel([np.mean(sublist) for sublist in mt_s_m],[np.mean(sublist) for sublist in mt_s_a],axis=None)
    print('mt s for means',p7_2)
    t8,p8=ttest_rel([item for sublist in mt_p_m for item in sublist],[item for sublist in mt_p_a for item in sublist],axis=None)
    print('mt p',p8)
    t8_2,p8_2=ttest_rel([np.mean(sublist) for sublist in mt_p_m],[np.mean(sublist) for sublist in mt_p_a],axis=None)
    print('mt p for means',p8_2)

    print('independent samples t-tests resuts: ')
    from scipy.stats.mstats import ttest_rel, ttest_ind
    #NOTES: we cannot user '..._filtered' arrays directly because of their structure
    #we need to flatten them to 1-D lists
    t,p=ttest_ind([item for sublist in PAs_s_m_filtered for item in sublist],[item for sublist in PAs_s_a_filtered for item in sublist],axis=None)
    print('PAS s',p)
    t2,p2=ttest_ind([item for sublist in PAs_p_m_filtered for item in sublist],[item for sublist in PAs_p_a_filtered for item in sublist],axis=None)
    print('PAS p',p2)
    t3,p3=ttest_ind([item for sublist in PAi_s_m_filtered for item in sublist],[item for sublist in PAi_s_a_filtered for item in sublist],axis=None)
    print('PAI s',p3)
    t4,p4=ttest_ind([item for sublist in PAi_p_m_filtered for item in sublist],[item for sublist in PAi_p_a_filtered for item in sublist],axis=None)
    print('PAI p',p4)
    t5,p5=ttest_ind([item for sublist in fl_s_m_filtered for item in sublist],[item for sublist in fl_s_a_filtered for item in sublist],axis=None)
    print('FL s',p5)
    t6,p6=ttest_ind([item for sublist in fl_p_m_filtered for item in sublist],[item for sublist in fl_p_a_filtered for item in sublist],axis=None)
    print('FL p',p6)
    t7,p7=ttest_ind([item for sublist in mt_s_m for item in sublist],[item for sublist in mt_s_a for item in sublist],axis=None)
    print('mt s',p7)
    t7_2,p7_2=ttest_ind([np.mean(sublist) for sublist in mt_s_m],[np.mean(sublist) for sublist in mt_s_a],axis=None)
    print('mt s for means',p7_2)
    t8,p8=ttest_ind([item for sublist in mt_p_m for item in sublist],[item for sublist in mt_p_a for item in sublist],axis=None)
    print('mt p',p8)
    t8_2,p8_2=ttest_ind([np.mean(sublist) for sublist in mt_p_m],[np.mean(sublist) for sublist in mt_p_a],axis=None)
    print('mt p for means',p8_2)

    #size effects
    s1=sizeEffect(PAs_s_m_filtered,PAs_s_a_filtered)
    s2=sizeEffect(PAs_p_m_filtered,PAs_p_a_filtered)
    s3=sizeEffect(PAi_s_m_filtered,PAi_s_a_filtered)
    s4=sizeEffect(PAi_p_m_filtered,PAi_p_a_filtered)
    s5=sizeEffect(fl_s_m_filtered,fl_s_a_filtered)
    s6=sizeEffect(fl_p_m_filtered,fl_p_a_filtered)
    s7=sizeEffect(mt_s_m,mt_s_a)
    s8=sizeEffect(mt_p_m,mt_p_a)
    print('Size effects: ')
    print('PAS s',s1)
    print('PAS p',s2)
    print('PAi s',s3)
    print('PAi p',s4)
    print('fl s',s5)
    print('fl p',s6)
    print('mt s',s7)
    print('mt p',s8)

    mt_s_a_filt = [[] for par in range(len(participants))]
    mt_s_m_filt = [[] for par in range(len(participants))]
    for p in range(len(mt_s_a)):
        for val in range(len(mt_s_a[p])):
            if p==9:
                if mt_s_a[p][val]>mt_s_m[p][val]+2.37 or mt_s_a[p][val]<mt_s_m[p][val]-2.08:
                    print('aberrante valeur: participant ', p, ' , place val ', val)
                else:
                    mt_s_a_filt[p].append(mt_s_a[p][val])
                    mt_s_m_filt[p].append(mt_s_m[p][val])
            else:
                mt_s_a_filt[p].append(mt_s_a[p][val])
                mt_s_m_filt[p].append(mt_s_m[p][val])

    print('apres avoir enleve les valeurs out of LoA: ')
    t7,p7=ttest_rel([item for sublist in mt_s_m_filt for item in sublist],[item for sublist in mt_s_a_filt for item in sublist],axis=None)
    print('mt s',p7)

max_stars_repo_name: Lisa-pa/SAMAE

max_issues_repo_name: Lisa-pa/SAMAE

max_forks_repo_name: Lisa-pa/AponeurosesDetection

reason: import numpy,from scipy


text:
# coding:utf-8
import math
import numpy as np
import pandas as pd

if __name__ == "__main__":
    data_set = pd.read_csv("sample_data.csv")
    print("================")
    print(data_set)
    print("================")
    math_list = np.array(data_set.get("数学"))
    literature = np.array(data_set.get("国語"))

    print("平均値,", np.average(math_list))
    print("中央値,", np.median(math_list))

    count = np.bincount(math_list.astype("int32"))
    print("最頻値,", np.argmax(count))
    print("================")
    print("分散,", np.var(math_list))
    print("標準偏差", math.floor(np.std(math_list)))


text:
# Can we create a standalone executable?
# # target = "llvm -link-params"
# https://discuss.tvm.apache.org/t/can-we-create-a-standalone-executable/8773

import numpy as np
import tvm
import tvm.testing
import tvm.topi.testing
from tvm import relay
from tvm.contrib import graph_executor

import onnx
from onnx import TensorProto, helper, mapping, numpy_helper
import pdb

def get_input_data_shape_dict(graph_def, input_data):
    if isinstance(input_data, list):
        input_names = {}
        shape_dict = {}
        for i, _ in enumerate(input_data):
            input_names[i] = graph_def.graph.input[i].name
            shape_dict[input_names[i]] = input_data[i].shape
    else:
        input_names = graph_def.graph.input[0].name
        shape_dict = {input_names: input_data.shape}

    return input_names, shape_dict

def get_tvm_output_with_vm(
    graph_def, input_data, target, device, opset=None, freeze_params=False, convert_to_static=False
):
    """Generic function to execute and get tvm output with vm executor"""
    if not isinstance(input_data, list):
        input_data = [input_data]
    _, shape_dict = get_input_data_shape_dict(graph_def, input_data)

    mod, params = relay.frontend.from_onnx(
        graph_def, shape_dict, opset=opset, freeze_params=freeze_params
    )

    if convert_to_static:
        mod = relay.transform.DynamicToStatic()(mod)

    ex = relay.create_executor("vm", mod=mod, device=device, target=target)
    result = ex.evaluate()(*input_data, **params)
    if isinstance(result, tvm.runtime.NDArray):
        return result.numpy()
    return [r.numpy() for r in result]

def verify_simple_dynamic_model(a_shape, b_shape, target, dev):
    def verify_model(ex, a_shape, b_shape):
        a_array = np.random.uniform(size=a_shape).astype("float32")
        b_array = np.random.uniform(size=b_shape).astype("float32")
        # matmul
        out_np = np.matmul(a_array, b_array)
        # relu
        out_np[out_np < 0] = 0

        tvm_out = ex.evaluate()(a_array, b_array).numpy()
        tvm.testing.assert_allclose(out_np, tvm_out, rtol=1e-5, atol=1e-5)

    mul_node = helper.make_node("MatMul", ["a", "b"], ["out"])
    relu_node = helper.make_node("Relu", ["out"], ["relu"])

    a_array = np.random.uniform(size=a_shape).astype("float32")
    b_array = np.random.uniform(size=b_shape).astype("float32")
    # matmul
    out_np = np.matmul(a_array, b_array)

    graph = helper.make_graph(
        [mul_node, relu_node],
        "matmul_test",
        inputs=[
            helper.make_tensor_value_info("a", TensorProto.FLOAT, list(a_shape)),
            helper.make_tensor_value_info("b", TensorProto.FLOAT, list(b_shape)),
        ],
        outputs=[helper.make_tensor_value_info("relu", TensorProto.FLOAT, list(out_np.shape))],
    )

    model = helper.make_model(graph, producer_name="matmul_test")

    a_anys = [relay.Any()] * len(a_shape)
    b_anys = [relay.Any()] * len(b_shape)

    mod, params = relay.frontend.from_onnx(model, {"a": a_anys, "b": b_anys})
    # https://discuss.tvm.apache.org/t/relay-frontend-can-relay-take-none-include-shape/5772/2

    # # xxxx8888
    # # opt_level = 3
    # # with tvm.transform.PassContext(opt_level=opt_level):
    # executable = tvm.relay.backend.vm.compile(mod, target, params=params)
    # # code, lib = executable.save()
    # # Examples
    # # --------------------------------------------
    # #     import numpy as np
    # #     import tvm
    # #     from tvm import te
    # #     from tvm import relay
    # #     # define a simple network.
    # #     x = relay.var('x', shape=(10, 10))
    # #     f = relay.Function([x], x + x)
    # #     mod = tvm.IRModule({"main": f})
    # #     # create a Relay VM.
    # #     dev = tvm.cpu()
    # #     target = "llvm"
    # #     executable = relay.vm.compile(mod, target)
    # #     code, lib = executable.save()

    # #     # save and load the code and lib file.
    # #     tmp = tvm.contrib.utils.tempdir()
    # #     path_lib = tmp.relpath("lib.so")
    # #     lib.export_library(path_lib)

    # #     with open(tmp.relpath("code.ro"), "wb") as fo:
    # #         fo.write(code)

    # #     loaded_lib = tvm.runtime.load_module(path_lib)
    # #     loaded_code = bytearray(open(tmp.relpath("code.ro"), "rb").read())

    # #     # deserialize.
    # #     des_exec = tvm.runtime.vm.Executable.load_exec(loaded_code, loaded_lib)
    # #     # execute the deserialized executable.

    # #     x_data = np.random.rand(10, 10).astype('float32')
    # #     des_vm = tvm.runtime.vm.VirtualMachine(des_exec, dev)
    # #     res = des_vm.run(x_data)
    # #     print(res.numpy())
    # # pdb.set_trace()

    # ex = relay.create_executor("vm", mod=mod, device=dev, target=target)
    verify_model(ex, a_shape, b_shape)
    verify_model(ex, [a * 2 for a in a_shape], [b * 2 for b in b_shape])
    verify_model(ex, [a * 3 for a in a_shape], [b * 3 for b in b_shape])

    pdb.set_trace()

# TODO(mbrookhart, electriclilies): Add CUDA as a target once batch matmul is fixed
@tvm.testing.parametrize_targets("llvm")
def test_batch_matmul_dynamic_model(target, dev):
    verify_simple_dynamic_model((2, 3, 4, 3), (2, 3, 3, 4), target, dev)
    verify_simple_dynamic_model((2, 4, 3), (3, 4), target, dev)
    verify_simple_dynamic_model((2, 3, 4, 3), (3, 4), target, dev)

if __name__ == "__main__":
    target = "llvm"
    device = tvm.cpu(0)
    test_batch_matmul_dynamic_model(target, device)

max_stars_repo_name: delldu/Study

max_issues_repo_name: delldu/Study

max_forks_repo_name: delldu/Study


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

# Copyright 2017, IBM.
#
# This source code is licensed under the Apache License, Version 2.0 found in
# the LICENSE.txt file in the root directory of this source tree.

# pylint: disable=invalid-name,missing-docstring,broad-except

"""IBMQJob Test."""

import unittest
import time
from concurrent import futures
import numpy
from scipy.stats import chi2_contingency

from qiskit import (ClassicalRegister, QuantumCircuit, QuantumRegister,
                    QuantumJob)
import qiskit._compiler
from qiskit.backends.ibmq import IBMQProvider
from qiskit.backends.ibmq.ibmqjob import IBMQJob, IBMQJobError
from qiskit.backends.basejob import JobStatus
from .common import requires_qe_access, QiskitTestCase, slow_test

def lowest_pending_jobs(backends):
    """Returns the backend with lowest pending jobs."""
    backends = filter(lambda x: x.status.get('available', False), backends)
    by_pending_jobs = sorted(backends,
                             key=lambda x: x.status['pending_jobs'])
    return by_pending_jobs[0]

class TestIBMQJob(QiskitTestCase):
    """
    Test ibmqjob module.
    """

    @classmethod
    @requires_qe_access
    def setUpClass(cls, QE_TOKEN, QE_URL, hub=None, group=None, project=None):
        # pylint: disable=arguments-differ
        super().setUpClass()
        # create QuantumCircuit
        qr = QuantumRegister(2, 'q')
        cr = ClassicalRegister(2, 'c')
        qc = QuantumCircuit(qr, cr)
        qc.h(qr[0])
        qc.cx(qr[0], qr[1])
        qc.measure(qr, cr)
        cls._qc = qc
        cls._provider = IBMQProvider(QE_TOKEN, QE_URL, hub, group, project)
        cls._using_hub = bool(hub and group and project)

    def test_run_simulator(self):
        backend = self._provider.get_backend('ibmq_qasm_simulator')
        qobj = qiskit._compiler.compile(self._qc, backend)
        shots = qobj['config']['shots']
        quantum_job = QuantumJob(qobj, backend, preformatted=True)
        job = backend.run(quantum_job)
        result = job.result()
        counts_qx = result.get_counts(result.get_names()[0])
        counts_ex = {'00': shots/2, '11': shots/2}
        states = counts_qx.keys() | counts_ex.keys()
        # contingency table
        ctable = numpy.array([[counts_qx.get(key, 0) for key in states],
                              [counts_ex.get(key, 0) for key in states]])
        self.log.info('states: %s', str(states))
        self.log.info('ctable: %s', str(ctable))
        contingency = chi2_contingency(ctable)
        self.log.info('chi2_contingency: %s', str(contingency))
        self.assertGreater(contingency[1], 0.01)

    @slow_test
    def test_run_device(self):
        backends = self._provider.available_backends({'simulator': False})
        self.log.info('devices: %s', [b.name for b in backends])
        backend = lowest_pending_jobs(backends)
        self.log.info('using backend: %s', backend.name)
        qobj = qiskit._compiler.compile(self._qc, backend)
        shots = qobj['config']['shots']
        quantum_job = QuantumJob(qobj, backend, preformatted=True)
        job = backend.run(quantum_job)
        while not (job.done or job.exception):
            self.log.info(job.status)
            time.sleep(4)
        if job.exception:
            raise job.exception
        self.log.info(job.status)
        result = job.result()
        counts_qx = result.get_counts(result.get_names()[0])
        counts_ex = {'00': shots/2, '11': shots/2}
        states = counts_qx.keys() | counts_ex.keys()
        # contingency table
        ctable = numpy.array([[counts_qx.get(key, 0) for key in states],
                              [counts_ex.get(key, 0) for key in states]])
        self.log.info('states: %s', str(states))
        self.log.info('ctable: %s', str(ctable))
        contingency = chi2_contingency(ctable)
        self.log.info('chi2_contingency: %s', str(contingency))
        self.assertDictAlmostEqual(counts_qx, counts_ex, shots*0.1)

    @slow_test
    def test_run_async_simulator(self):
        IBMQJob._executor = futures.ThreadPoolExecutor(max_workers=2)
        backend = self._provider.get_backend('ibmq_qasm_simulator')
        self.log.info('submitting to backend %s', backend.name)
        num_qubits = 16
        qr = QuantumRegister(num_qubits, 'qr')
        cr = ClassicalRegister(num_qubits, 'cr')
        qc = QuantumCircuit(qr, cr)
        for i in range(num_qubits-1):
            qc.cx(qr[i], qr[i+1])
        qc.measure(qr, cr)
        qobj = qiskit._compiler.compile([qc]*10, backend)
        quantum_job = QuantumJob(qobj, backend, preformatted=True)
        num_jobs = 5
        job_array = [backend.run(quantum_job) for _ in range(num_jobs)]
        found_async_jobs = False
        timeout = 30
        start_time = time.time()
        while not found_async_jobs:
            check = sum([job.running for job in job_array])
            if check >= 2:
                self.log.info('found %d simultaneous jobs', check)
                break
            if all([job.done for job in job_array]):
                # done too soon? don't generate error
                self.log.warning('all jobs completed before simultaneous jobs '
                                 'could be detected')
                break
            for job in job_array:
                self.log.info('%s %s %s %s', job.status['status'], job.running,
                              check, job.job_id)
            self.log.info('-'*20 + ' ' + str(time.time()-start_time))
            if time.time() - start_time > timeout:
                raise TimeoutError('failed to see multiple running jobs after '
                                   '{0} s'.format(timeout))
            time.sleep(0.2)

        result_array = [job.result() for job in job_array]
        self.log.info('got back all job results')
        # Ensure all jobs have finished.
        self.assertTrue(all([job.done for job in job_array]))
        self.assertTrue(all([result.get_status() == 'COMPLETED' for result in result_array]))

        # Ensure job ids are unique.
        job_ids = [job.job_id for job in job_array]
        self.assertEqual(sorted(job_ids), sorted(list(set(job_ids))))

    @slow_test
    def test_run_async_device(self):
        backends = self._provider.available_backends({'simulator': False})
        backend = lowest_pending_jobs(backends)
        self.log.info('submitting to backend %s', backend.name)
        num_qubits = 5
        qr = QuantumRegister(num_qubits, 'qr')
        cr = ClassicalRegister(num_qubits, 'cr')
        qc = QuantumCircuit(qr, cr)
        for i in range(num_qubits-1):
            qc.cx(qr[i], qr[i+1])
        qc.measure(qr, cr)
        qobj = qiskit._compiler.compile(qc, backend)
        quantum_job = QuantumJob(qobj, backend, preformatted=True)
        num_jobs = 3
        job_array = [backend.run(quantum_job) for _ in range(num_jobs)]
        time.sleep(3)  # give time for jobs to start (better way?)
        job_status = [job.status['status'] for job in job_array]
        num_init = sum([status == JobStatus.INITIALIZING for status in job_status])
        num_queued = sum([status == JobStatus.QUEUED for status in job_status])
        num_running = sum([status == JobStatus.RUNNING for status in job_status])
        num_done = sum([status == JobStatus.DONE for status in job_status])
        num_error = sum([status == JobStatus.ERROR for status in job_status])
        self.log.info('number of currently initializing jobs: %d/%d',
                      num_init, num_jobs)
        self.log.info('number of currently queued jobs: %d/%d',
                      num_queued, num_jobs)
        self.log.info('number of currently running jobs: %d/%d',
                      num_running, num_jobs)
        self.log.info('number of currently done jobs: %d/%d',
                      num_done, num_jobs)
        self.log.info('number of errored jobs: %d/%d',
                      num_error, num_jobs)
        self.assertTrue(num_jobs - num_error - num_done > 0)

        # Wait for all the results.
        result_array = [job.result() for job in job_array]

        # Ensure all jobs have finished.
        self.assertTrue(all([job.done for job in job_array]))
        self.assertTrue(all([result.get_status() == 'COMPLETED' for result in result_array]))

        # Ensure job ids are unique.
        job_ids = [job.job_id for job in job_array]
        self.assertEqual(sorted(job_ids), sorted(list(set(job_ids))))

    def test_cancel(self):
        if not self._using_hub:
            self.skipTest('job cancellation currently only available on hubs')
        backends = self._provider.available_backends({'simulator': False})
        self.log.info('devices: %s', [b.name for b in backends])
        backend = backends[0]
        self.log.info('using backend: %s', backend.name)
        num_qubits = 5
        qr = QuantumRegister(num_qubits, 'qr')
        cr = ClassicalRegister(num_qubits, 'cr')
        qc = QuantumCircuit(qr, cr)
        for i in range(num_qubits-1):
            qc.cx(qr[i], qr[i+1])
        qc.measure(qr, cr)
        qobj = qiskit._compiler.compile(qc, backend)
        quantum_job = QuantumJob(qobj, backend, preformatted=True)
        num_jobs = 3
        job_array = [backend.run(quantum_job) for _ in range(num_jobs)]
        success = False
        self.log.info('jobs submitted: %s', num_jobs)
        while any([job.status['status'] == JobStatus.INITIALIZING for job in job_array]):
            self.log.info('jobs initializing')
            time.sleep(1)
        for job in job_array:
            job.cancel()
        while not success:
            job_status = [job.status for job in job_array]
            for status in job_status:
                self.log.info(status)
            if any([status['status'] == JobStatus.CANCELLED for status in job_status]):
                success = True
            if all([status['status'] == JobStatus.DONE for status in job_status]):
                raise IBMQJobError('all jobs completed before any could be cancelled')
            self.log.info('-' * 20)
            time.sleep(2)
        self.assertTrue(success)

    def test_job_id(self):
        backend = self._provider.get_backend('ibmq_qasm_simulator')
        qobj = qiskit._compiler.compile(self._qc, backend)
        quantum_job = QuantumJob(qobj, backend, preformatted=True)
        job = backend.run(quantum_job)
        while job.status['status'] == JobStatus.INITIALIZING:
            time.sleep(0.1)
        self.log.info('job_id: %s', job.job_id)
        self.assertTrue(job.job_id is not None)

    def test_get_backend_name(self):
        backend_name = 'ibmq_qasm_simulator'
        backend = self._provider.get_backend(backend_name)
        qobj = qiskit._compiler.compile(self._qc, backend)
        quantum_job = QuantumJob(qobj, backend, preformatted=True)
        job = backend.run(quantum_job)
        self.assertTrue(job.backend_name == backend_name)

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

max_stars_repo_name: Hosseinyeganeh/qiskit-core

max_issues_repo_name: Hosseinyeganeh/qiskit-core

max_forks_repo_name: Hosseinyeganeh/qiskit-core