from SimPEG import *



class FieldsDC_CC(Problem.Fields):
    knownFields = {'phi_sol':'CC'}
    aliasFields = {
                    'phi' : ['phi_sol','CC','_phi'],
                    'e' : ['phi_sol','F','_e'],
                    'j' : ['phi_sol','F','_j']
                  }

    def __init__(self,mesh,survey,**kwargs):
        super(FieldsDC_CC, self).__init__(mesh, survey, **kwargs)

    def startup(self):
        self._cellGrad = self.survey.prob.mesh.cellGrad
        self._Mfinv = self.survey.prob.mesh.getFaceInnerProduct(invMat=True)

    def _phi(self, phi_sol, srcList):
        phi = phi_sol
        # for i, src in enumerate(srcList):
        #     phi_p = src.phi_p(self.survey.prob)
        #     if phi_p is not None:
        #         phi[:,i] += phi_p     
        return phi

    def _e(self, phi_sol, srcList):
        e = -self._cellGrad*phi_sol
        # for i, src in enumerate(srcList):
        #     e_p = src.e_p(self.survey.prob)
        #     if e_p is not None:
        #         e[:,i] += e_p     
        return e
    
    def _j(self, phi_sol, srcList):

        j = -self._Mfinv*self.survey.prob.Msig*self._cellGrad*phi_sol
        # for i, src in enumerate(srcList):
        #     j_p = src.j_p(self.survey.prob)
        #     if j_p is not None:
        #         j[:,i] += j_p
        return j



class SrcDipole(Survey.BaseSrc):
    """A dipole source, locA and locB are moved to the closest cell-centers"""

    current = 1
    loc  = None
    # _rhsDict = None

    def __init__(self, rxList, locA, locB, **kwargs):
        self.loc = (locA, locB)
        super(SrcDipole, self).__init__(rxList, **kwargs)

    def eval(self, prob):
        # Recompute rhs 
        # if getattr(self, '_rhsDict', None) is None:
        #     self._rhsDict = {}
        # if mesh not in self._rhsDict:
        pts = [self.loc[0], self.loc[1]]
        inds = Utils.closestPoints(prob.mesh, pts)
        q = np.zeros(prob.mesh.nC)
        q[inds] = - self.current * ( np.r_[1., -1.] / prob.mesh.vol[inds] )
        # self._rhsDict[mesh] = q
        # return self._rhsDict[mesh]
        return q


class RxDipole(Survey.BaseRx):
    """A dipole source, locA and locB are moved to the closest cell-centers"""
    def __init__(self, locsM, locsN, **kwargs):
        locs = (locsM, locsN)
        assert locsM.shape == locsN.shape, 'locs must be the same shape.'
        super(RxDipole, self).__init__(locs, 'dipole', storeProjections=False, **kwargs)

    @property
    def nD(self):
        """Number of data in the receiver."""
        return self.locs[0].shape[0]

    def getP(self, mesh):
        P0 = mesh.getInterpolationMat(self.locs[0], self.projGLoc)
        P1 = mesh.getInterpolationMat(self.locs[1], self.projGLoc)
        return P0 - P1


class SurveyDC(Survey.BaseSurvey):
    """
        **SurveyDC**

        Geophysical DC resistivity data.

    """

    def __init__(self, srcList, **kwargs):
        self.srcList = srcList
        Survey.BaseSurvey.__init__(self, **kwargs)
        # self._rhsDict = {}
        self._Ps = {}

    def projectFields(self, u):
        """
            Predicted data.

            .. math::
                d_\\text{pred} = Pu(m)
        """
        P = self.getP(self.prob.mesh)
        return P*mkvc(u[self.srcList, 'phi_sol'])

    def getP(self, mesh):
        if mesh in self._Ps:
            return self._Ps[mesh]

        P_src = [sp.vstack([rx.getP(mesh) for rx in src.rxList]) for src in self.srcList]

        self._Ps[mesh] = sp.block_diag(P_src)
        return self._Ps[mesh]


