text:
from typing import Set, Dict, Any
from metagraph import ConcreteType, dtypes
from ..core.types import Matrix, EdgeSet, EdgeMap, Graph
from ..core.wrappers import EdgeSetWrapper, EdgeMapWrapper, GraphWrapper
from .. import has_scipy
import numpy as np

if has_scipy:
    import scipy.sparse as ss

    class ScipyEdgeSet(EdgeSetWrapper, abstract=EdgeSet):
        """
        scipy.sparse matrix is the minimal size to contain all edges.
        If nodes are not sequential, a node_list must be provided to map the matrix index to NodeId.
        Nodes which are present in the matrix but have no edges are not allowed as
            they will not survive a roundtrip translation.
        The actual values in the matrix are unused.
        """

        def __init__(self, data, node_list=None, *, aprops=None):
            super().__init__(aprops=aprops)
            self._assert_instance(data, ss.spmatrix)
            nrows, ncols = data.shape
            self._assert(nrows == ncols, "Adjacency Matrix must be square")
            self.value = data
            if node_list is None:
                node_list = np.arange(nrows)
            else:
                self._assert_instance(node_list, (np.ndarray, list, tuple))
                if not isinstance(node_list, np.ndarray):
                    node_list = np.array(node_list)
            self._assert(
                nrows == len(node_list),
                f"node list size ({len(node_list)}) and data matrix size ({nrows}) don't match.",
            )
            self.node_list = node_list

        # def copy(self):
        #     return ScipyEdgeSet(self.value.copy(), node_list=self.node_list.copy())

        class TypeMixin:
            @classmethod
            def _compute_abstract_properties(
                cls, obj, props: Set[str], known_props: Dict[str, Any]
            ) -> Dict[str, Any]:
                ret = known_props.copy()

                # slow properties, only compute if asked
                for prop in props - ret.keys():
                    if prop == "is_directed":
                        ret[prop] = (obj.value.T != obj.value).nnz > 0

                return ret

            @classmethod
            def assert_equal(
                cls,
                obj1,
                obj2,
                aprops1,
                aprops2,
                cprops1,
                cprops2,
                *,
                rel_tol=None,
                abs_tol=None,
            ):
                m1, m2 = obj1.value, obj2.value
                assert (
                    m1.shape[0] == m2.shape[0]
                ), f"size mismatch: {m1.shape[0]} != {m2.shape[0]}"
                assert m1.nnz == m2.nnz, f"num edges mismatch: {m1.nnz} != {m2.nnz}"
                if not (obj1.node_list == obj2.node_list).all():
                    nl1 = set(obj1.node_list)
                    nl2 = set(obj2.node_list)
                    assert nl1 == nl2, f"node list mismatch: {nl1 ^ nl2}"
                assert aprops1 == aprops2, f"property mismatch: {aprops1} != {aprops2}"
                # Convert to COO format, apply node ids, then sort to allow comparison of indices
                d1 = m1.tocoo()
                d2 = m2.tocoo()
                r1, c1 = obj1.node_list[d1.row], obj1.node_list[d1.col]
                r2, c2 = obj2.node_list[d2.row], obj2.node_list[d2.col]
                sort1 = np.lexsort((c1, r1))
                sort2 = np.lexsort((c2, r2))
                r1, c1 = r1[sort1], c1[sort1]
                r2, c2 = r2[sort2], c2[sort2]
                assert (r1 == r2).all(), f"rows mismatch {r1} != {r2}"
                assert (c1 == c2).all(), f"cols mismatch {c1} != {c2}"

    class ScipyEdgeMap(EdgeMapWrapper, abstract=EdgeMap):
        """
        scipy.sparse matrix is the minimal size to contain all edges.
        If nodes are not sequential, a node_list must be provided to map the matrix index to NodeId.
        Nodes which are present in the matrix but have no edges are not allowed as
            they will not survive a roundtrip translation.
        """

        def __init__(self, data, node_list=None, *, aprops=None):
            super().__init__(aprops=aprops)
            self._assert_instance(data, ss.spmatrix)
            nrows, ncols = data.shape
            self._assert(nrows == ncols, "Adjacency Matrix must be square")
            self.value = data
            if node_list is None:
                node_list = np.arange(nrows)
            else:
                self._assert_instance(node_list, (np.ndarray, list, tuple))
                if not isinstance(node_list, np.ndarray):
                    node_list = np.array(node_list)
            self._assert(
                nrows == len(node_list),
                f"node list size ({len(node_list)}) and data matrix size ({nrows}) don't match.",
            )
            self.node_list = node_list

        # def copy(self):
        #     node_list = (
        #         self.node_list if self.node_list is None else self.node_list.copy()
        #     )
        #     return ScipyEdgeMap(self.value.copy(), node_list=node_list)

        class TypeMixin:
            @classmethod
            def _compute_abstract_properties(
                cls, obj, props: Set[str], known_props: Dict[str, Any]
            ) -> Dict[str, Any]:
                ret = known_props.copy()

                # fast properties
                for prop in {"dtype"} - ret.keys():
                    if prop == "dtype":
                        ret[prop] = dtypes.dtypes_simplified[obj.value.dtype]

                # slow properties, only compute if asked
                for prop in props - ret.keys():
                    if prop == "is_directed":
                        ret[prop] = (obj.value.T != obj.value).nnz > 0
                    if prop == "has_negative_weights":
                        if ret["dtype"] in {"bool", "str"}:
                            neg_weights = None
                        else:
                            min_val = obj.value.data.min()
                            if min_val < 0:
                                neg_weights = True
                            else:
                                neg_weights = False
                        ret[prop] = neg_weights

                return ret

            @classmethod
            def assert_equal(
                cls,
                obj1,
                obj2,
                aprops1,
                aprops2,
                cprops1,
                cprops2,
                *,
                rel_tol=1e-9,
                abs_tol=0.0,
            ):
                m1, m2 = obj1.value, obj2.value
                assert (
                    m1.shape[0] == m2.shape[0]
                ), f"size mismatch: {m1.shape[0]} != {m2.shape[0]}"
                assert m1.nnz == m2.nnz, f"num edges mismatch: {m1.nnz} != {m2.nnz}"
                if not (obj1.node_list == obj2.node_list).all():
                    nl1 = set(obj1.node_list)
                    nl2 = set(obj2.node_list)
                    assert nl1 == nl2, f"node list mismatch: {nl1 ^ nl2}"
                assert aprops1 == aprops2, f"property mismatch: {aprops1} != {aprops2}"
                # Convert to COO format, apply node ids, then sort to allow comparison of indices and values
                d1 = m1.tocoo()
                d2 = m2.tocoo()
                r1, c1, v1 = obj1.node_list[d1.row], obj1.node_list[d1.col], d1.data
                r2, c2, v2 = obj2.node_list[d2.row], obj2.node_list[d2.col], d2.data
                sort1 = np.lexsort((c1, r1))
                sort2 = np.lexsort((c2, r2))
                r1, c1, v1 = r1[sort1], c1[sort1], v1[sort1]
                r2, c2, v2 = r2[sort2], c2[sort2], v2[sort2]
                assert (r1 == r2).all(), f"rows mismatch {r1} != {r2}"
                assert (c1 == c2).all(), f"cols mismatch {c1} != {c2}"
                if issubclass(d1.dtype.type, np.floating):
                    assert np.isclose(v1, v2, rtol=rel_tol, atol=abs_tol).all()
                else:
                    assert (v1 == v2).all()

    class ScipyGraph(GraphWrapper, abstract=Graph):
        """
        scipy.sparse matrix is the minimal size to contain all nodes in the graph.
        If nodes are not sequential, a node_list must be provided to map the matrix index to NodeId.
        node_vals (if populated) contains node weights
        """

        def __init__(self, matrix, node_list=None, node_vals=None, *, aprops=None):
            super().__init__(aprops=aprops)
            self._assert_instance(matrix, ss.spmatrix)
            nrows, ncols = matrix.shape
            self._assert(
                nrows == ncols, f"adjacency matrix must be square, not {nrows}x{ncols}"
            )
            if node_list is None:
                node_list = np.arange(nrows)
            else:
                self._assert_instance(node_list, (np.ndarray, list, tuple))
                if not isinstance(node_list, np.ndarray):
                    node_list = np.array(node_list)
            self._assert(
                nrows == len(node_list),
                f"node list size ({len(node_list)}) and data matrix size ({nrows}) don't match.",
            )
            if node_vals is not None:
                self._assert_instance(node_vals, (np.ndarray, list, tuple))
                if not isinstance(node_vals, np.ndarray):
                    node_vals = np.array(node_vals)
                self._assert(
                    nrows == len(node_vals),
                    f"node vals size ({len(node_vals)}) and data matrix size ({nrows}) don't match",
                )
            self.value = matrix
            self.node_list: np.ndarray = node_list
            self.node_vals: np.ndarray = node_vals

        # def copy(self):
        #     node_vals = None if self.node_vals is None else self.node_vals.copy()
        #     return ScipyGraph(self.value.copy(), self.node_list.copy(), node_vals)

        class TypeMixin:
            # Both forward and reverse lookup
            _edge_prop_map = {
                "is_directed": "is_directed",
                "edge_dtype": "dtype",
                "edge_has_negative_weights": "has_negative_weights",
                "dtype": "edge_dtype",
                "has_negative_weights": "edge_has_negative_weights",
            }

            @classmethod
            def _compute_abstract_properties(
                cls, obj, props: Set[str], known_props: Dict[str, Any]
            ) -> Dict[str, Any]:
                ret = known_props.copy()

                # fast properties
                for prop in {"node_type", "edge_type"} - ret.keys():
                    if prop == "node_type":
                        ret[prop] = "set" if obj.node_vals is None else "map"
                    elif prop == "edge_type":
                        ret[prop] = "set" if obj.value.dtype == bool else "map"

                # Delegate to ScipyEdge{Set/Map} to compute edge properties
                if ret["edge_type"] == "set":
                    ret["edge_dtype"] = None
                    ret["edge_has_negative_weights"] = None
                    edgeclass = ScipyEdgeSet
                else:
                    edgeclass = ScipyEdgeMap
                edge_props = {
                    cls._edge_prop_map[p] for p in props if p in cls._edge_prop_map
                }
                known_edge_props = {
                    cls._edge_prop_map[p]: v
                    for p, v in known_props.items()
                    if p in cls._edge_prop_map
                }
                edges = edgeclass(obj.value, obj.node_list)
                edge_computed_props = edgeclass.Type._compute_abstract_properties(
                    edges, edge_props, known_edge_props
                )
                ret.update(
                    {cls._edge_prop_map[p]: v for p, v in edge_computed_props.items()}
                )

                # slow properties, only compute if asked
                for prop in props - ret.keys():
                    if prop == "node_dtype":
                        if ret["node_type"] == "set":
                            ret[prop] = None
                        else:
                            ret[prop] = dtypes.dtypes_simplified[obj.node_vals.dtype]

                return ret

            @classmethod
            def assert_equal(
                cls,
                obj1,
                obj2,
                aprops1,
                aprops2,
                cprops1,
                cprops2,
                *,
                rel_tol=1e-9,
                abs_tol=0.0,
            ):
                assert aprops1 == aprops2, f"property mismatch: {aprops1} != {aprops2}"
                subprops1 = {
                    cls._edge_prop_map[p]: v
                    for p, v in aprops1.items()
                    if p in cls._edge_prop_map
                }
                subprops2 = {
                    cls._edge_prop_map[p]: v
                    for p, v in aprops2.items()
                    if p in cls._edge_prop_map
                }
                if aprops1["edge_type"] == "set":
                    edgeset1 = ScipyEdgeSet(obj1.value, obj1.node_list)
                    edgeset2 = ScipyEdgeSet(obj2.value, obj2.node_list)
                    ScipyEdgeSet.Type.assert_equal(
                        edgeset1,
                        edgeset2,
                        subprops1,
                        subprops2,
                        {},
                        {},
                        rel_tol=rel_tol,
                        abs_tol=abs_tol,
                    )
                else:
                    edgemap1 = ScipyEdgeMap(obj1.value, obj1.node_list)
                    edgemap2 = ScipyEdgeMap(obj2.value, obj2.node_list)
                    ScipyEdgeMap.Type.assert_equal(
                        edgemap1,
                        edgemap2,
                        subprops1,
                        subprops2,
                        {},
                        {},
                        rel_tol=rel_tol,
                        abs_tol=abs_tol,
                    )
                if aprops1["node_type"] == "map":
                    sort1 = np.argsort(obj1.node_list)
                    sort2 = np.argsort(obj2.node_list)
                    vals1 = obj1.node_vals[sort1]
                    vals2 = obj2.node_vals[sort2]
                    if issubclass(vals1.dtype.type, np.floating):
                        assert np.isclose(
                            vals1, vals2, rtol=rel_tol, atol=abs_tol
                        ).all()
                    else:
                        assert (vals1 == vals2).all()

