def __init__(self, X_train, y_train, X_test, y_text, batch_size=100, max_iter=1000, M=20, n_hidden=50, a0=1, b0=0.1, master_stepsize=0.001, auto_corr=0.9, h=(- 1), alpha=0.9, method='none', m=5, cf=False, uStat=True, regCoeff=0.1, adver=False, adverMaxIter=5, maxTime=20, numTimeSteps=20):
    self.n_hidden = n_hidden
    self.d = X_train.shape[1]
    self.M = M
    num_vars = (((self.d * n_hidden) + (n_hidden * 2)) + 3)
    self.theta = np.zeros([self.M, num_vars])
    '\n            We keep the last 10% (maximum 500) of training data points for model developing\n        '
    size_dev = min(int(np.round((0.1 * X_train.shape[0]))), 500)
    (X_dev, y_dev) = (X_train[(- size_dev):], y_train[(- size_dev):])
    (X_train, y_train) = (X_train[:(- size_dev)], y_train[:(- size_dev)])
    '\n            The data sets are normalized so that the input features and the targets have zero mean and unit variance\n        '
    self.std_X_train = np.std(X_train, 0)
    self.std_X_train[(self.std_X_train == 0)] = 1
    self.mean_X_train = np.mean(X_train, 0)
    self.mean_y_train = np.mean(y_train)
    self.std_y_train = np.std(y_train)
    '\n            Theano symbolic variables\n            Define the neural network here\n        '
    X = T.matrix('X')
    y = T.vector('y')
    w_1 = T.matrix('w_1')
    b_1 = T.vector('b_1')
    w_2 = T.vector('w_2')
    b_2 = T.scalar('b_2')
    N = T.scalar('N')
    log_gamma = T.scalar('log_gamma')
    log_lambda = T.scalar('log_lambda')
    prediction = (T.dot(T.nnet.relu((T.dot(X, w_1) + b_1)), w_2) + b_2)
    ' define the log posterior distribution '
    log_lik_data = ((((- 0.5) * X.shape[0]) * (T.log((2 * np.pi)) - log_gamma)) - ((T.exp(log_gamma) / 2) * T.sum(T.power((prediction - y), 2))))
    log_prior_data = ((((a0 - 1) * log_gamma) - (b0 * T.exp(log_gamma))) + log_gamma)
    log_prior_w = (((((((- 0.5) * (num_vars - 2)) * (T.log((2 * np.pi)) - log_lambda)) - ((T.exp(log_lambda) / 2) * ((((w_1 ** 2).sum() + (w_2 ** 2).sum()) + (b_1 ** 2).sum()) + (b_2 ** 2)))) + ((a0 - 1) * log_lambda)) - (b0 * T.exp(log_lambda))) + log_lambda)
    log_posterior = ((((log_lik_data * N) / X.shape[0]) + log_prior_data) + log_prior_w)
    (dw_1, db_1, dw_2, db_2, d_log_gamma, d_log_lambda) = T.grad(log_posterior, [w_1, b_1, w_2, b_2, log_gamma, log_lambda])
    logp_gradient = theano.function(inputs=[X, y, w_1, b_1, w_2, b_2, log_gamma, log_lambda, N], outputs=[dw_1, db_1, dw_2, db_2, d_log_gamma, d_log_lambda])
    self.nn_predict = theano.function(inputs=[X, w_1, b_1, w_2, b_2], outputs=prediction)
    '\n            Training with SVGD\n        '
    (X_train, y_train) = self.normalization(X_train, y_train)
    N0 = X_train.shape[0]
    ' initializing all particles '
    for i in range(self.M):
        (w1, b1, w2, b2, loggamma, loglambda) = self.init_weights(a0, b0)
        ridx = np.random.choice(range(X_train.shape[0]), np.min([X_train.shape[0], 1000]), replace=False)
        y_hat = self.nn_predict(X_train[ridx, :], w1, b1, w2, b2)
        loggamma = (- np.log(np.mean(np.power((y_hat - y_train[ridx]), 2))))
        self.theta[i, :] = self.pack_weights(w1, b1, w2, b2, loggamma, loglambda)
    grad_theta = np.zeros([self.M, num_vars])
    fudge_factor = 1e-06
    historical_grad = 0
    self.y_historical_grad = 0
    self.h_historical_grad = 0
    self.rmse_overTime = np.zeros(numTimeSteps)
    self.llh_overTime = np.zeros(numTimeSteps)
    self.iter_overTime = np.zeros(numTimeSteps)
    timeStepUnit = (maxTime / numTimeSteps)
    timeInd = 0
    start_time = time.time()
    for iter in range(max_iter):
        if (method == 'subparticles'):
            self.Sqy = np.zeros([m, num_vars])
        elif ((method == 'inducedPoints') or (method == 'none')):
            self.Sqx = np.zeros([self.M, num_vars])
        h = (- 1)
        batch = [(i % N0) for i in range((iter * batch_size), ((iter + 1) * batch_size))]
        if ((method == 'none') or (method == 'inducedPoints')):
            for i in range(self.M):
                (w1, b1, w2, b2, loggamma, loglambda) = self.unpack_weights(self.theta[i, :])
                (dw1, db1, dw2, db2, dloggamma, dloglambda) = logp_gradient(X_train[batch, :], y_train[batch], w1, b1, w2, b2, loggamma, loglambda, N0)
                self.Sqx[i, :] = self.pack_weights(dw1, db1, dw2, db2, dloggamma, dloglambda)
            if (method == 'none'):
                grad_theta = self.svgd_kernel(h=h)
            elif (method == 'inducedPoints'):
                self.yInd = np.random.choice(self.theta.shape[0], m, replace=False)
                self.y = self.theta[self.yInd]
                grad_theta = self.svgd_kernel_inducedPoints(h=h, uStat=uStat, regCoeff=regCoeff, adver=adver, adverMaxIter=adverMaxIter)
        elif (method == 'subparticles'):
            self.yInd = np.random.choice(self.theta.shape[0], m, replace=False)
            self.y = self.theta[self.yInd]
            for i in range(m):
                (w1, b1, w2, b2, loggamma, loglambda) = self.unpack_weights(self.y[i, :])
                (dw1, db1, dw2, db2, dloggamma, dloglambda) = logp_gradient(X_train[batch, :], y_train[batch], w1, b1, w2, b2, loggamma, loglambda, N0)
                self.Sqy[i, :] = self.pack_weights(dw1, db1, dw2, db2, dloggamma, dloglambda)
            grad_theta = self.svgd_kernel_subset(h=(- 1), cf=cf)
        [adj_grad, historical_grad] = self.get_adamUpdate(iter, grad_theta, historical_grad, master_stepsize, alpha, fudge_factor)
        self.theta = (self.theta + adj_grad)
        elapsed_time = (time.time() - start_time)
        if (elapsed_time > timeStepUnit):
            self.thetaCopy = np.copy(self.theta)
            '\n                    Model selection by using a development set\n                '
            X_dev = self.normalization(X_dev)
            for i in range(self.M):
                (w1, b1, w2, b2, loggamma, loglambda) = self.unpack_weights(self.thetaCopy[i, :])
                pred_y_dev = ((self.nn_predict(X_dev, w1, b1, w2, b2) * self.std_y_train) + self.mean_y_train)

                def f_log_lik(loggamma):
                    return np.sum(np.log(((np.sqrt(np.exp(loggamma)) / np.sqrt((2 * np.pi))) * np.exp((((- 1) * (np.power((pred_y_dev - y_dev), 2) / 2)) * np.exp(loggamma))))))
                lik1 = f_log_lik(loggamma)
                loggamma = (- np.log(np.mean(np.power((pred_y_dev - y_dev), 2))))
                lik2 = f_log_lik(loggamma)
                if (lik2 > lik1):
                    self.thetaCopy[(i, (- 2))] = loggamma
            (svgd_rmse, svgd_ll) = self.evaluation(X_test, y_test)
            self.rmse_overTime[timeInd] = svgd_rmse
            self.llh_overTime[timeInd] = svgd_ll
            self.iter_overTime[timeInd] = iter
            start_time = time.time()
            timeInd = (timeInd + 1)
            if (timeInd >= numTimeSteps):
                print('Reached ', iter, 'iterations\n')
                break
