Python numpy.piecewise() Examples

def quintic_spline(dx, dtype=np.float64):
    def inner(x):
        return 1 + x ** 3 / 12 * (-95 + 138 * x - 55 * x ** 2)

    def middle(x):
        return (x - 1) * (x - 2) / 24 * (-138 + 348 * x - 249 * x ** 2 + 55 * x ** 3)

    def outer(x):
        return (x - 2) * (x - 3) ** 2 / 24 * (-54 + 50 * x - 11 * x ** 2)

    window = np.arange(-3, 4)
    x = np.abs(dx - window)
    result = np.piecewise(
        x,
        [x <= 1, (x > 1) & (x <= 2), (x > 2) & (x <= 3)],
        [lambda x: inner(x), lambda x: middle(x), lambda x: outer(x)],
    )
    return result, window 

def _SetCPUCount(self):
    """Set cpu count.

    GCP: ratings from
    https://cloud.google.com/compute/docs/disks/performance#ssd-pd-performance
    AWS: to achieve good performance on EBS, one needs to use an
    EBS-optimized VM instance, and the smallest VM instance that can be EBS
    optimized is *.large VM types (e.g., c4.large), those comes with 2 cores.
    ratings from
    http://docs.aws.amazon.com/AWSEC2/latest/UserGuide/EBSOptimized.html

    """
    if self._provider == GCP:
      value = self._iops
      self._cpu_count = int(
          numpy.piecewise([value], [[value <= 15000], [
              (value > 15000) & (value <= 25000)
          ], [value > 25000]], [lambda x: 1, lambda x: 16, lambda x: 32]))
    elif self._provider == AWS:
      self._cpu_count = 2 

def _SetSize(self):
    """Set minimum size (GB) necessary to achieve _iops level.

    Rating performance levels as of May 2017, sources found below.
    GCP: ratings from https://cloud.google.com/compute/docs/disks/. Storage can
          go as high as 64TB per disk but IOPS maxes out at 30,000 iops/disk
          or (1000GB).
    AWS: ratings from
          http://docs.aws.amazon.com/AWSEC2/latest/
          UserGuide/EBSVolumeTypes.html#EBSVolumeTypes_gp2. Storage can go as
          high as 16TiB per disk but IOPS maxes out at 10000 IOPS/volume size.
          Reference gives volume size in GiB so converted to GB below.
    """
    if self._provider == GCP:
      self._size = min(int(math.ceil(self._iops / 30.0)), 30000 / 30)
    elif self._provider == AWS:
      value = self._iops
      value = numpy.array(value)
      self._size = int(
          numpy.piecewise([value], [[value <= 100], [(value > 100) & (
              value <= 9999)], [value > 9999]], [
                  lambda x: int(math.ceil(1.07374)),
                  lambda x: int(math.ceil(3 * value)),
                  lambda x: int(math.ceil(3579.855))])) 

def get_frame_negative(self,frame_data):

		# print(frame_data.size)
		frame = np.rec.array(None, dtype=[('value', np.float16),('valid', np.bool_)], shape=(self.height, self.width))
		frame.valid.fill(False)
		frame.value.fill(0.)
		# print(frame.size)

		for datum in np.nditer(frame_data):
			# print(datum['y'])
			ts_val = datum['ts']
			f_data = frame[datum['y'], datum['x']]
			if(datum['p']== False):
				f_data.value += 1

		img = frame.value/20*255
		img = img.astype('uint8')
		# img = np.piecewise(img, [img <= 0, (img > 0) & (img < 255), img >= 255], [0, lambda x: x, 255])
		cv2.normalize(img,img,0,255,cv2.NORM_MINMAX)
		img = cv2.flip(img, 1)
		img = np.rot90(img)
		# cv2.imshow('img_neg', img)
		# cv2.waitKey(0)
		return img 

