Reads and returns the energy bin edges, fluxes and npreds from
    a FITs HDU
179920
def read_spectral_data(hdu):
    
    ebins = read_energy_bounds(hdu)
    fluxes = np.ndarray((len(ebins)))
    try:
        fluxes[0:-1] = hdu.data.field("E_MIN_FL")
        fluxes[-1] = hdu.data.field("E_MAX_FL")[-1]
        npreds = hdu.data.field("NPRED")
    except:
        fluxes = np.ones((len(ebins)))
        npreds = np.ones((len(ebins)))
    return ebins, fluxes, npreds
995945
def create_from_flux_points(cls, txtfile):
        """"""

        tab = Table.read(txtfile, format='ascii.ecsv')
        dnde_unit = u.ph / (u.MeV * u.cm ** 2 * u.s)
        loge = np.log10(np.array(tab['e_ref'].to(u.MeV)))
        norm = np.array(tab['norm'].to(dnde_unit))
        norm_errp = np.array(tab['norm_errp'].to(dnde_unit))
        norm_errn = np.array(tab['norm_errn'].to(dnde_unit))
        norm_err = 0.5 * (norm_errp + norm_errn)
        dloge = loge[1:] - loge[:-1]
        dloge = np.insert(dloge, 0, dloge[0])
        emin = 10 ** (loge - dloge * 0.5)
        emax = 10 ** (loge + dloge * 0.5)
        ectr = 10 ** loge
        deltae = emax - emin
        flux = norm * deltae
        eflux = norm * deltae * ectr

        ones = np.ones(flux.shape)
        ref_spec = ReferenceSpec(emin, emax, ones, ones, ones, ones)

        spec_data = SpecData(ref_spec, norm, flux, eflux, norm_err)

        stephi = np.linspace(0, 1, 11)
        steplo = -np.linspace(0, 1, 11)[1:][::-1]

        loscale = 3 * norm_err
        hiscale = 3 * norm_err
        loscale[loscale > norm] = norm[loscale > norm]

        norm_vals_hi = norm[:, np.newaxis] + \
            stephi[np.newaxis, :] * hiscale[:, np.newaxis]
        norm_vals_lo = norm[:, np.newaxis] + \
            steplo[np.newaxis, :] * loscale[:, np.newaxis]

        norm_vals = np.hstack((norm_vals_lo, norm_vals_hi))
        nll_vals = 0.5 * \
            (norm_vals - norm[:, np.newaxis]) ** 2 / \
            norm_err[:, np.newaxis] ** 2

        norm_vals *= flux[:, np.newaxis] / norm[:, np.newaxis]

        return cls(norm_vals, nll_vals, spec_data, 'flux')
1127378
def get_shape(bin_edges, sid):
    
    mag_bins, dist_bins, lon_bins, lat_bins, eps_bins = bin_edges
    return (len(mag_bins) - 1, len(dist_bins) - 1,
            len(lon_bins[sid]) - 1, len(lat_bins[sid]) - 1, len(eps_bins) - 1)
339008
def compute_ck2004_ldints(self):
        

        if 'ck2004_all' not in self.content:
            print('Castelli & Kurucz (2004) intensities are not computed yet. Please compute those first.')
            return None

        ldaxes = self._ck2004_intensity_axes
        ldtable = self._ck2004_Imu_energy_grid
        pldtable = self._ck2004_Imu_photon_grid

        self._ck2004_ldint_energy_grid = np.nan*np.ones((len(ldaxes[0]), len(ldaxes[1]), len(ldaxes[2]), 1))
        self._ck2004_ldint_photon_grid = np.nan*np.ones((len(ldaxes[0]), len(ldaxes[1]), len(ldaxes[2]), 1))

        mu = ldaxes[3]
        Imu = 10**ldtable[:,:,:,:]/10**ldtable[:,:,:,-1:]
        pImu = 10**pldtable[:,:,:,:]/10**pldtable[:,:,:,-1:]

        # To compute the fluxes, we need to evaluate \int_0^1 2pi Imu mu dmu.

        for a in range(len(ldaxes[0])):
            for b in range(len(ldaxes[1])):
                for c in range(len(ldaxes[2])):

                    ldint = 0.0
                    pldint = 0.0
                    for i in range(len(mu)-1):
                        ki = (Imu[a,b,c,i+1]-Imu[a,b,c,i])/(mu[i+1]-mu[i])
                        ni = Imu[a,b,c,i]-ki*mu[i]
                        ldint += ki/3*(mu[i+1]**3-mu[i]**3) + ni/2*(mu[i+1]**2-mu[i]**2)

                        pki = (pImu[a,b,c,i+1]-pImu[a,b,c,i])/(mu[i+1]-mu[i])
                        pni = pImu[a,b,c,i]-pki*mu[i]
                        pldint += pki/3*(mu[i+1]**3-mu[i]**3) + pni/2*(mu[i+1]**2-mu[i]**2)

                    self._ck2004_ldint_energy_grid[a,b,c] = 2*ldint
                    self._ck2004_ldint_photon_grid[a,b,c] = 2*pldint

        self.content.append('ck2004_ldint')
