halo.py

"""

halo
====

Implements halo catalogue functions. If you have a supported halo
catalogue on disk or a halo finder installed and correctly configured,
you can access a halo catalogue through f.halos() where f is a
SimSnap.

See the `halo tutorial
<http://pynbody.github.io/pynbody/tutorials/halos.html>`_ for some
examples.

"""

import numpy as np
import weakref
import os.path
import glob
import re
import copy
import sys
from array import SimArray
import gzip
import logging
from . import snapshot, util, config, config_parser, units
from .snapshot import gadget

logger = logging.getLogger("pynbody.halo")

class DummyHalo(object):

    def __init__(self):
        self.properties = {}


class Halo(snapshot.IndexedSubSnap):

    """
    Generic class representing a halo.
    """

    def __init__(self, halo_id, halo_catalogue, *args):
        super(Halo, self).__init__(*args)
        self._halo_catalogue = halo_catalogue
        self._halo_id = halo_id
        self._descriptor = "halo_" + str(halo_id)
        self.properties = copy.copy(self.properties)
        self.properties['halo_id'] = halo_id

    def is_subhalo(self, otherhalo):
        """
        Convenience function that calls the corresponding function in
        a halo catalogue.
        """

        return self._halo_catalogue.is_subhalo(self._halo_id, otherhalo._halo_id)


# ----------------------------#
# General HaloCatalogue class #
#-----------------------------#

class HaloCatalogue(object):

    """
    Generic halo catalogue object.
    """

    def __init__(self, sim):
        self._base = weakref.ref(sim)
        self._halos = {}
        self.lazy_off = util.ExecutionControl()

    def calc_item(self, i):
        if i in self._halos:  # and self._halos[i]() is not None :
            return self._halos[i]  # ()
        else:
            h = self._get_halo(i)
            self._halos[i] = h  # weakref.ref(h)
            return h

    def __len__(self):
        return len(self._halos)

    def __iter__(self):
        return self._halo_generator()

    def __getitem__(self, item):
        if isinstance(item, slice):
            for x in self._halo_generator(item.start,item.stop) : pass
            indices = item.indices(len(self._halos))
            res = [self.calc_item(i) for i in range(*indices)]
            return res
        else:
            return self.calc_item(item)

    @property
    def base(self):
        return self._base()

    def _halo_generator(self, i_start=None, i_stop=None) :
        if len(self) == 0 : return
        if i_start is None :
            try :
                self[0]
                i = 0
            except KeyError :
                i = 1
        else :
            i = i_start

        if i_stop is None :
            i_stop = len(self)

        while True:
            try:
                yield self[i]
                i+=1
                if len(self[i]) == 0: continue
            except RuntimeError:
                break
            if i == i_stop: raise StopIteration

    def _init_iord_to_fpos(self):
        if not hasattr(self, "_iord_to_fpos"):
            self._iord_to_fpos = np.empty(self.base['iord'].max()+1,dtype=np.int64)
            self._iord_to_fpos[self.base['iord']] = np.arange(len(self.base))

    def is_subhalo(self, childid, parentid):
        """Checks whether the specified 'childid' halo is a subhalo
        of 'parentid' halo.
        """
        if (childid in self._halos[parentid].properties['children']):
            return True
        else:
            return False

    def contains(self, haloid):
        if (haloid in self._halos):
            return True
        else:
            return False

    def __contains__(self, haloid):
        return self.contains(haloid)

    def get_group_array(self):
        """Return an array with an integer for each particle in the simulation
        indicating which halo that particle is associated with. If there are multiple
        levels (i.e. subhalos), the number returned corresponds to the lowest level, i.e.
        the smallest subhalo."""
        raise NotImplementedError

    @staticmethod
    def _can_load(self):
        return False

    @staticmethod
    def _can_run(self):
        return False

class RockstarIntermediateCatalogue(HaloCatalogue):
    """Reader for Rockstar intermediate catalogues as generated by
    Michael Tremmel's tool"""

    _halo_type = np.dtype([('id',np.int64),('num_p',np.int64),('indstart',np.int64)])
    _part_type = np.dtype('int64')

    def __init__(self, sim, sort=True, correct=False):
        assert isinstance(sim,snapshot.SimSnap)
        self._correct=correct
        HaloCatalogue.__init__(self, sim)

        self._halos = {}

        self._index_filename = sim.filename+".rockstar.halos"
        self._particles_filename = sim.filename+".rockstar.halo_particles_fpos"

        self._init_index()
        if sort:
            self._sort_index()

    def __len__(self):
        return self._nhalos

    def _get_particles_for_halo(self, num):
        halo_info = self._halo_info[num]
        with util.open_(self._particles_filename) as f:
            f.seek(halo_info['indstart']*self._part_type.itemsize)
            halo_ptcls=np.fromfile(f,dtype=self._part_type,count=halo_info['num_p']-1)
            halo_ptcls.sort()


        return halo_ptcls

    def _add_halo_id(self, halo, num):
        halo.properties['rockstar_halo_id']=self._halo_info[num]['id']
        if self._correct is True:
            with util.open_(self._particles_filename) as f:
                f.seek(self._halo_info['indstart'][num]*self._part_type.itemsize)
                one_ptcl = np.fromfile(f,dtype=self._part_type,count=1)
                from . import load
                onep = load(self.base.filename, take = one_ptcl)
                while not onep['grpid']:
                    one_ptcl = np.fromfile(f,dtype=self._part_type,count=1)
                    onep = load(self.base.filename, take = one_ptcl)
                halo.properties['rockstar_halo_id'] = onep['grpid'][0]

    def _get_halo(self, i):
        if self.base is None:
            raise RuntimeError("Parent SimSnap has been deleted")

        halo_ptcls = self._get_particles_for_halo(i)
        h = Halo(i, self, self.base, halo_ptcls)
        self._add_halo_id(h,i)
        return h

    def load_copy(self, i):
        """Load a fresh SimSnap with only the particles in halo i"""
        if i>=len(self):
            raise KeyError, "No such halo"

        from . import load
        halo = load(self.base.filename, take=self._get_particles_for_halo(i))
        self._add_halo_id(halo, i)
        return halo

    @staticmethod
    def _can_load(sim,*args,**kwargs):
        return os.path.exists(sim.filename+".rockstar.halos") and \
               os.path.exists(sim.filename+".rockstar.halo_particles_fpos")

    def _init_index(self):
        with open(self._index_filename,"rb") as f:
            self._nhalos = np.fromfile(f,np.int64,1)[0]
            #self._halo_info = np.fromfile(f,self._halo_type,self._nhalos)
            self._halo_info = np.fromfile(f,self._halo_type,-1)
            ok, = np.where(self._halo_info['indstart']>=0)
            self._halo_info = self._halo_info[ok]

    def _sort_index(self):
        self._halo_info[::-1].sort(order='num_p')

    def get_fam_group_array(self, family='star'):
        if family == 'star':
            target = self.base.star
        if family == 'gas':
            target = self.base.gas
        if family == 'dark':
            target = self.base.dark
        if family == 'BH':
            temptarget = self.base.star
            target = temptarget[(temptarget['tform']<0)]
        if family not in ['star','gas','dark', 'BH']:
            raise TypeError("family input not understood. Must be star, gas, or dark")
        ar = -1 * np.ones(len(target))
        if family != 'dark':
            fpos_ar = target.get_index_list(self.base)
        ngas = len(self.base.gas)
        ndark = len(self.base.dark)
        with util.open_(self._particles_filename) as f:
            for i in range(len(self._halo_info)):
                #print i, len(self._halo_info)
                #if i%100 == 0: print float(i)/float(len(self._halo_info))*100, "% done"
                f.seek(self._halo_info[i]['indstart']*self._part_type.itemsize)
                halo_ptcls=np.fromfile(f,dtype=self._part_type,count=self._halo_info[i]['num_p']-1)
                if family == 'gas':
                    halo_ptcls = halo_ptcls[(halo_ptcls>=ndark)&(halo_ptcls<ndark+ngas)]
                if family == 'dark':
                    halo_ptcls = halo_ptcls[(halo_ptcls<ndark)]
                if family == 'star' or family == 'BH':
                    halo_ptcls = halo_ptcls[(halo_ptcls>=ndark+ngas)]
                if family != 'dark':
                    match, = np.where(np.in1d(fpos_ar, halo_ptcls))
                    ar[match] = i
                else:
                    ar[halo_ptcls] = i
        return ar.astype(np.int)