class ProblemDC_CC(Problem.BaseProblem):
    """
        **ProblemDC**

        Geophysical DC resistivity problem.

    """

    surveyPair = SurveyDC
    Solver     = Solver
    fieldsPair = FieldsDC_CC
    Ainv = None

    def __init__(self, mesh, **kwargs):
        Problem.BaseProblem.__init__(self, mesh)
        self.mesh.setCellGradBC('neumann')
        Utils.setKwargs(self, **kwargs)


    deleteTheseOnModelUpdate = ['_A', '_Msig', '_dMdsig']

    @property
    def Msig(self):
        if getattr(self, '_Msig', None) is None:
            sigma = self.curModel.transform
            Av = self.mesh.aveF2CC
            self._Msig = Utils.sdiag(1/(self.mesh.dim * Av.T * (1/sigma)))
        return self._Msig

    @property
    def dMdsig(self):
        if getattr(self, '_dMdsig', None) is None:
            sigma = self.curModel.transform
            Av = self.mesh.aveF2CC
            dMdprop = self.mesh.dim * Utils.sdiag(self.Msig.diagonal()**2) * Av.T * Utils.sdiag(1./sigma**2)
            self._dMdsig = lambda Gu: Utils.sdiag(Gu) * dMdprop
        return self._dMdsig

    @property
    def A(self):
        """
            Makes the matrix A(m) for the DC resistivity problem.

            :param numpy.array m: model
            :rtype: scipy.csc_matrix
            :return: A(m)

            .. math::
                c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0

            Where M() is the mass matrix and mT is the model transform.
        """
        if getattr(self, '_A', None) is None:
            D = self.mesh.faceDiv
            G = self.mesh.cellGrad
            self._A = D*self.Msig*G
            # Remove the null space from the matrix.
            self._A[0,0] /= self.mesh.vol[0]
            self._A = self._A.tocsc()
        return self._A

    def getRHS(self):
        # if self.mesh not in self._rhsDict:
        RHS = np.array([src.eval(self) for src in self.survey.srcList]).T
        # self._rhsDict[mesh] = RHS
        # return self._rhsDict[mesh]
        return RHS

    def fields(self, m):

        F = self.fieldsPair(self.mesh, self.survey)        
        self.curModel = m
        A    = self.A
        self.Ainv = self.Solver(A, **self.solverOpts)
        RHS    = self.getRHS()
        Phi  = self.Ainv * RHS   
        Srcs = self.survey.srcList
        F[Srcs, 'phi_sol'] = Phi

        return F

    def Jvec(self, m, v, u=None):
        """
            :param numpy.array m: model
            :param numpy.array v: vector to multiply
            :param numpy.array u: fields
            :rtype: numpy.array
            :return: Jv

            .. math::
                c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0

                \\nabla_u (A(m)u - q) = A(m)

                \\nabla_m (A(m)u - q) = G\\text{sdiag}(Du)\\nabla_m(M(mT(m)))

            Where M() is the mass matrix and mT is the model transform.

            .. math::
                J = - P \left( \\nabla_u c(m, u) \\right)^{-1} \\nabla_m c(m, u)

                J(v) = - P ( A(m)^{-1} ( G\\text{sdiag}(Du)\\nabla_m(M(mT(m))) v ) )
        """
        # Set current model; clear dependent property $\mathbf{A(m)}$
        self.curModel = m
        sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
        if u is None:
            # Run forward simulation if $u$ not provided
            u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
        else:
            u = u[self.survey.srcList, 'phi_sol']

        D = self.mesh.faceDiv
        G = self.mesh.cellGrad
        # Derivative of model transform, $\deriv{\sigma}{\m}$
        dsigdm_x_v = self.curModel.transformDeriv * v

        # Take derivative of $C(m,u)$ w.r.t. $m$
        dCdm_x_v = np.empty_like(u)
        # loop over fields for each source
        for i in range(self.survey.nSrc):
            # Derivative of inner product, $\left(\mathbf{M}_{1/\sigma}^f\right)^{-1}$
            dAdsig         = D * self.dMdsig( G * u[:,i] )
            dCdm_x_v[:, i] = dAdsig *  dsigdm_x_v

        # Take derivative of $C(m,u)$ w.r.t. $u$
        dA_du = self.A
        # Solve for $\deriv{u}{m}$
        # dCdu_inv = self.Solver(dCdu, **self.solverOpts)
        if self.Ainv is None:
            self.Ainv = self.Solver(dA_du, **self.solverOpts)

        P        = self.survey.getP(self.mesh)
        Jv    = - P * mkvc( self.Ainv * dCdm_x_v )
        return Jv

    def Jtvec(self, m, v, u=None):

        self.curModel = m
        sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
        if u is None:
            # Run forward simulation if $u$ not provided
            u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
        else:
            u = u[self.survey.srcList, 'phi_sol']

        shp = u.shape
        P = self.survey.getP(self.mesh)
        PT_x_v = (P.T*v).reshape(shp, order='F')

        D = self.mesh.faceDiv
        G = self.mesh.cellGrad
        dA_du = self.A
        mT_dm = self.mapping.deriv(m)

        # We probably always need this due to the linesearch .. (?)
        self.Ainv = self.Solver(dA_du.T, **self.solverOpts)
        # if self.Ainv is None:
        #     self.Ainv = self.Solver(dCdu, **self.solverOpts)

        w = self.Ainv * PT_x_v

        Jtv = 0
        for i, ui in enumerate(u.T):  # loop over each column
            Jtv += self.dMdsig( G * ui ).T * ( D.T * w[:,i] )

        Jtv = - mT_dm.T * ( Jtv )
        return Jtv

