Python numpy.sqrt() Examples

def great_circle_dist(p1, p2):
    """Return the distance (in km) between two points in
    geographical coordinates.
    """
    lon0, lat0 = p1
    lon1, lat1 = p2
    EARTH_R = 6372.8
    lat0 = np.radians(float(lat0))
    lon0 = np.radians(float(lon0))
    lat1 = np.radians(float(lat1))
    lon1 = np.radians(float(lon1))
    dlon = lon0 - lon1
    y = np.sqrt(
        (np.cos(lat1) * np.sin(dlon)) ** 2
        + (np.cos(lat0) * np.sin(lat1)
           - np.sin(lat0) * np.cos(lat1) * np.cos(dlon)) ** 2)
    x = np.sin(lat0) * np.sin(lat1) + \
        np.cos(lat0) * np.cos(lat1) * np.cos(dlon)
    c = np.arctan2(y, x)
    return EARTH_R * c 

def update_core_cpu(self, param):
        grad = param.grad
        if grad is None:
            return
        hp = self.hyperparam
        eps = grad.dtype.type(hp.eps)
        if hp.eps != 0 and eps == 0:
            raise ValueError(
                'eps of Adam optimizer is too small for {} ({})'.format(
                    grad.dtype.name, hp.eps))
        m, v = self.state['m'], self.state['v']

        m += (1 - hp.beta1) * (grad - m)
        v += (1 - hp.beta2) * (grad * grad - v)

        vhat = v
        # This adam multipies schduled adaptive learning rate
        # with both main term and weight decay.
        # Normal Adam: param.data -= hp.eta * (self.lr * m / (numpy.sqrt(vhat) + hp.eps) +
        #                                      hp.weight_decay_rate * param.data)
        param.data -= hp.eta * self.lr * (m / (numpy.sqrt(vhat) + hp.eps) +
                                          hp.weight_decay_rate * param.data) 

def update_core_gpu(self, param):
        grad = param.grad
        if grad is None:
            return

        hp = self.hyperparam
        eps = grad.dtype.type(hp.eps)
        if hp.eps != 0 and eps == 0:
            raise ValueError(
                'eps of Adam optimizer is too small for {} ({})'.format(
                    grad.dtype.name, hp.eps))

        cuda.elementwise(
            'T grad, T lr, T one_minus_beta1, T one_minus_beta2, T eps, \
             T eta, T weight_decay_rate',
            'T param, T m, T v',
            '''m += one_minus_beta1 * (grad - m);
               v += one_minus_beta2 * (grad * grad - v);
               param -= eta * lr * (m / (sqrt(v) + eps) +
                               weight_decay_rate * param);''',
            'adam')(grad, self.lr, 1 - hp.beta1,
                    1 - hp.beta2, hp.eps,
                    hp.eta, hp.weight_decay_rate,
                    param.data, self.state['m'], self.state['v']) 

def calculate_diff_stress(self, x, u, nu, side=1):
        """
        Calculate the derivative of the Von Mises stress given the densities x,
        displacements u, and young modulus nu. Optionally, provide the side
        length (default: 1).
        """
        rho = self.penalized_densities(x)
        EB = self.E(nu).dot(self.B(side))
        EBu = sum([EB.dot(u[:, i][self.edofMat]) for i in range(u.shape[1])])
        s11, s22, s12 = numpy.hsplit((EBu * rho / float(u.shape[1])).T, 3)
        drho = self.diff_penalized_densities(x)
        ds11, ds22, ds12 = numpy.hsplit(
            ((1 - rho) * drho * EBu / float(u.shape[1])).T, 3)
        vm_stress = numpy.sqrt(s11**2 - s11 * s22 + s22**2 + 3 * s12**2)
        if abs(vm_stress).sum() > 1e-8:
            dvm_stress = (0.5 * (1. / vm_stress) * (2 * s11 * ds11 -
                ds11 * s22 - s11 * ds22 + 2 * s22 * ds22 + 6 * s12 * ds12))
            return dvm_stress
        return 0 

def test_arcsinh_ad_results():
	# positive real numbers
	x = AutoDiff(1, 2)
	f = ef.arcsinh(x)
	assert f.val == np.arcsinh(1)
	assert f.der == np.array([[((2)/np.sqrt((1)**2 + 1))]])
	assert f.jacobian == np.array([[((1)/np.sqrt((1)**2 + 1))]])
	# negative real numbers
	y = AutoDiff(-1, 2)
	f = ef.arcsinh(y)
	assert f.val == np.arcsinh(-1)
	assert f.der == np.array([[((2)/np.sqrt((-1)**2 + 1))]])
	assert f.jacobian == np.array([[((1)/np.sqrt((-1)**2 + 1))]])
	# zero
	z = AutoDiff(0, 2)
	f = ef.arcsinh(z)
	assert f.val == np.arcsinh(0)
	assert f.der == np.array([[((2)/np.sqrt((0)**2 + 1))]])
	assert f.jacobian == np.array([[((1)/np.sqrt((0)**2 + 1))]]) 

