text:
import unittest
import sys
import os
sys.path.insert(0,os.path.abspath('../'))
import circle
import circlecover
import line
import pdb
import random
import math
import logging
import excessarea
import numpy as np
from line import Line
import printcover
from shapely.geometry import Polygon
from shapely.geometry import Point
from shapely.geometry import LineString

# set up logging to file - see previous section for more details
logging.basicConfig(level=logging.DEBUG,
                    format='%(asctime)s %(name)-12s %(levelname)-8s %(message)s',
                    datefmt='%m-%d %H:%M',
                    filename='debug.txt',
                    filemode='w')
logger = logging.getLogger('circlecover')
# create file handler which logs even debug messages
sh = logging.StreamHandler()
fh = logging.FileHandler('debug.txt')
fh.setLevel(logging.DEBUG)
sh.setLevel(logging.DEBUG)
logger.addHandler(fh)
logger.addHandler(sh)

FIXED_RADIUS = "FixedRadiusPointCover"
VAR_RADIUS = "VariableRadiusLineCover"
AREA_COVER = "AreaCover"

def m_to_km(km):
    return [m/1000.0 for m in km]

def total_area(circle_collection):
    total_area = 0
    for c in circle_collection:
        total_area = total_area + c.area()
    return total_area

class CircleCoverTest(unittest.TestCase):

    def setUp(self):
        pass

    def testMinimumCircleSetCoverForLineSetGreedy(self):
        print("testMinimumAreaCircleSetCoverForLineSetGreedy")
        line_endpoints = [[20,80],[50,70],[60,60],[80,50],[60,40],[80,30],[100,30],[90,20]]
        centers = [(10,70),(30,60),(50,55),(50,30),(60,20)]
        testName =  "testMinimumCircleSetCoverForLineSetGreedy"
        circ,covered_segments = circlecover.min_line_cover_greedy(centers,line_endpoints)
        bp = excessarea.generate_bounding_polygon(line_endpoints,centers)
        printcover.printCover(bp,circ,centers,0,covered_segments,testName, VAR_RADIUS)
        self.assertTrue(len(circ) == len(covered_segments))
        # check that for every line segment, both endpoints are covered by
        # at least one circle in our solution.
        for point in line_endpoints:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)

        # Check if the segments in the cover are distinct subsets.
        for i in range(0,len(covered_segments)):
            for k in range(i+1,len(covered_segments) -1 ):
                for l in covered_segments[i]:
                    for m in covered_segments[k]:
                        self.assertFalse(l == m)

        circ,covered_points = circlecover.min_point_cover_greedy_with_fixed_discs(centers,line_endpoints)
        printcover.printCover(bp,circ,centers,0,covered_segments,testName, FIXED_RADIUS)
        # check that for every line segment, both endpoints are covered by
        # at least one circle in our solution.
        for point in line_endpoints:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)

    def testMinimumCircleSetCoverForLineSetGreedy2(self):
        print("testMinimumCircleSetCoverForLineSetGreedy2")
        line_endpoints = [[20,80],[50,70]]
        centers = [(10,70),(40,60)]
        savedCentrs = list(centers)
        circ,segments = circlecover.min_line_cover_greedy(centers,line_endpoints)
        # check that for every line segment, both endpoints are covered by
        # at least one circle in our solution.
        for point in line_endpoints:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)

        testName = "testMinimumCircleSetCoverForLineSetGreedy2"
        bp = excessarea.generate_bounding_polygon(line_endpoints,centers)
        printcover.printCover(bp,circ,savedCentrs,0,segments,testName, VAR_RADIUS)

    def testMinimumCircleSetCoverForLineSetGreedy3(self):
        print("testMinimumCircleSetCoverForLineSetGreedy3")
        line_endpoints = [[20,80],[50,70]]
        centers = [(10,70),(30,60)]
        polygon_to_cover = excessarea.generate_bounding_polygin(centers,line_endpoints)
        circ = circlecover.min_area_cover_greedy(centers,polygon_to_cover)
        # check that for every line segment, both endpoints are covered by
        # at least one circle in our solution.
        testName = "testMinimumCircleSetCoverForLineSetGreedy3"
        bp = excessarea.generate_bounding_polygon(line_endpoints,centers)
        printcover.printCover(bp,circ,centers,0,[],testName, AREA_COVER)

        for point in line_endpoints:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)

    def testMinimumCircleSetCoverForLineSetGreedyRandom(self):
        print("testMinimumCircleSetCoverForLineSetGreedyRandom")
        line_endpoints = []
        line_segments = []
        random.seed(0)
        start = [100,100]
        end = [200,120]

        p1 = start[0]
        nsteps = 200
        for i in range(1,nsteps):
            p2 = random.randint(start[1],end[1])
            line_endpoints.append((p1,p2))
            p1 = start[0] + float(end[0] - start[0])/float(nsteps) * i

        lineString = LineString(line_endpoints)
        centers = []
        x_coords= [ random.randint(100,200) for i in range(1,150)]
        sorted_xcoords = sorted(x_coords)

        for p1 in sorted_xcoords:
            p2 = random.randint(90,100)
            centers.append((p1,p2))

        for i in range( len(line_endpoints) - 1 ):
            line_segments.append(line.Line(line_endpoints[i],line_endpoints[i+1]))
        circ,included = circlecover.min_line_cover_greedy(centers,line_endpoints,min_center_distance = 20)
        testName = "LineCoverGreedyRandom"
        bp = excessarea.generate_bounding_polygon(line_endpoints,centers)
        printcover.printCover(bp,circ,centers,20,included,testName, VAR_RADIUS)

        for point in line_endpoints:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)

        u = Polygon()

        for c in circ:
            u  = u.union(c.get_geometry())

        self.assertTrue(u.buffer(1).contains(lineString))

        # Check if the segments in the cover are distinct subsets.
        for i in range(0,len(included)):
            for k in range(i+1,len(included) -1 ):
                for l in included[i]:
                    for m in included[k]:
                        self.assertFalse(l == m)

        circ,included = circlecover.min_point_cover_greedy_with_fixed_discs(centers,line_endpoints,min_center_distance=20)
        for point in line_endpoints:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)
        bp = excessarea.generate_bounding_polygon(line_endpoints,centers)
        printcover.printCover(bp,circ,centers,20,[],testName, FIXED_RADIUS)

    def testVB(self):
        print("testEscCoverVB")
        esc_loc_x = m_to_km([1771380,1769310,1769790,1768380,1767390,1764690,1762020,1759920,1753110,1741950,1752210,1757010,1761870,1768230,1772820,1777110,1781610,1786920,1793220])
        esc_loc_y = m_to_km([1827030,1817070,1806990,1797090,1787100,1776840,1767270,1756950,1746690,1735050,1727220,1717290,1707360,1697370,1687320,1677450,1667400,1657350,1647360])
        ship_loc_x = m_to_km([1847012,1844913,1845660,1834150,1823280,1811715,1807512,1806671,1810710,1807769,1817910,1822503,1827218,1823623,1828432,1842183,1846928,1852378,1858591])
        ship_loc_y = m_to_km([1843636,1833617,1823583,1811442,1799284,1787072,1777140,1767066,1759078,1749183,1741311,1731358,1721401,1709309,1699318,1691518,1681523,1671542,1661589])

        centers = []
        for i in range(0,len(esc_loc_x)):
            center = (esc_loc_x[i],esc_loc_y[i])
            centers.append(center)

        line_endpoints = []
        for i in range(0,len(ship_loc_x)):
            p = (ship_loc_x[i],ship_loc_y[i])
            line_endpoints.append(p)

        circ,included = circlecover.min_line_cover_greedy(centers,line_endpoints,min_center_distance = 60)
        testName = "VirginiaBeach"
        printcover.printCover(line_endpoints,circ,centers,60,included,testName,VAR_RADIUS)
        for point in line_endpoints:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)

        circ,included = circlecover.min_point_cover_greedy_with_fixed_discs(centers,line_endpoints,min_center_distance = 60)
        printcover.printCover(line_endpoints,circ,centers,60,[],testName,FIXED_RADIUS)
        for point in line_endpoints:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)

        polygon_to_cover = excessarea.generate_bounding_polygin(centers,line_endpoints)
        circ = circlecover.min_area_cover_greedy(centers,polygon_to_cover,min_center_distance = 60)
        printcover.printCover(polygon_to_cover,circ,centers,60,[],testName,AREA_COVER)
        for point in line_endpoints:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)

    def testSF(self):
        esc_loc_x = m_to_km([-2300850,-2297160,-2284680,-2283390,-2284800,-2289540,-2287620,-2287740,-2287620,-2291760,-2289540,-2283720,-2279730,-2254320,-2252430,-2253120,-2256900,-2273160,-2273970,-2273910])
        esc_loc_y = m_to_km([1986840,1977120,1966620,1957680,1947570,1937730,1926720,1917720,1907880,1897830,1887360,1876560,1867620,1852470,1843620,1833720,1824660,1817640,1807710,1797930])
        ship_loc_x = m_to_km([-2414875,-2401190,-2405002,-2406089,-2407670,-2402495,-2402056,-2400759,-2390693,-2394883,-2393517,-2394351,-2393072,-2396507,-2393492,-2394127,-2399780,-2396297,-2387368,-2387021])
        ship_loc_y = m_to_km([2019979,2007266,2001267,1992909,1982826,1970152,1959565,1950068,1937330,1927338,1917048,1908072,1899767,1892623,1883306,1873335,1864772,1852304,1839577,1829662])

        centers = []
        for i in range(0,len(esc_loc_x)):
            center = (esc_loc_x[i],esc_loc_y[i])
            centers.append(center)
        line_endpoints = []
        for i in range(0,len(ship_loc_x)):
            p = (ship_loc_x[i],ship_loc_y[i])
            line_endpoints.append(p)

        line_segments = []
        for i in range( len(line_endpoints) - 1 ):
            line_segments.append(line.Line(line_endpoints[i],line_endpoints[i+1]))
        testName = "SanFrancisco"
        circ,included = circlecover.min_line_cover_greedy(centers,line_endpoints,min_center_distance = 60)
        printcover.printCover(line_endpoints,circ,centers,60,included,testName,VAR_RADIUS)
        for point in line_endpoints:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)
        circ,included = circlecover.min_point_cover_greedy_with_fixed_discs(centers,line_endpoints,min_center_distance = 60)
        printcover.printCover(line_endpoints,circ,centers,60,[],testName,FIXED_RADIUS)
        for point in line_endpoints:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)

        polygon = excessarea.generate_bounding_polygin(centers,line_endpoints)
        circ = circlecover.min_area_cover_greedy(centers,polygon,min_center_distance = 60)
        printcover.printCover(polygon,circ,centers,60,[],testName,AREA_COVER)
        for point in line_endpoints:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)

    def testEstuary(self):
        """
        Test for a deep estuary.
        """
        interference_contour = [(20,55),(35,65),(40,60),(45,65),(50,55)]
        centers = [(20,46),(25,30),(30,20),(40,15),(50,30),(60,50)]

        polygon = excessarea.generate_bounding_polygon(centers,interference_contour)
        circ = circlecover.min_area_cover_greedy(centers,polygon,min_center_distance=0)
        testName = "Estuary"
        printcover.printCover(interference_contour,circ,centers,0,[],testName,AREA_COVER)
        for point in interference_contour:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)

        circ,included = circlecover.min_point_cover_greedy_with_fixed_discs(centers,interference_contour,min_center_distance=0)
        printcover.printCover(interference_contour,circ,centers,0,[],testName,FIXED_RADIUS)
        for point in interference_contour:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)

        circ,included = circlecover.min_line_cover_greedy(centers,interference_contour,min_center_distance=0)
        polygon = excessarea.generate_bounding_polygin(centers,line_endpoints)
        printcover.printCover(polygon,circ,centers,0,[],testName,VAR_RADIUS)
        for point in interference_contour:
            flag = False
            for c in circ:
                if c.inside(point):
                    flag = True
                    break
            self.assertTrue(flag)

max_stars_repo_name: usnistgov/esc-antenna-cover

max_issues_repo_name: usnistgov/esc-antenna-cover

max_forks_repo_name: usnistgov/esc-antenna-cover


text:
#!/usr/bin/env python

# ase.py
# Copyright (c) 2020. All rights reserved.

from typing import Sequence, TypeVar, Union, Dict
import os
import sys
import numpy as np

from scipy.linalg import orthogonal_procrustes

