def __call__(self, gradient, error):
        products = self.previous_gradient * gradient
        signs = sign(gradient)

        # For positive gradient parts.
        positive = (products > 0).astype('int8')
        pos_step = self.step * self.upfactor * positive
        clip(pos_step, -self.bound, self.bound)
        pos_update = self.values - signs * pos_step

        # For negative gradient parts.
        negative = (products < 0).astype('int8')
        neg_step = self.step * self.downfactor * negative
        clip(neg_step, -self.bound, self.bound)
        if error <= self.previous_error:
            # If the error has decreased, do nothing.
            neg_update = zeros(gradient.shape)
        else:
            # If it has increased, move back 2 steps.
            neg_update = self.more_prev_values
        # Set all negative gradients to zero for the next step.
        gradient *= positive

        # Bookkeeping.
        self.previous_gradient = gradient
        self.more_prev_values = self.prev_values
        self.prev_values = self.values.copy()
        self.previous_error = error

        # Updates.
        self.step[:] = pos_step + neg_step
        self.values[:] = positive * pos_update + negative * neg_update

        return self.values

 def normalize(self, sensors):
     """ The function scales the parameters to be between -1 and 1. e.g. [(-pi, pi), (0, 1), (-0.001, 0.001)] """
     assert(len(self.sensor_limits) == len(sensors))
     result = []
     for l, s in zip(self.sensor_limits, sensors):
         if not l:
             result.append(s)
         else:
             result.append((s - l[0]) / (l[1] - l[0]) * 2 - 1.0)
     if self.clipping:
         clip(result, -1, 1)
     return asarray(result)

def getBeamFluxSpline(beam,plasma,t,lim1,lim2,points = 1000):
    """ generates a spline off of the beampath.  Assumes
    that the change in flux is MONOTONIC"""

    lim = beam.norm.s
    beam.norm.s = scipy.linspace(0,lim[-1],points)
    h = time.time()
    psi = plasma.eq.rz2rmid(beam.r()[0],beam.r()[2],t) #evaluates all psi's at once
    print(time.time()-h)
    outspline = len(t)*[0]
    inspline = len(t)*[0]
    for i in xrange(t.size):
        temp = lim1
        mask = scipy.logical_and(scipy.isfinite(psi[i]),psi[i] < lim2+.02)

        try:
            minpos = scipy.argmin(psi[i][mask])
            test = psi[i][mask][minpos]
        except ValueError:
            test = lim2+.03
            
        #plt.plot(beam.x()[0][mask],psi[i][mask])
        #plt.show()
        sizer = psi[i][mask].size
        if not test > lim2:

        #plt.plot(beam.x()[0][mask][0:minpos],psi[i][mask][0:minpos],beam.x()[0][mask][minpos:],psi[i][mask][minpos:])
        #plt.show()
        #limout = scipy.insert(lim,(2,2),(beam.norm.s[mask][minpos],beam.norm.s[mask][minpos]))  # add minimum flux s for bound testing
            if lim1 < test:
                temp = test

            try:
                temp1 = scipy.clip(scipy.digitize((lim1,lim2),psi[i][mask][minpos::-1]),0,minpos)
                outspline[i] = beam.norm.s[mask][minpos::-1][temp1]
            
            except ValueError:
                tempmask = (psi[i][mask] < lim2)[0]
                outspline[i] = scipy.array([beam.norm.s[mask][minpos],beam.norm.s[mask][tempmask]])

            try:
                temp2 = scipy.clip(scipy.digitize((lim1,lim2),psi[i][mask][minpos:]),0,sizer-minpos-1)
                inspline[i] = beam.norm.s[mask][minpos:][temp2]
                
            except ValueError:
                inspline[i] = scipy.array([beam.norm.s[mask][minpos],beam.norm.s[mask][-1]])

        else:
            outspline[i] = scipy.array([[],[]])
            inspline[i] = scipy.array([[],[]])

    return (outspline,inspline)

 def normalize(self, sensors):
     """ limits is a list of 2-tuples, one tuple per parameter, giving min and max for that parameter.
         The function scales the parameters to be between -1 and 1. e.g. [(-pi, pi), (0, 1), (-0.001, 0.001)] """
     assert len(self.sensor_limits) == len(sensors)
     result = []
     for l, s in zip(self.sensor_limits, sensors):
         if not l:
             result.append(s)
         else:
             result.append((s - l[0]) / (l[1] - l[0]) * 2 - 1.0)
     if self.clipping:
         clip(result, -1, 1)
     return result

    def keyPressEvent(self, event): # reimplementation
        if event.key() == 16777234:
