text:
# This file contains information about the different table-
# based file types and methods for validating them
import logging
import re
import os
from functools import reduce

import pandas as pd
import numpy as np

from django.conf import settings
from django.core.paginator import Paginator, Page
from rest_framework.pagination import PageNumberPagination

from .base import DataResource, ParseException, UnexpectedTypeValidationException
from api.data_structures import Feature, \
    FeatureSet, \
    Observation, \
    ObservationSet, \
    create_attribute, \
    convert_dtype, \
    numeric_attribute_typenames
from api.utilities.basic_utils import alert_admins
from api.serializers.feature_set import FeatureSetSerializer
from api.serializers.observation_set import ObservationSetSerializer

logger = logging.getLogger(__name__)

# acceptable file extensions which give us a
# clue about how to parse the files.
CSV = 'csv'
TSV = 'tsv'
TAB = 'tab'
BED = 'bed'
VCF = 'vcf'
XLS = 'xls'
XLSX = 'xlsx'
TAB_DELIMITED_EXTENSIONS = [TSV, TAB, BED, VCF]
COMMA_DELIMITED_EXTENSIONS = [CSV]
EXCEL_EXTENSIONS = [XLS, XLSX]

class ParserNotFoundException(Exception):
    '''
    For raising exceptions when a proper
    parser cannot be found.
    '''
    pass

# Some error messages:
PARSE_ERROR = ('There was an unexpected problem when'
    ' parsing and validating the file.')

PARSER_NOT_FOUND_ERROR = ('Could not find an appropriate parser'
    ' for the resource.  Please check the instructions.')

NON_NUMERIC_ERROR = 'Some of your columns contained non-numeric entries. Please check the file.'

NON_INTEGER_ERROR = 'Some of your columns contained non-integer entries. Please check the file.'

TRIVIAL_TABLE_ERROR = ('The file contained only a single column'
    ' which provided an index.  No data was provided in additional columns.')

BED_FORMAT_ERROR = ('When parsing the BED file, we detected issues with'
    ' column(s): {cols}. Note that BED files must NOT have column headers and can'
    ' contain only integers in the second and third columns, which correspond'
    ' to the start and end of a genomic location.  Please check your entries and'
    ' ensure your file does not have a header line.')

NUMBERED_COLUMN_NAMES_ERROR = ('All the column names were numbers.  Often this is due'
    ' to a missing column header.  In that case, there will be a missing row in'
    ' your table. If you named your columns with numbers, please change them'
    ' to something else to avoid incorrect parsing of the file.')

NUMBERED_ROW_NAMES_ERROR = ('All the row names were numbers.  We use the first column'
    ' to uniquely identify the rows for filtering purposes. If you named your rows with'
    ' numbers, please change them to something else (e.g. add "x" to the beginning'
    ' to avoid incorrect parsing of the file.')

NONUNIQUE_ROW_NAMES_ERROR = ('Your row names were not unique, which could cause'
    ' unexpected behavior.')

MISSING_HEADER_WARNING = ('One of your column names matched the values in the'
    ' corresponding problem.  This is not an error, but may indicate that a'
    ' proper header line was missing.  Please check to ensure the file was'
    ' parsed correctly.')

EMPTY_TABLE_ERROR = ('The parsed table was empty. This can happen if you are'
    ' importing a table with only a single column. Alternatively, if you are'
    ' trying to import an Excel spreadsheet, please ensure that the data is'
    ' contained in the first sheet of the workbook.')

def col_str_formatter(x):
    '''
    x is a tuple with the column number
    and column name
    '''
    return '%s (column %d)' % (x[0],x[1])

class TableResourcePage(Page):
    '''
    Overrides some methods of the django.core.paginator.Page
    class. The default version does a cast to a list, which is a problem
    when the "query set" is a dataframe instead of a QuerySet instance
    '''
    def __getitem__(self, index):
        if not isinstance(index, (int, slice)):
            raise TypeError(
                'Page indices must be integers or slices, not %s.'
                % type(index).__name__
            )

        return self.object_list[index]

class TableResourcePaginator(Paginator):
    '''
    This overrides the default Django paginator, which
    typically expects to deal with a database query (a QuerySet).

    Dataframes are list-like (sort of...) but the problem is that
    the default django paginator ultimately does a cast to a list. When you
    cast a pandas dataframe as a list, you only get a list of the
    column names, NOT a list of the "records" (rows).

    The django default Paginator mentions using the _get_page method as a
    way to override the django.core.paginator.Page class (which is where the
    cast occurs). By overriding that method, we can make our own Page-like
    class.

    Finally, this class can be provided to a subclass of DRF's Paginator
    classes.
    '''
    def _get_page(self, *args, **kwargs):
        '''
        Method to override the default behavior of the django core
        Paginator class. Should return a Page-like class.
        '''
        return TableResourcePage(*args, **kwargs)

class TableResourcePageNumberPagination(PageNumberPagination):
    django_paginator_class = TableResourcePaginator
    page_size_query_param = settings.PAGE_SIZE_PARAM

