text:
"""
Try out this finite-horizon q learning.
"""

import gym
from gym.envs.mujoco import HopperEnv
import numpy as np

import rlkit.torch.pytorch_util as ptu
from rlkit.envs.mujoco.discrete_reacher import DiscreteReacherEnv
from rlkit.envs.mujoco.discrete_swimmer import DiscreteSwimmerEnv
from rlkit.envs.mujoco.reacher_env import ReacherEnv
from rlkit.envs.wrappers import DiscretizeEnv
from rlkit.finite_q_learning.finite_horizon_dqn import FiniteHorizonDQN
from rlkit.launchers.launcher_util import run_experiment
from rlkit.torch.networks import Mlp
import rlkit.misc.hyperparameter as hyp
from rlkit.envs.mujoco.hopper_env import HopperEnv

def experiment(variant):
    env = variant['env_class'](**variant['env_kwargs'])
    env = DiscretizeEnv(env, variant['num_bins'])
    # env = DiscreteReacherEnv(**variant['env_kwargs'])

    qf = Mlp(
        input_size=int(np.prod(env.observation_space.shape)),
        output_size=env.action_space.n,
        **variant['qf_kwargs']
    )
    algorithm = FiniteHorizonDQN(
        env,
        qf,
        **variant['algo_kwargs']
    )
    algorithm.to(ptu.device)
    algorithm.train()

if __name__ == "__main__":
    n_seeds = 1
    mode = 'local'
    exp_prefix = 'dev'

    n_seeds = 3
    mode = 'ec2'
    exp_prefix = 'fhql-vs-ddqn-hooper-H1000-1000-networks-2'

    # noinspection PyTypeChecker
    variant = dict(
        algo_kwargs=dict(
            num_epochs=200,
            num_steps_per_epoch=1000,
            num_steps_per_eval=1000,
            batch_size=128,
            max_path_length=1000,
            max_horizon=1000,
            discount=1.,
            random_action_prob=0.05,
        ),
        env_kwargs=dict(
        ),
        algorithm="Finite-Horizon-DQN",
        num_bins=5,
    )

    search_space = {
        # 'algo_kwargs.discount': [0.99],
        # 'algo_kwargs.random_action_prob': [0.05, 0.2],
        'qf_kwargs.hidden_sizes': [[32, 32]],
        # 'env_kwargs.frame_skip': [2, 5],
        # 'env_kwargs.num_bins': [5],
        'env_class': [HopperEnv],
        'num_bins': [3],
    }
    sweeper = hyp.DeterministicHyperparameterSweeper(
        search_space, default_parameters=variant,
    )
    for exp_id, variant in enumerate(sweeper.iterate_hyperparameters()):
        for _ in range(n_seeds):
            run_experiment(
                experiment,
                exp_prefix=exp_prefix,
                mode=mode,
                variant=variant,
                exp_id=exp_id,
            )


text:
import numpy as np
import streamlit as st

embedding = {}
glove_vocab = []
def get_glove_vocab(embedding_dim, vocab_size, glove_vocab_size=6):
    with open(f"glove.{glove_vocab_size}B.{embedding_dim}d.txt", encoding='utf-8') as f:
            for line in f:
                values = line.split()
                glove_vocab.append(values[0])
    vocab = " ".join(w for w in glove_vocab)
    st.write(vocab)

def load_embed(embedding_dim, vocab_size, glove_vocab_size=6):
    try:
        with open(f"glove.{glove_vocab_size}B.{embedding_dim}d.txt", encoding='utf-8') as f:
            for line in f:
                values = line.split()
                embedding[values[0]] = np.asarray(values[1:], dtype = 'float32')
    except IOError:
        pass
    embedding_matrix = np.zeros((len(embedding.keys()), embedding_dim))
    return

load_embed(100, 263)

ext: py

max_stars_repo_name: streamlit-badge-bot/Light

max_issues_repo_name: streamlit-badge-bot/Light

max_forks_repo_name: streamlit-badge-bot/Light

reason: import numpy


text:
import os
import re
import time
import cv2
import numpy as np
from evaluate import evaluate

MAX_FEATURES = 500
GOOD_MATCH_PERCENT = 0.15

def getImageList(dir):
    """
    Searches through a directory and fetches all files
    """
    fileList = os.listdir(dir)

    rgbList = []
    for file in fileList:
        rgbList.append(cv2.cvtColor(cv2.imread(
            dir + file), cv2.COLOR_BGR2RGB))

    return rgbList

def alignImageList(originalList):
    """
    Take all images and get aligned versions in reference to the first image
    """
    alignedList = []

    for image in originalList:
        imReg, h = alignImages(image, originalList[0])
        alignedList.append(imReg)

    return alignedList

def alignImages(im1, im2):
    """
    Align Images based on features
    https://www.learnopencv.com/image-alignment-feature-based-using-opencv-c-python/
    """

    # Convert images to grayscale
    im1Gray = cv2.cvtColor(im1, cv2.COLOR_BGR2GRAY)
    im2Gray = cv2.cvtColor(im2, cv2.COLOR_BGR2GRAY)

    # Detect ORB features and compute descriptors.
    orb = cv2.ORB_create(MAX_FEATURES)
    keypoints1, descriptors1 = orb.detectAndCompute(im1Gray, None)
    keypoints2, descriptors2 = orb.detectAndCompute(im2Gray, None)

    # Match features.
    matcher = cv2.DescriptorMatcher_create(
        cv2.DESCRIPTOR_MATCHER_BRUTEFORCE_HAMMING)
    matches = matcher.match(descriptors1, descriptors2, None)

    # Sort matches by score
    matches.sort(key=lambda x: x.distance, reverse=False)

    # Remove not so good matches
    numGoodMatches = int(len(matches) * GOOD_MATCH_PERCENT)
    matches = matches[:numGoodMatches]

    # Draw top matches
    imMatches = cv2.drawMatches(
        im1, keypoints1, im2, keypoints2, matches, None)
    # cv2.imwrite("matches.jpg", imMatches)

    # Extract location of good matches
    points1 = np.zeros((len(matches), 2), dtype=np.float32)
    points2 = np.zeros((len(matches), 2), dtype=np.float32)

    for i, match in enumerate(matches):
        points1[i, :] = keypoints1[match.queryIdx].pt
        points2[i, :] = keypoints2[match.trainIdx].pt

    # Find homography
    h, mask = cv2.findHomography(points1, points2, cv2.RANSAC)

    # Use homography
    height, width, channels = im2.shape
    im1Reg = cv2.warpPerspective(im1, h, (width, height))

    return im1Reg, h

def blend(rgbList):
    """
    Median blending function
    """
    h = rgbList[0].shape[0]
    w = rgbList[0].shape[1]
    blended = np.copy(rgbList[0])

    # Initialize lists
    r = []
    g = []
    b = []

    start = time.time()
    for y in range(0, h):
        for x in range(0, w):
            # Loop and append to coresponding list
            for rgbImg in rgbList:
                # Don't include pure black pixels since they are usually from cropped/transformed images
                if all(pixelVal != 0 for pixelVal in rgbImg[y][x]):
                    r.append(rgbImg[y][x][0])
                    g.append(rgbImg[y][x][1])
                    b.append(rgbImg[y][x][2])

            # If there was a value that was appended, else don't change orignal
            if not (len(r) == 0 and len(g) == 0 and len(b) == 0):
                # Set current pixel to the median
                blended[y][x] = [np.median(r), np.median(g), np.median(b)]
                r.clear()
                b.clear()
                g.clear()

    end = time.time()

    # Convert back to BGR
    blendedBGR = cv2.cvtColor(blended, cv2.COLOR_RGB2BGR)

    # Write image
    return blendedBGR, end-start

if __name__ == '__main__':
    # Change to desired directory
    dir = "./blur/"

    # Retrieve list
    print("Retrieving list from %s" % dir)
    originalList = getImageList(dir)

    # Align images
    print("Aligning images...")
    rgbList = alignImageList(originalList)

    # Blend images
    print("Blending images...")
    print("(Warning: algorithm is VERY slow, please wait)")
    blendedBGR, timeElap = blend(rgbList)
    print("Process completed in %d seconds." % timeElap)

    # Output file
    filename = re.sub('[./]', '', dir)
    cv2.imwrite("%s_blended.jpg" % filename, blendedBGR)
    print("File generated.")

    # Check if a ground truth directory exists
    grndTDir = dir.replace('./', '').replace('/', '_grndT')
    if os.path.isdir(grndTDir):
        grndT = cv2.imread(grndTDir + "/grndtruth.jpg")
        evaluate(blendedBGR, grndT)  # Run test (in evaluate.py)
    else:
        print("> There was no ground truth folder found")

    # Display result and a input file
    cv2.imshow("Reference Image", cv2.cvtColor(rgbList[0], cv2.COLOR_RGB2BGR))
    cv2.imshow("Result", blendedBGR)
    print("Images displayed, press any key to close...")
    cv2.waitKey(0)
    cv2.destroyAllWindows()

max_stars_repo_name: tjagda/thanos-snap

max_issues_repo_name: tjagda/thanos-snap

max_forks_repo_name: tjagda/thanos-snap


text:
"""Simulation coupled rossler system with different network topology"""
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import random
import network_topology
from mpl_toolkits.mplot3d import Axes3D

def runge_kutta(f, tspan, Y0, Nt):
    '''f deviation function;
    tspan 2*1 vector tspan(1)start time tspan(2) end time;
    Y0 is a columb vector, Y0(i) is the initial value of i-th variable;
    Nt is the number of t'''
    Nvar = Y0.shape[1]  #number of variables
    dt = (tspan[1]-tspan[0])/(Nt-1)
    Y = np.zeros((Nvar, Nt))   #initial value of Y
    Y[:,0] = Y0
    t = list(np.linspace(tspan[0], tspan[1], Nt))
    for i in range(0, Nt-1):
        K1 = f(Y[:,i], t[i])
        K2 = f(Y[:,i]+(K1*dt/2).flatten('F'), t[i]+dt/2.0)
        K3 = f(Y[:,i]+(K2*dt/2).flatten('F'), t[i]+dt/2.0)
        K4 = f(Y[:,i]+(K3*dt/2).flatten('F'), t[i]+dt)
        Y[:,i+1] = Y[:,i] + dt/6*(K1+2*K2+2*K3+K4).flatten('F')
    return Y

def function_g(n, A, sub_n, H, X):
    """coupling function g(x,y)
    n--> node n
    sub_n --> dimension x11
    A--> adjacent matrix of node n x1"""
    arrow_n = A[n]
    neighbours = np.nonzero(arrow_n)[0] # neighbours of node n
    degree_n = np.sum(A,axis=1)[n]
    g1 = sum([np.dot(H[sub_n], X[i]) for i in list(neighbours)])
    g = g1 - degree_n*(np.dot(H[sub_n], X[n]))
    return g

def f_func(Y, t, A):
    """A -- adjacent"""
    YF= Y.flatten('F')
    step = 3
    X = [YF[i:i + step] for i in range(0, len(YF), step)]  # variable
    a = 0.2
    b = 0.2
    c = 6.0
    coup =0.03
    Nvar1 = len(list(YF)) #number of variables
    Y1 = np.zeros((Nvar1, 1))
    H = np.array([[1, 0, 0],
                  [0, 0, 0],
                  [0, 0, 0]])
    for i in range(Nvar1):
        a1 = i//step  # X1
        b1 = i%step   # X11
        coupling_g = function_g(a1, A, b1, H, X)
        if b1 == 0:
            Y1[i, 0] = -X[a1][1] - X[a1][2] + coup*(coupling_g)
        elif b1 == 1:
            Y1[i, 0] = X[a1][0] + a*X[a1][1] + coup*(coupling_g)
        elif b1 == 2:
            Y1[i, 0] = b + X[a1][2]*(X[a1][0]-c) + coup*(coupling_g)
    return Y1