#            print " left arrow "
            self.Data_Display.Frame_Visualizer.frame = self.Data_Display.Frame_Visualizer.frame - 1
            self.Data_Display.Frame_Visualizer.frame = sp.clip(self.Data_Display.Frame_Visualizer.frame,0,self.Main.Data.nFrames-1)

        if event.key() == 16777236:
#            print " right arrow "
            self.Data_Display.Frame_Visualizer.frame = self.Data_Display.Frame_Visualizer.frame + 1
            self.Data_Display.Frame_Visualizer.frame = sp.clip(self.Data_Display.Frame_Visualizer.frame,0,self.Main.Data.nFrames-1)

        self.Data_Display.Frame_Visualizer.update_frame()
        self.Data_Display.Traces_Visualizer.update_vline(self.Data_Display.Frame_Visualizer.frame) # one call is enougth because this one calls the other as well

def ojf(x,s,d,override=False):
    #print "called ojf: "+str(x)
    try:
        x=x.flatten(0)
    except:
        pass
    
    xlow = [-2.,-2.]
    xupp = [2.,2.]
    
    xthis = [xlow[i]+0.5*(xin+1)*(xupp[i]-xlow[i]) for i,xin in enumerate(x)]
    hyp = [10**i for i in xthis]
    
    print hyp
    t0=time.clock()
    llk = sp.clip(GPdc.GP_LKonly(X,Y,S,D,GPdc.kernel(GPdc.MAT52,1,sp.array(hyp))).plk(pm,ps),-1e60,1e60)
    
    t1=time.clock()
    if llk<-1.:
        out = sp.log(-llk)+1.
    else:
        out = -llk
    print "--->llk: {0} {1}    t: {2}".format(llk,out,t1-t0)
    
    return [out,t1-t0]

def two_channel_to_color(im):
    """Converts a two-channel microarray image to a color image, as described in the paper associated with this 
    codebase"""
    lower = sp.percentile(im, 5)
    upper = sp.percentile(im, 98)   
    
    channel_0 = sp.clip((im[:, :, 0] - lower)/(upper - lower), 0, 1)
    channel_2 = sp.clip((im[:, :, 1] - lower)/(upper - lower), 0, 1)
    channel_1 = ((channel_0 + channel_2)/2.)
    
    im = sp.array((channel_0, channel_1, channel_2))
    im = sp.rollaxis(im, 0, 3)
    
    im = (255*im).astype(sp.uint8)    
    
    return im

def plotHeatmap(fwrap, aclass, algoparams, trials, maxsteps):
    """ Visualizing performance across trials and across time 
    (iterations in powers of 2) """
    psteps = int(log2(maxsteps)) + 1
    storesteps = [0] + [2 ** x  for x in range(psteps)]
    ls = lossTraces(fwrap, aclass, dim=trials, maxsteps=maxsteps,
                    storesteps=storesteps, algoparams=algoparams,
                    minLoss=1e-10)
            
    initv = mean(ls[0])
    maxgain = exp(fwrap.stochfun.maxLogGain(maxsteps) + 1)
    maxneggain = (sqrt(maxgain))
    
    M = zeros((psteps, trials))
    for sid in range(psteps):
        # skip the initial values
        winfactors = clip(initv / ls[sid+1], 1. / maxneggain, maxgain)
        winfactors[isnan(winfactors)] = 1. / maxneggain
        M[sid, :] = log10(sorted(winfactors))
        
    pylab.imshow(M.T, interpolation='nearest', cmap=cm.RdBu, #@UndefinedVariable
                 aspect=psteps / float(trials) / 1,  
                 vmin= -log10(maxgain), vmax=log10(maxgain),
                 )   
    pylab.xticks([])
    pylab.yticks([])
    return ls

