Python pygimli.cat函数代码示例

    def block1dInversionNew(self, nlay=2, lam=100.0, verbose=True):
        """invert all data together by a 1D model (more general solution)"""
        data, error = pg.RVector(), pg.RVector()
        for mrs in self.mrs:
            data = pg.cat(data, mrs.data)
            error = pg.cat(error, mrs.error)

        f = JointMRSModelling(self.mrs, nlay)
        mrsobj = self.mrs[0]
        for i in range(3):
            f.region(i).setParameters(mrsobj.startval[i], mrsobj.lowerBound[i], mrsobj.upperBound[i])
        #        f.region(0).setStartValue(mrsobj.startval[0])
        #        f.region(1).setStartValue(mrsobj.startval[1])
        #        f.region(2).setStartValue(mrsobj.startval[2])
        #        # Model transformation instances saved in class
        #        transTH = pg.RTransLogLU(mrsobj.lowerBound[0], mrsobj.upperBound[0])
        #        transWC = pg.RTransLogLU(mrsobj.lowerBound[1], mrsobj.upperBound[1])
        #        transT2 = pg.RTransLogLU(mrsobj.lowerBound[2], mrsobj.upperBound[2])
        #        f.region(0).setTransModel(transTH)
        #        f.region(1).setTransModel(transWC)
        #        f.region(2).setTransModel(transT2)
        INV = pg.RInversion(data, f, verbose)
        INV.setLambda(lam)
        INV.setMarquardtScheme(0.8)
        INV.stopAtChi1(False)  # now in MarquardtScheme
        INV.setDeltaPhiAbortPercent(0.5)
        INV.setAbsoluteError(error)
        model = INV.run()
        return model

 def setData(self, data):
     """
     """
     if type(data) is list:
         if len(data) == len(self.managers):
             self.tD.clear()
             self.dataVals.clear()
             self.dataErrs.clear()
             
             self.fop.setData(data)
             
             for i, mgr in enumerate(self.managers):
                 t = mgr.tD
                 self.tD.add(t, data[i].size())
 
                 self.dataVals = pg.cat(self.dataVals, data[i](mgr.dataToken()))
         
                 if mgr.errIsAbsolute:
                     self.dataErrs = pg.cat(self.dataErrs, data[i]('err')/data[i](mgr.dataToken()))
                 else:
                     self.dataErrs = pg.cat(self.dataErrs, data[i]('err'))
             
             self.data = data
             
             self.inv.setTransData(self.tD)
             self.inv.setTransModel(self.tM)
         else:
             raise BaseException("To few datacontainer given")

    def runEA(self,nlay=None,type='GA',pop_size=100,max_evaluations=10000,**kwargs):
        import inspyred
        import random
        
        def mygenerate( random, args ):
            """ generate a random vector of model size """
            return [random.random() for i in range( nlay*3 - 1 )]
        
        def my_observer(population, num_generations, num_evaluations, args):
            best = min(population)
            print('{0:6} -- {1}'.format(num_generations,best.fitness))
        
        @inspyred.ec.evaluators.evaluator
        def datafit( individual, args ):
            misfit = (self.data-self.f.response(self.genMod(individual)))/self.error
            return np.mean(misfit**2)
        
        # prepare forward operator
        if self.f is None or (nlay is not None and nlay is not self.nlay): self.createFOP(nlay)
        
        lowerBound = pg.cat( pg.cat( pg.RVector(self.nlay-1,self.lowerBound[0]), 
            pg.RVector(self.nlay,self.lowerBound[1])), pg.RVector(self.nlay,self.lowerBound[2]) )
        upperBound = pg.cat( pg.cat( pg.RVector(self.nlay-1,self.upperBound[0]), 
            pg.RVector(self.nlay,self.upperBound[1])), pg.RVector(self.nlay,self.upperBound[2]) )
        if self.logpar:
            self.lLB, self.lUB = pg.log(lowerBound), pg.log(upperBound) # ready mapping functions
        else:
            self.lLB, self.lUB = lowerBound, upperBound
        
