Python matplotlib.pyplot 模块，Normalize() 实例源码

def draw_heatmaps(screen, surf, hm, thr=0.5, vmin=-15, vmax=10):
    hm_idx = [
            ( 8, 0*hmsurf_size[0], 0*hmsurf_size[1]),   # R. wrist
            ( 9, 1*hmsurf_size[0], 0*hmsurf_size[1]),   # L. wrist
            ( 6, 0*hmsurf_size[0], 1*hmsurf_size[1]),   # R. elbow
            ( 7, 1*hmsurf_size[0], 1*hmsurf_size[1]),   # L. elbow
            ( 3, 0*hmsurf_size[0], 2*hmsurf_size[1]),   # Head
            ( 0, 1*hmsurf_size[0], 2*hmsurf_size[1]),   # Pelvis
            (12, 0*hmsurf_size[0], 3*hmsurf_size[1]),   # R. knee
            (13, 1*hmsurf_size[0], 3*hmsurf_size[1])]   # L. knee

    for idx in hm_idx:
        h = np.transpose(hm[:,:,idx[0]].copy(), (1, 0))
        h[h < vmin] = vmin
        h[h > vmax] = vmax
        cmap = plt.cm.jet
        norm = plt.Normalize(vmin=vmin, vmax=vmax)
        cm = np.zeros((34, 34, 3))
        cm[1:33, 1:33, :] = cmap(norm(h))[:,:,0:3]
        cm = scipy.ndimage.zoom(cm, (5, 5, 1), order=1)
        pygame.surfarray.pixels3d(surf)[:,:,:] = np.array(255.*cm, dtype=int)
        screen.blit(surf, (idx[1] + img_size[0], idx[2]))

def colorline(x, y, z=None, cmap=plt.get_cmap('copper'),
              norm=plt.Normalize(0.0, 1.0), linewidth=3, alpha=1.0):
    """
    Plot a colored line with coordinates x and y
    Optionally specify colors in the array z
    Optionally specify a colormap, a norm function and a line width
    """

    # Default colors equally spaced on [0,1]:
    if z is None:
        z = np.linspace(0.0, 1.0, len(x))

    z = np.asarray(z)

    segments = make_segments(x, y)
    lc = LineCollection(segments, array=z, cmap=cmap, norm=norm, \
                        linewidth=linewidth, alpha=alpha)
    ax = plt.gca()
    ax.add_collection(lc)
    return lc

def colorline(x, y, z=None, cmap=plt.get_cmap('Spectral_r'), cmin=None, cmax=None, lw=3):
    '''
    Plot a colored line with coordinates x and y
    Optionally specify colors in the array z
    Optionally specify a colormap, a norm function and a line width
    '''

    cmap = copy(cmap)
    cmap.set_over('k')
    cmap.set_under('k')

    # Default colors equally spaced on [0,1]:
    if z is None:
        z = np.linspace(0.0, 1.0, len(x))

    # Special case if a single number:
    if not hasattr(z, "__iter__"):  # to check for numerical input -- this is a hack
        z = np.array([z])

    z = np.asarray(z)

    segments = make_segments(x, y)
    return LineCollection(segments, array=z, cmap=cmap, norm=plt.Normalize(vmin=cmin, vmax=cmax), linewidth=lw)

def animate(self, j):
        '''
        Play frame `j` of the animation.

        '''

        if not self.pause:

            # Normalize the time index
            j = int(j * len(self.t) / 100.)

            # Time tracker
            self.tracker.set_xdata(self.bodies[0].time_hr[self.t[j]])

            # Occultor images
            x0 = self.body.x_hr[self.t[j]]
            y0 = self.body.y_hr[self.t[j]]
            for k, occultor in enumerate(self.occultors):
                if self.xy is None:
                    xo, yo = occultor.x_hr[self.t[j]] - x0, \
                             occultor.y_hr[self.t[j]] - y0
                else:
                    xo, yo = self.xy(occultor.x_hr[self.t[j]] - x0,
                                     occultor.y_hr[self.t[j]] - y0)
                self.occ[k].center = (xo / self.body._r, yo / self.body._r)

                # BODY orbits
                for k, b in enumerate(self.bodies):
                    self.pt_xz[k].set_xdata(b.x_hr[self.t[j]])
                    self.pt_xz[k].set_ydata(b.z_hr[self.t[j]])
                    self.pt_zy[k].set_xdata(b.z_hr[self.t[j]])
                    self.pt_zy[k].set_ydata(b.y_hr[self.t[j]])

