def __init__(self, which_case, LUT, RandomSamples, interp_type):
     print 'SciPy Interpolating ', which_case
     
     select = {\
     "rhoe":('Density','StaticEnergy'),\
     "PT":('Pressure','Temperature'),\
     "Prho":('Pressure','Density'),\
     "rhoT":('Density','Temperature'),\
     "Ps":('Pressure','Entropy'),\
     "hs":('Enthalpy','Entropy')\
     }
     
     thermo1, thermo2, = select[which_case]
     x =getattr(LUT,thermo1)
     y =getattr(LUT,thermo2)
     samples_x = getattr(RandomSamples,thermo1)
     samples_y = getattr(RandomSamples,thermo2)
     setattr(self,thermo1, samples_x)
     setattr(self,thermo2, samples_y)
     
     variables = sp.array(['Temperature','Density','Enthalpy','StaticEnergy',\
     'Entropy','Pressure','SoundSpeed2','dPdrho_e','dPde_rho',\
     'dTdrho_e','dTde_rho','Cp','Mu','Kt']);
     
     for var in variables[sp.where((variables!=thermo1) * (variables!=thermo2))]:
         z = getattr(LUT,var)            
         interp_func = sp.interpolate.griddata((x,y),z,sp.column_stack((samples_x,samples_y)),\
         method=interp_type) 
         nan_index = sp.where(sp.isnan(interp_func))
         interp_func[nan_index]= sp.interpolate.griddata((x,y),z,\
         sp.column_stack((samples_x[nan_index],samples_y[nan_index])),\
         method='nearest') 
         setattr(self,var,interp_func)
         
     return  

    def _packData(self, G1, indices2select, effect):
        if effect == 'fixed':
            if G1 is None and self.G0 is None:
                data = self.X[self.data_permutation][indices2select]
            elif G1 is None:
                data = sp.column_stack((self.G0[self.data_permutation][indices2select],
                                        self.X[self.data_permutation][indices2select]))
            elif self.G0 is None:
                data = sp.column_stack((G1[self.data_permutation][indices2select],
                                        self.X[self.data_permutation][indices2select]))
            else:
                data = sp.column_stack((self.G0[self.data_permutation][indices2select],
                                        G1[self.data_permutation][indices2select],
                                        self.X[self.data_permutation][indices2select]))
        elif effect == 'mixed':
            X = self.X[self.data_permutation]
            if self.G0 is not None:
                G0 = self.G0[self.data_permutation]
            if G1 is not None:
                G1 = G1[self.data_permutation]

            data = []
            for i in range(len(indices2select)):

                lis = [X[indices2select[i]]]
                if G0 is not None:
                    lis.append( G0[indices2select[i]] )
                if G1 is not None:
                    lis.append( G1[indices2select[i]] )

                data.append( lis )
            else:
                assert False, 'Unkown effect type.'

        return (data, self.Y[self.data_permutation][indices2select])

def GetFermiSigns(filename,refstate=None,channel=None):
    attr=GetAttr(filename)
    filetype=''
    try:
        filetype=attr['type']
    except KeyError:
        mo=re.match('.*/?[0-9]+-(.*)\.h5',filename)
        filetype=mo.groups()[0]
    L=attr['L']
    if refstate==None:
        refstate=sc.zeros(2*L*L)
        refstate[0::2]=1
    if filetype=='WaveFunction':
        hfile=h5py.File(filename,'r')
        if 'states_up' in hfile.keys():
            states=sc.column_stack([hfile['states_up'],hfile['states_do']])
        else:
            states=sc.column_stack([hfile['states_0'],hfile['states_1']])
        hfile.close()
        return sf.fermisigns(states,refstate)
    else:
        if channel==None:
            channel=attr['channel']
        L=int(attr['L'])
        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]
        q=[float(attr['qx'])/L,float(attr['qy'])/L]
        if channel=='trans':
            return sf.transfermisigns(L,L,shift,q)
        elif channel=='long':
            return sf.longfermisigns(L,L,shift,q)
        else:
            raise KeyError('\"channel\" must be either \"trans\" or \"long\".')

