Output: m lpc coefficients, excitation energy
4181
float lpc_from_curve(float[] curve, float[] lpc){
    int n=ln;
    float[] work=new float[n+n];
    float fscale=(float)(.5/n);
    int i, j;

    // input is a real curve. make it complex-real
    // This mixes phase, but the LPC generation doesn't care.
    for(i=0; i<n; i++){
      work[i*2]=curve[i]*fscale;
      work[i*2+1]=0;
    }
    work[n*2-1]=curve[n-1]*fscale;

    n*=2;
    fft.backward(work);

    // The autocorrelation will not be circular.  Shift, else we lose
    // most of the power in the edges.

    for(i=0, j=n/2; i<n/2;){
      float temp=work[i];
      work[i++]=work[j];
      work[j++]=temp;
    }

    return (lpc_from_data(work, lpc, n, m));
  }
172570
def LML(self,params=None):
        
        if params is not None:
            self.setParams(params)

        KV = self._update_cache()
        alpha = KV['alpha']
        L = KV['L']

        lml_quad = 0.5 * (alpha*self.Y).sum()
        lml_det = self.t *SP.log(SP.diag(L)).sum()
        lml_const = 0.5*self.n*self.t*SP.log(2*SP.pi)
        LML = lml_quad + lml_det + lml_const
        return LML
40060
def _conversion(self, input_signal, vad_file):
    
    
    e = bob.ap.Energy(rate_wavsample[0], self.win_length_ms, self.win_shift_ms)
    energy_array = e(rate_wavsample[1])
    labels = self.use_existing_vad(energy_array, vad_file)
    
    return labels
1159939
def generate_synthObs(self, bases_wave, bases_flux, basesCoeff, Av_star, z_star, sigma_star, resample_range = None, resample_int = 1):
        
        ''''''
        nbases              = basesCoeff.shape[0]        
        bases_wave_resam    = arange(int(resample_range[0]), int(resample_range[-1]), resample_int, dtype=float)
        npix_resample       = len(bases_wave_resam)
        
        #Resampling the range
        bases_flux_resam = empty((nbases, npix_resample))
        for i in range(nbases):
#             print bases_wave[i][0], bases_wave[i][-1]
#             print bases_wave_resam[0], bases_wave_resam[-1]
            bases_flux_resam[i,:] = interp1d(bases_wave[i], bases_flux[i], bounds_error=True)(bases_wave_resam)            
        
        #Display physical parameters
        synth_wave                  = bases_wave_resam * (1 + z_star)
        
        Av_vector                   = Av_star * ones(nbases)
        Xx_redd                     = CCM89_Bal07(3.4, bases_wave_resam)   
        r_sigma                     = sigma_star/(synth_wave[1] - synth_wave[0])

        #Defining empty kernel
        box                         = int(3 * r_sigma) if int(3 * r_sigma) < 3 else 3
        kernel_len                  = 2 * box + 1
        kernel                      = zeros((1, kernel_len)) 
        kernel_range                = arange(0, 2 * box + 1)
        
        #Generating the kernel with sigma (the norm factor is the sum of the gaussian)
        kernel[0,:]                 = exp(-0.5 * ((square(kernel_range-box)/r_sigma)))
        norm                        = np_sum(kernel[0,:])        
        kernel                      = kernel / norm

        #Convove bases with respect to kernel for dispersion velocity calculation
        bases_grid_convolve         = convolve2d(bases_flux_resam, kernel, mode='same', boundary='symm')  
        
        #Interpolate bases to wavelength range
        interBases_matrix           = (interp1d(bases_wave_resam, bases_grid_convolve, axis=1, bounds_error=True)(bases_wave_resam)).T       

        #Generate final flux model including dust        
        dust_attenuation            = power(10, -0.4 * outer(Xx_redd, Av_vector))
        bases_grid_model            = interBases_matrix * dust_attenuation
                    
        #Generate combined flux
        synth_flux                  = np_sum(basesCoeff.T * bases_grid_model, axis=1)
                                      
        return synth_wave, synth_flux
300399
def _calc_ML(sampler, modelidx=0, e_range=None, e_npoints=100):
    

    ML, MLp, MLerr, ML_model = find_ML(sampler, modelidx)

    if e_range is not None:
        # prepare bogus data for calculation
        e_range = validate_array(
            "e_range", u.Quantity(e_range), physical_type="energy"
        )
        e_unit = e_range.unit
        energy = (
            np.logspace(
                np.log10(e_range[0].value),
                np.log10(e_range[1].value),
                e_npoints,
            )
            * e_unit
        )
        data = {
            "energy": energy,
            "flux": np.zeros(energy.shape) * sampler.data["flux"].unit,
        }
        modelout = sampler.modelfn(MLp, data)

        if isinstance(modelout, np.ndarray):
            blob = modelout
        else:
            blob = modelout[modelidx]

        if isinstance(blob, u.Quantity):
            modelx = data["energy"].copy()
            model_ML = blob.copy()
        elif len(blob) == 2:
            modelx = blob[0].copy()
            model_ML = blob[1].copy()
        else:
            raise TypeError("Model {0} has wrong blob format".format(modelidx))

        ML_model = (modelx, model_ML)

    return ML, MLp, MLerr, ML_model
