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 query(self, coords, order=1):
        """
        Returns the map value at the specified location(s) on the sky.

        Args:
            coords (`astropy.coordinates.SkyCoord`): The coordinates to query.
            order (Optional[int]): Interpolation order to use. Defaults to `1`,
                for linear interpolation.

        Returns:
            A float array containing the map value at every input coordinate.
            The shape of the output will be the same as the shape of the
            coordinates stored by `coords`.
        """
        out = np.full(len(coords.l.deg), np.nan, dtype='f4')

        for pole in self.poles:
            m = (coords.b.deg >= 0) if pole == 'ngp' else (coords.b.deg < 0)

            if np.any(m):
                data, w = self._data[pole]
                x, y = w.wcs_world2pix(coords.l.deg[m], coords.b.deg[m], 0)
                out[m] = map_coordinates(data, [y, x], order=order, mode='nearest')

        return out 

def test_select_confounds_error(confounds_file, tmp_path):
    import pandas as pd
    import numpy as np

    confounds_df = pd.read_csv(str(confounds_file), sep='\t', na_values='n/a')

    confounds_df['white_matter'][0] = np.nan

    conf_file = tmp_path / "confounds.tsv"

    confounds_df.to_csv(str(conf_file), index=False, sep='\t', na_rep='n/a')

    with pytest.raises(ValueError) as val_err:
        _select_confounds(str(conf_file), ['white_matter', 'csf'])

    assert "The selected confounds contain nans" in str(val_err.value) 

def wnba_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
    """

    try:
        if row["etype"] == 6:
            try:
                return foul_dict[row["mtype"]]
            except KeyError:
                return np.nan
        return np.nan
    except KeyError:
        return np.nan 

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
    """

    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 wnba_shot_types(row):
    """
    function to parse what type of shot is being taken

    Inputs:
    row - pandas row of play by play dataframe

    Outputs:
    shot_type - returns a shot type of the values hook, jump, layup, dunk, tip
    """
    try:
        if row["etype"] in [1, 2, 3]:
            return SHOT_DICT[row["etype"]][row["mtype"]]
        else:
            return np.nan
    except KeyError:
        return np.nan 

def parse_shot_types(row):
    """
    function to parse what type of shot is being taken

    Inputs:
    row - pandas row of play by play dataframe

    Outputs:
    shot_type - returns a shot type of the values hook, jump, layup, dunk, tip
    """
    try:
        if row["eventmsgtype"] in [1, 2, 3]:
            return SHOT_DICT[row["eventmsgtype"]][row["eventmsgactiontype"]]
        else:
            return np.nan
    except KeyError:
        return np.nan 

def apply_cmap(zs, cmap, vmin=None, vmax=None, unit=None, logrescale=False):
    '''
    apply_cmap(z, cmap) applies the given cmap to the values in z; if vmin and/or vmax are passed,
      they are used to scale z.

    Note that this function can automatically rescale data into log-space if the colormap is a
    neuropythy log-space colormap such as log_eccentricity. To enable this behaviour use the
    optional argument logrescale=True.
    '''
    zs = pimms.mag(zs) if unit is None else pimms.mag(zs, unit)
    zs = np.asarray(zs, dtype='float')
    if pimms.is_str(cmap): cmap = matplotlib.cm.get_cmap(cmap)
    if logrescale:
        if vmin is None: vmin = np.log(np.nanmin(zs))
        if vmax is None: vmax = np.log(np.nanmax(zs))
        mn = np.exp(vmin)
        u = zdivide(nanlog(zs + mn) - vmin, vmax - vmin, null=np.nan)
    else:        
        if vmin is None: vmin = np.nanmin(zs)
        if vmax is None: vmax = np.nanmax(zs)
        u = zdivide(zs - vmin, vmax - vmin, null=np.nan)
    u[np.isnan(u)] = -np.inf
    return cmap(u) 

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 plot_results(start=0, stop=0):  # from utils.utils import *; plot_results()
    # Plot training results files 'results*.txt'
    fig, ax = plt.subplots(2, 5, figsize=(14, 7))
    ax = ax.ravel()
    s = ['GIoU', 'Objectness', 'Classification', 'Precision', 'Recall',
         'val GIoU', 'val Objectness', 'val Classification', 'mAP', 'F1']
    for f in sorted(glob.glob('results*.txt') + glob.glob('../../Downloads/results*.txt')):
        results = np.loadtxt(f, usecols=[2, 3, 4, 8, 9, 12, 13, 14, 10, 11], ndmin=2).T
        n = results.shape[1]  # number of rows
        x = range(start, min(stop, n) if stop else n)
        for i in range(10):
            y = results[i, x]
            if i in [0, 1, 2, 5, 6, 7]:
                y[y == 0] = np.nan  # dont show zero loss values
            ax[i].plot(x, y, marker='.', label=f.replace('.txt', ''))
            ax[i].set_title(s[i])
            if i in [5, 6, 7]:  # share train and val loss y axes
                ax[i].get_shared_y_axes().join(ax[i], ax[i - 5])

    fig.tight_layout()
    ax[1].legend()
    fig.savefig('results.png', dpi=200) 

