text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,py:light
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.1.7
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# +
import cv2
import matplotlib.pyplot as plt
import tensorflow as tf
import os
import numpy as np
import time
import sys
if 'lib\\models\\research\\' not in sys.path:
    sys.path.append('lib\\models\\research\\')

import object_detection.utils.visualization_utils
from google.protobuf import text_format
from object_detection.protos import string_int_label_map_pb2
import pickle

# +
with open('label_map.pbtxt', 'rb') as f:
    label_map = string_int_label_map_pb2.StringIntLabelMap()
    text_format.Merge(f.read(), label_map)

label_map_dict = {i.id: {'name' : i.name} for i in label_map.item}

# -

model = tf.saved_model.load('data/frcnn_model/saved_model')
model = model.signatures['serving_default']

out_fn = 'data/out1.avi'

fn = 'data/eval/Videos/видео-20200124-100159-d6ad16e5.MOV'
vid = cv2.VideoCapture(fn)
if not vid.isOpened(): raise ValueError("Error opening file '" + fn + "'")
try:
    with open(out_fn + '.boxes1', 'wb') as out_boxes:
        while True:
            ret, frame = vid.read()
            if not ret: break
            frame = np.flip(np.transpose(frame, [1,0,2]), 1)
            im = tf.convert_to_tensor(frame)[tf.newaxis,...]
            res = model(im)
            for nm in ['detection_scores', 'detection_boxes', 'detection_classes']:
                pickle.dump(res[nm].numpy(), out_boxes)
                print('frame')
finally:
    vid.release()

fn = 'data/eval/Videos/видео-20200124-100159-d6ad16e5.MOV'
vid = cv2.VideoCapture(fn)
if not vid.isOpened(): raise ValueError("Error opening file '" + fn + "'")
try:
    width = int(vid.get(cv2.CAP_PROP_FRAME_WIDTH))
    height = int(vid.get(cv2.CAP_PROP_FRAME_HEIGHT))
    if os.path.exists(out_fn): os.remove(out_fn)
    out = cv2.VideoWriter(
        out_fn,
        fourcc=cv2.VideoWriter_fourcc(*'MP42'),
        fps=float(vid.get(cv2.CAP_PROP_FPS)),
#         frameSize = (width, height),
        frameSize = (height, width),
    )
    try:
        with open(out_fn + '.boxes', 'rb') as in_boxes:
            while True:
                ret, frame = vid.read()
                if not ret: break
                frame = np.flip(np.transpose(frame, [1,0,2]), 1)

                detection_scores, detection_boxes, detection_classes = [pickle.load(in_boxes) for _ in range(3)]
                flt = detection_scores > .8
                im = object_detection.utils.visualization_utils.visualize_boxes_and_labels_on_image_array(
                      frame,
                      detection_boxes[flt],
                      detection_classes[flt].astype(np.int),
                      detection_scores[flt],
                      label_map_dict,
                      instance_masks=None,
                      use_normalized_coordinates=True,
                      line_thickness=8
                )
                out.write(im)
                print("frame")
    finally:
        out.release()
finally:
    vid.release()


text:
import numpy as np

class RescaleVectors:
	def __init__(self):
		pass

	def Rescale(self,X,new_min,new_max):
	#Rescaling:
        #This function rescales values so that the vector in a matrix contains values in
        #the range new_min <= xx <= new_max

		SizeX = X.shape
                X_resc = np.ones(SizeX)
                pdim = SizeX[1]
		MinVec = []
		MaxVec = []
		MultFacVec = []
		for i in range(pdim):
			old_min = np.min(X[:,i])
			old_max = np.max(X[:,i])
			MinVec.append(old_min)
			MaxVec.append(old_max)
			mult_fac = (new_max-new_min)/(old_max-old_min)
			MultFacVec.append(mult_fac)
			X_resc[:,i] = mult_fac*(X[:,i]-old_min) + new_min;
		self.X = X_resc
		self.MIN = MinVec
		self.MAX = MaxVec
		self.MULTFACT = MultFacVec

	def Standardization(self,X):
	#Standardization: difference with mean value divided by the
        #standard deviation

		mean_p = []
    		std_p = []
		SizeX = X.shape
    		X_stand = np.ones(SizeX)
		pdim = SizeX[1]
    		for i in range(pdim):
        		avgp = np.mean(X[:,i])
        		stdp = np.std(X[:,i])
        		mean_p.append(avgp)
        		std_p.append(stdp)
        		X_stand[:,i] = (X[:,i] - avgp)/stdp

		self.X = X_stand
		self.MEAN = mean_p
		self.STD = std_p

ext: py

max_stars_repo_name: L-F-A/Changing-matrix-vector

max_issues_repo_name: L-F-A/Changing-matrix-vector

max_forks_repo_name: L-F-A/Changing-matrix-vector

reason: import numpy


text:
# -*- coding: utf-8 -*-
"""User-friendly top-level functions which allow the user to handle

#. Call-background segmentation
#. CF-FM call part segmentation
#. Measurement of CF-FM audio parts

Let's take a look at an example where we [TO BE COMPLETED!!!]

.. code-block:: python

    import scipy.signal as signal
    from itsfm.user_interface import segment_and_measure_call
    from itsfm.view_horseshoebat_call import *
    from itsfm.simulate_calls import make_cffm_call

    # create synthetic call
    call_parameters = {'cf':(100000, 0.01),
                        'upfm':(80000, 0.002),
                        'downfm':(60000, 0.003),
                        }

    fs = 500*10**3 # 500kHz sampling rate
    synthetic_call, freq_profile = make_cffm_call(call_parameters, fs)

    # window and reduce overall signal level
    synthetic_call *= signal.tukey(synthetic_call.size, 0.1)
    synthetic_call *= 0.75

    # measuring a well-selected call (without silent background)

    # measuing a call with a silent background

    # and add 2ms of additional background_noise of ~ -60dBrms
    samples_1ms = int(fs*0.001)
    final_size = synthetic_call.size + samples_1ms*2
    call_with_noise = np.random.normal(0,10**(-60/20.0),final_size)
    call_with_noise[samples_1ms:-samples_1ms] += synthetic_call

    #

    seg_and_msmts = segment_and_measure_call(call_with_noise, fs,
                                             segment_from_background=True)
    call_segmentation, call_parts, measurements, backg_segment = seg_and_msmts

"""

import matplotlib
import matplotlib.pyplot as plt
import matplotlib.backends.backend_pdf
import itsfm.segment
from itsfm.segment import segment_call_from_background
from itsfm.segment import segment_call_into_cf_fm
from itsfm.measure import measure_hbc_call
from itsfm.measure import parse_cffm_segments

def segment_and_measure_call(main_call, fs,
                             segment_from_background=False,
                             **kwargs):
    '''Segments the CF and FM parts of a call and then
    proceeds to measure their characteristics. If required,
    also segments call from background.

    Parameters
    ----------
    main_call : np.array
    fs : float>0
        sampling rate in Hz
    segment_from_background : boolean
        Whether to segment the call in the main_call audio.
        Defaults to False.

    Keyword Arguments
    -----------------
    For further keyword arguments see segment_call_from_background,
    segment_call_into_cf_fm and measure_hbc_call

    Returns
    -------
    segmentation_outputs : tuple
        The outputs of segment_call_into_cf_fm in a tuple
    call_parts_audio : dictionary
        Dictionary with numbered entries. If a sound has the following
        order of Cf and FM: FM-CF-FM, then the keys will be
        'fm1','cf1','fm2'. The numbering is according to the chronological
        order.
    measurements : pd.DataFrame
        All the measurements from the FM and CF parts.

    Example
    -------
    Let's simulate a call to demonstrate how the measurement+segmentation
    works.

    >>> import scipy.signal  as signal
    >>> from itsfm.simulate_calls import make_cffm_call
    >>> call_properties = {'cf':(80000, 0.01), 'upfm':(70000, 0.002),
                           'downfm':(50000, 0.002)}
    >>> fs = 500000
    >>> call, profile = make_cffm_call(call_properties, fs)
    >>> call *= signal.tukey(call.size, 0.1)
    >>> plt.figure()
    >>> plot1 = plt.subplot(211)
    >>> plt.plot(profile)
    >>> #segment the CF and FM parts with the default 'peak percentage' method.
    >>> segm_out, call_parts, measures, _ = segment_and_measure_call(call,
                                                                    fs,
                                                                    segment_method='peak_percentage',
                                                                    peak_percentage=0.999,
                                                                    window_size=int(fs*0.5*10**-3))
    >>> print(measures)

    Now segment with frequency tracking implemented with the Pseudo Wigner Ville
    Distribution, and the set the fmrate threshold to 10 kHz/ms

    >>> segm_out, call_parts, measures, _ = segment_and_measure_call(call,
                                                                     fs,
                                                                     segment_method='pwvd',
                                                                     fmrate_threshold=10,
                                                                     medianfilter_length=0.5*10**-3,
                                                                     )
    >>> plt.subplot(212, sharex=plot1)
    >>> plt.plot(segm_out[-1]['fmrate'])
    >>> print(measures)

    See Also
    --------
    itsfm.measure

    '''
    cf, fm, info = segment_call_into_cf_fm(main_call, fs,
                                                           **kwargs)

    measurements = measure_hbc_call(main_call, fs, cf, fm,
                                                        **kwargs)
    call_segments = parse_cffm_segments(cf,fm)

    call_parts_audio = {}
    for (part_name, location) in call_segments:
        call_parts_audio[part_name] = main_call[location]

    return (cf, fm, info), call_parts_audio, measurements

def save_overview_graphs(all_subplots, analysis_name, file_name, index,
                          **kwargs):
    '''Saves overview graphs.

    Parameters
    ----------
    all_subplots : list
        List with plt.subplot objects in them.
        For each figure to be saved, one subplot object is enough.
    analysis_name : str
        The name of the analysis. If this funciton is called
        through a batchfile, then it becomes the name of the
        batchfile
    file_name : str
    index : int, optional
        A numeric identifier for each graph. This is especially relevant
        for analyses driven by batch files as there may be cases where the
        calls are selected from the same audio file but in different parts.

    Returns
    -------
    None

    Notes
    -----
    This function has the main side effect of saving all the input figures
    into a pdf file with >1 pages (one page per plot) for the user to inspect
    the results.

    Example
    ---------
    import numpy as np

    # 1st plot
    plt.figure()
    a = plt.subplot(211)
    plt.plot([1,2,3])
    b = plt.subplot(212)
    plt.plot([5,4,3])

    #2nd plot
    plt.figure()
    c = plt.subplot(121)
    plt.plot(np.random.normal(0,1,100))
    d = plt.subplot(122)
    plt.plot(np.random.normal(0,1,10))

    save_overview_graphs([a,c], 'example_plots', 'example_file',0)
    '''

    final_file_name = analysis_name+'_'+file_name+'_'+str(index)

    # thanks to J0e3gan : https://stackoverflow.com/a/17788764
    pdf = matplotlib.backends.backend_pdf.PdfPages(final_file_name+".pdf")

    for one_subplot in all_subplots:
        pdf.savefig(one_subplot.figure)
    pdf.close()

max_stars_repo_name: thejasvibr/itsfm

reason: import numpy,import scipy


text:
import torch
import torchvision
import torchvision.transforms as transforms
import numpy as np
import urllib.request

import os

def create_train_dataset(batch_size = 128, root = '../data'):
    """
    Create different training dataset
    """

    transform_train = transforms.Compose([
    transforms.ToTensor(),
    ])
    trainset = torchvision.datasets.MNIST(root=root, train=True, download=True, transform=transform_train)
    trainloader = torch.utils.data.DataLoader(trainset, batch_size=batch_size, shuffle=True, num_workers=2)

    return trainloader