def get_frame(self,frame_data):

		# print(frame_data.size)
		frame = np.rec.array(None, dtype=[('value', np.float16),('valid', np.bool_)], shape=(self.height, self.width))
		frame.valid.fill(False)
		frame.value.fill(0.)
		# print(frame.size)

		for datum in np.nditer(frame_data, flags=['zerosize_ok']):
			# print(datum['y'])
			ts_val = datum['ts']
			f_data = frame[datum['y'], datum['x']]
			f_data.value += 1

		img = frame.value/20*255
		img = img.astype('uint8')
		# img = np.piecewise(img, [img <= 0, (img > 0) & (img < 255), img >= 255], [0, lambda x: x, 255])
		# cv2.normalize(img,img,0,255,cv2.NORM_L1)
		cv2.normalize(img,img,0,255,cv2.NORM_MINMAX)
		img = cv2.flip(img, 1)
		img = np.rot90(img)
		# cv2.imshow('img_f', img)
		# cv2.waitKey(0)
		return img 

def get_random(self, size=None):
        """Draw a random number from the distribution.

        If size is not None but an integer N, return an array of N numbers.

        For the MultivariateNumericalDistribution, the PDF from which the
        random numbers are drawn is approximated to be piecewise constant in
        hypercubes around the points of the lattice spanned by the `xi`. A finer
        lattice spacing will lead to a smoother distribution of random numbers
        (but will also be slower).
        """

        if size is None:
            return self._get_random()
        else:
            return np.array([self._get_random() for i in range(size)]) 

def minimum(self, ndim):
        return -2.903534 * np.ones(ndim)


# class Simionescu(OptimizationTestFunction):
#
#     def __init__(self):
#         OptimizationTestFunction.__init__(self, mindim=2, maxdim=2, domain=np.array([-1.25, 1.25]))
#
#     @staticmethod
#     def function(x):
#         rt = 1
#         rs = 0.2
#         n = 8
#         return np.piecewise(x,
#                      [x[0]**2 + x[1]**2 <= (rt + rs*np.cos(n*np.arctan(x[0]/x[1])))**2,
#                       x[0]**2 + x[1]**2 > (rt + rs*np.cos(n*np.arctan(x[0]/x[1])))**2], [0.1*x[0]*x[1], 1])
#
#
#     def minimum(self, ndim):
#         assert ndim == 2
#         return -0.84852813*np.ones(ndim) 

def B_0123(x, der=0):
    """A quadratic B-spline function B(x | 0, 1, 2, 3)."""
    x = np.atleast_1d(x)
    conds = [x < 1, (x > 1) & (x < 2), x > 2]
    if der == 0:
        funcs = [lambda x: x*x/2.,
                 lambda x: 3./4 - (x-3./2)**2,
                 lambda x: (3.-x)**2 / 2]
    elif der == 2:
        funcs = [lambda x: 1.,
                 lambda x: -2.,
                 lambda x: 1.]
    else:
        raise ValueError('never be here: der=%s' % der)
    pieces = np.piecewise(x, conds, funcs)
    return pieces 

def spherical_covariance_model(self, d):
        '''Spherical covariance model, effective range equals range parameter, valid in R3'''
        psill = self.sill - self.nugget
        d = d.astype(float)
        gamma = np.piecewise(d, [d <= self.range_, d > self.range_],
                             [lambda d:
                              psill * (1 - ((3. * d) / (2. * self.range_)
                                            - (d ** 3.) / (2. * self.range_ ** 3.))),
                              lambda d: 0])
        return gamma

    # TODO: Make this better and nicer and everything
    # option for covariance
    # display range, sill, nugget, practical range etc. 

def test_piecewise(self):
        x = np.linspace(-2.5, 2.5, 6) * u.m
        out = np.piecewise(x, [x < 0, x >= 0], [-1*u.s, 1*u.day])
        expected = np.piecewise(x.value, [x.value < 0, x.value >= 0],
                                [-1, 24*3600]) * u.s
        assert out.unit == expected.unit
        assert np.all(out == expected)

        out2 = np.piecewise(x, [x < 1 * u.m, x >= 0],
                            [-1*u.s, 1*u.day, lambda x: 1*u.hour])
        expected2 = np.piecewise(x.value, [x.value < 1, x.value >= 0],
                                 [-1, 24*3600, 3600]) * u.s
        assert out2.unit == expected2.unit
        assert np.all(out2 == expected2)

        out3 = np.piecewise(x, [x < 1 * u.m, x >= 0],
                            [0, 1*u.percent, lambda x: 1*u.one])
        expected3 = np.piecewise(x.value, [x.value < 1, x.value >= 0],
                                 [0, 0.01, 1]) * u.one
        assert out3.unit == expected3.unit
        assert np.all(out3 == expected3)

        with pytest.raises(TypeError):  # no Quantity in condlist.
            np.piecewise(x, [x], [0.])

        with pytest.raises(TypeError):  # no Quantity in condlist.
            np.piecewise(x.value, [x], [0.]) 