def pso(func, nswarm, lbound, ubound, vmax, args=(), maxiter=1000, cp=2.0, cg=2.0):
    ndim = len(lbound)
    lbound = sp.asarray(lbound)
    ubound = sp.asarray(ubound)
    vmax = sp.asarray(vmax)

    # initialize the swarm
    swarm = lbound + sp.rand(nswarm, ndim)*(ubound-lbound)

    # initialize the "personal best" values
    pbestv = sp.zeros(nswarm, sp.Float)
    for i in sp.arange(nswarm):
        pbestv[i] = func(swarm[i])
    pbest = sp.array(swarm)

    # initialize the "global best" values
    gbesti = sp.argmin(pbestv)
    gbestv = pbestv[gbesti]
    gbest = pbest[gbesti]

    # initialize velocities
    velocities = 2*vmax*sp.randn(nswarm, ndim) - vmax

    for i in sp.arange(maxiter):

        values = sp.zeros(nswarm, sp.Float)
        for j in sp.arange(nswarm):
            values[j] = func(swarm[j])

        mask = values < pbestv
        mask2d = sp.repeat(mask, ndim)
        mask2d.shape = (nswarm, ndim)
        pbestv = sp.where(mask, values, pbestv)
        pbest = sp.where(mask2d, swarm, pbest)

        if sp.minimum.reduce(pbestv) < gbestv:
            gbesti = sp.argmin(pbestv)
            gbestv = pbestv[gbesti]
            gbest = pbest[gbesti]

        velocities += (cp*sp.rand()*(pbest - swarm) +
                       cg*sp.rand()*(gbest - swarm))
        velocities = sp.clip(velocities, -vmax, vmax)
        swarm += velocities
        swarm = sp.clip(swarm, lbound, ubound)
        yield gbest

    def __call__(self,x_new):
        """Find linearly interpolated y_new = <name>(x_new).

        Inputs:
          x_new -- New independent variables.

        Outputs:
          y_new -- Linearly interpolated values corresponding to x_new.
        """
        # 1. Handle values in x_new that are outside of x.  Throw error,
        #    or return a list of mask array indicating the outofbounds values.
        #    The behavior is set by the bounds_error variable.
        ## RHC -- was   x_new = atleast_1d(x_new)
        x_new_1d = atleast_1d(x_new)
        out_of_bounds = self._check_bounds(x_new_1d)
        # 2. Find where in the orignal data, the values to interpolate
        #    would be inserted.
        #    Note: If x_new[n] = x[m], then m is returned by searchsorted.
        x_new_indices = searchsorted(self.x,x_new_1d)
        # 3. Clip x_new_indices so that they are within the range of 
        #    self.x indices and at least 1.  Removes mis-interpolation
        #    of x_new[n] = x[0]
        # RHC -- changed Int to Numeric_Int to avoid name clash with numarray
        x_new_indices = clip(x_new_indices,1,len(self.x)-1).astype(Numeric_Int)
        # 4. Calculate the slope of regions that each x_new value falls in.
        lo = x_new_indices - 1; hi = x_new_indices
        
        # !! take() should default to the last axis (IMHO) and remove
        # !! the extra argument.
        x_lo = take(self.x,lo,axis=self.interp_axis)
        x_hi = take(self.x,hi,axis=self.interp_axis)
        y_lo = take(self.y,lo,axis=self.interp_axis)
        y_hi = take(self.y,hi,axis=self.interp_axis)
        slope = (y_hi-y_lo)/(x_hi-x_lo)
        # 5. Calculate the actual value for each entry in x_new.
        y_new = slope*(x_new_1d-x_lo) + y_lo 
        # 6. Fill any values that were out of bounds with NaN
        # !! Need to think about how to do this efficiently for 
        # !! mutli-dimensional Cases.
        yshape = y_new.shape
        y_new = y_new.flat
        new_shape = list(yshape)
        new_shape[self.interp_axis] = 1
        sec_shape = [1]*len(new_shape)
        sec_shape[self.interp_axis] = len(out_of_bounds)
        out_of_bounds.shape = sec_shape
        new_out = ones(new_shape)*out_of_bounds
        putmask(y_new, new_out.flat, self.fill_value)
        y_new.shape = yshape
        # Rotate the values of y_new back so that they correspond to the
        # correct x_new values.
        result = swapaxes(y_new,self.interp_axis,self.axis)
        try:
            len(x_new)
            return result
        except TypeError:
            return result[0]
        return result

 def _boltzmannProbs(qvalues, temperature=1.):
     if temperature == 0:
         tmp = zeros(len(qvalues))        
         tmp[r_argmax(qvalues)] = 1.
     else:
         tmp = qvalues / temperature            
         tmp -= max(tmp)        
         tmp = exp(clip(tmp, -20, 0))
     return tmp / sum(tmp)

 def performAction(self, action):
     """ a filtered mapping towards performAction of the underlying environment. """
     # scaling
     self.incStep()
     action = (action + 1.0) / 2.0 * self.dif + self.env.fraktMin * self.env.dists[0]
     #Clipping the maximal change in actions (max force clipping)
     action = clip(action, self.action - self.maxSpeed, self.action + self.maxSpeed)
     EpisodicTask.performAction(self, action)
     self.action = action.copy()

 def getReward(self):
     # calculate reward and return reward
     if self.count < 800:
         return 0.0
     else:
         reward = self.env.getSensorByName('SpecificBodyPositionSensor8')[1] / float(self.epiLen - 800) #reward is hight of head
         #to prevent jumping reward can't get bigger than head position while standing absolut upright
         reward = clip(reward, -14.0, 4.0)
         return reward

