simulator.py

# -*- coding: utf-8 -*-
import logging
import math as m
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import random
from scipy.special import comb
from scipy.stats import gamma

from wbia_lca import cluster_tools as ct
from wbia_lca import exp_scores as es
from wbia_lca import weighter as wgtr


logger = logging.getLogger('wbia_lca')


class simulator(object):  # NOQA
    """
    Important note (2019-08-01): when human-weighted edges are added, only the
    latest edge is currently "remembered" in this graph.  This could
    diverge from the actual graph when that graph combines weights
    additively.  One solution is to make this a multigraph, but not now.
    """

    def __init__(self, params, wgtr, seed=None):
        self.params = params
        self.wgtr = wgtr
        self.gt_clustering = dict()  # cid -> at first a list of nodes, then a set
        self.gt_node2cid = dict()  # node -> cid
        self.ranker_matches = dict()  # node -> set of nodes
        self.G = nx.Graph()
        self.G_orig = nx.Graph()  # the generated graph without algorithmic additions
        self.verify_edges = []
        self.human_edges = []
        if seed is not None:
            random.seed(seed)
        self.prev_num_human = 0
        self.max_delay_steps = 10
        self.human_steps_until_return = -1
        self.verify_steps_until_return = -1

    def generate(self):
        expected_nodes = 1 + self.params['gamma_shape'] * self.params['gamma_scale']

        digits_per_node = 2 + int(m.log10(expected_nodes))
        next_index = 0
        nodes = []  # list of node ids
        edges = []  # list of edge 3-tuples

        num_correct_positive = 0
        num_correct_negative = 0
        num_correct_zero = 0
        num_incorrect_positive = 0
        num_incorrect_negative = 0
        num_incorrect_zero = 0

        """
        Step 0:
        """
        samples = np.random.gamma(
            self.params['gamma_shape'],
            self.params['gamma_scale'],
            self.params['num_clusters'],
        )
        samples = 1 + np.round(samples)
        samples = samples.astype(int)

        """
        Step 1:
        Generate the clusters, the nodes within the cluster, and the
        "correct" inter-cluster edges.  Note that since we are
        assuming an imperfect ranking algorithm, this does not ensure
        that each cluster is connected.
        """
        num_from_ranker = self.params['num_from_ranker']
        cids = ct.cids_from_range(len(samples), prefix='ct')
        for i, cid in enumerate(cids):
            self.gt_clustering[cid] = list()

            n = samples[i]

            # Create the nodes in the cluster
            for i in range(n):
                node_id = 'n' + str(next_index).zfill(digits_per_node)
                next_index += 1
                nodes.append(node_id)
                self.gt_clustering[cid].append(node_id)
                self.gt_node2cid[node_id] = cid
                self.ranker_matches[node_id] = set()

            #  Create the positive edges between nodes in a cluster.
            #  These are symmetric. Don't allow more than num_from_ranker
            #  matches / edges for any node.
            for i, ith_node in enumerate(self.gt_clustering[cid]):
                for j in range(i + 1, len(self.gt_clustering[cid])):
                    prob = random.uniform(0, 1)
                    jth_node = self.gt_clustering[cid][j]
                    if (
                        prob < self.params['p_ranker_correct']
                        and len(self.ranker_matches[ith_node]) < num_from_ranker
                        and len(self.ranker_matches[jth_node]) < num_from_ranker
                    ):
                        self.ranker_matches[ith_node].add(jth_node)
                        self.ranker_matches[jth_node].add(ith_node)
                        is_match_correct = True
                        wgt = self.wgtr.random_wgt(is_match_correct)
                        if wgt > 0:
                            num_correct_positive += 1
                        elif wgt == 0:
                            num_correct_zero += 1
                        else:
                            num_correct_negative += 1

                        e = (ith_node, jth_node, wgt)
                        edges.append(e)

        assert num_from_ranker > 0
        num_nodes = len(nodes)

        # Change the list to a set
        self.gt_clustering = {
            cid: set(cluster) for cid, cluster in self.gt_clustering.items()
        }

        """
        Step 2:
        Generate "incorrect" match edges, sufficient to have the required
        number of edges generated by the ranking algorithm.
        """
        for i, ith_node in enumerate(nodes):
            matches = self.ranker_matches[ith_node]
            cid = self.gt_node2cid[ith_node]
            cluster = set(self.gt_clustering[cid])

            """
            Generate (incorrect) edges between clusters
            """
            is_match_correct = False
            while len(matches) < num_from_ranker:
                j = random.randint(0, num_nodes - 1)
                jth_node = nodes[j]
                if jth_node not in matches and jth_node not in cluster:
                    matches.add(jth_node)
                    wgt = self.wgtr.random_wgt(is_match_correct)
                    if wgt > 0:
                        num_incorrect_positive += 1
                    elif wgt == 0:
                        num_incorrect_zero += 1
                    else:
                        num_incorrect_negative += 1

                    if ith_node < jth_node:
                        e = (ith_node, jth_node, wgt)
                    else:
                        e = (jth_node, ith_node, wgt)
                    edges.append(e)

        self.G.add_weighted_edges_from(edges)
        logger.info('simulator::generate: adding %d edges' % len(edges))
        logger.info('%d correct match edges have positive weight' % num_correct_positive)
        logger.info('%d correct match edges have zero weight' % num_correct_zero)
        logger.info('%d correct match edges have negative weight' % num_correct_negative)
        logger.info(
            '%d incorrect match edges have positive weight' % num_incorrect_positive
        )
        logger.info('%d incorrect match edges have zero weight' % num_incorrect_zero)
        logger.info(
            '%d incorrect match edges have negative weight' % num_incorrect_negative
        )

        self.G_orig.add_nodes_from(self.G)
        self.G_orig.add_weighted_edges_from(edges)

        """
        Step 3: Generate the "reachable" ground truth, the obtainable
        result given simulated failures to match that could disconnect
        a correct match.
        """
        self.r_clustering = dict()
        k = 0
        for cc in self.gt_clustering.values():
            H = self.G.subgraph(cc)
            prev_k = k
            for new_cc in nx.connected_components(H):
                self.r_clustering[k] = new_cc
                k += 1
            if k - prev_k > 1:
                logger.info('GT cluster %a split into %a ...' % (cc, k - prev_k))
                for i in range(prev_k, k):
                    logger.info('   %a' % self.r_clustering[i])
            else:
                logger.info('GT cluster %a is intact' % cc)
        self.r_node2cid = ct.build_node_to_cluster_mapping(self.r_clustering)

        """
        Step 4: Reconfigure edges to maks the expected input to the
        graph algorithm weight manager.
        """
        aug_names = ['verifier', 'human']
        edges = [(n0, n1, w, aug_names[0]) for n0, n1, w in edges]

        return edges, aug_names

    def print_clusters(self):
        logger.info('Ground truth clusters:')
        for cid in self.gt_clustering:
            logger.info('%a: %a' % (cid, self.gt_clustering[cid]))

    def save_graph(self):
        assert False

    def restore_graph(self):
        assert False

    def augmentation_request(self, triples):
        pairs_for_verify = []
        pairs_for_human = []
        for n0, n1, a_name in triples:
            if a_name == 'verifier':
                pairs_for_verify.append((n0, n1))
            elif a_name == 'human':
                pairs_for_human.append((n0, n1))
            else:
                assert False
        self.verify_request(pairs_for_verify)
        self.human_request(pairs_for_human)

    def augmentation_result(self, node_set=None):
        """
        Ignore the node_set...
        """
        v_results = [(n0, n1, w, 'verifier') for n0, n1, w in self.verify_result()]
        h_results = [(n0, n1, w, 'human') for n0, n1, w in self.human_result()]
        return v_results + h_results

    def verify_request(self, node_pairs):
        """
        The result of verification is deterministic and symmetric.
        So we need to keep and return the same result if called multiple
        times. To do this, we keep it in the graph.  By contrast, the
        result of each human decision is generated randomly for each
        request.
        """
        if self.verify_steps_until_return < 0:
            self.verify_steps_until_return = random.randint(0, self.max_delay_steps)
            """
            logger.info(
                'sim::verify_request: first delay will be %d steps'
                % self.verify_steps_until_return
            )
            """

        new_edges = []
        for pr in set(node_pairs):  # make sure no repeats
            if pr[1] not in self.G[pr[0]]:  # no edge already
                assert pr[0] < pr[1]
                cid0 = self.gt_node2cid[pr[0]]
                cid1 = self.gt_node2cid[pr[1]]
                is_match_correct = cid0 == cid1
                wgt = self.wgtr.random_wgt(is_match_correct)
                e = tuple([pr[0], pr[1], wgt])
                new_edges.append(e)
                self.verify_edges.append(e)

        self.G.add_weighted_edges_from(new_edges)

    def verify_result(self):
        if self.verify_steps_until_return == 0:
            edges_returned = self.verify_edges.copy()
            self.verify_edges.clear()
            self.verify_steps_until_return = random.randint(0, self.max_delay_steps)
            logger.info(
                'sim::verify_request: next delay will be %d steps'
                % self.verify_steps_until_return
            )
            return edges_returned
        else:
            self.verify_steps_until_return -= 1
            return []

    def gen_human_wgt(self, pr):
        cid0 = self.gt_node2cid[pr[0]]
        cid1 = self.gt_node2cid[pr[1]]
        is_match_correct = cid0 == cid1
        wgt = self.wgtr.human_random_wgt(is_match_correct)
        return wgt

    def human_request(self, node_pairs):
        if self.human_steps_until_return < 0:
            self.human_steps_until_return = random.randint(0, self.max_delay_steps)
            logger.info(
                'sim::human_request: delay will be %d steps'
                % self.human_steps_until_return
            )

        new_edges = []
        for pr in set(node_pairs):
            if pr[0] > pr[1]:
                pr = (pr[1], pr[0])
            wgt = self.gen_human_wgt(pr)
            e = tuple([pr[0], pr[1], wgt])
            new_edges.append(e)
            self.human_edges.append(e)
        self.G.add_weighted_edges_from(new_edges)

    def human_result(self):
        if self.human_steps_until_return == 0:
            edges_returned = self.human_edges.copy()
            self.human_edges.clear()
            self.human_steps_until_return = random.randint(0, self.max_delay_steps)
            logger.info(
                'sim::human_result: next delay will be %d steps'
                % self.human_steps_until_return
            )
            return edges_returned
        else:
            self.human_steps_until_return -= 1
            return []

    @staticmethod
    def incremental_stats(
        num_human, clustering, node2cid, true_clustering, true_node2cid
    ):
        frac, prec, rec = ct.percent_and_PR(
            clustering, node2cid, true_clustering, true_node2cid
        )
        result = {
            'num human': num_human,
            'num clusters': len(clustering),
            'num true clusters': len(true_clustering),
            'frac correct': frac,
            'precision': prec,
            'recall': rec,
        }
        return result

    def trace_start_human(self, clustering, node2cid):
        """
        Beging to record information about the number of human decisions
        vs. the accuracy of the current clustering.  The comparison is
        made against both the ground truth clustering and the "reachable"
        ground truth clustering. For each new number of
        human decisions, we record (1) this number, (2) the number of
        ground truth clusters (fixed value), (3) the number of current
        clusters, (4) the fraction of current clusters that are
        exactly correct, (5) the precision and (6) the recall.  The
        same thing will be done for the "reachable" clusters.
        """
        result = self.incremental_stats(
            0, clustering, node2cid, self.gt_clustering, self.gt_node2cid
        )
        self.gt_results = [result]
        result = self.incremental_stats(
            0, clustering, node2cid, self.r_clustering, self.r_node2cid
        )
        self.r_results = [result]
        self.prev_num_human = 0

    def trace_iter_compare_to_gt(self, clustering, node2cid, num_human):
        if num_human <= self.prev_num_human:
            return
        result = self.incremental_stats(
            num_human, clustering, node2cid, self.gt_clustering, self.gt_node2cid
        )
        self.gt_results.append(result)
        result = self.incremental_stats(
            num_human, clustering, node2cid, self.r_clustering, self.r_node2cid
        )
        self.r_results.append(result)
        self.prev_num_human = num_human

    @staticmethod
    def csv_output(out_file, results):
        with open(out_file, 'w') as f:
            f.write(
                'human decisions, num clusters, num true clusters, '
                + 'frac eq true, precision, recall\n'
            )
            for res in results:
                f.write(
                    '%d, %d, %d, %.4f, %.4f %.4f\n'
                    % (
                        res['num human'],
                        res['num clusters'],
                        res['num true clusters'],
                        res['frac correct'],
                        res['precision'],
                        res['recall'],
                    )
                )

    @staticmethod
    def plot_convergence(results, filename=None):
        num_human_decisions = [r['num human'] for r in results]
        num_actual_clusters = [r['num clusters'] for r in results]
        num_true_clusters = [r['num true clusters'] for r in results]
        # frac_correct = [r['frac correct'] for r in results]

        fig, ax1 = plt.subplots()

        ax1.set_xlabel('Number of human decisions')
        ax1.set_ylabel('Number of clusters')

        color = 'blue'
        ax1.plot(num_human_decisions, num_actual_clusters, color=color)

        # ADD FRAC CORRECT TO PLOT ON RIGHT SIDE!!!!

        color = 'green'
        ax1.plot(num_human_decisions, num_true_clusters, color=color)

        if filename is None:
            plt.show()
        else:
            plt.savefig(filename)
            logger.info('Saved plot to %s' % filename)
            plt.close()

    def generate_plots(self, out_prefix):
        out_name = out_prefix + '_gt.pdf'
        self.plot_convergence(self.gt_results, out_name)
        out_name = out_prefix + '_reach_gt.pdf'
        self.plot_convergence(self.r_results, out_name)