max_stars_repo_name: paul-tqh-nguyen/metagraph

max_issues_repo_name: paul-tqh-nguyen/metagraph

max_forks_repo_name: paul-tqh-nguyen/metagraph

reason: import numpy,import scipy


text:
"""
Network Plot
===============================================================================

>>> from techminer2 import *
>>> directory = "/workspaces/techminer2/data/"
>>> file_name = "/workspaces/techminer2/sphinx/images/network_plot.png"
>>> coc_matrix = co_occurrence_matrix(
...     column='author_keywords',
...     min_occ=7,
...     directory=directory,
... )
>>> from techminer2.network_api.network import network
>>> network_ = network(coc_matrix)
>>> from techminer2.network_api.network_plot import network_plot
>>> network_plot(network_).savefig(file_name)

.. image:: images/network_plot.png
    :width: 700px
    :align: center

"""
from operator import itemgetter

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np

group_colors = [
    "tab:blue",
    "tab:orange",
    "tab:green",
    "tab:red",
    "tab:purple",
    "tab:brown",
    "tab:pink",
    "tab:gray",
    "tab:olive",
    "tab:cyan",
    "cornflowerblue",
    "lightsalmon",
    "limegreen",
    "tomato",
    "mediumvioletred",
    "darkgoldenrod",
    "lightcoral",
    "silver",
    "darkkhaki",
    "skyblue",
] * 5

def network_plot(
    network,
    figsize=(7, 7),
    k=0.20,
    iterations=50,
    max_labels=50,
):

    fig = plt.figure(figsize=figsize)
    ax = fig.subplots()

    options = {
        "width": 1,
        "with_labels": False,
        "font_size": 7,
        "font_weight": "regular",
        "alpha": 0.7,
    }

    # --------------------------------------------------------------------------------
    G = network["G"]
    nodes = network["nodes"]
    edges = network["edges"]

    # --------------------------------------------------------------------------------
    pos = nx.spring_layout(G, k=k, iterations=iterations)

    colors = [node[1]["group"] for node in nodes]
    colors = [group_colors[color] for color in colors]

    sizes = [node[1]["size"] for node in G.nodes.data()]
    sizes = [size / max(sizes) for size in sizes]
    max_size = max(sizes)
    min_size = min(sizes)
    if max_size == min_size:
        sizes = [500 for size in sizes]
    else:
        sizes = [
            (size - min_size) / (max_size - min_size) * 1400 + 100 for size in sizes
        ]

    edge_colors = ["silver"] * len(edges)

    # draws the network
    nx.draw(
        G,
        node_color=colors,
        node_size=sizes,
        edge_color=edge_colors,
        pos=pos,
        **options,
    )

    # edge color of nodes
    ax.collections[0].set_edgecolor("k")

    # plot centers as black dots
    x_points = [value[0] for value in pos.values()]
    y_points = [value[1] for value in pos.values()]

    size = [
        (node[0], node[1]["size"], pos[node[0]][0], pos[node[0]][1])
        for node in G.nodes.data()
    ]
    size = sorted(size, key=itemgetter(1), reverse=True)
    x_points_marked = [value[2] for value in size[:max_labels]]
    y_points_marked = [value[3] for value in size[:max_labels]]

    ax.scatter(
        x_points_marked,
        y_points_marked,
        marker="o",
        s=20,
        c="k",
        alpha=1.0,
        zorder=10,
    )

    #  Center of the plot
    x_mean = sum(x_points) / len(x_points)
    y_mean = sum(y_points) / len(y_points)

    xlim = ax.get_xlim()
    ylim = ax.get_ylim()

    factor = 0.05
    rx = factor * (xlim[1] - xlim[0])
    ry = factor * (ylim[1] - ylim[0])
    radious = np.sqrt(rx ** 2 + ry ** 2)

    for label in size[:max_labels]:

        label = label[0]

        x_point, y_point = pos[label]

        x_c = x_point - x_mean
        y_c = y_point - y_mean
        angle = np.arctan(np.abs(y_c / x_c))
        x_label = x_point + np.copysign(radious * np.cos(angle), x_c)
        y_label = y_point + np.copysign(radious * np.sin(angle), y_c)

        ha = "left" if x_point > x_mean else "right"
        va = "center"

        ax.text(
            x_label,
            y_label,
            s=label,
            fontsize=7,
            bbox=dict(
                facecolor="w",
                alpha=1.0,
                edgecolor="gray",
                boxstyle="round,pad=0.5",
            ),
            horizontalalignment=ha,
            verticalalignment=va,
            alpha=0.9,
            zorder=13,
        )

        ax.plot(
            [x_point, x_label],
            [y_point, y_label],
            lw=1,
            ls="-",
            c="k",
            zorder=13,
        )

    fig.set_tight_layout(True)

    return fig

ext: py

max_stars_repo_name: jdvelasq/techminer-api

max_issues_repo_name: jdvelasq/techminer-api

max_forks_repo_name: jdvelasq/techminer-api

reason: import numpy,import networkx


text:
import srez_demo
import srez_input
import srez_model
import srez_train

import os.path
import random
import numpy as np
import numpy.random

import tensorflow as tf

FLAGS = tf.app.flags.FLAGS

# Configuration (alphabetically)
tf.app.flags.DEFINE_integer('batch_size', 16,
                            "Number of samples per batch.")

tf.app.flags.DEFINE_string('checkpoint_dir', 'checkpoint',
                           "Output folder where checkpoints are dumped.")

tf.app.flags.DEFINE_integer('checkpoint_period', 10000,
                            "Number of batches in between checkpoints")

tf.app.flags.DEFINE_string('summary_dir', 'summary',
                            "Diretory to save TensorBoard summaries")

tf.app.flags.DEFINE_string('dataset', 'dataset',
                           "Path to the dataset directory.")

tf.app.flags.DEFINE_float('epsilon', 1e-8,
                          "Fuzz term to avoid numerical instability")

tf.app.flags.DEFINE_string('run', 'demo',
                            "Which operation to run. [demo|train]")

tf.app.flags.DEFINE_float('gene_l1_factor', .90,
                          "Multiplier for generator L1 loss term")

tf.app.flags.DEFINE_float('learning_beta1', 0.5,
                          "Beta1 parameter used for AdamOptimizer")

tf.app.flags.DEFINE_float('learning_beta2', 0.9,
                          "Beta1 parameter used for AdamOptimizer")

tf.app.flags.DEFINE_float('learning_rate_start', 0.0001,
                          "Starting learning rate used for AdamOptimizer")

tf.app.flags.DEFINE_float('train_noise', 0.03,
                          "level of Gaussian noise added to training images")

tf.app.flags.DEFINE_integer('learning_rate_half_life', 5000,
                            "Number of batches until learning rate is halved")

tf.app.flags.DEFINE_bool('log_device_placement', False,
                         "Log the device where variables are placed.")

tf.app.flags.DEFINE_integer('sample_size', 64,
                            "Image sample size in pixels. Range [64,128]")

tf.app.flags.DEFINE_integer('summary_period', 200,
                            "Number of batches between summary data dumps")

tf.app.flags.DEFINE_integer('random_seed', 0,
                            "Seed used to initialize rng.")

tf.app.flags.DEFINE_integer('test_vectors', 16,
                            """Number of features to use for testing""")

tf.app.flags.DEFINE_string('train_dir', 'train',
                           "Output folder where training logs are dumped.")

tf.app.flags.DEFINE_integer('train_time', 200,
                            "Time in minutes to train the model")

tf.app.flags.DEFINE_string('optimizer', 'adam',
                            "Which optimizer to use: [rmsprop | adam]?")

tf.app.flags.DEFINE_string('loss_func', 'dcgan',
                            "Which loss function to use: [dcgan | wgan | wgangp]?")

tf.app.flags.DEFINE_integer('weigh_clip', 0,
                            "Perform weigh clipping(1) or not(0)?")

tf.app.flags.DEFINE_string('critic_norm', 'batch',
                            "Normalization in the discriminator: [batch|layer|none]?")

tf.app.flags.DEFINE_integer('output_image', 0,
                            "Present images from all methods (0) or only present images from the generator (1)")

tf.app.flags.DEFINE_integer('LAMBDA', 0,
                            "Gradient penalty lambda hyperparameter in improved WGAN")

