*! version 1.0 P.MILLAR 22Jul2012 *! version 1.1 P.MILLAR 02Aug2012 added probabilistic shapes *! This software can be used for non-commercial purposes only. *! The copyright is retained by the developer. *! Copyright 2005-2008 Paul Millar program define fractal, rclass byable(recall) version 12.0 syntax [anything] , [HShape(numlist sort min=3) VShape(numlist min=3) HRange(numlist sort min=2) VRange(numlist sort min=2) ITer(integer 1) KEEPVars SAVEGraph HSHAPE2(numlist sort min=3) VSHAPE2 (numlist min=3) PROB2(real 0) HSHAPE3(numlist sort min=3) VSHAPE3(numlist min=3) PROB3(real 0) ] tempvar x y if `iter' <1 { di as error "Number of iterations must be at least 1" exit } /* default Hrange and Vrange */ if "`hrange'" == "" { local hrange="0 100" } if "`vrange'" == "" { local vrange="0 100" } /* ----------------------------------------- */ /* deal with the probabilities of each shape */ /* ----------------------------------------- */ local nshapes=1 local sprob2=0 if "`prob2'" != "" { local sprob2=`prob2' if `sprob2' <0 | `sprob2' > 1 { di as error "PROB2 must be between 0 and 1" exit } local nshapes =`nshapes' + 1 } local sprob3=0 if "`prob3'" != "" { local sprob3=`prob3' if `sprob3' <0 | `sprob3' > 1 { di as error "PROB3 must be between 0 and 1" exit } local nshapes =`nshapes' + 1 } local prob23=`sprob2' + `sprob3' if `prob23' > 1 { di as error "The sum of PROB2 and PROB3 must be less than 1" exit } local sprob1=1-`prob23' local plevel1=`sprob1' local plevel2=`plevel1'+`sprob2' local plevel3=`plevel2'+`sprob3' /* find number of horizontal points (initial shape) */ tokenize `hshape' local word="`1'" local nhpoints=1 while "`word'" !="" { local hpoint`nhpoints'="`word'" local nhpoints=`nhpoints'+1 local word="``nhpoints'+1'" } local nhpoints=`nhpoints'-1 if `nhpoints' < 2 { di as error "No shape specified" exit } /* find number of vertical points (initial shape) */ tokenize `vshape' local word="`1'" local nvpoints=1 while "`word'" !="" { local vpoint`nvpoints'="`word'" local nvpoints=`nvpoints'+1 local word="``nvpoints'+1'" } local nvpoints=`nvpoints'-1 local nsegs=`nvpoints'-1 /* standardize vertical scale */ forvalues i=1/`nvpoints' { local vpoint`i'=`vpoint`i''/`vpoint`nvpoints'' } if `nhpoints' != `nvpoints' { di as error "Number of points in the horizontal shape must be the same as the number of points in the vertical shape" exit } local nsegs1=`nsegs' local nh1points=`nhpoints' forvalues i=1/`nh1points' { local h1point`i'=`hpoint`i'' local v1point`i'=`vpoint`i'' } /* --------------- */ /* process shape 2 */ /* --------------- */ tokenize `hshape2' local word="`1'" local nh2points=1 while "`word'" !="" { local h2point`nh2points'="`word'" local nh2points=`nh2points'+1 local word="``nh2points'+1'" } local nh2points=`nh2points'-1 tokenize `vshape2' local word="`1'" local nv2points=1 while "`word'" !="" { local v2point`nv2points'="`word'" local nv2points=`nv2points'+1 local word="``nv2points'+1'" } local nv2points=`nv2points'-1 local nsegs2=`nv2points'-1 /* standardize vertical scale */ forvalues i=1/`nv2points' { local v2point`i'=`v2point`i''/`v2point`nv2points'' } /* --------------- */ /* process shape 3 */ /* --------------- */ tokenize `hshape3' local word="`1'" local nh3points=1 while "`word'" !="" { local h3point`nh3points'="`word'" local nh3points=`nh3points'+1 local word="``nh3points'+1'" } local nh3points=`nh3points'-1 tokenize `vshape3' local word="`1'" local nv3points=1 while "`word'" !="" { local v3point`nv3points'="`word'" local nv3points=`nv3points'+1 local word="``nv3points'+1'" } local nv3points=`nv3points'-1 local nsegs3=`nv3points'-1 forvalues i=1/`nv3points' { local v3point`i'=`v3point`i''/`v3point`nv3points'' } /* decode the Hrange (same for all shapes) */ tokenize `hrange' local hmin="`1'" local hmax="`2'" local hlength=`hmax' - `hmin' /* decode the Vrange */ tokenize `vrange' local vmin="`1'" local vmax="`2'" local vlength=`vmax' - `vmin' qui gen `x'=. qui gen `y'=. /* ------------------------------------------------------------ */ /* make sure there are enough cases to cover all the iterations */ /* ------------------------------------------------------------ */ local maxpoints=0 forvalues i=1/`nshapes' { if `nh`i'points' > `maxpoints' { local maxpoints=`nh`i'points' } } local newobs=((`maxpoints'-1)^(`iter'))+1 local oldobs=_N // di "maxpoints=`maxpoints', newobs=`newobs', oldobs=`oldobs'" if `newobs' > _N { qui set obs `newobs' } /* --------------- */ /* first iteration */ /* --------------- */ local nsegs=`nhpoints'-1 local xlength=`hlength' local ylength=`vlength' forvalues i=1/`nhpoints' { qui replace `x'=`hmin' + `hpoint`i'' * `xlength' in `i' qui replace `y'=`vmin' + `vpoint`i'' * `ylength' in `i' } /* save the results of first iteration if asked to do so */ if "`keepvars'"=="keepvars" { capture qui gen _frctlx1=. if _rc !=0 { qui drop _frctlx1 qui gen _frctlx1=. } capture qui gen _frctly1=. if _rc != 0 { qui drop _frctly1 qui gen _frctly1=. } qui replace _frctlx1=`x' qui replace _frctly1=`y' } if "`savegraph'" == "savegraph" { line `y' `x', ytitle("Y") xtitle("X") qui graph save _frctl1,replace } local i=`nsegs'+1 /* -------------------------------- */ /* second and subsequent iterations */ /* -------------------------------- */ forvalues iter8=2/`iter' { local nsegs= `i'-1 // di "--------->iter8=`iter8', nsegs=`nsegs', i=`i'" forvalues seg=1/`nsegs' { // di "=====>iter8=`iter8', seg=`seg'" local j=`seg' local k=`seg'+1 local curxmin=`x' in `j' local curxmax=`x' in `k' local curxlength=`curxmax' - `curxmin' local curymin=`y' in `j' local curymax=`y' in `k' local curylength=`curymax' - `curymin' // di " curxmin=`curxmin', curxmax=`curxmax', curxlength=`curxlength'" // di " curymin=`curymin', curymax=`curymax', curylength=`curylength'" /* choose the shape to insert */ local rand=runiform() if `rand' < `plevel1' { local shapeno = 1 } else if `rand' < `plevel2' { local shapeno = 2 } else { local shapeno = 3 } // di "rand=`rand', plevel1=`plevel1', plevel2=`plevel2', plevel3=`plevel3', shape=`shapeno' local ncalcs=`nh`shapeno'points'-1 forvalues point=2/`ncalcs' { local i=`i' +1 // local newx=`curxmin' + `h`shapeno'point`point'' * `curxlength' // local newy=`curymin' + `v`shapeno'point`point'' * `curylength' // di " i=`i', newx=`newx', newy=`newy'" qui replace `x'=`curxmin' + `h`shapeno'point`point'' * `curxlength' in `i' qui replace `y'=`curymin' + `v`shapeno'point`point'' * `curylength' in `i' } } sort `x' /* save the results of this iteration if asked to do so */ if "`keepvars'"=="keepvars" { capture qui gen _frctlx`iter8'=`x' if _rc!=0 { qui drop _frctlx`iter8' qui gen _frctlx`iter8'=`x' } capture qui gen _frctly`iter8'=`y' if _rc!=0 { qui drop _frctly`iter8' qui gen _frctly`iter8'=`y' } } if "`savegraph'" == "savegraph" { line `y' `x' , ytitle("Y") xtitle("X") qui graph save _frctl`iter8', replace } local nsegs=(`nh`shapeno'points'-1)*`nsegs' // di " end of loop, nsegs=`nsegs'" } /* end of loop */ capture qui gen _frctlx=`x' if _rc !=0 { qui drop _frctlx qui gen _frctlx=`x' } qui label var _frctlx "X values of fractal" capture qui gen _frctly=`y' if _rc != 0 { qui drop _frctly qui gen _frctly=`y' } qui label var _frctly "Y values of fractal" /* ---------------------------- */ /* Remove unneeded observations */ /* ---------------------------- */ qui summ _frctlx local curobs=r(N) local bign=_N if `curobs' < _N { if `oldobs' > `curobs' { qui drop in `oldobs'/`bign' } else { qui drop in `curobs'/`bign' } } di " " di as result "Fractal generated with `iter' iterations" if "`savegraph'" == "savegraph" { di "To display list of graphs type -graph dir-" } di " " end