def test_arccosh_ad_results():
	# value defined at positive real numbers x >= 1
	# derivative defined at positive real numbers x > 1
	x = AutoDiff(1.1, 2)
	f = ef.arccosh(x)
	assert f.val == np.arccosh(1.1)
	assert f.der == np.array([[((2)/np.sqrt((1.1)**2 - 1))]])
	assert f.jacobian == np.array([[((1)/np.sqrt((1.1)**2 - 1))]])
	# value defined at x = 1, derivative not defined
	with pytest.warns(RuntimeWarning):
		y = AutoDiff(1, 2)
		f = ef.arccosh(y)
		assert np.isinf(f.der)
		assert np.isinf(f.jacobian)
	# neither value nor derivative defined at x < 1
	with pytest.warns(RuntimeWarning):
		z = AutoDiff(0, 2)
		f = ef.arccosh(z)
		assert np.isnan(f.val)
		assert np.isnan(f.der)
		assert np.isnan(f.jacobian) 

def test_sqrt_ad_results():
	# Positive reals
	x = AutoDiff(0.5, 2.0)
	f = ef.sqrt(x)
	assert f.val == np.array([[np.sqrt(0.5)]])
	assert f.der == np.array([[0.5 * 0.5 ** (-0.5) * 2.0]])
	assert f.jacobian == np.array([[0.5 * 0.5 ** (-0.5) * 1]])
	# Value defined but derivative undefined when x == 0
	with pytest.warns(RuntimeWarning):
		y = AutoDiff(0, 2)
		f = ef.sqrt(y)
		assert f.val == np.array([[0]])
		assert np.isinf(f.der[0][0])
		assert np.isinf(f.jacobian[0][0])
	# Value and derivative undefined when x < 0
	with pytest.warns(RuntimeWarning):
		z = AutoDiff(-0.5, 2)
		f = ef.sqrt(z)
		assert np.isnan(f.val[0][0])
		assert np.isnan(f.der[0][0])
		assert np.isnan(f.jacobian[0][0]) 

def test_arcsinh_ad_results():
	# positive real numbers
	x = Dual(1, 2)
	f = ef.arcsinh(x)
	assert f.Real == np.arcsinh(1)
	assert f.Dual == np.array([[((2)/np.sqrt((1)**2 + 1))]])
	
	# negative real numbers
	y = Dual(-1, 2)
	f = ef.arcsinh(y)
	assert f.Real == np.arcsinh(-1)
	assert f.Dual == np.array([[((2)/np.sqrt((-1)**2 + 1))]])
	
	# zero
	z = Dual(0, 2)
	f = ef.arcsinh(z)
	assert f.Real == np.arcsinh(0)
	assert f.Dual == np.array([[((2)/np.sqrt((0)**2 + 1))]]) 

def test_arccosh_ad_results():
	# Realue defined at positive real numbers x >= 1
	# Dualivative defined at positive real numbers x > 1
	x = Dual(1.1, 2)
	f = ef.arccosh(x)
	assert f.Real == np.arccosh(1.1)
	assert f.Dual == np.array([[((2)/np.sqrt((1.1)**2 - 1))]])
	
	# Realue defined at x = 1, Dualivative not defined
	with pytest.warns(RuntimeWarning):
		y = Dual(1, 2)
		f = ef.arccosh(y)
		assert np.isinf(f.Dual)
	
	# neither Realue nor Dualivative defined at x < 1
	with pytest.warns(RuntimeWarning):
		z = Dual(0, 2)
		f = ef.arccosh(z)
		assert np.isnan(f.Real)
		assert np.isnan(f.Dual) 

def test_sqrt_ad_results():
	# Positive reals
	x = Dual(0.5, 2.0)
	f = ef.sqrt(x)
	assert f.Real == np.array([[np.sqrt(0.5)]])
	assert f.Dual == np.array([[0.5 * 0.5 ** (-0.5) * 2.0]])

	# Realue defined but Dualivative undefined when x == 0
	with pytest.warns(RuntimeWarning):
		y = Dual(0, 2)
		f = ef.sqrt(y)
		assert f.Real == np.array([[0]])
		assert np.isinf(f.Dual)

	# Realue and Dualivative undefined when x < 0
	with pytest.warns(RuntimeWarning):
		z = Dual(-0.5, 2)
		f = ef.sqrt(z)
		assert np.isnan(f.Real)
		assert np.isnan(f.Dual) 

def get_fans(shape, dim_ordering='th'):
    if len(shape) == 2:
        fan_in = shape[0]
        fan_out = shape[1]
    elif len(shape) == 4 or len(shape) == 5:
        # assuming convolution kernels (2D or 3D).
        # TH kernel shape: (depth, input_depth, ...)
        # TF kernel shape: (..., input_depth, depth)
        if dim_ordering == 'th':
            receptive_field_size = np.prod(shape[2:])
            fan_in = shape[1] * receptive_field_size
            fan_out = shape[0] * receptive_field_size
        elif dim_ordering == 'tf':
            receptive_field_size = np.prod(shape[:2])
            fan_in = shape[-2] * receptive_field_size
            fan_out = shape[-1] * receptive_field_size
        else:
            raise ValueError('Invalid dim_ordering: ' + dim_ordering)
    else:
        # no specific assumptions
        fan_in = np.sqrt(np.prod(shape))
        fan_out = np.sqrt(np.prod(shape))
    return fan_in, fan_out 

