Python scipy.squeeze函数代码示例

def interpolateBetweenBinaryObjects(obj1, obj2, slices):
    """
    Takes two binary objects and puts slices slices in-between them, each of which
    contains a smooth binary transition between the objects.
    """
    # constants
    temporal_dimension = 3
    # flip second returned binary objects along temporal axis
    slicer = [slice(None) for _ in range(obj1.ndim + 1)]
    slicer[temporal_dimension] = slice(None, None, -1)
    # logical-and combination
    ret = (
        __interpolateBetweenBinaryObjects(obj1, obj2, slices)
        | __interpolateBetweenBinaryObjects(obj2, obj1, slices)[slicer]
    )
    # control step: see if last volume corresponds to obj2
    slicer[temporal_dimension] = slice(-1, None)
    if not scipy.all(scipy.squeeze(ret[slicer]) == obj2.astype(scipy.bool_)):
        raise Exception(
            "Last created object does not correspond to obj2. Difference of {} voxels.".format(
                len(scipy.nonzero(scipy.squeeze(ret[slicer]) & obj2.astype(scipy.bool_))[0])
            )
        )

    return ret

def Au(U,GF,EpsArr,NX,NY,NZ):
    """Returns the result of matrix-vector multiplication
       by the system matrix A=I-GX
    """
    # reshaping input vector into 4-D array
    Uarr=sci.reshape(U,(NX,NY,NZ,3))
    # extended zero-padded arrays
    Uext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    Vext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    Jext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    JFext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    Uext[0:NX,0:NY,0:NZ,:]=Uarr
    # contrast current array
    s=0
    while s<=2:
        Jext[0:NX,0:NY,0:NZ,s]=Uext[0:NX,0:NY,0:NZ,s]*(EpsArr[0:NX,0:NY,0:NZ]-1.0)
        JFext[:,:,:,s]=fft.fftn(sci.squeeze(Jext[:,:,:,s]))
        s=s+1
    Vext[:,:,:,0]=Uext[:,:,:,0]-\
    fft.ifftn(sci.squeeze(sci.multiply(GF[:,:,:,0,0],JFext[:,:,:,0])+\
                          sci.multiply(GF[:,:,:,0,1],JFext[:,:,:,1])+\
                          sci.multiply(GF[:,:,:,0,2],JFext[:,:,:,2])))
    Vext[:,:,:,1]=Uext[:,:,:,1]-\
    fft.ifftn(sci.squeeze(sci.multiply(GF[:,:,:,1,0],JFext[:,:,:,0])+\
                          sci.multiply(GF[:,:,:,1,1],JFext[:,:,:,1])+\
                          sci.multiply(GF[:,:,:,1,2],JFext[:,:,:,2])))
    Vext[:,:,:,2]=Uext[:,:,:,2]-\
    fft.ifftn(sci.squeeze(sci.multiply(GF[:,:,:,2,0],JFext[:,:,:,0])+\
                          sci.multiply(GF[:,:,:,2,1],JFext[:,:,:,1])+\
                          sci.multiply(GF[:,:,:,2,2],JFext[:,:,:,2])))
    # reshaping output into column vector
    V=sci.reshape(Vext[0:NX,0:NY,0:NZ,:],(NX*NY*NZ*3,1))

    return V

 def dB_Phase_Error_vect(self,other):
     """Subtract two bodes dBmag and phase and return a Bode with
     dBmag and phase."""
     outbode=rwkbode()
     outbode.phase=squeeze(self.phase)-squeeze(other.phase)
     outbode.dBmag=self.dBmag()-other.dBmag()
     return outbode

 def __mul__(self,other):
     """Multiply two bodes."""
     if type(other)==float or type(other)==int:
         myoutput=copy.deepcopy(self)
         myoutput.mag=myoutput.mag*other
         return myoutput
     myin='in'
     myout='out'
     match=1
     if self.output==other.input:
         first=self
         second=other
     elif self.input==other.output:
         first=other
         second=self
     else:
         warnme=1
         if (self.input=='in' and self.output=='out') or (other.input=='in' and other.output=='out'):
             warnme=0
         if warnme:
             print('Warning: multiplying Bodes without a matching input/output pair:\n'+self.output+'/'+self.input+ ' * ' +other.output+'/'+other.input)
         match=0
         first=self
         second=other
     if match:
         myin=first.input
         myout=first.output
     myoutput=copy.deepcopy(self)
     myoutput.input=myin
     myoutput.output=myout
     myoutput.mag=squeeze(colwise(first.mag)*colwise(second.mag))
     myoutput.phase=squeeze(colwise(first.phase)+colwise(second.phase))
     return myoutput

