text:
# Author : skjang@mnc.ai
# Date : 2020-12-03

import logging
import numpy as np
import pandas as pd
import os, glob, time, datetime
import pickle
import gzip
import copy
import json

import cv2
import random
import torch
from sklearn.model_selection import train_test_split
from sklearn.model_selection import StratifiedKFold, StratifiedShuffleSplit
import torch.distributed as dist
from torch.utils.data.dataset import Dataset
from torch.utils.data.distributed import DistributedSampler
from torch.nn.parallel import DistributedDataParallel
from multiprocessing import Pool

from torch.optim import lr_scheduler
import torch.optim as optim
import torch.nn as nn
import torchvision.transforms as A
from PIL import Image
from efficientnet_pytorch import EfficientNet

import matplotlib.pyplot as plt
import pydicom as dicom
from pydicom.pixel_data_handlers.numpy_handler import unpack_bits
from torch.utils.tensorboard import SummaryWriter

disease = ['Sinusitis','Oral_cancer'][0]

formatter = logging.Formatter('[%(asctime)s] [%(name)s] [%(levelname)s] %(message)s')

def get_logger(name, level=logging.DEBUG, resetlogfile=False, path='log'):
    fname = os.path.join(path, name+'.log')
    os.makedirs(path, exist_ok=True)
    if resetlogfile :
        if os.path.exists(fname):
            os.remove(fname)
    logger = logging.getLogger(name)
    logger.handlers.clear()
    logger.setLevel(level)

    ch = logging.StreamHandler()
    ch.setLevel(level)
    ch.setFormatter(formatter)
    logger.addHandler(ch)

    fh = logging.FileHandler(fname)
    fh.setLevel(level)
    fh.setFormatter(formatter)
    logger.addHandler(fh)
    return logger

today_datever = datetime.datetime.now().strftime("%y%m%d")
logger = get_logger(f'{disease}_EfficientNet_{today_datever}', resetlogfile=True)
logger.setLevel(logging.INFO)

def prepare_metatable(filenames):
    cvdf = pd.DataFrame(filenames)
    cvdf.columns = ['filename']
    cvdf['patientID'] = cvdf['filename'].apply(lambda x : x.split('/')[-1].split('_')[0])
    #cvdf['year'] = cvdf['filename'].apply(lambda x : x.split('/')[-2].split('_')[1])

    cvdf['left_label_org'] = cvdf['filename'].apply(lambda x : x.replace('.dcm', '').replace('.jpg', '').split('_')[-2])
    cvdf['right_label_org'] = cvdf['filename'].apply(lambda x : x.replace('.dcm', '').replace('.jpg', '').split('_')[-1])
    #print(pd.crosstab(cvdf['left_label_org'], cvdf['right_label_org'], margins=True))

    cvdf['left_label'] = '1'
    cvdf.at[cvdf['left_label_org']=='0','left_label'] = '0'
    cvdf.at[cvdf['left_label_org']=='x','left_label'] = 'x'
    cvdf['right_label'] = '1'
    cvdf.at[cvdf['right_label_org']=='0','right_label'] = '0'
    cvdf.at[cvdf['right_label_org']=='x','right_label'] = 'x'
    #print(pd.crosstab(cvdf['left_label'], cvdf['right_label']))

    cvdf['FOLD'] = np.nan

    oldcolumns = cvdf.columns.tolist()
    cvdf['index'] = cvdf.index
    cvdf = cvdf[['index']+oldcolumns]
    return cvdf

def save_validationlist(root='.'):
    # list up filenames of valid data
    # totalfiles = glob.glob(os.path.join(root,"test_20??_withUPID","*.dcm"))
    # filenames = glob.glob(os.path.join(root,"test_20??_withUPID","*_[0-3]_[0-3].dcm"))
    data_dir = ["final_dcm","final_crop"][0]
    logger.info('['*10+' '*20 +'START ANALYSIS'+' '*20+ ']'*10)
    filenames = glob.glob(os.path.join(root,data_dir,"*" + (".dcm" if data_dir=='final_dcm' else '.jpg')))
    logger.info(f'No. of total datasets : {len(filenames)} patients') # 6516

    rmfn = glob.glob(os.path.join(root,data_dir,"*_x_x"+(".dcm" if data_dir=='final_dcm' else '.jpg')))
    if len(rmfn)>1:
        logger.info('  x_x.dcm  :')
        logger.info(rmfn)
        filenames.remove(rmfn)
    logger.info(f'No. of valid datasets : {len(filenames)} patients (excluded x_x.dcm )') #2980 (20.10.7 ver)

    cvdf = prepare_metatable(filenames)

    n_folds = 10
    plen = len(filenames)
    logger.info(f'----- Split patients for {n_folds} Cross-validation')

    skf = StratifiedKFold(n_splits=n_folds, random_state=42, shuffle=True)
    for ii, (train_pindex, test_pindex) in enumerate(skf.split(range(plen),cvdf['left_label'])):
        # record fold index
        cvdf.at[test_pindex,'FOLD']= ii
        cvdf[f'FOLD{ii}_testset'] = 0
        cvdf.at[test_pindex,f'FOLD{ii}_testset'] = 1

    # save metadata
    filelist_dir = os.path.join(root,'inputlist')
    os.makedirs(filelist_dir, exist_ok=True)

    cvdf.to_csv(os.path.join(filelist_dir,"input_metadata_table.csv"),index=False)
    cvdf[['index','filename']].to_csv(os.path.join(filelist_dir,"input_filenames_total.csv"),index=False)
    for i in range(n_folds):
        cvdf.loc[cvdf[f'FOLD{i}_testset']==1,'filename'].to_csv(os.path.join(filelist_dir,f"input_filenames_fold{i}.csv"),index=False)

    # statistics
    logger.info(f'----- Data statistics by fold',cvdf['FOLD'].value_counts())

    logger.info(cvdf['FOLD'].value_counts())

    labelfreq_left = pd.crosstab(cvdf['FOLD'], cvdf['left_label'], margins=True)
    labelfreq_left_ratio = pd.crosstab(cvdf['FOLD'], cvdf['left_label'], margins=True, normalize='index')
    labelfreq_right = pd.crosstab(cvdf['FOLD'], cvdf['right_label'], margins=True)
    labelfreq_right_ratio = pd.crosstab(cvdf['FOLD'], cvdf['right_label'], margins=True, normalize='index')

    labelfreq = pd.concat([labelfreq_left, labelfreq_right], axis=1, keys=['left_sinus', 'right_sinus'], names=[' ','label'])
    labelfreq_ratio = pd.concat([labelfreq_left_ratio, labelfreq_right_ratio], axis=1, keys=['left_sinus', 'right_sinus'], names=[' ','label (ratio)'])

    labelfreq.to_csv(os.path.join(filelist_dir,f"label_freq_byfold.csv"))
    labelfreq_ratio.to_csv(os.path.join(filelist_dir,f"label_freq_ratio_byfold.csv"),float_format = '%.2f')
    logger.info(f'----- Label frequency by fold')
    logger.info(labelfreq)
    logger.info(f'----- Label frequency (ratio) by fold')
    logger.info(labelfreq_ratio)

import multiprocessing

class ImageDataset(Dataset):
    def __init__(self, root='.', input_csv='inputlist/input_filenames_fold', annotation_path=None,
    fold_num=0, data_type='train', carrydata=True, transform=None, savejpg=True):
        super(ImageDataset, self).__init__()
        self.root = root
        self.input_csv = input_csv
        self.annotation_path = annotation_path
        self.fold_num = fold_num
        self.data_type = data_type # train, val, test
        self.carrydata = carrydata
        self.transform = transform
        self.savejpg = savejpg

        if self.annotation_path is not None:
            json_file = open(self.annotation_path)
            roi_annotation = json.load(json_file) #coco
            json_file.close()

            self.roi_dict = dict()
            for segmentation in roi_annotation:
                image_name = list(segmentation.keys())[0].replace('.dcm','')
                bbox_dict = list(segmentation.values())[0][-1]
                assert bbox_dict['name']=='bounding box'
                self.roi_dict[image_name] = bbox_dict['points'] # {'00001334_0_0' : [128.91, 230, 920.48, 786.83]}

        logger.info('--'*20)
        logger.info(f"- Build {self.data_type} dataset")
        logger.info('-- Transform')
        logger.info(self.transform)
        logger.info('-- ')

        n_folds = 10

        train_fold = list(range(n_folds))
        val_fold = train_fold[fold_num-1]
        train_fold.remove(fold_num) # test set
        train_fold.remove(val_fold) # validation set

        if data_type=="train":
            self.filenames = []
            for i in train_fold:
                fl_i = pd.read_csv(f'{input_csv}{i}.csv')['filename'].tolist()
                self.filenames.extend(fl_i)
        elif data_type=="val":
            self.filenames = pd.read_csv(f'{input_csv}{val_fold}.csv')['filename'].tolist()
        elif data_type=="test":
            self.filenames = pd.read_csv(f'{input_csv}{fold_num}.csv')['filename'].tolist()
        # self.filenames = glob.glob(os.path.join(self.root,"test_20??_withUPID","*_[0-3]_[0-3].dcm"))

        logger.info(f" Read {len(self.filenames)} files : '{self.filenames[0]}', etc.") # tmp

        self.imagedata = []
        self.targetdata = [] # binary class (0 to 1)
        self.orginal_targetdata = [] # multiclass (0 to 3)

        save_dir = 'preprocessed_data'
        if os.path.exists(save_dir):
            logger.info(f"preprocessed_data folder already exists. Use data in this folder...")
        else:
            logger.info(f"Create preprocessed_data folder...")
            os.makedirs(save_dir, exist_ok=True)

        if self.carrydata:
            if self.savejpg:
                os.makedirs('./png_tmp', exist_ok=True)

            n_cpu = multiprocessing.cpu_count()
            used_th = 8*6
            logger.info(f'no. cpu existed : {n_cpu}, use {used_th} threads')
            pool = multiprocessing.Pool(processes=used_th)
            result = pool.map(self.read_data, self.filenames)
            pool.close()
            pool.join()

            for i in result :
                self.imagedata.extend(i[1])
                self.orginal_targetdata.extend(i[2])

            self.targetdata = [0 if x==0 else 1 for x in self.orginal_targetdata]
            logger.info(f' No. of datasets loaded  : {len(self.imagedata)} images')

    def __len__(self):
        return len(self.imagedata)

    def __getitem__(self, index):
        if self.carrydata:
            img = self.imagedata[index]
            target = self.targetdata[index]
        else:
            _, img, target = self.read_data(self.filenames[index])
        if self.transform is not None:
            img = Image.fromarray(img)
            img = self.transform(img)
        return img, target

    def patients_length(self):
        return len(self.filenames)

    def read_data(self, fn, pkl_dir='preprocessed_data'):
        targets = fn.split('.')[-2].split('_')[-2:]
        patient_id = fn.split('/')[-1].split('_')[0]

        pkl_fname = os.path.join(pkl_dir, f"{fn.replace('.dcm', '').replace('./', '').replace('/', '__')}.pkl")

        if os.path.exists(pkl_fname):
            # load pickle file
            logger.info(f' > Load preprocessed data : patient # {patient_id}')
            with gzip.open(pkl_fname,'rb') as f:
                [patient_id, imgs, targets] = pickle.load(f)
        else:
            logger.info(f' > Load & preprocess data : patient # {patient_id}')
            img_format = fn.split('.')[-1]
            if img_format == 'dcm' or img_format == 'DCM' :
                d = dicom.read_file(fn)
                org_img = d.pixel_array.astype('float')

                logger.info(f' >> min of intensity : {np.min(org_img)}, max: {np.max(org_img)} (patient # {patient_id}, {d.PhotometricInterpretation})')

                if d.PhotometricInterpretation == 'MONOCHROME1':
                    org_img = np.max(org_img) - org_img

                if self.annotation_path is not None:
                    # crop roi
                    img_name = fn.split('/')[-1].split('.')[0]
                    roi_location = self.roi_dict[img_name]
                    roi_location = [max(x, 0) for x in roi_location] # for avoiding negative location
                    x1, y1, x2, y2 = roi_location
                    cropped_img = org_img[int(y1):int(y2),int(x1):int(x2)]

                    plt.imsave(f'png_tmp/{img_name}_original.jpg',org_img, cmap=plt.cm.gray)
                    org_img = cropped_img.copy()

                    os.makedirs('final_crop_fromDicom',exist_ok=True)
                        #cv2.imwrite(f'final_crop_fromDicom/{img_name}.jpg', org_img)
                    plt.imsave(f'final_crop_fromDicom/{img_name}.jpg',cropped_img, cmap=plt.cm.gray)
            elif img_format == 'jpg' or img_format == 'jpeg' or img_format == 'JPG' or \
            img_format == 'png' or img_format == 'PNG':
                org_img = Image.open(fn).convert('L')
                org_img = np.array(org_img)
            if (org_img.shape[0]==0) or (org_img.shape[1]==0):
                logger.warn(f'error in {fn}')
                logger.warn(org_img.shape) #(1174, 1144)
                logger.error(cropped_img.shape) #(574, 0)

            ##tmp : mimic crooped image
            # h, w = org_img.shape
            # h__, w__ = h//10, w//10
            # org_img = org_img[2*h__:8*h__, 2*w__:8*w__] # cropped
            #plt.imsave('./png_tmp/'+patient_id +'-original.png',d.pixel_array, cmap=plt.cm.bone)

            norm_img = self.standardize_image(org_img) ##float64
            imgs = self.split_leftright(norm_img)

            logger.info(f' >> Original image shape : {org_img.shape[0]} x {org_img.shape[1]}') # 1101 x 1134 -> 696 x 636
            logger.info(f' >> Patch shape after split_leftright : {imgs.shape[0]} x {imgs.shape[1]} x {imgs.shape[2]}') # expected : 2 x 318 x 318 # tmp

            # visualize processed image
            if self.savejpg:
                plt.imsave('./png_tmp/'+patient_id +'-cropped.png', org_img, cmap=plt.cm.bone) # tmp

                if targets[0]!='x':
                    plt.imsave('./png_tmp/'+patient_id +'-left.png',imgs[0,:,:], cmap=plt.cm.bone) # tmp
                if targets[1]!='x':
                    plt.imsave('./png_tmp/'+patient_id +'-right.png',imgs[1,:,:], cmap=plt.cm.bone) # tmp

            # drop data if label=='x'
            notx_index = [i for i in range(len(targets)) if targets[i]!='x']
            targets = [targets[i] for i in notx_index]
            targets = [int(x) for x in targets]

            imgs = imgs[notx_index,:,:]

            # save pickle file
            objs = [patient_id, imgs, targets]
            with gzip.open(pkl_fname, 'wb') as f:
                pickle.dump(objs, f)
        if 'x' in targets:
                logger.error(fn)
                logger.error(np.where(targets=='x'))
                logger.error('target contains unvalid value')
        return patient_id, imgs, targets

    def split_leftright(self, img_arr):
        def make_square(img_l):
            if img_l.shape[0]>img_l.shape[1]:
                img_l = img_l[1:,:]
            if img_l.shape[1]>img_l.shape[0]:
                img_l = img_l[:,1:]
            return img_l
        h, w = img_arr.shape
        h_, w_ = h//2, w//2
        img_l = img_arr[:, w_:]
        img_r = img_arr[:, :w_]
        img_l = cv2.resize(img_l, (300,300))
        img_r = cv2.resize(img_r, (300,300))
        imgs = np.stack((img_l,img_r),axis=0)
        return imgs

    def image_minmax(self, img):
        img_minmax = img.copy()

        if np.ndim(img) == 3:
            for d in range(3):
                img_minmax[d] = image_minmax(img[d])
        else:
            img_minmax = ((img - np.min(img)) / (np.max(img) - np.min(img))).copy()

        return img_minmax

    def standardize_image(self, img):
        normalized_img = img.copy()

        if np.ndim(img) == 3:
            for d in range(3):
                normalized_img[d] = standardize_image(img[d])
        else:
            normalized_img = ((img - np.mean(img)) / np.std(img)).copy()

        return normalized_img

