numpy.mat

Here are the examples of the python api numpy.mat taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.

69 Examples 7

Page 1 Selected Page 2

Example 1

Project: python-control Source File: matlab_test.py
    def testSS2cont(self):
        sys = ss(
            np.mat("-3 4 2; -1 -3 0; 2 5 3"),
            np.mat("1 4 ; -3 -3; -2 1"),
            np.mat("4 2 -3; 1 4 3"),
            np.mat("-2 4; 0 1"))
        sysd = c2d(sys, 0.1)
        np.testing.assert_array_almost_equal(
            np.mat(
                """0.742840837331905  0.342242024293711  0.203124211149560;
                  -0.074130792143890  0.724553295044645 -0.009143771143630;
                   0.180264783290485  0.544385612448419  1.370501013067845"""),
            sysd.A)
        np.testing.assert_array_almost_equal(
            np.mat(""" 0.012362066084719   0.301932197918268;
                      -0.260952977031384  -0.274201791021713;
                      -0.304617775734327   0.075182622718853"""), sysd.B)

Example 2

Project: nimfa Source File: linalg.py
Function: trace
def trace(X):
    """
    Return trace of sparse or dense square matrix X.

    :param X: Target matrix.
    :type X: :class:`scipy.sparse` of format csr, csc, coo, bsr, dok, lil,
    dia or :class:`numpy.matrix`
    """
    assert X.shape[0] == X.shape[1], "X should be square matrix."
    if sp.isspmatrix(X):
        return sum(X[i, i] for i in range(X.shape[0]))
    else:
        return np.trace(np.mat(X))

Example 3

Project: RLScore Source File: sq_mprank_measure.py
def sqmprank_multitask(Y, Y_predicted):  
    Y = np.mat(Y)
    Y_predicted = np.mat(Y_predicted)
    vlen = Y.shape[0]
    centeredsqerror = Y - Y_predicted
    onevec = np.mat(np.ones((vlen, 1)))
    tempvec = onevec.T * centeredsqerror
    np.multiply(vlen, centeredsqerror, centeredsqerror)
    np.subtract(centeredsqerror, tempvec, centeredsqerror)
    np.multiply(centeredsqerror, centeredsqerror, centeredsqerror)
    performances = np.mean(centeredsqerror, axis = 0) / ((vlen ** 2 - vlen) / 2)
    performances = np.array(performances)[0]
    return performances

Example 4

Project: RLScore Source File: array_tools.py
def as_dense_matrix(A):
    """Returns the input as matrix

    Parameters
    ----------
    A: {array-like, sparse matrix}, shape = 2D
    
    Returns
    -------
    A : np.matrix
    """
    if sp.issparse(A):
        return A.todense()
    else:
        return np.mat(A)

Example 5

Project: RLScore Source File: query_rankrls.py
    def __init__(self, X, Y, qids, regparam = 1.0, kernel='LinearKernel', basis_vectors = None, **kwargs):
        kwargs["bias"] = 0.
        kwargs['kernel'] =  kernel
        kwargs['X'] = X
        if basis_vectors is not None:
            kwargs['basis_vectors'] = basis_vectors
        self.svdad = adapter.createSVDAdapter(**kwargs)
        self.Y = np.mat(array_tools.as_2d_array(Y))
        self.regparam = regparam
        self.svals = np.mat(self.svdad.svals)
        self.svecs = self.svdad.rsvecs
        self.size = self.Y.shape[0]
        self.size = self.Y.shape[0]
        self.qids = map_qids(qids)
        self.qidlist = qids_to_splits(self.qids)
        self.solve(self.regparam)

Example 6

Project: RLScore Source File: global_rankrls.py
Function: init
    def __init__(self, X, Y, regparam = 1.0, kernel='LinearKernel', basis_vectors = None, **kwargs):
        Y = array_tools.as_2d_array(Y)
        self.Y = np.mat(Y)
        if X.shape[0] != Y.shape[0]:
            raise Exception("First dimension of X and Y must be the same")
        if basis_vectors is not None:
            if X.shape[1] != basis_vectors.shape[1]:
                raise Exception("Number of columns for X and basis_vectors must be the same")
        kwargs["bias"] = 0.
        kwargs['kernel'] =  kernel
        kwargs['X'] = X
        if basis_vectors is not None:
            kwargs['basis_vectors'] = basis_vectors
        self.svdad = adapter.createSVDAdapter(**kwargs)
        self.regparam = regparam
        self.svals = np.mat(self.svdad.svals)
        self.svecs = np.mat(self.svdad.rsvecs)
        self.size = self.Y.shape[0]
        self.solve(self.regparam)

Example 7

Project: nmrglue Source File: proc_lp.py
def find_lpc_svd(D, d):
    """
    Find linear prediction filter using single value decomposition.
    """
    L = D.shape[0]
    m = D.shape[1]
    U, s, Vh = scipy.linalg.svd(D)  # SVD decomposition
    U, Vh = np.mat(U), np.mat(Vh)   # make U and Vh matrices
    Si = pinv_diagsvd(s, m, L)      # construct the pseudo-inverse sigma matrix
    return np.array(Vh.H * Si * U.H * d)

Example 8

Project: RLScore Source File: accuracy_measure.py
def accuracy_multitask(Y, P):
    Y = np.mat(Y)
    P = np.mat(P)
    if not np.all((Y==1) + (Y==-1)):
        raise UndefinedPerformance("binary classification accuracy accepts as Y-values only 1 and -1")
    vlen = float(Y.shape[0])
    performances = np.sum(np.sign(np.multiply(Y, P)) + 1., axis = 0) / (2 * vlen)
    performances = np.array(performances)[0]
    return performances

Example 9

Project: RHEAS Source File: kalman.py
Function: init
    def __init__(self, A, HA, d, E):
        """Initialize Ensemble Kalman Filter object using a state matrix *A*,
        a predicted measurement matrix *HA*, an observation vector *d*,
        and an observation error covariance matrix *R*."""
        if HA is not None:
            self.HA = HA
        else:
            self.HA = A
        self.A = np.mat(A)
        self.HA = np.mat(self.HA)
        self.d = np.mat(d)
        self.E = np.mat(E)
        self.ndim, self.nens = A.shape
        self.nobs = len(d)
        self.R = np.mat(
            np.diag(np.diag((1.0 / (self.nens - 1)) * self.E * self.E.T)))

Example 10

Project: RLScore Source File: sq_mprank_measure.py
def sqmprank_singletask(Y, P):
    correct = Y
    predictions = P
    correct = np.mat(correct)
    predictions = np.mat(predictions)
    vlen = correct.shape[0]
    diff = correct - predictions
    onevec = np.mat(np.ones((vlen, 1)))
    centereddiff = vlen * diff - onevec * (onevec.T * diff)
    sqerror = (centereddiff.T * diff)[0, 0] / ((len(correct) ** 2 - len(correct)) / 2)
    return sqerror

Example 11

Project: RLScore Source File: global_rankrls.py
Function: leave_one_out
    def leave_one_out(self):
        """Computes leave-one-out predictions for a trained RankRLS
        
        Returns
        -------
        F : array, shape = [n_samples, n_labels]
            leave-one-out predictions

        Notes
        -----
    
        Provided for reference, usually you should not call this, but
        rather use leave_pair_out.

        """        
        LOO = mat(zeros((self.size, self.ysize), dtype=float64))
        for i in range(self.size):
            LOO[i,:] = self.holdout([i])
        return LOO

Example 12

Project: RLScore Source File: rls.py
Function: init
    def __init__(self, X, Y, regparam = 1.0, kernel='LinearKernel', basis_vectors = None, **kwargs):
        self.Y = array_tools.as_2d_array(Y)
        if X.shape[0] != Y.shape[0]:
            raise Exception("First dimension of X and Y must be the same")
        if basis_vectors is not None:
            if X.shape[1] != basis_vectors.shape[1]:
                raise Exception("Number of columns for X and basis_vectors must be the same")
        kwargs['X'] = X
        kwargs['kernel'] = kernel
        if basis_vectors is not None:
            kwargs['basis_vectors'] = basis_vectors
        self.svdad = adapter.createSVDAdapter(**kwargs)
        self.regparam = regparam
        self.svals = np.mat(self.svdad.svals)
        self.svecs = np.mat(self.svdad.rsvecs)
        self.size = self.Y.shape[0]
        self.solve(self.regparam)   

