def get_pivots(observations, pivot_interval, start_date=None, end_date=None, pivot_interval_units="months"):
    """Calculate pivots for a given list of floating point observation dates and
    interval between pivots.

    Start and end pivots will be based on the range of given observed dates,
    unless a start or end date are provided to override these defaults.

    Parameters
    ----------
    observations : list
        a list of observed floating point dates per sample
    pivot_interval : int
        number of months (or weeks) between pivots
    start_date : float
        optional start of the pivots interval
    end_date : float
        optional end of the pivots interval
    pivot_interval : str
        whether pivots are measured in "months" or in "weeks"

    Returns
    -------
    pivots : ndarray
        floating point pivots spanning the given the dates

    """
    # Convert months between pivots to pivot frequency.
    pivot_frequency = pivot_interval / 12.0
    if pivot_interval_units == "weeks":
        pivot_frequency = pivot_interval / 52.1429

    pivot_start = start_date if start_date else np.floor(np.min(observations) / pivot_frequency) * pivot_frequency
    pivot_end = end_date if end_date else np.ceil(np.max(observations) / pivot_frequency) * pivot_frequency

    if pivot_interval_units == "months":
        offset = "%sMS" % pivot_interval
    elif pivot_interval_units == "weeks":
        offset = "%sW" % pivot_interval
    else:
        raise ValueError(f"The given interval unit '{pivot_interval_units}' is not supported.")

    datetime_pivots = pd.date_range(
        float_to_datestring(pivot_start),
        float_to_datestring(pivot_end),
        freq = offset
    )
    pivots = np.array([timestamp_to_float(pivot) for pivot in datetime_pivots])

    return np.around(pivots, 4) class alignment_frequencies(object):
    '''
    calculates frequencies of mutations in an alignment. uses nested frequencies
    for mutations at a particular site such that mutations are forced to add
    up to one.

    '''
    def __init__(self, aln, tps, pivots, **kwargs):
        """Create an instance of the alignment frequency estimate

        Parameters
        ----------
        aln : Bio.Align.MultipleSeqAlignment
            alignment
        tps : np.array(float)
            Array of numerical dates, one for each sequence in the
            alignment in the SAME ORDER!
        pivots : np.array(float)
            pivot values for which frequencies are estimated
        **kwargs
            Description
        """
        self.aln = np.array(aln)
        self.tps = np.array(tps)
        self.kwargs = kwargs
        self.pivots = pivots
        self.counts = count_observations(self.pivots, self.tps)


    def estimate_genotype_frequency(self, gt):
        '''
        slice an alignment at possibly multiple positions and calculate the
        frequency trajectory of this multi-locus genotype

        Parameters
        ----------
        gt : list
            a list of (position, state) tuples specifying the genotype
            whose frequency is to be estimated

        Returns
        -------
        np.array
            frequency trajectory
        '''
        match = []
        for pos, state in gt:
            match.append(aln[:,pos]==state)
        obs = np.array(match).all(axis=0)

        fe = frequency_estimator(zip(self.tps, obs),self.pivots, **kwargs)
        fe.learn()
        return fe.frequency_estimate


    def mutation_frequencies(self, min_freq=0.01, include_set=None, ignore_char=''):
        '''
        estimate frequencies of single site mutations for each alignment column.
        This function populates a dictionary class.frequencies with the frequency
        trajectories.

        Parameters
        ----------
        min_freq : float, optional
            minimal all-time frequency for an aligment column to be considered
        include_set : list/set, optional
            set of alignment column that will be used regardless of variation
        ignore_char : str, optional
            ignore this character in an alignment column (missing data)
        '''
        if include_set is None:
            include_set=[]
        alphabet = np.unique(self.aln)
        # af: rows correspond to letters of the alphabet, columns positions in the alignemtn, values: frequencies
        af = np.zeros((len(alphabet), self.aln.shape[1]))
        for ni, n in enumerate(alphabet):
            af[ni] = (self.aln==n).mean(axis=0)

        # minor_freqs: frequencies of all alleles not the major one! ignore_char is excluded
        if ignore_char:
            minor_freqs = af[alphabet!=ignore_char].sum(axis=0) - af[alphabet!=ignore_char].max(axis=0)
        else:
            minor_freqs = 1.0 - af.max(axis=0)
        self.frequencies = {}
        for pos in sorted(set.union(set(np.where(minor_freqs>min_freq)[0]), include_set)):
            # indices determine the indicies of mutations in descending order by frequency
            nis = np.argsort(af[:,pos])[::-1]
            nis = nis[af[nis,pos]>0]
            column = self.aln[:,pos]
            if ignore_char: # subset sequences, time points and alphabet to non-gapped and non X
                good_seq = (column!=ignore_char)&(column!='X') #array of bools, length = # seqs in aln
                tps=self.tps[good_seq]
                column = column[good_seq]
                nis = nis[(alphabet[nis]!=ignore_char)&(alphabet[nis]!='X')]
            else:
                tps = self.tps

            if len(nis)==0:
                # no characters with any frequencies (column of ignore chars)
                continue

            muts = alphabet[nis]
            obs = {}
            for ni, mut in zip(nis, muts):
                if (af[ni,pos]>min_freq and af[ni,pos]<1-min_freq) or ni==0 or (pos in include_set):
                    obs[(pos, mut)] = column==mut
                else:
                    break

            # if multiple states are below frequency threshold, glob them together as 'other'
            if len(obs)==0:
                obs[(pos, muts[0])] = column==muts[0]
            elif len(obs)!=len(nis):
                tmp = ~np.any(list(obs.values()), axis=0) # pull out sequences not yet assigned

                if any(tmp):
                    if len(obs)==len(nis)-1: # if only category left, assign it
                        obs[(pos, muts[-1])] = tmp
                    else: #other wise, call that mutation 'other'
                        obs[(pos, 'other')] = tmp

            print("Estimating frequencies of position: {}".format(pos), end="\r")
            sys.stdout.flush()
            # print("Variants found at frequency:", [(k,o.mean()) for k,o in obs.items()])

            # calculate frequencies, which will be added to the frequencies dict with (pos, mut) as key
            ne = nested_frequencies(tps, obs, self.pivots, **self.kwargs)
            self.frequencies.update(ne.calc_freqs())


    def calc_confidence(self):
        """calculate a crude binomial sampling confidence interval
        of the frequency estimate. This ignores autocorrelation of the
        trajectory and returns one standard deviation (68%).

        Returns
        -------
        dict
            dictionary of standard deviations for each estimated trajectory
        """
        self.confidence = {}
        for key, freq in self.frequencies.items():
            # add a pseudo count 1/(n+1) and normalize to n+1
            self.confidence[key] = np.sqrt((1.0/(1+self.counts)+freq*(1-freq))/(1.0+self.counts))

        return self.confidence class tree_frequencies(object):
    '''
    class that estimates frequencies for nodes in the tree. each internal node is assumed
    to be named with an attribute clade, of root doesn't have such an attribute, clades
    will be numbered in preorder. Each node is assumed to have an attribute `attr` with a
    key "num_date".
    '''
    def __init__(self, tree, pivots, node_filter=None, min_clades=10, verbose=0, pc=1e-4, **kwargs):
        '''
        set up the internal tree, the pivots and cutoffs

        Parameters
        ----------
        tree : Bio.Phylo.calde
            Biopython tree
        pivots : int/array
            number or list of pivots
        node_filter : callable, optional
            function that evaluates to true/false to filter nodes
        min_clades : int, optional
            minimal size of clades to estimate proper frequencies
        verbose : int, optional
            verbosity
        pc : float, optional
            pseudo-counts/lower cutoff to drive estimates away from 0/1
        **kwargs
            Description
        '''
        self.tree = tree
        self.min_clades = min_clades
        self.pivots = pivots
        self.kwargs = kwargs
        self.verbose = verbose
        self.pc = pc
        if node_filter is None:
            self.node_filter = lambda x:True
        else:
            self.node_filter = node_filter
        self.prepare()


    def prepare(self):
        # name nodes if they aren't already named
        if not hasattr(self.tree.root, 'clade'):
            for ni,node in enumerate(self.tree.find_clades(order='preorder')):
                node.clade = ni

        # extract all observations
        tps = []
        leaf_count = 0
        for node in self.tree.find_clades(order='postorder'):
            if node.is_terminal():
                if self.node_filter(node):
                    tps.append(node.attr["num_date"])
                    node.leafs = np.array([leaf_count], dtype=int)
                    leaf_count+=1
                else:
                    node.leafs = np.array([], dtype=int)
            else:
                node.leafs = np.concatenate([c.leafs for c in node.clades])
        self.tps = np.array(tps)

        if np.isscalar(self.pivots):
            self.frequencies = {self.tree.root.clade:np.ones(self.pivots)}
        else:
            self.frequencies = {self.tree.root.clade:np.ones_like(self.pivots)}

        self.counts = count_observations(self.pivots, self.tps)


    def estimate_clade_frequencies(self):
        for node in self.tree.get_nonterminals(order='preorder'):
            if self.verbose: print("Estimating frequencies of children of node ",node.clade)
            node_tps = self.tps[node.leafs]
            obs_to_estimate = {}
            small_clades = []
            if len(node.clades)==1:
                self.frequencies[node.clades[0].clade] = self.frequencies[node.clade]
                continue

            for c in node.clades:
                if len(c.leafs)>self.min_clades:
                    obs_to_estimate[c.clade] = np.in1d(node.leafs, c.leafs)
                else:
                    # Only include internal nodes or tips that pass the node
                    # filter as small clades.
                    if not c.is_terminal() or self.node_filter(c):
                        small_clades.append(c)
            if len(obs_to_estimate):
                if len(small_clades):
                    remainder = {}
                    for c in small_clades:
                        remainder[c.clade] = np.in1d(node.leafs, c.leafs)
                    if len(small_clades)==1:
                        obs_to_estimate.update(remainder)
                    else:
                        obs_to_estimate['other'] = np.any(list(remainder.values()), axis=0)

                ne = nested_frequencies(node_tps, obs_to_estimate, self.pivots, pc=self.pc, **self.kwargs)
                freq_est = ne.calc_freqs()
                for clade, tmp_freq in freq_est.items():
                    if clade != "other":
                        self.frequencies[clade] = self.frequencies[node.clade] * tmp_freq

                if len(small_clades) > 1:
                    total_leaves_in_small_clades = 0
                    for clade in small_clades:
                        total_leaves_in_small_clades += len(clade.leafs)

                    for clade in small_clades:
                        if total_leaves_in_small_clades > 0:
                            frac = len(clade.leafs) / total_leaves_in_small_clades
                        else:
                            frac = 0.0

                        self.frequencies[clade.clade] = frac * freq_est["other"] * self.frequencies[node.clade]
            else:
                for clade in small_clades:
                    if len(node.leafs):
                        frac = 1.0*len(clade.leafs)/len(node.leafs)
                    else:
                        frac = 0.0
                    self.frequencies[clade.clade] = frac*self.frequencies[node.clade]

        # Assign zero frequencies to tips that did not pass the node_filter.
        for tip in self.tree.get_terminals():
            if not self.node_filter(tip):
                self.frequencies[tip.clade] = np.zeros(len(self.pivots))


    def calc_confidence(self):
        '''
        for each frequency trajectory, calculate the bernouilli sampling error
        -- in reality we should look at the inverse hessian, but this is a
        useful approximation in most cases

        Returns
        -------
        dict
            dictionary with estimated confidence intervals
        '''
        self.confidence = {}
        for key, freq in self.frequencies.items():
            # add a pseudo count 1/(n+1) and normalize to n+1
            self.confidence[key] = np.sqrt((1.0/(1+self.counts)+freq*(1-freq))/(1.0+self.counts))
        return self.confidence class AlignmentKdeFrequencies(KdeFrequencies):
    """KDE frequency estimator for multiple sequence alignments

    Estimates frequencies for samples provided as sequences in a multiple
    sequence alignment and corresponding observation dates for each sequence.
    """
    def estimate(self, alignment, observations):
        """Estimate frequencies of bases/residues at each site in a given multiple sequence
        alignment based on the observation dates associated with each sequence
        in the alignment.
        """
        # Calculate pivots for the given tree.
        self.pivots = get_pivots(
            observations,
            self.pivot_frequency,
            start_date=self.start_date,
            end_date=self.end_date,
            pivot_interval_units=self.pivot_interval_units
        )

        # Pair alignment sequence indices with observation dates.
        samples = zip(range(len(alignment)), observations)

        # Sort samples by date and extract sequence indices and dates.
        samples = np.array(sorted(samples, key=lambda row: row[1]))
        sample_indices = samples[:, 0]
        sample_dates = samples[:, 1].astype(float)

        # Map sample alignment indices to their corresponding frequency matrix row indices.
        sample_to_index = {sample_indices[i]: i for i in range(len(samples))}

        # Calculate tip frequencies and normalize.
        normalized_freq_matrix = self.estimate_frequencies(
            sample_dates,
            self.pivots,
            sigma_narrow=self.sigma_narrow,
            sigma_wide=self.sigma_wide,
            proportion_wide=self.proportion_wide,
            max_date=self.max_date,
            censored=self.censored
        )

        # Calculate frequencies per site and base/residue from the estimated
        # frequencies per sample.
        frequencies = {}
        for position in range(alignment.get_alignment_length()):
            # Only estimate frequencies at sites with variation.
            if len(set(alignment[:, position])) > 1:
                for sample, base in enumerate(alignment[:, position]):
                    # Estimate frequencies by current position and base/residue
                    # at the position for each sample.
                    key = "%s:%s" % (position + 1, base.upper())
                    if key not in frequencies:
                        frequencies[key] = np.zeros_like(self.pivots)

                    frequencies[key] += normalized_freq_matrix[sample_to_index[sample], :]

        self.frequencies = frequencies

        return frequencies class TreeKdeFrequencies(KdeFrequencies):
    """KDE frequency estimator for phylogenetic trees

    Estimates frequencies for samples provided as annotated tips on a given
    tree and optionally reports clade-level frequencies based on the sum of each
    clade's respective tips.
    """
    def tip_passes_filters(self, tip):
        """Returns a boolean indicating whether a given tip passes the node filters
        defined for the current instance.

        If no filters are defined, returns True.

        Args:
            tip (Bio.Phylo): tip from a Bio.Phylo tree annotated with attributes in `tip.attr`

        Returns:
            bool: whether the given tip passes the defined filters or not
        """
        return (self.node_filters is None or
                all([tip.attr[key] in values for key, values in self.node_filters.items()]))

    def estimate_tip_frequencies_to_proportion(self, tips, proportion):
        """Estimate frequencies for a given set of tips up to a given proportion of total frequencies.

        Args:
            tips (list): a list of Bio.Phylo terminals annotated with attributes in `tip.attr`
            proportion (float): the proportion of the total frequency that the given tips should represent

        Returns:
            dict: frequencies of given tips indexed by tip name
        """
        clade_frequencies = {}

        # If none of the tips pass the given node filters, do not try to
        # normalize tip frequencies.
        if len(tips) == 0:
            return clade_frequencies

        tips = np.array(sorted(tips, key=lambda row: row[1]))
        clades = tips[:, 0]
        tip_dates = tips[:, 1].astype(float)

        # Map clade ids to their corresponding frequency matrix row index.
        clade_to_index = {clades[i]: i for i in range(len(clades))}

        # Calculate tip frequencies and normalize.
        normalized_freq_matrix = self.estimate_frequencies(
            tip_dates,
            self.pivots,
            normalize_to=proportion,
            sigma_narrow=self.sigma_narrow,
            sigma_wide=self.sigma_wide,
            proportion_wide=self.proportion_wide,
            max_date=self.max_date,
            censored=self.censored
        )

        for clade in clades:
            clade_frequencies[clade] = normalized_freq_matrix[clade_to_index[clade]]

        return clade_frequencies

    def estimate(self, tree):
        """Estimate frequencies for a given tree using the parameters defined for this instance.

        If weights are defined, frequencies represent a weighted mean across the
        values in attribute defined by `self.weights_attribute`.

        Args:
            tree (Bio.Phylo): annotated tree whose nodes all have an `attr` attribute with at least  "num_date" key

        Returns:
            frequencies (dict): node frequencies by clade

        """
        # Calculate pivots for the given tree.
        observations = [tip.attr["num_date"] for tip in tree.get_terminals()]
        self.pivots = get_pivots(
            observations,
            self.pivot_frequency,
            start_date=self.start_date,
            end_date=self.end_date,
            pivot_interval_units=self.pivot_interval_units
        )

        # If weights are defined, calculate frequencies for all tips by their
        # weight attribute values such that the total frequency of all tips with
        # a given value sums to the proportion of tips that have that value.
        # If weights are not defined, estimate frequencies such that they sum to 1.
        if self.weights:
            # Determine which weight attributes are represented in the
            # tree. Drop any unrepresented attributes and renormalize the
            # remaining proportions per attribute to sum to one before
            # estimating frequencies.
            weights_represented = set([
                tip.attr[self.weights_attribute]
                for tip in tree.find_clades(terminal=True)
            ])

            if len(self.weights) != len(weights_represented):
                # Remove unrepresented weights.
                weights_unrepresented = set(self.weights.keys()) - weights_represented
                for weight in weights_unrepresented:
                    del self.weights[weight]

                # Renormalize the remaining weights.
                weight_total = sum(self.weights.values())
                for key, value in self.weights.items():
                    self.weights[key] = value / weight_total

            # Confirm that one or more weights are represented by tips in the
            # tree. If there are no more weights, raise an exception because
            # this likely represents a data error (either in the tree
            # annotations or the weight definitions).
            if len(self.weights) == 0:
                raise TreeKdeFrequenciesError("None of the provided weight attributes were represented by tips in the given tree. Doublecheck weight attribute definitions and their representations in the tree.")

            # Estimate frequencies for all tips within each weight attribute
            # group.
            weight_keys, weight_values = zip(*sorted(self.weights.items()))
            proportions = np.array(weight_values) / np.array(weight_values).sum(axis=0)
            frequencies = {}

            for (weight_key, proportion) in zip(weight_keys, proportions):
                # Find tips with the current weight attribute.
                tips = [(tip.name, tip.attr["num_date"])
                        for tip in tree.get_terminals()
                        if tip.attr[self.weights_attribute] == weight_key and self.tip_passes_filters(tip)]
                frequencies.update(self.estimate_tip_frequencies_to_proportion(tips, proportion))
        else:
            tips = [(tip.name, tip.attr["num_date"])
                    for tip in tree.get_terminals()
                    if self.tip_passes_filters(tip)]
            frequencies = self.estimate_tip_frequencies_to_proportion(tips, proportion=1.0)

        # Assign zero frequencies to any tips that were filtered out of the frequency estimation.
        for tip in tree.get_terminals():
            if not tip.name in frequencies:
                frequencies[tip.name] = np.zeros_like(self.pivots)

        if self.include_internal_nodes:
            for node in tree.find_clades(order="postorder"):
                if not node.is_terminal():
                    frequencies[node.name] = np.array(
                        [frequencies[child.name] for child in node.clades]
                    ).sum(axis=0)

        # Store frequencies in the current instance for simplified exporting.
        self.frequencies = frequencies

        return frequencies class TreeKdeFrequenciesError(Exception):
    """Represents an error estimating KDE frequencies for a tree.
    """
    pass def read_metadata(fname, query=None, as_data_frame=False):
    return MetadataFile(fname, query, as_data_frame).read() def read_node_data(fnames, tree=None):
    return NodeDataReader(fnames, tree).read() def write_json(data, file_name, indent=(None if os.environ.get("AUGUR_MINIFY_JSON") else 2), include_version=True):
    """
    Write ``data`` as JSON to the given ``file_name``, creating parent directories
    if necessary. The augur version is included as a top-level key "augur_version".

    Parameters
    ----------
    data : dict
        data to write out to JSON
    file_name : str
        file name to write to
    indent : int or None, optional
        JSON indentation level. Default is `None` if the environment variable `AUGUR_MINIFY_JSON`
        is truthy, else 1
    include_version : bool, optional
        Include the augur version. Default: `True`.

    Raises
    ------
    OSError
    """
    #in case parent folder does not exist yet
    parent_directory = os.path.dirname(file_name)
    if parent_directory and not os.path.exists(parent_directory):
        try:
            os.makedirs(parent_directory)
        except OSError: #Guard against race condition
            if not os.path.isdir(parent_directory):
                raise

    if include_version:
        data["generated_by"] = {"program": "augur", "version": get_augur_version()}

    with open(file_name, 'w', encoding='utf-8') as handle:
        json.dump(data, handle, indent=indent, sort_keys=True) def get_numerical_dates(meta_dict, name_col = None, date_col='date', fmt=None, min_max_year=None):
    if fmt:
        numerical_dates = {}

        if isinstance(meta_dict, dict):
            for k,m in meta_dict.items():
                v = m[date_col]
                numerical_dates[k] = get_numerical_date_from_value(
                    v,
                    fmt,
                    min_max_year
                )
        elif isinstance(meta_dict, pd.DataFrame):
            strains = meta_dict.index.values
            dates = meta_dict[date_col].apply(
                lambda date: get_numerical_date_from_value(
                    date,
                    fmt,
                    min_max_year
                )
            ).values
            numerical_dates = dict(zip(strains, dates))
    else:
        if isinstance(meta_dict, dict):
            numerical_dates = {k:float(v) for k,v in meta_dict.items()}
        elif isinstance(meta_dict, pd.DataFrame):
            strains = meta_dict.index.values
            dates = meta_dict[date_col].astype(float)
            numerical_dates = dict(zip(strains, dates))

    return numerical_dates