def writeCSVOutput(hdf5File=None,cout=None):
    f = h5py.File(hdf5File,'r')

    csv_filename = os.path.join(cout,f['phenotype_name'].value.replace(" ","_") + ".csv") 
    csv_header = None
    csv_matrix = None
    
    if "betas" in f.keys():
        csv_header = sp.array(["CHR","Positions","P-Value","TestStatistic","Q-Value","Benjamini-Hochberg-P-Value","Benjamini-Hochberg-Yekutieli-P-Value","Beta0","SEBeta0","Beta1","SEBeta1","MAF","SNP-Hash"])
    else:
        csv_header = sp.array(["CHR","Positions","P-Value","TestStatistic","Q-Value","Benjamini-Hochberg-P-Value","Benjamini-Hochberg-Yekutieli-P-Value","MAF","SNP-Hash"])
    tmp_matrix = None
    if "betas" in f.keys():
        tmp_matrix = sp.column_stack([sp.array(f["chromosomes"],dtype="S50"),
                                          f['positions'],
                                          f['p_values'],
                                          f['scores'],
                                          f['q_values'],
                                          f['bh_p_values'],
                                          f['bhy_p_values'],
                                          f['betas'][:,0],
                                          f['betas_se'][:,0],
                                          f['betas'][:,1],
                                          f['betas_se'][:,1],
                                          f['maf'],
                                          f['snp_hash']])
    else:
        tmp_matrix = sp.column_stack([sp.array(f["chromosomes"],dtype="S50"),
                              f['positions'],
                              f['p_values'],
                              f['scores'],
                              f['q_values'],
                              f['bh_p_values'],
                              f['bhy_p_values'],
                              f['maf'],
                              f['snp_hash']])
     
    csv_matrix = tmp_matrix
    mf = open(csv_filename,'w')
    string = ""
    for i in xrange(csv_header.shape[0]):
        string += str(csv_header[i]) + ","
    mf.write(string[:-1] + "\n")
    for i in xrange(csv_matrix.shape[0]):
        string = ""
        for j in xrange(csv_matrix.shape[1]):
            string += str(csv_matrix[i,j]) + ","
        mf.write(string[:-1] + "\n")
    mf.close()
    f.close()

 def bm29write(self, filename= None,  newdata=None, comment=True, newdatacol_head=None, columns=None):
     """function to write a bm29 file, define a series of array with the same name
     of column defined in the file
     input \n
     filename =>string
     newdata =>numpy array
     comment =>boolean or new comments line
     col_head =>string heders of column
     columns =>list of strings with the columns to write 
     P.S or newdata and col_head or columns
     """
     if not(filename):
         try:
             filename =getattr(self,"fullfilename")
         except:
             raise FileFormatError('no filename specified and self.fullfilenames not defined')
     if os.path.exists(filename):
             filename += ".1"
     outFile = open(filename, 'w')
     
     if comment is True:
         if hasattr(self, "comments"):
             outFile.writelines(self.comments[:-2])
     elif comment is False:
           pass
     else: outFile.writelines(comment)
     if newdata is None:
         if columns:
             outFile.writelines("#N "+str(len(columns))+ "\n")   
             col_head= "#L "+" ".join(columns)
             outFile.writelines(col_head+"\n")
             f=lambda x: getattr(self,x)
             newdata= scipy.column_stack(map(f,columns))
         else:
             if self.All_Column == False:
                 outFile.writelines("#N 2\n")                    
                 outFile.writelines("#L   E     Mu\n")
                 newdata= scipy.column_stack((self.E, self.Mu))
             elif self.All_Column == True: 
                 try: 
                     outFile.writelines(self.comments[-2])                        
                     outFile.write("#L  "+" ".join(getattr(self, "col_head")))
                     outFile.write("\n")
                     newdata = self.data
                 except:
                     pass    
     scipy.savetxt(outFile, newdata, fmt= '%1.10f')
     outFile.close
     return

