text:
import numpy as np

def mask_center_2d(center_x, center_y, r, x, y):
    """

    :param center: 2D coordinate of center position of circular mask
    :param r: radius of mask in pixel values
    :param data: data image
    :return:
    """
    x_shift = x - center_x
    y_shift = y - center_y
    R = np.sqrt(x_shift*x_shift + y_shift*y_shift)
    mask = np.empty_like(R)
    mask[R > r] = 1
    mask[R <= r] = 0
    n = int(np.sqrt(len(x)))
    mask_2d = mask.reshape(n, n)
    return mask_2d

def mask_sphere(x, y, center_x, center_y, r):
    """

    :param center: 2D coordinate of center position of circular mask
    :param r: radius of mask in pixel values
    :param data: data image
    :return:
    """
    x_shift = x - center_x
    y_shift = y - center_y
    R = np.sqrt(x_shift*x_shift + y_shift*y_shift)
    mask = np.empty_like(R)
    mask[R > r] = 0
    mask[R <= r] = 1
    return mask

def mask_half_moon(x, y, center_x, center_y, r_in, r_out, phi0=0, delta_phi=2*np.pi):
    """

    :param x:
    :param y:
    :param center_x:
    :param center_y:
    :param r_in:
    :param r_out:
    :param phi:
    :param delta_phi:
    :return:
    """
    x_shift = x - center_x
    y_shift = y - center_y
    R = np.sqrt(x_shift*x_shift + y_shift*y_shift)
    phi = np.arctan2(x_shift, y_shift)
    #phi = np.abs(phi)
    phi_min = phi0 - delta_phi/2.
    phi_max = phi0 + delta_phi/2.
    mask = np.zeros_like(x)
    if phi_max > phi_min:
        mask[(R < r_out) & (R > r_in) & (phi > phi_min) & (phi < phi_max)] = 1
    else:
        mask[(R < r_out) & (R > r_in) & (phi > phi_max)] = 1
        mask[(R < r_out) & (R > r_in) & (phi < phi_min)] = 1
    return mask

max_stars_repo_name: guoxiaowhu/lenstronomy

max_issues_repo_name: guoxiaowhu/lenstronomy

max_forks_repo_name: guoxiaowhu/lenstronomy


text:
import torch
from torch.autograd import Variable
import matplotlib.pyplot as plt
import numpy as np
import shutil
import pickle
import json
import os
import pdb
from lib.coarsening import perm_data_torch2
import visdom
import collections
from tqdm import tqdm

if torch.cuda.is_available():
    print('cuda available')
    dtypeFloat = torch.cuda.FloatTensor
    dtypeLong = torch.cuda.LongTensor
    torch.cuda.manual_seed(1)
else:
    print('cuda not available')
    dtypeFloat = torch.FloatTensor
    dtypeLong = torch.LongTensor
    torch.manual_seed(1)

def save_checkpoint(state, is_best, filename='checkpoint.pth.tar'):
    torch.save(state, filename)
    if is_best:
        shutil.copyfile(filename, 'model_best.pth.tar')

class ecg_trainer:
    def __init__(self, nn_model, train_loader, test_loader):
        self.nn_model = nn_model
        self.train_loader = train_loader
        self.test_loader = test_loader

    def train_model(self, optimizer, global_lr, L, lmax, num_epoch=15, l2_reg=5e-4, plt_fig=False):
        running_loss = 0.0
        running_accuray = 0
        running_total = 0
        lr = global_lr
        for epoch in range(num_epoch):
            running_loss = 0
            tqdm_len = len(self.train_loader)
            with tqdm(total=tqdm_len) as pbar:
                for ind, sample in enumerate(self.train_loader):
                    instance = Variable(sample['sig'].cuda())
                    label = sample['label']
                    label = torch.LongTensor(label).type(dtypeLong)
                    label = Variable(label, requires_grad=False)
                    # label = Variable(sample['label'])
                    # self.optimizer.zero_grad()
                    # pdb.set_trace()
                    y_pred = self.nn_model.forward(instance, 0, L, lmax)
                    # y_pred = y_pred.view(-1, 2)
                    # pdb.set_trace()
                    loss = self.nn_model.loss(y_pred, label, l2_reg)
                    loss_train = loss.data[0]
                    # print(loss.data[0])
                    acc_train = self.nn_model.evaluation(y_pred, label.data)

                    loss.backward()

                    optimizer.step()
                    optimizer.zero_grad()

                    # loss, accuracy
                    running_loss += loss_train
                    running_accuray += acc_train
                    running_total += 1
                    pbar.update(len(sample['label']))
                    if not running_total % 100:  # print every x mini-batches
                        print('epoch= %d, i= %4d, loss(batch)= %.4f, accuray(batch)= %.2f' %
                              (epoch+1, running_total, loss_train, acc_train))
                print('epoch= %d, loss(train)= %.3f, accuracy(train)= %.3f, lr= %.5f' %
                      (epoch+1, running_loss/running_total, running_accuray/running_total, lr))

            lr = global_lr / (epoch+1)
            optimizer = self.nn_model.update_learning_rate(optimizer, lr)

            running_accuray_test = 0
            running_total_test = 0

            for ind, sample in enumerate(self.test_loader):
                test_x = Variable(torch.FloatTensor(sample[0]).type(dtypeFloat),
                                  requires_grad=False)
                y_pred = self.nn_model.forward(test_x, 0.0, L, lmax)
                test_y = sample[1]
                test_y = torch.LongTensor(test_y).type(dtypeLong)
                test_y = Variable(test_y, requires_grad=False)
                acc_test = self.nn_model.evaluation(y_pred, test_y.data)
                running_accuray_test += acc_test
                running_total_test += 1

            print('  accuracy(test) = %.3f %%' % (running_accuray_test / running_total_test))