def keplerRK4(A, AA):
    Dimension = 3  # 3 dimension
    N_node = A.shape[1]  # number of nodes
    Y0 = np.zeros((1, Dimension*N_node))
    for i in range(Dimension*N_node):
        random.seed(i)
        Y0[0, i] = i/100.0
    tspan = [0, 100]
    Nt = 1000
    f = lambda Y,t:f_func(Y,t,A)
    Y = runge_kutta(f, tspan, Y0, Nt)
    f1 = lambda Y1, t: f_func(Y1, t, AA)
    Y1 = runge_kutta(f1, tspan, Y0, Nt)

    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot(Y1[0], Y1[1], Y1[2])
    ax.plot(Y1[3], Y1[4], Y1[5])
    plt.show()

    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot(Y[0], Y[1], Y[2])
    ax.plot(Y[3], Y[4], Y[5])
    plt.show()

    plt.plot(Y[0], Y[1], label='node1')
    plt.plot(Y[3], Y[4], label='node2')
    plt.legend()
    plt.xlabel(r'$x^{(1)}$')
    plt.ylabel(r'$x^{(2)}$')
    plt.savefig('attractor2.pdf')
    plt.show()

    fig3 = plt.figure(figsize=(20, 10))
    #fig3 = plt.figure()
    plt.subplot(2, 3, 1)
    c1 = [Y[3][i] - Y[0][i] for i in range(0, len(Y[5]))]
    c2 = [Y[4][i] - Y[1][i] for i in range(0, len(Y[5]))]
    c3 = [Y[5][i] - Y[2][i] for i in range(0, len(Y[5]))]
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c1, label=r'$x_{2}^{(1)}-x_{1}^{(1)}$')
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c2, label=r'$x_{2}^{(2)}-x_{1}^{(2)}$')
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c3, label=r'$x_{2}^{(3)}-x_{1}^{(3)}$')
    plt.xlabel('t')
    plt.legend()

    plt.subplot(2, 3, 2)
    c1 = [Y[9][i] - Y[6][i] for i in range(0, len(Y[5]))]
    c2 = [Y[10][i] - Y[7][i] for i in range(0, len(Y[5]))]
    c3 = [Y[11][i] - Y[8][i] for i in range(0, len(Y[5]))]
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c1, label=r'$x_{4}^{(1)}-x_{3}^{(1)}$')
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c2, label=r'$x_{4}^{(2)}-x_{3}^{(2)}$')
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c3, label=r'$x_{4}^{(3)}-x_{3}^{(3)}$')
    plt.xlabel('t')
    plt.legend()

    plt.subplot(2, 3, 3)
    c1 = [Y[15][i] - Y[12][i] for i in range(0, len(Y[5]))]
    c2 = [Y[16][i] - Y[13][i] for i in range(0, len(Y[5]))]
    c3 = [Y[17][i] - Y[14][i] for i in range(0, len(Y[5]))]
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c1, label=r'$x_{6}^{(1)}-x_{5}^{(1)}$')
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c2, label=r'$x_{6}^{(2)}-x_{5}^{(2)}$')
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c3, label=r'$x_{6}^{(3)}-x_{5}^{(3)}$')
    plt.xlabel('t')
    plt.legend()

    plt.subplot(2, 3, 4)
    c1 = [Y1[3][i] - Y1[0][i] for i in range(0, len(Y1[5]))]
    c2 = [Y1[4][i] - Y1[1][i] for i in range(0, len(Y1[5]))]
    c3 = [Y1[5][i] - Y1[2][i] for i in range(0, len(Y1[5]))]
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c1, label=r'$x_{2}^{(1)}-x_{1}^{(1)}$')
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c2, label=r'$x_{2}^{(2)}-x_{1}^{(2)}$')
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c3, label=r'$x_{2}^{(3)}-x_{1}^{(3)}$')
    plt.xlabel('t')
    plt.legend()

    plt.subplot(2, 3, 5)
    c1 = [Y1[9][i] - Y1[6][i] for i in range(0, len(Y1[5]))]
    c2 = [Y1[10][i] - Y1[7][i] for i in range(0, len(Y1[5]))]
    c3 = [Y1[11][i] - Y1[8][i] for i in range(0, len(Y1[5]))]
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c1, label=r'$x_{4}^{(1)}-x_{3}^{(1)}$')
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c2, label=r'$x_{4}^{(2)}-x_{3}^{(2)}$')
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c3, label=r'$x_{4}^{(3)}-x_{3}^{(3)}$')
    plt.xlabel('t')
    plt.legend()

    plt.subplot(2, 3, 6)
    c1 = [Y1[15][i] - Y1[12][i] for i in range(0, len(Y1[5]))]
    c2 = [Y1[16][i] - Y1[13][i] for i in range(0, len(Y1[5]))]
    c3 = [Y1[17][i] - Y1[14][i] for i in range(0, len(Y1[5]))]
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c1, label=r'$x_{6}^{(1)}-x_{5}^{(1)}$')
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c2, label=r'$x_{6}^{(2)}-x_{5}^{(2)}$')
    plt.plot(list(np.linspace(tspan[0], tspan[1], Nt)), c3, label=r'$x_{6}^{(3)}-x_{5}^{(3)}$')
    plt.xlabel('t')
    plt.legend()

    plt.savefig('node12desynch.pdf')
    plt.show()

A = network_topology.A # dessyn
AA = network_topology.AA # syn
keplerRK4(A, AA)

max_stars_repo_name: Martinezou/synchronization

max_issues_repo_name: Martinezou/synchronization

max_forks_repo_name: Martinezou/synchronization

reason: import numpy,import networkx


text:
#!/usr/bin/env python
"""
Simple Python interface to pupil_function.c

Hazen 10/17
"""
import ctypes
import numpy
from numpy.ctypeslib import ndpointer
import os

import storm_analysis.sa_library.loadclib as loadclib

import storm_analysis.simulator.pupil_math as pupilMath

#import storm_analysis.sa_library.recenter_psf as recenterPSF

pupil_fn = loadclib.loadCLibrary("storm_analysis.sa_library", "pupil_function")

pupil_fn.pfnCleanup.argtypes = [ctypes.c_void_p]

pupil_fn.pfnGetPSF.argtypes = [ctypes.c_void_p,
                               ndpointer(dtype = numpy.float64),
                               ndpointer(dtype = numpy.float64)]

pupil_fn.pfnGetPSFdx.argtypes = [ctypes.c_void_p,
                                 ndpointer(dtype = numpy.float64),
                                 ndpointer(dtype = numpy.float64)]

pupil_fn.pfnGetPSFdy.argtypes = [ctypes.c_void_p,
                                 ndpointer(dtype = numpy.float64),
                                 ndpointer(dtype = numpy.float64)]

pupil_fn.pfnGetPSFdz.argtypes = [ctypes.c_void_p,
                                 ndpointer(dtype = numpy.float64),
                                 ndpointer(dtype = numpy.float64)]

pupil_fn.pfnInitialize.argtypes = [ndpointer(dtype = numpy.float64),
                                   ndpointer(dtype = numpy.float64),
                                   ndpointer(dtype = numpy.float64),
                                   ctypes.c_int]
pupil_fn.pfnInitialize.restype = ctypes.c_void_p

pupil_fn.pfnSetPF.argtypes = [ctypes.c_void_p,
                              ndpointer(dtype = numpy.float64),
                              ndpointer(dtype = numpy.float64)]

pupil_fn.pfnTranslate.argtypes = [ctypes.c_void_p,
                                  ctypes.c_double,
                                  ctypes.c_double,
                                  ctypes.c_double]

class PupilFunction(object):

    def __init__(self, geometry = None, **kwds):
        """
        geometry is a simulation.pupil_math.Geometry object.

        Note: All of the getXY functions return complex numbers. In order
              to convert the output of getPSF() to a magnitude you need
              to do mag = psf*numpy.conj(psf).
        """
        super(PupilFunction, self).__init__(**kwds)
        assert isinstance(geometry, pupilMath.Geometry)

        self.size = geometry.size
        # Size must be an even number.
        assert ((self.size%2)==0)

        # geometry.kz will be a complex number, but the magnitude of the
        # imaginary component is zero so we just ignore it.
        self.pfn = pupil_fn.pfnInitialize(numpy.ascontiguousarray(geometry.kx, dtype = numpy.float64),
                                          numpy.ascontiguousarray(geometry.ky, dtype = numpy.float64),
                                          numpy.ascontiguousarray(numpy.real(geometry.kz), dtype = numpy.float64),
                                          geometry.size)

    def cleanup(self):
        pupil_fn.pfnCleanup(self.pfn)
        self.pfn = None

    def getCPointer(self):
        return self.pfn

    def getPSF(self):
        return self.getXX(pupil_fn.pfnGetPSF)

    def getPSFdx(self):
        return self.getXX(pupil_fn.pfnGetPSFdx)

    def getPSFdy(self):
        return self.getXX(pupil_fn.pfnGetPSFdy)

    def getPSFdz(self):
        return self.getXX(pupil_fn.pfnGetPSFdz)

    def getXX(self, fn):
        r = numpy.zeros((self.size, self.size), dtype = numpy.float64)
        c = numpy.zeros((self.size, self.size), dtype = numpy.float64)
        fn(self.pfn, r, c)
        return r + 1j*c

    def setPF(self, pf):
        pupil_fn.pfnSetPF(self.pfn,
                          numpy.ascontiguousarray(numpy.real(pf), dtype = numpy.float64),
                          numpy.ascontiguousarray(numpy.imag(pf), dtype = numpy.float64))

    def translate(self, dx, dy, dz):
        pupil_fn.pfnTranslate(self.pfn, dx, dy, dz)

max_stars_repo_name: oxfordni/storm-analysis

max_issues_repo_name: oxfordni/storm-analysis

max_forks_repo_name: oxfordni/storm-analysis

reason: import numpy,from numpy


text:
"""

"""

import unittest

import numpy as np

from probnum.filtsmooth.statespace.discrete import discretegaussianmodel

TEST_NDIM = 4

class TestDiscreteGaussianModel(unittest.TestCase):
    """
    """

    def setUp(self):
        """
        """
        dynamat = np.random.rand(TEST_NDIM, TEST_NDIM)
        diffmat = dynamat @ dynamat.T + np.eye(TEST_NDIM)
        self.nl = discretegaussianmodel.DiscreteGaussianModel(
            lambda t, x: dynamat @ x, lambda t: diffmat
        )

    def test_dynamics(self):
        """
        """
        some_input = np.random.rand(TEST_NDIM)
        val = self.nl.dynamics(0.0, some_input)
        self.assertEqual(val.ndim, 1)
        self.assertEqual(val.shape[0], TEST_NDIM)

    def test_diffusionmatrix(self):
        """
        """
        val = self.nl.diffusionmatrix(0.0)
        self.assertEqual(val.ndim, 2)
        self.assertEqual(val.shape[0], TEST_NDIM)
        self.assertEqual(val.shape[1], TEST_NDIM)

    def test_jacobian(self):
        """
        """
        some_input = np.random.rand(TEST_NDIM)
        with self.assertRaises(NotImplementedError):
            self.nl.jacobian(0.0, some_input)

    def test_sample(self):
        """
        """
        some_input = np.random.rand(TEST_NDIM)
        samp = self.nl.sample(0.0, some_input)
        self.assertEqual(samp.ndim, 1)
        self.assertEqual(samp.shape[0], TEST_NDIM)

    def test_ndim(self):
        """
        """
        self.assertEqual(self.nl.ndim, TEST_NDIM)

    def test_pdf(self):
        """
        """
        some_state = np.random.rand(TEST_NDIM)
        evl = self.nl.pdf(some_state, 0.0, some_state)
        self.assertEqual(np.isscalar(evl), True)

class TestLinear(unittest.TestCase):
    """
    """

    def setUp(self):
        """
        """
        dynamat = np.random.rand(TEST_NDIM, TEST_NDIM)
        diffmat = dynamat @ dynamat.T + np.eye(TEST_NDIM)
        self.lin = discretegaussianmodel.DiscreteGaussianLinearModel(
            lambda t: dynamat, lambda t: np.random.rand(TEST_NDIM), lambda t: diffmat
        )

    def test_dynamics(self):
        """
        """
        some_input = np.random.rand(TEST_NDIM)
        val = self.lin.dynamics(0.0, some_input)
        self.assertEqual(val.ndim, 1)
        self.assertEqual(val.shape[0], TEST_NDIM)

    def test_diffusionmatrix(self):
        """
        """
        val = self.lin.diffusionmatrix(0.0)
        self.assertEqual(val.ndim, 2)
        self.assertEqual(val.shape[0], TEST_NDIM)
        self.assertEqual(val.shape[1], TEST_NDIM)

    def test_jacobian(self):
        """
        """
        some_input = np.random.rand(TEST_NDIM)
        jac = self.lin.jacobian(0.0, some_input)
        self.assertEqual(jac.ndim, 2)
        self.assertEqual(jac.shape[0], TEST_NDIM)
        self.assertEqual(jac.shape[1], TEST_NDIM)

    def test_dynamicsmatrix(self):
        """
        """
        dyna = self.lin.dynamicsmatrix(0.0)
        self.assertEqual(dyna.ndim, 2)
        self.assertEqual(dyna.shape[0], TEST_NDIM)
        self.assertEqual(dyna.shape[1], TEST_NDIM)

    def test_force(self):
        """
        """
        force = self.lin.force(0.0)
        self.assertEqual(force.ndim, 1)
        self.assertEqual(force.shape[0], TEST_NDIM)

    def test_sample(self):
        """
        """
        some_input = np.random.rand(TEST_NDIM)
        samp = self.lin.sample(0.0, some_input)
        self.assertEqual(samp.ndim, 1)
        self.assertEqual(samp.shape[0], TEST_NDIM)

    def test_ndim(self):
        """
        """
        self.assertEqual(self.lin.ndim, TEST_NDIM)

    def test_pdf(self):
        """
        """
        some_state = np.random.rand(TEST_NDIM)
        evl = self.lin.pdf(some_state, 0.0, some_state)
        self.assertEqual(np.isscalar(evl), True)