def readMahalih5(filename,des_site):
    """ This function will read the mahali GPS data into a GeoData data structure.
        The user only has to give a filename and name of the desired site.
        Input
            filename - A string that holds the file name.
            des_site - The site name. Should be listed in the h5 file in the
                table sites.
    """
    h5fn = Path(filename).expanduser()

    with h5py.File(str(h5fn), "r", libver='latest') as f:

        despnts = sp.where(f['data']['site']==des_site)[0]
        # TODO: hard coded for now
        doy =  doy= f['data']['time'][despnts]
        year = 2015*sp.ones_like(doy,dtype=int)

        TEC = f['data']['los_tec'][despnts]

        nTEC = f['data']['err_los_tec'][despnts]

        vTEC = f['data']['vtec'][despnts]
        az2sat = f['data']['az'][despnts]
        el2sat = f['data']['az'][despnts]

        piercelat = f['data']['pplat'][despnts]
        piercelong = f['data']['pplon'][despnts]
        satnum= f['data']['prn'][despnts]
        recBias = f['data']['rec_bias'][despnts]
        nrecBias = f['data']['err_rec_bias'][despnts]

    # Make the integration time on the order of 15 seconds.
    if (year==year[1]).all():
        unixyear =(datetime(year[0],1,1,0,0,0,tzinfo=UTC) - EPOCH).total_seconds()
        uttime = unixyear + sp.round_(24*3600*sp.column_stack((doy,doy+15./24./3600.))) # Making the difference in time to be a minute
    else:
        (y_u,y_iv) = np.unique(year,return_inverse=True)
        unixyearu = sp.array([(datetime(iy,1,1,0,0,0,tzinfo=UTC) - EPOCH).total_seconds() for iy in y_u])
        unixyear = unixyearu[y_iv]
        uttime = unixyear + 24*3600*sp.column_stack((doy,doy+15./24./3600.))


    data = {'TEC':TEC,'nTEC':nTEC,'vTEC':vTEC,'recBias':recBias,'nrecBias':nrecBias,'satnum':satnum,'az2sat':az2sat,'el2sat':el2sat}
    coordnames = 'WGS84'
    sensorloc = sp.nan*sp.ones(3)
    dataloc = sp.column_stack((piercelat,piercelong, 350e3*sp.ones_like(piercelat)))

    return (data,coordnames,dataloc,sensorloc,uttime)

    def kmeanspp_initialisation( self, X ):
        """Initialise means using K-Means++"""
        N, _ = X.shape
        k, d = self.k, self.d
        M = []

        # Choose one center amongst the X at random
        m = sc.random.randint( N )
        M.append( X[m] )

        # Choose k centers
        while( len( M ) < self.k ):
            # Create a probability distribution D^2 from the previous mean
            D = cdist( X, M ).min( 1 )**2
            assert( D.shape == (N,) )

            # Normalise and sample a new point
            D /= D.sum()

            m = sc.random.multinomial( 1, D ).argmax()
            M.append( X[m] )

        M = sc.column_stack( M )
        sigma = sc.sqrt(cdist( X, M.T, 'sqeuclidean').sum(0)/(N))
        w = ones( k )/float(k)

        return M, sigma, w

