def solve(self):
    """
    """
    s    = "::: solving density, overburden stress, and grain radius :::"
    text = colored(s, 'cyan')
    print text
    
    firn   = self.firn
    config = self.config

    # newton's iterative method :
    solve(self.delta == 0, firn.U, bcs=self.bcs, J=self.J, 
          solver_parameters=config['enthalpy']['solver_params'])
    firn.rho, firn.sigma, firn.r = firn.U.split(True)

    rhop = firn.rho.vector().array()

    # update kc term in drhodt :
    # if rho >  550, kc = kcHigh
    # if rho <= 550, kc = kcLow
    # with parameterizations given by ligtenberg et all 2011
    A                   = firn.rhoi/firn.rhow * 1e3 * firn.adot
    rhoCoefNew          = ones(firn.n)
    rhoHigh             = where(rhop >  550)[0]
    rhoLow              = where(rhop <= 550)[0]
    rhoCoefNew[rhoHigh] = firn.kcHh * (2.366 - 0.293*ln(A))
    rhoCoefNew[rhoLow]  = firn.kcLw * (1.435 - 0.151*ln(A))
    firn.assign_variable(firn.rhoCoef, rhoCoefNew)
    firn.assign_variable(firn.drhodt,  project(self.drhodt))
    
    #firn.assign_variable(firn.rho, rhop)
    firn.print_min_max(firn.rho,   'rho')
    firn.print_min_max(firn.sigma, 'sigma')
    firn.print_min_max(firn.r,     'r^2')

 def v(current_data):
     from pylab import where,sqrt
     q = current_data.q
     h = q[0,:]
     dry_tol = 0.001
     u = where(h>dry_tol, q[1,:]/h, 0.)
     v = where(h>dry_tol, q[2,:]/h, 0.)
     return v

def fwhm(x, y):
	hm = pl.amax(y/2.0);
	y_diff = pl.absolute(y-hm);
	y_diff_sorted = pl.sort(y_diff);
	i1 = pl.where(y_diff==y_diff_sorted[0]);
	i2 = pl.where(y_diff==y_diff_sorted[1]);
	fwhm = pl.absolute(x[i1]-x[i2]);
	return hm, fwhm

 def speed(current_data):
     from pylab import where,sqrt
     q = current_data.q
     h = q[0,:]
     dry_tol = 0.001
     u = where(h>dry_tol, q[1,:]/h, 0.)
     v = where(h>dry_tol, q[2,:]/h, 0.)
     s = sqrt(u**2 + v**2)
     return s

 def shadowing(self):
     "Select the shadowed antennas from the FLAG column and return the index of the shadowed measurement and the percentage of shadowing "
     
     indexFlag=pl.concatenate((pl.where(self.f==1)[0],pl.where(self.ff[0,0,])[0],pl.where(self.ff[1,0,])[0]))
     indexNoFlag=pl.concatenate((pl.where(self.f==0)[0],pl.where(self.ff[0,0,]==False)[0],pl.where(self.ff[1,0,]==False)[0]))
     
     Ntot=len(indexFlag)+len(indexNoFlag)
     fractionShadow=100.*len(indexFlag)/Ntot
     
     return(indexFlag,fractionShadow)

 def gauge_speed(current_data):
     # different than speed function because q is function of time, not
     # x,y at the gauges.
     from numpy import where, sqrt
     h = q[0,:]
     dry_tol = 0.001
     u = where(h>dry_tol, q[1,:]/h, 0.)
     v = where(h>dry_tol, q[2,:]/h, 0.)
     s = sqrt(u**2 + v**2)
     return s

