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

def tod2map(tod,pointing,instrument_in,detector_list=False,disp=True,kmax=2,displaytime=False):
    t0=time.time()
    #### Detectors
    mask_packed = np.ones(len(instrument_in.detector.packed), bool)
    if detector_list:
        mask_packed[detector_list] = False
        mask_unpacked = instrument_in.unpack(mask_packed)
        instrument = QubicInstrument('monochromatic', removed=mask_unpacked,nside=instrument_in.sky.nside)
    else:
        instrument = instrument_in
        
    #### Observations
    obs = QubicAcquisition(instrument, pointing)
    projection = obs.get_projection_peak_operator(kmax=kmax)
    coverage = projection.pT1()
    mask = coverage == 0
    projection = pack_projection_inplace(projection, mask)
    hwp = obs.get_hwp_operator()
    polgrid = DenseOperator([[0.5, 0.5, 0],
                             [0.5,-0.5, 0]])
    H = polgrid * hwp * projection
    preconditioner = DiagonalOperator(1/coverage[~mask], broadcast='rightward')
    solution = pcg(H.T * H, H.T(tod), M=preconditioner, disp=disp, tol=1e-3)
    output_map = unpack(solution['x'], mask)
    t1=time.time()
    if displaytime: print(' Map done in {0:.4f} seconds'.format(t1-t0))
    return output_map,coverage

def get_coverage_onedet(idetector=231, angspeed=1., delta_az=15., angspeed_psi=0., maxpsi=15., nsweeps_el=100, duration=24, ts=1., decrange=2., decspeed=2., recenter=False):
    print('##### Getting Coverage for: ') 
    print('## idetector = '+str(idetector))
    print('## duration = '+str(duration))
    print('## ts = '+str(ts))
    print('## recenter = '+str(recenter))
    print('## angspeed = '+str(angspeed))
    print('## delta_az = '+str(delta_az))
    print('## nsweeps_el = '+str(nsweeps_el))
    print('## decrange = '+str(decrange))
    print('## decspeed = '+str(decspeed))
    print('## angspeed_psi = '+str(angspeed_psi))
    print('## maxpsi = '+str(maxpsi))
    print('##########################')
    nside = 256
    racenter = 0.0
    deccenter = -57.0
    pointing = pointings_modbyJC.create_sweeping_pointings(
        [racenter, deccenter], duration, ts, angspeed, delta_az, nsweeps_el,
        angspeed_psi, maxpsi, decrange=decrange, decspeed=decspeed, recenter=recenter)
    pointing.angle_hwp = np.random.random_integers(0, 7, pointing.size) * 22.5
    ntimes = len(pointing)

    # get instrument model with only one detector
    instrument = QubicInstrument('monochromatic')
    mask_packed = np.ones(len(instrument.detector.packed), bool)
    mask_packed[idetector] = False
    mask_unpacked = instrument.unpack(mask_packed)
    instrument = QubicInstrument('monochromatic', removed=mask_unpacked)
    
    obs = QubicAcquisition(instrument, pointing)
    convolution = obs.get_convolution_peak_operator()
    projection = obs.get_projection_peak_operator(kmax=0)
    coverage = projection.pT1()
    return coverage

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

def test_add_subtract_grid_operator():
    acq = QubicAcquisition(150, sampling, detector_ngrids=2)
    tod = map2tod(acq, input_map, convolution=False)
    add = acq.get_add_grids_operator()
    sub = acq.get_subtract_grids_operator()
    assert_same(add(tod), tod[:992] + tod[992:])
    assert_same(sub(tod), tod[:992] - tod[992:])

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)    

def map2tod(maps,pointing,instrument_in,detector_list=False,kmax=2):
    #### Detectors
    mask_packed = np.ones(len(instrument_in.detector.packed), bool)
    if detector_list:
        mask_packed[detector_list] = False
        mask_unpacked = instrument_in.unpack(mask_packed)
        instrument = QubicInstrument('monochromatic', removed=mask_unpacked,nside=instrument_in.sky.nside)
    else:
        instrument = instrument_in
        
    #### Observations
    obs = QubicAcquisition(instrument, pointing)
    #C = obs.get_convolution_peak_operator()
    #convmaps=np.transpose(np.array([C(maps[:,0]),C(maps[:,1]),C(maps[:,2])]))
    projection = obs.get_projection_peak_operator(kmax=kmax)
    coverage = projection.pT1()
    mask = coverage == 0
    projection = pack_projection_inplace(projection, mask)
    hwp = obs.get_hwp_operator()
    polgrid = DenseOperator([[0.5, 0.5, 0],
                             [0.5,-0.5, 0]])
    #H = polgrid * hwp * projection * C
    H = polgrid * hwp * projection
    x1 = pack(maps, mask)
    y = H(x1)
    #return y,convmaps
    return y

