def test_quasi_newton_broyden_bad():
    op = p.OptimizationProblem(c.chebyquad)
    guess=linspace(0,1,4)
    cn  = p.QuasiNewtonBroydenBad(op)
    assert near(sol4,(cn.optimize(guess)))
    guess=linspace(0,1,8)
    assert near(sol8,(cn.optimize(guess)))

def extract_seg_plot(data, output_filename_base):
    subclone_plot_file_name = output_filename_base + '.MixClone.segplot.png'
    
    subclone_prev_lst = []
    copynumber_lst = []
    seg_num = data.seg_num

    print "Extracting segments plot file..."
    sys.stdout.flush()

    for j in range(0, seg_num):
        if data.segments[j].baseline_label == True or data.segments[j].allele_type == 'PM':
            continue
        subclone_prev_lst.append(data.segments[j].subclone_prev)
        copynumber_lst.append(data.segments[j].copy_number)

    X = len(subclone_prev_lst)
    
    plt.figure(figsize=(8,8), dpi = 150)
    plt.plot(range(1, X+1), subclone_prev_lst, 'o')
    plt.xlim(0, X+1)
    plt.ylim(0, 1)
    plt.xlabel('Copy number')
    plt.ylabel('Subclonal cellular prevalence')
    plt.xticks(sp.linspace(1, X, X), copynumber_lst)
    plt.yticks(sp.linspace(0, 1, 11), ['0%', '10%', '20%', '30%', '40%', '50%', '60%', '70%', '80%', '90%', '100%'])
    plt.savefig(subclone_plot_file_name, bbox_inches='tight')

    def split(self, sagi, meri):
        """ utilizes geometry.grid to change the rectangle into a generalized surface,
        it is specified with a single set of basis vectors to describe the meridonial,
        normal, and sagittal planes."""
        ins = float((sagi - 1))/sagi
        inm = float((meri - 1))/meri
        stemp = self.norm.s/sagi
        mtemp = self.meri.s/meri

        z,theta = scipy.meshgrid(scipy.linspace(-self.norm.s*ins,
                                                self.norm.s*ins,
                                                sagi),
                                 scipy.linspace(-self.meri.s*inm,
                                                self.meri.s*inm,
                                                meri))

        vecin =geometry.Vecr((self.sagi.s*scipy.ones(theta.shape),
                              theta+scipy.pi/2,
                              scipy.zeros(theta.shape))) #this produces an artificial
        # meri vector, which is in the 'y_hat' direction in the space of the cylinder
        # This is a definite patch over the larger problem, where norm is not normal
        # to the cylinder surface, but is instead the axis of rotation.  This was
        # done to match the Vecr input, which works better with norm in the z direction
               
        pt1 = geometry.Point(geometry.Vecr((scipy.zeros(theta.shape),
                                            theta,
                                            z)),
                             self)

        pt1.redefine(self._origin)

        vecin = vecin.split()

        x_hat = self + pt1 #creates a vector which includes all the centers of the subsurface

        out = []
        #this for loop makes me cringe super hard
        for i in xrange(meri):
            try:
                temp = []
                for j in xrange(sagi):
                    inp = self.rot(vecin[i][j])
                    temp += [Cyl(geometry.Vecx(x_hat.x()[:,i,j]),
                                 self._origin,
                                 [2*stemp,2*mtemp],
                                 self.sagi.s,
                                 vec=[inp, self.norm.copy()],
                                 flag=self.flag)]
                out += [temp]
            except IndexError:
                inp = self.rot(vecin[i])
                out += [Cyl(geometry.Vecx(x_hat.x()[:,i]),
                            self._origin,
                            [2*stemp,2*mtemp],
                            self.norm.s,
                            vec=[inp, self.norm.copy()],
                            flag=self.flag)]
                

        return out

 def loadImageAsGreyScale(self, subscriber=0):
     im = Image.open(self.paramFilename.value)
     if im.mode == "I;16":
         im = im.convert("I")
         data = scipy.misc.fromimage(im).astype("int16")
     else:
         data = scipy.misc.fromimage(im, flatten=True)
     Ny, Nx = data.shape
     xUnit = Quantity(self.paramXScale.value.encode("utf-8"))
     xAxis = DataContainer.FieldContainer(
         scipy.linspace(0.0, xUnit.value, Nx, True), xUnit / xUnit.value, longname="x-coordinate", shortname="x"
     )
     if self.paramYScale.value == "link2X":
         yUnit = xUnit * float(Ny) / Nx
     else:
         yUnit = Quantity(self.paramYScale.value.encode("utf-8"))
     yAxis = DataContainer.FieldContainer(
         scipy.linspace(0.0, yUnit.value, Ny, True), yUnit / yUnit.value, longname="y-coordinate", shortname="y"
     )
     try:
         FieldUnit = Quantity(self.paramFieldUnit.value.encode("utf-8"))
     except AttributeError:
         FieldUnit = self.paramFieldUnit.value
     result = DataContainer.FieldContainer(
         data, FieldUnit, longname="Image", shortname="I", dimensions=[yAxis, xAxis]
     )
     result.seal()
     return result