def create_test_dataset(batch_size = 128, root = '../data'):
    transform_test = transforms.Compose([
    transforms.ToTensor(),
    ])
    testset = torchvision.datasets.MNIST(root=root, train=False, download=True, transform=transform_test)
    testloader = torch.utils.data.DataLoader(testset, batch_size=batch_size, shuffle=False, num_workers=2)
    return testloader

def download_model(url, file):
    print('Dowloading from {} to {}'.format(url, file))
    try:
        urllib.request.urlretrieve(url, file)
    except:
        raise Exception("Download failed! Make sure you have stable Internet connection and enter the right name")

def save_checkpoint(now_epoch, net, optimizer, lr_scheduler, file_name):
    checkpoint = {'epoch': now_epoch,
                  'state_dict': net.state_dict(),
                  'optimizer_state_dict': optimizer.state_dict(),
                  'lr_scheduler_state_dict':lr_scheduler.state_dict()}
    if os.path.exists(file_name):
        print('Overwriting {}'.format(file_name))
    torch.save(checkpoint, file_name)
    # link_name = os.path.join(*file_name.split(os.path.sep)[:-1], 'last.checkpoint')
    # #print(link_name)
    # make_symlink(source = file_name, link_name=link_name)

def load_checkpoint(file_name, net = None, optimizer = None, lr_scheduler = None):
    if os.path.isfile(file_name):
        print("=> loading checkpoint '{}'".format(file_name))
        check_point = torch.load(file_name)
        if net is not None:
            print('Loading network state dict')
            net.load_state_dict(check_point['state_dict'])
        if optimizer is not None:
            print('Loading optimizer state dict')
            optimizer.load_state_dict(check_point['optimizer_state_dict'])
        if lr_scheduler is not None:
            print('Loading lr_scheduler state dict')
            lr_scheduler.load_state_dict(check_point['lr_scheduler_state_dict'])

        return check_point['epoch']
    else:
        print("=> no checkpoint found at '{}'".format(file_name))

def make_symlink(source, link_name):
    """
    Note: overwriting enabled!
    """

    if os.path.exists(link_name):
        print("Link name already exist! Removing '{}' and overwriting".format(link_name))
        os.remove(link_name)
    if os.path.exists(source):
        os.symlink(source, link_name)
        return
    else:
        print('Source path not exists')

from texttable import Texttable
def tab_printer(args):
    """
    Function to print the logs in a nice tabular format.
    input:
        param args: Parameters used for the model.
    """
    args = vars(args)
    keys = sorted(args.keys())
    t = Texttable()
    t.add_rows([["Parameter", "Value"]] +  [[k.replace("_"," ").capitalize(), args[k]] for k in keys])
    print(t.draw())

def onehot_like(a, index, value=1):
    """Creates an array like a, with all values
    set to 0 except one.
    Parameters
    ----------
    a : array_like
        The returned one-hot array will have the same shape
        and dtype as this array
    index : int
        The index that should be set to `value`
    value : single value compatible with a.dtype
        The value to set at the given index
    Returns
    -------
    `numpy.ndarray`
        One-hot array with the given value at the given
        location and zeros everywhere else.
    """
    #TODO: change the note here.
    x = np.zeros_like(a)
    x[index] = value
    return x

def reduce_sum(x, keepdim=True):
    # silly PyTorch, when will you get proper reducing sums/means?
    for a in reversed(range(1, x.dim())):
        x = x.sum(a, keepdim=keepdim)
    return x

def arctanh(x, eps=1e-6):
    """
    Calculate arctanh(x)
    """
    x *= (1. - eps)
    return (np.log((1 + x) / (1 - x))) * 0.5

def l2r_dist(x, y, keepdim=True, eps=1e-8):
    d = (x - y)**2
    d = reduce_sum(d, keepdim=keepdim)
    d += eps  # to prevent infinite gradient at 0
    return d.sqrt()

def l2_dist(x, y, keepdim=True):
    d = (x - y)**2
    return reduce_sum(d, keepdim=keepdim)

def l1_dist(x, y, keepdim=True):
    d = torch.abs(x - y)
    return reduce_sum(d, keepdim=keepdim)

def l2_norm(x, keepdim=True):
    norm = reduce_sum(x*x, keepdim=keepdim)
    return norm.sqrt()

def l1_norm(x, keepdim=True):
    return reduce_sum(x.abs(), keepdim=keepdim)

def adjust_learning_rate(optimizer, epoch, learning_rate):
    """decrease the learning rate"""
    lr = learning_rate
    if epoch >= 55:
        lr = learning_rate * 0.1
    if epoch >= 75:
        lr = learning_rate * 0.01
    if epoch >= 90:
        lr = learning_rate * 0.001
    for param_group in optimizer.param_groups:
        param_group['lr'] = lr

    return optimizer

def progress_bar(current, total, msg=None):
    global last_time, begin_time
    if current == 0:
        begin_time = time.time()  # Reset for new bar.

    cur_len = int(TOTAL_BAR_LENGTH*current/total)
    rest_len = int(TOTAL_BAR_LENGTH - cur_len) - 1

    sys.stdout.write(' [')
    for i in range(cur_len):
        sys.stdout.write('=')
    sys.stdout.write('>')
    for i in range(rest_len):
        sys.stdout.write('.')
    sys.stdout.write(']')

    cur_time = time.time()
    step_time = cur_time - last_time
    last_time = cur_time
    tot_time = cur_time - begin_time

    L = []
    L.append('  Step: %s' % format_time(step_time))
    L.append(' | Tot: %s' % format_time(tot_time))
    if msg:
        L.append(' | ' + msg)

    msg = ''.join(L)
    sys.stdout.write(msg)
    for i in range(term_width-int(TOTAL_BAR_LENGTH)-len(msg)-3):
        sys.stdout.write(' ')

    # Go back to the center of the bar.
    for i in range(term_width-int(TOTAL_BAR_LENGTH/2)+2):
        sys.stdout.write('\b')
    sys.stdout.write(' %d/%d ' % (current+1, total))

    if current < total-1:
        sys.stdout.write('\r')
    else:
        sys.stdout.write('\n')
    sys.stdout.flush()

max_stars_repo_name: smaranjitghose/DeepRobust

max_issues_repo_name: qilong-zhang/DeepRobust

max_forks_repo_name: qilong-zhang/DeepRobust


text:
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
import os
import copy
import sys
from scipy.stats import norm
import time
import seaborn as sns
import pandas as pd
from baseclass import BaseClass
from mpi4py import MPI

class SMC(BaseClass):
    """Class for different sampling algorithms"""

    def __init__(self,model,qstart,qlims,nsamples,Vstart=None,output_file="samples.txt",
                     method="SMC",adapt_interval=100,printoutput=True,plotfigures=True,
                     saveoutputfile=True):
        """
        Inhereted from BaseClass
        """
        BaseClass.__init__(self,model=model,
                                qstart=qstart,
                                qlims=qlims,
                                nsamples=nsamples,
                                output_file=output_file,
                                printoutput=printoutput,
                                plotfigures=plotfigures,
                                saveoutputfile=saveoutputfile)
        if Vstart==None:
            self.Vstart=np.eye(len(qstart))*1.0/len(qstart)
        else:
            self.Vstart=Vstart
        self.adapt_interval=adapt_interval
        self.method=method

        self.ESSfract=0.5
        self.tauThres=1e6
        self.NSMC_MCMC=3
        # self.nsamples=self.nsamples

        self.comm = MPI.COMM_WORLD
        self.myrank = self.comm.Get_rank()
        self.size=self.comm.size

    def sample(self):
        """
        Function to call the sampling method of choice and plot results
        """
        if self.myrank==0:
            print("Sampling using sequential monte carlo")

        start_time=time.time()
        qsamples,ntot,naccept=self.SMCsampler()
        comp_time=time.time()-start_time

        if self.myrank==0:
            super(SMC, self).printoutput(qsamples,ntot,naccept,naccept/ntot,self.size,comp_time)
            # super(SMC, self).plotsamples(qsamples,self.method)
            super(SMC, self).savesamples(qsamples,self.method)

        return naccept,naccept/ntot

    def SMCsampler(self):
        """
        Function for sampling using sequential monter carlo

        Return:
            Q_MCMC: final samples
        """
        tau=0;
        tcount=0
        allaccepts_qs=[]

        if self.myrank==0:
            qprtls=np.random.rand(self.qlims.shape[0],self.nsamples)
            qprtls=np.reshape(self.qlims[:,0],(-1,1))+(np.reshape(self.qlims[:,1],(-1,1))-np.reshape(self.qlims[:,0],(-1,1)))*qprtls[:,:]

            Var=copy.deepcopy(self.Vstart)
            VarR=copy.deepcopy(Var)
            W=np.ones(self.nsamples)*1/self.nsamples

            pbar = tqdm(total = 1000)
            count=0
            countold=0
        else:
            qprtls=None

        start_time=time.time()

        while tau < self.tauThres:

            tcount=tcount+1
            partitions=self.distributor(qprtls)
            qprtls_div=self.comm.scatter(partitions,root=0)

            tau,wnt,ESSeff=self.find_tau(qprtls_div,tauL=tau,tauU=self.tauThres)
            wnt_all=self.comm.gather(wnt,root=0)

            if self.myrank==0:
                wnt_all=np.hstack(wnt_all)
                W=W*wnt_all;    W=W/np.sum(W)
                qprtls,W=self.resampling(qprtls,W)
                Var,VarR=self.calSMCvar(qprtls,W,Var,VarR)
            else:
                Var=None; VarR=None

            partitions=self.distributor(qprtls)
            qprtls_div=self.comm.scatter(partitions,root=0)
            Var=self.comm.bcast(Var,root=0)
            VarR=self.comm.bcast(VarR,root=0)

            iacceptall=0
            if tau>=self.tauThres:
                for i in range(0,qprtls_div.shape[1]):
                    # qprtls[:,i],iaccept,allaccepts_dum=self.SMCRW(tau,qprtls[:,i],Var,VarR)
                    qprtls_div[:,i],iaccept,allaccepts_dum=self.GibbsSampler(tau,qprtls_div[:,i],Var,VarR)
                    iacceptall=iacceptall+iaccept
                    allaccepts_qs=allaccepts_qs+allaccepts_dum
                allaccepts_qs=np.asarray(allaccepts_qs)
                allaccepts_qs=self.comm.gather(allaccepts_qs,root=0)

            else:
                for i in range(0,qprtls_div.shape[1]):
                    # qprtls[:,i],iaccept,allaccepts_dum=self.SMCRW(tau,qprtls[:,i],Var,VarR)
                    qprtls_div[:,i],iaccept,_=self.GibbsSampler(tau,qprtls_div[:,i],Var,VarR)
                    iacceptall=iacceptall+iaccept

            qprtls=self.comm.gather(qprtls_div,root=0)
            if self.myrank==0:
                qprtls=np.hstack(qprtls)
                count=np.abs(int(np.log10(tau)*1000/6))
                pbar.update(count-countold)
                countold=copy.deepcopy(count)

            if tcount>200:
                print("Number of iterations are more than ",tcount)

                for i in range(0,qprtls_div.shape[1]):
                    # qprtls[:,i],iaccept,allaccepts_dum=self.SMCRW(tau,qprtls[:,i],Var,VarR)
                    qprtls_div[:,i],iaccept,allaccepts_dum=self.GibbsSampler(tau,qprtls_div[:,i],Var,VarR)
                    iacceptall=iacceptall+iaccept
                    allaccepts_qs=allaccepts_qs+allaccepts_dum
                allaccepts_qs=np.asarray(allaccepts_qs)
                allaccepts_qs=self.comm.gather(allaccepts_qs,root=0)

                break

        if self.myrank==0:
            allaccepts_qs=np.vstack(allaccepts_qs)
            allaccepts_qs=np.unique(allaccepts_qs,axis=0)
            EFF_prtls=np.unique(allaccepts_qs,axis=0).shape[0]
            print("Accepted samples:",allaccepts_qs.shape)
        else:
            EFF_prtls=1

        ntot=(self.nsamples*self.NSMC_MCMC)*tcount
        return allaccepts_qs,ntot,EFF_prtls

    def distributor(self,qprtls):

        if self.myrank==0:
            partitions=[]
            for prtls in np.array_split(np.arange(0,self.nsamples),self.size,axis=0):
                partitions.append(qprtls[:,prtls])
        else:
            partitions=[]

        return partitions

    def find_tau(self,qprtls,tauL,tauU,deltatau=100.0,Niter=20):
        """
        Function to find tau using ESS

        Return:
            tau: returns
        """