class RockstarCatalogue(HaloCatalogue):
    def __init__(self, sim, dummy=False, filename=None, sort=False):
        """Initialize a RockstarCatalogue.

        **kwargs** :


        *dummy*: if True, the particle file is not loaded, and all
                 halos returned are just dummies (with the correct
                 properties dictionary loaded). Use load_copy to get
                 the actual data in this case.
                 
        *sort*: if True, resort the halos into descending order of
                particle number. Otherwise, leave in RockStar output order.

        """
        HaloCatalogue.__init__(self, sim)

        self._dummy = dummy

        if filename is not None:
            self._files = filename
        else:
            self._files = glob.glob(os.path.join(os.path.dirname(sim.filename),'halos*.bin'))

        self._cpus = [RockstarCatalogueOneCpu(sim,dummy,file_i) for file_i in self._files]
        self._prune_files_from_wrong_scalefactor()
        self._cpus[0]._init_iord_to_fpos()
        for cpu in self._cpus:
            cpu._iord_to_fpos = self._cpus[0]._iord_to_fpos
        
        self._init_index_ar()
        if sort:
            self._sort_index_ar()

    def _prune_files_from_wrong_scalefactor(self):
        new_cpus = []
        for cpu in self._cpus:
            if abs(self.base.properties['a']-cpu._head['scale'])<1e-6:
                new_cpus.append(cpu)
        self._cpus = new_cpus

    def _pass_on(self, function, *args, **kwargs):
        if self._index_ar is not None:
            cpui, hi = self._cpus[self._index_ar[args[0],0]], self._index_ar[args[0],1]
            return function(cpui,hi,*args[1:],**kwargs)
        for i in self._cpus:
            try:
                return function(i,*args,**kwargs)
            except KeyError:
                pass

    def __getitem__(self, k):
        return self._pass_on(RockstarCatalogueOneCpu.__getitem__, k)

    def load_copy(self, k):
        return self._pass_on(RockstarCatalogueOneCpu.load_copy, k)

    def __len__(self):
        return sum([len(x) for x in self._cpus])
    
    def _run_rockstar(self, sim):
        import pynbody
        fileformat = 'TIPSY'
        if (sim is pynbody.gadget.GadgetSnap):
            fileformat = 'GADGET'
        import pynbody.units as units

        # find AHFstep
        groupfinder = config_parser.get('RockstarCatalogue', 'Path')

        if groupfinder == 'None':
            for directory in os.environ["PATH"].split(os.pathsep):
                ahfs = glob.glob(os.path.join(directory, "rockstar-galaxies"))
                for iahf, ahf in enumerate(ahfs):
                    if ((len(ahfs) > 1) & (iahf != len(ahfs)-1) &
                            (os.path.basename(ahf) == 'rockstar-galaxies')):
                        continue
                    else:
                        groupfinder = ahf
                        break

        if not os.path.exists(groupfinder):
            raise RuntimeError("Path to Rockstar (%s) is invalid" % groupfinder)

        f = open('quickstart.cfg', 'w')
        print >>f, config_parser.get('RockstarCatalogue', 'Config', vars={
                'format': fileformat,
                'partmass': sim.d['mass'].in_units('Msol h^-1',**sim.conversion_context()).min(),
                'expfac': sim.properties['a'],
                'hub': sim.properties['h'],
                'omega0': sim.properties['omegaM0'],
                'lambda0': sim.properties['omegaL0'],
                'boxsize': sim['pos'].units.ratio('Mpc a h^-1', **sim.conversion_context()),
                'vunit': sim['vel'].units.ratio('km s^-1 a', **sim.conversion_context()),
                'munit': sim['mass'].units.ratio('Msol h^-1', **sim.conversion_context()),
                'softening': sim.s['eps'].in_units('Mpc a h^-1', **sim.conversion_context()).min()
            })

        f.close()

        if (not os.path.exists(sim._filename)):
            os.system("gunzip "+sim._filename+".gz")
        # determine parallel possibilities

        if os.path.exists(groupfinder):
            # run it
            if config['verbose']:
                print "RockstarCatalogue: running %s"%groupfinder
            os.system(groupfinder+" -c quickstart.cfg "+sim._filename)
            return

    def _init_index_ar(self):
        index_ar = np.empty((len(self),2),dtype=np.int32)

        for cpu_id, cpu in enumerate(self._cpus):
            i0 = cpu._halo_min
            i1 = cpu._halo_max
            index_ar[i0:i1,0]=cpu_id
            index_ar[i0:i1,1]=np.arange(i0,i1,dtype=np.int32)

        self._index_ar = index_ar
    

    def _sort_index_ar(self):
        num_ar = np.empty(len(self))

        for cpu_id, cpu in enumerate(self._cpus):
            i0 = cpu._halo_min
            i1 = cpu._halo_max
            num_ar[i0:i1]=self._cpus[cpu_id]._halo_lens

        num_ar = np.argsort(num_ar)[::-1]
        self._index_ar = self._index_ar[num_ar]


    @staticmethod
    def _can_run(sim):
        if config_parser.getboolean('RockstarCatalogue', 'AutoRun'):
            if config_parser.get('RockstarCatalogue', 'Path') == 'None':
                for directory in os.environ["PATH"].split(os.pathsep):
                    if (len(glob.glob(os.path.join(directory, "rockstar-galaxies"))) > 0):
                        return True
            else:
                path = config_parser.get('RockstarCatalogue', 'Path')
                return os.path.exists(path)
        return False



        
class RockstarCatalogueOneCpu(HaloCatalogue):
    """
    Class to handle catalogues produced by Rockstar
    """

    head_type = np.dtype([('magic',np.uint64),('snap',np.int64),
                          ('chunk',np.int64),('scale','f'),
                          ('Om','f'),('Ol','f'),('h0','f'),
                          ('bounds','f',6),('num_halos',np.int64),
                          ('num_particles',np.int64),('box_size','f'),
                          ('particle_mass','f'),('particle_type',np.int64),
                          ('format_revision',np.int32),
                          ('rockstar_version',np.str_,12)])

    halo_type = np.dtype([('id',np.int64),('pos','f',3),('vel','f',3),
                          ('corevel','f',3),('bulkvel','f',3),('m','f'),
                          ('r','f'),
                          ('child_r','f'),('vmax_r','f'),('mgrav','f'),
                          ('vmax','f'),('rvmax','f'),('rs','f'),
                          ('klypin_rs','f'),('vrms','f'),('J','f',3),
                          ('energy','f'),('spin','f'),('alt_m','f',4),
                          ('Xoff','f'),('Voff','f'),('b_to_a','f'),
                          ('c_to_a','f'),('A','f',3),('b_to_a2','f'),
                          ('c_to_a2','f'),('A2','f',3),('bullock_spin','f'),
                          ('kin_to_pot','f'),('m_pe_b','f'),('m_pe_d','f'),
                          ('dum',np.str_,4),
                          ('num_p',np.int64),('num_child_particles',np.int64),
                          ('p_start',np.int64),('desc',np.int64),
                          ('flags',np.int64),('n_core',np.int64),
                          ('min_pos_err','f'),('min_vel_err','f'),
                          ('min_bulkvel_err','f'),('type',np.int32),
                          ('sm','f'),('gas','f'),('bh','f'),
                          ('peak_density','f'),('av_density','f'),
                          ('odum',np.str_,4)])


    def __init__(self, sim, dummy=False, filename=None):
        """Initialize a RockstarCatalogue.

        **kwargs** :


        *dummy*: if True, the particle file is not loaded, and all
                 halos returned are just dummies (with the correct
                 properties dictionary loaded). Use load_copy to get
                 the actual data in this case.


        """

        import os.path
        if not self._can_load(sim):
            self._run_rockstar(sim)
        HaloCatalogue.__init__(self,sim)

        self._dummy = dummy

        if filename is not None: self._rsFilename = filename
        else:
            self._rsFilename = util.cutgz(glob.glob('halos*.bin')[0])

        try :
            f = util.open_(self._rsFilename)
        except IOError:
            raise IOError("Halo catalogue not found -- check the file name of catalogue data or try specifying a catalogue using the filename keyword")

        with f:
            self._head = np.fromstring(f.read(self.head_type.itemsize),
                                       dtype=self.head_type)
            unused = f.read(256 - self._head.itemsize)

            self._nhalos = self._head['num_halos'][0]

            self._load_rs_halos(f,sim)
            
        


    def __len__(self):
        return len(self._halo_lens)

    def calc_item(self, i):
        if self.base is None:
            raise RuntimeError("Parent SimSnap has been deleted")

        self._load_rs_halo_if_required(i)
        self._load_rs_particles_for_halo_if_required(i)

        return self._halos[i]

    def load_copy(self, i):
        """Load a fresh SimSnap with only the particles in halo i"""
        if i<self._halo_min or i>=self._halo_max:
            raise KeyError, "No such halo"

        from . import load
        return load(self.base.filename, take=self._get_particles_for_halo(i))

    @staticmethod
    def _can_load(sim,**kwargs):
        for file in glob.glob('halos*.bin'):
            if os.path.exists(file):
                return True
        return False



    def make_grp(self):
        """
        Creates a 'grp' array which labels each particle according to
        its parent halo.
        """
        try:
            self.base['grp']
        except:
            self.base['grp'] = np.zeros(len(self.base),dtype='i')

        for halo in self._halos.values():
            halo[name][:] = halo._halo_id

        if config['verbose']:  print "writing %s"%(self._base().filename+'.grp')
        self._base().write_array('grp',overwrite=True,binary=False)
        
    def _setup_children(self):
        """
        Creates a 'children' array inside each halo's 'properties'
        listing the halo IDs of its children. Used in case the reading
        of substructure data from the AHF-supplied _substructure file
        fails for some reason.
        """

        for i in xrange(self._nhalos):
            self._halos[i+1].properties['children'] = []

        for i in xrange(self._nhalos):
            host = self._halos[i+1].properties.get('hostHalo', -2)
            if host > -1:
                try:
                    self._halos[host+1].properties['children'].append(i+1)
                except KeyError:
                    pass





    def _load_rs_halo_if_required(self, i):
        if i not in self._halos:
            self._halos[i] = self._get_dummy_for_halo(i)
            
    def _get_dummy_for_halo(self, n):
        if n<self._halo_min or n>=self._halo_max:
            raise KeyError, "No such halo"
        
        with util.open_(self._rsFilename) as f:
            f.seek(self._haloprops_offset+(n-self._halo_min)*self.halo_type.itemsize)
            halo_data = np.fromfile(f, dtype=self.halo_type, count=1)

        hn = DummyHalo()
        # TODO: properties are in Msun / h, Mpc / h
        hn.properties = dict(zip(halo_data.dtype.names,halo_data[0]))
        return hn  
                                                
        
    def _load_rs_halos(self, f, sim):
        self._haloprops_offset = f.tell()
        self._halo_offsets = np.empty(self._head['num_halos'],dtype=np.int64)
        self._halo_lens = np.empty(self._head['num_halos'],dtype=np.int64)
        offset = self._haloprops_offset+self.halo_type.itemsize*self._head['num_halos']
        
        
        self._halo_min = int(np.fromfile(f, dtype=self.halo_type, count=1)['id'])
        self._halo_max = self._halo_min+self._head['num_halos']

        f.seek(self._haloprops_offset)
        
        this_id = self._halo_min
        
        for h in xrange(self._head['num_halos']):
            halo_data =np.fromfile(f, dtype=self.halo_type, count=1)
            self._halo_offsets[this_id-self._halo_min] = int(offset)
            self._halo_lens[this_id-self._halo_min] = int(halo_data['num_p'])
            offset+=halo_data['num_p']*np.dtype('int64').itemsize
            this_id+=1
            

    def _get_particles_for_halo(self, num):
        self._init_iord_to_fpos()
        with util.open_(self._rsFilename) as f:
            f.seek(self._halo_offsets[num-self._halo_min])
            halo_ptcls=np.fromfile(f,dtype=np.int64,count=self._halo_lens[num-self._halo_min])
            halo_ptcls = self._iord_to_fpos[halo_ptcls]
            halo_ptcls.sort()
            
        return halo_ptcls
    
    def _load_rs_particles_for_halo(self, num):
        halo_ptcls = self._get_particles_for_halo(num)

        properties_from_proxy = self._halos[num].properties

        self._halos[num]=Halo(num, self, self.base, halo_ptcls)
        self._halos[num]._descriptor = "halo_"+str(num)

        self._halos[num].properties.update(properties_from_proxy)

    def _load_rs_particles_for_halo_if_required(self, num):
        if isinstance(self._halos[num], DummyHalo) and not self._dummy:
            self._load_rs_particles_for_halo(num)



    def _load_ahf_substructure(self, filename):
        f = util.open_(filename)

        for i in xrange(len(self._halos)):

            haloid, nsubhalos = [int(x) for x in f.readline().split()]
            self._halos[haloid+1].properties['children'] = [
                int(x)+1 for x in f.readline().split()]

        f.close()




    def writestat(self, outfile=None, hubble=None):
        """
        write a condensed skid.stat style ascii file from ahf_halos
        file.  header + 1 halo per line. should reproduce `Alyson's
        idl script' except does not do last 2 columns (Is it a
        satellite?) and (Is central halo is `false'ly split?).  output
        units are set to Mpc Msun, km/s.

        user can specify own hubble constant hubble=(H0/(100
        km/s/Mpc)), ignoring the snaphot arg for hubble constant
        (which sometimes has a large roundoff error).
        """
        s = self._base()
        mindarkmass = min(s.dark['mass'])

        if hubble is None:
            hubble = s.properties['h']

        if outfile is None: outfile = self._base().filename+'.stat'
        print "write stat file to ", outfile
        fpout = open(outfile, "w")
        header = "#Grp  N_tot     N_gas      N_star    N_dark    Mvir(M_sol)       Rvir(kpc)       GasMass(M_sol) StarMass(M_sol)  DarkMass(M_sol)  V_max  R@V_max  VelDisp    Xc   Yc   Zc   VXc   VYc   VZc   Contam   Satellite?   False?   ID_A"
        print >> fpout, header
        for ii in np.arange(self._nhalos)+1:
            print '%d '%ii,
            sys.stdout.flush()
            h = self[ii].properties  # halo index starts with 1 not 0
