Python numpy 模块，c_() 实例源码

def array2PIL(arr, size):
    mode = 'RGBA'
    arr = arr.reshape(arr.shape[0]*arr.shape[1], arr.shape[2])
    if len(arr[0]) == 3:
        arr = numpy.c_[arr, 255*numpy.ones((len(arr),1), numpy.uint8)]
    return Image.frombuffer(mode, size, arr.tostring(), 'raw', mode, 0, 1)

def best_fit_plane(self):
        """Fits a plane to the point cloud using least squares.

        Returns
        -------
        :obj:`tuple` of :obj:`numpy.ndarray` of float
            A normal vector to and point in the fitted plane.
        """
        X = np.c_[self.x_coords, self.y_coords, np.ones(self.num_points)]
        y = self.z_coords
        A = X.T.dot(X)
        b = X.T.dot(y)
        w = np.linalg.inv(A).dot(b)
        n = np.array([w[0], w[1], -1])
        n = n / np.linalg.norm(n)
        n = Direction(n, self._frame)
        x0 = self.mean()
        return n, x0

def predict_bbox_regressor(model, feat, ex_boxes):
  if ex_boxes.size == 0:
    return np.array([]).reshape(-1, 4)

  # predict regression targets
  Y = np.dot(feat, model.Beta[:-1]) + model.Beta[-1]

  # invert transformation
  Y = dot(Y, model.T_inv)

  # read out prediction
  dst_size = Y[:, 2:]
  dst_ctr = Y[:, 2:]

  src_size = ex_boxes[:, 2:]
  src_ctr = ex_boxes[:, :2] + 0.5 * src_size

  pred_size = np.exp(dst_size) * src_size
  pred_ctr = dst_ctr * src_ctr + src_ctr

  pred = np.c_[pred_ctr - 0.5 * pred_size, pred_size]

  return pred

def get_examples(bbox, gt):
  # compute overlap ratio
  O = overlap_ratio(bbox, gt)

  # compute answer
  src_size = bbox[:, 2:]
  src_ctr = bbox[:, :2] + 0.5 * src_size

  gt_size = gt[2:]
  gt_ctr = gt[:2] + 0.5 * gt_size

  dst_size = np.log(gt_size / src_size)
  dst_ctr = (gt_ctr - src_ctr) * 1.0 / src_ctr

  Y = np.c_[dst_ctr, dst_size]

  return Y, O

def computeRotMatrix(self,Phi=False):

        #######################################
        # COMPUTE ROTATION MATRIX SUCH THAT m(t) = A*L(t)*A'*Hp
        # Default set such that phi1,phi2 = 0 is UXO pointed towards North

        if Phi is False:
            phi1 = np.radians(self.phi[0])
            phi2 = np.radians(self.phi[1])
            phi3 = np.radians(self.phi[2])
        else:
            phi1 = np.radians(Phi[0])           # Roll (CCW)
            phi2 = np.radians(Phi[1])           # Inclination (+ve is nose pointing down)
            phi3 = np.radians(Phi[2])           # Declination (degrees CW from North)

        # A1 = np.r_[np.c_[np.cos(phi1),-np.sin(phi1),0.],np.c_[np.sin(phi1),np.cos(phi1),0.],np.c_[0.,0.,1.]] # CCW Rotation about z-axis
        # A2 = np.r_[np.c_[1.,0.,0.],np.c_[0.,np.cos(phi2),np.sin(phi2)],np.c_[0.,-np.sin(phi2),np.cos(phi2)]] # CW Rotation about x-axis (rotates towards North)
        # A3 = np.r_[np.c_[np.cos(phi3),-np.sin(phi3),0.],np.c_[np.sin(phi3),np.cos(phi3),0.],np.c_[0.,0.,1.]] # CCW Rotation about z-axis (direction of head of object)

        A1 = np.r_[np.c_[np.cos(phi1),np.sin(phi1),0.],np.c_[-np.sin(phi1),np.cos(phi1),0.],np.c_[0.,0.,1.]] # CW Rotation about z-axis
        A2 = np.r_[np.c_[1.,0.,0.],np.c_[0.,np.cos(phi2),np.sin(phi2)],np.c_[0.,-np.sin(phi2),np.cos(phi2)]] # CW Rotation about x-axis (rotates towards North)
        A3 = np.r_[np.c_[np.cos(phi3),np.sin(phi3),0.],np.c_[-np.sin(phi3),np.cos(phi3),0.],np.c_[0.,0.,1.]] # CW Rotation about z-axis (direction of head of object)

        return np.dot(A3,np.dot(A2,A1))

