;+ ; NAME: ; mem ; PURPOSE: ; Implement a Maximum Entropy Algorithm within OVSA's ; Solar Microwave Imaging Software ; CATEGORY: ; OVSA APC imaging ; CALLING SEQUENCE: ; pro mem, cmap, cln, uv, vis, wgt, nmap, misc ; INPUTS: ; cmap CLEAN map array of [m, m] dimension, ; where m is the number of pixels in one direction. ; Temperatures are in [MK]. ; cln Structure returned by clean.pro. ; uv 2D complex array of uv coordinate in [time, baseline] form. ; The real part is the west component and the imaginary part is ; the north component. ; vis 2D complex visibility array in [time, baseline] form ; wgt 2D complex array of the weight for each uv point ; in [time, baseline] form. The real part is for RCP ; and the imaginary part is for LCP. ; OPTIONAL (KEYWORD) INPUT PARAMETERS: ; model Model map array of the same size as cmap. Temperatures are in [MK]. ; Default models are cmap truncated at tmin. ; flux Flux array of each frequency map in Solar Flux Unit (S.F.U.). ; Default values are model map fluxes. ; ferr Array of acceptable flux errors in S.F.U. ; Default values are tol * flux. ; tmin Array of the minimum temperature acceptable as model in [MK]. ; Dafault values are 0.01. ; tol Small number used in convergence criteria. Default value is 0.03 ; savefile File name to save the result. ; Saved file can be read using readmou.pro ; result = readmou(savefile). ; result is anonymous structure ; f_ghz float[1] Frequency in [GHz]. ; tb_xy float[m, m] Temperature map in [MK]. ; xyint float[1] Pixel size in [arcsec]. ; bmin float Minor axis of the CLEAN beam in pixels. ; bmaj float Major axis of the CLEAN beam in pixels. ; pa float Position angle of the CLEAN beam in degree. ; Counter clockwise from the y axis. ; alp float[1] Lagrange multiplier alpha. ; bet float[1] Lagrange multiplier beta. ; continue Set this keyword to resumed the aborted process. ; Maps and Lagrange multipliers are restored from savefile ; but other inputs are assumed to be the same as ; that of the aborted precess. ; ROUTINES CALLED: ; mem_get_input.pro ; readmou.pro ; ssmem_uvgrid.pro ; ssmem_grd_chi2.pro ; ssmem_show_uv.pro ; ssmem_chi2chk.pro ; ssmsm_control.pro ; ssmem_optimiz.pro ; fft_xy2uv.pro ; plot_3maps.pro ; plot_vis2vis_simpl.pro ; savemou.pro ; OUTPUTS: ; nmap MEM map of the same size as cmap. Temperatures are in [MK]. ; misc String array containing informative values. ; misc = ; ['Remark : Maximum Entropy Method', $ ; 'Remark : MEM, frq[GHz] = (frequency in [GHz])', $ ; 'Remark : MEM, iterat = (iterations spent in optimization)', $ ; 'Remark : MEM, chi2_nrm = (Actual / Expected)', $ ; 'Remark : MEM, flux_nrm = (Actual / Expected)', $ ; 'Remark : MEM, alpha = (Lagrange multipliers)', $ ; 'Remark : MEM, beta = (Lagrange multipliers)'] ; COMMENTS ; This is a modified program of ssmem.pro that is currently compatible with imagr.pro ; The internal process is essentially same as SSMEM. ; To use in imagr.pro, remove mem_get_input procedure at the beginning. ; Otherwise, all the inputs are retrieved from ssmem_get_input.pro. ; SIDE EFFECTS: ; ; RESTRICTIONS: ; MODIFICATION HISTORY: ; Written 2001 by Su-Chan Bong ;- pro mem, cmap, cln, uv, vis, wgt, nmap, misc, $ model = mmap, flux = flux, ferr = ferr, tmin = tmin_, $ tol = tol_, savefile = sfile, continue = cont_ mem_get_input, cmap, cln, uv, vis, wgt, nmap, misc, $ model = mmap, flux = flux, ferr = ferr, tmin = tmin_, chk = chk if chk eq 0 then return ;backup current display setting device, get_decomposed = dc0 device, decomposed = 0 tvlct, r0, g0, b0, /get multi0 = !p.multi ch0 = !p.charsize f0 = !p.font xm0 = !x.margin ym0 = !y.margin win0 = !d.window ;set parameters if keyword_set(tol) then tol = tol_ else tol = 3.e-2 sz = size(cmap) m = sz[1] ;map size at the highest frequency m2 = m / 2l f = [cln.f_ghz] ;frequency array nf = n_elements(f) nf1 = nf - 1l xyint = [cln.xyint] tb2flux = 3.6e-5 * (xyint * f)^2 if keyword_set(tmin_) then tmin = tmin_ else tmin = replicate(0.01, nf) spwgt = 0. spset = 0b cbm = [cln.bmin / xyint[nf1], cln.bmaj / xyint[nf1], cln.pa] sfset = keyword_set(sfile) cont = keyword_set(cont_) pvisobs = ptrarr(nf) ptmap = ptrarr(nf) pmmap = ptrarr(nf) nul = fltarr(nf) chi2 = fltarr(nf) chi2mmap = fltarr(nf) info = {iter : 0, rchi2 : nul, rflux : nul, a : nul, b : nul, $ perror : ptr_new([0l]), control : 0l, tol : tol, init : 1b} if keyword_set(mmap) $ then for k = 0l, nf1 do pmmap[k] = ptr_new(mmap[*, *, k] > tmin[k]) $ else for k = 0l, nf1 do pmmap[k] = ptr_new(cmap[*, *, k] > tmin[k]) if keyword_set(flux) then begin tflux = flux / tb2flux tferr = tflux * tol if keyword_set(ferr) then tferr = ferr / TB2flux > tferr endif else begin tflux = fltarr(nf) for k = 0l, nf1 do tflux[k] = total(*pmmap[k], /double) tferr = tflux * tol endelse if cont and sfset then begin mou = readmou(sfile) for k = 0l, nf1 do ptmap[k] = ptr_new(mou.tb_xy[*, *, k]) info.a = mou.alp info.b = mou.bet info.init = 0b mou = 0. endif else for k = 0l, nf1 do ptmap[k] = ptr_new(*pmmap[k]) for k = 0l, nf1 do begin visobsk = ssmem_uvgrid(uv[*, *, k], vis[*, *, k] / tb2flux[k], $ wgt[*, *, k] * tb2flux[k]^2, m, xyint[k]) chi2[k] = n_elements(visobsk.iu) pvisobs[k] = ptr_new(visobsk) dummy = ssmem_grd_chi2(*pmmap[k], visobsk, chi2mem = chi2memk) chi2mmap[k] = chi2memk endfor cherr = chi2 * tol > sqrt(2. * chi2) ;check validity chk = ssmem_show_uv(m, xyint[nf1], uv[*, *, nf1], *pvisobs[nf1]) if not chk then begin ptr_free, pmmap, ptmap, pvisobs return endif chi2chk = bytarr(nf) for k = 0l, nf1 do chi2chk[k] = ssmem_chi2chk(*pvisobs[k], *ptmap[k], tol) bad = where(chi2mmap le chi2, nbad) if nbad gt 0l then begin error = [3l, nbad, bad, 0l] ;error code 3, invalid input *info.perror = error print, 'Error: Too large data error.' goto, skip endif bad = where(chi2chk, nbad) if nbad gt 0l then begin error = [3l, nbad, bad, 0l] ;error code 3, invalid input *info.perror = error print, 'Error: Too small map size.' goto, skip endif entparm = {pmmap : pmmap, spset : spset, spwgt : spwgt, $ f : f, nf : nf, fqwgt : 0.} ;make control panel control = ssmem_control(init = {f : f, m : m}, /single) info.control = control.tlb ;main algorithm ssmem_optimiz, ptmap, chi2, cherr, tflux, tferr, pvisobs, entparm, info ;destroy control panel status = ssmem_control(control.tlb, /destroy) skip: ;compare observation and map visibility sz = size(vis) ntb = sz[1] * sz[2] cvis = complexarr(ntb, nf) + complex(!values.f_nan, !values.f_nan) uvint = 206265. / (xyint * m) u = reform(float(uv), ntb, nf) v = reform(imaginary(uv), ntb, nf) good = where(finite(vis[*, *, 0])) cvis[good, 0] = fft_xy2uv(*ptmap[0], tb2flux[0], $ u[good, 0] / uvint[0] + m2, v[good, 0] / uvint[0] + m2) nmap = *ptmap[0l] ptr_free, pmmap, ptmap, pvisobs big_tit ='I M A G R / M E M (SCBong)' tit3 = ['CLEAN MAP','MEM MAP','Parameters'] map3 = fltarr(m, m, 3) map3[*, *, 0] = cmap map3[*, *, 1] = nmap winsiz = (2.5 * max(abs(cln.box)) < m2) * xyint[nf1] ; [arcsec] txt = ['Maximum Entropy Method', $ 'Itera = ' + strcompress(string(info.iter), /remove_all), $ 'Rel Chi^2 = ' + strcompress(string(info.rchi2), /remove_all), $ 'Rel Flux = ' + strcompress(string(info.rflux), /remove_all), $ 'Alpha = ' + strcompress(string(info.a), /remove_all), $ 'Beta = ' + strcompress(string(info.b), /remove_all)] plot_3maps, map3, tit3, big_tit, xyint[nf1], winsiz, txt = txt plot_vis2vis_simpl, vis, cvis ;organize outputs: if keyword_set(sfile) then $ savemou, sfile, f, nmap, xyint, cbm, info.a, info.b error = *info.perror ;error information. [code, nbad, [bad....]]. case error[0] of 0l: print, 'Successful end: MEM' 1l: print, 'Maximization failure: MEM' 2l: print, 'Aborted by user: MEM' 3l: print, 'Inconsistent constraint: MEM' 4l: print, 'Optimization failure: MEM' endcase i0 = 1l nsgl = error[i0] if nsgl gt 0l then begin print, 'Error: Inconsistency detected.' i0 = i0 + 2l endif else i0 = i0 + 1l nbad = error[i0] if nbad gt 0l then begin print, 'Warning: Constraint not satisfied.' endif misc = $ ['Remark : Maximum Entropy Method', $ 'Remark : MEM, iterat = ' + string(info.iter), $ 'Remark : MEM, chi2_nrm = ' + string(info.rchi2), $ 'Remark : MEM, flux_nrm = ' + string(info.rflux), $ 'Remark : MEM, alpha = ' + string(info.a), $ 'Remark : MEM, beta = ' + string(info.b)] ptr_free, info.perror ;restore previous display setting wset, win0 !y.margin = ym0 !x.margin = xm0 !p.font = f0 !p.charsize = ch0 !p.multi = multi0 tvlct, r0, g0, b0 device, decomposed = dc0 end