Python QubicAcquisition.get_noise方法代码示例

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_noise [as 别名]
def get_qubic_map(instrument, sampling, scene, input_maps, withplanck=True, covlim=0.1):
    acq = QubicAcquisition(instrument, sampling, scene, photon_noise=True, effective_duration=1)
    C = acq.get_convolution_peak_operator()
    coverage = acq.get_coverage()
    observed = coverage > covlim * np.max(coverage)
    acq_restricted = acq[:, :, observed]
    H = acq_restricted.get_operator()
    x0_convolved = C(input_maps)
    if not withplanck:
        pack = PackOperator(observed, broadcast='rightward')
        y_noiseless = H(pack(x0_convolved))
        noise = acq.get_noise()
        y = y_noiseless + noise
        invntt = acq.get_invntt_operator()
        A = H.T * invntt * H
        b = (H.T * invntt)(y)
        preconditioner = DiagonalOperator(1 / coverage[observed], broadcast='rightward')
        solution_qubic = pcg(A, b, M=preconditioner, disp=True, tol=1e-3, maxiter=1000)
        maps = pack.T(solution_qubic['x'])
        maps[~observed] = 0
    else:
        acq_planck = PlanckAcquisition(150, acq.scene, true_sky=x0_convolved)#, fix_seed=True)
        acq_fusion = QubicPlanckAcquisition(acq, acq_planck)
        map_planck_obs=acq_planck.get_observation()
        H = acq_fusion.get_operator()
        invntt = acq_fusion.get_invntt_operator()
        y = acq_fusion.get_observation()
        A = H.T * invntt * H
        b = H.T * invntt * y
        solution_fusion = pcg(A, b, disp=True, maxiter=1000, tol=1e-3)
        maps = solution_fusion['x']
        maps[~observed] = 0
    x0_convolved[~observed] = 0
    return(maps, x0_convolved, observed)    

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_noise [as 别名]
def make_a_map(x0, pointing, instrument, nside, coverage_threshold=0.01, todnoise=None, fits_string=None, noiseless=False):
    ############# Make TODs ###########################################################
    acquisition = QubicAcquisition(instrument, pointing,
                                 nside=nside,
                                 synthbeam_fraction=0.99)
    tod, x0_convolved = map2tod(acquisition, x0, convolution=True)
    if todnoise is None:
        todnoise = acquisition.get_noise()
    factnoise=1
    if noiseless:
        factnoise=0
    ##################################################################################

    
    ############# Make mapss ###########################################################
    print('Making map')
    maps, cov = tod2map_all(acquisition, tod + todnoise * factnoise, tol=1e-4, coverage_threshold=coverage_threshold)
    if MPI.COMM_WORLD.rank == 0:
        fitsmapname = 'maps_'+fits_string+'.fits'
        fitscovname = 'cov_'+fits_string+'.fits'
        print('Saving the map: '+fitsmapname)
        qubic.io.write_map(fitsmapname,maps)
        print('Saving the coverage: '+fitsmapname)
        qubic.io.write_map(fitscovname,cov)
    ##################################################################################

    return maps, cov, todnoise

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_noise [as 别名]
def get_maps(spectra, inst, sampling, nside, x0, coverage_threshold=0.01,savefile=None, savefile_noiseless=None, noI=False):
  #if x0 is None:
  #  print("Running Synfast")
  #  x0 = np.array(hp.synfast(spectra[1:],nside,fwhm=0,pixwin=True,new=True)).T
  #  if noI: x0[:,0]*=0
  acquisition = QubicAcquisition(inst, sampling,
                                nside=nside,
                                synthbeam_fraction=0.99)
                                #max_nbytes=6e9)
  # simulate the timeline
  print('Now doing MAP2TOD (simulate data)')
  tod, x0_convolved = map2tod(acquisition, x0, convolution=True)
  print('TOD Done, now adding noise')
  bla = acquisition.get_noise()
  tod_noisy = tod + bla 
  # reconstruct using all available bolometers
  print('Now doing TOD2MAP (make map)')
  map_all, cov_all = tod2map_all(acquisition, tod_noisy, tol=1e-4, coverage_threshold=coverage_threshold)
  print('Map done')
  print('Now doing TOD2MAP (make map) on noiseless data')
  map_all_noiseless, cov_all_noiseless = tod2map_all(acquisition, tod, tol=1e-4, coverage_threshold=coverage_threshold)
  print('Noiseless Map done')
  mask = map_all_noiseless != 0
  x0_convolved[~mask,:] = 0
  rank = MPI.COMM_WORLD.rank
  if rank == 0:
    if savefile is not None:
      print('I am rank='+str(rank)+' and I try to save the file '+savefile)
      FitsArray(np.array([map_all[:,0], map_all[:,1], map_all[:,2], cov_all, x0_convolved[:,0], x0_convolved[:,1], x0_convolved[:,2]]), copy=False).save(savefile)
      print('I am rank='+str(rank)+' and I just saved the file '+savefile)
      print('I am rank='+str(rank)+' and I try to save the file '+savefile_noiseless)
      FitsArray(np.array([map_all_noiseless[:,0], map_all_noiseless[:,1], map_all_noiseless[:,2], cov_all_noiseless, x0_convolved[:,0], x0_convolved[:,1], x0_convolved[:,2]]), copy=False).save(savefile_noiseless)
      print('I am rank='+str(rank)+' and I just saved the file '+savefile_noiseless)
  return map_all, x0_convolved, cov_all, x0

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_noise [as 别名]
def get_tod(instrument, sampling, scene, input_maps, withplanck=True, covlim=0.1, photon_noise=True): 
    acq = QubicAcquisition(instrument, sampling, scene, photon_noise=photon_noise)
    C = acq.get_convolution_peak_operator()
    coverage = acq.get_coverage()
    observed = coverage > covlim * np.max(coverage)
    acq_restricted = acq[:, :, observed]
    H = acq_restricted.get_operator()
    x0_convolved = C(input_maps)
    pack = PackOperator(observed, broadcast='rightward')
    y_noiseless = H(pack(x0_convolved))
    noise = acq.get_noise()
    y = y_noiseless + noise
    return (y_noiseless, noise, y)

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_noise [as 别名]
def TOD(subnu_min,subnu_max,subdelta_nu,cmb,dust,sampling,scene,effective_duration,verbose=True, photon_noise=True):

	sh = cmb.shape
	Nbpixels = sh[0]

	###################
	### Frequencies ###
	###################

	Nbsubfreq=int(floor(log(subnu_max/subnu_min)/log(1+subdelta_nu)))+1
	sub_nus_edge=subnu_min*np.logspace(0,log(subnu_max/subnu_min)/log(1+subdelta_nu),Nbsubfreq,endpoint=True,base=subdelta_nu+1)
	sub_nus=np.array([(sub_nus_edge[i]+sub_nus_edge[i-1])/2 for i in range(1,Nbsubfreq)])
	sub_deltas=np.array([(sub_nus_edge[i]-sub_nus_edge[i-1]) for i in range(1,Nbsubfreq)])
	Delta=subnu_max-subnu_min
	Nbsubbands=len(sub_nus) ## = Nbsubfreq-1

	if verbose:
		print('Nombre de bandes utilisées pour la construction : '+str(Nbsubbands))
		print('Sous fréquences centrales utilisées pour la construction : '+str(sub_nus))
		print('Edges : '+str(sub_nus_edge))

	################
	### Coverage ###
	################
	sub_instruments=[QubicInstrument(filter_nu=sub_nus[i]*10**9,filter_relative_bandwidth=sub_deltas[i]/sub_nus[i],detector_nep=2.7e-17) for i in range(Nbsubbands)]
	sub_acqs=[QubicAcquisition(sub_instruments[i], sampling, scene,photon_noise=photon_noise, effective_duration=effective_duration) for i in range(Nbsubbands)]
	covlim=0.1
	coverage = np.array([sub_acqs[i].get_coverage() for i in range(Nbsubbands)])
	observed = [(coverage[i] > covlim*np.max(coverage[i])) for i in range(Nbsubbands)]
	obs=reduce(logical_and,tuple(observed[i] for i in range(Nbsubbands)))
	pack = PackOperator(obs, broadcast='rightward')

	#################
	### Input map ###
	#################

	x0=np.zeros((Nbsubbands,Nbpixels,3))
	for i in range(Nbsubbands):
		#x0[i,:,0]=cmb.T[0]+dust.T[0]*scaling_dust(150,sub_nus[i]e-9,1.59)
		x0[i,:,1]=cmb.T[1]+dust.T[1]*scaling_dust(150,sub_nus[i],1.59)
		x0[i,:,2]=cmb.T[2]+dust.T[2]*scaling_dust(150,sub_nus[i],1.59)


	###############################
	### Construction of the TOD ###
	###############################
	dnu=sub_instruments[0].filter.bandwidth
	Y=0
	Y_noisy=0


	# Noiseless TOD
	for i in range(Nbsubbands):
		sub_acqs_restricted=sub_acqs[i][:,:,obs]
		operator=sub_acqs_restricted.get_operator()
		C=HealpixConvolutionGaussianOperator(fwhm=sub_instruments[i].synthbeam.peak150.fwhm * (150 / (sub_nus[i])))
		Y=Y+operator*pack*C*x0[i]


	# Global instrument creqted to get the noise over the entire instrument bandwidth
	Global_instrument=QubicInstrument(filter_nu=(subnu_max+subnu_min)/2,filter_relative_bandwidth=Delta/((subnu_max+subnu_min)/2),detector_nep=2.7e-17)
	Global_acquisition=QubicAcquisition(Global_instrument, sampling, scene,photon_noise=photon_noise, effective_duration=effective_duration)


	noise_instrument=Global_acquisition.get_noise()
	#sigma=np.std(noise_instrument)
	#mean=np.mean(noise_instrument)
	#white_noise=np.random.normal(mean,sigma,shape(Y))
	#Y_noisy=Y+white_noise

	#noise=sub_acqs[0].get_noise()*np.sum(sub_deltas)*10**9/dnu
	Y_noisy=Y + noise_instrument

	return Y_noisy,obs

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_noise [as 别名]
# acquisition model
acq = QubicAcquisition(150, sampling, kind='I', synthbeam_fraction=0.99,
                       detector_sigma=sigma, detector_fknee=fknee,
                       detector_fslope=fslope, detector_ncorr=ncorr)
