Python numpy.cbrt() Examples

def from_xyz100(self, xyz):
        def f(t):
            delta = 6.0 / 29.0
            out = numpy.array(t, dtype=float)
            is_greater = out > delta ** 3
            out[is_greater] = 116 * numpy.cbrt(out[is_greater]) - 16
            out[~is_greater] = out[~is_greater] / (delta / 2) ** 3
            return out

        L = f(xyz[1] / self.whitepoint[1])

        x, y, z = xyz
        p = x + 15 * y + 3 * z
        u = 4 * x / p
        v = 9 * y / p

        wx, wy, wz = self.whitepoint
        q = wx + 15 * wy + 3 * wz
        un = 4 * wx / q
        vn = 9 * wy / q
        return numpy.array([L, 13 * L * (u - un), 13 * L * (v - vn)]) 

def _subsample_array(self, array):
        original_shape = array.shape
        spacing = self.spacing
        extent = tuple((o_s - 1) * s for o_s, s in zip(original_shape, spacing))
        dim_ratio = np.cbrt((array.nbytes / 1e6) / self.max_data_size)
        max_shape = tuple(int(o_s / dim_ratio) for o_s in original_shape)
        dowsnscale_factor = [max(o_s, m_s) / m_s for m_s, o_s in zip(max_shape, original_shape)]

        if any([d_f > 1 for d_f in dowsnscale_factor]):
            try:
                import scipy.ndimage as nd
                sub_array = nd.interpolation.zoom(array, zoom=[1 / d_f for d_f in dowsnscale_factor], order=0)
            except ImportError:
                sub_array = array[::int(np.ceil(dowsnscale_factor[0])),
                                  ::int(np.ceil(dowsnscale_factor[1])),
                                  ::int(np.ceil(dowsnscale_factor[2]))]
            self._sub_spacing = tuple(e / (s - 1) for e, s in zip(extent, sub_array.shape))
        else:
            sub_array = array
            self._sub_spacing = self.spacing
        return sub_array 

def symmetrized_stress(self, stress):
        #Make the matrix lower-diagonal
        for (i,j) in [(1,0),(2,1),(2,2)]:
            stress[i][j] = stress[i][j] + stress[j][i]
            stress[j][i] = 0
        #Normalize the matrix
        snm = self.struc.lattice.stress_normalization_matrix
        m2 = np.multiply(stress, snm)
        #Normalize the on-diagonal elements
        indices = self.struc.lattice.stress_indices
        if len(indices) == 2:
            total = 0
            for index in indices:
                total += stress[index]        
            value = total**0.5
            for index in inices:
                m2[index] = value
        elif len(indices) == 3:
            total = 0
            for index in indices:
                total += stress[index]        
            value = np.cbrt(total)
            for index in inices:
                m2[index] = value
        return m2 

def average_velocity(self):
        """
        This property calculates and returns the cube root of the
        mean cubed velocity in the turbine's rotor swept area (m/s).

        Returns:
            float: The average velocity across a rotor.

        Examples:
            To get the average velocity for a turbine:

            >>> avg_vel = floris.farm.turbines[0].average_velocity()
        """
        # remove all invalid numbers from interpolation
        data = self.velocities[np.where(~np.isnan(self.velocities))]
        avg_vel = np.cbrt(np.mean(data ** 3))
        if np.isnan(avg_vel):
            avg_vel = 0
        elif np.isinf(avg_vel):
            avg_vel = 0

        return avg_vel 

def symmetrized_stress(self, stress):
        #Make the matrix lower-diagonal
        for (i,j) in [(1,0),(2,1),(2,2)]:
            stress[i][j] = stress[i][j] + stress[j][i]
            stress[j][i] = 0
        #Normalize the matrix
        snm = self.struc.lattice.stress_normalization_matrix
        m2 = np.multiply(stress, snm)
        #Normalize the on-diagonal elements
        indices = self.struc.lattice.stress_indices
        if len(indices) == 2:
            total = 0
            for index in indices:
                total += stress[index]        
            value = np.sqrt(total)
            for index in inices:
                m2[index] = value
        elif len(indices) == 3:
            total = 0
            for index in indices:
                total += stress[index]        
            value = np.cbrt(total)
            for index in inices:
                m2[index] = value
        return m2 

