text:
''' Script to check the correctness of clustering for masked pixels.

    Masked pixels are pixels that are disabled or noisy.
'''

import unittest
import os

import numpy as np

from pixel_clusterizer.clusterizer import HitClusterizer, default_hits_dtype, default_clusters_dtype, default_cluster_hits_dtype

def create_hits(n_hits, max_column, max_row, max_frame, max_charge, hit_dtype=default_hits_dtype, hit_fields=None):
    hits = np.zeros(shape=(n_hits, ), dtype=hit_dtype)
    if not hit_fields:
        for i in range(n_hits):
            hits[i]['event_number'], hits[i]['frame'], hits[i]['column'], hits[i]['row'], hits[i]['charge'] = i / 3, i % max_frame, i % max_column + 1, 2 * i % max_row + 1, i % max_charge
    else:
        hit_fields_inverse = dict((v, k) for k, v in hit_fields.items())
        for i in range(n_hits):
            hits[i][hit_fields_inverse['event_number']], hits[i][hit_fields_inverse['frame']], hits[i][hit_fields_inverse['column']], hits[i][hit_fields_inverse['row']], hits[i][hit_fields_inverse['charge']] = i / 3, i % max_frame, i % max_column + 1, 2 * i % max_row + 1, i % max_charge
    return hits