class TestLTI(unittest.TestCase):
    """
    """

    def setUp(self):
        """
        """
        dynamat = np.random.rand(TEST_NDIM, TEST_NDIM)
        diffmat = dynamat @ dynamat.T + np.eye(TEST_NDIM)
        self.lti = discretegaussianmodel.DiscreteGaussianLTIModel(
            dynamat, dynamat[0], diffmat
        )

    def test_dynamics(self):
        """
        """
        some_input = np.random.rand(TEST_NDIM)
        val = self.lti.dynamics(0.0, some_input)
        self.assertEqual(val.ndim, 1)
        self.assertEqual(val.shape[0], TEST_NDIM)

    def test_dynamicsmatrix(self):
        """
        """
        some_input = np.random.rand(TEST_NDIM)
        dyna = self.lti.dynamicsmatrix(0.0)
        self.assertEqual(dyna.ndim, 2)
        self.assertEqual(dyna.shape[0], TEST_NDIM)
        self.assertEqual(dyna.shape[1], TEST_NDIM)

    def test_diffusionmatrix(self):
        """
        """
        val = self.lti.diffusionmatrix(0.0)
        self.assertEqual(val.ndim, 2)
        self.assertEqual(val.shape[0], TEST_NDIM)
        self.assertEqual(val.shape[1], TEST_NDIM)

    def test_jacobian(self):
        """
        """
        some_input = np.random.rand(TEST_NDIM)
        jac = self.lti.jacobian(0.0, some_input)
        self.assertEqual(jac.ndim, 2)
        self.assertEqual(jac.shape[0], TEST_NDIM)
        self.assertEqual(jac.shape[1], TEST_NDIM)

    def test_force(self):
        """
        """
        force = self.lti.force(0.0)
        self.assertEqual(force.ndim, 1)
        self.assertEqual(force.shape[0], TEST_NDIM)

    def test_sample(self):
        """
        """
        some_input = np.random.rand(TEST_NDIM)
        samp = self.lti.sample(0.0, some_input)
        self.assertEqual(samp.ndim, 1)
        self.assertEqual(samp.shape[0], TEST_NDIM)

    def test_ndim(self):
        """
        """
        self.assertEqual(self.lti.ndim, TEST_NDIM)

    def test_pdf(self):
        """
        """
        some_state = np.random.rand(TEST_NDIM)
        evl = self.lti.pdf(some_state, 0.0, some_state)
        self.assertEqual(np.isscalar(evl), True)

max_stars_repo_name: coldfix/probnum

max_issues_repo_name: coldfix/probnum

max_forks_repo_name: coldfix/probnum


text:
from PySide6.QtWidgets import QMainWindow, QMessageBox
from PySide6.QtCore import Signal, QObject, QThread

from ui_mainwindowtrabajofinal import Ui_MainWindowTrabajoFinal

from abbengine import ABBEngine

import sim
import numpy as np
import time
import csv

class MainWindowTrabajoFinal(QMainWindow):
    # Eventos
    closed = Signal()

    # Comunicacion con el hilo
    sceneStart = Signal()
    sceneStop = Signal()

    def __init__(self, parent=None):
        super(MainWindowTrabajoFinal, self).__init__(parent)

        self.ui = Ui_MainWindowTrabajoFinal()
        self.ui.setupUi(self)

        # Controlador de escena
        self.sceneThread = QThread(self)
        self.scene = TrabajoFinal()

        self.scene.moveToThread(self.sceneThread)

        self.sceneThread.start()

        # Eventos de botones
        self.ui.pushButtonIniciar.clicked.connect(self.pushButtonIniciar_onCliked)
        self.ui.pushButtonDetener.clicked.connect(self.pushButtonDetener_onCliked)

        # Conexiones
        self.sceneStart.connect(self.scene.start)
        self.sceneStop.connect(self.scene.stop)

        self.scene.sceneStarted.connect(self.sceneStarted)
        self.scene.sceneError.connect(self.sceneError)
        self.scene.sceneStoped.connect(self.sceneStop)

        self.scene.sceneCarObsSensorChanged.connect(self.sceneCarObsSensorChanged)
        self.scene.sceneCarSensorChanged.connect(self.sceneCarSensorChanged)

    def closeEvent(self, event):
        self.sceneStop.emit()

        self.sceneThread.exit()
        self.sceneThread.wait(1000)

        self.closed.emit()

    def pushButtonIniciar_onCliked(self):
        self.sceneStart.emit()

    def pushButtonDetener_onCliked(self):
        self.sceneStop.emit()

    def sceneStarted(self):
        QMessageBox(QMessageBox.Information, "Escena", "La escena se inicio correctamente",
                    QMessageBox.Ok, self).exec()

    def sceneError(self, error):
        if error == 0:
            QMessageBox(QMessageBox.Critical, "Escena", "Error de conexión",
                        QMessageBox.Ok, self).exec()

        elif error == 1:
            QMessageBox(QMessageBox.Critical, "Escena", "Error de configuracion",
                        QMessageBox.Ok, self).exec()

    def sceneStoped(self):
        QMessageBox(QMessageBox.Warning, "Escena", "La escena se detuvo correctamente",
                    QMessageBox.Ok, self).exec()

    def sceneCarSensorChanged(self, values):
        self.ui.labelCarSensor0.setText(str("%.4f" % values[0]))
        self.ui.labelCarSensor1.setText(str("%.4f" % values[1]))
        self.ui.labelCarSensor2.setText(str("%.4f" % values[2]))
        self.ui.labelCarSensor3.setText(str("%.4f" % values[3]))
        self.ui.labelCarSensor4.setText(str("%.4f" % values[4]))
        self.ui.labelCarSensor5.setText(str("%.4f" % values[5]))

    def sceneCarObsSensorChanged(self, values):
        self.ui.labelCarObsSensor1.setText(str("%.4f" % values[0]))
        self.ui.labelCarObsSensor8.setText(str("%.4f" % values[7]))
        self.ui.labelCarObsSensor9.setText(str("%.4f" % values[8]))
        self.ui.labelCarObsSensor10.setText(str("%.4f" % values[9]))
        self.ui.labelCarObsSensor11.setText(str("%.4f" % values[10]))
        self.ui.labelCarObsSensor12.setText(str("%.4f" % values[11]))
        self.ui.labelCarObsSensor13.setText(str("%.4f" % values[12]))
        self.ui.labelCarObsSensor14.setText(str("%.4f" % values[13]))
        self.ui.labelCarObsSensor15.setText(str("%.4f" % values[14]))
        self.ui.labelCarObsSensor16.setText(str("%.4f" % values[15]))

