def scaleSignal(img, fitParams=None, backgroundToZero=False, reference=None):
    '''
    scale the image between...
    
    backgroundToZero=True -> 0 (average background) and 1 (maximum signal)
    backgroundToZero=False -> signal+-3std
    
    reference -> reference image -- scale image to fit this one
    
    returns:
    scaled image
    '''
    img = imread(img)
    if reference is not None: 
        low, high = signalRange(img, fitParams)
        low2, high2 = signalRange(reference)
        img = np.asfarray(img)
        ampl = (high2-low2)/(high-low)
        img-=low
        img *= ampl
        img += low2
        return img
    else:
        offs, div = scaleParams(img, fitParams, backgroundToZero)
        img = np.asfarray(img)  - offs 
        img /= div 
        print 'offset: %s, divident: %s' %(offs, div)
        return img

    def SetInitialPoints(self, x0, radius=0.05):
        """Set Initial Points with Guess (x0)

input::
    - x0: must be a sequence of length self.nDim
    - radius: generate random points within [-radius*x0, radius*x0]
        for i!=0 when a simplex-type initial guess in required"""
        x0 = asfarray(x0)
        rank = len(x0.shape)
        if rank is 0:
            x0 = asfarray([x0])
            rank = 1
        if not -1 < rank < 2:
            raise ValueError, "Initial guess must be a scalar or rank-1 sequence."
        if len(x0) != self.nDim:
            raise ValueError, "Initial guess must be length %s" % self.nDim

        #slightly alter initial values for solvers that depend on randomness
        min = x0*(1-radius)
        max = x0*(1+radius)
        numzeros = len(x0[x0==0])
        min[min==0] = asarray([-radius for i in range(numzeros)])
        max[max==0] = asarray([radius for i in range(numzeros)])
        self.SetRandomInitialPoints(min,max)
        #stick initial values in population[i], i=0
        self.population[0] = x0.tolist()

    def testCopyConstructorBP(self, level=1):
	""" test if a copied bandpass ESN generates the same result """
        
	# set bandpass parameters
	self.net.setSimAlgorithm(SIM_BP)
	f1 = N.linspace(0.1, 1., self.net.getSize())
	f2 = N.linspace(0.0001, 0.5, self.net.getSize())
	self.net.init()
	self.net.setBPCutoff(f1,f2)
	
	# set output weight matrix
	trainin = N.random.rand(self.ins,self.train_size) * 2 - 1
	trainout = N.random.rand(self.outs,self.train_size) * 2 - 1
	trainin = N.asfarray(trainin, self.dtype)
	trainout = N.asfarray(trainout, self.dtype)
	self.net.train(trainin,trainout,1)
	
	# copy network
	# ATTENTION: operator= is shallow copy !
	if self.dtype is 'float32':
		netA = SingleESN(self.net)
	else:
		netA = DoubleESN(self.net)
	
	# simulate both networks separate and test result
	indata = N.random.rand(self.ins,self.sim_size)*2-1
	indata = N.asfarray(indata, self.dtype)
	outdata = N.empty((self.outs,self.sim_size),self.dtype)
	outdataA = N.empty((self.outs,self.sim_size),self.dtype)
	self.net.simulate( indata, outdata )
	netA.simulate( indata, outdataA )
	assert_array_almost_equal(outdata,outdataA)

    def testPI(self, level=1):
	""" test TRAIN_PI with zero input and feedback """
        
	# init network
	self.net.setSimAlgorithm(SIM_STD)
	self.net.setTrainAlgorithm(TRAIN_PI)
	self.net.init()
	
	# train network
	washout = 2
	# test with zero input:
	indata = N.zeros((self.ins,self.train_size),self.dtype)
	outdata = N.random.rand(self.outs,self.train_size) * 2 - 1
	indata = N.asfarray( indata, self.dtype )
	outdata = N.asfarray( outdata, self.dtype )
	self.net.train( indata, outdata, washout )
	wout_target = self.net.getWout().copy()
	
	# teacher forcing, collect states
	X = self._teacherForcing(indata,outdata)
	
	# restructure data
	M = N.r_[X,indata]
	M = M[:,washout:self.train_size].T
	T = outdata[:,washout:self.train_size].T
	
	# calc pseudo inverse: wout = pinv(M) * T
	wout = ( N.dot(pinv(M),T) ).T
	
	# normalize result for comparison
	wout = wout / abs(wout).max()
	wout_target = wout_target / abs(wout_target).max()
	assert_array_almost_equal(wout_target,wout,2)