class TableResource(DataResource):
    '''
    The `TableResource` is the most generic form of a delimited file.  Any
    type of data that can be represented as rows and columns.

    This or any of the more specific subclasses can be contained in files
    saved in CSV, TSV, or Excel (xls/xlsx) format.  If in Excel format, the
    data of interest must reside in the first sheet of the workbook.

    Special tab-delimited files like BED or VCF files are recognized by
    their canonical extension (e.g. ".bed" or ".vcf").

    Note that unless you create a "specialized" implementation (e.g. like
    for a BED file), then we assume you have features as rows and observables
    as columns.
    '''
    ACCEPTABLE_EXTENSIONS = [
        CSV,
        TSV,
        TAB,
        BED,
        VCF,
        XLS,
        XLSX
    ]

    # the "standardized" format we will save all table-based files as:
    STANDARD_FORMAT = TSV

    # Create a list of query params that are "reserved" and we ignore when
    # attempting to filter on the actual table content (e.g. the column/rows)
    IGNORED_QUERY_PARAMS = [settings.PAGE_SIZE_PARAM, settings.PAGE_PARAM, settings.SORT_PARAM]

    def __init__(self):
        self.table = None

    @staticmethod
    def get_paginator():
        return TableResourcePageNumberPagination()

    @staticmethod
    def get_reader(resource_path):
        '''
        By using the file extension, we infer the delimiter
        Returns a pandas "reader" (e.g. `read_csv` or `read_table`)
        '''
        file_extension = DataResource.get_extension(resource_path)

        if file_extension in COMMA_DELIMITED_EXTENSIONS:
            return pd.read_csv
        elif file_extension in TAB_DELIMITED_EXTENSIONS:
            return pd.read_table
        elif file_extension in EXCEL_EXTENSIONS:
            return pd.read_excel
        else:
            logger.error('Could not infer the file format from the file'
            ' extension of {ext}.  Full resource path was {path}'.format(
                ext = file_extension,
                path = resource_path
            ))
            return None

    @staticmethod
    def index_all_numbers(names):
        '''
        Works for both row and column indexes.  Returns
        True if all the index labels are numbers.
        '''
        if all([re.match('\d+', str(x)) for x in names]):
            return True
        else:
            return False

    def read_resource(self, resource_path):
        '''
        One common spot to define how the file is read
        '''
        reader = TableResource.get_reader(resource_path)
        if reader is None:
            raise ParserNotFoundException('')
        else:
            try:
                # read the table using the appropriate parser:
                self.table = reader(resource_path, index_col=0, comment='#')

                # drop extra/empty cols and rows
                self.table.dropna(axis=0, how='all', inplace=True)
                self.table.dropna(axis=1, how='all', inplace=True)

                # call a method to
            except Exception as ex:
                logger.info('Could not use {reader} to parse the file'
                ' at {path}'.format(
                    reader = reader,
                    path = resource_path
                ))
                raise ParseException('Failed when parsing the table-based resource.')

    def performs_validation(self):
        '''
        Since we have methods to validate table-based DataResource types, we
        override this method and return True, which indicates that we CAN
        perform validation.
        '''
        return True

    def validate_type(self, resource_path):
        '''
        In this base method, we determine attempt to parse the file.
        If there are no restrictions on content, succesful parsing
        of the file is good enough.

        More specific constraints on the file content are handled in child
        classes.  This method, however, fills in the `self.table` member
        which is then accessible to children.
        '''
        try:
            self.read_resource(resource_path)
            if self.table.shape == (0,0):
                return (False, EMPTY_TABLE_ERROR )

            if self.table.shape[1] == 0:
                return (False, TRIVIAL_TABLE_ERROR)

            # check if all the column names are numbers-- which would USUALLY
            # indicate a missing header
            columns_all_numbers = TableResource.index_all_numbers(self.table.columns)
            if columns_all_numbers:
                return (False, NUMBERED_COLUMN_NAMES_ERROR)

            # check if all the rownames are numbers, which would usually
            # indicate missing row names (i.e. a column of data is read
            # as the index)
            rows_all_numbers = TableResource.index_all_numbers(self.table.index)

            if rows_all_numbers:
                return (False, NUMBERED_ROW_NAMES_ERROR)

            # check for duplicate row names
            if self.table.index.has_duplicates:
                return (False, NONUNIQUE_ROW_NAMES_ERROR)
            return (True, None)

        except ParserNotFoundException as ex:
            return (False, PARSER_NOT_FOUND_ERROR)

        except ParseException as ex:
            return (False, PARSE_ERROR)

    def do_type_cast(self, v, typename):
        '''
        Used for casting the type when query params are provided.
        '''
        if typename in numeric_attribute_typenames:
            try:
                val = float(v)
            except ValueError as ex:
                raise ParseException('Could not parse "{v}" as a number.'.format(v=v))
        else:
            val = v
        return val

    def perform_sorting(self, query_params):
        '''
        Sorts the table as requested. Sorts self.table in place

        A general sort string (after the equals sign) would be:
        <sort kw>:<col>,<sort kw>:<col>,...
        e.g. for sorting q-value (padj) ascending followed by fold-change descending,
        [asc]:padj,[desc]:log2Foldchange
        '''
        if settings.SORT_PARAM in query_params:
            sort_strings = query_params[settings.SORT_PARAM].split(',')
            sort_order_list = []
            column_list = []
            for s in sort_strings:
                try:
                    sort_order, col = s.split(settings.QUERY_PARAM_DELIMITER)
                except ValueError:
                    raise ParseException('The sorting request was not properly formatted. '
                        'Please use "<ordering keyword>:<column name>"')
                if sort_order in settings.SORTING_OPTIONS:
                    sort_order_list.append(sort_order)
                else:
                    raise ParseException('The sort order "{s}" is not an available option. Choose from: {opts}'.format(
                        s = sort_order,
                        opts = ','.join(settings.SORTING_OPTIONS)
                    ))

                if col in self.table.columns:
                    column_list.append(col)
                else:
                    raise ParseException('The column identifier "{s}" does not exist in this resource. Options are: {opts}'.format(
                        s = col,
                        opts = ','.join(self.table.columns)
                    ))
            # at this point, all the sort orders and cols were OK. Now perform the sorting:
            # Need to convert our strings (e.g. "[asc]") to bools for the pandas sort_values method.
            order_bool = [True if x==settings.ASCENDING else False for x in sort_order_list]
            self.table.sort_values(by=column_list, ascending=order_bool, inplace=True)

    def filter_against_query_params(self, query_params):
        '''
        Looks through the query params to subset the table
        '''
        table_cols = self.table.columns

        # guard against some edge case where the table we are filtering happens to have
        # columns that conflict with the pagination parameters. We simply inform the admins
        # and ignore that conflict by not using that filter

        if any([x in self.IGNORED_QUERY_PARAMS for x in table_cols]):
            logger.warning('One of the column names conflicted with the pagination query params.')
            alert_admins('Edge-case error: when filtering on a column, one of the column names'
                ' conflicted with the pagination query params. Query params was: {p} and column'
                ' names were: {c}'.format(
                    p = query_params,
                    c = ','.join(table_cols)
                )
            )
        filters = []

        # used to map the pandas native type to a MEV-type so we can do type casting consistently
        type_dict = self.get_type_dict()
        for k,v in query_params.items():
            split_v = v.split(settings.QUERY_PARAM_DELIMITER)
            if (not k in self.IGNORED_QUERY_PARAMS) and (k in table_cols):
                # v is either a value (in the case of strict equality)
                # or a delimited string which will dictate the comparison.
                # For example, to filter on the 'pval' column for values less than or equal to 0.01,
                # v would be "[lte]:0.01". The "[lte]" string is set in our general settings file.
                column_type = type_dict[k] # gets a type name (as a string, e.g. "Float")
                if len(split_v) == 1:
                    # strict equality
                    val = self.do_type_cast(v, column_type)
                    try:
                        filters.append(self.table[k] == val)
                    except Exception as ex:
                        logger.error('Encountered exception!!')
                elif len(split_v) == 2:
                    val = self.do_type_cast(split_v[1], column_type)
                    try:
                        op = settings.OPERATOR_MAPPING[split_v[0]]
                    except KeyError as ex:
                        raise ParseException('The operator string ("{s}") was not understood. Choose'
                            ' from among: {vals}'.format(
                                s = split_v[0],
                                vals = ','.join(settings.OPERATOR_MAPPING.keys())
                            )
                        )
                    filters.append(
                        self.table[k].apply(lambda x: op(x, val))
                    )
                else:
                    raise ParseException('The query param string ({v}) for filtering on'
                        ' the {col} column was not formatted properly.'.format(
                            v = v,
                            col = k
                        )
                    )
            elif k in self.IGNORED_QUERY_PARAMS:
                 pass
            elif k == settings.ROWNAME_FILTER:
                if len(split_v) != 2:
                    raise ParseException('The query for filtering on the rows'
                        ' was not properly formatted. It should be [<op>]:<value>')
                # we don't allow indexes that are all numbers, so don't worry about casting
                # the filter value from a string
                val = split_v[1]
                try:
                    op = settings.OPERATOR_MAPPING[split_v[0]]
                except KeyError as ex:
                    raise ParseException('The operator string ("{s}") was not understood. Choose'
                        ' from among: {vals}'.format(
                            s = split_v[0],
                            vals = ','.join(settings.OPERATOR_MAPPING.keys())
                        )
                    )
                try:
                    rowname_filter = self.table.index.to_series().apply(lambda x: op(x, val))
                    filters.append(rowname_filter)
                except Exception as ex:
                    raise ParseException('Error encountered with filter on rows.'
                        ' Admin has been notified.'
                    )
                    alert_admins('Error when attempting to perform a row filter. Exception was: {x}'.format(x=ex))
            else:
                raise ParseException('The column "{c}" is not available for filtering.'.format(c=k))
        if len(filters) > 1:
            combined_filter = reduce(lambda x,y: x & y, filters)
            self.table = self.table.loc[combined_filter]
        elif len(filters) == 1:
            self.table = self.table.loc[filters[0]]

    def get_type_dict(self, **kwargs):
        '''
        Gets a mapping from the pandas native dtype to a MEV-compatible
        type.
        '''
        type_dict = {}
        for c in self.table.dtypes.index:
            # the convert_dtype function takes the native pandas dtype
            # and returns an attribute "type" that MEV understands.
            type_dict[c] = convert_dtype(str(self.table.dtypes[c]), **kwargs)
        return type_dict

    def replace_special_values(self):
        '''
        NaN and Inf values cause issues when we are serializing into JSON. If we have
        an action/event that requires serialization of the table-based resource, then we
        should use this method to safely convert those values

        If there is any filtering involved (e.g. filtering for values greater than x)
        ensure you don't call this method first, as it will replace infinity values with
        strings and the filter won't work properly.
        '''
        self.table = self.table.replace({
            -np.infty: settings.NEGATIVE_INF_MARKER,
            np.infty: settings.POSITIVE_INF_MARKER
        })
        self.table = self.table.mask(pd.isnull, None)

    def _resource_specific_modifications(self):
        '''
        This is a hook where derived classes can implement
        special filtering/behavior/etc. (if necessary)
        when the resource contents are requested.

        Called by the `get_contents` method.
        '''
        pass

    def main_contents_converter(self, row):
        '''
        A method that takes the resource contents to
        something that can be serialized.

        Input arg is a row of a pandas dataframe
        '''

        return {'rowname': row.name, 'values': row.to_dict()}

    def extra_contents_converter(self, row):
        '''
        A method that takes the resource contents to
        something that can be serialized.

        Input arg is a row of a pandas dataframe

        Note that this is different since we don't create a 'values' key
        as in the other contents converter method
        '''
        d = row.to_dict()
        d['rowname'] = row.name
        return d

    def get_contents(self, resource_path, query_params={}):
        '''
        Returns a dataframe of the table contents

        The dataframe allows the caller to subset as needed to 'paginate'
        the rows of the table
        '''

        try:
            logger.info('Read resource at {p}'.format(p=resource_path))
            self.read_resource(resource_path)
            self.additional_exported_cols = []

            # if there were any filtering params requested, apply those
            self.filter_against_query_params(query_params)
            self._resource_specific_modifications()
            self.perform_sorting(query_params)
            self.replace_special_values()
            if self.table.shape[0] > 0:
                standard_cols = [x for x in self.table.columns if not x in self.additional_exported_cols]
                content = self.table[standard_cols].apply(self.main_contents_converter, axis=1).tolist()
                if self.additional_exported_cols:
                    additional_content = self.table[self.additional_exported_cols].apply(self.extra_contents_converter, axis=1).tolist()
                    for x,y in zip(content, additional_content):
                        x.update(y)
                return content
            else:
                return []
        # for these first two exceptions, we already have logged
        # any problems when we called the `read_resource` method
        except ParserNotFoundException as ex:
            raise ex
        except ParseException as ex:
            raise ex
        # catch any other types of exceptions that we did not anticipate.
        except Exception as ex:
            logger.error('An unexpected error occurred when preparing'
                ' a resource preview for the resource at {path}. Exception'
                ' was: {ex}'.format(
                    path=resource_path,
                    ex=ex
                ))
            raise ex

    def extract_metadata(self, resource_path, parent_op_pk=None):
        '''
        This method extracts metadata from the Resource in question and
        saves it to the database.

        In the case of new Resources being added, the `parent_op` is None
        since no MEV-based analyses were responsible for the creation of the
        Resource.  If the Resource is created by some MEV-based analysis,
        the primary-key for that ExecutedOperation will be passed.

        '''
        logger.info('Extracting metadata from resource with path ({path}).'.format(
            path = resource_path
        ))

        # If the self.table field was not already filled, we need to
        # read the data
        if self.table is None:
            logger.info('Resource with path ({path}) was not '
                'previously parsed.  Do that now.'.format(
                    path=resource_path
                )
            )
            is_valid, message = self.validate_type(resource_path)
            if not is_valid:
                raise UnexpectedTypeValidationException(message)
        else:
            logger.info('Resource with path ({path}) was '
                'previously parsed.'.format(
                    path=resource_path
                )
            )        # now we have a table loaded at self.table.

        # call the super method to initialize the self.metadata
        # dictionary
        logger.info('Pre-populate/setup empty metadata.')
        super().setup_metadata()
        logger.info('Done with setting up the metadata.')
        # now add the information to self.metadata:
        if parent_op_pk:
            self.metadata[DataResource.PARENT_OP] = parent_op_pk

    def save_in_standardized_format(self, resource_path, resource_name):
        '''
        To avoid all the analyses having to concern themselves with data formats
        like csv, tsv, xlsx, etc. we just save table-based formats as a TSV.
        '''
        logger.info('Saving resource with path ({path}) to the standard format.'
            ' for a table-based resource'.format(
            path = resource_path
        ))

        # If the self.table field was not already filled, we need to
        # read the data
        if self.table is None:
            logger.info('Resource with path ({path}) was not '
                'previously parsed.  Do that now.'.format(
                    path=resource_path
                )
            )
            try:
                # this call will set the self.table member
                self.read_resource(resource_path)
            except Exception as ex:
                logger.error('Failed when trying to save in standard format.'
                    ' Specifically, failed when reading the resource, as self.table'
                    ' was None. Path of resource was: {path}'.format(
                        path=resource_path
                    )
                )
        # ok, self.table is set-- save it.
        current_file_extension = DataResource.get_extension(resource_path)
        file_dir =  os.path.dirname(resource_path)
        basename = os.path.basename(resource_path)
        basename_contents = basename.split('.')
        basename_contents[-1] = self.STANDARD_FORMAT
        new_basename = '.'.join(basename_contents)
        new_path = os.path.join(file_dir, new_basename)

        name_contents = resource_name.split('.')
        name_contents[-1] = self.STANDARD_FORMAT
        new_name = '.'.join(name_contents)

        logger.info('Writing the reformatted table-based resource to: {p}.'
            ' The new name is {n}.'.format(
            p = new_path,
            n = new_name
        ))
        self.table.to_csv(new_path, sep='\t')
        return new_path, new_name

