text:
import pennylane as qml
from pennylane import numpy as np

dev = qml.device('default.qubit', wires = 5)

def basis_state(features, wires):
	qml.templates.state_preparations.basis.BasisStatePreparation(features, wires = range(4))
	# features(grad_required = False), they are never trained and are fixed.

@qml.qnode(dev)
def variational_layers():
	pass

def circuit(features, weights):
	# wires = qml.wires.Wires(wires)
	basis_state(features, wires = range(5))

	# variational_layers()

	return qml.expval(qml.PauliZ(wires = 0))

def variational_classifier(parameter, features):
	weights = parameter[0]
	bias = parameter[1]
	return circuit(features, weights) + bias

def square_loss(labels, predictions):
	loss = 0
	for l, p in zip(labels, predictions):
		loss += (l - p) ** 2

	loss = loss / len(labels)

	return loss

def accuracy():
	# Calculate how many inputs are predicted correctly
	accuracy = 0
	threshold = 1e-5

	for l, p in zip(labels, predictions):
		if (l - p) < threshold:
			accuracy += 1

	return accuracy

# def

features = np.array([0, 1, 1, 0])
parameters = [[1, 4, 5, 8], 0.0]

print(variational_classifier(parameters, features))

drawer = qml.draw(circuit)
print(drawer(features))

max_stars_repo_name: SaashaJoshi/qml

max_issues_repo_name: SaashaJoshi/quantum-machine-learning

max_forks_repo_name: SaashaJoshi/quantum-machine-learning


text:
#!/usr/bin/env python
__author__ = 'heroico'
import os
import io
import gzip
import logging
import numpy

from metax import WeightDBUtilities
from metax import PrediXcanFormatUtilities
from metax import ThousandGenomesUtilities
from metax import Logging
from metax import Utilities
from metax import Formats
from timeit import default_timer as timer

def pathLeaf(path):
    head, tail = os.path.split(path)
    return tail or os.path.basename(head)

class ProcessWeightDB(object):
    def __init__(self, args):
        self.weight_db = pathLeaf(args.weight_db)
        self.db_path = args.weight_db
        self.data_folder = args.input_folder
        self.correlation_output = args.correlation_output
        self.covariance_output = args.covariance_output
        if args.covariance_output is None:
            comp = os.path.splitext(self.weight_db)[0]
            name = comp + ".cov.txt.gz"
            path = os.path.join("intermediate", "cov")
            path = os.path.join(path, name)
            self.covariance_output = path

        self.input_format = args.input_format

        self.found_genes_for_covariance = {}
        self.found_genes_for_correlation = {}

        self.min_maf_filter = float(args.min_maf_filter) if args.min_maf_filter else None
        self.max_maf_filter = float(args.max_maf_filter) if args.max_maf_filter else None

        self.max_snps_in_gene = int(args.max_snps_in_gene) if args.max_snps_in_gene else None

        self.delimiter = args.delimiter

    def run(self):
        start = timer()

        if not self.correlation_output and not self.covariance_output:
            logging.info("Provide --correlation_output or --covariance_output or both")
            return

        logging.info("Loading Weights")
        weight_db_logic = WeightDBUtilities.WeightDBEntryLogic(self.db_path)

        logging.info("Building files")
        self.buildFiles(weight_db_logic)

        end = timer()
        logging.info("Ran successfully in %s seconds",  str(end-start))

    def buildFiles(self, weight_db_logic):
        do_correlations = self.correlation_output is not None
        if do_correlations:
            if os.path.exists(self.correlation_output):
                logging.info("%s already exists, delete it if you want it figured out again", self.correlation_output)
                do_correlations = False
            else:
                correlation_dir = os.path.dirname(self.correlation_output)
                if not os.path.exists(correlation_dir):
                    os.makedirs(correlation_dir)
                self.writeFileHeader(self.correlation_output)

        do_covariances = self.covariance_output is not None
        if do_covariances:
            if os.path.exists(self.covariance_output):
                logging.info("%s already exists, delete it if you want it figured out again", self.covariance_output)
                do_covariances = False
            else:
                covariance_dir = os.path.dirname(self.covariance_output)
                if not os.path.exists(covariance_dir):
                    os.makedirs(covariance_dir)
                self.writeFileHeader(self.covariance_output)

        if not do_covariances and not do_correlations:
            return

        names = Utilities.dosageNamesFromFolder(self.data_folder)
        for name in names:
            snps, snps_by_rsid = self.getSNPS(name, weight_db_logic)
            if do_correlations:
                self.addToCorrelationFile(weight_db_logic, name, snps, snps_by_rsid)

            if do_covariances:
                self.addToCovarianceFile(weight_db_logic, name, snps, snps_by_rsid)

    def writeFileHeader(self,path):
        with io.TextIOWrapper(gzip.open(path, "ab"), newline="") as file:
            file.write("GENE RSID1 RSID2 VALUE\n")

    def getSNPS(self, name, weight_db_logic):
        dosageLoader = None
        if self.input_format == Formats.IMPUTE:
            #TODO: update
            dosageLoader = ThousandGenomesUtilities.IMPUTEDosageLoader(self.data_folder, name) #outdated code
        elif self.input_format == Formats.PrediXcan:
            #dosageName = Utilities.dosageName(name)
            path = os.path.join(self.data_folder, name)
            dosageLoader = PrediXcanFormatUtilities.PrediXcanFormatDosageLoader(path, weight_db_logic, self.delimiter)
        else:
            logging.info("Invalid input format: %s", self.input_format)
            return
        snps, snps_by_rsid = dosageLoader.load()
        return snps, snps_by_rsid

    def addToCovarianceFile(self, weight_db_logic, name, snps, snps_by_rsid):
        logging.info("Adding to covariance for %s-%s", name, self.weight_db)

        genes = weight_db_logic.weights_by_gene.keys()
        total_genes = len(genes)
        last_reported_percent = 0
        processed = 0
        for gene in genes:
            processed += 1
            percent = int(processed*100.0 / total_genes)
            if percent == last_reported_percent+10:
                logging.info("%d percent genes processed", percent)
                last_reported_percent = percent

            entries = self.buildCovarianceEntries(name, gene, weight_db_logic, snps_by_rsid)

            if len(entries) == 0:
                logging.log(6,"Gene %s has no snps in current file", gene)
                continue

            self.addToFile(self.covariance_output, gene, entries)

    def addToFile(self, path, gene, entries):
        with io.TextIOWrapper(gzip.open(path, "ab"), newline="") as file:
            for entry in entries:
                line = " ".join([gene, entry[0], entry[1], entry[2]])+"\n"
                file.write(line)

    def buildCovarianceEntries(self, name, gene, weight_db_logic, snps_by_rsid):
        weights_in_gene = weight_db_logic.weights_by_gene[gene]
        rsids_from_genes = list(weights_in_gene.keys())

        #gather as much data as we can work on
        related_rsids, related_data = self.buildRelatedData(rsids_from_genes, snps_by_rsid, weights_in_gene)

        if len(related_rsids) == 0:
            return []

        self.updateFoundCovariance(gene, name)

        #covariance matrix of related SNP's data
        array = numpy.array(related_data)
        cov = numpy.cov(array)

        #translate into sql entries
        entries = self.buildMatrixOutputEntries(cov, rsids_from_genes, related_rsids, snps_by_rsid)
        if not len(entries):
            raise NameError("Couldn not build covariance entries for (%s,%s)" %(name,gene))
        return entries

    def updateFoundCovariance(self, gene, name):
        found = None
        if gene in self.found_genes_for_covariance:
            found = self.found_genes_for_covariance[gene]
            logging.info("Gene %s found again for %s", gene, name)
        else:
            found = []
            self.found_genes_for_covariance[gene] = found
        found.append(name)

    def buildRelatedData(self, rsids_from_genes, snps_by_rsid, weights_in_gene):
        related_rsids = []
        related_data = []

        l = len(rsids_from_genes)
        if self.max_snps_in_gene and l > self.max_snps_in_gene:
            logging.info("Skipping covariance too large: %d", l)
            return related_data, related_rsids

        for rsid in rsids_from_genes:
            if not rsid in snps_by_rsid:
                logging.log(5, "related rsid %s not present in genotype data", rsid)
                continue

            related_snp = snps_by_rsid[rsid]
            freq = sum(related_snp.data)*1.0/(2*len(related_snp.data))
            if self.min_maf_filter and self.min_maf_filter > freq:
                logging.log(6, "related rsid %s below min maf: %s", rsid, freq)
                continue

            if self.max_maf_filter and self.max_maf_filter < freq:
                logging.log(6, "related rsid %s  above max maf: %s", rsid, freq)
                continue

            data = related_snp.data
            weight = weights_in_gene[rsid]
            if weight.ref_allele == related_snp.eff_allele and\
                weight.eff_allele == related_snp.ref_allele:
                logging.log(7, "related rsid %s has alleles flipped compared to model, transforming dosage", rsid)
                data = [2-x for x in data]

            related_data.append(data)
            related_rsids.append(rsid)
        return related_rsids, related_data

    def buildMatrixOutputEntries(self, matrix, rsids_from_genes, related_rsids, snps_by_rsid):
        entries = []

        #special case: we might have a single rsid!
        if matrix.ndim == 0:
            c = str(float(matrix))
            r = rsids_from_genes[0]
            entries.append((r,r,c))
            return entries

        for i in range(0, len(rsids_from_genes)):
            rsid_i = rsids_from_genes[i]
            related_i = -1
            if rsid_i in related_rsids:
                related_i = related_rsids.index(rsid_i)

            for j in range(0, len(rsids_from_genes)):
                rsid_j = rsids_from_genes[j]

                related_j = -1
                if rsid_j in related_rsids:
                    related_j = related_rsids.index(rsid_j)

                value = "NA"
                if related_i > -1 and related_j > -1:
                    value = str(matrix[related_i][related_j])

                if i == j:
                    entries.append((rsid_i, rsid_i, value))
                else:
                    if value == "NA":
                        if rsid_i < rsid_j:
                            entries.append((rsid_i, rsid_j, value))
                    else:
                        snp_i = snps_by_rsid[rsid_i]
                        snp_j = snps_by_rsid[rsid_j]

                        if snp_i.position < snp_j.position:
                            entries.append((rsid_i, rsid_j, value))
        return entries

    def addToCorrelationFile(self, weight_db_logic, name, snps, snps_by_rsid):
        logging.info("Building correlation database for %s-%s", name, self.weight_db)
        genes = list(weight_db_logic.weights_by_gene.keys())
        total_genes = len(genes)
        last_reported_percent = 0
        processed = 0
        for gene in genes:
            processed += 1
            percent = int(processed*100.0 / total_genes)
            if percent == last_reported_percent+10:
                logging.info("%d percent genes processed", percent)
                last_reported_percent = percent

            entries = self.buildCorrelationEntries(name, gene, weight_db_logic, snps_by_rsid)

            if len(entries) == 0:
                logging.log(6,"Gene %s has no snps in current file", gene)
                continue

            self.addToFile(self.correlation_output, gene, entries)

    def buildCorrelationEntries(self, name, gene, weight_db_logic, snps_by_rsid):
        weights_in_gene = weight_db_logic.weights_by_gene[gene]
        rsids_from_genes = list(weights_in_gene.keys())

        #gather as much data as we can work on
        related_rsids, related_data = self.buildRelatedData(rsids_from_genes, snps_by_rsid, weights_in_gene)

        if len(related_rsids) == 0:
            return []

        self.updateFoundCorrelation(gene, name)

        #correlation matrix of related SNP's data
        array = numpy.array(related_data)
        cor = numpy.corrcoef(array)

        #translate into sql entries
        entries = self.buildMatrixOutputEntries(cor, rsids_from_genes, related_rsids, snps_by_rsid)
        if not len(entries):
            raise NameError("Couldn not build correlation entries for (%s,%s)" %(name,gene))
        return entries

    def updateFoundCorrelation(self, gene, name):
        found = None
        if gene in self.found_genes_for_correlation:
            found = self.found_genes_for_correlation[gene]
            logging.info("Gene %s found again for %s", gene, name)
        else:
            found = []
            self.found_genes_for_correlation[gene] = found
        found.append(name)

