def solve(self, pixel=None, solver='cvxopt.coneqp'):
        '''
        Solves minimization problem using quadratic program formulation.
        '''
        n, m = self.image.shape
        if not pixel:
            i, j = floor(n/2), floor(m/2)
        else:
            i, j = pixel
        r, p, q = self.psf_tensor.shape
        psf_bbox = (max(0, p/2-i), max(0, q/2-j),
                    min(p, p/2-i+n), min(q, q/2-j+m))
        img_bbox = (max(0, i-n/2), max(0, j-q/2),
                    min(n, i+n/2), min(m, j+m/2))
        self.A = self._compute_A(psf_bbox)
        self.normalization_factor = self._compute_normalization(psf_bbox)
        self.y = self._compute_y(img_bbox)

        Q = self._compute_Q()
        p = self._compute_p()
        G = self._compute_G()
        h = self._compute_h()
        result = cvxopt.solvers.qp(Q, p, G, h)

        self.result = result
        self.x = result['x']

        return result['status']

def svd_example():
    a = np.floor(np.random.rand(4, 4)*20-6)
    logger.info("Matrix A:\n %s", a)
    b = np.floor(np.random.rand(4, 1)*20-6)
    logger.info("Matrix B:\n %s", b)

    u, s, v_t = np.linalg.svd(a) # SVD decomposition of A
    logger.info("Matrix U:\n %s", u)
    logger.info("Matrix S:\n %s", s)
    logger.info("Matrix V(transpose:\n %s", u)

    logger.info("Computing inverse using linalg.pinv")
    # Computing the inverse using pinv
    inv_pinv = np.linalg.pinv(a)
    logger.info("pinv:\n %s", inv_pinv)

    # Computing inverse using matrix decomposition
    logger.info("Computing inverse using svd matrix decomposition")
    inv_svd = np.dot(np.dot(v_t.T, np.linalg.inv(np.diag(s))), u.T)
    logger.info("svd inverse:\n %s", inv_svd)
    logger.info("comparing the results from pinv and svd_inverse:\n %s",
                np.allclose(inv_pinv, inv_svd))

    logger.info("Sol1: Solving x using pinv matrix... x=A^-1 x b")
    result_pinv_x = np.dot(inv_pinv, b)

    logger.info("Sol2: Solving x using svd_inverse matrix... x=A^-1 x b")
    result_svd_x = np.dot(inv_svd, b)

    if not np.allclose(result_pinv_x, result_svd_x):
        raise ValueError('Should have been True')

    def __init__(self, filepath, qually, mode='444'):
        '''
        '''
        imOrig = cv2.imread(filepath,1)
        self.filepath = filepath
        self.mode = mode
        #Taxa de compressão e Redundancia
        self.CRate = 0; self.Redunc = 0
        self.avgBits = 0
        #Qualidade
        self.qually = qually
        #Dimensões da imagem original
        self.Mo, self.No, self.Do = imOrig.shape
        self.r, self.c = [8, 8]       #DIMENSAO DOS BLOCOS
        #TRATA AS DIMENSOES DA IMAGEM
        (self.M, self.N, self.D), self.img = h.adjImg(imOrig)
        #NUMERO DE BLOCOS NA VERTICAL E HORIZONTAL
        self.nBlkRows = int(np.floor(self.M/self.r))
        self.nBlkCols = int(np.floor(self.N/self.c))
        #Gera Tabela de Qunatizaçao
        self.Z = h.genQntb(self.qually)
        #TRANSFORMA DE RGB PARA YCbCr
        self.Ymg = cv2.cvtColor(self.img, cv2.COLOR_BGR2YCR_CB)
        self.NumBits = 0
        if self.Do == 2:
            self.NCHNL = 1
        elif self.Do == 3:
            self.NCHNL = 3
            