def _azimuthal_wvnum(k, l, N, nfactor):
    k = k.values
    l = l.values
    K = np.sqrt(k[np.newaxis,:]**2 + l[:,np.newaxis]**2)
    nbins = int(N/nfactor)
    if k.max() > l.max():
        ki = np.linspace(0., l.max(), nbins)
    else:
        ki = np.linspace(0., k.max(), nbins)

    kidx = np.digitize(np.ravel(K), ki)
    area = np.bincount(kidx)

    kr = np.bincount(kidx, weights=K.ravel()) / area

    return kidx, area, kr 

def classical_mds(self, D):
        ''' 
        Classical multidimensional scaling

        Parameters
        ----------
        D : square 2D ndarray
            Euclidean Distance Matrix (matrix containing squared distances between points
        '''

        # Apply MDS algorithm for denoising
        n = D.shape[0]
        J = np.eye(n) - np.ones((n,n))/float(n)
        G = -0.5*np.dot(J, np.dot(D, J))

        s, U = np.linalg.eig(G)

        # we need to sort the eigenvalues in decreasing order
        s = np.real(s)
        o = np.argsort(s)
        s = s[o[::-1]]
        U = U[:,o[::-1]]

        S = np.diag(s)[0:self.dim,:]
        self.X = np.dot(np.sqrt(S),U.T) 

def trilateration(self, D):
        '''
        Find the location of points based on their distance matrix using trilateration

        Parameters
        ----------
        D : square 2D ndarray
            Euclidean Distance Matrix (matrix containing squared distances between points
        '''

        dist = np.sqrt(D)

        # Simpler algorithm (no denoising)
        self.X = np.zeros((self.dim, self.m))

        self.X[:,1] = np.array([0, dist[0,1]])
        for i in xrange(2,m):
            self.X[:,i] = self.trilateration_single_point(self.X[1,1],
                    dist[0,i], dist[1,i]) 

def mtx_freq2visi(M, p_mic_x, p_mic_y):
    """
    build the matrix that maps the Fourier series to the visibility
    :param M: the Fourier series expansion is limited from -M to M
    :param p_mic_x: a vector that constains microphones x coordinates
    :param p_mic_y: a vector that constains microphones y coordinates
    :return:
    """
    num_mic = p_mic_x.size
    ms = np.reshape(np.arange(-M, M + 1, step=1), (1, -1), order='F')
    G = np.zeros((num_mic * (num_mic - 1), 2 * M + 1), dtype=complex, order='C')
    count_G = 0
    for q in range(num_mic):
        p_x_outer = p_mic_x[q]
        p_y_outer = p_mic_y[q]
        for qp in range(num_mic):
            if not q == qp:
                p_x_qqp = p_x_outer - p_mic_x[qp]
                p_y_qqp = p_y_outer - p_mic_y[qp]
                norm_p_qqp = np.sqrt(p_x_qqp ** 2 + p_y_qqp ** 2)
                phi_qqp = np.arctan2(p_y_qqp, p_x_qqp)
                G[count_G, :] = (-1j) ** ms * sp.special.jv(ms, norm_p_qqp) * \
                                np.exp(1j * ms * phi_qqp)
                count_G += 1
    return G 

def draw(self, ax, color, line_width=1, fillcolor=None, name=None, arrow=True, alpha=0.2, scale=50):
        ax.add_patch(PolygonPatch(self.contour, alpha=alpha, fc=fillcolor, ec=color, linewidth=line_width))

        vertices = np.array(self.contour.exterior.coords)[1:]

        if arrow:
            arrow_center = np.mean(vertices, axis=0)
            arrow_direction = (vertices[2] - vertices[1]) / 1.5
            arrow_tail = arrow_center - arrow_direction / 2
            arrow_head = arrow_center + arrow_direction / 2
            style = plt_patches.ArrowStyle.Simple(head_length=.4, head_width=.6, tail_width=.1)
            x = np.array(ax.axis())
            scale_factor = np.sqrt(np.prod(np.abs(x[::2] - x[1::2])) / (60 * 60))
            arrow_patch = plt_patches.FancyArrowPatch(posA=arrow_tail, posB=arrow_head, arrowstyle=style,
                                                      color='w', mutation_scale= scale / scale_factor, alpha=0.4)
            ax.add_patch(arrow_patch)
        elif name is None:
            name = 'front'

        if name is not None:
            text_location = np.mean(vertices[[0, -1]], axis=0)
            ax.text(text_location[0], text_location[1], name, ha='center', va='top', color='w') 