tf.app.flags.DEFINE_string('gen_architect', 'resnet',
                            "Artchitect of the generator: [resnet | deconv]")

tf.app.flags.DEFINE_integer('gen_norm', 1, "Normalization in the generator or not")

tf.app.flags.DEFINE_string('input', 'noise', "Input image: [scaled | noise]")

tf.app.flags.DEFINE_integer('noise_dimension', 4, "Dimension of noise input")

def prepare_dirs(delete_train_dir=False):
    # Create checkpoint dir (do not delete anything)
    if not tf.gfile.Exists(FLAGS.checkpoint_dir):
        tf.gfile.MakeDirs(FLAGS.checkpoint_dir)

    # Create summary dir
    if not tf.gfile.Exists(FLAGS.summary_dir):
        tf.gfile.MakeDirs(FLAGS.summary_dir)

    # Cleanup train dir
    if delete_train_dir:
        if tf.gfile.Exists(FLAGS.train_dir):
            tf.gfile.DeleteRecursively(FLAGS.train_dir)
        tf.gfile.MakeDirs(FLAGS.train_dir)

    # Return names of training files
    if not tf.gfile.Exists(FLAGS.dataset) or \
       not tf.gfile.IsDirectory(FLAGS.dataset):
        raise FileNotFoundError("Could not find folder `%s'" % (FLAGS.dataset,))

    filenames = tf.gfile.ListDirectory(FLAGS.dataset)
    filenames = sorted(filenames)
    filenames = [os.path.join(FLAGS.dataset, f) for f in filenames]

    return filenames

def setup_tensorflow():
    # Create session
    config = tf.ConfigProto(log_device_placement=FLAGS.log_device_placement)
    sess = tf.Session(config=config)

    # Initialize rng with a deterministic seed
    with sess.graph.as_default():
        tf.set_random_seed(FLAGS.random_seed)

    random.seed(FLAGS.random_seed)
    np.random.seed(FLAGS.random_seed)

    summary_writer = tf.summary.FileWriter(FLAGS.summary_dir, sess.graph)

    return sess, summary_writer

def _demo():
    # Load checkpoint
    if not tf.gfile.IsDirectory(FLAGS.checkpoint_dir):
        raise FileNotFoundError("Could not find folder `%s'" % (FLAGS.checkpoint_dir,))

    # Setup global tensorflow state
    sess, summary_writer = setup_tensorflow()

    # Prepare directories
    filenames = prepare_dirs(delete_train_dir=False)

    # Setup async input queues
    features, labels = srez_input.setup_inputs(sess, filenames)

    # Create and initialize model
    [gene_minput, gene_moutput,
     gene_output, gene_var_list,
     disc_real_output, disc_fake_output, disc_var_list] = \
            srez_model.create_model(sess, features, labels)

    # Restore variables from checkpoint
    saver = tf.train.Saver()
    filename = 'checkpoint_new.txt'
    filename = os.path.join(FLAGS.checkpoint_dir, filename)
    saver.restore(sess, filename)

    # Execute demo
    srez_demo.demo1(sess)

class TrainData(object):
    def __init__(self, dictionary):
        self.__dict__.update(dictionary)

def _train():
    # Setup global tensorflow state
    sess, summary_writer = setup_tensorflow()

    # Prepare directories
    all_filenames = prepare_dirs(delete_train_dir=True)

    # Separate training and test sets
    # train_filenames = all_filenames[:-FLAGS.test_vectors]
    # test_filenames  = all_filenames[-FLAGS.test_vectors:]

    # We chose a pre-determined set of faces for the convenience of comparing results across models
    determined_test = [73883-1, 110251-1, 36510-1, 132301-1, 57264-1, 152931-1, 93861-1,
    124938-1, 79512-1, 106152-1, 127384-1, 134028-1, 67874-1,
    10613-1, 198694-1, 100990-1]
    all_filenames = np.array(all_filenames)
    train_filenames = list(np.delete(all_filenames, determined_test))

#     test_filenames = list(all_filenames[determined_test])

    # Setup async input queues
    train_features, train_labels = srez_input.setup_inputs(sess, train_filenames)

    # test_features,  test_labels  = srez_input.setup_inputs(sess, test_filenames)

    # Test sets are stored in 'testset_label.npy'
    test_labels = np.load('testset_label.npy')
    test_labels = tf.convert_to_tensor(test_labels, dtype = tf.float32)

    if FLAGS.input == 'scaled':
        test_features = tf.image.resize_area(test_labels, [16, 16])
    elif FLAGS.input == 'noise':
        test_features = tf.random_uniform(shape=[16, FLAGS.noise_dimension, FLAGS.noise_dimension, 3],minval= -1., maxval=1.)

    # Add some noise during training (think denoising autoencoders)
    noise_level = FLAGS.train_noise
    noisy_train_features = train_features + \
                           tf.random_normal(train_features.get_shape(), stddev=noise_level)

    [gene_minput, gene_moutput, gene_output, gene_var_list,
     disc_real_output, disc_fake_output, gradients, disc_var_list] = \
            srez_model.create_model(sess, noisy_train_features, train_labels)

    # >>> add summary scalars for test set
    max_samples = 10 # output 10 test images
    gene_output_clipped = tf.maximum(tf.minimum(gene_moutput, 1.0), 0.)

    # Calculate the L1 error between output samples and labels as a objective measure of image quality
    if FLAGS.input != 'noise':
      l1_quality  = tf.reduce_sum(tf.abs(gene_output_clipped - test_labels), [1,2,3])
      l1_quality = tf.reduce_mean(l1_quality[:max_samples])
      mse_quality  = tf.reduce_sum(tf.square(gene_output_clipped - test_labels), [1,2,3])
      mse_quality = tf.reduce_mean(mse_quality[:max_samples])
      tf.summary.scalar('l1_quality', l1_quality, collections=['test_scalars'])
      tf.summary.scalar('mse_quality', mse_quality, collections=['test_scalars'])

    gene_loss = srez_model.create_generator_loss(disc_fake_output, gene_output, train_features)
    disc_real_loss, disc_fake_loss = \
                     srez_model.create_discriminator_loss(disc_real_output, disc_fake_output)

    # Different training objectives
    if FLAGS.loss_func == 'dcgan':
        # for DCGAN
        disc_loss = tf.add(disc_real_loss, disc_fake_loss, name='disc_loss')
    elif FLAGS.loss_func == 'wgan':
        # for WGAN
        disc_loss = tf.subtract(disc_real_loss, disc_fake_loss, name='disc_loss')
    elif FLAGS.loss_func == 'wgangp':
        # for WGANGP
        disc_loss = tf.subtract(disc_real_loss, disc_fake_loss)
        slopes = tf.sqrt(tf.reduce_sum(tf.square(gradients), reduction_indices=[1]))
        gradient_penalty = tf.reduce_mean((slopes-1.)**2)
        disc_loss = tf.add(disc_loss, FLAGS.LAMBDA*gradient_penalty, name='disc_loss')

    (global_step, learning_rate, gene_minimize, disc_minimize, d_clip) = \
            srez_model.create_optimizers(gene_loss, gene_var_list, disc_loss, disc_var_list)

    # For tensorboard
    tf.summary.scalar('generator_loss', gene_loss)
    tf.summary.scalar('discriminator_real_loss', disc_real_loss)
    tf.summary.scalar('discriminator_fake_loss', disc_fake_loss)
    tf.summary.scalar('discriminator_tot_loss', disc_loss)

    # Train model
    train_data = TrainData(locals())
    srez_train.train_model(train_data)

def main(argv=None):
    # Training or showing off?

    if FLAGS.run == 'demo':
        _demo()
    elif FLAGS.run == 'train':
        _train()

if __name__ == '__main__':
  tf.app.run()

max_stars_repo_name: MandyZChen/srez

reason: import numpy


text:
#!/usr/bin/env python3
import logging
import os
import unittest
from collections import OrderedDict
import joblib as par
import numpy as np
from rdkit.Chem import (MolFromMolBlock, GetMolFrags)
from rdkit.Chem.AllChem import Compute2DCoords
from rdkit.Chem.Draw import MolToFile
from scipy.sparse import dok_matrix
from common_wrangler.common import (InvalidDataError, capture_stdout, silent_remove, diff_lines, capture_stderr)
from ligninkmc.create_lignin import (DEF_TEMP, calc_rates, create_initial_monomers, create_initial_events,
                                     degree, create_initial_state, overall_branching_coefficient,
                                     adj_analysis_to_stdout, get_bond_type_v_time_dict)
from ligninkmc.kmc_common import (Event, Monomer, G, S, H, C, C5O4, OX, C5C5, B5, BB, BO4, AO4, B1, DEF_RXN_RATES,
                                  MON_OLI, MONOMER, GROW, TIME, MONO_LIST, ADJ_MATRIX, CHAIN_LEN, BONDS,
                                  RCF_YIELDS, RCF_BONDS, B1_ALT, DEF_E_BARRIER_KCAL_MOL, MAX_NUM_DECIMAL)
from ligninkmc.kmc_functions import (run_kmc, generate_mol, gen_tcl, find_fragments, fragment_size, break_bond_type,
                                     analyze_adj_matrix, count_oligomer_yields, count_bonds)

# logging.basicConfig(level=logging.DEBUG)
logger = logging.getLogger(__name__)
DISABLE_REMOVE = logger.isEnabledFor(logging.DEBUG)

# Constants #
DATA_DIR = os.path.join(os.path.dirname(__file__), 'test_data')
SUB_DATA_DIR = os.path.join(DATA_DIR, 'run_kmc')

# Output files #
TEST_PNG = os.path.join(SUB_DATA_DIR, 'test.png')
C_LIGNIN_MOL_OUT = os.path.join(SUB_DATA_DIR, 'c_lignin_molfile.txt')
GOOD_C_LIGNIN_MOL_OUT = os.path.join(SUB_DATA_DIR, 'c_lignin_molfile_good.txt')

TCL_FNAME = "psfgen.tcl"
TCL_FILE_LOC = os.path.join(SUB_DATA_DIR, TCL_FNAME)
GOOD_TCL_OUT = os.path.join(SUB_DATA_DIR, "good_psfgen.tcl")
GOOD_TCL_C_LIGNIN_OUT = os.path.join(SUB_DATA_DIR, "good_c_lignin.tcl")
GOOD_TCL_SHORT_SIM_OUT = os.path.join(SUB_DATA_DIR, "good_short_sim.tcl")
GOOD_TCL_NO_GROW_OUT = os.path.join(SUB_DATA_DIR, "good_no_grow.tcl")
GOOD_TCL_SHORT = os.path.join(SUB_DATA_DIR, "good_short.tcl")