def get_dataloaders(args, dataset_class, batch, root, fold_num=0, multinode=False, augment=False):

    logger.info(f'----- Get data loader : Validation No. {fold_num}')

    # augmentation
    transform_test = A.Compose(
        [
            A.ToTensor()
        ]
    )

    if augment:
        transform_train = A.Compose(
            [
                A.ShiftScaleRotate(shift_limit=0.05, scale_limit=0.05, rotate_limit=10, p=0.5),
                A.RandomBrightnessContrast(p=0.5),
                #A.Resize(300,300),
                ToTensorV2()
            ]
        )
    else:
        transform_train = transform_test

    # prepare training dataset using transform_train
    if args.trained_model is not None:
        train_dataset = None
    else:
        train_dataset = dataset_class(data_type='train', fold_num=fold_num, transform=transform_train, annotation_path='selected_list(6516).json')

    # prepare validation(for hyperparam. tuning), test sets using transform_test
    if args.trained_model is not None:
        train_dataset = None
    else:
        val_dataset = dataset_class(data_type='val', fold_num=fold_num, transform=transform_test, annotation_path='selected_list(6516).json')
    test_dataset = dataset_class(data_type='test', fold_num=fold_num, transform=transform_test, annotation_path='selected_list(6516).json')

    train_sampler, valid_sampler, test_sampler = None, None, None
    trainloader, validloader, testloader = None, None, None

    #define sampler for imbalanced data
    if args.trained_model is None :

        target_list = torch.tensor(train_dataset.targetdata)
        target_list = target_list[torch.randperm(len(target_list))]
        class_sample_count = np.array([len(np.where(target_list==t)[0]) for t in np.unique(target_list)])
        class_weights = 1./class_sample_count
        class_weights = torch.from_numpy(class_weights.copy()).type('torch.DoubleTensor')

        logger.info('class_weights : ')
        logger.info(class_weights)

        class_weights_all = class_weights[target_list]

        train_sampler = torch.utils.data.sampler.WeightedRandomSampler(
            weights=class_weights_all,
            num_samples=len(class_weights_all),
            replacement=True
        )

        # data loader
        trainloader = torch.utils.data.DataLoader(
            train_dataset, batch_size=batch, shuffle=True if train_sampler==None else False, pin_memory=True,
            sampler=train_sampler, drop_last=True)

        validloader = torch.utils.data.DataLoader(
            val_dataset, batch_size=batch, shuffle=False, pin_memory=True,
            sampler=valid_sampler, drop_last=False)

    testloader = torch.utils.data.DataLoader(
        test_dataset, batch_size=batch, shuffle=False, pin_memory=True,
        sampler=test_sampler, drop_last=False)

    return train_dataset, trainloader, validloader, testloader


text:
""" Tracking Utilities

Classes:

* :py:class:`Link`: The main track object to store timeseries of cells

Functions:

* :py:func:`find_flat_regions`: Smooth tracks to find persistent regions
* :py:func:`position_to_velocity`: Convert a track of coordinates to velocities
* :py:func:`smooth_velocity`: Smooth the velocity field
* :py:func:`rolling_window`: Create strided windows over numpy arrays

I/O Functions:

* :py:func:`load_simulation_tsvfile`: Load tracks from Demarcus's simulations
* :py:func:`load_track_csvfile`: Load tracks from the tracking algorithm
* :py:func:`save_track_csvfile`: Save tracks for the tracking algorithm

API Documentation
-----------------

"""

# Imports
import pathlib
from typing import Tuple, List, Optional

# 3rd party
import numpy as np

from scipy.ndimage.measurements import label
from scipy.interpolate import UnivariateSpline, interp1d

from skimage.morphology import remove_small_holes, remove_small_objects

# Our own imports
from . import rolling_interp, rolling_slope
from ..utils import BBox

# Classes