def plot(samples):
    width = min(12,int(np.sqrt(len(samples))))
    fig = plt.figure(figsize=(width, width))
    gs = gridspec.GridSpec(width, width)
    gs.update(wspace=0.05, hspace=0.05)

    for ind, sample in enumerate(samples):
        if ind >= width*width:
            break
        ax = plt.subplot(gs[ind])
        plt.axis('off')
        ax.set_xticklabels([])
        ax.set_yticklabels([])
        ax.set_aspect('equal')
        sample = sample * 0.5 + 0.5
        sample = np.transpose(sample, (1, 2, 0))
        plt.imshow(sample)

    return fig 

def solve_modal(model,k:int):
    """
    Solve eigen mode of the MDOF system
    
    params:
        model: FEModel.
        k: number of modes to extract.
    """
    K_,M_=model.K_,model.M_
    if k>model.DOF:
        logger.info('Warning: the modal number to extract is larger than the system DOFs, only %d modes are available'%model.DOF)
        k=model.DOF
    omega2s,modes = sl.eigsh(K_,k,M_,sigma=0,which='LM')
    delta = modes/np.sum(modes,axis=0)
    model.is_solved=True
    model.mode_=delta
    model.omega_=np.sqrt(omega2s).reshape((k,1)) 

def tile_image(x_gen, tiles=None):
    """Tiled image representations.

    Args:
      x_gen: 4D array of images (n x w x h x 3)
      tiles (int pair, optional): number of rows and columns

    Returns:
      Array of tiled images (1 x W x H x 3)
    """
    n_images = x_gen.shape[0]
    if tiles is None:
        for i in range(int(np.sqrt(n_images)), 0, -1):
            if n_images % i == 0: break
        n_rows = i; n_cols = n_images // i
    else:
        n_rows, n_cols = tiles
    full = [np.hstack(x_gen[c * n_rows:(c + 1) * n_rows]) for c in range(n_cols)]
    return np.expand_dims(np.vstack(full), 0) 

def Givens_rotation_double(a, b):
    # Working with actual trigonometric functions
    # angle = np.arctan2(-a, b)
    # c = np.cos(angle)
    # s = np.sin(angle)

    # Using naive definitions
    # root = np.sqrt(a ** 2 + b ** 2)
    # c = a / root
    # s = -b / root

    # Using Matrix Computations solution
    if b == 0:
        c = 1.0
        s = 0.0
    else:
        if np.abs(b) > np.abs(a):
            tau = - a / b
            s = 1 / (np.sqrt(1 + tau ** 2))
            c = s * tau
        else:
            tau = - b / a
            c = 1 / (np.sqrt(1 + tau ** 2))
            s = c * tau

    return c, s 

def conic_to_general_reference(conic_coeffs):
    """Transform from conic section format to general format.

    :param conic_coeffs: The six coefficients defining the ellipse as a conic shape.
    :type conic_coeffs: :py:class:`numpy.ndarray` or tuple
    :param verbose: If debug printout is desired.
    :type verbose: bool
    :return:
    :rtype: tuple

    """
    a, b, c, d, e, f = conic_coeffs

    angle = np.arctan2(b, a - c) / 2

    cos_theta = np.cos(angle)  # np.sqrt((1 + np.cos(2*angle)) / 2)
    sin_theta = np.sin(angle)  # np.sqrt((1 - np.cos(2*angle)) / 2)

    a_prime = a * (cos_theta ** 2) + (b * cos_theta * sin_theta) + c * (sin_theta ** 2)
    c_prime = a * (sin_theta ** 2) - (b * cos_theta * sin_theta) + c * (cos_theta ** 2)
    d_prime = (d * cos_theta) + (e * sin_theta)
    e_prime = (-(d * sin_theta)) + (e * cos_theta)
    f_prime = f

    x_prime = (-d_prime) / (2 * a_prime)
    y_prime = (-e_prime) / (2 * c_prime)
    major_axis = np.sqrt(
        ((-4 * f_prime * a_prime * c_prime) + (c_prime * (d_prime ** 2)) + (a_prime * (e_prime ** 2))) /
        (4 * a_prime * (c_prime ** 2)))
    minor_axis = np.sqrt(
        ((-4 * f_prime * a_prime * c_prime) + (c_prime * (d_prime ** 2)) + (a_prime * (e_prime ** 2))) /
        (4 * (a_prime ** 2) * c_prime))

    if a_prime > c_prime:
        angle += np.pi / 2

    x = x_prime * cos_theta - y_prime * sin_theta
    y = x_prime * sin_theta + y_prime * cos_theta

    return [x, y], [major_axis, minor_axis], angle 

