Python numpy.zeros() Examples

def polygonize(self, n=73):
        """Gets a approximate polygon array representing the ellipse.

        Note that the last point is the same as the first point, creating a closed
        polygon.

        :param n: The number of points to generate. Default is 73 (one vertex every 5 degrees).
        :type n: int
        :return: An [n  x 2] numpy array, describing the boundary vertices of
                 the polygonized ellipse.
        :rtype: :py:class:`numpy.ndarray`

        """
        t = np.linspace(0, 2 * np.pi, num=n, endpoint=True)
        out = np.zeros((len(t), 2), dtype='float')
        out[:, 0] = (self.center_point[0] +
                     self.radii[0] * np.cos(t) * np.cos(self.rotation_angle) -
                     self.radii[1] * np.sin(t) * np.sin(self.rotation_angle))
        out[:, 1] = (self.center_point[1] +
                     self.radii[0] * np.cos(t) * np.sin(self.rotation_angle) +
                     self.radii[1] * np.sin(t) * np.cos(self.rotation_angle))
        return out 

def get_a_datum(self):
        if self._compressed:
            datum = extract_sample(
                self._data[self._cur], self._mean, self._resize)
        else:
            datum = self._data[self._cur]
        # start parsing labels
        label_elems = parse_label(self._label[self._cur])
        label = np.zeros(self._label_dim)
        if not self._multilabel:
            label[0] = label_elems[0]
        else:
            for i in label_elems:
                label[i] = 1
        self._cur = (self._cur + 1) % self._sample_count
        return datum, label 

def fit_B2AC(points):
    """Ellipse fitting in Python with numerically unstable algorithm.

    Described `here <http://research.microsoft.com/pubs/67845/ellipse-pami.pdf>`_.

    N_POLYPOINTS.B. Do not use, since it works with almost singular matrix.

    :param points: The [Nx2] array of points to fit ellipse to.
    :type points: :py:class:`numpy.ndarray`
    :return: The conic section array defining the fitted ellipse.
    :rtype: :py:class:`numpy.ndarray`

    """
    import scipy.linalg as scla

    constraint_matrix = np.zeros((6, 6))
    constraint_matrix[0, 2] = 2
    constraint_matrix[1, 1] = -1
    constraint_matrix[2, 0] = 2

    S = _calculate_scatter_matrix_py(points[:, 0], points[:, 1])

    evals, evect = scla.eig(S, constraint_matrix)
    ind = np.where(evals == (evals[evals > 0].min()))[0][0]
    return evect[:, ind] 

def get_bounding_box(self):
        """Get a four point Polygon describing the bounding box of the current Polygon.

        :return: A new polygon, describing this one's bounding box.
        :rtype: :py:class:`b2ac.geometry.polygon.B2ACPolygon`

        """
        mins = self.polygon_points.min(axis=0)
        maxs = self.polygon_points.max(axis=0)
        bb_pnts = np.zeros((4, 2), dtype=self.polygon_points.dtype)
        bb_pnts[0, :] = mins
        bb_pnts[1, :] = [mins[0], maxs[1]]
        bb_pnts[2, :] = maxs
        bb_pnts[3, :] = [maxs[0], mins[1]]
        out = B2ACPolygon(bb_pnts)
        return out 

def query(self, x_0):
        if self.contains(x_0):
            return x_0, []

        x_star = self.project(x_0) 
        if self._unqueried:
            # This is an abuse of the word hyperplane; this function
            # actually returns a set of hyperplanes that exactly identifies
            # the zero set. If added to an outer approximation, the
            # interpretation is that we are doing a presolve by forcing
            # x to be 0.
            assert np.prod(self._x.size) == x_0.shape[0]
            A = scipy.sparse.eye(x_0.shape[0])
            # Note that we signal to the dynamic polyhedron that this
            # hyperplane should be pinned. The reasoning is that this
            # "hyperplane" is simply a linear algebraic presolve, and enforcing
            # it should not increase our computational complexity if our
            # implementation is sound.
            h = [Hyperplane(self._x, A, np.zeros(x_0.shape[0]), pin=True)]
            self._unqueried = False
            self._info.append(h)
            return x_star, h
        else:
            return x_star, [] 

def solve(self, problem):
        left_set = problem.sets[0]
        right_set = problem.sets[1]

        iterate = (self._initial_iterate if
            self._initial_iterate is not None else np.ones(problem.dimension))

        # (p_n), (q_n) are the auxiliary sequences, (a_n), (b_n) the main
        # sequences, defined in Bauschke's 98 paper
        # (Dykstra's Alternating Projection Algorithm for Two Sets).
        self.p = [None] * (self.max_iters + 1)
        self.q = [None] * (self.max_iters + 1)
        self.a = [None] * (self.max_iters + 1)
        self.b = [None] * (self.max_iters + 1)
        zero_vector = np.zeros(problem.dimension)
        self.p[0] = self.q[0] = zero_vector
        self.b[0] = self.a[0] = iterate
        residuals = []

        status = Optimizer.Status.INACCURATE
        for n in xrange(1, self.max_iters + 1):
            if self.verbose:
                print 'iteration %d' % n
            # TODO(akshayka): Robust stopping criterion
            residuals.append(self._compute_residual(
                self.b[n-1], left_set, right_set))
            if self.verbose:
                print '\tresidual: %e' % sum(residuals[-1])
            if self._is_optimal(residuals[-1]):
                status = Optimizer.Status.OPTIMAL
                if not self.do_all_iters:
                    break

            self.a[n] = left_set.project(self.b[n-1] + self.p[n-1])
            self.b[n] = right_set.project(self.a[n] + self.q[n-1])
            self.p[n] = self.b[n-1] + self.p[n-1] - self.a[n]
            self.q[n] = self.a[n] + self.q[n-1] - self.b[n]

        # TODO(akshayka): does it matter if I return self.b vs self.a?
        # the first implementation returned self.a ...
        return self.b, residuals, status 

def load_keypoints(image_filepath, image_height, image_width):
    """Load facial keypoints of one image."""
    fp_keypoints = "%s.cat" % (image_filepath,)
    if not os.path.isfile(fp_keypoints):
        raise Exception("Could not find keypoint coordinates for image '%s'." \
                        % (image_filepath,))
    else:
        coords_raw = open(fp_keypoints, "r").readlines()[0].strip().split(" ")
        coords_raw = [abs(int(coord)) for coord in coords_raw]
        keypoints = []
        #keypoints_arr = np.zeros((9*2,), dtype=np.int32)
        for i in range(1, len(coords_raw), 2): # first element is the number of coords
            x = np.clip(coords_raw[i], 0, image_width-1)
            y = np.clip(coords_raw[i+1], 0, image_height-1)
            keypoints.append((x, y))

        return keypoints 

def load_text_and_labels(path, lowercase = True, multilingual = False, distinct_labels_index = None):
	"""Loads text instances from files (one text one line), splits the data into words and generates labels (as one-hot vectors).
		Returns split sentences and labels.
	"""
    # Load data from files
	lines = [(s.lower() if lowercase else s).strip().split() for s in list(codecs.open(path,'r',encoding='utf8', errors='replace').readlines())]
	x_instances = [l[1:-1] for l in lines] if multilingual else [l[:-1] for l in lines]

	if multilingual: 
		langs = [l[0] for l in lines]
	labels = [l[-1] for l in lines]
	
	dist_labels = list(set(labels)) if distinct_labels_index is None else distinct_labels_index
	y_instances = [np.zeros(len(dist_labels)) for l in labels]
	for i in range(len(y_instances)):
		y_instances[i][dist_labels.index(labels[i])] = 1

	return [x_instances, y_instances, langs, dist_labels] if multilingual else [x_instances, y_instances, dist_labels] 

def wer(self, r, h):
        # initialisation
        d = np.zeros((len(r)+1)*(len(h)+1), dtype=np.uint8)
        d = d.reshape((len(r)+1, len(h)+1))
        for i in range(len(r)+1):
            for j in range(len(h)+1):
                if i == 0:
                    d[0][j] = j
                elif j == 0:
                    d[i][0] = i

        # computation
        for i in range(1, len(r)+1):
            for j in range(1, len(h)+1):
                if r[i-1] == h[j-1]:
                    d[i][j] = d[i-1][j-1]
                else:
                    substitution = d[i-1][j-1] + 1
                    insertion    = d[i][j-1] + 1
                    deletion     = d[i-1][j] + 1
                    d[i][j] = min(substitution, insertion, deletion)

        return d[len(r)][len(h)] 

def compute_edge_angles(edges):
    """
    Compute the angles between the edges and the x-axis.

    Parameters
    ----------
    edges : (Mx2x2) array
        The coordinates of the sets of points that define the edges.

    Returns
    -------
    edge_angles : (Mx1) array
        The angles between the edges and the x-axis.
    """
    edges_count = len(edges)
    edge_angles = np.zeros(edges_count)
    for i in range(edges_count):
        edge_x = edges[i][1][0] - edges[i][0][0]
        edge_y = edges[i][1][1] - edges[i][0][1]
        edge_angles[i] = math.atan2(edge_y, edge_x)

    return np.unique(edge_angles) 

def compliance_function_fdiff(self, x, dc):
        obj = self.compliance_function(x, dc)

        x0 = x.copy()
        dc0 = dc.copy()
        dcf = np.zeros(dc.shape)
        for i, v in enumerate(x):
            x = x0.copy()
            x[i] += 1e-6
            o1 = self.compliance_function(x, dc)
            x[i] = x0[i] - 1e-6
            o2 = self.compliance_function(x, dc)
            dcf[i] = (o1 - o2) / (2e-6)
        print("finite differences: {:g}".format(np.linalg.norm(dcf - dc0)))
        dc[:] = dc0

        return obj 