##  'Contaminated'? means multiple dark matter particle masses in halo)"
            icontam = np.where(self[ii].dark['mass'] > mindarkmass)
            if (len(icontam[0]) > 0):
                contam = "contam"
            else:
                contam = "clean"
## may want to add implement satellite test and false central breakup test.
            ss = "     "  # can adjust column spacing
            outstring = str(ii)+ss
            outstring += str(len(self[ii]))+ss+str(len(self[ii].g))+ss
            outstring += str(len(self[ii].s)) + ss+str(len(self[ii].dark))+ss
            outstring += str(h['m']/hubble)+ss+str(h['r']/hubble)+ss
            outstring += str(self[ii].g['mass'].in_units('Msol').sum())+ss
            outstring += str(self[ii].s['mass'].in_units('Msol').sum())+ss
            outstring += str(self[ii].d['mass'].in_units('Msol').sum())+ss
            outstring += str(h['vmax'])+ss+str(h['vmax_r']/hubble)+ss
            outstring += str(h['vrms'])+ss
        ## pos: convert kpc/h to mpc (no h).
            outstring += str(h['pos'][0][0]/hubble)+ss
            outstring += str(h['pos'][0][1]/hubble)+ss
            outstring += str(h['pos'][0][2]/hubble)+ss
            outstring += str(h['vel'][0][0])+ss+str(h['vel'][0][1])+ss
            outstring += str(h['vel'][0][2])+ss
            outstring += contam+ss
            outstring += "unknown" + \
                ss  # unknown means sat. test not implemented.
            outstring += "unknown"+ss  # false central breakup.
            print >> fpout, outstring
        fpout.close()

    def writetipsy(self, outfile=None, hubble=None):
        """
        write halos to tipsy file (write as stars) from ahf_halos
        file.  returns a shapshot where each halo is a star particle.

        user can specify own hubble constant hubble=(H0/(100
        km/s/Mpc)), ignoring the snaphot arg for hubble constant
        (which sometimes has a large roundoff error).
        """
        from . import analysis
        from . import tipsy
        from .analysis import cosmology
        from snapshot import _new as new
        import math
        s = self._base()
        if outfile is None: outfile = s.filename+'.gtp'
        print "write tipsy file to ", outfile
        sout = new(star=self._nhalos)  # create new tipsy snapshot written as halos.
        sout.properties['a'] = s.properties['a']
        sout.properties['z'] = s.properties['z']
        sout.properties['boxsize'] = s.properties['boxsize']
        if hubble is None: hubble = s.properties['h']
        sout.properties['h'] = hubble
    ### ! dangerous -- rho_crit function and unit conversions needs simplifying
        rhocrithhco = cosmology.rho_crit(s, z=0, unit="Msol Mpc^-3 h^2")
        lboxkpc = sout.properties['boxsize'].ratio("kpc a")
        lboxkpch = lboxkpc*sout.properties['h']
        lboxmpch = lboxkpc*sout.properties['h']/1000.
        tipsyvunitkms = lboxmpch * 100. / (math.pi * 8./3.)**.5
        tipsymunitmsun = rhocrithhco * lboxmpch**3 / sout.properties['h']

        print "transforming ", self._nhalos, " halos into tipsy star particles"
        for ii in xrange(self._nhalos):
            h = self[ii+1].properties
            sout.star[ii]['mass'] = h['m']/hubble / tipsymunitmsun
            ## tipsy units: box centered at 0. (assume 0<=x<=1)
            sout.star[ii]['x'] = h['pos'][0][0]/lboxmpch - 0.5
            sout.star[ii]['y'] = h['pos'][0][1]/lboxmpch - 0.5
            sout.star[ii]['z'] = h['pos'][0][2]/lboxmpch - 0.5
            sout.star[ii]['vx'] = h['vel'][0][0]/tipsyvunitkms
            sout.star[ii]['vy'] = h['vel'][0][1]/tipsyvunitkms
            sout.star[ii]['vz'] = h['vel'][0][2]/tipsyvunitkms
            sout.star[ii]['eps'] = h['r']/lboxkpch
            sout.star[ii]['metals'] = 0.
            sout.star[ii]['phi'] = 0.
            sout.star[ii]['tform'] = 0.
        print "writing tipsy outfile %s"%outfile
        sout.write(fmt=tipsy.TipsySnap, filename=outfile)
        return sout

    def writehalos(self, hubble=None, outfile=None):
        s = self._base()
        if outfile is None:
            statoutfile = s.filename+".rockstar.stat"
            tipsyoutfile = s.filename+".rockstar.gtp"
        else:
            statoutfile = outfile+'.stat'
            gtpoutfile = outfile+'.gtp'
        self.make_grp()
        self.writestat(statoutfile, hubble=hubble)
        shalos = self.writetipsy(gtpoutfile, hubble=hubble)
        return shalos