class TrabajoFinal(QObject):
    # Señales
    sceneStarted = Signal()
    sceneError = Signal(int)
    sceneStoped = Signal()

    sceneCarObsSensorChanged = Signal(list)
    sceneCarSensorChanged = Signal(list)

    def __init__(self, parent=None):
        super(TrabajoFinal, self).__init__(parent)

        # Bandera que detiene la ejecucion
        self.programInterrupt = False

        # Posiciones de interes
        self.abb1_p_home = [374.0 + 67.0, 0.0, 630.0, np.pi, np.pi / 2.0, 0.0]
        self.abb2_p_home = [200.0, 0.0, 630.0, np.pi, np.pi / 2.0, 0.0]

        # Valores iniciales
        self.coppeliaIdClient = -1

        # Brazo 1
        # Valores actuales del robot
        self.abb1_q1 = 0.0
        self.abb1_q2 = 0.0
        self.abb1_q3 = 0.0
        self.abb1_q4 = 0.0
        self.abb1_q5 = 0.0
        self.abb1_q6 = 0.0

        self.abb1_x = self.abb1_p_home[0]
        self.abb1_y = self.abb1_p_home[1]
        self.abb1_z = self.abb1_p_home[2]

        self.abb1_a = self.abb1_p_home[3]
        self.abb1_b = self.abb1_p_home[4]
        self.abb1_c = self.abb1_p_home[5]

        # Juntas
        self.coppeliaABB1JointHandle = [None, None, None, None, None, None]

        # Vector de trayectorias
        self.abb1Trajectory = []

        # Motor de cinematica
        self.abb1Engine = ABBEngine()

        # Herramienta
        self.abb1Engine.tool = np.array(
            [
                [1.0, 0.0, 0.0, 0.0],
                [0.0, 1.0, 0.0, 0.0],
                [0.0, 0.0, 1.0, 67.0],
                [0.0, 0.0, 0.0, 1.0]
            ]
        )

        # Brazo 2
        # Valores actuales del robot
        self.abb2_q1 = 0.0
        self.abb2_q2 = 0.0
        self.abb2_q3 = 0.0
        self.abb2_q4 = 0.0
        self.abb2_q5 = 0.0
        self.abb2_q6 = 0.0

        self.abb2_x = self.abb1_p_home[0]
        self.abb2_y = self.abb1_p_home[1]
        self.abb2_z = self.abb1_p_home[2]

        self.abb2_a = self.abb1_p_home[3]
        self.abb2_b = self.abb1_p_home[4]
        self.abb2_c = self.abb1_p_home[5]

        # Juntas
        self.coppeliaABB2JointHandle = [None, None, None, None, None, None]

        # Vector de trayectorias
        self.abb2Trajectory = []

        # Motor de cinematica
        self.abb2Engine = ABBEngine()

        # Herramienta
        self.abb2Engine.tool = np.array(
            [
                [1.0, 0.0, 0.0, 0.0],
                [0.0, 1.0, 0.0, 0.0],
                [0.0, 0.0, 1.0, 67.0],
                [0.0, 0.0, 0.0, 1.0]
            ]
        )

        # Autito
        # Parametros del autito
        self.carWheelsAxisDistanceL = 1
        self.carWheelsRadius = 1

        # Objetos de la escena de Coppelia
        # Ruedas
        self.carLeftJoint = None
        self.carRigthJoint = None

        # Sensores
        self.carSensor = [None, None, None, None, None, None, None, None]
        self.carObsSensor = []

        for i in range(0, 16):
            self.carObsSensor.append(None)

    def start(self):
        if self.coppeliaConnect("127.0.0.1", 19999):
            if self.coppeliaBegin():
                self.sceneStarted.emit()

                self.exec()

            else:
                self.sceneError.emit(1)

        else:
            self.sceneError.emit(0)

    def stop(self):
        self.coppeliaDisconnect()

        self.sceneStoped.emit()

    def exec(self):
        # Establecemos las trayectorias
        # Brazo 1
        if not self.abb1LoadTrajectory():
            return

        # Brazo 2
        if not self.abb2LoadTrajectory():
            return

        # Movimientos al home
        self.abb1ExecLinTrajectory(self.abb1_p_home[0], self.abb1_p_home[1], self.abb1_p_home[2],
                                   self.abb1_p_home[3], self.abb1_p_home[4], self.abb1_p_home[5])

        self.abb2ExecLinTrajectory(self.abb2_p_home[0], self.abb2_p_home[1], self.abb2_p_home[2],
                                   self.abb2_p_home[3], self.abb2_p_home[4], self.abb2_p_home[5])

        # Movimientos del brazo 1
        for i in range(0, len(self.abb1Trajectory)):
            self.abb1ExecLinTrajectory(self.abb1Trajectory[i][0], self.abb1Trajectory[i][1],
                                       self.abb1Trajectory[i][2], self.abb1Trajectory[i][3],
                                       self.abb1Trajectory[i][4], self.abb1Trajectory[i][5])

            if (i == 3) or (i == 10) or (i == 18):
                self.coppeliaABB1SetVacuumGripper(1)

            if (i == 6) or (i == 14) or (i == 22):
                self.coppeliaABB1SetVacuumGripper(0)

        # Autito
        car_complete = False

        # Seguidor de linea
        # Ponderacion
        lin_pond = [2.0, 1.0, 0.5, 0.3, 0.3, 0.5, 1.0, 2.0]

        lin_kp = 2.5
        lin_kd = 2.5

        lin_error_prev = 0.0

        # Esquive de obstaculos
        box_kp = 2.0

        # Maquina de estados para la solución del problema
        car_state = 0

        while not car_complete:
            # Leemos los valores del ultrasonico
            car_obs_sensor_value = self.carReadObsSensors()

            # Enviamos los valores
            send_obs_val = []

            for i in car_obs_sensor_value:
                send_obs_val.append(i[0])

            self.sceneCarObsSensorChanged.emit(send_obs_val)

            # Leemos el valor de los sensores de vision
            car_sensor_value = self.carReadSensors()

            if car_sensor_value:
                self.sceneCarSensorChanged.emit(car_sensor_value)

            # Maquina de estados
            # El auto esta siguiendo la linea
            if car_state == 0:
                if car_sensor_value:
                    # PID
                    lin_error = 0.0

                    lin_error += car_sensor_value[0] * lin_pond[0] - car_sensor_value[7] * lin_pond[7]
                    lin_error += car_sensor_value[1] * lin_pond[1] - car_sensor_value[6] * lin_pond[6]
                    lin_error += car_sensor_value[2] * lin_pond[2] - car_sensor_value[5] * lin_pond[5]
                    lin_error += car_sensor_value[3] * lin_pond[3] - car_sensor_value[4] * lin_pond[4]

                    lin_error /= 8

                    lin_correction = lin_error * lin_kp + lin_error_prev * lin_kd

                    self.carSetVelocity(-0.5, -lin_correction)

                    # Valores del PID
                    lin_error_prev = lin_error

                # El auto detecto un obtaculo en frente
                for i in range(10, 16):
                    if car_obs_sensor_value[i][1]:
                        if car_obs_sensor_value[i][0] < 0.15:
                            self.carSetVelocity(0.0, 0.0)

                            car_state = 1

            # El auto tiene que analizar hacia que lado evitara el obstaculo
            elif car_state == 1:
                direction_pond = [-2.0, -1.0, -0.5, 0.5, 1.0, 2.0]

                # Direccion
                direction = 0.0

                direction += car_obs_sensor_value[10][0] * direction_pond[0] + \
                             car_obs_sensor_value[15][0] * direction_pond[5]
                direction += car_obs_sensor_value[11][0] * direction_pond[1] + \
                             car_obs_sensor_value[14][0] * direction_pond[4]
                direction += car_obs_sensor_value[12][0] * direction_pond[2] + \
                             car_obs_sensor_value[13][0] * direction_pond[3]

                # El auto deberia esquivar por la izquierda
                if direction < 0.0:
                    car_state = 10

                else:
                    car_state = 20

            # Esquivando el auto por la derecha
            elif car_state == 10:
                self.carSetVelocity(0.0, 0.5)
                time.sleep(5.0)

                car_state = 11

            elif car_state == 11:
                # P
                box_error = 0.125 - car_obs_sensor_value[7][0]
                box_correction = box_error * box_kp
                self.carSetVelocity(-0.5, box_correction)

                # Medicion de la linea
                for i in range(0, 8):
                    if car_sensor_value[i] < 0.2:
                        car_state = 0

            # Esquivando el auto por la izquierda
            elif car_state == 20:
                self.carSetVelocity(0.0, -0.5)
                time.sleep(5.0)

                car_state = 21

            elif car_state == 21:
                # P
                box_error = 0.125 - car_obs_sensor_value[7][0]
                box_correction = box_error * box_kp
                self.carSetVelocity(-0.5, -box_correction)

                # Medicion de la linea
                for i in range(0, 8):
                    if car_sensor_value[i] < 0.2:
                        car_state = 0

            car_sim = sim.simxGetInt32Signal(self.coppeliaIdClient, "car_finished", sim.simx_opmode_oneshot)

            if car_sim[0] == sim.simx_return_ok:
                if car_sim[1] == 1:
                    self.carSetVelocity(0.0, 0.0)
                    car_complete = True

            time.sleep(0.1)

        # Brazo 2
        # Posicionamos las cajas como corresponde
        sim.simxSetInt32Signal(self.coppeliaIdClient, "box_teleport", 1, sim.simx_opmode_oneshot)

        # Movimientos del brazo 2
        for i in range(0, len(self.abb2Trajectory)):
            self.abb2ExecLinTrajectory(self.abb2Trajectory[i][0], self.abb2Trajectory[i][1],
                                       self.abb2Trajectory[i][2], self.abb2Trajectory[i][3],
                                       self.abb2Trajectory[i][4], self.abb2Trajectory[i][5])

            if (i == 0) or (i == 13) or (i == 28):
                self.coppeliaABB2SetVacuumGripper(1)

            if (i == 8) or (i == 22) or (i == 34):
                self.coppeliaABB2SetVacuumGripper(0)

    def abb1ExecTrajectory(self, traj):
        for i in traj:
            self.coppeliaABB1SetJointsRotations(i[0], i[1], i[2], i[3], i[4], i[5])
            self.abb1_q1, self.abb1_q2, self.abb1_q3, self.abb1_q4, self.abb1_q5, self.abb1_q6 = \
                i[0], i[1], i[2], i[3], i[4], i[5]
            time.sleep(0.01)

    def abb1ExecPtpTrajectory(self, x, y, z, a, b, c):
        traj = self.abb1Engine.ptpTrajectory(self.abb1_x, self.abb1_y, self.abb1_z,
                                             self.abb1_a, self.abb1_b, self.abb1_c,
                                             x, y, z, a, b, c,
                                             100)

        if traj:
            self.abb1ExecTrajectory(traj)
            self.abb1_x, self.abb1_y, self.abb1_z, self.abb1_a, self.abb1_b, self.abb1_c = x, y, z, a, b, c

            return True

        else:
            return False

    def abb1ExecLinTrajectory(self, x, y, z, a, b, c):
        traj = self.abb1Engine.linTrajectory(self.abb1_x, self.abb1_y, self.abb1_z,
                                             self.abb1_a, self.abb1_b, self.abb1_c,
                                             x, y, z, a, b, c,
                                             100)

        if traj:
            self.abb1ExecTrajectory(traj)
            self.abb1_x, self.abb1_y, self.abb1_z, self.abb1_a, self.abb1_b, self.abb1_c = x, y, z, a, b, c

            return True

        else:
            return False

    def abb1LoadTrajectory(self):
        file = open("trajs/ABB1.traj", "r")

        if file:
            a = list(csv.reader(file, delimiter=','))

            self.abb1Trajectory.clear()

            for row in a:
                _row = []

                for i in range(0, 6):
                    _row.append(float(row[i]))

                self.abb1Trajectory.append(_row)

            file.close()

            return True

        else:
            return False

    def abb2ExecTrajectory(self, traj):
        for i in traj:
            self.coppeliaABB2SetJointsRotations(i[0], i[1], i[2], i[3], i[4], i[5])
            self.abb2_q1, self.abb2_q2, self.abb2_q3, self.abb2_q4, self.abb2_q5, self.abb2_q6 = \
                i[0], i[1], i[2], i[3], i[4], i[5]
            time.sleep(0.01)

    def abb2ExecPtpTrajectory(self, x, y, z, a, b, c):
        traj = self.abb2Engine.ptpTrajectory(self.abb2_x, self.abb2_y, self.abb2_z,
                                             self.abb2_a, self.abb2_b, self.abb2_c,
                                             x, y, z, a, b, c,
                                             100)

        if traj:
            self.abb2ExecTrajectory(traj)
            self.abb2_x, self.abb2_y, self.abb2_z, self.abb2_a, self.abb2_b, self.abb2_c = x, y, z, a, b, c

            return True

        else:
            return False

    def abb2ExecLinTrajectory(self, x, y, z, a, b, c):
        traj = self.abb2Engine.linTrajectory(self.abb2_x, self.abb2_y, self.abb2_z,
                                             self.abb2_a, self.abb2_b, self.abb2_c,
                                             x, y, z, a, b, c,
                                             100)

        if traj:
            self.abb2ExecTrajectory(traj)
            self.abb2_x, self.abb2_y, self.abb2_z, self.abb2_a, self.abb2_b, self.abb2_c = x, y, z, a, b, c

            return True

        else:
            return False

    def abb2LoadTrajectory(self):
        file = open("trajs/ABB2.traj", "r")

        if file:
            a = list(csv.reader(file, delimiter=','))

            self.abb2Trajectory.clear()

            for row in a:
                _row = []

                for i in range(0, 6):
                    _row.append(float(row[i]))

                self.abb2Trajectory.append(_row)

            file.close()

            return True

        else:
            return False

    def carReadSensors(self):
        values = []

        # Valores de los sensores
        for i in range(0, 8):
            ret_code, res, value = sim.simxReadVisionSensor(self.coppeliaIdClient,
                                                            self.carSensor[i],
                                                            sim.simx_opmode_streaming)
            if value:
                _val = np.mean(value[0])

                values.append(_val)

            else:
                values.append(0.0)

        return values

    def carReadObsSensors(self):
        values = []

        # Valores de los sensores
        for i in range(0, 16):
            ret = sim.simxReadProximitySensor(
                self.coppeliaIdClient, self.carObsSensor[i], sim.simx_opmode_oneshot
            )

            val = np.sqrt(np.power(ret[2][0], 2) + np.power(ret[2][1], 2) + np.power(ret[2][2], 2))

            if val > 1.0:
                val = 1.0

            values.append([val, ret[1]])

        return values

    def carSetVelocity(self, vel, omega):
        # Obtenemos la velocidad de cada rueda para las velocidades
        omega_left = (vel - self.carWheelsAxisDistanceL * omega) / self.carWheelsRadius
        omega_right = (vel + self.carWheelsAxisDistanceL * omega) / self.carWheelsRadius

        self.coppeliaCarSetWheelsVelocity(omega_left, omega_right)

    def coppeliaConnect(self, ip, port):
        # Cerramos todas las conexiones
        self.coppeliaDisconnect()

        # Guardamos el ID de conexión
        self.coppeliaIdClient = sim.simxStart(ip, port, True, True, 1000, 5)

        return self.coppeliaIdClient != -1

    def coppeliaDisconnect(self):
        sim.simxFinish(-1)

    def coppeliaIsConnected(self):
        return sim.simxGetConnectionId(self.coppeliaIdClient) != -1

    def coppeliaBegin(self):
        # Para el brazo 1
        for i in range(0, 6):
            ret_code, self.coppeliaABB1JointHandle[i] = sim.simxGetObjectHandle(
                self.coppeliaIdClient, "IRB120_joint_" + str(i + 1), sim.simx_opmode_blocking
            )

            if ret_code != sim.simx_return_ok:
                return False

        # Para el brazo 2
        for i in range(0, 6):
            ret_code, self.coppeliaABB2JointHandle[i] = sim.simxGetObjectHandle(
                self.coppeliaIdClient, "IRB120_joint_" + str(i + 1) + "_2", sim.simx_opmode_blocking
            )

            if ret_code != sim.simx_return_ok:
                return False

        # Handles del autito
        ret_code, self.carLeftJoint = sim.simxGetObjectHandle(self.coppeliaIdClient,
                                                              "Pioneer_p3dx_leftMotor",
                                                              sim.simx_opmode_blocking)

        if ret_code != sim.simx_return_ok:
            return False

        ret_code, self.carRigthJoint = sim.simxGetObjectHandle(self.coppeliaIdClient,
                                                               "Pioneer_p3dx_rightMotor",
                                                               sim.simx_opmode_blocking)

        if ret_code != sim.simx_return_ok:
            return False

        # Handles de los sensores de linea
        for i in range(0, 8):
            ret_code, self.carSensor[i] = sim.simxGetObjectHandle(self.coppeliaIdClient,
                                                                  "Sensor" + str(i),
                                                                  sim.simx_opmode_blocking)

            if ret_code != sim.simx_return_ok:
                return False

        # Handles de los sensores de parede
        for i in range(0, 16):
            ret_code, self.carObsSensor[i] = sim.simxGetObjectHandle(
                self.coppeliaIdClient,
                "Pioneer_p3dx_ultrasonicSensor" + str(i + 1),
                sim.simx_opmode_blocking
            )

            if ret_code != sim.simx_return_ok:
                return False

        return True

    def coppeliaABB1SetJointsRotations(self, q1, q2, q3, q4, q5, q6):
        sim.simxSetJointPosition(self.coppeliaIdClient, self.coppeliaABB1JointHandle[0],
                                 q1, sim.simx_opmode_oneshot)
        sim.simxSetJointPosition(self.coppeliaIdClient, self.coppeliaABB1JointHandle[1],
                                 q2, sim.simx_opmode_oneshot)
        sim.simxSetJointPosition(self.coppeliaIdClient, self.coppeliaABB1JointHandle[2],
                                 q3, sim.simx_opmode_oneshot)
        sim.simxSetJointPosition(self.coppeliaIdClient, self.coppeliaABB1JointHandle[3],
                                 q4, sim.simx_opmode_oneshot)
        sim.simxSetJointPosition(self.coppeliaIdClient, self.coppeliaABB1JointHandle[4],
                                 q5, sim.simx_opmode_oneshot)
        sim.simxSetJointPosition(self.coppeliaIdClient, self.coppeliaABB1JointHandle[5],
                                 q6, sim.simx_opmode_oneshot)

    def coppeliaABB2SetJointsRotations(self, q1, q2, q3, q4, q5, q6):
        sim.simxSetJointPosition(self.coppeliaIdClient, self.coppeliaABB2JointHandle[0],
                                 q1, sim.simx_opmode_oneshot)
        sim.simxSetJointPosition(self.coppeliaIdClient, self.coppeliaABB2JointHandle[1],
                                 q2, sim.simx_opmode_oneshot)
        sim.simxSetJointPosition(self.coppeliaIdClient, self.coppeliaABB2JointHandle[2],
                                 q3, sim.simx_opmode_oneshot)
        sim.simxSetJointPosition(self.coppeliaIdClient, self.coppeliaABB2JointHandle[3],
                                 q4, sim.simx_opmode_oneshot)
        sim.simxSetJointPosition(self.coppeliaIdClient, self.coppeliaABB2JointHandle[4],
                                 q5, sim.simx_opmode_oneshot)
        sim.simxSetJointPosition(self.coppeliaIdClient, self.coppeliaABB2JointHandle[5],
                                 q6, sim.simx_opmode_oneshot)

    def coppeliaABB1SetVacuumGripper(self, active):
        sim.simxSetInt32Signal(self.coppeliaIdClient, "VacuumGripper_active", active, sim.simx_opmode_oneshot)

    def coppeliaABB2SetVacuumGripper(self, active):
        sim.simxSetInt32Signal(self.coppeliaIdClient, "VacuumGripperB_active", active, sim.simx_opmode_oneshot)

    def coppeliaCarSetWheelsVelocity(self, wheel_left, wheel_right):
        sim.simxSetJointTargetVelocity(self.coppeliaIdClient, self.carLeftJoint, wheel_left,
                                       sim.simx_opmode_oneshot)
        sim.simxSetJointTargetVelocity(self.coppeliaIdClient, self.carRigthJoint, wheel_right,
                                       sim.simx_opmode_oneshot)