def __load_images(label_image_n, mask_image_n, boundary_image_n, fg_image_n = False, bg_image_n = False):
    """
    Load and return all image data in preprocessed ndarrays.
    The label image will be relabeled to start from 1.
    @return label, ground-truth-mask, boundary-result-mask[, fg-markers, bg-markers]
    """
    # load images
    label_image = load(label_image_n)
    mask_image = load(mask_image_n)
    boundary_image = load(boundary_image_n)
    if fg_image_n: fg_image = load(fg_image_n)
    if bg_image_n: bg_image = load(bg_image_n)
    
    # extract image data
    label_image_d = scipy.squeeze(label_image.get_data())
    mask_image_d = scipy.squeeze(mask_image.get_data()).astype(scipy.bool_)
    boundary_image_d = scipy.squeeze(boundary_image.get_data()).astype(scipy.bool_)
    if fg_image_n: fg_image_d = scipy.squeeze(fg_image.get_data()).astype(scipy.bool_)
    if bg_image_n: bg_image_d = scipy.squeeze(bg_image.get_data()).astype(scipy.bool_)
    
    # check if images are of same dimensionality
    if label_image_d.shape != mask_image_d.shape: raise argparse.ArgumentError('The mask image {} must be of the same dimensionality as the label image {}.'.format(mask_image_d.shape, label_image_d.shape))
    if label_image_d.shape != boundary_image_d.shape: raise argparse.ArgumentError('The boundary term image {} must be of the same dimensionality as the label image {}.'.format(boundary_image_d.shape, label_image_d.shape))
    if fg_image_n:
        if label_image_d.shape != fg_image_d.shape: raise argparse.ArgumentError('The foreground markers image {} must be of the same dimensionality as the label image {}.'.format(fg_image_d.shape, label_image_d.shape))
    if bg_image_n:
        if label_image_d.shape != bg_image_d.shape: raise argparse.ArgumentError('The background markers image {} must be of the same dimensionality as the label image {}.'.format(bg_image_d.shape, label_image_d.shape))
    
    # relabel the label image to start from 1
    label_image_d = medpy.filter.relabel(label_image_d, 1)
    
    if fg_image_n:
        return label_image_d, mask_image_d, boundary_image_d, fg_image_d, bg_image_d
    else:
        return label_image_d, mask_image_d, boundary_image_d

def CombinedBodes(bodelist):
    """This function seeks to make one combined rwkbode instance from
    a list of them.  This may be useful in the case of having several
    experimental Bode data sets with various targeted frequency ranges
    that need to be treated as one data set."""
    docoh=True
    for cb in bodelist:
        if cb.coh is None:
            docoh=False
            break
    first=1
