Python numpy.log10() Examples

def test_mimo(self):
      # MIMO
      B = np.matrix('1,0;0,1')
      D = np.matrix('0,0')
      sysMIMO = ss(self.A,B,self.C,D)

      frqMIMO = sysMIMO.freqresp(self.omega)
      tfMIMO = tf(sysMIMO)

      #bode(sysMIMO) # - should throw not implemented exception
      #bode(tfMIMO) # - should throw not implemented exception

      #plt.figure(3)
      #plt.semilogx(self.omega,20*np.log10(np.squeeze(frq[0])))

      #plt.figure(4)
      #bode(sysMIMO,self.omega) 

def maxdmag(self, s):
        """Calculates the maximum value of dMag for projected separation
        
        Args:
            s (float):
                Projected separation (AU)
        
        Returns:
            maxdmag (float):
                Maximum planet delta magnitude
        
        """
        
        if s == 0.0:
            maxdmag = self.cdmax - 2.5*np.log10(self.Phi(np.pi))
        elif s < self.rmax:
            maxdmag = self.cdmax - 2.5*np.log10(np.abs(self.Phi(np.pi-np.arcsin(s/self.rmax))))
        else:
            maxdmag = self.cdmax - 2.5*np.log10(self.Phi(np.pi/2.0))

        return maxdmag 

def calclogf(self):
        """
        # wavelength dependence, from Table 19 in Leinert et al 1998
        # interpolated w/ a quadratic in log-log space
        Returns:
            interpolant (object):
                a 1D quadratic interpolant of intensity vs wavelength

        """
        self.zodi_lam = np.array([0.2, 0.3, 0.4, 0.5, 0.7, 0.9, 1.0, 1.2, 2.2, 3.5,
                4.8, 12, 25, 60, 100, 140]) # um
        self.zodi_Blam = np.array([2.5e-8, 5.3e-7, 2.2e-6, 2.6e-6, 2.0e-6, 1.3e-6,
                1.2e-6, 8.1e-7, 1.7e-7, 5.2e-8, 1.2e-7, 7.5e-7, 3.2e-7, 1.8e-8,
                3.2e-9, 6.9e-10]) # W/m2/sr/um
        x = np.log10(self.zodi_lam)
        y = np.log10(self.zodi_Blam)
        return interp1d(x, y, kind='quadratic') 

def test_dist_sma_radius(self):
        """
        Test that sma and radius values outside of the range have zero probability
        """
        
        for mod in self.allmods:
            if 'dist_sma_radius' in mod.__dict__:
                with RedirectStreams(stdout=self.dev_null):
                    pp = mod(**self.spec)
                
                a = np.logspace(np.log10(pp.arange[0].value/10.),np.log10(pp.arange[1].value*100),100)
                Rp = np.logspace(np.log10(pp.Rprange[0].value/10.),np.log10(pp.Rprange[1].value*100),100)
                
                aa, RR = np.meshgrid(a,Rp)
                
                fr = pp.dist_sma_radius(aa,RR)
                self.assertTrue(np.all(fr[aa < pp.arange[0].value] == 0),'dist_sma_radius low bound failed on sma for %s'%mod.__name__)
                self.assertTrue(np.all(fr[aa > pp.arange[1].value] == 0),'dist_sma_radius high bound failed on sma for %s'%mod.__name__)
                self.assertTrue(np.all(fr[RR < pp.Rprange[0].value] == 0),'dist_sma_radius low bound failed on radius for %s'%mod.__name__)
                self.assertTrue(np.all(fr[RR > pp.Rprange[1].value] == 0),'dist_sma_radius high bound failed on radius for %s'%mod.__name__)
                self.assertTrue(np.all(fr[(aa > pp.arange[0].value) & (aa < pp.arange[1].value) & (RR > pp.Rprange[0].value) & (RR < pp.Rprange[1].value)] > 0),'dist_sma_radius is improper pdf for %s'%mod.__name__) 