class simple_trainer:
    def __init__(self, nn_model, graph_model, train_loader=None,
                 train_loader2=None, train_loader_graph=None, test_loader=None, test_loader2=None,
                 loss_fn=None, optimizer='adam', graph_optimizer='adam'):
        self.nn_model = nn_model
        self.graph_model = graph_model
        # self.nn_model.cuda()
        self.train_loader = train_loader
        self.train_loader2 = train_loader2
        self.train_loader_graph = train_loader_graph
        self.test_loader = test_loader
        self.test_loader2 = test_loader2
        self.loss_fn = loss_fn
        self.start_epoch = 0
        # self.dropout_value = 0.5
        if optimizer == 'adam':
            self.optimizer = torch.optim.Adam(self.nn_model.parameters())
        if graph_optimizer == 'adam':
            self.graph_optimizer = torch.optim.Adam(self.graph_model.parameters())
        # self.valid_acc_list = list()

    def load_model(self, load_path='model_best.pth.tar'):
        if os.path.isfile(load_path):
            print("=> loading checkpoint '{}'".format(load_path))
            checkpoint = torch.load(load_path)
            self.start_epoch = checkpoint['epoch']
            # best_prec1 = checkpoint['best_prec1']
            self.nn_model.load_state_dict(checkpoint['state_dict'])
            self.optimizer.load_state_dict(checkpoint['optimizer'])
            print("=> loaded checkpoint '{}' (epoch {})"
                  .format(load_path, checkpoint['epoch']))
        else:
            print("=> no checkpoint found at '{}'".format(load_path))

    def train_model(self, num_epoch=15, plt_fig=False):
        print('TrainSet :', len(self.train_loader))
        # self.valid_acc_list = np.arange(15)
        # plt.axis([0, num_epoch, 0, 1])
        # plt.ion()
        epoch_list = 0
        val_acc = 0
        if plt_fig:
            fig = plt.figure(1)
            ax = fig.add_subplot(111)
            ax.set_xlim(0, num_epoch)
            ax.set_ylim(0, 1)
            line, = ax.plot(epoch_list, val_acc, 'ko-')
        curr_acc = 0
        best_acc = 0
        for epoch in range(num_epoch):
            running_loss = 0
            num_tr_iter = 0
            for ind, sample in enumerate(self.train_loader):
                instance = Variable(sample['sig'].cuda())
                label = Variable(sample['label'].cuda())
                # label = Variable(sample['label'])
                self.optimizer.zero_grad()
                # pdb.set_trace()
                y_pred, pred_layer_outs = self.nn_model(instance)
                # pdb.set_trace()
                y_pred = y_pred.view(-1, 2)
                # pdb.set_trace()
                loss = self.loss_fn(y_pred, label)
                # print(loss.data[0])

                loss.backward()
                self.optimizer.step()
                running_loss += loss.data[0]
                # print(self.nn_model.parameters())
                # if ind % 100 == 0:
                # if True:
                # print(epoch, running_loss/num_tr_points)
                num_tr_iter += 1
            self.curr_epoch = self.start_epoch + epoch
            print('epoch', self.curr_epoch, running_loss/num_tr_iter)
            # pdb.set_trace()
            if plt_fig:
                epoch_list = np.concatenate((line.get_xdata(), [epoch]))
                val_acc = np.concatenate((line.get_ydata(), [self.test_model()]))
                # plt.plot(epoch, val_acc, '.r-')
                line.set_data(epoch_list, val_acc)
                plt.pause(0.01)

            # self.valid_acc_list.append(val_acc)
            else:               #
                curr_acc = self.test_model()
                is_best = False
                # if curr_acc > best_acc:
                if False:
                    best_acc = curr_acc
                    is_best = True
                    save_checkpoint({
                        'epoch': self.curr_epoch + 1,
                        # 'arch': args.arch,
                        'state_dict': self.nn_model.state_dict(),
                        # 'best_prec1': best_prec1,
                        'optimizer': self.optimizer.state_dict(),
                    }, is_best)
        if plt_fig:
            return fig
        else:
            return

    # def get_nn_features(self, instance):
    #     # for ind, sample in enumerate(self.train_loader):
    #         # instance = Variable(sample['sig'].cuda())
    #     y_pred, pred_layer_outs = self.nn_model(instance)
    #     return pred_layer_outs

    # def retesting_model(self):
    #     print('TestSet :', len(self.test_loader))
    #     num_corr = 0
    #     tot_num = 0
    #     for sample in self.test_loader:
    #         instance = Variable(sample['sig'].cuda())

    def test_model(self):
        print('TestSet :', len(self.test_loader))
        num_corr = 0
        tot_num = 0
        # self.nn_model.eval()
        for sample in self.test_loader:
            instance = Variable(sample['sig'].cuda())
            y_pred, pred_layer_outs = self.nn_model(instance)
            # pdb.set_trace()
            y_pred = y_pred.view(-1, 2)
            _, label_pred = torch.max(y_pred.data, 1)
            # if (label_pred == sample['label'].cuda()).cpu().numpy():
            # if (label_pred.cpu() == sample['label']).numpy():

            num_corr += (label_pred.cpu() == sample['label']).any()
            # tot_num += label_pred.shape[0]
            tot_num += 1
        print(num_corr, tot_num, num_corr/tot_num)
        # self.nn_model.train()
        return num_corr/tot_num

    def cnn_features_save(self, train_bool=False, fname='../data/cnn_features.pkl'):
        out_train_list = dict()
        sample_count = 0
        loader = self.test_loader2
        if train_bool:
            loader = self.train_loader2

        for ind, sample in enumerate(loader):
            # for ind, sample in enumerate(self.test_loader2):
            out_train = dict()
            instance = Variable(sample['sig'].cuda())
            # label = Variable(sample['label'].cuda())
            idx = Variable(sample['idx'])
            pidx = Variable(sample['pidx'])
            y_pred, channel_layer_outs = self.nn_model(instance)
            out_train['idx'] = idx
            out_train['pidx'] = pidx
            out_train['label'] = sample['label']
            out_train['channel_layer_outs'] = channel_layer_outs
            # out_train_list.append(out_train)
            # out_train_list['index'] = sample_count
            out_train_list[sample_count] = out_train
            sample_count += 1
            print('Ind', ind, 'done')
        # print(out_train_list[0])
        with open(fname, 'wb') as f:
            pickle.dump(out_train_list, f)

    def graph_nn_train(self, L, lmax, perm, l2_reg=5e-4, num_epoch=10, dr_val=0.5):
        running_loss = 0
        running_accuray = 0
        running_total = 0

        for epoch in range(num_epoch):
            # running_loss = 0
            # num_tr_iter = 0
            for ind, sample in enumerate(self.train_loader_graph):
                instance = Variable(sample['sig'].cuda())
                label = Variable(sample['label'].cuda(), requires_grad=False)
                # last_layer_features = self.get_nn_features(instance)
                y_pred, channel_layer_outs = self.nn_model(instance)
                inp_fin_outs1 = [channel_layer_outs[i]['fc1'] for i in range(6)]
                # pdb.set_trace()
                inp_fin_outs = torch.stack(inp_fin_outs1)
                inp_fin_outs = inp_fin_outs.detach()
                inp_fin_outs = perm_data_torch2(inp_fin_outs, perm)
                # inp_fin_outs.cuda()
                inp_fin_outs = inp_fin_outs.permute(1, 0, 2)
                # pdb.set_trace()
                self.graph_optimizer.zero_grad()
                y_pred_graph = self.graph_model(inp_fin_outs, dr_val, L, lmax)
                # y_pred_graph = y_pred.view(-1, 2)
                y_pred_graph = y_pred_graph.view(-1, 2)
                # pdb.set_trace()
                loss = self.graph_model.loss(y_pred_graph, label, l2_reg)
                loss_train = loss.data[0]

                acc_train = self.graph_model.evaluation(y_pred_graph, label.data)
                # pdb.set_trace()
                loss.backward()

                # global_step += self.batch_size
                self.graph_optimizer.step()

                running_loss += loss_train
                running_accuray += acc_train
                running_total += 1
                # print('Run Total', running_total)

            print('epoch= %d, loss(train)= %.3f, accuracy(train)= %.3f' %
                  (epoch+1, running_loss/running_total, running_accuray/running_total))

        # self.graph_model(channel_layer_outs, )
        # pdb.set_trace()

class end_to_end_trainer:
    def __init__(self, nn_model, train_loader=None, test_loader=None,
                 loss_fn=None, optimizer='adam', tovis=False):
        self.nn_model = nn_model
        if torch.cuda.is_available():
            self.nn_model.cuda()
        self.train_loader = train_loader
        self.test_loader = test_loader
        self.loss_fn = loss_fn
        self.start_epoch = 0
        # Note: may need to change the below line
        self.curr_epoch = 0
        if optimizer == 'adam':
            self.optimizer = torch.optim.Adam(self.nn_model.parameters())
        elif optimizer == 'sgd':
            self.optimizer = torch.optim.SGD(self.nn_model.parameters(), lr=0.17)
        self.tovis = tovis
        if self.tovis:
            self.vis = visdom.Visdom()
            self.loss_window = self.vis.line(X=torch.zeros((1)).cpu(), Y=torch.zeros((1)).cpu(),
                                             opts=dict(xlabel='minibatches', ylabel='Loss',
                                                       title='Training Loss', legend=['Loss']))

    def loss_l2(self, y_pred, label_gt):
        loss = self.loss_fn(y_pred, label_gt)
        for params in self.nn_model.parameters():
            data = params * params
            loss += data.sum()
        return loss

    def save_checkpoint1(state, is_best, filename='checkpoint_e2e.pth.tar'):
        # f = open(filename, 'w')
        torch.save(state, filename)
        if is_best:
            shutil.copyfile(filename, 'e2e_model_best.pth.tar')
        # f.close()

    def load_model(self, load_path='e2e_model_best.pth.tar'):
        if os.path.isfile(load_path):
            print("=> loading checkpoint '{}'".format(load_path))
            checkpoint = torch.load(load_path)
            self.start_epoch = checkpoint['epoch']
            # best_prec1 = checkpoint['best_prec1']
            self.nn_model.load_state_dict(checkpoint['state_dict'])
            self.optimizer.load_state_dict(checkpoint['optimizer'])
            print("=> loaded checkpoint '{}' (epoch {})"
                  .format(load_path, checkpoint['epoch']))
        else:
            print("=> no checkpoint found at '{}'".format(load_path))
        return

    def get_output(self, instance, d, L, lmax, perm, ep):
        return self.nn_model(instance, d, L, lmax, perm, ep)

    def proc_inp(self, inp):
        return inp

    def proc_pred(self, pred):
        return pred

    def proc_pred_lab(self, plab):
        return plab

    def train_model(self, d=None, L=None, lmax=None, perm=None, num_epoch=15):
        # pdb.set_trace()
        print('TrainSet: ', len(self.train_loader))
        self.tot_ind = len(self.train_loader)
        self.nn_model.train()
        curr_acc = 0
        best_acc = 0
        for epoch in range(num_epoch):
            running_loss = 0
            num_tr_iter = 0
            for ind, sample in enumerate(tqdm(self.train_loader)):