def fbank(signal,samplerate=16000,winlen=0.025,winstep=0.01,
          nfilt=26,nfft=512,lowfreq=0,highfreq=None,preemph=0.97):
    """Compute Mel-filterbank energy features from an audio signal.

    :param signal: the audio signal from which to compute features. Should be an N*1 array
    :param samplerate: the samplerate of the signal we are working with.
    :param winlen: the length of the analysis window in seconds. Default is 0.025s (25 milliseconds)    
    :param winstep: the step between successive windows in seconds. Default is 0.01s (10 milliseconds)    
    :param nfilt: the number of filters in the filterbank, default 26.
    :param nfft: the FFT size. Default is 512.
    :param lowfreq: lowest band edge of mel filters. In Hz, default is 0.
    :param highfreq: highest band edge of mel filters. In Hz, default is samplerate/2
    :param preemph: apply preemphasis filter with preemph as coefficient. 0 is no filter. Default is 0.97. 
    :returns: 2 values. The first is a numpy array of size (NUMFRAMES by nfilt) containing features. Each row holds 1 feature vector. The
        second return value is the energy in each frame (total energy, unwindowed)
    """          
    highfreq= highfreq or samplerate/2
    print "preemph %s"%(preemph)
    signal = sigproc.preemphasis(signal,preemph)
    frames = sigproc.framesig(signal, winlen*samplerate, winstep*samplerate)
    matchframes(frames[0], frames[1])
    pspec = sigproc.powspec(frames,nfft)
    energy = pylab.sum(pspec,1) # this stores the total energy in each frame
    energy = pylab.where(energy == 0, pylab.finfo(float).eps, energy) # if energy is zero, we get problems with log
    fb = get_filterbanks(nfilt, nfft, samplerate, lowfreq, highfreq)
    print "len(fb) %s"%(len(fb))
    colour = "k-"
    for i in range(len(fb)):
        if colour == "k-":
            colour = "r-"
        else:
            colour = "k-"
        startedplot = False
        midpoint = 0
        for j in range(len(fb[i])):
            if fb[i][j] > 0:
                if startedplot == False:
                    startedplot = j
                if j > 0:
                    pylab.plot([j-1, j], [fb[i][j-1], fb[i][j]], colour)
                    if fb[i][j] == 1.0:
                        midpoint = j
            else:
                if not startedplot == False:
                    pylab.plot([j-1, j], [fb[i][j-1], 0], colour)
                    try:
                        print "slope to midpoint %.3f, slope from midpoint %.3f"%(1.0/float(midpoint-startedplot), 1.0/float(midpoint-j+1))
                    except:
                        pass
                    break
    pylab.show()
    feat = pylab.dot(pspec, fb.T) # compute the filterbank energies
    feat = pylab.where(feat == 0, pylab.finfo(float).eps, feat) # if feat is zero, we get problems with log
    return feat, energy

 def int_peak(self,fitrange=None, intrange=None, normalize=False, plot=False, npoints=10):
     """
     Fits a linear background, subtracts the background, and integrates. Intended to be used for integrating peaks.
     
     wavelen : list
         list of wavelengths in nm. Can be sorted from low to high or high to low
     lum : list
         list of luminescence
     fitrange : 2-element list, optional
         Defaults to the span of the data. Input: [low nm, high nm]
     intrange : 2-element list, optional
         Defaults to the span of the data or fitrange (if given). Input: [low nm, high nm]
     normalize : boolean, optional
         Default is False
     plot : boolean, optional
         Default is False. Plots the original data, the linear background, and the data with the background subtracted
     npoints : int
         Default is 10. Number of points above and below the given fitrange point to average over.
     """
     if fitrange is None:
         fitindex=[0+npoints/2, len(self._wavelen)-1-npoints/2]
     else:
         fitindex=[0, 0]
         fitindex[0]=py.where(self._wavelen>fitrange[0])[0][0]
         fitindex[1]=py.where(self._wavelen>fitrange[1])[0][0]
     
     wavelenfit=py.concatenate((self._wavelen[fitindex[0]-npoints/2:fitindex[0]+npoints/2], 
                            self._wavelen[fitindex[1]-npoints/2:fitindex[1]+npoints/2]))
     lumfit=py.concatenate((self._lum[fitindex[0]-npoints/2:fitindex[0]+npoints/2], 
                         self._lum[fitindex[1]-npoints/2:fitindex[1]+npoints/2]))
     linearfit = py.polyfit(wavelenfit, lumfit, 1)
     linear_bg = py.polyval( linearfit, self._wavelen[fitindex[0]:fitindex[1]+1] )
     wavelen_bg = self._wavelen[fitindex[0]:fitindex[1]+1].copy()
     lum_bg = self._lum[fitindex[0]:fitindex[1]+1].copy()
     lum_bg -= linear_bg
     
     if plot is True:
         py.plot(self._wavelen,self._lum,'k')
         py.plot(wavelen_bg,linear_bg,'k:')
         py.plot(wavelen_bg,lum_bg,'r')
         py.show()
     
     intindex=[0,0]
     if intrange is None:
         wavelen_int = wavelen_bg
         lum_int = lum_bg  
     else:
         intindex[0]=py.where(wavelen_bg>intrange[0])[0][0]
         intindex[1]=py.where(wavelen_bg>intrange[1])[0][0]    
         wavelen_int = wavelen_bg[intindex[0]:intindex[1]+1]
         lum_int = lum_bg[intindex[0]:intindex[1]+1]
     
     peak_area = py.trapz(lum_int, x=wavelen_int)
     return peak_area

 def speed(current_data):
     from pylab import where,sqrt
     q = current_data.q
     h = q[0,:]
     dry_tol = 0.001
     u = where(h>dry_tol, q[1,:]/h, 0.)
     v = where(h>dry_tol, q[2,:]/h, 0.)
     s = sqrt(u**2 + v**2)
     #s = s / sqrt(9.81/0.97)  # so comparable to eta
     s = s / 10.  # as in Figure 5 of paper
     return s

 def save_sonogram(
     self,
     replace=False,
     n_fft=settings.N_FFT,
     min_freq=settings.MIN_FREQ,
     max_freq=settings.MAX_FREQ,
     dpi=100,
     width=1000,
     height=350,
     max_framerate=settings.MAX_FRAMERATE,
 ):
     filename = self.get_sonogram_name()
     name = os.path.join(settings.SONOGRAM_DIR, filename)
     path = os.path.join(settings.MEDIA_ROOT, name)
     try:
         if not os.path.exists(path):
             replace = True
     except (ValueError, SuspiciousOperation, AttributeError):
         replace = True
     if replace:
         audio, framerate = self.get_audio(max_framerate=max_framerate)
         Pxx, freqs, bins, im = specgram(audio, NFFT=n_fft, Fs=framerate)
         f = where(logical_and(freqs > min_freq, freqs <= max_freq))[0]
         Pxx[where(Pxx > percentile(Pxx[f].flatten(), 99.99))] = percentile(Pxx[f].flatten(), 99.99)
         Pxx[where(Pxx < percentile(Pxx[f].flatten(), 0.01))] = percentile(Pxx[f].flatten(), 0.01)
         clf()
         fig = figure(figsize=(float(width) / dpi, float(height) / dpi), dpi=dpi)
         imshow(
             flipud(10 * log10(Pxx[f,])),
             extent=(bins[0], bins[-1], freqs[f][0], freqs[f][-1]),
             aspect="auto",
             cmap=cm.gray,
         )
         gca().set_ylabel("Frequency (Hz)")
         gca().set_xlabel("Time (s)")
         axis_pixels = gca().transData.transform(np.array((gca().get_xlim(), gca().get_ylim())).T)
         st, created = SonogramTransform.objects.get_or_create(
             n_fft=n_fft,
             framerate=framerate,
             min_freq=min_freq,
             max_freq=max_freq,
             duration=self.duration,
             width=width,
             height=height,
             dpi=dpi,
             top_px=max(axis_pixels[:, 1]),
             bottom_px=min(axis_pixels[:, 1]),
             left_px=min(axis_pixels[:, 0]),
             right_px=max(axis_pixels[:, 0]),
         )
         savefig(open(path, "wb"), format="jpg", dpi=dpi)
         sonogram, created = Sonogram.objects.get_or_create(snippet=self, transform=st, path=name)
         close()

