Python QubicAcquisition.get_operator方法代碼示例

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_operator [as 別名]
 def func(scene, max_nbytes, ref1, ref2, ref3, ref4, ref5, ref6):
     acq = QubicAcquisition(instrument, sampling, scene,
                            max_nbytes=max_nbytes)
     sky = np.ones(scene.shape)
     H = acq.get_operator()
     actual1 = H(sky)
     assert_same(actual1, ref1, atol=10)
     actual2 = H.T(actual1)
     assert_same(actual2, ref2, atol=10)
     actual2 = (H.T * H)(sky)
     assert_same(actual2, ref2, atol=10)
     actual3, actual4 = tod2map_all(acq, ref1, disp=False)
     assert_same(actual3, ref3, atol=10000)
     assert_same(actual4, ref4)
     actual5, actual6 = tod2map_each(acq, ref1, disp=False)
     assert_same(actual5, ref5, atol=1000)
     assert_same(actual6, ref6)

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_operator [as 別名]
 def func(sampling, kind, sky, ref1, ref2, ref3, ref4, ref5, ref6):
     nprocs_instrument = max(size // 2, 1)
     acq = QubicAcquisition(instrument, sampling, kind=kind,
                            nprocs_instrument=nprocs_instrument)
     assert_equal(acq.comm.size, size)
     assert_equal(acq.instrument.detector.comm.size, nprocs_instrument)
     assert_equal(acq.sampling.comm.size, size / nprocs_instrument)
     H = acq.get_operator()
     invntt = acq.get_invntt_operator()
     tod = H(sky)
     #actual1 = acq.unpack(H(sky))
     #assert_same(actual1, ref1, atol=20)
     actual2 = H.T(invntt(tod))
     assert_same(actual2, ref2, atol=20)
     actual2 = (H.T * invntt * H)(sky)
     assert_same(actual2, ref2, atol=20)
     actual3, actual4 = tod2map_all(acq, tod, disp=False, maxiter=2)
     assert_same(actual3, ref3, atol=20)
     assert_same(actual4, ref4, atol=20)

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_operator [as 別名]
def test():
    kinds = 'I', 'IQU'
    instrument = QubicInstrument(synthbeam_dtype=float)[:400]
    np.random.seed(0)
    sampling = create_random_pointings([0, 90], 30, 5)
    skies = np.ones(12 * 256**2), np.ones((12 * 256**2, 3))

    def func(sampling, kind, sky, ref1, ref2, ref3, ref4, ref5, ref6):
        nprocs_instrument = max(size // 2, 1)
        acq = QubicAcquisition(instrument, sampling, kind=kind,
                               nprocs_instrument=nprocs_instrument)
        assert_equal(acq.comm.size, size)
        assert_equal(acq.instrument.detector.comm.size, nprocs_instrument)
        assert_equal(acq.sampling.comm.size, size / nprocs_instrument)
        H = acq.get_operator()
        invntt = acq.get_invntt_operator()
        tod = H(sky)
        #actual1 = acq.unpack(H(sky))
        #assert_same(actual1, ref1, atol=20)
        actual2 = H.T(invntt(tod))
        assert_same(actual2, ref2, atol=20)
        actual2 = (H.T * invntt * H)(sky)
        assert_same(actual2, ref2, atol=20)
        actual3, actual4 = tod2map_all(acq, tod, disp=False, maxiter=2)
        assert_same(actual3, ref3, atol=20)
        assert_same(actual4, ref4, atol=20)
        #actual5, actual6 = tod2map_each(acq, tod, disp=False)
        #assert_same(actual5, ref5, atol=1000)
        #assert_same(actual6, ref6)

    for kind, sky in zip(kinds, skies):
        acq = QubicAcquisition(instrument, sampling, kind=kind,
                               comm=MPI.COMM_SELF)
        assert_equal(acq.comm.size, 1)
        H = acq.get_operator()
        invntt = acq.get_invntt_operator()
        tod = H(sky)
        ref1 = acq.unpack(tod)
        ref2 = H.T(invntt(tod))
        ref3, ref4 = tod2map_all(acq, tod, disp=False, maxiter=2)
        ref5, ref6 = None, None #tod2map_each(acq, tod, disp=False)
        yield (func, sampling, kind, sky, ref1, ref2, ref3, ref4, ref5, ref6)

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_operator [as 別名]
def test():
    instrument = QubicInstrument()[:10]
    sampling = create_random_pointings([0, 90], 30, 5)
    nside = 64
    scenes = QubicScene(nside, kind='I'), QubicScene(nside)

    def func(scene, max_nbytes, ref1, ref2, ref3, ref4, ref5, ref6):
        acq = QubicAcquisition(instrument, sampling, scene,
                               max_nbytes=max_nbytes)
        sky = np.ones(scene.shape)
        H = acq.get_operator()
        actual1 = H(sky)
        assert_same(actual1, ref1, atol=10)
        actual2 = H.T(actual1)
        assert_same(actual2, ref2, atol=10)
        actual2 = (H.T * H)(sky)
        assert_same(actual2, ref2, atol=10)
        actual3, actual4 = tod2map_all(acq, ref1, disp=False)
        assert_same(actual3, ref3, atol=10000)
        assert_same(actual4, ref4)
        actual5, actual6 = tod2map_each(acq, ref1, disp=False)
        assert_same(actual5, ref5, atol=1000)
        assert_same(actual6, ref6)

    for scene in scenes:
        acq = QubicAcquisition(instrument, sampling, scene)
        sky = np.ones(scene.shape)
        nbytes_per_sampling = acq.get_operator_nbytes() // len(acq.sampling)
        H = acq.get_operator()
        ref1 = H(sky)
        ref2 = H.T(ref1)
        ref3, ref4 = tod2map_all(acq, ref1, disp=False)
        ref5, ref6 = tod2map_each(acq, ref1, disp=False)
        for max_sampling in 10, 29, 30:
            max_nbytes = None if max_sampling is None \
                              else max_sampling * nbytes_per_sampling
            yield (func, scene, max_nbytes, ref1, ref2, ref3, ref4, ref5, ref6)

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_operator [as 別名]
delta_az = 15.
angspeed_psi = 0.1
maxpsi = 45.
nsweeps_el = 300
duration = 24   # hours
ts = 20         # seconds
sampling = create_sweeping_pointings(
    [racenter, deccenter], duration, ts, angspeed, delta_az, nsweeps_el,
    angspeed_psi, maxpsi)
scene = QubicScene(nside, kind='I')

# acquisition model
acq = QubicAcquisition(150, sampling, scene, synthbeam_fraction=0.99,
                       detector_fknee=fknee, detector_fslope=fslope,
                       detector_ncorr=ncorr)
H_ga = acq.get_operator()
C = acq.get_convolution_peak_operator()
H = H_ga * C

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

# noise
sigma = acq.instrument.detector.nep / np.sqrt(2 * sampling.period)
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()
noise[...] = 0

# map-making

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_operator [as 別名]
acq_qubic = QubicAcquisition(150, sampling, scene, effective_duration=1)
convolved_sky = acq_qubic.instrument.get_convolution_peak_operator()(sky)
acq_planck = PlanckAcquisition(150, acq_qubic.scene, true_sky=convolved_sky)
acq_fusion = QubicPlanckAcquisition(acq_qubic, acq_planck)

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=maxiter, tol=tol)

acq_qubic = QubicAcquisition(150, sampling, scene, effective_duration=1)
H = acq_qubic.get_operator()
invntt = acq_qubic.get_invntt_operator()
y, sky_convolved = acq_qubic.get_observation(sky, convolution=True)

A = H.T * invntt * H
b = H.T * invntt * y

solution_qubic = pcg(A, b, disp=True, maxiter=maxiter, tol=tol)


# some display
def display(input, msg, iplot=1):
    out = []
    for i, (kind, lim) in enumerate(zip('IQU', [50, 5, 5])):
        map = input[..., i]
        out += [hp.gnomview(map, rot=center, reso=5, xsize=800, min=-lim,

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_operator [as 別名]
center = gal2equ(center_gal[0], center_gal[1])

# sampling model
np.random.seed(0)
sampling = create_random_pointings(center, 1000, 10)

# scene model
scene = QubicScene(hp.npix2nside(x0.size), kind='I')

# instrument model
instrument = QubicInstrument(filter_nu=150e9)

# acquisition model
acq = QubicAcquisition(instrument, sampling, scene)
x0_convolved = acq.get_convolution_peak_operator()(x0)
H = acq.get_operator()
coverage = H.T(np.ones(H.shapeout))
mask = coverage > 0

# restrict the scene to the observed pixels
acq_restricted = acq[..., mask]
H_restricted = acq_restricted.get_operator()
x0_restricted = x0[mask]
y = H_restricted(x0_restricted)
invntt = acq_restricted.get_invntt_operator()

# solve for x
A = H_restricted.T * invntt * H_restricted
b = H_restricted.T(invntt(y))
solution = pcg(
    A, b, M=DiagonalOperator(1/coverage[mask]), disp=True, tol=1e-4)