def calcZ01andZ10(Y, MPS):
    try:
        U, S, V = spla.svd(Y, full_matrices=True)
    except spla.LinAlgError as err:
        if 'empty' in err.message:
            row, col = Y.shape
            Z01 = np.array([], dtype=Y.dtype).reshape(row, 0)
            Z10 = np.array([], dtype=Y.dtype).reshape(0, col)
            print "Empty", Z01.shape, Z10.shape
        else:
            print >> sys.stderr, "calcZ01andZ10: Error", I, err
            raise
    else:
        print "S", S, "\nU", U, "\nV", V
        __, chi, __ = MPS.shape
        mask = (S > expS) #np.array([True] * S.shape[0])
        mask[xiTilde - chi:] = False
        U = np.compress(mask, U, 1)
        S = np.compress(mask, S, 0)
        V = np.compress(mask, V, 0)

        Ssq = np.diag(np.sqrt(S))
        Z01 = np.dot(U, Ssq)
        Z10 = np.dot(Ssq, V)
        print "Fill ", U.shape, V.shape, "mask", mask

    eps = np.linalg.norm(np.dot(Z01, Z10))
    print "eps", I, eps
    print "Z01", Z01.shape, "\n", Z01, "\nZ10", Z10.shape, "\n", Z10

    return Z01, Z10

    def create_projection_as_numeric_array_3D(self, attr_indices, **settings_dict):
        valid_data = settings_dict.get("valid_data")
        class_list = settings_dict.get("class_list")
        jitter_size = settings_dict.get("jitter_size", 0.0)

        if valid_data == None:
            valid_data = self.get_valid_list(attr_indices)
        if sum(valid_data) == 0:
            return None

        if class_list == None and self.data_has_class:
            class_list = self.original_data[self.data_class_index]

        xarray = self.no_jittering_scaled_data[attr_indices[0]]
        yarray = self.no_jittering_scaled_data[attr_indices[1]]
        zarray = self.no_jittering_scaled_data[attr_indices[2]]
        if jitter_size > 0.0:
            xarray += (np.random.random(len(xarray))-0.5)*jitter_size
            yarray += (np.random.random(len(yarray))-0.5)*jitter_size
            zarray += (np.random.random(len(zarray))-0.5)*jitter_size
        if class_list != None:
            data = np.compress(valid_data, np.array((xarray, yarray, zarray, class_list)), axis = 1)
        else:
            data = np.compress(valid_data, np.array((xarray, yarray, zarray)), axis = 1)
        data = np.transpose(data)
        return data

def calculate_realexptime(id_arr, utc_arr, dsec_arr, diff_arr, req_texp, utc_list):
    """Calculates the real exposure time.
    This makes the following assumptions:
    #. That the measurement after the turn of the second is a fiducial
    #. That there is an integer number of frames between each fiducial exposure
    #. We then set up a metric which is Y=np.sum(i-int(i)) where i=dt/t_exp
    #. Then the minimum of Y is found between the requested exposure time and the median time difference
    #. And the best exposure time is the time at that minimum

    returns median exposure time and real exposure time
    """
    t_exp=0

    # calculate the median time
    try:
        t_wrong=np.median(diff_arr)
    except:
        raise SaltError('Unable to calculate median time difference')

    # Compress the arrays to find those closest to the second mark
    mask=(dsec_arr<t_wrong)*(diff_arr>0)
    t=np.compress(mask,utc_arr)
    s=np.compress(mask,dsec_arr)
    id=np.compress(mask,id_arr)

    # Now set up the components in the equation
    try:
        t_start=t[0]
        dt=t[1:]-t[0]
    except Exception, e:
        msg='Unable to set up necessary arrays because %s' % e
        raise SaltError(msg)