if __name__ == "__main__":
    import argparse
    parser = argparse.ArgumentParser(description='Build correlations and/or covariances from PHASE3 data and weights database.')

    parser.add_argument("--verbosity",
                        help="Log verbosity level. 1 is everything being logged. 10 is only high level messages, above 10 will hardly log anything",
                        default = "10")

    parser.add_argument("--weight_db",
                        help="name of weight db in data folder",
                        default="data/DGN-WB_0.5.db")

    parser.add_argument("--input_folder",
                        help="name of folder containing PHASE 3 data",
                        default="intermediate/TGF_EUR")

    parser.add_argument("--correlation_output",
                        help="Name of file to dump correlation results in.",
                        default=None)
                        #default="intermediate/1000GP_Phase3_chr_cor_dgnwb_cor.txt.gz")

    parser.add_argument("--covariance_output",
                        help="Name of file to dump covariance results in. Defaults to 'intermediate/cov/' + file name prefix from '--weight_db' argument",
                        default=None)

    parser.add_argument('--input_format',
                   help='Input dosage files format. Valid options are: IMPUTE, PrediXcan',
                   default=Formats.PrediXcan)

    parser.add_argument('--min_maf_filter',
                   help="Filter snps according to this maf",
                   default=None)

    parser.add_argument('--max_maf_filter',
                   help="Filter snps according to this maf",
                   default=None)

    parser.add_argument("--max_snps_in_gene",
                        help="Ignore any gene that has snps above this value",
                        type=int,
                        default=None)

    parser.add_argument("--delimiter", help="dosage delimiter", default=" ")

    args = parser.parse_args()

    Logging.configureLogging(int(args.verbosity))

    work = ProcessWeightDB(args)
    work.run()

ext: py

max_stars_repo_name: adellanno/MetaXcan

max_issues_repo_name: adellanno/MetaXcan

max_forks_repo_name: adellanno/MetaXcan

reason: import numpy


text:
import sys
import os
sys.path.append(os.pardir)

import numpy as np
import unittest

import coropy.utils as cu

class TestUtils(unittest.TestCase):
    """Class for `utils.py` tests."""

    def test_normalize(self):
        self.assertIsNone(
            np.testing.assert_array_equal(cu.normalize([10, 50]), [0., 1.]))
        self.assertIsNone(
            np.testing.assert_array_almost_equal(
                cu.normalize(range(100)), np.linspace(0, 1, 100), decimal=2))
        self.assertIsNone(
            np.testing.assert_array_equal(cu.normalize([-1, 1]), [0., 1.]))
        self.assertRaises(ValueError, cu.normalize, 5.)
        self.assertRaises(AssertionError, cu.normalize, [[1, 2, 3]])
        self.assertRaises(AssertionError, cu.normalize, [5])
        self.assertRaises(ValueError, cu.normalize, [-np.nan, 50])
        self.assertWarns(RuntimeWarning, cu.normalize, [0, 0, 0])

    def test_restore(self):
        self.assertIsNone(
            np.testing.assert_array_equal(
                cu.restore([0., 1.], [10, 50]), [10., 50.]))
        self.assertIsNone(
            np.testing.assert_array_almost_equal(
                cu.restore(np.linspace(0, 1, 100), range(100)),
                range(100),
                decimal=2))
        self.assertRaises(ValueError, cu.restore, [5.], 0)
        self.assertRaises(AssertionError, cu.restore, [[1, 2, 3]], [0, 3])
        self.assertRaises(ValueError, cu.restore, [0, 1], [-np.nan, 50])

    def test_moving_average(self):
        a = np.array([1, 2, 3, 4, 5])
        a_2 = np.array([1.5, 2.5, 3.5, 4.5])
        a_5 = np.array([3.])
        self.assertIsNone(
            np.testing.assert_array_equal(cu.moving_average(a, 1), a))
        self.assertIsNone(
            np.testing.assert_array_equal(cu.moving_average(a, 2), a_2))
        self.assertIsNone(
            np.testing.assert_array_equal(cu.moving_average(a, 5), a_5))
        self.assertRaises(ValueError, cu.moving_average, a, 1.5)
        self.assertRaises(ValueError, cu.moving_average, 1)
        self.assertRaises(AssertionError, cu.moving_average, a.reshape(-1, 1))
        self.assertRaises(AssertionError, cu.moving_average, a, 6)

    def test_train_test_split(self):
        a_list = [range(10)]
        a = np.array(a_list)
        a_train, a_test = a[:8], a[8:]
        self.assertIsNone(
            np.testing.assert_array_equal(
                cu.train_test_split(a, 0.8)[0], a_train))
        self.assertIsNone(
            np.testing.assert_array_equal(
                cu.train_test_split(a, 0.8)[1], a_test))
        self.assertRaises(ValueError, cu.train_test_split, a_list, .5)
        self.assertRaises(ValueError, cu.train_test_split, a, 1.5)

    def test_mse(self):
        y_true_list = [0, 1, 2, 3, 4, 5]
        y_true = np.array(y_true_list)
        y_pred_list = [1, 2, 3, 4, 5, 6]
        y_pred = np.array(y_pred_list)
        mse = 1.0
        self.assertEqual(cu.mse(y_true, y_pred), mse)
        self.assertRaises(ValueError, cu.mse, y_true_list, y_pred_list)

    def test_rmse(self):
        y_true_list = [0, 1, 2, 3, 4, 5]
        y_true = np.array(y_true_list)
        y_pred_list = [1, 2, 3, 4, 5, 6]
        y_pred = np.array(y_pred_list)
        rmse = 1.0
        self.assertEqual(cu.rmse(y_true, y_pred), rmse)
        self.assertRaises(ValueError, cu.rmse, y_true_list, y_pred_list)

    def test_msle(self):
        y_true_list = [0, 1, 2, 3, 4, 5]
        y_true = np.array(y_true_list)
        y_pred_list = [1, 2, 3, 4, 5, 6]
        y_pred = np.array(y_pred_list)
        msle = 0.13906
        self.assertAlmostEqual(round(cu.msle(y_true, y_pred), 5), msle, places=4)
        self.assertRaises(ValueError, cu.msle, y_true_list, y_pred_list)

    def test_mae(self):
        y_true_list = [0, 1, 2, 3, 4, 5]
        y_true = np.array(y_true_list)
        y_pred_list = [1, 2, 3, 4, 5, 6]
        y_pred = np.array(y_pred_list)
        mae = 1.0
        self.assertEqual(cu.mse(y_true, y_pred), mae)
        self.assertRaises(ValueError, cu.mae, y_true_list, y_pred_list)

if __name__ == "__main__":
    unittest.main()

max_issues_repo_name: antelk/coropy

max_forks_repo_name: antelk/coropy


text:
import unittest

import numpy as np
import tensorflow as tf
from mock import Mock

from tfsnippet.dataflow import DataFlow
from tfsnippet.scaffold import TrainLoop
from tfsnippet.trainer import *

class AutoLossWeightTestCase(unittest.TestCase):

    def test_auto_loss_weight(self):
        self.assertEquals(5., auto_loss_weight(np.arange(5)))
        self.assertEquals(7., auto_loss_weight(np.arange(7), np.arange(5)))
        self.assertEquals(1., auto_loss_weight(None))

class ValidatorTestCase(tf.test.TestCase):

    def test_props(self):
        loop = Mock(valid_metric_name='valid_loss')
        df = Mock()

        v = Validator(loop, 12, [34, 56], df)
        self.assertIs(loop, v.loop)
        self.assertEquals(12, v.loss)
        self.assertEquals([34, 56], v.inputs)
        self.assertIs(v.data_flow, df)
        self.assertEquals({}, v.feed_dict)
        self.assertEquals('valid_time', v.time_metric_name)
        self.assertEquals('valid_loss', v.loss_metric_name)
        self.assertIs(auto_loss_weight, v.loss_weight_func)

        loss_weight_func = Mock(return_value=123.)
        v = Validator(loop, 12, [34, 56], df,
                      feed_dict={'a': 1},
                      time_metric_name='valid_time_x',
                      loss_metric_name='valid_loss_x',
                      loss_weight_func=loss_weight_func)
        self.assertEquals('valid_time_x', v.time_metric_name)
        self.assertEquals('valid_loss_x', v.loss_metric_name)
        self.assertIs(loss_weight_func, v.loss_weight_func)

    def test_run(self):
        with self.test_session() as session:
            df = DataFlow.arrays([np.arange(6, dtype=np.float32)], batch_size=4)
            ph = tf.placeholder(tf.float32, shape=[None])
            ph2 = tf.placeholder(tf.float32, shape=[])
            ph3 = tf.placeholder(tf.float32, shape=[])

            # test default loss weight and merged feed dict
            with TrainLoop([], max_epoch=1) as loop:
                v = Validator(loop, tf.reduce_mean(ph), [ph], df,
                              feed_dict={ph2: 34})
                v._run_batch = Mock(wraps=v._run_batch)

                for epoch in loop.iter_epochs():
                    v.run({ph3: 56})
                    np.testing.assert_almost_equal(
                        2.5, loop._epoch_metrics._metrics['valid_loss'].mean)

                self.assertEquals(2, len(v._run_batch.call_args_list))
                for i, call_args in enumerate(v._run_batch.call_args_list):
                    call_session, call_feed_dict = call_args[0]
                    self.assertIs(session, call_session)
                    np.testing.assert_equal(
                        np.arange(6, dtype=np.float32)[i*4: (i+1)*4],
                        call_feed_dict[ph]
                    )
                    self.assertEquals(34, call_feed_dict[ph2])
                    self.assertEquals(56, call_feed_dict[ph3])

            # test None loss weight and override feed dict
            with TrainLoop([], max_epoch=1) as loop:
                v = Validator(loop, tf.reduce_mean(ph), [ph], df,
                              feed_dict={ph2: 34},
                              loss_weight_func=None,
                              loss_metric_name='valid_loss_x')
                v._run_batch = Mock(wraps=v._run_batch)

                for epoch in loop.iter_epochs():
                    v.run({ph2: 56})
                    np.testing.assert_almost_equal(
                        3.0, loop._epoch_metrics._metrics['valid_loss_x'].mean)

                for i, call_args in enumerate(v._run_batch.call_args_list):
                    call_session, call_feed_dict = call_args[0]
                    self.assertEquals(56, call_feed_dict[ph2])
                    self.assertNotIn(ph3, call_feed_dict)