def test_primary_beam():
    def primary_beam(theta, phi):
        import numpy as np
        with settingerr(invalid='ignore'):
            return theta <= np.radians(9)
    a = QubicAcquisition(150, s, primary_beam=primary_beam)
    p = a.get_projection_operator()
    assert_equal(p.matrix.ncolmax, 5)

def write_reference():
    np.random.seed(0)
    p = create_random_pointings([0, 90], NPTG, ANGLE)
    FitsArray([p.azimuth, p.elevation, p.pitch, p.angle_hwp]).save(FILEPTG)
    for kind, map, filename in zip(['I', 'IQU'], [MAPIQU[..., 0], MAPIQU],
                                   [FILETODI, FILETODIQU]):
        acq = QubicAcquisition(INSTRUMENT, p, nside=256, kind=kind)
        tod = acq.get_observation(map, noiseless=True)
        FitsArray(tod).save(filename)

def map2TOD(input_map, pointing,kmax=2):
    ns=hp.npix2nside(input_map.size)
    qubic = QubicInstrument('monochromatic,nopol',nside=ns)
    #### configure observation
    obs = QubicAcquisition(qubic, pointing)
    C = obs.get_convolution_peak_operator()
    P = obs.get_projection_peak_operator(kmax=kmax)
    H = P * C
    # Produce the Time-Ordered data
    tod = H(input_map)
    input_map_conv = C(input_map)
    return(tod,input_map_conv,P)

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)

 def func_observation(obs, tod, info):
     filename_ = os.path.join(outpath, 'obs-' + str(uuid1()))
     obs._save_observation(filename_, tod, info)
     obs2, tod2, info2 = QubicAcquisition._load_observation(filename_)
     assert_equal(str(obs), str(obs2))
     assert_equal(tod, tod2)
     assert_equal(info, info2)

def TOD2map(tod,pointing,nside,kmax=2,disp=True,P=False,covmin=10):
    qubic = QubicInstrument('monochromatic,nopol',nside=nside)
    #### configure observation
    obs = QubicAcquisition(qubic, pointing)
    if not P:
        P = obs.get_projection_peak_operator(kmax=kmax)
    coverage = P.T(np.ones_like(tod))
    mask=coverage < covmin
    P.matrix.pack(mask)
    P_packed = ProjectionInMemoryOperator(P.matrix)
    unpack = UnpackOperator(mask)
    # data
    solution = pcg(P_packed.T * P_packed, P_packed.T(tod), M=DiagonalOperator(1/coverage[~mask]), disp=disp)
    output_map = unpack(solution['x'])
    output_map[mask] = np.nan
    coverage[mask]=np.nan
    return(output_map,mask,coverage)

 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)

 def func_simulation(obs, input_map, tod, info):
     filename_ = os.path.join(outpath, 'simul-' + str(uuid1()))
     obs.save_simulation(filename_, input_map, tod, info)
     obs2, input_map2, tod2, info2 = QubicAcquisition.load_simulation(
         filename_)
     assert_equal(str(obs), str(obs2))
     assert_equal(input_map, input_map2)
     assert_equal(tod, tod2)
     assert_equal(info, info2)

 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)

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)

def test_load_save():
    info = 'test\nconfig'

    def func_acquisition(obs, info):
        filename_ = os.path.join(outpath, 'config-' + str(uuid1()))
        obs.save(filename_, info)
        obs2, info2 = QubicAcquisition.load(filename_)
        assert_equal(str(obs), str(obs2))
        assert_equal(info, info2)

    def func_observation(obs, tod, info):
        filename_ = os.path.join(outpath, 'obs-' + str(uuid1()))
        obs._save_observation(filename_, tod, info)
        obs2, tod2, info2 = QubicAcquisition._load_observation(filename_)
        assert_equal(str(obs), str(obs2))
        assert_equal(tod, tod2)
        assert_equal(info, info2)

    def func_simulation(obs, input_map, tod, info):
        filename_ = os.path.join(outpath, 'simul-' + str(uuid1()))
        obs.save_simulation(filename_, input_map, tod, info)
        obs2, input_map2, tod2, info2 = QubicAcquisition.load_simulation(
            filename_)
        assert_equal(str(obs), str(obs2))
        assert_equal(input_map, input_map2)
        assert_equal(tod, tod2)
        assert_equal(info, info2)

    for p in ptgs:
        p = np.asarray(p)
        smp = QubicSampling(azimuth=p[..., 0], elevation=p[..., 1],
                            pitch=p[..., 2])
        acq = QubicAcquisition(qubic, smp)
        P = acq.get_projection_peak_operator()
        tod = P(input_map)
        yield func_acquisition, acq, info
        yield func_observation, acq, tod, info
        yield func_simulation, acq, input_map, tod, info

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)

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