import json, os, sys
import datetime
import numpy as np
import treetime.utils
from collections import defaultdict
from Bio import Phylo, AlignIO
from Bio.Align import MultipleSeqAlignment
from .frequency_estimators import get_pivots, alignment_frequencies, tree_frequencies
from .frequency_estimators import AlignmentKdeFrequencies, TreeKdeFrequencies, TreeKdeFrequenciesError
from .utils import read_metadata, read_node_data, write_json, get_numerical_dates
    parser.add_argument("--include-internal-nodes", action="store_true",
                        help="calculate frequencies for internal nodes as well as tips")

    # Alignment-specific arguments
    parser.add_argument('--alignments', type=str, nargs='+',
                        help="alignments to estimate mutations frequencies for")
    parser.add_argument('--gene-names', nargs='+', type=str,
                        help="names of the sequences in the alignment, same order assumed")
    parser.add_argument('--ignore-char', type=str, default='',
                        help="character to be ignored in frequency calculations")
    parser.add_argument('--minimal-frequency', type=float, default=0.05,
                        help="minimal all-time frequencies for a trajectory to be estimates")

    # KDE-specific arguments
    parser.add_argument("--narrow-bandwidth", type=float, default=1 / 12.0, help="the bandwidth for the narrow KDE")
    parser.add_argument("--wide-bandwidth", type=float, default=3 / 12.0, help="the bandwidth for the wide KDE")
    parser.add_argument("--proportion-wide", type=float, default=0.2, help="the proportion of the wide bandwidth to use in the KDE mixture model")
    parser.add_argument("--weights", help="a dictionary of key/value mappings in JSON format used to weight KDE tip frequencies")
    parser.add_argument("--weights-attribute", help="name of the attribute on each tip whose values map to the given weights dictionary")
    parser.add_argument("--censored", action="store_true", help="calculate censored frequencies at each pivot")

    # Diffusion frequency specific arguments
    parser.add_argument('--minimal-clade-size', type=int, default=0,
                        help="minimal number of tips a clade must have for its diffusion frequencies to be reported")
    parser.add_argument('--minimal-clade-size-to-estimate', type=int, default=10,
                        help="""minimal number of tips a clade must have for its diffusion frequencies to be estimated
                                by the diffusion likelihood; all smaller clades will inherit frequencies from their
                                parents""")
    parser.add_argument("--stiffness", type=float, default=10.0, help="parameter penalizing curvature of the frequency trajectory")
    parser.add_argument("--inertia", type=float, default=0.0, help="determines how frequencies continue "
                        "in absense of data (inertia=0 -> go flat, inertia=1.0 -> continue current trend)")

    # Output arguments
    parser.add_argument('--output-format', default='auspice', choices=['auspice', 'nextflu'],
                        help="format to export frequencies JSON depending on the viewing interface")
    parser.add_argument('--output', '-o', type=str,
                        help='JSON file to save estimated frequencies to')