def index_approx(v, value):
    """Find index of element closest to value"""
    if isinstance(v, (list, tuple)):
        a = min([abs(i-value) for i in v])
        try:
            return v.index(value + a)
        except ValueError:
            return v.index(value - a)
    elif isinstance(v, pylab.ndarray):
        a = min(abs(v-value))
        if pylab.where(v == value + a)[0].size:
            return pylab.where(v == value + a)[0][0]
        else:
            return pylab.where(v == value - a)[0][0]

 def __setitem__(self, wavelength, intensity):
     index, = pylab.where(self.wavelengths == wavelength)
     if pylab.any(index.shape):
         self.intensities[index] = intensity
     else:
         index, = pylab.where(self.wavelengths < wavelength)
         if pylab.any(index.shape):
             self.wavelengths = pylab.insert(self.wavelengths, index[-1] + 1,
                                             wavelength)
             self.intensities = pylab.insert(self.intensities, index[-1] + 1,
                                             intensity)
         else:
             self.wavelengths = pylab.insert(self.wavelengths, 0, wavelength)
             self.intensities = pylab.insert(self.intensities, 0, intensity)

def compare(features, targets, columns):
    q0mask = pyl.where(targets == 0)[0]
    q1mask = pyl.where(targets == 1)[0]
    q2mask = pyl.where(targets == 2)[0]

    pyl.scatter(features[:, columns[0]][q0mask],
                features[:, columns[1]][q0mask], c='r', label='0')
    pyl.scatter(features[:, columns[0]][q1mask],
                features[:, columns[1]][q1mask], c='g', label='1')
    pyl.scatter(features[:, columns[0]][q2mask],
                features[:, columns[1]][q2mask], c='b', label='2')

    pyl.legend(loc='best')
    pyl.show()

 def gs(current_data):
     q = current_data.q
     # different than speed function because q is function of time, not
     # x,y at the gauges.
     from numpy import where, sqrt
     h = q[0,:]
     #print('shape of h ' +  str(h.shape))
     dry_tol = 0.001
     u = where(h>dry_tol, q[1,:]/h, 0.)
     v = where(h>dry_tol, q[2,:]/h, 0.)
     ssq = sqrt(u*u+v*v)
     #s = sqrt(u**2 + v**2)
     s = sqrt(ssq)
     return ssq