def readUBC_DC2DModel(fileName):

    from SimPEG import np, mkvc
    """
        Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg

        Input:
        :param fileName, path to the UBC GIF 2D model file

        Output:
        :param SimPEG TensorMesh 2D object
        :return
        
        Created on Thu Nov 12 13:14:10 2015

        @author: dominiquef

    """
    
    # Open fileand skip header... assume that we know the mesh already

    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
    
    dim = np.array(obsfile[0].split(),dtype=float)
    
    temp = np.array(obsfile[1].split(),dtype=float)
    
    if len(temp) > 1:
        model = np.zeros(dim)
        
        for ii in range(len(obsfile)-1):
            mm = np.array(obsfile[ii+1].split(),dtype=float)
            model[:,ii] = mm
            
        model = model[:,::-1]
        
    else:
        
        if len(obsfile[1:])==1:
            mm = np.array(obsfile[1:].split(),dtype=float)
            
        else:
            mm = np.array(obsfile[1:],dtype=float)
            
        # Permute the second dimension to flip the order
        model = mm.reshape(dim[1],dim[0])
    
        model = model[::-1,:]
        model = np.transpose(model, (1, 0))
        
    model = mkvc(model)


    return model

def plot_pseudoSection(Tx,Rx,data,z0, stype):
    
    from SimPEG import np, mkvc
    from scipy.interpolate import griddata
    from matplotlib.colors import LogNorm
    import pylab as plt
    import re
    """
        Read list of 2D tx-rx location and plot a speudo-section of apparent
        resistivity.
        
        Assumes flat topo for now...
    
        Input:
        :param d2D, z0
        :switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole)
    
        Output:
        :figure scatter plot overlayed on image
        
        Created on Mon December 7th, 2015
    
        @author: dominiquef
    
    """
    #d2D = np.asarray(d2D)
    
    midl = []
    midz = []
    rho = []
    
    for ii in range(len(Tx)):
        # Get distances between each poles
        rC1P1 = np.abs(Tx[ii][0] - Rx[ii][:,0]) 
        rC2P1 = np.abs(Tx[ii][1] - Rx[ii][:,0])
        rC1P2 = np.abs(Tx[ii][1] - Rx[ii][:,1])
        rC2P2 = np.abs(Tx[ii][0] - Rx[ii][:,1])
        rP1P2 = np.abs(Rx[ii][:,1] - Rx[ii][:,0])    
    
        # Compute apparent resistivity
        if re.match(stype,'pdp'):
            rho = np.hstack([rho, data[ii] * 2*np.pi  * rC1P1 * ( rC1P1 + rP1P2 ) / rP1P2] )
            
        elif re.match(stype,'dpdp'):
            rho = np.hstack([rho, data[ii] * 2*np.pi / ( 1/rC1P1 - 1/rC2P1 - 1/rC1P2 + 1/rC2P2 ) ])
    
        Cmid = (Tx[ii][0] + Tx[ii][1])/2
        Pmid = (Rx[ii][:,0] + Rx[ii][:,1])/2
    
        midl = np.hstack([midl, ( Cmid + Pmid )/2 ])
        midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + z0 ])
    
   
    # Grid points
    grid_x, grid_z = np.mgrid[np.min(midl):np.max(midl), np.min(midz):np.max(midz)]
    grid_rho = griddata(np.c_[midl,midz], np.log10(abs(1/rho.T)), (grid_x, grid_z), method='linear')
    
    
    #plt.subplot(2,1,2)
    plt.imshow(grid_rho.T, extent = (np.min(midl),np.max(midl),np.min(midz),np.max(midz)), origin='lower', alpha=0.8)
    cbar = plt.colorbar(format = '%.2f',fraction=0.02)
    cmin,cmax = cbar.get_clim()
    ticks = np.linspace(cmin,cmax,3)
    cbar.set_ticks(ticks)
    
    # Plot apparent resistivity
    plt.scatter(midl,midz,s=50,c=np.log10(abs(1/rho.T)))

