Return a SpinTaylorF2 waveform using CUDA to generate the phase and amplitude
345021
def spintaylorf2(**kwds):
    
    #####Pull out the input arguments#####
    f_lower = double(kwds['f_lower'])
    delta_f = double(kwds['delta_f'])
    distance = double(kwds['distance'])
    mass1 = double(kwds['mass1'])
    mass2 = double(kwds['mass2'])
    spin1x = double(kwds['spin1x'])
    spin1y = double(kwds['spin1y'])
    spin1z = double(kwds['spin1z'])
    phi0 = double(kwds['coa_phase'])               #Orbital Phase at coalescence
    phase_order = int(kwds['phase_order'])
    amplitude_order = int(kwds['amplitude_order'])
    inclination = double(kwds['inclination'])
    lnhatx = sin(inclination)
    lnhaty = 0.
    lnhatz = cos(inclination)
    psi = 0.

    tC= -1.0 / delta_f
    M = mass1 + mass2
    eta = mass1 * mass2 / (M * M)
    m_sec = M * lal.MTSUN_SI
    piM = lal.PI * m_sec

    vISCO = 1. / sqrt(6.)
    fISCO = vISCO * vISCO * vISCO / piM
    f_max = ceilpow2(fISCO)
    n = int(f_max / delta_f + 1)
    kmax = int(fISCO / delta_f)
    kmin = int(numpy.ceil(f_lower / delta_f))
    kmax = kmax if (kmax<n) else n

    #####Calculate the Orientation#####
    v0 = pow(piM *  kmin * delta_f,1./3)
    chi = sqrt(spin1x**2+spin1y**2+spin1z**2)
    kappa = (lnhatx*spin1x+lnhaty*spin1y+lnhatz*spin1z)/chi if (chi > 0.)  else 1.
    Jx0 = mass1*mass2*lnhatx/v0 + mass1*mass1*spin1x
    Jy0 = mass1*mass2*lnhaty/v0 + mass1*mass1*spin1y
    Jz0 = mass1*mass2*lnhatz/v0 + mass1*mass1*spin1z
    thetaJ = acos(Jz0 / sqrt(Jx0**2+Jy0**2+Jz0**2))
    psiJ = atan2(Jy0, -Jx0) # FIXME: check that Jy0 and Jx0 are not both 0
    # Rotate Lnhat back to frame where J is along z, to figure out initial alpha
    rotLx = lnhatx*cos(thetaJ)*cos(psiJ) - lnhaty*cos(thetaJ)*sin(psiJ) + lnhatz*sin(thetaJ)
    rotLy = lnhatx*sin(psiJ) + lnhaty*cos(psiJ)
    alpha0 = atan2(rotLy, rotLx) # FIXME: check that rotLy and rotLx are not both 0
    psiJ_P =psiJ + psi
    psiJ_C =psiJ + psi + lal.PI/4.

    #####Calculate the Coefficients#####
    #quadparam = 1.
    gamma0 = mass1*chi/mass2
    #Calculate the spin corrections
    # FIXME should use pycbc's function, but sigma has different expression
    # in Andy's code, double check
    # pn_beta, pn_sigma, pn_gamma = pycbc.pnutils.mass1_mass2_spin1z_spin2z_to_beta_sigma_gamma(
    #                               mass1, mass2, chi*kappa, 0) # FIXME: spin2 is taken to be 0
    pn_beta = (113.*mass1/(12.*M) - 19.*eta/6.)*chi*kappa
    pn_sigma = (  (5.*(3.*kappa*kappa-1.)/2.) + (7. - kappa*kappa)/96.  ) * (mass1*mass1*chi*chi/M/M)
    pn_gamma = (5.*(146597. + 7056.*eta)*mass1/(2268.*M) - 10.*eta*(1276. + 153.*eta)/81.)*chi*kappa
    prec_fac0 = 5.*(4. + 3.*mass2/mass1)/64.
    dtdv2 = 743./336. + 11.*eta/4.
    dtdv3 = -4.*lal.PI + pn_beta
    dtdv4 = 3058673./1016064. + 5429.*eta/1008. + 617.*eta*eta/144. - pn_sigma
    dtdv5 = (-7729./672.+13.*eta/8.)*lal.PI + 9.*pn_gamma/40.

    #####Calculate the Initial Euler Angles alpha_ref, beta_ref=0 and zeta_ref#####
    gam = gamma0*v0
    sqrtfac = sqrt(1. + 2.*kappa*gam + gam*gam)
    logv0 = log(v0)
    logfac1 = log(1. + kappa*gam + sqrtfac)
    logfac2 = log(kappa + gam + sqrtfac)
    v02 = v0 * v0
    v03 = v0 * v02
    kappa2 = kappa * kappa
    kappa3 = kappa2 * kappa
    gamma02 = gamma0 * gamma0
    gamma03 = gamma02 *gamma0

    alpha_ref = prec_fac0*(  logfac2 *( dtdv2*gamma0 + dtdv3*kappa - dtdv5*kappa/(2.*gamma02) + dtdv4/(2.*gamma0) - dtdv4*kappa2/(2.*gamma0) + (dtdv5*kappa3)/(2.*gamma02) )  +  logfac1*( - dtdv2*gamma0*kappa - dtdv3 + kappa*gamma03/2. - gamma03*kappa3/2. ) + logv0 *( dtdv2*gamma0*kappa + dtdv3 - kappa*gamma03/2. + gamma03*kappa3/2. ) + sqrtfac *( dtdv3 + dtdv4*v0/2. + dtdv5/gamma02/3. + dtdv4*kappa/(2.*gamma0) + dtdv5*kappa*v0/(6.*gamma0) - dtdv5*kappa2/(2.*gamma02) - 1/(3.*v03) - gamma0*kappa/(6.*v02) - dtdv2/v0 - gamma02/(3.*v0) + gamma02*kappa2/(2.*v0) + dtdv5*v02/3. ))  - alpha0

    zeta_ref = prec_fac0*( dtdv3*gamma0*kappa*v0 + dtdv4*v0 + logfac2 *(-dtdv2*gamma0 - dtdv3*kappa + dtdv5*kappa/(2.*gamma02) - dtdv4/(2.*gamma0) + dtdv4*kappa2/(2.*gamma0) - dtdv5*kappa3/(2.*gamma02) ) + logv0 *( kappa*gamma03/2. - gamma03*kappa3/2. ) + logfac1 *( dtdv2*gamma0*kappa + dtdv3 - kappa*gamma03/2. + gamma03*kappa3/2. ) - 1/(3.*v03) - gamma0*kappa/(2.*v02) - dtdv2/v0 + dtdv4*gamma0*kappa*v02/2. + dtdv5*v02/2. + sqrtfac *( -dtdv3 - dtdv4*v0/2. - dtdv5/(3.*gamma02) - dtdv4*kappa/(2.*gamma0) - dtdv5*kappa*v0/(6.*gamma0) + dtdv5*kappa2/(2.*gamma02) + 1/(3.*v03) + gamma0*kappa/(6.*v02) + dtdv2/v0 + gamma02/(3.*v0) - gamma02*kappa2/(2.*v0) - dtdv5*v02/3. ) + dtdv5*gamma0*kappa*v03/3. )

    #####Calculate the Complex sideband factors, mm=2 is first entry#####
    RE_SBfac0= (1.+cos(thetaJ)**2)/2.
    RE_SBfac1= sin(2.*thetaJ)
    RE_SBfac2= 3.*sin(thetaJ)**2
    RE_SBfac3= -sin(2.*thetaJ)
    RE_SBfac4= (1.+cos(thetaJ)**2)/2.
    IM_SBfac0= -cos(thetaJ)
    IM_SBfac1= -2.*sin(thetaJ)
    IM_SBfac2= 0.
    IM_SBfac3= -2.*sin(thetaJ)
    IM_SBfac4= cos(thetaJ)

    #####Calculate the PN terms # FIXME replace with functions in lalsimulation #####
    theta = -11831./9240.
    lambdaa = -1987./3080.0
    pfaN = 3.0/(128.0 * eta)
    pfa2 = 5.0*(743.0/84 + 11.0 * eta)/9.0
    pfa3 = -16.0*lal.PI + 4.0*pn_beta
    pfa4 = 5.0*(3058.673/7.056 + 5429.0/7.0 * eta + 617.0 * eta*eta)/72.0 - \
            10.0*pn_sigma
    pfa5 = 5.0/9.0 * (7729.0/84.0 - 13.0 * eta) * lal.PI - pn_gamma
    pfl5 = 5.0/3.0 * (7729.0/84.0 - 13.0 * eta) * lal.PI - pn_gamma * 3
    pfa6 = (11583.231236531/4.694215680 - 640.0/3.0 * lal.PI * lal.PI- \
            6848.0/21.0*lal.GAMMA) + \
            eta * (-15335.597827/3.048192 + 2255./12. * lal.PI * \
            lal.PI - 1760./3.*theta +12320./9.*lambdaa) + \
            eta*eta * 76055.0/1728.0 - \
            eta*eta*eta*  127825.0/1296.0
    pfl6 = -6848.0/21.0
    pfa7 = lal.PI * 5.0/756.0 * ( 15419335.0/336.0 + 75703.0/2.0 * eta - \
            14809.0 * eta*eta)

    FTaN = 32.0 * eta*eta / 5.0
    FTa2 = -(12.47/3.36 + 3.5/1.2 * eta)
    FTa3 = 4.0 * lal.PI
    FTa4 = -(44.711/9.072 - 92.71/5.04 * eta - 6.5/1.8 * eta*eta)
    FTa5 = -(81.91/6.72 + 58.3/2.4 * eta) * lal.PI
    FTa6 = (664.3739519/6.9854400 + 16.0/3.0 * lal.PI*lal.PI -
            17.12/1.05 * lal.GAMMA +
         (4.1/4.8 * lal.PI*lal.PI - 134.543/7.776) * eta -
         94.403/3.024 * eta*eta - 7.75/3.24 * eta*eta*eta)
    FTl6 = -8.56/1.05
    FTa7 = -(162.85/5.04 - 214.745/1.728 * eta - 193.385/3.024 * eta*eta) \
            * lal.PI

    dETaN = 2 * -eta/2.0
    dETa1 = 2 * -(3.0/4.0 + 1.0/12.0 * eta)
    dETa2 = 3 * -(27.0/8.0 - 19.0/8.0 * eta + 1./24.0 * eta*eta)
    dETa3 = 4 * -(67.5/6.4 - (344.45/5.76 - 20.5/9.6 * lal.PI*lal.PI) *
                             eta + 15.5/9.6 * eta*eta + 3.5/518.4 * eta*eta*eta)

    amp0 = -4. * mass1 * mass2 / (1.0e+06 * distance * lal.PC_SI ) * \
                    lal.MRSUN_SI * lal.MTSUN_SI * sqrt(lal.PI/12.0)

    htildeP = FrequencySeries(zeros(n,dtype=complex128), delta_f=delta_f, copy=False)
    htildeC = FrequencySeries(zeros(n,dtype=complex128), delta_f=delta_f, copy=False)
    spintaylorf2_kernel(htildeP.data[kmin:kmax], htildeC.data[kmin:kmax],
                        kmin, phase_order, amplitude_order, delta_f, piM, pfaN,
                        pfa2, pfa3, pfa4, pfa5, pfl5,
                        pfa6, pfl6, pfa7, FTaN, FTa2,
                        FTa3, FTa4, FTa5, FTa6,
                        FTl6, FTa7, dETaN, dETa1, dETa2,  dETa3,
                        amp0, tC, phi0,
                        kappa, prec_fac0, alpha_ref, zeta_ref,
                        dtdv2, dtdv3, dtdv4, dtdv5,
                        RE_SBfac0, RE_SBfac1, RE_SBfac2, RE_SBfac3, RE_SBfac4,
                        IM_SBfac0, IM_SBfac1, IM_SBfac2, IM_SBfac3, IM_SBfac4,
                        psiJ_P, psiJ_C, gamma0)

    return htildeP, htildeC
