text:
import torch
import numpy as np
from time import time

class FIM:
  """Fisher information matrix details"""
  def __init__(self):
    ##Following should be set by subclasses if needed
    self.npars = None ## Positive integer: used in makeA (always needed)
    self.nsamples = None ## Positive integer: used in estimate_expected_FIM
    self.prior = None ## PyTorch distribution: used in estimate_expected_FIM

  def estimate_FIM(self, theta, design):
    """Return estimates of Fisher information matrices

    `theta` is a `self.nsamples` x `self.npars` tensor of parameters
    `design` is a tensor of designs, batched along dimension 0

    The output is a `design.shape[0]` x `self.nsamples` x `self.npars` x `self.npars` tensor
    """
    raise NotImplementedError

  def estimate_expected_FIM(self, design):
    """Return estimates of the expected Fisher information

    `design` is a tensor of designs, batched along dimension 0

    The output is a `design.shape[0]` x `self.npars` x `self.npars` tensor
    """
    ## Default code performs Monte Carlo estimation.
    ## Subclasses can override and use other approaches if desired.
    theta = self.prior.sample((self.nsamples,))
    fim = self.estimate_FIM(theta, design)
    return torch.mean(fim, dim=1)

  def estimateK(self, design, A):
    """Return unbiased estimates of K, the adversarial design objective

    `design` is a tensor of designs, batched along dimension 0
    `A` is a tensor of A matrices, of dimension `design.shape[0]` x `self.npars` x `self.npars`

    The output is a tensor of shape `design.shape[0]`
    """
    temp = self.estimate_expected_FIM(design)
    temp = torch.matmul(A, temp)
    temp = torch.matmul(temp, A.transpose(1,2))
    return -temp.diagonal(dim1=1, dim2=2).sum(dim=1)

  def initialise_J_estimation(self):
    """Initialise variables needed for estimating J"""
    self.thetas_for_J = self.prior.sample((1000,))

  def estimateJ(self, design, adv=True):
    """Return an estimate of J objective

    `design` is a tensor of designs, batched along dimension 0
    `adv` - whether to calculate J_ADV (if True) or J_FIG (if False)

    The output is a tensor of shape `design.shape[0]`
    """
    temp = self.estimate_FIM(self.thetas_for_J, design)
    temp = torch.mean(temp, dim=1)
    if adv:
      temp = torch.det(temp)
    else:
      temp = temp.diagonal(dim1=1,dim2=2).sum(dim=1)
    return temp

  def eta_dim(self):
    """Number of entries required in eta vector

    i.e. count of non-zero entries in `npars` x `npars` triangular matrix minus 1
    """
    return self.npars * (self.npars + 1) // 2 - 1

  def makeA(self, x):
    """Convert tensor `x` to tensor of `A` matrices

    The final dimension of `x` should be `npars(npars+1)/2 - 1`

    Returns `A`, a tensor with shape `x.shape[:-1] + (npars, npars)`
    consisting of lower triangular matrices with positive diagonal
    and bottom right element set to give determinant 1.
    The elements of `x` fill the remainder of the matrix in an
    anticlockwise spiral.
    See `math.fill_triangular` from tensorflow probability for more details:
    the code here uses the same mathematical approach.
    """
    z = torch.zeros(x.shape[:-1] + (1,), dtype=x.dtype)
    x = torch.cat([z,x], dim=-1)
    x_list = [x[..., self.npars:], x.flip(-1)]
    A = torch.reshape(torch.cat(x_list, dim=-1), x.shape[:-1] + (self.npars, self.npars))
    A = torch.tril(A)
    d = A.diagonal(dim1=-2, dim2=-1)
    A[...,-1,-1] = -d.sum(-1)
    d.exp_()
    return A

class AdvOpt:
  """Optimises a Bayesian experimental design problem"""

  def __init__(self, fim,
               optimizers, schedulers,
               init_design_raw, init_A_raw,
               make_design=lambda x:x,
               penalty = lambda x:0., # Not batched but works via broadcasting
               report_every=500, text_progress=True,
               track_J=None):
    """
    `fim` - object of `FIM` class
    `optimizers` - dict of optimizers for `experimenter` and `adversary`
    (with no params)
    `schedulers` - dict of schedulers for `experimenter` and `adversary`
    `init_design_raw` - batched initial values for `design_raw`
    `init_A_raw` - batched initial values for variables controlling A matrix
    `make_design` - function mapping a batched `design_raw` tensor to designs
    `penalty` - function mapping a batched `design` tensor to penalties
    `report_every` - how often to record/report progress
    `text_progress` - report text summaries of progress if `True`
    `track_J` - can be "ADV", "FIG" or None
    """
    self.start_time = time()
    self.report_every = report_every
    self.track_J = track_J
    self.text_progress = text_progress
    if text_progress == True:
      np.set_printoptions(precision=3)

    self.fim = fim
    self.make_design = make_design
    self.penalty = penalty

    self.design_raw = torch.tensor(init_design_raw, dtype=torch.float32,
                                   requires_grad=True)
    eta_dim = fim.eta_dim()
    if init_A_raw.shape[1] != eta_dim:
      raise ValueError("Last init_A_raw dimension wrong, should be " + str(eta_dim))
    self.A_raw = torch.tensor(init_A_raw, dtype=torch.float32,
                              requires_grad=True)
    self.optimizers = optimizers
    self.schedulers = schedulers
    self.optimizers['experimenter'].add_param_group({'params': self.design_raw})
    self.optimizers['adversary'].add_param_group({'params': self.A_raw})

    self.output = {'iterations':[], 'time':[], 'design':[], 'A':[],
                   'objectiveK':[]}
    if track_J is not None:
      self.fim.initialise_J_estimation()
      self.output['objectiveJ'] = []
    self.iteration = 0

  def optimize(self, iterations=10**5):
    """Find optimal design numerically"""
    for i in range(iterations):
      for opt in self.optimizers.values():
        opt.zero_grad()
      A = self.fim.makeA(self.A_raw)
      design = self.make_design(self.design_raw)
      objective = self.fim.estimateK(design, A) + self.penalty(design)
      objective.sum().backward()
      self.A_raw.grad *= -1. # Do ascent instead of descent
      self.optimizers['experimenter'].step()
      self.schedulers['experimenter'].step()
      self.optimizers['adversary'].step()
      self.schedulers['adversary'].step()
      self.iteration += 1
      if (self.iteration % self.report_every) == 0:
        elapsed = time() - self.start_time
        self.output['time'].append(elapsed)
        self.output['iterations'].append(self.iteration)
        design_np = design.detach().numpy().copy()
        A_np = A.detach().numpy().copy()
        self.output['design'].append(design_np)
        self.output['A'].append(A_np)
        self.output['objectiveK'].append(objective.detach().numpy().copy())
        if self.track_J == "ADV":
          self.output['objectiveJ'].append(self.fim.estimateJ(design, adv=True).detach().numpy().copy())
        elif self.track_J == "FIG":
          self.output['objectiveJ'].append(self.fim.estimateJ(design, adv=False).detach().numpy().copy())

        if self.text_progress:
          print("Iteration {:d}, time (mins) {:.1f}, K objective {:.2f}".\
                format(i+1, elapsed/60, float(objective[0])))
          print("Design:\n", design_np[0, ...])
          print("A matrix:\n", A_np[0, ...])
          print("Learning rate: experimenter {:.6f} adversary {:.6f}\n".\
                format(self.optimizers['experimenter'].param_groups[0]['lr'],
                       self.optimizers['adversary'].param_groups[0]['lr']))

  def pointExchange(self, maxits=10, adv=True):
    """Point exchange optimisation to fine tune design

    This requires `self.design` to have dimensions `batches` x `design points` x ...

    `maxits` is maximum number of iterations to use
    `adv` should be True for the J_ADV objective or False for J_FIG objective
    """
    designs = self.make_design(self.design_raw).detach().numpy().copy()
    for rep in range(designs.shape[0]):
      design = np.expand_dims(designs[rep, ...], 0)
      bestJ = self.fim.estimateJ(torch.tensor(design), adv)[0].detach().numpy()
      best_new_design = design.copy()
      for outer_counter in range(maxits):
        improvement = False
        for i in range(design.shape[1]):
          for j in range(design.shape[1]):
            if i==j:
              continue
            proposed_design = design.copy()
            proposed_design[0,i,...] = design[0,j,...]
            newJ = self.fim.estimateJ(torch.tensor(proposed_design), adv)[0].detach().numpy()
            if (newJ > bestJ):
              bestJ = newJ
              best_new_design = proposed_design
              improvement = True
          ## End j for loop
        ## End i for loop
        design = best_new_design
        if self.text_progress==True:
          print("Replicate {:d}".format(rep+1))
          print("Point exchange iteration {:d}".format(outer_counter+1))
          print("Design:\n", np.sort(proposed_design[0,...]))
          print("Objective J estimate: {:.3f}\n".format(bestJ))
        if improvement == False:
          break
      ## End outer_counter for loop
      designs[rep, ...] = design
    ## End rep for loop
    return designs

  def stacked_output(self):
    """Return optimisation output stacked into arrays"""
    ##This always puts iteration as first dimension
    ##Not intuitive for arrays which also have a batch dimension
    stacked_output = {}
    for key in self.output.keys():
      stacked_output[key] = np.stack(self.output[key])

    return stacked_output

max_stars_repo_name: dennisprangle/AdversarialDesignCode

max_issues_repo_name: dennisprangle/AdversarialDesignCode

max_forks_repo_name: dennisprangle/AdversarialDesignCode


text:
"""
    THIS CODE IS UNDER THE BSD 2-Clause LICENSE. YOU CAN FIND THE COMPLETE
                        FILE AT THE SOURCE DIRECTORY.

    Copyright (C) 2017 V.-M. Katsageorgiou - All rights reserved

    @author : vasiliki.katsageorgiou@gmail.com

                        Publication:
    A Novel Unsupervised Analysis of Electrophysiological
        Signals Reveals New Sleep Sub-stages in Mice

*****************************************************************************

Script written for analysing the per strain daily profiles in the inferred
(from a trained RBM model) latent states for the study of sleep in mice:

here we test the homostatic response of the observed latent states in the
recovery mode (i.e. after sleep deprivation), comparing it to the baseline.

*********************************  OUTPUT    *********************************

Output : a folder named as "homostaticResponse" including the per latent
         state homostatic response profiles of the diferent mouse groups.

*****************************************************************************

<vkatsageorgiou@vassia-PC>
"""

from __future__ import division

import sys
import os
#from pylab import *
import numpy as np
import math
from numpy.random import RandomState
from numpy.polynomial import Polynomial
from scipy.io import loadmat, savemat
from configparser import *
import datetime
import PIL.Image

import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
from matplotlib.collections import LineCollection
from matplotlib import use
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d import proj3d
from matplotlib.ticker import MaxNLocator

from scipy.special import expit
from scipy.optimize import curve_fit
from scipy.interpolate import interp1d, UnivariateSpline
from scipy import stats
import _pickle as cPickle
import shutil