def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
    
    from SimPEG import np
    import re
    """
        Load in endpoints and survey specifications to generate Tx, Rx location
        stations.
        
        Assumes flat topo for now...
    
        Input:
        :param endl -> input endpoints [x1, y1, z1, x2, y2, z2]
        :object mesh -> SimPEG mesh object
        :switch stype -> "dpdp" (dipole-dipole) | "pdp" (pole-dipole) | 'gradient'
        : param a, n -> pole seperation, number of rx dipoles per tx
        
        Output:
        :param Tx, Rx -> List objects for each tx location
            Lines: P1x, P1y, P1z, P2x, P2y, P2z
        
        Created on Wed December 9th, 2015
    
        @author: dominiquef
    
    """
    def xy_2_r(x1,x2,y1,y2):
        r = np.sqrt( np.sum((x2 - x1)**2 + (y2 - y1)**2) )
        return r 
        
    ## Evenly distribute electrodes and put on surface
    # Mesure survey length and direction
    dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1])
    
    dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
    dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
       
    nstn = np.floor( dl_len / a )
    
    # Compute discrete pole location along line
    stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a
    stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a
    
    # Create line of P1 locations
    M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
    
    # Create line of P2 locations
    N = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
    
    ## Build list of Tx-Rx locations depending on survey type
    # Dipole-dipole: Moving tx with [a] spacing -> [AB a MN1 a MN2 ... a MNn]
    # Pole-dipole: Moving pole on one end -> [A a MN1 a MN2 ... MNn a B]
    Tx = []
    Rx = []
    
    if not re.match(stype,'gradient'):
        
        for ii in range(0, int(nstn)-1): 
            
            
            if re.match(stype,'dpdp'):
                tx = np.c_[M[ii,:],N[ii,:]]
            elif re.match(stype,'pdp'):
                tx = np.c_[M[ii,:],M[ii,:]]
                
            #Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]])
            
            # Current elctrode seperation
            AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
            
            # Number of receivers to fit
            nstn = np.min([np.floor( (AB - b) / a ) , n])
            
            # Check if there is enough space, else break the loop
            if nstn <= 0:
                continue
            
            # Compute discrete pole location along line
            stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a
            stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a
            
            # Create receiver poles
            # Create line of P1 locations
            P1 = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
            
            # Create line of P2 locations
            P2 = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
            
            Rx.append(np.c_[P1,P2])
            Tx.append(tx)            
            