def test8():
    global L0, N

    L = deepcopy(L0)
    rho = zeros(N, 'double')

    rho[random.sample(xrange(N), N/2)] = 1

    print rho

    LI = linalg.inv(L)
    #print L
    #print LI
    #I = numpy.dot(L,LI)
    #I[abs(I)<0.001] = 0
    #print I

    t = numpy.greater(rho, 0)
    X = numpy.zeros((N,N))
    for i in xrange(N):
        X[0][i] = i
    print X
    LIC = numpy.compress(t, LI, 1)
    print LIC
    LIC = numpy.compress(t, LIC, 0)
    print LIC
    LICI = linalg.inv(LIC)
    print LICI

 def myzpk2tf(self, z, p, k):
         z = np.atleast_1d(z)
         k = np.atleast_1d(k)
         if len(z.shape) > 1:
                 temp = np.poly(z[0])
                 b = np.zeros((z.shape[0], z.shape[1] + 1), temp.dtype.char)
                 if len(k) == 1:
                         k = [k[0]] * z.shape[0]
                 for i in range(z.shape[0]):
                         b[i] = k[i] * poly(z[i])
         else:
                 b = k * np.poly(z)
         a = np.atleast_1d(np.poly(p))
         # Use real output if possible. Copied from numpy.poly, since
         # we can't depend on a specific version of numpy.
         if issubclass(b.dtype.type, np.complexfloating):
                 # if complex roots are all complex conjugates, the roots are real.
                 roots = np.asarray(z, complex)
                 pos_roots = np.compress(roots.imag > 0, roots)
                 neg_roots = np.conjugate(np.compress(roots.imag < 0, roots))
                 if len(pos_roots) == len(neg_roots):
                         if np.all(np.sort_complex(neg_roots) == np.sort_complex(pos_roots)):
                                 b = b.real.copy()
         if issubclass(a.dtype.type, np.complexfloating):
                 # if complex roots are all complex conjugates, the roots are real.
                 roots = np.asarray(p, complex)
                 pos_roots = np.compress(roots.imag > 0, roots)
                 neg_roots = np.conjugate(np.compress(roots.imag < 0, roots))
                 if len(pos_roots) == len(neg_roots):
                         if np.all(np.sort_complex(neg_roots) == np.sort_complex(pos_roots)):
                                 a = a.real.copy()
         return b, a

 def estimateState(self):
     """ Updates the estimate of the state """
     best = numpy.argmax(self.Weights)
     beststate = self.States[best,:]
     
     #print "Best State:", beststate
     
     cond = (numpy.sum(numpy.fabs(self.States - beststate), axis=1) < 1)
     beststates = numpy.compress(cond, self.States, axis=0)
     bestweights = numpy.compress(cond, self.Weights)
     
     #print "States", self.States
     #print "States within window:", cond
     #print "States close to best", len(beststates), beststates
     #print "Weights close to best", bestweights
     
     #print "Product:", (bestweights*beststates.T).T
     bestweights /= numpy.sum(bestweights)
     self.State = numpy.sum((bestweights*beststates.T).T, axis=0)
     #print "Estimate:", self.State
     
     #print numpy.fabs(numpy.arctan2(self.State[Localisation.YDOT], self.State[Localisation.XDOT]) - self.State[Localisation.THETA]) -  self.__controlToVelocityVector()
     
     if numpy.isnan(self.State[0]):
         print "FAIL"
     self.__updateAttributesFromState()

