def produce_heatmap(model, every = True, save = False):
    col_label = range(28)
    row_label = range(28)
    if every:
        for i in range(10):
            plt.pcolor(np.flipud(model[i]))
            plt.xticks(col_label)
            plt.yticks(row_label)
            plt.axis('off')
            plt.title("HeatMap for %d" % (i))
            cb = plt.colorbar()
            cb.set_label("Frequency")
            if save:
                plt.savefig('imgs/%d.png' % (i), bbox_inches='tight')
            else:
                plt.show()
            plt.close()
    else:
        plt.pcolor(np.flipud(model))
        plt.xticks(col_label)
        plt.yticks(row_label)
        plt.axis('off')
        cb = plt.colorbar()
        cb.set_label("Frequency")
        if save:
            plt.savefig('imgs/temp.png', bbox_inches='tight')
        else:
            plt.show()
        plt.close()

    def get_output(self, idx):
        img_id=idx/self.pertnum
        pert_id=idx%self.pertnum
        rot_id=pert_id%self.param['rotate']
        off_id=pert_id/self.param['rotate']
        [h, w]=self.output[img_id].shape
        [dy, dx]=self.get_offset(h, w, off_id)
        dy+=self.param['mrgsize']
        dx+=self.param['mrgsize']
        res=self.output[img_id][dy:dy+self.param['outsize'], dx:dx+self.param['outsize']]
        #res=np.rot90(res) #rotate 90
        if rot_id==1:
            res=np.fliplr(res)
        elif rot_id==2:
            res=np.flipud(res).T
        elif rot_id==3:
            res=res.T
        elif rot_id==4:
            res=np.fliplr(res).T
        elif rot_id==5:
            res=np.flipud(res)
        elif rot_id==6:
            res=np.rot90(res,2)
        elif rot_id==7:
            res=np.rot90(res,2).T

        return res

def rotate_data(bg, overlay, slices_list, axis_name, shape):
    # Rotate the data as required
    # Return the rotated data, and an updated slice list if necessary
    if axis_name == 'axial':
        # Align so that right is right
        overlay = np.rot90(overlay)
        overlay = np.fliplr(overlay)
        bg = np.rot90(bg)
        bg = np.fliplr(bg)
    
    elif axis_name == 'coronal':
        overlay = np.rot90(overlay)
        bg = np.rot90(bg)
        overlay = np.flipud(np.swapaxes(overlay, 0, 2))
        bg = np.flipud(np.swapaxes(bg, 0, 2))
        slices_list[1] = [ shape - n - 3 for n in slices_list[1] ] 
        
    elif axis_name == 'sagittal':
        overlay = np.flipud(np.swapaxes(overlay, 0, 2))
        bg = np.flipud(np.swapaxes(bg, 0, 2))
    
    else:
        print '\n************************'
        print 'ERROR: data could not be rotated\n'
        parser.print_help()
        sys.exit()
    
    return bg, overlay, slices_list

def get_image_quadrants(IM, reorient=False):
    """
    Given an image (m,n) return its 4 quadrants Q0, Q1, Q2, Q3
    as defined in abel.hansenlaw.iabel_hansenlaw

    Parameters:
      - IM: 1D or 2D array
      - reorient: reorient image as required by abel.hansenlaw.iabel_hansenlaw
    """
    IM = np.atleast_2d(IM)

    n, m = IM.shape

    n_c = n//2 + n%2
    m_c = m//2 + m%2

    # define 4 quadrants of the image
    # see definition in abel.hansenlaw.iabel_hansenlaw
    Q1 = IM[:n_c, :m_c]
    Q2 = IM[-n_c:, :m_c]
    Q0 = IM[:n_c, -m_c:]
    Q3 = IM[-n_c:, -m_c:]

    if reorient:
        Q1 = np.fliplr(Q1)
        Q3 = np.flipud(Q3)
        Q2 = np.fliplr(np.flipud(Q2))

    return Q0, Q1, Q2, Q3