def find_np_ratio(gamma_shape, gamma_scale, ranker_per_node, prob_match):
    logger.info('=====')
    logger.info(
        'gamma_shape %a gamma_scale %a ranker_per_node %d, prob_match %1.3f'
        % (gamma_shape, gamma_scale, ranker_per_node, prob_match)
    )
    logger.info('from these we get (requiring at least one per node)')
    mean = 1 + gamma_shape * gamma_scale
    mode = 1 + (gamma_shape - 1) * gamma_scale
    std_dev = m.sqrt(gamma_shape) * gamma_scale

    logger.info(
        'mean %s mode %s std_dev %.2f'
        % (
            mean,
            mode,
            std_dev,
        )
    )
    limit = int(mean + 10 * std_dev)

    pos_per_node = 0
    for k in range(2, limit + 1):
        prob_k = gamma.cdf(k - 0.5, gamma_shape, scale=gamma_scale) - gamma.cdf(
            k - 1.5, gamma_shape, scale=gamma_scale
        )

        max_correct = min(ranker_per_node, k - 1)
        exp_correct = 0
        sum_prob_i = 0
        for i in range(max_correct + 1):
            prob_i = comb(k - 1, i) * prob_match ** i * (1 - prob_match) ** (k - 1 - i)
            term = i * prob_i
            sum_prob_i += prob_i
            exp_correct += term
        exp_correct /= sum_prob_i

        pos_per_node += prob_k * exp_correct

    neg_per_node = ranker_per_node - pos_per_node
    ratio = neg_per_node / pos_per_node
    logger.info(
        'neg_per_node %1.2f,  pos_per_node %1.2f, np_ratio %1.2f'
        % (neg_per_node, pos_per_node, ratio)
    )
    return ratio