class Matrix(TableResource):
    '''
    A `Matrix` is a delimited table-based file that has only numeric types.
    These types can be mixed, like floats and integers
    '''
    ACCEPTABLE_EXTENSIONS = [
        CSV,
        TSV,
        TAB,
        XLS,
        XLSX
    ]

    DESCRIPTION = 'A table of where all the entries are numbers'\
        ' except the first column (which names the rows) and the' \
        ' first line (which gives the column names). The cell at the' \
        ' first row and column may be left blank. Acceptable file extensions' \
        ' include: {s}.'.format(s = ', '.join(ACCEPTABLE_EXTENSIONS))

    EXAMPLE = [
        {
            "rowname":"gene1",
            "values": {
                "sampleA": 1.1,
                "sampleB": 2.2,
                "sampleC": 3.3,
            }
        },
        {
            "rowname":"gene2",
            "values": {
                "sampleA": 11.1,
                "sampleB": 22.2,
                "sampleC": 33.3,
            }
        },
        {
            "rowname":"gene3",
            "values": {
                "sampleA": 111.1,
                "sampleB": 222.2,
                "sampleC": 333.3,
            }
        }
    ]

    # looking for integers OR floats.  Both are acceptable
    TARGET_PATTERN = '(float|int)\d{0,2}'

    # Define some additional filtering/sorting params that
    # are specific to this type and its children.
    # A special filter for allowing filtering of numeric tables based on
    # the row means
    ROWMEAN_KEYWORD = '__rowmean__'
    INCLUDE_ROWMEANS = '__incl_rowmeans__'
    EXTRA_MATRIX_QUERY_PARAMS = [
        ROWMEAN_KEYWORD,
        INCLUDE_ROWMEANS
    ]

    # Copy the ignored params from the parent (don't want to modify that)
    IGNORED_QUERY_PARAMS = [x for x in TableResource.IGNORED_QUERY_PARAMS]
    IGNORED_QUERY_PARAMS.extend(EXTRA_MATRIX_QUERY_PARAMS)

    def check_column_types(self, target_pattern):
        '''
        Checks each column against a specific numpy/pandas dtype.
        The specific dtype comes from the class member.
        '''
        problem_columns = []
        for i,col in enumerate(self.table.dtypes):
            if not re.match(target_pattern, str(col)):
                colname = self.table.columns[i]
                problem_columns.append(
                    (colname, i+1)
                )
        return problem_columns

    def validate_type(self, resource_path):
        is_valid, error_msg = super().validate_type(resource_path)
        if not is_valid:
            return (False, error_msg)

        # was able to at least open/parse the file.
        # now check for numeric types
        problem_columns = self.check_column_types(Matrix.TARGET_PATTERN)

        if len(problem_columns) > 0:
            error_message = NON_NUMERIC_ERROR
            return (False, error_message)

        return (True, None)

    def extract_metadata(self, resource_path, parent_op_pk=None):

        super().extract_metadata(resource_path, parent_op_pk)

        # Note: removed the addition of FeatureSets to the metadata as it was causing
        # issues with large json objects being inserted into the database.
        # the FeatureSet comes from the rows:
        # f_set = FeatureSet([Feature(x) for x in self.table.index])
        # self.metadata[DataResource.FEATURE_SET] = FeatureSetSerializer(f_set).data

        # the ObservationSet comes from the cols:
        o_set = ObservationSet([Observation(x) for x in self.table.columns])
        self.metadata[DataResource.OBSERVATION_SET] = ObservationSetSerializer(o_set).data
        return self.metadata

    def _resource_specific_modifications(self):
        if not self.extra_query_params:
            return

        filters = []
        if self.ROWMEAN_KEYWORD in self.extra_query_params:
            self.additional_exported_cols.append(self.ROWMEAN_KEYWORD)
            self.table[self.ROWMEAN_KEYWORD] = self.table.mean(axis=1)
            filter_string = self.extra_query_params[self.ROWMEAN_KEYWORD]
            split_str = filter_string.split(settings.QUERY_PARAM_DELIMITER)
            if len(split_str) == 1:
                # strict equality
                try:
                    val = float(split_str[0])
                except ValueError as ex:
                    raise ParseException('Could not parse the request'
                        ' for filtering on the row means. The value'
                        ' could not be interpreted as a number.'
                    )

                filters.append(self.table[self.ROWMEAN_KEYWORD] == val)

            elif len(split_str) == 2:
                try:
                    val = float(split_str[1])
                except ValueError as ex:
                    raise ParseException('Could not parse the request'
                        ' for filtering on the row means. The value'
                        ' could not be interpreted as a number.'
                    )
                try:
                    op = settings.OPERATOR_MAPPING[split_str[0]]
                except KeyError as ex:
                    raise ParseException('The operator string ("{s}") was not understood. Choose'
                        ' from among: {vals}'.format(
                            s = split_str[0],
                            vals = ','.join(settings.OPERATOR_MAPPING.keys())
                        )
                    )
                filters.append(self.table[self.ROWMEAN_KEYWORD].apply(lambda x: op(x, val)))
            else:
                raise ParseException('The query param string ({v}) for filtering on'
                    ' the mean values was not formatted properly.'.format(
                        v = filter_string
                    )
                )

        if self.INCLUDE_ROWMEANS in self.extra_query_params:
            # if the rowmean filter was also applied, we can skip--
            # This query param is requesting the row means, which is inferred by the rowmean filter
            if not self.ROWMEAN_KEYWORD in self.extra_query_params:
                # if here, no filtering was requested on the row means. Simply calculate.
                # Note that we use the rowmean_keyword as the column name and in the
                # export columns list so that the final converter knows to send this column
                # in the response.
                self.additional_exported_cols.append(self.ROWMEAN_KEYWORD)
                self.table[self.ROWMEAN_KEYWORD] = self.table.mean(axis=1)

        # apply filters (if any)
        if len(filters) > 1:
            combined_filter = reduce(lambda x,y: x & y, filters)
            self.table = self.table.loc[combined_filter]
        elif len(filters) == 1:
            self.table = self.table.loc[filters[0]]

    def get_contents(self, resource_path, query_params={}):
        '''
        This method allows us to add on additional content that is
        allowable for matrix types, as they are all numeric.
        Examples include providing the row means, etc.
        '''
        # pull out any matrix-specific query params so they don't interfere
        # with the general query params passed to the parent.
        # If the request doesn't contain any of these params, then
        # the self.current_query_params dict is empty and no further
        # modifications to the returned content are made.
        self.extra_query_params = {}
        for p in self.EXTRA_MATRIX_QUERY_PARAMS:
            if p in query_params:
                self.extra_query_params[p] = query_params[p]

        # additional filtering/behavior specific to a Matrix (if requested)
        # is handled in the _resource_specific_modifications method
        return super().get_contents(resource_path, query_params)