#                 instance = Variable(sample['sig'].cuda())
#                 label = Variable(sample['label'].type(dtypeLong))
                instance = Variable(sample[0])
                label = Variable(sample[1])
                instance = self.proc_inp(instance)
                # pdb.set_trace()
                self.optimizer.zero_grad()
                y_pred = self.get_output(instance, d, L, lmax, perm,
                                         self.start_epoch + epoch)
                y_pred = y_pred.view(-1, 2)
                y_pred = self.proc_pred(y_pred)
                # pdb.set_trace()
                loss = self.loss_fn(y_pred, label)
                loss.backward()
                self.optimizer.step()
                running_loss += loss.data[0]
                num_tr_iter += 1

                if self.tovis:
                    # pdb.set_trace()
                    self.vis.line(
                        X=torch.ones((1)).cpu() * (ind + (self.start_epoch + epoch)
                                                   * self.tot_ind),
                        Y=torch.Tensor([loss.data[0]]).cpu(),
                        win=self.loss_window,
                        update='append')

            self.curr_epoch = self.start_epoch + epoch
            print('epoch', self.curr_epoch, running_loss/num_tr_iter)
            curr_acc = self.test_model(d, L, lmax, perm)
            is_best = False

            if curr_acc > best_acc:
                best_acc = curr_acc
                is_best = True
                save_checkpoint({
                    'epoch': self.curr_epoch + 1,
                    # 'arch': args.arch,
                    'state_dict': self.nn_model.state_dict(),
                    # 'best_prec1': best_prec1,
                    'optimizer': self.optimizer.state_dict(),
                }, is_best, 'e2e_checkpoint.pth.tar')

    def test_model(self, d=None, L=None, lmax=None, perm=None):
        print('TestSet :', len(self.test_loader))
        num_corr = 0
        tot_num = 0
        false_0 = 0
        false_1 = 0
        true_0 = 0
        true_1 = 0
        pt_dis_pred_dict = {}
        pt_dis_actual_dict = {}
        self.nn_model.eval()
        for sample in tqdm(self.test_loader):
#             instance = Variable(sample['sig'].cuda())
            instance = Variable(sample[0]).cuda()
            instance = self.proc_inp(instance)

            y_pred = self.get_output(instance, d, L, lmax, perm, self.curr_epoch)
            y_pred = y_pred.view(-1, 2)
            _, label_pred = torch.max(y_pred.data, 1)
            label_pred = self.proc_pred_lab(label_pred)
            for ind1, sidx in enumerate(sample['idx']):
                if pt_dis_actual_dict.get(sidx) is None:
                    pt_dis_actual_dict[sidx] = sample['label'][ind1]
                    pt_dis_pred_dict[sidx] = [label_pred[ind1]]
                else:
                    assert pt_dis_actual_dict[sidx] == sample['label'][ind1]
                    pt_dis_pred_dict[sidx].append(label_pred[ind1])

        for k, v in pt_dis_actual_dict.items():
            pred_list = pt_dis_pred_dict[k]

            pred_list_counter = collections.Counter(pred_list)
            fin_label_pred, fin_label_count = pred_list_counter.most_common()[0]

            if fin_label_pred == v:
                num_corr += 1
            true_0 += (fin_label_pred == v) and (fin_label_pred == 0)
            true_1 += (fin_label_pred == v) and (fin_label_pred == 1)
            false_0 += (fin_label_pred != v) and (fin_label_pred == 0)
            false_1 += (fin_label_pred != v) and (fin_label_pred == 1)
            tot_num += 1
        self.nn_model.train()
        print(true_0, true_1, false_0, false_1)
        prec = true_1 / (false_1 + true_1)
        rec = true_1 / (true_1 + false_0)
        f1s = 2 * prec * rec / (prec + rec)

        print(prec, rec, f1s)
        print(num_corr, tot_num, num_corr/tot_num)
        return num_corr / tot_num

    def test_model2(self, d, L, lmax, perm):
        print('TestSet :', len(self.test_loader))
        num_corr = 0
        tot_num = 0
        self.nn_model.eval()
        for sample in self.test_loader:
            instance = Variable(sample['sig'].cuda())
            y_pred = self.nn_model(instance, d, L, lmax, perm, self.curr_epoch)
            y_pred = y_pred.view(-1, 2)
            _, label_pred = torch.max(y_pred.data, 1)

            num_corr += (label_pred.cpu() == sample['label']).any()

            # tot_num += label_pred.shape[0]
            tot_num += 1
        print(num_corr, tot_num, num_corr/tot_num)
        self.nn_model.train()
        return num_corr / tot_num

class ml_cnn_trainer(end_to_end_trainer):
    def __init__(self, config, train_loader, test_loader, nn_model,
                 loss_fn=None, optimizer='adam', tovis=False):
        self.config = config
        super(ml_cnn_trainer, self).__init__(nn_model, train_loader, test_loader,
                                             loss_fn, optimizer, tovis)

    def get_output(self, instance, d, L, lmax, perm, ep):
        # pdb.set_trace()
        return self.nn_model(instance)

class pe2e_trainer(ml_cnn_trainer):
    def get_output(self, instance, d, L, lmax, perm, ep):
        return self.nn_model(instance, L[0])

class ml_cnn_consensus_trainer(ml_cnn_trainer):
    def proc_inp(self, inp):
        self.B, self.cons, Din, num_c = inp.shape
        return inp.view(self.B*self.cons, Din, num_c)

    def proc_pred(self, pred):
        return pred.view(self.B, self.cons, 2)

    def proc_pred_lab(self, plab):
        label_pred_cons = plab.view(self.B, self.cons)
        label_pred, _ = label_pred_cons.mode(dim=1)
        return label_pred

ext: py

reason: import numpy


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Calculate image shifts between images.

Example
-------
To calculate the image shifts between the reference image (``ref_image``) and another image (``image``):

	>>> imk = ImageMotionKernel(ref_image=ref_image, warpmode='translation')
	>>> kernel = img.calc_kernel(image)
	>>> print(kernel)