# #------------------------manual tau
#         pitallminus=self.pitallfun(tauL,qprtls)
#         taunew=10**(np.log10(tauL+1)+0.1/2)
#         wnt=self.ESSfun(taunew,qprtls,pitallminus)
#         wnt_all=self.comm.gather(wnt,root=0)
#         if self.myrank==0:
#             wnt_all=np.hstack(wnt_all)
#             ESS=(np.sum(wnt_all))**2/(np.sum(wnt_all**2)+1e-16)
#         else:
#             ESS=None
#         ESS=self.comm.bcast(ESS,root=0)
#         return taunew,wnt,ESS/self.nsamples

#-------------------------------------------

        pitallminus=self.pitallfun(tauL,qprtls)
        funL=1.0-self.ESSfract
        wnt=self.ESSfun(tauU,qprtls,pitallminus)

        wnt_all=self.comm.gather(wnt,root=0)
        if self.myrank==0:
            wnt_all=np.hstack(wnt_all)
            ESS=(np.sum(wnt_all))**2/(np.sum(wnt_all**2)+1e-16)
        else:
            ESS=None
        ESS=self.comm.bcast(ESS,root=0)

        funU=ESS/self.nsamples-self.ESSfract
        fact=1.0

        if funU>=0.0:
            return tauU,wnt,ESS/self.nsamples

        for iter in np.arange(0,Niter):

            taunew=tauL+fact*deltatau
            wnt=self.ESSfun(taunew,qprtls,pitallminus)
            wnt_all=self.comm.gather(wnt,root=0)
            if self.myrank==0:
                wnt_all=np.hstack(wnt_all)
                ESS=(np.sum(wnt_all))**2/(np.sum(wnt_all**2)+1e-16)
            else:
                ESS=None
            ESS=self.comm.bcast(ESS,root=0)

            funnew=ESS/self.nsamples-self.ESSfract
            if funnew<0:
                fact=fact/2.0
                tauU=copy.deepcopy(taunew)
                funU=copy.deepcopy(funnew)
            elif funnew>0:
                fact=fact*2.0
                tauL=copy.deepcopy(taunew)
                funL=copy.deepcopy(funnew)

            if np.abs(funnew)<1e-3:
                break

        return taunew,wnt,ESS/self.nsamples

    def ESSfun(self,tau,qprtls,pitallminus):
        """
        Calcualtes effective samples

        Return:
            ESS: returns
        """
        pitall=self.pitallfun(tau,qprtls)
        wnt=pitall/(pitallminus+1e-16)
        # ESS=(np.sum(wnt))**2/(np.sum(wnt**2)+1e-16)
        # return ESS,wnt
        return wnt

    def pitallfun(self,tau,qprtls):
        """
        Function for evaluatinng the target distribution for sequential monte carlo
        𝜋(𝑥) = 𝜙(𝜏𝐶1(𝑥)) * 𝜙(𝜏𝐶2(𝑥)) * .... * 𝜙(𝜏𝐶k(𝑥))
        for all particles

        Return:
            pit: returns
        """
        pitall=[]
        for i in range(0,qprtls.shape[1]):
            pitall.append(self.pitfun(tau,qprtls[:,i]))
        pitall=np.asarray(pitall)

        return pitall

    def pitfun(self,tau,qval):
        """
        Function for evaluatinng the target distribution for sequential monte carlo
        𝜋(𝑥) = 𝜙(𝜏𝐶1(𝑥)) * 𝜙(𝜏𝐶2(𝑥)) * .... * 𝜙(𝜏𝐶k(𝑥))
        for single particle

        Return:
            pit: returns
        """
        Call= self.model.apply_eval(list(qval))
        Call=Call+list(qval-self.qlims[:,0])
        Call=Call+list(-qval+self.qlims[:,1])

        pit=tau*np.asarray(Call)
        pit=norm.cdf(pit)
        pit=np.prod(pit)

        return pit

    def resampling(self,qprtls,W):
        """
        Resampling using systematic resampling

        Return:
            pit: returns
        """
        qprtlsnew=copy.deepcopy(qprtls)

        j=0
        isavejend=0
        sumW=W[j]
        u=np.random.rand()/self.nsamples

        for i in range(0,self.nsamples):
            while sumW<u:
                j=j+1
                isavejend=copy.deepcopy(i)
                sumW=sumW+W[j]
            qprtlsnew[:,i]=qprtls[:,j]
            u=u+1.0/self.nsamples
        W=np.ones(self.nsamples)/self.nsamples
        return qprtlsnew,W

    def calSMCvar(self,qprtls,W,Varold,VarRold):

        qmean=np.reshape(np.mean(qprtls,axis=1),(-1,1))
        Var=(qprtls-qmean)
        Var=np.dot(Var,Var.T)/qprtls.shape[1]
        try:
            VarR = np.linalg.cholesky(Var)
        except:
            VarR=copy.deepcopy(VarRold)
            Var=copy.deepcopy(Varold)

        return Var,VarR

    def SMCRW(self,tau,qval,Var,VarR):
        """
        Function for sampling using random walk Metropolis algorithm for SMC

        Return:
            Q_MCMC: Accepted samples
        """
        qval=np.reshape(np.asarray(qval),(-1,1))
        qmean=copy.deepcopy(qval)
        qmean_old=copy.deepcopy(qval)
        Var=copy.deepcopy(Var)
        Varold=copy.deepcopy(Var)

        iaccept=0
        output0=self.pitfun(tau,qval[:,0])+1e-16
        accepted_qs=[qval[:,0]]
        for isamples in range(0,self.NSMC_MCMC):

            z = np.random.randn(qval.shape[0],1)
            qnew = qval + np.dot(VarR.T,z)

            output=self.pitfun(tau,qnew[:,0])
            ratio=output/output0

            if ratio>np.random.rand():
                qval=copy.deepcopy(qnew)
                output0=copy.deepcopy(output0)
                iaccept=iaccept+1
                accepted_qs.append(list(qval[:,0]))

        return qval[:,0],iaccept,accepted_qs

    def GibbsSampler(self,tau,qval,Var,VarR):
        """
        Function for sampling using gibs sampler for SMC

        Return:
            qval: Accepted samples
        """
        qval=np.reshape(np.asarray(qval),(-1,1))

        iaccept=0
        output0=self.pitfun(tau,qval[:,0])+1e-16
        accepted_qs=[qval[:,0]]

        qindG=np.random.randint(qval.shape[0],size=self.NSMC_MCMC)

        for ind in qindG:
            # print(ind)
            qnew=copy.deepcopy(qval)
            qnew[ind,0]=qval[ind,0]+np.sqrt(Var[ind,ind])*np.random.randn()

            output=self.pitfun(tau,qnew[:,0])
            ratio=output/output0

            if ratio>np.random.rand():
                qval=copy.deepcopy(qnew)
                output0=copy.deepcopy(output0)
                iaccept=iaccept+1
                accepted_qs.append(list(qval[:,0]))

        return qval[:,0],iaccept,accepted_qs

reason: import numpy,from scipy


text:
from __future__ import division
import refnx
from refnx.reflect import structure, ReflectModel, SLD
from refnx.dataset import ReflectDataset
from refnx.analysis import (
    Transform,
    CurveFitter,
    Objective,
    GlobalObjective,
    Parameter,
)

import numpy as np
import scipy
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns

sns.set(palette="colorblind")
import corner

import sys

sys.path.insert(0, "../reports/code_blocks")
import mol_vol as mv
import ref_help as rh

import matplotlib as mpl

mpl.rcParams["xtick.labelsize"] = 8
mpl.rcParams["ytick.labelsize"] = 8
mpl.rcParams["axes.facecolor"] = "w"
mpl.rcParams["lines.linewidth"] = 2
mpl.rcParams["xtick.top"] = False
mpl.rcParams["xtick.bottom"] = True
mpl.rcParams["ytick.left"] = True
mpl.rcParams["grid.linestyle"] = "--"
mpl.rcParams["legend.fontsize"] = 8
mpl.rcParams["legend.facecolor"] = [1, 1, 1]
mpl.rcParams["legend.framealpha"] = 0.75
mpl.rcParams["axes.labelsize"] = 8
mpl.rcParams["axes.linewidth"] = 1
mpl.rcParams["axes.edgecolor"] = "k"

# The lipid to be investigated
lipid = sys.argv[1]
# Number of carbon atoms in tail group
t_length = int(sys.argv[2])
# The type of radiation that is used
neutron = int(sys.argv[7])

# The surface pressures to probe
if lipid == 'dspc':
    sp = sys.argv[3]
    cont = ['d13acmw', 'd13d2o', 'hd2o', 'd70acmw',
            'd70d2o', 'd83acmw', 'd83d2o']
    apm = sys.argv[4]
    neutron = 1
else:
    if neutron:
        sp = sys.argv[3]
        cont = [sys.argv[4], sys.argv[5]]
    else:
        sp = np.array(
            [sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6]]
        )

# A label for the output files
if lipid == 'dspc':
    label = "{}_{}".format(lipid, sp)
else:
    if neutron:
        radiation = "neutron"
        label = "{}_{}_{}".format(lipid, radiation, sp)

    else:
        radiation = "xray"
        label = "{}_{}".format(lipid, radiation)

# Relative directory locations
if lipid == 'dspc':
    data_dir = "../data/reflectometry2/{}".format(label)
    fig_dir = "../reports/figures/reflectometry2/{}".format(label)
    anal_dir = "../output/reflectometry2/{}".format(label)
else:
    if neutron:
        data_dir = "../data/reflectometry1/{}".format(label[:-3])
        fig_dir = "../reports/figures/reflectometry1/{}".format(
            label[:-3]
        )
        anal_dir = "../output/reflectometry1/{}".format(label[:-3])
        xanal_dir = "../output/reflectometry1/{}_xray".format(lipid)
    else:
        data_dir = "../data/reflectometry1/{}".format(label)
        fig_dir = "../reports/figures/reflectometry1/{}".format(label)
        anal_dir = "../output/reflectometry1/{}".format(label)

# For the analysis to be reproduced exactly, the same versions of `refnx`, `scipy`, and `numpy`, must be present.
# The versions used in the original version are:
#
# refnx.version.full_version == 0.1.2
# scipy.version.version == 1.1.0
# np.__version__ == 1.15.4

refnx.version.full_version, scipy.version.version, np.__version__