def format_frequencies(freq):
    return [round(x,6) for x in freq]


def run(args):
    metadata, columns = read_metadata(args.metadata)
    dates = get_numerical_dates(metadata, fmt='%Y-%m-%d')
    stiffness = args.stiffness
    inertia = args.inertia

    if args.method == "kde":
        # Load weights if they have been provided.
        if args.weights:
            with open(args.weights, "r", encoding='utf-8') as fh:
                weights = json.load(fh)

            weights_attribute = args.weights_attribute
        else:
            weights = None
            weights_attribute = None

    if args.tree:
        tree = Phylo.read(args.tree, 'newick')
        tps = []
        for tip in tree.get_terminals():
            tip.attr = {"num_date": np.mean(dates[tip.name])}
            tps.append(tip.attr["num_date"])

            # Annotate tips with metadata to enable filtering and weighting of
            # frequencies by metadata attributes.
            for key, value in metadata[tip.name].items():
                tip.attr[key] = value

        if args.method == "diffusion":
            # estimate tree frequencies
            pivots = get_pivots(tps, args.pivot_interval, args.min_date, args.max_date, args.pivot_interval_units)
            frequency_dict = {"pivots":format_frequencies(pivots)}
            frequency_dict["counts"] = {}

            for region in args.regions:
                # Omit strains sampled prior to the first pivot from frequency calculations.
                # (these tend to be reference strains included for phylogenetic context)
                if region=='global':
                    node_filter_func = lambda node: node.attr["num_date"] >= pivots[0]
                else:
                    node_filter_func = lambda node: (node.attr["region"] == region
                                                    and node.attr["num_date"] >= pivots[0])

                tree_freqs = tree_frequencies(tree, pivots, method='SLSQP',
                                              node_filter = node_filter_func,
                                              ws = max(2, tree.count_terminals()//10),
                                              stiffness = stiffness, inertia=inertia,
                                              min_clades=args.minimal_clade_size_to_estimate)

                tree_freqs.estimate_clade_frequencies()

                frequency_dict["counts"][region] = [int(x) for x in tree_freqs.counts]
                if args.output_format == "nextflu":
                    # Export frequencies in nextflu-format by region and clade id.
                    for clade_id, clade_frequencies in tree_freqs.frequencies.items():
                        frequency_dict["%s_clade:%d" % (region, clade_id)] = format_frequencies(clade_frequencies)

                else:
                    # Export frequencies in auspice-format by strain name.
                    for node in tree.find_clades(order='postorder'):
                        if node.is_terminal():
                            node.tipcount=1
                        else:
                            node.tipcount = np.sum([c.tipcount for c in node])

                        if (node.is_terminal() or args.include_internal_nodes) and node.tipcount>args.minimal_clade_size:
                            if node.name not in frequency_dict:
                                frequency_dict[node.name] = {}
                            frequency_dict[node.name][region] = format_frequencies(tree_freqs.frequencies[node.clade])

        elif args.method == "kde":
            if args.output_format == "nextflu":
                print("ERROR: nextflu format is not supported for KDE frequencies", file=sys.stderr)
                return 1

            # Estimate frequencies.

            kde_frequencies = TreeKdeFrequencies(