C = acq.get_convolution_peak_operator()
P = acq.get_projection_operator()
H = P * C

# produce the Time-Ordered data
y = H(x0)

# noise
psd = _gaussian_psd_1f(len(acq.sampling), sigma=sigma, fknee=fknee,
                       fslope=fslope, sampling_frequency=1/ts)
invntt = acq.get_invntt_operator()
noise = acq.get_noise()

# map-making
coverage = P.pT1()
mask = coverage > 10
P = P.restrict(mask, inplace=True)
unpack = UnpackOperator(mask)

# map without covariance matrix
solution1 = pcg(P.T * P, P.T(y + noise),
                M=DiagonalOperator(1/coverage[mask]), disp=True)
x1 = unpack(solution1['x'])

# map with covariance matrix
solution2 = pcg(P.T * invntt * P, (P.T * invntt)(y + noise),
                M=DiagonalOperator(1/coverage[mask]), disp=True)

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_noise [as 别名]
  print('new I map RMS is : '+str(np.std(x0_noI[:,0])))
  print('new Q map RMS is : '+str(np.std(x0_noI[:,1])))
  print('new U map RMS is : '+str(np.std(x0_noI[:,2])))
  ##################################################################################





  ############# Make TODs ###########################################################
  acquisition = QubicAcquisition(instFull, new_sampling,
                                 nside=nside,
                                 synthbeam_fraction=0.99)
  tod, x0_convolved = map2tod(acquisition, x0, convolution=True)
  tod_noI, x0_noI_convolved = map2tod(acquisition, x0_noI, convolution=True)
  todnoise = acquisition.get_noise()
  ##################################################################################




  ############# Make maps ###########################################################
  th_acquisition = QubicAcquisition(instFull, sampling,
                                 	nside=nside,
                                 	synthbeam_fraction=0.99)

  strrnd = random_string(10)

  fits_noI_input = 'maps_ns'+str(nside)+'_noise'+str(noise)+'_sigptg'+str(sigptg)+'_noI_input_'+strrnd+'.fits'
  fits_input = 'maps_ns'+str(nside)+'_noise'+str(noise)+'_sigptg'+str(sigptg)+'_input_'+strrnd+'.fits'

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_noise [as 别名]

