Python numpy.nan() Examples

def trix(df, n):
    """Calculate TRIX for given data.
    
    :param df: pandas.DataFrame
    :param n: 
    :return: pandas.DataFrame
    """
    EX1 = df['Close'].ewm(span=n, min_periods=n).mean()
    EX2 = EX1.ewm(span=n, min_periods=n).mean()
    EX3 = EX2.ewm(span=n, min_periods=n).mean()
    i = 0
    ROC_l = [np.nan]
    while i + 1 <= df.index[-1]:
        ROC = (EX3[i + 1] - EX3[i]) / EX3[i]
        ROC_l.append(ROC)
        i = i + 1
    Trix = pd.Series(ROC_l, name='Trix_' + str(n))
    df = df.join(Trix)
    return df 

def __init__(self, gamma, sinDec_range=(-1., 1.), sinDec_bandwidth=0.1,
                 e_range=(0., np.inf), E0=1., GeV=1., seed=None):
        super(PointSourceInjector, self).__init__()

        self.gamma = gamma

        self.sinDec_range = sinDec_range
        self.sinDec_bandwidth = sinDec_bandwidth
        self.src_dec = np.nan

        self.e_range = e_range

        self.E0 = E0
        self.GeV = GeV

        self.random = np.random.RandomState(seed) 

def get(self):
        """Get the current evaluation result.

        Returns
        -------
        name : str
           Name of the metric.
        value : float
           Value of the evaluation.
        """
        self._update()  # update metric at this time
        if self.num is None:
            if self.num_inst == 0:
                return (self.name, float('nan'))
            else:
                return (self.name, self.sum_metric / self.num_inst)
        else:
            names = ['%s'%(self.name[i]) for i in range(self.num)]
            values = [x / y if y != 0 else float('nan') \
                for x, y in zip(self.sum_metric, self.num_inst)]
            return (names, values)

    # pylint: disable=arguments-differ, too-many-nested-blocks 

def _recall_prec(self):
        """ get recall and precision from internal records """
        n_fg_class = max(self._n_pos.keys()) + 1
        prec = [None] * n_fg_class
        rec = [None] * n_fg_class

        for l in self._n_pos.keys():
            score_l = np.array(self._score[l])
            match_l = np.array(self._match[l], dtype=np.int32)

            order = score_l.argsort()[::-1]
            match_l = match_l[order]

            tp = np.cumsum(match_l == 1)
            fp = np.cumsum(match_l == 0)

            # If an element of fp + tp is 0,
            # the corresponding element of prec[l] is nan.
            with np.errstate(divide='ignore', invalid='ignore'):
                prec[l] = tp / (fp + tp)
            # If n_pos[l] is 0, rec[l] is None.
            if self._n_pos[l] > 0:
                rec[l] = tp / self._n_pos[l]

        return rec, prec 

def _average_precision(self, rec, prec):
        """
        calculate average precision, override the default one,
        special 11-point metric

        Params:
        ----------
        rec : numpy.array
            cumulated recall
        prec : numpy.array
            cumulated precision
        Returns:
        ----------
        ap as float
        """
        if rec is None or prec is None:
            return np.nan
        ap = 0.
        for t in np.arange(0., 1.1, 0.1):
            if np.sum(rec >= t) == 0:
                p = 0
            else:
                p = np.max(np.nan_to_num(prec)[rec >= t])
            ap += p / 11.
        return ap 

def compute_cor_loc(num_gt_imgs_per_class,
                    num_images_correctly_detected_per_class):
  """Compute CorLoc according to the definition in the following paper.

  https://www.robots.ox.ac.uk/~vgg/rg/papers/deselaers-eccv10.pdf

  Returns nans if there are no ground truth images for a class.

  Args:
    num_gt_imgs_per_class: 1D array, representing number of images containing
        at least one object instance of a particular class
    num_images_correctly_detected_per_class: 1D array, representing number of
        images that are correctly detected at least one object instance of a
        particular class

  Returns:
    corloc_per_class: A float numpy array represents the corloc score of each
      class
  """
  return np.where(
      num_gt_imgs_per_class == 0,
      np.nan,
      num_images_correctly_detected_per_class / num_gt_imgs_per_class) 

def __init__(self,
               num_groundtruth_classes,
               matching_iou_threshold=0.5,
               nms_iou_threshold=1.0,
               nms_max_output_boxes=10000):
    self.per_image_eval = per_image_evaluation.PerImageEvaluation(
        num_groundtruth_classes, matching_iou_threshold, nms_iou_threshold,
        nms_max_output_boxes)
    self.num_class = num_groundtruth_classes

    self.groundtruth_boxes = {}
    self.groundtruth_class_labels = {}
    self.groundtruth_is_difficult_list = {}
    self.num_gt_instances_per_class = np.zeros(self.num_class, dtype=int)
    self.num_gt_imgs_per_class = np.zeros(self.num_class, dtype=int)

    self.detection_keys = set()
    self.scores_per_class = [[] for _ in range(self.num_class)]
    self.tp_fp_labels_per_class = [[] for _ in range(self.num_class)]
    self.num_images_correctly_detected_per_class = np.zeros(self.num_class)
    self.average_precision_per_class = np.empty(self.num_class, dtype=float)
    self.average_precision_per_class.fill(np.nan)
    self.precisions_per_class = []
    self.recalls_per_class = []
    self.corloc_per_class = np.ones(self.num_class, dtype=float) 

def testReturnsCorrectNanLoss(self):
    batch_size = 3
    num_anchors = 10
    code_size = 4
    prediction_tensor = tf.ones([batch_size, num_anchors, code_size])
    target_tensor = tf.concat([
        tf.zeros([batch_size, num_anchors, code_size / 2]),
        tf.ones([batch_size, num_anchors, code_size / 2]) * np.nan
    ], axis=2)
    weights = tf.ones([batch_size, num_anchors])
    loss_op = losses.WeightedL2LocalizationLoss()
    loss = loss_op(prediction_tensor, target_tensor, weights=weights,
                   ignore_nan_targets=True)

    expected_loss = (3 * 5 * 4) / 2.0
    with self.test_session() as sess:
      loss_output = sess.run(loss)
      self.assertAllClose(loss_output, expected_loss) 

def lnlike(self, pars):
        # Pull theta out of pars
        theta = pars[:self.Nbins]

        # Generate the inner summation
        gamma = np.ones_like(self.bin_idx) * np.nan
        good = (self.bin_idx < self.Nbins) & (self.bin_idx >= 0)  # nans in q get put in nonexistent bins
        gamma[good] = self.Nobs * self.censoring_fcn(self.mcmc_samples[good]) * theta[self.bin_idx[good]]
        summation = np.nanmean(gamma, axis=1)

        # Calculate the integral
        I = self._integral_fcn(theta)

        # Generate the log-likelihood
        ll = -I + np.nansum(np.log(summation))
        return ll 

def ada2MonPix(self, xyCoor = tuple):
        
        # check argument values
        if xyCoor is None:
            raise ValueError("No coordinate values have been specified.")
        elif not isinstance(xyCoor, tuple):
            raise TypeError("XY coordinates must be given as tuple.")
        elif isinstance(xyCoor, tuple) and len(xyCoor) is not 2: 
            raise ValueError("Wrong number of coordinate dimensions")

        if np.isnan(xyCoor[0]) and np.isnan(xyCoor[1]):
            monPix = (np.nan, np.nan)
            return monPix

        # convert so point of gaze on monitor is accurate
        monPix = (int(xyCoor[0] * self.win.getSizePix()[0]),
                  int(xyCoor[1] * self.win.getSizePix()[1]))                                
        return monPix        
    
        
# ----- Functions for collecting eye and gaze data -----
      
    # function for collecting gaze coordinates in psychopy pixel coordinate 
    # system. currently written to return the average (x, y) position of both 
    # eyes, but can be easily rewritten to return data from one or both eyes 

def DistMAXBAS(FPL,SPMAXBAS, q_uz):
    
    MAXBAS = np.nanmax(SPMAXBAS)
    FPLArray = FPL.ReadAsArray()
    rows = FPL.RasterYSize
    cols = FPL.RasterXSize
    NoDataValue = np.float32(FPL.GetRasterBand(1).GetNoDataValue())
    FPLArray[FPLArray == NoDataValue] = np.nan # replace novalue cells by nan
    
    MaxFPL = np.nanmax(FPLArray)
    MinFPL = np.nanmin(FPLArray)
#    resize_fun = lambda x: np.round(((((x - min_dist)/(max_dist - min_dist))*(1*maxbas - 1)) + 1), 0)
    resize_fun = lambda x: ((((x - MinFPL)/(MaxFPL - MinFPL))*(1*MAXBAS - 1)) + 1)
    
    NormalizedFPL = resize_fun(FPLArray)
    
    for x in range(rows):
        for y in range(cols):
            if not np.isnan(FPLArray[x,y]):# FPLArray[x,y] != np.nan: #NoDataValue:
                q_uz[x,y,:] = Routing.TriangularRouting(q_uz[x,y,:], NormalizedFPL[x,y]) 
    
    return q_uz 

def test_linear_sum_assignment_input_validation():
    assert_raises(ValueError, linear_sum_assignment, [1, 2, 3])

    C = [[1, 2, 3], [4, 5, 6]]
    assert_array_equal(linear_sum_assignment(C), linear_sum_assignment(np.asarray(C)))
    # assert_array_equal(linear_sum_assignment(C),
    #                    linear_sum_assignment(matrix(C)))

    I = np.identity(3)
    assert_array_equal(linear_sum_assignment(I.astype(np.bool)), linear_sum_assignment(I))
    assert_raises(ValueError, linear_sum_assignment, I.astype(str))

    I[0][0] = np.nan
    assert_raises(ValueError, linear_sum_assignment, I)

    I = np.identity(3)
    I[1][1] = np.inf
    assert_raises(ValueError, linear_sum_assignment, I) 

def __init__(self, graph: MultiplexNet):

        if not isinstance(graph, MultiplexNet):
            raise ValueError("Unrecognized graph object")

        self.graph = graph

        # create average adjacency matrix
        self.avg_adj = np.zeros((self.graph.nodes, self.graph.nodes))
        connections = self.graph.connections
        connections[connections == 0] = np.nan

        for a in self.graph.adj_matrices:
            self.avg_adj = np.nansum(np.dstack((a, self.avg_adj)), 2)
        self.avg_adj = self.avg_adj / self.graph.connections

        # impute NaN values in average adjacency with SVD-EM algorithm
        if np.sum(np.isnan(self.avg_adj)) > 0:
            self.avg_adj = svdEM(self.avg_adj)

        # layer impact balancing parameters
        self.lambdas = [(1. / samples.shape[0]) ** 2 for samples in self.graph.samples] 