def defineSensorLoc(self,XYZ):
        #############################################
        # DEFINE TRANSMITTER AND RECEIVER LOCATIONS
        #
        # XYZ: N X 3 array containing transmitter center locations
        # **NOTE** for this sensor, we know where the receivers are relative to transmitters
        self.TxLoc = XYZ

        dx,dy = np.meshgrid([-0.8,-0.4,0.,0.4,0.8],[-0.8,-0.4,0.,0.4,0.8])
        dx = mkvc(dx)
        dy = mkvc(dy)

        N = np.shape(XYZ)[0]

        X = np.kron(XYZ[:,0],np.ones((25))) + np.kron(np.ones((N)),dx)
        Y = np.kron(XYZ[:,1],np.ones((25))) + np.kron(np.ones((N)),dy)
        Z = np.kron(XYZ[:,2],np.ones((25)))

        self.RxLoc = np.c_[X,Y,Z]

        self.TxID = np.kron(np.arange(1,np.shape(XYZ)[0]+1),np.ones((25)))
        self.RxComp = np.kron(3*np.ones(np.shape(XYZ)[0]),np.ones((25)))

def updatePolarizations(self,r0,UB):

        # Set operator and solution array
        Hp = self.computeHp(r0=r0)
        Brx = self.computeBrx(r0=r0)
        P = self.computeP(Hp,Brx)
        dunc = self.dunc
        dobs = self.dobs

        K = np.shape(dobs)[1]
        q = np.zeros((6,K))

        lb = np.zeros(6)
        ub = UB*np.ones(6)

        for kk in range(0,K):

            LHS = P/np.c_[dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk]]
            RHS = dobs[:,kk]/dunc[:,kk]
            Sol = op.lsq_linear(LHS,RHS,bounds=(lb,ub),tol=1e-5)
            q[:,kk] = Sol.x

        self.q = q

def updatePolarizations(self,r0,UB):

        # Set operator and solution array
        Hp = self.computeHp(r0=r0)
        Brx = self.computeBrx(r0=r0)
        P = self.computeP(Hp,Brx)
        dunc = self.dunc
        dobs = self.dobs

        K = np.shape(dobs)[1]
        q = np.zeros((6,K))

        lb = np.zeros(6)
        ub = UB*np.ones(6)

        for kk in range(0,K):

            LHS = P/np.c_[dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk]]
            RHS = dobs[:,kk]/dunc[:,kk]
            Sol = op.lsq_linear(LHS,RHS,bounds=(lb,ub),tol=1e-7)
            q[:,kk] = Sol.x

        self.q = q

def analytic_infinite_wire(obsloc,wireloc,orientation,I=1.):
    """
    Compute the response of an infinite wire with orientation 'orientation'
    and current I at the obsvervation locations obsloc

    Output:
    B: magnetic field [Bx,By,Bz]
    """

    n,d = obsloc.shape
    t,d = wireloc.shape
    d = np.sqrt(np.dot(obsloc**2.,np.ones([d,t]))+np.dot(np.ones([n,d]),(wireloc.T)**2.)
    - 2.*np.dot(obsloc,wireloc.T))
    distr = np.amin(d, axis=1, keepdims = True)
    idxmind = d.argmin(axis=1)
    r = obsloc - wireloc[idxmind]

    orient = np.c_[[orientation for i in range(obsloc.shape[0])]]
    B = (mu_0*I)/(2*np.pi*(distr**2.))*np.cross(orientation,r)

    return B

def predict_with_glm(X, y, model):
    """ Predict number of mutation with GLM.

    Args:
        X (np.array): feature matrix.
        y (pd.df): response.
        model (dict): model meta-data.

    Returns:
        np.array: array of predictions.

    """
    # Add const. to X
    X = np.c_[X, np.ones(X.shape[0])]
    if model['model_name'] == 'Binomial':
        pred = np.array(model['model'].predict(X) * y.length * y.N)
    elif model['model_name'] == 'NegativeBinomial':
        pred = np.array(model['model'].predict(X, exposure=(y.length * y.N).values + 1))
    else:
        sys.stderr.write('Wrong model name in model info: {}. Need Binomial or NegativeBinomial.'.format(model['model_name']))
        sys.exit(1)
    return pred

def normals(self):
        """Face Normals

        :rtype: numpy.array
        :return: normals, (sum(nF), dim)
        """
        if self.dim == 2:
            nX = np.c_[
                np.ones(self.nFx), np.zeros(self.nFx)
            ]
            nY = np.c_[
                np.zeros(self.nFy), np.ones(self.nFy)
            ]
            return np.r_[nX, nY]
        elif self.dim == 3:
            nX = np.c_[
                np.ones(self.nFx), np.zeros(self.nFx), np.zeros(self.nFx)
            ]
            nY = np.c_[
                np.zeros(self.nFy), np.ones(self.nFy), np.zeros(self.nFy)
            ]
            nZ = np.c_[
                np.zeros(self.nFz), np.zeros(self.nFz), np.ones(self.nFz)
            ]
            return np.r_[nX, nY, nZ]

