text:
from unittest import TestCase, main as unittest_main
import os
import pandas as pd
import numpy as np

from q2_winnowing.step1_3.pipeline import main as step1_3_main

class Step1_3Tests( TestCase ):
    # <><><> Testing class for Step 1 to 3 <><><>

    package = 'q2_winnowing.tests'

    # <><> read input values and expected output values for testing <><>
    testing_data = {
        "betweenness": (
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_in_data.csv",
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_out_features_bw.csv",
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_out_important_bw.csv",
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_out_abundances_bw.csv",
        ),
        "degree": (
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_in_data.csv",
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_out_features_dg.csv",
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_out_important_dg.csv",
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_out_abundances_dg.csv",
        ),
        "closeness": (
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_in_data.csv",
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_out_features_cl.csv",
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_out_important_cl.csv",
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_out_abundances_cl.csv",
        ),
        "eigenvector": (
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_in_data.csv",
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_out_features_ei.csv",
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_out_important_ei.csv",
            f"{os.path.dirname(os.path.realpath(__file__))}/sample_data/step1_3/test_out_abundances_ei.csv",
        )
    }

    def test_step1_3_main_graph_bw(self):

        data_path, features_path, important_path, adundances_path = self.testing_data["betweenness"]
        data_in = pd.read_csv( data_path )
        data_in.name = "TESTING_FRAME"
        features_out = pd.read_csv( features_path )
        important_out = pd.read_csv( important_path, index_col=0 )
        abundances_out = pd.read_csv( adundances_path )

        features_result, important_result, abundances_result = \
            step1_3_main( False, data_in, None, metric_name="graph_centrality", c_type="add_one", min_count=3,
                          total_select=25, iteration_select=128, pca_components=4, smooth_type="sliding_window",
                          window_size=3, centrality_type="betweenness", keep_threshold=0.5, correlation="spearman",
                          weighted=True, corr_prop="both", evaluation_type="kl_divergence", min_connected=0,
                          detailed=False )

        run_out = features_out.pop("run time")
        run_result = features_result.pop("run time")
        self.assertLess( abs( run_result.tail(1).iloc[0] - run_out.tail(1).iloc[0] ), 100 )
            # test that time for step is reasonably close to expected ( within range of 100 sec )

        # Remove names column
        features_out.pop("dataframe1")
        features_result.pop("dataframe1")

        np.testing.assert_array_equal(
            features_result.to_numpy().astype(object),
            features_out.to_numpy().astype(object),
        )
        np.testing.assert_array_almost_equal(
            important_result.to_numpy().astype(object),
            important_out.to_numpy().astype(object),
            15
        )
        np.testing.assert_array_equal(
            abundances_result.to_numpy().astype(object),
            abundances_out.to_numpy().astype(object),
        )

    def test_step1_3_main_graph_dg(self):

        data_path, features_path, important_path, adundances_path = self.testing_data["degree"]
        data_in = pd.read_csv(data_path)
        data_in.name = "TESTING_FRAME"
        features_out = pd.read_csv(features_path)
        important_out = pd.read_csv(important_path, index_col=0)
        abundances_out = pd.read_csv(adundances_path)

        features_result, important_result, abundances_result = \
            step1_3_main(False, data_in, None, metric_name="graph_centrality", c_type="add_one", min_count=3,
                         total_select=25, iteration_select=128, pca_components=4, smooth_type="sliding_window",
                         window_size=3, centrality_type="degree", keep_threshold=0.5,
                         correlation="spearman",
                         weighted=True, corr_prop="both", evaluation_type="kl_divergence", min_connected=0,
                         detailed=False )

        run_out = features_out.pop("run time")
        run_result = features_result.pop("run time")
        self.assertLess(abs(run_result.tail(1).iloc[0] - run_out.tail(1).iloc[0]), 100)
        # test that time for step is reasonably close to expected ( within range of 100 sec )

        # Remove names column
        features_out.pop("dataframe1")
        features_result.pop("dataframe1")

        np.testing.assert_array_equal(
            features_result.to_numpy().astype(object),
            features_out.to_numpy().astype(object),
        )
        np.testing.assert_array_almost_equal(
            important_result.to_numpy().astype(object),
            important_out.to_numpy().astype(object),
            15
        )
        np.testing.assert_array_equal(
            abundances_result.to_numpy().astype(object),
            abundances_out.to_numpy().astype(object),
        )

    def test_step1_3_main_graph_ei(self):

        data_path, features_path, important_path, adundances_path = self.testing_data["eigenvector"]
        data_in = pd.read_csv(data_path)
        data_in.name = "TESTING_FRAME"
        features_out = pd.read_csv(features_path)
        important_out = pd.read_csv(important_path, index_col=0)
        abundances_out = pd.read_csv(adundances_path)

        features_result, important_result, abundances_result = \
            step1_3_main(False, data_in, None, metric_name="graph_centrality", c_type="add_one", min_count=3,
                         total_select=25, iteration_select=128, pca_components=4, smooth_type="sliding_window",
                         window_size=3, centrality_type="eigenvector", keep_threshold=0.5,
                         correlation="spearman",
                         weighted=True, corr_prop="both", evaluation_type="kl_divergence", min_connected=0,
                         detailed=False )

        run_out = features_out.pop("run time")
        run_result = features_result.pop("run time")
        self.assertLess(abs(run_result.tail(1).iloc[0] - run_out.tail(1).iloc[0]), 100)
        # test that time for step is reasonably close to expected ( within range of 100 sec )

        # Remove names column
        features_out.pop("dataframe1")
        features_result.pop("dataframe1")

        np.testing.assert_array_equal(
            features_result.to_numpy().astype(object),
            features_out.to_numpy().astype(object),
        )
        np.testing.assert_array_almost_equal(
            important_result.to_numpy().astype(object),
            important_out.to_numpy().astype(object),
            15
        )
        np.testing.assert_array_equal(
            abundances_result.to_numpy().astype(object),
            abundances_out.to_numpy().astype(object),
        )

    def test_step1_3_main_graph_cl(self):

        data_path, features_path, important_path, adundances_path = self.testing_data["closeness"]
        data_in = pd.read_csv(data_path)
        data_in.name = "TESTING_FRAME"
        features_out = pd.read_csv(features_path)
        important_out = pd.read_csv(important_path, index_col=0)
        abundances_out = pd.read_csv(adundances_path)

        features_result, important_result, abundances_result = \
            step1_3_main(False, data_in, None, metric_name="graph_centrality", c_type="add_one", min_count=3,
                         total_select=25, iteration_select=128, pca_components=4, smooth_type="sliding_window",
                         window_size=3, centrality_type="closeness", keep_threshold=0.5,
                         correlation="spearman",
                         weighted=True, corr_prop="both", evaluation_type="kl_divergence", min_connected=0,
                         detailed=False )

        run_out = features_out.pop("run time")
        run_result = features_result.pop("run time")
        self.assertLess(abs(run_result.tail(1).iloc[0] - run_out.tail(1).iloc[0]), 100)
        # test that time for step is reasonably close to expected ( within range of 100 sec )

        # Remove names column
        features_out.pop("dataframe1")
        features_result.pop("dataframe1")

        np.testing.assert_array_equal(
            features_result.to_numpy().astype(object),
            features_out.to_numpy().astype(object),
        )
        np.testing.assert_array_almost_equal(
            important_result.to_numpy().astype(object),
            important_out.to_numpy().astype(object),
            15
        )
        np.testing.assert_array_equal(
            abundances_result.to_numpy().astype(object),
            abundances_out.to_numpy().astype(object),
        )

if __name__ == '__main__':
    unittest_main()


text:
from amuse.community import *
from amuse.test.amusetest import TestWithMPI

from amuse.community.bonsai2.interface import Bonsai2Interface, Bonsai2

import numpy

from amuse.units import nbody_system
from amuse.units import units
from amuse import datamodel
from amuse.rfi import channel
from amuse.ic.plummer import new_plummer_model
from amuse.support.exceptions import AmuseException

default_options = dict()
#default_options = dict(redirection="none")
#default_options = dict(redirection="none", debugger="gdb")

class TestBonsaiInterface(TestWithMPI):

    def test1(self):
        plummer_size = 500
        plummer =  new_plummer_model(plummer_size)
        mass=plummer.mass.number
        radius=plummer.radius.number
        x=plummer.x.number
        y=plummer.y.number
        z=plummer.z.number
        vx=plummer.vx.number
        vy=plummer.vy.number
        vz=plummer.vz.number

        instance = self.new_instance_of_an_optional_code(Bonsai2Interface, **default_options)
        self.assertEqual(0, instance.initialize_code())
        self.assertEqual([0, 0], list(instance.get_number_of_particles().values()))
        ids, errors = instance.new_particle(mass,x,y,z,vx,vy,vz,radius)
        self.assertEqual(0, errors)
        self.assertEqual(list(range(plummer_size)), ids)
        self.assertEqual(0, instance.commit_particles())

        self.assertEqual([500, 0], list(instance.get_number_of_particles().values()))
        masses, errors = instance.get_mass(list(range(500)))
        self.assertEqual(0, errors)
        self.assertAlmostEqual(0.002, masses)
        masses,xs,ys,zs,vxs,vys,vzs,radii, errors = instance.get_state(list(range(500)))
        self.assertEqual(0, errors)
        self.assertAlmostRelativeEquals(xs, x, 6)

        self.assertEqual(0, instance.evolve_model(0.00001))
        energy_total_init = instance.get_potential_energy()["potential_energy"] + instance.get_kinetic_energy()["kinetic_energy"]
        self.assertEqual(0, instance.evolve_model(1))
        energy_total_final = instance.get_potential_energy()["potential_energy"] + instance.get_kinetic_energy()["kinetic_energy"]
        self.assertAlmostRelativeEqual(energy_total_init, energy_total_final, 2)
        instance.stop()

class TestBonsai(TestWithMPI):

    def test1(self):
        print("Testing Bonsai initialization")
        instance = self.new_instance_of_an_optional_code(Bonsai2, **default_options)
        instance.initialize_code()
        print("cleaning the code")
        instance.cleanup_code()
        print("done")
        instance.stop()

    def test2(self):
        print("Testing Bonsai parameters")
        instance = self.new_instance_of_an_optional_code(Bonsai2,  **default_options)
        instance.initialize_code()
        self.assertAlmostEqual(instance.parameters.epsilon_squared, 0.0025 | nbody_system.length**2)
        self.assertAlmostEqual(instance.parameters.timestep, 1.0 / 64 | nbody_system.time)
        instance.parameters.epsilon_squared = 0.01 | nbody_system.length**2
        instance.parameters.timestep = 0.001 | nbody_system.time
        instance.commit_parameters()
        print((instance.parameters.epsilon_squared))
        print((instance.parameters.timestep))
        self.assertAlmostEqual(instance.parameters.epsilon_squared, 0.01 | nbody_system.length**2)
        self.assertAlmostEqual(instance.parameters.timestep, 0.001 | nbody_system.time)
        instance.stop()

    def test3(self):
        print("Testing Bonsai, N-body units")
        instance = self.new_instance_of_an_optional_code(Bonsai2, **default_options)
        instance.initialize_code()
        plummer = new_plummer_model(500)
        instance.particles.add_particles(plummer)
        self.assertAlmostEqual(instance.particles.mass, 0.002 | nbody_system.mass)
        instance.evolve_model(1.0 | nbody_system.time)
        self.assertAlmostEqual(instance.model_time, 1.0 | nbody_system.time)
        instance.stop()

    def test4(self):
        print("Testing Bonsai, SI units")
        convert_nbody = nbody_system.nbody_to_si(100.0 | units.MSun, 1.0 | units.parsec)
        instance = self.new_instance_of_an_optional_code(Bonsai2, convert_nbody, **default_options)
        instance.initialize_code()
        plummer = new_plummer_model(500, convert_nbody = convert_nbody)
        instance.particles.add_particles(plummer)
        instance.commit_particles()
        instance.evolve_model(1 | nbody_system.time)
        instance.stop()

    def test5(self):
        print("Testing Bonsai remove_particle")
        convert_nbody = nbody_system.nbody_to_si(100.0 | units.MSun, 1.0 | units.parsec)
        instance = self.new_instance_of_an_optional_code(Bonsai2, convert_nbody, **default_options)
        instance.initialize_code()
        plummer = new_plummer_model(500, convert_nbody = convert_nbody)
        instance.particles.add_particles(plummer)
        instance.commit_particles()
        instance.particles.remove_particle(plummer[0])