def generateWedge(a, b, n, type, plotShape=False):
    # Define x and y coordinate
    xFrontTop_ = np.linspace(0, a / 2.0, n);                    xFrontTop = xFrontTop_
    xBackTop_ = np.linspace(a / 2.0, a, n);                     xBackTop = xBackTop_[1:]
    xBackBottom_ = np.flipud(np.linspace(a / 2.0, a, n));       xBackBottom = xBackBottom_[1:]
    xFrontBottom_ = np.flipud(np.linspace(0, a / 2.0, n));      xFrontBottom = xFrontBottom_[1:-1]

    yFrontTop_ = b / a * xFrontTop_;                            yFrontTop = yFrontTop_
    yBackTop_ = -b / a * (xBackTop_ - a / 2.0) + b / 2;         yBackTop = yBackTop_[1:]
    yBackBottom_ = b / a * (xBackBottom_ - a / 2.0) - b / 2;    yBackBottom = yBackBottom_[1:]
    yFrontBottom_ = -b / a * xFrontBottom_;                     yFrontBottom = yFrontBottom_[1:-1]

    if type == 'twoSided':
        # x = np.concatenate((xFrontTop, xBackTop, xBackBottom, xFrontBottom))
        # y = np.concatenate((yFrontTop, yBackTop, yBackBottom, yFrontBottom))
        xTop = np.concatenate((xFrontTop, xBackTop))
        xBottom = np.concatenate((xBackBottom, xFrontBottom))
        yTop = np.concatenate((yFrontTop, yBackTop))
        yBottom = np.concatenate((yBackBottom, yFrontBottom))
    elif type == 'oneSided':
        x = np.concatenate((xFrontTop, xBackTop))
        y = np.concatenate((yFrontTop, yBackTop))

    if plotShape:
        plt.figure()
        plt.plot(x, y, 'k')
        plt.xlim([-5, 14])
        plt.ylim([-5, 5])
        plt.show()

    return (x, y)

    def update_lambda(self, sstats, word_list, opt_o):
        self.m_status_up_to_date = False
        # rhot will be between 0 and 1, and says how much to weight
        # the information we got from this mini-chunk.
        rhot = self.m_scale * pow(self.m_tau + self.m_updatect, -self.m_kappa)
        if rhot < rhot_bound:
            rhot = rhot_bound
        self.m_rhot = rhot

        # Update appropriate columns of lambda based on documents.
        self.m_lambda[:, word_list] = self.m_lambda[:, word_list] * (1 - rhot) + \
            rhot * self.m_D * sstats.m_var_beta_ss / sstats.m_chunksize
        self.m_lambda_sum = (1 - rhot) * self.m_lambda_sum + \
            rhot * self.m_D * np.sum(sstats.m_var_beta_ss, axis=1) / sstats.m_chunksize

        self.m_updatect += 1
        self.m_timestamp[word_list] = self.m_updatect
        self.m_r.append(self.m_r[-1] + np.log(1 - rhot))

        self.m_varphi_ss = (1.0 - rhot) * self.m_varphi_ss + rhot * \
            sstats.m_var_sticks_ss * self.m_D / sstats.m_chunksize

        if opt_o:
            self.optimal_ordering()

        ## update top level sticks
        self.m_var_sticks[0] = self.m_varphi_ss[:self.m_T - 1] + 1.0
        var_phi_sum = np.flipud(self.m_varphi_ss[1:])
        self.m_var_sticks[1] = np.flipud(np.cumsum(var_phi_sum)) + self.m_gamma

 def setUp(self):
     """
     This generates a minimum data-set to be used for the regression.
     """
     # Test A: Generates a data set assuming b=1 and N(m=4.0)=10.0 events
     self.dmag = 0.1
     mext = np.arange(4.0, 7.01, 0.1)
     self.mval = mext[0:-1] + self.dmag / 2.0
     self.bval = 1.0
     self.numobs = np.flipud(
         np.diff(np.flipud(10.0 ** (-self.bval * mext + 8.0))))
     # Test B: Generate a completely artificial catalogue using the
     # Gutenberg-Richter distribution defined above
     numobs = np.around(self.numobs)
     size = int(np.sum(self.numobs))
     magnitude = np.zeros(size)
     lidx = 0
     for mag, nobs in zip(self.mval, numobs):
         uidx = int(lidx + nobs)
         magnitude[lidx:uidx] = mag + 0.01
         lidx = uidx
     year = np.ones(size) * 1999
     self.catalogue = Catalogue.make_from_dict(
         {'magnitude': magnitude, 'year': year})
     # Create the seismicity occurrence calculator
     self.aki_ml = AkiMaxLikelihood()