def create_color_bar(f, cmap, colorbar_axis, bar_font, epochsInds, title):
    sm = plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=0, vmax=1))
    sm._A = []
    cbar_ax = f.add_axes(colorbar_axis)
    cbar = f.colorbar(sm, ticks=[], cax=cbar_ax)
    cbar.ax.tick_params(labelsize=bar_font)
    cbar.set_label(title, size=bar_font)
    cbar.ax.text(0.5, -0.01, epochsInds[0], transform=cbar.ax.transAxes,
                 va='top', ha='center', size=bar_font)
    cbar.ax.text(0.5, 1.0, str(epochsInds[-1]), transform=cbar.ax.transAxes,
                 va='bottom', ha='center', size=bar_font)

def cm_bounds_to_norm(cm_bounds, data=None):
    cmin = float(cm_bounds[0]) if cm_bounds[0] is not None \
               else (np.min(data) if data is not None else 0)
    cmax = float(cm_bounds[1]) if cm_bounds[1] is not None \
               else (np.max(data) if data is not None else cmin+1)
    return matplotlibpyplotNormalize(cmin, cmax)

def colorbar(cmap="jet", cm_bounds=(0, 1)):
    """
    cmap, norm, mappable = colorbar('jet', min, max)
    plt.scatter(x, y, c, cmap=cmap, norm=norm)
    cb = plt.colorbar(mappable)
    if arr is given, forces cm_min and cm_max to min-max of the arr
    """
    if isinstance(cmap, str):
        cmap = mat.cm.get_cmap(cmap)
    cmin, cmax = list(map(float, cm_bounds[:2]))
    norm = matplotlibpyplotNormalize(cmin, cmax)
    mappable = mat.cm.ScalarMappable(cmap=cmap, norm=norm)
    mappable._A = []
    return cmap, norm, mappable

def compute_node_colors(self):
        """Compute the node colors. Also computes the colorbar."""
        data = [self.graph.node[n][self.node_color] for n in self.nodes]
        data_reduced = sorted(list(set(data)))
        dtype = infer_data_type(data)
        n_grps = num_discrete_groups(data)

        if dtype == 'categorical' or dtype == 'ordinal':
            cmap = get_cmap(cmaps['Accent_{0}'.format(n_grps)].mpl_colormap)
        elif dtype == 'continuous' and not is_data_diverging(data):
            cmap = get_cmap(cmaps['continuous'].mpl_colormap)
        elif dtype == 'continuous' and is_data_diverging(data):
            cmap = get_cmap(cmaps['diverging'].mpl_colormap)

        for d in data:
            idx = data_reduced.index(d) / n_grps
            self.node_colors.append(cmap(idx))

        # Add colorbar if required.
        logging.debug('length of data_reduced: {0}'.format(len(data_reduced)))
        logging.debug('dtype: {0}'.format(dtype))
        if len(data_reduced) > 1 and dtype == 'continuous':
            self.sm = plt.cm.ScalarMappable(cmap=cmap,
                                            norm=plt.Normalize(vmin=min(data_reduced),  # noqa
                                                               vmax=max(data_reduced)   # noqa
                                                               )
                                            )
            self.sm._A = []

