IDL code to create the enormous correlation matrix plot from Sinha et al 2018

compute_joint_wp_and_gmf_correlation_matrix_from_theory_with_ngal_from_mocks.pro

function return_left_edge,  axis,   index,   value
  common tick_block,  fixed_value_left_x, fixed_value_right_x
  if value lt n_elements(fixed_value_left_x) then begin
     tickmark =  string_logexp(fixed_value_left_x[value])
  endif else begin
     if value eq n_elements(fixed_right_x) then begin
        tickmark =  string_logexp(fixed_value_right_x[-1])
     endif else begin
        tickmark =   'N!Dgal!N'
     endelse
  endelse

  RETURN,  tickmark[0]
end
  

compile_opt idl2, strictarrsubs
common tick_block
lum_sample = 21
nseeds = 50
wp_nbins = 14
nrots = 4
generate_ps = 1

;;; [-0.2,1.0] for Mr21 and [-0.35,1.0] for Mr19 to compare with subsampled. [-0.2,1.0] to compare between lum. samples
;; if lum_sample eq 21 then begin 
;;    cbrange =  [-0.2,  1.0]
;; endif else begin
;;    cbrange =  [-0.35,  1.0]
;; endelse
cbrange = [-0.4, 1.0]
ncolors = 256

cameron_mocks = 0
subsampled_mocks = 0

plot_mock_measurements = 1
highlight_stats = 1
plot_stat_extent_arrows =  1
plot_vertical_lines_with_arrows = 0 ;;; only takes effect when arrows are being plotted 

case lum_sample of
   19: BEGIN
      gmf_nbins = 8
      start_seed = 4001

      ;;; Cameron made the wp measurements for the -19
      if cameron_mocks eq 1 then begin
         cameron_wp = 1
         wp_datadir = '/data0/cameron/data/lasdamas/main19_manodeep/rppi_s0.2t20_pimax40by10/proj/'
         gmf_filebase = '/fs1/sinham/lasdamas/Consuelo/Mr19/best_box_fits/with_covariance/out2_SEED_irot_IROT_gmf_geometry.dat'
         ngal_file = '/data2/cameron/data/lasdamas/main19_manodeep/out2_ncount.dat'
         
      ;;; these are the subsampled files from the correlation matrix
      ;;; plots
      endif else begin
         cameron_wp = 0
         if subsampled_mocks eq 1 then begin
            dir =  '/fs1/sinham/lasdamas/Consuelo/Mr19/best_box_fits/with_covariance/subsampled/'
            wp_filebase  = dir + 'Consuelo_SEED_irot_IROT_geometry_subsampled.wp'
            gmf_filebase = dir + 'Consuelo_SEED_irot_IROT_mock_geometry_subsampled.grps'
            ngal_file = dir + 'makemock_ngal.dat'
         endif else begin
            dir =  '/fs1/sinham/lasdamas/Consuelo/Mr19/best_box_fits/joint_correlation_matrix/'
            wp_filebase  = dir + 'Consuelo_SEED_irot_IROT_geometry.wp'
            gmf_filebase = dir + 'Consuelo_SEED_irot_IROT_gmf_geometry.dat'
            ngal_file = dir + 'makemock_ngal.dat'
         endelse
      endelse
      
      sdss_wp_file = '/home/sinham/code/andreas/data/sdss_Mr19_wp_northonly.dat'
      sdss_gmf_file = '/home/sinham/code/andreas/data/sdss_Mr19_fib0_gmf_northonly.dat'
      sdss_ngal_file = '/home/sinham/code/andreas/data/sdss_Mr19_ngal_northonly.dat'
      volume = 5.555221d6 ;;approx volume of the mocks -- the real volume is (from Cameron through email) is 5.555221e+06
      gmf_tickv = [5, 10, 20, 50, 200]
   END

   21: BEGIN
      gmf_nbins = 6
      start_seed = 2001

      cameron_wp = 0
      if subsampled_mocks eq 0 then begin
         dir = "/fs1/sinham/lasdamas/Carmen/mocks_for_joint_covariance_matrix/"
         wp_filebase  =  dir + 'Carmen_SEED_irot_IROT_geometry.wp'
         gmf_filebase =  dir + 'Carmen_SEED_irot_IROT_gmf_geometry.dat'
      endif else begin
         dir = "/fs1/sinham/lasdamas/Carmen/mocks_for_joint_covariance_matrix/subsampled/"
         wp_filebase  =  dir + 'Carmen_SEED_irot_IROT_mock_geometry_subsampled.wp'
         gmf_filebase =  dir + 'Carmen_SEED_irot_IROT_mock_geometry_subsampled.grps'
      endelse
      
      ngal_file    =  dir + 'makemock_ngal.dat'
      sdss_wp_file   =  '/home/sinham/code/andreas/data/sdss_Mr21_wp_northonly.dat'
      sdss_gmf_file  =  '/home/sinham/code/andreas/data/sdss_Mr21_fib0_gmf_northonly.dat'
      sdss_ngal_file =  '/home/sinham/code/andreas/data/sdss_Mr21_ngal_northonly.dat'
      volume =  7.837407d7
      gmf_tickv =   [5,   10,   20,   40]
   END
   else:BEGIN
      print, 'ERROR: luminosity sample not supported'
      stop
   END
