pro measles_model ; Program for numerical solution of PDE of measles transmission ; Andrew Noymer & Katelyn Corey ; Sept 2014-May 2016, Irvine ; USE FONT: Andale Mono(space) 38pt ; "tic" and "toc" are timer calls in IDL ; uncomment this "tic" and the "toc" below ; to get run-times for the program ; ; TIC Omega=50 ; the oldest age to which anyone lives ipy=146 ; iterations per year a/k/a 2.5 day time/age step ypi=double(1.0/(1.0*ipy)) ; years per iteration (floating-point) num_steps=Omega*ipy ; 7300 if ipy and Omega remain the same dps=double(365.0/ipy) ; "days per step" ; the following is /12 to convert months to years... ; median age (post-expiration of maternal Ab... see notes ; "lambda.pdf"... six months is subtracted [already]). This is ; drawn from the literature and may vary from place to place etc. median_age_seroconvert=ipy*(24.0/12.0) two=double(2.0) ; double-precision 2.0 mean_age_sero=median_age_seroconvert/ALOG(two) mean_contagious=7.0 ; in days cdr = double(1.0/mean_contagious) ; "contagious decay rate" ncp = EXP(-cdr*dps) ; "net contagious persistence" mean_latent=10.0 ; in days ldr=double(1.0/mean_latent) ; "latent decay rate" nlp=EXP(-ldr*dps) ; "net latent persistence" ; population classes: general=dblarr(num_steps) ; population before apportionment ; "general" is just a placeholder ; for the population while the model ; does some set-up work. one can think ; of it as the "total population" but it ; is un-used once the setup work is done lifetable=dblarr(num_steps) ; life table l(x) column M=dblarr(num_steps) ; maternal S=dblarr(num_steps) ; susceptible L=dblarr(num_steps) ; latent C=dblarr(num_steps) ; contagious Z=dblarr(num_steps) ; immune ; rates: mu=dblarr(num_steps) ; read data: ; mu is a file containing survival probabilities ; produced separately: calc_lifetable_v16.pro get_lun, LUN openr, LUN, 'mu.txt' readf, LUN, mu close, LUN free_lun, LUN ;; Run quick simulation to convert mu[·] into ;; a lifetable l(x) column (which will be needed ;; to pre-populate a stable population, but which ;; is not used in the model per se) lifetable[0]=1.0 for time=1,num_steps-1 do begin lifetable[time]=$ lifetable[time-1]*mu[time-1] endfor ;; ;; start by populating the model ;; TotalN=ulong(300000) ; starting population of model district PropUnderFifty=0.75 ; from Coale and Demeny (e.g. p.183) N=TotalN*PropUnderFifty ; de facto population growth=double(2.6/100.0); population growth rate in percent (DHS estimate) CBR=47.0/1000.0 ; crude birth rate (DHS country report, p.4 [PDF p.29], "Taux de natalité") ;; HEATMAP - comment-out this section and its associated loop ;; we are tagging the heatmap code with ";%;" -- these are the ;; lines that need to be un-commented to do the heatmap. ;%; num_growth_rates=100 ;%; min_growth_rate=double(0.5/100.0) ;%; max_growth_rate=double(2.8/100.0) ;%; growth_rate_range=max_growth_rate-min_growth_rate ;%; growth_rate_step=growth_rate_range/float(num_growth_rates) ;%; growth_rate=min_growth_rate+$ ;%; growth_rate_step*findgen(num_growth_rates) ;%; num_immuniz_rates=100 ;%; min_immuniz_rate=double(0.60) ;%; max_immuniz_rate=double(1.00) ;%; immuniz_rate_range=max_immuniz_rate-min_immuniz_rate ;%; immuniz_rate_step=immuniz_rate_range/float(num_immuniz_rates) ;%; immuniz_rate=min_immuniz_rate+$ ;%; immuniz_rate_step*findgen(num_immuniz_rates) ; heatmap is a 2D-array of waiting TIMES so is an integer-type array ;%; heatmap=UINTARR(num_growth_rates,num_immuniz_rates) ;%; imm_map=DBLARR(num_growth_rates,num_immuniz_rates) ;%; r_map=DBLARR(num_growth_rates,num_immuniz_rates) ;%; FOR growth_loop=0,num_growth_rates-1 DO BEGIN ;%; growth=growth_rate[growth_loop] ;; POPULATION INITIALIZATION ;; step 1: populate the "general" array ;; with a stable population structure based ;; off a calcuation, which should be v. close ;; to the correct stable pop structure. This ;; saves a ton of execution time vs. starting ;; with arbitrary (e.g. uniform) pop structure ;; and simulating-out the equilibrium ;; ;; ;; calibrate the CURRENT birth cohort size ;; based on the empirical crude birth rate ;; (recall that ypi=1/146<<1) starting_birth_cohort=TotalN*CBR*ypi age=Omega*(findgen(num_steps)/num_steps) ; see Keyfitz, AMD (Springer Edn 1985), p.78 ; note that "dx" is not included below but ; is implicit because "ypi"==dx and is ; in the above defn of starting_birth_cohort general=starting_birth_cohort*$ EXP(-growth*age)*lifetable ;; ;; then apportion into M,S,L,C,Z classes ;; ;; first part: rough ; scale pop: ; this just restores the pop to N=TotalN*PropUnderFifty ; because the stable-pop calculation inflated it a bit ; this does not change the structure, just the total totalpop=total(general,/double) paf=(N/totalpop) ; population adjustment factor general=paf*general ; this variable just control the number of steps in which ; to sprinkle seed cases steps_latent=1 ;; pre-assignment of classes ; assign M zeta=(ipy/2)-1 M[0:zeta]=general[0:zeta] ; assign S ; kappa (below) is a fudge-factor to adjust lambda so that ; (BY TRIAL AND ERROR) the "model-empirical" age at 50% ; seropositivity equals the true-empirical value indicated ; above see repeat/until loop below for final calculation ; ; kappa = 1.0 is what one expects from theory but that ; assumes a rectangular population whereas our model uses ; a Lotka population, i.e. demographic stability. In other ; words, there are design effects. So we make a small ; adjustment to lambda to get out from the simulation the ; correct median age of serology (etc.) ; ; kappa is a design-effect paramter and should not deviate ; from 1.0 by more than a few percent. we align kappa by ; trial and error but it could also be automated such that ; everything below is a repeat..until loop. It didn't take ; a lot of iteration to assign kappa and it's easier to do ; a few iterations and then lock-in the value. in any case, ; if we are iterating we would search in a neighborhood of kappa=1.0 kappa = 1.03 lambda=kappa*1.0/mean_age_sero age_dummy=findgen(num_steps-zeta) S[zeta+1:*]=$ general[zeta+1:*]*(exp(-lambda*age_dummy)) ;*ypi ; assign L ; L&C don't matter much since they don't affect the force ; of transmission which is fixed (for now only...) at lambda ; we can sprinkle a number of L&C into the population just ; for fun, but in the end the model has to run until a stable ; equilibrum is achieved... ; ; number of seeded cases to be inserted: seeds=4.0 ; make sure it doesn't exceed the population: num_seeds=min([seeds,general[zeta+1]]) ; insert that many latents: L[zeta+1]=num_seeds ; deduct the seeds from S: S[zeta+1]=S[zeta+1]-num_seeds ; assign Z Z=general-C-L-S-M ;; NOW we run a simulation (until equilibrium) OF measles ;; transmission with force of infection fixed at lambda ; run the model for 50 years (=Omega) to achieve a full ; demographic/epidemiologic equilibrium ...trial and error ; shows that it takes this long to achieve a full equilibrium max_init_inteations=num_steps ; pre-calculate certain quantities to save execution time: exp_lambda=exp(-lambda) for time=0,max_init_inteations-1 do begin S_old=S L_old=L C_old=C ; note that the left hand side of the array assignments below ; are single array cells, but what is being put there is a ; multi-element array subset. IDL "knows" how to do this and ; makes sure that the overflow (so to say) goes to the right place ; ; for more details, see this URL: ; http://www.harrisgeospatial.com/Company/PressRoom/Blogs/IDLDataPointDetail/TabId/902/ArtMID/2926/ArticleID/14505/Increasing-Performance-when-Assigning-Subsets-of-Arrays.aspx ; if this link is broken, see array_assign.pdf ; why not just do S[zeta+1:*] instead of S[zeta+1] (e.g. in ; the LHS, first line of code below)? it is clearer code. ; A: yes, it is clearer code, but due to IDL internals (described ; at the above URL), it executes slower. We are making a small ; decrease in code readability in favor of speed. And given that ; IDL is designed to work this way, the decrease in readability ; is debatable (i.e., if the reader knows IDL and is not just ; treating it all as pseudocode, then it's clear enough what's ; going on). How much is the speed increase? About 4% in ; our tests, so not large but certainly not negligible. ; And, are we sure it works? Yes, tested. S[zeta+1]=$ mu[zeta:-2]*S_old[zeta:-2]*exp_lambda L[zeta+1]=$ mu[zeta:-2]*(L_old[zeta:-2]*nlp+S_old[zeta:-2]*(1.0-exp_lambda)) C[zeta+1]=$ mu[zeta:-2]*((C_old[zeta:-2]*ncp)+L_old[zeta:-2]*(1.0-nlp)) Z[zeta+1]=$ mu[zeta:-2]*(Z[zeta:-2]+C_old[zeta:-2]*(1.0-ncp)) ; The value of S, just above zeta, is now zero, since it was ; "aged" from s[zeta] in the above loop. S[zeta+1] in fact comes ; from M[zeta]. i.e. when M[zeta] age, they age into S[zeta+1]. ; The next line does this. ; transition from M to S: S[zeta+1]=mu[zeta]*M[zeta] M[1:zeta]=mu[0:zeta-1]*(M[0:zeta-1]) ; birth cohort function: M[0]=M[0]*EXP(growth*ypi) endfor totalpop=total(M+S+L+C+Z,/double) beta = lambda * totalpop/total(C,/double) ; print, 'beta', beta ; re-apportion the population (again) paf=(N/totalpop) ; population adjustment factor M=M*paf S=S*paf L=L*paf C=C*paf Z=Z*paf G=M+S+L+C+Z ; save initial configuration of starting population... ; we will need it later. (This is only needed for the ; heatmap operations, or similar.) ;%; M_orig=M ;%; S_orig=S ;%; L_orig=L ;%; C_orig=C ;%; Z_orig=Z ;%; G_orig=G ; median serology: ; locate youngest age at which >=50% of ; pop has measles Ab (i.e. is recovered) ; ; THIS LOOP MAY BE COMMENTED ONCE kappa is set ;; x=zeta-1 ;; REPEAT BEGIN ;; ++x ;; ENDREP UNTIL Z[x]/G[x] GE 0.5 ;; IMMUNIZATION LOOP ;%; FOR IMM_count=0,num_immuniz_rates-1 DO BEGIN ;; ;; vaccination issues: ;; max_model_iterations=25*ipy ; 25 years ; 0-2 years (of time): no vaccination, at all ; (but model is running in its perturbed state) ; ; at 2 years (of time: vaccination starts, ; with a campaign. as described below. ; ; starting immed. after the campaign, there is ; programmatic vaccination, as described below ; ; (i) evolution of vaccination coverage over time. vaccination ; coverage was 30% for the first two years, then 35% for two years. ; Then 55% for two years. Then maxes out at 65% and remains. ; program vaccine efficacy is 75% ; ; (ii) at 2 years into the model, there is a campaign. ; the campaign looks like this: ; vaccinate kids, 9 to 23 months of age. ; coverage is 70% of pop ; vaccine efficacy is 80% ; ; In this version of the model, there are more frequent campaigns. ; It's governed by the variable camp_time, which is an array with ; entries at each time step, 1 if campaign, 0 otherwise. See below. ; ; A note about implementation of vaccine camaign: the most ; runtime-efficient way to do this would be to run the model ; for 2 years, then simulate the vaccine campaign, then run ; the loop some more. this would require two loops... which is not ; runtime-costly, but adds approx 50 lines of code (incl blank ; lines + comments), which breaks up the flow of the program from ; a human-reading perspective. what we do instead is add a boolean ; eval inside the loop, checking for whether it's campaign-time or ; not. the disadvantage of this is that this boolean eval will ; be executed in every loop iteration, so is costly for runtime, ; but shouldn't be too costly. at this point we prefer to have ; readable code. ; ; the PROGRAM vaccination is implemented as an vector (over time) ; of the immunization rate, with the pre-vax era set to 0.0. HOWEVER, ; we cannot implement campaign vaccination this way since campaign ; vaccination occurs over an age RANGE, not a single-point of age, ; so the analogue would be a matrix (age*time) not a vector, ; and that would be more costly than a boolean eval. ; net immunization is: coverage * efficacy. ; for program vaccination we make an array of net immunization X time: immuniz=dblarr(max_model_iterations) vax_begins=2*ipy ; is there a campaign now? camp_time=bytarr(max_model_iterations) ; 2 years vax_campaign_time = vax_begins camp_time[vax_campaign_time]=1 campaign_frequency=1.7 ; years c_f_steps=FLOOR(ipy*campaign_frequency) FOR vct = vax_campaign_time+1,max_model_iterations-1 do begin if (vct MOD c_f_steps) EQ 0 then camp_time[vct]=1 endfor ;print, 'number of campaigns: ', total(camp_time,/integer) ; vaccine coverage: v1=0.3 v2=v1 v3=0.35 v4=v3 v5=0.55 v6=v5 v7=0.65 v8=0.65 ;v8=0.85 ;v8=immuniz_rate[IMM_count] ; remove or comment this block to return to BAU: ; the following is a jump right into ; the final vaccination coverage ;; v1=v8 ;; v2=v8 ;; v3=v8 ;; v4=v8 ;; v5=v8 ;; v6=v8 ;; v7=v8 ; "simulated Muyinga" was only 65% coverage ; with a 75% effective vaccine vax_efficacy=0.75 ;vax_efficacy=1.00 ;;; this represents PROGRAM vax only, not campaign immuniz[0:vax_begins-1]=0.0 immuniz[vax_begins]=0.0 ; campaign vax not program immuniz[vax_begins+1:3*ipy-1]=v1*vax_efficacy immuniz[3*ipy:4*ipy-1]=v2*vax_efficacy immuniz[4*ipy:5*ipy-1]=v3*vax_efficacy immuniz[5*ipy:6*ipy-1]=v4*vax_efficacy immuniz[6*ipy:7*ipy-1]=v5*vax_efficacy immuniz[7*ipy:8*ipy-1]=v6*vax_efficacy immuniz[8*ipy:9*ipy-1]=v7*vax_efficacy immuniz[9*ipy:*]=v8*vax_efficacy age_vax=CEIL((9.0/12.0)*ipy) ; 9 months to 2 years of age... ;; note that we are changing this now... some of the ;; above comments may indicate that camp. vax is 9mo-2yr ;; but in reality it's whatever is immed. below: max_age_campaign_vax=(6*ipy)-1 ; ; since age_vax is used for program vaccination as well, ; define the following: min_age_campaign_vax=(2*ipy)-1 ;; DATE of vaccine campaign ; age-structure of campaign vaccination: campaign_vax_by_age=dblarr(num_steps) campaign_coverage=double(.7) campaign_efficacy=double(.8) net_immun=campaign_coverage*campaign_efficacy number_steps_campaign_vax=$ 1+max_age_campaign_vax-min_age_campaign_vax imminz_profile=replicate(net_immun,number_steps_campaign_vax) dummy_vec_0=dblarr(min_age_campaign_vax-1) dummy_vec_1=$ dblarr(num_steps-((min_age_campaign_vax-1)+number_steps_campaign_vax)) success_campaign_vax= $ [ dummy_vec_0, imminz_profile, dummy_vec_1] ; so, we have an array, num_steps in length, which ; gives the net immunizing effect (by age) of the ; CAMPAIGN vaccination. It is mostly zero, but ; from 9 to 23 months it is "net_immun" (see definition ; above of that constant. Multiply the susceptibles by ; this (element-wise) and we get the net immunized. Then subtract ; these from the susceptibles and add them to the immunized, ; and that's the campaign, done. One array multiplication ; (element-wise), and two array additions/subtractions ;; ;; final model: ;; agevector_detail=findgen(num_steps)/ipy ; "R" for reporting statistics R_mean=DBLARR(max_model_iterations) R_SD=DBLARR(max_model_iterations) measles_over_time=dblarr(max_model_iterations) measles_surface=dblarr(num_steps,max_model_iterations) print, 'perturbing the model' perturbation=.1*S S=S-perturbation Z=Z+perturbation age_x_section=DBLARR(num_steps,max_model_iterations/ipy) age_x_counter=0 timestamp=LONARR(max_model_iterations/ipy) ; cbi: cumulative (time-)binned incidence cbi=dblarr(1+CEIL(max_model_iterations/12)) cbi_count=LONG(-1) ; data structure to record the average of incident ; cases from an arbitrary start to an arbitrary ; stop point (e.g. a trough-to-trough epidemic) ; start+end are determined by trial and error ; (see "diagnostic" code in graph) epid_cases=dblarr(num_steps) epid_start=455 epid_end=epid_start+480 ;; variable to keep track of ascending vs descending epidemic epi_peak=double(0.0) old_epi_peak=epi_peak time_peak=long(0) ;;; BEGIN MAIN LOOP ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; for time=0,max_model_iterations-1 do begin ; start of every iteration: ; do PROGRAM vaccination: immunized=S[age_vax]*immuniz[time] S[age_vax]=S[age_vax]-immunized Z[age_vax]=Z[age_vax]+immunized ; now we implement CAMPAIGN vaccination: if camp_time[time] EQ 1 then begin camp_immunized=S*success_campaign_vax S=S-camp_immunized Z=Z+camp_immunized endif totC=total(C,/double) measles_over_time[time]=totC measles_surface[*,time]=C ; this loop is basically like that above but ; lambda is now a function that changes each ; iteration lambda=beta*totC/total(M+S+L+C+Z) ; to save execution time: exp_lambda=exp(-lambda) S_old=S L_old=L C_old=C S[zeta+1:*]=$ mu[zeta:-2]*S_old[zeta:-2]*exp_lambda L[zeta+1:*]=$ mu[zeta:-2]*(L_old[zeta:-2]*nlp+S_old[zeta:-2]*(1.0-exp_lambda)) C[zeta+1:*]=$ mu[zeta:-2]*((C_old[zeta:-2]*ncp)+L_old[zeta:-2]*(1.0-nlp)) new_cases=$ total(mu[zeta:-2]*(L_old[zeta:-2]*(1.0-nlp)),/double) Z[zeta+1:*]=$ mu[zeta:-2]*(Z[zeta:-2]+C_old[zeta:-2]*(1.0-ncp)) ; transition from M to S: S[zeta+1]=mu[zeta]*M[zeta] M[1:zeta]=mu[0:zeta-1]*(M[0:zeta-1]) ; birth cohort function: M[0]=M[0]*EXP(growth*ypi) ;; REPORTING: ; mean age is weighted average of age and C vecotr R_mean[time]=total(C*agevector_detail,/double)/total(C,/double) ; geomean & harmonic mean are not informative in add to arithmetic ;R_geomean[time]=EXP(total(C[1:*]*ALOG(agevector_detail[1:*]))/total(C[1:*])) ; ;R_harmean[time]=total(C[1:*])/total(C[1:*]/agevector_detail[1:*]) R_SD[time]=$ SQRT(TOTAL(C*(agevector_detail-R_mean[time])^2,/double)/total(C,/double)) if (time MOD ipy) EQ 0 THEN BEGIN timestamp[age_x_counter]=time age_x_section[*,age_x_counter++]=C ENDIF ; monthly binning so 12 time steps makes 30 days IF (time mod 12) EQ 0 THEN cbi_count++ cbi[cbi_count]=cbi[cbi_count]+new_cases ; some small gains in runtime can be achieved ; by commenting-out the following, when the ; age distribution data (as in "table 3" of the original ; paper) are not needed IF (time LT epid_end) AND (time GE epid_start) THEN BEGIN epid_cases=epid_cases+(mu*(L_old*(1.0-nlp))) ENDIF ; avoid roundoff problems at sustained low incidence: ; 0.0000150 produces 5 bad 'pixels' at 100x100 ; 0.00000750 produces 2 bad 'pixels' at 100x100 ;%; tolerance=double(0.00000100) ;; ***************************************************** ;; comment-out this to turn OFF short-circuting the model ;; looking for triangle pattern ;%; ;; IF (time GT vax_campaign_time) AND $ ;%; ;; (new_cases LT epi_peak) AND $ ;%; ;; ((epi_peak-tolerance) GT old_epi_peak) THEN BEGIN ;%; ;; time_peak=time-1 ;%; ;; GOTO, LOOP_END ;%; ;; ENDIF old_epi_peak=epi_peak epi_peak=new_cases time_peak=time endfor ;%; LOOP_END: ;print, 'peak ', time_peak ;; heatmap[growth_loop,imm_count]=time_peak ;; imm_map[growth_loop,imm_count]=v8 ;; r_map[growth_loop,imm_count]=growth ; we are now at the bottom of the ; IMMUNIZ loop. Before we start again, ; reset M,S,L,... to starting value ;%; M=M_orig ;%; S=S_orig ;%; L=L_orig ;%; C=C_orig ;%; Z=Z_orig ;%; G=G_orig ;%; ;ENDFOR ; IMM_count ;%; ;ENDFOR ; growth_loop ;%; openw, lun, 'heatmap.txt', /get_lun, width=10*num_growth_rates ;%; printf, lun, heatmap ;%; free_LUN, lun ;; print, heatmap ;; print ;; print, imm_map ;; print ;; print, r_map ; calculate ages for "table 3" (rel. to orig paper) ; age groups (months) were as follows: ;; 0-5 ;; 6-11 ;; 12-23 ;; 24-35 ;; 36-59 ;; 60-600 table_3=dblarr(6) ; bin the cases by age: table_3[0]=total(epid_cases[0:72],/double) table_3[1]=total(epid_cases[73:145],/double) table_3[2]=total(epid_cases[146:291],/double) table_3[3]=total(epid_cases[292:437],/double) table_3[4]=total(epid_cases[438:583],/double) table_3[5]=total(epid_cases[584:*],/double) ; make it proportional: table_3=table_3/total(epid_cases,/double) ;; print, 'table 3' ;; print, transpose(table_3) ;; alternate version of table 3 where we ignore ;; cases >5 years (like the paper) table_3[0]=total(epid_cases[0:72],/double) table_3[1]=total(epid_cases[73:145],/double) table_3[2]=total(epid_cases[146:291],/double) table_3[3]=total(epid_cases[292:437],/double) table_3[4]=total(epid_cases[438:583],/double) ; make it proportional: table_3=table_3/total(epid_cases[0:583],/double) print, 'table 3' print, transpose(table_3[0:4]) ; check growth rate: ; print, ALOG(POP_ONE/POP_ZERO)/10.0 ;;;;;;; downsample the result ;;;;;;;;;;;;;;;;; ; ; for the wireframe surface we want ; less than the full 25 years of ; data: surf_subset=min([10*ipy,max_model_iterations]) ;years new_cols=num_steps/25 new_rows=surf_subset/25 measles_surface_downsamp=$ congrid(measles_surface[*,0:surf_subset-1],new_cols,new_rows) ; REBIN == congrid and is is faster but the ; new dimensions must be an integer multple or ; factor of the old dimensions ; drop over 25 age25=floor(25*ipy/25) measles_surface_downsamp3=measles_surface_downsamp[0:age25,*] ;;;;;;; make graphs ;;;;;;;;;;;;;;;;; ; make graph ;agevector_complete=indgen(51) ; print, agevector_detail[x] SET_PLOT, 'PS' ; use PostScript fonts: !P.FONT=0 ; ; thickness (for line weight, etc, 1.0 is IDL default): CH=4.5 CS=2 ;CS=3.5 T=7 U=T ;U=T+2 ; LOADCT, 0 color_Measles=x11colors('dark_olive_green') shade=x11colors('dark_grey');x11colors('gainsboro');x11colors('light_grey');x11colors('gainsboro') ; make plot: ;%; DEVICE, XSize=9, YSize=9, /INCHES, /ENCAPSULATED, /DECOMPOSED, $ ;%; /COLOR, Bits_Per_Pixel=8, File='heat_map.eps' ;%; !P.MULTI = [0, 1, 2] ;%; ymax=num_immuniz_rates ;%; ymin=0 ;%; xmin=0 ;%; xmax=num_growth_rates ;%; num_x_ticks=6 ;%; xtickmarks=100.0*(min_growth_rate+$ ;%; findgen(num_x_ticks)*(max_growth_rate-min_growth_rate)/float(num_x_ticks-1)) ;%; num_y_ticks=3 ;%; ytickmarks=100.0*(min_immuniz_rate+$ ;%; findgen(num_y_ticks)*(max_immuniz_rate-min_immuniz_rate)/float(num_y_ticks-1)) ;%; PLOT, findgen(num_growth_rates), findgen(num_immuniz_rates), YRANGE=[ymin,ymax], $ ;%; XRANGE=[xmin,xmax], ticklen=-0.02, $ ;%; /NODATA, $ ;%; XTHICK=T, YTHICK=T, YSTYLE=1, XSTYLE=1, $ ;%; XMINOR=1, YMINOR=2, YTICKS=num_y_ticks-1, XTICKS=num_x_ticks-1, $ ;%; XTITLE='Population growth rate, r (%)', $ ;%; YTITLE='Immunization rate (%)', $ ;%; ytickname=string(ytickmarks,format='(F5.1)'), $ ;%; xtickname=string(xtickmarks,format='(F4.2)'), $ ;%; CHARTHICK=T, $ ;%; CHARSIZE=CS, $ ;%; POSITION=[0.2,0.3,0.95,0.9] ;%; print, 'fastest epidemic: ', min(heatmap)/float(ipy) ;%; print, 'longest time till (no) epidemic: ', max(heatmap)/float(ipy) ; Flip the color because highest value means longest ; wait between epidemics, which means "coolest" epi ; sitation. But naturally, higher values are coded ; red, or "hottest". ;%; heat_map=byte(255)-bytscl(heatmap) ;%; for col=0,num_growth_rates-1 do begin ;%; for row=0,num_immuniz_rates-1 do begin ;%; polyfill, [col,col,col+.99999,col+.99999,col],$ ;%; [row,row+.99999,row+.99999,row,row], $ ;%; color=blue_red_256(heat_map[col,row]) ;%; ; now foolproof way to know what is "white: " ;%; IF heat_map[col,row] EQ 127 THEN BEGIN ;%; print, 'white: ', heatmap[col,row] ;%; ENDIF ;%; endfor ;%; endfor ; re-draw axes where POLYFILL has rubbed up against ;%;AXIS, XAX=0, XRANGE=[XMIN,XMAX], Xstyle=1, Xticks=1, XTHICK=T, $ ;%; XTICKNAME=replicate(' ',2), XTICKLEN=-0.0001, /DATA ;%; ;%;AXIS, XAX=1, XRANGE=[XMIN,XMAX], Xstyle=1, Xticks=1, XTHICK=T, $ ;%; XTICKNAME=replicate(' ',2), XTICKLEN=-0.0001, /DATA ;%; ;%;AXIS, YAX=0, XRANGE=[YMIN,YMAX], Ystyle=1, Yticks=1, YTHICK=T, $ ;%; YTICKNAME=replicate(' ',2), YTICKLEN=-0.0001, /DATA ;%; ;%;AXIS, YAX=1, XRANGE=[YMIN,YMAX], Ystyle=1, Yticks=1, YTHICK=T, $ ;%; YTICKNAME=replicate(' ',2), YTICKLEN=-0.0001, /DATA ; legend: ;%;PLOT, indgen(256), indgen(256), YRANGE=[0,1], $ ;%; XRANGE=[0,255], ticklen=-0.02, $ ;%; /NODATA, YTICKLEN=-0.000001, $ ;%; XTHICK=T, YTHICK=T, YSTYLE=1, XSTYLE=1, $ ;%; XMINOR=1, YMINOR=1, YTICKS=1, XTICKS=2, $ ;%;; XTITLE='Waiting time to epidemic', $ ;%; ytickname=replicate(' ',2), $ ;%; xtickname=['>25 years','15 years','5 years'], $ ;%; CHARTHICK=T, $ ;%; CHARSIZE=1.25, $ ;%; POSITION=[0.19,0.07,0.94,0.12] ;%;x_pos=indgen(256) ;%;legend_colors=bytscl(x_pos) ;%; ;%;for col=0,255 do begin ;%; PLOTS, replicate(x_pos[col],2), $ ;%; [0,1], color=blue_red_256(legend_colors[col]), THICK=T ;%;endfor ;%; ;%;XYOUTS, 0.45, 0.013, 'Waiting time to epidemic', /NORM, $ ;%; charsize=1.25 ;%; ; re-draw axes where PLOTS has rubbed up against ;%;AXIS, XAX=0, XRANGE=[0,255], Xstyle=1, Xticks=1, XTHICK=T, $ ;%; XTICKNAME=replicate(' ',2), XTICKLEN=-0.0001, /DATA ;%; ;%;AXIS, XAX=1, XRANGE=[0,255], Xstyle=1, Xticks=1, XTHICK=T, $ ;%; XTICKNAME=replicate(' ',2), XTICKLEN=-0.0001, /DATA ;%; ;%;AXIS, YAX=0, XRANGE=[0,1], Ystyle=1, Yticks=1, YTHICK=T, $ ;%; YTICKNAME=replicate(' ',2), YTICKLEN=-0.0001, /DATA ;%; ;%;AXIS, YAX=1, XRANGE=[0,1], Ystyle=1, Yticks=1, YTHICK=T, $ ;%; YTICKNAME=replicate(' ',2), YTICKLEN=-0.0001, /DATA ;%; DEVICE, /CLOSE_file ; reset !P.MULTI = 0 ; make plot: DEVICE, XSize=9, YSize=7, /INCHES, /ENCAPSULATED, /DECOMPOSED, $ /COLOR, Bits_Per_Pixel=8, File='timeseries.eps' monthvector=indgen(CEIL(max_model_iterations/12.0)) ymax=3000;ceil(max(cbi)) ymin=30 xmin=0-1.5 xmax=max(monthvector)-1 PLOT, monthvector, cbi, YRANGE=[ymin,ymax], $ XRANGE=[xmin,xmax], ticklen=-0.02, $ /NODATA, /ylog, $ XTHICK=T, YTHICK=T, YSTYLE=1, XSTYLE=1, $ XMINOR=2, YMINOR=9, YTICKS=2, XTICKS=13, $ XTITLE='Time (years)', $ YTITLE='Monthly incidence', $ CHARTHICK=T, $ CHARSIZE=CS, $ xtickv=24*indgen(13), $ xtickname=['0',' ','4',' ','8',' ','12',' ','16',' ','20',' ','24',' '], $ POSITION=[0.2,0.2,0.95,0.9] polyfill, [ XMIN , XMIN, XMAX, XMAX ], [ YMIN, ymax , ymax , YMIN ], color=shade PLOTS, [ XMIn, XMAX ] , [ 300, 300 ], /data, thick=t-2, color=x11colors('white') ;; add a white rule to show vax campaign onset ;; this will be commented-out in many versions ;; of the results PLOTS, [ vax_campaign_time/12.0, vax_campaign_time/12.0 ], $ [ ymin, ymax ], /data, thick=2, color=x11colors('white') ;; add a white rule to confirm peak-timing of 1st post-H epi ;; PLOTS, [ 1155.0/12.0, 1155.0/12.0 ], $ ;; [ ymin, ymax ], /data, thick=2, color=x11colors('blue') ; add epid_start & epid_end (diagnostic only) ;;plots, [epid_start/12.0, epid_start/12.0 ], [ymin, ymax], thick=T ;;plots, [epid_end/12.0, epid_end/12.0 ], [ymin, ymax], thick=T ; plot to the penultimate cbi value because final ; bin may not be 4 full weeks OPLOT, monthvector, cbi[0:-2], THICK=T,$ color=color_measles ; output the data for plotting multiple timeseries on ; one set of axes. timeseries_data=[transpose(monthvector),transpose(cbi)] openw, lun, 'timeseries.txt', /get_lun, width=100 ; reduce scope for re-naming error by storing the ; final vax coverage as the 1st line of the output data printf, lun, vax_efficacy, v8, FORMAT='(F0,"'+string(9b)+'",F0)' printf, lun, timeseries_data[0:1,*], $ FORMAT='(%"%I\t%22.15F")' ; this is tab-separated-file output. the \t is tab and the % causes ; IDL to switch to C-style format codes, which support \t. the ; default format codes do not support \t. For whitespace-delimited, ; just use: FORMAT='(I0," ",F0)' free_LUN, lun ; then manually re-name this file before trying a ; different value for "v8" ; re-draw axes where POLYFILL has rubbed up against AXIS, XAX=0, XRANGE=[XMIN,XMAX], Xstyle=1, Xticks=1, XTHICK=T, $ XTICKNAME=replicate(' ',2), XTICKLEN=-0.0001, /DATA AXIS, XAX=1, XRANGE=[XMIN,XMAX], Xstyle=1, Xticks=1, XTHICK=T, $ XTICKNAME=replicate(' ',2), XTICKLEN=-0.0001, /DATA AXIS, YAX=0, XRANGE=[YMIN,YMAX], Ystyle=1, Yticks=1, YTHICK=T, $ YTICKNAME=replicate(' ',2), YTICKLEN=-0.0001, /DATA AXIS, YAX=1, XRANGE=[YMIN,YMAX], Ystyle=1, Yticks=1, YTHICK=T, $ YTICKNAME=replicate(' ',2), YTICKLEN=-0.0001, /DATA DEVICE, /CLOSE_file ; make plot: DEVICE, XSize=9, YSize=7, /INCHES, /ENCAPSULATED, /DECOMPOSED, $ /COLOR, Bits_Per_Pixel=8, File='timeseries_nolog.eps' monthvector=indgen(CEIL(max_model_iterations/12.0)) ymax=3000;ceil(max(cbi)) ymin=0 xmin=0-1.5 xmax=max(monthvector)-1 PLOT, monthvector, cbi, YRANGE=[ymin,ymax], $ XRANGE=[xmin,xmax], ticklen=-0.02, $ /NODATA, $ XTHICK=T, YTHICK=T, YSTYLE=1, XSTYLE=1, $ XMINOR=2, YMINOR=9, YTICKS=2, XTICKS=13, $ XTITLE='Time (years)', $ YTITLE='Monthly incidence', $ CHARTHICK=T, $ CHARSIZE=CS, $ xtickv=24*indgen(13), $ xtickname=['0',' ','4',' ','8',' ','12',' ','16',' ','20',' ','24',' '], $ POSITION=[0.2,0.2,0.95,0.9] polyfill, [ XMIN , XMIN, XMAX, XMAX ], [ YMIN, ymax , ymax , YMIN ], color=shade PLOTS, [ XMIn, XMAX ] , [ 300, 300 ], /data, thick=t-2, color=x11colors('white') ;; add a white rule to show vax campaign onset ;; this will be commented-out in many versions ;; of the results PLOTS, [ vax_campaign_time/12.0, vax_campaign_time/12.0 ], $ [ ymin, ymax ], /data, thick=2, color=x11colors('white') ;; add a white rule to confirm peak-timing of 1st post-H epi ;; PLOTS, [ 1155.0/12.0, 1155.0/12.0 ], $ ;; [ ymin, ymax ], /data, thick=2, color=x11colors('blue') ; add epid_start & epid_end (diagnostic only) ;;plots, [epid_start/12.0, epid_start/12.0 ], [ymin, ymax], thick=T ;;plots, [epid_end/12.0, epid_end/12.0 ], [ymin, ymax], thick=T ; plot to the penultimate cbi value because final ; bin may not be 4 full weeks OPLOT, monthvector, cbi[0:-2], THICK=T,$ color=color_measles ; output the data for plotting multiple timeseries on ; one set of axes. timeseries_data=[transpose(monthvector),transpose(cbi)] openw, lun, 'timeseries.txt', /get_lun, width=100 ; reduce scope for re-naming error by storing the ; final vax coverage as the 1st line of the output data printf, lun, vax_efficacy, v8, FORMAT='(F0,"'+string(9b)+'",F0)' printf, lun, timeseries_data[0:1,*], $ FORMAT='(%"%I\t%22.15F")' ; this is tab-separated-file output. the \t is tab and the % causes ; IDL to switch to C-style format codes, which support \t. the ; default format codes do not support \t. For whitespace-delimited, ; just use: FORMAT='(I0," ",F0)' free_LUN, lun ; then manually re-name this file before trying a ; different value for "v8" ; re-draw axes where POLYFILL has rubbed up against AXIS, XAX=0, XRANGE=[XMIN,XMAX], Xstyle=1, Xticks=1, XTHICK=T, $ XTICKNAME=replicate(' ',2), XTICKLEN=-0.0001, /DATA AXIS, XAX=1, XRANGE=[XMIN,XMAX], Xstyle=1, Xticks=1, XTHICK=T, $ XTICKNAME=replicate(' ',2), XTICKLEN=-0.0001, /DATA AXIS, YAX=0, XRANGE=[YMIN,YMAX], Ystyle=1, Yticks=1, YTHICK=T, $ YTICKNAME=replicate(' ',2), YTICKLEN=-0.0001, /DATA AXIS, YAX=1, XRANGE=[YMIN,YMAX], Ystyle=1, Yticks=1, YTHICK=T, $ YTICKNAME=replicate(' ',2), YTICKLEN=-0.0001, /DATA DEVICE, /CLOSE_file ymax=CEIL(max(measles_over_time)) xmin=0 xmax=25 sd_color=x11colors('blue') ; cf. measles_v041.pro for the final version that also ; includes an alternate version of the Mean/SD graph ; wherein mean is on one Y-axis and SD is on another. ; ... same thing goes for mean-only + sd-only graphs ; ; make plot: DEVICE, XSize=9, YSize=7, /INCHES, /ENCAPSULATED, /DECOMPOSED, $ /COLOR, Bits_Per_Pixel=8, File='mean_sd.eps' AA=xmin-0.4 BB=xmax+0.4 CC=2.3 DD=5.1 ygridlines=[3,4,5] xgridlines=[0,10,20] PLOT, agevector_detail, R_SD, YRANGE=[CC,DD], $ XRANGE=[AA,BB], ticklen=-0.02, $ /NODATA, $ XTHICK=T, YTHICK=T, YSTYLE=1, XSTYLE=1, $ XMINOR=2, YMINOR=2, YTICKS=2, XTICKS=5, $ XTITLE='time (years)', $ YTITLE='age of infected (years)', $ CHARTHICK=T, $ CHARSIZE=CS, $ ytickv=ygridlines,$ xtickv=[0,5,10,15,20,25] POSITION=[0.2,0.2,0.9,0.9] polyfill, [ AA, AA, BB, BB, AA ], $ [ CC, DD, DD, CC, CC ], color=x11colors('dark_grey') for i=0,N_elements(ygridlines)-1 do begin plots, [ AA, BB ], $ replicate(ygridlines[i], 2), color=x11colors('white'), THICK=T plots, replicate(xgridlines[i], 2), $ [ CC, DD], color=x11colors('white'), THICK=T endfor OPLOT, agevector_detail[0:max_model_iterations-1], $ R_SD, THICK=T,$ color=sd_color OPLOT, agevector_detail[0:max_model_iterations-1], $ R_mean, THICK=T XYOUTS, 20.2,4.25,'SD',/data,color=SD_color, charsize=2.0, charthick=1.5 XYOUTS, 15.8, 4.6, 'Mean', /data, charsize=2.0, charthick=1.5 AXIS, YAX=0, YTHICK=T, YSTYLE=1, $ YMINOR=1, YTICKS=1, YTICKLEN=-0.0000000005, $ YTICKNAME=[' ',' '] AXIS, YAX=1, YTHICK=T, YSTYLE=1, $ YMINOR=1, YTICKS=1, YTICKLEN=-0.0000000005, $ YTICKNAME=[' ',' '] AXIS, XAX=1, XTHICK=T, XSTYLE=1, $ XMINOR=1, XTICKS=1, XTICKLEN=-0.0000000005, $ XTICKNAME=[' ',' '] AXIS, XAX=0, XTHICK=T, XSTYLE=1, $ XMINOR=1, XTICKS=1, XTICKLEN=-0.0000000005, $ XTICKNAME=[' ',' '] DEVICE, /CLOSE_file max_age_plot=12 ; make plot: DEVICE, XSize=9, YSize=7, /INCHES, /ENCAPSULATED, /DECOMPOSED, $ /COLOR, Bits_Per_Pixel=8, File='age_x_section.eps' PLOT, age, age, YRANGE=[0,0.915], $ XRANGE=[0,max_age_plot], ticklen=-0.02, $ /NODATA, $ XTHICK=T, YTHICK=T, YSTYLE=3, XSTYLE=3, $ XMINOR=2, YMINOR=2, YTICKS=3, XTICKS=4, $ XTITLE='age (years)', $ YTITLE='N infected per age bin', $ CHARTHICK=T, $ CHARSIZE=CS, $ YTICKV=[0,.3,.6,.9],$ POSITION=[0.2,0.2,0.95,0.9] xsection_color=$ [x11colors('blue'),$ ;0 x11colors('green'),$ ;1 x11colors('red'),$ ;2 x11colors('yellow'),$ ;3 x11colors('SeaGreen4'),$ ;4 x11colors('brown'),$ ;5 x11colors('violet'),$ ;6 x11colors('medium_blue'),$;7 x11colors('orange'),$ ;8 x11colors('olive_drab')] ;9 max_plot_age=max_age_plot*ipy-1 yleg=.6 FOR i=0,N_ELEMENTS(timestamp)-1 do begin j=(N_ELEMENTS(timestamp)-1)-i if (j EQ 0) OR (j EQ 2) OR (j EQ 4) OR (j EQ 8) then begin OPLOT, age[0:max_plot_age], $ age_x_section[0:max_plot_age,j], THICK=T,$ color=xsection_color[j] PLOTS, [.7,.8],[yleg,yleg], /norm, thick=T, color=xsection_color[j] XYOUTS,.807, yleg-.01, STRING(timestamp(j)/(1.0*ipy), FORMAT='(F3.1)'), $ charsize=1.5, charthick=1.2, /norm yleg=yleg+.05 endif endfor XYOUTS,.70, yleg, 'Time (years):', $ charsize=1.5, charthick=1.2, /norm PLOTS, .675, yleg+.05, /norm PLOTS, .875, yleg+.05, /norm, /CONTINUE, thick=T-2 PLOTS, .875, .56, /norm, /CONTINUE, thick=T-2 PLOTS, .675, .56, /norm, /CONTINUE, thick=T-2 PLOTS, .675, yleg+.05, /norm, /CONTINUE, thick=T-2 DEVICE, /CLOSE_file ; make plot: !P.FONT=1 ; for surfaces... ;; surface, zcoord , xcoord, ycoord, $ ;; YTITLE='age (years)', $ ;; ZTITLE='number affeced by rumor', $ ;; XTITLE='time (years)';, AZ=90;, ZRANGE=[0,200] ;;; a new attempt at the surface ;; how many cols and rows are there? numrows=N_ELEMENTS(measles_surface_downsamp3[0,*]) numcols=N_ELEMENTS(measles_surface_downsamp3[*,0]) agevec=25*(findgen(numcols)/numcols-1) timevec=10*(findgen(numrows)/numrows-1) DEVICE, XSize=9, YSize=7, /INCHES, /ENCAPSULATED, /DECOMPOSED, $ /COLOR, Bits_Per_Pixel=8, File='new_surf.eps' ; note that AZ=230 is the "optimal" setting ; for the data. however, the axes are drawn ; "behind" the data. ugly. so we do two ; things in tandem. (i) we get the axes right ; by reversing the data (see xrange yrange etc). ; (ii) the remedy in (i) fixes the axes but ; flips the data --- unflip it by setting ; AZ={orig_AZ}-180 (=50 in other words). done. SURFACE, measles_surface_downsamp3, $ agevec, timevec, $ AZ=50, AX=45, $ XSTYLE=3, YSTYLE=3, ZSTYLE=1, $ ZAXIS=1, $ XTHICK=T, YTHICK=T, ZTHICK=T, $ XRANGE = [MAX(agevec), MIN(agevec)], $ YRANGE=[MAX(timevec), MIN(timevec)], $ XTICKLEN=-0.04, YTICKLEN=-0.04, $ ZTICKLEN=-0.04, $ CHARSIZE=CS+2, /SAVE, $ yminor=2, $ ytickname=reverse(['0','2','4','6','8']), $ xminor=2, $ xtickname=reverse(['0','5','10','15','20']), $ XTITLE='age (years)', $ zminor=2, $ ztickname=replicate(' ',10), $ ztitle='incidence', zcharsize=CS-0.25 ; to display the y-axis label in such a way that ; the viewer does not need to turn her head to ; read the text, we use XYOUTS in place of an ; axis label. position is +/- by trial and error. ; NOTE: /T3d must be set to make the text look ; like it's in the same plane as the graph; otherwise ; it looks strange. Note also that /save must be called ; by the surface command or else the 3d transformation is ; not sved. XYOUTS, 7.5, -8 , 'time (years)', CHARSIZE=CS+2, /T3d, $ orientation=-90 DEVICE, /CLOSE_file ; ; ; now try in color ; DEVICE, XSize=9, YSize=7, /INCHES, /ENCAPSULATED, $ /COLOR, Bits_Per_Pixel=8, File='new_surf_color.eps' DEVICE, DECOMPOSED=0 ;; loadcolors ;; bottom=16 ;; loadct, 13, bottom=bottom bottom=byte(30); !d.table_size/2 loadct, 10, bottom=bottom top=byte(!d.table_size-70);- 100) ;; bottom=byte(20); !d.table_size/2 ;; loadct, 2, bottom=bottom ;; top=byte(!d.table_size-150);- 100) ;; bottom=byte(150); !d.table_size/2 ;; loadct, 11, bottom=bottom ;; top=byte(!d.table_size-1);- 100) ;; bottom=byte(1); !d.table_size/2 ;; loadct, 40, bottom=bottom ;; top=byte(!d.table_size-170);- 100) ;; TVLCT, r, g, b, /Get ;; TVLCT, Reverse(r), Reverse(g), Reverse(b) ceiling=percentiles(measles_surface_downsamp3,value=[0.925]) floor=percentiles(measles_surface_downsamp3,value=[0.1]) shades=BytScl(measles_surface_downsamp3, max=ceiling, top=top) $ + byte(bottom) shades=byte(255)-shades ; print, !d.table_size SURFACE, measles_surface_downsamp3, $ agevec, timevec, $ AZ=50, AX=45, $ ; COLOR=[0,0,0], $ ; color experiment shades=shades, $ ; end color experiment XSTYLE=3, YSTYLE=3, ZSTYLE=1, $ ZAXIS=1, $ XTHICK=T, YTHICK=T, ZTHICK=T, $ XRANGE = [MAX(agevec), MIN(agevec)], $ YRANGE=[MAX(timevec), MIN(timevec)], $ XTICKLEN=-0.04, YTICKLEN=-0.04, $ ZTICKLEN=-0.04, $ CHARSIZE=CS+2, /SAVE, $ yminor=2, $ ytickname=reverse(['0','2','4','6','8']), $ xminor=2, $ xtickname=reverse(['0','5','10','15','20']), $ XTITLE='age (years)', $ zminor=2, $ ztickname=replicate(' ',10), $ ztitle='incidence', zcharsize=CS-0.25 XYOUTS, 7.5, -8 , 'time (years)', CHARSIZE=CS+2, /T3d, $ orientation=-90 DEVICE, /CLOSE_file spawn, 'sh embed_labels.sh' ; TOC end