def colorline(x, y, cmap=None, cm_range=(0, 0.7), **kwargs):
    """Colorline plots a trajectory of (x,y) points with a colormap"""

    # plt.plot(x, y, '-k', zorder=1)
    # plt.scatter(x, y, s=40, c=plt.cm.RdBu(np.linspace(0,1,40)), zorder=2, edgecolor='k')

    assert len(cm_range)==2, "cm_range must have (min, max)"
    assert len(x) == len(y), "x and y must have the same number of elements!"

    ax = kwargs.get('ax', plt.gca())
    lw = kwargs.get('lw', 2)
    if cmap is None:
        cmap=plt.cm.Blues_r

    t = np.linspace(cm_range[0], cm_range[1], len(x))

    points = np.array([x, y]).T.reshape(-1, 1, 2)
    segments = np.concatenate([points[:-1], points[1:]], axis=1)

    lc = LineCollection(segments, cmap=cmap, norm=plt.Normalize(0, 1),
                        zorder=50)
    lc.set_array(t)
    lc.set_linewidth(lw)

    ax.add_collection(lc)

    return lc

def plotmap(mapname, axlines, segments, coords, figtype="png"):
    """plotmap(mapname, axlines, segments, coords)
    plot axial map with integration data as colour values
    """
    import matplotlib.pyplot as plt
    import matplotlib.cm as cm

    fig,axes=plt.subplots()

    valary=[]
    for ax in axlines:
        I = axlines[ax]["Integration"]
        valary.append(I)

    cmap = cm.gnuplot
    norm=plt.Normalize(vmin=max(0, min(valary)), vmax=max(valary))
    sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
    sm.set_array(valary)

    for i,a in enumerate(axlines):
        cval=float(axlines[a]["Integration"])
        for seg in axlines[a]["segments"]:
            n1, n2 = segments[seg]
            x1, y1 = coords[n1]
            x2, y2 = coords[n2]
            axes.plot([x1,x2], [y1,y2], color=sm.to_rgba(cval), linewidth=3)
    axes.set_xlabel("Longitude")
    axes.set_ylabel("Latitude")
    cb=plt.colorbar(sm, ax=axes)
    cb.set_label("Integration")

    plt.savefig(mapname+"-integration."+figtype)

def plot_fancy(nodes, elems, phi, charge, u=None, charge_max=None,
               show=False, save=None):
    """ Plots fancily. """
    fig = Figure(colorbar=False, tight_layout=True, show=show,
                 xlabel="", ylabel="", save=save, ticks=False)

    if charge_max is None:
        charge_max = max(np.max(np.abs(charge)), 1e-10)

    cmap = plt.cm.get_cmap('Greys')
    cmap._init()
    cmap._lut[:, :] = 0.
    length = len(cmap._lut[:, -1])
    # cmap._lut[:, -1] = np.linspace(0., 1.0, length)
    cmap._lut[:length/2, -1] = 0.
    cmap._lut[length/2:, -1] = 1.

    phi[phi > 1.] = 1.
    phi[phi < -1.] = -1.

    plt.tripcolor(nodes[:, 0], nodes[:, 1], elems, charge,
                  cmap=plt.get_cmap("coolwarm"), shading="gouraud",
                  vmin=-charge_max, vmax=charge_max)
    plt.tricontourf(nodes[:, 0], nodes[:, 1], elems, phi,
                    cmap=cmap, levels=[-2.0, 0., 2.0], antialiased=True)

    if u is not None:
        u_norm = np.sqrt(u[:, 0]**2 + u[:, 1]**2) + 1e-10
        colors = phi
        norm = plt.Normalize()
        norm.autoscale(colors)
        colormap = cmap  # plt.cm.get_cmap('inferno')
        cmap._lut[:, -1] = 0.5
        cmap._lut[length/2:, :-1] = 1.

        fig.ax.quiver(nodes[:, 0], nodes[:, 1],
                      u[:, 0]/u_norm, u[:, 1]/u_norm,
                      color=colormap(norm(colors)))

    return fig

