text:
import tensorflow as tf

import numpy as np
import json

import keras
from keras.callbacks import Callback
from keras.models import Model
from keras.layers import *
from keras import losses
import keras.backend as K

import sys
sys.path.append(sys.path[0]+"/LatentSpacePhysics/src/")

from LatentSpacePhysics.src.util.filesystem import make_dir

from itertools import chain
from collections import defaultdict
from ops import *

# --------------------------------------------------------------------------------------------------------------------------------------------------
# Network Helper Functions -------------------------------------------------------------------------------------------------------------------------
# --------------------------------------------------------------------------------------------------------------------------------------------------
# Freeze weights in all layers of a network
def make_layers_trainable(net, val, recurse=False, prefix=None):
    net.trainable = val
    for l in net.layers:
        if prefix is None:
            l.trainable = val
        if prefix and l.name.startswith(prefix):
            l.trainable = val
        if recurse and isinstance(l, Model):
            make_layers_trainable(l, val, recurse, prefix)

# --------------------------------------------------------------------------------------------------------------------------------------------------
# Copyright Ferry Boender, released under the MIT license.
# https://www.electricmonk.nl/log/2017/05/07/merging-two-python-dictionaries-by-deep-updating/
def deepupdate(target, src):
    for k, v in src.items():
        if type(v) == list:
            if not k in target:
                target[k] = copy.deepcopy(v)
            else:
                target[k].extend(v)
        elif type(v) == dict:
            if not k in target:
                target[k] = copy.deepcopy(v)
            else:
                deepupdate(target[k], v)
        elif type(v) == set:
            if not k in target:
                target[k] = v.copy()
            else:
                target[k].update(v.copy())
        else:
            target[k] = copy.copy(v)

# --------------------------------------------------------------------------------------------------------------------------------------------------
def merge_dicts(dict1, dict2):
    dict_res = defaultdict(list)
    for k, v in chain(dict1.items(), dict2.items()):
        dict_res[k].append(v)
    return  dict_res

# --------------------------------------------------------------------------------------------------------------------------------------------------
def merge_histories(hist_1, hist_2):
    hist_1.epoch += hist_2.epoch
    deepupdate(hist_1.history, hist_2.history)
    return hist_1

# --------------------------------------------------------------------------------------------------------------------------------------------------
# serialize model to JSON
def model_to_json(model, path):
    model_json = model.to_json()
    with open(path, "w") as json_file:
        json_file.write(model_json)

# --------------------------------------------------------------------------------------------------------------------------------------------------
# Callbacks ----------------------------------------------------------------------------------------------------------------------------------------
# --------------------------------------------------------------------------------------------------------------------------------------------------
class SaveCheckpoint(Callback):
    """See ModelCheckpoint documentation"""

    def __init__(self, filepath, network, monitor='val_loss', verbose=0,
                 save_best_only=False, mode='auto', period=1):
        super(SaveCheckpoint, self).__init__()
        self.network = network
        self.monitor = monitor
        self.verbose = verbose
        self.filepath = filepath
        make_dir(self.filepath)
        self.save_best_only = save_best_only
        self.period = period
        self.epochs_since_last_save = 0
        self.history = []

        if mode not in ['auto', 'min', 'max']:
            warnings.warn('SaveCheckpoint mode %s is unknown, '
                          'fallback to auto mode.' % (mode),
                          RuntimeWarning)
            mode = 'auto'

        if mode == 'min':
            self.monitor_op = np.less
            self.best = np.Inf
        elif mode == 'max':
            self.monitor_op = np.greater
            self.best = -np.Inf
        else:
            if 'acc' in self.monitor or self.monitor.startswith('fmeasure'):
                self.monitor_op = np.greater
                self.best = -np.Inf
            else:
                self.monitor_op = np.less
                self.best = np.Inf

    def on_epoch_end(self, epoch, logs=None):
        logs = logs or {}
        self.epochs_since_last_save += 1
        if self.epochs_since_last_save >= self.period:
            self.epochs_since_last_save = 0
            filepath = self.filepath #.format(epoch=epoch + 1, **logs)
            if self.save_best_only:
                current = logs.get(self.monitor)
                if current is None:
                    warnings.warn('Can save best model only with %s available, '
                                  'skipping.' % (self.monitor), RuntimeWarning)
                else:
                    if self.monitor_op(current, self.best):
                        if self.verbose > 0:
                            print('\nEpoch %05d: %s improved from %0.5f to %0.5f,'
                                  ' saving model to %s'
                                  % (epoch + 1, self.monitor, self.best,
                                     current, filepath))
                        self.best = current
                        self.network.save_model(filepath)
                        self.history.append("{}: {}".format(epoch, self.best))
                    else:
                        if self.verbose > 0:
                            print('\nEpoch %05d: %s did not improve from %0.5f' %
                                  (epoch + 1, self.monitor, self.best))
            else:
                if self.verbose > 0:
                    print('\nEpoch %05d: saving model to %s' % (epoch + 1, filepath))
                self.network.save_model(filepath)
                self.history.append("{}: {}".format(epoch, logs.get(self.monitor)))
        with open(self.filepath+"/history.json", "w") as outfile:
            json.dump(self.history, outfile)

# --------------------------------------------------------------------------------------------------------------------------------------------------
class PlotAEFields(Callback):
    def __init__(self, ae_func, x, func, path, batch_manager, name="AE_EncDec"):
        self._ae_func = ae_func
        self._x = x
        self._func = func
        self._counter = 0
        self._path = path
        self._batch_manager = batch_manager
        self.name = name

    def on_epoch_end(self, acc, loss):
        self._y = self._ae_func(self._x)[0]
        self._func(self._y, self._counter, self._path, self._batch_manager, self.name)
        self._counter += 1

# --------------------------------------------------------------------------------------------------------------------------------------------------
class StatefulResetCallback(Callback):
    def __init__(self, model):
        self.model = model

    def on_batch_end(self, batch, logs={}):
        self.model.reset_states()

# --------------------------------------------------------------------------------------------------------------------------------------------------
# Layer --------------------------------------------------------------------------------------------------------------------------------------------

# --------------------------------------------------------------------------------------------------------------------------------------------------
def upscale_layer(x, data_format='NHWC', is_3d=False):
    if is_3d:
        return Lambda(upscale3, arguments={'scale': 2})(x)
    else:
        return Lambda(upscale, arguments={'scale': 2, 'data_format': 'NHWC'})(x)

# --------------------------------------------------------------------------------------------------------------------------------------------------
def conv_layer( x,
                filters,
                kernel_size,
                stride=1,
                padding='valid',
                activation=None,
                kernel_initializer='glorot_uniform',
                kernel_regularizer=None,
                is_3d=False):
    if is_3d:
        return Conv3D(filters, kernel_size, strides=(stride,stride, stride), activation=activation, padding=padding, kernel_initializer=kernel_initializer, kernel_regularizer=kernel_regularizer)(x)
    else:
        return Conv2D(filters, kernel_size, strides=(stride,stride), activation=activation, padding=padding, kernel_initializer=kernel_initializer, kernel_regularizer=kernel_regularizer)(x)

# --------------------------------------------------------------------------------------------------------------------------------------------------
def jacobian_layer(x, data_format='NHWC', is_3d=False):
    if is_3d:
        dudx = x[:,:,:,1:,0] - x[:,:,:,:-1,0]
        dvdx = x[:,:,:,1:,1] - x[:,:,:,:-1,1]
        dwdx = x[:,:,:,1:,2] - x[:,:,:,:-1,2]
        dudy = x[:,:,:-1,:,0] - x[:,:,1:,:,0] # horizontally flipped
        dvdy = x[:,:,:-1,:,1] - x[:,:,1:,:,1] # horizontally flipped
        dwdy = x[:,:,:-1,:,2] - x[:,:,1:,:,2] # horizontally flipped
        dudz = x[:,1:,:,:,0] - x[:,:-1,:,:,0]
        dvdz = x[:,1:,:,:,1] - x[:,:-1,:,:,1]
        dwdz = x[:,1:,:,:,2] - x[:,:-1,:,:,2]

        dudx = K.concatenate([dudx, K.expand_dims(dudx[:,:,:,-1], axis=3)], axis=3)
        dvdx = K.concatenate([dvdx, K.expand_dims(dvdx[:,:,:,-1], axis=3)], axis=3)
        dwdx = K.concatenate([dwdx, K.expand_dims(dwdx[:,:,:,-1], axis=3)], axis=3)

        dudy = K.concatenate([K.expand_dims(dudy[:,:,0,:], axis=2), dudy], axis=2)
        dvdy = K.concatenate([K.expand_dims(dvdy[:,:,0,:], axis=2), dvdy], axis=2)
        dwdy = K.concatenate([K.expand_dims(dwdy[:,:,0,:], axis=2), dwdy], axis=2)

        dudz = K.concatenate([dudz, K.expand_dims(dudz[:,-1,:,:], axis=1)], axis=1)
        dvdz = K.concatenate([dvdz, K.expand_dims(dvdz[:,-1,:,:], axis=1)], axis=1)
        dwdz = K.concatenate([dwdz, K.expand_dims(dwdz[:,-1,:,:], axis=1)], axis=1)

        u = dwdy - dvdz
        v = dudz - dwdx
        w = dvdx - dudy

        j = K.stack([dudx,dudy,dudz,dvdx,dvdy,dvdz,dwdx,dwdy,dwdz], axis=-1)
        c = K.stack([u,v,w], axis=-1)

        return j, c
    else:
        if data_format == 'NCHW':
            x = nchw_to_nhwc(x)

        dudx = x[:,:,1:,0] - x[:,:,:-1,0]
        dudy = x[:,:-1,:,0] - x[:,1:,:,0] # horizontally flipped
        dvdx = x[:,:,1:,1] - x[:,:,:-1,1]
        dvdy = x[:,:-1,:,1] - x[:,1:,:,1] # horizontally flipped

        dudx = K.concatenate([dudx, K.expand_dims(dudx[:,:,-1], axis=2)], axis=2)
        dvdx = K.concatenate([dvdx, K.expand_dims(dvdx[:,:,-1], axis=2)], axis=2)
        dudy = K.concatenate([K.expand_dims(dudy[:,0,:], axis=1), dudy], axis=1)
        dvdy = K.concatenate([K.expand_dims(dvdy[:,0,:], axis=1), dvdy], axis=1)

        j = K.stack([dudx,dudy,dvdx,dvdy], axis=-1)
        w = K.expand_dims(dvdx - dudy, axis=-1) # vorticity (for visualization)

        if data_format == 'NCHW':
            j = nhwc_to_nchw(j)
            w = nhwc_to_nchw(w)

        return j, w

