text:
# Reproduce Fig 4 of the paper

import sys
sys.path.append('../code')
import numpy as np
from scipy import ndimage
import scipy as sp
import matplotlib.pylab as pl
import ot
import scipy.io

#%%

from sklearn.datasets import fetch_olivetti_faces

dataset = fetch_olivetti_faces()
Xtot1 = dataset.data

y=dataset.target

def get_data(x,n):

    idx=np.random.permutation(x.shape[0])

    xr =x[idx[:n],:]
    return xr

#%%
n=400
x=get_data(Xtot1,n)
x=x/x.max()

#%%
from cot import cot_numpy

np.random.seed=1986

def cot_clustering(X,ns,nv,niter_cluster=10, niter=10, algo='emd'):

    Xc=np.random.randn(ns,nv)
    old_cost=0

    for i in range(niter_cluster):
        Ts,Tv,cost=cot_numpy(X, Xc, niter=10, algo=algo,verbose=False)

        Xc=Ts.T.dot(X).dot(Tv)*ns*nv

        print(cost)

        if abs(old_cost-cost)==0:
            break

        old_cost=cost

    return Ts,Tv,Xc

ns=9

nv=40

ot.tic()
Ts,Tv,Xc= cot_clustering(x,ns,nv,niter_cluster=50, niter=20, algo='emd')
ot.toc()

#%% plot var pos

mxc=Xc.mean(0)
isort=np.argsort(mxc)
#isort=np.arange(nv)
Iclass=np.zeros((64,64))

#pl.figure(4,(8,5))

for i in range(nv):
    #pl.subplot(5,8,i+1)
    #pl.imshow(Tv[:,isort[i]].reshape((64,64))*mxc[isort[i]]/np.max(mxc),vmax=Tv.max(),cmap='jet')
    #pl.title('V={:1.2f}'.format(Xc.mean(0)[isort[i]]))
    Iclass[Tv[:,isort[i]].reshape((64,64))>0]=i#*(Tv[:,isort[i]].reshape((64,64))>0)
    #pl.axis('off')

pl.figure(5)
pl.imshow(Iclass,cmap='jet')
pl.title('Pixel (variable) clustering')

#%%

pl.figure(6,(16,10))

clustrs=Tv.dot(Xc.T).T

for i in range(ns):
    pl.subplot(3,4,i+1)
    pl.imshow(clustrs[i,:].reshape((64,64)),cmap='gray')
    pl.axis('off')

#%%

data_rec=Ts.dot(Xc).dot(Tv.T)

idplot=np.random.permutation(y.shape[0])
nbv=100

pl.figure(7,(10,10))
for i in range(nbv):
    pl.subplot(10,10,i+1)
    pl.imshow(x[idplot[i],:].reshape((64,64)),cmap='gray')

pl.figure(8,(10,10))
for i in range(nbv):
    pl.subplot(10,10,i+1)
    pl.imshow(data_rec[idplot[i],:].reshape((64,64)),cmap='gray')

#%%

nx=3
ny=3
delta=1
wimg=64
Icluster=np.zeros((wimg*ny+delta*(ny-1),wimg*nx+delta*(nx-1)))+clustrs.min()
for i in range(ny):
    for j in range(nx):
        Icluster[i*(wimg+delta):i*(wimg+delta)+wimg,j*(wimg+delta):j*(wimg+delta)+wimg]=clustrs[i+ny*j,:].reshape((64,64))

pl.figure(1)
pl.imshow(Icluster,cmap='gray')
pl.title('Sample clusters')
pl.axis('off')

#%%

isel=range(80,400,10)

nx=3
ny=3
delta=1
wimg=64
Idata=np.zeros((wimg*ny+delta*(ny-1),wimg*nx+delta*(nx-1)))+clustrs.min()
for i in range(ny):
    for j in range(nx):
        Idata[i*(wimg+delta):i*(wimg+delta)+wimg,j*(wimg+delta):j*(wimg+delta)+wimg]=x[isel[i+ny*j],:].reshape((64,64))

pl.figure(1)
pl.imshow(Idata,cmap='gray')
pl.title('Face dataset')
pl.axis('off')

#%%

pl.figure(10,(12,4))
pl.clf()
pl.subplot(1,3,1)
pl.imshow(Idata,cmap='gray')
pl.title('Face dataset')
pl.axis('off')

pl.subplot(1,3,2)
pl.imshow(Icluster,cmap='gray')
pl.title('Centroids for sample clustering')
pl.axis('off')

pl.subplot(1,3,3)
pl.imshow(Iclass+1,cmap='nipy_spectral')
pl.title('Pixel (feature) clustering')
pl.axis('off')

left, bottom, width, height = pl.gca().get_position().bounds
cax = pl.gcf().add_axes([left+1.05*width, height*0.23, width*0.05, height])
pl.colorbar( cax=cax)

pl.savefig('./clootclustering_faces.png')
pl.savefig('./clootclustering_faces.pdf',bbox_inches='tight')

max_stars_repo_name: ievred/COOT

max_issues_repo_name: ievred/COOT

max_forks_repo_name: ievred/COOT

reason: import numpy,import scipy,from scipy


text:
from genericpath import isfile
import time
import cv2
from dataset import DatasetRegistry, register_coco, register_balloon
import numpy as np
import os
os.environ['KMP_WARNINGS'] = '0'
os.environ['ENFLAME_LOG_LEVEL'] = 'FATAL'
os.environ['SDK_LOG_LEVEL'] = '3'
import tensorflow as tf
from config import config as cfg
from eval import _paste_mask, DetectionResult
from common import CustomResize, clip_boxes
import tensorpack.utils.viz as tpviz
from tensorpack.utils import logger
from viz import (
    draw_annotation, draw_final_outputs, draw_predictions,
    draw_proposal_recall, draw_final_outputs_blackwhite)