def random_shear_matrix(width=1.0, unitary=False):
    """
    Generate a random symmetric shear matrix with Gaussian elements. If unitary
    is True, normalize to determinant 1

    Args:
        width: the width of the normal distribution to use when choosing values.
            Passed to np.random.normal
        unitary: whether or not to normalize the matrix to determinant 1
    
    Returns:
        a 3x3 numpy array of floats
    """
    mat = np.zeros([3,3])
    determinant = 0
    while determinant == 0:
        a, b, c = np.random.normal(scale=width), np.random.normal(scale=width), np.random.normal(scale=width)
        mat = np.array([[1,a,b],[a,1,c],[b,c,1]])
        determinant = np.linalg.det(mat)
    if unitary:
        new = mat / np.cbrt(np.linalg.det(mat))
        return new
    else: return mat 

def _no_extrapolation_hessian(internal_func, params_value, method):
    finite_diff = getattr(aux, method)
    hess = np.empty((len(params_value), len(params_value)))
    for i, val_1 in enumerate(params_value):
        h_1 = (1.0 + abs(val_1)) * np.cbrt(np.finfo(float).eps)
        for j, val_2 in enumerate(params_value):
            h_2 = (1.0 + abs(val_2)) * np.cbrt(np.finfo(float).eps)
            params_r = params_value.copy()
            params_r[j] += h_2
            # Calculate the first derivative w.r.t. var_1 at (params + h_2) with
            # the central method. This is not the right f_x0, but the real one
            # isn't needed for the central method.
            f_plus = finite_diff(internal_func, None, params_r, i, h_1)
            params_l = params_value.copy()
            params_l[j] -= h_2
            # Calculate the first derivative w.r.t. var_1 at (params - h_2) with
            # the central method. This is not the right f_x0, but the real one
            # isn't needed for the central method.
            f_minus = finite_diff(internal_func, None, params_l, i, h_1)
            f_diff = (f_plus - f_minus) / (2.0 * h_1 * h_2)
            hess[i, j] = f_diff
            hess[i, j] = f_diff
    return hess 

def test_cbrt_scalar(self):
        assert_almost_equal((np.cbrt(np.float32(-2.5)**3)), -2.5) 

def sample(center, radius, mode):
        """
        Samples a point according to the mode, the center and the radius.

       :param center, A list of three values, describing the x, y and z coordinate of the center of the sphere. Type: mathutils.Vector
       :param radius, The radius of the sphere. Type: float
       :param mode, Mode of sampling. SURFACE - sampling from the 2-sphere, INTERIOR - sampling from the 3-ball. Type: str
        """
        # Sample
        direction = np.random.normal(size=3)
        
        if np.count_nonzero(direction) == 0:  # Check no division by zero
            direction[0] = 1e-5

        # For normalization
        norm = np.sqrt(direction.dot(direction))

        # If sampling from the surface set magnitude to radius of the sphere
        if mode == "SURFACE":
            magnitude = radius
        # If sampling from the interior set it to scaled radius
        elif mode == "INTERIOR":
            magnitude = radius * np.cbrt(np.random.unform())
        else:
            raise Exception("Unknown sampling mode: " + mode)
        
        # Normalize
        sampled_point = list(map(lambda x: magnitude*x/norm, direction))
        
        # Add center
        location = mathutils.Vector(np.array(sampled_point) + center)

        return location 

def test_cbrt(self):
        x = np.array([1., 2., -3., np.inf, -np.inf])
        assert_almost_equal(np.cbrt(x**3), x)

        assert_(np.isnan(np.cbrt(np.nan)))
        assert_equal(np.cbrt(np.inf), np.inf)
        assert_equal(np.cbrt(-np.inf), -np.inf) 

def test_cbrt_scalar(self):
        assert_almost_equal((np.cbrt(np.float32(-2.5)**3)), -2.5) 