def conic_to_general_2(conic_coeffs):
    """Transform from conic section format to general format.

    :param conic_coeffs: The six coefficients defining the ellipse as a conic shape.
    :type conic_coeffs: :py:class:`numpy.ndarray` or tuple
    :return: The general form for the ellipse. Returns tuple :math:`(x_c,\ y_c),\\ (a,\\ b),\\ \\theta`
        that fulfills the equation

        .. math::

            \\frac{((x-x_c)\\cos(\\theta) + (y-y_c)\\sin(\\theta))^2}{a^2} +
            \\frac{((x-x_c)\\sin(\\theta) - (y-y_c)\\sin(\\theta))^2}{b^2} = 1

    :rtype: tuple

    """
    a, b, c, d, e, f = conic_coeffs
    denom = 2 * ((b ** 2) - (a * c))
    x = (c * d - b * e) / denom
    y = (a * e - b * d) / denom
    mu = 1 / ((a * (x ** 2)) + (2 * b * x * y) + (c * (y ** 2)) - f)

    sqrt_expr = np.sqrt(((mu * a - mu * c) ** 2) + (4 * ((mu * b) ** 2)))
    min_axis = 1 / np.sqrt((mu * a + mu * c + sqrt_expr) / 2)
    maj_axis = 1 / np.sqrt((mu * a + mu * c - sqrt_expr) / 2)
    angle = np.arctan2(-2 * b, c - a)

    return [x, y], [maj_axis, min_axis], angle 

def inverse_iteration_for_eigenvector_int(A, eigenvalue):
    """Performs a series of inverse iteration steps with a known
    eigenvalue to produce its eigenvector.

    :param A: The 3x3 matrix to which the eigenvalue belongs.
    :type A: :py:class:`numpy.ndarray`
    :param eigenvalue: One approximate eigenvalue of the matrix A.
    :type eigenvalue: int
    :return: The eigenvector of this matrix and eigenvalue combination.
    :rtype: :py:class:`numpy.ndarray`

    """
    A = np.array(A, 'int64')

    # Subtract the eigenvalue from the diagonal entries of the matrix.
    for k in xrange(A.shape[0]):
        A[k, k] -= eigenvalue
    A, scale = fp.scale_64bit_matrix(A)

    # Obtain the inverse of the matrix.
    adj_A = mo.inverse_3by3_int64(A.flatten(), False)
    eigenvector = adj_A.reshape((3, 3)).sum(1)
    eigenvector, scale = fp.scale_64bit_vector(eigenvector)

    e_norm = int(np.sqrt((eigenvector ** 2).sum()))
    if (eigenvector[0] < 0) and (eigenvector[2] < 0):
        eigenvector = -eigenvector
    return eigenvector, e_norm 

def _scheduled_learning_rate(schedule, hp, t):
    if t == 0:
        raise RuntimeError(
            'Can\'t determine the learning rate of Adam optimizer '
            'because the update steps have not been started.')
    fix1 = 1. - math.pow(hp.beta1, t)
    fix2 = 1. - math.pow(hp.beta2, t)
    lrt = hp.alpha * math.sqrt(fix2) / fix1
    lrt *= schedule(t / hp.t_total, hp.warmup)
    return lrt 

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 calculate_stress(self, x, u, nu, side=1):
        """
        Calculate the Von Mises stress given the densities x, displacements u,
        and young modulus nu.
        """
        s11, s22, s12 =  self.calculate_principle_stresses(x, u, nu, side)
        vm_stress = numpy.sqrt(s11**2 - s11 * s22 + s22**2 + 3 * s12**2)
        return vm_stress 

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 calculate_stress(self, x, u, nu, side=1):
        """
        Calculate the Von Mises stress given the densities x, displacements u,
        and young modulus nu.
        """
        s11, s22, s12 =  self.calculate_principle_stresses(x, u, nu, side)
        vm_stress = numpy.sqrt(s11**2 - s11 * s22 + s22**2 + 3 * s12**2)
        return vm_stress 

def arccos(X):
	''' Compute the arccos of an AutoDiff object and its derivative.

		INPUTS
		======
		X: an AutoDiff object or constant

		RETURNS
		=======
		A new AutoDiff object or scalar with calculated value and derivative.

		EXAMPLES
		========
		>>> X = AutoDiff(0.5, 2)
		>>> arccosAutoDiff = arccos(X)
		>>> arccosAutoDiff.val
		1.0471975511965976
		>>> arccosAutoDiff.der
		-2.3094010767585034
		>>> arccosAutoDiff.jacobian
		-1.1547005383792517
		'''
	try:
		# Is another ADT
		new_val = np.arccos(X.val) #if (-1 <= X.val and X.val <= 1) else np.nan
		new_der = (-1/np.sqrt(1-X.val**2))*X.der #if (-1 < X.val and X.val < 1) else np.nan
		new_jacobian = (-1/np.sqrt(1-X.val**2))*X.jacobian #if (-1 < X.val and X.val < 1) else np.nan

		return AutoDiff(new_val, new_der, X.n, 0, new_jacobian)
	except AttributeError:
		try:
			return Dual(np.arccos(X.Real), -X.Dual/np.sqrt(1-X.Real**2))		
		except AttributeError:
			try:
				return Dual(arccos(X.Real), -X.Dual/sqrt(1-X.Real**2))
			except AttributeError:
			# Constant
				return_val = np.arccos(X)
				return return_val

# arc tangent 