def calculate_fdiff_stress(self, x, u, nu, side=1, dx=1e-6):
        """
        Calculate the derivative of the Von Mises stress using finite
        differences given the densities x, displacements u, and young modulus
        nu. Optionally, provide the side length (default: 1) and delta x
        (default: 1e-6).
        """
        ds = self.calculate_diff_stress(x, u, nu, side)
        dsf = numpy.zeros(x.shape)
        x = numpy.expand_dims(x, -1)
        for i in range(x.shape[0]):
            delta = scipy.sparse.coo_matrix(([dx], [[i], [0]]), shape=x.shape)
            s1 = self.calculate_stress((x + delta.A).squeeze(), u, nu, side)
            s2 = self.calculate_stress((x - delta.A).squeeze(), u, nu, side)
            dsf[i] = ((s1 - s2) / (2. * dx))[i]
        print("finite differences: {:g}".format(numpy.linalg.norm(dsf - ds)))
        return dsf 

def __init__(self, nelx, nely, volfrac, penal, rmin, ft, gui, bc):
        self.n = nelx * nely
        self.opt = nlopt.opt(nlopt.LD_MMA, self.n)
        self.passive = bc.get_passive_elements()
        self.xPhys = np.ones(self.n)
        if self.passive is not None:
            self.xPhys[self.passive] = 0

        # set bounds
        ub = np.ones(self.n, dtype=float)
        self.opt.set_upper_bounds(ub)
        lb = np.zeros(self.n, dtype=float)
        self.opt.set_lower_bounds(lb)

        # set stopping criteria
        self.opt.set_maxeval(2000)
        self.opt.set_ftol_rel(0.001)

        # set objective and constraint functions
        self.opt.set_min_objective(self.compliance_function)
        self.opt.add_inequality_constraint(self.volume_function, 0)

        # setup filter
        self.ft = ft
        self.filtering = Filter(nelx, nely, rmin)

        # setup problem def
        self.init_problem(nelx, nely, penal, bc)
        self.volfrac = volfrac

        # set GUI callback
        self.init_gui(gui) 

def compliance_function_fdiff(self, x, dc):
        obj = self.compliance_function(x, dc)

        x0 = x.copy()
        dc0 = dc.copy()
        dcf = np.zeros(dc.shape)
        for i, v in enumerate(x):
            x = x0.copy()
            x[i] += 1e-6
            o1 = self.compliance_function(x, dc)
            x[i] = x0[i] - 1e-6
            o2 = self.compliance_function(x, dc)
            dcf[i] = (o1 - o2) / (2e-6)
        print("finite differences: {:g}".format(np.linalg.norm(dcf - dc0)))
        dc[:] = dc0

        return obj 

def calculate_fdiff_stress(self, x, u, nu, side=1, dx=1e-6):
        """
        Calculate the derivative of the Von Mises stress using finite
        differences given the densities x, displacements u, and young modulus
        nu. Optionally, provide the side length (default: 1) and delta x
        (default: 1e-6).
        """
        ds = self.calculate_diff_stress(x, u, nu, side)
        dsf = numpy.zeros(x.shape)
        x = numpy.expand_dims(x, -1)
        for i in range(x.shape[0]):
            delta = scipy.sparse.coo_matrix(([dx], [[i], [0]]), shape=x.shape)
            s1 = self.calculate_stress((x + delta.A).squeeze(), u, nu, side)
            s2 = self.calculate_stress((x - delta.A).squeeze(), u, nu, side)
            dsf[i] = ((s1 - s2) / (2. * dx))[i]
        print("finite differences: {:g}".format(numpy.linalg.norm(dsf - ds)))
        return dsf 

def vectorizeDual(functions, order):
	n = len(functions)
	m = order + 1
	vector = np.zeros([n, m])
	for i, f in enumerate(functions):
   	 f.buildCoefficients(order)
   	 for j, coeff in enumerate(f.coefficients):
   		 vector[i,j] = coeff

	# Create a copy of one function
	f_everything = functions[0]
	f_everything.coefficients = np.zeros([n, m])
	for i in range(m):
   	 f_everything.coefficients[:, i] = vector[:, i]

	return f_everything 

def setup_buffer(width, height):
    """Set up the internal pixel buffer.

    :param width: width of buffer, ideally in multiples of 16
    :param height: height of buffer, ideally in multiples of 16

    """
    global _buffer_width, _buffer_height, _buf

    _buffer_width = width
    _buffer_height = height
    _buf = numpy.zeros((_buffer_width, _buffer_height, 3), dtype=int) 

def compute_mfcc(audio, **kwargs):
    """
    Compute the MFCC for a given audio waveform. This is
    identical to how DeepSpeech does it, but does it all in
    TensorFlow so that we can differentiate through it.
    """

    batch_size, size = audio.get_shape().as_list()
    audio = tf.cast(audio, tf.float32)

    # 1. Pre-emphasizer, a high-pass filter
    audio = tf.concat((audio[:, :1], audio[:, 1:] - 0.97*audio[:, :-1], np.zeros((batch_size,1000),dtype=np.float32)), 1)

    # 2. windowing into frames of 320 samples, overlapping
    windowed = tf.stack([audio[:, i:i+400] for i in range(0,size-320,160)],1)

    # 3. Take the FFT to convert to frequency space
    ffted = tf.spectral.rfft(windowed, [512])
    ffted = 1.0 / 512 * tf.square(tf.abs(ffted))

    # 4. Compute the Mel windowing of the FFT
    energy = tf.reduce_sum(ffted,axis=2)+1e-30
    filters = np.load("filterbanks.npy").T
    feat = tf.matmul(ffted, np.array([filters]*batch_size,dtype=np.float32))+1e-30

    # 5. Take the DCT again, because why not
    feat = tf.log(feat)
    feat = tf.spectral.dct(feat, type=2, norm='ortho')[:,:,:26]

    # 6. Amplify high frequencies for some reason
    _,nframes,ncoeff = feat.get_shape().as_list()
    n = np.arange(ncoeff)
    lift = 1 + (22/2.)*np.sin(np.pi*n/22)
    feat = lift*feat
    width = feat.get_shape().as_list()[1]

    # 7. And now stick the energy next to the features
    feat = tf.concat((tf.reshape(tf.log(energy),(-1,width,1)), feat[:, :, 1:]), axis=2)
    
    return feat 

def get_logits(new_input, length, first=[]):
    """
    Compute the logits for a given waveform.

    First, preprocess with the TF version of MFC above,
    and then call DeepSpeech on the features.
    """
    # new_input = tf.Print(new_input, [tf.shape(new_input)])

    # We need to init DeepSpeech the first time we're called
    if first == []:
        first.append(False)
        # Okay, so this is ugly again.
        # We just want it to not crash.
        tf.app.flags.FLAGS.alphabet_config_path = "DeepSpeech/data/alphabet.txt"
        DeepSpeech.initialize_globals()
        print('initialized deepspeech globals')

    batch_size = new_input.get_shape()[0]

    # 1. Compute the MFCCs for the input audio
    # (this is differentable with our implementation above)
    empty_context = np.zeros((batch_size, 9, 26), dtype=np.float32)
    new_input_to_mfcc = compute_mfcc(new_input)[:, ::2]
    features = tf.concat((empty_context, new_input_to_mfcc, empty_context), 1)

    # 2. We get to see 9 frames at a time to make our decision,
    # so concatenate them together.
    features = tf.reshape(features, [new_input.get_shape()[0], -1])
    features = tf.stack([features[:, i:i+19*26] for i in range(0,features.shape[1]-19*26+1,26)],1)
    features = tf.reshape(features, [batch_size, -1, 19*26])

    # 3. Whiten the data
    mean, var = tf.nn.moments(features, axes=[0,1,2])
    features = (features-mean)/(var**.5)

    # 4. Finally we process it with DeepSpeech
    logits = DeepSpeech.BiRNN(features, length, [0]*10)

    return logits 

def _split_into_chunks(signal, chunks_per_second=24):
    # TODO currently broken
    raise NotImplemented("Splitting to chunks is currently broken.")
    window_length_ms = 1/chunks_per_second * 1000
    intervals = np.arange(window_length_ms, signal.shape[0], window_length_ms, dtype=np.int32)
    chunks = np.array_split(signal, intervals, axis=0)
    pad_to = _next_power_of_two(np.max([chunk.shape[0] for chunk in chunks]))
    padded_chunks = np.stack(np.concatenate([chunk, np.zeros((pad_to - chunk.shape[0],))]) for chunk in chunks)
    return padded_chunks 

def display(self, image_generator):
        im = plt.imshow(np.zeros((1, 1, 3)), animated=True)
        def update(frame, *_):
            im.set_array(frame)
            return im,

        ani = FuncAnimation(self.fig, func=update, frames=image_generator, interval=self.pause_time, blit=True)
        plt.show() 

def random_start_img(img_size, batch_size, num_channels=3, num_ones_offset=None):
    zeros = np.zeros((*img_size, batch_size * num_channels))
    num_ones_offset = np.random.choice([1, 0], size=batch_size * num_channels)
    zeros += num_ones_offset * np.random.random(size=batch_size * num_channels)
    img = np.transpose(zeros.reshape((*img_size, batch_size, num_channels)), axes=[2, 0, 1, 3])
    return img 

def batch_norm(x, is_train_tensor=None):

    if is_train_tensor is None:
        is_training_tensor = tf.get_default_graph().get_tensor_by_name('is_training:0')

    normed, _, _ = tf.nn.fused_batch_norm(x, scale=tf.ones(shape=(x.get_shape()[-1],)),
                                          offset=tf.zeros(shape=(x.get_shape()[-1],)),
                                          is_training=is_train_tensor, name=None)
    return normed 

def bilinear_weight_initializer(filter_size, add_noise=True):

    def _initializer(shape, dtype=tf.float32, partition_info=None):
        if shape:
            # second last dimension is input, last dimension is output
            fan_in = float(shape[-2]) if len(shape) > 1 else float(shape[-1])
            fan_out = float(shape[-1])
        else:
            fan_in = 1.0
            fan_out = 1.0

        # define weight matrix (set dtype always to float32)
        weights = np.zeros((filter_size[0], filter_size[1], int(fan_in), int(fan_out)), dtype=dtype.as_numpy_dtype())

        # get bilinear kernel
        bilinear = bilinear_filt(filter_size=filter_size)
        bilinear = bilinear / fan_out  # normalize by number of channels

        # set filter in weight matrix (also allow channel mixing)
        for i in range(weights.shape[2]):
            for j in range(weights.shape[3]):
                weights[:, :, i, j] = bilinear

        # add small noise for symmetry breaking
        if add_noise:
            # define standard deviation so that it is equal to 1/2 of the smallest weight entry
            std = np.min(bilinear) / 2
            noise = truncated_normal(shape=shape, mean=0.0,
                                                              stddev=std, seed=None, dtype=dtype)
            weights += noise

        return weights

    return _initializer 