Example 13

Project: nimfa Source File: linalg.py
Function: svd
def svd(X):
    """
    Compute standard SVD on matrix X.

    :param X: The input matrix.
    :type X: :class:`scipy.sparse` of format csr, csc, coo, bsr, dok, lil,
    dia or :class:`numpy.matrix`
    """
    if sp.isspmatrix(X):
        if X.shape[0] <= X.shape[1]:
            U, S, V = _svd_left(X)
        else:
            U, S, V = _svd_right(X)
    else:
        U, S, V = nla.svd(np.mat(X), full_matrices=False)
        S = np.mat(np.diag(S))
    return U, S, V

Example 14

Project: python-control Source File: matlab_test.py
Function: test_connect
    def testConnect(self):
        sys1 = ss("1. -2; 3. -4", "5.; 7", "6, 8", "9.")
        sys2 = ss("-1.", "1.", "1.", "0.")
        sys = append(sys1, sys2)
        Q= np.mat([ [ 1, 2], [2, -1] ]) # basically feedback, output 2 in 1
        sysc = connect(sys, Q, [2], [1, 2])
        # print(sysc)
        np.testing.assert_array_almost_equal(
            sysc.A, np.mat('1 -2 5; 3 -4 7; -6 -8 -10'))
        np.testing.assert_array_almost_equal(
            sysc.B, np.mat('0; 0; 1'))
        np.testing.assert_array_almost_equal(
            sysc.C, np.mat('6 8 9; 0 0 1'))
        np.testing.assert_array_almost_equal(
            sysc.D, np.mat('0; 0'))

Example 15

Project: scipy Source File: test_ltisys.py
    def test_jordan_block(self):
        # Non-diagonalizable A matrix
        #   x1' + x1 = x2
        #   x2' + x2 = u
        #   y = x1
        # Exact solution with u = 0 is y(t) = t exp(-t)
        A = np.mat("-1. 1.; 0. -1.")
        B = np.mat("0.; 1.")
        C = np.mat("1. 0.")
        system = self.lti_nowarn(A, B, C, 0.)
        t = np.linspace(0,5)
        u = np.zeros_like(t)
        tout, y, x = lsim(system, u, t, X0=[0.0, 1.0])
        expected_y = tout * np.exp(-tout)
        assert_almost_equal(y, expected_y)

Example 16

Project: scipy Source File: test_ltisys.py
    def test_double_integrator(self):
        # double integrator: y'' = 2u
        A = np.mat("0. 1.; 0. 0.")
        B = np.mat("0.; 1.")
        C = np.mat("2. 0.")
        system = self.lti_nowarn(A, B, C, 0.)
        t = np.linspace(0,5)
        u = np.ones_like(t)
        tout, y, x = lsim(system, u, t)
        expected_x = np.transpose(np.array([0.5 * tout**2, tout]))
        expected_y = tout**2
        assert_almost_equal(x, expected_x)
        assert_almost_equal(y, expected_y)

Example 17

Project: RLScore Source File: sqerror_measure.py
def sqerror_singletask(correct, predictions):
    correct = np.mat(correct)
    predictions = np.mat(predictions)
    diff = correct - predictions
    sqerror = np.dot(diff.T, diff)[0,0]
    sqerror /= correct.shape[0]
    return sqerror

Example 18

Project: nmrglue Source File: proc_lp.py
def find_lpc_cholesky(D, d):
    """
    Find linear prediction filter using a Cholesky decomposition.
    """
    # form the normal equation (D.H*D)*a = D.H*d
    # SPEED
    # this can be improved by using the Hankel nature of D
    D = np.mat(D)
    DhD = np.mat(np.dot(D.H, D))
    Dhd = np.mat(np.dot(D.H, d))

    c, lower = scipy.linalg.cho_factor(DhD)     # Compute Cholesky decomp.
    return scipy.linalg.cho_solve((c, lower), Dhd)  # solve normal equation

Example 19

Project: RLScore Source File: array_tools.py
Function: as_matrix
def as_matrix(A):
    """Returns the input as matrix or sparse matrix

    Parameters
    ----------
    A: {array-like, sparse matrix}, shape = 2D
    
    Returns
    -------
    A : {matrix, sparse matrix}
    """
    if sp.issparse(A):
        return A
    else:
        return np.mat(A)

Example 20

Project: nmrglue Source File: proc_lp.py
def find_lpc_qr(D, d):
    """
    Find linear prediction filter using QR decomposition.
    """
    L = D.shape[0]
    m = D.shape[1]
    q, r = scipy.linalg.qr(D)
    q, r = np.mat(q), np.mat(r)

    # SPEED
    # the next line is slow and the use of pinv2 should be avoided as
    # pseudo inversion of r involves a computationally expensive SVD
    # decomposition which is not needed.  Rather r*x = q.H*d should be
    # solved for x using LAPACK's ZTRTRS function (or similar function with
    # different prefix).  This is not currently available in scipy/numpy and
    # therefore is not used here.
    return scipy.linalg.pinv2(r) * q.H * d

Example 21

Project: RLScore Source File: sqerror_measure.py
def sqerror_multitask(Y, Y_predicted):  
    Y = np.mat(Y)
    Y_predicted = np.mat(Y_predicted)
    sqerror = Y - Y_predicted
    multiply(sqerror, sqerror, sqerror)
    performances = mean(sqerror, axis = 0)
    performances = np.array(performances)[0]
    return performances

Example 22

Project: universal-portfolios Source File: corn.py
    def optimal_weights(self, X):
        X = np.mat(X)

        n,m = X.shape
        P = 2 * matrix(X.T * X)
        q = -3 * matrix(np.ones((1,n)) * X).T

        G = matrix(-np.eye(m))
        h = matrix(np.zeros(m))
        A = matrix(np.ones(m)).T
        b = matrix(1.)

        sol = solvers.qp(P, q, G, h, A, b)
        return np.squeeze(sol['x'])

Example 23

Project: msmbuilder-legacy Source File: lprmsd.py
Function: computeoverlap
def ComputeOverlap(theta, elem, xyz1, xyz2):
    """
    Computes an 'overlap' between two molecules based on some
    fictitious density.  Good for fine-tuning alignment but gets stuck
    in local minima.
    """

    xyz2R = np.array(EulerMatrix(theta[0], theta[1], theta[2]) * np.mat(xyz2.T)).T
    Obj = 0.0
    for i in set(elem):
        for j in np.where(elem == i)[0]:
            for k in np.where(elem == i)[0]:
                dx = xyz1[j] - xyz2R[k]
                dx2 = np.dot(dx, dx)
                Obj -= np.exp(-0.5 * dx2)
    return Obj

Example 24

Project: python-control Source File: matlab_test.py
Function: test_damp
    def testDamp(self):
        A = np.mat('''-0.2  0.06 0    -1;
               0    0    1     0;
             -17    0   -3.8   1;
               9.4  0   -0.4  -0.6''')
        B = np.mat('''-0.01  0.06;
               0     0;
             -32     5.4;
               2.6  -7''')
        C = np.eye(4)
        D = np.zeros((4,2))
        sys = ss(A, B, C, D)
        wn, Z, p = damp(sys, False)
        # print (wn)
        np.testing.assert_array_almost_equal(
            wn, np.array([4.07381994,   3.28874827,   3.28874827,
                          1.08937685e-03]))
        np.testing.assert_array_almost_equal(
            Z, np.array([1.0, 0.07983139,  0.07983139, 1.0]))

Example 25

