<jupyter_start><jupyter_text>This notebook times the execution of the galsim and scarlet resampling routines<jupyter_code>%pylab inline
data_dir=galsim.meta_data.share_dir
#Reference catalog for simulations
cat = galsim.COSMOSCatalog(dir=data_dir, file_name = 'real_galaxy_catalog_23.5_example.fits')
#Survey dictionaries
HST, EUCLID, ROMAN, HSC, RUBIN = load_surveys()
#CHannels
channel_hr = ['hr']
channel_lr = ['lr']
channels = channel_lr+channel_hr

# Shapes tested for timing
shapes = [200, 150, 100, 70, 50, 30, 20]
#List of surveys
surveys = [HST, EUCLID, ROMAN, HSC, RUBIN]
#Size of the psf (pixels on-a-side)
npsf = 41

matplotlib.rc('image', cmap='gist_stern')
matplotlib.rc('image', interpolation='none')
matplotlib.rc('xtick', labelsize=20) 
# Run galsim's sinc interpolation
sinc = 0

# Dictionary for storing results
structure = {'survey_hr':[], 'survey_lr': [], 'n_hr': [], 'n_lr': [], 's_mean': [], 's_std': [], 
             'g_mean': [], 'g_std': [], 'gs_mean': [], 'gs_std': []}


#Timing for all combinations of high-low resolution images and all shapes in `shapes`
for nhr in shapes:
    for k, surveyhr in enumerate(surveys):
        print('survey hr  ' + str(nhr), surveyhr)
        for g, surveylr in enumerate(surveys[k+1:]):
            print(surveylr)
            # Making sure that the lr image spans roughly the same area on the sky as the hr image.
            nlr = int(nhr*surveyhr['pixel']/surveylr['pixel'])
            
            # Make simulations, picks a galaxy at random
            r = np.random.rand(1)*90
            pic_hr, pic_lr = gct.mk_scene(surveyhr, 
                                          surveylr, 
                                          cat,
                                          (nhr, nhr), 
                                          (nlr, nlr), 
                                          1,
                                          'real',
                                          noise = False,
                                          random_seds=False,
                                          index = 1,
                                          use_cat = False)
            
            wcs_hr = pic_hr.wcs
            wcs_lr = pic_lr.wcs
        
            data_hr = pic_hr.cube[0]/np.sum(pic_hr.cube[0])
            data_lr = pic_lr.cube[0]/np.sum(pic_lr.cube[0])
            
        
            psf_hr = scarlet.ImagePSF(np.array(pic_hr.psfs)[0][None, :,:])
            psf_lr = scarlet.ImagePSF(np.array(pic_lr.psfs)[0][None, :,:])
        
            ref = wcs_hr.wcs.crval.reshape((1,2))
            angle = 0
            
            if g == 0:
                # scarlet setup
                obs, frame = gct.setup_scarlet(data_hr, data_lr, wcs_hr, wcs_lr, psf_hr, psf_lr, channels, 'intersection')
                obs_lr, obs_hr = obs
                renderer = ResolutionRenderer(obs_lr, frame)
             
            # Galsim setup:
            ## GSO from psf_hr for galsim
            psf_hr_galsim = galsim.InterpolatedImage(galsim.Image(psf_hr.get_model()[0]), 
                                               scale = surveyhr['pixel'], use_true_center = False)
            ## deconvolution kernel for diff kernel
            deconv = galsim.Deconvolve(psf_hr_galsim)
            ## Interpolation of low resolution psf at high resolution
            psf_lr_hr = galsim.InterpolatedImage(galsim.Image(psf_lr.get_model()[0]), 
                                               scale = surveylr['pixel'], use_true_center = False)
            ## Difference kernel from galsim
            diff_gal = galsim.Convolve(deconv, psf_lr_hr)
            
            if sinc == 0:
                # Scarlet timing
                t_s = %timeit -n 1000 -r 2 -o renderer(data_hr[None, :,:])
                s_mean = np.array(t_s.all_runs).mean()/t_s.loops
                s_std = np.array(t_s.all_runs).std()/t_s.loops 
                # Galsim timing
                t_g = %timeit -n 1000 -r 2 -o gct.interp_galsim(data_hr, data_lr, diff_gal, angle, surveyhr['pixel'], surveylr['pixel'])
                g_mean = np.array(t_g.all_runs).mean()/t_g.loops 
                g_std = np.array(t_g.all_runs).std()/t_g.loops
                
                #Storage
                structure['survey_hr'].append(surveyhr) 
                structure['survey_lr'].append(surveylr) 
                structure['n_hr'].append(nhr) 
                structure['n_lr'].append(nlr)
                structure['s_mean'].append(s_mean)
                structure['s_std'].append(s_std)
                structure['g_mean'].append(g_mean) 
                structure['g_std'].append(g_std)
            else:
                t_gs = %timeit -n 1000 -r 2 -o gct.interp_galsim_sinc(data_hr, data_lr, diff_gal, angle, surveyhr['pixel'], surveylr['pixel'])
                gs_mean = np.array(t_gs.all_runs).mean()/t_gs.loops 
                gs_std = np.array(t_gs.all_runs).std()/t_gs.loops
            
                #Storage
                structure['survey_hr'].append(surveyhr) 
                structure['survey_lr'].append(surveylr) 
                structure['n_hr'].append(nhr) 
                structure['n_lr'].append(nlr)
                structure['gs_mean'].append(gs_mean) 
                structure['gs_std'].append(gs_std)
            
            
    #Saving results
    if sinc == 0:
        afile = open('Timings_npsf='+str(npsf)+'.pkl', 'wb')
        pickle.dump(structure, afile)
        afile.close()
    else:
        afile = open('Timings_sinc_npsf='+str(npsf)+'.pkl', 'wb')
        pickle.dump(structure, afile)
        afile.close()
