text:
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import os
import sys
import random
import numpy as np
import logging
import tensorflow as tf
from tensorflow.contrib.keras import utils

from bigfile import BigFile
from constant import *
import utility
from dataset import TestDataBatchIterator
from mlp import MLPPredictor

logger = logging.getLogger(__file__)
formatter_log = "[%(asctime)s - %(filename)s:line %(lineno)s] %(message)s"
logging.basicConfig(
    format = formatter_log,
    datefmt = '%d %b %H:%M:%S')
logger.setLevel(logging.INFO)

flags = tf.app.flags

# GPU option
flags.DEFINE_integer("ses_threads", 2, "Tensorflow CPU session threads to use")
flags.DEFINE_float("gpu_memory_fraction", 0.05, "Fraction of GPU memory to use")
flags.DEFINE_integer("gpu", 1, "select a GPU, 0:gpu0 1:gpu1 (default: 1)")

flags.DEFINE_string("rootpath", ROOT_PATH, "rootdir of the data and model (default: %s)"%ROOT_PATH)
flags.DEFINE_string('train_collection', 'msrvtt10ktrain', 'collection dataset for model training')
flags.DEFINE_string('val_collection', 'msrvtt10kval', 'collection dataset for model validation')
flags.DEFINE_string('test_collection', 'msrvtt10ktest', 'collection dataset for evaluating the performance')
flags.DEFINE_string("model_name", DEFAULT_MODEL_NAME, "model configuration (default: %s)" % DEFAULT_MODEL_NAME)
flags.DEFINE_string('annotation_name', 'conceptsbigram288.txt', 'annotation / lablel name list')
flags.DEFINE_string('vf_name', 'mean_pyresnet-101_rbps13k_step_3_aug_1_dist_2x2,resnet-imagenet-101-0,39,flatten0_output', 'name of the visual feature')
flags.DEFINE_integer('overwrite', 0, 'overwrite existing file (default: 0)')
flags.DEFINE_integer('multi_task', 0, 'multi task learning (default: 0)')

flags.DEFINE_integer('batch_size', 256, 'size of test batch')
flags.DEFINE_integer('top_k', 10, 'the number of top predicted labels')
flags.DEFINE_string('eval_k', '1-5-10', 'the list of top k for evaluation')

flags.DEFINE_string('aux_train_collection', '', 'auxiliary collection dataset for model training')
flags.DEFINE_string('aux_annotation_name', '', 'auxiliary annotation / lablel name list')

FLAGS = flags.FLAGS

def test():
    rootpath = FLAGS.rootpath
    overwrite = FLAGS.overwrite
    train_collection = FLAGS.train_collection
    aux_train_collection = FLAGS.aux_train_collection
    val_collection = FLAGS.val_collection
    test_collection = FLAGS.test_collection
    vf_name = FLAGS.vf_name
    annotation_name = FLAGS.annotation_name
    aux_annotation_name = FLAGS.aux_annotation_name
    top_k = FLAGS.top_k

    config_path = os.path.join(os.path.dirname(__file__), 'model_conf', FLAGS.model_name + '.py')
    config = utility.load_config(config_path)

    model_dir = utility.get_model_dir(FLAGS)
    val_perf_file = os.path.join(model_dir, '%s.txt' % val_collection)
    epoch, val_perf = open(val_perf_file).readline().strip().split()
    epoch = int(epoch.rsplit('_', 1)[1])
    model_path = os.path.join(model_dir, 'epoch_%d.ckpt' % epoch)

    test_data_iter = TestDataBatchIterator(test_collection, vf_name, FLAGS.batch_size, rootpath=rootpath)

    res_dir = utility.get_pred_dir(FLAGS)
    res_file = os.path.join(res_dir, 'id.tagvotes.txt')

    if os.path.exists(res_file) and not overwrite:
        logger.info('%s exists. quit', res_file)
        return

    concept_file = os.path.join(rootpath, train_collection, 'Annotations', annotation_name)
    concepts = map(str.strip, open(concept_file).readlines())

    if FLAGS.multi_task:
        aux_concept_file = os.path.join(rootpath, aux_train_collection, 'Annotations', aux_annotation_name)
        aux_concepts = map(str.strip, open(aux_concept_file).readlines())
    else:
        aux_concepts = []

    mlp_predictor = MLPPredictor(n_layers=[test_data_iter.feat_dim, config.hidden_size, len(concepts)],
                                  last_nl_fun=config.output,
                                  loss_fun=config.loss,
                                  top_k=min(len(concepts),FLAGS.top_k),
                                  ses_threads=FLAGS.ses_threads,
                                  aux_classes=len(aux_concepts),
                                  gpu_memory_fraction=FLAGS.gpu_memory_fraction,
                                  gpu=FLAGS.gpu)

    mlp_predictor.load_model(model_path)

    progbar = utils.Progbar(test_data_iter.num_batches)
    if not os.path.exists(res_dir):
        os.makedirs(res_dir)
    fout = open(res_file, "w")

    for minibatch_index, (ids_test, X_test, _, _) in enumerate(test_data_iter):
        pred_top_values_batch, pred_top_indices_batch = mlp_predictor.predict(X_test)

        for idx, img_id in enumerate(ids_test):
            tagvotes = zip(pred_top_indices_batch[idx], pred_top_values_batch[idx])
            res = ' '.join(['%s %g' % (concepts[x[0]], x[1]) for x in tagvotes])
            fout.write('%s %s\n' % (img_id, res))
        progbar.add(1)
    fout.close()

def main(unused_args):
    test()

if __name__ == '__main__':
    tf.app.run()

ext: py

max_stars_repo_name: li-xirong/coco-cn

max_issues_repo_name: li-xirong/coco-cn

max_forks_repo_name: li-xirong/coco-cn

reason: import numpy


text:
from typing import Any
import numpy as np
from parallellinear.calculations.CalculationsManager import CalculationsManager
from parallellinear.calculations.NumpyLinear import NumpyLinear

class Matrix:

    CALCULATIONS_MANAGER=CalculationsManager(NumpyLinear.getLinearCalculator(), np.float32)

    def __init__(self, rows:int, data:np.ndarray, calcManager=None):
        self.rows = rows
        self.columns=int(len(data)/rows)
        self.data=data
        if calcManager==None:
            self.calcManager = Matrix.CALCULATIONS_MANAGER
        else:
            self.calcManager = calcManager

    @classmethod
    def setCalculationsManager(cls, calculationsManager:CalculationsManager):
        Matrix.CALCULATIONS_MANAGER = calculationsManager

    @classmethod
    def loadCustomFunctionToCalculationsManager(cls, function_name:str, func:str):
        Matrix.CALCULATIONS_MANAGER.getCalculator().loadCustomFunction(function_name, func)

    @classmethod
    def random(cls, rows:int, columns:int, calcManager=None, random_low=0, random_high=1):
        if calcManager==None:
            calcManager = Matrix.CALCULATIONS_MANAGER
        out = cls(rows=rows, data=np.random.rand(rows*columns).astype(calcManager.getPrecision()), calcManager=calcManager)
        if random_low != 0 or random_high != 1:
            out.scale(random_high-random_low)
            out.addScaler(random_low)
        return out

    @classmethod
    def fromFlatListGivenRowNumber(cls, rows:int, data:list, calcManager=None,):
        if calcManager==None:
            calcManager = Matrix.CALCULATIONS_MANAGER
        if len(data) % rows != 0:
                raise ValueError("Matrix from list requires the following assertion to be true len(second parameter) % first parameter == 0")
        return cls(rows=rows, data=np.array(data, dtype=calcManager.getPrecision()), calcManager=calcManager)

    @classmethod
    def zeros(cls, rows:int, columns:int, calcManager=None):
        if calcManager==None:
            calcManager = Matrix.CALCULATIONS_MANAGER
        return cls(rows=rows, data=np.zeros(rows*columns, dtype=calcManager.getPrecision()), calcManager=calcManager)

    @classmethod
    def filledWithValue(cls, rows:int, columns:int, value, calcManager=None):
        if calcManager==None:
            calcManager = Matrix.CALCULATIONS_MANAGER
        return cls(rows=rows, data=np.empty(rows*columns, dtype=calcManager.getPrecision()).fill(value), calcManager=calcManager)

    def getNumberOfColumns(self):
        return self.columns

    def getNumberOfRows(self):
        return self.rows

    def setAtPos(self, x, y, val):
        self.data[x*self.columns+y] = val

    def getAtPos(self, x, y):
        return self.data[x*self.columns+y]

    def __str__(self):
        return str(self.data)

    def __len__(self) -> int:
        return len(self.data)

    def getData(self):
        return self.data

    def exportToListRowNoAtLast(self):
        out=self.exportToList()
        out.append(self.rows)
        return out

    def exportToList(self):
        return self.data.tolist()

    def add(self, a, in_place = True) -> Any:
        if in_place:
            self.calcManager.getCalculator()._addInPlace(self.data, a.getData())
        else:
            return Matrix(self.rows, self.calcManager.getCalculator()._add(self.data, a.getData()), self.calcManager)

    def sub(self, a, in_place = True) -> Any:
        if in_place:
            self.calcManager.getCalculator()._subInPlace(self.data, a.getData())
        else:
            return Matrix(self.rows, self.calcManager.getCalculator()._sub(self, a), self.calcManager)

    def addScaler(self, scaler, in_place = True) -> Any:
        if in_place:
            self.calcManager.getCalculator()._addScalerInPlace(self.data, scaler)
        else:
            return Matrix(self.rows, self.calcManager.getCalculator()._addScaler(self.data, scaler), self.calcManager)

    def subScaler(self, scaler, in_place = True) -> Any:
        if in_place:
            self.calcManager.getCalculator()._subScalerInPlace(self.data, scaler)
        else:
            return Matrix(self.rows, self.calcManager.getCalculator()._subScaler(self.data, scaler), self.calcManager)

    def subScalerFrom(self, scaler, in_place = True) -> Any:
        if in_place:
            self.calcManager.getCalculator()._subScalerFromInPlace(self.data, scaler)
        else:
            return Matrix(self.rows, self.calcManager.getCalculator()._subScalerFrom(self.data, scaler), self.calcManager)

    def scale(self, scaler, in_place = True) -> Any:
        if in_place:
            self.calcManager.getCalculator()._scaleInPlace(self.data, scaler)
        else:
            return Matrix(self.rows, self.calcManager.getCalculator()._scale(self.data, scaler), self.calcManager)

    def descale(self, scaler, in_place = True) -> Any:
        if in_place:
            self.calcManager.getCalculator()._descale(self.data, scaler)
        else:
            return Matrix(self.rows, self.calcManager.getCalculator()._descale(self.data, scaler), self.calcManager)

    def __getitem__(self, indicies):
        return self.data[self.columns * indicies:self.columns*(indicies +1)]

    def transpose(self, in_place = True):

        bufferList = []
        for i in range(0, self.columns):
            for j in range(0, self.rows):
                bufferList.append(self.data[(j*self.columns)+i])

        if in_place:
            self.data = np.array(bufferList, dtype=self.calcManager.getPrecision())
            tmp=self.columns
            self.columns = self.rows
            self.rows = tmp
        else:
            return Matrix.fromFlatListGivenRowNumber(self.columns, bufferList)

    def multiply(self, a):
        return Matrix(self.rows, self.calcManager.getCalculator()._multiply(self.data, a.getData(), self.rows, self.columns, a.getNumberOfRows(), a.getNumberOfColumns()), self.calcManager)

    def applyCustomFunction(self, func_name, in_place = True):
        if in_place:
            self.calcManager.getCalculator()._applyCustomFunctionInPlace(self.data, func_name)
        else:
            return Matrix(self.rows, self.calcManager.getCalculator()._applyCustomFunction(self.data, func_name), self.calcManager)

    def sum(self):
        return self.calcManager.getCalculator()._sum(self.data)

    def elementWiseMultiply(self, a, in_place = True):
        if in_place:
            self.calcManager.getCalculator()._elementWiseMultiplyInPlace(self.data, a.getData())
        else:
            return Matrix(self.rows, self.calcManager.getCalculator()._elementWiseMultiply(self.data, a.getData()), self.calcManager)


text:
import numpy as np
import cv2

# cap = cv2.VideoCapture('MARKfun\IRdrive\data\ir_drive_camcal.avi')
cap = cv2.VideoCapture('data\ir_drive_nowhitecal.avi')
out = cv2.VideoWriter('drive.avi',cv2.VideoWriter_fourcc(*'DIVX'), 15, (240,320))

blackThreshold = 60