def _get_peak_coverage(self):
        norm_coverage_folder = "{}/normalized_coverage".format(
            self._project_folder)
        coverage_files = [join(norm_coverage_folder, f) for f in listdir(
            norm_coverage_folder) if isfile(join(norm_coverage_folder, f))]
        wiggle_parser = WiggleParser()
        cons_value_dict = defaultdict(dict)
        for coverage_file in coverage_files:
            cons_values = np.zeros(self._consensus_length)
            with open(coverage_file, 'r') as cov_fh:
                for wiggle_entry in wiggle_parser.entries(cov_fh):
                    lib_name_and_strand = wiggle_entry.track_name
                    lib_name = '_'.join(lib_name_and_strand.split('_')[:-1])
                    lib_strand = '+' if lib_name_and_strand.split(
                        '_')[-1] == "forward" else '-'
                    replicon = wiggle_entry.replicon
                    pos_value_pairs = dict(wiggle_entry.pos_value_pairs)
                    self._get_coverage_for_replicon_peaks(
                        replicon, lib_strand, pos_value_pairs, cons_values)
            cons_value_dict[lib_name][lib_strand] = cons_values
        # combine strands
        comb_cons_value_dict = {}
        for lib in cons_value_dict:
            comb_cons_value_dict[lib] = np.zeros(self._consensus_length)
            for strand in cons_value_dict[lib]:
                comb_cons_value_dict[lib] += cons_value_dict[lib][strand]
        return comb_cons_value_dict 

def comenzarPartida():

    global juegoComenzado
    global jugador1
    global jugador2
    global texto
    global canvas
    global matrix
    global fin
    global usados
    
    if juegoComenzado == False:
        jugador1 = True

    #Reiniciamos la partida
    elif juegoComenzado == True:
        
        jugador1 = True
        jugador2 = False
        canvas.delete("all")
        #Canvas donde se dibuja todo
        canvas.create_line(0,100,400,100) #Lineas que dibujan el tablero 4x4
        canvas.create_line(0,200,400,200)
        canvas.create_line(0,300,400,300)
        canvas.create_line(0,400,400,400)
        canvas.create_line(100,0,100,400)
        canvas.create_line(200,0,200,400)
        canvas.create_line(300,0,300,400)
        canvas.create_line(400,0,400,400)
        matrix = numpy.zeros((4,4)) 
        fin = False
        usados=0
        
    texto.set("Estado del juego: Turno jugador 1.")
    juegoComenzado = True 

def inverse_3by3_double(M):
    """C style inverse of a flattened 3 by 3 matrix, returning full matrix in
    floating-point, also in flattened format.

    :param M: The matrix to find inverse to.
    :type M: :py:class:`numpy.ndarray`
    :return: The inverse matrix.
    :rtype: :py:class:`numpy.ndarray`

    """
    if len(M.shape) > 1:
        M = M.flatten()

    M = np.array(M, 'float')

    determinant = 0.
    adj_M = np.zeros((9,), 'float')

    # First row of adjugate matrix
    adj_M[0] = (M[4] * M[8] - M[7] * M[5])  # Det #0
    adj_M[1] = -(M[1] * M[8] - M[7] * M[2])
    adj_M[2] = (M[1] * M[5] - M[4] * M[2])

    # Second row of adjugate matrix
    adj_M[3] = -(M[3] * M[8] - M[6] * M[5])  # Det #1
    adj_M[4] = (M[0] * M[8] - M[6] * M[2])
    adj_M[5] = -(M[0] * M[5] - M[3] * M[2])

    # Third row of adjugate matrix
    adj_M[6] = (M[3] * M[7] - M[6] * M[4])  # Det #2
    adj_M[7] = -(M[0] * M[7] - M[6] * M[1])
    adj_M[8] = (M[0] * M[4] - M[3] * M[1])

    determinant += M[0] * adj_M[0]
    determinant += M[1] * adj_M[3]  # Using addition since minus is integrated in adjugate matrix.
    determinant += M[2] * adj_M[6]

    return (adj_M / determinant) 

def QR_algorithm_shift_Givens_double(A):
    """The QR algorithm with largest value shift for finding eigenvalues.

    Using Givens rotations for finding RQ.

    :param A: The square matrix to find eigenvalues of.
    :type A: :py:class:`numpy.ndarray`
    :return: The eigenvalues.
    :rtype: list

    """

    # First, tridiagonalize the input matrix to a Hessenberg matrix.
    T = ma.convert_to_Hessenberg_Givens_double(A.copy())

    m, n = T.shape
    if n != m:
        raise np.linalg.LinAlgError("Array must be square.")

    convergence_measure = []
    eigenvalues = np.zeros((m, ), dtype='float')
    m -= 1

    while m > 0:
        # Obtain the shift from the lower right corner of the matrix.
        mu_matrix = np.eye(T.shape[0]) * T[m, m]
        # Perform Givens QR step (which returns RQ) on the shifted
        # matrix and then shift it back.
        T = Givens_QR_step_double(T - mu_matrix) + mu_matrix
        # Add convergence information and extract eigenvalue if close enough.
        convergence_measure.append(np.abs(T[m, m - 1]))
        if convergence_measure[-1] < QR_ALGORITHM_TOLERANCE:
            eigenvalues[m] = T[m, m]
            T = T[:m, :m]
            m -= 1

    eigenvalues[0] = T
    return eigenvalues, convergence_measure 

def QR_algorithm_shift_Givens_int(A):
    """The QR algorithm with largest value shift for finding eigenvalues, in integer precision.

    Using Givens rotations for finding RQ.

    :param A: The square matrix to find eigenvalues of.
    :type A: :py:class:`numpy.ndarray`
    :return: The eigenvalues.
    :rtype: list

    """

    A, scale = fp.scale_64bit_matrix(A.copy())

    # First, tridiagonalize the input matrix to a Hessenberg matrix.
    T = ma.convert_to_Hessenberg_Givens_int(A)

    m, n = T.shape
    if n != m:
        raise np.linalg.LinAlgError("Array must be square.")

    convergence_measure = []
    eigenvalues = np.zeros((m, ), dtype='int64')
    m -= 1

    while m > 0:
        # Obtain the shift from the lower right corner of the matrix.
        mu_matrix = np.eye(T.shape[0], dtype='int64') * T[m, m]
        # Perform Givens QR step (which returns RQ) on the shifted
        # matrix and then shift it back.
        T = Givens_QR_step_int(T - mu_matrix) + mu_matrix
        # Add convergence information and extract eigenvalue if close enough.
        convergence_measure.append(np.abs(T[m, m - 1]))
        if convergence_measure[-1] <= 1:
            eigenvalues[m] = T[m, m]
            T = T[:m, :m]
            m -= 1

    eigenvalues[0] = T
    return eigenvalues << scale, convergence_measure 

def overlap(x,y):
    # You can't overlap sdrs
    # with different shapes
    if x.shape != y.shape:
        return None
    result = np.zeros(x.shape,dtype=np.bool_)
    for i in range(x.shape[0]):
        result[i] = x[i] and y[i]
    return result 

def run(self,input_sdr):
        self.current = 0
        output = np.zeros(len(self.collumns),dtype=np.bool_)
        for i in range(len(self.collumns)):
            output[i] = self.spatial_pooler(input_sdr)
            self.current += 1
        self.current = 0
        return output 

def stsb_label_encoding(labels, nclass=6):
    """
    Label encoding from Tree LSTM paper (Tai, Socher, Manning)
    """
    Y = np.zeros((len(labels), nclass)).astype(np.float32)
    for j, y in enumerate(labels):
        for i in range(nclass):
            if i == np.floor(y) + 1:
                Y[j, i] = y - np.floor(y)
            if i == np.floor(y):
                Y[j, i] = np.floor(y) - y + 1
    return Y 

def transform_roc(X1, X2, X3):
    n_batch = len(X1)
    xmb = np.zeros((n_batch, 2, n_ctx, 2), dtype=np.int32)
    mmb = np.zeros((n_batch, 2, n_ctx), dtype=np.float32)
    start = encoder['_start_']
    delimiter = encoder['_delimiter_']
    for i, (x1, x2, x3), in enumerate(zip(X1, X2, X3)):
        x12 = [start] + x1[:max_len] + [delimiter] + x2[:max_len] + [clf_token]
        x13 = [start] + x1[:max_len] + [delimiter] + x3[:max_len] + [clf_token]
        l12 = len(x12)
        l13 = len(x13)
        xmb[i, 0, :l12, 0] = x12
        xmb[i, 1, :l13, 0] = x13
        mmb[i, 0, :l12] = 1
        mmb[i, 1, :l13] = 1
    xmb[:, :, :, 1] = np.arange(
        n_vocab + n_special, n_vocab + n_special + n_ctx)
    return xmb, mmb 

def transform_sst(X1):
    n_batch = len(X1)
    xmb = np.zeros((n_batch, 1, n_ctx, 2), dtype=np.int32)
    mmb = np.zeros((n_batch, 1, n_ctx), dtype=np.float32)
    start = encoder['_start_']
    delimiter = encoder['_delimiter_']
    for i, x1, in enumerate(X1):
        x1 = [start] + x1[:max_len] + [clf_token]
        l1 = len(x1)
        xmb[i, 0, :l1, 0] = x1
        mmb[i, 0, :l1] = 1
    xmb[:, :, :, 1] = np.arange(
        n_vocab + n_special, n_vocab + n_special + n_ctx)
    return xmb, mmb 