# --------------------------------------------------------------------------------------------------------------------------------------------------
def jacobian_with_time_layer(x, is_3d=False):
    j = None
    if is_3d:
        dudx = x[:,:,:,:,1:,0] - x[:,:,:,:,:-1,0]
        dvdx = x[:,:,:,:,1:,1] - x[:,:,:,:,:-1,1]
        dwdx = x[:,:,:,:,1:,2] - x[:,:,:,:,:-1,2]
        dudy = x[:,:,:,:-1,:,0] - x[:,:,:,1:,:,0] # horizontally flipped
        dvdy = x[:,:,:,:-1,:,1] - x[:,:,:,1:,:,1] # horizontally flipped
        dwdy = x[:,:,:,:-1,:,2] - x[:,:,:,1:,:,2] # horizontally flipped
        dudz = x[:,:,1:,:,:,0] - x[:,:,:-1,:,:,0]
        dvdz = x[:,:,1:,:,:,1] - x[:,:,:-1,:,:,1]
        dwdz = x[:,:,1:,:,:,2] - x[:,:,:-1,:,:,2]

        dudx = K.concatenate([dudx, K.expand_dims(dudx[:,:,:,:,-1], axis=4)], axis=4)
        dvdx = K.concatenate([dvdx, K.expand_dims(dvdx[:,:,:,:,-1], axis=4)], axis=4)
        dwdx = K.concatenate([dwdx, K.expand_dims(dwdx[:,:,:,:,-1], axis=4)], axis=4)

        dudy = K.concatenate([K.expand_dims(dudy[:,:,:,0,:], axis=3), dudy], axis=3)
        dvdy = K.concatenate([K.expand_dims(dvdy[:,:,:,0,:], axis=3), dvdy], axis=3)
        dwdy = K.concatenate([K.expand_dims(dwdy[:,:,:,0,:], axis=3), dwdy], axis=3)

        dudz = K.concatenate([dudz, K.expand_dims(dudz[:,:,-1,:,:], axis=2)], axis=2)
        dvdz = K.concatenate([dvdz, K.expand_dims(dvdz[:,:,-1,:,:], axis=2)], axis=2)
        dwdz = K.concatenate([dwdz, K.expand_dims(dwdz[:,:,-1,:,:], axis=2)], axis=2)

        j = K.stack([dudx,dudy,dudz,dvdx,dvdy,dvdz,dwdx,dwdy,dwdz], axis=-1)
    else:
        dudx = x[:, :,:,1:,0] - x[:, :,:,:-1,0]
        dudy = x[:, :,:-1,:,0] - x[:, :,1:,:,0] # horizontally flipped
        dvdx = x[:, :,:,1:,1] - x[:, :,:,:-1,1]
        dvdy = x[:, :,:-1,:,1] - x[:, :,1:,:,1] # horizontally flipped

        dudx = K.concatenate([dudx, K.expand_dims(dudx[:, :,:,-1], axis=3)], axis=3)
        dvdx = K.concatenate([dvdx, K.expand_dims(dvdx[:, :,:,-1], axis=3)], axis=3)
        dudy = K.concatenate([K.expand_dims(dudy[:, :,0,:], axis=2), dudy], axis=2)
        dvdy = K.concatenate([K.expand_dims(dvdy[:, :,0,:], axis=2), dvdy], axis=2)

        j = K.stack([dudx,dudy,dvdx,dvdy], axis=-1)

    return j

# --------------------------------------------------------------------------------------------------------------------------------------------------
def pf_advect(src, v, dt, mac_adv, name):
    import phi.flow as pf
    box_shape = [d-1 for d in pf.math.staticshape(v)[1:-1]]
    v_sg = pf.StaggeredGrid.from_tensors(name, pf.Box(0, box_shape), pf.unstack_staggered_tensor(v))
    if mac_adv: # e.g. self advection
        src_sg = pf.StaggeredGrid.from_tensors(name, pf.Box(0, box_shape), pf.unstack_staggered_tensor(src))
        src_advected = pf.advect.semi_lagrangian(src_sg, v_sg, dt=dt).staggered_tensor()
    else:
        box_shape = [d for d in pf.math.staticshape(src)[1:-1]]
        src_cg = pf.CenteredGrid(name, pf.Box(0, box_shape), src)
        src_advected = pf.advect.semi_lagrangian(src_cg, v_sg, dt=dt).data
    return src_advected
# --------------------------------------------------------------------------------------------------------------------------------------------------
# src_in = [d, v]
def advect(src_in, dt, mac_adv=False, name="velocity"):
    #src = np.zeros([1, 80, 64, 2])  # 2D: (batch, y, x, 2)
    #src = np.zeros([1, 80, 64, 64, 3])  # 3D  (batch, z, y, x, 3)
    # phiflow dimension order zyx <-> manta dim order xyz
    src = src_in[0]
    v = src_in[1]
    src = src[...,::-1]
    v = v[...,::-1]
    if len(src.shape) == 4: # 2D
        src = src[:,::-1]
        v = v[:,::-1]
    else: # 3D
        src = src[:,:,::-1]
        v = v[:,:,::-1]
    src_advected = pf_advect(src, v, dt, mac_adv=mac_adv, name=name)
    src_advected = src_advected[...,::-1]
    if len(src.shape) == 4: # 2D
        src_advected = src_advected[:,::-1]
    else: # 3D
        src_advected = src_advected[:,:,::-1]
    return src_advected

# --------------------------------------------------------------------------------------------------------------------------------------------------
def grad_density(x, data_format='NHWC', is_3d=False):
    print("grad_density")
    print("\tx shape: {}".format(x.shape))

    if is_3d:
        dudx = x[:,:,:,1:,3] - x[:,:,:,:-1,3]
        dudy = x[:,:,:-1,:,3] - x[:,:,1:,:,3] # horizontally flipped
        dudz = x[:,1:,:,:,3] - x[:,:-1,:,:,3]

        dudx = K.concatenate([dudx, K.expand_dims(dudx[:,:,:,-1], axis=3)], axis=3)
        dudy = K.concatenate([K.expand_dims(dudy[:,:,0,:], axis=2), dudy], axis=2)
        dudz = K.concatenate([dudz, K.expand_dims(dudz[:,-1,:,:], axis=1)], axis=1)

        j = K.stack([dudx,dudy,dudz], axis=-1)
        return j
    else:
        if data_format == 'NCHW':
            x = nchw_to_nhwc(x)

        dudx = x[:,:,1:,2] - x[:,:,:-1,2]
        dudy = x[:,:-1,:,2] - x[:,1:,:,2] # horizontally flipped

        dudx = K.concatenate([dudx, K.expand_dims(dudx[:,:,-1], axis=2)], axis=2)
        dudy = K.concatenate([K.expand_dims(dudy[:,0,:], axis=1), dudy], axis=1)

        j = K.stack([dudx,dudy], axis=-1)

        if data_format == 'NCHW':
            j = nhwc_to_nchw(j)
        return j

# --------------------------------------------------------------------------------------------------------------------------------------------------
# returns the normalized [-1,1] spherical coordinates
def cartesian_to_spherical(x, is_3d=False):
    # length of 3d vectors
    if is_3d:
        r = tf.norm(x[:,:,:,:,:3], axis=4, keepdims=True)
        phi = tf.atan2(x[:,:,:,:,1], x[:,:,:,:,0]) / math.pi
        theta = tf.acos(x[:,:,:,:,2] / r) / math.pi
        return K.concatenate([K.expand_dims(phi, axis=4), K.expand_dims(theta, axis=4)], axis=4)
    else:
        phi = tf.atan2(x[:,:,:,1], x[:,:,:,0]) / math.pi
        return K.expand_dims(phi, axis=3)

# --------------------------------------------------------------------------------------------------------------------------------------------------
def vector_length(x, is_3d=False):
    if is_3d:
        return tf.norm(x[:,:,:,:,:3], axis=4, keepdims=True)
    else:
        return tf.norm(x[:,:,:,:2], axis=3, keepdims=True)

# --------------------------------------------------------------------------------------------------------------------------------------------------
# Losses -------------------------------------------------------------------------------------------------------------------------------------------

# --------------------------------------------------------------------------------------------------------------------------------------------------
class SquaredDifferenceLoss(object):
    #---------------------------------------------------------------------------------
    def __init__(self):
        self.__name__ = "SquaredDifferenceLoss"
        print("Spawned {}".format(self.__name__))
    #---------------------------------------------------------------------------------
    def __call__(self, y_true, y_pred):
        return 0.5 * K.abs(K.sum(K.square(y_pred) - K.square(y_true))) * 0.00005

# --------------------------------------------------------------------------------------------------------------------------------------------------
class AE_Loss(object):
    #---------------------------------------------------------------------------------
    def __init__(self, sqrd_diff_loss=False, density=False, is_3d=False):
        self.__name__ = "AE_Loss"
        self.sqrd_diff_loss = sqrd_diff_loss
        self.density = density
        self.is_3d = is_3d
        print("Using density loss")
        if self.sqrd_diff_loss:
            self.sqrd_diff_loss_op = SquaredDifferenceLoss()
        print("Spawned {}".format(self.__name__))

    #---------------------------------------------------------------------------------
    def __call__(self, y_true, y_pred):
        ae_pred = y_pred
        ae_grad_pred, _ = jacobian_layer(y_pred, is_3d=self.is_3d)

        ae_true = y_true
        ae_grad_true, _ = jacobian_layer(y_true, is_3d=self.is_3d)

        velo_dim = 3 if self.is_3d else 2

        # MAE on velo
        result = losses.mean_absolute_error(ae_pred[...,0:velo_dim], ae_true[...,0:velo_dim])
        # Gradient Loss
        result += losses.mean_absolute_error(ae_grad_pred, ae_grad_true)
        if self.sqrd_diff_loss:
            result += self.sqrd_diff_loss_op(ae_true[...,0:velo_dim], ae_pred[...,0:velo_dim])
        if self.density:
            result += losses.mean_squared_error(ae_pred[...,velo_dim:velo_dim+1], ae_true[...,velo_dim:velo_dim+1])
        return result

# --------------------------------------------------------------------------------------------------------------------------------------------------
class AE_Loss_Multitile(object):
    #---------------------------------------------------------------------------------
    def __init__(self, sqrd_diff_loss=False, density=False, is_3d=False, tile_loss=False, vort_loss=False):
        self.__name__ = "AE_Loss"
        self.sqrd_diff_loss = sqrd_diff_loss
        self.vort_loss = vort_loss
        self.density = density
        self.is_3d = is_3d
        self.loss_history = [None] * 20 # adjust number if more losses arise
        self.loss_data_scale = [1.0] * 20 # adjust number if more losses arise
        self.loss_weights = [1.0] * 20
        self.loss_weights[2] = 0.5
        self.loss_weights[3] = 0.5
        print("Using density loss")
        if self.sqrd_diff_loss:
            self.sqrd_diff_loss_op = SquaredDifferenceLoss()
        print("Spawned {}".format(self.__name__))

    #---------------------------------------------------------------------------------
    def add_to_history(self, idx, val):
        if self.loss_history[idx] is not None:
            if len(self.loss_history[idx] > 100): # 100 is hyper parameter, eval different values
                self.loss_data_scale[idx] = 1.0 / statistics.mean(self.loss_history[idx])
                self.loss_history[idx] = []
            self.loss_history[idx].append(val)
        else:
            self.loss_history[idx] = [val]

    #---------------------------------------------------------------------------------
    def __call__(self, y_true, y_pred):
        ae_pred = y_pred
        ae_grad_pred, w_pred = jacobian_layer(y_pred, is_3d=self.is_3d)

        ae_true = y_true
        ae_grad_true, w_true = jacobian_layer(y_true, is_3d=self.is_3d)

        velo_dim = 3 if self.is_3d else 2

        result = 0

        # MAE on velo
        cur_result = losses.mean_absolute_error(ae_pred[...,0:velo_dim], ae_true[...,0:velo_dim])
        self.add_to_history(0, cur_result)
        result += self.loss_weights[0] * self.loss_data_scale[0] * cur_result

        # Gradient Loss
        cur_result = losses.mean_absolute_error(ae_grad_pred, ae_grad_true)
        self.add_to_history(1, cur_result)
        result += self.loss_weights[1] * self.loss_data_scale[1] * cur_result

        # # Vector length loss
        # ae_pred_length = vector_length(y_pred, is_3d=self.is_3d)
        # ae_true_length = vector_length(y_true, is_3d=self.is_3d)
        # cur_result = losses.mean_squared_error(ae_pred_length, ae_true_length)
        # self.add_to_history(2, cur_result)
        # result += self.loss_weights[2] * self.loss_data_scale[2] * cur_result

        # # Direction loss
        # ae_pred_spherical = cartesian_to_spherical(y_pred, is_3d=self.is_3d)
        # ae_true_spherical = cartesian_to_spherical(y_true, is_3d=self.is_3d)
        # cur_result = losses.mean_squared_error(ae_pred_spherical, ae_true_spherical)
        # self.add_to_history(3, cur_result)
        # result += self.loss_weights[3] * self.loss_data_scale[3] * cur_result

        if self.sqrd_diff_loss:
            cur_result = self.sqrd_diff_loss_op(ae_true[...,0:velo_dim], ae_pred[...,0:velo_dim])
            self.add_to_history(4, cur_result)
            result += self.loss_weights[4] * self.loss_data_scale[4] * cur_result
        if self.density:
            ae_grad_den_pred = grad_density(y_pred, is_3d=self.is_3d)
            ae_grad_den_true = grad_density(y_true, is_3d=self.is_3d)
            cur_result = losses.mean_absolute_error(ae_grad_den_pred, ae_grad_den_true)
            self.add_to_history(5, cur_result)
            result += self.loss_weights[5] * self.loss_data_scale[5] * cur_result
            cur_result = losses.mean_squared_error(ae_pred[...,velo_dim:velo_dim+1], ae_true[...,velo_dim:velo_dim+1])
            self.add_to_history(6, cur_result)
            result += self.loss_weights[6] * self.loss_data_scale[6] * cur_result

        if self.vort_loss:
            # Vort Loss
            cur_result = losses.mean_absolute_error(w_pred, w_true)
            self.add_to_history(7, cur_result)
            result += self.loss_weights[7] * self.loss_data_scale[7] * cur_result

        return result