# Data #
SHORT_TIME = 0.00001
ADJ_ZEROS = dok_matrix([[0, 0, 0, 0, 0],
                        [0, 0, 0, 0, 0],
                        [0, 0, 0, 0, 0],
                        [0, 0, 0, 0, 0],
                        [0, 0, 0, 0, 0]])

MONO_DRAW_3 = [0.48772, 0.15174, 0.7886]
MONO_DRAW_20 = [0.48772, 0.15174, 0.7886, 0.48772, 0.15174, 0.7886, 0.48772, 0.15174, 0.7886, 0.48772, 0.15174, 0.7886,
                0.48772, 0.15174, 0.7886, 0.48772, 0.15174, 0.7886, 0.48772, 0.15174]

ADJ2 = dok_matrix([[0, 4, 0, 0, 0],
                   [8, 0, 0, 0, 0],
                   [0, 0, 0, 8, 0],
                   [0, 0, 5, 0, 0],
                   [0, 0, 0, 0, 0]])

ADJ3 = dok_matrix([[0, 4, 8, 0, 0],
                   [8, 0, 0, 0, 0],
                   [5, 0, 0, 0, 0],
                   [0, 0, 0, 0, 0],
                   [0, 0, 0, 0, 0]])

def get_avg_num_bonds(bond_type, num_opts, result_list, num_repeats):
    analysis = []
    for i in range(num_opts):
        opt_results = result_list[i]
        cur_adjs = [opt_results[j][ADJ_MATRIX] for j in range(num_repeats)]
        analysis.append([analyze_adj_matrix(cur_adjs[j]) for j in range(num_repeats)])

    num_bonds = [[analysis[j][i][BONDS][bond_type]/sum(analysis[j][i][BONDS].values())
                  for i in range(num_repeats)] for j in range(num_opts)]
    av_num_bonds = [np.mean(bond_pcts) for bond_pcts in num_bonds]
    std_num_bonds = [np.sqrt(np.var(bond_pcts)) for bond_pcts in num_bonds]
    return av_num_bonds, std_num_bonds

def create_sample_kmc_result(max_time=1., num_initial_monos=3, max_monos=10, sg_ratio=0.75, seed=10):
    # The set lists are to minimize randomness in testing (adding while debugging source of randomness in some tests;
    #     leaving because it doesn't hurt a thing; also leaving option to make a monomer_draw of arbitrary length
    #     using a seed, but rounding those numbers because the machine precision differences in floats was the bug
    np.random.seed(seed)
    if num_initial_monos == 3:
        monomer_draw = MONO_DRAW_3
    elif num_initial_monos == 20:
        monomer_draw = [0.77132064, 0.02075195, 0.63364823, 0.74880388, 0.49850701, 0.22479665, 0.19806286,
                        0.76053071, 0.16911084, 0.08833981, 0.68535982, 0.95339335, 0.00394827, 0.51219226,
                        0.81262096, 0.61252607, 0.72175532, 0.29187607, 0.91777412, 0.71457578]
    else:
        monomer_draw = np.around(np.random.rand(num_initial_monos), MAX_NUM_DECIMAL)

    # these are tested separately elsewhere
    initial_monomers = create_initial_monomers(sg_ratio, monomer_draw)
    initial_events = create_initial_events(initial_monomers, DEF_RXN_RATES)
    initial_state = OrderedDict(create_initial_state(initial_events, initial_monomers))
    initial_events.append(Event(GROW, [], rate=1e4))
    result = run_kmc(DEF_RXN_RATES, initial_state, initial_events, n_max=max_monos, t_max=max_time,
                     random_seed=10, sg_ratio=sg_ratio)
    return result

def create_sample_kmc_result_c_lignin(num_monos=2, max_monos=12, seed=10):
    initial_monomers = [Monomer(C, i) for i in range(num_monos)]
    # noinspection PyTypeChecker
    initial_events = create_initial_events(initial_monomers, DEF_RXN_RATES)
    initial_state = create_initial_state(initial_events, initial_monomers)
    initial_events.append(Event(GROW, [], rate=1e4))
    result = run_kmc(DEF_RXN_RATES, initial_state, sorted(initial_events), n_max=max_monos, t_max=2, random_seed=seed)
    return result

def get_avg_num_bonds_single_option(bond_type, result_list, num_repeats):
    analysis = []
    for j in range(num_repeats):
        run_results = result_list[j]
        cur_adjs = run_results[ADJ_MATRIX]
        analysis.append(analyze_adj_matrix(cur_adjs))

    num_bonds = [analysis[j][BONDS][bond_type]/sum(analysis[j][BONDS].values())
                 for j in range(num_repeats)]
    av_num_bonds = np.mean(num_bonds)
    std_num_bonds = np.sqrt(np.var(num_bonds))
    return av_num_bonds, std_num_bonds

# Tests #

class TestCalcRates(unittest.TestCase):
    """
    Tests calculation of rate coefficients by the Eyring equation.
    """
    def testCalcRatesFromKcalMol(self):
        rxn_rates = calc_rates(DEF_TEMP, ea_kcal_mol_dict=DEF_E_BARRIER_KCAL_MOL)
        self.assertTrue(len(rxn_rates) == len(DEF_RXN_RATES))
        rxn_type, substrate, substrate_type = None, None, None  # to make IDE happy
        try:
            for rxn_type in DEF_RXN_RATES:
                for substrate in DEF_RXN_RATES[rxn_type]:
                    for substrate_type in DEF_RXN_RATES[rxn_type][substrate]:
                        self.assertAlmostEqual(DEF_RXN_RATES[rxn_type][substrate][substrate_type],
                                               rxn_rates[rxn_type][substrate][substrate_type])
        except (TypeError, IndexError) as e:
            print(f'{e}\nError when looking at rxn_type: {rxn_type} substrate: {substrate}    '
                  f'substrate_type:    {substrate_type}')

class TestMonomers(unittest.TestCase):
    def testCreateConiferyl(self):
        mon = Monomer(G, 0)  # Makes a guaiacol type monomer with ID = 0
        self.assertTrue(mon.open == {8, 4, 5})
        self.assertTrue(str(mon) == '0: coniferyl alcohol is connected to {0} and active at position 0')

    def testCreateSyringol(self):
        mon = Monomer(S, 2)  # Makes a syringol type monomer with ID = 2
        self.assertTrue(mon.open == {4, 8})
        self.assertTrue(mon.connectedTo == {2})
        self.assertTrue(str(mon) == '2: sinapyl alcohol is connected to {2} and active at position 0')
        self.assertTrue(repr(mon) == '2: sinapyl alcohol \n')

    def testHUnit(self):
        # todo: update once H is added
        try:
            mon = Monomer(H, 2)
            # H type is not currently implemented
            self.assertFalse(mon)  # should not be reached
        except InvalidDataError as e:
            self.assertTrue("only the following" in e.args[0])

    def testUnknownUnit(self):
        try:
            mon = Monomer("@", 2)
            # not a real type
            self.assertFalse(mon)  # should not be reached
        except InvalidDataError as e:
            self.assertTrue("only the following" in e.args[0])

    def testHash(self):
        mon1 = Monomer(S, 5)
        mon2 = Monomer(S, 5)
        check_set = {mon1, mon2}
        self.assertTrue(len(check_set) == 1)

class TestEvent(unittest.TestCase):
    def testIDRepr(self):
        rxn = OX
        # noinspection PyTypeChecker
        event1 = Event(rxn, [2], DEF_RXN_RATES[rxn][G][MONOMER])
        self.assertTrue(str(event1) == "Performing oxidation on index 2")

    def testIDReprBond(self):
        rxn = BO4
        # noinspection PyTypeChecker
        event1 = Event(rxn, [1, 2], DEF_RXN_RATES[rxn][(G, S)][MON_OLI], (4, 8))
        good_str = "Forming bo4 bond between indices [1, 2] (adjacency_matrix update (4, 8))"
        self.assertTrue(str(event1) == good_str)
        self.assertTrue(repr(event1) == good_str)

    def testEventIDHash(self):
        monomer_a = Monomer(S, 4)
        monomer_b = Monomer(S, 4)
        events_a = create_initial_events([monomer_a], DEF_RXN_RATES)
        events_b = create_initial_events([monomer_b], DEF_RXN_RATES)
        self.assertTrue(events_a == events_b)
        check_set = {events_a[0], events_b[0]}
        self.assertTrue(len(check_set) == 1)

class TestCreateInitialMonomers(unittest.TestCase):
    def testCreate3Monomers(self):
        initial_monomers = create_initial_monomers(0.75, [0.48772, 0.15174, 0.7886])
        self.assertTrue(len(initial_monomers) == 3)
        self.assertTrue(initial_monomers[0].type == S)
        self.assertTrue(initial_monomers[1].type == S)
        self.assertTrue(initial_monomers[2].type == G)
        self.assertTrue(initial_monomers[1] < initial_monomers[2])
        self.assertFalse(initial_monomers[0] == initial_monomers[1])

class TestState(unittest.TestCase):
    def testCreateInitialState(self):
        sg_ratio = 0.75
        monomer_draw = [0.48772, 0.15174, 0.7886]
        initial_monomers = create_initial_monomers(sg_ratio, monomer_draw)
        initial_events = create_initial_events(initial_monomers, DEF_RXN_RATES)
        initial_state = create_initial_state(initial_events, initial_monomers)
        self.assertTrue(len(initial_state) == 3)
        self.assertTrue(str(initial_monomers[0]) == str(initial_state[0][MONOMER]))