#    pdb.set_trace()
    bigcoh=[]
    for cb in bodelist:
        if first:
            first=0
            bigmag=colwise(cb.mag)
            bigphase=colwise(cb.phase)
            if docoh:
                bigcoh=colwise(cb.coh)
        else:
            bigmag=r_[bigmag,colwise(cb.mag)]
            bigphase=r_[bigphase,colwise(cb.phase)]
            if docoh:
                bigcoh=r_[bigcoh,colwise(cb.coh)]
    myoutbode=copy.deepcopy(bodelist[0])
    myoutbode.mag=squeeze(bigmag)
    myoutbode.phase=squeeze(bigphase)
    myoutbode.coh=squeeze(bigcoh)
    myoutbode.freqlim=[]
    return myoutbode

 def __init__(self, output='out', input='in', \
              mag=None, phase=None, coh=None, \
              freqlim=[], maglim=[], phaselim=[], \
              averaged='not specified', \
              seedfreq=-1, seedphase=0,
              labels=[], legloc=-1, compin=[]):
     self.output = output
     self.input = input
     if len(compin) > 0:
        if mag is None:
           self.mag = squeeze(colwise(abs(compin)))
        if phase is None:
           self.phase = squeeze(colwise(arctan2(imag(compin),real(compin))*180.0/pi))
     else:
         self.mag = squeeze(mag)
         self.phase = squeeze(phase)
     self.coh = coh
     self.averaged = averaged
     self.seedfreq = seedfreq
     self.seedphase = seedphase
     self.freqlim = freqlim
     self.maglim = maglim
     self.phaselim = phaselim
     self.labels = labels
     self.legloc = legloc

 def dB_Phase_Error_sum(self,other, phaseweight=0.1):
     """Subtract two bodes dBmag and phase and return the sum of
     the squared error."""
     phaseerror = squeeze(self.phase)-squeeze(other.phase)
     e_phase = (phaseerror**2).sum()
     dBmagerror = self.dBmag()-other.dBmag()
     e_dB = (dBmagerror**2).sum()
     return e_dB + e_phase*phaseweight

 def ToComp(self,ind=None):
     if ind!=None:
         x=self.mag[ind]*cos(pi/180.0*self.phase[ind])
         y=self.mag[ind]*sin(pi/180.0*self.phase[ind])
         return x+y*1.0j
     else:
         x=squeeze(self.mag)*squeeze(cos(pi/180.0*self.phase))
         y=squeeze(self.mag)*squeeze(sin(pi/180.0*self.phase))
     return squeeze(x+y*1.0j)

def AveBodeFromSpectra(specin):
    Gxxave=squeeze(specin.Gxx.mean(1))
    Gxyave=squeeze(specin.Gxy.mean(1))
    Gyyave=squeeze(specin.Gyy.mean(1))
    Have=Gxyave/Gxxave
    cohnum=norm2(Gxyave)
    cohden=Gxxave*Gyyave
    mybode=rwkbode(input=specin.input, output=specin.output, \
                   compin=squeeze(Have), coh=squeeze(cohnum/cohden))
    return mybode