def plot_results_overlay(start=0, stop=0):  # from utils.utils import *; plot_results_overlay()
    # Plot training results files 'results*.txt', overlaying train and val losses
    s = ['train', 'train', 'train', 'Precision', 'mAP', 'val', 'val', 'val', 'Recall', 'F1']  # legends
    t = ['GIoU', 'Objectness', 'Classification', 'P-R', 'mAP-F1']  # titles
    for f in sorted(glob.glob('results*.txt') + glob.glob('../../Downloads/results*.txt')):
        results = np.loadtxt(f, usecols=[2, 3, 4, 8, 9, 12, 13, 14, 10, 11], ndmin=2).T
        n = results.shape[1]  # number of rows
        x = range(start, min(stop, n) if stop else n)
        fig, ax = plt.subplots(1, 5, figsize=(14, 3.5))
        ax = ax.ravel()
        for i in range(5):
            for j in [i, i + 5]:
                y = results[j, x]
                if i in [0, 1, 2]:
                    y[y == 0] = np.nan  # dont show zero loss values
                ax[i].plot(x, y, marker='.', label=s[j])
            ax[i].set_title(t[i])
            ax[i].legend()
            ax[i].set_ylabel(f) if i == 0 else None  # add filename
        fig.tight_layout()
        fig.savefig(f.replace('.txt', '.png'), dpi=200) 

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 OutlierDetection(CMat, s):
    n = np.amax(s)
    _, N = CMat.shape
    OutlierIndx = list()
    FailCnt = 0
    Fail = False

    for i in range(0, N):
        c = CMat[:, i]
        if np.sum(np.isnan(c)) >= 1:
            OutlierIndx.append(i)
            FailCnt += 1
    sc = s.astype(float)
    sc[OutlierIndx] = np.nan
    CMatC = CMat.astype(float)
    CMatC[OutlierIndx, :] = np.nan
    CMatC[:, OutlierIndx] = np.nan
    OutlierIndx = OutlierIndx

    if FailCnt > (N - n):
        CMatC = np.nan
        sc = np.nan
        Fail = True
    return CMatC, sc, OutlierIndx, Fail 

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 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 query(self, coords, component='mean'):
        """
        Returns the extinction density (in e-foldings / kpc, in Gaia G-band)
        at the given coordinates.

        Args:
            coords (:obj:`astropy.coordinates.SkyCoord`): Coordinates at which
                to query the extinction. Must be 3D (i.e., include distance
                information).
            component (str): Which component to return. Allowable values are
                'mean' (for the mean extinction density) and 'std' (for the
                standard deviation of extinction density). Defaults to 'mean'.

        Returns:
            The extinction density, in units of e-foldings / pc, as either a
            numpy array or float, with the same shape as the input
            :obj:`coords`.
        """
        idx,mask = self._coords2idx(coords)

        v = self._data[component][idx[0], idx[1], idx[2]]
        
        if np.any(mask):
            # Set extinction to NaN for out-of-bounds (x, y, z)
            v[mask] = np.nan

        return v 

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 confounds_file(deriv_dir, preproc_file,
                   deriv_regressor_fname=deriv_regressor_fname):
    confounds_file = deriv_dir.ensure(deriv_regressor_fname)
    confound_dict = {}
    tp = nib.load(str(preproc_file)).shape[-1]
    ix = np.arange(tp)
    # csf
    confound_dict['csf'] = np.cos(2*np.pi*ix*(50/tp)) * 0.1
    # white matter
    confound_dict['white_matter'] = np.sin(2*np.pi*ix*(22/tp)) * 0.1
    # framewise_displacement
    confound_dict['framewise_displacement'] = np.random.random_sample(tp)
    confound_dict['framewise_displacement'][0] = np.nan
    # motion outliers
    for motion_outlier in range(0, 5):
        mo_name = 'motion_outlier0{}'.format(motion_outlier)
        confound_dict[mo_name] = np.zeros(tp)
        confound_dict[mo_name][motion_outlier] = 1
    # derivatives
    derive1 = [
        'csf_derivative1',
        'csf_derivative1_power2',
        'global_signal_derivative1_power2',
        'trans_x_derivative1',
        'trans_y_derivative1',
        'trans_z_derivative1',
        'trans_x_derivative1_power2',
        'trans_y_derivative1_power2',
        'trans_z_derivative1_power2',
    ]
    for d in derive1:
        confound_dict[d] = np.random.random_sample(tp)
        confound_dict[d][0] = np.nan

    # transformations
    for dir in ["trans_x", "trans_y", "trans_z"]:
        confound_dict[dir] = np.random.random_sample(tp)

    confounds_df = pd.DataFrame(confound_dict)
    confounds_df.to_csv(str(confounds_file), index=False, sep='\t', na_rep='n/a')
    return confounds_file 

