text:
import numpy as np
from collections import deque

class Time_calc:
    def __init__(self, parameters, episode):
        self.parmeters = parameters

        """Random Seed for random numbers"""
        np.random.seed(parameters['SEED'] + episode)
        self.randomStreams = {}
        self.randomStreams["process_time"] = [np.random.RandomState(np.random.randint(100)) for i in range(parameters['NUM_MACHINES'])]
        self.randomStreams["machine_failure"] = [np.random.RandomState(np.random.randint(100)) for i in range(parameters['NUM_MACHINES'])]
        self.randomStreams["repair_time"] = [np.random.RandomState(np.random.randint(100)) for i in range(parameters['NUM_MACHINES'])]
        self.randomStreams["order_generation"] = [np.random.RandomState(np.random.randint(100)) for i in range(parameters['NUM_SOURCES'])]
        self.randomStreams["order_sequence"] = np.random.RandomState(np.random.randint(100))
        self.randomStreams["transp_agent"] = [np.random.RandomState(np.random.randint(100)) for i in range(parameters['NUM_TRANSP_AGENTS'])]
        self.randomStreams["filled_initial_system"] = np.random.RandomState(np.random.randint(100))

    """
    Utility procedures
    """
    def get_inventory_level(self, statistics):
        return statistics['stat_inv_episode'][-1][1]

    def processing_time(self, machine, statistics, parameters, order):
        """Return actual processing time for a concrete part."""
        result_time = min(parameters['MAX_PROCESS_TIME'][machine.id], max(parameters['MIN_PROCESS_TIME'][machine.id], self.randomStreams["process_time"][machine.id].exponential(scale=parameters['AVERAGE_PROCESS_TIME'][machine.id])))
        statistics['stat_machines_working'][machine.id] += result_time
        return result_time

    def transp_time(self, start, end, transp, statistics, parameters):
        """Return actual processing time for a concrete part."""
        result_time = parameters['TRANSP_TIME'][start.id][end.id]
        statistics['stat_transp_walking'][transp.id] += result_time
        statistics['stat_transp_working'][transp.id] += result_time
        return result_time

    def handling_time(self, MachineOrSource, LoadOrUnload, transp, statistics, parameters):
        """Return actual handling time for a concrete part."""
        result_time = 0
        if MachineOrSource == "machine":
            if LoadOrUnload == "load":
                result_time += parameters['TIME_TO_LOAD_MACHINE']
            elif LoadOrUnload == "unload":
                result_time += parameters['TIME_TO_UNLOAD_MACHINE']
        elif MachineOrSource == "source":
            if LoadOrUnload == "load":
                result_time += parameters['TIME_TO_LOAD_SOURCE']
            elif LoadOrUnload == "unload":
                result_time += parameters['TIME_TO_UNLOAD_SOURCE']
        statistics['stat_transp_handling'][transp.id] += result_time
        statistics['stat_transp_working'][transp.id] += result_time
        return result_time

    def changeover_time(self, machine, current_variant, next_variant, statistics, parameters):
        """Return actual changeover time for a machine."""
        result_time = parameters['CHANGEOVER_TIME']
        statistics['stat_machines_changeover'][machine.id] += result_time
        return result_time

    def time_to_failure(self, machine, statistics, parameters):
        """Return time until next failure for a machine."""
        result_time = self.randomStreams["machine_failure"][machine.id].exponential(scale=parameters['MTBF'][machine.id])
        return result_time

    def repair_time(self, machine, statistics, parameters):
        """Return time until next failure for a machine."""
        result_time = self.randomStreams["repair_time"][machine.id].exponential(scale=parameters['MTOL'][machine.id])
        statistics['stat_machines_broken'][machine.id] += result_time
        return result_time

    def time_to_order_generation(self, source, statistics, parameters):
        """Return time until next failure for a machine."""
        result_time = self.randomStreams["order_generation"][source.id - self.parmeters['NUM_MACHINES']].exponential(scale=parameters['MTOG'][source.id])
        return result_time

    def create_intermediate_production_steps_and_variant(self, statistics, parameters, resources, at_resource):
        """Return random sequence of production steps based on the distribution that is defined over all machines."""
        if at_resource.type == "source":
            result_prod_steps = [at_resource]
            while True:
                first_machine = self.randomStreams["order_sequence"].choice(resources['machines'], 1, p=parameters['ORDER_DISTRIBUTION'])[0]
                if first_machine.id in at_resource.resp_area:
                    break
            result_prod_steps.append(first_machine)
            result_prod_steps.extend(self.randomStreams["order_sequence"].choice(resources['machines'], parameters['NUM_PROD_STEPS'] - 1, p=parameters['ORDER_DISTRIBUTION']))
        else:
            first_machine = at_resource
            result_prod_steps = [x for x in resources['sources'] if first_machine.id in x.resp_area]
            result_prod_steps.append(first_machine)
            actual_step = self.randomStreams["filled_initial_system"].randint(1, high=parameters['NUM_PROD_STEPS'])
            result_prod_steps.extend(self.randomStreams["order_sequence"].choice(resources['machines'], parameters['NUM_PROD_STEPS'] - 1 - actual_step, p=parameters['ORDER_DISTRIBUTION']))
        if result_prod_steps[-1].id < 2:
            result_prod_steps.append(resources['sinks'][0])
        elif result_prod_steps[-1].id < 5:
            result_prod_steps.append(resources['sinks'][1])
        elif result_prod_steps[-1].id < 8:
            result_prod_steps.append(resources['sinks'][2])
        result_variant = self.randomStreams["order_sequence"].choice(parameters['NUM_PROD_VARIANTS'], 1, p=parameters['VARIANT_DISTRIBUTION'])

        return result_prod_steps, result_variant

def update_mov_avg(**kwargs):
    """ Function to iteratively calculate moving average with a given window"""
    kwargs['cont'].appendleft(kwargs['value'])
    new_mean = sum(kwargs['cont']) / len(kwargs['cont'])
    return new_mean

def update_mov_std(**kwargs):
    """ Function to iteratively calculate moving std with a given window"""
    kwargs['cont_sq'].appendleft(kwargs['value'] ** 2)
    len_wd = len(kwargs['cont_sq'])
    new_var = ((1 / len_wd) * sum(kwargs['cont_sq'])) - ((1 / len_wd) * sum(kwargs['cont'])) ** 2

    if new_var < 0:
        new_var = 0
    return np.sqrt(new_var)

def update_exp_weighted_mean(**kwargs):
    """ Function to iteratively calculate exponentially weighted mean"""
    new_mean = (1 - kwargs['alpha']) * kwargs['oldMean'] + kwargs['alpha'] * kwargs['value']
    return new_mean

def update_exp_weightes_std(**kwargs):
    """ Function to iteratively calculate exponentially weighted std"""
    diff = kwargs['value'] - kwargs['oldMean']
    incr = kwargs['alpha'] * diff
    new_var = (1 - kwargs['alpha']) * (kwargs['oldStd'] ** 2 + diff * incr)
    return np.sqrt(new_var)

class ZScoreNormalization(object):
    running_mean = dict(
        exp=update_exp_weighted_mean,
        mov=update_mov_avg
    )
    running_std = dict(
        exp=update_exp_weightes_std,
        mov=update_mov_std
    )

    def __init__(self, type, **config):
        """ Calculates the z score normalization for a given method of calculating the mean and std """
        self.update_mean = ZScoreNormalization.running_mean[type]
        self.update_std = ZScoreNormalization.running_std[type]
        self.config = dict(**config)
        self.type = type
        self.counter = 0
        self.mean = 0
        self.std = 0
        self.setup()

    def __call__(self, value):
        """ By calling this function calculates the mean, std and z score normalization of the next iteration with given value """
        self.counter += 1
        update = dict(
            oldMean=self.mean,
            value=value,
            oldStd=self.std,
            counter=self.counter,
            **self.attr_alg
        )
        self.mean = self.update_mean(**update)
        self.std = self.update_std(**update)

    def setup(self):
        """ Setup function to set all parameters for calculations"""
        if self.type is 'mov':
            size_wd = self.config['window']
            self.attr_alg = dict(cont_sq=deque([], maxlen=size_wd), cont=deque([], maxlen=size_wd), **self.config)
        elif self.type is 'exp':
            self.attr_alg = dict(**self.config)

    def get_z_score_normalization(self, value):
        if self.std != 0 and value != None:
            normalized = (value - self.mean) / self.std
            return normalized
        elif value == None:
            return -1
        else:
            return 0

    def reset(self):
        self.counter = 0
        self.mean = 0
        self.std = 0
        self.setup()

