res_reg_lmnt_awikner.csc_mult

  1import numpy as np
  2from scipy.sparse import csc_matrix, diags, csr_matrix, coo_matrix
  3from numba import jit, objmode
  4
  5#def create_csc_matrix(data, indices, indptr, shape):
  6#    return csc_matrix((data, indices, indptr), shape = (shape[0], shape[1]))
  7
  8@jit(nopython = True, fastmath = True)
  9def matrix_dot_left_T(data, indices, indptr, shape, mat):
 10    with objmode(out = 'double[:,:]'):
 11        csc = csc_matrix((data, indices, indptr), shape = (shape[0], shape[1]))
 12        out = csc.T.dot(mat).T
 13    return out
 14
 15@jit(nopython = True, fastmath = True)
 16def matrix_diag_mult(dmat, b):
 17    with objmode(out = 'double[:,:]'):
 18        out = diags(dmat, format = 'csc').dot(b)
 19    return out
 20
 21@jit(nopython = True, fastmath = True)
 22def matrix_sparse_mult(sp, mat):
 23    with objmode(out = 'double[:,:]'):
 24        out = csc_matrix(sp).dot(mat)
 25    return out
 26
 27@jit(nopython = True, fastmath = True)
 28def matrix_diag_sparse_mult_add(diag_mat, data, indices, indptr, shape, add_mat):
 29    with objmode(out = 'double[:,:]'):
 30        out_sp = diags(diag_mat, format = 'csc').dot(csc_matrix((data, indices, indptr), shape = (shape[0], shape[1]))) + csc_matrix(add_mat)
 31        out    = out_sp.toarray()
 32    return out
 33
 34@jit(nopython = True, fastmath = True)
 35def matrix_diag_sparse_mult_add_sparse(diag_mat, data, indices, indptr, shape, add_mat_data, add_mat_indices,\
 36    add_mat_indptr, add_mat_shape):
 37    with objmode(data = 'double[:]', indices = 'int32[:]', indptr = 'int32[:]', shape = 'int32[:]'):
 38        out = diags(diag_mat, format = 'csc').dot(csc_matrix((data, indices, indptr), shape = (shape[0], shape[1]))) + csc_matrix((add_mat_data, add_mat_indices, add_mat_indptr), shape = (shape[0], shape[1]))
 39        data = out.data
 40        indices = out.indices
 41        indptr = out.indptr
 42        shape = np.array(list(out.shape), dtype = np.int32)
 43
 44    return data, indices, indptr, shape
 45
 46@jit(nopython = True, fastmath = True)
 47def matrix_diag_sparse_diag_mult(dmat, b_data, b_indices, b_indptr, b_shape, dmat_inv):
 48    with objmode(data = 'double[:]', indices = 'int32[:]', indptr = 'int32[:]', shape = 'int32[:]'):
 49        out = diags(dmat, format = 'csc').dot(csc_matrix((b_data, b_indices, b_indptr), shape = (b_shape[0], b_shape[1]))).dot(diags(dmag_inv, format = 'csc'))
 50        data = out.data
 51        indices = out.indices
 52        indptr = out.indptr
 53        shape = np.array(list(out.shape), dtype = np.int32)
 54
 55    return data, indices, indptr, shape
 56
 57@jit(nopython = True, fastmath = True)
 58def matrix_diag_sparse_mult(dmat, b_data, b_indices, b_indptr, b_shape):
 59    with objmode(data = 'double[:]', indices = 'int32[:]', indptr = 'int32[:]', shape = 'int32[:]'):
 60        out = diags(dmat, format = 'csc').dot(csc_matrix((b_data, b_indices, b_indptr), shape = (b_shape[0], b_shape[1])))
 61        data = out.data
 62        indices = out.indices
 63        indptr = out.indptr
 64        shape = np.array(list(out.shape), dtype = np.int32)
 65
 66    return data, indices, indptr, shape
 67
 68@jit(nopython = True, fastmath = True)
 69def dense_to_sparse(W):
 70    with objmode(data = 'double[:]', indices = 'int32[:]', indptr = 'int32[:]', shape = 'int32[:]'):
 71        out = csc_matrix(W)
 72        data = out.data
 73        indices = out.indices
 74        indptr = out.indptr
 75        shape = np.array(list(out.shape), dtype = np.int32)
 76
 77    return data, indices, indptr, shape
 78
 79
 80@jit(nopython = True, fastmath = True)
 81def mult_vec(data, indices, indptr, shape, mat):
 82    assert(shape[1] == mat.size)
 83    out = np.zeros(shape[0])
 84    for i in range(mat.size):
 85        for k in range(indptr[i], indptr[i+1]):
 86            out[indices[k]] += data[k] * mat[i]
 87    return out
 88
 89@jit(nopython = True, fastmath = True)
 90def construct_jac_mat_csc(Win, data_in, indices_in, indptr_in, shape_in, rsvr_size, n, squarenodes = False):
 91    with objmode(data = 'double[:]', indices = 'int32[:]', indptr = 'int32[:]', shape = 'int32[:]'):
 92        W_conv = csc_matrix((data_in, indices_in, indptr_in), shape = (shape_in[0], shape_in[1])).toarray()
 93        if squarenodes:
 94            mat    = csc_matrix(np.concatenate((Win[:,0].reshape(-1,1), W_conv, np.zeros((rsvr_size,n+rsvr_size))), axis = 1))
 95        else:
 96            mat    = csc_matrix(np.concatenate((Win[:,0].reshape(-1,1), W_conv, np.zeros((rsvr_size,n))), axis = 1))
 97        data   = mat.data
 98        indices = mat.indices
 99        indptr = mat.indptr