structure = pickle.load(open('Timings_npsf='+str(npsf)+'.pkl', 'rb'))
structure_sinc = pickle.load(open('Timings_sinc_npsf='+str(npsf)+'.pkl', 'rb'))

#Reading results
n_hrs = np.array(structure['n_hr'])
n_lrs = np.array(structure['n_lr'])
s_mean = (np.array(structure['s_mean']))
s_std = np.array(structure['s_std'])
g_mean = np.array(structure['g_mean'])
g_std = np.array(structure['g_std'])
gs_mean = np.array(structure_sinc['gs_mean'])
gs_std = np.array(structure_sinc['gs_std'])

matplotlib.rc('xtick', labelsize=20) 
matplotlib.rc('ytick', labelsize=20) 
plt.close()
#Plotting results. Mean timing of resampling runs as a function of N
plt.figure(0, figsize = (15,15))
plt.subplots(figsize=(10, 6))
#ax2.errorbar(n_lrs, s_mean, yerr = s_std, capsize = 2, fmt = 'o', label = 'scarlet')
sc1 = plt.scatter(n_hrs, s_mean,c = n_lrs, marker = 'o', cmap = 'winter')
#ax2.errorbar(n_lrs, g_mean, yerr = g_std, capsize = 2, fmt = 's', label = 'galsim')
sc2 = plt.scatter(n_hrs, g_mean,c = n_lrs, marker = 's', cmap = 'copper')

sc3 = plt.scatter(n_hrs, gs_mean,c = n_lrs, marker = 'x', cmap = 'summer')
plt.yscale('log')
cbar = plt.colorbar(sc1)
plt.colorbar(sc2)
plt.colorbar(sc3)
cbar.set_label('M (# of pixels)', rotation=270, labelpad = 20, fontsize = 20)
# Adding plotting parameters
plt.title('Mean timing of resampling runs', fontsize=20)
plt.xlabel('N (# of pixels)', fontsize=20)
plt.ylabel('t (s)', fontsize=20)
plt.ylim([0.001,1])
plt.show()

#Plotting results. Mean timing of resampling runs as a function of M
plt.figure(1, figsize = (15,15))
plt.subplots(figsize=(10, 6))
#ax2.errorbar(n_lrs, s_mean, yerr = s_std, capsize = 2, fmt = 'o', label = 'scarlet')
sc1 = plt.scatter(n_lrs, s_mean,c = n_hrs, marker = 'o', cmap = 'winter')
#ax2.errorbar(n_lrs, g_mean, yerr = g_std, capsize = 2, fmt = 's', label = 'galsim')
sc2 = plt.scatter(n_lrs, g_mean,c = n_hrs, marker = 's', cmap = 'copper')