def test_norm(self):
        """
        Checking that the norm function works properly for various powers.
        Checked both by
        """
        self.assertAlmostEqual(np.sqrt(41) ,mb.norm(5,4,2))
        self.assertAlmostEqual(np.cbrt(72.248104), mb.norm(2.5,3.84,3)) 

def test_cbrt(self):
        x = np.array([1., 2., -3., np.inf, -np.inf])
        assert_almost_equal(np.cbrt(x**3), x)

        assert_(np.isnan(np.cbrt(np.nan)))
        assert_equal(np.cbrt(np.inf), np.inf)
        assert_equal(np.cbrt(-np.inf), -np.inf) 

def test_cbrt_scalar(self):
        assert_almost_equal((np.cbrt(np.float32(-2.5)**3)), -2.5) 

def test_cbrt_scalar(self):
        assert_almost_equal((np.cbrt(np.float32(-2.5)**3)), -2.5) 

def test_cbrt(self):
        x = np.array([1., 2., -3., np.inf, -np.inf])
        assert_almost_equal(np.cbrt(x**3), x)

        assert_(np.isnan(np.cbrt(np.nan)))
        assert_equal(np.cbrt(np.inf), np.inf)
        assert_equal(np.cbrt(-np.inf), -np.inf) 

def test_cbrt(self):
        x = np.array([1., 2., -3., np.inf, -np.inf])
        assert_almost_equal(np.cbrt(x**3), x)

        assert_(np.isnan(np.cbrt(np.nan)))
        assert_equal(np.cbrt(np.inf), np.inf)
        assert_equal(np.cbrt(-np.inf), -np.inf) 

def _build_placeholder(self):

        # Create placeholders for the input to the siamese network
        self.anchor_input = tf.placeholder(dtype=tf.float32, shape=[None, int(np.cbrt(self.config.input_dim)),
                                                                    int(np.cbrt(self.config.input_dim)),
                                                                    int(np.cbrt(self.config.input_dim)), 1],
                                           name='X_reference')

        self.positive_input = tf.placeholder(dtype=tf.float32, shape=[None, int(np.cbrt(self.config.input_dim)),
                                                                      int(np.cbrt(self.config.input_dim)),
                                                                      int(np.cbrt(self.config.input_dim)), 1],
                                             name='X_positive')

        # Global step for optimization
        self.global_step = tf.Variable(0, trainable=False) 

def test_cbrt_scalar(self):
        assert_almost_equal((np.cbrt(np.float32(-2.5)**3)), -2.5) 

def test_cbrt(self):
        x = np.array([1., 2., -3., np.inf, -np.inf])
        assert_almost_equal(np.cbrt(x**3), x)

        assert_(np.isnan(np.cbrt(np.nan)))
        assert_equal(np.cbrt(np.inf), np.inf)
        assert_equal(np.cbrt(-np.inf), -np.inf) 

def chirp_mass_and_primary_mass_to_mass_ratio(chirp_mass, mass_1):
    """
    Convert chirp mass and mass ratio of a binary to its total mass.

    Rearranging the relation for chirp mass (as a function of mass_1 and
    mass_2) and q = mass_2 / mass_1, it can be shown that

        (chirp_mass/mass_1)^5 = q^3 / (1 + q)

    Solving for q, we find the releation expressed in python below for q.

    Parameters
    ----------
    chirp_mass: float
        Chirp mass of the binary
    mass_1: float
        The primary mass

    Return
    ------
    mass_ratio: float
        Mass ratio (mass_2/mass_1) of the binary
    """
    a = (chirp_mass / mass_1) ** 5
    t0 = np.cbrt(9 * a + np.sqrt(3) * np.sqrt(27 * a ** 2 - 4 * a ** 3))
    t1 = np.cbrt(2) * 3 ** (2 / 3)
    t2 = np.cbrt(2 / 3) * a
    return t2 / t0 + t0 / t1 