datasets = []
if lipid == 'dspc':
    for i in range(len(cont)):
        datasets.append(
            ReflectDataset(
                "{}/{}{}.dat".format(
                    data_dir, cont[i], sp
                )
            )
        )
else:
    if neutron:
        for i in range(len(cont)):
            datasets.append(
                ReflectDataset(
                    "{}/{}_{}_sp_{}.dat".format(
                        data_dir, label[:-3], cont[i], sp
                    )
                )
            )
            datasets[i].mask = datasets[i].x <= 0.6
    else:
        for i in range(sp.size):
            datasets.append(
                ReflectDataset(
                    "{}/{}_sp_{}.dat".format(data_dir, label, sp[i])
                )
            )
            datasets[i].mask = datasets[i].x <= 0.6

if lipid == 'dspc':
    b_head = []
    b_tail = []
    for i in range(len(cont)):
        if cont[i][:3] == 'd13':
            head = {"C": 10, "D": 13, "H": 5, "O": 8, "N": 1, "P": 1}
            tail = {"C": t_length * 2, "H": t_length * 4 + 2}
        elif cont[i][:3] == 'd70':
            head = {"C": 10, "D": 5, "H": 13, "O": 8, "N": 1, "P": 1}
            tail = {"C": t_length * 2, "D": t_length * 4 + 2}
        elif cont[i][:3] == 'd83':
            head = {"C": 10, "D": 18, "O": 8, "N": 1, "P": 1}
            tail = {"C": t_length * 2, "D": t_length * 4 + 2}
        elif cont[i][:1] == 'h':
            head = {"C": 10, "H": 18, "O": 8, "N": 1, "P": 1}
            tail = {"C": t_length * 2, "H": t_length * 4 + 2}
        b_head.append(rh.get_scattering_length(head, 1))
        b_tail.append(rh.get_scattering_length(tail, 1))
else:
    if lipid == "dmpg":
        head = {"C": 8, "H": 12, "O": 10, "Na": 1, "P": 1}
    else:
        head = {"C": 10, "H": 18, "O": 8, "N": 1, "P": 1}
    if neutron:
        tail = {"C": t_length * 2, "D": t_length * 4 + 2}
    else:
        tail = {"C": t_length * 2, "H": t_length * 4 + 2}
    b_head = rh.get_scattering_length(head, neutron)
    b_tail = rh.get_scattering_length(tail, neutron)

def get_value(file):
    f = open(file, "r")
    for line in f:
        k = line
    if "^" in k:
        l = k.split("$")[1].split("^")[0]
    else:
        l = k.split("$")[1].split("\\pm")[0]
    return float(l)

if lipid == 'dspc':
    d_h = 8.5
    V_h = 330.0
    V_t = 1100.0
    min_d_t = 9
else:
    if neutron:
        d_h = get_value(
            "{}/{}.tex".format(xanal_dir, "{}_xray-d_h".format(lipid))
        )
        V_h = get_value(
            "{}/{}.tex".format(xanal_dir, "{}_xray-V_h".format(lipid))
        )
        if lipid[:2] == "dl":
            V_t = get_value(
                "{}/{}.tex".format(xanal_dir, "{}_xray-V_t".format(lipid))
            )
            min_d_t = 5
        elif lipid[:2] == "dm":
            V_t = get_value(
                "{}/{}.tex".format(xanal_dir, "{}_xray-V_t".format(lipid))
            )
            min_d_t = 7
        elif lipid[:2] == "dp":
            V_t = get_value(
                "{}/{}.tex".format(xanal_dir, "{}_xray-V_t".format(lipid))
            )
            min_d_t = 9
    else:
        d_h = 12.0
        V_h = 330.0
        if lipid[:2] == "dl":
            V_t = 667.0
            min_d_t = 5
        elif lipid[:2] == "dm":
            V_t = 779.0
            if lipid[2:] == "pg":
                min_d_t = 5
            else:
                min_d_t = 7
        elif lipid[:2] == "dp":
            V_t = 891.0
            min_d_t = 9

if lipid == 'dspc':
    lipids = []
    for i in range(len(cont)):
        lipids.append(
            mv.VolMono(
                [V_h, V_t], [b_head[i], b_tail[i]], d_h, t_length, name=label
            )
        )
else:
    if neutron:
        t = len(cont)
    else:
        t = sp.size
    lipids = []
    for i in range(t):
        lipids.append(
            mv.VolMono(
                [V_h, V_t], [b_head, b_tail], d_h, t_length, name=label
            )
        )

air = SLD(0, "air")
if lipid == 'dspc':
    structures = []
    for i in range(len(cont)):
        if cont[i][-3:] == 'd2o':
            water = SLD(6.35, "d2o")
        elif cont[i][-4:] == 'acmw':
            water = SLD(0.0, "acmw")
        structures.append(air(0, 0) | lipids[i] | water(0, 3.3))
else:
    if neutron:
        des_1 = SLD(0.43, "hdes")
        des_2 = SLD(3.15, "hddes")

        structures = []
        structure_lipid_1 = air(0, 0) | lipids[0] | des_1(0, 0)
        structure_lipid_2 = air(0, 0) | lipids[1] | des_2(0, 0)
        structures.append(structure_lipid_1)
        structures.append(structure_lipid_2)

    else:
        des = SLD(10.8, "des")

        structures = []
        for i in range(sp.size):
            structures.append(air(0, 0) | lipids[i] | des(0, 3.3))

if lipid == 'dspc':
    for i in range(len(cont)):
        lipids[i].vol[0].setp(
            vary=True, bounds=(V_h * 0.8, V_h * 1.2)
        )
        lipids[i].vol[1].setp(
            vary=True, bounds=(V_t * 0.8, V_t * 1.2)
        )
        lipids[i].d[0].setp(vary=True, bounds=(5, 15))
        max_d_t = 1.54 + 1.265 * t_length
        lipids[i].d[1].setp(vary=True, bounds=(min_d_t, max_d_t))
        lipids[i].vol[1].constraint = lipids[i].d[1] * float(apm)
        structures[i][-1].rough.setp(vary=True, bounds=(3.3, 6))
else:
    if neutron:
        t = len(cont)
        for i in range(t):
            lipids[i].vol[0].setp(vary=False)
            lipids[i].vol[1].setp(vary=False)
            lipids[i].d[0].setp(vary=False)
            max_d_t = 1.54 + 1.265 * t_length
            lipids[i].d[1].setp(vary=True, bounds=(min_d_t, max_d_t))
            structures[i][-1].rough.setp(vary=True, bounds=(3.3, 6))
    else:
        t = sp.size
        for i in range(t):
            lipids[i].vol[0].setp(
                vary=True, bounds=(V_h * 0.8, V_h * 1.2)
            )
            lipids[i].vol[1].setp(
                vary=True, bounds=(V_t * 0.8, V_t * 1.2)
            )
            lipids[i].d[0].setp(vary=True, bounds=(10, 17))
            max_d_t = 1.54 + 1.265 * t_length
            lipids[i].d[1].setp(vary=True, bounds=(min_d_t, max_d_t))
            structures[i][-1].rough.setp(vary=True, bounds=(3.3, 6))

if neutron:
    lipids = rh.set_constraints(
        lipids,
        structures,
        hold_tails=True,
        hold_rough=True,
        hold_phih=True,
    )
else:
    lipids = rh.set_constraints(lipids, structures, hold_rough=False)

models = []
if neutron:
    t = len(cont)
else:
    t = sp.size

for i in range(t):
    models.append(ReflectModel(structures[i]))
    models[i].scale.setp(vary=True, bounds=(0.005, 10))
    models[i].bkg.setp(datasets[i].y[-1], vary=False)

objectives = []
if neutron:
    t = len(cont)
else:
    t = sp.size
for i in range(t):
    objectives.append(
        Objective(models[i], datasets[i], transform=Transform("YX4"))
    )

global_objective = GlobalObjective(objectives)

fitter = CurveFitter(global_objective)
np.random.seed(1)
res = fitter.fit("differential_evolution")

print(global_objective)

fitter.sample(200, random_state=1)
fitter.sampler.reset()
res = fitter.sample(
    1000,
    nthin=1,
    random_state=1,
    f="{}/{}_chain.txt".format(anal_dir, label),
)


text:
import pyglet
import random
from time import sleep
import pyaudio
import numpy as np

RATE    = 16000
CHUNK   = 256

p = pyaudio.PyAudio()

player = p.open(format=pyaudio.paInt16, channels=1, rate=RATE, output=True,
frames_per_buffer=CHUNK)
stream = p.open(format=pyaudio.paInt16, channels=1, rate=RATE, input=True, frames_per_buffer=CHUNK)
maxx = 1000

Width = 1000
Height = 500
window = pyglet.window.Window(Width, Height)
batch = pyglet.graphics.Batch()
snow = pyglet.image.load("snow.png").texture
x_offset = 5
y_offset = 2

class Graph:
    def __init__(self):
        self.size = 256
        self.points = []
        self.objects = []

    def clear(self):
        self.points = []
        self.objects = []

    def set_size(self, size):
        self.size = size

    def generate(self):
        self.objects = []
        data = np.fromstring(stream.read(CHUNK), dtype=np.int16)
        for x in range(0, self.size):
            y = data[x]/300
            self.objects.append(pyglet.sprite.Sprite(snow, batch=batch, x=x*x_offset, y=y*y_offset+250))

    def next_point(self):
        self.points.pop(0)
        self.objects.pop(0)
        temp_list = []
        for point in self.points:
            temp_list.append((point[0]-1, point[1]))
        for obj in self.objects:
            obj.x -= x_offset
        self.points = temp_list
        self.pointer.next_point()
        y = self.pointer.get_point()
        self.points.append((self.size-1, y))

        self.objects.append(pyglet.sprite.Sprite(snow, batch=batch, x=self.size*x_offset-x_offset, y=y*y_offset+250))

    def get_points(self):
        return self.points

class Display:
    def __init__(self):
        self.graph = Graph()
        self.color = (255, 255, 255)

    def draw(self):
        #for point in self.graph.get_points():
        #    draw_pixel(point, self.color)
        batch.draw()

    def next_graph_point(self):
        self.graph.generate()

    def new_graph(self, size):
        self.graph.clear()
        self.graph.generate()

def draw_pixel(cord=(0, 0), color=(0, 0, 0)):
    x = round(cord[0] * 5)
    y = round(cord[1]) + int(Height/2)
    pyglet.graphics.draw(1, pyglet.gl.GL_POINTS, ('v2i', (x, y)), ('c3B', (color[0], color[1], color[2])))

displayable = []
first_run = True

def update(dt):
    global first_run
    if first_run:
        display = Display()
        display.new_graph(200)
        displayable.append(display)
        first_run = False
    if len(displayable) != 0:
        for item in displayable:
            item.next_graph_point()
    on_draw()

def on_draw():
    window.clear()
    if len(displayable) != 0:
        for item in displayable:
            item.draw()

pyglet.clock.schedule_interval(update, 1/120)
pyglet.app.run()


text:
# ## Copyright (c) 2014, Zhenwen Dai
# Licensed under the BSD 3-clause license (see LICENSE.txt)
from __future__ import print_function

import numpy as np
import sys

try:
    #In Python 2, cPickle is faster. It does not exist in Python 3 but the underlying code is always used
    #if available
    import cPickle as pickle
except ImportError:
    import pickle