class TestRunKMC(unittest.TestCase):
    def testMissingRequiredSGRatio(self):
        # set up variable to allow running run_kmc without specifying sg_ratio
        initial_sg_ratio = 0.75
        num_initial_monos = 3
        monomer_draw = np.around(np.random.rand(num_initial_monos), MAX_NUM_DECIMAL)
        # these are tested separately
        initial_monomers = create_initial_monomers(initial_sg_ratio, monomer_draw)
        initial_events = create_initial_events(initial_monomers, DEF_RXN_RATES)
        initial_state = create_initial_state(initial_events, initial_monomers)
        events = {initial_events[i] for i in range(num_initial_monos)}
        events.add(Event(GROW, [], rate=1e4))
        try:
            run_kmc(DEF_RXN_RATES, initial_state, sorted(events), n_max=20, t_max=1, random_seed=10)
            self.assertFalse("Should not arrive here; An error should have be raised")
        except InvalidDataError as e:
            self.assertTrue("A numeric sg_ratio" in e.args[0])

    def testSampleRunKMC(self):
        result = create_sample_kmc_result()
        self.assertTrue(len(result[TIME]) == 41)
        self.assertAlmostEqual(result[TIME][-1], 0.002627864716612143)
        self.assertTrue(len(result[MONO_LIST]) == 10)
        self.assertTrue(str(result[MONO_LIST][-1]) == '9: sinapyl alcohol is connected to '
                                                      '{0, 1, 2, 3, 4, 5, 6, 7, 8, 9} and active at position 4')
        good_dok_keys = [(0, 1), (1, 0), (1, 2), (2, 1), (2, 3), (3, 2), (3, 4), (4, 3), (4, 5), (5, 4), (5, 6),
                         (6, 5), (7, 8), (8, 7), (8, 6), (6, 8), (7, 9), (9, 7)]
        good_dok_vals = [8.0, 8.0, 4.0, 8.0, 4.0, 8.0, 5.0, 8.0, 5.0, 8.0, 4.0, 8.0, 8.0, 8.0, 5.0, 4.0, 4.0, 8.0]

        self.assertTrue(list(result[ADJ_MATRIX].keys()) == good_dok_keys)
        self.assertTrue(list(result[ADJ_MATRIX].values()) == good_dok_vals)

    def testSampleRunKMCCLignin(self):
        result = create_sample_kmc_result_c_lignin()
        self.assertTrue(len(result[TIME]) == 35)
        self.assertAlmostEqual(result[TIME][-1], 0.001749379102515669)
        self.assertTrue(len(result[MONO_LIST]) == 12)
        self.assertTrue(str(result[MONO_LIST][-1]) == '11: caffeoyl alcohol is connected to '
                                                      '{9, 10, 11} and active at position 0')
        good_dok_keys = [(1, 0), (0, 1), (0, 2), (2, 0), (3, 2), (2, 3), (5, 4), (4, 5), (6, 7), (7, 6), (8, 7),
                         (7, 8), (10, 9), (9, 10), (11, 9), (9, 11)]
        good_dok_vals = [5.0, 8.0, 4.0, 8.0, 8.0, 4.0, 8.0, 5.0, 5.0, 8.0, 8.0, 4.0, 5.0, 8.0, 8.0, 4.0]
        self.assertTrue(list(result[ADJ_MATRIX].keys()) == good_dok_keys)
        self.assertTrue(list(result[ADJ_MATRIX].values()) == good_dok_vals)

# noinspection DuplicatedCode
class TestAnalyzeKMCParts(unittest.TestCase):
    def testFindOneFragment(self):
        a = dok_matrix((2, 2))
        frags, branches = find_fragments(a)
        good_frags = [{0}, {1}]
        self.assertEqual(frags, good_frags)
        good_branches = [0, 0]
        self.assertEqual(branches, good_branches)

    def testFindTwoFragments(self):
        a_array = [[0., 1., 1., 0., 0.],
                   [1., 0., 0., 0., 0.],
                   [1., 0., 0., 0., 0.],
                   [0., 0., 0., 0., 1.],
                   [0., 0., 0., 1., 0.]]
        a = dok_matrix(a_array)
        frags, branches = find_fragments(a)
        good_frags = [{0, 1, 2}, {3, 4}]
        self.assertEqual(frags, good_frags)
        good_branches = [0, 0]
        self.assertEqual(branches, good_branches)

    def testFindThreeFragments(self):
        # does not increase coverage, but that's okay
        a = dok_matrix((5, 5))
        a[0, 4] = 1
        a[4, 0] = 1
        frags, branches = find_fragments(a)
        good_frags = [{0, 4}, {1}, {2}, {3}]
        self.assertEqual(frags, good_frags)
        good_branches = [0, 0, 0, 0]
        self.assertEqual(branches, good_branches)

    def testFragmentSize1(self):
        frags = [{0}, {1}]
        result = fragment_size(frags)
        good_result = {0: 1, 1: 1}
        self.assertEqual(result, good_result)

    def testFragmentSize2(self):
        # Does not increase coverage; keep anyway
        frags = [{0, 4, 2}, {1, 3}]
        result = fragment_size(frags)
        good_result = {0: 3, 2: 3, 4: 3, 1: 2, 3: 2}
        self.assertEqual(result, good_result)

    def testFragmentSize3(self):
        frags = [{0, 1, 2, 3, 4}]
        result = fragment_size(frags)
        good_result = {0: 5, 1: 5, 2: 5, 3: 5, 4: 5}
        self.assertEqual(result, good_result)

    def testCountYieldsAllMonomers(self):
        good_olig_len_dict = {1: 5}
        good_olig_monos_dict = {1: 5}
        good_olig_branch_dict = {1: 0}
        good_olig_branch_coeff_dict = {1: 0}
        olig_len_dict, olig_monos_dict, olig_branch_dict, olig_branch_coeff_dict = count_oligomer_yields(ADJ_ZEROS)
        self.assertTrue(olig_len_dict == good_olig_len_dict)
        self.assertTrue(olig_monos_dict == good_olig_monos_dict)
        self.assertTrue(olig_branch_dict == good_olig_branch_dict)
        self.assertTrue(olig_branch_coeff_dict == good_olig_branch_coeff_dict)

    def testCountYields1(self):
        adj_1 = dok_matrix([[0, 4, 0, 0, 0],
                            [8, 0, 0, 0, 0],
                            [0, 0, 0, 0, 0],
                            [0, 0, 0, 0, 0],
                            [0, 0, 0, 0, 0]])
        good_olig_len_dict = {1: 3, 2: 1}
        good_olig_monos_dict = {1: 3, 2: 2}
        good_olig_branch_dict = {1: 0, 2: 0}
        good_olig_branch_coeff_dict = {1: 0, 2: 0}
        olig_len_dict, olig_monos_dict, olig_branch_dict, olig_branch_coeff_dict = count_oligomer_yields(adj_1)
        self.assertTrue(olig_len_dict == good_olig_len_dict)
        self.assertTrue(olig_monos_dict == good_olig_monos_dict)
        self.assertTrue(olig_branch_dict == good_olig_branch_dict)
        self.assertTrue(olig_branch_coeff_dict == good_olig_branch_coeff_dict)

    def testCountYields2(self):
        olig_len_dict, olig_monos_dict, olig_branch_dict, olig_branch_coeff_dict = count_oligomer_yields(ADJ2)
        good_olig_len_dict = {1: 1, 2: 2}
        good_olig_monos_dict = {1: 1, 2: 4}
        good_olig_branch_dict = {1: 0, 2: 0}
        good_olig_branch_coeff_dict = {1: 0, 2: 0}
        self.assertTrue(olig_len_dict == good_olig_len_dict)
        self.assertTrue(olig_monos_dict == good_olig_monos_dict)
        self.assertTrue(olig_branch_dict == good_olig_branch_dict)
        self.assertTrue(olig_branch_coeff_dict == good_olig_branch_coeff_dict)

    def testCountYields3(self):
        olig_len_dict, olig_monos_dict, olig_branch_dict, olig_branch_coeff_dict = count_oligomer_yields(ADJ3)
        good_olig_len_dict = {1: 2, 3: 1}
        good_olig_monos_dict = {1: 2, 3: 3}
        good_olig_branch_dict = {1: 0, 3: 0}
        good_olig_branch_coeff_dict = {1: 0, 3: 0}
        self.assertTrue(olig_len_dict == good_olig_len_dict)
        self.assertTrue(olig_monos_dict == good_olig_monos_dict)
        self.assertTrue(olig_branch_dict == good_olig_branch_dict)
        self.assertTrue(olig_branch_coeff_dict == good_olig_branch_coeff_dict)

    def testCountYields4(self):
        adj = dok_matrix((10, 10), dtype=np.float32)
        adj_dict = {(0, 1): 8.0, (1, 0): 8.0, (1, 2): 4.0, (2, 1): 8.0, (2, 3): 4.0, (3, 2): 8.0, (3, 4): 5.0,
                    (4, 3): 8.0, (5, 4): 8.0, (4, 5): 5.0, (5, 6): 4.0, (6, 5): 8.0, (7, 8): 8.0, (8, 7): 8.0,
                    (0, 8): 4.0, (8, 0): 5.0, (8, 9): 4.0, (9, 8): 8.0}
        for key, val in adj_dict.items():
            adj[key] = val
        olig_len_dict, olig_monos_dict, olig_branch_dict, olig_branch_coeff_dict = count_oligomer_yields(adj)
        good_olig_len_dict = {10: 1}
        good_olig_monos_dict = {10: 10}
        good_olig_branch_dict = {10: 1}
        good_olig_branch_coeff_dict = {10: 0.1}
        self.assertTrue(olig_len_dict == good_olig_len_dict)
        self.assertTrue(olig_monos_dict == good_olig_monos_dict)
        self.assertTrue(olig_branch_dict == good_olig_branch_dict)
        self.assertTrue(olig_branch_coeff_dict == good_olig_branch_coeff_dict)

    def testCountBonds(self):
        good_bond_dict = {BO4: 2, B1: 0, BB: 1, B5: 1, C5C5: 0, AO4: 0, C5O4: 0}
        adj_a = dok_matrix([[0, 8, 0, 0, 0],
                            [4, 0, 8, 0, 0],
                            [0, 5, 0, 8, 0],
                            [0, 0, 8, 0, 4],
                            [0, 0, 0, 8, 0]])
        adj_bonds = count_bonds(adj_a)
        self.assertTrue(adj_bonds == good_bond_dict)

    def testBreakBO4Bonds(self):
        good_broken_adj = np.asarray([[0, 0, 0, 0, 0],
                                      [0, 0, 0, 0, 0],
                                      [0, 0, 0, 8, 0],
                                      [0, 0, 8, 0, 0],
                                      [0, 0, 0, 0, 0]])
        a = dok_matrix((5, 5))
        a[1, 0] = 4
        a[0, 1] = 8
        a[2, 3] = 8
        a[3, 2] = 8
        broken_adj = break_bond_type(a, BO4).toarray()
        self.assertTrue(np.array_equal(broken_adj, good_broken_adj))

    def testBreakB1_ALTBonds(self):
        good_broken_adj = np.asarray([[0, 0, 0, 0, 0],
                                      [0, 0, 1, 0, 0],
                                      [0, 8, 0, 0, 0],
                                      [0, 0, 0, 0, 0],
                                      [0, 0, 0, 0, 0]])
        a = dok_matrix([[0, 4, 0, 0, 0],
                        [8, 0, 1, 0, 0],
                        [0, 8, 0, 0, 0],
                        [0, 0, 0, 0, 0],
                        [0, 0, 0, 0, 0]])
        broken_adj = break_bond_type(a, B1_ALT).toarray()
        self.assertTrue(np.array_equal(broken_adj, good_broken_adj))

    def testDegree0(self):
        # Testing all monomers
        adj = ADJ_ZEROS
        mono_degree = list(degree(adj))
        good_mono_degree = []
        self.assertTrue(mono_degree == good_mono_degree)
        branch_coeff = overall_branching_coefficient(adj)
        self.assertTrue(branch_coeff == 0)

    def testBranchCoeff1(self):
        # testing dimers and a monomer
        adj = ADJ2
        mono_degree = list(degree(adj))
        good_mono_degree = [1, 1, 1, 1]
        self.assertTrue(mono_degree == good_mono_degree)
        branch_coeff = overall_branching_coefficient(adj)
        self.assertTrue(branch_coeff == 0)

    def testBranchCoeff2(self):
        # testing timer and dimers
        adj = ADJ3
        branch_degree = list(degree(adj))
        good_branch_degree = [2, 1, 1]
        self.assertTrue(branch_degree == good_branch_degree)
        branch_coeff = overall_branching_coefficient(adj)
        self.assertTrue(branch_coeff == 0)

    def testBranchCoeff3(self):
        adj = dok_matrix((10, 10), dtype=np.float32)
        adj_dict = {(0, 1): 8.0, (1, 0): 8.0, (1, 2): 4.0, (2, 1): 8.0, (2, 3): 4.0, (3, 2): 8.0, (3, 4): 5.0,
                    (4, 3): 8.0, (5, 4): 8.0, (4, 5): 5.0, (5, 6): 4.0, (6, 5): 8.0, (7, 8): 8.0, (8, 7): 8.0,
                    (0, 8): 4.0, (8, 0): 5.0, (8, 9): 4.0, (9, 8): 8.0}
        for key, val in adj_dict.items():
            adj[key] = val
        branch_degree = list(degree(adj))
        good_branch_degree = [2, 2, 2, 2, 2, 2, 1, 1, 3, 1]
        self.assertTrue(branch_degree == good_branch_degree)
        branch_coeff = overall_branching_coefficient(adj)
        good_branch_coeff = 0.1
        self.assertAlmostEqual(branch_coeff, good_branch_coeff)

    def testBranchCoeff4(self):
        # this has 3 fragments: 2 dimers and a pentamer, which is branched
        adj = dok_matrix((9, 9), dtype=np.float32)
        adj_dict = {(4, 8): 8.0, (8, 4): 8.0, (7, 3): 8.0, (3, 7): 8.0, (0, 2): 5.0, (2, 0): 8.0, (4, 2): 4.0,
                    (2, 4): 5.0, (5, 6): 8.0, (6, 5): 5.0, (1, 2): 8.0, (2, 1): 4.0}
        for key, val in adj_dict.items():
            adj[key] = val
        branch_degree = list(degree(adj))
        good_branch_degree = [1, 1, 3, 1, 2, 1, 1, 1, 1]
        self.assertTrue(branch_degree == good_branch_degree)
        branch_coeff = overall_branching_coefficient(adj)
        good_branch_coeff = 1/9
        self.assertAlmostEqual(branch_coeff, good_branch_coeff)
        # Uncomment below to visually check output
        # mol = MolFromMolBlock(block)
        # Compute2DCoords(mol)
        # MolToFile(mol, TEST_PNG, size=(2000, 1000))