#        self.f = MRS1dBlockQTModelling(nlay, self.K, self.z, self.t)
        # setup random generator
        rand = random.Random()
        rand.seed(int(time.time()))
        # choose among different evolution algorithms
        if type == 'GA': 
            ea = inspyred.ec.GA(rand)
            ea.variator = [inspyred.ec.variators.blend_crossover, inspyred.ec.variators.gaussian_mutation]
            ea.selector = inspyred.ec.selectors.tournament_selection
            ea.replacer = inspyred.ec.replacers.generational_replacement
        if type == 'SA': ea = inspyred.ec.SA(rand)
        if type == 'DEA': ea = inspyred.ec.DEA(rand)
        if type == 'PSO': ea = inspyred.swarm.PSO(rand)
        if type == 'ACS': ea = inspyred.swarm.ACS(rand,[])
        if type == 'ES': 
            ea = inspyred.ec.ES(rand)
            ea.terminator = [inspyred.ec.terminators.evaluation_termination, 
                             inspyred.ec.terminators.diversity_termination]            
        else:
            ea.terminator = inspyred.ec.terminators.evaluation_termination                     

        #ea.observer = my_observer
        ea.observer = [inspyred.ec.observers.stats_observer, inspyred.ec.observers.file_observer]
        self.pop = ea.evolve(evaluator=datafit,generator=mygenerate,maximize=False,
                             pop_size=pop_size,max_evaluations=max_evaluations,num_elites=1,
                             bounder=inspyred.ec.Bounder(0.,1.),**kwargs)
        self.pop.sort(reverse=True)
        self.fits=[ind.fitness for ind in self.pop]

def test_VESManager(showProgress=False):
    """
        run from console with: python -c 'import pygimli.physics.ert.ves as pg; pg.test_VESManager(1)'
    """
    thicks = [2., 10.]
    res = [100., 5., 30]
    phi = [0., 20., 0.]

    # model fails
    thicks = [2., 6., 10.]
    res = [100., 500., 20., 800.]
    phi = [0., 20., 50., 0]

    synthModel = pg.cat(thicks, res)
    ab2 = np.logspace(np.log10(1.5), np.log10(100.), 25)

    mgr = VESManager(verbose=True, debug=False)
    mgr.fop.setRegionProperties(0, limits=[0.5, 200], trans='log')
    ra, err = mgr.simulate(synthModel, ab2=ab2, mn2=1.0, noiseLevel=0.01)
    mgr.exportData('synth.ves', ra, err)

    mgr.invert(ra, err, nLayer=4, lam=100,
               showProgress=showProgress)

    pg.wait()
    ### Test -- reinit with new parameter count
    mgr.invert(ra, err, nLayer=3,
               showProgress=showProgress)

    #np.testing.assert_array_less(mgr.inv.inv.chi2(), 1)

    ### Test -- reinit with new data basis
    ab2 = np.logspace(np.log10(1.5), np.log10(50.), 10)
    ra, err = mgr.simulate(synthModel, ab2=ab2, mn2=1.0, noiseLevel=0.01)

    mgr2 = VESManager(verbose=False, debug=False)
    mgr2.invert(ra, err, nLayer=3, ab2=ab2, mn2=1.0,
                showProgress=showProgress)

    #np.testing.assert_array_less(mgr2.inv.inv.chi2(), 1)

    pg.wait()
    ### Test -- reinit with complex resistivies
    mgr.complex = True
    synthModel =  pg.cat(synthModel, phi)

    ra, err = mgr.simulate(synthModel, ab2=ab2, mn2=1.0, noiseLevel=0.01)
    mgr.exportData('synthc.ves', ra, err)
    mgr.invert(ra, err,
               showProgress=showProgress)

    np.testing.assert_array_less(mgr.inv.inv.chi2(), 1)

    if showProgress:
        print("test done");
        pg.wait()

    def setData(self, data):
        self.fop.setData(data)
        self.dataVals = pg.Vector(0)
        self.dataErrs = pg.Vector(0)

        for i, mgr in enumerate(self.mgrs):
            self.tD.add(mgr.tD, data[i].size())
            self.dataVals = pg.cat(self.dataVals, mgr.dataVals(data[i]))
            self.dataErrs = pg.cat(self.dataErrs, mgr.relErrorVals(data[i]))
        self.inv.setTransData(self.tD)

 def response(self, par):
     """Response vector as combined in-phase and out-phase data."""
     thk = np.asarray(par[:self.nlay-1], dtype=np.float)
     res = np.asarray(par[self.nlay-1:2*self.nlay-1], dtype=np.float)
     mur = np.asarray(par[2*self.nlay-1:3*self.nlay-1], dtype=np.float) + 1
     ip, op = self.vmd_hem(self.height, rho=res, d=thk, mur=mur)
     return pg.cat(ip, op)

 def drawModel(self, ax, model):
     nLay = (len(model)+1) // 3
     super(VESCModelling, self).drawModel(ax, model[0:nLay*2-1])
     pg.mplviewer.drawModel1D(ax=ax,
                              model=pg.cat(model[0:nLay-1], model[nLay*2-1::]),
                              plot='plot',
                              xlabel='Phase [mrad]')