if __name__ == '__main__':
    unittest.main()


text:
# Importing NumPy Library
import numpy as np
import sys

# Reading number of unknowns
n = int(input('Enter number of unknowns: '))

# Making numpy array of n x n+1 size and initializing
# to zero for storing augmented matrix
a = np.zeros((n,n+1))

# Making numpy array of n size and initializing
# to zero for storing solution vector
x = np.zeros(n)

# Reading augmented matrix coefficients
print('Enter Augmented Matrix Coefficients:')
for i in range(n):
    for j in range(n+1):
        a[i][j] = float(input( 'a['+str(i)+']['+ str(j)+']='))

# Applying Gauss Elimination
for i in range(n):
    if a[i][i] == 0.0:
        sys.exit('Divide by zero detected!')

    for j in range(i+1, n):
        ratio = a[j][i]/a[i][i]

        for k in range(n+1):
            a[j][k] = a[j][k] - ratio * a[i][k]

# Back Substitution
x[n-1] = a[n-1][n]/a[n-1][n-1]

for i in range(n-2,-1,-1):
    x[i] = a[i][n]

    for j in range(i+1,n):
        x[i] = x[i] - a[i][j]*x[j]

    x[i] = x[i]/a[i][i]

# Displaying solution
print('\nRequired solution is: ')
for i in range(n):
    print('X%d = %0.2f' %(i,x[i]), end = '\t')

max_stars_repo_name: nikitazviahin/UniversityProjects

max_issues_repo_name: nikitazviahin/UniversityProjects

max_forks_repo_name: nikitazviahin/UniversityProjects


text:
import numpy as np
import cv2

# Identify pixels above the threshold
# Threshold of RGB > 160 does a nice job of identifying ground pixels only
def above_color_thresh(img, rgb_thresh=(160, 160, 160)):
    # Create an array of zeros same xy size as img, but single channel
    color_select = np.zeros_like(img[:,:,0])
    # Require that each pixel be above all three threshold values in RGB
    # above_thresh will now contain a boolean array with "True"
    # where threshold was met
    above_thresh = (img[:,:,0] > rgb_thresh[0]) \
                & (img[:,:,1] > rgb_thresh[1]) \
                & (img[:,:,2] > rgb_thresh[2])
    # Index the array of zeros with the boolean array and set to 1
    color_select[above_thresh] = 1
    # Return the binary image
    return color_select

# Identify pixels below threshold
def below_color_thresh(img, rgb_thresh=(110, 110, 110)):
    # Create an array of zeros same xy size as img, but single channel
    color_select = np.zeros_like(img[:,:,0])
    # Require that each pixel be BELOW all three threshold values in RGB
    # below_thresh will now contain a boolean array with "True"
    # where threshold was met
    below_thresh = (img[:,:,0] < rgb_thresh[0]) \
                   & (img[:,:,1] < rgb_thresh[1]) \
                   & (img[:,:,2] < rgb_thresh[2])
    # Index the array of zeros with the boolean array and set to 1
    color_select[below_thresh] = 1
    # Return the binary image

    return color_select

# Identify pixels between two threshold values
def between_color_thresh(img, lower_rgb_thresh=(130, 120, 30), upper_rgb_thresh=(220,180,90)):
    # Create an array of zeros same xy size as img, but single channel
    color_select = np.zeros_like(img[:,:,0])
    # Require that each pixel be BETWEEN all three lower threshold values in RGB
    # and all three upper threshold values in RGB
    # below_thresh will now contain a boolean array with "True"
    # where both conditions were met (sample was above lower threshold and below upper threshold)

    between_thresh = ((img[:,:,0] > lower_rgb_thresh[0]) == ((img[:,:,0] < upper_rgb_thresh[0]))) \
                     & ((img[:,:,1] > lower_rgb_thresh[1]) == (img[:,:,1] < upper_rgb_thresh[1])) \
                     & ((img[:,:,2] > lower_rgb_thresh[2]) == (img[:,:,2] < upper_rgb_thresh[2]))
    # Index the array of zeros with the boolean array and set to 1
    color_select[between_thresh] = 1

    # Return the binary image
    return color_select

# Define a function to convert from image coords to rover coords
def rover_coords(binary_img):
    # Identify nonzero pixels
    ypos, xpos = binary_img.nonzero()
    # Calculate pixel positions with reference to the rover position being at the
    # center bottom of the image.
    x_pixel = -(ypos - binary_img.shape[0]).astype(np.float)
    y_pixel = -(xpos - binary_img.shape[1]/2 ).astype(np.float)
    return x_pixel, y_pixel

# Define a function to convert to radial coords in rover space
def to_polar_coords(x_pixel, y_pixel):
    # Convert (x_pixel, y_pixel) to (distance, angle)
    # in polar coordinates in rover space
    # Calculate distance to each pixel
    dist = np.sqrt(x_pixel**2 + y_pixel**2)
    # Calculate angle away from vertical for each pixel
    angles = np.arctan2(y_pixel, x_pixel)
    return dist, angles

# Define a function to map rover space pixels to world space
def rotate_pix(xpix, ypix, yaw):
    # Convert yaw to radians
    yaw_rad = yaw * np.pi / 180
    xpix_rotated = (xpix * np.cos(yaw_rad)) - (ypix * np.sin(yaw_rad))

    ypix_rotated = (xpix * np.sin(yaw_rad)) + (ypix * np.cos(yaw_rad))
    # Return the result
    return xpix_rotated, ypix_rotated

def translate_pix(xpix_rot, ypix_rot, xpos, ypos, scale):
    # Apply a scaling and a translation
    xpix_translated = (xpix_rot / scale) + xpos
    ypix_translated = (ypix_rot / scale) + ypos
    # Return the result
    return xpix_translated, ypix_translated

# Define a function to apply rotation and translation (and clipping)
# Once you define the two functions above this function should work
def pix_to_world(xpix, ypix, xpos, ypos, yaw, world_size, scale):
    # Apply rotation
    xpix_rot, ypix_rot = rotate_pix(xpix, ypix, yaw)
    # Apply translation
    xpix_tran, ypix_tran = translate_pix(xpix_rot, ypix_rot, xpos, ypos, scale)
    # Perform rotation, translation and clipping all at once
    x_pix_world = np.clip(np.int_(xpix_tran), 0, world_size - 1)
    y_pix_world = np.clip(np.int_(ypix_tran), 0, world_size - 1)
    # Return the result
    return x_pix_world, y_pix_world

# Define a function to perform a perspective transform
def perspect_transform(img, src, dst):

    M = cv2.getPerspectiveTransform(src, dst)
    warped = cv2.warpPerspective(img, M, (img.shape[1], img.shape[0]))# keep same size as input image

    return warped

# Apply the above functions in succession and update the Rover state accordingly
def perception_step(Rover):
    # Perform perception steps to update Rover()
    # NOTE: camera image is coming to you in Rover.img
    # 1) Define source and destination points for perspective transform

    # Define calibration box in source (actual) and destination (desired) coordinates
    # These source and destination points are defined to warp the image
    # to a grid where each 10x10 pixel square represents 1 square meter
    # The destination box will be 2*dst_size on each side
    dst_size = 5
    # Set a bottom offset to account for the fact that the bottom of the image
    # is not the position of the rover but a bit in front of it
    bottom_offset = 6
    source = np.float32([[14, 140], [301 ,140],[200, 96], [118, 96]])
    destination = np.float32([[Rover.img.shape[1]/2 - dst_size, Rover.img.shape[0] - bottom_offset],
                              [Rover.img.shape[1]/2 + dst_size, Rover.img.shape[0] - bottom_offset],
                              [Rover.img.shape[1]/2 + dst_size, Rover.img.shape[0] - 2*dst_size - bottom_offset],
                              [Rover.img.shape[1]/2 - dst_size, Rover.img.shape[0] - 2*dst_size - bottom_offset],
                              ])

    # 2) Apply perspective transform
    warped = perspect_transform(Rover.img, source, destination)

    # 3) Apply color threshold to identify navigable terrain/obstacles/rock samples
    navigable = above_color_thresh(warped, (160, 160, 160))
    obstacle = below_color_thresh(warped, (100, 100, 100))
    rock = between_color_thresh(warped, (140, 155, 50), (190, 190, 130))

    # 4) Update Rover.vision_image (this will be displayed on left side of screen)
    # Example: Rover.vision_image[:,:,0] = obstacle color-thresholded binary image
    #          Rover.vision_image[:,:,1] = rock_sample color-thresholded binary image
    #          Rover.vision_image[:,:,2] = navigable terrain color-thresholded binary image
    # For now, do not update vision_image for rock. That will be handled later.
    Rover.vision_image[:, :, 0] = obstacle * 255
    Rover.vision_image[:, :, 2] = navigable * 255

    # 5) Convert map image pixel values to rover-centric coords
    nav_xpix, nav_ypix = rover_coords(navigable)
    obs_xpix, obs_ypix = rover_coords(obstacle)
    rock_xpix, rock_ypix = rover_coords(rock)

    # 6) Convert rover-centric pixel values to world coordinates
    # First set the scale to 2 x the destination size.
    scale = 2 * dst_size

    # Extract coordinates and yaw from telemetry provided
    rover_xpos = Rover.pos[0]
    rover_ypos = Rover.pos[1]
    rover_yaw = Rover.yaw

    # Now use telemetry data to get navigable pixel positions in world coords
    # Handle navigable pixels
    nav_x_world, nav_y_world = pix_to_world(nav_xpix, nav_ypix, rover_xpos,
                                            rover_ypos, rover_yaw,
                                            Rover.worldmap.shape[0], scale)

    # Now convert obstacle pixels
    obs_x_world, obs_y_world = pix_to_world(obs_xpix, obs_ypix, rover_xpos,
                                            rover_ypos, rover_yaw,
                                            Rover.worldmap.shape[0], scale)

    # And finally convert rock pixels
    rock_x_world, rock_y_world = pix_to_world(rock_xpix, rock_ypix, rover_xpos,
                                              rover_ypos, rover_yaw,
                                              Rover.worldmap.shape[0], scale)

    # For rock pixels, the perspective-transformed version of the camera image will
    # have what amounts to a "smear" where the rock was found. Rather than map all
    # of those pixels, which would be wrong, map only the closest pixel (which essentially
    # represents the face of the rock).
    if rock.any(): # Checks the image to see if there are even any rocks in the image
        # Found rocks, so convert to polar coordinates
        rock_dist, rock_angles = to_polar_coords(rock_xpix, rock_ypix)
        # Now get the index of the nearest pixel - this will basically be the face of the rock
        # and the true position of it. These steps ignore the smear that would be seen on the map
        rock_idx = np.argmin(rock_dist)
        rock_xcen = rock_x_world[rock_idx]
        rock_ycen = rock_y_world[rock_idx]
        # Now that we have identified the rock location in the code just above, set the vision_image
        # to show the whole rock map.
        Rover.vision_image[:, :, 1] = rock * 255
    else: # No rocks found so zero out the channel
        Rover.vision_image[:, :, 1] = 0

    # 7) Update Rover worldmap (to be displayed on right side of screen)
        # Example: Rover.worldmap[obstacle_y_world, obstacle_x_world, 0] += 1
        #          Rover.worldmap[rock_y_world, rock_x_world, 1] += 1
        #          Rover.worldmap[navigable_y_world, navigable_x_world, 2] += 1

    # Check roll and pitch and only add to map if these values are very close to zero, otherwise
    # accuracy of the map is compromised
    if (Rover.roll > 359.5 or Rover.roll < .5) and (Rover.pitch > 359.5 or Rover.pitch < .5):
        # if any rock pixels are present, add the xcen,ycen coordinates to the worldmap. By
        # adding only the xcen,ycen coordinates, this effectively maps only the rock location
        # and not the smear.
        if rock.any():
            Rover.worldmap[rock_ycen, rock_xcen, 1] = 255

        # Now add the nav and obstacle pixels to the worldmap
        Rover.worldmap[nav_y_world, nav_x_world, 2] = 255
        Rover.worldmap[obs_y_world, obs_x_world, 0] = 255

        # There will be overlap between obstacle pixels and navigable pixels. Assume that a navigable
        # pixel "over-rides" an obstacle pixel. In that case, set the corresponding obstacle pixel to
        # 0.
        # This next line extracts the navigable pixels by examining the appropriate channel of the
        # worldmap image. nav_pix will contain essentially a copy of that channel containing all
        # the nav pix
        nav_pix = Rover.worldmap[:,:,2] > 0
        # Now use the nav pixels as a "masK" to set to 0 the corresponding pixels in the obstacle channel
        # The following iterates through the nav_pix array, which is an array of x,y coordinates in
        # world view that have nav pixels. If there is a nav pixel present, then set the corresponding
        # pixel in the obstacle map (channel 0) to 0. This is where the obstacle map is being over-ridden
        # with the navigable pixel data
        Rover.worldmap[nav_pix, 0] = 0

    # 8) Convert rover-centric pixel positions to polar coordinates
    nav_dist, nav_angles = to_polar_coords(nav_xpix, nav_ypix)

    # Update Rover pixel distances and angles
    Rover.nav_dists = np.mean(nav_angles)
    Rover.nav_angles = nav_angles

    return Rover