250631
def to_fock_representation(self, index_symbol='n', max_terms=None):
        """"""
        phase_factor = sympy.exp(
            sympy.Rational(-1, 2) * self.ampl * self.ampl.conjugate())
        if not isinstance(index_symbol, IdxSym):
            index_symbol = IdxSym(index_symbol)
        n = index_symbol
        if max_terms is None:
            index_range = IndexOverFockSpace(n, hs=self._hs)
        else:
            index_range = IndexOverRange(n, 0, max_terms-1)
        term = (
            (self.ampl**n / sympy.sqrt(sympy.factorial(n))) *
            BasisKet(FockIndex(n), hs=self._hs))
        return phase_factor * KetIndexedSum(term, index_range)
567326
def get_waveform_amplitude(mean_masks, mean_waveforms):
    """"""

    assert mean_waveforms.ndim == 2
    n_samples, n_channels = mean_waveforms.shape

    assert mean_masks.ndim == 1
    assert mean_masks.shape == (n_channels,)

    mean_waveforms = mean_waveforms * mean_masks
    assert mean_waveforms.shape == (n_samples, n_channels)

    # Amplitudes.
    m, M = mean_waveforms.min(axis=0), mean_waveforms.max(axis=0)
    return M - m
85828
def write_ndarray(self, result, dst_paths, nodata=None, compress='lzw'):
        """"""

        assert len(dst_paths) == result.shape[2]
        assert result.shape[0] == self._height
        assert result.shape[1] == self._width
        assert result.shape[2] == len(dst_paths)
        with rasterio.open(self._mrio._get_template_for_given_resolution(self._mrio.dst_res, "path")) as src_layer:
            pass # later we need src_layer for src_layer.window_transform(win)
        for i, pth in enumerate(dst_paths):
            dst_path_chunk = self.get_chunk_path_from_layer_path(pth, self.ji)

            result_layer_i = result[:, :, [i]]
            assert result_layer_i.shape[2] == 1
            kwargs = self._mrio._get_template_for_given_resolution(
                res=self._mrio.dst_res, return_="meta").copy()
            kwargs.update({"driver": "GTiff",
                           "compress": compress,
                           "nodata": nodata,
                           "height": self._height,
                           "width": self._width,
                           "dtype": result_layer_i.dtype,
                           "transform": src_layer.window_transform(self._window)})
            with rasterio.open(dst_path_chunk, "w", **kwargs) as dst:
                dst.write(result_layer_i[:, :, 0], 1)