#--------------------------#
# AHF Halo Catalogue class #
#--------------------------#

class AHFCatalogue(HaloCatalogue):

    """
    Class to handle catalogues produced by Amiga Halo Finder (AHF).
    """

    def __init__(self, sim, make_grp=None, dummy=False, use_iord=None, ahf_basename=None):
        """Initialize an AHFCatalogue.

        **kwargs** :

        *make_grp*: if True a 'grp' array is created in the underlying
                    snapshot specifying the lowest level halo that any
                    given particle belongs to. If it is False, no such
                    array is created; if None, the behaviour is
                    determined by the configuration system.

        *dummy*: if True, the particle file is not loaded, and all
                 halos returned are just dummies (with the correct
                 properties dictionary loaded). Use load_copy to get
                 the actual data in this case.

        *use_iord*: if True, the particle IDs in the Amiga catalogue
                    are taken to refer to the iord array. If False,
                    they are the particle offsets within the file. If
                    None, the parameter defaults to True for
                    GadgetSnap, False otherwise.

        *ahf_basename*: specify the basename of the AHF halo catalog
                        files - the code will append 'halos',
                        'particles', and 'substructure' to this
                        basename to load the catalog data.

        """

        import os.path
        if not self._can_load(sim):
            self._run_ahf(sim)

        HaloCatalogue.__init__(self,sim)

        if use_iord is None:
            use_iord = isinstance(sim.ancestor, gadget.GadgetSnap)

        self._use_iord = use_iord

        self._dummy = dummy

        if ahf_basename is not None:
            self._ahfBasename = ahf_basename
        else:
            self._ahfBasename = util.cutgz(
                glob.glob(sim._filename + '*z*AHF_halos*')[0])[:-5]

        try:
            f = util.open_(self._ahfBasename + 'halos')
        except IOError:
            raise IOError(
                "Halo catalogue not found -- check the base name of catalogue data or try specifying a catalogue using the ahf_basename keyword")

        for i, l in enumerate(f):
            pass
        self._nhalos = i
        f.close()

        logger.info("AHFCatalogue loading particles")

        self._load_ahf_particles(self._ahfBasename + 'particles')

        logger.info("AHFCatalogue loading halos")

        self._load_ahf_halos(self._ahfBasename + 'halos')

        if os.path.isfile(self._ahfBasename + 'substructure'):
            logger.info("AHFCatalogue loading substructure")

            self._load_ahf_substructure(self._ahfBasename + 'substructure')
        else:
            self._setup_children()

        if make_grp is None:
            make_grp = config_parser.getboolean('AHFCatalogue', 'AutoGrp')

        if make_grp:
            self.make_grp()

        if config_parser.getboolean('AHFCatalogue', 'AutoPid'):
            sim['pid'] = np.arange(0, len(sim))

        logger.info("AHFCatalogue loaded")

    def make_grp(self, name='grp'):
        """
        Creates a 'grp' array which labels each particle according to
        its parent halo.
        """
        self.base[name] = self.get_group_array()

    def get_group_array(self, top_level=False):
        ar = np.zeros(len(self.base), dtype=int)
        if top_level is False:
            for halo in self._halos.values():
                ar[halo.get_index_list(self.base)] = halo._halo_id
        else:
            for halo in self._halos.values()[::-1]:
                ar[halo.get_index_list(self.base)] = halo._halo_id
        return ar

    def _setup_children(self):
        """
        Creates a 'children' array inside each halo's 'properties'
        listing the halo IDs of its children. Used in case the reading
        of substructure data from the AHF-supplied _substructure file
        fails for some reason.
        """

        for i in xrange(self._nhalos):
            self._halos[i + 1].properties['children'] = []

        for i in xrange(self._nhalos):
            host = self._halos[i + 1].properties.get('hostHalo', -2)
            if host > -1:
                try:
                    self._halos[host + 1].properties['children'].append(i + 1)
                except KeyError:
                    pass

    def _get_halo(self, i):
        if self.base is None:
            raise RuntimeError("Parent SimSnap has been deleted")

        return self._halos[i]

    @property
    def base(self):
        return self._base()

    def load_copy(self, i):
        """Load the a fresh SimSnap with only the particle in halo i"""

        from . import load

        f = util.open_(self._ahfBasename + 'particles')

        if self.isnew:
            nhalos = int(f.readline())
        else:
            nhalos = self._nhalos

        for h in xrange(i):
            ids = self._load_ahf_particle_block(f)

        f.close()

        return load(self.base.filename, take=ids)

    def _load_ahf_particle_block(self, f):
        """Load the particles for the next halo described in particle file f"""
        ng = len(self.base.gas)
        nds = len(self.base.dark) + len(self.base.star)
        nparts = int(f.readline().split()[0])

        if self.isnew:
            if not isinstance(f, gzip.GzipFile):
                data = (np.fromfile(
                    f, dtype=int, sep=" ", count=nparts * 2).reshape(nparts, 2))[:, 0]
            else:
                # unfortunately with gzipped files there does not
                # seem to be an efficient way to load nparts lines
                data = np.zeros(nparts, dtype=int)
                for i in xrange(nparts):
                    data[i] = int(f.readline().split()[0])

            if self._use_iord:
                data = self._iord_to_fpos[data]
            else:
                hi_mask = data >= nds
                data[np.where(hi_mask)] -= nds
                data[np.where(~hi_mask)] += ng
        else:
            if not isinstance(f, gzip.GzipFile):
                data = np.fromfile(f, dtype=int, sep=" ", count=nparts)
            else:
                # see comment above on gzipped files
                data = np.zeros(nparts, dtype=int)
                for i in xrange(nparts):
                    data[i] = int(f.readline())
        data.sort()
        return data

    def _load_ahf_particles(self, filename):
        if self._use_iord:
            self._iord_to_fpos = np.zeros(self.base['iord'].max()+1,dtype=int)
            self._iord_to_fpos[self.base['iord']] = np.arange(len(self._base()))

        f = util.open_(filename)
        if filename.split("z")[-2][-1] is ".":
            self.isnew = True
        else:
            self.isnew = False

        if self.isnew:
            nhalos = int(f.readline())
        else:
            nhalos = self._nhalos

        if not self._dummy:
            for h in xrange(nhalos):
                self._halos[h + 1] = Halo(
                    h + 1, self, self.base, self._load_ahf_particle_block(f))
                self._halos[h + 1]._descriptor = "halo_" + str(h + 1)
        else:
            for h in xrange(nhalos):
                self._halos[h + 1] = DummyHalo()

        f.close()

    def _load_ahf_halos(self, filename):
        f = util.open_(filename,"rt")
        # get all the property names from the first, commented line
        # remove (#)
        keys = [re.sub('\([0-9]*\)', '', field)
                for field in f.readline().split()]
        # provide translations
        for i, key in enumerate(keys):
            if self.isnew:
                if(key == '#npart'):
                    keys[i] = 'npart'
            else:
                if(key == '#'):
                    keys[i] = 'dumb'
            if(key == 'a'):
                keys[i] = 'a_axis'
            if(key == 'b'):
                keys[i] = 'b_axis'
            if(key == 'c'):
                keys[i] = 'c_axis'
            if(key == 'Mvir'):
                keys[i] = 'mass'

        if self.isnew:
            # fix for column 0 being a non-column in some versions of the AHF
            # output
            if keys[0] == '#':
                keys = keys[1:]

        for h, line in enumerate(f):
            values = [float(x) if '.' in x or 'e' in x or 'nan' in x else int(
                x) for x in line.split()]
            # XXX Unit issues!  AHF uses distances in Mpc/h, possibly masses as
            # well
            for i, key in enumerate(keys):
                if self.isnew:
                    self._halos[h + 1].properties[key] = values[i]
                else:
                    self._halos[h + 1].properties[key] = values[i - 1]
        f.close()

    def _load_ahf_substructure(self, filename):
        f = util.open_(filename)
        # nhalos = int(f.readline())  # number of halos?  no, some crazy number
        # that we will ignore
        for i in xrange(len(self._halos)):
            try:
                haloid, nsubhalos = [int(x) for x in f.readline().split()]
                self._halos[haloid + 1].properties['children'] = [
                    int(x) + 1 for x in f.readline().split()]
                for ichild in self._halos[haloid + 1].properties['children']:
                    self._halos[ichild].properties['parentid'] = haloid+1
            except KeyError:
                pass
            except ValueError:
                break
        f.close()

    def writegrp(self, grpoutfile=False):
        """
        simply write a skid style .grp file from ahf_particles
        file. header = total number of particles, then each line is
        the halo id for each particle (0 means free).
        """
        snapshot = self[1].ancestor
        try:
            snapshot['grp']
        except:
            self.make_grp()
        if not grpoutfile:
            grpoutfile = snapshot.filename + '.grp'
        logger.info("Writing grp file to %s" % grpoutfile)
        fpout = open(grpoutfile, "w")
        print >> fpout, len(snapshot['grp'])

        # writing 1st to a string sacrifices memory for speed.
        # but this is much faster than numpy.savetxt (could make an option).
        # it is assumed that max halo id <= nhalos (i.e.length of string is set
        # len(str(nhalos))
        stringarray = snapshot['grp'].astype(
            '|S' + str(len(str(self._nhalos))))
        outstring = "\n".join(stringarray)
        print >> fpout, outstring
        fpout.close()

    def writestat(self, snapshot, halos, statoutfile, hubble=None):
        """
        write a condensed skid.stat style ascii file from ahf_halos
        file.  header + 1 halo per line. should reproduce `Alyson's
        idl script' except does not do last 2 columns (Is it a
        satellite?) and (Is central halo is `false'ly split?).  output
        units are set to Mpc Msun, km/s.

        user can specify own hubble constant hubble=(H0/(100
        km/s/Mpc)), ignoring the snaphot arg for hubble constant
        (which sometimes has a large roundoff error).
        """
        s = snapshot
        mindarkmass = min(s.dark['mass'])

        if hubble is None:
            hubble = s.properties['h']

        outfile = statoutfile
        logger.info("Writing stat file to %s" % statoutfile)
        fpout = open(outfile, "w")
        header = "#Grp  N_tot     N_gas      N_star    N_dark    Mvir(M_sol)       Rvir(kpc)       GasMass(M_sol) StarMass(M_sol)  DarkMass(M_sol)  V_max  R@V_max  VelDisp    Xc   Yc   Zc   VXc   VYc   VZc   Contam   Satellite?   False?   ID_A"
        print >> fpout, header
        nhalos = halos._nhalos
        for ii in xrange(nhalos):
            h = halos[ii + 1].properties  # halo index starts with 1 not 0
            # 'Contaminated'? means multiple dark matter particle masses in halo)"
            icontam = np.where(halos[ii + 1].dark['mass'] > mindarkmass)
            if (len(icontam[0]) > 0):
                contam = "contam"
            else:
                contam = "clean"
            # may want to add implement satellite test and false central
            # breakup test.

            n_dark = h['npart'] - h['n_gas'] - h['n_star']
            M_dark = h['mass'] - h['M_gas'] - h['M_star']
            ss = "     "  # can adjust column spacing
            outstring = str(int(h['halo_id'])) + ss
            outstring += str(int(h['npart'])) + ss + str(int(h['n_gas'])) + ss
            outstring += str(int(h['n_star'])) + ss + str(int(n_dark)) + ss
            outstring += str(h['mass'] / hubble) + ss + \
                str(h['Rvir'] / hubble) + ss
            outstring += str(h['M_gas'] / hubble) + ss + \
                str(h['M_star'] / hubble) + ss
            outstring += str(M_dark / hubble) + ss
            outstring += str(h['Vmax']) + ss + str(h['Rmax'] / hubble) + ss
            outstring += str(h['sigV']) + ss
            # pos: convert kpc/h to mpc (no h).
            outstring += str(h['Xc'] / hubble / 1000.) + ss
            outstring += str(h['Yc'] / hubble / 1000.) + ss
            outstring += str(h['Zc'] / hubble / 1000.) + ss
            outstring += str(h['VXc']) + ss + \
                str(h['VYc']) + ss + str(h['VZc']) + ss
            outstring += contam + ss
            outstring += "unknown" + \
                ss  # unknown means sat. test not implemented.
            outstring += "unknown" + ss  # false central breakup.
            print >> fpout, outstring
        fpout.close()
        return 1

    def writetipsy(self, snapshot, halos, tipsyoutfile, hubble=None):
        """
        write halos to tipsy file (write as stars) from ahf_halos
        file.  returns a shapshot where each halo is a star particle.

        user can specify own hubble constant hubble=(H0/(100
        km/s/Mpc)), ignoring the snaphot arg for hubble constant
        (which sometimes has a large roundoff error).
        """
        from . import analysis
        from . import tipsy
        from .analysis import cosmology
        from snapshot import _new as new
        import math
        s = snapshot
        outfile = tipsyoutfile
        nhalos = halos._nhalos
        nstar = nhalos
        sout = new(star=nstar)  # create new tipsy snapshot written as halos.
        sout.properties['a'] = s.properties['a']
        sout.properties['z'] = s.properties['z']
        sout.properties['boxsize'] = s.properties['boxsize']
        if hubble is None:
            hubble = s.properties['h']
        sout.properties['h'] = hubble
    # ! dangerous -- rho_crit function and unit conversions needs simplifying
        rhocrithhco = cosmology.rho_crit(s, z=0, unit="Msol Mpc^-3 h^2")
        lboxkpc = sout.properties['boxsize'].ratio("kpc a")
        lboxkpch = lboxkpc * sout.properties['h']
        lboxmpch = lboxkpc * sout.properties['h'] / 1000.
        tipsyvunitkms = lboxmpch * 100. / (math.pi * 8. / 3.) ** .5
        tipsymunitmsun = rhocrithhco * lboxmpch ** 3 / sout.properties['h']

        for ii in xrange(nhalos):
            h = halos[ii + 1].properties
            sout.star[ii]['mass'] = h['mass'] / hubble / tipsymunitmsun
            # tipsy units: box centered at 0. (assume 0<=x<=1)
            sout.star[ii]['x'] = h['Xc'] / lboxkpch - 0.5
            sout.star[ii]['y'] = h['Yc'] / lboxkpch - 0.5
            sout.star[ii]['z'] = h['Zc'] / lboxkpch - 0.5
            sout.star[ii]['vx'] = h['VXc'] / tipsyvunitkms
            sout.star[ii]['vy'] = h['VYc'] / tipsyvunitkms
            sout.star[ii]['vz'] = h['VZc'] / tipsyvunitkms
            sout.star[ii]['eps'] = h['Rvir'] / lboxkpch
            sout.star[ii]['metals'] = 0.
            sout.star[ii]['phi'] = 0.
            sout.star[ii]['tform'] = 0.

        sout.write(fmt=tipsy.TipsySnap, filename=outfile)
        return sout

    def writehalos(self, snapshot, halos, hubble=None, outfile=None):
        """ Write the (ahf) halo catalog to disk.  This is really a
        wrapper that calls writegrp, writetipsy, writestat.  Writes
        .amiga.grp file (ascii group ids), .amiga.stat file (ascii
        halo catalog) and .amiga.gtp file (tipsy halo catalog).
        default outfile base simulation is same as snapshot s.
        function returns simsnap of halo catalog.
        """
        s = snapshot
        grpoutfile = s.filename + ".amiga.grp"
        statoutfile = s.filename + ".amiga.stat"
        tipsyoutfile = s.filename + ".amiga.gtp"
        halos.writegrp(s, halos, grpoutfile)
        halos.writestat(s, halos, statoutfile, hubble=hubble)
        shalos = halos.writetipsy(s, halos, tipsyoutfile, hubble=hubble)
        return shalos

    @staticmethod
    def _can_load(sim,**kwargs):
        for file in glob.glob(sim._filename + '*z*particles*'):
            if os.path.exists(file):
                return True
        return False

    def _run_ahf(self, sim):
        # if (sim is pynbody.tipsy.TipsySnap) :
        typecode = 90
        # elif (sim is pynbody.gadget.GadgetSnap):
        #   typecode = '60' or '61'
        import pynbody.units as units
        # find AHFstep

        groupfinder = config_parser.get('AHFCatalogue', 'Path')

        if groupfinder == 'None':
            for directory in os.environ["PATH"].split(os.pathsep):
                ahfs = glob.glob(os.path.join(directory, "AHF*"))
                for iahf, ahf in enumerate(ahfs):
                    # if there are more AHF*'s than 1, it's not the last one, and
                    # it's AHFstep, then continue, otherwise it's OK.
                    if ((len(ahfs) > 1) & (iahf != len(ahfs) - 1) &
                            (os.path.basename(ahf) == 'AHFstep')):
                        continue
                    else:
                        groupfinder = ahf
                        break

        if not os.path.exists(groupfinder):
            raise RuntimeError("Path to AHF (%s) is invalid" % groupfinder)

        if (os.path.basename(groupfinder) == 'AHFstep'):
            isAHFstep = True
        else:
            isAHFstep = False
        # build units file
        if isAHFstep:
            f = open('tipsy.info', 'w')
            f.write(str(sim.properties['omegaM0']) + "\n")
            f.write(str(sim.properties['omegaL0']) + "\n")
            f.write(str(sim['pos'].units.ratio(
                units.kpc, a=1) / 1000.0 * sim.properties['h']) + "\n")
            f.write(
                str(sim['vel'].units.ratio(units.km / units.s, a=1)) + "\n")
            f.write(str(sim['mass'].units.ratio(units.Msol)) + "\n")
            f.close()
            # make input file
            f = open('AHF.in', 'w')
            f.write(sim._filename + " " + str(typecode) + " 1\n")
            f.write(sim._filename + "\n256\n5\n5\n0\n0\n0\n0\n")
            f.close()
        else:
            # make input file
            f = open('AHF.in', 'w')

            lgmax = np.min([int(2 ** np.floor(np.log2(
                1.0 / np.min(sim['eps'])))), 131072])
            # hardcoded maximum 131072 might not be necessary

            print >>f, config_parser.get('AHFCatalogue', 'Config', vars={
                'filename': sim._filename,
                'typecode': typecode,
                'gridmax': lgmax
            })

            print >>f, config_parser.get('AHFCatalogue', 'ConfigTipsy', vars={
                'omega0': sim.properties['omegaM0'],
                'lambda0': sim.properties['omegaL0'],
                'boxsize': sim['pos'].units.ratio('Mpc a h^-1', **sim.conversion_context()),
                'vunit': sim['vel'].units.ratio('km s^-1 a', **sim.conversion_context()),
                'munit': sim['mass'].units.ratio('Msol h^-1', **sim.conversion_context()),
                'eunit': 0.03  # surely this can't be right?
            })

            f.close()

        if (not os.path.exists(sim._filename)):
            os.system("gunzip " + sim._filename + ".gz")
        # determine parallel possibilities

        if os.path.exists(groupfinder):
            # run it
            os.system(groupfinder + " AHF.in")
            return

    @staticmethod
    def _can_run(sim):
        if config_parser.getboolean('AHFCatalogue', 'AutoRun'):
            if config_parser.get('AHFCatalogue', 'Path') == 'None':
                for directory in os.environ["PATH"].split(os.pathsep):
                    if (len(glob.glob(os.path.join(directory, "AHF*"))) > 0):
                        return True
            else:
                path = config_parser.get('AHFCatalogue', 'Path')
                return os.path.exists(path)
        return False