100        shape = np.array(list(mat.shape), dtype = np.int32)
101    return data, indices, indptr, shape
102
103@jit(nopython = True, fastmath = True)
104def construct_jac_mat_csc_csc(Win_data, Win_indices, Win_indptr, Win_shape,\
105        data_in, indices_in, indptr_in, shape_in, rsvr_size,  n, squarenodes = False):
106    with objmode(data = 'double[:]', indices = 'int32[:]', indptr = 'int32[:]', shape = 'int32[:]'):
107        W_conv = csc_matrix((data_in, indices_in, indptr_in), shape = (shape_in[0],  shape_in[1])).toarray()
108        Win_conv = csc_matrix((Win_data, Win_indices, Win_indptr), shape = (Win_shape[0],  Win_shape[1])).toarray()
109        if squarenodes:
110            mat    = csc_matrix(np.concatenate((Win_conv[:,0].reshape(-1,1), W_conv, np.zeros((rsvr_size,n+rsvr_size))), axis = 1))
111        else:
112            mat    = csc_matrix(np.concatenate((Win_conv[:,0].reshape(-1,1), W_conv, np.zeros((rsvr_size,n))), axis = 1))
113        data   = mat.data
114        indices = mat.indices
115        indptr = mat.indptr
116        shape = np.array(list(mat.shape), dtype = np.int32)
117    return data, indices, indptr, shape
118
119@jit(nopython = True, fastmath = True)
120def construct_leakage_mat(rsvr_size, n, leakage, squarenodes):
121    if squarenodes:
122        leakage_mat = np.concatenate((np.zeros((rsvr_size,1)), (1-leakage) *
123                            np.identity(rsvr_size), np.zeros((rsvr_size, n+rsvr_size))), axis=1)
124    else:
125        leakage_mat = np.concatenate((np.zeros((rsvr_size,1)), (1-leakage) *
126                             np.identity(rsvr_size), np.zeros((rsvr_size, n))),   axis=1)
127    with objmode(data = 'double[:]', indices = 'int32[:]', indptr = 'int32[:]', shape = 'int32[:]'):
128        mat = csc_matrix(leakage_mat)
129        data   = mat.data
130        indices = mat.indices
131        indptr = mat.indptr
132        shape = np.array(list(mat.shape), dtype = np.int32)
133
134    return data, indices, indptr, shape
135
136@jit(nopython = True, fastmath = True)
137def construct_leakage_mat_mlonly(rsvr_size, n, leakage, squarenodes):
138    leakage_mat = (1-leakage) * np.identity(rsvr_size)
139    with objmode(data = 'double[:]', indices = 'int32[:]', indptr = 'int32[:]', shape = 'int32[:]'):
140        mat = csc_matrix(leakage_mat)
141        data   = mat.data
142        indices = mat.indices
143        indptr = mat.indptr
144        shape = np.array(list(mat.shape), dtype = np.int32)
145
146    return data, indices, indptr, shape
147
148@jit(nopython = True, fastmath = True)
149def get_Win_nobias(data_in, indices_in, indptr_in, shape_in):
150    with objmode(data = 'double[:]', indices = 'int32[:]', indptr = 'int32[:]', shape = 'int32[:]'):
151        Win = csc_matrix((data_in, indices_in, indptr_in), shape = (shape_in[0],  shape_in[1])).toarray()
152        Win_nobias = csc_matrix(np.ascontiguousarray(Win[:, 1:]))
153        data = Win_nobias.data
154        indices = Win_nobias.indices
155        indptr = Win_nobias.indptr
156        shape = np.array(list(Win_nobias.shape), dtype = np.int32)
157    return data, indices, indptr, shape
158
159
160@jit(nopython = True, fastmath = True)
161def coo_to_csc(data_in, row_in, col_in, shape_in):
162    with objmode(data = 'double[:]', indices = 'int32[:]', indptr = 'int32[:]'):
163        mat     = coo_matrix((data_in, (row_in, col_in)), shape = (shape_in[0], shape_in[1])).tocsc()
164        data    = mat.data
165        indices = mat.indices
166        indptr  = mat.indptr
167    return data, indices, indptr
168
169@jit(nopython = True, fastmath = True)
170def csc_to_coo(data_in, indices_in, indptr_in, shape_in):
171    with objmode(data = 'double[:]', rows = 'int32[:]', cols = 'int32[:]'):
172        mat    = csc_matrix((data_in, indices_in, indptr_in), shape = (shape_in[0], shape_in[1])).tocoo()
173        data   = mat.data
174        rows   = mat.row
175        cols   = mat.col
176    return data, rows, cols
177
178@jit(nopython = True, fastmath = True)
179def matrix_diag_sparse_mult_sparse_add(diag_mat, data, indices, indptr, shape, \
180        add_mat_data, add_mat_indices, add_mat_indptr, add_mat_shape):
181    with objmode(data = 'double[:]', indices = 'int32[:]', indptr = 'int32[:]', shape = 'int32[:]'):
182        mat = diags(diag_mat, format = "csc").dot(csc_matrix((data, indices, indptr), shape = (shape[0], shape[1])))\
183            + csc_matrix((add_mat_data, add_mat_indices, add_mat_indptr), shape = (add_mat_shape[0], add_mat_shape[1]))
184        data   = mat.data
185        indices = mat.indices
186        indptr = mat.indptr
187        shape = np.array(list(mat.shape), dtype = np.int32)
188    return data, indices, indptr, shape
189
190@jit(nopython = True, fastmath = True)
191def matrix_sparse_sparseT_mult(mat1_data, mat1_indices, mat1_indptr, mat1_shape):
192    with objmode(out = 'double[:,:]'):
193        mat_sp = csc_matrix((mat1_data, mat1_indices, mat1_indptr), shape = (mat1_shape[0], mat1_shape[1]))
194        mat_sp2 = mat_sp @ mat_sp.T
195        out = mat_sp2.toarray()
196    return out
197
198@jit(nopython = True, fastmath = True)
199def matrix_sparse_sparse_mult(mat1_data, mat1_indices, mat1_indptr, mat1_shape,\
200                              mat2_data, mat2_indices, mat2_indptr, mat2_shape):
201    with objmode(data = 'double[:]', indices = 'int32[:]', indptr = 'int32[:]', shape = 'int32[:]'):
202        mat = csc_matrix((mat1_data, mat1_indices, mat1_indptr), shape = (mat1_shape[0], mat1_shape[1])) @ \
203              csc_matrix((mat2_data, mat2_indices, mat2_indptr), shape = (mat2_shape[0], mat2_shape[1]))
204        data     = mat.data
205        indices  = mat.indices
206        indptr   = mat.indptr
207        shape    = np.array(list(mat.shape), dtype = np.int32)
208    return data, indices, indptr, shape
209
210@jit(nopython = True, fastmath = True)
211def matrix_sparse_sparse_conv_mult(mat1_data, mat1_indices, mat1_indptr, mat1_shape,\
212                              mat2_data, mat2_indices, mat2_indptr, mat2_shape):
213    with objmode(data = 'double[:]', indices = 'int32[:]', indptr = 'int32[:]', shape = 'int32[:]'):
214        mat = csc_matrix(csc_matrix((mat1_data, mat1_indices, mat1_indptr), shape = (mat1_shape[0], mat1_shape[1])) @ \
215              csc_matrix((mat2_data, mat2_indices, mat2_indptr), shape = (mat2_shape[0], mat2_shape[1])).toarray())
216        data     = mat.data
217        indices  = mat.indices
218        indptr   = mat.indptr
219        shape    = np.array(list(mat.shape), dtype = np.int32)
220    return data, indices, indptr, shape
221
222@jit(nopython = True, fastmath = True)
223def matrix_sparse_sparseT_conv_mult(mat1_data, mat1_indices, mat1_indptr, mat1_shape):
224    with objmode(out = 'double[:,:]'):
225        mat = csc_matrix((mat1_data, mat1_indices, mat1_indptr), shape = (mat1_shape[0], mat1_shape[1])).toarray()
226        out = mat @ mat.T
227    return out
228
229@jit(fastmath = True, nopython = True)
230def get_D_n(D, X, Win_nobias_data, Win_nobias_indices, Win_nobias_indptr, Win_nobias_shape, \
231        D_n_shape, rsvr_size, res_feature_size, n, squarenodes):
232    D_n_base_data, D_n_base_indices, D_n_base_indptr, D_n_base_shape = matrix_diag_sparse_mult(D, \
233        Win_nobias_data, Win_nobias_indices, Win_nobias_indptr, Win_nobias_shape)
234    D_n_data_coo, D_n_rows, D_n_cols = csc_to_coo(\
235        D_n_base_data, D_n_base_indices, D_n_base_indptr, D_n_base_shape)
236    D_n_rows += 1
237    if squarenodes:
238        D_n_2_data, D_n_2_indices, D_n_2_indptr, D_n_2_shape = matrix_diag_sparse_mult(2*X, D_n_base_data,\
239            D_n_base_indices, D_n_base_indptr, D_n_base_shape)
240        #D_n[1+rsvr_size:1+res_feature_size,:,i] = matrix_diag_mult(2*X[:, i], D_n_base)
241        D_n_2_data_coo, D_n_2_rows, D_n_2_cols = csc_to_coo(\
242            D_n_2_data, D_n_2_indices, D_n_2_indptr, D_n_2_shape)
243        D_n_2_rows += 1+rsvr_size
244        D_n_data_coo, D_n_rows, D_n_cols = np.concatenate((D_n_data_coo, D_n_2_data_coo, np.ones(n))),\
245                np.concatenate((D_n_rows, D_n_2_rows, np.arange(res_feature_size+1,D_n_shape[0]))),\
246                np.concatenate((D_n_cols, D_n_2_cols, np.arange(D_n_shape[1])))
247    else:
248        #D_n[1:rsvr_size+1, :, i] = D_n_base
249        D_n_data_coo, D_n_rows, D_n_cols = np.append(D_n_data_coo, np.ones(n)),\
250                np.append(D_n_rows, np.arange(res_feature_size+1,D_n_shape[0])),\
251                np.append(D_n_cols, np.arange(D_n_shape[1]))
252    D_n_data, D_n_indices, D_n_indptr = coo_to_csc(D_n_data_coo, D_n_rows, D_n_cols, D_n_shape)
253        #D_n[res_feature_size+1:, :, i] = np.identity(n)
254    return D_n_data, D_n_indices, D_n_indptr
255
256@jit(fastmath = True, nopython = True)
257def get_D_n_mlonly(D, Win_nobias_data, Win_nobias_indices, Win_nobias_indptr, Win_nobias_shape, \
258        D_n_shape, rsvr_size, res_feature_size, n, squarenodes):
259    D_n_data, D_n_indices, D_n_indptr, tmp = matrix_diag_sparse_mult(D, Win_nobias_data, Win_nobias_indices,\
260        Win_nobias_indptr, Win_nobias_shape)
261    return D_n_data,D_n_indices,D_n_indptr
262
263@jit(fastmath = True, nopython = True)
264def get_E_n(D, X, E_n_shape, rsvr_size, W_mat_data, W_mat_indices, W_mat_indptr, W_mat_shape, \
265            leakage_mat_data, leakage_mat_indices, leakage_mat_indptr, leakage_mat_shape, squarenodes):
266    E_n_base_data, E_n_base_indices, E_n_base_indptr, E_n_base_shape = matrix_diag_sparse_mult_sparse_add(
267        D, W_mat_data, W_mat_indices, W_mat_indptr, W_mat_shape, leakage_mat_data, \
268        leakage_mat_indices, leakage_mat_indptr, leakage_mat_shape)
269    E_n_base_coo, E_n_base_rows, E_n_base_cols = csc_to_coo(E_n_base_data,\
270        E_n_base_indices, E_n_base_indptr, E_n_base_shape)
271    E_n_base_rows += 1
272    if squarenodes:
273        E_n_2_coo, E_n_2_rows, E_n_2_cols = csc_to_coo(*matrix_diag_sparse_mult(2*X,\
274            E_n_base_data, E_n_base_indices, E_n_base_indptr, E_n_base_shape))
275        E_n_2_rows += 1+rsvr_size
276        #E_n[1+rsvr_size:1+res_feature_size,:,i-1] = matrix_diag_mult(2*X[:,  i], E_n_base)
277        E_n_data_coo, E_n_rows, E_n_cols = np.append(E_n_base_coo, E_n_2_coo),\
278            np.append(E_n_base_rows, E_n_2_rows), np.append(E_n_base_cols, E_n_2_cols)
279    else:
280        E_n_data_coo, E_n_rows, E_n_cols = E_n_base_coo, E_n_base_rows, E_n_base_cols
281    E_n_data, E_n_indices, E_n_indptr = coo_to_csc(E_n_data_coo, E_n_rows, E_n_cols, E_n_shape)
282    return E_n_data, E_n_indices, E_n_indptr
283
284@jit(fastmath = True, nopython = True)
285def get_E_n_mlonly(D, E_n_shape, rsvr_size, W_mat_data, W_mat_indices, W_mat_indptr, W_mat_shape, \
286            leakage_mat_data, leakage_mat_indices, leakage_mat_indptr, leakage_mat_shape, squarenodes):
287    E_n_data, E_n_indices, E_n_indptr, tmp = matrix_diag_sparse_mult_add_sparse(D, W_mat_data, W_mat_indices, W_mat_indptr, W_mat_shape, leakage_mat_data, leakage_mat_indices, leakage_mat_indptr, leakage_mat_shape)
288    return E_n_data, E_n_indices, E_n_indptr