def test_vtkvol_pane_from_np_array(document, comm):
    # Test empty initialisation
    pane = VTKVolume()
    model = pane.get_root(document, comm=comm)

    pane.object = np.ones((10,10,10))
    from operator import eq
    # Create pane
    assert isinstance(model, VTKVolumePlot)
    assert pane._models[model.ref['id']][0] is model
    assert np.all(np.frombuffer(base64.b64decode(model.data['buffer'].encode())) == 1)
    assert all([eq(getattr(pane, k), getattr(model, k))
                for k in ['slice_i', 'slice_j', 'slice_k']])

    # Test update data
    pane.object = 2*np.ones((10,10,10))
    assert np.all(np.frombuffer(base64.b64decode(model.data['buffer'].encode())) == 2)

    # Test size limitation of date sent
    pane.max_data_size = 0.1 # limit data size to 0.1MB
    # with uint8
    data = (255*np.random.rand(50,50,50)).astype(np.uint8)
    assert data.nbytes/1e6 > 0.1
    pane.object = data
    data_model = np.frombuffer(base64.b64decode(model.data['buffer'].encode()))
    assert data_model.nbytes/1e6 <= 0.1
    # with float64
    data = np.random.rand(50,50,50)
    assert data.nbytes/1e6 > 0.1
    pane.object = data
    data_model = np.frombuffer(base64.b64decode(model.data['buffer'].encode()), dtype=np.float64)
    assert data_model.nbytes/1e6 <= 0.1

    # Test conversion of the slice_i number with subsample array
    param = pane._process_property_change({'slice_i': (np.cbrt(data_model.size)-1)//2})
    assert param == {'slice_i': (50-1)//2}

    # Cleanup
    pane._cleanup(model)
    assert pane._models == {} 

def from_blocks_numpy(a, output_shape):
    """Combine 4D array of non-overlapping blocks `a` into 3D array of shape
    `output_shape`.

    For the reverse of this function, see `to_blocks_numpy`.

    Parameters
    ----------
    a: array-like, 4D array of blocks with shape (N, *block_shape), where N is
        the number of blocks.
    output_shape: tuple of len 3, shape of the combined array.

    Returns
    -------
    Rank 3 array with shape `output_shape`.
    """
    a = np.asarray(a)

    if a.ndim != 4:
        raise ValueError("This function only works for 4D arrays.")
    if len(output_shape) != 3:
        raise ValueError("output_shape must have three values.")

    n_blocks = a.shape[0]
    block_shape = a.shape[1:]
    ncbrt = np.cbrt(n_blocks).round(6)
    if not ncbrt.is_integer():
        raise ValueError("Cubed root of number of blocks is not an integer")
    ncbrt = int(ncbrt)
    intershape = (ncbrt, ncbrt, ncbrt, *block_shape)

    return a.reshape(intershape).transpose((0, 3, 1, 4, 2, 5)).reshape(output_shape) 

def from_blocks(x, output_shape):
    """Combine 4D array of non-overlapping sub-volumes `x` into 3D tensor of
    shape `output_shape`.

    For the reverse of this function, see `to_blocks`.

    Parameters
    ----------
    x: tensor, 4D tensor with shape `(N, *block_shape)`, where `N` is the number
        of sub-volumes.
    output_shape: tuple of length 3, shape of resulting volumes.

    Returns
    -------
    Tensor with shape `output_shape`.
    """
    x = tf.convert_to_tensor(x)
    n_blocks = x.shape[0]
    block_shape = x.shape[1:]
    ncbrt = np.cbrt(n_blocks).round(6)
    if not ncbrt.is_integer():
        raise ValueError("Cubed root of number of blocks is not an integer")
    ncbrt = int(ncbrt)
    intershape = (ncbrt, ncbrt, ncbrt, *block_shape)
    perm = (0, 3, 1, 4, 2, 5)  # 3D only

    return tf.reshape(
        tf.transpose(tf.reshape(x, shape=intershape), perm=perm), shape=output_shape
    )


# Below this line, we implement methods similar to those above but using Numpy.
# This is particularly useful when we use models to predict, because it is
# usually more pleasant to predict on Numpy arrays. 

def random_matrix(width=1.0, unitary=False):
    mat = np.zeros([3,3])
    for x in range(3):
        for y in range(3):
            mat[x][y] = normal(scale=width)
    if unitary:
        new = mat / cbrt(det(mat))
        return new
    else: return mat 

def strain_matrix(a, b, c):
    a = mat_a(a)
    b = mat_b(b)
    c = mat_c(c)
    raw = a+b+c-2*identity
    return raw/cbrt(det(raw)) 

def from_xyz100(self, xyz100):
        X, Y, Z = xyz100
        s = numpy.sum(xyz100, axis=0)

        # Avoid division by s, could be 0.
        YKs2 = (
            4.4934 * Y * X ** 2
            + 4.3034 * Y ** 3
            - 4.276 * X * Y ** 2
            - 1.3744 * X * Y * s
            - 2.5643 * Y ** 2 * s
            + 1.8103 * Y * s ** 2
        )
        Y0 = numpy.zeros_like(Y)
        idx = numpy.abs(s) > 1.0e-15
        Y0[idx] = YKs2[idx] / s[idx] ** 2

        #  L' is L in original article
        L_prime = 5.9 * (numpy.cbrt(Y0) - 2 / 3 + 0.042 * numpy.cbrt(Y0 - 30))

        C = L_prime / (5.9 * (numpy.cbrt(Y0) - 2 / 3))
        R, G, B = dot(self.M, xyz100)

        a = -13.7 * numpy.cbrt(R) + 17.7 * numpy.cbrt(G) - 4 * numpy.cbrt(B)
        b = 1.7 * numpy.cbrt(R) + 8 * numpy.cbrt(G) - 9.7 * numpy.cbrt(B)

        L = (L_prime - 14.3993) / numpy.sqrt(2)
        g = C * a
        j = C * b

        return numpy.array([L, g, j]) 

def from_xyz100(self, xyz):
        def f(t):
            delta = 6 / 29
            out = numpy.array(t, dtype=float)
            is_greater = out > delta ** 3
            out[is_greater] = 116 * numpy.cbrt(out[is_greater]) - 16
            out[~is_greater] = out[~is_greater] / (delta / 2) ** 3
            return out

        fx, fy, fz = f((xyz.T / self.whitepoint).T)
        return numpy.array([fy, 125 / 29 * (fx - fy), 50 / 29 * (fy - fz)]) 

def test_speed(N=2):
    numpy.random.seed(1)
    osa = colorio.OsaUcs()
    cielab = colorio.CIELAB()
    # cam16 = colorio.CAM16(0.69, 20, L_A=64 / numpy.pi / 5)
    ciecam02 = colorio.CIECAM02(0.69, 20, L_A=64 / numpy.pi / 5)

    # This close probably means that another figure hasn't been properly closed.
    import matplotlib.pyplot as plt

    plt.close()

    perfplot.show(
        # Don't use numpy.random.rand(3, n) to avoid the CIECAM breakdown
        setup=lambda n: numpy.outer(numpy.random.rand(3), numpy.ones(n)) * 10,
        equality_check=None,
        kernels=[
            osa.to_xyz100,
            cielab.to_xyz100,
            # cam16.to_xyz100,
            lambda Jsh: ciecam02.to_xyz100(Jsh, "Jsh"),
            numpy.cbrt,
        ],
        labels=["OSA-UCS", "CIELAB", "CIECAM02", "cbrt"],
        n_range=[2 ** n for n in range(N)],
        logx=True,
        logy=True,
        # relative_to=3
    )
    # import tikzplotlib as tpl
    # tpl.save("out.tex") 

def _split_epsilon(self, dims, total_iters, rho=0.225):
        """Split epsilon between sum perturbation and count perturbation, as proposed by Su et al.

        Parameters
        ----------
        dims : int
            Number of dimensions to split `epsilon` across.

        total_iters : int
            Total number of iterations to split `epsilon` across.

        rho : float, default: 0.225
            Coordinate normalisation factor.

        Returns
        -------
        epsilon_0 : float
            The epsilon value for satisfying differential privacy on the count of a cluster.

        epsilon_i : float
            The epsilon value for satisfying differential privacy on each dimension of the center of a cluster.

        """
        epsilon_i = 1
        epsilon_0 = np.cbrt(4 * dims * rho ** 2)

        normaliser = self.epsilon / total_iters / (epsilon_i * dims + epsilon_0)

        return epsilon_i * normaliser, epsilon_0 * normaliser