def jacobianMatrix(diffeo, resol=[1.,1.,1.], periodic=False):
    if diffeo.ndim > 4:
        print 'No jacobian in dimension larget than 3'
        return

    if diffeo.ndim == 4:
        if periodic == True:
            w = np.mgrid[0:diffeo.shape[1], 0:diffeo.shape[2], 0:diffeo.shape[3]]
            dw = diffeo-w
            for k in range(3):
                diffeo[k,:,:,:] -= np.rint(dw[k,:,:,:]/diffeo.shape[k+1])*diffeo.shape[k+1]
        grad = np.zeros([3,3,diffeo.shape[1], diffeo.shape[2], diffeo.shape[3]])
        grad[0,:,:,:,:] = gradient(np.squeeze(diffeo[0,:,:,:]), resol=resol)
        grad[1,:,:,:,:] = gradient(np.squeeze(diffeo[1,:,:,:]), resol=resol)
        grad[2,:,:,:,:] = gradient(np.squeeze(diffeo[2,:,:,:]), resol=resol)
    elif diffeo.ndim == 3:
        if periodic == True:
            w = np.mgrid[0:diffeo.shape[1], 0:diffeo.shape[2]]
            dw = diffeo-w
            for k in range(2):
                diffeo[k,:,:] -= np.rint(dw[k,:,:]/diffeo.shape[k+1])*diffeo.shape[k+1]
        grad = np.zeros([2,2,diffeo.shape[1], diffeo.shape[2]])
        grad[0,:,:,:] = gradient(np.squeeze(diffeo[0,:,:]), resol=resol)
        grad[1,:,:,:] = gradient(np.squeeze(diffeo[1,:,:]), resol=resol)
    else:
        if periodic == True:
            w = np.mgrid[0:diffeo.shape[0]]
            dw = diffeo-w
            diffeo -= np.rint(dw/diffeo.shape[0])*diffeo.shape[0]
        grad =  np.fabs(gradient(np.squeeze(diffeo)), resol=resol)
    return grad

    def solution(self):
        if self.workers == 1:
            # there is just 1 worker at the bottom, then the tree is linear
            # height moves up in powers of 2

            if self.height == 1: # just a single worker
                return 0, 1
            else:
                k = np.rint(np.log(self.height) / np.log(2))
                return k, 2 ** (k + 1) - 1

        ratio = np.log(self.height) / np.log(self.workers)
        N = self.binary_search(ratio, 1, self.workers)

        if N is None:
            raise ValueError, "No solution exists for [%d, %d]!" \
                              % (self.height, self.workers)

        k = np.rint(np.log(self.workers) / np.log(N))
        assert self.height == (N+1) ** k, \
            "This should not happen!"

        inodes = (self.workers - 1) / (N - 1)
        ipl = self.height * (N+1) - self.workers * N

        return inodes, ipl

    def ll2xy(self, lat, lon):
        if not self.swnavdisp:
            # RADAR mode:
            # Convert lat/lon to pixel x,y

            # Normal case
            if self.lon1 > self.lon0:
                x = self.width * (lon - self.lon0) / (self.lon1 - self.lon0)

            # Wrap around:
            else:
                dellon = 180. - self.lon0 + self.lon1 + 180.
                xlon = lon + (lon < 0.) * 360.
                x = (xlon - self.lon0) / dellon * self.width

            y = self.height * (self.lat1 - lat) / (self.lat1 - self.lat0)
        else:
            # NAVDISP mode:
            qdr, dist = geo.qdrdist(self.ndlat, self.ndlon, lat, lon)
            alpha = np.radians(qdr - self.ndcrs)
            base = 30. * (self.lat1 - self.lat0)
            x = dist * np.sin(alpha) / base * self.height + self.width / 2
            y = -dist * np.cos(alpha) / base * self.height + self.height / 2

        return np.rint(x), np.rint(y)