class Link(object):
    """ Link object in a chain

    :param \\*args:
        The initial link coordinates (index, tp, x, y)
    """

    def __init__(self, *args):
        if len(args) == 0:
            self.index = None
            self.line_t = []
            self.line_x = []
            self.line_y = []
        else:
            index, tp, x, y = args
            self.index = index
            self.line_t = [tp]
            self.line_x = [x]
            self.line_y = [y]

    @property
    def first(self):
        return self.line_t[0], self.line_x[0], self.line_y[0]

    @property
    def last(self):
        return self.line_t[-1], self.line_x[-1], self.line_y[-1]

    @classmethod
    def merge_mean(cls, *args):
        """ Take the average over several tracks """
        if not all([len(t) == len(args[0]) for t in args]):
            raise ValueError('All tracks must have the same length')

        merge_tt = None
        merge_xx = []
        merge_yy = []
        for track in args:
            tt, xx, yy = track.to_arrays()
            if merge_tt is None:
                merge_tt = tt
            else:
                if not np.allclose(tt, merge_tt):
                    raise ValueError('Time vectors must be aligned')
            merge_xx.append(xx)
            merge_yy.append(yy)
        return cls.from_arrays(tt, np.mean(merge_xx, axis=0), np.mean(merge_yy, axis=0))

    @classmethod
    def join(cls, *args, **kwargs):
        """ Merge two or more links into a single link """
        interp = kwargs.get('interp', None)

        link = cls()
        for a in args:
            if len(link.line_t) >= 2 and interp == 'linear':
                tstep = link.line_t[-1] - link.line_t[-2]
                t0 = link.line_t[-1]
                x0 = link.line_x[-1]
                y0 = link.line_y[-1]
                dt = a.line_t[0] - t0
                dx = a.line_x[0] - x0
                dy = a.line_y[0] - y0

                num_steps = int(round(dt / tstep))

                for i in range(1, num_steps):
                    link.line_t.append(i * dt / num_steps + t0)
                    link.line_x.append(i * dx / num_steps + x0)
                    link.line_y.append(i * dy / num_steps + y0)

            link.line_t += a.line_t
            link.line_x += a.line_x
            link.line_y += a.line_y

        return link

    @classmethod
    def from_tuples(cls, tuples):
        """ Create a link from a list of tuples """

        chain = cls()
        for t, x, y in tuples:
            chain.add(-1, int(t), float(x), float(y))
        return chain

    @classmethod
    def from_arrays(cls, tt, xx, yy):
        """ Create a link from a set of arrays

        :param ndarray tt:
            The time numpy array
        :param ndarry xx:
            The x coordinate array
        :param ndarray yy:
            The y coordinate array
        :returns:
            The link class from the numpy arrays
        """

        assert tt.shape[0] == xx.shape[0]
        assert tt.shape[0] == yy.shape[0]

        chain = cls()
        chain.line_t = list(tt)
        chain.line_x = list(xx)
        chain.line_y = list(yy)
        return chain

    def get_track_length(self, space_scale: float = 1.0) -> float:
        """ Get the arc length along the track

        :returns:
            A float giving the total distance along the track
        """
        _, x, y = self.to_arrays()
        dx = x[1:] - x[:-1]
        dy = y[1:] - y[:-1]
        return np.sum(np.sqrt(dx**2 + dy**2)) * space_scale

    def get_track_displacement(self, space_scale: float = 1.0) -> float:
        """ Get the distance from track start to track end

        :returns:
            A float giving the total displacement of the track
        """
        dx = self.line_x[-1] - self.line_x[0]
        dy = self.line_y[-1] - self.line_y[0]
        return np.sqrt(dx**2 + dy**2) * space_scale

    def get_bbox(self) -> BBox:
        """ Get the bounding box for this track

        :returns:
            The bounding box for this track
        """
        min_x = np.min(self.line_x)
        max_x = np.max(self.line_x)
        min_y = np.min(self.line_y)
        max_y = np.max(self.line_y)
        return BBox(min_x, min_y, max_x, max_y)

    def should_add(self, index, tp, x, y) -> bool:
        """ See if this set of points matches

        :param index:
            The unique key that should match this track
        :param int tp:
            The timepoint for this set
        :param float x:
            The x coordinate for this set
        :param float y:
            The y coordinate for this set
        :returns:
            True if this point fits with the current index
        """
        if self.index is None:
            return True
        dx = self.line_x[-1] - x
        dy = self.line_y[-1] - y
        if index == self.index and np.sqrt(dx**2 + dy**2) < 1e-5:
            return True
        return False

    def add(self, index, tp: int, x: float, y: float):
        """ Add a timepoint

        :param index:
            A unique key for this track
        :param int tp:
            The timepoint for this set
        :param float x:
            The x coordinate for this set
        :param float y:
            The y coordinate for this set
        """
        self.index = index
        self.line_t.append(tp)
        self.line_x.append(x)
        self.line_y.append(y)

    def find_timepoint(self, timepoint: int):
        """ Find the x,y coordinates for a timepoint

        :param int timepoint:
            The timepoint to find
        :returns:
            x, y for the given timepoint or None, None if no match
        """
        try:
            idx = self.line_t.index(timepoint)
        except ValueError:
            return None, None
        else:
            return self.line_x[idx], self.line_y[idx]

    def to_padded_arrays(self,
                         min_t: int = 0,
                         max_t: int = -1,
                         extrapolate: bool = False) -> Tuple[np.ndarray]:
        """ Pad out the arrays in time

        :param int min_t:
            The minimum time index or 0
        :param int max_t:
            The maximum time index or -1 for the last value in our array
        :param bool interpolate:
            If True, interpolate holes in the samples
        :returns:
            The t, x, y numpy arrays, padded with nans
        """
        line_t, line_x, line_y = self.to_arrays()
        if min_t < 0:
            min_t = line_t[min_t]
        if max_t < 0:
            max_t = line_t[max_t]
        pad_t = np.arange(min_t, max_t, 1)

        if extrapolate:
            fill_value_x = (line_x[0], line_x[-1])
            fill_value_y = (line_y[0], line_y[-1])
        else:
            fill_value_x = fill_value_y = np.nan
        fxn_x = interp1d(line_t, line_x,
                         kind='linear',
                         bounds_error=False,
                         fill_value=fill_value_x)
        fxn_y = interp1d(line_t, line_y,
                         kind='linear',
                         bounds_error=False,
                         fill_value=fill_value_y)
        return pad_t, fxn_x(pad_t), fxn_y(pad_t)

    def to_arrays(self) -> Tuple[np.ndarray]:
        """ Convert the internal pointset to numpy arrays

        :returns:
            The t, x, y numpy arrays
        """
        return (np.array(self.line_t),
                np.array(self.line_x),
                np.array(self.line_y))

    def vel_x(self) -> np.ndarray:
        """ X-coordinate velocity

        :returns:
            An (n-1) x 1 array of velocities in the x-axis
        """
        x = np.array(self.line_x)
        t = np.array(self.line_t)
        if x.shape[0] < 2:
            return np.array([0.0])

        dt = t[1:] - t[:-1]
        dx = x[1:] - x[:-1]
        return dx / dt

    def vel_y(self) -> np.ndarray:
        """ Y-coordinate velocity

        :returns:
            An (n-1) x 1 array of velocities in the y-axis
        """
        y = np.array(self.line_y)
        t = np.array(self.line_t)
        if y.shape[0] < 2:
            return np.array([0.0])

        dt = t[1:] - t[:-1]
        dy = y[1:] - y[:-1]
        return dy / dt

    def vel_mag(self) -> np.ndarray:
        """ Velocity magnitude for the point set

        :returns:
            An (n-1) x 1 array of velocity magnitudes
        """
        return np.sqrt(self.vel_x()**2 + self.vel_y()**2)

    def scale_points(self,
                     x_scale: float = 1.0,
                     y_scale: float = 1.0,
                     x_center: Optional[float] = None,
                     y_center: Optional[float] = None):
        """ Scale the track by a factor around its center of mass

        :param float x_scale:
            Scale x by this amount
        :param float y_scale:
            Scale y by this amount
        """
        x = np.array(self.line_x)
        y = np.array(self.line_y)
        if x_center is None:
            x_center = np.mean(x)
        if y_center is None:
            y_center = np.mean(y)

        scaled_x = (x - x_center) * x_scale
        scaled_y = (y - y_center) * y_scale

        self.line_x = list(scaled_x + x_center)
        self.line_y = list(scaled_y + y_center)

    def interpolate_points(self,
                           step: int = 1,
                           smoothing: Optional[int] = None):
        """ Interpolate the points back onto a regular grid

        :param int step:
            The timestep to interpolate over
        :param int smoothing:
            How strongly to smooth the line, s=0 is linear, s=3 is default
        """

        out_t = np.arange(self.line_t[0], self.line_t[-1]+step, step)

        # ext - 3 - return the boundary value
        # Default is ext=0, extrapolate, which does bad things
        x_fit = UnivariateSpline(self.line_t, self.line_x, s=smoothing, ext=3)
        out_x = x_fit(out_t)

        y_fit = UnivariateSpline(self.line_t, self.line_y, s=smoothing, ext=3)
        out_y = y_fit(out_t)

        self.line_t = list(out_t)
        self.line_x = list(out_x)
        self.line_y = list(out_y)

    def __getitem__(self, idx):
        return self.line_t[idx], self.line_x[idx], self.line_y[idx]

    def __eq__(self, other):
        if len(self) != len(other):
            return False
        return all([np.allclose(self.line_t, other.line_t),
                    np.allclose(self.line_x, other.line_x),
                    np.allclose(self.line_y, other.line_y)])

    def __len__(self):
        return len(self.line_x)

    def __repr__(self):
        return f'L({self.line_x},{self.line_y})'

# Functions

def find_flat_regions(tt: np.ndarray,
                      yy: np.ndarray,
                      interp_points: int = 100,
                      cutoff: int = 10,
                      noise_points: int = 5) -> List[np.ndarray]:
    """ Look for regions in of directional persistence or change

    :param ndarray tt:
        The time variable to search over
    :param ndarray yy:
        The space variable to search
    :param int interp_points:
        The number of points to use for the rolling slope fit
    :param float cutoff:
        The largest allowable slope magnitude
    :param int noise_points:
        Number of points to treat a connected region as noise
    :returns:
        A list of masks for each flat region
    """
    slopes = rolling_slope(tt, tt, yy, interp_points)

    # Find all slopes smaller than the cutoff and then convert to individual masks
    mask = (np.abs(slopes) < cutoff).astype(np.bool)
    if np.any(mask) and not np.all(mask):
        mask = remove_small_objects(mask, noise_points)
        mask = remove_small_holes(mask, noise_points)
        rois = label(mask)[0]
    else:
        rois = mask

    # Split up the connected region labels
    regions = []
    for roi_idx in np.unique(rois):
        if roi_idx < 1:
            continue
        regions.append(rois == roi_idx)
    return regions

def smooth_velocity(tt: np.ndarray, xx: np.ndarray, yy: np.ndarray,
                    resample_factor: int = 10,
                    interp_points: int = 7,
                    smooth_points: int = 5) -> Tuple[np.ndarray]:
    """ Rolling window velocity smoothing

    :param ndarray tt:
        The 1D array of times
    :param ndarray xx:
        The 1D array of x positions
    :param ndarray yy:
        The 1D array of y positions
    :param int resample_factor:
        How many subsamples to take between samples
    :param int interp_points:
        How many support points to use in the linear interpolation
    :param int smooth_points:
        How many support points to use in the rolling average smoothing
    :returns:
        A new set of time, x, y coordinate arrays for the smoothed track
    """

    # Interpolate the points onto a larger space with linear smoothing
    sm_tt = np.linspace(np.min(tt), np.max(tt), resample_factor*tt.shape[0])
    sm_xx = rolling_interp(sm_tt, tt, xx, interp_points, order=1)
    sm_yy = rolling_interp(sm_tt, tt, yy, interp_points, order=1)

    # Remove bumps in the points with rolling averaging
    if smooth_points > 0:
        sm_xx = np.mean(rolling_window(sm_xx, smooth_points*resample_factor, pad='same'), 1)
        sm_yy = np.mean(rolling_window(sm_yy, smooth_points*resample_factor, pad='same'), 1)
    return sm_tt, sm_xx, sm_yy

def position_to_velocity(tt: np.ndarray,
                         xx: np.ndarray,
                         yy: np.ndarray,
                         pad: str = 'valid') -> Tuple[np.ndarray]:
    """ Do the basic position to velocity metric conversion

    :param ndarray tt:
        The 1D array of times
    :param ndarray xx:
        The 1D array of x positions
    :param ndarray yy:
        The 1D array of y positions
    :returns:
        The deltas for time, postion, velocity, and angle
    """
    dt = (tt[1:] - tt[:-1])
    dx = (xx[1:] - xx[:-1])
    dy = (yy[1:] - yy[:-1])
    if pad == 'same':
        dt = np.concatenate([dt, dt[-1]])
        dx = np.concatenate([dx, dx[-1]])
        dy = np.concatenate([dy, dy[-1]])

    ds = np.sqrt((dx**2 + dy**2))
    dv = ds / dt
    dtheta = np.arctan2(dy, dx)
    return dt, dx, dy, ds, dv, dtheta

def load_simulation_tsvfile(track_tsvfile: pathlib.Path) -> List[Link]:
    """ Load the track data from Demarcus' simulations

    :param Path track_tsvfile:
        The tab-separated track data to load
    :returns:
        The list of tracks in the file
    """
    cells = {}

    with track_tsvfile.open('rt') as fp:
        for i, line in enumerate(fp):
            if i == 0:
                continue
            line = line.strip()
            if line == '':
                continue
            rec = [r.strip() for r in line.split('\t')]
            if len(rec) != 7:
                raise ValueError(f'Invalid line: {rec}')

            # Unpack the record
            cell_tp = float(rec[0])
            cell_id = int(rec[1])
            cell_x = float(rec[2])
            cell_y = float(rec[3])

            cells.setdefault(cell_id, Link()).add(cell_id, cell_tp, cell_x, cell_y)

    return list(cells.values())

def save_track_csvfile(track_csvfile: pathlib.Path,
                       final_chains: List[Link]):
    """ Save the track CSV file

    :param Path track_csvfile:
        The track file to save
    :param list[Link] final_chains:
        The tracks to save
    """

    track_csvfile.parent.mkdir(exist_ok=True, parents=True)

    with track_csvfile.open('wt') as fp:
        for i, chain in enumerate(final_chains):
            fp.write(f'Track{i:03d}\n')
            chain_t = ','.join(f'{int(t):d}' for t in chain.line_t)
            fp.write(chain_t + '\n')
            chain_x = ','.join(f'{x:0.4f}' for x in chain.line_x)
            fp.write(chain_x + '\n')
            chain_y = ','.join(f'{y:0.4f}' for y in chain.line_y)
            fp.write(chain_y + '\n')

def load_track_csvfile(track_csvfile: pathlib.Path,
                       min_len: int = 0) -> List[Link]:
    """ Load the track csvfile back

    :param Path track_csvfile:
        The file to load
    :param int min_len:
        The minimum length of a track
    :returns:
        A list of Link objects, one for each track
    """

    final_chains = []
    cur_chain_t = None
    cur_chain_x = None
    cur_chain_y = None

    num_tracks = -1

    def append_last_chain():
        if all([c is not None for c in [cur_chain_t, cur_chain_x, cur_chain_y]]):
            assert len(cur_chain_t) == len(cur_chain_x)
            assert len(cur_chain_t) == len(cur_chain_y)

            if min_len <= 0 or len(cur_chain_t) >= min_len:
                chain = Link()
                chain.line_t = list(cur_chain_t)
                chain.line_x = list(cur_chain_x)
                chain.line_y = list(cur_chain_y)
                final_chains.append(chain)

    with track_csvfile.open('rt') as fp:
        for i, line in enumerate(fp):
            line = line.split('#', 1)[0].strip()
            if line == '':
                continue
            if line.startswith('Track'):
                num_tracks += 1
                append_last_chain()
                cur_chain_t = cur_chain_x = cur_chain_y = None
            elif cur_chain_t is None:
                cur_chain_t = [int(t.strip()) for t in line.split(',')]
            elif cur_chain_x is None:
                cur_chain_x = [float(x.strip()) for x in line.split(',')]
            elif cur_chain_y is None:
                cur_chain_y = [float(y.strip()) for y in line.split(',')]
            else:
                raise ValueError(f'Invalid line #{i}: {line}')
    append_last_chain()
    return final_chains