#        self.OUTCOMES = self._run_()
        self._run_()

    def connect_composition_III(self):
        synapse_list = []
        Mt3v_list = self.non_columnar_neurons['Mt3v']
        Mt3h_list = self.non_columnar_neurons['Mt3h']
        
        for neuron in Mt3v_list:
            neuron.assign_pos(0., 0.)
        for neuron in Mt3h_list:
            neuron.assign_pos(0., 0.)
        
        rule3synapsesv = self.other_synapse_dict[self.other_synapse_dict['postname'] == 'Mt3v']
        rule3synapsesh = self.other_synapse_dict[self.other_synapse_dict['postname'] == 'Mt3h']
        
        dtnames = rule3synapsesv.dtype.names
        for cartridge in self.cartridges:
            synapse = Synapse(dict(zip(dtnames, [np.asscalar(p) for p in rule3synapsesv[0]])))
            mtn = int(np.floor(cartridge.neurons['L2'].ypos / ((self.hexarray.Y[-1][-1]+1)/4)))
            synapse.link(cartridge.neurons['L2'], Mt3v_list[mtn])
            synapse_list.append(synapse)
            synapse = Synapse(dict(zip(dtnames, [np.asscalar(p) for p in rule3synapsesh[0]])))
            mtn = int(np.floor(cartridge.neurons['L2'].xpos / ((self.hexarray.X[-1][-1]+1)/4)))
            synapse.link(cartridge.neurons['L2'], Mt3h_list[mtn])
            synapse_list.append(synapse)

        self.composition_rules.append({'synapses': synapse_list})

    def _find_tails(self, mag, rounding=True, half=5, full=10, flag=50):
        '''
        Find how many of each of the tail pieces is necessary.  Flag
        specifies the increment for a flag, barb for a full barb, and half for
        half a barb. Mag should be the magnitude of a vector (ie. >= 0).

        This returns a tuple of:

            (*number of flags*, *number of barbs*, *half_flag*, *empty_flag*)

        *half_flag* is a boolean whether half of a barb is needed,
        since there should only ever be one half on a given
        barb. *empty_flag* flag is an array of flags to easily tell if
        a barb is empty (too low to plot any barbs/flags.
        '''

        #If rounding, round to the nearest multiple of half, the smallest
        #increment
        if rounding:
            mag = half * (mag / half + 0.5).astype(np.int)

        num_flags = np.floor(mag / flag).astype(np.int)
        mag = np.mod(mag, flag)

        num_barb = np.floor(mag / full).astype(np.int)
        mag = np.mod(mag, full)

        half_flag = mag >= half
        empty_flag = ~(half_flag | (num_flags > 0) | (num_barb > 0))

        return num_flags, num_barb, half_flag, empty_flag

def interp(pic,flow):
    ys=np.arange(pic.shape[0]*pic.shape[1])/pic.shape[1]
    ud=(flow[:,:,0].reshape(-1)+ys)%pic.shape[0]
    xs=np.arange(pic.shape[0]*pic.shape[1])%pic.shape[1]
    lr=(flow[:,:,1].reshape(-1)+xs)%pic.shape[1]

    u=np.int32(np.floor(ud))
    d=np.int32(np.ceil(ud))%pic.shape[0]
    udiffs=ud-u
    udiffs=np.dstack((udiffs,udiffs,udiffs))
    l=np.int32(np.floor(lr))
    r=np.int32(np.ceil(lr))%pic.shape[1]
    ldiffs=lr-l
    ldiffs=np.dstack((ldiffs,ldiffs,ldiffs))

    ul=pic[u,l,:]
    ur=pic[u,r,:]
    dl=pic[d,l,:]
    dr=pic[d,r,:]


    udl=ul*(1-udiffs)+dl*udiffs
    udr=ur*(1-udiffs)+dr*udiffs
    ans=np.zeros(pic.shape)
    ans[ys,xs,:]=udl*(1-ldiffs)+udr*ldiffs
    return ans

 def getRowCol(self,lat,lon,returnFloat=False):
     """Return data row and column from given geographic coordinates (lat/lon decimal degrees).
     
     :param lat: 
        Input latitude.
     :param lon: 
        Input longitude.
     :param returnFloat: 
        Boolean indicating whether floating point row/col coordinates should be returned.
     :returns: 
        Tuple of row and column.
     """
     ulx = self._geodict['xmin']
     uly = self._geodict['ymax']
     xdim = self._geodict['xdim']
     ydim = self._geodict['ydim']
     #check to see if we're in a scenario where the grid crosses the meridian
     if self._geodict['xmax'] < ulx and lon < ulx:
         lon += 360
     col = (lon-ulx)/xdim
     row = (uly-lat)/ydim
     if returnFloat:
         return (row,col)
     
     return (np.floor(row).astype(int),np.floor(col).astype(int))