def main():
    saved_handler = sp.seterrcall(err_handler)
    saved_err = sp.seterr(all='call')

    print('============ Part 1: Plotting =============================')
    x, y = load_data('ex2/ex2data1.txt')
    plot_data(x, y)
    pl.show()

    print('============ Part 2: Compute Cost and Gradient ============')
    m, n = x.shape
    x = sp.column_stack((sp.ones((m, 1)), x))
    init_theta = sp.asmatrix(sp.zeros((n + 1, 1)))
    cost, grad = cost_function(init_theta, x, y)
    print('Cost at initial theta: %s' % cost)
    print('Gradient at initial theta:\n %s' % grad)

    print('============ Part 3: Optimizing minimize ====================')
    # res = op.minimize(cost_function, init_theta, args=(x, y), jac=True, method='Newton-CG')
    res = op.minimize(cost_function_without_grad, init_theta, args=(x, y), method='Powell')
    # print('Cost at theta found by fmin: %s' % cost)
    print('Result by minimize:\n%s' % res)
    plot_decision_boundary(res.x, x, y)
    pl.show()

    print('============ Part 4: Optimizing fmin ====================')
    res = op.fmin(cost_function_without_grad, init_theta, args=(x, y))
    # print('Cost at theta found by fmin: %s' % cost)
    print('Result by fmin:\n%s' % res)
    plot_decision_boundary(res, x, y)
    pl.show()

    sp.seterrcall(saved_handler)
    sp.seterr(**saved_err)

 def apply_flow(self,flowrate):
     r'''
     Convert the invaded sequence into an invaded time for a given flow rate
     considering the volume of invaded pores and throats.
     
     Parameters
     ----------
     flowrate : float
         The flow rate of the injected fluid
         
     Returns
     -------
     Creates a throat array called 'invasion_time' in the Algorithm 
     dictionary
     
     '''
     P12 = self._net['throat.conns']  # List of throats conns
     a = self['throat.invasion_sequence']  # Invasion sequence
     b = sp.argsort(self['throat.invasion_sequence'])
     P12_inv = self['pore.invasion_sequence'][P12]  # Pore invasion sequence
     # Find if the connected pores were invaded with or before each throat
     P1_inv = P12_inv[:,0] == a
     P2_inv = P12_inv[:,1] == a
     c = sp.column_stack((P1_inv,P2_inv))  
     d = sp.sum(c,axis=1,dtype=bool)  # List of Pores invaded with each throat
     # Find volume of these pores
     P12_vol = sp.zeros((self.Nt,))
     P12_vol[d] = self._net['pore.volume'][P12[c]]
     # Add invaded throat volume to pore volume (if invaded)
     T_vol = P12_vol + self._net['throat.volume']
     # Cumulative sum on the sorted throats gives cumulated inject volume
     e = sp.cumsum(T_vol[b]/flowrate)
     t = sp.zeros((self.Nt,))
     t[b] = e  # Convert back to original order
     self._phase['throat.invasion_time'] = t

 def run(self):
 
     # Parameters passed are current data array, along with time step
     # between current data points
     self.times = sp.arange(0,self.Tfinal,self.dt)
     self.sim = odeint(self.eqns,self.init,self.times,(self.inj,self.injdt))
     sp.savetxt('simulation.txt',sp.column_stack((self.times,self.sim)))

    def backsolve(self, T=None, vterm=None):
        """Solve finite system by backward recursion

        Parameters
        -------------
        T : int, optional
            Number of periods of time.

        Returns
        ----------
        X : array, shape (n, T)
            Optimal controls. An optimal policy for each starting state
        V : array, shape (n, T + 1)
            Value function.             

        """
        if T is None:
            if self.T is not None:
                T = self.T
            else:
                print ("Not a finite time model")
                return
        if vterm is None and self.vterm is None:
            vterm = sp.zeros(self.n)
        else:
            vterm = self.vterm
        x = sp.zeros((self.n, T), dtype=int)
        v = sp.column_stack((sp.zeros((self.n, T)), vterm))
        pstar = sp.zeros((self.n, self.n, T))
        for t in sp.arange(T - 1, -1, -1):
            v[ :, t] , x[ :, t]  = self.valmax(v[ : , t + 1])
            pstar[..., t] = self.valpol(x[:, t])[0]
        return (x, v, pstar)

def recover_topics( P, T, k, a0 ):
    """Recover the k components given input Pairs and Triples and
    $\\alpha_0$"""

    # Consider the k rank approximation of P,
    P = approxk( P, k )
    # Get the whitening matrix and coloring matrices
    W, Wt = get_whitener( P, k )

    # Whiten the third moment
    Tw = lambda theta: W.T.dot( T( W.dot(theta) ) ).dot( W )

    # Project Tw onto a matrix
    theta = orthogonal( k ).T[0] 

    U, S, _ = svd( Tw( theta ) )
    assert( (S > 1e-10).all() ) # Make sure it is non-singular

    O = []
    for i in xrange( k ):
        v = U.T[i]
        Zinv = (a0 + 2)/2 * (v.T.dot(Tw(v)).dot(v))
        O.append( Zinv * Wt.T.dot( v ) )

    O = sc.column_stack( O )

    return abs( O )

