def __call__(self, value, clip=None):
     if clip is None:
         clip = self.clip
     if cbook.iterable(value):
         vtype = 'array'
         val = ma.asarray(value).astype(np.float)
     else:
         vtype = 'scalar'
         val = ma.array([value]).astype(np.float)
     self.autoscale_None(val)
     vmin, vmax = self.vmin, self.vmax
     cmin, cmax = self.cmin * vmin, self.cmax * vmax
     if vmin > vmax:
         raise ValueError("minvalue must be less than or equal to maxvalue")
     elif vmin == vmax:
         result = 0.0 * val
     else:
         if clip:
             mask = ma.getmask(val)
             val = ma.array(np.clip(val.filled(vmax), vmin, vmax),
                             mask=mask)
         result = 0. * val + 0.5
         result[val > cmax] = (ma.log10(val[val > cmax]) - ma.log10(cmax)) / (np.log10(vmax) - np.log10(cmax)) / 2. + 0.5
         result[val < cmin] = -(ma.log10(-val[val < cmin]) - ma.log10(-cmin)) / (np.log10(-vmin) - np.log10(-cmin)) / 2. + 0.5
     if vtype == 'scalar':
         result = result[0]
     return result

 def test_testUfuncs1(self):
     # Test various functions such as sin, cos.
     (x, y, a10, m1, m2, xm, ym, z, zm, xf, s) = self.d
     assert_(eq(np.cos(x), cos(xm)))
     assert_(eq(np.cosh(x), cosh(xm)))
     assert_(eq(np.sin(x), sin(xm)))
     assert_(eq(np.sinh(x), sinh(xm)))
     assert_(eq(np.tan(x), tan(xm)))
     assert_(eq(np.tanh(x), tanh(xm)))
     with np.errstate(divide='ignore', invalid='ignore'):
         assert_(eq(np.sqrt(abs(x)), sqrt(xm)))
         assert_(eq(np.log(abs(x)), log(xm)))
         assert_(eq(np.log10(abs(x)), log10(xm)))
     assert_(eq(np.exp(x), exp(xm)))
     assert_(eq(np.arcsin(z), arcsin(zm)))
     assert_(eq(np.arccos(z), arccos(zm)))
     assert_(eq(np.arctan(z), arctan(zm)))
     assert_(eq(np.arctan2(x, y), arctan2(xm, ym)))
     assert_(eq(np.absolute(x), absolute(xm)))
     assert_(eq(np.equal(x, y), equal(xm, ym)))
     assert_(eq(np.not_equal(x, y), not_equal(xm, ym)))
     assert_(eq(np.less(x, y), less(xm, ym)))
     assert_(eq(np.greater(x, y), greater(xm, ym)))
     assert_(eq(np.less_equal(x, y), less_equal(xm, ym)))
     assert_(eq(np.greater_equal(x, y), greater_equal(xm, ym)))
     assert_(eq(np.conjugate(x), conjugate(xm)))
     assert_(eq(np.concatenate((x, y)), concatenate((xm, ym))))
     assert_(eq(np.concatenate((x, y)), concatenate((x, y))))
     assert_(eq(np.concatenate((x, y)), concatenate((xm, y))))
     assert_(eq(np.concatenate((x, y, x)), concatenate((x, ym, x))))