1217436
def _iso_ic_on_monochromatic(
        electron_energy, seed_energy, seed_edensity, gamma_energy
    ):
        
        photE0 = (seed_energy / mec2).decompose().value
        phn = seed_edensity

        # electron_energy = electron_energy[:, None]
        gamma_energy = gamma_energy[:, None]
        photE0 = photE0[:, None, None]
        phn = phn[:, None, None]

        b = 4 * photE0 * electron_energy
        w = gamma_energy / electron_energy
        q = w / (b * (1 - w))
        fic = (
            2 * q * np.log(q)
            + (1 + 2 * q) * (1 - q)
            + (1.0 / 2.0) * (b * q) ** 2 * (1 - q) / (1 + b * q)
        )

        gamint = (
            fic
            * heaviside(1 - q)
            * heaviside(q - 1.0 / (4 * electron_energy ** 2))
        )
        gamint[np.isnan(gamint)] = 0.0

        if phn.size > 1:
            phn = phn.to(1 / (mec2_unit * u.cm ** 3)).value
            gamint = trapz_loglog(gamint * phn / photE0, photE0, axis=0)  # 1/s
        else:
            phn = phn.to(mec2_unit / u.cm ** 3).value
            gamint *= phn / photE0 ** 2
            gamint = gamint.squeeze()

        # gamint /= mec2.to('erg').value

        # r0 = (e**2 / m_e / c**2).to('cm')
        # sigt = ((8 * np.pi) / 3 * r0**2).cgs
        sigt = 6.652458734983284e-25
        c = 29979245800.0

        gamint *= (3.0 / 4.0) * sigt * c / electron_energy ** 2

        return gamint
1022972
def calc_signal_power(params, fname, transient=200, Df=None, mlab=True, NFFT=1000,
                      noverlap=0, window=plt.mlab.window_hanning):
    '''
    

    '''


    if type(fname) is str and os.path.isfile(fname):
        #open file
        f = h5py.File(fname)
        data = f['data'].value
        srate = f['srate'].value 
        tvec = np.arange(data.shape[1]) * 1000. / srate
        f.close()
    elif type(fname) is np.ndarray:
        data = fname
        srate = 1000.
        tvec = np.arange(data.shape[1]) * 1000. / srate
    else:
        raise Exception, '{} not a file or array'.format(fname)
    
    # slice
    slica = (tvec >= transient)
    data = data[:,slica]

    # subtract mean
    dataT = data.T - data.mean(axis=1)
    data = dataT.T
    
    #extract PSD
    PSD=[]
    for i in np.arange(len(params.electrodeParams['z'])):
        if mlab:
            Pxx, freqs=plt.mlab.psd(data[i], NFFT=NFFT,
                                    Fs=srate, noverlap=noverlap, window=window)
        else:
            [freqs, Pxx] = helpers.powerspec([data[i,]], tbin= 1.,
                                        Df=Df, pointProcess=False)
            mask = np.where(freqs >= 0.)
            freqs = freqs[mask]
            Pxx = Pxx.flatten()
            Pxx = Pxx[mask]
            Pxx = Pxx/tvec[tvec >= transient].size**2
        PSD +=[Pxx.flatten()]
        
    PSD=np.array(PSD)

    return freqs, PSD
977177
def defValueFuncs(self):
        '''
        
        '''
        MPCnvrs      = self.MPC**(-self.CRRA/(1.0-self.CRRA))
        vFuncNvrs    = LinearInterp(np.array([self.mNrmMin, self.mNrmMin+1.0]),np.array([0.0, MPCnvrs]))
        self.vFunc   = ValueFunc(vFuncNvrs,self.CRRA)
        self.vPfunc  = MargValueFunc(self.cFunc,self.CRRA)
348905
def stSpectralEntropy(X, n_short_blocks=10):
    """"""
    L = len(X)                         # number of frame samples
    Eol = numpy.sum(X ** 2)            # total spectral energy

    sub_win_len = int(numpy.floor(L / n_short_blocks))   # length of sub-frame
    if L != sub_win_len * n_short_blocks:
        X = X[0:sub_win_len * n_short_blocks]

    sub_wins = X.reshape(sub_win_len, n_short_blocks, order='F').copy()  # define sub-frames (using matrix reshape)
    s = numpy.sum(sub_wins ** 2, axis=0) / (Eol + eps)                      # compute spectral sub-energies
    En = -numpy.sum(s*numpy.log2(s + eps))                                    # compute spectral entropy

    return En