class homostaticResponse(object):
    '''
    Object for analyzing the daily profiles of the observed latent states,
    testing the response of the latent states in the Recovery mode
    (i.e. after sleep deprivation).
    '''

    def __init__(self, refDir, epochID, polydeg):
        # directory containing all the configuration files for the experiment
        self.refDir = refDir
        # id of the epoch to be analysed
        self.epochID = int(epochID)
        # degree of polynomial
        self.polydeg = int(polydeg)

        np.random.seed(124)
        self.prng =  RandomState(123)

    def loadData(self):
        '''
        Method loading the visible data.
        '''
        visData = 'visData.npz'
        dataFile = np.load(self.refDir + visData)
        self.d = dataFile['data']
        self.epochIDs = dataFile['obsKeys'].astype(int)
        self.epochTime = dataFile['epochTime']#.astype(int)

        """
        Move in final analysed model folder
        """
        os.chdir(self.refDir + 'analysis/epoch%d' %self.epochID)

        """
        Load analysed data: array with epochIDs & the corresponding
        latent-states' IDs.
        """
        self.obsKeys = np.load('obsKeys.npz')['obsKeys'].astype(int)

        """
        Associate each epoch with the strain ID
        """
        self.obsKeys = np.insert(self.obsKeys, self.obsKeys.shape[1], 0, axis=1)
        for i in range(self.obsKeys.shape[0]):
            self.obsKeys[i, self.obsKeys.shape[1]-1] = int(str(self.obsKeys[i, self.obsKeys.shape[1]-2])[0])

        self.lstates = np.unique(self.obsKeys[:, 1])
        self.subjIDs = np.unique(self.obsKeys[:, self.obsKeys.shape[1]-2])
        self.strainIDs = np.unique(self.obsKeys[:, self.obsKeys.shape[1]-1])

        #self.epochTime = self.epochTime[:self.obsKeys.shape[0],:]

        self.utc_transform()

    def utc_transform(self):
        '''
        Method for transforming the serial date-time to integer.
        '''

        # add the column with the mode id: 1-baseline, 2-recovery
        self.obsKeys = np.hstack((self.obsKeys, self.epochTime[:, self.epochTime.shape[1]-2].reshape(-1,1))).astype(int)

        self.epochTime = np.insert(self.epochTime, self.epochTime.shape[1], 0, axis=1)

        """
        Iterate through epochs:
        """
        for i in range(len(self.epochTime)):
            self.epochTime[i, self.epochTime.shape[1]-1] = int(datetime.datetime.utcfromtimestamp((self.epochTime[i, self.epochTime.shape[1]-2]-25569) * 86400.0).strftime("%H"))

        """
        Create unique array: epochsIDs - lstatesIDs - subjects - mode - dayTime
        """
        self.obsKeys = np.hstack((self.obsKeys, self.epochTime[:, self.epochTime.shape[1]-1].reshape(-1,1))).astype(int)
        savemat('obsKeysTime.mat', mdict={'obsKeys':self.obsKeys})

    def combinedGroupHistograms(self):
        '''
        Method for visualizing the distribution of each latent state in time
        per group (strain).
        '''
        if not os.path.isdir('homostaticResponse'):
            os.makedirs('homostaticResponse')

        self.baselineHours = [7, 7, 7, 7, 7, 7, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0, 1, 2, 3, 4, 5, 6, 6, 6, 6, 6, 6, 6]
        self.recoveryHours = [13, 13, 13, 13, 13, 13, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0, 1, 2, 3, 4, 5, 6, 6, 6, 6, 6, 6, 6]
        self.baselineInd = np.arange(len(self.baselineHours))
        self.recoveryInd = np.arange(len(self.recoveryHours))

        plt.style.use('bmh')
        palete = plt.rcParams['axes.color_cycle']
        colors = ['#b2182b', palete[9], palete[8], '#1a1a1a', '#004529']

        width = 0.4

        """ Iterate through latent states """
        for self.lstate in self.lstates:

            idx = np.where(self.obsKeys[:,1]==self.lstate)[0]

            if len(idx) >= 100:

                fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(60,35))

                ax1.spines['top'].set_visible(False)
                ax1.spines['right'].set_visible(False)

                ax1.spines['bottom'].set_linewidth(20)
                ax1.spines['bottom'].set_color('k')
                ax1.spines['left'].set_linewidth(20)
                ax1.spines['left'].set_color('k')

                ax1.yaxis.set_ticks_position('left')
                ax1.xaxis.set_ticks_position('bottom')
                ax1.grid(False)

                ax2.spines['top'].set_visible(False)
                ax2.spines['right'].set_visible(False)

                ax2.spines['bottom'].set_linewidth(20)
                ax2.spines['bottom'].set_color('k')
                ax2.spines['left'].set_linewidth(20)
                ax2.spines['left'].set_color('k')
                ax2.yaxis.set_ticks_position('left')
                ax2.xaxis.set_ticks_position('bottom')
                ax2.grid(False)

                currObsKeys = self.obsKeys[idx,:]
                epochsW = round((len(np.where((currObsKeys[:, currObsKeys.shape[1]-5]==1))[0]) / len(currObsKeys)), 3)
                epochsNR = round((len(np.where((currObsKeys[:, currObsKeys.shape[1]-5]==2))[0]) / len(currObsKeys)), 3)
                epochsR = round((len(np.where((currObsKeys[:, currObsKeys.shape[1]-5]==3))[0]) / len(currObsKeys)), 3)

                plt.text(0.01, 0.98, 'LS: ' + str(self.lstate), transform=ax1.transAxes, horizontalalignment='left', verticalalignment='top', color='#252525', fontsize=150, fontweight='bold')

                behaviorDistribution = [epochsW, epochsNR, epochsR]
                max_idx = behaviorDistribution.index(max(behaviorDistribution))

                if max_idx == 0:
                    plt.text(0.5, 1., 'W: ' + str(epochsW*100) + '%', transform=ax1.transAxes, horizontalalignment='left', verticalalignment='top', color='#b2182b', fontsize=150, fontweight='bold')
                    plt.text(0.5, 0.93, 'NR: ' + str(epochsNR*100) + '%', transform=ax1.transAxes, horizontalalignment='left', verticalalignment='top', color='#1a1a1a', fontsize=150)
                    plt.text(0.5, 0.85, 'R: ' + str(epochsR*100) + '%', transform=ax1.transAxes, horizontalalignment='left', verticalalignment='top', color='#1a1a1a', fontsize=150)
                elif max_idx == 1:
                    plt.text(0.5, 1., 'W: ' + str(epochsW*100) + '%', transform=ax1.transAxes, horizontalalignment='left', verticalalignment='top', color='#1a1a1a', fontsize=150)
                    plt.text(0.5, 0.93, 'NR: ' + str(epochsNR*100) + '%', transform=ax1.transAxes, horizontalalignment='left', verticalalignment='top', color='#b2182b', fontsize=150, fontweight='bold')
                    plt.text(0.5, 0.85, 'R: ' + str(epochsR*100) + '%', transform=ax1.transAxes, horizontalalignment='left', verticalalignment='top', color='#1a1a1a', fontsize=150)
                else:
                    plt.text(0.5, 1., 'W: ' + str(epochsW*100) + '%', transform=ax1.transAxes, horizontalalignment='left', verticalalignment='top', color='#1a1a1a', fontsize=150)
                    plt.text(0.5, 0.93, 'NR: ' + str(epochsNR*100) + '%', transform=ax1.transAxes, horizontalalignment='left', verticalalignment='top', color='#1a1a1a', fontsize=150)
                    plt.text(0.5, 0.85, 'R: ' + str(epochsR*100) + '%', transform=ax1.transAxes, horizontalalignment='left', verticalalignment='top', color='#b2182b', fontsize=150, fontweight='bold')

                """ Iterate through mouse strains """
                ci = 0
                maxStrainID = max(self.strainIDs)
                max_graph = 0.

                for self.strain in self.strainIDs:
                    idxStr = np.where(currObsKeys[:, currObsKeys.shape[1]-3]==self.strain)[0]

                    if len(idxStr) != 0:
                        strainObsKeys = currObsKeys[idxStr, :]

                        # spit into modes
                        baseline = strainObsKeys[np.where( strainObsKeys[:, strainObsKeys.shape[1]-2] == 1 )[0], :]
                        recovery = strainObsKeys[np.where( strainObsKeys[:, strainObsKeys.shape[1]-2] == 2 )[0], :]

                        if (len(baseline) > 0. and len(recovery) > 0.):

                            count_baseline, p_baseline = self.computeDistribution(baseline, 'baseline')
                            count_recovery, p_recovery = self.computeDistribution(recovery, 'recovery')

                            idx_b_start = np.where( p_baseline.linspace()[0] >= 6. )[0][0]
                            idx_b_end = np.where( p_baseline.linspace()[0] >= len(model.baselineHours)-6 )[0][0]

                            idx_r_start = np.where( p_recovery.linspace()[0] >= 6. )[0][0]
                            idx_r_end = np.where( p_recovery.linspace()[0] >= len(model.recoveryHours)-6 )[0][0]

                            if max(max(p_recovery.linspace()[1][idx_r_start:idx_r_end]), max(p_baseline.linspace()[1][idx_b_start:idx_b_end])) > max_graph:
                                max_graph = max(max(p_recovery.linspace()[1][idx_r_start:idx_r_end]), max(p_baseline.linspace()[1][idx_b_start:idx_b_end]))

                            """ Compute curve difference """

                            if self.strain<maxStrainID:

                                l = '$\mathregular{Zfhx3^{Sci{/}{+}}}$'

                                ax2.plot(*p_baseline.linspace(), c=colors[ci+1], linewidth=20.0, label='baseline')
                                ax2.plot(p_recovery.linspace()[0][idx_r_start:idx_r_end] + 6., p_recovery.linspace()[1][idx_r_start:idx_r_end], c=colors[ci], linewidth=20.0, label='recovery')

                                t2 = ax2.set_title(l, fontsize=150, fontweight='bold', loc='left')
                                t2.set_position([0.0, 1.05])
                                ax2.axvspan(0., p_recovery.linspace()[0][idx_r_start]+6, facecolor='#b2182b', alpha=0.3)
                                ax2.text(1.21, .08, 'total: ' + str(len(strainObsKeys)) + ' epochs', transform=ax1.transAxes, horizontalalignment='left', verticalalignment='top', color='#252525', fontsize=150)
                            else:
                                l = '$\mathregular{Zfhx3^{{+}{/}{+}}}$'

                                ax1.plot(*p_baseline.linspace(), c=colors[ci+1], linewidth=20.0, label='baseline')
                                ax1.plot(p_recovery.linspace()[0][idx_r_start:idx_r_end] + 6., p_recovery.linspace()[1][idx_r_start:idx_r_end], c=colors[ci], linewidth=20.0, label='recovery')

                                t1 = ax1.set_title(l, fontsize=150, fontweight='bold', loc='left')
                                t1.set_position([0.0, 1.05])
                                ax1.axvspan(0., p_recovery.linspace()[0][idx_r_start]+6, facecolor='#b2182b', alpha=0.3)
                                ax1.text(0.01, .08, 'total: ' + str(len(strainObsKeys)) + ' epochs', transform=ax1.transAxes, horizontalalignment='left', verticalalignment='top', color='#252525', fontsize=150)

                        else:
                            continue

                    else:
                        continue

                ax1.set_ylabel('Number of epochs', fontweight='bold', fontsize=150, labelpad=30)

                xTickMarks = []
                for j in np.arange(-6, 31):
                    if j in [0, 6, 12]:
                        xTickMarks.append('%s' %str(j))
                    elif j == 23:
                        xTickMarks.append('%s' %str(24))
                    else:
                        xTickMarks.append('')

                ax1.set_xticks(self.baselineInd+width/2)
                xtickNames = ax1.set_xticklabels(xTickMarks, fontweight='bold', fontsize=130)

                ax1.xaxis.set_ticks_position('none')
                ax1.yaxis.set_ticks_position('none')

                ax1.set_ylim([0, int(math.ceil(max_graph / 10.0)) * 10])

                yint = []
                locs, labels = plt.yticks()
                for each in locs:
                    yint.append(int(each))
                plt.yticks(yint)

                ax1.set_yticklabels(ax1.get_yticks(), fontweight='bold', fontsize=130)

                ax1.set_xlim([6, len(self.baselineHours)-6])
                ax1.set_xlabel('ZT', fontweight='bold', fontsize=150, labelpad=40)

                ax2.xaxis.set_ticks_position('none')
                ax2.yaxis.set_ticks_position('none')

                ax2.set_xticks(self.baselineInd+width/2)
                xtickNames = ax2.set_xticklabels(xTickMarks, fontweight='bold', fontsize=130)

                ax2.set_xlim([6, len(self.baselineHours)-6])
                ax2.set_xlabel('ZT', fontweight='bold', fontsize=150, labelpad=40)

                legend_properties = {'weight':'bold', 'size':130}
                legend = plt.legend(bbox_to_anchor=(.5, 1.), loc='upper left', borderaxespad=0., prop=legend_properties)
                frame = legend.get_frame().set_alpha(0)

                fname = 'lstate%d.tiff' %self.lstate
                fname = os.path.join('./homostaticResponse/', fname)
                fig.savefig(fname, format='tiff', transparent=True, dpi=100)
                #plt.show()
                plt.close(fig)

    def computeDistribution(self, d, mode):
        """
        Function computing the discrete distribution (histogram) of each
        latent state over the 24h in one hour bin
        """

        if mode == 'baseline' :

            """    Iterate through hours and compute histogram    """

            count = np.zeros((len(np.arange(24)), 2), dtype=float)
            for t in np.arange(24):
                count[t, 0] = t
                count[t, 1] = len( np.where(d[:, d.shape[1]-1]==t)[0] )

            idx_baseHours = []
            for i in self.baselineHours:
                i = float(i)
                idx_baseHours.append( np.where(count[:, 0] == i)[0][0] )
            count = count[idx_baseHours, :]

            """    Fit curve on histogram : Polynomial fit    """

            p = Polynomial.fit(self.baselineInd, count[self.baselineHours, 1], self.polydeg)
            pConv = p.convert(domain=[-1, 1])

        else:
            """    Iterate through hours and compute histogram """
            count = np.zeros((len(np.arange(18)), 2), dtype=float)
            i = 0
            for t in np.arange(24):
                if t in self.recoveryHours:
                    count[i, 0] = t
                    count[i, 1] = len( np.where(d[:, d.shape[1]-1]==t)[0] )
                    i += 1

            idx_recHours = []
            for i in self.recoveryHours:
                i = float(i)
                idx_recHours.append( np.where(count[:, 0] == i)[0][0] )
            count = count[idx_recHours, :]

            """    Fit curve on histogram : Polynomial fit    """
            p = Polynomial.fit(self.recoveryInd, count[:, 1], self.polydeg)
            pConv = p.convert(domain=[-1, 1])

        return count, p