def RenormalizeFactor(excfile,gsfile,channel=None,Nsamp=1,O=None,q=None):
    if not type(excfile)==list:
        excfile=[excfile]
    if not type(gsfile)==list:
        gsfile=[gsfile]
    exat=GetAttr(excfile[0])
    gsat=GetAttr(gsfile[0])
    L=exat['L']
    if q==None:
        q=sc.array([exat['qx'],exat['qy']])
    if 'phasex' in exat.keys():
        shift=sc.array([exat['phasex']/2.0,exat['phasey']/2.0])
    else:
        shift=sc.array([exat['phase_shift_x']/2.0,exat['phase_shift_y']/2.0])
    phi=exat['phi']
    neel=exat['neel']
    qx,qy,Sq=GetSq(gsfile)
    kx,ky=sf.fermisea(L,L,shift)
    qidx=ml.find((qx==q[0])*(qy==q[1]))
    if O==None:
        _,O,_,_=GetEigSys(excfile,Nsamp)
    pk=None
    sqq=None
    if channel==None:
        channel=exat['channel']
    if channel=='trans':
        pk=sc.squeeze(sf.phiktrans(kx,ky,q[0]/L,q[1]/L,[phi,neel]))
        sqq=sc.real(0.5*(Sq[0,1,qidx]+Sq[0,2,qidx]))
    elif channel=='long':
        pkup=sc.squeeze(sf.phiklong(kx,ky,q[0]/L,q[1]/L,1,[phi,neel]))
        pkdo=sc.squeeze(sf.phiklong(kx,ky,q[0]/L,q[1]/L,-1,[phi,neel]))
        if (q[0]/L==0.5 and q[1]/L==0.5) or (q[0]/L==0 and q[1]/L==0):
            pk=sc.zeros(2*sc.shape(pkup)[0]+1,complex)
        else:
            pk=sc.zeros(2*sc.shape(pkup)[0],complex)
        pk[0:2*sc.shape(pkup)[0]:2]=pkup
        pk[1:2*sc.shape(pkdo)[0]:2]=pkdo
        if (qx[0]/L==0.5 and q[1]/L==0.5) or (q[0]/L==0 and q[1]/L==0):
            if neel==0:
                pk[-1]=0
            else:
                pk[-1]=sum(neel/sf.omega(kx,ky,[phi,neel]))
        sqq=Sq[0,0,qidx]
    else:
        raise(InputFileError('In file \''+excfile+'\', channel=\''+str(channel)+'\'. Should be \'trans\' or \'long\''))
    sqe=sc.einsum('i,jik,k->j',sc.conj(pk),O,pk)
    out=sc.zeros(Nsamp)
    for n in range(Nsamp):
        if abs(sqq)<1e-6 or abs(sqe[n])<1e-6:
            warnings.warn('Probably ill-defined renormalization, returns 1 for sample {0} out of {1}'.format(n,Nsamp),UserWarning)
            out[n]=1
        else:
            out[n]=sc.real(sqq/sqe[n])
    return out

def GetEigSys(filename,gsfile=None,Nsamp=1,channel=None,wavefile=None,q=None):
    if type(filename)==str:
        filename=[filename]
    hfile=h5py.File(filename[0],'r')
    attr=GetAttr(filename[0])
    if channel==None:
        channel=attr['channel']
    dpath,args=GetStat(filename,Nsamp)
    dat=sc.array(hfile["/rank-1/data-0"])
    hfile.close()
    N=int(sc.shape(dat)[0]/2)
    L=attr['L']
    shift=None
    if 'phasex' in attr.keys():
        shift=[attr['phasex']/2.0,attr['phasey']/2.0]
    else:
        shift=[attr['phase_shift_x']/2.0,attr['phase_shift_y']/2.0]
    H=sc.zeros([Nsamp,N,N],complex)
    O=sc.zeros([Nsamp,N,N],complex)
    E=sc.zeros([Nsamp,N])
    V=sc.zeros([Nsamp,N,N],complex)
    for sample,b in enumerate(args):
        for d in b:
            hfile=h5py.File(dpath[d][0],'r')
            dat=hfile[dpath[d][1]]
            H[sample,:,:]+=dat[0:N,0:2*N:2]+1j*dat[0:N,1:2*N:2]
            O[sample,:,:]+=dat[N:2*N,0:2*N:2]+1j*dat[N:2*N,1:2*N:2]
            hfile.close()
        H[sample,:,:]=0.5*(H[sample,:,:]+sc.conj(H[sample,:,:].T))/len(b)
        O[sample,:,:]=0.5*(O[sample,:,:]+sc.conj(O[sample,:,:].T))/len(b)
    if channel=='groundstate':
        return H
    fs=None
    refstate=sc.zeros(2*L*L)
    refstate[0::2]=1
    if wavefile==None:
        fs=GetFermiSigns(filename[0],refstate,channel=channel)
    else:
        fs=GetFermiSigns(wavefile,refstate,channel=channel)
    for s in range(sc.shape(H)[0]):
        H[s,:,:]=sc.dot(sc.diag(fs),sc.dot(H[s,:,:],sc.diag(fs)))
        O[s,:,:]=sc.dot(sc.diag(fs),sc.dot(O[s,:,:],sc.diag(fs)))
    ren=sc.ones(Nsamp)
    if gsfile!=None:
        ren=RenormalizeFactor(filename,gsfile,Nsamp=1,channel=channel,O=O,q=q)
    print('{0} pair of (H,O) matrices loaded, now diagonalize'.format(sc.shape(H)[0]))
    H=sc.einsum('ijk,i->ijk',H,ren)
    O=sc.einsum('ijk,i->ijk',O,ren)
    for s in range(sc.shape(H)[0]):
        E[s,:],V[s,:,:]=vln.geneigh(sc.squeeze(H[s,:,:]),sc.squeeze(O[s,:,:]))
    print('diagonalization finished')
    return H,O,E,V

    def return_results(self, pores=None, throats=None, **kwargs):
        r"""
        Send results of simulation out the the appropriate locations.

        This is a basic version of the update that simply sends out the main
        result (quantity). More elaborate updates should be subclassed.
        """
        if pores is None:
            pores = self.Ps
        if throats is None:
            throats = self.Ts

        phase_quantity = self._quantity.replace(self._phase.name + '_', '')
        if phase_quantity not in self._phase.props():
            self._phase[phase_quantity] = sp.nan
        self._phase[phase_quantity][pores] = self[self._quantity][pores]
        conn_arr = self._net.find_connected_pores(self.Ts)
        dx = sp.squeeze(sp.diff(self[self._quantity][conn_arr], n=1, axis=1))
        g = self['throat.conductance']
        rate = sp.absolute(g * dx)
        if 'throat.rate' not in self._phase.props():
            self._phase['throat.rate'] = sp.nan
        self._phase['throat.rate'][throats] = rate[throats]
        logger.debug('Results of ' + self.name +
                     ' algorithm have been added to ' + self._phase.name)