def rolling_window(a, window, pad='valid'):
    """ Rolling window over an array

    :param ndarray a:
        The array to roll a window over
    :param int window:
        The window size
    :returns:
        A window stride object
    """
    orig_shape = a.shape[-1]
    if pad == 'same':
        window_left = window//2
        window_right = window_left + orig_shape
        na = np.empty(a.shape[:-1] + (orig_shape + window, ))
        na[..., window_left:window_right] = a
        na[..., :window_left] = a[..., 0]
        na[..., window_right:] = a[..., -1]
        a = na

    shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
    strides = a.strides + (a.strides[-1],)
    res = np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)
    if pad == 'same':
        res = res[:orig_shape]
    return res

ext: py

max_stars_repo_name: JackToppen/deep-hipsc-tracking

max_issues_repo_name: JackToppen/deep-hipsc-tracking

max_forks_repo_name: JackToppen/deep-hipsc-tracking

reason: import numpy,from scipy


text:
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from os.path import join
from pathlib import Path

# Root, data, and result directories
project_root = str(Path(__file__).parent.parent)
data_dir = join(project_root, 'data')
raw_data_dir = join(data_dir, 'raw')
result_dir = join(project_root, 'results')

# Raw migration data
gcd_sample_raw_file_path = join(raw_data_dir, 'Mig_Mat_Nordic_SME_LC_RE.xlsx')
gcd_raw_file_path = join(raw_data_dir, 'gcd_nordic_migration.csv')
seb_raw_file_path = join(raw_data_dir, 'seb_raw_data.xlsx')

# Raw indicator data
indicator_small_raw_path = join(raw_data_dir, 'systemic_indicators_small_2.xlsx')
indicator_large_raw_path = join(raw_data_dir, 'Swe_Fin_Nor_Den_EconomicData_WorldBank_2280395.xlsx')

# Processed migration data
gcd_sample_path = join(data_dir, 'gcd_sample.csv')
gcd_data_path = join(data_dir, 'gcd_data.csv')
seb_data_path = join(data_dir, 'seb_data.csv')

# Processed indicator data
indicator_small_path = join(data_dir, 'indicator_small.csv')
indicator_large_path = join(data_dir, 'indicator_large.csv')
indicator_pca_path = join(data_dir, 'indicator_pca.csv')

# Rating and industry code data
nace_gcd_raw_path = join(raw_data_dir, 'GCD_nace.csv')
rating_table_path = join(data_dir, 'rating_table.csv')
nace_seb_table_path = join(data_dir, 'nace_seb_table.csv')
nace_gcd_table_path = join(data_dir, 'nace_gcd_table.csv')

# Custom industry code data
code_gcd_path = join(data_dir, 'code_table_gcd.csv')
code_seb_path = join(data_dir, 'code_table_seb.csv')

# Result files
portfolio_result_path = join(result_dir, 'portfolio_rmse.csv')
portfolio_kf_result_path = join(result_dir, 'portfolio_kf_rmse.csv')
sector_result_path = join(result_dir, 'sector_rmse.csv')
sector_grouped_result_path = join(result_dir, 'sector_grouped_rmse.csv')
country_sector_result_path = join(result_dir, 'country_sector_rmse.csv')

# Path dictionaries
raw_path_dict = {'gcd_sample': gcd_sample_raw_file_path, 'gcd': gcd_raw_file_path,  'seb': seb_raw_file_path, 'nace_gcd': nace_gcd_raw_path,
                    'indicator_small': indicator_small_raw_path, 'indicator_large': indicator_large_raw_path}
processed_path_dict = {'gcd_sample': gcd_sample_path, 'gcd':gcd_data_path, 'seb': seb_data_path,
                        'rating': rating_table_path, 'nace_seb': nace_seb_table_path, 'nace_gcd': nace_gcd_table_path,
                        'code_gcd': code_gcd_path, 'code_seb': code_seb_path,
                        'indicator_small': indicator_small_path, 'indicator_large': indicator_large_path, 'indicator_pca': indicator_pca_path, 'indicator_pca': indicator_pca_path,
                        'portfolio_result': portfolio_result_path, 'portfolio_kf': portfolio_kf_result_path,
                        'sector_result': sector_result_path, 'sector_grouped_result': sector_grouped_result_path, 'country_sector_result': country_sector_result_path}

code_types = ['0', '1', '2', 'nace']

def load_raw_data_file(dataset):
    '''Loading raw data files in pandas ExcelFile
    '''
    return pd.ExcelFile(raw_path_dict[dataset])

def process_data(dataset=[], rating='seb'):
    '''Main function for processing the data, the following was done:
        removing unnecessary columns,
        renaming columns,
        mapping ratings and industries,
        filtering and normalizing indicators,
        applying PCA on the large indicator set.
    The function body is removed to comply with NDAs.
    '''
    pass

def load_data(dataset):
    return pd.read_csv(processed_path_dict[dataset], keep_default_na=False)

def get_gcd_transition_matrix(year=[], quarter=[4], asset_class=[], country=[], industry=[], code=0, rating='seb', count=True, sample=False):
    '''
    Function for creating a transition matrix for given years, quarter, asset class, country, industry from GCD data set
    '''

    if sample:
        df = load_data('gcd_sample')
    else:
        df = load_data('gcd')
    code_dict = {0: 'industry_0', 1: 'industry_gcd'}

    # Filtering
    if year:
        df = df[df['year'].isin(year)]
    if quarter:
        df = df[df['quarter'].isin(quarter)]
    if asset_class:
        df = df[df['asset_class'].isin(asset_class)]
    if country:
        df = df[df['country'].isin(country)]
    if industry:
        df = df[df[code_dict[code]].isin(industry)]

    # Creating the transition matrix
    if rating == 'seb':
        matrix_size = (15,16)
    else:
        matrix_size = (24,24)

    m = np.zeros(matrix_size, dtype=int)

    for i in range(1, matrix_size[0] + 1):
        for j in range(1, matrix_size[1] + 1):
            df_ = df[df['t1']==i]
            df_ = df_[df_['t2']==j]
            m[i-1,j-1] = int(df_['count'].sum())

    if not count:
        m = (m.T / np.sum(m, axis=1)).T

    return m

def get_seb_transition_matrix(year=[], industry=[], country=[], code=0, count=True):
    '''
    Function for creating a transition matrix for given years, quarter, asset class, country, industry from SEB data set
    '''
    df = load_data('seb')
    code_dict = {0: 'industry_0', 1: 'industry_1', 2: 'industry_2'}

    # Filtering
    if year:
        df = df[df['year'].isin(year)]
    if country:
        df = df[df['country'].isin(country)]
    if industry:
        df = df[df[code_dict[code]].isin(industry)]

    # Creating the transition matrix
    matrix_size = (15,16)
    m = np.zeros(matrix_size, dtype=int)

    for i in range(1, matrix_size[0] + 1):
        for j in range(1, matrix_size[1] + 1):
            df_ = df[df['t1']==i]
            df_ = df_[df_['t2']==j]
            m[i-1,j-1] = int(df_['count'].sum())

    if not count:
        m = (m.T / np.sum(m, axis=1)).T

    return m

def get_code_name(code, table, code_type=0):
    '''
    Helper function for getting industry names from the industry code table
    '''

    code_col = 'code_' + code_types[code_type]
    name_col = 'name_' + code_types[code_type]
    code_table = load_data(table)[[code_col, name_col]]
    code_table = code_table[code_table[code_col]==code]
    code_table.drop_duplicates(inplace=True)
    return list(code_table[name_col])

def get_px(m):
    '''
    Proportion of ratings, p(X=x)
    '''

    m = m + 1e-6
    return np.sum(m, axis=1) / np.sum(m)

def get_pdx(m):
    '''
    Conditional default probability knowing rating, p(D=1|X=x)
    '''

    m = m + 1e-10
    return m[:,-1] / (np.sum(m, axis=1) + 1e-4)

def get_pd(m):
    '''
    Marginal probability of default p(D=1)
    '''

    m = m + 1e-10
    return np.sum(get_px(m) * get_pdx(m)) + 1e-10

def get_transitions(data, years=[2008, 2020], country=[], quarter=[4], industry=[], code=0):
    '''
    Getting list of transition matrices of portfolio, country, or industry in chronological order
    '''

    M = []

    for year in list(range(*years)):
        if data == 'seb':
            M.append(get_seb_transition_matrix(year=[year], country=country, industry=industry, code=code))
        elif data == 'gcd':
            M.append(get_gcd_transition_matrix(year=[year], quarter=quarter, country=country, industry=industry, code=code))
        elif dataset == 'gcd_sample':
            M.append(get_gcd_transition_matrix(year=[year], quarter=quarter, country=country, industry=industry, code=code, sample=True))

    return np.array(M, dtype=int)

def get_z(data, years=[2008, 2020], country='SE', pca_n=0):
    '''
    Getting the Z factor matrix (shape = years x num_indicators)
    '''

    years = list(range(*years))

    if data == 'large':
        df = load_data('indicator_large')
    elif data == 'pca':
        df  = load_data('indicator_pca')
    else:
        df = load_data('indicator_small')
        if data == 'single':
            df = df[['country', 'year', 'GDP']]

    df = df[df['country'] == country]
    df = df[df['year'].isin(years)]
    df = df.sort_values('year')
    df = df.drop(columns=['country', 'year'])
    Z = df.to_numpy()

    if data == 'pca' and pca_n > 0:
        Z = Z[:, :pca_n]
    return Z

max_stars_repo_name: trananhthong/sebcreditrisk

max_issues_repo_name: trananhthong/sebcreditrisk

max_forks_repo_name: trananhthong/sebcreditrisk

reason: import numpy


text:
from tkinter import Image

import cv2
import pytesseract
import mysql.connector
import imutils
import numpy

# Create Database Connection
from PIL.ImageOps import grayscale

mydb = mysql.connector.connect(host="localhost", user="root", passwd="");
mycursor = mydb.cursor()
mycursor.execute("use attendance")

# import pytesseract
pytesseract.pytesseract.tesseract_cmd = r"C:\\Program Files\\Tesseract\\tesseract.exe";

# Load image
x = int(input("Enter the day : "));

if (x == 1):
    img = cv2.imread('1.png')
    query = 'UPDATE attend set day1=%s where id=%s'
    mycursor.execute("ALTER TABLE `attend` ADD `day1` VARCHAR(50) NOT NULL AFTER `name`;")
if (x == 2):
    img = cv2.imread('2.png')
    query = 'UPDATE attend set day2=%s where id=%s'
    mycursor.execute("ALTER TABLE `attend` ADD `day2` VARCHAR(50) NOT NULL AFTER `day1`;")
if (x == 3):
    img = cv2.imread('3.png')
    query = 'UPDATE attend set day3=%s where id=%s'
    mycursor.execute("ALTER TABLE `attend` ADD `day3` VARCHAR(50) NOT NULL AFTER `day2`;")
if (x == 4):
    img = cv2.imread('4.png')
    query = 'UPDATE attend set day4=%s where id=%s'
    mycursor.execute("ALTER TABLE `attend` ADD `day4` VARCHAR(50) NOT NULL AFTER `day3`;")