if __name__ == "__main__":

    import argparse

    parser = argparse.ArgumentParser()
    parser.add_argument('-refDir', help='Experiment path')
    parser.add_argument('-epoch', help='Epoch analysed')
    parser.add_argument('-deg', help='Degrees of polynomial')
    args = parser.parse_args()

    print ('Initialization..')
    model = homostaticResponse(args.refDir, args.epoch, args.deg)

    print ('Loading data..')
    model.loadData()

    print ('Computing circadian profiles..')
    model.combinedGroupHistograms()

    with open(args.refDir + 'doneHomoResp', 'w') as doneFile:
        doneFile.write(datetime.datetime.strftime(datetime.datetime.now(), '%d/%m/%Y %H:%M:%S'))

ext: py

max_stars_repo_name: davidmam/EEG-sleep-analysis

max_issues_repo_name: davidmam/EEG-sleep-analysis

max_forks_repo_name: davidmam/EEG-sleep-analysis

reason: import numpy,from numpy,from scipy


text:
import os

from xcs_soxs.background.spectra import BackgroundSpectrum, \
    ConvolvedBackgroundSpectrum
from xcs_soxs.background.events import make_diffuse_background
from xcs_soxs.utils import soxs_files_path, parse_prng, mylog
import numpy as np

# X-ray foreground from Hickox & Markevitch 2007
# (http://adsabs.harvard.edu/abs/2007ApJ...661L.117H)
hm_bkgnd_file = os.path.join(soxs_files_path, "hm_cxb_bkgnd.h5")
hm_astro_bkgnd = BackgroundSpectrum.from_file(hm_bkgnd_file)

def make_foreground(event_params, arf, rmf, prng=None):
    import pyregion._region_filter as rfilter

    prng = parse_prng(prng)

    conv_frgnd_spec = ConvolvedBackgroundSpectrum(hm_astro_bkgnd, arf)

    energy = conv_frgnd_spec.generate_energies(event_params["exposure_time"],
                                               event_params["fov"], prng=prng,
                                               quiet=True).value

    prng = parse_prng(prng)

    bkg_events = {}

    n_events = energy.size

    nx = event_params["num_pixels"]
    bkg_events["detx"] = prng.uniform(low=-0.5*nx, high=0.5*nx, size=n_events)
    bkg_events["dety"] = prng.uniform(low=-0.5*nx, high=0.5*nx, size=n_events)
    bkg_events["energy"] = energy

    if event_params["chips"] is None:
        bkg_events["chip_id"] = np.zeros(n_events, dtype='int')
    else:
        bkg_events["chip_id"] = -np.ones(n_events, dtype='int')
        for i, chip in enumerate(event_params["chips"]):
            thisc = np.ones(n_events, dtype='bool')
            rtype = chip[0]
            args = chip[1:]
            r = getattr(rfilter, rtype)(*args)
            inside = r.inside(bkg_events["detx"], bkg_events["dety"])
            thisc = np.logical_and(thisc, inside)
            bkg_events["chip_id"][thisc] = i

    keep = bkg_events["chip_id"] > -1

    if keep.sum() == 0:
        raise RuntimeError("No astrophysical foreground events were detected!!!")
    else:
        mylog.info("Making %d events from the astrophysical foreground." % keep.sum())

    for key in bkg_events:
        bkg_events[key] = bkg_events[key][keep]

    return make_diffuse_background(bkg_events, event_params, rmf, prng=prng)

reason: import numpy


text:
# Copyright (c) 2021, NVIDIA CORPORATION.  All rights reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import os
import pickle as pkl
import random
from argparse import ArgumentParser

import h5py
import numpy as np
import torch
from omegaconf import DictConfig, OmegaConf
from sklearn.decomposition import PCA
from tqdm import tqdm

from nemo.collections.nlp.models import EntityLinkingModel
from nemo.utils import logging

try:
    import faiss
except ModuleNotFoundError:
    logging.warning("Faiss is required for building the index. Please install faiss-gpu")

device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

def build_index(cfg: DictConfig, model: object):
    """
    Builds faiss index from index dataset specified in the config.

    Args:
        cfg (DictConfig): Config file specifying index parameters
        model (object): Encoder model
    """

    # Get index dataset embeddings
    # PCA model exists and index embeddings have already been PCAed, no need to re-extract/PCA them
    if cfg.apply_pca and os.path.isfile(cfg.pca.pca_save_name) and os.path.isfile(cfg.pca_embeddings_save_name):
        logging.info("Loading reduced dimensionality embeddings")
        embeddings = h5py.File(cfg.pca_embeddings_save_name, "r")
        embeddings = embeddings[cfg.index_ds.name][:]

    elif os.path.isfile(cfg.embedding_save_name):
        logging.info("Loading previously extracted index dataset embeddings")
        embeddings = h5py.File(cfg.embedding_save_name, "r")
        embeddings = embeddings[cfg.index_ds.name][:]

    else:
        logging.info("Encoding index dataset, this may take a while")
        index_dataloader = model.setup_dataloader(cfg.index_ds, is_index_data=True)
        embeddings, concept_ids = get_index_embeddings(cfg, index_dataloader, model)

    # Create pca model to reduce dimensionality of index dataset and decrease memory footprint
    if cfg.apply_pca:

        # Need to train PCA model and apply PCA transformation with newly trained model
        if not os.path.isfile(cfg.pca.pca_save_name):
            logging.info("Fitting PCA model for embedding dimensionality reduction")
            pca_train_set = random.sample(list(embeddings), k=int(len(embeddings) * cfg.pca.sample_fraction))
            pca = PCA(n_components=cfg.pca.output_dim)
            pca.fit(pca_train_set)
            pkl.dump(pca, open(cfg.pca.pca_save_name, "wb"))
            embeddings = reduce_embedding_dim(pca, embeddings, cfg)

        # PCA model already trained, just need to reduce dimensionality of all embeddings
        elif not os.path.isfile(cfg.pca_embeddings_save_name):
            pca = pkl.load(open(cfg.pca.pca_save_name, "rb"))
            embeddings = reduce_embedding_dim(pca, embeddings, cfg)

    # Build faiss index from embeddings
    logging.info(f"Training index with embedding dim size {cfg.dims} using {faiss.get_num_gpus()} gpus")
    quantizer = faiss.IndexFlatL2(cfg.dims)
    index = faiss.IndexIVFFlat(quantizer, cfg.dims, cfg.nlist)
    index = faiss.index_cpu_to_all_gpus(index)
    index.train(embeddings)

    logging.info("Adding dataset embeddings to index")
    for i in tqdm(range(0, embeddings.shape[0], cfg.index_batch_size)):
        index.add(embeddings[i : i + cfg.index_batch_size])

    logging.info("Saving index")
    faiss.write_index(faiss.index_gpu_to_cpu(index), cfg.index_save_name)
    logging.info("Index built and saved")

def reduce_embedding_dim(pca, embeddings, cfg):
    """Apply PCA transformation to index dataset embeddings"""

    logging.info("Applying PCA transformation to entire index dataset")
    embeddings = np.array(pca.transform(embeddings), dtype=np.float32)
    emb_file = h5py.File(cfg.pca_embeddings_save_name, "w")
    emb_file.create_dataset(cfg.index_ds.name, data=embeddings)
    emb_file.close()

    return embeddings

def get_index_embeddings(cfg: DictConfig, dataloader: object, model: object):
    """Use entity linking encoder to get embeddings for full index dataset"""
    embeddings = []
    concept_ids = []

    with torch.no_grad():
        for batch in tqdm(dataloader):
            input_ids, token_type_ids, input_mask, batch_concept_ids = batch
            input_ids = input_ids.to(device)
            token_type_ids = token_type_ids.to(device)
            input_mask = input_mask.to(device)
            batch_embeddings = model.forward(
                input_ids=input_ids, token_type_ids=token_type_ids, attention_mask=input_mask
            )

            embeddings.extend(batch_embeddings.detach().cpu().numpy())
            concept_ids.extend(batch_concept_ids.numpy())

    emb_file = h5py.File(cfg.embedding_save_name, "w")
    emb_file.create_dataset(cfg.index_ds.name, data=embeddings)
    emb_file.close()

    pkl.dump(concept_ids, open(cfg.concept_id_save_name, "wb"))

    return embeddings, concept_ids