def numerical_partials(f, p, f0=None, pmin=None, pmax=None, prel=1.e-6,
                       pabs=1.e-9, pmask=None, args=(), kwargs=None):
    """Compute partial derivatives of f(p) wrt p by finite differences."""
    if kwargs is None:
        kwargs = {}
    f0, p = f(p) if f0 is None else asfarray(f0), asfarray(p)
    dp = zeros_like(p)
    prel, pabs = prel + dp, pabs + dp
    dp = maximum(prel*absolute(p), pabs)
    pfull = p.copy()
    if pmask is not None:
        p, dp = p[pmask], dp[pmask]
    else:
        pmask = ones(p.shape, dtype=bool)
    # Assume that pmin <= p <= pmax, but check p+dp.
    if pmax is not None:
        mask = p+dp > pmax
        dp[mask] *= -1
        if mask.any():
            if pmin is not None and (p+dp < pmin).any():
                raise ValueError("pmin and pmax too close together")
    dfdp = []
    for dp, p1 in zip(dp, p+diag(dp)):
        if not dp:
            raise ValueError("zero step size, check prel and pabs")
        pfull[pmask] = p1
        dfdp.append((f(pfull, *args, **kwargs) - f0)/dp)
    return array(dfdp)

    def __init__(self,ad):
        """Load parameters from a FITS header into the fitting function such that
           we can evaluate the function.
           The header has keywords of the form COEFF_A, COEFF_B..., from low
           to high order.
           :param xmin,xmax:  Min,max range where to eveluate
           :param fitname:  Fitting function name used.
           :type fitname: {'polynomial','legendre',or 'chebyshev'}
           :param order: polynomial order used to evaluate the coeffs
           :param coeff: Coefficients from the fitting
           :type coeff: List
           :return: An array of evaluated values

           Example:
             ef = gfit.Pix2coord(ad)
             ef(1300)   # Evaluate wavelength at pixel 1300 in the middle row
             ef(1300,400)   # Evaluate wavelength at pixel 1300 at row 400
        """

        tbh = ad['WAVECAL',1].header
        pixsample = np.asfarray(tbh['pixrange'].split())
        self.fitname = tbh['fitname']
        self.order = tbh['fitorder']
        self.coeff = np.asfarray(tbh['fitcoeff'].split())
        # Set up the pix2wavelength evaluator function 
        xmin, xmax = pixsample
        self.xnorm = lambda x: 2.*(x - xmin) / (xmax-xmin) - 1

        f3dcoeff = []
        for k in ['A', 'B', 'C', 'D', 'E', 'F']:
            f3dcoeff.append(tbh['COEFF_'+k])
        self.f3dcoeff = np.asfarray(f3dcoeff)

def relu(x,deriv = False):
    if not deriv:
        return np.asfarray(np.maximum(0,x))
    if deriv:
        der = np.asfarray(np.maximum(0,x))
        der[der > 0] = 1
        return der

def snv_on_chromosome(axis, variants, segments, genes,
                      do_trend, do_boost=False):
    # XXX only set x-limits if not already done for probes/segments
    # setup_chromosome(axis, None, segments, variants,
    #                  0.0, 1.0, "VAF")
    axis.set_ylim(0.0, 1.0)
    axis.set_ylabel("VAF")
    axis.set_xlabel("Position (Mb)")
    axis.get_yaxis().tick_left()
    axis.get_xaxis().tick_top()
    axis.tick_params(which='both', direction='out',
                     labelbottom=False, labeltop=False)

    x_mb = variants["start"] * MB
    if do_boost:
        y = variants.tumor_boost()
    else:
        y = np.asfarray(variants["alt_freq"])
    axis.scatter(x_mb, y, color=POINT_COLOR, alpha=0.3)
    # TODO - highlight genes/selection
    if segments:
        # Draw average VAF within each segment
        posns = np.asfarray(variants["start"]) # * MB
        y = np.asfarray(y)
        for v_start, v_end, v_freq in group_snvs_by_segments(posns, y,
                                                             segments):
            # ENH: color by segment gain/loss
            axis.plot([v_start * MB, v_end * MB], [v_freq, v_freq],
                      color='#C0C0C0', linewidth=2, #zorder=1,
                      solid_capstyle='round')