class IntegerMatrix(Matrix):
    '''
    An `IntegerMatrix` further specializes the `Matrix`
    to admit only integers.
    '''
    # looking for only integers.
    TARGET_PATTERN = 'int\d{0,2}'

    DESCRIPTION = 'A table where all the entries are integers'\
        ' except the first column (which names the rows) and the' \
        ' first line (which gives the column names). The cell at the' \
        ' first row and column may be left blank. Acceptable extensions' \
        ' include: {s}'.format(s = ', '.join(Matrix.ACCEPTABLE_EXTENSIONS))

    EXAMPLE = [
        {
            "rowname":"gene1",
            "values": {
                "sampleA": 1,
                "sampleB": 2,
                "sampleC": 3,
            }
        },
        {
            "rowname":"gene2",
            "values": {
                "sampleA": 11,
                "sampleB": 22,
                "sampleC": 33,
            }
        },
        {
            "rowname":"gene3",
            "values": {
                "sampleA": 111,
                "sampleB": 222,
                "sampleC": 333,
            }
        }
    ]

    def validate_type(self, resource_path):

        # first check that it has all numeric types.  If that fails
        # immediately return--
        is_valid, error_message = super().validate_type(resource_path)
        if not is_valid:
            return (False, error_message)

        # was valid for numeric types.  Now check for integer
        problem_columns = self.check_column_types(IntegerMatrix.TARGET_PATTERN)
        if len(problem_columns) > 0:

            # one problem with pandas is that NaN values cause a column
            # to be parsed as a float, even if all other values in the
            # column are integers.  We can do a secondary check, however,
            # to see if the remaining values (non-NaN) are basically
            # integers.  We check if the numbers look like "2.0".
            # If that is the case, we remove that column from the
            # "problem columns".
            for i,c in enumerate(problem_columns):
                # recall c is a tuple of (colname, col number)
                if all([
                        re.match('\d+\.0', str(x))
                        for x in self.table[c[0]].dropna()]):
                    problem_columns.pop(i)

            # if there are still any remaining problematic cols,
            # we now issue an error
            if len(problem_columns) > 0:
                error_message = NON_INTEGER_ERROR
                return (False, error_message)

        return (True, None)

class RnaSeqCountMatrix(IntegerMatrix):
    '''
    A very-explicit class (mainly for making things user-friendly)
    where we provide specialized behavior/messages specific to count matrices
    generated from RNA-seq data. The same as an integer matrix, but named to be
    suggestive for WebMEV users.
    '''
    DESCRIPTION = 'A table of integer-based counts corresponding to'\
        ' the number of sequencing reads associated with a particular' \
        ' gene or transcript. Acceptable file extensions' \
        ' include: {s}'.format(s = ', '.join(Matrix.ACCEPTABLE_EXTENSIONS))

    EXAMPLE = IntegerMatrix.EXAMPLE

class ElementTable(TableResource):
    '''
    An ElementTable captures common behavior of tables which
    annotate `Observation`s (AnnotationTable) or `Feature`s (FeatureTable)

    It's effectively an abstract class. See the derived classes which implement
    the specific behavior for `Observation`s or `Feature`s.
    '''

    ACCEPTABLE_EXTENSIONS = [
        CSV,
        TSV,
        TAB,
        XLS,
        XLSX
    ]

    def validate_type(self, resource_path):

        # check that file can be parsed:
        is_valid, error_message = super().validate_type(resource_path)
        if not is_valid:
            return (False, error_message)

        # check that the file is "useful" in that it has
        # more than one column.  It's not REALLY an error, but it does not
        # provide any information.  This can also be caught earlier, but
        # we provide it here just as a secondary guard.
        if self.table.shape[1] == 0:
            return (False, TRIVIAL_TABLE_ERROR)
        return (True, None)

    def prep_metadata(self, element_class):
        '''
        When we extract the metadata from an ElementTable, we
        expect the Element instances (Observations or Features)
        to be contained in the rows.

        Additional columns specify attributes which we incorporate.

        The `element_class` arg is a class which implements the specific
        type we want (i.e. Observation or Feature)
        '''
        logger.info('For element type {s},'
            ' extract out metadata'.format(s=type(element_class)))
        # Go through the columns and find out the primitive types
        # for each column/covariate.
        # Note that we can't determine specific types (e.g. bounded integers)
        # from general annotations.  We basically allow floats, integers, and
        # "other" types, which get converted to strings.
        # Also note that we pass the `allow_unrestricted_strings` kwarg
        # which will allow the attribute values to take on arbitrary string
        # values.
        type_dict = self.get_type_dict(allow_unrestricted_strings=True)

        # convert NaN and infs to our special marker values
        self.replace_special_values()

        element_list = []
        for id, row_series in self.table.iterrows():
            d = row_series.to_dict()
            attr_dict = {}
            for key, val in d.items():
                # Note the 'allow_null=True', so that attributes can be properly serialized
                # if they are missing a value. This happens, for instance, in FeatureTable
                # instances where p-values were not assigned.
                attr = create_attribute(key,
                    {
                        'attribute_type': type_dict[key],
                        'value': val
                    },
                    allow_null=True
                )
                attr_dict[key] = attr
            element_list.append(element_class(id, attr_dict))
        return element_list

class AnnotationTable(ElementTable):
    '''
    An `AnnotationTable` is a special type of table that will be responsible
    for annotating Observations/samples (e.g. adding sample names and
    associated attributes like experimental group or other covariates)

    The first column will give the sample names and the remaining columns will
    each individually represent different covariates associated with that sample.
    '''

    DESCRIPTION = 'This type of file is used to add metadata to your samples.' \
        ' The first column has the sample name and the remaining columns contain' \
        ' metadata about each sample (for instance, experimental group,'\
        ' treatment, or similar. Acceptable file extensions' \
        ' include: {s}.'.format(s = ', '.join(ElementTable.ACCEPTABLE_EXTENSIONS))

    EXAMPLE = [
        {
            "rowname":"SampleA",
            "values": {
                "cell_type": 'CD4',
                "treatment": 'Y'
            }
        },
        {
            "rowname":"SampleB",
            "values": {
                "cell_type": 'CD8',
                "treatment": 'Y'
            }
        },
        {
            "rowname":"SampleC",
            "values": {
                "cell_type": 'Monocyte',
                "treatment": 'N'
            }
        }
    ]

    def validate_type(self, resource_path):

        # check that file can be parsed:
        is_valid, error_message = super().validate_type(resource_path)
        if not is_valid:
            return (False, error_message)

        # it is hard to check for proper headers for annotation
        # files since they have relatively free format.  However,
        # if the column name matches any values in its column, the
        # annotation was likely missing a header line.  For example,
        # >>> df
        # SW1_Control   CTRL
        # 0  SW2_Control   CTRL
        # 1  SW3_Control   CTRL
        # 2  SW4_Treated  TREAT
        # 3  SW5_Treated  TREAT
        # 4  SW6_Treated  TREAT
        #
        # Here, "CTRL" becomes the header, but it's clearly
        # just due to a missing header.  We don't issue an
        # error, but we do warn the user by adding a comment.
        flagged_columns = []
        for c in self.table.columns:
            if np.sum(c == self.table[c]) > 0:
                flagged_columns.append(c)
        if len(flagged_columns) > 0:
            return (True, MISSING_HEADER_WARNING)

        return (True, None)

    def extract_metadata(self, resource_path, parent_op_pk=None):
        '''
        When we extract the metadata from an AnnotationTable, we
        expect the Observation instances to be the rows.

        Additional columns specify attributes of each Observation,
        which we incorporate
        '''
        super().extract_metadata(resource_path, parent_op_pk)

        observation_list = super().prep_metadata(Observation)
        o_set = ObservationSet(observation_list)
        self.metadata[DataResource.OBSERVATION_SET] = ObservationSetSerializer(o_set).data
        return self.metadata

class FeatureTable(ElementTable):
    '''
    A `FeatureTable` is a type of table that has aggregate information about
    the features, but does not have any "observations" in the columns.  An example
    would be the results of a differential expression analysis.  Each row corresponds
    to a gene (feature) and the columns are information about that gene (such as p-value).

    Another example could be a table of metadata about genes (e.g. pathways or perhaps a
    mapping to a different gene identifier).

    The first column will give the feature/gene identifiers and the remaining columns will
    have information about that gene
    '''

    DESCRIPTION = 'This type of file describes the "features" of your data.  In the genomics' \
        ' context, one example of a feature is a gene.  Therefore, you could use this table' \
        ' to give additional information about each gene, such as alternative symbols,' \
        ' oncogene status, or similar.  Each row contains information about a single gene.' \
        ' Note, however, that this concept is completely general and not restricted' \
        ' to information about genes or transcripts. Acceptable file extensions' \
        ' include: {s}'.format(s = ', '.join(ElementTable.ACCEPTABLE_EXTENSIONS))

    EXAMPLE = [
        {
            "rowname":"gene1",
            "values": {
                "logFoldChange": 1.1,
                "pvalue": 0.03,
                "FDR": 0.2,
            }
        },
        {
            "rowname":"gene2",
            "values": {
                "logFoldChange": 2.2,
                "pvalue": 0.01,
                "FDR": 0.03,
            }
        },
        {
            "rowname":"gene3",
            "values": {
                "logFoldChange": 3.3,
                "pvalue": 0.000003,
                "FDR": 0.001,
            }
        }
    ]

    def extract_metadata(self, resource_path, parent_op_pk=None):
        '''
        When we extract the metadata from a FeatureTable, we
        expect the Feature instances to be the rows.

        Additional columns specify attributes of each Feature,
        which we incorporate
        '''
        logger.info('Extract metadata from a FeatureTable')
        super().extract_metadata(resource_path, parent_op_pk)

        # Note: removed the addition of FeatureSets to the metadata as it was causing
        # issues with large json objects being inserted into the database.
        # feature_list = super().prep_metadata(Feature)
        # f_set = FeatureSet(feature_list)
        # self.metadata[DataResource.FEATURE_SET] = FeatureSetSerializer(f_set).data
        return self.metadata