class Metropolis_Hastings(object):
    def __init__(self,model,cov=None):
        """Metropolis Hastings, with tunings according to Gelman et al. """
        self.model = model
        current = self.model.optimizer_array
        self.D = current.size
        self.chains = []
        if cov is None:
            self.cov = np.eye(self.D)
        else:
            self.cov = cov
        self.scale = 2.4/np.sqrt(self.D)
        self.new_chain(current)

    def new_chain(self, start=None):
        self.chains.append([])
        if start is None:
            self.model.randomize()
        else:
            self.model.optimizer_array = start

    def sample(self, Ntotal=10000, Nburn=1000, Nthin=10, tune=True, tune_throughout=False, tune_interval=400):
        current = self.model.optimizer_array
        fcurrent = self.model.log_likelihood() + self.model.log_prior()
        accepted = np.zeros(Ntotal,dtype=np.bool)
        for it in range(Ntotal):
            print("sample %d of %d\r"%(it+1,Ntotal),end="")
            sys.stdout.flush()
            prop = np.random.multivariate_normal(current, self.cov*self.scale*self.scale)
            self.model.optimizer_array = prop
            fprop = self.model.log_likelihood() + self.model.log_prior()

            if fprop>fcurrent:#sample accepted, going 'uphill'
                accepted[it] = True
                current = prop
                fcurrent = fprop
            else:
                u = np.random.rand()
                if np.exp(fprop-fcurrent)>u:#sample accepted downhill
                    accepted[it] = True
                    current = prop
                    fcurrent = fprop

            #store current value
            if (it > Nburn) & ((it%Nthin)==0):
                self.chains[-1].append(current)

            #tuning!
            if it & ((it%tune_interval)==0) & tune & ((it<Nburn) | (tune_throughout)):
                pc = np.mean(accepted[it-tune_interval:it])
                self.cov = np.cov(np.vstack(self.chains[-1][-tune_interval:]).T)
                if pc > .25:
                    self.scale *= 1.1
                if pc < .15:
                    self.scale /= 1.1

    def predict(self,function,args):
        """Make a prediction for the function, to which we will pass the additional arguments"""
        param = self.model.param_array
        fs = []
        for p in self.chain:
            self.model.param_array = p
            fs.append(function(*args))
        # reset model to starting state
        self.model.param_array = param
        return fs

max_stars_repo_name: ekalosak/GPy

max_issues_repo_name: ekalosak/GPy

max_forks_repo_name: ekalosak/GPy


text:
import logging
import pandas as pd
import numpy as np
import time
from scipy.optimize import linear_sum_assignment
from sklearn.preprocessing import OrdinalEncoder
from collections import Counter
from sklearn import metrics

__all__ = ['general_statistics', 'cluster_event_match', 'mcminn_eval', 'cluster_acc']

def vizualize(vectors, data):
    labels = data.label.unique()
    print(len(labels))
    inter_dist = []
    intra_dist = []
    max_dist = []
    # avg_dist = np.zeros((labels.size, labels.size))
    # max_dist = np.zeros((labels.size, labels.size))

    for i, ilabel in enumerate(labels):
        t0 = time.time()
        for j, jlabel in enumerate(labels):
            if i <= j:
                pairwise_distance = metrics.pairwise_distances(
                    vectors[(data.label == ilabel)],
                    vectors[(data.label == jlabel)],
                    metric="cosine"
                )
                mean_pairwise_distance = pairwise_distance.mean()
                max_pairwise_distance = pairwise_distance.max()
                # avg_dist[i, j] = mean_pairwise_distance
                # max_dist[i, j] = max_pairwise_distance
                max_dist.append(max_pairwise_distance)
                if i < j:
                    inter_dist.append(mean_pairwise_distance)
                elif i == j:
                    # max_pairwise_distance = metrics.pairwise_distances(
                    #     vectors[(data.label == ilabel).values],
                    #     vectors[(data.label == jlabel).values],
                    #     metric=metric
                    # ).max()
                    intra_dist.append(mean_pairwise_distance)
                    logging.info("mean intra_distance event {}: {}".format(i, mean_pairwise_distance))
                    logging.info("max intra_distance event {}: {}".format(i, max_pairwise_distance))
    logging.info("max inter distance: {}".format(np.array(max_dist).max()))
    logging.info("mean inter distance: {}".format(np.array(inter_dist).mean()))
    logging.info("mean intra distance: {}".format(np.array(intra_dist).mean()))

def cluster_acc(data, pred):
    """
    Calculate clustering accuracy.
    (Taken from https://github.com/XifengGuo/IDEC-toy/blob/master/DEC.py)
    # Arguments
        data: pd.DataFrame with shape `(n_samples,n_columns)` with a "label" column containing true labels
        pred: predicted labels, numpy.array with shape `(n_samples,)`
    # Return
        accuracy, in [0,1]
    """
    data["pred"] = pd.Series(pred, dtype=data.label.dtype)
    data = data[data.label.notna()].as_matrix(columns=["label", "pred"])
    enc = OrdinalEncoder()
    data = enc.fit_transform(data)
    y_true = data[:,0].astype(np.int64)
    y_pred = data[:,1].astype(np.int64)
    assert y_pred.size == y_true.size
    D = max(y_pred.max(), y_true.max()) + 1
    w = np.zeros((D, D), dtype=np.int64)
    for i in range(y_pred.size):
        w[y_pred[i], y_true[i]] += 1
    ind = linear_sum_assignment(w.max() - w) # Optimal label mapping based on the Hungarian algorithm

    return sum([w[i, j] for i, j in ind]) * 1.0 / y_pred.size

def general_statistics(pred):
    s = pd.Series(pred)
    stats = s.groupby(by=s.values).size().describe(percentiles=[])
    return dict(stats)

def cluster_event_match(data, pred):
    data["pred"] = pd.Series(pred, dtype=data.label.dtype)
    data = data[data.label.notna()]
    logging.info("{} labels, {} preds".format(len(data.label.unique()), len(data.pred.unique())))
    t0 = time.time()

    match = data.groupby(["label", "pred"], sort=False).size().reset_index(name="a")
    b, c = [], []
    for idx, row in match.iterrows():
        b_ = ((data["label"] != row["label"]) & (data["pred"] == row["pred"]))
        b.append(b_.sum())
        c_ = ((data["label"] == row["label"]) & (data["pred"] != row["pred"]))
        c.append(c_.sum())
    logging.info("match clusters with events took: {} seconds".format(time.time() - t0))
    match["b"] = pd.Series(b)
    match["c"] = pd.Series(c)
    # recall = nb true positive / (nb true positive + nb false negative)
    match["r"] = match["a"] / (match["a"] + match["c"])
    # precision = nb true positive / (nb true positive + nb false positive)
    match["p"] = match["a"] / (match["a"] + match["b"])
    match["f1"] = 2 * match["r"] * match["p"] / (match["r"] + match["p"])
    match = match.sort_values("f1", ascending=False)
    macro_average_f1 = match.drop_duplicates("label").f1.mean()
    macro_average_precision = match.drop_duplicates("label").p.mean()
    macro_average_recall = match.drop_duplicates("label").r.mean()
    return macro_average_precision, macro_average_recall, macro_average_f1

def mcminn_eval(data, pred, nb_tweets=5, share_tweets=0.8):
    data["pred"] = pd.Series(pred, dtype=data.label.dtype)
    count_label = data.label.value_counts()
    count_pred = data.pred.value_counts()
    data = data[data.label.isin(count_label[count_label >= nb_tweets].index)]
    labels = data.label.unique().tolist()
    all_labels = len(labels)
    data = data[data.pred != -1]
    data = data[data.pred.isin(count_pred[count_pred >= nb_tweets].index)]
    match = data.groupby(["pred", "label"], sort=False).size().reset_index(name="a")
    a, b = 0, 0
    for p in data.pred.unique():
        candidates = match[match.pred == p]
        main = candidates.a.max()
        total = candidates.a.sum()
        if main/total >= share_tweets:
            a += 1
            try:
                labels.remove(candidates[candidates.a == main].label.values[0])
            except ValueError:
                pass
        else:
            b += 1
    precision = a / (a + b)
    recall = (all_labels - len(labels))/all_labels
    f1 = 2*precision*recall/(precision+recall)
    return precision, recall, f1

max_stars_repo_name: Yomguithereal/twembeddings

max_issues_repo_name: Yomguithereal/twembeddings

max_forks_repo_name: Yomguithereal/twembeddings


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

import trisicell as tsc

def _colors(alist):
    return mcolors.LinearSegmentedColormap.from_list("", alist)

def heatmap(
    adata, color_attrs=None, layer="X", figsize=(12, 7), vmin=0, vmax=2, rvb=None
):
    """Plot HeatMap.

    For flips: rvb=["#FFFFFF", "#D92347", "#A9D0F5"]
    For CNV: rvb="RdBu_r", vmin=-0.5, vmax=0.5

    Parameters
    ----------
    adata : [type]
        [description]
    color_attrs : [type], optional
        [description], by default None
    layer : str, optional
        [description], by default "X"
    figsize : tuple, optional
        [description], by default (12, 7)
    vmin : [type], optional
        [description], by default None
    vmax : [type], optional
        [description], by default None
    rvb : [type], optional
        [description], by default None
    """

    if color_attrs is not None:
        row_colors = []
        if isinstance(color_attrs, list):
            for attr in color_attrs:
                row_colors.append(adata.obs[attr])
        else:
            row_colors.append(adata.obs[color_attrs])
    else:
        row_colors = None

    if layer == "X":
        rvb = _colors(["#A9D0F5", "#000000", "#FFFFFF"])
        adatac = adata[:, adata.var.sort_values(["CHROM", "POS"]).index].copy()
        df = adatac.to_df().copy()
        df[df == 3] = 2
        df.index.name = "cells"
        df.columns.name = "mutations"
        chromosoms = {}
        for i in list(range(1, 23, 2)) + ["Y"]:
            chromosoms[f"chr{i}"] = "#969696"
        for i in list(range(2, 23, 2)) + ["X"]:
            chromosoms[f"chr{i}"] = "#252525"
        adatac.var["chrom_color"] = adata.var["CHROM"].map(chromosoms)
        tsc.logg.info(adatac.var.CHROM.value_counts().sort_index().to_frame())
        adatac.var["chrom_color"].name = ""
        column_colors = adatac.var["chrom_color"]
    else:
        df = adata.obsm[layer].copy()
        if isinstance(rvb, list):
            rvb = _colors(rvb)
        column_colors = None

    sns.clustermap(
        df,
        vmin=vmin,
        vmax=vmax,
        metric="euclidean",
        cmap=rvb,
        row_cluster=False,
        col_cluster=False,
        row_colors=row_colors,
        col_colors=column_colors,
        cbar_pos=None,
        figsize=figsize,
        xticklabels=False,
        yticklabels=False,
        colors_ratio=(0.02, 0.02),
        dendrogram_ratio=0,
    )
    # plt.savefig(filepath, bbox_inches="tight", pad_inches=0)

def plot_dist(adata, attr):
    x = adata.obsp[attr][np.triu_indices(adata.shape[0], 1)]
    plt.hist(x, bins=50)

max_stars_repo_name: faridrashidi/trisicell

max_issues_repo_name: faridrashidi/trisicell

max_forks_repo_name: faridrashidi/trisicell


text:
# Copyright 2020 Alibaba Group Holding Limited. 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.
# =============================================================================
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import sys
import uuid
import numpy as np
import warnings

from graphlearn import pywrap_graphlearn as pywrap
from graphlearn.python.data.values import Nodes
from graphlearn.python.gsl.dag_edge import get_dag_edge, get_eid
from graphlearn.python.utils import strategy2op