def _det(xvert, yvert):
    """
    Compute twice the area of the triangle defined by points using the
    determinant formula.

    Parameters
    ----------
    xvert : array
        A vector of nodal x-coords.
    yvert : array
        A vector of nodal y-coords.

    Returns
    -------
    area : float
        Twice the area of the triangle defined by the points:
            area is positive if points define polygon in anticlockwise order.
            area is negative if points define polygon in clockwise order.
            area is zero if at least two of the points are concident or if
            all points are collinear.

    """
    xvert = np.asfarray(xvert)
    yvert = np.asfarray(yvert)
    x_prev = np.concatenate(([xvert[-1]], xvert[:-1]))
    y_prev = np.concatenate(([yvert[-1]], yvert[:-1]))
    return np.sum(yvert * x_prev - xvert * y_prev, axis=0)

  def calcNLL(self, doc, state):
    """Calculates the negative log likelihood of the document, given the relevant information. This is the DocState object again, but this time with the entire state object as well. Probability (Expressed as negative log likelihood.) is specificly calculated using all terms that contain a variable in the document, but none that would be identical for all documents. That is, it contains the probability of the cluster, the probability of the dp given the cluster, and the probability of the samples, which factors in both the drawing of the topic and the drawing of the word. The ordering of the samples is considered irrelevant, with both the topic and word defining uniqueness. Some subtle approximation is made - see if you can spot it in the code!"""
    self.nll = 0.0

    # Probability of drawing the cluster...
    self.nll -= math.log(state.clusterUse[doc.cluster])
    self.nll += math.log(state.clusterUse.sum()+state.clusterConc)


    # Probability of drawing the documents dp from its cluster, taking into account the abnormal entrys...
    cl = state.cluster[doc.cluster]
    logBMN = numpy.log(cl[2] / (cl[2]*numpy.asfarray(doc.behFlags)).sum())

    behInstCounts = numpy.zeros(doc.behFlags.shape[0], dtype=numpy.int32)
    instCounts = numpy.zeros(cl[0].shape[0], dtype=numpy.int32)
    for ii in xrange(doc.use.shape[0]):
      behInstCounts[doc.use[ii,0]] += 1
      if doc.use[ii,0]==0: instCounts[doc.use[ii,1]] += 1

    self.nll -= (logBMN * behInstCounts).sum()
    self.nll -= scipy.special.gammaln(behInstCounts.sum() + 1.0)
    self.nll += scipy.special.gammaln(behInstCounts + 1.0).sum()

    norm = cl[0][:,1].sum() + cl[1]
    self.nll -= (numpy.log(numpy.asfarray(cl[0][:,1])/norm)*instCounts).sum()
    self.nll -= scipy.special.gammaln(instCounts.sum() + 1.0) # Cancels with a term from the above - can optimise, but would rather have neat code.
    self.nll += scipy.special.gammaln(instCounts + 1.0).sum()


    # Count the numbers of word/topic instance pairs in the data structure - sum using a dictionary...
    samp_count = collections.defaultdict(int) # [instance,word]
    for s in xrange(doc.samples.shape[0]):
      samp_count[doc.samples[s,0],doc.samples[s,1]] += 1

    # Calculate the probability distribution of drawing each topic instance and the probability of drawing each word/topic assignment...
    inst = numpy.asfarray(doc.use[:,2])
    inst /= inst.sum() + doc.conc
    
    topicWord = numpy.asfarray(state.topicWord) + state.beta
    topicWord = (topicWord.T/topicWord.sum(axis=1)).T

    abnormTopicWord = numpy.asfarray(state.abnormTopicWord) + state.beta
    abnormTopicWord = (abnormTopicWord.T/abnormTopicWord.sum(axis=1)).T

    instLog = numpy.log(inst)
    wordLog = numpy.log(topicWord)
    abnormLog = numpy.log(abnormTopicWord)


    # Now sum into nll the probability of drawing the samples that have been drawn - gets a tad complex as includes the probability of drawing the topic from the documents dp and then the probability of drawing the word from the topic, except I've merged them such that it doesn't look like that is what is happening...
    self.nll -= scipy.special.gammaln(doc.samples.shape[0]+1.0)
    for pair, count in samp_count.iteritems():
      inst, word = pair
      beh = doc.use[inst,0]
      if beh==0:
        topic = cl[0][doc.use[inst,1],0]
        self.nll -= count * (wordLog[topic,word] + instLog[inst])
      else:
        self.nll -= count * (abnormLog[beh,word] + instLog[inst])
      self.nll += scipy.special.gammaln(count+1.0)