from d3m import utils
from d3m import container
from d3m import exceptions
from d3m.primitive_interfaces.unsupervised_learning import UnsupervisedLearnerPrimitiveBase
from d3m.metadata import hyperparams, base as metadata_module, params
from d3m.primitive_interfaces import base
from d3m.primitive_interfaces.base import CallResult

Inputs = container.pandas.DataFrame
Outputs = container.pandas.DataFrame

class Params(params.Params):
    W: container.ndarray

class Hyperparams(hyperparams.Hyperparams):
    pass

class PartialProcrustes(UnsupervisedLearnerPrimitiveBase[Inputs, Outputs, Params, Hyperparams]):
    """
    Creates a similarity matrix from pairwise distances and nominates one-to-one
    smallest distance vertex match.
    """
    # This should contain only metadata which cannot be automatically determined from the code.
    metadata = metadata_module.PrimitiveMetadata({
        # Simply an UUID generated once and fixed forever. Generated using "uuid.uuid4()".
        'id': 'd830b3c3-e8f1-480d-9970-c404e5064edd',
        'version': "0.1.0",
        'name': "jhu.partial_procrustes",
        # Keywords do not have a controlled vocabulary. Authors can put here whatever they find suitable.
        'keywords': ['procrustes', 'matching', 'alignment'],
        'source': {
            'name': "JHU",
            'uris': [
                # Unstructured URIs. Link to file and link to repo in this case.
                'https://github.com/neurodata/primitives-interfaces/blob/master/jhu_primitives/partial_procrustes/partial_procrustes.py',
                'https://github.com/neurodata/primitives-interfaces',
            ],
            'contact': 'mailto:asaadel1@jhu.edu'
        },
        'description': 'Aligns two datasets based on the rotation for a seeded group of entries',
        'hyperparams_configuration': {
            # 'max_dimension': 'The maximum dimension that can be used for eigendecomposition',
            # 'which_elbow': 'The scree plot "elbow" to use for dimensionality reduction. High values leads to more dimensions selected.',
            # 'use_attributes': 'Boolean which indicates whether to use the attributes of the nodes.'
        },
        # A list of dependencies in order. These can be Python packages, system
        # packages, or Docker images. Of course Python packages can also have
        # their own dependencies, but sometimes it is necessary to install a
        # Python package first to be even able to run setup.py of another
        # package. Or you have a dependency which is not on PyPi.
        'installation': [
            {
            'type': 'UBUNTU',
            'package': 'libxml2-dev',
            'version': '2.9.4'
            },
            {
            'type': 'UBUNTU',
            'package': 'libpcre3-dev',
            'version': '2.9.4'
            },
            {
            'type': 'PIP',
            'package_uri': 'git+https://github.com/neurodata/primitives-interfaces.git@{git_commit}#egg=jhu_primitives'.format(
                git_commit=utils.current_git_commit(os.path.dirname(__file__)),
             ),
        }],
        # URIs at which one can obtain code for the primitive, if available.
        # 'location_uris': [
        #     'https://gitlab.com/datadrivendiscovery/tests-data/raw/{git_commit}/primitives/test_primitives/monomial.py'.format(
        #         git_commit=utils.current_git_commit(os.path.dirname(__file__)),
        #     ),
        # ],
        # The same path the primitive is registered with entry points in setup.py.
        'python_path': 'd3m.primitives.graph_matching.partial_procrustes.JHU',
        # Choose these from a controlled vocabulary in the schema. If anything
        # is missing which would best describe the primitive, make a merge
        # request.
        'algorithm_types': [
            "RANDOM_GRAPH"
        ],
        'primitive_family':
            'GRAPH_MATCHING',
        'preconditions': [
            'NO_MISSING_VALUES'
        ]
    })

    def __init__(self, *, hyperparams: Hyperparams, random_seed: int = 0,
                 docker_containers: Dict[str, base.DockerContainer] = None) -> None:
        super().__init__(hyperparams=hyperparams,
                         random_seed=random_seed,
                         docker_containers=docker_containers)
        self._fitted: bool = False
        self._training_inputs: Inputs = None
        self._w: container.ndarray = None

    def fit(self, *, timeout: float = None, iterations: int = None) -> CallResult[None]:
        if self._fitted:
            return base.CallResult(None)

        xhat = self._inputs_1
        yhat = self._inputs_2

        temp_train = self._reference.merge(xhat, how='left', on='e_nodeID')
        temp_train = temp_train.merge(yhat, how='left', on='g_nodeID')

        temp_train = temp_train[temp_train['match'].astype(int).astype(bool)]

        xhat_train = temp_train.values[:, 4:-300].astype(np.float32)
        yhat_train = temp_train.values[:, -300:].astype(np.float32)

        w, _ = orthogonal_procrustes(yhat_train, xhat_train)
        self._w = container.ndarray(w)
        self._fitted = True

        return CallResult(None)

    def produce(self, *,
                inputs_1: Inputs,
                inputs_2: Inputs,
                reference: Inputs,
                timeout: float = None,
                iterations: int = None) -> CallResult[Outputs]:
        xhat = inputs_1
        yhat = inputs_2

        yhat_embedding_align = yhat.values[:,1:].astype(np.float32) @ self._w
        yhat_align = yhat.copy()
        yhat_align[yhat.columns[1:]] = yhat_embedding_align

        return base.CallResult(yhat_align,
                               has_finished = True, iterations_done = 1)

        # return base.CallResult(reference, #results,
        #                        has_finished=True,
        #                        iterations_done=1)

    def multi_produce(self, *,
                      produce_methods: Sequence[str],
                      inputs_1: Inputs,
                      inputs_2: Inputs,
                      reference: Inputs,
                      timeout: float = None,
                      iterations: int = None) -> base.MultiCallResult:  # type: ignore
        return self._multi_produce(produce_methods=produce_methods,
                                   timeout=timeout,
                                   iterations=iterations,
                                   inputs_1=inputs_1,
                                   inputs_2=inputs_2,
                                   reference=reference)

    def fit_multi_produce(self, *,
                          produce_methods: Sequence[str],
                          inputs_1: Inputs,
                          inputs_2: Inputs,
                          reference: Inputs,
                          timeout: float = None,
                          iterations: int = None) -> base.MultiCallResult:  # type: ignore
        return self._fit_multi_produce(produce_methods=produce_methods,
                                       timeout=timeout,
                                       iterations=iterations,
                                       inputs_1=inputs_1,
                                       inputs_2=inputs_2,
                                       reference=reference)

    def set_training_data(self, *,
                    inputs_1: Inputs,
                    inputs_2: Inputs,
                    reference: Inputs) -> None:
        self._inputs_1 = inputs_1
        self._inputs_2 = inputs_2
        self._reference = reference

    def get_params(self) -> Params:
        if not self._fitted:
            raise ValueError("Fit not performed.")

        return Params(
            W = self._w)

    def set_params(self, *, params: Params) -> None:
        self._fitted = True
        self._w = params['W']

ext: py

reason: import numpy,from scipy


text:
# Copyright (c) 2022 NVIDIA CORPORATION.  All rights reserved.
# NVIDIA CORPORATION and its licensors retain all intellectual property
# and proprietary rights in and to this software, related documentation
# and any modifications thereto.  Any use, reproduction, disclosure or
# distribution of this software and related documentation without an express
# license agreement from NVIDIA CORPORATION is strictly prohibited.

# include parent path
import numpy as np
import math

import warp as wp
from warp.tests.test_base import *

np.random.seed(42)

wp.init()

@wp.kernel
def sample_mesh_query(mesh: wp.uint64,
                query_points: wp.array(dtype=wp.vec3),
                query_faces: wp.array(dtype=int),
                query_signs: wp.array(dtype=float),
                query_dist: wp.array(dtype=float)):

    tid = wp.tid()

    face_index = int(0)
    face_u = float(0.0)
    face_v = float(0.0)
    sign = float(0.0)

    max_dist = 10012.0

    p = query_points[tid]

    wp.mesh_query_point(mesh, p, max_dist, sign, face_index, face_u, face_v)

    cp = wp.mesh_eval_position(mesh, face_index, face_u, face_v)

    query_signs[tid] = sign
    query_faces[tid] = face_index
    query_dist[tid] = wp.length(cp-p)

@wp.func
def triangle_closest_point(a: wp.vec3, b: wp.vec3, c: wp.vec3, p: wp.vec3):
    ab = b - a
    ac = c - a
    ap = p - a

    d1 = wp.dot(ab, ap)
    d2 = wp.dot(ac, ap)

    if (d1 <= 0.0 and d2 <= 0.0):
        return wp.vec2(1.0, 0.0)

    bp = p - b
    d3 = wp.dot(ab, bp)
    d4 = wp.dot(ac, bp)

    if (d3 >= 0.0 and d4 <= d3):
        return wp.vec2(0.0, 1.0)

    vc = d1 * d4 - d3 * d2
    v = d1 / (d1 - d3)
    if (vc <= 0.0 and d1 >= 0.0 and d3 <= 0.0):
        return wp.vec2(1.0 - v, v)

    cp = p - c
    d5 = wp.dot(ab, cp)
    d6 = wp.dot(ac, cp)

    if (d6 >= 0.0 and d5 <= d6):
        return wp.vec2(0.0, 0.0)

    vb = d5 * d2 - d1 * d6
    w = d2 / (d2 - d6)
    if (vb <= 0.0 and d2 >= 0.0 and d6 <= 0.0):
        return wp.vec2(1.0 - w, 0.0)

    va = d3 * d6 - d5 * d4
    w = (d4 - d3) / ((d4 - d3) + (d5 - d6))
    if (va <= 0.0 and (d4 - d3) >= 0.0 and (d5 - d6) >= 0.0):
        return wp.vec2(0.0, 1.0 - w)

    denom = 1.0 / (va + vb + vc)
    v = vb * denom
    w = vc * denom
    u = 1.0 - v - w

    return wp.vec2(u, v)

@wp.func
def solid_angle(v0: wp.vec3, v1: wp.vec3, v2: wp.vec3, p: wp.vec3):
    a = v0 - p
    b = v1 - p
    c = v2 - p

    a_len = wp.length(a)
    b_len = wp.length(b)
    c_len = wp.length(c)

    det = wp.dot(a, wp.cross(b, c))
    den = a_len*b_len*c_len + wp.dot(a, b)*c_len + wp.dot(b, c)*a_len + wp.dot(c, a)*b_len

    return wp.atan2(det, den)

@wp.kernel
def sample_mesh_brute(
                tri_points: wp.array(dtype=wp.vec3),
                tri_indices: wp.array(dtype=int),
                tri_count: int,
                query_points: wp.array(dtype=wp.vec3),
                query_faces: wp.array(dtype=int),
                query_signs: wp.array(dtype=float),
                query_dist: wp.array(dtype=float)):

    tid = wp.tid()

    min_face = int(0)
    min_dist = float(1.e+6)

    winding_angle = float(0.0)

    p = query_points[tid]

    for i in range(0, tri_count):

        a = tri_points[tri_indices[i*3+0]]
        b = tri_points[tri_indices[i*3+1]]
        c = tri_points[tri_indices[i*3+2]]

        winding_angle += solid_angle(a, b, c, p)

        bary = triangle_closest_point(a, b, c, p)
        u = bary[0]
        v = bary[1]

        cp = u*a + v*b + (1.0-u-v)*c
        cp_dist = wp.length(cp - p)

        if (cp_dist < min_dist):
            min_dist = cp_dist
            min_face = i

    query_faces[tid] = min_face
    query_signs[tid] = winding_angle
    query_dist[tid] = min_dist

# constructs a grid of evenly spaced particles
def particle_grid(dim_x, dim_y, dim_z, lower, radius, jitter):
    points = np.meshgrid(np.linspace(0, dim_x, dim_x), np.linspace(0, dim_y, dim_y), np.linspace(0, dim_z, dim_z))
    points_t = np.array((points[0], points[1], points[2])).T*radius*2.0 + np.array(lower)
    points_t = points_t + np.random.rand(*points_t.shape)*radius*jitter

    return points_t.reshape((-1, 3))