def load_model(cfg: DictConfig, restore: bool):
    """
    Loads encoder model.

    Args:
        cfg: Config file specifying model parameters
        restore: Whether to restore model weights trained
                 by the user. Otherwise will load weights
                 used before self alignment pretraining.
    """

    if restore:
        model = EntityLinkingModel.restore_from(cfg.nemo_path)
    else:
        cfg.train_ds = None
        cfg.validation_ds = None
        cfg.test_ds = None
        model = EntityLinkingModel(cfg)

    model = model.to(device)

    return model

def main(cfg: DictConfig, restore: bool):
    """
    Builds new index if one hasn't been built yet.

    Args:
        cfg: Config file specifying index parameters
        restore: Whether to restore model weights trained
                 by the user. Otherwise will load weights
                 used before self alignment pretraining.
    """

    logging.info("Loading entity linking encoder model")
    model = load_model(cfg.model, restore)

    if not os.path.isfile(cfg.index.index_save_name) or (
        cfg.apply_pca and not os.path.isfile(cfg.index.pca.pca_save_name)
    ):
        logging.info("Building index")
        build_index(cfg.index, model)
    else:
        logging.info("Index and pca model (if required) already exists. Skipping build index step.")

    if not os.path.isfile(cfg.index.idx_to_id):
        logging.info("Mapping entity index postions to ids")
        map_idx_to_ids(cfg.index)
    else:
        logging.info("Map from concept index to id already exists. Skipping mapping step.")

if __name__ == '__main__':
    parser = ArgumentParser()
    parser.add_argument(
        "--restore", action="store_true", help="Whether to restore encoder model weights from nemo path"
    )
    parser.add_argument("--project_dir", required=False, type=str, default=".")
    parser.add_argument("--cfg", required=False, type=str, default="./conf/umls_medical_entity_linking_config.yaml")
    args = parser.parse_args()

    cfg = OmegaConf.load(args.cfg)
    cfg.project_dir = args.project_dir

    main(cfg, args.restore)

max_stars_repo_name: madhukarkm/NeMo

max_issues_repo_name: madhukarkm/NeMo

max_forks_repo_name: madhukarkm/NeMo


text:
# coding: utf-8

import os
import datetime
os.environ["OPENBLAS_NUM_THREADS"] = "1"

'''
Created on 27 марта 2016 г.

@author: keen
'''
import CO2_tree as co2
from multiprocessing import Pool
from functools import partial

from numpy import argmax
from numpy import multiply
from numpy import asarray

from joblib import Parallel, delayed
import os
from scipy.sparse.csr import csr_matrix

from scipy.special import kl_div

from numpy import load
from numpy import save
import uuid
import numpy

from scipy.spatial.distance import jensenshannon
from scipy.stats import entropy

#import cProfile

#def profile(func):
#    """Decorator for run function profile"""
#    def wrapper(*args, **kwargs):
#        profile_filename = func.__name__ + '.prof'
#        profiler = cProfile.Profile()
#        result = profiler.runcall(func, *args, **kwargs)
#        profiler.dump_stats(profile_filename)
#        return result
#    return wrapper

def fitter(uuids,forest,shapex,seed_):
    dataX = load('/dev/shm/' + uuids + 'DataX.npy',mmap_mode='r')
    indX = load('/dev/shm/' + uuids + "IndX.npy",mmap_mode='r')
    ptrX = load('/dev/shm/' + uuids + "PtrX.npy",mmap_mode='r')
    x = csr_matrix((dataX,indX,ptrX), shape=shapex,dtype=numpy.float32,copy=False)
    Y = load('/dev/shm/' + uuids + "DataY.npy",mmap_mode='r')

    tree = co2.CO2Tree(C=forest.C , kernel=forest.kernel,\
    tol=forest.tol, max_iter=forest.max_iter,max_deth = forest.max_deth,\
     min_samples_split = forest.min_samples_split,dual=forest.dual,\
    min_samples_leaf = forest.min_samples_leaf, seed = None,\
     sample_ratio = forest.sample_ratio, feature_ratio = forest.feature_ratio, \
     gamma=forest.gamma,intercept_scaling=forest.intercept_scaling,dropout_low=forest.dropout_low,dropout_high=forest.dropout_high,noise=forest.noise,cov_dr=forest.cov_dr, criteria = forest.criteria)
    tree.fit(x,Y, preprocess = False)
    return tree

def probber(uuids, shapex,tree,stat_only,use_weight = True,withY = False):
        dataX = load('/dev/shm/' + uuids + 'DataX.npy',mmap_mode='r')
        indX = load('/dev/shm/' + uuids + "IndX.npy",mmap_mode='r')
        ptrX = load('/dev/shm/' + uuids + "PtrX.npy",mmap_mode='r')

        if withY:
            Y = load('/dev/shm/' + uuids + "DataY.npy",mmap_mode='r')
        else:
            Y = None
        x = csr_matrix((dataX,indX,ptrX), shape=shapex,dtype=numpy.float32, copy=False)
        return tree.predict_proba(x,Y,preprocess = False,stat_only=stat_only,use_weight=use_weight)

def statter(tree):
    return tree.getWeights()

def jsd(norm,i):
    sum_ = entropy(norm[i]) * norm.shape[0]

    for j in range(norm.shape[0]):
        if i != j:
            sum_ += entropy(norm[i],norm[j]) + entropy(norm[j],norm[i]) + entropy(norm[j])
    return float(sum_) / norm.shape[0]

def weighter(tree,forest,norm,diff,min_nom,uniform):
    t = forest.trees[tree]
    t.estimateChunkWeights(norm,diff,min_nom,uniform)
    return t

class CO2_forest:
    def stat(self):
        return Parallel(n_jobs=self.n_jobs)(delayed(statter)(t) for t in self.trees)

    def sequential_fit(self,x,Y):
        self.trees = []
        forest = self
        for i in range(self.n_estimators):
            tree = co2.CO2Tree(C=forest.C , kernel=forest.kernel,\
            tol=forest.tol, max_iter=forest.max_iter,max_deth = forest.max_deth,\
            min_samples_split = forest.min_samples_split,dual=forest.dual,\
            min_samples_leaf = forest.min_samples_leaf, seed = None,\
            sample_ratio = forest.sample_ratio, feature_ratio = forest.feature_ratio, \
            gamma=forest.gamma,intercept_scaling=forest.intercept_scaling,dropout_low=forest.dropout_low,dropout_high=forest.dropout_high,noise=forest.noise,cov_dr=forest.cov_dr, criteria = forest.criteria)

            tree.fit(x,Y, preprocess = False)
            self.trees.append(tree)

    #@profile
    def fit(self,x,Y):
        uuids = str(uuid.uuid4())

        save('/dev/shm/'+ uuids + "DataX",x.data)
        save('/dev/shm/'+ uuids + "IndX",x.indices)
        save('/dev/shm/'+ uuids + "PtrX",x.indptr)
        save('/dev/shm/'+ uuids + "DataY",Y)

        #pool = Pool(self.n_jobs)
        self.trees = Parallel(n_jobs=self.n_jobs,backend="multiprocessing")(delayed(fitter)(uuids,self,x.shape,i) for i in range(self.n_estimators))

        #self.trees = pool.map(func,range(self.n_estimators))
        #pool.close()
        #pool.join()

        os.remove('/dev/shm/'+ uuids + "DataX.npy")
        os.remove('/dev/shm/'+ uuids + "IndX.npy")
        os.remove('/dev/shm/'+ uuids + "PtrX.npy")
        os.remove('/dev/shm/'+ uuids + "DataY.npy")

        if self.cov_dr > 0:

            counts = []
            for t in self.trees:
                counts.append(t.getCounts())

            norm = numpy.concatenate(counts)
            #norm = asarray(norm,float) / (norm.sum(axis=1).reshape(-1,1)) + 0.0001

            self.covs = numpy.zeros((norm.shape[1],norm.shape[1]))

            norm = norm.sum(axis=1)

            #means = norm.mean(axis=1)
            #for i in range(norm.shape[1]):
                #nm = norm[:,i].mean()
            #    for j in range(norm.shape[1]):
                    #self.covs[i,j] =  ((norm[:,i] - nm)*(norm[:,j] - norm[:,j].mean())).mean()
            #        cvs = (norm[:,i] - means[i])*(norm[:,j] - means[j])
                    #cvs[cvs < 0.] = 0.
            #        self.covs[i,j] =  cvs.mean()

            #print  (datetime.datetime.now())
            #chain_covs = Parallel(n_jobs=self.n_jobs)(delayed(jsd)(norm,i) for i in range(norm.shape[0]))
            #div_arr = 0.5*kl_div(norm,avg)
            #sm = div_arr[i].sum()*norm.shape[0] + div_arr[:i].sum() + div_arr[i+1:].sum()

            #print  (datetime.datetime.now())
            #chain_covs = numpy.asarray(chain_covs)
            #chain_covs = norm.dot(numpy.transpose(1. - norm)).sum(axis=1)

            min_cov =  norm.min()
            max_cov =  norm.max()
            #avg_cov = (max_cov + min_cov) / 2.0
            max_diff = max_cov - min_cov

            if not self.cov_dr == 0.:
                global_counter = 0
                for t in self.trees:
                    for lidx in t.leaves:
                        dst = float(norm[global_counter] - min_cov) / max_diff
                        #dst = 1. / (numpy.exp(-2* dst + 1) + 1)
                        #dst = chain_covs[global_counter]
                        t.nodes[lidx].chunk_weight = dst
                        global_counter += 1
            #self.trees = Parallel(n_jobs=self.n_jobs,backend="multiprocessing")(delayed( weighter)(i,self,norm, diff,min_nom,uniform)for i in range(self.n_estimators))

    def sequential_predict(self,x,use_weight=True):
        probas = []
        for tree in self.trees:
            probas.append(tree.predict_proba(x,None,preprocess = False,stat_only=False,use_weight=use_weight))

        proba =  multiply(asarray(probas).sum(axis=0), 1. / self.n_estimators)
        res =  argmax(proba, axis = 1)
        zr = res == 0
        res[zr] = 1
        return res

    def predict(self,x,Y=None,use_weight=True):
        if Y is not None:
            proba, cmp = self.predict_proba(x,Y,use_weight=use_weight)
        else:
            proba = self.predict_proba(x,Y,use_weight=use_weight)
        res =  argmax(proba, axis = 1)
        zr = res == 0
        res[zr] = 1
        if Y is not None:
            return res,cmp
        else:
            return res

    def predict_proba(self,x,Y=None,avg='macro',use_weight=True):
        uuids = str(uuid.uuid4())

        save('/dev/shm/'+ uuids + "DataX",x.data)
        save('/dev/shm/'+ uuids + "IndX",x.indices)
        save('/dev/shm/'+ uuids + "PtrX",x.indptr)

        if Y is not None:
            save('/dev/shm/'+ uuids + "DataY",Y)
        res = Parallel(n_jobs=self.n_jobs)(delayed(probber)(uuids,x.shape,t,False,use_weight,Y is not None) for t in self.trees)

        os.remove('/dev/shm/'+ uuids + "DataX.npy")
        os.remove('/dev/shm/'+ uuids + "IndX.npy")
        os.remove('/dev/shm/'+ uuids + "PtrX.npy")

        rr = []
        rs = []
        if Y is not None:
            os.remove('/dev/shm/'+ uuids + "DataY.npy")
            for r in res:
                rr.append(r[0])
                rs += r[1]
            res = rr
        if Y is not None:
            if avg == 'macro':
                return multiply(asarray(res).sum(axis=0), 1. / self.n_estimators),rs
            else:
                return asarray(res), rs
        else:
            if avg == 'macro':
                return multiply(asarray(res).sum(axis=0), 1. / self.n_estimators)
            else:
                return asarray(res)

    def __init__(self,C, kernel = 'linear', max_deth = None, tol = 0.001, min_samples_split = 2, \
                 dual=True,max_iter=1000000,
                 min_samples_leaf = 1, n_jobs=1, n_estimators = 10,sample_ratio = 1.0,feature_ratio=1.0,gamma=1000.,intercept_scaling=1.,dropout_low=0.,dropout_high=1.0,noise=0.,cov_dr=0., criteria='gini'):
        self.criteria = criteria
        self.C = C
        self.min_samples_split = min_samples_split
        self.min_samples_leaf = min_samples_leaf
        self.kernel = kernel
        self.max_deth = max_deth
        self.n_estimators = n_estimators
        self.n_jobs = n_jobs
        self.trees = []
        self.tol = tol
        self.sample_ratio = sample_ratio
        self.gamma = gamma
        self.dual = dual
        self.max_iter = max_iter
        self.feature_ratio=feature_ratio
        self.intercept_scaling = intercept_scaling
        self.dropout_low = dropout_low
        self.dropout_high = dropout_high
        self.noise = noise
        os.environ["OPENBLAS_NUM_THREADS"] = "1"
        self.cov_dr = cov_dr