.. codeauthor:: Mikkel N. Lund
.. codeauthor:: Rasmus Handberg <rasmush@phys.au.dk>
"""

import numpy as np
import logging
import cv2
#from skimage.transform import estimate_transform, warp, AffineTransform, EuclideanTransform
import math
from skimage.filters import scharr
from scipy.interpolate import interp1d
from astropy.io import fits
from astropy.wcs import WCS, NoConvergence
import warnings

class ImageMovementKernel(object):

	N_PARAMS = {
		'unchanged': 0,
		'translation': 2,
		'euclidian': 3,
		'affine': 6,
		'wcs': 1
	}

	#----------------------------------------------------------------------------------------------
	def __init__(self, warpmode='euclidian', image_ref=None, wcs_ref=None):
		"""
		Initialize ImageMovementKernel.

		Parameters:
			warpmode (str): Options are ``'wcs'``, ``'unchanged'``, ``'translation'``,
				``'euclidian'`` and ``'affine'``. Default is ``'euclidian'``.
			image_ref (2D ndarray): Reference image used.
			wcs_ref (:class:`astropy.wcs.WCS`): Reference WCS when using warpmode='wcs'.

		.. codeauthor:: Rasmus Handberg <rasmush@phys.au.dk>
		"""

		if warpmode not in ('wcs', 'unchanged', 'translation', 'euclidian', 'affine'):
			raise ValueError("Invalid warpmode")

		self.warpmode = warpmode
		self.image_ref = image_ref
		self.wcs_ref = wcs_ref
		self.n_params = ImageMovementKernel.N_PARAMS[self.warpmode]

		if self.image_ref is not None:
			self.image_ref = self._prepare_flux(self.image_ref)

		if self.wcs_ref is not None and not isinstance(self.wcs_ref, WCS):
			if not isinstance(self.wcs_ref, str):
				self.wcs_ref = self.wcs_ref.decode("utf-8") # For Python 3
			self.wcs_ref = WCS(header=fits.Header().fromstring(self.wcs_ref))

		self._interpolator = None

	#----------------------------------------------------------------------------------------------
	def __call__(self, *args, **kwargs):
		return self.apply_kernel(*args, **kwargs)

	#----------------------------------------------------------------------------------------------
	def _prepare_flux(self, flux):
		"""
		Preparation of images for Enhanced Correlation Coefficient (ECC) Maximization
		estimation of movement - used for estimation of jitter.

		Parameters:
			flux (array): flux pixel image

		Returns:
			array: Gradient (using Scharr method) of image in logarithmic units.

		.. codeauthor:: Mikkel N. Lund
		.. codeauthor:: Rasmus Handberg <rasmush@phys.au.dk>
		"""

		# Convert to logarithmic units, avoiding taking log if zero:
		flux = np.asarray(flux)
		flux = np.log10(flux - np.nanmin(flux) + 1.0)

		# Convert image to flux in range -1 to 1 (for gradient determination)
		fmax = np.nanmax(flux)
		fmin = np.nanmin(flux)
		ran = np.abs(fmax - fmin)
		flux1 = -1 + 2*((flux - fmin)/ran)

		# Calculate Scharr gradient
		flux1 = scharr(flux1)

		# Remove potential NaNs in gradient image
		flux1[np.isnan(flux1)] = 0

		# Make sure image is in proper units for ECC routine
		return np.asarray(flux1, dtype='float32')

	#----------------------------------------------------------------------------------------------
	def apply_kernel(self, xy, kernel):
		"""
		Application of warp matrix to pixel coordinates

		Parameters:
			xy (2D ndarray): 2D array of image positions to be transformed.
			kernel (1D ndarray): The kernel to transform against.

		Returns:
			ndarray: Change in positions compared to reference.

		.. codeauthor:: Mikkel N. Lund
		.. codeauthor:: Rasmus Handberg <rasmush@phys.au.dk>
		"""

		xy = np.atleast_2d(xy)
		delta_pos = np.empty_like(xy)

		if self.warpmode == 'wcs':
			# Calculate RA and DEC of target in the reference image:
			radec = self.wcs_ref.all_pix2world(xy, 0, ra_dec_order=True)
			# Use RA and DEC to find the position in the kernel image:
			# TODO: Better handling of NoConvergence exception, which is currently silenced
			delta_pos = kernel.all_world2pix(radec, 0, ra_dec_order=True, maxiter=50, quiet=True)
			# Calculate the difference in pixel-position:
			delta_pos -= xy

		elif self.warpmode == 'unchanged':
			delta_pos.fill(0)

		elif self.warpmode == 'translation':
			delta_pos[:, 0] = kernel[0]
			delta_pos[:, 1] = kernel[1]

		elif self.warpmode == 'euclidian':
			dx = kernel[0]
			dy = kernel[1]
			theta = kernel[2]

			# Set up warp matrix:
			c = np.cos(theta)
			s = np.sin(theta)
			R = np.array([[c, -s, dx], [s, c, dy]])

			# Apply warp to all positions:
			for i in range(xy.shape[0]):
				x = xy[i, 0]
				y = xy[i, 1]
				delta_pos[i, :] = np.dot(R, [x, y, 1])

			# Subtract the reference positions to return the change in positions:
			delta_pos -= xy

		elif self.warpmode == 'affine':
			# In the affine case, the kernel is simply the full matrix:
			warp_matrix = kernel.reshape(2, 3)

			# Apply warp to all positions:
			for i in range(xy.shape[0]):
				x = xy[i, 0]
				y = xy[i, 1]
				delta_pos[i, :] = np.dot(warp_matrix, [x, y, 1])

			# Subtract the reference positions to return the change in positions:
			delta_pos -= xy

		return delta_pos

	#----------------------------------------------------------------------------------------------
	def calc_kernel(self, image, number_of_iterations=10000, termination_eps=1e-6):
		"""
		Calculates the position movement kernel for a given image. This kernel is
		a set of numbers that can be passed to `apply_kernel` to calculate the movement
		of a star at specific coordinates.

		Calculation of Enhanced Correlation Coefficient (ECC) Maximization using OpenCV.

		Parameters:
			image (ndarray): Image to calculate kernel for.
			number_of_iterations (integer, optional): Specify the number of iterations.
			termination_eps (float, optional): Specify the threshold of the increment in the
				correlation coefficient between two iterations.

		Returns:
			ndarray:

		.. codeauthor:: Mikkel N. Lund
		.. codeauthor:: Rasmus Handberg <rasmush@phys.au.dk>
		"""

		logger = logging.getLogger(__name__)

		if self.warpmode == 'unchanged':
			return []

		# Check that reference image was actually given:
		if self.image_ref is None:
			raise Exception("Reference image not defined")

		# Define the motion model
		if self.warpmode == 'euclidian':
			warp_mode = cv2.MOTION_EUCLIDEAN
		elif self.warpmode == 'translation':
			warp_mode = cv2.MOTION_TRANSLATION
		elif self.warpmode == 'affine':
			warp_mode = cv2.MOTION_AFFINE

		# Prepare comparison image for estimation of motion
		image = self._prepare_flux(image)

		# Mask of good pixels in input image:
		mask = np.isfinite(image)
		mask = np.asarray(mask, dtype='uint8') # For some reason OpenCV doesn't like boolean arrays

		# Define 2x3 warp matrix and initialize the matrix to identity
		warp_matrix = np.eye(2, 3, dtype='float32')

		# Define termination criteria
		criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, number_of_iterations, termination_eps)

		# Run the ECC algorithm. The results are stored in warp_matrix.
		try:
			inputMask = cv2.UMat(mask)
			cc, warp_matrix = cv2.findTransformECC(self.image_ref, image, warp_matrix, warp_mode, criteria, inputMask, 5)
			warp_matrix = np.asarray(warp_matrix.get(), dtype='float64')
		except: # noqa: E722
			logger.exception("Could not find transform")
			return np.full(self.n_params, np.NaN)

		# The affine transformation uses six parameters,
		# so let's just store a flattened version of the matrix:
		if self.warpmode == 'affine':
			return warp_matrix.flatten()

		# Extract movement in pixel units in x- and y direction
		dx = warp_matrix[0,2]
		dy = warp_matrix[1,2]

		if self.warpmode == 'euclidian':
			# Estimate rotation angle in radians
			theta = math.atan2(warp_matrix[1,0], warp_matrix[0,0])
			return [dx, dy, theta]
		else:
			# Translation only:
			return [dx, dy]

	#----------------------------------------------------------------------------------------------
	def load_series(self, times, kernels):
		"""
		Load time-series of kernels and create interpolator.

		The interpolator (:py:func:`interpolator`) can be used to obtain movements
		at a arbitrery timestamp within the timestamps provided in ``times``.

		Parameters:
			time (1D array): Timestamps to be interpolated against. Timestamps must be sorted.
			kernels (2D array): List of kernels.

		Raises:
			ValueError: If kernels have the wrong shape.
		"""

		self.series_times = np.asarray(times)
		self.series_kernels = kernels

		if self.warpmode == 'wcs':
			# Check the lenghts of the provided vectors:
			if len(kernels) != len(times):
				raise ValueError("Wrong shape of kernels.")

			good_series = np.ones_like(self.series_times, dtype='bool')
			for k in range(len(kernels)):
				if not isinstance(self.series_kernels[k], WCS):
					# Assuming that is is a string then:
					hdr_string = self.series_kernels[k]
					if not isinstance(hdr_string, str):
						hdr_string = hdr_string.decode("utf-8") # For Python 3

					# If the string is empty, remove the point from the series:
					if hdr_string.strip() == '':
						good_series[k] = False
						continue

					# Create a WCS object from the header string:
					self.series_kernels[k] = WCS(header=fits.Header().fromstring(hdr_string), relax=True)

				# Try if the WCS can return pixel coordinates for the test-coordinates:
				# If it can't we will remove that timestamp from the series, and the
				# pixel coordinates will therefore be interpolated when calculation jitter.
				# Using a (ra, dec) coordinates from the actual footprint of the WCS, and
				# using axes=(2,2), since we here don't nessacerily know the size of the image,
				# and we are only using the first corner anyway.
				fp = self.series_kernels[k].calc_footprint(axes=(2, 2))
				test_coords = np.atleast_2d(fp[0, :])
				try:
					self.series_kernels[k].all_world2pix(test_coords, 0, ra_dec_order=True, maxiter=50)
				except (NoConvergence, ValueError):
					good_series[k] = False

			# Remove any bad series points from the lists:
			self.series_kernels = np.asarray(self.series_kernels)
			self.series_times = self.series_times[good_series]
			self.series_kernels = self.series_kernels[good_series]

		else:
			# For these warpmodes, the kernels should be 2D arrays:
			self.series_kernels = np.atleast_2d(self.series_kernels)

			# Check shape of the input:
			if self.series_kernels.shape != (len(self.series_times), self.n_params):
				raise ValueError("Wrong shape of kernels. Anticipated ({0},{1}), but got {2}".format(
					len(self.series_times),
					self.n_params,
					self.series_kernels.shape
				))

			# Only take the kernels that are well-defined:
			# TODO: Should we raise a warning if there are many undefined?
			indx = np.isfinite(times) & np.all(np.isfinite(kernels), axis=1)

			# Create interpolator:
			self._interpolator = interp1d(times[indx], kernels[indx, :],
				axis=0,
				assume_sorted=True,
				bounds_error=False,
				fill_value=(kernels[0, :], kernels[-1, :]))

	#----------------------------------------------------------------------------------------------
	def interpolate(self, time, xy):
		"""
		Interpolate in the kernel time-series provided in :py:func:`load_series`
		to obtain movment a arbitrery time.

		Parameters:
			time (float): Timestamp to return movement for.
			xy (2D array): row and column positions to be modified.

		Returns:
			``numpy.ndarray``: Array with the same size as `xy` containing the
				changes to rows and columns. These can be added
				to `xy` to yield the new positions.

		Raises:
			ValueError: If timeseries has not been provided.
		"""

		if self.warpmode == 'wcs':
			# Methods where the kernel is complex (non-numeric)
			# Handle the case where we are requesting a timestamp outside the
			# range of the loaded kernel timeseries:
			if time < self.series_times[0] or time > self.series_times[-1]:
				# Allow for a bit of a margin before and after the ends of the
				# timeseries, to account for e.g. round-off errors in the timestamps:
				dt = np.median(np.diff(self.series_times))
				if np.abs(time - self.series_times[0]) < dt:
					return self.apply_kernel(xy, self.series_kernels[0])
				elif np.abs(time - self.series_times[-1]) < dt:
					return self.apply_kernel(xy, self.series_kernels[-1])
				else:
					raise ValueError("Timestamp outside timeseries interval")

			# Find the point in the series where the timestamp falls:
			k = np.searchsorted(self.series_times, time, side='right')
			t1 = self.series_times[k-1]
			# Find the jitter in that kernel:
			jitter_1 = self.apply_kernel(xy, self.series_kernels[k-1])
			if t1 == time:
				# We actually hit spot on, so let's just return the jitter:
				return jitter_1
			else:
				#
				t2 = self.series_times[k]
				jitter_2 = self.apply_kernel(xy, self.series_kernels[k])

				int_time = [t1, t2]
				jitter_row = interp1d(int_time, np.column_stack((jitter_1[:,0], jitter_2[:,0])), axis=1, assume_sorted=True)
				jitter_col = interp1d(int_time, np.column_stack((jitter_1[:,1], jitter_2[:,1])), axis=1, assume_sorted=True)

				return np.column_stack((jitter_row(time), jitter_col(time)))
		else:
			#
			if self._interpolator is None:
				raise ValueError("Interpolator is not defined. ")

			# Get the kernel parameters for the timestamp:
			with warnings.catch_warnings():
				warnings.filterwarnings('ignore', category=RuntimeWarning, module='scipy')
				kernel = self._interpolator(time)

			return self.apply_kernel(xy, kernel)

	#----------------------------------------------------------------------------------------------
	def jitter(self, time, column, row):
		"""
		Calculate the change to a given position as a function of time.

		Parameters:
			time (ndarray): Array of timestamps to calculate position changes for.
			column (float): Column position at reference time.
			row (float): Row position at reference time.

		Returns:
			ndarray: 2D array with changes in column and row for each timestamp.
		"""

		xy = np.array([column, row])
		jtr = np.empty((len(time), 2), dtype='float64')
		for k in range(len(time)):
			jtr[k, :] = self.interpolate(time[k], xy)

		return jtr

max_stars_repo_name: aditya-sengupta/tesscomp-prototyping

max_issues_repo_name: aditya-sengupta/tesscomp-prototyping

max_forks_repo_name: aditya-sengupta/tesscomp-prototyping

reason: import numpy,from scipy,from astropy


text:
import numpy as np
from numpy import testing
import pytest

from .basis import Basis
from .grid import Grid
from .cell import Cell

def _parabolic(n, d2=False):
    grid = [np.linspace(-.5, .5, n, endpoint=False)] * 2
    if not d2:
        grid.append(np.array((0,)))
    mg = np.meshgrid(*grid, indexing='ij')
    data = sum(i ** 2 for i in mg) ** .5

    return Grid(Basis.orthorhombic((1, 1) if d2 else (1, 1, 1)), grid, data)

grid1 = _parabolic(1)
cell1 = grid1.as_cell()
grid50 = _parabolic(50)
cell50 = grid50.as_cell()
grid50_2d = _parabolic(50, True)
cell50_2d = grid50_2d.as_cell()
cell = Cell(Basis.orthorhombic((1, 1)), [[.2, .3], [.6, .1], [0, .5]], [1, 2, 3])

@pytest.mark.parametrize("grid", (grid1, cell1, grid50, cell50, grid50_2d, cell50_2d, cell))
def test_triangulation(grid):
    tri = grid.compute_triangulation()
    testing.assert_allclose(tri.weights.sum(), 1)

@pytest.mark.parametrize("grid", (grid50, cell50))
def test_default(grid):
    d = grid.tetrahedron_density((-.1, 0, .1, .2))
    testing.assert_allclose(d, (0, 0, 2 * np.pi * 0.1, 2 * np.pi * 0.2), rtol=1e-2)

def test_resolved():
    d = grid50.tetrahedron_density((-.1, 0, .1, .2), resolved="grid")
    testing.assert_equal(d.values.shape, (*grid50.grid_shape, 4))
    testing.assert_allclose(d.values.sum(axis=(0, 1, 2)), (0, 0, 2 * np.pi * 0.1, 2 * np.pi * 0.2), rtol=1e-2)

@pytest.mark.parametrize("grid", (grid50, cell50))
def test_weighted(grid):
    d = grid.tetrahedron_density((-.1, 0, .1, .2), weights=np.ones_like(grid.values) * .5)
    testing.assert_allclose(d, (0, 0, np.pi * 0.1, np.pi * 0.2), rtol=1e-2)

def test_fail():
    g = grid50.copy(
        vectors=grid50.vectors[:2, :2],
        coordinates=grid50.coordinates[:2],
        values=grid50.values[:, :, 0, ...],
    )
    with pytest.raises(AssertionError):
        g.tetrahedron_density((-.1, 0, .1, .2))

max_stars_repo_name: pulkin/pycoordinates

max_issues_repo_name: pulkin/pycoordinates

max_forks_repo_name: pulkin/pycoordinates

reason: import numpy,from numpy


text:
import networkx as nx
import numpy as np
import math
import matplotlib.pyplot as plt

class HashFamily():
    """Universal hash family
    Args:
        bins (int): Number of bins to hash to. Better if a prime number.
        moduler (int,optional): Temporary hashing. Has to be a prime number.
    """
    def __init__(self, bins, moduler=None):
        if moduler and moduler <= bins:
            raise ValueError("p (moduler) should be >> m (buckets)")

        self.bins = bins
        self.moduler = moduler if moduler else self._next_prime(np.random.randint(self.bins + 1, 2**32))

        # do not allow same a and b, as it could mean shifted hashes
        self.sampled_a = set()
        self.sampled_b = set()

    def _is_prime(self, x):
        """Naive is prime test."""
        for i in range(2, int(np.sqrt(x))):
            if x % i == 0:
                return False
        return True

    def _next_prime(self, n):
        """Naively gets the next prime larger than n."""
        while not self._is_prime(n):
            n += 1

        return n

    def draw_hash(self, a=None, b=None):
        """Draws a single hash function from the family."""
        if a is None:
            while a is None or a in self.sampled_a:
                a = np.random.randint(1, self.moduler - 1)
                assert len(self.sampled_a) < self.moduler - 2, "please give a bigger moduler"

            self.sampled_a.add(a)
        if b is None:
            while b is None or b in self.sampled_b:
                b = np.random.randint(0, self.moduler - 1)
                assert len(self.sampled_b) < self.moduler - 1, "please give a bigger moduler"

            self.sampled_b.add(b)

        return lambda x: ((a * x + b) % self.moduler) % self.bins

    def draw_hashes(self, n, **kwargs):
        """Draws n hash function from the family."""
        return [self.draw_hash() for i in range(n)]

def power2(color):
    return 2**(color)

def dp_helper(dp_mat, cur_color, G, colors, k, khop_set, start_node, subgraph='path'):
    new_khop_set = set()
    if len(khop_set)==0:
        return []
    for col in range(dp_mat.shape[1]):
#         find the correct columns representing current number of color
        if bin(col).count("1")!=cur_color-1:
            continue
        for v in khop_set:
            if dp_mat[v][col]>0:
                for neighbor in G[v]:
                    aggregate_color = (power2(colors[neighbor])|col)
                    if bin(aggregate_color).count("1")==cur_color:
                        new_khop_set.add(neighbor)
                        dp_mat[neighbor][aggregate_color] = 1
    if cur_color==k:
        if subgraph=='path':
            return np.where(dp_mat[:,-1]==1)[0]
        elif subgraph=='cycle':
            return set(np.where(dp_mat[:,-1]==1)[0]) & set(G[start_node])
    return dp_helper(dp_mat, cur_color+1, G, colors, k, new_khop_set, start_node, subgraph)

def backtrack(dp_mat, G, colors, target, s, prev_path, prev_colors, verbose):
    if target==s:
        if verbose:
            print(prev_path)
        global_path.add(','.join(prev_path))
        prev_path.pop()
        return 1
    count = 0
    for neigh in G[target]:
        not_visited = dp_mat.shape[1]-1^prev_colors
        if is_set(prev_colors, colors[neigh]):
            continue
        if dp_mat[neigh][not_visited]>0:
            count += backtrack(dp_mat, G, colors, neigh, s, prev_path+[str(neigh)], set_bit(prev_colors, colors[neigh]), verbose)
    prev_path.pop()
    return count

def colorful_subgraph(G, colors, k, s, subgraph_type, get_detail=True, verbose=False):
    """Detect if colorful path exists fom s to any node by dynamic programming.
    Args:
        G (nx.Graph): with n nodes and m edges
        colors (list): list of integers represents node colors
        k (int): number of colors
        s (int): starting node
        subgraph (str): 'path' or 'cycle'
    Return:
        List: nodes connect to s with at least one colorful path
    """
    n = G.number_of_nodes()
    col = 2**k
    dp_mat = np.zeros((n, col))
    dp_mat[s][power2(colors[s])] = 1
    targets = dp_helper(dp_mat, 2, G, colors, k, set([s]), s, subgraph_type)
    if not get_detail:
        return targets
    else:
        empty_color = 0
        total_count = 0
        for target in targets:
            total_count += backtrack(dp_mat, G, colors, target, s, [str(target)], set_bit(empty_color, colors[target]), verbose)
        if verbose:
            print('from node', s, 'find in total', total_count, 'colorful paths of length', k)
        return total_count

def color_by_hashes(G, k, repeat, modulor=None):
    hashes = HashFamily(k, modulor)
    hash_list = hashes.draw_hashes(repeat)
    color_list = []
    for hashF in hash_list:
        color_list.append([hashF(i) for i in G.nodes])
    return color_list

def colorful_cycle(G, colors, k, s):
    n = G.number_of_nodes()
    col = 2**k
    dp_mat = np.zeros((n, col))
    dp_mat[s][power2(colors[s])] = 1
    return dp_helper(dp_mat, 2, G, colors, k, set([s]), s, 'cycle')

def is_set(x, n):
    '''
    Check if bit n of integer x is set.
    '''
    return x & 1 << n != 0

def set_bit(x, n):
    return x | (1<<n)

def backtrack(dp_mat, G, colors, target, s, prev_path, prev_colors, verbose):
    if target==s:
        if verbose:
            print(prev_path)
        global_path.add(','.join(prev_path))
        prev_path.pop()
        return 1
    count = 0
    for neigh in G[target]:
        not_visited = dp_mat.shape[1]-1^prev_colors
        if is_set(prev_colors, colors[neigh]):
            continue
        if dp_mat[neigh][not_visited]>0:
            count += backtrack(dp_mat, G, colors, neigh, s, prev_path+[str(neigh)], set_bit(prev_colors, colors[neigh]), verbose)
    prev_path.pop()
    return count

def find_subgraph(G, k, subgraph, color_repeat, get_detail=True, verbose=False, plot=True):
    count = 0
    global global_path
    global_path = set()

    repeat = math.ceil(math.exp(k))
    colors = color_by_hashes(G, k, repeat, modulor=97)
    for color in colors:
        if plot:
            nx.draw(G, node_color=color, with_labels=True)
            plt.show()
        for i in G.nodes():
            tmp_res = colorful_subgraph(G, color, k, i, subgraph, get_detail, verbose)
    #         print(i, tmp_res)
        count += 1
        if count >= color_repeat:
            break
    print('set of', subgraph, 'found:', global_path)

reason: import numpy,import networkx


text:
import numpy as np
import pandas as pd
import argparse
import sys
import os

sys.path.append(os.path.dirname(os.path.abspath(__file__)) + '/../')
from db import engine, Session
from common.time_utils import get_local_datetime
from config.settings import MIN_LAT, MIN_LON, DELTA_LAT, DELTA_LON, MAP_WIDTH, MAP_HEIGHT,\
    GLOBAL_STATE_UPDATE_CYCLE, DESTINATION_PROFILE_TEMPORAL_AGGREGATION, DESTINATION_PROFILE_SPATIAL_AGGREGATION
# from dqn.settings import FLAGS

# def create_predicted_demand(source_table, profile_table, variance_rate=None):
#     query = "SELECT * FROM {}".format(source_table)
#     df = pd.read_sql(query, engine, index_col="id")
#     print("# of rows {}".format(len(df)))
#
#     df["datetime_obj"] = df.request_datetime.apply(lambda x: get_local_datetime(x))
#     df["dayofweek"] = df.datetime_obj.apply(lambda x: x.weekday())
#     df["hour"] = df.datetime_obj.apply(lambda x: x.hour)
#     df['day'] = df.datetime_obj.apply(lambda x: x.day)
#     df['month'] = df.datetime_obj.apply(lambda x: x.month)
#     df['x'] = ((df['origin_lon'] - MIN_LON) / DELTA_LON).astype(int)
#     df['y'] = ((df['origin_lat'] - MIN_LAT) / DELTA_LAT).astype(int)
#
#     df = df.groupby(['month', 'day', 'hour', 'x', 'y']).size()
#     if variance_rate:
#         df += np.maximum(-df.values, np.minimum(df.values, np.random.normal(0, np.sqrt(variance_rate * df.values))))
#     df.name = 'demand'
#     df = df.reset_index()
#
#     drop_table = """
#     DROP TABLE IF EXISTS {};
#     """.format(profile_table)
#     Session.execute(drop_table)
#     Session.commit()
#     df.to_sql(profile_table, engine, flavor=None, schema=None, if_exists='fail',
#                index=True, index_label=None, chunksize=None, dtype=None)
#
#     create_index = """
#     CREATE INDEX index_{table} ON {table} (month, day, hour);
#     """.format(table=profile_table)
#     Session.execute(create_index)
#     Session.commit()