1317462
def gpu_a_trous():
    

    ker1 = SourceModule("""
                        __global__ void gpu_a_trous_row_kernel(float *in1, float *in2, float *wfil, int *scale)
                        {
                            const int len = gridDim.x*blockDim.x;
                            const int col = (blockDim.x * blockIdx.x + threadIdx.x);
                            const int i = col;
                            const int row = (blockDim.y * blockIdx.y + threadIdx.y);
                            const int j = row*len;
                            const int tid2 = i + j;
                            const int lstp = exp2(float(scale[0] + 1));
                            const int sstp = exp2(float(scale[0]));

                            in2[tid2] = wfil[2]*in1[tid2];

                            if (row < lstp)
                                { in2[tid2] += wfil[0]*in1[col + len*(lstp - row - 1)]; }
                            else
                                { in2[tid2] += wfil[0]*in1[tid2 - lstp*len]; }

                            if (row < sstp)
                                { in2[tid2] += wfil[1]*in1[col + len*(sstp - row - 1)]; }
                            else
                                { in2[tid2] += wfil[1]*in1[tid2 - sstp*len]; }

                            if (row >= (len - sstp))
                                { in2[tid2] += wfil[3]*in1[col + len*(2*len - row - sstp - 1)]; }
                            else
                                { in2[tid2] += wfil[3]*in1[tid2 + sstp*len]; }

                            if (row >= (len - lstp))
                                { in2[tid2] += wfil[4]*in1[col + len*(2*len - row - lstp - 1)]; }
                            else
                                { in2[tid2] += wfil[4]*in1[tid2 + lstp*len]; }
                        }
                        """, keep=True)

    ker2 = SourceModule("""
                        __global__ void gpu_a_trous_col_kernel(float *in1, float *in2, float *wfil, int *scale)
                        {
                            const int len = gridDim.x*blockDim.x;
                            const int col = (blockDim.x * blockIdx.x + threadIdx.x);
                            const int i = col;
                            const int row = (blockDim.y * blockIdx.y + threadIdx.y);
                            const int j = row*len;
                            const int tid2 = i + j;
                            const int lstp = exp2(float(scale[0] + 1));
                            const int sstp = exp2(float(scale[0]));

                            in2[tid2] = wfil[2]*in1[tid2];

                            if (col < lstp)
                                { in2[tid2] += wfil[0]*in1[j - col + lstp - 1]; }
                            else
                                { in2[tid2] += wfil[0]*in1[tid2 - lstp]; }

                            if (col < sstp)
                                { in2[tid2] += wfil[1]*in1[j - col + sstp - 1]; }
                            else
                                { in2[tid2] += wfil[1]*in1[tid2 - sstp]; }

                            if (col >= (len - sstp))
                                { in2[tid2] += wfil[3]*in1[j + 2*len - sstp - col - 1]; }
                            else
                                { in2[tid2] += wfil[3]*in1[tid2 + sstp]; }

                            if (col >= (len - lstp))
                                { in2[tid2] += wfil[4]*in1[j + 2*len - lstp - col - 1]; }
                            else
                                { in2[tid2] += wfil[4]*in1[tid2 + lstp]; }

                        }
                       """, keep=True)

    return ker1.get_function("gpu_a_trous_row_kernel"), ker2.get_function("gpu_a_trous_col_kernel")