def process_coordinates_from_data(fit_parameters, frames, dt):
    """(center_x, center_y, resting_x, resting_y, extended_x, extended_y)
    finds the resting and extended position for each dot, using the data."""
    X, Y = 0, 1
    center_x = ((fit_parameters[-1][0].offset - fit_parameters[-1][0].amplitude) -
                (fit_parameters[0][0].offset + fit_parameters[0][0].amplitude)) / 2 + fit_parameters[0][0].offset
    center_y = ((fit_parameters[-1][1].offset - fit_parameters[-1][1].amplitude) -
                (fit_parameters[0][1].offset + fit_parameters[0][1].amplitude)) + fit_parameters[0][1].offset

    # resting y positions fall when y is maximized, at t = period (pi/2 - phase) / (2 pi)
    # (this is because of our choice of coordinate system, where extension is up, towards 0)
    N = len(frames)
    y_max_t = [np.arange(start=yfit.period * (pi/2-yfit.phase) / (2*pi),
                         stop=N*dt,
                         step=yfit.period) for _, yfit in fit_parameters]
    y_max_t = [a[a > 0] for a in y_max_t]

    # extended y positions fall when y is minimized, at t = period (3 pi / 2 - phase) / (2 pi)
    y_min_t = [np.arange(start=yfit.period * (3*pi/2-yfit.phase) / (2*pi),
                         stop=N*dt,
                         step=yfit.period) for _, yfit in fit_parameters]
    y_min_t = [a[a > 0] for a in y_min_t]

    y_max_i = [np.rint(yt / dt).astype(int) for yt in y_max_t]
    y_min_i = [np.rint(yt / dt).astype(int) for yt in y_min_t]

    resting_x, resting_y = [], []
    extended_x, extended_y = [], []
    n_dots = len(frames[0])
    for dot in range(n_dots):
        extended_x.append(mean([frames[i][dot].xpos for i in y_min_i[dot]]))
        extended_y.append(mean([frames[i][dot].ypos for i in y_min_i[dot]]))
        resting_x.append(mean([frames[i][dot].xpos for i in y_max_i[dot]]))
        resting_y.append(mean([frames[i][dot].ypos for i in y_max_i[dot]]))
    return (center_x, center_y, resting_x, resting_y, extended_x, extended_y)

    def makeMaskCircle(self):
        self.initMask()
        if self.parent.data is not None and self.maskingMode > 0:
            (radiusX, radiusY) = self.mask_circle.size()
            (cornerX, cornerY) = self.mask_circle.pos()
            i0, j0 = np.meshgrid(range(int(radiusY)),
                                 range(int(radiusX)), indexing='ij')
            r = np.sqrt(np.square((i0 - radiusY / 2).astype(np.float)) +
                        np.square((j0 - radiusX / 2).astype(np.float)))
            i0 = np.rint(i0[np.where(r < radiusY / 2.)] + cornerY).astype(np.int)
            j0 = np.rint(j0[np.where(r < radiusX / 2.)] + cornerX).astype(np.int)
            i01 = i0[(i0 >= 0) & (i0 < self.parent.data.shape[1]) & (j0 >= 0) & (j0 < self.parent.data.shape[0])]
            j01 = j0[(i0 >= 0) & (i0 < self.parent.data.shape[1]) & (j0 >= 0) & (j0 < self.parent.data.shape[0])]

            _mask = np.ones_like(self.parent.data)
            _mask[j01, i01] = 0
            if self.maskingMode == 1:  # masking mode
                self.userMaskAssem *= _mask
            elif self.maskingMode == 2:  # unmasking mode
                self.userMaskAssem[j01, i01] = 1
            elif self.maskingMode == 3:  # toggle mode
                self.userMaskAssem[j01, i01] = (1 - self.userMaskAssem[j01, i01])

            # update userMask
            self.userMask = self.parent.det.ndarray_from_image(self.parent.evt, self.userMaskAssem, pix_scale_size_um=None,
                                                        xy0_off_pix=None)

            self.displayMask()
            self.parent.pk.algInitDone = False
            self.parent.pk.updateClassification()
        if self.parent.args.v >= 1: print "done makeMaskCircle!!!!!!"

