Python QubicAcquisition.get_invntt_operator方法代码示例

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_invntt_operator [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_invntt_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_invntt_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_invntt_operator [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),

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_invntt_operator [as 别名]
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,
                            max=lim, title=msg + ' ' + kind,