763641
def get_samples(self, sample_count):
        
        if self.amplitude.value <= 0:
            return None
        # Build samples by rolling the period cache through the buffer
        rolled_array = numpy.roll(self.wave_cache,
                                  -1 * self.last_played_sample)
        # Append remaining partial period
        full_count, remainder = divmod(sample_count, self.cache_length)
        final_subarray = rolled_array[:int(remainder)]
        return_array = numpy.concatenate((numpy.tile(rolled_array, full_count),
                                          final_subarray))
        # Keep track of where we left off to prevent popping between chunks
        self.last_played_sample = int(((self.last_played_sample + remainder) %
                                       self.cache_length))
        # Multiply output by amplitude
        return return_array * (self.amplitude.value *
                               self.amplitude_multiplier)
267653
def spa_tmplt(**kwds):
    
    # Pull out the input arguments
    f_lower = kwds['f_lower']
    delta_f = kwds['delta_f']
    distance = kwds['distance']
    mass1 = kwds['mass1']
    mass2 = kwds['mass2']
    s1z = kwds['spin1z']
    s2z = kwds['spin2z']
    phase_order = int(kwds['phase_order'])
    #amplitude_order = int(kwds['amplitude_order'])
    spin_order = int(kwds['spin_order'])

    if 'out' in kwds:
        out = kwds['out']
    else:
        out = None

    amp_factor = spa_amplitude_factor(mass1=mass1, mass2=mass2) / distance

    lal_pars = lal.CreateDict()
    if phase_order != -1:
        lalsimulation.SimInspiralWaveformParamsInsertPNPhaseOrder(
            lal_pars, phase_order)

    if spin_order != -1:
        lalsimulation.SimInspiralWaveformParamsInsertPNSpinOrder(
            lal_pars, spin_order)

    #Calculate the PN terms
    phasing = lalsimulation.SimInspiralTaylorF2AlignedPhasing(
                                    float(mass1), float(mass2),
                                    float(s1z), float(s2z),
                                    lal_pars)

    pfaN = phasing.v[0]
    pfa2 = phasing.v[2] / pfaN
    pfa3 = phasing.v[3] / pfaN
    pfa4 = phasing.v[4] / pfaN
    pfa5 = phasing.v[5] / pfaN
    pfa6 = (phasing.v[6] - phasing.vlogv[6] * log(4)) / pfaN
    pfa7 = phasing.v[7] / pfaN

    pfl5 = phasing.vlogv[5] / pfaN
    pfl6 = phasing.vlogv[6] / pfaN

    piM = lal.PI * (mass1 + mass2) * lal.MTSUN_SI

    kmin = int(f_lower / float(delta_f))

    vISCO = 1. / sqrt(6.)
    fISCO = vISCO * vISCO * vISCO / piM
    kmax = int(fISCO / delta_f)
    f_max = ceilpow2(fISCO)
    n = int(f_max / delta_f) + 1

    if not out:
        htilde = FrequencySeries(zeros(n, dtype=numpy.complex64), delta_f=delta_f, copy=False)
    else:
        if type(out) is not Array:
            raise TypeError("Output must be an instance of Array")
        if len(out) < kmax:
            kmax = len(out)
        if out.dtype != complex64:
            raise TypeError("Output array is the wrong dtype")
        htilde = FrequencySeries(out, delta_f=delta_f, copy=False)

    spa_tmplt_engine(htilde[kmin:kmax], kmin, phase_order, delta_f, piM, pfaN,
                     pfa2, pfa3, pfa4, pfa5, pfl5,
                     pfa6, pfl6, pfa7, amp_factor)
    return htilde