def symmetryError(latt, parentlatt):
    """check that the lattice obeys all symmetry operations of a parent lattice:  R.latt.inv(R) will give an integer matrix"""
    symmerr = 0.0
    for iop in range(parentlatt.nops):
        lmat = array(latt)
        #        if det(lmat) == 0:
        #            print 'Determinant zero'
        #            print lmat
        #            return symmerror
        mmat = trimSmall(dot(dot(inv(lmat), parentlatt.symops[:, :, iop]), lmat))
        #        print 'mmat', iop
        #        print trimSmall(mmat)
        operr = 0.0
        for i in range(3):
            for j in range(3):
                if abs(rint(mmat[i, j]) - mmat[i, j]) > 1.0e-4:
                    operr += abs(rint(mmat[i, j]) - mmat[i, j])
        #                    print iop, 'Symmetry failed for mmat[i,j]',mmat[i,j]
        #                    print 'Cartesian operator'
        #                    print parentlatt.symops[:,:,iop]
        #                    print 'Cartesian Lattice'
        #                    print lmat
        if operr > 2.0e-4:
            symmerr += operr
    #            print 'Noninteger operator in superlattice for operation %s, with error %s.' % (iop,str(operr))
    #    if operr < 2.0e-4:
    #        return 0.0
    #    else:
    return symmerr

 def diff(self, i):
     """Determine the reciprocal space difference between the calculated 
     (hkl) and the closest integer (hkl) of the specified peak"""
     h, k, l = self.hkl(i)
     Q = np.matrix((h, k, l)).T
     Q0 = np.matrix((np.rint(h), np.rint(k), np.rint(l))).T
     return norm(self.Bmat * (Q - Q0))

  def get_map(self):
    import numpy as np
    map_data = "{:d} 0 {:d} {:d} {:d} 0 0\n".format(self.elements.shape[0], 
                                             self.n[0], self.n[1], self.n[2])

    my_n = np.copy(self.n)
    if not self.boundaries[0] == 'P':
      my_n[1] += 1
    if not self.boundaries[1] == 'P':
      my_n[0] += 1
    if not self.boundaries[5] == 'P':
      my_n[2] += 1

    for e in range(self.map.shape[0]):
      ix = np.rint((self.elements[e,0] - self.root[0])/self.delta[0])
      iy = np.rint((self.elements[e,4] - self.root[1])/self.delta[1])
      iz = np.rint((self.elements[e,8] - self.root[2])/self.delta[2])

      map_data += "{:d} {:d} {:d} {:d} {:d} {:d} {:d} {:d} {:d}\n".format(
                   self.map[e], 
                   get_ind(ix,   iy  , iz  , my_n),
                   get_ind(ix+1, iy  , iz  , my_n),
                   get_ind(ix,   iy+1, iz  , my_n),
                   get_ind(ix+1, iy+1, iz  , my_n),
                   get_ind(ix,   iy  , iz+1, my_n),
                   get_ind(ix+1, iy  , iz+1, my_n),
                   get_ind(ix,   iy+1, iz+1, my_n),
                   get_ind(ix+1, iy+1, iz+1, my_n))
     
    return map_data

 def __call__(self, t):
     if self.clock is None:
         raise ValueError('Can only call timed arrays if they are based on a clock.')
     else:
         if isinstance(t, (list, tuple)):
             t = numpy.array(t)
         if isinstance(t, neurongroup.TArray):
             # In this case, we know that t = ones(N)*t so we just use the first value
             t = t[0]
         elif isinstance(t, numpy.ndarray):
             if len(self.shape) > 2:
                 raise ValueError('Calling TimedArray with array valued t only supported for 1D or 2D TimedArray.')
             if len(self.shape) == 2 and len(t) != self.shape[1]:
                 raise ValueError('Calling TimedArray with array valued t on 2D TimedArray requires len(t)=arr.shape[1]')
             t = numpy.array(numpy.rint((t - self._t_init) / self._dt), dtype=int)
             t[t < 0] = 0
             t[t >= len(self.times)] = len(self.times) - 1
             if len(self.shape) == 1:
                 return numpy.asarray(self)[t]
             return numpy.asarray(self)[t, numpy.arange(len(t))]
         t = float(t)
         ot = t
         t = int(numpy.rint((t - self._t_init) / self._dt))
         if t < 0: t = 0
         if t >= len(self.times): t = len(self.times) - 1
         return numpy.asarray(self)[t]