class DagNode(object):
  def __init__(self, dag, op_name="", params={}):
    self._dag = dag
    self._op_name = op_name

    # Lazy init.
    self._alias = None
    self._params = {}  # keep params for the node
    self._output_field = None
    self._shape = None
    self._sparse = False
    self._path = None  # Upstream edge type
    self._type = None  # Edge type or node type of current DagNode
    self._strategy = "random"

    # Add by GSL.
    # All the DagEdges expect link to sink node are mapped to DagEdgeDef.
    self._out_edges = []  # downstream DagEdges
    self._in_edges = []  # upstream DagEdges

    # Indicates setting which field to src_input of the downstream edges
    self._graph = self._dag.graph
    self._decoder = None

    self._lookup_node = None  # Each traverse node has a LookupDagNode
    self._degree_nodes = []  # Traverse node may have several DegreeDagNodes

    # Init until Dag is ready.
    self._node_def = None
    self._nid = -1

    self._params.update(params)
    # Downstream Traverse DagNodes
    self._pos_downstreams = []
    # Downstream Negative Sampled Traverse DagNodes
    self._neg_downstreams = []

  def _add_degree_node(self, edge_type, node_from):
    edge = self._new_edge(dst_input=pywrap.kNodeIds)
    # Add an edge from upstream node to degree node.
    self._add_out_edge(edge)
    node =  DegreeDagNode(
      "GetDegree", self, [edge],
      {pywrap.kEdgeType: edge_type,
       pywrap.kNodeFrom: int(node_from),
       pywrap.kPartitionKey: pywrap.kNodeIds})
    self._degree_nodes.append(node)

  @property
  def spec(self):
    return self._decoder

  @property
  def nid(self):
    return self._nid

  @property
  def op_name(self):
    return self._op_name

  @property
  def in_edges(self):
    return self._in_edges

  @property
  def out_edges(self):
    return self._out_edges

  @property
  def pos_downstreams(self):
    return self._pos_downstreams

  @property
  def neg_downstreams(self):
    return self._neg_downstreams

  def set_path(self, path, node_from):
    assert isinstance(path, str)
    assert isinstance(node_from, pywrap.NodeFrom)
    topo = self._graph.get_topology()
    path_to_type = {pywrap.NodeFrom.NODE: lambda x: x,
                    pywrap.NodeFrom.EDGE_SRC: topo.get_src_type,
                    pywrap.NodeFrom.EDGE_DST: topo.get_dst_type}
    self._path = path
    self._type = path_to_type.get(node_from)(path)
    if self._type in self._graph.get_node_decoders().keys():
      self._decoder = self._graph.get_node_decoder(self._type)
    else:
      self._decoder = self._graph.get_edge_decoder(self._type)

  def _add_param(self, name, value):
    self._params[name] = value

  def _add_in_edge(self, edge):
    self._in_edges.append(edge)

  def _add_out_edge(self, edge):
    self._out_edges.append(edge)

  def set_output_field(self, field):
    self._output_field = field

  @property
  def output_field(self):
    return self._output_field

  @property
  def type(self):
    return self._type

  @property
  def shape(self):
    return self._shape

  @property
  def sparse(self):
    return self._sparse

  @property
  def node_def(self):
    return self._node_def

  def get_alias(self):
    return self._alias

  def get_lookup_node(self):
    return self._lookup_node

  def get_degree_nodes(self):
    return self._degree_nodes

  """ GSL APIs """
  def alias(self, alias):
    self._set_alias(alias, temp=False)
    return self

  def batch(self, batch_size):
    assert isinstance(batch_size, int) and batch_size > 0
    self._shape = (batch_size,)
    self._add_param(pywrap.kBatchSize, batch_size)
    self._add_param(pywrap.kEpoch, sys.maxsize >> 32)
    self._add_param(pywrap.kStrategy, "by_order")
    return self

  def shuffle(self, traverse=False):
    strategy_map = {True: "shuffle", False: "random"}
    self._add_param(pywrap.kStrategy, strategy_map.get(traverse, False))
    return self

  def sample(self, count):
    assert isinstance(count, int)
    self._add_param(pywrap.kNeighborCount, count)
    self._add_param(pywrap.kPartitionKey, pywrap.kSrcIds)
    self._shape = (np.prod(self._shape), count)
    return self

  def by(self, strategy):
    self._strategy = strategy
    if self._op_name == "NegativeSampler":
      assert strategy in ["random", "in_degree", "conditional", "node_weight"]
    elif self._op_name == "Sampler":
      assert strategy in \
        ["random", "topk", "in_degree", "edge_weight", "full"]
    else:
      raise ValueError("`by(strategy)` can only be used after`sample(count)`")

    self._sparse = (strategy == "full")
    self._op_name = strategy2op(self._strategy, self._op_name)
    self._add_param(pywrap.kStrategy, self._op_name)
    return self

  def filter(self, target):
    """Filter the samples that are not equal to target ids.

    Args:
      target (string): Alias of upstream TraverseVertexDagNode.

    Raises:
      ValueError: target upstream is not existed.
    """
    if isinstance(target, str):
      target = self._dag.get_node(target)
    if not isinstance(target, TraverseVertexDagNode):
      raise ValueError("filter only accepts upstream Nodes.")
    edge = self._new_edge(src_output=target.output_field,
                          dst_input=pywrap.kFilterIds)
    target._add_out_edge(edge)
    self._add_in_edge(edge)
    self._add_param(pywrap.kFilterType, 1)
    return self

  def where(self, target, condition={}):
    """ Add condition for negative samlpler. Used after `by`.

    Args:
      target (string): Alias of upstream TraverseVertexDagNode which is the
        postive sample that condition should match.
      condition (dict, optional): Keys are as following.
        "batch_share" (bool, optional): Whether sampled negative samples are
          shared by this batch. Defaults to False.
        "unique" (bool, optional): Whether sampled negtive samples are unique.
          Defaults to False.
        "int_cols" (int list, optional): int columns as condition.
          Defaults to [].
        "int_props" (float list, optional) : proportions of int columns.
          Defaults to [].
        "float_cols" (int list, optional): float columns as condition.
          Defaults to [].
        "float_props" (float list, optional): proportions of float columns.
          Defaults to [].
        "str_cols" (int list, optional): string columns as condition.
          Defaults to [].
        "str_props" (float list, optional): proportions of string columns.
          Defaults to [].

    Raises:
      ValueError: target upstream is not existed.
    """
    if isinstance(target, str):
      target = self._dag.get_node(target)
    if not isinstance(target, TraverseVertexDagNode):
      raise ValueError("where only accepts upstream Nodes.")
    edge = self._new_edge(src_output=target.output_field,
                          dst_input=pywrap.kDstIds)
    target._add_out_edge(edge)
    self._add_in_edge(edge)

    default_kvs = {
      "batch_share": (pywrap.kBatchShare, False),
      "unique": (pywrap.kUnique, False),
      "int_cols": (pywrap.kIntCols, None),
      "int_props": (pywrap.kIntProps, None),
      "float_cols": (pywrap.kFloatCols, None),
      "float_props": (pywrap.kFloatProps, None),
      "str_cols": (pywrap.kStrCols, None),
      "str_props": (pywrap.kStrProps, None)
    }

    for k in condition.keys():
      if k not in default_kvs.keys():
        raise ValueError("condition {} is not supported.".format(k))
    for k, v in default_kvs.items():
      param_key, default_value = v
      value = condition.get(k, default_value)
      if value is not None:
        self._add_param(param_key, value)

    self._op_name = "ConditionalNegativeSampler"
    self._add_param(pywrap.kStrategy, self._strategy)
    self._add_param(pywrap.kDstType, target.type)
    return self

  def each(self, func):
    func(self)
    return self

  def values(self, func=lambda x: x):
    self._dag.set_ready(func)
    return self._dag

  """ GSL Apis """

  def _set_alias(self, alias=None, temp=False):
    if self._alias:
      return
    if not alias:
      alias = str(uuid.uuid1())
    self._alias = alias
    self._dag.add_node(alias, self, temp=temp)

    self._lookup_node = self._lookup()
    self._link_to_sink()

  def _get_shape_and_degrees(self, dag_values):
    shape = self._shape
    degrees = None
    if self._sparse:
      assert isinstance(shape, tuple) and len(shape) == 2
      degrees = pywrap.get_dag_value(dag_values, self._nid, pywrap.kDegreeKey)
      shape = (shape[0], shape[1] if shape[1] else max(degrees))
    return shape, degrees

  def feed_values(self, dag_values):
    pass

  def _reverse_edge(self, edge_type, force=True):
    reverse_mask = "_reverse"
    if edge_type.endswith(reverse_mask):
      return edge_type[: -len(reverse_mask)]
    elif force:
      return edge_type + reverse_mask
    return edge_type

  def _new_edge(self, src_output=None, dst_input=None):
    # add an edge for cur node
    eid = get_eid()
    cur_edge = get_dag_edge(eid)
    default_field = "fake"
    cur_edge.src_output = src_output or self._output_field or default_field
    cur_edge.dst_input = dst_input or default_field
    return cur_edge

  def _new_edge_node(self, op_name, edge_type, in_edge):
    assert edge_type is not None and isinstance(edge_type, str)

    self._add_out_edge(in_edge)
    shape = self._shape

    next_node = TraverseEdgeDagNode(self._dag, op_name=op_name,
                                    upstream=self)
    next_node._shape = shape
    next_node._add_param(pywrap.kEdgeType, edge_type)
    next_node._add_in_edge(in_edge)

    next_node.set_path(edge_type, pywrap.NodeFrom.NODE)
    next_node.set_output_field(pywrap.kEdgeIds)
    return next_node

  def _new_vertex_node(self, op_name, edge_type, in_edge,
                       node_from=pywrap.NodeFrom.EDGE_DST):
    assert edge_type is not None and isinstance(edge_type, str)

    self._add_out_edge(in_edge)
    shape = self._shape

    next_node = TraverseVertexDagNode(self._dag, op_name=op_name)
    next_node._shape = shape
    next_node._add_param(pywrap.kEdgeType, edge_type)
    next_node._add_in_edge(in_edge)

    next_node.set_path(edge_type, node_from)
    next_node.set_output_field(pywrap.kNodeIds)
    return next_node

  def set_ready(self, node_id):
    """ Set dag_node ready and format the DagNodeDef proto.
    """
    add_param_map = {
      int: pywrap.add_dag_node_int_params,
      str: pywrap.add_dag_node_string_params,
      unicode: pywrap.add_dag_node_string_params}
    add_vector_param_map = {
      int: pywrap.add_dag_node_int_vector_params,
      float: pywrap.add_dag_node_float_vector_params}

    self._nid = node_id
    node_def = pywrap.new_dag_node()
    pywrap.set_dag_node_id(node_def, node_id)
    pywrap.set_dag_node_op_name(node_def, self._op_name)
    for in_edge in self._in_edges:
      in_edge.dst = self
      pywrap.add_dag_node_in_edge(
        node_def, in_edge.dag_edge_def)
    for out_edge in self._out_edges:
      out_edge.src = self
      pywrap.add_dag_node_out_edge(
        node_def, out_edge.dag_edge_def)

    for k, v in self._params.items():
      if isinstance(v, bool):
        v = int(v)
      if not isinstance(v, list):
        add_param_map[type(v)](node_def, k, v)
        continue
      if len(v) == 0:
        continue
      add_vector_param_map[type(v[0])](node_def, k, v)
    self._node_def = node_def
    return True

  def _lookup(self):
    return None

  def _link_to_sink(self):
    edge = self._new_edge()
    self._add_out_edge(edge)
    self._dag.sink_node._add_in_edge(edge)