def uvToELz_grid(ulinspace,vlinspace,R=1.,t=-4.,pot='bar',beta=0.,
            potparams=(0.9,0.01,25.*_degtorad,.8,None)):
    """
    NAME:
       uvToELz_grid
    PURPOSE:
       calculate uvToLz on a grid in (u,v)
    INPUT:
       ulinspace, vlinspace - build the grid using scipy's linspace with
                              these arguments
       R - Galactocentric Radius
       t - time to integrate backwards for 
           (interpretation depends on potential)
       pot - type of non-axisymmetric, time-dependent potential
       beta - power-law index of rotation curve
       potparams - parameters for this potential
    OUTPUT:
       final (E,Lz) on grid [nus,nvs,2]
       E=E/vo^2; Lz= Lz/Ro/vo
    HISTORY:
       2010-03-01 - Written - Bovy (NYU)
    """
    us= sc.linspace(*ulinspace)
    vs= sc.linspace(*vlinspace)
    nus= len(us)
    nvs= len(vs)
    out= sc.zeros((nus,nvs,2))
    for ii in range(nus):
        for jj in range(nvs):
            tmp_out= uvToELz(UV=(us[ii],vs[jj]),R=R,t=t,pot=pot,beta=beta,potparams=potparams)
            out[ii,jj,0]= tmp_out[0]
            out[ii,jj,1]= tmp_out[1]
    return out

def volweight(numsplit=(3,3), factor=1, fact2=None, eq='/home/ian/python/g1120824019.01400'):

    b =  TRIPPy.Tokamak(eqtools.EqdskReader(gfile=eq))
    
    rgrid = b.eq.getRGrid()
    zgrid = b.eq.getZGrid()
    rgrid = scipy.linspace(rgrid[0],rgrid[-1],len(rgrid)*factor)
    zgrid = scipy.linspace(zgrid[0],zgrid[-1],len(zgrid)*factor)
    
    twopi2 = twopi(b)
    
    surfs = twopi2[0].split(numsplit[0],numsplit[1])

    out = scipy.zeros((len(rgrid)-1,len(zgrid)-1))

    for i in surfs:
        for j in i:
            surf2 = j
            if fact2 is None:
                surf2 = j
            else:
                surf2 = j.split(fact2,fact2)

            beam = TRIPPy.beam.multiBeam(surf2,twopi2[1])
            b.trace(beam)
            #TRIPPy.plot.mayaplot.plotLine(beam)
            out += TRIPPy.beam.volWeightBeam(beam,rgrid,zgrid)

    
    return out