endcase
if ((cameron_wp ne 0) &&  (lum_sample ne 19)) then stop

statistics_label = ['n!Dgal!N', 'w!Dp!N(r!Dp!N)', 'n(N)']

seeds =  lindgen(nseeds) + start_seed
totnbins = wp_nbins+gmf_nbins+1
totnmocks = nseeds*nrots
if n_elements(model_data) eq 0 then $
   model_data = dblarr(totnbins, totnmocks)
sdss_data = dblarr(totnbins)

get_lun, lun
idat = lonarr(5)

stat_bin = 0
if n_elements(all_ngal) eq 0 then begin
   readcol, ngal_file, ngal, format = 'D', /silent
   if n_elements(ngal) ne nseeds*nrots then stop
   all_ngal =  ngal
   all_ngal /= volume
   model_data[stat_bin, 0:totnmocks-1L] = all_ngal
endif
stat_bin++

if n_elements(all_wp) eq 0 then begin
   all_wp = dblarr(wp_nbins, nseeds*nrots)
   index = 0
   for iseed = 0, nseeds-1 do begin
      for irot = 0, nrots-1 do begin
         if max(all_wp[*, index]) eq 0.0 then begin
            case lum_sample of 
               19: BEGIN
                  if cameron_wp eq 1 then begin
                     file = wp_datadir + 'rppi_s0.2t20_pimax40by10.out2_'+strn(seeds[iseed]) + '_irot_' + strn(irot) + '.wp.dat'
                     readcol, file, wp, format = 'D', /silent
                  endif else begin
                     file =  str_replace(wp_filebase,  'SEED',  string(seeds[iseed],  format =  '(I0)'))
                     file =  str_replace(file,  'IROT',  string(irot,  format =  '(I0)'))
                     readcol,  file,  wp,  format =  'X,X,X,D',  /silent ;;; rpmin,rpmax,rpavg,wp,npairs -> format. Skip rpmin,rpmax and rpavg
                  endelse
                  if n_elements(wp) ne wp_nbins then stop
               END
               21: BEGIN
                  file =  str_replace(wp_filebase,  'SEED',  string(seeds[iseed],  format =  '(I0)'))
                  file =  str_replace(file,  'IROT',  string(irot,  format =  '(I0)'))
                  readcol,  file,  wp,  format =  'X,X,X,D',  /silent ;;; rpmin,rpmax,rpavg,wp,npairs -> format. Skip rpmin,rpmax and rpavg
                  if n_elements(wp) ne wp_nbins then stop
               END
               ELSE: BEGIN
                  print, 'ERROR: luminosity sample not supported'
                  sto
               END
            endcase
            all_wp[0:wp_nbins-1L,  index] =   wp               
            model_data[stat_bin:stat_bin+wp_nbins-1L, index] = wp
         endif
         index++
      endfor
   endfor