Project: python-control Source File: matlab_test.py
Function: test_connect2
    def testConnect2(self):
        sys = append(ss([[-5, -2.25], [4, 0]], [[2], [0]],
                          [[0, 1.125]], [[0]]),
                       ss([[-1.6667, 0], [1, 0]], [[2], [0]],
                          [[0, 3.3333]], [[0]]),
                       1)
        Q = [ [ 1, 3], [2, 1], [3, -2]]
        sysc = connect(sys, Q, [3], [3, 1, 2])
        np.testing.assert_array_almost_equal(
            sysc.A, np.mat([[-5, -2.25, 0, -6.6666],
                            [4, 0, 0, 0],
                            [0, 2.25, -1.6667, 0],
                            [0, 0, 1, 0]]))
        np.testing.assert_array_almost_equal(
            sysc.B, np.mat([[2], [0], [0], [0]]))
        np.testing.assert_array_almost_equal(
            sysc.C, np.mat([[0, 0, 0, -3.3333],
                            [0, 1.125, 0, 0],
                            [0, 0, 0, 3.3333]]))
        np.testing.assert_array_almost_equal(
            sysc.D, np.mat([[1], [0], [0]]))

Example 26

Project: RLScore Source File: two_step_rls.py
Function: init
    def __init__(self, **kwargs):
        Y = kwargs["Y"]
        Y = np.mat(array_tools.as_2d_array(Y))
        if kwargs.has_key('K1'):
            K1 = np.mat(kwargs['K1'])
            K2 = np.mat(kwargs['K2'])
            Y = Y.reshape((K1.shape[0], K2.shape[0]), order = 'F')
            self.K1, self.K2 = K1, K2
            self.kernelmode = True
        else:
            X1 = np.mat(kwargs['X1'])
            X2 = np.mat(kwargs['X2'])
            Y = Y.reshape((X1.shape[0], X2.shape[0]), order = 'F')
            self.X1, self.X2 = X1, X2
            self.kernelmode = False
        self.Y = Y
        self.regparam1 = kwargs["regparam1"]
        self.regparam2 = kwargs["regparam2"]
        self.trained = False
        self.solve(self.regparam1, self.regparam2)

Example 27

Project: python-control Source File: timeresp_test.py
    def test_lsim_double_integrator(self):
        # Note: scipy.signal.lsim fails if A is not invertible
        A = np.mat("0. 1.;0. 0.")
        B = np.mat("0.; 1.")
        C = np.mat("1. 0.")
        D = 0.
        sys = StateSpace(A, B, C, D)

        def check(u, x0, xtrue):
            _t, yout, xout = forced_response(sys, t, u, x0)
            np.testing.assert_array_almost_equal(xout, xtrue, decimal=6)
            ytrue = np.squeeze(np.asarray(C.dot(xtrue)))
            np.testing.assert_array_almost_equal(yout, ytrue, decimal=6)

        # test with zero input
        npts = 10
        t = np.linspace(0, 1, npts)
        u = np.zeros_like(t)
        x0 = np.array([2., 3.])
        xtrue = np.zeros((2, npts))
        xtrue[0, :] = x0[0] + t * x0[1]
        xtrue[1, :] = x0[1]
        check(u, x0, xtrue)

        # test with step input
        u = np.ones_like(t)
        xtrue = np.array([0.5 * t**2, t])
        x0 = np.array([0., 0.])
        check(u, x0, xtrue)

        # test with linear input
        u = t
        xtrue = np.array([1./6. * t**3, 0.5 * t**2])
        check(u, x0, xtrue)

Example 28

Project: RLScore Source File: steepest_descent_mmc.py
Function: init
    def __init__(self, X, regparam=1.0, number_of_clusters=2, kernel='LinearKernel', basis_vectors=None, Y = None, fixed_indices=None, callback=None,  **kwargs):
        kwargs['X'] = X
        kwargs['kernel'] = kernel
        if basis_vectors is not None:
            kwargs['basis_vectors'] = basis_vectors
        self.svdad = adapter.createSVDAdapter(**kwargs)
        self.svals = np.mat(self.svdad.svals)
        self.svecs = np.mat(self.svdad.rsvecs)
        self.callbackfun = callback
        self.regparam = regparam
        self.constraint = 0
        self.labelcount = int(number_of_clusters)
         
        if self.labelcount == 2:
            self.oneclass = True
        else:
            self.oneclass = False
        if Y is not None:
            Y_orig = array_tools.as_array(Y)
            #if Y_orig.shape[1] == 1:
            if len(Y_orig.shape) == 1:
                self.Y = np.zeros((Y_orig.shape[0], 2))
                self.Y[:, 0] = Y_orig
                self.Y[:, 1] = - Y_orig
                self.oneclass = True
            else:
                self.Y = Y_orig.copy()
                self.oneclass = False
            for i in range(self.Y.shape[0]):
                largestind = 0
                largestval = self.Y[i, 0]
                for j in range(self.Y.shape[1]):
                    if self.Y[i, j] > largestval:
                        largestind = j
                        largestval = self.Y[i, j]
                    self.Y[i, j] = -1.
                self.Y[i, largestind] = 1.
        else:
            size = self.svecs.shape[0]
            ysize = self.labelcount
            if self.labelcount is None: self.labelcount = 2
            self.Y = RandomLabelSource(size, ysize).readLabels()
        self.size = self.Y.shape[0]
        self.labelcount = self.Y.shape[1]
        self.classvec = - np.ones((self.size), dtype = np.int32)
        self.classcounts = np.zeros((self.labelcount), dtype = np.int32)
        for i in range(self.size):
            clazzind = 0
            largestlabel = self.Y[i, 0]
            for j in range(self.labelcount):
                if self.Y[i, j] > largestlabel:
                    largestlabel = self.Y[i, j]
                    clazzind = j
            self.classvec[i] = clazzind
            self.classcounts[clazzind] = self.classcounts[clazzind] + 1
        
        self.svecs_list = []
        for i in range(self.size):
            self.svecs_list.append(self.svecs[i].T)
        self.fixedindices = []
        if fixed_indices is not None:
            self.fixedindices = fixed_indices
        else:
            self.fixedindices = []
        self.results = {}
        self.solve(self.regparam)

Example 29

Project: python-control Source File: yottalab.py
def dlqr(*args, **keywords):
    """Linear quadratic regulator design for discrete systems

    Usage
    =====
    [K, S, E] = dlqr(A, B, Q, R, [N])
    [K, S, E] = dlqr(sys, Q, R, [N])

    The dlqr() function computes the optimal state feedback controller
    that minimizes the quadratic cost

        J = \sum_0^\infty x' Q x + u' R u + 2 x' N u

    Inputs
    ------
    A, B: 2-d arrays with dynamics and input matrices
    sys: linear I/O system 
    Q, R: 2-d array with state and input weight matrices
    N: optional 2-d array with cross weight matrix

    Outputs
    -------
    K: 2-d array with state feedback gains
    S: 2-d array with solution to Riccati equation
    E: 1-d array with eigenvalues of the closed loop system
    """

    # 
    # Process the arguments and figure out what inputs we received
    #
    
    # Get the system description
    if (len(args) < 3):
        raise ControlArgument("not enough input arguments")

    elif (ctrlutil.issys(args[0])):
        # We were passed a system as the first argument; extract A and B
        A = array(args[0].A, ndmin=2, dtype=float);
        B = array(args[0].B, ndmin=2, dtype=float);
        index = 1;
        if args[0].dt==0.0:
            print "dlqr works only for discrete systems!"
            return
    else:
        # Arguments should be A and B matrices
        A = array(args[0], ndmin=2, dtype=float);
        B = array(args[1], ndmin=2, dtype=float);
        index = 2;

    # Get the weighting matrices (converting to matrices, if needed)
    Q = array(args[index], ndmin=2, dtype=float);
    R = array(args[index+1], ndmin=2, dtype=float);
    if (len(args) > index + 2): 
        N = array(args[index+2], ndmin=2, dtype=float);
        Nflag = 1;
    else:
        N = zeros((Q.shape[0], R.shape[1]));
        Nflag = 0;

    # Check dimensions for consistency
    nstates = B.shape[0];
    ninputs = B.shape[1];
    if (A.shape[0] != nstates or A.shape[1] != nstates):
        raise ControlDimension("inconsistent system dimensions")

    elif (Q.shape[0] != nstates or Q.shape[1] != nstates or
          R.shape[0] != ninputs or R.shape[1] != ninputs or
          N.shape[0] != nstates or N.shape[1] != ninputs):
        raise ControlDimension("incorrect weighting matrix dimensions")

    if Nflag==1:
        Ao=A-B*inv(R)*N.T
        Qo=Q-N*inv(R)*N.T
    else:
        Ao=A
        Qo=Q
    
    #Solve the riccati equation
    (X,L,G) = dare(Ao,B,Qo,R)
