Ellipsoid test objective function
183271
def elli(self, x, rot=0, xoffset=0, cond=1e6, actuator_noise=0.0, both=False):
        """"""
        if not isscalar(x[0]):  # parallel evaluation
            return [self.elli(xi, rot) for xi in x]  # could save 20% overall
        if rot:
            x = rotate(x)
        N = len(x)
        if actuator_noise:
            x = x + actuator_noise * np.random.randn(N)

        ftrue = sum(cond**(np.arange(N) / (N - 1.)) * (x + xoffset)**2)

        alpha = 0.49 + 1. / N
        beta = 1
        felli = np.random.rand(1)[0]**beta * ftrue * \
                max(1, (10.**9 / (ftrue + 1e-99))**(alpha * np.random.rand(1)[0]))
        # felli = ftrue + 1*np.random.randn(1)[0] / (1e-30 +
        #                                           np.abs(np.random.randn(1)[0]))**0
        if both:
            return (felli, ftrue)
        else:
            # return felli  # possibly noisy value
            return ftrue  # + np.random.randn()
352938
def xstep(self):
        
        """

        ngsit = 0
        gsrrs = np.inf
        YU = self.Y - self.U
        SYU = self.S + YU[..., -1]
        YU[..., -1] = 0.0
        ATYU = self.cnst_AT(YU)
        while gsrrs > self.opt['GSTol'] and ngsit < self.opt['MaxGSIter']:
            self.X = self.GaussSeidelStep(
                SYU, self.X, ATYU, 1.0, self.lcw, 1.0)
            gsrrs = sl.rrs(
                self.cnst_AT(self.cnst_A(self.X)),
                self.cnst_AT(self.cnst_c() - self.cnst_B(self.Y) - self.U)
            )
            ngsit += 1

        self.xs = (ngsit, gsrrs)
1199606
def cylinders_from_ellipses(ellipses):
    """"""
    ellipses = np.asarray(ellipses)
    ellipsoids = np.zeros((ellipses.shape[0], 10))
    ellipsoids[:, [0, 1, 2, 4, 5, 7]] = ellipses
    ellipsoids[:, 3] = 100000.0

    return ellipsoids
650799
def _distance_covariance_sqr_fast_generic(
        x, y, unbiased=False):  # pylint:disable=too-many-locals
    # This function has many locals so it can be compared
    # with the original algorithm.
    """"""
    x = np.asarray(x)
    y = np.asarray(y)

    x = np.ravel(x)
    y = np.ravel(y)

    n = x.shape[0]
    assert n > 3
    assert n == y.shape[0]
    temp = range(n)

    # Step 1
    ix0 = np.argsort(x)
    vx = x[ix0]

    ix = np.zeros(n, dtype=int)
    ix[ix0] = temp

    iy0 = np.argsort(y)
    vy = y[iy0]

    iy = np.zeros(n, dtype=int)
    iy[iy0] = temp

    # Step 2
    sx = np.cumsum(vx)
    sy = np.cumsum(vy)

    # Step 3
    alpha_x = ix
    alpha_y = iy

    beta_x = sx[ix] - vx[ix]
    beta_y = sy[iy] - vy[iy]

    # Step 4
    x_dot = np.sum(x)
    y_dot = np.sum(y)

    # Step 5
    a_i_dot = x_dot + (2 * alpha_x - n) * x - 2 * beta_x
    b_i_dot = y_dot + (2 * alpha_y - n) * y - 2 * beta_y

    sum_ab = np.sum(a_i_dot * b_i_dot)

    # Step 6
    a_dot_dot = 2 * np.sum(alpha_x * x) - 2 * np.sum(beta_x)
    b_dot_dot = 2 * np.sum(alpha_y * y) - 2 * np.sum(beta_y)

    # Step 7
    gamma_1 = _partial_sum_2d(x, y, np.ones(n, dtype=x.dtype))
    gamma_x = _partial_sum_2d(x, y, x)
    gamma_y = _partial_sum_2d(x, y, y)
    gamma_xy = _partial_sum_2d(x, y, x * y)

    # Step 8
    aijbij = np.sum(x * y * gamma_1 + gamma_xy - x * gamma_y - y * gamma_x)

    if unbiased:
        d3 = (n - 3)
        d2 = (n - 2)
        d1 = (n - 1)
    else:
        d3 = d2 = d1 = n

    # Step 9
    d_cov = (aijbij / n / d3 - 2 * sum_ab / n / d2 / d3 +
             a_dot_dot / n * b_dot_dot / d1 / d2 / d3)

    return d_cov
58625
function(x, y, elem) {
        Garnish.hitTest._$elem = $(elem);
        Garnish.hitTest._offset = Garnish.hitTest._$elem.offset();
        Garnish.hitTest._x1 = Garnish.hitTest._offset.left;
        Garnish.hitTest._y1 = Garnish.hitTest._offset.top;
        Garnish.hitTest._x2 = Garnish.hitTest._x1 + Garnish.hitTest._$elem.outerWidth();
        Garnish.hitTest._y2 = Garnish.hitTest._y1 + Garnish.hitTest._$elem.outerHeight();

        return (x >= Garnish.hitTest._x1 && x < Garnish.hitTest._x2 && y >= Garnish.hitTest._y1 && y < Garnish.hitTest._y2);
    }
35901
def set_differentiable_objective(self):
    """"""
    # Checking if graphs are already created
    if self.vector_g is not None:
      return

    # Computing the scalar term
    bias_sum = 0
    for i in range(0, self.nn_params.num_hidden_layers):
      bias_sum = bias_sum + tf.reduce_sum(
          tf.multiply(self.nn_params.biases[i], self.lambda_pos[i + 1]))
    lu_sum = 0
    for i in range(0, self.nn_params.num_hidden_layers + 1):
      lu_sum = lu_sum + tf.reduce_sum(
          tf.multiply(tf.multiply(self.lower[i], self.upper[i]),
                      self.lambda_lu[i]))

    self.scalar_f = -bias_sum - lu_sum + self.final_constant

    # Computing the vector term
    g_rows = []
    for i in range(0, self.nn_params.num_hidden_layers):
      if i > 0:
        current_row = (self.lambda_neg[i] + self.lambda_pos[i] -
                       self.nn_params.forward_pass(self.lambda_pos[i+1],
                                                   i, is_transpose=True) +
                       tf.multiply(self.lower[i]+self.upper[i],
                                   self.lambda_lu[i]) +
                       tf.multiply(self.lambda_quad[i],
                                   self.nn_params.biases[i-1]))
      else:
        current_row = (-self.nn_params.forward_pass(self.lambda_pos[i+1],
                                                    i, is_transpose=True)
                       + tf.multiply(self.lower[i]+self.upper[i],
                                     self.lambda_lu[i]))
      g_rows.append(current_row)

    # Term for final linear term
    g_rows.append((self.lambda_pos[self.nn_params.num_hidden_layers] +
                   self.lambda_neg[self.nn_params.num_hidden_layers] +
                   self.final_linear +
                   tf.multiply((self.lower[self.nn_params.num_hidden_layers]+
                                self.upper[self.nn_params.num_hidden_layers]),
                               self.lambda_lu[self.nn_params.num_hidden_layers])
                   + tf.multiply(
                       self.lambda_quad[self.nn_params.num_hidden_layers],
                       self.nn_params.biases[
                           self.nn_params.num_hidden_layers-1])))
    self.vector_g = tf.concat(g_rows, axis=0)
    self.unconstrained_objective = self.scalar_f + 0.5 * self.nu
430190
def reconstruct_from_shape(self, shape, optimize=False):
        
        axes_lengths = list(self.elementary_axes_lengths)
        if self.ellipsis_positions != (math.inf, math.inf):
            if len(shape) < len(self.input_composite_axes) - 1:
                raise EinopsError('Expected at least {} dimensions, got {}'.format(
                    len(self.input_composite_axes) - 1, len(shape)))
        else:
            if len(shape) != len(self.input_composite_axes):
                raise EinopsError('Expected {} dimensions, got {}'.format(len(self.input_composite_axes), len(shape)))
        for input_axis, (known_axes, unknown_axes) in enumerate(self.input_composite_axes):
            before_ellipsis = input_axis
            after_ellipsis = input_axis + len(shape) - len(self.input_composite_axes)
            if input_axis == self.ellipsis_positions[0]:
                assert len(known_axes) == 0 and len(unknown_axes) == 1
                unknown_axis, = unknown_axes
                ellipsis_shape = shape[before_ellipsis:after_ellipsis + 1]
                if any(d is None for d in ellipsis_shape):
                    raise EinopsError("Couldn't infer shape for one or more axes represented by ellipsis")
                axes_lengths[unknown_axis] = _product(ellipsis_shape)
            else:
                if input_axis < self.ellipsis_positions[0]:
                    length = shape[before_ellipsis]
                else:
                    length = shape[after_ellipsis]
                known_product = 1
                for axis in known_axes:
                    known_product *= axes_lengths[axis]

                if len(unknown_axes) == 0:
                    if isinstance(length, int) and isinstance(known_product, int) and length != known_product:
                        raise EinopsError('Shape mismatch, {} != {}'.format(length, known_product))
                else:
                    if isinstance(length, int) and isinstance(known_product, int) and length % known_product != 0:
                        raise EinopsError("Shape mismatch, can't divide axis of length {} in chunks of {}".format(
                            length, known_product))

                    unknown_axis, = unknown_axes
                    axes_lengths[unknown_axis] = length // known_product

        init_shapes = axes_lengths
        reduced_axes_lengths = [dim for i, dim in enumerate(axes_lengths) if i not in self.reduced_elementary_axes]
        final_shapes = []
        for output_axis, grouping in enumerate(self.output_composite_axes):
            if output_axis == self.ellipsis_positions[1]:
                final_shapes.extend(ellipsis_shape)
            else:
                lengths = [reduced_axes_lengths[elementary_axis] for elementary_axis in grouping]
                if any(l is None for l in lengths):
                    final_shapes.append(None)
                else:
                    final_shapes.append(_product(lengths))
        reduced_axes = self.reduced_elementary_axes
        axes_reordering = self.final_axes_grouping_flat
        if optimize:
            return _optimize_transformation(init_shapes, reduced_axes, axes_reordering, final_shapes)
        else:
            return init_shapes, reduced_axes, axes_reordering, final_shapes
1176092
def Ellipsoid(pos=(0, 0, 0), axis1=(1, 0, 0), axis2=(0, 2, 0), axis3=(0, 0, 3),
              c="c", alpha=1, res=24):
    
    elliSource = vtk.vtkSphereSource()
    elliSource.SetThetaResolution(res)
    elliSource.SetPhiResolution(res)
    elliSource.Update()
    l1 = np.linalg.norm(axis1)
    l2 = np.linalg.norm(axis2)
    l3 = np.linalg.norm(axis3)
    axis1 = np.array(axis1) / l1
    axis2 = np.array(axis2) / l2
    axis3 = np.array(axis3) / l3
    angle = np.arcsin(np.dot(axis1, axis2))
    theta = np.arccos(axis3[2])
    phi = np.arctan2(axis3[1], axis3[0])

    t = vtk.vtkTransform()
    t.PostMultiply()
    t.Scale(l1, l2, l3)
    t.RotateX(angle * 57.3)
    t.RotateY(theta * 57.3)
    t.RotateZ(phi * 57.3)
    tf = vtk.vtkTransformPolyDataFilter()
    tf.SetInputData(elliSource.GetOutput())
    tf.SetTransform(t)
    tf.Update()
    pd = tf.GetOutput()

    actor = Actor(pd, c=c, alpha=alpha)
    actor.GetProperty().BackfaceCullingOn()
    actor.GetProperty().SetInterpolationToPhong()
    actor.SetPosition(pos)
    actor.base = -np.array(axis1) / 2 + pos
    actor.top = np.array(axis1) / 2 + pos
    settings.collectable_actors.append(actor)
    return actor
265867
def is_functional(cls):
        
        if not cls._tested:
            cls._tested = True
            np.random.seed(SEED)
            test_problem_dimension = 10
            mat = np.random.randn(test_problem_dimension, test_problem_dimension)
            posmat = mat.dot(mat.T)
            posvar = cvxpy.Variable(test_problem_dimension, test_problem_dimension)
            prob = cvxpy.Problem(cvxpy.Minimize((cvxpy.trace(posmat * posvar)
                                                 + cvxpy.norm(posvar))),
                                 [posvar >> 0, cvxpy.trace(posvar) >= 1.])

            try:
                prob.solve(SOLVER)
                cls._functional = True
            except cvxpy.SolverError:  # pragma no coverage
                _log.warning("No convex SDP solver found. You will not be able to solve"
                             " tomography problems with matrix positivity constraints.")
        return cls._functional
1249036
def cart2pol(x, y):
    """"""
    theta = np.arctan2(y, x)
    rho = np.hypot(x, y)
    return theta, rho
1231775
def poisson_objective(X, m, w):
    
    clusters, cells = w.shape
    genes = X.shape[0]
    #m = m.reshape((X.shape[0], w.shape[0]))
    d = m.dot(w)+eps
    #temp = X/d
    #w_sum = w.sum(1)
    #w2 = w.dot(temp.T)
    #deriv = w_sum - w2.T
    return np.sum(d - X*np.log(d))/genes
1012078
def _modified_shepp_logan_ellipsoids(ellipsoids):
    
    intensities = [1.0, -0.8, -0.2, -0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]

    # Add minimal numbers to ensure that the result is nowhere negative.
    # This is needed due to numerical issues.
    intensities[2] += 5e-17
    intensities[3] += 5e-17

    assert len(ellipsoids) == len(intensities)

    for ellipsoid, intensity in zip(ellipsoids, intensities):
        ellipsoid[0] = intensity
659844
def hinit(func, x, t, pos_neg, f0, iord, hmax, rtol, atol, args):
    
    sk = atol + rtol * np.fabs(x)
    dnf = np.sum(np.square(f0 / sk), axis=0)
    dny = np.sum(np.square(x / sk), axis=0)

    h = np.sqrt(dny / dnf) * 0.01

    h = np.min([h, np.fabs(hmax)])
    h = custom_sign(h, pos_neg)

    # perform an explicit Euler step
    xx1 = x + h * f0
    f1 = np.array(func(xx1, t[0] + h, *args))

    # estimate the second derivative of the solution
    der2 = np.sum(np.square((f1 - f0) / sk), axis=0)
    der2 = np.sqrt(der2) / h

    # step size is computed such that h ** iord * max_d(norm(f0), norm(der2)) = 0.01
    der12 = np.max([np.fabs(der2), np.sqrt(dnf)])
    h1 = np.power(0.01 / der12, 1.0 / iord)

    h = np.min([100.0 * np.fabs(h), np.min([np.fabs(h1), np.fabs(hmax)])])

    return custom_sign(h, pos_neg), f0, f1, xx1
670567
def _ellipsoid_bootstrap_expand(args):
    """"""

    # Unzipping.
    points, pointvol = args
    rstate = np.random

    # Resampling.
    npoints, ndim = points.shape
    idxs = rstate.randint(npoints, size=npoints)  # resample
    idx_in = np.unique(idxs)  # selected objects
    sel = np.ones(npoints, dtype='bool')
    sel[idx_in] = False
    idx_out = np.arange(npoints)[sel]  # "missing" objects
    if len(idx_out) < 2:  # edge case
        idx_out = np.append(idx_out, [0, 1])
    points_in, points_out = points[idx_in], points[idx_out]

    # Compute bounding ellipsoid.
    ell = bounding_ellipsoid(points_in, pointvol=pointvol)

    # Compute normalized distances to missing points.
    dists = [ell.distance(p) for p in points_out]

    # Compute expansion factor.
    expand = max(1., max(dists))

    return expand
339175
def add_ellipsoid(self, x0, radii, lcar=None, with_volume=True, holes=None):
        
        if holes is None:
            holes = []

        if holes:
            assert with_volume

        # Add points.
        p = [
            self.add_point(x0, lcar=lcar),
            self.add_point([x0[0] + radii[0], x0[1], x0[2]], lcar=lcar),
            self.add_point([x0[0], x0[1] + radii[1], x0[2]], lcar=lcar),
            self.add_point([x0[0], x0[1], x0[2] + radii[2]], lcar=lcar),
            self.add_point([x0[0] - radii[0], x0[1], x0[2]], lcar=lcar),
            self.add_point([x0[0], x0[1] - radii[1], x0[2]], lcar=lcar),
            self.add_point([x0[0], x0[1], x0[2] - radii[2]], lcar=lcar),
        ]
        # Add skeleton.
        # Alternative for circles:
        # `self.add_circle_arc(a, b, c)`
        c = [
            self.add_ellipse_arc(p[1], p[0], p[6], p[6]),
            self.add_ellipse_arc(p[6], p[0], p[4], p[4]),
            self.add_ellipse_arc(p[4], p[0], p[3], p[3]),
            self.add_ellipse_arc(p[3], p[0], p[1], p[1]),
            self.add_ellipse_arc(p[1], p[0], p[2], p[2]),
            self.add_ellipse_arc(p[2], p[0], p[4], p[4]),
            self.add_ellipse_arc(p[4], p[0], p[5], p[5]),
            self.add_ellipse_arc(p[5], p[0], p[1], p[1]),
            self.add_ellipse_arc(p[6], p[0], p[2], p[2]),
            self.add_ellipse_arc(p[2], p[0], p[3], p[3]),
            self.add_ellipse_arc(p[3], p[0], p[5], p[5]),
            self.add_ellipse_arc(p[5], p[0], p[6], p[6]),
        ]
        # Add surfaces (1/8th of the ball surface).
        ll = [
            # one half
            self.add_line_loop([c[4], c[9], c[3]]),
            self.add_line_loop([c[8], -c[4], c[0]]),
            self.add_line_loop([-c[9], c[5], c[2]]),
            self.add_line_loop([-c[5], -c[8], c[1]]),
            # the other half
            self.add_line_loop([c[7], -c[3], c[10]]),
            self.add_line_loop([c[11], -c[7], -c[0]]),
            self.add_line_loop([-c[10], -c[2], c[6]]),
            self.add_line_loop([-c[1], -c[6], -c[11]]),
        ]
        # Create a surface for each line loop.
        s = [self.add_surface(l) for l in ll]

        # Combine the surfaces to avoid seams
        if self._gmsh_major() == 3:
            s = [self.add_compound_surface(s[:4]), self.add_compound_surface(s[4:])]
        else:
            assert self._gmsh_major() == 4
            # <https://gitlab.onelab.info/gmsh/gmsh/issues/507>
            self.add_raw_code(
                "Compound Surface{{{}}};".format(",".join([surf.id for surf in s[:4]]))
            )
            self.add_raw_code(
                "Compound Surface{{{}}};".format(",".join([surf.id for surf in s[4:]]))
            )

        # Create the surface loop.
        surface_loop = self.add_surface_loop(s)
        # if holes:
        #     # Create an array of surface loops; the first entry is the outer
        #     # surface loop, the following ones are holes.
        #     surface_loop = self.add_array([surface_loop] + holes)
        # Create volume.
        volume = self.add_volume(surface_loop, holes) if with_volume else None

        class Ellipsoid(object):
            dimension = 3

            def __init__(self, x0, radii, surface_loop, volume, lcar=None):
                self.x0 = x0
                self.lcar = lcar
                self.radii = radii
                self.surface_loop = surface_loop
                self.volume = volume
                return

        return Ellipsoid(x0, radii, surface_loop, volume, lcar=lcar)
1089007
def _derY(self,x,y):
        '''
        
        '''
        x_pos, y_pos = self.findSector(x,y)
        alpha, beta = self.findCoords(x,y,x_pos,y_pos)

        # Get four corners data for each point
        xA = self.x_values[x_pos,y_pos]
        xB = self.x_values[x_pos+1,y_pos]
        xC = self.x_values[x_pos,y_pos+1]
        xD = self.x_values[x_pos+1,y_pos+1]
        yA = self.y_values[x_pos,y_pos]
        yB = self.y_values[x_pos+1,y_pos]
        yC = self.y_values[x_pos,y_pos+1]
        yD = self.y_values[x_pos+1,y_pos+1]
        fA = self.f_values[x_pos,y_pos]
        fB = self.f_values[x_pos+1,y_pos]
        fC = self.f_values[x_pos,y_pos+1]
        fD = self.f_values[x_pos+1,y_pos+1]

        # Calculate components of the alpha,beta --> x,y delta translation matrix
        alpha_x = (1-beta)*(xB-xA) + beta*(xD-xC)
        alpha_y = (1-beta)*(yB-yA) + beta*(yD-yC)
        beta_x  = (1-alpha)*(xC-xA) + alpha*(xD-xB)
        beta_y  = (1-alpha)*(yC-yA) + alpha*(yD-yB)

        # Invert the delta translation matrix into x,y --> alpha,beta
        det = alpha_x*beta_y - beta_x*alpha_y
        y_alpha = -beta_x/det
        y_beta  = alpha_x/det

        # Calculate the derivative of f w.r.t. alpha and beta
        dfda = (1-beta)*(fB-fA) + beta*(fD-fC)
        dfdb = (1-alpha)*(fC-fA) + alpha*(fD-fB)

        # Calculate the derivative with respect to x (and return it)
        dfdy = y_alpha*dfda + y_beta*dfdb
        return dfdy
830793
def ellipse(center,covariance_matrix,level=1, n=1000):
    """"""
    # singular value decomposition
    U, s, rotation_matrix = N.linalg.svd(covariance_matrix)
    # semi-axes (largest first)

    saxes = N.sqrt(s)*level ## If the _area_ of a 2s ellipse is twice that of a 1s ellipse
    # If the _axes_ are supposed to be twice as long, then it should be N.sqrt(s)*width

    u = N.linspace(0, 2*N.pi, n)
    data = N.column_stack((saxes[0]*N.cos(u), saxes[1]*N.sin(u)))
    # rotate data
    return N.dot(data, rotation_matrix)+ center