850972
def _load_fits(self, h5file):
        
        fits = {}
        for key in ['mf']:
            fits[key] = self._load_scalar_fit(fit_key=key, h5file=h5file)
        for key in ['chif', 'vf']:
            fits[key] = self._load_vector_fit(key, h5file)
        return fits
1090141
def get_kernel_data(self):
        """"""
        return {
            'scratch_mot_float_type': LocalMemory(
                'mot_float_type', 8 +
                                  2 * self._var_replace_dict['NMR_OBSERVATIONS'] +
                                  5 * self._var_replace_dict['NMR_PARAMS'] +
                                  self._var_replace_dict['NMR_PARAMS'] * self._var_replace_dict['NMR_OBSERVATIONS']),
            'scratch_int': LocalMemory('int', self._var_replace_dict['NMR_PARAMS'])
        }
1127151
def get_wulff_shape(self, material_id):
        
        from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
        from pymatgen.analysis.wulff import WulffShape, hkl_tuple_to_str

        structure = self.get_structure_by_material_id(material_id)
        surfaces = self.get_surface_data(material_id)["surfaces"]
        lattice = (SpacegroupAnalyzer(structure)
                   .get_conventional_standard_structure().lattice)
        miller_energy_map = {}
        for surf in surfaces:
            miller = tuple(surf["miller_index"])
            # Prefer reconstructed surfaces, which have lower surface energies.
            if (miller not in miller_energy_map) or surf["is_reconstructed"]:
                miller_energy_map[miller] = surf["surface_energy"]
        millers, energies = zip(*miller_energy_map.items())
        return WulffShape(lattice, millers, energies)
541739
def get_total_flux(self, energies):
        

        fluxes = []

        for src in self._point_sources:

            fluxes.append(self._point_sources[src](energies))

        return np.sum(fluxes, axis=0)
230201
def extract_images_from_tscube(infile, outfile):
    
    inhdulist = fits.open(infile)
    wcs = pywcs.WCS(inhdulist[0].header)
    map_shape = inhdulist[0].data.shape

    t_eng = Table.read(infile, "EBOUNDS")
    t_scan = Table.read(infile, "SCANDATA")
    t_fit = Table.read(infile, "FITDATA")

    n_ebin = len(t_eng)
    energies = np.ndarray((n_ebin + 1))
    energies[0:-1] = t_eng["E_MIN"]
    energies[-1] = t_eng["E_MAX"][-1]

    cube_shape = (n_ebin, map_shape[1], map_shape[0])

    wcs_cube = wcs_utils.wcs_add_energy_axis(wcs, energies)

    outhdulist = [inhdulist[0], inhdulist["EBOUNDS"]]

    FIT_COLNAMES = ['FIT_TS', 'FIT_STATUS', 'FIT_NORM',
                    'FIT_NORM_ERR', 'FIT_NORM_ERRP', 'FIT_NORM_ERRN']
    SCAN_COLNAMES = ['TS', 'BIN_STATUS', 'NORM', 'NORM_UL',
                     'NORM_ERR', 'NORM_ERRP', 'NORM_ERRN', 'LOGLIKE']

    for c in FIT_COLNAMES:
        data = t_fit[c].data.reshape(map_shape)
        hdu = fits.ImageHDU(data, wcs.to_header(), name=c)
        outhdulist.append(hdu)

    for c in SCAN_COLNAMES:
        data = t_scan[c].data.swapaxes(0, 1).reshape(cube_shape)
        hdu = fits.ImageHDU(data, wcs_cube.to_header(), name=c)
        outhdulist.append(hdu)

    hdulist = fits.HDUList(outhdulist)
    hdulist.writeto(outfile, clobber=True)
    return hdulist
1233406
def read(self, line, f, data):
        """"""
        data["energy"] = float(f.readline().split()[1])
        N = len(data["symbols"])
        # if the data are already read before, just overwrite them
        gradient = data.get("gradient")
        if gradient is None:
            gradient = np.zeros((N,3), float)
            data["gradient"] = gradient
        for i in range(N):
            words = f.readline().split()
            gradient[i,0] = float(words[2])
            gradient[i,1] = float(words[3])
            gradient[i,2] = float(words[4])