# --------------------------------------------------------------------------------------------------------------------------------------------------
class GradLoss(object):
    #---------------------------------------------------------------------------------
    def __init__(self, use_mse=False, is_3d=False):
        self.__name__ = "GradLoss"
        self.use_mse = use_mse
        self.is_3d = is_3d
        print("Spawned {}".format(self.__name__))

    #---------------------------------------------------------------------------------
    def __call__(self, y_true, y_pred):
        ae_pred = y_pred
        ae_grad_pred = jacobian_with_time_layer(y_pred, is_3d=self.is_3d)

        ae_true = y_true
        ae_grad_true = jacobian_with_time_layer(y_true, is_3d=self.is_3d)

        result = 0.0
        if self.use_mse:
            result = losses.mean_squared_error(ae_pred, ae_true) + losses.mean_squared_error(ae_grad_pred, ae_grad_true)
        else:
            result = losses.mean_absolute_error(ae_pred, ae_true) + losses.mean_absolute_error(ae_grad_pred, ae_grad_true)

        return result

# --------------------------------------------------------------------------------------------------------------------------------------------------
class Pred_Loss(object):
    #---------------------------------------------------------------------------------
    def __init__(self, GT_split_idx, skip_steps=False, gradient_loss=False, sqrd_diff_loss=False, gradient_loss_mse=False, density=False, is_3d=False):
        self.__name__ = "Pred_Loss"
        self.split_idx = GT_split_idx
        self.skip_steps = skip_steps
        self.gradient_loss = gradient_loss
        self.sqrd_diff_loss = sqrd_diff_loss
        self.density = density
        self.is_3d = is_3d
        print("Using density loss")
        if self.gradient_loss:
            self.gradient_loss_op = GradLoss(gradient_loss_mse, is_3d=self.is_3d)
        if self.sqrd_diff_loss:
            self.sqrd_diff_loss_op = SquaredDifferenceLoss()

        print("Spawned {}".format(self.__name__))
    #---------------------------------------------------------------------------------
    def __call__(self, y_true, y_pred):
        y_true = y_pred[:, self.split_idx:]
        y_pred = y_pred[:, :self.split_idx]
        if self.skip_steps:
            idx = tf.random_uniform([1], 0, self.split_idx, dtype = tf.int32)[0]
        else:
            idx = 0
        loss = losses.mean_squared_error(y_true[:, idx:], y_pred[:, idx:])
        return loss
# --------------------------------------------------------------------------------------------------------------------------------------------------
class Pred_Decoded_Loss(object):
    #---------------------------------------------------------------------------------
    def __init__(self, skip_steps=False, gradient_loss=False, sqrd_diff_loss=False, gradient_loss_mse=False, density=False, is_3d=False):
        self.__name__ = "Pred_Decoded_Loss"
        self.skip_steps = skip_steps
        self.gradient_loss = gradient_loss
        self.sqrd_diff_loss = sqrd_diff_loss
        self.density = density
        self.is_3d = is_3d
        print("Using density loss")
        if self.gradient_loss:
            self.gradient_loss_op = GradLoss(gradient_loss_mse, is_3d=self.is_3d)
        if self.sqrd_diff_loss:
            self.sqrd_diff_loss_op = SquaredDifferenceLoss()

        print("Spawned {}".format(self.__name__))
    #---------------------------------------------------------------------------------
    def __call__(self, y_true, y_pred):
        if self.skip_steps:
            idx = tf.random_uniform([1], 0, y_pred.shape[1], dtype = tf.int32)[0]
        else:
            idx = 0

        velo_dim = 3 if self.is_3d else 2

        if self.gradient_loss:
            loss = self.gradient_loss_op(y_true[:, idx:, ..., 0:velo_dim], y_pred[:, idx:, ..., 0:velo_dim])
        else:
            loss = losses.mean_absolute_error(y_true[:, idx:, ..., 0:velo_dim], y_pred[:, idx:, ..., 0:velo_dim])
        if self.sqrd_diff_loss:
            loss += self.sqrd_diff_loss_op(y_true[:, idx:, ..., 0:velo_dim], y_pred[:, idx:, ..., 0:velo_dim])
        if self.density:
            loss += losses.mean_squared_error(y_true[:, idx:, ..., velo_dim:velo_dim+1], y_pred[:, idx:, ..., velo_dim:velo_dim+1])

        return loss

# --------------------------------------------------------------------------------------------------------------------------------------------------
class Split_Loss(object):
    #---------------------------------------------------------------------------------
    def __init__(self, first_idx, last_idx):
        self.__name__ = "Split_Loss"
        self.first_idx = first_idx
        self.last_idx = last_idx
        print("Spawned {}".format(self.__name__))
    #---------------------------------------------------------------------------------
    def __call__(self, y_true, y_pred):
        loss = losses.mean_absolute_error(y_true[..., self.first_idx:self.last_idx], y_pred[..., self.first_idx:self.last_idx])
        return loss

max_stars_repo_name: wiewel/LatentSpaceSubdivision

max_issues_repo_name: wiewel/LatentSpaceSubdivision

max_forks_repo_name: wiewel/LatentSpaceSubdivision


text:
import json
import logging
import time
from collections import deque
from typing import Tuple, Union

import numpy as np

MAX_TRIALS = 3
TRANSITION_DELAY = 2
REENTER_DELAY = 3
INIT_PENALTY_DELAY = 10

class Processor(object):
    def __init__(self):
        self.stages = self._load_stages()
        self.users = self._load_users()
        self.autofills = self._load_autofills()
        self.current_user = self.users["test"]

        self.current_stage = self.stages["entry"]
        self.next_stage = self._get_next_stage()
        self.delay = 2
        self.ts = time.time()

        self.configs = None
        self.custom_rule = None

        self.inputs = deque()
        self.connections = deque()
        self.remaining_trials = MAX_TRIALS
        self.penalty_delay = INIT_PENALTY_DELAY

        logging.info(" Processor initialized")
        self._logging_stage_info()

    def _load_stages(self) -> dict:
        with open('stages.json') as f:
            return json.load(f)

    def _load_users(self) -> dict:
        with open('users.json') as f:
            return json.load(f)

    def _save_users(self) -> None:
        with open('users.json', 'w') as f:
            f.write(json.dumps(self.users, indent=4))

    def _load_autofills(self) -> dict:
        with open('autofills.json') as f:
            return json.load(f)

    def _get_next_stage(self, action: str = None):
        if self.stages[self.current_stage]["multi-branching"] and action:
            return self.stages[self.current_stage]["next"][action]
        elif not self.stages[self.current_stage]["multi-branching"]:
            return self.stages[self.current_stage]["next"]
        else:
            raise NotImplementedError(
                'Trying to move to another branch without action')

    def _is_done(self):
        # Always running
        return False
        # return (self.current_stage == 'complete' or self.next_stage == 'complete')

    def _update_stage(self) -> None:
        if not self.next_stage or (time.time() - self.ts) < self.delay:
            return
        self.current_stage = self.next_stage
        self.next_stage = None
        self.ts = None

        try:
            self.custom_rule = self.stages[self.current_stage]["custom-rule"]
        except KeyError:
            self.custom_rule = False
        self.configs = self._load_config()
        input_len = self.stages[self.current_stage]["input-length"]
        self.inputs = deque(maxlen=input_len)
        if not self.remaining_trials:
            self.remaining_trials = MAX_TRIALS

        # Used in password related stages only
        self.connections = deque()

        logging.info(f" Stage is updated to {self.current_stage}")
        logging.info(f" Current stage config {self.configs}")

    def _load_config(self) -> Union[dict, None]:
        config_file = self.stages[self.current_stage]["config-file"]
        if config_file:
            with open(config_file) as f:
                return json.load(f)
        else:
            return None

    def process(self, position: list) -> Tuple[bool, bool]:
        if self.next_stage or self._is_done():
            return False, self._is_done()
        if self._is_reset_area(position):
            logging.info(" - RESET - Detected\n")
            self.inputs.clear()
            self.connections.clear()
            return True, self._is_done()

        button = self._get_button(position)

        if button and button not in self.inputs:
            # TODO: Replace above with next line to enable nonduplication inputs
            # if button:
            self.inputs.append(button)
            logging.info(f" - {button} - Detected\n" +
                         f"\tCurrent inputs: {self.inputs}")

            if self.custom_rule:
                self._handle_custom_rules()
            else:
                self._handle_common_rule()
        return True, self._is_done()

    def _is_reset_area(self, position: list) -> bool:
        if self.current_stage in ['register_input_phase_one',
                                  'register_input_phase_two',
                                  'login_input']:
            if position[0] < 0.21875 or position[0] > 0.78125:
                return True
        return False

    def _get_button(self, position: list) -> Union[str, None]:
        for button_name, config in self.configs.items():
            center, radius = config["center"], config["range"]
            if np.linalg.norm(np.array(position) - np.array(center)) < radius:
                return button_name
        return None

    def _handle_common_rule(self):
        self.next_stage = self._get_next_stage(self.inputs[0])
        self.ts = time.time()

        self._logging_stage_info()

    def _get_latest_two_inputs(self):
        return ''.join(list(self.inputs)[-2:])

    def _insert_autofill(self):
        latest_two = self._get_latest_two_inputs()
        if latest_two in self.autofills:
            autofill = self.autofills[latest_two]

            if autofill in self.inputs:
                conn1 = ''.join([autofill, self.inputs[-2], autofill])
                conn2 = ''.join([autofill, self.inputs[-1], autofill])
                if (conn1[:2] not in self.connections and conn1[1:] not in self.connections):
                    self.connections.append(conn1[:2])
                if (conn2[:2] not in self.connections and conn2[1:] not in self.connections):
                    self.connections.append(conn2[:2])
                return

            logging.info(f" Current input: {self.inputs}\n" +
                         f" Autofilling - {autofill} - ...")
            if len(self.inputs) == self.inputs.maxlen:
                self.inputs.pop()
                self.inputs.append(autofill)
                self.connections.append(self._get_latest_two_inputs())
            else:
                latest_one = self.inputs.pop()
                self.inputs.append(autofill)
                self.connections.append(self._get_latest_two_inputs())
                self.inputs.append(latest_one)
                self.connections.append(self._get_latest_two_inputs())

            logging.info(" Autofill complete\n" +
                         f" Current inputs: {self.inputs}")
        else:
            self.connections.append(self._get_latest_two_inputs())

    def _handle_custom_rules(self) -> None:
        if self.current_stage == 'register_input_phase_one':
            if len(self.inputs) > 1:
                self._insert_autofill()
            if len(self.inputs) == self.inputs.maxlen:
                self.current_user["password"] = ''.join(self.inputs)
                self.next_stage = self._get_next_stage()
                self.ts = time.time()
                self.delay = 5

                self._logging_stage_info()
        elif self.current_stage == 'register_input_phase_two':
            if len(self.inputs) > 1:
                self._insert_autofill()
            if len(self.inputs) == self.inputs.maxlen:
                match = self.current_user["password"] == ''.join(self.inputs)
                if match:
                    self._save_users()
                self.next_stage = self._get_next_stage(str(match))
                self.ts = time.time()

                self._logging_stage_info()
        elif self.current_stage == 'login_input':
            if len(self.inputs) > 1:
                self._insert_autofill()
            if len(self.inputs) == self.inputs.maxlen:
                match = self.current_user["password"] == ''.join(self.inputs)
                if not match:
                    self.remaining_trials -= 1
                    if not self.remaining_trials:
                        self.delay = self.penalty_delay
                        self.penalty_delay *= 2
                    else:
                        self.delay = REENTER_DELAY
                else:
                    self.delay = TRANSITION_DELAY
                self.next_stage = self._get_next_stage(str(match))
                self.ts = time.time()

                self._logging_stage_info()
        else:
            logging.error(" Custom rule not specified")

    def get_stage_info(self) -> dict:
        self._update_stage()
        transition_delay = 0 if not self.ts else (self.ts + self.delay -
                                                  time.time())

        return {
            "currentStage": self.current_stage,
            "nextStage": self.next_stage,
            "transitionDelay": transition_delay,
            "params": {
                "buttons": list(self.inputs),
                "links": list(self.connections),
                "remainingTrials": self.remaining_trials,
                "nextPenaltyTime": self.penalty_delay
            }
        }

    def _logging_stage_info(self) -> None:
        logging.info(f" Changing stage: \n" +
                     f"\tfrom: {self.current_stage}\n" +
                     f"\tto: {self.next_stage}\n" +
                     f"\tafter: {self.delay}s\n")