text:
"""Runge-Kutta methods
"""

import numpy as np

class RungeKutta(object):
    """Base class for Runge-Kutta methods.
    """
    nstage: int = NotImplemented
    butcherTable: np.ndarray = NotImplemented

    def __init__(self):
        self.beta = self.butcherTable[:-1, 1:]
        self.w = self.butcherTable[-1, 1:]
        self.alpha = self.butcherTable[:-1,0]

    def step(self):
        """Perform one integration step.
        """
        raise NotImplementedError

    def nstep(self, fun, u, t0=0, dt=1., n=1):
        """Perform n times integration steps.

        Args:
            fun (callable): Right-hand side of the equation. du/dt = fun(t, u).
            u (np.array): Initial state.
            t0 (float): Initial time.
            dt (float): Timestep size.
            n (int, optional): Number of steps. Defaults to 1.

        Returns:
            np.array: u(t0 + dt * n)
        """
        self.t = t0
        for i in range(1, n+1):
            u = self.step(fun, u, t0, dt)
            self.t = t0 + dt*i
        return u

class RungeKutta2(RungeKutta):
    """Base class for Runge-Kutta methods of order 2.
    """
    def step(self, fun, u, t0=0, dt=1.):
        """Perform one integration step.

        Args:
            fun (callable): Right-hand side of the system.
            u (np.array): Initial state.
            t0 (float): Initial time.
            dt (float): Timestep size.

        Returns:
            np.array: u(t + dt)
        """
        k0 = dt * fun(t0, u)
        k1 = dt * fun(t0 + dt*self.alpha[1], u + self.beta[1,0]*k0)
        return u + self.w[0]*k0 + self.w[1]*k1

class RungeKutta3(RungeKutta):
    """Base class for Runge-Kutta methods of order 3.
    """
    def step(self, fun, u, t0=0, dt=1.):
        """Perform one integration step.

        Args:
            fun (callable): Right-hand side of the system.
            u (np.array): Initial state.
            t0 (float): Initial time.
            dt (float): Timestep size.

        Returns:
            np.array: u(t + dt)
        """
        k0 = dt * fun(t0, u)
        k1 = dt * fun(t0 + dt*self.alpha[1], u + self.beta[1,0]*k0)
        k2 = dt * fun(t0 + dt*self.alpha[2], u + self.beta[2,0]*k0 + self.beta[2,1]*k1)
        return u + self.w[0]*k0 + self.w[1]*k1 + self.w[2]*k2

class RK2(RungeKutta2):
    """Runge-Kutta methods of order 2.

    The Butcher Table is from Table I Case no.2 [1]_.

    References:
    .. [1] J.H Williamson, "Low-storage Runge-Kutta schemes", 1980.
    """
    butcherTable = np.array([[0, 0, 0],
                             [2/3, 2/3, 0],
                             [0, 1/4, 3/4]])

class RK3(RungeKutta3):
    """Classic Runge-Kutta methods of order 3.

    The Butcher Table is from Table I Case no.11 [1]_.

    References:
    .. [1] J.H Williamson, "Low-storage Runge-Kutta schemes", 1980.
    """
    butcherTable = np.array([[0, 0, 0, 0],
                             [1/2, 1/2, 0, 0],
                             [1, -1, 2, 0],
                             [0, 1/6, 2/3, 1/6]])

max_stars_repo_name: mkihara/cfdpy

max_issues_repo_name: mkihara/cfdpy

max_forks_repo_name: mkihara/cfdpy


text:
from typing import Union, Optional, Any, Mapping, Callable, NamedTuple, Generator, Tuple

import numpy as np
import scipy
from anndata import AnnData
from numpy.random import RandomState
from scipy.sparse import issparse, coo_matrix, csr_matrix
from sklearn.metrics import pairwise_distances
from sklearn.utils import check_random_state

from .. import logging as logg
from .. import utils
from ..utils import doc_params
from ..tools._utils import choose_representation, doc_use_rep, doc_n_pcs

N_DCS = 15  # default number of diffusion components
N_PCS = 50  # default number of PCs

Metric = Callable[[np.ndarray, np.ndarray], float]

@doc_params(n_pcs=doc_n_pcs, use_rep=doc_use_rep)
def neighbors(
    adata: AnnData,
    n_neighbors: int = 15,
    n_pcs: Optional[int] = None,
    use_rep: Optional[str] = None,
    knn: bool = True,
    random_state: Optional[Union[int, RandomState]] = 0,
    method: str = 'umap',
    metric: Union[str, Metric] = 'euclidean',
    metric_kwds: Mapping[str, Any] = {},
    copy: bool = False
) -> Optional[AnnData]:
    """\
    Compute a neighborhood graph of observations [McInnes18]_.

    The neighbor search efficiency of this heavily relies on UMAP [McInnes18]_,
    which also provides a method for estimating connectivities of data points -
    the connectivity of the manifold (`method=='umap'`). If `method=='gauss'`,
    connectivities are computed according to [Coifman05]_, in the adaption of
    [Haghverdi16]_.

    Parameters
    ----------
    adata
        Annotated data matrix.
    n_neighbors
        The size of local neighborhood (in terms of number of neighboring data
        points) used for manifold approximation. Larger values result in more
        global views of the manifold, while smaller values result in more local
        data being preserved. In general values should be in the range 2 to 100.
        If `knn` is `True`, number of nearest neighbors to be searched. If `knn`
        is `False`, a Gaussian kernel width is set to the distance of the
        `n_neighbors` neighbor.
    {n_pcs}
    {use_rep}
    knn
        If `True`, use a hard threshold to restrict the number of neighbors to
        `n_neighbors`, that is, consider a knn graph. Otherwise, use a Gaussian
        Kernel to assign low weights to neighbors more distant than the
        `n_neighbors` nearest neighbor.
    random_state
        A numpy random seed.
    method : {{'umap', 'gauss', `None`}}  (default: `'umap'`)
        Use 'umap' [McInnes18]_ or 'gauss' (Gauss kernel following [Coifman05]_
        with adaptive width [Haghverdi16]_) for computing connectivities.
    metric
        A known metric’s name or a callable that returns a distance.
    metric_kwds
        Options for the metric.
    copy
        Return a copy instead of writing to adata.

    Returns
    -------
    Depending on `copy`, updates or returns `adata` with the following:

    **connectivities** : sparse matrix (`.uns['neighbors']`, dtype `float32`)
        Weighted adjacency matrix of the neighborhood graph of data
        points. Weights should be interpreted as connectivities.
    **distances** : sparse matrix (`.uns['neighbors']`, dtype `float32`)
        Instead of decaying weights, this stores distances for each pair of
        neighbors.
    """
    start = logg.info('computing neighbors')
    adata = adata.copy() if copy else adata
    if adata.isview:  # we shouldn't need this here...
        adata._init_as_actual(adata.copy())
    neighbors = Neighbors(adata)
    neighbors.compute_neighbors(
        n_neighbors=n_neighbors, knn=knn, n_pcs=n_pcs, use_rep=use_rep,
        method=method, metric=metric, metric_kwds=metric_kwds,
        random_state=random_state,
    )
    adata.uns['neighbors'] = {}
    adata.uns['neighbors']['params'] = {'n_neighbors': n_neighbors, 'method': method}
    adata.uns['neighbors']['params']['metric'] = metric
    if metric_kwds:
        adata.uns['neighbors']['params']['metric_kwds'] = metric_kwds
    if use_rep is not None:
        adata.uns['neighbors']['params']['use_rep'] = use_rep
    if n_pcs is not None:
        adata.uns['neighbors']['params']['n_pcs'] = n_pcs
    adata.uns['neighbors']['distances'] = neighbors.distances
    adata.uns['neighbors']['connectivities'] = neighbors.connectivities
    if neighbors.rp_forest is not None:
        adata.uns['neighbors']['rp_forest'] = neighbors.rp_forest
    logg.info(
        '    finished',
        time=start,
        deep=(
            'added to `.uns[\'neighbors\']`\n'
            '    \'distances\', distances for each pair of neighbors\n'
            '    \'connectivities\', weighted adjacency matrix'
        ),
    )
    return adata if copy else None