def imshow2d(data, ax=None, cmap2d='brightwheel', huenorm=None, huevmin=None,
             huevmax=None, lightnorm=None, lightvmin=None, lightvmax=None,
             **kwargs):
    """
    Plot 2 parameter 2D data array to current axis.

    :param data: numpy array with shape (2, nwidth, nheight). The first index
                 corresponds to the hue and the second to the lightness of the
                 colors.
    :param ax: a matplotlib axis instance.
    :param cmap: either:
                 numpy array with shape (nwidth, nheight, 4) that contains
                 the 4 rgba values in hue (width) and lightness (height).
                 Can be obtained by a call to get_cmap2d(name).
                 or:
                 name where name is one of the following strings:
                 'brightwheel', 'darkwheel', 'hardwheel', 'newwheel',
                 'smoothwheel', 'wheel'
    :param huenorm: a plt.Normalize() instance that normalizes the hue values.
    :param huevmin: the minimum of the huevalues. Only used if huenorm=None.
    :param huevmax: the maximum of the huevalues. Only used if huenorm=None.
    :param lightnorm: a plt.Normalize() instance that normalizes the lightness
                      values.
    :param lightvmin: the minimum of the lightness values.
                      Only used if lightnorm=None.
    :param lightvmax: the maximum of the lightness values.
                      Only used if lightnorm=None.
    :param **kwargs: remaining kwargs are passed to plt.imshow()
    """
    if ax is None:
        ax = plt.gca()
    rgb_data = data2d_to_rgb(data, cmap2d=cmap2d,
                             huenorm=huenorm, huevmin=huevmin,
                             huevmax=huevmax, lightnorm=lightnorm,
                             lightvmin=lightvmin, lightvmax=lightvmax)
    im = ax.imshow(rgb_data, **kwargs)
    return im

def colorbar(cmap="jet", cm_min=0, cm_max=1):
    if isinstance(cmap, str):
        cmap = cmget_cmap(cmap)
    norm = matplotlibpyplotNormalize(cm_min, cm_max)
    mappable = cmScalarMappable(cmap=cmap, norm=norm)
    mappable._A = []
    return cmap, norm, mappable

def normalize_colors(vmin, vmax, clip=False):
    """Helper to handle matplotlib API"""
    import matplotlib.pyplot as plt
    try:
        return plt.Normalize(vmin, vmax, clip=clip)
    except AttributeError:
        return plt.normalize(vmin, vmax, clip=clip)

def _colorline(ax, x, y, color = (0, 0, 0), **kwargs):
    '''
    Plots the curve `y(x)` with linearly increasing alpha.
    Adapted from `http://nbviewer.jupyter.org/github/dpsanders/matplotlib-examples/blob/master/colorline.ipynb`_.

    '''

    # A bit hacky... But there doesn't seem to be
    # an easy way to get the hex code for a named color...
    if isinstance(color, string_types):
        if color.startswith("#"):
            hex = color[1:]
        else:
            if len(color) == 1:
                if color == 'k':
                    color = 'black'
                elif color == 'r':
                    color = 'red'
                elif color == 'b':
                    color = 'blue'
                elif color == 'g':
                    color = 'green'
                elif color == 'y':
                    color = 'yellow'
                elif color == 'w':
                    color = 'white'
                else:
                    # ?!
                    color = 'black'
            hex = matplotlib.colors.cnames[color.lower()][1:]
        r, g, b = tuple(int(hex[i:i+2], 16) / 255. for i in (0, 2, 4))
    else:
        r, g, b = color

    colors = [(r, g, b, i) for i in np.linspace(0, 1, 3)]
    cmap = LinearSegmentedColormap.from_list('alphacmap', colors, N = 1000)
    points = np.array([x, y]).T.reshape(-1, 1, 2)
    segments = np.concatenate([points[:-1], points[1:]], axis=1)
    lc = LineCollection(segments, array = np.linspace(0.0, 1.0, len(x)),
                        cmap = cmap, norm = pl.Normalize(0.0, 1.0), **kwargs)
    ax.add_collection(lc)

    return lc