class BEDFile(TableResource):
    '''
    A file format that corresponds to the BED format.  This is
    the minimal BED format, which has:

    - chromosome
    - start position
    - end position

    Additional columns are ignored.

    By default, BED files do NOT contain headers and we enforce that here.
    '''

    ACCEPTABLE_EXTENSIONS = [BED,]

    DESCRIPTION = 'A three-column BED-format file. https://ensembl.org/info/website/upload/bed.html'\
        ' BED files do NOT have column headers.' \
        ' Acceptable file extensions include: {s}'.format(s=', '.join(ACCEPTABLE_EXTENSIONS))

    def validate_type(self, resource_path):
        reader = TableResource.get_reader(resource_path)

        # if the BED file has a header, the reader below will incorporate
        # that into the columns and the 2nd and 3rd columns will no longer have
        # the proper integer type.
        table = reader(resource_path,
            names=['chrom','start','stop'],
            usecols=[0,1,2])
        start_col_int = re.match('int\d{0,2}', str(table['start'].dtype))
        stop_col_int = re.match('int\d{0,2}', str(table['stop'].dtype))
        if start_col_int and stop_col_int:
            return (True, None)
        else:
            problem_columns = []
            if start_col_int is None:
                problem_columns.append(2)
            if stop_col_int is None:
                problem_columns.append(3)

            cols = ','.join([str(x) for x in problem_columns])
            error_message = BED_FORMAT_ERROR.format(cols=cols)
            return (False, error_message)

    def extract_metadata(self, resource_path, parent_op_pk=None):
        super().extract_metadata(resource_path, parent_op_pk)
        return self.metadata

max_stars_repo_name: hsph-qbrc/mev-backend

max_issues_repo_name: hsph-qbrc/mev-backend

max_forks_repo_name: hsph-qbrc/mev-backend


text:
"""
With this class new feedback can be added to the database.
"""
import os
from Product.Database.DatabaseManager.Insert.InsertFeedback import InsertFeedback
from Product.RecommendationManager.model import generate_model as generate_model
from Product.Database.DatabaseManager.Retrieve.RetrieveMovie import RetrieveMovie
from Product.Database.DatabaseManager.Retrieve.RetrieveUser import RetrieveUser
import scipy.sparse as sp

class Feedback(object):
    """
    Author: Alexander Dahl, Marten Bolin
    Date: 2017-11-17
    Last update:
    Purpose: Feedback class that can be used to create feedback to the database
    """
    @staticmethod
    def insert_feedback(user_id, movie_id, watched=None, rating=None):
        """
        Author: Alexander Dahl, Marten Bolin
        Date: 2017-11-17
        Last update:
        Purpose: Make Insert of Feedback to the database.
        if a feedback has a rating this will be used to evolve the lightFM model.
        :param user_id : The id of the user that has made the rating or watched
        :type int
        :param movie_id : The id of the movie that has been rated or watched
        :type int
        :param watched : 1 if has watched (optional)
        :type int
        :param rating : the rating that was made, 1-5 (optional)
        :type float
        """
        InsertFeedback().insert_feedback(user_id, movie_id, watched, rating)

        if rating:
            model = generate_model.load_model(generate_model.get_path())
            # Converting to lists because otherwise it is not possible to convert them to
            # coo matrices
            rating_list = [rating]
            user_list = [user_id]
            movie_list = [movie_id]

            # Need to define shape so that it is coherent with the previous model
            # TODO Change dimensions to greater than 1 if problem with dimensions
            user_matrix = sp.coo_matrix((rating_list, (user_list, movie_list)),
                                        shape=(RetrieveUser().retrieve_largest_user_id()+1,
                                               RetrieveMovie().retrieve_largest_movie_id()+1))
            generate_model.evolve_model(generate_model.get_path(), model, user_matrix)

ext: py

max_stars_repo_name: VincentDehaye/recommender-system-liu

max_issues_repo_name: VincentDehaye/recommender-system-liu

max_forks_repo_name: VincentDehaye/recommender-system-liu

reason: import scipy


text:
import numpy as np
import torch

# import matplotlib.pyplot as plt

IIWA_JOINT_MIN_LIMITS = torch.deg2rad(torch.tensor([-170., -120., -170., -120., -170., -120., -175.]))
IIWA_JOINT_MAX_LIMITS = torch.deg2rad(torch.tensor([170., 120., 170., 120., 170., 120., 175.]))

class Rotation:
    def __init__(self, quat, normalize=True):
        self._quat = quat if isinstance(quat, torch.Tensor) else torch.tensor(quat).double()
        if normalize:
            self._quat /= torch.norm(quat)

    @classmethod
    def from_quat(cls, quat):
        return cls(quat, normalize=True)

    @classmethod
    def from_matrix(cls, matrix):
        trace = torch.trace(matrix)
        quat = torch.zeros(4).double()
        if trace > 1e-5:
            s = torch.sqrt(trace + 1) * 2
            quat[0] = s / 4
            quat[1] = (matrix[2, 1] - matrix[1, 2]) / s
            quat[2] = (matrix[0, 2] - matrix[2, 0]) / s
            quat[3] = (matrix[1, 0] - matrix[0, 1]) / s
        elif (matrix[0, 0] > matrix[1, 1]) and (matrix[0, 0] > matrix[2, 2]):
            s = torch.sqrt(1.0 + matrix[0, 0] - matrix[1, 1] - matrix[2, 2]) * 2
            quat[0] = (matrix[2, 1] - matrix[1, 2]) / s
            quat[1] = 0.25 * s
            quat[2] = (matrix[0, 1] + matrix[1, 0]) / s
            quat[3] = (matrix[0, 2] + matrix[2, 0]) / s
        elif (matrix[1, 1] > matrix[2, 2]):
            s = torch.sqrt(1.0 + matrix[1, 1] - matrix[0, 0] - matrix[2, 2]) * 2
            quat[0] = (matrix[0, 2] - matrix[2, 0]) / s
            quat[1] = (matrix[0, 1] + matrix[1, 0]) / s
            quat[2] = 0.25 * s
            quat[3] = (matrix[1, 2] + matrix[2, 1]) / s
        else:
            s = torch.sqrt(1.0 + matrix[2, 2] - matrix[0, 0] - matrix[1, 1]) * 2
            quat[0] = (matrix[1, 0] - matrix[0, 1]) / s
            quat[1] = (matrix[0, 2] + matrix[2, 0]) / s
            quat[2] = (matrix[1, 2] + matrix[2, 1]) / s
            quat[3] = 0.25 * s
        return cls(quat, normalize=False)

    @classmethod
    def from_rotvec(cls, rotvec):
        quat = torch.zeros(4).double()
        norm = torch.norm(rotvec)
        small_angle = (norm <= 1e-3)

        if small_angle:
            scale = (0.5 - norm ** 2 / 48 + norm ** 4 / 3840)
        else:
            scale = (torch.sin(norm / 2) / norm)

        quat[0] = torch.cos(norm / 2)
        quat[1:] = scale * rotvec
        return cls(quat, normalize=True)

    @classmethod
    def from_euler_yzy(cls, yzy):
        a = yzy[0]
        b = yzy[1]
        c = yzy[2]
        rot_mat = torch.tensor(
            [[-torch.sin(a) * torch.sin(c) + torch.cos(a) * torch.cos(b) * torch.cos(c),
              -torch.sin(b) * torch.cos(c),
              torch.sin(a) * torch.cos(b) * torch.cos(c) + torch.sin(c) * torch.cos(a)],
             [torch.sin(b) * torch.cos(a),
              torch.cos(b),
              torch.sin(a) * torch.sin(b)],
             [-torch.sin(a) * torch.cos(c) - torch.sin(c) * torch.cos(a) * torch.cos(b),
              torch.sin(b) * torch.sin(c),
              -torch.sin(a) * torch.sin(c) * torch.cos(b) + torch.cos(a) * torch.cos(c)]])
        return cls.from_matrix(rot_mat)

    def as_quat(self):
        return self._quat

    def as_matrix(self):
        quat = self._quat
        R = torch.eye(3).double()
        R[0, 0] = 1 - 2 * (quat[2] ** 2 + quat[3] ** 2)
        R[0, 1] = 2 * (quat[1] * quat[2] - quat[3] * quat[0])
        R[0, 2] = 2 * (quat[1] * quat[3] + quat[2] * quat[0])

        R[1, 0] = 2 * (quat[1] * quat[2] + quat[3] * quat[0])
        R[1, 1] = 1 - 2 * (quat[1] ** 2 + quat[3] ** 2)
        R[1, 2] = 2 * (quat[2] * quat[3] - quat[1] * quat[0])

        R[2, 0] = 2 * (quat[1] * quat[3] - quat[2] * quat[0])
        R[2, 1] = 2 * (quat[2] * quat[3] + quat[1] * quat[0])
        R[2, 2] = 1 - 2 * (quat[1] ** 2 + quat[2] ** 2)
        return R

    def as_rotvec(self):
        angle = 2 * np.arctan2(np.linalg.norm(self._quat[1:]), self._quat[0])

        small_angle = angle <= 1e-3

        if small_angle:
            scale = (2 + angle[small_angle] ** 2 / 12 + 7 * angle[small_angle] ** 4 / 2880)
        else:
            scale = (angle / np.sin(angle / 2))

        rotvec = scale * self._quat[1:]
        return rotvec

    def as_euler_yzy(self):
        rot_mat = self.as_matrix()
        yzy = torch.zeros(3)
        if torch.isclose(torch.abs(rot_mat[1, 1]), torch.ones(1).double()):
            yzy[0] = torch.atan2(rot_mat[0, 2], rot_mat[0, 0])
            yzy[1] = torch.acos(rot_mat[1, 1])
            yzy[2] = 0
        else:
            yzy[0] = torch.atan2(rot_mat[1, 2], rot_mat[1, 0])
            yzy[1] = torch.atan2(torch.sqrt(rot_mat[1, 2] ** 2 + rot_mat[1, 0] ** 2), rot_mat[1, 1])
            yzy[2] = torch.atan2(rot_mat[2, 1], -rot_mat[0, 1])
        return yzy