while(cap.isOpened()):
    ret, frame = cap.read()
    # print("ret "+str(ret)+" type "+str(type(frame)))
    if not ret or not isinstance(frame,np.ndarray):
        break
    frame = cv2.rotate(frame,cv2.ROTATE_90_CLOCKWISE)
    h,w,_= frame.shape
    img2gray = cv2.cvtColor(frame,cv2.COLOR_BGR2GRAY)
    # img = cv2.medianBlur(img2gray,5)
    # img = cv2.GaussianBlur(img2gray,(15,15),0)
    mmatch = np.zeros((h,256,3))
    # histogram of grey shades
    shades, bin_edges = np.histogram(img2gray, bins=np.arange(256),density=True)
    shademax = shades.max()
    shades = 100 - shades * 100 / shademax;
    # Count black pixels
    ret, mask = cv2.threshold(img2gray, blackThreshold, 255, cv2.THRESH_BINARY_INV)
    maskview = mask.sum(axis=0)
    masksum = mask.sum(axis=0)
    # Matlab type position, get two thirds of maximum to bin
    masksum[masksum < masksum.max()*(2/3)] = 0
    ind = masksum > masksum.max()*(2/3)
    maskpossum = np.sum(np.arange(w)[ind])
    maskposnum = np.sum(np.ones(w)[ind])
    maskposf = maskpossum/maskposnum
    # print("maskpos "+str(maskposf))
    masksum=319-(masksum/masksum.max())*100
    maskview=319-(maskview/maskview.max())*100
    # Draw overlay graphs
    pts = np.vstack((np.arange(w),masksum)).astype(np.int32).T
    ptsx = np.vstack((np.arange(w),maskview)).astype(np.int32).T
    pts2 = np.vstack((np.arange(240),shades[0:240])).astype(np.int32).T
    cv2.polylines(mmatch, [ptsx], isClosed=False, color=(255,0,0))
    cv2.polylines(frame,  [pts], isClosed=False, color=(255,0,0))
    cv2.polylines(frame,  [pts2], isClosed=False, color=(0,256,0))
    cv2.polylines(mmatch, [pts2], isClosed=False, color=(0,256,0))
    pts5 = np.vstack((np.array([blackThreshold,blackThreshold]),np.array([1,100]))).astype(np.int32).T
    cv2.polylines(mmatch, [pts5], isClosed=False, color=(255,0,0))

    pts4 = np.vstack((np.array([w/2,w/2]),np.array([200,300]))).astype(np.int32).T
    cv2.polylines(frame, [pts4], isClosed=True, color=(0,256,0))

    if maskposnum != 0:
        maskpos=int(maskposf)
        pts3 = np.vstack((np.array([maskpos,maskpos]),np.array([200,300]))).astype(np.int32).T
        # cv2.polylines(mmatch, [pts3], isClosed=True, color=(0,256,0))
        cv2.polylines(frame, [pts3], isClosed=True, color=(256,256,0))

    # Try other thresholding
    img = cv2.GaussianBlur(img2gray,(13,13),0)
    maskam = cv2.adaptiveThreshold(img,255,cv2.ADAPTIVE_THRESH_MEAN_C,cv2.THRESH_BINARY_INV,11,2)
    maskag = cv2.adaptiveThreshold(img,255,cv2.ADAPTIVE_THRESH_GAUSSIAN_C,cv2.THRESH_BINARY_INV,11,2)
    ret, masko = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)

    out.write(frame)

    cv2.imshow('frame',frame)
    cv2.imshow('mask',mask)
    cv2.imshow('adaptive mean',maskam)
    cv2.imshow('adaptive gaussian',maskag)
    cv2.imshow('otsu thresholding',masko)
    cv2.imshow('histo',mmatch)

    keyOut = cv2.waitKey(25) & 0xFF
    if keyOut == ord('q'):
        break
    if keyOut == ord('p'):
        # Wait p unpress
        while keyOut == ord('p'):
            keyOut = cv2.waitKey(25) & 0xFF
        while keyOut != ord('p'):
            keyOut = cv2.waitKey(25) & 0xFF

cap.release()
out.release()
cv2.destroyAllWindows()

max_stars_repo_name: masterhapero/MARKfun

max_issues_repo_name: masterhapero/MARKfun

max_forks_repo_name: masterhapero/MARKfun


text:
import numpy as np
import pandas as pd
import os
import argparse
from keras.models import Model
from keras.layers import Input,Conv2D, MaxPooling2D, Flatten, Dense
from keras_preprocessing.image import ImageDataGenerator
from keras.callbacks import TensorBoard, LearningRateScheduler, ModelCheckpoint
from sklearn.metrics import confusion_matrix, classification_report

# input_shape = 64
filters = 32
kernel_size = 3
strides = (3,3)
iterations = 782
# resize_shape = (64,64)

def createFolder(directory):

    if not os.path.exists(directory):
        os.makedirs(directory)

def arguments():
    parser = argparse.ArgumentParser(description = "binary image classification")

    parser.add_argument('--train_path' ,
                        required = True,
                        help = 'train data directory path')
    parser.add_argument('--resize_x',
                        required = False,
                        type = int,
                        default = 32,
                       help = 'resizing x')
    parser.add_argument('--epoch',
                       required = False,
                       type= int,
                        default = 50,
                       help = 'epoch size')
    parser.add_argument('--sample_count',
                       required = False,
                        type = int,
                        default = 5000,
                       help = 'tarining sampe count')
    args  = parser.parse_args()

    return args

def working_model(input_dim, filters , kernel_size , strides):
    inputs = Input(shape = (input_dim, input_dim, 3))
    cnn = Conv2D(filters = filters, kernel_size=kernel_size, strides=strides, activation = 'relu')(inputs)
    cnn = MaxPooling2D(pool_size=(2,2))(cnn)
    cnn = Flatten()(cnn)
    cnn = Dense(128, activation='relu')(cnn)
    output = Dense(1, activation='sigmoid')(cnn)
    classifier = Model(inputs,output)

    return classifier

def image_generator(train_path,  resize_shape):
    train_datagen = ImageDataGenerator(rescale=1./255,
                                      shear_range = 0.2,
                                      zoom_range = 0.2,
                                      horizontal_flip = True,
                                      validation_split = 0.2)
    # test_datagen = ImageDataGenerator(rescale= 1./255)
    training_set = train_datagen.flow_from_directory(train_path,
                                                    target_size = resize_shape,
                                                    batch_size = 32,
                                                    class_mode = 'binary',
                                                    subset = 'training')
    validation_set = train_datagen.flow_from_directory(train_path,
                                                    target_size = resize_shape,
                                                    batch_size = 1,
                                                    class_mode = 'binary',
                                                     subset = 'validation')

    # test_set = test_datagen.flow_from_directory(test_path,
    #                                                 target_size = resize_shape,
    #                                                 batch_size = 32,
    #                                                 class_mode = 'binary')
    return training_set, validation_set

def train(train_set,validation_set ,model, epoch, sample_size):

    model.compile(optimizer = 'adam', loss = 'binary_crossentropy', metrics = ['accuracy'])

#     tb_hist = TensorBoard(log_dir='./graph/image_classification/binary',
#                           histogram_freq=0, write_graph=True, write_images=True)

    tb_cb = TensorBoard(log_dir = './graph/', histogram_freq=0,
                        write_graph=True, write_images=True)

    ckpt = ModelCheckpoint('./weight/ckpt.h5',
                           save_best_only = True, mode = 'auto', period = 10)
    cbks = [tb_cb, ckpt]

    model.fit_generator(train_set,
#                         steps_per_epoch = 15,
#                         epochs = 5,
                       samples_per_epoch = sample_size,
                       epochs = epoch,
                        # steps_per_epoch = iterations,
                        steps_per_epoch = train_set.samples// 32,
                       validation_data = validation_set,
                        # validation_steps=1,
                        validation_steps=validation_set.samples,
                       # nb_val_samples = 2000,
                       verbose = 2,
                       callbacks = cbks,
                       workers = os.cpu_count(),
                       use_multiprocessing = True)
    return model

def score(validation_set, model):
    predict_res = model.predict_generator(validation_set,
                                         verbose = 2,
                                         workers = os.cpu_count(),
                                         use_multiprocessing = True,
                                         steps = validation_set.samples)
    print(validation_set.class_indices)
    def label(x):
        if x>0.5:
            return 1
        else:
            return 0
    label_encode = np.vectorize(label)
    res = label_encode(predict_res)
    # print(np.unique(res))
    # print(res.shape)
    # print(confusion_matrix(validation_set.classes,res))
    print(classification_report(validation_set.classes, res))

if __name__ == '__main__':

    args = arguments()
    createFolder('./graph')
    createFolder('./weight')

    model = working_model(args.resize_x, filters, kernel_size, strides)
    print(model.summary())

    print('Data argumentation for binary classification')
    resize_shape = (args.resize_x, args.resize_x)
    train_datagen, validation_datagen = image_generator(args.train_path, resize_shape)
    model = train(train_datagen ,
                  validation_datagen,
                  model,
                  args.epoch,
                  args.sample_count)
    model.save('./weight/model_weight.h5')

    # predict_res = model.predict_generator()

    score(validation_datagen, model)


text:
#!/usr/bin/env python

import sys
import struct
import numpy as np
import IO

import ISCV

def rle_compress(data):
	'''compress a binary image by run-length encoding.'''
	shape = data.shape
	data = data.reshape(-1)
	tmp = data - np.concatenate(([0],data[:-1]))
	return np.array(np.where(tmp != 0)[0],dtype=np.int32),shape

def rle_decompress(rle, shape):
	'''decompress a run-length encoding into a binary image.'''
	ret = np.zeros(shape,dtype=np.bool)
	data = ret.reshape(-1)
	for c0,c1 in zip(rle[::2],rle[1::2]): data[c0:c1] = True
	return ret

def read_float(data,s):
	return struct.unpack('<d',data[s:s+8])[0],s+8

def read_int(data,s):
	return struct.unpack('<i',data[s:s+4])[0],s+4

def read_half(data,s):
	return struct.unpack('<h',data[s:s+2])[0],s+2

def read_bool(data,s):
	return [False,True][ord(data[s])],s+1

def read_mask(data,s):
	assert(data[s] == chr(2)); s+= 1
	assert(data[s] == chr(1)); s+= 1
	height_a,s = read_int(data,s)
	width_a,s = read_int(data,s)
	chans_a,s = read_int(data,s)
	chans_b,s = read_int(data,s)
	size_bytes,s = read_int(data,s)
	end_s = s + size_bytes
	assert(data[s] == chr(1)); s+= 1
	height,s = read_int(data,s)
	width,s = read_int(data,s)
	chans,s = read_int(data,s)
	size,s = read_int(data,s)
	assert(width*height == size*8)
	assert(chans == 1)
	block = []
	while data[s] != chr(0xef) or data[s+1] != chr(0xbe):
		v = ord(data[s]); s += 1
		count = 1
		if v == 0 or v == 255: count,s = read_int(data,s)
		v = [v&1,v&2,v&4,v&8,v&0x10,v&0x20,v&0x40,v&0x80]
		block.extend(v*count)
	s += 4
	assert(s == end_s)
	assert(len(block) == size*8)
	return np.array(block,dtype=np.bool).reshape(height,width),s

def read_addendum(data,s):
	#assert(data[s] == chr(0x0a))
	ret = []
	func = {'h':read_half,'i':read_int,'d':read_float,'b':read_bool}
	for t in 'hibddddiiiddbddddddiiidiiiiibiiiddiddihiiibiiihibiibibddidiiibiididdiiiid':
		v,s=func[t](data,s)
		print t,v
		ret.append(v)
	return ret,s

def read_CAL(filename):
	'''Reader for .cal format (optitrack).'''
	data = open(filename,'rb').read()
	s = 0
	fourcc,s = read_int(data,s)
	assert(fourcc == 844258806)
	version,s = read_int(data,s)
	assert(version == 20)
	num_cameras,s = read_int(data,s)
	#print num_cameras
	cameras = []
	for ci in xrange(num_cameras):
		#print ' '.join(map(lambda x:hex(x)[-2:].replace('x','0'),map(ord,data[s:s+10])))
		camera_tag,s = read_half(data,s)
		assert(camera_tag == 0x4)
		cam = {}
		# NB because the images are square, we can't be sure which is width and which is height.
		# NB Undecoded CAMERA_0.77 and CAMERA_0.93 parameters are not constants
		for t,k in zip('iiiiiiiiiibiiiiibbbbdbibbiiidiiiiddddddddddddrmiiiid',['INDEX','ID','i11','i0_a','i0_b','HEIGHT','WIDTH',\
					   '?i225_600','?i167_200','strobe_output_Watts','b0_a','i4','fps','i0_c','i0_d','i100','b1_a','b0_b','b0_c','b0_d','d0','b0_e','i3_b',\
					   'b1_b','b1_c','i2','?i1_4_a','i2000','d0.6000000238418579','i48','i64','i4','i2','OX','OY','FX','FY',\
					   'K1','K2','K3','P1','P2','TX','TY','TZ','R','MASK','i64','i1_b','i1_c','i1_d','d0.5']):
			if t == 'm':
				v,s = read_mask(data,s)
				rle,shape = rle_compress(v)
				v2 = rle_decompress(rle, shape)
				assert(np.all(v == v2))
				v = rle,shape
			elif t == 'r':
				v = np.fromstring(data[s:s+9*8],dtype=np.float64).reshape(3,3)
				s += 9*8
			else:
				assert(data[s] == chr(1)); s += 1
				if t == 'i': v,s = read_int(data,s)
				elif t == 'b':v,s = read_bool(data,s)
				elif t == 'd':v,s = read_float(data,s)
				else: raise
			# if we get a file with different values, cause a crash
			# we want to know what these values might mean!
			if k.startswith('?'): print k,v
			if k.startswith('i'): assert(v == int(k[1:].partition('_')[0])); continue
			if k.startswith('d'): assert(v == float(k[1:].partition('_')[0])); continue
			if k.startswith('b0'): assert(v is False); continue
			if k.startswith('b1'): assert(v is True); continue
			cam[k] = v
		cameras.append(cam)
	# no idea what this data is
	addenda = []
	num_addenda,s = read_int(data,s)
	print 'num_addenda',num_addenda
	if 0:
		for ai in xrange(num_addenda):
			addendum,s = read_addendum(data,s)
			addenda.append(addendum)
	#assert(s == len(data))
	return {'filename':filename,'cameras':cameras,'addenda':addenda}