def _distance_testing(func, val):
    assert func in [similarity.euclidean_dist, similarity.chi_squared_dist, similarity.agreement_dist]

    assert np.allclose(func(np.array([1, 2]), np.array([1, 2])), 0)
    assert np.allclose(func(np.arange(16), np.arange(16)), 0)

    v = np.array([1 if i != 0 else 2 for i in range(100)])
    assert np.allclose(func(v, np.ones(100)), val / 100)

    assert np.allclose(func(np.array([1, 1, 1]), np.array([1, 1, 2])), val / 3)
    assert np.allclose(func(np.array([1, np.nan, 1]), np.array([1, np.nan, 2])), val / 2)
    assert np.allclose(func(np.array([np.nan, 1, 1]), np.array([1, np.nan, 2])), val)

    assert np.allclose(func(np.array([np.nan, 1, 1]), np.array([1, np.nan, 2]), missing=0.5), np.nan, equal_nan=True)
    assert np.allclose(func(np.array([np.nan, 1]), np.array([1, np.nan])), np.nan, equal_nan=True)

    with pytest.raises(AssertionError):
        func(np.array([1]), np.array([1, 1]))

    with pytest.raises(AttributeError):
        func([2], [3]) 

def test_init():
    a0 = np.random.random((10, 10))
    a0 = (a0 * a0.T) / 2
    a1 = np.random.random((10, 10))
    a1 = (a1 * a1.T) / 2
    sample_names = ['sample_{}'.format(i) for i in range(10)]

    with pytest.raises(ValueError):
        SumoNMF(a0)

    graph = MultiplexNet(adj_matrices=[a0, a1], node_labels=np.array(sample_names))
    nmf = SumoNMF(graph)
    assert np.array_equal((a0 + a1) / 2, nmf.avg_adj)

    # missing values in one layer
    a0[0, 1], a0[1, 0] = np.nan, np.nan
    graph = MultiplexNet(adj_matrices=[a0, a1], node_labels=np.array(sample_names))
    nmf = SumoNMF(graph)
    assert np.allclose(nmf.avg_adj[0, 1], a1[0, 1])

    # missing values in both layers
    a1[0, 1], a1[1, 0] = np.nan, np.nan
    graph = MultiplexNet(adj_matrices=[a0, a1], node_labels=np.array(sample_names))
    nmf = SumoNMF(graph)
    assert not np.isnan(nmf.avg_adj[0, 1]) 

def arccos(X):
	''' Compute the arccos of an AutoDiff object and its derivative.

		INPUTS
		======
		X: an AutoDiff object or constant

		RETURNS
		=======
		A new AutoDiff object or scalar with calculated value and derivative.

		EXAMPLES
		========
		>>> X = AutoDiff(0.5, 2)
		>>> arccosAutoDiff = arccos(X)
		>>> arccosAutoDiff.val
		1.0471975511965976
		>>> arccosAutoDiff.der
		-2.3094010767585034
		>>> arccosAutoDiff.jacobian
		-1.1547005383792517
		'''
	try:
		# Is another ADT
		new_val = np.arccos(X.val) #if (-1 <= X.val and X.val <= 1) else np.nan
		new_der = (-1/np.sqrt(1-X.val**2))*X.der #if (-1 < X.val and X.val < 1) else np.nan
		new_jacobian = (-1/np.sqrt(1-X.val**2))*X.jacobian #if (-1 < X.val and X.val < 1) else np.nan

		return AutoDiff(new_val, new_der, X.n, 0, new_jacobian)
	except AttributeError:
		try:
			return Dual(np.arccos(X.Real), -X.Dual/np.sqrt(1-X.Real**2))		
		except AttributeError:
			try:
				return Dual(arccos(X.Real), -X.Dual/sqrt(1-X.Real**2))
			except AttributeError:
			# Constant
				return_val = np.arccos(X)
				return return_val

# arc tangent 

def values_for_reg_arr():
    return np.array([[-2., 1.],
                     [np.nan, 5.],
                     [3., 3.],
                     [4., 4.2]]) 

def test_mask_var(data_for_reg_calcs, region):
    # Test region masks first row and second column of test data.  Note that
    # first element of second row is np.nan in initial dataset.
    expected_data = [[np.nan, np.nan],
                     [np.nan, np.nan],
                     [3., np.nan],
                     [4., np.nan]]
    expected = data_for_reg_calcs.copy(deep=True)
    expected.values = expected_data
    result = region.mask_var(data_for_reg_calcs)
    xr.testing.assert_identical(result, expected) 

def test_mask_var_non_aospy_names(data_reg_alt_names, region):
    # Test region masks first row and second column of test data.  Note that
    # first element of second row is np.nan in initial dataset.
    expected_data = [[np.nan, np.nan],
                     [np.nan, np.nan],
                     [3., np.nan],
                     [4., np.nan]]
    expected = data_reg_alt_names.copy(deep=True)
    expected.values = expected_data
    result = region.mask_var(data_reg_alt_names, lon_str=_alt_names[LON_STR],
                             lat_str=_alt_names[LAT_STR])
    xr.testing.assert_identical(result, expected) 

def decrease_wavelengths(original_dataset, A_set = 2, specific_subset = None):
    '''
    Takes in the original, full dataset and removes specific wavelengths, or only keeps every
    multipl of A_set. Returns a new, smaller dataset that should be easier to solve
    
    Args:
        original_dataset (DataFrame):   the data to be processed
        A_set (float, optional):  optional user-provided multiple of wavelengths to keep. i.e. if
                                    3, every third value is kept. Default is 2.
        specific_subset (list or dict, optional): If the user already knows which wavelengths they would like to
                                    remove, then a list containing these can be included.
        
    Returns:
        DataFrame with the smaller dataset
    
    '''
    if specific_subset != None:
        if not isinstance(specific_subset, (list, dict)):
            raise RuntimeError("subset must be of type list or dict!")
             
        if isinstance(specific_subset, dict):
            lists1 = sorted(specific_subset.items())
            x1, y1 = zip(*lists1)
            specific_subset = list(x1)
            
        new_D = pd.DataFrame(np.nan,index=original_dataset.index, columns = specific_subset)
        for t in original_dataset.index:
            for l in original_dataset.columns.values:
                if l in subset:
                    new_D.at[t,l] = self.model.D[t,l]           
    else:
        count=0
        for l in original_dataset.columns.values:
            remcount = count%A_set
            if remcount==0:
                original_dataset.drop(columns=[l],axis = 1)
            count+=1
        new_D = original_dataset[original_dataset.columns[::A_set]]     
    return new_D 

def calc_mean_diff_conf(zipper, feat_subset):
    """
    Calculates the confidence intervals for numerical features.

    :param zipper: Zipper created from a DataCollection.
    :param feat_subset: An iterable of features for which the confidence intervals shall be calculated.
    :return:
    """
    # initialize dictionary which stores the conf_invs
    conf_invs = dict()

    for feat in feat_subset:  # run through all variables

        # initiate dict in dict for d1 vs d2, d2 vs d3 etc. per feature
        conf_invs[feat] = dict()

        for i in range(len(zipper[feat]) - 1):  # select dataset1
            for j in range(i + 1, len(zipper[feat])):  # select dataset2

                # only calculate score if there are values in each dataset
                if zipper[feat][i] and zipper[feat][j]:
                    interval = np.round(calc_mean_diff(zipper[feat][i], zipper[feat][j]), 2)

                    # indicator = True if 0 is not in the interval
                    if interval[0] >= 0 or interval[1] <= 0:
                        flag = True
                    else:
                        flag = False

                    conf_invs[feat][i + 1, j + 1] = (interval, flag)

                # if one or both sets are empty
                else:
                    conf_invs[feat][i + 1, j + 1] = (np.nan, np.nan)

    return conf_invs 

def p_correction(p_values):
    """
    Corrects p_values for multiple testing.

    :param p_values: Dictionary storing p_values with corresponding feature names as keys.
    :return: DataFrame which shows the results of the analysis; p-value, corrected p-value and boolean indicating \
    significance.
    """

    p_trans = _transform_p_dict(p_values)

    # get and drop features which are NaN to skip them in multitest correction
    nan_features = p_trans[pd.isnull(p_trans[0])]
    p_trans = p_trans.dropna(axis=0, subset=[0])

    # extract p_value column to pass into multiple testing correction
    p_val_col = p_trans[0].sort_values()

    # add NaN features back to p_trans to include them into result table later on
    p_trans = pd.concat([p_trans, nan_features])

    # raise Error if no p_values where calculated that can be passed into multiple test correction
    if p_val_col.values.size == 0:
        # unpack the p_values which are stored in 2 layer nested dicts.
        nested_values = []
        for value in p_values.values():
            nested_values.append(*value.values())

        # if all p_values are nan, return an all nan result table
        if pd.isnull(nested_values).all():
            result_table = _create_result_table(None, p_val_col, p_trans, conf_invs, counts)
            return result_table.sort_index()

        raise ValueError("No p_values have been submitted into multiple test correction.")

    # correct p-values
    result = multipletests(p_val_col.values)

    return result, p_val_col, p_trans 

def __init__(self, df_list, df_names, categorical_feats, numerical_feats=None, exclude_nan_feats=True):
        """
        Initialize DataCollection object.

        :param df_list: List containing pandas.DataFrames
        :param df_names: List of strings containing the names of the datasets in df_list
        :param categorical_feats: List of categorical features
        """
        super().__init__(df_list)

        # check if there is an empty dataframe in the datacollection
        if any(df.empty for df in self):
            warnings.warn("One of the dataframes in the DataCollection is empty!", UserWarning)

        self.df_names = df_names
        self.categorical_feats = categorical_feats

        # exclude all nan features
        if exclude_nan_feats:
            self.exclude_empty_feats()

            # update categorical feature name list to only include the ones present after feature elimination
            self.update_categorical_feats()

        # get numerical features if not given
        if numerical_feats is None:
            self.numerical_feats = self.get_numerical_features(self.categorical_feats)
        else:
            self.numerical_feats = numerical_feats 

def test_transformer(rbp):
    """Test the transformer function
    """
    posTransf = PositionTransformer(ALL_LANDMARKS, "{0}/dataloader_files/position_transformer.pkl".format(rbp))
    x = np.arange(len(ALL_LANDMARKS)).reshape((1, -1)).astype(float)
    x[0, 1] = np.nan
    res = posTransf.transform(x)
    assert set(res.keys()) == set(["dist_" + x for x in ALL_LANDMARKS])
    print("Test passed!") 

def __getitem__(self, idx):
        if self.fasta is None:
            from fasta_utils import FastaFile
            self.fasta = FastaFile(self.fasta_file)

        out = {}
        if self.MISO_AS:
            gene = self.genes[idx]
            inputs, ranges = self.get_seq(gene)
            out['inputs'] = inputs
            if self.Y is not None:
                out['targets'] = self.Y.get_target(gene.geneName)
            else:
                out['targets'] = np.nan
            out['metadata'] = {}
            out['metadata']['geneName'] = gene.geneName
            out['metadata']['chrom'] = gene.chrom
            out['metadata']['strand'] = gene.strand
            out['metadata']['start'] = gene.start
            out['metadata']['stop'] = gene.stop
            out['metadata']['extracted_regions'] = ranges

        else:
            spliceSite = self.spliceSites[idx]
            out['inputs'] = spliceSite.get_seq(self.fasta)
            out['metadata'] = {}
            out['metadata']['geneID'] = spliceSite.geneID
            out['metadata']['transcriptID'] = spliceSite.transcriptID
            out['metadata']['biotype'] = spliceSite.biotype
            out['metadata']['order'] = spliceSite.order
            out['metadata']['ranges'] = GenomicRanges(spliceSite.chrom,
                                                      spliceSite.grange[0] - 1,  # use 0-base indexing
                                                      spliceSite.grange[1],
                                                      spliceSite.geneID,
                                                      spliceSite.strand)

        return out 