def create_demand_profile(df, profile_table, n_weeks):
    df["dayofweek"] = df.datetime_obj.apply(lambda x: x.weekday())
    df["hour"] = df.datetime_obj.apply(lambda x: x.hour)
    df['x'] = ((df['origin_lon'] - MIN_LON) / DELTA_LON).astype(int)
    df['y'] = ((df['origin_lat'] - MIN_LAT) / DELTA_LAT).astype(int)

    profile_df = df.groupby(['dayofweek', 'hour', 'x', 'y']).size() / float(n_weeks)
    profile_df.name = 'demand'
    drop_table = """
    DROP TABLE IF EXISTS {};
    """.format(profile_table)
    Session.execute(drop_table)
    Session.commit()
    profile_df.reset_index().to_sql(profile_table, engine, flavor=None, schema=None, if_exists='fail',
               index=True, index_label=None, chunksize=None, dtype=None)
    create_index = """
    CREATE INDEX index_{table} ON {table} (dayofweek, hour);
    """.format(table=profile_table)
    Session.execute(create_index)
    Session.commit()

def create_od_profile(df, profile_table, n_weeks):
    hours_bin = DESTINATION_PROFILE_TEMPORAL_AGGREGATION
    n_agg = DESTINATION_PROFILE_SPATIAL_AGGREGATION
    df["dayofweek"] = df.datetime_obj.apply(lambda x: x.weekday())
    df["hours_bin"] = (df.datetime_obj.apply(lambda x: x.hour) / hours_bin).astype(int)
    df['origin_x'] = ((df['origin_lon'] - MIN_LON) / (DELTA_LON * n_agg)).astype(int)
    df['origin_y'] = ((df['origin_lat'] - MIN_LAT) / (DELTA_LAT * n_agg)).astype(int)
    df['destination_x'] = ((df['destination_lon'] - MIN_LON) / (DELTA_LON * n_agg)).astype(int)
    df['destination_y'] = ((df['destination_lat'] - MIN_LAT) / (DELTA_LAT * n_agg)).astype(int)

    od_df = df.groupby(['dayofweek', 'hours_bin', 'origin_x', 'origin_y', 'destination_x', 'destination_y']
                        ).trip_time.agg(['count', 'mean'])
    od_df = od_df.rename(columns={'count' : 'demand', 'mean' : 'trip_time'}).reset_index()
    drop_table = """
    DROP TABLE IF EXISTS {};
    """.format(profile_table)
    Session.execute(drop_table)
    Session.commit()
    od_df.to_sql(profile_table, engine, flavor=None, schema=None, if_exists='fail',
               index=True, index_label=None, chunksize=None, dtype=None)

    create_index = """
    CREATE INDEX index_{table} ON {table} (dayofweek, hours_bin);
    """.format(table=profile_table)
    Session.execute(create_index)
    Session.commit()

