Add clockpr option to specify prior on substitution rate

Author: 4ment

phylostan/generate_script.py

@@ -1010,8 +1010,11 @@ def get_model(params):
 		if params.estimate_rate:
 			if params.clock == 'strict':
 				parameters_block.append('real <lower=0> rate;')
-				functions_block.append(ctmc_scale_prior())
-				model_priors.append('rate ~ ctmc_scale(blensUnscaled);')
+				if params.clockpr == 'ctmcscale':
+					functions_block.append(ctmc_scale_prior())
+					model_priors.append('rate ~ ctmc_scale(blensUnscaled);')
+				else:
+					model_priors.append('rate ~ exponential(1000);')
 			elif params.clock.endswith('mrf'):
 				transformed_parameters_declarations.append('real substrates[bcount];')
 				parameters_block.append('real deltas[2*S-3];')
@@ -1189,7 +1192,7 @@ def get_model(params):
 		raise ValueError('Supports JC69, HKY and GTR only.')
 
 	# Tree likelihood
-	model_block.append(likelihood(params.categories > 1 or params.invariant, params.clock is not None))
+	model_block.append(likelihood2(params.categories > 1 or params.invariant, params.clock is not None))
 
 	if params.clock is not None:
 		model_block.append(jacobian(params.heterochronous))

phylostan/phylostan.py

@@ -82,6 +82,10 @@ def create_build_parser(subprasers, prog, help):
 	parser.add_argument('--clock', required=False,
 						choices=['strict', 'ace', 'acln', 'acg', 'aoup', 'ucln', 'uced', 'gmrf', 'hsmrf'], default=None,
 						help="""Type of clock""")
+	parser.add_argument('--clockpr', default='ctmcscale',
+						type=lambda x: distribution_type(x, ('exponential', 'ctmcscale')),
+						help="""prior on substitution rate [default: %(default)s]""",
+	)
 	parser.add_argument('--estimate_rate', action='store_true', help="""Estimate substitution rate""")
 	parser.add_argument('-c', '--coalescent', choices=['constant', 'skyride', 'skygrid'], default=None,
 						help="""Type of coalescent (constant or skyride)""")
@@ -105,6 +109,22 @@ def create_compile_parser(subprasers, prog, help):
 	return parser
 
 
+def distribution_type(arg, choices):
+    """Used by argparse for specifying distributions with optional
+    parameters."""
+    res = arg.split('(')
+    if (isinstance(choices, tuple) and res[0] in choices) or res[0] == choices:
+        return arg
+    else:
+        if isinstance(choices, tuple):
+            message = "'" + "','".join(choices) + '"'
+        else:
+            message = "'" + choices + "'"
+        raise argparse.ArgumentTypeError(
+            'invalid choice (choose from a number or ' + message + ')'
+        )
+
+
 def main():
 	parser_epilog = """To get some information for each sub-command:\n
 	phylostan build --help