- test:

- name: 'test with python ${{ matrix.python-version }}'

- runs-on: ubuntu-20.04

- strategy:

- fail-fast: false

- matrix:

- python-version:

- - '3.7'

- - '3.8'

- - '3.9'

- - '3.10'

- - '3.11'

-

- steps:

- - uses: actions/checkout@v3

- - name: Set up Python ${{ matrix.python-version }}

- uses: actions/setup-python@v4

- with:

- python-version: ${{ matrix.python-version }}

- - name: Install dependencies

- run: |

- python -m pip install --upgrade pip

- pip install pytest wheel

- - name: Install treetime

- run: python -m pip install .

- - name: Run tests

- run: bash test.sh

-ignore=.git .idea .input .output .temp .venv .vscode

- bad-continuation,

- bad-whitespace,

- too-many-branches,

- useless-object-inheritance,

-

- duplicate-code,

-

-

- simplifiable-if-expression,

- simplifiable-if-statement,

- singleton-comparison,

-

- attribute-defined-outside-init,

- multiple-statements,

- redefined-outer-name,

-

- cyclic-import,

- import-error,

- import-outside-toplevel,

- reimported,

-extension-pkg-whitelist=numpy

-export UID=$(shell id -u)

-export GID=$(shell id -g)

-

-export DOCS_CONTAINER_NAME=treetime-docs

-.PHONY: docs docker-docs

- set -euox

- from treetime.gtr import GTR

- dummy_prof = np.random.random(size=(10000,5))

-

- dates = parse_dates(base_name+'.metadata.csv')

- tt = TreeTime(gtr='Jukes-Cantor', tree = base_name+'.nwk', use_fft=True,

- aln = base_name+'.fasta', verbose = 3, dates = dates, debug=True)

- tt.run(root="best", branch_length_mode='input', max_iter=2, time_marginal=True)

- msgs_to_multiply = [child.joint_pos_Lx for child in tt.tree.root.clades

- if child.joint_pos_Lx is not None]

- res, res_t = NodeInterpolator.convolve(subtree_distribution, tt.tree.root.clades[1].branch_length_interpolator, max_or_integral='max', inverse_time=True, n_grid_points = tt.node_grid_points, n_integral=tt.n_integral, rel_tol=tt.rel_tol_refine)

-

- msgs_to_multiply = [child.marginal_pos_Lx for child in tt.tree.root.clades

- if child.marginal_pos_Lx is not None]

- res, res_t = NodeInterpolator.convolve(subtree_distribution, tt.tree.root.clades[1].branch_length_interpolator, max_or_integral='integral', inverse_time=True, n_grid_points = tt.node_grid_points, n_integral=tt.n_integral, rel_tol=tt.rel_tol_refine)

- res, res_t = NodeInterpolator.convolve(subtree_distribution, tt.tree.root.clades[1].branch_length_interpolator, max_or_integral='integral', inverse_time=False, n_grid_points = tt.node_grid_points, n_integral=tt.n_integral, rel_tol=tt.rel_tol_refine)

- res = NodeInterpolator.convolve_fft(subtree_distribution, tt.tree.root.clades[1].branch_length_interpolator, inverse_time=True)

- res = NodeInterpolator.convolve_fft(subtree_distribution, tt.tree.root.clades[1].branch_length_interpolator, inverse_time=False)

-

-

-