endif
stat_bin += wp_nbins

if n_elements(all_gmf) eq 0 then begin
   all_gmf =  dblarr(gmf_nbins, nseeds*nrots)
   index = 0
   readcol, sdss_gmf_file, gmf_bin_left, gmf_bin_right, sdss_gmf, format = 'D,D,D', /silent
   if n_elements(gmf_bin_left) ne gmf_nbins or n_elements(gmf_bin_right) ne gmf_nbins then stop
   bin_width =  gmf_bin_right-gmf_bin_left
   for iseed = 0, nseeds-1 do begin
      for irot = 0, nrots-1 do begin
         if max(all_gmf[*, index]) eq 0.0 then begin
            file =  str_replace(gmf_filebase,  'SEED',  strn(seeds[iseed]))
            file =  str_replace(file,  'IROT',  strn(irot))
            readcol,  file,  ngalaxies, numgroups,  format =  'D,D',  /silent
            gmf =  dblarr(gmf_nbins)
            for ibin = 0, gmf_nbins-1 do begin
               ind =  where(ngalaxies ge gmf_bin_left[ibin] and ngalaxies lt gmf_bin_right[ibin], cnt)
               if cnt gt 0 then begin
                  gmf[ibin] =  total(numgroups[ind], /double)/volume/bin_width[ibin]
               endif
            endfor
            all_gmf[0:gmf_nbins-1L, index] =  gmf
            model_data[stat_bin:stat_bin+gmf_nbins-1L, index] = gmf
         endif
         index++
      endfor
      print, iseed, format = '("Done reading groups for seed # ",I3)'
   endfor
endif
correlation_matrix = make_correlation_matrix(model_data)
mincov =  min([correlation_matrix], max = maxcov)
print, mincov, maxcov
wait, 0.5

stat_bin = 0
readcol, sdss_ngal_file, sdss_ngal, sdss_ngal_error, format = 'D,D', /silent
sdss_data[stat_bin] = sdss_ngal
stat_bin++

readcol, sdss_wp_file, fixed_value_left_x, fixed_value_right_x,sdss_wp, format = 'D,D,x,D', /silent
sdss_data[stat_bin:stat_bin+wp_nbins-1L] = sdss_wp
stat_bin += wp_nbins

fixed_left_x  = lindgen(wp_nbins) + 1
fixed_right_x = fixed_left_x + 1
sdss_rp_left = fixed_value_left_x
sdss_rp_right = fixed_value_right_x

wp_tickv = [0.2,0.5, 1.0,2.0, 5.0, 10.0, 20.0]
wp_ticknames = string(wp_tickv, format = '(F0.1)')
quadterp, double(fixed_value_left_x), double(fixed_left_x), wp_tickv, wp_xtickv

readcol, sdss_gmf_file, ngal_left, ngal_right,sdss_gmf, format = 'D,D,D', /silent
sdss_data[stat_bin:stat_bin+gmf_nbins-1L] = sdss_gmf
stat_bin += gmf_nbins

fixed_left_x = lindgen(gmf_nbins) + 1 + wp_nbins
fixed_right_x = fixed_left_x + 1

fixed_value_left_x =  ngal_left
gmf_ticknames = string(gmf_tickv, format = '(I0)')
quadterp, double(fixed_value_left_x), fixed_left_x, gmf_tickv, gmf_xtickv

xtickv = [wp_xtickv, gmf_xtickv]
ticknames = [wp_ticknames, gmf_ticknames]

fixed_left_x  = lindgen(totnbins)
fixed_right_x = fixed_left_x + 1