def test_dist_mass(self):
        """
        Test that masses outside of the range have zero probability

        """

        for mod in self.allmods:
            if 'dist_mass' in mod.__dict__:
                with RedirectStreams(stdout=self.dev_null):
                    pp = mod(**self.spec)

                Mp = np.logspace(np.log10(pp.Mprange[0].value/10.),np.log10(pp.Mprange[1].value*100.),100) 

                fr = pp.dist_mass(Mp)
                self.assertTrue(np.all(fr[Mp < pp.Mprange[0].value] == 0),'dist_mass high bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fr[Mp > pp.Mprange[1].value] == 0),'dist_mass high bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fr[(Mp >= pp.Mprange[0].value) & (Mp <= pp.Mprange[1].value)] > 0)) 

def test_dist_radius(self):
        """
        Test that radii outside of the range have zero probability

        """
        for mod in self.allmods:
            if 'dist_radius' in mod.__dict__:
                with RedirectStreams(stdout=self.dev_null):
                    pp = mod(**self.spec)

                Rp = np.logspace(np.log10(pp.Rprange[0].to('earthRad').value/10.),np.log10(pp.Rprange[1].to('earthRad').value*100.),100) 

                fr = pp.dist_radius(Rp)
                self.assertTrue(np.all(fr[Rp < pp.Rprange[0].to('earthRad').value] == 0),'dist_radius high bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fr[Rp > pp.Rprange[1].to('earthRad').value] == 0),'dist_radius high bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fr[(Rp >= pp.Rprange[0].to('earthRad').value) & (Rp <= pp.Rprange[1].to('earthRad').value)] > 0),'dist_radius generates zero probabilities within range for %s'%mod.__name__) 

def test_dist_sma(self):
        """
        Test that smas outside of the range have zero probability

        """

        for mod in self.allmods:
            if 'dist_sma' in mod.__dict__:
                with RedirectStreams(stdout=self.dev_null):
                    pp = mod(**self.spec)

                a = np.logspace(np.log10(pp.arange[0].to('AU').value/10.),np.log10(pp.arange[1].to('AU').value*10.),100)

                fa = pp.dist_sma(a)
                self.assertTrue(np.all(fa[a < pp.arange[0].to('AU').value] == 0),'dist_sma high bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fa[a > pp.arange[1].to('AU').value] == 0),'dist_sma low bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fa[(a >= pp.arange[0].to('AU').value) & (a <= pp.arange[1].to('AU').value)] >= 0.),'dist_sma generates negative densities within range for %s'%mod.__name__) 

def zdr(z_h, z_v):
    """
    Differential reflectivity [dB].

    From Rinehart (1997), Eqn 10.3 and Seliga and Bringi (1976)

    Parameters
    ----------
    z_h : float or array
        Horizontal reflectivity [mm^6/m^3]
    z_v : float or array
        Vertical reflectivity [mm^6/m^3]

    Notes
    -----
    Ensure that both powers have the same units!

    Alternating horizontally and linearly polarized pulses are averaged.

    Notes
    -----
    If arrays are used, either for one okay, or both must be the same length.
    """
    return 10. * np.log10(np.asarray(z_h) / np.asarray(z_v)) 

def gain_Pratio(p1, p2):
    """
    Antenna gain [dB] via power ratio.

    From Rinehart (1997), Eqn 2.1

    Parameters
    ----------
    p1 : float or array
        Power on the beam axis [W]
    p2 : float or array
        Power from an isotropic antenna [W]

    Notes
    -----
    Ensure that both powers have the same units!
    If arrays are used, either for one okay, or both must be the same length.
    """
    return 10. * np.log10(np.asarray(p1) / np.asarray(p2)) 

def deltaMag(p,Rp,d,Phi):
    """ Calculates delta magnitudes for a set of planets, based on their albedo, 
    radius, and position with respect to host star.
    
    Args:
        p (ndarray):
            Planet albedo
        Rp (astropy Quantity array):
            Planet radius in units of km
        d (astropy Quantity array):
            Planet-star distance in units of AU
        Phi (ndarray):
            Planet phase function
    
    Returns:
        dMag (ndarray):
            Planet delta magnitudes
    
    """
    dMag = -2.5*np.log10(p*(Rp/d).decompose()**2*Phi).value
    
    return dMag 