# triangulate a list of polygon face indices
def triangulate(face_counts, face_indices):
    num_tris = np.sum(np.subtract(face_counts, 2))
    num_tri_vtx = num_tris * 3
    tri_indices = np.zeros(num_tri_vtx, dtype=int)
    ctr = 0
    wedgeIdx = 0

    for nb in face_counts:
        for i in range(nb-2):
            tri_indices[ctr] = face_indices[wedgeIdx]
            tri_indices[ctr + 1] = face_indices[wedgeIdx + i + 1]
            tri_indices[ctr + 2] = face_indices[wedgeIdx + i + 2]
            ctr+=3
        wedgeIdx+=nb

    return tri_indices

def test_mesh_query_point(test, device):

    from pxr import Usd, UsdGeom, Gf, Sdf

    mesh = Usd.Stage.Open(os.path.abspath(os.path.join(os.path.dirname(__file__), "assets/torus.usda")))
    mesh_geom = UsdGeom.Mesh(mesh.GetPrimAtPath("/World/Torus"))

    mesh_counts = mesh_geom.GetFaceVertexCountsAttr().Get()
    mesh_indices = mesh_geom.GetFaceVertexIndicesAttr().Get()

    tri_indices = triangulate(mesh_counts, mesh_indices)

    mesh_points = wp.array(np.array(mesh_geom.GetPointsAttr().Get()), dtype=wp.vec3, device=device)
    mesh_indices = wp.array(np.array(tri_indices), dtype=int, device=device)

    # create mesh
    mesh = wp.Mesh(
        points=mesh_points,
        velocities=None,
        indices=mesh_indices)

    p = particle_grid(32, 32, 32, np.array([-5.0, -5.0, -5.0]), 0.1, 0.1)*100.0

    query_count = len(p)
    query_points = wp.array(p, dtype=wp.vec3, device=device)

    signs_query = wp.zeros(query_count, dtype=float, device=device)
    faces_query = wp.zeros(query_count, dtype=int, device=device)
    dist_query = wp.zeros(query_count, dtype=float, device=device)

    signs_brute = wp.zeros(query_count, dtype=float, device=device)
    faces_brute = wp.zeros(query_count, dtype=int, device=device)
    dist_brute = wp.zeros(query_count, dtype=float, device=device)

    wp.launch(kernel=sample_mesh_query, dim=query_count, inputs=[mesh.id, query_points, faces_query, signs_query, dist_query], device=device)
    wp.launch(kernel=sample_mesh_brute, dim=query_count, inputs=[mesh_points, mesh_indices, int(len(mesh_indices)/3), query_points, faces_brute, signs_brute, dist_brute], device=device)

    signs_query = signs_query.numpy()
    faces_query = faces_query.numpy()
    dist_query = dist_query.numpy()

    signs_brute = signs_brute.numpy()
    faces_brute = faces_brute.numpy()
    dist_brute = dist_brute.numpy()

    query_points = query_points.numpy()

    inside_query = []
    inside_brute = []

    for i in range(query_count):

        if (signs_query[i] < 0.0):
            inside_query.append(query_points[i].tolist())

        if (signs_brute[i] > 6.0):
            inside_brute.append(query_points[i].tolist())

    inside_query = np.array(inside_query)
    inside_brute = np.array(inside_brute)

    dist_error = np.max(np.abs(dist_query - dist_brute))
    sign_error = np.max(np.abs(inside_query - inside_brute))

    tolerance = 1.5e-4
    test.assertTrue(dist_error < tolerance, f"dist_error is {dist_error} which is >= {tolerance}")
    test.assertTrue(sign_error < tolerance, f"sign_error is {sign_error} which is >= {tolerance}")

    # import warp.render

    # stage = warp.render.UsdRenderer("tests/outputs/test_mesh_query_point.usd")

    # radius = 10.0
    # stage.begin_frame(0.0)
    # stage.render_mesh(points=mesh_points.numpy(), indices=mesh_indices.numpy(), name="mesh")
    # stage.render_points(points=inside_query, radius=radius, name="query")
    # stage.render_points(points=inside_brute, radius=radius, name="brute")
    # stage.render_points(points=query_points, radius=radius, name="all")
    # stage.end_frame()

    # stage.save()

@wp.kernel
def mesh_query_point_loss(mesh: wp.uint64,
                query_points: wp.array(dtype=wp.vec3),
                projected_points: wp.array(dtype=wp.vec3),
                loss: wp.array(dtype=float)):

    tid = wp.tid()

    face_index = int(0)
    face_u = float(0.0)
    face_v = float(0.0)
    sign = float(0.0)

    max_dist = 10012.0

    p = query_points[tid]

    wp.mesh_query_point(mesh, p, max_dist, sign, face_index, face_u, face_v)
    q = wp.mesh_eval_position(mesh, face_index, face_u, face_v)

    projected_points[tid] = q

    dist = wp.length(wp.sub(p, q))
    loss[tid] = dist

def test_adj_mesh_query_point(test, device):

    from pxr import Usd, UsdGeom, Gf, Sdf

    mesh = Usd.Stage.Open(os.path.abspath(os.path.join(os.path.dirname(__file__), "assets/torus.usda")))
    mesh_geom = UsdGeom.Mesh(mesh.GetPrimAtPath("/World/Torus"))

    mesh_counts = mesh_geom.GetFaceVertexCountsAttr().Get()
    mesh_indices = mesh_geom.GetFaceVertexIndicesAttr().Get()

    tri_indices = triangulate(mesh_counts, mesh_indices)

    mesh_points = wp.array(np.array(mesh_geom.GetPointsAttr().Get()), dtype=wp.vec3, device=device)
    mesh_indices = wp.array(np.array(tri_indices), dtype=int, device=device)

    # test tri
    # print("Testing Single Triangle")
    # mesh_points = wp.array(np.array([[0.0, 0.0, 0.0], [2.0, 0.0, 0.0], [0.0, 2.0, 0.0]]), dtype=wp.vec3, device=device)
    # mesh_indices = wp.array(np.array([0,1,2]), dtype=int, device=device)

    # create mesh
    mesh = wp.Mesh(
        points=mesh_points,
        velocities=None,
        indices=mesh_indices)

    # p = particle_grid(32, 32, 32, np.array([-5.0, -5.0, -5.0]), 0.1, 0.1)*100.0
    p = wp.vec3(50.0, 50.0, 50.0)

    tape = wp.Tape()

    # analytic gradients
    with tape:

        query_points = wp.array(p, dtype=wp.vec3, device=device, requires_grad=True)
        projected_points = wp.zeros(n=1, dtype=wp.vec3, device=device)
        loss = wp.zeros(n=1, dtype=float, device=device, requires_grad=True)

        wp.launch(kernel=mesh_query_point_loss, dim=1, inputs=[mesh.id, query_points, projected_points, loss], device=device)

    tape.backward(loss=loss)
    analytic = tape.gradients[query_points].numpy().flatten()

    # numeric gradients
    eps = 1.e-3
    loss_values = []
    numeric = np.zeros(3)

    offset_query_points = [
        wp.vec3(p[0] - eps, p[1], p[2]), wp.vec3(p[0] + eps, p[1], p[2]),
        wp.vec3(p[0], p[1] - eps, p[2]), wp.vec3(p[0], p[1] + eps, p[2]),
        wp.vec3(p[0], p[1], p[2] - eps), wp.vec3(p[0], p[1], p[2] + eps)]

    for i in range(6):
        q = offset_query_points[i]

        query_points = wp.array(q, dtype=wp.vec3, device=device)
        projected_points = wp.zeros(n=1, dtype=wp.vec3, device=device)
        loss = wp.zeros(n=1, dtype=float, device=device)

        wp.launch(kernel=mesh_query_point_loss, dim=1, inputs=[mesh.id, query_points, projected_points, loss], device=device)

        loss_values.append(loss.numpy()[0])

    for i in range(3):
        l_0 = loss_values[i*2]
        l_1 = loss_values[i*2+1]
        gradient = (l_1 - l_0) / (2.0*eps)
        numeric[i] = gradient

    error = ((analytic - numeric) * (analytic - numeric)).sum(axis=0)

    tolerance = 1.e-3
    test.assertTrue(error < tolerance, f"error is {error} which is >= {tolerance}")

def register(parent):

    devices = wp.get_devices()

    class TestMeshQuery(parent):
        pass

    add_function_test(TestMeshQuery, "test_mesh_query_point", test_mesh_query_point, devices=devices)
    add_function_test(TestMeshQuery, "test_adj_mesh_query_point", test_adj_mesh_query_point, devices=devices)

    return TestMeshQuery

if __name__ == '__main__':
    c = register(unittest.TestCase)
    unittest.main(verbosity=2)

max_stars_repo_name: NVIDIA/warp

max_issues_repo_name: NVIDIA/warp

max_forks_repo_name: NVIDIA/warp

reason: import numpy


text:
import cv2
import numpy as np
import math
from skimage import transform as trans

def get_norm(image, plandmarks, cropSize=(128, 128), ec_mc_y=48, ec_y=40):
    '''
    image: input image
    plandmarks: ((x0, y0), (x1, y1), (x2, y2), (x3, y3), (x4, y4))
    cropSize:  crop face size, default:(128, 128)
    ec_mc_y: distance between two eyes' center and mouth corner center, default:48
    ec_y: distance between two eyes, default:40
    '''

    midEye = ((plandmarks[0][0] + plandmarks[1][0])/2.0, (plandmarks[0][1] + plandmarks[1][1])/2.0)
    midMouth = ((plandmarks[3][0] + plandmarks[4][0])/2.0, (plandmarks[3][1] + plandmarks[4][1])/2.0)
    distance_ec_mc = ((midEye[0] - midMouth[0])**2 + (midEye[1] - midMouth[1])**2)**0.5
    scale = ec_mc_y / distance_ec_mc
    angle = math.atan2((plandmarks[1][1] - plandmarks[0][1]), (plandmarks[1][0] - plandmarks[0][0]))

    M = np.zeros((2, 3))
    M[0][0] = scale*math.cos(angle)
    M[0][1] = scale*math.sin(angle)
    M[1][0] = -scale*math.sin(angle)
    M[1][1] = scale*math.cos(angle)
    M[0][2] = -(midEye[0]*M[0][0] + midEye[1]*M[0][1]) + cropSize[0] / 2
    M[1][2] = -(midEye[0]*M[1][0] + midEye[1]*M[1][1]) + ec_y
    norm = cv2.warpAffine(image, M, cropSize)
    return norm

def get_norm2(img, bbox=None, landmark=None, **kwargs):
    if isinstance(img, str):
        img = read_image(img, **kwargs)
    M = None
    image_size = []
    str_image_size = kwargs.get('image_size', '')
    if len(str_image_size)>0:
        image_size = [int(x) for x in str_image_size.split(',')]
        if len(image_size)==1:
            image_size = [image_size[0], image_size[0]]
        # assert len(image_size)==2
        # assert image_size[0]==112
        # assert image_size[0]==112 or image_size[1]==96
    if landmark is not None:
        assert len(image_size)==2
        src = np.array([
        [30.2946, 51.6963],
        [65.5318, 51.5014],
        [48.0252, 71.7366],
        [33.5493, 92.3655],
        [62.7299, 92.2041] ], dtype=np.float32 )
        if image_size[1]==112:
            src[:,0] += 8.0
        dst = landmark.astype(np.float32)

        tform = trans.SimilarityTransform()
        tform.estimate(dst, src)
        M = tform.params[0:2,:]
        #M = cv2.estimateRigidTransform( dst.reshape(1,5,2), src.reshape(1,5,2), False)

    if M is None:
        if bbox is None: #use center crop
            det = np.zeros(4, dtype=np.int32)
            det[0] = int(img.shape[1]*0.0625)
            det[1] = int(img.shape[0]*0.0625)
            det[2] = img.shape[1] - det[0]
            det[3] = img.shape[0] - det[1]
        else:
            det = bbox
        margin = kwargs.get('margin', 44)
        bb = np.zeros(4, dtype=np.int32)
        bb[0] = np.maximum(det[0]-margin/2, 0)
        bb[1] = np.maximum(det[1]-margin/2, 0)
        bb[2] = np.minimum(det[2]+margin/2, img.shape[1])
        bb[3] = np.minimum(det[3]+margin/2, img.shape[0])
        ret = img[bb[1]:bb[3],bb[0]:bb[2],:]
        if len(image_size)>0:
            ret = cv2.resize(ret, (image_size[1], image_size[0]))
        return ret
    else: #do align using landmark
        assert len(image_size)==2

        #src = src[0:3,:]
        #dst = dst[0:3,:]

        #print(src.shape, dst.shape)
        #print(src)
        #print(dst)
        #print(M)
        warped = cv2.warpAffine(img,M,(image_size[1],image_size[0]), borderValue = 0.0)

        #tform3 = trans.ProjectiveTransform()
        #tform3.estimate(src, dst)
        #warped = trans.warp(img, tform3, output_shape=_shape)
        return warped