def tangents(self):
        """Edge Tangents

        :rtype: numpy.array
        :return: normals, (sum(nE), dim)
        """
        if self.dim == 2:
            tX = np.c_[
                np.ones(self.nEx), np.zeros(self.nEx)
            ]
            tY = np.c_[
                np.zeros(self.nEy), np.ones(self.nEy)
            ]
            return np.r_[tX, tY]
        elif self.dim == 3:
            tX = np.c_[
                np.ones(self.nEx), np.zeros(self.nEx), np.zeros(self.nEx)
            ]
            tY = np.c_[
                np.zeros(self.nEy), np.ones(self.nEy), np.zeros(self.nEy)
            ]
            tZ = np.c_[
                np.zeros(self.nEz), np.zeros(self.nEz), np.ones(self.nEz)
            ]
            return np.r_[tX, tY, tZ]

def test_invPropertyTensor2D(self):
        M = discretize.TensorMesh([6, 6])
        a1 = np.random.rand(M.nC)
        a2 = np.random.rand(M.nC)
        a3 = np.random.rand(M.nC)
        prop1 = a1
        prop2 = np.c_[a1, a2]
        prop3 = np.c_[a1, a2, a3]

        for prop in [4, prop1, prop2, prop3]:
            b = invPropertyTensor(M, prop)
            A = makePropertyTensor(M, prop)
            B1 = makePropertyTensor(M, b)
            B2 = invPropertyTensor(M, prop, returnMatrix=True)

            Z = B1*A - sp.identity(M.nC*2)
            self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
            Z = B2*A - sp.identity(M.nC*2)
            self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)

def test_TensorType3D(self):
        M = discretize.TensorMesh([6, 6, 7])
        a1 = np.random.rand(M.nC)
        a2 = np.random.rand(M.nC)
        a3 = np.random.rand(M.nC)
        a4 = np.random.rand(M.nC)
        a5 = np.random.rand(M.nC)
        a6 = np.random.rand(M.nC)
        prop1 = a1
        prop2 = np.c_[a1, a2, a3]
        prop3 = np.c_[a1, a2, a3, a4, a5, a6]

        for ii, prop in enumerate([4, prop1, prop2, prop3]):
            self.assertTrue(TensorType(M, prop) == ii)

        self.assertRaises(Exception, TensorType, M, np.c_[a1, a2, a3, a3])
        self.assertTrue(TensorType(M, None) == -1)

def test_invPropertyTensor3D(self):
        M = discretize.TensorMesh([6, 6, 6])
        a1 = np.random.rand(M.nC)
        a2 = np.random.rand(M.nC)
        a3 = np.random.rand(M.nC)
        a4 = np.random.rand(M.nC)
        a5 = np.random.rand(M.nC)
        a6 = np.random.rand(M.nC)
        prop1 = a1
        prop2 = np.c_[a1, a2, a3]
        prop3 = np.c_[a1, a2, a3, a4, a5, a6]

        for prop in [4, prop1, prop2, prop3]:
            b = invPropertyTensor(M, prop)
            A = makePropertyTensor(M, prop)
            B1 = makePropertyTensor(M, b)
            B2 = invPropertyTensor(M, prop, returnMatrix=True)

            Z = B1*A - sp.identity(M.nC*3)
            self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
            Z = B2*A - sp.identity(M.nC*3)
            self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)

def getError(self):

        funR = lambda r, z: np.sin(2.*np.pi*r)
        funZ = lambda r, z: np.sin(2.*np.pi*z)

        sol = lambda r, t, z: (2*np.pi*r*np.cos(2*np.pi*r) + np.sin(2*np.pi*r))/r + 2*np.pi*np.cos(2*np.pi*z)

        Fc = cylF2(self.M, funR, funZ)
        Fc = np.c_[Fc[:, 0], np.zeros(self.M.nF), Fc[:, 1]]
        F = self.M.projectFaceVector(Fc)

        divF = self.M.faceDiv.dot(F)
        divF_ana = call3(sol, self.M.gridCC)

        err = np.linalg.norm((divF-divF_ana), np.inf)
        return err

def getError(self):

        funR = lambda r, z: np.sin(2.*np.pi*z) * np.cos(np.pi*r)
        funZ = lambda r, z: np.sin(3.*np.pi*z) * np.cos(2.*np.pi*r)

        Fc = cylF2(self.M, funR, funZ)
        Fc = np.c_[Fc[:, 0], np.zeros(self.M.nF), Fc[:, 1]]
        F = self.M.projectFaceVector(Fc)

        aveF = self.M.aveF2CCV * F

        aveF_anaR = funR(self.M.gridCC[:, 0], self.M.gridCC[:, 2])
        aveF_anaZ = funZ(self.M.gridCC[:, 0], self.M.gridCC[:, 2])

        aveF_ana = np.hstack([aveF_anaR, aveF_anaZ])

        err = np.linalg.norm((aveF-aveF_ana), np.inf)
        return err

def surface_points(self, grid_basis=True):
        """Returns the points on the surface.

        Parameters
        ----------
        grid_basis : bool
            If False, the surface points are transformed to the world frame.
            If True (default), the surface points are left in grid coordinates.

        Returns
        -------
        :obj:`tuple` of :obj:`numpy.ndarray` of int, :obj:`numpy.ndarray` of float
            The points on the surface and the signed distances at those points.
        """
        surface_points = np.where(np.abs(self.data_) < self.surface_thresh_)
        x = surface_points[0]
        y = surface_points[1]
        z = surface_points[2]
        surface_points = np.c_[x, np.c_[y, z]]
        surface_vals = self.data_[surface_points[:,0], surface_points[:,1], surface_points[:,2]]
        if not grid_basis:
            surface_points = self.transform_pt_grid_to_obj(surface_points.T)
            surface_points = surface_points.T

        return surface_points, surface_vals