def find_peaks(lpix, indl, indr):
    """
      Given the left and right edges of a line list,
      calculates the peak center by simply centroid algorithm
    """

    centres = []
    max_flux = []
    wth = []
    for i0,i1 in zip(indl,indr):
        fl = lpix[i0:i1]
        wa = np.arange(i0,i1)
        if not len(wa): continue
        try:
            ew = len(wa)*fl
            ewtot = np.sum(ew)
            wa_ewtot = np.sum(wa * ew)
            center = wa_ewtot / ewtot
        except:
            center = (i1-i0)/2.

        centres.append(center)
        try:
           if i0==i1:
              print 'FNDPK00:',i0
           max_flux.append(max(fl))
        except:
           print 'FNDPK:',i0,i1
        wth.append(abs(i1-i0))

    return np.asfarray(centres),np.asfarray(max_flux),np.asfarray(wth)

    def __init__(self, X=None, Y=None, Z=None, clipboard=False,
                 x_label=None, x_unit=None,
                 y_label=None, y_unit=None,
                 label=None, unit=None):
                    
        self.sp = None

        if isinstance(X, mesh1d):
            self.X = X
        else:
            self.X = mesh1d(np.asfarray(X, dtype='float64'))

        if isinstance(Y, mesh1d):
            self.Y = Y
        else:
            self.Y = mesh1d(np.asfarray(Y, dtype='float64'))

        self.label = label
        self.unit = unit
        self.Z = mesh1d(np.asfarray(Z, dtype='float64'), self.label, self.unit)  # np.array(Z)

        if clipboard is True:
            self.read_clipboard()
            self.__init_sp()
        elif self.Z.shape[0] > 1:
            self.__init_sp()

 def __solver__(self, p):
     
     p.xk = p.x0.copy()
     p.fk = asfarray((p.f(p.x0)) ** 2).sum().flatten()
         
     p.iterfcn()
     if p.istop:
         p.xf, p.ff = p.xk, p.fk
         return 
     
     if p.userProvided.df:
         xf, cov_x, infodict, mesg, ier = leastsq(p.f, p.x0, Dfun=p.df, xtol = p.xtol, ftol = p.ftol, maxfev = p.maxFunEvals, full_output = 1)
     else:
         xf, cov_x, infodict, mesg, ier = leastsq(p.f, p.x0, xtol = p.xtol, maxfev = p.maxFunEvals, epsfcn = p.diffInt, ftol = p.ftol, full_output = 1)
     
     if ier == 1: p.istop = 1000
     else: p.istop = -1000
     p.msg = mesg
         
     ff = asfarray((p.f(xf)) ** 2).sum().flatten()
     p.xk = xf
     p.fk = ff
     
     p.xf = xf
     p.ff = ff        
     p.iterfcn()

    def _cdf(self, xloc, left, right, cache):
        """
        Cumulative distribution function.

        Example:
            >>> print(chaospy.Uniform().fwd([-0.5, 0.5, 1.5, 2.5]))
            [0.  0.5 1.  1. ]
            >>> print(chaospy.Add(chaospy.Uniform(), 1).fwd([-0.5, 0.5, 1.5, 2.5]))
            [0.  0.  0.5 1. ]
            >>> print(chaospy.Add(1, chaospy.Uniform()).fwd([-0.5, 0.5, 1.5, 2.5]))
            [0.  0.  0.5 1. ]
            >>> print(chaospy.Add(1, 1).fwd([-0.5, 0.5, 1.5, 2.5]))
            [0. 0. 0. 1.]
        """
        left = evaluation.get_forward_cache(left, cache)
        right = evaluation.get_forward_cache(right, cache)

        if isinstance(left, Dist):
            if isinstance(right, Dist):
                raise evaluation.DependencyError(
                    "under-defined distribution {} or {}".format(left, right))
        elif not isinstance(right, Dist):
            return numpy.asfarray(left+right <= xloc)
        else:
            left, right = right, left
        xloc = (xloc.T-numpy.asfarray(right).T).T
        output = evaluation.evaluate_forward(left, xloc, cache=cache)
        assert output.shape == xloc.shape
        return output

