*! Date : 12 Jan 2005 *! Version : 1.33 *! Author : Adrian Mander *! Email : adrian.p.mander@gsk.com program define surface7 preserve version 7.0 local varlist "ex min(3) max(3)" local options "SAVING(string) Box(string) EYE(string) ROUND(int 5) LABELRND(real 0.01) ORIENT(string) XLABel(string) YLABel(string) ZLAB(string) NOWIRE XTITLE(string) YTITLE(string) ZTITLE(string)" parse "`*'" parse "`varlist'", parse(" ") local x `1' local y `2' local z `3' if "`xtitle'"=="" { local xtitle "X-axis" } if "`ytitle'"=="" { local ytitle "Y-axis" } if "`ztitle'"=="" { local ztitle "Z-axis" } if "`orient'"=="" { local orient = "xyz" } else { if "`orient'"=="xyz" | "`orient'"=="xzy" | "`orient'"=="yxz" | "`orient'"=="yzx" | "`orient'"=="zxy" | "`orient'"=="zyx" { } else { di in red "orient must contain one of the following strings xyz, xzy, yxz, yzx,zxy,zyx" di in red "Re-setting orient to xyz" local orient = "xyz" } } /* Set up graph box variables leftx rightx topy ------------------------- | | | | height | | | | | | | | boty ------------------------- <- width -> */ global max_by=23063 global max_rx=32000 global topy=2000 global boty=21000 global leftx=3500 global rightx=30000 global width=$rightx-$leftx global height=$boty-$topy if "`box'" ~= "" { parse "`box'", parse(",") if(`1'<2000) { di "WARNING Bounding box wrong top y `1'<2000" exit(666) } if(`3'>21000) { di "WARNING Bounding box wrong bottom y `3'>21000" exit(666) } if(`5'<3500) { di "WARNING Bounding box wrong left x `5'<3500" exit(666) } if(`7'>30000) { di "WARNING Bounding box wrong right x `7'>30000" exit(666) } if(`1'>=`3' | `5'>=`7') { di "WARNING Bounding box wrong Is `1'>=`3'? OR is `5'>=`7'?" exit(666) } else { global topy=`1' global boty=`3' global leftx=`5' global rightx=`7' global width=$rightx-$leftx global height=$boty-$topy } } /* GET data in a square data set ready for drawing */ sq_dat `x' `y' `z', round(`round') /* * Calculate all the x y z data lengths */ mat origin = J(1,2,0) mat origin[1,1]=16000 mat origin[1,2]=12000 qui summ `y' global dly = _result(5) global dry = _result(6) qui summ `x' global dlx = _result(5) global drx = _result(6) qui summ `z' global dlz = _result(5) global drz = _result(6) mat xdir = J(1,2,0) mat ydir = J(1,2,0) mat zdir = J(1,2,0) /* Just initialise x/y/z dir matrices */ eyeball, orient(`orient') global scalex = 1 global scaley = 1 chek_box 1 /* di "Scaling $scaley $scalex" */ if $scaley~=1 | $scalex ~= 1 { mat xdir[1,1] = $scalex*xdir[1,1] mat ydir[1,1] = $scalex*ydir[1,1] mat zdir[1,1] = $scalex*zdir[1,1] mat xdir[1,2] = $scaley*xdir[1,2] mat ydir[1,2] = $scaley*ydir[1,2] mat zdir[1,2] = $scaley*zdir[1,2] } chek_box 2 /***************************************** * Saving the file or not. * NB if file exists it is deleted!! *****************************************/ if ("`saving'"~="") { gph open, saving(`saving') } else { gph open } draw_axe global labrnd= `labelrnd' lab_axe , xtitle(`xtitle') ytitle(`ytitle') ztitle(`ztitle') gph pen 2 if "`nowire'"~="" { local i 1 while `i' <= _N { draw_pt `x'[`i'] `y'[`i'] `z'[`i'] local i=`i'+1 } } else { if "`orient'"=="zxy" { draw_wis `x' `y' `z', orient(`orient') } else { draw_wis `x' `y' `z' } } gph close restore end /* This works out the perspective.. used to be that you could change the eyeball position * to anywhere but it morphed the axes box */ program define eyeball local options "ORIENT(string) " parse "`*'" local string1 = substr("`orient'",-3,1) local string2 = substr("`orient'",-2,1) local string3 = substr("`orient'",-1,1) mat `string1'dir[1,1] = 1 mat `string1'dir[1,2] =0 mat `string2'dir[1,1] = -0.5 mat `string2'dir[1,2] = 0.5 mat `string3'dir[1,1] = 0 mat `string3'dir[1,2] = -1 end program define draw_axe /* horizontal boxe lines */ gph pen 1 gph line $pt1y $pt1x $pt2y $pt2x gph line $pt1y $pt1x $pt3y $pt3x gph line $pt1y $pt1x $pt5y $pt5x gph line $pt2y $pt2x $pt4y $pt4x gph line $pt2y $pt2x $pt6y $pt6x gph line $pt3y $pt3x $pt7y $pt7x gph line $pt3y $pt3x $pt4y $pt4x gph line $pt4y $pt4x $pt8y $pt8x gph line $pt5y $pt5x $pt6y $pt6x gph line $pt5y $pt5x $pt7y $pt7x gph line $pt6y $pt6x $pt8y $pt8x gph line $pt7y $pt7x $pt8y $pt8x end /* Translate the x y z to the 2-D dimension and plot that point */ program define draw_pt local x1= origin[1,1] + (`1'-$dlx)*xdir[1,1] +(`2'-$dly)*ydir[1,1] local y1 = origin[1,2] + (`2'-$dly)*ydir[1,2] + (`1'-$dlx)*xdir[1,2] local x2 = `x1' + (`3'-$dlz)*zdir[1,1] local y2 = `y1'+ (`3'-$dlz)*zdir[1,2] gph line `y1' `x1' `y2' `x2' gph point `y2' `x2' 150 4 end /* Draw the wirefram translating the Z-axis to the 2-D image */ program define draw_wis local varlist "ex min(3) max(3)" local options "ORIENT(string) " parse "`*'" parse "`varlist'", parse(" ") local x `1' local y `2' local z `3' sort `x' `y' tempvar yy1 yy2 xx1 gen `yy1' =0 gen `yy2'=0 qui gen `xx1' = 0 local i 1 while `i'<=rowsof(xlines) { /* hold x and draw along y */ qui replace `yy1' = cond(`x'[_n]==xlines[`i',1], origin[1,2] + (`y'[_n]-$dly)*ydir[1,2] + (`x'[_n]-$dlx)*xdir[1,2] +(`z'[_n]-$dlz)*zdir[1,2], 0) qui replace `xx1' = cond(`x'[_n]==xlines[`i',1], origin[1,1] + (`y'[_n]-$dly)*ydir[1,1] + (`x'[_n]-$dlx)*xdir[1,1]+(`z'[_n]-$dlz)*zdir[1,1], 0) gph vline `yy1' `xx1' if `yy1'~=0 local i = `i'+1 } local i 1 while `i'<=rowsof(ylines) { /* hold x and draw along y */ qui replace `yy1' = cond(`y'[_n]==ylines[`i',1], origin[1,2] + (`y'[_n]-$dly)*ydir[1,2] + (`x'[_n]-$dlx)*xdir[1,2] +(`z'[_n]-$dlz)*zdir[1,2], 0) qui replace `xx1' = cond(`y'[_n]==ylines[`i',1], origin[1,1] + (`y'[_n]-$dly)*ydir[1,1] + (`x'[_n]-$dlx)*xdir[1,1]+(`z'[_n]-$dlz)*zdir[1,1], 0) gph vline `yy1' `xx1' if `yy1'~=0 local i = `i'+1 } end /************************************ * To scale up the axes ************************************/ program define chek_box /* need more options 8 pts to a cube take 1 as origin 3 4 ________ / / / / / / --------- 1 2 */ local i 1 while `i'<9 { matrix pt`i' = J(1,2,0) local i=`i'+1 } /* Fix lengths of x,y and z directions */ if "`1'"=="1" { local x1 = ($drx-$dlx)*xdir[1,1] local x2 = ($drx-$dlx)*xdir[1,2] local z1 = ($drz-$dlz)*zdir[1,1] local z2 = ($drz-$dlz)*zdir[1,2] local y1 = ($dry-$dly)*ydir[1,1] local y2 = ($dry-$dly)*ydir[1,2] if `x1'~=0 { mat xdir[1,1] = xdir[1,1]*abs(16000/`x1') } if `x2'~=0 { mat xdir[1,2] = xdir[1,2]*abs(11000/`x2') } if `y1'~=0 { mat ydir[1,1] = ydir[1,1]*abs(14000/`y1') } if `y2'~=0 { mat ydir[1,2] = ydir[1,2]*abs(10000/`y2') } if `z2'~=0 { mat zdir[1,2] = zdir[1,2]*abs(12000/`z2') } if `z1'~=0 { mat zdir[1,1] = zdir[1,1]*abs(15000/`z1') } } local x1 = ($drx-$dlx)*xdir[1,1] local x2 = ($drx-$dlx)*xdir[1,2] local z1 = ($drz-$dlz)*zdir[1,1] local z2 = ($drz-$dlz)*zdir[1,2] local y1 = ($dry-$dly)*ydir[1,1] local y2 = ($dry-$dly)*ydir[1,2] mat pt1[1,1]=origin[1,1] mat pt1[1,2]=origin[1,2] mat pt2[1,1]=origin[1,1]+ `x1' mat pt2[1,2]=origin[1,2]+ `x2' mat pt3[1,1]=origin[1,1]+ `z1' mat pt3[1,2]=origin[1,2]+ `z2' mat pt4[1,1]=origin[1,1]+ `z1'+ `x1' mat pt4[1,2]=origin[1,2]+ `z2'+ `x2' mat pt5[1,1]=origin[1,1]+ `y1' mat pt5[1,2]=origin[1,2]+ `y2' mat pt6[1,1]=origin[1,1]+ `x1' + `y1' mat pt6[1,2]=origin[1,2]+ `x2' + `y2' mat pt7[1,1]=origin[1,1]+ `z1' + `y1' mat pt7[1,2]=origin[1,2]+ `z2' + `y2' mat pt8[1,1]=origin[1,1]+ `z1' + `y1' + `x1' mat pt8[1,2]=origin[1,2]+ `z2' + `y2' + `x2' qui minax_8 1 global rbx= $S_max global lbx = $S_min qui minax_8 2 global tby= $S_min global bby = $S_max if "`1'"=="1" { local frame=4000 global scalex = ($max_rx-`frame')/($rbx-$lbx) global scaley = ($max_by-`frame')/($bby-$tby) mat origin[1,1]=int((origin[1,1]-$lbx)*$scalex)+3*`frame'/4 mat origin[1,2]=int((origin[1,2]-$tby)*$scaley)+`frame'/2 } end /************************************************ * Label axes ************************************************/ program define lab_axe syntax [varlist] [,XTITLE(string) YTITLE(string) ZTITLE(string)] minax_9 gph pen 1 /* X-Axis */ local labopt " xtitle(`xtitle') ytitle(`ytitle') ztitle(`ztitle')" if xdir[1,2]== 0 { labxyz, dirx(x) `labopt' } if ydir[1,2]== 0 { labxyz, dirx(y) `labopt' } if zdir[1,2]== 0 { labxyz, dirx(z) `labopt' } if xdir[1,1]== 0 { labxyz, dirz(x) `labopt' if xdir[1,2]~= 0 & xdir[1,1]~= 0{ labxyz, diry(x x) `labopt' } if ydir[1,2]~= 0 & ydir[1,1]~= 0{ labxyz, diry(x y) `labopt' } if zdir[1,2]~= 0 & zdir[1,1]~= 0{ labxyz, diry(x z) `labopt'} } if ydir[1,1]== 0 { labxyz, dirz(y) `labopt' if xdir[1,2]~= 0 & xdir[1,1]~= 0{ labxyz, diry(y x) `labopt'} if ydir[1,2]~= 0 & ydir[1,1]~= 0{ labxyz, diry(y y) `labopt'} if zdir[1,2]~= 0 & zdir[1,1]~= 0{ labxyz, diry(y z) `labopt'} } if zdir[1,1]== 0 { labxyz, dirz(z) `labopt' if xdir[1,2]~= 0 & xdir[1,1]~= 0{ labxyz, diry(z x) `labopt' } if ydir[1,2]~= 0 & ydir[1,1]~= 0{ labxyz, diry(z y) `labopt' } if zdir[1,2]~= 0 & zdir[1,1]~= 0{ labxyz, diry(z z) `labopt' } } end /* Do the axis labelling */ program define labxyz local options "DIRX(string) DIRZ(string) DIRY(string) XTITLE(string) YTITLE(string) ZTITLE(string)" parse "`*'" if "`dirx'"~="" { local dr = "dr`dirx'" local dl = "dl`dirx'" local y1 = pt$blcpt[1,2]+($`dr'-$`dl')*`dirx'dir[1,2] + 1000 local x1 = pt$blcpt[1,1]+($`dr'-$`dl')*`dirx'dir[1,1]/2 /* local text = "`dirx'-axis" */ local text "``dirx'title'" gph text `y1' `x1' 0 0 `text' local sty = pt$blcpt[1,2] local stx = pt$blcpt[1,1] local endy = pt$blcpt[1,2]+($`dr'-$`dl')*`dirx'dir[1,2] local endx = pt$blcpt[1,1]+($`dr'-$`dl')*`dirx'dir[1,1] ticks, line(`sty' `stx' `endy' `endx') dirx(`dirx') } if "`dirz'"~="" { local dr = "dr`dirz'" local dl = "dl`dirz'" local y1 = pt$blcpt[1,2]+($`dr'-$`dl')*`dirz'dir[1,2]/2 + 1000 local x1 = pt$blcpt[1,1]+($`dr'-$`dl')*`dirz'dir[1,1]- 1000 /* local text = "`dirz'-axis"*/ local text "``dirz'title'" gph text `y1' `x1' 0 0 `text' local sty = pt$blcpt[1,2] local stx = pt$blcpt[1,1] local endy = pt$blcpt[1,2]+($`dr'-$`dl')*`dirz'dir[1,2] local endx = pt$blcpt[1,1]+($`dr'-$`dl')*`dirz'dir[1,1] ticks, line(`sty' `stx' `endy' `endx') dirz(`dirz') } if "`diry'"~="" { parse "`diry'", parse(" ") local dr1 = "dr`1'" local dl1 = "dl`1'" local dr2 = "dr`2'" local dl2 = "dl`2'" local y2 = pt$blcpt[1,2]+($`dr1'-$`dl1')*`1'dir[1,2]-($`dr2'-$`dl2')*`2'dir[1,2]/2-1000 local x2 = pt$blcpt[1,1]+($`dr1'-$`dl1')*`1'dir[1,1]-($`dr2'-$`dl2')*`2'dir[1,1]/2-2000 /* local text = "`2'-axis" */ local text = "``2'title'" gph text `y2' `x2' 0 0 `text' local sty = pt$blcpt[1,2]+($`dr1'-$`dl1')*`1'dir[1,2] local stx = pt$blcpt[1,1]+($`dr1'-$`dl1')*`1'dir[1,1] local endy = pt$blcpt[1,2]+($`dr1'-$`dl1')*`1'dir[1,2]-($`dr2'-$`dl2')*`2'dir[1,2] local endx = pt$blcpt[1,1]+($`dr1'-$`dl1')*`1'dir[1,1]-($`dr2'-$`dl2')*`2'dir[1,1] ticks, line(`sty' `stx' `endy' `endx') diry(`2') } end /************************************************ * Draw in all the ticks..... ************************************************/ program define ticks local options "DIRX(string) DIRZ(string) DIRY(string) LINE(string)" parse "`*'" parse "`line'", parse(" ") if "`dirx'"~="" { local dr = "dr`dirx'" local dl = "dl`dirx'" local y1 = `1' local x1 = `2' local y2 = `1'+300 local text : di %9.3f $`dl' gph line `y1' `x1' `y2' `x1' local y2 =`y2'+700 gph text `y2' `x1' 0 0 `text' local dr = "dr`dirx'" local dl = "dl`dirx'" local y1 = `3' local x1 = `4' local y2 = `3'+300 local text : di %9.3f $`dr' gph line `y1' `x1' `y2' `x1' local y2 =`y2'+700 gph text `y2' `x1' 0 0 `text' } if "`dirz'"~="" { local dr = "dr`dirz'" local dl = "dl`dirz'" local y1 = `1' local x1 = `2' local x2 = `2'-400 local text :di %9.3f $`dl' gph line `y1' `x1' `y1' `x2' local x2 =`x2'-400 gph text `y1' `x2' 0 1 `text' local dr = "dr`dirz'" local dl = "dl`dirz'" local y1 = `3' local x1 = `4' local x2 = `4'-400 local text : di %9.3f $`dr' gph line `y1' `x1' `y1' `x2' local x2 =`x2'-400 gph text `y1' `x2' 0 1 `text' } if "`diry'"~="" { local dr = "dr`diry'" local dl = "dl`diry'" local y1 = `1' local x1 = `2' local y2 = `1'-300 local text : di %9.3f $`dr' gph line `y1' `x1' `y2' `x1' local y2 =`y2'-700 gph text `y2' `x1' 0 0 `text' local dr = "dr`diry'" local dl = "dl`diry'" local y1 = `3' local x1 = `4' local y2 = `3'-300 local text : di %9.3f $`dl' gph line `y1' `x1' `y2' `x1' local y2 =`y2'-600 gph text `y2' `x1' 0 0 `text' } end /************************************************ * find the number of macros ************************************************/ program define macno global macn=0 global last=0 local tmp_str "" local tmp_str "`2'" parse "`1'", parse("`2'") while "`1'"~="" { while "`1'"=="`tmp_str'" { mac shift global last = $last+1 } global macn = $macn+1 mac shift global last= $last+1 } end /************************************************ * This checks the state of the data and sees if * it can fit it in a squarer grid. * ************************************************/ program define sq_dat version 5.0 local varlist "ex min(3) max(3)" local options "ROUND(int 5)" parse "`*'" parse "`varlist'", parse(" ") tempname xx yy local error 0 sort `1' cap unique `1' if _rc==908 { di in red "ERROR set matsize higher" exit(908) } if _rc~=0 { di "problem with variable `1'" error _rc } if _rc==0 { mat xlines = resp } else { local error=_rc } sort `2' cap unique `2' if _rc==908 { di in red "ERROR set matsize higher" exit(908) } if _rc~=0 { di "problem with variable `2'" error _rc } if _rc==0 { mat ylines = resp } else { local error=_rc+`error' } sort `1' `2' local manyx = sqrt(2*_N) if rowsof(xlines)>`manyx' | `error'>0 { di in red "WARNING the x-variable contains too many unique values attempting to round" qui summ `1' local rnd = `round'*(_result(6)-_result(5))/_N qui gen `xx' = round(`1',`rnd') drop `1' rename `xx' `1' cap unique `1' if _rc==908 { di in red "ERROR set matsize higher OR increase the round() value" exit(908) } mat xlines = resp } local manyy = sqrt(2*_N) if rowsof(ylines)>`manyy' | `error'>0 { di in red "WARNING the y-variable contains too many unique values attempting to round" qui summ `2' local rnd = `round'*(_result(6)-_result(5))/_N gen `yy' = round(`2',`rnd') drop `2' rename `yy' `2' cap unique `2' if _rc==908 { di in red "ERROR set matsize higher OR increase the round() value" exit(908) } mat ylines = resp } end /************************************************ * Take the 8 vertices of the cube and try and * find the extremes in each direction. ************************************************/ program define minax_8 local ind `1' local inde = mod(`ind',2)+1 local min = 99999999 local max = -99999999 local maxpti`ind' = 0 local minpti`ind' = 0 local i=1 while `i'<9 { if `min'>=pt`i'[1,`ind'] { if `min'==pt`i'[1,`ind'] { if pt`i'[1,`inde'] < pt`minpti`ind''[1,`inde'] & `ind'==2 { local minpti`ind' = `i' local min=pt`i'[1,`ind'] } if pt`i'[1,`inde'] < pt`minpti`ind''[1,`inde'] & `ind'==1 { local minpti`ind' = `i' local min=pt`i'[1,`ind'] } } if abs(`min'-pt`i'[1,`ind'])>0.0001 { local min=pt`i'[1,`ind'] local minpti`ind' = `i' } } if `max'<=pt`i'[1,`ind'] { if `max'==pt`i'[1,`ind'] { if pt`i'[1,`inde'] < pt`maxpti`ind''[1,`inde'] { local maxpti`ind' = `i' } } else { local max=pt`i'[1,`ind'] local maxpti`ind' = `i' } } if `ind'==1 { global pt`i'x = pt`i'[1,`ind'] } else { global pt`i'y = pt`i'[1,`ind'] } local i=`i'+1 } if `ind'==1 { global maxpti1 = `maxpti`ind'' global minpti1 = `minpti`ind'' } else { global maxpti2 = `maxpti`ind'' global minpti2 = `minpti`ind'' } global S_min=`min' global S_max=`max' end /* Botttom left corner and work from there */ program define minax_9 local minx = pt1[1,1] local minpt = 1 global blcpt =1 local i=2 while `i'<9 { if abs(`minx'-pt`i'[1,1])<0.00001 | `minx'> pt`i'[1,1]{ if abs(`minx'-pt`i'[1,1])<0.00001 { if pt`i'[1,2] > pt`minpt'[1,2] { local minpt = `i' local minx=pt`i'[1,1] } } else { local minpt = `i' local minx=pt`i'[1,1] } } local i=`i'+1 } global blcpt = `minpt' end /************************************************ * Take the 3 directions and try and * find the extremes in each direction. ************************************************/ program define minax_3 local ind `1' local min = 99999999 local max = -99999999 global maxdir`ind' = "q" global mindir`ind' = "q" local dir = "x" local i=1 while `i'<=3 { local dr = "dr`dir'" local dl = "dl`dir'" if `min'>($`dr'-$`dl')*`dir'dir[1,`ind'] { local min = ($`dr'-$`dl')*`dir'dir[1,`ind'] global mindir`ind' = "`dir'" } if `max'< ($`dr'-$`dl')*`dir'dir[1,`ind'] { local max = ($`dr'-$`dl')*`dir'dir[1,`ind'] global maxdir`ind' = "`dir'" } if "`dir'"=="x" { local dir = "y" } else { local dir = "z" } local i=`i'+1 } end program define unique version 5.0 local varlist "req ex min(1) max(1)" parse "`*'" parse "`varlist'", parse(" ") local var `1' preserve sort `var' qui by `var': keep if _n==1 global S_1=_N qui drop if `var'==. global S_2=_N mkmat `var', matrix(resp) restore end