def test_LGST_1overSqrtN_dependence(self):
        my_datagen_gateset = self.model.depolarize(op_noise=0.05, spam_noise=0)
        # !!don't depolarize spam or 1/sqrt(N) dependence saturates!!

        nSamplesList = np.array([ 16, 128, 1024, 8192 ])
        diffs = []
        for nSamples in nSamplesList:
            ds = pygsti.construction.generate_fake_data(my_datagen_gateset, self.lgstStrings, nSamples,
                                                        sampleError='binomial', seed=100)
            mdl_lgst = pygsti.do_lgst(ds, self.fiducials, self.fiducials, self.model, svdTruncateTo=4, verbosity=0)
            mdl_lgst_go = pygsti.gaugeopt_to_target(mdl_lgst, my_datagen_gateset, {'spam':1.0, 'gate': 1.0}, checkJac=True)
            diffs.append( my_datagen_gateset.frobeniusdist(mdl_lgst_go) )

        diffs = np.array(diffs, 'd')
        a, b = polyfit(np.log10(nSamplesList), np.log10(diffs), deg=1)
        #print "\n",nSamplesList; print diffs; print a #DEBUG
        self.assertLess( a+0.5, 0.05 ) 

def speech_aug_scheduler(step, s_r, s_i, s_f, peak_lr):
    # Starting from 0, ramp-up to set LR and  converge to 0.01*LR, w/ exp. decay
    final_lr_ratio = 0.01
    exp_decay_lambda = -np.log10(final_lr_ratio)/(s_f-s_i) # Approx. w/ 10-based
    cur_step = step+1

    if cur_step<s_r:
        # Ramp-up
        return peak_lr*float(cur_step)/s_r
    elif cur_step<s_i:
        # Hold
        return peak_lr
    elif cur_step<=s_f:
        # Decay
        return peak_lr*np.power(10,-exp_decay_lambda*(cur_step-s_i))
    else:
        # Converge
        return peak_lr*final_lr_ratio 

def logSE(x1, x2=-1):
    """
    10 * log10(e**2)    
    This function accepts two series of data or directly
    one series with error.

    **Args:**

    * `x1` - first data series or error (1d array)

    **Kwargs:**

    * `x2` - second series (1d array) if first series was not error directly,\\
        then this should be the second series

    **Returns:**

    * `e` - logSE of error (1d array) obtained directly from `x1`, \\
        or as a difference of `x1` and `x2`. The values are in dB!

    """
    e = get_valid_error(x1, x2)
    return 10*np.log10(e**2) 

def db2mag(db):
    """Convert a gain in decibels (dB) to a magnitude

    If A is magnitude,

        db = 20 * log10(A)

    Parameters
    ----------
    db : float or ndarray
        input value or array of values, given in decibels

    Returns
    -------
    mag : float or ndarray
        corresponding magnitudes

    """
    return 10. ** (db / 20.) 

def mag2db(mag):
    """Convert a magnitude to decibels (dB)

    If A is magnitude,

        db = 20 * log10(A)

    Parameters
    ----------
    mag : float or ndarray
        input magnitude or array of magnitudes

    Returns
    -------
    db : float or ndarray
        corresponding values in decibels

    """
    return 20. * np.log10(mag) 

def plotstft(audiopath, binsize=2**10, plotpath=None, colormap="gray", channel=0, name='tmp.png', alpha=1, offset=0):
    samplerate, samples = wav.read(audiopath)
    samples = samples[:, channel]
    s = stft(samples, binsize)

    sshow, freq = logscale_spec(s, factor=1, sr=samplerate, alpha=alpha)
    sshow = sshow[2:, :]
    ims = 20.*np.log10(np.abs(sshow)/10e-6) # amplitude to decibel
    timebins, freqbins = np.shape(ims)
    
    ims = np.transpose(ims)
    # ims = ims[0:256, offset:offset+768] # 0-11khz, ~9s interval
    ims = ims[0:256, :] # 0-11khz, ~10s interval
    #print "ims.shape", ims.shape
    
    image = Image.fromarray(ims) 
    image = image.convert('L')
    image.save(name) 

def _interp_ebv(self, datum, dist):
        """
        Calculate samples of E(B-V) at an arbitrary distance (in kpc) for one
        test coordinate.
        """
        dm = 5. * (np.log10(dist) + 2.)
        idx_ceil = np.searchsorted(datum['DM_bin_edges'], dm)
        if idx_ceil == 0:
            dist_0 = 10.**(datum['DM_bin_edges'][0]/5. - 2.)
            return dist/dist_0 * datum['samples'][:,0]
        elif idx_ceil == len(datum['DM_bin_edges']):
            return datum['samples'][:,-1]
        else:
            dm_ceil = datum['DM_bin_edges'][idx_ceil]
            dm_floor = datum['DM_bin_edges'][idx_ceil-1]
            a = (dm_ceil - dm) / (dm_ceil - dm_floor)
            return (
                (1.-a) * datum['samples'][:,idx_ceil]
                +    a * datum['samples'][:,idx_ceil-1]
            ) 