def save_pareto_fitness_plot():
    """
    Saves a plot of the current fitness for a pareto front.

    :return: Nothing
    """

    from algorithm.parameters import params

    # Initialise up figure instance.
    fig = plt.figure()
    ax1 = fig.add_subplot(1, 1, 1)

    # Set up iterator for color plotting.
    color = iter(plt.cm.jet(np.linspace(0, 1, len(first_pareto_list))))

    # Get labels for individual fitnesses.
    ffs = params['FITNESS_FUNCTION'].fitness_functions

    # Find the direction for step lines to "bend"
    step_dir = 'pre' if ffs[0].maximise else 'post'

    # Plot data.
    for i, gen in enumerate(first_pareto_list):
        c = next(color)
        ax1.step(gen[0], gen[1], linestyle='--',
                 where=step_dir, color=c, lw=0.35, alpha=0.25)
        ax1.plot(gen[0], gen[1], 'o', color=c, ms=1)

    # Set labels with class names.
    ax1.set_xlabel(ffs[0].__class__.__name__, fontsize=14)
    ax1.set_ylabel(ffs[1].__class__.__name__, fontsize=14)

    # Plot title and legend.
    plt.title("First pareto fronts by generation")

    # Set up colorbar instead of legend. Normalise axis to scale of data.
    sm = plt.cm.ScalarMappable(cmap="jet",
                   norm=plt.Normalize(vmin=0, vmax=len(first_pareto_list) - 1))

    # Fake up the array of the scalar mappable.
    sm._A = []

    # Plot the colorbar.
    cbar = plt.colorbar(sm, ticks=[0, len(first_pareto_list) - 1])

    # Set label of colorbar.
    # cbar.ax.get_yaxis().labelpad = 15
    cbar.ax.set_ylabel('Generation', rotation=90)

    # Save plot and close.
    plt.savefig(path.join(params['FILE_PATH'], "fitness.pdf"))
    plt.close()

def main():
    # some parameters
    npts_real = 200
    npts_imag = 100
    realaxis = np.linspace(-2 * np.pi, 2 * np.pi, npts_real)
    imagaxis = np.linspace(-1., 1., npts_imag)
    regrid, imgrid = np.meshgrid(realaxis, imagaxis)

    # put complex magnitude and argument in 2d array
    complex_sine = np.sin(regrid + 1j * imgrid)
    data = np.empty((2, npts_imag, npts_real))
    data[0] = np.abs(complex_sine)
    data[1] = np.angle(complex_sine)

    # normalize array
    norm0 = plt.Normalize(0, data[0].max())
    norm1 = plt.Normalize(-np.pi, np.pi)
    data[0] = norm0(data[0])
    data[1] = norm1(data[1])

    paths_cmap = glob.glob('colormap2d/colormaps/*.npy')
    ncmaps = len(paths_cmap) + 1  # one extra cmap for hsv
    fig, axes = plt.subplots(ncmaps, 2, figsize=(10, ncmaps * 3))
    for path_cmap, (col1, col2) in zip(paths_cmap, axes):
        dirname, fname = os.path.split(path_cmap)
        cmap = np.load(path_cmap).transpose((1, 0, 2))
        col1.set(title='{}'.format(fname))
        col2.set(title='complex sine')

        rgb_colors = cmap_file2d(data, cmap)

        col1.imshow(cmap, aspect='auto')
        col2.imshow(rgb_colors, aspect='auto')

    # hsv colormap
    ax20, ax21 = axes[-1, 0], axes[-1, 1]
    ax20.set(title='HSV colormap')
    ax21.set(title='complex sine')

    xx, yy = np.meshgrid(np.linspace(0., 1., 100), np.linspace(0., 1., 100))
    cmap_grid = np.array([yy, xx])
    cmap = cmap_multidim_hsv(cmap_grid)
    cmap = np.roll(cmap, 48, axis=1)
    rgb_colors = cmap_file2d(data, cmap)
    ax20.imshow(cmap, aspect='auto')
    ax21.imshow(rgb_colors, aspect='auto')

    fig.tight_layout(pad=0.5)
    plt.show()