def _check_freq(f):
    """Check the frequency definition."""
    f = np.atleast_2d(np.asarray(f))
    #
    if len(f.reshape(-1)) == 1:
        raise ValueError("The length of f should at least be 2.")
    elif 2 in f.shape:  # f of shape (N, 2) or (2, N)
        if f.shape[1] is not 2:
            f = f.T
    elif np.squeeze(f).shape == (4,):  # (fstart, fend, fwidth, fstep)
        f = _pair_vectors(*tuple(np.squeeze(f)))
    else:  # Sequential
        f = f.reshape(-1)
        f.sort()
        f = np.c_[f[0:-1], f[1::]]
    return f

def convert_mask_to_locations(mask):
    """ Return the converted Cartesian mask as sampling locations.

    Parameters
    ----------
    mask: np.ndarray, {0,1}
        2D matrix, not necessarly a square matrix.

    Returns
    -------
    samples_locations: np.ndarray
        list of the samples between [-0.5, 0.5[.
    """
    row, col = np.where(mask == 1)
    row = row.astype("float") / mask.shape[0] - 0.5
    col = col.astype("float") / mask.shape[1] - 0.5
    return np.c_[row, col]

def _test_double_optimization():
    """Test double optimization on a simple example."""
    # A simple sparse-sum function
    X = [[1, 2], [3, 4], [5, 6]]
    y = [sum(x) for x in X]
    T = [[7, 8], [9, 10], [2, 1]]

    # noisy variables
    np.random.seed(0)
    X = np.c_[X, np.random.random((3, 100))]
    T = np.c_[T, np.random.random((3, 100))]    

    # Select the first 2 variables and calculate a linear model on them
    dstep = DoubleStepEstimator(Lasso(tau=1.0), RidgeRegression(mu=0.0)).train(X, y)

    # Coefficients
    lasso = dstep.selector
    ridge = dstep.estimator
    assert_array_almost_equal([0.90635646, 0.90635646], lasso.beta[:2])
    assert_array_almost_equal([1.0, 1.0], ridge.beta)
    assert_array_almost_equal([1.0, 1.0], dstep.beta[:2])

    # Prediction
    y_ = dstep.predict(T)
    assert_array_almost_equal([15., 19., 3.], y_)

def generate_hills(width, height, nhills):
    '''
    @param width float, terrain width
    @param height float, terrain height
    @param nhills int, #hills to gen. #hills actually generted is sqrt(nhills)^2
    '''
    # setup coordinate grid
    xmin, xmax = -width/2.0, width/2.0
    ymin, ymax = -height/2.0, height/2.0
    x, y = np.mgrid[xmin:xmax:STEP, ymin:ymax:STEP]
    pos = np.empty(x.shape + (2,))
    pos[:, :, 0] = x; pos[:, :, 1] = y

    # generate hilltops
    xm, ym = np.mgrid[xmin:xmax:width/np.sqrt(nhills), ymin:ymax:height/np.sqrt(nhills)]
    mu = np.c_[xm.flat, ym.flat]
    sigma = float(width*height)/(nhills*8)
    for i in range(mu.shape[0]):
        mu[i] = multivariate_normal.rvs(mean=mu[i], cov=sigma)

    # generate hills
    sigma = sigma + sigma*np.random.rand(mu.shape[0])
    rvs = [ multivariate_normal(mu[i,:], cov=sigma[i]) for i in range(mu.shape[0]) ]
    hfield = np.max([ rv.pdf(pos) for rv in rvs ], axis=0)
    return x, y, hfield

def generate_data(sample_size=200, pd=[[0.4, 0.4], [0.1, 0.1]]):
    pd = np.array(pd)
    pd /= pd.sum()
    offset = 50
    bins = np.r_[np.zeros((1,)), np.cumsum(pd)]
    bin_counts = np.histogram(np.random.rand(sample_size), bins)[0]
    data = np.empty((0, 2))
    targets = []
    for ((i, j), p), count in zip(np.ndenumerate(pd), bin_counts):
        xs = np.random.uniform(low=0.0, high=50.0, size=count) + j * offset
        ys = np.random.uniform(low=0.0, high=50.0, size=count) + -i * offset
        data = np.vstack((data, np.c_[xs, ys]))
        if i == j:
            targets.extend([1] * count)
        else:
            targets.extend([-1] * count)
    return np.c_[data, targets]