def late_filling(target, pressure='pore.pressure',
                 Pc_star='pore.pc_star',
                 Swp_star=0.2, eta=3):
    r"""
    Calculates the fraction of a pore or throat filled with invading fluid
    based on the capillary pressure in the invading phase.  The invading phase
    volume is calculated from:

        .. math::
            S_{nwp} = 1 - S_{wp}^{*} (P^{*}/P_{c})^{\eta}

    Parameters
    ----------
    pressure : string
        The capillary pressure in the non-wetting phase (Pc > 0).

    Pc_star : string
        The minimum pressure required to create an interface within the pore
        body or throat.  Typically this would be calculated using the Washburn
        equation.

    Swp_star : float
        The residual wetting phase in an invaded pore or throat at a pressure
        of ``pc_star``.

    eta : float
        Exponent controlling the rate at which wetting phase is displaced with
        increasing pressure.

    Returns
    -------
    An array containing the fraction of each pore or throat that would be
    filled with non-wetting phase at the given phase pressure.  This does not
    account for whether or not the element is actually invaded, which requires
    a percolation algorithm of some sort.

    """
    element = pressure.split('.')[0]
    network = target.project.network
    phase = target.project.find_phase(target)
    pc_star = phase[Pc_star]
    Pc = phase[pressure]
    # Remove any 0's from the Pc array to prevent numpy div by 0 warning
    Pc = sp.maximum(Pc, 1e-9)
    Swp = Swp_star*((pc_star/Pc)**eta)
    values = sp.clip(1 - Swp, 0.0, 1.0)
    # Now map element onto target object
    if element == 'throat':
        Ts = network.map_throats(throats=target.Ts, origin=target)
        values = values[Ts]
    else:
        Ps = network.map_pores(pores=target.Ps, origin=target)
        values = values[Ps]
    return values