sc3 = plt.scatter(n_lrs, gs_mean,c = n_hrs, marker = 'x', cmap = 'summer')
plt.yscale('log')
cbar = plt.colorbar(sc1)
plt.colorbar(sc2)
plt.colorbar(sc3)
cbar.set_label('N (# ofpixels)', rotation=270, fontsize = 20, labelpad = 20)
# Adding plotting parameters
plt.title('Mean timing of resampling runs', fontsize=20)
plt.xlabel('M (# of pixels)', fontsize=20)
plt.ylabel('t (s)', fontsize=20)
plt.ylim(0.001,1)
plt.show()

#Plotting results. Mean timing of resampling runs as a function of M
plt.figure(1, figsize = (15,15))
plt.subplots(figsize=(10, 6))
for n in np.unique(n_hrs):
    plt.errorbar(n_lrs[n_hrs == n], s_mean[n_hrs == n],yerr = s_std[n_hrs == n], fmt = 'o', label = 'N = '+str(n))
plt.yscale('log')
plt.legend()
# Adding plotting parameters
plt.title('Mean timing of resampling runs', fontsize=20)
plt.xlabel('M (# of pixels)', fontsize=20)
plt.ylabel('t (s)', fontsize=20)
    
plt.show()
#Timing ratios as a function of M
plt.figure(2, figsize = (15,15))
fig2, ax2 = plt.subplots(figsize=(10, 6))
sc1 = ax2.scatter(n_lrs, s_mean/g_mean,c = n_hrs, marker = 'o', cmap = 'winter')
#plot((np.arange(5,200)), (np.arange(5,200))/30, label = 'M/30')
plt.legend()
yscale('log')
cbar = colorbar(sc1)
plot([0,220], [1,1], '--k')
cbar.set_label('N (# of pixels)', rotation=270, fontsize = 20, labelpad = 20)
# Adding plotting parameters
ax2.set_title('Timing ratios of resampling runs', fontsize=20)
ax2.set_xlabel('M (# of pixels)', fontsize=20)
ax2.set_ylabel('time ratios', fontsize=20)
savefig('galsim_scarlet_ratio.png')
show()

#Timing ratios as a function of N
figure(3, figsize = (15,15))
fig2, ax2 = plt.subplots(figsize=(10, 6))
sc1 = ax2.scatter(n_hrs, s_mean/g_mean,c = n_lrs, marker = 'o', cmap = 'winter')
yscale('log')
cbar = colorbar(sc1)
cbar.set_label('M (# of pixels)', rotation=270, fontsize = 20, labelpad = 20)
plot([0,220], [1,1], '--k')
# Adding plotting parameters
ax2.set_title('Timing ratios of resampling runs', fontsize=14)
plt.legend()
ax2.set_xlabel('N (# of pixels)', fontsize=20)
ax2.set_ylabel('time ratios', fontsize=20)
show()

#Colour representation of the colour ratios
fig2, ax2 = plt.subplots(figsize=(10, 8))
sc1 = ax2.scatter(n_hrs, n_lrs,c = np.log10(s_mean/g_mean), vmax = 2, vmin = -2, marker = 'o', cmap = 'seismic')
plt.ylim([0,220])
colorbar(sc1)
# Adding plotting parameters
ax2.set_title('Timing ratios of resampling runs', fontsize=14)
ax2.set_xlabel('N (# of pixels)', fontsize=20)
ax2.set_ylabel('M (# of pixels)', fontsize=20)
plt.savefig('measured_ratios.png')
show()
#Comparison with ggalsim sinc
#Timing ratios as a function of M
plt.figure(2, figsize = (15,15))
fig2, ax2 = plt.subplots(figsize=(10, 6))
sc1 = ax2.scatter(n_lrs,gs_mean/g_mean,c = n_hrs, marker = 'o', cmap = 'winter')
#plot((np.arange(5,200)), (np.arange(5,200))/30, label = 'M/30')
plt.legend()
yscale('log')
cbar = colorbar(sc1)
plot([0,220], [1,1], '--k')
cbar.set_label('N (# of pixels)', rotation=270, fontsize = 20, labelpad = 20)
# Adding plotting parameters
ax2.set_title('Timing ratios of resampling runs', fontsize=20)
ax2.set_xlabel('M (# of pixels)', fontsize=20)
ax2.set_ylabel('time ratios', fontsize=20)
savefig('galsim_sinc_scarlet_ratio.png')
show()