def test_find_np_ratio():
    gamma_shape = 2
    gamma_scale = 3
    ranker_per_node = 10
    prob_match = 0.85
    ratio = find_np_ratio(gamma_shape, gamma_scale, ranker_per_node, prob_match)
    logger.info('Final np ratio is %1.4f' % ratio)

    gamma_shape = 1
    gamma_scale = 3
    ranker_per_node = 15
    prob_match = 0.85
    ratio = find_np_ratio(gamma_shape, gamma_scale, ranker_per_node, prob_match)
    logger.info('Final np ratio is %1.4f' % ratio)


def ensure_after_gen(sim, params):
    logger.info(
        '\n================================\n'
        'Testing the generated simulator:\n'
        '================================'
    )
    logger.info(
        'Number of clusters:  should be %d is %d'
        % (sim.params['num_clusters'], len(sim.gt_clustering))
    )

    num_nodes_in_graph = len(sim.G.nodes())
    num_nodes_in_node2cid = len(sim.gt_node2cid)
    num_nodes_in_clusters = sum([len(c) for c in sim.gt_clustering.values()])

    logger.info(
        'All three numbers of nodes should be equal:\n '
        + '    %d in graph\n' % (num_nodes_in_graph,)
        + '    %d in node2cid\n' % (num_nodes_in_node2cid,)
        + '    %d in cluster' % (num_nodes_in_clusters,)
    )

    """
    All node_matches sets are the same length, and that length is the params values.
    """
    num_from_ranker = params['num_from_ranker']
    logger.info('Each node should have %d matches from the ranker' % num_from_ranker)
    all_same = True
    for node, matches in sim.ranker_matches.items():
        if len(matches) != num_from_ranker:
            all_same = False
            logger.info('Error: node %s has %d matches' % (node, len(matches)))
    if all_same:
        logger.info('Yes. All have the correct number of ranker matches')

    """
    Examine the distribution of true and false edges
    """
    tp = fp = tn = fn = zero_pos = zero_neg = 0
    for n0 in sim.G.nodes():
        cid0 = sim.gt_node2cid[n0]
        for n1 in sim.G[n0]:
            if n0 < n1:
                cid1 = sim.gt_node2cid[n1]
                wgt = sim.G[n0][n1]['weight']
                if cid0 == cid1:
                    if wgt > 0:
                        tp += 1
                    elif wgt < 0:
                        fn += 1
                    else:
                        zero_pos += 1
                else:
                    if wgt > 0:
                        fp += 1
                    elif wgt < 0:
                        tn += 1
                    else:
                        zero_neg += 1

    logger.info(
        'tp %d    fp %d\nfn %d    tn %d\nzero_pos %d  zero_neg %d'
        % (tp, fp, fn, tn, zero_pos, zero_neg)
    )
    logger.info('actual ratio = %1.2f' % ((fp + tn + zero_neg) / (tp + fn + zero_pos)))
    logger.info('actual number of xx positive pairs %s' % (tp + fn + zero_pos,))
    logger.info('actual number of xx negative pairs %s' % (tn + fp + zero_neg,))
    """
    exp_neg, exp_pos = sim.wgtr.eval_overlap()
    logger.info("Correct match fraction with positive weight = %5.3f, expected = %5.3f"
          % ((tp / (tp+fn)), exp_pos))
    logger.info("Incorrect match fraction with negative weight = %5.3f, expected = %5.3f"
          % ((tn / (tn+fp)), exp_neg))
    """

    """
    Did we get the right fraction of correct matches becoming real matches?
    """
    num_correct_edges = tp + fn + zero_pos
    num_possible = 0
    for c in sim.gt_clustering.values():
        nc = len(c)
        num_possible += nc * (nc - 1) // 2
    logger.info(
        'Fraction of true edges generated = %1.3f, expected = %1.3f'
        % (num_correct_edges / num_possible, sim.params['p_ranker_correct'])
    )

    avg_per_cluster = len(sim.G.nodes()) / len(sim.gt_clustering)
    exp_per_cluster = 1 + sim.params['gamma_shape'] * sim.params['gamma_scale']
    logger.info('Average cluster size is %1.2f' % avg_per_cluster)
    logger.info('Expected cluster size is %1.2f' % exp_per_cluster)