def scale_for_cmap(cmap, x, vmin=Ellipsis, vmax=Ellipsis, unit=Ellipsis):
    '''
    scale_for_cmap(cmap, x) yields the values in x rescaled to be appropriate for the given
      colormap cmap. The cmap must be the name of a colormap or a colormap object.

    For a given cmap argument, if the object is a colormap itself, it is treated as cmap.name.
    If the cmap names a colormap known to neuropythy, neuropythy will rescale the values in x
    according to a heuristic.
    '''
    import matplotlib as mpl
    if isinstance(cmap, mpl.colors.Colormap): cmap = cmap.name
    (name, cm) = (None, None)
    if cmap not in colormaps:
        for (k,v) in six.iteritems(colormaps):
            if cmap in k:
                (name, cm) = (k, v)
                break
    else: (name, cm) = (cmap, colormaps[cmap])
    if cm is not None:
        cm = cm if len(cm) == 3 else (cm + (None,))
        (cm, (mn,mx), uu) = cm
        if vmin is Ellipsis: vmin = mn
        if vmax is Ellipsis: vmax = mx
        if unit is Ellipsis: unit = uu
    if vmin is Ellipsis: vmin = None
    if vmax is Ellipsis: vmax = None
    if unit is Ellipsis: unit = None
    x = pimms.mag(x) if unit is None else pimms.mag(x, unit)
    if name is not None and name.startswith('log_'):
        emn = np.exp(vmin)
        x = np.log(x + emn)
    vmin = np.nanmin(x) if vmin is None else vmin
    vmax = np.nanmax(x) if vmax is None else vmax
    return zdivide(x - vmin, vmax - vmin, null=np.nan) 

def cifti_split(cii, label=('lh', 'rh', 'rest'), subject=None, hemi=None, null=np.nan):
    '''
    cifti_split(cii, label) yields the rows or columns of the given cifti file that correspond to
      the given label (see below).
    cifti_split(cii) is equivalent to cifti_split(cii, ('lh', 'rh', 'rest')).

    The label argument may be any of the following:
      * a valid CIFTI label name such as 'CIFTI_STRUCTURE_CEREBELLUM' or
        'CIFTI_STRUCTURE_CORTEX_LEFT';
      * an abbreviated name such as 'cerebellum' for 'CIFTI_STRUCTURE_CEREBELLUM'.
      * the abbreviations 'lh' and 'rh' which stand for 'CIFTI_STRUCTURE_CORTEX_LEFT' and 
        'CIFTI_STRUCTURE_CORTEX_RIGHT';
      * the special keyword 'rest', which represents all the rows/columns not collected by any other
        instruction ('rest', by itself, results in the whole matrix being returned); or
      * A tuple of the above, indicating that each of the items listed should be returned
        sequentially in a tuple.

    The following optional arguments may be given:
      * subject (default: None) may specify the subject
      * hemi (default: None) can specify the hemisphere object that 
    '''
    dat = np.asanyarray(cii.dataobj if is_image(cii) else cii)
    n = dat.shape[-1]
    atlas = cifti_split._size_data.get(n, None)
    if atlas is None: raise ValueError('cannot split cifti with size %d' % n)
    if atlas not in cifti_split._atlas_cache:
        patt = os.path.join('data', 'fs_LR', '%s.atlasroi.%dk_fs_LR.shape.gii')
        lgii = nib.load(os.path.join(library_path(), patt % ('lh', atlas)))
        rgii = nib.load(os.path.join(library_path(), patt % ('rh', atlas)))
        cifti_split._atlas_cache[atlas] = tuple([pimms.imm_array(gii.darrays[0].data.astype('bool'))
                                                 for gii in (lgii, rgii)])
    (lroi,rroi) = cifti_split._atlas_cache[atlas]
    (ln,lN) = (np.sum(lroi), len(lroi))
    (rn,rN) = (np.sum(rroi), len(rroi))
    (ldat,rdat,sdat) = [np.full(dat.shape[:-1] + (k,), null) for k in [lN, rN, n - ln - rn]]
    ldat[..., lroi] = dat[..., :ln]
    rdat[..., rroi] = dat[..., ln:(ln+rn)]
    sdat[...] = dat[..., (ln+rn):]
    if ln + rn >= n: sdat = None
    return (ldat, rdat, sdat) 