def __init__(self, x):
        """
        Args:
            x (cvxpy.Variable): a symbolic representation of
                members of the set
        """
        constr = [x == np.zeros(np.prod(x.size))]
        self._unqueried = True
        super(Zeros, self).__init__(x, constr) 

def project(self, x_0):
        return x_0 if self.contains(x_0) else np.zeros(x_0.shape) 

def project(self, x_0):
        if self.contains(x_0):
            return x_0
        z = x_0[:-1]
        t = x_0[-1]
        norm_z = np.linalg.norm(z, 2)

        # As given in [BV04, Chapter 8.Exercises]
        if norm_z <= -t:
            # this certainly will never happen when t > 0
            return np.zeros(np.shape(x_0))
        elif self._contains(norm_z, t):
            return x_0
        else:
            return 0.5 * (1 + t/norm_z) * np.append(z, norm_z) 

def project(self, x_0):
        assert self._shape == x_0.shape, \
            'cone shape (%s) != x_0 shape (%s)' % (
            str(self._shape), str(x_0.shape))
        x_star = np.zeros(x_0.shape)
        for s, slx in zip(self.sets, self.slices):
            x_star[slx] = s.project(x_0[slx])
        return x_star 

def test_projection(self):
        """Test projection, query methods of the AffineSet oracle."""
        m = 150
        n = 200
        x = cvxpy.Variable(n)

        sol = np.random.randn(n)
        A = np.random.randn(m, n)
        b = A.dot(sol)

        affine = AffineSet(x, A, b)
        constr = affine._constr

        # Test idempotency
        x_0 = sol
        x_star = affine.project(x_0)
        self.assertTrue(affine.contains(x_star))
        self.assertTrue(np.array_equal(x_0, x_star), "projection not "
            "idemptotent")
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.pnorm(x - x_0, 2)), constr)
        p.solve()
        self.assertTrue(np.isclose(np.array(x.value).flatten(), x_star,
            atol=1e-3).all())
        utils.query_helper(self, x_0, x_star, affine, idempotent=True)
        
        # Test random x_0
        x_0 = np.random.randn(n)
        while (np.isclose(A.dot(x_0), np.zeros(m)).all()):
            x_0 = np.random.randn(n)
        x_star = affine.project(x_0)
        self.assertTrue(np.isclose(A.dot(x_star), b).all())
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.pnorm(x - x_0, 2)), constr)
        p.solve()
        self.assertTrue(np.isclose(np.array(x.value).flatten(), x_star,
            atol=1e-3).all())
        utils.query_helper(self, x_0, x_star, affine, idempotent=False) 

def build_text_and_labels(texts, class_labels, lowercase = True, multilingual = False, langs = None, distinct_labels_index = None):
	# Load data from files
	lines = [(text.lower() if lowercase else text).strip().split() for text in texts]
	x_instances = [l[1:-1] for l in lines] if multilingual else [l[:-1] for l in lines]
	
	dist_labels = list(set(class_labels)) if distinct_labels_index is None else distinct_labels_index
	y_instances = [np.zeros(len(dist_labels)) for l in class_labels]
	for i in range(len(y_instances)):
		y_instances[i][dist_labels.index(class_labels[i])] = 1

	return [x_instances, y_instances, langs, dist_labels] if multilingual else [x_instances, y_instances, dist_labels] 

def aggregate_phrase_embedding(words, stopwords, embs, emb_size, l2_norm_vec = True, lang = 'en'):
	vec_res = np.zeros(emb_size)
	fit_words = [w.lower() for w in words if w.lower() not in stopwords and w.lower() in embs.lang_vocabularies[lang]]
	if len(fit_words) == 0:
		return None

	for w in fit_words:
		vec_res += embs.get_vector(lang, w) 	
	res = np.multiply(1.0 / (float(len(fit_words))), vec_res)
	if l2_norm_vec:
		res = np.multiply(1.0 / np.linalg.norm(res), res)
	return res 

def __init__(self, labels = [], predictions = [], gold = [], one_hot_encoding = False, class_indices = False):
		# rows are true labels, columns predictions
		self.matrix = np.zeros(shape = (len(labels), len(labels)))
		self.labels = labels

		if len(predictions) != len(gold):
			raise ValueError("Predictions and gold labels do not have the same count.")
		for i in range(len(predictions)):
			index_pred = np.argmax(predictions[i]) if one_hot_encoding else (predictions[i] if class_indices else labels.index(predictions[i]))
			index_gold = np.argmax(gold[i]) if one_hot_encoding else (gold[i] if class_indices else labels.index(gold[i]))
			self.matrix[index_gold][index_pred] += 1

		if len(predictions) > 0: 
			self.compute_all_scores() 

def covariance_matrix(first, second):
	if first.shape[0] != second.shape[0]:
		raise ValueError("Input matrices must have the same number of row vectors (same first dimensions)")
	
	cov_mat = np.zeros((first.shape[1], second.shape[1]))
	for i in range(first.shape[0]):
		cov_mat = np.add(cov_mat, np.matmul(np.transpose([first[i]]), np.array([second[i]])))
	cov_mat = np.multiply(1.0 / float(first.shape[0]), cov_mat)
	return cov_mat 

def edge_property_as_matrix(g, prop_name):
    """edge_property_as_matrix
    Returns an edge property as a matrix. If the graph is undirected, a
    symmetric graph is returned.
    
    Parameters
    ----------
        g : :obj:`graph_tool.Graph` 
            A graph-tool type graph.
        prop_name : :obj:`str`
            The name of an edge property belonging to Graph, g.

    Returns
    -------
     mat : :obj:`np.array` 
        The edge property in matrix form, with each
    """
    n_vertices = len([_ for _ in g.vertices()])
    mat = np.zeros((n_vertices, n_vertices))
    prop = g.ep[prop_name]
    edges = sorted(g.edges())
    if g.is_directed():
        for e in edges:
            s = int(e.source())
            t = int(e.target())            
            mat[s][t] = prop[e]
    else:
        for e in edges:
            s = int(e.source())
            t = int(e.target())            
            mat[s][t] = prop[e]
            mat[t][s] = prop[e]
    return mat 

def specification(self, specification):
        if isinstance(specification, (int)):
            if np.abs(specification) > self._triangsamples.vertlist.shape[0]:
                raise ValueError("""The Number of selected basic functions is
                too large.""")
            else:
                if specification == 0:
                    self._specification = \
                        np.ones(self._triangsamples.vertlist.shape[0])
                else:
                    self._specification = \
                        np.zeros(self._triangsamples.vertlist.shape[0])
                    if specification > 0:
                        self._specification[:specification] = 1
                    else:
                        self._specification[specification:] = 1
        elif isinstance(specification, (list, tuple, np.ndarray)):
            specification = np.asarray(specification)
            if specification.shape[0] != self._triangsamples.vertlist.shape[1]:
                raise IndexError("""The length of the specification vector
                does not match the number of spatial sample points. """)
            else:
                self._specification = specification
        else:
            raise TypeError("""The parameter specification has to be
            int or a vecor""") 

def from_q_to_rad_axis(q):
        """
        :param q: quaternion in w-i-j-k format
        :return: rad in format x-y-z
        """
        norm_axis = math.sqrt(q[1] ** 2 + q[2] ** 2 + q[3] ** 2)
        rad = 2.0 * math.atan2(norm_axis, q[0])
        if norm_axis < 1e-5:
            return np.zeros((3,), dtype=np.float32)
        else:
            return rad * np.asarray(q[1:], dtype=np.float32) / norm_axis 

def fhs(self, n_scenarios=250, start_date=None, end_date=None):
        x = sliding_window(n_scenarios+1, range(len(self.ts.index)))
        scenarios = np.zeros((len(self.ts.index), n_scenarios+1))
        for i, el in enumerate(x):
            l = list(el)
            cur_idx, hist_idx = l[-1], l[:-1]
            neutral = self.ts.Value.values[cur_idx]
            ret = self.ts.DevolLogReturns.values[hist_idx]
            vol = self.ts.Vola.values[cur_idx]
            scenarios[cur_idx, 1:] = self.scenario_values(ret, neutral, vol)
            scenarios[cur_idx, 0] = neutral
        return scenarios 

def get_forces(self):
        # Return the force vector for the problem
        ndof = 2 * (self.nelx + 1) * (self.nely + 1)
        f = np.zeros((ndof, 1))
        f[1, 0] = -1
        return f 

def __init__(self, nelx, nely, rmin):
        """
        Filter: Build (and assemble) the index+data vectors for the coo matrix
        format.
        """
        nfilter = int(nelx * nely * ((2 * (np.ceil(rmin) - 1) + 1)**2))
        iH = np.zeros(nfilter)
        jH = np.zeros(nfilter)
        sH = np.zeros(nfilter)
        cc = 0
        for i in range(nelx):
            for j in range(nely):
                row = i * nely + j
                kk1 = int(np.maximum(i - (np.ceil(rmin) - 1), 0))
                kk2 = int(np.minimum(i + np.ceil(rmin), nelx))
                ll1 = int(np.maximum(j - (np.ceil(rmin) - 1), 0))
                ll2 = int(np.minimum(j + np.ceil(rmin), nely))
                for k in range(kk1, kk2):
                    for l in range(ll1, ll2):
                        col = k * nely + l
                        fac = rmin - np.sqrt(
                            ((i - k) * (i - k) + (j - l) * (j - l)))
                        iH[cc] = row
                        jH[cc] = col
                        sH[cc] = np.maximum(0.0, fac)
                        cc = cc + 1
        # Finalize assembly and convert to csc format
        self.H = scipy.sparse.coo_matrix((sH, (iH, jH)),
            shape=(nelx * nely, nelx * nely)).tocsc()
        self.Hs = self.H.sum(1) 