def create_latest_demand(source_table, latest_table):
    query = "SELECT * FROM {}".format(source_table)
    df = pd.read_sql(query, engine, index_col="id")
    print("# of rows {}".format(len(df)))

    unit_min = GLOBAL_STATE_UPDATE_CYCLE / 60
    df["t"] = (df.request_datetime / unit_min).astype(int) * unit_min
    df['x'] = ((df['origin_lon'] - MIN_LON) / DELTA_LON).astype(int)
    df['y'] = ((df['origin_lat'] - MIN_LAT) / DELTA_LAT).astype(int)

    latest_df = df.groupby(['t', 'x', 'y']).size()
    latest_df.name = 'demand'
    drop_table = """
    DROP TABLE IF EXISTS {};
    """.format(latest_table)
    Session.execute(drop_table)
    Session.commit()
    latest_df.reset_index().to_sql(latest_table, engine, flavor=None, schema=None, if_exists='fail',
               index=True, index_label=None, chunksize=None, dtype=None)
    create_index = """
    CREATE INDEX index_{table} ON {table} (t);
    """.format(table=latest_table)
    Session.execute(create_index)
    Session.commit()

def create_training_dataset(df, n_weeks):
    t_start = df.request_datetime.min()
    t_end = t_start + 3600 * 24 * 7 * n_weeks
    df = df[(df.request_datetime >= t_start) & (df.request_datetime < t_end)]
    df["datetime_obj"] = df.request_datetime.apply(lambda x: get_local_datetime(x))
    return df

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument("input_file", help = "path to the training data file")
    args = parser.parse_args()

    df = pd.read_csv(args.input_file, index_col='id')
    print("load {} rows".format(len(df)))

    n_weeks = 4
    df = create_training_dataset(df, n_weeks)
    print("created training dataset")

    create_demand_profile(df, "demand_profile", n_weeks)
    print("created demand_profile table")

    create_od_profile(df, "od_profile", n_weeks)
    print("created od_profile table")

    create_latest_demand("request_backlog", "demand_latest")
    print("created demand_latest table")
    # t_start = 1462075200 # 160501
    # sim_datetime = 1464753600 # 160601