#-----------------------------#
# General Grp Catalogue class #
#-----------------------------#

class GrpCatalogue(HaloCatalogue):
    """
    A generic catalogue using a .grp file to specify which particles
    belong to which group.
    """
    def __init__(self, sim, array='grp'):
        sim[array]
    # trigger lazy-loading and/or kick up a fuss if unavailable
        self._halos = {}
        self._array = array
        self._sorted = None
        HaloCatalogue.__init__(self,sim)

    def __len__(self):
        return self.base[self._array].max()

    def precalculate(self):
        """Speed up future operations by precalculating the indices
        for all halos in one operation. This is slow compared to
        getting a single halo, however."""
        self._sorted = np.argsort(
            self.base[self._array], kind='mergesort')  # mergesort for stability
        self._boundaries = util.find_boundaries(
            self.base[self._array][self._sorted])

    def get_group_array(self):
        return self.base[self._array]

    def _get_halo(self, i):
        if self.base is None:
            raise RuntimeError("Parent SimSnap has been deleted")

        no_exist = ValueError("Halo %s does not exist" % (str(i)))

        if self._sorted is None:
            # one-off selection
            x = Halo(i, self, self.base, np.where(self.base[self._array] == i))
            if len(x) == 0:
                raise no_exist
            x._descriptor = "halo_" + str(i)
            return x
        else:
            # pre-calculated
            if i >= len(self._boundaries) or i < 0:
                raise no_exist
            if self._boundaries[i] < 0:
                raise no_exist

            start = self._boundaries[i]
            if start is None:
                raise no_exist

            end = None
            j = i + 1
            while j < len(self._boundaries) and end is None:
                end = self._boundaries[j]
                j += 1

            x = Halo(i, self, self.base, self._sorted[start:end])
            x._descriptor = "halo_" + str(i)

            return x

    @property
    def base(self):
        return self._base()

    @staticmethod
    def _can_load(sim, arr_name='grp'):
        if (arr_name in sim.loadable_keys()) or (arr_name in sim.keys()) :
            return True
        else:
            return False