from config import finalize_configs
# import argparse
import TopsInference
model_path = "/tmp/maskrcnn_1600.pb"
if __name__ == '__main__':
    register_coco(cfg.DATA.BASEDIR)
    finalize_configs(is_training=False)
    img = cv2.imread("./sq.jpg", cv2.IMREAD_COLOR)
    orig_shape = img.shape[:2]
    resizer = CustomResize(cfg.PREPROC.TEST_SHORT_EDGE_SIZE, cfg.PREPROC.MAX_SIZE)
    resized_img = resizer.augment(img).astype(np.float32)
    scale = np.sqrt(resized_img.shape[0] * 1.0 / img.shape[0] * resized_img.shape[1] / img.shape[1])
    # resized_img1 = np.load("/tmp/input.npy").astype(np.float32)
    with TopsInference.device(0, 0):
        os.environ['DTU_UMD_FLAGS'] = 'ib_pool_size=209715200'

        ## BACKBONE
        if(os.path.isfile("./engines1600/backbone.exec")):
             engine_backbone = TopsInference.load("./engines1600/backbone.exec")
        else:
            tf_parser_backbone = TopsInference.create_parser(TopsInference.TF_MODEL)
            tf_parser_backbone.set_input_names(['image'])
            tf_parser_backbone.set_input_shapes([[1600, 1600, 3]])
            tf_parser_backbone.set_output_names([
                                'fpn/posthoc_3x3_p2/output',
                                'fpn/posthoc_3x3_p3/output',
                                'fpn/posthoc_3x3_p4/output',
                                'fpn/posthoc_3x3_p5/output',
                                'generate_fpn_proposals/Lvl2/Reshape_1',
                                'generate_fpn_proposals/Lvl3/Reshape_1',
                                'generate_fpn_proposals/Lvl4/Reshape_1',
                                'generate_fpn_proposals/Lvl5/Reshape_1',
                                'generate_fpn_proposals/Lvl6/Reshape_1',
                                'rpn/transpose_1',
                                'rpn_1/transpose_1',
                                'rpn_2/transpose_1',
                                'rpn_3/transpose_1',
                                'rpn_4/transpose_1'
                                ])
            network_backbone = tf_parser_backbone.read(model_path)
            optimizer = TopsInference.create_optimizer()
            engine_backbone = optimizer.build(network_backbone)
            engine_backbone.save_executable("./engines1600/backbone.exec")

        #FC cascade_rcnn_stage1
        if(os.path.isfile("./engines1600/fc_1.exec")):
             engine_fc_1 = TopsInference.load("./engines1600/fc_1.exec")
        else:
            tf_parser_fc_1 = TopsInference.create_parser(TopsInference.TF_MODEL)
            tf_parser_fc_1.set_input_names(['cascade_rcnn_stage1/multilevel_roi_align/output'])
            tf_parser_fc_1.set_input_shapes([[1000, 256, 7, 7]])
            tf_parser_fc_1.set_output_names(['cascade_rcnn_stage1/head/fc7/output'])
            network_fc_1 = tf_parser_fc_1.read(model_path)
            optimizer = TopsInference.create_optimizer()
            engine_fc_1 = optimizer.build(network_fc_1)
            engine_fc_1.save_executable("./engines1600/fc_1.exec")

        #FC cascade_rcnn_stage2
        if(os.path.isfile("./engines1600/fc_2.exec")):
             engine_fc_2 = TopsInference.load("./engines1600/fc_2.exec")
        else:
            tf_parser_fc_2 = TopsInference.create_parser(TopsInference.TF_MODEL)
            tf_parser_fc_2.set_input_names(['cascade_rcnn_stage2/multilevel_roi_align/output'])
            tf_parser_fc_2.set_input_shapes([[1000, 256, 7, 7]])
            tf_parser_fc_2.set_output_names(['cascade_rcnn_stage2/head/fc7/output'])
            network_fc_2 = tf_parser_fc_2.read(model_path)
            optimizer = TopsInference.create_optimizer()
            engine_fc_2 = optimizer.build(network_fc_2)
            engine_fc_2.save_executable("./engines1600/fc_2.exec")

        #FC cascade_rcnn_stage3
        if(os.path.isfile("./engines1600/fc_3.exec")):
             engine_fc_3 = TopsInference.load("./engines1600/fc_3.exec")
        else:
            tf_parser_fc_3 = TopsInference.create_parser(TopsInference.TF_MODEL)
            tf_parser_fc_3.set_input_names(['cascade_rcnn_stage3/multilevel_roi_align/output'])
            tf_parser_fc_3.set_input_shapes([[1000, 256, 7, 7]])
            tf_parser_fc_3.set_output_names(['cascade_rcnn_stage3/head/fc7/output'])
            network_fc_3 = tf_parser_fc_3.read(model_path)
            optimizer = TopsInference.create_optimizer()
            engine_fc_3 = optimizer.build(network_fc_3)
            engine_fc_3.save_executable("./engines1600/fc_0.exec")

        #FC maskrcnn
        if(os.path.isfile("./engines1600/fc_maskrcnn.exec")):
             engine_fc_maskrcnn = TopsInference.load("./engines1600/fc_maskrcnn.exec")
        else:
            tf_parser_fc_maskrcnn = TopsInference.create_parser(TopsInference.TF_MODEL)
            tf_parser_fc_maskrcnn.set_input_names(['multilevel_roi_align/output'])
            tf_parser_fc_maskrcnn.set_input_shapes([[12, 256, 14, 14]])
            tf_parser_fc_maskrcnn.set_output_names(['maskrcnn/fcn3/output'])
            network_fc_maskrcnn = tf_parser_fc_maskrcnn.read(model_path)
            optimizer = TopsInference.create_optimizer()
            engine_fc_maskrcnn = optimizer.build(network_fc_maskrcnn)
            engine_fc_maskrcnn.save_executable("./engines1600/fc_maskrcnn.exec")

        with tf.device("/device:CPU:0"):
            tf.compat.v1.reset_default_graph()
            with tf.compat.v1.Session() as sess:
                with tf.gfile.FastGFile(model_path, 'rb') as f:
                    graph_def = tf.compat.v1.GraphDef()
                    graph_def.ParseFromString(f.read())
                    sess.graph.as_default()
                    tf.import_graph_def(graph_def, name='')
                    sess.run(tf.compat.v1.global_variables_initializer())
                    backbone_time = []
                    total_time = []
                    for i in range(10):
                        print("="*88)
                        print("CURRENT {} ROUND".format(i + 1))
                        print("="*88)
                        time_total_begin = time.time()
                        time_begin = time.time()
                        # # standard input
                        # input = sess.graph.get_tensor_by_name('image:0')

                        # standard output
                        boxes = sess.graph.get_tensor_by_name('output/boxes:0')
                        scores = sess.graph.get_tensor_by_name('output/scores:0')
                        labels = sess.graph.get_tensor_by_name('output/labels:0')
                        masks = sess.graph.get_tensor_by_name('output/masks:0')

                        # backbone output
                        input_fpn_out_2 = sess.graph.get_tensor_by_name(
                            'fpn/posthoc_3x3_p2/output:0')
                        input_fpn_out_3 = sess.graph.get_tensor_by_name(
                            'fpn/posthoc_3x3_p3/output:0')
                        input_fpn_out_4 = sess.graph.get_tensor_by_name(
                            'fpn/posthoc_3x3_p4/output:0')
                        input_fpn_out_5 = sess.graph.get_tensor_by_name(
                            'fpn/posthoc_3x3_p5/output:0')

                        Lvl6_Reshape_1 = sess.graph.get_tensor_by_name(
                            'generate_fpn_proposals/Lvl6/Reshape_1:0')
                        Lvl5_Reshape_1 = sess.graph.get_tensor_by_name(
                            'generate_fpn_proposals/Lvl5/Reshape_1:0')
                        Lvl4_Reshape_1 = sess.graph.get_tensor_by_name(
                            'generate_fpn_proposals/Lvl4/Reshape_1:0')
                        Lvl3_Reshape_1 = sess.graph.get_tensor_by_name(
                            'generate_fpn_proposals/Lvl3/Reshape_1:0')
                        Lvl2_Reshape_1 = sess.graph.get_tensor_by_name(
                            'generate_fpn_proposals/Lvl2/Reshape_1:0')
                        rpn_trans = sess.graph.get_tensor_by_name('rpn/transpose_1:0')
                        rpn1_trans = sess.graph.get_tensor_by_name('rpn_1/transpose_1:0')
                        rpn2_trans = sess.graph.get_tensor_by_name('rpn_2/transpose_1:0')
                        rpn3_trans = sess.graph.get_tensor_by_name('rpn_3/transpose_1:0')
                        rpn4_trans = sess.graph.get_tensor_by_name('rpn_4/transpose_1:0')
                        # crop_resize_1 input

                        # crop_resize_1 output

                        #fc_1
                        fpn_box = sess.graph.get_tensor_by_name('generate_fpn_proposals/boxes:0')
                        stage_1_box = sess.graph.get_tensor_by_name('cascade_rcnn_stage1/output_boxes:0')
                        stage_2_box = sess.graph.get_tensor_by_name('cascade_rcnn_stage2/output_boxes:0')
                        fc_1_in = sess.graph.get_tensor_by_name('cascade_rcnn_stage1/multilevel_roi_align/output:0')
                        fc_1_out = sess.graph.get_tensor_by_name('cascade_rcnn_stage1/head/fc7/output:0')
                        fc_2_in = sess.graph.get_tensor_by_name('cascade_rcnn_stage2/multilevel_roi_align/output:0')
                        fc_2_out = sess.graph.get_tensor_by_name('cascade_rcnn_stage2/head/fc7/output:0')
                        fc_3_in = sess.graph.get_tensor_by_name('cascade_rcnn_stage3/multilevel_roi_align/output:0')
                        fc_3_out = sess.graph.get_tensor_by_name('cascade_rcnn_stage3/head/fc7/output:0')
                        #fc_maskrcnn
                        fc_maskrcnn_in = sess.graph.get_tensor_by_name('multilevel_roi_align/output:0')
                        fc_maskrcnn_out = sess.graph.get_tensor_by_name('maskrcnn/fcn3/output:0')
                        time_end = time.time()
                        duration_tf_find_tensor = time_end - time_begin
                        print("[TIME][CPU] find tensor time-consuming: {:.3f}ms"
                            .format(duration_tf_find_tensor * 1000))
                        # DTU RUN backbone
                        time_begin = time.time()
                        dtu_outputs_backbone = []
                        engine_backbone.run(
                            [resized_img], dtu_outputs_backbone,
                            TopsInference.TIF_ENGINE_RSC_IN_HOST_OUT_HOST)
                        time_end = time.time()
                        duration_dtu_backbone = time_end - time_begin
                        print("[TIME][DTU] backbone time-consuming: {:.3f}ms"
                            .format(duration_dtu_backbone * 1000))
                        backbone_time.append(duration_dtu_backbone * 1000)

                        # CPU RUN fc_1_in
                        time_begin = time.time()
                        foutputs_fc_1 = sess.run(
                                            [
                                                fc_1_in, fpn_box
                                            ],
                                            feed_dict={
                                                # input: resized_img
                                                input_fpn_out_2:dtu_outputs_backbone[0],
                                                input_fpn_out_3:dtu_outputs_backbone[1],
                                                input_fpn_out_4:dtu_outputs_backbone[2],
                                                input_fpn_out_5:dtu_outputs_backbone[3],
                                                Lvl2_Reshape_1:dtu_outputs_backbone[4],
                                                Lvl3_Reshape_1:dtu_outputs_backbone[5],
                                                Lvl4_Reshape_1:dtu_outputs_backbone[6],
                                                Lvl5_Reshape_1:dtu_outputs_backbone[7],
                                                Lvl6_Reshape_1:dtu_outputs_backbone[8],
                                                rpn_trans:dtu_outputs_backbone[9],
                                                rpn1_trans:dtu_outputs_backbone[10],
                                                rpn2_trans:dtu_outputs_backbone[11],
                                                rpn3_trans:dtu_outputs_backbone[12],
                                                rpn4_trans:dtu_outputs_backbone[13],
                                            })
                        time_end = time.time()
                        duration_cpu_fc_1_in = time_end - time_begin
                        print("[TIME][CPU] calc fc_1_in time-consuming: {:.3f}ms"
                            .format(duration_cpu_fc_1_in * 1000))

                        # DTU RUN fc_1_out
                        time_begin = time.time()
                        dtu_outputs_fc_1 = []
                        engine_fc_1.run(
                            [foutputs_fc_1[0]], dtu_outputs_fc_1,
                            TopsInference.TIF_ENGINE_RSC_IN_HOST_OUT_HOST)
                        time_end = time.time()
                        duration_dtu_fc_1 = time_end - time_begin
                        print("[TIME][DTU] fc_1 time-consuming: {:.3f}ms"
                            .format(duration_dtu_fc_1 * 1000))

                        # CPU RUN fc_2_in
                        time_begin = time.time()
                        foutputs_fc_2 = sess.run(
                                            [
                                                fc_2_in, stage_1_box
                                            ],
                                            feed_dict={
                                                # input: resized_img
                                                input_fpn_out_2:dtu_outputs_backbone[0],
                                                input_fpn_out_3:dtu_outputs_backbone[1],
                                                input_fpn_out_4:dtu_outputs_backbone[2],
                                                input_fpn_out_5:dtu_outputs_backbone[3],
                                                fc_1_out:dtu_outputs_fc_1[0],
                                                fpn_box:foutputs_fc_1[1]
                                            })
                        time_end = time.time()
                        duration_cpu_fc_2_in = time_end - time_begin
                        print("[TIME][CPU] calc fc_2_in time-consuming: {:.3f}ms"
                            .format(duration_cpu_fc_2_in * 1000))

                        # DTU RUN fc_2_out
                        time_begin = time.time()
                        dtu_outputs_fc_2 = []
                        engine_fc_2.run(
                            [foutputs_fc_2[0]], dtu_outputs_fc_2,
                            TopsInference.TIF_ENGINE_RSC_IN_HOST_OUT_HOST)
                        time_end = time.time()
                        duration_dtu_fc_2 = time_end - time_begin
                        print("[TIME][DTU] fc_2 time-consuming: {:.3f}ms"
                            .format(duration_dtu_fc_2 * 1000))

                        # CPU RUN fc_3_in
                        time_begin = time.time()
                        foutputs_fc_3 = sess.run(
                                            [
                                                fc_3_in, stage_2_box
                                            ],
                                            feed_dict={
                                                # input: resized_img
                                                input_fpn_out_2:dtu_outputs_backbone[0],
                                                input_fpn_out_3:dtu_outputs_backbone[1],
                                                input_fpn_out_4:dtu_outputs_backbone[2],
                                                input_fpn_out_5:dtu_outputs_backbone[3],
                                                fc_2_out:dtu_outputs_fc_2[0],
                                                stage_1_box:foutputs_fc_2[1]
                                            })
                        time_end = time.time()
                        duration_cpu_fc_3_in = time_end - time_begin
                        print("[TIME][CPU] calc fc_3_in time-consuming: {:.3f}ms"
                            .format(duration_cpu_fc_3_in * 1000))

                        # DTU RUN fc_3_out
                        time_begin = time.time()
                        dtu_outputs_fc_3 = []
                        engine_fc_3.run(
                            [foutputs_fc_3[0]], dtu_outputs_fc_3,
                            TopsInference.TIF_ENGINE_RSC_IN_HOST_OUT_HOST)
                        time_end = time.time()
                        duration_dtu_fc_3 = time_end - time_begin
                        print("[TIME][DTU] fc_3 time-consuming: {:.3f}ms"
                            .format(duration_dtu_fc_3 * 1000))

                        # CPU RUN fc_maskrcnn_in
                        time_begin = time.time()
                        foutputs_fc_maskrcnn = sess.run(
                                            [
                                                fc_maskrcnn_in, boxes, scores, labels,
                                            ],
                                            feed_dict={
                                                # input: resized_img
                                                input_fpn_out_2:dtu_outputs_backbone[0],
                                                input_fpn_out_3:dtu_outputs_backbone[1],
                                                input_fpn_out_4:dtu_outputs_backbone[2],
                                                input_fpn_out_5:dtu_outputs_backbone[3],
                                                fc_1_out:dtu_outputs_fc_1[0],
                                                fc_2_out:dtu_outputs_fc_2[0],
                                                fc_3_out:dtu_outputs_fc_3[0],
                                                stage_2_box:foutputs_fc_3[1]
                                            })
                        time_end = time.time()
                        duration_cpu_fc_maskrcnn_in = time_end - time_begin
                        print("[TIME][CPU] calc fc_maskrcnn_in time-consuming: {:.3f}ms"
                            .format(duration_cpu_fc_maskrcnn_in * 1000))

                        # DTU RUN fc_maskrcnn_out
                        time_begin = time.time()
                        dtu_outputs_fc_maskrcnn = []
                        engine_fc_maskrcnn.run(
                            [foutputs_fc_maskrcnn[0]], dtu_outputs_fc_maskrcnn,
                            TopsInference.TIF_ENGINE_RSC_IN_HOST_OUT_HOST)
                        time_end = time.time()
                        duration_dtu_fc_maskrcnn = time_end - time_begin
                        print("[TIME][DTU] fc_maskrcnn time-consuming: {:.3f}ms"
                            .format(duration_dtu_fc_maskrcnn * 1000))

                        # CPU RUN endding
                        time_begin = time.time()
                        foutputs_end = sess.run(
                                            [
                                                masks,
                                            ],
                                            feed_dict={
                                                fc_maskrcnn_out:dtu_outputs_fc_maskrcnn[0],
                                                labels:foutputs_fc_maskrcnn[3]
                                            })
                        time_end = time.time()
                        duration_cpu_endding = time_end - time_begin
                        print("[TIME][CPU] calc endding time-consuming: {:.3f}ms"
                            .format(duration_cpu_endding * 1000))
                        time_total_end = time.time()
                        print("[TIME] TOTAL time-consuming: {:.3f}ms"
                        .format((time_total_end - time_total_begin) * 1000))
                        total_time.append((time_total_end - time_total_begin) * 1000)

                        boxes = foutputs_fc_maskrcnn[1]
                        scores = foutputs_fc_maskrcnn[2]
                        labels = foutputs_fc_maskrcnn[3]
                        masks = foutputs_end[0]
                    print(np.mean(backbone_time[1:]))
                    print(np.mean(total_time[1:]))
                    # for i in range(4):
                    #     print("$$$$"*30)
                    #     print(foutputs[4 + i].shape)

                    # Some slow numpy postprocessing:
                    boxes = boxes / scale
                    # boxes are already clipped inside the graph, but after the floating point scaling, this may not be true any more.
                    boxes = clip_boxes(boxes, orig_shape)
                    full_masks = [_paste_mask(box, mask, orig_shape)
                                for box, mask in zip(boxes, masks)]
                    masks = full_masks
                    results = [DetectionResult(*args) for args in zip(boxes, scores, labels.tolist(), masks)]
                    final = draw_final_outputs_blackwhite(img, results)
                    viz = np.concatenate((img, final), axis=1)
                    cv2.imwrite("output.png", viz)
                    # logger.info("Inference output for {} written to output.png".format(input_file))
                    tpviz.interactive_imshow(viz)

