Label Prop

import numpy as np
from scipy.sparse import csr_matrix
import scipy.sparse
from Kaczmarz_heap import runKaczmarz_heap
import Kaczmarz as KZ
import sample_descent as sd
import plot_experiment as plt
import legal_selection as LE
import heap_selection as HE
from pretty_plot_ import pretty_plot
from scipy.io import loadmat
#import FastKaczmark as FKZ
import scipy.sparse as sp
import scipy.io as sio
from numba import autojit
from numba import float64, int64



def save_sparse_csr(filename,array):
    np.savez(filename,data = array.data ,indices=array.indices,
             indptr =array.indptr, shape=array.shape )

def load_sparse_csr(filename):
    loader = np.load(filename)
    return csr_matrix((  loader['data'], loader['indices'], loader['indptr']),
                         shape = loader['shape'])

def get_leastsq_formulation(X,y):
    n, m = X.shape
    I = np.eye(n)
    Z = np.zeros((m,m))
    A = np.vstack([np.hstack([X,-I]),np.hstack([Z,X.T])])
    z = np.zeros((1,3)).flatten()
    b = np.hstack([y,z])
    return A, b

@autojit("uint8[:,:], float64[:,:],int64[:],int64,int64", nopython=True)
def get_label_prop_A(W,A,S,s,n):
    for i in range(s):
        for j in range(s):
            A[i,j] = -(W[S[j],S[i]]+W[S[i],S[j]])
        for j in range(n):
            A[i,i] += W[j,S[i]] + W[S[i],j]
    return


@autojit("uint8[:,:], uint8[:,:],float64[:],int64[:],int64,int64[:],int64,int64", nopython=True)
def get_label_prop_B(W,y,b,S,s,S_k,s_k,n):
    for i in range(s):
        for j in range(s_k):
            b[i] += (W[S[i],S_k[j]] + W[S_k[j],S[i]])*y[S_k[j]]
    return

def main():
    M = sio.loadmat("/home/alimv/Downloads/Sara/labelprop.mat")
    y = M["y"]
    W = M["W"]
    rows, column = W.shape
    s = np.sum(y==0)
    S = np.nonzero(y == 0)[0]
    assert S.shape[0] == s
    A = np.zeros((s,s))
    get_label_prop_A(W,A,S,s,rows)
    s_k = np.sum(y!=0)
    S_k = np.nonzero(y != 0)[0]
    assert S_k.shape[0] == s_k
    b = np.zeros(s)
    get_label_prop_B(W,y,b,S,s,S_k,s_k,rows)
    print b.shape
    print A.shape
    return

    n_samples = 5000
    n_features = 1000
    #A = csr_matrix(scipy.sparse.rand(m, n,(np.log(m)*np.log(n))/(n*m), random_state = np.random.RandomState(seed = 6666)))

    # O(log n) each column, stacking columns
    sparse_1 = sp.rand(n_samples, 1, density=np.log(n_samples)/(n_features*n_samples))
    for i in range(n_features-1):
        sparse_1 = sp.hstack([sparse_1, sp.rand(n_samples, 1, density=np.log(n_samples)/(2*n_samples))])

    # O(log m) each row, stacking rows
    sparse_2 = sp.rand(1, n_features, density=np.log(n_features)/(n_features*n_samples))
    for j in range(n_samples-1):
        sparse_2 = sp.vstack([sparse_2, sp.rand(1,n_features, density=np.log(n_features)/(2*n_features))])

    A  = csr_matrix(sparse_1 + sparse_2)

    np.random.seed(6666)
    x_true = np.random.randn(n_features)
    b = A.dot(x_true.T)

    #import scipy.io as sio
    #M = sio.loadmat("/home/alimv/Downloads/Sara/sara.mat")
    #A = M["A"].tocsr()






    #for algorithm in ["uniform","non_uniform","greedy_res","greedy_dis","legal_rule", "legal_rule_uni"]:
    # ["cycle", "uniform", "heap_rule_distance", "heap_rule_residual", "non_uniform"]:
    results = []

    for algorithm in ["non_uniform"]:
        print algorithm
        model = sd.Kaczmarz(n_epochs=10000, algorithm=algorithm)
        model.fit(A, b)

        results.append((model.results, algorithm))


    #for (y_axis, x_axis) in [("Distance", "Iteration"), ("Distance", "Time"),
    #                        ("Squared Error", "Iteration"), ("Squared Error", "Time")]:
    for (y_axis, x_axis) in [("Squared Error", "Time")]:
        pp = pretty_plot(title="samples="+str(n_samples)+", features="+str(n_features),
                         xlabel=x_axis, ylabel=y_axis,
                         n_plots=len(results))
        plt.plotkaczmarz(results, x_axis, y_axis, pp)
    #print "getting true x"
    #x_true = KZ.get_true_x(A,b)
    #print "have true x"
    """
    results = []
    #FKZ.runKaczmarz(A,b,x_true,"SADSRFdghjkgfdsadfgh")

    results.append((KZ.runKaczmarz(A,b,x_true, HE.getNextExample),"Heap"))

    #results.append((KZ.runKaczmarz(A,b,x_true, KZ.getNextExample_greedy),"Naive Greedy"))
    results.append((runKaczmarz_heap(A,b,x_true),"Heap"))
    #results.append((KZ.runKaczmarz(A,b,x_true, KZ.getNextExample_uni),"Uni"))
    #results.append((KZ.runKaczmarz(A,b,x_true, KZ.getNextExample_non_uni),"Non_Uni"))
    #results.append((KZ.runKaczmarz(A,b,x_true, LE.getNextExample),"Adapt"))
    #results.append((KZ.runKaczmarz(A,b,x_true, LE.getNextExample_aux),"Adapt-Uni"))
    for (y_axis, x_axis) in [("Distance", "Iteration"), ("Distance", "Time"),
                             ("Squared Error", "Iteration"), ("Squared Error", "Time")]:
        pp = pretty_plot(title="(less sparse) samples="+str(m)+", features="+str(n), xlabel=x_axis, ylabel=y_axis, n_plots=len(results))
        plt.plotkaczmarz(results, x_axis, y_axis, pp)
    """
if __name__ == "__main__":
    main()