text:
import math
import numpy as np

def absolute(column):
    i = 0
    column= convert_num_col(column)
    result = list()
    while i < len(column):
        result.insert(i+1,abs(column[i]))
        i+=1

    return result

def cbrt(column):
    i = 0
    column= convert_num_col(column)
    result = list()
    while i < len(column):
        result.insert(i+1,round(column[i]**(1/3.),2))
        i+=1

    return result

def ceil(column):

    i = 0
    column= convert_num_col(column)
    result = list()
    while i < len(column):
        result.insert(i+1,math.ceil(column[i]))
        i+=1

    return result

def ceiling(column):
    return ceil(column)

def degrees(column):
    i = 0
    column= convert_num_col(column)
    result = list()
    while i < len(column):
        result.insert(i+1,math.degrees(column[i]))
        i+=1

    return result

def div(column1,column2):
   return div_columns(column1,column2)

def exp(column):
    i = 0
    column= convert_num_col(column)
    result = list()
    while i < len(column):
        result.insert(i+1,math.exp(column[i]))
        i+=1

    return result

def factorial(column):
    i = 0
    column= convert_num_col(column)
    result = list()
    while i < len(column):
        result.insert(i+1,math.factorial(column[i]))
        i+=1

    return result

def floor(column):
    i = 0
    column= convert_num_col(column)
    result = list()
    while i < len(column):
        result.insert(i+1,math.floor(column[i]))
        i+=1

    return result

def gcd(column1, column2):
    i = 0
    column= convert_num_col(column)
    result = list()
    if (len(column1)==len(column2)):
        while i < len(column1):
                result.insert(i+1,math.gcd(column1[i],column2[i]))
                i+=1

    return result