class FlatTree(NamedTuple):
    hyperplanes: None
    offsets: None
    children: None
    indices: None

RPForestDict = Mapping[str, Mapping[str, np.ndarray]]

def _rp_forest_generate(rp_forest_dict: RPForestDict) -> Generator[FlatTree, None, None]:
    props = FlatTree._fields[0]
    num_trees = len(rp_forest_dict[props]['start'])-1

    for i in range(num_trees):
        tree = []
        for prop in props:
            start = rp_forest_dict[prop]['start'][i]
            end = rp_forest_dict[prop]['start'][i+1]
            tree.append(rp_forest_dict[prop]['data'][start:end])
        yield FlatTree(*tree)

    tree = []
    for prop in props:
        start = rp_forest_dict[prop]['start'][num_trees]
        tree.append(rp_forest_dict[prop]['data'][start:])
    yield FlatTree(*tree)

def neighbors_update(adata, adata_new, k=10, queue_size=5, random_state=0):
    # only with use_rep='X' for now
    from umap.nndescent import make_initialisations, make_initialized_nnd_search, initialise_search
    from umap.umap_ import INT32_MAX, INT32_MIN
    from umap.utils import deheap_sort
    import umap.distances as dist

    if 'metric_kwds' in adata.uns['neighbors']['params']:
        dist_args = tuple(adata.uns['neighbors']['params']['metric_kwds'].values())
    else:
        dist_args = ()
    dist_func = dist.named_distances[adata.uns['neighbors']['params']['metric']]

    random_init, tree_init = make_initialisations(dist_func, dist_args)
    search = make_initialized_nnd_search(dist_func, dist_args)

    search_graph = adata.uns['neighbors']['distances'].copy()
    search_graph.data = (search_graph.data > 0).astype(np.int8)
    search_graph = search_graph.maximum(search_graph.transpose())
    # prune it?

    random_state = check_random_state(random_state)
    rng_state = random_state.randint(INT32_MIN, INT32_MAX, 3).astype(np.int64)

    if 'rp_forest' in adata.uns['neighbors']:
        rp_forest = _rp_forest_generate(adata.uns['neighbors']['rp_forest'])
    else:
        rp_forest = None
    train = adata.X
    test = adata_new.X

    init = initialise_search(rp_forest, train, test, int(k * queue_size), random_init, tree_init, rng_state)
    result = search(train, search_graph.indptr, search_graph.indices, init, test)

    indices, dists = deheap_sort(result)
    return indices[:, :k], dists[:, :k]

def compute_neighbors_umap(
    X: Union[np.ndarray, csr_matrix],
    n_neighbors: int,
    random_state: Optional[Union[int, RandomState]] = None,
    metric: Union[str, Metric] = 'euclidean',
    metric_kwds: Mapping[str, Any] = {},
    angular: bool = False,
    verbose: bool = False,
):
    """This is from umap.fuzzy_simplicial_set [McInnes18]_.

    Given a set of data X, a neighborhood size, and a measure of distance
    compute the fuzzy simplicial set (here represented as a fuzzy graph in
    the form of a sparse matrix) associated to the data. This is done by
    locally approximating geodesic distance at each point, creating a fuzzy
    simplicial set for each such point, and then combining all the local
    fuzzy simplicial sets into a global one via a fuzzy union.

    Parameters
    ----------
    X: array of shape (n_samples, n_features)
        The data to be modelled as a fuzzy simplicial set.
    n_neighbors
        The number of neighbors to use to approximate geodesic distance.
        Larger numbers induce more global estimates of the manifold that can
        miss finer detail, while smaller values will focus on fine manifold
        structure to the detriment of the larger picture.
    random_state
        A state capable being used as a numpy random state.
    metric
        The metric to use to compute distances in high dimensional space.
        If a string is passed it must match a valid predefined metric. If
        a general metric is required a function that takes two 1d arrays and
        returns a float can be provided. For performance purposes it is
        required that this be a numba jit'd function. Valid string metrics
        include:
            * euclidean
            * manhattan
            * chebyshev
            * minkowski
            * canberra
            * braycurtis
            * mahalanobis
            * wminkowski
            * seuclidean
            * cosine
            * correlation
            * haversine
            * hamming
            * jaccard
            * dice
            * russelrao
            * kulsinski
            * rogerstanimoto
            * sokalmichener
            * sokalsneath
            * yule
        Metrics that take arguments (such as minkowski, mahalanobis etc.)
        can have arguments passed via the metric_kwds dictionary. At this
        time care must be taken and dictionary elements must be ordered
        appropriately; this will hopefully be fixed in the future.
    metric_kwds
        Arguments to pass on to the metric, such as the ``p`` value for
        Minkowski distance.
    angular
        Whether to use angular/cosine distance for the random projection
        forest for seeding NN-descent to determine approximate nearest
        neighbors.
    verbose
        Whether to report information on the current progress of the algorithm.

    Returns
    -------
    **knn_indices**, **knn_dists** : np.arrays of shape (n_observations, n_neighbors)
    """
    from umap.umap_ import nearest_neighbors

    random_state = check_random_state(random_state)

    knn_indices, knn_dists, forest = nearest_neighbors(
        X, n_neighbors, random_state=random_state,
        metric=metric, metric_kwds=metric_kwds,
        angular=angular, verbose=verbose,
    )

    return knn_indices, knn_dists, forest

def get_sparse_matrix_from_indices_distances_umap(knn_indices, knn_dists, n_obs, n_neighbors):
    rows = np.zeros((n_obs * n_neighbors), dtype=np.int64)
    cols = np.zeros((n_obs * n_neighbors), dtype=np.int64)
    vals = np.zeros((n_obs * n_neighbors), dtype=np.float64)

    for i in range(knn_indices.shape[0]):
        for j in range(n_neighbors):
            if knn_indices[i, j] == -1:
                continue  # We didn't get the full knn for i
            if knn_indices[i, j] == i:
                val = 0.0
            else:
                val = knn_dists[i, j]

            rows[i * n_neighbors + j] = i
            cols[i * n_neighbors + j] = knn_indices[i, j]
            vals[i * n_neighbors + j] = val

    result = coo_matrix((vals, (rows, cols)),
                                      shape=(n_obs, n_obs))
    result.eliminate_zeros()
    return result.tocsr()

def compute_connectivities_umap(
    knn_indices, knn_dists,
    n_obs, n_neighbors, set_op_mix_ratio=1.0,
    local_connectivity=1.0,
):
    """This is from umap.fuzzy_simplicial_set [McInnes18]_.

    Given a set of data X, a neighborhood size, and a measure of distance
    compute the fuzzy simplicial set (here represented as a fuzzy graph in
    the form of a sparse matrix) associated to the data. This is done by
    locally approximating geodesic distance at each point, creating a fuzzy
    simplicial set for each such point, and then combining all the local
    fuzzy simplicial sets into a global one via a fuzzy union.
    """
    from umap.umap_ import fuzzy_simplicial_set

    X = coo_matrix(([], ([], [])), shape=(n_obs, 1))
    connectivities = fuzzy_simplicial_set(X, n_neighbors, None, None,
                                          knn_indices=knn_indices, knn_dists=knn_dists,
                                          set_op_mix_ratio=set_op_mix_ratio,
                                          local_connectivity=local_connectivity)
    distances = get_sparse_matrix_from_indices_distances_umap(knn_indices, knn_dists, n_obs, n_neighbors)

    return distances, connectivities.tocsr()

def get_sparse_matrix_from_indices_distances_numpy(indices, distances, n_obs, n_neighbors):
    n_nonzero = n_obs * n_neighbors
    indptr = np.arange(0, n_nonzero + 1, n_neighbors)
    D = csr_matrix((
        distances.copy().ravel(),  # copy the data, otherwise strange behavior here
        indices.copy().ravel(),
        indptr,
    ), shape=(n_obs, n_obs))
    D.eliminate_zeros()
    return D

def get_indices_distances_from_sparse_matrix(D, n_neighbors: int):
    indices = np.zeros((D.shape[0], n_neighbors), dtype=int)
    distances = np.zeros((D.shape[0], n_neighbors), dtype=D.dtype)
    n_neighbors_m1 = n_neighbors - 1
    for i in range(indices.shape[0]):
        neighbors = D[i].nonzero()  # 'true' and 'spurious' zeros
        indices[i, 0] = i
        distances[i, 0] = 0
        # account for the fact that there might be more than n_neighbors
        # due to an approximate search
        # [the point itself was not detected as its own neighbor during the search]
        if len(neighbors[1]) > n_neighbors_m1:
            sorted_indices = np.argsort(D[i][neighbors].A1)[:n_neighbors_m1]
            indices[i, 1:] = neighbors[1][sorted_indices]
            distances[i, 1:] = D[i][
                neighbors[0][sorted_indices], neighbors[1][sorted_indices]]
        else:
            indices[i, 1:] = neighbors[1]
            distances[i, 1:] = D[i][neighbors]
    return indices, distances

def get_indices_distances_from_dense_matrix(D, n_neighbors: int):
    sample_range = np.arange(D.shape[0])[:, None]
    indices = np.argpartition(D, n_neighbors-1, axis=1)[:, :n_neighbors]
    indices = indices[sample_range, np.argsort(D[sample_range, indices])]
    distances = D[sample_range, indices]
    return indices, distances

def _backwards_compat_get_full_X_diffmap(adata: AnnData) -> np.ndarray:
    if 'X_diffmap0' in adata.obs:
        return np.c_[adata.obs['X_diffmap0'].values[:, None],
                     adata.obsm['X_diffmap']]
    else:
        return adata.obsm['X_diffmap']

def _backwards_compat_get_full_eval(adata: AnnData):
    if 'X_diffmap0' in adata.obs:
        return np.r_[1, adata.uns['diffmap_evals']]
    else:
        return adata.uns['diffmap_evals']

def _make_forest_dict(forest):
    d = {}
    props = ('hyperplanes', 'offsets', 'children', 'indices')
    for prop in props:
        d[prop] = {}
        sizes = np.fromiter((getattr(tree, prop).shape[0] for tree in forest), dtype=int)
        d[prop]['start'] = np.zeros_like(sizes)
        if prop == 'offsets':
            dims = sizes.sum()
        else:
            dims = (sizes.sum(), getattr(forest[0], prop).shape[1])
        dtype = getattr(forest[0], prop).dtype
        dat = np.empty(dims, dtype=dtype)
        start = 0
        for i, size in enumerate(sizes):
            d[prop]['start'][i] = start
            end = start+size
            dat[start:end] = getattr(forest[i], prop)
            start = end
        d[prop]['data'] = dat
    return d