def unpack_data(path, delimiter, filtr=False, split_column=-1):
    """Measurements and errors are assumed to be alternating. The last
    pair of columns corresponds to the dependent variable
    while the preceeding are independent.

    If filtr is True, values larger than the error are removed.

    If split_column is given, the data is split into lumps with a column
    value in that column, e.g if split_column=(n-1) [n.b we count from 0] and
    the nth column contains trial number, chemical type etc. this value will
    be used to categorise the rest of the data and the other procedures
    will run sequentially on each category, as if they were in different files."""

    raw = np.loadtxt(path, delimiter=delimiter, skiprows=1)
    data_name = os.path.splitext(os.path.basename(path))[0]

    if split_column != -1:
        raws = split_file(raw, split_column, data_name)
    else:
        # Needed to generalise following iterative step.
        raws = [(data_name, raw)]
    for (name, raw) in raws:
        meas = raw[:, ::2].transpose()
        err = raw[:, 1::2].transpose()
        if filtr:
            test = (abs(meas) >= err).prod(axis=0)
            meas = np.compress(test, meas, axis=1)
            err = np.compress(test, err, axis=1)

        if meas.shape[0] == 2:
            A = (meas[:-1].ravel(), err[:-1].ravel())
            yield name, (A, (meas[-1], err[-1]))
        else:
            yield name, ((meas[:-1], err[:-1]), (meas[-1], err[-1]))

    def utest( self, score ):
        """
        Gives the Mann-Withney U test probability that the score is
        random.  See:

        Mason & Graham (2002) Areas beneath the relative operating
        characteristics (ROC) and relative operating levels (ROL)
        curves: Statistical significance and interpretation

        Note (1): P-values below ~1e-16 are reported as 0.0.
        See zprob() in Biskit.Statistics.stats!

        Note (2): the P-value does not distinguish between positive
        and negative deviations from random -- a ROC area of 0.1 will
        get the same P-value as a ROC area of 0.9.

        @param score: the score predicted for each item
        @type  score: [ float ]

        @return: 1-tailed P-value
        @rtype: float
        """
        sample1 = N.compress( self.positives, score )
        sample1 = sample1[-1::-1]  # invert order

        sample2 = N.compress( N.logical_not( self.positives ), score )
        sample2 = sample2[-1::-1]  # invert order

        sample1 = sample1.tolist()
        sample2 = sample2.tolist()

        p = stats.mannwhitneyu( sample1, sample2 )
        return p[1]

def whiskers_and_fliers(x, q1, q3, transformout=None):
    wnf = {}
    if transformout is None:
        transformout = lambda x: x

    iqr = q3 - q1
    # get low extreme
    loval = q1 - (1.5 * iqr)
    whislo = np.compress(x >= loval, x)
    if len(whislo) == 0 or np.min(whislo) > q1:
        whislo = q1
    else:
        whislo = np.min(whislo)

    # get high extreme
    hival = q3 + (1.5 * iqr)
    whishi = np.compress(x <= hival, x)
    if len(whishi) == 0 or np.max(whishi) < q3:
        whishi = q3
    else:
        whishi = np.max(whishi)

    wnf['fliers'] = np.hstack([
        transformout(np.compress(x < whislo, x)),
        transformout(np.compress(x > whishi, x))
    ])
    wnf['whishi'] = transformout(whishi)
    wnf['whislo'] = transformout(whislo)

    return wnf

def doit(input_file, output_file, regularization, wants_normalization):
	# Read the entire file.
	data= tuple(tuple(map(float, line.split(','))) for line in input_file)
	if len(data) == 0:
		print("no data", file=sys.stderr)
		return

	# Create X and Y indices.  Assume the last column contains the output and
	# the rest contain the inputs.
	y_index= len(data[0]) - 1
	x_indices= tuple(range(y_index))

	# Create and print the model parameters, normalizing the data if requested.
	data= np.array(data)
	x= np.compress(as_bools(x_indices), data, 1)
	mu= list(it.repeat(0.0, x.shape[1]))
	sigma= list(it.repeat(1.0, x.shape[1]))
	if wants_normalization:
		for i in range(x.shape[1]):
			mu[i]= np.mean(x[:,i])
			sigma[i]= np.std(x[:,i])
			if sigma[i] == 0.0:
				sigma[i]= 1.0
			x[:,i]= (x[:,i] - mu[i]) / sigma[i]
	y= np.compress(as_bools(y_index), data, 1).squeeze()
	model= MinimizedModel(x, y, regularization, mu, sigma)
	print(model, file=output_file)

 def setFlaggedImageRange(self):
   (nx,ny) = self.raw_image.shape
   num_elements = nx * ny
   if self._flags_array is None:
     if not self._nan_flags_array is None:
       flags_array = self._nan_flags_array.copy()
     else:
       flags_array = numpy.zeros((nx,ny),int);
   else:
     flags_array = self._flags_array.copy()
     if not self._nan_flags_array is None:
       flags_array = flags_array + self._nan_flags_array
   flattened_flags = numpy.reshape(flags_array,(num_elements,))
   if self.raw_image.dtype == numpy.complex64 or self.raw_image.dtype == numpy.complex128:
     real_array =  self.raw_image.real
     imag_array =  self.raw_image.imag
     flattened_real_array = numpy.reshape(real_array.copy(),(num_elements,))
     flattened_imag_array = numpy.reshape(imag_array.copy(),(num_elements,))
     real_flagged_array = numpy.compress(flattened_flags == 0, flattened_real_array)
     imag_flagged_array = numpy.compress(flattened_flags == 0, flattened_imag_array)
     flagged_image = numpy.zeros(shape=real_flagged_array.shape,dtype=self.raw_image.dtype)
     flagged_image.real = real_flagged_array
     flagged_image.imag = imag_flagged_array
   else:
     flattened_array = numpy.reshape(self.raw_image.copy(),(num_elements,))
     flagged_image = numpy.compress(flattened_flags == 0, flattened_array)
   self.setImageRange(flagged_image)