class AmigaGrpCatalogue(GrpCatalogue):
    def __init__(self, sim, arr_name='amiga.grp'):
        GrpCatalogue.__init__(self, sim, arr_name)

    @staticmethod
    def _can_load(sim,arr_name='amiga.grp'):
        return GrpCatalogue._can_load(sim, arr_name)


#-----------------------------------------------------------------------#
# SubFind Catalogue classes -- including classes for handling HDF format #
#-----------------------------------------------------------------------#

class SubfindCatalogue(HaloCatalogue):

    """
        Class to handle catalogues produced by the SubFind halo finder. 
        Can import FoF groups (default) or subhalos (by setting subs=True).
        Groups are sorted by mass (descending), most massive group is halo[0]. Subhalos within each group are also sorted by mass.
        IDs within groups are sorted by binding energy, UNLESS the IDs in the snapshots are not ordered starting from 0. In that case,
        setting order=False will generate the correct halo catalogue but the ordering by binding energy within each halo will be lost.
        Additional properties calculated by SubFind can be accessed via e.g. 'halo[0].properties'.
        make_grp=True generates a particle array with the Fof group number, and -1 if the particles is not in a group. *This will take some time! v=True enables a crude progress report.
    """


    def __init__(self, sim, subs=False, order=True, make_grp=None, v=False):
        self._base = weakref.ref(sim)
        self._order=order
        self._subs=subs
        if self._order is False:
            if (self.base['iord'][0] != 0 or self.base['iord'][1] != 1):
                print 'ID[0] and ID[1]:', self.base['iord'][0:2]
                raise Exception("IDs are not ordered, this won't work! Use f.halos(order=False).")
        
        self._halos = {}
        HaloCatalogue.__init__(self,sim)
        self.dtype_int = sim['iord'].dtype
        self.dtype_flt='float32' #SUBFIND data is apparently always single precision???
        self.halodir = self._name_of_catalogue(sim)
        self.header = self._readheader()
        if subs is True:
            if self.header[6]==0:
                raise Exception("This files does not contain subhalos.")
            if make_grp is True:
                raise Exception("subs=True and make_grp=True are not compatible.")
        self._tasks = self.header[4]
        self.ids = self._read_ids()
        self._keys={}
        self._halodat, self._subhalodat=self._read_groups()
        #self.data_len, self.data_off = self._read_groups()
        if make_grp:
            self.make_grp(v=v)

    def make_grp(self, name='grp', v=False):
        """
        #Creates a 'grp' array which labels each particle according to
        #its parent halo. This can take quite some time!
        Option: v=True prints out 'progress' in terms of total number of groups.
        #"""
        self.base[name] = self.get_group_array(v=v) #np.zeros(len(self.base), dtype=int)#self.get_group_array()

    def get_group_array(self, v=False):
        ar = np.zeros(len(self.base), dtype=int)-1
        for i in range(0, self.header[1]): #total number of groups
            if v:
                print "Halo #", i , "of", self.header[1]
            halo=self[i]
            ar[halo.get_index_list(self.base)] = halo._halo_id
        return ar

    def _get_halo(self,i):
        """This also works if the IDs are not sorted, but it will break the ordering by binding energy which is not desirable. We do however save the group's mostboundID"""
        if self._order is False:
            if self._subs is True:
                #this needs to be tested again on a snapshot that is not ordered!
                x = Halo(i, self, self.base, np.where(np.in1d(self.base['iord'], self.ids[self._subhalodat['sub_off'][i]:self._subhalodat['sub_off'][i]+self._subhalodat['sub_len'][i]]  )))
            else:
                x = Halo(i, self, self.base, np.where(np.in1d(self.base['iord'], self.ids[self._halodat['group_off'][i]:self._halodat['group_off'][i]+self._halodat['group_len'][i]]  )))
      
        else:
            if self._subs is False: #to use groups as halos:
                x = Halo(i, self, self.base, self.ids[self._halodat['group_off'][i]:self._halodat['group_off'][i]+self._halodat['group_len'][i]] )        
            else:
                x=Halo(i, self, self.base, self.ids[self._subhalodat['sub_off'][i]:self._subhalodat['sub_off'][i]+self._subhalodat['sub_len'][i]] )
            
        x._descriptor = "halo_"+str(i)
            
        if self._subs is False:
            for key in self._keys:
                x.properties[key]=self._halodat[key][i]
            if self.header[6]>0:
                x.properties['children']=np.where( self._subhalodat['sub_groupNr']==i )[0] #this is the FIRST level of substructure, sub-subhalos (etc) can be accessed via the subs=True output (below)
        else:
            for key in self._keys:
                x.properties[key]=self._subhalodat[key][i]
            x.properties['children']=np.where(self._subhalodat['sub_parent']==i)[0] #this goes down one level in the hierarchy, i.e. a subhalo will have all its sub-subhalos listed, but not its sub-sub-subhalos (those will be listed in each sub-subhalo)
        return x

    def _readheader(self):
        header = np.array([], dtype='int32')
        filename = self.halodir + "/subhalo_tab_" + \
            self.halodir.split("_")[-1] + ".0"
        fd = open(filename, "rb")
        # read header: this is strange but it works: there is an extra value in
        # header which we delete in the next step
        header1 = np.fromfile(fd, dtype='int32', sep="", count=8)
        header = np.delete(header1, 4)
        fd.close()
        return header  # [4]

    def _read_ids(self):
        data_ids = np.array([], dtype=self.dtype_int)
        for n in range(0, self._tasks):
            filename = self.halodir + "/subhalo_ids_" + \
                self.halodir.split("_")[-1] + "." + str(n)
            fd = open(filename, "rb")
            # for some reason there is an extra value in header which we delete
            # in the next step
            header1 = np.fromfile(fd, dtype='int32', sep="", count=7)
            header = np.delete(header1, 4)
            # TODO: include a check if both headers agree (they better)
            ids = np.fromfile(fd, dtype=self.dtype_int, sep="", count=-1)
            fd.close()
            data_ids = np.append(data_ids, ids)
        return data_ids

    def _read_groups(self):
        halodat={}
        keys_flt=['mass', 'pos', 'mmean_200', 'rmean_200', 'mcrit_200', 'rcrit_200', 'mtop_200', 'rtop_200', 'contmass']
        keys_int=['group_len', 'group_off',  'first_sub', 'Nsubs', 'cont_count', 'mostboundID']
        for key in keys_flt:
            halodat[key]=np.array([], dtype=self.dtype_flt)
        for key in keys_int:
            halodat[key]=np.array([], dtype='int32')

        subhalodat={}
        subkeys_int=['sub_len', 'sub_off', 'sub_parent', 'sub_mostboundID', 'sub_groupNr']
        subkeys_flt=['sub_pos', 'sub_vel', 'sub_CM', 'sub_mass', 'sub_spin', 'sub_veldisp', 'sub_VMax', 'sub_VMaxRad', 'sub_HalfMassRad', ]
        for key in subkeys_int:
            subhalodat[key]=np.array([], dtype='int32')
        subhalodat['sub_mostboundID']=np.array([], dtype=self.dtype_int)
        #subhalodat['sub_groupNr']=np.array([], dtype=self.dtype_int) #these are special
        for key in subkeys_flt:
            subhalodat[key]=np.array([], dtype=self.dtype_flt)

        self._keys=keys_flt+keys_int
        if self._subs is True:
            self._keys=subkeys_flt+subkeys_int

        for n in xrange(0,self._tasks):
            filename=self.halodir+"/subhalo_tab_"+self.halodir.split("_")[-1]+"." +str(n)
            fd=open(filename, "rb")
            header1=np.fromfile(fd, dtype='int32', sep="", count=8)
            header=np.delete(header1,4)
            #read groups
            if header[0]>0:
                halodat['group_len']=np.append(halodat['group_len'], np.fromfile(fd, dtype='int32', sep="", count=header[0]))
                halodat['group_off']=np.append(halodat['group_off'], np.fromfile(fd, dtype='int32', sep="", count=header[0]))
                halodat['mass']=np.append(halodat['mass'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=header[0]))
                halodat['pos']=np.append(halodat['pos'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=3*header[0]) )
                halodat['mmean_200']=np.append(halodat['mmean_200'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=header[0]))
                halodat['rmean_200']=np.append(halodat['rmean_200'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=header[0]))
                halodat['mcrit_200']=np.append(halodat['mcrit_200'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=header[0]))
                halodat['rcrit_200']=np.append(halodat['rcrit_200'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=header[0]))
                halodat['mtop_200']=np.append(halodat['mtop_200'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=header[0]))
                halodat['rtop_200']=np.append(halodat['rtop_200'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=header[0]))
                halodat['cont_count']=np.append(halodat['cont_count'], np.fromfile(fd, dtype='int32', sep="", count=header[0]))
                halodat['contmass']=np.append(halodat['contmass'],np.fromfile(fd, dtype=self.dtype_flt, sep="", count=header[0]))
                halodat['Nsubs']=np.append(halodat['Nsubs'],np.fromfile(fd, dtype='int32', sep="", count=header[0]))
                halodat['first_sub']=np.append(halodat['first_sub'],np.fromfile(fd, dtype='int32', sep="", count=header[0]))
            #read subhalos only if expected to exist from header
            if header[5]>0:
                subhalodat['sub_len']=np.append(subhalodat['sub_len'], np.fromfile(fd, dtype='int32', sep="", count=header[5]))
                subhalodat['sub_off']=np.append(subhalodat['sub_off'], np.fromfile(fd, dtype='int32', sep="", count=header[5]))
                subhalodat['sub_parent']=np.append(subhalodat['sub_parent'], np.fromfile(fd, dtype='int32', sep="", count=header[5]))
                subhalodat['sub_mass']=np.append(subhalodat['sub_mass'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=header[5]))
                subhalodat['sub_pos']=np.append(subhalodat['sub_pos'],np.fromfile(fd, dtype=self.dtype_flt, sep="", count=3*header[5]))
                subhalodat['sub_vel']=np.append(subhalodat['sub_vel'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=3*header[5]))
                subhalodat['sub_CM']=np.append(subhalodat['sub_CM'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=3*header[5]))
                subhalodat['sub_spin']=np.append(subhalodat['sub_spin'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=3*header[5]))
                subhalodat['sub_veldisp']=np.append(subhalodat['sub_veldisp'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=header[5]))
                subhalodat['sub_VMax']=np.append(subhalodat['sub_VMax'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=header[5]))
                subhalodat['sub_VMaxRad']=np.append(subhalodat['sub_VMaxRad'],np.fromfile(fd, dtype=self.dtype_flt, sep="", count=header[5]))
                subhalodat['sub_HalfMassRad']=np.append(subhalodat['sub_HalfMassRad'], np.fromfile(fd, dtype=self.dtype_flt, sep="", count=header[5]))
                subhalodat['sub_mostboundID']=np.append(subhalodat['sub_mostboundID'], np.fromfile(fd, dtype=self.dtype_int, sep="", count=header[5]))
                subhalodat['sub_groupNr']=np.append(subhalodat['sub_groupNr'], np.fromfile(fd, dtype='int32', sep="", count=header[5]))
            fd.close()

        halodat['pos']=np.reshape(halodat['pos'], (header[1],3))
        if header[6]>0:
            #some voodoo because some SubFind files may have (at least?) one extra entry which is not really a subhalo
            real_ones=np.where(halodat['first_sub']<header[6])[0]
            fake_ones=np.where(halodat['first_sub']>=header[6])[0]
            halodat['mostboundID']=np.zeros(len(halodat['Nsubs']),dtype=self.dtype_int)-1
            halodat['mostboundID'][real_ones]=subhalodat['sub_mostboundID'][halodat['first_sub'][real_ones]]  #useful for the case of unordered snapshot IDs

            subhalodat['sub_pos']=np.reshape(subhalodat['sub_pos'], (header[6],3))
            subhalodat['sub_vel']=np.reshape(subhalodat['sub_vel'], (header[6],3))
            subhalodat['sub_CM']=np.reshape(subhalodat['sub_CM'], (header[6],3))
            subhalodat['sub_spin']=np.reshape(subhalodat['sub_spin'], (header[6],3))

        return halodat, subhalodat


    @staticmethod
    def _name_of_catalogue(sim):
        # standard path for multiple snapshot files
        snapnum = os.path.basename(
            os.path.dirname(sim.filename)).split("_")[-1]
        parent_dir = os.path.dirname(os.path.dirname(sim.filename))
        dir_path=os.path.join(parent_dir,"groups_" + snapnum)

        if os.path.exists(dir_path):
            return dir_path
        # alternative path if snapshot is single file
        else:
            snapnum = os.path.basename(sim.filename).split("_")[-1]
            parent_dir = os.path.dirname(sim.filename)
            return os.path.join(parent_dir,"groups_" + snapnum)

    @property
    def base(self):
        return self._base()

    @staticmethod
    def _can_load(sim, **kwargs):
        if os.path.exists(SubfindCatalogue._name_of_catalogue(sim)):
            return True
        else:
            return False