class OnFlySymMatrix:
    """Emulate a matrix where elements are calculated on the fly.
    """
    def __init__(
        self,
        get_row: Callable[[Any], np.ndarray],
        shape: Tuple[int, int],
        DC_start: int = 0,
        DC_end: int = -1,
        rows: Optional[Mapping[Any, np.ndarray]] = None,
        restrict_array: Optional[np.ndarray] = None,
    ):
        self.get_row = get_row
        self.shape = shape
        self.DC_start = DC_start
        self.DC_end = DC_end
        self.rows = {} if rows is None else rows
        self.restrict_array = restrict_array  # restrict the array to a subset

    def __getitem__(self, index):
        if isinstance(index, (int, np.integer)):
            if self.restrict_array is None:
                glob_index = index
            else:
                # map the index back to the global index
                glob_index = self.restrict_array[index]
            if glob_index not in self.rows:
                self.rows[glob_index] = self.get_row(glob_index)
            row = self.rows[glob_index]
            if self.restrict_array is None:
                return row
            else:
                return row[self.restrict_array]
        else:
            if self.restrict_array is None:
                glob_index_0, glob_index_1 = index
            else:
                glob_index_0 = self.restrict_array[index[0]]
                glob_index_1 = self.restrict_array[index[1]]
            if glob_index_0 not in self.rows:
                self.rows[glob_index_0] = self.get_row(glob_index_0)
            return self.rows[glob_index_0][glob_index_1]

    def restrict(self, index_array):
        """Generate a view restricted to a subset of indices.
        """
        new_shape = index_array.shape[0], index_array.shape[0]
        return OnFlySymMatrix(
            self.get_row, new_shape, DC_start=self.DC_start,
            DC_end=self.DC_end,
            rows=self.rows, restrict_array=index_array,
        )

class Neighbors:
    """Data represented as graph of nearest neighbors.

    Represent a data matrix as a graph of nearest neighbor relations (edges)
    among data points (nodes).

    Parameters
    ----------
    adata
        Annotated data object.
    n_dcs
        Number of diffusion components to use.
    """

    def __init__(self, adata: AnnData, n_dcs: Optional[int] = None):
        self._adata = adata
        self._init_iroot()
        # use the graph in adata
        info_str = ''
        self.knn: Optional[bool] = None
        self._distances: Union[np.ndarray, csr_matrix, None] = None
        self._connectivities: Union[np.ndarray, csr_matrix, None] = None
        self._transitions_sym: Union[np.ndarray, csr_matrix, None] = None
        self._number_connected_components: Optional[int] = None
        self._rp_forest: Optional[RPForestDict] = None
        if 'neighbors' in adata.uns:
            if 'distances' in adata.uns['neighbors']:
                self.knn = issparse(adata.uns['neighbors']['distances'])
                self._distances = adata.uns['neighbors']['distances']
            if 'connectivities' in adata.uns['neighbors']:
                self.knn = issparse(adata.uns['neighbors']['connectivities'])
                self._connectivities = adata.uns['neighbors']['connectivities']
            if 'rp_forest' in adata.uns['neighbors']:
                self._rp_forest = adata.uns['neighbors']['rp_forest']
            if 'params' in adata.uns['neighbors']:
                self.n_neighbors = adata.uns['neighbors']['params']['n_neighbors']
            else:
                def count_nonzero(a: Union[np.ndarray, csr_matrix]) -> int:
                    return a.count_nonzero() if issparse(a) else np.count_nonzero(a)
                # estimating n_neighbors
                if self._connectivities is None:
                    self.n_neighbors = int(count_nonzero(self._distances) / self._distances.shape[0])
                else:
                    self.n_neighbors = int(count_nonzero(self._connectivities) / self._connectivities.shape[0] / 2)
            info_str += '`.distances` `.connectivities` '
            self._number_connected_components = 1
            if issparse(self._connectivities):
                from scipy.sparse.csgraph import connected_components
                self._connected_components = connected_components(self._connectivities)
                self._number_connected_components = self._connected_components[0]
        if 'X_diffmap' in adata.obsm_keys():
            self._eigen_values = _backwards_compat_get_full_eval(adata)
            self._eigen_basis = _backwards_compat_get_full_X_diffmap(adata)
            if n_dcs is not None:
                if n_dcs > len(self._eigen_values):
                    raise ValueError(
                        'Cannot instantiate using `n_dcs`={}. '
                        'Compute diffmap/spectrum with more components first.'
                        .format(n_dcs))
                self._eigen_values = self._eigen_values[:n_dcs]
                self._eigen_basis = self._eigen_basis[:, :n_dcs]
            self.n_dcs = len(self._eigen_values)
            info_str += '`.eigen_values` `.eigen_basis` `.distances_dpt`'
        else:
            self._eigen_values = None
            self._eigen_basis = None
            self.n_dcs = None
        if info_str != '':
            logg.debug(f'    initialized {info_str}')

    @property
    def rp_forest(self) -> Optional[RPForestDict]:
        return self._rp_forest

    @property
    def distances(self) -> Union[np.ndarray, csr_matrix, None]:
        """Distances between data points (sparse matrix).
        """
        return self._distances

    @property
    def connectivities(self) -> Union[np.ndarray, csr_matrix, None]:
        """Connectivities between data points (sparse matrix).
        """
        return self._connectivities

    @property
    def transitions(self) -> Union[np.ndarray, csr_matrix]:
        """Transition matrix (sparse matrix).

        Is conjugate to the symmetrized transition matrix via::

            self.transitions = self.Z *  self.transitions_sym / self.Z

        where ``self.Z`` is the diagonal matrix storing the normalization of the
        underlying kernel matrix.

        Notes
        -----
        This has not been tested, in contrast to `transitions_sym`.
        """
        if issparse(self.Z):
            Zinv = self.Z.power(-1)
        else:
            Zinv = np.diag(1./np.diag(self.Z))
        return self.Z @ self.transitions_sym @ Zinv

    @property
    def transitions_sym(self) -> Union[np.ndarray, csr_matrix, None]:
        """Symmetrized transition matrix (sparse matrix).

        Is conjugate to the transition matrix via::

            self.transitions_sym = self.Z /  self.transitions * self.Z

        where ``self.Z`` is the diagonal matrix storing the normalization of the
        underlying kernel matrix.
        """
        return self._transitions_sym

    @property
    def eigen_values(self):
        """Eigen values of transition matrix (numpy array).
        """
        return self._eigen_values

    @property
    def eigen_basis(self):
        """Eigen basis of transition matrix (numpy array).
        """
        return self._eigen_basis

    @property
    def distances_dpt(self):
        """DPT distances (on-fly matrix).

        This is yields [Haghverdi16]_, Eq. 15 from the supplement with the
        extensions of [Wolf19]_, supplement on random-walk based distance
        measures.
        """
        return OnFlySymMatrix(self._get_dpt_row, shape=self._adata.shape)

    def to_igraph(self):
        """Generate igraph from connectiviies.
        """
        return utils.get_igraph_from_adjacency(self.connectivities)

    @doc_params(n_pcs=doc_n_pcs, use_rep=doc_use_rep)
    def compute_neighbors(
        self,
        n_neighbors: int = 30,
        knn: bool = True,
        n_pcs: Optional[int] = None,
        use_rep: Optional[str] = None,
        method: str = 'umap',
        random_state: Optional[Union[RandomState, int]] = 0,
        write_knn_indices: bool = False,
        metric: str = 'euclidean',
        metric_kwds: Mapping[str, Any] = {}
    ) -> None:
        """\
        Compute distances and connectivities of neighbors.

        Parameters
        ----------
        n_neighbors
             Use this number of nearest neighbors.
        knn
             Restrict result to `n_neighbors` nearest neighbors.
        {n_pcs}
        {use_rep}

        Returns
        -------
        Writes sparse graph attributes `.distances` and `.connectivities`.
        Also writes `.knn_indices` and `.knn_distances` if
        `write_knn_indices==True`.
        """
        start_neighbors = logg.debug('computing neighbors')
        if n_neighbors > self._adata.shape[0]:  # very small datasets
            n_neighbors = 1 + int(0.5*self._adata.shape[0])
            logg.warning(f'n_obs too small: adjusting to `n_neighbors = {n_neighbors}`')
        if method == 'umap' and not knn:
            raise ValueError('`method = \'umap\' only with `knn = True`.')
        if method not in {'umap', 'gauss'}:
            raise ValueError('`method` needs to be \'umap\' or \'gauss\'.')
        if self._adata.shape[0] >= 10000 and not knn:
            logg.warning('Using high n_obs without `knn=True` takes a lot of memory...')
        self.n_neighbors = n_neighbors
        self.knn = knn
        X = choose_representation(self._adata, use_rep=use_rep, n_pcs=n_pcs)
        # neighbor search
        use_dense_distances = (metric == 'euclidean' and X.shape[0] < 8192) or knn == False
        if use_dense_distances:
            _distances = pairwise_distances(X, metric=metric, **metric_kwds)
            knn_indices, knn_distances = get_indices_distances_from_dense_matrix(
                _distances, n_neighbors)
            if knn:
                self._distances = get_sparse_matrix_from_indices_distances_numpy(
                    knn_indices, knn_distances, X.shape[0], n_neighbors)
            else:
                self._distances = _distances
        else:
            # non-euclidean case and approx nearest neighbors
            if X.shape[0] < 4096:
                X = pairwise_distances(X, metric=metric, **metric_kwds)
                metric = 'precomputed'
            knn_indices, knn_distances, _ = compute_neighbors_umap(
                X, n_neighbors, random_state, metric=metric, metric_kwds=metric_kwds)
            #self._rp_forest = _make_forest_dict(forest)
        # write indices as attributes
        if write_knn_indices:
            self.knn_indices = knn_indices
            self.knn_distances = knn_distances
        start_connect = logg.debug('computed neighbors', time=start_neighbors)
        if not use_dense_distances or method == 'umap':
            # we need self._distances also for method == 'gauss' if we didn't
            # use dense distances
            self._distances, self._connectivities = compute_connectivities_umap(
                knn_indices,
                knn_distances,
                self._adata.shape[0],
                self.n_neighbors,
            )
        # overwrite the umap connectivities if method is 'gauss'
        # self._distances is unaffected by this
        if method == 'gauss':
            self._compute_connectivities_diffmap()
        logg.debug('computed connectivities', time=start_connect)
        self._number_connected_components = 1
        if issparse(self._connectivities):
            from scipy.sparse.csgraph import connected_components
            self._connected_components = connected_components(self._connectivities)
            self._number_connected_components = self._connected_components[0]

    def _compute_connectivities_diffmap(self, density_normalize=True):
        # init distances
        if self.knn:
            Dsq = self._distances.power(2)
            indices, distances_sq = get_indices_distances_from_sparse_matrix(
                Dsq, self.n_neighbors)
        else:
            Dsq = np.power(self._distances, 2)
            indices, distances_sq = get_indices_distances_from_dense_matrix(
                Dsq, self.n_neighbors)

        # exclude the first point, the 0th neighbor
        indices = indices[:, 1:]
        distances_sq = distances_sq[:, 1:]

        # choose sigma, the heuristic here doesn't seem to make much of a difference,
        # but is used to reproduce the figures of Haghverdi et al. (2016)
        if self.knn:
            # as the distances are not sorted
            # we have decay within the n_neighbors first neighbors
            sigmas_sq = np.median(distances_sq, axis=1)
        else:
            # the last item is already in its sorted position through argpartition
            # we have decay beyond the n_neighbors neighbors
            sigmas_sq = distances_sq[:, -1]/4
        sigmas = np.sqrt(sigmas_sq)

        # compute the symmetric weight matrix
        if not issparse(self._distances):
            Num = 2 * np.multiply.outer(sigmas, sigmas)
            Den = np.add.outer(sigmas_sq, sigmas_sq)
            W = np.sqrt(Num/Den) * np.exp(-Dsq/Den)
            # make the weight matrix sparse
            if not self.knn:
                mask = W > 1e-14
                W[mask == False] = 0
            else:
                # restrict number of neighbors to ~k
                # build a symmetric mask
                mask = np.zeros(Dsq.shape, dtype=bool)
                for i, row in enumerate(indices):
                    mask[i, row] = True
                    for j in row:
                        if i not in set(indices[j]):
                            W[j, i] = W[i, j]
                            mask[j, i] = True
                # set all entries that are not nearest neighbors to zero
                W[mask == False] = 0
        else:
            W = Dsq.copy()  # need to copy the distance matrix here; what follows is inplace
            for i in range(len(Dsq.indptr[:-1])):
                row = Dsq.indices[Dsq.indptr[i]:Dsq.indptr[i+1]]
                num = 2 * sigmas[i] * sigmas[row]
                den = sigmas_sq[i] + sigmas_sq[row]
                W.data[Dsq.indptr[i]:Dsq.indptr[i+1]] = np.sqrt(num/den) * np.exp(
                    -Dsq.data[Dsq.indptr[i]: Dsq.indptr[i+1]] / den)
            W = W.tolil()
            for i, row in enumerate(indices):
                for j in row:
                    if i not in set(indices[j]):
                        W[j, i] = W[i, j]
            W = W.tocsr()

        self._connectivities = W

    def compute_transitions(self, density_normalize=True):
        """Compute transition matrix.

        Parameters
        ----------
        density_normalize : `bool`
            The density rescaling of Coifman and Lafon (2006): Then only the
            geometry of the data matters, not the sampled density.

        Returns
        -------
        Makes attributes `.transitions_sym` and `.transitions` available.
        """
        start = logg.info('computing transitions')
        W = self._connectivities
        # density normalization as of Coifman et al. (2005)
        # ensures that kernel matrix is independent of sampling density
        if density_normalize:
            # q[i] is an estimate for the sampling density at point i
            # it's also the degree of the underlying graph
            q = np.asarray(W.sum(axis=0))
            if not issparse(W):
                Q = np.diag(1.0/q)
            else:
                Q = scipy.sparse.spdiags(1.0/q, 0, W.shape[0], W.shape[0])
            K = Q @ W @ Q
        else:
            K = W

        # z[i] is the square root of the row sum of K
        z = np.sqrt(np.asarray(K.sum(axis=0)))
        if not issparse(K):
            self.Z = np.diag(1.0/z)
        else:
            self.Z = scipy.sparse.spdiags(1.0/z, 0, K.shape[0], K.shape[0])
        self._transitions_sym = self.Z @ K @ self.Z
        logg.info('    finished', time=start)

    def compute_eigen(self, n_comps=15, sym=None, sort='decrease'):
        """Compute eigen decomposition of transition matrix.

        Parameters
        ----------
        n_comps : `int`
            Number of eigenvalues/vectors to be computed, set `n_comps = 0` if
            you need all eigenvectors.
        sym : `bool`
            Instead of computing the eigendecomposition of the assymetric
            transition matrix, computed the eigendecomposition of the symmetric
            Ktilde matrix.
        matrix : sparse matrix, np.ndarray, optional (default: `.connectivities`)
            Matrix to diagonalize. Merely for testing and comparison purposes.

        Returns
        -------
        Writes the following attributes.

        eigen_values : numpy.ndarray
            Eigenvalues of transition matrix.
        eigen_basis : numpy.ndarray
             Matrix of eigenvectors (stored in columns).  `.eigen_basis` is
             projection of data matrix on right eigenvectors, that is, the
             projection on the diffusion components.  these are simply the
             components of the right eigenvectors and can directly be used for
             plotting.
        """
        np.set_printoptions(precision=10)
        if self._transitions_sym is None:
            raise ValueError('Run `.compute_transitions` first.')
        matrix = self._transitions_sym
        # compute the spectrum
        if n_comps == 0:
            evals, evecs = scipy.linalg.eigh(matrix)
        else:
            n_comps = min(matrix.shape[0]-1, n_comps)
            # ncv = max(2 * n_comps + 1, int(np.sqrt(matrix.shape[0])))
            ncv = None
            which = 'LM' if sort == 'decrease' else 'SM'
            # it pays off to increase the stability with a bit more precision
            matrix = matrix.astype(np.float64)
            evals, evecs = scipy.sparse.linalg.eigsh(matrix, k=n_comps,
                                                  which=which, ncv=ncv)
            evals, evecs = evals.astype(np.float32), evecs.astype(np.float32)
        if sort == 'decrease':
            evals = evals[::-1]
            evecs = evecs[:, ::-1]
        logg.info('    eigenvalues of transition matrix\n'
                  '    {}'.format(str(evals).replace('\n', '\n    ')))
        if self._number_connected_components > len(evals)/2:
            logg.warning('Transition matrix has many disconnected components!')
        self._eigen_values = evals
        self._eigen_basis = evecs

    def _init_iroot(self):
        self.iroot = None
        # set iroot directly
        if 'iroot' in self._adata.uns:
            if self._adata.uns['iroot'] >= self._adata.n_obs:
                logg.warning(
                    f'Root cell index {self._adata.uns["iroot"]} does not '
                    f'exist for {self._adata.n_obs} samples. It’s ignored.'
                )
            else:
                self.iroot = self._adata.uns['iroot']
            return
        # set iroot via xroot
        xroot = None
        if 'xroot' in self._adata.uns: xroot = self._adata.uns['xroot']
        elif 'xroot' in self._adata.var: xroot = self._adata.var['xroot']
        # see whether we can set self.iroot using the full data matrix
        if xroot is not None and xroot.size == self._adata.shape[1]:
            self._set_iroot_via_xroot(xroot)

    def _get_dpt_row(self, i):
        mask = None
        if self._number_connected_components > 1:
            label = self._connected_components[1][i]
            mask = self._connected_components[1] == label
        row = sum(
            (
                self.eigen_values[l] / (1-self.eigen_values[l])
                * (self.eigen_basis[i, l] - self.eigen_basis[:, l])
            )**2
            # account for float32 precision
            for l in range(0, self.eigen_values.size)
            if self.eigen_values[l] < 0.9994
        )
        # thanks to Marius Lange for pointing Alex to this:
        # we will likely remove the contributions from the stationary state below when making
        # backwards compat breaking changes, they originate from an early implementation in 2015
        # they never seem to have deteriorated results, but also other distance measures (see e.g.
        # PAGA paper) don't have it, which makes sense
        row += sum(
            (self.eigen_basis[i, l] - self.eigen_basis[:, l])**2
            for l in range(0, self.eigen_values.size)
            if self.eigen_values[l] >= 0.9994
        )
        if mask is not None:
            row[~mask] = np.inf
        return np.sqrt(row)

    def _set_pseudotime(self):
        """Return pseudotime with respect to root point.
        """
        self.pseudotime = self.distances_dpt[self.iroot].copy()
        self.pseudotime /= np.max(self.pseudotime[self.pseudotime < np.inf])

    def _set_iroot_via_xroot(self, xroot):
        """Determine the index of the root cell.

        Given an expression vector, find the observation index that is closest
        to this vector.

        Parameters
        ----------
        xroot : np.ndarray
            Vector that marks the root cell, the vector storing the initial
            condition, only relevant for computing pseudotime.
        """
        if self._adata.shape[1] != xroot.size:
            raise ValueError(
                'The root vector you provided does not have the '
                'correct dimension.')
        # this is the squared distance
        dsqroot = 1e10
        iroot = 0
        for i in range(self._adata.shape[0]):
            diff = self._adata.X[i, :] - xroot
            dsq = diff @ diff
            if dsq < dsqroot:
                dsqroot = dsq
                iroot = i
                if np.sqrt(dsqroot) < 1e-10: break
        logg.debug(f'setting root index to {iroot}')
        if self.iroot is not None and iroot != self.iroot:
            logg.warning(f'Changing index of iroot from {self.iroot} to {iroot}.')
        self.iroot = iroot