class KinematicsTorch:
    def __init__(self, tcp_pos=None, tcp_quat=None):
        self.d_bs = torch.tensor(0.36)
        self.d_se = torch.tensor(0.42)
        self.d_ew = torch.tensor(0.4)
        self.d_wf = torch.tensor(0.151)

        # DH_paramters a_i, alpha_i, d_i
        self.dh_a = torch.tensor([0., 0., 0., 0., 0., 0., 0.])
        self.dh_alpha = torch.tensor([-np.pi / 2, np.pi / 2, np.pi / 2, -np.pi / 2, -np.pi / 2, np.pi / 2, 0.])
        self.dh_d = torch.tensor([self.d_bs, 0., self.d_se, 0., self.d_ew, 0., self.d_wf])

        self.joint_limits = torch.stack((IIWA_JOINT_MIN_LIMITS, IIWA_JOINT_MAX_LIMITS), dim=0).T

        self.joint_vel_limits = torch.deg2rad(torch.tensor([[-85., 85.],
                                                            [-85., 85.],
                                                            [-100., 100.],
                                                            [-75., 75.],
                                                            [-130., 130.],
                                                            [-135., 135.],
                                                            [-135., 135.]]))

        self.singularity_eps = 0.1

        self.T_ee = torch.eye(4).double()
        if tcp_pos is not None:
            self.T_ee[:3, 3] = tcp_pos
        if tcp_quat is not None:
            self.T_ee[:3, :3] = Rotation.from_quat(tcp_quat).as_matrix()

    def _transform_i(self, index, theta_i):
        alpha_i = self.dh_alpha[index]
        a_i = self.dh_a[index]
        d_i = self.dh_d[index]
        T_i = torch.eye(4).double()
        T_i[0, 0] = torch.cos(theta_i)
        T_i[0, 1] = -torch.sin(theta_i) * torch.cos(alpha_i)
        T_i[0, 2] = torch.sin(theta_i) * torch.sin(alpha_i)
        T_i[0, 3] = a_i * torch.cos(theta_i)
        T_i[1, 0] = torch.sin(theta_i)
        T_i[1, 1] = torch.cos(theta_i) * torch.cos(alpha_i)
        T_i[1, 2] = -torch.cos(theta_i) * torch.sin(alpha_i)
        T_i[1, 3] = a_i * torch.sin(theta_i)
        T_i[2, 0] = 0.
        T_i[2, 1] = torch.sin(alpha_i)
        T_i[2, 2] = torch.cos(alpha_i)
        T_i[2, 3] = d_i
        return T_i

    def _transform(self, q):
        T = torch.eye(4).double()
        for i, theta_i in enumerate(q):
            T_i = self._transform_i(i, theta_i)
            T = T @ T_i
        T = T @ self.T_ee
        return T

    def forward_kinematics(self, q):
        T = self._transform(q)
        position = T[:-1, 3]
        rot_mat = T[:-1, :-1]
        return torch.cat([position, Rotation.from_matrix(rot_mat).as_quat()], dim=0)

    def get_jacobian(self, q):
        q.requires_grad_(True)
        y = self.forward_kinematics(q)
        jac = torch.zeros(y.shape[0], q.shape[0]).double()
        for i in range(q.shape[0]):
            jac[i] = torch.autograd.grad(y[i], q, create_graph=True, retain_graph=True)[0]
        return jac

    def get_redundancy(self, q):
        pose = self.forward_kinematics(q)
        gc_4 = torch.sign(q[3])
        p = pose[:3]
        rot_mat = Rotation.from_quat(pose[3:]).as_matrix()

        rot_0_e = rot_mat @ self.T_ee[:3, :3].T
        p_0_e = p - rot_0_e @ self.T_ee[:3, 3]

        p_0_2 = torch.tensor([0., 0., self.d_bs]).double()
        p_6_7 = torch.tensor([0., 0., self.d_wf]).double()

        # get q4 and p2_6 Shoulder-Elbow-Wrist (SEW) Plane
        p_2_6 = p_0_e - p_0_2 - rot_0_e @ p_6_7

        cosq_4_v = (torch.norm(p_2_6, dim=0) ** 2 - self.d_se ** 2 - self.d_ew ** 2) / (2 * self.d_se * self.d_ew)
        q_4_v = gc_4 * torch.acos(cosq_4_v)

        R_0_1_z = torch.tensor([0., 0., 1.]).double()
        if torch.norm(torch.cross(p_2_6, R_0_1_z)) < 1e-3:
            q_1_v = 0
        else:
            q_1_v = torch.atan2(p_2_6[1], p_2_6[0])

        phi = torch.acos(
            (self.d_se ** 2 + torch.norm(p_2_6) ** 2 - self.d_ew ** 2) / (2 * self.d_se * torch.norm(p_2_6)))

        q_2_v = torch.atan2(torch.sqrt(p_2_6[0] ** 2 + p_2_6[1] ** 2), p_2_6[2]) + gc_4 * phi

        q_3_v = torch.tensor(0.).double()

        p_0_2_v = (self._transform_i(0, q_1_v) @ self._transform_i(1, q_2_v))[:3, 3]
        p_0_4_v = (self._transform_i(0, q_1_v) @ self._transform_i(1, q_2_v) @
                   self._transform_i(2, q_3_v) @ self._transform_i(3, q_4_v))[:3, 3]

        p_0_4 = (self._transform_i(0, q[0]) @ self._transform_i(1, q[1]) @
                 self._transform_i(2, q[2]) @ self._transform_i(3, q[3]))[:3, 3]

        vec_2_4_v = (p_0_4_v - p_0_2_v) / torch.norm(p_0_4_v - p_0_2_v)
        vec_2_6_v = (p_2_6) / torch.norm(p_2_6)
        vec_2_4 = (p_0_4 - p_0_2) / torch.norm(p_0_4 - p_0_2)
        vec_2_6 = (p_2_6) / torch.norm(p_2_6)

        vec_sew_v = torch.cross(vec_2_4_v, vec_2_6_v)
        vec_sew = torch.cross(vec_2_4, vec_2_6)
        sign_psi = torch.sign(torch.cross(vec_sew_v, vec_sew).dot(p_2_6))
        psi = sign_psi * torch.acos(vec_sew_v.dot(vec_sew) / (torch.norm(vec_sew) * torch.norm(vec_sew_v)))
        return psi

    def _get_auxiliary_parameter(self, p, rot_mat, gc):
        p0_2 = torch.tensor([0., 0., self.d_bs]).double()
        p6_7 = torch.tensor([0., 0., self.d_wf]).double()
        R0_7 = rot_mat

        self.gc = np.array([gc[0], gc[0], gc[0], gc[1], gc[2], gc[2], gc[2]])

        # get q4 and p2_6 Shoulder-Elbow-Wrist (SEW) Plane
        p2_6 = p - p0_2 - rot_mat @ p6_7
        cosq4_v = (torch.norm(p2_6, dim=0) ** 2 - self.d_se ** 2 - self.d_ew ** 2) / (2 * self.d_se * self.d_ew)
        if abs(cosq4_v) > 1 + 1e-9:
            self.q4_v = None
            self.a = torch.zeros((7, 2))
            self.b = torch.zeros((7, 2))
            self.c = torch.zeros((7, 2))
            return False
        else:
            cosq4_v = torch.clamp(cosq4_v, -1., 1.)
            self.q4_v = gc[1] * torch.acos(cosq4_v)
        p2_6_norm = torch.norm(p2_6)

        q1_v = torch.atan2(p2_6[1], p2_6[0])
        phi = torch.acos((self.d_se ** 2 + p2_6_norm ** 2 - self.d_ew ** 2) / (2 * self.d_se * p2_6_norm))
        q2_v = torch.atan2(torch.sqrt(p2_6[0] ** 2 + p2_6[1] ** 2), p2_6[2]) + gc[1] * phi
        q3_v = torch.tensor(0.0)
        T_0_3_v = self._transform_i(0, q1_v) @ self._transform_i(1, q2_v) @ self._transform_i(2, q3_v)
        R_0_3_v = T_0_3_v[:3, :3]

        # cross product matrixq_true
        p2_6_hat = p2_6 / p2_6_norm
        cpm = torch.tensor([[0., -p2_6_hat[2], p2_6_hat[1]],
                            [p2_6_hat[2], 0., -p2_6_hat[0]],
                            [-p2_6_hat[1], p2_6_hat[0], 0.]])
        A_s = cpm @ R_0_3_v
        B_s = - (cpm @ cpm) @ R_0_3_v
        C_s = torch.ger(p2_6_hat, p2_6_hat) @ R_0_3_v

        T_3_4 = self._transform_i(3, self.q4_v)
        R_3_4 = T_3_4[:3, :3]

        A_w = R_3_4.T @ A_s.T @ R0_7
        B_w = R_3_4.T @ B_s.T @ R0_7
        C_w = R_3_4.T @ C_s.T @ R0_7

        # x[3, :]=0 For joint 4 is zero
        # x[1 5, 1]=0 a_d, b_d, c_d =0
        self.a = torch.zeros((7, 2))
        self.b = torch.zeros((7, 2))
        self.c = torch.zeros((7, 2))

        self.a[0, 0] = A_s[1, 1]
        self.b[0, 0] = B_s[1, 1]
        self.c[0, 0] = C_s[1, 1]
        self.a[0, 1] = A_s[0, 1]
        self.b[0, 1] = B_s[0, 1]
        self.c[0, 1] = C_s[0, 1]

        self.a[1, 0] = A_s[2, 1]
        self.b[1, 0] = B_s[2, 1]
        self.c[1, 0] = C_s[2, 1]

        self.a[2, 0] = -A_s[2, 2]
        self.b[2, 0] = -B_s[2, 2]
        self.c[2, 0] = -C_s[2, 2]
        self.a[2, 1] = -A_s[2, 0]
        self.b[2, 1] = -B_s[2, 0]
        self.c[2, 1] = -C_s[2, 0]

        self.a[4, 0] = A_w[1, 2]
        self.a[4, 1] = A_w[0, 2]
        self.b[4, 0] = B_w[1, 2]
        self.b[4, 1] = B_w[0, 2]
        self.c[4, 0] = C_w[1, 2]
        self.c[4, 1] = C_w[0, 2]

        self.a[5, 0] = A_w[2, 2]
        self.b[5, 0] = B_w[2, 2]
        self.c[5, 0] = C_w[2, 2]

        self.a[6, 0] = A_w[2, 1]
        self.a[6, 1] = -A_w[2, 0]
        self.b[6, 0] = B_w[2, 1]
        self.b[6, 1] = -B_w[2, 0]
        self.c[6, 0] = C_w[2, 1]
        self.c[6, 1] = -C_w[2, 0]
        return True

    def _transform_tcp2ee(self, pose):
        # R_0_tcp = R_0_e @ R_e_tcp -> R_0_e = R_0_tcp @ R_e_tcp.T
        rot_0_e = Rotation.from_quat(pose[3:]).as_matrix() @ self.T_ee[:3, :3].T
        # p_0_tcp = p_0_e + R_0_e * p_e_tcp
        p_0_e = pose[:3] - rot_0_e @ self.T_ee[:3, 3]
        return p_0_e, rot_0_e

    def inverse_kinematics(self, pose, psi, gc):
        q = torch.zeros((7,))
        # get auxiliary parameter a, b, c: [7, 2],

        p_0_e, rot_0_e = self._transform_tcp2ee(pose)

        ret = self._get_auxiliary_parameter(p_0_e, rot_0_e, gc)

        q[:] = np.nan
        if not ret:
            return False, q

        cos_q2 = self.a[1, 0] * torch.sin(psi) + self.b[1, 0] * torch.cos(psi) + self.c[1, 0]
        if torch.abs(cos_q2) > 1 + 1e-9:
            return False, q
        else:
            cos_q2 = torch.clamp(cos_q2, -1, 1)

        q[0] = torch.atan2(self.gc[0] * (self.a[0, 0] * torch.sin(psi) + self.b[0, 0] * torch.cos(psi) + self.c[0, 0]),
                           self.gc[0] * (self.a[0, 1] * torch.sin(psi) + self.b[0, 1] * torch.cos(psi) + self.c[0, 1]))
        q[1] = self.gc[1] * torch.acos(cos_q2)
        q[2] = torch.atan2(self.gc[2] * (self.a[2, 0] * torch.sin(psi) + self.b[2, 0] * torch.cos(psi) + self.c[2, 0]),
                           self.gc[2] * (self.a[2, 1] * torch.sin(psi) + self.b[2, 1] * torch.cos(psi) + self.c[2, 1]))

        q[3] = self.q4_v

        cos_q6 = self.a[5, 0] * torch.sin(psi) + self.b[5, 0] * torch.cos(psi) + self.c[5, 0]
        if torch.abs(cos_q6) > 1 + 1e-9:
            return False, q
        else:
            cos_q6 = torch.clamp(cos_q6, -1, 1)

        q[4] = torch.atan2(self.gc[4] * (self.a[4, 0] * torch.sin(psi) + self.b[4, 0] * torch.cos(psi) + self.c[4, 0]),
                           self.gc[4] * (self.a[4, 1] * torch.sin(psi) + self.b[4, 1] * torch.cos(psi) + self.c[4, 1]))
        q[5] = self.gc[5] * torch.acos(cos_q6)
        q[6] = torch.atan2(self.gc[6] * (self.a[6, 0] * torch.sin(psi) + self.b[6, 0] * torch.cos(psi) + self.c[6, 0]),
                           self.gc[6] * (self.a[6, 1] * torch.sin(psi) + self.b[6, 1] * torch.cos(psi) + self.c[6, 1]))
        return True, q

    def _get_psi_i(self, theta_i, joint_id):
        if joint_id in [0, 2, 4, 6]:
            a_p = self.gc[joint_id] * ((self.c[joint_id, 1] - self.b[joint_id, 1]) *
                                       torch.tan(theta_i) + (self.b[joint_id, 0] - self.c[joint_id, 0]))
            b_p = 2 * self.gc[joint_id] * (self.a[joint_id, 1] * torch.tan(theta_i) - self.a[joint_id, 0])
            c_p = self.gc[joint_id] * ((self.c[joint_id, 1] + self.b[joint_id, 1]) *
                                       torch.tan(theta_i) - (self.b[joint_id, 0] + self.c[joint_id, 0]))
            square_p = b_p ** 2 - 4 * a_p * c_p
            psi_ = []
            if square_p < 0:
                return psi_
            else:
                psi_1 = 2 * torch.atan((-b_p - torch.sqrt(square_p)) / (2 * a_p))
                theta_1 = torch.atan2(self.gc[joint_id] * (self.a[joint_id, 0] * torch.sin(psi_1) +
                                                           self.b[joint_id, 0] * torch.cos(psi_1) +
                                                           self.c[joint_id, 0]),
                                      self.gc[joint_id] * (self.a[joint_id, 1] * torch.sin(psi_1) +
                                                           self.b[joint_id, 1] * torch.cos(psi_1) +
                                                           self.c[joint_id, 1]))

                psi_2 = 2 * torch.atan((-b_p + torch.sqrt(square_p)) / (2 * a_p))
                theta_2 = torch.atan2(self.gc[joint_id] * (self.a[joint_id, 0] * torch.sin(psi_2) +
                                                           self.b[joint_id, 0] * torch.cos(psi_2) +
                                                           self.c[joint_id, 0]),
                                      self.gc[joint_id] * (self.a[joint_id, 1] * torch.sin(psi_2) +
                                                           self.b[joint_id, 1] * torch.cos(psi_2) +
                                                           self.c[joint_id, 1]))

                if torch.isclose(theta_i, theta_1) and torch.isclose(theta_i, theta_2):
                    psi_.append(psi_1)
                    psi_.append(psi_2)

                return psi_

        elif joint_id in [1, 5]:
            square_p = self.a[joint_id, 0] ** 2 + self.b[joint_id, 0] ** 2 - \
                       (self.c[joint_id, 0] - torch.cos(theta_i)) ** 2
            psi_ = []
            if square_p < 0:
                return psi_
            else:
                psi_1 = 2 * torch.atan((self.a[joint_id, 0] - torch.sqrt(square_p)) /
                                       (torch.cos(theta_i) + self.b[joint_id, 0] - self.c[joint_id, 0]))
                theta_1 = torch.acos(self.a[joint_id, 0] * torch.sin(psi_1) +
                                     self.b[joint_id, 0] * torch.cos(psi_1) + self.c[joint_id, 0])

                psi_2 = 2 * torch.atan((self.a[joint_id, 0] + torch.sqrt(square_p)) /
                                       (torch.cos(theta_i) + self.b[joint_id, 0] - self.c[joint_id, 0]))
                theta_2 = torch.acos(self.a[joint_id, 0] * torch.sin(psi_2) +
                                     self.b[joint_id, 0] * torch.cos(psi_2) +
                                     self.c[joint_id, 0])

                if torch.isclose(theta_i, theta_1) and torch.isclose(theta_i, theta_2):
                    psi_.append(psi_1)
                    psi_.append(psi_2)
                return psi_

        else:
            raise NotImplementedError()

    def _check_stationary(self, joint_id):
        # check whether contains stationary point
        if joint_id in [0, 2, 4, 6]:
            self.a_t = self.c[joint_id, 0] * self.b[joint_id, 1] - self.b[joint_id, 0] * self.c[joint_id, 1]
            self.b_t = self.a[joint_id, 0] * self.c[joint_id, 1] - self.c[joint_id, 0] * self.a[joint_id, 1]
            self.c_t = self.a[joint_id, 0] * self.b[joint_id, 1] - self.b[joint_id, 0] * self.a[joint_id, 1]

            square_t = self.a_t ** 2 + self.b_t ** 2 - self.c_t ** 2
            if torch.abs(square_t) > 0 + 1e-3:
                return 0
            else:
                return 1

    def _get_interval_i(self, joint_id):
        u_bound = []
        l_bound = []
        psi_singularity = []
        theta_l = self.joint_limits[joint_id, 0]
        theta_u = self.joint_limits[joint_id, 1]
        # Check for pivot joints
        if joint_id in [0, 2, 4, 6]:
            # check whether stationary point
            has_stationary = self._check_stationary(joint_id)
            if has_stationary:
                if torch.abs(self.b_t - self.c_t) < 1e-6:
                    psi_singularity.append(np.pi - self.singularity_eps)
                    psi_singularity.append(-np.pi + self.singularity_eps)
                    u_bound.append(np.pi - self.singularity_eps)
                    l_bound.append(-np.pi + self.singularity_eps)
                else:
                    psi_singularity.append(2 * torch.atan(self.a_t / (self.b_t - self.c_t)) - self.singularity_eps)
                    psi_singularity.append(2 * torch.atan(self.a_t / (self.b_t - self.c_t)) + self.singularity_eps)
                    u_bound.append(2 * torch.atan(self.a_t / (self.b_t - self.c_t)) - self.singularity_eps)
                    l_bound.append(2 * torch.atan(self.a_t / (self.b_t - self.c_t)) + self.singularity_eps)
            for theta_i in [theta_l, theta_u]:
                psi_i = self._get_psi_i(theta_i, joint_id)
                for psi_tmp in psi_i:
                    dtheta = self.a_t * torch.sin(psi_tmp) + self.b_t * torch.cos(psi_tmp) + self.c_t
                    if theta_i * dtheta > 0:
                        u_bound.append(psi_tmp)
                    elif theta_i * dtheta < 0:
                        l_bound.append(psi_tmp)
                    else:
                        pass

        elif joint_id in [1, 5]:
            # psi_test = torch.linspace(-np.pi, np.pi, 1000)
            # theta_test = self.gc[joint_id] * torch.acos(self.a[joint_id, 0] * torch.sin(psi_test) + self.b[joint_id, 0] * torch.cos(psi_test) + self.c[joint_id, 0])

            # plt.figure(1)
            # plt.clf()
            # plt.plot(psi_test, theta_test, c='b')
            # plt.plot(psi_test, torch.ones_like(psi_test) * theta_l, c='r')
            # plt.plot(psi_test, torch.ones_like(psi_test) * theta_u, c='r')
            # plt.plot(psi_test, torch.zeros_like(psi_test), linestyle = '--',  c='r', )
            # plt.xlim(-np.pi, np.pi)
            # plt.ylim(-np.pi, np.pi)
            # plt.draw()
            # plt.pause(0.1)

            # symmetry at theta=0, only check upper part
            theta_i = theta_u
            psi_i = self._get_psi_i(theta_i, joint_id)

            for psi_tmp in psi_i:
                dtheta = -self.a[joint_id, 0] * torch.cos(psi_tmp) + self.b[joint_id, 0] * torch.sin(psi_tmp)
                if dtheta > 0:
                    u_bound.append(psi_tmp)
                elif dtheta < 0:
                    l_bound.append(psi_tmp)
                else:
                    pass

        l_bound.append(torch.tensor(-np.pi))
        u_bound.append(torch.tensor(np.pi))

        if not (len(l_bound) == len(u_bound)):
            raise ValueError("Number of upper and lower bound are not the same")
        else:
            num_bound = len(l_bound)

        l_bound_torch = torch.cat((torch.tensor(l_bound).unsqueeze(1), torch.zeros((num_bound, 1))), dim=1)
        u_bound_torch = torch.cat((torch.tensor(u_bound).unsqueeze(1), torch.ones((num_bound, 1))), dim=1)

        bound = torch.cat((l_bound_torch, u_bound_torch), dim=0)
        bound = bound[bound[:, 0].sort(dim=0).indices]

        index = 0
        interval_ = []
        while index < bound.shape[0]:
            if bound[index, 1] == 0:
                iter_low = bound[index, 0]
                index += 1
                if bound[index, 1] == 1:
                    iter_high = bound[index, 0]
                    # interval.extend([[iter_low, 0.], [iter_high, 1.]])
                    interval_.append([iter_low, iter_high])
                    index += 1
                else:
                    continue
            else:
                index += 1
        interval_ = torch.tensor(interval_)
        return interval_, psi_singularity

    def _get_intersection_interval(self, interval_1, interval_2):
        intersections = []
        for int_1 in interval_1:
            for int_2 in interval_2:
                low = torch.tensor((int_1[0], int_2[0])).max()
                high = torch.tensor((int_1[1], int_2[1])).min()
                if high >= low:
                    intersections.append([low, high])
        return torch.tensor(intersections)

    def get_feasible_interval(self, pose, gc):
        p, rot = self._transform_tcp2ee(pose)
        if not self._get_auxiliary_parameter(p, rot, gc):
            return None
        interval_prev = torch.tensor([[-np.pi, np.pi]])
        for i in range(7):
            interval_i, singularity = self._get_interval_i(i)
            interval_prev = self._get_intersection_interval(interval_prev, interval_i)

        return interval_prev