#Timing ratios as a function of N
figure(3, figsize = (15,15))
fig2, ax2 = plt.subplots(figsize=(10, 6))
sc1 = ax2.scatter(n_hrs, gs_mean/g_mean,c = n_lrs, marker = 'o', cmap = 'winter')
yscale('log')
cbar = colorbar(sc1)
cbar.set_label('M (# of pixels)', rotation=270, fontsize = 20, labelpad = 20)
plot([0,220], [1,1], '--k')
# Adding plotting parameters
ax2.set_title('Timing ratios of resampling runs', fontsize=14)
plt.legend()
ax2.set_xlabel('N (# of pixels)', fontsize=20)
ax2.set_ylabel('time ratios', fontsize=20)
show()

#Colour representation of the colour ratios
fig2, ax2 = plt.subplots(figsize=(10, 8))
sc1 = ax2.scatter(n_hrs, n_lrs,c = np.log10(gs_mean/g_mean), vmax = 2, vmin = -2, marker = 'o', cmap = 'seismic')
plt.ylim([0,220])
colorbar(sc1)
# Adding plotting parameters
ax2.set_title('Timing ratios of resampling runs', fontsize=14)
ax2.set_xlabel('N (# of pixels)', fontsize=20)
ax2.set_ylabel('M (# of pixels)', fontsize=20)
plt.savefig('measured_sinc_ratios.png')
show()
#theoretical version:
def galsim_timing(M, N):
    return N**2*(6**2+4**2*np.log(x*N))+M**2*(3+np.log(M))
def scarlet_timing(M, N):
    return N**2*((M+1)*np.log(N) + M*(M*1))
mnms = arange(20,201,10)

x,y = np.meshgrid(mnms, mnms)

plt.figure(figsize = (10,8))
plt.ylim([0,220])
plt.scatter(x,y, c = np.log10(scarlet_timing(x,y)/galsim_timing(x,y)), vmax = 2, vmin = -2, cmap = 'seismic')
#scatter(n_hrs, n_lrs,c = np.log10(s_mean/g_mean), marker = 'x', vmax = 1, vmin = -1, cmap = 'seismic')
plt.xlabel('M', fontsize = '30')
plt.ylabel('N', fontsize = '30')
plt.title('theoretical timing ratios', fontsize = '30')
cbar = plt.colorbar()
cbar.set_label('time ratios', rotation=270, fontsize = 25, labelpad = 20)
plt.savefig('theory_timing.png')
plt.show()
#Theoretical timing ratio for dominant components 
plt.figure(figsize = (10,8))
plt.ylim([0,220])
plt.scatter(x,y, c = np.log10(scarlet_timing(x,y)), cmap = 'seismic')
#scatter(n_hrs, n_lrs,c = np.log10(s_mean/g_mean), marker = 'x', vmax = 1, vmin = -1, cmap = 'seismic')
plt.xlabel('M', fontsize = '30')
plt.ylabel('N', fontsize = '30')
plt.title('dominant timing ratios', fontsize = '30')
cbar = plt.colorbar()
cbar.set_label('time ratios', rotation=270, fontsize = 25, labelpad = 20)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15949834577210126, 'fraction_numerical': 0.03346926213741633, 'mean_word_length': 3.0746081504702194, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '425310', 'license': 'MIT', 'n_tokens_mistral': 5280, 'n_tokens_neox': 4875, 'n_words': 1025, 'path': 'Timing_test/Timing_comparison.ipynb', 'repo_name': 'herjy/Multi_resolution_comparisons', 'size': '12997'}