ext: py

max_stars_repo_name: rlvc/tensorpack

max_issues_repo_name: rlvc/tensorpack

max_forks_repo_name: rlvc/tensorpack

reason: import numpy


text:
# -*- coding: utf-8 -*-
"""Some specialization to the trans functions."""

from .core import pgcore
import numpy as np

__TransCumulative_addForGC__ = pgcore.RTransCumulative.add

def __TransCumulative_addForGC_MP__(self, T, *args):
    """Don't use directly.

    Monkey patch to keep the GC happy until redesign.
    """
    if not hasattr(self, '__trans__'):
        self.__trans__ = []
    self.__trans__.append(T)

    if len(args) == 1:
        # be sure avoid auto conversion from int to IndexArray
        if isinstance(args[0], int):
            return __TransCumulative_addForGC__(self, T, size=args[0])
    return __TransCumulative_addForGC__(self, T, *args)

pgcore.RTransCumulative.add = __TransCumulative_addForGC_MP__

# Aliases
Trans = pgcore.RTrans
TransLinear = pgcore.RTransLinear
TransLin = pgcore.RTransLin
TransPower = pgcore.RTransPower
TransLog = pgcore.RTransLog
TransLogLU = pgcore.RTransLogLU
TransCotLU = pgcore.RTransCotLU
TransCumulative = pgcore.RTransCumulative

class TransSymLog(pgcore.RTrans):
    """Transformation using a bilogarithmic scaling."""

    def __init__(self, tol=1e-12):
        """Forward transformation."""
        super().__init__()
        self.tol = tol

    def trans(self, x):
        """Forward transformation."""
        return pgcore.log(1 + np.abs(x) / self.tol) * np.sign(x)

    def invTrans(self, y):
        """Inverse transformation."""
        return (pgcore.exp(np.abs(y)) - 1.) * self.tol * np.sign(y)

    def deriv(self, x):
        """Derivative of the transformation."""
        return 1. / (np.abs(x) / self.tol + 1) / self.tol

max_stars_repo_name: JuliusHen/gimli

max_issues_repo_name: JuliusHen/gimli

max_forks_repo_name: JuliusHen/gimli


text:
import numpy
from ..default_metrics import DefaultMetrics
from ...utils.helpers import round_value

class ClassificationMetrics(DefaultMetrics):

    def __init__(self, environment):
        self.environment_ = environment

    def metrics_for_team(self, team):
        msg = ""
        if team.extra_metrics_:
            msg += "\n\n### Classification-specific metrics for the best team:"
            msg += "\nrecall per action: "+str(team.extra_metrics_['recall_per_action'])
            msg += "\n\naccuracy: "+str(round_value(team.extra_metrics_['accuracy']))
            msg += "\n\nconfusion matrix:\n"+str(team.extra_metrics_['confusion_matrix'])
            if team.diversity_:
                msg += "\n\ndiversities:"
                for key, value in team.diversity_.iteritems():
                    msg += "\n - "+str(key)+": "+str(value)
        return msg

    def initialize_attributes_for_run_info(self, run_info):
        run_info.recall_per_validation_ = []

    def generate_output_for_attributes_for_run_info(self, run_info):
        msg = ""
        msg += "\n\n\n\n#################### Classification-specific Metrics:"
        msg += "\n\nChampion Recall per Action per Validation: "+str(run_info.recall_per_validation_)
        return msg

    def quick_metrics(self):
        msg = ""
        msg += "\n### Dataset Info:"
        msg += "\ntotal inputs: "+str(self.environment_.total_inputs_)
        msg += "\ntotal actions: "+str(self.environment_.total_actions_)
        msg += "\nactions mapping: "+str(self.environment_.action_mapping_)
        msg += ("\nclass distribution (train set, "+str(len(self.environment_.train_population_))+" "
            "samples): "+str(self.environment_.trainset_class_distribution_))
        msg += ("\nclass distribution (test set, "+str(len(self.environment_.test_population_))+" "
            "samples): "+str(self.environment_.testset_class_distribution_))
        return msg

    def store_per_validation_metrics(self, run_info, best_team, teams_population, programs_population,
        current_generation):
        super(ClassificationMetrics, self).store_per_validation_metrics(run_info, best_team,
            teams_population, programs_population, current_generation)

        run_info.recall_per_validation_.append(best_team.extra_metrics_['recall_per_action'])

    def generate_overall_metrics_output(self, run_infos):
        msg = super(ClassificationMetrics, self).generate_overall_metrics_output(run_infos)
        msg += self._generate_summary(run_infos)
        return msg

    def _generate_summary(self, run_infos): # Obs.: Everything here is duplicated code
        msg = "\n\n\n\n######### SUMMARY:"

        msg += "\n\n\n### Summary Per Run"

        best_scores = [run.global_mean_fitness_per_generation_[-1] for run in run_infos]
        msg += "\n\nGlobal Mean Fitness Score per Training per Run:"
        msg += "\n"+str(best_scores)
        msg += "\nmean: "+str(round_value(numpy.mean(best_scores)))

        best_scores = [run.global_max_fitness_per_generation_[-1] for run in run_infos]
        msg += "\n\nGlobal Max. Fitness Score per Training per Run:"
        msg += "\n"+str(best_scores)
        msg += "\nmean: "+str(round_value(numpy.mean(best_scores)))

        best_scores = [round_value(run.best_team_.score_champion_) for run in run_infos]
        msg += "\n\nChampion Score per Validation per Run:"
        msg += "\n"+str(best_scores)
        msg += "\nmean: "+str(round_value(numpy.mean(best_scores)))
        return msg

max_stars_repo_name: jpbonson/SBBReinforcementLearning

max_issues_repo_name: jpbonson/SBBReinforcementLearner

max_forks_repo_name: jpbonson/SBBReinforcementLearner


text:
"""
Trains Gaussian process emulators.

When run as a script, allows retraining emulators, specifying the number of
principal components, and other options (however it is not necessary to do this
explicitly --- the emulators will be trained automatically when needed).  Run
``python -m src.emulator --help`` for usage information.

Uses the `scikit-learn <http://scikit-learn.org>`_ implementations of
`principal component analysis (PCA)
<http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html>`_
and `Gaussian process regression
<http://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.GaussianProcessRegressor.html>`_.
"""

import logging
import pickle

import numpy as np
from sklearn.decomposition import PCA
from sklearn.externals import joblib
from sklearn.gaussian_process import GaussianProcessRegressor as GPR
from sklearn.gaussian_process import kernels
from sklearn.preprocessing import StandardScaler

from . import cachedir, lazydict, model
from .design import Design

class _Covariance:
    """
    Proxy object to extract observable sub-blocks from a covariance array.
    Returned by Emulator.predict().

    """
    def __init__(self, array, slices):
        self.array = array
        self._slices = slices

    def __getitem__(self, key):
        (obs1, subobs1), (obs2, subobs2) = key
        return self.array[
            ...,
            self._slices[obs1][subobs1],
            self._slices[obs2][subobs2]
        ]

class Emulator:
    """
    Multidimensional Gaussian process emulator using principal component
    analysis.

    The model training data are standardized (subtract mean and scale to unit
    variance), then transformed through PCA.  The first `npc` principal
    components (PCs) are emulated by independent Gaussian processes (GPs).  The
    remaining components are neglected, which is equivalent to assuming they
    are standard zero-mean unit-variance GPs.

    This class has become a bit messy but it still does the job.  It would
    probably be better to refactor some of the data transformations /
    preprocessing into modular classes, to be used with an sklearn pipeline.
    The classes would also need to handle transforming uncertainties, which
    could be tricky.

    """
    #: Observables to emulate as a list of 2-tuples
    #: ``(obs, [list of subobs])``.
    observables = [
        ('dNch_deta', [None]),
        ('dET_deta', [None]),
        ('dN_dy', ['pion', 'kaon', 'proton']),
        ('mean_pT', ['pion', 'kaon', 'proton']),
        ('pT_fluct', [None]),
        ('vnk', [(2, 2), (3, 2), (4, 2)]),
    ]

    def __init__(self, system, npc=10, nrestarts=0):
        logging.info(
            'training emulator for system %s (%d PC, %d restarts)',
            system, npc, nrestarts
        )

        Y = []
        self._slices = {}

        # Build an array of all observables to emulate.
        nobs = 0
        for obs, subobslist in self.observables:
            self._slices[obs] = {}
            for subobs in subobslist:
                Y.append(model.data[system][obs][subobs]['Y'])
                n = Y[-1].shape[1]
                self._slices[obs][subobs] = slice(nobs, nobs + n)
                nobs += n

        Y = np.concatenate(Y, axis=1)

        self.npc = npc
        self.nobs = nobs
        self.scaler = StandardScaler(copy=False)
        self.pca = PCA(copy=False, whiten=True, svd_solver='full')

        # Standardize observables and transform through PCA.  Use the first
        # `npc` components but save the full PC transformation for later.
        Z = self.pca.fit_transform(self.scaler.fit_transform(Y))[:, :npc]

        # Define kernel (covariance function):
        # Gaussian correlation (RBF) plus a noise term.
        design = Design(system)
        ptp = design.max - design.min
        kernel = (
            1. * kernels.RBF(
                length_scale=ptp,
                length_scale_bounds=np.outer(ptp, (.1, 10))
            ) +
            kernels.WhiteKernel(
                noise_level=.1**2,
                noise_level_bounds=(.01**2, 1)
            )
        )

        # Fit a GP (optimize the kernel hyperparameters) to each PC.
        self.gps = [
            GPR(
                kernel=kernel, alpha=0,
                n_restarts_optimizer=nrestarts,
                copy_X_train=False
            ).fit(design, z)
            for z in Z.T
        ]

        # Construct the full linear transformation matrix, which is just the PC
        # matrix with the first axis multiplied by the explained standard
        # deviation of each PC and the second axis multiplied by the
        # standardization scale factor of each observable.
        self._trans_matrix = (
            self.pca.components_
            * np.sqrt(self.pca.explained_variance_[:, np.newaxis])
            * self.scaler.scale_
        )

        # Pre-calculate some arrays for inverse transforming the predictive
        # variance (from PC space to physical space).

        # Assuming the PCs are uncorrelated, the transformation is
        #
        #   cov_ij = sum_k A_ki var_k A_kj
        #
        # where A is the trans matrix and var_k is the variance of the kth PC.
        # https://en.wikipedia.org/wiki/Propagation_of_uncertainty

        # Compute the partial transformation for the first `npc` components
        # that are actually emulated.
        A = self._trans_matrix[:npc]
        self._var_trans = np.einsum(
            'ki,kj->kij', A, A, optimize=False).reshape(npc, nobs**2)

        # Compute the covariance matrix for the remaining neglected PCs
        # (truncation error).  These components always have variance == 1.
        B = self._trans_matrix[npc:]
        self._cov_trunc = np.dot(B.T, B)

        # Add small term to diagonal for numerical stability.
        self._cov_trunc.flat[::nobs + 1] += 1e-4 * self.scaler.var_

    @classmethod
    def from_cache(cls, system, retrain=False, **kwargs):
        """
        Load the emulator for `system` from the cache if available, otherwise
        train and cache a new instance.

        """
        cachefile = cachedir / 'emulator' / '{}.pkl'.format(system)

        # cache the __dict__ rather than the Emulator instance itself
        # this way the __name__ doesn't matter, e.g. a pickled
        # __main__.Emulator can be unpickled as a src.emulator.Emulator
        if not retrain and cachefile.exists():
            logging.debug('loading emulator for system %s from cache', system)
            emu = cls.__new__(cls)
            emu.__dict__ = joblib.load(cachefile)
            return emu

        emu = cls(system, **kwargs)

        logging.info('writing cache file %s', cachefile)
        cachefile.parent.mkdir(exist_ok=True)
        joblib.dump(emu.__dict__, cachefile, protocol=pickle.HIGHEST_PROTOCOL)

        return emu

    def _inverse_transform(self, Z):
        """
        Inverse transform principal components to observables.

        Returns a nested dict of arrays.

        """
        # Z shape (..., npc)
        # Y shape (..., nobs)
        Y = np.dot(Z, self._trans_matrix[:Z.shape[-1]])
        Y += self.scaler.mean_

        return {
            obs: {
                subobs: Y[..., s]
                for subobs, s in slices.items()
            } for obs, slices in self._slices.items()
        }

    def predict(self, X, return_cov=False, extra_std=0):
        """
        Predict model output at `X`.

        X must be a 2D array-like with shape ``(nsamples, ndim)``.  It is passed
        directly to sklearn :meth:`GaussianProcessRegressor.predict`.

        If `return_cov` is true, return a tuple ``(mean, cov)``, otherwise only
        return the mean.

        The mean is returned as a nested dict of observable arrays, each with
        shape ``(nsamples, n_cent_bins)``.

        The covariance is returned as a proxy object which extracts observable
        sub-blocks using a dict-like interface:

        >>> mean, cov = emulator.predict(X, return_cov=True)

        >>> mean['dN_dy']['pion']
        <mean prediction of pion dN/dy>

        >>> cov[('dN_dy', 'pion'), ('dN_dy', 'pion')]
        <covariance matrix of pion dN/dy>

        >>> cov[('dN_dy', 'pion'), ('mean_pT', 'kaon')]
        <covariance matrix between pion dN/dy and kaon mean pT>

        The shape of the extracted covariance blocks are
        ``(nsamples, n_cent_bins_1, n_cent_bins_2)``.

        NB: the covariance is only computed between observables and centrality
        bins, not between sample points.

        `extra_std` is additional uncertainty which is added to each GP's
        predictive uncertainty, e.g. to account for model systematic error.  It
        may either be a scalar or an array-like of length nsamples.

        """
        gp_mean = [gp.predict(X, return_cov=return_cov) for gp in self.gps]

        if return_cov:
            gp_mean, gp_cov = zip(*gp_mean)

        mean = self._inverse_transform(
            np.concatenate([m[:, np.newaxis] for m in gp_mean], axis=1)
        )

        if return_cov:
            # Build array of the GP predictive variances at each sample point.
            # shape: (nsamples, npc)
            gp_var = np.concatenate([
                c.diagonal()[:, np.newaxis] for c in gp_cov
            ], axis=1)

            # Add extra uncertainty to predictive variance.
            extra_std = np.array(extra_std, copy=False).reshape(-1, 1)
            gp_var += extra_std**2

            # Compute the covariance at each sample point using the
            # pre-calculated arrays (see constructor).
            cov = np.dot(gp_var, self._var_trans).reshape(
                X.shape[0], self.nobs, self.nobs
            )
            cov += self._cov_trunc

            return mean, _Covariance(cov, self._slices)
        else:
            return mean

    def sample_y(self, X, n_samples=1, random_state=None):
        """
        Sample model output at `X`.

        Returns a nested dict of observable arrays, each with shape
        ``(n_samples_X, n_samples, n_cent_bins)``.

        """
        # Sample the GP for each emulated PC.  The remaining components are
        # assumed to have a standard normal distribution.
        return self._inverse_transform(
            np.concatenate([
                gp.sample_y(
                    X, n_samples=n_samples, random_state=random_state
                )[:, :, np.newaxis]
                for gp in self.gps
            ] + [
                np.random.standard_normal(
                    (X.shape[0], n_samples, self.pca.n_components_ - self.npc)
                )
            ], axis=2)
        )