def show_classification_areas(X, Y, lr):
    x_min, x_max = X[:, 0].min() - .5, X[:, 0].max() + .5
    y_min, y_max = X[:, 1].min() - .5, X[:, 1].max() + .5
    xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.02), np.arange(y_min, y_max, 0.02))
    Z = lr.predict(np.c_[xx.ravel(), yy.ravel()])

    Z = Z.reshape(xx.shape)
    plt.figure(1, figsize=(30, 25))
    plt.pcolormesh(xx, yy, Z, cmap=plt.cm.Pastel1)

    # Plot also the training points
    plt.scatter(X[:, 0], X[:, 1], c=np.abs(Y - 1), edgecolors='k', cmap=plt.cm.coolwarm)
    plt.xlabel('X')
    plt.ylabel('Y')

    plt.xlim(xx.min(), xx.max())
    plt.ylim(yy.min(), yy.max())
    plt.xticks(())
    plt.yticks(())

    plt.show()

def plot_decision_boundary(X, Y, model):
    # X - some data in 2dimensional np.array
    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.01),
                         np.arange(y_min, y_max, 0.01))

    # here "model" is your model's prediction (classification) function
    Z = model(np.c_[xx.ravel(), yy.ravel()])

    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    plt.contourf(xx, yy, Z, cmap=plt.cm.Paired)
    plt.axis('off')

    for i in x:
        print i

    # Plot also the training points
    plt.scatter(X[:, 0], X[:, 1], c=Y, cmap=plt.cm.Paired)


#???????

def get_double_integrator(m=1000, b=50, d=1):
    N = 2
    sys = LTISystem(
        np.c_[[0, 1], [1, -b/m]],  # A
        np.r_[0, d/m],  # B
        # np.r_[0, 1],  # C
    )

    def ref_func(*args):
        if len(args) == 1:
            x = np.zeros(N)
        else:
            x = args[1]
        return np.r_[d/m, 0]-x
    ref = SystemFromCallable(ref_func, N, N)

    return sys, ref

def get_electromechanical(b=1, R=1, L=1, K=np.pi/5, M=1):
    # TODO: determine good reference and/or initial_condition
    # TODO: determine good default values for b, R, L, M
    N = 3
    sys = LTISystem(
        np.c_[  # A
            [0, 0, 0],
            [1, -b/M, -K/L],
            [0, K/M, -R/L]
        ],
        np.r_[0, 0, 1/L],  # B
        # np.r_[1, 0, 0],  # C
    )
    sys.initial_condition = np.ones(N)

    def ref_func(*args):
        if len(args) == 1:
            x = np.zeros(N)
        else:
            x = args[1]
        return np.r_[0, 0, 0]-x
    ref = SystemFromCallable(ref_func, N, N)

    return sys, ref

def get_cart_pendulum(m=1, M=3, L=0.5, g=9.81, pedant=False):
    N = 4
    sys = LTISystem(
        np.c_[  # A
            [0, 0, 0, 0],
            [1, 0, 0, 0],
            [0, m*g/M, 0, (-1)**(pedant)*(m+M)*g/(M*L)],
            [0, 0, 1, 0]
        ],
        np.r_[0, 1/M, 0, 1/(M*L)],  # B
        # np.r_[1, 0, 1, 0]  # C
    )
    sys.initial_condition = np.r_[0, 0, np.pi/3, 0]

    def ref_func(*args):
        if len(args) == 1:
            x = np.zeros(N)
        else:
            x = args[1]
        return np.zeros(N)-x
    ref = SystemFromCallable(ref_func, N, N)

    return sys, ref

def update_equation_function(self, *args):
        event_var = self.event_variable_equation_function(*args)
        if self.condition_idx is None:
            self.condition_idx = np.where(np.all(np.r_[
                    np.c_[[[True]], event_var >= self.event_bounds],
                    np.c_[event_var <= self.event_bounds, [[True]]]
                    ], axis=0))[0][0]
        else:
            sq_dist = (event_var - self.event_bounds)**2
            crossed_root_idx = np.where(sq_dist == np.min(sq_dist))[1][0]
            if crossed_root_idx == self.condition_idx:
                self.condition_idx += 1
            elif crossed_root_idx == self.condition_idx-1:
                self.condition_idx -= 1
            else:
                warnings.warn("SwitchedSystem did not cross a neighboring " +
                              "boundary. This may indicate an integration " +
                              "error. Continuing without updating " +
                              "condition_idx", UserWarning)
        return self.state_update_equation_function(*args)

def plane2xyz(center, ij, plane):
        """
        converts image pixel indices to xyz on the PLANE.

        center : 2-tuple
        ij : nx2 int array
        plane : 4-tuple

        return nx3 array.
        """
        ij = np.atleast_2d(ij)
        n = ij.shape[0]
        ij = ij.astype('float')
        xy_ray = (ij-center[None,:]) / DepthCamera.f
        z = -plane[2]/(xy_ray.dot(plane[:2])+plane[3])
        xyz = np.c_[xy_ray, np.ones(n)] * z[:,None]
        return xyz