def __init__(self, nelx, nely, volfrac, penal, rmin, ft, gui, bc):
        self.n = nelx * nely
        self.opt = nlopt.opt(nlopt.LD_MMA, self.n)
        self.passive = bc.get_passive_elements()
        self.xPhys = np.ones(self.n)
        if self.passive is not None:
            self.xPhys[self.passive] = 0

        # set bounds
        ub = np.ones(self.n, dtype=float)
        self.opt.set_upper_bounds(ub)
        lb = np.zeros(self.n, dtype=float)
        self.opt.set_lower_bounds(lb)

        # set stopping criteria
        self.opt.set_maxeval(2000)
        self.opt.set_ftol_rel(0.001)

        # set objective and constraint functions
        self.opt.set_min_objective(self.compliance_function)
        self.opt.add_inequality_constraint(self.volume_function, 0)

        # setup filter
        self.ft = ft
        self.filtering = Filter(nelx, nely, rmin)

        # setup problem def
        self.init_problem(nelx, nely, penal, bc)
        self.volfrac = volfrac

        # set GUI callback
        self.init_gui(gui) 

def build_indices(self, nelx, nely):
        """ FE: Build the index vectors for the for coo matrix format. """
        self.KE = self.lk()
        self.edofMat = np.zeros((nelx * nely, 8), dtype=int)
        for elx in range(nelx):
            for ely in range(nely):
                el = ely + elx * nely
                n1 = (nely + 1) * elx + ely
                n2 = (nely + 1) * (elx + 1) + ely
                self.edofMat[el, :] = np.array([2 * n1 + 2, 2 * n1 + 3,
                    2 * n2 + 2, 2 * n2 + 3, 2 * n2, 2 * n2 + 1, 2 * n1,
                    2 * n1 + 1])
        # Construct the index pointers for the coo format
        self.iK = np.kron(self.edofMat, np.ones((8, 1))).flatten()
        self.jK = np.kron(self.edofMat, np.ones((1, 8))).flatten() 

def __init__(self, nelx, nely, title=""):
        """Initialize plot and plot the initial design"""
        plt.ion()  # Ensure that redrawing is possible
        self.fig, self.ax = plt.subplots()
        self.im = self.ax.imshow(-np.zeros((nelx, nely)).T, cmap='gray',
            interpolation='none', norm=colors.Normalize(vmin=-1, vmax=0))
        plt.xlabel(title)
        # self.fig.tight_layout()
        self.fig.show()
        self.nelx, self.nely = nelx, nely 

def build_indices(self):
        edofMat = numpy.zeros((8, self.nelx * self.nely), dtype=int)
        for elx in range(self.nelx):
            for ely in range(self.nely):
                el = ely + elx * self.nely # Element index
                n1 = (self.nely + 1) * elx + ely # Left nodes
                n2 = (self.nely + 1) * (elx + 1) + ely # Right nodes
                edofMat[:, el] = numpy.array([2 * n1 + 2, 2 * n1 + 3,
                    2 * n2 + 2, 2 * n2 + 3, 2 * n2, 2 * n2 + 1, 2 * n1,
                    2 * n1 + 1])
        return edofMat 