max_stars_repo_name: zeehio/scanpy

max_issues_repo_name: waldoPHD/scanpy

max_forks_repo_name: waldoPHD/scanpy

reason: import numpy,from numpy,import scipy,from scipy


text:
import tarmac.tarmac as tm
import numpy as np
import pytest
import matplotlib.pyplot as plt

def test_corner_plot():
    """Test the corner plot function.

    Construct a corner plot using the example data.

    """

    data = np.load("extra_data.npy")

    fig = plt.figure(figsize=(10, 10))

    tm.corner_plot(fig,
                   data,
                   labels=['a','b','c','d'],
                   cmap='viridis'
                   )
    plt.show()
    return

def test_walker_trace():
    """Test the walker trace function.

    Construct a walker trace plot using the example data.

    """

    data = np.load("extra_data.npy")

    fig = plt.figure(figsize=(10, 10))
    tm.walker_trace(fig, data, labels=["a", "b", "c", "d"])
    plt.show()
    return

max_stars_repo_name: peytondmurray/tarmac

max_issues_repo_name: peytondmurray/tarmac

max_forks_repo_name: peytondmurray/tarmac


text:
import numpy as np
import math
import cv2

import extramath as em

class Joint:
	def __init__(self, label, degMin, degMax, index, end, base = 0):
		self.label = label
		self.degMin = degMin
		self.degMax = degMax
		self.index = index
		self.end = end
		self.base = base

	def getCurlAngle(self, lm):
		return em.angleByPoints(
			lm[self.index],
			lm[self.end],
			lm[self.base])

	def getCurlPercentage(self, lm):
		return self.angleToPercent(self.getCurlAngle(lm))

	def angleToPercent(self, angle):
		p = em.map(angle, self.degMin, self.degMax, 100, 0)
		return np.clip(p, 0, 100)

fingers = [
Joint('Thumb', 120, 150, 1, 4),
Joint('Index', 50, 170, 5, 8),
Joint('Middle', 30, 170, 9, 12),
Joint('Ring', 30, 170, 13, 16),
Joint('Pinky', 60, 170, 17, 20),
]