def getImageDescriptor(model, im, conf):
	im = standardizeImage(im)
	height, width = im.shape[:2]
	numWords = model.vocab.shape[1]
	frames, descrs = getPhowFeatures(im, conf.phowOpts)
	# quantize appearance
	if model.quantizer == 'vq':
		binsa, _ = vq(descrs.T, model.vocab.T)
	elif model.quantizer == 'kdtree':
		raise ValueError('quantizer kdtree not implemented')
	else:
		raise ValueError('quantizer {0} not known or understood'.format(model.quantizer))
	hist = []
	for n_spatial_bins_x, n_spatial_bins_y in zip(model.numSpatialX, model.numSpatialX):
		binsx, distsx = vq(frames[0, :], linspace(0, width, n_spatial_bins_x))
		binsy, distsy = vq(frames[1, :], linspace(0, height, n_spatial_bins_y))
		# binsx and binsy list to what spatial bin each feature point belongs to
		if (numpy.any(distsx < 0)) | (numpy.any(distsx > (width/n_spatial_bins_x+0.5))):
			print ("something went wrong")
			import pdb; pdb.set_trace()
		if (numpy.any(distsy < 0)) | (numpy.any(distsy > (height/n_spatial_bins_y+0.5))):
			print ("something went wrong")
			import pdb; pdb.set_trace()
		# combined quantization
		number_of_bins = n_spatial_bins_x * n_spatial_bins_y * numWords
		temp = arange(number_of_bins)
		# update using this: http://stackoverflow.com/questions/15230179/how-to-get-the-linear-index-for-a-numpy-array-sub2ind
		temp = temp.reshape([n_spatial_bins_x, n_spatial_bins_y, numWords])
		bin_comb = temp[binsx, binsy, binsa]
		hist_temp, _ = histogram(bin_comb, bins=range(number_of_bins+1), density=True)
		hist.append(hist_temp)
	
	hist = hstack(hist)
	hist = array(hist, 'float32') / sum(hist)
	return hist

def test_interpolateArray():
	grid_x = scipy.linspace(1, 5, 20)
	grid_y = scipy.linspace(-1, 1, 10)
	def fn(x):
		return scipy.sin(x[0] + x[1])
	((xlist, ylist), f) = applyGrid([grid_x, grid_y], fn)
	fig = plt.figure()
	ax = Axes3D(fig)	
	ax.scatter(xlist, ylist, f.ravel())
	
	grid2_x = scipy.linspace(1, 5, 40)
	grid2_y = scipy.linspace(-1, 1, 20)
	f2 = interpolateArray([grid_x, grid_y], [grid2_x, grid2_y], f)
	xy_list = itertools.product(grid2_x, grid2_y)
	(xlist2, ylist2) = zip(*xy_list)	
	fig = plt.figure()	
	ax = Axes3D(fig)
	ax.scatter(xlist2, ylist2, f2.ravel())
	
	grid3_x = grid_x
	grid3_y = grid_y
	f3 = interpolateArray([grid2_x, grid2_y], [grid3_x, grid3_y], f2)
	xy_list = itertools.product(grid3_x, grid3_y)
	(xlist3, ylist3) = zip(*xy_list)	
	fig = plt.figure()	
	ax = Axes3D(fig)	
	ax.scatter(xlist3, ylist3, f3.ravel())	
	
	
	
	
	
	
	