def plot_kmeans(dist_m, shape_port, dat_port, dat_star, dat_kclass, ft, humfile, sonpath, base, p):

   Zdist = dist_m[shape_port[-1]*p:shape_port[-1]*(p+1)]
   extent = shape_port[1] 

   levels = [0.5,0.75,1.25,1.5,1.75,2,3]

   fig = plt.figure()
   plt.subplot(2,1,1)
   ax = plt.gca()
   plt.imshow(np.vstack((np.flipud(dat_port), dat_star)), cmap='gray',extent=[min(Zdist), max(Zdist), -extent*(1/ft), extent*(1/ft)],origin='upper')
   
   CS = plt.contourf(np.flipud(dat_kclass), levels, alpha=0.4, extent=[min(Zdist), max(Zdist), -extent*(1/ft), extent*(1/ft)],origin='upper', cmap='YlOrRd', vmin=0.5, vmax=3)   
   plt.ylabel('Horizontal distance (m)')
   plt.xlabel('Distance along track (m)')
   plt.axis('tight')

   try:
      divider = make_axes_locatable(ax)
      cax = divider.append_axes("right", size="5%", pad=0.05)
      plt.colorbar(CS, cax=cax)
   except:
      plt.colorbar()

   custom_save(sonpath,base+'class_kmeans'+str(p))
   del fig

def data_reorient_udrot(d, pl, ud=None, rot90=None):
    """
    """
    print >>sys.stderr, "Reorienting...",
    if (not ud is None) and (not rot90 is None):
        if ud > 0:
            d = np.flipud(d)
        if rot90 > 0:
            d = np.rot90(d, rot90)
        print >>sys.stderr, "done."
        return d
        
    ud = pl.get_val('REORIENT_CCD_FLIPUD')
    if ud == None:
        ud = pl.get_val('CORRECTION_FLIPUD')
    ud = int(ud)
    rot90 = pl.get_val('REORIENT_CCD_ROT90')
    if rot90 == None:
        rot90 = pl.get_val('CORRECTION_ROTATION')
    rot90 = int(rot90)
    if ud > 0:
        d = np.flipud(d)
    if rot90 > 0:
        d = np.rot90(d, rot90)
    print >>sys.stderr, "done."
    return d

    def to_unit(self, unit):
        """Return wavelength and transmission in new wavelength units.

        If the requested units are the same as the current units, self is
        returned.

        Parameters
        ----------
        unit : `~astropy.units.Unit` or str
            Target wavelength unit.

        Returns
        -------
        wave : `~numpy.ndarray`
        trans : `~numpy.ndarray`
        """

        if unit is u.AA:
            return self.wave, self.trans

        d = u.AA.to(unit, self.wave, u.spectral())
        t = self.trans
        if d[0] > d[-1]:
            d = np.flipud(d)
            t = np.flipud(t)
        return d, t

def schedule_jobs(job_matrix, mode='ratio'):
    job_weight = job_matrix[:, 0]
    job_length = job_matrix[:, 1]
    if mode == 'difference':
        job_metric = numpy.subtract(job_weight, job_length)
        job_matrix = numpy.vstack((job_weight, job_length, job_metric)).T
        job_matrix = job_matrix[job_matrix[:, 2].argsort()]
        job_matrix = numpy.flipud(job_matrix)

        metric_value_freq = Counter(job_matrix[:, 2])

        # Break ties between jobs using weight (i.e. higher weight goes first)
        i = 0
        for key in sorted(metric_value_freq.keys(), reverse=True):
            subarray = job_matrix[i:i + metric_value_freq[key]]
            sorted_subarray = subarray[subarray[:, 0].argsort()]
            sorted_subarray = numpy.flipud(sorted_subarray)
            job_matrix[i:i + metric_value_freq[key]] = sorted_subarray
            i += metric_value_freq[key]

    elif mode == 'ratio':
        job_metric = numpy.divide(job_weight, job_length)

        job_matrix = numpy.vstack((job_weight, job_length, job_metric)).T
        job_matrix = job_matrix[job_matrix[:, 2].argsort()]
        job_matrix = numpy.flipud(job_matrix)

    return job_matrix

    def testCircular(self):
        g = Graph()

        op = OpLazyCC(graph=g)
        op.ChunkShape.setValue((3, 3, 1))

        vol = np.asarray(
            [
                [0, 0, 0, 0, 0, 0, 0, 0, 0],
                [0, 1, 1, 1, 0, 1, 1, 1, 0],
                [0, 1, 0, 0, 0, 0, 0, 1, 0],
                [0, 1, 0, 0, 0, 0, 0, 1, 0],
                [0, 1, 0, 0, 0, 0, 0, 1, 0],
                [0, 1, 0, 0, 0, 0, 0, 1, 0],
                [0, 1, 0, 0, 0, 0, 0, 1, 0],
                [0, 1, 1, 1, 1, 1, 1, 1, 0],
                [0, 0, 0, 0, 0, 0, 0, 0, 0],
            ],
            dtype=np.uint8,
        )
        vol1 = vigra.taggedView(vol, axistags="yx")
        vol2 = vigra.taggedView(vol, axistags="xy")
        vol3 = vigra.taggedView(np.flipud(vol), axistags="yx")
        vol4 = vigra.taggedView(np.flipud(vol), axistags="xy")

        for v in (vol1, vol2, vol3, vol4):
            op.Input.setValue(v)
            for x in [0, 3, 6]:
                for y in [0, 3, 6]:
                    if x == 3 and y == 3:
                        continue
                    op.Input.setDirty(slice(None))
                    out = op.Output[x : x + 3, y : y + 3].wait()
                    print(out.squeeze())
                    assert out.max() == 1

