text:
# -*- coding: utf-8 -*-
"""ReNS experiments - CIFAR10

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1byZ4xTfCK2x1Rhkxpl-Vv4sqA-bo4bis

# SETUP
"""

#@title Insatlling Pyorch

# !pip install torch
# !pip install torchvision

#@title Import Dependencies
import numpy as np
import torch
import torch.nn as nn
import torchvision.datasets as dsets
import torchvision.transforms as transforms
from torch.autograd import Variable
from tqdm import tqdm
from typing import Optional, Union, Tuple, List, Sequence, Iterable
import math
from scipy.spatial.distance import euclidean
from torch.nn.modules.utils import _pair
from torchvision import models
from sklearn.metrics import jaccard_score
import matplotlib.pyplot as plt
from models.models import RegularAutoEncoder, ModulatedAutoEncoder, PseudoRecAutoEncoder

"""# TRAINING"""

batch_size = 32
num_epochs = 5

transform = transforms.Compose(
    [transforms.ToTensor(),
     transforms.Normalize((0.5, 0.5, 0.5), (0.5, 0.5, 0.5))])

# Load MNIST data.
train_data = dsets.CIFAR10(root = './data', train = True,
                        transform = transform, download = True)

test_data = dsets.CIFAR10(root = './data', train = False,
                       transform = transform)

train_gen = torch.utils.data.DataLoader(dataset = train_data,
                                             batch_size = batch_size,
                                             shuffle = True)

test_gen = torch.utils.data.DataLoader(dataset = test_data,
                                      batch_size = batch_size,
                                      shuffle = False)

reflexor_size = 500
image_size = 32
channels = 3

# net = recurrentLayer(784, 784, 10, 5, 10, 0)
net1 = RegularAutoEncoder(channels * image_size ** 2, channels * image_size ** 2, reflexor_size)
net2 = ModulatedAutoEncoder(channels * image_size ** 2, channels * image_size ** 2, reflexor_size)
net3 = PseudoRecAutoEncoder(channels * image_size ** 2, channels * image_size ** 2, reflexor_size)

lr = .0001 # size of step
loss_function = nn.MSELoss()

# Unnormalize the image to display it
def img_fix(img):
  return np.transpose((img / 2 + 0.5).numpy(), (1, 2, 0))

# Commented out IPython magic to ensure Python compatibility.
train_losses = [[],[],[]]
test_losses = [[],[],[]]

real_imgs = [[],[],[]]
reconstructed_imgs = [[],[],[]]

param_counts = np.ones(3)

steps = [[],[],[]]