text:
import abc
import numpy as np

from latent_with_splitseqs.memory.replay_buffer_for_latent import LatentReplayBuffer

import self_supervised.utils.typed_dicts as td

from my_utils.np_utils.take_per_row import take_per_row
from my_utils.np_utils.np_array_equality import np_array_equality

class LatentReplayBufferSplitSeqSamplingBase(LatentReplayBuffer,
                                             metaclass=abc.ABCMeta):
    def __init__(
            self,
            *args,
            min_sample_seqlen: int = 2,
            **kwargs
    ):
        super().__init__(*args, **kwargs)
        self._min_sample_seqlen = min_sample_seqlen

    @abc.abstractmethod
    def _get_sample_seqlen(self) -> int:
        raise NotImplementedError

    @property
    def horizon_len(self):
        return self._seq_len

    @staticmethod
    def _add_left_zero_padding(
            arr: np.ndarray,
            padding_len: int,
            padding_dim: int
    ) -> np.ndarray:
        size_padding = list(arr.shape)
        size_padding[padding_dim] = padding_len
        padding_arr = np.zeros(size_padding)
        padded = np.concatenate([padding_arr, arr], axis=padding_dim)
        return padded

    def _take_elements_from_batch(
            self,
            batch: dict,
            batch_size: int,
            sample_seq_len: int,
    ) -> dict:
        """
        Prepends zero padding and samples sequences of seq_len from batch
        Args:
            batch               : dict with (N, data_dim, S) numpy arrays
            batch_size          : N
            sample_seq_len      : length of sampling that are taken along dimension S
        Return:
            sampled_batch       : dict with (N, data_dim, sample_seq_len) numpy arrays
        """
        batch_dim = 0
        data_dim = 1
        seq_dim = 2
        skill_key = 'mode'
        next_obs_key = 'next_obs'
        horizon_len = batch[next_obs_key].shape[seq_dim]
        transition_mode_mapping_kwargs = {}
        start_idx = np.random.randint(
            low=0,
            high=self.horizon_len, # padding is needed
            size=(batch_size,)
        )

        if self._padding:
            for key, el in batch.items():
                if isinstance(el, np.ndarray) and key != skill_key:
                    el_padded = self._add_left_zero_padding(
                        arr=el,
                        padding_len=sample_seq_len,
                        padding_dim=seq_dim,
                    )

                elif isinstance(el, np.ndarray) and key == skill_key:
                    assert np_array_equality(
                        el,
                        np.stack([el[..., 0]] * horizon_len, axis=seq_dim)
                    )
                    padding_arr = np.stack([el[..., 0]] * sample_seq_len, axis=seq_dim)
                    el_padded = np.concatenate([padding_arr, el], axis=seq_dim)

                else:
                    el_padded = None

                # Sample from seqs
                if isinstance(el_padded, np.ndarray):
                    el_bsd = np.swapaxes(el_padded, axis1=data_dim, axis2=seq_dim)
                    el_slices_bsd = take_per_row(el_bsd, start_idx, num_elem=sample_seq_len)
                    transition_mode_mapping_kwargs[key] = \
                        np.swapaxes(el_slices_bsd, axis1=1, axis2=2)

        else:
            end_idx_bool = start_idx > self.horizon_len - sample_seq_len
            start_idx[end_idx_bool] = np.random.randint(
                self.horizon_len - sample_seq_len,
                size=(batch_size,)
            )

            for key, el in batch.items():
                if isinstance(el, np.ndarray):
                    el_bsd = np.swapaxes(el, axis1=data_dim, axis2=seq_dim)
                    el_slices_bsd = take_per_row(el_bsd, start_idx, num_elem=sample_seq_len)
                    transition_mode_mapping_kwargs[key] \
                        = np.swapaxes(el_slices_bsd, axis1=1, axis2=2)

        return transition_mode_mapping_kwargs

    def random_batch(self,
                     batch_size: int) -> td.TransitionModeMapping:
        """
        Sample batches with random sequence length
        Args:
            batch_size                 : N
        Return:
            TransitionModeMapping      : consisting of (N, data_dim, S) tensors
        """
        sample_seq_len = self._get_sample_seqlen()
        batch_horizon = super().random_batch(
            batch_size=batch_size,
        )
        transition_mode_mapping_kwargs = self._take_elements_from_batch(
            batch=batch_horizon,
            batch_size=batch_size,
            sample_seq_len=sample_seq_len,
        )

        return td.TransitionModeMapping(
            **transition_mode_mapping_kwargs
        )

    def random_batch_latent_training(self,
                                     batch_size: int) -> dict:
        """
        Sample batches with random sequence length
        Args:
            batch_size              : N
        """
        sample_seq_len = self._get_sample_seqlen()
        batch_horizon = super().random_batch_latent_training(
            batch_size=batch_size,
        )
        return_dict = self._take_elements_from_batch(
            batch=batch_horizon,
            sample_seq_len=sample_seq_len,
            batch_size=batch_size,
        )

        return return_dict

    def _add_sample_if(self, path_len) -> bool:
        """
        No terminal handling in this class
        """
        return True