def load_CAL(filename):
	ret = read_CAL(filename)['cameras']
	#print ret[0]; exit()
	mats = []
	for d in ret:
		hsx,hsy = d['WIDTH']*0.5,d['HEIGHT']*0.5
		sensor = np.array([d['WIDTH'],d['HEIGHT']],dtype=np.int32)
		fx = d['FX'] / hsx
		fy = d['FY'] / hsy
		skew = 0. #d['SKEW'] / hsx
		ox,oy = d['OX'],d['OY']
		ox = (ox - hsx) / hsx
		oy = (oy - hsy) / hsx
		K = np.array([[fx, skew, -ox],[0, fy, oy], [0, 0, 1]], dtype=np.float32)
		RT = np.zeros((3,4),dtype=np.float32)
		RT[:3,:3] = d['R'].T
		pos = np.array([d['TX'],d['TY'],d['TZ']],dtype=np.float32)*1000.
		RT[:,3] = -np.dot(RT[:3,:3], pos)
		k1 = d['K1'] * 0.25 # their image goes from -0.5 to +0.5 whereas ours goes from -1 to +1
		k2 = d['K2'] * (0.25*0.25)
		mats.append([K,RT,np.dot(K,RT),np.array([k1,k2],dtype=np.float32),pos,sensor]) # 9+12+12+2+3+2 = 40 numbers
	order = np.argsort([c['INDEX'] for c in ret])
	print 'order',order#np.array([c['INDEX'] for c in ret])[order]
	return [mats[oi] for oi in order],[ret[oi] for oi in order]

if __name__ == '__main__':
	import os
	import sys
	import IO
	import ASFReader
	filename = sys.argv[1]
	cal = load_CAL(filename)
	cams = cal['cameras']
	keys = cams[0].keys()
	#for key in keys:
	#	vals = [c[key] for c in cams]
	#	print key,vals
	IO.save('cal.imscal',cal)
	import pylab as pl
	mask = cal['cameras'][0]['MASK']
	mask = rle_decompress(mask[0],mask[1])
	print mask.shape
	pl.imshow(mask)
	pl.hold()
	pl.show()

max_stars_repo_name: davidsoncolin/IMS

max_issues_repo_name: davidsoncolin/IMS

max_forks_repo_name: davidsoncolin/IMS


text:
import networkx as nx
import json
import codecs
from collections import Counter

def createRelationNetwork():

    # struct node{
    # key = author's name
    # type = 'V'/'T' ;; V means author isn't in training.txt
    # interest = []
    # }

    with codecs.open("author_cooperators.json","r","utf-8") as fid:
        author_cooperator = json.load(fid)
    with codecs.open("author_interest.json","r","utf-8") as fid:
        author_interest = json.load(fid)

    relation_net = nx.Graph()
    # graph to record the realtions between the authors

    for key,value in author_cooperator.items():
        if key in author_interest.keys():
            relation_net.add_node(key,type='T',interest = author_interest[key])
        else:
            relation_net.add_node(key,type='V',interest = [])
        value_count = Counter(value)
        # create graph , set authors' name as node, set the times that they
        # cooperated as the weighted edge

        for k,v in value_count.items():
            relation_net.add_weighted_edges_from([(key,k,v)])
    return relation_net

def expandInterest(name,relation_net):
    # input authors' name
    # expand interest
    expand_interest = []
    cooperators = nx.all_neighbors(relation_net,name)
    for v in cooperators:
        expand_interest.extend(relation_net.node[v]['interest']*relation_net[name][v]['weight'])
    expand_interest = Counter(expand_interest).most_common(5)
    interest = []
    for i in range (len(expand_interest)):
        interest.append(expand_interest[i][0])
    relation_net.node[name]['interest'].extend(interest)
    return interest

def predict(relation_net):
    '''
    predict the validations' interest
    '''

    fid_result = codecs.open('predict2.txt','w','utf-8')

    with codecs.open('validation.txt','r','utf-8') as fid:
        for line in fid:
            if line == '\n':
                continue
            line = line.strip()
            fid_result.write(line)
            interest = []
            for v in nx.all_neighbors(relation_net,line):
                if len(relation_net.node[v]['interest']) <= 0:
                    interest.extend(expandInterest(v,relation_net)*relation_net[line][v]['weight'])
                else:
                    interest.extend(relation_net.node[v]['interest']*relation_net[line][v]['weight'])
            interest = Counter(interest).most_common(5)

            for i in range(len(interest)):
                fid_result.write('\t%s' %interest[i][0])

            for i in range(5-len(interest)):
                fid_result.write('\t0')
            fid_result.write('\n')

    fid_result.close()

if __name__ == "__main__":
    print('begin to create relation network ...')
    relation_net = createRelationNetwork()
    print('succeed creating network!\n')
    print('begin to predict ...')
    predict(relation_net)

max_stars_repo_name: luolidong/aca

max_issues_repo_name: luolidong/aca

max_forks_repo_name: luolidong/aca

reason: import networkx


text:
import numpy as np
import argparse
import cv2

ap = argparse.ArgumentParser()
ap.add_argument("-i", "--image", type=str, required=True,
                help="path to input black and white image")
ap.add_argument("-p", "--prototxt", type=str, required=True,
                help="path to Caffe prototxt file")
ap.add_argument("-m", "--model", type=str, required=True,
                help="path to Caffe pre-trained model")
ap.add_argument("-c", "--points", type=str, required=True,
                help="path to cluster center points")
args = vars(ap.parse_args())

print("[INFO] loading model...")
net = cv2.dnn.readNetFromCaffe(args["prototxt"], args["model"])
pts = np.load(args["points"])
class8 = net.getLayerId("class8_ab")
conv8 = net.getLayerId("conv8_313_rh")
pts = pts.transpose().reshape(2, 313, 1, 1)
net.getLayer(class8).blobs = [pts.astype("float32")]
net.getLayer(conv8).blobs = [np.full([1, 313], 2.606, dtype="float32")]

image = cv2.imread(args["image"])
scaled = image.astype("float32") / 255.0
lab = cv2.cvtColor(scaled, cv2.COLOR_BGR2LAB)

resized = cv2.resize(lab, (224, 224))
L = cv2.split(resized)[0]
L -= 50

'print("[INFO] colorizing image...")'
net.setInput(cv2.dnn.blobFromImage(L))
ab = net.forward()[0, :, :, :].transpose((1, 2, 0))

ab = cv2.resize(ab, (image.shape[1], image.shape[0]))

L = cv2.split(lab)[0]
colorized = np.concatenate((L[:, :, np.newaxis], ab), axis=2)

colorized = cv2.cvtColor(colorized, cv2.COLOR_LAB2BGR)
colorized = np.clip(colorized, 0, 1)

colorized = (255 * colorized).astype("uint8")
cv2.imshow("Original", image)
cv2.imshow("Colorized", colorized)
cv2.waitKey(0)

max_stars_repo_name: prashanth-exe/ReColour


text:
import glob, os, cv2
import numpy as np
import json,tqdm
import copy
import numpy as np                                 # (pip install numpy)
from skimage import measure                        # (pip install scikit-image)
from shapely.geometry import Polygon, MultiPolygon # (pip install Shapely)