size = 1200
ps_size = 16
if generate_ps eq 1 then begin
   set_plot,  'ps'
   ;; psfile =   'joint_count_in_bins_Mr'+strn(lum_sample)+'sample_bestHOD_from_amoeba_wp_with_ngal_from_mocks_gaussian_overlay.eps'
   ;; psfile =   'correlation_matrix_Mr'+strn(lum_sample)+'sample_bestHOD_from_amoeba_wp_with_ngal_from_mocks.eps'
   ;; psfile = 'joint_count_in_bins_Mr'+strn(lum_sample)+'sample_bestHOD_from_amoeba_wp_with_ngal_from_mocks_gaussian_overlay_subsampled.eps'
   psfile = '../figures/joint_CorrMatrix_Mr' + strn(lum_sample)
   if subsampled_mocks eq 1 then psfile =  psfile + '_subsampled'
   psfile = psfile + '.eps'
   device,    /decomposed,    bits_per_pixel =    8,    /color,    /tt_font,    /bold,    file =    psfile,    xsize =    ps_size,    ysize =    ps_size,    /inches,    /helvetica,   /encap
endif else begin
   window, xsize = size, ysize = size
endelse

tickinterval = 4
minor = tickinterval
xthick = 3
ythick = 3
thick = 6
charsize = 2
xticklen = 0.02
yticklen = 0.02
xrange = [0, totnbins]
yrange = xrange
cb_title = '[R!Di j!N]'
cb_title = tex2idl("[ $R_{ij}$ ]")
cb_title = ''
cb_charsize = 3
pos =   [0.15,   0.15,   0.85,   0.85]
cb_pos =  [0.15,  0.88,  0.85,  0.91]
if n_elements(ncolors) eq 0 then $
   ncolors = 256

cgloadct,   25,   /brewer,   /reverse,   ncolors =   ncolors
if n_elements(cbrange) eq 0 then $
   cbrange = [mincov, maxcov]

ms_plotimage_varbins, correlation_matrix, fixed_left_x, fixed_right_x, fixed_left_x, fixed_right_x, cbrange, logcolor = 0, $
                      xrange = yrange, yrange = yrange, $
                      xminor =   minor,   yminor =   minor,   xtickinterval =   tickinterval,   ytickinterval =   tickinterval,  $
                      xthick =   xthick,   ythick =   ythick,   xticklen =   xticklen,   yticklen =   yticklen,  $
                      xstyle =  5,  ystyle =  5,  xtickformat =  '(A1)',  ytickformat =  '(A1)', $
                      position = pos, ncolors = ncolors, charsize = charsize

;;; highlight the bins if requested
if ((n_elements(highlight_stats) ne 0) &&  (highlight_stats eq 1)) then begin

   ;;; First draw a thick square around ngal
   xarr = [0, 1, 1, 0, 0]
   yarr = [0, 0, 1, 1, 0]
   cgplots, xarr, yarr, thick = 2*thick, line = 0, color = 'black'

   xarr = [1, 1+wp_nbins, 1+wp_nbins, 1, 1]
   yarr = [1, 1, 1+wp_nbins, 1+wp_nbins, 1]
   cgplots, xarr, yarr, thick = 2*thick, line = 0, color = 'black'

   xarr = [1+wp_nbins, totnbins, totnbins, 1+wp_nbins, 1+wp_nbins]
   yarr = [1+wp_nbins, 1+wp_nbins, totnbins, totnbins, 1+wp_nbins]
   cgplots, xarr, yarr, thick = 2*thick, line = 0, color = 'black'

endif