text:
import os
import os.path
import sys
from glob import glob
from pathlib import Path
from typing import Optional
import random

import adaptive
import matplotlib
import matplotlib.cm
import matplotlib.colors as colors
import matplotlib.font_manager as fm
import matplotlib.patheffects as patheffects
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
import numpy as np
from matplotlib import pyplot as plt

def get_cmap(cmap, min_clip=0.0, max_clip=1.0, exp=1.0):
    fcmap = getattr(plt.cm, cmap)
    return colors.LinearSegmentedColormap.from_list(
        "my_colormap", fcmap(np.linspace(min_clip, max_clip, 256) ** exp),
    )

def learner_till(till, learner, data):
    new_learner = adaptive.Learner2D(None, bounds=learner.bounds)
    new_learner.data = {k: v for k, v in data[:till]}
    for x, y in learner._bounds_points:
        # always include the bounds
        new_learner.tell((x, y), learner.data[x, y])
    return new_learner

def plot_tri(learner, ax, xy_size):
    ip = learner.ip()
    tri = ip.tri
    xs, ys = tri.points.T
    x_size, y_size = xy_size
    triang = mtri.Triangulation(x_size * xs, y_size * ys, triangles=tri.vertices)
    return ax.triplot(triang, c="k", lw=0.3, alpha=1, zorder=2), (ip.values, triang)

def to_gradient(data, horizontal, cmap, spread=20, mid=0.5):
    n, m = data.shape if horizontal else data.shape[::-1]
    x = np.linspace(1, 0, n)
    x = 1 / (np.exp((x - mid) * spread) + 1)  # Fermi-Dirac like
    gradient = x.reshape(1, -1).repeat(m, 0)
    if not horizontal:
        gradient = gradient.T
    gradient_rgb = cmap(data)
    gradient_rgb[:, :, -1] = gradient
    return gradient_rgb

def get_new_artists(npoints_tri, learner, data, ax, xy_size, npoints_interp, cmap):
    new_learner = learner_till(npoints_tri, learner, data)
    (line1, line2), (zs, triang) = plot_tri(new_learner, ax, xy_size)
    data = learner.interpolated_on_grid(npoints_interp)[
        -1
    ]  # This uses the original learner!
    x_size, y_size = xy_size
    im = ax.imshow(
        to_gradient(np.rot90(data), horizontal=False, cmap=cmap),
        extent=(-0.5 * x_size, 0.5 * x_size, -0.5 * y_size, 0.5 * y_size),
        zorder=3,
    )
    ax.tripcolor(triang, zs.flatten(), zorder=0, cmap=cmap)
    return im, line1, line2

def generate_cover(
    learner,
    save_fname: Optional[str] = "thesis-cover.pdf",
    with_lines=False,
    npoints_interp=1000,
    dpi=300,
    cmap=None,
    personal_text=None,
    edition=None,
    with_text=True,
):
    data = list(learner.data.items())

    # Measured from proefdruk
    x_total = 34.95  # cm total sides + back
    y_total = 24  # cm top to bottom

    inch_per_cm = 2.54
    margin = 0.5  # add 5 mm margin on each side

    x_size = (x_total + margin) / inch_per_cm
    y_size = (y_total + margin) / inch_per_cm
    xy_size = x_size, y_size

    spine_size = 1.1 / inch_per_cm

    fig, ax = plt.subplots(figsize=(x_size, y_size))
    fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=None, hspace=None)

    ax.set_xticks([])
    ax.set_yticks([])

    cmap = cmap or get_cmap("inferno", 0.15, 0.95, 1.15)
    npoints_tri = len(data) // 4
    if len(data) > 4000:
        npoints_tri = max(npoints_tri, 4000)

    im, line1, line2 = get_new_artists(
        npoints_tri, learner, data, ax, xy_size, npoints_interp, cmap
    )

    title = "Towards realistic numerical simulations \n of Majorana devices"
    title2 = "Towards realistic numerical simulations of Majorana devices"
    author = "Bas Nijholt"

    text_color = "white"

    ax.axis("off")
    if with_text:
        font = "proxima_ssv/ProximaNova-Regular.otf"
        text_kwargs = dict(
            path_effects=[
                patheffects.withStroke(
                    linewidth=0.7, foreground="black", capstyle="round", alpha=1
                )
            ],
            zorder=4,
            verticalalignment="center",
            fontproperties=fm.FontProperties(fname=font),
        )
        for pos, text in zip([-0.8, 0.7], [author, title]):
            ax.text(
                x_size / 4,
                pos * (y_size - margin) / 2,
                text.upper(),
                color=text_color,
                weight="bold",
                **text_kwargs,
                horizontalalignment="center",
                fontsize=18,
            )

        ax.text(
            -0.09,
            y_size / 4 - 0.9,
            title2,
            color=text_color,
            weight="bold",
            rotation=-90,
            **text_kwargs,
            fontsize=12,
            horizontalalignment="left",
        )
        ax.text(
            -0.09,
            -y_size / 4 - 1,
            author,
            color=text_color,
            weight="bold",
            rotation=-90,
            **text_kwargs,
            fontsize=12,
            horizontalalignment="left",
        )

        lower_text_back = "Casimir PhD series 2020-11\nISBN 978-90-8593-438-7"
        if edition is not None:
            lower_text_back += f"\nedition {edition} of 120"
        ax.text(
            -x_size / 4,
            -0.8 * (y_size - margin) / 2,
            lower_text_back,
            color=text_color,
            weight="bold",
            horizontalalignment="center",
            **text_kwargs,
            fontsize=11,
        )
        if personal_text is not None:
            ax.text(
                -x_size / 4,
                0.4 * (y_size - margin) / 2,
                personal_text,
                color=text_color,
                weight="bold",
                horizontalalignment="center",
                **text_kwargs,
                fontsize=16,
            )

        if with_lines:
            for i in [-1, +1]:
                line_kwargs = dict(color="cyan", zorder=10, linestyles=":")
                ax.vlines(i * spine_size / 2, -y_size / 2, y_size / 2, **line_kwargs)
                ax.vlines(
                    -i * x_size / 2 + i * margin, -y_size / 2, y_size / 2, **line_kwargs
                )
                ax.hlines(
                    -i * y_size / 2 + i * margin, -x_size / 2, x_size / 2, **line_kwargs
                )

    ax.set_xlim(-x_size / 2, x_size / 2)
    ax.set_ylim(-y_size / 2, y_size / 2)
    print(f"Saving {save_fname}")
    if save_fname is not None:
        fig.savefig(
            save_fname,
            format=save_fname.suffix[1:],
            pad_inches=0,
            dpi=dpi,
        )
    plt.close(fig)

def bounds_from_saved_learner(fname):
    learner = adaptive.Learner2D(None, [(-1, 1), (-1, 1)])
    learner.load(fname)
    xs, ys = np.array(list(learner.data.keys())).T
    bounds = [(xs.min(), xs.max()), (ys.min(), ys.max())]
    return bounds

def load_learner(fname="data/mu-sweep2/data_learner_0246.pickle"):
    learner = adaptive.Learner2D(None, bounds_from_saved_learner(fname))
    learner.load(fname)
    return learner

def save(fname):
    print(f"Opening {fname}")
    f = fname.replace("/", "__")[:-7]
    pdf_fname = f"covers/{f}.pdf"
    print(pdf_fname)
    if os.path.exists(pdf_fname):
        print("exists, exit!")
        sys.exit(0)

    learner = load_learner(fname)
    generate_cover(learner, pdf_fname, with_lines=False, npoints_interp=2000)

def fname_out(folder, fname):
    fname_friendly = str(fname).replace("/", "__")
    return folder / f"{fname_friendly}.pdf"

max_stars_repo_name: basnijholt/thesis-cover

max_issues_repo_name: basnijholt/thesis-cover

max_forks_repo_name: basnijholt/thesis-cover


text:
# $Header: /opt/cvs/python/packages/share1.5/Pmv/crystalCommands.py,v 1.4 2009/05/22 18:40:03 vareille Exp $
#
# $Id: crystalCommands.py,v 1.4 2009/05/22 18:40:03 vareille Exp $
#
"""
Displays Unit Cell and Packing when Crystal Info is available in input cif file.
"""
from symserv.spaceGroups import spaceGroups
import numpy
from mglutil.math.crystal import Crystal
from Pmv.mvCommand import MVCommand
from ViewerFramework.VFCommand import CommandGUI
from mglutil.gui.InputForm.Tk.gui import InputFormDescr
import Pmw, Tkinter
import tkMessageBox
from DejaVu.Box import Box
from mglutil.util.callback import CallBackFunction
def instanceMatricesFromGroup(molecule):
    returnMatrices = [numpy.eye(4,4)]
    crystal = Crystal(molecule.cellLength, molecule.cellAngles)
    matrices = spaceGroups[molecule.spaceGroup]
    for matrix in matrices:
        tmpMatix = numpy.eye(4,4)
        tmpMatix[:3, :3] = matrix[0]
        tmpMatix[:3, 3] = crystal.toCartesian(matrix[1])
        returnMatrices.append(tmpMatix)
    molecule.crystal = crystal
    return returnMatrices