def makeinputh5(Iono,basedir):
    """This will make a h5 file for the IonoContainer that can be used as starting
    points for the fitter. The ionocontainer taken will be average over the x and y dimensions
    of space to make an average value of the parameters for each altitude.
    Inputs
    Iono - An instance of the Ionocontainer class that will be averaged over so it can
    be used for fitter starting points.
    basdir - A string that holds the directory that the file will be saved to.
    """
    # Get the parameters from the original data
    Param_List = Iono.Param_List
    dataloc = Iono.Cart_Coords
    times = Iono.Time_Vector
    velocity = Iono.Velocity
    zlist,idx = sp.unique(dataloc[:,2],return_inverse=True)
    siz = list(Param_List.shape[1:])
    vsiz = list(velocity.shape[1:])

    datalocsave = sp.column_stack((sp.zeros_like(zlist),sp.zeros_like(zlist),zlist))
    outdata = sp.zeros([len(zlist)]+siz)
    outvel = sp.zeros([len(zlist)]+vsiz)
    #  Do the averaging across space
    for izn,iz in enumerate(zlist):
        arr = sp.argwhere(idx==izn)
        outdata[izn] = sp.mean(Param_List[arr],axis=0)
        outvel[izn] = sp.mean(velocity[arr],axis=0)

    Ionoout = IonoContainer(datalocsave,outdata,times,Iono.Sensor_loc,ver=0,
                            paramnames=Iono.Param_Names, species=Iono.Species,velocity=outvel)
    Ionoout.saveh5(basedir/'startdata.h5')

    def train(self):
        if self.__algo_model=="MWUrt" or self.__algo_model=="WCrt":
            data = asso.MatrixData(x=self.__x,y=self.__y,covariates=self.__cov)
            self.__ass.setData(data)
            self.__ass.train()
        else:
            self.__ass.setPhenotype(self.__y)
            self.__ass.setGenotype(self.__x)

            if not self.__cov is None:
                self.__ass.setCovariates(sp.column_stack([sp.ones(self.__y.shape),self.__cov]))
            if self.__permutation==True:
                if self.__x.shape[1]<1000:
                    self.__perms = 1000000
                elif self.__x.shape[1]<10000:
                    self.__perms = 100000
                elif self.__x.shape[1]<100000:
                    self.__perms = 10000
                elif self.__x.shape[1]<1000000:
                    self.__perms = 1000
                elif self.__x.shape[1]<10000000:
                    self.__perms = 100
                else:
                    self.__perms = 10
                print "SNPS: ", self.__x.shape[1]
                print "Perms: ", self.__perms
                self.__ass.permutations(self.__perms)
            else:
                self.__ass.test_associations()

def makeinputh5(Iono,basedir):
    basedir = Path(basedir).expanduser()

    Param_List = Iono.Param_List
    dataloc = Iono.Cart_Coords
    times = Iono.Time_Vector
    velocity = Iono.Velocity
    zlist,idx = sp.unique(dataloc[:,2],return_inverse=True)
    siz = list(Param_List.shape[1:])
    vsiz = list(velocity.shape[1:])

    datalocsave = sp.column_stack((sp.zeros_like(zlist),sp.zeros_like(zlist),zlist))
    outdata = sp.zeros([len(zlist)]+siz)
    outvel = sp.zeros([len(zlist)]+vsiz)

    for izn,iz in enumerate(zlist):
        arr = sp.argwhere(idx==izn)
        outdata[izn]=sp.mean(Param_List[arr],axis=0)
        outvel[izn]=sp.mean(velocity[arr],axis=0)

    Ionoout = IonoContainer(datalocsave,outdata,times,Iono.Sensor_loc,ver=0,
                            paramnames=Iono.Param_Names, species=Iono.Species,velocity=outvel)


    ofn = basedir/'startdata.h5'
    print('writing {}'.format(ofn))
    Ionoout.saveh5(str(ofn))