if (x == 5):
    img = cv2.imread('5.png')
    query = 'UPDATE attend set day5=%s where id=%s'
    mycursor.execute("ALTER TABLE `attend` ADD `day5` VARCHAR(50) NOT NULL AFTER `day4`;")

# First Iteration
# Crop Image and get index number
x = 80
y = 275
w = 120
h = 20
cropImg1 = img[y:y + h, x:x + w]

img1 = cv2.cvtColor(cropImg1, cv2.COLOR_BGR2RGB);
image = imutils.resize(img1, width=612, height=800)
id = (pytesseract.image_to_string(image));
id = id.strip()
cv2.imshow("Segmentation",cropImg1)
cv2.imshow("BGR to RGB",img1)
cv2.imshow("Resacle Image",image)

# crop image and get Signature
sigx = 530
sigy = 275
sigw = 135
sigh = 20
cropSigImg = img[sigy:sigy + sigh, sigx:sigx + sigw]
cv2.imshow("Segmentation Signature",cropSigImg)
grayImage1 = cv2.cvtColor(cropSigImg, cv2.COLOR_BGR2GRAY)
cv2.imshow("Grays Scale",grayImage1)
(thresh, blackAndWhiteImage) = cv2.threshold(grayImage1, 127, 255, cv2.THRESH_BINARY)
cv2.imshow("Segmentation Signature",cropSigImg)
cv2.imshow("Grays Scale",grayImage1)
cv2.imshow("black & white",blackAndWhiteImage)
whitepix = numpy.sum(blackAndWhiteImage == 0)
# print(whitepix);
if (whitepix > 100):
    status = str("present")
else:
    status = str("absence")
print(id)
print(status)

mycursor.execute(query, (status, id));

# Second Iteration
# Crop Image and get index number
x = 80
y = 300
w = 120
h = 20
cropImg1 = img[y:y + h, x:x + w]
img1 = cv2.cvtColor(cropImg1, cv2.COLOR_BGR2RGB);
image = imutils.resize(img1, width=612, height=800)
id = (pytesseract.image_to_string(image));
id = id.strip()
cv2.imshow("Segmentation",cropImg1)
cv2.imshow("BGR to RGB",img1)
cv2.imshow("Resacle Image",image)

# crop image and get signature
sigx = 530
sigy = 305
sigw = 135
sigh = 20
cropSigImg = img[sigy:sigy + sigh, sigx:sigx + sigw]
grayImage1 = cv2.cvtColor(cropSigImg, cv2.COLOR_BGR2GRAY)
(thresh, blackAndWhiteImage) = cv2.threshold(grayImage1, 127, 255, cv2.THRESH_BINARY)
cv2.imshow("Segmentation Signature",cropSigImg)
cv2.imshow("Grays Scale",grayImage1)
cv2.imshow("black & white",blackAndWhiteImage)
whitepix = numpy.sum(blackAndWhiteImage == 0)
# print(whitepix);
if (whitepix > 100):
    status = str("present")
else:
    status = str("absence")
print(id)
print(status)

mycursor.execute(query, (status, id));

# Third Iteration
# Crop Image and get index number
x = 80
y = 330
w = 120
h = 20
cropImg1 = img[y:y + h, x:x + w]
img1 = cv2.cvtColor(cropImg1, cv2.COLOR_BGR2RGB);
image = imutils.resize(img1, width=612, height=800)
id = (pytesseract.image_to_string(image));
id = id.strip()
cv2.imshow("Segmentation",cropImg1)
cv2.imshow("BGR to RGB",img1)
cv2.imshow("Resacle Image",image)

# crop image and get signature
sigx = 530
sigy = 330
sigw = 135
sigh = 20
cropSigImg = img[sigy:sigy + sigh, sigx:sigx + sigw]
grayImage1 = cv2.cvtColor(cropSigImg, cv2.COLOR_BGR2GRAY)
(thresh, blackAndWhiteImage) = cv2.threshold(grayImage1, 127, 255, cv2.THRESH_BINARY)
cv2.imshow("Segmentation Signature",cropSigImg)
cv2.imshow("Grays Scale",grayImage1)
cv2.imshow("black & white",blackAndWhiteImage)
whitepix = numpy.sum(blackAndWhiteImage == 0)
# print(whitepix);
if (whitepix > 100):
    status = str("present")
else:
    status = str("absence")
print(id)
print(status)

mycursor.execute(query, (status, id));

# Fourth Iteration
# Crop Image and get index number
x = 80
y = 355
w = 120
h = 20
cropImg1 = img[y:y + h, x:x + w]
img1 = cv2.cvtColor(cropImg1, cv2.COLOR_BGR2RGB);
image = imutils.resize(img1, width=612, height=800)
id = (pytesseract.image_to_string(image));
id = id.strip()
cv2.imshow("Segmentation",cropImg1)
cv2.imshow("BGR to RGB",img1)
cv2.imshow("Resacle Image",image)

# crop image and get signature
sigx = 530
sigy = 360
sigw = 135
sigh = 20
cropSigImg = img[sigy:sigy + sigh, sigx:sigx + sigw]
grayImage1 = cv2.cvtColor(cropSigImg, cv2.COLOR_BGR2GRAY)
(thresh, blackAndWhiteImage) = cv2.threshold(grayImage1, 127, 255, cv2.THRESH_BINARY)
cv2.imshow("Segmentation Signature",cropSigImg)
cv2.imshow("Grays Scale",grayImage1)
cv2.imshow("black & white",blackAndWhiteImage)
whitepix = numpy.sum(blackAndWhiteImage == 0)
# print(whitepix);
if (whitepix > 100):
    status = str("present")
else:
    status = str("absence")
print(id)
print(status)

mycursor.execute(query, (status, id));

# Fifth Iteration
# Crop Image and get index number
x = 80
y = 383
w = 120
h = 20
cropImg1 = img[y:y + h, x:x + w]
img1 = cv2.cvtColor(cropImg1, cv2.COLOR_BGR2RGB);
image = imutils.resize(img1, width=612, height=800)
id = (pytesseract.image_to_string(image));
id = id.strip()
cv2.imshow("Segmentation",cropImg1)
cv2.imshow("BGR to RGB",img1)
cv2.imshow("Resacle Image",image)

# crop image and get signature
sigx = 530
sigy = 383
sigw = 135
sigh = 20
cropSigImg = img[sigy:sigy + sigh, sigx:sigx + sigw]
grayImage1 = cv2.cvtColor(cropSigImg, cv2.COLOR_BGR2GRAY)
(thresh, blackAndWhiteImage) = cv2.threshold(grayImage1, 127, 255, cv2.THRESH_BINARY)
cv2.imshow("Segmentation Signature",cropSigImg)
cv2.imshow("Grays Scale",grayImage1)
cv2.imshow("black & white",blackAndWhiteImage)
whitepix = numpy.sum(blackAndWhiteImage == 0)
# print(whitepix);
if (whitepix > 100):
    status = str("present")
else:
    status = str("absence")
print(id)
print(status)

mycursor.execute(query, (status, id));

# sixth Iteration
# Crop Image and get index number
x = 80
y = 409
w = 120
h = 20
cropImg1 = img[y:y + h, x:x + w]
img1 = cv2.cvtColor(cropImg1, cv2.COLOR_BGR2RGB)
image = imutils.resize(img1, width=612, height=800)
id = (pytesseract.image_to_string(image));
id = id.strip()
cv2.imshow("Segmentation",cropImg1)
cv2.imshow("BGR to RGB",img1)
cv2.imshow("Resacle Image",image)

# crop image and get signature
sigx = 530
sigy = 409
sigw = 135
sigh = 20
cropSigImg = img[sigy:sigy + sigh, sigx:sigx + sigw]
grayImage1 = cv2.cvtColor(cropSigImg, cv2.COLOR_BGR2GRAY)
(thresh, blackAndWhiteImage) = cv2.threshold(grayImage1, 127, 255, cv2.THRESH_BINARY)
cv2.imshow("Segmentation Signature",cropSigImg)
cv2.imshow("Grays Scale",grayImage1)
cv2.imshow("black & white",blackAndWhiteImage)
whitepix = numpy.sum(blackAndWhiteImage == 0)
# print(whitepix);
if (whitepix > 100):
    status = str("present")
else:
    status = str("absence")
print(id)
print(status)

mycursor.execute(query, (status, id));

cv2.imshow('Result',blackAndWhiteImage);
cv2.waitKey(0);


text:
import numpy as np
import schafkopf.players.data.load_data as ld
import schafkopf.players.data.encodings as enc
from schafkopf.game_modes import PARTNER_MODE
from schafkopf.ranks import KING, TEN, OBER, NINE, EIGHT, SEVEN
from schafkopf.suits import BELLS, ACORNS, LEAVES, HEARTS

filepath = 'test_data.p'
game_data_dic = {'player_hands': [[(4, 0), (7, 1), (6, 1), (5, 1), (2, 1), (5, 3), (5, 2), (6, 0)],
                                  [(3, 2), (3, 1), (0, 1), (7, 3), (6, 3), (2, 2), (0, 2), (7, 0)],
                                  [(4, 3), (4, 2), (4, 1), (3, 0), (1, 1), (2, 3), (0, 3), (6, 2)],
                                  [(3, 3), (1, 3), (7, 2), (1, 2), (5, 0), (2, 0), (1, 0), (0, 0)]],
                 'game_mode': (PARTNER_MODE, LEAVES),
                 'played_cards': [((5, 2), 0), ((0, 2), 1), ((6, 2), 2), ((7, 2), 3), ((1, 2), 3), ((7, 1), 0),
                                  ((2, 2), 1), ((3, 0), 2), ((4, 3), 2), ((3, 3), 3), ((2, 1), 0), ((0, 1), 1),
                                  ((4, 2), 2), ((5, 0), 3), ((6, 1), 0), ((3, 1), 1), ((0, 3), 2), ((1, 3), 3),
                                  ((5, 3), 0), ((7, 3), 1), ((6, 3), 1), ((2, 3), 2), ((2, 0), 3), ((6, 0), 0),
                                  ((7, 0), 1), ((4, 1), 2), ((0, 0), 3), ((5, 1), 0), ((1, 1), 2), ((1, 0), 3),
                                  ((4, 0), 0), ((3, 2), 1)],
                 'declaring_player': 2}

def test_num_games_in_file():
    assert ld.num_games_in_file(filepath) == 2751

def test_num_augmented_examples_in_file():
    assert ld.num_augmented_games_in_file(filepath) == 28260

def test_prepare_data_bidding():
    data_list = ld.prepare_data_bidding(game_data_dic)
    assert len(data_list) == 4
    x, y = data_list[0]
    assert x.shape == (8, 32)
    assert y.shape == (9, )
    data_list_augmented = ld.prepare_data_bidding(game_data_dic, augment_data=True)
    assert len(data_list_augmented) == 24

def test_prepare_data_inference():
    card_sequences, hands_to_predict = ld.prepare_data_inference(game_data_dic, num_samples=26)
    assert len(card_sequences) == 26
    assert len(hands_to_predict) == 26
    seq = card_sequences[0]
    hands = hands_to_predict[0]
    assert hands.shape == (128, )
    assert seq.shape == (28, 36)
    shortest_seq = None
    pl_hands = None
    for s, i in zip(card_sequences, range(len(hands_to_predict))):
        if np.array_equal(s[2][:32], np.zeros(32)):
            shortest_seq = s
            pl_hands = hands_to_predict[i]
    played_card = shortest_seq[1][:32]
    assert pl_hands.shape == (128, )
    assert enc.decode_one_hot_card(played_card) == (KING, LEAVES)
    goal = np.zeros(128)
    indices = [16, 29, 25, 21, 9, 23, 24, 46, 45, 33, 63, 59, 42, 34, 60, 83,
               82, 81, 76, 69, 75, 67, 90, 111, 103, 126, 102, 116, 104, 100, 96]
    for index in indices:
        goal[index] = 1
    for i in range(128):
        assert goal[i] == pl_hands[i], 'failed at card {}'.format(i)
    assert np.array_equal(goal, pl_hands)