emulators = lazydict(Emulator.from_cache)

if __name__ == '__main__':
    import argparse
    from . import systems

    def arg_to_system(arg):
        if arg not in systems:
            raise argparse.ArgumentTypeError(arg)
        return arg

    parser = argparse.ArgumentParser(
        description='train emulators for each collision system',
        argument_default=argparse.SUPPRESS
    )

    parser.add_argument(
        '--npc', type=int,
        help='number of principal components'
    )
    parser.add_argument(
        '--nrestarts', type=int,
        help='number of optimizer restarts'
    )

    parser.add_argument(
        '--retrain', action='store_true',
        help='retrain even if emulator is cached'
    )
    parser.add_argument(
        'systems', nargs='*', type=arg_to_system,
        default=systems, metavar='SYSTEM',
        help='system(s) to train'
    )

    args = parser.parse_args()
    kwargs = vars(args)

    for s in kwargs.pop('systems'):
        emu = Emulator.from_cache(s, **kwargs)

        print(s)
        print('{} PCs explain {:.5f} of variance'.format(
            emu.npc,
            emu.pca.explained_variance_ratio_[:emu.npc].sum()
        ))

        for n, (evr, gp) in enumerate(zip(
                emu.pca.explained_variance_ratio_, emu.gps
        )):
            print(
                'GP {}: {:.5f} of variance, LML = {:.5g}, kernel: {}'
                .format(n, evr, gp.log_marginal_likelihood_value_, gp.kernel_)
            )


text:
import numpy as np

import aesara
import aesara.tensor as tt
from aesara.gof import Apply, Op
from aesara.gradient import DisconnectedType

class PdbBreakpoint(Op):
    """
    This is an identity-like op with the side effect of enforcing a
    conditional breakpoint, inside a aesara function, based on a symbolic
    scalar condition. It automatically detects available debuggers and uses
    the first available in the following order: `pudb`, `ipdb`, or `pdb`.

    :type name: String
    :param name: name of the conditional breakpoint. To be printed when the
                 breakpoint is activated.

    :note: WARNING. At least one of the outputs of the op must be used
                    otherwise the op will be removed from the Aesara graph
                    due to its outputs being unused

    :note: WARNING. Employing the function inside a aesara graph can prevent
                    Aesara from applying certain optimizations to improve
                    performance, reduce memory consumption and/or reduce
                    numerical instability.

            Detailed explanation:
            As of 2014-12-01 the PdbBreakpoint op is not known by any
            optimization. Setting a PdbBreakpoint op in the middle of a
            pattern that is usually optimized out will block the optimization.

    Example:

    .. code-block:: python

        import aesara
        import aesara.tensor as tt
        from aesara.breakpoint import PdbBreakpoint

        input = tt.fvector()
        target = tt.fvector()

        # Mean squared error between input and target
        mse = (input - target) ** 2

        # Conditional breakpoint to be activated if the total MSE is higher
        # than 100. The breakpoint will monitor the inputs, targets as well
        # as the individual error values
        breakpointOp = PdbBreakpoint("MSE too high")
        condition = tt.gt(mse.sum(), 100)
        mse, monitored_input, monitored_target = breakpointOp(condition, mse,
                                                              input, target)

        # Compile the aesara function
        fct = aesara.function([input, target], mse)

        # Use the function
        print fct([10, 0], [10, 5]) # Will NOT activate the breakpoint
        print fct([0, 0], [10, 5]) # Will activate the breakpoint

    """

    __props__ = ("name",)

    def __init__(self, name):
        self.name = name

    def make_node(self, condition, *monitored_vars):

        # Ensure that condition is a aesara tensor
        if not isinstance(condition, aesara.Variable):
            condition = tt.as_tensor_variable(condition)

        # Validate that the condition is a scalar (else it is not obvious how
        # is should be evaluated)
        assert condition.ndim == 0

        # Because the user might be tempted to instantiate PdbBreakpoint only
        # once and apply it many times on different number of inputs, we must
        # create a new instance of the op here, define the instance attributes
        # (view_map and var_types) in that instance and then apply it on the
        # inputs.
        new_op = PdbBreakpoint(name=self.name)
        new_op.view_map = {}
        new_op.inp_types = []
        for i in range(len(monitored_vars)):
            # Every output i is a view of the input i+1 because of the input
            # condition.
            new_op.view_map[i] = [i + 1]
            new_op.inp_types.append(monitored_vars[i].type)

        # Build the Apply node
        inputs = [condition] + list(monitored_vars)
        outputs = [inp.type() for inp in monitored_vars]
        return Apply(op=new_op, inputs=inputs, outputs=outputs)

    def perform(self, node, inputs, output_storage):
        condition = inputs[0]

        if condition:
            try:
                monitored = [np.asarray(inp) for inp in inputs[1:]]
            except Exception:
                raise ValueError(
                    "Some of the inputs to the PdbBreakpoint op "
                    "'%s' could not be casted to NumPy arrays" % self.name
                )

            print("\n")
            print("-------------------------------------------------")
            print("Conditional breakpoint '%s' activated\n" % self.name)
            print("The monitored variables are stored, in order,")
            print("in the list variable 'monitored' as NumPy arrays.\n")
            print("Their contents can be altered and, when execution")
            print("resumes, the updated values will be used.")
            print("-------------------------------------------------")

            try:
                import pudb

                pudb.set_trace()
            except ImportError:
                try:
                    import ipdb

                    ipdb.set_trace()
                except ImportError:
                    import pdb

                    pdb.set_trace()

            # Take the new values in monitored, cast them back to their
            # original type and store them in the output_storage
            for i in range(len(output_storage)):
                output_storage[i][0] = self.inp_types[i].filter(monitored[i])

        else:
            # Simply return views on the monitored variables
            for i in range(len(output_storage)):
                output_storage[i][0] = inputs[i + 1]

    def grad(self, inputs, output_gradients):
        return [DisconnectedType()()] + output_gradients

    def infer_shape(self, inputs, input_shapes):
        # Return the shape of every input but the condition (first input)
        return input_shapes[1:]

    def connection_pattern(self, node):

        nb_inp = len(node.inputs)
        nb_out = nb_inp - 1

        # First input is connected to no output and every other input n is
        # connected to input n-1
        connections = [
            [out_idx == inp_idx - 1 for out_idx in range(nb_out)]
            for inp_idx in range(nb_inp)
        ]
        return connections

max_stars_repo_name: abdalazizrashid/Theano-PyMC

max_issues_repo_name: abdalazizrashid/Theano-PyMC

max_forks_repo_name: abdalazizrashid/Theano-PyMC


text:
# -*- coding: utf-8 -*-
"""
Created on Thu Dec 10 13:39:40 2015

@author: abzooba
"""

import pandas as pd
from textblob import TextBlob
import numpy as np
from textblob.classifiers import NaiveBayesClassifier
from textblob.classifiers import DecisionTreeClassifier
from textblob.classifiers import MaxEntClassifier
from sklearn.cross_validation import StratifiedKFold
from sklearn.metrics import accuracy_score, classification_report, confusion_matrix
import pickle
import sys

train = pd.read_csv('train.csv', index_col = 'id', encoding = 'utf8')
test = pd.read_csv('test.csv', index_col = 'id', encoding = 'utf8')

def packageDataForClassification(sample, feature, target, withIndex = False):
    packaged = []

    for index, row in sample.iterrows():
        features_text = row[feature]
        if withIndex:
            packaged.append((index, features_text, row[target]))
        else:
            print row[target] + ' : ' + features_text
            packaged.append((features_text, row[target]))
    return packaged

k_folds = 5
target = 'cuisine'
y = train[target]
feature = 'ingredients'

fold = 1
skf = StratifiedKFold(y, n_folds=k_folds, shuffle=True)

#for train_index, validation_index in skf:
validation_index, train_index = list(skf)[0]
train_data = train.loc[ train.iloc[train_index].index ]
test_data = train.loc[ train.iloc[validation_index].index ]
#print train_data.cuisine.value_counts()
#print test_data.cuisine.value_counts()
training = packageDataForClassification(train_data, feature, target)
testing = packageDataForClassification(test_data, feature, target, True)
testing_ = packageDataForClassification(test_data, feature, target)
#print training
# Naive Bayes Classifier
naive_clf = NaiveBayesClassifier(training)
#naive_accuracy = naive_clf.accuracy(testing_)

## Decision Tree Classifier
#tree_clf = DecisionTreeClassifier(train)
#tree_accuracy = tree_clf.accuracy(test_)

#for index, feature_x, y in test:
#    self.predictions.loc[ index, 'naive bayes classifier' ] = naive_clf.classify(feature_x)
#    self.predictions.loc[ index, 'decision tree classifier' ] = tree_clf.classify(feature_x)