-Upon pushing a commit, travis will run a few simple tests. These use data available in the [neherlab/treetime_examples](https://github.com/neherlab/treetime_examples) repository.

-Uploaders: Andreas Tille <tille@debian.org>,

- Étienne Mollier <emollier@debian.org>

-Priority: optional

-Build-Depends: debhelper-compat (= 13),

- dh-sequence-python3,

- python3-all,

- python3-biopython,

- python3-numpy,

- python3-pandas,

- python3-scipy,

- python3-setuptools,

- debhelper

-Standards-Version: 4.7.0

-Rules-Requires-Root: no

-Depends: ${python3:Depends},

- ${misc:Depends},

- python3-numpy,

- python3-scipy,

- python3-biopython,

- python3-pandas

-Recommends: python3-matplotlib

-syntaxwarning.patch

-Description: fix SyntaxWarnings in treetime/merger_models.py.

- Embedded TeX mathematical snippets are causing the below warnings:

- .

- /usr/lib/python3/dist-packages/treetime/merger_models.py:187:

- SyntaxWarning: invalid escape sequence '\k'

- '''

- /usr/lib/python3/dist-packages/treetime/merger_models.py:197:

- SyntaxWarning: invalid escape sequence '\l'

- '''

- /usr/lib/python3/dist-packages/treetime/merger_models.py:209:

- SyntaxWarning: invalid escape sequence '\l'

- '''

- .

- Switching the affected strings to raw strings with prefix r resolves

- the issue first, and second allows getting back our τ characters with

- single slash; it is assumed that the duplicate slash was necessary

- because Python's regular string would otherwise throw a <Tab>au instead

- of a \tau.

-Author: Étienne Mollier <emollier@debian.org>

-Bug-Debian: https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=1086965

-Forwarded: https://github.com/neherlab/treetime/pull/285

-Last-Update: 2024-11-09

----

-This patch header follows DEP-3: http://dep.debian.net/deps/dep3/

---- python-treetime.orig/treetime/merger_models.py

-+++ python-treetime/treetime/merger_models.py

-@@ -184,7 +184,7 @@

- np.concatenate(([cost[0]], cost,[cost[-1]])), kind='linear')

-

- def branch_merger_rate(self, t):

-- '''

-+ r'''

- rate at which one particular branch merges with any other branch at time t,

- in the Kingman model this is: :math:`\kappa(t) = (k(t)-1)/(2Tc(t))`

- '''

-@@ -194,7 +194,7 @@

- return 0.5*np.maximum(0.5,self.nbranches(t)-1.0)/self.Tc(t)

-

- def total_merger_rate(self, t):

-- '''

-+ r'''

- rate at which any branch merges with any other branch at time t,

- in the Kingman model this is: :math:`\lambda(t) = k(t)(k(t)-1)/(2Tc(t))`

- '''

-@@ -206,12 +206,12 @@

-

-

- def cost(self, t_node, branch_length, multiplicity=2.0):

-- '''

-+ r'''

- returns the cost associated with a branch starting with divergence time t_node (:math:`t_n`)

-- having a branch length :math:`\\tau`.

-+ having a branch length :math:`\tau`.

- This is equivalent to the probability of there being no merger on that branch and a merger at the end of the branch,

- calculated in the negative log

-- :math:`-log(\lambda(t_n+ \\tau)^{(m-1)/m}) + \int_{t_n}^{t_n+ \\tau} \kappa(t) dt`, where m is the multiplicity

-+ :math:`-log(\lambda(t_n+ \tau)^{(m-1)/m}) + \int_{t_n}^{t_n+ \tau} \kappa(t) dt`, where m is the multiplicity

-

- Parameters

- ----------

-override_dh_install:

- dh_install

-override_dh_auto_clean:

- dh_auto_clean

- rm -rf phylo_treetime.egg-info

-version=4

-opts="uversionmangle=s/-beta/~beta/,filenamemangle=s%(?:.*?)?v?(\d[\d.]*)\.tar\.gz%@PACKAGE@-$1.tar.gz%" \

-https://github.com/neherlab/treetime/tags (?:.*?/)?v?@ANY_VERSION@@ARCHIVE_EXT@

-#needs_sphinx = '1.0'

- 'sphinxarg.ext'

-#source_encoding = 'utf-8-sig'

-project = u'TreeTime'

-copyright = u'2017-2021, Pavel Sagulenko and Richard Neher'

-author = u'Pavel Sagulenko and Richard Neher'

-#today = ''

-#today_fmt = '%B %d, %Y'

-#default_role = None

-#add_function_parentheses = True

-#add_module_names = True

-#show_authors = False

-#modindex_common_prefix = []

-#keep_warnings = False

-#html_theme_options = {}

-#html_theme_path = []

-#html_title = u'TreeTime v1.0'

-#html_short_title = None

-#html_logo = None

-#html_favicon = None

-#html_extra_path = []

-#html_last_updated_fmt = None

-#html_use_smartypants = True

-#html_sidebars = {}

-#html_additional_pages = {}

-#html_domain_indices = True

-#html_use_index = True

-#html_split_index = False

-#html_show_sourcelink = True

-#html_show_sphinx = True

-#html_show_copyright = True

-#html_use_opensearch = ''

-#html_file_suffix = None

-#html_search_language = 'en'

-#html_search_options = {'type': 'default'}

-#html_search_scorer = 'scorer.js'

-# The paper size ('letterpaper' or 'a4paper').

-#'papersize': 'letterpaper',

-

-# The font size ('10pt', '11pt' or '12pt').

-#'pointsize': '10pt',

-

-# Additional stuff for the LaTeX preamble.

-#'preamble': '',

-

-# Latex figure (float) alignment

-#'figure_align': 'htbp',

- (master_doc, 'TreeTime.tex', u'TreeTime Documentation',

- u'Pavel Sagulenko and Richard Neher', 'manual'),

-#latex_logo = None

-#latex_use_parts = False

-#latex_show_pagerefs = False

-#latex_show_urls = False

-#latex_appendices = []

-#latex_domain_indices = True

-man_pages = [

- (master_doc, 'treetime', u'TreeTime Documentation',

- [author], 1)

-]

-#man_show_urls = False

- (master_doc, 'TreeTime', u'TreeTime Documentation',

- author, 'TreeTime', 'One line description of project.',

- 'Miscellaneous'),

-#texinfo_appendices = []

-#texinfo_domain_indices = True

-#texinfo_show_urls = 'footnote'

-#texinfo_no_detailmenu = False

-.. image:: https://anaconda.org/bioconda/treetime/badges/installer/conda.svg

- v = "0.0.0"

- if line[:7]=='version':

-with open("README.md", "r") as fh:

- name = "phylo-treetime",

- version = get_version(),

- author = "Pavel Sagulenko, Emma Hodcroft, and Richard Neher",

- author_email = "richard.neher@unibas.ch",

- description = ("Maximum-likelihood phylodynamic inference"),

- long_description = long_description,

- long_description_content_type="text/markdown",

- license = "MIT",

- keywords = "Time-stamped phylogenies, phylogeography, virus evolution",

- url = "https://github.com/neherlab/treetime",

- packages=['treetime'],

- install_requires = [

- 'biopython>=1.67,!=1.77,!=1.78',

- 'numpy>=1.10.4',

- 'pandas>=0.17.1',

- 'scipy>=0.16.1'

- ],

- extras_require = {

- ':python_version >= "3.6"':['matplotlib>=2.0'],

- },

- classifiers=[

- "Development Status :: 5 - Production/Stable",

- "Topic :: Scientific/Engineering :: Bio-Informatics",

- "License :: OSI Approved :: MIT License",

- "Programming Language :: Python :: 3.7",

- "Programming Language :: Python :: 3.8",

- "Programming Language :: Python :: 3.9",

- "Programming Language :: Python :: 3.10",

- ],

- entry_points = {

- "console_scripts": [

- "treetime = treetime.__main__:main",

- ]

- }

- )

-

-def write_to_file(times_and_names, name = 'master.txt'):

- f.write("time \t name \t bad_branch\n")

- f.write("\t")

- f.write("\t")

- f.write("\n")

- lines = line.split("\n")[0].split("\t")

-

- ebola=True

- master = False ##should be True for master branch

- location_master = '../../TreeTimeMaster/treetime/' ##only needed for test branch

- seq_kwargs = {"marginal_sequences":True,

- "branch_length_mode": 'input',

- "sample_from_profile":"root",

- "reconstruct_tip_states":False}

- tt_kwargs = {'clock_rate':clock_rate,

- 'time_marginal':'assign'}

- coal_kwargs ={'Tc':10000,

- 'time_marginal':'assign'}

-

- dates = parse_dates(base_name+'.metadata.csv')

- tt = TreeTime(gtr='Jukes-Cantor', tree = base_name+'.nwk', use_fft=False,

- aln = base_name+'.fasta', verbose = 1, dates = dates, precision=3, debug=True)

-

- tt._set_branch_length_mode(seq_kwargs["branch_length_mode"])

- tt.infer_ancestral_sequences(infer_gtr=False, marginal=seq_kwargs["marginal_sequences"])

- tt.clock_filter(reroot='least-squares', n_iqd=1, plot=False, fixed_clock_rate=tt_kwargs["clock_rate"])

- tt.reroot(root='least-squares', clock_rate=tt_kwargs["clock_rate"])

- tt.make_time_tree(clock_rate=tt_kwargs["clock_rate"], time_marginal=tt_kwargs["time_marginal"])

- output_comparison = compare(tree_events_tt, read_from_file(location_master + "master.txt"))

- plt.plot(output_comparison[1][output_comparison[1]['bad_branch']==0].time, output_comparison[1][output_comparison[1]['bad_branch']==0].difference, 'o')

- tt.add_coalescent_model(coal_kwargs ["Tc"])

- tt.make_time_tree(clock_rate=tt_kwargs ["clock_rate"], time_marginal=coal_kwargs ["time_marginal"])

- write_to_file(tree_events_tt, "fft_branch"+ ".txt")

- output_comparison = compare(tree_events_tt_post_coal, read_from_file(location_master + "master.txt"))

- plt.plot(output_comparison[1].time, output_comparison[1].difference, marker='o', color=color[int(name)], linestyle='', ms=1, label=name)

- plt.xlabel("nodes ranging from root at 0 to most recent")

- plt.ylabel("difference time_before_present coalescent branch - master branch")

- plt.savefig("time_before_present-differences-ebola.png")

- plt.savefig("time_before_present-differences-h3n2_na.png")

- plt.plot(output_comparison[1][output_comparison[1]['bad_branch']==0].time, output_comparison[1][output_comparison[1]['bad_branch']==0].difference, 'o')

- Phylo.draw(tt.tree, label_func=lambda x:"")

- if coal_kwargs['time_marginal']=='assign':

- t = np.linspace(test_node.time_before_present-0.0001,test_node.time_before_present+0.0001,1000)

- test_node= test_node.up \ No newline at end of file

- "x = np.linspace(start, end, num=1000, endpoint=False) #create an x vector of points that should be sampled\n",

- "y_1 = gamma.pdf(x, alpha_1, scale=1/beta)\n",

- "y_1[y_1 == 0] = TINY_NUMBER ##remove -inf to prevent numerical errors\n",

- "fig= plt.figure()\n",

- "plt.plot(x, gamma.pdf(x, alpha_1, scale=1/beta), label=\"PDF\")\n",

- "plt.ylabel(\"PDF Gamma(\"+ str(alpha_1) +\",\"+str(beta)+ \")\")\n",

- "plt.title(\"Distribution 1\")\n",

- "y_2 = gamma.pdf(x, alpha_2, scale=1/beta)\n",

- "y_2[y_2 == 0] = TINY_NUMBER ##remove -inf to prevent numerical errors\n",

- "plt.plot(x, gamma.pdf(x, alpha_2, scale=1/beta), color=\"blue\")\n",

- "plt.plot(x, np.exp(-dist_2.y), color=\"green\") ##just to check that distribution object is initialized accurately \n",

- "plt.ylabel(\"PDF Gamma(\"+ str(alpha_2) +\",\"+str(beta)+ \")\")\n",

- "plt.title(\"Distribution 2\")\n",

- "y_sol = gamma.pdf(x, alpha_sol, scale=1/beta)\n",

- "y_sol[y_sol == 0] = TINY_NUMBER ##remove -inf to prevent numerical errors\n",

- "plt.plot(x, gamma.pdf(x, alpha_sol, scale=1/beta))\n",

- "plt.ylabel(\"PDF Gamma(\"+ str(alpha_sol) +\",\"+str(beta)+ \")\")\n",

- "plt.title(\"Convolution Distribution 1 and Distribution 2\")\n",

- " '''\n",

- " '''\n",

- " if calc_type==\"fft\":\n",

- " if alpha_1==1 or alpha_2==1:\n",

- " res = NodeInterpolator.convolve_fft(dist_1, dist_2, fft_grid_size=int(grid_size/factor))\n",

- " if calc_type==\"num\":\n",

- " res = NodeInterpolator.convolve(dist_1, dist_2, n_grid_points=grid_size, n_integral=grid_size)[0] ##second output is grid scale, which is also given in the x parameter\n",

- " '''\n",

- " '''\n",

- " return y\n",

- "\n",

- " "

- "grid_size = [100, 200, 300, 500, 1000] #more desired grid_sizes, true grid size is given as output when the accuracy_and_time function is called\n",

- "time_fft = np.empty((len(grid_size),2))\n",

- "time_num = np.empty((len(grid_size),2))\n",

- " #print(t)\n",

- " time_fft_full = accuracy_and_time(dist_1, dist_2, \"fft\", grid_size[t])\n",

- " time_num_full = accuracy_and_time(dist_1, dist_2, \"num\", grid_size[t])\n",

- " accuracy_num.append(time_num_full[0])\n"

- " print(\"grid points NUM (output):\" +str(time_list_num[t][1]))\n",

- " norm_y_est_num = get_P(x, accuracy_list_num[t] )\n",

- " print(\"grid points FFT (output):\" +str(time_list_fft[t][1]))\n",

- " norm_y_est_fft = get_P(x, accuracy_list_fft[t] )\n",

- " #plt.autoscale()\n",

- " #plt.subplot(1,2,1)\n",

- " plt.plot(x, norm_y_sol, color= \"green\", label=\"true dist\")\n",

- " plt.plot(x, norm_y_est_num, color=\"blue\", label=\"est num dist\")\n",

- " plt.plot(x, norm_y_est_fft, color=\"red\", label=\"est fft dist\")\n",

- " plt.axvline(x=eff_support[0], label= \"effective support start\")\n",

- " plt.axvline(x=eff_support[1], label= \"effective support end\")\n",

- " plt.ylabel(\"P\", fontsize=7)\n",

- " plt.xlim((0,1))\n",

- " #plt.subplot(1,2,2)\n",

- " plt.plot(x, norm_y_sol, color= \"green\", label=\"true dist\")\n",

- " plt.plot(x, norm_y_est_num, color=\"blue\", label=\"est num dist\")\n",

- " plt.plot(x, norm_y_est_fft, color=\"red\", label=\"est fft dist\")\n",

- " plt.axvline(x=eff_support[0], label= \"effective support start\")\n",

- " plt.axvline(x=eff_support[1], label= \"effective support end\")\n",

- " plt.xlim((0,1))\n",

- " plt.yscale(\"log\")\n",

- " plt.ylabel(\"log(P)\", fontsize=7)\n",

- " print(\"FFT grid points (output):\" +str(time_fft[t][1]))\n",

- " print(\"NUM grid points (output):\" +str(time_num[t][1]))\n",

- " plt.plot(x[2:], ((norm_y_est_fft -norm_y_sol)/(norm_y_sol+epsilon))[2:], 'o', markersize=1, color=\"blue\", label=\"normalized diff FFT\")\n",

- " plt.plot(x, (norm_y_est_num -norm_y_sol)/(norm_y_sol+epsilon), 'o', markersize=1, color=\"red\", label=\"normalized diff NUM\")\n",

- " plt.axvline(x=eff_support[0], label= \"effective support start\")\n",

- " plt.axvline(x=eff_support[1], label= \"effective support end\")\n",

- " plt.legend(prop={\"size\":10})\n",

- " plt.xlim((0,1))\n",

- " #plt.yscale(\"log\")\n",

- " plt.ylabel(\"(calculation - analytical_sol)/(analytical_sol + epsilon)\", fontsize=7)\n",

- " plt.xlabel(\"Time\", fontsize=7)\n",

- " plt.title(\"Normalized Difference Estimated and True Solution\", fontsize=12)\n",

- " x_support = np.linspace(eff_support[1]-0.1, eff_support[1]+0.1, num=1000, endpoint=False) #create an x vector of points that should be sampled\n",

- " norm_y_est_num = get_P(x_support, accuracy_num[t]) \n",

- " #print(accuracy_fft[t].x)\n",

- " plt.plot(x_support[2:], ((norm_y_est_fft -norm_y_sol_sup)/(norm_y_sol_sup+epsilon))[2:], 'o', markersize=1, color=\"blue\", label=\"normalized diff FFT\")\n",

- " #plt.plot(x_support, norm_y_est_num, 'o', markersize=1, color=\"red\", label=\"normalized diff NUM\")\n",

- " #plt.plot(x_support[2:], norm_y_est_fft[2:], 'o', markersize=1, color=\"blue\", label=\"normalized diff FFT\")\n",

- " #plt.plot(x_support, norm_y_sol_sup, 'o', markersize=1, color=\"green\", label=\"normalized diff NUM\")\n",

- " plt.plot(x_support, (norm_y_est_num -norm_y_sol_sup)/(norm_y_sol_sup+epsilon), 'o', markersize=1, color=\"red\", label=\"normalized diff NUM\")\n",

- " plt.ylabel(\"(calculation - analytical_sol)/(analytical_sol + epsilon)\", fontsize=7)\n",

- " plt.xlabel(\"Time\", fontsize=7)\n",

- " plt.title(\"Normalized Difference Estimated and True Solution\", fontsize=12)\n",

- " plt.axvline(x=eff_support[1], label= \"effective support end\")\n",

- " plt.legend(prop={\"size\":10})\n",

- " #plt.yscale(\"log\")\n",

- " print(\"Absolute difference at end of effective support: \" +str(get_P(eff_support[1] , accuracy_fft[t])- get_P(eff_support[1], dist_sol)))\n",

- " print(\"Absolute difference at end plus epsilon of effective support: \" +str(get_P(eff_support[1] + epsilon, accuracy_fft[t])- get_P(eff_support[1]+epsilon, dist_sol)))\n",

- " print(\"Absolute difference at start of effective support: \" +str(get_P(eff_support[0], accuracy_fft[t])- get_P(eff_support[0], dist_sol)))\n",

- " \n",

- "print_relative_diff_grids(plot_support_closeup=True)\n"

- " '''\n",

- " get the location of the start and end points of the linearly spaced output grid \n",

- " '''\n",

- " est_new_peak_pos = dist_1.peak_pos + dist_2.peak_pos ##as we have assumed inverse time\n",

- " est_joint_fwhm = (dist_1.fwhm + dist_2.fwhm)\n",

- " center_width = 3*est_joint_fwhm\n",

- " print(\"The bounds for the linear grid should be:\" + str(est_new_peak_pos - center_width) + \"to \" + str(est_new_peak_pos + center_width))\n",

- " plt.axvline(x=end_conv_linear, color=\"green\", label= \"linear grid end\")\n",

- " plt.legend(prop={\"size\":8})\n",

- "\n",

- " x_linear = np.linspace(end_conv_linear-0.1, end_conv_linear+0.1, num=1000, endpoint=False) #create an x vector of points that should be sampled\n",

- " norm_y_est_num = get_P(x_linear, accuracy_num[t]) \n",

- " plt.plot(x_linear[2:], ((norm_y_est_fft -norm_y_sol_sup)/(norm_y_sol_sup+epsilon))[2:], 'o', markersize=1, color=\"blue\", label=\"normalized diff FFT\")\n",

- " plt.plot(x_linear, (norm_y_est_num -norm_y_sol_sup)/(norm_y_sol_sup+epsilon), 'o', markersize=1, color=\"red\", label=\"normalized diff NUM\")\n",

- " plt.ylabel(\"(calculation - analytical_sol)/(analytical_sol + epsilon)\", fontsize=7)\n",

- " plt.xlabel(\"Time\", fontsize=7)\n",

- " plt.title(\"Normalized Difference Estimated and True Solution\", fontsize=12)\n",

- " plt.axvline(x=end_conv_linear, color=\"green\", label= \"linear grid end\")\n",

- " if (end_conv_linear+0.1 > eff_support[1]):\n",

- " plt.axvline(x=eff_support[1], color=\"blue\", label= \"effective support end\")\n",

- " plt.legend(prop={\"size\":10})\n",

- " print(\"FFT grid points (output):\" +str(time_fft[t][1]))\n",

- " print(\"NUM grid points (output):\" +str(time_num[t][1]))\n",

- " print(\"FFT grid points, after adjusting grid (output):\" +str(len(accuracy_fft[t].y)))\n",

- " print(\"NUM grid points, after adjusting grid (output):\" +str(len(accuracy_num[t].y)))\n",

- " x_support = np.linspace(eff_support[1]-0.1, eff_support[1]+0.1, num=1000, endpoint=False) #create an x vector of points that should be sampled\n",

- " plt.plot(x_support[2:], ((norm_y_est_fft -norm_y_sol_sup)/(norm_y_sol_sup+epsilon))[2:], 'o', markersize=1, color=\"blue\", label=\"normalized diff FFT\")\n",

- " plt.plot(x_support, (norm_y_est_num -norm_y_sol_sup)/(norm_y_sol_sup+epsilon), 'o', markersize=1, color=\"red\", label=\"normalized diff NUM\")\n",

- " plt.ylabel(\"(calculation - analytical_sol)/(analytical_sol + epsilon)\", fontsize=7)\n",

- " plt.xlabel(\"Time\", fontsize=7)\n",

- " plt.title(\"Normalized Difference Estimated and True Solution\", fontsize=12)\n",

- " plt.axvline(x=eff_support[1], label= \"effective support end\")\n",

- " plt.legend(prop={\"size\":10})\n",

- "plot_effects_adjust_grid() "

- "print(MAX_BRANCH_LENGTH) \n",

- " grid_left = sigma * (1 - np.linspace(1, 0.0, n_grid_points//3)**2.0)\n",

- " grid_zero = grid_left[1]*np.logspace(-20,0,6)[:5]\n",

- " grid_zero2 = grid_left[1]*np.linspace(0,1,10)[1:-1]\n",

- " grid_right = sigma + (3*sigma*(np.linspace(0, 1, n_grid_points//3)**2))\n",

- " far_grid = grid_right.max() + MAX_BRANCH_LENGTH*np.linspace(0, 1, n_grid_points//3)**2\n",

- " grid = np.concatenate((grid_zero,grid_zero2, grid_left,grid_right[1:],far_grid[1:]))\n",

- " grid.sort() # just for safety\n",

- "y_1_bl = gamma.pdf(x_1, alpha_1, scale=1/beta)\n",

- "y_1_bl[y_1_bl == 0] = TINY_NUMBER ##remove -inf to prevent numerical errors\n",

- "y_2_bl = gamma.pdf(x_2, alpha_2, scale=1/beta)\n",

- "y_2_bl[y_2_bl == 0] = TINY_NUMBER ##remove -inf to prevent numerical errors\n",

- "plt.plot(x_2, gamma.pdf(x_2, alpha_2, scale=1/beta), color=\"blue\")\n",

- "plt.plot(x_2, np.exp(-dist_2_bl.y), color=\"green\") ##just to check that distribution object is initialized accurately \n",

- "plt.ylabel(\"PDF Gamma(\"+ str(alpha_2) +\",\"+str(beta)+ \")\")\n",

- "plt.title(\"Distribution 2\")\n",

- "y_sol = gamma.pdf(x, alpha_sol, scale=1/beta)\n",

- "y_sol[y_sol == 0] = TINY_NUMBER ##remove -inf to prevent numerical errors\n",

- "dist_sol = Distribution(x, y_sol, kind='linear', is_log=False)\n"

- "grid_size = [100, 200, 300, 500, 1000] #more desired grid_sizes, true grid size is given as output when the accuracy_and_time function is called\n",

- "time_fft = np.empty((len(grid_size),2))\n",

- "time_num = np.empty((len(grid_size),2))\n",

- " #print(t)\n",

- " time_fft_full = accuracy_and_time(dist_1_bl, dist_2_bl, \"fft\", grid_size[t])\n",

- " time_num_full = accuracy_and_time(dist_1_bl, dist_2_bl, \"num\", grid_size[t])\n",

- " \n",

- tree_events_tt = sorted([(n.time_before_present, n.name, int(n.bad_branch)) for n in tt.tree.find_clades()], key=lambda x:-x[0])

- bad_branches = np.float_(old_times_and_names.iloc[:,1]) + 2*np.float_(new_times_and_names.iloc[:,1])

- df = pd.DataFrame(dict(time=np.float_(new_times_and_names.iloc[:,0]),

- difference=(np.float_(new_times_and_names.iloc[:,0])-np.float_(old_times_and_names.iloc[:,0])),

- bad_branch=bad_branches), index=new_times_and_names.index)

- return np.all(new_times_and_names.iloc[:,0]==old_times_and_names.iloc[:,0]), df

- large_differences = df[abs(df.difference) > abs(np.mean(df.difference)) + 2*np.std(df.difference)]

- plt.xlabel("nodes ranging from root at 0 to most recent")

- plt.ylabel("difference time_before_present fft - numerical")

- plt.plot(dist_node_num.x, dist_node_num.prob_relative(dist_node_num.x), marker='o', linestyle='', ms=2, label='numerical')

- if xlimits!='empty':

- plt.xlim((x_min-0.0001, x_max+0.0001))

- plt.title(node_num.name + "_" + dist_name)

- ebola=True

- seq_kwargs = {"marginal_sequences":True,

- "branch_length_mode": 'input',

- "sample_from_profile":"root",

- "reconstruct_tip_states":False}

- tt_kwargs = {'clock_rate':clock_rate,

- 'time_marginal':'assign'}

- coal_kwargs ={'Tc':10000,

- 'time_marginal':'assign'}

-

- dates = parse_dates(base_name+'.metadata.csv')

- tt = TreeTime(gtr='Jukes-Cantor', tree = base_name+'.nwk', use_fft=False,

- aln = base_name+'.fasta', verbose = 1, dates = dates, precision=3, debug=True)

- tt_fft = TreeTime(gtr='Jukes-Cantor', tree = base_name+'.nwk', use_fft=True, precision_fft=200,

- aln = base_name+'.fasta', verbose = 1, dates = dates, precision=3, debug=True)

- tree._set_branch_length_mode(seq_kwargs["branch_length_mode"])

- tree.infer_ancestral_sequences(infer_gtr=False, marginal=seq_kwargs["marginal_sequences"])

- tree.clock_filter(reroot='least-squares', n_iqd=1, plot=False, fixed_clock_rate=tt_kwargs["clock_rate"])

- tree.reroot(root='least-squares', clock_rate=tt_kwargs["clock_rate"])

- tree.make_time_tree(clock_rate=tt_kwargs["clock_rate"], time_marginal=tt_kwargs["time_marginal"])

- #compare marginal_pos_LH

- compare_dist(test_node, test_node_fft, test_node.marginal_pos_LH, test_node_fft.marginal_pos_LH, 'pos_LH')

- #compare msg_from_parent

- compare_dist(test_node, test_node_fft, test_node.msg_from_parent, test_node_fft.msg_from_parent, 'msg_from_parent')

- #compare subtree_distribution

- compare_dist(test_node, test_node_fft, test_node.subtree_distribution, test_node_fft.subtree_distribution, 'subtree_dist')

- #compare branch_length_interpolator

- compare_dist(test_node, test_node_fft, test_node.branch_length_interpolator, test_node_fft.branch_length_interpolator, 'branch length dist', xlimits='empty')

- #compare marginal_pos_Lx

- xlimits = (test_node.up.time_before_present-0.01,test_node.up.time_before_present+0.01)

- compare_dist(test_node, test_node_fft, test_node.marginal_pos_Lx, test_node_fft.marginal_pos_Lx, 'marginal_pos_Lx', xlimits=xlimits)

-

- #compare marginal_pos_LH of parent

- print("The parent of numerical is: " + str(test_node.up.name) + " the parent of fft is: " + str(test_node_fft.up.name))

- compare_dist(test_node.up, test_node_fft.up, test_node.up.marginal_pos_LH, test_node_fft.up.marginal_pos_LH, 'pos_LH parent')

-

- #compare marginal_pos_LH of children

- print("Now looking at children of node of interest")

- print("Time difference is " + str(test_node.time_before_present - test_node_fft.time_before_present ))

- print("Time of numerical:" + str(test_node.time_before_present) + " time of fft:"+ str(test_node_fft.time_before_present))

- print("delta")

- compare_dist(test_node, test_node_fft, test_node.marginal_pos_LH, test_node_fft.marginal_pos_LH, 'pos_LH')

- tree.add_coalescent_model(coal_kwargs ["Tc"])

- tree.make_time_tree(clock_rate=tt_kwargs ["clock_rate"], time_marginal=coal_kwargs ["time_marginal"])

- #plot_differences(output_comparison[1])

- msgs_to_multiply_tt_fft = [child.marginal_pos_Lx for child in test_node_fft.clades

- if child.marginal_pos_Lx is not None]

- msgs_to_multiply_tt = [child.marginal_pos_Lx for child in test_node.clades

- if child.marginal_pos_Lx is not None]

- compare_dist(test_node, test_node_fft, subtree_distribution_tt, subtree_distribution_tt_fft, 'multiplied Lx of children')

- # compare_multiply_functions(nodes_to_look_at) \ No newline at end of file

- df_list[i]= df_list[i].reindex(index=df_list[0].index)

- tree._set_branch_length_mode(kwargs["branch_length_mode"])

- tree.infer_ancestral_sequences(infer_gtr=False, marginal=kwargs["marginal_sequences"])

- tree.clock_filter(reroot='least-squares', n_iqd=1, plot=False, fixed_clock_rate=kwargs["clock_rate"])

- tree.reroot(root='least-squares', clock_rate=kwargs["clock_rate"])

- tree.make_time_tree(clock_rate=kwargs["clock_rate"], time_marginal=kwargs["time_marginal"], divide=False)

- par2.axis["right"] = new_fixed_axis(loc="right", axes=par2, offset=(offset, 0))

- par2.axis["right"].toggle(all=True)

- par1.axis["right"].toggle(all=True)

- host.set_xlabel("Branch Grid Size")

- par1.set_ylabel("Time [sec]")

- par2.set_ylabel("Memory Consumption [MB]")

- if len(accuracy)==(len(branch_grid_size)-1):

- p1, = host.plot(branch_grid_size[1:], accuracy, label=title)

- p2, = par1.plot(branch_grid_size, np.array(time_branch), label='time')

- p3, = par2.plot(branch_grid_size, np.array(memory_branch), label='memory consumption')

- p1, = host.plot(branch_grid_size, accuracy, label=title)

- p2, = par1.plot(branch_grid_size, np.array(time_branch), label='time')

- p3, = par2.plot(branch_grid_size, np.array(memory_branch), label='memory consumption')

- host.axis["left"].label.set_color(p1.get_color())

- par1.axis["right"].label.set_color(p2.get_color())

- par2.axis["right"].label.set_color(p3.get_color())

- plt.title(title + " - FFT Grid Size=" + str(prec_fft))

- plt.savefig("./results/Memory_Time_" + str(title)+ ".png")

-

- ebola=True

- dates = parse_dates(base_name+'.metadata.csv')

- tt_numerical = TreeTime(gtr='Jukes-Cantor', tree = base_name+'.nwk', use_fft=False,

- aln = base_name+'.fasta', verbose = 1, dates = dates, precision=3, debug=True)

- kwargs = {"marginal_sequences":True,

- "branch_length_mode": 'input',

- "sample_from_profile":"root",

- "reconstruct_tip_states":False,

- "max_iter":1,

- 'clock_rate':clock_rate,

- 'time_marginal':'assign'}

- print("Time to make tree using ultra fine numerical grid: " + str(tree_make_time))

- numerical_LH= tt_numerical.tree.positional_LH

- print("likelihood of tree under numerical:" + str(numerical_LH))

- '''

- '''

- time_fft = [] #time needed to make fft

- pp = PdfPages('./results/DifferencesComparison.pdf') #output the difference plots to visually assess if there is an issue

- pre_divergence_times_fft_list =[]

- tt_fft = TreeTime(gtr='Jukes-Cantor', tree = base_name+'.nwk', precision_fft=prec, use_fft=True,

- aln = base_name+'.fasta', verbose = 1, dates = dates, precision=3, precision_branch=b, debug=True)

- print("Time to make tree using fft grid of size " + str(prec) +" :"+str(tree_make_time_fft))

- if len([t[1] for t in tree_events_fft if np.isnan(t[0])])>0:

- print("error in tree! a clade has no time before present!")

- fig = plot_differences(output_comparison[1], title="Differences branch grid " + str(b) + " fft grid " + str(prec))

- if len([t for t in sorted_treelist[i][0] if np.isnan(t)])>0:

- divergence_time_diff.append(np.mean(np.abs(divergence_times[i] - divergence_times[i-1])/(precision_fft[i] - precision_fft[i-1])))

- #for each possible branch length grid size plot the divergence time differences for different fft grid sizes

- plt.plot(precision_fft[1:len(divergence_times_fft_list[i][0])], divergence_time_diff_list[i][0], 'o-', label=branch_grid_size[i])

- plt.legend(title="branch grid size")

- plt.ylabel("log(average rate of change)")

- plt.xlabel("FFT grid size")

- plt.title("Change in Infered Divergence Times - Average Rate of Change")

- plt.savefig("./results/Optimal_Grid_Size.png")

- time_fft, divergence_times_fft_list, divergence_time_diff_list = find_optimal_branch_fft_gridsizes(precision_fft, branch_grid_size)

- pre_divergence_times_fft_list =[]

- tt_fft = TreeTime(gtr='Jukes-Cantor', tree = base_name+'.nwk', precision_fft=prec_fft, use_fft=True,

- aln = base_name+'.fasta', verbose = 1, dates = dates, precision=3, precision_branch = b, debug=True)

- print("Time to make tree using fft grid of size " + str(150) +" :"+str(tree_make_time_fft))

- if len([t[1] for t in tree_events_fft if np.isnan(t[0])])>0:

- print("error in tree! a clade has no time before present!")

- plot_differences(output_comparison[1], title="Differences branch grid " + str(b) + " fft grid " + str(prec_fft))

- print("likelihood of tree under fft:" + str(tt_fft.tree.positional_LH))

-

- if len([t for t in sorted_treelist[i][0] if np.isnan(t)])>0:

- divergence_time_diff.append(np.nanmean(np.abs(divergence_times[i] - divergence_times[i-1])/(branch_grid_size[i] - branch_grid_size[i-1])))

-

- memory_branch = np.array(memory_branch)/(10**6)

- accuracy_rate = np.array(divergence_time_diff)/1e-8

- plot_accuracy_memory_time(branch_grid_size, accuracy_rate, time_branch, memory_branch, "Average Rate of Change Inferred Time [1e-8]", "average_rate_of_change")

- plot_accuracy_memory_time(branch_grid_size, np.array(fft_LH), time_branch, memory_branch, "tree log(LH) ", "log(LH)")

-

-

-

- return test_node.branch_length_interpolator.prob(t)*np.exp(tt.merger_model.cost(test_node.time_before_present, t))

-

- ebola=False

- dates = parse_dates(base_name+'.metadata.csv')

-

- tt_old= TreeTime(gtr='Jukes-Cantor', tree = base_name+'.nwk',

- aln = base_name+'.fasta', verbose = 1, dates = dates, precision=3, debug=True)

- tt_smooth= TreeTime(gtr='Jukes-Cantor', tree = base_name+'.nwk',

- aln = base_name+'.fasta', verbose = 1, dates = dates, precision=3, debug=True)

- #tt._ml_t_marginal(assign_dates=True)

- print("Time for smooth nbranches (using posterior dist):" + str(time.process_time()-start))

- print("Time for discrete nbranches:" + str(time.process_time()-start))

- node_pattern= 'V517'

- plt.plot(tt_old.merger_model.nbranches.x/tt_old.date2dist.clock_rate, tt_old.merger_model.nbranches.y, label="old nbranches function")

- plt.plot(tt_smooth.merger_model.nbranches.x/tt_smooth.date2dist.clock_rate, tt_smooth.merger_model.nbranches.y, label="new smooth nbranches function")

- plt.xlabel("time before present")

- plt.ylabel("nbranches")

- plt.xlim((0,40))

-

- t = np.linspace(0,4*fixed_clock_rate,1000)

- test_nodes =[test_node.up for test_node in test_nodes]

-

- fig, axs = plt.subplots(1,2, sharey=True, figsize=(12,8))

- Phylo.draw(tt_old.tree, label_func=lambda x:"", axes=axs[0])

- axs[0].set_title("old time tree")

- Phylo.draw(tt_smooth.tree, label_func=lambda x:"", axes=axs[1])

- axs[1].set_title("smooth time tree")

- next=0

- fig, axs = plt.subplots(2,2, sharey=True, figsize=(12,8))

- fig.suptitle("Effect of smooth nbranches on LH distribution")

- if next==0:

- i, j=0, 0

- elif next==1:

- i, j=1, 0

- elif next==2:

- i, j=0, 1

- elif next==3:

- i, j=1, 1

- t = np.linspace(test_nodes[0].time_before_present-0.005,test_nodes[0].time_before_present+0.003,2000)

- axs[i,j].plot(t, test_nodes[0].marginal_pos_LH.prob_relative(t), label='old', ls='-')

- axs[i,j].plot(t, test_nodes[1].marginal_pos_LH.prob_relative(t), label='smooth', ls='-')

- axs[i,j].set_xlabel("time before present")

- axs[i,j].set_ylabel("marginal LH")

- axs[i,j].legend()

- if next==4:

- next=0

- fig, axs = plt.subplots(2,2, sharey=True, figsize=(12,8))

- fig.suptitle("Effect of smooth nbranches on LH distribution")

- test_nodes =[test_node.up for test_node in test_nodes]

-

- print("finished")

-

- print("testing short imports")

- print("testing assign gamma")

- fasta = root_dir + "/treetime_examples/data/h3n2_na/h3n2_na_20.fasta"

- nwk = root_dir + "/treetime_examples/data/h3n2_na/h3n2_na_20.nwk"

- dates = parse_dates(root_dir + "/treetime_examples/data/h3n2_na/h3n2_na_20.metadata.csv")

- seq_kwargs = {"marginal_sequences":True,

- "branch_length_mode": 'input',

- "reconstruct_tip_states":False}

- tt_kwargs = {'clock_rate': 0.0001,

- 'time_marginal':'assign'}

- myTree = TreeTime(gtr='Jukes-Cantor', tree = nwk, use_fft=False,

- aln = fasta, verbose = 1, dates = dates, precision=3, debug=True, rng_seed=1234)

- custom_gtr = root_dir + "/test_sequence_evolution_model.txt"

- assert (gtr.Pi.sum() - 1.0)**2<1e-14

- assert np.allclose(gtr.W, np.array([[0, 0.7003, 3.0669, 0.2651, 0.9742],

- [0.7003, 0, 0.3354, 3.399, 0.999],

- [3.0669, 0.3354, 0, 0.4258, 0.9892],

- [0.2651, 3.399, 0.4258, 0, 0.9848],

- [0.9742, 0.999, 0.9892, 0.9848, 0]]), atol=1e-4)

- print('testing GTR, model:',model)

- assert (myGTR.Pi.sum() - 1.0)**2<1e-14

- assert np.abs(myGTR.v.sum()) > 1e-10 # **and** v is not zero

- tiny_tree = Phylo.read(StringIO("((A:0.60100000009,B:0.3010000009):0.1,C:0.2):0.001;"), 'newick')

- "A": "?",

- "B": "North America",

- "C": "West Asia",

- missing_data="?",

- assert "A" in letter_to_state

- assert "B" in letter_to_state

- assert letter_to_state["C"] == "?"

- print('ancestral reconstruction' + ("marginal" if marginal else "joint"))

- assert t.data.compressed_to_full_sequence(t.tree.root.cseq, as_string=True) == 'ATGAATCCAAATCAAAAGATAATAACGATTGGCTCTGTTTCTCTCACCATTTCCACAATATGCTTCTTCATGCAAATTGCCATCTTGATAACTACTGTAACATTGCATTTCAAGCAATATGAATTCAACTCCCCCCCAAACAACCAAGTGATGCTGTGTGAACCAACAATAATAGAAAGAAACATAACAGAGATAGTGTATCTGACCAACACCACCATAGAGAAGGAAATATGCCCCAAACCAGCAGAATACAGAAATTGGTCAAAACCGCAATGTGGCATTACAGGATTTGCACCTTTCTCTAAGGACAATTCGATTAGGCTTTCCGCTGGTGGGGACATCTGGGTGACAAGAGAACCTTATGTGTCATGCGATCCTGACAAGTGTTATCAATTTGCCCTTGGACAGGGAACAACACTAAACAACGTGCATTCAAATAACACAGTACGTGATAGGACCCCTTATCGGACTCTATTGATGAATGAGTTGGGTGTTCCTTTTCATCTGGGGACCAAGCAAGTGTGCATAGCATGGTCCAGCTCAAGTTGTCACGATGGAAAAGCATGGCTGCATGTTTGTATAACGGGGGATGATAAAAATGCAACTGCTAGCTTCATTTACAATGGGAGGCTTGTAGATAGTGTTGTTTCATGGTCCAAAGAAATTCTCAGGACCCAGGAGTCAGAATGCGTTTGTATCAATGGAACTTGTACAGTAGTAATGACTGATGGAAGTGCTTCAGGAAAAGCTGATACTAAAATACTATTCATTGAGGAGGGGAAAATCGTTCATACTAGCACATTGTCAGGAAGTGCTCAGCATGTCGAAGAGTGCTCTTGCTATCCTCGATATCCTGGTGTCAGATGTGTCTGCAGAGACAACTGGAAAGGCTCCAATCGGCCCATCGTAGATATAAACATAAAGGATCATAGCATTGTTTCCAGTTATGTGTGTTCAGGACTTGTTGGAGACACACCCAGAAAAAACGACAGCTCCAGCAGTAGCCATTGTTTGGATCCTAACAATGAAGAAGGTGGTCATGGAGTGAAAGGCTGGGCCTTTGATGATGGAAATGACGTGTGGATGGGAAGAACAATCAACGAGACGTCACGCTTAGGGTATGAAACCTTCAAAGTCATTGAAGGCTGGTCCAACCCTAAGTCCAAATTGCAGATAAATAGGCAAGTCATAGTTGACAGAGGTGATAGGTCCGGTTATTCTGGTATTTTCTCTGTTGAAGGCAAAAGCTGCATCAATCGGTGCTTTTATGTGGAGTTGATTAGGGGAAGAAAAGAGGAAACTGAAGTCTTGTGGACCTCAAACAGTATTGTTGTGTTTTGTGGCACCTCAGGTACATATGGAACAGGCTCATGGCCTGATGGGGCGGACCTCAATCTCATGCCTATA'

- from Bio import Phylo,AlignIO

- tiny_tree = Phylo.read(StringIO("((A:0.60100000009,B:0.3010000009):0.1,C:0.2):0.001;"), 'newick')

- tiny_aln = AlignIO.read(StringIO(">A\nAAAAAAAAAAAAAAAACCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTT\n"

- ">B\nAAAACCCCGGGGTTTTAAAACCCCGGGGTTTTAAAACCCCGGGGTTTTAAAACCCCGGGGTTTT\n"

- ">C\nACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT\n"), 'fasta')

- mygtr = GTR.custom(alphabet = np.array(['A', 'C', 'G', 'T']), pi = np.array([0.9, 0.06, 0.02, 0.02]), W=np.ones((4,4)))

- lhsum = np.exp(t.sequence_LH(pos=np.arange(4**3))).sum()

- print (lhsum)

- assert(np.abs(lhsum-1.0)<1e-6)

- tiny_tree = Phylo.read(StringIO("((A:.060,B:.01200)C:.020,D:.0050)E:.004;"), 'newick')

- mygtr = GTR.custom(alphabet = np.array(['A', 'C', 'G', 'T']),

- pi = np.array([0.15, 0.95, 0.05, 0.3]),

- W=np.ones((4,4)))

-

- print ("Root sequence: " + ''.join(tree.root.ref_seq.astype('U')))

- p = mygtr.evolve( seq_utils.seq2prof(node.up.ref_seq, mygtr.profile_map), t)

- p=(p.T/p.sum(axis=1)).T

- node.ref_mutations = [(anc, pos, der) for pos, (anc, der) in

- enumerate(zip(node.up.ref_seq, node.ref_seq)) if anc!=der]

- print (node.name, len(node.ref_mutations), node.ref_mutations)

- alnstr = ""

- alnstr += ">" + leaf.name + "\n" + ''.join(leaf.ref_seq.astype('U')) + '\n'

- print (alnstr)

- print("%s: True sequence does not equal inferred sequence. parent %s"%(node.name, node.up.name))

- print("%s: True sequence equals inferred sequence. parent %s"%(node.name, node.up.name))

- assert diff_count/seq_len<2*(1.0*mut_count/seq_len)**2

- LH_p = (myTree.tree.sequence_LH)

- print ("Difference between reference and inferred LH:", (LH - LH_p).sum())

- assert ((LH - LH_p).sum())<1e-9

- mygtr = GTR.custom(alphabet = np.array(['A', 'C', 'G', 'T']), pi = np.array([0.91, 0.05, 0.02, 0.02]), W=np.ones((4,4)))

- tiny_tree = Phylo.read(StringIO("((A:.0060,B:.30)C:.030,D:.020)E:.004;"), 'newick')

- #terminal node sequences (single nuc)

- A_seq = ''.join(np.repeat(A_char,16))

- B_seq = ''.join(np.repeat(B_char,16))

- D_seq = ''.join(np.repeat(D_char,16))

- tiny_aln = AlignIO.read(StringIO(">A\n" + A_seq + "\n"

- ">B\n" + B_seq + "\n"

- ">D\n" + D_seq + "\n"

- ">C\nAAAACCCCGGGGTTTT\n"

- ">E\nACGTACGTACGTACGT\n"), 'fasta')

- myTree = TreeAnc(gtr=mygtr, tree = tiny_tree,

- aln =tiny_aln, verbose = 4, rng_seed=1234)

- tiny_aln_1 = AlignIO.read(StringIO(">A\n"+A_char+"\n"

- ">B\n"+B_char+"\n"

- ">D\n"+D_char+"\n"), 'fasta')

- myTree_1 = TreeAnc(gtr=mygtr, tree = tiny_tree,

- aln=tiny_aln_1, verbose = 4, rng_seed=1234)

- real = real_lh()

- print(abs(ref.max() - real) )

- assert(abs(ref.max() - real) < 1e-10)

- fname = d / "no-header.vcf"

-class TestMalformedVcf:

- read_vcf(create_vcf(tmp_path, """\

- """))

- read_vcf(create_vcf(tmp_path, """\

- """))

-

- read_vcf(create_vcf(tmp_path, """\

- """))

- read_vcf(create_vcf(tmp_path, """\

- """))

- read_vcf(create_vcf(tmp_path, """\

- """))

- read_vcf(create_vcf(tmp_path, """\

- """))

- return idx-1

-def write_data(tmp_path_factory, sample_names, data_lines, reference_seq, meta_lines=None, reference_name='reference_name'):

- vcf = "\n".join(meta_lines)+"\n"

- vcf += "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t" + "\t".join(sample_names) + "\n"

- vcf += "\t".join(line) + "\n"

- reference=dedent(f"""\

- dir = tmp_path_factory.mktemp("data")

- vcf_fn = dir / "snps.vcf"

- ref_fn = dir / "reference.fasta"

- @pytest.fixture(scope="class")

- sample_names = ["sample_A", "sample_B", "sample_C"]

- ["1", "5", ".", "T", "C", ".", ".", ".", "GT", "1", "1", "1"],

- ["1", "7", ".", "T", "G", ".", ".", ".", "GT", "1", "1", "0"],

- ["1", "14", ".", "C", "T", ".", ".", ".", "GT", "1", "1", "0"],

- ["1", "18", ".", "C", "T", ".", ".", ".", "GT", "0", "0", "1"],

- ["1", "28", ".", "A", "N", ".", ".", ".", "GT", "0", "0", "1"],

- ["1", "29", ".", "A", "N", ".", ".", ".", "GT", "0", "0", "1"],

- ["1", "30", ".", "A", "N", ".", ".", ".", "GT", "0", "0", "1"],

- ["1", "33", ".", "A", "C,G", ".", ".", ".", "GT", "1", "2", "0"],

- ["1", "39", ".", "C", "T", ".", ".", ".", "GT", "1", "0", "0"],

- ["1", "42", ".", "G", "A", ".", ".", ".", "GT", "0", "1", "0"],

- ["1", "43", ".", "A", "T", ".", ".", ".", "GT", "1", "1", "0"],

- "filenames": filenames,

- "positions": positions,

- "data_lines": data_lines,

- "sample_names": sample_names,

- "reference_seq": reference_seq,

- assert(set(data['sample_names'])==set(vcf_data['sequences'].keys()))

- assert(set(data['sample_names'])==set(vcf_data['insertions'].keys()))

- return ",".join([f"{pos+1}{alt}" for pos,alt in vcf_data['sequences'][sample].items()])

- assert("5C,7G,14T,33C,39T,43T"==summarise('sample_A'))

- assert("5C,7G,14T,33G,42A,43T"==summarise('sample_B'))

- assert("5C,18T,28N,29N,30N"==summarise('sample_C'))

- assert(ins=={})

- assert(data['reference_seq']==vcf_data['reference'])

- assert(data['positions'] == [pos+1 for pos in vcf_data['positions']])

- @pytest.fixture(scope="class")

- reference="ATCGA"

- sample_names = ["sample_A", "sample_B", "sample_C", "sample_D"]

- ["1", "2", ".", "T", "A", ".", ".", ".", "GT", ".", "1", "0", "0"], # sample A has T2N

- ["1", "4", ".", "GA", "C", ".", ".", ".", "GT", ".", "1", "0", "0"], # sample A has GA -> NN

- ["1", "3", ".", "C", "*,G", ".", ".", ".", "GT", "1", "2", "0", "0"], # sample A has C->N

- ["1", "3", ".", "C", "T,.", ".", ".", ".", "GT", "0", "0", "1", "2"], # sample D has C->N

- ["1", "4", ".", "GA", ".,*", ".", ".", ".", "GT", "0", "0", "1", "2"], # both samples C & D have G->N and A->N

- return {"filenames": filenames}

- assert(vcf_data['sequences']['sample_A'][zero_based(2)]=='N')

- assert(vcf_data['sequences']['sample_A'][zero_based(4)]=='N')

- assert(vcf_data['sequences']['sample_A'][zero_based(5)]=='N')

- assert(vcf_data['sequences']['sample_B'][zero_based(4)]=='C')

- assert(vcf_data['sequences']['sample_B'][zero_based(5)]=='-')

- assert(vcf_data['sequences']['sample_A'][zero_based(3)]=='N')

- assert(vcf_data['sequences']['sample_B'][zero_based(3)]=='G')

- assert(vcf_data['sequences']['sample_D'][zero_based(4)]=='N')

- assert(vcf_data['sequences']['sample_D'][zero_based(5)]=='N')

- assert(vcf_data['sequences']['sample_C'][zero_based(3)]=='T')

- assert(vcf_data['sequences']['sample_D'][zero_based(3)]=='N')

-

- assert(vcf_data['sequences']['sample_C'][zero_based(4)]=='N')

- assert(vcf_data['sequences']['sample_C'][zero_based(5)]=='N')

- read_vcf(create_vcf(tmp_path, """\

- """))

- read_vcf(create_vcf(tmp_path, """\

- """))

- lines = ["##fileformat=VCFv4.2", # Actual version may differ, but we don't parse this so it doesn't matter (yet)

- "##contig=<ID=1,length=50>",

- "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">",

- "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t" + "\t".join(input['samples']),

- lines.append(f"{d['chrom']}\t{d['pos']}\t.\t{d['ref']}\t{d['alt']}\t.\t.\t.\tGT\t" + "\t".join([str(x) for x in d['values']]))

- return create_vcf(tmp_path, "\n".join(lines))

- assert(sample_name in vcf['sequences'])

- assert(len(snp)==2)

- assert(vcf['sequences'][sample_name][zero_based(int(snp[0]))] == snp[1])

- assert(len(vcf['sequences'][sample_name].keys()) == len(sample_data))

- assert(sample_name in vcf['insertions'])

- assert(vcf['insertions'][sample_name][zero_based(int(ins[0]))] == ins[1:])

- assert(len(vcf['insertions'][sample_name].keys()) == len(sample_data))

- {'chrom': 20, 'pos': 2, 'ref': 'TC', 'alt': 'T,TCA', 'values': [0, 1, 2]}

- 'expect_insertions': {'C': ['3CA']} # see comment above. Actually only a single base "A" insertion

- 'data': [

- {'chrom': 20, 'pos': 2, 'ref': 'TC', 'alt': 'TG,T', 'values': [1,2]}

- ],

- 'expect_insertions': {}

- 'data': [

- {'chrom': 20, 'pos': 2, 'ref': 'TCG', 'alt': 'TG,T,TCAG', 'values': [1,2,3]}

- ],

- 'expect_insertions': {'C': ['4AG']} # See comment above

- 'data': [

- {'chrom': 20, 'pos': 3, 'ref': 'C', 'alt': 'T', 'values': [1]}

- ],

- 'expect_insertions': {}

- 'data': [

- {'chrom': 20, 'pos': 3, 'ref': 'C', 'alt': 'CTAG', 'values': [1]}

- ],

- 'expect_insertions': {'A': ['3CTAG']} # See comment above

- 'data': [

- {'chrom': 20, 'pos': 2, 'ref': 'TCG', 'alt': 'T', 'values': [1]}

- ],

- 'expect_insertions': {}

- setattr(TestVcfSpecExamples, f"test_{spec_version}_example_{example_key}", TestVcfSpecExamples.add_test(example_data))

- @pytest.fixture(scope="class")

- reference="ATCGA"

- sample_names = ["sample_A", "sample_B"]

- ["1", "2", ".", "T", "GA", ".", ".", ".", "GT", "1", "0"], # sample A has both a C->G mutation _and_ a subsequent "A" insertion

- ["1", "3", ".", "CGA", "NTTTA,CCGT", ".", ".", ".", "GT", "1", "2"], # complex! Both samples have multiple mutations + an insertion

- return {"filenames": filenames}

- assert(vcf_data['sequences']['sample_A'][zero_based(2)]=='G')

- assert(vcf_data['insertions']['sample_A'][zero_based(2)]=='GA') # see comment above re: insertion encoding

- assert(vcf_data['sequences']['sample_A'][zero_based(3)]=='N')

- assert(vcf_data['sequences']['sample_A'][zero_based(4)]=='T')

- assert(vcf_data['sequences']['sample_A'][zero_based(5)]=='T')

- assert(vcf_data['insertions']['sample_A'][zero_based(5)]=='TTA') # see comment above re: insertion encoding

- assert(vcf_data['sequences']['sample_B'][zero_based(4)]=='C')

- assert(vcf_data['sequences']['sample_B'][zero_based(5)]=='G')

- assert(vcf_data['insertions']['sample_B'][zero_based(5)]=='GT') # see comment above re: insertion encoding

-

- data = read_vcf(create_vcf(tmp_path, """\

- """))

- assert(data['metadata']['chrom']=='1')

- assert(data['metadata']['ploidy']==1)

- assert(data['metadata']['meta_lines']==['##fileformat=VCFv4.3'])

- data = read_vcf(create_vcf(tmp_path, """\

- """))

- assert(data['metadata']['chrom']=='foo')

- assert(data['metadata']['ploidy']==2)

- assert(data['metadata']['meta_lines']==['##fileformat=VCFv4.3', '##contig=<ID=1,length=50>', '##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">'])

-

-

- vcf_filename_a = str(dir / "a.vcf")

- vcf_filename_b = str(dir / "b.vcf")

- print(dedent(f"""\

- """), file=fh)

- vcf_lines = meta_lines[:] + \

- ["#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t" + "\t".join(sample_names)] + \

- ["\t".join(data) for data in data_lines]

- print("\n".join(vcf_lines), file=fh)

- [vcf_a, vcf_b] = roundtrip(tmp_path,

- reference = "ATCGACC",

- meta_lines = [

- "##fileformat=VCFv4.3",

- "##contig=<ID=1,length=7>",

- "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">"

- sample_names = ["sample_A", "sample_B"],

- data_lines = [

- ["foo", "2", ".", "T", "G", ".", ".", ".", "GT", "1", "0"],

- ["foo", "3", ".", "C", "GA,*", ".", ".", ".", "GT", "1", "2"], # "A" insertion will be lost on roundtrip!

- ["foo", "5", ".", "A", "TC", ".", ".", ".", "GT", ".", "1"],

- ["foo", "6", ".", "CC", "C", ".", ".", ".", "GT", "1", "0"], # del of (1-based) pos 7 in sample A

- pass_metadata=False

- assert(vcf_a['reference']==vcf_a['reference'])

- assert(vcf_a['metadata']['chrom'] == "foo" and vcf_b['metadata']['chrom'] == "1")

- assert(vcf_a['metadata']['ploidy'] == 1 and vcf_b['metadata']['ploidy'] == 2)

- assert(vcf_a['positions']==vcf_a['positions'])

- assert(vcf_a['sequences']==vcf_b['sequences'])

- [vcf_a, vcf_b] = roundtrip(tmp_path,

- reference = "ATCGACC",

- meta_lines = [

- "##fileformat=VCFv4.3",

- "##contig=<ID=1,length=7>",

- "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">"

- sample_names = ["sample_A", "sample_B"],

- data_lines = [

- ["1", "2", ".", "T", "G", ".", ".", ".", "GT", "1", "0"],

- ["1", "3", ".", "C", "GA,*", ".", ".", ".", "GT", "1", "2"], # "A" insertion will be lost on roundtrip!

- ["1", "5", ".", "A", "TC", ".", ".", ".", "GT", ".", "1"],

- ["1", "6", ".", "CC", "C", ".", ".", ".", "GT", "1", "0"], # del of (1-based) pos 7 in sample A

- pass_metadata=True

- assert(vcf_a['reference']==vcf_a['reference'])

- assert(vcf_a['metadata']['chrom'] == vcf_b['metadata']['chrom'])

- assert(vcf_a['metadata']['ploidy'] == vcf_b['metadata']['ploidy'])

- assert(vcf_a['positions']==vcf_a['positions'])

- assert(vcf_a['sequences']==vcf_b['sequences'])

- print("designated output location %s is not a directory"%params.outdir, file=sys.stderr)

- outdir = outdir_stem + suffix.rstrip('/')+'/'

- outdir = outdir_stem + '-%04d'%count + suffix.rstrip('/')+'/'

- pos = m.GENOMIC_POSITION-1+offset #put in correct numbering

- DRM_info = {'DRMs': DRMs,

- 'drmPositions': drmPositions}

- print("ERROR: a reference Fasta is required with VCF-format alignments")

- if not hasattr(params, 'gtr') or params.gtr=="infer": #if not specified, set it:

- fixed_pi = [ref.count(base)/len(ref) for base in alpha]

- fixed_pi = [v-0.01 for v in fixed_pi]

- plt.savefig(fname)

- print("--- root-to-tip plot saved to \n\t"+fname)

-

-def export_sequences_and_tree(tt, basename, is_vcf=False, zero_based=False,

- report_ambiguous=False, timetree=False, confidence=False,

- reconstruct_tip_states=False, tree_suffix=''):

- AlignIO.write(tt.get_reconstructed_alignment(reconstruct_tip_states=reconstruct_tip_states), outaln_name, 'fasta')

- print("\n--- alignment including ancestral nodes saved as \n\t %s\n"%outaln_name)

- mutations_out = open(basename + "branch_mutations.txt", "w")

- mutations_out.write("node\tstate1\tpos\tstate2\n")

- fh_dates.write('%s\t--\t--\t--\t--\n'%(n.name))

- fh_dates.write('%s\t%s\t%f\t%f\t%f\n'%(n.name, n.date, n.numdate,conf[0], conf[1]))

- fh_dates.write('%s\t--\t--\n'%(n.name))

- fh_dates.write('%s\t%s\t%f\n'%(n.name, n.date, n.numdate))

- n.confidence=None

- if n.is_terminal() and len(n.name)>40 and bioversion<"1.69":

- n.name = n.name[:35]+'_%03d'%terminal_count

- terminal_count+=1

- n.comment=''

- n.comment= '&mutations="' + ','.join([a+str(pos + offset)+d for (a,pos, d) in n.mutations])+'"'

- n.comment= '&mutations="' + ','.join([a+str(pos + offset)+d for (a,pos, d) in n.mutations

- if tt.gtr.ambiguous not in [a,d]])+'"'

- n.comment= '&mutations="' + ','.join([a+str(pos + offset)+d for (a,pos, d) in n.mutations if n.mask[pos]>0])+f'",mcc="{n.mcc}"'

- n.comment= '&mutations="' + ','.join([a+str(pos + offset)+d for (a,pos, d) in n.mutations

- if tt.gtr.ambiguous not in [a,d] and n.mask[pos]>0])+f'",mcc="{n.mcc}"'

-

- for (a, pos, d) in n.mutations:

- if tt.gtr.ambiguous not in [a,d] or report_ambiguous:

- mutations_out.write("%s\t%s\t%s\t%s\n" %(n.name, a, pos + 1, d))

- n.comment+=(',' if n.comment else '&') + 'date=%1.2f'%n.numdate

- fmt_bl = "%1.7f" if tt.data.full_length<1e6 else "%1.9e"

- print("--- saved divergence times in \n\t %s\n"%dates_fname)

- print("--- tree saved in nexus format as \n\t %s\n"%outtree_name)

- print("--- tree saved in auspice json format as \n\t %s\n"%outtree_name_json)

- print("--- divergence tree saved in nexus format as \n\t %s\n"%outtree_name)

- print("--- saved detected outliers as " + basename + 'outliers.tsv')

- if save: fh = open(save, 'w', encoding='utf-8')

- header1 = "Skyline assuming "+ str(gen)+" gen/year and approximate confidence bounds (+/- %f standard deviations of the LH)\n"%n_std

- header2 = "date \tN_e \tlower \tupper"

- if screen: print('\t'+header1+'\t'+header2)

- if save: fh.write("#"+ header1+'#'+header2+'\n')

- for (x,y, y1, y2) in zip(skyline.x, skyline.y, conf[0], conf[1]):

- if screen: print("\t%1.3f\t%1.3e\t%1.3e\t%1.3e"%(x,y, y1, y2))

- if save: fh.write("%1.3f\t%1.3e\t%1.3e\t%1.3e\n"%(x,y, y1, y2))

- print("\n --- written skyline to %s\n"%save)

- plt.ticklabel_format(axis='x',useOffset=False)

-

- "title": f"Auspice visualization of TreeTime (v{treetime_version}) analysis",

- "build_url": "https://github.com/neherlab/treetime",

- "last_updated": datetime.now().strftime("%Y-%m-%d"),

- "treetime_version": treetime_version,

- "genome_annotations": {

- "nuc":{"start":1, "end":int(tt.data.full_length), "type":"source", "strand":"+:"}

- },

- "panels":["tree", "entropy"],

- "colorings": [

- "title": "Date",

- "type": "continuous",

- "key": "num_date",

- "title": "Genotype",

- "type": "categorical",

- "key": "gt",

- {

- "title": "Excluded",

- "type": "categorical",

- "key": "bad_branch"

- },

- {

- "title": "Branch Support",

- "type": "continuous",

- "key": "confidence"

- }

- "display_defaults": {"color_by":"bad_branch"},

- "filters": ["bad_branch"]

- j = {"name":n.name, "node_attrs":{}, "branch_attrs":{}}

- j["children"] = []

- j["node_attrs"]["num_date"] = {"value":float(n.numdate)}

- j["node_attrs"]["num_date"]["confidence"] = (float(conf[0]), float(conf[1]))

- j["node_attrs"]["div"] = float(pdiv + n.mutation_length)

- j["node_attrs"]["bad_branch"] = {"value": "Yes" if n.bad_branch else "No"}

- if seq_info: # only add mutations to the json if run with sequence data (fasta or vcf)

- j["branch_attrs"]["mutations"] = {"nuc": [f"{a}{pos+1}{d}" for a,pos,d in n.mutations if d in "ACGT-"]}

- j["node_attrs"]["confidence"] = {"value":round(1-np.exp(-len([pos for a,pos,d in n.mutations if d in "ACGT"])),3)

- if not n.is_terminal() else 1.0}

- n_json["children"].append(node_to_json(c, n_json["node_attrs"]["div"]))

- node_lookup[c.name] = n_json["children"][-1]

- return {"meta":meta, "tree":tree}

-version="0.11.4"

-recursion_limit = os.environ.get("TREETIME_RECURSION_LIMIT")

-

-

- matplotlib.use("AGG")

-if __name__ == "__main__":

-from __future__ import division, print_function, absolute_import, absolute_import

- 0.06600504])

- Spis = np.sqrt(pis[None, :] / pis[:,None])

- gtr = GTR(alphabet=alphabets['aa_nogap'])

-

- a1 = {s.id:s for s in AlignIO.read(aln1, 'fasta')}

- a2 = {s.id:s for s in AlignIO.read(aln2, 'fasta')}

- for aln in [a1,a2]:

- for s,seq in aln.items():

- seqstr = "".join(seq2array(seq, fill_overhangs=fill_overhangs))

- l1 = len(a1[leaf])

- l2 = len(a2[leaf])

- return {"MCCs": MCCs, "trees":[t1,t2], "alignment":MultipleSeqAlignment(aln_combined),

- "masks":[mask1,mask2], "combined_mask":combined_mask}

-

-def setup_arg(T, aln, total_mask, segment_mask, dates, MCCs, gtr='JC69',

- verbose=0, fill_overhangs=True, reroot=True, fixed_clock_rate=None, alphabet='nuc', **kwargs):

- tt = TreeTime(dates=dates, tree=T,

- aln=aln, gtr=gtr, alphabet=alphabet, verbose=verbose,

- fill_overhangs=fill_overhangs, keep_node_order=True,

- compress=False, **kwargs)

-

- tt.reroot("least-squares", force_positive=True, clock_rate=fixed_clock_rate)

- for mi,mcc in enumerate(MCCs):

- if (n.mcc is not None) and n.up and n.up.mcc==n.mcc:

- leaf.branch_length = max(0.5*one_mutation, leaf.branch_length)

- n.branch_length = max(0.5*one_mutation, n.branch_length)

- if n==tree.root:

- if n.up.mcc in n.child_mccs: # parent MCC part of children -> that is the MCC

- elif len(n.child_mccs)==1: # child is an MCC

- else: # no unique child MCC and no match with parent -> not part of an MCCs

-from .wrappers import ancestral_reconstruction, mugration, scan_homoplasies,\

- timetree, estimate_clock_model, arg_time_trees

- if len(sys.argv)==1:

-treetime_description = \

- "TreeTime: Maximum Likelihood Phylodynamics\n\n"

-subcommand_description = \

- "In addition, TreeTime implements several sub-commands:\n\n"\

- "\t ancestral\tinfer ancestral sequences maximizing the joint or marginal likelihood.\n"\

- "\t homoplasy\tanalyze patterns of recurrent mutations aka homoplasies.\n"\

- "\t clock\t\testimate molecular clock parameters and reroot the tree.\n"\

- "\t mugration\tmap discrete character such as host or country to the tree.\n\n"\

- "(note that 'tt' is a default subcommand in python2 that doesn't need to be specified).\n"\

- "To print a description and argument list of the individual sub-commands, type:\n\n"\

- "\t treetime <subcommand> -h\n\n"

-

-ref_msg = \

- "If you use results from treetime in a publication, please cite:"\

- "\n\n\tSagulenko et al. TreeTime: Maximum-likelihood phylodynamic analysis"\

- "\n\tVirus Evolution, vol 4, https://academic.oup.com/ve/article/4/1/vex042/4794731\n"

-

-timetree_description=\

- "TreeTime infers a time scaled phylogeny given a tree topology, an alignment, "\

- "and tip dates. Reconstructs ancestral sequences and infers a molecular clock tree. "\

- "TreeTime will reroot the tree and resolve polytomies by default. "\

- "In addition, treetime will infer ancestral sequences and a GTR substitution model. "\

- "Inferred mutations are included as comments in the output tree.\n\n"

-

-gtr_description = "GTR model to use. '--gtr infer' will infer a model "\

- "from the data. Alternatively, specify the model type. If the specified model "\

-gtr_params_description = "GTR parameters for the model specified by "\

- "the --gtr argument. The parameters should be feed as 'key=value' "\

- "list of parameters. Example: '--gtr K80 --gtr-params kappa=0.2 "\

- "pis=0.25,0.25,0.25,0.25'. See the exact definitions of the "\

- "parameters in the GTR creation methods in treetime/nuc_models.py "\

- "or treetime/aa_models.py"

-

-reroot_description = "Reroot the tree using root-to-tip regression. Valid choices are "\

- "'min_dev', 'least-squares', and 'oldest'. 'least-squares' adjusts the root to "\

- "minimize residuals of the root-to-tip vs sampling time regression, " \

- "'min_dev' minimizes variance of root-to-tip distances. "\

- "'least-squares' can be combined with --covariation to account for shared ancestry. "\

- "Alternatively, you can specify a node name or a list of node names "\

- "to be used as outgroup or use 'oldest' to reroot to the oldest node. "\

- "By default, TreeTime will reroot using 'least-squares'. "\

- "Use --keep-root to keep the current root."

-

-tree_description = "Name of file containing the tree in newick, nexus, or phylip format, "\

- "the branch length of the tree should be in units of average number of nucleotide or protein "\

- "substitutions per site. If no file is provided, treetime will attempt "\

- "to build a tree from the alignment using fasttree, iqtree, or raxml "\

- "(assuming they are installed). "

-

-aln_description = "alignment file (fasta)"

-

-dates_description = "csv file with dates for nodes with 'node_name, date' where date is float (as in 2012.15) or in ISO-format (YYYY-MM-DD). "\

- "Imprecisely known dates can be specified as '2023-XX-XX' or [2013.2:2013.7]"

-

-coalescent_description = \

- "coalescent time scale -- sensible values are on the order of the average "\

- "hamming distance of contemporaneous sequences. In addition, 'opt' "\

- "'skyline' are valid options and estimate a constant coalescent rate "\

- "or a piecewise linear coalescent rate history"

-

-ancestral_description = \

- "Reconstructs ancestral sequences and maps mutations to the tree. "\

- "The output consists of a file 'ancestral.fasta' with ancestral sequences "\

- "and a tree 'annotated_tree.nexus' with mutations added as comments "\

- "like A45G,G136T,..., number in SNPs used 1-based index by default. "\

- "The inferred GTR model is written to stdout."

-

-homoplasy_description = \

- "Reconstructs ancestral sequences and maps mutations to the tree. "\

- "The tree is then scanned for homoplasies. An excess number of homoplasies "\

- "might suggest contamination, recombination, culture adaptation or similar."

-

-mugration_description = \

- "Reconstructs discrete ancestral states, for example "\

- "geographic location, host, or similar. In addition to ancestral states, "\

- "a GTR model of state transitions is inferred."

- parser.add_argument('--sequence-length', type=int, help="length of the sequence, "

- "used to calculate expected variation in branch length. "

- "Not required if alignment is provided.")

- parser.add_argument('--vcf-reference', type=str, help='only for vcf input: fasta file of the sequence the VCF was mapped to.')

- parser.add_argument('--clock-filter', type=float, default=4.0,

- help="ignore tips that don't follow a loose clock, "

- "'clock-filter=number of interquartile ranges from regression (method=`residual`)' "

- "or z-score of local clock deviation (method=`local`). "

- "Default=4.0, set to 0 to switch off.")

- parser.add_argument('--clock-filter-method', choices=['residual', 'local'], default='residual',

- help="Use residuals from global clock (`residual`, default) or local clock deviation (`clock`) "

- "to filter out tips that don't follow the clock")

- reroot_group.add_argument('--keep-root', required = False, action="store_true", default=False,

- help ="don't reroot the tree. Otherwise, reroot to minimize the "

- "the residual of the regression of root-to-tip distance and sampling time")

- parser.add_argument('--tip-slack', type=float, default=3,

- help="excess variance associated with terminal nodes accounting for "

- " overdispersion of the molecular clock")

- parser.add_argument('--covariation', action='store_true', help="Account for covariation when estimating rates "

- "or rerooting using root-to-tip regression, default False.")

- parser.add_argument('--aa', action='store_true', help="use aminoacid alphabet")

- parser.add_argument('--custom-gtr', default = None, type=str, help="filename of pre-defined custom GTR model in standard TreeTime format")

- parser.add_argument('--name-column', type=str, help="label of the column to be used as taxon name")

- parser.add_argument('--date-column', type=str, help="label of the column to be used as sampling date")

- parser.add_argument('--keep-overhangs', default = False, action='store_true', help='do not fill terminal gaps')

- parser.add_argument('--zero-based', default = False, action='store_true', help='zero based mutation indexing')

- parser.add_argument('--reconstruct-tip-states', default = False, action='store_true', help='overwrite ambiguous states on tips with the most likely inferred state')

- parser.add_argument('--report-ambiguous', default=False, action="store_true", help='include transitions involving ambiguous states')

- parser.add_argument('--method-anc', default='probabilistic', type=str, choices = ['parsimony', 'fitch', 'probabilistic', 'ml'],

- help="method used for reconstructing ancestral sequences, default is 'probabilistic'")

- parser.add_argument('--verbose', default=1, type=int, help='verbosity of output 0-6')

- parser.add_argument('--outdir', type=str, help='directory to write the output to')

- parser.add_argument('--clock-rate', type=float, help="if specified, the rate of the molecular clock won't be optimized.")

- parser.add_argument('--clock-std-dev', type=float, help="standard deviation of the provided clock rate estimate")

- parser.add_argument('--branch-length-mode', default='auto', type=str, choices=['auto', 'input', 'joint', 'marginal'],

- help="If set to 'input', the provided branch length will be used without modification. "

- "Note that branch lengths optimized by treetime are only accurate at short evolutionary distances.")

- parser.add_argument('--confidence', action='store_true', help="estimate confidence intervals of divergence times using the marginal"

- " posterior distribution, if `--time-marginal` is False (default) inferred divergence"

- " times will still be calculated using the jointly most likely tree configuration.")

- parser.add_argument('--time-marginal', default='false', choices = ['false', 'true', 'assign', 'always', 'only-final', 'never'],

- help="For 'false' or 'never', TreeTime uses the jointly most likely values for the divergence times. "

- "For 'true' and 'always', it uses the marginal inference mode at every round of optimization, for 'only-final' "

- "(or 'assign' for compatibility with previous versions) only uses the marginal "

- "distribution in the final round.")

- parser.add_argument('--keep-polytomies', default=False, action='store_true',

- help="Don't resolve polytomies using temporal information.")

- parser.add_argument('--stochastic-resolve', default=False, action='store_true',

- help="Resolve polytomies using a random coalescent tree.")

- parser.add_argument('--greedy-resolve', action='store_false', dest='stochastic_resolve',

- help="Resolve polytomies greedily. Currently default, but will "

- "switched to `stochastic-resolve` in future versions.")

- parser.add_argument('--relax',nargs=2, type=float,

- help='use an autocorrelated molecular clock. Strength of the gaussian priors on'

- ' branch specific rate deviation and the coupling of parent and offspring'

- ' rates can be specified e.g. as --relax 1.0 0.5. Values around 1.0 correspond'

- ' to weak priors, larger values constrain rate deviations more strongly.'

- ' Coupling 0 (--relax 1.0 0) corresponds to an un-correlated clock.')

- parser.add_argument('--max-iter', default=2, type=int,

- help='maximal number of iterations the inference cycle is run. Note that for polytomy resolution and coalescence models max_iter should be at least 2')

- parser.add_argument('--coalescent', default="0.0", type=str,

- help=coalescent_description)

- parser.add_argument('--n-skyline', default="20", type=int,

- help="number of grid points in skyline coalescent model")

- parser.add_argument('--gen-per-year', default="50.0", type=float,

- help="number of generations per year - used for estimating N_e in coalescent models")

- parser.add_argument('--n-branches-posterior', default=False, action='store_true',

- help= "add posterior LH to coalescent model: use the posterior probability distributions of "

- "divergence times for estimating the number of branches when calculating the coalescent merger"

- "rate or use inferred time before present (default)." )

- parser.add_argument('--plot-tree', default="timetree.pdf",

- help = "filename to save the plot to. Suffix will determine format"

- " (choices pdf, png, svg, default=pdf)")

- parser.add_argument('--plot-rtt', default="root_to_tip_regression.pdf",

- help = "filename to save the plot to. Suffix will determine format"

- " (choices pdf, png, svg, default=pdf)")

- parser.add_argument('--tip-labels', action='store_true',

- help = "add tip labels (default for small trees with <30 leaves)")

- parser.add_argument('--no-tip-labels', action='store_true',

- help = "don't show tip labels (default for small trees with >=30 leaves)")

- parser = argparse.ArgumentParser(description = "",

- usage=treetime_description)

- t_parser.add_argument('--rng-seed', type=int, help="random number generator seed for treetime")

- t_parser.add_argument("--version", action="version", version="%(prog)s " + version)

- print(treetime_description+timetree_description+subcommand_description+

- "'--dates' and '--aln' or '--tree' are REQUIRED inputs, type 'treetime -h' for a full list of arguments.\n")

-

- h_parser.add_argument('--tree', type = str, help=tree_description)

- h_parser.add_argument('--rng-seed', type=int, help="random number generator seed for treetime")

- h_parser.add_argument('--const', type = int, default=0, help ="number of constant sites not included in alignment")

- h_parser.add_argument('--rescale', type = float, default=1.0, help ="rescale branch lengths")

- h_parser.add_argument('--detailed', required = False, action="store_true", help ="generate a more detailed report")

- h_parser.add_argument('--zero-based', default = False, action='store_true', help='zero based mutation indexing')

- h_parser.add_argument('-n', default = 10, type=int, help='number of mutations/nodes that are printed to screen')

- h_parser.add_argument('--drms', type=str, help='TSV file containing DRM info. columns headers: GENOMIC_POSITION, ALT_BASE, DRUG, GENE, SUBSTITUTION')

- a_parser.add_argument('--tree', type=str, help=tree_description)

- a_parser.add_argument('--rng-seed', type=int, help="random number generator seed for treetime")

- a_parser.add_argument('--marginal', default=False, action="store_true", help ="marginal reconstruction of ancestral sequences")

- m_parser.add_argument('--tree', required = True, type=str, help=tree_description)

- m_parser.add_argument('--rng-seed', type=int, help="random number generator seed for treetime")

- m_parser.add_argument('--name-column', type=str, help="label of the column to be used as taxon name")

- m_parser.add_argument('--attribute', type=str, help ="attribute to reconstruct, e.g. country")

- m_parser.add_argument('--states', required = True, type=str, help ="csv or tsv file with discrete characters."

- "\n#name,country,continent\ntaxon1,micronesia,oceania\n...")

- m_parser.add_argument('--weights', type=str, help="csv or tsv file with probabilities of that a randomly sampled "

- "sequence at equilibrium has a particular state. E.g. population of different continents or countries. E.g.:"

- "\n#country,weight\nmicronesia,0.1\n...")

- m_parser.add_argument('--confidence', action="store_true", help="output confidence of mugration inference")

- m_parser.add_argument('--pc', type=float, default=1.0, help ="pseudo-counts higher numbers will results in 'flatter' models")

- m_parser.add_argument('--missing-data', type=str, default='?', help ="string indicating missing data")

- m_parser.add_argument('--sampling-bias-correction', type=float,

- help='a rough estimate of how many more events would have been observed'

- ' if sequences represented an even sample. This should be'

- ' roughly the (1-sum_i p_i^2)/(1-sum_i t_i^2), where p_i'

- ' are the equilibrium frequencies and t_i are apparent ones.'

- '(or rather the time spent in a particular state on the tree)')

-

- c_parser = subparsers.add_parser('clock',

- description="Calculates the root-to-tip regression and quantifies the 'clock-i-ness' of the tree. "

- "It will reroot the tree to maximize the clock-like "

- "signal and recalculate branch length unless run with --keep-root.")

- c_parser.add_argument('--tree', required=True, type=str, help=tree_description)

- c_parser.add_argument('--rng-seed', type=int, help="random number generator seed for treetime")

- c_parser.add_argument('--prune-outliers', action='store_true', default=False,

- help="remove detected outliers from the output tree")

- c_parser.add_argument('--allow-negative-rate', required = False, action="store_true", default=False,

- help="By default, rates are forced to be positive. For trees with little temporal "

- "signal it is advisable to remove this restriction to achieve essentially mid-point rooting.")

- c_parser.add_argument('--plot-rtt', default="root_to_tip_regression.pdf",

- help = "filename to save the plot to. Suffix will determine format"

- " (choices pdf, png, svg, default=pdf)")

- arg_parser = subparsers.add_parser('arg',

- description="Calculates the root-to-tip regression and quantifies the 'clock-i-ness' of the tree. "

- "It will reroot the tree to maximize the clock-like "

- "signal and recalculate branch length unless run with --keep_root.")

- arg_parser.add_argument('--rng-seed', type=int, help="random number generator seed for treetime")

-

- v_parser.set_defaults(func=lambda x: print("treetime "+version))

-class BranchLenInterpolator (Distribution):

- def __init__(self, node, gtr, one_mutation=None, min_width=ttconf.MIN_INTEGRATION_PEAK,

- branch_length_mode = 'joint', pattern_multiplicity = None,

- n_grid_points = ttconf.BRANCH_GRID_SIZE, ignore_gaps=True):

-

- raise Exception("Cannot create branch length interpolator for the root node.")

- one_mutation = 1.0/L

- if mutation_length < np.min((1e-5, 0.1*one_mutation)): # zero-length

- short_range = 10*one_mutation

- grid = np.concatenate([short_range*(np.linspace(0, 1.0 , n_grid_points//2)[:-1]),

- (short_range + (ttconf.MAX_BRANCH_LENGTH - short_range)*(np.linspace(0, 1.0 , n_grid_points//2+1)**2))])

-

- else: # branch length is not zero

- sigma = mutation_length #np.max([self.average_branch_len, mutation_length])

- grid_left = mutation_length * (1 - np.linspace(1, 0.0, n_grid_points//3)**2.0)

- grid_zero = grid_left[1]*np.logspace(-20,0,6)[:5]

- grid_zero2 = grid_left[1]*np.linspace(0,1,10)[1:-1]

- grid_right = mutation_length + (3*sigma*(np.linspace(0, 1, n_grid_points//3)**2))

- far_grid = grid_right.max() + ttconf.MAX_BRANCH_LENGTH*np.linspace(0, 1, n_grid_points//3)**2

- grid = np.concatenate((grid_zero,grid_zero2, grid_left,grid_right[1:],far_grid[1:]))

- grid.sort() # just for safety

- if branch_length_mode=='input':

- p0 = 1.0-np.sum(self.gtr.Pi**2)

- if mutation_length<0.05:

- log_prob = np.array([ k - mutation_length*np.log(k+ttconf.MIN_BRANCH_LENGTH*one_mutation) for k in grid])/one_mutation

- #sigma_sq = (mutation_length+one_mutation)*variance_scale

- l = (mutation_length+one_mutation)

- nm_inv = np.exp(l/p0)

- sigma_sq = p0*(nm_inv-1)*(nm_inv - p0*(nm_inv-1))*one_mutation

- sigma = np.sqrt(sigma_sq+ttconf.MIN_BRANCH_LENGTH*one_mutation)

- log_prob = np.array(np.min([[ 0.5*(mutation_length-k)**2/sigma_sq for k in grid],

- 100 + np.abs([(mutation_length-k)/sigma for k in grid])], axis=0))

- elif branch_length_mode=='marginal':

- log_prob = np.array([-self.gtr.prob_t_profiles(node.profile_pair,

- pattern_multiplicity,

- k,

- return_log=True)

- for k in grid])

- raise Exception("profile pairs need to be assigned to node")

-

- elif branch_length_mode=='joint':

- raise Exception("branch state pairs need to be assigned to nodes")

-

- log_prob = np.array([-self.gtr.prob_t_compressed(node.branch_state['pair'],

- node.branch_state['multiplicity'],

- k,

- return_log=True)

- for k in grid])

- raise Exception("unknown branch length mode! "+branch_length_mode)

- super(BranchLenInterpolator, self).__init__(grid, log_prob, is_log=True,

- kind='linear', min_width=min_width)

-

- ratio = self._gamma/new_gamma

-

-

- if hasattr(node, 'raw_date_constraint') and (node.raw_date_constraint is not None):

- res[node] = node.dist2root - clock_rate*np.mean(node.raw_date_constraint) - icpt

- iqd = np.percentile(residuals,75) - np.percentile(residuals,25)

- for node,r in res.items():

- if abs(r)>n_iqd*iqd and node.up.up is not None:

- node.bad_branch=True

- outliers[node.name] = {'tau':(node.dist2root - icpt)/clock_rate, 'avg_date': np.mean(node.raw_date_constraint),

- 'exact_date': node.raw_date_constraint if type(node) is float else None,

- 'residual': r/iqd}

- node.bad_branch=False

- tt.outliers=None

- outlier_df = pd.DataFrame(outliers).T.loc[:,['avg_date', 'tau', 'residual']]\

- .rename(columns={'avg_date':'given_date', 'tau':'apparent_date'})

- tt.logger("Clock_filter.residual_filter marked the following outliers:", 2, warn=True)

- if tt.verbose>=2:

- tt.logger(f"TreeTime.ClockFilter: starting local_outlier_detection", 2)

- tt.logger(f"TreeTime.ClockFilter: z-scale {z_scale:1.2f}", 2)

- tt.outliers=None

- outlier_df = pd.DataFrame(outliers).T.loc[:,['avg_date', 'tau', 'z', 'diagnosis']]\

- .rename(columns={'avg_date':'given_date', 'tau':'apparent_date', 'z':'z-score'})

- tt.logger("Clock_filter.local_filter marked the following outliers", 2, warn=True)

- if tt.verbose>=2:

- diagnosis=''

- if z<0:

- if np.abs(n_info['avg_date']-parent_tau) > n_info["nmuts"]/mu:

- diagnosis='date_too_early'

- mu = tt.clock_model['slope']*tt.data.full_length

- continue # do not label children of the root as bad -- typically a problem with root choice that will be fixed anyway

- z = (n_info['avg_date'] - n_info['tau'])/z_scale

- zs = [(n_info['tau'] - x)/z_scale for x in n.raw_date_constraint]

- if zs[0]*zs[1]>0 and np.min(np.abs(zs))>z_score_threshold:

- outliers[n.name] = add_outlier_info(zs[0] if np.abs(zs[0])<np.abs(zs[1]) else zs[1],

- n, n_info, parent_tau, mu)

- mu = tt.clock_model['slope']*tt.data.full_length

- sigma_sq = (3/mu)**2

- tt.logger(f"Clockfilter.calculate_node_timings: mu={mu:1.3e}/y, sigma={3/mu:1.3e}y", 2)

- prefactor = (p["observations"]/sigma_sq + mu**2/(p["nmuts"]+eps))

- p["a"] = (p["avg_date"]/sigma_sq + mu*p["nmuts"]/(p["nmuts"]+eps))/prefactor

- children = [node_info[c.name] for c in n if (not node_info[c.name]['skip']) and node_info[c.name]['exact_date']]

- if n==tt.tree.root:

- tmp_children_1 = mu*np.sum([(mu*c["a"]-c["nmuts"])/(eps+c["nmuts"]) for c in children])

- tmp_children_2 = mu**2*np.sum([(1-c["b"])/(eps+c["nmuts"]) for c in children])

- prefactor = (p["observations"]/sigma_sq + tmp_children_2)

- p["a"] = (p["observations"]*p["avg_date"]/sigma_sq + tmp_children_1)/prefactor

- tmp_children_1 = mu*np.sum([(mu*c["a"]-c["nmuts"])/(eps+c["nmuts"]) for c in children])

- tmp_children_2 = mu**2*np.sum([(1-c["b"])/(eps+c["nmuts"]) for c in children])

- prefactor = (p["observations"]/sigma_sq + mu**2/(p["nmuts"]+eps) + tmp_children_2)

- p["a"] = (p["observations"]*p["avg_date"]/sigma_sq + mu*p["nmuts"]/(p["nmuts"]+eps)+tmp_children_1)/prefactor

- p["b"] = mu**2/(p["nmuts"]+eps)/prefactor

- node_info[tt.tree.root.name]["tau"] = node_info[tt.tree.root.name]["a"]

- c_info['tau']=p['tau']

- c_info["tau"] = c_info["a"] + c_info["b"]*p["tau"]

- c_info["tau"] = p["tau"] + c_info['nmuts']/mu

- dev.append(c_info['avg_date']-c_info['tau'])

- date_uncertainty = [np.abs(n.raw_date_constraint[1]-n.raw_date_constraint[0])

- if type(n.raw_date_constraint)!=float else 0.0

- for n in tt.tree.get_terminals()

- if n.raw_date_constraint is not None]

- uncertainty_cutoff = scoreatpercentile(date_uncertainty, percentile_for_exact_date)*1.01

- parent = {"dates": [], "tips": {}, "skip":False}

- child = {'skip':False}

- child["nmuts"] = len([m for m in c.mutations if m[-1] in 'ACGT']) if aln \

- else np.round(c.branch_length*L)

- elif type(c.raw_date_constraint)==float:

- child['exact_date'] = np.abs(c.raw_date_constraint[1]-c.raw_date_constraint[0])<=uncertainty_cutoff

- if child["nmuts"]==0:

- parent["tips"][c.name]={'date': np.mean(c.raw_date_constraint),

- 'exact_date':child['exact_date']}

- child["avg_date"] = np.mean(c.raw_date_constraint)

- parent['exact_date'] = exact_dates>0

- parent["nmuts"] = len([m for m in n.mutations if m[-1] in 'ACGT']) if aln else np.round(n.branch_length*L)

- parent["observations"] = len(d)

- parent["avg_date"] = np.mean(d) if len(d) else 0.0

- def __init__(self, *args, dates=None, debug=False, real_dates=True, precision_fft = 'auto',

- precision='auto', precision_branch='auto', branch_length_mode='joint', use_covariation=False,

- use_fft=True, **kwargs):

-

- raise MissingDataError("ClockTree requires date constraints!")

- self.debug=debug

- # to terminal branches in covariance calculation

- self.clock_model=None

- self.use_covariation=use_covariation # if false, covariation will be ignored in rate estimates.

-

- self.logger("WARNING: ClockTree.init: node %s has a bad date: %s"%(node.name, str(tmp_date)), 2, warn=True)

- self.logger("WARNING: ClockTree.init: node %s has a bad date: %s"%(node.name, str(tmp_date)), 2, warn=True)

- else: # nodes without date contraints

-

- if bad_branch_counter>self.tree.count_terminals()-3:

- raise MissingDataError("ERROR: ALMOST NO VALID DATE CONSTRAINTS")

- self.logger("ClockTree._assign_dates: assigned date constraints to {} out of {} tips.".format(self.tree.count_terminals()-bad_branch_counter, self.tree.count_terminals()), 1)

-

- '''

- '''

- self.min_width = 10*self.one_mutation

- if precision in [0,1,2,3]:

- self.precision=precision

- if self.one_mutation and self.one_mutation<1e-4 and precision<2:

- self.logger("ClockTree._set_precision: FOR LONG SEQUENCES (>1e4) precision>=2 IS RECOMMENDED."

- " precision %d was specified by the user"%precision, level=0)

- if self.one_mutation>1e-4:

- self.precision=1

- self.precision=2

- self.precision=1

- self.logger("ClockTree: Setting precision to level %s"%self.precision, 2)

- if self.precision==0:

- elif self.precision==2:

- elif self.precision==3:

- '''

- function to set the number of grid points for the minimal FWHM window and branch grid

- when calculating the marginal distribution using the FFT-based approach

- The default parameters ttconf.FFT_FWHM_GRID_SIZE and branch_grid_points determined in set_precision

- are used unless an integer value is specified using precision_fft and precision_branch

- '''

-

- if type(precision_fft) is int:

- self.logger("ClockTree.init._set_precision_fft: setting fft grid size explicitly,"

- " fft_grid_points=%.3e"%(precision_fft), 2)

- self.fft_grid_size = precision_fft

- elif precision_fft!='auto':

- raise UnknownMethodError(f"ClockTree: precision_fft needs to be either 'auto' or an integer, got '{precision_fft}'.")

- else:

- self.fft_grid_size = ttconf.FFT_FWHM_GRID_SIZE

- if type(precision_branch) is int:

- self.logger("ClockTree.init._set_precision_fft: setting branch grid size explicitly,"

- " branch_grid_points=%.3e"%(precision_branch), 2)

- self.branch_grid_points = precision_branch

- self.logger("ClockTree.date2dist: Setting new molecular clock."

- " rate=%.3e, R^2=%.4f"%(val.clock_rate, val.r_val**2), 2)

-

- tip_value = lambda x:np.mean(x.raw_date_constraint) if (x.is_terminal() and (x.bad_branch is False)) else None

- branch_value = lambda x:x.mutation_length

- branch_variance = lambda x:((max(0,x.clock_length) if hasattr(x,'clock_length') else x.mutation_length)

- +(self.tip_slack**2*om if x.is_terminal() else 0.0))*om

- branch_variance = lambda x:1.0 if x.is_terminal() else 0.0

- Treg = TreeRegression(self.tree, tip_value=tip_value,

- branch_value=branch_value, branch_variance=branch_variance)

-

- self.logger(f'ClockTree.get_clock_model: estimating clock model with covariation={covariation}',3)

- raise ValueError("Clock rate estimation failed. If your data lacks temporal signal, please specify the rate explicitly!")

- self.clock_model['valid_confidence']=False

- self.clock_model['valid_confidence']=True

-

- self.logger("ClockTree.init_date_constraints...",2)

- self.infer_ancestral_sequences('probabilistic', marginal=self.branch_length_mode=='marginal',

- sample_from_profile='root',**kwarks)

- self.logger('ClockTree.init_date_constraints: Initializing branch length interpolation objects...',2)

- if hasattr(node,'branch_length_interpolator') and node.branch_length_interpolator is not None:

- if self.branch_length_mode=='marginal':

- elif self.branch_length_mode=='joint' and (not hasattr(node, 'branch_state')):

- node.branch_length_interpolator = BranchLenInterpolator(node, self.gtr,

- pattern_multiplicity = self.data.multiplicity(mask=node.mask), min_width=self.min_width,

- one_mutation=self.one_mutation, branch_length_mode=self.branch_length_mode,

- n_grid_points = self.branch_grid_points)

- use_cov = (np.sum(has_clock_length) > len(has_clock_length)*0.7) and self.use_covariation

- for node in self.tree.find_clades(order="postorder"):

- self.logger("ClockTree.init_date_constraints -- WARNING: Branch is marked as bad"

- ", excluding it from the optimization process."

- " Date constraint will be ignored!", 4, warn=True)

- else: # node without sampling date set

-

- '''

- '''

- self.logger("ClockTree: Maximum likelihood tree optimization with temporal constraints",1)

-

-

- self.logger("ClockTree - Joint reconstruction: Propagating leaves -> root...", 2)

- else: # all other nodes

- if node.date_constraint is not None and node.date_constraint.is_delta: # there is a strict time constraint

- node.joint_pos_Lx = Distribution(x, -self.merger_model.integral_merger_rate(node.date_constraint.peak_pos)

- + node.branch_length_interpolator(bl), min_width=self.min_width, is_log=True)

- node.joint_pos_Lx = Distribution(x, node.branch_length_interpolator(bl), min_width=self.min_width, is_log=True)

- node.joint_pos_Cx = Distribution(x, bl, min_width=self.min_width) # map back to the branch length

- else: # all nodes without precise constraint but positional information

- msgs_to_multiply.append(Distribution(time_points, -self.merger_model.integral_merger_rate(time_points), is_log=True))

- if len(msgs_to_multiply) == 0: # there are no constraints

- elif len(msgs_to_multiply)>1: # combine the different msgs and constraints

- else: # there is exactly one constraint.

- if node.up is None: # this is the root, set dates

- subtree_distribution = Distribution.multiply([subtree_distribution, Distribution(subtree_distribution.x,

- self.merger_model.integral_merger_rate(subtree_distribution.x), is_log=True)])

- else: # otherwise propagate to parent

- res, res_t = NodeInterpolator.convolve(subtree_distribution,

- node.branch_length_interpolator,

- max_or_integral='max',

- inverse_time=True,

- n_grid_points = self.node_grid_points,

- n_integral=self.n_integral,

- rel_tol=self.rel_tol_refine)

- node.joint_inverse_cdf=interp1d([0,1], node.date_constraint.peak_pos*np.ones(2), kind="linear")

- int_y = np.concatenate(([0], np.cumsum(dt*(y[1:]+y[:-1])/2.0)))

- int_y/=int_y[-1]

- node.joint_inverse_cdf = interp1d(int_y, subtree_distribution.x, kind="linear")

- #node.joint_cdf = interp1d(subtree_distribution.x, int_y, kind="linear")

-

- self.logger("ClockTree - Joint reconstruction: Propagating root -> leaves...", 2)

- if node.up is None: # root node

- continue # the position was already set on the previous step

-

- if node.joint_pos_Cx is None: # no constraints or branch is bad - reconstruct from the branch len interpolator

- merger_cost = Distribution(node.branch_length_interpolator.x, self.merger_model.cost(node.time_before_present,

- node.branch_length_interpolator.x, multiplicity=len(node.up.clades)), is_log=True)

- node.branch_length = node.joint_pos_Cx(max(node.joint_pos_Cx.xmin,

- node.up.time_before_present))

- if node.branch_length<0:

- if node.branch_length>-2*ttconf.TINY_NUMBER:

- self.logger(f"ClockTree - Joint reconstruction: correcting rounding error of {node.name} bl={node.branch_length:1.2e}", 4)

- self.logger(f"ClockTree - Joint reconstruction: NEGATIVE BRANCH LENGTH {node.name} bl={node.branch_length:1.2e}", 2, warn=True)

-

- '''

- '''

- LH = self.tree.root.marginal_pos_LH.integrate(return_log=True, a=self.tree.root.marginal_pos_LH.xmin, b=self.tree.root.marginal_pos_LH.xmax, n=1000)

- if node.up is None: # root node

-

- #del node.marginal_pos_LH

- self.logger(f"ClockTree - Marginal reconstruction using {method} convolution: Propagating leaves -> root...", 2)

- else: # all other nodes

- if node.date_constraint is not None and node.date_constraint.is_delta: # there is a hard time constraint

- node.marginal_pos_Lx = Distribution(x, -self.merger_model.integral_merger_rate(node.date_constraint.peak_pos)

- +node.branch_length_interpolator(bl), min_width=self.min_width, is_log=True)

- node.marginal_pos_Lx = Distribution(x, node.branch_length_interpolator(bl), min_width=self.min_width, is_log=True)

- else: # all nodes without precise constraint but positional information

- # subtree likelihood given the node's constraint and child msg:

- msgs_to_multiply.extend([child.marginal_pos_Lx for child in node.clades

- if child.marginal_pos_Lx is not None])

- if len(msgs_to_multiply)==0:

- elif len(msgs_to_multiply)==1:

- else: # combine the different msgs and constraints

- merger_contribution = Distribution(time_points, -self.merger_model.integral_merger_rate(time_points), is_log=True)

- node.subtree_distribution = Distribution.multiply([merger_contribution, node.product_of_child_messages])

- if node.up is None: # this is the root, set dates

- node.marginal_pos_LH = Distribution.multiply([node.subtree_distribution, Distribution(node.subtree_distribution.x,

- self.merger_model.integral_merger_rate(node.subtree_distribution.x), is_log=True)])

- else: # otherwise propagate to parent

- res, res_t = NodeInterpolator.convolve_fft(node.subtree_distribution,

- node.branch_length_interpolator, self.fft_grid_size), None

- res, res_t = NodeInterpolator.convolve(node.subtree_distribution,

- node.branch_length_interpolator,

- max_or_integral='integral',

- n_grid_points = self.node_grid_points,

- n_integral=self.n_integral,

- rel_tol=self.rel_tol_refine)

- self.logger("ClockTree - Marginal reconstruction: Propagating root -> leaves...", 2)

- node.msg_from_parent = None #if internal node has a delta constraint no previous information is passed on

- node.msg_from_parent = None # nothing beyond the root

- msg_parent_to_node =None

- if len(parent.clades)<5:

- complementary_msgs.extend([sister.marginal_pos_Lx for sister in parent.clades

- if (sister != node) and (sister.marginal_pos_Lx is not None)])

- complementary_msgs = [Distribution.divide(parent.product_of_child_messages, node.marginal_pos_Lx)]

- if len(time_points)<5:

- time_points = np.unique(np.concatenate([

- time_points,

- np.linspace(np.min([x.xmin for x in complementary_msgs]),

- np.max([x.xmax for x in complementary_msgs]), 50),

- np.linspace(np.min([x.effective_support[0] for x in complementary_msgs]),

- np.max([x.effective_support[1] for x in complementary_msgs]), 50),

- ]))

- complementary_msgs.append(Distribution(time_points, self.merger_model.integral_merger_rate(time_points), is_log=True))