Merge pull request #13 from cumc/master

sync

multivariate-fine-mapping/mixture_prior.ipynb

@@ -71,6 +71,12 @@
  " --datadir $m --name `basename $m` &> bovy_$m.log\n",
  "```\n",
  "\n",
+ "### Plot results\n",
+ "\n",
+ "```\n",
+ "sos run mixture_prior.ipynb plot_U --model-data /path/to/mixture_model.rds --cwd output\n",
+ "```\n",
+ "\n",
  "Notice that for production use, each `sos run` command should be submitted to the cluster as a job."
  ]
  },
@@ -452,6 +458,8 @@
  "# A list of models where we only update the scales and not the matrices\n",
  "# A typical choice is to estimate scales only for canonical components\n",
  "parameter: scale_only = []\n",
+ "# Tolerance for change in likelihood\n",
+ "parameter: ud_tol_lik = 1e-3\n",
  "input: [f\"{cwd}/{name}.rds\"] + [f\"{cwd}/{name}.{m}.rds\" for m in mixture_components]\n",
  "output: f'{cwd}/{name}.{ud_method}{\"_unconstrained\" if len(scale_only) == 0 else \"\"}{(\".\" + os.path.basename(resid_cor)[:-4]) if resid_cor.is_file() else \"\"}.rds'\n",
  "task: trunk_workers = 1, walltime = '36h', trunk_size = 1, mem = '10G', cores = 4, tags = f'{_output:bn}'\n",
@@ -481,7 +489,7 @@
  " library(udr)\n",
  " message(paste(\"Running ${ud_method.upper()} via udr package for\", length(U), \"mixture components\"))\n",
  " f0 = ud_init(X = as.matrix(dat$strong.z), V = V, U_scaled = U_scaled, U_unconstrained = U, n_rank1=0)\n",
- " res = ud_fit(f0, X = na.omit(f0$X), control = list(unconstrained.update = \"${ud_method}\", resid.update = 'none', scaled.update = \"fa\", maxiter=5000, tol = 1e-06), verbose=TRUE)\n",
+ " res = ud_fit(f0, X = na.omit(f0$X), control = list(unconstrained.update = \"${ud_method}\", resid.update = 'none', scaled.update = \"fa\", maxiter=5000, tol.lik = ${ud_tol_lik}), verbose=TRUE)\n",
  " # add back col and row names to U\n",
  " # https://github.com/stephenslab/udr/issues/9\n",
  " for (i in 1:length(res$U)) {\n",
@@ -499,6 +507,7 @@
  "outputs": [],
  "source": [
  "[ed]\n",
+ "parameter: ed_tol = 1e-6\n",
  "input: [f\"{cwd}/{name}.rds\"] + [f\"{cwd}/{name}.{m}.rds\" for m in mixture_components]\n",
  "output: f'{cwd}/{name}.ed_bovy{(\".\" + os.path.basename(resid_cor)[:-4]) if resid_cor.is_file() else \"\"}.rds'\n",
  "task: trunk_workers = 1, walltime = '36h', trunk_size = 1, mem = '10G', cores = 4, tags = f'{_output:bn}'\n",
@@ -515,7 +524,7 @@
  " # Fit mixture model using ED code by J. Bovy\n",
  " mash_data = mashr::mash_set_data(dat$strong.z, V=V)\n",
  " message(paste(\"Running ED via J. Bovy's code for\", length(U), \"mixture components\"))\n",
- " res = mashr:::bovy_wrapper(mash_data, U, logfile=${_output:nr}, tol = 1e-06)\n",
+ " res = mashr:::bovy_wrapper(mash_data, U, logfile=${_output:nr}, tol = ${ed_tol})\n",
  " saveRDS(list(U=res$Ulist, w=res$pi, loglik=scan(\"${_output:n}_loglike.log\")), ${_output:r})"
  ]
  },