class CrystalCommand(MVCommand):
    def guiCallback(self):
        molNames = []
        for mol in self.vf.Mols:
            if hasattr(mol, 'spaceGroup'):
                molNames.append(mol.name)
        if not molNames:
            tkMessageBox.showinfo("Crystal Info is Needed", "No Molecule in the Viewer has Crystal Info.")
            return
        ifd = InputFormDescr(title='Crystal Info')
        ifd.append({'name':'moleculeList',
        'widgetType':Pmw.ScrolledListBox,
        'tooltip':'Select a molecule with Crystal Info.',
        'wcfg':{'label_text':'Select Molecule: ',
            'labelpos':'nw',
            'items':molNames,
            'listbox_selectmode':'single',
            'listbox_exportselection':0,
            'usehullsize': 1,
            'hull_width':100,'hull_height':150,
            'listbox_height':5},
        'gridcfg':{'sticky':'nsew', 'row':1, 'column':0}})
        val = self.vf.getUserInput(ifd, modal=1, blocking=1)
        if val:
            molecule = self.vf.getMolFromName(val['moleculeList'][0])
            matrices = instanceMatricesFromGroup(molecule)
            geom = molecule.geomContainer.geoms['master']
            geom.Set(instanceMatrices=matrices)
            if not molecule.geomContainer.geoms.has_key('Unit Cell'):
                fractCoords=((1,1,0),(0,1,0),(0,0,0),(1,0,0),(1,1,1),(0,1,1),
                                                                (0,0,1),(1,0,1))
                coords = []
                coords = molecule.crystal.toCartesian(fractCoords)
                box=Box('Unit Cell', vertices=coords)
                self.vf.GUI.VIEWER.AddObject(box, parent=geom)
                molecule.geomContainer.geoms['Unit Cell'] = box
            ifd = InputFormDescr(title='Crystal Options')
            visible = molecule.geomContainer.geoms['Unit Cell'].visible
            if visible:
                showState = 'active'
            else:
                showState = 'normal'
            ifd.append({'name': 'Show Cell',
            'widgetType':Tkinter.Checkbutton,
            'text': 'Hide Unit Cell',
            'state':showState,
            'gridcfg':{'sticky':Tkinter.W},
            'command': CallBackFunction(self.showUnitCell, molecule)})

            ifd.append({'name': 'Show Packing',
            'widgetType':Tkinter.Checkbutton,
            'text': 'Hide Packing',
            'state':'active',
            'gridcfg':{'sticky':Tkinter.W},
            'command': CallBackFunction(self.showPacking, molecule)})

            val = self.vf.getUserInput(ifd, modal=0, blocking=1)
            if not val:
                geom.Set(instanceMatrices=[numpy.eye(4,4)])
                molecule.geomContainer.geoms['Unit Cell'].Set(visible=False)

    def showUnitCell(self, molecule):
        visible = not molecule.geomContainer.geoms['Unit Cell'].visible
        molecule.geomContainer.geoms['Unit Cell'].Set(visible=visible)

    def showPacking(self, molecule):
        geom = molecule.geomContainer.geoms['master']
        if len(geom.instanceMatricesFortran) >= 2:
            geom.Set(instanceMatrices=[numpy.eye(4,4)])
        else:
            matrices = instanceMatricesFromGroup(molecule)
            geom.Set(instanceMatrices=matrices)

CrystalCommandGUI = CommandGUI()
CrystalCommandGUI.addMenuCommand('menuRoot', 'Display', 'Crystal')

commandList  = [{'name':'crystalCommand','cmd':CrystalCommand(),'gui':CrystalCommandGUI}]
def initModule(viewer):
    for _dict in commandList:
        viewer.addCommand(_dict['cmd'],_dict['name'],_dict['gui'])

max_stars_repo_name: J-E-J-S/aaRS-Pipeline

max_issues_repo_name: J-E-J-S/aaRS-Pipeline

max_forks_repo_name: J-E-J-S/aaRS-Pipeline


text:
import random
import cv2
from PIL import Image
import requests
from io import BytesIO
import numpy as np

def addPreset(im, union_mask):
	green = random.choice([40, 100, 240, 200, 80])
	temp = im
	h1_mask = union_mask
	WHITE_BORDER_FRACTION = 0.07
	for i in range(h1_mask.shape[0]):
		for j in range(h1_mask.shape[1]):

			if(h1_mask[i][j]==0):
				if(i>int(h1_mask.shape[0]*WHITE_BORDER_FRACTION) and i<int((1-WHITE_BORDER_FRACTION)*h1_mask.shape[0]) and j>int(h1_mask.shape[1]*WHITE_BORDER_FRACTION) and j<int((1-WHITE_BORDER_FRACTION)*h1_mask.shape[1])):
					red = 255*(1-(j/h1_mask.shape[1]))
					temp[i][j] = [red, green, 255]
				else:
					temp[i][j] = [255, 255, 255]
	return temp

def createThumbnail(filename):
	'''
		createThumbnail('me.jpg')
		creates lower resolution png format image
	'''
	from PIL import Image
	path = './'
	size = 1080, 1296
	im = Image.open(path+filename)
	im.thumbnail(size, Image.ANTIALIAS)
	im.save(path+filename.split('.')[0]+'.png', "PNG")

def loadImage(URL):
	# with urllib.request.urlopen(URL) as url:
	#	 img = tf.keras.preprocessing.image.load_img(BytesIO(url.read()), target_size=(w, h))
	response = requests.get(URL)
	img = Image.open(BytesIO(response.content))
	return np.array(img)

def feature_1(detector, filename):
	# Remove bg and output image with bg-color-preset
	im, human_masks, boxes = detector.return_attributes(filename)

	# createThumbnail('me_winter.jpg')
	final_mask = detector.return_union_mask(im, human_masks)
	im = addPreset(im, final_mask)
	return im
	# cv2.imwrite('./x_'+filename.split('./')[1], im)

def feature_2(detector, filename, catagory):
	'''
		Add unsplash image as background.
		TO-DO: run detector and loadImage in parallel (using threads or processes)
	'''
	im, human_masks, boxes = detector.return_attributes(filename)
	final_mask = detector.return_union_mask(im, human_masks)
	print('Segment mask calculated')

	bg = loadImage('https://source.unsplash.com/1500x1500/?'+catagory)
	# bg = cv2.imread('./background.jpeg')
	print('Unsplash image loaded')

	WHITE_BORDER_FRACTION = 0.07
	for i in range(final_mask.shape[0]):
		for j in range(final_mask.shape[1]):

			if(final_mask[i][j]==0):
				if(i>int(final_mask.shape[0]*WHITE_BORDER_FRACTION) and i<int((1-WHITE_BORDER_FRACTION)*final_mask.shape[0]) and j>int(final_mask.shape[1]*WHITE_BORDER_FRACTION) and j<int((1-WHITE_BORDER_FRACTION)*final_mask.shape[1])):
					red = 255*(1-(j/final_mask.shape[1]))
					im[i][j] = bg[i,j,:]
				else:
					im[i][j] = [255, 255, 255]
	return im

def feature_3(segment_detector, face_detector, smile_detector, filename, catagory):
	'''
		Add smile related caption to the image
	'''
	WHITE_BORDER_FRACTION = 0.07

	im, human_masks, boxes = segment_detector.return_attributes(filename)
	final_mask = segment_detector.return_union_mask(im, human_masks)

	print('Segment mask calculated')

	# bg = loadImage('https://source.unsplash.com/1500x1500/?'+catagory)

	bg = cv2.imread('./background.jfif')

	print('Unsplash image loaded')

	############################### ###############################
	# pred_boxes = face_detector.return_pred_boxes(im)

	# print('Face boxes calculated')

	# num_faces = pred_boxes.shape[0]
	# prob_smiling_faces = []

	# for i in pred_boxes:
	# 	temp = im[i[1]:i[3],i[0]:i[2],:]
	# 	proba = smile_detector._return_smile_prob(temp)
	# 	prob_smiling_faces.append(proba[0][0]) # Smiling probability of each person

	# prob_smiling_faces = np.array(prob_smiling_faces)
	# index_smile_most = np.argmax(prob_smiling_faces) # Person smiling most
	# smiling_prob_final = prob_smiling_faces[index_smile_most] # Smiling probability of person smiling most

	# face_final = pred_boxes[index_smile_most]

	print('Removing background and adding preset')
	for i in range(final_mask.shape[0]):
		for j in range(final_mask.shape[1]):

			if(final_mask[i][j]==0):
				if(i>int(final_mask.shape[0]*WHITE_BORDER_FRACTION) and i<int((1-WHITE_BORDER_FRACTION)*final_mask.shape[0]) and j>int(final_mask.shape[1]*WHITE_BORDER_FRACTION) and j<int((1-WHITE_BORDER_FRACTION)*final_mask.shape[1])):
					red = 255*(1-(j/final_mask.shape[1]))
					im[i][j] = bg[i,j,:]
				else:
					im[i][j] = [255, 255, 255]

	print(im.shape[1]*WHITE_BORDER_FRACTION)

	smiling_prob_final = smile_detector._return_smile_prob(im)[0][1]

	print('Smiling prediction: ', smiling_prob_final)

	isSmiling = smiling_prob_final > 0.5
	if isSmiling:
		textOnImage = 'Smiling: ' + str('%.2f' % smiling_prob_final)
		textOnImage = '#KEEPSMILING ;)'
	else:
		textOnImage = 'Not Smiling: ' + str('%.2f' % (1-smiling_prob_final))
		textOnImage = '#OFFMOOOOD'
	############################### ###############################

	font				   = cv2.FONT_HERSHEY_SIMPLEX
	bottomLeftCornerOfText = (100,700)
	fontScale			  = 3
	fontColor			  = (10,0,250)
	lineType			   = 3

	cv2.putText(
		im,
		textOnImage,
		bottomLeftCornerOfText,
		font,
		fontScale,
		fontColor,
		lineType
	)

	return im


text:
import numpy as np
import pandas as pd
import random
import Scripts.PopUpManager as PopUpManager
import time