def subdivide(g, rmin, rmax, m):
    """
    Subdivide the given interval by equal integrals of g.

    **example**
     >>> g = lambda r : r**2
     >>> rmin = 0.0
     >>> rmax = 1.0
     >>> m = 2
     >>> r = subdivide(g, rmin, rmax, m)
     >>> abs(r - array([0.0, .5**(1./3.), 1.0])).sum() <= 1e-4
     True
    """
    r = sp.linspace(rmin, rmax, m+1)

    # total integral value
    tot = quad(g, rmin, rmax)[0]
    cuts = sp.linspace(0.0, tot, m+1)[1:-1]
    
    # define intervals r_n to r_n+1 to give equal area under g(r)
    r[1:-1] = [brentq((lambda r_: quad(g, rmin, r_)[0]-cut), rmin, rmax,
                      xtol=1.0e-4)
                for cut in cuts]
    # return the sequence of subinterval boundaries
    return r

	def draw_plot(self):
		self.axes.clear()
		self.histax.clear()
		global mode
		if mode == 'hsv':
			self.axes.set_xlabel('Hue')
			self.axes.set_ylabel('Saturation')
			self.axes.set_zlabel('Value')
			tubepixels = np.array(list(self.hsv[greenloc[0], greenloc[1]] for greenloc in coloredpixels(self.mask, (0,255,0))))
			nontubepixels = np.array(list(self.hsv[redloc[0], redloc[1]] for redloc in coloredpixels(self.mask, (0,0,255))))
		else:
			self.axes.set_xlabel('Blue')
			self.axes.set_ylabel('Green')
			self.axes.set_zlabel('Red')
			tubepixels = np.array(list(self.smallimg[greenloc[0], greenloc[1]] for greenloc in coloredpixels(self.mask, (0,255,0))))
			nontubepixels = np.array(list(self.smallimg[redloc[0], redloc[1]] for redloc in coloredpixels(self.mask, (0,0,255))))
		onehundredtube = [int(x) for x in scipy.linspace(0,len(tubepixels)-1, 500)]
		onehundrednon = [int(x) for x in scipy.linspace(0,len(nontubepixels)-1, 500)]
		if len(tubepixels):
			self.axes.scatter(tubepixels[onehundredtube,0], tubepixels[onehundredtube,1], tubepixels[onehundredtube,2], c='g')
			self.histax.hist(tubepixels[onehundredtube, 1], 100, normed=1, facecolor='green')
		if len(nontubepixels):
			self.axes.scatter(nontubepixels[onehundrednon,0], nontubepixels[onehundrednon,1], nontubepixels[onehundrednon,2], c='r')
			self.histax.hist(nontubepixels[onehundrednon, 1], 100, normed=1, facecolor='red')
		self.figure.canvas.draw()
		self.histfig.canvas.draw()

    def testSnrFuncs(self):
        """test for signal to noise ratio functions"""

        # trivial
        data_triv = sp.ones((3, 10))
        snr_triv_test = sp.ones(3)
        assert_equal(
            snr_peak(data_triv, 1.0),
            snr_triv_test)
        assert_equal(
            snr_power(data_triv, 1.0),
            snr_triv_test)
        assert_equal(
            snr_maha(data_triv, sp.eye(data_triv.shape[1])),
            snr_triv_test)

        # application
        data = sp.array([
            sp.sin(sp.linspace(0.0, 2 * sp.pi, 100)),
            sp.sin(sp.linspace(0.0, 2 * sp.pi, 100)) * 2,
            sp.sin(sp.linspace(0.0, 2 * sp.pi, 100)) * 5,
        ])
        assert_equal(
            snr_peak(data, 1.0),
            sp.absolute(data).max(axis=1))
        assert_equal(
            snr_power(data, 1.0),
            sp.sqrt((data * data).sum(axis=1) / data.shape[1]))
        assert_almost_equal(
            snr_maha(data, sp.eye(data.shape[1])),
            sp.sqrt((data * data).sum(axis=1) / data.shape[1]))

def test_newton_exact_line_search():
    op = p.OptimizationProblem(c.chebyquad)
    guess=linspace(0,1,4)
    cn  = p.NewtonExactLine(op)
    assert near(sol4,(cn.optimize(guess)))
    guess=linspace(0,1,8)
    assert near(sol8,(cn.optimize(guess)))

def test_classic_newton():
    op = p.OptimizationProblem(c.chebyquad)
    guess=linspace(0,1,4)
    cn  = p.ClassicNewton(op)
    assert near(sol4,(cn.optimize(guess)))
    guess=linspace(0,1,8)
    assert near(sol8,(cn.optimize(guess)))

 def get_znodes(self):
     " Compute a nodes for a log-lower and linear-upper grid. "
     zlower = sp.exp(sp.linspace(sp.log(self.zmin), sp.log(self.zmid), self.Nlo))
     zupper = sp.linspace(self.zmid, self.zmax, self.Nhi)
     znodes = sp.concatenate([zlower, 
                              zupper[1:]])
     return znodes