#print 'fold : ' + str(fold) + ':: Naive Bayes = ' + str(round(naive_accuracy, 4)) + '; Decision Tree = ' + str(round(tree_accuracy, 4))
#    fold += 1

max_stars_repo_name: Bolaka/kaggle-whats-cooking

max_issues_repo_name: Bolaka/kaggle-whats-cooking

max_forks_repo_name: Bolaka/kaggle-whats-cooking


text:
from collections import Counter
from typing import List, Optional, Sequence, Union, Counter as CounterType

import numpy as np
import pandas as pd
from tqdm import tqdm
from scipy import sparse
from sklearn.preprocessing import RobustScaler
from transformers import PreTrainedTokenizerFast
from sklearn.base import BaseEstimator, TransformerMixin

from transformers_domain_adaptation.type import Corpus, Token
from transformers_domain_adaptation.data_selection.metrics import (
    SIMILARITY_FEATURES,
    DIVERSITY_FEATURES,
    similarity_func_factory,
    diversity_func_factory,
)

class DataSelector(BaseEstimator, TransformerMixin):
    """Select subset of data that is likely to be beneficial for domain pre-training.

    This class is sklearn-compatible and implements the sklearn Transformers interface.
    """

    def __init__(
        self,
        keep: Union[int, float],
        tokenizer: PreTrainedTokenizerFast,
        similarity_metrics: Optional[Sequence[str]] = None,
        diversity_metrics: Optional[Sequence[str]] = None,
    ):
        """
        Args:
            keep: Quantity of documents from corpus to keep.
                  To specify number of documents, use :obj:`int`.
                  To specify percentage of documents in corpus, use :obj:`float`.
            tokenizer: A Rust-based 🤗 Tokenizer
            similarity_metrics: An optional list of similarity metrics
            diversity_metrics: An optional list of diversity metrics

        Note:
            For a list of similarity and diversity metrics, refer to :ref:`data-selection-metrics`

        Note:
            At least one similarity/diversity metric must be provided.
        """
        if isinstance(keep, int) and keep <= 0:
            raise ValueError(f"Int value for `keep` must be strictly positive.")
        if isinstance(keep, float) and not 0 <= keep <= 1:
            raise ValueError(
                f"Float value for `keep` must be between 0 and 1 (inclusive)."
            )
        if similarity_metrics is not None:
            _invalid_sim_metrics = set(similarity_metrics) - SIMILARITY_FEATURES
            if _invalid_sim_metrics:
                raise ValueError(
                    f"Invalid similarity metric(s) {_invalid_sim_metrics} found"
                )
        if diversity_metrics is not None:
            _invalid_div_metrics = set(diversity_metrics) - DIVERSITY_FEATURES
            if _invalid_div_metrics:
                raise ValueError(
                    f"Invalid diversity metric(s) {_invalid_div_metrics} found"
                )
        if similarity_metrics is None and diversity_metrics is None:
            raise ValueError(
                f"No metrics provided. Please provide at least one similarity or diversity metric."
            )

        self.keep = keep
        self.tokenizer = tokenizer
        self.similarity_metrics = similarity_metrics
        self.diversity_metrics = diversity_metrics

    def to_term_dist(self, text: str) -> np.ndarray:
        if not len(text.strip()):
            raise ValueError(f"A non-empty string must be provided.")

        tokenized: List[Token] = self.tokenizer.tokenize(text)
        term_counts = Counter(tokenized)

        vocab = self.tokenizer.vocab

        # Create a term distribution
        term_dist: np.ndarray = np.zeros(len(vocab))
        for term, count in term_counts.items():
            term_dist[vocab[term]] = count
        term_dist /= term_dist.sum()

        return term_dist

    def to_term_dist_batch(self, texts: Sequence[str]) -> np.ndarray:
        # * Assumption: Token ID 0 is a special token id and never appears in tokenization with `add_special_tokens=False`

        # Tokenize all documents using Rust tokenizer
        counters: List[CounterType[int]] = [
            Counter(enc.ids)
            for enc in self.tokenizer.backend_tokenizer.encode_batch(
                texts, add_special_tokens=False
            )
        ]

        rows = np.array(
            [val for i, counter in enumerate(counters) for val in [i] * len(counter)]
        )
        cols = np.array([key for counter in counters for key in counter.keys()])
        data = np.array([value for counter in counters for value in counter.values()])
        term_counts = sparse.csr_matrix(
            (data, (rows, cols)),
            shape=(len(counters), len(self.tokenizer)),
            dtype=np.uint16 if len(self.tokenizer) < 2 ** 16 else np.uint32,
        )

        term_dist = term_counts / term_counts.sum(axis=1)
        return np.array(term_dist)

    def fit(self, ft_corpus: Corpus):
        """Compute corpus-level term distribution of :obj:`ft_corpus`.

        A new fitted attribute ``.ft_term_dist_`` of shape (:math:`V`,) is created,
        where :math:`V` is the size of the :obj:`tokenizer` vocabulary.

        Args:
            ft_corpus: The fine-tuning corpus. Not to be confused with
                       the domain pre-training corpus (which is used in :meth:`transform`)

        Note:
            The :obj:`ft_corpus` is treated as a single "document", which will be compared
            against other documents in the in-domain corpus in :meth:`transform`
        """
        self.ft_term_dist_ = self.to_term_dist(" ".join(ft_corpus))
        return self

    def transform(self, docs: Corpus) -> Corpus:
        """Create a relevant subset of documents from the training corpus based on the provided data selection metrics.

        Args:
            docs: The training corpus

        Returns:
            A subset of relevant :obj:`docs` for domain pre-training
        """
        scores = self.compute_metrics(docs)
        composite_scores = scores["composite"].sort_values(ascending=False)

        n_select = (
            self.keep if isinstance(self.keep, int) else int(self.keep * len(docs))
        )
        selection_index = composite_scores.index[:n_select]
        subset_corpus = pd.Series(docs)[selection_index]

        return subset_corpus.tolist()

    def compute_metrics(self, docs: Corpus) -> pd.DataFrame:
        scores = pd.concat(
            [
                self.compute_similarities(docs),
                self.compute_diversities(docs),
            ],
            axis=1,
        )

        # Ensure metrics are normalized, before combining them into a composite score
        scores = pd.DataFrame(
            RobustScaler().fit_transform(scores), columns=scores.columns
        )
        scores["composite"] = scores.sum(axis=1)
        return scores

    def compute_similarities(self, docs: Corpus) -> pd.DataFrame:
        similarities = pd.DataFrame()  # of shape (n_docs, n_metrics)
        if (
            self.similarity_metrics is None
        ):  # Short-circuit function to avoid unnecessary computations
            return similarities

        term_dists = self.to_term_dist_batch(docs)

        pbar = tqdm(
            self.similarity_metrics,
            desc="computing similarity",
            unit="metric",
            dynamic_ncols=True,
        )
        for metric in pbar:
            sim_func = similarity_func_factory(metric)
            similarities[metric] = sim_func(
                term_dists, self.ft_term_dist_.reshape(1, -1)
            )

        return similarities

    def compute_diversities(self, docs: Corpus) -> pd.DataFrame:
        diversities = pd.DataFrame()  # of shape (n_docs, n_metrics)
        if self.diversity_metrics is None:
            return diversities

        tokenized_docs: List[List[Token]] = [
            enc.tokens for enc in self.tokenizer.backend_tokenizer.encode_batch(docs)
        ]

        pbar = tqdm(
            self.diversity_metrics,
            desc="computing diversity",
            unit="metric",
            dynamic_ncols=True,
        )
        for metric in pbar:
            div_func = diversity_func_factory(
                metric,
                train_term_dist=self.ft_term_dist_,
                vocab2id=self.tokenizer.vocab,
            )
            diversities[metric] = pd.Series(
                (div_func(tokenized_doc) for tokenized_doc in tokenized_docs)
            )

        return diversities

reason: import numpy,from scipy


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

from __future__ import (absolute_import, division, print_function,
    generators, nested_scopes, unicode_literals, with_statement)

"""
A basic example of a Linear Programming Problem solved using PuLP (and NumPy).

In this case, decision variables of types LpInteger or LpContinuous yields the
same results.
"""

import numpy as np
from pulp import *

prob = LpProblem("Giapetto's Problem", LpMaximize)
x = [[LpVariable(str('x_{}_{}').format(i, j), cat=LpContinuous)
        for j in range(2) ]
            for i in range(1)]

c   = [3, 2]
A, b = zip(*[
    ( [2, 1], 100 ),
    ( [1, 1],  80 ),
    ( [1, 0],  40 ),
])

A, b, c, x = map(np.array, (A, b, c, x))

obj = lpSum(c * x)
Ax = A * x

for cn in [ lpSum(Axi)<=bi for Axi, bi in zip(Ax,b) ]:
    prob += cn
prob+= obj

status = prob.solve()

# Printing results:

print(prob)
print('-'*50)

print('status:', status)
print('LpStatus[status]:', LpStatus[status])

optimal_x = tuple(zip(x.ravel(), map(value,x.ravel())))

print('Dec. vars:')
for name, val in optimal_x:
    print("  {}: {}".format(name, val))

print('Objective:', value(prob.objective))

# Example inspired by: http://lipas.uwasa.fi/~tsottine/or_with_octave/or_with_octave.pdf

max_stars_repo_name: coppe-or/algorithms

max_issues_repo_name: coppe-or/algorithms

max_forks_repo_name: coppe-or/algorithms


text:
"""
author: Michał Szlupowicz

project is a part of "Dzikie Koty Finansów" cooperation

file contans utilities for gathering the data
"""

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import io
import requests
import os

from random import randint
from time import sleep
from tqdm import tqdm

def recurent_request(quiry, max_depth = 1):
    """
    recurently repeats request in case of fail
    """
    try:
        data = requests.get(quiry).content
    except:
        max_depth -= 1

        if max_depth <= 0:
            print("request fail!!!")
            return None
        else:
            data = recurent_request(quiry, max_depth = max_depth)
    return data

def download_stock_data(name, interval = "d", skip = "1111111", save = False, save_dir = "stock_data"):
    """
    Function for downloading stock data for given users.
    INPUT:
        name - stirng with shortcut of stock name
        interval - stirng with type of interval {"d","w","m","q","y"}
        skip - logical value. If to skip aciditonal data
    """
    adress = "https://stooq.pl/q/d/l/"
    quiry_name = "s=" + str(name)
    quiry_interval = "i=" + str(interval)
    quiry_skip = "o=" + str(skip)
    quiry = adress + "?" + quiry_name + "&" + quiry_interval + "&" + quiry_skip

    data = recurent_request(quiry, 5)#requests.get(quiry).content

    try:
        data_frame = pd.read_csv(io.StringIO(data.decode('utf-8')))
        data_frame["Data"]
        if save:
            directory = save_dir
            file_name = str(name).upper()  + ".tsv"
            path = os.path.join(directory, file_name)

            cwd = os.getcwd()
            dir_path = os.path.join(cwd,directory)

            ### try to create directory
            try:
                os.makedirs(dir_path)
            except:
                pass

            data_frame.to_csv(path, sep = "\t" )
        return data_frame
    except:
        print(f"{name} is broken!!! Check it")
        return 0

def get_stock_names(file = "main_names.csv", directory = "names_data"):
    """
    reads stocks names from file
    """
    path = os.path.join(directory, file)
    df = pd.read_csv(path, sep = "\t")
    df = df.loc[:,["Nazwa", "Skrót"]]
    return df.dropna()

def download_multi_stock(names, save_stats = False, wait = True, save_dir = "stock_data"):
    """
    function downloads and saves stock data for given names
    """

    stats = {"Data": [],
             "Name": []}

    for name in tqdm(names):
        stock = download_stock_data(name, save = True, save_dir = save_dir)
        try:
            stats["Data"].append(stock["Data"][0])
            stats["Name"].append(name)

            ### wait random time
            if wait:
                sleep(randint(10,50)/100)
        except:
            pass

    stats = pd.DataFrame(stats)

    if save_stats:
        stats.to_csv("stats.csv", sep = "\t")

    return stats