#    X = bb_dare(Ao,B,Qo,R)

    # Now compute the return value
    Phi=mat(A)
    H=mat(B)
    K=inv(H.T*X*H+R)*(H.T*X*Phi+N.T)
    L=eig(Phi-H*K)
    return K,X,L

Example 30

Project: RLScore Source File: greedy_rls.py
    def _solve_new(self, regparam, floattype):
        
        #Legacy code. Works only with a single output but can work with given performance measures and is faster than _solve_bu
        
        if not self.Y.shape[1] == 1:
            raise Exception('This variation of GreedyRLS supports only one output at a time. The output matrix is now of shape ' + str(self.Y.shape) + '.')
        self.regparam = regparam
        X = self.X
        Y = np.mat(self.Y, dtype=floattype)
        
        bias_slice = np.sqrt(self.bias)*np.mat(np.ones((1,X.shape[1]), dtype=floattype))
        
        tsize = self.size
        fsize = X.shape[0]
        assert X.shape[1] == tsize
        self.A = np.mat(np.zeros((fsize,1),dtype=floattype))
        
        rp = regparam
        rpinv = 1. / rp
        
        
        #Biaz
        cv = np.sqrt(self.bias)*np.mat(np.ones((1, tsize), dtype=floattype))
        ca = np.mat(rpinv * (1. / (1. + cv * rpinv * cv.T)) * (cv * rpinv), dtype=floattype)
        
        
        self.dualvec = rpinv * Y - cv.T * rpinv * (1. / (1. + cv * rpinv * cv.T)) * (cv * rpinv * Y)
        
        GXT = cv.T * np.mat((rpinv * (1. / (1. + cv * rpinv * cv.T)) * (cv * rpinv)) * X.T, dtype=floattype)
        tempmatrix = np.mat(np.zeros(X.T.shape, dtype=floattype))
        np.multiply(X.T, rpinv, tempmatrix)
        #tempmatrix = rpinv * X.T
        np.subtract(tempmatrix, GXT, GXT)
        
        diagG = []
        for i in range(tsize):
            diagGi = rpinv - cv.T[i, 0] * ca[0, i]
            diagG.append(diagGi)
        diagG = np.mat(diagG, dtype=floattype).T
        
        self.selected = []
        self.performances = []
        currentfcount = 0
        
        temp2 = np.mat(np.zeros(tempmatrix.shape, dtype=floattype))
        
        while currentfcount < self.desiredfcount:
            
            np.multiply(X.T, GXT, tempmatrix)
            XGXTdiag = np.sum(tempmatrix, axis = 0)
            
            XGXTdiag = 1. / (1. + XGXTdiag)
            np.multiply(GXT, XGXTdiag, tempmatrix)
            
            tempvec1 = np.multiply((X * self.dualvec).T, XGXTdiag)
            np.multiply(GXT, tempvec1, temp2)
            np.subtract(self.dualvec, temp2, temp2)
            
            np.multiply(tempmatrix, GXT, tempmatrix)
            np.subtract(diagG, tempmatrix, tempmatrix)
            np.divide(1, tempmatrix, tempmatrix)
            np.multiply(tempmatrix, temp2, tempmatrix)
            
            
            if not self.measure is None:
                np.subtract(Y, tempmatrix, tempmatrix)
                np.multiply(temp2, 0, temp2)
                np.add(temp2, Y, temp2)
                looperf = self.measure.multiTaskPerformance(temp2, tempmatrix)
                looperf = np.mat(looperf, dtype=floattype)
                if self.measure.iserror:
                    looperf[0, self.selected] = float('inf')
                    bestcind = np.argmin(looperf)
                    self.bestlooperf = np.amin(looperf)
                else:
                    looperf[0, self.selected] = - float('inf')
                    bestcind = np.argmax(looperf)
                    self.bestlooperf = np.amax(looperf)
            else:
                np.multiply(tempmatrix, tempmatrix, temp2)
                looperf = np.sum(temp2, axis = 0)
                looperf[0, self.selected] = float('inf')
                bestcind = np.argmin(looperf)
                self.bestlooperf = np.amin(looperf)
                self.loo_predictions = Y - tempmatrix[:, bestcind]
            
            self.looperf = looperf   #Needed in test_GreedyRLS module
            
            self.performances.append(self.bestlooperf)
            cv = X[bestcind]
            GXT_bci = GXT[:, bestcind]
            ca = GXT_bci * (1. / (1. + cv * GXT_bci))
            self.dualvec = self.dualvec - ca * (cv * self.dualvec)
            diagG = diagG - np.multiply(ca, GXT_bci)
            np.multiply(tempmatrix, 0, tempmatrix)
            np.add(tempmatrix, ca, tempmatrix)
            tempvec1 = cv * GXT
            np.multiply(tempmatrix, tempvec1, tempmatrix)
            np.subtract(GXT, tempmatrix, GXT)
            self.selected.append(bestcind)
            currentfcount += 1
            
            #Linear predictor with bias
            self.A[self.selected] = X[self.selected] * self.dualvec
            self.b = bias_slice * self.dualvec * np.sqrt(self.bias)
            self.predictor = predictor.LinearPredictor(self.A, self.b)
            
            if not self.callbackfun is None:
                self.callbackfun.callback(self)
        if not self.callbackfun is None:
            self.callbackfun.finished(self)
        self.A[self.selected] = X[self.selected] * self.dualvec
        self.b = bias_slice * self.dualvec * np.sqrt(self.bias)
        self.results[SELECTED_FEATURES] = self.selected
        self.results[GREEDYRLS_LOO_PERFORMANCES] = self.performances
        #self.results['predictor'] = self.getModel()
        self.predictor = predictor.LinearPredictor(self.A, self.b)

Example 31

Project: python-control Source File: yottalab.py
def d2c(sys,method='zoh'):
    """Continous to discrete conversion with ZOH method

    Call:
    sysc=c2d(sys,method='log')

    Parameters
    ----------
    sys :   System in statespace or Tf form 
    method: 'zoh' or 'bi'

    Returns
    -------
    sysc: continous system ss or tf
    

    """
    flag = 0
    if isinstance(sys, TransferFunction):
        sys=tf2ss(sys)
        flag=1

    a=sys.A
    b=sys.B
    c=sys.C
    d=sys.D
    Ts=sys.dt
    n=shape(a)[0]
    nb=shape(b)[1]
    nc=shape(c)[0]
    tol=1e-12
    
    if method=='zoh':
        if n==1:
            if b[0,0]==1:
                A=0
                B=b/sys.dt
                C=c
                D=d
        else:
            tmp1=hstack((a,b))
            tmp2=hstack((zeros((nb,n)),eye(nb)))
            tmp=vstack((tmp1,tmp2))
            s=logm(tmp)
            s=s/Ts
            if norm(imag(s),inf) > sqrt(sp.finfo(float).eps):
                print "Warning: accuracy may be poor"
            s=real(s)
            A=s[0:n,0:n]
            B=s[0:n,n:n+nb]
            C=c
            D=d
    elif method=='foh':
        a=mat(a)
        b=mat(b)
        c=mat(c)
        d=mat(d)
        Id = mat(eye(n))
        A = logm(a)/Ts
        A = real(around(A,12))
        Amat = mat(A)
        B = (a-Id)**(-2)*Amat**2*b*Ts
        B = real(around(B,12))
        Bmat = mat(B)
        C = c
        D = d - C*(Amat**(-2)/Ts*(a-Id)-Amat**(-1))*Bmat
        D = real(around(D,12))
    elif method=='bi':
        a=mat(a)
        b=mat(b)
        c=mat(c)
        d=mat(d)
        poles=eigvals(a)
        if any(abs(poles-1)<200*sp.finfo(float).eps):
            print "d2c: some poles very close to one. May get bad results."
        
        I=mat(eye(n,n))
        tk = 2 / sqrt (Ts)
        A = (2/Ts)*(a-I)*inv(a+I)
        iab = inv(I+a)*b
        B = tk*iab
        C = tk*(c*inv(I+a))
        D = d- (c*iab)
    else:
        print "Method not supported"
        return
    
    sysc=StateSpace(A,B,C,D)
    if flag==1:
        sysc=ss2tf(sysc)
    return sysc