if __name__ == "__main__":
    tol = 1e-4
    import quaternion

    for i in range(100):
        quat_test = np.random.rand(4)
        quat_test = quat_test / np.linalg.norm(quat_test)
        rot_true = quaternion.from_float_array(quat_test)
        rot_torch = Rotation.from_quat(torch.tensor(quat_test))
        if np.linalg.norm(rot_torch.as_matrix().detach().numpy() - quaternion.as_rotation_matrix(rot_true)) > 1e-8:
            print("Rotation to Matrix error")

        if np.linalg.norm(rot_torch.from_matrix(rot_torch.as_matrix()).as_quat().detach().numpy() - quat_test) > 1e-8:
            print("Matrix to Quaternion Error")

        if np.linalg.norm(rot_torch.as_rotvec().detach().numpy() - quaternion.as_rotation_vector(rot_true)) > 1e-8:
            print("Quaternion to Rotation Vector Error")

        if np.linalg.norm(rot_torch.from_rotvec(rot_torch.as_rotvec()).as_quat().detach().numpy() - quat_test) > 1e-8:
            print("Rotation Vector to Quaternion Error")

    for i in range(10):
        yzy_test = torch.rand(3) * torch.tensor([2 * np.pi, np.pi, 2 * np.pi]) - torch.tensor([np.pi, 0, np.pi])
        kine_test = KinematicsTorch(tcp_pos=torch.rand(3),
                                    tcp_quat=Rotation.from_euler_yzy(yzy_test).as_quat())
        q_test = torch.deg2rad(torch.rand((7,)) * 360 - 180).requires_grad_()
        pose_test = kine_test.forward_kinematics(q_test)

        psi_test = torch.rand(1)
        gc_test = torch.tensor([1, 1, 1])
        res, q_inv = kine_test.inverse_kinematics(pose_test, psi_test, gc_test)
        if res:
            pose_inv = kine_test.forward_kinematics(q_inv)
            position_error = torch.norm(pose_test[:3] - pose_inv[:3])
            if position_error > tol:
                print("Position Error:", position_error.item())
            rotation_error = torch.norm(pose_test[3:] - pose_inv[3:])
            if rotation_error > tol:
                print("Rotation Error", rotation_error.item())

        # interval_test = kine_test.get_feasible_interval(pose_test.requires_grad_(True), gc_test)
        # print("Interval:", interval_test)

        jac = kine_test.get_jacobian(q_test)
        dq_test = torch.ones(7).double() * 1e-3
        pose_dq_test = kine_test.forward_kinematics(dq_test + q_test) - pose_test
        pose_jac_test = jac @ dq_test
        pose_jac_error = torch.norm(pose_dq_test - pose_jac_test)
        if pose_jac_error > tol:
            print("Jacobian Error:", pose_jac_error)

    # Test for calculating redundancy
    for i in range(10):
        # yzy_test = torch.rand(3) * torch.tensor([2 * np.pi, np.pi, 2 * np.pi]) - torch.tensor([np.pi, 0, np.pi])
        # kine_test = KinematicsTorch(tcp_pos=torch.rand(3),
        #                             tcp_quat=Rotation.from_euler_yzy(yzy_test).as_quat())
        # q_test = torch.rand(7) * (kine_test.joint_limits[:, 1]) * 2 + kine_test.joint_limits[:, 0]

        yzy_test = torch.rand(3) * torch.tensor([2 * np.pi, np.pi, 2 * np.pi]) - torch.tensor([np.pi, 0, np.pi])
        kine_test = KinematicsTorch(tcp_pos=torch.tensor([0.1, 0.2, 0.5]),
                                tcp_quat=torch.tensor([1., 0.5, 0.1, 0.6])/torch.norm(torch.tensor([1., 0.5, 0.1, 0.6])))
        q_test = torch.tensor([-0.1, 0.2, -0.3, -0.4, -0.5, 0.6, -0.7])

        pose_test = kine_test.forward_kinematics(q_test)
        gc = torch.sign(q_test[[1, 3, 5]])
        psi_test = kine_test.get_redundancy(q_test)
        res, q_test_out = kine_test.inverse_kinematics(pose_test, psi_test, gc)

        q_error = torch.norm(q_test - q_test_out)
        if q_error > tol:
            print("psi is not correct")

    print("Finish!")

reason: import numpy