def to_one_csv(directory = "stock_data", save = False):
        """
        glues together time series to one file
        """
        files_names = sorted([f.name for f in os.scandir("./"+directory) if f.is_file()])

        data_frames = {}
        names = []

        ### reed all files in directory
        for file_name in tqdm(files_names):
            name = file_name[:-4]
            path = os.path.join(directory,file_name)
            tmp = pd.read_csv(path, sep = "\t", index_col = 0)
            tmp["Data"] = pd.to_datetime(tmp["Data"], dayfirst=True)
            tmp = tmp.set_index('Data')
            data_frames[name] = tmp

        ### get start and end of each series
        firsts = []
        lasts = []
        for key in data_frames.keys():
            tmp = data_frames[key].index
            firsts.append(tmp[0])
            lasts.append(tmp[len(tmp)-1])
        first = sorted(firsts)[0]
        last = sorted(lasts)[-1]

        ### create data frame based on longest time_series
        time_range = pd.date_range(first, last,  freq='24H')
        final_df = pd.DataFrame({"Data": time_range})
        final_df = final_df.set_index('Data')

        ### adds series to data frame by date
        for ii, name in enumerate(data_frames.keys()):
            final_df[name] = data_frames[name]["Zamkniecie"]

        ### drops rows full of NaNs
        final_df = final_df.dropna(how='all')

        ### adds zeros on beginning of series
        for ii, name in enumerate(data_frames.keys()):
            if first != firsts[ii]:
                final_df.loc[first:firsts[ii]-1,name] = 0

        ### interpolate Nas by filling with previous value in time series
        final_df = final_df.fillna(method= "ffill")

        ### save to file
        if save:
            final_df.to_csv("final_df.tsv", sep = "\t")

        return  final_df


text:
# Generic skeleton script for later use

# import the opencv library
import cv2
import numpy as np
import sys
from cv2 import aruco

def customAruco():
    # define an empty custom dictionary with
    aruco_dict = cv2.aruco.custom_dictionary(0, 5, 1)
    # add empty bytesList array to fill with 3 markers later
    aruco_dict.bytesList = np.empty(shape = (5, 4, 4), dtype = np.uint8)
    # add new marker(s)
    mybits = np.array([[1,0,0,0,0],[1,0,0,0,0],[1,0,0,0,0],[1,0,0,0,0],[1,0,1,1,1]], dtype = np.uint8)
    aruco_dict.bytesList[0] = cv2.aruco.Dictionary_getByteListFromBits(mybits)
    mybits = np.array([[1,0,0,0,0],[1,0,0,0,0],[1,0,0,0,0],[1,0,0,0,0],[0,1,0,0,1]], dtype = np.uint8)
    aruco_dict.bytesList[1] = cv2.aruco.Dictionary_getByteListFromBits(mybits)
    mybits = np.array([[1,0,0,0,0],[1,0,0,0,0],[1,0,0,0,0],[1,0,0,0,0],[0,1,1,1,0]], dtype = np.uint8)
    aruco_dict.bytesList[2] = cv2.aruco.Dictionary_getByteListFromBits(mybits)
    mybits = np.array([[1,0,0,0,0],[1,0,0,0,0],[1,0,0,0,0],[1,0,1,1,1],[1,0,0,0,0]], dtype = np.uint8)
    aruco_dict.bytesList[3] = cv2.aruco.Dictionary_getByteListFromBits(mybits)
    mybits = np.array([[1,0,0,0,0],[1,0,0,0,0],[1,0,0,0,0],[1,0,1,1,1],[1,0,1,1,1]], dtype = np.uint8)
    aruco_dict.bytesList[4] = cv2.aruco.Dictionary_getByteListFromBits(mybits)
    # adjust dictionary parameters for better marker detection
    parameters =  cv2.aruco.DetectorParameters_create()
    parameters.cornerRefinementMethod = 3
    parameters.errorCorrectionRate = 0.2
    return aruco_dict, parameters

# define a video capture object
vid = cv2.VideoCapture(0)
aruco_dict, arucoParams = customAruco()

with np.load('res/cal_out/cam_params.npz') as X:
        camera_matrix, dist_coeffs = [X[i] for i in ('mtx','dist')]

while(True):

    # Capture the video frame
    # by frame
    ret, frame = vid.read()

    # Rotate frame 180 degrees
    frame = cv2.rotate(frame, cv2.ROTATE_180)
    gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
    corners, ids, rejectedImgPoints = cv2.aruco.detectMarkers(gray, aruco_dict,
    parameters=arucoParams)

    if len(corners) > 0:
        for i in range(0, len(ids)):  # Iterate in markers
            # Estimate pose of each marker and return the values rvec and tvec---different from camera coefficients
            rvec, tvec, markerPoints = cv2.aruco.estimatePoseSingleMarkers(corners[i], 0.015, camera_matrix,dist_coeffs)
            print("Rvec", rvec)
            print("Tvec",  tvec)
            #cv2.aruco.drawDetectedMarkers(frame, corners, ids)
            #cv2.aruco.drawAxis(frame, camera_matrix, dist_coeffs, rvec, tvec, 0.01)

    # Display the resulting frame
    cv2.imshow('frame', frame)

    # the 'q' button is set as the
    # quitting button you may use any
    # desired button of your choice
    if cv2.waitKey(1) & 0xFF == ord('q'):
        break

# After the loop release the cap object
vid.release()
# Destroy all the windows
cv2.destroyAllWindows()


text:
import unittest
import functools
import numpy as np

from log_uniform import LogUniformSampler

def EXPECT_NEAR(x, y, epsilon):
  return abs(x - y) <= epsilon

class LogUniformSamplerTestCase(unittest.TestCase):
    def test_probability(self):
      N = 1000000
      sampler = LogUniformSampler(N)
      idx = 100
      while idx < N:
        ratio = sampler.probability(idx) / sampler.probability(idx / 2)
        EXPECT_NEAR(ratio, 0.5, 0.1);
        idx *= 2

    def test_unique(self):
      N = 100
      nsampled = 50
      RND = 100

      sampler = LogUniformSampler(N)
      histogram = [0 for idx in range(N)]

      all_values = [idx for idx in range(N)]
      sample_ids, expected, sample_freq = sampler.sample(nsampled, np.asarray(all_values, dtype=np.int32))

      sample_set = set(sample_ids)
      self.assertEqual(len(sample_set), nsampled)

      for rnd in range(RND):
        sample_ids, true_freq, sample_freq = sampler.sample(nsampled, np.asarray(all_values, dtype=np.int32))

        for idx in range(N):
          self.assertTrue(EXPECT_NEAR(expected[idx], true_freq[idx], expected[idx] * 0.5))

        for idx in range(nsampled):
          histogram[sample_ids[idx]] += 1

      for idx in range(N):
        average_count = histogram[idx] / RND
        self.assertTrue(EXPECT_NEAR(expected[idx], average_count, 0.2))

    def test_avoid(self):
      N = 100
      nsampled = 98
      sampler = LogUniformSampler(N)
      labels = [17, 23]
      sample_ids = sampler.sample_unique(nsampled, np.asarray(labels, dtype=np.int32))

      total = functools.reduce(lambda x,y: x+y, sample_ids, 0.0)
      expected_sum = 100 * 99 / 2 - labels[0] - labels[1]
      self.assertEqual(expected_sum, total)

if __name__ == '__main__':
    unittest.main()


text:
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import ImageGrid
import numpy as np
from os import listdir, getcwd
from os import chdir
from PIL import Image
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.image as mimage
from matplotlib.backends.backend_pdf import PdfPages

files = listdir('CNN_run2/Visualisations_w_folders/max_pooling_3')
chdir('CNN_run2/Visualisations_w_folders/max_pooling_3')

images = [Image.open(f).convert('LA') for f in files]

"""
    fig = plt.figure()

    grid = ImageGrid(fig, 111, # similar to subplot(111)
    nrows_ncols = (2, 5), # creates 2x2 grid of axes
    axes_pad=0.1, # pad between axes in inch.
    )
    """

num_rows = 1
num_cols = 128

fig = plt.figure()
gs = gridspec.GridSpec(num_rows, num_cols, wspace=0.0)

i = 0
for g in gs:
    ax = plt.subplot(g)
    ax.imshow(images[i])
    ax.set_xticks([])
    ax.set_yticks([])
    i = i + 1
#    ax.set_aspect('auto')

plt.axis('off')
plt.show()


text:
""" Contains all models specific functions for CUDA Python XVA solver

This file acts as the container for all the question specific function for CUDA
Python SGBM program, namely the forward simulation for the stock models,
thespecific sorting function for bundling and regression basis with its expectation.

.. todo:: Adding the geometric basket payoff and basis
"""
import time
import math
import numpy as np
from numba import cuda, int32, float64
from pyculib import rand

from mathfunctions import star_and_bin_array
from cudafunctions import cuda_jit_step_euler

from example_xva_v2 import Example

# Example related parameters are imported here as golbal variables.
# The usage of golbal variable is due to local variable cannot be used to assign GPU memory.
# The rest of the data is imported in the same way for consistency.
example = Example
# The following parameters are related to stock model.
INITIAL_VALUE = example.initial_value
NUM_OF_ASSETS = example.num_of_assets
NUM_OF_BROWNIAN_MOTION = example.num_of_brownian_motion
MU_BAR = example.mu_bar
SIGMA_BAR = example.sigma_bar
CHOLESKY_DECOMPOSITION = example.cholesky_decomposition
# The following parameters are related to the financial product.
WEIGHT = example.weight
STRIKE = example.strike
BUY_SELL = example.buy_sell
# The following parameters are related to the regression basis.
BASIS_ORDER = example.basis_order

# The followings are stock-related function.
@cuda.jit
def cuda_jit_black_scholes_drift(d_stock, mu_bar, d_drift):
    """Calculate Black-Scholes drift in parallel with GPU"""
    i = cuda.grid(1)
    if i < d_drift.shape[0]:
        for j in range(d_drift.shape[1]):
            d_drift[i, j] = d_stock[i, j] * mu_bar[j]

@cuda.jit
def cuda_jit_black_scholes_volatility(d_stock, sigma_bar, cholesky_decomposition, d_volatility):
    """Calculate Black-Scholes volatility in parallel with GPU"""
    i = cuda.grid(1)
    if i < d_volatility.shape[0]:
        for j in range(d_volatility.shape[1]):
            for k in range(d_volatility.shape[2]):
                d_volatility[i, j, k] = sigma_bar[j]*d_stock[i, j]*cholesky_decomposition[j, k]

# The following are basis-related functions.
@cuda.jit
def cuda_intrinsicvalue(weight, x, res):
    """Calculate the arithmetic weight average of stock prices in parallel with GPU"""
    i = cuda.grid(1)
    if i < x.shape[0]:
        res[i] = 0
        for j in range(x.shape[1]):
            res[i] += weight[j] * x[i, j]

@cuda.jit
def cuda_jit_intrisic_value_basis(weight_stock, basis_partition, degree, basis):
    """Calculate the regression basis based on arithmetic weighted average in parallel with GPU

    This function calculates the arithmetic weighted average of stocks to power 1
    or 2. The resulting values are used as the basis for linear regression.

    :param weight_stock: Array of stock value multipied by weight.
    :param basis_partition: Array of all possible distribution of a given degree in over
    the stock dimension
    :param degree: The power to raise the weighted sum to
    :param basis: Array to hold the resulting value
    :type weight_stock: N x q array
    :type basis_partition: (q + degree -1 C degree) x q array
    :type degree: int
    :type basis: N array
    :returns: basis
    :rtype: N array
    """
    i = cuda.grid(1)
    if i < weight_stock.shape[0]:
        temp = cuda.local.array(1, dtype=float64)
        basis[i] = 0
        for j in range(basis_partition.shape[0]):
            temp[0] = math.gamma(degree + 1.)
            for k in range(basis_partition.shape[1]):
                temp[0] *= weight_stock[i, k] ** basis_partition[j, k] \
                /math.gamma(basis_partition[j, k] +1.)
            basis[i] += temp[0]

@cuda.jit
def cuda_weighted_stock(weight, x, res):
    """Calculate the stock multiplied by weight in parallel with GPU"""
    i = cuda.grid(1)
    if i < x.shape[0]:
        for j in range(x.shape[1]):
            res[i, j] = weight[j] * x[i, j]