class SubFindHDFHaloCatalogue(HaloCatalogue) :
    """
    Gadget's SubFind Halo catalogue -- used in concert with :class:`~SubFindHDFSnap`
    """


    def __init__(self, sim) :
        super(SubFindHDFHaloCatalogue,self).__init__(sim)

        if not isinstance(sim, snapshot.gadgethdf.SubFindHDFSnap):
            raise ValueError, "SubFindHDFHaloCatalogue can only work with a SubFindHDFSnap simulation"

        self.__init_halo_offset_data()
        self.__init_subhalo_relationships()
        self.__init_halo_properties()
        self.__reshape_multidimensional_properties()
        self.__reassign_properties_from_sub_to_fof()

    def __init_ignorable_keys(self):
        self.fof_ignore = map(str.strip,config_parser.get("SubfindHDF","FoF-ignore").split(","))
        self.sub_ignore = map(str.strip,config_parser.get("SubfindHDF","Sub-ignore").split(","))

        for t in self.base._family_to_group_map.values():
            # Don't add SubFind particles ever as this list is actually spherical overdensity
            self.sub_ignore.append(t[0])
            self.fof_ignore.append(t[0])

    def __init_halo_properties(self):
        self.__init_ignorable_keys()
        self._fof_properties = self.__get_property_dictionary_from_hdf('FOF')
        self._sub_properties = self.__get_property_dictionary_from_hdf('SUBFIND')


    def __get_property_dictionary_from_hdf(self, hdf_key):
        sim = self.base
        hdf0 = sim._hdf_files.get_file0_root()

        props = {}
        for property_key in hdf0[hdf_key].keys():
            if property_key not in self.fof_ignore:
                props[property_key] = np.array([])

        for h in sim._hdf_files.iterroot():
            for property_key in props.keys():
                props[property_key] = np.append(props[property_key], h[hdf_key][property_key].value)

        for property_key in props.keys():
            arr_units = sim._get_units_from_hdf_attr(hdf0[hdf_key][property_key].attrs)
            if property_key in props:
                props[property_key] = props[property_key].view(SimArray)
                props[property_key].units = arr_units
                props[property_key].sim = sim

        return props



    def __reshape_multidimensional_properties(self):
        sub_properties = self._sub_properties
        fof_properties = self._fof_properties

        for key in sub_properties.keys():
            # Test if there are no remainders, i.e. array is multiple of halo length
            # then solve for the case where this is 1, 2 or 3 dimension
            if len(sub_properties[key]) % self.nsubhalos == 0:
                ndim = len(sub_properties[key]) / self.nsubhalos
                if ndim > 1:
                    sub_properties[key] = sub_properties[key].reshape(self.nsubhalos, ndim)

            try:
                # The case fof FOF
                if len(fof_properties[key]) % self.ngroups == 0:
                    ndim = len(fof_properties[key]) / self.ngroups
                    if ndim > 1:
                        fof_properties[key] = fof_properties[key].reshape(self.ngroups, ndim)
            except KeyError:
                pass

    def __reassign_properties_from_sub_to_fof(self):
        reassign = []
        for k,v in self._sub_properties.iteritems():
            if v.shape[0]==self.ngroups:
                reassign.append(k)

        for reassign_i in reassign:
            self._fof_properties[reassign_i] = self._sub_properties[reassign_i]
            del self._sub_properties[reassign_i]


    def __init_subhalo_relationships(self):

        nsub = 0
        nfof = 0
        for h in self.base._hdf_files.iterroot():
            parent_groups = h['SUBFIND']['GrNr']
            self._subfind_halo_parent_groups[nsub:nsub + len(parent_groups)] = parent_groups
            nsub += len(parent_groups)

            first_groups = h['SUBFIND']['FirstSubOfHalo']
            self._fof_group_first_subhalo[nfof:nfof + len(first_groups)] = first_groups
            nfof += len(first_groups)

    def __init_halo_offset_data(self):

        hdf0 = self.base._hdf_files.get_file0_root()

        self._fof_group_offsets = {}
        self._fof_group_lengths = {}
        self._subfind_halo_offsets = {}
        self._subfind_halo_lengths = {}

        self.ngroups = hdf0['FOF'].attrs['Total_Number_of_groups']
        self.nsubhalos = hdf0['FOF'].attrs['Total_Number_of_subgroups']
        self._subfind_halo_parent_groups = np.empty(self.nsubhalos, dtype=int)
        self._fof_group_first_subhalo = np.empty(self.ngroups, dtype=int)
        for ptype in self.base._family_to_group_map.values():
            ptype = ptype[0]
            self._fof_group_offsets[ptype] = np.empty(self.ngroups, dtype='int64')
            self._fof_group_lengths[ptype] = np.empty(self.ngroups, dtype='int64')
            self._subfind_halo_offsets[ptype] = np.empty(self.ngroups, dtype='int64')
            self._subfind_halo_lengths[ptype] = np.empty(self.ngroups, dtype='int64')

            curr_groups = 0
            curr_subhalos = 0

            for h in self.base._hdf_files:
                # fof groups
                offset = h[ptype]['Offset']
                length = h[ptype]['Length']
                self._fof_group_offsets[ptype][curr_groups:curr_groups + len(offset)] = offset
                self._fof_group_lengths[ptype][curr_groups:curr_groups + len(offset)] = length
                curr_groups += len(offset)

                # subfind subhalos
                offset = h[ptype]['SUB_Offset']
                length = h[ptype]['SUB_Length']
                self._subfind_halo_offsets[ptype][curr_subhalos:curr_subhalos + len(offset)] = offset
                self._subfind_halo_lengths[ptype][curr_subhalos:curr_subhalos + len(offset)] = length
                curr_subhalos += len(offset)


    def _get_halo(self, i) :
        if self.base is None :
            raise RuntimeError("Parent SimSnap has been deleted")

        if i > len(self)-1 :
            raise RuntimeError("Group %d does not exist"%i)

        type_map = self.base._family_to_group_map

        # create the particle lists
        tot_len = 0
        for g_ptype in type_map.values() :
            g_ptype = g_ptype[0]
            tot_len += self._fof_group_lengths[g_ptype][i]

        plist = np.zeros(tot_len,dtype='int64')

        npart = 0
        for ptype in type_map.keys() :
            # family slice in the SubFindHDFSnap
            sl = self.base._family_slice[ptype]

            # gadget ptype
            g_ptype = type_map[ptype][0]

            # add the particle indices to the particle list
            offset = self._fof_group_offsets[g_ptype][i]
            length = self._fof_group_lengths[g_ptype][i]
            ind = np.arange(sl.start + offset, sl.start + offset + length)
            plist[npart:npart+length] = ind
            npart += length

        return SubFindFOFGroup(i, self, self.base, plist)


    def __len__(self) :
        return self.base._hdf_files[0].attrs['Total_Number_of_groups']


    @property
    def base(self):
        return self._base()