def _interpolate_for_irf(w_orig, w_interp, mat_in):
  '''
  Interpolate matrices for the IRF calculations
  '''
  mat_interp = np.zeros([mat_in.shape[0], mat_in.shape[1], w_interp.size])

  flip = False

  if w_orig[0] > w_orig[1]:

    w_tmp = np.flipud(w_orig)
    flip = True

  else:

    w_tmp = w_orig

  for i in xrange(mat_in.shape[0]):

    for j in xrange(mat_in.shape[1]):

      if flip is True:

        rdTmp = np.flipud(mat_in[i, j, :])

      else:
        rdTmp = mat_in[i, j, :]

      f = interpolate.interp1d(x=w_tmp, y=rdTmp)
      mat_interp[i, j, :] = f(w_interp)

  return mat_interp

def var_border(v, di=1, dj=1):
    """
  Border of 2d numpy array
  di,dj is the interval between points along columns and lines
  Corner points are kept even with di and dj not 1
  """

    j, i = v.shape
    if (di, dj) == (1, 1):
        xb = np.arange(2 * i + 2 * j, dtype=v.dtype)
        yb = np.arange(2 * i + 2 * j, dtype=v.dtype)

        xb[0:j] = v[:, 0]
        xb[j : j + i] = v[-1, :]
        xb[j + i : j + i + j] = np.flipud(v[:, -1])
        xb[j + i + j :] = np.flipud(v[0, :])
    else:
        # ensure corner points are kept!!
        tmp1 = v[::dj, 0]
        tmp2 = v[-1, ::di]
        tmp3 = np.flipud(v[:, -1])[::dj]
        tmp4 = np.flipud(v[0, :])[::di]
        xb = np.concatenate((tmp1, tmp2, tmp3, tmp4))

    return xb

	def _northup( self, latitude='lat' ):
		''' this works only for global grids to be downscaled flips it northup '''
		if self.ds[ latitude ][0].data < 0: # meaning that south is north globally
			self.ds[ latitude ] = np.flipud( self.ds[ latitude ] )
			# flip each slice of the array and make a new one
			flipped = np.array( [ np.flipud( arr ) for arr in self.ds[ self.variable ].data ] )
			self.ds[ self.variable ] = (('time', 'lat', 'lon' ), flipped )

    def zwriter(self, imageslice, fname):
        if "OpenEXR" in self.backends:
            imageslice = numpy.flipud(imageslice)
            exr.save_depth(imageslice, fname)

        elif "VTK" in self.backends:
            height = imageslice.shape[1]
            width = imageslice.shape[0]

            file = open(fname, mode='w')
            file.write("Image type: L 32F image\r\n")
            file.write("Name: A cinema depth image\r\n")
            file.write("Image size (x*y): "+str(height) + "*" + str(width) + "\r\n")
            file.write("File size (no of images): 1\r\n")
            file.write(chr(26))
            imageslice.tofile(file)
            file.close()

        elif "PIL" in self.backends:
            imageslice = numpy.flipud(imageslice)
            pimg = PIL.Image.fromarray(imageslice)
            #TODO:
            # don't let ImImagePlugin.py insert the Name: filename in line two
            # why? because ImImagePlugin.py reader has a 100 character limit
            pimg.save(fname)

        else:
            print "Warning: need OpenEXR or PIL or VTK to write to " + fname

def polysmooth(x,y,z,NI,NJ):

    # size of the incoming array
    Nx, Ny = np.shape(z)
    x1d = x[:,0]
    y1d = y[0,:]

    # Get the C coefficients
    #NI = 7
    CIj = np.zeros((NI,Ny))
    for j in range (Ny):
        CIj[:,j] = np.flipud(np.polyfit(x1d,z[:,j],NI-1))

    # Get the D coefficients
    #NJ = 7
    DIJ = np.zeros((NI,NJ))
    for I in range (NI):
        DIJ[I,:] = np.flipud(np.polyfit(y1d,CIj[I,:],NJ-1))
    
    # Reconstruct the entire surface
    zsmooth = np.zeros((Nx,Ny))
    for I in range(NI):
        for J in range(NJ):
            zsmooth += DIJ[I,J]*x**I*y**J

    return zsmooth