#    self.assertRaises(AmuseException, instance.particles.remove_particle, plummer[0],
#            expected_message = "Error when calling 'delete_particle' of a 'Bonsai', errorcode is -2")
        print(instance.particles[0].mass)
        instance.stop()

    def test6(self):
        print("Testing Bonsai states")
        plummer = new_plummer_model(500)

        print("First do everything manually:")
        instance = self.new_instance_of_an_optional_code(Bonsai2, **default_options)
        self.assertEqual(instance.get_name_of_current_state(), 'UNINITIALIZED')
        instance.initialize_code()
        self.assertEqual(instance.get_name_of_current_state(), 'INITIALIZED')
        instance.commit_parameters()
        self.assertEqual(instance.get_name_of_current_state(), 'EDIT')
        instance.particles.add_particles(plummer)
        instance.commit_particles()
        self.assertEqual(instance.get_name_of_current_state(), 'RUN')
        mass = instance.particles[0].mass
        instance.evolve_model(0.1 | nbody_system.time)
        self.assertEqual(instance.get_name_of_current_state(), 'EVOLVED')
        instance.cleanup_code()
        self.assertEqual(instance.get_name_of_current_state(), 'END')
        instance.stop()

        print(("initialize_code(), commit_parameters(), commit_particles(), and cleanup_code() should be called " +
            "automatically before editing parameters, new_particle(), get_xx(), and stop():"))
        instance = self.new_instance_of_an_optional_code(Bonsai2, **default_options)
        self.assertEqual(instance.get_name_of_current_state(), 'UNINITIALIZED')
        timestep = instance.parameters.timestep
        self.assertEqual(instance.get_name_of_current_state(), 'INITIALIZED')
        instance.particles.add_particles(plummer)
        self.assertEqual(instance.get_name_of_current_state(), 'EDIT')
        mass = instance.particles[0].mass
        self.assertEqual(instance.get_name_of_current_state(), 'RUN')
        instance.evolve_model(0.1 | nbody_system.time)
        self.assertEqual(instance.get_name_of_current_state(), 'EVOLVED')
        instance.stop()
        self.assertEqual(instance.get_name_of_current_state(), 'STOPPED')

    def test7(self):
        print("Testing Bonsai properties")
        numpy.random.seed(12345)
        plummer = new_plummer_model(500)
        instance = self.new_instance_of_an_optional_code(Bonsai2, **default_options)
        instance.particles.add_particles(plummer)

        self.assertEqual(instance.model_time, 0.0 | nbody_system.time)
        instance.evolve_model(0.0001 | nbody_system.time)
        print(instance.model_time)
        self.assertTrue(instance.model_time >= 0.0001 | nbody_system.time)
        self.assertEqual(instance.model_time, 1.0 * instance.parameters.timestep)

        self.assertAlmostEqual(instance.potential_energy, -0.50625962019 | nbody_system.energy)
        self.assertAlmostEqual(instance.kinetic_energy,    0.244611829519 | nbody_system.energy)
        self.assertAlmostEqual(instance.total_mass, 1.0 | nbody_system.mass)

        E0 = instance.kinetic_energy + instance.potential_energy

        self.assertRaises(AmuseException, getattr, instance, "total_radius", expected_message =
            "Error when calling 'get_total_radius' of a 'Bonsai2', errorcode is -2, "
            "error is 'Called function is not implemented.'")
        self.assertRaises(AmuseException, getattr, instance, "center_of_mass_position", expected_message =
            "Error when calling 'get_center_of_mass_position' of a 'Bonsai2', "
            "errorcode is -2, error is 'Called function is not implemented.'")
        self.assertRaises(AmuseException, getattr, instance, "center_of_mass_velocity", expected_message =
            "Error when calling 'get_center_of_mass_velocity' of a 'Bonsai2', "
            "errorcode is -2, error is 'Called function is not implemented.'")

        instance.evolve_model(1.0 | nbody_system.time)
        self.assertAlmostEqual(instance.model_time, 1.0 | nbody_system.time)
        self.assertAlmostEqual(instance.potential_energy, -0.5115 | nbody_system.energy, 4)
        self.assertAlmostEqual(instance.kinetic_energy,    0.24986 | nbody_system.energy, 4)
        self.assertAlmostEqual(
            instance.kinetic_energy + instance.potential_energy,
            E0, 4)
        instance.particles.remove_particle(plummer[2])
        instance.evolve_model(2.0 | nbody_system.time)
        instance.particles.remove_particle(plummer[20])
        instance.particles.remove_particle(plummer[30])
        instance.evolve_model(3.0 | nbody_system.time)
        instance.particles.remove_particle(plummer[35])
        instance.stop()

"""
    def test8(self):
        print "Testing Bonsai2 collision_detection"
        particles = datamodel.Particles(7)
        particles.mass = 0.001 | nbody_system.mass
        particles.radius = 0.01 | nbody_system.length
        particles.x = [-101.0, -100.0, -0.5, 0.5, 100.0, 101.0, 104.0] | nbody_system.length
        particles.y = 0 | nbody_system.length
        particles.z = 0 | nbody_system.length
        particles.velocity = [[2, 0, 0], [-2, 0, 0]]*3 + [[-4, 0, 0]] | nbody_system.speed

        instance = self.new_instance_of_an_optional_code(Bonsai2)
        instance.initialize_code()
        instance.parameters.set_defaults()
        instance.particles.add_particles(particles)
        collisions = instance.stopping_conditions.collision_detection
        collisions.enable()
        instance.evolve_model(1.0 | nbody_system.time)

        self.assertTrue(collisions.is_set())
        self.assertTrue(instance.model_time < 0.5 | nbody_system.time)
        self.assertEquals(len(collisions.particles(0)), 3)
        self.assertEquals(len(collisions.particles(1)), 3)
        self.assertEquals(len(particles - collisions.particles(0) - collisions.particles(1)), 1)
        self.assertEquals(abs(collisions.particles(0).x - collisions.particles(1).x) <
                (collisions.particles(0).radius + collisions.particles(1).radius),
                [True, True, True])

        sticky_merged = datamodel.Particles(len(collisions.particles(0)))
        sticky_merged.mass = collisions.particles(0).mass + collisions.particles(1).mass
        sticky_merged.radius = collisions.particles(0).radius
        for p1, p2, merged in zip(collisions.particles(0), collisions.particles(1), sticky_merged):
            merged.position = (p1 + p2).center_of_mass()
            merged.velocity = (p1 + p2).center_of_mass_velocity()

        print instance.model_time
        print instance.particles
        instance.particles.remove_particles(collisions.particles(0) + collisions.particles(1))
        instance.particles.add_particles(sticky_merged)

        instance.evolve_model(1.0 | nbody_system.time)
        print
        print instance.model_time
        print instance.particles
        self.assertTrue(collisions.is_set())
        self.assertTrue(instance.model_time < 1.0 | nbody_system.time)
        self.assertEquals(len(collisions.particles(0)), 1)
        self.assertEquals(len(collisions.particles(1)), 1)
        self.assertEquals(len(instance.particles - collisions.particles(0) - collisions.particles(1)), 2)
        self.assertEquals(abs(collisions.particles(0).x - collisions.particles(1).x) <
                (collisions.particles(0).radius + collisions.particles(1).radius),
                [True])
        instance.stop()

    def test9(self):
        print "Testing Bonsai2 tree build exception"
        plummer = new_plummer_model(50)
        #instance = self.new_instance_of_an_optional_code(Bonsai2, **default_options)
        instance = self.new_instance_of_an_optional_code(Bonsai2,debugger="gdb", **default_options)
        instance.particles.add_particles(plummer)

        instance.particles[0].position -= [1e9, 0, 0] | nbody_system.length
        self.assertRaises(AmuseException, instance.evolve_model, 1.0 | nbody_system.time, expected_message =
            "Error when calling 'evolve_model' of a 'Bonsai2', errorcode is -4, error is "
            "'The tree has become too deep, consider the removal of far away particles to prevent a too large box.'")

        instance.particles.remove_particle(instance.particles[0])
        instance.evolve_model(0.1 | nbody_system.time)
        instance.stop()
"""

ext: py

max_stars_repo_name: sibonyves/amuse

max_issues_repo_name: sibonyves/amuse

max_forks_repo_name: sibonyves/amuse

reason: import numpy


text:
from abc import ABC, abstractmethod
from typing import Union
import numpy as np