def fit_gauss_to_hist(binheights,
                     binedges,
                     binerrors,
                     fitmin=0,
                     fitmax=None,
                     p0=None,  # guesses for norm, mu, sigma
                     fitcolor="r"):

    left_binedges = binedges[:-1]
    
    if fitmax == None:
        fitmax = np.max(binedges)
    
    # cut data to values needed for fitting
    cut_mask = (left_binedges>fitmin)*(left_binedges<fitmax)
    fitx = np.compress(cut_mask, left_binedges)
    fity = np.compress(cut_mask, binheights)
    cut_binerrors = np.compress(cut_mask, binerrors)
    
    # p0 = [200.,meanguess,10.] 
    popt, pcov = scipy.optimize.curve_fit(gauss, fitx, fity,
                                          sigma=cut_binerrors,
                                          absolute_sigma=True,
                                          p0=p0)
    perr = np.sqrt(np.diag(pcov))
    # draw fitfunction
    xbase = np.linspace(fitmin,fitmax,1000)
    plt.plot(xbase, gauss(xbase, popt[0], popt[1], popt[2]),
             color=fitcolor,
             linewidth=2.) 
    print "optimized norm, mu, sigma:\n", popt
    print "corresponding errors\n", np.sqrt(np.diag(pcov))
    print "corresponding covariance matrix:\n", pcov
    return popt, pcov

 def _addChildren(self, parent_node_num, I, cur_depth, right_mask, left_mask, y):
     """Modifies self.nodes_dict, self.stack
     """
     # do the right branch
     r_tmp = numpy.compress(right_mask, y)
     if (r_tmp.shape[0] > 0): 
         # then there is a reason to add a right child
         r_node_num = self.num_nodes
         r_child = TreeNode()
         r_child.parent = parent_node_num
         r_child.constval = numpy.average(r_tmp)
         self.nodes_dict[parent_node_num].Rchild = r_node_num
         self.nodes_dict[r_node_num] = r_child
         self.stack.append( self.StackEntry(r_node_num, cur_depth+1,\
                                            numpy.compress(right_mask, I)))
         self.num_nodes += 1
    
     # do the left branch
     l_tmp = numpy.compress(left_mask, y)
     if (l_tmp.shape[0] > 0): 
         l_node_num = self.num_nodes
         l_child = TreeNode()
         l_child.parent = parent_node_num
         l_child.constval = numpy.average(l_tmp)
         self.nodes_dict[parent_node_num].Lchild = l_node_num
         self.nodes_dict[l_node_num] = l_child
         self.stack.append( self.StackEntry(l_node_num, cur_depth+1,\
                                            numpy.compress(left_mask, I)))
         self.num_nodes += 1