1193044
def read_projection_from_fits(fitsfile, extname=None):
    
    f = fits.open(fitsfile)
    nhdu = len(f)
    # Try and get the energy bounds
    try:
        ebins = find_and_read_ebins(f)
    except:
        ebins = None

    if extname is None:
        # If there is an image in the Primary HDU we can return a WCS-based
        # projection
        if f[0].header['NAXIS'] != 0:
            proj = WCS(f[0].header)
            return proj, f, f[0]
    else:
        if f[extname].header['XTENSION'] == 'IMAGE':
            proj = WCS(f[extname].header)
            return proj, f, f[extname]
        elif extname in ['SKYMAP', 'SKYMAP2']:
            proj = HPX.create_from_hdu(f[extname], ebins)
            return proj, f, f[extname]
        elif f[extname].header['XTENSION'] == 'BINTABLE':
            try:
                if f[extname].header['PIXTYPE'] == 'HEALPIX':
                    proj = HPX.create_from_hdu(f[extname], ebins)
                    return proj, f, f[extname]
            except:
                pass
        return None, f, None

    # Loop on HDU and look for either an image or a table with HEALPix data
    for i in range(1, nhdu):
        # if there is an image we can return a WCS-based projection
        if f[i].header['XTENSION'] == 'IMAGE':
            proj = WCS(f[i].header)
            return proj, f, f[i]
        elif f[i].header['XTENSION'] == 'BINTABLE':
            if f[i].name in ['SKYMAP', 'SKYMAP2']:
                proj = HPX.create_from_hdu(f[i], ebins)
                return proj, f, f[i]
            try:
                if f[i].header['PIXTYPE'] == 'HEALPIX':
                    proj = HPX.create_from_hdu(f[i], ebins)
                    return proj, f, f[i]
            except:
                pass

    return None, f, None
61604
def energy_data():
    
    cur = db.cursor().execute("""SELECT timestamp, current FROM Readings""")
    original = TimeSeries()
    original.initialize_from_sql_cursor(cur)
    original.normalize("day", fusionMethod = "sum")
    return itty.Response(json.dumps(original, cls=PycastEncoder), content_type='application/json')
1305820
def energy_profile(self,ixaxis):
        

        mass = self.get('mass')
        radius = self.get('radius') * ast.rsun_cm
        eps_nuc = self.get('eps_nuc')
        eps_neu = self.get('non_nuc_neu')

        if ixaxis == 'mass':
            xaxis = mass
            xlab = 'Mass / M$_\odot$'
        else:
            xaxis = old_div(radius, 1.e8) # Mm
            xlab = 'Radius (Mm)'

        pl.plot(xaxis, np.log10(eps_nuc),
                'k-',
                label='$\epsilon_\mathrm{nuc}>0$')
        pl.plot(xaxis, np.log10(-eps_nuc),
                'k--',
                label='$\epsilon_\mathrm{nuc}<0$')
        pl.plot(xaxis, np.log10(eps_neu),
                'r-',
                label='$\epsilon_\\nu$')

        pl.xlabel(xlab)
        pl.ylabel('$\log(\epsilon_\mathrm{nuc},\epsilon_\\nu)$')
        pl.legend(loc='best').draw_frame(False)
1125515
def fit(self, X):
        
        self.centers_, self.labels_, self.sse_arr_, self.n_iter_ = \
              _kmedoids(X, self.n_clusters, self.distance, self.max_iter, self.n_trials, self.tol, self.rng)
69830
def read_excitation_energies(self):
        

        transitions = list()

        # read in file
        with zopen(self.filename, "r") as f:
            line = f.readline()
            td = False
            while line != "":
                if re.search(r"^\sExcitation energies and oscillator strengths:", line):
                    td = True

                if td:
                    if re.search(r"^\sExcited State\s*\d", line):
                        val = [float(v) for v in float_patt.findall(line)]
                        transitions.append(tuple(val[0:3]))
                line = f.readline()
        return transitions
1278048
def _read_or_calc_samples(
    sampler,
    modelidx=0,
    n_samples=100,
    last_step=False,
    e_range=None,
    e_npoints=100,
    threads=None,
):
    

    if not e_range:
        # return the results saved in blobs
        modelx, model = _process_blob(sampler, modelidx, last_step=last_step)
    else:
        # prepare bogus data for calculation
        e_range = validate_array(
            "e_range", u.Quantity(e_range), physical_type="energy"
        )
        e_unit = e_range.unit
        energy = (
            np.logspace(
                np.log10(e_range[0].value),
                np.log10(e_range[1].value),
                e_npoints,
            )
            * e_unit
        )
        data = {
            "energy": energy,
            "flux": np.zeros(energy.shape) * sampler.data["flux"].unit,
        }
        # init pool and select parameters
        chain = sampler.chain[-1] if last_step else sampler.flatchain
        pars = chain[np.random.randint(len(chain), size=n_samples)]
        blobs = []

        p = Pool(threads)
        modelouts = p.map(partial(sampler.modelfn, data=data), pars)
        p.close()
        p.terminate()

        for modelout in modelouts:
            if isinstance(modelout, np.ndarray):
                blobs.append([modelout])
            else:
                blobs.append(modelout)

        modelx, model = _process_blob(
            blobs, modelidx=modelidx, energy=data["energy"]
        )

    return modelx, model
1206852
def _generate_energy_edges_single(ene):
    """"""
    midene = np.sqrt((ene[1:] * ene[:-1]))
    elo, ehi = np.zeros(len(ene)) * ene.unit, np.zeros(len(ene)) * ene.unit
    elo[1:] = ene[1:] - midene
    ehi[:-1] = midene - ene[:-1]
    elo[0] = ene[0] * (1 - ene[0] / (ene[0] + ehi[0]))
    ehi[-1] = elo[-1]
    return u.Quantity([elo, ehi])