ext: py

reason: import numpy


text:
#!/usr/bin/python
"""
FlowCal Python API example, using an input Excel UI file.

This script is divided in two parts. Part one uses FlowCal's Excel UI
to process calibration beads data and cell samples according to an input
Excel file. The exact operations performed are identical to when normally
using the Excel UI. However, instead of generating an output Excel file, a
list of objects representing gated and transformed flow cytometry samples
is obtained.

Part two exemplifies how to use the processed cell sample data with
FlowCal's plotting and statistics modules, in order to produce interesting
plots.

For details about the experiment, samples, and instrument used, please
consult readme.txt.

"""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import FlowCal

if __name__ == "__main__":

    ###
    # Part 1: Obtaining processed flow cytometry data from an input Excel file
    ###

    # ``FlowCal.excel_ui.read_table() reads a table from an Excel file, and
    # returns the contents as a pandas dataframe.
    # We need tables describing the instruments used, the calibration beads
    # files, and the cell samples.
    instruments_table = FlowCal.excel_ui.read_table(filename='experiment.xlsx',
                                                    sheetname='Instruments',
                                                    index_col='ID')
    beads_table = FlowCal.excel_ui.read_table(filename='experiment.xlsx',
                                              sheetname='Beads',
                                              index_col='ID')
    samples_table = FlowCal.excel_ui.read_table(filename='experiment.xlsx',
                                                sheetname='Samples',
                                                index_col='ID')

    # Process beads samples
    # ``FlowCal.excel_ui.process_beads_table()`` reads a properly formatted
    # table describing calibration beads samples (``beads_table``), performs
    # density gating as specified in this table, and generates MEF
    # transformation functions that can later be applied to cell samples.
    # To do so, it requires a table describing the flow cytometer used
    # (``instruments_table``). Here, we also use verbose mode, and indicate that
    # plots describing individual steps should be generated in the folder
    # "plot_beads". The result is a list of ``FCSData`` objects representing
    # gated and transformed calibration beads samples (``beads_samples``), and
    # a dictionary containing MEF transformation functions
    # (``mef_transform_fxns``). This will be used later to process cell samples.
    beads_samples, mef_transform_fxns = FlowCal.excel_ui.process_beads_table(
        beads_table=beads_table,
        instruments_table=instruments_table,
        verbose=True,
        plot=True,
        plot_dir='plot_beads')

    # Process cell samples
    # ``FlowCal.excel_ui.process_samples_table()`` reads a properly formatted
    # table describing cell samples (``samples_table``), and performs density
    # gating and unit transformations as specified in this table. To do so, it
    # requires a table describing the flow cytometer used
    # (``instruments_table``) and a corresponding dictionary with MEF
    # transformation functions (``mef_transform_fxns``). Here, we also use
    # verbose mode, and indicate that plots of each sample should be generated
    # in the folder "plot_samples". The result is a list of ``FCSData`` objects
    # (``samples``) representing gated and transformed flow cytometry cell
    # samples.
    samples = FlowCal.excel_ui.process_samples_table(
        samples_table=samples_table,
        instruments_table=instruments_table,
        mef_transform_fxns=mef_transform_fxns,
        verbose=True,
        plot=True,
        plot_dir='plot_samples')

    ###
    # Part 2: Examples on how to use processed cell sample data
    ###

    # We will read IPTG concentrations for each sample from the Excel file.
    # ``samples_table`` contains all the data from sheet "Samples", including
    # data not directly used by ``FlowCal.excel_ui.process_samples_table()``.
    # ``samples_table`` is a pandas dataframe, with each column having the
    # same name as the corresponding header in the Excel file.
    # We multiply the IPTG concentration by 1000 to convert to micromolar.
    iptg = samples_table['IPTG (mM)']*1000

    # Histogram of all samples
    # Here, we plot the fluorescence histograms of all five samples in the same
    # figure, using ``FlowCal.plot.hist1d``. Note how this function can be used
    # in the context of accessory matplotlib functions to modify the axes
    # limits and labels and add a legend, among others.
    plt.figure(figsize=(6,3.5))
    FlowCal.plot.hist1d(samples,
                        channel='FL1',
                        histtype='step',
                        bins=128)
    plt.ylim([0, 2000])
    plt.xlabel('FL1  (Molecules of Equivalent Fluorescein, MEFL)')
    plt.legend(['{} $\mu M$ IPTG'.format(i) for i in iptg],
               loc='upper left',
               fontsize='small')
    plt.tight_layout()
    plt.savefig('histograms.png', dpi=200)
    plt.close()

    # Here we illustrate how to obtain statistics from the fluorescence of each
    # sample, and how to use them in a plot.
    # The stats module contains functions to calculate different statistics
    # such as mean, median, and standard deviation. Here, we calculate the
    # geometric mean from channel FL1 of each sample, and plot them against the
    # corresponding IPTG concentrations.
    samples_fluorescence = [FlowCal.stats.gmean(s, channels='FL1')
                            for s in samples]
    plt.figure(figsize=(5.5, 3.5))
    plt.plot(iptg,
             samples_fluorescence,
             marker='o',
             color=(0, 0.4, 0.7))
    plt.xlabel('IPTG Concentration ($\mu M$)')
    plt.ylabel('FL1 Fluorescence (MEFL)')
    plt.tight_layout()
    plt.savefig('dose_response.png', dpi=200)
    plt.close()

    print("\nDone.")

max_stars_repo_name: castillohair/FlowCal

max_issues_repo_name: castillohair/FlowCal

max_forks_repo_name: castillohair/FlowCal


text:
#Author
#Felice De Luca
#https://github.com/felicedeluca

import sys
import os

import networkx as nx
from networkx.drawing.nx_agraph import write_dot
from networkx.drawing.nx_agraph import read_dot as nx_read_dot

import math

import numpy as np # to compute inner product

def compute_edges_vects(G):
    '''
    Computes the vectors for each edge in G and resturns them in a list.
    '''

    vects_list = []

    pos_dict = nx.get_node_attributes(G, 'pos')

    for e in nx.edges(G):

        (source, target) = e

        x_source = float(pos_dict[source].split(",")[0])
        x_target = float(pos_dict[target].split(",")[0])

        y_source = float(pos_dict[source].split(",")[1])
        y_target = float(pos_dict[target].split(",")[1])

        # The following check is to force the upwardness
        # (if the edges vects) are inverted
        if y_target < y_source:
            x_temp = x_target
            x_target = x_source
            x_source = x_temp

            y_temp = y_target
            y_target = y_source
            y_source = y_temp

        len = [x_source - x_target, y_source - y_target]
        norm = math.sqrt(len[0] ** 2 + len[1] ** 2)
        dir = [len[0]/norm, len[1]/norm]
        vect = [dir[0], dir[1]]

        vects_list.append(vect)

    return vects_list

def compute_upwardflow(G):
    '''
    Computes the Upward flow of the given upward drawing of the graph G.

    This metric determines the proportion of edge segments of a drawing which have a consistent direction.
    Edge segments are used rather than edges, as edges with edge segments of alternating direction are generally considered undesirable.

    The desired direction is usually upwards or downwards with respect to a vertical axis, but the metric makes no assumptions about its orientation. It is assumed that G is a directed graph.

    The notation ei indicates the vector corresponding to the ith directed edge in the drawing:The inner product of two vectors is denoted <v1, v2> : The unit vector parallel with the desired direction is denoted 1 and is considered to point in the direction of desired flow.

    This metric is 0 for non directed graphs.

    <tt>https://pdfs.semanticscholar.org/be7e/4c447ea27e0891397ae36d8957d3cbcea613.pdf</tt>

    '''

    e_vects = compute_edges_vects(G)

    prods_sum = 0

    for e in e_vects:
        ones = [1, 1]
        inn_prod = np.inner(e, ones)

        if inn_prod > 0:
            prods_sum += 1

    flow = prods_sum/len(e_vects)
    flow = round(flow, 3)

    return flow

reason: import numpy,import networkx,from networkx


text:
# coding: utf-8

# In[ ]:

# Awesome image patch finder,
# Due to http://stackoverflow.com/questions/9525313

import os
import sys
import numpy
import random

import numpy as np
import scipy.ndimage as ndimage
import scipy.spatial as spatial
import scipy.misc as misc
import matplotlib.pyplot as plt
import matplotlib.patches as patches

import cv2
class BBox(object):
    def __init__(self, x1, y1, x2, y2):
        '''
        (x1, y1) is the upper left corner,
        (x2, y2) is the lower right corner,
        with (0, 0) being in the upper left corner.
        '''
        if x1 > x2: x1, x2 = x2, x1
        if y1 > y2: y1, y2 = y2, y1
        self.x1 = x1
        self.y1 = y1
        self.x2 = x2
        self.y2 = y2
    def taxicab_diagonal(self):
        '''
        Return the taxicab distance from (x1,y1) to (x2,y2)
        '''
        return self.x2 - self.x1 + self.y2 - self.y1
    def overlaps(self, other):
        '''
        Return True iff self and other overlap.
        '''
        return not ((self.x1 > other.x2)
                    or (self.x2 < other.x1)
                    or (self.y1 > other.y2)
                    or (self.y2 < other.y1))
    def __eq__(self, other):
        return (self.x1 == other.x1
                and self.y1 == other.y1
                and self.x2 == other.x2
                and self.y2 == other.y2)