319803
def get_lps(self):
        
        if self.header is not None:
            for linguisticProcessor in self.header:
                for lp in linguisticProcessor:
                    yield lp
1109380
def _compute_Lp_matrix(self):
        
        self.Lp = sum([1/self.eigen_values[i]
                      * np.outer(self.eigen_basis[:, i], self.eigen_basis[:, i])
                      for i in range(1, self.eigen_values.size)])
398697
def PartialDynamicSystem(self, ieq, variable):
        

        if ieq == 0:

            #            if variable==self.physical_nodes[0].variable:
            # print('1')
            #                # U1 is output
            #                # U1=i1/pC+U2
            #                Uc=Variable(hidden=True)
            #                block1=ODE(self.variables[0],Uc,[1],[0,self.C])
            #                sub1=Sum([self.physical_nodes[1].variable,Uc],variable)
            #                return [block1,sub1]
            #            elif variable==self.physical_nodes[1].variable:
            #                print('2')
            #                # U2 is output
            #                # U2=U1-i1/pC
            #                Uc=Variable(hidden=True)
            #                block1=ODE(self.variables[0],Uc,[-1],[0,self.C])
            #                sum1=Sum([self.physical_nodes[0].variable,Uc],variable)
            #                return [block1,sum1]
            if variable == self.variables[0]:  # i1=(u1-u2)/Lp
                print('3')
                # i1 is output
                # i1=pC(U1-U2)
                Uc = Variable(hidden=True)
                subs1 = Subtraction(
                    self.physical_nodes[0].variable, self.physical_nodes[1].variable, Uc)
                block1 = ODE(Uc, variable, [1], [0, self.L])
                return [block1, subs1]
        elif ieq == 1:
            # i1=-i2
            if variable == self.variables[0]:
                # i1 as output
                return [Gain(self.variables[1], self.variables[0], -1)]
            elif variable == self.variables[1]:
                # i2 as output
                return [Gain(self.variables[0], self.variables[1], -1)]
281789
def _compute_anelastic_attenuation_term(self, C, rrup, mag):
        
        r = (rrup**2. + (C['c5'] * np.exp(C['c6'] * mag +
                                          C['c7'] * (8.5 - mag)**2.5))**2.)**.5
        f3 = ((C['c4'] + C['c13'] * mag) * np.log(r) +
              (C['c8'] + C['c12'] * mag) * r)
        return f3
813744
def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, target_mu, target_S):
        
        self._psi_computations(Z, mu, S)
        tmp = (self.inv_lengthscale2 * self._psi2[:, :, :, None]) / self._psi2_denom
        target_mu += -2.*(dL_dpsi2[:, :, :, None] * tmp * self._psi2_mudist).sum(1).sum(1)
        target_S += (dL_dpsi2[:, :, :, None] * tmp * (2.*self._psi2_mudist_sq - 1)).sum(1).sum(1)
29677
def energy_value(h: List[Union[int, float]],
                 J: Dict[Tuple[int, int], Union[int, float]],
                 sol: List[int]) -> Union[int, float]:
    
    ener_ising = 0
    for elm in J.keys():
        if elm[0] == elm[1]:
            raise TypeError("""Interaction term must connect two different variables""")
        else:
            ener_ising += J[elm] * int(sol[elm[0]]) * int(sol[elm[1]])
    for i in range(len(h)):
        ener_ising += h[i] * int(sol[i])
    return ener_ising
892613
def PartialDynamicSystem(self, ieq, variable):
        
#        print(ieq,variable.name)
        if ieq == 0:
            # U1-U2=R(i1)
            if variable == self.physical_nodes[0].variable:
                # U1 is output
                # U1=R(i1)+U2
                return [WeightedSum([self.physical_nodes[1].variable, self.variables[0]], variable, [1, self.R])]
            elif variable == self.physical_nodes[1].variable:
                # U2 is output
                # U2=-R(i1)+U2
                return [WeightedSum([self.physical_nodes[0].variable, self.variables[0]], variable, [1, -self.R])]
            elif variable == self.variables[0]:
                # i1 is output
                # i1=(U1-U2)/R
                return [WeightedSum([self.physical_nodes[0].variable, self.physical_nodes[1].variable], variable, [1/self.R, -1/self.R])]
        elif ieq == 1:
            # i1=-i2
            if variable == self.variables[0]:
                # i1 as output
                return [Gain(self.variables[1], self.variables[0], -1)]
            elif variable == self.variables[1]:
                # i2 as output
                return [Gain(self.variables[0], self.variables[1], -1)]
1370914
def lpgamma(p, a):
    
    res = p * (p - 1.0) / 4.0 * log(pi)
    for ii in range(1, p+1):
        res += gammaln(a + (1.0 - ii) / 2.0)

    return res