def solver_dt(t0, y0, f, dt, tk, method='expl_RK4', tuning='', v=0, outform='arr'):
	
	step=methods[method]
	#we calculate number of steps
	if np.abs(t0-tk/dt)-np.floor(np.abs(t0-tk/dt)) >= 0.5:
		nos = int(np.ceil(np.abs(t0-tk)/dt))
	else:
		nos = int(np.floor(np.abs(t0-tk)/dt))
		
	if v==2: print "Number of steps I`ll take is "+str(nos) 
	
	#in solver_nos above code will be changed!
	
	
	soln = [[t0, y0]]
	
	for i in range(nos-1):
		tmp1, tmp2 = step(f, soln[i][1], soln[i][0], dt, tuning=tuning)
		soln.append([tmp1, tmp2])
	
	if v>=1: print "solved. Took "+str(nos)+" steps to solve on interval ["+str(t0)+";"+str(tk)+"]"
	if v==2: print "solve method was "+method+". Additional tweaks for this solver were:"+tuning
	if v==2: print "output form was "+outform+"(arr- numpy array such that arr[i] = [t, y1, ..., yn], list - nested list list[i]=[t, array(y)])"
	
	if outform=='arr':
		return np.array([[i[0]]+ ([a for a in i[1]] if type(i[1])==np.ndarray else [i[1]])  for i in soln])
	elif outform=='list':
		return soln

def get_border_to_nucleus_properties(segment):


	pts = np.array(segment. isotropic_border_coords)
	res = min([segment.xres, segment.yres, segment.zres])
	centroid = np.floor(np.mean(pts,0))

	centroid_with_offset = np.floor(np.mean(segment.border_coords,0) + np.array([segment.bounding_box.xmin, segment.bounding_box.ymin, segment.bounding_box.zmin]))

	centroid_by_res = centroid_with_offset * np.array([segment.xres, segment.yres, segment.zres])

#	if centroid_with_offset[2] == 26:
#		pdb.set_trace()

	dists =  ((pts[:,0] - centroid[0])**2 + (pts[:,1] - centroid[1])**2 + (pts[:,2] - centroid[2])**2 ) **0.5
	max_dist = dists.max()
	dists = dists/max_dist	

	bins = np.arange(0, 1.05, 0.05)
	dist_hist = histogram(dists, bins = bins) [0]

	
	segment.add_feature("centroid_res", tuple(centroid_by_res))
	segment.add_feature("border_to_nucleus_dist_hist", dist_hist)
	segment.add_feature("border_to_nucleus_dist_mean", np.mean(dists))
	segment.add_feature("border_to_nucleus_dist_std", np.std(dists))
	segment.add_feature("distance_to_border_scale_factor", max_dist)


	pts = np.array(segment.border_coords)
	segment.add_feature("centroid", tuple(centroid_with_offset.astype("int")))

def calc_slit_box_aps_1id(slit_box_corners, inclip=(1, 10, 1, 10)):
    """
    Calculate the clip box based on given slip corners.

    Parameters
    ----------
    slit_box_corners : np.ndarray
        Four corners of the slit box as a 4x2 matrix
    inclip : tuple, optional
        Extra inclipping to avoid clipping artifacts

    Returns
    -------
    Tuple:
        Cliping indices as a tuple of four
        (clipFromTop, clipToBottom, clipFromLeft, clipToRight)

    """
    return (
        np.floor(slit_box_corners[:, 0].min()).astype(
            int) + inclip[0],  # clip top    row
        np.ceil(slit_box_corners[:, 0].max()).astype(
            int) - inclip[1],  # clip bottom row
        np.floor(slit_box_corners[:, 1].min()).astype(
            int) + inclip[2],  # clip left   col
        np.ceil(slit_box_corners[:, 1].max()).astype(
            int) - inclip[3],  # clip right  col
    )