max_issues_repo_name: AndreasKuhnle/SimPyRLFab

max_forks_repo_name: AndreasKuhnle/SimPyRLFab


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Oct 22 11:29:57 2018

@author: virati
File for patient specific readout
"""

#DBSpace libraries and sublibraries here
from DBSpace.readout import BR_DataFrame as BRDF
from DBSpace.readout import ClinVect, DSV
from DBSpace.readout.BR_DataFrame import BR_Data_Tree

# General python libraries
import scipy.signal as sig
import numpy as np
import pickle

# Plotting Libraries
import matplotlib.pyplot as plt
import seaborn as sns
#Do some cleanup of the plotting space
plt.close('all')
sns.set_context('paper')
sns.set(font_scale=5)

#%%

test_scale = 'HDRS17'
do_detrend = 'None'

#%%

ClinFrame = ClinVect.CFrame(norm_scales=True)
#BRFrame = BRDF.BR_Data_Tree(preFrame='Chronic_Frame.pickle')
BRFrame = pickle.load(open('/home/virati/Chronic_Frame.pickle',"rb"))

# The below has been folded into the class
# Run our main sequence for populating the BRFrame
#BRFrame.full_sequence(data_path='/home/virati/Chronic_Frame_july.npy')
#BRFrame.check_empty_phases() # Check to see if there are any empty phases. This should be folded into full_sequence soon TODO

#%%
# This sets up the regression
analysis = DSV.ORegress(BRFrame,ClinFrame)

#%%

model = 'ENR_Osc'
for pt in ['901','903','905','906','907','908']:
    print('Patient ' + pt + ' specific model')

    analysis.O_regress(method=model,doplot=False,avgweeks=True,ignore_flags=False,circ='day',scale=test_scale,lindetrend=do_detrend,train_pts=pt,pt_specific=True)
    coeff_runs = analysis.O_models(plot=False,models=[model])
    summ_stats_runs = analysis.Clinical_Summary(model,plot_indiv=False,ranson=True,doplot=False)

    analysis.plot_model_coeffs(model=model,pt=pt)
    # need to now work with the output model but in a reasonable way
    #want to plot the coefficients learned

ext: py

reason: import numpy,import scipy


text:
# https://deeplearningcourses.com/c/unsupervised-machine-learning-hidden-markov-models-in-python
# https://udemy.com/unsupervised-machine-learning-hidden-markov-models-in-python
# http://lazyprogrammer.me
# Continuous-observation HMM in Theano using gradient descent.
from __future__ import print_function, division
from builtins import range
# Note: you may need to update your version of future
# sudo pip install -U future

import wave
import theano
import theano.tensor as T
import numpy as np
import matplotlib.pyplot as plt

# from theano.sandbox import solve # does not have gradient functionality
from generate_c import get_signals, big_init

def random_normalized(d1, d2):
    x = np.random.random((d1, d2))
    return x / x.sum(axis=1, keepdims=True)

class HMM:
    def __init__(self, M, K):
        self.M = M # number of hidden states
        self.K = K # number of Gaussians

    def fit(self, X, learning_rate=1e-2, max_iter=10):
        # train the HMM model using gradient descent

        N = len(X)
        D = X[0].shape[1] # assume each x is organized (T, D)

        pi0 = np.ones(self.M) / self.M # initial state distribution
        A0 = random_normalized(self.M, self.M) # state transition matrix
        R0 = np.ones((self.M, self.K)) / self.K # mixture proportions
        mu0 = np.zeros((self.M, self.K, D))
        for i in range(self.M):
            for k in range(self.K):
                random_idx = np.random.choice(N)
                x = X[random_idx]
                random_time_idx = np.random.choice(len(x))
                mu0[i,k] = x[random_time_idx]
        sigma0 = np.zeros((self.M, self.K, D, D))
        for j in range(self.M):
            for k in range(self.K):
                sigma0[j,k] = np.eye(D)

        thx, cost = self.set(pi0, A0, R0, mu0, sigma0)

        pi_update = self.pi - learning_rate*T.grad(cost, self.pi)
        pi_update = pi_update / pi_update.sum()

        A_update = self.A - learning_rate*T.grad(cost, self.A)
        A_update = A_update / A_update.sum(axis=1).dimshuffle(0, 'x')

        R_update = self.R - learning_rate*T.grad(cost, self.R)
        R_update = R_update / R_update.sum(axis=1).dimshuffle(0, 'x')

        updates = [
            (self.pi, pi_update),
            (self.A, A_update),
            (self.R, R_update),
            (self.mu, self.mu - learning_rate*T.grad(cost, self.mu)),
            (self.sigma, self.sigma - learning_rate*T.grad(cost, self.sigma)),
        ]

        train_op = theano.function(
            inputs=[thx],
            updates=updates,
        )

        costs = []
        for it in range(max_iter):
            print("it:", it)

            for n in range(N):
                c = self.log_likelihood_multi(X).sum()
                print("c:", c)
                costs.append(c)
                train_op(X[n])

        print("A:", self.A.get_value())
        print("mu:", self.mu.get_value())
        print("sigma:", self.sigma.get_value())
        print("R:", self.R.get_value())
        print("pi:", self.pi.get_value())

        plt.plot(costs)
        plt.show()

    def set(self, pi, A, R, mu, sigma):
        self.pi = theano.shared(pi)
        self.A = theano.shared(A)
        self.R = theano.shared(R)
        self.mu = theano.shared(mu)
        self.sigma = theano.shared(sigma)
        M, K = R.shape
        self.M = M
        self.K = K

        D = self.mu.shape[2]
        twopiD = (2*np.pi)**D

        # set up theano variables and functions
        thx = T.matrix('X') # represents a TxD matrix of sequential observations
        def mvn_pdf(x, mu, sigma):
            k = 1 / T.sqrt(twopiD * T.nlinalg.det(sigma))
            e = T.exp(-0.5*(x - mu).T.dot(T.nlinalg.matrix_inverse(sigma).dot(x - mu)))
            return k*e

        def gmm_pdf(x):
            def state_pdfs(xt):
                def component_pdf(j, xt):
                    Bj_t = 0
                    for k in range(self.K):
                        Bj_t += self.R[j,k] * mvn_pdf(xt, self.mu[j,k], self.sigma[j,k])
                    return Bj_t

                Bt, _ = theano.scan(
                    fn=component_pdf,
                    sequences=T.arange(self.M),
                    n_steps=self.M,
                    outputs_info=None,
                    non_sequences=[xt],
                )
                return Bt

            B, _ = theano.scan(
                fn=state_pdfs,
                sequences=x,
                n_steps=x.shape[0],
                outputs_info=None,
            )
            return B.T

        B = gmm_pdf(thx)
        # scale = T.zeros((thx.shape[0], 1), dtype=theano.config.floatX)
        # scale[0] = (self.pi*B[:,0]).sum()

        def recurrence(t, old_a, B):
            a = old_a.dot(self.A) * B[:, t]
            s = a.sum()
            return (a / s), s

        [alpha, scale], _ = theano.scan(
            fn=recurrence,
            sequences=T.arange(1, thx.shape[0]),
            outputs_info=[self.pi*B[:,0], None],
            n_steps=thx.shape[0]-1,
            non_sequences=[B],
        )

        cost = -T.log(scale).sum()
        self.cost_op = theano.function(
            inputs=[thx],
            outputs=cost,
        )
        return thx, cost

    def log_likelihood_multi(self, X):
        return np.array([self.cost_op(x) for x in X])

def real_signal():
    spf = wave.open('helloworld.wav', 'r')

    #Extract Raw Audio from Wav File
    # If you right-click on the file and go to "Get Info", you can see:
    # sampling rate = 16000 Hz
    # bits per sample = 16
    # The first is quantization in time
    # The second is quantization in amplitude
    # We also do this for images!
    # 2^16 = 65536 is how many different sound levels we have
    signal = spf.readframes(-1)
    signal = np.fromstring(signal, 'Int16')
    T = len(signal)
    signal = (signal - signal.mean()) / signal.std()

    hmm = HMM(5, 3)
    # signal needs to be of shape N x T(n) x D
    hmm.fit(signal.reshape(1, T, 1), learning_rate=1e-5, max_iter=20)

def fake_signal():
    signals = get_signals()
    hmm = HMM(5, 3)
    hmm.fit(signals)
    L = hmm.log_likelihood_multi(signals).sum()
    print("LL for fitted params:", L)

    # test in actual params
    _, _, _, pi, A, R, mu, sigma = big_init()
    hmm.set(pi, A, R, mu, sigma)
    L = hmm.log_likelihood_multi(signals).sum()
    print("LL for actual params:", L)

if __name__ == '__main__':
    # real_signal()
    fake_signal()

max_stars_repo_name: ruslanmv/Machine-Learning-Codes

max_issues_repo_name: ruslanmv/Machine-Learning-Codes

max_forks_repo_name: ruslanmv/Machine-Learning-Codes

reason: import numpy,import theano,from theano


text:
import cv2
import numpy as np
import os, glob, shutil
import json
import re
import fnmatch
from PIL import Image
from pycococreatortools import pycococreatortools

def rgb2masks(label_img, save_dir):
    lbl_name = os.path.split(label_img)[-1].split('.')[0]
    lbl = cv2.imread(label_img, 1)
    h, w = lbl.shape[:2]
    lbl_dict = {}
    idx = 0
    white_mask = np.ones((h, w, 3), dtype=np.uint8) * 255
    for i in range(h):
        for j in range(w):
            if tuple(lbl[i][j]) in lbl_dict or \
                    tuple(lbl[i][j]) == (0, 0, 0) or \
                    tuple(lbl[i][j]) == (255, 255, 255):
                continue
            lbl_dict[tuple(lbl[i][j])] = idx
            mask = (lbl == lbl[i][j]).all(-1)
            binary_mask = np.where(mask[..., None], white_mask, 0)
            mask_name = lbl_name + '_' + str(idx) + '.png'
            cv2.imwrite(os.path.join(save_dir, mask_name), binary_mask)
            idx += 1
    # print(lbl_dict)

def filter_for_jpeg(root, files):
    file_types = ['*.jpeg', '*.jpg', '*.png']
    file_types = r'|'.join([fnmatch.translate(x) for x in file_types])
    files = [os.path.join(root, f) for f in files]
    files = [f for f in files if re.match(file_types, f)]
    return files

def filter_for_annotations(root, files, image_filename):
    file_types = ['*.png']
    file_types = r'|'.join([fnmatch.translate(x) for x in file_types])
    basename_no_extension = os.path.splitext(os.path.basename(image_filename))[0]
    file_name_prefix = basename_no_extension + '.*'
    files = [os.path.join(root, f) for f in files]
    files = [f for f in files if re.match(file_types, f)]
    files = [f for f in files if re.match(file_name_prefix, os.path.splitext(os.path.basename(f))[0])]
    return files

if __name__ == '__main__':
    root_dir = '/Users/dyy/Desktop/human_seg_datasets'
    img_dir = '/Users/dyy/Desktop/human/images'
    lbl_dir = '/Users/dyy/Desktop/human/labels'
    xml_dir = '/Users/dyy/Desktop/human/xmls'
    for dir in [img_dir, lbl_dir, xml_dir]:
        if not os.path.exists(dir):
            os.makedirs(dir)

    for root, dirs, files in os.walk(root_dir):
        for dir in dirs:
            if root.split('/')[-1].startswith('2019'):
                if dir == 'segmentationObject':
                    labels = glob.glob(os.path.join(root, dir, '*.png'))
                    for label in labels:
                        name = label.split('/')[-3][:12] + '_' + label.split('/')[-1]
                        shutil.copy(label, os.path.join(lbl_dir, name))

                if dir == 'JPEGimages':
                    imgs = glob.glob(os.path.join(root, dir, '*.jpg'))
                    for img in imgs:
                        name = img.split('/')[-3][:12] + '_' + img.split('/')[-1]
                        shutil.copy(img, os.path.join(img_dir, name))

                if dir == 'VOC':
                    xmls = glob.glob(os.path.join(root, dir, '*.xml'))
                    for xml in xmls:
                        name = xml.split('/')[-3][:12] + '_' + xml.split('/')[-1]
                        shutil.copy(xml, os.path.join(xml_dir, name))

    label_dir = '/Users/dyy/Desktop/human/labels/'
    save_dir = '/Users/dyy/Desktop/human/masks/'
    if not os.path.exists(save_dir):
        os.makedirs(save_dir)
    label_list = glob.glob(os.path.join(label_dir, '*.png'))[:1]
    for i, label_name in enumerate(label_list):
        print(i, label_name)
        rgb2masks(label_name, save_dir)

    ROOT_DIR = '/Users/dyy/Desktop/human'
    IMAGE_DIR = os.path.join(ROOT_DIR, "images")
    ANNOTATION_DIR = os.path.join(ROOT_DIR, "masks")

    CATEGORIES = [
        {
            'supercategory': 'none',
            'name': '0',
            'id': 0
        },
        {
            'supercategory': 'none',
            'name': '1',
            'id': 1
        },
        {
            'supercategory': 'none',
            'name': '2',
            'id': 2
        }
    ]

    coco_output = {
        "images": [],
        "annotations": [],
        "categories": CATEGORIES
    }

    image_id = 1
    segmentation_id = 1

    # filter for jpeg images
    for root, _, files in os.walk(IMAGE_DIR):
        image_files = filter_for_jpeg(root, files)

        # go through each image
        for image_filename in image_files:
            image = Image.open(image_filename)
            image_info = pycococreatortools.create_image_info(
                image_id, os.path.basename(image_filename), image.size)
            coco_output["images"].append(image_info)

            # filter for associated png annotations
            for root, _, files in os.walk(ANNOTATION_DIR):
                annotation_files = filter_for_annotations(root, files, image_filename)

                # go through each associated annotation
                for annotation_filename in annotation_files:

                    class_id = [x['id'] for x in CATEGORIES if x['name'] in annotation_filename][0]

                    category_info = {'id': class_id, 'is_crowd': 'crowd' in image_filename}
                    binary_mask = np.asarray(Image.open(annotation_filename)
                                             .convert('1')).astype(np.uint8)

                    annotation_info = pycococreatortools.create_annotation_info(
                        segmentation_id, image_id, category_info, binary_mask,
                        image.size, tolerance=2)

                    if annotation_info is not None:
                        coco_output["annotations"].append(annotation_info)

                    segmentation_id = segmentation_id + 1

            image_id = image_id + 1

    with open('{}/test.json'.format(ROOT_DIR), 'w') as output_json_file:
        json.dump(coco_output, output_json_file)

max_issues_repo_name: hukefei/SOLO-master

max_forks_repo_name: hukefei/SOLO-master

reason: import numpy


text:
import sys
import numpy as np

from mspasspy.ccore.seismic import Seismogram, SlownessVector
from mspasspy.ccore.utility import SphericalCoordinate
from mspasspy.algorithms.basic import ator, rtoa, rotate, rotate_to_standard, free_surface_transformation, transform, linear_taper, cosine_taper, vector_taper

# module to test
sys.path.append("python/tests")

from helper import (get_live_seismogram,
                    get_live_timeseries,
                    get_sin_timeseries)

def test_ator_rtoa():
    ts = get_live_timeseries()
    original_t0 = ts.t0
    ts_new = ator(ts, 1)
    assert ts_new.time_is_relative()
    assert ts_new.t0 == original_t0 - 1

    ts_new2 = rtoa(ts_new)
    assert ts_new2.time_is_UTC()
    assert ts_new2.t0 == original_t0

def test_rotate():
    seis = Seismogram()
    seis.npts = 100
    seis.t0 = 0.0
    seis.dt = 0.001
    seis.live = 1
    for i in range(3):
        for j in range(100):
            if i == 0:
                seis.data[i, j] = 1.0
            else:
                seis.data[i, j] = 0.0
    seis.data[0, 1] = 1.0
    seis.data[0, 2] = 1.0
    seis.data[0, 3] = 0.0
    seis.data[1, 1] = 1.0
    seis.data[1, 2] = 1.0
    seis.data[1, 3] = 0.0
    seis.data[2, 1] = 1.0
    seis.data[2, 2] = 0.0
    seis.data[2, 3] = 1.0

    sc = SphericalCoordinate()
    sc.phi = 0.0
    sc.theta = np.pi/4
    seis2 = rotate(seis, sc)
    assert all(np.isclose(seis2.data[:, 3], [0, -0.707107, 0.707107]))
    seis3 = rotate_to_standard(seis2)
    assert all(np.isclose(seis3.data[:, 3], [0, 0, 1]))

def test_transform():
    seis = Seismogram()
    seis.npts = 100
    seis.t0 = 0.0
    seis.dt = 0.001
    seis.live = 1
    for i in range(3):
        for j in range(100):
            if i == 0:
                seis.data[i, j] = 1.0
            else:
                seis.data[i, j] = 0.0
    seis.data[0, 1] = 1.0
    seis.data[0, 2] = 1.0
    seis.data[0, 3] = 0.0
    seis.data[1, 1] = 1.0
    seis.data[1, 2] = 1.0
    seis.data[1, 3] = 0.0
    seis.data[2, 1] = 1.0
    seis.data[2, 2] = 0.0
    seis.data[2, 3] = 1.0

    a = np.zeros((3, 3))
    a[0][0] = 1.0
    a[0][1] = 1.0
    a[0][2] = 1.0
    a[1][0] = -1.0
    a[1][1] = 1.0
    a[1][2] = 1.0
    a[2][0] = 0.0
    a[2][1] = -1.0
    a[2][2] = 0.0
    seis1 = transform(seis, a)
    assert all(np.isclose(seis1.data[:, 0], [1, -1,  0]))
    assert all(np.isclose(seis1.data[:, 1], [3,  1, -1]))
    assert all(np.isclose(seis1.data[:, 2], [2,  0, -1]))
    assert all(np.isclose(seis1.data[:, 3], [1,  1,  0]))
    seis2 = rotate_to_standard(seis1)
    assert all(np.isclose(seis2.data[:, 0], [1, 0, 0]))
    assert all(np.isclose(seis2.data[:, 1], [1, 1, 1]))
    assert all(np.isclose(seis2.data[:, 2], [1, 1, 0]))
    assert all(np.isclose(seis2.data[:, 3], [0, 0, 1]))

    uvec = SlownessVector()
    uvec.ux = 0.17085  # cos(-20deg)/5.5
    uvec.uy = -0.062185  # sin(-20deg)/5.5
    seis3 = free_surface_transformation(seis2, uvec, 5.0, 3.5)
    assert (np.isclose(seis3.tmatrix,
                       np.array([[-0.171012, -0.469846,  0],
                                 [0.115793, -0.0421458, 0.445447],
                                 [-0.597975,  0.217647,  0.228152]]))).all()

    # test with invalid uvec, but inplace return
    seis4 = free_surface_transformation(seis2, SlownessVector(1.0, 1.0, 0.0), 5.0, 3.5)
    assert seis4

def test_taper():
    ts = get_live_timeseries()
    ts.t0 = 0
    ts.dt = 1
    ts.npts = 200
    ts.data += 1

    ts_l = linear_taper(ts, 4, 14, 170, 180)
    assert ts_l.data[4] == 0
    assert ts_l.data[9] == 0.5
    assert ts_l.data[14] == 1
    assert ts_l.data[170] == 1
    assert ts_l.data[175] == 0.5
    assert ts_l.data[180] == 0

    ts.npts = 200
    ts.data += 1

    ts_c = cosine_taper(ts, 4, 14, 170, 180)
    assert ts_c.data[4] == 0
    assert ts_c.data[9] == 0.5
    assert ts_c.data[14] == 1
    assert ts_c.data[170] == 1
    assert ts_c.data[175] == 0.5
    assert ts_c.data[180] == 0

    ts.npts = 200
    ts.data += 1

    vtaper = np.zeros(200)
    vtaper += 0.5
    ts_v = vector_taper(ts, vtaper)
    assert ts_v.data[4] == 0.5
    assert ts_v.data[9] == 0.5
    assert ts_v.data[14] == 0.5
    assert ts_v.data[170] == 0.5
    assert ts_v.data[175] == 0.5
    assert ts_v.data[180] == 0.5

max_stars_repo_name: Yangzhengtang/mspass

max_issues_repo_name: Yangzhengtang/mspass

max_forks_repo_name: Yangzhengtang/mspass


text:
import cv2
import numpy as np

def bgr2rgb(image):
    return cv2.cvtColor(image, cv2.COLOR_BGR2RGB)

def flipimg(image):
    return cv2.flip(image, 1)

def cropimg(image):
    cropped = image[60:130, :]
    return cropped

def resize(image, shape=(160, 70)):
    return cv2.resize(image, shape)

def crop_and_resize(image):
    cropped = cropimg(image)
    resized = resize(cropped)
    return resized


text:
"""Tests for module utils for timing and parallel computation """

# Author: Remi Flamary <remi.flamary@unice.fr>
#
# License: MIT License

import ot
import numpy as np
import sys

def test_parmap():

    n = 10

    def f(i):
        return 1.0 * i * i

    a = np.arange(n)

    l1 = list(map(f, a))

    l2 = list(ot.utils.parmap(f, a))

    np.testing.assert_allclose(l1, l2)

def test_tic_toc():

    import time

    ot.tic()
    time.sleep(0.5)
    t = ot.toc()
    t2 = ot.toq()

    # test timing
    np.testing.assert_allclose(0.5, t, rtol=1e-1, atol=1e-1)

    # test toc vs toq
    np.testing.assert_allclose(t, t2, rtol=1e-1, atol=1e-1)

def test_kernel():

    n = 100

    x = np.random.randn(n, 2)

    K = ot.utils.kernel(x, x)

    # gaussian kernel  has ones on the diagonal
    np.testing.assert_allclose(np.diag(K), np.ones(n))

def test_unif():

    n = 100

    u = ot.unif(n)

    np.testing.assert_allclose(1, np.sum(u))

def test_dist():

    n = 100

    x = np.random.randn(n, 2)

    D = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            D[i, j] = np.sum(np.square(x[i, :] - x[j, :]))

    D2 = ot.dist(x, x)
    D3 = ot.dist(x)

    # dist shoul return squared euclidean
    np.testing.assert_allclose(D, D2)
    np.testing.assert_allclose(D, D3)

def test_dist0():

    n = 100
    M = ot.utils.dist0(n, method='lin_square')

    # dist0 default to linear sampling with quadratic loss
    np.testing.assert_allclose(M[0, -1], (n - 1) * (n - 1))

def test_dots():

    n1, n2, n3, n4 = 100, 50, 200, 100

    A = np.random.randn(n1, n2)
    B = np.random.randn(n2, n3)
    C = np.random.randn(n3, n4)

    X1 = ot.utils.dots(A, B, C)

    X2 = A.dot(B.dot(C))

    np.testing.assert_allclose(X1, X2)

def test_clean_zeros():

    n = 100
    nz = 50
    nz2 = 20
    u1 = ot.unif(n)
    u1[:nz] = 0
    u1 = u1 / u1.sum()
    u2 = ot.unif(n)
    u2[:nz2] = 0
    u2 = u2 / u2.sum()

    M = ot.utils.dist0(n)

    a, b, M2 = ot.utils.clean_zeros(u1, u2, M)

    assert len(a) == n - nz
    assert len(b) == n - nz2

def test_cost_normalization():

    C = np.random.rand(10, 10)

    # does nothing
    M0 = ot.utils.cost_normalization(C)
    np.testing.assert_allclose(C, M0)

    M = ot.utils.cost_normalization(C, 'median')
    np.testing.assert_allclose(np.median(M), 1)

    M = ot.utils.cost_normalization(C, 'max')
    np.testing.assert_allclose(M.max(), 1)

    M = ot.utils.cost_normalization(C, 'log')
    np.testing.assert_allclose(M.max(), np.log(1 + C).max())

    M = ot.utils.cost_normalization(C, 'loglog')
    np.testing.assert_allclose(M.max(), np.log(1 + np.log(1 + C)).max())

def test_check_params():

    res1 = ot.utils.check_params(first='OK', second=20)
    assert res1 is True

    res0 = ot.utils.check_params(first='OK', second=None)
    assert res0 is False

def test_deprecated_func():

    @ot.utils.deprecated('deprecated text for fun')
    def fun():
        pass

    def fun2():
        pass

    @ot.utils.deprecated('deprecated text for class')
    class Class():
        pass

    if sys.version_info < (3, 5):
        print('Not tested')
    else:
        assert ot.utils._is_deprecated(fun) is True

        assert ot.utils._is_deprecated(fun2) is False

def test_BaseEstimator():

    class Class(ot.utils.BaseEstimator):

        def __init__(self, first='spam', second='eggs'):

            self.first = first
            self.second = second

    cl = Class()

    names = cl._get_param_names()
    assert 'first' in names
    assert 'second' in names

    params = cl.get_params()
    assert 'first' in params
    assert 'second' in params

    params['first'] = 'spam again'
    cl.set_params(**params)

    assert cl.first == 'spam again'


text:
import numpy as np
import tensorflow as tf
from tensorflow.keras.initializers import glorot_uniform, zeros
from tensorflow.keras.optimizers import Adam
from tensorflow.keras.activations import softmax, relu
from tensorflow.keras.losses import SparseCategoricalCrossentropy

from graphadv import is_binary
from graphadv.attack.untargeted.untargeted_attacker import UntargetedAttacker
from graphadv.utils.surrogate_utils import train_a_surrogate
from graphadv.utils.graph_utils import likelihood_ratio_filter
from graphgallery.nn.models import DenseGCN
from graphgallery import tqdm, asintarr, normalize_adj_tensor, astensor

# cora lr=0.1, citeseer lr=0.01, lambda_=1. reaches best result
class BaseMeta(UntargetedAttacker):
    '''Base model for Mettack.'''

    def __init__(self, adj, x, labels,
                 idx_train, idx_unlabeled,
                 hidden_layers, use_relu, self_training_labels=None,
                 seed=None, name=None, device='CPU:0', **kwargs):

        super().__init__(adj=adj, x=x, labels=labels, seed=seed, name=name, device=device, **kwargs)
        adj, x, labels = self.adj, self.x, self.labels

        idx_train = asintarr(idx_train)
        idx_unlabeled = asintarr(idx_unlabeled)

        if self_training_labels is None:
            surrogate = DenseGCN(adj, x, labels, device='GPU', norm_x=None, seed=None)
            surrogate.build(16, activations='relu' if use_relu else None)
            his = surrogate.train(idx_train, verbose=0, epochs=200, save_best=False)
            self_training_labels = surrogate.predict(idx_unlabeled).argmax(1)

        self.ll_ratio = None
        # mettack can also conduct feature attack
        self.allow_feature_attack = True

        with tf.device(self.device):
            self.idx_train = astensor(idx_train, dtype=self.intx)
            self.idx_unlabeled = astensor(idx_unlabeled, dtype=self.intx)
            self.labels_train = astensor(self.labels[idx_train], dtype=self.floatx)
            self.self_training_labels = astensor(self_training_labels, dtype=self.floatx)
            self.tf_adj = astensor(adj.A, dtype=self.floatx)
            self.tf_x = astensor(x, dtype=self.floatx)
            self.build(hidden_layers=hidden_layers)
            self.use_relu = use_relu
            self.loss_fn = SparseCategoricalCrossentropy(from_logits=True)

            self.adj_changes = tf.Variable(tf.zeros_like(self.tf_adj))
            self.x_changes = tf.Variable(tf.zeros_like(self.tf_x))

    def reset(self):
        super().reset()
        self.structure_flips = []
        self.feature_flips = []

        with tf.device(self.device):
            self.adj_changes.assign(tf.zeros_like(self.tf_adj))
            self.x_changes.assign(tf.zeros_like(self.tf_x))

    def filter_potential_singletons(self, modified_adj):
        """
        Computes a mask for entries potentially leading to singleton nodes, i.e. one of the two nodes corresponding to
        the entry have degree 1 and there is an edge between the two nodes.
        Returns
        -------
        tf.Tensor shape [N, N], float with ones everywhere except the entries of potential singleton nodes,
        where the returned tensor has value 0.
        """
        N = self.n_nodes
        degrees = tf.reduce_sum(modified_adj, axis=1)
        degree_one = tf.equal(degrees, 1)
        resh = tf.reshape(tf.tile(degree_one, [N]), [N, N])
        l_and = tf.logical_and(resh, tf.equal(modified_adj, 1))
        logical_and_symmetric = tf.logical_or(l_and, tf.transpose(l_and))
        flat_mask = 1. - tf.cast(logical_and_symmetric, self.floatx)
        return flat_mask

    def log_likelihood_constraint(self, adj, modified_adj, ll_cutoff):
        """
        Computes a mask for entries that, if the edge corresponding to the entry is added/removed, would lead to the
        log likelihood constraint to be violated.

        """
        t_d_min = tf.constant(2., dtype=self.floatx)
        t_possible_edges = tf.constant(np.array(np.triu(np.ones([self.n_nodes, self.n_nodes]), k=1).nonzero()).T,
                                       dtype=self.intx)
        allowed_mask, current_ratio = likelihood_ratio_filter(t_possible_edges, modified_adj,
                                                              adj, t_d_min, ll_cutoff)

        return allowed_mask, current_ratio

    @tf.function
    def get_perturbed_adj(self, adj, adj_changes):
        adj_changes_square = adj_changes - tf.linalg.band_part(adj_changes, 0, 0)
#         adj_changes_sym = adj_changes_square + tf.transpose(adj_changes_square)
        clipped_adj_changes = self.clip(adj_changes_square)
        return adj + clipped_adj_changes

    @tf.function
    def get_perturbed_x(self, x, x_changes):
        return x + self.clip(x_changes)

    def forward(self, x, adj):
        h = x
        for w in self.weights:
            h = adj @ h @ w
            if self.use_relu:
                h = relu(h)
        return h

    def structure_score(self, modified_adj, adj_grad, ll_constraint=None, ll_cutoff=None):
        adj_meta_grad = adj_grad * (-2. * modified_adj + 1.)
        # Make sure that the minimum entry is 0.
        adj_meta_grad -= tf.reduce_min(adj_meta_grad)

        if not self.allow_singleton:
            # Set entries to 0 that could lead to singleton nodes.
            singleton_mask = self.filter_potential_singletons(modified_adj)
            adj_meta_grad *= singleton_mask

        if ll_constraint:
            allowed_mask, self.ll_ratio = self.log_likelihood_constraint(modified_adj, self.tf_adj, ll_cutoff)
            adj_meta_grad = adj_meta_grad * allowed_mask

        return tf.reshape(adj_meta_grad, [-1])

    def feature_score(self, modified_x, x_grad):
        x_meta_grad = x_grad * (-2. * modified_x + 1.)
        x_meta_grad -= tf.reduce_min(x_meta_grad)
        return tf.reshape(x_meta_grad, [-1])

    def clip(self, matrix):
        clipped_matrix = tf.clip_by_value(matrix, -1., 1.)
        return clipped_matrix

class Metattack(BaseMeta):

    def __init__(self,  adj, x, labels,
                 idx_train, idx_unlabeled,
                 lr=0.1, epochs=100,
                 momentum=0.9, lambda_=0.,
                 hidden_layers=[16], use_relu=True, self_training_labels=None,
                 seed=None, name=None, device='CPU:0', **kwargs):

        super().__init__(adj, x, labels,
                         idx_train, idx_unlabeled,
                         hidden_layers=hidden_layers, use_relu=use_relu,
                         self_training_labels=self_training_labels,
                         seed=seed, name=name, device=device, **kwargs)

        self.lr = lr
        self.epochs = epochs
        self.momentum = momentum
        self.lambda_ = lambda_

        if lambda_ not in (0., 0.5, 1.):
            raise ValueError('Invalid value of `lanbda_`, allowed values [0: (meta-self), 1: (meta-train), 0.5: (meta-both)].')

    def build(self, hidden_layers):

        weights, velocities = [], []
        zeros_initializer = zeros()

        pre_hid = self.n_attrs
        for hid in hidden_layers + [self.n_classes]:
            shape = (pre_hid, hid)
            # use zeros_initializer temporary to save time
            weight = tf.Variable(zeros_initializer(shape=shape, dtype=self.floatx))
            w_velocity = tf.Variable(zeros_initializer(shape=shape, dtype=self.floatx))

            weights.append(weight)
            velocities.append(w_velocity)

            pre_hid = hid

        self.weights, self.velocities = weights, velocities

    def initialize(self):
        w_initializer = glorot_uniform()
        zeros_initializer = zeros()

        for w, wv in zip(self.weights, self.velocities):
            w.assign(w_initializer(w.shape, dtype=self.floatx))
            wv.assign(zeros_initializer(wv.shape, dtype=self.floatx))

    @tf.function
    def train_step(self, x, adj, index, labels):
        with tf.GradientTape() as tape:
            output = self.forward(x, adj)
            logit = tf.gather(output, index)
            loss = self.loss_fn(labels, logit)

        weight_grads = tape.gradient(loss, self.weights)
        return weight_grads

    def inner_train(self, adj, x):

        self.initialize()
        adj_norm = normalize_adj_tensor(adj)

        for _ in range(self.epochs):
            weight_grads = self.train_step(x, adj_norm, self.idx_train, self.labels_train)

            for v, g in zip(self.velocities, weight_grads):
                v.assign(self.momentum * v + g)

            for w, v in zip(self.weights, self.velocities):
                w.assign_sub(self.lr * v)

    @tf.function
    def meta_grad(self):

        modified_adj, modified_x = self.tf_adj, self.tf_x
        persistent = self.structure_attack and self.feature_attack

        with tf.GradientTape(persistent=persistent) as tape:
            if self.structure_attack:
                modified_adj = self.get_perturbed_adj(self.tf_adj, self.adj_changes)

            if self.feature_attack:
                modified_x = self.get_perturbed_x(self.tf_x, self.x_changes)

            adj_norm = normalize_adj_tensor(modified_adj)
            output = self.forward(modified_x, adj_norm) / 5.0
            logit_labeled = tf.gather(output, self.idx_train)
            logit_unlabeled = tf.gather(output, self.idx_unlabeled)

            loss_labeled = self.loss_fn(self.labels_train, logit_labeled)
            loss_unlabeled = self.loss_fn(self.self_training_labels, logit_unlabeled)

            attack_loss = self.lambda_ * loss_labeled + (1 - self.lambda_) * loss_unlabeled

        adj_grad, x_grad = None, None

        if self.structure_attack:
            adj_grad = tape.gradient(attack_loss, self.adj_changes)

        if self.feature_attack:
            x_grad = tape.gradient(attack_loss, self.x_changes)

        return adj_grad, x_grad

    def attack(self, n_perturbations=0.05, structure_attack=True, feature_attack=False,
               ll_constraint=False, ll_cutoff=0.004, disable=False):

        super().attack(n_perturbations, structure_attack, feature_attack)

        if ll_constraint:
            raise NotImplementedError('`log_likelihood_constraint` has not been well tested.'
                                      ' Please set `ll_constraint=False` to achieve a better performance.')

        if feature_attack and not is_binary(self.x):
            raise ValueError("Attacks on the node features are currently only supported for binary attributes.")

        with tf.device(self.device):
            modified_adj, modified_x = self.tf_adj, self.tf_x
            adj_changes, x_changes = self.adj_changes, self.x_changes
            structure_flips, feature_flips = self.structure_flips, self.feature_flips

            for _ in tqdm(range(self.n_perturbations), desc='Peturbing Graph', disable=disable):

                if structure_attack:
                    modified_adj = self.get_perturbed_adj(self.tf_adj, adj_changes)

                if feature_attack:
                    modified_x = self.get_perturbed_x(self.tf_x, x_changes)

                self.inner_train(modified_adj, modified_x)

                adj_grad, x_grad = self.meta_grad()

                adj_meta_score = tf.constant(0.0)
                x_meta_score = tf.constant(0.0)

                if structure_attack:
                    adj_meta_score = self.structure_score(modified_adj, adj_grad, ll_constraint, ll_cutoff)

                if feature_attack:
                    x_meta_score = self.feature_score(modified_x, x_grad)

                if tf.reduce_max(adj_meta_score) >= tf.reduce_max(x_meta_score):
                    adj_meta_argmax = tf.argmax(adj_meta_score)
                    row, col = divmod(adj_meta_argmax.numpy(), self.n_nodes)
                    adj_changes[row, col].assign(-2. * modified_adj[row, col] + 1.)
                    adj_changes[col, row].assign(-2. * modified_adj[col, row] + 1.)
                    structure_flips.append((row, col))
                else:
                    x_meta_argmax = tf.argmax(x_meta_score)
                    row, col = divmod(x_meta_argmax.numpy(), self.n_attrs)
                    x_changes[row, col].assign(-2 * modified_x[row, col] + 1)
                    feature_flips.append((row, col))

class MetaApprox(BaseMeta):

    def __init__(self,  adj, x, labels,
                 idx_train, idx_unlabeled,
                 lr=0.1, epochs=100, lambda_=0.,
                 hidden_layers=[16], use_relu=True, self_training_labels=None,
                 seed=None, name=None, device='CPU:0', **kwargs):

        self.lr = lr
        self.epochs = epochs
        self.lambda_ = lambda_

        if lambda_ not in (0., 0.5, 1.):
            raise ValueError('Invalid value of `lanbda_`, allowed values [0: (meta-self), 1: (meta-train), 0.5: (meta-both)].')

        super().__init__(adj, x, labels,
                         idx_train, idx_unlabeled,
                         hidden_layers=hidden_layers, use_relu=use_relu,
                         self_training_labels=self_training_labels,
                         seed=seed, name=name, device=device, **kwargs)

    def build(self, hidden_layers):

        weights = []
        zeros_initializer = zeros()

        pre_hid = self.n_attrs
        for hid in hidden_layers + [self.n_classes]:
            shape = (pre_hid, hid)
            # use zeros_initializer temporary to save time
            weight = tf.Variable(zeros_initializer(shape=shape, dtype=self.floatx))
            weights.append(weight)
            pre_hid = hid

        self.weights = weights
        self.adj_grad_sum = tf.Variable(tf.zeros_like(self.tf_adj))
        self.x_grad_sum = tf.Variable(tf.zeros_like(self.tf_x))
        self.optimizer = Adam(self.lr, epsilon=1e-8)

    def initialize(self):

        w_initializer = glorot_uniform()
        zeros_initializer = zeros()

        for w in self.weights:
            w.assign(w_initializer(w.shape, dtype=self.floatx))

        if self.structure_attack:
            self.adj_grad_sum.assign(zeros_initializer(self.adj_grad_sum.shape, dtype=self.floatx))

        if self.feature_attack:
            self.x_grad_sum.assign(zeros_initializer(self.x_grad_sum.shape, dtype=self.floatx))

        # reset optimizer
        for var in self.optimizer.variables():
            var.assign(tf.zeros_like(var))

    @tf.function
    def meta_grad(self):
        self.initialize()

        modified_adj, modified_x = self.tf_adj, self.tf_x
        adj_grad_sum, x_grad_sum = self.adj_grad_sum, self.x_grad_sum
        optimizer = self.optimizer

        for _ in tf.range(self.epochs):

            with tf.GradientTape(persistent=True) as tape:
                if self.structure_attack:
                    modified_adj = self.get_perturbed_adj(self.tf_adj, self.adj_changes)

                if self.feature_attack:
                    modified_x = self.get_perturbed_x(self.tf_x, self.x_changes)

                adj_norm = normalize_adj_tensor(modified_adj)
                output = self.forward(modified_x, adj_norm) / 5.0
                logit_labeled = tf.gather(output, self.idx_train)
                logit_unlabeled = tf.gather(output, self.idx_unlabeled)

                loss_labeled = self.loss_fn(self.labels_train, logit_labeled)
                loss_unlabeled = self.loss_fn(self.self_training_labels, logit_unlabeled)

                attack_loss = self.lambda_ * loss_labeled + (1 - self.lambda_) * loss_unlabeled

            adj_grad, x_grad = None, None

            gradients = tape.gradient(loss_labeled, self.weights)
            optimizer.apply_gradients(zip(gradients, self.weights))

            if self.structure_attack:
                adj_grad = tape.gradient(attack_loss, self.adj_changes)
                adj_grad_sum.assign_add(adj_grad)

            if self.feature_attack:
                x_grad = tape.gradient(attack_loss, self.x_changes)
                x_grad_sum.assign_add(x_grad)

            del tape

        return adj_grad_sum, x_grad_sum

    def attack(self, n_perturbations=0.05, structure_attack=True, feature_attack=False,
               ll_constraint=False, ll_cutoff=0.004, disable=False):

        super().attack(n_perturbations, structure_attack, feature_attack)

        if ll_constraint:
            raise NotImplementedError('`log_likelihood_constraint` has not been well tested.'
                                      ' Please set `ll_constraint=False` to achieve a better performance.')

        if feature_attack and not is_binary(self.x):
            raise ValueError("Attacks on the node features are currently only supported for binary attributes.")

        with tf.device(self.device):
            modified_adj, modified_x = self.tf_adj, self.tf_x
            adj_changes, x_changes = self.adj_changes, self.x_changes
            structure_flips, feature_flips = self.structure_flips, self.feature_flips

            for _ in tqdm(range(self.n_perturbations), desc='Peturbing Graph', disable=disable):

                adj_grad, x_grad = self.meta_grad()

                adj_meta_score = tf.constant(0.0)
                x_meta_score = tf.constant(0.0)

                if structure_attack:
                    modified_adj = self.get_perturbed_adj(self.tf_adj, adj_changes)
                    adj_meta_score = self.structure_score(modified_adj, adj_grad, ll_constraint, ll_cutoff)

                if feature_attack:
                    modified_x = self.get_perturbed_x(self.tf_x, x_changes)
                    x_meta_score = self.feature_score(modified_x, feature_grad)

                if tf.reduce_max(adj_meta_score) >= tf.reduce_max(x_meta_score):
                    adj_meta_argmax = tf.argmax(adj_meta_score)
                    row, col = divmod(adj_meta_argmax.numpy(), self.n_nodes)
                    adj_changes[row, col].assign(-2. * modified_adj[row, col] + 1.)
                    adj_changes[col, row].assign(-2. * modified_adj[col, row] + 1.)
                    structure_flips.append((row, col))
                else:
                    x_meta_argmax = tf.argmax(x_meta_score)
                    row, col = divmod(x_meta_argmax.numpy(), self.n_attrs)
                    x_changes[row, col].assign(-2 * modified_x[row, col] + 1)
                    feature_flips.append((row, col))

max_stars_repo_name: EdisonLeeeee/graphadv

max_issues_repo_name: EdisonLeeeee/graphadv

max_forks_repo_name: EdisonLeeeee/graphadv


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

# ### Calculating Null Distribution
#
# - Null distribution - is generated by getting the median correlation score of randomly combined compounds that do not share/come from the same MOAs.
#
#
# ### The goal here:
# - is to compute the p-value for each MOA per dose by evaluating the probability of random combinations of compounds (from different MOAs) having greater median correlation score than compounds of the same MOA.
#
#
#
#
# - In our case, we generated 1000 median correlation scores from randomly combined compounds as the **null distribution** for each MOA_SIZE class ***i.e. for a moa_size class - we have 1000 medians scores from randomly combined compounds of different MOAs.***
#
#
#
#
# - Moa_size is the number of compounds in a specific MOA and moa_size class is a specific group of MOAs that have the same number of compounds ***e.g all MOAs with just 2 compounds in them are in the same moa_size class.***
#
#
# ### Note:
#
# To generate the null distribution for all the four consensus data, you will have to execute this notebook four times, one for each.

# In[1]:

import os
import pathlib
import pandas as pd
import numpy as np
from collections import defaultdict
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
import seaborn as sns
from pycytominer import feature_select
from statistics import median
import random
sns.set_style("darkgrid")
from scipy import stats
import pickle
from io import BytesIO
from urllib.request import urlopen
from zipfile import ZipFile

# In[2]:

np.random.seed(42)

# #### - Load in the datasets required,
#
# - They were generated from the `cell_painting_moa_median_scores_calculation notebook`

# In[3]:

data_moa = pd.read_csv(os.path.join('moa_sizes_consensus_datasets', 'modz_consensus_data.csv'))
data_moa_vals = pd.read_csv(os.path.join('moa_sizes_consensus_datasets', 'modz_moa_median_scores.csv'))
data_moa_cpds = pd.read_csv(os.path.join('moa_sizes_consensus_datasets', 'cellpainting_moa_compounds.csv'))

# In[4]:

data_moa.shape

# In[5]:

data_moa_vals.shape

# In[6]:

data_moa_cpds.shape

# In[7]:

def conv_cols_to_list(df_moa_cpds):
    """This function convert string values in compound dataframe to lists"""

    moa_cpd_cols = [col for col in df_moa_cpds.columns.tolist()
                 if (col.startswith("moa_cpds_"))]
    for col in moa_cpd_cols:
        df_moa_cpds[col] = df_moa_cpds[col].apply(lambda row: row.split(';'))
    return df_moa_cpds

# In[8]:

data_moa_cpds = conv_cols_to_list(data_moa_cpds)

# In[9]:

def get_cpd_agg(data_moa, dose_number):
    """
    This function aggregate values for a particular
    dose by taking the mean value of distinct compounds in the dose
    """

    df_dose = data_moa[data_moa['Metadata_dose_recode'] == dose_number].copy()
    meta_cols = [col for col in df_dose.columns.tolist()
                 if (col.startswith("Metadata_"))]
    df_dose.drop(meta_cols, axis = 1, inplace = True)
    df_compound_agg = df_dose.groupby(['pert_iname']).agg(['mean'])
    df_compound_agg.columns  = df_compound_agg.columns.droplevel(1)
    df_compound_agg.rename_axis(None, axis=0, inplace = True)

    return df_compound_agg

# In[10]:

def cpds_found_in_all_doses(data_moa):
    """This function return a list of compounds found in all doses (1 - 6)"""
    cpds_fd = []
    for num in range(1,7):
        df_cpd_agg = get_cpd_agg(data_moa, num)
        all_cpds = df_cpd_agg.index.tolist()
        cpds_fd.append(all_cpds)

    cpds_fd_in_all = [cpd for list_cpds in cpds_fd
                      for cpd in list_cpds
                      if all(cpd in list_of_cpds for list_of_cpds in cpds_fd)]
    cpds_fd_in_all = list(set(cpds_fd_in_all))

    return cpds_fd_in_all

# In[11]:

cpds_fd_in_all = cpds_found_in_all_doses(data_moa)

# In[12]:

len(cpds_fd_in_all)

# In[13]:

all_moa_list = data_moa['moa'].unique().tolist()

# In[14]:

len(all_moa_list)

# In[15]:

all_moa_dict = {moa: [cpd for cpd in data_moa['pert_iname'][data_moa['moa']== moa].unique().tolist()
                      if cpd in cpds_fd_in_all]
                for moa in all_moa_list}
all_moa_dict = {kys:all_moa_dict[kys] for kys in all_moa_dict if all_moa_dict[kys]}

# In[16]:

len(all_moa_dict)

# In[17]:

def generate_moa_size_dict(df_moa_cpds):
    """
    Generates a dictionary with distinct moa_sizes
    (moa_size == number of compounds that is present in each MOA)
    as the keys and all compounds of MOAs with that particular size as the values
    """
    moa_size_dict = {}
    for size in df_moa_cpds['moa_size'].unique():
        size_df = df_moa_cpds[df_moa_cpds['moa_size'] == size].drop(['moa_size', 'moa'], axis = 1)
        size_df_values = size_df.values.tolist()
        size_df_vals = list(set([cpd for size_list in size_df_values for sizes in size_list for cpd in sizes]))
        moa_size_dict[size] = size_df_vals
    return moa_size_dict

# In[18]:

moa_sizes_dict = generate_moa_size_dict(data_moa_cpds)

# In[19]:

len(data_moa_cpds['moa_size'].unique())

# In[20]:

len(moa_sizes_dict)

# In[21]:

def check_similar_cpds(cpds, moa_dict):
    """This function checks if two compounds are found in the same moa"""

    for x in range(len(cpds)):
        for y in range(x+1, len(cpds)):
            for kys in moa_dict:
                if all(i in moa_dict[kys] for i in [cpds[x], cpds[y]]):
                    return True
    return False

# In[22]:

def get_random_cpds(all_cpds, moa_size, moa_cpds, all_moa_cpds):
    """
    This function return a list of random cpds that are not of the same moas
    or found in the current moa cpd's list
    """
    while (True):
        random_cpds = random.sample(all_cpds, moa_size)
        if not (any(cpds in moa_cpds for cpds in random_cpds) & (check_similar_cpds(random_cpds, all_moa_cpds))):
            break
    return random_cpds

# #### - You only need to generate the null distribution once (i.e. you can re-use the pickled null distribution for other consensus data), since the 1000 lists of randomly generated compounds combinations  for each MOA are found in all doses and all consensus datasets

# In[23]:

def get_null_distribution_cpds(moa_size_dict, cpds_list, all_moa_dict, rand_num = 1000):

    """
    This function returns the null distribution dictionary, with MOAs as the keys and
    1000 lists of randomly selected compounds combinations as the values for each MOA
    """
    null_distribution_moa = {}
    for size in moa_size_dict:
        moa_cpds = moa_size_dict[size]
        moa_cpds_list = []
        for idx in range(rand_num):
            start_again = True
            while (start_again):
                rand_cpds = get_random_cpds(cpds_list, size, moa_cpds, all_moa_dict)
                if rand_cpds not in moa_cpds_list:
                    start_again = False
            moa_cpds_list.append(rand_cpds)
        null_distribution_moa[size] = moa_cpds_list

    return null_distribution_moa

# In[24]:

null_distribution_moa = get_null_distribution_cpds(moa_sizes_dict, cpds_fd_in_all, all_moa_dict)

# In[25]:

#save the null_distribution_moa to pickle, you only need to run the code once
with open(os.path.join('moa_sizes_consensus_datasets', 'null_distribution.pickle'), 'wb') as handle:
    pickle.dump(null_distribution_moa, handle, protocol=pickle.HIGHEST_PROTOCOL)

# In[26]:

##load the null_distribution_moa from pickle
with open(os.path.join('moa_sizes_consensus_datasets', 'null_distribution.pickle'), 'rb') as handle:
    null_distribution_moa = pickle.load(handle)

# In[27]:

print('moa_size', '\tnumber of generated lists of randomly combined compounds')
for keys in null_distribution_moa:
    print(keys, '\t\t', len(null_distribution_moa[keys]))

# In[28]:

def assert_null_distribution(null_distribution_moa):

    """
    This function assert that each of the list in the 1000 lists of
    random compounds combination for each MOA are distinct with no duplicates
    """

    duplicates_moa = {}
    for keys in null_distribution_moa:
        null_dist = null_distribution_moa[keys]
        for cpds_moa in null_dist:
            cpds_duplicates = []
            new_list = list(filter(lambda cpds_list: cpds_list != cpds_moa, null_dist))
            if (len(new_list) != len(null_dist) - 1):
                cpds_duplicates.append(cpds_moa)
        if cpds_duplicates:
            duplicates_moa[keys] = cpds_duplicates
    return duplicates_moa

# In[29]:

duplicates_cpds_list = assert_null_distribution(null_distribution_moa)

# In[30]:

duplicates_cpds_list ##no duplicate found

# In[31]:

def calc_null_dist_median_scores(data_moa, dose_num, moa_cpds_list):
    """
    This function calculate the median of the correlation
    values for each of the list in the 1000 lists of
    random compounds combination for each MOA
    """
    df_cpd_agg = get_cpd_agg(data_moa, dose_num)
    median_corr_list = []
    for list_of_cpds in moa_cpds_list:
        df_cpds = df_cpd_agg.loc[list_of_cpds]
        cpds_corr = df_cpds.T.corr(method = 'spearman').values
        median_corr_val = median(list(cpds_corr[np.triu_indices(len(cpds_corr), k = 1)]))
        median_corr_list.append(median_corr_val)
    return median_corr_list

# **A P value can be computed nonparametrically by evaluating the probability of random compounds of different MOAs having greater median similarity value than compounds of the same MOAs.**

# In[32]:

def get_p_value(median_scores_list, df_moa_values, dose_name, moa_name):
    """
    This function calculate the p-value from the
    null_distribution median scores for each MOA
    """
    actual_med = df_moa_values.loc[moa_name, dose_name]
    p_value = np.sum(median_scores_list >= actual_med) / len(median_scores_list)
    return p_value

# In[33]:

def get_null_dist_median_scores(null_distribution_moa, df_moa):
    """
    This function calculate the median correlation scores for all
    1000 lists of randomly combined compounds for each moa_size class
    across all doses (1-6)
    """
    null_distribution_medians = {}
    for key in null_distribution_moa:
        median_score_list = []
        for num in range(1,7):
            moa_size_median_scores = calc_null_dist_median_scores(df_moa, num, null_distribution_moa[key])
            median_score_list.append(moa_size_median_scores)
        null_distribution_medians[key] = median_score_list
    return null_distribution_medians

# In[34]:

null_distribution_medns = get_null_dist_median_scores(null_distribution_moa, data_moa)

# In[35]:

def transform_dataframe(rep, rep_name):
    """
    Transforms replicate correlation dataframe to have 3 columns:
    dose, correlation_values and type of replicates
    """
    df_reps = pd.DataFrame.from_dict(rep, orient='index').T
    rep_melt = df_reps.melt(var_name="dose", value_name="correlation_values")
    rep_melt['type'] = rep_name
    return rep_melt

# In[36]:

threshold_df = []
for n_replicate in null_distribution_medns.keys():
    matched_null = pd.DataFrame(null_distribution_medns[n_replicate])
    for dose, dose_row in matched_null.iterrows():
        thresh = dose_row.quantile(0.95)
        dose_id = dose + 1
        threshold_df.append([n_replicate, dose_id, thresh])

threshold_df = pd.DataFrame(threshold_df, columns=["n_replicates", "dose", "95th_threshold"])
threshold_df

# In[37]:

def get_moa_p_vals(null_dist_median, df_moa_values):
    """
    This function returns a dict, with MOAs as the keys and the MOA's
    p-values for each dose (1-6) as the values
    """
    null_p_vals = {}
    df_moa_values = df_moa_values.set_index('moa').rename_axis(None, axis=0)
    for key in null_dist_median:
        df_moa_size = df_moa_values[df_moa_values['moa_size'] == key]
        for moa in df_moa_size.index:
            dose_p_values = []
            for num in range(1,7):
                dose_name = 'dose_' + str(num)
                moa_p_value = get_p_value(null_dist_median[key][num-1], df_moa_size, dose_name, moa)
                dose_p_values.append(moa_p_value)
            null_p_vals[moa] = dose_p_values
    sorted_null_p_vals = {key:value for key, value in sorted(null_p_vals.items(), key=lambda item: item[0])}
    return sorted_null_p_vals

# In[38]:

null_p_vals = get_moa_p_vals(null_distribution_medns, data_moa_vals)

# In[39]:

df_null_p_vals = pd.DataFrame.from_dict(null_p_vals, orient='index',
                                        columns = ['dose_' + str(x)
                                                   for x in range(1,7)]).reset_index().rename(columns={"index": "moa"})

# In[40]:

df_null_p_vals['moa_size'] = data_moa_vals['moa_size']

# In[41]:

df_null_p_vals.shape

# In[42]:

df_null_p_vals.head(10)

# In[43]:

def save_to_csv(df, path, file_name):
    """saves moa dataframes to csv"""

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

    df.to_csv(os.path.join(path, file_name), index=False)

# In[44]:

save_to_csv(df_null_p_vals, 'moa_sizes_consensus_datasets', 'modz_null_p_values.csv')

# In[45]:

# Output files for visualization
cpd_summary_file = pathlib.Path("moa_sizes_consensus_datasets/matching_score_per_MOA_CellPainting_p_values_compared_to_nonparametric_null.tsv.gz")

dose_recode_info = {
    'dose_1': '0.04 uM', 'dose_2':'0.12 uM', 'dose_3':'0.37 uM',
    'dose_4': '1.11 uM', 'dose_5':'3.33 uM', 'dose_6':'10 uM'
}

cpd_score_summary_df = (
    df_null_p_vals
    .rename(columns={"moa_size": "no_of_replicates"})
    .melt(
        id_vars=["moa", "no_of_replicates"],
        value_vars=["dose_1", "dose_2", "dose_3", "dose_4", "dose_5", "dose_6"],
        var_name="dose",
        value_name="p_value"
    )
)

cpd_score_summary_df.dose = cpd_score_summary_df.dose.replace(dose_recode_info)

cpd_score_summary_df.to_csv(cpd_summary_file, sep="\t", index=False)
cpd_score_summary_df.head()

reason: import numpy,from scipy