#==============================================================================
#     elif re.match(stype,'dpdp'):
#         
#         for ii in range(0, int(nstn)-2):  
#             
#             indx = np.min([ii+n+1,nstn])
#             Tx.append(np.c_[M[ii,:],N[ii,:]])
#             Rx.append(np.c_[M[ii+2:indx,:],N[ii+2:indx,:]])
#==============================================================================
            
    elif re.match(stype,'gradient'):
        
        # Gradient survey only requires Tx at end of line and creates a square
        # grid of receivers at in the middle at a pre-set minimum distance
        Tx.append(np.c_[M[0,:],N[-1,:]])
              
        # Get the edge limit of survey area
        min_x = endl[0,0] + dl_x * b
        min_y = endl[0,1] + dl_y * b
            
        max_x = endl[1,0] - dl_x * b
        max_y = endl[1,1] - dl_y * b
        
        box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
        box_w = box_l/2.
        
        nstn = np.floor( box_l / a )
        
        # Compute discrete pole location along line
        stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
        stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
        
        # Define number of cross lines
        nlin = int(np.floor( box_w / a ))
        lind = range(-nlin,nlin+1) 
        
        ngrad = nstn * len(lind)
        
        rx = np.zeros([ngrad,6])
        for ii in range( len(lind) ):
            
            # Move line in perpendicular direction by dipole spacing
            lxx = stn_x - lind[ii]*a*dl_y
            lyy = stn_y + lind[ii]*a*dl_x
            
            
            M = np.c_[ lxx, lyy , np.ones(nstn).T*mesh.vectorNz[-1]]
            N = np.c_[ lxx+a*dl_x, lyy+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
            
            rx[(ii*nstn):((ii+1)*nstn),:] = np.c_[M,N]
            
        Rx.append(rx)
        
    else:
        print """stype must be either 'pdp', 'dpdp' or 'gradient'. """

        
   
    return Tx, Rx  

def writeUBC_DCobs(fileName,Tx,Rx,d,wd, dtype):
    
    from SimPEG import np, mkvc
    import re
    """
        Read UBC GIF DCIP 3D observation file and generate arrays for tx-rx location
    
        Input:
        :param fileName, path to the UBC GIF 3D obs file
    
        Output:
        :param rx, tx, d, wd
        :return
        
        Created on Mon December 7th, 2015
    
        @author: dominiquef
    
    """
    fid = open(fileName,'w')
    fid.write('! GENERAL FORMAT\n')  
    
    for ii in range(len(Tx)):
        
        tx = np.asarray(Tx[ii])
        rx = np.asarray(Rx[ii])
        nrx = rx.shape[0]
        
        fid.write('\n')
        
        if re.match(dtype,'2D'):
            
            for jj in range(nrx):
                
                fid.writelines("%e " % ii for ii in mkvc(tx))
                fid.writelines("%e " % ii for ii in mkvc(rx[jj]))
                fid.write('%e %e\n'% (d[ii][jj],wd[ii][jj]))
                #np.savetxt(fid, np.c_[ rx ,np.asarray(d[ii]), np.asarray(wd[ii]) ], fmt='%e',delimiter=' ',newline='\n') 
        
        elif re.match(dtype,'3D'):
                                
            fid.write('\n')
            fid.writelines("%e " % ii for ii in mkvc(tx))
            fid.write('%i\n'% nrx)
            np.savetxt(fid, np.c_[ rx ,np.asarray(d[ii]), np.asarray(wd[ii]) ], fmt='%e',delimiter=' ',newline='\n') 
        
        
    fid.close()  

def convertObs_DC3D_to_2D(Tx,Rx):
    
    from SimPEG import np
    import numpy.matlib as npm
    """
        Read list of 3D Tx Rx location and change coordinate system to distance
        along line assuming all data is acquired along line
        First transmitter pole is assumed to be at the origin

        Assumes flat topo for now...
    
        Input:
        :param Tx, Rx
        
        Output:
        :figure Tx2d, Rx2d
        
        Created on Mon December 7th, 2015
    
        @author: dominiquef
    
    """
    
                
    Tx2d = []
    Rx2d = []

    for ii in range(len(Tx)):
        
        if ii == 0:
            endp = Tx[0][0:2,0]
        
        nrx = Rx[ii].shape[0]
                  
        rP1 = np.sqrt( np.sum( ( endp - Tx[ii][0:2,0] )**2 , axis=0))
        rP2 = np.sqrt( np.sum( ( endp - Tx[ii][0:2,1] )**2 , axis=0))
        rC1 = np.sqrt( np.sum( ( npm.repmat(endp.T,nrx,1) - Rx[ii][:,0:2] )**2 , axis=1))
        rC2 = np.sqrt( np.sum( ( npm.repmat(endp.T,nrx,1) - Rx[ii][:,3:5] )**2 , axis=1))
        
        Tx2d.append( np.r_[rP1, rP2] )
        Rx2d.append( np.c_[rC1, rC2] )
            #np.savetxt(fid, data, fmt='%e',delimiter=' ',newline='\n')

    return Tx2d, Rx2d

def readUBC_DC3Dobs(fileName):
    
    from SimPEG import np
    """
        Read UBC GIF DCIP 3D observation file and generate arrays for tx-rx location
    
        Input:
        :param fileName, path to the UBC GIF 3D obs file
    
        Output:
        :param rx, tx, d, wd
        :return
        
        Created on Mon December 7th, 2015
    
        @author: dominiquef
    
    """
       
    # Load file
    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
    
    # Pre-allocate
    Tx = []
    Rx = []
    d = []
    wd = []
    
    # Countdown for number of obs/tx
    count = 0
    for ii in range(obsfile.shape[0]):
        
        if not obsfile[ii]:
            continue
        
        # First line is transmitter with number of receivers
        if count==0:
    
            temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T)
            count = int(temp[-1])
            temp = np.reshape(temp[0:-1],[2,3]).T
            
            Tx.append(temp)
            rx = []
            continue
        
        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ')
        
            
        rx.append(temp)          
        
        count = count -1        
        
        # Reach the end of  
        if count == 0:
            temp = np.asarray(rx)
            Rx.append(temp[:,0:6])
            
            # Check for data + uncertainties
            if temp.shape[1]==8:
                d.append(temp[:,6])
                wd.append(temp[:,7])
                
            # Check for data only    
            elif temp.shape[1]==7:
                d.append(temp[:,6])
            
    return Tx, Rx, d, wd