class TestAnalyzeKMCSummary(unittest.TestCase):
    def testKMCResultSummary(self):
        result = create_sample_kmc_result()
        summary = analyze_adj_matrix(result[ADJ_MATRIX], break_co_bonds=True)
        self.assertTrue(summary[CHAIN_LEN] == {10: 1})
        self.assertTrue(summary[BONDS] == {BO4: 4, BB: 2, B5: 2, B1: 0, C5O4: 1, AO4: 0, C5C5: 0})
        self.assertTrue(summary[RCF_YIELDS] == {1: 3, 2: 2, 3: 1})
        self.assertTrue(summary[RCF_BONDS] == {BO4: 0, BB: 2, B5: 2, B1: 0, C5O4: 0, AO4: 0, C5C5: 0})

    def testKMCResultSummaryDescription(self):
        break_co_bonds = True
        result = create_sample_kmc_result()
        summary = analyze_adj_matrix(result[ADJ_MATRIX], break_co_bonds=break_co_bonds)
        good_chain_summary = "Lignin KMC created 10 monomers, which formed:\n" \
                             "       1 oligomer(s) of chain length 10, with branching coefficient 0.0"
        good_bond_summary = "composed of the following bond types and number:\n    BO4:    4     BB:    2" \
                            "     B5:    2     B1:    0    5O4:    1    AO4:    0     55:    0"
        good_rcf_chain_summary = "Breaking C-O bonds to simulate RCF results in:\n       3 monomer(s) (chain length " \
                                 "1)\n       2 dimer(s) (chain length 2)\n       1 trimer(s) (chain length 3)"
        good_rcf_bond_summary = "with the following remaining bond types and number:\n    BO4:    0     BB:    2    " \
                                " B5:    2     B1:    0    5O4:    0    AO4:    0     55:    0"
        with capture_stdout(adj_analysis_to_stdout, summary, break_co_bonds=break_co_bonds) as output:
            self.assertTrue(good_chain_summary in output)
            self.assertTrue(good_bond_summary in output)
            self.assertTrue(good_rcf_chain_summary in output)
            self.assertTrue(good_rcf_bond_summary in output)

    def testKMCShortSimResultSummaryDescription(self):
        break_co_bonds = True
        result = create_sample_kmc_result(max_time=SHORT_TIME)
        summary = analyze_adj_matrix(result[ADJ_MATRIX], break_co_bonds=break_co_bonds)
        # adj_analysis_to_stdout(summary)
        good_chain_summary = "Lignin KMC created 3 monomers, which formed:\n" \
                             "       1 trimer(s) (chain length 3)"
        good_bond_summary = "composed of the following bond types and number:\n    BO4:    1     BB:    1" \
                            "     B5:    0     B1:    0    5O4:    0    AO4:    0     55:    0"
        good_rcf_olig_summary = "Breaking C-O bonds to simulate RCF results in:\n       1 monomer(s) (chain " \
                                "length 1)\n       1 dimer(s) (chain length 2)"
        good_rcf_bond_summary = "with the following remaining bond types and number:\n    BO4:    0     BB:    1    " \
                                " B5:    0     B1:    0    5O4:    0    AO4:    0     55:    0"
        with capture_stdout(adj_analysis_to_stdout, summary, break_co_bonds=break_co_bonds) as output:
            self.assertTrue(good_chain_summary in output)
            self.assertTrue(good_bond_summary in output)
            self.assertTrue(good_rcf_olig_summary in output)
            self.assertTrue(good_rcf_bond_summary in output)

    def testKMCShortSimMoreMonosResultSummaryDescription(self):
        break_co_bonds = False
        result = create_sample_kmc_result(max_time=SHORT_TIME, num_initial_monos=20, max_monos=40)
        summary = analyze_adj_matrix(result[ADJ_MATRIX], break_co_bonds=break_co_bonds)
        # adj_analysis_to_stdout(summary)
        good_chain_summary = "Lignin KMC created 21 monomers, which formed:\n       2 monomer(s) (chain length 1)\n" \
                             "       8 dimer(s) (chain length 2)\n       1 trimer(s) (chain length 3)"
        good_bond_summary = "composed of the following bond types and number:\n    BO4:    1     BB:    6" \
                            "     B5:    3     B1:    0    5O4:    0    AO4:    0     55:    0"
        with capture_stdout(adj_analysis_to_stdout, summary, break_co_bonds=break_co_bonds) as output:
            self.assertTrue(good_chain_summary in output)
            self.assertTrue(good_bond_summary in output)