def analysisdump(maindir, configfile, suptitle=None):
    """ This function will perform all of the plotting functions in this module
        given the main directory that all of the files live.
        Inputs
            maindir - The directory for the simulation.
            configfile - The name of the configuration file used.
            suptitle - The supertitle used on the files.
    """
    maindir = Path(maindir)
    plotdir = maindir.joinpath("AnalysisPlots")
    if not plotdir.is_dir():
        plotdir.mkdir()

    # plot spectrums
    filetemplate1 = str(maindir.joinpath("AnalysisPlots", "Spec"))
    filetemplate3 = str(maindir.joinpath("AnalysisPlots", "ACF"))
    filetemplate4 = str(maindir.joinpath("AnalysisPlots", "AltvTime"))
    (sensdict, simparams) = readconfigfile(configfile)
    angles = simparams["angles"]
    ang_data = sp.array([[iout[0], iout[1]] for iout in angles])
    if not sensdict["Name"].lower() in ["risr", "pfisr"]:
        ang_data_temp = ang_data.copy()
        beamlistlist = sp.array(simparams["outangles"]).astype(int)
        ang_data = sp.array([ang_data_temp[i].mean(axis=0) for i in beamlistlist])

    zenang = ang_data[sp.argmax(ang_data[:, 1])]
    rnggates = simparams["Rangegatesfinal"]
    rngchoices = sp.linspace(sp.amin(rnggates), sp.amax(rnggates), 4)
    angtile = sp.tile(zenang, (len(rngchoices), 1))
    coords = sp.column_stack((sp.transpose(rngchoices), angtile))
    times = simparams["Timevec"]

    filetemplate2 = str(maindir.joinpath("AnalysisPlots", "Params"))
    if simparams["Pulsetype"].lower() == "barker":
        params = ["Ne"]
        if suptitle is None:
            plotbeamparametersv2(times, configfile, maindir, params=params, filetemplate=filetemplate2, werrors=True)
        else:
            plotbeamparametersv2(
                times, configfile, maindir, params=params, filetemplate=filetemplate2, suptitle=suptitle, werrors=True
            )
    else:
        params = ["Ne", "Nepow", "Te", "Ti", "Vi"]
        if suptitle is None:
            plotspecs(coords, times, configfile, maindir, cartcoordsys=False, filetemplate=filetemplate1)
            plotacfs(coords, times, configfile, maindir, cartcoordsys=False, filetemplate=filetemplate3)
            plotbeamparametersv2(times, configfile, maindir, params=params, filetemplate=filetemplate2, werrors=True)
            beamvstime(configfile, maindir, params=params, filetemplate=filetemplate4)
        else:
            plotspecs(
                coords, times, configfile, maindir, cartcoordsys=False, filetemplate=filetemplate1, suptitle=suptitle
            )
            plotacfs(
                coords, times, configfile, maindir, cartcoordsys=False, filetemplate=filetemplate3, suptitle=suptitle
            )
            plotbeamparametersv2(
                times, configfile, maindir, params=params, filetemplate=filetemplate2, suptitle=suptitle, werrors=True
            )
            beamvstime(configfile, maindir, params=params, filetemplate=filetemplate4, suptitle=suptitle)

def SRIparams2iono(filename):

    fullfile = h5file(filename)
    fullfiledict = fullfile.readWholeh5file()

    #Size = Nrecords x Nbeams x Nranges x Nions+1 x 4 (fraction, temperature, collision frequency, LOS speed)
    fits = fullfiledict['/FittedParams']['Fits']
    (nt,nbeams,nrng,nspecs,nstuff) = fits.shape
    nlocs = nbeams*nrng
    fits = fits.transpose((1,2,0,3,4))
    fits = fits.reshape((nlocs,nt,nspecs,nstuff))
    #  Nrecords x Nbeams x Nranges
    Ne = fullfiledict['/FittedParams']['Ne']
    Ne = Ne.transpose((1,2,0))
    Ne = Ne.reshape((nlocs,nt))
    param_lists =sp.zeros((nlocs,nt,nspecs,2))
    param_lists[:,:,:,0] = fits[:,:,:,0]
    param_lists[:,:,:,1] = fits[:,:,:,1]
    param_lists[:,:,-1,0]=Ne
    Velocity = fits[:,:,0,3]


    if fullfiledict['/FittedParams']['IonMass']==16:
        species = ['O+','e-']
        pnames = sp.array([['Ni','Ti'],['Ne','Te']])

    time= fullfiledict['/Time']['UnixTime']
    time = time
    rng = fullfiledict['/FittedParams']['Range']
    bco = fullfiledict['/']['BeamCodes']
    angles = bco[:,1:3]
    (nang,nrg) = rng.shape

    allang = sp.tile(angles[:,sp.newaxis],(1,nrg,1))
    all_loc = sp.column_stack((rng.flatten(),allang.reshape(nang*nrg,2)))
    lkeep = ~ sp.any(sp.isnan(all_loc),1)
    all_loc = all_loc[lkeep]
    Velocity = Velocity[lkeep]
    param_lists = param_lists[lkeep]
    all_loc[:,0]=all_loc[:,0]*1e-3
    iono1 = IonoContainer(all_loc,param_lists,times=time,ver = 1,coordvecs = ['r','theta','phi'],
                          paramnames = pnames,species=species,velocity=Velocity)
                          
                          
                          
    # MSIS
    tn = fullfiledict['/MSIS']['Tn']
    tn = tn.transpose((1,2,0))
    tn = tn.reshape((nlocs,nt))
    
    
    startparams = sp.ones((nlocs,nt,2,2))
    startparams[:,:,0,1] = tn
    startparams[:,:,1,1] = tn
    startparams = startparams[lkeep]
    ionoS = IonoContainer(all_loc,startparams,times=time,ver = 1,coordvecs = ['r','theta','phi'],
                          paramnames = pnames,species=species)
    return iono1,ionoS