if plot_mock_measurements eq 1 then begin
   for i =  0,  totnbins-1 do begin
      x1data = model_data[i, 0:totnmocks-1L]
      sdss_x1data = sdss_data[i]

      for j = 0, totnbins-1 do begin
         x2data = model_data[j, 0:totnmocks-1L]
         sdss_x2data = sdss_data[j]
         color =  cgcolor('black')
         newpos =  [i,  j,  i+1,  j+1]
         newpos_normal =   double(newpos)
         newpos_normal[0:1] =  (convert_coord(newpos[0:1],  /data,  /to_normal))[0:1]
         newpos_normal[2:3] =  (convert_coord(newpos[2:3],  /data,  /to_normal))[0:1]
         oldx =  !x
         oldy =  !y
         oldp =  !p

         
         if i eq j then begin
            if min(x1data) ne max(x1data) then begin
         ;;; Find the yrange
               histogauss, x1data, A, XX, YY,  GX, GY, /noplot, /double
               plot,  [0],  xrange =  minmax([XX, GX]),  yrange =  [0, max([A[0], YY])*1.02],  $
                      /noerase,  position =  newpos_normal,  xstyle =  5,  ystyle =  5,  $
                      /nodata
               cgplots, XX, YY, line = 0, color = !black, thick = thick, noclip = 0
               cgplots, GX, GY, line = 0, color = !white, thick = thick, noclip = 0
               quadterp, gx, gy, sdss_x1data, yloc
               cgplots, sdss_x1data, yloc, psym = 1, color = 'white', noclip = 0, symsize = 1.5
               print, i, A[4], format = '("Bin # ",I2," gaussianity = ",F0.2)'
            endif
         endif else begin
            plot,  [0],  xrange =  minmax(x1data),  yrange =  minmax(x2data),  $
                   /noerase,  position =  newpos_normal,  xstyle =  5,  ystyle =  5,  $
                   /nodata
            cgplots,  x1data,  x2data,  psym =  16,  color =  color,  noclip =  0, symsize = 0.4
            cgplots, sdss_x1data, sdss_x2data, psym = 1, color = 'white', noclip = 0, symsize = 1.5
         endelse
         !x =  oldx
         !y =  oldy
         !p =  oldp
      endfor
   endfor
endif
cgcolorbar,   range =   cbrange,  ncolors =  ncolors,  position = cb_pos, title = cb_title,  /top, charsize = cb_charsize
axis, xaxis = 0, xticks = n_elements(xtickv)-1, xtickv = xtickv, xthick = xthick, xticklen = xticklen, xtickname = ticknames, charsize = charsize
axis, yaxis = 0, yticks = n_elements(xtickv)-1, ytickv = xtickv, ythick = ythick, yticklen = yticklen, ytickname = ticknames, charsize = charsize
axis, xaxis = 1, xticks = n_elements(xtickv)-1, xtickv = xtickv, xthick = xthick, xticklen = xticklen, xtickname = ticknames, charsize = charsize, xtickformat = '(A1)'
axis, yaxis = 1, yticks = n_elements(xtickv)-1, ytickv = xtickv, ythick = ythick, yticklen = yticklen, ytickname = ticknames, charsize = charsize, ytickformat = '(A1)'

axis_label_pos =  -2.0
axis_label_charsize = 1.5*charsize
arrow_to_label_gap = 1
arrow_to_marker_gap = 0.05
if generate_ps eq 1 then begin
   arrow_hsize = !D.X_SIZE / 64.
endif else begin
   arrow_hsize = 10
endelse

arrow_thick = 0.3*thick
arrow_color = 'darkgray'

ngal_label = 'n!Dgal!N'
ngal_label_pos =  0.5*(0 + 1) ;;; mid-point of the entire ngal index range
ngal_align = 0.5
xyouts, ngal_label_pos, axis_label_pos, ngal_label, /data, alignment = ngal_align, charsize = axis_label_charsize
xyouts, axis_label_pos, ngal_label_pos,  ngal_label, /data, alignment = ngal_align, charsize = axis_label_charsize,  orientation = 90

wp_label =  'w!Dp!N(r!Dp!N)'
wp_start =  1  ;;; ngal is at index 0
wp_end = wp_start + wp_nbins
wp_label_pos =  0.5*(wp_start + wp_end) ;;; mid-point of the entire wp index range (wp[0] is at index 1, since ngal is at index 0)
xyouts, wp_label_pos, axis_label_pos, wp_label, /data, alignment = 0.5, charsize = axis_label_charsize
xyouts, axis_label_pos, wp_label_pos,  wp_label, /data, alignment = 0.5, charsize = axis_label_charsize,  orientation = 90