def getJointDataString(landmarks):
	result = "### Joint Curl Data: ###\n"
	resultRaw = ""

	for joint in fingers:
		angle = np.rint(joint.getCurlAngle(landmarks))
		percentage = np.rint(joint.angleToPercent(angle))
		start = f"| {joint.label}:"

		result += start + f"{percentage}%".center(8) + '|'
		resultRaw += start + f"{angle}º".center(8) + '|'

	return result + '\n' + resultRaw


text:
from typing import List, Optional

from .baseLinearReward import BaseLinearReward
import numpy as np
import tensorflow.compat.v1 as tf
import math

from proseco.utility.constants import GRAVITY, ACTION_DURATION

tf.disable_v2_behavior()

class FeatureNames:
    """
    Enum defining the feature indexes in the feature vec
    """

    DIFF_VEL_VEL_DES = 0
    ABS_LANE_DIFF = 1
    DIFF_DES_LANE_CENT = 2
    COLLIDED = 3
    INVALID_STATE = 4
    INVALID_ACTION = 5
    ACC_Y = 6

class LinearIrlReward(BaseLinearReward):
    """
    Class which specifies the LinearIRL function (in particular its features) which corresponds to the same class in the MCTS

    Instance variables:
        feature_names list -- list of the names of all used features in this reward function
        expert_cost_params np.array -- parameters of the reward function for the expert trajectories
        number_of_features int -- number of features in the reward function
        expert_trajectories_pickle_path string -- path for input and output of the expert trajectory messages, dependend if experts are generated or training is proceeded
        name string -- name of the reward
        mcts_cost_model_name -- name of the cost class inside the MCTS which mirrors the features
        T float -- parameter of the reward function that is normally equal to the episode length of the trajectories inside the IRL procedure (must be set by the irl model before training starts - this is also asserted in create_cost_message)
    """

    def __init__(
        self,
        expert_trajectories_pickle_path: str,
        expert_cost_params: Optional[List[float]],
    ):
        self.name = "LinearIrl"
        self.mcts_cost_model_name = "costLinear"
        self.expert_trajectories_pickle_path = expert_trajectories_pickle_path
        self.T = None
        self.feature_names = [
            "diff_vel_vel_des",
            "abs_lane_diff",
            "diff_des_lane_cent",
            "collided",
            "invalidState",
            "invalidAction",
            "accY",
        ]
        self.number_of_features = len(self.feature_names)

        if expert_cost_params:
            self.expert_cost_params = np.array(
                expert_cost_params, ndmin=2
            ).T  # Transpose so that we get the same shape as linearIrl.initial_parameters
        else:
            self.expert_cost_params = None

    def get_name(self):
        """
        return name of the reward model

        Return:
            string -- name of the reward model
        """
        return self.name

    def is_cooperative(self):
        return False

    def get_feature_names(self):
        """
        returns the names of the features

        Return:
            list -- names of features
        """
        return self.feature_names

    def get_number_of_features(self):
        """
        return number of features == number of input neurons in the reward neural net

        Return:
            int -- number of input features
        """
        return self.number_of_features

    def get_expert_parameters(self):
        """
        return the parameters of the expert reward model to create experts

        Return:
            np.array -- parameters create expert trajectories with
        """
        return self.expert_cost_params

    def get_expert_pickle_folder(self):
        """
        return the path to the pickle files where the expert trajectory messages are saved or will be written to (important for training + expert creation)

        Return:
            string -- path where the pickle lie or should put to
        """
        return self.expert_trajectories_pickle_path

    ##### Main methods #############

    def features_to_vec(
        self, features, features_prev, meta_info, feature_vec_other_agents
    ):
        """Extracts features from the feature message of one agent in one timstep

        Arguments:
            features ros_proseco_planning.msg.Features -- feature message for one agent in timestep t
            features_prev ros_proseco_planning.msg.Features -- feature message for one agent in timestep t-1
            meta_info ros_proseco_planning.msg.ScenarioInfo -- scenarioInfo - necessary to get for example lanewidth

        Returns:
            np.array -- feature valus in list
        """

        feature_vec = [0] * self.number_of_features
        lane_width = meta_info.laneWidth

        feature_vec[FeatureNames.DIFF_VEL_VEL_DES] = self.feature_vel_vel_des(
            features.diff_vel_vel_des, features.desired_vel
        )
        feature_vec[FeatureNames.ABS_LANE_DIFF] = self.feature_abs_lane_diff(
            features.abs_lane_diff
        )
        feature_vec[FeatureNames.DIFF_DES_LANE_CENT] = self.feature_diff_des_lane_cent(
            features.diff_des_lane_cent, lane_width
        )
        feature_vec[FeatureNames.COLLIDED] = int(features.collided)
        feature_vec[FeatureNames.INVALID_STATE] = int(features.invalidState)
        feature_vec[FeatureNames.INVALID_ACTION] = self.feature_invalid_action(
            features.invalidAction
        )
        feature_vec[FeatureNames.ACC_Y] = self.feature_acc_y(features.accY)
        return np.array(feature_vec)

    def create_cost_message(self, cost_params):
        """Set current cost message which is retrieved to the slaves via the service

        Arguments:
            params {np.array} -- cost parameter vector

        Return:
            json -- cost model
        """

        return {
            "name": [self.mcts_cost_model_name],
            "w_acceleration_x": [0.0],
            "w_acceleration_y": cost_params[FeatureNames.ACC_Y],
            "cost_collision": cost_params[FeatureNames.COLLIDED],
            "cost_invalid_action": cost_params[FeatureNames.INVALID_ACTION],
            "cost_invalid_state": cost_params[FeatureNames.INVALID_STATE],
            "w_lane_center_deviation": cost_params[FeatureNames.DIFF_DES_LANE_CENT],
            "w_lane_change": [0.0],
            "w_lane_deviation": cost_params[FeatureNames.ABS_LANE_DIFF],
            "w_velocity_deviation": cost_params[FeatureNames.DIFF_VEL_VEL_DES],
            "cost_enter_safe_range": [0.0],
            "reward_terminal": [0.0],
            "T": [self.T],
        }

    def feature_acc_y(self, cumulated_squared_acceleration: float) -> float:
        """Calculates the acceleration y feature.

        Parameters
        ----------
        cumulated_squared_acceleration : float
            The Integral of the squared acceleration in y direction.

        Returns
        -------
        float
            The value of the acceleration y feature.
        """
        assert cumulated_squared_acceleration >= 0
        g_y = math.sqrt(cumulated_squared_acceleration / ACTION_DURATION) / GRAVITY
        return max(1.0 - 4 * g_y, -1.0) / self.T

    def feature_vel_vel_des(
        self, velocity_deviation: float, desired_vel: float
    ) -> float:
        """Calculates the velocity deviation feature.

        Parameters
        ----------
        velocity_deviation : float
            The difference between the current velocity and the desired velocity.
        desired_vel : float
            The desired velocity.

        Returns
        -------
        float
            The value of the velocity deviation feature.
        """
        return (
            max(1.0 - abs(velocity_deviation) / (abs(desired_vel) / 10.0), -1.0)
            / self.T
        )

    def feature_abs_lane_diff(self, lane_deviation: int) -> float:
        """Calculates the lane deviation feature.

        Parameters
        ----------
        lane_deviation : int
            The difference between the current lane and the desired lane.

        Returns
        -------
        float
            The value of the lane deviation feature.
        """
        return max(1.0 - abs(lane_deviation), -1.0) / self.T

    def feature_diff_des_lane_cent(self, d: float, lane_width: float) -> float:
        """Calculates the desired lane center feature. The feature is positive as long as the deviation is smaller than a quarter of the lane width.

        Parameters
        ----------
        d : float
            The difference between the y position and the center of the lane.
        lane_width : [type]
            The width of the lane.

        Returns
        -------
        float
            The value of the desired lane center feature.
        """
        return max(1.0 - abs(d) / (lane_width / 4.0), -1.0) / self.T

    def feature_invalid_action(self, is_invalid) -> float:
        """Calculates the invalid action feature.

        Parameters
        ----------
        is_invalid : bool
            Whether the action is invalid or not.

        Returns
        -------
        float
            The value of the invalid action feature.
        """
        if is_invalid:
            return 1.0 / self.T
        else:
            return 0.0


text:
import matplotlib.lines as mlines
import matplotlib.patches as mpatches
from utils import generate_u3
from estimators import QiskitAnsatz, QiskitFlammiaAnsatz, QiskitEstimator, QiskitFlammiaEstimator
from probability_distributions import ChiProbDist, FlammiaProbDist
from utils import GateObj, U_from_hamiltonian, get_filename
from qutip import Qobj, cnot, tensor, qeye
from qiskit.aqua import QuantumInstance
from qiskit import IBMQ, Aer, execute, QuantumRegister, ClassicalRegister, QuantumCircuit
import pickle
from os import path, getcwd, makedirs
import sys
import matplotlib.pyplot as plt
zero_fidelsimport numpy as np

savepath = path.join(getcwd(), 'data', 'projective_estimations', 'real_machine',
                     '3q_random')

with open(path.join(savepath, 'process_fidelity_estimates.pickle'), 'rb') as f:
    proc_fidels = pickle.load(f)

with open(path.join(savepath, 'zero_fidelity_estimates.pickle'), 'rb') as f:
    zero_fidesl = pickle.load(f)

point_alpha = 0.8
fill_alpha = 0.12

x = np.linspace(0, 1, 11)
y = np.linspace(0, 1, 11)

plt.scatter([i for i in range(len(proc_fidels))], [
            f for f in proc_fidels], marker='x', s=15, alpha=point_alpha)
plt.fill_between([-1 + i for i in range(len(proc_fidels)+2)], f_mean-np.std(proc_fidels),
                 f_mean+np.std(proc_fidels), alpha=fill_alpha, color=col_cyc[0])
plt.hlines(np.mean(proc_fidels), -1, len(proc_fidels), lw=0.7, color=col_cyc[0])
p1 = mlines.Line2D([], [], color=col_cyc[0], marker='x', markersize=np.sqrt(15), lw=0.7)
p2 = mpatches.Patch(color=col_cyc[0], alpha=fill_alpha, linewidth=0)

plt.scatter([i for i in range(len(zero_fidels))], [
            f for f in zero_fidels], marker='^', s=15, alpha=point_alpha)
plt.fill_between([-1 + i for i in range(len(zero_fidels)+2)], fom_mean-np.std(zero_fidels),
                 fom_mean+np.std(zero_fidels), alpha=fill_alpha, color=col_cyc[1])
plt.hlines(np.mean(zero_fidels), -1, len(zero_fidels), lw=0.7, color=col_cyc[1])
p3 = mlines.Line2D([], [], color=col_cyc[1], marker='^', markersize=np.sqrt(15), lw=0.7)
p4 = mpatches.Patch(color=col_cyc[1], alpha=fill_alpha, linewidth=0)

plt.xlabel("Evaluation")
plt.ylabel("Estimated value")

plt.legend(((p1, p2), (p3, p4),), ('Process Fidelity', r"$0-$Fidelity",), frameon=True)
plt.savefig(path.join(savepath, "toronto_50evals_129024_exps_mean_stds.png"))