def apply_filter(image, img_filter, horizontally=True):
    filter_size = img_filter.shape
    filter_size = filter_size[0]

    pad = filter_size // 2
    height, width = image.shape

    new_image = np.zeros((height, width), dtype=np.uint8)

    if horizontally:
        for h in range(pad, height - pad):
            for w in range(pad, width - pad):
                new_image[h, w] = np.clip(np.rint(
                    np.sum(np.multiply(image[h - pad: h + pad + 1, w - pad: w + pad + 1], img_filter))),
                    MIN_PIXEL_VALUE, MAX_PIXEL_VALUE)

                aa = np.sum(np.multiply(image[h - pad: h + pad + 1, w - pad: w + pad + 1], img_filter))
    else:
        for w in range(pad, width - pad):
            for h in range(pad, height - pad):
                new_image[h, w] = np.clip(np.rint(
                    np.sum(np.multiply(image[h - pad: h + pad + 1, w - pad: w + pad + 1], img_filter))),
                    MIN_PIXEL_VALUE, MAX_PIXEL_VALUE)

    return new_image

    def createFromRandomVector(vector, alphaDivisor, randomVector, alphaStage):
        if alphaDivisor is None:
            new = Solution(vector + numpy.rint(randomVector).astype(int))
        else:
            # Calculate new requests component
            newVector = numpy.zeros(vector.shape, dtype=object)
            if alphaStage == 0:
                newVector[0] = vector[0] \
                    + (Solution.sizeDomainRequestComponent * (alphaDivisor.numerator * int(round(randomVector[0] * Solution.randomPrecisionMult)))) \
                        // (alphaDivisor.denominator * Solution.randomPrecisionMult)
            else:
                newVector[0] = vector[0] + int(round(randomVector[0]))

            candidate = Solution(newVector)

            # Set the rest of the vector if the first component is valid
            if alphaStage != 0 and candidate.canApplyTimeAdjust():
                newVector[1:] = vector[1:] \
                    + (candidate.getSizeDomainEachBus() * (alphaDivisor.numerator * numpy.rint(randomVector[1:] * Solution.randomPrecisionMult).astype(int).astype(object))) \
                        // (alphaDivisor.denominator * Solution.randomPrecisionMult)
            else:
                newVector[1:] = vector[1:]
            new = Solution(newVector)

        return new