def arcsinh(x):
	''' Compute the hyperbolic arc sine of an AutoDiff object and its derivative.
	
	INPUTS
	======
	x: an AutoDiff object
	
	RETURNS
	=======
	A new AutoDiff object with calculated value and derivative.
	
	EXAMPLES
	========
	>>> x = AutoDiff(0.5, 2, 1)
	>>> myAutoDiff = arcsinh(x)
	>>> myAutoDiff.val
	2.3124383412727525
	>>> myAutoDiff.der
	0.39223227027
	>>> myAutoDiff.jacobian
	0.19611613513818404
	
	'''
	try:
		new_val = np.arcsinh(x.val)
		new_der = ((1)/np.sqrt(x.val**2 + 1))*x.der
		new_jacobian = ((1)/np.sqrt(x.val**2 + 1))*x.jacobian
		return AutoDiff(new_val, new_der, x.n, 0, new_jacobian)
	except AttributeError:
		try:
			return Dual(np.arcsinh(x.Real), x.Dual/np.sqrt((x.Real**2)+1))		
		except AttributeError:
			try:
				return Dual(arcsinh(x.Real), (x.Dual*(1+x.Real**2)**-0.5))
			except AttributeError:
			# Constant
				return_val = np.arcsinh(x)
				return return_val

# hyperbolic arc cosine 

def arccosh(x):
	''' Compute the hyperbolic arc cosine of an AutoDiff object and its derivative.
	
	INPUTS
	======
	x: an AutoDiff object
	
	RETURNS
	=======
	A new AutoDiff object with calculated value and derivative.
	
	EXAMPLES
	========
	>>> x = AutoDiff(1.1, 2)
	>>> myAutoDiff = arccosh(x)
	>>> myAutoDiff.val
	0.4435682543851154
	>>> myAutoDiff.der
	(2/np.sqrt(1.1**2 - 1))
	>>> myAutoDiff.jacobian
	(1/np.sqrt(1.1**2 - 1))
	
	'''
	try:
		new_val = np.arccosh(x.val)
		# Derivative of arccosh is only defined when x > 1
		new_der = ((1)/np.sqrt(x.val**2 - 1))*x.der  # if x.val > 1 else None
		new_jacobian = ((1)/np.sqrt(x.val**2 - 1))*x.jacobian  # if x.val > 1 else None
		return AutoDiff(new_val, new_der, x.n, 0, new_jacobian)
	except AttributeError:
		try:
			return Dual(np.arccosh(x.Real), x.Dual/np.sqrt((x.Real**2)-1))		
		except AttributeError:
			try:
				return Dual(arccosh(x.Real), (x.Dual*((x.Real**2)-1)**-0.5))
			except AttributeError:
			# Constant
				return_val = np.arccosh(x)
				return return_val

# hyperbolic arc tangent 

def test_arccos_ad_results():
	# positive real numbers
	x = AutoDiff(0.5, 2)
	f = ef.arccos(x)
	assert f.val == np.array([[np.arccos(0.5)]])
	assert f.der == np.array([[-2/np.sqrt(1-0.5**2)]])
	assert f.jacobian == np.array([[-1/np.sqrt(1-0.5**2)]])

	# out of bounds - negative sqrt
	with pytest.warns(RuntimeWarning):
		y = AutoDiff(-2, 2)
		f = ef.arccos(y)
		assert np.isnan(f.val[0][0])
		assert np.isnan(f.der[0][0])
		assert np.isnan(f.jacobian[0][0])

	# out of bounds - divide by 0
	with pytest.warns(RuntimeWarning):
		y = AutoDiff(1, 2)
		f = ef.arccos(y)
		assert f.val == np.array([[np.arccos(1)]])
		assert np.isneginf(f.der[0][0])
		assert np.isneginf(f.jacobian[0][0])

	# zero
	z = AutoDiff(0, 2)
	f = ef.arccos(z)
	assert f.val == np.array([[np.arccos(0)]])
	assert f.der == np.array([[-2/np.sqrt(1-0**2)]])
	assert f.jacobian == np.array([[-1/np.sqrt(1-0**2)]]) 

def test_sqrt_constant_results():
	a = ef.sqrt(5)
	assert a == np.sqrt(5)
	b = ef.sqrt(0)
	assert b == np.sqrt(0)
	# Value undefined when x < 0
	with pytest.warns(RuntimeWarning):
		c = ef.sqrt(-5)
		assert np.isnan(c) 

def test_arcsin_ad_results():
	# positive real numbers
	x = Dual(0.5, 2)
	f = ef.arcsin(x)
	assert f.Real == np.array([[np.arcsin(0.5)]])
	assert f.Dual == np.array([[2/np.sqrt(1-0.5**2)]])

	# out of bounds - undefined sqrt
	with pytest.warns(RuntimeWarning):
		y = Dual(-2, 2)
		f = ef.arcsin(y)
		assert np.isnan(f.Real)
		assert np.isnan(f.Dual)

	# out of bounds - div by zero
	with pytest.warns(RuntimeWarning):
		y = Dual(1, 2)
		f = ef.arcsin(y)
		assert f.Real == np.array([[np.arcsin(1)]])
		assert np.isinf(f.Dual)

	# zero
	z = Dual(0, 2)
	f = ef.arcsin(z)
	assert f.Real == np.array([[0.0]])
	assert f.Dual == np.array([[2.0]]) 