class FacultyAssigner:
    def __init__(self, input_file_path, destination, root, progress):
        self.function = progress
        self.root = root
        self.input_file_path = input_file_path
        self.destination = destination
        self.df = pd.read_csv(input_file_path)
        self.faculty_list = list(self.df.Faculty)
        df = self.df.dropna()
        self.dates = list(df.Dates)
        self.slot1 = list(df.Slot1)
        self.slot2 = list(df.Slot2)
        self.slot3 = list(df.Slot3)
        self.max_occurrence = ((sum(self.slot1) + sum(self.slot2)) // len(self.faculty_list)) + 1
        self.faculty_key = {}
        for i in self.faculty_list:
            self.faculty_key[i] = 0

    def get_slot(self, slot):
        df = pd.DataFrame({})
        for i in range(len(self.dates)):
            data = []
            for _ in range(int(slot[i])):
                choice = None
                while choice is None:
                    faculty = random.choice(self.faculty_list)
                    if (self.faculty_key[faculty] < self.max_occurrence) and (faculty not in data):
                        self.faculty_key[faculty] += 1
                        choice = faculty
                    else:
                        continue
                if choice is None:
                    PopUpManager.error_popup(self.root, 'Some unexpected error occurred')
                else:
                    data.append(choice)
            for y in range(len(self.faculty_list)-len(data)):
                data.append(np.nan)
            df[self.dates[i]] = data
        time.sleep(1)

        return df

    def generate_check_file(self, slot):
        df = pd.DataFrame({}, columns=self.dates, index=self.faculty_list)
        for i in self.dates:
            for j in self.faculty_list:
                if j in list(slot[i]):
                    df[i][j] = '✓'
                else:
                    df[i][j] = np.nan
        return df

    def runner(self):
        slot1 = self.get_slot(self.slot1)
        slot2 = self.get_slot(self.slot2)
        file1 = self.generate_check_file(slot1)
        file2 = self.generate_check_file(slot2)
        file1.to_excel(self.destination+'/Slot1.xlsx')
        file2.to_excel(self.destination+'/Slot2.xlsx')
        self.function()


text:
#!/usr/bin/env python
# @Copyright 2007 Kristjan Haule
from scipy import *
from scipy import linalg
import struct, sys, re
import rdU, utils
import findEF

if len(sys.argv)<=1:
    print 'Give dictionary of coh_orb'
    sys.exit(0)

#coh_orb = {0:2,1:2}
coh_orb = eval(sys.argv[1])
print 'coh_orb=', coh_orb

fhp = 233
cpuID = 0

dmfe = utils.DmftEnvironment()  # DMFT paths
w2k = utils.W2kEnvironment()    # W2k filenames and paths
case = w2k.case
(EF,NOE) = findEF.FindChemicalPotential(case,'')

# Read Udmft
filename = 'Udmft.'+str(cpuID)
(nkp, nsymop, norbitals) = rdU.fopen(fhp, filename)
nindo = rdU.read1(fhp, norbitals)

print 'nsymop=', nsymop

fc = open('cohfactors.dat', 'w')

print >> fc, '# Coherence factors for KS states'

for ikp in range(nkp):

    (iikp, nbands, tmaxdim2, tnorbitals, nemin) = rdU.read2(fhp)

    print 'ikp=', ikp,  'nbands=', nbands, 'maxdim2=', tmaxdim2, 'norbitals=', tnorbitals, 'nemin=', nemin

    for isym in range(nsymop):

        iisym = rdU.read3(fhp)
        if iisym!=isym+1:
            print 'ERROR: isym!=iisym', isym, iisym

        gs = zeros((nbands,nbands),dtype=complex)

        for iorb in range(norbitals):
            # Reading DMFT transformation UDMFT
            U = []
            for ind in range(nindo[iorb]):
                U.append( rdU.read4(fhp, nbands) )
            U = array(U)
            # print ikp, isym, iorb, shape(U)

            if (isym==0):
                # Prepares coherence factors gs(i,j)=U(orb,i)*U(orb,j).conj()
                if coh_orb.has_key(iorb):
                    xorb = coh_orb[iorb]
                    for iband in range(nbands):
                        gs[iband,:] += U[xorb,iband]*conj(U[xorb,:])

        if isym==0:
            # Printing the coherence factors
            print >> fc, iikp, nemin, nbands
            for i in range(len(gs)):
                for j in range(len(gs)):
                    print >> fc, "%15.10f " % abs(gs[j,i]),
                print >> fc

rdU.fclose(fhp)

reason: from scipy


text:
import numpy as np
import math
from functools import partial
import torch

class LRSchedulerStep(object):
    def __init__(self, fai_optimizer, total_step, lr_phases, mom_phases):
        self.optimizer = fai_optimizer
        self.total_step = total_step
        self.lr_phases = []

        for i, (start, lambda_func) in enumerate(lr_phases):
            if len(self.lr_phases) != 0:
                assert self.lr_phases[-1][0] < int(start * total_step)
            if isinstance(lambda_func, str):
                lambda_func = eval(lambda_func)
            if i < len(lr_phases) - 1:
                self.lr_phases.append((int(start * total_step),
                                       int(lr_phases[i + 1][0] * total_step),
                                       lambda_func))
            else:
                self.lr_phases.append((int(start * total_step), total_step,
                                       lambda_func))
        assert self.lr_phases[0][0] == 0
        self.mom_phases = []

        for i, (start, lambda_func) in enumerate(mom_phases):
            if len(self.mom_phases) != 0:
                assert self.mom_phases[-1][0] < int(start * total_step)
            if isinstance(lambda_func, str):
                lambda_func = eval(lambda_func)
            if i < len(mom_phases) - 1:
                self.mom_phases.append((int(start * total_step),
                                        int(mom_phases[i + 1][0] * total_step),
                                        lambda_func))
            else:
                self.mom_phases.append((int(start * total_step), total_step,
                                        lambda_func))
        if len(mom_phases) > 0:
            assert self.mom_phases[0][0] == 0

    def step(self, step):
        lrs = []
        moms = []
        for start, end, func in self.lr_phases:
            if step >= start:
                lrs.append(func((step - start) / (end - start)))
        if len(lrs) > 0:
            self.optimizer.lr = lrs[-1]
        for start, end, func in self.mom_phases:
            if step >= start:
                moms.append(func((step - start) / (end - start)))
                self.optimizer.mom = func((step - start) / (end - start))
        if len(moms) > 0:
            self.optimizer.mom = moms[-1]

    @property
    def learning_rate(self):
        return self.optimizer.lr

def annealing_cos(start, end, pct):
    # print(pct, start, end)
    "Cosine anneal from `start` to `end` as pct goes from 0.0 to 1.0."
    cos_out = np.cos(np.pi * pct) + 1
    return end + (start - end) / 2 * cos_out

class OneCycle(LRSchedulerStep):
    def __init__(self, fai_optimizer, total_step, lr_max, moms, div_factor,
                 pct_start):
        self.lr_max = lr_max
        self.moms = moms
        self.div_factor = div_factor
        self.pct_start = pct_start
        a1 = int(total_step * self.pct_start)
        a2 = total_step - a1
        low_lr = self.lr_max / self.div_factor
        lr_phases = ((0, partial(annealing_cos, low_lr, self.lr_max)),
                     (self.pct_start,
                      partial(annealing_cos, self.lr_max, low_lr / 1e4)))
        mom_phases = ((0, partial(annealing_cos, *self.moms)),
                      (self.pct_start, partial(annealing_cos,
                                               *self.moms[::-1])))
        fai_optimizer.lr, fai_optimizer.mom = low_lr, self.moms[0]
        super().__init__(fai_optimizer, total_step, lr_phases, mom_phases)

class ExponentialDecay(LRSchedulerStep):
    def __init__(self,
                 fai_optimizer,
                 total_step,
                 initial_learning_rate,
                 decay_length,
                 decay_factor,
                 staircase=True):
        """
        Args:
            decay_length: must in (0, 1)
        """
        assert decay_length > 0
        assert decay_length < 1
        self._decay_steps_unified = decay_length
        self._decay_factor = decay_factor
        self._staircase = staircase
        step = 0
        stage = 1
        lr_phases = []
        if staircase:
            while step <= total_step:
                func = lambda p, _d=initial_learning_rate * stage: _d
                lr_phases.append((step / total_step, func))
                stage *= decay_factor
                step += int(decay_length * total_step)
        else:
            func = lambda p: pow(decay_factor, (p / decay_length))
            lr_phases.append((0, func))
        super().__init__(fai_optimizer, total_step, lr_phases, [])

class ManualStepping(LRSchedulerStep):
    def __init__(self, fai_optimizer, total_step, boundaries, rates):
        assert all([b > 0 and b < 1 for b in boundaries])
        assert len(boundaries) + 1 == len(rates)
        boundaries.insert(0, 0.0)
        lr_phases = []
        for start, rate in zip(boundaries, rates):
            func = lambda p, _d=rate: _d
            lr_phases.append((start, func))
        super().__init__(fai_optimizer, total_step, lr_phases, [])

class FakeOptim:
    def __init__(self):
        self.lr = 0
        self.mom = 0

if __name__ == "__main__":
    import matplotlib.pyplot as plt
    opt = FakeOptim()  # 3e-3, wd=0.4, div_factor=10
    # schd = OneCycle(opt, 100, 3e-3, (0.95, 0.85), 10.0, 0.4)
    schd = ExponentialDecay(opt, 100, 3e-4, 0.1, 0.8, staircase=True)
    schd = ManualStepping(opt, 100, [0.8, 0.9], [0.001, 0.0001, 0.00005])
    lrs = []
    moms = []
    for i in range(100):
        schd.step(i)
        lrs.append(opt.lr)
        moms.append(opt.mom)

    plt.plot(lrs)
    # plt.plot(moms)
    # plt.show()
    # plt.plot(moms)
    plt.show()


text:
from __future__ import division
import sys
import os
import numpy as np
import matplotlib.pyplot as plt
from math import ceil
from statistics import mode
from statsmodels.distributions.empirical_distribution import ECDF

try:
    from xml.etree import cElementTree as ElementTree
except ImportError:
    from xml.etree import ElementTree

def parse_time_ns(tm):
    if tm.endswith('ns'):
        return int(tm[:-4])
    raise ValueError(tm)

## FiveTuple
class FiveTuple(object):
    ## class variables
    ## @var sourceAddress
    #  source address
    ## @var destinationAddress
    #  destination address
    ## @var protocol
    #  network protocol
    ## @var sourcePort
    #  source port
    ## @var destinationPort
    #  destination port
    ## @var __slots__
    #  class variable list
    __slots__ = ['sourceAddress', 'destinationAddress', 'protocol', 'sourcePort', 'destinationPort']

    def __init__(self, el):
        '''The initializer.
        @param self The object pointer.
        @param el The element.
        '''
        self.sourceAddress = el.get('sourceAddress')
        self.destinationAddress = el.get('destinationAddress')
        self.sourcePort = int(el.get('sourcePort'))
        self.destinationPort = int(el.get('destinationPort'))
        self.protocol = int(el.get('protocol'))

## Histogram
class Histogram(object):
    ## class variables
    ## @var bins
    #  histogram bins
    ## @var nbins
    #  number of bins
    ## @var number_of_flows
    #  number of flows
    ## @var __slots__
    #  class variable list
    __slots__ = 'bins', 'nbins', 'number_of_flows'

    def __init__(self, el=None):
        ''' The initializer.
        @param self The object pointer.
        @param el The element.
        '''
        self.bins = []
        if el is not None:
            # self.nbins = int(el.get('nBins'))
            for bin in el.findall('bin'):
                self.bins.append((float(bin.get("start")), float(bin.get("width")), int(bin.get("count"))))

## Flow
class Flow(object):
    ## class variables
    ## @var flowId
    #  delay ID
    ## @var delayMean
    #  mean delay
    ## @var packetLossRatio
    #  packet loss ratio
    ## @var rxBitrate
    #  receive bit rate
    ## @var txBitrate
    #  transmit bit rate
    ## @var fiveTuple
    #  five tuple
    ## @var packetSizeMean
    #  packet size mean
    ## @var probe_stats_unsorted
    #  unsirted probe stats
    ## @var hopCount
    #  hop count
    ## @var flowInterruptionsHistogram
    #  flow histogram
    ## @var rx_duration
    #  receive duration
    ## @var __slots__
    #  class variable list
    __slots__ = ['flowId', 'delayMean', 'packetLossRatio', 'rxBitrate', 'txBitrate',
                 'fiveTuple', 'packetSizeMean', 'probe_stats_unsorted',
                 'hopCount', 'flowInterruptionsHistogram', 'rx_duration', 'rxPackets']

    def __init__(self, flow_el):
        ''' The initializer.
        @param self The object pointer.
        @param flow_el The element.
        '''
        self.flowId = int(flow_el.get('flowId'))
        rxPackets = int(flow_el.get('rxPackets'))
        txPackets = int(flow_el.get('txPackets'))
        tx_duration = float(
            int(flow_el.get('timeLastRxPacket')[:-4]) - int(flow_el.get('timeFirstTxPacket')[:-4])) * 1e-9
        rx_duration = float(
            int(flow_el.get('timeLastRxPacket')[:-4]) - int(flow_el.get('timeFirstRxPacket')[:-4])) * 1e-9
        self.rx_duration = rx_duration
        self.probe_stats_unsorted = []
        if rxPackets:
            self.hopCount = float(flow_el.get('timesForwarded')) / rxPackets + 1
            self.rxPackets = rxPackets
        else:
            self.hopCount = -1000
        if rxPackets:
            self.delayMean = float(flow_el.get('delaySum')[:-4]) / rxPackets * 1e-9
            self.packetSizeMean = float(flow_el.get('rxBytes')) / rxPackets
        else:
            self.delayMean = None
            self.packetSizeMean = None
        if rx_duration > 0:
            self.rxBitrate = int(flow_el.get('rxBytes')) * 8 / rx_duration
        else:
            self.rxBitrate = None
        if tx_duration > 0:
            self.txBitrate = int(flow_el.get('txBytes')) * 8 / tx_duration
        else:
            self.txBitrate = None
        lost = float(flow_el.get('lostPackets'))
        # print "rxBytes: %s; txPackets: %s; rxPackets: %s; lostPackets: %s" % (flow_el.get('rxBytes'), txPackets, rxPackets, lost)
        if rxPackets == 0:
            self.packetLossRatio = None
        else:
            self.packetLossRatio = (lost / (rxPackets + lost))

        interrupt_hist_elem = flow_el.find("flowInterruptionsHistogram")
        if interrupt_hist_elem is None:
            self.flowInterruptionsHistogram = None
        else:
            self.flowInterruptionsHistogram = Histogram(interrupt_hist_elem)

## ProbeFlowStats
class ProbeFlowStats(object):
    ## class variables
    ## @var probeId
    #  probe ID
    ## @var packets
    #  network packets
    ## @var bytes
    #  bytes
    ## @var delayFromFirstProbe
    #  delay from first probe
    ## @var __slots__
    #  class variable list
    __slots__ = ['probeId', 'packets', 'bytes', 'delayFromFirstProbe']

'''
/home/zwa/xq/cce-scripts/venv/bin/python /home/zwa/xq/cce-scripts/src/cce/parser/flow_parser.py
{
    'name': 'FlowMonitor',
    'count': 1,
    'children': [
        {
            'name': 'FlowStats',
            'count': 1,
            'attribute': {},
            'children': [
                {
                    'name': 'Flow',
                    'count': 10000,
                    'attribute': {
                        'flowId': '1',
                        'timeFirstTxPacket': '+0.0ns',
                        'timeFirstRxPacket': '+8472.0ns',
                        'timeLastTxPacket': '+135.0ns',
                        'timeLastRxPacket': '+11768.0ns',
                        'delaySum': '+20105.0ns',
                        'jitterSum': '+3161.0ns',
                        'lastDelay': '+11633.0ns',
                        'txBytes': '2056',
                        'rxBytes': '2056',
                        'txPackets': '2',
                        'rxPackets': '2',
                        'lostPackets': '0',
                        'timesForwarded': '10'
                    },
                    'children': [
                        {
                            'name': 'delayHistogram',
                            'count': 1,
                            'attribute': {'nBins': '1'},
                            'children': [
                                {
                                    'name': 'bin',
                                    'count': 1,
                                    'attribute': {
                                        'index': '0',
                                        'start': '0',
                                        'width': '0.001',
                                        'count': '2'
                                    },
                                    'children': []
                                }
                            ]
                        },
                        {
                            'name': 'jitterHistogram',
                            'count': 1,
                            'attribute': {'nBins': '1'},
                            'children': [
                                {
                                    'name': 'bin',
                                    'count': 1,
                                    'attribute': {
                                        'index': '0',
                                        'start': '0',
                                        'width': '0.001',
                                        'count': '1'
                                    },
                                    'children': []
                                }
                            ]
                        },
                        {
                            'name': 'packetSizeHistogram',
                            'count': 1,
                            'attribute': {'nBins': '52'},
                            'children': [
                                {
                                    'name': 'bin',
                                    'count': 1,
                                    'attribute': {
                                        'index': '51',
                                        'start': '1020',
                                        'width': '20',
                                        'count': '2'
                                    },
                                    'children': []
                                }
                            ]
                        },
                        {
                            'name': 'flowInterruptionsHistogram',
                            'count': 1,
                            'attribute': {'nBins': '0'},
                            'children': []
                        }
                    ]
                }
            ]
        },
        {
            'name': 'Ipv4FlowClassifier',
            'count': 1,
            'attribute': {},
            'children': [
                {
                    'name': 'Flow',
                    'count': 10000,
                    'attribute': {
                        'flowId': '1',
                        'sourceAddress': '1.1',
                        'destinationAddress': '1.161',
                        'protocol': '17',
                        'sourcePort': '49153',
                        'destinationPort': '9199'
                    },
                    'children': []
                }
            ]
        },
        {
            'name': 'FlowProbes',
            'count': 1,
            'attribute': {},
            'children': [
                {
                    'name': 'FlowProbe',
                    'count': 288,
                    'attribute': {'index': '0'},
                    'children': [
                        {
                            'name': 'FlowStats',
                            'count': 1000,
                            'attribute': {
                                'flowId': '1',
                                'packets': '2',
                                'bytes': '2056',
                                'delayFromFirstProbeSum': '+0.0ns'
                            },
                            'children': []
                        }
                    ]
                }
            ]
        }
    ]
}

Process finished with exit code 0
'''
## Simulation
class Simulation(object):
    ## class variables
    ## @var flows
    #  list of flows
    def __init__(self, simulation_el):
        ''' The initializer.
        @param self The object pointer.
        @param simulation_el The element.
        '''
        self.flows = []
        FlowClassifier_el, = simulation_el.findall("Ipv4FlowClassifier")
        flow_map = {}
        for flow_el in simulation_el.findall("FlowStats/Flow"):
            flow = Flow(flow_el)
            flow_map[flow.flowId] = flow
            self.flows.append(flow)
        for flow_cls in FlowClassifier_el.findall("Flow"):
            flowId = int(flow_cls.get('flowId'))
            flow_map[flowId].fiveTuple = FiveTuple(flow_cls)

        for probe_elem in simulation_el.findall("FlowProbes/FlowProbe"):
            probeId = int(probe_elem.get('index'))
            for stats in probe_elem.findall("FlowStats"):
                flowId = int(stats.get('flowId'))
                s = ProbeFlowStats()
                s.packets = int(stats.get('packets'))
                s.bytes = int(stats.get('bytes'))
                s.probeId = probeId
                if s.packets > 0:
                    s.delayFromFirstProbe = parse_time_ns(stats.get('delayFromFirstProbeSum')) / float(s.packets)
                else:
                    s.delayFromFirstProbe = 0
                flow_map[flowId].probe_stats_unsorted.append(s)

def ProcessXMLFile(filename):
    file_obj = open(filename)
    # print ("Reading XML file ",filename)

    sys.stdout.flush()
    level = 0
    sim_list = []
    for event, elem in ElementTree.iterparse(file_obj, events=("start", "end")):
        if event == "start":
            level += 1
        if event == "end":
            level -= 1
            if level == 0 and elem.tag == 'FlowMonitor':
                sim = Simulation(elem)
                sim_list.append(sim)
                elem.clear()  # won't need this any more
                # sys.stdout.write(".")
                sys.stdout.flush()
    # print (" done.")

    th = []
    rth = []
    wth = []

    for sim in sim_list:
        sum = 0
        # th=[]
        # print()
        for flow in sim.flows:
            t = flow.fiveTuple
            proto = {6: 'TCP', 17: 'UDP'}[t.protocol]
            if flow.txBitrate is None:
                print("\tTX bitrate: None")
            else:
                th.append(flow.txBitrate)
                if t.destinationAddress == '1.161':
                    wth.append(flow.txBitrate)
                else:
                    rth.append(flow.txBitrate)

    return th, rth, wth

def checkSize(filename):
    file_obj = open(filename)
    # print ("Reading XML file ",filename)

    sys.stdout.flush()
    level = 0
    sim_list = []
    for event, elem in ElementTree.iterparse(file_obj, events=("start", "end")):
        if event == "start":
            level += 1
        if event == "end":
            level -= 1
            if level == 0 and elem.tag == 'FlowMonitor':
                sim = Simulation(elem)
                sim_list.append(sim)
                elem.clear()  # won't need this any more
                # sys.stdout.write(".")
                sys.stdout.flush()
    # print (" done.")
    print(len(sim_list))
    size = []
    rsize = []
    wsize = []

    for sim in sim_list:
        sum = 0
        # th=[]
        # print()
        for flow in sim.flows:
            t = flow.fiveTuple
            proto = {6: 'TCP', 17: 'UDP'}[t.protocol]
            if flow.rxPackets is None:
                print("\tTX bitrate: None")
            else:
                size.append(flow.rxPackets)
                if t.destinationAddress == '1.161':
                    wsize.append(flow.rxPackets)
                else:
                    rsize.append(flow.rxPackets)

    return size, rsize, wsize

def printThroughputs(th):
    print(len(th), np.mean(th), np.median(th), np.percentile(th, 10))

def Compare2(th, th1):
    print(len(th), np.mean(th), np.median(th), np.percentile(th, 10))  # mode(th), th.count(mode(th)))
    print(len(th1), np.mean(th1), np.median(th1), np.percentile(th1, 10))  # mode(th1), th1.count(mode(th1)))
    # bin_size=5000
    # mx=int(max(max(th),max(th1)))
    # mn=int(min(min(th),min(th1)))
    # bin_num=ceil((mx-mn) / bin_size)
    # ecdf=ECDF(th)
    # ecdf1=ECDF(th1)
    # plt.plot(ecdf.x, ecdf.y,'r',ecdf1.x, ecdf1.y,'b')
    # plt.show()

def DrawCDF(th1, bin_size):
    f, axs = plt.subplots(3, 1, sharex=True)
    for x, th in enumerate(th1):
        print(len(th), np.mean(th), mode(th), th.count(mode(th)))
        mx = int(max(th))
        mn = int(min(th))
        bin_num = ceil((mx - mn) / bin_size)
        ecdf = ECDF(th)
        axs[x].plot(ecdf.x, ecdf.y)

    plt.xlim(0, 100)
    plt.show()

def main(argv):
    s, rs, ws = checkSize(argv[1])
    # DrawCDF([s,rs,ws],1)

    th, rth, wth = ProcessXMLFile(argv[1])
    # th1,rth1,wth1 = ProcessXMLFile(argv[2])

    # Compare2(th,th1)
    printThroughputs(th)

    # Compare2(rth,rth1)
    # Compare2(wth,wth1)
    # print ("\tTX bitrate: %.2f kbit/s" % (flow.txBitrate*1e-3,))
    # if flow.rxBitrate is None:
    #    print ("\tRX bitrate: None")
    # else:
    #    print ("\tRX bitrate: %.2f kbit/s" % (flow.rxBitrate*1e-3,))
    # if flow.delayMean is None:
    #    print ("\tMean Delay: None")
    # else:
    #    print ("\tMean Delay: %.2f ms" % (flow.delayMean*1e3,))
    # if flow.packetLossRatio is None:
    #    print ("\tPacket Loss Ratio: None")
    # else:
    #    print ("\tPacket Loss Ratio: %.2f %%" % (flow.packetLossRatio*100))
    # sum=sum+flow.txBitrate*1e-3
    # print(sum/len(sim.flows))

if __name__ == '__main__':
    # main(sys.argv)
    checkSize('../../../data/flow_monitor.xml')

max_stars_repo_name: TiffanyXQZ/cce-scripts

max_issues_repo_name: TiffanyXQZ/cce-scripts

max_forks_repo_name: TiffanyXQZ/cce-scripts

reason: import numpy,from statsmodels