def scatterColor(x0, y, w):
    """Creates scatter plot with points colored by variable.
    All input arrays must have matching lengths

    :param x0: x values to plot
    :type x0: list
    :param y: y values to plot
    :type y: list
    :param w: z values to plot

    :returns: plot; slope and intercept of the RLM best fit line shown on the plot
    .. warning:: all input arrays must have matching lengths and scalar values
    .. note:: See documentation at http://statsmodels.sourceforge.net/0.6.0/generated/statsmodels.robust.robust_linear_model.RLM.html
    for the RLM line
    """
    import matplotlib as mpl
    import matplotlib.cm as cm
    import statsmodels.api as sm
    from scipy.stats import linregress
    cmap = plt.cm.get_cmap('RdYlBu')
    norm = mpl.colors.Normalize(vmin=w.min(), vmax=w.max())
    m = cm.ScalarMappable(norm=norm, cmap=cmap)
    m.set_array(w)
    sc = plt.scatter(x0, y, label='', color=m.to_rgba(w))

    xa = sm.add_constant(x0)

    est = sm.RLM(y, xa).fit()
    r2 = sm.WLS(y, xa, weights=est.weights).fit().rsquared
    slope = est.params[1]

    x_prime = np.linspace(np.min(x0), np.max(x0), 100)[:, np.newaxis]
    x_prime = sm.add_constant(x_prime)
    y_hat = est.predict(x_prime)

    const = est.params[0]
    y2 = [i * slope + const for i in x0]

    lin = linregress(x0, y)
    x1 = np.arange(np.min(x0), np.max(x0), 0.1)
    y1 = [i * lin[0] + lin[1] for i in x1]
    y2 = [i * slope + const for i in x1]
    plt.plot(x1, y1, c='g',
             label='simple linear regression m = {:.2f} b = {:.0f}, r^2 = {:.2f}'.format(lin[0], lin[1], lin[2] ** 2))
    plt.plot(x1, y2, c='r', label='rlm regression m = {:.2f} b = {:.0f}, r2 = {:.2f}'.format(slope, const, r2))
    plt.legend()
    cbar = plt.colorbar(m)

    cbar.set_label('Julian Date')

    return slope, const