class SubFindFOFGroup(Halo) :
    """
    SubFind FOF group class
    """

    def __init__(self, group_id, *args) :
        super(SubFindFOFGroup,self).__init__(group_id, *args)

        self._subhalo_catalogue = SubFindHDFSubhaloCatalogue(group_id, self._halo_catalogue)

        self._descriptor = "fof_group_"+str(group_id)

        # load properties
        for key in self._halo_catalogue._fof_properties.keys() :
            self.properties[key] = SimArray(self._halo_catalogue._fof_properties[key][group_id],
                                            self._halo_catalogue._fof_properties[key].units)
            self.properties[key].sim = self.base


    def __getattr__(self, name):
        if name == 'sub':
            return self._subhalo_catalogue
        else :
            return super(SubFindFOFGroup,self).__getattr__(name)


class SubFindHDFSubhaloCatalogue(HaloCatalogue) :
    """
    Gadget's SubFind HDF Subhalo catalogue.

    Initialized with the parent FOF group catalogue and created
    automatically when an fof group is created
    """

    def __init__(self, group_id, group_catalogue) :
        super(SubFindHDFSubhaloCatalogue,self).__init__(group_catalogue.base)

        self._group_id = group_id
        self._group_catalogue = group_catalogue



    def __len__(self):
        if self._group_id == (len(self._group_catalogue._fof_group_first_subhalo)-1) :
            return self._group_catalogue.nsubhalos - self._group_catalogue._fof_group_first_subhalo[self._group_id]
        else:
            return (self._group_catalogue._fof_group_first_subhalo[self._group_id + 1] -
                    self._group_catalogue._fof_group_first_subhalo[self._group_id])

    def _get_halo(self, i):
        if self.base is None :
            raise RuntimeError("Parent SimSnap has been deleted")

        if i > len(self)-1 :
            raise RuntimeError("FOF group %d does not have subhalo %d"%(self._group_id, i))

        # need this to index the global offset and length arrays
        absolute_id = self._group_catalogue._fof_group_first_subhalo[self._group_id] + i

        # now form the particle IDs needed for this subhalo
        type_map = self.base._family_to_group_map

        halo_lengths = self._group_catalogue._subfind_halo_lengths
        halo_offsets = self._group_catalogue._subfind_halo_offsets

        # create the particle lists
        tot_len = 0
        for g_ptype in type_map.values() :
            g_ptype = g_ptype[0]
            tot_len += halo_lengths[g_ptype][absolute_id]

        plist = np.zeros(tot_len,dtype='int64')

        npart = 0
        for ptype in type_map.keys() :
            # family slice in the SubFindHDFSnap
            sl = self.base._family_slice[ptype]

            # gadget ptype
            g_ptype = type_map[ptype][0]

            # add the particle indices to the particle list
            offset = halo_offsets[g_ptype][absolute_id]
            length = halo_lengths[g_ptype][absolute_id]
            ind = np.arange(sl.start + offset, sl.start + offset + length)
            plist[npart:npart+length] = ind
            npart += length

        return SubFindHDFSubHalo(i, self._group_id, self, self.base, plist)


    @property
    def base(self) :
        return self._base()

class SubFindHDFSubHalo(Halo) :
    """
    SubFind subhalo class
    """

    def __init__(self,halo_id, group_id, *args) :
        super(SubFindHDFSubHalo,self).__init__(halo_id, *args)

        self._group_id = group_id
        self._descriptor = "fof_group_%d_subhalo_%d"%(group_id,halo_id)

        # need this to index the global offset and length arrays
        absolute_id = self._halo_catalogue._group_catalogue._fof_group_first_subhalo[self._group_id] + halo_id

        # load properties
        sub_props = self._halo_catalogue._group_catalogue._sub_properties
        for key in sub_props :
            self.properties[key] = SimArray(sub_props[key][absolute_id], sub_props[key].units)
            self.properties[key].sim = self.base



def _get_halo_classes():
    # AmigaGrpCatalogue MUST be scanned first, because if it exists we probably
    # want to use it, but an AHFCatalogue will probably be on-disk too.
    _halo_classes = [GrpCatalogue, AmigaGrpCatalogue, AHFCatalogue,
                     RockstarCatalogue, SubfindCatalogue, SubFindHDFHaloCatalogue,
                     RockstarIntermediateCatalogue]

    return _halo_classes