# noinspection DuplicatedCode
class TestVisualization(unittest.TestCase):
    def testMakePNG(self):
        # smoke test only--that it doesn't fail, not that it looks correct (that's outside the scope of this package)
        # The choices shown resulted (at last check) in 3 fragments, one of which has a branch
        try:
            silent_remove(TEST_PNG)
            result = create_sample_kmc_result(num_initial_monos=24, max_monos=24, seed=1, max_time=SHORT_TIME)
            summary = analyze_adj_matrix(result[ADJ_MATRIX])
            adj_analysis_to_stdout(summary, break_co_bonds=False)
            nodes = result[MONO_LIST]
            adj = result[ADJ_MATRIX]
            block = generate_mol(adj, nodes)
            mol = MolFromMolBlock(block)
            Compute2DCoords(mol)
            MolToFile(mol, TEST_PNG, size=(2000, 1200))
            self.assertTrue(os.path.isfile(TEST_PNG))
        finally:
            silent_remove(TEST_PNG, disable=DISABLE_REMOVE)
            pass

    def testMakeTCL(self):
        try:
            silent_remove(TCL_FILE_LOC)
            result = create_sample_kmc_result()
            gen_tcl(result[ADJ_MATRIX], result[MONO_LIST], tcl_fname=TCL_FNAME, chain_id="L",
                    toppar_dir='toppar', out_dir=SUB_DATA_DIR)
            self.assertFalse(diff_lines(TCL_FILE_LOC, GOOD_TCL_OUT))
        finally:
            silent_remove(TCL_FILE_LOC, disable=DISABLE_REMOVE)
            pass

    def testMakeTCLCLignin(self):
        # Only adds 3 lines to coverage... oh well! At least it's quick.
        try:
            seed = 1
            monos = 7
            silent_remove(TCL_FILE_LOC)
            result = create_sample_kmc_result_c_lignin(num_monos=monos, max_monos=monos*2, seed=seed)
            good_last_time = 0.004474718478040326
            self.assertAlmostEqual(result[TIME][-1], good_last_time)
            gen_tcl(result[ADJ_MATRIX], result[MONO_LIST], tcl_fname=TCL_FNAME, chain_id="L", toppar_dir=None,
                    out_dir=SUB_DATA_DIR)
            self.assertFalse(diff_lines(TCL_FILE_LOC, GOOD_TCL_C_LIGNIN_OUT))
        finally:
            silent_remove(TCL_FILE_LOC, disable=DISABLE_REMOVE)
            pass

    def testGenMolCLignin(self):
        # Only testing parts of this package; printing block to ease checking output
        try:
            result = create_sample_kmc_result_c_lignin(num_monos=12, max_monos=14)
            block = generate_mol(result[ADJ_MATRIX], result[MONO_LIST])
            with open(C_LIGNIN_MOL_OUT, "w") as f:
                f.write(block)
            self.assertFalse(diff_lines(C_LIGNIN_MOL_OUT, GOOD_C_LIGNIN_MOL_OUT))
            # # Uncomment below to visually check output
            # mol = MolFromMolBlock(block)
            # Compute2DCoords(mol)
            # MolToFile(mol, TEST_PNG, size=(2000, 1000))
        finally:
            # silent_remove(TEST_PNG, disable=DISABLE_REMOVE)
            silent_remove(C_LIGNIN_MOL_OUT, disable=DISABLE_REMOVE)
            pass

    def testB1BondGenMol(self):
        ini_mono_type_list = [S, S, S, G, S]
        sg_ratio = 1.0
        max_monos = 12
        random_num = 55
        initial_monomers = [Monomer(mono_type, i) for i, mono_type in enumerate(ini_mono_type_list)]
        initial_events = create_initial_events(initial_monomers, DEF_RXN_RATES)
        initial_events.append(Event(GROW, [], rate=1e4))
        initial_state = create_initial_state(initial_events, initial_monomers)
        result = run_kmc(DEF_RXN_RATES, initial_state, initial_events, n_max=max_monos, t_max=2,
                         random_seed=random_num, sg_ratio=sg_ratio)
        nodes = result[MONO_LIST]
        adj = result[ADJ_MATRIX]
        # generate_mol(adj, nodes)
        with capture_stderr(generate_mol, adj, nodes) as output:
            self.assertFalse(output)

        mol = MolFromMolBlock(generate_mol(adj, nodes))
        mols = GetMolFrags(mol)

        analysis = analyze_adj_matrix(adj)
        frag_sizes = analysis[CHAIN_LEN]

        # Make sure there are the same number of separate fragments calculated by RDKIT
        # as we get from just separating the alternate B1
        self.assertEqual(np.sum(list(frag_sizes.values())), len(mols))

    def testDynamics(self):
        # Tests procedures in the Dynamics.ipynb
        # minimize number of random calls during testing (here, set monomer type distribution)
        monomer_type_list = [G, S, G, G, S, S, S, G, S, S, G, G, S, G, G, G, G, S, G, G, G, S, S, G, S, S, G, G, ]
        num_monos = len(monomer_type_list)
        initial_monomers = [Monomer(mono_type, i) for i, mono_type in enumerate(monomer_type_list)]
        initial_events = create_initial_events(initial_monomers, DEF_RXN_RATES)
        initial_state = create_initial_state(initial_events, initial_monomers)
        # since GROW is not added to event_dict, no additional monomers will be added (sg_ratio is thus not needed)
        result = run_kmc(DEF_RXN_RATES, initial_state, sorted(initial_events), random_seed=10, dynamics=True)
        # With dynamics, the MONO_LIST will be a list of monomer lists:
        #    the inner list is the usual MONO_LIST, but here is it saved for every time step
        t_steps = result[TIME]
        expected_num_t_steps = 61
        self.assertEqual(len(t_steps), expected_num_t_steps)
        self.assertTrue(len(result[MONO_LIST]) == expected_num_t_steps)
        self.assertTrue(len(result[MONO_LIST][-1]) == num_monos)
        # want dict[key: [], ...] where the inner list is values by timestep
        #                         instead of list of time steps with [[key: val, ...], ... ]
        adj_list = result[ADJ_MATRIX]
        (bond_type_dict, olig_len_dict, sum_list, olig_count_dict,
            sum_count_list) = get_bond_type_v_time_dict(adj_list, sum_len_larger_than=10)

        # test results by checking sums
        good_bond_type_sum_dict = {BO4: 16, BB: 171, B1: 0, B5: 119, C5C5: 0, AO4: 0, C5O4: 23}
        bond_type_sum_dict = {}
        for bond_type, val_list in bond_type_dict.items():
            self.assertEqual(len(val_list), expected_num_t_steps)
            bond_type_sum_dict[bond_type] = sum(val_list)
        self.assertEqual(bond_type_sum_dict, good_bond_type_sum_dict)

        good_olig_len_sum_dict = {1: 1112, 2: 474, 3: 21, 4: 56, 5: 45}
        olig_len_sum_dict = {}
        for olig_len, val_list in olig_len_dict.items():
            self.assertEqual(len(val_list), expected_num_t_steps)
            olig_len_sum_dict[olig_len] = sum(val_list)
        self.assertEqual(olig_len_sum_dict, good_olig_len_sum_dict)

    def testIniRates(self):
        # Note: this test did not increase coverage. Added to help debug notebook.
        run_multi = False
        num_repeats = 4
        if run_multi:
            fun = par.delayed(run_kmc)
            num_jobs = num_repeats
        else:
            fun = None
            num_jobs = None

        sg_ratio = 1.1

        # minimize random calls
        monomer_type_list = [S, G]
        initial_monomers = [Monomer(mono_type, i) for i, mono_type in enumerate(monomer_type_list)]
        max_monos = 12
        initial_events = create_initial_events(initial_monomers, DEF_RXN_RATES)
        # FYI: np.logspace(start, stop, num=50, endpoint=True, base=10.0, dtype=None, axis=0)[source]
        num_rates = 3
        add_rates = np.logspace(4, 12, num_rates)
        add_rates_result_list = []

        # will add to random seed in the iterations to insure using a different seed for each repeat
        random_seed = 2

        for add_rate in add_rates:
            initial_state = create_initial_state(initial_events, initial_monomers)
            initial_events.append(Event(GROW, [], rate=add_rate))
            if run_multi:
                results = par.Parallel(n_jobs=num_jobs)([fun(DEF_RXN_RATES, initial_state, initial_events,
                                                             n_max=max_monos, t_max=1, sg_ratio=sg_ratio,
                                                             random_seed=(random_seed + i))
                                                         for i in range(num_repeats)])
            else:
                results = [run_kmc(DEF_RXN_RATES, initial_state, initial_events, n_max=max_monos, t_max=1,
                                   sg_ratio=sg_ratio, random_seed=(random_seed + i)) for i in range(num_repeats)]
            add_rates_result_list.append(results)

        av_bo4_bonds, std_bo4_bonds = get_avg_num_bonds(BO4, num_rates, add_rates_result_list, num_repeats)

        good_av_bo4 = [0.3680555555555555, 0.2863636363636364, 0.03125]
        good_std_bo4 = [0.08187379251771941, 0.013636363636363641, 0.05412658773652741]
        self.assertTrue(np.allclose(av_bo4_bonds, good_av_bo4))
        self.assertTrue(np.allclose(std_bo4_bonds, good_std_bo4))

    def testNoGrowth(self):
        # Here, all the monomers are available at the beginning of the simulation
        try:
            # minimize random calls by providing set list of monomer types
            mono_type_list = [S, S, S, S, G, S, S, S, S, S, S, G, S, S, S, S, S, S, S, S, S, S, S, S]
            random_num = 24
            initial_monomers = [Monomer(mono_type, i) for i, mono_type in
                                enumerate(mono_type_list)]
            initial_events = create_initial_events(initial_monomers, DEF_RXN_RATES)
            initial_state = create_initial_state(initial_events, initial_monomers)
            # since GROW is not added to event_dict, no additional monomers will be added
            result = run_kmc(DEF_RXN_RATES, initial_state, sorted(initial_events), t_max=0.0001,
                             random_seed=random_num)
            gen_tcl(result[ADJ_MATRIX], result[MONO_LIST], tcl_fname=TCL_FNAME, chain_id="L",
                    out_dir=SUB_DATA_DIR)
            self.assertFalse(diff_lines(TCL_FILE_LOC, GOOD_TCL_NO_GROW_OUT))
        finally:
            silent_remove(TCL_FILE_LOC, disable=DISABLE_REMOVE)
            pass

    def testTCLTruncateSegname(self):
        # Tests providing a chain_id that is longer than one character
        try:
            # easier to run_kmc to create monomer_list than recreate it here (adj easier) so doing so
            # minimize random calls by providing set list of monomer types
            initial_mono_type_list = [S, S, G, S, S, S, G, S]
            num_monos = len(initial_mono_type_list)
            initial_monomers = [Monomer(mono_type, i) for i, mono_type in enumerate(initial_mono_type_list)]
            initial_events = create_initial_events(initial_monomers, DEF_RXN_RATES)
            initial_state = create_initial_state(initial_events, initial_monomers)
            # since GROW is not added to event_dict, no additional monomers will be added
            result = run_kmc(DEF_RXN_RATES, initial_state, sorted(initial_events), t_max=2, random_seed=8)
            # quick tests to make sure run_kmc gives expected results (not what we want to test here)
            # self.assertAlmostEqual(result[TIME][-1], 0.000766574526703574)
            self.assertTrue(len(result[MONO_LIST]) == num_monos)
            # the function we want to test here is below
            with capture_stderr(gen_tcl, result[ADJ_MATRIX], result[MONO_LIST], chain_id="lignin",
                                out_dir=SUB_DATA_DIR) as output:
                self.assertTrue("should be one character" in output)
            self.assertFalse(diff_lines(TCL_FILE_LOC, GOOD_TCL_SHORT))
        finally:
            silent_remove(TCL_FILE_LOC, disable=DISABLE_REMOVE)
            pass

    def testTCLGenEmptySegname(self):
        # tcl_fname="psfgen.tcl", psf_fname='lignin', chain_id="L", toppar_dir="toppar/"
        # Here, all the monomers are available at the beginning of the simulation
        # Increases coverage of gen_tcl
        try:
            # easier to run_kmc to create monomer_list than recreate it here (adj easier) so doing so
            # minimize random calls by providing set list of monomer types
            initial_mono_type_list = [S, S, G, S, S, S, G, S]
            num_monos = len(initial_mono_type_list)
            initial_monomers = [Monomer(mono_type, i) for i, mono_type in enumerate(initial_mono_type_list)]
            initial_events = create_initial_events(initial_monomers, DEF_RXN_RATES)
            initial_state = create_initial_state(initial_events, initial_monomers)
            # since GROW is not added to event_dict, no additional monomers will be added
            result = run_kmc(DEF_RXN_RATES, initial_state, sorted(initial_events), t_max=2, random_seed=8)
            # quick tests to make sure run_kmc gives expected results (not what we want to test here)
            self.assertAlmostEqual(result[TIME][-1], 0.00015059250794459398)
            self.assertTrue(len(result[MONO_LIST]) == num_monos)
            # the function we want to test here is below
            with capture_stderr(gen_tcl, result[ADJ_MATRIX], result[MONO_LIST], chain_id=" ",
                                out_dir=SUB_DATA_DIR) as output:
                self.assertTrue("should be one character" in output)
            self.assertFalse(diff_lines(TCL_FILE_LOC, GOOD_TCL_SHORT))
        finally:
            silent_remove(TCL_FILE_LOC, disable=DISABLE_REMOVE)
            pass

    def testCheckBO4Fraction(self):
        # similar to a test above; was useful for comparing output from different versions
        monomer_types = [[G, S, G, G, S, S, S, G, S, S, G, G, S, G, G, G, G, S, G, G, G, S, G, S, S, S, G, S, S, G, G],
                         [S, S, G, G, S, G, S, G, G, G, G, S, S, S, S, S, G, S, S, S, G, G, S, G, S, G, S, S, G, S, S],
                         [S, S, S, S, G, S, S, G, G, S, G, S, G, G, G, G, S, S, S, S, S, S, S, G, S, S, G, S, G, S, G]]
        num_repeats = len(monomer_types)
        sg_result_list = []

        # will add to random seed in the iterations to insure using a different seed for each repeat
        random_seed = 32
        for i in range(num_repeats):
            # Initialize the monomers, event_dict, and state
            initial_monomers = [Monomer(mono_type, m) for m, mono_type in enumerate(monomer_types[i])]
            num_monos = len(initial_monomers)
            initial_events = create_initial_events(initial_monomers, DEF_RXN_RATES)
            initial_state = create_initial_state(initial_events, initial_monomers)
            results = run_kmc(DEF_RXN_RATES, initial_state, initial_events,
                              n_max=num_monos, t_max=2, random_seed=random_seed + i)
            sg_result_list.append(results)

        av_bo4_bonds, std_bo4_bonds = get_avg_num_bonds_single_option(BO4, sg_result_list, num_repeats)
        print("Average fraction BO4 bonds: {:.3f}".format(av_bo4_bonds))
        print("Std dev fraction BO4 bonds: {:.3f}".format(std_bo4_bonds))
        self.assertTrue(np.allclose(av_bo4_bonds, 0.21020733652312598))
        self.assertTrue(np.allclose(std_bo4_bonds, 0.04743254939825481))