def test_prepare_extended_data_inference():
    card_sequences, player_hands, hands_to_predict = ld.prepare_extended_data_inference(game_data_dic, num_samples=26)
    assert len(player_hands) == 26
    opp_hands = None
    chosen_seq = None
    pl_hand = None
    for seq, hand, opp in zip(card_sequences, player_hands, range(len(hands_to_predict))):
        if np.array_equal(seq[21][:32], np.zeros(32)) and not np.array_equal(seq[20][:32], np.zeros(32)):
            chosen_seq = seq
            pl_hand = hand
            opp_hands = hands_to_predict[opp]

    assert opp_hands is not None
    assert chosen_seq is not None
    assert pl_hand is not None
    opp_hand_1 = enc.decode_hand_inference(opp_hands[:32])
    opp_hand_2 = enc.decode_hand_inference(opp_hands[32:64])
    opp_hand_3 = enc.decode_hand_inference(opp_hands[64:])
    assert set(opp_hand_1) == {(TEN, BELLS), (KING, HEARTS), (OBER, BELLS)}
    assert set(opp_hand_2) == {(NINE, ACORNS), (OBER, HEARTS), (EIGHT, HEARTS)}
    assert set(opp_hand_3) == {(NINE, BELLS), (EIGHT, BELLS), (SEVEN, BELLS)}
    assert set(enc.decode_on_hot_hand(pl_hand)) == {(3, 2), (3, 1), (0, 1), (7, 3), (6, 3), (2, 2), (0, 2), (7, 0)}
    assert chosen_seq[20][31] == 1
    assert chosen_seq[20][33] == 1

max_stars_repo_name: Taschee/schafkopf

max_forks_repo_name: Taschee/schafkopf


text:
import sys
import numpy as np
import cv2 as cv
from scipy import signal

Alpha_F = 0.1
Alpha_L = 1.0
Alpha_T = 0.3

V_F = 0.5
V_L = 0.2
V_T = 20.0

Num = 10
Beta = 0.1

W = [0.5, 1, 0.5, 1, 0, 1, 0.5, 1, 0.5,]
W = np.array(W, np.float).reshape((3, 3))
M = W

global F
global L
global Y
global T
global Y_AC

def main(argv):
    default_file = 'images/3.jpeg'
    filename = argv[0] if len(argv) > 0 else default_file

    # Loads an image
    src = cv.imread(cv.samples.findFile(filename), cv.IMREAD_GRAYSCALE)

    # Check if image is loaded fine
    if src is None:
        print ('Error opening image!')
        print ('Usage: hough_lines.py [image_name -- default ' + default_file + '] \n')
        return -1

    cv.imshow("Original", src)

    dim = src.shape

    F = np.zeros( dim, np.float)
    L = np.zeros( dim, np.float)
    Y = np.zeros( dim, np.float)
    T = np.ones( dim, np.float)
    Y_AC = np.zeros( dim, np.float)

    #normalize image
    S = cv.normalize(src.astype('float'), None, 0.0, 1.0, cv.NORM_MINMAX)

    for cont in range(Num):
        #numpy.convolve(W, Y, mode='same')
        F = np.exp(-Alpha_F) * F + V_F * signal.convolve2d(Y, W, mode='same') + S
        L = np.exp(-Alpha_L) * L + V_L * signal.convolve2d(Y, M, mode='same')
        U = F * (1 + Beta * L)
        T = np.exp(-Alpha_T) * T + V_T * Y
        Y = (U>T).astype(np.float)
        Y_AC = Y_AC + Y

    cv.imshow("Result", Y)
    Y_AC = cv.normalize(Y_AC.astype('float'), None, 0.0, 1.0, cv.NORM_MINMAX)
    cv.imshow("Result Acumulated", Y_AC)

    #Y = (Y*255).astype(np.uint8)
    cv.imwrite('result.jpg', Y*255)
    cv.imwrite('result Acumulated.jpg', Y_AC*255)

    cv.waitKey()
    return 0

if __name__ == "__main__":
    main(sys.argv[1:])


text:
import numpy as np
from PIL import Image
from matplotlib import pyplot as plt
from numpy.fft import fft2, fftshift, ifft2

# %% read images
img_raw = np.asarray(Image.open('Q5_2.tif'))
row, col = img_raw.shape

# %% fft
img_fourier = fft2(np.pad(img_raw, ((0, row), (0, col))))
img_view = np.log10(fftshift(np.abs(img_fourier)) + 1).astype(np.uint8)

plt.subplot(121)
plt.imshow(img_raw, cmap='gray')
plt.title('Raw Q5_2.tif')

plt.subplot(122)
plt.imshow(img_view, cmap='gray')
plt.title('Fourier Transformed Q5_2.tif')
# plt.savefig('Q5_2.png')
plt.show()

# %% low pass filtering
def gaussian(d0: int):
    def d(mu: int, nu: int):
        return np.sqrt((mu - row) ** 2 + (nu - col) ** 2)

    gaussian_lpf = np.fromfunction(lambda mu, nu: np.exp(-d(mu, nu) ** 2 / (2 * d0 ** 2)),
                                   (2 * row, 2 * col))
    img_freq_filtered = img_fourier * fftshift(gaussian_lpf)
    img_freq_filtered_view = np.log10(fftshift(np.abs(img_freq_filtered)) + 1).astype(np.uint8)
    img_filtered = np.real(ifft2(img_freq_filtered))[0:row, 0:col].astype(np.uint8)

    # gaussian_lpf = np.ones((row * 2, col * 2)) - gaussian_lpf

    return img_filtered, gaussian_lpf

filtered10, filter10 = gaussian(10)
filtered30, filter30 = gaussian(30)
filtered60, filter60 = gaussian(60)
filtered160, filter160 = gaussian(160)

plt.subplot(121)
plt.imshow(img_raw, cmap='gray')
plt.title('Raw Q5_2.tif')

plt.subplot(122)
plt.imshow(filtered10, cmap='gray')
# plt.title('Gaussian low pass filtered Q5_2.tif\n$D_0=30$')
plt.title('Gaussian high-pass filtered Q5_2.tif\n$D_0=10$')
# plt.savefig('Q5_2_10_high.png')
plt.show()

# %%
plt.figure(figsize=(12.8, 4.8 * 2))
plt.subplot(231)
plt.imshow(filtered30, cmap='gray')
plt.title('Gaussian low-pass filtered image\n$D_0=30$')
# plt.title('Gaussian high-pass filtered image\n$D_0=30$')

plt.subplot(232)
plt.imshow(filtered60, cmap='gray')
plt.title('Gaussian low-pass filtered image\n$D_0=60$')
# plt.title('Gaussian high-pass filtered image\n$D_0=60$')

plt.subplot(233)
plt.imshow(filtered160, cmap='gray')
plt.title('Gaussian low-pass filtered image, $D_0=160$')
# plt.title('Gaussian high-pass filtered image\n$D_0=160$')

plt.subplot(234)
plt.imshow(filter30, cmap='gray')
plt.title('Gaussian low-pass filter, $D_0=30$')
# plt.title('Gaussian high-pass filter\n$D_0=30$')

plt.subplot(235)
plt.imshow(filter60, cmap='gray')
plt.title('Gaussian low-pass filter, $D_0=60$')
# plt.title('Gaussian high-pass filter\n$D_0=60$')

plt.subplot(236)
plt.imshow(filter160, cmap='gray')
plt.title('Gaussian low-pass filter, $D_0=160$')
# plt.title('Gaussian high-pass filter\n$D_0=160$')
plt.savefig('Q5_2_filtered_low.png')
# plt.savefig('Q5_2_filtered_high.png')

plt.show()

# # %%
# plt.subplot(121)
# plt.imshow(img_raw, cmap='gray')
# plt.title('Raw Q5_2.tif')
#
# plt.subplot(122)
# plt.imshow(filtered30, cmap='gray')
# plt.title('Gaussian low pass filtered Q5_2.tif\n$D_0=30$')
# plt.savefig('Q5_2_30.png')
# plt.show()
#
# plt.subplot(121)
# plt.imshow(img_raw, cmap='gray')
# plt.title('Raw Q5_2.tif')
#
# plt.subplot(122)
# plt.imshow(filtered60, cmap='gray')
# plt.title('Gaussian low pass filtered Q5_2.tif\n$D_0=60$')
# plt.savefig('Q5_2_60.png')
# plt.show()
#
# plt.subplot(121)
# plt.imshow(img_raw, cmap='gray')
# plt.title('Raw Q5_2.tif')
#
# plt.subplot(122)
# plt.imshow(filtered30, cmap='gray')
# plt.title('Gaussian low pass filtered Q5_2.tif\n$D_0=160$')
# plt.savefig('Q5_2_160.png')
# plt.show()

max_stars_repo_name: kommunium/dip-lab

max_issues_repo_name: kommunium/dip-lab

max_forks_repo_name: kommunium/dip-lab

reason: import numpy,from numpy


text:
#!/usr/bin/env python3

import rospy
import numpy as np
import utils
import threading, time, collections, heapq, itertools, math, os
from itertools import count
from utils import Circle, Path, SearchNode, SearchNodeTree, TreeNode

from geometry_msgs.msg import PoseStamped, PolygonStamped
from nav_msgs.msg import Odometry, OccupancyGrid
from visualization_msgs.msg import MarkerArray, Marker
from std_msgs.msg import ColorRGBA
from nav_msgs.srv import GetMap

import matplotlib.pyplot as plt