def B_012(x):
    """ A linear B-spline function B(x | 0, 1, 2)."""
    x = np.atleast_1d(x)
    return np.piecewise(x, [(x < 0) | (x > 2),
                            (x >= 0) & (x < 1),
                            (x >= 1) & (x <= 2)],
                           [lambda x: 0., lambda x: x, lambda x: 2.-x]) 

def stabilize(td):
    """Compensate for motion of the ATIS sensor during recording of the Neuromorphic datasets
    Applies to the N-MNIST and N-Caltech101 datasets.
    The image motion is originally induced by egorotation of the ATIS sensor
    td: eventvision.Events
    """
    assert isinstance(td, ev.Events)

    def correct_saccade1(data):
        data.x -= np.rint(3.5 * data.ts / 105e3).astype(np.uint16)
        data.y -= np.rint(7 * data.ts / 105e3).astype(np.uint16)
        return data

    def correct_saccade2(data):
        data.x -= np.rint(3.5 + 3.5 * (data.ts - 105e3) / 105e3).astype(np.uint16)
        data.y -= np.rint(7 - 7 * (data.ts - 105e3) / 105e3).astype(np.uint16)
        return data

    def correct_saccade3(data):
        data.x -= np.rint(7 - 7 * (data.ts - 210e3) / 105e3).astype(np.uint16)
        return data

    copy = np.piecewise(td.data,\
        [td.data.ts <= 105e3, (td.data.ts > 105e3) & (td.data.ts <= 210e3), (td.data.ts > 210e3)],\
        [correct_saccade1, correct_saccade2, correct_saccade3]).view(np.recarray)

    # after saccades, we might end up with invalid x and y values, have to
    # correct these
    x_vals = copy.x
    y_vals = copy.y
    copy.x = np.piecewise(x_vals,\
        [x_vals >= 65000, (x_vals < 65000) & (x_vals >= td.width), x_vals < td.width],\
        [0, td.width - 1, lambda x: x])
    copy.y = np.piecewise(y_vals,\
        [y_vals >= 65000, (y_vals < 65000) & (y_vals >= td.height), y_vals < td.height],\
        [0, td.height - 1, lambda y: y])

    return copy 

def show_td(self, wait_delay=1):
        """Displays the TD events (change detection ATIS or DVS events)
        waitDelay: milliseconds
        """
        frame_length = 24e3
        t_max = self.data.ts[-1]
        frame_start = self.data[0].ts
        frame_end = self.data[0].ts + frame_length
        td_img = np.ones((self.height, self.width), dtype=np.uint8)
        while frame_start < t_max:
            frame_data = self.data[(self.data.ts >= frame_start) & (self.data.ts < frame_end)]
            
            if frame_data.size > 0:
                td_img.fill(128)

                #with timer.Timer() as em_playback_timer:
                for datum in np.nditer(frame_data):
                    td_img[datum['y'].item(0), datum['x'].item(0)] = datum['p'].item(0)
                #print 'prepare td frame by iterating events took %s seconds'
                #%em_playback_timer.secs

                td_img = np.piecewise(td_img, [td_img == 0, td_img == 1, td_img == 128], [0, 255, 128])
                cv2.imshow('img', td_img)
                cv2.waitKey(wait_delay)

            frame_start = frame_end + 1
            frame_end = frame_end + frame_length + 1

        cv2.destroyAllWindows()
        return 