def average_in_flux(mag, dmag, axis=None):
    flux = 10**(mag / -2.5)
    dflux = np.log(10) / 2.5 * flux * dmag
    avg_dflux = np.power(np.sum(np.power(dflux, -2), axis), -0.5)
    avg_flux = np.sum(flux * np.power(dflux, -2), axis) * avg_dflux**2
    avg_mag = -2.5 * np.log10(avg_flux)
    avg_dmag = 2.5 / np.log(10) * np.divide(avg_dflux, avg_flux)
    return avg_mag, avg_dmag

    def inverse(self, value):

        # ORIGINAL MATPLOTLIB CODE

        if not self.scaled():
            raise ValueError("Not invertible until scaled")

        vmin, vmax = self.vmin, self.vmax

        # CUSTOM APLPY CODE

        if cbook.iterable(value):
            val = ma.asarray(value)
        else:
            val = value

        if self.stretch == 'linear':

            pass

        elif self.stretch == 'log':

            val = (ma.power(10., val * ma.log10(self.midpoint)) - 1.) / (self.midpoint - 1.)

        elif self.stretch == 'sqrt':

            val = val * val

        elif self.stretch == 'arcsinh':

            val = self.midpoint * \
                  ma.sinh(val * ma.arcsinh(1. / self.midpoint))

        elif self.stretch == 'square':

            val = ma.power(val, (1. / 2))

        elif self.stretch == 'power':

            val = ma.power(val, (1. / self.exponent))

        else:

            raise Exception("Unknown stretch in APLpyNormalize: %s" %
                            self.stretch)

        return vmin + val * (vmax - vmin)

 def transform_non_affine(self, a):
     """
     This transform takes an Nx1 ``numpy`` array and returns a
     transformed copy.  Since the range of the Mercator scale
     is limited by the user-specified threshold, the input
     array must be masked to contain only valid values.
     ``matplotlib`` will handle masked arrays and remove the
     out-of-range data from the plot.  Importantly, the
     ``transform`` method *must* return an array that is the
     same shape as the input array, since these values need to
     remain synchronized with values in the other dimension.
     """
     masked = ma.masked_where(((a - self.zero) * self.sign <= 0.0), (a - self.zero) * self.sign)
     if masked.mask.any():
         return ma.log10(masked) / np.log10 (self.base)
     else:
         return np.log10((a - self.zero) * self.sign) / np.log10 (self.base)

 def transform_non_affine(self, a):
     lower = a[np.where(a<=change)]
     greater = a[np.where(a> change)]
     if lower.size:
         lower = self._handle_nonpos(lower * 10.0)/10.0
         if isinstance(lower, ma.MaskedArray):
             lower = ma.log10(lower)
         else:
             lower = np.log10(lower)
         lower = factor*lower
     if greater.size:
         greater = (factor*np.log10(change) + (greater-change))
     # Only low
     if not(greater.size):
         return lower
     # Only high
     if not(lower.size):
         return greater
     return np.concatenate((lower, greater))

    def __call__(self, value, clip=None):

        #read in parameters
        method = self.stretch
        exponent = self.exponent
        midpoint = self.midpoint

        # ORIGINAL MATPLOTLIB CODE

        if clip is None:
            clip = self.clip

        if cbook.iterable(value):
            vtype = 'array'
            val = ma.asarray(value).astype(np.float)
        else:
            vtype = 'scalar'
            val = ma.array([value]).astype(np.float)

        self.autoscale_None(val)
        vmin, vmax = self.vmin, self.vmax
        if vmin > vmax:
            raise ValueError("minvalue must be less than or equal to maxvalue")
        elif vmin==vmax:
            return 0.0 * val
        else:
            if clip:
                mask = ma.getmask(val)
                val = ma.array(np.clip(val.filled(vmax), vmin, vmax),
                                mask=mask)
            result = (val-vmin) * (1.0/(vmax-vmin))

            # CUSTOM APLPY CODE

            # Keep track of negative values
            negative = result < 0.

            if self.stretch == 'linear':

                pass

            elif self.stretch == 'log':

                result = ma.log10(result * (self.midpoint - 1.) + 1.) \
                       / ma.log10(self.midpoint)

            elif self.stretch == 'sqrt':

                result = ma.sqrt(result)

            elif self.stretch == 'arcsinh':

                result = ma.arcsinh(result/self.midpoint) \
                       / ma.arcsinh(1./self.midpoint)

            elif self.stretch == 'power':

                result = ma.power(result, exponent)

            else:

                raise Exception("Unknown stretch in APLpyNormalize: %s" %
                                self.stretch)

            # Now set previously negative values to 0, as these are
            # different from true NaN values in the FITS image
            result[negative] = -np.inf

        if vtype == 'scalar':
            result = result[0]

        return result

 def transform_non_affine(self, a):
     a = self._handle_nonpos(a * 10.0)
     if isinstance(a, ma.MaskedArray):
         return ma.log10(a)
     return np.log10(a)

def find_max(array_like):
    no_zeros = array_like[nonzero(array_like)]
    with errstate(all='ignore'):
        max = nanmax(10 ** log10(no_zeros))
    return max

where_are_NaNs = isnan(img_masked1)
img_masked1[where_are_NaNs] = 0

plt.imshow(img_masked1, origin = 'lower',aspect='auto') 
plt.pcolor(img_masked1,norm=LogNorm()) 
plt.set_cmap('seismic')
cbar=plt.colorbar() 


with np.errstate(divide='ignore', invalid='ignore'):
    ratio = np.true_divide(img_masked,img_masked1)
    ratio[ratio == np.inf] = 0
    ratio= np.nan_to_num(ratio)

 ratio_norm=ratio/2.76
  balmer = 2.5*(ma.log10(ratio_norm))

plt.imshow(balmer, origin = 'lower',aspect='auto') 
plt.pcolor(balmer,norm=LogNorm()) 
plt.set_cmap('jet')
cbar=plt.colorbar() 


ratio_norm=ratio/2.76
balmer = 2.5*(ma.log(ratio_norm))
print balmer.filled(0)
outfile = 'balmer.fits.gz'

hdu = fits.PrimaryHDU(balmer)
hdu.writeto(outfile, clobber=True)