def get_forces(self):
        # Return the force vector for the problem
        topx_to_id = np.vectorize(
            lambda x: xy_to_id(x, 0, self.nelx, self.nely))
        topx = 2 * topx_to_id(np.arange((self.nelx + 1) // 2)) + 1
        f = np.zeros((2 * (self.nelx + 1) * (self.nely + 1), 1))
        f[topx, 0] = -100
        return f 

def get_forces(self):
        # Return the force vector for the problem
        ndof = 2 * (self.nelx + 1) * (self.nely + 1)
        f = np.zeros((ndof, 1))
        f[1, 0] = -1
        return f 

def build_indices(self, nelx, nely):
        """ FE: Build the index vectors for the for coo matrix format. """
        self.KE = self.lk()
        self.edofMat = np.zeros((nelx * nely, 8), dtype=int)
        for elx in range(nelx):
            for ely in range(nely):
                el = ely + elx * nely
                n1 = (nely + 1) * elx + ely
                n2 = (nely + 1) * (elx + 1) + ely
                self.edofMat[el, :] = np.array([2 * n1 + 2, 2 * n1 + 3,
                    2 * n2 + 2, 2 * n2 + 3, 2 * n2, 2 * n2 + 1, 2 * n1,
                    2 * n1 + 1])
        # Construct the index pointers for the coo format
        self.iK = np.kron(self.edofMat, np.ones((8, 1))).flatten()
        self.jK = np.kron(self.edofMat, np.ones((1, 8))).flatten() 

def __init__(self, nelx, nely, penal, bc):
        # Problem size
        self.nelx = nelx
        self.nely = nely

        # Max and min stiffness
        self.Emin = 1e-9
        self.Emax = 1.0

        # SIMP penalty
        self.penal = penal

        # dofs:
        self.ndof = 2 * (nelx + 1) * (nely + 1)

        # FE: Build the index vectors for the for coo matrix format.
        self.build_indices(nelx, nely)

        # BC's and support (half MBB-beam)
        dofs = np.arange(2 * (nelx + 1) * (nely + 1))
        self.fixed = bc.get_fixed_nodes()
        self.free = np.setdiff1d(dofs, self.fixed)

        # Solution and RHS vectors
        self.f = bc.get_forces()
        self.u = np.zeros(self.f.shape)

        # Per element compliance
        self.ce = np.zeros(nely * nelx) 

def __init__(self, nelx, nely, title=""):
        """Initialize plot and plot the initial design"""
        plt.ion()  # Ensure that redrawing is possible
        self.fig, self.ax = plt.subplots()
        self.im = self.ax.imshow(-np.zeros((nelx, nely)).T, cmap='gray',
            interpolation='none', norm=colors.Normalize(vmin=-1, vmax=0))
        plt.xlabel(title)
        # self.fig.tight_layout()
        self.fig.show()
        self.nelx, self.nely = nelx, nely 

def get_forces(self):
        # Return the force vector for the problem
        topx_to_id = np.vectorize(
            lambda x: xy_to_id(x, 0, self.nelx, self.nely))
        topx = 2 * topx_to_id(np.arange((self.nelx + 1) // 2)) + 1
        f = np.zeros((2 * (self.nelx + 1) * (self.nely + 1), 1))
        f[topx, 0] = -100
        return f 

def get_forces(self):
        """ Return the force vector for the problem. """
        ndof = 2 * (self.nelx + 1) * (self.nely + 1)
        f = np.zeros((ndof, 1))
        fx = self.nelx
        # fy = (self.nely - self.passive_max_y) // 2 + self.passive_max_y
        for i in range(1, 2):
            fy = self.passive_max_y - 1 + 2 * i
            id = xy_to_id(fx, fy, self.nelx, self.nely)
            f[2 * id + 1, 0] = -1
        return f 

def mesh_from_img(img):
    nv = (img.shape[0] + 1) * (img.shape[1] + 1)
    nf = img.size * 2

    v_count = 0
    f_count = 0
    V_dict = {}

    V = np.zeros([nv, 2])
    F = np.zeros([nf, 3], dtype=np.int)

    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            val = img[i, j]
            if val == 255.0:
                continue

            v_idx = []
            for v_i in [(i, j), (i + 1, j), (i, j + 1), (i + 1, j + 1)]:
                if v_i in V_dict:
                    v_idx.append(V_dict[v_i])
                else:
                    V_dict[v_i] = v_count
                    V[v_count, :] = np.array((v_i[1], -v_i[0]))
                    v_count += 1
                    v_idx.append(v_count - 1)

            v1, v2, v3, v4 = v_idx

            F[f_count, :] = np.array([v1, v2, v4])
            F[f_count + 1, :] = np.array([v1, v4, v3])
            f_count += 2

    V = np.resize(V, [v_count, 2])
    F = np.resize(F, [f_count, 3])

    return V, F 

def _calcJacobian(self, k, n, jacobian_value):
		if np.all(np.equal(jacobian_value, None)):
			if k != 0:
				rows = self.val.shape[0]
				seed = np.zeros([rows, n])
				seed[:, k-1] = 1
				return seed
		else:
			return jacobian_value 

def to_phalf_from_pfull(arr, val_toa=0, val_sfc=0):
    """Compute data at half pressure levels from values at full levels.

    Could be the pressure array itself, but it could also be any other data
    defined at pressure levels.  Requires specification of values at surface
    and top of atmosphere.
    """
    phalf = np.zeros((arr.shape[0] + 1, arr.shape[1], arr.shape[2]))
    phalf[0] = val_toa
    phalf[-1] = val_sfc
    phalf[1:-1] = 0.5*(arr[:-1] + arr[1:])
    return phalf 

def inverse_symmetric_3by3_double(M):
    """C style inverse of a symmetric, flattened 3 by 3 matrix.

    Returns the adjoint matrix and the determinant, s.t.

    .. math::

        M^{-1} = \\frac{1}{\\det(M)} \\cdot \\text{adj}(M).


    For integer matrices, this then returns exact results by avoiding
    to apply the division.

    :param M: The matrix to find inverse to. Assumes array with shape (6,).
    :type M: :py:class:`numpy.ndarray`
    :return: The inverse matrix, flattened.
    :rtype: :py:class:`numpy.ndarray`

    """

    determinant = 0
    inverse = np.zeros((9,), dtype='float')

    # First row of adjunct matrix
    inverse[0] = (M[3] * M[5] - (M[4] ** 2))  # Det #0
    inverse[1] = -(M[1] * M[5] - M[4] * M[2])  # Det #1
    inverse[2] = (M[1] * M[4] - M[3] * M[2])  # Det #2

    # Second row of adjunct matrix
    inverse[3] = inverse[1]
    inverse[4] = (M[0] * M[5] - (M[2] ** 2))
    inverse[5] = -(M[0] * M[4] - M[1] * M[2])

    # Third row of adjunct matrix
    inverse[6] = inverse[2]
    inverse[7] = inverse[5]
    inverse[8] = (M[0] * M[3] - (M[1] ** 2))

    determinant += M[0] * inverse[0]
    determinant += M[1] * inverse[1]  # Using addition since minus is integrated in adjunct matrix.
    determinant += M[2] * inverse[2]

    for i in xrange(len(inverse)):
        inverse[i] /= determinant

    return inverse 

def inverse_symmetric_3by3_int(M):
    """C style inverse of a symmetric, flattened 3 by 3 matrix.

    Returns the adjoint matrix and the determinant, s.t.

    .. math::

        M^{-1} = \\frac{1}{\\det(M)} \\cdot \\text{adj}(M).


    For integer matrices, this then returns exact results by avoiding
    to apply the division.

    :param M: The matrix to find inverse to. Assumes array with shape (6,).
    :type M: :py:class:`numpy.ndarray`
    :return: The adjoint flattened matrix and the determinant.
    :rtype: tuple

    """

    determinant = 0
    adj_M = np.zeros((9,), dtype='int32')

    # First row of adjunct matrix
    adj_M[0] = (M[3] * M[5] - (M[4] ** 2))  # Det #0
    adj_M[1] = -(M[1] * M[5] - M[4] * M[2])  # Det #1
    adj_M[2] = (M[1] * M[4] - M[3] * M[2])  # Det #2

    # Second row of adjunct matrix
    adj_M[3] = adj_M[1]
    adj_M[4] = (M[0] * M[5] - (M[2] ** 2))
    adj_M[5] = -(M[0] * M[4] - M[1] * M[2])

    # Third row of adjunct matrix
    adj_M[6] = adj_M[2]
    adj_M[7] = adj_M[5]
    adj_M[8] = (M[0] * M[3] - (M[1] ** 2))

    determinant += np.int64(M[0]) * np.int64(adj_M[0])
    determinant += np.int64(M[1]) * np.int64(adj_M[1])  # Using addition since minus is integrated in adjunct matrix.
    determinant += np.int64(M[2]) * np.int64(adj_M[2])

    return adj_M, determinant 

def inverse_3by3_int(M):
    """C style inverse of a flattened 3 by 3 matrix.

    Returns the adjoint matrix and the determinant, s.t.

    .. math::

        M^{-1} = \\frac{1}{\\det(M)} \\cdot \\text{adj}(M).


    For integer matrices, this then returns exact results by avoiding
    to apply the division.

    :param M: The matrix to find inverse to.
    :type M: :py:class:`numpy.ndarray`
    :return: The adjoint flattened matrix and the determinant.
    :rtype: tuple

    """
    if len(M.shape) > 1:
        M = M.flatten()

    determinant = 0
    adj_M = np.zeros((9,), 'int')

    # First row of adjunct matrix
    adj_M[0] = (M[4] * M[8] - M[7] * M[5])  # Det #0
    adj_M[1] = -(M[1] * M[8] - M[7] * M[2])
    adj_M[2] = (M[1] * M[5] - M[4] * M[2])

    # Second row of adjunct matrix
    adj_M[3] = -(M[3] * M[8] - M[6] * M[5])  # Det #1
    adj_M[4] = (M[0] * M[8] - M[6] * M[2])
    adj_M[5] = -(M[0] * M[5] - M[3] * M[2])

    # Third row of adjunct matrix
    adj_M[6] = (M[3] * M[7] - M[6] * M[4])  # Det #2
    adj_M[7] = -(M[0] * M[7] - M[6] * M[1])
    adj_M[8] = (M[0] * M[4] - M[3] * M[1])

    determinant += np.int64(M[0]) * np.int64(adj_M[0])
    determinant += np.int64(M[1]) * np.int64(adj_M[3])  # Using addition since minus is integrated in adjunct matrix.
    determinant += np.int64(M[2]) * np.int64(adj_M[6])

    return adj_M, determinant 

def inverse_3by3_double(M):
    """C style inverse of a flattened 3 by 3 matrix.

    .. math::

        M^{-1} = \\frac{1}{\\det(M)} \\cdot \\text{adj}(M).


    For integer matrices, this then returns exact results by avoiding
    to apply the division.

    :param M: The matrix to find inverse to.
    :type M: :py:class:`numpy.ndarray`
    :return: The inverse matrix.
    :rtype: :py:class:`numpy.ndarray`

    """
    if len(M.shape) > 1:
        M = M.flatten()

    determinant = 0
    adj_M = np.zeros((9,), 'int')

    # First row of adjunct matrix
    adj_M[0] = (M[4] * M[8] - M[7] * M[5])  # Det #0
    adj_M[1] = -(M[1] * M[8] - M[7] * M[2])
    adj_M[2] = (M[1] * M[5] - M[4] * M[2])

    # Second row of adjunct matrix
    adj_M[3] = -(M[3] * M[8] - M[6] * M[5])  # Det #1
    adj_M[4] = (M[0] * M[8] - M[6] * M[2])
    adj_M[5] = -(M[0] * M[5] - M[3] * M[2])

    # Third row of adjunct matrix
    adj_M[6] = (M[3] * M[7] - M[6] * M[4])  # Det #2
    adj_M[7] = -(M[0] * M[7] - M[6] * M[1])
    adj_M[8] = (M[0] * M[4] - M[3] * M[1])

    determinant += M[0] * adj_M[0]
    determinant += M[1] * adj_M[3]  # Using addition since minus is integrated in adjunct matrix.
    determinant += M[2] * adj_M[6]

    return adj_M / determinant 

def convert_to_Hessenberg_symmetric_double(A):
    """Tridiagonalize a symmetric square matric to Hessenberg form using Householder reflections.

    Costs 4/3 * n^3 + O(n^2) operations.

    .. code-block:: matlab

        function A = mytridiagturbo(A)

        [m,n] = size(A);
        if (m ~= n)
            error('Input matrix is not square.')
        end

        for k = 1:(m - 2)
            vk = A((k+1):m,k);
            beta = sign(vk(1))*norm(vk);
            vk(1) = vk(1) + beta;
            vk = vk / norm(vk);

            % Calculate A_hat for current iteration.
            u = -2*(A((k+1:m),(k+1:m)) * vk);
            d = -vk' * u;
            w = u + d * vk;
            wv = w * vk';
            A((k+1:m),(k+1:m)) = A((k+1:m),(k+1:m)) + wv + wv';

            % Adjust tridiagonality.
            A((k+1):m,k) = [-beta; zeros(m-k-1,1)];
            A(k,(k+1):m) = A((k+1):m,k)';
         end

    :param A: The matrix to convert.
    :type A: :py:class:`numpy.ndarray`
    :return: The Hessenberg matrix.
    :rtype: :py:class:`numpy.ndarray`

    """

    m, n = A.shape
    A = np.array(A, 'float')

    for k in xrange(m - 1):
        vk = A[(k + 1):, k].copy()
        beta = np.sign(vk[0]) * np.linalg.norm(vk)
        vk[0] += beta
        vk = vk / np.linalg.norm(vk)

        # Calculate A_hat for current iteration.
        u = -2 * np.dot(A[(k + 1):, (k + 1):], vk)
        d = np.dot(-vk, u)
        w = u + d * vk
        wv = np.outer(w, vk)
        A[(k + 1):, (k + 1):] += wv + wv.T

        # Adjust tridiagonality.
        t = np.zeros((m - (k + 1), ), dtype='float')
        t[0] = -beta
        A[(k + 1):, k] = t
        A[k, (k + 1):] = t
    return A 

def inverse_symmetric_3by3_int64(M):
    """C style inverse of a symmetric, flattened 3 by 3 matrix, represented
    by the six unique values

    .. code-block:: python

        np.array([S[0, 0], S[0, 1], S[0, 2], S[1, 1], S[1, 2], S[2, 2]])

    Returns the flattened full adjugate matrix and the determinant, s.t.

    .. math::

        M^{-1} = \\frac{1}{\\det(M)} \\cdot \\text{adj}(M).


    For integer matrices, this then returns exact results by avoiding
    to apply the division.

    :param M: The matrix to find inverse to. Assumes array with shape (6,).
    :type M: :py:class:`numpy.ndarray`
    :return: The adjugate flattened matrix and the determinant.
    :rtype: tuple

    """

    determinant = 0
    adj_M = np.zeros((9,), dtype='int64')

    # First row of adjugate matrix
    adj_M[0] = (M[3] * M[5] - (M[4] ** 2))  # Det #0
    adj_M[1] = -(M[1] * M[5] - M[4] * M[2])  # Det #1
    adj_M[2] = (M[1] * M[4] - M[3] * M[2])  # Det #2

    # Second row of adjugate matrix
    adj_M[3] = adj_M[1]
    adj_M[4] = (M[0] * M[5] - (M[2] ** 2))
    adj_M[5] = -(M[0] * M[4] - M[1] * M[2])

    # Third row of adjugate matrix
    adj_M[6] = adj_M[2]
    adj_M[7] = adj_M[5]
    adj_M[8] = (M[0] * M[3] - (M[1] ** 2))

    determinant += np.int64(M[0]) * np.int64(adj_M[0])
    determinant += np.int64(M[1]) * np.int64(adj_M[1])  # Using addition since minus is integrated in adjugate matrix.
    determinant += np.int64(M[2]) * np.int64(adj_M[2])

    return adj_M, determinant 

def inverse_3by3_int64(M, return_determinant=True):
    """C style inverse of a flattened 3 by 3 matrix.

    Returns the full adjugate matrix and the determinant, s.t.

    .. math::

        M^{-1} = \\frac{1}{\\det(M)} \\cdot \\text{adj}(M).


    For integer matrices, this then returns exact results by avoiding
    to apply the division.

    :param M: The matrix to find inverse to.
    :type M: :py:class:`numpy.ndarray`
    :return: The adjugate flattened matrix and the determinant.
    :rtype: tuple

    """
    if len(M.shape) > 1:
        M = M.flatten()

    determinant = np.int64(0)
    adj_M = np.zeros((9,), 'int64')

    # First row of adjugate matrix
    adj_M[0] = (M[4] * M[8] - M[7] * M[5])  # Det #0
    adj_M[1] = -(M[1] * M[8] - M[7] * M[2])
    adj_M[2] = (M[1] * M[5] - M[4] * M[2])

    # Second row of adjugate matrix
    adj_M[3] = -(M[3] * M[8] - M[6] * M[5])  # Det #1
    adj_M[4] = (M[0] * M[8] - M[6] * M[2])
    adj_M[5] = -(M[0] * M[5] - M[3] * M[2])

    # Third row of adjugate matrix
    adj_M[6] = (M[3] * M[7] - M[6] * M[4])  # Det #2
    adj_M[7] = -(M[0] * M[7] - M[6] * M[1])
    adj_M[8] = (M[0] * M[4] - M[3] * M[1])

    if return_determinant:
        if ((np.log2(np.abs(M[0])) + np.log2(np.abs(adj_M[0]))) > 63 or
                (np.log2(np.abs(M[1])) + np.log2(np.abs(adj_M[1]))) > 63 or
                (np.log2(np.abs(M[2])) + np.log2(np.abs(adj_M[6]))) > 63):
            print("inverse_3by3_int64: Overflow in determinant calculation!")
            determinant += int(M[0]) * int(adj_M[0])
            determinant += int(M[1]) * int(adj_M[3])  # Using addition since minus is integrated in adjugate matrix.
            determinant += int(M[2]) * int(adj_M[6])
        else:
            determinant += np.int64(M[0]) * np.int64(adj_M[0])
            determinant += np.int64(M[1]) * np.int64(adj_M[3])  # Using addition since minus is integrated in adjugate matrix.
            determinant += np.int64(M[2]) * np.int64(adj_M[6])
        return adj_M, determinant
    else:
        return adj_M 

def _calculate_scatter_matrix_c(x, y):
    """Calculates the upper triangular scatter matrix for the input coordinates.

    :param x: The x coordinates.
    :type x: :py:class:`numpy.ndarray`
    :param y: The y coordinates.
    :type y: :py:class:`numpy.ndarray`
    :return: The upper triangular scatter matrix.
    :rtype: :py:class:`numpy.ndarray`

    """
    S = np.zeros((6, 6), 'int32')

    for i in xrange(len(x)):
        tmp_x2 = x[i] ** 2
        tmp_x3 = tmp_x2 * x[i]
        tmp_y2 = y[i] ** 2
        tmp_y3 = tmp_y2 * y[i]

        S[0, 0] += tmp_x2 * tmp_x2
        S[0, 1] += tmp_x3 * y[i]
        S[0, 2] += tmp_x2 * tmp_y2
        S[0, 3] += tmp_x3
        S[0, 4] += tmp_x2 * y[i]
        S[0, 5] += tmp_x2
        S[1, 2] += tmp_y3 * x[i]
        S[1, 4] += tmp_y2 * x[i]
        S[1, 5] += x[i] * y[i]
        S[2, 2] += tmp_y2 * tmp_y2
        S[2, 4] += tmp_y3
        S[2, 5] += tmp_y2
        S[3, 5] += x[i]
        S[4, 5] += y[i]

    S[5, 5] = len(x)

    # Doubles
    S[1, 1] = S[0, 2]
    S[1, 3] = S[0, 4]
    S[2, 3] = S[1, 4]
    S[3, 3] = S[0, 5]
    S[3, 4] = S[1, 5]
    S[4, 4] = S[2, 5]

    return S 

def _calculate_scatter_matrix_c(x, y):
    """Calculates the upper triangular scatter matrix for the input coordinates.

    :param x: The x coordinates.
    :type x: :py:class:`numpy.ndarray`
    :param y: The y coordinates.
    :type y: :py:class:`numpy.ndarray`
    :return: The upper triangular scatter matrix.
    :rtype: :py:class:`numpy.ndarray`

    """
    S = np.zeros((6, 6), 'int64')

    for i in xrange(len(x)):
        tmp_x2 = x[i] ** 2
        tmp_x3 = tmp_x2 * x[i]
        tmp_y2 = y[i] ** 2
        tmp_y3 = tmp_y2 * y[i]

        S[0, 0] += tmp_x2 * tmp_x2
        S[0, 1] += tmp_x3 * y[i]
        S[0, 2] += tmp_x2 * tmp_y2
        S[0, 3] += tmp_x3
        S[0, 4] += tmp_x2 * y[i]
        S[0, 5] += tmp_x2
        S[1, 2] += tmp_y3 * x[i]
        S[1, 4] += tmp_y2 * x[i]
        S[1, 5] += x[i] * y[i]
        S[2, 2] += tmp_y2 * tmp_y2
        S[2, 4] += tmp_y3
        S[2, 5] += tmp_y2
        S[3, 5] += x[i]
        S[4, 5] += y[i]

    S[5, 5] = len(x)

    # Doubles
    S[1, 1] = S[0, 2]
    S[1, 3] = S[0, 4]
    S[2, 3] = S[1, 4]
    S[3, 3] = S[0, 5]
    S[3, 4] = S[1, 5]
    S[4, 4] = S[2, 5]

    return S 

def test_projection(self):
        """Test projection, query methods of the SOC oracle."""
        x = cvxpy.Variable(1000)
        x_0 = np.ones(1000)
        soc = SOC(x)
        constr = soc._constr

        # Test idempotency
        z = x_0[:-1]
        norm_z = np.linalg.norm(x_0[:-1], 2)
        x_0[-1] = norm_z + 10
        x_star = soc.project(x_0)
        self.assertTrue(soc.contains(x_star))
        self.assertTrue(np.array_equal(x_0, x_star), "projection not "
            "idemptotent")
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.pnorm(x - x_0, 2)), constr)
        p.solve()
        self.assertTrue(np.isclose(np.array(x.value).flatten(), x_star,
            atol=1e-3).all())
        utils.query_helper(self, x_0, x_star, soc, idempotent=True)
        

        # Test the case in which norm_z <= -t
        x_0[-1] = -1 * norm_z
        x_star = soc.project(x_0)
        self.assertTrue(np.array_equal(x_star, np.zeros(1000)),
            "projection not  zero")
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.pnorm(x - x_0, 2)), constr)
        p.solve()
        self.assertTrue(np.isclose(np.array(x.value).flatten(), x_star,
            atol=1e-3).all())
        utils.query_helper(self, x_0, x_star, soc, idempotent=False)

        # Test the case in which norm_z > abs(t)
        x_0[-1] = norm_z - 10
        x_star = soc.project(x_0)
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.pnorm(x - x_0, 2)), constr)
        p.solve()
        self.assertTrue(np.isclose(np.array(x.value).flatten(), x_star,
            atol=1e-3).all())
        utils.query_helper(self, x_0, x_star, soc, idempotent=False)

        # Test random projection
        x_0 = np.random.randn(1000)
        x_star = soc.project(x_0)
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.pnorm(x - x_0, 2)), constr)
        p.solve()
        self.assertTrue(np.isclose(np.array(x.value).flatten(), x_star,
            atol=1e-3).all()) 