def build_X(bids, bot_or_human):
    X = bot_or_human
    X = X.drop("payment_account", 1)
    X = X.drop("address", 1)

    bids["day"] = np.floor((bids["time"] - startt) / one_day)
    bids["week"] = np.floor(bids["day"] / 7.0)

    # print bids.bidder_id[0]
    print "starting ips"
    X = ip(X, bids)
    print "starting bid order"
    X = bid_order(X, bids)
    print "starting dt"
    X = dt(X, bids)
    print "startin day"
    X = day(X, bids)
    print "starting n_bids"
    X = n_bids(X, bids)
    print "starting urls"
    X = urls(X, bids)
    # print 'starting bid order'
    # X = bid_order(X, bids)
    print "starting countries"
    X = user_countries_per_auction(X, bids)
    print "starting merch"
    X = merch(X, bids)

    return X

def check_orbits(p1, t1, p2, t2, tmn, tmx, tol):
    n1 = t1 + p1 * np.arange(np.floor((tmn-t1)/p1), np.ceil((tmx-t1)/p1))
    n1 = n1[(tmn <= n1) * (n1 <= tmx)]
    n2 = t2 + p2 * np.arange(np.floor((tmn-t2)/p2), np.ceil((tmx-t2)/p2))
    n2 = n2[(tmn <= n2) * (n2 <= tmx)]
    delta = np.fabs(n1[:, None] - n2[None, :])
    return max(len(n1), len(n2)) == np.sum(delta < tol)

