text:
import os
os.environ["CUDA_DEVICE_ORDER"]="PCI_BUS_ID"
os.environ["CUDA_VISIBLE_DEVICES"]="0"

import tensorflow as tf
from keras import backend as K
K.set_image_dim_ordering('th')

import numpy as np
import pandas as pd
import cv2
import zarr
import glob
import time
import matplotlib.pyplot as plt

from keras.models import Sequential, load_model, Model
from keras.layers import Dense, Dropout, Activation, Flatten
from keras.layers import Input
from keras.optimizers import Adam, SGD, RMSprop, Nadam

from keras.layers.convolutional import Convolution2D, MaxPooling2D
from keras.layers.advanced_activations import PReLU
from keras.layers import BatchNormalization, GlobalAveragePooling2D, GlobalMaxPooling2D
from keras.callbacks import EarlyStopping, ModelCheckpoint
from keras.layers.core import SpatialDropout3D
from keras.utils.np_utils import to_categorical

# In[2]:

zarr_dir = '/home/w/DS_Projects/Kaggle/DS Bowl 2017/input_data/stage1_zarr/DSBowl.zarr/'
zarr_store = zarr.DirectoryStore(zarr_dir)
zarr_load_group = zarr.hierarchy.open_group(store=zarr_store, mode='r')
dsb_pats = os.listdir('/home/w/DS_Projects/Kaggle/DS Bowl 2017/input_data/stage1_zarr/DSBowl.zarr/nodules_candidates/')

def load_zarr(patient_id):
    lung_cand_zarr = zarr_load_group['nodules_candidates'][patient_id]
    return np.array(lung_cand_zarr).astype('float32')

def load_data(start, end):
    print('Loading 2D full-size candidates.')
    df = pd.read_csv('/home/w/DS_Projects/Kaggle/DS Bowl 2017/input_data/stage1_labels_full2.csv')[start:end]
    df = df[df['id'].isin(dsb_pats)]
    t = time.time()
    masks = np.zeros((0, 1, 512, 512))
    labels = np.zeros(0)
    for i in range(len(df)):
        mask = load_zarr('{}'.format(df.iloc[i, 0]))
        mask[mask <= 1.0] = 0.
        z_nonzero = np.unique(np.nonzero(mask)[0])
        mask = mask[z_nonzero[0]:z_nonzero[-1], :, :]
        print('Nonzero mask shape:', mask.shape[0])
        mask = mask/255.
        masks = np.concatenate((masks, mask), 0)
        if df.iloc[i, 1] == 1:
            label = np.ones(mask.shape[0])
        if df.iloc[i, 1] == 0:
            label = np.zeros(mask.shape[0])
        labels = np.concatenate((labels, label), 0)
    print('Data shape:', masks.shape)
    print('Time it took to load the data:', time.time() - t)
    return masks, df, labels

def check_shapes(start, end):
    print('Loading 2D full-size candidates.')
    df = pd.read_csv('/home/w/DS_Projects/Kaggle/DS Bowl 2017/input_data/stage1_labels_full2.csv')[start:end]
    df = df[df['id'].isin(dsb_pats)]
    full_shape = 0
    for i in range(len(df)):
        mask = load_zarr('{}'.format(df.iloc[i, 0]))
        full_shape += mask.shape[0]
    return full_shape