542356
def reduce_cuda(g, x, axes, dtype):
    """"""
    if axes != 0:
        raise NotImplementedError("'axes' keyword is not implemented for CUDA")

    return g(x, dtype=dtype)
1112935
def _print_neural_compute_device(cuda_gpus, use_mps, cuda_mem_req=None, has_mps_impl=True):
    
    num_cuda_gpus = len(cuda_gpus)
    if num_cuda_gpus >= 1:
        gpu_names = ', '.join(gpu['name'] for gpu in cuda_gpus)

    if use_mps:
        from ._mps_utils import mps_device_name
        print('Using GPU to create model ({})'.format(mps_device_name()))
    elif num_cuda_gpus >= 1:
        from . import _mxnet_utils
        plural = 's' if num_cuda_gpus >= 2 else ''
        print('Using GPU{} to create model ({})'.format(plural, gpu_names))
        if cuda_mem_req is not None:
            _mxnet_utils._warn_if_less_than_cuda_free_memory(cuda_mem_req, max_devices=num_cuda_gpus)
    else:
        import sys
        print('Using CPU to create model')
        if sys.platform == 'darwin' and _mac_ver() < (10, 14) and has_mps_impl:
            print('NOTE: If available, an AMD GPU can be leveraged on macOS 10.14+ for faster model creation')