def _get_extreme_points(pts):
  l, t = min(pts[:, 0]), min(pts[:, 1])
  r, b = max(pts[:, 0]), max(pts[:, 1])
  # 3 degrees
  thresh = 0.02
  w = r - l + 1
  h = b - t + 1

  pts = np.concatenate([pts[-1:], pts, pts[:1]], axis=0)
  t_idx = np.argmin(pts[:, 1])
  t_idxs = [t_idx]
  tmp = t_idx + 1
  while tmp < pts.shape[0] and pts[tmp, 1] - pts[t_idx, 1] <= thresh * h:
    t_idxs.append(tmp)
    tmp += 1
  tmp = t_idx - 1
  while tmp >= 0 and pts[tmp, 1] - pts[t_idx, 1] <= thresh * h:
    t_idxs.append(tmp)
    tmp -= 1
  tt = [(max(pts[t_idxs, 0]) + min(pts[t_idxs, 0])) // 2, t]

  b_idx = np.argmax(pts[:, 1])
  b_idxs = [b_idx]
  tmp = b_idx + 1
  while tmp < pts.shape[0] and pts[b_idx, 1] - pts[tmp, 1] <= thresh * h:
    b_idxs.append(tmp)
    tmp += 1
  tmp = b_idx - 1
  while tmp >= 0 and pts[b_idx, 1] - pts[tmp, 1] <= thresh * h:
    b_idxs.append(tmp)
    tmp -= 1
  bb = [(max(pts[b_idxs, 0]) + min(pts[b_idxs, 0])) // 2, b]

  l_idx = np.argmin(pts[:, 0])
  l_idxs = [l_idx]
  tmp = l_idx + 1
  while tmp < pts.shape[0] and pts[tmp, 0] - pts[l_idx, 0] <= thresh * w:
    l_idxs.append(tmp)
    tmp += 1
  tmp = l_idx - 1
  while tmp >= 0 and pts[tmp, 0] - pts[l_idx, 0] <= thresh * w:
    l_idxs.append(tmp)
    tmp -= 1
  ll = [l, (max(pts[l_idxs, 1]) + min(pts[l_idxs, 1])) // 2]

  r_idx = np.argmax(pts[:, 0])
  r_idxs = [r_idx]
  tmp = r_idx + 1
  while tmp < pts.shape[0] and pts[r_idx, 0] - pts[tmp, 0] <= thresh * w:
    r_idxs.append(tmp)
    tmp += 1
  tmp = r_idx - 1
  while tmp >= 0 and pts[r_idx, 0] - pts[tmp, 0] <= thresh * w:
    r_idxs.append(tmp)
    tmp -= 1
  rr = [r, (max(pts[r_idxs, 1]) + min(pts[r_idxs, 1])) // 2]

  return np.array([tt, ll, bb, rr])

#http://www.immersivelimit.com/tutorials/create-coco-annotations-from-scratch/#create-custom-coco-dataset
def create_sub_mask_annotation(sub_mask, image_id, category_id, annotation_id, is_crowd):
    # Find contours (boundary lines) around each sub-mask
    # Note: there could be multiple contours if the object
    # is partially occluded. (E.g. an elephant behind a tree)
    contours = measure.find_contours(sub_mask, 0.5, positive_orientation='low')

    segmentations = []
    polygons = []
    for contour in contours:
        # Flip from (row, col) representation to (x, y)
        # and subtract the padding pixel
        for i in range(len(contour)):
            row, col = contour[i]
            contour[i] = (col - 1, row - 1)

        # Make a polygon and simplify it
        poly = Polygon(contour)
        # print("[create_sub_mask_annotation] poly", poly, "contour", contour)
        poly = poly.simplify(1.0, preserve_topology=True) #False
        polygons.append(poly)
        segmentation = np.array(poly.exterior.coords).ravel().tolist()
        segmentations.append(segmentation)

    # Combine the polygons to calculate the bounding box and area
    multi_poly = MultiPolygon(polygons)
    x, y, max_x, max_y = multi_poly.bounds
    width = max_x - x
    height = max_y - y
    bbox = (x, y, width, height)
    area = multi_poly.area

    annotation = {
        'segmentation': segmentations,
        'iscrowd': is_crowd,
        'image_id': image_id,
        'category_id': category_id,
        'id': annotation_id,
        'bbox': bbox,
        'area': area
    }

    return annotation

if '__main__' == __name__:
	output_file  = "../data/medical_img/annotations"
	if not os.path.exists(output_file):
		print("NOT EXIST, to create one")
		os.makedirs(output_file)
	output_file_extreme = os.path.join(output_file,"instances_extreme_train2017.json")
	output_file_init = os.path.join(output_file,"instances_train2017.json")
	del output_file

	labels_list_to_check = [4,5,6]#[3,4]
	data_extreme = {"annotations":[],"images":[],"categories": []}
	min_base_index = min(labels_list_to_check) - 1
	# for label in labels_list_to_check:
  #========
	# data_extreme["categories"].append({"id":3-min_base_index,"name":"spin_cord", "supercategory":"spin_cord"})
	# data_extreme["categories"].append({"id":4-min_base_index,"name":"probe_right", "supercategory":"probe"})
	#========
	data_extreme["categories"].append({"id":4-min_base_index,"name":"probe_right", "supercategory":"probe"})
	data_extreme["categories"].append({"id":6-min_base_index,"name":"probe_left", "supercategory":"probe"})
	data_extreme["categories"].append({"id":5-min_base_index,"name":"scissor", "supercategory":"scissor"})
  #========
	mask_list = glob.glob('../../../medical_img/data/test/mask_label/*.png')
	# print("mask_list", mask_list)
	mask_list = [m for m in mask_list if "corrected" not in m]
	mask_list.sort()
	annotation_id = 1
	data_init = copy.deepcopy(data_extreme)
	for i in tqdm.tqdm(range(len(mask_list))):
		mask_name = mask_list[i]
		splits = mask_name.split('/')[-1].split('_')[0]#[:-4].split('_')[0].split('e')
		index = int(''.join([s for s in splits if s.isdigit()]))
		# mask = cv2.imread(mask_name,0)#.astype(np.float64)

		img_origin_directory = "/home/yezheng/medical_img/data/test/img"
		img_origin = cv2.imread(os.path.join(img_origin_directory,
			"{:06d}.jpg".format(index) ) )
		if None is img_origin:
			img_origin_directory = "/Users/yezheng/medical_img/data/test/img"
			img_origin = cv2.imread(os.path.join(img_origin_directory,
			"{:06d}.jpg".format(index) ) )
		data_extreme["images"].append({"file_name":os.path.join(img_origin_directory,
			"{:06d}.jpg".format(index) ),"id":index})
		data_init["images"].append({"file_name":os.path.join(img_origin_directory,
			"{:06d}.jpg".format(index) ),"id":index})

		# mask_corrected = np.zeros(mask.shape)
		for label in labels_list_to_check:

			mask_corrected_name = '../../../medical_img/data/test/mask_label/_corrected{}/correctedFrame{:04d}_ordered_{}.png'.format(label, index, label)

			mask_corrected = cv2.imread(mask_corrected_name,0)
			# print(mask_corrected.shape)
			if None is mask_corrected:
				print("mask_corrected_name NOT EXIST", mask_corrected_name)
				continue
			mask_corrected = (mask_corrected >0).astype(np.int32)
			ann = create_sub_mask_annotation(sub_mask = mask_corrected, image_id = index, category_id = label-min_base_index, annotation_id = annotation_id , is_crowd = 0)
			# print("[gen_medical_extreme_points] ann", ann)
			annotation_id += 1 # automatically increase
			tmp = np.where(mask_corrected > 0)
			pts = np.asarray(tmp).transpose()[:, ::-1].astype(np.int32)
			extreme_points = _get_extreme_points(pts).astype(np.int32)
			data_init['annotations'].append(ann)
			ann['extreme_points'] = extreme_points.copy().tolist()
			data_extreme['annotations'].append(ann)
			del ann

	json.dump(data_extreme, open(output_file_extreme, 'w'))
	json.dump(data_init, open(output_file_init, 'w'))


text:
import numpy as np
import jax
import jax.numpy as jnp
from jax import jit, device_put, vmap
import optax
from sklearn.cluster import KMeans
from tqdm.auto import tqdm
from .linear_operator_base import LinearOperator, Lazy
from .linear_operators import ConcatLazy, I, lazify, densify, LazyJVP
import logging
import matplotlib.pyplot as plt
from functools import reduce
from emlp.utils import export, memory
import torch

from plum import dispatch
import emlp.reps

# TODO: add rep,v = flatten({'Scalar':..., 'Vector':...,}), to_dict(rep,vector) returns {'Scalar':..., 'Vector':...,}
# TODO and simpler rep = flatten({Scalar:2,Vector:10,...}),
# Do we even want + operator to implement non canonical orderings?

__all__ = ["V", "Vector", "Scalar"]

@export
class Rep(object):
    r""" The base Representation class. Representation objects formalize the vector space V
        on which the group acts, the group representation matrix ρ(g), and the Lie Algebra
        representation dρ(A) in a single object. Representations act as types for vectors coming
        from V. These types can be manipulated and transformed with the built in operators
        ⊕,⊗,dual, as well as incorporating custom representations. Rep objects should
        be immutable.

        At minimum, new representations need to implement ``rho``, ``__str__``."""

    is_permutation = False

    def __init__(self):
        self.equivariant_projector = memory.cache(self.equivariant_projector)
        self.equivariant_basis = memory.cache(self.equivariant_basis)

    def rho(self, M):
        """ Group representation of the matrix M of shape (d,d)"""
        raise NotImplementedError

    def drho(self, A):
        """ Lie Algebra representation of the matrix A of shape (d,d)"""
        In = jnp.eye(A.shape[0])
        return LazyJVP(self.rho, In, A)

    def __call__(self, G):
        """ Instantiate (non concrete) representation with a given symmetry group"""
        raise NotImplementedError

    def __str__(self):
        raise NotImplementedError

    # TODO: separate __repr__ and __str__?
    def __repr__(self):
        return str(self)

    def __eq__(self, other):
        if type(self) != type(other): return False
        d1 = tuple(
            [(k, v) for k, v in self.__dict__.items() if (k not in ['_size', 'is_permutation', 'is_orthogonal'])])
        d2 = tuple(
            [(k, v) for k, v in other.__dict__.items() if (k not in ['_size', 'is_permutation', 'is_orthogonal'])])
        return d1 == d2

    def __hash__(self):
        state_dict = self.__dict__.copy()

        # Turn lists into tuples for hashing purposes
        for k in list(state_dict.keys()):
            if isinstance(state_dict[k], list):
                state_dict[k] = tuple(state_dict[k])

        d1 = tuple([(k, v) for k, v in state_dict.items() if (k not in ['_size', 'is_permutation', 'is_orthogonal'])])
        return hash((type(self), d1))

    def size(self):
        """ Dimension dim(V) of the representation """
        if hasattr(self, '_size'):
            return self._size
        elif self.concrete and hasattr(self, "G"):
            self._size = self.rho(self.G.sample()).shape[-1]
            return self._size
        else:
            raise NotImplementedError

    def canonicalize(self):
        """ An optional method to convert the representation into a canonical form
            in order to reuse equivalent solutions in the solver. Should return
            both the canonically ordered representation, along with a permutation
            which can be applied to vectors of the current representation to achieve
            that ordering. """
        return self, np.arange(self.size())  # return canonicalized rep

    def rho_dense(self, M):
        """ A convenience function which returns rho(M) as a dense matrix."""
        return densify(self.rho(M))

    def drho_dense(self, A):
        """ A convenience function which returns drho(A) as a dense matrix."""
        return densify(self.drho(A))

    def constraint_matrix(self):
        """ Constructs the equivariance constrant matrix (lazily) by concatenating
        the constraints (ρ(hᵢ)-I) for i=1,...M and dρ(Aₖ) for k=1,..,D from the generators
        of the symmetry group. """
        n = self.size()
        constraints = []
        constraints.extend([lazify(self.rho(h)) - I(n) for h in self.G.discrete_generators])
        constraints.extend([lazify(self.drho(A)) for A in self.G.lie_algebra])
        return ConcatLazy(constraints) if constraints else lazify(jnp.zeros((1, n)))

    solcache = {}

    def equivariant_basis(self):
        """ Computes the equivariant solution basis for the given representation of size N.
            Canonicalizes problems and caches solutions for reuse. Output [Q (N,r)] """
        if self == Scalar: return jnp.ones((1, 1))
        canon_rep, perm = self.canonicalize()
        invperm = np.argsort(perm)
        if canon_rep not in self.solcache:
            logging.info(f"{canon_rep} cache miss")
            logging.info(f"Solving basis for {self}" + (f", for G={self.G}" if hasattr(self, "G") else ""))
            # if isinstance(group,Trivial): return np.eye(size(rank,group.d))
            C_lazy = canon_rep.constraint_matrix()
            if C_lazy.shape[0] * C_lazy.shape[1] > 3e7:  # Too large to use SVD
                result = krylov_constraint_solve(C_lazy)
            else:
                C_dense = C_lazy.to_dense()
                result = orthogonal_complement(C_dense)
            self.solcache[canon_rep] = result
        return self.solcache[canon_rep][invperm]

    def equivariant_projector(self):
        """ Computes the (lazy) projection matrix P=QQᵀ that projects to the equivariant basis."""
        Q = self.equivariant_basis()
        Q_lazy = lazify(Q)
        P = Q_lazy @ Q_lazy.H
        # print(f'For rep {self.__str__()}')
        # print('P = ', P.to_dense().shape)
        # print('Q = ', Q.shape)
        return P

    @property
    def concrete(self):
        return hasattr(self, "G") and self.G is not None
        # if hasattr(self,"_concrete"): return self._concrete
        # else:
        #     return hasattr(self,"G") and self.G is not None

    def __add__(self, other):
        """ Direct sum (⊕) of representations. """
        if isinstance(other, int):
            if other == 0:
                return self
            else:
                return self + other * Scalar
        elif emlp.reps.product_sum_reps.both_concrete(self, other):
            return emlp.reps.product_sum_reps.SumRep(self, other)
        else:
            return emlp.reps.product_sum_reps.DeferredSumRep(self, other)

    def __radd__(self, other):
        if isinstance(other, int):
            if other == 0:
                return self
            else:
                return other * Scalar + self
        else:
            return NotImplemented

    def __mul__(self, other):
        """ Tensor sum (⊗) of representations. """
        return mul_reps(self, other)

    def __rmul__(self, other):
        return mul_reps(other, self)

    def __pow__(self, other):
        """ Iterated tensor product. """
        assert isinstance(other, int), f"Power only supported for integers, not {type(other)}"
        assert other >= 0, f"Negative powers {other} not supported"
        return reduce(lambda a, b: a * b, other * [self], Scalar)

    def __rshift__(self, other):
        """ Linear maps from self -> other """
        return other * self.T

    def __lshift__(self, other):
        """ Linear maps from other -> self """
        return self * other.T

    def __lt__(self, other):
        """ less than defined to disambiguate ordering multiple different representations.
            Canonical ordering is determined first by Group, then by size, then by hash"""
        if other == Scalar: return False
        try:
            if self.G < other.G: return True
            if self.G > other.G: return False
        except (AttributeError, TypeError):
            pass
        if self.size() < other.size(): return True
        if self.size() > other.size(): return False
        return hash(self) < hash(other)  # For sorting purposes only

    def __mod__(self, other):  # Wreath product
        """ Wreath product of representations (Not yet implemented)"""
        raise NotImplementedError

    @property
    def T(self):
        """ Dual representation V*, rho*, drho*."""
        if hasattr(self, "G") and (self.G is not None) and self.G.is_orthogonal: return self
        return Dual(self)

@dispatch
def mul_reps(ra, rb: int):
    if rb == 1: return ra
    if rb == 0: return 0
    if (not hasattr(ra, 'concrete')) or ra.concrete:
        return emlp.reps.product_sum_reps.SumRep(*(rb * [ra]))
    else:
        return emlp.reps.product_sum_reps.DeferredSumRep(*(rb * [ra]))

@dispatch
def mul_reps(ra: int, rb):
    return mul_reps(rb, ra)

# Continued with non int cases in product_sum_reps.py

# A possible
class ScalarRep(Rep):
    def __init__(self, G=None):
        super().__init__()
        self.G = G
        self.is_permutation = True

    def __call__(self, G):
        self.G = G
        return self

    def size(self):
        return 1

    def __repr__(self):
        return str(self)  # f"T{self.rank+(self.G,)}"

    def __str__(self):
        return "V⁰"

    @property
    def T(self):
        return self

    def rho(self, M):
        return jnp.eye(1)

    def drho(self, M):
        return 0 * jnp.eye(1)

    def __hash__(self):
        return 0

    def __eq__(self, other):
        return isinstance(other, ScalarRep)

    def __mul__(self, other):
        if isinstance(other, int): return super().__mul__(other)
        return other

    def __rmul__(self, other):
        if isinstance(other, int): return super().__rmul__(other)
        return other

    @property
    def concrete(self):
        return True

class Base(Rep):
    """ Base representation V of a group."""

    def __init__(self, G=None):
        super().__init__()
        self.G = G
        if G is not None: self.is_permutation = G.is_permutation

    def __call__(self, G):
        return self.__class__(G)

    def rho(self, M):
        if hasattr(self, 'G') and isinstance(M, dict): M = M[self.G]
        return M

    def drho(self, A):
        if hasattr(self, 'G') and isinstance(A, dict): A = A[self.G]
        return A

    def size(self):
        assert self.G is not None, f"must know G to find size for rep={self}"
        return self.G.d

    def __repr__(self):
        return str(self)  # f"T{self.rank+(self.G,)}"

    def __str__(self):
        return "V"  # +(f"_{self.G}" if self.G is not None else "")

    def __hash__(self):
        return hash((type(self), self.G))

    def __eq__(self, other):
        return type(other) == type(self) and self.G == other.G

    def __lt__(self, other):
        if isinstance(other, Dual): return True
        return super().__lt__(other)
    # @property
    # def T(self):
    #     return Dual(self.G)

class Dual(Rep):
    def __init__(self, rep):
        super().__init__()
        self.rep = rep
        self.G = rep.G
        if hasattr(rep, "is_permutation"): self.is_permutation = rep.is_permutation

    def __call__(self, G):
        return self.rep(G).T

    def rho(self, M):
        rho = self.rep.rho(M)
        rhoinvT = rho.invT() if isinstance(rho, LinearOperator) else jnp.linalg.inv(rho).T
        return rhoinvT

    def drho(self, A):
        return -self.rep.drho(A).T

    def __str__(self):
        return str(self.rep) + "*"

    def __repr__(self):
        return str(self)

    @property
    def T(self):
        return self.rep

    def __eq__(self, other):
        return type(other) == type(self) and self.rep == other.rep

    def __hash__(self):
        return hash((type(self), self.rep))

    def __lt__(self, other):
        if other == self.rep: return False
        return super().__lt__(other)

    def size(self):
        return self.rep.size()

V = Vector = Base()  #: Alias V or Vector for an instance of the Base representation of a group

Scalar = ScalarRep()  #: An instance of the Scalar representation, equivalent to V**0

@export
def T(p, q=0, G=None):
    """ A convenience function for creating rank (p,q) tensors."""
    return (V ** p * V.T ** q)(G)

def orthogonal_complement(proj):
    """ Computes the orthogonal complement to a given matrix proj"""
    U, S, VH = jnp.linalg.svd(proj, full_matrices=True)
    rank = (S > 1e-5).sum()
    return VH[rank:].conj().T

def krylov_constraint_solve(C, tol=1e-5):
    """ Computes the solution basis Q for the linear constraint CQ=0  and QᵀQ=I
        up to specified tolerance with C expressed as a LinearOperator. """
    r = 5
    if C.shape[0] * r * 2 > 2e9: raise Exception(f"Solns for contraints {C.shape} too large to fit in memory")
    found_rank = 5
    while found_rank == r:
        r *= 2  # Iterative doubling of rank until large enough to include the full solution space
        if C.shape[0] * r > 2e9:
            logging.error(f"Hit memory limits, switching to sample equivariant subspace of size {found_rank}")
            break
        Q = krylov_constraint_solve_upto_r(C, r, tol)
        found_rank = Q.shape[-1]
    return Q

def krylov_constraint_solve_upto_r(C, r, tol=1e-5, lr=1e-2):  # ,W0=None):
    """ Iterative routine to compute the solution basis to the constraint CQ=0 and QᵀQ=I
        up to the rank r, with given tolerance. Uses gradient descent (+ momentum) on the
        objective |CQ|^2, which provably converges at an exponential rate."""
    W = np.random.randn(C.shape[-1], r) / np.sqrt(C.shape[-1])  # if W0 is None else W0
    W = device_put(W)
    opt_init, opt_update = optax.sgd(lr, .9)
    opt_state = opt_init(W)  # init stats

    def loss(W):
        return (jnp.absolute(C @ W) ** 2).sum() / 2  # added absolute for complex support

    loss_and_grad = jit(jax.value_and_grad(loss))
    # setup progress bar
    pbar = tqdm(total=100, desc=f'Krylov Solving for Equivariant Subspace r<={r}',
                bar_format="{l_bar}{bar}| {n:.3g}/{total_fmt} [{elapsed}<{remaining}, {rate_fmt}{postfix}]")
    prog_val = 0
    lstart, _ = loss_and_grad(W)

    for i in range(20000):

        lossval, grad = loss_and_grad(W)
        updates, opt_state = opt_update(grad, opt_state, W)
        W = optax.apply_updates(W, updates)
        # update progress bar
        progress = max(100 * np.log(lossval / lstart) / np.log(tol ** 2 / lstart) - prog_val, 0)
        progress = min(100 - prog_val, progress)
        if progress > 0:
            prog_val += progress
            pbar.update(progress)

        if jnp.sqrt(lossval) < tol:  # check convergence condition
            pbar.close()
            break  # has converged
        if lossval > 2e3 and i > 100:  # Solve diverged due to too high learning rate
            logging.warning(f"Constraint solving diverged, trying lower learning rate {lr / 3:.2e}")
            if lr < 1e-4: raise ConvergenceError(f"Failed to converge even with smaller learning rate {lr:.2e}")
            return krylov_constraint_solve_upto_r(C, r, tol, lr=lr / 3)
    else:
        raise ConvergenceError("Failed to converge.")
    # Orthogonalize solution at the end
    U, S, VT = np.linalg.svd(np.array(W), full_matrices=False)
    # Would like to do economy SVD here (to not have the unecessary O(n^2) memory cost)
    # but this is not supported in numpy (or Jax) unfortunately.
    rank = (S > 10 * tol).sum()
    Q = device_put(U[:, :rank])
    # final_L
    final_L = loss_and_grad(Q)[0]
    if final_L > tol: logging.warning(f"Normalized basis has too high error {final_L:.2e} for tol {tol:.2e}")
    scutoff = (S[rank] if r > rank else 0)
    assert rank == 0 or scutoff < S[rank - 1] / 100, f"Singular value gap too small: {S[rank - 1]:.2e} \
        above cutoff {scutoff:.2e} below cutoff. Final L {final_L:.2e}, earlier {S[rank - 5:rank]}"
    # logging.debug(f"found Rank {r}, above cutoff {S[rank-1]:.3e} after {S[rank] if r>rank else np.inf:.3e}. Loss {final_L:.1e}")
    return Q

class ConvergenceError(Exception): pass

@export
def sparsify_basis(Q, lr=1e-2):  # (n,r)
    """ Convenience function to attempt to sparsify a given basis by applying an orthogonal transformation
        W, Q' = QW where Q' has only 1s, 0s and -1s. Notably this method does not have the same convergence
        gauruntees of krylov_constraint_solve and can fail (even silently). Intended to be used only for
        visualization purposes, use at your own risk. """
    W = np.random.randn(Q.shape[-1], Q.shape[-1])
    W, _ = np.linalg.qr(W)
    W = device_put(W.astype(jnp.float32))
    opt_init, opt_update = optax.adam(lr)  # optax.sgd(1e2,.9)#optax.adam(lr)#optax.sgd(3e-3,.9)#optax.adam(lr)
    opt_update = jit(opt_update)
    opt_state = opt_init(W)  # init stats

    def loss(W):
        return jnp.abs(Q @ W.T).mean() + .1 * (jnp.abs(W.T @ W - jnp.eye(W.shape[0]))).mean() + .01 * \
               jax.numpy.linalg.slogdet(W)[1] ** 2

    loss_and_grad = jit(jax.value_and_grad(loss))

    for i in tqdm(range(3000), desc=f'sparsifying basis'):
        lossval, grad = loss_and_grad(W)
        updates, opt_state = opt_update(grad, opt_state, W)
        W = optax.apply_updates(W, updates)
        # W,_ = np.linalg.qr(W)
        if lossval > 1e2 and i > 100:  # Solve diverged due to too high learning rate
            logging.warning(f"basis sparsification diverged, trying lower learning rate {lr / 3:.2e}")
            return sparsify_basis(Q, lr=lr / 3)
    Q = np.copy(Q @ W.T)
    Q[np.abs(Q) < 1e-2] = 0
    Q[np.abs(Q) > 1e-2] /= np.abs(Q[np.abs(Q) > 1e-2])
    A = Q @ (1 + np.arange(Q.shape[-1]))
    if len(np.unique(np.abs(A))) != Q.shape[-1] + 1 and len(np.unique(np.abs(A))) != Q.shape[-1]:
        logging.error(f"Basis elems did not separate: found only {len(np.unique(np.abs(A)))}/{Q.shape[-1]}")
        # raise ConvergenceError(f"Basis elems did not separate: found only {len(np.unique(A))}/{Q.shape[-1]}")
    return Q

class TorchLazyProjection:
    def __init__(self, W_multiplicities, x_multiplicities, nelems, reduced_indices_dict, inv_perm, mat_shape):
        self.W_multiplicities = W_multiplicities
        self.x_multiplicities = x_multiplicities
        self.reduced_indices_dict = reduced_indices_dict
        self.inv_perm = inv_perm
        self.mat_shape = mat_shape

    def __call__(self, params, x):
        bshape = x.shape[:-1]
        x = x.reshape(-1, x.shape[-1])
        bs = x.shape[0]
        i = 0
        Ws = []
        for rep, W_mult in self.W_multiplicities.items():
            if rep not in self.x_multiplicities:
                Ws.append(torch.zeros((bs, W_mult * rep.size()), device=x.device))
                continue
            x_mult = self.x_multiplicities[rep]
            n = nelems(x_mult, rep)
            i_end = i + W_mult * n
            bids = self.reduced_indices_dict[rep]
            bilinear_params = params[i:i_end].reshape(W_mult, n)  # bs,nK-> (nK,bs)
            i = i_end  # (bs,W_mult,d^r) = (W_mult,n)@(n,d^r,bs)
            bilinear_elems = bilinear_params @ x[..., bids].T.reshape(n, rep.size() * bs)
            bilinear_elems = bilinear_elems.reshape(W_mult * rep.size(), bs).T
            Ws.append(bilinear_elems)
        Ws = torch.cat(Ws, dim=-1)  # concatenate over rep axis
        return Ws[..., self.inv_perm].reshape(*bshape, *self.mat_shape)  # reorder to original rank ordering

def nelems(nx, rep):
    return min(nx, rep.size())

# @partial(jit,static_argnums=(0,1))
@export
def torch_bilinear_weights(out_rep, in_rep):
    # TODO: replace lazy_projection function with LazyDirectSum LinearOperator
    W_rep, W_perm = (in_rep >> out_rep).canonicalize()
    # possible bug when in_rep and out_rep are both non sumreps? #TODO: investigate
    inv_perm = np.argsort(W_perm)
    mat_shape = out_rep.size(), in_rep.size()
    x_rep = in_rep
    W_multiplicities = W_rep.reps
    x_multiplicities = x_rep.reps
    x_multiplicities = {rep: n for rep, n in x_multiplicities.items() if rep != Scalar}
    active_dims = sum([W_multiplicities.get(rep, 0) * nelems(n, rep) for rep, n in x_multiplicities.items()])
    reduced_indices_dict = {rep: ids[np.random.choice(len(ids), nelems(len(ids), rep))].reshape(-1) \
                            for rep, ids in x_rep.as_dict(np.arange(x_rep.size())).items()}

    # Apply the projections for each rank, concatenate, and permute back to orig rank order
    # @torch.jit.script
    # def lazy_projection(params, x):  # (r,), (*c) #TODO: find out why backwards of this function is so slow
    #     bshape = x.shape[:-1]
    #     x = x.reshape(-1, x.shape[-1])
    #     bs = x.shape[0]
    #     i = 0
    #     Ws = []
    #     for rep, W_mult in W_multiplicities.items():
    #         if rep not in x_multiplicities:
    #             Ws.append(torch.zeros((bs, W_mult * rep.size()), device=x.device))
    #             continue
    #         x_mult = x_multiplicities[rep]
    #         n = nelems(x_mult, rep)
    #         i_end = i + W_mult * n
    #         bids = reduced_indices_dict[rep]
    #         bilinear_params = params[i:i_end].reshape(W_mult, n)  # bs,nK-> (nK,bs)
    #         i = i_end  # (bs,W_mult,d^r) = (W_mult,n)@(n,d^r,bs)
    #         bilinear_elems = bilinear_params @ x[..., bids].T.reshape(n, rep.size() * bs)
    #         bilinear_elems = bilinear_elems.reshape(W_mult * rep.size(), bs).T
    #         Ws.append(bilinear_elems)
    #     Ws = torch.cat(Ws, dim=-1)  # concatenate over rep axis
    #     return Ws[..., inv_perm].reshape(*bshape, *mat_shape)  # reorder to original rank ordering
    lazy_projection = TorchLazyProjection(W_multiplicities, x_multiplicities, nelems, reduced_indices_dict, inv_perm, mat_shape)

    return active_dims, lazy_projection

@export
def bilinear_weights(out_rep, in_rep):
    # TODO: replace lazy_projection function with LazyDirectSum LinearOperator
    W_rep, W_perm = (in_rep >> out_rep).canonicalize()
    # possible bug when in_rep and out_rep are both non sumreps? #TODO: investigate
    inv_perm = np.argsort(W_perm)
    mat_shape = out_rep.size(), in_rep.size()
    x_rep = in_rep
    W_multiplicities = W_rep.reps
    x_multiplicities = x_rep.reps
    x_multiplicities = {rep: n for rep, n in x_multiplicities.items() if rep != Scalar}
    nelems = lambda nx, rep: min(nx, rep.size())
    active_dims = sum([W_multiplicities.get(rep, 0) * nelems(n, rep) for rep, n in x_multiplicities.items()])
    reduced_indices_dict = {rep: ids[np.random.choice(len(ids), nelems(len(ids), rep))].reshape(-1) \
                            for rep, ids in x_rep.as_dict(np.arange(x_rep.size())).items()}

    # Apply the projections for each rank, concatenate, and permute back to orig rank order
    @jit
    def lazy_projection(params, x):  # (r,), (*c) #TODO: find out why backwards of this function is so slow
        bshape = x.shape[:-1]
        x = x.reshape(-1, x.shape[-1])
        bs = x.shape[0]
        i = 0
        Ws = []
        for rep, W_mult in W_multiplicities.items():
            if rep not in x_multiplicities:
                Ws.append(jnp.zeros((bs, W_mult * rep.size())))
                continue
            x_mult = x_multiplicities[rep]
            n = nelems(x_mult, rep)
            i_end = i + W_mult * n
            bids = reduced_indices_dict[rep]
            bilinear_params = params[i:i_end].reshape(W_mult, n)  # bs,nK-> (nK,bs)
            i = i_end  # (bs,W_mult,d^r) = (W_mult,n)@(n,d^r,bs)
            bilinear_elems = bilinear_params @ x[..., bids].T.reshape(n, rep.size() * bs)
            bilinear_elems = bilinear_elems.reshape(W_mult * rep.size(), bs).T
            Ws.append(bilinear_elems)
        Ws = jnp.concatenate(Ws, axis=-1)  # concatenate over rep axis
        return Ws[..., inv_perm].reshape(*bshape, *mat_shape)  # reorder to original rank ordering

    return active_dims, lazy_projection

# @jit
# def mul_part(bparams,x,bids):
#     b = prod(x.shape[:-1])
#     return (bparams@x[...,bids].T.reshape(bparams.shape[-1],-1)).reshape(-1,b).T

@export
def vis(repin, repout, cluster=True):
    """ A function to visualize the basis of equivariant maps repin>>repout
        as an image. Only use cluster=True if you know Pv will only have
        r distinct values (true for G<S(n) but not true for many continuous groups)."""
    rep = (repin >> repout)
    P = rep.equivariant_projector()  # compute the equivariant basis
    Q = rep.equivariant_basis()
    v = np.random.randn(P.shape[1])  # sample random vector
    v = np.round(P @ v, decimals=4)  # project onto equivariant subspace (and round)
    if cluster:  # cluster nearby values for better color separation in plot
        v = KMeans(n_clusters=Q.shape[-1]).fit(v.reshape(-1, 1)).labels_
    plt.imshow(v.reshape(repout.size(), repin.size()))
    plt.axis('off')

def scale_adjusted_rel_error(t1, t2, g):
    error = jnp.sqrt(jnp.mean(jnp.abs(t1 - t2) ** 2))
    tscale = jnp.sqrt(jnp.mean(jnp.abs(t1) ** 2)) + jnp.sqrt(jnp.mean(jnp.abs(t2) ** 2))
    gscale = jnp.sqrt(jnp.mean(jnp.abs(g - jnp.eye(g.shape[-1])) ** 2))
    scale = jnp.maximum(tscale, gscale)
    return error / jnp.maximum(scale, 1e-7)

@export
def equivariance_error(W, repin, repout, G):
    """ Computes the equivariance relative error rel_err(Wρ₁(g),ρ₂(g)W)
        of the matrix W (dim(repout),dim(repin)) [or basis Q: (dim(repout)xdim(repin), r)]
        according to the input and output representations and group G. """
    W = W.reshape(repout.size(), repin.size(), -1).transpose((2, 0, 1))[None]

    # Sample 5 group elements and verify the equivariance for each
    gs = G.samples(5)
    ring = vmap(repin.rho_dense)(gs)[:, None]
    routg = vmap(repout.rho_dense)(gs)[:, None]
    equiv_err = scale_adjusted_rel_error(W @ ring, routg @ W, gs)
    return equiv_err

import emlp.groups  # Why is this necessary to avoid circular import?

max_stars_repo_name: marcelroed/equivariant-MLP

max_issues_repo_name: marcelroed/equivariant-MLP

max_forks_repo_name: marcelroed/equivariant-MLP

reason: import numpy,import jax,from jax


text:
import os
import argparse
import sys

sys.path.append('../')
from load_paths import load_box_paths
from processing_helpers import *
import pandas as pd
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np

mpl.rcParams['pdf.fonttype'] = 42
sns.set_style('whitegrid', {'axes.linewidth': 0.5})

def parse_args():
    description = "Simulation run for modeling Covid-19"
    parser = argparse.ArgumentParser(description=description)

    parser.add_argument(
        "-s",
        "--stem",
        type=str,
        help="Name of simulation experiment"
    )
    parser.add_argument(
        "-loc",
        "--Location",
        type=str,
        help="Local or NUCLUSTER",
        default="Local"
    )
    parser.add_argument(
        "-perc",
        "--overflow_threshold_percents",
        type=float,
        nargs='+',
        help="Calculate probability for specified percent of capacity limit",
        default=99
    )
    return parser.parse_args()

def get_latest_filedate(file_path=os.path.join(datapath, 'covid_IDPH', 'Corona virus reports',
                                               'hospital_capacity_thresholds'), extraThresholds=False):
    files = os.listdir(file_path)
    files = sorted(files, key=len)
    if extraThresholds == False:
        files = [name for name in files if not 'extra_thresholds' in name]
    if extraThresholds == True:
        files = [name for name in files if 'extra_thresholds' in name]

    filedates = [item.replace('capacity_weekday_average_', '') for item in files]
    filedates = [item.replace('.csv', '') for item in filedates]
    latest_filedate = max([int(x) for x in filedates])
    fname = f'capacity_weekday_average_{latest_filedate}.csv'
    if extraThresholds == True:
        fname = f'capacity_weekday_average_{latest_filedate}__extra_thresholds.csv'
    return fname

def get_probs(ems_nr, channels=['total_hosp_census', 'crit_det', 'ventilators'], overflow_threshold_percents=[1, 0.8],
              param=None, save_csv=False, plot=True):
    """Define columns and labels"""
    if ems_nr == 0:
        region_suffix = "_All"
        region_label = 'Illinois'
    else:
        region_suffix = "_EMS-" + str(ems_nr)
        region_label = region_suffix.replace('_EMS-', 'COVID-19 Region ')

    column_list = ['scen_num', 'sample_num', 'time', 'startdate']
    grp_channels = ['date']
    if param is not None:
        column_list = column_list + param
        grp_channels = ['date'] + param

    column_list_t = column_list
    for channel in ['hosp_det', 'crit_det']:
        column_list_t.append(channel + region_suffix)

    """ Load dataframes"""
    df = load_sim_data(exp_name, region_suffix=region_suffix, column_list=column_list, add_incidence=False)
    df['total_hosp_census'] = df['hosp_det'] + df['crit_det']
    df['ventilators'] = get_vents(df['crit_det'])

    capacity_df = load_capacity(ems_nr)
    len(df['scen_num'].unique())
    df['N_scen_num'] = df.groupby(grp_channels)['scen_num'].transform('count')

    df_all = pd.DataFrame()
    for channel in channels:
        if channel == "crit_det": channel_label = 'icu_availforcovid'
        if channel == "hosp_det": channel_label = 'hb_availforcovid'
        if channel == "total_hosp_census": channel_label = 'hb_availforcovid'
        if channel == "ventilators": channel_label = 'vent_availforcovid'

        for overflow_threshold_percent in overflow_threshold_percents:
            thresh = capacity_df[f'{channel}'] * overflow_threshold_percent

            mdf = df.copy()
            mdf.loc[df[f'{channel}'] < thresh, 'above_yn'] = 0
            mdf.loc[df[f'{channel}'] >= thresh, 'above_yn'] = 1

            mdf = mdf.groupby(grp_channels)['above_yn'].agg(['sum', 'count', 'nunique']).rename_axis(None, axis=0)
            mdf = mdf.reset_index()
            mdf['prob'] = mdf['sum'] / mdf['count']

            mdf = mdf.rename(columns={'sum': 'n_above', 'count': 'N_scen_num', 'index': 'date'})
            mdf['overflow_threshold_percent'] = overflow_threshold_percent
            mdf['capacity_channel'] = channel_label
            mdf['availforcovid'] = capacity_df[f'{channel}']
            mdf['region'] = ems_nr
            del thresh

            if df_all.empty:
                df_all = mdf
            else:
                df_all = pd.concat([df_all, mdf])
            del mdf
    if plot:
        plot_probs(df=df_all, region_label=region_label)

    if save_csv:
        filename = f'overflow_probabilities_over_time_region_{ems_nr}.csv'
        df_all.to_csv(os.path.join(sim_output_path, filename), index=False, date_format='%Y-%m-%d')
    return df_all

def plot_probs(df, region_label):
    fig = plt.figure(figsize=(12, 4))
    fig.suptitle(region_label, y=0.97, fontsize=14)
    fig.subplots_adjust(right=0.98, wspace=0.2, left=0.05, hspace=0.4, top=0.84, bottom=0.13)
    palette = sns.color_palette('Set1', 12)
    axes = [fig.add_subplot(1, 3, x + 1) for x in range(3)]

    linestyles = ['solid', 'dashed']
    for c, channel in enumerate(df.capacity_channel.unique()):

        mdf = df[df['capacity_channel'] == channel]
        ax = axes[c]
        ax.set_ylim(0, 1)
        ax.set_title(channel)
        ax.set_ylabel(f'Probability of overflow')

        for e, t in enumerate(list(df.overflow_threshold_percent.unique())):
            line_label = f'{channel} ({t * 100})%'
            adf = mdf[mdf['overflow_threshold_percent'] == t]
            ax.plot(adf['date'], adf['prob'], linestyle=linestyles[e], color=palette[c], label=line_label)

        ax.xaxis.set_major_formatter(mdates.DateFormatter('%d\n%b\n%Y'))

    axes[-1].legend()

    plotname = f'overflow_probabilities_{region_label}'
    plt.savefig(os.path.join(plot_path, f'{plotname}.png'))
    plt.savefig(os.path.join(plot_path, 'pdf', f'{plotname}.pdf'))

def write_probs_to_template(df, plot=True):
    fname_capacity = get_latest_filedate()
    civis_template = pd.read_csv(
        os.path.join(datapath, 'covid_IDPH', 'Corona virus reports', 'hospital_capacity_thresholds', fname_capacity))
    civis_template = civis_template.drop_duplicates()
    civis_template['date_window_upper_bound'] = pd.to_datetime(civis_template['date_window_upper_bound'])

    civis_template_all = pd.DataFrame()
    for index, row in civis_template.iterrows():
        upper_limit = row['date_window_upper_bound']
        lower_limit = upper_limit - pd.Timedelta(7, 'days')

        df_sub = df[df['date'].between(lower_limit, upper_limit)]
        df_sub = df_sub[df_sub['region'] == int(row['geography_modeled'].replace("covidregion_", ""))]
        df_sub = df_sub[df_sub['capacity_channel'] == row['resource_type']]
        df_sub = df_sub[df_sub['overflow_threshold_percent'] == row['overflow_threshold_percent']]

        """Take maximum of previous 7 days"""
        civis_template.loc[index, 'percent_of_simulations_that_exceed'] = df_sub['prob'].max()

    if civis_template_all.empty:
        civis_template_all = civis_template
    else:
        civis_template_all = pd.concat([civis_template_all, civis_template])

    """Replace NAs with zero """
    civis_template_all['percent_of_simulations_that_exceed'] = civis_template_all[
        'percent_of_simulations_that_exceed'].fillna(0)
    """Scenario name of simulation - here hardcoded to baseline!!"""
    civis_template_all['scenario_name'] = 'baseline'
    file_str = 'nu_hospitaloverflow_' + str(exp_name[:8]) + '.csv'
    civis_template_all.to_csv(os.path.join(sim_output_path, file_str), index=False)

    if plot:
        plot_probs_from_template(df=civis_template_all)

def plot_probs_from_template(df=None, show_75=True):
    if df is None:
        file_str = 'nu_hospitaloverflow_' + str(exp_name[:8]) + '.csv'
        df = pd.read_csv(os.path.join(sim_output_path, file_str))

    regionlist = df['geography_modeled'].unique()
    df['date_md'] = df['date_window_upper_bound'].dt.strftime('%m-%d\n%Y')
    df['region'] = df['geography_modeled'].str.replace('covidregion_', '')

    fig = plt.figure(figsize=(14, 12))
    fig.suptitle('Overflow probability per week dates by COVID-19 Region', y=0.97, fontsize=14)
    fig.subplots_adjust(right=0.98, wspace=0.4, left=0.05, hspace=0.4, top=0.90, bottom=0.07)
    palette = sns.color_palette('Set1', len(df.resource_type.unique()))
    axes = [fig.add_subplot(4, 3, x + 1) for x in range(len(regionlist))]

    for c, reg_nr in enumerate(regionlist):
        reg_label = reg_nr.replace('covidregion_', 'COVID-19 Region ')
        mdf = df[df['geography_modeled'] == reg_nr]
        ax = axes[c]
        ax.set_ylim(0, 1)
        ax.set_title(reg_label)
        ax.set_ylabel(f'Probability of overflow')

        for e, t in enumerate(list(df.resource_type.unique())):
            adf = mdf[mdf['resource_type'] == t]
            adf1 = adf[adf['overflow_threshold_percent'] == 1]
            adf2 = adf[adf['overflow_threshold_percent'] != 1]
            ax.plot(adf1['date_md'], adf1['percent_of_simulations_that_exceed'], color=palette[e], label=t)
            if show_75:
                ax.plot(adf2['date_md'], adf2['percent_of_simulations_that_exceed'], color=palette[e], label='',
                        alpha=0.5)

    axes[-1].legend()

    plt.savefig(os.path.join(plot_path, 'overflow_probabilities.png'))
    plt.savefig(os.path.join(plot_path, 'pdf', 'overflow_probabilities.pdf'))

if __name__ == '__main__':

    args = parse_args()
    stem = args.stem
    Location = args.Location
    overflow_threshold_percents = args.overflow_threshold_percents
    datapath, projectpath, wdir, exe_dir, git_dir = load_box_paths(Location=Location)

    first_plot_day = pd.Timestamp.today() - pd.Timedelta(14, 'days')
    last_plot_day = pd.Timestamp.today() + pd.Timedelta(90, 'days')

    if overflow_threshold_percents == 99:
        fname_capacity = get_latest_filedate()
        civis_template = pd.read_csv(
            os.path.join(datapath, 'covid_IDPH', 'Corona virus reports', 'hospital_capacity_thresholds',
                         fname_capacity))
        overflow_threshold_percents = civis_template.overflow_threshold_percent.unique()
    print(overflow_threshold_percents)

    exp_names = [x for x in os.listdir(os.path.join(wdir, 'simulation_output')) if stem in x]
    for exp_name in exp_names:
        print(exp_name)
        sim_output_path = os.path.join(wdir, 'simulation_output', exp_name)
        plot_path = os.path.join(sim_output_path, '_plots')

        df_all = pd.DataFrame()
        for ems_nr in range(1, 12):
            print("Start processing region " + str(ems_nr))
            df = get_probs(ems_nr, overflow_threshold_percents=overflow_threshold_percents, save_csv=False, plot=False)
            df = df[df['date'].between(first_plot_day, last_plot_day)]
            if df_all.empty:
                df_all = df
            else:
                df_all = pd.concat([df_all, df])
        df_all.to_csv(os.path.join(sim_output_path, 'overflow_probabilities.csv'), index=False)
        write_probs_to_template(df=df_all)

max_stars_repo_name: msingam/covid-chicago

max_issues_repo_name: msingam/covid-chicago

max_forks_repo_name: msingam/covid-chicago


text:
import matplotlib.pylab as plt
import pandas as pd
import json
import seaborn as sns
import numpy as np

def pbc(positions,lBox):
    dims=len(lBox)
    positions=positions.reshape(positions.size//dims,dims)
    for i in range(dims):
        xpos= positions[:,i] > 0
        positions[xpos]=positions[xpos] - ( (positions[xpos] + lBox[i]/2.)//lBox[i]) * lBox[i]

        xpos= positions[:,i] < 0
        positions[xpos]=positions[xpos] + ( (-(positions[xpos] - lBox[i]/2.))//lBox[i]) * lBox[i]

    return positions

def configurationTable(jSonData):
    jSonConfs=jSonData["configurations"];

    walkers = len(jSonConfs)
    species= len(jSonConfs[0])

    dfs=[]
    for iWalker,walkerData in enumerate(jSonConfs):

        for iSpecies,speciesData in enumerate(walkerData):
            df=pd.DataFrame(speciesData)
            dfs.append( df)
            df["walker"]=iWalker;
            df["species"]=iSpecies

    if ( len(dfs)==0):
        return None
    else:
        data=pd.concat(dfs)
        x= pbc ( np.array(data["x"]).reshape(len(data["x"]),1) ,[10.] )

        data["x"]=x

        return data

def plotConfigurationTable(data):

    fig = plt.figure()
    for species,df in data.groupby("species"):
        print(len(df))
        plt.plot(df["x"],df["x"] * 0  ,marker="o",alpha=0.7,linestyle="none",markeredgewidth=0)
    plt.ylim(-0.1,0.1)

if __name__ == "__main__":

    with open("walkers-Rank0.json") as f:
        data = json.load(f)

    data=configurationTable(data) ;

    data=data[data["walker"]==10 ]
    plotConfigurationTable(data)
    plt.show()


text:
#!/usr/bin/python3

import unittest
from unittest.mock import patch

import numpy as np

from tmc import points

from tmc.utils import load, get_out, patch_helper
module_name='src.almost_meeting_lines'
almost_meeting_lines = load(module_name, 'almost_meeting_lines')
ph = patch_helper(module_name)

@points('p03-08.1')
class AlmostMeetingLines(unittest.TestCase):

    def test_types(self):
        a1=1
        b1=4
        a2=3
        b2=2
        p=(a1,b1,a2,b2)
        system="(a1=%i, b1=%i, a2=%i, b2=%i)" % p
        (x, y), exact = almost_meeting_lines(*p)
        self.assertIsInstance(exact, bool, f"Expected the exact value to be a bool. Got {type(exact)}.")
        self.assertIsInstance(y, float, f"Expected the y coordinate to have type float! Got {type(y)}.")
        self.assertIsInstance(x, float, f"Expected the x coordinate to have type float! Got {type(x)}.")

    def test_first(self):
        a1=1
        b1=4
        a2=3
        b2=2
        p=(a1,b1,a2,b2)
        system="(a1=%i, b1=%i, a2=%i, b2=%i)" % p
        (x, y), exact = almost_meeting_lines(*p)
        self.assertEqual(exact, True, msg="Expected exact solution for system %s!" % system)
        self.assertAlmostEqual(y, a1*x + b1,
                               msg="Solution %f does not satisfy the first equation of system %s!"%(y,system))
        self.assertAlmostEqual(y, a2*x + b2,
                               msg="Solution %f does not satisfy the second equation of system %s!"%(y,system))

    def test_underdetermined(self):
        a1=1
        b1=4
        p=(a1,b1,a1,b1)
        system="(a1=%i, b1=%i, a2=%i, b2=%i)" % p
        (x, y), exact = almost_meeting_lines(*p)
        self.assertEqual(exact, False,
                         msg="Did not expect exact solution for underdetermined system %s!" % system)

    def test_inconsistent(self):
        a1=1
        b1=4
        p=(a1,b1,a1,b1+1)
        system="(a1=%i, b1=%i, a2=%i, b2=%i)" % p
        (x, y), exact = almost_meeting_lines(*p)
        self.assertEqual(exact, False, msg="Did not expect exact solution for inconsistent system %s!" % system)

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


text:
from class_gto import *
from scipy import *
from scipy import linalg
from scipy import interpolate
import sys
import time

def extract_sig(Sigma,oms_f_log,oms_f,outputname="SigGW.out"):
	data = zeros((len(oms_f),3),float)
	data[:,0] = oms_f
	Sig = Spline_Complex(Sigma, oms_f_log,oms_f)
	data[:,1] = real(Sig)
	data[:,2] = imag(Sig)
	savetxt(outputname,data)

def log_mesh(om,nom, ntail_,del_nom = 1):
    """Creates logarithmic mesh on Matsubara axis
       Takes first istart points from mesh om and
       the rest of om mesh is replaced by ntail poinst
       redistribued logarithmically.
       Input:
           om      -- original long mesh
           istart  -- first istart points unchanged
           ntail   -- tail replaced by ntail points only
       Output:
           som             -- smaller mesh created from big om mesh
           sSig[nom,nc]    -- Sig on small mesh
       Also computed but not returned:
           ind_om  -- index array which conatins index to
                      kept Matsubara points
    """

    istart = min(nom, len(om))
    ntail = min(ntail_, len(om)-istart)

    ind_om=[]
    alpha = log((len(om)-1.)/istart)/(ntail-1.)
    for i in range(istart)[::del_nom]:
        ind_om.append(i)
    for i in range(ntail):
        t = int(istart*exp(alpha*i)+0.5)
        if (t != ind_om[-1]):
            ind_om.append(t)

    ind_oms_equal = [[0]]
    for it in range(1,len(ind_om)-1):
        istart = int(0.5*(ind_om[it-1]+ind_om[it])+0.51)
        iend = int(0.5*(ind_om[it]+ind_om[it+1])-0.01)
        equal = [i for i in range(istart,iend+1)]
        ind_oms_equal.append(equal)
    istart = int(0.5*(ind_om[-2]+ind_om[-1])+0.51)
    equal = [i for i in range(istart,ind_om[-1]+1)]
    ind_oms_equal.append(equal)

    oms_equal=[]
    for ind in ind_oms_equal:
        oms_equal.append( array([om[i] for i in ind]) )

    return (ind_om,oms_equal)

def Spline_Real(Arr,mesh_log,mesh_reg):
	#Arr = 1D real array
	#Arr = Arr[:10]
	#mesh_log = mesh_log[:10]
	#print Arr, mesh_log

	Re = interpolate.UnivariateSpline(mesh_log,Arr,s=0)
	Arr_large = Re(mesh_reg)
	#Nmesh = len(mesh_log)
	#print mesh_log[Nmesh/2]
	#print Re(0.05)
	#sys.exit()
	return Arr_large

def Spline_Complex(Arr,mesh_log,mesh_reg):
	#Arr = 1D complex array
	Re = interpolate.UnivariateSpline(mesh_log,real(Arr),s=0)
	Im = interpolate.UnivariateSpline(mesh_log,imag(Arr),s=0)
	Arr_large = Re(mesh_reg)+1j*Im(mesh_reg)
	return Arr_large

def	Gw_to_Gtau(Gw ,oms_f_log,oms_f,tau_mesh,beta):
	if len(shape(Gw)) <= 1:
		Gw_large = Spline_Complex( Gw,oms_f_log,oms_f)
		Gtau = zeros((len(tau_mesh)),float)
		corr2 = -0.5*beta
		for itau,tau in enumerate(tau_mesh):
			phase = exp(-1j*oms_f*tau)
			om_inv = 1/(1j*oms_f)
			G_error = dot(phase, Gw_large)
			corr1 = dot(phase,om_inv)
			Gtau[itau] = 2*real(G_error - corr1) + corr2
		Gtau *= 1./beta
		return Gtau

	else:
		t1 = time.clock()
		M = len(Gw[0,0])
		ist, dict_ist = index_double(M)
		N_ist = len(ist)
		Gw_large = zeros((len(oms_f),M,M),complex)

		for i,(i1,i2) in enumerate(ist):
			Gw_large[:,i1,i2] = Spline_Complex( Gw[:,i1,i2],oms_f_log,oms_f)

		# sum_{iw_n 0 to infty} 1/iw exp(-1j*iw*tau)
		phase = exp(-1j*dot(tau_mesh[newaxis].T,oms_f[newaxis]))
		G_error = tensordot(phase,Gw_large,([1],[0]))
		om_inv = tensordot(1/(1j*oms_f)[newaxis].T,\
				identity(M)[newaxis],axes=([1],[0]))
		corr1 = tensordot(phase,om_inv,([1],[0]))
		corr2 = -0.5*beta*identity(M)[newaxis]
		G2 = 2*real(G_error-corr1)+corr2
		Gtau = reshape(G2,(len(G2),M*M))
		Gtau *= 1./beta
		t2 = time.clock()
		#print "time for Gw_to_Gtau = %2.3f" %(t2-t1)

		return Gtau

def cot(x): return 1./tan(x)

def corr(x):
	"""
	correction function (imaginary part)
	when transforming from tau to freq
	"""
	#Corr_im = 0.0 											#no correction
	#Corr_im = - 1j*x/12. 							#first order
	#Corr_im = 1j*(- x/12. -x**3/720*)  #third order
	if x == 0.: Corr_im = 0.0
	else:
		Corr_im = 1j*((0.5*cot(0.5*x))-1/x) #almost exact

	return Corr_im

def Gtau_to_Gw(Gtau, tau_mesh, tau_large, oms_f_log):
	dtau = tau_large[1]-tau_large[0]
	if len(shape(Gtau)) <= 1:
		Gtau_large = Spline_Real(Gtau,tau_mesh, tau_large)
		Gw = zeros((len(oms_f_log)),complex)
		for iw, w in enumerate(oms_f_log):
			x = w*dtau
			phase = exp(1j*w*tau_large)
			Gw_error = dtau*dot(phase, Gtau_large)
			Corr_im = dtau*corr(x)
			Gw[iw] = Gw_error + (Corr_im)
		#Corr_re = .0
		Corr_re = -real(Gw[-1])
		Gw = Gw + Corr_re
		return Gw

	else:
		M = int(sqrt(len(Gtau[0]))+0.1)
		ist, dict_ist = index_double(M)
		Gtau_large = zeros((len(tau_large),M,M),float)
		for i, (i1,i2) in enumerate(ist):
			Gtau_large[:,i1,i2] = Spline_Real(Gtau[:,i],tau_mesh, tau_large)
		Gw = zeros((len(oms_f_log),M,M),complex)
		for iw, w in enumerate(oms_f_log):
			x = w*dtau
			phase = exp(1j*w*tau_large)
			Gw_error = dtau*tensordot(phase, Gtau_large, axes= ([0],[0]))
			Corr_im = dtau*corr(x)
			Gw[iw] = Gw_error + (Corr_im)*identity(M)
		#Corr_re = 0.0
		Corr_re = -real(Gw[-1])
		Gw = Gw + Corr_re
		return Gw

def Gtau_to_Gw2(Gtau, tau_mesh,tau_large, oms_f_log):
	M = int(sqrt(len(Gtau[0]))+0.1)
	ist, dict_ist = index_double(M)
	Gtau_large = zeros((len(tau_large),M,M),float)
	for i, (i1,i2) in enumerate(ist):
		Gtau_large[:,i1,i2] = Spline_Real(Gtau[:,i],tau_mesh, tau_large)
	return Atau_to_Aw(Gtau_large,tau_large,oms_f_log)

def Atau_to_Aw(Atau, tau_mesh, oms_f_log,option=0):
	exp_wt = exp(+1j*dot(oms_f_log[newaxis].T,tau_mesh[newaxis]))
	if not option:
		Aw = FT_matrix(tau_mesh,exp_wt,Atau)
	else:
		dtau = tau_mesh[1]-tau_mesh[0]
		Aw = tensordot(exp_wt,Atau,axes=([1],[0]))*dtau
	return Aw

def Gtau_to_Pw(Gtau, tau_mesh,tau_large,oms_b_log):
	"""
	P(t) = G(t)*G(-t) = -G(t)*G(beta-t)
			 = -G(t)*G(t)[::-1]
	"""

	if len(shape(Gtau)) <= 1:
		Gtau_large = Spline_Real(Gtau,tau_mesh, tau_large)
		Gtau_large_inv = Gtau_large[::-1] #G(-tau)
		Ptau_large = -2.*Gtau_large*Gtau_large_inv
		option = 0 #0: simps, 1: naive integration
		Pw = real(Atau_to_Aw(Ptau_large,tau_large,oms_b_log,option))
		return Pw

	else:
		t1 = time.clock()
		M = int(sqrt(len(Gtau[0]))+0.1)
		ist, dict_ist = index_double(M)

		Gtau_large = zeros((len(tau_large),M*M),float)
		for i, (i1,i2) in enumerate(ist):
			Gtau_large[:,i] = Spline_Real(Gtau[:,i],tau_mesh, tau_large)

		Nist = len(ist)
		Ntau = len(tau_large)
		Ptau_large = zeros((Ntau, Nist, Nist),float)
		Gtau_large_inv = Gtau_large[::-1] #G(-tau)
		for p,(i,l) in enumerate(ist):
			#print p
			for q,(k,j) in enumerate(ist): #note the order kj, not jk
				#P(ijkl) = P(il)(kj)
				#	=  2*G(ki)(tau)G(lj)(-tau)
				# = -2*G(ki)(tau)G(lj)(beta-tau)
				# = -2*G(ki)(tau)G(lj)[::-1]
				p1 = dict_ist[k,i]
				q1 = dict_ist[l,j]
				Ptau_large[:,p,q] = -2.*Gtau_large[:,p1]*Gtau_large_inv[:,q1]
		t2 = time.clock()
		#print "t Gtau_to_Ptau = %2.2f sec" %(t2-t1)
		option = 0 #0: simps, 1: naive integration
		Pw = real(Atau_to_Aw(Ptau_large,tau_large,oms_b_log,option))
		#Pw = Atau_to_Aw(Ptau_large,tau_large,oms_b_log,option)
		t3 = time.clock()
		#print "t Ptau_to_Pw = %2.2f sec" %(t3-t2)
		return Pw

def FT_matrix(tau,phase,A):
	"""
	Fourier transformation:
			F(iw) = int dtau exp(i*iw*tau) F(tau)
	integration done by simpson's rule
	shape(A) = (Ntau,p,q)
	shape(phase) = (Nw,Ntau)
	"""
	if len(tau)!=len(A):
		print "the # of samplings of X(%2d)\
		and Y(%2d) are different" %(len(tau), len(A))
		sys.exit()

	N = len(tau)
	res = 0.0
	#print N
	if N%2 != 0:
		t1 = tau[1:-1][::2]-tau[:-2][::2]
		t2 = tau[2:][::2]-tau[1:-1][::2]
		T  = (t1+t2)/6.
		T1 = (2.-t2/t1)*T
		T2 = (t1+t2)**2/(t1*t2)*T
		T3 = (2.-t1/t2)*T

		A1 = A[:-2][::2]
		A2 = A[1:-1][::2]
		A3 = A[2:][::2]

		TP1 = T1[newaxis]*phase[:,:-2][:,::2] #t_m * exp(iwn*tm)
		TP2 = T2[newaxis]*phase[:,1:-1][:,::2] #t_m * exp(iwn*tm)
		TP3 = T3[newaxis]*phase[:,2:][:,::2] #t_m * exp(iwn*tm)

		res += tensordot(TP1,A[:-2][::2],axes=([1],[0]))
		res += tensordot(TP2,A[1:-1][::2],axes=([1],[0]))
		res += tensordot(TP3,A[2:][::2],axes=([1],[0]))

	else:
		print "Nmesh for Simpsons is even"
		t1 = tau[1:-1][::2]-tau[:-2][::2]
		t2 = tau[2:][::2]-tau[1:-1][::2]
		#print t1
		#print t2
		T  = (t1+t2)/6.
		T1 = (2.-t2/t1)*T
		T2 = (t1+t2)**2/(t1*t2)*T
		T3 = (2.-t1/t2)*T

		A1 = A[:-2][::2]
		A2 = A[1:-1][::2]
		A3 = A[2:][::2]

		TP1 = T1[newaxis]*phase[:,:-2][:,::2] #t_m * exp(iwn*tm)
		TP2 = T2[newaxis]*phase[:,1:-1][:,::2] #t_m * exp(iwn*tm)
		TP3 = T3[newaxis]*phase[:,2:][:,::2] #t_m * exp(iwn*tm)

		res += tensordot(TP1,A[:-2][::2],axes=([1],[0]))
		res += tensordot(TP2,A[1:-1][::2],axes=([1],[0]))
		res += tensordot(TP3,A[2:][::2],axes=([1],[0]))

		#res += 0.5*(tau[-1]-tau[-2])*(P[-1]+P[-2])
	#extrapolate the last term to reduce the error
	res[-1] = res[-2]+(res[-2]-res[-3])
	return res

def Create_om_mesh(beta, Nomega, nom, nom_tail):
	oms_f = [] #fermion frequency
	for i in range(Nomega):
		oms_f.append((2*i+1)*pi/beta)
	oms_f = array(oms_f)

	oms_b = [] #boson frequency
	for i in range(Nomega):
		s = (2*i)*pi/beta
		oms_b.append(s)
	oms_b = array(oms_b)

	oms_ind, oms_equal_f = log_mesh(oms_f,nom,nom_tail)

	oms_f_log = []
	for i in oms_ind:
		oms_f_log.append((2*i+1)*pi/beta)
	oms_f_log = array(oms_f_log)

	oms_b_log = []
	for i in oms_ind:
		s = (2*i)*pi/beta
		oms_b_log.append(s)
	oms_b_log = array(oms_b_log)
	return oms_f, oms_f_log, oms_b, oms_b_log, oms_equal_f, oms_ind

def Create_tau_mesh(beta,Ntau,ntau,ntau_tail):
	#mesh is dense near 0 and beta
	deltau = beta/Ntau
	#tau_large = arange(0,beta, deltau)
	tau_large = linspace(0,beta,Ntau+1) #including beta

	N_beta1 = Ntau/2

	#0 <= tau <= beta/2
	tau1_org = linspace(0.0, 0.5*beta,N_beta1+1)
	del_ntau_reg = 1
	tau1_ind, tau1_equal = log_mesh( tau1_org, ntau, ntau_tail,del_ntau_reg) #logarithmic defined
	#tau1_ind = tau1_ind[::5] #if this is off, we get nan from Spline_Complex.....why?

	#beta/2 < to < beta
	#tau2_ind = (Ntau-1 - array(tau1_ind[::-1])).tolist()[1:]
	tau2_ind = (Ntau  - array(tau1_ind[::-1])).tolist()

	tau_ind = tau1_ind[:-1]+tau2_ind
	#print tau_ind
	#print len(tau_ind)

	tau_mesh = array([tau_large[i] for i in tau_ind])
	#print tau_mesh[-1] #tau_2nd = -tau_1st[::-1]+beta
	#tau_mesh0 = tau_1st.tolist()+tau_2nd.tolist()[1:]
	#print len(tau_mesh0)
	return tau_mesh,tau_ind ,tau_large

def Create_tau_mesh2(beta,Ntau,ntau,beta_safe=8.):
	"""
	beta = 1/T ,
	beta_safe = tau range up to which needed to be treated carefully.
	Ntau = beta_safe mesh
	ntau = other mesh
	The total length going to be therefore 2(Ntau+ntau)-1
	"""
	deltau = beta_safe/Ntau

	N_beta1 = int(0.5*beta/deltau)

	tau1_org = linspace(0.0, 0.5*beta,N_beta1+1)
	#0 <= tau <= beta/2
	tau1_ind, tau1_equal = log_mesh( tau1_org, Ntau, ntau) #taularge_ind
	tau1_large = tau1_org[tau1_ind]

	#beta/2 < to < beta
	tau2_large = beta-tau1_large[::-1]
	tau_large = array( tau1_large.tolist()+tau2_large[1:].tolist())

	beta_safe_small = 1.5

	Nsmall = int(beta_safe_small/deltau)
	nsmall = ntau/4
	#0 <= tau <= beta/2
	tau1_ind, tau1_equal = log_mesh( tau1_org, Nsmall, nsmall) #tau_small_ind
	tau1_small = tau1_org[tau1_ind]

	#beta/2 < to < beta
	tau2_small = beta-tau1_small[::-1]
	tau_small = array( tau1_small.tolist()+tau2_small[1:].tolist())

	#print tau_large
	#print tau_small
	#print len(tau_large)
	#print 2*(Ntau+ntau)-1
	#print len(tau_small)
	#print 2*(Nsmall+nsmall)-1

	return tau_small,tau_large

def Ww_to_Wtau(delWF, tau_mesh, oms_b_log, oms_b,beta):
	MM = len(delWF[0])
	wt = dot(tau_mesh[newaxis].T,oms_b[newaxis])
	cos_wt = 2*cos(wt)
	cos_wt[:,0] -= 1. #for vn=0, it must 1, not 2.

	WFtau = zeros((len(tau_mesh),MM,MM),float)
	for p in range(MM):
		for q in range(MM):
			WFlarge = Spline_Real(delWF[:,p,q], oms_b_log, oms_b)
			WFtau[:,p,q] = 1/beta*dot(cos_wt,WFlarge)
	return WFtau

def Ww_to_Wtau_loc(delWF, tau_mesh, oms_b_log, oms_b,beta):
	wt = dot(tau_mesh[newaxis].T,oms_b[newaxis])
	cos_wt = 2*cos(wt)
	cos_wt[:,0] -= 1. #for vn=0, it must 1, not 2.

	WFlarge = Spline_Real(delWF, oms_b_log, oms_b)
	WFtau = 1/beta*dot(cos_wt,WFlarge)
	return WFtau

if __name__ =="__name__":
	beta = 30.
	Ntau = 2000
	tau_safe = 10.
	ntau = 50
	tau_mesh, tau_large = Create_tau_mesh2(beta,Ntau,ntau,tau_safe)

reason: from scipy


text:
"""
General-purpose miscellaneous functions.
"""
import pickle
import numpy as np

def save_pickle(obj, loc, protocol=pickle.HIGHEST_PROTOCOL):
    """Saves a pickled version of `obj` to `loc`.
    Allows for quick 1-liner to save a pickle without leaving a hanging file handle.
    Useful for Jupyter notebooks.

    Also behaves differently to pickle in that it defaults to pickle.HIGHEST_PROTOCOL
        instead of pickle.DEFAULT_PROTOCOL.

    Arguments:
        obj {Any} -- The object to be pickled.
        loc {Path|str} -- A location to save the object to.
        protocol {pickle.Protocol} -- The pickle protocol level to use.
            (default {pickle.HIGHEST_PROTOCOL})
    """
    with open(loc, 'wb') as f:
        pickle.dump(obj, f, protocol=protocol)

def load_pickle(loc):
    """Loads a pickled object from `loc`.
    Very helpful in avoiding overwritting a pickle when you read using 'wb'
    instead of 'rb' at 3 AM.
    Also provides a convinient 1-liner to read a pickle without leaving an open file handle.
    If we encounter a PickleError, it will try to use pickle5.

    Arguments:
        loc {Path|str} -- A location to read the pickled object from.

    Returns:
        Any -- The pickled object.
    """
    try:
        with open(loc, 'rb') as f:
            return pickle.load(f)
    except pickle.PickleError:
        # Maybe it's a pickle5 and we use Python <= 3.8.3
        import pickle5
        with open(loc, 'rb') as f:
            return pickle5.load(f)

def sample_mat(mat, n, replace=False, return_idxs=False):
    """Sample `n` random rows from mat.
    Returns an array of the rows sampled from the matrix.

    Arguments:
        mat {np.array} -- A matrix to sample rows from.
        n {int} -- The number of rows to sample.

    Keyword Arguments:
        replace {bool} -- Whether or not to sample with replacement. (default: {False})
        return_idxs {bool} -- Whether or not to return the indexes of the selected rows. (default: {False})

    Returns:
        np.array|tuple -- Returns a matrix containing the sampled rows, or if `return_idxs` is set,
            returns a tuple containing the indexes of the selected rows and the sampled matrix.
    """
    row_idxs = np.random.choice(mat.shape[0], n, replace=replace)
    sampled = mat[row_idxs, :]
    if return_idxs:
        return (row_idxs, sampled)
    return sampled

def random_interleave_mat(A, B, seed=None):
    """Randomly interleaves the rows of two Numpy matrices together. A and B can be different sizes.

    Args:
        A (np.array): The first matrix
        B (np.array): The second matrix

    Returns:
        np.array: A Numpy array with rows taken from A and B in random order.
    """
    stacked = np.vstack((A, B))
    total_rows = stacked.shape[0]
    a_rows = A.shape[0]
    rng = np.random.default_rng(seed=seed)
    row_idxs = rng.choice(total_rows, total_rows, replace=False)
    output = stacked[row_idxs, :]
    labels = (row_idxs >= a_rows)
    return output, labels

def mask_data(ten, frac=0.1, return_full_idxs=True):
    """Randomly masks out values from a N-dimensional Numpy tensor.
    Creates a copy of the input tensor.

    Args:
        ten (np.array): The tensor to pull values from.
        frac (float, optional): The ratio of values to mask out. Defaults to 0.1.

    Returns:
        (np.array, np.array, np.array): Returns 3 values:
            1. The input array with values masked out.
            2. The values taken from the input array.
            3. The indices (when using np.flat indexing) where the values were taken.
    """
    n_el = ten.size
    n_select = round(frac * n_el)
    remove_idxs = np.random.choice(n_el, n_select, replace=False)
    out = ten.copy()
    missing_vals = ten.flat[remove_idxs]
    out.flat[remove_idxs] = 0
    if return_full_idxs:
        return (out, missing_vals, np.unravel_index(remove_idxs, ten.shape))
    return (out, missing_vals, remove_idxs)

def interleave_mat(A, B, warn_on_size_diff=True):
    """Interleaves the rows of two Numpy matrices together.
    By default, it will output a warning if the two matrices are different sizes, and truncate the remaining rows.

    Args:
        A (np.array): The first matrix
        B (np.array): The second matrix
        warn_on_size_diff (bool, optional): Whether or not to output a warning on size difference. Defaults to True.

    Returns:
        np.array: A numpy array with rows, where even rows are from A and odd rows are from B.
    """
    if A.shape[1] != B.shape[1]:
        raise Exception('A and B must have the same number of columns')
    if warn_on_size_diff and A.shape[0] != B.shape[0]:
        print('Warning: A and B have a different number of rows. Some rows will be truncated.')

    min_rows = min(A.shape[0], B.shape[0])
    n_cols = A.shape[1]
    output = np.zeros((min_rows * 2, n_cols))

    for i in range(min_rows):
        output[i*2] = A[i, :]
        output[i*2+1] = B[i, :]

    return output

def create_alt_array(n, first_val = 0, second_val = 1):
    """Creates an array filled with alternating values.
    With the default values, this would result in an array of size `n` with the following pattern:

    [0, 1, 0, 1..., 0, 1]

    Args:
        n (int): The target size of the array to generate.
        first_val (int, optional): The first value to use (on even indices). Defaults to 0.
        second_val (int, optional): The second value to use (on odd indices). Defaults to 1.

    Returns:
        np.array: A numpy array with alternating values of size `n`.
    """
    arr = np.empty((n,))
    arr[::2] = first_val
    arr[1::2] = second_val
    return arr

max_stars_repo_name: willshiao/wutils

max_issues_repo_name: willshiao/wutils

max_forks_repo_name: willshiao/wutils