def get_positive_missing(sim, num_requested):
    pos_missing = []
    nodes = list(sim.gt_node2cid.keys())
    random.shuffle(nodes)
    for n0 in nodes:
        cid0 = sim.gt_node2cid[n0]
        for n1 in sim.gt_clustering[cid0]:
            if n0 != n1 and n1 not in sim.G[n0]:
                pr = (min(n0, n1), max(n0, n1))
                pos_missing.append(pr)
                break
        if len(pos_missing) == num_requested:
            break
    return pos_missing


def get_negative_missing(sim, num_requested):
    neg_missing = []
    max_tries = 50  # avoid infinite loop
    tries = 0
    nodes = list(sim.gt_node2cid.keys())
    num_nodes = len(nodes)
    while len(neg_missing) < num_requested and tries < max_tries:
        tries += 1
        n0 = nodes[random.randint(0, num_nodes - 1)]
        n1 = nodes[random.randint(0, num_nodes - 1)]
        cid0 = sim.gt_node2cid[n0]
        cid1 = sim.gt_node2cid[n1]
        if n0 != n1 and cid0 != cid1 and n1 not in sim.G[n0]:
            pr = (min(n0, n1), max(n0, n1))
            neg_missing.append(pr)
    return neg_missing


def ensure_verify(sim, params):
    """
    1. Test abilty to ignore edges already in the graph.
    """
    pairs = []
    nodes = list(sim.gt_node2cid.keys())
    num_nodes = len(nodes)
    num_already_there = 2
    for _ in range(num_already_there):

        n0 = nodes[random.randint(0, num_nodes - 1)]
        n1 = next(iter(sim.G[n0].keys()))  # first nbr in n1 dictionary
        pr = (min(n0, n1), max(n0, n1))
        pairs.append(pr)
    sim.verify_request(pairs)
    edges = sim.verify_result()
    logger.info('test_verify:')
    logger.info(
        'After adding edges already in the graph, num returned should be 0. It is %s'
        % (len(edges),)
    )

    """
    2. Test ability to gather missing edges.  Do it 2x
    """
    for tries in range(2):
        # Get missing edges
        pos_missing = get_positive_missing(sim, tries + 1)
        neg_missing = get_negative_missing(sim, tries + 2)
        missing = pos_missing + neg_missing
        logger.info('Verify request / results try %s' % (tries,))
        logger.info('Num pos_missing (should be 3) is %s' % (len(pos_missing),))
        logger.info('Num neg_missing (should be 2) is %s' % (len(neg_missing),))

        #  Make the request for weights for these edges. This does not get them
        #  yet.  Make sure though that it has set the delay.
        sim.verify_request(missing)
        steps_until = sim.verify_steps_until_return
        logger.info(
            'Delay steps until return should be non-negative and is %s' % (steps_until,)
        )

        # For the given number of steps, the result should be empty
        has_error = False
        while steps_until > 0 and not has_error:
            steps_until -= 1
            edges = sim.verify_result()
            if len(edges) > 0:
                has_error = True
                logger.info(
                    'Error: with %s iterations remaining returned %s edges.'
                    % (
                        steps_until,
                        len(edges),
                    )
                )

        if not has_error:
            logger.info(
                'Successfully iterated through the delay without returning edges.'
            )

        #  This time it should return edges.
        edges = sim.verify_result()
        logger.info(
            'Requested missing edges: positive %s negative: %s'
            % (
                pos_missing,
                neg_missing,
            )
        )
        logger.info('returned are %s' % (edges,))
        if len(missing) != len(edges):
            logger.info('Error: incorrect number returned')
        else:
            logger.info('Correct number returned')
        logger.info(
            'sim.verify_edges should be len 0, and it is %s' % (len(sim.verify_edges),)
        )