def correlation(VN,r,maxs,dmaxs,bins,N,ld):
    R = np.zeros(shape=(N,N))
    for i in range(N):
        for j in range(i):
            dx=r[i,0]-r[j,0]
            dy=r[i,1]-r[j,1]
            dz=r[i,2]-r[j,2]
            dx -= np.rint(dx / ld) * ld
            dy -= np.rint(dy / ld) * ld
            dz -= np.rint(dz / ld) * ld        
            drr = math.sqrt( dx * dx + dy * dy + dz * dz)
            R[i,j] = drr
            R[j,i] = drr

    #bins = 30
    #maxs = np.linspace(0,ld,num=bins)
    #dmaxs = maxs[1] - maxs[0]
    g = np.zeros(shape=(bins,))
    for i in range(1,bins):
        maxx = maxs[i]
        for j in range(N):
            for k in range(j):
                dist = 1/(4*math.pi*R[j,k]*R[j,k]*dmaxs)
                x = R[j,k] - maxx
                if (-1*dmaxs < (x) <0):
                    g[i-1] += 2*dist
            
    g = g*VN
    return g

    def configure(self, bin_width_s, record_length_s, number_of_gates=0):
        """ Configuration of the fast counter.

        @param float bin_width_s: Length of a single time bin in the time trace
                                  histogram in seconds.
        @param float record_length_s: Total length of the timetrace/each single
                                      gate in seconds.
        @param int number_of_gates: optional, number of gates in the pulse
                                    sequence. Ignore for not gated counter.

        @return tuple(binwidth_s, gate_length_s, number_of_gates):
                    binwidth_s: float the actual set binwidth in seconds
                    gate_length_s: the actual set gate length in seconds
                    number_of_gates: the number of gated, which are accepted
        """
        # Do nothing if fast counter is running
        if self.statusvar >= 2:
            binwidth_s = self._binwidth / self._internal_clock_hz
            gate_length_s = self._gate_length_bins * binwidth_s
            return binwidth_s, gate_length_s, self._number_of_gates

        # set class variables
        self._binwidth = int(np.rint(bin_width_s * self._internal_clock_hz))

        # calculate the actual binwidth depending on the internal clock:
        binwidth_s = self._binwidth / self._internal_clock_hz

        self._gate_length_bins = int(np.rint(record_length_s / bin_width_s))
        gate_length_s = self._gate_length_bins * binwidth_s

        self._number_of_gates = number_of_gates

        self.statusvar = 1
        return binwidth_s, gate_length_s, number_of_gates

def main():
    """Main Function"""
    # Read dat file and save pinned pixels
    if os.path.exists( datFile ):
        # Make output directory
        if not os.path.exists( resultsDir ):
            os.system( 'mkdir %s' % resultsDir )
            
        img = ReadDat ( datFile )
        np.array( img.pinned , dtype=bool ).tofile('%s/pinned.dat' % resultsDir )
        
        # Save details of image size and miniblock sizes
        f = open( '%s/array.log' % resultsDir , 'w' )
        f.write ( '%s\t%s\t%s\t%s' % ( img.rows , img.cols , img.miniR , img.miniC ) )
        f.close ( )
        
        # Empty chip beadfind analysis
        bf  = BeadFind ( img )
        np.array( bf['actpix'] , dtype=bool ).tofile('%s/actpix.dat' % resultsDir )
        
        # Note that these are only useful in the miniblock sizes 
        # 12 x 14 blocks (miniR = 111, miniC = 92)
        # For thumbnail, it's different, of course.  8 x 12 blocks of (miniR = 100 , miniC = 100)
        np.array( np.rint( bf['bfmat'] ) , dtype=np.dtype('i2') ).tofile('%s/ebfvals.dat' % resultsDir )
        np.array( np.rint( 10000 * bf['gains'] ) , dtype=np.dtype('i2') ).tofile('%s/gaincorr.dat' % resultsDir )
        
        # Buffering analysis
        bt  = BufferTest( img )
        
        np.array( np.rint( bt['slopes'] ) , dtype=np.dtype('i2') ).tofile('%s/slopes.dat' % resultsDir )
        np.array( bt['t0'] , dtype=np.dtype('i1') ).tofile('%s/t0.dat' % resultsDir )
        
    else:
        print('Error!  Acquisition file not found.  Please do not skip the calibration before loading.')

 def __init__(self, n=None, azim_max=np.pi/2,
              diameter=speed_of_sound_in_air*650*usecond,
              itd=None, samplerate=None, fractional_itds=False):
     if itd is None:
         azim = np.linspace(-azim_max, azim_max, n)
         itd = diameter*np.sin(azim)/speed_of_sound_in_air
         coords = make_coordinates(azim=azim, itd=itd)
     else:
         coords = make_coordinates(itd=itd)
     self.itd = itd
     samplerate = self.samplerate = get_samplerate(samplerate)
     if not fractional_itds:
         dl = itd.copy()
         dr = -itd
         dl[dl<0] = 0
         dr[dr<0] = 0
         dl = np.array(np.rint(dl*samplerate), dtype=int)
         dr = np.array(np.rint(dr*samplerate), dtype=int)
         idxmax = max(np.amax(dl), np.amax(dr))
         data = np.zeros((2, len(itd), idxmax+1))
         data[0, np.arange(len(itd)), dl] = 1
         data[1, np.arange(len(itd)), dr] = 1
     else:
         delays = np.hstack((itd/2, -itd/2))
         fd = FractionalDelay(silence(1*ms, samplerate=samplerate), delays)
         ir = fd.impulse_response
         data = np.zeros((2, len(itd), fd.filter_length))
         data[0, :, :] = ir[:len(itd), :]
         data[1, :, :] = ir[len(itd):, :]
         self.delay_offset = fd.delay_offset
     self.hrtfset = HRTFSet(data, samplerate, coords)
     self.hrtfset.name = 'ITDDatabaseSubject'
     self.subjects = ['0']

    def output(self):
        """Return the drawn line and the resulting scan.

        Returns
        -------
        line_image : (M, N) uint8 array, same shape as image
            An array of 0s with the scanned line set to 255.
            If the linewidth of the line tool is greater than 1,
            sets the values within the profiled polygon to 128.
        scan : (P,) or (P, 3) array of int or float
            The line scan values across the image.
        """
        end_points = self.line_tool.end_points
        line_image = np.zeros(self.image_viewer.image.shape[:2],
                              np.uint8)
        width = self.line_tool.linewidth
        if width > 1:
            rp, cp = measure.profile._line_profile_coordinates(
                *end_points[:, ::-1], linewidth=width)
            # the points are aliased, so create a polygon using the corners
            yp = np.rint(rp[[0, 0, -1, -1],[0, -1, -1, 0]]).astype(int)
            xp = np.rint(cp[[0, 0, -1, -1],[0, -1, -1, 0]]).astype(int)
            rp, cp = draw.polygon(yp, xp, line_image.shape)
            line_image[rp, cp] = 128
        (x1, y1), (x2, y2) = end_points.astype(int)
        rr, cc = draw.line(y1, x1, y2, x2)
        line_image[rr, cc] = 255
        return line_image, self.scan_data