178210
def tdot_blas(mat, out=None):
    """"""
    if (mat.dtype != 'float64') or (len(mat.shape) != 2):
        return np.dot(mat, mat.T)
    nn = mat.shape[0]
    if out is None:
        out = np.zeros((nn, nn))
    else:
        assert(out.dtype == 'float64')
        assert(out.shape == (nn, nn))
        # FIXME: should allow non-contiguous out, and copy output into it:
        assert(8 in out.strides)
        # zeroing needed because of dumb way I copy across triangular answer
        out[:] = 0.0

    # # Call to DSYRK from BLAS
    mat = np.asfortranarray(mat)
    out = blas.dsyrk(alpha=1.0, a=mat, beta=0.0, c=out, overwrite_c=1,
                     trans=0, lower=0)

    symmetrify(out, upper=True)
    return np.ascontiguousarray(out)
640617
def corr2(X):
    ''' 
    '''
    X = np.atleast_3d(X)
    N = X.shape[0]
    D = X.shape[2]

    if D == 1:
        return np.zeros(N, dtype=np.float)

    trii = np.triu_indices(D, k=1)
    DD = len(trii[0])
    r = np.zeros((N, DD))
    for i in np.arange(N):
        rmat = np.corrcoef(X[i])  # get the ith window from each signal, result will be DxD
        r[i] = rmat[trii]
    return r