Example 32

Project: RLScore Source File: interactive_rls_classifier.py
Function: init
    def __init__(self, X, regparam=1.0, number_of_clusters=2, kernel='LinearKernel', basis_vectors=None, Y = None, fixed_indices=None, callback=None,  **kwargs):
        kwargs['X'] = X 
        kwargs['kernel'] = kernel
        if basis_vectors is not None:
            kwargs['basis_vectors'] = basis_vectors
        self.svdad = adapter.createSVDAdapter(**kwargs)
        self.svals = np.mat(self.svdad.svals)
        self.svecs = np.mat(self.svdad.rsvecs)
        self.callbackfun = callback
        self.regparam = regparam
        self.constraint = 0
        self.labelcount = number_of_clusters
        self.size = X.shape[0] 
        #if self.labelcount == 2:
        #    self.oneclass = True
        #else:
        #    self.oneclass = False
        
        if Y is None:
            self.classvec = np.zeros(self.size, np.int)
        else:
            self.classvec = Y
        #self.size = self.classvec.shape[0]
        self.Y = -np.ones((self.size, self.labelcount))
        self.classcounts = np.zeros((self.labelcount), dtype = np.int32)
        for i in range(self.size):
            clazzind = self.classvec[i]
            self.Y[i, clazzind] = 1
            self.classcounts[clazzind] = self.classcounts[clazzind] + 1
        
        self.fixedindices = []
        if fixed_indices is not None:
            self.fixedindices = fixed_indices
        self.train()

Example 33

Project: RLScore Source File: measure_utilities.py
Function: wrapper
def wrapper(measure, Y, Y_predicted, qids):
    Y = np.mat(Y)
    Y_predicted = np.mat(Y_predicted)
    qids_perfs = []
    for inds in qids:
        Y_sub = Y[inds]
        Y_predicted_sub = Y_predicted[inds]
        perfs = measure.getPerformance(Y_sub, Y_predicted_sub)
        qids_perfs.append(perfs)
    #quite a bit juggling follows to handle the fact, that nans encode
    #queries for which performance is undefined (happens sometimes
    #in ranking
    #
    #count the number of non-nan values in each column
    perfs = np.vstack(qids_perfs)
    normalizers = np.isnan(perfs)
    normalizers = np.logical_not(normalizers)
    normalizers = np.sum(normalizers, axis=0)
    normalizers = np.where(normalizers>0,normalizers,np.nan)
    #turns nans into zeroes
    perfs = np.nan_to_num(perfs)
    perfs = np.sum(perfs, axis=0)
    perfs = perfs/normalizers
    return perfs

Example 34

Project: python-control Source File: yottalab.py
def red_obs(sys,T,poles):
    """Reduced order observer of the system sys

    Call:
    obs=red_obs(sys,T,poles)

    Parameters
    ----------
    sys : System in State Space form
    T: Complement matrix
    poles: desired observer poles

    Returns
    -------
    obs: ss
    Reduced order Observer

    """
    if isinstance(sys, TransferFunction):
        "System must be in state space form"
        return
    a=mat(sys.A)
    b=mat(sys.B)
    c=mat(sys.C)
    d=mat(sys.D)
    T=mat(T)
    P=mat(vstack((c,T)))
    invP=inv(P)
    AA=P*a*invP
    ny=shape(c)[0]
    nx=shape(a)[0]
    nu=shape(b)[1]

    A11=AA[0:ny,0:ny]
    A12=AA[0:ny,ny:nx]
    A21=AA[ny:nx,0:ny]
    A22=AA[ny:nx,ny:nx]

    L1=placep(A22.T,A12.T,poles)
    L1=mat(L1).T

    nn=nx-ny

    tmp1=mat(hstack((-L1,eye(nn,nn))))
    tmp2=mat(vstack((zeros((ny,nn)),eye(nn,nn))))
    Ar=tmp1*P*a*invP*tmp2
 
    tmp3=vstack((eye(ny,ny),L1))
    tmp3=mat(hstack((P*b,P*a*invP*tmp3)))
    tmp4=hstack((eye(nu,nu),zeros((nu,ny))))
    tmp5=hstack((-d,eye(ny,ny)))
    tmp4=mat(vstack((tmp4,tmp5)))

    Br=tmp1*tmp3*tmp4

    Cr=invP*tmp2

    tmp5=hstack((zeros((ny,nu)),eye(ny,ny)))
    tmp6=hstack((zeros((nn,nu)),L1))
    tmp5=mat(vstack((tmp5,tmp6)))
    Dr=invP*tmp5*tmp4
    
    obs=StateSpace(Ar,Br,Cr,Dr,sys.dt)
    return obs

Example 35

Project: python-control Source File: yottalab.py
def full_obs(sys,poles):
    """Full order observer of the system sys

    Call:
    obs=full_obs(sys,poles)

    Parameters
    ----------
    sys : System in State Space form
    poles: desired observer poles

    Returns
    -------
    obs: ss
    Observer

    """
    if isinstance(sys, TransferFunction):
        "System must be in state space form"
        return
    a=mat(sys.A)
    b=mat(sys.B)
    c=mat(sys.C)
    d=mat(sys.D)
    L=placep(a.T,c.T,poles)
    L=mat(L).T
    Ao=a-L*c
    Bo=hstack((b-L*d,L))
    n=shape(Ao)
    m=shape(Bo)
    Co=eye(n[0],n[1])
    Do=zeros((n[0],m[1]))
    obs=StateSpace(Ao,Bo,Co,Do,sys.dt)
    return obs

Example 36