def viz_textbb(fignum,text_im, bb_list,alpha=1.0):
    """
    text_im : image containing text
    bb_list : list of 2x4xn_i boundinb-box matrices
    """
    plt.close(fignum)
    plt.figure(fignum)
    plt.imshow(text_im)
    plt.hold(True)
    H,W = text_im.shape[:2]
    for i in xrange(len(bb_list)):
        bbs = bb_list[i]
        ni = bbs.shape[-1]
        for j in xrange(ni):
            bb = bbs[:,:,j]
            bb = np.c_[bb,bb[:,0]]
            plt.plot(bb[0,:], bb[1,:], 'r', linewidth=2, alpha=alpha)
    plt.gca().set_xlim([0,W-1])
    plt.gca().set_ylim([H-1,0])
    plt.show(block=False)

def bbox(self, image_coords=True):
        """
        Return a 3 by 2 matrix, corresponding to the bounding box of the 
        annotation within the scan. If `scan_slice` is a numpy array 
        containing aslice of the scan, each slice of the annotation is 
        contained within the box:

            bbox[0,0]:bbox[0,1]+1, bbox[1,0]:bbox[1,1]+1

        If `image_coords` is `False` then each annotation slice is 
        instead contained within:

            bbox[1,0]:bbox[1,1]+1, bbox[0,0]:bbox[0,1]+1

        The last row of `bbox` give the inclusive lower and upper 
        bounds of the `image_z_position`.
        """
        matrix = self.contours_to_matrix()
        bbox   = np.c_[matrix.min(axis=0), matrix.max(axis=0)]
        return bbox if not image_coords else bbox[[1,0,2]]

def plot_and_save(x, func, xvv_inst, plot=False, fname=False):
    mat = x.T
    for m in range(xvv_inst.nsites):
        for n in range(m+1):
            if fname:
                mat = np.c_[mat, func[:,m,n].T]
            if plot:
                plt.plot(x, func[:,m,n], 
                         label='{}-{}'.format(xvv_inst.atom_names[m], 
                                              xvv_inst.atom_names[n]))
    if fname:
        np.savetxt(fname, mat)
    if plot:
        plt.legend()
        plt.savefig('graph.png', dpi=300)
        plt.show()

def read_csv(filename, skip_lines=0):
    csvfile = file(filename, 'rb')
    reader = csv.reader(csvfile)
    data = np.empty(0, dtype=object)
    last_count = np.NAN
    for line in reader:
        if skip_lines > 0:
            skip_lines = skip_lines - 1
            continue
        if data.size > 0:
            if len(line) != last_count:
                raise Exception('unequal columes found')
            data = np.c_[data, line]
            last_count = len(line)
        else:
            data = np.array(line, dtype=object)
            data = data.reshape(len(data), 1)
            last_count = len(line)
    csvfile.close()
    return data.T

def print_evaluation_result(clf, bags_test, args):
  pred_score = np.array([clf(B.data()) for B in bags_test])
  pred_label = np.array([1 if score >= 0 else -1 for score in pred_score])
  true_label = np.array([B.y for B in bags_test])
  a = accuracy (pred_label, true_label)  # accuracy
  p = precision(pred_label, true_label)  # precision
  r = recall   (pred_label, true_label)  # recall
  f = f_score  (pred_label, true_label)  # F-score
  auc = metrics.roc_auc_score((true_label+1)/2, pred_score)

  if not args.aucplot:
    sys.stdout.write("""# accuracy,precision,recall,f-score,ROC-AUC
{:.3f},{:.3f},{:.3f},{:.3f},{:.3f}\n""".format(a, p, r, f, auc))
    sys.stdout.flush()

  else:
    sys.stdout.write("""# accuracy,precision,recall,f-score,ROC-AUC
# {:.3f},{:.3f},{:.3f},{:.3f},{:.3f}\n""".format(a, p, r, f, auc))
    sys.stdout.flush()
    np.savetxt(sys.stdout.buffer, np.c_[pred_score, true_label])

def trotx(theta, unit="rad", xyz=[0, 0, 0]):
    """
    TROTX Rotation about X axis

    :param theta: rotation in radians or degrees
    :param unit: "rad" or "deg" to indicate unit being used
    :param xyz: the xyz translation, if blank defaults to [0,0,0]
    :return: homogeneous transform matrix

    trotx(THETA) is a homogeneous transformation (4x4) representing a rotation
    of THETA radians about the x-axis.
    trotx(THETA, 'deg') as above but THETA is in degrees
    trotx(THETA, 'rad', [x,y,z]) as above with translation of [x,y,z]
    """
    check_args.unit_check(unit)
    tm = rotx(theta, unit)
    tm = np.r_[tm, np.zeros((1, 3))]
    mat = np.c_[tm, np.array([[xyz[0]], [xyz[1]], [xyz[2]], [1]])]
    mat = np.asmatrix(mat.round(15))
    return mat


# ---------------------------------------------------------------------------------------#