def find_paws(data, smooth_radius = 5, threshold = 0.0001):
    # http://stackoverflow.com/questions/4087919/how-can-i-improve-my-paw-detection
    """Detects and isolates contiguous regions in the input array"""
    # Blur the input data a bit so the paws have a continous footprint
    data = ndimage.uniform_filter(data, smooth_radius)
    # Threshold the blurred data (this needs to be a bit > 0 due to the blur)
    thresh = data > threshold
    # Fill any interior holes in the paws to get cleaner regions...
    filled = ndimage.morphology.binary_fill_holes(thresh)
    # Label each contiguous paw
    coded_paws, num_paws = ndimage.label(filled)
    # Isolate the extent of each paw
    # find_objects returns a list of 2-tuples: (slice(...), slice(...))
    # which represents a rectangular box around the object
    data_slices = ndimage.find_objects(coded_paws)
    return data_slices

def slice_to_bbox(slices):
    for s in slices:
        dy, dx = s[:2]
        yield BBox(dx.start, dy.start, dx.stop+1, dy.stop+1)

def remove_overlaps(bboxes):
    '''
    Return a set of BBoxes which contain the given BBoxes.
    When two BBoxes overlap, replace both with the minimal BBox that contains both.
    '''
    # list upper left and lower right corners of the Bboxes
    corners = []

    # list upper left corners of the Bboxes
    ulcorners = []

    # dict mapping corners to Bboxes.
    bbox_map = {}

    for bbox in bboxes:
        ul = (bbox.x1, bbox.y1)
        lr = (bbox.x2, bbox.y2)
        bbox_map[ul] = bbox
        bbox_map[lr] = bbox
        ulcorners.append(ul)
        corners.append(ul)
        corners.append(lr)

    # Use a KDTree so we can find corners that are nearby efficiently.
    tree = spatial.KDTree(corners)
    new_corners = []
    for corner in ulcorners:
        bbox = bbox_map[corner]
        # Find all points which are within a taxicab distance of corner
        indices = tree.query_ball_point(
            corner, bbox_map[corner].taxicab_diagonal(), p = 1)
        for near_corner in tree.data[indices]:
            near_bbox = bbox_map[tuple(near_corner)]
            if bbox != near_bbox and bbox.overlaps(near_bbox):
                # Expand both bboxes.
                # Since we mutate the bbox, all references to this bbox in
                # bbox_map are updated simultaneously.
                bbox.x1 = near_bbox.x1 = min(bbox.x1, near_bbox.x1)
                bbox.y1 = near_bbox.y1 = min(bbox.y1, near_bbox.y1)
                bbox.x2 = near_bbox.x2 = max(bbox.x2, near_bbox.x2)
                bbox.y2 = near_bbox.y2 = max(bbox.y2, near_bbox.y2)
    return set(bbox_map.values())

if False:
    fig = plt.figure()
    ax = fig.add_subplot(111)

    data = misc.imread('/data/vision/torralba/gigaSUN/www/unit_annotation/result_segments_iterations/caffeNet_imagenet2places_iter_1/html/image/conv5-0000.jpg')
    data = data[:160,3:163,:]

    im = ax.imshow(data)
    data_slices = find_paws(data, smooth_radius = 20, threshold = 52)

    bboxes = remove_overlaps(slice_to_bbox(data_slices))
    for bbox in bboxes:
        xwidth = bbox.x2 - bbox.x1
        ywidth = bbox.y2 - bbox.y1
        p = patches.Rectangle((bbox.x1, bbox.y1), xwidth, ywidth,
                              fc = 'none', ec = 'red')
        ax.add_patch(p)

    plt.show()

def input_image_filename(basename, iter, zunit):
#    return ('/data/vision/torralba/gigaSUN/www/unit_annotation/result_segments_iterations' +
    return ('/data/vision/torralba/scratch2/davidbau/iccv' +
        '/%s_iter_%d/html/image/conv5-%04d.jpg' % (basename, iter, zunit))

def output_image_filename(basename, iter, zunit):
#    return ('/data/vision/torralba/gigaSUN/www/unit_annotation/result_segments_iterations' +
    return ('/data/vision/torralba/scratch2/davidbau/iccv' +
        '/%s_iter_%d/html/image/conv5-%04d_crop.jpg' % (basename, iter, zunit))