class SinkNode(DagNode):
  def __init__(self, dag):
    super(SinkNode, self).__init__(dag, "Sink")
    self._dag.sink_node = self

  def alias(self, alias):
    self._set_alias(alias, temp=True)
    return self

  def _lookup(self):
    # Override
    return None

  def _link_to_sink(self):
    # Override
    pass

class TraverseVertexDagNode(DagNode):
  def __init__(self, dag, op_name="", params={}):
    super(TraverseVertexDagNode, self).__init__(dag, op_name, params)

  def outV(self, edge_type=None):
    self._set_alias()
    in_edge = self._new_edge(dst_input=pywrap.kSrcIds)
    next_node = self._new_vertex_node("Sampler", edge_type, in_edge)
    self._pos_downstreams.append(next_node)
    self._add_degree_node(edge_type, pywrap.NodeFrom.EDGE_SRC)
    return next_node

  def inV(self, edge_type=None):
    self._set_alias()
    in_edge = self._new_edge(dst_input=pywrap.kSrcIds)
    next_node = self._new_vertex_node("Sampler",
                                      self._reverse_edge(edge_type), in_edge)
    self._pos_downstreams.append(next_node)
    self._add_degree_node(self._reverse_edge(edge_type), pywrap.NodeFrom.EDGE_SRC)
    return next_node

  def outE(self, edge_type):
    self._set_alias()
    in_edge = self._new_edge(dst_input=pywrap.kSrcIds)
    next_node = self._new_edge_node("Sampler", edge_type, in_edge)
    self._pos_downstreams.append(next_node)
    return next_node

  def inE(self, edge_type):
    self._set_alias()
    in_edge = self._new_edge(dst_input=pywrap.kSrcIds)
    next_node = self._new_edge_node("Sampler",
                                    self._reverse_edge(edge_type), in_edge)
    self._pos_downstreams.append(next_node)
    return next_node

  def outNeg(self, edge_type):
    self._set_alias()
    in_edge = self._new_edge(dst_input=pywrap.kSrcIds)
    next_node = self._new_vertex_node("NegativeSampler", edge_type, in_edge)
    self._neg_downstreams.append(next_node)
    return next_node

  def inNeg(self, edge_type):
    self._set_alias()
    in_edge = self._new_edge(dst_input=pywrap.kSrcIds)
    next_node = self._new_vertex_node("NegativeSampler",
                                      self._reverse_edge(edge_type), in_edge)
    self._neg_downstreams.append(next_node)
    return next_node

  def Neg(self, node_type):
    self._set_alias()
    in_edge = self._new_edge(dst_input=pywrap.kSrcIds)
    next_node = self._new_vertex_node("NegativeSampler", node_type, in_edge,
                                      node_from=pywrap.NodeFrom.NODE)
    self._neg_downstreams.append(next_node)
    return next_node

  def _lookup(self):
    # Override
    # Generate an edge from traverse node to it's lookup node.
    edge = self._new_edge(dst_input=pywrap.kNodeIds)
    self._add_out_edge(edge)
    return LookupDagNode("LookupNodes", self, [edge],
                         {pywrap.kNodeType: self._type,
                          pywrap.kPartitionKey: pywrap.kNodeIds})

  def feed_values(self, dag_values):
    shape, degrees = self._get_shape_and_degrees(dag_values)
    return self._graph.get_nodes(
      self._type, pywrap.get_dag_value(
        dag_values, self._nid, pywrap.kNodeIds), degrees, shape)

class TraverseNegVertexDagNode(TraverseVertexDagNode):
  def __init__(self, dag, op_name="", params={}):
    super(TraverseNegVertexDagNode, self).__init__(dag, op_name, params)

class TraverseEdgeDagNode(DagNode):
  def __init__(self, dag, op_name="", params={}, upstream=None):
    super(TraverseEdgeDagNode, self).__init__(dag, op_name, params)
    self._upstream = upstream

  def inV(self):
    self._set_alias()
    next_node = FakeNode(self)
    next_node.set_path(self._path, pywrap.NodeFrom.EDGE_DST)
    next_node.set_output_field(pywrap.kNodeIds)
    return next_node

  def outV(self):
    raise warnings.warn("outV is just the upstream Nodes.")

  def _lookup(self):
    # Override
    # Generate an edge from current node to it's lookup node.
    edge = self._new_edge(dst_input=pywrap.kEdgeIds)
    # Add an edge from upstream node to lookup node.
    extra_edge = self._new_edge(
      src_output=self._upstream.output_field, dst_input=pywrap.kSrcIds)
    self._add_out_edge(edge)
    self._upstream._add_out_edge(extra_edge)

    return LookupDagNode(
        "LookupEdges", self, [edge, extra_edge],
        {pywrap.kEdgeType: self._reverse_edge(self._type, False),
         pywrap.kPartitionKey: pywrap.kSrcIds})

  def feed_values(self, dag_values):
    shape, degrees = self._get_shape_and_degrees(dag_values)
    edge_ids = pywrap.get_dag_value(dag_values, self._nid, pywrap.kEdgeIds)
    assert isinstance(shape, tuple) and len(shape) == 2
    nbr_counts = degrees if self._sparse else [shape[1]] * shape[0]
    src_ids = pywrap.get_dag_value(dag_values,
                                   self._upstream.nid,
                                   pywrap.kNodeIds)
    src_ids = src_ids.reshape(self._upstream.shape)
    src_ids = np.concatenate(
      [src_ids[idx].repeat(d) for idx, d in enumerate(nbr_counts)])
    dst_ids = pywrap.get_dag_value(dag_values, self._nid, pywrap.kNodeIds)
    return self._graph.get_edges(
        self._type, src_ids, dst_ids, edge_ids, degrees, shape)

class TraverseSourceEdgeDagNode(TraverseEdgeDagNode):
  def __init__(self, dag, op_name="", params={}):
    super(TraverseSourceEdgeDagNode, self).__init__(
      dag, op_name=op_name, params=params)

  def outV(self):
    self._set_alias()
    next_node = FakeNode(self)
    next_node.set_path(self._path, pywrap.NodeFrom.EDGE_SRC)
    next_node.set_output_field(pywrap.kSrcIds)
    return next_node

  def inV(self, edge_type=None):
    self._set_alias()
    next_node = FakeNode(self)
    next_node.set_path(self._path, pywrap.NodeFrom.EDGE_DST)
    next_node.set_output_field(pywrap.kDstIds)
    return next_node

  def _lookup(self):
    # Override
    edge = self._new_edge(dst_input=pywrap.kEdgeIds)
    extra_edge = self._new_edge(src_output=pywrap.kSrcIds,
                                dst_input=pywrap.kSrcIds)
    self._add_out_edge(edge)
    self._add_out_edge(extra_edge)

    return LookupDagNode(
        "LookupEdges", self, [edge, extra_edge],
        {pywrap.kEdgeType: self._reverse_edge(self._type, False),
         pywrap.kPartitionKey: pywrap.kSrcIds})

  def feed_values(self, dag_values):
    shape, degrees = self._get_shape_and_degrees(dag_values)
    edge_ids = pywrap.get_dag_value(dag_values, self._nid, pywrap.kEdgeIds)
    src_ids = pywrap.get_dag_value(dag_values, self._nid, pywrap.kSrcIds)
    dst_ids = pywrap.get_dag_value(dag_values, self._nid, pywrap.kDstIds)
    return self._graph.get_edges(
        self._type, src_ids, dst_ids, edge_ids, degrees, shape)

class LookupDagNode(DagNode):
  def __init__(self, op_name="", upstream=None, in_edges=[], params={}):
    super(LookupDagNode, self).__init__(upstream._dag, op_name, params)
    self._upstream = upstream
    self._shape = upstream._shape
    for edge in in_edges:
      self._add_in_edge(edge)
    self._set_alias(temp=True)
    self.set_output_field("properties")

class DegreeDagNode(DagNode):
  def __init__(self, op_name="", upstream=None, in_edges=[], params={}):
    super(DegreeDagNode, self).__init__(upstream._dag, op_name, params)
    self._upstream = upstream
    self._shape = upstream._shape
    for edge in in_edges:
      self._add_in_edge(edge)
    self._set_alias(temp=True)
    self.set_output_field(pywrap.kDegrees)

  @property
  def edge_type(self):
    return self._params[pywrap.kEdgeType]

  @property
  def node_from(self):
    return self._params[pywrap.kNodeFrom]

class FakeNode(TraverseVertexDagNode):
  """ FakeNode is used for adding corresponding DagNode of E.outV()/inV()
  to Dag. E.outV()/inV() doesn't raise any operator, but only changes the
  field which downstream absorbs.
  """
  def __init__(self, dag_node):
    super(FakeNode, self).__init__(dag_node._dag)
    if not isinstance(dag_node, TraverseEdgeDagNode):
      raise ValueError("FakeNode is a fake for TraverseEdgeDagNode, not {}"
                       .format(type(dag_node)))
    self._upstream = dag_node
    self._shape = dag_node._shape

  @property
  def nid(self):
    return self._upstream.nid

  def set_ready(self, node_id):
    # Override
    return False

  def _add_out_edge(self, edge):
    # Override
    self._upstream._add_out_edge(edge)

  def _link_to_sink(self):
    # Override
    pass

  def feed_values(self, dag_values):
    edges = self._upstream.feed_values(dag_values)
    return (edges.dst_nodes, edges.src_nodes)[int(self._output_field == pywrap.kSrcIds)]

max_stars_repo_name: gasdaf/graph-learn

max_issues_repo_name: gasdaf/graph-learn

max_forks_repo_name: gasdaf/graph-learn


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri May  1 17:48:08 2020