def troty(theta, unit="rad", xyz=[0, 0, 0]):
    """
    TROTY Rotation about Y axis

    :param theta: rotation in radians or degrees
    :param unit: "rad" or "deg" to indicate unit being used
    :param xyz: the xyz translation, if blank defaults to [0,0,0]
    :return: homogeneous transform matrix

    troty(THETA) is a homogeneous transformation (4x4) representing a rotation
    of THETA radians about the y-axis.
    troty(THETA, 'deg') as above but THETA is in degrees
    troty(THETA, 'rad', [x,y,z]) as above with translation of [x,y,z]
    """
    check_args.unit_check(unit)
    tm = roty(theta, unit)
    tm = np.r_[tm, np.zeros((1, 3))]
    mat = np.c_[tm, np.array([[xyz[0]], [xyz[1]], [xyz[2]], [1]])]
    mat = np.asmatrix(mat.round(15))
    return mat


# ---------------------------------------------------------------------------------------#

def trotz(theta, unit="rad", xyz=[0, 0, 0]):
    """
    TROTZ Rotation about Z axis

    :param theta: rotation in radians or degrees
    :param unit: "rad" or "deg" to indicate unit being used
    :param xyz: the xyz translation, if blank defaults to [0,0,0]
    :return: homogeneous transform matrix

    trotz(THETA) is a homogeneous transformation (4x4) representing a rotation
    of THETA radians about the z-axis.
    trotz(THETA, 'deg') as above but THETA is in degrees
    trotz(THETA, 'rad', [x,y,z]) as above with translation of [x,y,z]
    """
    check_args.unit_check(unit)
    tm = rotz(theta, unit)
    tm = np.r_[tm, np.zeros((1, 3))]
    mat = np.c_[tm, np.array([[xyz[0]], [xyz[1]], [xyz[2]], [1]])]
    mat = np.asmatrix(mat.round(15))
    return mat


# ---------------------------------------------------------------------------------------#

def trot2(theta, unit='rad'):
    """
    TROT2 SE2 rotation matrix

    :param theta: rotation in radians or degrees
    :param unit: "rad" or "deg" to indicate unit being used
    :return: homogeneous transform matrix (3x3)

    TROT2(THETA) is a homogeneous transformation (3x3) representing a rotation of
    THETA radians.
    TROT2(THETA, 'deg') as above but THETA is in degrees.
    Notes::
    - Translational component is zero.
    """
    tm = rot2(theta, unit)
    tm = np.r_[tm, np.zeros((1, 2))]
    mat = np.c_[tm, np.array([[0], [0], [1]])]
    return mat


# ---------------------------------------------------------------------------------------#

def _scipy_bivariate_kde(x, y, bw, gridsize, cut, clip):
    """Compute a bivariate kde using scipy."""
    data = np.c_[x, y]
    kde = stats.gaussian_kde(data.T)
    data_std = data.std(axis=0, ddof=1)
    if isinstance(bw, string_types):
        bw = "scotts" if bw == "scott" else bw
        bw_x = getattr(kde, "%s_factor" % bw)() * data_std[0]
        bw_y = getattr(kde, "%s_factor" % bw)() * data_std[1]
    elif np.isscalar(bw):
        bw_x, bw_y = bw, bw
    else:
        msg = ("Cannot specify a different bandwidth for each dimension "
               "with the scipy backend. You should install statsmodels.")
        raise ValueError(msg)
    x_support = _kde_support(data[:, 0], bw_x, gridsize, cut, clip[0])
    y_support = _kde_support(data[:, 1], bw_y, gridsize, cut, clip[1])
    xx, yy = np.meshgrid(x_support, y_support)
    z = kde([xx.ravel(), yy.ravel()]).reshape(xx.shape)
    return xx, yy, z

def spiral(num_cls, dim, point_per_cls, rnd_state=1024):
    np.random.seed(rnd_state)
    points_per_cls = 100  # number of points per class
    dim = 2  # dimensionality
    num_cls = 3  # number of classes
    X_data = np.zeros((points_per_cls * num_cls, dim))
    y_data = np.zeros(points_per_cls * num_cls, dtype='uint8')
    for j in range(num_cls):
        ix = range(points_per_cls * j, points_per_cls * (j + 1))
        r = np.linspace(0.0, 1, points_per_cls)
        t = np.linspace(j * 4, (j + 1) * 4, points_per_cls) + np.random.randn(points_per_cls) * 0.2  # theta
        X_data[ix] = np.c_[r * np.sin(t), r * np.cos(t)]
        y_data[ix] = j

    y_data_encoded = np.zeros((points_per_cls * num_cls, num_cls))
    y_data_encoded[range(points_per_cls * num_cls), y_data] = 1
    return X_data, y_data, y_data_encoded

def gradient(x0, X, y, alpha):
    # gradient of the logistic loss

    w, c = x0[1:137], x0[0]

    #print("c is " + str(c))
    z = X.dot(w) + c
    z = phi(y * z)
    z0 = (z - 1) * y
    grad_w = np.matmul(z0,X) / X.shape[0] + alpha * w
    grad_c = z0.sum() / X.shape[0]

    grad_c = np.array(grad_c)
    #print(grad_w[0,1:5])
    return np.c_[([grad_c], grad_w)]


##### Stochastic Gradient Descent Optimiser ######