def multilook_attributes(atr_dict,lks_az,lks_rg):
    #####
    atr = dict()
    for key, value in atr_dict.iteritems():  atr[key] = str(value)
  
    ##### calculate new data size
    length = int(atr['FILE_LENGTH'])
    width  = int(atr['WIDTH'])
    length_mli = int(np.floor(length/lks_az))
    width_mli  = int(np.floor(width/lks_rg))
  
    ##### Update attributes
    atr['FILE_LENGTH'] = str(length_mli)
    atr['WIDTH']       = str(width_mli)
    try:
        atr['Y_STEP'] = str(lks_az*float(atr['Y_STEP']))
        atr['X_STEP'] = str(lks_rg*float(atr['X_STEP']))
    except: pass
    try:
        atr['AZIMUTH_PIXEL_SIZE'] = str(lks_az*float(atr['AZIMUTH_PIXEL_SIZE']))
        atr['RANGE_PIXEL_SIZE']   = str(lks_rg*float(atr['RANGE_PIXEL_SIZE']))
    except: pass
  
    try:
        atr['ref_y'] = str(int(int(atr['ref_y'])/lks_az))
        atr['ref_x'] = str(int(int(atr['ref_x'])/lks_rg))
    except: pass
    try:
        atr['subset_y0'] = str(int(int(atr['subset_y0'])/lks_az))
        atr['subset_y1'] = str(int(int(atr['subset_y1'])/lks_az))
        atr['subset_x0'] = str(int(int(atr['subset_x0'])/lks_rg))
        atr['subset_x1'] = str(int(int(atr['subset_x1'])/lks_rg))
    except: pass
  
    return atr

 def __init__(self, hipparcos, symbad, intervals):
     self.hipparcos = hipparcos
     self.symbad = symbad
     self.ra_min = min(hipparcos['RA_J2000'])
     self.ra_max = max(hipparcos['RA_J2000'])
     self.de_min = min(hipparcos['DE_J2000'])
     self.de_max = max(hipparcos['DE_J2000'])
     self.intervals = intervals
     self.ra_dif = self.ra_max - self.ra_min
     self.de_dif = self.de_max - self.de_min
     self.ra_sp = self.ra_dif/(self.intervals-1)
     self.de_sp = self.de_dif/(self.intervals-1)
     
     # set the position at the space grid for hipparcos catalog
     self.stars = pandas.DataFrame(
                numpy.zeros((len(self.hipparcos),2)), columns=['rec','dec'])
     self.stars['rec'] = numpy.floor(
         (self.hipparcos['RA_J2000'] - self.ra_min) / self.ra_sp) + 1
     self.stars['dec'] = numpy.floor(
         (self.hipparcos['DE_J2000'] - self.de_min) / self.de_sp) + 1
     
     # set the position at the space grid for symbad catalog
     self.symbad['rac'] = numpy.floor(
         (self.symbad['RA_J2000'] - self.ra_min) / self.ra_sp) + 1
     self.symbad['dec'] = numpy.floor(
         (self.symbad['DE_J2000'] - self.de_min) / self.de_sp) + 1

 def _get_interv_graticule(self,pmin,pmax,dpar,mmin,mmax,dmer,verbose=True):
     def set_prec(d,n,nn=2):
         arcmin=False
         if d/n < 1.:
             d *= 60
             arcmin = True
             nn = 1
         x = d/n
         y = nn*x
         ex = np.floor(np.log10(y))
         z = np.around(y/10**ex)*10**ex/nn
         if arcmin:
             z = 1./np.around(60./z)
         return z
     max_n_par = 18
     max_n_mer = 36
     n_par = (pmax-pmin)/dpar
     n_mer = (mmax-mmin)/dmer
     if n_par > max_n_par:
         dpar = set_prec((pmax-pmin)/dtor,max_n_par/2)*dtor
     if n_mer > max_n_mer:
         dmer = set_prec((mmax-mmin)/dtor,max_n_mer/2,nn=1)*dtor
     if dmer/dpar < 0.2 or dmer/dpar > 5.:
         dmer = dpar = max(dmer,dpar)
     vdeg = np.floor(np.around(dpar/dtor,10))
     varcmin = (dpar/dtor-vdeg)*60.
     if verbose: print "The interval between parallels is %d deg %.2f'."%(vdeg,varcmin)
     vdeg = np.floor(np.around(dmer/dtor,10))
     varcmin = (dmer/dtor-vdeg)*60.
     if verbose: print "The interval between meridians is %d deg %.2f'."%(vdeg,varcmin)
     return dpar,dmer

def plot_scatter_with_histograms(xvals, yvals, colour='k', oneToOneLine=True, xlabel=None, ylabel=None, title=None):
    gs = gridspec.GridSpec(5, 5)
    xmin = np.floor(min(xvals))
    xmax = np.ceil(max(xvals))
    ymin = np.floor(min(yvals))
    ymax = np.ceil(max(yvals))
    plt.subplot(gs[1:, 0:4])
    plt.plot(xvals, yvals, 'o', color=colour)
    if xlabel is not None:
        plt.xlabel(xlabel)
    if ylabel is not None:
        plt.ylabel(ylabel)
    if oneToOneLine:
        oneToOneMax = max([max(xvals),max(yvals)])
        plt.plot([0,oneToOneMax],[0,oneToOneMax],'b--')
    plt.xlim(xmin,xmax)
    plt.ylim(ymin,ymax)
    plt.subplot(gs[0, 0:4])
    plt.hist(xvals, np.linspace(xmin,xmax,50))
    plt.axis('off')
    plt.subplot(gs[1:,4])
    plt.hist(yvals, np.linspace(ymin,ymax,50), orientation='horizontal')
    plt.axis('off')
    if title is not None:
        plt.suptitle(title)

def qwtCanvasClip(canvas, canvasRect):
    x1 = np.ceil(canvasRect.left())
    x2 = np.floor(canvasRect.right())
    y1 = np.ceil(canvasRect.top())
    y2 = np.floor(canvasRect.bottom())
    r = QRect(x1, y1, x2-x1-1, y2-y1-1)
    return canvas.borderPath(r)