def filtfilt(b,a,x):
    """
    What does this function do?  Alberto??
    """
    
    #For now only accepting 1d arrays
    ntaps=max(len(a),len(b))
    edge=ntaps*3
        
    if x.ndim != 1:
        raise ValueError, "Filiflit is only accepting 1 dimension arrays."

    #x must be bigger than edge
    if x.size < edge:
        raise ValueError, "Input vector needs to be bigger than 3 * max(len(a),len(b)."
        
    if len(a) < ntaps:
	a=np.r_[a,zeros(len(b)-len(a))]

    if len(b) < ntaps:
	b=np.r_[b,zeros(len(a)-len(b))]         
    
    zi=lfilter_zi(b,a)
    
    #Grow the signal to have edges for stabilizing 
    #the filter with inverted replicas of the signal
    s=np.r_[2*x[0]-x[edge:1:-1],x,2*x[-1]-x[-1:-edge:-1]]
    #in the case of one go we only need one of the extrems 
    # both are needed for filtfilt
    
    (y,zf)=lfilter(b,a,s,-1,zi*s[0])

    (y,zf)=lfilter(b,a,np.flipud(y),-1,zi*y[-1])
    
    return np.flipud(y[edge-1:-edge+1])

def apply_Gumbel_original(gumbelFile, trgFolder, prefix, return_periods, cellArea, logger):
    # read gumbel file
    nc_src = nc4.Dataset(gumbelFile, 'r')
    # read axes and revert the y-axis
    x = nc_src.variables['lon'][:]
    y = np.flipud(nc_src.variables['lat'][:])
    
    # read different variables
    loc   = nc_src.variables['flvol_location'][0,:,:]#; loc    = gumbel_loc.data;  loc[loc==gumbel_loc._FillValue]       = np.nan
    scale = nc_src.variables['flvol_scale'][0,:,:]#;    scale  = gumbel_scale.data;scale[scale==gumbel_scale._FillValue] = np.nan
    p_zero    = nc_src.variables['flvol_zero_prob'][0,:,:]#;p_zero = zero_prob.data;   p_zero[p_zero==zero_prob._FillValue]  = np.nan
    # loop over all return periods
    for return_period in return_periods:
        logger.info('Preparing return period %05.f' % return_period)
        flvol = inv_gumbel(p_zero, loc, scale, return_period)
        # any area with cell > 0, fill in a zero. This may occur because:
            # a) dynRout produces missing values (occuring in some pixels in the Sahara)
            # b) the forcing data is not exactly overlapping the cell area mask (e.g. EU-WATCH land cells are slightly different from PCR-GLOBWB mask)
            # c) the probability of zero flooding is 100%. This causes a division by zero in the inv_gumbel function
        test = logical_and(flvol.mask, cellArea > 0)
        flvol[test] = 0.
        test = logical_and(np.isnan(flvol), cellArea > 0)
        flvol[test] = 0.
        # if any values become negative due to the statistical extrapolation, fix them to zero (may occur if the sample size for fitting was small and a small return period is requested)
        flvol = np.maximum(flvol, 0.)
        # write to a PCRaster file
        flvol_data = flvol.data
        # finally mask the real not-a-number cells
        flvol_data[flvol.mask] = -9999.
        fileName = os.path.join(trgFolder, '%s_RP_%05.f.map') % (prefix, return_period)
        writeMap(fileName, 'PCRaster', x, y, np.flipud(flvol_data), -9999.)
    nc_src.close()

    def MatForChannel(signals, picHeight, picWidth):
        # signals is a vector! not a mat

        newSignals = signals * ((picHeight - 1) / 2)
        width = len(newSignals)
        stride = int(np.ceil(float(width) / picWidth))

        # pad
        temp = np.zeros(picWidth * stride)
        temp[0:width] = newSignals

        # 1 - reshape according to stride
        maxVec = temp.reshape((picWidth, stride))
        # 2 - get max vector of each stride
        maxVec = maxVec.max(axis=1).astype(int)
        # 3 - get max mat
        maxMat = np.zeros((picWidth, picHeight / 2))
        maxMat[np.arange(picWidth), maxVec] = 1

        # flip it when it's little...
        # TODO: maybe we can cheat  this part (but it's really fast)
        maxMat = maxMat.T
        maxMat = np.flipud(maxMat)

        # 4 - create fill matrix
        # TODO: can do a static mat of this
        fillMat = np.ones((picHeight / 2, picHeight / 2))
        fillMat = np.tril(fillMat, 0)

        # result
        endMat = np.dot(fillMat, maxMat)
        endMat = np.vstack((endMat, np.flipud(endMat)))

        return endMat