def test_projection(self):
        """Test projection, query methods of the NonNeg oracle."""
        x = cvxpy.Variable(1000)
        nonneg = NonNeg(x)
        constr = nonneg._constr

        # Test idempotency
        x_0 = np.abs(np.random.randn(1000))
        x_star = nonneg.project(x_0)
        self.assertTrue(nonneg.contains(x_star))
        self.assertTrue(np.array_equal(x_0, x_star), "projection not "
            "idemptotent")
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.pnorm(x - x_0, 2)), constr)
        p.solve()
        self.assertTrue(np.isclose(np.array(x.value).flatten(), x_star,
            atol=1e-3).all())
        utils.query_helper(self, x_0, x_star, nonneg, idempotent=True)
        

        # Test the case in which x_0 <= 0
        x_0 = -1 * np.abs(np.random.randn(1000))
        x_star = nonneg.project(x_0)
        self.assertTrue(np.array_equal(x_star, np.zeros(1000)),
            "projection not zero")
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.pnorm(x - x_0, 2)), constr)
        p.solve()
        self.assertTrue(np.isclose(np.array(x.value).flatten(), x_star,
            atol=1e-3).all())
        utils.query_helper(self, x_0, x_star, nonneg, idempotent=False)

        # Test the case in which x_0 == 0 (idempotent as well)
        x_0 = np.zeros(1000)
        x_star = nonneg.project(x_0)
        self.assertTrue(np.array_equal(x_star, np.zeros(1000)),
            "projection not zero")
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.pnorm(x - x_0, 2)), constr)
        p.solve()
        self.assertTrue(np.isclose(np.array(x.value).flatten(), x_star,
            atol=1e-3).all())
        utils.query_helper(self, x_0, x_star, nonneg, idempotent=True)

        # Test random projection
        x_0 = np.random.randn(1000)
        x_star = nonneg.project(x_0)
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.pnorm(x - x_0, 2)), constr)
        p.solve()
        self.assertTrue(np.isclose(np.array(x.value).flatten(), x_star,
            atol=1e-3).all()) 

def test_projection(self):
        """Test projection, query methods of the Zeros oracle."""
        x = cvxpy.Variable(1000)
        zeros = Zeros(x)
        constr = zeros._constr

        # Test idempotency
        x_0 = np.zeros(1000)
        x_star = zeros.project(x_0)
        self.assertTrue(zeros.contains(x_star))
        self.assertTrue(np.array_equal(x_0, x_star), "projection not "
            "idemptotent")
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.pnorm(x - x_0, 2)), constr)
        p.solve()
        self.assertTrue(np.isclose(np.array(x.value).flatten(), x_star,
            atol=1e-3).all())
        utils.query_helper(self, x_0, x_star, zeros, idempotent=True)
        

        # Test the case in which x_0 >= 0
        x_0 = np.abs(np.random.randn(1000))
        x_star = zeros.project(x_0)
        self.assertTrue(np.array_equal(x_star, np.zeros(1000)),
            "projection not  zero")
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.pnorm(x - x_0, 2)), constr)
        p.solve()
        self.assertTrue(np.isclose(np.array(x.value).flatten(), x_star,
            atol=1e-3).all())
        utils.query_helper(self, x_0, x_star, zeros, idempotent=False)

        # Test the case in which x_0 <= 0
        x_0 = -1 * np.abs(np.random.randn(1000))
        x_star = zeros.project(x_0)
        self.assertTrue(np.array_equal(x_star, np.zeros(1000)),
            "projection not  zero")
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.pnorm(x - x_0, 2)), constr)
        p.solve()
        self.assertTrue(np.isclose(np.array(x.value).flatten(), x_star,
            atol=1e-3).all())
        utils.query_helper(self, x_0, x_star, zeros, idempotent=False)

        # Test random projection
        x_0 = np.random.randn(1000)
        x_star = zeros.project(x_0)
        self.assertTrue(np.array_equal(x_star, np.zeros(1000)),
            "projection not  zero")
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.pnorm(x - x_0, 2)), constr)
        p.solve()
        self.assertTrue(np.isclose(np.array(x.value).flatten(), x_star,
            atol=1e-3).all()) 