def preprocessing(input_image, scale, factor):
    """Given an image object, return the difference of gaussian matrix
    
    """
    var1 = scale / sqrt(2)
    var2 = scale * sqrt(2)
    image_array = (asarray(input_image)).astype(float64)/255. * scale
    threshold = image_array.max() * factor
    image_array[where(image_array < threshold)] = 0
    image_array[where(image_array >= threshold)] = 255
    
    image_gaussian = filters.gaussian_filter(image_array, sigma=(var1,var1)) - \
                     filters.gaussian_filter(image_array, sigma=(var2,var2))
    
    return image_gaussian    

def fwhm_2gauss(x, y, dx=0.001):
	'''
	Finds the FWHM for the profile y(x), with accuracy dx=0.001
	Uses a 2-Gauss 1D fit.
	'''
	popt, pcov = curve_fit(gauss2, x, y);
	xx = pl.arange(pl.amin(x), pl.amax(x)+dx, dx);
	ym = gauss2(xx, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
	hm = pl.amax(ym/2.0);
	y_diff = pl.absolute(ym-hm);
	y_diff_sorted = pl.sort(y_diff);
	i1 = pl.where(y_diff==y_diff_sorted[0]);
	i2 = pl.where(y_diff==y_diff_sorted[1]);
	fwhm = pl.absolute(xx[i1]-xx[i2]);
	return hm, fwhm, xx, ym

 def vvel(current_data):
     from numpy import where, sqrt
     q = current_data.q
     h = q[0,:]
     dry_tol = 0.001
     v = where(h>dry_tol, q[2,:]/h, 0.)
     return v

 def activate(self, inputs):
   assert len(inputs) == NN.ni, 'incorrect number of inputs'
   self.ai[:] = np.tanh(inputs).copy()
   self.ah[:] = np.tanh(np.sum(self.wi.T * self.ai, axis=1))
   self.ao[:] = np.sum(self.wo.T * self.ah, axis=1) + 0.5
   self.ao[:] = where(self.ao == max(self.ao), 1.0, 0.0)
   return self.ao.copy()

def specgram_demo():
   '''
   the demo in matplotlib. But calls
   interactive.specgram
   '''
   from pylab import arange, sin, where, logical_and, randn, pi

   dt = 0.0005
   t = arange(0.0, 20.0, dt)
   s1 = sin(2*pi*100*t)
   s2 = 2*sin(2*pi*400*t)

   # create a transient "chirp"
   mask = where(logical_and(t>10, t<12), 1.0, 0.0)
   s2 = s2 * mask

   # add some noise into the mix
   nse = 0.01*randn(len(t))

   x = s1 + s2 + nse # the signal
   NFFT = 1024       # the length of the windowing segments
   Fs = int(1.0/dt)  # the sampling frequency

   from ifigure.interactive import figure, specgram, nsec, plot, isec, clog, hold

   figure()
   hold(True)
   nsec(2)
   isec(0)
   plot(t, x)
   isec(1)
   specgram(x, NFFT=NFFT, Fs=Fs, noverlap=900)
   clog()

    def addDataVectorAccessor(self, data_vector_accessor):
        self.__data_vectors_accessors__.append(data_vector_accessor)

        _sum = pl.sum(data_vector_accessor.signal)
        _min = pl.amin(data_vector_accessor.signal)
        _max = pl.amax(data_vector_accessor.signal)

        if self.__minimal_signal__ == None:
            self.__minimal_signal__ = _sum
            self.__minimal_data_vector_accessor__ = data_vector_accessor

            self.__min_signal__ = _min
            self.__max_signal__ = _max

        if _sum < self.__minimal_signal__:
            self.__minimal_data_vector_accessor__ = data_vector_accessor
            self.__minimal_signal__ = _sum

        if _min < self.__min_signal__:
            self.__min_signal__ = _min

        if _max > self.__max_signal__:
            self.__max_signal__ = _max

        #collects unique annotations (>0) as a set
        if not data_vector_accessor.annotation == None:
            unique_annotations = pl.unique(data_vector_accessor.annotation[
                                pl.where(data_vector_accessor.annotation > 0)])
            if len(unique_annotations) > 0:
                #union of sets
                self.__unique_annotations__ |= set(unique_annotations)