def gaussian_square(times: np.ndarray, amp: complex, center: float, square_width: float,
                    sigma: float, zeroed_width: Optional[float] = None) -> np.ndarray:
    r"""Continuous gaussian square pulse.

    Args:
        times: Times to output pulse for.
        amp: Pulse amplitude.
        center: Center of the square pulse component.
        square_width: Width of the square pulse component.
        sigma: Standard deviation of Gaussian rise/fall portion of the pulse.
        zeroed_width: Subtract baseline of gaussian square pulse
            to enforce $\OmegaSquare(center \pm zeroed_width/2)=0$.

    Raises:
        PulseError: if zeroed_width is not compatible with square_width.
    """
    square_start = center-square_width/2
    square_stop = center+square_width/2
    if zeroed_width:
        if zeroed_width < square_width:
            raise PulseError("zeroed_width cannot be smaller than square_width.")
        gaussian_zeroed_width = zeroed_width-square_width
    else:
        gaussian_zeroed_width = None

    funclist = [functools.partial(gaussian, amp=amp, center=square_start, sigma=sigma,
                                  zeroed_width=gaussian_zeroed_width, rescale_amp=True),
                functools.partial(gaussian, amp=amp, center=square_stop, sigma=sigma,
                                  zeroed_width=gaussian_zeroed_width, rescale_amp=True),
                functools.partial(constant, amp=amp)]
    condlist = [times <= square_start, times >= square_stop]
    return np.piecewise(times.astype(np.complex_), condlist, funclist) 

def __call__(self, a):
        self.variables = (a,)
        return np.piecewise(
            a.data, [a.data == 0, a.data != 0], [np.pi / 2, lambda x: np.arctan(1 / x)]
        ) 

def spherical_variogram_model(self, d):
        '''Spherical variogram model, effective range equals range parameter, valid in R3'''
        psill = self.sill - self.nugget
        d = d.astype(float)
        gamma = np.piecewise(d, [d <= self.range_, d > self.range_],
                             [lambda d:
                              psill * ((3. * d) / (2. * self.range_)
                                       - (d ** 3.) / (2. * self.range_ ** 3.)) + self.nugget,
                              lambda d: self.sill])
        return gamma 

def backward_var(self, grad, index, **kwargs):
        a = self.variables[index]
        return grad * np.piecewise(
            a.data, [a.data < 0, a.data == 0, a.data > 0], [-1, np.nan, 1]
        ) 

def logpdf(self, x):
        # return -inf for negative x values
        inf0 = np.piecewise(np.asarray(x, dtype=float), [x<0, x>=0], [-np.inf, 0.])
        return inf0 + self.scipy_dist.logpdf(x) + np.log(self._pdf_scale) 

def cdf(self, x):
        cdf0 = self.scipy_dist.cdf(0)
        cdf = (self.scipy_dist.cdf(x) - cdf0)/(1-cdf0)
        return np.piecewise(
                    np.asarray(x, dtype=float),
                    [x<0, x>=0],
                    [0., cdf]) # return 0 for negative x 

def get_generator(self):
        return lambda x: numpy.piecewise(x, self.domain(x), [self._generators[j] for j in range(len(self._generators))] + [self.finals[-1]]) 

def backward_var(self, grad, index, **kwargs):
        a = self.variables[index]
        x = a.data
        return np.pi * grad * np.piecewise(x, [x == 0, x != 0], [np.zeros_like, _dsinc]) 

def aflare(t, p):
    """
    This is the Analytic Flare Model from the flare-morphology paper.
    Reference Davenport et al. (2014) http://arxiv.org/abs/1411.3723
    Note: this model assumes the flux before the flare is zero centered
    Note: many sub-flares can be modeled by this method by changing the
    number of parameters in "p". As a result, this routine may not work
    for fitting with methods like scipy.optimize.curve_fit, which require
    a fixed number of free parameters. Instead, for fitting a single peak
    use the aflare1 method.
    Parameters
    ----------
    t : 1-d array
        The time array to evaluate the flare over
    p : 1-d array
        p == [tpeak, fwhm (units of time), amplitude (units of flux)] x N
    Returns
    -------
    flare : 1-d array
        The flux of the flare model evaluated at each time
    """
    _fr = [1.00000, 1.94053, -0.175084, -2.24588, -1.12498]
    _fd = [0.689008, -1.60053, 0.302963, -0.278318]

    Nflare = int( np.floor( (len(p)/3.0) ) )

    flare = np.zeros_like(t)
    # compute the flare model for each flare
    for i in range(Nflare):
        outm = np.piecewise(t, [(t<= p[0+i*3]) * (t-p[0+i*3])/p[1+i*3] > -1.,
                                (t > p[0+i*3])],
                            [lambda x: (_fr[0]+                             # 0th order
                                        _fr[1]*((x-p[0+i*3])/p[1+i*3])+     # 1st order
                                        _fr[2]*((x-p[0+i*3])/p[1+i*3])**2.+  # 2nd order
                                        _fr[3]*((x-p[0+i*3])/p[1+i*3])**3.+  # 3rd order
                                        _fr[4]*((x-p[0+i*3])/p[1+i*3])**4. ),# 4th order
                             lambda x: (_fd[0]*np.exp( ((x-p[0+i*3])/p[1+i*3])*_fd[1] ) +
                                        _fd[2]*np.exp( ((x-p[0+i*3])/p[1+i*3])*_fd[3] ))]
                            ) * p[2+i*3] # amplitude
        flare = flare + outm
        
    return flare 