def value_diff(I, mu):
    '''we pass in a numpy array: with the R, G, B components being the innermost'''
    if len(mu.shape) == 1:
        mu = np.rint(mu[np.newaxis,np.newaxis,:]).astype(np.uint8)
    if len(I.shape) == 1:
        I = np.rint(I[np.newaxis,np.newaxis,:]).astype(np.uint8)
    return(np.squeeze(cv2.subtract(cv2.cvtColor(I, cv2.COLOR_BGR2HSV)[...,2], cv2.cvtColor(mu, cv2.COLOR_BGR2HSV)[...,2], dtype=cv2.CV_32S)))

def writeImageFile(xrayImageFilename,BMDprojected,imageSize,imageFormat='bmp',smooth=True):
    """Create an image from array and write to file"""
    
    # Re-import PIL.Image and PIL.ImageFilter into current function namespace
    from PIL import Image, ImageFilter
    
    # Convert to 8-bit
    xray = BMDprojected.copy()
    xray = np.asarray(xray,dtype=np.int8)

    # Create image from array.
    xray = xray[:,::-1]       
    xrayImage = Image.fromarray(xray.transpose(),mode='L')
    
    # Resize image
    xsize,ysize = xrayImage.size    
    bigSide = np.argmax(xrayImage.size)
    if bigSide==0: scale = float(imageSize)/xsize
    else:          scale = float(imageSize)/ysize
    xsize = int(np.rint(scale*xsize))
    ysize = int(np.rint(scale*ysize))  
    xrayImage = xrayImage.resize((xsize,ysize),Image.BILINEAR)
    if smooth: xrayImage = xrayImage.filter(ImageFilter.SMOOTH)
    
    # Save xray image to file
    if xrayImageFilename.split('.')[-1]!=imageFormat: xrayImageFilename+='.%s' % imageFormat
    xrayImage.save(xrayImageFilename,imageFormat) 

    return