def generate_train(start, end, batch_size, shuffle = True):
    step = 10
    subsets = np.arange(start, end, step)
    while True:
        for s in range(len(subsets)-1):
            print(subsets[s], subsets[s+1])
            masks, _, labels = load_data(subsets[s], subsets[s+1])
            y = to_categorical(labels, 2)
            #y = y.reshape((-1,1))
            sample_index = np.arange(0, masks.shape[0])

            if shuffle:
                print('Shuffling data')
                perm = np.random.permutation(len(sample_index))
                masks = masks[perm]
                y = y[perm]
            if len(sample_index) % batch_size != 0:
                mod = len(sample_index) % batch_size
                sample_index = sample_index[:-mod]
                assert len(sample_index) % batch_size == 0
            for j in range(len(sample_index)//batch_size):
                batch_index = sample_index[j*batch_size:j*batch_size+batch_size]
                lung_batch = masks[batch_index, :, :, :]
                y_batch = y[batch_index]
                yield(lung_batch, y_batch)

def generate_val(start, end, batch_size):
    step = 10
    subsets = np.arange(start, end, step)
    while True:
        for s in range(len(subsets)-1):
            print(subsets[s], subsets[s+1])
            masks, _, labels = load_data(subsets[s], subsets[s+1])
            y = to_categorical(labels, 2)
            #y = y.reshape((-1,1))
            sample_index = np.arange(0, masks.shape[0])

            if len(sample_index) % batch_size != 0:
                mod = len(sample_index) % batch_size
                sample_index = sample_index[:-mod]
                assert len(sample_index) % batch_size == 0
            for j in range(len(sample_index)//batch_size):
                batch_index = sample_index[j*batch_size:j*batch_size+batch_size]
                lung_batch = masks[batch_index, :, :, :]
                y_batch = y[batch_index]
                yield(lung_batch, y_batch)

def cnn2d():
    width = 4
    inputs = Input((1, 512, 512))
    conv1 = Convolution2D(width*2, 3, 3, activation='relu', border_mode='same')(inputs)
    #conv1 = BatchNormalization(axis = 1)(conv1)
    conv1 = Dropout(0.2)(conv1)
    conv1 = Convolution2D(width*2, 3, 3, activation='relu', border_mode='same')(conv1)
    #conv1 = BatchNormalization(axis = 1)(conv1)
    pool1 = MaxPooling2D(pool_size=(2, 2))(conv1)

    conv2 = Convolution2D(width*4, 3, 3, activation='relu', border_mode='same')(pool1)
    #conv2 = BatchNormalization(axis = 1)(conv2)
    conv2 = Dropout(0.2)(conv2)
    conv2 = Convolution2D(width*4, 3, 3, activation='relu', border_mode='same')(conv2)
    #conv2 = BatchNormalization(axis = 1)(conv2)
    pool2 = MaxPooling2D(pool_size=(2, 2))(conv2)

    conv3 = Convolution2D(width*8, 3, 3, activation='relu', border_mode='same')(pool2)
    #conv3 = BatchNormalization(axis = 1)(conv3)
    conv3 = Dropout(0.2)(conv3)
    conv3 = Convolution2D(width*8, 3, 3, activation='relu', border_mode='same')(conv3)
    #conv3 = BatchNormalization(axis = 1)(conv3)
    pool3 = MaxPooling2D(pool_size=(2, 2))(conv3)

    conv4 = Convolution2D(width*16, 3, 3, activation='relu', border_mode='same')(pool3)
    #conv4 = BatchNormalization(axis = 1)(conv4)
    conv4 = Dropout(0.2)(conv4)
    conv4 = Convolution2D(width*16, 3, 3, activation='relu', border_mode='same')(conv4)
    #conv4 = BatchNormalization(axis = 1)(conv4)
    pool4 = MaxPooling2D(pool_size=(2, 2), border_mode='same')(conv4)

    output = Flatten(name='flatten')(pool4)
    output = Dropout(0.2)(output)
    output = Dense(256)(output)
    output = PReLU()(output)
    output = BatchNormalization()(output)
    output = Dropout(0.2)(output)
    output = Dense(128)(output)
    output = PReLU()(output)
    output = BatchNormalization()(output)
    output = Dropout(0.3)(output)
    output = Dense(2, activation='softmax', name = 'predictions')(output)

    model = Model(input=inputs, output=output)
    #model.summary()
    model.compile(loss='categorical_crossentropy', optimizer = Adam(lr = 1e-4), metrics = ['accuracy'])
    return model

def cnn_genfit(name, batch_size, cnn, samples_tr, samples_val,
               start_tr, end_tr, start_val, end_val):

    callbacks = [EarlyStopping(monitor='val_loss', patience = 3,
                                   verbose = 1),
    ModelCheckpoint('/home/w/DS_Projects/Kaggle/DS Bowl 2017/Scripts/LUNA/CNN/Checkpoints/{}.h5'.format(name),
                        monitor='val_loss',
                        verbose = 0, save_best_only = True)]

    model = cnn()
    model.fit_generator(generate_train(start_tr, end_tr, batch_size),
                        nb_epoch = 25, verbose = 1, callbacks = callbacks,
                        samples_per_epoch = samples_tr,
                        validation_data = generate_val(start_val, end_val, batch_size),
                        nb_val_samples = samples_val)
    return

#samplestr = check_shapes(0, 1398)
#samplesval = check_shapes(1398, 1594)

samplestr = 443092
samplesval = 61992

cnn_genfit('1stgentry_CNN2DClassifier', 8, cnn2d, samplestr, samplesval, 0, 1402, 1398, 1598)

max_stars_repo_name: Wrosinski/Kaggle-DSB

max_issues_repo_name: Wrosinski/Kaggle-DSB

max_forks_repo_name: Wrosinski/Kaggle-DSB


text:
import torch
import numpy as np

def normalize(input, p=2, dim=1, eps=1e-12):
    return input / input.norm(p, dim).clamp(min=eps).unsqueeze(dim).expand_as(input)

def calculate_l2_loss(x, y):
    assert x.size() == y.size()
    ret = torch.mean(torch.pow((x - y), 2))
    return ret

def get_edge_loss(vertices,faces):
    # bs*2562*3 bs*5120*3
    edges = torch.cat((faces[:,:,:2],faces[:,:,[0,2]],faces[:,:,1:]),1) # bs * (3*5120) *2
    edges_vertices = vertices.index_select(1,edges.view(-1)).\
        view(vertices.size(0)*edges.size(0),edges.size(1),edges.size(2),vertices.size(2)) #
    indices = (torch.arange(0,vertices.size(0))*(1+vertices.size(0))).type(torch.cuda.LongTensor)
    edges_vertices = edges_vertices.index_select(0,indices) # bs * (3*5120) *2 *3
    edges_len = torch.norm((edges_vertices[:,:,0]-edges_vertices[:,:,1]),2,2)
    edges_len = torch.pow(edges_len,2)
    nonzero = len(edges_len.nonzero())
    edge_loss = torch.sum(edges_len)/nonzero
    return edge_loss

def get_edge_loss_stage1(vertices,edge):
    # vertices bs*points_number*3 edge edge_number*2
    vertices_edge = vertices.index_select(1,edge.view(-1)).\
        view(vertices.size(0),edge.size(0),edge.size(1),vertices.size(2))
    vertices_edge_vector = vertices_edge[:,:,0] - vertices_edge[:,:,1]
    vertices_edge_len = torch.pow(vertices_edge_vector.norm(2,2),2)
    edge_loss = torch.mean(vertices_edge_len)
    return edge_loss

def get_normal_loss(vertices, faces, gt_normals, idx2):
    idx2 = idx2.type(torch.cuda.LongTensor).detach()
    edges = torch.cat((faces[:,:,:2],faces[:,:,[0,2]],faces[:,:,1:]),1)
    edges_vertices = vertices.index_select(1,edges.view(-1)).\
        view(vertices.size(0)*edges.size(0),edges.size(1),edges.size(2),vertices.size(2))
    indices = (torch.arange(0,vertices.size(0))*(1+vertices.size(0))).type(torch.cuda.LongTensor)
    edges_vertices = edges_vertices.index_select(0,indices)
    edges_len1 = edges_vertices[:,:,0] - edges_vertices[:,:,1]
    edges_len2 = edges_vertices[:,:,1] - edges_vertices[:,:,0]
    edges_vector = torch.stack((edges_len1,edges_len2),2)
    gt_normals = gt_normals.index_select(1, idx2.contiguous().view(-1)).contiguous().view(gt_normals.size(0) * idx2.size(0),
                                                                                idx2.size(1), gt_normals.size(2))
    gt_normals = gt_normals.index_select(0, indices)
    gt_normals_edges = gt_normals.index_select(1, edges.view(-1)).view(gt_normals.size(0) * edges.size(0),
                                                                       edges.size(1), edges.size(2), gt_normals.size(2))
    gt_normals_edges = gt_normals_edges.index_select(0, indices)
    gt_normals_edges = normalize(gt_normals_edges, p=2, dim=3)
    edges_vector = normalize(edges_vector,p=2,dim=3)
    cosine = torch.abs(torch.sum(torch.mul(edges_vector, gt_normals_edges), 3))
    nonzero = len(cosine.nonzero())
    normal_loss = torch.sum(cosine)/nonzero

    return normal_loss

def smoothness_loss_parameters(faces):
    # faces faces_number*3(array)
    print('calculating the smoothness loss parameters, gonna take a few moments')
    if hasattr(faces, 'get'):
        faces = faces.get()
    vertices = list(set([tuple(v) for v in np.sort(np.concatenate((faces[:, 0:2], faces[:, 1:3]), axis=0))])) # edges

    v0s = np.array([v[0] for v in vertices], 'int32')
    v1s = np.array([v[1] for v in vertices], 'int32')
    v2s = []
    v3s = []
    for v0, v1 in zip(v0s, v1s):
        count = 0
        for face in faces:
            if v0 in face and v1 in face:
                v = np.copy(face)
                v = v[v != v0]
                v = v[v != v1]
                if count == 0:
                    v2s.append(int(v[0]))
                    count += 1
                else:
                    v3s.append(int(v[0]))
        if len(v3s) < len(v2s):
            v3s.append(0)

    v2s = np.array(v2s, 'int32')
    v3s = np.array(v3s, 'int32')
    print('calculated')
    return v0s, v1s, v2s, v3s

def get_smoothness_loss_stage1(vertices, parameters, eps=1e-6):
    # make v0s, v1s, v2s, v3s
    # vertices (bs*num_points*3)
    v0s, v1s, v2s, v3s = parameters
    batch_size = vertices.size(0)

    v0s = torch.from_numpy(v0s).type(torch.cuda.LongTensor)
    v1s = torch.from_numpy(v1s).type(torch.cuda.LongTensor)
    v2s = torch.from_numpy(v2s).type(torch.cuda.LongTensor)
    v3s = torch.from_numpy(v3s).type(torch.cuda.LongTensor)

    v0s = vertices.index_select(1, v0s)
    v1s = vertices.index_select(1, v1s)
    v2s = vertices.index_select(1, v2s)
    v3s = vertices.index_select(1, v3s)

    a1 = v1s - v0s
    b1 = v2s - v0s
    a1l2 = torch.sum(a1.pow(2),dim=2)
    b1l2 = torch.sum(b1.pow(2),dim=2)
    a1l1 = torch.sqrt(a1l2 + eps)
    b1l1 = torch.sqrt(b1l2 + eps)
    ab1 = torch.sum(a1*b1,dim=2)

    cos1 = ab1 / (a1l1 * b1l1 + eps)
    sin1 = torch.sqrt(1 - cos1.pow(2) + eps)
    c1 = a1 * (((ab1/(a1l2+eps)).unsqueeze(2)).expand_as(a1))

    cb1 = b1 - c1
    cb1l1 = b1l1 * sin1

    a2 = v1s - v0s
    b2 = v3s - v0s
    a2l2 = torch.sum(a2.pow(2),dim=2)
    b2l2 = torch.sum(b2.pow(2),dim=2)
    a2l1 = torch.sqrt(a2l2 + eps)
    b2l1 = torch.sqrt(b2l2 + eps)
    ab2 = torch.sum(a2*b2,dim=2)

    cos2 = ab2 / (a2l1 * b2l1 + eps)
    sin2 = torch.sqrt(1 - cos2.pow(2) + eps)
    c2 = a2 * (((ab2 / (a2l2 + eps)).unsqueeze(2)).expand_as(a2))

    cb2 = b2 - c2
    cb2l1 = b2l1 * sin2

    cos = torch.sum(cb1*cb2, dim=2) / (cb1l1 * cb2l1 + eps)
    loss = torch.sum((cos+1).pow(2)) / batch_size

    return loss

def get_smoothness_loss(vertices, parameters,faces_bn,eps=1e-6):
    # make v0s, v1s, v2s, v3s
    # vertices (bs*num_points*3)
    v0s, v1s, v2s, v3s = parameters
    batch_size = vertices.size(0)

    v0s = torch.from_numpy(v0s).type(torch.cuda.LongTensor)
    v1s = torch.from_numpy(v1s).type(torch.cuda.LongTensor)
    v2s = torch.from_numpy(v2s).type(torch.cuda.LongTensor)
    v3s = torch.from_numpy(v3s).type(torch.cuda.LongTensor)

    vs = torch.stack((v0s,v1s,v2s,v3s),1)

    faces_bn_view = faces_bn.view(faces_bn.size(0),-1)
    faces_bn_view = faces_bn_view.sort(1)[0]

    count = torch.ones(1).expand_as(faces_bn_view).type(torch.cuda.LongTensor)
    count_sum = torch.zeros(1).expand((faces_bn.size(0),vertices.shape[1])).type(torch.cuda.LongTensor)

    count_sum = count_sum.scatter_add(1, faces_bn_view, count)
    count_sum = (count_sum > 0).type(torch.cuda.ByteTensor)

    b1 = count_sum.ne(1).type(torch.cuda.LongTensor)

    i2 = vs.expand([faces_bn.size(0),vs.size(0),vs.size(1)])
    i2_unrolled = i2.view(i2.size()[0],-1)
    out_mask = torch.gather(b1,1,i2_unrolled).resize_as_(i2)
    zero_mask = out_mask.sum(2,keepdim=True).long().eq(0).long().expand_as(out_mask)
    final = i2 * zero_mask

    v0s_bn = final[:,:,0]
    v1s_bn = final[:,:,1]
    v2s_bn = final[:,:,2]
    v3s_bn = final[:,:,3]

    v0s = vertices.index_select(1, v0s_bn.view(-1)).view(vertices.size(0)*v0s_bn.size(0),v0s_bn.size(1),vertices.size(2))
    v1s = vertices.index_select(1, v1s_bn.view(-1)).view(vertices.size(0)*v1s_bn.size(0),v1s_bn.size(1),vertices.size(2))
    v2s = vertices.index_select(1, v2s_bn.view(-1)).view(vertices.size(0)*v2s_bn.size(0),v2s_bn.size(1),vertices.size(2))
    v3s = vertices.index_select(1, v3s_bn.view(-1)).view(vertices.size(0)*v3s_bn.size(0),v3s_bn.size(1),vertices.size(2))

    indices = (torch.arange(0, vertices.size(0)) * (1 + vertices.size(0))).type(torch.cuda.LongTensor)

    v0s = v0s.index_select(0,indices)
    v1s = v1s.index_select(0,indices)
    v2s = v2s.index_select(0,indices)
    v3s = v3s.index_select(0,indices)

    a1 = v1s - v0s
    b1 = v2s - v0s
    a1l2 = torch.sum(a1.pow(2),dim=2)
    b1l2 = torch.sum(b1.pow(2),dim=2)
    a1l1 = torch.sqrt(a1l2 + eps)
    b1l1 = torch.sqrt(b1l2 + eps)
    ab1 = torch.sum(a1*b1,dim=2)

    cos1 = ab1 / (a1l1 * b1l1 + eps)
    sin1 = torch.sqrt(1 - cos1.pow(2) + eps)
    c1 = a1 * (((ab1/(a1l2+eps)).unsqueeze(2)).expand_as(a1))
    cb1 = b1 - c1
    cb1l1 = b1l1 * sin1

    a2 = v1s - v0s
    b2 = v3s - v0s
    a2l2 = torch.sum(a2.pow(2),dim=2)
    b2l2 = torch.sum(b2.pow(2),dim=2)
    a2l1 = torch.sqrt(a2l2 + eps)
    b2l1 = torch.sqrt(b2l2 + eps)
    ab2 = torch.sum(a2*b2,dim=2)
    cos2 = ab2 / (a2l1 * b2l1 + eps)
    sin2 = torch.sqrt(1 - cos2.pow(2) + eps)
    c2 = a2 * (((ab2 / (a2l2 + eps)).unsqueeze(2)).expand_as(a2))

    cb2 = b2 - c2
    cb2l1 = b2l1 * sin2

    cos = torch.sum(cb1*cb2, dim=2) / (cb1l1 * cb2l1 + eps)
    loss = torch.sum((cos+1).pow(2)) / batch_size
    return loss

ext: py

max_stars_repo_name: mulinlanmao/TMNet

max_issues_repo_name: mulinlanmao/TMNet

max_forks_repo_name: mulinlanmao/TMNet

reason: import numpy


text:
#!/usr/bin/env python

from optparse import OptionParser
import sys, copy, os, pdb, random, shutil, subprocess, time
import h5py
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import numpy as np
import pandas as pd
from scipy.stats import spearmanr
import seaborn as sns
from sklearn import preprocessing
import dna_io
from itertools import product, combinations
import scipy.stats as stats

os.environ['BASSETDIR']=""
################################################################################
# basset_motifs.py
#
# Collect statistics and make plots to explore the first convolution layer
# of the given model using the given sequences.
################################################################################

#sns.set()

################################################################################
# main
################################################################################
def main():
    usage = 'usage: %prog [options] <model_file> <input_file>'
    parser = OptionParser(usage)
    parser.add_option('-a', dest='act_t', default=0.0, type='float', help='Activation threshold (as proportion of max) to consider for PWM [Default: %default]')
    parser.add_option('-c', dest='conv_hdf5_file', default=None, help='Pre-computed model output as HDF5.')
    parser.add_option('-o', dest='out_dir', default='.')
    parser.add_option('-p', dest='plot_logo', default=False, action='store_true', help='Plot LOGO (in probability) for weights of each filter [Default: %default]')
    parser.add_option('-t', dest='test_set', default=False, action='store_true', help='Generate positional convolution heatmap for each sequence in the test set [Default: %default]')
    (options,args) = parser.parse_args()

    if (options.conv_hdf5_file is None) and (len(args) != 2):
        parser.error('Must provide Basset model file and test data in HDF5 format.')
    elif options.conv_hdf5_file is None:
        model_file = args[0]
        input_file = args[1]

        if not os.path.isdir(options.out_dir):
            os.mkdir(options.out_dir)

        #################################################################
        # load HDF5 file
        #################################################################
        try:
            # input_file is HDF5
            input_hdf5 = input_file
            test_hdf5_in = h5py.File(input_hdf5, 'r')
            seqs_1hot = np.array(test_hdf5_in['test_in'])
            seq_class = np.array(test_hdf5_in['test_out'])
            test_hdf5_in.close()

        except:
            # input_file is not HDF5
            print "Input file has to be a HDF5 file."
            sys.exit()

        #################################################################
        # Torch Get Model Layers
        #################################################################
        options.conv_hdf5_file = '%s/convolution_out.h5' % options.out_dir
        #torch_cmd = '%s/src/basset_motifs_predict.lua %s %s %s' % (os.environ['BASSETDIR'],model_file, input_hdf5, options.conv_hdf5_file)
        torch_cmd = '%s/src/basset_motifs_positional_infl.lua -batch_size %d %s %s %s' % (os.environ['BASSETDIR'], seq_class.shape[0], model_file, input_hdf5, options.conv_hdf5_file)
        print torch_cmd
        subprocess.call(torch_cmd, shell=True)

    # load convolution output
    conv_hdf5_in = h5py.File(options.conv_hdf5_file, 'r')
    filter_weights = np.array(conv_hdf5_in['filter_weights'])    #shape (50, 4, 5)
    relu_outs = np.array(conv_hdf5_in['relu_outs'])              #shape (131, 50, 400) this is ReLU outputs
    norm_beta = np.array(conv_hdf5_in['norm_beta'])       #shape (50, )
    norm_gamma = np.array(conv_hdf5_in['norm_gamma'])     #shape (50, )
    norm_mean = np.array(conv_hdf5_in['norm_mean'])       #shape (50, )
    norm_var = np.array(conv_hdf5_in['norm_var'])         #shape (50, )
    pos_infl_entropy = np.array(conv_hdf5_in['filter_infl_entropy'])   #shape (50, 400)
    pos_infl_targets = np.array(conv_hdf5_in['filter_infl_targets'])   #shape (50, 3, 400)
    original_entropy = np.array(conv_hdf5_in['original_entropy'])[0]   #single value
    conv_hdf5_in.close()

    # store useful variables
    num_filters = filter_weights.shape[0]
    filter_size = filter_weights.shape[2]

    #################################################################
    # individual filter plots
    #################################################################
    # also save information contents
    filters_ic = []
    #meme_out = meme_intro('%s/filters_meme.txt'%options.out_dir, seqs)

    # generate all possible kmers
    kmers=np.array(list(product([[0,0,0,1], [0,0,1,0], [0,1,0,0],[1,0,0,0]], repeat=filter_size)))
    kmers=np.swapaxes(kmers,1,2) #shape (50,4,5)

    #For the paper
    #plot_pos_infl_entropy_heatmap(pos_infl_entropy, filter_size, [1,9,10,18,21,22,25,27,29,37,47,49], "RdBu_r", '%s/positional_infl_entropy_heatmap.pdf' % options.out_dir)

    #For the patent
    motif_index1 = [18, 25, 49, 9, 10, 21, 29, 47, 1, 22, 27, 37]
    plot_pos_infl_entropy_heatmap(pos_infl_entropy, filter_size, motif_index1, "RdBu_r", '%s/positional_infl_entropy_heatmap_inclusion_all.pdf' % options.out_dir)
    #motif_index2 = [41,5,20,2,3,40,43,42,44,38,47,29,12,31,21,9,10]
    #plot_pos_infl_entropy_heatmap(pos_infl_entropy, filter_size, motif_index2, "RdBu_r", '%s/positional_infl_entropy_heatmap_exclusion_all.pdf' % options.out_dir)
    #motif_index3 = [37,16,1,22,48,17,27,33,6,46]
    #plot_pos_infl_entropy_heatmap(pos_infl_entropy, filter_size, motif_index3, "RdBu_r", '%s/positional_infl_entropy_heatmap_unchanged_all.pdf' % options.out_dir)

    #plot_pos_pval_heatmap(relu_outs, filter_size, seq_class, [18,25,49], ["Inclusion","Stable"], "Reds", '%s/positional_IvU_heatmap.pdf' % options.out_dir)
    #plot_pos_pval_heatmap(relu_outs, filter_size, seq_class, [18,25,49], ["Exclusion","Stable"], "Blues", '%s/positional_EvU_heatmap.pdf' % options.out_dir)

    #for f in range(num_filters) :
    #    print 'Filter %d' % f

        # plot LOGO (in probability) for weights of this filter
        #plot_weight_logo(filter_weights[f,:,:], kmers, norm_beta[f], norm_gamma[f], norm_mean[f], norm_var[f], '%s/filter%d'%(options.out_dir,f), maxpct_t=options.act_t)

        # plot positional heatmap
        #set seaborn global params

    #    plot_pos_heatmap_single(relu_outs[:,f,:], filter_size, seq_class, '%s/filter%d_positional_heatmap.pdf' % (options.out_dir,f))

        # write possum motif file
        #filter_possum(filter_weights[f,:,:], 'filter%d'%f, '%s/filter%d_possum.txt'%(options.out_dir,f), options.trim_filters)

        # make a PWM for the filter
        #filter_pwm, nsites = make_filter_pwm('%s/filter%d_logo.fa'%(options.out_dir,f))

        #if nsites < 10:
            # no information
        #    filters_ic.append(0)
        #else:
            # compute and save information content
        #    filters_ic.append(info_content(filter_pwm))

            # add to the meme motif file
        #    meme_add(meme_out, f, filter_pwm, nsites, options.trim_filters)

    #meme_out.close()

################################################################################
# plot_pos_pval_heatmap
#
# Plot positional heatmap for a single filter based on ReLU output, grouped by treatment response
#
# Input
#  ReLU_out_matrix: np.array of the ReLU output
#  out_pdf:
################################################################################
def plot_pos_infl_entropy_heatmap(mat, filter_size, filter_index, color, out_prefix, seg=4):
    #figure out gap size
    gap = filter_size - 1

    #figure out padding size
    pad_side = int(gap // 2)

    #figure out interval between two sub-regions
    subinterval = int(mat.shape[-1] / seg)

    #figure out sub-region length
    sublen = subinterval - gap

    #figure out range of each sub-region
    subs=[]
    for iseg in range(seg):
        subs.append(range(pad_side+subinterval*iseg,pad_side+subinterval*iseg+sublen))

    p_tensor = mat[filter_index,:]

    #plot heatmap
    yticks=["Motif%s" % x for x in filter_index]
    fig = plt.figure(figsize=(8.27, 11.69/3))
    gs =  gridspec.GridSpec(3, seg)
    gs.update(left=0.02, right=0.98, top=0.98, bottom=0.02, wspace=0.02, hspace=0.02)
    for sub, g in zip(subs,range(seg)):
        if g==0 or g==2:
            xticks=list(range(-24,-24+sublen))
        else:
            xticks=list(range(-75,-75+sublen))
        if g>0:
            yticks=False
        ax1=plt.subplot(gs[0,g])
        sns.heatmap(p_tensor[:,sub], linecolor="black", cmap=color, linewidths=0.1, xticklabels=xticks, yticklabels=yticks, cbar=True, cbar_kws={"orientation": "horizontal"}, ax=ax1,vmin=-0.01, vmax=0.01)

    plt.savefig(out_prefix,dpi=300)
    plt.close(fig)

################################################################################
# plot_pos_pval_heatmap
#
# Plot positional heatmap for a single filter based on ReLU output, grouped by treatment response
#
# Input
#  ReLU_out_matrix: np.array of the ReLU output
#  out_pdf:
################################################################################
def plot_pos_pval_heatmap(relu_outs_3D, filter_size, class_vec, filter_index, comparison, color, out_prefix, seg=4):
    #check if the comparison is and only is between two levels
    if len(comparison)!=2:
        print "'comparison' argument has to be a list of length two, where the second one defines the baseline level."
        sys.exit()

    #figure out gap size
    gap = filter_size - 1

    #figure out padding size
    pad_side = int(gap // 2)

    #figure out interval between two sub-regions
    subinterval = int(relu_outs_3D.shape[-1] / seg)

    #figure out sub-region length
    sublen = subinterval - gap

    #figure out range of each sub-region
    subs=[]
    for iseg in range(seg):
        subs.append(range(pad_side+subinterval*iseg,pad_side+subinterval*iseg+sublen))

    #define the total number of test, which is required to correct p values
    cate=list(set(class_vec))
    total_cmps = list(combinations(range(len(cate)),2))
    total_test = len(total_cmps)

    #figure out which sequences are in the corresponding groups
    seq_index1=[i for i in range(len(class_vec)) if class_vec[i] == comparison[0]]
    seq_index2=[i for i in range(len(class_vec)) if class_vec[i] == comparison[1]]

    #ReLU 3D tensor for sequences in two classes
    relu_outs1=relu_outs_3D[seq_index1,:,:][:,filter_index,:]
    relu_outs2=relu_outs_3D[seq_index2,:,:][:,filter_index,:]

    #Initiate p value tensor
    p_tensor=np.zeros((len(filter_index),relu_outs_3D.shape[-1]),dtype=np.float32)

    #1-way ANOVA for each position
    print comparison[0],"vs",comparison[1]
    for mi in range(len(filter_index)):
        for ri in range(p_tensor.shape[-1]):
            t,p=stats.ttest_ind(relu_outs1[:,mi,ri],relu_outs2[:,mi,ri])
            q = p*total_test
            if np.isnan(p):
                q=1
            elif q >= 0.05:
                q=1
            p_tensor[mi,ri]= np.log10(q)*(-1)

    #plot heatmap
    yticks=["Motif%s" % x for x in filter_index]

    fig = plt.figure(figsize=(8.27, 11.69/12))
    gs =  gridspec.GridSpec(1, seg)
    gs.update(left=0.02, right=0.98, top=0.98, bottom=0.02, wspace=0.02, hspace=0.02)
    for sub, g in zip(subs,range(seg)):
        if g==0 or g==2:
            xticks=list(range(-24,-24+sublen))
        else:
            xticks=list(range(-75,-75+sublen))
        if g>1:
            yticks=False
        ax1=plt.subplot(gs[0:2,g])
        sns.heatmap(p_tensor[:,sub], linecolor="black", cmap=color, linewidths=0.1, xticklabels=xticks, yticklabels=yticks, cbar=True, cbar_kws={"orientation": "horizontal"}, ax=ax1)

    plt.savefig(out_prefix,dpi=300)
    plt.close(fig)

################################################################################
# plot_weight_logo
#
# Plot a weblogo of the filter's weight of the model
#
# Input
#  param_matrix: np.array of the filter's parameter matrix
#  out_pdf:
################################################################################
def plot_weight_logo(filter_weights, seq_tensor, norm_beta, norm_gamma, norm_mean, norm_var, out_prefix, maxpct_t=0.0):
    #Get convolution
    kmer_conv=np.tensordot(filter_weights,seq_tensor,axes=((0,1),(1,2))) #result in a 1D array

    #normalization using trained parameters
    kmer_conv_norm = (kmer_conv - norm_mean) * norm_beta / np.sqrt(norm_var + 0.00001) + norm_gamma

    #ReLU
    kmer_conv_relu = [0 if val <0 else val for val in kmer_conv_norm]

    #calculate a ReLU cutoff to plot weblogo
    if maxpct_t == 0:
        relu_act = 0
    else:
        all_outs = np.ravel(kmer_conv_relu)
        all_outs_mean = all_outs.mean()
        all_outs_norm = all_outs - all_outs_mean
        relu_act = maxpct_t * all_outs_norm.max() + all_outs_mean

    #Find which kmer pass the cutoff
    kmer_ok_index=[i for i in range(len(kmer_conv_relu)) if kmer_conv_relu[i]>relu_act]
    if len(kmer_ok_index)>0:
        kmer_ok=np.array([seq_tensor[i,:,:] for i in kmer_ok_index]) #shape (50,4,5)
        kmer_ok_4D=np.expand_dims(kmer_ok, axis=2) #shape (50,4,1,5)
        kmer_seq = dna_io.vecs2dna(kmer_ok_4D)
        kmer_fasta=[">"+str(i)+"\n"+seq+"\n" for i,seq in zip(range(len(kmer_seq)),kmer_seq)]
        kmer_fasta_out=open("%s_activated_kmers.fa" % out_prefix,"w")
        kmer_fasta_out.writelines(kmer_fasta)
        kmer_fasta_out.close()
        if relu_act > 0:
            subprocess.call("weblogo -X NO -Y NO -F pdf --resolution 300 --errorbars NO --fineprint '' -C '#CB2026' A A -C '#34459C' C C -C '#FBB116' G G -C '#0C8040' T T < %s_activated_kmers.fa > %s_weight_LOGO.pdf"%(out_prefix, out_prefix), shell=True)
        else:
            subprocess.call("weblogo -U probability -F pdf --resolution 300 --errorbars NO --fineprint '' -C '#CB2026' A A -C '#34459C' C C -C '#FBB116' G G -C '#0C8040' T T < %s_activated_kmers.fa > %s_weight_LOGO.pdf"%(out_prefix, out_prefix), shell=True)

################################################################################
# plot_pos_heatmap_single
#
# Plot positional heatmap for a single filter based on ReLU output, grouped by treatment response
#
# Input
#  ReLU_out_matrix: np.array of the ReLU output
#  out_pdf:
################################################################################
def plot_pos_heatmap_single(relu_outs, filter_size, class_vec, out_prefix, byclass=True, seg=4):
    #figure out gap size
    gap = filter_size - 1

    #figure out padding size
    pad_side = gap // 2

    #figure out interval between two sub-regions
    subinterval = relu_outs.shape[-1] / seg

    #figure out sub-region length
    sublen = subinterval - gap

    #figure out range of each sub-region
    subs=[]
    for iseg in range(seg):
        subs.append(range(pad_side+subinterval*iseg,pad_side+subinterval*iseg+sublen))

    #define class
    cate=list(set(class_vec))

    #Initiate plot
    fig = plt.figure(figsize=(8.27, 11.69))
    gs =  gridspec.GridSpec(6, seg)
    gs.update(left=0.02, right=0.98, top=0.58, bottom=0.02, wspace=0.1, hspace=0.2)

    conv_vec_byclass=[]
    #group relu vectors into classes
    for aclass, i in zip(cate,range(len(cate))):
        conv_vec=[]
        for cand, j in zip(class_vec,range(relu_outs.shape[0])):
            if cand==aclass:
                conv_vec.append(relu_outs[j,:])
        conv_vec=np.array(conv_vec)

        #copy into a final list of each class
        conv_vec_byclass.append(conv_vec)
        #conv_vec_colmean=np.mean(conv_vec,axis=0)
        #conv_vec_colmean=np.expand_dims(conv_vec_colmean,0)

    #1-way ANOVA
    cmps = list(combinations(range(len(cate)),2))
    #print cmps
    p_tensor=np.zeros((len(cmps),relu_outs.shape[1]),dtype=np.float32)
    for ci in range(len(cmps)):
        print cate[cmps[ci][0]],"vs",cate[cmps[ci][1]]
        p_vec=np.zeros((1,relu_outs.shape[1]),dtype=np.float32)
        for ri in range(p_tensor.shape[1]):
            t,p=stats.ttest_ind(conv_vec_byclass[cmps[ci][0]][:,ri],conv_vec_byclass[cmps[ci][1]][:,ri])
            q = p*len(cmps)
            if np.isnan(p):
                q=1
            elif q >= 0.05:
                q=1
            p_vec[0,ri]= np.log10(q)*(-1)
        p_tensor[ci,:]=p_vec

    #plot
    plot_cbar=False
    for sub, g in zip(subs,range(seg)):
        #lineplot
        ax1=plt.subplot(gs[0:2,g])
        for data in conv_vec_byclass:
            if g==0 or g==2:
                junc=list(range(-24,-24+sublen))*data.shape[0]
            else:
                junc=list(range(-75,-75+sublen))*data.shape[0]
            sns.lineplot(x=junc, y=np.ravel(data[:,sub]), ax=ax1)
            if g>0:
                ax1.tick_params(axis="y",labelleft=False)
        sns.despine()

        #heatmap
        ax2=plt.subplot(gs[3,g])
        if g==3:
            plot_cbar=False
        sns.heatmap(p_tensor[:,sub], linecolor="white", linewidths=.2, cbar=plot_cbar, xticklabels=False, yticklabels=False, ax=ax2)

    plt.savefig(out_prefix,dpi=300)
    plt.close(fig)

################################################################################
# __main__
################################################################################
if __name__ == '__main__':
    main()

reason: import numpy,import scipy,from scipy


text:
#Copyright (c) 2016 PaddlePaddle Authors. All Rights Reserved
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import os
import sys
import time
import argparse
import logging
import paddle.fluid as fluid
import paddle
import utils
import numpy as np
from nets import SequenceSemanticRetrieval

logging.basicConfig(format="%(asctime)s - %(levelname)s - %(message)s")
logger = logging.getLogger("fluid")
logger.setLevel(logging.INFO)

def parse_args():
    parser = argparse.ArgumentParser("sequence semantic retrieval")
    parser.add_argument(
        "--train_dir", type=str, default='train_data', help="Training file")
    parser.add_argument(
        "--base_lr", type=float, default=0.01, help="learning rate")
    parser.add_argument(
        '--vocab_path', type=str, default='vocab.txt', help='vocab file')
    parser.add_argument(
        "--epochs", type=int, default=10, help="Number of epochs")
    parser.add_argument(
        '--parallel', type=int, default=0, help='whether parallel')
    parser.add_argument(
        '--use_cuda', type=int, default=0, help='whether use gpu')
    parser.add_argument(
        '--print_batch', type=int, default=10, help='num of print batch')
    parser.add_argument(
        '--model_dir', type=str, default='model_output', help='model dir')
    parser.add_argument(
        "--hidden_size", type=int, default=128, help="hidden size")
    parser.add_argument(
        "--batch_size", type=int, default=50, help="number of batch")
    parser.add_argument(
        "--embedding_dim", type=int, default=128, help="embedding dim")
    parser.add_argument(
        '--num_devices', type=int, default=1, help='Number of GPU devices')
    parser.add_argument(
        '--step_num', type=int, default=1000, help='Number of steps')
    parser.add_argument(
        '--enable_ce',
        action='store_true',
        help='If set, run the task with continuous evaluation logs.')
    parser.add_argument(
        '--role', type=str, default='pserver', help='trainer or pserver')
    parser.add_argument(
        '--endpoints',
        type=str,
        default='127.0.0.1:6000',
        help='The pserver endpoints, like: 127.0.0.1:6000, 127.0.0.1:6001')
    parser.add_argument(
        '--current_endpoint',
        type=str,
        default='127.0.0.1:6000',
        help='The current_endpoint')
    parser.add_argument(
        '--trainer_id',
        type=int,
        default=0,
        help='trainer id ,only trainer_id=0 save model')
    parser.add_argument(
        '--trainers',
        type=int,
        default=1,
        help='The num of trianers, (default: 1)')
    return parser.parse_args()

def get_cards(args):
    return args.num_devices

def train_loop(main_program, avg_cost, acc, train_input_data, place, args,
               train_reader):
    data_list = [var.name for var in train_input_data]
    feeder = fluid.DataFeeder(feed_list=data_list, place=place)
    exe = fluid.Executor(place)
    exe.run(fluid.default_startup_program())
    train_exe = exe

    total_time = 0.0
    ce_info = []
    for pass_id in range(args.epochs):
        epoch_idx = pass_id + 1
        print("epoch_%d start" % epoch_idx)
        t0 = time.time()
        i = 0
        for batch_id, data in enumerate(train_reader()):
            i += 1
            loss_val, correct_val = train_exe.run(
                feed=feeder.feed(data), fetch_list=[avg_cost.name, acc.name])
            ce_info.append(float(np.mean(correct_val)) / args.batch_size)
            if i % args.print_batch == 0:
                logger.info(
                    "Train --> pass: {} batch_id: {} avg_cost: {}, acc: {}".
                    format(pass_id, batch_id,
                           np.mean(loss_val),
                           float(np.mean(correct_val)) / args.batch_size))
            if args.enable_ce and i > args.step_num:
                break
        t1 = time.time()
        total_time += t1 - t0
        print("epoch:%d num_steps:%d time_cost(s):%f" %
              (epoch_idx, i, total_time / epoch_idx))
        save_dir = "%s/epoch_%d" % (args.model_dir, epoch_idx)
        fluid.save(fluid.default_main_program(), save_dir)
        print("model saved in %s" % save_dir)

    # only for ce
    if args.enable_ce:
        ce_acc = 0
        try:
            ce_acc = ce_info[-2]
        except:
            print("ce info error")
        epoch_idx = args.epochs
        device = get_device(args)
        if args.use_cuda:
            gpu_num = device[1]
            print("kpis\teach_pass_duration_gpu%s\t%s" %
                  (gpu_num, total_time / epoch_idx))
            print("kpis\ttrain_acc_gpu%s\t%s" % (gpu_num, ce_acc))
        else:
            cpu_num = device[1]
            threads_num = device[2]
            print("kpis\teach_pass_duration_cpu%s_thread%s\t%s" %
                  (cpu_num, threads_num, total_time / epoch_idx))
            print("kpis\ttrain_acc_cpu%s_thread%s\t%s" %
                  (cpu_num, threads_num, ce_acc))

def train(args):
    if args.enable_ce:
        SEED = 102
        fluid.default_startup_program().random_seed = SEED
        fluid.default_main_program().random_seed = SEED
    use_cuda = True if args.use_cuda else False
    parallel = True if args.parallel else False
    print("use_cuda:", use_cuda, "parallel:", parallel)
    train_reader, vocab_size = utils.construct_train_data(
        args.train_dir, args.vocab_path, args.batch_size * get_cards(args))
    place = fluid.CUDAPlace(0) if use_cuda else fluid.CPUPlace()
    ssr = SequenceSemanticRetrieval(vocab_size, args.embedding_dim,
                                    args.hidden_size)
    # Train program
    train_input_data, cos_pos, avg_cost, acc = ssr.train()

    # Optimization to minimize lost
    optimizer = fluid.optimizer.Adagrad(learning_rate=args.base_lr)
    optimizer.minimize(avg_cost)

    print("run distribute training")
    t = fluid.DistributeTranspiler()
    t.transpile(
        args.trainer_id, pservers=args.endpoints, trainers=args.trainers)
    if args.role == "pserver":
        print("run psever")
        pserver_prog = t.get_pserver_program(args.current_endpoint)
        pserver_startup = t.get_startup_program(args.current_endpoint,
                                                pserver_prog)
        exe = fluid.Executor(fluid.CPUPlace())
        exe.run(pserver_startup)
        exe.run(pserver_prog)
    elif args.role == "trainer":
        print("run trainer")
        train_loop(t.get_trainer_program(), avg_cost, acc, train_input_data,
                   place, args, train_reader)

def get_device(args):
    if args.use_cuda:
        gpus = os.environ.get("CUDA_VISIBLE_DEVICES", 1)
        gpu_num = len(gpus.split(','))
        return "gpu", gpu_num
    else:
        threads_num = os.environ.get('NUM_THREADS', 1)
        cpu_num = os.environ.get('CPU_NUM', 1)
        return "cpu", int(cpu_num), int(threads_num)

def main():
    args = parse_args()
    train(args)

if __name__ == "__main__":
    main()

max_stars_repo_name: suytingwan/models

max_issues_repo_name: suytingwan/models

max_forks_repo_name: suytingwan/models


text:
import math
import numpy as np
np.set_printoptions(suppress=True)

print(np.random.rand(100,100))


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Oct  8 12:19:47 2021

@author: Hendrik Timm
"""

import numpy as np
from code.feature_extraction.feature_extractor import FeatureExtractor

# class for extracting the number of hashtags in a tweet as a feature
class Hashtags(FeatureExtractor):

    # constructor
    def __init__(self, input_column):
        super().__init__([input_column], "{0}_#hashtags".format(input_column))

    def _get_values(self, inputs):

        num_hashtags = np.zeros((len(inputs[0]),1))
        for i in range(len(inputs[0])):

            num_hashtags[i][0] = inputs[0][i].count(",")

            if(len(inputs[0][i])>2):
                num_hashtags[i][0] +=1

        return num_hashtags

max_stars_repo_name: Jackylinelein/MLinPractice

max_issues_repo_name: Jackylinelein/MLinPractice

max_forks_repo_name: Jackylinelein/MLinPractice


text:
"""
Computation of the measure 'AUC Mu'. This measure requires installation of the
numpy, Apache Spark and pyspark interpreter.
This code corresponds to the paper: Kleiman, R., Page, D. ``AUC Mu: A
Performance Metric for Multi-Class Machine Learning Models``, Proceedings of the
2019 International Conference on Machine Learning (ICML).
"""

__author__ = "PoWeiHuang"
__copyright__ = "Copyright 2019"
__credits__ = ["PoWeiHuang"]
__license__ = "MIT"
__version__ = "1.0"
__maintainer__ = "PoWeiHuang"
__email__ = "poweihuang@cathayholdings.com.tw"
__status__ = "Production"

import numpy as np
from pyspark.mllib.linalg import Vector, VectorUDT, DenseVector
from pyspark.ml.evaluation import BinaryClassificationEvaluator
from pyspark.sql.types import IntegerType
from pyspark.ml.linalg import Vectors
from pyspark.sql.functions import col, udf, DataFrame, avg

#----------------------------------------------------------------------
def pyspark_auc_mu(data, y_true, y_score, A=None, W=None):
    """
    Compute the multi-class measure AUC Mu from prediction scores and labels from Spark dataframe.

    Parameters
    ----------
    data : Spark dataframe
        The prediction output table from pyspark.ml

    y_true : double type column in data, shape = [n_samples]
        The true class labels column from spark dataframe in the range [0, n_samples-1]

    y_score : vector type column in data, shape = [n_samples, n_classes]
        Target scores, where each element in a vector is a categorical distribution over the
        n_classes.

    A : array, shape = [n_classes, n_classes], optional
        The partition (or misclassification cost) matrix. If ``None`` A is the
        argmax partition matrix. Entry A_{i,j} is the cost of classifying an
        instance as class i when the true class is j. It is expected that
        diagonal entries in A are zero and off-diagonal entries are positive.

    W : array, shape = [n_classes, n_classes], optional
        The weight matrix for incorporating class skew into AUC Mu. If ``None``,
        the standard AUC Mu is calculated. If W is specified, it is expected to
        be a lower triangular matrix where entrix W_{i,j} is a positive float
        from 0 to 1 for the partial score between classes i and j. Entries not
        in the lower triangular portion of W must be 0 and the sum of all
        entries in W must be 1.

    Returns
    -------
    auc_mu : float

    References
    ----------
    .. [1] Kleiman, R., Page, D. ``AUC Mu: A Performance Metric for Multi-Class
           Machine Learning Models``, Proceedings of the 2019 International
           Conference on Machine Learning (ICML).
       [2] https://github.com/kleimanr/auc_mu
    """

    n_classes = data.select(y_true).distinct().count()
    n_samples = data.select(y_score).count()
    # Validate input arguments
    if not isinstance(data, DataFrame):
        raise TypeError("Expected data to be DataFrame, got: %s"
                        % type(data))
    if not data.select(y_true).dtypes[0][1] == 'double':
        raise TypeError("Expected column y_true to be double, got: %s"
                        % data.select(y_true).dtypes[0][1])
    if not data.select(y_score).dtypes[0][1] == 'vector':
        raise TypeError("Expected column y_score to be vector, got: %s"
                        % data.select(y_true).dtypes[0][1])
    if not data.select(y_true).count() == n_samples:
        raise ValueError("Expected y_true to be shape %s, got: %s"
                        %(str(data.select(y_score).count()), str(data.select(y_true).count())))

    slen = udf(lambda s: len(s), IntegerType())
    if not data.select(slen(col(y_score))).groupBy().avg().collect()[0][0] == n_classes:
        raise ValueError("Expected y_true values in range 0..%i, got: %s"
                        %(n_classes-1, str(data.select(slen(data.y_score)).groupBy().avg().collect()[0][0])))
    if A is None:
        A = np.ones((n_classes, n_classes)) - np.eye(n_classes)
    if not isinstance(A, np.ndarray):
        raise TypeError("Expected A to be np.ndarray, got: %s"
                        % type(A))
    if not A.ndim == 2:
        raise ValueError("Expected A to be 2 dimensional, got: %s"
                         % A.ndim)
    if not A.shape == (n_classes, n_classes):
        raise ValueError("Expected A to be shape (%i, %i), got: %s"
                         %(n_classes, n_classes, str(A.shape)))
    if not np.all(A.diagonal() == np.zeros(n_classes)):
        raise ValueError("Expected A to be zero on the diagonals")
    if not np.all(A >= 0):
        raise ValueError("Expected A to be non-negative")

    if W is None:
        W = np.tri(n_classes, k=-1)
        W /= W.sum()
    if not isinstance(W, np.ndarray):
        raise TypeError("Expected W to be np.ndarray, got: %s"
                        % type(W))
    if not W.ndim == 2:
        raise ValueError("Expected W to be 2 dimensional, got: %s"
                         % W.ndim)
    if not W.shape == (n_classes, n_classes):
        raise ValueError("Expected W to be shape (%i, %i), got: %s"
                         %(n_classes, n_classes, str(W.shape)))

    auc_total = 0.0
    for class_i in xrange(n_classes):
        preds_i = data.select(y_score).where(col(y_true) == class_i)
        n_i = preds_i.count()

        for class_j in xrange(class_i):
            preds_j = data.select(y_score).where(col(y_true) == class_j)
            temp_preds = preds_i.union(preds_j)

            n_j = preds_j.count()
            n = n_i+n_j
            #temp_preds: concat prob vectors which class = i and j
            temp_preds = DenseVector(temp_preds.select("probability").rdd.map(lambda x: x[0]).collect())

            #temp_labels: convert the two selected classes to a binary class vector
            temp_labels = np.zeros((n), dtype=int)
            temp_labels[n_i:n] = 1

            # v: differencing by vector_{i,.} and vector_{j,.} in partition matrix
            v = A[class_i, :] - A[class_j, :]
            score = temp_preds.dot(v)
            df = np.column_stack([score, temp_labels])
            concat_df = map(lambda x: (float(x[0]), float(x[1:])), df)
            auc_mu_df = spark.createDataFrame(concat_df,schema=["score", "temp_label"])

            evaluator = BinaryClassificationEvaluator(labelCol="temp_label", rawPredictionCol="score", metricName='areaUnderROC')
            score_i_j = evaluator.evaluate(auc_mu_df)
            auc_total += W[class_i, class_j]*score_i_j

    return auc_total

max_stars_repo_name: poweihuang/spark-aucmu

max_issues_repo_name: poweihuang/spark-aucmu

max_forks_repo_name: poweihuang/spark-aucmu


text:
import scipy
import scipy.integrate
import scipy.interpolate
from scipy.interpolate import make_interp_spline as mkspl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.axes import Axes
from matplotlib.figure import Figure
import numpy as np
from numpy import array, linalg, matrix
from scipy.special import comb as nOk
import deepracing_models.math_utils as math_utils
import torch
def lsqfit(points,M):
    M_ = linalg.pinv(M)
    return M_ * points
pi = np.pi
num_points = 1000
tmin = -1.0
tmax = 1.0
d = 4
kbezier = d - 0
numcurves = d
Mtk = lambda i, n, t: t**(i)*(1-t)**(n-i)*nOk(n,i)
bezierFactors = lambda ts, n: matrix([[Mtk(k_,n,t) for k_ in range(n+1)] for t in ts])
P = np.zeros((numcurves,num_points,2))
Pprime = np.zeros((numcurves,num_points,2))
Pprimeprime = np.zeros((numcurves,num_points,2))
tnp = np.zeros((numcurves,num_points))
t = np.linspace(tmin, tmax, num_points )
tfunc = t
for order in range(1,numcurves+1):
    i = order-1
   #P[i] = np.zeros((num_points,2))
   # tfunc = t-t[int(len(t)/2)]
    tnp[i] = tfunc.copy()
    px = 1.0*np.random.randn(order)
    py = 1.0*np.random.randn(order)
    P[i,:,0] = np.polyval(px,tnp[i])
    P[i,:,1] = np.polyval(py,tnp[i])
    P[i,:,0] = P[i,:,0] - P[i,0,0]
    P[i,:,1] = P[i,:,1] - P[i,0,1]
    Pprime[i,:,0] = np.polyval(np.polyder(px, m=1),tnp[i])
    Pprime[i,:,1] = np.polyval(np.polyder(py, m=1),tnp[i])
    Pprimeprime[i,:,0] = np.polyval(np.polyder(px, m=2),tnp[i])
    Pprimeprime[i,:,1] = np.polyval(np.polyder(py, m=2),tnp[i])

Ptorch = torch.from_numpy(P.copy()).double()#.unsqueeze(0)
Pprimetorch = torch.from_numpy(Pprime.copy()).double()#.unsqueeze(0)
Pprimeprimetorch = torch.from_numpy(Pprimeprime.copy()).double()#.unsqueeze(0)
ttorch = torch.from_numpy(tnp.copy()).double()#.unsqueeze(0)
dt = ttorch[:,-1]-ttorch[:,0]
storch = (ttorch - ttorch[:,0,None])/dt[:,None]
#ttorch = ttorch.repeat(numcurves,1)
print(storch)
print(Ptorch.shape)
M , bezier_control_points = math_utils.bezierLsqfit(Ptorch, kbezier, t = storch )
print(bezier_control_points.shape)
Pbeziertorch = torch.matmul(M,bezier_control_points)

Mderiv, Pdot_s = math_utils.bezierDerivative(bezier_control_points,storch, order=1)
Pdot_t= Pdot_s/dt[:,None,None]

Mdotdot, Pdotdot_s = math_utils.bezierDerivative(bezier_control_points, storch, order=2)
Pdotdot_t= Pdotdot_s/((dt**2)[:,None,None])

fig : Figure = plt.figure()
for i in range(numcurves):
   #ax.plot(P[i,:,0],P[i,:,1],'r-')
    plt.title("A Bézier Curve and its Control Points")
    plt.plot(Pbeziertorch[i,:,0].numpy(),Pbeziertorch[i,:,1].numpy(), color='blue')#, label="A Bézier Curve")
   # ax.scatter(Pbeziertorch[i,:,0].numpy(),Pbeziertorch[i,:,1].numpy(), facecolors='none', edgecolors='b')
    skipn = 50
    #ax.plot(bezier_control_points[i,:,0].numpy(),bezier_control_points[i,:,1].numpy(),'go')
    plt.scatter(bezier_control_points[i,:,0].numpy(),bezier_control_points[i,:,1].numpy(), facecolors='none', edgecolors='g')#, label="Control Points")
    plt.quiver(Pbeziertorch[i,::skipn,0].numpy(),Pbeziertorch[i,::skipn,1].numpy(),Pdot_t[i,::skipn,0].numpy(),Pdot_t[i,::skipn,1].numpy(), angles='xy', color='black')#, label="Velocity Vectors")

    print("Mean distance: %f" %(torch.mean(torch.norm(Pbeziertorch[i]-Ptorch[i],dim=1,p=2)).item()))
    print("Mean velocity diff: %f" %(torch.mean(torch.norm(Pdot_t[i]-Pprimetorch[i],dim=1,p=2)).item()))
    print("Mean acceleration diff: %f" %(torch.mean(torch.norm(Pdotdot_t[i]-Pprimeprimetorch[i],dim=1,p=2)).item()))
plt.legend(loc="best")
plt.show()

max_stars_repo_name: linklab-uva/deepracing

max_issues_repo_name: linklab-uva/deepracing

max_forks_repo_name: linklab-uva/deepracing

reason: import numpy,from numpy,import scipy,from scipy


text:
import ampligraph.datasets
import ampligraph.latent_features
from ampligraph.evaluation import hits_at_n_score, mr_score, evaluate_performance, mrr_score

import argparse
import os
import json
import sys
import yaml
import logging

import numpy as np
from beautifultable import BeautifulTable
from tqdm import tqdm
import warnings

os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3"

warnings.simplefilter(action="ignore", category=Warning)

SUPPORT_DATASETS = ["fb15k", "fb15k-237", "wn18", "wn18rr"]
SUPPOORT_MODELS = ["complex", "transe", "distmult", "hole"]

def display_scores(scores):
    output_readme = {}
    output_rst = {}

    for obj in scores:
        output_rst[obj["dataset"]] = BeautifulTable()
        output_rst[obj["dataset"]].set_style(BeautifulTable.STYLE_RST)
        output_rst[obj["dataset"]].column_headers = \
            ["Model", "MR", "MRR", "Hits@1", \
             "Hits@3", "Hits@10", "Hyperparameters"]

    for obj in scores:
        try:
            output_rst[obj["dataset"]] \
                .append_row([obj["model"],
                             "{0:.2f}".format(obj["mr"]),
                             "{0:.2f}".format(obj["mrr"]),
                             "{0:.2f}".format(obj["H@1"]),
                             "{0:.2f}".format(obj["H@3"]),
                             "{0:.2f}".format(obj["H@10"]),
                             yaml.dump(obj["hyperparams"],
                                       default_flow_style=False)])
        except:
            output_rst[obj["dataset"]] \
                .append_row([obj["model"],
                             ".??",
                             ".??",
                             ".??",
                             ".??",
                             ".??",
                             ".??"])

    for key, value in output_rst.items():
        print(key)
        print(value)

def run_single_exp(config, dataset, model):
    hyperparams = config["hyperparams"][dataset][model]
    if hyperparams is None:
        logging.info("dataset {0}...model {1} \
                      experiment is not conducted yet..." \
                     .format(dataset, config["model_name_map"][model]))
        return {
            "hyperparams": ".??"
        }
    logging.info("dataset {0}...model {1}...\
                  best hyperparameter:...{2}" \
                 .format(dataset,
                         config["model_name_map"][model],
                         hyperparams))

    es_code = "{0}_{1}".format(dataset, model)

    load_func = getattr(ampligraph.datasets,
                        config["load_function_map"][dataset])
    X = load_func()
    # logging.debug("Loaded...{0}...".format(dataset))

    # load model
    model_class = getattr(ampligraph.latent_features,
                          config["model_name_map"][model])
    model = model_class(**hyperparams)
    # Fit the model on training and validation set
    # The entire dataset will be used to filter out false positives statements
    # created by the corruption procedure:
    filter = np.concatenate((X['train'], X['valid'], X['test']))

    if es_code in config["no_early_stopping"]:
        logging.debug("Fit without early stopping...")
        model.fit(X["train"])
    else:
        logging.debug("Fit with early stopping...")
        model.fit(X["train"], True,
                  {
                      'x_valid': X['valid'][::10],
                      'criteria': 'mrr',
                      'x_filter': filter,
                      'stop_interval': 2,
                      'burn_in': 0,
                      'check_interval': 100
                  })

    # Run the evaluation procedure on the test set. Will create filtered rankings.
    # To disable filtering: filter_triples=None
    ranks = evaluate_performance(X['test'],
                                 model,
                                 filter,
                                 verbose=False)

    # compute and print metrics:
    mr = mr_score(ranks)
    mrr = mrr_score(ranks)
    hits_1 = hits_at_n_score(ranks, n=1)
    hits_3 = hits_at_n_score(ranks, n=3)
    hits_10 = hits_at_n_score(ranks, n=10)

    return {
        "mr": mr,
        "mrr": mrr,
        "H@1": hits_1,
        "H@3": hits_3,
        "H@10": hits_10,
        "hyperparams": hyperparams
    }

def run_all(config):
    obj = []
    for dataset in tqdm(config["hyperparams"].keys(),
                        desc="evaluation done: "):
        for model in config["hyperparams"][dataset].keys():
            result = run_single_exp(config, dataset, model)
            obj.append({
                **result,
                "dataset": dataset,
                "model": config["model_name_map"][model]
            })
    return obj

def run_single_dataset(config, dataset):
    obj = []
    for model in tqdm(config["hyperparams"][dataset].keys(),
                      desc="evaluation done: "):
        result = run_single_exp(config, dataset, model)
        obj.append({
            **result,
            "dataset": dataset,
            "model": config["model_name_map"][model]
        })
    return obj

def run_single_model(config, model):
    obj = []
    for dataset in tqdm(config["hyperparams"].keys(),
                        desc="evaluation done: "):
        result = run_single_exp(config, dataset, model)
        obj.append({
            **result,
            "dataset": dataset,
            "model": config["model_name_map"][model]
        })
    return obj

def main():
    with open("config.json", "r") as fi:
        config = json.load(fi)

    # set GPU id to run
    os.environ["CUDA_VISIBLE_DEVICES"] = config["CUDA_VISIBLE_DEVICES"]
    logging.debug("Will use gpu number...{0}" \
                  .format(config["CUDA_VISIBLE_DEVICES"]))

    parser = argparse.ArgumentParser()
    parser.add_argument("-d", "--dataset",
                        type=str,
                        choices=SUPPORT_DATASETS)
    parser.add_argument("-m", "--model",
                        type=str,
                        choices=SUPPOORT_MODELS)

    args = parser.parse_args()
    logging.debug("Input dataset...{0}...input model...{1}..." \
                  .format(args.dataset, args.model))

    if args.dataset is None:
        if args.model is None:
            display_scores(run_all(config))
        else:
            if args.model.upper() not in config["model_name_map"]:
                sys.exit("Input model is not valid...")
            display_scores(run_single_model(config, args.model.upper()))
    else:
        if args.model is not None:
            if args.model.upper() not in config["model_name_map"] \
                    or args.dataset.upper() \
                            not in config["load_function_map"]:
                sys.exit("Input model or dataset is not valid...")

            result = run_single_exp(config,
                                    args.dataset.upper(),
                                    args.model.upper())
            display_scores([{
                **result,
                "dataset": args.dataset,
                "model": config["model_name_map"][args.model.upper()]
            }])
        else:
            if args.dataset.upper() not in config["load_function_map"]:
                sys.exit("Input dataset is not supported yet...")
            display_scores(run_single_dataset(config,
                                              args.dataset.upper()))

if __name__ == "__main__":
    main()

max_stars_repo_name: pyvandenbussche/AmpliGraph

max_issues_repo_name: rezacsedu/AmpliGraph

max_forks_repo_name: rezacsedu/AmpliGraph


text:
import numpy as np

from benchml.pipeline import FitTransform

class KernelBase(object):
    def __init__(self, **kwargs):
        self.name = "base"
        self.X_fit = None
        self.K_fit = None

    def evaluate(self, X1, X2, symmetric, **kwargs):
        raise NotImplementedError("<evaluate> not defined")

    def evaluateFit(self, X, y="not_used", **kwargs):
        self.X_fit = X
        self.K_fit = self.evaluate(X, X, symmetric=True, **kwargs)
        return self.K_fit

    def evaluatePredict(self, X, **kwargs):
        return self.evaluate(X, self.X_fit, symmetric=False, **kwargs)

class KernelDotDeprecated(KernelBase):
    def __init__(self, **kwargs):
        KernelBase.__init__(self, **kwargs)
        self.power = kwargs["power"]
        self.name = "dot"

    def evaluate(self, X1, X2, symmetric="not_used", **kwargs):
        k = X1.dot(X2.T)
        k = k ** self.power
        return k

class KernelDot(FitTransform):
    default_args = {"power": 1, "self_kernel": False}
    req_inputs = ("X",)
    allow_params = {"X"}
    allow_stream = {"K", "K_self"}
    stream_kernel = ("K",)
    stream_samples = ("K_self",)
    precompute = True

    def evaluate(self, x1, x2=None, diagonal_only=False):
        if x2 is None:
            x2 = x1
        if diagonal_only:
            return np.einsum("ai,ai->a", x1, x2, optimize="greedy") ** self.args["power"]
        else:
            return x1.dot(x2.T) ** self.args["power"]

    def _fit(self, inputs, stream, params):
        K = self.evaluate(inputs["X"])
        params.put("X", np.copy(inputs["X"]))
        stream.put("K", K)
        if self.args["self_kernel"]:
            stream.put("K_self", K.diagonal())

    def _map(self, inputs, stream):
        K = self.evaluate(inputs["X"], self.params().get("X"))
        stream.put("K", K)
        if self.args["self_kernel"]:
            K_self = self.evaluate(inputs["X"], inputs["X"], diagonal_only=True)
            stream.put("K_self", K_self)

class KernelGaussian(FitTransform):
    default_args = {"scale": 1, "self_kernel": False, "epsilon": 1e-10}
    req_inputs = ("X",)
    allow_params = {"X", "sigma"}
    allow_stream = {"K", "K_self"}
    stream_kernel = ("K",)
    stream_samples = ("K_self",)
    precompute = True

    def evaluate(self, x1, x2=None, sigma=None, diagonal_only=False):
        x1s = x1 / (sigma + self.args["epsilon"])
        z1 = np.sum(x1s ** 2, axis=1)
        if x2 is None:
            x2s = x1s
            z2 = z1
        else:
            x2s = x2 / sigma
            z2 = np.sum(x2s ** 2, axis=1)
        if diagonal_only:
            zz = -0.5 * (z1 + z2)
            xx = np.einsum("ai,ai->a", x1s, x2s, optimize="greedy")
        else:
            zz = -0.5 * np.add.outer(z1, z2)
            xx = x1s.dot(x2s.T)
        return np.exp(zz + xx)

    def _fit(self, inputs, stream, params):
        X = inputs["X"]
        sigma = self.args["scale"] * np.std(X, axis=0)
        K = self.evaluate(x1=inputs["X"], sigma=sigma)
        params.put("sigma", sigma)
        params.put("X", np.copy(inputs["X"]))
        stream.put("K", K)
        if self.args["self_kernel"]:
            stream.put("K_self", K.diagonal())

    def _map(self, inputs, stream):
        K = self.evaluate(
            x1=inputs["X"], x2=self.params().get("X"), sigma=self.params().get("sigma")
        )
        stream.put("K", K)
        if self.args["self_kernel"]:
            K_self = self.evaluate(
                x1=inputs["X"], x2=inputs["X"], sigma=self.params().get("sigma"), diagonal_only=True
            )
            stream.put("K_self", K_self)

max_stars_repo_name: rudolfspetrovs/benchml

max_issues_repo_name: rudolfspetrovs/benchml

max_forks_repo_name: rudolfspetrovs/benchml


text:
import numpy as np

def Standardization(data, n_features,axis = 0):
    """
    Performs a data standardization

    Parameters:
    data (np.array): Dataset
    n_features (int): Number of features.
    axis (int): Direction

    Returns:
    np.array: Standardized data
    """
    X = data.copy()

    #Mean subtraction
    mean = np.mean(X, axis = axis)
    for i in range(n_features):
        X[:,i] -= mean[i]

    # Normalize the Variance
    std = np.std(X, axis =axis)
    for i in range(n_features):
        X[:,i] /= std[i]

    return X


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

import math
import numpy as np
import random
import datetime
import subprocess
import pickle
import sys
import os
import argparse
from options import add_neurosat_options
from neurosat import NeuroSAT

parser = argparse.ArgumentParser()
add_neurosat_options(parser)

parser.add_argument('test_dir', action='store', type=str, help='Directory with directories of testation data')
parser.add_argument('restore_id', action='store', type=int)
parser.add_argument('restore_epoch', action='store', type=int)
parser.add_argument('n_rounds', action='store', type=int)

opts = parser.parse_args()
setattr(opts, 'run_id', None)
setattr(opts, 'n_saves_to_keep', 1)

print(opts)

g = NeuroSAT(opts)
g.restore()

results = g.test(opts.test_dir)
for (test_filename, etest_cost, etest_mat) in results:
    print("%s %.6f (%.6f, %.6f, %.6f, %.6f) accuracy %.6f" % (
        test_filename, etest_cost,
        etest_mat.ff, etest_mat.ft, etest_mat.tf, etest_mat.tt,
        etest_mat.ff + etest_mat.tt))

max_stars_repo_name: motus/neurosat

max_issues_repo_name: motus/neurosat

max_forks_repo_name: motus/neurosat


text:
import os
import torch
import torch.nn.functional as F
import glob
import imageio
import numpy as np
import pickle
from util import get_image_to_tensor_balanced, get_mask_to_tensor
import random
from .FluidShakeDataset import FluidShakeDataset
from .WaterPourDataset import FluidPourDataset

class MixedDataset(torch.utils.data.Dataset):
    def __init__(self, stage='train', image_size=(128, 128), img_format='png', seed=10):
        super(MixedDataset, self).__init__()
        self.shake_dataset = FluidShakeDataset(stage=stage, image_size=image_size, img_format=img_format)
        self.pour_dataset = FluidPourDataset(stage=stage, image_size=image_size, img_format=img_format)
        self.dataset_size = self.shake_dataset.__len__() + self.pour_dataset.__len__()

        random.seed(seed)

        order = [i for i in range(self.dataset_size)]
        random.shuffle(order)

        self.order = order

    def __len__(self):
        return self.dataset_size

    def __getitem__(self, item):
        item = self.order[item]
        if item < self.shake_dataset.__len__():
            return self.shake_dataset.__getitem__(item)
        else:
            return self.pour_dataset.__getitem__(item - self.shake_dataset.__len__())

max_stars_repo_name: xavihart/pixel-nerf

max_issues_repo_name: xavihart/pixel-nerf

max_forks_repo_name: xavihart/pixel-nerf


text:
"""
BSD 3-Clause License

Copyright (c) 2018, Jerrad Genson
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

* Neither the name of the copyright holder nor the names of its
  contributors may be used to endorse or promote products derived from
  this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

"""

import numpy as np

from .trainer import Trainer
from learners.learner import LearnerType, InvalidLearnerTypeError

class Backpropagation(Trainer):
    """
    Backpropagation Trainer for the MultilayerPerceptron Learner.

    iterations: Number of iterations to run the training algorithm.
                If this is set too low, the algorithm might not converge
                on a solution. If set too high, it might take too long to
                run and/or overfit the data. Defaults to 1000.
    learning_rate: A float between 0 and 1 that determines the magnitude
                   of updates to the network's weights. A high learning
                   rate will cause the network to converge faster, but
                   might negatively impact the precision/solution quality.
                   Defaults to 0.25.
    randomize: A flag that indicates whether or not to randomize inputs
               and targets. This can improve the speed at which the
               training algorithm converges. Default value is False.
    momentum: The amount of "momentum" to conserve during training as a float
              between 0 and 1. Defaults to 0.9.

    """

    def __init__(self, iterations=1000, learning_rate=0.25, randomize=False,
                 momentum=0.9):

        self.learning_rate = learning_rate
        self.randomize = randomize
        self.momentum = momentum
        return super().__init__(iterations)

    @staticmethod
    def _error_gradient(beta, activations, previous_layer, previous_delta):
        """
        Compute the error gradient for a hidden layer.

        Args
          beta: Learning rate of the network.
          activations: Output values for the current layer.
          previous_layer: Weights of the previous layer.
          previous_delta: Delta value for the previous layer's error gradient.

        Returns
          A delta value for the current layer's error gradient.

        """

        return (activations * beta * (1.0 - activations) *
                (np.dot(previous_delta, np.transpose(previous_layer))))

    def _backpropagate(self, learner, activations, layers, layers_updates,
                       previous_delta):

        """
        Recursively perform backpropagation on each network layer.

        Args
          learner: An instance of learners.learner.Learner.
          activations: A sequence of all the network's activations.
          layers: A sequence of the network's layer weights.
          layers_updates: A sequence of the layer weight updates computed by
                          previous runs of backpropagation.
          previous_delta: The delta values computed for the previous layer.

        Returns
          A tuple of (layer_weights, updates) where layers_weights is a list of
          new network layer weights and updates is a list of new layer weight
          updates.

        """

        layer_activations = activations[-1]
        layer_inputs = activations[-2]
        layer_weights = layers[-2]
        previous_weights = layers[-1]
        layer_delta = self._error_gradient(learner.beta,
                                           layer_activations,
                                           previous_weights,
                                           previous_delta)

        if len(layers) == 2 and learner.bias:
            layer_delta = layer_delta[:, :-1]

        updates = (self.learning_rate *
                   np.dot(np.transpose(layer_inputs), layer_delta) +
                   self.momentum * layers_updates[-1])

        layer_weights -= updates
        if len(layers) == 2:
            # Basis step
            return [layer_weights], [updates]

        else:
            # Inductive step
            next_layer, next_updates = self._backpropagate(learner,
                                                           activations[:-1],
                                                           layers[:-1],
                                                           layers_updates[:-1],
                                                           layer_delta)

            return next_layer + [layer_weights], next_updates + [updates]

    def _train(self, learner, inputs, targets):
        """
        Train the neural network using backpropagation.
        Training happens en batch, which means all the training data is fed to
        the algorithm at once. Mutates self.hidden_weights and
        self.output_weights.

        Args
          learner: An instance of learners.learner.Learner.
          inputs: Training inputs to the network as a numpy array of arrays,
                  where each inner array is one set of inputs.
          targets: Target outputs for the network as a numpy array of arrays,
                   where each inner array is one set of target outputs. Target
                   arrays must match the order of input arrays.

        Returns
          Sum of squares error of the last network recall on the input data.

        """

        # Add the inputs that match the bias node
        training_dataset_rows = np.shape(inputs)[0]
        inputs = learner._concat_bias(inputs)

        # Compute the initial order of input and target nodes so we can
        # randomize them if we so choose.
        node_order = list(range(training_dataset_rows))

        # Create arrays to store update values for weight vectors.
        layers_updates = [np.zeros(layer.shape) for layer in learner.layers]
        for iteration in range(self.iterations):
            outputs = learner._recall(inputs, False)
            error = 0.5 * np.sum((outputs - targets)**2)

            # Compute the output layer error gradient for different activation functions.
            if learner.learner_type == LearnerType.REGRESSION:
                deltao = (outputs - targets) / training_dataset_rows

            elif learner.learner_type == LearnerType.CLASSIFICATION:
                deltao = learner.beta * (outputs - targets) * \
                    outputs * (1.0 - outputs)

            elif learner.learner_type == LearnerType.ONE_OF_N:
                deltao = ((outputs - targets) *
                          (outputs * -outputs + outputs) /
                          training_dataset_rows)

            else:
                raise InvalidLearnerTypeError(
                    'learner_type not member of LearnerType')

            output_updates = (self.learning_rate *
                              np.dot(np.transpose(learner.activations[-2]), deltao) +
                              self.momentum * layers_updates[-1])

            new_output_layer = learner.layers[-1] - output_updates
            new_layers, new_updates = self._backpropagate(learner,
                                                          learner.activations[:-1],
                                                          learner.layers,
                                                          layers_updates[:-1],
                                                          deltao)

            learner.layers = new_layers + [new_output_layer]
            layers_updates = new_updates + [output_updates]
            if self.randomize:
                # Randomize order of input vector and update target vector correspondingly.
                np.random.shuffle(node_order)
                inputs = inputs[node_order, :]
                targets = targets[node_order, :]

        return error

max_stars_repo_name: jerradgenson/rapid-learn

max_issues_repo_name: jerradgenson/rapid-learn

max_forks_repo_name: jerradgenson/rapid-learn


text:
"""This module is for gamma spectrometry analysis. Currently defines a
GammaSpectrum class, reads a .spe file Will in the future have functions
for activity calculations.
"""

from pyne.utils import QA_warn

import numpy as np

from pyne import spectanalysis

QA_warn(__name__)

class GammaSpectrum(spectanalysis.PhSpectrum):
    """GammaSpec class includes gamma specific variables"""

    def __init__(self, real_time=0, live_time=0, det_id='',
                 det_descp='', start_date='', start_time='',
                 calib_e_fit=None, calib_fwhm_fit=None, eres=0, file_name=''):
        """Define gamma spectrum specific variables """
        super(GammaSpectrum, self).__init__()
        self.real_time = real_time
        self.live_time = live_time
        self.dead_time = real_time - live_time
        self.det_id = det_id
        self.det_descp = det_descp
        self.start_date = start_date
        self.start_time = start_time
        self.calib_e_fit = [] if calib_e_fit is None else calib_e_fit
        self.calib_fwhm_fit = [] if calib_fwhm_fit is None else calib_fwhm_fit
        self.eres = eres
        self.file_name = file_name

    def calc_ebins(self):
        """Calculate the energy value for each channel."""
        channels = self.channels = np.asarray(self.channels, float)
        self.ebin = self.calib_e_fit[0] + (self.calib_e_fit[1] * channels) + \
                    (self.calib_e_fit[2] * channels ** 2)

    def __str__(self):
        """Print debug information"""
        print_string = ('Debug print of all header variables\n'
                        'The real time is: {x.real_time}\n'
                        'The live time is: {x.live_time}\n'
                        'The dead time is: {x.dead_time}\n'
                        'Detector ID: {x.det_id}\n'
                        'Detector description: {x.det_descp}\n'
                        'Start date: {x.start_date}\n'
                        'Start time: {x.start_time}\n'
                        'Start channel number: {x.start_chan_num}\n'
                        'Number of channels: {x.num_channels}\n'
                        'Energy calibration fit: {x.calib_e_fit}\n'
                        'FWHM calibration fit: {x.calib_fwhm_fit}\n'
                        'Spectrum: {x.counts}\n'
                        'File name: {x.file_name}').format(x=self)
        return print_string

def read_dollar_spe_file(spec_file_path):
    """Reads a .spe file with the $format
    """

    with open(spec_file_path, "r") as spec_file:
        full_file_text = spec_file.read()
    file_split = full_file_text.splitlines()
    spec_file.close()

    # check version of .spe file matches this functions format
    if (file_split[0] != '$SPEC_ID:'):
        raise RuntimeError('spe file format not supported by this function')

    spectrum = GammaSpectrum()
    # descriptive variables
    spectrum.file_name = spec_file_path
    spectrum.spec_name = file_split[file_split.index("$SPEC_ID:") + 1]
    spectrum.spec_name=spectrum.spec_name.strip()
    tmp = file_split[file_split.index("$SPEC_REM:") + 1]
    tmp = tmp.split(" ")
    spectrum.det_id = tmp[1]
    tmp = file_split[file_split.index("$SPEC_REM:") + 2]
    tmp = tmp.split(" ")
    spectrum.det_descp = tmp[1]
    # time variables
    tmp = file_split[file_split.index("$DATE_MEA:") + 1]
    tmp = tmp.split(" ")
    spectrum.start_date = tmp[0]
    spectrum.start_time = tmp[1]
    tmp = file_split[file_split.index("$MEAS_TIM:") + 1]
    tmp = tmp.split(" ")
    spectrum.real_time = float(tmp[1])
    spectrum.live_time = float(tmp[0])
    # pulse height
    tmp = file_split[file_split.index("$DATA:") + 1]
    tmp = tmp.split(" ")
    spectrum.start_chan_num = int(tmp[0])
    spectrum.num_channels = int(tmp[1])+1
    tmp = file_split[file_split.index("$DATA:") + 2:
                                 file_split.index("$DATA:") + 2
                                 + int(spectrum.num_channels) ]
    for c in tmp:
        val=c.strip()
        spectrum.counts.append(float(val))

    tmp = file_split[file_split.index("$MCA_CAL:") + 2]
    tmp = tmp.split(" ")
    spectrum.calib_e_fit.append(float(tmp[0]))
    spectrum.calib_e_fit.append(float(tmp[1]))
    spectrum.calib_e_fit.append(float(tmp[2]))

    tmp = file_split[file_split.index("$SHAPE_CAL:") + 2]
    tmp = tmp.split(" ")
    spectrum.calib_fwhm_fit.append(float(tmp[0]))
    spectrum.calib_fwhm_fit.append(float(tmp[1]))
    spectrum.calib_fwhm_fit.append(float(tmp[2]))

    # calculate additional parameters based on .spe file
    spectrum.dead_time = spectrum.real_time - spectrum.live_time
    spectrum.channels = np.arange(0, len(spectrum.counts))
    spectrum.calc_ebins()

    return spectrum

def read_spe_file(spec_file_path):
    """Reads a .spe file

    Parameters
    ----------
    spec_file_path : str
        Path to spe file

    Returns
    -------
    spectrum : a GammaSpec object
        Contains all information from .spe file

    """

    spectrum = GammaSpectrum()
    spectrum.file_name = spec_file_path

    with open(spec_file_path, "r") as spec_file:
        full_file_text = spec_file.read()
    file_split = full_file_text.splitlines()
    spec_file.close()
    inspec = False

    # check version of .spe file matches currently supported version
    if (file_split[0] == '$SPEC_ID:'):
        raise RuntimeError('Spe file format is not supported by this function')

    for item in file_split:
        line = item.split(":")
        # processes the spectrum into 2 lists 1 for channel numbers
        # the other for counts
        if (inspec):
            if (len(line) > 1):
                spectrum.channels.append(int(line[0]))
                temp = line[1].strip()
                spectrum.counts.append(float(temp))

        if (line[0] == "Spectrum name"):
            spectrum.spec_name = line[1]
            spectrum.spec_name=spectrum.spec_name.strip()
        elif (line[0] == "Detector ID"):
            spectrum.det_id = line[1].strip()
        elif (line[0] == "Detector description"):
            spectrum.det_descp = line[1].strip()
        elif (line[0] == "Real Time"):
            spectrum.real_time = float(line[1])
        elif (line[0] == "Live Time"):
            spectrum.live_time = float(line[1])
        elif (line[0] == "Acquisition start date"):
            spectrum.start_date = line[1].strip()
        elif (line[0] == "Acquisition start time"):
            spectrum.start_time = line[1].strip() + ":" + line[2] + ":" + line[3]
            spectrum.start_time=spectrum.start_time.strip()
        elif (line[0] == "Starting channel number"):
            spectrum.start_chan_num = int(line[1].strip())
        elif (line[0] == "Number of channels"):
            spectrum.num_channels = int(line[1].strip())
        elif (line[0] == "Energy Fit"):
            temp = line[1].strip()
            temp = temp.split(" ")
            spectrum.calib_e_fit.append(float(temp[0]))
            spectrum.calib_e_fit.append(float(temp[2]))
            spectrum.calib_e_fit.append(float(temp[4]))
        elif (line[0] == "FWHM Fit"):
            temp = line[1].strip()
            temp = temp.split(" ")
            spectrum.calib_fwhm_fit.append(float(temp[0]))
            spectrum.calib_fwhm_fit.append(float(temp[2]))
            spectrum.calib_fwhm_fit.append(float(temp[4]))
        elif (line[0] == "SPECTRUM"):
            inspec = True

    spectrum.counts = np.array(spectrum.counts)
    spectrum.channels = np.array(spectrum.channels)
    # calculate additional parameters based on .spe file
    spectrum.dead_time = spectrum.real_time - spectrum.live_time
    spectrum.calc_ebins()
    return spectrum

def calc_e_eff(energy, eff_coeff, eff_fit=1):
    """Detector efficiency calculation

    Parameters
    ----------
    energy : float
        Energy to calcuate det eff
    eff_coeff : arr
        An array with the coefficients for the energy fit
        the length is not fixed, the length of the array determines the
        number of terms in the expansion
    eff_fit : int
        Determines what type of fit to use

    Returns
    -------
    eff : float
        Value of efficiency for the input energy using the selected fitting eqn

    """
    # eff_fit used to choose between calibration fit eqns
    # energy to be in MeV

    if eff_fit == 1:
        # eff_fit 1 uses series ao + a1(lnE)^1+ a2(lnE)^2+ ....
        log_eff = eff_coeff[0]
        i = 1
        while i < len(eff_coeff):
            log_eff = log_eff + (eff_coeff[i] *
                                 (np.power(np.log(energy), i)))
            i = i + 1
        eff = np.exp(log_eff)
    elif eff_fit == 2:
        # eff_fit 2 uses series a0 + a1(1/E)^1 + a2(1/E)^2+...
        log_eff = eff_coeff[0]
        i = 1
        while i < len(eff_coeff):
            log_eff = log_eff + (eff_coeff[i] * ((1 / energy) ** i))
            i = i + 1
        eff = np.exp(log_eff)
    else:
        raise ValueError('The selected eff_fit is not valid')
        eff = 0

    return eff


text:
# coding=utf-8

import sys
import numpy as np
from pyspark import SparkContext, SparkConf

conf = SparkConf()
sc = SparkContext(conf=conf)

# https://blog.cambridgespark.com/nowadays-recommender-systems-are-used-to-personalize-your-experience-on-the-web-telling-you-what-120f39b89c3c
# https://towardsdatascience.com/various-implementations-of-collaborative-filtering-100385c6dfe0

# aws s3 cp s3://afrikana-collective/1M_Movie_Rating.py ./
# aws s3 cp s3://afrikana-collective/movie-reviews/movies.dat ./
# spark-submit --executor-memory 1g 1M_Movie_Rating.py 260 1493
#
# Item-Item Collaborative Filtering: “Users who liked this movie also liked …”
# User-Item Collaborative Filtering: “Users who are similar to you also liked …”
#

def loadMovieNames():
    movieNames = {}
    with open("movies.dat") as f:
        for line in f:
            fields = line.split("::")
            movieNames[int(fields[0])] = fields[1].decode('ascii','ignore')
    return movieNames

print("\nLoading movie names...\n")
nameDict = loadMovieNames()

movieRatingsRdd = sc.textFile("s3n://afrikana-collective/movie-reviews/ratings.dat")

# Remove duplicate ratings (hint: one half of the n*n matrix)
# Key1,Key2 can be either user or item, based on filtering technique
def isDuplicate(pairedUserRatings):
    ratingsPair = pairedUserRatings[1]
    (key1,rating1) = ratingsPair[0]
    (key2,rating2) = ratingsPair[1]
    return key1 < key2

# Item-based filtering
# (User, Movie => Rating)
movieRatings = movieRatingsRdd\
                .map(lambda line:line.split("::"))\
                .map(lambda line:(int(line[0]),(int(line[1]),float(line[2]))))

# Combine every movie rated by same user
# ratings consists of userID => ((movieID, rating), (movieID, rating))
movieRatingsPartitioned = movieRatings.partitionBy(100)
ratings = movieRatings.join(movieRatingsPartitioned)
uniqueItemRatings = ratings.filter(isDuplicate)

# User-based filtering
# Record: (Movie, User => Rating)
userRatings = movieRatingsRdd\
                .map(lambda line:line.split("::"))\
                .map(lambda line:(int(line[1]),(int(line[0]),float(line[2]))))

# Combine every user rating by same movie
# Ratings consists of movieID => ((userID, rating), (userID, rating))
ratings = userRatings.join(userRatings)
uniqueUserRatings = ratings.filter(isDuplicate)

# Memory based filtering
#
# Item-based filtering: convert ratings from userID => ((movieID, rating), (movieID, rating))
# to userID => ((movieID1, movieID2), (rating1, rating2))
#
# User-based filtering: convert ratings from movieID => ((userID, rating), (userID, rating))
# to movieID => ((userID1, userID2), (rating1, rating2))
#
def pairwiseMovieRatings(pairedRating):
    ratingsPair = pairedRating[1]
    (key1, rating1) = ratingsPair[0]
    (key2, rating2) = ratingsPair[1]
    return (key1,key2),(rating1,rating2)

moviePairs = uniqueItemRatings.map(pairwiseMovieRatings).partitionBy(100)
userPairs = uniqueUserRatings.map(pairwiseMovieRatings).partitionBy(100)

userPairingRatings = userPairs.groupByKey()  # Convert to (user1, user2) => (rating1, rating2),(rating1, rating2)...
moviePairingRatings = moviePairs.groupByKey()  # Convert to (movie1, movie2) => (rating1, rating2),(rating1, rating2)...

def flattenPairings(pairsRatings):
    ratings1, ratings2 = map(list, zip(*pairsRatings))
    return ratings1,ratings2

movieRatingSimilarities = moviePairingRatings.mapValues(flattenPairings)
userRatingSimilarities = userPairingRatings.mapValues(flattenPairings)

# Input as (movie1, movie2) => [rating1.1, rating1.2, ...],[rating2.1, rating2.2, ...]
# Input as (user1, user2) => [rating1.1, rating1.2, ...],[rating2.1, rating2.2, ...]
def getScores(moviePairRatingsTuple):
    numPairs = len(moviePairRatingsTuple[0])
    score = np.dot(moviePairRatingsTuple[0], moviePairRatingsTuple[1]) / (
            np.linalg.norm(moviePairRatingsTuple[0]) * np.linalg.norm(moviePairRatingsTuple[1]))

    return numPairs, score

# Compute similarities in scoring
movieRatingSimilaritieScores = movieRatingSimilarities.mapValues(getScores).persist()
userRatingSimilaritieScores = userRatingSimilarities.mapValues(getScores).persist()

klScoreThreshold = 0.975
coOccurenceMoviesThreshold = 1000 # Similar movie ratings have at least 50 co-occurrences
coOccurenceUserRatingThreshold = 10 #Similar users ratings have at least 10 co-occurrences

# Filter for movies with similarities that are "good" as defined by
# our quality thresholds above

if len(sys.argv) > 1:
    movieID = int(sys.argv[1])
    userID = int(sys.argv[2])

    filteredResultsByMovie = movieRatingSimilaritieScores.filter(lambda result:
                                                               (result[0][0] == movieID or result[0][1] == movieID) and
                                                               result[1][0] > coOccurenceMoviesThreshold and
                                                               result[1][1] > klScoreThreshold)

    filteredResultsByUser = userRatingSimilaritieScores.filter(lambda result:
                                                             (result[0][0] == userID or result[0][1] == userID) and
                                                             result[1][0] > coOccurenceUserRatingThreshold and
                                                             result[1][1] > klScoreThreshold)

    resultsMovies = filteredResultsByMovie\
        .map(lambda result:(result[1][1],(result[0][1] if result[0][1]!=movieID else result[0][0],result[1][0])))\
        .sortByKey(ascending=False)\
        .take(5)

    resultsUsers = filteredResultsByUser \
        .map(lambda result: (result[1][1], (result[0][1] if result[0][1] != userID else result[0][0], result[1][0]))) \
        .sortByKey(ascending=False) \
        .take(5)

    print("\n\nTop 5 similar movies for {}\n".format(nameDict[movieID]))
    for result in resultsMovies:
        print("{}\t\tscore: {:.4f}\t\tvotes: {}".format(nameDict[result[1][0]],float(result[0]),result[1][1]))

    print("\n\nTop 5 similar users for user {}\n".format(userID))
    for result in resultsUsers:
        print("{}\t\tscore: {:.4f}\t\tvotes: {}".format(result[1][0],float(result[0]),result[1][1]))