def test_arccos_ad_results():
	# positive real numbers
	x = Dual(0.5, 2)
	f = ef.arccos(x)
	assert f.Real == np.array([[np.arccos(0.5)]])
	assert f.Dual == np.array([[-2/np.sqrt(1-0.5**2)]])

	# out of bounds - negative sqrt
	with pytest.warns(RuntimeWarning):
		y = Dual(-2, 2)
		f = ef.arccos(y)
		assert np.isnan(f.Real)
		assert np.isnan(f.Dual)

	# out of bounds - divide by 0
	with pytest.warns(RuntimeWarning):
		y = Dual(1, 2)
		f = ef.arccos(y)
		assert f.Real == np.array([[np.arccos(1)]])
		assert np.isneginf(f.Dual)

	# zero
	z = Dual(0, 2)
	f = ef.arccos(z)
	assert f.Real == np.array([[np.arccos(0)]])
	assert f.Dual == np.array([[-2/np.sqrt(1-0**2)]]) 

def test_sqrt_types():
	with pytest.raises(TypeError):
		ef.sqrt('x')
	with pytest.raises(TypeError):
		ef.sqrt("1234")


# ---------------LOGISTIC FUNC----------------# 

def __call__(self, shape, dtype=None):
        dtype = dtype or K.floatx()

        kernel_height, kernel_width, _, out_filters = shape
        fan_out = int(kernel_height * kernel_width * out_filters)
        return tf.random_normal(
            shape, mean=0.0, stddev=np.sqrt(2.0 / fan_out), dtype=dtype)


# Obtained from https://github.com/tensorflow/tpu/blob/master/models/official/efficientnet/efficientnet_model.py 

def __call__(self, shape, dtype=None):
        dtype = dtype or K.floatx()

        init_range = 1.0 / np.sqrt(shape[1])
        return tf.random_uniform(shape, -init_range, init_range, dtype=dtype)


# Obtained from https://github.com/tensorflow/tpu/blob/master/models/official/efficientnet/efficientnet_model.py 

def shared(shape, name):
#{{{
    """
    Create a shared object of a numpy array.
    """ 
    init=initializations.get('glorot_uniform');
    if len(shape) == 1:
        value = np.zeros(shape)  # bias are initialized with zeros
        return theano.shared(value=value.astype(theano.config.floatX), name=name)
    else:
        drange = np.sqrt(6. / (np.sum(shape)))
        value = drange * np.random.uniform(low=-1.0, high=1.0, size=shape)
        return init(shape=shape,name=name);
#}}} 

def lecun_uniform(shape, name=None, dim_ordering='th'):
    ''' Reference: LeCun 98, Efficient Backprop
        http://yann.lecun.com/exdb/publis/pdf/lecun-98b.pdf
    '''
    fan_in, fan_out = get_fans(shape, dim_ordering=dim_ordering)
    scale = np.sqrt(3. / fan_in)
    return uniform(shape, scale, name=name) 

def glorot_normal(shape, name=None, dim_ordering='th'):
    ''' Reference: Glorot & Bengio, AISTATS 2010
    '''
    fan_in, fan_out = get_fans(shape, dim_ordering=dim_ordering)
    s = np.sqrt(2. / (fan_in + fan_out))
    return normal(shape, s, name=name) 

def glorot_uniform(shape, name=None, dim_ordering='th'):
    fan_in, fan_out = get_fans(shape, dim_ordering=dim_ordering)
    s = np.sqrt(6. / (fan_in + fan_out))
    return uniform(shape, s, name=name) 

def he_uniform(shape, name=None, dim_ordering='th'):
    fan_in, fan_out = get_fans(shape, dim_ordering=dim_ordering)
    s = np.sqrt(6. / fan_in)
    return uniform(shape, s, name=name) 

def calc_l2_norm_of_vector_per_entry(grad_vec):
    return np.sqrt(np.sum(np.square(grad_vec))) / float(grad_vec.size) 

def trilateration_single_point(self, c, Dx, Dy):
        '''
        Given x at origin (0,0) and y at (0,c) the distances from a point
        at unknown location Dx, Dy to x, y, respectively, finds the position of the point.
        '''

        z = (c**2 - (Dy**2 - Dx**2)) / (2*c)
        t = np.sqrt(Dx**2 - z**2)

        return np.array([t,z]) 