def BodeFromSpectra(specin,calccoh=False):
    H=specin.Gxy/specin.Gxx
    if calccoh:
       cohnum=squeeze(norm2(specin.Gxy))
       cohden=squeeze(specin.Gxx*specin.Gyy)
       coh=cohnum/cohden
       ave=True
    else:
       coh=[]
       ave=False
    mybode=rwkbode(input=specin.input, output=specin.output, \
                   compin=squeeze(H), coh=coh, averaged=ave)
    mybode.Gxx=specin.Gxx
    mybode.Gyy=specin.Gyy
    mybode.Gxy=specin.Gxy
    return mybode

def plotVol(volume,pts=15,**kwargs):
    fluxGrid = scipy.squeeze(volume.getFluxGrid()).T

    datain = scipy.zeros((fluxGrid.shape[0],
                          pts,
                          fluxGrid.shape[1]))

    for idx in range(pts):
        datain[:,idx,:] = fluxGrid

        temp = genCylGrid(plasma.eq.getRGrid(),
                      scipy.linspace(0,2*scipy.pi,pts),
                      plasma.eq.getZGrid())
    verticies = genVertsFromPixel(temp)

    hex_type = tvtk.api.tvtk.Hexahedron().cell_type
    temp = temp.reshape((temp.size/3,3))
    verticies = verticies.reshape((verticies.size/8,8))

    sg = tvtk.api.tvtk.UnstructuredGrid(points=temp)
    
    sg.set_cells(hex_type,verticies)
    sg.point_data.scalars = datain.flatten()
    sg.point_data.scalars.name = 'temp'
    
    psi_0 = plasma.eq.getFluxAxis()
    psi_LCFS = plasma.eq.getFluxLCFS()
    psi_min = fluxGrid.min()
    psi_max = fluxGrid.max()
    v1 = (psi_0 - psi_min)/(psi_max - psi_min)
    v2 = (psi_LCFS - psi_min)/(psi_max - psi_min)
    mlab.pipeline.volume(sg)