def get_fasta_stats(cnarr, fa_fname):
    """Calculate GC and RepeatMasker content of each bin in the FASTA genome."""
    logging.info("Calculating GC and RepeatMasker content in %s ...", fa_fname)
    gc_rm_vals = [calculate_gc_lo(subseq)
                  for subseq in fasta_extract_regions(fa_fname, cnarr)]
    gc_vals, rm_vals = zip(*gc_rm_vals)
    return np.asfarray(gc_vals), np.asfarray(rm_vals)

    def generateTrajectory(self,strokes):
	# Loop over each stroke pair
	current = self.arm.configuration
	configurations = []
	for stroke in strokes:
	    wp1 = stroke[0] # Start point of stroke
	    wp2 = stroke[1] # End point of stroke
	    # Check if current location is current waypoint (i.e. don't lift off)
	    #currentHover = np.copy(self.arm.configuration)	# Set hover over curret point
	    currentHover = np.asfarray([current[0], current[1], current[2] + 2-0])
	    print 'Current:',self.arm.configuration, currentHover
	    target1 = self.paper.paperToWorld(wp1) 	# Compute world coordinates of target
	    target2 = self.paper.paperToWorld(wp2) 	# Compute world coordinates of endpoint
	    print 'WORLD PTS:',target1,target2
	    config1 = self.arm.fullIK(target1) 		# Perform full IK on 1st target
	    config1Hover = config1 + np.asfarray([0,0,-20])	# Hover configuration over 1st target
	    config1Hover[2] = config1Hover[2] - 20	# Set height to hover
	    config2 = self.arm.fullIK(target2) 		# Perform full IK on 2nd target
	    
	    print 'Configs:', str(config1), str(config1Hover), str(currentHover)
	    current = current.reshape((3,))
	    currentHover = currentHover.reshape((3,))
	    config1Hover = config1Hover.reshape((3,))
	    config1 = config1.reshape((3,))
	    config2 = config2.reshape((3,))
	    print 'CONFIG2\n' , str(config2)
	    curHoverPts = interpolateLinear(current, currentHover, 10)
	    current = config2
	    configurations = configurations + [currentHover, config1Hover, config1, config2]
	    
	    configurations = configurations + interpolateLinear(currentHover, config1Hover, 50)
	    configurations = configurations + interpolateLinear(config1Hover, config1, 10)
	    configurations = configurations + interpolateLinear(config1, config2, 50)
	    
	return configurations

def _det(xvert, yvert):
    '''Compute twice the area of the triangle defined by points with using
    determinant formula.

    Input parameters:

    xvert -- A vector of nodal x-coords (array-like).
    yvert -- A vector of nodal y-coords (array-like).

    Output parameters:

    Twice the area of the triangle defined by the points.

    Notes:

    _det is positive if points define polygon in anticlockwise order.
    _det is negative if points define polygon in clockwise order.
    _det is zero if at least two of the points are concident or if
        all points are collinear.

    '''
    xvert = np.asfarray(xvert)
    yvert = np.asfarray(yvert)
    x_prev = np.concatenate(([xvert[-1]], xvert[:-1]))
    y_prev = np.concatenate(([yvert[-1]], yvert[:-1]))
    return np.sum(yvert * x_prev - xvert * y_prev, axis=0)