def get_audio_data(self):
        frames = self.rec.get_frames()
        result = [0] * self.bins
        if len(frames) > 0:
            # keeps only the last frame
            current_frame = frames[-1]
            # plots the time signal
            # self.line_top.set_data(self.time_vect, current_frame)
            # computes and plots the fft signal
            fft_frame = np.fft.rfft(current_frame)
            if self.auto_gain:
                fft_frame /= np.abs(fft_frame).max()
            else:
                fft_frame *= (1 + self.gain) / 5000000.

            fft_frame = np.abs(fft_frame)
            if self.log_scale:
                fft_frame = np.log10(np.add(1, np.multiply(10, fft_frame)))

            result = [min(int(max(i, 0.) * 1023), 1023) for i in fft_frame][0:self.bins]

        return result 

def get_audio_data(self):
        frames = self.rec.get_frames()
        result = [0] * self.bins
        if len(frames) > 0:
            # keeps only the last frame
            current_frame = frames[-1]
            # plots the time signal
            # self.line_top.set_data(self.time_vect, current_frame)
            # computes and plots the fft signal
            fft_frame = np.fft.rfft(current_frame)
            if self.auto_gain:
                fft_frame /= np.abs(fft_frame).max()
            else:
                fft_frame *= (1 + self.gain) / 5000000.

            fft_frame = np.abs(fft_frame)
            if self.log_scale:
                fft_frame = np.log10(np.add(1, np.multiply(10, fft_frame)))

            result = [min(int(max(i, 0.) * 1023), 1023) for i in fft_frame][0:self.bins]

        return result 

def test_amplitude_to_decibel():
    """test for backend_keras.amplitude_to_decibel"""
    from kapre.backend_keras import amplitude_to_decibel

    x = np.array([[1e-20, 1e-5, 1e-3, 5e-2], [0.3, 1.0, 20.5, 9999]])  # random positive numbers

    amin = 1e-5
    dynamic_range = 80.0

    x_decibel = 10 * np.log10(np.maximum(x, amin))
    x_decibel = x_decibel - np.max(x_decibel, axis=(1,), keepdims=True)
    x_decibel_ref = np.maximum(x_decibel, -1 * dynamic_range)

    x_var = K.variable(x)
    x_decibel_kapre = amplitude_to_decibel(x_var, amin, dynamic_range)

    assert np.allclose(K.eval(x_decibel_kapre), x_decibel_ref, atol=TOL) 

def show_plot(
    key, segment_times, sample_freqs, spec, duration, wav_data, vad_feat
):
    """This function plots the vad against the signal and the spectrogram.

    Args:
        segment_times: the time intervals acting as the x axis
        sample_freqs: the frequency bins acting as the y axis
        spec: the spectrogram
        duration: duration of the wave file
        wav_data: the wave data
        vad_feat: VAD features
    """

    import matplotlib.pyplot as plt
    import matplotlib.mlab as mlb

    plt.subplot(3, 1, 1)
    plt.pcolormesh(segment_times, sample_freqs, 10 * np.log10(spec), cmap="jet")
    plt.ylabel("Frequency [Hz]")
    plt.xlabel("Time [sec]")

    plt.subplot(3, 1, 2)
    axes = plt.gca()
    axes.set_xlim([0, duration])
    tmp_axis = np.linspace(0, duration, wav_data.shape[0])
    plt.plot(tmp_axis, wav_data / np.abs(np.max(wav_data)))
    plt.xlabel("Time [sec]")

    plt.subplot(3, 1, 3)
    axes = plt.gca()
    axes.set_xlim([0, duration])
    tmp_axis = np.linspace(0, duration, vad_feat.shape[0])
    plt.plot(tmp_axis, vad_feat)
    plt.xlabel("Time [sec]")

    plt.savefig("plots/" + key, bbox_inches="tight") 