#.........这里部分代码省略.........
    Nbbands = len(nus)

    # frequencies assumed to have been used for construction of TOD
    subnu_min = nu_min
    subnu_max = nu_max
    Nbsubfreq = int(floor(log(subnu_max / subnu_min) / log(1 + subdelta_nu))) + 1
    sub_nus_edge = subnu_min * np.logspace(
        0, log(subnu_max / subnu_min) / log(1 + subdelta_nu), Nbsubfreq, endpoint=True, base=subdelta_nu + 1
    )
    sub_nus = np.array([(sub_nus_edge[i] + sub_nus_edge[i - 1]) / 2 for i in range(1, Nbsubfreq)])
    sub_deltas = np.array([(sub_nus_edge[i] - sub_nus_edge[i - 1]) for i in range(1, Nbsubfreq)])
    Nbsubbands = len(sub_nus)

    # Bands
    bands = [sub_nus[reduce(logical_and, (sub_nus <= nus_edge[i + 1], sub_nus >= nus_edge[i]))] for i in range(Nbbands)]
    numbers = np.cumsum(np.array([len(bands[i]) for i in range(Nbbands)]))
    numbers = np.append(0, numbers)
    bands_numbers = np.array([(np.arange(numbers[i], numbers[i + 1])) for i in range(Nbbands)])

    if verbose:
        print ("Nombre de bandes utilisées pour la reconstruction : " + str(Nbsubbands))
        print ("Nombre de bandes reconstruites : " + str(Nbbands))
        print ("Résolution spectrale : " + str(delta_nu))
        print ("Bandes reconstruites : " + str(bands))

    ################
    ### Coverage ###
    ################
    sub_instruments = [
        QubicInstrument(
            filter_nu=sub_nus[i] * 10 ** 9, filter_relative_bandwidth=sub_deltas[i] / sub_nus[i], detector_nep=2.7e-17
        )
        for i in range(Nbsubbands)
    ]
    sub_acqs = [
        QubicAcquisition(sub_instruments[i], sampling, scene, photon_noise=True, effective_duration=effective_duration)
        for i in range(Nbsubbands)
    ]
    covlim = 0.1
    coverage = np.array([sub_acqs[i].get_coverage() for i in range(Nbsubbands)])
    observed = [(coverage[i] > covlim * np.max(coverage[i])) for i in range(Nbsubbands)]
    obs = reduce(logical_and, tuple(observed[i] for i in range(Nbsubbands)))
    pack = PackOperator(obs, broadcast="rightward")

    ######################
    ### Reconstruction ###
    ######################
    sub_acqs_restricted = [sub_acqs[i][:, :, obs] for i in range(Nbsubbands)]
    operators = np.array([sub_acqs_restricted[i].get_operator() for i in range(Nbsubbands)])
    K = np.array([np.sum([operators[j] for j in bands_numbers[i]], axis=0) for i in range(Nbbands)])
    H = BlockRowOperator([K[i] for i in range(Nbbands)], new_axisin=0)
    invntt = sub_acqs[0].get_invntt_operator()
    A = H.T * invntt * H
    b = (H.T * invntt)(Y)
    preconditioner = BlockDiagonalOperator(
        [DiagonalOperator(1 / coverage[0][obs], broadcast="rightward") for i in range(Nbbands)], new_axisin=0
    )
    solution_qubic = pcg(A, b, M=preconditioner, disp=True, tol=1e-3, maxiter=1000)
    blockpack = BlockDiagonalOperator([pack for i in range(Nbbands)], new_axisin=0)
    maps = blockpack.T(solution_qubic["x"])
    maps[:, ~obs] = 0

    ####################################
    ### Monochromatic reconstruction ###
    ####################################

    x0 = np.zeros((Nbsubbands, Nbpixels, 3))
    for i in range(Nbsubbands):
        # x0[i,:,0]=cmb.T[0]+dust.T[0]*scaling_dust(150,sub_nus[i]e-9,1.59)
        x0[i, :, 1] = cmb.T[1] + dust.T[1] * scaling_dust(150, sub_nus[i], 1.59)
        x0[i, :, 2] = cmb.T[2] + dust.T[2] * scaling_dust(150, sub_nus[i], 1.59)

    maps_mono = np.zeros((Nbbands, Nbpixels, 3))
    if return_mono:
        (m, n) = shape(Y)
        Y_mono = np.zeros((Nbbands, m, n))
        for i in range(Nbbands):
            for j in bands_numbers[i]:
                C = HealpixConvolutionGaussianOperator(
                    fwhm=sub_instruments[j].synthbeam.peak150.fwhm * (150 / (sub_nus[j]))
                )
                Y_mono[i] = Y_mono[i] + operators[j] * pack * C * x0[j]
            Global_instrument = QubicInstrument(
                filter_nu=nus[i], filter_relative_bandwidth=deltas[i] / nus[i], detector_nep=2.7e-17
            )
            Global_acquisition = QubicAcquisition(
                Global_instrument, sampling, scene, photon_noise=True, effective_duration=effective_duration
            )
            noise = Global_acquisition.get_noise()
            Y_mono[i] = Y_mono[i] + noise
            H_mono = np.sum([operators[j] for j in bands_numbers[i]], axis=0)
            A_mono = H_mono.T * invntt * H_mono
            b_mono = (H_mono.T * invntt)(Y_mono[i])
            preconditioner_mono = DiagonalOperator(1 / coverage[0][obs], broadcast="rightward")
            solution_qubic_mono = pcg(A_mono, b_mono, M=preconditioner_mono, disp=True, tol=1e-3, maxiter=1000)
            maps_mono[i] = pack.T(solution_qubic_mono["x"])
            maps_mono[i, ~obs] = 0
        return maps, maps_mono, bands, deltas

    return maps, bands, deltas

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_noise [as 别名]
    [racenter, deccenter], duration, ts, angspeed, delta_az, nsweeps_el,
    angspeed_psi, maxpsi)