text:
import gensim

import numpy as np

from sklearn.base import BaseEstimator, ClassifierMixin

class MeanW2VEmbeddingVectorizer(BaseEstimator, ClassifierMixin):
    def __init__(self):
        pass

    def fit(self, X, y=None):
        X = [val.split() for val in X.to_list()]
        self.model = gensim.models.Word2Vec(X, size=100)
        self.word2vec = dict(zip(self.model.wv.index2word, self.model.wv.vectors))

        if self.word2vec:
            self.dim = next(iter(self.word2vec.values())).shape[0]
        else:
            self.dim = 0
        return self

    def get_params(self, **params):
        return {}

    def transform(self, X):
        return np.array([
            np.mean([self.word2vec[w] for w in words if w in self.word2vec]
                    or [np.zeros(self.dim)], axis=0)
            for words in X
        ])

max_stars_repo_name: castargo/SarcDetectionRusModels

max_issues_repo_name: castargo/SarcDetectionRusModels

max_forks_repo_name: castargo/SarcDetectionRusModels


text:
import numpy as np
from torchstain.normalizers.he_normalizer import HENormalizer
import scipy.linalg

"""
Source code adapted from: https://github.com/schaugf/HEnorm_python
Original implementation: https://github.com/mitkovetta/staining-normalization
"""
class NumpyMacenkoNormalizer(HENormalizer):
    def __init__(self):
        super().__init__()

        self.HERef = np.array([[0.5626, 0.2159],
                          [0.7201, 0.8012],
                          [0.4062, 0.5581]])
        self.maxCRef = np.array([1.9705, 1.0308])

    def __convert_rgb2od(self, I, Io=240, beta=0.15):
        # calculate optical density
        OD = -np.log((I.astype(np.float)+1)/Io)

        # remove transparent pixels
        ODhat = OD[~np.any(OD < beta, axis=1)]

        return OD, ODhat

    def __find_HE(self, ODhat, eigvecs, alpha):
        #project on the plane spanned by the eigenvectors corresponding to the two
        # largest eigenvalues
        That = ODhat.dot(eigvecs[:,1:3])

        phi = np.arctan2(That[:,1],That[:,0])

        minPhi = np.percentile(phi, alpha)
        maxPhi = np.percentile(phi, 100-alpha)

        vMin = eigvecs[:,1:3].dot(np.array([(np.cos(minPhi), np.sin(minPhi))]).T)
        vMax = eigvecs[:,1:3].dot(np.array([(np.cos(maxPhi), np.sin(maxPhi))]).T)

        # a heuristic to make the vector corresponding to hematoxylin first and the
        # one corresponding to eosin second
        if vMin[0] > vMax[0]:
            HE = np.array((vMin[:,0], vMax[:,0])).T
        else:
            HE = np.array((vMax[:,0], vMin[:,0])).T

        return HE

    def __find_concentration(self, OD, HE):
        # rows correspond to channels (RGB), columns to OD values
        Y = np.reshape(OD, (-1, 3)).T

        # determine concentrations of the individual stains
        C = np.linalg.lstsq(HE, Y, rcond=None)[0]

        return C

    def compute_matrices(self, I, Io, alpha, beta):
        I = I.reshape((-1,3))

        OD, ODhat = self.__convert_rgb2od(I, Io=Io, beta=beta)

        # compute eigenvectors
        _, eigvecs = np.linalg.eigh(np.cov(ODhat.T))

        HE = self.__find_HE(ODhat, eigvecs, alpha)

        C = self.__find_concentration(OD, HE)

        # normalize stain concentrations
        maxC = np.array([np.percentile(C[0,:], 99), np.percentile(C[1,:],99)])

        return HE, C, maxC

    # my mod
    def compute_OD_HE(self, I, Io, alpha, beta):
        I = I.reshape((-1,3))
        OD, ODhat = self.__convert_rgb2od(I, Io=Io, beta=beta)
        # compute eigenvectors
        _, eigvecs = scipy.linalg.eigh(np.cov(ODhat.T))
        HE = self.__find_HE(ODhat, eigvecs, alpha)

        return (OD, HE)

    def fit(self, I, Io=240, alpha=1, beta=0.15):
        HE, _, maxC = self.compute_matrices(I, Io, alpha, beta)

        self.HERef = HE
        self.maxCRef = maxC

    def normalize(self, I, precalc_HE=None, precalc_maxC=None, Io=240, alpha=1, beta=0.15, stains=True):
        ''' Normalize staining appearence of H&E stained images

        Example use:
            see test.py

        Input:
            I: RGB input image
            Io: (optional) transmitted light intensity

        Output:
            Inorm: normalized image
            H: hematoxylin image
            E: eosin image

        Reference:
            A method for normalizing histology slides for quantitative analysis. M.
            Macenko et al., ISBI 2009
        '''
        h, w, c = I.shape
        I = I.reshape((-1,3))

        # HE, C, maxC = self.compute_matrices(I, Io, alpha, beta)
        OD, HE = self.compute_OD_HE(I, Io, alpha, beta)

        if precalc_HE is not None:
            HE = precalc_HE
        # if bool_mask is not None:
        #     bool_mask = bool_mask.reshape((-1))
        #     valid_pixels = I[bool_mask]
        #     _, HE = self.compute_OD_HE(valid_pixels, Io, alpha, beta)
        #     # HE, _, maxC = self.compute_matrices(valid_pixels, Io, alpha, beta)

        C = self.__find_concentration(OD, HE)

        if precalc_maxC is not None:
            maxC = precalc_maxC
        else:
            maxC = np.array([np.percentile(C[0,:], 99), np.percentile(C[1,:],99)])

        maxC = np.divide(maxC, self.maxCRef)
        C2 = np.divide(C, maxC[:, np.newaxis])

        # recreate the image using reference mixing matrix
        Inorm = np.multiply(Io, np.exp(-self.HERef.dot(C2)))
        Inorm[Inorm > 255] = 255
        Inorm = np.reshape(Inorm.T, (h, w, c)).astype(np.uint8)

        H, E = None, None

        if stains:
            # unmix hematoxylin and eosin
            H = np.multiply(Io, np.exp(np.expand_dims(-self.HERef[:,0], axis=1).dot(np.expand_dims(C2[0,:], axis=0))))
            H[H > 255] = 255
            H = np.reshape(H.T, (h, w, c)).astype(np.uint8)

            E = np.multiply(Io, np.exp(np.expand_dims(-self.HERef[:,1], axis=1).dot(np.expand_dims(C2[1,:], axis=0))))
            E[E > 255] = 255
            E = np.reshape(E.T, (h, w, c)).astype(np.uint8)

        return Inorm, H, E

max_stars_repo_name: juanigp/torchstain

max_issues_repo_name: juanigp/torchstain

max_forks_repo_name: juanigp/torchstain

reason: import numpy,import scipy


text:
# Greedy layerwise training of a 2 layer autoencoder (MLP) on Fashion MNIST

# Code is based on
# https://github.com/ageron/handson-ml2/blob/master/17_autoencoders_and_gans.ipynb

import superimport

import numpy as np
import matplotlib.pyplot as plt

import os
figdir = "../figures"
def save_fig(fname): plt.savefig(os.path.join(figdir, fname))

import tensorflow as tf
from tensorflow import keras

(X_train_full, y_train_full), (X_test, y_test) = keras.datasets.fashion_mnist.load_data()
X_train_full = X_train_full.astype(np.float32) / 255
X_test = X_test.astype(np.float32) / 255
X_train, X_valid = X_train_full[:-5000], X_train_full[-5000:]
y_train, y_valid = y_train_full[:-5000], y_train_full[-5000:]

def rounded_accuracy(y_true, y_pred):
    return keras.metrics.binary_accuracy(tf.round(y_true), tf.round(y_pred))

def plot_image(image):
    plt.imshow(image, cmap="binary")
    plt.axis("off")

def show_reconstructions(model, images=X_valid, n_images=5):
    reconstructions = model.predict(images[:n_images])
    plt.figure(figsize=(n_images * 1.5, 3))
    for image_index in range(n_images):
        plt.subplot(2, n_images, 1 + image_index)
        plot_image(images[image_index])
        plt.subplot(2, n_images, 1 + n_images + image_index)
        plot_image(reconstructions[image_index])

def train_autoencoder(n_neurons, X_train, X_valid, loss, optimizer,
                      n_epochs=10, output_activation=None, metrics=None):
    n_inputs = X_train.shape[-1]
    encoder = keras.models.Sequential([
        keras.layers.Dense(n_neurons, activation="selu", input_shape=[n_inputs])
    ])
    decoder = keras.models.Sequential([
        keras.layers.Dense(n_inputs, activation=output_activation),
    ])
    autoencoder = keras.models.Sequential([encoder, decoder])
    autoencoder.compile(optimizer, loss, metrics=metrics)
    autoencoder.fit(X_train, X_train, epochs=n_epochs,
                    validation_data=[X_valid, X_valid])
    return encoder, decoder, encoder(X_train), encoder(X_valid)

tf.random.set_seed(42)
np.random.seed(42)

K = keras.backend
X_train_flat = K.batch_flatten(X_train) # equivalent to .reshape(-1, 28 * 28)
X_valid_flat = K.batch_flatten(X_valid)

# Reconstruct binary image
enc1, dec1, X_train_enc1, X_valid_enc1 = train_autoencoder(
    100, X_train_flat, X_valid_flat, "binary_crossentropy",
    keras.optimizers.SGD(lr=1.5), output_activation="sigmoid",
    metrics=[rounded_accuracy])

# Reconstruct real-valued codes
enc2, dec2, _, _ = train_autoencoder(
    30, X_train_enc1, X_valid_enc1, "mse", keras.optimizers.SGD(lr=0.05),
    output_activation="selu")

# Stack models, no fine tuning
stacked_ae = keras.models.Sequential([
    keras.layers.Flatten(input_shape=[28, 28]),
    enc1, enc2, dec2, dec1,
    keras.layers.Reshape([28, 28])
])

show_reconstructions(stacked_ae)
plt.show()

# Fine tune stacked model end to end
stacked_ae.compile(loss="binary_crossentropy",
                          optimizer=keras.optimizers.SGD(lr=0.1), metrics=[rounded_accuracy])
history = stacked_ae.fit(X_train, X_train, epochs=10,
                                validation_data=[X_valid, X_valid])
show_reconstructions(stacked_ae)
plt.show()

max_stars_repo_name: vipavlovic/pyprobml

max_issues_repo_name: vipavlovic/pyprobml

max_forks_repo_name: vipavlovic/pyprobml


text:
import numpy as np
import pickle
import torch
from torch.utils.data import Dataset

class YelpDataset(Dataset):
    def __init__(self, X_file, Y_file, num_train_lines, vocab, max_seq_len, update_vocab):
        self.raw_X = []
        self.X = []
        with open(X_file) as inf:
            cnt = 0
            for line in inf:
                words = line.rstrip().split()
                if max_seq_len > 0:
                    words = words[:max_seq_len]
                self.raw_X.append(words)
                if update_vocab:
                    for w in words:
                        vocab.add_word(w)
                # save lengths
                self.X.append(([vocab.lookup(w) for w in words], len(words)))
                cnt += 1
                if num_train_lines > 0 and cnt >= num_train_lines:
                    break

        self.max_seq_len = max_seq_len
        if isinstance(Y_file, str):
            self.Y = (torch.from_numpy(np.loadtxt(Y_file)) - 1).long()
        else:
            self.Y = Y_file
        if num_train_lines > 0:
            self.X = self.X[:num_train_lines]
            self.Y = self.Y[:num_train_lines]
        self.num_labels = 5
        # self.Y = self.Y.to(opt.device)
        assert len(self.X) == len(self.Y), 'X and Y have different lengths'
        print('Loaded Yelp dataset of {} samples'.format(len(self.X)))

    def __len__(self):
        return len(self.Y)

    def __getitem__(self, idx):
        return (self.X[idx], self.Y[idx])

    def set_max_seq_len(self, max_seq_len):
        self.X = [(x[0][:max_seq_len], min(x[1], max_seq_len)) for x in self.X]
        self.max_seq_len = max_seq_len

    def get_max_seq_len(self):
        if not hasattr(self, 'max_seq_len'):
            self.max_seq_len = max([x[1] for x in self.X])
        return self.max_seq_len