def average_ndcg(labels, query_ids, predicted_labels):
    ndcg_list = np.zeros(len(set(query_ids)))
    k = 0
    for i in set(query_ids):
        idx = [query_ids == i]
        orders = np.c_[labels[idx],predicted_labels[idx]]

        sorted_orders = orders[orders[:,1].argsort()[::-1]][:,0]
        ndcg_list[k] = ndcg(sorted_orders)

        k +=1
        if k%2000 == 0:
            print(str(k) + " queries calculated")
            print("mean ndcg so far: " + str(np.mean(ndcg_list[0:k])))
    return np.mean(ndcg_list)


# average ndcg is 0.26333

def average_ndcg(labels, query_ids, predicted_labels):
    ndcg_list = np.zeros(len(set(query_ids)))
    k = 0
    for i in set(query_ids):
        idx = [query_ids == i]
        orders = np.c_[labels[idx],predicted_labels[idx]]

        sorted_orders = orders[orders[:,1].argsort()[::-1]][:,0]
        ndcg_list[k] = ndcg(sorted_orders)

        k +=1
        if k%2000 == 0:
            print(str(k) + " queries calculated")
            print("mean ndcg so far: " + str(np.mean(ndcg_list[0:k])))
    return np.mean(ndcg_list)


# average ndcg is 0.26333

def computeGaussianWidthCandidates(self, referenceSamples=None, testSamples=None) :
        """
        Compute a candidate list of Gaussian kernel widths. The best width will be
        selected via cross-validation
        """
        allSamples     = numpy.c_[referenceSamples, testSamples]
        medianDistance = self.getMedianDistanceBetweenSamples(allSamples.T)

        return medianDistance * numpy.array([0.6, 0.8, 1, 1.2, 1.4])

def _fit(self,X,y=None):
        if isinstance(X, pyisc.DataObject) and y is None:
            assert y is None # Contained in the data object
            self.class_column = X.class_column
            if self.class_column >= 0:
                self.classes_ = X.classes_

            self._anomaly_detector._SetParams(
                0,
                -1 if X.class_column is None else X.class_column,
                self.anomaly_threshold,
                1 if self.is_clustering else 0
            )
            self._anomaly_detector._TrainData(X)
            return self
        if isinstance(X, ndarray):
            class_column = -1
            data_object = None
            assert X.ndim <= 2
            if X.ndim == 2:
                max_class_column = X.shape[1]
            else:
                max_class_column = 1
            if isinstance(y, list) or isinstance(y, ndarray):
                assert len(X) == len(y)
                class_column = max_class_column
                data_object = pyisc.DataObject(numpy.c_[X, y], class_column=class_column)
            elif y is None or int(y) == y and y > -1 and y <= max_class_column:
                self.class_column = y
                data_object = pyisc.DataObject(X,class_column=y)

            if data_object is not None:
                return self._fit(data_object)

        raise ValueError("Unknown type of data to fit X, y:", type(X), type(y))

def _convert_to_data_object_in_scoring(self, X, y):
        data_object = None
        if isinstance(y, list) or isinstance(y, ndarray):
            assert X.ndim == 2 and self.class_column == X.shape[1] or X.ndim == 1 and self.class_column == 1
            data_object = pyisc.DataObject(numpy.c_[X, y], class_column=self.class_column,classes=self.classes_)
        else:
            assert self.class_column == y
            data_object = pyisc.DataObject(X, class_column=self.class_column,classes=self.classes_ if y is not None else None)
        return data_object

def loglikelihood(self,X,y=None):
        assert isinstance(X, ndarray) and (self.class_column is None and y is None or len(y) == len(X))

        if y is not None:
            return self._anomaly_detector._LogProbabilityOfData(pyisc.DataObject(c_[X,y], class_column=len(X[0])), len(X)).sum()
        else:
            return self._anomaly_detector._LogProbabilityOfData(pyisc.DataObject(X), len(X)).sum()

def test_dataobject_set_column_values(self):
        X = array([norm(1.0).rvs(10) for _ in range(1000)])
        y = [None] * 1000

        DO = DataObject(c_[X,y], class_column=len(X[0]))
        assert_equal(len(X[0]), DO.class_column)
        assert_equal(unique(y), DO.classes_)

        classes=[None] + ['1', '2', '3', '4', '5']
        DO = DataObject(c_[X,y], class_column=len(X[0]), classes=classes)
        assert_equal(len(X[0]), DO.class_column)
        assert_equal(classes, DO.classes_)

        X2 = DO.as_2d_array()
        assert_allclose(X2.T[:-1].T.astype(float), X)
        assert_equal(X2.T[-1],y)

        new_y = ["%i"%(divmod(i,5)[1]+1) for i in range(len(X))]
        DO.set_column_values(len(X[0]), new_y)

        assert_equal(len(X[0]), DO.class_column)
        assert_equal([None]+list(unique(new_y)), DO.classes_)

        X2 = DO.as_2d_array()
        assert_allclose(X2.T[:-1].T.astype(float), X)
        assert_equal(X2.T[-1], new_y)