def affine(self,
             num_out_channels,
             input_layer=None,
             num_channels_in=None,
             bias=0.0,
             stddev=None,
             activation='relu'):
    if input_layer is None:
      input_layer = self.top_layer
    if num_channels_in is None:
      num_channels_in = self.top_size
    name = 'affine' + str(self.counts['affine'])
    self.counts['affine'] += 1
    with tf.variable_scope(name):
      init_factor = 2. if activation == 'relu' else 1.
      stddev = stddev or np.sqrt(init_factor / num_channels_in)
      kernel = self.get_variable(
          'weights', [num_channels_in, num_out_channels],
          self.variable_dtype, self.dtype,
          initializer=tf.truncated_normal_initializer(stddev=stddev))
      biases = self.get_variable('biases', [num_out_channels],
                                 self.variable_dtype, self.dtype,
                                 initializer=tf.constant_initializer(bias))
      logits = tf.nn.xw_plus_b(input_layer, kernel, biases)
      if activation == 'relu':
        affine1 = tf.nn.relu(logits, name=name)
      elif activation == 'linear' or activation is None:
        affine1 = logits
      else:
        raise KeyError('Invalid activation type \'%s\'' % activation)
      self.top_layer = affine1
      self.top_size = num_out_channels
      return affine1 

def affine(self,
             num_out_channels,
             input_layer=None,
             num_channels_in=None,
             bias=0.0,
             stddev=None,
             activation='relu'):
    if input_layer is None:
      input_layer = self.top_layer
    if num_channels_in is None:
      num_channels_in = self.top_size
    name = 'affine' + str(self.counts['affine'])
    self.counts['affine'] += 1
    with tf.variable_scope(name):
      init_factor = 2. if activation == 'relu' else 1.
      stddev = stddev or np.sqrt(init_factor / num_channels_in)
      kernel = self.get_variable(
          'weights', [num_channels_in, num_out_channels],
          self.variable_dtype, self.dtype,
          initializer=tf.truncated_normal_initializer(stddev=stddev))
      biases = self.get_variable('biases', [num_out_channels],
                                 self.variable_dtype, self.dtype,
                                 initializer=tf.constant_initializer(bias))
      logits = tf.nn.xw_plus_b(input_layer, kernel, biases)
      if activation == 'relu':
        affine1 = tf.nn.relu(logits, name=name)
      elif activation == 'linear' or activation is None:
        affine1 = logits
      else:
        raise KeyError('Invalid activation type \'%s\'' % activation)
      self.top_layer = affine1
      self.top_size = num_out_channels
      return affine1 

def batchnorm(self, layer, inp):
        if not self.var:
            temp = (inp - layer.w['moving_mean'])
            temp /= (np.sqrt(layer.w['moving_variance']) + 1e-5)
            temp *= layer.w['gamma']
            return temp
        else:
            args = dict({
                'center' : False, 'scale' : True,
                'epsilon': 1e-5, 'scope' : self.scope,
                'updates_collections' : None,
                'is_training': layer.h['is_training'],
                'param_initializers': layer.w
                })
            return slim.batch_norm(inp, **args) 

def weight_variable(shape, name):
    #xavier initialisation
    in_dim = shape[0]
    xavier_stddev = 1. / tf.sqrt(in_dim / 2.)
    return tf.Variable(tf.random_normal(shape=shape, stddev=xavier_stddev), name=name) 

def train_test_validation(self, M, train_idx, test_idx, valid_idx, n_steps=100000, result_path='result/'):
        nonzero_user_idx = M.nonzero()[0]
        nonzero_item_idx = M.nonzero()[1]

        trainM = np.zeros(M.shape)
        trainM[nonzero_user_idx[train_idx], nonzero_item_idx[train_idx]] = M[nonzero_user_idx[train_idx], nonzero_item_idx[train_idx]]

        validM = np.zeros(M.shape)
        validM[nonzero_user_idx[valid_idx], nonzero_item_idx[valid_idx]] = M[nonzero_user_idx[valid_idx], nonzero_item_idx[valid_idx]]

        testM = np.zeros(M.shape)
        testM[nonzero_user_idx[test_idx], nonzero_item_idx[test_idx]] = M[nonzero_user_idx[test_idx], nonzero_item_idx[test_idx]]

        for i in range(self.num_user):
            if np.sum(trainM[i]) == 0:
                testM[i] = 0
                validM[i] = 0

        train_writer = tf.summary.FileWriter(
            result_path + '/train', graph=self.sess.graph)

        best_val_rmse = np.inf
        best_test_rmse = 0

        self.sess.run(tf.global_variables_initializer())
        for step in range(1, n_steps):
            feed_dict = {self.user: trainM, self.valid_rating:validM, self.test_rating:testM}

            _, mse, mae, valid_rmse, test_rmse,  summary_str = self.sess.run(
                [self.train_step, self.MSE, self.MAE, self.valid_RMSE, self.test_RMSE, self.summary_op], feed_dict=feed_dict)
            train_writer.add_summary(summary_str, step)
            print("Iter {0} Train RMSE:{1}, Valid RMSE:{2}, Test RMSE:{3}".format(step, np.sqrt(mse), valid_rmse, test_rmse))

            if best_val_rmse > valid_rmse:
                best_val_rmse = valid_rmse
                best_test_rmse = test_rmse

        self.saver.save(self.sess, result_path + "/model.ckpt")
        return best_test_rmse