951604
private void buildSpectrogram() {

        short[] amplitudes = wave.getSampleAmplitudes();
        int numSamples = amplitudes.length;

        int pointer = 0;
        // overlapping
        if (overlapFactor > 1) {
            int numOverlappedSamples = numSamples * overlapFactor;
            int backSamples = fftSampleSize * (overlapFactor - 1) / overlapFactor;
            short[] overlapAmp = new short[numOverlappedSamples];
            pointer = 0;
            for (int i = 0; i < amplitudes.length; i++) {
                overlapAmp[pointer++] = amplitudes[i];
                if (pointer % fftSampleSize == 0) {
                    // overlap
                    i -= backSamples;
                }
            }
            numSamples = numOverlappedSamples;
            amplitudes = overlapAmp;
        }
        // end overlapping

        numFrames = numSamples / fftSampleSize;
        framesPerSecond = (int) (numFrames / wave.length());

        // set signals for fft
        WindowFunction window = new WindowFunction();
        window.setWindowType("Hamming");
        double[] win = window.generate(fftSampleSize);

        double[][] signals = new double[numFrames][];
        for (int f = 0; f < numFrames; f++) {
            signals[f] = new double[fftSampleSize];
            int startSample = f * fftSampleSize;
            for (int n = 0; n < fftSampleSize; n++) {
                signals[f][n] = amplitudes[startSample + n] * win[n];
            }
        }
        // end set signals for fft

        absoluteSpectrogram = new double[numFrames][];
        // for each frame in signals, do fft on it
        FastFourierTransform fft = new FastFourierTransform();
        for (int i = 0; i < numFrames; i++) {
            absoluteSpectrogram[i] = fft.getMagnitudes(signals[i], false);
        }

        if (absoluteSpectrogram.length > 0) {

            numFrequencyUnit = absoluteSpectrogram[0].length;
            unitFrequency = (double) wave.getWaveHeader().getSampleRate() / 2 / numFrequencyUnit; // frequency could be caught within the half of nSamples according to Nyquist theory

            // normalization of absoultSpectrogram
            spectrogram = new double[numFrames][numFrequencyUnit];

            // set max and min amplitudes
            double maxAmp = Double.MIN_VALUE;
            double minAmp = Double.MAX_VALUE;
            for (int i = 0; i < numFrames; i++) {
                for (int j = 0; j < numFrequencyUnit; j++) {
                    if (absoluteSpectrogram[i][j] > maxAmp) {
                        maxAmp = absoluteSpectrogram[i][j];
                    } else if (absoluteSpectrogram[i][j] < minAmp) {
                        minAmp = absoluteSpectrogram[i][j];
                    }
                }
            }
            // end set max and min amplitudes

            // normalization
            // avoiding divided by zero 
            double minValidAmp = 0.00000000001F;
            if (minAmp == 0) {
                minAmp = minValidAmp;
            }

            double diff = Math.log10(maxAmp / minAmp); // perceptual difference
            for (int i = 0; i < numFrames; i++) {
                for (int j = 0; j < numFrequencyUnit; j++) {
                    if (absoluteSpectrogram[i][j] < minValidAmp) {
                        spectrogram[i][j] = 0;
                    } else {
                        spectrogram[i][j] = (Math.log10(absoluteSpectrogram[i][j] / minAmp)) / diff;
                    }
                }
            }
            // end normalization
        }
    }
124553
def cublasCsyr(handle, uplo, n, alpha, x, incx, A, lda):
    

    status = _libcublas.cublasCsyr_v2(handle,
                                        _CUBLAS_FILL_MODE[uplo], n, 
                                        ctypes.byref(cuda.cuFloatComplex(alpha.real,
                                                                        alpha.imag)),
                                        int(x), incx, int(A), lda)
    cublasCheckStatus(status)