Project: RLScore Source File: global_rankrls.py
Function: hold_out
    def holdout(self, indices):
        """Computes hold-out predictions for a trained RankRLS
        
        Parameters
        ----------
        indices : list of indices, shape = [n_hsamples]
            list of indices of training examples belonging to the set for which the
            hold-out predictions are calculated. The list can not be empty.

        Returns
        -------
        F : array, shape = [n_hsamples, n_labels]
            holdout predictions
            
        Notes
        -----
    
        The algorithm is a modification of the ones published in [1,2] for the regular RLS method.
        
        References
        ----------
        
        [1] Tapio Pahikkala, Jorma Boberg, and Tapio Salakoski.
        Fast n-Fold Cross-Validation for Regularized Least-Squares.
        Proceedings of the Ninth Scandinavian Conference on Artificial Intelligence,
        83-90, Otamedia Oy, 2006.
        
        [2] Tapio Pahikkala, Hanna Suominen, and Jorma Boberg.
        Efficient cross-validation for kernelized least-squares regression with sparse basis expansions.
        Machine Learning, 87(3):381--407, June 2012.     
        """
        
        indices = array_tools.as_index_list(indices, self.Y.shape[0])
        
        if len(indices) != len(np.unique(indices)):
            raise IndexError('Hold-out can have each index only once.')
        
        Y = self.Y
        m = self.size
        
        evals, V = self.evals, self.svecs
        
        #results = []
        
        C = np.mat(np.zeros((self.size, 3), dtype = np.float64))
        onevec = np.mat(np.ones((self.size, 1), dtype = np.float64))
        for i in range(self.size):
            C[i, 0] = 1.
        
        
        VTY = V.T * Y
        VTC = V.T * onevec
        CTY = onevec.T * Y
        
        holen = len(indices)
        
        newevals = multiply(evals, 1. / ((m - holen) * evals + self.regparam))
        
        R = np.mat(np.zeros((holen, holen + 1), dtype = np.float64))
        for i in range(len(indices)):
            R[i, 0] = -1.
            R[i, i + 1] = sqrt(self.size - float(holen))
        
        Vho = V[indices]
        Vhov = multiply(Vho, newevals)
        Ghoho = Vhov * Vho.T
        GCho = Vhov * VTC
        GBho = Ghoho * R
        for i in range(len(indices)):
            GBho[i, 0] += GCho[i, 0]
        
        CTGC = multiply(VTC.T, newevals) * VTC
        RTGCho = R.T * GCho
        
        BTGB = R.T * Ghoho * R
        for i in range(len(indices) + 1):
            BTGB[i, 0] += RTGCho[i, 0]
            BTGB[0, i] += RTGCho[i, 0]
        BTGB[0, 0] += CTGC[0, 0]
        
        BTY = R.T * Y[indices]
        BTY[0] = BTY[0] + CTY[0]
        
        GDYho = Vhov * (self.size - holen) * VTY
        GLYho = GDYho - GBho * BTY
        
        CTGDY = multiply(VTC.T, newevals) * (self.size - holen) * VTY
        BTGLY = R.T * GDYho - BTGB * BTY
        BTGLY[0] = BTGLY[0] + CTGDY[0]
        
        F = GLYho - GBho * la.inv(-mat(eye(holen + 1)) + BTGB) * BTGLY
        
        #results.append(F)
        #return results
        F = np.squeeze(np.array(F))
        return F

Example 37

Project: RLScore Source File: kron_rls.py
Function: init
    def __init__(self, **kwargs):
        Y = kwargs["Y"]
        Y = array_tools.as_2d_array(Y)
        Y = np.mat(Y)
        if kwargs.has_key('K1'):
            K1 = np.mat(kwargs['K1'])
            K2 = np.mat(kwargs['K2'])
            Y = Y.reshape((K1.shape[0], K2.shape[0]), order = 'F')
            self.K1, self.K2 = K1, K2
            self.kernelmode = True
        else:
            X1 = np.mat(kwargs['X1'])
            X2 = np.mat(kwargs['X2'])
            Y = Y.reshape((X1.shape[0], X2.shape[0]), order = 'F')
            self.X1, self.X2 = X1, X2
            self.kernelmode = False
        self.Y = Y
        if kwargs.has_key("regparam"):
            self.regparam = kwargs["regparam"]
        else:
            self.regparam = 1.
        self.trained = False
        self.solve(self.regparam)

Example 38

Project: python-control Source File: yottalab.py
def comp_form(sys,obs,K):
    """Compact form Conroller+Observer

    Call:
    contr=comp_form(sys,obs,K)

    Parameters
    ----------
    sys : System in State Space form
    obs : Observer in State Space form
    K: State feedback gains

    Returns
    -------
    contr: ss
    Controller

    """
    nx=shape(sys.A)[0]
    ny=shape(sys.C)[0]
    nu=shape(sys.B)[1]
    no=shape(obs.A)[0]

    Bu=mat(obs.B[:,0:nu])
    By=mat(obs.B[:,nu:])
    Du=mat(obs.D[:,0:nu])
    Dy=mat(obs.D[:,nu:])

    X=inv(eye(nu,nu)+K*Du)

    Ac = mat(obs.A)-Bu*X*K*mat(obs.C);
    Bc = hstack((Bu*X,By-Bu*X*K*Dy))
    Cc = -X*K*mat(obs.C);
    Dc = hstack((X,-X*K*Dy))
    contr = StateSpace(Ac,Bc,Cc,Dc,sys.dt)
    return contr

Example 39

Project: RLScore Source File: kron_rls.py
    def solve_linear_conditional_ranking(self, regparam):
        """Trains conditional ranking KronRLS, that ranks objects from
        domain 2 against objects from domain 1.
               
        Parameters
        ----------
        regparam : float, optional
            regularization parameter, regparam > 0 (default=1.0)
            
        Notes
        -----
        Minimizes RankRLS type of loss. Currently only linear kernel
        supported. Including the code here is a hack, this should
        probably be implemented as an independent learner.
        """
        self.regparam = regparam
        X1, X2 = self.X1, self.X2
        Y = self.Y.reshape((X1.shape[0], X2.shape[0]), order = 'F')
        
        V, svals1, rsvecs1 = linalg.svd_economy_sized(X1)
        svals1 = np.mat(svals1)
        self.svals1 = svals1.T
        self.evals1 = np.multiply(self.svals1, self.svals1)
        self.V = V
        self.rsvecs1 = np.mat(rsvecs1)
        
        qlen = X2.shape[0]
        onevec = (1. / np.math.sqrt(qlen)) * np.mat(np.ones((qlen, 1)))
        C = np.mat(np.eye(qlen)) - onevec * onevec.T
        
        U, svals2, rsvecs2 = linalg.svd_economy_sized(C * X2)
        svals2 = np.mat(svals2)
        self.svals2 = svals2.T
        self.evals2 = np.multiply(self.svals2, self.svals2)
        self.U = U
        self.rsvecs2 = np.mat(rsvecs2)
        
        self.VTYU = V.T * Y * C * U
        
        kronsvals = self.svals1 * self.svals2.T
        
        newevals = np.divide(kronsvals, np.multiply(kronsvals, kronsvals) + regparam)
        self.W = np.multiply(self.VTYU, newevals)
        self.W = self.rsvecs1.T * self.W * self.rsvecs2
        self.predictor = LinearPairwisePredictor(np.array(self.W))

Example 40

Project: RLScore Source File: mmc.py
Function: init
    def __init__(self, X, regparam=1.0, number_of_clusters=2, kernel='LinearKernel', basis_vectors=None, Y = None, fixed_indices=None, callback=None,  **kwargs):
        kwargs['X'] = X 
        kwargs['kernel'] = kernel
        if basis_vectors is not None:
            kwargs['basis_vectors'] = basis_vectors
        self.svdad = adapter.createSVDAdapter(**kwargs)
        self.svals = np.mat(self.svdad.svals)
        self.svecs = self.svdad.rsvecs
        self.regparam = regparam
        self.constraint = 0
        #if not kwargs.has_key('number_of_clusters'):
        #    raise Exception("Parameter 'number_of_clusters' must be given.")
        self.labelcount = number_of_clusters
        if self.labelcount == 2:
            self.oneclass = True
        else:
            self.oneclass = False
        self.callbackfun = callback
        if Y is not None:
            Y_orig = array_tools.as_array(Y)
            #if Y_orig.shape[1] == 1:
            if len(Y_orig.shape) == 1:
                self.Y = np.zeros((Y_orig.shape[0], 2))
                self.Y[:, 0] = Y_orig
                self.Y[:, 1] = - Y_orig
                self.oneclass = True
            else:
                self.Y = Y_orig.copy()
                self.oneclass = False
            for i in range(self.Y.shape[0]):
                largestind = 0
                largestval = self.Y[i, 0]
                for j in range(self.Y.shape[1]):
                    if self.Y[i, j] > largestval:
                        largestind = j
                        largestval = self.Y[i, j]
                    self.Y[i, j] = -1.
                self.Y[i, largestind] = 1.
        else:
            size = self.svecs.shape[0]
            ysize = self.labelcount
            if self.labelcount is None: self.labelcount = 2
            self.Y = RandomLabelSource(size, ysize).readLabels()
        self.size = self.Y.shape[0]
        self.labelcount = self.Y.shape[1]
        self.classvec = - np.ones((self.size), dtype = np.int32)
        self.classcounts = np.zeros((self.labelcount), dtype = np.int32)
        for i in range(self.size):
            clazzind = 0
            largestlabel = self.Y[i, 0]
            for j in range(self.labelcount):
                if self.Y[i, j] > largestlabel:
                    largestlabel = self.Y[i, j]
                    clazzind = j
            self.classvec[i] = clazzind
            self.classcounts[clazzind] = self.classcounts[clazzind] + 1
        
        self.svecs_list = []
        for i in range(self.size):
            self.svecs_list.append(self.svecs[i].T)
        self.fixedindices = []
        if fixed_indices is not None:
            self.fixedindices = fixed_indices
        else:
            self.fixedindices = []
        self.results = {}
        self.solve(self.regparam)