def AnalyzeSamind(str_f_samind, str_f_sg, str_of_iost, str_of_iosg, int_cutoff=0):
    dfm_samind = pd.read_csv(str_f_samind, sep='\t', header=None)
    dfm_samind.columns = ['readID', 'batchID', 'chrom', 'sbeg', 'send', 'sgID', 'sgstrand', 'sgbeg', 'sgend', 'sgseq', 'c_site', 'CIGAR', 'idtype', 'idbeg', 'idend', 'idlen', 'fm_status', 'factor', 'count', 'freq']

    dfm_sg = pd.read_csv(str_f_sg, header=None, index_col=None, sep='\t')
    dfm_sg.columns = ['sgID', 'chrom', 'strand', 'sbeg', 'send', 'qseq', 'c_site']

    func_samind = lambda x: pd.DataFrame(dict(batchID=x.batchID.unique(),
                                           sgID=x.sgID.unique(),
                                           allrc = np.int(np.rint(x.freq.sum())),
                                           none=np.int(np.rint(x.freq[x.fm_status=='None'].sum())),
                                              r_ind = x.freq[x.fm_status!='None'].sum()/x.freq.sum() if x.freq.sum()!=0 else 0,
                                           inf=np.int(np.rint(x.freq[x.fm_status=='INF'].sum())),
                                           otf=np.int(np.rint(x.freq[x.fm_status=='OTF'].sum())),
                                           r_reffect=x.freq[x.fm_status=='OTF'].sum()/x.freq[x.fm_status!='None'].sum() if x.freq[x.fm_status!='None'].sum()!=0 else 0,
                                              r_effect=x.freq[x.fm_status=='OTF'].sum()/x.freq.sum() if x.freq.sum()!=0 else 0))
    
    dfm_iost = dfm_samind.groupby(['batchID', 'sgID'], group_keys=False).apply(func_samind).reset_index(drop=True).reindex(columns=['batchID', 'sgID', 'allrc', 'none', 'r_ind', 'inf', 'otf', 'r_reffect', 'r_effect'])

    func_io = lambda x: pd.DataFrame(dict(sgID=x.sgID.unique(),
                                          allrc=x.allrc.sum(),
                                          none=x.none.sum(),
                                          inf=x.inf.sum(),
                                          otf=x.otf.sum(),
                                          r_ind=(x.inf.sum()+x.otf.sum())/x.allrc.sum() if x.allrc.sum()!=0 else 0,
                                          r_reffect=x.otf.sum()/(x.inf.sum()+x.otf.sum()) if x.inf.sum()+x.otf.sum()!=0 else 0,
                                          r_effect=x.otf.sum()/x.allrc.sum() if x.allrc.sum()!=0 else 0))
    dfm_iost = dfm_iost.groupby(['sgID'], group_keys=False).apply(func_io).reset_index(drop=True).reindex(columns=['sgID', 'allrc', 'none', 'inf', 'otf', 'r_ind', 'r_reffect', 'r_effect'])
    dfm_iost = dfm_iost.ix[dfm_iost.allrc>=int_cutoff, ]
    dfm_iost.to_csv(str_of_iost, sep='\t', index=None)

    dfm_sg = pd.merge(dfm_sg, dfm_iost[['sgID']], on='sgID')
    dfm_sg.to_csv(str_of_iosg, sep='\t', header=None, index=None)

 def angle_diff(self, i):
     """Determine the polar angle difference between the calculated 
     (hkl) and the closest integer (hkl) of the specified peak"""
     h, k, l = self.hkl(i)
     (h0, k0, l0) = (np.rint(h), np.rint(k), np.rint(l))
     polar0 = 2 * np.arcsin(self.unitcell.ds((h0,k0,l0))*self.wavelength/2)
     return np.abs(self.polar(i) - polar0)