def get_yelp_datasets(vocab,
                      X_train_filename,
                      Y_train_filename,
                      num_train_lines,
                      X_test_filename,
                      Y_test_filename,
                      max_seq_len):
    train_dataset = YelpDataset(X_train_filename, Y_train_filename,
            num_train_lines, vocab, max_seq_len, update_vocab=True)
    valid_dataset = YelpDataset(X_test_filename, Y_test_filename,
            0, vocab, max_seq_len, update_vocab=True)
    return train_dataset, valid_dataset

max_stars_repo_name: forest-snow/adan

max_issues_repo_name: forest-snow/adan

max_forks_repo_name: forest-snow/adan


text:
"""
utility functions that support scBasset.
"""

import tensorflow as tf
import numpy as np
import anndata
from scipy import sparse
import h5py
from Bio import SeqIO
from scbasset.basenji_utils import *

###############################
# function for pre-processing #
###############################
def make_h5(
    input_ad,
    input_bed,
    input_fasta,
    out_file,
    seq_len=1344,
    train_ratio=0.9,
):
    """Preprocess to generate h5 for scBasset training.
    Args:
        input_ad:       anndata. the peak by cell matrix.
        input_bed:      bed file. genomic range of peaks.
        input_fasta:    fasta file. genome fasta. (hg19, hg38, mm10 etc.)
        out_file:       output file name.
        seq_len:        peak size to train on. default to 1344.
        train_ratio:    fraction of data used for training. default to 0.9.
    Returns:
        None.           Save a h5 file to 'out_file'. X: 1-hot encoded feature matrix.
                        Y: peak by cell matrix. train_ids: data indices used for train.
                        val_ids: data indices used for val. test_ids: data indices unused
                        during training, can be used for test.
    """
    # generate 1-hot-encoding from bed
    (seqs_dna, seqs_coords,) = make_bed_seqs(
        input_bed,
        fasta_file=input_fasta,
        seq_len=seq_len,
    )
    dna_array = [dna_1hot(x) for x in seqs_dna]
    dna_array = np.array(dna_array)
    ids = np.arange(dna_array.shape[0])
    np.random.seed(10)
    test_val_ids = np.random.choice(
        ids,
        int(len(ids) * (1 - train_ratio)),
        replace=False,
    )
    train_ids = np.setdiff1d(ids, test_val_ids)
    val_ids = np.random.choice(
        test_val_ids,
        int(len(test_val_ids) / 2),
        replace=False,
    )
    test_ids = np.setdiff1d(test_val_ids, val_ids)

    # generate binary peak*cell matrix
    ad = anndata.read_h5ad(input_ad)
    if sparse.issparse(ad.X):
        m = (np.array(ad.X.todense()).transpose() != 0) * 1
    else:
        m = (np.array(ad.X).transpose() != 0) * 1

    # save train_test_val splits
    f = h5py.File(out_file, "w")
    f.create_dataset(
        "X",
        data=dna_array,
        dtype="bool",
    )
    f.create_dataset("Y", data=m, dtype="int8")
    f.create_dataset(
        "train_ids",
        data=train_ids,
        dtype="int",
    )
    f.create_dataset(
        "val_ids",
        data=val_ids,
        dtype="int",
    )
    f.create_dataset(
        "test_ids",
        data=test_ids,
        dtype="int",
    )
    f.close()

################
# create model #
################
def make_model(
    bottleneck_size,
    n_cells,
    seq_len=1344,
    show_summary=True,
):
    """create keras CNN model.
    Args:
        bottleneck_size:int. size of the bottleneck layer.
        n_cells:        int. number of cells in the dataset. Defined the number of tasks.
        seq_len:        int. peak size used to train. Default to 1344.
        show_summary:   logical. Whether to print the model summary. Default to True.
    Returns:
        keras model:    object of keras model class.
    """
    sequence = tf.keras.Input(
        shape=(seq_len, 4),
        name="sequence",
    )
    current = sequence
    (current, reverse_bool,) = StochasticReverseComplement()(
        current
    )  # enable random rv
    current = StochasticShift(3)(current)  # enable random shift
    current = conv_block(
        current,
        filters=288,
        kernel_size=17,
        pool_size=3,
    )  # 1 cnn block
    current = conv_tower(
        current,
        filters_init=288,
        filters_mult=1.122,
        repeat=6,
        kernel_size=5,
        pool_size=2,
    )  # cnn tower
    current = conv_block(
        current,
        filters=256,
        kernel_size=1,
    )
    current = dense_block(
        current,
        flatten=True,
        units=bottleneck_size,
        dropout=0.2,
    )
    current = GELU()(current)
    current = final(
        current,
        units=n_cells,
        activation="sigmoid",
    )
    current = SwitchReverse()(
        [current, reverse_bool]
    )  # switch back, # this doesn't matter
    current = tf.keras.layers.Flatten()(current)
    model = tf.keras.Model(inputs=sequence, outputs=current)
    if show_summary:
        model.summary()
    return model

################################
# function for post-processing #
################################
def get_cell_embedding(model):
    """get cell embeddings from trained model"""
    return model.layers[-3].get_weights()[0].transpose()

def get_intercept(model):
    """get intercept from trained model"""
    return model.layers[-3].get_weights()[1]

def imputation_Y(X, model):
    """Perform imputation. Don't normalize for depth.
    Args:
        X:              feature matrix from h5.
        model:          a trained scBasset model.
    Returns:
        array:          a peak*cell imputed accessibility matrix. Sequencing depth
                        isn't corrected for.
    """
    Y_impute = model.predict(X)
    return Y_impute

# perform imputation. Depth normalized.
def imputation_Y_normalize(X, model, scale_method=None):
    """Perform imputation. Normalize for depth.
    Args:
        X:              feature matrix from h5.
        model:          a trained scBasset model.
    Returns:
        array:          a peak*cell imputed accessibility matrix. Sequencing depth corrected
                        for. scale_method=None, don't do any scaling of output. The raw
                        normalized output would have both positive and negative values.
                        scale_method="all_positive" scales the output by subtracting minimum value.
                        scale_method="sigmoid" scales the output by sigmoid transform.
    """
    new_model = tf.keras.Model(
        inputs=model.layers[0].input,
        outputs=model.layers[-4].output,
    )
    Y_pred = new_model.predict(X)
    w = model.layers[-3].get_weights()[0]
    accessibility_norm = np.dot(Y_pred.squeeze(), w)

    if scale_method == "all_positive":
        accessibility_norm = accessibility_norm - np.min(accessibility_norm)
    if scale_method == "sigmoid":
        accessibility_norm = np.divide(
            1,
            1 + np.exp(-accessibility_norm),
        )

    return accessibility_norm

def pred_on_fasta(fa, model):
    """Run a trained model on a fasta file.
    Args:
        fa:             fasta file to run on. Need to have a fixed size of 1344. Default
                        sequence size of trained model.
        model:          a trained scBasset model.
    Returns:
        array:          a peak*cell imputed accessibility matrix. Sequencing depth corrected for.
    """
    records = list(SeqIO.parse(fa, "fasta"))
    seqs = [str(i.seq) for i in records]
    seqs_1hot = np.array([dna_1hot(i) for i in seqs])
    pred = imputation_Y_normalize(seqs_1hot, model)
    return pred

def motif_score(tf, model, motif_fasta_folder):
    """score motifs for any given TF.
    Args:
        tf:             TF of interest. By default we only provide TFs to score in
                        https://storage.googleapis.com/scbasset_tutorial_data/Homo_sapiens_motif_fasta.tar.gz.
                        To score on additional motifs, follow make_fasta.R in the tarball
                        to create dinucleotide shuffled sequences with and without motifs of
                        interest.
        model:          a trained scBasset model.
        motif_fasta_folder: folder for dinucleotide shuffled sequences with and without any motif.
                        We provided motifs from CIS-BP/Homo_sapiens.meme downloaded from the
                        MEME Suite (https://meme-suite.org/meme/) in
                        https://storage.googleapis.com/scbasset_tutorial_data/Homo_sapiens_motif_fasta.tar.gz.
    Returns:
        array:          a vector for motif activity per cell. (cell order is the
                        same order as the model.)
    """
    fasta_motif = "%s/shuffled_peaks_motifs/%s.fasta" % (motif_fasta_folder, tf)
    fasta_bg = "%s/shuffled_peaks.fasta" % motif_fasta_folder

    pred_motif = pred_on_fasta(fasta_motif, model)
    pred_bg = pred_on_fasta(fasta_bg, model)
    tf_score = pred_motif.mean(axis=0) - pred_bg.mean(axis=0)
    tf_score = (tf_score - tf_score.mean()) / tf_score.std()
    return tf_score

max_stars_repo_name: calico/scBasset

max_issues_repo_name: calico/scBasset

max_forks_repo_name: calico/scBasset


text:
# -*- coding: utf-8 -*-

import glob
import numpy as np
import cv2
import matplotlib.image as mpimg
import matplotlib.pyplot as plt

def abs_sobel_thresh(img, orient='x', sobel_kernel=3, thresh=(0, 255)):
    # Convert to grayscale
    gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
    # Apply x or y gradient with the OpenCV Sobel() function
    # and take the absolute value
    if orient == 'x':
        abs_sobel = np.absolute(cv2.Sobel(gray, cv2.CV_64F, 1, 0, ksize=sobel_kernel))
    if orient == 'y':
        abs_sobel = np.absolute(cv2.Sobel(gray, cv2.CV_64F, 0, 1, ksize=sobel_kernel))
    # Rescale back to 8 bit integer
    scaled_sobel = np.uint8(255 * abs_sobel / np.max(abs_sobel))
    # Create a copy and apply the threshold
    binary_output = np.zeros_like(scaled_sobel)
    # Here I'm using inclusive (>=, <=) thresholds, but exclusive is ok too
    binary_output[(scaled_sobel >= thresh[0]) & (scaled_sobel <= thresh[1])] = 1

    # Return the result
    return binary_output

def mag_thresh(img, sobel_kernel=3, mag_thresh=(0, 255)):
    # Convert to grayscale
    gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
    # Take both Sobel x and y gradients
    sobelx = cv2.Sobel(gray, cv2.CV_64F, 1, 0, ksize=sobel_kernel)
    sobely = cv2.Sobel(gray, cv2.CV_64F, 0, 1, ksize=sobel_kernel)
    # Calculate the gradient magnitude
    gradmag = np.sqrt(sobelx ** 2 + sobely ** 2)
    # Rescale to 8 bit
    scale_factor = np.max(gradmag) / 255
    gradmag = (gradmag / scale_factor).astype(np.uint8)
    # Create a binary image of ones where threshold is met, zeros otherwise
    binary_output = np.zeros_like(gradmag)
    binary_output[(gradmag >= mag_thresh[0]) & (gradmag <= mag_thresh[1])] = 1

    # Return the binary image
    return binary_output

def dir_threshold(img, sobel_kernel=3, thresh=(0, np.pi / 2)):
    # Grayscale
    gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
    # Calculate the x and y gradients
    sobelx = cv2.Sobel(gray, cv2.CV_64F, 1, 0, ksize=sobel_kernel)
    sobely = cv2.Sobel(gray, cv2.CV_64F, 0, 1, ksize=sobel_kernel)
    # Take the absolute value of the gradient direction,
    # apply a threshold, and create a binary image result
    absgraddir = np.arctan2(np.absolute(sobely), np.absolute(sobelx))
    binary_output = np.zeros_like(absgraddir)
    binary_output[(absgraddir >= thresh[0]) & (absgraddir <= thresh[1])] = 1

    # Return the binary image
    return binary_output

if __name__ == "__main__":
    images = glob.glob('./output_images/undistorted/*.jpg')
    for fname in images:
        image = mpimg.imread(fname)
        binary_output = abs_sobel_thresh(image, orient='x', sobel_kernel=3, thresh=(10, 150))
        plt.imshow(binary_output, cmap='gray')
        plt.show()

    images = glob.glob('./output_images/undistorted/*.jpg')
    for fname in images:
        image = mpimg.imread(fname)
        binary_output = abs_sobel_thresh(image, orient='y', sobel_kernel=3, thresh=(0, 25))
        plt.imshow(binary_output, cmap='gray')
        plt.show()