def main():
    # parse command-line arguments
    parser = argparse.ArgumentParser(
        formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    ## targets to fit
    parser.add_argument("--name", type=str, default=None,
        help="base name of combined skim file")
    parser.add_argument("--subsample-step", type=int, default=1000,
        help="step size used for subsampling observations")
    parser.add_argument("--dont-save", action="store_true",
        help="dont save delta field (just do all the preprocessing)")
    args = parser.parse_args()

    # import data
    skim = h5py.File(args.name+'.hdf5', 'r')
    norm = skim['norm'][:][:,np.newaxis]
    loglam = skim['loglam'][:]
    wave = np.power(10.0, loglam)

    quasar_redshifts = skim['z'][:]

    linear_continuum = h5py.File(args.name+'-linear-continuum.hdf5', 'r')
    params_a = linear_continuum['params_a'].value
    params_b = linear_continuum['params_b'].value
    continuum = linear_continuum['continuum'].value
    continuum_wave = linear_continuum['continuum_wave'].value
    continuum_interp = scipy.interpolate.UnivariateSpline(continuum_wave, continuum, s=0, ext=1)

    wave_lya = linear_continuum.attrs['wave_lya']
    abs_alpha = linear_continuum.attrs['abs_alpha']
    abs_beta = linear_continuum.attrs['abs_beta']
    forest_wave_ref = linear_continuum.attrs['forest_wave_ref']

    print 'Adjusting weights for pipeline variance and LSS variance...'

    forest_min_z = linear_continuum.attrs['forest_min_z']
    forest_max_z = linear_continuum.attrs['forest_max_z']
    forest_dz = 0.1
    forest_z_bins = np.arange(forest_min_z, forest_max_z + forest_dz, forest_dz)
    var_lss = scipy.interpolate.UnivariateSpline(forest_z_bins, 0.05 + 0.06*(forest_z_bins - 2.0)**2, s=0)
    var_pipe_scale = scipy.interpolate.UnivariateSpline(forest_z_bins, 0.7 + 0.2*(forest_z_bins - 2.0)**2, s=0)

    forest_pixel_redshifts = wave/wave_lya - 1
    abs_coefs = abs_alpha*np.power(1+forest_pixel_redshifts, abs_beta)

    forest_wave_refs = forest_wave_ref*(1+quasar_redshifts)
    def model_flux(a, b):
        return a*np.power(wave/forest_wave_refs[:,np.newaxis], b)*continuum_interp(wave/(1+quasar_redshifts[:,np.newaxis]))*np.exp(-abs_coefs)

    mflux = model_flux(params_a[:,np.newaxis],params_b[:,np.newaxis])

    # (1.0 + quasar_redshifts[:,np.newaxis])*forest_wave/args.wave_lya - 1.0

    print forest_pixel_redshifts.shape

    pixel_mask = skim['mask'][:]

    print pixel_mask.shape

    flux = np.ma.MaskedArray(skim['flux'][:], mask=pixel_mask)
    ivar = np.ma.MaskedArray(skim['ivar'][:], mask=pixel_mask)

    delta_flux = flux/mflux - 1.0
    delta_ivar = ivar*mflux*mflux

    delta_weight = delta_ivar*var_pipe_scale(forest_pixel_redshifts)
    delta_weight = delta_weight/(1 + delta_weight*var_lss(forest_pixel_redshifts))

    redshift_order = np.argsort(quasar_redshifts)
    export_exact_image(args.name+'-delta-flux.png', delta_flux[redshift_order][::args.subsample_step], dpi=100,
        vmin=-5, vmax=5, cmap=plt.get_cmap('bwr'), origin='lower')
    export_exact_image(args.name+'-delta-weight.png', ma.log10(delta_flux[redshift_order][::args.subsample_step]), dpi=100,
        vmin=-5, vmax=2, cmap=plt.get_cmap('Purples'), origin='lower')
    export_exact_image(args.name+'-delta-mask.png', pixel_mask[redshift_order][::args.subsample_step], dpi=100,
        origin='lower')

    print 'Computing mean delta...'

    mask_params = (params_a > .1) & (params_a < 10) & (params_b > -10) & (params_b < 10)

    delta_mean = ma.average(delta_flux[mask_params], axis=0)
    delta_mean_weighted = ma.average(delta_flux[mask_params], weights=delta_weight[mask_params], axis=0)
    delta_mean_ivar_weighted = ma.average(delta_flux[mask_params], weights=delta_ivar[mask_params], axis=0)

    plt.figure(figsize=(12,9))
    plt.plot(wave, delta_mean, label='Unweighted Mean')
    plt.plot(wave, delta_mean_weighted, label='LSS weighted Mean')
    plt.plot(wave, delta_mean_ivar_weighted, label='Ivar weighted Mean')
    # plt.ylim(0.06*np.array([-1,1]))
    plt.xlabel(r'Observed Wavelength ($\AA$)')
    plt.ylabel(r'Delta Mean')
    plt.grid()
    plt.legend()
    plt.savefig(args.name+'-lssweighted-delta-mean.png', dpi=100, bbox_inches='tight')
    plt.close()

    if args.dont_save:
        return -1

    outfile = h5py.File(args.name+'-delta.hdf5', 'w')
#.........这里部分代码省略.........

 def transform_non_affine(self, a):
   masked = ma.masked_where(a > 1-10**(-1-self.nines), a)
   if masked.mask.any():
     return -ma.log10(1-a)
   else:
     return -np.log10(1-a)