max_stars_repo_name: fgitmichael/SelfSupevisedSkillDiscovery

max_issues_repo_name: fgitmichael/SelfSupevisedSkillDiscovery

max_forks_repo_name: fgitmichael/SelfSupevisedSkillDiscovery


text:
import numpy
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import math

SIZE_X = 128

u = numpy.zeros(SIZE_X)
b = numpy.zeros(SIZE_X)

for x in range(40,50):
    b[x] = 1.0

# b[8,8] = 1

# def gauss_seidel(U,b,theta):
#     rows = U.shape[0]
#     T = numpy.zeros_like(U)
#     for x in range(theta, (rows - theta),theta):
#         T[x] = (U[x+theta] +
#                     U[x-theta] + b[x])/2.0
#     numpy.copyto(U,T)
#     return T

def residual(u, f):
    """
    f - A u
    """
    n = len(f)
    r = numpy.zeros(len(u))
    r[1:-1] = f[1:-1] - ((n+1)**2) * (2 * u[1:-1] - u[2:] - u[:-2])
    r[0]    = f[0]    - ((n+1)**2) * (2 * u[0]    - u[1])
    r[-1]   = f[-1]   - ((n+1)**2) * (2 * u[-1]   - u[-2])
    return r

# def residual(phi,f):
#     r = numpy.zeros_like(phi)
#     rows = phi.shape[0]
#
#     for x in range(1, (rows - 1),1): #(((rows**2.0) / 1.0)*
#         r[x] = (phi[x+1] + phi[x-1] - 2.0*phi[x])/2.0
#     return r

# def construct_T(U):
#     rows = U.shape[0]
#     T = numpy.zeros((rows,rows))
#
#     for x in range(1, (rows - 1)):
#         T[x,x-1] = 1
#         T[x,x] = -2
#         T[x,x+1] = 1
#
#
#
#     return T

#
def restriction(X):
    rows = X.shape[0]

    O = numpy.zeros((int(rows/2)))

    for x in range(0,int(rows/2)-1):
        O[x] = X[x*2]

    return O

def prolongate(X):
    rows = X.shape[0]

    O = numpy.zeros((rows*2))
    for x in range(1,rows-1):
        O[2*x] = X[x]
        O[2*x+1] = (X[x]+X[x+1])/2.0

    return O

# it's not the root residual you idiot
# it's each residual recursively

def V_Cycle(u,f,h):
    phi = u.copy()

    gauss_seidel(phi,f,1)

    r = residual(phi,f)

    rhs = restriction(r)
    #
    eps = numpy.zeros_like(rhs)
    #
    if(h == 32):
        for i in range(0,10):
            gauss_seidel(eps,rhs,1) #this is the error or correction equation.
    else:
        eps = V_Cycle(eps,rhs,2*h)

    T = prolongate(eps)

    if(h==1):
        print("Residual: {} Step: {}".format(numpy.linalg.norm(r)/(len(phi)**2),h))

        plt.subplot(2, 3, 1)
        plt.gca().set_title('Potentials')
        plt.plot(phi)

        plt.subplot(2, 3, 2)
        plt.gca().set_title('Residuals')
        plt.plot(r)

        plt.subplot(2, 3, 3)
        plt.gca().set_title('Restricted')
        plt.plot(rhs)

        plt.subplot(2, 3, 4)
        plt.gca().set_title('Correction')
        plt.plot(eps)

        plt.subplot(2, 3, 5)
        plt.gca().set_title('Prolongated')
        plt.plot(T)

        plt.draw()
        plt.pause(0.1)
        plt.clf()
        plt.cla()

    return phi + T

convergence = []
ims = []
t = 0
c1 = 1

v1 = u.copy()

r1 = u-v1

while True:
# for i in range(0,5):
    # u = restriction(b)
    # u = restriction(u)
    #
    # plt.subplot(2, 2, 1)
    # plt.gca().set_title('Potentials')
    # plt.imshow(b)
    #
    # plt.subplot(2, 2, 2)
    # plt.gca().set_title('Potentials')
    # plt.imshow(u)
    #
    # T = prolongate(u)
    #
    # plt.subplot(2, 2, 3)
    # plt.gca().set_title('Potentials')
    # plt.imshow(T)

    # u = V_Cycle(u,b,1)

    # gauss_seidel(u,b,1)
    # r = residual(u,b)
    # r = restriction(r)
    #
    # v = numpy.zeros_like(r)
    # gauss_seidel(v,r,1)
    # print(numpy.linalg.norm(r)/(len(u)**2))
    # print(numpy.linalg.norm(v))

    T = numpy.zeros_like(u)
    for x in range(1, (u.shape[0] - 2)):
        T[x] = (u[x+1] + u[x-1] + b[x])/2.0

    u = T.copy()
    plt.subplot(2, 3, 1)
    plt.gca().set_title('Potentials')
    plt.plot(u)

    # plt.subplot(2, 3, 2)
    # plt.gca().set_title('Correction')
    # plt.plot(r)
    #
    # plt.subplot(2, 3, 3)
    # plt.gca().set_title('Prolongated')
    # plt.plot(v)

    plt.draw()
    plt.pause(0.1)
    plt.clf()
    plt.cla()

    # u += prolongate(v)

    #
    # convergence.append(numpy.linalg.norm(r))
    #

    #
    # u = u + v
    #
    #
    # plt.subplot(2, 2, 2)
    # plt.gca().set_title('Residual')
    # plt.imshow(r)
    # plt.subplot(2, 2, 3)
    # plt.gca().set_title('Correction')
    # plt.imshow(v)
    # plt.subplot(2, 2, 4)
    # plt.yscale('log')
    # plt.gca().set_title('Convergence')
    # plt.plot(convergence)
    # # plt.savefig(str(t) + '.png')
    # t+=1
    # print("Residual: {} convergence factor: {} Step: {}".format(numpy.linalg.norm(r),numpy.linalg.norm(r)/c1,t))
    # c1 = numpy.linalg.norm(r)
    #

    # plt.close()