for num, net in enumerate([net1, net2, net3]):

  optimizer = torch.optim.Adam( net.parameters(), lr=lr)
  param_counts[num] = (sum(p.numel() for p in net.parameters() if p.requires_grad))

  for epoch in range(num_epochs):
    for i ,(images,labels) in enumerate(train_gen):
      #images = Variable(images.view(-1,28*28))
      labels = Variable(images.view(-1,3 * image_size ** 2))

      optimizer.zero_grad()
      outputs = net(images)
      loss = loss_function(outputs, labels)
      loss.backward()
      optimizer.step()

      if (i+1) % 300 == 0:
        temp_loss = loss.item()
        print('Epoch [%d/%d], Step [%d/%d], Loss: %.4f'
                   %(epoch+1, num_epochs, i+1, len(train_data)//batch_size, temp_loss))
        dupe = Variable(outputs[0].data, requires_grad=False)
        # plt.imshow(img_fix(images[0]))
        # plt.show()
        # plt.imshow(img_fix(dupe.view(3, image_size, image_size)))
        # plt.show()
        train_losses[num].append(temp_loss)
        steps[num].append((50000 * epoch) + ((i + 1) * batch_size))

        real_imgs[num].append(img_fix(images[0]))
        reconstructed_imgs[num].append(img_fix(dupe.view(3, image_size, image_size)))

        # Test Data

        score = 0
        total = 0
        for images,labels in test_gen:
          #images = Variable(images.view(-1,784))

          output = net(images)
          score += loss_function(output, images.view(-1, 3 * image_size ** 2)).item()
        test_losses[num].append((score))

plt.plot(steps[0], train_losses[0], label= "Baseline")
plt.plot(steps[1], train_losses[1], label= "Modulated")
plt.plot(steps[2], train_losses[2], label= "Recurrent with Modulation")
plt.xlabel('Iteration')
plt.ylabel('Loss')
plt.title('Training loss history')
plt.legend()
plt.show()

plt.plot(steps[0], test_losses[0], label= "Baseline")
plt.plot(steps[1], test_losses[1], label= "Modulated")
plt.plot(steps[2], test_losses[2], label= "Recurrent with Modulation")
plt.xlabel('Iteration')
plt.ylabel('Loss')
plt.title('Testing loss history')
plt.legend()
plt.show()

for num,count in enumerate(param_counts):
  param_counts[num] /= 1000

plt.bar(["Base", "Modulated", "ReNS"], param_counts)
plt.xlabel('Model')
plt.ylabel('# of thousands of Parameters')
plt.show()

from mpl_toolkits.axes_grid1 import ImageGrid

num_smaples = len(real_imgs[0])

for num in [0,1,2]:
  fig = plt.figure(figsize=(20.,20.))
  grid = ImageGrid(fig, 111,  # similar to subplot(111)
                  nrows_ncols=(2, num_smaples),  # creates 2x2 grid of axes
                  axes_pad=0.1,  # pad between axes in inch.
                  )

  for ax, im in zip(grid, real_imgs[num]+reconstructed_imgs[num]):
      # Iterating over the grid returns the Axes.
      ax.imshow(im)
      ax.axis("off")

  plt.show()

max_stars_repo_name: coopersigrist/RecurrentNeuralSystem-

max_issues_repo_name: coopersigrist/RecurrentNeuralSystem-

max_forks_repo_name: coopersigrist/RecurrentNeuralSystem-

reason: import numpy,from scipy


text:
import os
import sys
sys.path.insert(1, os.path.join(sys.path[0], '../'))
import numpy as np
import h5py
import argparse
import time
import logging
from sklearn import metrics
from utils import utilities, data_generator
import core
import inference_rt
import transfer_learning

# import keras
# from keras.models import Model
# from keras.layers import (Input, Dense, BatchNormalization, Dropout, Lambda,
#                           Activation, Concatenate, LSTM, GRU, Reshape)
# import keras.backend as K
# from keras.optimizers import Adam

try:
    import cPickle
except BaseException:
    import _pickle as cPickle

import vim_params as vp

def train(args):
    # model_type = args.model_type
    # time_steps = 10
    # freq_bins = 128
    # classes_num = 527

    # # Hyper parameters
    # hidden_units = 1024
    # drop_rate = 0.5

    args.batch_size = vp.BATCH_SIZE    #128, 192

    # Train
    core.train(args)

# Main
if __name__ == '__main__':

    parser = argparse.ArgumentParser(description="")

    parser.add_argument('--data_dir', type=str, required=True)

    parser.add_argument('--workspace', type=str, required=True)

    parser.add_argument('--mini_data', action='store_true',
                        default=False)

    parser.add_argument('--quantize', action='store_true',
                        default=False)  # quantize the model

    parser.add_argument('--balance_type', type=str,
                        default='balance_in_batch',
                        choices=['no_balance', 'balance_in_batch'])

    parser.add_argument('--model_type', type=str, required=True,
                        choices=['decision_level_max_pooling',
                                 'decision_level_average_pooling',
                                 'decision_level_single_attention',
                                 'decision_level_multi_attention',
                                 'feature_level_attention'])

    parser.add_argument('--learning_rate', type=float, default=1e-3)

    subparsers = parser.add_subparsers(dest='mode')
    parser_train = subparsers.add_parser('train')
    parser_get_avg_stats = subparsers.add_parser('get_avg_stats')
    parser_inference = subparsers.add_parser('inference_rt')
    parser_train = subparsers.add_parser('transfer')

    args = parser.parse_args()

    args.filename = utilities.get_filename(__file__)

    # Logs
    logs_dir = os.path.join(args.workspace, 'logs', args.filename)
    utilities.create_folder(logs_dir)
    logging = utilities.create_logging(logs_dir, filemode='w')

    logging.info(os.path.abspath(__file__))
    logging.info(args)

    if args.mode == "train":
        train(args)

    elif args.mode == 'transfer':
        args.batch_size = vp.BATCH_SIZE
        transfer_learning.train(args)

    elif args.mode == 'get_avg_stats':
        args.bgn_iteration = 10000
        args.fin_iteration = 50001
        args.interval_iteration = 5000
        utilities.get_avg_stats(args)

    elif args.mode == 'inference_rt':
        inference_rt.core(args)

    else:
        raise Exception("Error!")

ext: py

reason: import numpy


text:
"""
Assignment 1 Question 2 Part C script

@author         Ian Tongs, 27765369
@since          21 April 2020
"""

# Imports:

import matplotlib.pyplot as plt
import numpy as np
from matplotlib import rc

###### Functions:
# Cobwebbing:

# Set in function:
func = lambda x, m, c: m*x+c

# Actual Plotter:
def plot_cobweb(f, l, b, x0, xmin, xmax, title, nmax=40):
    # Set up plot:
    x = np.linspace(xmin, xmax, 1000)
    fig = plt.figure()
    ax = fig.add_subplot(111)

    # Plot y = f(x) and y = x as required for the web
    ax.plot(x, x, c='#444444', lw=2, label='y=x')
    ax.plot(x, f(x, l, b), c='#FF0000', lw=2, label='y=f(x)')

    # Create Numpy arrays for x and y (the webbing)
    px, py = np.empty((2, nmax+1, 2))

    # Iterate x = f(x) for each allowed step, starting at (x0, 0), for the webbing itself
    px[0], py[0] = x0, 0
    for n in range(1, nmax, 2):
        px[n] = px[n-1]
        py[n] = f(px[n-1], l, b)
        px[n+1] = py[n]
        py[n+1] = py[n]

    # Plot the path traced out by the iteration.
    ax.plot(px, py, c='b', alpha=0.7)

    # Add dashes and axis labels
    ax.minorticks_on()
    ax.grid(which='minor', alpha=0.5)
    ax.grid(which='major', alpha=0.5)
    ax.set_aspect('equal')
    ax.set_xlabel('$x_t$')
    ax.set_ylabel('$x_{t+1}$')
    ax.set_title(title)

    # Show the graph and add a legend
    plt.legend()
    plt.show()

# |m| < 1, with c = -2
# plot_cobweb(func, 0.5, -2, 35, -10, 40, 'Absolute value of m less than 1:', 20)

# |m| > 1, with c = 2
# plot_cobweb(func, 2, 2, 1, -10, 40, 'Absolute value of m greater than 1:', 8)

# m = 1, with c = 0
# plot_cobweb(func, 1, 0, 15, -10, 40, 'm = 1 and c = 0', 8)

# m = 1, with c = 3
# plot_cobweb(func, 1, 3, 15, -10, 40, 'm = 1 and c = 3', 16)

# m = -1, with c = 0
# plot_cobweb(func, -1, 0, 6, -10, 40, 'm = -1', 8)


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Author:
    Lennart van Sluijs
Date:
    Dec 21 2018
Description:
    Library of usefull auxillary functions and variables
    used by multiple scripts.
"""
#%%
import os
from astropy.io import fits
import numpy as np
from matplotlib import pyplot as plt
from scipy.ndimage import affine_transform

#%%
def bin2d(x, binsize):
    """
    Description:
        Bin 2D data to specified binsize.
    Input:
        x - 2D numpy array
        binsize - integer specifying binsize
    Output:
        xb - binned list or numpy array
    """

    # bin the rows of 2D array x by the binsize
    xb = np.mean(x.reshape(int(x.shape[0]/binsize), binsize, x.shape[1]),
               axis = 1)

    return xb

def bin1d(x, binsize):
    """
    Description:
        Bin 1D data to specified binsize.
    Input:
        x - 1D numpy array
        binsize - integer specifying binsize
    Output:
        xb - binned list or numpy array
    """

    # bin the rows of 1D array x by the binsize
    xb = np.mean(x.reshape(int(len(x)/binsize), binsize), axis = 1)

    return xb

def get_flist(dirname):
    """
    Description:
        Get a sorted list of all files in the directory
    Input:
        dirname - name of the input directory
    """

    # get file list and sort
    flist = [fname for fname in os.listdir(dirname) if \
             os.path.isfile(os.path.join(dirname, fname))]
    flist = sorted(flist)

    return flist

def findall_containing_string(list, string, inelement = False):
    """
    Description:
        Find all elements containing a substring and return those elements.
    Input:
        list - list object
        string - string to look for
        invert - if True, find all that do not contain string
    """

    # return updated list (not) containing the substring
    if inelement:
        return [element for element in list if string in element]
    else:
        return [element for element in list if string not in element]

def wfits(fname, im):
    """
    Description:
        Write image data to a .fits file.
    Input:
        fname - file name
        im - image data
    """

    # wfits - write im to file fname, automatically overwriting any old file
    hea = fits.PrimaryHDU(im)
    hea.writeto(fname, overwrite=True)

def cen_rot2(im, rot, dim_out, offset1=(0,0), offset2=(0,0), order=2):
    """
    Description:
        takes a cube of images im, and a set of rotation angles in rot,
        and translates the middle of the frame with a size dim_out to the middle of
        a new output frame with an additional rotation of rot.
    Input:
        im - image data
        rot - rotation angle
        dim_out - image output dimension
        offset1 - first offset
        offset2 - second offset
        order - used order
    dst - rotated image
        """

    # convert to radian
    a = rot * np.pi / 180.

    # make a rotation matrix
    transform=np.array([[np.cos(a),-np.sin(a)],[np.sin(a),np.cos(a)]])

    # -0.5 is there for fencepost counting error
    c_in = np.array(offset1) - 0.5
    c_out = 0.5 * np.array(dim_out) - 0.5

    # c_out has to be pre-rotated to make offset correct
    offset = c_in - c_out.dot(transform) - np.array(offset2).dot(transform)

    # perform the transformation
    dst=affine_transform( \
        im,transform.T, order=order,offset=offset, \
        output_shape=dim_out, cval=np.nan)

    return(dst)

def discrete_cmap(N, base_cmap=None):
    """
    Description:
        Create an N-bin discrete colormap from the specified input map.
    Input:
        N - integer, number of colors to use
        base_cmap - base colormap to use
    """

    # get base of discrete cmap
    base = plt.cm.get_cmap(base_cmap, N)

    return base

#%%
"""
Initizalize tableau20 color scheme used by many of the scripts.
"""

# define color scheme
tableau20 = [(31, 119, 180), (174, 199, 232), (255, 127, 14), (255, 187, 120),
             (44, 160, 44), (152, 223, 138), (214, 39, 40), (255, 152, 150),
             (148, 103, 189), (197, 176, 213), (140, 86, 75), (196, 156, 148),
             (227, 119, 194), (247, 182, 210), (127, 127, 127),
             (199, 199, 199), (188, 189, 34), (219, 219, 141), (23, 190, 207),
             (158, 218, 229)]

# scale the RGB values to the [0, 1] range
for i in range(len(tableau20)):
    r, g, b = tableau20[i]
    tableau20[i] = (r / 255., g / 255., b / 255.)

reason: import numpy,from scipy,from astropy


text:
import time
import numpy as np
from scipy.constants import (physical_constants as phys, c, m_u, hbar, pi, e, m_e)

from nbody import NBody

class Ion:
    name = 0
    def __init__(self, r=np.zeros(3), v=np.zeros(3), m=m_u, k=0, Z=1):
        self.r = r
        self.v = v
        self.m = m
        self.Z = Z
        self.k = k
        self.name = Ion.name
        Ion.name += 1
        return

    def __eq__(self, other):
        return self.name == other.name

    def __hash__(self):
        return hash(self.name)

    def __str__(self):
        return "{:06.2f} {:06.2f}".format(self.r[0], self.v[0])

class Trap(NBody):
    """
    linear quadrupole trap
    fill it with ions:
    positions
    velocities
    weights (weights of ions in the sum used to calculate field
             i.e. E=k*Z*e*r/r**3, so weight = k*Z*e)
    f2a (field to acceleration for each ion, i.e. Z*e/m)
    viscosities (a = -kv, k is viscosity)
    """
    micro = True # calculate micromotion
    V = 350 # V, AC
    w = 2*pi*18.0e6 # 2*pi*Hz
    U = 0 # V, DC
    Uz = 100 # V, DC
    xySize = 1.14e-3 # m, radius
    zSize = 7e-3 #m, half of length

    beamR = 1e-3 # m
    # use self.config(Z=1, m=171) to change any of the above
    # or during self.init(1e4, sigmaV=1e-2, Z=1, m=171)

    # load ions, or use 'ions' random ones if 'ions' is a number
    def __init__(self, ions, sigmaV=1, viscosity=1e3, **kwargs):
        self.config(**kwargs)
        # assume Thorium trap by default
        E = e/(4*pi*phys['electric constant'][0])
        Z = 1
        try:
            N = int(ions)
            rnd = True
        except TypeError:
            rnd = False
        if rnd:
            np.random.seed(42)
            positions = (2*np.random.rand(N,3)-1)*min(self.xySize, self.zSize)/100
            velocities = np.random.randn(N,3)*sigmaV
            weights = np.array([Z*E]*N)

            # TODO
            # this should really be of the shape (N,3)
            viscosities = np.array([viscosity]*N)

            f2a = np.array([Z*e/(25*m_u)]*N)
        else:
            positions = np.array(ion.r for ion in ions)
            velocities = np.array(ion.v for ion in ions)
            weights = np.array(E*ion.Z for ion in ions)
            viscosities = np.array(ion.k for ion in ions)
            f2a = np.array(e*ion.Z/ion.m for ion in ions)
        self.N = len(positions)
        return super().__init__(positions, velocities, weights,
                                viscosities, f2a)

    def config(self, **kwargs):
        # careful not to change something you are not supposed to
        # I won't fool-proof this
        for key,value in kwargs.items():
            setattr(self, key, value)
        self.u = self.U/(self.xySize*self.xySize)
        self.v = self.V/(self.xySize*self.xySize)
        self.uz = self.Uz/(self.zSize*self.zSize)
        return

    def confinement(self):
        if self.micro:
            xConf = self.u - self.v*np.cos(self.w*self.T)
            yConf = -xConf
            if self.T*self.w > 42*2*pi:
                self.T -= 42*2*pi/self.w
        else:
            xConf = self.v/1.414
            yConf = xConf
        zConf = self.uz/1.414
        return (xConf, yConf, zConf)

    def size(self):
        return self.xySize, self.xySize, self.zSize

    def viscosityR(self):
        return self.beamR

    def clean(self):
        # leave everyone who is out of range of our trap
        newPos = np.empty([self.N,3], dtype=np.float32)
        newVel = np.empty([self.N,3], dtype=np.float32)
        newWei = np.empty(self.N, dtype=np.float32)
        newVis = np.empty(self.N, dtype=np.float32)
        newf2a = np.empty(self.N, dtype=np.float32)
        j = 0
        for i in range(self.N):
            if (abs(self.positions[i][0]) < self.xySize and
                abs(self.positions[i][1]) < self.xySize and
                abs(self.positions[i][2]) < self.zSize):
                newPos[j] = self.positions[i]
                newVel[j] = self.velocities[i]
                newWei[j] = self.weights[i]
                newVis[j] = self.viscosities[i]
                newf2a[j] = self.f2a[i]
                j += 1
        if j < 1:
            print("WE'VE LOST THEM ALL", self.velocities[0])
            return
        self.positions = newPos[:j,:]
        self.velocities = newVel[:j,:]
        self.weights = newWei[:j]
        self.viscosities = newVis[:j]
        self.f2a = newf2a[:j]
        self.N = len(self.weights)

    def __call__(self, dt=1e-8, steps=1, debug=False):
        before = time.perf_counter()

        super().__call__(dt, steps)
        self.clean()

        return ('{0:10.2e}'.format(self.t*1000),
                '{0:05d}'.format(self.N),
                '{:10.2e}'.format(time.perf_counter() - before))


text:
#!/usr/bin/env python
# @Copyright 2007 Kristjan Haule
import optparse
import os, shutil
import subprocess
from numpy import zeros, log
import indmffile
from utils import DmftEnvironment, W2kEnvironment
import sys

def print_prompt(text):
    print "----->", text

inpt={}  # input from command line
if len(sys.argv)>1:
    if sys.argv[1] in ['-h', '--help']:
        print """Program : inid_dmft.py

Initializes the EDMFTF run by producing case.indmfl and case indmfi files.
The first connects DMFT problem with the Kohn-Sham orbitals,
      and the second connects DMFT with the impurity solver.

Usage init_dmft.py  [-option val]

 where option's are:
   -h, --help     -- displays this messgae
   -ca  1,2,3     -- list of correlated atoms (which atoms from
                                 case.struct are treated by DMFT)
   -ot  d,d,d     -- orbital type being correlated
                     (which orbital is treated by DMFT, i.e, d or f)
   -qs  7,7,7     -- q-split described the symmetry of the orbital.
                      (7 cubic, 2 real harmonics,..).
                      For help execute -h qs
   -p  5          -- projector type (For help execute -h p )
   -hr -10,10     -- hybridization window in eV
   -a  i          -- real (r) or imaginary (i) axis calculation
   -us "1,2 3,4"  -- unique set of orbitals, i.e., atom 1 and 2 map to
                     an equivalent imputy problem, and 3,4 map to the
                     second equivalent impurity problem.
   -cl "1,2 3,4"  -- grouping of atoms into cluster for cluster-DMFT
                     by default atoms are treated by single-site DMFT
"""
        if len(sys.argv)>2 and sys.argv[2] == 'qs':
            print indmffile.qsplit_doc
        if len(sys.argv)>2 and sys.argv[2] == 'p':
            print indmffile.projector_doc
        sys.exit(0)
    else:
        for i in range(len(sys.argv)/2):
            inpt[sys.argv[2*i+1][1:]] = sys.argv[2*i+2]
        print 'given options are', inpt

dmfenv = DmftEnvironment()
w2kenv = W2kEnvironment()

indmf = indmffile.Indmf(w2kenv.case)

# First, check if case.indmf exists
create_indmf = True
if indmf.file_exists():
    prompt = "File `%s` already exists.  Do you want to create new inputs? (y/n): " % indmf.filename()
    userin = raw_input(prompt).lower().strip()
    if userin.lower() == 'n':
        create_indmf = False

# create new case.indmf (if desired)
if create_indmf:
    indmf.user_input(inpt)
    indmf.write()
    fname = indmf.filename()
    fnametemp = fname + '_st'

    if inpt:
        finished_editing = True  # no editing when using command line
    else:
        finished_editing = False
    while not finished_editing:
        shutil.copy(fname, fnametemp)  # copy to temp file for user to edit
        print_prompt("Edit %s file." % (fnametemp,))
        subprocess.call(w2kenv.EDITOR.split() + [fnametemp])

        # read user-edited file in just to check for any syntax errors
        indmf_st = indmffile.Indmf(w2kenv.case)
        indmf_st.extn = 'indmf_st'
        try:
            indmf_st.read()
            finished_editing = True
            shutil.copy(fnametemp, fname)  # move it back
        except:
            print_prompt("Your edits have introduced syntax errors in %s." % (fnametemp,))
            prompt = "Do you want to edit the file again? (y/n): "
            userin = raw_input(prompt).lower().strip()
            if userin.lower() == 'n':
                finished_editing = True

# ask user if this is a SO run
so_run = False
if os.path.isfile(w2kenv.case+".inso") and os.path.getsize(w2kenv.case+".inso")>0 :
    print 'Found '+w2kenv.case+'.inso file, hence assuming so-coupling exists. Switching -so switch!'
    so_run = True

#prompt = "Is this a spin-orbit run? (y/n): "
#userin = raw_input(prompt).lower().strip()
#if userin.lower() == 'y':
#    so_run = True

# run sigen to calculate siginds, cftrans and write output to case.indmfl file
sigen = os.path.join( dmfenv.ROOT, 'sigen.py' )
#sigen = './sigen.py'
args = ["--so"] if so_run else []
ireturn = subprocess.call([sigen] + args, shell=True,stdout=sys.stdout,stderr=sys.stdout)
if ireturn:
    print "Error executing sigen.py."
    exit(1)

# copy case.indmfl --> case.indmfl_st to allow user to make changes
findmfl = w2kenv.case + '.indmfl'
findmfltemp = findmfl + '_st'

if inpt:
    finished_editing = True  # no editing when using command line
else:
    finished_editing = False
while not finished_editing:
    shutil.copy(findmfl, findmfltemp)
    print_prompt("Edit %s file." % (findmfltemp,))
    subprocess.call(w2kenv.EDITOR.split() + [findmfltemp])

    # read user-edited file in just to check for any syntax errors
    inl = indmffile.Indmfl(w2kenv.case)
    inl.extn = 'indmfl_st'
    try:
        inl.read()
        finished_editing = True
        shutil.copy(findmfltemp, findmfl)  # move it back
    except:
        print_prompt("Your edits have introduced syntax errors in %s." % (findmfltemp,))
        prompt = "Do you want to edit the file again? (y/n): "
        userin = raw_input(prompt).lower().strip()
        if userin.lower() == 'n':
            finished_editing = True

# Generate case.indmfi file -----------------------------------

indmf.read()
inl = indmffile.Indmfl(w2kenv.case)
inl.read()

f = open(w2kenv.case + '.indmfi', 'w')
print >> f, len(indmf.ucps), "  # number of sigind blocks"
for iucp,icps in indmf.ucps.iteritems():
    icp = icps[0]  # just take sigind of first correlated problem  (TODO: modify appropriately for broken symmetry runs)
    sigind = inl.siginds[icp]
    print >> f, sigind.shape[0], "  # dimension of this sigind block"
    max_sigfig = 1 + int(log(max(sigind.flat))/log(10))
    format = '%' + str(max_sigfig) + 'd'
    for row in sigind:
        print >> f, ' '.join([format % elem for elem in row])
f.close()

# Generate blank sig.inp --------------------------------------

#prompt = "Generate blank sig.inp? (y/n): "
#userin = raw_input(prompt).lower().strip()
#if userin.lower() == 'y':
#
#    cmd = os.path.join( dmfenv.ROOT, 'szero.py' )
#    stdoutdata, stderrdata = subprocess.Popen(cmd, stdout=subprocess.PIPE).communicate()

    #
    #
    ## use siginds to determine number of columns needed in sig.inp
    #inl = indmffile.Indmfl(w2kenv.case)
    #inl.read()
    #nsigs = max([max(sigind.flat) for sigind in inl.siginds.values()])  # assumes siginds start with 1
    #ncols = 2*nsigs   # real/imag for each sigma column
    #
    ## use gaumesh.py to generate standard mesh
    #gaumesh = os.path.join( dmfenv.ROOT, 'gaumesh.py' )
    ## gaumesh = './gaumesh.py'
    #args = [
    #    'x0=[0,0]',
    #    'dx0=[0.25,0.01]',
    #    'fwhm=[30, 3]',
    #    'xmin=-20',
    #    'xmax=20'
    #    ]
    #stdoutdata, stderrdata = subprocess.Popen([gaumesh] + args, stdout=subprocess.PIPE).communicate()
    #omegas = [x for x in stdoutdata.split('\n') if x.strip() and not x.strip().startswith('#')]
    #
    ## if sig.inp already exists, make backup
    #if os.path.isfile('sig.inp'):
    #    print "File `sig.inp` already exists, backing up to 'sig.inp.bak'"
    #    shutil.move('sig.inp', 'sig.inp.bak')
    #
    ## now generate zeros and write to file
    #f = open('sig.inp', 'w')
    #nsig_zeros = ', '.join(['0.0']*nsigs)       # first generate s_oo and Edc
    #print >> f, '# s_oo = [' + nsig_zeros + ']'
    #print >> f, '# Edc  = [' + nsig_zeros + ']'
    #ncol_zeros = ' '.join(['0.0']*ncols)       # now write out main file
    #for omega in omegas:
    #    print >> f, omega, ncol_zeros
    #f.close()

#    print 'Blank self-energy written to `sig.inp`.'

reason: from numpy


text:
# Copyright 2020 MIT Probabilistic Computing Project.
# See LICENSE.txt

from math import log

import pytest

import numpy

from sppl.distributions import choice
from sppl.distributions import gamma
from sppl.distributions import norm
from sppl.dnf import dnf_to_disjoint_union
from sppl.math_util import allclose
from sppl.math_util import isinf_neg
from sppl.math_util import logdiffexp
from sppl.math_util import lognorm
from sppl.math_util import logsumexp
from sppl.sets import Interval
from sppl.sets import inf as oo
from sppl.spe import LeafSPE
from sppl.spe import ProductSPE
from sppl.spe import SumSPE
from sppl.transforms import Exp
from sppl.transforms import Id
from sppl.transforms import Log

def test_product_distribution_normal_gamma_basic():
    X1 = Id('X1')
    X2 = Id('X2')
    X3 = Id('X3')
    X4 = Id('X4')
    children = [
        ProductSPE([
            X1 >> norm(loc=0, scale=1),
            X4 >> norm(loc=10, scale=1),
        ]),
        X2 >> gamma(loc=0, a=1),
        X3 >> norm(loc=2, scale=3)
    ]
    spe = ProductSPE(children)
    assert spe.children == (
        children[0].children[0],
        children[0].children[1],
        children[1],
        children[2],)
    assert spe.get_symbols() == frozenset([X1, X2, X3, X4])
    assert spe.size() == 5

    samples = spe.sample(2)
    assert len(samples) == 2
    for sample in samples:
        assert len(sample) == 4
        assert all([X in sample for X in (X1, X2, X3, X4)])

    samples = spe.sample_subset((X1, X2), 10)
    assert len(samples) == 10
    for sample in samples:
        assert len(sample) == 2
        assert X1 in sample
        assert X2 in sample

    samples = spe.sample_func(lambda X1, X2, X3: (X1, (X2**2, X3)), 1)
    assert len(samples) == 1
    assert len(samples[0]) == 2
    assert len(samples[0][1]) == 2

    with pytest.raises(ValueError):
        spe.sample_func(lambda X1, X5: X1 + X4, 1)

def test_product_inclusion_exclusion_basic():
    X = Id('X')
    Y = Id('Y')
    spe = ProductSPE([X >> norm(loc=0, scale=1), Y >> gamma(a=1)])

    a = spe.logprob(X > 0.1)
    b = spe.logprob(Y < 0.5)
    c = spe.logprob((X > 0.1) & (Y < 0.5))
    d = spe.logprob((X > 0.1) | (Y < 0.5))
    e = spe.logprob((X > 0.1) | ((Y < 0.5) & ~(X > 0.1)))
    f = spe.logprob(~(X > 0.1))
    g = spe.logprob((Y < 0.5) & ~(X > 0.1))

    assert allclose(a, spe.children[0].logprob(X > 0.1))
    assert allclose(b, spe.children[1].logprob(Y < 0.5))

    # Pr[A and B]  = Pr[A] * Pr[B]
    assert allclose(c, a + b)
     # Pr[A or B] = Pr[A] + Pr[B] - Pr[AB]
    assert allclose(d, logdiffexp(logsumexp([a, b]), c))
    # Pr[A or B] = Pr[A] + Pr[B & ~A]
    assert allclose(e, d)
    # Pr[A and B]  = Pr[A] * Pr[B]
    assert allclose(g, b + f)
    # Pr[A or (B & ~A)] = Pr[A] + Pr[B & ~A]
    assert allclose(e, logsumexp([a, b+f]))

    # (A => B) => Pr[A or B] = Pr[B]
    # i.e.,k (X > 1) => (X > 0).
    assert allclose(spe.logprob((X > 0) | (X > 1)), spe.logprob(X > 0))

    # Positive probability event.
    # Pr[A] = 1 - Pr[~A]
    event = ((0 < X) < 0.5) | ((Y < 0) & (1 < X))
    assert allclose(
        spe.logprob(event),
        logdiffexp(0, spe.logprob(~event)))

    # Probability zero event.
    event = ((0 < X) < 0.5) & ((Y < 0) | (1 < X))
    assert isinf_neg(spe.logprob(event))
    assert allclose(spe.logprob(~event), 0)

def test_product_condition_basic():
    X = Id('X')
    Y = Id('Y')
    spe = ProductSPE([X >> norm(loc=0, scale=1), Y >> gamma(a=1)])

    # Condition on (X > 0) and ((X > 0) | (Y < 0))
    # where the second clause reduces to first as Y < 0
    # has probability zero.
    for event in [(X > 0), (X  > 0) | (Y < 0)]:
        dX = spe.condition(event)
        assert isinstance(dX, ProductSPE)
        assert dX.children[0].symbol == Id('X')
        assert dX.children[0].conditioned
        assert dX.children[0].support == Interval.open(0, oo)
        assert dX.children[1].symbol == Id('Y')
        assert not dX.children[1].conditioned
        assert dX.children[1].Fl == 0
        assert dX.children[1].Fu == 1

    # Condition on (Y < 0.5)
    dY = spe.condition(Y < 0.5)
    assert isinstance(dY, ProductSPE)
    assert dY.children[0].symbol == Id('X')
    assert not dY.children[0].conditioned
    assert dY.children[1].symbol == Id('Y')
    assert dY.children[1].conditioned
    assert dY.children[1].support == Interval.Ropen(0, 0.5)

    # Condition on (X > 0) & (Y < 0.5)
    dXY_and = spe.condition((X > 0) & (Y < 0.5))
    assert isinstance(dXY_and, ProductSPE)
    assert dXY_and.children[0].symbol == Id('X')
    assert dXY_and.children[0].conditioned
    assert dXY_and.children[0].support == Interval.open(0, oo)
    assert dXY_and.children[1].symbol == Id('Y')
    assert dXY_and.children[1].conditioned
    assert dXY_and.children[1].support == Interval.Ropen(0, 0.5)

    # Condition on (X > 0) | (Y < 0.5)
    event = (X > 0) | (Y < 0.5)
    dXY_or = spe.condition((X > 0) | (Y < 0.5))
    assert isinstance(dXY_or, SumSPE)
    assert all(isinstance(d, ProductSPE) for d in dXY_or.children)
    assert allclose(dXY_or.logprob(X > 0),dXY_or.weights[0])
    samples = dXY_or.sample(100, prng=numpy.random.RandomState(1))
    assert all(event.evaluate(sample) for sample in samples)

    # Condition on a disjoint union with one term in second clause.
    dXY_disjoint_one = spe.condition((X > 0) & (Y < 0.5) | (X <= 0))
    assert isinstance(dXY_disjoint_one, SumSPE)
    component_0 = dXY_disjoint_one.children[0]
    assert component_0.children[0].symbol == Id('X')
    assert component_0.children[0].conditioned
    assert component_0.children[0].support == Interval.open(0, oo)
    assert component_0.children[1].symbol == Id('Y')
    assert component_0.children[1].conditioned
    assert component_0.children[1].support == Interval.Ropen(0, 0.5)
    component_1 = dXY_disjoint_one.children[1]
    assert component_1.children[0].symbol == Id('X')
    assert component_1.children[0].conditioned
    assert component_1.children[0].support == Interval(-oo, 0)
    assert component_1.children[1].symbol == Id('Y')
    assert not component_1.children[1].conditioned

    # Condition on a disjoint union with two terms in each clause
    dXY_disjoint_two = spe.condition((X > 0) & (Y < 0.5) | ((X <= 0) & ~(Y < 3)))
    assert isinstance(dXY_disjoint_two, SumSPE)
    component_0 = dXY_disjoint_two.children[0]
    assert component_0.children[0].symbol == Id('X')
    assert component_0.children[0].conditioned
    assert component_0.children[0].support == Interval.open(0, oo)
    assert component_0.children[1].symbol == Id('Y')
    assert component_0.children[1].conditioned
    assert component_0.children[1].support == Interval.Ropen(0, 0.5)
    component_1 = dXY_disjoint_two.children[1]
    assert component_1.children[0].symbol == Id('X')
    assert component_1.children[0].conditioned
    assert component_1.children[0].support == Interval(-oo, 0)
    assert component_1.children[1].symbol == Id('Y')
    assert component_1.children[1].conditioned
    assert component_1.children[1].support == Interval(3, oo)

    # Some various conditioning.
    spe.condition((X > 0) & (Y < 0.5) | ((X <= 1) | ~(Y < 3)))
    spe.condition((X > 0) & (Y < 0.5) | ((X <= 1) & (Y < 3)))

def test_product_condition_or_probabilithy_zero():
    X = Id('X')
    Y = Id('Y')
    spe = ProductSPE([X >> norm(loc=0, scale=1), Y >> gamma(a=1)])

    # Condition on event which has probability zero.
    event = (X > 2) & (X < 2)
    with pytest.raises(ValueError):
        spe.condition(event)
    assert spe.logprob(event) == -float('inf')

    # Condition on event which has probability zero.
    event = (Y < 0) | (Y < -1)
    with pytest.raises(ValueError):
        spe.condition(event)
    assert spe.logprob(event) == -float('inf')
    # Condition on an event where one clause has probability
    # zero, yielding a single product.
    spe_condition = spe.condition((Y < 0) | ((Log(X) >= 0) & (1 <= Y)))
    assert isinstance(spe_condition, ProductSPE)
    assert spe_condition.children[0].symbol == X
    assert spe_condition.children[0].conditioned
    assert spe_condition.children[0].support == Interval(1, oo)
    assert spe_condition.children[1].symbol == Y
    assert spe_condition.children[1].conditioned
    assert spe_condition.children[0].support == Interval(1, oo)

    # We have (X < 2) & ~(1 < exp(|3X**2|) is empty.
    # Thus Y remains unconditioned,
    #   and X is partitioned into (-oo, 0) U (0, oo) with equal weight.
    event = (Exp(abs(3*X**2)) > 1) | ((Log(Y) < 0.5) & (X < 2))
    spe_condition = spe.condition(event)
    #
    # The most concise representation of spe_condition is:
    #   (Product (Sum [.5 .5] X|X<0 X|X>0) Y)
    assert isinstance(spe_condition, ProductSPE)
    assert isinstance(spe_condition.children[0], SumSPE)
    assert spe_condition.children[0].weights == (-log(2), -log(2))
    assert spe_condition.children[0].children[0].conditioned
    assert spe_condition.children[0].children[1].conditioned
    assert spe_condition.children[0].children[0].support \
        in [Interval.Ropen(-oo, 0), Interval.Lopen(0, oo)]
    assert spe_condition.children[0].children[1].support \
        in [Interval.Ropen(-oo, 0), Interval.Lopen(0, oo)]
    assert spe_condition.children[0].children[0].support \
        != spe_condition.children[0].children[1].support
    assert spe_condition.children[1].symbol == Y
    assert not spe_condition.children[1].conditioned
    #
    # However, now we solve the constraint and factor the query into
    # a disjoint union at the root of the SPE, not once per node, we
    # no longer have guarantees of generating the smallest network.
    # In this case, the answer is the more verbose:
    #   (Sum [.5 .5] (Product X|X<0 Y) (Product X|X>0 Y)
    # assert isinstance(spe_condition, SumSPE)
    # assert spe_condition.weights == (-log(2), -log(2))
    # assert isinstance(spe_condition.children[0], ProductSPE)
    # assert X \
    #     == spe_condition.children[0].children[0].symbol \
    #     == spe_condition.children[1].children[0].symbol
    # assert spe_condition.children[0].children[0].conditioned
    # assert spe_condition.children[1].children[0].conditioned
    # assert spe_condition.children[0].children[0].support \
    #     == Interval.Ropen(-oo, 0)
    # assert spe_condition.children[1].children[0].support \
    #     == Interval.Lopen(0, oo)
    # assert Y \
    #     == spe_condition.children[0].children[1].symbol \
    #     == spe_condition.children[1].children[1].symbol
    # assert not spe_condition.children[0].children[1].conditioned
    # assert not spe_condition.children[1].children[1].conditioned

def test_product_disjoint_union_numerical():
    X = Id('X')
    Y = Id('Y')
    Z = Id('Z')
    spe = ProductSPE([
        X >> norm(loc=0, scale=1),
        Y >> norm(loc=0, scale=2),
        Z >> norm(loc=0, scale=2),
    ])

    for event in [
        (1/X < 4) | (X > 7),
        (2*X-3 > 0) | (Log(Y) < 3),
        ((X > 0) & (Y < 1)) | ((X < 1) & (Y < 3)) | ~(X<<{1, 2}),
        ((X > 0) & (Y < 1)) | ((X < 1) & (Y < 3)) | (Z < 0),
        ((X > 0) & (Y < 1)) | ((X < 1) & (Y < 3)) | (Z < 0) | ~(X <<{1, 3}),
    ]:
        clauses = dnf_to_disjoint_union(event)
        logps = [spe.logprob(s) for s in clauses.subexprs]
        assert allclose(logsumexp(logps), spe.logprob(event))

def test_product_disjoint_union_nominal():
    N = Id('N')
    P = Id('P')

    nationality = N >> choice({'India': 0.5, 'USA': 0.5})
    perfect = P >> choice({'Imperfect': 0.99, 'Perfect': 0.01})
    student = nationality & perfect

    condition_1 = (N << {'India'}) & (P << {'Imperfect'})
    condition_2 = (N << {'India'}) & (P << {'Perfect'})
    condition_3 = (N << {'USA'}) & (P << {'Imperfect'})
    condition_4 = (N << {'USA'}) & (P << {'Perfect'})

    event_1 = condition_1
    event_2 = condition_2 & ~condition_1
    event_3 = condition_3 & ~condition_2 & ~condition_1
    event_4 = condition_4 & ~condition_3 & ~condition_2 & ~condition_1

    assert allclose(student.prob(event_1), 0.5*0.99)
    assert allclose(student.prob(event_2), 0.5*0.01)
    assert allclose(student.prob(event_3), 0.5*0.99)
    assert allclose(student.prob(event_4), 0.5*0.01)

A = Id('A')
B = Id('B')
C = Id('C')
D = Id('D')
spe_abcd \
    = norm(loc=0, scale=1)(A) \
    & norm(loc=0, scale=1)(B) \
    & norm(loc=0, scale=1)(C) \
    & norm(loc=0, scale=1)(D)
def test_product_condition_simplify_a():
    spe = spe_abcd.condition((A > 1) | (A < -1))
    assert isinstance(spe, ProductSPE)
    assert spe_abcd.children[1] in spe.children
    assert spe_abcd.children[2] in spe.children
    assert spe_abcd.children[3] in spe.children
    idx_sum = [i for i, c in enumerate(spe.children) if isinstance(c, SumSPE)]
    assert len(idx_sum) == 1
    assert allclose(spe.children[idx_sum[0]].weights[0], -log(2))
    assert allclose(spe.children[idx_sum[0]].weights[1], -log(2))
    assert spe.children[idx_sum[0]].children[0].conditioned
    assert spe.children[idx_sum[0]].children[1].conditioned

def test_product_condition_simplify_ab():
    spe = spe_abcd.condition((A > 1) | (B < 0))
    assert isinstance(spe, ProductSPE)
    assert spe_abcd.children[2] in spe.children
    assert spe_abcd.children[2] in spe.children
    idx_sum = [i for i, c in enumerate(spe.children) if isinstance(c, SumSPE)]
    assert len(idx_sum) == 1
    spe_sum = spe.children[idx_sum[0]]
    assert isinstance(spe_sum.children[0], ProductSPE)
    assert isinstance(spe_sum.children[1], ProductSPE)
    lp0 = spe_abcd.logprob(A > 1)
    lp1 = spe_abcd.logprob((B < 0) & ~(A > 1))
    weights = lognorm([lp0, lp1])
    assert allclose(spe_sum.weights[0], weights[0])
    assert allclose(spe_sum.weights[1], weights[1])

def test_product_condition_simplify_abc():
    spe = spe_abcd.condition((A > 1) | (B < 0) | (C > 2))
    assert isinstance(spe, ProductSPE)
    assert len(spe.children) == 2
    assert spe_abcd.children[3] in spe.children
    idx_sum = [i for i, c in enumerate(spe.children) if isinstance(c, SumSPE)]
    assert len(idx_sum) == 1
    spe_sum = spe.children[idx_sum[0]]
    assert len(spe_sum.children) == 2
    lp0 = spe_abcd.logprob(A > 1)
    lp1 = spe_abcd.logprob((B < 0) & ~(A > 1))
    lp2 = spe_abcd.logprob((C > 2) & ~(B < 0) & ~(A > 1))
    weights = lognorm([lp0, lp1, lp2])
    assert allclose(spe_sum.weights[1], weights[-1])
    assert isinstance(spe_sum.children[0], ProductSPE)
    assert isinstance(spe_sum.children[0].children[0], SumSPE)
    assert isinstance(spe_sum.children[0].children[1], LeafSPE)
    assert isinstance(spe_sum.children[1], ProductSPE)
    assert isinstance(spe_sum.children[1].children[0], LeafSPE)
    assert isinstance(spe_sum.children[1].children[1], LeafSPE)
    assert isinstance(spe_sum.children[1].children[2], LeafSPE)

max_stars_repo_name: SEICS/sppl

max_issues_repo_name: SEICS/sppl

max_forks_repo_name: SEICS/sppl


text:
#!/usr/bin/env python3
# need first: pip install numba
# python fibn.py

from numba import jit
import time

@jit(nopython=True) # Set "nopython" mode for best performance, equivalent to @njit
def fib_recursive(n):
	if n <= 1: return n
	return fib_recursive(n - 1) + fib_recursive(n - 2)

def fib_iterative(n):
	a = 0
	b = 1
	for x in range(n):
		c = a + b
		a, b = b, c
	return a

t1 = time.time()
print(fib_recursive(42))
t2 = time.time()
print('Recursive took %0.3fsec.' % (t2-t1))

t1 = time.time()
print(fib_iterative(42))
t2 = time.time()
print('Iterative took %0.3fsec.' % (t2-t1))

max_stars_repo_name: glella/fib

max_issues_repo_name: glella/fib

max_forks_repo_name: glella/fib

reason: from numba


text:
import numpy as np
import numpy.linalg as la
import scipy.stats as stats

from .bayesian_single_target_tracker import BayesianSingleTargetTracker
from utils import maharanobis_distance

class NNKF(BayesianSingleTargetTracker):
    """Nearest neighbor kalman filter."""

    def __init__(self, model, covar_w, covar_v, x, covar_x, z, PG):

        self.PG = PG  # gate threshold probability
        self.gate_thresh = None

        # linear state space model of the plant
        # x_{t+1} = F_t * x_t + G_t * u_t + L * w_t
        # y_{t} = H_t * x_t + v_t
        self.model = model

        # covariance
        self.Q = covar_w
        self.R = covar_v

        # expectaion to estimate.
        self.x_pred = x
        self.x_est = None

        self.z_pred = z

        # covariance to estimate.
        self.P_pred = covar_x
        self.P_est = None

        # kalman gain
        self.Gain = None

        # innvation: z - z_pred
        self.innovation = None
        self.valid_innovations = None
        self.covar_innovation = None

    def extract_valid_measurments(self, z, z_hat, inv_covar):
        """Extract measurments in the validation region."""

        if z.shape[0] == 0:
            # shape: (0, nz)
            return np.empty_like(z), np.empty_like(z)

        distances = maharanobis_distance(z, z_hat, inv_covar)

        if not self.gate_thresh:
            self.gate_thresh = stats.chi2.ppf(q=self.PG, df=z_hat.shape[0])

        mask = distances <= self.gate_thresh
        return z[mask], distances[mask]

    def associate(self, t, z):
        """Associate the measurements with the predicted sate variable."""

        P = self.P_pred
        H = self.model.H
        R = self.R

        self.covar_innovation = H @ P @ H.T + R
        inv_covar = la.inv(self.covar_innovation)
        z, d = self.extract_valid_measurments(z, self.z_pred, inv_covar)
        print(f"step {t} : the number of valid measurments is {z.shape[0]}")

        if z.shape[0] == 0:
            self.valid_innovations = None
            self.innovation = None
        else:
            self.valid_innovations = z - self.z_pred
            self.innovation = self.valid_innovations[np.argmin(d)]

    def filtering(self):
        """Estimate the state variable."""

        # extract variables.
        x, P = self.x_pred, self.P_pred
        S = self.covar_innovation
        H = self.model.H

        # update the gain.
        W = P @ H.T @ np.linalg.inv(S)
        self.Gain = W  # store

        if self.valid_innovations is None:
            self.x_est = x
            self.P_est = P
            return

        # x(t|t), P(t|t)
        self.x_est = x + W @ self.innovation
        self.P_est = (np.eye(W.shape[0]) - W @ H) @ P

    def predict(self, t, u):
        """Predict the sate variable."""

        x, P = self.x_est, self.P_est
        F = self.model.F
        L = self.model.L
        Q = self.Q

        # x(t+1|t), P(t+1|t)
        self.x_pred = self.model.state_equation(t, x, u)
        self.P_pred = F @ P @ F.T + L @ Q @ L.T

        # z(t+1|t)
        self.z_pred = self.model.observation_equation(t+1, self.x_pred)

max_stars_repo_name: azimuth-san/pdaf-tracking

max_issues_repo_name: azimuth-san/pdaf-tracking

max_forks_repo_name: azimuth-san/pdaf-tracking

reason: import numpy,import scipy


text:
from pgfutils import save, setup_figure

setup_figure(width=1, height=0.4, columns=2)

from matplotlib import pyplot as plt
import numpy as np

t = np.linspace(0, 10, 400)
s = 0.3 * t + 2.5 * np.cos(2 * np.pi * 0.85 * t) - 0.8

plt.plot(t, s)
plt.xlim(0, 10)
plt.grid(True)

save()

max_stars_repo_name: bcbnz/matplotlib-pgfutils

max_issues_repo_name: bcbnz/matplotlib-pgfutils

max_forks_repo_name: bcbnz/matplotlib-pgfutils


text:
# performs a simple device inquiry, followed by a remote name request of each
# discovered device

import datetime as dt
import struct
import sys
import time
import logging
import warnings
import binascii

import bluetooth
import bluetooth._bluetooth as bluez
import numpy as np
from scipy.interpolate import interp1d
from tabulate import tabulate
from .utils import RealTimePlot, spin, rssi_to_colour_str
from ._base import show_header, term, \
    MAJOR_CLASSES, MINOR_CLASSES, SERVICES
from ._exceptions import AdapterUnaccessibleError

DEVICE_ID = 0
LOADING_HANDLER = None
VALUES_PER_FRAME = 50
CATEGORY_VALUES = [0, -10, -30, -50, -70]
OUT_OF_RANGE = (-300, -200)
NAME_DICT = dict()
EXTRA_INFO_DICT = dict()
CLASS_DICT = dict()

def printpacket(pkt):
    for c in pkt:
        sys.stdout.write("%02x " % struct.unpack("B", c)[0])

def get_device_extra(addr):
    """
    return the extra device info from the global dict
    """
    try:
        _info_dict = EXTRA_INFO_DICT[addr]
        major_device = _info_dict.get("major_device", "XXXX")
        minor_device = _info_dict.get("minor_device", "XXXX")
        services = _info_dict.get("services", "XXXX")
        return [major_device, minor_device, services]
    except:
        EXTRA_INFO_DICT[addr] = {
            "major_device": "",
            "minor_device": "",
            "services": "",
        }
        return ["XXXX", "XXXX", "XXXX"]

def populate_info_dict():
    """
    call to populate the global info dictionary
    """

    # extract hex value dictionary
    hex_dict = dict()
    for i in CLASS_DICT:
        hex_dict[i] = "%X" % CLASS_DICT[i]
    # check against odd length hex values
    for i in hex_dict:
        if len(hex_dict[i]) % 2 != 0:
            hex_dict[i] = "0{}".format(hex_dict[i])

    # initialize entries in EXTRA_INFO_DICT using vars
    for i in hex_dict:
        if i not in EXTRA_INFO_DICT:
            EXTRA_INFO_DICT[i] = {
                "major_device": "",
                "minor_device": "",
                "services": "",
            }
    # extract byte dictionary
    byte_dict = {i: binascii.unhexlify(hex_dict[i]) for i in hex_dict}
    # extract bit dictionary using big byte-decode
    bit_dict = {i: bin(int.from_bytes(byte_dict[i], 'big')) for i in byte_dict}

    # service bits
    serv_numbers = [13, 16, 17, 18, 19, 20, 21, 22, 23]
    services = ""
    for i in bit_dict:
        bit_stream = bit_dict[i]
        major_stream = bit_stream[-13:-8]
        major_class = MAJOR_CLASSES.get(major_stream, "XXXX")
        minor_class = "XXXX"
        if major_class in ("Miscellaneous", "Device code not specified"):
            minor_class = "XXXX"
            services = "XXXX"
        elif major_class == "Computer":
            minor_class = MINOR_CLASSES[major_class].get(bit_stream[-8:-2], "XXXX")
        elif major_class == "Phone":
            minor_class = MINOR_CLASSES[major_class].get(bit_stream[-8:-2], "XXXX")
        elif major_class == "LAN/Network Access Point":
            minor_class = MINOR_CLASSES[major_class].get(bit_stream[-8:-5], "XXXX")
        elif major_class == "Audio/Video":
            minor_class = MINOR_CLASSES[major_class].get(bit_stream[-8:-2], "XXXX")
        elif major_class == "Peripheral":
            minor_class = MINOR_CLASSES[major_class].get(bit_stream[-8:-6], "XXXX")
        elif major_class == "Imaging":
            minor_class = MINOR_CLASSES[major_class].get(bit_stream[-8:-4], "XXXX")
        elif major_class == "Wearable":
            minor_class = MINOR_CLASSES[major_class].get(bit_stream[-8:-2], "XXXX")
        elif major_class == "Toy":
            minor_class = MINOR_CLASSES[major_class].get(bit_stream[-8:-2], "XXXX")
        elif major_class == "Health":
            minor_class = MINOR_CLASSES[major_class].get(bit_stream[-8:-2], "XXXX")

        # parse services logic, appending each available service
        serv = ""
        if services != "XXXX":
            # run from 13 to 23 excluding (14, 15)
            for x in serv_numbers:
                # appending 0 before odd-numbered hex values sometimes causes
                # the bit_stream[-x] query to hit the 'b' flag of the bit stream
                try:
                    # if bit at position -x is 1, append service
                    if bool(int(bit_stream[-x])):
                        serv += "{}|".format(SERVICES[str(x)])
                except Exception as e:
                    pass
            services = serv

        EXTRA_INFO_DICT[i] = {
            "major_device": major_class,
            "minor_device": minor_class,
            "services": services,
        }

def read_inquiry_mode(sock):
    """returns the current mode, or -1 on failure"""
    # save current filter
    old_filter = sock.getsockopt(bluez.SOL_HCI, bluez.HCI_FILTER, 14)

    # Setup socket filter to receive only events related to the
    # read_inquiry_mode command
    flt = bluez.hci_filter_new()
    opcode = bluez.cmd_opcode_pack(bluez.OGF_HOST_CTL,
                                   bluez.OCF_READ_INQUIRY_MODE)
    bluez.hci_filter_set_ptype(flt, bluez.HCI_EVENT_PKT)
    bluez.hci_filter_set_event(flt, bluez.EVT_CMD_COMPLETE)
    bluez.hci_filter_set_opcode(flt, opcode)
    sock.setsockopt(bluez.SOL_HCI, bluez.HCI_FILTER, flt)

    try:
        # first read the current inquiry mode.
        bluez.hci_send_cmd(sock, bluez.OGF_HOST_CTL, bluez.OCF_READ_INQUIRY_MODE)
    except bluez.error as e:
        raise AdapterUnaccessibleError("Are you sure this a bluetooth 1.2 device? \nTurn On Your Bluetooth")

    pkt = sock.recv(255)

    status, mode = struct.unpack("xxxxxxBB", pkt)
    if status != 0: mode = -1

    # restore old filter
    sock.setsockopt(bluez.SOL_HCI, bluez.HCI_FILTER, old_filter)
    return mode

def write_inquiry_mode(sock, mode):
    """returns 0 on success, -1 on failure"""
    # save current filter
    old_filter = sock.getsockopt(bluez.SOL_HCI, bluez.HCI_FILTER, 14)

    # Setup socket filter to receive only events related to the
    # write_inquiry_mode command
    flt = bluez.hci_filter_new()
    opcode = bluez.cmd_opcode_pack(bluez.OGF_HOST_CTL,
                                   bluez.OCF_WRITE_INQUIRY_MODE)
    bluez.hci_filter_set_ptype(flt, bluez.HCI_EVENT_PKT)
    bluez.hci_filter_set_event(flt, bluez.EVT_CMD_COMPLETE)
    bluez.hci_filter_set_opcode(flt, opcode)
    sock.setsockopt(bluez.SOL_HCI, bluez.HCI_FILTER, flt)

    # send the command!
    bluez.hci_send_cmd(sock, bluez.OGF_HOST_CTL,
                       bluez.OCF_WRITE_INQUIRY_MODE, struct.pack("B", mode))

    pkt = sock.recv(255)

    status = struct.unpack("xxxxxxB", pkt)[0]

    # restore old filter
    sock.setsockopt(bluez.SOL_HCI, bluez.HCI_FILTER, old_filter)
    if status != 0: return -1
    return 0

def device_inquiry_with_with_rssi(sock, show_name=False, show_extra_info=False, color=True, ret_table=False):
    global LOADING_HANDLER

    # save current filter
    old_filter = sock.getsockopt(bluez.SOL_HCI, bluez.HCI_FILTER, 14)

    # perform a device inquiry on bluetooth device #0
    # The inquiry should last 8 * 1.28 = 10.24 seconds
    # before the inquiry is performed, bluez should flush its cache of
    # previously discovered devices
    flt = bluez.hci_filter_new()
    bluez.hci_filter_all_events(flt)
    bluez.hci_filter_set_ptype(flt, bluez.HCI_EVENT_PKT)
    sock.setsockopt(bluez.SOL_HCI, bluez.HCI_FILTER, flt)

    duration = 1
    max_responses = 255
    cmd_pkt = struct.pack("BBBBB", 0x33, 0x8b, 0x9e, duration, max_responses)
    # TODO Optimize code for performance
    # update the global device name dictionary before sending hci cmd(which changes mode)
    headers = ["Name", "MAC Address", "RSSI"]
    data = []
    results = []
    if show_extra_info or show_name:
        devices = bluetooth.discover_devices(lookup_names=True, lookup_class=True)
        if show_name:
            update_dict = {i[0]: i[1] for i in devices}
            NAME_DICT.update(update_dict)
        if show_extra_info:
            update_dict = {i[0]: i[2] for i in devices}
            CLASS_DICT.update(update_dict)
            headers.extend(["Major Dev Class", "Minor Dev Class", "Services"])
            populate_info_dict()
    bluez.hci_send_cmd(sock, bluez.OGF_LINK_CTL, bluez.OCF_INQUIRY, cmd_pkt)

    done = False
    while not done:
        pkt = sock.recv(255)
        ptype, event, plen = struct.unpack("BBB", pkt[:3])
        if event == bluez.EVT_INQUIRY_RESULT_WITH_RSSI:
            pkt = pkt[3:]
            nrsp = bluetooth.get_byte(pkt[0])
            for i in range(nrsp):
                # get human readable addr
                addr = bluez.ba2str(pkt[1 + 6 * i:1 + 6 * i + 6])
                rssi = bluetooth.byte_to_signed_int(
                    bluetooth.get_byte(pkt[1 + 13 * nrsp + i]))
                # retrieve device name, or assign address as name
                try:
                    name = NAME_DICT[addr]
                except:
                    name = addr
                results.append((addr, rssi, name))
                if color:
                    data.append([name, addr, rssi_to_colour_str(rssi)])
                else:
                    data.append([name, addr, rssi])
                if show_extra_info:
                    extra_info = get_device_extra(addr)
                    # extend last data list with extra info
                    data[-1].extend(extra_info)
        elif event == bluez.EVT_INQUIRY_COMPLETE:
            done = True
        elif event == bluez.EVT_CMD_STATUS:
            status, ncmd, opcode = struct.unpack("BBH", pkt[3:7])
            if status != 0:
                print("uh oh...")
                printpacket(pkt[3:7])
                done = True
        elif event == bluez.EVT_INQUIRY_RESULT:
            pkt = pkt[3:]
            nrsp = bluetooth.get_byte(pkt[0])
            for i in range(nrsp):
                addr = bluez.ba2str(pkt[1 + 6 * i:1 + 6 * i + 6])
                results.append((addr, -1, "UNK"))
                print("[%s] (no RRSI)" % addr)
        else:
            logging.debug("unrecognized packet type 0x%02x" % ptype)
        logging.debug("event %s", event)

    # restore old filter
    sock.setsockopt(bluez.SOL_HCI, bluez.HCI_FILTER, old_filter)
    # if ordered to return a table by analyze_all, ignore other sequence
    if ret_table:
        if len(results) < 1:
            return ((None, headers))
        return ((data, headers))
    else:
        # print all the data at once since blessings clears the screen just before
        if len(results) >= 1:
            # terminate concurrent loading handler
            if bool(LOADING_HANDLER):
                LOADING_HANDLER.terminate()
            show_header("BLUETOOTH")
            print(tabulate(data, headers=headers, disable_numparse=True))
        else:
            # LOADING_HANDLER = spin(before="Searching",
            #                    after="\nNo devices found in nearby range")
            LOADING_HANDLER.terminate()
            LOADING_HANDLER = spin(before="No BT devices in nearby range")
        return results

def animate(i, ax, plt, val_dict, xs, sock, show_name=False, show_extra_info=False):
    """
    Instance function to create matplotlib graph

    """
    # TODO Hide/cutout devices with rssi < -200
    results = device_inquiry_with_with_rssi(sock, show_name=show_name)
    # append datetime string as a float to represent time axis
    xs.append(float(dt.datetime.now().strftime("%H.%M%S")))
    NAME_DICT.update({i[0]: i[2] for i in results})
    for i in results:
        try:
            # check for dict key if it exists
            affect_list = val_dict[i[0]]
            affect_list.append(i[1])
        except Exception as e:
            # create new list with prior values out of range
            val_dict[i[0]] = list()
            val_dict[i[0]].extend([np.random.random_integers(*OUT_OF_RANGE) \
                                   for i in range(len(xs))])

    ax.clear()
    # limit both axis to VALUES_PER_FRAME values at a time maximum
    xs = xs[-VALUES_PER_FRAME:]
    for i in val_dict:
        device_name = NAME_DICT.get(i, "XXXX")
        val_dict[i] = val_dict[i][-VALUES_PER_FRAME:]
        # if device has dissapeared, append zeros to make up length
        if len(val_dict[i]) < len(xs):
            val_dict[i].extend([np.random.random_integers(*OUT_OF_RANGE) \
                                for i in range(len(xs) - len(val_dict[i]))])
        # smoothen out x axis before display
        x = np.array(xs)
        y = np.array(val_dict[i])
        x_new = np.linspace(x.min(), x.max(), 500)
        # check if points are enough to interpolate on and use box(nearest) interpolation
        # to display levels to this
        if len(x) > 2:
            f = interp1d(x, y, kind='nearest')
            y_smooth = f(x_new)
            # plot smooth plot with scatter point plots
            ax.plot(x_new, y_smooth, label=device_name)
        else:
            ax.plot(xs, y, label=device_name)
        # ax.scatter(xs, y)
    # display legend, attempt to supress warnings
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        ax.legend()

    plt.xticks([])
    plt.ylim(-100, 0)
    plt.title("Bluetooth Devices RSSI against time")
    plt.ylabel("BT RSSI")
    plt.xlabel("Time")

def bluelyze(**kwargs):
    global LOADING_HANDLER
    show_graph = kwargs.pop("graph")
    show_name = kwargs.pop("show_name")
    show_extra_info = kwargs.pop("show_extra_info")
    analyze_all = kwargs.pop("analyze_all")
    _color = kwargs.get("color", True)

    try:
        sock = bluez.hci_open_dev(DEVICE_ID)
    except:
        print("Error accessing bluetooth device...\n"
              "Confirm if your bluetooth device is correctly installed and try again")
        sys.exit(1)

    try:
        mode = read_inquiry_mode(sock)
        logging.debug("current inquiry mode is %d" % mode)

        if mode != 1:
            logging.debug("writing inquiry mode...")
            try:
                result = write_inquiry_mode(sock, 1)
            except Exception as e:
                print("error writing inquiry mode.  Are you sure you're root?")
                print(e)
                sys.exit(1)
            if result != 0:
                print("error while setting inquiry mode")
            logging.debug("result: %d" % result)

        if analyze_all:
            return device_inquiry_with_with_rssi(sock, show_name, show_extra_info, _color, ret_table=True)
        else:
            print(term.clear())
            show_header("BLUETOOTH")
            LOADING_HANDLER = spin(before="Initializing...")

            if show_graph:
                # create general figure object
                xs = []
                results = device_inquiry_with_with_rssi(sock, show_name, show_extra_info, _color)
                # initialize dictionary to store real time values of devices
                val_dict = {key: list() for key, value, name in results}
                realtimeplot = RealTimePlot(
                    func=animate,
                    func_args=(val_dict, xs, sock, show_name, show_extra_info, _color),
                )
                realtimeplot.animate()

            else:
                while True:
                    device_inquiry_with_with_rssi(sock, show_name, show_extra_info, _color)

    except (Exception, bluez.error) as e:
        if LOADING_HANDLER:
            LOADING_HANDLER.terminate()
        # Analyze implements its own error handler
        if analyze_all:
            raise (e)
        else:
            logging.debug("error reading inquiry mode.  ")
            show_header("BLUETOOTH")
            print("Are you sure this a bluetooth 1.2 device? \nTurn On Your Bluetooth")
            logging.debug(e)
            sys.exit(1)

max_stars_repo_name: bisoncorps/signalum

max_issues_repo_name: bisoncorps/signalum

max_forks_repo_name: bisoncorps/signalum


text:
from __future__ import print_function
"""
    Model evaluation on VOC for vggSSD separately
    Execute: python3 eval_voc_vrmSSD.py --trained_model weights/_your_trained_SSD_model_.pth

    Model evaluation on VOC for resnetSSD separately
    Execute: python3 eval_voc_vrmSSD.py --use_res --trained_model weights/_your_trained_SSD_model_.pth

    Model evaluation on VOC for mobileSSD v1 separately
    Execute: python3 eval_voc_vrmSSD.py --use_m1 --trained_model weights/_your_trained_SSD_model_.pth

    Model evaluation on VOC for mobileSSD v2 separately
    Execute: python3 eval_voc_vrmSSD.py --use_m2 --trained_model weights/_your_trained_SSD_model_.pth

    (Take care of different versions of .pth file, can be solved by changing state_dict)

    Author: xuhuahuang as intern in YouTu 07/2018
    Status: checked
"""

import torch
import torch.nn as nn
import torch.backends.cudnn as cudnn
from torch.autograd import Variable
from data import *
import torch.utils.data as data

from models.SSD_vggres import build_ssd
from models.SSD_mobile import build_mssd

import sys
import os
import time
import argparse
import numpy as np
import pickle
import cv2
cv2.setNumThreads(0) # pytorch issue 1355: possible deadlock in DataLoader
# OpenCL may be enabled by default in OpenCV3;
# disable it because it because it's not thread safe and causes unwanted GPU memory allocations
cv2.ocl.setUseOpenCL(False)

if sys.version_info[0] == 2:
    import xml.etree.cElementTree as ET
else:
    import xml.etree.ElementTree as ET

def str2bool(v):
    return v.lower() in ("yes", "true", "t", "1")

parser = argparse.ArgumentParser(
    description='Single Shot MultiBox Detector Evaluation')
parser.add_argument('--trained_model',
                    default='weights/ssd300_mAP_77.43_v2.pth', type=str,
                    help='Trained state_dict file path to open')
parser.add_argument('--save_folder', default='eval/', type=str,
                    help='File path to save results')
parser.add_argument('--confidence_threshold', default=0.01, type=float,
                    help='Detection confidence threshold')
# 200 in SSD paper, 200 for COCO, 300 for VOC
parser.add_argument('--max_per_image', default=200, type=int,
                    help='Top number of detections kept per image, further restrict the number of predictions to parse')
parser.add_argument('--cuda', default=True, type=str2bool,
                    help='Use cuda to train model')
parser.add_argument('--voc_root', default= VOC_ROOT,# XL_ROOT, for VOC_xlab_products dataset
                    help='Location of XL root directory')
parser.add_argument('--cleanup', default=True, type=str2bool,
                    help='Cleanup and remove results files following eval')
# for resnet backbone
parser.add_argument("--use_res", dest="use_res", action="store_true")
parser.set_defaults(use_res=False)
# for mobilev1 backbone
parser.add_argument("--use_m1", dest="use_m1", action="store_true")
parser.set_defaults(use_m1=False)
# for mobilev2 backbone
parser.add_argument("--use_m2", dest="use_m2", action="store_true")
parser.set_defaults(use_m2=False)

args = parser.parse_args()

if not os.path.exists(args.save_folder):
    os.mkdir(args.save_folder)

if torch.cuda.is_available():
    if args.cuda:
        torch.set_default_tensor_type('torch.cuda.FloatTensor')
    if not args.cuda:
        print("WARNING: It looks like you have a CUDA device, but aren't using \
              CUDA.  Run with --cuda for optimal eval speed.")
        torch.set_default_tensor_type('torch.FloatTensor')
else:
    torch.set_default_tensor_type('torch.FloatTensor')

set_type = 'test'
cfg = voc #xl, for VOC_xlab_products dataset

# test function for vggSSD
"""
    Args:
        save_folder: the eval results saving folder
        net: test-type ssd net
        dataset: validation dataset
        transform: BaseTransform - not used here
"""
def test_net(save_folder, net, cuda,
             testset, transform, max_per_image=300, thresh=0.05):

    if not os.path.exists(save_folder):
        os.mkdir(save_folder)

    num_images = len(testset)
    num_classes = cfg['num_classes']
    # all detections are collected into:
    #    all_boxes[cls][image] = N x 5 array of detections in
    #    (x1, y1, x2, y2, score)
    all_boxes = [[[] for _ in range(num_images)]
                 for _ in range(num_classes)]

    # timers
    _t = {'im_detect': Timer(), 'misc': Timer()}
    #output_dir = get_output_dir('ssd300_120000', set_type) #directory storing output results
    #det_file = os.path.join(output_dir, 'detections.pkl') #file storing output result under output_dir
    det_file = os.path.join(save_folder, 'detections.pkl')

    for i in range(num_images):
        im, gt, h, w = testset.pull_item(i) # include BaseTransform inside

        x = Variable(im.unsqueeze(0)) #insert a dimension of size one at the dim 0
        if cuda:
            x = x.cuda()

        _t['im_detect'].tic()
        detections = net(x=x, test=True).data # get the detection results
        detect_time = _t['im_detect'].toc(average=False) #store the detection time

        # skip j = 0, because it's the background class
        for j in range(1, detections.size(1)): # for every class
            dets = detections[0, j, :]#size( ** , 5)
            mask = dets[:, 0].gt(0.).expand(5, dets.size(0)).t()
            dets = torch.masked_select(dets, mask).view(-1, 5)
            if dets.dim() == 0:
                continue
            #if dets.size(0) == 0:
            #    continue
            boxes = dets[:, 1:]
            boxes[:, 0] *= w
            boxes[:, 2] *= w
            boxes[:, 1] *= h
            boxes[:, 3] *= h
            scores = dets[:, 0].cpu().numpy()
            cls_dets = np.hstack((boxes.cpu().numpy(),
                                  scores[:, np.newaxis])).astype(np.float32,
                                                                 copy=False)
            all_boxes[j][i] = cls_dets #[class][imageID] = 1 x 5 where 5 is box_coord + score
        if (i + 1) % 100 == 0:
            print('im_detect: {:d}/{:d} {:.3f}s'.format(i + 1,
                                                        num_images, detect_time))

    #write the detection results into det_file
    with open(det_file, 'wb') as f:
        pickle.dump(all_boxes, f, pickle.HIGHEST_PROTOCOL)

    print('Evaluating detections')
    APs,mAP = testset.evaluate_detections(all_boxes, save_folder)

if __name__ == '__main__':
    # load net
    num_classes = cfg['num_classes']
    if args.use_res:
        net = build_ssd('test', cfg, 300, num_classes, base='resnet', max_per_image = args.max_per_image) # initialize SSD (resnet)
    elif args.use_m1:
        net = build_mssd('test', cfg, 300, num_classes, base='m1', max_per_image = args.max_per_image) # backbone network is m1
    elif args.use_m2:
        net = build_mssd('test', cfg, 300, num_classes, base='m2', max_per_image = args.max_per_image) # backbone network is m2
    else:
        net = build_ssd('test', cfg, 300, num_classes, base='vgg', max_per_image = args.max_per_image) # initialize SSD (vgg)
    # if you want to eval SSD from original version ssd.pytorch because self.vgg was changed to self.base
    '''
    # load resume SSD network
    state_dict = torch.load(args.trained_model)
    from collections import OrderedDict
    new_state_dict = OrderedDict()
    for k, v in state_dict.items():
        # change from self.vgg to self.base
        head = k[:4]
        if head == 'vgg.':
            name = 'base.' + k[4:]
        else:
            name = k
        new_state_dict[name] = v
    net.load_state_dict(new_state_dict)
    '''
    net.load_state_dict(torch.load(args.trained_model))
    net.eval()
    print('Finished loading model!')
    # load data
#    dataset = XLDetection(args.voc_root, [set_type], # for VOC_xlab_products dataset
#                           BaseTransform(300, cfg['dataset_mean']),
#                           XLAnnotationTransform())
    dataset = VOCDetection(args.voc_root, [('2007', set_type)],
                           BaseTransform(300, cfg['dataset_mean']),
                           VOCAnnotationTransform())
    if args.cuda:
        net = net.cuda()
        cudnn.benchmark = True
    # evaluation
    test_net(args.save_folder, net, args.cuda, dataset,
             BaseTransform(net.size, cfg['dataset_mean']), args.max_per_image,
             thresh=args.confidence_threshold)


text:
import pyglet
import numpy as np

class DmControlViewer:
    def __init__(self, width, height, depth=False):
        self.window = pyglet.window.Window(width=width, height=height, display=None)
        self.width = width
        self.height = height

        self.depth = depth

        if depth:
            self.format = 'RGB'
            self.pitch = self.width * -3
        else:
            self.format = 'RGB'
            self.pitch = self.width * -3

    def update(self, pixel):
        self.window.clear()
        self.window.switch_to()
        self.window.dispatch_events()
        if self.depth:
            pixel = np.dstack([pixel.astype(np.uint8)] * 3)
        pyglet.image.ImageData(self.width, self.height, self.format, pixel.tobytes(), pitch=self.pitch).blit(0, 0)
        self.window.flip()

    def close(self):
        self.window.close()


text:
#!/usr/bin/env python

# import sys,os
from math import floor
from PyQt5 import QtWidgets, QtCore, uic
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
from matplotlib.backends.backend_qt5agg import NavigationToolbar2QT as NavigationToolbar
from matplotlib.figure import Figure
import numpy as np
from . artists import TemplateSourcePlot

class _Figure(FigureCanvas):
	def __init__(self, parent=None):
		self.figure = Figure(dpi=100)
		self.ax     = self.figure.add_axes([0,0,1,1])
		super().__init__(self.figure)
		self.setParent(parent)
		super().setSizePolicy(QtWidgets.QSizePolicy.Expanding, QtWidgets.QSizePolicy.Expanding)
		super().updateGeometry()
		self.setFocusPolicy( QtCore.Qt.ClickFocus )
		self.setFocus()

		policy = self.sizePolicy()
		policy.setRetainSizeWhenHidden(True)
		self.setSizePolicy(policy)

	def reset(self):
		self.ax.figure.clf()
		self.ax     = self.figure.add_axes([0,0,1,1])
		self.ax.figure.canvas.draw()
		self.tsplot = None

	def set_visible(self, visible):
		self.setVisible(visible)

class FigureLandmarksRegistered(_Figure):

	color_source_selected    = (1, 0.56, 0)
	color_source_unselected  = '0.7'

	def _init(self, y0=None, y=None):
		self.y0           = np.random.randn(101) if (y0 is None) else y0
		self.y            = np.random.randn(8,101) if (y is None) else y
		self.h0           = self.ax.plot(self.y0,  lw=4,   color='k')[0]
		self.h            = self.ax.plot(self.y.T, lw=0.5, color=self.color_source_unselected)
		self.ax.figure.canvas.draw()

	def clear_source_selection(self, update=True):
		for h in self.h:
			h.set_color( self.color_source_unselected )
			h.set_linewidth(0.5)
		if update:
			self.ax.figure.canvas.draw()

	def plot(self, y):
		self.ax.plot( y.T )

	def get_warped_sources_as_array(self):
		return np.array([hh.get_ydata()  for hh in self.h])

	def set_source_selected(self, ind):
		self.clear_source_selection(update=False)
		self.h[ind].set_color( self.color_source_selected )
		self.h[ind].set_linewidth( 1 )
		self.ax.figure.canvas.draw()

	def set_sources_visible(self, visible):
		[h.set_visible(visible) for h in self.h]

	def set_sources_visible_bool_list(self, visible):
		for h,v in zip(self.h, visible):
			h.set_visible(v)

	def set_template_visible(self, visible):
		self.h0.set_visible( visible )

class FigureLandmarks(_Figure):

	landmark_deleted                     = QtCore.pyqtSignal(int)

	def __init__(self, parent=None):
		super().__init__(parent)
		self.J              = None
		self.Q              = None
		self.tsplot         = None

	def _init(self, template, sources, lmtemplate=None, lmsources=None):
		J,Q               = sources.shape
		self.J            = J
		self.Q            = Q
		ytemplate         = np.random.randn(101) if (template is None) else template
		ysources          = np.random.randn(8,101) if (sources is None) else sources
		lmtemplate        = [25, 50, 75] if (lmtemplate is None) else lmtemplate
		lmsources         = np.vstack([[25]*J, [50]*J, [75]*J]).T if (lmsources is None) else lmsources

		self.tsplot       = TemplateSourcePlot(self.ax, self.Q)
		self.tsplot.set_template( ytemplate, lmtemplate )
		self.tsplot.set_sources( ysources, lmsources )
		self.tsplot.init_legend()

		self.tsplot.point_deleted.connect(self.on_point_deleted)

		self.curve_left_click_selected    = self.tsplot.curve_left_click_selected
		self.curve_arrow_selected         = self.tsplot.curve_arrow_selected
		self.dragging_stopped             = self.tsplot.dragging_stopped
		self.template_point_dragged       = self.tsplot.template_point_dragged
		self.source_point_dragged         = self.tsplot.source_point_dragged
		self.template_landmark_added      = self.tsplot.template_landmark_added
		self.toggle_unselected_visible    = self.tsplot.toggle_unselected_visible
		self.update_flagged_source_colors = self.tsplot.update_flagged_source_colors

		self.maxlandmarksreached       = self.tsplot.maxlandmarksreached

		self.get_active_source         = self.tsplot.get_active_source
		self.reset_source_landmarks    = self.tsplot.reset_source_landmarks
		self.set_template_active       = self.tsplot.set_template_active
		self.set_template_landmark_adding_enabled        = self.tsplot.set_template_landmark_adding_enabled
		self.set_template_landmark_dragging_enabled      = self.tsplot.set_template_landmark_dragging_enabled

		self.ax.figure.canvas.draw()

	def get_source_landmarks_objects(self):
		return [s.landmarks for s in self.tsplot.sources]

	def get_template_landmarks_object(self):
		return self.tsplot.template.landmarks

	def on_point_deleted(self, obj, ind):
		self.landmark_deleted.emit(ind)

	def set_next_selected(self):
		self.tsplot.set_next_selected()
		# self.curve_selected.emit()
	def set_previous_selected(self):
		self.tsplot.set_previous_selected()
	def select_source(self, ind):
		self.tsplot.set_source_active(ind)
	def select_source_by_table_header_click(self, ind):
		self.tsplot.toggle_unselected_visible()
		self.tsplot.set_source_active(ind)
		self.tsplot.toggle_unselected_visible()

	def set_source_landmarks_visible(self, visible):
		[source.landmarks.set_visible(visible) for source in self.tsplot.sources]

	def set_sources_visible(self, visible):
		[source.set_visible(visible) for source in self.tsplot.sources]

	def set_template_locked(self, locked):
		self.set_template_landmark_adding_enabled(not locked)
		self.set_template_landmark_dragging_enabled(not locked)
		self.set_sources_visible(locked)
		self.set_source_landmarks_visible(False)
		self.ax.figure.canvas.draw()

	def set_template_selected(self):
		self.tsplot.set_template_active()

	def set_template_visible(self, visible):
		self.tsplot.set_template_visible(visible)


text:
from deepar.dataset import Dataset
import numpy as np
import pandas as pd
import logging

logger = logging.getLogger('deepar')

class MockTs(Dataset):
    """
    This class generates 'mock' time series data of the form (y = t * np.sin(t/6) / 3 +np.sin(t*2))
    Created mainly for showcase/testing purpose
    """
    def __init__(self, t_min=0, t_max=30, resolution=.1):
        self.t_min = t_min
        self.t_max = t_max
        self.resolution = resolution
        self.data = True

    @staticmethod
    def _time_series(t):
        return t * np.sin(t/6) / 3 + np.sin(t*2)

    def next_batch(self, batch_size, n_steps):
        """
        Generate next batch (x, y), generate y by lagging x (1 step)
        """
        t0 = np.random.rand(batch_size, 1) * (self.t_max - self.t_min - n_steps * self.resolution)
        Ts = t0 + np.arange(0., n_steps + 1) * self.resolution
        ys = self._time_series(Ts)
        return ys[:, :-1].reshape(-1, n_steps, 1), ys[:, 1:].reshape(-1, n_steps, 1)

    @property
    def mock_ts(self):
        """
        Return the data used for training (ranging from self.t_min and self.t_max, with resolution self.resolution)
        :return: a Numpy array
        """
        t_list = [self.t_min]
        results = [self._time_series(t_list[0])]
        while t_list[-1] < self.t_max:
            t_list.append(t_list[-1] + self.resolution)
            results.append(self._time_series(t_list[-1]))
        return results

    def generate_test_data(self, n_steps):
        """
        Generate test data outside of the training set (t > self.t_max)
        :param n_steps:
        :return:
        """
        t_list = [self.t_max + self.resolution]
        results = [self._time_series(t_list[0])]
        for i in range(1, n_steps):
            t_list.append(t_list[-1] + self.resolution)
            results.append(self._time_series(t_list[-1]))
        return results

class TimeSeries(Dataset):
    def __init__(self, pandas_df, one_hot_root_list=None, grouping_variable='category', scaler=None):
        super().__init__()
        self.data = pandas_df
        self.one_hot_root_list = one_hot_root_list
        self.grouping_variable = grouping_variable
        if self.data is None:
            raise ValueError('Must provide a Pandas df to instantiate this class')
        self.scaler = scaler

    def _one_hot_padding(self, pandas_df, padding_df):
        """
        Util padding function
        :param padding_df:
        :param one_hot_root_list:
        :return:
        """
        for one_hot_root in self.one_hot_root_list:
            one_hot_columns = [i for i in pandas_df.columns   # select columns equal to 1
                               if i.startswith(one_hot_root) and pandas_df[i].values[0] == 1]
            for col in one_hot_columns:
                padding_df[col] = 1
        return padding_df

    def _pad_ts(self, pandas_df, desired_len, padding_val=0):
        """
        Add padding int to the time series
        :param pandas_df:
        :param desired_len: (int)
        :param padding_val: (int)
        :return: X (feature_space), y
        """
        pad_length = desired_len - pandas_df.shape[0]
        padding_df = pd.concat([pd.DataFrame({col: padding_val for col in pandas_df.columns},
                                             index=[i for i in range(pad_length)])])

        if self.one_hot_root_list:
            padding_df = self._one_hot_padding(pandas_df, padding_df)

        return pd.concat([padding_df, pandas_df]).reset_index(drop=True)

    @staticmethod
    def _sample_ts(pandas_df, desired_len):
        """

        :param pandas_df: input pandas df with 'target' columns e features
        :param desired_len: desired sample length (number of rows)
        :param padding_val: default is 0
        :param initial_obs: how many observations to skip at the beginning
        :return: a pandas df (sample)
        """
        if pandas_df.shape[0] < desired_len:
            raise ValueError('Desired sample length is greater than df row len')
        if pandas_df.shape[0] == desired_len:
            return pandas_df

        start_index = np.random.choice([i for i in range(0, pandas_df.shape[0] - desired_len + 1)])
        return pandas_df.iloc[start_index: start_index+desired_len, ]

    def next_batch(self, batch_size, n_steps,
                   target_var='target', verbose=False,
                   padding_value=0):
        """
        :param batch_size: how many time series to be sampled in this batch (int)
        :param n_steps: how many RNN cells (int)
        :param target_var: (str)
        :param verbose: (boolean)
        :param padding_value: (float)
        :return: X (feature space), y
        """

        # Select n_batch time series
        groups_list = self.data[self.grouping_variable].unique()
        np.random.shuffle(groups_list)
        selected_groups = groups_list[:batch_size]
        input_data = self.data[self.data[self.grouping_variable].isin(set(selected_groups))]

        # Initial padding for each selected time series to reach n_steps
        sampled = []
        for cat, cat_data in input_data.groupby(self.grouping_variable):
                if cat_data.shape[0] < n_steps:
                    sampled_cat_data = self._pad_ts(pandas_df=cat_data,
                                                    desired_len=n_steps,
                                                    padding_val=padding_value)
                else:
                    sampled_cat_data = self._sample_ts(pandas_df=cat_data,
                                                       desired_len=n_steps)
                sampled.append(sampled_cat_data)
                if verbose:
                    logger.debug('Sampled data for {}'.format(cat))
                    logger.debug(sampled_cat_data)
        rnn_output = pd.concat(sampled).drop(columns=self.grouping_variable).reset_index(drop=True)

        if self.scaler:
            batch_scaler = self.scaler()
            n_rows = rnn_output.shape[0]
            # Scaling will have to be extended to handle multiple variables!
            rnn_output['feature_1'] = rnn_output.feature_1.astype('float')
            rnn_output[target_var] = rnn_output[target_var].astype('float')

            rnn_output['feature_1'] = batch_scaler.fit_transform(rnn_output.feature_1.values.reshape(n_rows, 1)).reshape(n_rows)
            rnn_output[target_var] = batch_scaler.fit_transform(rnn_output[target_var].values.reshape(n_rows, 1)).reshape(n_rows)

        return rnn_output.drop(target_var, 1).values.reshape(batch_size, n_steps, -1), \
               rnn_output[target_var].values.reshape(batch_size, n_steps, 1)

max_stars_repo_name: mannimm/deepar

max_issues_repo_name: mannimm/deepar

max_forks_repo_name: mannimm/deepar