Major changes

inst/pydevil/.virtual_documents/notebook_test/batch_vs_full.ipynb

@@ -0,0 +1,132 @@
+import scanpy as sc
+import patsy as ptsy
+import torch
+import matplotlib.pyplot as plt
+from scipy.optimize import curve_fit
+import numpy as np
+import pandas as pd
+
+
+data_path = "../../../../../../../Dropbox/2021. CLL CD49/scRNA/data/SevenBridges_counts/Novaseq_SampleTag01_hs_RM238_2/Novaseq_SampleTag01_hs_RM238_2_DBEC_MolsPerCell.csv"
+labels_path = "../../../../../../../Dropbox/2021. CLL CD49/scRNA/data/kmeans_labels/Novaseq_SampleTag01_hs_RM238_2_Ab_Labels.tsv"
+labels = pd.read_csv(labels_path, sep="\t")
+adata = sc.read_csv(data)
+
+
+adata = adata[:, 5:]
+adata.obs["cell_type"] = list(labels['labels'])
+
+
+adata = adata[:, adata.X.sum(axis=0) > 0]
+
+
+import sys
+sys.path.append("../")
+%reload_ext autoreload
+%autoreload 2
+
+import pydevil
+
+
+covariates = ptsy.dmatrix("~ 1", adata.obs)
+UMI = np.array(X.sum(1)).squeeze() / np.mean(np.array(X.sum(1)).squeeze())
+gene_names = adata.var_names
+obs_names = adata.obs_names
+X = torch.tensor(adata.X).int()
+
+
+res  = pydevil.run_SVDE(X, covariates, UMI,
+                        gene_names = gene_names,
+                        cell_names = obs_names, 
+                        #group_matrix = random_effects,
+                        jit_compile=False,
+                        optimizer_name = "ClippedAdam", lr = 0.5, gamma_lr=0.0001, steps = 20, 
+                        full_cov = True, cuda = True, batch_size=X.shape[0])
+
+
+res_batch  = pydevil.run_SVDE(X, covariates, UMI,
+                        gene_names = gene_names,
+                        cell_names = obs_names, 
+                        #group_matrix = random_effects,
+                        jit_compile=False,
+                        optimizer_name = "ClippedAdam", lr = 0.5, gamma_lr=0.0001, steps = 20, 
+                        full_cov = False, cuda = True, batch_size=int(X.shape[0] / 4.0))
+
+
+contrast = np.array([1,-1])
+
+
+betas = res_nocov['params']['beta']
+var = res_nocov['params']['variance']
+
+
+(betas * contrast.reshape([-1,1])).sum(axis = 0), (var**2 * abs(contrast)).sum(axis = 1)
+
+
+batch_factor = 2
+res_no_cov_batch  = pydevil.run_SVDE(X, covariates, UMI,
+                        gene_names = gene_names,
+                        cell_names = obs_names, 
+                        #group_matrix = random_effects,
+                        jit_compile=False,
+                        optimizer_name = "ClippedAdam", lr = 0.5 / batch_factor, gamma_lr=0.0001, steps = 500, 
+                        full_cov = True, cuda = True, batch_size=int(X.shape[1] / batch_factor))
+
+
+pydevil.test_posterior_null(res, contrast)
+
+
+pydevil.test_posterior_null(res_batch, contrast)
+
+
+pydevil.test_posterior_null(res, np.array([1]))
+pydevil.test_posterior_null(res_nocov, np.array([1]))
+pydevil.test_posterior_null(res_no_cov_batch, np.array([1]))
+
+
+theta = res['params']['theta']
+beta = res['params']['beta'].squeeze(0)
+beta_loc = res['params']['variance']
+
+theta_nocov = res_nocov['params']['theta']
+beta_nocov = res_nocov['params']['beta'].squeeze(0)
+beta_loc_nocov = res_nocov['params']['variance']
+
+theta_batch = res_no_cov_batch['params']['theta']
+beta_batch = res_no_cov_batch['params']['beta'].squeeze(0)
+beta_loc_bathc = res_no_cov_batch['params']['variance']
+
+
+fig, (ax1, ax2) = plt.subplots(1, 2)
+fig.suptitle('Plot for overdispersion')
+ax1.scatter(theta, theta_batch, s = .5)
+ax2.scatter([i for i in range(len(theta))], theta / theta_batch, s = .5)
+plt.show()
+
+
+fig, (ax1, ax2) = plt.subplots(1, 2)
+fig.suptitle('Plot for overdispersion')
+ax1.scatter(theta, theta_batch_nc, s = .5)
+ax2.scatter([i for i in range(len(theta))], theta / theta_batch_nc, s = .5)
+plt.show()
+
+
+fig, (ax1, ax2) = plt.subplots(1, 2)
+fig.suptitle('Plot for beta')
+ax1.scatter(beta, beta_batch, s = .5)
+ax2.scatter([i for i in range(len(beta))], beta / beta_batch, s = .5)
+
+
+fig, (ax1, ax2) = plt.subplots(1, 2)
+fig.suptitle('Plot for beta')
+ax1.scatter(beta_loc, beta_loc_batch, s = .5)
+ax2.scatter([i for i in range(len(beta_loc))], beta_loc_batch / beta_loc, s = .5)
+
+
+beta_loc_bathc_nc.shape, 
+
+
+fig, (ax1, ax2) = plt.subplots(1, 2)
+fig.suptitle('Plot for beta')
+ax1.scatter(beta_loc, beta_loc_bathc_nc, s = .5)
+ax2.scatter([i for i in range(len(beta_loc))], beta_loc_bathc_nc / beta_loc, s = .5)