def get_positive_for_human(sim, num_pos):
    """
    get several that are known positive and several that are known negative
    """
    pos = []
    cluster_ids = list(sim.gt_clustering.keys())
    random.shuffle(cluster_ids)
    for cid in cluster_ids:
        c = sim.gt_clustering[cid]
        if len(c) == 1:
            continue
        c = list(c)
        pr = (min(c[0], c[1]), max(c[0], c[1]))
        pos.append(pr)
        if len(pos) >= num_pos:
            break
    return pos


def get_negative_for_human(sim, num_neg):
    nodes = list(sim.gt_node2cid.keys())
    num_nodes = len(nodes)
    neg = []
    max_tries = 30  # avoid infinite loop
    tries = 0
    while len(neg) < num_neg and tries < max_tries:
        tries += 1
        n0 = nodes[random.randint(0, num_nodes - 1)]
        n1 = nodes[random.randint(0, num_nodes - 1)]
        cid0 = sim.gt_node2cid[n0]
        cid1 = sim.gt_node2cid[n1]
        if n0 != n1 and cid0 != cid1:
            pr = (min(n0, n1), max(n0, n1))
            neg.append(pr)
    return neg


def ensure_human(sim, params):
    logger.info('Test simulation of human')
    for tries in range(2):
        logger.info('try %s' % (tries,))
        logger.info('Test human request / result try %s' % (tries,))
        pos_for_human = get_positive_for_human(sim, tries + 2)
        neg_for_human = get_negative_for_human(sim, tries + 1)
        for_human = pos_for_human + neg_for_human
        logger.info(for_human)
        sim.human_request(for_human)

        logger.info(
            'Requested edges from human: positive %s negative: %s'
            % (
                pos_for_human,
                neg_for_human,
            )
        )
        steps_until = sim.human_steps_until_return
        logger.info(
            'Delay steps until return should be non-negative and is %s' % (steps_until,)
        )

        # For the given number of steps, the result should be empty
        has_error = False
        while steps_until > 0 and not has_error:
            steps_until -= 1
            edges = sim.human_result()
            if len(edges) > 0:
                has_error = True
                logger.info(
                    'Error: with %s iterations remaining returned %s edges.'
                    % (
                        steps_until,
                        len(edges),
                    )
                )

        if not has_error:
            logger.info(
                'Successfully iterated through the delay without returning edges.'
            )

        #  This time it should return edges.
        edges = sim.human_result()
        logger.info(
            'Requested human edges: positive %s negative: %s'
            % (
                pos_for_human,
                neg_for_human,
            )
        )
        logger.info('returned are %s' % (edges,))
        if len(for_human) != len(edges):
            logger.info('Error: incorrect number returned')
        else:
            logger.info('Correct number returned')
        logger.info(
            'sim.human_edges should be len 0, and it is %s' % (len(sim.human_edges),)
        )