Example 41

Project: RLScore Source File: test_reduced_set_approximation.py
Function: testrls
    def testRLS(self):
        
        print
        print
        print
        print
        print "Testing the cross-validation routines of the RLS module."
        print
        print
        
        m, n = 100, 300
        Xtrain = random.rand(m, n)
        Y = mat(random.rand(m, 1))
        basis_vectors = [0,3,7,8]
        
        def complement(indices, m):
            compl = range(m)
            for ind in indices:
                compl.remove(ind)
            return compl
        
        #hoindices = [45, 50, 55]
        hoindices = [0, 1, 2]
        hocompl = complement(hoindices, m)
        
        bk = GaussianKernel(**{'X':Xtrain[basis_vectors], 'gamma':0.001})
        
        rpool = {}
        rpool['X'] = Xtrain
        bk2 = GaussianKernel(**{'X':Xtrain, 'gamma':0.001})
        K = np.mat(bk2.getKM(Xtrain))
        
        Yho = Y[hocompl]
        
        
        rpool = {}
        rpool['Y'] = Y
        rpool['X'] = Xtrain
        rpool['basis_vectors'] = Xtrain[basis_vectors]
        
        Xhocompl = Xtrain[hocompl]
        testX = Xtrain[hoindices]
        
        rpool = {}
        rpool['Y'] = Yho
        rpool['X'] = Xhocompl
        rpool["kernel"] = "RsetKernel"
        rpool["base_kernel"] = bk
        rpool["basis_features"] = Xtrain[basis_vectors]
        #rk = RsetKernel(**{'base_kernel':bk, 'basis_features':Xtrain[basis_vectors], 'X':Xhocompl})
        dualrls_naive = RLS(**rpool)
        
        rpool = {}
        rpool['Y'] = Yho
        rpool['X'] = Xhocompl
        
        rsaK = K[:, basis_vectors] * la.inv(K[ix_(basis_vectors, basis_vectors)]) * K[basis_vectors]
        rsaKho = rsaK[ix_(hocompl, hocompl)]
        rsa_testkm = rsaK[ix_(hocompl, hoindices)]
        loglambdas = range(-5, 5)
        for j in range(0, len(loglambdas)):
            regparam = 2. ** loglambdas[j]
            print
            print "Regparam 2^%1d" % loglambdas[j]
            
            print (rsa_testkm.T * la.inv(rsaKho + regparam * eye(rsaKho.shape[0])) * Yho).T, 'Dumb HO (dual)'
            dumbho = np.squeeze(np.array(rsa_testkm.T * la.inv(rsaKho + regparam * eye(rsaKho.shape[0])) * Yho))
            
            dualrls_naive.solve(regparam)
            predho1 = np.squeeze(dualrls_naive.predictor.predict(testX))
            print predho1.T, 'Naive HO (dual)'
            
            #dualrls.solve(regparam)
            #predho2 = np.squeeze(dualrls.computeHO(hoindices))
            #print predho2.T, 'Fast HO (dual)'
            
            for predho in [dumbho, predho1]:#, predho2]:
                self.assertEqual(dumbho.shape, predho.shape)
                for row in range(predho.shape[0]):
                    #for col in range(predho.shape[1]):
                    #    self.assertAlmostEqual(dumbho[row,col],predho[row,col])
                        self.assertAlmostEqual(dumbho[row],predho[row])

Example 42

Project: RLScore Source File: global_rankrls.py
Function: reference
    def _reference(self, pairs):
        
        evals, evecs = self.evals, self.svecs
        Y = self.Y
        m = self.size
        
        
        results = []
        
        D = mat(zeros((self.size, 1), dtype=float64))
        C = mat(zeros((self.size, 3), dtype=float64))
        for i in range(self.size):
            D[i, 0] = self.size - 2.
            C[i, 0] = 1.
        lamb = self.regparam
        
        G = evecs * multiply(multiply(evals, 1. / ((m - 2.) * evals + lamb)).T, evecs.T)
        
        
        DY = multiply(D, Y)
        GDY = G * DY
        GC = G * C
        CTG = GC.T
        CT = C.T
        CTGC = CT * GC
        CTY = CT * Y
        CTGDY = CT * GDY
        
        minusI3 = -mat(eye(3))
        for i, j in pairs:
            hoinds = [i, j]
            
            R = mat(zeros((2, 3), dtype=float64))
            R[0, 0] = -1.
            R[1, 0] = -1.
            R[0, 1] = sqrt(self.size - 2.)
            R[1, 2] = sqrt(self.size - 2.)
            
            GBho = GC[hoinds] + G[np.ix_(hoinds, hoinds)] * R
            
            BTGB = CTGC \
                + R.T * GC[hoinds] \
                + CTG[:, hoinds] * R \
                + R.T * G[np.ix_(hoinds, hoinds)] * R
            
            BTY = CTY + R.T * Y[hoinds]
            
            GLYho = GDY[hoinds] - GBho * BTY
            BTGLY = CTGDY + R.T * GDY[hoinds] - BTGB * BTY
            
            F = GLYho - GBho * la.inv(minusI3 + BTGB) * BTGLY
            
            results.append(F)
        return results

Example 43

Project: python-control Source File: modelsimp.py
Function: markov
def markov(Y, U, M):
    """
    Calculate the first `M` Markov parameters [D CB CAB ...]
    from input `U`, output `Y`.

    Parameters
    ----------
    Y: array_like
        Output data
    U: array_like
        Input data
    M: integer
        Number of Markov parameters to output

    Returns
    -------
    H: matrix
        First M Markov parameters

    Notes
    -----
    Currently only works for SISO

    Examples
    --------
    >>> H = markov(Y, U, M)
    """

    # Convert input parameters to matrices (if they aren't already)
    Ymat = np.mat(Y)
    Umat = np.mat(U)
    n = np.size(U)

    # Construct a matrix of control inputs to invert
    UU = Umat
    for i in range(1, M-1):
        newCol = np.vstack((0, UU[0:n-1,i-2]))
        UU = np.hstack((UU, newCol))
    Ulast = np.vstack((0, UU[0:n-1,M-2]))
    for i in range(n-1,0,-1):
        Ulast[i] = np.sum(Ulast[0:i-1])
    UU = np.hstack((UU, Ulast))

    # Invert and solve for Markov parameters
    H = UU.I
    H = np.dot(H, Y)

    return H

Example 44

Project: python-control Source File: statefbk.py
def acker(A, B, poles):
    """Pole placement using Ackermann method

    Call:
    K = acker(A, B, poles)

    Parameters
    ----------
    A, B : 2-d arrays
        State and input matrix of the system
    poles: 1-d list
        Desired eigenvalue locations

    Returns
    -------
    K: matrix
        Gains such that A - B K has given eigenvalues

    """
    # Convert the inputs to matrices
    a = np.mat(A)
    b = np.mat(B)

    # Make sure the system is controllable
    ct = ctrb(A, B)
    if sp.linalg.det(ct) == 0:
        raise ValueError("System not reachable; pole placement invalid")

    # Compute the desired characteristic polynomial
    p = np.real(np.poly(poles))

    # Place the poles using Ackermann's method
    n = np.size(p)
    pmat = p[n-1]*a**0
    for i in np.arange(1,n):
        pmat = pmat + p[n-i-1]*a**i
    K = sp.linalg.inv(ct) * pmat

    K = K[-1][:]                # Extract the last row
    return K