@cuda.jit
def cuda_intrisic_value_basis_expect(WEIGHT, x, time, mu, sigma, delta_t, order,
                                     combination, basis_expectation):
    """Calculate the expectation of arithmetic weight average basis in parallel with GPU"""
    i = cuda.grid(1)
    if i < x.shape[0]:
        basis_expectation[i] = arithemticbasketmoment(WEIGHT, x[i], time, mu[i], sigma[i],
                         delta_t, order, combination, basis_expectation[i])

@cuda.jit
def cuda_intrisic_value_basis_brownian_expect(
        weight, x, time, mu, sigma, delta_t, order, dimension, combination,
        basis_brownian_expectation):
    """Calculate the expectation of arithmetic weight average multiplied by brownian motion"""
    i = cuda.grid(1)
    if i < x.shape[0]:
        temp = cuda.local.array(2, dtype=float64)
        temp[0] = 0
        for j in range(weight.shape[0]):
            temp[0] += weight[j] * sigma[i, order, j] # From the derivative and Ito of the basis
        temp[1] = arithemticbasketmoment(weight, x[i], time, mu[i], sigma[i], \
            delta_t, order-1, combination, temp[1])
        basis_brownian_expectation[i] = order * temp[1] * temp[0]

@cuda.jit(device=True)
def arithemticbasketmoment(WEIGHT, x, time, mu, sigma, delta_t, order, combination, rec):
    """Calculate the expectation of arithmetic weight average using the geometric mean

    This functions calculates the expectation of the arithmetrtic weight average by using
    the multinomial theorem. For each possible distribution of degree in dimensions, the
    code calculates the cooresponding multinomial coefficient and calling the geometric
    mean function. Finaaly the result is aggregated.
    """
    if order == 0:
        return 0
    rec = 0.
    coeff = cuda.local.array(1, dtype=float64)
    for i in range(combination.shape[0]):
        coeff[0] = math.gamma(order + 1.)
        for j in range(combination.shape[1]):
            coeff[0] *= WEIGHT[j] ** combination[i, j] /math.gamma(combination[i, j] + 1.)
            # multinomial coefficients times weight
        temp = cuda.local.array(1, dtype=float64)
        temp = cuda_diffusion_product(mu, sigma, order, combination[i], x, time, delta_t, temp)
        rec += coeff[0] * temp[0]
    return rec

# This function may be moved to cudafunctions
@cuda.jit(device=True)
def cuda_diffusion_product(mu, sigma, order, d, x, time, delta_t, temp):
    """ Calculate the expectations of product of a diffusion process

    This function calculates the expectations of product (up to second order) of
    the Euler approximation of a diffusion process in one time step.
    """
    temp[0] = 1
    if order == 1:
        # Result for first order
        for i in range(d.shape[0]):
            temp[0] *= (x[i] + mu[i] * delta_t) ** d[i]
        return temp
    elif order == 2:
        iterator = cuda.local.array(1, dtype=int32)
        iterator[0] = 0
        data = cuda.local.array(2, dtype=int32)
        # Here the active dimensions are determinated.
        for i in range(d.shape[0]):
            if d[i] > 0:
                for j in range(d[i]):
                    data[iterator[0]] = i
                    iterator[0] += 1
                if iterator[0] > 1:
                    break
        # Formula for second order
        temp[0] = x[data[0]] * x[data[1]] + x[data[1]] * mu[data[0]] * delta_t\
        +  x[data[0]] * mu[data[1]] * delta_t + mu[data[0]]  * mu[data[1]] * delta_t **2
        for k in range(sigma.shape[1]):
            temp[0] += sigma[data[0], k] * sigma[data[1], k] * delta_t
        return temp

# Financial product/Payoff-related functions
@cuda.jit
def cuda_put_payoff(x, weight, strike, res):
    """ Calculate the payoff of an arithmetic basket put option"""
    i = cuda.grid(1)
    if i < x.shape[0]:
        temp = cuda.local.array(1, dtype=float64)
        temp[0] = 0.
        for j in range(weight.shape[0]):
            temp[0] += weight[j] * x[i, j]
        res[i, 0] = max(strike - temp[0], 0.0)

@cuda.jit
def cuda_put_terminal_delta(x, weight, strike, volatility, res):
    """ Calculate the terminal Z of the payoff of an arithmetic basket put option"""
    i = cuda.grid(1)
    if i < x.shape[0]:
        temp = cuda.local.array(2, dtype=float64)
        temp[0] = 0.
        temp[1] = 0.
        for k in range(weight.shape[0]):
            temp[1] += weight[k] * x[i, k]
        for j in range(volatility.shape[2]):
            res[i, j] = 0
            for l in range(volatility.shape[1]):
                res[i, j] -= weight[l] * (max((strike - temp[1], 0.0)) > 0) * volatility[i, l, j]

def cuda_jit_montecarlo(timespan, n):
    """Compute the Monte Carlo scenarios

    :param timespan: Array of the time discretization
    :param n: Number of Monte Carlo paths
    :type timespan: M+1 array
    :type n: int
    :returns: (sorting_value, basis, expect_basis, expect_brownian_basis, riskfree_price,
              riskfree_delta, adjusted_price, adjusted_delta)
    :rtypes: (N x (M+1) array, N x Q x M array, N x Q x M array, N x Q x q x M array,
    N x 1 array, N x q array, Nx 1 array, N x q array)
    """
    no_of_timesteps = timespan.shape[0] # M+1
    stock_data = np.ones((n, 1)) * INITIAL_VALUE
    sorting_value = np.empty((n, no_of_timesteps))
    basis = np.empty((n, BASIS_ORDER, no_of_timesteps-1))
    expect_basis = np.empty((n, BASIS_ORDER, no_of_timesteps-1))
    expect_brownian_basis = np.empty((n, BASIS_ORDER, NUM_OF_BROWNIAN_MOTION, no_of_timesteps-1))
    riskfree_price = np.zeros((n, 1))
    riskfree_delta = np.zeros((n, NUM_OF_BROWNIAN_MOTION))
    adjusted_price = np.zeros((n, 1))
    adjusted_delta = np.zeros((n, NUM_OF_BROWNIAN_MOTION))

    blksz = 256
    gridsz = int(math.ceil(float(n) / blksz))

    # instantiate a CUDA stream for queueing async CUDA cmds
    stream = cuda.stream()

    # instantiate a cuRAND PRNG
    prng = rand.PRNG(rndtype=rand.PRNG.MRG32K3A, seed=int(time.time()), stream=stream)

    # Allocate device side array
    d_normdist = cuda.device_array(n * NUM_OF_BROWNIAN_MOTION, dtype=np.double, stream=stream)
    d_drift = cuda.device_array((n, NUM_OF_ASSETS), dtype=np.double, stream=stream)
    d_volatility = \
    cuda.device_array((n, NUM_OF_ASSETS, NUM_OF_BROWNIAN_MOTION), dtype=np.double, stream=stream)
    d_sorting_value = cuda.device_array_like(sorting_value[:, 0], stream=stream)
    d_weighted_stock = cuda.device_array_like(stock_data, stream=stream)
    d_basis = cuda.device_array(n, dtype=np.double, stream=stream)
    d_basis_expectation = cuda.device_array(n, dtype=np.double, stream=stream)
    d_basis_brownian_expectation = cuda.device_array(n, dtype=np.double, stream=stream)
    d_price = cuda.device_array((n, 1), dtype=np.double, stream=stream)
    d_delta = cuda.device_array((n, NUM_OF_BROWNIAN_MOTION), dtype=np.double, stream=stream)

    step_cfg = cuda_jit_step_euler[gridsz, blksz, stream] # The forward simulation scheme

    for k in range(no_of_timesteps):
        # Calculate the time interval
        current_time = timespan[k]
        if k != 0:
            previous_time = timespan[k-1]
            dt = current_time - previous_time
        if k != no_of_timesteps - 1:
            next_time = timespan[k+1]
            forward_time_interval = next_time - current_time

        if k == 0:
            # transfer the initial prices
            d_last = cuda.to_device(stock_data, stream=stream)
        else:
            # call cuRAND to populate d_normdist with gaussian noises
            prng.normal(d_normdist, mean=0, sigma=1)
            # invoke step kernel asynchronously
            step_cfg(d_last, dt, d_drift, d_volatility, d_normdist)

        # Calculate the stock model parameters.
        cuda_jit_drift[gridsz, blksz, stream](d_last, MU_BAR, d_drift)
        cuda_jit_volatility[gridsz, blksz, stream](d_last, SIGMA_BAR,
                           CHOLESKY_DECOMPOSITION, d_volatility)

        # Calculate the sorting value
        cuda_sorting_method[gridsz, blksz, stream](WEIGHT, d_last, d_sorting_value)
        sorting_value[:, k] = d_sorting_value.copy_to_host()

        # Calculate the basis function and its two types of expectations
        if k != 0:
            cuda_weighted_stock[gridsz, blksz, stream](WEIGHT, d_last, d_weighted_stock)
            for j in range(BASIS_ORDER):
                basis_partition = star_and_bin_array(j, NUM_OF_ASSETS)
                cuda_jit_basis[gridsz, blksz, stream](d_weighted_stock, basis_partition, j, d_basis)
                basis[:, j, k-1] = d_basis.copy_to_host()
        if k != no_of_timesteps-1:
            expect_basis[:, 0, k] = 1.
            expect_brownian_basis[:, 0, :, k] = 0.
            for j in range(1, BASIS_ORDER):
                basis_partition = star_and_bin_array(j, NUM_OF_ASSETS)
                basis_gradient_partition = star_and_bin_array(j-1, NUM_OF_ASSETS)
                cuda_basis_expect[gridsz, blksz, stream](WEIGHT, d_last, current_time, d_drift, d_volatility, forward_time_interval, j, basis_partition, d_basis_expectation)
                expect_basis[:, j, k] = d_basis_expectation.copy_to_host()
                for l in range(d_volatility.shape[2]):
                    cuda_basis_brownian_expect[gridsz, blksz, stream](WEIGHT, d_last, current_time, d_drift, d_volatility, forward_time_interval, j, l, basis_gradient_partition, d_basis_brownian_expectation)
                    expect_brownian_basis[:, j, l, k] = d_basis_brownian_expectation.copy_to_host()

        stream.synchronize()   # This ensure all GPU work is completed

    #  Calculate the termianl conditions
    cuda_payoff[gridsz, blksz, stream](d_last, WEIGHT, STRIKE, d_price)
    cuda_terminal_delta[gridsz, blksz, stream](d_last, WEIGHT, STRIKE, d_volatility, d_delta)
    riskfree_price[:, :] = d_price.copy_to_host()
    riskfree_price *= BUY_SELL
    adjusted_price[:,:] = d_price.copy_to_host()
    adjusted_price *= BUY_SELL
    riskfree_delta[:, :] = d_delta.copy_to_host()
    riskfree_delta *= BUY_SELL
    adjusted_delta[:,:] = d_delta.copy_to_host()
    adjusted_delta *= BUY_SELL

    return sorting_value, basis, expect_basis, expect_brownian_basis, riskfree_price, riskfree_delta, adjusted_price, adjusted_delta

# The following is the functions currently under development for the geometric basket option.

# The following functions are related to the regression basis.
@cuda.jit
def cuda_geometric_value(weight, x, res):
    """Calculate the weighted geometric average of stock prices in parallel with GPU"""
    i = cuda.grid(1)
    if i < x.shape[0]:
        res[i] = 1.
        for j in range(x.shape[1]):
            res[i] *= x[i, j] ** weight[j]

@cuda.jit
def cuda_geometric_value_basis(stock, basis_partition, degree, basis):
    """Calculate the geometric average of stock prices to a given order in parallel with GPU"""
    i = cuda.grid(1)
    if i < stock.shape[0]:
        basis[i] = 1
        for j in range(stock.shape[1]):
            basis[i] *= stock[j] ** (degree)

