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)

 def genMods( individual ):
     """ generate MRS and VES models from unit vector """
     model = pg.asvector( individual ) * ( self.lUB - self.lLB ) + self.lLB
     if self.logpar:
         model = pg.exp( model )
     
     modMRS = model(0,nlay*3-1)
     modVES = pg.cat(model(0,nlay-1),model(nlay*3-1,nlay*4-1))
     return modMRS, modVES

    def fitDebyeModel(self, ePhi=0.001, lam=1e3, lamFactor=0.8,
                      mint=None, maxt=None, nt=None, new=True,
                      showFit=False, cType=1):
        """fit a (smooth) continuous Debye model (Debye decomposition)"""
        nf = len(self.f)
        if mint is None:
            mint = .1 / max(self.f)
        if maxt is None:
            maxt = .5 / min(self.f)
        if nt is None:
            nt = nf*2
        # %% discretize tau, setup DD and perform DD inversion
        self.tau = np.logspace(log10(mint), log10(maxt), nt)
        phi = self.phi
        tLin, tLog, tM = pg.RTrans(), pg.RTransLog(), pg.RTransLogLU(0., 1.)
        if new:
            reNorm, imNorm = self.zNorm()
            fDD = DebyeComplex(self.f, self.tau)
            Znorm = pg.cat(reNorm, imNorm)
            IDD = pg.RInversion(Znorm, fDD, tLog, tM, False)
            IDD.setAbsoluteError(max(Znorm)*0.003+0.01)
        else:
            fDD = DebyePhi(self.f, self.tau)
            IDD = pg.RInversion(phi, fDD, tLin, tM, True)
            IDD.setAbsoluteError(ePhi)  # 1 mrad

        fDD.regionManager().setConstraintType(cType)
        IDD.stopAtChi1(False)
        startModel = pg.RVector(nt, 0.01)
        IDD.setModel(startModel)
        IDD.setLambda(lam)
        IDD.setLambdaFactor(lamFactor)
        self.mDD = IDD.run()
        IDD.echoStatus()
        if new:
            resp = np.array(IDD.response())
            respRe = resp[:nf]
            respIm = resp[nf:]
            respC = ((1 - respRe) + respIm * 1j) * max(self.amp)
            self.phiDD = np.angle(respC)
            self.ampDD = np.abs(respC)
            if showFit:
                fig, ax = self.showData(znorm=True, nrows=3)
                ax[0].plot(self.f, respRe, 'r-')
                ax[1].plot(self.f, respIm, 'r-')
                ax[2].semilogx(self.tau, self.mDD, 'r-')
                ax[2].set_xlim(max(self.tau), min(self.tau))
                ax[2].set_ylim(0., max(self.mDD))
                ax[2].grid(True)
                ax[2].set_xlabel(r'$\tau$ [s]')
                ax[2].set_xlabel('$m$ [-]')
        else:
            self.phiDD = IDD.response()
            if showFit:
                fig, ax = self.showData(nrows=3)
                ax[2].semilogx(self.tau, self.mDD, 'r-')

    def inv2D(self, nlay, lam=100., resL=1., resU=1000., thkL=1.,
              thkU=100., minErr=1.0):
        """
            2d LCI inversion class
        """

        if isinstance(nlay, int):
            modVec = pg.RVector(nlay * 2 - 1, 30.)
            cType = 0  # no reference model
        else:
            modVec = nlay
            cType = 10  # use this as referencemodel
            nlay = (len(modVec) + 1) / 2

        # init forward operator
        self.f2d = self.FOP2d(nlay)

        # transformations
        self.tD = pg.RTrans()
        self.tThk = pg.RTransLogLU(thkL, thkU)
        self.tRes = pg.RTransLogLU(resL, resU)

        for i in range(nlay - 1):
            self.f2d.region(i).setTransModel(self.tThk)

        for i in range(nlay - 1, nlay * 2 - 1):
            self.f2d.region(i).setTransModel(self.tRes)

        # set constraints
        self.f2d.region(0).setConstraintType(cType)
        self.f2d.region(1).setConstraintType(cType)

        # collect data vector
        datvec = pg.RVector(0)

        for i in range(len(self.x)):
            datvec = pg.cat(datvec, self.datavec(i))

        # collect error vector
        if self.ERR is None:
            error = 1.0
        else:
            error = []
            for i in range(len(self.x)):
                err = np.maximum(self.ERR[i][self.activeFreq] * 0.701, minErr)
                error.extend(err)

        # generate starting model by repetition
        model = pg.asvector(np.repeat(modVec, len(self.x)))
        INV = pg.RInversion(datvec, self.f2d, self.tD)
        INV.setAbsoluteError(error)
        INV.setLambda(lam)
        INV.setModel(model)
        INV.setReferenceModel(model)

        return INV

    def blockLCInversion(self, nlay=2, startModel=None, **kwargs):
        """Laterally constrained (piece-wise 1D) block inversion."""
        data, error, self.nData = pg.RVector(), pg.RVector(), []
        for mrs in self.mrs:
            data = pg.cat(data, mrs.data)
            error = pg.cat(error, mrs.error)
            self.nData.append(len(mrs.data))

        fop = MRSLCI(self.mrs, nlay=nlay)
        fop.region(0).setZWeight(kwargs.pop('zWeight', 0))
        fop.region(0).setConstraintType(kwargs.pop('cType', 1))
        transData, transMod = pg.RTrans(), pg.RTransLog()  # LU(1., 500.)
        if startModel is None:
            startModel = self.block1dInversion(nlay, verbose=False)
        model = kwargs.pop('startvec', np.tile(startModel, len(self.mrs)))
        INV = pg.RInversion(data, fop, transData, transMod, True, False)
        INV.setModel(model)
        INV.setReferenceModel(model)
        INV.setAbsoluteError(error)
        INV.setLambda(kwargs.pop('lam', 100))
        INV.setMaxIter(kwargs.pop('maxIter', 20))
#        INV.stopAtChi1(False)
        INV.setLambdaFactor(0.9)
        INV.setDeltaPhiAbortPercent(0.1)
        model = INV.run()
        self.WMOD, self.TMOD = [], []
        for par in np.reshape(model, (len(self.mrs), 3*nlay-1)):
            thk = par[0:nlay-1]
            self.WMOD.append(np.hstack((thk, par[nlay-1:2*nlay-1])))
            self.TMOD.append(np.hstack((thk, par[2*nlay-1:3*nlay-1])))

        ind = np.hstack((0, np.cumsum(self.nData)))
        resp = INV.response()
        misfit = data - resp
        emisfit = misfit / error
        misfit *= 1e9
        self.totalChi2 = INV.chi2()
        self.totalRMS = INV.absrms()*1e9
        self.RMSvec, self.Chi2vec = [], []
        for i in range(len(self.mrs)):
            self.RMSvec.append(np.sqrt(np.mean(misfit[ind[i]:ind[i+1]]**2)))
            self.Chi2vec.append(np.mean(emisfit[ind[i]:ind[i+1]]**2))

    def block1dInversionNew(self, nlay=2, lam=100., 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])

        INV = pg.RInversion(data, f, verbose)
        INV.setLambda(lam)
        INV.setMarquardtScheme(0.8)
#        INV.stopAtChi1(False)  # should be already in MarquardtScheme
        INV.setDeltaPhiAbortPercent(0.5)
        INV.setAbsoluteError(error)
        model = INV.run()
        return model