def composed(cls, terms_to_compose, magnitude):
        logmag = _math.log10(magnitude) if magnitude > 0 else -LARGE
        first = terms_to_compose[0]
        coeffrep = first.coeff
        pre_ops = first.pre_ops[:]
        post_ops = first.post_ops[:]
        for t in terms_to_compose[1:]:
            coeffrep = coeffrep.mult(t.coeff)
            pre_ops += t.pre_ops
            post_ops += t.post_ops
        return SVTermRep(coeffrep, magnitude, logmag, first.pre_state, first.post_state,
                         first.pre_effect, first.post_effect, pre_ops, post_ops) 

def set_magnitude(self, mag):
        self.magnitude = mag
        self.logmagnitude = _math.log10(mag) if mag > 0 else -LARGE 

def db(audio):
    if len(audio.shape) > 1:
        maxx = np.max(np.abs(audio), axis=1)
        return 20 * np.log10(maxx) if np.any(maxx != 0) else np.array([0])
    maxx = np.max(np.abs(audio))
    return 20 * np.log10(maxx) if maxx != 0 else np.array([0]) 

def length_over_time(dfs, path, figformat, title, log_length=False, plot_settings={}):
    if log_length:
        time_length = Plot(path=path + "TimeLogLengthViolinPlot." + figformat,
                           title="Violin plot of log read lengths over time")
    else:
        time_length = Plot(path=path + "TimeLengthViolinPlot." + figformat,
                           title="Violin plot of read lengths over time")
    sns.set(style="white", **plot_settings)
    if log_length:
        length_column = "log_lengths"
    else:
        length_column = "lengths"

    if "length_filter" in dfs:  # produced by NanoPlot filtering of too long reads
        temp_dfs = dfs[dfs["length_filter"]]
    else:
        temp_dfs = dfs

    ax = sns.violinplot(x="timebin",
                        y=length_column,
                        data=temp_dfs,
                        inner=None,
                        cut=0,
                        linewidth=0)
    ax.set(xlabel='Interval (hours)',
           ylabel="Read length",
           title=title or time_length.title)
    if log_length:
        ticks = [10**i for i in range(10) if not 10**i > 10 * np.amax(dfs["lengths"])]
        ax.set(yticks=np.log10(ticks),
               yticklabels=ticks)
    plt.xticks(rotation=45, ha='center', fontsize=8)
    time_length.fig = ax.get_figure()
    time_length.save(format=figformat)
    plt.close("all")
    return time_length 

def _hrv_nonlinear_poincare(rri, out):
    """Compute SD1 and SD2
    - Do existing measures of Poincare plot geometry reflect nonlinear features of heart rate
    variability? - Brennan (2001)
    """

    # HRV and hrvanalysis
    rri_n = rri[:-1]
    rri_plus = rri[1:]
    x1 = (rri_n - rri_plus) / np.sqrt(2)  # Eq.7
    x2 = (rri_n + rri_plus) / np.sqrt(2)
    sd1 = np.std(x1, ddof=1)
    sd2 = np.std(x2, ddof=1)

    out["SD1"] = sd1
    out["SD2"] = sd2

    # SD1 / SD2
    out["SD1SD2"] = sd1 / sd2

    # Area of ellipse described by SD1 and SD2
    out["S"] = np.pi * out["SD1"] * out["SD2"]

    # CSI / CVI
    T = 4 * out["SD1"]
    L = 4 * out["SD2"]
    out["CSI"] = L / T
    out["CVI"] = np.log10(L * T)
    out["CSI_Modified"] = L ** 2 / T

    return out 

def _rsp_rrv_nonlinear(bbi):
    diff_bbi = np.diff(bbi)
    out = {}

    # Poincaré plot
    out["SD1"] = np.sqrt(np.std(diff_bbi, ddof=1) ** 2 * 0.5)
    out["SD2"] = np.sqrt(2 * np.std(bbi, ddof=1) ** 2 - 0.5 * np.std(diff_bbi, ddof=1) ** 2)
    out["SD2SD1"] = out["SD2"] / out["SD1"]

    # CSI / CVI
    #    T = 4 * out["SD1"]
    #    L = 4 * out["SD2"]
    #    out["CSI"] = L / T
    #    out["CVI"] = np.log10(L * T)
    #    out["CSI_Modified"] = L ** 2 / T

    # Entropy
    out["ApEn"] = entropy_approximate(bbi, dimension=2)
    out["SampEn"] = entropy_sample(bbi, dimension=2, r=0.2 * np.std(bbi, ddof=1))

    # DFA
    if len(bbi) / 10 > 16:
        out["DFA_1"] = fractal_dfa(bbi, windows=np.arange(4, 17))
    if len(bbi) > 65:
        out["DFA_2"] = fractal_dfa(bbi, windows=np.arange(16, 65))

    return out