def PlotErrBoxPlot(x, y, delta, ax, showXTicks):
  if x.size < 1:
    return
  ids = np.floor((x)/delta).astype(np.int)
  data = []
  for i in range(ids.min(), ids.max()+1):
    if (ids==i).any():
      data.append(y[ids==i])
  bp = plt.boxplot(data)
#  plt.plot(x,y,'.', color=c1, alpha=0.3)
  # set xticks
  if showXTicks:
    ticks = np.floor((np.arange(ids.min(), ids.max()+1)+0.5)*delta).astype(np.int)
    if np.unique(ticks).size < ticks.size:
      ticks = np.floor((np.arange(ids.min(), ids.max()+1)+0.5)*delta*10.)/10.
    xtickNames = plt.setp(ax, xticklabels=ticks)
    plt.setp(xtickNames, rotation=45)
  else:
    plt.setp(ax.get_xticklabels(), visible=False) 
  for box in bp["boxes"]:
    box.set(color=c1)
    #box.set(facecolor=c1)
  for whisker in bp["whiskers"]:
    whisker.set(color=c1)
  for cap in bp["caps"]:
    cap.set(color=c1)
  for median in bp["medians"]:
    median.set(color=c2)
  for flier in bp["fliers"]:
    flier.set(color=c3, marker=".", alpha=0.15) #,s=6)

def _search_fine(sino, srad, step, init_cen, ratio, drop):
    """
    Fine search for finding the rotation center.
    """
    Nrow, Ncol = sino.shape
    centerfliplr = (Ncol + 1.0) / 2.0 - 1.0

    # Use to shift the sinogram 2 to the raw CoR.
    shiftsino = np.int16(2 * (init_cen - centerfliplr))
    _copy_sino = np.roll(np.fliplr(sino[1:]), shiftsino, axis=1)
    lefttake = 0
    righttake = Ncol - 1
    if init_cen <= centerfliplr:
        lefttake = np.ceil(srad + 1)
        righttake = np.floor(2 * init_cen - srad - 1)
    else:
        lefttake = np.ceil(
            init_cen - (Ncol - 1 - init_cen) + srad + 1)
        righttake = np.floor(Ncol - 1 - srad - 1)
    Ncol1 = righttake - lefttake + 1
    mask = _create_mask(2 * Nrow - 1, Ncol1, 0.5 * ratio * Ncol, drop)
    numshift = np.int16((2 * srad + 1.0) / step)
    listshift = np.linspace(-srad, srad, num=numshift)
    listmetric = np.zeros(len(listshift), dtype='float32')
    num1 = 0
    for i in listshift:
        _sino = ndimage.interpolation.shift(
            _copy_sino, (0, i), prefilter=False)
        sinojoin = np.vstack((sino, _sino))
        listmetric[num1] = np.sum(np.abs(np.fft.fftshift(
            pyfftw.interfaces.numpy_fft.fft2(
                sinojoin[:, lefttake:righttake + 1]))) * mask)
        num1 = num1 + 1
    minpos = np.argmin(listmetric)
    return init_cen + listshift[minpos] / 2.0

def errorScalingFactor(observable, beta):
  """
  Look up the numerical factors to apply to the sky averaged parallax error in order to obtain error
  values for a given astrometric parameter, taking the Ecliptic latitude and the number of transits into
  account.

  Parameters
  ----------

  observable - Name of astrometric observable (one of: alphaStar, delta, parallax, muAlphaStar, muDelta)
  beta       - Values(s) of the Ecliptic latitude.

  Returns
  -------

  Numerical factors to apply to the errors of the given observable.
  """
  if isscalar(beta):
    index=int(floor(abs(sin(beta))*_numStepsSinBeta))
    if index == _numStepsSinBeta:
      return _astrometricErrorFactors[observable][_numStepsSinBeta-1]
    else:
      return _astrometricErrorFactors[observable][index]
  else:
    indices = array(floor(abs(sin(beta))*_numStepsSinBeta), dtype=int)
    indices[(indices==_numStepsSinBeta)] = _numStepsSinBeta-1
    return _astrometricErrorFactors[observable][indices]