sampling.angle_hwp = np.random.random_integers(0, 7, len(sampling)) * 11.25

# get the acquisition model
acquisition = QubicAcquisition(150, sampling,
                               nside=nside,
                               synthbeam_fraction=0.99,
                               detector_tau=0.01,
                               detector_sigma=1.,
                               detector_fknee=1.,
                               detector_fslope=1)

# simulate the timeline
tod, x0_convolved = map2tod(acquisition, x0, convolution=True)
tod_noisy = tod + acquisition.get_noise()

# reconstruct using two methods
map_all, cov_all = tod2map_all(acquisition, tod, tol=1e-2)
map_each, cov_each = tod2map_each(acquisition, tod, tol=1e-2)


# some display
def display(map, cov, msg, sub):
    for i, (kind, lim) in enumerate(zip('IQU', [200, 10, 10])):
        map_ = map[..., i].copy()
        mask = cov == 0
        map_[mask] = np.nan
        hp.gnomview(map_, rot=center, reso=5, xsize=400, min=-lim, max=lim,
                    title=msg + ' ' + kind, sub=(3, 3, 3 * (sub-1) + i+1))

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_noise [as 别名]
        acq = QubicAcquisition(instrument_one, sampling, scene, photon_noise=False)
        C = acq.get_convolution_peak_operator()
        coverage = acq.get_coverage()
        observed = coverage > covlim * np.max(coverage)
        acq_restricted = acq[:, :, observed]
        H = acq_restricted.get_operator()
        x0_convolved = C(input_maps)
        pack = PackOperator(observed, broadcast='rightward')
        y_noiseless = H(pack(x0_convolved))
        invntt = acq.get_invntt_operator()
        A = H.T * invntt * H
        preconditioner = DiagonalOperator(1 / coverage[observed], broadcast='rightward')
        
        for nn in xrange(nbmc):
            noise = acq.get_noise()/np.sqrt(real_duration*ndetectors)
            y = y_noiseless + noise
            b = (H.T * invntt)(y)
            solution_qubic = pcg(A, b, M=preconditioner, disp=True, tol=1e-3, maxiter=1000)
            maps = pack.T(solution_qubic['x'])
            maps[~observed] = 0
            x0_convolved[~observed,:]=0    
        

            residuals = maps-x0_convolved
            freq, ps = powspec_inst(ts, y_noiseless+noise)

            pxsize_arcmin2 = 4*pi*(180/pi)**2 / (12*nside**2) * 60**2
            sigs= np.std(residuals[observed,:], axis=0)*sqrt(pxsize_arcmin2)
            print(s, f, nn, sigs)
            allsigs[i,j,nn,:] = sigs

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_noise [as 别名]
    angspeed_psi, maxpsi)
    
### Prepare input / output data
nside = 64
scene = QubicScene(nside, kind='IQU')

noise_vals = logspace(-5,0,10)
valspix = []
valsstd = []
for n in noise_vals:
    instrument = QubicInstrument(filter_nu=150e9,
                    detector_nep=n*4.7e-17)

    acq = QubicAcquisition(instrument, sampling, scene, photon_noise=False, 
                            effective_duration=effective_duration)
    noise = acq.get_noise()
    valspix.append(std(noise, axis=1))
    valsstd.append(std(noise))

clf()
plot(noise_vals,valsstd)
xscale('log')

########## Important: mettre photon_noise = False dans QubicAcquisition quand on veut baisser le niveau de bruit !
########## Maintenant il y a le keyword effective_duration (en années) qui scale le bruit (valable uniquement sous l'hypothèse de bruit blanc hélas)




sampling = create_sweeping_pointings(
    [racenter, deccenter], duration, ts, angspeed, delta_az, nsweeps_el,