class TestClusterizer(unittest.TestCase):

    @classmethod
    def setUpClass(cls):
        cls.pure_python = os.getenv('PURE_PYTHON', False)

    def test_disabled_pixels(self):
        # Create some fake data
        hits = np.ones(shape=(7, ), dtype=default_hits_dtype)
        hits[0]['column'], hits[0]['row'], hits[0]['charge'], hits[0]['event_number'], hits[0]['frame'] = 1, 2, 4, 0, 0
        hits[1]['column'], hits[1]['row'], hits[1]['charge'], hits[1]['event_number'], hits[1]['frame'] = 2, 2, 4, 0, 0
        hits[2]['column'], hits[2]['row'], hits[2]['charge'], hits[2]['event_number'], hits[2]['frame'] = 2, 2, 5, 1, 10
        hits[3]['column'], hits[3]['row'], hits[3]['charge'], hits[3]['event_number'], hits[3]['frame'] = 2, 2, 6, 2, 0
        hits[4]['column'], hits[4]['row'], hits[4]['charge'], hits[4]['event_number'], hits[4]['frame'] = 2, 3, 6, 2, 0
        hits[5]['column'], hits[5]['row'], hits[5]['charge'], hits[5]['event_number'], hits[5]['frame'] = 3, 3, 6, 2, 0
        hits[6]['column'], hits[6]['row'], hits[6]['charge'], hits[6]['event_number'], hits[6]['frame'] = 3, 3, 7, 3, 11

        # Create clusterizer object
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=2, row_cluster_distance=2, frame_cluster_distance=4, ignore_same_hits=True)

        # Case 1: Test max hit charge cut, accept all hits
        cluster_hits, clusters = clusterizer.cluster_hits(hits, disabled_pixels=[[2, 2], [3, 3]])  # cluster hits

        # Check cluster
        expected_clusters = np.zeros(shape=(2, ), dtype=default_clusters_dtype)
        expected_clusters['event_number'] = [0, 2]
        expected_clusters['ID'] = [0, 0]
        expected_clusters['n_hits'] = [1, 1]
        expected_clusters['charge'] = [4, 6]
        expected_clusters['seed_column'] = [1, 2]
        expected_clusters['seed_row'] = [2, 3]
        expected_clusters['mean_column'] = [1.0, 2.0]
        expected_clusters['mean_row'] = [2.0, 3.0]

        expected_cluster_hits = np.zeros(shape=(7, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [0, -1, -1, -1, 0, -1, -1]
        expected_cluster_hits['is_seed'] = [1, 0, 0, 0, 1, 0, 0]
        expected_cluster_hits['cluster_size'] = [1, 0, 0, 0, 1, 0, 0]
        expected_cluster_hits['n_cluster'] = [1, 1, 0, 1, 1, 1, 0]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

    def test_noisy_pixels(self):
        # Create some fake data
        hits = np.ones(shape=(9, ), dtype=default_hits_dtype)
        hits[0]['column'], hits[0]['row'], hits[0]['charge'], hits[0]['event_number'], hits[0]['frame'] = 1, 2, 8, 0, 0
        hits[1]['column'], hits[1]['row'], hits[1]['charge'], hits[1]['event_number'], hits[1]['frame'] = 2, 2, 4, 0, 0
        hits[2]['column'], hits[2]['row'], hits[2]['charge'], hits[2]['event_number'], hits[2]['frame'] = 2, 2, 5, 1, 10
        hits[3]['column'], hits[3]['row'], hits[3]['charge'], hits[3]['event_number'], hits[3]['frame'] = 2, 2, 12, 2, 0
        hits[4]['column'], hits[4]['row'], hits[4]['charge'], hits[4]['event_number'], hits[4]['frame'] = 2, 3, 6, 2, 0
        hits[5]['column'], hits[5]['row'], hits[5]['charge'], hits[5]['event_number'], hits[5]['frame'] = 3, 3, 3, 2, 0
        hits[6]['column'], hits[6]['row'], hits[6]['charge'], hits[6]['event_number'], hits[6]['frame'] = 3, 3, 7, 3, 11
        hits[7]['column'], hits[7]['row'], hits[7]['charge'], hits[7]['event_number'], hits[7]['frame'] = 3, 15, 1, 4, 0
        hits[8]['column'], hits[8]['row'], hits[8]['charge'], hits[8]['event_number'], hits[8]['frame'] = 20, 15, 1, 5, 0

        # Create clusterizer object
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=2, row_cluster_distance=2, frame_cluster_distance=4, ignore_same_hits=True)

        # Case 1: Test max hit charge cut, accept all hits
        cluster_hits, clusters = clusterizer.cluster_hits(hits, noisy_pixels=[[2, 2], [3, 3], [3, 15]])  # cluster hits

        # Check cluster
        expected_clusters = np.zeros(shape=(3, ), dtype=default_clusters_dtype)
        expected_clusters['event_number'] = [0, 2, 5]
        expected_clusters['ID'] = [0, 0, 0]
        expected_clusters['n_hits'] = [2, 3, 1]
        expected_clusters['charge'] = [8 + 4, 12 + 6 + 3, 1]
        expected_clusters['seed_column'] = [1, 2, 20]
        expected_clusters['seed_row'] = [2, 2, 15]
        expected_clusters['mean_column'] = [(9 * 1 + 5 * 2) / float(9 + 5), (13 * 2 + 7 * 2 + 4 * 3) / float(13 + 7 + 4), 20]
        expected_clusters['mean_row'] = [(9 * 2 + 5 * 2) / float(9 + 5), (13 * 2 + 7 * 3 + 4 * 3) / float(13 + 7 + 4), 15]

        expected_cluster_hits = np.zeros(shape=(9, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [0, 0, -1, 0, 0, 0, -1, -1, 0]
        expected_cluster_hits['is_seed'] = [1, 0, 0, 1, 0, 0, 0, 0, 1]
        expected_cluster_hits['cluster_size'] = [2, 2, 0, 3, 3, 3, 0, 0, 1]
        expected_cluster_hits['n_cluster'] = [1, 1, 0, 1, 1, 1, 0, 0, 1]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

    def test_noisy_and_disabled_pixels(self):
        # Create single hit data
        hits = np.zeros(shape=(1, ), dtype=default_hits_dtype)
        hits[0]['column'], hits[0]['row'] = 1, 1

        # Case 1: Test single noisy pixel
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=1, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, noisy_pixels=[[1, 1]])

        expected_clusters = np.zeros(shape=(0, ), dtype=default_clusters_dtype)

        expected_cluster_hits = np.zeros(shape=(1, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [-1]
        expected_cluster_hits['is_seed'] = [0]
        expected_cluster_hits['cluster_size'] = [0]
        expected_cluster_hits['n_cluster'] = [0]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Case 2: Test single disabled pixel
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=1, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, disabled_pixels=[[1, 1]])

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Create double hit data
        hits = np.zeros(shape=(2, ), dtype=default_hits_dtype)
        hits[0]['column'], hits[0]['row'] = 1, 1
        hits[1]['column'], hits[1]['row'] = 1, 2

        # Case 3: Test double noisy pixel
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=1, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, noisy_pixels=[[1, 1], [1, 2]])

        expected_clusters = np.zeros(shape=(1, ), dtype=default_clusters_dtype)
        expected_clusters['event_number'] = [0]
        expected_clusters['ID'] = [0]
        expected_clusters['n_hits'] = [2]
        expected_clusters['charge'] = [0]
        expected_clusters['seed_column'] = [1]
        expected_clusters['seed_row'] = [1]
        expected_clusters['mean_column'] = [1.0]
        expected_clusters['mean_row'] = [1.5]

        expected_cluster_hits = np.zeros(shape=(2, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [0, 0]
        expected_cluster_hits['is_seed'] = [1, 0]
        expected_cluster_hits['cluster_size'] = [2, 2]
        expected_cluster_hits['n_cluster'] = [1, 1]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Case 4: Test double noisy pixel
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=1, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, disabled_pixels=[[1, 1], [1, 2]])

        expected_clusters = np.zeros(shape=(0, ), dtype=default_clusters_dtype)

        expected_cluster_hits = np.zeros(shape=(2, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [-1, -1]
        expected_cluster_hits['is_seed'] = [0, 0]
        expected_cluster_hits['cluster_size'] = [0, 0]
        expected_cluster_hits['n_cluster'] = [0, 0]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Case 5: Test noisy and disabled pixel
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=1, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, noisy_pixels=[[1, 1]], disabled_pixels=[[1, 2]])

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Create triple hit data
        hits = np.zeros(shape=(3, ), dtype=default_hits_dtype)
        hits[0]['column'], hits[0]['row'] = 1, 1
        hits[1]['column'], hits[1]['row'] = 1, 2
        hits[2]['column'], hits[2]['row'] = 1, 3

        # Case 6: Test triple pixel
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=1, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, disabled_pixels=[[1, 1]], noisy_pixels=[[1, 2], [1, 3]])

        expected_clusters = np.zeros(shape=(1, ), dtype=default_clusters_dtype)
        expected_clusters['event_number'] = [0]
        expected_clusters['ID'] = [0]
        expected_clusters['n_hits'] = [2]
        expected_clusters['charge'] = [0]
        expected_clusters['seed_column'] = [1]
        expected_clusters['seed_row'] = [2]
        expected_clusters['mean_column'] = [1.0]
        expected_clusters['mean_row'] = [2.5]

        expected_cluster_hits = np.zeros(shape=(3, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [-1, 0, 0]
        expected_cluster_hits['is_seed'] = [0, 1, 0]
        expected_cluster_hits['cluster_size'] = [0, 2, 2]
        expected_cluster_hits['n_cluster'] = [1, 1, 1]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Create quadruple hit data
        hits = np.zeros(shape=(4, ), dtype=default_hits_dtype)
        hits[0]['column'], hits[0]['row'] = 1, 1
        hits[1]['column'], hits[1]['row'] = 1, 2
        hits[2]['column'], hits[2]['row'] = 1, 3
        hits[3]['column'], hits[3]['row'] = 1, 4

        # Case 7: Test quadruple pixel with single disabled pixel
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=1, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, disabled_pixels=[[1, 3]])

        expected_clusters = np.zeros(shape=(2, ), dtype=default_clusters_dtype)
        expected_clusters['event_number'] = [0, 0]
        expected_clusters['ID'] = [0, 1]
        expected_clusters['n_hits'] = [2, 1]
        expected_clusters['charge'] = [0, 0]
        expected_clusters['seed_column'] = [1, 1]
        expected_clusters['seed_row'] = [1, 4]
        expected_clusters['mean_column'] = [1.0, 1.0]
        expected_clusters['mean_row'] = [1.5, 4.0]

        expected_cluster_hits = np.zeros(shape=(4, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [0, 0, -1, 1]
        expected_cluster_hits['is_seed'] = [1, 0, 0, 1]
        expected_cluster_hits['cluster_size'] = [2, 2, 0, 1]
        expected_cluster_hits['n_cluster'] = [2, 2, 2, 2]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Case 8: Test quadruple pixel with single noisy pixel
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=1, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, noisy_pixels=[[1, 3]])

        expected_clusters = np.zeros(shape=(1, ), dtype=default_clusters_dtype)
        expected_clusters['event_number'] = [0]
        expected_clusters['ID'] = [0]
        expected_clusters['n_hits'] = [4]
        expected_clusters['charge'] = [0]
        expected_clusters['seed_column'] = [1]
        expected_clusters['seed_row'] = [1]
        expected_clusters['mean_column'] = [1.0]
        expected_clusters['mean_row'] = [2.5]

        expected_cluster_hits = np.zeros(shape=(4, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [0, 0, 0, 0]
        expected_cluster_hits['is_seed'] = [1, 0, 0, 0]
        expected_cluster_hits['cluster_size'] = [4, 4, 4, 4]
        expected_cluster_hits['n_cluster'] = [1, 1, 1, 1]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Create quintuple hit data
        hits = np.zeros(shape=(5, ), dtype=default_hits_dtype)
        hits[0]['column'], hits[0]['row'] = 1, 1
        hits[1]['column'], hits[1]['row'] = 1, 2
        hits[2]['column'], hits[2]['row'] = 1, 3
        hits[3]['column'], hits[3]['row'] = 1, 4
        hits[4]['column'], hits[4]['row'] = 1, 5

        # Case 9: Test quintuple pixel with 2 disabled pixels
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=1, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, disabled_pixels=[[1, 3], [1, 4]])

        expected_clusters = np.zeros(shape=(2, ), dtype=default_clusters_dtype)
        expected_clusters['event_number'] = [0, 0]
        expected_clusters['ID'] = [0, 1]
        expected_clusters['n_hits'] = [2, 1]
        expected_clusters['charge'] = [0, 0]
        expected_clusters['seed_column'] = [1, 1]
        expected_clusters['seed_row'] = [1, 5]
        expected_clusters['mean_column'] = [1.0, 1.0]
        expected_clusters['mean_row'] = [1.5, 5.0]

        expected_cluster_hits = np.zeros(shape=(5, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [0, 0, -1, -1, 1]
        expected_cluster_hits['is_seed'] = [1, 0, 0, 0, 1]
        expected_cluster_hits['cluster_size'] = [2, 2, 0, 0, 1]
        expected_cluster_hits['n_cluster'] = [2, 2, 2, 2, 2]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Case 10: Test quintuple pixel with 2 noisy pixel
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=1, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, noisy_pixels=[[1, 3]])

        expected_clusters = np.zeros(shape=(1, ), dtype=default_clusters_dtype)
        expected_clusters['event_number'] = [0]
        expected_clusters['ID'] = [0]
        expected_clusters['n_hits'] = [5]
        expected_clusters['charge'] = [0]
        expected_clusters['seed_column'] = [1]
        expected_clusters['seed_row'] = [1]
        expected_clusters['mean_column'] = [1.0]
        expected_clusters['mean_row'] = [3.0]

        expected_cluster_hits = np.zeros(shape=(5, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [0, 0, 0, 0, 0]
        expected_cluster_hits['is_seed'] = [1, 0, 0, 0, 0]
        expected_cluster_hits['cluster_size'] = [5, 5, 5, 5, 5]
        expected_cluster_hits['n_cluster'] = [1, 1, 1, 1, 1]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Case 11: Test quintuple pixel with single noisy and disabled pixels
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=1, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, noisy_pixels=[[1, 3]], disabled_pixels=[[1, 4]])

        expected_clusters = np.zeros(shape=(2, ), dtype=default_clusters_dtype)
        expected_clusters['event_number'] = [0, 0]
        expected_clusters['ID'] = [0, 1]
        expected_clusters['n_hits'] = [3, 1]
        expected_clusters['charge'] = [0, 0]
        expected_clusters['seed_column'] = [1, 1]
        expected_clusters['seed_row'] = [1, 5]
        expected_clusters['mean_column'] = [1.0, 1.0]
        expected_clusters['mean_row'] = [2.0, 5.0]

        expected_cluster_hits = np.zeros(shape=(5, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [0, 0, 0, -1, 1]
        expected_cluster_hits['is_seed'] = [1, 0, 0, 0, 1]
        expected_cluster_hits['cluster_size'] = [3, 3, 3, 0, 1]
        expected_cluster_hits['n_cluster'] = [2, 2, 2, 2, 2]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Case 12: Test quintuple pixel with single noisy and disabled pixels
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=2, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, noisy_pixels=[[1, 3]], disabled_pixels=[[1, 4]])

        expected_clusters = np.zeros(shape=(1, ), dtype=default_clusters_dtype)
        expected_clusters['event_number'] = [0]
        expected_clusters['ID'] = [0]
        expected_clusters['n_hits'] = [4]
        expected_clusters['charge'] = [0]
        expected_clusters['seed_column'] = [1]
        expected_clusters['seed_row'] = [1]
        expected_clusters['mean_column'] = [1.0]
        expected_clusters['mean_row'] = [2.75]

        expected_cluster_hits = np.zeros(shape=(5, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [0, 0, 0, -1, 0]
        expected_cluster_hits['is_seed'] = [1, 0, 0, 0, 0]
        expected_cluster_hits['cluster_size'] = [4, 4, 4, 0, 4]
        expected_cluster_hits['n_cluster'] = [1, 1, 1, 1, 1]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Create sextuple hit data
        hits = np.zeros(shape=(6, ), dtype=default_hits_dtype)
        hits[0]['column'], hits[0]['row'] = 1, 1
        hits[1]['column'], hits[1]['row'] = 1, 2
        hits[2]['column'], hits[2]['row'] = 1, 3
        hits[3]['column'], hits[3]['row'] = 1, 4
        hits[4]['column'], hits[4]['row'] = 1, 5
        hits[5]['column'], hits[5]['row'] = 1, 6

        # Case 13: Test sextuple pixel with noisy and disabled pixels
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=1, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, noisy_pixels=[[1, 3], [1, 5]], disabled_pixels=[[1, 4]])

        expected_clusters = np.zeros(shape=(2, ), dtype=default_clusters_dtype)
        expected_clusters['event_number'] = [0, 0]
        expected_clusters['ID'] = [0, 1]
        expected_clusters['n_hits'] = [3, 2]
        expected_clusters['charge'] = [0, 0]
        expected_clusters['seed_column'] = [1, 1]
        expected_clusters['seed_row'] = [1, 5]
        expected_clusters['mean_column'] = [1.0, 1.0]
        expected_clusters['mean_row'] = [2.0, 5.5]

        expected_cluster_hits = np.zeros(shape=(6, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [0, 0, 0, -1, 1, 1]
        expected_cluster_hits['is_seed'] = [1, 0, 0, 0, 1, 0]
        expected_cluster_hits['cluster_size'] = [3, 3, 3, 0, 2, 2]
        expected_cluster_hits['n_cluster'] = [2, 2, 2, 2, 2, 2]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Case 14: Test sextuple pixel with noisy and disabled pixels
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=2, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, noisy_pixels=[[1, 3], [1, 5]], disabled_pixels=[[1, 4]])

        expected_clusters = np.zeros(shape=(1, ), dtype=default_clusters_dtype)
        expected_clusters['event_number'] = [0]
        expected_clusters['ID'] = [0]
        expected_clusters['n_hits'] = [5]
        expected_clusters['charge'] = [0]
        expected_clusters['seed_column'] = [1]
        expected_clusters['seed_row'] = [1]
        expected_clusters['mean_column'] = [1.0]
        expected_clusters['mean_row'] = [3.4]

        expected_cluster_hits = np.zeros(shape=(6, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [0, 0, 0, -1, 0, 0]
        expected_cluster_hits['is_seed'] = [1, 0, 0, 0, 0, 0]
        expected_cluster_hits['cluster_size'] = [5, 5, 5, 0, 5, 5]
        expected_cluster_hits['n_cluster'] = [1, 1, 1, 1, 1, 1]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Case 15: Test sextuple pixel with noisy and disabled pixels
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=1, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, disabled_pixels=[[1, 3], [1, 5]], noisy_pixels=[[1, 4]])

        expected_clusters = np.zeros(shape=(2, ), dtype=default_clusters_dtype)
        expected_clusters['event_number'] = [0, 0]
        expected_clusters['ID'] = [0, 1]
        expected_clusters['n_hits'] = [2, 1]
        expected_clusters['charge'] = [0, 0]
        expected_clusters['seed_column'] = [1, 1]
        expected_clusters['seed_row'] = [1, 6]
        expected_clusters['mean_column'] = [1.0, 1.0]
        expected_clusters['mean_row'] = [1.5, 6.0]

        expected_cluster_hits = np.zeros(shape=(6, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [0, 0, -1, -1, -1, 1]
        expected_cluster_hits['is_seed'] = [1, 0, 0, 0, 0, 1]
        expected_cluster_hits['cluster_size'] = [2, 2, 0, 0, 0, 1]
        expected_cluster_hits['n_cluster'] = [2, 2, 2, 2, 2, 2]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

        # Case 16: Test sextuple pixel with noisy and disabled pixels
        clusterizer = HitClusterizer(pure_python=self.pure_python, min_hit_charge=0, max_hit_charge=13, charge_correction=1, charge_weighted_clustering=True, column_cluster_distance=1, row_cluster_distance=2, frame_cluster_distance=1, ignore_same_hits=True)
        cluster_hits, clusters = clusterizer.cluster_hits(hits, disabled_pixels=[[1, 3], [1, 5]], noisy_pixels=[[1, 4]])

        expected_clusters = np.zeros(shape=(1, ), dtype=default_clusters_dtype)
        expected_clusters['event_number'] = [0]
        expected_clusters['ID'] = [0]
        expected_clusters['n_hits'] = [4]
        expected_clusters['charge'] = [0]
        expected_clusters['seed_column'] = [1]
        expected_clusters['seed_row'] = [1]
        expected_clusters['mean_column'] = [1.0]
        expected_clusters['mean_row'] = [3.25]

        expected_cluster_hits = np.zeros(shape=(6, ), dtype=default_cluster_hits_dtype)
        expected_cluster_hits['event_number'] = hits['event_number']
        expected_cluster_hits['frame'] = hits['frame']
        expected_cluster_hits['column'] = hits['column']
        expected_cluster_hits['row'] = hits['row']
        expected_cluster_hits['charge'] = hits['charge']
        expected_cluster_hits['cluster_ID'] = [0, 0, -1, 0, -1, 0]
        expected_cluster_hits['is_seed'] = [1, 0, 0, 0, 0, 0]
        expected_cluster_hits['cluster_size'] = [4, 4, 0, 4, 0, 4]
        expected_cluster_hits['n_cluster'] = [1, 1, 1, 1, 1, 1]

        # Test results
        self.assertTrue(np.array_equal(clusters, expected_clusters))
        self.assertTrue(np.array_equal(cluster_hits, expected_cluster_hits))

if __name__ == '__main__':
    suite = unittest.TestLoader().loadTestsFromTestCase(TestClusterizer)
    unittest.TextTestRunner(verbosity=2).run(suite)

max_stars_repo_name: SiLab-Bonn/pyPixelClusterizer

max_issues_repo_name: SiLab-Bonn/pyPixelClusterizer


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

from levenshtein import editDistance, distance
from scipy.cluster.hierarchy import ward, dendrogram
import matplotlib.pyplot as plt
import sys
import codecs

def HierarchicalCluster(datafile):
    sim_list = []
    sim_matrix = []
    data_list = [line.rstrip('\n') for line in codecs.open(datafile, "r", "utf-8")]
    data_list = [line.split('   ') for line in data_list]
    print(data_list)
    len_list = len(data_list)
    label_list = [data_list[l][0] for l in range(0,len(data_list))]
    text_list = [data_list[l][1] for l in range(0,len(data_list))]
    print(label_list)
    print(text_list)

    for i in range(0, len_list):
        pivot = text_list[i]
        for j in range(0, len_list):
            sim = distance(pivot, text_list[j]) # calcurate similarity(distance)
            print('n{}, n{} : {}'.format(i, j, sim))
            sim_list.append(sim)
            if j == len_list-1:
                sim_matrix.append(sim_list)
                sim_list = []

    print('-------------------------')
    print('matrix: {}'.format(sim_matrix))
    linkage_matrix = ward(sim_matrix)
    print('--------------------------')
    print(linkage_matrix)
    dendrogram(linkage_matrix, labels=label_list)
    plt.show()

HierarchicalCluster(sys.argv[1])

ext: py

max_stars_repo_name: gentom/sentences-similarity-cluster

max_issues_repo_name: gentom/sentences-similarity-cluster

max_forks_repo_name: gentom/sentences-similarity-cluster

reason: from scipy


text:
import argparse
import numpy as np
import tensorflow as tf

from src.pcns_v2.datasets.cifar10 import get_dataset, input_fn
from src.pcns_v2.architectures.resnet import ResNet
from src.pcns_v2.adversarial.helpers import advTrainModel
from src.pcns_v2.helpers.optimization_helpers import piecewise_scheduler

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description='ResNet adversarial training')
    parser.add_argument('--adv_method', type=str, default='None', help="Method to generate adversarial examples")
    parser.add_argument('--scale', type=float, default=1.0, help="Width of the network")

    args = parser.parse_args()
    adv_method = args.adv_method
    scale = args.scale
    print("Method to generate adversarial samples: %s" % adv_method)
    print("Width factor: %f" % scale)

    TOTAL_EPOCHS = 164
    BATCH_SIZE = 128
    VERBOSE = 1  # 1 = progress bar, 2 = one line per epoch

    # Get train/test data sets
    (train_x, train_y), (val_x, val_y), (test_x, test_y) = get_dataset()
    num_train = train_x.shape[0]
    num_test = test_x.shape[0]

    train_dataset = input_fn(train_x, train_y, num_epochs=TOTAL_EPOCHS, batch_size=BATCH_SIZE,
                                 is_training=True, normalize=False)
    test_dataset = input_fn(test_x, test_y, batch_size=BATCH_SIZE, normalize=False)

    # Create ResNet model
    resnet = ResNet(input_shape=(32, 32, 3), num_classes=10, bottleneck=False, num_filters_at_start=int(16*scale),
                    initial_kernel_size=3, initial_conv_strides=1, initial_pool_size=None, initial_pool_strides=None,
                    num_residual_blocks_per_stage=[3, 3, 3], first_block_strides_per_stage=[1, 2, 2], kernel_size=3,
                    project_first_residual=True, version='V1', data_format='channels_last')
    model, _ = resnet.get_model()
    model = advTrainModel(model, attack=adv_method)

    lr_schedule = piecewise_scheduler([82, 123], [1.0, 0.1, 0.01], base_rate=0.1, boundaries_as='epochs',
                                      num_images=num_train, batch_size=BATCH_SIZE)
    optimizer = tf.keras.optimizers.SGD(learning_rate=lr_schedule, momentum=0.9)
    loss = tf.keras.losses.CategoricalCrossentropy(from_logits=True)

    model.compile(loss=loss, optimizer=optimizer, metrics=['accuracy'])
    print(model.base_model.summary())
    resnet.print_ResNet(model.base_model)

    # Train model
    history = model.fit(train_dataset, epochs=TOTAL_EPOCHS, steps_per_epoch=np.floor(num_train/BATCH_SIZE),
                        validation_data=test_dataset, validation_steps=np.ceil(num_test/BATCH_SIZE), validation_freq=1,
                        verbose=VERBOSE)

    # Best model in terms of validation accuracy, would have to checkpoint this model during training to evaluate
    # all three natural/FGSM/PGD accuracies (see experiment_scripts/resnet_on_cifar for model checkpointing)
    best_epoch = np.argmax(history.history['val_accuracy'])
    best_acc = history.history['val_accuracy'][best_epoch]
    print("Best accuracy: {:.5f} after epoch {}".format(best_acc, best_epoch+1))

    # Evaluate fully trained model
    print("Evaluation of fully trained model:")

    print('Natural Loss and Acc')
    model.adv_attack = 'None'
    model.compile(loss=loss, optimizer=optimizer, metrics=['accuracy'])
    score = model.evaluate(test_dataset, steps=np.ceil(num_test/BATCH_SIZE), verbose=0)
    print("Test Loss {:.5f}, Test Acc: {:.5f}".format(score[0], score[1]))

    print('FGSM Loss and Acc')
    model.adv_attack = 'FGSM'
    model.compile(loss=loss, optimizer=optimizer, metrics=['accuracy'])
    score = model.evaluate(test_dataset, steps=np.ceil(num_test/BATCH_SIZE), verbose=0)
    print("Test Loss {:.5f}, Test Acc: {:.5f}".format(score[0], score[1]))

    print('PGD Loss and Acc')
    model.adv_attack = 'PGD'
    model.compile(loss=loss, optimizer=optimizer, metrics=['accuracy'])
    score = model.evaluate(test_dataset, steps=np.ceil(num_test/BATCH_SIZE), verbose=0)
    print("Test Loss {:.5f}, Test Acc: {:.5f}".format(score[0], score[1]))

    # # Something like this might be helpful:
    # loss = tf.keras.losses.CategoricalCrossentropy(reduction=tf.keras.losses.Reduction.NONE)
    # test_set_losses = [] # hold the loss of every example from the test set
    # for idx, batch in enumerate(test_dataset):
    #     x, y = batch
    #     # attack x
    #     y_pred = model.base_model(x, training=False)
    #     batch_losses = loss(y, y_pred)
    #     test_set_losses.append(batch_losses)
    #
    # test_set_losses = tf.concat(test_set_losses, axis=0).numpy() # array of length 10,000k

reason: import numpy


text:
"""
The point of this snakemake script is to let you programmatically view and save screenshots of fMRI data.

Created 12/7/2021 by Benjamin Velie.
veliebm@gmail.com
"""
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)

from os import PathLike
from typing import List, Tuple
import nilearn.plotting
import nilearn.image
import nibabel
import numpy
import matplotlib.pyplot as plt

def main():
    make_plot(
        underlay_path=snakemake.input.underlay,
        overlay_path=snakemake.input.overlay,
        save_to_path=snakemake.output.plot,
        threshold=snakemake.params.threshold,
        coordinates=snakemake.params.coordinates,
        title=snakemake.params.title,
    )

def make_plot(underlay_path: PathLike, overlay_path: PathLike, save_to_path: PathLike, threshold: Tuple[float, float], coordinates: List[int], title: str) -> None:
    """
    Plot an underlay and overlay.
    """
    view = plt.figure(figsize=(10, 5))
    view = nilearn.plotting.plot_anat(underlay_path, cut_coords=coordinates, figure=view, draw_cross=False, annotate=False)
    overlay_image = load_as_nifti(overlay_path)

    if isinstance(threshold, tuple):
        thresholded_overlay = apply_threshold(overlay_image, threshold[0], threshold[1])
        view.add_overlay(thresholded_overlay, colorbar=True, cmap="cold_hot")
    else:
        view.add_overlay(overlay_image, colorbar=True, cmap="cold_hot", threshold=threshold)

    view.title(title)

    view.savefig(save_to_path)
    view.close()

def apply_threshold(image: nibabel.Nifti1Image, bottom_threshold: float, top_threshold: float) -> nibabel.Nifti1Image:
    """
    Apply a bottom threshold and top threshold to an image.
    """
    image_array = numpy.array(image.dataobj)
    image_array[(image_array>bottom_threshold) & (image_array<top_threshold)] = numpy.nan
    thresholded_image = nibabel.Nifti1Image(image_array, image.affine, image.header)
    return thresholded_image

def load_as_nifti(path_to_image: PathLike) -> nibabel.Nifti1Image:
    """
    Load an image and convert it into NIfTI format.

    For some reason, our AFNI headers don't play nice with nilearn's plotting functions.
    By transforming them into NIFTI's, we circumvent that problem.
    The affine, databj, header, and filemap are preserved.

    Args:
        path_to_image (PathLike): Path to an AFNI image.

    Returns:
        nibabel.Nifti1Image: A nibabel NIfTI image object containing the image data.
    """
    raw_image = nibabel.load(path_to_image)
    nifti_image = nibabel.Nifti1Image(raw_image.dataobj, raw_image.affine, header=raw_image.header, file_map=raw_image.file_map)

    return nifti_image

if __name__ == "__main__":
    main()

max_stars_repo_name: veliebm/contrascan

max_issues_repo_name: veliebm/contrascan

max_forks_repo_name: veliebm/contrascan


text:
from io import BytesIO
from PIL import Image
import numpy as np
import cv2
import math

def get_image_extension(django_image):
    '''
    Returns image extension from a django image
    '''
    pil_image = Image.open(django_image)
    return pil_image.format

def opencv_image_to_django_image(opencv_image, ext):
    opencv_image = cv2.cvtColor(opencv_image, cv2.COLOR_BGR2RGB)
    django_image = BytesIO()

    pil_image = Image.fromarray(opencv_image)
    pil_image.save(django_image, format=ext)

    return django_image

def django_image_to_opencv_image(django_image):
    pil_image = Image.open(django_image)
    img = np.array(pil_image)
    img = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)
    return img

def remove_duplicate_lines(sorted_points):
    '''
    Serches for the lines that are drawn
    over each other in the image and returns
    a list of non duplicate line co-ordinates
    '''
    last_x1 = 0
    non_duplicate_points = []
    for point in sorted_points:
        ((x1, y1), (x2, y2)) = point
        if last_x1 == 0:
            non_duplicate_points.append(point)
            last_x1 = x1

        elif abs(last_x1 - x1) >= 25:
            non_duplicate_points.append(point)
            last_x1 = x1

    return non_duplicate_points

def get_points_in_x_and_y(hough_lines, max_y):
    '''
    Takes a list of trigonometric form of lines
    and returns their starting and ending
    co-ordinates
    '''
    points = []
    for line in hough_lines:
        rho, theta = line[0]
        a = np.cos(theta)
        b = np.sin(theta)
        x0 = a * rho
        y0 = b * rho
        x1 = int(x0 + (max_y + 100) * (-b))
        y1 = int(y0 + (max_y + 100) * (a))
        start = (x1, y1)

        x2 = int(x0 - (max_y + 100) * (-b))
        y2 = int(y0 - (max_y + 100) * (a))
        end = (x2, y2)

        points.append((start, end))

    # Add a line at the very end of the image
    points.append(((500, max_y), (500, 0)))

    return points

def shorten_line(points, y_max):
    '''
    Takes a list of starting and ending
    co-ordinates of different lines
    and returns their trimmed form matching
    the image height
    '''
    shortened_points = []
    for point in points:
        ((x1, y1), (x2, y2)) = point

        # Slope
        try:
            m = (y2 - y1) / (x2 - x1)
        except ZeroDivisionError:
            m = -1  # Infinite slope

        if m == -1:
            shortened_points.append(((x1, y_max), (x1, 0)))
            continue

        # From equation of line:
        # y-y1 = m (x-x1)
        # x = (y-y1)/m + x1
        # let y = y_max
        new_x1 = math.ceil(((y_max - y1) / m) + x1)
        start_point = (abs(new_x1), y_max)

        # Now let y = 0
        new_x2 = math.ceil(((0 - y1) / m) + x1)
        end_point = (abs(new_x2), 0)

        shortened_points.append((start_point, end_point))

    return shortened_points

def get_cropped_images(image, points):
    '''
    Takes a spine line drawn image and
    returns a list of opencv images splitted
    from the drawn lines
    '''
    image = image.copy()
    y_max, _, _ = image.shape
    last_x1 = 0
    last_x2 = 0
    cropped_images = []

    for point in points:
        ((x1, y1), (x2, y2)) = point

        crop_points = np.array([[last_x1, y_max],
                                [last_x2, 0],
                                [x2, y2],
                                [x1, y1]])

        # Crop the bounding rect
        rect = cv2.boundingRect(crop_points)
        x, y, w, h = rect
        cropped = image[y: y + h, x: x + w].copy()

        # make mask
        crop_points = crop_points - crop_points.min(axis=0)
        mask = np.zeros(cropped.shape[:2], np.uint8)
        cv2.drawContours(mask, [crop_points], -1, (255, 255, 255), -1, cv2.LINE_AA)

        # do bit-op
        dst = cv2.bitwise_and(cropped, cropped, mask=mask)
        cropped_images.append(dst)

        last_x1 = x1
        last_x2 = x2

    return cropped_images

def resize_img(img):
    img = img.copy()
    img_ht, img_wd, _ = img.shape
    ratio = img_wd / img_ht
    new_width = 500
    new_height = math.ceil(new_width / ratio)
    resized_image = cv2.resize(img, (new_width, new_height))

    return resized_image

def detect_spines(img):
    '''
    Returns a list of lines seperating
    the detected spines in the image
    '''
    img = img.copy()
    height, width, _ = img.shape

    blur = cv2.GaussianBlur(img, (5, 5), 0)

    gray = cv2.cvtColor(blur, cv2.COLOR_BGR2GRAY)

    edge = cv2.Canny(gray, 50, 70)

    # kernel = np.ones((4, 1), np.uint8)
    kernel = np.array([[0, 0, 0, 0, 1, 0, 0, 0, 0],
                       [0, 0, 0, 0, 1, 0, 0, 0, 0],
                       [0, 0, 0, 0, 1, 0, 0, 0, 0],
                       [0, 0, 0, 0, 1, 0, 0, 0, 0],
                       [0, 0, 0, 0, 1, 0, 0, 0, 0]], dtype=np.uint8)

    img_erosion = cv2.erode(edge, kernel, iterations=1)

    lines = cv2.HoughLines(img_erosion, 1, np.pi / 180, 100)
    if lines is None:
        return []
    points = get_points_in_x_and_y(lines, height)
    points.sort(key=lambda val: val[0][0])
    non_duplicate_points = remove_duplicate_lines(points)

    final_points = shorten_line(non_duplicate_points, height)

    return final_points

def get_spines(django_image):
    img = django_image_to_opencv_image(django_image)
    ext = get_image_extension(django_image)

    final_image = resize_img(img)
    final_points = detect_spines(final_image)
    cropped_images = get_cropped_images(final_image, final_points)

    django_cropped_images = []
    for cropped_image in cropped_images:
        django_cropped_images.append(
            opencv_image_to_django_image(
                cropped_image,
                ext
            )
        )

    return django_cropped_images

def draw_spine_lines(django_image):
    img = django_image_to_opencv_image(django_image)
    ext = get_image_extension(django_image)

    final_image = resize_img(img)
    final_points = detect_spines(final_image)

    for point in final_points:
        ((x1, y1), (x2, y2)) = point
        final_image = cv2.line(final_image, (x1, y1), (x2, y2), (0, 0, 255), 10)

    django_image = opencv_image_to_django_image(
        final_image,
        ext
    )
    return django_image, ext

max_stars_repo_name: LakshyaKhatri/Bookshelf-Reader-API

max_issues_repo_name: imranansari/Bookshelf-Reader-API

max_forks_repo_name: imranansari/Bookshelf-Reader-API


text:
# Copyright 2018 Adrien CHEVRIER, Florian HEPP, Xavier HERMAND,
#                Gauthier LEONARD, Audrey LY, Elliot MAINCOURT
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     http://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import json
import csv
import numpy as np
import os
import hashlib

from itertools import chain, combinations
from collections import defaultdict

from sklearn.cluster import DBSCAN, AffinityPropagation
from sklearn.preprocessing import StandardScaler
from random import randrange
from datetime import datetime, timedelta
from dateutil import parser
from statistics import mean

from adapt.intent import IntentBuilder
from mycroft.skills.core import MycroftSkill
from mycroft.skills.core import intent_handler
from mycroft.skills.context import removes_context
from mycroft.messagebus.message import Message
from mycroft.util.log import getLogger, LOG

__authors__ = 'Nuttymoon, adrienchevrier, florianhepp'

# Logger: used for debug lines, like "LOGGER.debug(xyz)". These
# statements will show up in the command line when running Mycroft.
LOGGER = getLogger(__name__)

SKILLS_FOLDERS = {
    "/opt/mycroft/skills/PFE1718-skill-listener": "skill listener",
    "/opt/mycroft/skills/PFE1718-habit-miner": "habit miner",
    "/opt/mycroft/skills/PFE1718-automation-handler": "automation handler"
}

class HabitsManager(object):
    """
    This class manages the reading and writting in the file habits.json

    Attributes:
        habits_file_path (str): path to the file habits.json
        habits (json): the json datastore corresponding to habits.json
    """

    def __init__(self):
        self.habits_file_path = os.path.expanduser(
            "~/.mycroft/skills/ListenerSkill/habits/habits.json")
        self.triggers_file_path = os.path.expanduser(
            "~/.mycroft/skills/ListenerSkill/habits/triggers.json")

    def load_files(self):
        self.habits = json.load(open(self.habits_file_path))
        self.triggers = json.load(open(self.triggers_file_path))

    def get_all_habits(self):
        """Return all the existing habits of the user"""
        return self.habits

    def get_habit_by_id(self, habit_id):
        """Return one particular habit of the user"""
        return self.habits[habit_id]

    def check_skill_habit(self, new_habit):
        """
        Verify if skill habit is already present
        """

        old_habits = self.get_all_habits()
        new_utterances = [str(utt['utterance']) for utt in new_habit]
        # Check utterances for each old habit
        for ohabit in old_habits:
            # Extract utterances from old habit
            if ohabit['trigger_type'] in "skill":
                old_utterances = [
                    str(outt['last_utterance']) for outt in ohabit['intents']]
                # Check if new utterances are in old utterances
                if set(new_utterances).issubset(set(old_utterances)):
                    LOG.info('skill habit already exists')
                    return True
        LOG.info('skill habit does not exist')
        return False

    def fusion_habits(self, intent, old_intent):
        old_intent.append(intent)
        return old_intent

    def check_habit_presence(self, X, time, days, interval_max):
        """returns true if a habit with same
        trigger_type,time and days alreasy exists,
        returns False if not"""

        # Init variables
        intent = {
            "name": str(X[0, 3]),
            "parameters": X[0, 4],
            "last_utterance": str(X[0, 5])}

        # used to compare time with previous habits
        s_time = datetime.strptime(time, '%H:%M')
        # check habit presence
        for old_habit in self.habits:
            # If habit with same dates and time exists
            if old_habit['trigger_type'] == "time":
                if days in old_habit['days']:
                    # Verify if time between old habit and new habit
                    # are in interval_max
                    o_time = datetime.strptime(old_habit['time'], '%H:%M')
                    diff = abs(
                        s_time.hour * 60 + s_time.minute -
                        o_time.hour * 60 - o_time.minute)
                    # if habit with same name exists
                    if ((
                            str(intent['last_utterance']) in map(
                                lambda x: x[
                                    'last_utterance'], old_habit['intents']
                            )) and diff <= interval_max):
                        LOG.info("old habit found, no habit written")
                        return 0
                    # fusion if habits dont have same utterance
                    elif diff <= interval_max:
                        old_habit['intents'] = self.fusion_habits(
                            intent,
                            old_habit['intents'])

                        # Write fusionned habit
                        with open(
                                self.habits_file_path, 'w') as habits_file:
                            json.dump(self.habits, habits_file)
                        return 0
                """
                    else:
                        LOG.info('no habit found same day, new habit created')
                else:
                    LOG.info('no habit found, new habit created')
            else:
                LOG.info('no habit found, new habit created')
            """

        # register new habit
        self.register_habit("time", [intent], time, [days], str(interval_max))
        return 1

    def register_habit(self, trigger_type, intents, time=None, days=None,
                       interval_max=None):
        """
        Register a new habit in habits.json

        Args:
            trigger_type (str): the habit trigger type ("time" or "skill")
            intents (datastore): the intents that are part of the habit
            time (str): the time of the habit (if time based)
            days (int[]): the days of the habit (if time based)
        """

        if trigger_type == "skill":
            self.habits += [
                {
                    "intents": intents,
                    "trigger_type": trigger_type,
                    "automatized": 0,
                    "user_choice": False,
                    "triggers": []
                }
            ]
        else:
            self.habits += [
                {
                    "intents": intents,
                    "trigger_type": trigger_type,
                    "automatized": 0,
                    "user_choice": False,
                    "time": time,
                    "days": days,
                    "interval_max": interval_max
                }
            ]
        with open(self.habits_file_path, 'w') as habits_file:
            json.dump(self.habits, habits_file)

    def automate_habit(self, habit_id, auto, new_triggers=None):
        """
        Register the automation of a habit in the habits.json

        Args:
            habit_id (str): the id of the habit to automate
            triggers (str[]): the intents to register as triggers of the habit
            auto (int): 1 for full automation, 2 for habit offer when triggered
        """
        habit = self.habits[habit_id]
        habit["user_choice"] = True
        habit["automatized"] = auto

        if habit["trigger_type"] == "skill":
            if not self.triggers:
                for i in new_triggers:
                    self.triggers += [
                        {
                            "intent": habit["intents"][i]["name"],
                            "parameters": habit["intents"][i]["parameters"],
                            "habit_id": habit_id
                        }
                    ]
            else:
                if not self.check_triggers(habit_id, habit, new_triggers):
                    return False

            habit["triggers"] = new_triggers
            with open(self.triggers_file_path, 'w+') as triggers_file:
                json.dump(self.triggers, triggers_file)

        self.habits[habit_id] = habit
        with open(self.habits_file_path, 'w+') as habits_file:
            json.dump(self.habits, habits_file)

        return True

    def check_triggers(self, habit_id, habit, new_triggers):
        to_add = []
        for known_trig in self.triggers:
            for i in new_triggers:
                LOGGER.info("Testing trigger" + str(habit["intents"][int(i)]))
                if (
                        habit[
                            "intents"][i]["name"] == known_trig[
                                "intent"]) and (
                    habit[
                        "intents"][i]["parameters"] == known_trig[
                        "parameters"]):
                    return False

                to_add += [
                    {
                        "intent": habit["intents"][i]["name"],
                        "parameters": habit["intents"][i]["parameters"],
                        "habit_id": habit_id
                    }
                ]
        self.triggers += to_add

        return True

    def not_automate_habit(self, habit_id):
        """
        Register the user choice of not automatizing a habit

        Args:
            habit_id (str): the id of the habit to not automate
        """
        self.habits[habit_id]["user_choice"] = True
        self.habits[habit_id]["automatized"] = 0
        with open(self.habits_file_path, 'w+') as habits_file:
            json.dump(self.habits, habits_file)

    def get_trigger_by_id(self, trigger_id):
        """Return one particular habit trigger"""
        return self.triggers[trigger_id]

class HabitMinerSkill(MycroftSkill):
    """
    This class implements the habit miner skill
    """

    def __init__(self):
        super(HabitMinerSkill, self).__init__(
            name="HabitMinerSkill")
        self.logs_file_path = os.path.expanduser(
            "~/.mycroft/skills/ListenerSkill/habits/logs.json")
        self.to_install = []

    @intent_handler(IntentBuilder("LaunchMiningIntent")
                    .require("LaunchMiningKeyword"))
    def handle_launch_mining(self, message):
        if not self.check_skills_intallation():
            return

        process_mining(self.logs_file_path)

    # region Dependent skills installation

    def check_skills_intallation(self):
        LOGGER.info("Checking for skills install...")
        ret = True
        self.to_install = []

        for folder, skill in SKILLS_FOLDERS.iteritems():
            if not os.path.isdir(folder):
                ret = False
                self.to_install += [skill]

        if not ret:
            self.set_context("InstallMissingContextMiner")
            dial = ("To use the skill automation handler, you also have to "
                    "install the skill")
            num_skill = "this skill"
            skills_list = ""
            for skill in self.to_install[:-1]:
                skills_list += skill + ", "
            if len(self.to_install) > 1:
                num_skill = "these {} skills".format(len(self.to_install))
                skills_list += "and "
                dial += "s"
            skills_list += self.to_install[-1]
            self.speak(dial + " " + skills_list +
                       ". Should I install {} for you?".format(num_skill),
                       expect_response=True)
        return ret

    @intent_handler(IntentBuilder("InstallMissingIntent")
                    .require("YesKeyword")
                    .require("InstallMissingContextMiner").build())
    @removes_context("InstallMissingContextMiner")
    def handle_install_missing(self):
        LOGGER.info(self.to_install)
        for skill in self.to_install:
            LOGGER.info("Installing " + skill)
            self.emitter.emit(
                Message("recognizer_loop:utterance",
                        {"utterances": ["install " + skill],
                         "lang": 'en-us'}))

    @intent_handler(IntentBuilder("NotInstallMissingIntent")
                    .require("NoKeyword")
                    .require("InstallMissingContextMiner").build())
    @removes_context("InstallMissingContextMiner")
    def handle_not_install_missing(self):
        pass

# endregion

    def stop(self):
        pass

def create_skill():
    return HabitMinerSkill()

# Read json logs

def read_json(logs_file_path):
    data = []
    # Open Json
    for line in open(logs_file_path, 'r+'):
        data.append(json.loads(line))
    return data

# Parse json and extract variables variables

def parse_json(data):
    days = []
    times = []
    ids = []
    intents = []
    params = []
    utt = []
    # parse each json variable
    for i in range(len(data)):
        datetime_object = parser.parse(data[i]["datetime"])
        # Extract day of week
        my_day = date_to_days(datetime_object)
        days.append(my_day)
        # Extract hour of day
        my_hour = time_to_hours(datetime_object)
        times.append(my_hour)
        # Extract intents
        my_id = str(data[i]["type"]) + (str(data[i]["parameters"]))
        m_param = str(data[i]["parameters"])
        ids.append(my_id)
        my_intent = str(data[i]["type"])
        intents.append(my_intent)
        my_param = data[i]["parameters"]
        params.append(my_param)
        my_utterance = str(data[i]["utterance"])
        utt.append(my_utterance)
    X = np.array((days, times, ids, intents, params, utt))
    X = X.transpose()
    return X

def extract_ids(data):
    unique_ids = np.unique(data)
    return unique_ids

def random_date(start, end):
    """
    This function will return a random datetime between two datetime
    objects.
    """

    delta = end - start
    int_delta = (delta.days * 24 * 60 * 60) + delta.seconds
    random_second = randrange(int_delta)
    return start + timedelta(seconds=random_second)

def date_to_days(date):
    """
    This function returns day of week in range 1-10
    """
    day_t = (float(date.weekday()) * 10.0) / 6.0
    return day_t

def time_to_seconds(date):
    """
    Function returns time of the day in seconds"""
    td = date.time()
    return ((td.hour * 3600 + td.second) * 10 ** 6 + td.microsecond) / 10 ** 6

def time_to_hours(date):
    """
    Function returns time of the day in hours in range 1-10"""
    td = date.time()
    td_t = (((float(td.hour) * 60 + float(td.minute)) / 60) * 10.0) / 24.0
    return float(td_t)

# Write new habit
def write_habit(X, labels, core_samples_mask_dbscan):
    n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)
    # Black removed and is used for noise instead.
    unique_labels = set(labels)
    num_clusters = 0

    for k in unique_labels:
        zero = ""
        # Check if cluster
        if k >= 0:

            class_member_mask = (labels == k)

            xy = X[class_member_mask & core_samples_mask_dbscan]
            if (len(xy[:, 1]) > 1):
                day = int(float(round(mean(xy[:, 0].astype(float)), 0)))
                hour = xy[:, 1].astype(float).astype(int)
                minute = (xy[:, 1].astype(float) - hour)
                minute = map(lambda x: x * 60, minute)
                # calculate mean hour and minute of the cluster
                hour_moy = int(float(mean(hour)))
                min_moy = int(float(mean(minute)))
                # Round time to 5 minutes
                min_moy_rounded = float(float(min_moy) / 60.0)
                min_moy_rounded = int(round(min_moy_rounded * 12) * 5)
                # Calculate inteval max to detect habit
                interval_max = np.ceil(
                    max(
                        np.absolute(
                            (hour_moy * 60 + min_moy) - hour * 60 - minute)))
                if min_moy_rounded < 10:
                    zero += "0"
                time = str(hour_moy) + ":" + zero + str(min_moy_rounded)
                # Register ID, params, intents, days, hours
                my_habit_manager = HabitsManager()
                my_habit_manager.load_files()
                # Add new clusters to cluster counter
                num_clusters = my_habit_manager.check_habit_presence(
                    X, time, day, interval_max) + num_clusters

    return num_clusters

def process_mining(logs_file_path):
    """
    MAIN STEPS
    take logs as input, process mining and set habits as output
    """

    # Try to read json logs
    try:
        raw_data = read_json(logs_file_path)
        X = parse_json(raw_data)

    except:
        LOGGER.info("Unexpected error, exiting habit miner")
        return -1

    # compute dbscan per id
    unique_ids = extract_ids(X[:, 2])
    nb_clusters = 0

    # Compute dbscan for each application ID
    for i in unique_ids:

        # Filter points
        X_mini = X[X[:, 2] == i]

        # Check if points >2
        if (X_mini[:, 0].shape[0] > 2):
            X_mini_cluster = X_mini[:, [0, 1]]

            # dbscan
            core_samples_mask, labels = compute_DBSCAN(X_mini_cluster)

            # Convert hours to 24 range
            def to_24(X): return map(lambda x: (x / 10) * 24, X)

            # Convert days to 7 range
            def to_7(X): return map(lambda x: (x / 10) * 6, X)

            X_mini[:, 0] = to_7(X_mini[:, 0].astype(float))
            X_mini[:, 1] = to_24(X_mini[:, 1].astype(float))

            # Write clean and write new habit
            nb_clusters = write_habit(
                X_mini, labels, core_samples_mask) + nb_clusters

    # Print infos
    LOG.info("nb nb_clusters:%d", nb_clusters)
    LOG.info("processing finished")

    run_apriori(logs_file_path)
    LOG.info("apriori finished")

# MODELS
# Compute Affinity Propagation
def compute_AP(features):
    af = AffinityPropagation(preference=-50, damping=0.9).fit(features)
    cluster_centers_indices = af.cluster_centers_indices_
    labels = af.labels_
    return cluster_centers_indices, labels

def compute_DBSCAN(features):
    features = StandardScaler().fit_transform(features.astype(float))
    # #############################################################################
    # Compute DBSCAN
    db = DBSCAN(eps=0.25, min_samples=4).fit(features)
    core_samples_mask = np.zeros_like(db.labels_, dtype=bool)
    core_samples_mask[db.core_sample_indices_] = True
    labels = db.labels_
    return core_samples_mask, labels

def subsets(arr):
    """ Returns non empty subsets of arr"""
    return chain(*[combinations(arr, i + 1) for i, a in enumerate(arr)])

def return_item_with_min_support(
        item_set, transaction_list, min_support, freq_set):
    """calculates the support for items in the itemSet and returns a subset
   of the itemSet each of whose elements satisfies the minimum support"""
    _itemSet = set()
    local_set = defaultdict(int)

    for item in item_set:
        for transaction in transaction_list:
            if item.issubset(transaction):
                freq_set[item] += 1
                local_set[item] += 1

    for item, count in local_set.items():
        support = float(count) / len(transaction_list)

        if support >= min_support:
            _itemSet.add(item)

    return _itemSet

def join_set(item_set, length):
    """Join a set with itself and returns the n-element itemsets"""
    return set(
        [
            i.union(
                j) for i in item_set for j in item_set if len(
                    i.union(j)) == length])

def get_item_set_transaction_list(data_iterator):
    transaction_list = list()
    item_set = set()
    for record in data_iterator:
        transaction = frozenset(record)
        transaction_list.append(transaction)
        for item in transaction:
            item_set.add(frozenset([item]))  # Generate 1-itemSets
    return item_set, transaction_list

def apriori(data_iter, min_support, min_confidence):
    """
    run the apriori algorithm. data_iter is a record iterator
    Return both:
     - items (tuple, support)
     - rules ((pretuple, posttuple), confidence)
    """
    item_set, transaction_list = get_item_set_transaction_list(data_iter)

    freq_set = defaultdict(int)
    large_set = dict()
    # Global dictionary which stores (key=n-itemSets,value=support)
    # which satisfy minSupport

    assocRules = dict()
    # Dictionary which stores Association Rules

    one_c_set = return_item_with_min_support(item_set,
                                             transaction_list,
                                             min_support,
                                             freq_set)

    current_l_set = one_c_set
    k = 2
    while current_l_set != set([]):
        large_set[k - 1] = current_l_set
        current_l_set = join_set(current_l_set, k)
        current_c_set = return_item_with_min_support(current_l_set,
                                                     transaction_list,
                                                     min_support,
                                                     freq_set)
        current_l_set = current_c_set
        k = k + 1

    def get_support(item):
        """local function which Returns the support of an item"""
        return float(freq_set[item]) / len(transaction_list)

    to_ret_items = []
    for key, value in large_set.items():
        to_ret_items.extend([(tuple(item), get_support(item))
                             for item in value])

    to_ret_rules = []
    for key, value in large_set.items()[1:]:
        for item in value:
            _subsets = map(frozenset, [x for x in subsets(item)])
            for element in _subsets:
                remain = item.difference(element)
                if len(remain) > 0:
                    confidence = get_support(item) / get_support(element)
                    if confidence >= min_confidence:
                        to_ret_rules.append(
                            ((tuple(element), tuple(remain)), confidence))
    return to_ret_items, to_ret_rules

def data_from_file(fname):
    """Function which reads from the file and yields a generator"""
    file_iter = open(fname, 'rU')
    for line in file_iter:
        line = line.strip().rstrip(',')  # Remove trailing comma
        record = frozenset(line.split(','))
        yield record

def run_apriori(logs_file_path, min_supp=0.05, min_confidence=0.8):
    hashes_temp = []
    table_csv = []
    csv_path = os.path.expanduser(
        "~/.mycroft/skills/HabitMinerSkill/inputApriori.csv")
    date_time_obj0 = datetime.strptime(
        '2018-01-01 00:00:00.0', '%Y-%m-%d %H:%M:%S.%f')
    habit_manager = HabitsManager()
    habit_manager.load_files()

    if not os.path.getsize(logs_file_path):
        return

    """
    Open logs and put them in a list,
    same line if consequent logs are within 5 minutes interval
    """
    with open(logs_file_path) as json_data:
        for i, line in enumerate(json_data):
            data = json.loads(line)
            date_time_obj1 = datetime.strptime(
                data['datetime'], '%Y-%m-%d %H:%M:%S.%f')
            delta = date_time_obj1 - date_time_obj0
            del data['datetime']
            hash = hashlib.md5(json.dumps(data)).hexdigest()

            if delta > timedelta(minutes=5):
                table_csv.append(hashes_temp)
                hashes_temp = []
                hashes_temp.append(hash)
                date_time_obj0 = date_time_obj1
            else:
                hashes_temp.append(hash)

    del table_csv[0]

    """
    Converts logs list to csv so that
    we can executre apriori algorithm on them
    """
    with open(csv_path, 'w+') as fp:
        LOG.info('opened')
        writer = csv.writer(fp, delimiter=',')
        for row in table_csv:
            writer.writerow(row)

    in_file = data_from_file(csv_path)

    _, rules = apriori(in_file, min_supp, min_confidence)

    # reformat the rules and sort the tuples in it
    hashes_temp = []

    for rule in rules:
        hash = rule[0][0] + rule[0][1]
        hashes_temp.append(sorted(hash))

    # this is to remove duplicates
    hashes = []

    for tuple in hashes_temp:
        if tuple not in hashes:
            hashes.append(tuple)

    # Hash to json data
    habits = []
    habit = []
    intents = []

    for hash in hashes:
        for intent in hash:
            with open(logs_file_path) as json_data:
                for line in json_data:
                    data = json.loads(line)
                    del data['datetime']
                    hash = hashlib.md5(json.dumps(data)).hexdigest()
                    if hash == intent:
                        habit.append(json.dumps(data))
                        intents.append(data)
                        break
        if not habit_manager.check_skill_habit(intents):
            habits.append(habit)
        habit = []

    # format habits as expected and register them
    intents = []

    for habit in habits:
        for intent in habit:
            intent = {
                'last_utterance': json.loads(intent)['utterance'],
                'name': json.loads(intent)['type'],
                'parameters': json.loads(intent)['parameters']
            }
            intents.append(intent)
        habit_manager.register_habit("skill", intents)
        intents = []

    os.remove(csv_path)


text:
from __future__ import absolute_import, division, print_function

import numpy as np
import tensorflow as tf
import random

# Enable font colors
class bcolors:
  """ For the purpose of print in terminal with colors """
  HEADER = '\033[95m'
  OKBLUE = '\033[94m'
  OKGREEN = '\033[92m'
  WARNING = '\033[93m'
  FAIL = '\033[91m'
  ENDC = '\033[0m'
  BOLD = '\033[1m'
  UNDERLINE = '\033[4m'

def obs2state(observation, low):
  """
  Helper function convert observation to discrete state
  We only use x and y to represent state, so the first 2 element in observation
  Args:
    observation: 0-d numpy array e.g. (-1.2, 3.3)
    low: lower bound of observation
  Return:
    state: a scalar
  """
  x = observation[0]
  y = observation[1]
  ind_x = int(x - low[0]) # index of x
  ind_y = int(y - low[1])
  state = ind_x*10 + ind_y

  return state

def generate_action_sequence(num_sequences, len_horizon, num_actions):
  """ Generate S random action sequences with H horizon
  """
  action_sequences = np.zeros((num_sequences, len_horizon))
  for s in range(num_sequences):
    for h in range(len_horizon):
      action_sequences[s,h] = random.randrange(num_actions)

  return action_sequences

def sample_to_batch(samples_list, num_states, num_actions):
  """ Create training batch from sampled memory
  """
  x_batch = np.zeros((len(samples_list), num_states+num_actions))
  y_batch = np.zeros((len(samples_list), num_states))
  for i, s in enumerate(samples_list):
    onehot_action = np.zeros(num_actions)
    onehot_action[s[1]] = 1
    x_batch[i] = np.concatenate((s[0], onehot_action))
    y_batch[i] = s[-1]

  return x_batch, y_batch

def create_dataset(input_features, output_labels, batch_size, shuffle=True, num_epochs=None):
  """ Create TF dataset from numpy arrays
  """
  dataset = tf.data.Dataset.from_tensor_slices((input_features, output_labels))
  if shuffle:
    dataset = dataset.shuffle(buffer_size = 1000)
  dataset = dataset.batch(batch_size)
  dataset = dataset.repeat(num_epochs)

  return dataset

def shoot_action(model, action_sequences, state, goal):
  """ Find an action with most reward using random shoot
  """
  sequence_rewards = np.zeros(action_sequences.shape[0])
  # Compute reward for every sequence
  for seq in range(action_sequences.shape[0]):
    old_state = np.array(state).reshape(1,-1).astype(np.float32)
    # print("old_state: {}".format(old_state)) # debug
    reward_in_horizon = 0
    for hor in range(action_sequences.shape[1]):
      action = np.array([[action_sequences[seq,hor]]])
      stac = np.concatenate((old_state, action), axis=1).astype(np.float32) # state-action pair
      new_state = model(stac)
      # print("new_state: {}".format(new_state)) # debug
      if np.linalg.norm(new_state[0,:2]-goal) < np.linalg.norm(old_state[0,:2]-goal):
        reward = 1
      else:
        reward = 0
      reward_in_horizon += reward
      old_state = new_state
      sequence_rewards[seq] = reward_in_horizon

    idx = np.argmax(sequence_rewards) # action sequence index with max reward
    optimal_action = int(action_sequences[idx,0]) # take first action of each sequence

    return optimal_action

reason: import numpy,from numpy


text:
'''简单的RNN model,预测效果一般,所以后面并没有使用RNN'''
import os
import time
import ast
import argparse
import copy
import gc
import numpy as np
import pandas as pd
import datetime as dt
from tqdm import tqdm

import torch
import torch.nn as nn

from torch.optim.lr_scheduler import ReduceLROnPlateau
from torch.utils.data import DataLoader
from torchvision import transforms

from data_loader import DoodleDataset
from utils import *

from visdom import Visdom

viz = Visdom(env = "rnn_model")
assert viz.check_connection()

best_model_path = "./rnn_model/debug/csv71-90-epoch-03-step-008000-map3-0.799841.pth"  # 加载最好的模型继续训练或者进行测试

csv_files_path = '../../input/shuffled_csv'
model_dir = "./rnn_model"

NCATS=340
NCSVS=100
CHANNEL=3

DEVICE = torch.device("cuda" if torch.cuda.is_available() else "cpu")
# if DEVICE == "cuda":
#     torch.cuda.set_device(device)

def get_HParams():
    parser = argparse.ArgumentParser(description='Get hyper params of the model.')
    parser.add_argument("-tag", type=str, default="debug")
    parser.add_argument("-max_seq_len", type=int, default=500)
    parser.add_argument("-hidden_size", type=int, default=512)
    parser.add_argument("-rnn_model", type=str, default='lstm', help="The encoder RNN model, lstm or gru")
    parser.add_argument("-num_layers", type=int, default=1)
    parser.add_argument("-batch_size", type=int, default=128)
    parser.add_argument("-learning_rate", type=float, default=0.001)
    parser.add_argument("-dropout_prob", type=float, default=0.0)  # 根据torch.nn.LSTM, 当num_layers=1时, 这个参数无效

#   parser.add_argument("-image_size", type=int, default=224)
    parser.add_argument("-epochs", type=int, default=20)
    parser.add_argument("-gpus", type=int, required=True)

    hp = parser.parse_args()
    return hp

class RNN_Model(nn.Module):
    def __init__(self, hp, batch_first=True):
        super(RNN_Model, self).__init__()
        self.hp = hp
        self.batch_first = batch_first
        # 下面的全连接层是用encoder的输出h求decoder的初始h0和C0, 乘2是因为是双向的
        self.fc_output = nn.Linear(2 * hp.hidden_size, NCATS)
        nn.init.normal_(self.fc_output.weight, std=0.001)
        self.relu = nn.ReLU(inplace=True)

        # bidirectional RNN:
        if hp.rnn_model == 'lstm':
            self.rnn = nn.LSTM(input_size=5, hidden_size=hp.hidden_size, num_layers=hp.num_layers,
                               dropout=hp.dropout_prob, bidirectional=True, batch_first=batch_first)
        elif hp.rnn_model == 'gru':
            self.rnn = nn.GRU(input_size=5, hidden_size=hp.hidden_size, num_layers=hp.num_layers,
                              dropout=hp.dropout_prob, bidirectional=True, batch_first=batch_first)
        else:
            print("ERROR:RNN model must be lstm or gru!!!")
            assert False

        # 激活 dropout:
        self.train()

    def feature(self, inputs, act_lens, h=None, c=None):
        """
        :param inputs shape (batch_first): (batch_size, seq_len, 5) or (not batch_first) (seq_len, batch_size, 5)
        :param h:
        :param c:
        :return: output, shape: (batch_size, seq_len, 2*enc_size)
        """
#         if h is None:
#             # then must init with zeros
#             h = torch.zeros(self.hp.batch_size, 2 * self.hp.num_layers,  self.hp.hidden_size, device=DEVICE)
#             c = torch.zeros(self.hp.batch_size, 2 * self.hp.num_layers, self.hp.hidden_size, device=DEVICE)

        if self.hp.rnn_model == 'lstm':
            # output包含了最后一层lstm的所有时间步的hidden_state, 而h和c分别包含了所有层的最后一个时间步的hidden_state和cell_state
            output, (h, c) = self.rnn(inputs.float())
        elif self.hp.rnn_model == 'gru':
            # output包含了最后一层gru的所有时间步的hidden_state, 而h包含了所有层的最后一个时间步的hidden_state
            output, h = self.rnn(inputs.float())
        else:
            assert False

        # output shape: (seq_len, batch_size, 2*hidden_size)
        # output = output.permute(1, 0, 2)  # 现在output shape: (batch_size, seq_len, 2*hidden_size)
        act_lens = act_lens - 1  # index 是从0开始的
        indices = act_lens.view(-1, 1, 1).repeat(1, 1, output.shape[-1]) # shape: (batch, 1, 2*hidden_size)
        act_last_h = torch.gather(output, 1, indices) # 取实际的最后一个h, shape: (batch, 1, 2*hidden_size)
        act_last_h = torch.squeeze(act_last_h, 1)  # shape: (batch, 2*hidden_size)

        return act_last_h

    def forward(self, inputs, act_lens, h=None, c=None):
        act_last_h = self.feature(inputs, act_lens, h, c)
        x = self.fc_output(act_last_h)
        logit = self.relu(x)

        # logit shape: (batch_size, 340)
        return logit

def plot_by_visdom(win, x, y, opts):
    if win is None:
        win = viz.line(X = x, Y = y, opts=opts)
    else:
        viz.line(X= x, Y = y, win=win, update='append', opts=opts)

    return win

def valid_model(model, valid_data_loader):
    model.eval()   # Set model to   mode

    correct_num = 0
    map3s = []
    top3_accs = []
    for samples_batch in tqdm(valid_data_loader):
        inputs = samples_batch["strokes"]
        labels = samples_batch["label"]
        act_lens = samples_batch["point_num"]

        inputs = inputs.to(DEVICE, dtype=torch.float)
        act_lens = act_lens.to(DEVICE, dtype=torch.long)
        labels = labels.to(DEVICE, dtype=torch.long)

        with torch.set_grad_enabled(False):
            # forward
            outputs = model(inputs, act_lens)
            top3_acc, map3 = mapk(outputs, labels, k=3)
            map3s.append(map3)
            top3_accs.append(top3_acc)

            _, preds = torch.max(outputs, 1)
            correct_num += torch.sum(preds.int() == labels.int())

    valid_acc = correct_num.double().item() / len(valid_data_loader.dataset)
    valid_map3 = sum(map3s) / float(len(map3s))
    valid_top3_acc = sum(top3_accs) / float(len(top3_accs))

    model.train()  # Set model to training mode

    return valid_acc, valid_map3, valid_top3_acc

def train_model(hps, model, train_data_loader, valid_data_loader, criterion, optimizer, num_epochs=20,
                print_every_step=20, save_every_step=None, batch_accumulate_size=1):
    train_start = time.time()
    train_steps = int(np.ceil(len(train_data_loader.dataset) / float(hps.batch_size)))
    if save_every_step is None:
        save_every_step = int(train_steps / 10)  # 1/10个epoch检查(保存)一次

    model = model.to(DEVICE) # to gpu

    scheduler = ReduceLROnPlateau(optimizer, mode='max', factor=0.5, patience=2, verbose=True)

    # best_model_wts = copy.deepcopy(model.state_dict())
    # best_valid_map3 = 0.0
    best_valid_map3 = float("0." + best_model_path.split(".")[-2])

    valid_win = None
    valid_opts = dict(title="valid", xlabel="step", ylabel="acc or map3",
                          legend=['valid_acc', "valid_map3", "valid_top3_acc"])

    step_count = 0
    for epoch in range(num_epochs):
        print('Epoch %d/%d' % (epoch+1, num_epochs))
        print('-' * 20)
        model.train()  # Set model to training mode

        # 每个epoch一个图
        train_win = None
        train_opts = dict(title="train-epoch-%d"%(epoch+1), xlabel="step", ylabel="acc or loss",
                          legend=['train_acc', "train_loss"])

        # Iterate over data.
        for step, samples_batch in enumerate(train_data_loader):
            start = time.time()
            inputs = samples_batch["strokes"]
            act_lens = samples_batch["point_num"]
            labels = samples_batch["label"]

            # loss, preds, step_acc = train_step(model, inputs, labels, criterion, optimizer)
            # 上面一行代码就是普通的训练过程,下面这部分使用了梯度累加的trick,相当于增大了batch_size
            # 参考https://discuss.pytorch.org/t/how-to-implement-accumulated-gradient/3822中Gopal_Sharma的系列回答
            #################将batch_accumulate_size个batch的梯度积累起来,只在最后一次更新网络参数###################
            inputs = inputs.to(DEVICE, dtype=torch.float)
            act_lens = act_lens.to(DEVICE, dtype=torch.long)
            labels = labels.to(DEVICE, dtype=torch.float)
            if step % batch_accumulate_size == 0:
                optimizer.zero_grad()
            with torch.set_grad_enabled(True):
                # forward
                outputs = model(inputs, act_lens)
                loss = criterion(outputs, labels.long()) / batch_accumulate_size # 一定要除以这个size,原因见上面链接的讨论
                loss.backward()

                _, preds = torch.max(outputs, 1)
            correct_num = torch.sum(preds == labels.long())
            step_acc = correct_num.double() / inputs.size(0)

            if (step + 1) % batch_accumulate_size == 0:
                optimizer.step()

            loss = batch_accumulate_size * loss.item() # 转换为数字方便后面用visdom画图
            step_acc = step_acc.item()
            ########################################################################################

            train_win = plot_by_visdom(train_win, np.array([step_count]),
                                       np.column_stack((np.array([step_acc]), np.array([loss]))), train_opts)

            if step % print_every_step == 0:
                print('Train RNN epoch %d/%d step %d/%d. Loss: %.4f Acc: %.4f.  time taken: %.2fs' % (
                    epoch+1, num_epochs, step, train_steps, loss, step_acc, time.time()-start))

            if (step > 0 and step % save_every_step == 0) or (step == train_steps - 1):
                valid_acc, valid_map3, valid_top3_acc = valid_model(model, valid_data_loader)
                scheduler.step(valid_map3) # 适当调整学习率
                print('Validating model ...\nbest_valid_map3:%.6f valid_acc:%.4f valid_map3:%.6f valid_top3_acc:%.6f' % (
                        best_valid_map3, valid_acc, valid_map3, valid_top3_acc))

                valid_win = plot_by_visdom(valid_win, np.array([step_count]),
                                           np.column_stack((np.array([valid_acc]),
                                                            np.array([valid_map3]),
                                                            np.array([valid_top3_acc]))),
                                           valid_opts)

                if best_valid_map3 < valid_map3:
                    best_valid_map3 = valid_map3
                    print("Better model. best_valid_map3 --> %.6f \nsaving model..." % best_valid_map3, end='')
                    save_name = "csv71-90-epoch-%02d-step-%06d-map3-%.6f.pth" % (epoch + 1, step, best_valid_map3)
                    save_model(model, model_save_dir=model_dir + "/" + hps.tag, save_name=save_name)

            del inputs, labels, outputs
            gc.collect()
            step_count += 1

    # valid_acc, valid_map3, valid_top3_acc = valid_model(model, valid_data_loader)
    save_name = "whole_model_epoch_%d.pth" % num_epochs
    save_model(model, model_save_dir=model_dir + "/" + hps.tag, save_name=save_name)

    time_elapsed = time.time() - train_start
    print('Training complete in {:.0f}m {:.0f}s'.format(time_elapsed // 60, time_elapsed % 60))
    print('Best valid map3:%.6f' % best_valid_map3)

def main():
    hps = get_HParams()

    hps_dict = vars(hps)
    for key, val in hps_dict.items():
        print('%s = %s' % (key, str(val)))

    # train_csv_files, valid_csv_files = get_csv_files_randomly(hps.train_files_num, hps.valid_files_num)
    train_ks = list(range(71, 91))
    if 12 in train_ks:
        train_ks.remove(12)  # 12是验证集
    assert 12 not in train_ks
    train_csv_files = [os.path.join(csv_files_path, 'train_%d_%d.csv.gz'%(k, NCSVS)) for k in train_ks]
    valid_csv_files = [os.path.join(csv_files_path, 'train_12_100.csv.gz')]

    model = RNN_Model(hp=hps)
    model.load_state_dict(torch.load(best_model_path)) # 将之前训练的模型load进来

    if hps.gpus > 1:
        print("------------use %d GPUs!------------" % hps.gpus)
        model = nn.DataParallel(model)
    # 如果使用了DataParallel,那么load也应该在它后面load,否则会报key对不上的错, 参考https://www.ptorch.com/news/74.html
    # model.load_state_dict(torch.load(best_model_path)) # 将之前训练的模型load进来

    train_data_loader = create_data_loader(csv_files=train_csv_files, channel=CHANNEL, hps=hps,
                                           transform=None, batch_num=None, for_cnn=False, for_rnn=True)
    valid_data_loader = create_data_loader(csv_files=valid_csv_files, channel=CHANNEL, hps=hps,
                                           transform=None, batch_num=400, for_cnn=False, for_rnn=True)

    optimizer = torch.optim.Adam(model.parameters(), lr=hps.learning_rate, amsgrad=True)
    criterion = nn.CrossEntropyLoss()

    train_model(hps, model, train_data_loader, valid_data_loader, criterion, optimizer,
                              num_epochs=hps.epochs, print_every_step=400, save_every_step=4000, batch_accumulate_size=8)

if __name__ == '__main__':
    main()

max_stars_repo_name: ShusenTang/kaggle-doodle-recognition

max_issues_repo_name: ShusenTang/kaggle-doodle-recognition

max_forks_repo_name: ShusenTang/kaggle-doodle-recognition


text:
import unittest
import numpy as np

from Core.Complex import Complex
from Core.Structure import StructureAgent
from TS.Edge import Edge
from TS.State import State, Vector, Memory
from TS.TransitionSystem import TransitionSystem

class TestTransitionSystem(unittest.TestCase):
    def setUp(self):
        self.str1 = StructureAgent("X", set())
        self.str2 = StructureAgent("Y", set())
        self.str3 = StructureAgent("Z", set())
        self.str4 = StructureAgent("W", set())

        self.c1 = Complex([self.str1], "rep")
        self.c2 = Complex([self.str2], "rep")
        self.c3 = Complex([self.str3], "rep")
        self.c4 = Complex([self.str4], "rep")

        ordering = (self.c1, self.c2, self.c3)
        ordering_wrong = (self.c1, self.c2, self.c3, self.c4),
        ordering_reordered = (self.c3, self.c1, self.c2)

        self.s1 = State(Vector(np.array((1, 2, 3))), Memory(0))
        self.s2 = State(Vector(np.array((1, 2, 4))), Memory(0))
        self.s3 = State(Vector(np.array((1, 2, 5))), Memory(0))
        self.s4 = State(Vector(np.array((4, 2, 3))), Memory(0))
        self.hell = State(Vector(np.array((np.inf, np.inf, np.inf))), Memory(0), True)

        self.s1_reordered = State(Vector(np.array((3, 1, 2))), Memory(0))
        self.s2_reordered = State(Vector(np.array((4, 1, 2))), Memory(0))
        self.s3_reordered = State(Vector(np.array((5, 1, 2))), Memory(0))

        self.ts = TransitionSystem(tuple(), 5)
        self.ts.states_encoding = {1: self.s1, 2: self.s2}

        self.edge_1 = Edge(0, 1, 0.5)

        # equality of TS

        self.ts_eq_1 = TransitionSystem(ordering, 5)
        self.ts_eq_1.states_encoding = {0: self.s1, 1: self.s2, 2: self.s3}
        self.ts_eq_1.edges = {Edge(0, 1, 0.3), Edge(0, 2, 0.7), Edge(1, 2, 1)}

        self.ts_eq_2 = TransitionSystem(ordering, 5)
        self.ts_eq_2.states_encoding = {1: self.s1, 2: self.s2, 0: self.s3}
        self.ts_eq_2.edges = {Edge(1, 2, 0.3), Edge(1, 0, 0.7), Edge(2, 0, 1)}

        self.ts_eq_3 = TransitionSystem(ordering_reordered, 5)
        self.ts_eq_3.states_encoding = {1: self.s1_reordered, 2: self.s2_reordered, 0: self.s3_reordered}
        self.ts_eq_3.edges = {Edge(1, 2, 0.3), Edge(1, 0, 0.7), Edge(2, 0, 1)}

        # inequality of TS

        self.ts_neq_1 = TransitionSystem(ordering, 5)
        self.ts_neq_1.states_encoding = {0: self.s1, 1: self.s2, 2: self.s3}
        self.ts_neq_1.edges = {Edge(0, 1, 0.3), Edge(0, 2, 0.7), Edge(1, 2, 1)}

        self.ts_neq_2 = TransitionSystem(ordering, 5)
        self.ts_neq_2.states_encoding = {1: self.s1, 2: self.s2, 0: self.s3}
        self.ts_neq_2.edges = {Edge(1, 2, 0.3), Edge(1, 0, 0.7), Edge(0, 2, 1)}

        self.ts_neq_3 = TransitionSystem(ordering_wrong, 5)
        self.ts_neq_3.states_encoding = {1: self.s1, 2: self.s2, 0: self.s3}
        self.ts_neq_3.edges = {Edge(1, 2, 0.3), Edge(1, 0, 0.7), Edge(0, 2, 1)}

        # test bigger TS

        ordering = (self.c1, self.c2, self.c3, self.c4)

        self.ts_bigger = TransitionSystem(ordering, 5)
        self.ts_bigger.states_encoding = {1: self.s1, 2: self.s2, 0: self.s3, 3: self.s4}
        self.ts_bigger.edges = {Edge(1, 2, 0.8), Edge(0, 1, 0.5), Edge(1, 0, 0.3),
                                Edge(3, 1, 0.9), Edge(3, 2, 0.1), Edge(1, 3, 0.2)}

    def test_eq(self):
        self.assertEqual(self.ts_eq_1, self.ts_eq_2)
        self.assertEqual(self.ts_eq_1, self.ts_eq_3)
        self.assertNotEqual(self.ts_neq_1, self.ts_neq_2)
        self.assertNotEqual(self.ts_neq_1, self.ts_neq_3)

    def test_get_state_encoding(self):
        self.assertEqual(self.ts.states_encoding[2], self.s2)
        self.assertEqual(self.ts.states_encoding[1], self.s1)
        self.assertEqual(len(self.ts.states_encoding), 2)

    def test_add_edge(self):
        self.assertEqual(Edge(0, 1, 0.5), self.edge_1)

    def test_hash_edge(self):
        self.assertNotEqual(hash(Edge(1, 2, 0.3)), hash(Edge(2, 1, 0.3)))

    def test_edge_comparision(self):
        self.assertTrue(Edge(1, 2, 0.3) < Edge(2, 2, 0.3))
        self.assertTrue(Edge(1, 2, 0.3) < Edge(1, 3, 0.3))
        self.assertFalse(Edge(1, 2, 0.3) > Edge(2, 2, 0.3))
        self.assertTrue(Edge(4, 0, 0.3) < Edge(5, 2, 0.3))

    def test_sort_edges(self):
        edges = {Edge(1, 2, 0.3), Edge(2, 2, 0.3), Edge(4, 0, 0.3), Edge(4, 1, 0.3), Edge(5, 2, 0.3), Edge(1, 3, 0.3)}
        edges_sorted = [Edge(1, 2, 0.3), Edge(1, 3, 0.3), Edge(2, 2, 0.3),
                        Edge(4, 0, 0.3), Edge(4, 1, 0.3), Edge(5, 2, 0.3)]
        self.assertEqual(sorted(edges), edges_sorted)

    def test_change_hell(self):
        ordering = (self.c1, self.c2, self.c3, self.c4)
        ts = TransitionSystem(ordering, 4)
        ts.states_encoding = {1: self.s1, 2: self.s2, 0: self.s3, 3: self.hell}
        ts.change_hell()
        new_hell = State(Vector(np.array([5, 5, 5])), Memory(0), True)
        new_encoding = {1: self.s1, 2: self.s2, 0: self.s3, 3: new_hell}
        self.assertEqual(ts.states_encoding, new_encoding)

max_stars_repo_name: sybila/eBCSgen

max_issues_repo_name: sybila/eBCSgen

max_forks_repo_name: sybila/eBCSgen


text:
import time, os
import numpy as np
import json

class Timer:
    def __init__(self, name, remove_start_msg=True):
        self.name = name
        self.remove_start_msg = remove_start_msg

    def __enter__(self):
        self.start_time = time.time()
        print('Run "%s".........' % self.name, end='\r' if self.remove_start_msg else '\n')

    def __exit__(self, exc_type, exc_val, exc_tb):
        time_diff = float(time.time() - self.start_time)
        time_str = '{:.1f}s'.format(time_diff) if time_diff >= 1 else '{:.0f}ms'.format(time_diff * 1000)

        print('Finish "{}" in {}'.format(self.name, time_str))

def output_csv(the_path, data_dict, order=None, delimiter=','):
    if the_path.endswith('.tsv'):
        delimiter = '\t'

    is_file_exists = os.path.exists(the_path)
    with open(the_path, 'a+') as op:
        keys = list(data_dict.keys())
        if order is not None:
            keys = order + [k for k in keys if k not in order]

        col_title = delimiter.join([str(k) for k in keys])
        if not is_file_exists:
            print(col_title, file=op)
        else:
            old_col_title = open(the_path, 'r').readline().strip()
            if col_title != old_col_title:
                old_order = old_col_title.split(delimiter)

                additional_keys = [k for k in keys if k not in old_order]
                if len(additional_keys) > 0:
                    print('WARNING! The data_dict has following additional keys %s' % (str(additional_keys)))

                no_key = [k for k in old_order if k not in keys]
                if len(no_key) > 0:
                    raise(RuntimeError('The data_dict does not have the following old keys: %s' % str(no_key)))

                keys = old_order + additional_keys

        print(delimiter.join([str(data_dict[k]) for k in keys]), file=op)

def vector_in(vec, names):
    is_kept = (vec == names[0])
    for m_name in names[1:]:
        is_kept = (is_kept | (vec == m_name))
    return is_kept