def drawSeismogramm(axes, mesh, u, ids, dt, i=None):
    r"""Extract and show time series from wave field

    Parameters
    ----------
    """
    axes.set_xlim(-20., 20.)
    axes.set_ylim(0., dt*len(u)*1000)
    axes.set_aspect(1)
    axes.set_ylabel('Time in ms')

    if i is None:
        i = len(u)-1

    t = np.linspace(0, i*dt*1000, i+1)

    for iw, n in enumerate(ids):
        pos = mesh.node(n).pos()
        print(pos)
        axes.plot(pos[0], 0.05, '^', color='black')

        trace = pg.cat(pg.RVector(0), u[:(i+1), n])
#        print(i+1, n)
#        print(trace, (max(pg.abs(trace))))

#        if max(pg.abs(trace)) > 1e-8:

        trace *= np.exp(0.5*t)
        trace /= (max(pg.abs(trace))*1.5)

        drawWiggle(axes, trace, t=t, xoffset=pos[0])
    axes.invert_yaxis()

    def response(self, model):
        """ yields forward model response """
        modA = np.asarray(model).reshape((self.nlay * 2 - 1, self.nx)).T
        resp = pg.RVector(0)
        for modi in modA:
            resp = pg.cat(resp, self.FOP1d.response(modi))

        return resp

    def response(self, par):
        ip, op = self.vmd_hem(self.height,
                              np.asarray(par)[self.nlay-1:self.nlay*2-1],
                              np.asarray(par)[:self.nlay-1])
#        ip, op = self.vmd_hem(self.height,
#                              np.asarray(par(self.nlay-1, self.nlay*2-1)),
#                              np.asarray(par(0, self.nlay-1)))
        return pg.cat(ip, op)

 def createDefaultStartModel(self):
     """
     """
     res = pb.getComplexData(self.data())
     parCount = self.regionManager().parameterCount()
     re = pg.RVector(parCount, pg.mean(pg.real(res)))
     im = pg.RVector(parCount, -pg.mean(pg.imag(res)))
     return pg.cat(re, im)

    def response(self, model):
        """
        """
        modA = np.asarray(model).reshape((self.nlay_*2-1,self.nx_)).T
        resp = pg.RVector(0)
        for modi in modA:
            resp = pg.cat(resp, self.FOP_.response(modi))

        return resp

    def getIntegrationWeights(self, rMin, rMax):
        """Retrieve Gauss-Legende/Laguerre integration weights."""
        nGauLegendre = max(int((6.0 * np.log10(rMax / rMin))), 4)
        nGauLaguerre = 4

        k = pg.RVector()
        w = pg.RVector()

        k0 = 1.0 / (2.0 * rMin)
        pg.GaussLegendre(0.0, 1.0, nGauLegendre, k, w)
        kLeg = k0 * k * k
        wLeg = 2.0 * k0 * k * w / np.pi

        pg.GaussLaguerre(nGauLaguerre, k, w)
        kLag = k0 * (k + 1.0)
        wLag = k0 * np.exp(k) * w / np.pi

        return pg.cat(kLeg, kLag), pg.cat(wLeg, wLag)

 def loadData(self, fileName, **kwargs):
     mat = np.loadtxt(fileName)
     if len(mat[0]) == 4:
         self.fop.setDataBasis(ab2=mat[:,0], mn2=mat[:,1])
         return mat.T
     if len(mat[0]) == 6:
         self.complex = True
         self.fop.setDataBasis(ab2=mat[:,0], mn2=mat[:,1])
         return mat[:,0], mat[:,1], np.array(pg.cat(mat[:,2], mat[:,4])), np.array(pg.cat(mat[:,3], mat[:,5]))

    def createInv(self,nlay,lam=100.,errVES=3,verbose=True):
        """ create marquardt type inversion instance with data transformation """
        self.createFOP(nlay)
        self.tMod = pg.RTransLog()
        self.tMRS = pg.RTrans()
        self.tVES = pg.RTransLog()
        self.transData = pg.RTransCumulative()
        self.transData.push_back( self.tMRS, len(self.data) )
        self.transData.push_back( self.tVES, len(self.rhoa) )
        data = pg.cat(self.data, self.rhoa)
        self.INV = pg.RInversion(data, self.f, self.transData, verbose)
        self.INV.setLambda(lam)
        self.INV.setMarquardtScheme(0.8)
        self.INV.stopAtChi1(False) # now in MarquardtScheme
        self.INV.setDeltaPhiAbortPercent(0.5)
#        self.INV.setMaxIter(1)
        error = pg.cat(self.error, self.rhoa*errVES/100.)
        self.INV.setAbsoluteError(error)