if plot_stat_extent_arrows eq 1 then begin
   cgarrow, wp_label_pos - arrow_to_label_gap,  axis_label_pos,  wp_start + arrow_to_marker_gap, axis_label_pos, /data, /solid, thick = arrow_thick,  hsize = arrow_hsize, color = arrow_color
   cgarrow, wp_label_pos + arrow_to_label_gap,  axis_label_pos,  wp_end - arrow_to_marker_gap, axis_label_pos, /data, /solid, thick = arrow_thick,  hsize = arrow_hsize, color = arrow_color
   cgarrow, axis_label_pos, wp_label_pos - arrow_to_label_gap,  axis_label_pos,  wp_start + arrow_to_marker_gap, /data, /solid, thick = arrow_thick,  hsize = arrow_hsize, color = arrow_color
   cgarrow, axis_label_pos, wp_label_pos + arrow_to_label_gap,  axis_label_pos, wp_end - arrow_to_marker_gap, /data, /solid, thick = arrow_thick,  hsize = arrow_hsize, color = arrow_color
endif

gmf_label = 'n(N)'
gmf_start = wp_end
gmf_end = gmf_start + gmf_nbins
gmf_label_pos = 0.5*(gmf_start + gmf_end) ;; mid-point of the gmf index range
xyouts, gmf_label_pos, axis_label_pos, gmf_label, /data, alignment = 0.5, charsize = axis_label_charsize
xyouts, axis_label_pos, gmf_label_pos,  gmf_label, /data, alignment = 0.5, charsize = axis_label_charsize, orientation = 90

if plot_stat_extent_arrows eq 1 then begin
   cgarrow, gmf_label_pos - arrow_to_label_gap,  axis_label_pos,  gmf_start + arrow_to_marker_gap, axis_label_pos, /data, /solid, thick = arrow_thick,  hsize = arrow_hsize, color = arrow_color
   cgarrow, gmf_label_pos + arrow_to_label_gap,  axis_label_pos,  gmf_end - arrow_to_marker_gap, axis_label_pos, /data, /solid, thick = arrow_thick,  hsize = arrow_hsize, color = arrow_color
   cgarrow, axis_label_pos, gmf_label_pos - arrow_to_label_gap,  axis_label_pos,  gmf_start + arrow_to_marker_gap, /data, /solid, thick = arrow_thick,  hsize = arrow_hsize, color = arrow_color
   cgarrow, axis_label_pos, gmf_label_pos + arrow_to_label_gap,  axis_label_pos, gmf_end - arrow_to_marker_gap, /data, /solid, thick = arrow_thick,  hsize = arrow_hsize, color = arrow_color
endif

;;; Now plot vertical lines marking the beginning and end of each stat
if (plot_stat_extent_arrows eq 1) && (plot_vertical_lines_with_arrows eq 1) then begin
   vertical_marker_size = 0.2
   vertical_marker_thick = 2*arrow_thick
   
   ypos = [axis_label_pos - vertical_marker_size, axis_label_pos + vertical_marker_size]
   cgplots,               [0, 0], ypos, /data, thick = vertical_marker_thick, color = arrow_color
   cgplots, [wp_start, wp_start], ypos, /data, thick = vertical_marker_thick, color = arrow_color
   cgplots,     [wp_end, wp_end], ypos, /data, thick = vertical_marker_thick, color = arrow_color
   cgplots,   [gmf_end, gmf_end], ypos, /data, thick = vertical_marker_thick, color = arrow_color
   
   ;;; on the Y-axis -> involves simply swapping the x-y from the previous 4 lines
   cgplots, ypos,               [0, 0],  /data, thick = vertical_marker_thick, color = arrow_color
   cgplots, ypos, [wp_start, wp_start],  /data, thick = vertical_marker_thick, color = arrow_color
   cgplots, ypos,     [wp_end, wp_end],  /data, thick = vertical_marker_thick, color = arrow_color
   cgplots, ypos,   [gmf_end, gmf_end],  /data, thick = vertical_marker_thick, color = arrow_color
endif

if !D.NAME eq 'PS' then begin
   device, /close
   set_plot, 'x'
   @idl_reset_graphics
endif

free_lun, lun
end