def getMaxPoints(arr):
    # [TODO] Work out for RGB rather than array, and maybe we don't need the filter, but hopefully speeds it up.
    # Reference http://scipy-cookbook.readthedocs.io/items/FiltFilt.html
    arra = filtfilt(b,a,arr)
    maxp = maxpoints(arra, order=(len(arra)/20), mode='wrap')
    minp = minpoints(arra, order=(len(arra)/20), mode='wrap')

    points = []

    for i in range(3):
        mas = np.equal(np.greater_equal(maxp,(i*(len(arra)/3))), np.less_equal(maxp,((i+1)*len(arra)/3)))
        k = np.compress(mas[0], maxp)
        if len(k)==0:
            continue
        points.append(sum(k)/len(k))

    if len(points) == 1:
        return points, []

    points = np.compress(np.greater_equal(arra[points],(max(arra)-min(arra))*0.40 + min(arra)),points)
    rifts = []
    for i in range(len(points)-1):
        mas = np.equal(np.greater_equal(minp, points[i]),np.less_equal(minp,points[i+1]))
        k = np.compress(mas[0], minp)
        rifts.append(k[arra[k].argmin()])

    return points, rifts

def gammaGunFilter(data, quiet=True):
    '''Filters gamma gun data set to clean it up some, removing some of the crap.'''
    if not quiet: print "Filtering"   
    data     = np.compress([len(event[2]) >= 2 for event in data], data, axis=0)                    #|Get at least two ROI's    
    data     = np.compress([np.max(event[2]['eta'])*np.min(event[2]['eta']) < 0 for event in data], 
                           data, axis=0)                                                            #|Require an eta separation of opposite endcaps to prevent high errors
    return data

 def _findRobots(self):
     """ Finds the robots amoung the edges found
     """
     ## for each right edge find the next closest left edge. This forms an edge pair that could be robot 
     self.Robots = list()
     if len(self.RightEdges) == 0 or len(self.LeftEdges) == 0:
         return
         
     for rightedge in self.RightEdges:
         leftedge = self.LeftEdges[0]
         i = 1
         while leftedge < rightedge:
             if i >= len(self.LeftEdges):
                 break
             leftedge = self.LeftEdges[i]
             i = i + 1
             
         ## now calculate the distance between the two edges
         distance = self.__calculateDistanceBetweenEdges(leftedge, rightedge)
         
         if distance > self.MINIMUM_NAO_WIDTH and distance < self.MAXIMUM_NAO_WIDTH:
             x = self.CartesianData[0,rightedge:leftedge+1]
             y = self.CartesianData[1,rightedge:leftedge+1]
             r = self.PolarData[0,rightedge:leftedge+1]
             c = numpy.less(r, 409.5)
             x = numpy.compress(c, x)
             y = numpy.compress(c, y)                
             robotx = self.__averageObjectDistance(x)
             roboty = self.__averageObjectDistance(y)
             c = numpy.logical_and(numpy.less(numpy.fabs(x - robotx), self.MAXIMUM_NAO_WIDTH), numpy.less(numpy.fabs(y - roboty), self.MAXIMUM_NAO_WIDTH))
             x = numpy.compress(c, x)
             y = numpy.compress(c, y)
             robotr = math.sqrt(robotx**2 + roboty**2)
             robotbearing = math.atan2(roboty, robotx)
             self.Robots.append(Robot(robotx, roboty, robotr, robotbearing, x, y))

def roiEnergyAnalysis(data):
    '''Troubleshooting function, compares the observed sum energy in an ROI to the genEnergy''' 
    genEnergies = [] 
    sumEnergies = []
    pbar = progressbar("Processing event &count&:", len(data)+1)
    pbar.start()
    count = 0
    for event in data:
        genEnergy = event[2]['getpt'] * np.cosh(event[2]['geneta'])         
        for i in range(len(genEnergy)):
            clustersIndices = np.compress(event[1]['ROI'] == i, event[1]['clusterID'], axis=0)      #|Only take clusters corresponding to right ROI
            clusterEnergies = []
            for clusterID in clustersIndices:                                                       #|Only take hits corresponding to correct cluster
                hits = np.compress(event[0]['clusterID'] == clusterID, event[0], axis=0) 
                energies = hits['en'] 
                for energy in energies: 
                    clusterEnergies.append(energy)                                                  #|Add the energy to the cluster energies
            ROIEnergy = np.sum(clusterEnergies)
            # Append to original lists
            genEnergies.append(genEnergy[i])
            sumEnergies.append(ROIEnergy)
        pbar.update(count)
        count += 1
    pbar.finish()
    # np.save("sums.npy", sumEnergies)
    # np.save("gens.npy", genEnergies)
    # Plot it
    Plotter.sumEnergyVsGenEnergy(sumEnergies, genEnergies) 