def CVXOPT_QP_Solver(p, solverName):
    if solverName == 'native_CVXOPT_QP_Solver': solverName = None
    cvxopt_solvers.options['maxiters'] = p.maxIter
    cvxopt_solvers.options['feastol'] = p.contol
    cvxopt_solvers.options['abstol'] = p.ftol
    cvxopt_solvers.options['reltol'] = 1e-16
    if p.iprint <= 0: 
        cvxopt_solvers.options['show_progress'] = False
        cvxopt_solvers.options['MSK_IPAR_LOG'] = 0
    xBounds2Matrix(p)
    
    f = copy(p.f).reshape(-1,1)
    
    sol = cvxopt_solvers.qp(Matrix(p.H), Matrix(p.f), Matrix(p.A), Matrix(p.b), Matrix(p.Aeq), Matrix(p.beq), solverName)
    
    p.msg = sol['status']
    if p.msg == 'optimal' :  p.istop = 1000
    else: p.istop = -100
    
    
    if sol['x'] is not None:
        p.xf = xf = asfarray(sol['x']).flatten()
        p.ff = asfarray(0.5*dot(xf, p.matMultVec(p.H, xf)) + p.dotmult(p.f, xf).sum()).flatten()
        p.duals = concatenate((asfarray(sol['y']).flatten(), asfarray(sol['z']).flatten()))
    else:
        p.ff = nan
        p.xf = nan*ones(p.n)

    def _central(self):
        ''' Return central difference function

        Member variables used
            n
            fun
            vectorized
        '''
        even_order = (np.remainder(self.n, 2) == 0)

        if self.vectorized:
            if even_order:
                f_del = lambda fun, f_x0i, x0i, h: (
                    fun(x0i + h) + fun(x0i - h)).ravel() / 2.0 - f_x0i
            else:
                f_del = lambda fun, f_x0i, x0i, h: (
                    fun(x0i + h) - fun(x0i - h)).ravel() / 2.0
        else:
            if even_order:
                f_del = lambda fun, f_x0i, x0i, h: np.asfarray(
                    [fun(x0i + h_j) + fun(x0i - h_j)
                                for h_j in h]).ravel() / 2.0 - f_x0i
            else:
                f_del = lambda fun, f_x0i, x0i, h: np.asfarray(
                    [fun(x0i + h_j) - fun(x0i - h_j)
                                        for h_j in h]).ravel() / 2.0
        return f_del

def MakeFrameMask(data,frame):
    pixelSize = data['pixelSize']
    scalex = pixelSize[0]/1000.
    scaley = pixelSize[1]/1000.
    blkSize = 512
    Nx,Ny = data['size']
    nXBlks = (Nx-1)/blkSize+1
    nYBlks = (Ny-1)/blkSize+1
    tam = ma.make_mask_none(data['size'])
    for iBlk in range(nXBlks):
        iBeg = iBlk*blkSize
        iFin = min(iBeg+blkSize,Nx)
        for jBlk in range(nYBlks):
            jBeg = jBlk*blkSize
            jFin = min(jBeg+blkSize,Ny)                
            nI = iFin-iBeg
            nJ = jFin-jBeg
            tax,tay = np.mgrid[iBeg+0.5:iFin+.5,jBeg+.5:jFin+.5]         #bin centers not corners
            tax = np.asfarray(tax*scalex,dtype=np.float32)
            tay = np.asfarray(tay*scaley,dtype=np.float32)
            tamp = ma.make_mask_none((1024*1024))
            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
                tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
            if tamp.shape:
                tamp = np.reshape(tamp[:nI*nJ],(nI,nJ))
                tam[iBeg:iFin,jBeg:jFin] = ma.mask_or(tamp[0:nI,0:nJ],tam[iBeg:iFin,jBeg:jFin])
            else:
                tam[iBeg:iFin,jBeg:jFin] = True
    return tam.T