max_stars_repo_name: runxuanjiang/lignin-kmc

max_issues_repo_name: runxuanjiang/lignin-kmc

max_forks_repo_name: runxuanjiang/lignin-kmc

reason: import numpy,from scipy


text:
#
# Object detector (by sequential file read from directory)
#
import os
import sys
import random
import math
import re
import time
import numpy as np
import tensorflow as tf
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.patches as patches

from glob import glob
import argparse
import skimage
import shutil

# Root directory of the project
ROOT_DIR = os.path.abspath("../../")

# Import Mask RCNN
sys.path.append(ROOT_DIR)  # To find local version of the library
from mrcnn import utils
from mrcnn import visualize
from mrcnn.visualize import display_images
import mrcnn.model as modellib
from mrcnn.model import log

#from samples.cats_dogs import cats_dogs
from samples.miyukiCamera import miyukiCamera

# Directory to save logs and trained model
MODEL_DIR = os.path.join(ROOT_DIR, "logs")

def process():

    class InferenceConfig(miyukiCamera.MiyukiCameraConfig):
        # Set batch size to 1 since we'll be running inference on
        # one image at a time. Batch size = GPU_COUNT * IMAGES_PER_GPU
        GPU_COUNT = 1
        IMAGES_PER_GPU = 1
    config = InferenceConfig()
    #config.display()

    # Device to load the neural network on.
    # Useful if you're training a model on the same
    # machine, in which case use CPU and leave the
    # GPU for training.
    DEVICE = "/gpu:0"  # /cpu:0 or /gpu:0

    # Inspect the model in training or inference modes
    # values: 'inference' or 'training'
    # TODO: code for 'training' test mode not ready yet
    TEST_MODE = "inference"

    # set model
    # Create model in inference mode
    with tf.device(DEVICE):
        model = modellib.MaskRCNN(mode="inference", config=config,
                                model_dir=MODEL_DIR)

    # Or, load the last model you trained
    weights_path = model.find_last()
    # Load weights
    print("Loading weights ", weights_path)
    model.load_weights(weights_path, by_name=True)

    return model, config

def detector(model,config, dataset, DATA_DIR):

    MRCNN_DATA_DIR = "/".join(  DATA_DIR.split('/')[:-1] )
    MRCNN_DATA_DIR = os.path.join( MRCNN_DATA_DIR, "mrcnn_image")
    print(MRCNN_DATA_DIR)

    images = glob( os.path.join(DATA_DIR, "*.jpg")   )

    print("* total length of images : ", len(images) )

    for f in images:

        print("Running on {}".format(f))
        # Read image
        image = skimage.io.imread(f)
        # Detect objects
        results = model.detect([image], verbose=1)
        r = results[0]

        print("- " * 40 )
        print("Scores --> ",  r['scores'])
        print("found Class Names --> ", [dataset.class_info[i]["name"]  for i in r['class_ids']] )

        classes = [dataset.class_info[i]["name"]  for i in r['class_ids']]

        if "prescription" in classes:
            print("found prescription on %s" % f.split("/")[-1])
            image_file = f.split("/")[-1]
            shutil.copy( f,  os.path.join( MRCNN_DATA_DIR, image_file )  )

def main():

    # Parse command line arguments
    parser = argparse.ArgumentParser(
        description='Sequential Reading File Object Detector.')
    parser.add_argument('--dataset', required=True,
                        metavar="/path/to/balloon/dataset",
                        help='Directory of the target dataset to detect')

    args = parser.parse_args()

    assert args.dataset ,\
            "Provide --image directory to apply detector"

    model, config = process()

    dataset = miyukiCamera.MiyukiCameraDataset()
    DATA_DIR = args.dataset

    detector(model, config, dataset, DATA_DIR)

if __name__ == "__main__":
    main()


text:
from .SpikeUnit import SpikeUnit
import numpy as np

def gaussian_kernel(sigma):
    """
The gaussian kernel.

.. math ::
    w(\\tau) = \\frac{1}{\\sqrt{2 \\pi} \\sigma_w}
    \\exp(-\\frac{\\tau^2}{2 \\sigma^2_w})

example:

.. code-block:: python

    >>> kernel('guassian', {'sigma': 0.4})
    """
    return lambda t: 1/(np.sqrt(2*np.pi)*sigma) * np.exp(-t**2/(2*sigma**2))

def causal_kernel(alpha):
    """
The causal kernel.

.. math::
    w(\\tau) = [\\alpha^2 \\tau \\exp(- \\alpha \\tau)]_+

example:

.. code-block:: python

    >>> kernel('causal', {'alpha': 0.4})
    """
    def causal(t):
        v = alpha**2 * t * np.exp(-alpha*t)
        v[v<0] = 0
        return v
    return causal

def rectangular_kernel(delta):
    """
Moving rectangular kernel.

.. math::
    w(\\tau) = \\begin{cases}
    \\frac{1}{\Delta t} &,\\ -\\frac{\Delta t}{2} \\leq t \\leq \\frac{\\Delta t}{2} \\\\
    0 &,\\ otherwise
    \\end{cases}

example:

.. code-block:: python

    >>> kernel('square', {'delta': 0.4})

    """
    return lambda t: np.ones_like(t[(t>=-delta/2)&(t<=delta/2)]) / delta

def kernel(name, **args):
    """
Get the kernel function.

Kernels:
    - gaussian, args: [sigma]
    - causal, args: [alpha]
    - square, args: [delta]

Args:
    - name:   name of the kernel
    - \*\*args: arguments for each kernel.

Return:
    - kernel function in lambda.

Example:

.. code-block:: python

    >>> kernel('guassian', {'sigma':0.4})

    """
    if name == "gaussian":
        return gaussian_kernel(**args)
    elif name == "causal":
        return causal_kernel(**args)
    elif name=="square":
        return rectangular_kernel(**args)
    else:
        raise NameError("unkown kernel: "+str(name))

def _check_and_convert_to_ndarray(subject):
    if isinstance(subject,SpikeUnit):
        target = subject.spike_train
    elif isinstance(subject,list):
        target = np.array(subject)
    elif isinstance(subject, np.ndarray):
        target = subject
    else:
        raise ValueError("type of to is not supported: "+str(type(subject)))
    return target

def generate_linear_filter(to, k):
    """
Convert the discrete spike train array into a linear filter funciton.

Args:
    - to: data in numpy.ndarray or SpikeUnit class
    - k: kernel

Returns:
    - kernel lambda function

    """
    target = _check_and_convert_to_ndarray(to)
    return lambda t: np.sum(k(target - t))

def apply_linear_filter(to, k, x_range=None, nbins=1000, returnX=True):
    """
Map linear filter into a ndarray.

Args:
    - to: data in numpy.ndarray or SpikeUnit class
    - k: kernel
    - x_range: time range tuple (starttime, endtime), default as None
    - nbins: number of steps, default as 1000
    - returnX: return timespan ndarray

Returns:
    - [X_range,] discrete_linear_filter

    """
    target = _check_and_convert_to_ndarray(to)
    linear_filter = generate_linear_filter(to, k)

    if x_range is None:
        x_range = np.linspace(0, target[-1], nbins)
    elif isinstance(x_range, tuple) or isinstance(x_range, list):
        x_range = np.linspace(x_range[0], x_range[1], nbins)
    elif isinstance(x_range, np.ndarray):
        returnX = False
        pass
    else:
        raise ValueError("invalid x_range")

    if returnX:
        return x_range, np.array(list(map(linear_filter, x_range)))  # XXX
    else:
        return np.array(list(map(linear_filter, x_range)))

def apply_linear_filter_withroi(train, k, starts, roi=(0,0), nbins=1000, pbar=None):
    """
Map linear filter into a ndarray with a region of interest.

Args:
    - to: data in numpy.ndarray or SpikeUnit class
    - k: kernel
    - starts: start time points in list or 1-d ndarray
    - roi: region of interest, in time tuple (starttime, endtime)
    - nbins: number of steps, default as 1000
    - pbar: progress bar in tqdm

Returns:
    - _mean_response: 1-d ndarray
    """
    _mean_response = None

    for each_start in starts:
        a = apply_linear_filter(train, k,
                                x_range=(each_start+roi[0], each_start+roi[1]),
                                nbins=nbins, returnX=False)
        if _mean_response is None:
            _mean_response = a
        else:
            _mean_response = np.vstack((_mean_response, a))

        if pbar:
            pbar.update(1)

    return _mean_response

max_stars_repo_name: ZaneMuir/NeuroAnalysis

max_issues_repo_name: ZaneMuir/NeuroAnalysis

max_forks_repo_name: ZaneMuir/NeuroAnalysis