def test_simulator():
    test_find_np_ratio()

    params = dict()
    params['pos_error_frac'] = 0.1
    params['num_clusters'] = 16

    # The following are parameters of the gamma distribution.
    # Recall the following properties:
    #      mean is shape*scale,
    #      mode is (shape-1)*scale
    #      variance is shape*scale**2 = mean*scale
    # So when these are both 2 the mean is 4, the mode is 2
    # and the variance is 4 (st dev 2).  And, when shape = 1,
    # we always have the mode at 0.
    #
    # The mean and the mode must be offset by 1 because every cluster
    # has at least one node.
    #
    params['gamma_shape'] = 1.5
    params['gamma_scale'] = 3
    # num_per_cluster = params['gamma_scale'] * params['gamma_shape'] + 1

    params['p_ranker_correct'] = 0.9
    params['p_human_correct'] = 0.98
    params['num_from_ranker'] = 10

    np_ratio = find_np_ratio(
        params['gamma_shape'],
        params['gamma_scale'],
        params['num_from_ranker'],
        params['p_ranker_correct'],
    )
    logger.info('np_ratio = %1.3f' % np_ratio)

    scorer = es.exp_scores.create_from_error_frac(params['pos_error_frac'], np_ratio)
    wgtr_ = wgtr.weighter(scorer, human_prob=params['p_human_correct'])

    sim = simulator(params, wgtr_)
    sim.generate()

    r_leng = len(sim.r_clustering)
    gt_leng = len(sim.gt_clustering)
    logger.info(
        'gt length %s reachable length %s'
        % (
            gt_leng,
            r_leng,
        )
    )

    """
    logger.info("\nClusters")
    for cid in sim.gt_clustering:
        logger.info("%3d: %a" % (cid, sim.gt_clusters[cid]))

    logger.info("Length gt_node2cid %d, should be equal to graph nodes %d" \
          % (len(sim.gt_node2cid), len(sim.G.nodes())))

    logger.info("\nRanker matches\n")
    for nid in sim.ranker_matches:
        logger.info("Node %s has matches:" % (nid, ), end=' ')
        for m in sim.ranker_matches[nid]:
            logger.info(m, end=' ')
        logger.info()

    logger.info("\nEdges for each node (ranker matches in each direction!)")
    for nid in sorted(sim.G.nodes()):
        logger.info("Node ", nid)
        for m in sim.G[nid]:
            n0, n1 = min(nid, m), max(nid, m)
            logger.info("(%a, %a, %1.4f)" % (n0, n1, sim.G[nid][m]['weight']))
    """

    ensure_after_gen(sim, params)
    ensure_verify(sim, params)
    ensure_human(sim, params)


if __name__ == '__main__':
    test_simulator()