def nans(shape, dtype=float):
    a = np.empty(shape, dtype)
    a.fill(np.nan)
    return a 

def __init__(self, eta, df, loc=0., scale=1.):
        self.eta = eta
        self.params = (df, loc, scale)
        self._chi2 = scipy.stats.chi2(df, loc, scale)
        self.eta_err = np.nan
        self.ks = np.nan 

def __init__(self, df=np.nan, **others):
        self.df = df
        self.others = others 

def __init__(self, eta, df, loc=0., scale=1.):
        params = np.empty(shape=(2, 3))
        params[:, 0] = df
        params[:, 1] = loc
        params[:, 2] = scale

        self.eta = eta
        self.params = params

        self._chi2 = tuple(scipy.stats.chi2(*p) for p in params)

        self.eta_err = np.nan
        self.ks = (np.nan, np.nan) 

def __init__(self, df=np.nan, left={}, right={}):
        self.df = np.empty(2)
        self.df[:] = df
        self.others = (left, right) 

def _calculateCoreStats(self):
        """
        @summary: This function calculates the standard PAM statistics
        """
        # Number of species at each site
        self.alpha = np.sum(self.pamData, axis=1)
        
        # Number of sites for each species
        self.omega = np.sum(self.pamData, axis=0)
        
        # Calculate the number of species by looking for columns that have any 
        #     presences.  This will let the stats ignore empty columns
        self.numSpecies = np.sum(np.any(self.pamData, axis=0))

        # Calculate the number of sites that have at least one species present
        self.numSites = np.sum(np.any(self.pamData, axis=1))
        
        # Site statistics
        self.alphaProp = self.alpha.astype(float) / self.numSpecies
        self.phi = self.pamData.dot(self.omega)
        # phiAvgProp can have np.nan values if empty rows and columns, 
        #     set to zero
        self.phiAvgProp = np.nan_to_num(
            self.phi.astype(float) / (self.numSites * self.alpha))
        
        # Species statistics
        self.omegaProp = self.omega.astype(float) / self.numSites
        self.psi = self.alpha.dot(self.pamData)
        # psiAvgProp can produce np.nan for empty row and columns, set to zero
        self.psiAvgProp = np.nan_to_num(
            self.psi.astype(float) / (self.numSpecies * self.omega))
    
    # ........................... 

def uniform_sampling_within_a_bin(self, bin_idx: int):
        num_bins = len(self.distribution_bins)
        if bin_idx == num_bins:
            return np.nan
        elif self.is_categorical:
            return self.distribution_bins[bin_idx]
        elif bin_idx < num_bins - 1:
            return uniform(self.distribution_bins[bin_idx], self.distribution_bins[bin_idx + 1])
        else:
            # sample from the last interval where the right edge is missing in self.distribution_bins
            neg_2, neg_1 = self.distribution_bins[-2:]
            return uniform(neg_1, 2 * neg_1 - neg_2) 

def classifer_metrics(label, pred):
    """
    computes f1, precision and recall on the entity class
    """
    prediction = np.argmax(pred, axis=1)
    label = label.astype(int)

    pred_is_entity = prediction != not_entity_index
    label_is_entity = label != not_entity_index

    corr_pred = (prediction == label) == (pred_is_entity == True)

    #how many entities are there?
    num_entities = np.sum(label_is_entity)
    entity_preds = np.sum(pred_is_entity)

    #how many times did we correctly predict an entity?
    correct_entitites = np.sum(corr_pred[pred_is_entity])

    #precision: when we predict entity, how often are we right?
    precision = correct_entitites/entity_preds
    if entity_preds == 0:
        precision = np.nan

    #recall: of the things that were an entity, how many did we catch?
    recall = correct_entitites / num_entities
    if num_entities == 0:
        recall = np.nan
    f1 = 2 * precision * recall / (precision + recall)
    return precision, recall, f1 