def add_intron_patch2(ax, start, stop, cnt, color='green'):
    ### compute a quadratic function through the three points

    ### we set the first root to 0 and shift only the plotting ...
    x2 = ((stop - start) / 2.0)
    x3 = float(stop - start)

    ### compute coefficients
    #z = float((x1*x1*x2 + x1*x3*x3 + x2*x2*x3) - (x3*x3*x2 + x2*x2*x1 + x1*x1*x3))
    z = float((x2*x2*x3) - (x3*x3*x2))
    if z == 0:
        return

    #a = float(cnt) * (x3 - x1) / z
    #b = float(cnt) * ((x1*x1) - (x3*x3)) / z
    #c = float(cnt) * ((x1*x3*x3) - (x1*x1*x3)) / z
    a = float(cnt) * x3 / z
    b = float(cnt) * (-1*(x3*x3)) / z

    ### get points
    #x = sp.linspace(start, stop, 100)
    x = sp.linspace(0, stop-start, 100)
    #y = (a*x*x) + (b*x) + c
    y = (a*x*x) + (b*x)
    ax.plot(sp.linspace(start, stop, 100), y, '-', color=color)

def plot_disc_norm():
    x = sp.linspace(-3,3,100);
    y = st.norm.pdf(x,0,1);
    fig, ax = plt.subplots()
    fig.canvas.draw()
    
    ax.plot(x,y)
    
    fill1_x = sp.linspace(-2,-1.5,100)
    fill1_y = st.norm.pdf(fill1_x,0,1)
    fill2_x = sp.linspace(-1.5,-1,100)
    fill2_y = st.norm.pdf(fill2_x,0,1)
    ax.fill_between(fill1_x,0,fill1_y,facecolor = 'blue', edgecolor = 'k',alpha = 0.75)
    ax.fill_between(fill2_x,0,fill2_y,facecolor = 'blue', edgecolor = 'k',alpha = 0.75)
    for label in ax.get_yticklabels():
        label.set_visible(False)
    for tick in ax.get_xticklines():
        tick.set_visible(False)
    for tick in ax.get_yticklines():
        tick.set_visible(False)
    
    plt.rc("font", size = 16)
    plt.xticks([-2,-1.5,-1])
    labels = [item.get_text() for item in ax.get_xticklabels()]
    labels[0] = r"$v_k$"
    labels[1] = r"$\varepsilon_k$"
    labels[2] = r"$v_{k+1}$"
    ax.set_xticklabels(labels)
    plt.ylim([0, .45])

    
    plt.savefig('discnorm.pdf')

 def loadImageAsGreyScale(self, subscriber=0):
     im = PIL.Image.open(self.paramFilename.value)
     if im.mode == "I;16":
         im = im.convert("I")
         data = scipy.misc.fromimage(im).astype("int16")
     else:
         data = scipy.misc.fromimage(im, flatten=True)
     Ny, Nx = data.shape
     xUnit = parseFCUnit(self.paramXScale.value)
     xAxis = DataContainer.FieldContainer(
         scipy.linspace(0.0, self._getScalar(xUnit), Nx, True),
         xUnit / self._getScalar(xUnit),
         longname="x-coordinate",
         shortname="x",
     )
     if self.paramYScale.value == "link2X":
         yUnit = xUnit * float(Ny) / Nx
     else:
         yUnit = parseFCUnit(self.paramYScale.value)
     yAxis = DataContainer.FieldContainer(
         scipy.linspace(0.0, self._getScalar(yUnit), Ny, True),
         yUnit / self._getScalar(yUnit),
         longname="y-coordinate",
         shortname="y",
     )
     fieldUnit = parseFCUnit(self.paramFieldUnit.value)
     result = DataContainer.FieldContainer(
         data, fieldUnit, longname="Image", shortname="I", dimensions=[yAxis, xAxis]
     )
     result.seal()
     return result