def bb_coords_to_grid(bb_coords_one, img_shape, grid_size):
    """Convert bounding box coordinates (corners) to ground truth heatmaps."""
    if isinstance(bb_coords_one, ia.KeypointsOnImage):
        bb_coords_one = bb_coords_one.keypoints

    # bb edges after augmentation
    x1b = min([kp.x for kp in bb_coords_one])
    x2b = max([kp.x for kp in bb_coords_one])
    y1b = min([kp.y for kp in bb_coords_one])
    y2b = max([kp.y for kp in bb_coords_one])

    # clip
    x1c = np.clip(x1b, 0, img_shape[1]-1)
    y1c = np.clip(y1b, 0, img_shape[0]-1)
    x2c = np.clip(x2b, 0, img_shape[1]-1)
    y2c = np.clip(y2b, 0, img_shape[0]-1)

    # project
    x1d = int((x1c / img_shape[1]) * grid_size)
    y1d = int((y1c / img_shape[0]) * grid_size)
    x2d = int((x2c / img_shape[1]) * grid_size)
    y2d = int((y2c / img_shape[0]) * grid_size)

    assert 0 <= x1d < grid_size
    assert 0 <= y1d < grid_size
    assert 0 <= x2d < grid_size
    assert 0 <= y2d < grid_size

    # output ground truth:
    # - 1 heatmap that is 1 everywhere where there is a bounding box
    # - 9 position sensitive heatmaps,
    #   e.g. the first one is 1 everywhere where there is the _top left corner_
    #        of a bounding box,
    #        the second one is 1 for the top center cell,
    #        the third one is 1 for the top right corner,
    #        ...
    grids = np.zeros((grid_size, grid_size, 1+9), dtype=np.float32)
    # first heatmap
    grids[y1d:y2d+1, x1d:x2d+1, 0] = 1
    # position sensitive heatmaps
    nb_cells_x = 3
    nb_cells_y = 3
    cell_width = (x2d - x1d) / nb_cells_x
    cell_height = (y2d - y1d) / nb_cells_y
    cell_counter = 0
    for j in range(nb_cells_y):
        cell_y1 = y1d + cell_height * j
        cell_y2 = cell_y1 + cell_height
        cell_y1_int = np.clip(int(math.floor(cell_y1)), 0, img_shape[0]-1)
        cell_y2_int = np.clip(int(math.floor(cell_y2)), 0, img_shape[0]-1)
        for i in range(nb_cells_x):
            cell_x1 = x1d + cell_width * i
            cell_x2 = cell_x1 + cell_width
            cell_x1_int = np.clip(int(math.floor(cell_x1)), 0, img_shape[1]-1)
            cell_x2_int = np.clip(int(math.floor(cell_x2)), 0, img_shape[1]-1)
            grids[cell_y1_int:cell_y2_int+1, cell_x1_int:cell_x2_int+1, 1+cell_counter] = 1
            cell_counter += 1
    return grids 

def stiffnessmatrix(self):
        """Stiffness matrix of a triangular mesh

        The method determines a stiffness matrix of a triangular mesh.

        Returns
        -------
        stiffmatrix : array, shape (n_points, n_points)
            Symmetric matrix, which contains the stiffness values for each edge
            and vertex for the FEM approch. The number of vertices of the
            triangular mesh is n_points.

        Examples
        --------

        >>> from spharapy import trimesh as tm
        >>> testtrimesh = tm.TriMesh([[0, 1, 2]], [[1., 0., 0.], [0., 2., 0.],
        ...                                        [0., 0., 3.]])
        >>> testtrimesh.stiffnessmatrix()
        array([[-0.92857143,  0.64285714,  0.28571429],
               [ 0.64285714, -0.71428571,  0.07142857],
               [ 0.28571429,  0.07142857, -0.35714286]])

        References
        ----------
        :cite:`vallet07`

        """

        # get the largest index from from triangle list
        maxindex = self.trilist.max()

        # fill the weight matrix with zeros, the size is (maxindex + 1)^2
        stiffmatrix = np.zeros((maxindex + 1, maxindex + 1), dtype=float)

        # compute the cot weight matrix
        weightmatrix = self.weightmatrix(mode='half_cotangent')

        # compute and return the stiffness matrix
        stiffmatrix = -np.diag(weightmatrix.sum(axis=0)) + weightmatrix

        return stiffmatrix 

def roll():
    parser = argparse.ArgumentParser(
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
        description=(
            """
            Parses a list of dice rolls and modifiers, and prints
            the total of one realization, the individual rolls that
            produced that total, and the mean/percentile of that roll
            relative to the overall distribution.

            Example Usage:
            roll 1d10 + 3d6 + 8

            
            """))
    parser.add_argument('rolls', nargs='+', help=(
        """The list of rolls and modifiers to simulate."""))
    
    args = parser.parse_args().rolls
    message1 = "= "
    message2 = ""
    total = 0
    means = 0
    distribution = np.zeros(n_trials)
    for arg in args:
        if arg == '+':
            message1 += ' + '
            message2 += ' + '
        else:
            the_sum, rolls, mean, realizations = parse_roll(arg)
            message1 += str(the_sum) + ' '*(len(rolls) - len(str(the_sum)))
            message2 += rolls
            total += the_sum
            means += mean
            distribution += realizations
    print('*************')
    print('Total: {:4d}'.format(total))
    print('*************')
    print(message1)
    print(message2)
    print('Mean: {:.1f}'.format(means))
    print('Percentile: {:.1f}%'.format((total > distribution).mean()*100.)) 

def estimate_from_dwis(data, axis=-2, return_mask=False, exclude_mask=None, ncores=-1, method='moments', verbose=0):
    '''Given the data, splits over each slice to compute parameters of the gamma distribution

    input
    ------
        data : array, input volume used to identify noise voxels and estimate the noise distribution

    optional
    --------
        axis : (0, 1 or 2) the axis to consider as a slab of uniform noise profile

        return_mask : bool, if True returns the identified noise voxels as a mask

        exclude_mask : array, mask indicating voxels to remove from all the computations, such as those containing huge artifacts.

        ncores : int, number of cores to use for multiprocessing

        method='moments' or method='maxlk' : which algorithm to use to estimate sigma and N

        verbose : print progress for joblib, can be an integer to increase verbosity

    output
    -------
    sigma, N, mask (optional)
    '''

    # guess a gross upper bound of sigma
    # if it's masked, median can be zero, so use the nonzero data
    median = np.median(data)

    if median == 0:
        median = np.median(data[data > 0])

    if axis < 0:
        axis = data.ndim + axis

    ranger = range(data.shape[axis])
    swapped_data = data.swapaxes(0, axis)

    if exclude_mask is None:
        exclude_mask = [None] * swapped_data.shape[0]
    else:
        exclude_mask = exclude_mask.swapaxes(0, axis)

    output = Parallel(n_jobs=ncores, verbose=verbose)(delayed(_inner)(swapped_data[i], median, exclude_mask[i], method) for i in ranger)

    # output is each slice we took along axis, so the mask might be reversed
    sigma = np.zeros(len(output))
    N = np.zeros(len(output))
    mask = np.zeros(data.shape[:-1], dtype=np.int16).swapaxes(0, axis)

    for i, s in enumerate(output):
        sigma[i] = s[0]
        N[i] = s[1]
        mask[i] = s[2]

    if return_mask:
        return sigma, N, mask.swapaxes(0, axis)
    return sigma, N 

def vectorize(ad_functions, n_inputs = 1, n_vectors = 1):
	'''
	INPUTS
	======
	ad_functions: 	a list or row vector of AutoDiff objects
	n_inputs:		number of input variables the final function will use
	n_vectors: 		length of vector for input variables

	RETURNS
	=======
	An AutoDiff object of vector functions with calculated value, derivative, and jacobian.

	EXAMPLES
	========
	>>> x = AutoDiff(3, np.array([[2, 0, 0]]), n=3, k=1)
	>>> y = AutoDiff(2, np.array([[0, 2, 0]]), n=3, k=2)
	>>> z = AutoDiff(-1, np.array([[0, 0, 2]]), n=3, k=3)
	>>> fs = [3*x + 2*y + 4*z, x - y + z, x/2, 2*x - 2*y]
	>>> f = vectorize(fs, 3)
	>>> f.val
	np.array([[9],[0],[1.5],[2]]
	>>> f.der
	np.array([[6, 4, 8], [2, -2, 2], [1, 0, 0], [4, -4, 0]])
	>>> f.jacobian
	np.array([[3, 2, 4], [1, -1, 1], [0.5, 0, 0], [2, -2, 0]])
	'''
	if n_vectors == 1:
		val = np.zeros([len(ad_functions), n_vectors])
		der = np.zeros([len(ad_functions), n_inputs])
		jacobian = np.zeros([len(ad_functions), n_inputs])
		for i, f in enumerate(ad_functions):
			val[i, :] = f.val
			der[i, :] = f.der
			jacobian[i, :] = f.jacobian
	else:
		val = np.zeros([len(ad_functions), n_vectors])
		der = np.zeros([n_vectors, len(ad_functions), n_inputs])
		jacobian = np.zeros([n_vectors, len(ad_functions), n_inputs])
		for j in range(n_vectors):
			der_j = der[j]
			jac_j = jacobian[j]
			for i, f in enumerate(ad_functions):
				val[i, :] = f.val.T
				der_j[i, :] = f.der[j]
				jac_j[i, :]	= f.jacobian[j]

	return AutoDiff(val, der, n_inputs, 0, jacobian)