def readUBC_DC2DLoc(fileName):

    from SimPEG import np
    """
        Read UBC GIF 2D observation file and generate arrays for tx-rx location

        Input:
        :param fileName, path to the UBC GIF 2D model file

        Output:
        :param rx, tx
        :return
        
        Created on Thu Nov 12 13:14:10 2015

        @author: dominiquef

    """
    
    # Open fileand skip header... assume that we know the mesh already
#==============================================================================
#     fopen = open(fileName,'r')
#     lines = fopen.readlines()
#     fopen.close()
#==============================================================================

    # Load file
    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
    
    # Check first line and figure out if 2D or 3D file format
    line = np.array(obsfile[0].split(),dtype=float)   
    
    tx_A  = []
    tx_B  = []
    rx_M  = []
    rx_N  = []
    d   = []
    wd  = []
    
    for ii in range(obsfile.shape[0]):
        
        # If len==3, then simple format where tx-rx is listed on each line
        if len(line) == 4:
        
            temp = np.fromstring(obsfile[ii], dtype=float,sep=' ')
            tx_A = np.hstack((tx_A,temp[0]))
            tx_B = np.hstack((tx_B,temp[1]))
            rx_M = np.hstack((rx_M,temp[2]))
            rx_N = np.hstack((rx_N,temp[3]))
            
        
    rx = np.transpose(np.array((rx_M,rx_N)))    
    tx = np.transpose(np.array((tx_A,tx_B)))
    
    return tx, rx, d, wd

def readUBC_DC2DMesh(fileName):

    from SimPEG import np
    """
        Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg

        Input:
        :param fileName, path to the UBC GIF mesh file

        Output:
        :param SimPEG TensorMesh 2D object
        :return
        
        Created on Thu Nov 12 13:14:10 2015

        @author: dominiquef

    """
    
    # Open file 
    fopen = open(fileName,'r')
    
    # Read down the file and unpack dx vector
    def unpackdx(fid,nrows):
        for ii in range(nrows):
            
            line = fid.readline()
            var = np.array(line.split(),dtype=float)
            
            if ii==0:
                x0= var[0]
                xvec = np.ones(int(var[2])) * (var[1] - var[0]) / int(var[2])
                xend = var[1]            
                
            else:
                xvec = np.hstack((xvec,np.ones(int(var[1])) * (var[0] - xend) / int(var[1])))
                xend = var[0] 
                
        return x0, xvec
        
    #%% Start with dx block
    # First line specifies the number of rows for x-cells
    line = fopen.readline()
    nl = np.array(line.split(),dtype=float)

    [x0, dx] = unpackdx(fopen,nl) 
    

    #%% Move down the file until reaching the z-block
    line = fopen.readline()
    if not line:
        line = fopen.readline()
        
    #%% End with dz block
    # First line specifies the number of rows for z-cells
    line = fopen.readline()
    nl = np.array(line.split(),dtype=float)

    [z0, dz] = unpackdx(fopen,nl) 
    
    # Flip z0 to be the bottom of the mesh for SimPEG
    z0 = z0 - sum(dz)
    dz = dz[::-1]
    #%% Make the mesh using SimPEG
    
    from SimPEG import Mesh
    tensMsh = Mesh.TensorMesh([dx,dz],(x0, z0))
    return tensMsh