@cuda.jit
def cuda_geometric_value_basis_expect(WEIGHT, x, time, mu, sigma, delta_t, order,
                                     combination, basis_expectation):
    """Calculate the expectation of geometric average basis in parallel with GPU"""
    i = cuda.grid(1)
    if i < x.shape[0]:
        #  Calculate the new drift and volitility of geometric average based on Euler
        #  approximation and Ito formula
        product = cuda.local.array(1, dtype=float64)
        new_mu = cuda.local.array(1, dtype=float64)
        new_sigma = cuda.local.array(NUM_OF_BROWNIAN_MOTION, dtype=float64)
        temp = cuda.local.array(1, dtype=float64)
        product[0] = 1.
        new_mu[0] = 0.
        for j in range(x.shape[1]):
            product *= x[i, j]
        for j in range(x.shape[1]):
            new_mu[0] += mu[i]/x[i, j]
            temp[0] = 0
            for k in range(sigma.shape[2]):
                temp[0] += sigma[i, j, k] * sigma[i, k, j]
            new_mu[0] -= 0.5 * temp/ (x[i, j] ** 2)
        temp[0] = 0.
        for j in range(x.shape[1]):
            for k in range(x.shape[2]):
                for l in range(sigma.shape[2]):
                    temp[0] += sigma[i, j, l] * sigma[i, l, k] / (x[i, j] * x[i, k])
        new_mu[0] += temp[0]
        new_mu[0] *= product[0]
        for j in range(sigma.shape[2]):
            new_sigma[j] = 0.
            for k in range(x.shape[1]):
                new_sigma[j]+= sigma[i, k, j] / x[i, k]
            new_sigma[j] *= product[0]
        #  Call the function to calculate the product of the new process (prod X)^k
        basis_expectation[i] = cuda_diffusion_product(new_mu, new_sigma, order,
                         np.array([order]), product, time, delta_t, basis_expectation[i])

@cuda.jit
def cuda_geometric_value_basis_brownian_expect(
        weight, x, time, mu, sigma, delta_t, order, dimension, combination,
        basis_brownian_expectation):
    """Calculate the expectation of the geometric average basis multiplied by a Brownian motion"""
    i = cuda.grid(1)
    if i < x.shape[0]:
        base_value = cuda.local.array(1, dtype=float64)
        temp = cuda.local.array(2, dtype=float64)
        gradient = cuda.local.array(1, dtype=float64)
        basis_brownian_expectation[i] = 0.
        for j in range(x.shaspe[1]):
            # Calculate the current geometric average to the given power
            base_value[0] = 1.
            for k in range(x.shape[1]):
                base_value[0] *= x[i, k] ** order
            base_value[0] /= x[i, j]
            # Calulation of the derivate of the basis wrt to dimension j
            gradient[0] = 1.
            temp[0] = 0
            for k in range(x.shape[1]):
                temp[0] += mu[i, k]/x[i, k]
                for l in range(sigma.shape[2]):
                    temp[0] -= 0.5 * sigma[i, k, l] * sigma[i, l, j] / x[i, k] **2
            for k in range(x.shape[1]):
                for l in range(x.shape[2]):
                    for m in range(sigma.shape[2]):
                        temp[0] += 0.5 * sigma[i, k, m] * sigma[i, m, l] / (x[i, k] * x[i, l])
            gradient[0] += order * temp[0] * delta_t
            gradient[0] -= mu[i,j] * delta_t/ x[i,j]
            temp[0] = 0.
            for k in range(sigma.shape[1]):
                temp[1] = 0.
                for l in range(sigma.shape[2]):
                    temp[1] += sigma[i, l, k]/x[i, l]
                temp[0] += temp[1]**2
            gradient[0] += order * (order - 1) * 0.5 * temp[0] * delta_t
            temp[0] = 0.
            for k in range(sigma.shape[1]):
                for l in range(sigma.shape[2]):
                    temp[0] += sigma[k,l]*sigma[j, l] / (x[i, k]*x[i, l])
            gradient[0] -= k * temp[0] * delta_t
            temp[0] = 0.
            for k in range(sigma.shape[2]):
                temp[0] += sigma[i, j, k]**2 / x[i, j]
            gradient[0] += temp[0] * delta_t
            gradient[0] *= base_value[0]
            basis_brownian_expectation[i] += gradient[0] * sigma[i, j, dimension]

# Switch for selecting functions for different models
if example.stock_model == "BS":
    cuda_jit_drift = cuda_jit_black_scholes_drift
    cuda_jit_volatility = cuda_jit_black_scholes_volatility
if example.sorting_method == "Intrinsic Value":
    cuda_sorting_method = cuda_intrinsicvalue
elif example.sorting_method == "Geometric Intrinsic Value":
    cuda_sorting_method = cuda_geometric_value
if example.basis == "Intrinsic Value":
    cuda_jit_basis = cuda_jit_intrisic_value_basis
    cuda_basis_expect = cuda_intrisic_value_basis_expect
    cuda_basis_brownian_expect = cuda_intrisic_value_basis_brownian_expect
elif example.basis == "Geometric Intrinsic Value":
    cuda_jit_basis =  cuda_geometric_value_basis
    cuda_basis_expect = cuda_geometric_value_basis_expect
    cuda_basis_brownian_expect = cuda_geometric_value_basis_brownian_expect
if example.put_call == "Put":
    cuda_payoff = cuda_put_payoff
    cuda_terminal_delta = cuda_put_terminal_delta

reason: import numpy,from numba


text:
import unittest
from collections import Counter

import numpy as np

from eskapade.analysis.histogram import Histogram, ValueCounts
from eskapade_python.observers import TestCaseObservable

class HistogramTest(unittest.TestCase, TestCaseObservable):
    def setUp(self):
        pass

    def test_simulate(self):
        size = 10000
        bins = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
        values = [0.03846154, 0.07692308, 0.11538462, 0.15384615, 0.19230769, 0.15384615, 0.11538462, 0.07692308,
                  0.03846154, 0.03846154]
        h = Histogram((values, bins), variable='x')
        data, _ = h.simulate(size)

        self.assertTrue(len(data) == size)
        v, _ = np.histogram(data, bins=bins)
        self.assertTrue(np.allclose(np.divide(v, v.sum()), values, rtol=0, atol=0.02))

    def test_to_normalized(self):
        bins = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
        values = [1, 2, 3, 4, 5, 4, 3, 2, 1, 1]
        h = Histogram((values, bins), variable='x')
        h_norm = h.to_normalized()
        values_norm, _ = h_norm.get_bin_vals()

        self.assertTrue(abs(values_norm.sum() - 1) < 0.01)

    def test_surface(self):
        bins = [0, 1, 3, 4, 6, 7]
        values = [3, 2, 4, 2, 5]
        h = Histogram((values, bins), variable='x')
        s = h.surface()

        self.assertTrue(s == 20)

    def test_constructor1(self):

        # constructor
        cnt = Counter()
        for i in range(10):
            cnt[i * 2] = i

        bin_specs = {'bin_width': 1, 'bin_offset': 0}

        h = Histogram(cnt, variable='x', bin_specs=bin_specs)

        self.assertIsInstance(h, Histogram)

    def test_constructor2(self):

        # constructor
        cnt = Counter()
        for i in range(10):
            cnt[i * 2] = i

        vc = ValueCounts(key='x', counts=cnt)
        bin_specs = {'bin_width': 1, 'bin_offset': 0}

        h = Histogram(vc, variable='x', bin_specs=bin_specs)

        self.assertIsInstance(h, Histogram)

    def test_contents(self):

        # constructor
        cnt = Counter()
        for i in range(10):
            cnt[i * 2] = i

        vc = ValueCounts(key='x', counts=cnt)
        bin_specs = {'bin_width': 1, 'bin_offset': 0}

        h = Histogram(vc, variable='x', bin_specs=bin_specs)

        self.assertEqual(h.n_dim, 1)
        self.assertEqual(h.n_bins, 10)
        self.assertEqual(h._val_counts.sum_counts, 45)

    def test_bin_range(self):

        # constructor
        cnt = Counter()
        for i in range(10):
            cnt[i * 2] = i

        vc = ValueCounts(key='x', counts=cnt)
        bin_specs = {'bin_width': 1, 'bin_offset': 0}

        h = Histogram(vc, variable='x', bin_specs=bin_specs)

        bin_range = (0, 19)
        self.assertTupleEqual(h.get_bin_range(), bin_range)

    def test_bin_edges(self):

        # constructor
        cnt = Counter()
        for i in range(10):
            cnt[i * 2] = i

        vc = ValueCounts(key='x', counts=cnt)
        bin_specs = {'bin_width': 1, 'bin_offset': 0}

        h = Histogram(vc, variable='x', bin_specs=bin_specs)

        # uniform
        bin_edges = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]
        self.assertListEqual(h.get_uniform_bin_edges(), bin_edges)

        # truncated uniform bin edges
        truncated_bin_edges = [5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0]
        self.assertListEqual(h.truncated_bin_edges([5.5, 12.5]), truncated_bin_edges)

        h_bin_edges = h.bin_edges()
        self.assertIsInstance(h_bin_edges, np.ndarray)
        self.assertListEqual(h_bin_edges.tolist(), bin_edges)

    def test_bin_centers(self):

        # constructor
        cnt = Counter()
        for i in range(10):
            cnt[i * 2] = i

        vc = ValueCounts(key='x', counts=cnt)
        bin_specs = {'bin_width': 1, 'bin_offset': 0}

        h = Histogram(vc, variable='x', bin_specs=bin_specs)

        bin_centers = [0.5, 2.5, 4.5, 6.5, 8.5, 10.5, 12.5, 14.5, 16.5, 18.5]
        h_bin_centers = h.bin_centers()
        self.assertIsInstance(h_bin_centers, np.ndarray)
        self.assertListEqual(h_bin_centers.tolist(), bin_centers)

    def test_bin_entries(self):

        # constructor
        cnt = Counter()
        for i in range(10):
            cnt[i * 2] = i

        vc = ValueCounts(key='x', counts=cnt)
        bin_specs = {'bin_width': 1, 'bin_offset': 0}

        h = Histogram(vc, variable='x', bin_specs=bin_specs)

        bin_entries = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
        h_bin_entries = h.bin_entries()
        self.assertIsInstance(h_bin_entries, np.ndarray)
        self.assertListEqual(h_bin_entries.tolist(), bin_entries)

    def test_bin_labels(self):

        # constructor
        cnt = Counter()
        for i in range(10):
            cnt[i * 2] = i

        vc = ValueCounts(key='x', counts=cnt)
        bin_specs = {'bin_width': 1, 'bin_offset': 0}

        h = Histogram(vc, variable='x', bin_specs=bin_specs)

        bin_labels = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18]
        h_bin_labels = h.bin_labels()
        self.assertIsInstance(h_bin_labels, np.ndarray)
        self.assertListEqual(h_bin_labels.tolist(), bin_labels)

    def test_bin_vals(self):

        # constructor
        cnt = Counter()
        for i in range(10):
            cnt[i * 2] = i

        vc = ValueCounts(key='x', counts=cnt)
        bin_specs = {'bin_width': 1, 'bin_offset': 0}

        h = Histogram(vc, variable='x', bin_specs=bin_specs)

        h_bin_vals = h.get_bin_vals()
        self.assertIsInstance(h_bin_vals, tuple)
        self.assertEqual(len(h_bin_vals), 2)
        h_bin_entries, h_bin_edges = h_bin_vals[0], h_bin_vals[1]

        bin_entries = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
        self.assertIsInstance(h_bin_entries, np.ndarray)
        self.assertListEqual(h_bin_entries.tolist(), bin_entries)
        bin_edges = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 19]
        self.assertIsInstance(h_bin_edges, np.ndarray)
        self.assertListEqual(h_bin_edges.tolist(), bin_edges)

        h_bin_vals = h.get_bin_vals(variable_range=[5.5, 12.5])
        self.assertIsInstance(h_bin_vals, tuple)
        self.assertEqual(len(h_bin_vals), 2)
        h_bin_entries, h_bin_edges = h_bin_vals[0], h_bin_vals[1]

        bin_entries = [0, 3, 0, 4, 0, 5, 0, 6]
        self.assertIsInstance(h_bin_entries, np.ndarray)
        self.assertListEqual(h_bin_entries.tolist(), bin_entries)
        bin_edges = [5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0]
        self.assertIsInstance(h_bin_edges, np.ndarray)
        self.assertListEqual(h_bin_edges.tolist(), bin_edges)

    def tearDown(self):
        pass

max_stars_repo_name: mbaak/Eskapade

max_issues_repo_name: mbaak/Eskapade

max_forks_repo_name: mbaak/Eskapade