def joint(dataA, dataB, g, minLimitA, maxLimitA, minLimitB, maxLimitB, alpha=2):
    '''Calculates a joint PDF based on the DEFT algorithm.

    Keyword arguments:
    dataA           data variable A NxP
                        N = elements
                        P = population
    dataB           data variable B NxP
                        N = elements
                        P = population
    minLimitA       Minimal possible value for variable A
    maxLimitA       Maximal possible value for variable A
    minLimitB       Minimal possible value for variable B
    maxLimitB       Maximal possible value for variable B
    g               Number of bins on which to compute the estimated density
    alpha           The smoothness parameter, specifying which derivative of the filed is constrained by the prior. May be any integer >= 1.
    Returns:    
                    pdf NxB
                        N = elements
                        BA = Bin A
                        BB = Bin B
    '''
    assert alpha >= 1, "Alpha must be greater or equal than 1"
    assert g >= 1, "Number of bins must be greater or equal than 1"
    result = np.zeros((dataA.shape[0], g, g))
    bbox = [minLimitA, maxLimitA, minLimitB, maxLimitB]
    xgridA = sp.linspace(bbox[0], bbox[1], g)
    xgridB = sp.linspace(bbox[2], bbox[3], g)

    for i in xrange(dataA.shape[0]):
        if len(dataA[i,:]) > 1 and len(dataB[i,:]) > 1:
            da = deft_2d(dataA[i,:], dataB[i,:], bbox, G=g, alpha=alpha)(xgridA, xgridB)
            result[i, :, :] = da

    return result

def create_grid(r_in, r_out, nshell, space = 'powerlaw1', end = True):
    # function to create grid
    if space == 'log10':
        from scipy import log10, logspace
        # get the exponent of the start- and
        # stop-radius in input units
        start = [log10(r_in), 0][r_in == 0]
        stop = log10(r_out)
        radii = logspace(start, stop, num=nshell, endpoint=end)
    elif space == "powerlaw1":
        from scipy import arange
        radii = r_in * (r_out/r_in)**(arange(nshell)/(nshell - 1.0))
    elif space == 'linear':
        from scipy import linspace
        # linearly spaced grid
        radii = linspace(r_in, r_out, num=nshell, endpoint=end)
    elif space == 'powerlaw2':
        from scipy import linspace
        # first check if coefficients to the power-law was given
        #~ if 'exp' in kwargs:
            #~ p_exp = kwargs['exp']
        #~ else: # if not, set it to 2, i.e. r^2
            #~ p_exp = 2
        radii = r_in + (r_out - r_in)*(linspace(r_in, r_out, num=nshell, endpoint=end)/(r_out))**2
        #pr_int('Not implemented yet.')
        #raise ParError(spaced)
    else:
        raise Exception(space)
    return radii

    def generateNodesAdaptive(self):
        innerDomainSize = self.innerDomainSize
        innerMeshSize   = self.innerMeshSize
        numberElementsInnerDomain = innerDomainSize/innerMeshSize
	assert(numberElementsInnerDomain < self.numberElements)
        domainCenter = (self.domainStart+self.domainEnd)/2
        nodes0 = np.linspace(domainCenter,innerDomainSize/2.0,(numberElementsInnerDomain/2.0)+1.0)
        nodes0 = np.delete(nodes0,-1)
        numberOuterIntervalsFromDomainCenter = (self.numberElements - numberElementsInnerDomain)/2.0
        const = np.log2(innerDomainSize/2.0)/0.5
        exp = np.linspace(const,np.log2(self.domainEnd*self.domainEnd),numberOuterIntervalsFromDomainCenter+1)
        nodes1 = np.power(np.sqrt(2),exp)
        nodesp = np.concatenate((nodes0,nodes1))
        nodesn = -nodesp[::-1]
        nodesn = np.delete(nodesn,-1)
        linNodalCoordinates = np.concatenate((nodesn,nodesp))
        nodalCoordinates = 0

        #Introduce higher order nodes
        if self.elementType == "quadratic" or self.elementType == "cubic":
           if self.elementType == "quadratic":
              numberNodesPerElement = 3 
           elif self.elementType == "cubic":
              numberNodesPerElement = 4

           for i in range(0,len(linNodalCoordinates)-1):
              newnodes = np.linspace(linNodalCoordinates[i],linNodalCoordinates[i+1],numberNodesPerElement)
              nodalCoordinates = np.delete(nodalCoordinates,-1)
              nodalCoordinates = np.concatenate((nodalCoordinates,newnodes))

        else:
           nodalCoordinates = linNodalCoordinates
    
        return nodalCoordinates