1112307
def cudaGetDevice():
    

    dev = ctypes.c_int()
    status = _libcudart.cudaGetDevice(ctypes.byref(dev))
    cudaCheckStatus(status)
    return dev.value
1137342
def _tv2(data,weight,Niter=50):
    

    if dev is None:
        dev = imgtools.__DEFAULT_OPENCL_DEVICE__

    if dev is None:
        raise ValueError("no OpenCLDevice found...")


    proc = OCLProcessor(dev,utils.absPath("kernels/tv_chambolle.cl"))

    if Ncut ==1:
        inImg = dev.createImage(data.shape[::-1],dtype = np.float32)

        pImgs = [ dev.createImage(data.shape[::-1],
                                  mem_flags = cl.mem_flags.READ_WRITE,
                                  dtype= np.float32,
                                  channel_order = cl.channel_order.RGBA)
                                  for i in range(2)]

        outImg = dev.createImage(data.shape[::-1],
                                 dtype = np.float32,
                                 mem_flags = cl.mem_flags.READ_WRITE)


        dev.writeImage(inImg,data.astype(np.float32));
        dev.writeImage(pImgs[0],np.zeros((4,)+data.shape,dtype=np.float32));
        dev.writeImage(pImgs[1],np.zeros((4,)+data.shape,dtype=np.float32));


        for i in range(Niter):
            proc.runKernel("div_step",inImg.shape,None,
                           inImg,pImgs[i%2],outImg)
            proc.runKernel("grad_step",inImg.shape,None,
                           outImg,pImgs[i%2],pImgs[1-i%2],
                           np.float32(weight))
        return dev.readImage(outImg,dtype=np.float32)

    else:
        res = np.empty_like(data,dtype=np.float32)
        Nz,Ny,Nx = data.shape
        # a heuristic guess: Npad = Niter means perfect
        Npad = 1+Niter/2
        for i0,(i,j,k) in enumerate(product(list(range(Ncut)),repeat=3)):
            logger.info("calculating box  %i/%i"%(i0+1,Ncut**3))
            sx = slice(i*Nx/Ncut,(i+1)*Nx/Ncut)
            sy = slice(j*Ny/Ncut,(j+1)*Ny/Ncut)
            sz = slice(k*Nz/Ncut,(k+1)*Nz/Ncut)
            sx1,sx2 = utils._extended_slice(sx,Nx,Npad)
            sy1,sy2 = utils._extended_slice(sy,Ny,Npad)
            sz1,sz2 = utils._extended_slice(sz,Nz,Npad)

            data_sliced = data[sz1,sy1,sx1].copy()
            _res = tv3_gpu(dev,data_sliced,weight,Niter,Ncut = 1)
            res[sz,sy,sx] = _res[sz2,sy2,sx2]

        return res
744112
def phaser(self,
               gain_in=0.9,
               gain_out=0.8,
               delay=1,
               decay=0.25,
               speed=2,
               triangular=False):
        """"""
        self.command.append("phaser")
        self.command.append(gain_in)
        self.command.append(gain_out)
        self.command.append(delay)
        self.command.append(decay)
        self.command.append(speed)
        if triangular:
            self.command.append('-t')
        else:
            self.command.append('-s')
        return self
719273
def cublasCgbmv(handle, trans, m, n, kl, ku, alpha, A, lda,
                x, incx, beta, y, incy):
    

    status = _libcublas.cublasCgbmv_v2(handle,
                                       trans, m, n, kl, ku,
                                       ctypes.byref(cuda.cuFloatComplex(alpha.real,
                                                                        alpha.imag)),
                                       int(A), lda, int(x), incx,
                                       ctypes.byref(cuda.cuFloatComplex(beta.real,
                                                                        beta.imag)),
                                       int(y), incy)
    cublasCheckStatus(status)