def calculate_switch_length(inheritance, positions, ignore_size=0,
                            index_only=False):
    assert inheritance.shape[0] == positions.size

    # only 1s and 2s are relevant
    exclude = np.any(inheritance < 3, axis=1)
    inh_copy = np.compress(exclude, inheritance.copy(), axis=0)

    forgiven = [forgive(col, ignore_size) for col in inh_copy.T]
    switches = [derive_position_switch_array(np.compress(fgv, col))
                for col, fgv in zip(inh_copy.T, forgiven)]

    filtered_pos = None
    if index_only:
        mean_length = [np.mean(s) for s in switches]
        medi_length = [np.median(s) for s in switches]
        maxi_length = [np.median(s) for s in switches]
    else:
        assert inheritance.shape[0] == positions.shape[0]
        pos = np.compress(exclude, positions)

        filtered_pos = [np.insert(np.take(np.compress(fgv, pos),
                                          sw.cumsum() - 1), 0, pos[0])
                        for fgv, sw in zip(forgiven, switches)]

        mean_length = np.array([np.mean(np.diff(f)) for f in filtered_pos])
        medi_length = np.array([np.median(np.diff(f)) for f in filtered_pos])
        maxi_length = np.array([np.max(np.diff(f)) for f in filtered_pos])

    return mean_length, medi_length, maxi_length, filtered_pos

def contributions(Ilength, Olength, scale, kernel,kernel_width):
    # Antialiasing for downsizing
    if scale < 1:
        h = lambda x: kernel(x,scale)
        kernel_width = kernel_width/scale
    else:
        h = kernel

    # output space coordinate
    x            = np.arange(Olength, dtype = float)
    x.shape     += (1,)
    # input space coord so that 0.5 in Out ~ 0.5 in In, and 0.5+scale in Out ~
    # 0.5 + 1 in In
    u            = x/scale + 0.5*(-1+1.0/scale)
    left         = np.floor(u-kernel_width/2)
    P            = math.ceil(kernel_width) + 2
    indices      = left + np.arange(P)
    weights      = h(u - indices)
    norm         = np.sum(weights,axis=1)
    norm.shape  += (1,)
    weights      = weights/norm
    indices      = np.minimum(np.maximum(0,indices),Ilength-1)
    indices      = np.array(indices,dtype                      = int)


    kill    = np.ma.any(weights,0)
    weights = np.compress(kill,weights,1)
    indices = np.compress(kill,indices,1)
    return (weights,indices)

 def lcylimits(self):
     """Determine the y-limts depending on what plots are selected """
     mask = (self.dtime > self.lcx1)*(self.dtime<self.lcx2)*(self.goodframes>0)
     if self.ratiovar.get():
         rarr=np.compress(mask,self.ratio)
         y1=rarr.min()
         y2=rarr.max()
         ylabel='Star1/Star2'
     else:
         if self.star2var.get() and self.star1var.get():
             cfarr=np.compress(mask,self.cflux).max()
             tfarr=np.compress(mask,self.tflux).max()
             y1=0
             y2=cfarr < tfarr and tfarr or cfarr
             ylabel='Star Flux'
         elif self.star2var.get():
             cfarr=np.compress(mask,self.cflux)
             y1=0
             y2=cfarr.max()
             ylabel='Star2 Flux'
         else:
             tfarr=np.compress(mask,self.tflux)
             y1=0
             y2=tfarr.max()
             ylabel='Star1 Flux'
     return y1, y2, ylabel