def lcm(column1, column2):
    i = 0
    column= convert_num_col(column)
    result = list()
    if (len(column1)==len(column2)):
        while i < len(column1):
                result.insert(i+1,abs(column1[i]*column2[i]) // math.gcd(column1[i],column2[i]))
                i+=1

    return result

def ln(column):
    i = 0
    column= convert_num_col(column)
    result = list()
    while i < len(column):
        result.insert(i+1,math.log(column[i]))
        i+=1

    return result

def log(column):

    return log10(column)

def log10(column):
    return log(column,10)

def log(column,base):
    i = 0
    column= convert_num_col(column)
    result = list()
    while i < len(column):
        result.insert(i+1,math.log(column[i],base))
        i+=1

    return result

def mod(column1, column2):
    return mod_columns(column1,column2)

def pi():
    return math.pi

def pow(column1, column2):
    i = 0
    column= convert_num_col(column)
    result = list()
    if (len(column1)==len(column2)):
        while i < len(column1):
            result.insert(i+1,math.pow(column1[i],column2[i]))
            i+=1

    return result

def radians(column):
    i = 0
    column= convert_num_col(column)
    result = list()
    while i < len(column):
        result.insert(i+1,math.radians(column[i]))
        i+=1

    return result

def random():
    return random()

def sign(column):
    return np.sign(column)

def round(column):
    i = 0
    column= convert_num_col(column)
    result = list()
    while i < len(column):
        result.insert(i+1,round(column[i]))
        i+=1

    return result

def sqrt(column):
    i = 0
    column= convert_num_col(column)
    result = list()
    while i < len(column):
        result.insert(i+1,math.sqrt(column[i]))
        i+=1

    return result

def with_bucket(expresion, rango_izq, rango_der, number_buckets):
    if(rango_izq!=rango_der):
        if(rango_izq==expresion):
            return 0
        if(rango_der==expresion):
            return number_buckets+1
        incremento = (rango_der - rango_izq)/number_buckets
        valor = ((expresion-rango_izq)/incremento)+1
        if(valor>number_buckets ):
            return number_buckets+1
        elif(valor<0):
            return 0
        else:
            return int(valor)

def truncate_col(column,decimals=0):
    i = 0
    column= convert_num_col(column)
    result = list()
    while i < len(column):
        result.insert(i+1,truncate(column[i],decimals))
        i+=1

    return result

def truncate(number, decimals=0):
    if decimals == 0:
        return math.trunc(number)

    factor = 10.0 ** decimals
    return math.trunc(number * factor) / factor

def sum_columns(column1, column2):
    i = 0
    column1= convert_num_col(column1)
    column2= convert_num_col(column2)
    result = list()
    if (len(column1)==len(column2)):
        while i < len(column1):
            result.insert(i+1,column1[i]+column2[i])
            i+=1

    return result

def rest_columns(column1, column2):
    i = 0
    column1= convert_num_col(column1)
    column2= convert_num_col(column2)
    result = list()
    if (len(column1)==len(column2)):
        while i < len(column1):
            result.insert(i+1,column1[i]-column2[i])
            i+=1

    return result

def mult_columns(column1, column2):
    i = 0
    column1= convert_num_col(column1)
    column2= convert_num_col(column2)
    result = list()
    if (len(column1)==len(column2)):
        while i < len(column1):
            result.insert(i+1,column1[i]*column2[i])
            i+=1

    return result

def div_columns(column1, column2):
    i = 0
    column1= convert_num_col(column1)
    column2= convert_num_col(column2)
    result = list()
    if (len(column1)==len(column2)):
        while i < len(column1):
            result.insert(i+1,column1[i]/column2[i])
            i+=1

    return result

def mod_columns(column1, column2):
    i = 0
    column1= convert_num_col(column1)
    column2= convert_num_col(column2)
    result = list()
    if (len(column1)==len(column2)):
        while i < len(column1):
            result.insert(i+1,column1[i]%column2[i])
            i+=1

    return result

def convert_num_col(num):
    if not isinstance(num,int):
       return num
    if isinstance(num,int):
        result =  [num]
        return result

def exp_columns(column1, column2):
    i = 0
    column1= convert_num_col(column1)
    column2= convert_num_col(column2)
    result = list()
    if (len(column1)==len(column2)):
        while i < len(column1):
            result.insert(i+1,column1[i]**column2[i])
            i+=1

    return result

max_stars_repo_name: susanliss/tytus

max_issues_repo_name: susanliss/tytus

max_forks_repo_name: susanliss/tytus


text:
# Process: Landsat 8 multi-spectral imagery

# Import required packages
import os, sys, subprocess
import shutil
import requests
import gzip
import pandas as pd
import datetime
import numpy as np
from math import floor, ceil
import matplotlib.pyplot as plt

# Import helper functions relevant to this script
sys.path.append('E:/mdm123/D/scripts/geo/')
from geo_helpers import extract_projection_info, get_geotiff_props, get_geotiff_projection, get_geotiff_nodatavalue, geotiff_to_array, array_to_geotiff, ls_subfolder_valid, ls8_spectral_index, create_bounded_geotiff

# List paths to GDAL scripts
gdal_warp = 'C:/Anaconda3/envs/geo/Library/bin/gdalwarp.exe'
gdal_calc = 'C:/Anaconda3/envs/geo/Scripts/gdal_calc.py'
gdal_merge = 'C:/Anaconda3/envs/geo/Scripts/gdal_merge.py'
gdal_buildvrt = 'C:/Anaconda3/envs/geo/Library/bin/gdalbuildvrt.exe'

# Define path to relevant folders
folder_srtm = 'E:/mdm123/D/data/DSM/SRTM/'
folder_dtm = 'E:/mdm123/D/data/DTM/proc/'
folder_ls8 = 'E:/mdm123/D/data/Landsat8/'
folder_archive = 'V:/mdm123/Landsat8/'
folder_fig = 'E:/mdm123/D/figures/'

# Define list of zones to be processed (separate LiDAR coverage areas)
zones = ['MRL18_WPE', 'MRL18_WVL', 'MRL18_WKW', 'MRL18_FGA', 'TSM17_STA', 'TSM17_LDM', 'TSM17_GLB', 'TSM16_ATG']

# Define dictionary to hold information relating to each zone covered by the Marlborough (2018) survey
dtm_dict = {'MRL18_WPE':{'label':'Wairau Plains East (Marlborough 2018)',
                         'wrs_prs':['073089'],
                         'survey_start':datetime.date(2018,7,28),
                         'survey_end':datetime.date(2018,8,9)},
            'MRL18_WVL':{'label':'Wairau Valley (Marlborough 2018)',
                         'wrs_prs':['073089','074089'],
                         'survey_start':datetime.date(2018,7,28),
                         'survey_end':datetime.date(2018,8,9)},
            'MRL18_WKW':{'label':'Picton - Waikawa (Marlborough 2018)',
                         'wrs_prs':['073089'],
                         'survey_start':datetime.date(2018,5,26),
                         'survey_end':datetime.date(2018,9,12)},
            'MRL18_FGA':{'label':'Flaxbourne, Grassmere & Lower Awatere (Marlborough 2018)',
                         'wrs_prs':['073089'],
                         'survey_start':datetime.date(2018,5,26),
                         'survey_end':datetime.date(2018,9,12)},
            'TSM17_STA':{'label':'St Arnaud (Tasman 2017)',
                         'wrs_prs':['073089','074089'],
                         'survey_start':datetime.date(2017,11,24),
                         'survey_end':datetime.date(2017,12,16)},
            'TSM17_LDM':{'label':'Lee Dam (Tasman 2017)',
                         'wrs_prs':['073089','074089'],
                         'survey_start':datetime.date(2017,11,24),
                         'survey_end':datetime.date(2017,12,16)},
            'TSM17_GLB':{'label':'Golden Bay & Farewell Spit (Tasman 2017)',
                         'wrs_prs':['074088'],
                         'survey_start':datetime.date(2017,11,24),
                         'survey_end':datetime.date(2017,12,16)},
            'TSM16_ATG':{'label':'Abel Tasman & Golden Bay (Tasman 2016)',
                         'wrs_prs':['074088','074089'],
                         'survey_start':datetime.date(2016,12,13),
                         'survey_end':datetime.date(2016,12,14)}}

# Define the range of survey windows to be tested (the duration of the period for which Ls8 data should be processed, centred on the SRTM survey period)
query_windows = [60, 90, 120, 150, 180]   # Number of days

# Define the number of cells of padding to add along each raster boundary
pad = 44

# Define the no_data value to be used for GeoTIFF creation
no_data_value = -9999

###############################################################################
# 1. Archive to Public folder any scenes for which cloud cover (land) > 85%   #
###############################################################################

# Define a threshold for the maximum cloud coverage (over land) to be allowed
cloud_cover_max = 85

# Download the up-to-date Landsat 8 metadata from USGS
ls8_metadata_url = 'https://landsat.usgs.gov/landsat/metadata_service/bulk_metadata_files/LANDSAT_8_C1.csv.gz'
ls8_metadata_read = requests.get(ls8_metadata_url)
ls8_metadata_path = '{}metadata/{}'.format(folder_ls8, ls8_metadata_url.split('/')[-1])
with open(ls8_metadata_path, 'wb') as ls8_metadata_write:
    ls8_metadata_write.write(ls8_metadata_read.content)

# Extract the compressed file (.gz) & read the CSV inside into a pandas dataframe
with gzip.open(ls8_metadata_path) as ls8_metadata_csv:
    ls8_metadata = pd.read_csv(ls8_metadata_csv)

# Initialise a count of valid & invalid scenes (wrt cloud cover over land)
valid_scenes = 0
invalid_scenes = 0

# Loop through each subfolder in the Level-2 folder, checking its cloud cover over land
print('Archiving scenes for which cloud cover (over land) > {}%:'.format(cloud_cover_max))
for local_subfolder in os.listdir('{}raw/Level2'.format(folder_ls8)):

    # Process further if path is to a folder (rather than a file)
    if os.path.isdir('{}raw/Level2/{}'.format(folder_ls8, local_subfolder)):

        # Check the cloud cover (over land) associated with that scene
        cloud_cover = ls8_metadata.loc[ls8_metadata['LANDSAT_PRODUCT_ID']==local_subfolder, 'CLOUD_COVER_LAND'].values[0]

        # Update scene counts based on whether cloud cover exceeds threshold or not
        if cloud_cover > cloud_cover_max:

            print(' - {} ({:.1f}%):'.format(local_subfolder, cloud_cover), end=' ')

            # Update count of invalid scenes
            invalid_scenes += 1

            # Level-2 subfolder: move from Local to Public archive folder
            l2_local = '{}raw/Level2/{}'.format(folder_ls8, local_subfolder)
            l2_archive = '{}raw/Level2/{}'.format(folder_archive, local_subfolder)
            shutil.move(l2_local, l2_archive)
            print('Level-2', end=' ')

            # Level-1: move the corresponding subfolder from Local to Public archive folder too
            l1_local = '{}raw/Level1/{}'.format(folder_ls8, local_subfolder)
            l1_archive = '{}raw/Level1/{}'.format(folder_archive, local_subfolder)
            shutil.move(l1_local, l1_archive)
            print('Level-1')

        else:
            valid_scenes += 1

print('\n{} scenes retained, with {} scenes rejected (cloud cover over land > {}%) & archived'.format(valid_scenes, invalid_scenes, cloud_cover_max))

###############################################################################
# 2. Archive redundant files to Public folder (to save space on local drive)  #
###############################################################################

# Bands 1-7 are available as Level-2 products (surface reflectance) so aren't required in Level-1 folders
for local_subfolder in os.listdir('{}raw/Level1'.format(folder_ls8)):

    # Process further if path is to a folder (rather than a file)
    if os.path.isdir('{}raw/Level1/{}'.format(folder_ls8, local_subfolder)):

        # If a corresponding subfolder doesn't exist in the Public folder, create one
        archive_path = '{}raw/Level1/{}'.format(folder_archive, local_subfolder)
        if not os.path.exists(archive_path): os.mkdir(archive_path)

        # Move rasters corresponding to Bands 1-7
        for f in os.listdir('{}raw/Level1/{}'.format(folder_ls8, local_subfolder)):
            if f.lower().endswith('b1.tif') or f.lower().endswith('b2.tif') or f.lower().endswith('b3.tif') or f.lower().endswith('b4.tif') or f.lower().endswith('b5.tif') or f.lower().endswith('b6.tif') or f.lower().endswith('b7.tif'):
                shutil.move('{}raw/Level1/{}/{}'.format(folder_ls8, local_subfolder, f), '{}/{}'.format(archive_path, f))

###############################################################################
# 3. Develop cloud-free composites (original grids) for all temporal windows  #
###############################################################################

# Define the QA pixel values considered acceptable (relevant to SR Bands 1-7 of Level 2 products). See https://www.usgs.gov/media/files/land-surface-reflectance-code-lasrc-product-guide
pqa_acceptable_values_hq = [322, 386, 834, 898, 1346, 324, 388, 836, 900, 1348]        # High quality filter: only entries for "Clear" & "Water"
pqa_acceptable_values_lq = pqa_acceptable_values_hq + [328, 352, 416, 480]             # Low quality filter: also allow low-confidence cloud codes

# Define the BQA pixel values considered acceptable (relevant to Bands 8-11 of Level 1 products). See https://www.usgs.gov/land-resources/nli/landsat/landsat-collection-1-level-1-quality-assessment-band?qt-science_support_page_related_con=0#qt-science_support_page_related_con
bqa_acceptable_values_hq = [2720, 2724, 2728, 2732]                                    # High quality filter: only entries for "Clear"
bqa_acceptable_values_lq = bqa_acceptable_values_hq + [2752, 2976, 3008, 3744, 3776]   # Low quality filter: also allow low/medium cloud/shadow/cirrus codes (with no radiation saturation)

# Loop through each LiDAR survey zone
for zone in zones:

    print('\nProcessing Landsat 8 imagery for {}...'.format(dtm_dict[zone]['label']))

    # Create necessary result folder if it doesn't yet exist
    folder_original_zone = '{}proc/cloudfree/Original/{}/'.format(folder_ls8, zone)
    if not os.path.exists(folder_original_zone): os.mkdir(folder_original_zone)

    # Define path to the padded GeoTIFF describing that zone's SRTM coverage
    srtm_tif = '{}proc/{}/SRTM_{}_Z_Pad44.tif'.format(folder_srtm, zone, zone)
    srtm_proj = get_geotiff_projection(srtm_tif)

    # Define which Ls8 WRS path-rows are relevant for that zone
    zone_wrs_prs = dtm_dict[zone]['wrs_prs']

    # Loop through the list of possible query windows
    for query_window in query_windows:

        print(' - Using a {} day search window'.format(query_window))

        # Define the period over which Landsat 8 data should be sought (centred on each survey's collection period)
        query_window = datetime.timedelta(days=query_window)

        # Determine start & end dates of that zone's LiDAR survey period
        survey_start = dtm_dict[zone]['survey_start']
        survey_end = dtm_dict[zone]['survey_end']
        survey_window = survey_end - survey_start

        # Determine search window for tiles, by buffering the start & end dates retrieved above
        query_start = survey_start - (query_window - survey_window)/2
        query_end = survey_end + (query_window - survey_window)/2

        # Loop through all bands, setting parameters as appropriate (based on whether band is from Level-1 or Level-2 collection)
        for band in range(1,12):

            # Define collection level & folder containing raw data for processing & lists of QA codes for processing
            level = 2 if band <= 7 else 1
            folder = '{}raw/Level{}/'.format(folder_ls8, level)
            print('   - Band {} (Level-{} data):'.format(band, level), end=' ')

            # Assign other raster processing parameters based on collection level
            if level == 2:
                # Level-2 collection
                band_prefix = 'sr_band'
                qa_suffix = 'pixel_qa'
                qa_acceptable_values_hq, qa_acceptable_values_lq = pqa_acceptable_values_hq, pqa_acceptable_values_lq
                valid_from = 0.
                valid_to = 10000.
            elif level == 1:
                # Level-1 collection
                band_prefix = 'B'
                qa_suffix = 'BQA'
                qa_acceptable_values_hq, qa_acceptable_values_lq = bqa_acceptable_values_hq, bqa_acceptable_values_lq
                valid_from = 0.
                valid_to = 65535.

            # Get list of rasters available in selected folder
            subfolders = os.listdir(folder)

            # Loop through subfolders, clipping each based on extent info extracted above, masking based on QA raster, and adding to a list of masked arrays

            # Initialise lists to hold masked rasters relating to that band
            band_masked_arrays_hq = []
            band_masked_arrays_lq = []

            # Initialise a count of the number of valid layers available
            n_valid_layers = 0

            # Loop through each available subfolder, extracting relevant data if it's within search window
            for subfolder in subfolders:

                # Process data if capture date is within search window
                if ls_subfolder_valid(subfolder, zone_wrs_prs, query_start, query_end):

                    # Update the count of valid layers
                    n_valid_layers += 1

                    # Define paths to the band raster & its associated QA raster
                    band_path = '{}{}/{}_{}{}.tif'.format(folder, subfolder, subfolder, band_prefix, band)
                    qa_path = '{}{}/{}_{}.tif'.format(folder, subfolder, subfolder, qa_suffix)

                    # Get the projection properties of the band raster
                    band_proj = get_geotiff_projection(band_path)
                    band_props = get_geotiff_props(band_path)

                    # Project the padded GeoTIFF describing that zone's SRTM coverage into the band's CRS
                    srtm_projected = '{}proc/{}/SRTM_{}_Z_Pad44_Projected.tif'.format(folder_srtm, zone, zone)
                    if os.path.exists(srtm_projected): os.remove(srtm_projected)
                    project_command = [gdal_warp, '-s_srs', srtm_proj, '-t_srs', band_proj, srtm_tif, srtm_projected]
                    project_result = subprocess.run(project_command, stdout=subprocess.PIPE)
                    if project_result.returncode != 0:
                        print('\nProcess failed, with error message: {}\n'.format(project_result.stdout))
                        break

                    # Get the projected GeoTIFF's extent (in the CRS used by the input band)
                    srtm_projected_props = get_geotiff_props(srtm_projected)
                    os.remove(srtm_projected)
                    AOI_x_min = srtm_projected_props['x_min']
                    AOI_x_max = srtm_projected_props['x_max']
                    AOI_y_min = srtm_projected_props['y_min']
                    AOI_y_max = srtm_projected_props['y_max']

                    # Snap these extents to match the band raster's grid, adding a buffer of 50 cells along each edge to be safe
                    buffer = 50
                    clip_x_min = band_props['x_min'] + ceil((AOI_x_min - band_props['x_min'])/band_props['res_x'])*band_props['res_x'] - buffer*band_props['res_x']
                    clip_x_max = band_props['x_max'] + floor((AOI_x_max - band_props['x_max'])/band_props['res_x'])*band_props['res_x'] + buffer*band_props['res_x']
                    clip_y_min = band_props['y_min'] + floor((AOI_y_min - band_props['y_min'])/-band_props['res_y'])*-band_props['res_y'] - buffer*-band_props['res_y']
                    clip_y_max = band_props['y_max'] + ceil((AOI_y_max - band_props['y_max'])/-band_props['res_y'])*-band_props['res_y'] + buffer*-band_props['res_y']

                    # Clip the band raster such that its extent matches AOI but its grid alignment matches original band raster
                    band_clip = '{}{}/{}_{}{}_{}.tif'.format(folder, subfolder, subfolder, band_prefix, band, zone)
                    warp_band_command = [gdal_warp, '-overwrite', band_path, band_clip, '-s_srs', band_props['proj'], '-t_srs', band_props['proj'], '-tr', str(band_props['res_x']), str(-band_props['res_y']), '-te', str(clip_x_min), str(clip_y_min), str(clip_x_max), str(clip_y_max), '-te_srs', band_props['proj'], '-r', 'near', '-dstnodata', '-9999']
                    subprocess.call(warp_band_command)

                    # Clip the QA raster too, to be used as the mask
                    qa_clip = '{}{}/{}_{}_{}.tif'.format(folder, subfolder, subfolder, qa_suffix, zone)
                    warp_qa_command = [gdal_warp, '-overwrite', qa_path, qa_clip, '-s_srs', band_props['proj'], '-t_srs', band_props['proj'], '-tr', str(band_props['res_x']), str(-band_props['res_y']), '-te', str(clip_x_min), str(clip_y_min), str(clip_x_max), str(clip_y_max), '-te_srs', band_props['proj'], '-r', 'near', '-dstnodata', '-9999']
                    subprocess.call(warp_qa_command)

                    # Read in arrays of the clipped band raster & the clipped QA raster (to use as a mask)
                    band_array = geotiff_to_array(band_clip)
                    qa_array = geotiff_to_array(qa_clip)

                    # Check the no_data_value used by the input band raster
                    band_nodatavalue = get_geotiff_nodatavalue(band_clip)

                    # Mask the band array wherever the band_nodatavalue is present
                    band_masked_array = np.ma.masked_where(band_array==band_nodatavalue, band_array)

                    # Restrict to valid range & rescale
                    # Replace values below lower bound with no_data_value
                    band_masked_array[band_masked_array < valid_from] = no_data_value
                    # Replace values above upper bound with no_data_value
                    band_masked_array[band_masked_array > valid_to] = no_data_value
                    # Re-mask, using the newly applied no_data_value
                    band_masked_array = np.ma.masked_where(band_masked_array==no_data_value, band_masked_array)
                    # Rescale values by appropriate factor
                    band_masked_array = band_masked_array / valid_to

                    # Mask the band array using the HQ version of the QA array
                    band_masked_array_hq = np.ma.masked_where(~np.isin(qa_array, qa_acceptable_values_hq), band_masked_array)
                    band_masked_arrays_hq.append(band_masked_array_hq)
                    del band_masked_array_hq

                    # Set up a second masked band array, using the LQ version of the QA array (for gap-filling later)
                    band_masked_array_lq = np.ma.masked_where(~np.isin(qa_array, qa_acceptable_values_lq), band_masked_array)
                    band_masked_arrays_lq.append(band_masked_array_lq)
                    del band_masked_array_lq

                    # Clean up arrays no longer needed
                    del band_array, qa_array, band_masked_array

            print('{} images available'.format(n_valid_layers), end=' ')

            # Calculate median of the HQ array list representing data for the selected band
            band_masked_arrays_hq_stack = np.ma.array(band_masked_arrays_hq)
            band_masked_median_hq = np.ma.median(band_masked_arrays_hq_stack, axis=0, overwrite_input=True)
            del band_masked_arrays_hq, band_masked_arrays_hq_stack

            # Calculate median of the LQ array list representing data for the selected band
            band_masked_arrays_lq_stack = np.ma.array(band_masked_arrays_lq)
            band_masked_median_lq = np.ma.median(band_masked_arrays_lq_stack, axis=0)
            del band_masked_arrays_lq, band_masked_arrays_lq_stack

            # Save results in which the LQ results are used to gap-fill the HQ results
            # Retrieve the mask derived from the HQ filter
            mask_hq = np.ma.getmask(band_masked_median_hq)
            # Make copies of the HQ median results
            band_masked_median_gapfill = np.ma.copy(band_masked_median_hq)
            # Try to fill in masked values using results from the LQ results, where possible
            band_masked_median_gapfill[mask_hq] = band_masked_median_lq[mask_hq]
            del band_masked_median_hq, band_masked_median_lq, mask_hq
            # Fill remaining gaps with the same no_data value used previously
            band_masked_median_gapfill = band_masked_median_gapfill.filled(fill_value=no_data_value)

            # Update the dictionary describing band raster properties
            band_props['x_min'] = clip_x_min
            band_props['x_max'] = clip_x_max
            band_props['y_min'] = clip_y_min
            band_props['y_max'] = clip_y_max
            band_props['width'] = int((clip_x_max - clip_x_min)/band_props['res_x'])
            band_props['height'] = int((clip_y_max - clip_y_min)/-band_props['res_y'])

            # Save these new arrays to TIF files
            tif_median_path_gapfill = '{}LS8_L{}_B{}_{}d_GapFill_Median_{}.tif'.format(folder_original_zone, level, band, str(query_window.days).zfill(3), zone)
            array_to_geotiff(band_masked_median_gapfill, tif_median_path_gapfill, no_data_value, band_props)
            del band_masked_median_gapfill, band_props
            print('DONE')

###############################################################################
# 4. Calculate spectral index products, using original composite rasters      #
###############################################################################

# Define the list of spectral index products to be generated
spectral_index_products = ['NDVI','EVI','AVI','SAVI','MSAVI','SI','BSI','NDMI','MNDWI','AWEInsh','AWEIsh','NDBI']

# Loop through all available survey zones
for zone in zones:

    print('\nCalculating Landsat 8 spectral indices for {}...'.format(dtm_dict[zone]['label']))
    folder_ls8_zone = '{}proc/cloudfree/Original/{}'.format(folder_ls8, zone)

    # Loop through all considered query windows
    for query_window in query_windows:

        print(' - Indices based on {}d window:'.format(query_window), end=' ')

        # Define paths to relevant band rasters
        band_2_path = '{}/LS8_L2_B2_{}d_GapFill_Median_{}.tif'.format(folder_ls8_zone, str(query_window).zfill(3), zone)
        band_3_path = '{}/LS8_L2_B3_{}d_GapFill_Median_{}.tif'.format(folder_ls8_zone, str(query_window).zfill(3), zone)
        band_4_path = '{}/LS8_L2_B4_{}d_GapFill_Median_{}.tif'.format(folder_ls8_zone, str(query_window).zfill(3), zone)
        band_5_path = '{}/LS8_L2_B5_{}d_GapFill_Median_{}.tif'.format(folder_ls8_zone, str(query_window).zfill(3), zone)
        band_6_path = '{}/LS8_L2_B6_{}d_GapFill_Median_{}.tif'.format(folder_ls8_zone, str(query_window).zfill(3), zone)
        band_7_path = '{}/LS8_L2_B7_{}d_GapFill_Median_{}.tif'.format(folder_ls8_zone, str(query_window).zfill(3), zone)

        # Read data arrays for each into memory
        b2, b3, b4, b5, b6, b7 = [geotiff_to_array(band_path) for band_path in [band_2_path, band_3_path, band_4_path, band_5_path, band_6_path, band_7_path]]
        b_props = get_geotiff_props(band_2_path)

        # Loop through all spectral index products
        for product in spectral_index_products:
            # Define output path
            product_path = '{}/LS8_{}_{}d_{}.tif'.format(folder_ls8_zone, product, str(query_window).zfill(3), zone)
            # Use Landsat 8 function to calculate the spectral index (including gap-filling)
            with np.errstate(divide='ignore'):  # Divide by zero issues are handled already - suppressing warnings
                ls8_spectral_index(b2, b3, b4, b5, b6, b7, product, b_props, product_path, no_data_value)
            print(product, end=' ')

        # Build a VRT (Virtual Dataset) containing the R,G,B bands - for visualisation as a True Colour Image (TCI)
        vrt_path = '{}/LS8_TCI_{}d_{}.vrt'.format(folder_ls8_zone, str(query_window).zfill(3), zone)
        vrt_command = [gdal_buildvrt, '-separate', vrt_path] + ['{}/LS8_L2_B{}_{}d_GapFill_Median_{}.tif'.format(folder_ls8_zone, band, str(query_window).zfill(3), zone) for band in ['4','3','2']]
        vrt_result = subprocess.run(vrt_command, stdout=subprocess.PIPE)
        if vrt_result.returncode != 0:
            print(vrt_result.stdout)
            break
        print('TCI')

###############################################################################
# 5. Resample cloud-free composites (all time windows) to match SRTM grids    #
###############################################################################

# The 'cubicspline' resampling method was found to provide the best results based on visual comparisons
resampling = 'cubicspline'

# Loop through all zones, resampling all available Landsat 8 imagery to arrays coincident with the SRTM arrays
for zone in zones:

    print('\nProcessing Landsat 8 data for {}...'.format(dtm_dict[zone]['label']))

    # 5a. Read the appropriate SRTM DSM raster into memory & retrieve its properties
    print(' - Analysing zonal SRTM raster to align grids...')
    srtm_filename = '{}proc/{}/SRTM_{}_Z.tif'.format(folder_srtm, zone, zone)
    srtm_proj, srtm_res_x, srtm_res_y, srtm_x_min, srtm_x_max, srtm_y_min, srtm_y_max, srtm_width, srtm_height = extract_projection_info(srtm_filename)

    # Define a new bounding box, including the padding required for the 2D convnet data pre-processing
    pad_x_min = srtm_x_min - pad*srtm_res_x
    pad_x_max = srtm_x_max + pad*srtm_res_x
    pad_y_min = srtm_y_min - pad*-srtm_res_y
    pad_y_max = srtm_y_max + pad*-srtm_res_y
    pad_width = srtm_width + 2*pad
    pad_height = srtm_height + 2*pad

    # Create new folder for resampled rasters, if it doesn't exist already
    folder_ls8_zone = '{}proc/cloudfree/Resampled/{}'.format(folder_ls8, zone)
    if not os.path.exists(folder_ls8_zone):
        os.makedirs(folder_ls8_zone)

    # 5b. Resample each band, setting parameters as appropriate (based on whether band is from Level-1 or Level-2 collection)
    for band in range(1,12):

        # Define collection level
        level = 2 if band <= 7 else 1
        print(' - Band {} (Level-{} data):'.format(band, level), end=' ')

        # Loop through the processed results, for the various time windows trialled
        for query_window in query_windows:

            # Open band raster and extract its coordinate reference system (CRS)
            ls8_tile_path = '{}proc/cloudfree/Original/{}/LS8_L{}_B{}_{}d_GapFill_Median_{}.tif'.format(folder_ls8, zone, level, band, str(query_window).zfill(3), zone)
            ls8_proj = get_geotiff_projection(ls8_tile_path)

            # Warp band raster to WGS84, aligning with padded SRTM raster
            ls8_resample_pad = '{}/LS8_B{}_{}d_{}_Pad44.tif'.format(folder_ls8_zone, band, str(query_window).zfill(3), zone)
            warp_command_pad = [gdal_warp, '-overwrite', ls8_tile_path, ls8_resample_pad, '-s_srs', ls8_proj, '-t_srs', srtm_proj, '-tr', str(srtm_res_x), str(-srtm_res_y), '-te', str(pad_x_min), str(pad_y_min), str(pad_x_max), str(pad_y_max), '-te_srs', srtm_proj, '-r', resampling, '-dstnodata', str(no_data_value)]
            warp_result_pad = subprocess.run(warp_command_pad, stdout=subprocess.PIPE)
            if warp_result_pad.returncode != 0:
                print(warp_result_pad.stdout)
                break

            # Now clip that padded raster to the unpadded extent (i.e. same as SRTM/DTM zone raster)
            ls8_resample = '{}/LS8_B{}_{}d_{}.tif'.format(folder_ls8_zone, band, str(query_window).zfill(3), zone)
            warp_command = [gdal_warp, '-overwrite', ls8_resample_pad, ls8_resample, '-s_srs', srtm_proj, '-t_srs', srtm_proj, '-tr', str(srtm_res_x), str(-srtm_res_y), '-te', str(srtm_x_min), str(srtm_y_min), str(srtm_x_max), str(srtm_y_max), '-te_srs', srtm_proj, '-r', 'near', '-dstnodata', str(no_data_value)]
            warp_result = subprocess.run(warp_command, stdout=subprocess.PIPE)
            if warp_result.returncode != 0:
                print(warp_result.stdout)
                break

        print('DONE')

    # 5c. Resample each index product
    for product in spectral_index_products:

        print(' - {} spectral index:'.format(product), end=' ')

        # Loop through the processed results, for the three time windows trialled
        for query_window in query_windows:

            # Open product raster and extract its coordinate reference system (CRS)
            ls8_tile_path = '{}proc/cloudfree/Original/{}/LS8_{}_{}d_{}.tif'.format(folder_ls8, zone, product, str(query_window).zfill(3), zone)
            ls8_proj = get_geotiff_projection(ls8_tile_path)

            # Warp band raster to WGS84, aligning with padded SRTM raster
            ls8_resample_pad = '{}/LS8_{}_{}d_{}_Pad44.tif'.format(folder_ls8_zone, product, str(query_window).zfill(3), zone)
            warp_command_pad = [gdal_warp, '-overwrite', ls8_tile_path, ls8_resample_pad, '-s_srs', ls8_proj, '-t_srs', srtm_proj, '-tr', str(srtm_res_x), str(-srtm_res_y), '-te', str(pad_x_min), str(pad_y_min), str(pad_x_max), str(pad_y_max), '-te_srs', srtm_proj, '-r', resampling, '-dstnodata', str(no_data_value)]
            warp_result_pad = subprocess.run(warp_command_pad, stdout=subprocess.PIPE)
            if warp_result_pad.returncode != 0:
                print(warp_result_pad.stdout)
                break

            # Now clip that padded raster to the unpadded extent (i.e. same as SRTM/DTM zone raster)
            ls8_resample = '{}/LS8_{}_{}d_{}.tif'.format(folder_ls8_zone, product, str(query_window).zfill(3), zone)
            warp_command = [gdal_warp, '-overwrite', ls8_resample_pad, ls8_resample, '-s_srs', srtm_proj, '-t_srs', srtm_proj, '-tr', str(srtm_res_x), str(-srtm_res_y), '-te', str(srtm_x_min), str(srtm_y_min), str(srtm_x_max), str(srtm_y_max), '-te_srs', srtm_proj, '-r', 'near', '-dstnodata', str(no_data_value)]
            warp_result = subprocess.run(warp_command, stdout=subprocess.PIPE)
            if warp_result.returncode != 0:
                print(warp_result.stdout)
                break

        print('DONE')

###############################################################################
# 6. Compare band/product value distribution across all available zones       #
###############################################################################

# Loop through all available resampled rasters (bands & products)
for raster in ['B{}'.format(b) for b in range(1,12)] + spectral_index_products:

    print('\nProcessing {} raster...'.format(raster), end=' ')

    # Loop through all available temporary search windows
    for query_window in query_windows:

        # Initialise a figure to show that raster's data distribution for all zones
        fig, axes = plt.subplots(nrows=len(zones), sharex=True, figsize=(9,15))

        # Loop through all available zones
        for i, zone in enumerate(zones):

            # Define path to appropriate raster & retrieve its values
            raster_path = '{}proc/cloudfree/Resampled/{}/LS8_{}_{}d_{}.tif'.format(folder_ls8, zone, raster, str(query_window).zfill(3), zone)
            raster_array = geotiff_to_array(raster_path)
            raster_values = [val for val in raster_array.flatten() if val != no_data_value]

            # Show that zone's data distribution using a histogram & label figure with zone name
            axes[i].hist(raster_values, bins=100)
            axes[i].annotate(zone, xy=(0.99, 0.96), xycoords='axes fraction', ha='right', va='top')

            # Show extreme values using dashed red lines
            axes[i].axvline(x=min(raster_values), color='red', linestyle='dashed', alpha=0.3)
            axes[i].axvline(x=max(raster_values), color='red', linestyle='dashed', alpha=0.3)

        # Finalise figure
        fig.suptitle('Distribution of {} values across all zones ({} day window)'.format(raster, query_window), fontweight='bold')
        fig.tight_layout()
        fig.subplots_adjust(top=0.96)
        fig.savefig('{}All/Distributions/Landsat8/{}_{}d.png'.format(folder_fig, raster, str(query_window).zfill(3)), dpi=300)
        plt.close()

        print('{}d'.format(query_window), end=' ')

###############################################################################
# 7. Bound all bands & appropriate spectral indices, to truncate artefacts    #
###############################################################################

# Define a dictionary of bounds for each spectral index product
bounds_dict = {'NDVI':(-1., 1.),       # Source: https://www.usgs.gov/land-resources/nli/landsat/landsat-normalized-difference-vegetation-index
               'EVI':(-1., 1.),        # Source: https://www.usgs.gov/land-resources/nli/landsat/landsat-enhanced-vegetation-index
               'AVI':(-1., 1.),
               'SAVI':(-1., 1.),       # Source: https://www.usgs.gov/land-resources/nli/landsat/landsat-soil-adjusted-vegetation-index
               'MSAVI':(-1., 1.),      # Source: https://www.usgs.gov/land-resources/nli/landsat/landsat-modified-soil-adjusted-vegetation-index
               'SI':(0., 1.),          # Source: based on equation used
               'BSI':(-1., 1.),        # Source: based on equation used
               'NDMI':(-1., 1.),       # Source: https://www.usgs.gov/land-resources/nli/landsat/normalized-difference-moisture-index
               'MNDWI':(-1., 1.),      # Source: Xu 2006
               'AWEInsh':(-7., 4.),    # Source: based on equation used
               'AWEIsh':(-3.25, 3.5),  # Source: based on equation used
               'NDBI':(-1., 1.)}       # Source: based on equation used

# Add bounds for all bands too
for band in range(1,12):
    bounds_dict['B{}'.format(band)] = (0., 1.)

# Loop through all bands & spectral index products
for raster in ['B{}'.format(b) for b in range(1,12)] + spectral_index_products:

    print('\nBounding {} raster...'.format(raster), end=' ')

    # Retrieve the bounds defined for that raster (if any)
    bounds = bounds_dict[raster]
    lower_bound = bounds[0]
    upper_bound = bounds[1]

    # Loop through all available temporary search windows
    for query_window in query_windows:

        # Loop through all available zones
        for i, zone in enumerate(zones):

            # Loop through both versions of each raster (padded & unpadded)
            for padding in ['', '_Pad44']:

                # Define path to that raster GeoTIFF
                raster_unbounded = '{}proc/cloudfree/Resampled/{}/LS8_{}_{}d_{}{}.tif'.format(folder_ls8, zone, raster, str(query_window).zfill(3), zone, padding)
                raster_bounded = '{}proc/cloudfree/Resampled/{}/LS8_{}_{}d_{}{}_Bounded.tif'.format(folder_ls8, zone, raster, str(query_window).zfill(3), zone, padding)

                # Create a bounded GeoTIFF, using the helper function defined
                with np.errstate(invalid='ignore'):
                    create_bounded_geotiff(raster_unbounded, raster_bounded, lower_bound, upper_bound, no_data_value)

        print('{}d'.format(query_window), end=' ')

###############################################################################
# 8. Compare BOUNDED band/product value distribution for all available zones  #
###############################################################################

# Loop through all available resampled rasters (bands & products)
for raster in ['B{}'.format(b) for b in range(1,12)] + spectral_index_products:

    print('\nProcessing {} raster...'.format(raster), end=' ')

    # Loop through all available temporary search windows
    for query_window in query_windows:

        # Initialise a figure to show that raster's data distribution for all zones
        fig, axes = plt.subplots(nrows=len(zones), sharex=True, figsize=(9,15))

        # Loop through all available zones
        for i, zone in enumerate(zones):

            # Define path to appropriate raster & retrieve its values
            raster_path = '{}proc/cloudfree/Resampled/{}/LS8_{}_{}d_{}_Bounded.tif'.format(folder_ls8, zone, raster, str(query_window).zfill(3), zone)
            raster_array = geotiff_to_array(raster_path)
            raster_values = raster_array[raster_array != no_data_value].flatten()

            # Show that zone's data distribution using a histogram & label figure with zone name
            axes[i].hist(raster_values, bins=100, color='green', alpha=0.5)
            axes[i].annotate(zone, xy=(0.99, 0.96), xycoords='axes fraction', ha='right', va='top')

            # Show extreme values using dashed red lines
            axes[i].axvline(x=min(raster_values), color='red', linestyle='dashed', alpha=0.3)
            axes[i].axvline(x=max(raster_values), color='red', linestyle='dashed', alpha=0.3)

        # Finalise figure
        fig.suptitle('Distribution of {} values (bounded) across all zones ({} day window)'.format(raster, query_window), fontweight='bold')
        fig.tight_layout()
        fig.subplots_adjust(top=0.96)
        fig.savefig('{}All/Distributions/Landsat8/{}_{}d_Bounded.png'.format(folder_fig, raster, str(query_window).zfill(3)), dpi=300)
        plt.close()

        print('{}d'.format(query_window), end=' ')

max_stars_repo_name: mdmeadows/DSM-to-DTM

max_issues_repo_name: sarathsrk/DSM-to-DTM

max_forks_repo_name: sarathsrk/DSM-to-DTM


text:
import argparse
import numpy as np
import cPickle as pickle

parser = argparse.ArgumentParser()
parser.add_argument('--vectors_file', default='../glove.6B.300d.txt', type=str)
args = parser.parse_args()

#with open(args.vocab_file, 'r') as f:
#    words = [x.rstrip().split(' ')[0] for x in f.readlines()]
with open(args.vectors_file, 'r') as f:
    vectors = {}
    for line in f:
        vals = line.rstrip().split(' ')
        vectors[vals[0]] = map(float, vals[1:])

embedding_dim = len(vals[1:])
with open('/scratch/costarendon.b/show-attend-and-tell-tensorflow/data/train/word_to_idx.pkl', 'rb') as f:
    word_to_idx = pickle.load(f)

embedding_matrix = np.empty((max(word_to_idx.values())+1, embedding_dim))

count = 0

for word, idx in word_to_idx.items():

    if word in vectors.keys():
        embedding_matrix[idx,:] = np.array(vectors[word])

    else:
        print("Word:", word, "not found")
        count = count + 1
        embedding_matrix[idx,:] = np.random.uniform(low=-1., high=1., size=embedding_dim)

print('Count of words not in glove:', count)
pickle.dump(embedding_matrix, open('embedding_matrix2.pkl', 'wb'), pickle.HIGHEST_PROTOCOL)

max_stars_repo_name: Mogbo/image-caption

max_issues_repo_name: Mogbo/image-caption

max_forks_repo_name: Mogbo/image-caption


text:
# -*- coding: utf-8 -*-

# This code is part of Qiskit.
#
# (C) Copyright IBM 2017, 2019.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

# pylint: disable=len-as-condition

"""
Predicates for operators.
"""

import numpy as np

ATOL_DEFAULT = 1e-8
RTOL_DEFAULT = 1e-5

def matrix_equal(mat1,
                 mat2,
                 ignore_phase=False,
                 rtol=RTOL_DEFAULT,
                 atol=ATOL_DEFAULT):
    """Test if two arrays are equal.

    The final comparison is implemented using Numpy.allclose. See its
    documentation for additional information on tolerance parameters.

    If ignore_phase is True both matrices will be multiplied by
    exp(-1j * theta) where `theta` is the first nphase for a
    first non-zero matrix element `|a| * exp(1j * theta)`.

    Args:
        mat1 (matrix_like): a matrix
        mat2 (matrix_like): a matrix
        ignore_phase (bool): ignore complex-phase differences between
            matrices [Default: False]
        rtol (double): the relative tolerance parameter [Default {}].
        atol (double): the absolute tolerance parameter [Default {}].

    Returns:
        bool: True if the matrices are equal or False otherwise.
    """.format(RTOL_DEFAULT, ATOL_DEFAULT)

    if atol is None:
        atol = ATOL_DEFAULT
    if rtol is None:
        rtol = RTOL_DEFAULT
    mat1 = np.array(mat1)
    mat2 = np.array(mat2)
    if mat1.shape != mat2.shape:
        return False
    if ignore_phase:
        # Get phase of first non-zero entry of mat1 and mat2
        # and multiply all entries by the conjugate
        phases1 = np.angle(mat1[mat1 != 0].ravel(order='F'))
        if len(phases1) > 0:
            mat1 = np.exp(-1j * phases1[0]) * mat1
        phases2 = np.angle(mat2[mat2 != 0].ravel(order='F'))
        if len(phases2) > 0:
            mat2 = np.exp(-1j * phases2[0]) * mat2
    return np.allclose(mat1, mat2, rtol=rtol, atol=atol)

def is_square_matrix(mat):
    """Test if an array is a square matrix."""
    mat = np.array(mat)
    if mat.ndim != 2:
        return False
    shape = mat.shape
    return shape[0] == shape[1]

def is_diagonal_matrix(mat, rtol=RTOL_DEFAULT, atol=ATOL_DEFAULT):
    """Test if an array is a diagonal matrix"""
    if atol is None:
        atol = ATOL_DEFAULT
    if rtol is None:
        rtol = RTOL_DEFAULT
    mat = np.array(mat)
    if mat.ndim != 2:
        return False
    return np.allclose(mat, np.diag(np.diagonal(mat)), rtol=rtol, atol=atol)

def is_symmetric_matrix(op, rtol=RTOL_DEFAULT, atol=ATOL_DEFAULT):
    """Test if an array is a symmetric matrix"""
    if atol is None:
        atol = ATOL_DEFAULT
    if rtol is None:
        rtol = RTOL_DEFAULT
    mat = np.array(op)
    if mat.ndim != 2:
        return False
    return np.allclose(mat, mat.T, rtol=rtol, atol=atol)

def is_hermitian_matrix(mat, rtol=RTOL_DEFAULT, atol=ATOL_DEFAULT):
    """Test if an array is a Hermitian matrix"""
    if atol is None:
        atol = ATOL_DEFAULT
    if rtol is None:
        rtol = RTOL_DEFAULT
    mat = np.array(mat)
    if mat.ndim != 2:
        return False
    return np.allclose(mat, np.conj(mat.T), rtol=rtol, atol=atol)

def is_positive_semidefinite_matrix(mat, rtol=RTOL_DEFAULT, atol=ATOL_DEFAULT):
    """Test if a matrix is positive semidefinite"""
    if atol is None:
        atol = ATOL_DEFAULT
    if rtol is None:
        rtol = RTOL_DEFAULT
    if not is_hermitian_matrix(mat, rtol=rtol, atol=atol):
        return False
    # Check eigenvalues are all positive
    vals = np.linalg.eigvalsh(mat)
    for v in vals:
        if v < -atol:
            return False
    return True

def is_identity_matrix(mat,
                       ignore_phase=False,
                       rtol=RTOL_DEFAULT,
                       atol=ATOL_DEFAULT):
    """Test if an array is an identity matrix."""
    if atol is None:
        atol = ATOL_DEFAULT
    if rtol is None:
        rtol = RTOL_DEFAULT
    mat = np.array(mat)
    if mat.ndim != 2:
        return False
    if ignore_phase:
        # If the matrix is equal to an identity up to a phase, we can
        # remove the phase by multiplying each entry by the complex
        # conjugate of the phase of the [0, 0] entry.
        theta = np.angle(mat[0, 0])
        mat = np.exp(-1j * theta) * mat
    # Check if square identity
    iden = np.eye(len(mat))
    return np.allclose(mat, iden, rtol=rtol, atol=atol)

def is_unitary_matrix(mat, rtol=RTOL_DEFAULT, atol=ATOL_DEFAULT):
    """Test if an array is a unitary matrix."""
    if atol is None:
        atol = ATOL_DEFAULT
    if rtol is None:
        rtol = RTOL_DEFAULT
    mat = np.array(mat)
    # Compute A^dagger.A and see if it is identity matrix
    mat = np.conj(mat.T).dot(mat)
    return is_identity_matrix(mat, ignore_phase=False, rtol=rtol, atol=atol)