def test_stn():
    np.set_printoptions(threshold=np.nan)
    num_filter = 2  # conv of loc net
    kernel = (3, 3)  # conv of loc net
    num_hidden = 6  # fc of loc net
    for n in [1, 2, 3, 4]:
        for c in [1, 2, 3, 4]:
            for h in [5, 9, 13, 17]:  # for convenience test, this third and forth input dim should be 4x + 1
                for w in [5, 9, 13, 17]:
                    data_shape = (n, c, h, w)
                    target_shape = (int((data_shape[2]+1)/2), int((data_shape[3]+1)/2))
                    data = mx.sym.Variable(name="data")
                    loc = mx.sym.Convolution(data=data, kernel=kernel, pad=(1, 1), num_filter=num_filter, name="loc_conv")
                    loc = mx.sym.Flatten(data=loc)
                    loc = mx.sym.FullyConnected(data=loc, num_hidden=num_hidden, name="loc_fc")
                    stn = mx.sym.SpatialTransformer(data=data, loc=loc, target_shape=target_shape,
                                                    transform_type="affine", sampler_type="bilinear")
                    arg_names = stn.list_arguments()
                    arg_shapes, out_shapes, _ = stn.infer_shape(data=data_shape)
                    # check shape
                    assert out_shapes[0] == (data_shape[0], data_shape[1], target_shape[0], target_shape[1])
                    dev = default_context()
                    #dev = mx.gpu(0)
                    args = {}
                    args['data'] = mx.random.normal(0, 1, data_shape, ctx=mx.cpu()).copyto(dev)
                    args['loc_conv_weight'] = mx.nd.zeros((num_filter, data_shape[1], kernel[0], kernel[1]), ctx=dev)
                    args['loc_conv_bias'] = mx.nd.zeros((num_filter,), ctx=dev)
                    args['loc_fc_weight'] = mx.nd.zeros((6, num_filter*data_shape[2]*data_shape[3]), ctx=dev)
                    args['loc_fc_bias'] = mx.nd.array([0.5, 0, 0, 0, 0.5, 0], ctx=dev)
                    grad_grad = [mx.nd.zeros(shape, ctx=dev) for shape in arg_shapes]
                    exe = stn.bind(dev, args=args, args_grad=grad_grad)
                    exe.forward(is_train=True)
                    out = exe.outputs[0].asnumpy()
                    # check forward
                    assert_almost_equal(out, args['data'].asnumpy()[:, :, h//4:h-h//4, w//4:w-w//4], rtol=1e-2, atol=1e-4)
                    out_grad = mx.nd.ones(out.shape, ctx=dev)
                    exe.backward([out_grad])
                    # check backward
                    assert_almost_equal(out_grad.asnumpy(), grad_grad[0].asnumpy()[:, :, h//4:h-h//4, w//4:w-w//4], rtol=1e-2, atol=1e-4) 

def _average_precision(self, rec, prec):
        """
        calculate average precision

        Params:
        ----------
        rec : numpy.array
            cumulated recall
        prec : numpy.array
            cumulated precision
        Returns:
        ----------
        ap as float
        """
        if rec is None or prec is None:
            return np.nan

        # append sentinel values at both ends
        mrec = np.concatenate(([0.], rec, [1.]))
        mpre = np.concatenate(([0.], np.nan_to_num(prec), [0.]))

        # compute precision integration ladder
        for i in range(mpre.size - 1, 0, -1):
            mpre[i - 1] = np.maximum(mpre[i - 1], mpre[i])

        # look for recall value changes
        i = np.where(mrec[1:] != mrec[:-1])[0]

        # sum (\delta recall) * prec
        ap = np.sum((mrec[i + 1] - mrec[i]) * mpre[i + 1])
        return ap 

def test_compute_cor_loc_nans(self):
    num_gt_imgs_per_class = np.array([100, 0, 0, 1, 1], dtype=int)
    num_images_correctly_detected_per_class = np.array([10, 0, 1, 0, 0],
                                                       dtype=int)
    corloc = metrics.compute_cor_loc(num_gt_imgs_per_class,
                                     num_images_correctly_detected_per_class)
    expected_corloc = np.array([0.1, np.nan, np.nan, 0, 0], dtype=float)
    self.assertAllClose(corloc, expected_corloc) 

def test_filter_groundtruth_with_nan_box_coordinates(self):
    input_tensors = {
        fields.InputDataFields.groundtruth_boxes:
        [[np.nan, np.nan, np.nan, np.nan], [0.2, 0.4, 0.1, 0.8]],
        fields.InputDataFields.groundtruth_classes:
        [1, 2],
        fields.InputDataFields.groundtruth_is_crowd:
        [False, True],
        fields.InputDataFields.groundtruth_area:
        [100.0, 238.7]
    }

    expected_tensors = {
        fields.InputDataFields.groundtruth_boxes:
        [[0.2, 0.4, 0.1, 0.8]],
        fields.InputDataFields.groundtruth_classes:
        [2],
        fields.InputDataFields.groundtruth_is_crowd:
        [True],
        fields.InputDataFields.groundtruth_area:
        [238.7]
    }

    output_tensors = ops.filter_groundtruth_with_nan_box_coordinates(
        input_tensors)
    with self.test_session() as sess:
      output_tensors = sess.run(output_tensors)
      for key in [fields.InputDataFields.groundtruth_boxes,
                  fields.InputDataFields.groundtruth_area]:
        self.assertAllClose(expected_tensors[key], output_tensors[key])
      for key in [fields.InputDataFields.groundtruth_classes,
                  fields.InputDataFields.groundtruth_is_crowd]:
        self.assertAllEqual(expected_tensors[key], output_tensors[key]) 

def test_prune_outside_window(self):
    keypoints = tf.constant([
        [[0.25, 0.5], [0.75, 0.75]],
        [[0.5, 0.0], [1.0, 1.0]]
    ])
    window = tf.constant([0.25, 0.25, 0.75, 0.75])

    expected_keypoints = tf.constant([[[0.25, 0.5], [0.75, 0.75]],
                                      [[np.nan, np.nan], [np.nan, np.nan]]])
    output = keypoint_ops.prune_outside_window(keypoints, window)

    with self.test_session() as sess:
      output_, expected_keypoints_ = sess.run([output, expected_keypoints])
      self.assertAllClose(output_, expected_keypoints_) 

def createTestLabelScoresWithMissingScore(self):
    return tf.constant([0.5, np.nan], dtype=tf.float32) 

def expectedLabelScoresAfterThresholdingWithMissingScore(self):
    return tf.constant([np.nan], dtype=tf.float32) 

def testStrictRandomCropImageWithKeypoints(self):
    image = self.createColorfulTestImage()[0]
    boxes = self.createTestBoxes()
    labels = self.createTestLabels()
    keypoints = self.createTestKeypoints()
    with mock.patch.object(
        tf.image,
        'sample_distorted_bounding_box'
    ) as mock_sample_distorted_bounding_box:
      mock_sample_distorted_bounding_box.return_value = (
          tf.constant([6, 143, 0], dtype=tf.int32),
          tf.constant([190, 237, -1], dtype=tf.int32),
          tf.constant([[[0.03, 0.3575, 0.98, 0.95]]], dtype=tf.float32))
      (new_image, new_boxes, new_labels,
       new_keypoints) = preprocessor._strict_random_crop_image(
           image, boxes, labels, keypoints=keypoints)
      with self.test_session() as sess:
        new_image, new_boxes, new_labels, new_keypoints = sess.run([
            new_image, new_boxes, new_labels, new_keypoints])

        expected_boxes = np.array([
            [0.0, 0.0, 0.75789469, 1.0],
            [0.23157893, 0.24050637, 0.75789469, 1.0],
        ], dtype=np.float32)
        expected_keypoints = np.array([
            [[np.nan, np.nan],
             [np.nan, np.nan],
             [np.nan, np.nan]],
            [[0.38947368, 0.07173],
             [0.49473682, 0.24050637],
             [0.60000002, 0.40928277]]
        ], dtype=np.float32)
        self.assertAllEqual(new_image.shape, [190, 237, 3])
        self.assertAllClose(
            new_boxes.flatten(), expected_boxes.flatten())
        self.assertAllClose(
            new_keypoints.flatten(), expected_keypoints.flatten()) 

def __init__(self, x, y, w=None, bbox=[None]*2, k=3, ext=0, fill_value=np.nan):
        """
        See docstring for InterpolatedUnivariateSpline.
        """
        self.bounds = (min(x), max(x))
        self.fill_value = fill_value
        super(ExtrapolatingUnivariateSpline, self).__init__(x, y, w=w, bbox=bbox, k=k, ext=ext) 

def read_observed_targets(target_filename=OBS_TARGET_FNAME):
    """
    Reads the observed targets excel file into a pandas dataframe
    :param target_filename: The filename to read. Has a very specific format!
    :return:
    """
    sample_names = ['identifier', 'RA/DEC (J2000)', 'plx', 'Vmag', 'Kmag', 'vsini', 'SpT', 'configuration',
                    'Instrument',
                    'Date',
                    'Temperature', 'Velocity', 'vsini_sec', '[Fe/H]', 'Significance', 'Sens_min', 'Sens_any',
                    'Comments',
                    'Rank', 'Keck', 'VLT', 'Gemini', 'Imaging_Detecton']

    def plx_convert(s):
        try:
            return float(s)
        except ValueError:
            return np.nan

    sample = pd.read_excel(target_filename, sheetname=0, na_values=['     ~'], names=sample_names,
                           converters=dict(plx=plx_convert))
    sample = sample.reset_index(drop=True)[1:]

    # Convert everything to floats
    for col in sample.columns:
        sample[col] = pd.to_numeric(sample[col], errors='ignore')

    return sample 

def fill_for_noncomputable_vals(input_data, result_data):
    non_computable_values = np.repeat(
        np.nan, len(input_data) - len(result_data)
        )
    filled_result_data = np.append(non_computable_values, result_data)
    return filled_result_data 

def record_tabular(self, key, val):
        assert (str(key) not in self._curr_recorded)
        self._curr_recorded.append(str(key))

        if key in self._tabular:
            self._tabular[key].append(val)
        else:
            self._tabular[key] = [np.nan] * self._num_dump_tabular_calls + [val] 

def dump_tabular(self, print_func=None):
        if len(self._curr_recorded) == 0:
            return ''

        ### reset
        self._curr_recorded = list()
        self._num_dump_tabular_calls += 1

        ### make sure all same length
        for k, v in self._tabular.items():
            if len(v) == self._num_dump_tabular_calls:
                pass
            elif len(v) == self._num_dump_tabular_calls - 1:
                self._tabular[k].append(np.nan)
            else:
                raise ValueError('key {0} should not have {1} items when {2} calls have been made'.format(
                    k, len(v), self._num_dump_tabular_calls))

        ### print
        if print_func is not None:
            log_str = tabulate(sorted([(k, v[-1]) for k, v in self._tabular.items()], key=lambda kv: kv[0]))
            for line in log_str.split('\n'):
                print_func(line)

        ### write to file
        tabular_pandas = pandas.DataFrame({k: pandas.Series(v) for k, v in self._tabular.items()})
        tabular_pandas.to_csv(self._csv_path) 

def ada2PsychoPix(self, xyCoor = tuple):
        
        # check argument values
        if xyCoor is None:
            raise ValueError("No coordinate values have been specified.")
        elif not isinstance(xyCoor, tuple):
            raise TypeError("XY coordinates must be given as tuple.")
        elif isinstance(xyCoor, tuple) and len(xyCoor) is not 2: 
            raise ValueError("Wrong number of coordinate dimensions")
            
        if np.isnan(xyCoor[0]) and np.isnan(xyCoor[1]):
            psychoPix = (np.nan, np.nan)
            return psychoPix

        # convert to pixels and correct for psychopy window coordinates
        monHW = (self.win.getSizePix()[0], 
                 self.win.getSizePix()[1])
        wShift, hShift = monHW[0] / 2 , monHW[1] / 2
        psychoPix = ((((xyCoor[0]* monHW[0]) - wShift)), 
                     (((xyCoor[1] * monHW[1]) - hShift) * -1))
        # return coordinates in psychowin 'pix' units
        return psychoPix


    # function for converting from tobiis active display coordinate system in 
    # normalized coordinates where (0,0) is the upper left corner, to monitor 
    # coordinates in pix, where (0,0) is the upper left corner 

def getAvgGazePos(self):
        
        # check to see if the eyetracker is connected and turned on
        if self.eyetracker is None:
            raise ValueError("There is no eyetracker.")
        if self.tracking is False:
            raise ValueError("The eyetracker is not turned on.")
            
        # while tracking
        while True:
            # access gaze data dictionary to get gaze position tuples
            lGazeXYZ = self.gazeData['left_gaze_point_on_display_area']
            rGazeXYZ = self.gazeData['right_gaze_point_on_display_area']       
            # get 2D gaze positions for left and right eye
            xs = (lGazeXYZ[0], rGazeXYZ[0])
            ys = (lGazeXYZ[1], rGazeXYZ[1])   
            
            # if all of the axes have data from at least one eye
            if all([x != -1.0 for x in xs]) and all([y != -1.0 for y in ys]):
                # take x and y averages
                avgGazePos = np.nanmean(xs), np.nanmean(ys)
            else:
                # or if no data, hide points by showing off screen
                avgGazePos = (np.nan, np.nan)
            return self.ada2PsychoPix(avgGazePos)

                
    # function for finding the avg 3d position of subject's eyes, so that they
    # can be drawn in the virtual track box before calibration. The x and y 
    # coordinates are returned in normalized "tobii track box" units. 

def explained_variance(ypred,y):
    """
    Computes fraction of variance that ypred explains about y.
    Returns 1 - Var[y-ypred] / Var[y]

    interpretation:
        ev=0  =>  might as well have predicted zero
        ev=1  =>  perfect prediction
        ev<0  =>  worse than just predicting zero

    """
    assert y.ndim == 1 and ypred.ndim == 1
    vary = np.var(y)
    return np.nan if vary==0 else 1 - np.var(y-ypred)/vary 

def read_tb(path):
    """
    path : a tensorboard file OR a directory, where we will find all TB files
           of the form events.*
    """
    import pandas
    import numpy as np
    from glob import glob
    from collections import defaultdict
    import tensorflow as tf
    if osp.isdir(path):
        fnames = glob(osp.join(path, "events.*"))
    elif osp.basename(path).startswith("events."):
        fnames = [path]
    else:
        raise NotImplementedError("Expected tensorboard file or directory containing them. Got %s"%path)
    tag2pairs = defaultdict(list)
    maxstep = 0
    for fname in fnames:
        for summary in tf.train.summary_iterator(fname):
            if summary.step > 0:
                for v in summary.summary.value:
                    pair = (summary.step, v.simple_value)
                    tag2pairs[v.tag].append(pair)
                maxstep = max(summary.step, maxstep)
    data = np.empty((maxstep, len(tag2pairs)))
    data[:] = np.nan
    tags = sorted(tag2pairs.keys())
    for (colidx,tag) in enumerate(tags):
        pairs = tag2pairs[tag]
        for (step, value) in pairs:
            data[step-1, colidx] = value
    return pandas.DataFrame(data, columns=tags) 

def cal_er(tokenizer, pred, truth, mode='wer', ctc=False):
    # Calculate error rate of a batch
    if pred is None:
        return np.nan
    elif len(pred.shape) >= 3:
        pred = pred.argmax(dim=-1)
    er = []
    for p, t in zip(pred, truth):
        p = tokenizer.decode(p.tolist(), ignore_repeat=ctc)
        t = tokenizer.decode(t.tolist())
        if mode == 'wer':
            p = p.split(' ')
            t = t.split(' ')
        er.append(float(ed.eval(p, t))/len(t))
    return sum(er)/len(er) 

def muskingum(inflow,Qinitial,k,x,dt):
    """
    ===========================================================
     muskingum(inflow,Qinitial,k,x,dt)
    ===========================================================
    
    inputs:
    ----------
            1- inflow
            2- Qinitial initial value for outflow
            3- k travelling time (hours)
            4- x surface nonlinearity coefficient (0,0.5)
            5- dt delta t 
    Outputs:
    ----------
   """ 
    c1=(dt-2*k*x)/(2*k*(1-x)+dt)
    c2=(dt+2*k*x)/(2*k*(1-x)+dt)
    c3=(2*k*(1-x)-dt)/(2*k*(1-x)+dt)
    
#    if c1+c2+c3!=1:
#        raise("sim of c1,c2 & c3 is not 1")
    
    outflow=np.ones_like(inflow)*np.nan    
    outflow[0]=Qinitial
    
    for i in range(1,len(inflow)):
        outflow[i]=c1*inflow[i]+c2*inflow[i-1]+c3*outflow[i-1]
    
#    outflow[outflow<0]=0
    outflow=np.round(outflow,4)
    return outflow 

def euclidean_dist(a: np.ndarray, b: np.ndarray, missing=0.1):
    """ Calculate euclidean distance between two vectors of continuous variables """
    assert a.shape == b.shape
    threshold = a.shape[0] * missing
    eps = np.spacing(1)
    vec = np.power(a - b, 2)
    common = ~np.isnan(vec)  # common / not missing values
    return (sqrt(np.sum(vec[common])) + eps) / (np.sum(common) + eps) if np.sum(common) >= threshold else np.nan


# TODO update chi_squared_dist and agreement_dist to better deal with missing values 

def feature_corr_similarity(f: np.ndarray, missing: float = 0.1, method="pearson"):
    """ Generate similarity matrix from genomic assay based on correlation between samples

        Args:
            f (Numpy.ndarray): Feature matrix (n x k, where 'n' - samples, 'k' - measurements)
            missing (float): acceptable fraction of values for assessment of distance/similarity between two samples \
            (default of 0.1, means that up to 90 % of missing values is acceptable)
            method (str): either "pearson" or "spearman"
        Returns:
            sim (Numpy.ndarray): symmetric matrix describing similarity between samples (n x n)

    """
    assert method in CORR_METHODS
    # filter missing samples
    samples = np.array([i for i in range(f.shape[0]) if not np.all(np.isnan(f[i, :]))])

    # euclidean distance matrix
    sim = np.zeros((f.shape[0], f.shape[0]))
    sim[:] = np.nan

    for i in [sample for sample in range(f.shape[0]) if sample in samples]:
        for j in [sample for sample in range(i, f.shape[0]) if sample in samples]:
            sim[i, j] = corr(f[i, :], f[j, :], method=method, missing=missing)
            if i != j:
                sim[j, i] = sim[i, j]

    vals = sim[~np.isnan(sim)]
    vals[vals < 0] = 0
    sim[~np.isnan(sim)] = vals
    return sim 

def test_init():
    a0 = _create_test_adj(10)
    a1 = _create_test_adj(10)
    sample_names = ['sample_{}'.format(i) for i in range(10)]

    graph = MultiplexNet(adj_matrices=[a0, a1], node_labels=np.array(sample_names))
    assert graph.layers == 2 and graph.nodes == 10
    assert np.sum(graph.connections) == a0.size * 2  # all values in both layers

    # missing samples
    a0[:, 0], a0[0, :] = np.nan, np.nan
    graph = MultiplexNet(adj_matrices=[a0, a1], node_labels=np.array(sample_names))
    assert 0 not in graph.samples[0] and 0 in graph.samples[1]

    # incorrect adjacency matrices
    with pytest.raises(TypeError):
        MultiplexNet(adj_matrices=[[2], [4]], node_labels=np.array(sample_names))

    a2 = a0.copy()
    a2[2, 3] = np.nan
    with pytest.raises(ValueError):
        MultiplexNet(adj_matrices=[a1, a2], node_labels=np.array(sample_names))

    with pytest.raises(ValueError):
        MultiplexNet(adj_matrices=[], node_labels=np.array(sample_names))

    # different number of nodes in layers
    with pytest.raises(ValueError):
        MultiplexNet(adj_matrices=[a0[:5, :5], a1], node_labels=np.array(sample_names))

    # not normalized data
    a2 = a0.copy()
    a2[0, 1], a2[1, 0] = -1, -1
    with pytest.raises(ValueError):
        MultiplexNet(adj_matrices=[a1, a2], node_labels=np.array(sample_names))

    a2[0, 1], a2[1, 0] = 2, 2
    with pytest.raises(ValueError):
        MultiplexNet(adj_matrices=[a1, a2], node_labels=np.array(sample_names)) 

def test_feature_to_adjacency():
    f = np.random.random((10, 20))

    with pytest.raises(ValueError):
        similarity.feature_to_adjacency(f, variable_type="random")

    a = similarity.feature_to_adjacency(f, variable_type="continuous")
    assert check_matrix_symmetry(a)
    assert np.all(np.diag(a) == 1)

    # incorrect hyperparameter
    with pytest.raises(ValueError):
        similarity.feature_to_adjacency(f, variable_type="continuous", alpha=0)

    # missing samples
    f2 = f.copy()
    f2[:2, :] = np.nan
    a = similarity.feature_to_adjacency(f2, variable_type="continuous")
    assert np.all(np.isnan(a[:2, :])) and np.all(np.isnan(a[:, :2]))

    # very similar samples
    f[:5, :] = f[5:, :]
    with pytest.raises(ValueError):
        similarity.feature_to_adjacency(f, variable_type="continuous")

    a = similarity.feature_to_adjacency(f, variable_type="continuous", n=0.2)
    assert check_matrix_symmetry(a)
    assert np.all(np.diag(a) == 1) 

def test_filter_features_and_samples():
    data_vals = np.random.random((10, 20))
    data = DataFrame(data_vals.T, columns=['sample_{}'.format(i) for i in range(data_vals.shape[0])],
                     index=['feature_{}'.format(i) for i in range(data_vals.shape[1])])

    filtered = filter_features_and_samples(data)
    assert filtered.values.shape == (20, 10)

    # missing samples and features
    new_data = data.copy()
    new_data['sample_0'] = np.nan
    new_data.loc['feature_0'] = np.nan

    filtered = filter_features_and_samples(new_data)
    assert 'feature_0' not in filtered.index
    assert 'sample_0' not in filtered.columns

    # missing values sample filtering
    new_data = data.copy()
    new_data['sample_0'][1:3] = np.nan

    filtered = filter_features_and_samples(new_data)
    assert 'sample_0' in filtered.columns

    filtered = filter_features_and_samples(new_data, drop_samples=0.05)
    assert 'sample_0' not in filtered.columns

    # missing values feature filtering
    new_data = data.copy()
    new_data.loc['feature_0'][1] = np.nan

    filtered = filter_features_and_samples(new_data)
    assert 'feature_0' in filtered.index

    filtered = filter_features_and_samples(new_data, drop_features=0.05)
    assert 'feature_0' not in filtered.index 

def test_check_matrix_symmetry():
    assert utils.check_matrix_symmetry(np.array([[1, 2], [2, 1]])) is True
    assert utils.check_matrix_symmetry(np.array([[1, 2], [3, 1]])) is False
    assert utils.check_matrix_symmetry(np.array([[1, 2, 3], [1, 2, 3]])) is False

    m = np.random.rand(5, 5)
    m = (m + m.T) / 2
    assert utils.check_matrix_symmetry(m) is True

    m[3, 2] = np.nan
    m[2, 3] = np.nan
    assert utils.check_matrix_symmetry(m) is True 

def safemean(xs):
    return np.nan if len(xs) == 0 else np.mean(xs) 

def tan(x):
	''' Compute the tangent of an AutoDiff object and its derivative.

	INPUTS
	======
	x: an AutoDiff object

	'''
	try:
		# Value and derivative undefined when divisible by pi/2 but not pi
		# To make sure the asymptotes are undefined:
		if x.val%(np.pi/2)==0 and x.val%np.pi!=0:
			new_val = np.nan
			new_der = np.nan
			new_jacobian = np.nan
			warnings.warn('Undefined at value', RuntimeWarning)
		else:
			new_val = np.tan(x.val)
			new_der = x.der / (np.cos(x.val)**2.0)
			new_jacobian = x.jacobian / (np.cos(x.val)**2.0)
		return AutoDiff(new_val, new_der, x.n, 0, new_jacobian)
	except AttributeError:
		try:
			if x.Real%(np.pi/2)==0 and x.Real%np.pi!=0:
				ans = Dual(np.nan,np.nan)
				warnings.warn('Undefined at value', RuntimeWarning)
				return ans
			else:
				return Dual(np.tan(x.Real), x.Dual / (np.cos(x.Real))**2)
		except AttributeError:
			try:
				if x.Real%(np.pi/2)==0 and x.Real%np.pi!=0:
					ans = Dual(np.nan,np.nan)
					warnings.warn('Undefined at value', RuntimeWarning)
					return ans
				else:
					# return Dual(tan(x.Real), x.Dual / (cos(x.Real))**2)
					return sin(x)/cos(x)
			except AttributeError:
				if x%(np.pi/2)==0 and x%np.pi!=0:
					warnings.warn('Undefined at value', RuntimeWarning)
					return np.nan
				else:
					return np.tan(x)

#-------------------INVERSE TRIG FUNCTIONS-------------------#
# arc sin 

def sqrt(x):
	''' Compute the square root an AutoDiff object and its derivative.

	INPUTS
	======
	x: an AutoDiff object

	RETURNS
	=======
	A new AutoDiff object with calculated value and derivative.

	EXAMPLES
	========
	>>> x = AutoDiff(np.array([[5]]).T, np.array([[1]]), 1, 1)
	>>> myAutoDiff = sqrt(x)
	>>> myAutoDiff.val
	2.2360679775
	>>> myAutoDiff.der
	0.2236068

	'''
	try:
		new_val = np.sqrt(x.val)
		new_der = 0.5 * x.val ** (-0.5) * x.der
		new_jacobian = 0.5 * x.val ** (-0.5) * x.jacobian
		return AutoDiff(new_val, new_der, x.n, 0, new_jacobian)
	except AttributeError:
			try:
				if x.Real < 0.0:
					warnings.warn('Undefined at value', RuntimeWarning)
					dual=np.nan
				
				elif(x.Real==0):
					warnings.warn('Undefined at value', RuntimeWarning)
					dual = np.inf
				
				else:
					dual = 0.5 * x.Real ** (-0.5) * x.Dual
				
				real = np.sqrt(x.Real)
				return Dual(real, dual)
			except AttributeError:
				if x < 0.0:
					warnings.warn('Undefined at value', RuntimeWarning)
					return np.nan
				else:
					return np.sqrt(x)

# log base 

def run_param_est_with_subset_lambdas(self, builder_clone, end_time, subset, nfe, ncp, sigmas, solver='ipopt', ):
        """ Performs the parameter estimation with a specific subset of wavelengths.
            At the moment, this is performed as a totally new Pyomo model, based on the
            original estimation. Initialization strategies for this will be needed.

                Args:
                    builder_clone (TemplateBuidler): Template builder class of complete model
                                without the data added yet
                    end_time (float): the end time for the data and simulation
                    subset(list): list of selected wavelengths
                    nfe (int): number of finite elements
                    ncp (int): number of collocation points
                    sigmas(dict): dictionary containing the variances, as used in the ParameterEstimator class

                Returns:
                    results (Pyomo model solved): The solved pyomo model

        """
        if not isinstance(subset, (list, dict)):
            raise RuntimeError("subset must be of type list or dict!")

        if isinstance(subset, dict):
            lists1 = sorted(subset.items())
            x1, y1 = zip(*lists1)
            subset = list(x1)
        # should check whether the list contains wavelengths or not

        # This is the filter for creating the new data subset
        new_D = pd.DataFrame(np.nan, index=self._meas_times, columns=subset)
        for t in self._meas_times:
            for l in self._meas_lambdas:
                if l in subset:
                    new_D.at[t, l] = self.model.D[t, l]
        # print(new_D)
        # Now that we have a new DataFrame, we need to build the entire problem from this
        # An entire new ParameterEstimation problem should be set up, on the outside of
        # this function and class structure, from the model already developed by the user.
        new_template = construct_model_from_reduced_set(builder_clone, end_time, new_D)
        # need to put in an optional running of the variance estimator for the new
        # parameter estiamtion run, or just use the previous full model run to initialize...

        results, lof = run_param_est(new_template, nfe, ncp, sigmas, solver=solver)

        return results 

def test_num_feats(zipper, feat_subset=None, method=None):
    """
    Performs a hypothesis test to check if the value distributions of numerical features deviate signifcantly between
    the datasets. Currently t-test as a parametric and U-test as a non-parametric test are supported.

    :param zipper: Dictionary storing the feature values of the datasets in a list. Feature name is used as the key.
    :param feat_subset: A list containing feature names. If given, analysis will only be performed for the contained \
    features. If not given all features will be considered.
    :param method: Specify which statistical test should be used. "u" for Mann-Whitney-U-test, "t" for t-test and \
    "wilcoxon" for a Wilcoxon signed rank test.
    :return: dictionary storing the p_values of the analysis. Feature names are used as keys.
    """

    # if no method is specified used Mann-Whitney-U-test as standard
    if method is None:
        method = "u"

    # initialize dictionary which stores the p_values
    p_values = dict()

    if feat_subset is None:
        feat_subset = zipper.keys()

    for feat in feat_subset:  # run through all variables

        # initiate dict in dict for d1 vs d2, d2 vs d3 etc. per feature
        p_values[feat] = dict()

        for i in range(len(zipper[feat]) - 1):  # select dataset1
            for j in range(i + 1, len(zipper[feat])):  # select dataset2

                # handle the case that all values of current feature are equal across current datasets
                if _test_if_all_vals_equal(zipper[feat][i], zipper[feat][j]):
                    warnings.warn(
                        "Values of \"{}\" are the identical across the two datasets. It will be skipped.".format(feat),
                        UserWarning)
                    # delete already created dict for i, j in p_values and continue with next feature
                    del p_values[feat]
                    continue

                # only calculate score if there are values in each dataset
                if zipper[feat][i] and zipper[feat][j]:
                    # calculate u-test and return p-value
                    if method == "u":
                        stat_test_result = mannwhitneyu(zipper[feat][i], zipper[feat][j], alternative="two-sided")
                    # calculate t-test and return p-value
                    elif method == "t":
                        stat_test_result = ttest_ind(zipper[feat][i], zipper[feat][j])
                    elif method == "wilcoxon":
                        stat_test_result = wilcoxon(zipper[feat][i], zipper[feat][j])

                    p_values[feat][i + 1, j + 1] = stat_test_result.pvalue

                # if one or both sets are empty
                else:
                    p_values[feat][i + 1, j + 1] = np.nan

    return p_values 

def calc_prop_diff_conf(zipper, feat_subset):
    """
    Calculates the confidence intervals for numerical features.

    :param zipper: Zipper created from a DataCollection.
    :param feat_subset: An iterable of features for which the confidence intervals shall be calculated.
    :return:
    """
    # initialize dictionary which stores the conf_invs
    conf_invs = dict()

    for feat in feat_subset:  # run through all variables

        # initiate dict in dict for d1 vs d2, d2 vs d3 etc. per feature
        conf_invs[feat] = dict()

        for i in range(len(zipper[feat]) - 1):  # select dataset1
            for j in range(i + 1, len(zipper[feat])):  # select dataset2

                # only calculate score if there are values in each dataset
                if zipper[feat][i] and zipper[feat][j]:

                    invs = calc_prop_diff(zipper[feat][i], zipper[feat][j], feat)

                    # check for each factor of the categorical feature if 0 is in the CI or not and append it to dict
                    for factor in invs:
                        inv = invs[factor]

                        # TODO: pass the counts over to the next function such that in can be included into the result table
                        #count1 = invs[factor][0]
                        #count2 = invs[factor][1]

                        # indicator = True if 0 is not in the interval
                        flag = inv[0] >= 0.00 or inv[1] <= 0.00

                        # save interval in dict
                        try:
                            conf_invs[factor][i + 1, j + 1] = (inv, flag)
                        except KeyError:
                            conf_invs[factor] = dict()
                            conf_invs[factor][i + 1, j + 1] = (inv, flag)

                # if one or both sets are empty
                else:
                    conf_invs[feat][i + 1, j + 1] = (np.nan, np.nan)

    return conf_invs 

def parse_foul(row):
    """
    function to determine what type of foul is being commited by the player

    Input:
    row - row of nba play by play

    Output:
    foul_type - the foul type of the fould commited by the player
    """

    foul_dict = {
        0: "no foul",
        1: "personal",
        2: "shooting",
        3: "loose ball",
        4: "offensive",
        5: "inbound",
        6: "away from play",
        7: "elbow",
        8: "punching",
        9: "clear path",
        10: "double personal",
        11: "technical",
        12: "non-unsportsmanlike technical",
        13: "hanging technical",
        14: "flagrant type 1",
        15: "flagrant type 2",
        16: "double technical",
        17: "defense 3 second",
        18: "delay technical",
        19: "taunting technical",
        20: "player_control",
        21: "indirect technical",
        25: "excess timeout technical",
        26: "offensive charge",
        27: "personal block",
        28: "shooting block",
        30: "too many players technical",
    }
    try:
        if row["eventmsgtype"] == 6:
            try:
                return foul_dict[row["eventmsgactiontype"]]
            except KeyError:
                return np.nan
        return np.nan
    except KeyError:
        return np.nan 

def testRunRandomCropImageWithKeypointsOutsideCrop(self):
    image = self.createColorfulTestImage()
    boxes = self.createTestBoxes()
    labels = self.createTestLabels()
    keypoints = self.createTestKeypointsOutsideCrop()

    tensor_dict = {
        fields.InputDataFields.image: image,
        fields.InputDataFields.groundtruth_boxes: boxes,
        fields.InputDataFields.groundtruth_classes: labels,
        fields.InputDataFields.groundtruth_keypoints: keypoints
    }

    preprocessor_arg_map = preprocessor.get_default_func_arg_map(
        include_keypoints=True)

    preprocessing_options = [(preprocessor.random_crop_image, {})]

    with mock.patch.object(
        tf.image,
        'sample_distorted_bounding_box'
    ) as mock_sample_distorted_bounding_box:
      mock_sample_distorted_bounding_box.return_value = (
          tf.constant([6, 143, 0], dtype=tf.int32),
          tf.constant([190, 237, -1], dtype=tf.int32),
          tf.constant([[[0.03, 0.3575, 0.98, 0.95]]], dtype=tf.float32))
      distorted_tensor_dict = preprocessor.preprocess(
          tensor_dict, preprocessing_options, func_arg_map=preprocessor_arg_map)
      distorted_image = distorted_tensor_dict[fields.InputDataFields.image]
      distorted_boxes = distorted_tensor_dict[
          fields.InputDataFields.groundtruth_boxes]
      distorted_labels = distorted_tensor_dict[
          fields.InputDataFields.groundtruth_classes]
      distorted_keypoints = distorted_tensor_dict[
          fields.InputDataFields.groundtruth_keypoints]
      with self.test_session() as sess:
        (distorted_image_, distorted_boxes_, distorted_labels_,
         distorted_keypoints_) = sess.run(
             [distorted_image, distorted_boxes, distorted_labels,
              distorted_keypoints])

        expected_boxes = np.array([
            [0.0, 0.0, 0.75789469, 1.0],
            [0.23157893, 0.24050637, 0.75789469, 1.0],
        ], dtype=np.float32)
        expected_keypoints = np.array([
            [[np.nan, np.nan],
             [np.nan, np.nan],
             [np.nan, np.nan]],
            [[np.nan, np.nan],
             [np.nan, np.nan],
             [np.nan, np.nan]],
        ])
        self.assertAllEqual(distorted_image_.shape, [1, 190, 237, 3])
        self.assertAllEqual(distorted_labels_, [1, 2])
        self.assertAllClose(
            distorted_boxes_.flatten(), expected_boxes.flatten())
        self.assertAllClose(
            distorted_keypoints_.flatten(), expected_keypoints.flatten()) 

def _correct(self, df, cache=True):
        """
        This function is called by convert_measured_to_actual and is NOT meant to be called directly!
        It takes a pandas dataframe that all have the same star
        """

        # Group by instrument and addmode, and get the PDF for the actual temperature for each
        df = df.dropna(subset=['Tmeas'])
        if len(df) == 0:
            df['Corrected_Temperature'] = np.nan
            df['T_uperr'] = np.nan
            df['T_lowerr'] = np.nan
            return df[['Star', 'Corrected_Temperature', 'T_uperr', 'T_lowerr']].copy()
            
        groups = df.groupby(('Instrument', 'addmode'))
        Probabilities = []
        for (inst, addmode), group in groups:
            # Make a fitter instance
            d = {'instrument': inst,
                 'directory': self._caldir[inst],
                 'addmode': addmode}
            key = (inst, addmode)

            fitter = self._fitters[inst]()

            # get/set the cache
            if key in self._chainCache:
                chain, probs = self._chainCache[key]
                Tpredictions = self._predictionCache[key]
                fitter.spoof_sampler(chain, probs)
            else:
                chain = np.loadtxt(self._flatchain_format.format(**d))
                probs = np.loadtxt(self._flatlnprob_format.format(**d))
                fitter.spoof_sampler(chain, probs)

                Ta_arr = np.arange(2000, 12000, 2.0)
                Tmeas_pred = fitter.predict(Ta_arr, N=10000)
                Tpredictions = pd.DataFrame(Tmeas_pred, columns=Ta_arr)

                if cache:
                    self._chainCache[key] = (chain, probs)
                    self._predictionCache[key] = Tpredictions

            # Get the PDF (probability distribution function)
            Tmeas = group['Tmeas'].values
            Tmeas_err = group['Tmeas_err'].values
            for Tm, Tm_err in zip(Tmeas, Tmeas_err):
                temperature, probability = CCF_Systematics.get_actual_temperature(fitter, Tm, Tm_err,
                                                                                  cache=Tpredictions,
                                                                                  summarize=False)
                Probabilities.append(probability / probability.sum())

        # Multiply all the PDFs
        Prob = np.array(Probabilities).prod(axis=0)

        # Summarize the PDF (similar interface to np.percentile)
        l, m, h = integral(temperature, Prob, [0.16, 0.5, 0.84], k=0)

        # Save in a Pandas DataFrame and return
        return pd.DataFrame(data={'Star': df['Star'].values[0], 'Corrected_Temperature': m,
                                  'T_uperr': h - m, 'T_lowerr': m - l}, index=[0]) 

def weighted_mean_and_stddev(arr, weights=None, bad_value=np.nan):
    """
    Returns the weighted mean and standard deviation of an array.

    Parameters:
    ===========
    - arr:          Any list-like object (list, tuple, numpy array...)
                    The array of values.

    - weights:      A list-like object with the same shape as arr
                    The weights to apply to the values.

    - bad_value:    float
                    Return this value if the input array has zero length.

    Returns:
    ========
    The weighted mean and standard deviation of the array.
    """

    arr = np.atleast_1d(arr).astype(np.float)
    weights = np.atleast_1d(weights).astype(np.float)

    if len(arr) == 0:
        logging.warn('Zero-length array given! Mean and standard deviation are undefined')
        return bad_value, bad_value
    if weights is None:
        weights = np.ones_like(arr)

    avg = np.average(arr, weights=weights)
    if len(arr) > 1:
        var = np.average((arr - avg) ** 2, weights=weights)
        V1 = np.sum(weights)
        V2 = np.sum(weights ** 2)

        return avg, np.sqrt(var / (1 - V2 / V1 ** 2)) + np.nansum(1.0 / np.sqrt(weights))
    return avg, 1.0 / np.sqrt(weights[0])


# ############################################################
#  The following functions are really only useful for
#  my project.
############################################################# 

def multiGeomHandler(multi_geometry, coord_type, geom_type):
    """
    ===================================================================
       multiGeomHandler(multi_geometry, coord_type, geom_type)
    ===================================================================
    Function for handling multi-geometries. Can be MultiPoint, MultiLineString or MultiPolygon.
    Returns a list of coordinates where all parts of Multi-geometries are merged into a single list.
    Individual geometries are separated with np.nan which is how Bokeh wants them.
    # Bokeh documentation regarding the Multi-geometry issues can be found here (it is an open issue)
    # https://github.com/bokeh/bokeh/issues/2321
    
    inputs:
    ----------
        1- multi_geometry (geometry)
         the geometry of a shpefile
        2- coord_type (string)
         "string" either "x" or "y"
        3- geom_type (string)
            "MultiPoint" or "MultiLineString" or "MultiPolygon"
    outpus:
    ----------
        1-array: 
         contains x coordinates or y coordinates of all edges of the shapefile    
    """
    for i,part in enumerate(multi_geometry):
        # On the first part of the Multi-geometry initialize the coord_array (np.array)
        if i ==0:
            if geom_type=="MultiPoint":
                coord_arrays= getPointCoords(part, coord_type)#,np.nan)
            elif geom_type=="MultiLineString":
#                coord_arrays= np.append(getLineCoords(part,coord_type))#,np.nan)
                coord_arrays= getLineCoords(part,coord_type)
            elif geom_type=="MultiPolygon":
                coord_arrays= 999 #getPolyCoords(part,coord_type)#,np.nan)
        else:
            if geom_type=="MultiPoint":
                coord_arrays= np.concatenate([coord_arrays,getPointCoords(part, coord_type)]) #,np.nan
            elif geom_type=="MultiLineString":
                coord_arrays= np.concatenate([coord_arrays,getLineCoords(part,coord_type)]) #,np.nan
            elif geom_type=="MultiPolygon":
                coord_arrays= 999 #np.concatenate([coord_arrays,getPolyCoords(part,coord_type)]) #,np.nan
        # return the coordinates 
        return coord_arrays 

def multiGeomHandler(multi_geometry, coord_type, geom_type):
    """
    # =============================================================================
    #     multiGeomHandler(multi_geometry, coord_type, geom_type)
    # =============================================================================
    Function for handling multi-geometries. Can be MultiPoint, MultiLineString or MultiPolygon.
    Returns a list of coordinates where all parts of Multi-geometries are merged into a single list.
    Individual geometries are separated with np.nan which is how Bokeh wants them.
    # Bokeh documentation regarding the Multi-geometry issues can be found here (it is an open issue)
    # https://github.com/bokeh/bokeh/issues/2321
    
    inputs:
        1- multi_geometry (geometry)
         the geometry of a shpefile
        2- coord_type (string)
         "string" either "x" or "y"
        3- geom_type (string)
            "MultiPoint" or "MultiLineString" or "MultiPolygon"
    outpus:
        1-array: 
         contains x coordinates or y coordinates of all edges of the shapefile    
    """
    if geom_type=="MultiPoint" or geom_type== "MultiLineString":
        for i,part in enumerate(multi_geometry):
            # On the first part of the Multi-geometry initialize the coord_array (np.array)
            if i ==0:
                if geom_type=="MultiPoint":
                    coord_arrays= getPointCoords(part, coord_type)#,np.nan)
                elif geom_type=="MultiLineString":
    #                coord_arrays= np.append(getLineCoords(part,coord_type))#,np.nan)
                    coord_arrays= getLineCoords(part,coord_type)
            else:
                if geom_type=="MultiPoint":
                    coord_arrays= np.concatenate([coord_arrays,getPointCoords(part, coord_type)]) #,np.nan
                elif geom_type=="MultiLineString":
                    coord_arrays= np.concatenate([coord_arrays,getLineCoords(part,coord_type)]) #,np.nan
            
    elif geom_type=="MultiPolygon":
        if i ==0:
#            coord_arrays= getPolyCoords(part,coord_type)#,np.nan)
            multi_2_single=explode(multi_geometry)
            for j in range(len(multi_2_single)):
                if j ==0:
                    coord_arrays= getPolyCoords(multi_2_single[j],coord_type)#,np.nan)
                else: 
                    coord_arrays= np.concatenate([coord_arrays,getPolyCoords(multi_2_single[j],coord_type)]) #,np.nan
        else:
            # explode the multipolygon into polygons 
            multi_2_single=explode(part)
            for j in range(len(multi_2_single)):
                coord_arrays= np.concatenate([coord_arrays,getPolyCoords(multi_2_single[j],coord_type)]) #,np.nan
        # return the coordinates 
        return coord_arrays 

def get_targets(dem):
    """
    ===================================================================
        get_targets(dem)
    ===================================================================
    Returns the centres of the interpolation targets
    
    Parameters
    ----------
    dem : gdal_Dataset
        Get the data from the gdal datasetof the DEM
    
    Returns
    -------
    
    coords : nd_array [nxm - nan, 2]
        Array with a list of the coordinates to be interpolated, without the Nan
    
    mat_range : nd_array [n, m]
        Array with all the centres of cells in the domain of the DEM (rectangular)
    
    """
    # Getting data for the whole grid
    x_init, xx_span, xy_span, y_init, yy_span, yx_span = dem.GetGeoTransform()
    shape_dem = dem.ReadAsArray().shape
    
    # Getting data of the mask
    no_val = dem.GetRasterBand(1).GetNoDataValue()
    mask = dem.ReadAsArray()
    
    # Adding 0.5 to get the centre
    x = np.array([x_init + xx_span*(i+0.5) for i in range(shape_dem[0])])
    y = np.array([y_init + yy_span*(i+0.5) for i in range(shape_dem[1])])
    #mat_range = np.array([[(xi, yi) for xi in x] for yi in y])
    mat_range = [[(xi, yi) for xi in x] for yi in y]
    
    # applying the mask
    coords = []
    for i in range(len(x)):
        for j in range(len(y)):
            if mask[j, i] != no_val:
                coords.append(mat_range[j][i])
                #mat_range[j, i, :] = [np.nan, np.nan]

    return np.array(coords), np.array(mat_range) 

def reproject_points(lat,lng,from_epsg=4326,to_epsg=3857):
    """
    =====================================================================
      reproject_points(lat, lng, from_epsg=4326,to_epsg=3857)
    =====================================================================
    this function change the projection of the coordinates from a coordinate system
    to another (default from GCS to web mercator used by google maps)
    
    Inputs:
    ----------
        1- lat: 
            list of latitudes of the points 
        2- lng:
            list of longitude of the points
        3- from_epsg:
            integer reference number to the projection of the points (https://epsg.io/)
        4- to_epsg:
            integer reference number to the new projection of the points (https://epsg.io/)
    
    outputs:
    ----------
        1-x:
            list of x coordinates of the points 
        2-y:
            list of y coordinates of the points 
    
    Ex:
    ----------
        # from web mercator to GCS WGS64:
        x=[-8418583.96378159, -8404716.499972705], y=[529374.3212213353, 529374.3212213353]
        from_epsg = 3857, to_epsg = 4326
        longs, lats=reproject_points(y,x,from_epsg="3857", to_epsg="4326")
    """
    from pyproj import Proj, transform
    from_epsg="epsg:"+str(from_epsg)
    inproj = Proj(init=from_epsg) # GCS geographic coordinate system
    to_epsg="epsg:"+str(to_epsg)
    outproj=Proj(init=to_epsg) # WGS84 web mercator 
    
    x=np.ones(len(lat))*np.nan
    y=np.ones(len(lat))*np.nan
    
    for i in range(len(lat)):
        x[i],y[i]=transform(inproj,outproj,lng[i],lat[i])

    return x,y 

def _get_targets(dem):
    '''
    Returns the centres of the interpolation targets
    
    Parameters
    ----------
    dem : gdal_Dataset
        Get the data from the gdal datasetof the DEM
    
    Returns
    -------
    
    coords : nd_array [nxm - nan, 2]
        Array with a list of the coordinates to be interpolated, without the Nan
    
    mat_range : nd_array [n, m]
        Array with all the centres of cells in the domain of the DEM (rectangular)
    
    '''
    # Getting data for the whole grid
    x_init, xx_span, xy_span, y_init, yy_span, yx_span = dem.GetGeoTransform()
    shape_dem = dem.ReadAsArray().shape
    
    # Getting data of the mask
    no_val = dem.GetRasterBand(1).GetNoDataValue()
    mask = dem.ReadAsArray()
    
    # Adding 0.5 to get the centre
    x = np.array([x_init + xx_span*(i+0.5) for i in range(shape_dem[0])])
    y = np.array([y_init + yy_span*(i+0.5) for i in range(shape_dem[1])])
    #mat_range = np.array([[(xi, yi) for xi in x] for yi in y])
    mat_range = [[(xi, yi) for xi in x] for yi in y]
    
    # applying the mask
    coords = []
    for i in range(len(x)):
        for j in range(len(y)):
            if mask[j, i] != no_val:
                coords.append(mat_range[j][i])
                #mat_range[j, i, :] = [np.nan, np.nan]

    return np.array(coords), np.array(mat_range) 

def FD_array(self):
        """
        =============================================================
          FD_array(FD)
        =============================================================
        this function takes flow firection codes for the 8 directions (1,2,4,8,16,32,64,128)
        and return the indeces of the upstream cells
        
        inputs:
            1- FD:
                flow direction raster
        output:
            1-flow_direction_cell:
                [numpy array] with the same dimensions of the raster and 2 layers
                first layer for row index and second row for column index
        """
        FD=self.raster
        no_val = FD.GetRasterBand(1).GetNoDataValue()
        no_columns=FD.RasterXSize
        no_rows=FD.RasterYSize
    
        flow_direction=np.float32(FD.ReadAsArray())
        flow_direction[flow_direction==no_val]=np.nan
    
        flow_direction_cell=np.ones((no_rows,no_columns,2))*np.nan
    
        for i in range(no_rows):
            for j in range(no_columns):
                if flow_direction[i,j]==1:
                    flow_direction_cell[i,j,0]=i  # index of the row 
                    flow_direction_cell[i,j,1]=j+1 # index of the column
                elif flow_direction[i,j]==128:
                    flow_direction_cell[i,j,0]=i-1
                    flow_direction_cell[i,j,1]=j+1
                elif flow_direction[i,j]==64:
                    flow_direction_cell[i,j,0]=i-1
                    flow_direction_cell[i,j,1]=j   
                elif flow_direction[i,j]==32:
                    flow_direction_cell[i,j,0]=i-1
                    flow_direction_cell[i,j,1]=j-1
                elif flow_direction[i,j]==16:
                    flow_direction_cell[i,j,0]=i
                    flow_direction_cell[i,j,1]=j-1
                elif flow_direction[i,j]==8:
                    flow_direction_cell[i,j,0]=i+1
                    flow_direction_cell[i,j,1]=j-1
                elif flow_direction[i,j]==4:
                    flow_direction_cell[i,j,0]=i+1
                    flow_direction_cell[i,j,1]=j
                elif flow_direction[i,j]==2:
                    flow_direction_cell[i,j,0]=i+1
                    flow_direction_cell[i,j,1]=j+1
        
        return flow_direction_cell

#    def FlowAccTable() 

def multiGeomHandler(multi_geometry, coord_type, geom_type):
    """
    ===================================================================
       multiGeomHandler(multi_geometry, coord_type, geom_type)
    ===================================================================
    Function for handling multi-geometries. Can be MultiPoint, MultiLineString or MultiPolygon.
    Returns a list of coordinates where all parts of Multi-geometries are merged into a single list.
    Individual geometries are separated with np.nan which is how Bokeh wants them.
    # Bokeh documentation regarding the Multi-geometry issues can be found here (it is an open issue)
    # https://github.com/bokeh/bokeh/issues/2321
    
    inputs:
    ----------
        1- multi_geometry (geometry)
         the geometry of a shpefile
        2- coord_type (string)
         "string" either "x" or "y"
        3- geom_type (string)
            "MultiPoint" or "MultiLineString" or "MultiPolygon"
    outpus:
    ----------
        1-array: 
         contains x coordinates or y coordinates of all edges of the shapefile    
    """
    for i,part in enumerate(multi_geometry):
        # On the first part of the Multi-geometry initialize the coord_array (np.array)
        if i ==0:
            if geom_type=="MultiPoint":
                coord_arrays= getPointCoords(part, coord_type)#,np.nan)
            elif geom_type=="MultiLineString":
#                coord_arrays= np.append(getLineCoords(part,coord_type))#,np.nan)
                coord_arrays= getLineCoords(part,coord_type)
            elif geom_type=="MultiPolygon":
                coord_arrays= 999 #getPolyCoords(part,coord_type)#,np.nan)
        else:
            if geom_type=="MultiPoint":
                coord_arrays= np.concatenate([coord_arrays,getPointCoords(part, coord_type)]) #,np.nan
            elif geom_type=="MultiLineString":
                coord_arrays= np.concatenate([coord_arrays,getLineCoords(part,coord_type)]) #,np.nan
            elif geom_type=="MultiPolygon":
                coord_arrays= 999 #np.concatenate([coord_arrays,getPolyCoords(part,coord_type)]) #,np.nan
        # return the coordinates 
        return coord_arrays 

def muskingum(inflow,Qinitial,k,x,dt):
    """
    ===========================================================
     muskingum(inflow,Qinitial,k,x,dt)
    ===========================================================
    
    inputs:
    ----------
        1-inflow:
            [numpy array] time series of inflow hydrograph
        2-Qinitial:
            [numeric] initial value for outflow
        3-k:
            [numeric] travelling time (hours)
        4-x:
            [numeric] surface nonlinearity coefficient (0,0.5)
        5-dt:
            [numeric] delta t
        
    Outputs:
    ----------
        1-outflow:
            [numpy array] time series of routed hydrograph
    
    Examples:
    ----------
    pars[10]=k
    pars[11]=x
    p2[0]=1  # hourly time step
    q_routed = Routing.muskingum(q_uz,q_uz[0],pars[10],pars[11],p2[0]) 
   """ 
   
    c1=(dt-2*k*x)/(2*k*(1-x)+dt)
    c2=(dt+2*k*x)/(2*k*(1-x)+dt)
    c3=(2*k*(1-x)-dt)/(2*k*(1-x)+dt)
    
#    if c1+c2+c3!=1:
#        raise("sim of c1,c2 & c3 is not 1")
    
    outflow=np.ones_like(inflow)*np.nan    
    outflow[0]=Qinitial
    
    for i in range(1,len(inflow)):
        outflow[i]=c1*inflow[i]+c2*inflow[i-1]+c3*outflow[i-1]
    
    outflow=np.round(outflow,4)
    
    return outflow 

def multiGeomHandler(multi_geometry, coord_type, geom_type):
    """
    ===================================================================
       multiGeomHandler(multi_geometry, coord_type, geom_type)
    ===================================================================
    Function for handling multi-geometries. Can be MultiPoint, MultiLineString or MultiPolygon.
    Returns a list of coordinates where all parts of Multi-geometries are merged into a single list.
    Individual geometries are separated with np.nan which is how Bokeh wants them.
    # Bokeh documentation regarding the Multi-geometry issues can be found here (it is an open issue)
    # https://github.com/bokeh/bokeh/issues/2321
    
    inputs:
    ----------
        1- multi_geometry (geometry)
         the geometry of a shpefile
        2- coord_type (string)
         "string" either "x" or "y"
        3- geom_type (string)
            "MultiPoint" or "MultiLineString" or "MultiPolygon"
    outpus:
    ----------
        1-array: 
         contains x coordinates or y coordinates of all edges of the shapefile    
    """
    for i,part in enumerate(multi_geometry):
        # On the first part of the Multi-geometry initialize the coord_array (np.array)
        if i ==0:
            if geom_type=="MultiPoint":
                coord_arrays= getPointCoords(part, coord_type)#,np.nan)
            elif geom_type=="MultiLineString":
#                coord_arrays= np.append(getLineCoords(part,coord_type))#,np.nan)
                coord_arrays= getLineCoords(part,coord_type)
            elif geom_type=="MultiPolygon":
                coord_arrays= 999 #getPolyCoords(part,coord_type)#,np.nan)
        else:
            if geom_type=="MultiPoint":
                coord_arrays= np.concatenate([coord_arrays,getPointCoords(part, coord_type)]) #,np.nan
            elif geom_type=="MultiLineString":
                coord_arrays= np.concatenate([coord_arrays,getLineCoords(part,coord_type)]) #,np.nan
            elif geom_type=="MultiPolygon":
                coord_arrays= 999 #np.concatenate([coord_arrays,getPolyCoords(part,coord_type)]) #,np.nan
        # return the coordinates 
        return coord_arrays 

def ReprojectPoints(lat,lng,from_epsg=4326,to_epsg=3857):
    """
    =====================================================================
      reproject_points(lat, lng, from_epsg=4326,to_epsg=3857)
    =====================================================================
    this function change the projection of the coordinates from a coordinate system
    to another (default from GCS to web mercator used by google maps)
    
    Inputs:
    ----------
        1- lat: 
            list of latitudes of the points 
        2- lng:
            list of longitude of the points
        3- from_epsg:
            integer reference number to the projection of the points (https://epsg.io/)
        4- to_epsg:
            integer reference number to the new projection of the points (https://epsg.io/)
    
    outputs:
    ----------
        1-x:
            list of x coordinates of the points 
        2-y:
            list of y coordinates of the points 
    
    Ex:
    ----------
        # from web mercator to GCS WGS64:
        x=[-8418583.96378159, -8404716.499972705], y=[529374.3212213353, 529374.3212213353]
        from_epsg = 3857, to_epsg = 4326
        longs, lats=reproject_points(y,x,from_epsg="3857", to_epsg="4326")
    """
    from pyproj import Proj, transform
    from_epsg="epsg:"+str(from_epsg)
    inproj = Proj(init=from_epsg) # GCS geographic coordinate system
    to_epsg="epsg:"+str(to_epsg)
    outproj=Proj(init=to_epsg) # WGS84 web mercator 
    
    x=np.ones(len(lat))*np.nan
    y=np.ones(len(lat))*np.nan
    
    for i in range(len(lat)):
        x[i],y[i]=transform(inproj,outproj,lng[i],lat[i])

    return x,y 

def multiGeomHandler(multi_geometry, coord_type, geom_type):
    """
    # =============================================================================
    #     multiGeomHandler(multi_geometry, coord_type, geom_type)
    # =============================================================================
    Function for handling multi-geometries. Can be MultiPoint, MultiLineString or MultiPolygon.
    Returns a list of coordinates where all parts of Multi-geometries are merged into a single list.
    Individual geometries are separated with np.nan which is how Bokeh wants them.
    # Bokeh documentation regarding the Multi-geometry issues can be found here (it is an open issue)
    # https://github.com/bokeh/bokeh/issues/2321
    
    inputs:
        1- multi_geometry (geometry)
         the geometry of a shpefile
        2- coord_type (string)
         "string" either "x" or "y"
        3- geom_type (string)
            "MultiPoint" or "MultiLineString" or "MultiPolygon"
    outpus:
        1-array: 
         contains x coordinates or y coordinates of all edges of the shapefile    
    """
    if geom_type=="MultiPoint" or geom_type== "MultiLineString":
        for i,part in enumerate(multi_geometry):
            # On the first part of the Multi-geometry initialize the coord_array (np.array)
            if i ==0:
                if geom_type=="MultiPoint":
                    coord_arrays= GetPointCoords(part, coord_type)#,np.nan)
                elif geom_type=="MultiLineString":
    #                coord_arrays= np.append(getLineCoords(part,coord_type))#,np.nan)
                    coord_arrays= GetLineCoords(part,coord_type)
            else:
                if geom_type=="MultiPoint":
                    coord_arrays= np.concatenate([coord_arrays,GetPointCoords(part, coord_type)]) #,np.nan
                elif geom_type=="MultiLineString":
                    coord_arrays= np.concatenate([coord_arrays,GetLineCoords(part,coord_type)]) #,np.nan
            
    elif geom_type=="MultiPolygon":
        if i ==0:
#            coord_arrays= getPolyCoords(part,coord_type)#,np.nan)
            multi_2_single=Explode(multi_geometry)
            for j in range(len(multi_2_single)):
                if j ==0:
                    coord_arrays= GetPolyCoords(multi_2_single[j],coord_type)#,np.nan)
                else: 
                    coord_arrays= np.concatenate([coord_arrays,GetPolyCoords(multi_2_single[j],coord_type)]) #,np.nan
        else:
            # explode the multipolygon into polygons 
            multi_2_single=Explode(part)
            for j in range(len(multi_2_single)):
                coord_arrays= np.concatenate([coord_arrays,GetPolyCoords(multi_2_single[j],coord_type)]) #,np.nan
        # return the coordinates 
        return coord_arrays 

def muskingum(inflow,Qinitial,k,x,dt):
    """
    ===========================================================
     muskingum(inflow,Qinitial,k,x,dt)
    ===========================================================
    
    inputs:
    ----------
        1-inflow:
            [numpy array] time series of inflow hydrograph
        2-Qinitial:
            [numeric] initial value for outflow
        3-k:
            [numeric] travelling time (hours)
        4-x:
            [numeric] surface nonlinearity coefficient (0,0.5)
        5-dt:
            [numeric] delta t
        
    Outputs:
    ----------
        1-outflow:
            [numpy array] time series of routed hydrograph
    
    Examples:
    ----------
    pars[10]=k
    pars[11]=x
    p2[0]=1  # hourly time step
    q_routed = Routing.muskingum(q_uz,q_uz[0],pars[10],pars[11],p2[0]) 
   """ 
   
    c1=(dt-2*k*x)/(2*k*(1-x)+dt)
    c2=(dt+2*k*x)/(2*k*(1-x)+dt)
    c3=(2*k*(1-x)-dt)/(2*k*(1-x)+dt)
    
#    if c1+c2+c3!=1:
#        raise("sim of c1,c2 & c3 is not 1")
    
    outflow=np.ones_like(inflow)*np.nan    
    outflow[0]=Qinitial
    
    for i in range(1,len(inflow)):
        outflow[i]=c1*inflow[i]+c2*inflow[i-1]+c3*outflow[i-1]
    
    outflow=np.round(outflow,4)
    
    return outflow