def __init__(self, X_train, y_train, X_test, y_text, batch_size=100, max_iter=1000, M=20, n_hidden=50, a0=1, b0=0.1, master_stepsize=0.001, auto_corr=0.9, h=(- 1), alpha=0.9, method='none', m=5, cf=False, uStat=True, regCoeff=0.1, adver=False, adverMaxIter=5, maxTime=20, numTimeSteps=20):
    self.n_hidden = n_hidden
    self.d = X_train.shape[1]
    self.M = M
    num_vars = (((self.d * n_hidden) + (n_hidden * 2)) + 3)
    self.theta = np.zeros([self.M, num_vars])
    '\n            We keep the last 10% (maximum 500) of training data points for model developing\n        '
    size_dev = min(int(np.round((0.1 * X_train.shape[0]))), 500)
    (X_dev, y_dev) = (X_train[(- size_dev):], y_train[(- size_dev):])
    (X_train, y_train) = (X_train[:(- size_dev)], y_train[:(- size_dev)])
    '\n            The data sets are normalized so that the input features and the targets have zero mean and unit variance\n        '
    self.std_X_train = np.std(X_train, 0)
    self.std_X_train[(self.std_X_train == 0)] = 1
    self.mean_X_train = np.mean(X_train, 0)
    self.mean_y_train = np.mean(y_train)
    self.std_y_train = np.std(y_train)
    '\n            Theano symbolic variables\n            Define the neural network here\n        '
    X = T.matrix('X')
    y = T.vector('y')
    w_1 = T.matrix('w_1')
    b_1 = T.vector('b_1')
    w_2 = T.vector('w_2')
    b_2 = T.scalar('b_2')
    N = T.scalar('N')
    log_gamma = T.scalar('log_gamma')
    log_lambda = T.scalar('log_lambda')
    prediction = (T.dot(T.nnet.relu((T.dot(X, w_1) + b_1)), w_2) + b_2)
    ' define the log posterior distribution '
    log_lik_data = ((((- 0.5) * X.shape[0]) * (T.log((2 * np.pi)) - log_gamma)) - ((T.exp(log_gamma) / 2) * T.sum(T.power((prediction - y), 2))))
    log_prior_data = ((((a0 - 1) * log_gamma) - (b0 * T.exp(log_gamma))) + log_gamma)
    log_prior_w = (((((((- 0.5) * (num_vars - 2)) * (T.log((2 * np.pi)) - log_lambda)) - ((T.exp(log_lambda) / 2) * ((((w_1 ** 2).sum() + (w_2 ** 2).sum()) + (b_1 ** 2).sum()) + (b_2 ** 2)))) + ((a0 - 1) * log_lambda)) - (b0 * T.exp(log_lambda))) + log_lambda)
    log_posterior = ((((log_lik_data * N) / X.shape[0]) + log_prior_data) + log_prior_w)
    (dw_1, db_1, dw_2, db_2, d_log_gamma, d_log_lambda) = T.grad(log_posterior, [w_1, b_1, w_2, b_2, log_gamma, log_lambda])
    logp_gradient = theano.function(inputs=[X, y, w_1, b_1, w_2, b_2, log_gamma, log_lambda, N], outputs=[dw_1, db_1, dw_2, db_2, d_log_gamma, d_log_lambda])
    self.nn_predict = theano.function(inputs=[X, w_1, b_1, w_2, b_2], outputs=prediction)
    '\n            Training with SVGD\n        '
    (X_train, y_train) = self.normalization(X_train, y_train)
    N0 = X_train.shape[0]
    ' initializing all particles '
    for i in range(self.M):
        (w1, b1, w2, b2, loggamma, loglambda) = self.init_weights(a0, b0)
        ridx = np.random.choice(range(X_train.shape[0]), np.min([X_train.shape[0], 1000]), replace=False)
        y_hat = self.nn_predict(X_train[ridx, :], w1, b1, w2, b2)
        loggamma = (- np.log(np.mean(np.power((y_hat - y_train[ridx]), 2))))
        self.theta[i, :] = self.pack_weights(w1, b1, w2, b2, loggamma, loglambda)
    grad_theta = np.zeros([self.M, num_vars])
    fudge_factor = 1e-06
    historical_grad = 0
    self.y_historical_grad = 0
    self.h_historical_grad = 0
    self.rmse_overTime = np.zeros(numTimeSteps)
    self.llh_overTime = np.zeros(numTimeSteps)
    self.iter_overTime = np.zeros(numTimeSteps)
    timeStepUnit = (maxTime / numTimeSteps)
    timeInd = 0
    start_time = time.time()
    for iter in range(max_iter):
        if (method == 'subparticles'):
            self.Sqy = np.zeros([m, num_vars])
        elif ((method == 'inducedPoints') or (method == 'none')):
            self.Sqx = np.zeros([self.M, num_vars])
        h = (- 1)
        batch = [(i % N0) for i in range((iter * batch_size), ((iter + 1) * batch_size))]
        if ((method == 'none') or (method == 'inducedPoints')):
            for i in range(self.M):
                (w1, b1, w2, b2, loggamma, loglambda) = self.unpack_weights(self.theta[i, :])
                (dw1, db1, dw2, db2, dloggamma, dloglambda) = logp_gradient(X_train[batch, :], y_train[batch], w1, b1, w2, b2, loggamma, loglambda, N0)
                self.Sqx[i, :] = self.pack_weights(dw1, db1, dw2, db2, dloggamma, dloglambda)
            if (method == 'none'):
                grad_theta = self.svgd_kernel(h=h)
            elif (method == 'inducedPoints'):
                self.yInd = np.random.choice(self.theta.shape[0], m, replace=False)
                self.y = self.theta[self.yInd]
                grad_theta = self.svgd_kernel_inducedPoints(h=h, uStat=uStat, regCoeff=regCoeff, adver=adver, adverMaxIter=adverMaxIter)
        elif (method == 'subparticles'):
            self.yInd = np.random.choice(self.theta.shape[0], m, replace=False)
            self.y = self.theta[self.yInd]
            for i in range(m):
                (w1, b1, w2, b2, loggamma, loglambda) = self.unpack_weights(self.y[i, :])
                (dw1, db1, dw2, db2, dloggamma, dloglambda) = logp_gradient(X_train[batch, :], y_train[batch], w1, b1, w2, b2, loggamma, loglambda, N0)
                self.Sqy[i, :] = self.pack_weights(dw1, db1, dw2, db2, dloggamma, dloglambda)
            grad_theta = self.svgd_kernel_subset(h=(- 1), cf=cf)
        [adj_grad, historical_grad] = self.get_adamUpdate(iter, grad_theta, historical_grad, master_stepsize, alpha, fudge_factor)
        self.theta = (self.theta + adj_grad)
        elapsed_time = (time.time() - start_time)
        if (elapsed_time > timeStepUnit):
            self.thetaCopy = np.copy(self.theta)
            '\n                    Model selection by using a development set\n                '
            X_dev = self.normalization(X_dev)
            for i in range(self.M):
                (w1, b1, w2, b2, loggamma, loglambda) = self.unpack_weights(self.thetaCopy[i, :])
                pred_y_dev = ((self.nn_predict(X_dev, w1, b1, w2, b2) * self.std_y_train) + self.mean_y_train)

                def f_log_lik(loggamma):
                    return np.sum(np.log(((np.sqrt(np.exp(loggamma)) / np.sqrt((2 * np.pi))) * np.exp((((- 1) * (np.power((pred_y_dev - y_dev), 2) / 2)) * np.exp(loggamma))))))
                lik1 = f_log_lik(loggamma)
                loggamma = (- np.log(np.mean(np.power((pred_y_dev - y_dev), 2))))
                lik2 = f_log_lik(loggamma)
                if (lik2 > lik1):
                    self.thetaCopy[(i, (- 2))] = loggamma
            (svgd_rmse, svgd_ll) = self.evaluation(X_test, y_test)
            self.rmse_overTime[timeInd] = svgd_rmse
            self.llh_overTime[timeInd] = svgd_ll
            self.iter_overTime[timeInd] = iter
            start_time = time.time()
            timeInd = (timeInd + 1)
            if (timeInd >= numTimeSteps):
                print('Reached ', iter, 'iterations\n')
                break