max_stars_repo_name: jo-ny/CarND-Advanced-Lane-Lines

max_issues_repo_name: jo-ny/CarND-Advanced-Lane-Lines

max_forks_repo_name: jo-ny/CarND-Advanced-Lane-Lines


text:
# Utilities in this file are all adopted from the implementation of SimGNN: https://github.com/yunshengb/SimGNN

import os
import pickle
import random
import re
from glob import glob
from os.path import basename
from os.path import isfile

import networkx as nx
import numpy as np
from sklearn.preprocessing import OneHotEncoder

class SelfShuffleList(object):
    def __init__(self, li):
        assert (type(li) is list and list)
        self.li = li  # be careful! not a deep copy!
        self.idx = 0
        self._shuffle()  # shuffle at the beginning

    def get_next_one(self):
        if self.idx < len(self.li):
            rtn = self.li[self.idx]
            self.idx += 1
            return rtn
        else:
            self.idx = 0
            self._shuffle()
            return self.li[self.idx]

    def _shuffle(self):
        random.Random(123).shuffle(self.li)

class Data(object):
    def __init__(self, save_folder):
        print(save_folder)
        assert os.path.exists(save_folder)
        self.save_folder = save_folder

    def get_save_folder(self):
        return self.save_folder

    def get_dataset_folder_name(self):
        raise NotImplementedError

    def get_max_number_of_nodes(self):
        raise NotImplementedError

class AIDSData(Data):
    def __init__(self, save_folder, train_or_test):

        super(AIDSData, self).__init__(save_folder=save_folder)

        data_dir = '{}/{}/{}'.format(self.get_save_folder(), self.get_dataset_folder_name(), train_or_test)
        self.graphs = []
        self.graphs = iterate_get_graphs(data_dir)
        print('Loaded {} graphs from {}'.format(len(self.graphs), data_dir))

        number_of_nodes = []
        for g in self.graphs:
            number_of_nodes.append(g.number_of_nodes())
        self.max_number_of_nodes = max(number_of_nodes)
        print('The max number of nodes among all graphs of {} is # {}'.format(self.get_dataset_folder_name(),
                                                                              self.get_max_number_of_nodes()))

        if 'nef' in self.get_dataset_folder_name():
            print('Removing edge features -> valence feature')
            for g in self.graphs:
                for n1, n2, d in g.edges(data=True):
                    d.pop('valence', None)

    def get_dataset_folder_name(self):
        return 'AIDS700nef'

    def get_gids(self):
        return [g.graph['gid'] for g in self.graphs]

    def get_max_number_of_nodes(self):
        return self.max_number_of_nodes

class LinuxData(Data):
    def __init__(self, save_folder, train_or_test):
        super(LinuxData, self).__init__(save_folder=save_folder)
        # Load all graphs
        data_dir = '{}/{}/{}'.format(self.get_save_folder(), self.get_dataset_folder_name(), train_or_test)
        self.graphs = []
        self.graphs = iterate_get_graphs(data_dir)
        print('Loaded {} graphs from {}'.format(len(self.graphs), data_dir))

        number_of_nodes = []
        for g in self.graphs:
            number_of_nodes.append(g.number_of_nodes())
        self.max_number_of_nodes = max(number_of_nodes)
        print('The max number of nodes among all graphs of {} is # {}'.format(self.get_dataset_folder_name(),
                                                                              self.get_max_number_of_nodes()))

    def get_dataset_folder_name(self):
        return 'LINUX'

    def get_gids(self):
        return [g.graph['gid'] for g in self.graphs]

    def get_max_number_of_nodes(self):
        return self.max_number_of_nodes

class SiameseDataSet(object):
    def __init__(self, data_save_folder, data_set_name, validation_ratio, node_feat_name, node_encoder_name):

        self.data_save_folder = data_save_folder
        self.data_set_name = data_set_name
        self.validation_ratio = validation_ratio
        self.node_feat_name = node_feat_name
        self.node_encoder_name = node_encoder_name

        # Load train and validation graphs
        print('\nprepare the train/validation data set and split')
        orig_train_data = self.load_data_set(data=self.data_set_name, train_or_test='train')
        self.train_val_gs = orig_train_data.graphs

        self.train_gs, self.val_gs = self._train_val_split(orig_train_data)
        self.train_val_max_number_of_nodes = orig_train_data.get_max_number_of_nodes()
        print('\t\ttrain_gs.len={} and val_gs.len={}, and their max number of nodes = {}'.format(len(self.train_gs),
                                                                                len(self.val_gs),
                                                                                self.train_val_max_number_of_nodes))
        print('\nprepare the testing data set')
        orig_test_data = self.load_data_set(data=self.data_set_name, train_or_test='test')
        self.test_gs = orig_test_data.graphs
        self.test_max_number_of_nodes = orig_test_data.get_max_number_of_nodes()

        print('\t\ttest_gs.len={}, and their max number of nodes = {}'.format(len(self.test_gs),
                                                                              self.test_max_number_of_nodes))

        self.node_feat_encoder = self._get_node_feature_encoder(orig_train_data.graphs + self.test_gs)

    def load_data_set(self, data, train_or_test):
        if data.lower() == 'aids700nef':
            return AIDSData(save_folder=self.data_save_folder, train_or_test=train_or_test)
        elif data.lower() == 'linux':
            return LinuxData(save_folder=self.data_save_folder, train_or_test=train_or_test)
        else:
            raise RuntimeError('!!! Not recognized data set of {}'.format(data))

    def input_dim(self):
        return self.node_feat_encoder.input_dim()

    def _train_val_split(self, orig_train_data):
        if self.validation_ratio < 0 or self.validation_ratio > 1:
            raise RuntimeError('ratio of validation {} must be in [0, 1]'.format(self.validation_ratio))
        gs = orig_train_data.graphs
        sp = int(len(gs) * (1 - self.validation_ratio))
        train_graphs = gs[0:sp]
        valid_graphs = gs[sp:]
        self._check_graphs_num(train_graphs, 'train')
        self._check_graphs_num(valid_graphs, 'validation')
        return train_graphs, valid_graphs

    @staticmethod
    def _check_graphs_num(graphs, label):
        if len(graphs) <= 2:
            raise RuntimeError('Insufficient {} graphs {}'.format(label, len(graphs)))

    def _get_node_feature_encoder(self, graphs):
        if self.node_encoder_name == 'OneHot':
            return NodeFeatureOneHotEncoder(graphs, self.node_feat_name)
        elif 'constant' in self.node_encoder_name:
            return NodeFeatureConstantEncoder(graphs, self.node_feat_name, self.node_encoder_name)
        else:
            raise RuntimeError('Unknown node_feat_encoder {}'.format(self.node_feat_encoder))

class NodeFeatureEncoder(object):
    def encode(self, g):
        raise NotImplementedError()

    def input_dim(self):
        raise NotImplementedError()

class NodeFeatureOneHotEncoder(NodeFeatureEncoder):
    def __init__(self, gs, node_feat_name):
        self.node_feat_name = node_feat_name
        inputs_set = set()
        for g in gs:
            inputs_set = inputs_set | set(self._node_feat_dic(g).values())
        self.feat_idx_dic = {feat: idx for idx, feat in enumerate(inputs_set)}

        self.oe = OneHotEncoder().fit(np.array(list(self.feat_idx_dic.values())).reshape(-1, 1))

    def encode(self, g):
        node_feat_dic = self._node_feat_dic(g)
        temp = [self.feat_idx_dic[node_feat_dic[n]] for n in g.nodes()]
        return self.oe.transform(np.array(temp).reshape(-1, 1)).toarray()

    def input_dim(self):
        return self.oe.transform([[0]]).shape[1]

    def _node_feat_dic(self, g):
        return nx.get_node_attributes(g, self.node_feat_name)

class NodeFeatureConstantEncoder(NodeFeatureEncoder):

    def __init__(self, _, node_feat_name=None, node_encoder_name='constant_1_2'):
        assert (node_feat_name is None)
        self.input_dim_ = int(node_encoder_name.split('_')[1])
        self.const = float(node_encoder_name.split('_')[2])
        print('A constant feature encoder where the input_dim is {} and input feature is {}'.format(self.input_dim_,
                                                                                self.const))

    def encode(self, g):
        rnt = np.full((g.number_of_nodes(), self.input_dim_), self.const)
        return rnt

    def input_dim(self):
        return self.input_dim_

def sorted_nicely(l):
    def try_int(s):
        try:
            return int(s)
        except:
            return s

    def alphanum_key(s):
        return [try_int(c) for c in re.split('([0-9]+)', s)]

    return sorted(l, key=alphanum_key)

def iterate_get_graphs(directory):
    graphs = []
    for file in sorted_nicely(glob(directory + '/*.gexf')):

        gid = int(basename(file).split('.')[0])

        g = nx.read_gexf(file)
        g.graph['gid'] = gid
        graphs.append(g)
        if not nx.is_connected(g):
            raise RuntimeError('{} not connected'.format(gid))
    return graphs

def normalized_exp_kernel_ged(d, g1, g2, scale):
    normalized = 2 * d / (g1.number_of_nodes() + g2.number_of_nodes())
    return np.exp(-scale * normalized)

def proc_filepath(filepath):
    if type(filepath) is not str:
        raise RuntimeError('Did you pass a file path to this function?')
    ext = '.pickle'
    if ext not in filepath:
        filepath += ext
    return filepath

def load(filepath):
    filepath = proc_filepath(filepath)
    if isfile(filepath):
        with open(filepath, 'rb') as handle:
            return pickle.load(handle)
    else:
        return None

def top_k_ids(query, qid, k, inclusive, rm=0):
    sort_id_mat = np.argsort(query, kind='mergesort')[:, ::-1]
    _, n = sort_id_mat.shape
    if k < 0 or k >= n:
        raise RuntimeError('Invalid k {}'.format(k))
    if not inclusive:
        return sort_id_mat[qid][:k]

    while k < n:
        cid = sort_id_mat[qid][k - 1]
        nid = sort_id_mat[qid][k]
        if abs(query[qid][cid] - query[qid][nid]) <= rm:
            k += 1
        else:
            break
    return sort_id_mat[qid][:k], k

def computing_precision_ks(trues, predictions, ks, inclusive=True, rm=0):
    assert trues.shape == predictions.shape
    m, n = trues.shape

    precision_ks = np.zeros((m, len(ks)))
    inclusive_final_true_ks = np.zeros((m, len(ks)))
    inclusive_final_pred_ks = np.zeros((m, len(ks)))

    for i in range(m):

        for k_idx, k in enumerate(ks):
            assert (type(k) is int and 0 < k < n)
            true_ids, true_k = top_k_ids(trues, i, k, inclusive, rm)
            pred_ids, pred_k = top_k_ids(predictions, i, k, inclusive, rm)
            precision_ks[i][k_idx] = min(len(set(true_ids).intersection(set(pred_ids))), k) / k
            inclusive_final_true_ks[i][k_idx] = true_k
            inclusive_final_pred_ks[i][k_idx] = pred_k
    return np.mean(precision_ks, axis=0), np.mean(inclusive_final_true_ks, axis=0), np.mean(inclusive_final_pred_ks,
                                                                                axis=0)

def get_root_path():
    from os.path import dirname, abspath
    return dirname(dirname(abspath(__file__)))

def get_save_path():
    return get_root_path() + '/save'

def save_pkl(obj, handle):
    import pickle
    pickle.dump(obj, handle, protocol=pickle.HIGHEST_PROTOCOL)

def save(filepath, obj):
    with open(proc_filepath(filepath), 'wb') as handle:
        save_pkl(obj, handle)

class DistanceCalculator(object):

    def __init__(self, root_folder, data_name):

        self.sfn = '{}/{}_{}_{}_gidpair_dist_map'.format(root_folder, data_name.lower(), 'ged', 'astar')
        self.gid_pair_dist_map = load(self.sfn)
        print(type(self.gid_pair_dist_map))
        if not self.gid_pair_dist_map:
            raise RuntimeError('{} -> distance map does not exists!'.format(self.sfn))
        else:
            print('Loaded dist map from {} with {} entries'.format(self.sfn, len(self.gid_pair_dist_map)))

    def calculate_distance_btw_pairs(self, g1, g2, scale):
        gid1 = g1.graph['gid']
        gid2 = g2.graph['gid']
        pair = (gid1, gid2)
        d = self.gid_pair_dist_map.get(pair)
        if d is None:
            raise RuntimeWarning('{} distance is None'.format(d))
        return d, normalized_exp_kernel_ged(d, g1, g2, scale)