class HeuristicSearch(object):
	""" Abstract class for heuristic search on the provided set of cost/admissibility/heuristic/neighbor functions """
	def __init__(self, make_tree=False):
		self.make_tree = make_tree

	def reset(self, start_state):
		# Frontier is a priority queue.
		self.closed_set = []
		self.frontier = []
		self.found_paths = []
		self.goal_state = self.goal()
		self.step_count = 0
		if self.make_tree:
			self.tree_root = TreeNode(state=start_state, children=[])
			ss = SearchNodeTree(state=start_state, parent=None, cost=0, tree_node=self.tree_root, \
						heuristic=self.heuristic(start_state, self.goal_state))
		else:
			ss = SearchNode(state=start_state, parent=None, cost=0, \
						heuristic=self.heuristic(start_state, self.goal_state))

		self.neighbor_count = 0

		# a global
		self.tiebreaker = count()
		heapq.heappush(self.frontier, (ss.heuristic, next(self.tiebreaker), ss))

	def search(self, time_limit):
		if self.goal_state == None:
			return None

		start_time = time.time()
		# extend nodes until the time limit is reached
		while time.time() - start_time < time_limit:
			if len(self.frontier) == 0:
				print("Search failed, bailing early")
				return
			self._step()

	def _step(self, make_tree=False, neighbor_callback=None):
		""" Perform one iteration of heuristic search - extend a single node.
			Requires that self.frontier is not empty.
		"""
		# Pop the best path to split.
		parent_score, parent_count, parent_state = heapq.heappop(self.frontier)
		if self.make_tree:
			parent_tree_node = parent_state.tree_node
		# print "Expand node:", parent_state

		if neighbor_callback:
			neighbor_states = []

		for neighbor_state in self.neighbors(parent_state.state):
			# prune any path segments that are bound to fail
			if not self.should_bail(neighbor_state, self.goal_state) and self.is_admissible(neighbor_state):
				# print "   neighbor: ", neighbor_state
				# build the tree representation of the search
				if self.make_tree:
					ntn = TreeNode(state=neighbor_state, children=[])
					parent_tree_node.children.append(ntn)

					nss = SearchNodeTree(state=neighbor_state, parent=parent_state, tree_node=ntn, \
						cost=parent_state.cost + self.cost(neighbor_state),  \
						heuristic=self.heuristic(neighbor_state, self.goal_state))
				else:
					nss = SearchNode(state=neighbor_state, parent=parent_state, \
						cost=parent_state.cost + self.cost(neighbor_state),  \
						heuristic=self.heuristic(neighbor_state, self.goal_state))

				if neighbor_callback:
					neighbor_states.append(neighbor_state)

				self.neighbor_count += 1

				score = nss.cost + nss.heuristic
				if self.goal_met(neighbor_state, self.goal_state):
					# print "Met goal:", neighbor_state
					heapq.heappush(self.found_paths, (score, next(self.tiebreaker), nss))
				else:
					# print (score, len(self.frontier))
					heapq.heappush(self.frontier, (score, next(self.tiebreaker), nss))

		if neighbor_callback:
			neighbor_callback(neighbor_states)

		# self.map.add_circle_to_exploration_buffer(parent_state.state)
		self.closed_set.append(parent_state.state)
		self.step_count += 1

	def best(self):
		if len(self.found_paths) > 0:
			return self.make_path(self.found_paths[0][2], True, False)
		return None

	def make_path(self, end_node, add_goal=False, remove_start=True):
		# a path is a list of control states, in the order of traversal
		# builds it in reverse order by following parent pointers, then reverses result
		path = [end_node.state]
		while end_node.parent:
			path.append(end_node.parent.state)
			end_node = end_node.parent
		# don't execute the start state
		if remove_start:
			path.pop()
		path.reverse()

		if add_goal:
			path.append(self.goal_state)
		return Path(states=path)

	""" The following functions should be overridden by the parent class
	"""
	def cost(self, path):
		# estimate the cost of traversing a given path
		raise NotImplementedError("HeuristicSearch cost not specified")

	def heuristic(self, state, goal_state):
		# return an estimate for cost to go between the given state and the goal state
		raise NotImplementedError("HeuristicSearch heuristic not specified")

	def goal(self):
		# return the goal state
		raise NotImplementedError("HeuristicSearch goal not specified")

	def goal_met(self):
		# return true if the goal state is met
		raise NotImplementedError("HeuristicSearch goal not specified")

	def should_bail(self, state, goal_state):
		# return true if the node should be terminated
		return False

	def is_admissible(self, state):
		# return true if the given state is considered admissible - fulfills hard constraints
		raise NotImplementedError("HeuristicSearch is_admissible not specified")

	def neighbors(self, state):
		# return a discretized set of neighbors of the given state
		raise NotImplementedError("HeuristicSearch neighbors not specified")

class SpaceExploration(HeuristicSearch):
	""" This class implements lower dimensional search to provide a high quality
		heuristic for higher dimensional path planning. If this search fails to
		provide a solution, then there is a good chance that no solution exists
		and the car should just stop

		# TODO:
			put expanded nodes into a KD tree for fast nearest neighbors search.
				This can prevent getting stuck in dead ends

		See:
		    http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=7353741
		    http://www6.in.tum.de/Main/Publications/ChenChao2013a.pdf
	"""
	def __init__(self, omap):
		# fetch relevant parameters
		self.branch_factor    = int(rospy.get_param("~branch_factor"))
		self.min_turn_radius  = float(rospy.get_param("~minimum_turn_radius"))
		self.soft_min_radius  = float(rospy.get_param("~soft_min_radius"))
		self.soft_min_penalty = float(rospy.get_param("~soft_min_penalty"))
		self.hard_min_radius  = float(rospy.get_param("~hard_min_radius"))
		self.heuristic_bias   = float(rospy.get_param("~heuristic_bias"))
		self.min_goal_overlap = float(rospy.get_param("~min_goal_overlap"))
		self.half_space_theta = float(rospy.get_param("~half_space_theta", 0))

		# cache reused values
		self.map = omap
		self.thetas = np.zeros(self.branch_factor, dtype=float)
		self.radii  = np.zeros(self.branch_factor)
		self.euclidean_neighbors = np.zeros((self.branch_factor,2),dtype=np.float64)
		self.next_goal = None

		# initialize super
		super(SpaceExploration, self).__init__(make_tree=False)

	def add_states_to_exploration_buffer(self, neighbor_circles):
		if len(neighbor_circles):
			poses = np.array(map(lambda c: c.center, neighbor_circles))
			radii = np.array(map(lambda c: c.radius, neighbor_circles))
			self.map.add_circles_to_exploration_buffer(poses, radii)

	def reset(self, pose):
		# Frontier is a priority queue.
		start_state = Circle(center=pose[:2], radius=self.map.get_distances(np.array([pose]))[0], angle=pose[2], deflection=0)
		self.map.clear_exploration_buffer()
		super(SpaceExploration, self).reset(start_state)

	def cost(self, state):
		if state.radius < self.soft_min_radius:
			return state.radius * self.soft_min_penalty
		return state.radius

	def heuristic(self, state, goal_state):
		if goal_state == None:
			return np.inf
		return (np.linalg.norm(state.center - goal_state.center) - state.radius)*self.heuristic_bias

	def goal(self):
		return self.next_goal

	def set_goal(self, pose):
		r = self.map.get_distances(np.array([pose]))[0]
		self.next_goal = Circle(center=pose[:2], radius=r, angle=pose[2], deflection=0)

	def overlap(self, s1, s2, percentage=.15):
		# TODO/NOTE: this is a bit of a hack to go fast, percentage overlap not accurate
		if np.linalg.norm(s1.center - s2.center) > (s1.radius + s2.radius)*(1.0-percentage):
			return False
		else:
			return True

	def goal_met(self, state, goal_state):
		return self.overlap(state, goal_state, self.min_goal_overlap)

	def is_admissible(self, state):
		return state.radius > self.hard_min_radius

	def search(self, time_limit):
		# same as before, but this returns early if paths are found
		start_time = time.time()

		# extend nodes until the time limit is reached
		while time.time() - start_time < time_limit and not len(self.found_paths):
			if len(self.frontier) == 0:
				print("Search failed, bailing early")
				return
			self._step(neighbor_callback=self.add_states_to_exploration_buffer)

	def neighbors(self, state):
		max_angle = utils.max_angle(self.min_turn_radius, state.radius)
		self.thetas = np.linspace(-max_angle,max_angle, num=self.branch_factor)
		self.euclidean_neighbors[:,0] = state.radius * np.cos(self.thetas + state.angle) + state.center[0]
		self.euclidean_neighbors[:,1] = state.radius * np.sin(self.thetas + state.angle) + state.center[1]

		# we use this several times, so compute this icky math here to avoid unnecessary computation
		euc = self.euclidean_neighbors.copy()
		utils.world_to_map(euc, self.map.map_info)
		euc = np.round(euc).astype(int)
		radii = self.map.get_distances(euc, coord_convert=False, check_bounds=True)
		mask = np.logical_and(np.invert(self.map.get_explored(euc, coord_convert=False)), self.map.get_permissible(euc, coord_convert=False))

		neighbors = map(lambda c,r,t: Circle(center=c.copy(), radius=r, angle=state.angle + t, deflection=t),
			self.euclidean_neighbors[mask,:], radii[mask], self.thetas[mask])
		return neighbors

class PathPlanner(HeuristicSearch):
	''' This class is used to refine an initial trajectory from the space
	    exploration search. It uses the original path as a heuristic and
	    integrates a term in the cost function which penalizes excessive
	    turning. Generates suboptimal but visually smooth paths.
	'''
	def __init__(self, omap):
		# fetch relevant parameters
		self.branch_factor    = int(rospy.get_param("~fp_branch_factor"))
		self.min_turn_radius  = float(rospy.get_param("~fp_minimum_turn_radius"))
		self.soft_min_radius  = float(rospy.get_param("~soft_min_radius"))
		self.soft_min_penalty = float(rospy.get_param("~fp_soft_min_penalty"))
		self.hard_min_radius  = float(rospy.get_param("~hard_min_radius"))
		self.heuristic_bias   = float(rospy.get_param("~fp_heuristic_bias"))
		self.min_goal_overlap = float(rospy.get_param("~min_goal_overlap"))
		self.half_space_theta = float(rospy.get_param("~half_space_theta", 0))
		self.exploration_coeff= float(rospy.get_param("~fp_exploration_coeff"))
		self.max_circle_radius= float(rospy.get_param("~fp_max_circle_radius"))

		self.map = omap
		self.next_goal = None
		self.rough_trajectory = None

		# initialize super
		super(PathPlanner, self).__init__(make_tree=False)

	def set_heursitic_trajectory(self, trajectory):
		self.rough_trajectory = trajectory

	def reset(self, pose):
		start_state = Circle(center=pose[:2], radius=self.map.get_distances(np.array([pose]))[0], angle=pose[2], deflection=0)
		self.map.clear_exploration_buffer()
		super(PathPlanner, self).reset(start_state)

	def cost(self, state):
		''' This is a fairly ad hoc cost function which tries to minimize path length and
		    curvy paths.
		'''
		other_coeff = 0.0
		if state.radius < self.soft_min_radius:
			other_coeff += 0.5
			# print "softmax"
			# return state.radius * self.soft_min_penalty + 10.0*np.power(abs(state.deflection),10.0)
		# return state.radius
		# return np.power(state.radius, 0.75) + 1.0*np.power(abs(state.deflection),3.0)

		# if state.deflection == 0.0:
		# 	other_coeff = -0.15
		# return state.radius + 5.0*np.power(abs(state.deflection),2.0) + other_coeff

		# if state.deflection == 0.0:
		# 	other_coeff = -0.08
		# return state.radius + 3.0*np.power(abs(state.deflection),1.2) + other_coeff

		# if state.deflection == 0.0:
		# 	other_coeff = -0.02
		return state.radius + 1.5*np.power(abs(state.deflection),1.) + other_coeff

		# if state.deflection == 0.0:
		# 	other_coeff = -0.1
		# return state.radius + 3.0*np.power(abs(state.deflection),1.) + other_coeff
		# return state.radius + 1.8*np.power(abs(state.deflection),1.0) + other_coeff

	def heuristic(self, state, goal_state):
		''' Finds the nearest point along the rough trajectory and then uses path length from there to the end as heursitic
		'''
		nearest_point, nearest_dist, t, i = utils.nearest_point_on_trajectory(state.center, self.rough_trajectory.np_points)
		distance_to_goal_along_trajectory = self.rough_trajectory.distance_to_end(i+t)
		return (distance_to_goal_along_trajectory - state.radius*2.0) * self.heuristic_bias

	def goal(self):
		return self.next_goal

	def set_goal(self, pose):
		r = self.map.get_distances(np.array([pose]))[0]
		r = np.clip(r, 0.0, self.max_circle_radius)
		self.next_goal = Circle(center=pose[:2], radius=r, angle=pose[2], deflection=0)

	def overlap(self, s1, s2, percentage=.15):
		# TODO/NOTE: this is a bit of a hack to go fast, percentage overlap not accurate
		if np.linalg.norm(s1.center - s2.center) > (s1.radius + s2.radius)*(1.0-percentage):
			return False
		else:
			return True

	def goal_met(self, state, goal_state):
		return self.overlap(state, goal_state, self.min_goal_overlap)

	def is_admissible(self, state):
		return state.radius > self.hard_min_radius

	def search(self, time_limit):
		# same as before, but this returns early if paths are found
		start_time = time.time()

		# extend nodes until the time limit is reached
		while time.time() - start_time < time_limit and not len(self.found_paths):
			if len(self.frontier) == 0:
				print("Search failed, bailing early")
				return
			self._step(neighbor_callback=self.add_states_to_exploration_buffer)

	def add_states_to_exploration_buffer(self, neighbor_circles):
		if len(neighbor_circles):
			poses = np.array(map(lambda c: c.center, neighbor_circles))
			radii = np.array(map(lambda c: c.radius, neighbor_circles))
			self.map.add_circles_to_exploration_buffer(poses, radii, exp_coeff=self.exploration_coeff)

	def neighbors(self, state):
		''' This function determines the set of points along the perimeter of the given
		    circle state that are accessible according to ackermann geometry and steering
		    angle limitations. Then it samples the permissible theta space to generate neighbors.

		    TODO: computation of neighboring thetas is not accurate - it should actually be
		          the direction tangent to the path arc at the intersection with the circle radius.
		          For now this works, and is a strict underestimate of the reachable thetas
		'''
		max_angle = utils.max_angle(self.min_turn_radius, state.radius)
		percentage = np.power(2.0 * self.branch_factor / np.pi,0.9)
		actual_branch_factor = 2*int(max_angle * percentage / 2.0)+1
		thetas = np.linspace(-max_angle,max_angle, num=actual_branch_factor)

		euclidean_neighbors = np.zeros((actual_branch_factor,2))
		euclidean_neighbors[:,0] = state.radius * np.cos(thetas + state.angle) + state.center[0]
		euclidean_neighbors[:,1] = state.radius * np.sin(thetas + state.angle) + state.center[1]

		# perform coordinate space conversion here and then index into the exploration and permissible
		# buffers to prune states which are not permissible or already explored
		euc = euclidean_neighbors.copy()
		utils.world_to_map(euc, self.map.map_info)
		euc = np.round(euc).astype(int)
		radii = self.map.get_distances(euc, coord_convert=False, check_bounds=True) - 0.05
		radii = np.clip(radii, 0.0, self.max_circle_radius)
		mask = np.logical_and(np.invert(self.map.get_explored(euc, coord_convert=False)), self.map.get_permissible(euc, coord_convert=False))

		# generate a neighbor state for each of the permissible
		neighbors = map(lambda c,r,t: Circle(center=c.copy(), radius=r, angle=state.angle + t, deflection=t),
			euclidean_neighbors[mask,:], radii[mask], thetas[mask])
		return neighbors