def show_em(self):
        """Displays the EM events (grayscale ATIS events)"""
        frame_length = 24e3
        t_max = self.data.ts[-1]
        frame_start = self.data[0].ts
        frame_end = self.data[0].ts + frame_length
        max_val = 1.16e5
        min_val = 1.74e3
        val_range = max_val - min_val

        thr = np.rec.array(None, dtype=[('valid', np.bool_), ('low', np.uint64), ('high', np.uint64)], shape=(self.height, self.width))
        thr.valid.fill(False)
        thr.low.fill(frame_start)
        thr.high.fill(0)

        def show_em_frame(frame_data):
            """Prepare and show a single frame of em data to be shown"""
            for datum in np.nditer(frame_data):
                ts_val = datum['ts'].item(0)
                thr_data = thr[datum['y'].item(0), datum['x'].item(0)]

                if datum['p'].item(0) == 0:
                    thr_data.valid = 1
                    thr_data.low = ts_val
                elif thr_data.valid == 1:
                    thr_data.valid = 0
                    thr_data.high = ts_val - thr_data.low

            img = 255 * (1 - (thr.high - min_val) / (val_range))
            #thr_h = cv2.adaptiveThreshold(thr_h, 255,
            #cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, 3, 0)
            img = np.piecewise(img, [img <= 0, (img > 0) & (img < 255), img >= 255], [0, lambda x: x, 255])
            img = img.astype('uint8')
            cv2.imshow('img', img)
            cv2.waitKey(1)

        while frame_start < t_max:
            #with timer.Timer() as em_playback_timer:
            frame_data = self.data[(self.data.ts >= frame_start) & (self.data.ts < frame_end)]
            show_em_frame(frame_data)
            frame_start = frame_end + 1
            frame_end += frame_length + 1
            #print 'showing em frame took %s seconds' %em_playback_timer.secs

        cv2.destroyAllWindows()
        return 

def cubic_spline(dx, a=1, b=0):
    """Generate a cubix spline centered on `dx`.

    Parameters
    ----------
    dx: float
        Fractional amount that the kernel will be shifted
    a: float
        Cubic spline sharpness paremeter
    b: float
        Cubic spline shape parameter

    Returns
    -------
    result: array
        Cubic Spline kernel in a window from floor(dx)-1 to floor(dx) + 3
    window: array
        The pixel values for the window containing the kernel
    """
    if np.abs(dx) > 1:
        raise ValueError("The fractional shift dx must be between -1 and 1")

    def inner(x):
        """Cubic from 0<=abs(x)<=1
        """
        third = (-6 * a - 9 * b + 12) * x ** 3
        second = (6 * a + 12 * b - 18) * x ** 2
        zero = -2 * b + 6
        return (zero + second + third) / 6

    def outer(x):
        """Cubic from 1<=abs(x)<=2
        """
        third = (-6 * a - b) * x ** 3
        second = (30 * a + 6 * b) * x ** 2
        first = (-48 * a - 12 * b) * x
        zero = 24 * a + 8 * b
        return (zero + first + second + third) / 6

    window = np.arange(-1, 3) + np.floor(dx)
    x = np.abs(dx - window)
    result = np.piecewise(
        x, [x <= 1, (x > 1) & (x < 2)], [lambda x: inner(x), lambda x: outer(x)]
    )

    return result, np.array(window).astype(int)