def cifti_extract(cii, h, null=np.nan):
    '''
    cifti_extract(cii, h) yields the portion of the cifti vector or matrix that is associated
      with the given hemisphere h of the given subject. If h is None, then yields any subcortical
      voxels (or None if none are included).
    '''
    (l,r,s) = cifti_split(cii, null=null)
    if h is None:            return s
    elif h.startswith('lh'): return l
    else:                    return r

####################################################################################################
# The filemap instructions and related helper data/functions 

def line_intersection_2D(abarg, cdarg, atol=1e-8):
    '''
    line_intersection((a, b), (c, d)) yields the intersection point between the lines that pass
    through the given pairs of points. If any lines are parallel, (numpy.nan, numpy.nan) is
    returned; note that a, b, c, and d can all be 2 x n matrices of x and y coordinate row-vectors.
    '''
    ((x1,y1),(x2,y2)) = abarg
    ((x3,y3),(x4,y4)) = cdarg
    dx12 = (x1 - x2)
    dx34 = (x3 - x4)
    dy12 = (y1 - y2)
    dy34 = (y3 - y4)
    denom = dx12*dy34 - dy12*dx34
    unit = np.isclose(denom, 0, atol=atol)
    if unit is True: return (np.nan, np.nan)
    denom = unit + denom
    q12 = (x1*y2 - y1*x2) / denom
    q34 = (x3*y4 - y3*x4) / denom
    xi = q12*dx34 - q34*dx12
    yi = q12*dy34 - q34*dy12
    if   unit is False: return (xi, yi)
    elif unit is True:  return (np.nan, np.nan)
    else:
        xi = np.asarray(xi)
        yi = np.asarray(yi)
        xi[unit] = np.nan
        yi[unit] = np.nan
        return (xi, yi) 

def segment_intersection_2D(p12arg, p34arg, atol=1e-8):
    '''
    segment_intersection((a, b), (c, d)) yields the intersection point between the line segments
    that pass from point a to point b and from point c to point d. If there is no intersection
    point, then (numpy.nan, numpy.nan) is returned.
    '''
    (p1,p2) = p12arg
    (p3,p4) = p34arg
    pi = np.asarray(line_intersection_2D(p12arg, p34arg, atol=atol))
    p1 = np.asarray(p1)
    p2 = np.asarray(p2)
    p3 = np.asarray(p3)
    p4 = np.asarray(p4)
    u12 = p2 - p1
    u34 = p4 - p3
    cfn = lambda px,iis: (px if iis is None or len(px.shape) == 1 or px.shape[1] == len(iis) else
                          px[:,iis])
    dfn = lambda a,b:     a[0]*b[0] + a[1]*b[1]
    sfn = lambda a,b:     ((a-b)                 if len(a.shape) == len(b.shape) else
                           (np.transpose([a])-b) if len(a.shape) <  len(b.shape) else
                           (a - np.transpose([b])))
    fn  = lambda px,iis:  (1 - ((dfn(cfn(u12,iis), sfn(         px, cfn(p1,iis))) > 0) *
                                (dfn(cfn(u34,iis), sfn(         px, cfn(p3,iis))) > 0) *
                                (dfn(cfn(u12,iis), sfn(cfn(p2,iis),          px)) > 0) *
                                (dfn(cfn(u34,iis), sfn(cfn(p4,iis),          px)) > 0)))
    if len(pi.shape) == 1:
        if not np.isfinite(pi[0]): return (np.nan, np.nan)
        bad = fn(pi, None)
        return (np.nan, np.nan) if bad else pi
    else:
        nonpar = np.where(np.isfinite(pi[0]))[0]
        bad = fn(cfn(pi, nonpar), nonpar)
        (xi,yi) = pi
        bad = nonpar[np.where(bad)[0]]
        xi[bad] = np.nan
        yi[bad] = np.nan
        return (xi,yi)