class DiscountCurve(ABC):
    """
    Base class for discounting curves (for any discounting term structure: rates, dividends, etc.)

    Note: will add date based functionality, day counting, etc.
    For now, this class is time based, to support common modeling problems. In the future, especially when dealing
    with interest rate products, it will support date based arguments and day counting conventions.
    """

    def implied_rate(self, T: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
        """
        Implied continuous rate at a particular time
        :param T: float or np.ndarray, time from which we imply the continuous rate, over [0,T]
        :return: float or np.ndarray (matches shape of input), the implied rate
        """
        return -np.log(self.discount_T(T)) / T

    @abstractmethod
    def discount_T(self, T: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
        """
        Discount at time T in the future
        :param T: float or np.ndarray, time(s) in the future
        :return: float or np.ndarray, discounts(s) at time(s) in the future
        """
        raise NotImplementedError

    def __call__(self, T: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
        """
        Discount at time T in the future
        :param T: float or np.ndarray, time(s) in the future
        :return: float or np.ndarray, discounts(s) at time(s) in the future
        """
        return self.discount_T(T)

class DiscountCurve_ConstRate(DiscountCurve):
    def __init__(self, rate: float):
        """
        Constant rate discount curve, exp(-r*T)
        :param rate: float, rate of discounting (e.g interest rate, div yield, etc)
        """
        super().__init__()
        self._r = rate

    def discount_T(self, T: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
        """
        Discount at time T in the future
        :param T: float or np.ndarray, time(s) in the future
        :return: float or np.ndarray, discounts(s) at time(s) in the future
        """
        return np.exp(-self._r * T)


text:
import os,sys
from HiNT.corelib import *
from multiprocessing import Pool
import numpy as np
import math

def getDivisionMatrix(mat1,mat2,fname):
	observedMatrix = np.genfromtxt(mat1,delimiter="\t")
	expectedMatrix = np.genfromtxt(mat2,delimiter="\t")
	expectedFinitevalues = expectedMatrix[np.isfinite(expectedMatrix)]
	observedFinitevalues = observedMatrix[np.isfinite(observedMatrix)]
	expected = expectedMatrix/np.mean(expectedFinitevalues)
	observed = observedMatrix/np.mean(observedFinitevalues)
	expected[np.where(expected==0)]=1
	with np.errstate(all='ignore'): divisionMatrix = observed/expected
	np.savetxt(fname,divisionMatrix,fmt="%.3f",delimiter='\t')
	return

def runBPcaller(params):
	#print params
	RscriptBPcaller,outputsubdir,matrixfile,tempbpoutputfile = params
	#print RscriptBPcaller,outputsubdir,matrixfile,tempbpoutputfile
	command = "Rscript %s %s %s %s"%(RscriptBPcaller,outputsubdir,matrixfile,tempbpoutputfile)
	print(command)
	run_cmd(command)
	chrompair_bps = open(tempbpoutputfile).readlines()
	os.remove(tempbpoutputfile)
	return chrompair_bps

def getAllRoughBreakpoints(matrix100kbInfo,background100kbInfo,rpInfo,outdir,name,cutoff,RscriptBPcaller,threads):
	DivisionMatrixInfo = {}
	outputsubdir = os.path.join(outdir,"AdjustedMatrixFiles100kb")
	if not os.path.isdir(outputsubdir):
		os.mkdir(outputsubdir)

	for chrompair in matrix100kbInfo:
		if chrompair not in background100kbInfo:
			pass
		else:
			if float(rpInfo[chrompair]) <= float(cutoff):
				mat1 = matrix100kbInfo[chrompair]
				mat2 = background100kbInfo[chrompair]
				fname = os.path.join(outputsubdir,name + '_' + chrompair + '_DivisionMatrix.txt')
				getDivisionMatrix(mat1,mat2,fname)
				DivisionMatrixInfo[chrompair] = fname
			else:
				pass
	#print DivisionMatrixInfo
	bpoutputfile = os.path.join(outdir,name + '_roughBP_100kb.txt')
	matrixfiles = os.listdir(outputsubdir)
	allparamsInfo = []
	for matrixfile in matrixfiles:
		matrixfile_prefix = matrixfile.replace('_DivisionMatrix.txt','')
		tempbpoutputfile = os.path.join(outdir,matrixfile_prefix + '_roughBP_100kb.txt')
		allparamsInfo.append([RscriptBPcaller,outputsubdir,matrixfile,tempbpoutputfile])
	#print allparamsInfo[0]
	results = []
	p = Pool(threads)
	result = p.map_async(runBPcaller, allparamsInfo, callback=results.append)
	p.close()
	p.join()

	validResults = []
	for line in results[0]:
		if line[1].strip() == 'NA':
			pass
		else:
			validResults.append(line)
	if len(validResults) > 0:
		with open(bpoutputfile,'w') as outf:
			for i,line in enumerate(validResults):
				if i == 0:
					outf.write(line[0])
					outf.write(line[1])
				else:
					outf.write(line[1])
	else:
		Info("No valid breakpoints detected!")
		sys.exit(0)
	return bpoutputfile,DivisionMatrixInfo

##################################################################################################
##############    scripts below this line are doing the breakpoints filtering   ##################

def getchromsize(chromlengthf,resolution):
	chromSizeInfo = {}
	inf = open(chromlengthf)
	for line in inf:
		line = line.strip().split('\t')
		chrom,size = line
		chromSizeInfo[chrom] = [int(size),int(int(size)/int(resolution)+1)]
	return chromSizeInfo

def readBPs(breakpointsf,chromSizeInfo):
	bpInfo = {}
	inf = open(breakpointsf)
	for line in inf:
		if line.startswith('chrom'):
			pass
		else:
			line = line.strip().split('\t')
			chrom,bprow,maxrow,bpcol,maxcol = line
			chroms = chrom.split('_')
			chr1,chr2 = chroms
			size1 = chromSizeInfo[chr1][0]
			size2 = chromSizeInfo[chr2][0]
			if size1 > size2:
				bpInfo[chrom] = [(chr1,bpcol,maxcol),(chr2,bprow,maxrow)]
			else:
				bpInfo[chrom] = [(chr1,bprow,maxrow),(chr2,bpcol,maxcol)]
	return bpInfo

def mergeValidBPs(validBPs,matrixfile):
	matrix = np.loadtxt(matrixfile)
	rowsum = np.nansum(matrix,axis=1)
	colsum = np.nansum(matrix,axis=0)

	tempBPs = []
	for validbp in validBPs:
		if validbp not in tempBPs:
			tempBPs.append(validbp)
		else:
			pass
	bppairs = [(i[1],i[3]) for i in tempBPs]
	chrompair = list(set([(i[0],i[2]) for i in tempBPs]))
	if len(chrompair) != 1:
		print(chrompair)

	#print tempBPs,chrompair

	bppairs.sort()
	bp1s = [int(i[0]) for i in bppairs]
	bp2s = [int(i[1]) for i in bppairs]

	differences1 = ['NA']+[(bp1s[i+1]-bp1s[i]) for i in range(len(bp1s)-1)]
	differences2 = ['NA']+[(bp2s[i+1]-bp2s[i]) for i in range(len(bp2s)-1)]

	#print differences1,differences2
	tempMerged = [i for i in bppairs]
	#print tempMerged
	for j in range(1,len(differences1)):
		if differences1[j] == 0 and differences2[j] == 1:
			#print colsum[bp2s[j-1]],colsum[bp2s[j]]
			if colsum[bp2s[j-1]] > colsum[bp2s[j]]:
				tempMerged.remove(bppairs[j])
			else:
				tempMerged.remove(bppairs[j-1])
		elif differences1[j] == 1 and differences2[j] == 0:
			#print rowsum[bp2s[j-1]] > rowsum[bp2s[j]]
			if rowsum[bp2s[j-1]] > rowsum[bp2s[j]]:
				tempMerged.remove(bppairs[j])
			else:
				tempMerged.remove(bppairs[j-1])
		else:
			pass

	#print tempMerged

	mergedValidBPs = []
	for bp in tempMerged:
		allinfo = [chrompair[0][0],str(bp[0]),chrompair[0][1],str(bp[1])]
		mergedValidBPs.append(allinfo)

	return mergedValidBPs

def mergeBPs(bps,matrix,axis,windowsize=5):
	bps = [int(i)-1 for i in bps] #python index starts with 0
	bps.sort()
	mergedbps = [str(bps[0])]
	for i in range(1,len(bps)):
		a = int(bps[i])
		b = int(mergedbps[-1])
		#print a,b
		#a is always larger than b
		if (a - b) <= windowsize:
			if axis == "row":
				suma = np.nansum(matrix[a,:])
				sumar = np.nansum(matrix[a+1,:])
				sumb = np.nansum(matrix[b,:])
				sumbr = np.nansum(matrix[b+1,:])
			elif axis == 'column':
				suma = np.nansum(matrix[:,a])
				sumar = np.nansum(matrix[:,a+1])
				sumb = np.nansum(matrix[:,b])
				sumbr = np.nansum(matrix[:,b+1])
			else:
				print("give the correct axis, row or column")
			if suma > sumar and sumb < sumbr: # __|--|__
				fca = suma/sumar
				fcb = sumbr/sumb
				if suma > sumbr and fca > fcb:
					if axis == "row":
						print("merge",np.nansum(matrix[a-1,:]), suma, sumar, np.nansum(matrix[b-1,:]),sumb, sumbr, fca, fcb)
					else:
						print("merge",np.nansum(matrix[:,a-1]), suma, sumar, np.nansum(matrix[:,b-1]),sumb, sumbr, fca, fcb)
					mergedbps.append(str(a))
					mergedbps.remove(str(b))
				elif suma < sumbr and fca > fcb:
					mergedbps.append(str(a))
				else:
					continue
		else:
			mergedbps.append(str(a))
	return mergedbps

def relativefold(a,b):
	if float(b) != 0:
		rfold = (float(a) - float(b))/float(b)
	else:
		rfold = (float(a) - float(b))/(float(b)+1)
	return rfold

def filtering(bps,matrixfile,chromSizeInfo,windowsize):
	matrix = np.loadtxt(matrixfile)
	rowsum = np.nansum(matrix,axis=1)
	#print np.shape(rowsum)
	colsum = np.nansum(matrix,axis=0)
	#print np.shape(colsum)
	chrom1Info,chrom2Info = bps
	chr1,bp1,max1 = chrom1Info
	chr2,bp2,max2 = chrom2Info
	bp1 = bp1.split(',')
	bp2 = bp2.split(',')

	bp1s = bp1
	bp2s = bp2
	bp1s = list(set(bp1s))
	bp2s = list(set(bp2s))

	#print bp1s,bp2s
	mergedbp1s = mergeBPs(bp1s,matrix,'row')
	#print mergedbp1s
	if str(int(max1)-1) not in mergedbp1s:
		mergedbp1s = mergedbp1s + [str(int(max1)-1)]

	mergedbp2s = mergeBPs(bp2s,matrix,'column')
	#print mergedbp2s

	if str(int(max2)-1) not in mergedbp2s:
		mergedbp2s = mergedbp2s + [str(int(max2)-1)]

	print(mergedbp1s,mergedbp2s)

	validBPs = []
	binsize1 = chromSizeInfo[chr1][1]
	binsize2 = chromSizeInfo[chr2][1]

	cutoff = np.nanpercentile(matrix,99)
	#print "cutoff:", cutoff

	for bp1 in mergedbp1s:
		for bp2 in mergedbp2s:
			sm1 = np.nanmedian(matrix[int(bp1)-5:int(bp1),int(bp2)+1:int(bp2)+6])
			sm2 = np.nanmedian(matrix[int(bp1)+1:int(bp1)+6,int(bp2)+1:int(bp2)+6])
			sm3 = np.nanmedian(matrix[int(bp1)-5:int(bp1),int(bp2)-5:int(bp2)])
			sm4 = np.nanmedian(matrix[int(bp1)+1:int(bp1)+6,int(bp2)-5:int(bp2)])
			#print sm1,sm2,sm3,sm4
			sms = np.sort(np.asarray([sm1,sm2,sm3,sm4]))
			translocationType = "unbalanced"
			if sms[2] > min(np.nanpercentile(matrix,80),2) and sms[1] < max(np.nanpercentile(matrix,60),1.5):
				translocationType = "balanced"
			#print translocationType

			if translocationType == "NULL":
				continue
			else:

				if translocationType == "balanced":
					if sm3 > sm4: #bp1 belongs to the left
					#if rowsum[int(bp1)] > rowsum[int(bp1)+1]:
						x1 = int(bp1)-windowsize+1
						x2 = int(bp1)+1
					else:
						x1 = int(bp1)-windowsize
						x2 = int(bp1)
				if translocationType == "unbalanced":
					if np.nanmean(rowsum[int(bp1)-5:int(bp1)]) > np.nanmean(rowsum[int(bp1)+1:int(bp1)+6]):
						x1 = int(bp1)-windowsize+1
						x2 = int(bp1)+1
					else:
						x1 = int(bp1)-windowsize
						x2 = int(bp1)

				if translocationType == "balanced":
					if sm3 > sm1:
						#if colsum[int(bp2)] > colsum[int(bp2)+1]:
						y1 = int(bp2)-windowsize+1
						y2 = int(bp2)+1
					else:
						y1 = int(bp2)-windowsize
						y2 = int(bp2)

				if translocationType == "unbalanced":
					if np.nanmean(colsum[int(bp2)-5:int(bp2)]) > np.nanmean(colsum[int(bp2)+1:int(bp2)+6]):
						y1 = int(bp2)-windowsize+1
						y2 = int(bp2)+1
					else:
						y1 = int(bp2)-windowsize
						y2 = int(bp2)

				if x1 < 0:
					x1 = 0
				if y1 < 0:
					y1 = 0
				#print x1,x2,y1,y2
				r1 = matrix[x1:x2,y1:y2]
				mr1 = np.nanmean(r1)

				if translocationType == "balanced":
					if sm1 > sm2:
						x1 = int(bp1)+1
						x2 = int(bp1)+1+windowsize
					else:
						x1 = int(bp1)
						x2 = int(bp1)+windowsize
				if translocationType == "unbalanced":
					if np.nanmean(rowsum[int(bp1)-5:int(bp1)]) > np.nanmean(rowsum[int(bp1)+1:int(bp1)+6]):
						x1 = int(bp1)+1
						x2 = int(bp1)+1+windowsize
					else:
						x1 = int(bp1)
						x2 = int(bp1)+windowsize
				if translocationType == "balanced":
					if sm4 > sm2:
						y1 = int(bp2)+1
						y2 = int(bp2)+1+windowsize
					else:
						y1 = int(bp2)
						y2 = int(bp2)+windowsize
				if translocationType == "unbalanced":
					if np.nanmean(colsum[int(bp2)-5:int(bp2)]) > np.nanmean(colsum[int(bp2)+1:int(bp2)+6]):
						y1 = int(bp2)+1
						y2 = int(bp2)+1+windowsize
					else:
						y1 = int(bp2)
						y2 = int(bp2)+windowsize
				if x2 >= binsize1:
					x2 = binsize1
				if y2 >= binsize2:
					y1 = binsize2
				#print x1,x2,y1,y2
				r2 = matrix[x1:x2,y1:y2]
				mr2 = np.nanmean(r2)

				if translocationType == "balanced":
					if sm1 > sm2:
						x1 = int(bp1)-windowsize+1
						x2 = int(bp1)+1
					else:
						x1 = int(bp1)-windowsize
						x2 = int(bp1)
				if translocationType == "unbalanced":
					if np.nanmean(rowsum[int(bp1)-5:int(bp1)]) > np.nanmean(rowsum[int(bp1)+1:int(bp1)+6]):
						x1 = int(bp1)-windowsize+1
						x2 = int(bp1)+1
					else:
						x1 = int(bp1)-windowsize
						x2 = int(bp1)
				if translocationType == "balanced":
					if sm3 > sm1:
						y1 = int(bp2)+1
						y2 = int(bp2)+1+windowsize
					else:
						y1 = int(bp2)
						y2 = int(bp2)+windowsize
				if translocationType == "unbalanced":
					if np.nanmean(colsum[int(bp2)-5:int(bp2)]) > np.nanmean(colsum[int(bp2)+1:int(bp2)+6]):
						y1 = int(bp2)+1
						y2 = int(bp2)+1+windowsize
					else:
						y1 = int(bp2)
						y2 = int(bp2)+windowsize
				if x1 < 0:
					x1 = 0
				if y2 > binsize2:
					y1 = binsize2
				#print x1,x2,y1,y2
				r3 = matrix[x1:x2,y1:y2]
				mr3 = np.nanpercentile(r3,80,interpolation="lower")
				mr3 = np.nanmean(r3)

				if translocationType == "balanced":
					if sm3 > sm4:
						x1 = int(bp1)+1
						x2 = int(bp1)+1+windowsize
					else:
						x1 = int(bp1)
						x2 = int(bp1)+windowsize
				if translocationType == "unbalanced":
					if np.nanmean(rowsum[int(bp1)-5:int(bp1)]) > np.nanmean(rowsum[int(bp1)+1:int(bp1)+6]):
						x1 = int(bp1)+1
						x2 = int(bp1)+1+windowsize
					else:
						x1 = int(bp1)
						x2 = int(bp1)+windowsize
				if translocationType == "balanced":
					if sm4 > sm2:
						y1 = int(bp2)-windowsize+1
						y2 = int(bp2)+1
					else:
						y1 = int(bp2)-windowsize
						y2 = int(bp2)
				if translocationType == "unbalanced":
					if np.nanmean(colsum[int(bp2)-5:int(bp2)]) > np.nanmean(colsum[int(bp2)+1:int(bp2)+6]):
						y1 = int(bp2)-windowsize+1
						y2 = int(bp2)+1
					else:
						y1 = int(bp2)-windowsize
						y2 = int(bp2)

				if x2 >= binsize1:
					x2 = binsize1
				if y1 < 0:
					y1 = 0
				#print x1,x2,y1,y2
				r4 = matrix[x1:x2,y1:y2]
				#print(np.shape(r4))
				mr4 = np.nanpercentile(r4,80,interpolation="lower")
				mr4 = np.nanmean(r4)
				###   mr1 | mr3
				###   ----|-----
				###   mr4 | mr2

				#print mr1,mr2,mr3,mr4,cutoff
				count = 0
				mrs = [mr1,mr2,mr3,mr4]
				for mr in mrs:
					if np.isnan(mr):
						pass
					else:
						if mr > cutoff:
							count += 1
				#print count

				if count == 0:
					pass
				elif count >= 3:
					pass
				elif count == 2: #balanced translocation, need a submatrix to calculate the rowsums and colsums
					mrs.sort()
					if (mr1 > cutoff and mr2 > cutoff) or (mr3> cutoff and mr4 > cutoff):
						rfold1 = relativefold(mrs[-2],mrs[-3])
						#print rfold1
						if rfold1 > 5:
							#print 'TRUE',bp1,bp2
							submatrix = matrix[:,0:int(bp2)]
							subrowsum = np.nansum(submatrix,1)
							submatrix = matrix[0:int(bp1),:]
							subcolsum = np.nansum(submatrix,0)

							if subrowsum[int(bp1)] > subrowsum[int(bp1)+1]:
								correctbp1 = bp1
							else:
								correctbp1 = str(int(bp1)+1)
							if subcolsum[int(bp2)] > subcolsum[int(bp2)+1]:
								correctbp2 = bp2
							else:
								correctbp2 = str(int(bp2)+1)

							validBPs.append([chr1,str(int(correctbp1)+1),chr2,str(int(correctbp2)+1)])
					else:
						pass
				else:
					#print 'count: ', count
					sortedmrs = sorted(mrs, key = lambda x : float('-inf') if math.isnan(x) else x)
					#print sortedmrs
					rfold1 = 1
					rfold2 = 1
					if not math.isnan(sortedmrs[-2]) and sortedmrs[-2] != 0:
						rfold1 = relativefold(sortedmrs[-1],sortedmrs[-2])
						rfold2 = 1
						#print 'rfold1: ', rfold1
					else:
						rfold2 = relativefold(sortedmrs[-1],cutoff)
						#print 'rfold2: ', rfold2
						rfold1 = 1

					if rfold1 > 5 or rfold2 > 2:
						#print 'TRUE',bp1,bp2
						if rowsum[int(bp1)] > rowsum[int(bp1)+1]:
							correctbp1 = bp1
						else:
							correctbp1 = str(int(bp1)+1)
						if colsum[int(bp2)] > colsum[int(bp2)+1]:
							correctbp2 = bp2
						else:
							correctbp2 = str(int(bp2)+1)
						validBPs.append([chr1,str(int(correctbp1)+1),chr2,str(int(correctbp2)+1)])

	return validBPs

def getValidRoughBP(chromlengthf,matrix100kbInfo,background100kbInfo,rpInfo,outdir,name,cutoff,RscriptBPcaller,threads):
	bpoutputfile,DivisionMatrixInfo = getAllRoughBreakpoints(matrix100kbInfo,background100kbInfo,rpInfo,outdir,name,cutoff,RscriptBPcaller,threads)
	chromSizeInfo = getchromsize(chromlengthf,resolution = 100000)
	bpInfo = readBPs(bpoutputfile,chromSizeInfo)
	windowsize = 10
	validbpoutputfile = os.path.join(outdir,name + '_roughBP_100kb_filtered.txt')
	outf = open(validbpoutputfile,'w')
	#print DivisionMatrixInfo.keys()
	for chrompair in bpInfo:
		if chrompair in DivisionMatrixInfo:
			#print chrompair
			bps = bpInfo[chrompair]
			matrixfile = DivisionMatrixInfo[chrompair]
			validBPs = filtering(bps,matrixfile,chromSizeInfo,windowsize)
			if len(validBPs) == 1:
				mergedValidBPs = validBPs
			else:
				mergedValidBPs = mergeValidBPs(validBPs,matrixfile)
			for validbp in mergedValidBPs:
				res = [chrompair] + validbp
				outf.write('\t'.join(res) + '\n')
		else:
			print("This chromosomal pair (%s) is not valid"%chrompair)
	outf.close()
	return validbpoutputfile

def validbpSummary(validbpoutputfile,rpInfo):
	resolution = 100000
	inf = open(validbpoutputfile)
	output = validbpoutputfile.rstrip('.txt') + '_summarized.txt'
	outf = open(output,'w')
	title = ["chrompair","chrom1","region1","chrom2","region2","rankProduct"]
	outf.write('\t'.join(title) + '\n')
	allInfo = {}
	for line in inf:
		line = line.strip().split('\t')
		chrompair,chrom1,bp1,chrom2,bp2 = line
		rp = float(rpInfo[chrompair])
		allInfo[chrompair] = [chrom1,bp1,chrom2,bp2,rp]

	sortedInfo = sorted(allInfo.items(), key=lambda x: x[1][4], reverse=False)
	for item in sortedInfo:
		chrompair = item[0]
		chrom1,bp1,chrom2,bp2,rp = item[1]
		region1 = [str((int(bp1)-1)*resolution),str((int(bp1)+1)*resolution)]
		region2 = [str((int(bp2)-1)*resolution),str((int(bp2)+1)*resolution)]
		newline = [chrompair,chrom1,'-'.join(region1),chrom2,'-'.join(region2),str(rp)]
		outf.write('\t'.join(newline) + '\n')
	outf.close()
	return output

max_stars_repo_name: suwangbio/HiNT

max_issues_repo_name: suwangbio/HiNT

max_forks_repo_name: suwangbio/HiNT


text:
#!/usr/bin/env python
# coding: utf8
#
# Copyright (c) 2022 Centre National d'Etudes Spatiales (CNES).
#
# This file is part of Shareloc
# (see https://github.com/CNES/shareloc).
#
# 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.
#
"""
Image class to handle Image data.
Shareloc Reference raster image input based on rasterio.
"""
# pylint: disable=no-member

# Standard imports
import logging

# Third party imports
import numpy as np
import rasterio
from affine import Affine

# pylint: disable=too-many-instance-attributes
class Image:
    """
    class Image to handle image data
    Shareloc reference for raster image input based on rasterio.
    """

    def __init__(self, image_path, read_data=False, roi=None, roi_is_in_physical_space=False):
        """
        constructor

        :param image_path : image path
        :type image_path  : string or None
        :param read_data  : read image data
        :type read_data  : bool
        :param roi  : region of interest [row_min,col_min,row_max,col_max] or [ymin,xmin,ymax,xmax] if
             roi_is_in_physical_space activated
        :type roi  : list
        :param roi_is_in_physical_space  : roi value in physical space
        :type roi_is_in_physical_space  : bool
        """
        if image_path is not None:
            # Image path
            self.image_path = image_path

            # Rasterio dataset
            self.dataset = rasterio.open(image_path)

            # Geo-transform of type Affine with convention :
            # | pixel size col,   row rotation, origin col |
            # | col rotation  , pixel size row, origin row |
            self.transform = self.dataset.transform

            roi_window = None
            if roi is not None:
                if roi_is_in_physical_space:
                    row_off, col_off = self.transform_physical_point_to_index(roi[0], roi[1])
                    row_max, col_max = self.transform_physical_point_to_index(roi[2], roi[3])
                    # index is relative to pixel center, here the roi is defined with corners
                    row_off += 0.5
                    col_off += 0.5
                    row_max += 0.5
                    col_max += 0.5
                    # in case of negative pixel size y
                    if row_off > row_max:
                        row_max, row_off = row_off, row_max
                    row_off = np.floor(row_off)
                    col_off = np.floor(col_off)
                    width = int(np.ceil(col_max - col_off))
                    height = int(np.ceil(row_max - row_off))
                    # pylint: disable=logging-too-many-args
                    logging.info("roi in image , offset : %s %s size %s %s", col_off, row_off, width, height)
                else:
                    row_off = roi[0]
                    col_off = roi[1]
                    width = roi[3] - roi[1]
                    height = roi[2] - roi[0]
                roi_window = rasterio.windows.Window(col_off, row_off, width, height)
                self.transform = self.dataset.window_transform(roi_window)
                self.nb_rows = height
                self.nb_columns = width
            else:
                # Image size
                self.nb_rows = self.dataset.height
                self.nb_columns = self.dataset.width
                roi_window = None

            # Georeferenced coordinates of the upper-left origin
            self.origin_row = self.transform[5]
            self.origin_col = self.transform[2]

            # Pixel size
            self.pixel_size_row = self.transform[4]
            self.pixel_size_col = self.transform[0]

            # row/col rotation
            self.pixel_rotation_row = self.transform[3]
            self.pixel_rotation_col = self.transform[1]

            if self.dataset.crs is not None:
                self.epsg = self.dataset.crs.to_epsg()
            else:
                self.epsg = None

            self.nodata = self.dataset.nodata

            self.mask = None
            self.data = None
            if read_data:
                # Data of shape (nb band, nb row, nb col)
                self.data = np.squeeze(self.dataset.read(window=roi_window))
                if self.nodata is not None:
                    self.mask = np.squeeze(self.dataset.read_masks(window=roi_window))
                    # pylint: disable=logging-too-many-args
                    logging.info("image contains %d nodata values ", np.sum(self.mask == 0))

    def set_metadata(self, nb_row, nb_col, nb_band, transform, datatype=np.float32):
        """
        Set metadata and create data (np.array filled with zeros of shape (nb_band, nb_row, nb_col))

        :param nb_row: number of row
        :type nb_row: int
        :param nb_col: number of col
        :type nb_col: int
        :param nb_band: number of band
        :type nb_band: int
        :param transform: geo-transform
        :type transform: 1D numpy array with convention :
                | pixel size col, row rotation, origin col, col rotation, pixel size row, origin row |
        :param datatype: data type
        :type datatype: numpy type
        """
        self.transform = Affine(transform[0], transform[1], transform[2], transform[3], transform[4], transform[5])

        # Georeferenced coordinates of the upper-left origin
        self.origin_row = transform[5]
        self.origin_col = transform[2]

        # Image size
        self.nb_rows = nb_row
        self.nb_columns = nb_col

        # Pixel size
        self.pixel_size_row = transform[4]
        self.pixel_size_col = transform[0]

        self.data = np.zeros((nb_band, nb_row, nb_col), dtype=datatype)

    def transform_index_to_physical_point(self, row, col):
        """
        Transform index to physical point

        :param row: row index
        :type row: float or 1D np.array
        :param col: col index
        :type col: float or 1D np.array
        :return: Georeferenced coordinates (row, col)
        :rtype: Tuple(georeference row float or 1D np.array, georeference col float or 1D np.array)
        """

        trans = self.transform
        col_geo, row_geo = trans * (col + 0.5, row + 0.5)
        return row_geo, col_geo

    def transform_physical_point_to_index(self, row_geo, col_geo):
        """
        Transform physical point to index

        :param row_geo: physical point row
        :type row_geo: float or 1D np.array
        :param col_geo: physical point col
        :type col_geo: float or 1D np.array
        :return: index coordinates (row, col)
        :rtype: Tuple(row float or 1D np.array, col float or 1D np.array)
        """
        trans_inv = ~self.transform
        col, row = trans_inv * (col_geo, row_geo)
        return row - 0.5, col - 0.5

max_stars_repo_name: CNES/shareloc

max_issues_repo_name: CNES/shareloc

max_forks_repo_name: CNES/shareloc


text:
# Assorted convenient functions
#
# Kaiyu Zheng
import numpy as np
import random

def remap(oldvalue, oldmin, oldmax, newmin, newmax):
    if oldmax - oldmin == 0:
        print("Warning in remap: the old range has size 0")
        oldmax = oldmin + oldvalue
    return (((oldvalue - oldmin) * (newmax - newmin)) / (oldmax - oldmin)) + newmin

# Plotting
import matplotlib.pyplot as plt

def plot_points(xvals, yvals, color=None,
                size=1.5, label=None, connected=True, style="--", linewidth=1.5,
                xlabel='x', ylabel='f(x)', loc="lower right"):
    if not connected:
        plt.scatter(xvals, yvals, s=size, c=color, label=label)
    else:
        plt.plot(xvals, yvals, style, linewidth=linewidth, label=label)
    # plt.axhline(y=0, color='k')
    # plt.axvline(x=0, color='k')
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.legend(loc=loc)

def save_plot(path, bbox_inches='tight'):
    plt.savefig(path, bbox_inches=bbox_inches)
    plt.close()

# colors
def lighter(color, percent):
    '''assumes color is rgb between (0, 0, 0) and (255, 255, 255)'''
    color = np.array(color)
    white = np.array([255, 255, 255])
    vector = white-color
    return color + vector * percent

def rgb_to_hex(rgb):
    r,g,b = rgb
    return '#%02x%02x%02x' % (int(r), int(g), int(b))

def hex_to_rgb(hx):
    """hx is a string, begins with #. ASSUME len(hx)=7."""
    if len(hx) != 7:
        raise ValueError("Hex must be #------")
    hx = hx[1:]  # omit the '#'
    r = int('0x'+hx[:2], 16)
    g = int('0x'+hx[2:4], 16)
    b = int('0x'+hx[4:6], 16)
    return (r,g,b)

def inverse_color_rgb(rgb):
    r,g,b = rgb
    return (255-r, 255-g, 255-b)

def inverse_color_hex(hx):
    """hx is a string, begins with #. ASSUME len(hx)=7."""
    return inverse_color_rgb(hex_to_rgb(hx))

def random_unique_color(colors, ctype=1):
    """
    ctype=1: completely random
    ctype=2: red random
    ctype=3: blue random
    ctype=4: green random
    ctype=5: yellow random
    """
    if ctype == 1:
        color = "#%06x" % random.randint(0x444444, 0x999999)
        while color in colors:
            color = "#%06x" % random.randint(0x444444, 0x999999)
    elif ctype == 2:
        color = "#%02x0000" % random.randint(0xAA, 0xFF)
        while color in colors:
            color = "#%02x0000" % random.randint(0xAA, 0xFF)
    elif ctype == 4:  # green
        color = "#00%02x00" % random.randint(0xAA, 0xFF)
        while color in colors:
            color = "#00%02x00" % random.randint(0xAA, 0xFF)
    elif ctype == 3:  # blue
        color = "#0000%02x" % random.randint(0xAA, 0xFF)
        while color in colors:
            color = "#0000%02x" % random.randint(0xAA, 0xFF)
    elif ctype == 5:  # yellow
        h = random.randint(0xAA, 0xFF)
        color = "#%02x%02x00" % (h, h)
        while color in colors:
            h = random.randint(0xAA, 0xFF)
            color = "#%02x%02x00" % (h, h)
    else:
        raise ValueError("Unrecognized color type %s" % (str(ctype)))
    return color

# Printing
def json_safe(obj):
    if isinstance(obj, bool):
        return str(obj).lower()
    elif isinstance(obj, (list, tuple)):
        return [json_safe(item) for item in obj]
    elif isinstance(obj, dict):
        return {json_safe(key):json_safe(value) for key, value in obj.items()}
    else:
        return str(obj)
    return obj

# Math
def vec(p1, p2):
    """ vector from p1 to p2 """
    if type(p1) != np.ndarray:
        p1 = np.array(p1)
    if type(p2) != np.ndarray:
        p2 = np.array(p2)
    return p2 - p1

def proj(vec1, vec2, scalar=False):
    # Project vec1 onto vec2. Returns a vector in the direction of vec2.
    scale = np.dot(vec1, vec2) / np.linalg.norm(vec2)
    if scalar:
        return scale
    else:
        return vec2 * scale

def R_x(th):
    return np.array([
        1, 0, 0, 0,
        0, np.cos(th), -np.sin(th), 0,
        0, np.sin(th), np.cos(th), 0,
        0, 0, 0, 1
    ]).reshape(4,4)

def R_y(th):
    return np.array([
        np.cos(th), 0, np.sin(th), 0,
        0, 1, 0, 0,
        -np.sin(th), 0, np.cos(th), 0,
        0, 0, 0, 1
    ]).reshape(4,4)

def R_z(th):
    return np.array([
        np.cos(th), -np.sin(th), 0, 0,
        np.sin(th), np.cos(th), 0, 0,
        0, 0, 1, 0,
        0, 0, 0, 1
    ]).reshape(4,4)

def T(dx, dy, dz):
    return np.array([
        1, 0, 0, dx,
        0, 1, 0, dy,
        0, 0, 1, dz,
        0, 0, 0, 1
    ]).reshape(4,4)

def to_radians(th):
    return th*np.pi / 180

def R_between(v1, v2):
    if len(v1) != 3 or len(v2) != 3:
        raise ValueError("Only applicable to 3D vectors!")
    v = np.cross(v1, v2)
    c = np.dot(v1, v2)
    s = np.linalg.norm(v)
    I = np.identity(3)

    vX = np.array([
        0, -v[2], v[1],
        v[2], 0, -v[0],
        -v[1], v[0], 0
    ]).reshape(3,3)
    R = I + vX + np.matmul(vX,vX) * ((1-c)/(s**2))
    return R

def approx_equal(v1, v2, epsilon=1e-6):
    if len(v1) != len(v2):
        return False
    for i in range(len(v1)):
        if abs(v1[i] - v2[i]) > epsilon:
            return False
    return True

# Others
def safe_slice(arr, start, end):
    true_start = max(0, min(len(arr)-1, start))
    true_end = max(0, min(len(arr)-1, end))
    return arr[true_start:true_end]

# Colors on terminal https://stackoverflow.com/a/287944/2893053
class bcolors:
    WHITE = '\033[97m'
    CYAN = '\033[96m'
    MAGENTA = '\033[95m'
    BLUE = '\033[94m'
    GREEN = '\033[92m'
    YELLOW = '\033[93m'
    RED = '\033[91m'
    BOLD = '\033[1m'
    ENDC = '\033[0m'

    @staticmethod
    def disable():
        bcolors.WHITE   = ''
        bcolors.CYAN    = ''
        bcolors.MAGENTA = ''
        bcolors.BLUE    = ''
        bcolors.GREEN   = ''
        bcolors.YELLOW  = ''
        bcolors.RED     = ''
        bcolors.ENDC    = ''

    @staticmethod
    def s(color, content):
        """Returns a string with color when shown on terminal.
        `color` is a constant in `bcolors` class."""
        return color + content + bcolors.ENDC

def print_info(content):
    print(bcolors.s(bcolors.BLUE, content))

def print_note(content):
    print(bcolors.s(bcolors.YELLOW, content))

def print_error(content):
    print(bcolors.s(bcolors.RED, content))

def print_warning(content):
    print(bcolors.s(bcolors.YELLOW, content))

def print_success(content):
    print(bcolors.s(bcolors.GREEN, content))

def print_info_bold(content):
    print(bcolors.BOLD + bcolors.s(bcolors.BLUE, content))

def print_note_bold(content):
    print(bcolors.BOLD + bcolors.s(bcolors.GREEN, content))

def print_error_bold(content):
    print(bcolors.BOLD + bcolors.s(bcolors.RED, content))

def print_warning_bold(content):
    print(bcolors.BOLD + bcolors.s(bcolors.YELLOW, content))

def print_success_bold(content):
    print(bcolors.BOLD + bcolors.s(bcolors.GREEN, content))
# For your convenience:
# from pomdp_py.util import print_info, print_error, print_warning, print_success, print_info_bold, print_error_bold, print_warning_bold, print_success_bold

max_stars_repo_name: beckymathew/pomdp-py

max_issues_repo_name: beckymathew/pomdp-py

max_forks_repo_name: beckymathew/pomdp-py


text:
# Copyright 2019 Xilinx Inc.
#
# 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.

# This code is based on: https://github.com/nutonomy/second.pytorch.git
#
# MIT License
# Copyright (c) 2018
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

import time
from enum import Enum

import numpy as np
import torch
from torch import nn
from torch.nn import functional as F

from utils import box_torch_ops
from models.pointpillars import PillarFeatureNet, PointPillarsScatter, RPN
from utils.torchprofiler import Profiler

import os
if os.environ["W_QUANT"]=='1':
    from pytorch_nndct.nn import QuantStub, DeQuantStub

class LossNormType(Enum):
    NormByNumPositives = "norm_by_num_positives"
    NormByNumExamples = "norm_by_num_examples"
    NormByNumPosNeg = "norm_by_num_pos_neg"

def get_paddings_indicator(actual_num, max_num, axis=0):
    """Create boolean mask by actually number of a padded tensor.

    Args:
        actual_num ([type]): [description]
        max_num ([type]): [description]

    Returns:
        [type]: [description]
    """

    actual_num = torch.unsqueeze(actual_num, axis + 1)
    # tiled_actual_num: [N, M, 1]
    max_num_shape = [1] * len(actual_num.shape)
    max_num_shape[axis + 1] = -1
    max_num = torch.arange(
        max_num, dtype=torch.int, device=actual_num.device).view(max_num_shape)
    # tiled_actual_num: [[3,3,3,3,3], [4,4,4,4,4], [2,2,2,2,2]]
    # tiled_max_num: [[0,1,2,3,4], [0,1,2,3,4], [0,1,2,3,4]]
    paddings_indicator = actual_num.int() > max_num
    # paddings_indicator shape: [batch_size, max_num]
    return paddings_indicator

class PreProcess(object):
    def __init__(self,
                 with_distance=False,
                 voxel_size=(0.2, 0.2, 4),
                 pc_range=(0, -40, -3, 70.4, 40, 1)):
        # Need pillar (voxel) size and x/y offset in order to calculate pillar offset
        self.vx = voxel_size[0]
        self.vy = voxel_size[1]
        self.x_offset = self.vx / 2 + pc_range[0]
        self.y_offset = self.vy / 2 + pc_range[1]
        self._with_distance = with_distance

        # scale x,y,z according to the range
        self.vz = voxel_size[2]
        self.pc_range = pc_range

    def __call__(self, features, num_voxels, coors):
        device = features.device
        '''
        # Find distance of x, y, and z from cluster center
        points_mean = features[:, :, :3].sum(dim=1, keepdim=True) / num_voxels.type_as(features).view(-1, 1, 1)
        f_cluster = features[:, :, :3] - points_mean

        # Find distance of x, y, and z from pillar center
        f_center = torch.zeros_like(features[:, :, :2])
        f_center[:, :, 0] = features[:, :, 0] - (coors[:, 3].float().unsqueeze(1) * self.vx + self.x_offset)
        f_center[:, :, 1] = features[:, :, 1] - (coors[:, 2].float().unsqueeze(1) * self.vy + self.y_offset)

        # scale feature according to the range
        half_voxel_len = torch.Tensor([self.vx, self.vy, self.vz]).to(device) / 2.
        voxel_len = torch.Tensor([self.vx, self.vy, self.vz]).to(device)
        f_cluster = (f_cluster + half_voxel_len) / voxel_len
        f_center[:, :, 0] = (f_center[:, :, 0] + half_voxel_len[0]) / self.vx
        f_center[:, :, 1] = (f_center[:, :, 1] + half_voxel_len[1]) / self.vy
        '''
        pc_start = torch.from_numpy(self.pc_range[:3]).to(device)
        pc_len = torch.from_numpy(self.pc_range[3:6] - self.pc_range[:3]).to(device)
        features[:, :, :3] = features[:, :, :3] - pc_start
        features[:, :, :3] = features[:, :, :3] / pc_len

        # Combine together feature decorations
        # features_ls = [features, f_cluster, f_center]
        features_ls = [features]

        if self._with_distance:
            points_dist = torch.norm(features[:, :, :3], 2, 2, keepdim=True)
            features_ls.append(points_dist)
        features = torch.cat(features_ls, dim=-1)

        # The feature decorations were calculated without regard to whether pillar was empty. Need to ensure that
        # empty pillars remain set to zeros.
        voxel_count = features.shape[1]
        mask = get_paddings_indicator(num_voxels, voxel_count, axis=0)
        mask = torch.unsqueeze(mask, -1).type_as(features)
        features *= mask

        features = features.permute(2, 0, 1).contiguous().unsqueeze(0)

        return features

class PostProcess(object):

    def __init__(self,
                 num_class=2,
                 use_direction_classifier=True,
                 use_sigmoid_score=False,
                 encode_background_as_zeros=True,
                 use_rotate_nms=True,
                 multiclass_nms=False,
                 nms_score_threshold=0.5,
                 nms_pre_max_size=1000,
                 nms_post_max_size=20,
                 nms_iou_threshold=0.1,
                 target_assigner=None):
        super().__init__()
        self._box_coder = target_assigner.box_coder
        self._num_class = num_class
        self._encode_background_as_zeros = encode_background_as_zeros
        self._use_direction_classifier = use_direction_classifier
        self._use_sigmoid_score = use_sigmoid_score
        self._use_rotate_nms = use_rotate_nms
        self._multiclass_nms = multiclass_nms
        self._nms_pre_max_size = nms_pre_max_size
        self._nms_post_max_size = nms_post_max_size
        self._nms_iou_threshold = nms_iou_threshold
        self._nms_score_threshold = nms_score_threshold

    def __call__(self, example, preds_dict):
        t = time.time()
        batch_size = example['anchors'].shape[0]
        batch_anchors = example["anchors"].view(batch_size, -1, 7)

        # self._total_inference_count += batch_size
        batch_rect = example["rect"]
        batch_Trv2c = example["Trv2c"]
        batch_P2 = example["P2"]
        if "anchors_mask" not in example:
            batch_anchors_mask = [None] * batch_size
        else:
            batch_anchors_mask = example["anchors_mask"].view(batch_size, -1)
        batch_imgidx = example['image_idx']

        batch_box_preds = preds_dict["box_preds"]
        batch_cls_preds = preds_dict["cls_preds"]
        batch_box_preds = batch_box_preds.view(batch_size, -1,
                                               self._box_coder.code_size)
        num_class_with_bg = self._num_class
        if not self._encode_background_as_zeros:
            num_class_with_bg = self._num_class + 1

        batch_cls_preds = batch_cls_preds.view(batch_size, -1, num_class_with_bg)
        batch_box_preds = self._box_coder.decode_torch(batch_box_preds, batch_anchors)
        if self._use_direction_classifier:
            batch_dir_preds = preds_dict["dir_cls_preds"]
            batch_dir_preds = batch_dir_preds.view(batch_size, -1, 2)
        else:
            batch_dir_preds = [None] * batch_size

        predictions_dicts = []
        for box_preds, cls_preds, dir_preds, rect, Trv2c, P2, img_idx, a_mask in zip(
                batch_box_preds, batch_cls_preds, batch_dir_preds, batch_rect,
                batch_Trv2c, batch_P2, batch_imgidx, batch_anchors_mask
        ):
            if a_mask is not None:
                box_preds = box_preds[a_mask]
                cls_preds = cls_preds[a_mask]
            if self._use_direction_classifier:
                if a_mask is not None:
                    dir_preds = dir_preds[a_mask]
                # print(dir_preds.shape)
                dir_labels = torch.max(dir_preds, dim=-1)[1]
            if self._encode_background_as_zeros:
                # this don't support softmax
                assert self._use_sigmoid_score is True
                total_scores = torch.sigmoid(cls_preds)
            else:
                # encode background as first element in one-hot vector
                if self._use_sigmoid_score:
                    total_scores = torch.sigmoid(cls_preds)[..., 1:]
                else:
                    total_scores = F.softmax(cls_preds, dim=-1)[..., 1:]
            # Apply NMS in birdeye view
            if self._use_rotate_nms:
                nms_func = box_torch_ops.rotate_nms
            else:
                nms_func = box_torch_ops.nms
            selected_boxes = None
            selected_labels = None
            selected_scores = None
            selected_dir_labels = None

            if self._multiclass_nms:
                # curently only support class-agnostic boxes.
                boxes_for_nms = box_preds[:, [0, 1, 3, 4, 6]]
                if not self._use_rotate_nms:
                    box_preds_corners = box_torch_ops.center_to_corner_box2d(
                        boxes_for_nms[:, :2], boxes_for_nms[:, 2:4],
                        boxes_for_nms[:, 4])
                    boxes_for_nms = box_torch_ops.corner_to_standup_nd(
                        box_preds_corners)
                boxes_for_mcnms = boxes_for_nms.unsqueeze(1)
                selected_per_class = box_torch_ops.multiclass_nms(
                    nms_func=nms_func,
                    boxes=boxes_for_mcnms,
                    scores=total_scores,
                    num_class=self._num_class,
                    pre_max_size=self._nms_pre_max_size,
                    post_max_size=self._nms_post_max_size,
                    iou_threshold=self._nms_iou_threshold,
                    score_thresh=self._nms_score_threshold,
                )
                selected_boxes, selected_labels, selected_scores = [], [], []
                selected_dir_labels = []
                for i, selected in enumerate(selected_per_class):
                    if selected is not None:
                        num_dets = selected.shape[0]
                        selected_boxes.append(box_preds[selected])
                        selected_labels.append(
                            torch.full([num_dets], i, dtype=torch.int64))
                        if self._use_direction_classifier:
                            selected_dir_labels.append(dir_labels[selected])
                        selected_scores.append(total_scores[selected, i])
                if len(selected_boxes) > 0:
                    selected_boxes = torch.cat(selected_boxes, dim=0)
                    selected_labels = torch.cat(selected_labels, dim=0)
                    selected_scores = torch.cat(selected_scores, dim=0)
                    if self._use_direction_classifier:
                        selected_dir_labels = torch.cat(
                            selected_dir_labels, dim=0)
                else:
                    selected_boxes = None
                    selected_labels = None
                    selected_scores = None
                    selected_dir_labels = None
            else:
                # get highest score per prediction, than apply nms
                # to remove overlapped box.
                if num_class_with_bg == 1:
                    top_scores = total_scores.squeeze(-1)
                    top_labels = torch.zeros(
                        total_scores.shape[0],
                        device=total_scores.device,
                        dtype=torch.long)
                else:
                    top_scores, top_labels = torch.max(total_scores, dim=-1)

                if self._nms_score_threshold > 0.0:
                    thresh = torch.tensor(
                        [self._nms_score_threshold],
                        device=total_scores.device).type_as(total_scores)
                    top_scores_keep = (top_scores >= thresh)
                    top_scores = top_scores.masked_select(top_scores_keep)
                if top_scores.shape[0] != 0:
                    if self._nms_score_threshold > 0.0:
                        box_preds = box_preds[top_scores_keep]
                        if self._use_direction_classifier:
                            dir_labels = dir_labels[top_scores_keep]
                        top_labels = top_labels[top_scores_keep]
                    boxes_for_nms = box_preds[:, [0, 1, 3, 4, 6]]
                    if not self._use_rotate_nms:
                        box_preds_corners = box_torch_ops.center_to_corner_box2d(
                            boxes_for_nms[:, :2], boxes_for_nms[:, 2:4],
                            boxes_for_nms[:, 4])
                        boxes_for_nms = box_torch_ops.corner_to_standup_nd(
                            box_preds_corners)
                    # the nms in 3d detection just remove overlap boxes.
                    selected = nms_func(
                        boxes_for_nms,
                        top_scores,
                        pre_max_size=self._nms_pre_max_size,
                        post_max_size=self._nms_post_max_size,
                        iou_threshold=self._nms_iou_threshold,
                    )
                else:
                    selected = None
                if selected is not None:
                    selected_boxes = box_preds[selected]
                    if self._use_direction_classifier:
                        selected_dir_labels = dir_labels[selected]
                    selected_labels = top_labels[selected]
                    selected_scores = top_scores[selected]
            # finally generate predictions.

            if selected_boxes is not None:
                box_preds = selected_boxes
                scores = selected_scores
                label_preds = selected_labels
                if self._use_direction_classifier:
                    dir_labels = selected_dir_labels
                    # opp_labels = (box_preds[..., -1] > 0) ^ dir_labels.byte()
                    opp_labels = (box_preds[..., -1] > 0) ^ dir_labels.bool()
                    box_preds[..., -1] += torch.where(
                        opp_labels,
                        torch.tensor(np.pi).type_as(box_preds),
                        torch.tensor(0.0).type_as(box_preds))
                    # box_preds[..., -1] += (
                    #     ~(dir_labels.byte())).type_as(box_preds) * np.pi
                final_box_preds = box_preds
                final_scores = scores
                final_labels = label_preds
                final_box_preds_camera = box_torch_ops.box_lidar_to_camera(
                    final_box_preds, rect, Trv2c)
                locs = final_box_preds_camera[:, :3]
                dims = final_box_preds_camera[:, 3:6]
                angles = final_box_preds_camera[:, 6]
                camera_box_origin = [0.5, 1.0, 0.5]
                box_corners = box_torch_ops.center_to_corner_box3d(
                    locs, dims, angles, camera_box_origin, axis=1)
                box_corners_in_image = box_torch_ops.project_to_image(
                    box_corners, P2)
                # box_corners_in_image: [N, 8, 2]
                minxy = torch.min(box_corners_in_image, dim=1)[0]
                maxxy = torch.max(box_corners_in_image, dim=1)[0]
                # minx = torch.min(box_corners_in_image[..., 0], dim=1)[0]
                # maxx = torch.max(box_corners_in_image[..., 0], dim=1)[0]
                # miny = torch.min(box_corners_in_image[..., 1], dim=1)[0]
                # maxy = torch.max(box_corners_in_image[..., 1], dim=1)[0]
                # box_2d_preds = torch.stack([minx, miny, maxx, maxy], dim=1)
                box_2d_preds = torch.cat([minxy, maxxy], dim=1)
                # predictions
                predictions_dict = {
                    "bbox": box_2d_preds,
                    "box3d_camera": final_box_preds_camera,
                    "box3d_lidar": final_box_preds,
                    "scores": final_scores,
                    "label_preds": label_preds,
                    "image_idx": img_idx,
                }
            else:
                predictions_dict = {
                    "bbox": None,
                    "box3d_camera": None,
                    "box3d_lidar": None,
                    "scores": None,
                    "label_preds": None,
                    "image_idx": img_idx,
                }
            predictions_dicts.append(predictions_dict)
        return predictions_dicts

class VoxelNet(nn.Module):
    def __init__(self,
                 output_shape,
                 num_class=2,
                 num_input_features=4,
                 vfe_class_name="PillarFeatureNet",
                 vfe_num_filters=[32, 128],
                 with_distance=False,
                 middle_class_name="PointPillarsScatter",
                 middle_num_filters_d1=[64],
                 middle_num_filters_d2=[64, 64],
                 rpn_class_name="RPN",
                 rpn_layer_nums=[3, 5, 5],
                 rpn_layer_strides=[2, 2, 2],
                 rpn_num_filters=[128, 128, 256],
                 rpn_upsample_strides=[1, 2, 4],
                 rpn_num_upsample_filters=[256, 256, 256],
                 use_norm=True,
                 use_groupnorm=False,
                 num_groups=32,
                 use_sparse_rpn=False,
                 use_direction_classifier=True,
                 use_sigmoid_score=False,
                 encode_background_as_zeros=True,
                 use_rotate_nms=True,
                 multiclass_nms=False,
                 nms_score_threshold=0.5,
                 nms_pre_max_size=1000,
                 nms_post_max_size=20,
                 nms_iou_threshold=0.1,
                 target_assigner=None,
                 use_bev=False,
                 lidar_only=False,
                 cls_loss_weight=1.0,
                 loc_loss_weight=1.0,
                 pos_cls_weight=1.0,
                 neg_cls_weight=1.0,
                 direction_loss_weight=1.0,
                 loss_norm_type=LossNormType.NormByNumPositives,
                 encode_rad_error_by_sin=False,
                 loc_loss_ftor=None,
                 cls_loss_ftor=None,
                 voxel_size=(0.2, 0.2, 4),
                 pc_range=(0, -40, -3, 70.4, 40, 1),
                 name='voxelnet'):
        super().__init__()
        self.name = name
        self._num_class = num_class
        self._use_rotate_nms = use_rotate_nms
        self._multiclass_nms = multiclass_nms
        self._nms_score_threshold = nms_score_threshold
        self._nms_pre_max_size = nms_pre_max_size
        self._nms_post_max_size = nms_post_max_size
        self._nms_iou_threshold = nms_iou_threshold
        self._use_sigmoid_score = use_sigmoid_score
        self._encode_background_as_zeros = encode_background_as_zeros
        self._use_sparse_rpn = use_sparse_rpn
        self._use_direction_classifier = use_direction_classifier
        self._use_bev = use_bev
        self._total_forward_time = 0.0
        self._total_postprocess_time = 0.0
        self._total_inference_count = 0
        self._num_input_features = num_input_features
        self._box_coder = target_assigner.box_coder
        self._lidar_only = lidar_only
        self.target_assigner = target_assigner
        self._pos_cls_weight = pos_cls_weight
        self._neg_cls_weight = neg_cls_weight
        self._encode_rad_error_by_sin = encode_rad_error_by_sin
        self._loss_norm_type = loss_norm_type
        # self._dir_loss_ftor = WeightedSoftmaxClassificationLoss()
        self._dir_loss_ftor = None

        self._loc_loss_ftor = loc_loss_ftor
        self._cls_loss_ftor = cls_loss_ftor
        self._direction_loss_weight = direction_loss_weight
        self._cls_loss_weight = cls_loss_weight
        self._loc_loss_weight = loc_loss_weight

        if os.environ["W_QUANT"]=='1':
            self.quant = QuantStub()
            self.dequant = DeQuantStub()

        self.voxel_feature_extractor = PillarFeatureNet(
            num_input_features,
            use_norm,
            num_filters=vfe_num_filters,
            with_distance=with_distance,
            voxel_size=voxel_size,
            pc_range=pc_range
        )

        self.middle_feature_extractor = PointPillarsScatter(output_shape=output_shape,
                                                            num_input_features=vfe_num_filters[-1])
        num_rpn_input_filters = self.middle_feature_extractor.nchannels

        self.rpn = RPN(
            use_norm=True,
            num_class=num_class,
            layer_nums=rpn_layer_nums,
            layer_strides=rpn_layer_strides,
            num_filters=rpn_num_filters,
            upsample_strides=rpn_upsample_strides,
            num_upsample_filters=rpn_num_upsample_filters,
            num_input_filters=num_rpn_input_filters,
            num_anchor_per_loc=target_assigner.num_anchors_per_location,
            encode_background_as_zeros=encode_background_as_zeros,
            use_direction_classifier=use_direction_classifier,
            use_bev=use_bev,
            use_groupnorm=use_groupnorm,
            num_groups=num_groups,
            box_code_size=target_assigner.box_coder.code_size)

        self.register_buffer("global_step", torch.LongTensor(1).zero_())

    def update_global_step(self):
        self.global_step += 1

    def get_global_step(self):
        return int(self.global_step.cpu().numpy()[0])

    def forward(self, aug_voxels, coors):

        batch_size_dev = 1

        if os.environ["W_QUANT"]=='1':
            aug_voxels = self.quant(aug_voxels)
        voxel_features = self.voxel_feature_extractor(aug_voxels)
        if os.environ["W_QUANT"]=='1':
            voxel_features = self.dequant(voxel_features)

        spatial_features = self.middle_feature_extractor(
            voxel_features, coors, batch_size_dev)

        if os.environ["W_QUANT"]=='1':
            spatial_features = self.quant(spatial_features)
        if self._use_bev:
            preds_dict = self.rpn(spatial_features, example["bev_map"])
        else:
            preds_dict = self.rpn(spatial_features)
        if os.environ["W_QUANT"]=='1':
            preds_dict = self.dequant(preds_dict)

        # (1)
        print("aug_voxels.shape: ", aug_voxels.shape)
        profiler = Profiler(self.voxel_feature_extractor.cpu())
        profiler.run_with_input(aug_voxels.cpu())
        profiler.print_summary()
        # (2)
        profiler = Profiler(self.middle_feature_extractor, device="cpu")
        profiler.run_with_input(voxel_features, coors, batch_size_dev)
        profiler.print_summary()
        # (3)
        print("spatial_features shape: ", spatial_features.shape)
        profiler = Profiler(self.rpn.cpu())
        one_spatial_features = spatial_features[0,...].unsqueeze(0)
        profiler.run_with_input(one_spatial_features.cpu())
        profiler.print_summary()
        print("Done")
        import sys
        sys.exit()

        return preds_dict

    def clear_time_metrics(self):
        self._total_forward_time = 0.0
        self._total_postprocess_time = 0.0
        self._total_inference_count = 0

    @staticmethod
    def convert_norm_to_float(net):
        '''
        BatchNorm layers to have parameters in single precision.
        Find all layers and convert them back to float. This can't
        be done with built in .apply as that function will apply
        fn to all modules, parameters, and buffers. Thus we wouldn't
        be able to guard the float conversion based on the module type.
        '''
        if isinstance(net, torch.nn.modules.batchnorm._BatchNorm):
            net.float()
        for child in net.children():
            VoxelNet.convert_norm_to_float(net)
        return net

max_stars_repo_name: guochunhe/Vitis-AI

max_issues_repo_name: guochunhe/Vitis-AI

max_forks_repo_name: guochunhe/Vitis-AI


text:
from imutils.video import FileVideoStream
import numpy as np
import imutils
import time
import cv2
import posenet
import tensorflow.compat.v1 as tf
tf.disable_v2_behavior()

def loading_model(graph, text):
	'''load the soccer ball detection Tensorflow model graph and text files into OpenCV net'''

	net = cv2.dnn.readNetFromTensorflow(graph, text)

	return net

def initialize_constants():
	'''Initialize all constants, keys, and counters'''

	body_part_key = {0:"Head", 1:"Left Thigh", 2:"Right Thigh", 3:"Left Foot", 4:"Right Foot"}
	body_part_bounces = {"Head":0, "Left Thigh":0, "Right Thigh":0, "Left Foot":0, "Right Foot":0}
	body_part_key_simple = {0:"Head", 1:"Thigh", 2:"Thigh", 3:"Foot", 4:"Foot", 5:"Ground"}
	body_part_sequence = []
	confidence = 0.0
	(W, H) = (300, 300)
	bounces = 0
	frame_num = 0
	last_juggle_frame = -4
	first_detection = True
	cYs = [0,0] #list of cY for each frame
	cXs = [0,0] #list of cX for each frame
	skipped_frames = 100
	ground = False

	return (body_part_key, body_part_bounces, body_part_key_simple, body_part_sequence, confidence, W, H, bounces, frame_num, last_juggle_frame, first_detection, cYs, cXs, skipped_frames, ground)

def read_frame(video, frame_num):
	'''Read the next frame of the video'''

	frames_per_read = 2
	for i in range(frames_per_read):
		if video.more():
			frame = video.read()
			frame_num += 1
	return (frame, frame_num)

def process_frame(frame, net, W, H):
	'''Runs the current frame through the detection model'''

	blob = cv2.dnn.blobFromImage(frame, size = (W, H), swapRB=True)
	#blobFromImage can perform mean subtraction, scaling, and channel swapping
	net.setInput(blob)
	detection = net.forward()

	return detection

def check_bounce(cYs):
	'''Checks to see if a juggle has been made by looking at position and direction history'''

	if (cYs[-1] <= cYs[-2]) and (cYs[-2] >= cYs[-3]):
		return True
	else:
		return False

def get_centroid(frame, detection, W, H):
	'''Determine the centroid and diameter of the detected ball by looking at bounding box'''
	box = detection[0, 0, 0, 3:7] * np.array([W, H, W, H])

	(startX, startY, endX, endY) = box.astype("int")
	diameter = ((endX - startX) + (endY - startY)) / 2.0
	#cv2.rectangle(frame, (startX, startY), (endX, endY), (0, 255, 0), 2)

	cX = int((startX + endX) / 2.0)
	cY = int((startY + endY) / 2.0)

	return cX, cY, diameter

def get_pose(frame, sess, output_stride, model_outputs, scale_factor = 1):
	'''Passes frame through Posenet to determine the coordinates of key points of body'''

	input_image, display_image, output_scale = posenet.utils._process_input(frame, scale_factor, output_stride)

	heatmaps_result, offsets_result, displacement_fwd_result, displacement_bwd_result = sess.run(
		model_outputs,
		feed_dict={'image:0': input_image}
	)

	pose_scores, keypoint_scores, keypoint_coords = posenet.decode_multiple_poses(
		heatmaps_result.squeeze(axis=0),
		offsets_result.squeeze(axis=0),
		displacement_fwd_result.squeeze(axis=0),
		displacement_bwd_result.squeeze(axis=0),
		output_stride=output_stride,
		max_pose_detections=1,
		min_pose_score=0.15)

	keypoint_coords *= output_scale

	#drawing nodes and lines on frame
	draw_image = posenet.draw_skel_and_kp(
		display_image, pose_scores, keypoint_scores, keypoint_coords,
		min_pose_score=0.15, min_part_score=0.05)

	return (draw_image, keypoint_scores[0], keypoint_coords[0])

def get_closest_body_part(k_scores, k_coords, cX, cY):
	'''Determine which body part is closest to bell centroid (but not above)'''

	# PART_NAMES = ["nose", "leftEye", "rightEye", "leftEar", "rightEar", "leftShoulder",
 	#    "rightShoulder", "leftElbow", "rightElbow", "leftWrist", "rightWrist",
 	#    "leftHip", "rightHip", "leftKnee", "rightKnee", "leftAnkle", "rightAnkle"]

	#approximate length of thigh. This will be used to help determine if the touch is not part of the bady and instead is the ground or another object.
	knee_distance = ((k_coords[11, 1] - k_coords[13, 1])**2 + (k_coords[11, 0] - k_coords[13, 0])**2)**0.5

	#choosing hips instead of knees works better for detecting thigh vs foot
	parts = [0, 11, 12, 15, 16]
	distances = []

	#loop over five body parts that can make a touch
	for part in parts:
		#if the body part is above the ball, then penalize it so it cannot be chosen
		if cY > k_coords[part,0]:
			distance = 100000
		#distance from body part to ball
		else:
			distance = (cX - k_coords[part, 1])**2 + (cY - k_coords[part,0])**2

		distances.append(distance)

	body_part = distances.index(min(distances))

	#if detection is too far away from body, or if all body parts are penalized, then the ball hit something else, like the ground
	if abs(cX - k_coords[parts[body_part], 1]) > knee_distance * 2 or sum(distances) == 5 * 100000:
		return 5
	# 0 = head, 1 = left thigh, 2 = right thigh, 3 = left foot, 4 = right foot, 5 = ground
	return body_part

def run_video(path, net, sess, output_stride, model_outputs):
	'''Runs the uploaded video through the juggle counter algorithm
	'''

	#initialize constants and counters
	cv2.setUseOptimized(True)
	(body_part_key, body_part_bounces, body_part_key_simple, body_part_sequence, confidence, W, H, bounces, frame_num, last_juggle_frame, first_detection, cYs, cXs, skipped_frames, ground) = initialize_constants()

	# start the video stream
	print("Processing video...")
	vs = FileVideoStream(path).start()

	while vs.more() and not ground:

		(frame, frame_num) = read_frame(vs, frame_num)

		#resize the frame to fit the 300x300 mobilenet size
		try:
			frame = cv2.resize(frame, (W, H))
		except cv2.error as e:
			print("...end of video")
			break

		# pass the frame through the model and output detection
		detection = process_frame(frame, net, W, H)

		# get the centroid coordinate and size of ball detected
		cX, cY, diameter_temp = get_centroid(frame, detection, W, H)

		# Initialize ball: Determine the ball size from the detection on the first frame (most likely to be clear/stationary ball)
		if first_detection:
			diameter = diameter_temp
			first_detection = False

		# If the detection is not the same size as the known ball, skip it. Also only allow detections within a realistic distance from the previous detection (two diameters of the ball * # of consecutive skipped frames)
		skip_frame = True

		if (0.7 * diameter < diameter_temp < 2 * diameter) and (((cXs[-1] - cX)**2 + (cYs[-1] - cY)**2)**0.5 < skipped_frames * 2 * diameter):
			skip_frame = False
			skipped_frames = 1

		#if we've determined that there is a valide detection, show the detection of ball and add centroid to lists
		if not skip_frame:
			cv2.circle(frame, (cX, cY), int(W * 0.03), (0, 255, 0), -1)
			cXs.append(cX)
			cYs.append(cY)

			#check if there was a juggle, and make sure it has been long enough since last juggle (this is a dirty fix for some false detections)
			if (frame_num - last_juggle_frame > 8) and check_bounce(cYs):

				last_juggle_frame = frame_num

				#pose detect (pose detect first output is a frame with the pose drawn on. If we don't want to show the pose)
				(_, k_scores, k_coords) = get_pose(frame, sess, output_stride, model_outputs)

				# determine which body part is closes to the ball
				body_part = get_closest_body_part(k_scores, k_coords, cXs[-2], cYs[-2])

				# add one juggle to the correct body part and overall score
				if body_part == 5:
					ground = True
				else:
					body_part_bounces[body_part_key[body_part]] += 1
					bounces += 1
				body_part_sequence.append(body_part)

		#if there is no valid detection in the frame, increase skipped_frames counter
		else:
			skipped_frames += 0.7

		#print scores on the frame
		frame = cv2.resize(frame, (600, 600))
		cv2.putText(frame, str(bounces), (int(2 * W * 0.76), int(W * 0.3)), cv2.FONT_HERSHEY_SIMPLEX, int(W * 0.01), (0, 255, 0), 2)

		# Determining right vs left is not accurate enough, so don't differentiate between them
		cv2.putText(frame, "Head: " + str(body_part_bounces["Head"]), (int(2 * W * 0.76), int(W * 0.5)), cv2.FONT_HERSHEY_SIMPLEX, int(W * 0.0035), (0, 255, 0), 2)
		cv2.putText(frame, "Thigh: " + str(body_part_bounces["Left Thigh"]+body_part_bounces["Right Thigh"]), (int(2 * W * 0.76), int(W * 0.65)), cv2.FONT_HERSHEY_SIMPLEX, int(W * 0.0035), (0, 255, 0), 2)
		cv2.putText(frame, "Foot:  " + str(body_part_bounces["Left Foot"]+body_part_bounces["Right Foot"]), (int(2 * W * 0.76), int(W * 0.8)), cv2.FONT_HERSHEY_SIMPLEX, int(W * 0.0035), (0, 255, 0), 2)

		#Print body part of previous juggle, if a juggle has been made
		try:
			cv2.putText(frame, body_part_key_simple[body_part_sequence[-1]], (int(2 * W * 0.05), int(W * 0.2)), cv2.FONT_HERSHEY_SIMPLEX, int(W * 0.007), (0, 255, 0), 2)
		except:
			pass

		# Show the frame
		ret, img = cv2.imencode(".jpg", frame)
		img = img.tobytes()
		#video_bytes.append(img)
		yield (b'--frame\r\n'
               b'Content-Type: image/jpeg\r\n\r\n' + img + b'\r\n')
		#yield (bounces, body_part_bounces, body_part_sequence)

		#if run in terminal instead of app, can have OpenCV show image
		#cv2.imshow("Frame", frame)
		#key = cv2.waitKey(1) & 0xFF

		# Only for debuggin:
		# if the `q` key was pressed, break from the loop
		# if key == ord("q"):
		# 	break
	cv2.destroyAllWindows()
	vs.stop()

	#eturn (video_bytes, bounces, body_part_bounces, body_part_sequence)
	return (b'--frame\r\n'
               b'Content-Type: image/jpeg\r\n\r\n' + img + b'\r\n')

# Optional function to display the frames from a list like a video in real time
# def display_video(video_bytes):
# 	'''Takes a list of image bytes and yields them consecutively to play in browser like a video'''
# 	for frame in video_bytes:
# 		time.sleep(.07)
# 		yield (b'--frame\r\n'
#                b'Content-Type: image/jpeg\r\n\r\n' + frame + b'\r\n')