class Result(object):
    """
    The result object loads and stores the ranking result of a model
        for evaluation.
        Terminology:
            rtn: return value of a function.
            m: # of queries.
            n: # of database graphs.
    """

    def model(self):
        """
        :return: The model name.
        """
        return self.model_

    def m_n(self):
        return self.dist_sim_mat(norm=False).shape

    def dist_or_sim(self):
        raise NotImplementedError()

    def dist_sim_mat(self, norm):
        """
        Each result object stores either a distance matrix
            or a similarity matrix. It cannot store both.
        :param norm:
        :return: either the distance matrix or the similarity matrix.
        """
        raise NotImplementedError()

    def dist_sim(self, qid, gid, norm):
        """
        :param qid: query id (0-indexed).
        :param gid: database graph id (0-indexed) (NOT g.graph['gid']).
        :param norm:
        :return: (metric, dist or sim between qid and gid)
        """
        raise NotImplementedError()

    def sim_mat(self, sim_kernel, yeta, scale, norm):
        raise NotImplementedError()

    def top_k_ids(self, qid, k, norm, inclusive, rm):
        """
        :param qid: query id (0-indexed).
        :param k:
        :param norm:
        :param inclusive: whether to be tie inclusive or not.
            For example, the ranking may look like this:
            7 (sim_score=0.99), 5 (sim_score=0.99), 10 (sim_score=0.98), ...
            If tie inclusive, the top 1 results are [7, 9].
            Therefore, the number of returned results may be larger than k.
            In summary,
                len(rtn) == k if not tie inclusive;
                len(rtn) >= k if tie inclusive.
        :param rm:
        :return: for a query, the ids of the top k database graph
        ranked by this model.
        """
        sort_id_mat = self.sort_id_mat(norm)
        _, n = sort_id_mat.shape
        if k < 0 or k >= n:
            raise RuntimeError('Invalid k {}'.format(k))
        if not inclusive:
            return sort_id_mat[qid][:k]
        # Tie inclusive.
        dist_sim_mat = self.dist_sim_mat(norm)
        while k < n:
            cid = sort_id_mat[qid][k - 1]
            nid = sort_id_mat[qid][k]
            if abs(dist_sim_mat[qid][cid] - dist_sim_mat[qid][nid]) <= rm:
                k += 1
            else:
                break
        return sort_id_mat[qid][:k]

    def ranking(self, qid, gid, norm, one_based=True):
        """
        :param qid: query id (0-indexed).
        :param gid: database graph id (0-indexed) (NOT g.graph['gid']).
        :param norm:
        :param one_based: whether to return the 1-based or 0-based rank.
            True by default.
        :return: for a query, the rank of a database graph by this model.
        """
        # Assume self is ground truth.
        sort_id_mat = self.sort_id_mat(norm)
        finds = np.where(sort_id_mat[qid] == gid)
        assert (len(finds) == 1 and len(finds[0]) == 1)
        fid = finds[0][0]
        # Tie inclusive (always when find ranking).
        dist_sim_mat = self.dist_sim_mat(norm)
        while fid > 0:
            cid = sort_id_mat[qid][fid]
            pid = sort_id_mat[qid][fid - 1]
            if dist_sim_mat[qid][pid] == dist_sim_mat[qid][cid]:
                fid -= 1
            else:
                break
        if one_based:
            fid += 1
        return fid

    def classification_mat(self, thresh_pos, thresh_neg,
                           thresh_pos_sim, thresh_neg_sim, norm):
        raise NotImplementedError()

    def time(self, qid, gid):
        raise NotImplementedError()

    def time_mat(self):
        raise NotImplementedError()

    def mat(self, metric, norm):
        raise NotImplementedError()

    def sort_id_mat(self, norm):
        """
        :param norm:
        :return: a m by n matrix representing the ranking result.
            rtn[i][j]: For query i, the id of the j-th most similar
                       graph ranked by this model.
        """
        raise NotImplementedError()

    def ranking_mat(self, norm, one_based=True):
        """
        :param norm:
        :param one_based:
        :return: a m by n matrix representing the ranking result.
                 Note it is different from sort_id_mat.
            rtn[i][j]: For query i, the ranking of the graph j.
        """
        m, n = self.m_n()
        rtn = np.zeros((m, n))
        for i in range(m):
            for j in range(n):
                rtn[i][j] = self.ranking(i, j, norm, one_based=one_based)
        return rtn

class DistanceModelResult(Result):
    def __init__(self, dataset, model, result_folder, dist_mat, row_graphs, col_graphs, scale):
        self.dataset = dataset
        self.model_ = model
        self.result_folder = result_folder
        if dist_mat is not None:
            self.dist_mat_ = dist_mat
            self.time_mat_ = None
        else:
            self.dist_mat_ = self._load_result_mat(self.dist_metric(), self.model_, len(row_graphs), len(col_graphs))
            print('self.dist_mat_.shape = ', self.dist_mat_.shape)
        self.dist_norm_mat_ = np.copy(self.dist_mat_)

        m, n = self.dist_mat_.shape
        assert (m == len(row_graphs))
        assert (n == len(col_graphs))
        for i in range(m):
            for j in range(n):
                self.dist_norm_mat_[i][j] = normalized_exp_kernel_ged(self.dist_mat_[i][j], row_graphs[i],
                                                                      col_graphs[j], scale)
        self.sort_id_mat_ = np.argsort(self.dist_mat_, kind='mergesort')
        self.dist_norm_sort_id_mat_ = np.argsort(self.dist_norm_mat_, kind='mergesort')

    def dist_metric(self):
        raise NotImplementedError()

    def dist_mat(self, norm):
        return self._select_dist_mat(norm)

    def dist_sim_mat(self, norm):
        return self.dist_mat(norm)

    def dist_or_sim(self):
        return 'dist'

    def dist_sim(self, qid, gid, norm):
        return self.dist_metric(), self._select_dist_mat(norm)[qid][gid]

    def time(self, qid, gid):
        return self.time_mat_[qid][gid]

    def time_mat(self):
        return self.time_mat_

    def mat(self, metric, norm):
        if metric == self.dist_metric():
            return self._select_dist_mat(norm)
        elif metric == 'time':
            return self.time_mat_
        else:
            raise RuntimeError('Unknown metric {} for model {}'.format(metric, self.model_))

    def sort_id_mat(self, norm):
        return self._select_sort_id_mat(norm)

    def _load_result_mat(self, metric, model, m, n):
        file_p = self.result_folder + '/{}/{}/{}_{}_mat_{}_{}_*.npy'.format(self.dataset, metric, self.dist_metric(),
                                                                            metric, self.dataset, model)
        print('load_result_mat\nfile_p: ', file_p)
        li = glob(file_p)
        if not li:
            if 'astar' in model:
                if self.dataset != 'imdbmulti':
                    raise RuntimeError('Not imdbmulti and no astar results!')
                return self._load_merged_astar_from_other_three(metric, m, n)
            else:
                raise RuntimeError('No results found {}'.format(file_p))
        file = self._choose_result_file(li, m, n)
        # file = '../results' + '/{}/{}/{}'.format(self.dataset, metric, 'ged_ged_mat_aids700nef_astar_2018-06-27T19_38_42_qilin_10cpus.npy')
        return np.load(file)

    def _load_merged_astar_from_other_three(self, metric, m, n):
        print("_load_merged_astar_from_other_three: beam80, hungarian, vj")
        self.model_ = 'merged_astar'
        merge_models = ['beam80', 'hungarian', 'vj']
        dms = [self._load_result_mat(self.dist_metric(), model, m, n)
               for model in merge_models]
        tms = [self._load_result_mat('time', model, m, n)
               for model in merge_models]
        for i in range(len(merge_models)):
            assert (dms[i].shape == (m, n))
            assert (tms[i].shape == (m, n))
        rtn = np.zeros((m, n))
        for i in range(m):
            for j in range(n):
                d_li = [dm[i][j] for dm in dms]
                d_min = np.min(d_li)
                if metric == self.dist_metric():
                    rtn[i][j] = d_min
                else:
                    assert (metric == 'time')
                    d_argmin = np.argmin(d_li)
                    rtn[i][j] = tms[d_argmin][i][j]
        return rtn

    def _choose_result_file(self, files, m, n):
        cands = []
        for file in files:
            temp = np.load(file)
            if temp.shape == (m, n):
                cands.append(file)
                if 'qilin' in file:
                    print(file)
                    return file
        if cands:
            return cands[0]
        raise RuntimeError('No result files in {}'.format(files))  # TODO: smart choice and cross-checking

    def _select_dist_mat(self, norm):
        return self.dist_norm_mat_ if norm else self.dist_mat_

    def _select_sort_id_mat(self, norm):
        return self.dist_norm_sort_id_mat_ if norm else self.sort_id_mat_

class PairwiseGEDModelResult(DistanceModelResult):
    def dist_metric(self):
        return 'ged'

def load_result(dataset, model, result_folder="../GEDResults", sim=None, sim_mat=None, dist_mat=None, row_graphs=None,
                col_graphs=None, time_mat=None, model_info=None, scale=None):
    if 'beam' in model or model in ['astar', 'hungarian', 'vj']:
        return PairwiseGEDModelResult(dataset, model, result_folder, dist_mat, row_graphs, col_graphs, scale)
    elif 'siamese' in model:
        raise NotImplementedError
    else:
        raise RuntimeError('Unknown model {}'.format(model))

max_stars_repo_name: kleincup/MGMN

max_issues_repo_name: kleincup/MGMN

max_forks_repo_name: kleincup/MGMN

reason: import numpy,import networkx


text:
import numpy as np
import os, sys
import copy
import matplotlib.pyplot as plt
import palettable as pal
from matplotlib import cm

#from visualize import *
import mcmc

#set seed so we get predictable results
np.random.seed( 12 )
from mpl_toolkits.mplot3d import Axes3D

#Profiling / Timing
sys.path.insert(0, '../tools')
from timer import Timer

class Params:
    mins = None
    maxs = None
    lens = None

def randab(a, b):
    return a + (b-a)*np.random.rand()

#random 2d location
def randLoc(params):
    x = randab(params.mins[0], params.maxs[0])
    y = randab(params.mins[1], params.maxs[1])
    return (x,y)

#random 3D direction in terms of angles phi and theta
def randUnitSphere():
    vphi = randab( 0.0, 2.0*np.pi ) #azimuth
    vthe = randab(-np.pi/2.0, np.pi/2.0 ) #latitude
    return (vphi, vthe)

def sph2cart(vr, vphi, vthe):
    vx = vr * np.sin(vthe) * np.cos(vphi)
    vy = vr * np.sin(vthe) * np.sin(vphi)
    vz = vr * np.cos(vthe)
    return vx, vy, vz

def pol2cart(vr, vphi, vthe):
    vx = vr * np.cos(vphi)
    vy = vr * np.sin(vphi)
    return vx, vy, 0.0

def randVel(vabs):
    (vphi, vthe) = randUnitSphere()
    vx, vy, vz = sph2cart(vabs, vphi, vthe)
    #vx, vy, vz = pol2cart(vabs, vphi, vthe)
    return vx, vy, vz

def unitVecX(vec):
    return np.array([ 0.0, vec[0] ])
def unitVecY(vec):
    return np.array([ 0.0, vec[1] ])
def unitVecZ(vec):
    return np.array([ 0.0, vec[2] ])

def norm(vec):
    return np.sqrt( vec[0]**2 + vec[1]**2 + vec[2]**2 )

#unit vector into direction a x b
def uCross(vecA, vecB):
    vecC = np.cross(vecA, vecB)
    return vecC / norm(vecC)

def printVec(vec, name):
    print " {}: ({}, {}, {}) |v| = {} ".format(name, vec[0], vec[1], vec[2], norm(vec))

########################################

