add RSLDA code

HoganHPH · HoganHPH · commit 7196b1709ca2 · 2023-04-07T18:04:04.000+07:00
diff --git a/PCA1.py b/PCA1.py
@@ -0,0 +1,77 @@
+import os
+import pickle
+from mySVD import mySVD
+from scipy.sparse import issparse
+
+import numpy as np
+
+
+def PCA1(data, options=None):
+    # Principal Component Analysis
+
+    # Input:
+    # data - Data matrix. Each row vector of data is a data point.
+    # options.ReducedDim - The dimensionality of the reduced subspace. If 0, all the dimensions will be kept. Default is 0.
+    # Output:
+    # eigvector - Each column is an embedding function, for a new data point (row vector) x, y = x*eigvector will be the embedding result of x.
+    # eigvalue - The sorted eigvalue of PCA eigen-problem.
+
+    if options is None:
+        options = {}
+
+    ReducedDim = 0
+    if 'ReducedDim' in options:
+        ReducedDim = options['ReducedDim']
+
+    nSmp, nFea = data.shape
+    if (ReducedDim > nFea) or (ReducedDim <=0):
+        ReducedDim = nFea
+
+    if issparse(data):
+        data = data.toarray()
+    sampleMean = np.mean(data, axis=0)
+    data = (data - np.tile(sampleMean, (nSmp, 1)))
+
+    eigvector, eigvalue, _ = mySVD(data.T, ReducedDim)
+    eigvalue = np.square(eigvalue)
+
+    if 'PCARatio' in options:
+        sumEig = np.sum(eigvalue)
+        sumEig *= options['PCARatio']
+        sumNow = 0
+        for idx in range(len(eigvalue)):
+            sumNow += eigvalue[idx]
+            if sumNow >= sumEig:
+                break
+        eigvector = eigvector[:, :idx]
+
+    return eigvector, eigvalue
+
+
+if __name__ == '__main__':
+    dim = 100
+    options = {
+        "ReducedDim": dim
+    }
+
+    DATA_FOLDER = "P:/RESEARCH/DATA/CIFAR-100/CIFARDB/train/"
+    features = pickle.load(open("./database/features_YCrCb_CifarDataset.pkl", 'rb'))
+    paths = pickle.load(open("./database/paths_YCrCb_CifarDataset.pkl", 'rb'))
+    features = features
+    paths = paths
+    classes = os.listdir(DATA_FOLDER)
+    X = np.array(features)
+    y = []
+    class_check = 'apple'
+    for path in paths:
+        label = 0
+        if class_check in path:
+            label = 1
+        y.append(label)
+    y = np.array(y)
+    y = y.reshape(-1, 1)
+
+    eigvector, eigvalue = PCA1(X, options)
+    print("Shape of eigvector = ", eigvector.shape)
+    print("Shape of eigvalue = ", eigvalue.shape)
+
diff --git a/RSLDA.py b/RSLDA.py
@@ -0,0 +1,145 @@
+from PCA1 import PCA1
+from ScatterMat import ScatterMat
+
+import os
+import pickle
+import numpy as np
+
+from tqdm.auto import tqdm
+
+def RSLDA(X=None, label=None, lambda1 = 0.0002,
+                              lambda2 = 0.001,
+                              dim = 100,
+                              mu  = 0.1,
+                              rho = 1.01,
+                              max_iter = 100):
+    print("RUNNING RSLDA!")
+    m, n = X.shape
+    max_mu = 10**5
+
+    # Initialization
+    regu = 10**-5
+    Sw, Sb = ScatterMat(X, label)
+    options = {}
+    options['ReducedDim'] = dim
+    P1, _ = PCA1(X.T, options)
+    Q = np.ones((m, dim))
+    E = np.zeros((m, n))
+    Y = np.zeros((m, n))
+    v = np.sqrt(np.sum(Q*Q, axis=1) + np.finfo(float).eps)
+    D = np.diag(1./v)
+
+    # Main loop
+    for iter in tqdm(range(1, max_iter+1), total=max_iter):
+        
+        # Update P
+        if iter == 1:
+            P = P1
+        else:
+            M = X - E + Y/mu
+            U1, S1, V1 = np.linalg.svd(M @ X.T @ Q, full_matrices=False)
+            P = U1 @ V1
+            del M
+        
+        # Update Q
+        M = X - E + Y/mu
+        Q1 = 2*(Sw - regu*Sb) + lambda1*D + mu*X @ X.T
+        Q2 = mu*X @ M.T @ P
+        Q = np.linalg.solve(Q1, Q2)
+        v = np.sqrt(np.sum(Q*Q, axis=1) + np.finfo(float).eps)
+        D = np.diag(1./v)
+        
+        # Update E
+        eps1 = lambda2/mu
+        temp_E = X - P @ Q.T @ X + Y/mu
+        E = np.maximum(0, temp_E - eps1) + np.minimum(0, temp_E + eps1)
+        
+        # Update Y, mu
+        Y = Y + mu*(X - P @ Q.T @ X - E)
+        mu = min(rho*mu, max_mu)
+        leq = X - P @ Q.T @ X - E
+        EE = np.sum(np.abs(E), axis=1)
+        obj = np.trace(Q.T @ (Sw - regu*Sb) @ Q) + lambda1*np.sum(v) + lambda2*np.sum(EE)
+        
+        if iter > 2:
+            if np.linalg.norm(leq, np.inf) < 10**-7 and abs(obj - obj_prev) < 0.00001:
+                print(iter)
+                break
+        
+        obj_prev = obj
+    
+    return P, Q, E, obj
+
+
+def sort_power_of_features(Q):
+    row_norm = np.linalg.norm(Q, axis=1)
+    sorted_power = np.argsort(row_norm)[::-1] # DESC
+    return sorted_power
+
+
+if __name__ == '__main__':
+    # DATA_FOLDER = "P:/RESEARCH/DATA/CIFAR-100/CIFARDB/train/"
+
+    features = pickle.load(open("./database/features_handcrafted_Cifar.pkl", 'rb'))
+    paths = pickle.load(open("./database/paths_handcrafted_Cifar.pkl", 'rb'))
+
+    # X = np.array(features)
+    y = []
+
+    # class_check = 'apple'
+    # for i in range(len(paths)):
+    #     path = paths[i]
+    #     label =  0
+    #     if class_check in path:
+    #         label = 1
+    #     y.append(label)
+
+    X_train = []
+    n_pos = 0
+    n_neg = 0
+    
+    class_check = 'apple'
+    for i in range(len(paths)):
+        if n_pos >= 10 and n_neg >=100:
+            break
+        else:
+            path = paths[i]
+            if class_check in path and n_pos < 10:
+                label = 1
+                n_pos += 1
+                y.append(label)
+                X_train.append(features[i])
+            elif class_check not in path and n_neg < 50:
+                label = 0
+                n_neg += 1
+                y.append(label)
+                X_train.append(features[i])
+    
+    X_train = np.array(X_train)
+    X_train = X_train.T
+    # X = X.T
+    y = np.array(y)
+    y = y.reshape(-1, 1)
+
+    print(X_train.shape)
+    print(y.shape)
+    P,Q,E,obj = RSLDA(X=X_train,label=y)
+    print("Shape of P = ", P.shape)
+    print("Shape of Q = ", Q.shape)
+    print("Shape of E = ", E.shape)
+
+    # pickle.dump(Q, open("./database/Q.pkl", 'wb'))
+
+    # TEST
+    # Q = pickle.load(open("./database/Q.pkl", 'rb'))
+    # print(Q.shape)
+    # print(np.max(Q))
+    # print(np.min(Q))
+    # a = np.count_nonzero(Q)
+    # print(a)
+    # sort_power_of_features(Q=Q)
+    # print(Q.shape)
+
+
+
+
diff --git a/ScatterMat.py b/ScatterMat.py
@@ -0,0 +1,47 @@
+import os
+import numpy as np
+import pickle
+
+
+def ScatterMat(X, y):
+    dim, _ = X.shape
+    nclass = int(np.max(y)) + 1
+    
+    mean_X = np.mean(X, axis=1)
+    Sw = np.zeros((dim,dim))
+    Sb = np.zeros((dim,dim))
+    
+    for i in range(nclass):
+        inx_i = np.where(y==i)[0]
+        X_i = X[:,inx_i]
+        
+        mean_Xi = np.mean(X_i, axis=1)
+        Sw += np.cov(X_i, rowvar=True, bias=True)
+        Sb += len(inx_i)*(mean_Xi-mean_X)[:,None]*(mean_Xi-mean_X)[None,:]
+        
+    return Sw, Sb
+
+
+if __name__ == '__main__':
+    DATA_FOLDER = "P:/RESEARCH/DATA/CIFAR-100/CIFARDB/train/"
+    features = pickle.load(open("./database/features_YCrCb_CifarDataset.pkl", 'rb'))
+    paths = pickle.load(open("./database/paths_YCrCb_CifarDataset.pkl", 'rb'))
+    features = features
+    paths = paths
+    classes = os.listdir(DATA_FOLDER)
+    X = np.array(features)
+    y = []
+    class_check = 'apple'
+    for path in paths:
+        label = 0
+        if class_check in path:
+            label = 1
+        y.append(label)
+    y = np.array(y)
+    X = X.T
+    y = y.reshape(-1, 1)
+    print(X.shape)
+    print(y.shape)
+    Sw, Sb = ScatterMat(X, y)
+    print(Sw.shape, Sb.shape)
+    print(np.max(Sw), np.max(Sb))
diff --git a/mySVD.py b/mySVD.py
@@ -0,0 +1,110 @@
+import os
+import pickle
+import numpy as np
+import scipy as sp
+from scipy.sparse import issparse, spdiags
+from scipy.linalg import eig, eigh, svd
+
+
+def mySVD(X, ReducedDim=0):
+    MAX_MATRIX_SIZE = 1600
+    EIGVECTOR_RATIO = 0.1
+    
+    if not ReducedDim:
+        ReducedDim = 0
+    
+    mFea, nSmp = X.shape
+    
+    if nSmp/mFea > 1.0713:
+        ddata = X @ X.T
+        ddata = np.maximum(ddata, ddata.T)
+        
+        dimMatrix = ddata.shape[0]
+        if (ReducedDim > 0) and (dimMatrix > MAX_MATRIX_SIZE) and (ReducedDim < dimMatrix*EIGVECTOR_RATIO):
+            eigvalue, U = eigh(ddata, eigvals=(dimMatrix - ReducedDim, dimMatrix - 1))
+            eigvalue = eigvalue[::-1]
+            U = U[:, ::-1]
+        else:
+            if issparse(ddata):
+                ddata = ddata.toarray()
+                
+            eigvalue, U = eigh(ddata)
+            eigvalue = eigvalue[::-1]
+            U = U[:, ::-1]
+        
+        eigIdx = np.where(np.abs(eigvalue) / np.max(np.abs(eigvalue)) < 1e-10)[0]
+        eigvalue = np.delete(eigvalue, eigIdx)
+        U = np.delete(U, eigIdx, axis=1)
+        
+        if (ReducedDim > 0) and (ReducedDim < len(eigvalue)):
+            eigvalue = eigvalue[:ReducedDim]
+            U = U[:, :ReducedDim]
+        
+        eigvalue_Half = np.sqrt(np.abs(eigvalue))
+        S = spdiags(eigvalue_Half, 0, len(eigvalue_Half), len(eigvalue_Half))
+        
+        eigvalue_MinusHalf = eigvalue_Half ** -1
+        V = X.T @ (U * np.tile(eigvalue_MinusHalf, (U.shape[0], 1)))
+    else:
+        ddata = X.T @ X
+        ddata = np.maximum(ddata, ddata.T)
+        
+        dimMatrix = ddata.shape[0]
+        if (ReducedDim > 0) and (dimMatrix > MAX_MATRIX_SIZE) and (ReducedDim < dimMatrix*EIGVECTOR_RATIO):
+            eigvalue, V = eigh(ddata, eigvals=(dimMatrix - ReducedDim, dimMatrix - 1))
+            eigvalue = eigvalue[::-1]
+            V = V[:, ::-1]
+        else:
+            if issparse(ddata):
+                ddata = ddata.toarray()
+                
+            eigvalue, V = eigh(ddata)
+            eigvalue = eigvalue[::-1]
+            V = V[:, ::-1]
+        
+        eigIdx = np.where(np.abs(eigvalue) / np.max(np.abs(eigvalue)) < 1e-10)[0]
+        eigvalue = np.delete(eigvalue, eigIdx)
+        V = np.delete(V, eigIdx, axis=1)
+        
+        if (ReducedDim > 0) and (ReducedDim < len(eigvalue)):
+            eigvalue = eigvalue[:ReducedDim]
+            V = V[:, :ReducedDim]
+        
+        eigvalue_Half = eigvalue**0.5
+        S =  spdiags(eigvalue_Half,0,len(eigvalue_Half),len(eigvalue_Half))
+
+        eigvalue_MinusHalf = eigvalue_Half**-1
+
+        U = X.dot(V * np.tile(eigvalue_MinusHalf.T,(V.shape[0],1)))
+
+    return (U, S, V)
+
+
+
+if __name__ == '__main__':
+    DATA_FOLDER = "P:/RESEARCH/DATA/CIFAR-100/CIFARDB/train/"
+    features = pickle.load(open("./database/features_YCrCb_CifarDataset.pkl", 'rb'))
+    paths = pickle.load(open("./database/paths_YCrCb_CifarDataset.pkl", 'rb'))
+    features = features
+    paths = paths
+    classes = os.listdir(DATA_FOLDER)
+    X = np.array(features)
+    y = []
+    class_check = 'apple'
+    for path in paths:
+        label = 0
+        if class_check in path:
+            label = 1
+        y.append(label)
+    y = np.array(y)
+    X = X.T
+    y = y.reshape(-1, 1)
+    print(X.shape)
+    print(y.shape)
+    
+
+    # ====================================
+    U, S, V = mySVD(X, ReducedDim=100)
+    print("Size of U = ", U.shape)
+    print("Size of S = ", S.shape)
+    print("Size of V = ", V.shape)