def main():
    args = parse_args()
    _dataframe = DF()
    for spec in args.specid:
        if args.dark:
            if op.exists(op.join(args.outfolder, 'bias_DF_%s.csv' %spec)):
                _bias_dataframe = pd.read_csv(op.join(args.outfolder, 
                                          'bias_DF_%s.csv' %spec))
            else:
                print("Cannot find the dataframe for the biases.")

        ifuslot = CAM_IFUSLOT_DICT[spec]
        lower_folder_struct = op.join('virus', 'virus*', 'exp*', 'virus',
                                      '2*_%s*zro.fits' %ifuslot)
        progress = ProgressBar(len(args.cal_dirs), spec, fmt=ProgressBar.FULL)
        for date in args.cal_dirs:
            files = glob.glob(op.join(date,lower_folder_struct))
            for fn in files:
                with warnings.catch_warnings():
                    warnings.simplefilter("ignore")
                    _dataframe = build_dataframe(_dataframe, op.basename(date), fn, spec)
            progress.current+=1
            progress()
        progress.done()
    norm = plt.Normalize()
    colors = plt.cm.viridis_r(norm(np.arange(9+2)))
    if args.zero:
        _dataframe.to_csv('bias_DF_%s.csv' %spec)
    for amp in AMPS:
        fig1 = plt.figure(figsize=(8,6))
        fig2 = plt.figure(figsize=(8,6))
        fig3 = plt.figure(figsize=(8,6))
        df = _dataframe.query('AMP=="%s"'%amp)
        for i in xrange(9):
            strv = 'VAL' + str(i)
            df = df[(is_outlier(df['overscan'])<1)*(is_outlier(df[strv])<1)*
                    (is_outlier(df['temp'])<1)]
            plt.figure(fig1.number)
            plt.scatter(df['overscan'],df[strv]-df['overscan'], edgecolor='none',
                        s=25, color=colors[i,0:3], alpha=0.3)
            plt.xlabel('Overscan [ADU]')
            plt.ylabel('BIAS [ADU]')
            plt.figure(fig2.number)
            plt.scatter(df['temp'],df[strv]-df['overscan'], edgecolor='none',
                        s=25, color=colors[i,0:3], alpha=0.3)
            plt.xlabel('TEMP')
            plt.ylabel('BIAS [ADU]')
            plt.figure(fig3.number)
            plt.scatter(df['mjd'],df[strv]-df['overscan'], edgecolor='none',
                        s=25, color=colors[i,0:3], alpha=0.3)
            plt.xlabel('MJD')
            plt.ylabel('BIAS [ADU]')
        plt.figure(fig1.number)
        plt.savefig(op.join(args.output,'bias_struct_%s_%s_overscan.pdf' %(spec, amp)),dpi=150)
        plt.figure(fig2.number)
        plt.savefig(op.join(args.output,'bias_struct_%s_%s_temp.pdf' %(spec, amp)),dpi=150)
        plt.figure(fig3.number)
        plt.savefig(op.join(args.output,'bias_struct_%s_%s_mjd.pdf' %(spec, amp)),dpi=150)
        plt.close(fig1)
        plt.close(fig2)
        plt.close(fig3)

def throughput_fiberextract(Felist, args):
    nifu = len(Felist)
    nw = len(Felist[0][0].data[0,:])
    xp = np.linspace(0, 1, num=nw)
    nbspline = 12
    a = np.linspace(0, 1, nbspline)
    knots = np.hstack([0,0,np.vstack([a,a]).T.ravel(),1,1])
    b = Bspline(knots, 3)
    basis = np.array([b(xi) for xi in xp])
    B = np.zeros((nifu,nw))
    for i in xrange(nifu):
        spec = biweight_location(Felist[i][0].data,axis=(0,))
        mask = np.where((~np.isnan(spec))*(~np.isinf(spec))*(spec!=0))[0]
        sol = np.linalg.lstsq(basis[mask,:], spec[mask])[0]
        B[i,:] = np.dot(basis,sol)
#        if args.plot:  
#            pltfile = op.join(args.outfolder, 'spectrum_%i.pdf' %i)
#            fig = plt.figure(figsize=(8, 6))
#            plt.plot(xp, spec)
#            plt.plot(xp, B[i,:],'r--')
#            plt.xticks([])
#            plt.xlabel('Wavelength')
#            plt.ylabel('Arbitrary Units')
#            plt.xlim([0, 1])
#            fig.savefig(pltfile, dpi=150)
#            plt.close()
    if args.plot:
        norm = plt.Normalize()
        colors = plt.cm.viridis(norm(np.arange(nifu)+1))
        pltfile = op.join(args.outfolder, 'IFU_average_spectra.pdf')
        fig = plt.figure(figsize=(8, 6))
        avgB = biweight_location(B,axis=(0,))
        for i in xrange(nifu):
            with np.errstate(divide='ignore'):
                plt.plot(xp, B[i,:]/avgB, color=colors[i,0:3], alpha=0.9)
        plt.xticks([])
        plt.xlabel('Wavelength')
        plt.ylabel('Normalized Units')
        plt.xlim([0, 1])
        plt.ylim([0.5,1.5])
        fig.savefig(pltfile, dpi=150)
        plt.close()
    return B, avgB