@author: hovhannesstepanyan
"""

from Fuzzy import FS_plot_cat, FS_intersect
import numpy as np

class InteligentSystems:

    def FS_func(FSx, Uy, Func = lambda x, Ux, Uy: np.array([1]*len(Uy)), method='minimax'):
        masx = FSx.items()
        Ux = np.array(list(masx)).T[0]
        masmu = np.array([0] * len(Uy))
        if method == 'minimax':
            for el in Ux:
                masmu=np.max([masmu, np.min([Func(el, Ux, Uy), [FSx[el]] * len(Uy)], axis=0)], axis=0)
        elif method == 'probability':
            for el in Ux:
                mn = Func(el, masx, Uy) * FSx[el]
                masmu = masmu + mn - masmu * mn
        else:
            print("Unknown method", method)
        return dict(np.array([Uy, masmu]).T)

    def FS_f(x, Ux, Uy):
        minx = np.min(Ux)
        miny = np.min(Uy)
        maxy = np.max(Uy)
        k = float(x - minx) / float(Ux.max() - minx)
        gran = k * (maxy - miny) + miny
        return np.array([1. if y == gran else
                         float(y - miny) / (gran - miny) if y < gran else
                         float(maxy - y) / (maxy - gran) for y in Uy])

    def FClass(Ux, Func = lambda x, Ux: np.array([1]*len(Ux)), method='base'):
        massiv = np.array([Func(el, Ux) for el in Ux])
        if method == 'base':
            masmu = np.min(massiv / np.max([list(zip(el, massiv.T[nom]))
                                           for nom, el in enumerate(massiv)],
                                          axis = 2), axis = 1)
        elif method == 'alternative':
            masmu = np.min(massiv / massiv.T, axis = 1)
        else:
            masmu = np.array([0] * len(Ux))
            print('Unknown method', method)
        return dict(np.array([Ux, masmu]).T)

    def FSU(x, Fg):
        Fgt = Fg.T
        ress = 0.
        if x <= np.min(Fgt[0]):
            ress = Fg[np.argmin(Fgt[0])][2]
        elif x >= np.max(Fgt[1]):
            ress = Fg[np.argmax(Fgt[1])][3]
        else:
            kt = Fg[(Fgt[0] <= x) & (Fgt[1] >= x)][0]
            ress = (x - kt[0]) / (kt[1] - kt[0]) * (kt[3] - kt[2]) + kt[2]
        return ress

    def free_FS(U, Fg, Func=FSU):
        return dict([[el, Func(el, Fg)] for el in U])

    def free_FS_cat(X, Fg, Func=FSU):
        return dict([[nom, Func(el, Fg)] for nom, el in enumerate(X)])

    def dis_help_sys(datas):
        names = datas[0]
        names.pop()
        names.pop(0)
        Ca_cat = names
        arrayToPlot = []
        criterias = []
        for i in range(1, len(datas)):
            row = datas[i]
            T = row.pop()
            criteria = row.pop(0)
            FCS = InteligentSystems.make_FG(min(row), max(row), T)
            criterias.append(criteria)
            FS_cat = InteligentSystems.free_FS_cat(list(map(float, row)), FCS)
            arrayToPlot.append(FS_cat)

        FR1c = FS_intersect(arrayToPlot)
        FR2c = FS_intersect(arrayToPlot, method='probability')
        FS_plot_cat([FR1c, FR2c], U=Ca_cat,
                    labels=['minimax', 'probability'], title='Արդյունք',
                   colors=['red', 'green', 'blue', 'yellow', 'green', 'purple'],
                   name='result')

    def make_FG(x1, x2, T):
        if T == 'U':
            return np.array([[float(x1) - 10**len(x1), float(x2) + 10**len(x2), 0., 1.]])
        elif T == 'N':
            return np.array([[float(x1) - 10**len(x1), float(x2) + 10**len(x2), 1., 0.]])
        else:
            return np.array([[0., 1., 0., 1.]])

max_stars_repo_name: Hovo-Infinity/Expert-System

max_issues_repo_name: Hovo-Infinity/Expert-System

max_forks_repo_name: Hovo-Infinity/Expert-System


text:
#!/usr/bin/env python2
# -*- coding: utf-8 -*-

"""
@author: misiak

"""
import uproot
import numpy as np
import re
import os
import abc

class Artifact(object):
    """ Named object. Used to define namespace and later defining attributes.
    """
    def __init__(self, name):
        self.name = name

class Copse(object):
    """ Defining the cuts in a kinda clean/handy way.
    """

    def new_cut(self, name, thruth_array):
        setattr(self, name, thruth_array)

class Tree(object):
    """ Tree class, mimicking ROOT's Tree.
    """

    def __init__(self, root, tree_name):

        self._tree = root[tree_name]

        self._keys = list()
        for k in self._tree.keys():
            # sanitize keys and values
            key = k.decode('utf-8')
            (self._keys).append(key)

            value = np.array(self._tree[k].array())
            if value.ndim > 1 and value.shape[0]==1:
                value = value[0]
            setattr(self, key, value)

class Root_proto(object):
    """ Simpler root class for trigger and noise data.
    """
    def __init__(self, file_path):

        self._root = uproot.open(file_path)

        key_list = []
        for k in self._root.keys():
            key = k.decode('utf-8').split(';')[0]
            key_list.append(key)
        key_list = list(set(key_list))

        self._keys = key_list

        for key in self._keys:
            tree = Tree(
                    self._root,
                    key
            )
            setattr(self, key, tree)

class Guardian_proto(object):
    """ Concanetation of Root partitions.
    """
    def __init__(self, root_list):

        self.roots = root_list

        # self.run_tree for
        # stacking all the run_tree attributes
        # by listing then numpy.stack-ing
        # and
        # self.data_types.tree_labels for
        # concatenate all the data_types.tree_labels attributes
        # by listing then numpy.concatenate-ing

        # initializing these attributes with empty lists
        root0 = self.roots[0]

        self._keys = root0._keys

        for key in self._keys:

            tree = getattr(root0, key)
            setattr(self, key, Artifact(key))
            guard_tree = getattr(self, key)
            setattr(guard_tree, '_keys', tree._keys)

            for subkey in guard_tree._keys:
                setattr(guard_tree, subkey, list())
                guardian_arti = getattr(guard_tree, subkey)

                for root in self.roots:
                    root_tree = getattr(root, key)
                    root_arti = getattr(root_tree, subkey)

                    guardian_arti.append(root_arti)

                arti_array = np.concatenate(guardian_arti, axis=0)
                setattr(guard_tree, subkey, arti_array)

class Root_reader_proto(object):
    """ Analysis class.
    """

    def __init__(self, run, detector='RED71',
                 run_dir='/home/misiak/Data/data_run59',
                 data_type='trigger'):

        self.run = run
        self.detector = detector

        self.dir = '/'.join((run_dir, self.run, self.detector, 'RootFiles'))

        self.files = []
        for fname in os.listdir(self.dir):
            if data_type.lower() == 'trigger':
                if 'TriggerData' in fname:
                    self.files.append(fname)
            if data_type.lower() == 'noise':
                if 'NoiseData' in fname:
                    self.files.append(fname)

        (self.files).sort()

        self._part_num = [re.findall('S([0-9]+)_', f)[0] for f in self.files]

        self._part_label = ['part_'+ num for num in self._part_num]

        self.roots = list()

        for fname, label in zip(self.files, self._part_label):
            fpath = '/'.join((self.dir, fname))
            part_root = Root_proto(fpath)
            (self.roots).append(part_root)

        self.all = Guardian_proto(self.roots)

class Wisp(object):
    """ Preparing inheritance for Root and Thicket classes.
    """
    __metaclass__ = abc.ABCMeta

    def __init__(self):

        self.run_tree = Artifact('run_tree')

        self._event_types = ('trig', 'noise')
        self._processing_types = ('raw', 'filt', 'filt_decor')

        for etype in self._event_types:
            setattr(self, etype, Artifact(etype))

            for ptype in self._processing_types:
                arti = getattr(self, etype)
                setattr(arti, ptype, Artifact(ptype))

class Root(Wisp):
    """ Root class, mimicking ROOT's Root.
    """
    def __init__(self, file_path):

        Wisp.__init__(self)

        self._root = uproot.open(file_path)

        # general info about the run
        self.run_tree = Tree(self._root, 'RunTree_Normal')

        for etype in self._event_types:
            arti = getattr(self, etype)

            for ptype in self._processing_types:
                tree = Tree(self._root,
                            'EventTree_{}_Normal_{}'.format(etype, ptype))
                setattr(arti, ptype, tree)

class Guardian(Wisp):
    """ Concanetation of Root partitions.
    """
    def __init__(self, root_list):

        Wisp.__init__(self)

        self.roots = root_list

        # self.run_tree for
        # stacking all the run_tree attributes
        # by listing then numpy.stack-ing
        # and
        # self.data_types.tree_labels for
        # concatenate all the data_types.tree_labels attributes
        # by listing then numpy.concatenate-ing

        # initializing these attributes with empty lists
        root0 = self.roots[0]
        self.run_tree._keys = root0.run_tree._keys
        for key in self.run_tree._keys:
            setattr(self.run_tree, key, list())

        for etype in self._event_types:
            root_arti = getattr(root0, etype)
            guard_arti = getattr(self, etype)

            for ptype in self._processing_types:
                root_tree = getattr(root_arti, ptype)
                guard_tree = getattr(guard_arti, ptype)
                setattr(guard_tree, '_keys', root_tree._keys)

                for key in guard_tree._keys:
                    setattr(guard_tree, key, list())

        # appending the attributes of all root files to the empty list
        for root in self.roots:
            for key in self.run_tree._keys:
                root_attr = getattr(root.run_tree, key)
                attr_list = getattr(self.run_tree, key)
                attr_list.append(root_attr)

            for etype in self._event_types:
                root_arti = getattr(root, etype)
                guard_arti = getattr(self, etype)

                for ptype in self._processing_types:
                    root_tree = getattr(root_arti, ptype)
                    guard_tree = getattr(guard_arti, ptype)

                    for key in guard_tree._keys:
                        root_attr = getattr(root_tree, key)
                        guard_attr = getattr(guard_tree, key)
                        guard_attr.append(root_attr)

        # stacking for self.run_tree
        for key in self.run_tree._keys:
            attr_list = getattr(self.run_tree, key)
            attr_array = np.stack(attr_list, axis=0)
            setattr(self.run_tree, key, attr_array)

        # concatenating for self.data_types.tree_labels
        for etype in self._event_types:
            arti = getattr(self, etype)

            for ptype in self._processing_types:
                tree = getattr(arti, ptype)

                for key in tree._keys:
                    attr_list = getattr(tree, key)
                    attr_array = np.concatenate(attr_list, axis=0)
                    setattr(tree, key, attr_array)

            arti.nsamples = attr_array.shape[0]

class Root_reader(object):
    """ Analysis class.
    """

    def _temporal_data_extraction(self):
        """
        Extract data concerning time, freq, maintenance etc..
        Should be part of the original root. No arbitrary decision made here.
        """
        run_tree = self.all.run_tree

        # acquisition frequency
        freq_root = np.ravel((run_tree.f_max_heat, run_tree.f_max_ion))
        assert np.all(freq_root == freq_root[0])
        run_tree.freq = freq_root[0]

        # maintenance info
        maint_cycle = run_tree.MaintenanceCycle
        assert np.all(maint_cycle == maint_cycle[0])
        run_tree.maint_cycle = maint_cycle[0] / 3600 # in hours

        maint_duration = run_tree.MaintenanceDuration
        assert np.all(maint_duration == maint_duration[0])
        run_tree.maint_duration = maint_duration[0] / 3600 # in hours

    def _time_stamp_extraction(self):
        """ time_array in hours """

        for etype in self.all._event_types:
            etype_tree = getattr(self.all, etype)

            for ptype in self.all._processing_types:

                etype_attr = getattr(etype_tree, ptype)

                run_tree = self.all.run_tree
                secondes_array = etype_attr.MicroStp / run_tree.freq
                hours_array =  secondes_array / 3600 + etype_attr.NumPart

                setattr(etype_attr, 'time_stream', hours_array)

    def __init__(self, run, detector='RED80',
                 run_dir='/home/misiak/Data/data_run57'):

        self.run = run
        self.detector = detector

        self.dir = '/'.join((run_dir, self.run, self.detector))

        self.files = []
        for fname in os.listdir(self.dir):
            if 'ProcessedData' in fname:
                self.files.append(fname)

        (self.files).sort()

        self._part_num = [re.findall('S([0-9]+)_', f)[0] for f in self.files]

        self._part_label = ['part_'+ num for num in self._part_num]

        self.roots = list()

        for fname, label in zip(self.files, self._part_label):
            fpath = '/'.join((self.dir, fname))
            part_root = Root(fpath)
            (self.roots).append(part_root)

        self.all = Guardian(self.roots)

        self.all.trig.ref = self.all.trig.filt_decor

        self._temporal_data_extraction()
        self._time_stamp_extraction()

        # definng quality cuts
        self.all.trig.cut = Copse()
        self.all.noise.cut = Copse()

    def maintenance_cut(self, time_array):
        """ checking if the considered time is within a maintenance """
        run_tree = self.all.run_tree
        full_cycle = run_tree.maint_duration + run_tree.maint_cycle
        remainder_array = time_array % full_cycle
        # thruth array, True if time not in a maintenance
        return remainder_array > run_tree.maint_duration