def biggest_square_bbox(data, smooth_radius, threshold):
    data_slices = find_paws(data, smooth_radius=smooth_radius, threshold=threshold)
    if not data_slices:
        return BBox(0,0,0,0)
    bboxes = remove_overlaps(slice_to_bbox(data_slices))
    largest = max(bboxes, key=lambda b: (b.x2 - b.x1) * (b.y2 - b.y1))
    bwidth = largest.x2 - largest.x1 - 1
    bheight = largest.y2 - largest.y1 - 1

    if bwidth > bheight:
        # Height needs to expand, but keep it legal
        largest.y1 = min(data.shape[1] - bwidth, max(0,
                (largest.y1 + largest.y2 - bwidth) // 2))
        largest.y2 = largest.y1 + bwidth + 1
    else:
        # Width needs to expand, but keep it legal
        largest.x1 = min(data.shape[0] - bheight, max(0,
                (largest.x1 + largest.x2 - bheight) // 2))
        largest.x2 = largest.x1 + bheight + 1
    return largest

def best_tightcrop(imagedata):
    bbox = biggest_square_bbox(imagedata, smooth_radius = 10, threshold = 50)
    # print bbox.x1, bbox.x2, bbox.y1, bbox.y2
    # Don't zoom in to too tiny an area - just return the original
    if bbox.x2 - bbox.x1 < 10:
        return imagedata
    return cv2.resize(src=imagedata[bbox.y1:bbox.y2, bbox.x1:bbox.x2, :], dsize=imagedata.shape[:2])

def crop_tiled_image(data, border=3):
    outdata = data.copy()
    width = data.shape[0]
    for x in range(0, data.shape[1], width + border):
        onesquare = data[:,x:x+width,:]
        outdata[:,x:x+width,:] = best_tightcrop(onesquare)
    return outdata

def process_iteration(basename, iteration, show_things=False):
    for zunit in range(256):
        input_image = input_image_filename(basename, iteration, zunit)
        output_image = output_image_filename(basename, iteration, zunit)
        if os.path.exists(output_image)==False and os.path.exists(input_image)==True:
            print('processing ' + input_image)
            data = cv2.imread(input_image)
            cropped = crop_tiled_image(data)
            cv2.imwrite(output_image, cropped)
            if show_things:
                plt.imshow(data)
                plt.show()
                plt.imshow(crop_tiled_image(data))
                plt.show()

def process_snapshot():
    name = '/data/vision/oliva/scenedataset/modelzoo/list_iterations.txt'
    with open(name) as f:
        lines = f.readlines()
    for line in lines:
        print(line)
        model_name = line.rstrip()
        items = model_name.split('/')
        model_name = items[1]
        items = model_name.split('.')
        model_name = items[0]
        items = model_name.split('_')
        basename = '_'.join(items[:-2])
        iteration = int(items[-1])
        print(basename, iteration)
        process_iteration(basename, iteration)

def process_snapshot_iterations():
    basename = 'places'
    iters = [1,2,4,9,20,44,99,223,492,1108,2446,5509,12164,27396,60491,136238,300818,600818,1200818,2400818]
    for iteration in iters:
        process_iteration(basename, iteration)

max_stars_repo_name: zhoubolei/cnnvisualizer

max_issues_repo_name: Pandinosaurus/cnnvisualizer

max_forks_repo_name: Pandinosaurus/cnnvisualizer

reason: import numpy,import scipy


text:
import numpy as np
from scipy.optimize import linear_sum_assignment

def search_assignment(matrix,
                      row_assignment=False,
                      maximize=True,
                      inplace=False):
  r"""Solve the linear sum assignment problem.

  This function can also solve a generalization of the classic assignment
  problem where the cost matrix is rectangular. If it has more rows than
  columns, then not every row needs to be assigned to a column, and vice
  versa.

  Arguments:
    matrix : an Array.
      The cost matrix of the bipartite graph.
    inplace : a Boolean.
      If True, return a new matrix with the applied assignment

  Returns:
    row_ind, col_ind : array
        An array of row indices and one of corresponding column indices giving
        the optimal assignment. The cost of the assignment can be computed
        as ``cost_matrix[row_ind, col_ind].sum()``. The row indices will be
        sorted; in the case of a square cost matrix they will be equal to
        ``numpy.arange(cost_matrix.shape[0])``.
  """
  assert matrix.ndim == 2
  rows, cols = linear_sum_assignment(matrix if row_assignment else matrix.T,
                                     maximize=maximize)
  # select the right assignment and keep the order of other dimension intact
  if not row_assignment:
    rows, cols = (cols, rows)
    ids = np.argsort(rows)
  else:
    ids = np.argsort(cols)
  rows = rows[ids]
  cols = cols[ids]
  # inplace output matrix
  if inplace:
    # make sure all indices appear
    rows = rows.tolist()
    for i in range(matrix.shape[0]):
      if i not in rows:
        rows.append(i)
    cols = cols.tolist()
    for i in range(matrix.shape[0]):
      if i not in cols:
        cols.append(i)
    return matrix[rows] if row_assignment else matrix[:, cols]
  return rows, cols

max_stars_repo_name: tirkarthi/odin-ai

max_issues_repo_name: tirkarthi/odin-ai

max_forks_repo_name: tirkarthi/odin-ai

reason: import numpy,from scipy


text:
from imutils.video import VideoStream
import numpy as np
import imagezmq
import argparse
import socket
import time
import cv2
import threading

# construct the argument parser and parse the arguments
ap = argparse.ArgumentParser()

# ip = '192.168.1.143'
ip = 'localhost'
print('d')
sender = imagezmq.ImageSender(connect_to="tcp://{}:5555".format(ip))
thread_sender = imagezmq.ImageSender(connect_to="tcp://{}:5556".format(ip))

thread_frame = frame = None

rpiName = socket.gethostname()
# vs = VideoStream(usePiCamera=True).start()
# vs = VideoStream(src=0).start()
cap = cv2.VideoCapture(0)
time.sleep(2.0)
i = 0

# cap.set(cv2.CAP_PROP_FRAME_WIDTH, int(1024))
# cap.set(cv2.CAP_PROP_FRAME_HEIGHT, int(768))

def thread_send():
    global thread_sender, rpiName, thread_frame
    thread_sender.send_image(rpiName, thread_frame)

while True:
    # read the frame from the camera and send it to the server
    ret, frame = cap.read()
    frame = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)

    thread_frame, frame = np.array_split(frame, 2)

    # i+=1
    th = threading.Thread(target=thread_send)
    th.start()
    sender.send_image(rpiName, frame)
    th.join()
    # print(i)

max_stars_repo_name: lTheMightyMl/imagezmq

max_issues_repo_name: lTheMightyMl/imagezmq

max_forks_repo_name: lTheMightyMl/imagezmq


text:
#!/usr/bin/env python

"""General demonstrations /w set start and goal
"""
import os
import numpy as np
import yaml

# pylint: disable=invalid-name
# pylint: disable=import-error
# pylint: disable=wrong-import-position
parentdir = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
traj_path = os.path.join(parentdir, 'scripts')
os.sys.path.insert(0, traj_path)

from trajectory_solver import TrajectorySolver
from agent import Agent
# pylint: enable=invalid-name
# pylint: enable=import-error
# pylint: enable=wrong-import-position

# Read arguments from yaml file
PARENT_DIR = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
FILE_PATH = os.path.join(PARENT_DIR, 'swarm_manager/conf/swarm_conf.yaml')

with open(FILE_PATH) as f:
    YAML_CONF = yaml.load(f, Loader=yaml.FullLoader)

SOLVER_ARGS = YAML_CONF['trajectory_solver']

AGENT_ARGS = {'r_min': SOLVER_ARGS['r_min'],
              'col_radius_ratio': SOLVER_ARGS['col_radius_ratio'],
              'goal_dist_thres': SOLVER_ARGS['goal_dist_thres'],
              'goal_speed_thres': SOLVER_ARGS['goal_speed_thres'],
             }

def demo_two_agents():
    """Two agents trading spots, horizontal
    """
    a_1 = Agent(AGENT_ARGS, start_pos=[0.5, 2.0, 0.0], goal=[4.0, 2.0, 0.0])

    a_2 = Agent(AGENT_ARGS, start_pos=[4.0, 2.0, 0.0], goal=[0.5, 2.0, 0.0])

    return [a_1, a_2]

def demo_two_agents_vert():
    """Two agents trading spots, vertically
    """
    a_1 = Agent(AGENT_ARGS, start_pos=[2.5, 1.0, 0.0], goal=[2.5, 4.0, 0.0])

    a_2 = Agent(AGENT_ARGS, start_pos=[2.5, 4.0, 0.0], goal=[2.5, 1.0, 0.0])

    return [a_1, a_2]

def demo_wall(wall_coords=None):
    """Demo with a wall

    Args:
        wall_coords (list of tuple): [(x, y, z), (x, y, z)]

    Returns:
        list: [[agent_list], [wall_coords]]
    """
    a_1 = Agent(AGENT_ARGS, start_pos=[0.0, 2.0, 0.0], goal=[4.0, 2.0, 0.0])

    obs_coords = compute_obstacle([wall_coords], 15)

    return [a_1], obs_coords

def through_wall():
    """Wall"""
    a_1 = Agent(AGENT_ARGS, start_pos=[0.0, 2.0, 0.0], goal=[4.0, 2.0, 0.0])

    a_2 = Agent(AGENT_ARGS, start_pos=[5.0, 2.0, 0.0], goal=[0.0, 2.0, 0.0])

    obs_coords = compute_obstacle([[(2.0, -1.0, 0.), (2.0, 1.5, 0.0)],
                                   [(2.0, 2.5, 0.), (2.0, 5.0, 0.0)]], 30)

    return [a_1, a_2], obs_coords

def corners_2():
    """Two agents trading spots, starting from opposite corners
    """
    # a_1 = Agent(AGENT_ARGS, start_pos=[0.0, 4.0, 0.0], goal=[4.0, 0.0, 0.0])
    # a_2 = Agent(AGENT_ARGS, start_pos=[4.0, 0.0, 0.0], goal=[0.0, 4.0, 0.0])

    a_1 = Agent(AGENT_ARGS, start_pos=[0.0, 0.0, 0.0], goal=[4.0, 4.0, 0.0])
    a_2 = Agent(AGENT_ARGS, start_pos=[4.0, 4.0, 0.0], goal=[0.0, 0.0, 0.0])

    return [a_1, a_2]

def corners_2_2():
    """Corners 2, v2
    """
    a_1 = Agent(AGENT_ARGS, start_pos=[0.0, 4.0, 0.0], goal=[4.0, 0.0, 0.0])
    a_2 = Agent(AGENT_ARGS, start_pos=[4.0, 0.0, 0.0], goal=[0.0, 4.0, 0.0])

    return [a_1, a_2]

def corners_4():
    """Four agents, starting from opposite corners
    """
    a_1 = Agent(AGENT_ARGS, start_pos=[0.0, 4.0, 0.0], goal=[4.0, 0.0, 0.0])
    a_2 = Agent(AGENT_ARGS, start_pos=[4.0, 0.0, 0.0], goal=[0.0, 4.0, 0.0])
    a_3 = Agent(AGENT_ARGS, start_pos=[0.0, 0.0, 0.0], goal=[4.0, 4.0, 0.0])
    a_4 = Agent(AGENT_ARGS, start_pos=[4.0, 4.0, 0.0], goal=[0.0, 0.0, 0.0])

    return [a_1, a_2, a_3, a_4]

def six_agents():
    """Six agents
    """
    a_1 = Agent(AGENT_ARGS, start_pos=[0.0, 0.0, 0.0], goal=[1.5, 3.0, 0.0])
    a_2 = Agent(AGENT_ARGS, start_pos=[2.0, 0.0, 0.0], goal=[0.0, 4.0, 0.0])
    a_3 = Agent(AGENT_ARGS, start_pos=[1.0, 2.5, 0.0], goal=[4.0, 0.0, 0.0])
    a_4 = Agent(AGENT_ARGS, start_pos=[4.0, 4.0, 0.0], goal=[0.0, 0.0, 0.0])
    a_5 = Agent(AGENT_ARGS, start_pos=[2.5, 2.5, 0.0], goal=[4.0, 2.5, 0.0])
    a_6 = Agent(AGENT_ARGS, start_pos=[3.2, 3.2, 0.0], goal=[0.5, 0.0, 0.0])

    return [a_1, a_2, a_3, a_4, a_5, a_6]

def corners_6():
    """Six agents, starting from opposite corners
    """
    a_1 = Agent(AGENT_ARGS, start_pos=[0.0, 4.0, 0.0], goal=[4.0, 0.0, 0.0])
    a_2 = Agent(AGENT_ARGS, start_pos=[4.0, 0.0, 0.0], goal=[0.0, 4.0, 0.0])
    a_3 = Agent(AGENT_ARGS, start_pos=[0.0, 0.0, 0.0], goal=[4.0, 4.0, 0.0])
    a_4 = Agent(AGENT_ARGS, start_pos=[4.0, 4.0, 0.0], goal=[0.0, 0.0, 0.0])
    a_5 = Agent(AGENT_ARGS, start_pos=[2.0, 0.0, 0.0], goal=[2.0, 4.0, 0.0])
    a_6 = Agent(AGENT_ARGS, start_pos=[2.0, 4.0, 0.0], goal=[2.0, 0.0, 0.0])

    return [a_1, a_2, a_3, a_4, a_5, a_6]

def seven_agents():
    """Seven agents
    """
    a_1 = Agent(AGENT_ARGS, start_pos=[0.0, 0.0, 0.0], goal=[1.5, 3.0, 0.0])
    a_2 = Agent(AGENT_ARGS, start_pos=[2.0, 0.0, 0.0], goal=[0.0, 4.0, 0.0])
    a_3 = Agent(AGENT_ARGS, start_pos=[1.0, 2.5, 0.0], goal=[4.0, 0.0, 0.0])
    a_4 = Agent(AGENT_ARGS, start_pos=[4.0, 4.0, 0.0], goal=[0.0, 0.0, 0.0])
    a_5 = Agent(AGENT_ARGS, start_pos=[2.5, 2.5, 0.0], goal=[4.0, 2.5, 0.0])
    a_6 = Agent(AGENT_ARGS, start_pos=[3.2, 3.2, 0.0], goal=[0.5, 0.0, 0.0])

    a_7 = Agent(AGENT_ARGS, start_pos=[1.7, 0.5, 0.0], goal=[0.8, 3.8, 0.0])

    return [a_1, a_2, a_3, a_4, a_5, a_6, a_7]

def update_test():
    """To test when agents position are changed after solver initialization
    """
    a_list = ['a1', 'a2', 'a3']
    agents = {}
    for each_a in a_list:
        agents[each_a] = Agent(AGENT_ARGS, )

    agent_list = [agent for (_, agent) in agents.items()]

    solver = TrajectorySolver(agent_list)

    agents['a1'].set_starting_position([1.0, 1.0, 0.0])
    agents['a1'].set_goal([4.0, 4.0, 0.0])

    agents['a2'].set_starting_position([4.0, 4.0, 0.0])
    agents['a2'].set_goal([1.0, 1.0, 0.0])

    agents['a3'].set_starting_position([2.0, 0.5, 0.0])
    agents['a3'].set_goal([2.0, 4.0, 0.0])

    solver.update_agents_info()

    solver.solve_trajectories()
    solver.plot_trajectories()

def compute_obstacle(positions, n_pts):
    """Compute coordinates of a wall

    Args:
        position (list of list of list): [[obstacle1_start, obstacle1_end], [obstacle2]]; [x, y, z]
        n_pts (int): Number of pts

    Returns:
        list: All coords of wall
    """
    all_coords = []
    for each_obstacle in positions:
        obstacle_start = each_obstacle[0]
        obstacle_end = each_obstacle[1]

        coords = []
        obstacle_x = np.linspace(obstacle_start[0], obstacle_end[0], num=n_pts)
        obstacle_y = np.linspace(obstacle_start[1], obstacle_end[1], num=n_pts)

        for (x_coord, y_coord) in zip(obstacle_x, obstacle_y):
            coords.append([x_coord, y_coord, 0])

        all_coords.append(coords)

    return all_coords


text:
"""Test game reader"""
import functools
import io
import warnings

import numpy as np
import pytest

from gameanalysis import canongame
from gameanalysis import gambit
from gameanalysis import gamegen
from gameanalysis import gamereader
from gameanalysis import learning
from gameanalysis import matgame
from gameanalysis import rsgame

@functools.lru_cache()
def egame():
    """Empty game"""
    return rsgame.empty([3, 4], [4, 3])

@functools.lru_cache()
def game():
    """Base game"""
    return gamegen.gen_profiles(egame(), 0.5)

def sgame():
    """Sample game"""
    return gamegen.gen_noise(game())

@functools.lru_cache()
def agg():
    """Action graph game"""
    return gamegen.normal_aggfn([3, 4], [4, 3], 10)

def mat():
    """Matrix game"""
    return matgame.matgame(np.random.random((4, 3, 2, 3)))

@functools.lru_cache()
def rbf():
    """Rbf learning game"""
    with warnings.catch_warnings():
        warnings.filterwarnings(
            'ignore',
            'some lengths were at their bounds, this may indicate a poor fit',
            UserWarning)
        return learning.rbfgame_train(game())

def point():
    """Point learning game"""
    return learning.point(rbf())

def sample():
    """Sample learning game"""
    return learning.sample(rbf())

def neighbor():
    """Neighbor learning game"""
    return learning.neighbor(rbf())

def add():
    """Add game"""
    return canon() + agg()

def canon():
    """Canon game"""
    return canongame.canon(game())

def const():
    """Constant game"""
    return rsgame.const_replace(egame(), 0)

@pytest.mark.parametrize('base', [
    egame, game, sgame, agg, mat, rbf, point, sample, neighbor, add, canon,
    const, 'gambit'])
def test_automatic_deserialization(base):
    """Test that we can serialize and deserialize arbitrary games"""
    if base == 'gambit':
        base = mat()
        string = gambit.dumps(base)
        fil = io.StringIO(string)
    else:
        base = base()
        string = gamereader.dumps(base)
        fil = io.StringIO()
        gamereader.dump(base, fil)
        fil.seek(0)

    copy = gamereader.loads(string)
    assert base == copy

    copy = gamereader.load(fil)
    assert base == copy

def test_parse_fail():
    """Test invalid games"""
    with pytest.raises(ValueError):
        gamereader.loads('')
    with pytest.raises(ValueError):
        gamereader.loadj({'type': 'unknown.0'})

max_stars_repo_name: calebho/gameanalysis

max_issues_repo_name: calebho/gameanalysis

max_forks_repo_name: calebho/gameanalysis


text:
import numpy.random as npr
import numpy as np
from hypergrad.exact_rep import BitStore
npr.seed(1)

def test_pop_empty():
    store = BitStore(4)
    popempty = store.pop(3)
    assert all(popempty == [0, 0, 0, 0])
    store.push(popempty, 3)
    store.push([1, 4, 0, 5], 6)
    assert all(store.pop(6) == [1, 4, 0, 5])

def test_long_sequence():
    N_iters = 200
    vect_length = 10
    Ns = []
    Ms = []
    for i in range(N_iters):
        Ms.append(npr.randint(200) + 1)
        Ns.append(npr.randint(Ms[-1], size=vect_length))
    store = BitStore(vect_length)
    coinflips = npr.rand(N_iters)
    new_Ns = []
    for N, M, r in zip(Ns, Ms, coinflips):
        if r < 0.75:
            store.push(N, M)
        else:
            new_Ns.append(store.pop(M))

    for N, M, r in zip(Ns, Ms, coinflips)[::-1]:
        if r < 0.75:
            cur_N = store.pop(M)
            assert np.all(cur_N == N)
        else:
            store.push(new_Ns.pop(), M)

max_stars_repo_name: bigaidream-projects/drmad

max_issues_repo_name: LinZichuan/drmad

max_forks_repo_name: LinZichuan/drmad


text:
# -*- coding: utf-8 -*-
from types import SimpleNamespace
import numpy as np
from pyyeti import expmint

class SolveExp1(object):
    """
    1st order ODE time domain solver based on the matrix exponential.

    This class is for solving: ``yd - A y = f``

    This solver is exact assuming either piece-wise linear or
    piece-wise constant forces. See :class:`SolveExp2` for more
    details on how this algorithm solves the ODE.

    Examples
    --------
    Calculate impulse response of state-space system::

        xd = A @ x + B @ u
        y  = C @ x + D @ u

    where:
        - force = 0's
        - velocity(0) = B

    .. plot::
        :context: close-figs

        >>> from pyyeti import ode
        >>> from pyyeti.ssmodel import SSModel
        >>> import numpy as np
        >>> import matplotlib.pyplot as plt
        >>> f = 5           # 5 hz oscillator
        >>> w = 2*np.pi*f
        >>> w2 = w*w
        >>> zeta = .05
        >>> h = .01
        >>> nt = 500
        >>> A = np.array([[0, 1], [-w2, -2*w*zeta]])
        >>> B = np.array([[0], [3]])
        >>> C = np.array([[8, -5]])
        >>> D = np.array([[0]])
        >>> F = np.zeros((1, nt), float)
        >>> ts = ode.SolveExp1(A, h)
        >>> sol = ts.tsolve(B @ F, B[:, 0])
        >>> y = C @ sol.d
        >>> fig = plt.figure('Example')
        >>> fig.clf()
        >>> ax = plt.plot(sol.t, y.T,
        ...               label='SolveExp1')
        >>> ssmodel = SSModel(A, B, C, D)
        >>> z = ssmodel.c2d(h=h, method='zoh')
        >>> x = np.zeros((A.shape[1], nt+1), float)
        >>> y2 = np.zeros((C.shape[0], nt), float)
        >>> x[:, 0:1] = B
        >>> for k in range(nt):
        ...     x[:, k+1] = z.A @ x[:, k] + z.B @ F[:, k]
        ...     y2[:, k]  = z.C @ x[:, k] + z.D @ F[:, k]
        >>> ax = plt.plot(sol.t, y2.T, label='discrete')
        >>> leg = plt.legend(loc='best')
        >>> np.allclose(y, y2)
        True

        Compare against scipy:

        >>> from scipy import signal
        >>> ss = ssmodel.getlti()
        >>> tout, yout = ss.impulse(T=sol.t)
        >>> ax = plt.plot(tout, yout, label='scipy')
        >>> leg = plt.legend(loc='best')
        >>> np.allclose(yout, y.ravel())
        True
    """

    def __init__(self, A, h, order=1):
        """
        Instantiates a :class:`SolveExp1` solver.

        Parameters
        ----------
        A : 2d ndarray
            The state-space matrix: ``yd - A y = f``
        h : scalar or None
            Time step or None; if None, the `E`, `P`, `Q` members will
            not be computed.
        order : scalar, optional
            Specify which solver to use:

            - 0 for the zero order hold (force stays constant across
              time step)
            - 1 for the 1st order hold (force can vary linearly across
              time step)

        Notes
        -----
        The class instance is populated with the following members:

        =======   =================================================
        Member    Description
        =======   =================================================
        A         the input `A`
        h         time step
        n         number of total DOF (``A.shape[0]``)
        order     order of solver (0 or 1; see above)
        E, P, Q   output of :func:`pyyeti.expmint.getEPQ`; they are
                  matrices used to solve the ODE
        pc        True if E, P, and Q member have been calculated;
                  False otherwise
        =======   =================================================

        The E, P, and Q entries are used to solve the ODE::

            for j in range(1, nt):
                d[:, j] = E*d[:, j-1] + P*F[:, j-1] + Q*F[:, j]
        """
        if h:
            E, P, Q = expmint.getEPQ(A, h, order)
            self.E = E
            self.P = P
            self.Q = Q
        self.A = A
        self.h = h
        self.order = order
        self.n = A.shape[0]

    def tsolve(self, force, d0=None):
        """
        Solve time-domain 1st order ODE equations.

        Parameters
        ----------
        force : 2d ndarray
            The force matrix; ndof x time
        d0 : 1d ndarray; optional
            Displacement initial conditions; if None, zero ic's are
            used.

        Returns
        -------
        A record (SimpleNamespace class) with the members:

        d : 2d ndarray
            Displacement; ndof x time
        v : 2d ndarray
            Velocity; ndof x time
        h : scalar
            Time-step
        t : 1d ndarray
            Time vector: np.arange(d.shape[1])*h
        """
        force = np.atleast_2d(force)
        if force.shape[0] != self.n:
            raise ValueError(
                f"Force matrix has {force.shape[0]} rows; {self.n} rows are expected"
            )
        nt = force.shape[1]
        d = np.zeros((self.n, nt))  # , float, order='F')
        if d0 is not None:
            d[:, 0] = d0
        else:
            d0 = np.zeros(self.n, float)
        if nt > 1:
            if not self.h:
                raise RuntimeError("instantiate the class with a valid time step.")
            # calc force term outside loop:
            if self.order == 1:
                PQF = self.P @ force[:, :-1] + self.Q @ force[:, 1:]
            else:
                PQF = self.P @ force[:, :-1]
            E = self.E
            for j in range(1, nt):
                d0 = d[:, j] = E @ d0 + PQF[:, j - 1]
            t = self.h * np.arange(nt)
        else:
            t = np.array([0.0])
        return SimpleNamespace(d=d, v=force + self.A @ d, h=self.h, t=t)

max_stars_repo_name: twmacro/pyye

max_issues_repo_name: twmacro/pyye

max_forks_repo_name: twmacro/pyye


text:
import numpy as np
from geometer import (
    Point,
    Segment,
    Rectangle,
    Simplex,
    Triangle,
    Cuboid,
    Line,
    RegularPolygon,
    Polygon,
    SegmentCollection,
    PolygonCollection,
    PointCollection,
    dist,
    rotation,
    translation,
)

class TestSegment:
    def test_contains(self):
        # both points finite
        p = Point(0, 0)
        q = Point(2, 1)
        s = Segment(p, q)
        assert s.contains(p)
        assert s.contains(q)
        assert s.contains(0.5 * (p + q))
        assert not s.contains(p - q)
        assert not s.contains(Point([2, 1, 0]))

        # first point at infinity
        p = Point([-2, -1, 0])
        q = Point(2, 1)
        s = Segment(p, q)
        assert s.contains(p)
        assert s.contains(q)
        assert s.contains(0.5 * q)
        assert not s.contains(2 * q)

        # second point at infinity
        p = Point(0, 0)
        q = Point([2, 1, 0])
        s = Segment(p, q)
        assert s.contains(p)
        assert s.contains(q)
        assert s.contains(Point(2, 1))
        assert not s.contains(Point(-2, -1))

        # both points at infinity
        p = Point([-2, 1, 0])
        q = Point([2, 1, 0])
        s = Segment(p, q)
        assert s.contains(p)
        assert s.contains(q)
        assert s.contains(0.5 * (p + q))
        assert not s.contains(p - q)
        assert not s.contains(Point(0, 0))

    def test_equal(self):
        p = Point(0, 0)
        q = Point(2, 1)
        s = Segment(p, q)

        assert s == Segment(q, p)
        assert s == s
        assert s == Segment([(0, 0), (2, 1)], homogenize=True)
        assert s != Segment([(0, 0), (1, 2)], homogenize=True)

    def test_intersect(self):
        a = Point(0, 0)
        b = Point(0, 2)
        c = Point(2, 2)
        d = Point(2, 0)
        s1 = Segment(a, c)
        s2 = Segment(b, d)

        assert s1.intersect(s2) == [Point(1, 1)]

    def test_midpoint(self):
        p = Point(0, 0)
        q = Point(2, 2)
        s = Segment(p, q)
        assert s.midpoint == Point(1, 1)

        p = Point(0, 0, 0)
        q = Point(0, 2, 0)
        s = Segment(p, q)
        assert s.midpoint == Point(0, 1, 0)

    def test_transformation(self):
        p = Point(0, 0)
        q = Point(2, 2)
        s = Segment(p, q)

        r = rotation(np.pi / 2)
        assert r * s == Segment(p, Point(-2, 2))
        assert r.apply(s)._line == r.apply(s._line)

    def test_getitem(self):
        p = Point(0, 0)
        q = Point(2, 2)
        s = Segment(p, q)

        assert s[0] == p
        assert s[1] == q

class TestPolygon:
    def test_equal(self):
        points = np.random.rand(50, 3)

        p1 = Polygon(points)
        p2 = Polygon(*[Point(p) for p in points])

        assert p1 == p2

    def test_intersect(self):
        a = Point(0, 0)
        b = Point(0, 2)
        c = Point(2, 2)
        d = Point(2, 0)
        r = Rectangle(a, b, c, d)
        s = Segment(a, c)
        assert r.intersect(s) == [a, c]

        foo = Rectangle(
            Point(148.06094049635456, 10.151151779987144, 60.522099063951394),
            Point(129.78569335065157, -42.129870038015355, 60.54878245579997),
            Point(85.91668756014471, -26.79517716452499, 60.41723371984577),
            Point(104.19193470584759, 25.485844653477507, 60.390550327997126),
        )
        bar = Segment(
            Point(-38.9592826559563, -6.703132040294841, 64.78693707404751),
            Point(133.01711836447913, -6.633886165038485, 54.310634812542006),
        )

        assert len(foo.intersect(bar)) == 0

    def test_edges(self):
        a = Point(0, 0)
        b = Point(0, 2)
        c = Point(2, 2)
        d = Point(2, 0)
        r = Rectangle(a, b, c, d)
        assert r.edges == [Segment(a, b), Segment(b, c), Segment(c, d), Segment(d, a)]

    def test_contains(self):
        a = Point(0, 0)
        b = Point(0, 2)
        c = Point(2, 2)
        d = Point(2, 0)
        r = Rectangle(a, b, c, d)

        assert r.contains(a)
        assert r.contains(b)
        assert r.contains(c)
        assert r.contains(d)
        assert r.contains(Point(1, 1))
        assert not r.contains(Point([1, 1, 0]))

        a = Point(0, 0, 1)
        b = Point(1, 3, 1)
        c = Point(2, 0, 1)
        d = Point(1, 1, 1)
        p = Polygon(a, b, c, d)

        assert p.contains(Point(0.5, 1, 1))
        assert not p.contains(Point(0.5, 1, 0))
        assert np.all(p.contains(PointCollection([Point(0.5, 1, 1), Point(1.5, 1, 1)])))
        assert np.all(p.contains(PointCollection([a, c, d])))

        a = Point([1, 1, 2, 0])
        b = Point([-1, 1, 2, 0])
        c = Point([-1, -1, 2, 0])
        d = Point([1, -1, 2, 0])
        p = Polygon(a, b, c, d)

        assert all(p.contains(PointCollection([a, b, c, d])))
        assert p.contains(Point([0, 0, 1, 0]))
        assert not p.contains(Point([1, 1, 1, 0]))

    def test_transformation(self):
        a = Point(0, 0)
        b = Point(0, 1)
        c = Point(2, 1)
        d = Point(2, 0)
        r = Rectangle(a, b, c, d)
        r2 = rotation(np.pi / 2) * r

        assert r.area == r2.area
        assert r2.contains(Point(-0.5, 1.5))

        l = Line(Point(0, 0, -10), Point(0, 0, 10))
        r = Rectangle(
            Point(-10, -10, 0), Point(10, -10, 0), Point(10, 10, 0), Point(-10, 10, 0)
        )
        t = rotation(np.pi / 6, Point(1, 0, 0))

        assert r.intersect(l) == [Point(0, 0, 0)]
        assert (t * r).intersect(l) == [Point(0, 0, 0)]

    def test_copy(self):
        a = Point(0, 0)
        b = Point(0, 2)
        c = Point(2, 2)
        d = Point(2, 0)

        r1 = Rectangle(a, b, c, d)
        p1 = RegularPolygon(a, 1, 6)

        r2 = r1.copy()
        p2 = p1.copy()

        assert r1 == r2
        assert r1 is not r2
        assert r1.vertices == r2.vertices
        assert p1 == p2
        assert p1 is not p2
        assert p1.vertices == p2.vertices

    def test_centroid(self):
        a = Point(0, 0, 1)
        b = Point(2, 0, 1)
        c = Point(2, 2, 1)
        d = Point(0, 2, 1)
        r = Rectangle(a, b, c, d)

        assert r.centroid == Point(1, 1, 1)

    def test_getitem(self):
        a = Point(0, 0, 1)
        b = Point(2, 0, 1)
        c = Point(2, 2, 1)
        d = Point(0, 2, 1)
        r = Rectangle(a, b, c, d)

        assert r[0] == a
        assert r[1] == b
        assert r[2] == c
        assert r[3] == d

class TestTriangle:
    def test_contains(self):
        a = Point(0, 0)
        b = Point(0, 2)
        c = Point(2, 1)
        t = Triangle(a, b, c)

        assert t.contains(Point(1, 1))
        assert all(t.contains(PointCollection([a, b, c])))
        assert not t.contains(Point(-1, 1))
        assert not t.contains(Point([1, 1, 0]))

    def test_area(self):
        a = Point(0, 0)
        b = Point(2, 0)
        c = Point(0, 2)
        t = Triangle(a, b, c)
        assert np.isclose(t.area, 2)

    def test_centroid(self):
        a = Point(0, 0, 1)
        b = Point(2, 0, 1)
        c = Point(2, 2, 1)
        t = Triangle(a, b, c)

        s1, s2, s3 = t.edges
        l1 = s1.midpoint.join(c)
        l2 = s2.midpoint.join(a)

        assert t.centroid == (a + b + c) / 3
        assert t.centroid == l1.meet(l2)

    def test_circumcenter(self):
        a = Point(0, 0, 1)
        b = Point(2, 0, 1)
        c = Point(2, 2, 1)
        t = Triangle(a, b, c)

        assert t.circumcenter == Point(1, 1, 1)

class TestRegularPolygon:
    def test_init(self):
        a = Point(0, 0, 0)
        p = RegularPolygon(a, 1, 6, axis=Point(0, 0, 1))

        d = p.edges[0].length

        assert len(p.vertices) == 6
        assert np.isclose(dist(a, p.vertices[0]), 1)
        assert all(np.isclose(p.edges[1:].length, d))
        assert np.allclose(p.angles, np.pi / 3)
        assert p.center == a

    def test_radius(self):
        p = RegularPolygon(Point(0, 0, 0), 1, 6, axis=Point(0, 0, 1))

        assert np.isclose(p.radius, 1)
        assert np.isclose(p.inradius, np.cos(np.pi / 6))

    def test_transform(self):
        p = RegularPolygon(Point(0, 0, 0), 1, 6, axis=Point(0, 0, 1))
        t = translation(1, 1, 0)

        assert t * p == RegularPolygon(Point(1, 1, 0), 1, 6, axis=Point(0, 0, 1))
        assert isinstance(t * p, RegularPolygon)

class TestSimplex:
    def test_volume(self):
        a = Point(0, 0, 0)
        b = Point(1, 0, 0)
        c = Point(0, 1, 0)
        d = Point(0, 0, 1)
        s = Simplex(a, b, c, d)

        assert np.isclose(s.volume, 1 / 6)

        triangle = Simplex(a, b, c)
        assert np.isclose(triangle.volume, 1 / 2)

    def test_transform(self):
        a = Point(0, 0, 0)
        b = Point(1, 0, 0)
        c = Point(0, 1, 0)
        d = Point(0, 0, 1)
        s = Simplex(a, b, c, d)
        x = Point(1, 1, 1)
        t = translation(x)

        assert t * s == Simplex(a + x, b + x, c + x, d + x)

class TestCuboid:
    def test_intersect(self):
        a = Point(0, 0, 0)
        b = Point(1, 0, 0)
        c = Point(0, 1, 0)
        d = Point(0, 0, 1)
        cube = Cuboid(a, b, c, d)
        l = Line(Point(2, 0.5, 0.5), Point(-1, 0.5, 0.5))
        assert cube.intersect(l) == [Point(0, 0.5, 0.5), Point(1, 0.5, 0.5)]

    def test_edges(self):
        a = Point(0, 0, 0)
        b = Point(1, 0, 0)
        c = Point(0, 1, 0)
        d = Point(0, 0, 1)
        cube = Cuboid(a, b, c, d)

        assert len(cube.edges) == 12

    def test_area(self):
        a = Point(0, 0, 0)
        b = Point(1, 0, 0)
        c = Point(0, 1, 0)
        d = Point(0, 0, 1)
        cube = Cuboid(a, b, c, d)
        assert len(cube.faces) == 6
        assert len(cube.vertices) == 8
        assert cube.area == 6

    def test_transform(self):
        a = Point(0, 0, 0)
        b = Point(1, 0, 0)
        c = Point(0, 1, 0)
        d = Point(0, 0, 1)
        cube = Cuboid(a, b, c, d)
        x = Point(1, 1, 1)
        t = translation(x)

        assert t * cube == Cuboid(a + x, b + x, c + x, d + x)
        assert isinstance(t * cube, Cuboid)

    def test_add(self):
        a = Point(0, 0, 0)
        b = Point(1, 0, 0)
        c = Point(0, 1, 0)
        d = Point(0, 0, 1)
        cube = Cuboid(a, b, c, d)
        p = Point(1, 2, 3)

        assert cube + p == Cuboid(a + p, b + p, c + p, d + p)
        assert cube - p == Cuboid(a - p, b - p, c - p, d - p)

    def test_getitem(self):
        a = Point(0, 0, 0)
        b = Point(1, 0, 0)
        c = Point(0, 1, 0)
        d = Point(0, 0, 1)
        cube = Cuboid(a, b, c, d)

        x, y, z = b - a, c - a, d - a
        yz = Rectangle(a, a + z, a + y + z, a + y)
        xz = Rectangle(a, a + x, a + x + z, a + z)
        xy = Rectangle(a, a + x, a + x + y, a + y)

        assert isinstance(cube[0], Polygon)
        assert cube[0] == yz
        assert cube[1] == xz
        assert cube[2] == xy
        assert cube[0, 0] == a

class TestSegmentCollection:
    def test_contains(self):
        p = PointCollection([(0, 0), (1, 0)], homogenize=True)
        q = PointCollection([(2, 1), (3, 1)], homogenize=True)
        s = SegmentCollection(p, q)

        assert all(s.contains(p))
        assert all(s.contains(q))
        assert all(s.contains(0.5 * (p + q)))
        assert not any(s.contains(p - q))
        assert not any(s.contains(Point([2, 1, 0])))

    def test_midpoint(self):
        p = PointCollection([(0, 0), (2, 1)], homogenize=True)
        q = PointCollection([(2, 2), (4, 1)], homogenize=True)
        s = SegmentCollection(p, q)
        assert s.midpoint == PointCollection([(1, 1), (3, 1)], homogenize=True)

        p = PointCollection([(0, 0, 0), (1, 1, 1)], homogenize=True)
        q = PointCollection([(0, 2, 0), (3, 3, 3)], homogenize=True)
        s = SegmentCollection(p, q)
        assert s.midpoint == PointCollection([(0, 1, 0), (2, 2, 2)], homogenize=True)

    def test_intersect(self):
        a = PointCollection([(0, 0), (0, 0)], homogenize=True)
        b = PointCollection([(2, 0), (4, 0)], homogenize=True)
        c = PointCollection([(2, 2), (4, 4)], homogenize=True)
        d = PointCollection([(0, 2), (0, 4)], homogenize=True)
        s1 = SegmentCollection(a, c)
        s2 = SegmentCollection(b, d)

        assert s1.intersect(s2) == PointCollection([(1, 1), (2, 2)], homogenize=True)
        assert s1.intersect(Line(Point(0, 0), Point(1, 1))).size == 0


text:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec  # for unequal plot boxes
import scipy.optimize

# define function to calculate reduced chi-squared
def RedChiSqr(func, x, y, dy, params):
    resids = y - func(x, *params)
    chisq = ((resids/dy)**2).sum()
    return chisq/float(x.size-params.size)

# define fitting function
def SineGaussDecay(t, A, B, C, tau, omega):
    y = A * (1.0 + B*np.cos(omega*t)) * np.exp(-0.5*t*t/(tau*tau)) + C
    return y

# read in spectrum from data file
t, decay, unc = np.loadtxt("OscData.txt", skiprows=4, unpack=True)

# initial values for fitting parameters (guesses)
A0 = 15.0
B0 = 0.6
C0 = 1.2*A0
tau0 = 16.0
omega0 = 2.0 * (2.0*np.pi/tau0)
#omega0 = 2.34

# fit data using SciPy's Levenberg-Marquart method
nlfit, nlpcov = scipy.optimize.curve_fit(SineGaussDecay,
             t, decay, p0=[A0, B0, C0, tau0, omega0], sigma=unc)

# calculate reduced chi-squared
rchi = RedChiSqr(SineGaussDecay, t, decay, unc, nlfit)

# create fitting function from fitted parameters
A, B, C, tau, omega = nlfit
t_fit = np.linspace(0.0, 1.02*t[-1], 512)
d_fit = SineGaussDecay(t_fit, A, B, C, tau, omega)

# Create figure window to plot data
fig = plt.figure(1, figsize=(8,8))      # extra length for residuals
gs = gridspec.GridSpec(2, 1, height_ratios=[6, 2])

# Top plot: data and fit
ax1 = fig.add_subplot(gs[0])
ax1.plot(t_fit, d_fit)
ax1.errorbar(t, decay, yerr=unc, fmt='or', ecolor='black', ms=4)
ax1.set_xlabel('time (ms)')
ax1.set_ylabel('decay (arb units)')
ax1.text(0.55, 0.8, 'A = {0:0.1f}\nB = {1:0.3f}\nC = {2:0.1f}'.format(A, B, C),
        transform = ax1.transAxes)
ax1.text(0.75, 0.8, '$\\tau$ = {0:0.1f}\n$\omega$ = {1:0.3f}\n$\chi^2$ = {2:0.3f}'.format(tau, omega, rchi),
        transform = ax1.transAxes)
ax1.set_title('$d(t) = A (1+B\,\cos\,\omega t) e^{-t^2/2\\tau^2} + C$')

# Bottom plot: residuals
resids = decay - SineGaussDecay(t, A, B, C, tau, omega)
ax2 = fig.add_subplot(gs[1])
ax2.axhline(color="gray")
ax2.errorbar(t, resids, yerr = unc, ecolor="black", fmt="ro", ms=4)
ax2.set_xlabel('time (ms)')
ax2.set_ylabel('residuals')
ax2.set_ylim(-5, 5)
yticks = (-5, 0, 5)
ax2.set_yticks(yticks)

plt.savefig('FitOscDecay.pdf')

plt.show()

max_stars_repo_name: lorenghoh/pyman

max_issues_repo_name: lorenghoh/pyman

max_forks_repo_name: lorenghoh/pyman