import numpy as np

def rbf(X, gamma):
    assert(gamma > 0)
    n, d = np.shape(X)
    dotproduct = np.dot( X, X.T)
    diag_dotproduct = np.diag(dotproduct)
    allones = np.matrix(np.ones(n))
    offdiagonal = np.outer(diag_dotproduct, allones.T)
    distance_matrix = offdiagonal + offdiagonal.T - 2* dotproduct
    return np.exp( -1.0 * gamma * distance_matrix)

# p is dimension
# n is number of data points
def get_synthetic_data(n,p, gamma=0.2,holdfrac=0.2):
    nnz=p

    mu = np.zeros(p)
    D = np.matrix(np.reshape(np.random.normal(0, 1, p * p), (p, p)))
    Sigma = rbf(D, gamma)

    X = np.array(np.random.multivariate_normal(mu, Sigma, n))  # design matrix

    # [ss, uu] = np.linalg.eig(Sigma)

    b = np.zeros(p)
    nnzInd = np.random.permutation(p)
    b[nnzInd[:nnz]] = 1
    s = X.dot(b)
    yp = 1 / (1 + np.exp(-s))

    # y = (np.random.ranf(N) <= yp).astype(int)

    y = (yp >= 0.5).astype(int)

    numHold = np.floor(holdfrac * n).astype(int)
    y_test = y[:numHold]
    y_train = y[numHold:]
    X_test = np.hstack((X[:numHold, :], np.ones((numHold, 1))))
    X_train = np.hstack((X[numHold:, :], np.ones((n - numHold, 1))))

    return X_train, X_test, y_train, y_test