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
3
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)
3
Example 2
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))
3
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
3
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)
3
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)
3
Example 6
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)
3
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)
3
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
3
Example 9
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)))
3
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
3
Example 11
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
3
Example 12
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)
3
Example 13
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
3
Example 14
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'))
3
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)
3
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)
3
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
3
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
3
Example 19
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)
3
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
3
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
3
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'])
3
Example 23
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
3
Example 24
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]))
0
Example 25
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]]))
0
Example 26
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)
0
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)
0
Example 28
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)
0
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
0
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)
0
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
0
Example 32
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()
0
Example 33
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
0
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
0
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
0
Example 36
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
0
Example 37
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)
0
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
0
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))
0
Example 40
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)
0
Example 41
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])
0
Example 42
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
0
Example 43
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
0
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
0
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
0
Example 46
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
0
Example 47
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)
0
Example 48
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)
0
Example 49
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))
0
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