def planeFromPoints(points):
    """ Produces a plane from N points using least squares, and returns of the form:
        Z = aX + bY + c
        Follows logic from http://www.velocityreviews.com/forums/t368189-re-linear-regression-in-3-dimensions.html
    """
    x, y, z = zip(*points)
    A = column_stack([x, y, ones_like(x)])
    abc, residuals, rank, s = lstsq(A,z)
    return abc

def count_feature(X, tbl_lst=None, min_cnt=1):
    X_lst = [pd.Series(X[:, i]) for i in range(X.shape[1])]
    if tbl_lst is None:
        tbl_lst = [x.value_counts() for x in X_lst]
        if min_cnt > 1:
            tbl_lst = [s[s >= min_cnt] for s in tbl_lst]
    X = sp.column_stack([x.map(tbl).values for x, tbl in zip(X_lst, tbl_lst)])
    # NA(unseen values) to 0
    return np.nan_to_num(X), tbl_lst

    def __init__(self,ionoin,configfile):
        r2d = 180.0/sp.pi
        d2r = sp.pi/180.0
        (sensdict,simparams) = readconfigfile(configfile)
        nt = ionoin.Time_Vector.shape[0]
        nloc = ionoin.Sphere_Coords.shape[0]

        #Input location
        self.Cart_Coords_in = ionoin.Cart_Coords
        self.Sphere_Coords_In = ionoin.Sphere_Coords,
        self.Time_In = ionoin.Time_Vector
        self.Cart_Coords_In_Rep = sp.tile(ionoin.Cart_Coords,(nt,1))
        self.Sphere_Coords_In_Rep = sp.tile(ionoin.Sphere_Coords,(nt,1))
        self.Time_In_Rep  = sp.repeat(ionoin.Time_Vector,nloc,axis=0)

        #output locations
        rng_vec2 = simparams['Rangegatesfinal']
        nrgout = len(rng_vec2)

        angles = simparams['angles']
        nang =len(angles)

        ang_data = sp.array([[iout[0],iout[1]] for iout in angles])
        rng_all = sp.tile(rng_vec2,(nang))
        ang_all = sp.repeat(ang_data,nrgout,axis=0)
        nlocout = nang*nrgout

        ntout = len(simparams['Timevec'])
        self.Sphere_Coords_Out = sp.column_stack((rng_all,ang_all))
        (R_vec,Az_vec,El_vec) = (self.Sphere_Coords_Out[:,0],self.Sphere_Coords_Out[:,1],self.Sphere_Coords_Out[:,2])
        xvecmult = sp.cos(Az_vec*d2r)*sp.cos(El_vec*d2r)
        yvecmult = sp.sin(Az_vec*d2r)*sp.cos(El_vec*d2r)
        zvecmult = sp.sin(El_vec*d2r)
        X_vec = R_vec*xvecmult
        Y_vec = R_vec*yvecmult
        Z_vec = R_vec*zvecmult

        self.Cart_Coords_Out = sp.column_stack((X_vec,Y_vec,Z_vec))
        self.Time_Out = simparams['Timevec']
        self.Time_Out_Rep =sp.repeat(simparams['Timevec'],nlocout,axis=0)
        self.Sphere_Coords_Out_Rep =sp.tile(self.Sphere_Coords_Out,(ntout,1))
        self.RSTMat = makematPA(ionoin.Sphere_Coords,ionoin.Time_Vector)