;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;; GETJUL TAKES TWO SETS OF DATES AND ;;;;;; ;;;;;;;;; TIMES AND RETURNS THE CORRESPONDING ;;;;;; ;;;;;;;;; JULIAN DATES IN DOUBLE PRECISION ;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function getjul, date, time n = n_elements(date) ;NUMBER OF ELEMENTS IN INPUT juldates = dblarr(n) ;SET UP OUTPUT cent = '20' ;AVOID Y2.1K!!! for i=0,n-1 do begin ;PARSE THE DATE IN M/D/Y FORMAT pos1 = strpos(date[i],'/') ;FIRST '/' POSITION pos2 = strpos(date[i],'/',pos1+1) ;SECOND '/' POSITION month = fix(strmid(date[i],0,pos1)) day = fix(strmid(date[i],pos1+1,(pos2-1)-pos1)) year = fix(cent+strmid(date[i],pos2+1)) ;PARSE THE TIME IN HH:MM FORMAT pos3 = strpos(time[i],':') hour = fix(strmid(time[i],0,pos3)) min = fix(strmid(time[i],pos3+1)) ;THIS TESTS THE PARSE... USES NKDSTRIN WHICH IS JUST STRCOMPRESS(,/REMOVE_ALL) ;print, nkdstring(month)+'/'+nkdstring(day)+'/'+nkdstring(year) & $ ;print, nkdstring(hour)+':'+nkdstring(min) juldate, [year,month,day,hour,min], jd juldates[i] = double(jd) + 2400000.d endfor return, juldates end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;; CALCULATES THE CUSUM FOR A GIVEN VECTOR ;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function cusum, input, value=VALUE, median=MEDIAN vector = input if n_elements(value) le 0 then value = mean(vector) if n_elements(median) gt 0 then value = median(vector) vector = vector - value ;NEED THE INITIAL VALUES MINUS THE MEAN n = n_elements(vector) csum = dblarr(n) ;THE CUSUM ARRAY ;ADD ALL THE MEAN-SUBTRACTED POINTS BEFORE A GIVEN POINT for i=0,n-1 do begin csum[i]=total(vector[0:i]) endfor return, csum end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;; CALCULATES BLOOD PRESSURE STATISTICS ;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro bp restore, 'bp.template' bpstruct = read_ascii('bp.txt', template = bptemplate) error = 3.0 ;ERROR IN PRESSURE plotps = 1 ;PLOT TO .PS IF SET plotbmp = 0 ;PLOT TO .BMP IF SET @circsym ;DEFINES CIRCULAR PLOTTING SYMBOL n = n_elements(bpstruct.date) date = bpstruct.date time = bpstruct.time sys = bpstruct.dys dys = bpstruct.sys hrt = bpstruct.hr juldates = getjul(date,time) res = sys-dys ;SYSTOLIC/DIASTOLIC RESIDUALS ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;RUNNING MEAN IN +/- "WIND" DAYS rmnd = fltarr(n) & rmns = rmnd & rmnh = rmns & rmnr = rmnh wind = 2.0 ;NUMBER OF DAYS FOR RUNNING MEAN for i=0,n-1 do begin g = where(juldates ge juldates[i]-wind and juldates le juldates[i]+wind) rmnd[i] = mean(dys[g]) rmns[i] = mean(sys[g]) rmnh[i] = mean(hrt[g]) rmnr[i] = mean(res[g]) endfor ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;CALCULATE CUSUM STATISTICS ON RUNNING MEANS csumd = cusum(dys) csums = cusum(sys) csumh = cusum(hrt) csumr = cusum(res) ;csumd = cusum(rmnd) ;csums = cusum(rmns) ;csumh = cusum(rmnh) ;csumr = cusum(rmnr) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; bp_sys ;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if plotps eq 1 then begin a=PSWINDOW(/landscape, /cm, margin=0.025) set_plot, 'ps' outputname= '~/bp/bp_sys.ps' device, file=outputname, /landscape endif else begin plot, [0],[1] & wshow endelse temp = LABEL_DATE(DATE_FORMAT='%N/%D/%Z') ch = 0.750 ;charsize sysmn = mean(sys) & sysmd = median(sys) & syssd = stddev(sys) meanstring = string(sysmn,format='(f5.1)') medistring = string(sysmd,format='(f5.1)') stdvstring = string(syssd,format='(f4.1)') xpos = juldates[0]+0.075*(abs(juldates[n-1]-juldates[0])) tit = 'Mean = '+meanstring+' Median = '+medistring+' Standard Deviation = '+stdvstring rn=j_yrange(sys) jtxt = rn[1]-0.075*(abs(rn[1]-rn[0])) plot, juldates, sys, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Systolic Pressure (mmHg)', $ yrange=rn, thick=2, psym=-8, symsize=ch, /xstyle, /ystyle xyouts, xpos, jtxt, tit, alignment=0.0, charsize=1.0 if plotps eq 1 then begin device, /close set_plot, 'x' endif else stop ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; bp_dys ;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if plotps eq 1 then begin a=PSWINDOW(/landscape, /cm, margin=0.025) set_plot, 'ps' outputname= '~/bp/bp_dys.ps' device, file=outputname, /landscape endif else begin plot, [0],[1] & wshow endelse dysmn = mean(dys) & dysmd = median(dys) & dyssd = stddev(dys) meanstring = string(dysmn,format='(f4.1)') medistring = string(dysmd,format='(f4.1)') stdvstring = string(dyssd,format='(f3.1)') tit = 'Mean = '+meanstring+' Median = '+medistring+' Standard Deviation = '+stdvstring rn=j_yrange(dys) jtxt = rn[1]-0.075*(abs(rn[1]-rn[0])) plot, juldates, dys, psym=-8, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Diastolic Pressure (mmHg)',yrange=rn, $ thick=2, symsize=ch, /xstyle, /ystyle xyouts, xpos, jtxt, tit, alignment=0.0, charsize=1.0 if plotps eq 1 then begin device, /close set_plot, 'x' endif else stop ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; bp_hrt ;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if plotps eq 1 then begin a=PSWINDOW(/landscape, /cm, margin=0.025) set_plot, 'ps' outputname= '~/bp/bp_hrt.ps' device, file=outputname, /landscape endif else begin plot, [0],[1] & wshow endelse hrmn = mean(hrt) & hrmd = median(hrt) & hrsd = stddev(hrt) meanstring = string(hrmn,format='(f4.1)') medistring = string(hrmd,format='(f4.1)') stdvstring = string(hrsd ,format='(f4.1)') tit = 'Mean = '+meanstring+' Median = '+medistring+' Standard Deviation = '+stdvstring rn=j_yrange(hrt) jtxt = rn[1]-0.075*(abs(rn[1]-rn[0])) plot, juldates, hrt, psym=-8, yrange=rn, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Heart Rate (beats /min)', thick=2, symsize=ch, /xstyle, /ystyle xyouts, xpos, jtxt, tit, alignment=0.0, charsize=1.0 if plotps eq 1 then begin device, /close set_plot, 'x' endif else stop ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; bp_res ;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if plotps eq 1 then begin a=PSWINDOW(/landscape, /cm, margin=0.025) set_plot, 'ps' outputname= '~/bp/bp_res.ps' device, file=outputname, /landscape endif else begin plot, [0],[1] & wshow endelse resmn = mean(res) & resmd = median(res) & ressd = stddev(res) meanstring = string(resmn,format='(f4.1)') medistring = string(resmd,format='(f4.1)') stdvstring = string(ressd ,format='(f3.1)') tit = 'Mean = '+meanstring+' Median = '+medistring+' Standard Deviation = '+stdvstring rn=j_yrange(res) jtxt = rn[1]-0.075*(abs(rn[1]-rn[0])) plot, juldates, res, psym=-8, yrange=rn, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Difference between Systolic and Diastolic Pressures (mmHg)',$ thick=2, symsize=ch, /xstyle, /ystyle xyouts, xpos, jtxt, tit, alignment=0.0, charsize=1.0 if plotps eq 1 then begin device, /close set_plot, 'x' endif else stop ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; bp_com ;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if plotps eq 1 then begin a=PSWINDOW(/landscape, /cm, margin=0.025) set_plot, 'ps' outputname= '~/bp/bp_com.ps' device, file=outputname, /landscape endif else begin plot, [0],[1] & wshow endelse !p.multi = [0,2,2] temp = LABEL_DATE(DATE_FORMAT='%N/%D/%Z') ch = 0.750 ;charsize sy = 0.50 ;charsize sysmn = mean(sys) & sysmd = median(sys) & syssd = stddev(sys) meanstring = string(sysmn,format='(f5.1)') medistring = string(sysmd,format='(f5.1)') stdvstring = string(syssd,format='(f4.1)') tit = 'Mean = '+meanstring+' Median = '+medistring+' Standard Deviation = '+stdvstring rn=j_yrange(sys) jtxt = rn[1]-0.075*(abs(rn[1]-rn[0])) plot, juldates, sys, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Systolic Pressure (mmHg)', $ yrange=rn, thick=2, psym=-8, charsize=ch, symsize=sy, /xstyle, /ystyle xyouts, xpos, jtxt, tit, alignment=0.0, charsize=ch dysmn = mean(dys) & dysmd = median(dys) & dyssd = stddev(dys) meanstring = string(dysmn,format='(f4.1)') medistring = string(dysmd,format='(f4.1)') stdvstring = string(dyssd,format='(f3.1)') tit = 'Mean = '+meanstring+' Median = '+medistring+' Standard Deviation = '+stdvstring rn=j_yrange(dys) jtxt = rn[1]-0.075*(abs(rn[1]-rn[0])) plot, juldates, dys, psym=-8, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Diastolic Pressure (mmHg)',yrange=rn, $ thick=2, charsize=ch, symsize=sy, /xstyle, /ystyle xyouts, xpos, jtxt, tit, alignment=0.0, charsize=ch hrmn = mean(hrt) & hrmd = median(hrt) & hrsd = stddev(hrt) meanstring = string(hrmn,format='(f4.1)') medistring = string(hrmd,format='(f4.1)') stdvstring = string(hrsd ,format='(f4.1)') tit = 'Mean = '+meanstring+' Median = '+medistring+' Standard Deviation = '+stdvstring rn=j_yrange(hrt) jtxt = rn[1]-0.075*(abs(rn[1]-rn[0])) plot, juldates, hrt, psym=-8, yrange=rn, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Heart Rate (beats /min)', thick=2, charsize=ch, symsize=sy, /xstyle, /ystyle xyouts, xpos, jtxt, tit, alignment=0.0, charsize=ch resmn = mean(res) & resmd = median(res) & ressd = stddev(res) meanstring = string(resmn,format='(f4.1)') medistring = string(resmd,format='(f4.1)') stdvstring = string(ressd ,format='(f3.1)') tit = 'Mean = '+meanstring+' Median = '+medistring+' Standard Deviation = '+stdvstring rn=j_yrange(res) jtxt = rn[1]-0.075*(abs(rn[1]-rn[0])) plot, juldates, res, psym=-8, yrange=rn, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Difference between Systolic and Diastolic Pressures (mmHg)',$ thick=2, charsize=ch, symsize=sy, ycharsize=0.7, /xstyle, /ystyle xyouts, xpos, jtxt, tit, alignment=0.0, charsize=ch !p.multi = 0 if plotps eq 1 then begin device, /close set_plot, 'x' endif else stop if plotbmp eq 1 then begin wshow write_bmp, 'bp.bmp', tvrd() endif ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; bp_rmn ;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if plotps eq 1 then begin a=PSWINDOW(/landscape, /cm, margin=0.025) set_plot, 'ps' outputname= '~/bp/bp_rmn.ps' device, file=outputname, /landscape endif else begin plot, [0],[1] & wshow endelse !p.multi = [0,2,2] temp = LABEL_DATE(DATE_FORMAT='%N/%D/%Z') ch = 0.750 ;charsize sy = 0.50 ;charsize rmnsmn = mean(rmns) & rmnsmd = median(rmns) & rmnssd = stddev(rmns) meanstring = string(rmnsmn,format='(f5.1)') medistring = string(rmnsmd,format='(f5.1)') stdvstring = string(rmnssd,format='(f4.1)') tit = 'Mean = '+meanstring+' Median = '+medistring+' Standard Deviation = '+stdvstring rn=j_yrange(rmns,buff=[0.60,0.60]) jtxt = rn[1]-0.075*(abs(rn[1]-rn[0])) plot, juldates, rmns, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Systolic Pressure (mmHg)', $ yrange=rn, thick=2, psym=-8, charsize=ch, symsize=sy, /xstyle, /ystyle xyouts, xpos, jtxt, tit, alignment=0.0, charsize=ch rmndmn = mean(rmnd) & rmndmd = median(rmnd) & rmndsd = stddev(rmnd) meanstring = string(rmndmn,format='(f4.1)') medistring = string(rmndmd,format='(f4.1)') stdvstring = string(rmndsd,format='(f3.1)') tit = 'Mean = '+meanstring+' Median = '+medistring+' Standard Deviation = '+stdvstring rn=j_yrange(rmnd,buff=[0.60,0.60]) jtxt = rn[1]-0.075*(abs(rn[1]-rn[0])) plot, juldates, rmnd, psym=-8, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Diastolic Pressure (mmHg)',yrange=rn, $ thick=2, charsize=ch, symsize=sy, /xstyle, /ystyle xyouts, xpos, jtxt, tit, alignment=0.0, charsize=ch rmnhmn = mean(rmnh) & rmnhmd = median(rmnh) & rmnhsd = stddev(rmnh) meanstring = string(rmnhmn,format='(f4.1)') medistring = string(rmnhmd,format='(f4.1)') stdvstring = string(rmnhsd ,format='(f4.1)') tit = 'Mean = '+meanstring+' Median = '+medistring+' Standard Deviation = '+stdvstring rn=j_yrange(rmnh,buff=[0.60,0.60]) jtxt = rn[1]-0.075*(abs(rn[1]-rn[0])) plot, juldates, rmnh, psym=-8, yrange=rn, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Heart Rate (beats /min)', thick=2, charsize=ch, symsize=sy, /xstyle, /ystyle xyouts, xpos, jtxt, tit, alignment=0.0, charsize=ch rmnrmn = mean(rmnr) & rmnrmd = median(rmnr) & rmnrsd = stddev(rmnr) meanstring = string(rmnrmn,format='(f4.1)') medistring = string(rmnrmd,format='(f4.1)') stdvstring = string(rmnrsd ,format='(f3.1)') tit = 'Mean = '+meanstring+' Median = '+medistring+' Standard Deviation = '+stdvstring rn=j_yrange(rmnr,buff=[0.60,0.60]) jtxt = rn[1]-0.075*(abs(rn[1]-rn[0])) plot, juldates, rmnr, psym=-8, yrange=rn, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Difference between Systolic and Diastolic Pressures (mmHg)',$ thick=2, charsize=ch, symsize=sy, ycharsize=0.7, /xstyle, /ystyle xyouts, xpos, jtxt, tit, alignment=0.0, charsize=ch !p.multi = 0 if plotps eq 1 then begin device, /close set_plot, 'x' endif else stop ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; bp_csm ;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if plotps eq 1 then begin a=PSWINDOW(/landscape, /cm, margin=0.025) set_plot, 'ps' outputname= '~/bp/bp_csm.ps' device, file=outputname, /landscape endif else begin plot, [0],[1] & wshow endelse temp = LABEL_DATE(DATE_FORMAT='%N/%D/%Z') ch = 0.750 ;charsize sy = 0.50 ;charsize ;rn1 = [min(csums)-20.,max(csums)+20.] ;rn2 = [min(csumd)-20.,max(csumd)+20.] ;rn3 = [min(csumh)-20.,max(csumh)+20.] ;rn4 = [min(csumr)-20.,max(csumr)+20.] !p.multi = [0,2,2] plot, juldates, csums, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Systolic Pressure CUSUM (mmHg)', $ thick=2, psym=-8, charsize=ch, symsize=sy, $ yrange=j_yrange(csums,buff=[0.30,0.20]), /ystyle oplot, [juldates[0]-500,juldates[n-1]+500], [0.,0.], lines=2 legend, ['All plots are a mean CUSUM on raw data'], /bottom, $ charsize=ch*.8/.75 plot, juldates, csumd, psym=-8, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Diastolic Pressure CUSUM (mmHg)', $ thick=2, charsize=ch, symsize=sy, $ yrange=j_yrange(csumd), /ystyle oplot, [juldates[0]-500,juldates[n-1]+500], [0.,0.], lines=2 plot, juldates, csumh, psym=-8, xtickformat='LABEL_DATE', $ xtitle='Date', ytitle='Heart Rate CUSUM (beats /min)', thick=2, $ charsize=ch, symsize=sy, $ yrange=j_yrange(csumh), /ystyle oplot, [juldates[0]-500,juldates[n-1]+500], [0.,0.], lines=2 ;legend, ['All plots are a mean CUSUM on running means'], $ ; charsize=ch*.8/.75 plot, juldates, csumr, psym=-8, xtickformat='LABEL_DATE', $ xtitle='Date', $ ytitle='Difference between Systolic and Diastolic Pressures CUSUM (mmHg)',$ thick=2, charsize=ch, symsize=sy, ycharsize=ch*0.7/.75, $ yrange=j_yrange(csumr), /ystyle oplot, [juldates[0]-500,juldates[n-1]+500], [0.,0.], lines=2 !p.multi = 0 if plotps eq 1 then begin device, /close set_plot, 'x' endif else stop if plotbmp eq 1 then begin wshow write_bmp, 'bp.bmp', tvrd() endif @bptime stop end