class FindTrajectory(object):
	""" Performs two levels of heuristic search to plan paths between the robot pose and a clicked goal pose.
	"""
	def __init__(self):
		self.found_trajectory  = utils.LineTrajectory("/found_trajectory")
		self.rough_trajectory  = utils.LineTrajectory("/rough_trajectory")
		self.should_publish    = bool(rospy.get_param("~publish"))
		self.pub_topic         = rospy.get_param("~trajectory_topic")
		self.odom_topic        = rospy.get_param("~odom_topic")
		self.exploration_timeout  = rospy.get_param("~exploration_timeout")
		self.show_exploration_buffer = rospy.get_param("~show_exploration_buffer")
		self.save_trajectory   = rospy.get_param("~save_trajectory")
		self.save_path         = rospy.get_param("~save_path")

		self.should_refine_trajectory = bool(rospy.get_param("~refine_trajectory"))
		self.refining_timeout  = rospy.get_param("~refining_timeout")

		self.map = None
		self.map_initialized = False
		self.last_pose_time = 0.0
		self.circle_path = None

		# just a few initial poses for debugging purposes
		# self.last_pose = np.array([-1., 0., -np.pi/1.3])
		# self.last_pose = np.array([-1., 0., 0.0])
		self.last_pose = np.array([-3., -1., np.pi])
		# self.last_pose = np.array([-1.27, -2.6, np.pi])
		# self.last_pose = np.array([-6.27, -2.6, np.pi])
		# self.last_pose = np.array([15.3, 25.18, 0])
		# self.last_pose = np.array([-10,  33.8, 0])
		# self.last_pose = None

		self.get_omap()
		self.space_explorer = SpaceExploration(self.map)

		if self.should_publish:
			self.traj_pub = rospy.Publisher(self.pub_topic, PolygonStamped, queue_size=1)

		if self.should_refine_trajectory:
			self.path_planner = PathPlanner(self.map)
			self.fast_path = None
			self.fast_trajectory = utils.LineTrajectory("/fast_trajectory")

		# visualization publishers
		self.viz_namespace = "/circle_search"
		self.circle_pub = rospy.Publisher(self.viz_namespace + "/exploration_circles", MarkerArray, queue_size=1)
		self.fast_circle_pub = rospy.Publisher(self.viz_namespace + "/fast_circles", MarkerArray, queue_size=1)
		self.exploration_pub = rospy.Publisher(self.viz_namespace + "/exploration_buffer", OccupancyGrid, queue_size=1)

		self.nav_goal_sub = rospy.Subscriber("/move_base_simple/goal", PoseStamped, self.goal_point_callback, queue_size=1)
		self.odom_sum = rospy.Subscriber(self.odom_topic,  Odometry, self.odom_callback, queue_size=1)

		print("Initialized. Waiting on messages...")

	def get_omap(self):
		map_service_name = rospy.get_param("~static_map", "static_map")
		print("Fetching map from service: ", map_service_name)
		rospy.wait_for_service(map_service_name)
		map_msg = rospy.ServiceProxy(map_service_name, GetMap)().map
		self.map = utils.Map(map_msg)
		self.map_initialized = True
		print("Finished loading map")

	def publish_trajectory(self):
		print("Publishing trajectory to:", self.pub_topic)
		self.traj_pub.publish(self.rough_trajectory.toPolygon())

	def non_overlapping_paths(self, node):
		# given a tree data structure, this function will return list of lists of the
		# node "state" attributes
		# if visualized as a tree, these states will not contain overlapping segments
		if not node.children:
			return [node.state]
		else:
			paths = []

			for child in node.children:
				child_paths = self.non_overlapping_paths(child)

				if type(child_paths[0]) == list:
					# child is not a leaf node, add self to first path and store
					child_paths[0].insert(0, node.state)
					paths = paths + child_paths
				else:
					# this is a leaf node, add self to path and store
					child_paths.insert(0, node.state)
					paths.append(child_paths)

			return paths

	def visualize(self):
		if self.circle_pub.get_num_connections() > 0 and self.circle_path:
			print("CIRCLE PUB")
			# root = self.space_explorer.tree_root
			# markers = [utils.marker_clear_all("/map")]
			# explored = self.non_overlapping_paths(root)
			# if type(explored[0]) == list:
			# 	explored = reduce(lambda x,y: x+y, explored)
			# print "explored: ", len(explored)
			# markers += [utils.marker_from_circle(circle, index=i, linewidth=0.01, color=ColorRGBA(0, 1, 0, 0.1), lifetime=1.0) \
			#             for i, circle in enumerate(explored)]

			# marker_array = MarkerArray(markers=markers)
			# self.circle_pub.publish(marker_array)
			# print len(self.circle_path.states)
			# marker = utils.marker_from_circles(self.circle_path.states, 1, "Markers_NS2", "/map", 3.0, [1.,0.,0.])
			# self.circle_pub.publish(marker)

			markers = [utils.marker_clear_all("/map")]
			markers += [utils.marker_from_circle(circle, index=i, linewidth=0.05, color=ColorRGBA(1, 0, 0, 0.4), \
						lifetime=4.0)
						for i, circle in enumerate(self.circle_path.states)]
			# print len(markers)
			marker_array = MarkerArray(markers=markers)
			self.circle_pub.publish(marker_array)

		if self.should_refine_trajectory and self.fast_circle_pub.get_num_connections() > 0 and self.fast_path:
			print("FAST CIRCLE PUB")
			# root = self.path_planner.tree_root
			# markers = [utils.marker_clear_all("/map")]
			# explored = self.non_overlapping_paths(root)
			# if type(explored[0]) == list:
			# 	explored = reduce(lambda x,y: x+y, explored)
			# print "explored: ", len(explored)
			# markers += [utils.marker_from_circle(circle, index=i, linewidth=0.01, color=ColorRGBA(0, 0, 1, 0.1), lifetime=1.0) \
			#             for i, circle in enumerate(explored)]
			markers = [utils.marker_clear_all("/map")]
			markers += [utils.marker_from_circle(circle, index=i, linewidth=0.05, color=ColorRGBA(1, 0, 0, 0.4), \
						lifetime=4.0)
						for i, circle in enumerate(self.fast_path.states)]

			marker_array = MarkerArray(markers=markers)
			self.fast_circle_pub.publish(marker_array)

		if self.show_exploration_buffer and self.exploration_pub.get_num_connections() > 0:
			print("Publishing exploration grid")
			self.exploration_pub.publish(self.map.get_exploration_occupancy_grid())
		# plt.imshow(self.map.exploration_buffer)
		# plt.imshow(self.map.permissible_region)
		# plt.show()
		print("visualize")

	def goal_point_callback(self, msg):
		''' Initialize path search. First performs and exploration search, and then a second search which uses
		    the result from the first search as a heuristic. Roughly respects dynamic constraints of the car.
		'''
		self.goal = np.array([msg.pose.position.x, msg.pose.position.y, utils.quaternion_to_angle(msg.pose.orientation)])
		print("\nNew goal:", self.goal)

		if self.has_recent_pose():
			# perform the first search
			self.find_rough_trajectory()
			if self.circle_path:
				self.rough_trajectory.clear()
				self.rough_trajectory.points = map(lambda x: x.center, self.circle_path.states)
				self.rough_trajectory.update_distances()
				self.found_trajectory.points = self.rough_trajectory.points
				self.found_trajectory.publish_viz()

				if self.should_refine_trajectory:
					self.refine_trajectory()
					if self.fast_path:
						self.fast_trajectory.clear()
						self.fast_trajectory.points = map(lambda x: x.center, self.fast_path.states)
						self.fast_trajectory.update_distances()
						self.found_trajectory.points = self.fast_trajectory.points

				if self.publish_trajectory:
					self.traj_pub.publish(self.found_trajectory.toPolygon())

				if self.save_trajectory:
					self.found_trajectory.save(os.path.join(self.save_path, time.strftime("%Y-%m-%d-%H-%M-%S") + ".traj"))

				self.found_trajectory.update_distances()
				self.found_trajectory.publish_viz()
			self.visualize()
		else:
			print("No recent odometry, skipping search!")

	def find_rough_trajectory(self):
		print("Finding rough trajectory")
		self.space_explorer.set_goal(self.goal)
		self.space_explorer.reset(self.last_pose.copy())
		t = time.time()
		self.space_explorer.search(self.exploration_timeout)
		print("... rough search terminated in:", time.time() - t, "seconds")
		print("... expanded", self.space_explorer.step_count)
		print("... considered", self.space_explorer.neighbor_count)
		self.circle_path = self.space_explorer.best()
		if self.circle_path == None:
			print("...search failed")

	def refine_trajectory(self):
		print("Refining trajectory")
		self.rough_trajectory.make_np_array()
		self.path_planner.set_heursitic_trajectory(self.rough_trajectory)
		self.path_planner.set_goal(self.goal)
		self.path_planner.reset(self.last_pose.copy())

		t = time.time()
		self.path_planner.search(self.refining_timeout)
		print("... refine search terminated in:", time.time() - t, "seconds")
		print("... expanded", self.path_planner.step_count)
		print("... considered", self.path_planner.neighbor_count)
		self.fast_path = self.path_planner.best()
		if self.fast_path == None:
			print("... search failed")

	def has_recent_pose(self):
		# return True # this is useful for debugging
		return time.time() - self.last_pose_time < 0.5

	def odom_callback(self, msg):
		if not self.has_recent_pose():
			print("Recieved Odometry...")
		self.last_pose = np.array([msg.pose.pose.position.x, msg.pose.pose.position.y, utils.quaternion_to_angle(msg.pose.pose.orientation)])
		self.last_pose_time = time.time()

def make_flamegraph(filterx=None):
    import flamegraph
    perf_log_path = os.path.join(os.path.dirname(__file__), "../tmp/perf2.log")
    flamegraph.start_profile_thread(fd=open(perf_log_path, "w"),
                                    filter=filterx,
                                    interval=0.001)

if __name__=="__main__":
	rospy.init_node("trajectory_search")
	# make_flamegraph(r"goal_point_callback")
	pf = FindTrajectory()
	rospy.spin()