# Monte Carlo Compton scattering according to Sobol 1977
def comptonScatter(e, p, ax, plot=True):

    # notes about units
    # p.hv = hv/m_e c^2
    # e.v = v/c

    beta = np.array([ e.vx(), e.vy(), e.vz() ]) / e.beta() #electron unit vector
    omeg = np.array([ p.vx(), p.vy(), p.vz() ])         #photon unit vector
    theta = np.arccos( np.dot(beta, omeg) ) #angle between electron and photon

    #printVec(beta, "beta")
    #printVec(omeg, "omeg")
    #print "thet:", theta

    #rotate to scattering plane (i,j,k)
    kvec = np.array([ 0.0, -1.0, 0.0] ) #this points along r, here we select it arbitrarily
    #kvec = omeg #another choice so we are in direction of photon
    jvec = uCross(beta, omeg)
    ivec = uCross(kvec, jvec)
    M = np.array([ ivec, jvec, kvec ])

    #printVec(ivec, "i")
    #printVec(jvec, "j")
    #printVec(kvec, "k")

    ##################################################
    #print M
    #print np.linalg.inv(M)

    ##################################################

    #unit vector of electron in scattering coordinates
    v0  = np.matmul(M, np.array([ e.vx(), e.vy(), e.vz()] )) #rotate electron velocity to scattering plane (i,k)
    mu  = v0[0]*np.sin(theta) +  v0[2]*np.cos(theta)
    rho = np.sqrt( v0[0]**2 + v0[1]**2 )

    #printVec(v0, "v0")

    #Compton parameters
    x = 2.0 * p.hv() * e.gamma() * (1.0 - mu*e.beta() )
    y = x/2.0

    #print "Compton x: {}".format(x)

    #additional scattering angles (v0, w0, t0) define a frame of reference
    w0 = np.array([ v0[1], -v0[0], 0])/rho
    t0 = np.array([ v0[0]*v0[2], v0[1]*v0[2], -rho**2])/rho

    #printVec(w0, "w0")
    #printVec(t0, "t0")

    #ax.plot( unitVecX(v0), unitVecY(v0), unitVecZ(v0), label='v0' )
    #ax.plot( unitVecX(w0), unitVecY(w0), unitVecZ(w0), label='w0' )
    #ax.plot( unitVecX(t0), unitVecY(t0), unitVecZ(t0), label='t0' )

    #scatter
    done = False
    OmegaOmegap = 0.0
    while not(done):
        z1 = np.random.rand()
        z2 = np.random.rand()
        z3 = np.random.rand()

        mup  = (e.beta() + 2.0*z1 - 1.0)/(1.0 + e.beta()*(2.0*z1 - 1.0))
        phip = 2.0*np.pi*z2

        OmegaOmegap = mu*mup - np.sqrt( 1.0 - mup**2) * ( rho*np.sin(phip)*np.cos(theta)
                - (1.0/rho)*(v0[1]*np.cos(phip) + v0[0]*v0[2]*np.sin(phip))*np.sin(theta))

        yp = y / (1.0 + p.hv()*(1.0 - OmegaOmegap))/(e.gamma() * (1.0 - mup*e.beta()))

        Y = yp/y + (yp/y)**3 + (yp/y)**2 *( (1.0/yp - 1.0/y)**2 - 2.0*( 1.0/yp - 1.0/y) )

        if Y > 2.0*z3:
            done = True
    #now we have scattered successfully

    #new energy
    hvp = yp / ( e.gamma()*(1.0 - mup*e.beta()) )
    #print " energy shift: {}".format( hvp/p.hv() )

    #new direction
    Omegap_ijk = mup*v0 + np.sqrt(1.0-mup**2)*( w0*np.cos(phip) + t0*np.sin(phip) )
    Omegap = np.dot( np.linalg.inv(M), Omegap_ijk )
    Omegap = Omegap / np.sqrt( Omegap[0]**2 + Omegap[1]**2 + Omegap[2]**2 ) #normalize

    if plot:
        ax.plot( unitVecX(Omegap), unitVecY(Omegap), unitVecZ(Omegap), alpha=0.2, color='red')#, label='Omegap' )

    phs = mcmc.photon( hvp, Omegap[0], Omegap[1], Omegap[2] )

    return e, phs

def visualize(ax, slab, params):

    xs = slab.xloc
    ys = slab.yloc
    zs = slab.zloc

    ax.plot(xs, ys, zs,
            color='black',
            marker='.'
            )

    return

if __name__ == "__main__":

    #fig = plt.figure()
    #ax = fig.add_subplot(111, projection='3d')

    ##################################################
    # set up plotting and figure
    plt.fig = plt.figure(1, figsize=(8,4))
    plt.rc('font', family='serif', size=12)
    plt.rc('xtick')
    plt.rc('ytick')

    gs = plt.GridSpec(1, 2)
    gs.update(hspace = 0.5)

    axs = []
    axs.append( plt.subplot(gs[0]) )
    axs.append( plt.subplot(gs[1], projection='3d') )
    #axs.append( plt.subplot(gs[0], projection='3d') )

    ##################################################
    # Normalize 3d plot to make aspect ratio unity
    #MAX = 1.0
    #for direction in (-1, 1):
    #    for point in np.diag(direction * MAX * np.array([1,1,1])):
    #        #axs[1].plot([point[0]], [point[1]], [point[2]], 'w')
    #        #axs[2].plot([point[0]], [point[1]], [point[2]], 'w')
    #        #axs[0].plot([point[0]], [point[1]], [point[2]], 'w')

    ##################################################
    # set-up grid
    # xy
    params = Params()
    params.mins = [ 0.0, 0.0, 0.0 ]
    params.maxs = [ 1.0, 1.0, 1.0 ]
    params.lens = [ 1.0, 1.0, 1.0 ]

    ###################################################
    #test scattering for single case
    print "Compton scatter single interaction..."

    #create electron
    (vx, vy, vz) = (0.8, 0.0, 0.0) #x-dir electron
    e = mcmc.electron()
    e.loadVelComponents(vx, vy, vz)
    print "target electron with beta: {} gamma: {}".format(e.beta(), e.gamma() )

    #create photon
    (vx, vy, vz) = randVel(1.0) #direction on unit sphere
    ph = mcmc.photon( 1.0, vx, vy, vz )

    #visualize starting point of collision
    beta = np.array( [e.vx(), e.vy(), e.vz()]) / e.beta()  #electron unit vector
    print e.vx(), e.beta(), beta
    omeg = np.array([ ph.vx(), ph.vy(), ph.vz() ])      #photon unit vector
    axs[1].plot( unitVecX(-beta), unitVecY(-beta), unitVecZ(-beta), linestyle='dashed', color='black', linewidth=3.0)
    axs[1].plot( unitVecX(-omeg), unitVecY(-omeg), unitVecZ(-omeg), linestyle='dotted', color='black', linewidth=4.0)

    #es, ps = comptonScatter(e, ph, axs[1], plot=True)

    timer = Timer(["scatter", "cppCS", "bucket"])
    timer.start("scatter")

    for sc in range(500):
        es, ps = comptonScatter(e, ph, axs[1], plot=False)
        timer.lap("scatter")
    timer.stats("scatter")

    btmp = mcmc.photonBucket()
    slab = mcmc.Slab( btmp )

    timer.start("cppCS")
    for sc in range(500):
        (ph2, el2) = slab.comptonScatter(ph, e)
        timer.lap("cppCS")
    timer.stats("cppCS")

    #plt.savefig("comptonFan.png")
    #plt.show()
    #sys.exit()

    ##################################################
    # now Monte Carlo scatter the whole bucket
    bucket = mcmc.photonBucket()
    print "created bucket ({})".format( bucket.size() )

    #pour isotropic photons to the bucket
    timer.start("bucket")
    for i in range(1000000):
        (vx, vy, vz) = randVel(1.0) #direction on unit sphere
        #(vx, vy, vz) = (0.8, 0.0, 0.0) #x-dir electron
        E = 0.01 # E = h\nu

        ph = mcmc.photon( E, vx, vy, vz )
        bucket.push_back( ph )
    timer.stop("bucket")
    print "loaded bucket with {} photons".format( bucket.size() )
    timer.stats("bucket")

    ##################################################
    xmin = 0.5
    xmax = 100.0
    axs[0].set_xscale('log')
    axs[0].set_yscale('log')

    axs[0].set_xlim(xmin, xmax)
    axs[0].set_ylim(0.2, 1.1)
    axs[0].minorticks_on()

    e = mcmc.electron() # target electron

    #step in time
    xs = np.zeros( bucket.size() )
    for lap in range(3):
        print("lap : {}").format(lap)
        for i in range(bucket.size() ):

            (vx, vy, vz) = randVel( 0.8 ) #draw random electron
            #(vx, vy, vz) = (0.8, 0.0, 0.0) #x-dir electron

            e.loadVelComponents(vx, vy, vz)
            p = bucket.get(i)

            #el2, ph2 = comptonScatter(e, p, axs[1], plot=False)
            (ph2, el2) = slab.comptonScatter(p, e)

            #print p.hv()
            xs[i] = ph2.hv()
            bucket.replace(i, ph2)

        #hist, edges = np.histogram(xs, np.linspace(xmin, xmax, 20))
        hist, edges = np.histogram(xs, np.logspace(np.log10(xmin), np.log10(xmax), 20))
        hist = 1.0 * hist / hist.max()
        axs[0].plot(edges[:-1], hist )#, "k-")

    #visualize_data(axs[0], mesh, params)
    plt.savefig("compton.png")

max_stars_repo_name: Krissmedt/imprunko

max_issues_repo_name: Krissmedt/imprunko

max_forks_repo_name: Krissmedt/imprunko


text:
import scipy.spatial as sns

class AdaptiveTree():
  """
  """
  def __init__():
    pass

  def Load(self, dataset, embedding_function):
    """
    """
    pass

  def Build(self, sample_function):
    """
    """
    pass

  def Sample(self, location, n, r):
    """
    """
    pass

  def UpdateEmbeddings(self, embeddings, gradient):
    """
    """
    pass

  def Embeddings(self, embeddings):
    """
    """
    pass

max_stars_repo_name: LunarEngineer/adaptive-kd-tree

max_issues_repo_name: LunarEngineer/adaptive-information-tree

max_forks_repo_name: LunarEngineer/adaptive-information-tree

reason: import scipy


text:
"""
Standard Regression model
-------------------------
"""
import numpy as np
import pandas as pd

from typing import Union
from ..logging import get_logger
from .regression_model import RegressionModel
from sklearn.linear_model import LinearRegression

logger = get_logger(__name__)

class LinearRegressionModel(RegressionModel):

    def __init__(self,
                 lags: Union[int, list] = None,
                 lags_exog: Union[int, list, bool] = None,
                 **kwargs):
        """
        Simple wrapper for the linear regression model in scikit-learn, LinearRegression().

        Parameters
        ----------
        lags : Union[int, list]
            Number of lagged target values used to predict the next time step. If an integer is given
            the last `lags` lags are used (inclusive). Otherwise a list of integers with lags is required.
        lags_exog : Union[int, list, bool]
            Number of lagged exogenous values used to predict the next time step. If an integer is given
            the last `lags_exog` lags are used (inclusive). Otherwise a list of integers with lags is required.
            If True `lags` will be used to determine lags_exog. If False, the values of all exogenous variables
            at the current time `t`. This might lead to leakage if for predictions the values of the exogenous
            variables at time `t` are not known.
        **kwargs
            Additional keyword arguments passed to `sklearn.linear_model.LinearRegression`.
        """
        self.kwargs = kwargs
        super().__init__(
            lags=lags,
            lags_exog=lags_exog,
            model=LinearRegression(**kwargs)
        )

    def __str__(self):
        return 'LinearRegression(lags={}, lags_exog={})'.format(self.lags, self.lags_exog)


text:
#
# TODO: Use a more efficient algorithm for DLA.
#

import numpy as np
import matplotlib.pyplot as plt
from update import UpdateAll
from random_walk import RandWalk

side = 300
num_walkers = side**2//8
time_steps = 1500

agg = UpdateAll(side, 8)        # aggregation nbr
rw = RandWalk(side)             # random walk

init = rw.create_init(num_walkers, False)  # all walkers labeled as 1
init[side//2, side//2] = 2                 # sticky cell

def stick(cell, nbrs):
    if cell == 2:             # sticky cell
        b = nbrs == 1         # non-sticky nbrs
        if np.any(b):
            nbrs[b] = 2
    return nbrs

def run_once(init):
    ca1 = agg.update_nbrs(init, stick)
    return rw.run(ca1, 1, 1)

def run(init, n=1):
    out = [init]
    for _ in range(n):
        out.append(run_once(out[-1]))
    return out

def run_sans_memory(init, n=1):
    out = init
    for _ in range(n):
        out = run_once(out)
    return out

res = run_sans_memory(init, time_steps)
plt.axis('off')
plt.imshow(res)

plt.show()