Example 45

Project: python-control Source File: statefbk.py
def ctrb(A,B):
    """Controllabilty matrix

    Parameters
    ----------
    A, B: array_like or string
        Dynamics and input matrix of the system

    Returns
    -------
    C: matrix
        Controllability matrix

    Examples
    --------
    >>> C = ctrb(A, B)

    """

    # Convert input parameters to matrices (if they aren't already)
    amat = np.mat(A)
    bmat = np.mat(B)
    n = np.shape(amat)[0]
    # Construct the controllability matrix
    ctrb = bmat
    for i in range(1, n):
        ctrb = np.hstack((ctrb, amat**i*bmat))
    return ctrb

Example 46

Project: RLScore Source File: gaussian_kernel.py
Function: getkm
    def getKM(self, X):
        """Returns the kernel matrix between the basis vectors and X.
        
        Parameters
        ----------
        X : {array-like, sparse matrix}, shape = [n_samples, n_features]
        
        Returns
        -------
        K : array, shape = [n_samples, n_bvectors]
            kernel matrix
        """
        X = array_tools.as_2d_array(X, True)
        test_X = X 
        if sp.issparse(test_X):
            test_X = array_tools.spmat_resize(test_X, self.train_X.shape[1])
        else:
            test_X = array_tools.as_dense_matrix(test_X)
        gamma = self.gamma
        m = self.train_X.shape[0]
        n = test_X.shape[0]
        #The Gaussian kernel matrix is constructed from a linear kernel matrix
        linkm = self.train_X * test_X.T
        linkm = array_tools.as_dense_matrix(linkm)
        if sp.issparse(test_X):
            test_norms = ((test_X.T.multiply(test_X.T)).sum(axis=0)).T
        else:
            test_norms = (np.multiply(test_X.T, test_X.T).sum(axis=0)).T
        K = mat(np.ones((m, 1), dtype = float64)) * test_norms.T
        K = K + self.train_norms * mat(np.ones((1, n), dtype = float64))
        K = K - 2 * linkm
        K = - gamma * K
        K = np.exp(K)
        return K.A.T

Example 47

Project: RLScore Source File: global_rankrls.py
Function: solve
    def solve(self, regparam=1.0):
        """Re-trains RankRLS for the given regparam
               
        Parameters
        ----------
        regparam : float, optional
            regularization parameter, regparam > 0 (default=1.0)
            
        Notes
        -----
    
        Computational complexity of re-training:
        m = n_samples, d = n_features, l = n_labels, b = n_bvectors
        
        O(lm^2): basic case
        
        O(lmd): Linear Kernel, d < m
        
        O(lmb): Sparse approximation with basis vectors 
             
        """
        if not hasattr(self, "multiplyright"):
            
            #Eigenvalues of the kernel matrix
            self.evals = multiply(self.svals, self.svals)
            
            #Temporary variables
            ssvecs = multiply(self.svecs, self.svals)
            J = mat(ones((self.size, 1), dtype=float64))
            
            #These are cached for later use in solve function
            self.lowrankchange = ssvecs.T * J[range(ssvecs.shape[0])]
            self.multipleright = ssvecs.T * (self.size * self.Y - J * (J.T * self.Y))
        
        self.regparam = regparam
        
        #Compute the eigenvalues determined by the given regularization parameter
        neweigvals = 1. / (self.size * self.evals + regparam)
        
        P = self.lowrankchange
        nP = multiply(neweigvals.T, P)
        fastinv = 1. / (-1. + P.T * nP)
        self.A = self.svecs * multiply(1. / self.svals.T, \
            (multiply(neweigvals.T, self.multipleright) \
            - nP * (fastinv * (nP.T * self.multipleright))))
        self.predictor = self.svdad.createModel(self)

Example 48

Project: RLScore Source File: greedy_rls.py
Function: init
    def __init__(self, X, Y, subsetsize, regparam = 1.0, bias=1.0, callbackfun=None, **kwargs):
        self.callbackfun = callbackfun
        self.regparam = regparam
        if isinstance(X, sp.base.spmatrix):
            self.X = X.todense()
        else:
            self.X = X
        self.X = self.X.T
        self.X = self.X.astype("float64", copy=False)
        self.Y = np.mat(array_tools.as_2d_array(Y))
        #Number of training examples
        self.size = self.Y.shape[0]
        self.bias = bias
        self.measure = None
        fsize = X.shape[1]
        self.desiredfcount = subsetsize
        if not fsize >= self.desiredfcount:
            raise Exception('The overall number of features ' + str(fsize) + ' is smaller than the desired number ' + str(self.desiredfcount) + ' of features to be selected.')
        self.results = {}
        if 'use_default_callback' in kwargs and bool(kwargs['use_default_callback']):
            self.callbackfun = DefaultCallback(**kwargs)
        #The current version works only with the squared error measure
        self._solve_cython(self.regparam)

Example 49

Project: RLScore Source File: rls.py
Function: hold_out
    def holdout(self, indices):
        """Computes hold-out predictions
        
        Parameters
        ----------
        indices : list of indices, shape = [n_hsamples]
            list of indices of training examples belonging to the set for which the
            hold-out predictions are calculated. The list can not be empty.

        Returns
        -------
        F : array, shape = [n_hsamples, n_labels]
            holdout predictions
            
        Notes
        -----
        
        Computational complexity of holdout:
        m = n_samples, d = n_features, l = n_labels, b = n_bvectors, h=n_hsamples
        
        O(h^3 + lmh): basic case
        
        O(min(h^3 + lh^2, d^3 + ld^2) +ldh): Linear Kernel, d < m
        
        O(min(h^3 + lh^2, b^3 + lb^2) +lbh): Sparse approximation with basis vectors 
        
        The fast holdout algorithm is based on results presented in [1,2]. However, the removal of basis vectors decribed in [2] is currently not implemented.
            
        References
        ----------
        
        [1] Tapio Pahikkala, Jorma Boberg, and Tapio Salakoski.
        Fast n-Fold Cross-Validation for Regularized Least-Squares.
        Proceedings of the Ninth Scandinavian Conference on Artificial Intelligence,
        83-90, Otamedia Oy, 2006.
        
        [2] Tapio Pahikkala, Hanna Suominen, and Jorma Boberg.
        Efficient cross-validation for kernelized least-squares regression with sparse basis expansions.
        Machine Learning, 87(3):381--407, June 2012.
        """
        indices = array_tools.as_index_list(indices, self.Y.shape[0])
        
        if len(indices) != len(np.unique(indices)):
            raise IndexError('Hold-out can have each index only once.')
        
        bevals = multiply(self.evals, self.newevals)
        A = self.svecs[indices]
        right = self.svecsTY - A.T * self.Y[indices] #O(hrl)
        RQY = A * multiply(bevals.T, right) #O(hrl)
        B = multiply(bevals.T, A.T)
        if len(indices) <= A.shape[1]: #h < r
            I = mat(identity(len(indices)))
            result = la.inv(I - A * B) * RQY #O(h^3 + h^2 * l)
        else: #h > r
            I = mat(identity(A.shape[1]))
            result = RQY - A * (la.inv(B * A - I) * (B * RQY)) #O(r^3 + r^2 * l + h * r * l)
        return np.squeeze(np.array(result))

Example 50

Project: python-control Source File: statefbk.py
def obsv(A, C):
    """Observability matrix

    Parameters
    ----------
    A, C: array_like or string
        Dynamics and output matrix of the system

    Returns
    -------
    O: matrix
        Observability matrix

    Examples
    --------
    >>> O = obsv(A, C)

   """

    # Convert input parameters to matrices (if they aren't already)
    amat = np.mat(A)
    cmat = np.mat(C)
    n = np.shape(amat)[0]

    # Construct the controllability matrix
    obsv = cmat
    for i in range(1, n):
        obsv = np.vstack((obsv, cmat*amat**i))
    return obsv
See More Examples - Go to Next Page
Page 1 Selected Page 2