def pointchargePot(x,y,charge=1,scale=1):
    from scipy import sqrt,pi,clip
    size_x = x.max()
    size_y = y.max()
    Vbottom = 0
    x = x-size_x/2
    left_charge = 1/sqrt(x**2+y**2)
    right_charge = 1/sqrt(x**2+(y-size_y)**2)
    V = Vbottom +p.q*charge/(p.a*4*pi*p.eps0*p.epsr)*(left_charge+right_charge)
    V = clip(V,0,scale)
    return V

def calcProfilV(self,xy):
    """renvoie les valeurs des vitesses sur une section"""
    vxvy = self.getMfVitesse()
    grd  = self.parent.aquifere.getFullGrid()
    x0,y0,dx,dy,nx,ny = grd['x0'],grd['y0'],grd['dx'],grd['dy'],grd['nx'],grd['ny']
    x,y = zip(*xy)
    xl0, xl1 = x[:2]
    yl0, yl1 = y[:2]
    dd = min(dx,dy)*.95;dxp, dyp = xl1-xl0, yl1-yl0
    ld = max(ceil(abs(dxp/dx)),ceil(abs(dyp/dy)))
    ld = int(ld+1); ddx = dxp/ld; ddy = dyp/ld
    xp2 = xl0+arange(ld+1)*ddx
    yp2 = yl0+arange(ld+1)*ddy
    ix = floor((xp2-x0)/dx);ix=clip(ix.astype(int),0,nx-1)
    iy = floor((yp2-y0)/dy);iy=clip(iy.astype(int),0,ny-1)
    vx = take(ravel(vxvy[0]),iy*nx+ix)
    vy = take(ravel(vxvy[1]),iy*nx+ix)
    V = sqrt(vx**2+vy**2)
    cu = sqrt((xp2-xp2[0])**2+(yp2-yp2[0])**2)
    return [cu,V]

 def scroll_event(self,event):
     """ changes width of slice """
     if event.button == 'up':
         self.width += 2
     if event.button == 'down':
         self.width -= 2
         
     self.width = sp.clip(self.width,1,self.nPlaces)
     self.xs = self.calc_x(self.pos,self.width)
     self.Rect.set_xy((self.xs[0] ,0))
     self.Rect.set_width(self.width)
     self.update()

    def denormalize(self, actors):
        """ The function scales the parameters from -1 and 1 to the given interval (min, max) for each actor. """
        assert(len(self.actor_limits) == len(actors))
        result = []
        for l, a in zip(self.actor_limits, actors):
            if not l:
                result.append(a)
            else:
                r = (a + 1.0) / 2 * (l[1] - l[0]) + l[0]
                if self.clipping:
                    r = clip(r, l[0], l[1])
                result.append(r)

        return result

def get_color(listing, listings, f, cm):
    price = f(listing)
    prices = [f(l) for l in listings]

    lower = sp.percentile(prices, 10)
    upper = sp.percentile(prices, 90)

    relative_price = (price - lower)/(upper - lower)
    color = cm(sp.clip(relative_price, 0, 1))

    is_dark = sum(color[:3])/4 < 0.4
    background_color = tuple([int(255*c) for c in color[:3]])
    text_color = (230, 230, 230) if is_dark else (50, 50, 50)

    return background_color, text_color

    def denormalize(self, actors):
        """ limits is a list of 2-tuples, one tuple per parameter, giving min and max for that parameter.
            The function scales the parameters from -1 and 1 to the given interval (min, max) for each actor. """
        assert len(self.actor_limits) == len(actors)
        result = []
        for l, a in zip(self.actor_limits, actors):
            if not l:
                result.append(a)
            else:
                r = (a + 1.0) / 2 * (l[1] - l[0]) + l[0]
                if self.clipping:
                    r = clip(r, l[0], l[1])
                result.append(r)

        return result