# =============================================================================
# Internals
# ============================================================================= 

def plotstft(audiopath, binsize=2**10, plotpath=None, colormap="gray", channel=0, name='tmp.png', alpha=1, offset=0):
    samplerate, samples = wav.read(audiopath)
    samples = samples[:, channel]
    s = stft(samples, binsize)

    sshow, freq = logscale_spec(s, factor=1, sr=samplerate, alpha=alpha)
    sshow = sshow[2:, :]
    ims = 20.*np.log10(np.abs(sshow)/10e-6) # amplitude to decibel
    timebins, freqbins = np.shape(ims)
    
    ims = np.transpose(ims)
    ims = ims[0:256, offset:offset+768] # 0-11khz, ~9s interval
    #print "ims.shape", ims.shape
    
    image = Image.fromarray(ims) 
    image = image.convert('L')
    image.save(name) 

def search_c2x(self, zaxis, n):
        '''C2 axis which is perpendicular to z-axis'''
        decimals = int(-numpy.log10(TOLERANCE)) - 1
        for lst in self.group_atoms_by_distance:
            if len(lst) > 1:
                r0 = self.atoms[lst,1:]
                zcos = numpy.around(numpy.einsum('ij,j->i', r0, zaxis),
                                    decimals=decimals)
                uniq_zcos = numpy.unique(zcos)
                maybe_c2x = []
                for d in uniq_zcos:
                    if d > TOLERANCE:
                        mirrord = abs(zcos+d)<TOLERANCE
                        if mirrord.sum() == (zcos==d).sum():
                            above = r0[zcos==d]
                            below = r0[mirrord]
                            nelem = len(below)
                            maybe_c2x.extend([above[0] + below[i]
                                              for i in range(nelem)])
                    elif abs(d) < TOLERANCE: # plane which crosses the orig
                        r1 = r0[zcos==d][0]
                        maybe_c2x.append(r1)
                        r2 = numpy.dot(rotation_mat(zaxis, numpy.pi*2/n), r1)
                        if abs(r1+r2).sum() > TOLERANCE:
                            maybe_c2x.append(r1+r2)
                        else:
                            maybe_c2x.append(r2-r1)

                if len(maybe_c2x) > 0:
                    idx = norm(maybe_c2x, axis=1) > TOLERANCE
                    maybe_c2x = _normalize(maybe_c2x)[idx]
                    maybe_c2x = _remove_dupvec(maybe_c2x)
                    for c2x in maybe_c2x:
                        if (not parallel_vectors(c2x, zaxis) and
                            self.has_rotation(c2x, 2)):
                            return c2x 

def calc_radius_from_mass(self, Mp):
        """Forecast the Radius distribution given the mass distribution.
        
        Args:
            Mp (astropy Quantity array):
                Planet mass in units of Earth mass
        
        Returns:
            Rp (astropy Quantity array):
                Planet radius in units of Earth radius
        
        """
        
        mass = Mp.to('earthMass').value
        assert np.min(mass) > 3e-4 and np.max(mass) < 3e5, \
                "Mass range out of model expectation. Returning None."
        
        sample_size = len(mass)
        logm = np.log10(mass)
        prob = np.random.random(sample_size)
        logr = np.ones_like(logm)
        hyper_ind = np.random.randint(low=0, high=np.shape(self.all_hyper)[0], 
                size=sample_size)
        hyper = self.all_hyper[hyper_ind,:]
        
        for i in range(sample_size):
            logr[i] = self.piece_linear(hyper[i], logm[i], prob[i])
        
        Rp = 10.**logr*u.earthRad
        
        return Rp

#####################################################################################
# The following functions where adapted from the func.py file from the FORECASTER   #
# GitHub repository at https://github.com/chenjj2/forecaster (Chen & Kippling 2016) #
#####################################################################################