rm(list=ls()); aou='01720'; regi='ALA'; iyrf=1966; iyrl=2021; linear=1; yrver=2021 options(width=222) #trnhtm21 <- function(aou, regi, iyrf=1966, iyrl=2021, linear=1, yrver=2021) { # program to read in composite files and provide summary results jrs 2-22-2002 # modified for 2021 analysis by jeh = 4/21/2023 # *********Note..this one is modified to read 10000; NOT 10001 my screw-up in making the coda files in 2013. Gah. t1=Sys.time(); dbg=FALSE; cputime=codacnt=cpucoda=codagrp=st=ststr=sr=stnum=iregs=mflag=0; nreps=12000 if (iyrf<1900) { iyrf=iyrf+1900; dbg=TRUE; cat('
\n') } iyrf=iyrf-1965 iyrl=iyrl-1965 if(dbg)cat('aou=',aou,' regi=',regi,' yrf=',iyrf,iyrl,linear,yrver,'\n') yrs=yrver-1965 # 56 = 2021 - 1965 # ******************** input files needed: ******************** bugdir=paste0('/mnt/c/hines/csc1/bbs',yrver,'/bugs',yrver,'/') fin2=paste0(bugdir,'bug_input_2/scitab',yrver,'.txt') fin17=paste0(bugdir,'bug_input_2/ststr_firstyr_all.txt') # 17=old fortran unit number fin18=paste0(bugdir,'bug_input_2/r',aou,'.txt') if(yrver>=2021) fin18=paste0(bugdir,'bug_input_2/',aou,'r.txt') fin9=paste0(bugdir,'bug_outputs/code',aou,'.txt') fin10=paste0(bugdir,'bug_outputs/coda',aou,'.bin') #open(10,file=fin,status='OLD',form='unformatted',recl=nreps*4,err=1006) fin19=paste0(bugdir,'region/',regi) fout14=paste0(bugdir,'bug_input_p/',aou,'.txt') # ***************************************************************** # first, read files of weights and read of first years of data stnum=read.table(fin17) # open(17,file=fin,status='OLD',err=1001) a=readLines(fin2) # open(43,file=fin2,status='OLD',err=1002) iaou=substr(a,6,10); name=substr(a,14,53); sci=substr(a,55,94) mflag=which(aou==iaou); if (length(mflag)>0) { sci=sci[mflag]; name=name[mflag] # The approach is: # 1. Read in data (strata estimates) for each replicate # 2. Calculate composite indexes, trends for that replicate # 3. Store those puppies # 4. Sort them out, calculate summaries # 5. Bob's your uncle # Do not need to store state-stratum, just the states (or bcr's have to read in structure for each, then fill in from the read # tedium: structure for storing the results # now, read in region codes, to get the numbering of the regions iststr=1; a18=read.table(fin18) if (yrver < 2021) ststr=a18 else ststr=cbind(floor(a18[,1]/100), a18[,1]%%100, a18[,2], 0) ststr[,4]=0; iststr=nrow(a18) i=which((ststr[,1]*1000+ststr[,2]) %in% (stnum[,1]*1000+stnum[,2])) ststr[i,4]=stnum[i,3] # now, cut out regions that are not of interest 300 series suppressWarnings(a19 <- readLines(fin19)); lnum=0; a20=list() while(lnum0) { iflag=1; sr[jj,4]=ststr[jj,4] } if (iflag == 0) numin=0 # if iflag=0 go to next on list else { # otherwise, print it out to 20 iregs=iregs+1 a20[[iregs]]=list(numin=numin, statchr=statchr, sr=sr, stchr=stchr, st=st) # write(20,'(i3,a3)') numin,statchr } } # here is the read that repeats until list is read num=a20[[iregs]]$numin; statchr=a20[[iregs]]$statchr sr=a20[[iregs]]$sr; stchr=a20[[iregs]]$stchr; st=a20[[iregs]]$st; m=length(a20[[iregs]]$st) minyr = min(sr[sr[,4]!=0,4]) if (minyr > 3) minyr=3 iregs=iregs+1 # Here we calculate sample sizes by reading the raw data qqq rst,rstr,rcnt(50 rsamp,rsamp2,robs(50),rt2 a14 = read.table(fout14); n=nrow(a14) i=seq(1,n,2); j=seq(2,n,2) rst=floor(a14[j,1]/1000); rstr=a14[j,2]; rcnt=a14[i,-1:-2]; robs=a14[j,-1:-2] k = which((rst*1000+rstr) %in% (st*1000+sr[,1])) rsamp=sum(rcnt[k,iyrf:iyrl],na.rm=T) # so the list of states-strata are in. Match them to the indexes sumwts here, too. asumwt=0 for (j in 1:nrow(sr)) { i=which(st[j]==ststr[,1] & sr[j,1]==ststr[,2]) if (length(i)>0) { sr[j,3]=i; asumwt=asumwt+sr[j,2]} } a9 = readLines(fin9); i=grep('n\\[',a9); a9=a9[i] iiivec = as.numeric(gsub('.+\\[','',gsub(',.+','',a9))) jjjvec = as.numeric(gsub('.+,','',gsub('\\].+','',a9))) zz10=file(fin10, "rb"); codagrp=0 datain = matrix(0, nreps, yrs) # initialize the datain iflag = which(iiivec %in% sr[,3]) # this line checks to see if the stratum no (iii) matches # the no from the input reg list sr(i,3) # So, we sum by year for all the strata in the region. store the ones that match in datain by replicate and year # note that iflag= the region, to allow weight assignment t4=Sys.time() for (i in 1:length(iiivec)) { iflag = which(iiivec[i] == sr[,3]) # this line checks to see if the stratum no (iii) matches if (length(iflag)>0) { # the no from the input reg list sr(i,3) seek(zz10, (i-1)*(nreps+2)*4+4) cpucoda = cpucoda + system.time(xcoda <- readBin(zz10,what=double(), size=4, n=nreps)) # omit burn-in and omit strata not in list # sum to replicate (note that strata are summed within each replicate, with areas multiplied, # then divided afterwards by the total area. store these in datain indexed as i-nrepin # (nprein = burnin = 0 , so not used here.) datain[,jjjvec[i]] = datain[,jjjvec[i]] + xcoda * sr[iflag,2] codagrp=codagrp+1 #cat(i,iiivec[i],jjjvec[i],iflag); print(c(xcoda[1:2],sr[iflag,2],datain[1:2,jjjvec[i]])) } } close(zz10); codacnt=nreps*codagrp; cat('time reading coda files:',Sys.time()-t4,'\n') # this region has been read in, so now summarize by year. first divide each replicate by total region area asumwt datain = datain / asumwt; if (iyrf < minyr) iyrf=minyr # yrs+1 stores the trend firstyr is minyr*** datain1=(datain[,iyrl]/datain[,iyrf])^(1/(iyrl-iyrf)) datain2=(datain[,yrs]/datain[,minyr])^(1/(yrs-minyr)) datain=cbind(datain,datain1,datain2) amean=colMeans(datain); var=apply(datain,2,sd)^2 # sort datain..find min, switch with current value ilow=floor(nreps*0.025) # find the 5, 50 and 95 % points imed=round(nreps*0.5) iupp=nreps-ilow # note that all years are output# cat(' North American Breeding Bird Survey
\nSummary of Population Change (Ver',yrver,')
\n') cat('',name,'',sci,'
\n',regi,' trend results
\n') cat('Species aou:',aou,' Region:',regi,'
\n') if (regi == 'A24') cat('N of Survey Routes: ',rsamp,'
\n') cat('Hierarchical Model Results\n') cat('\nTrend period',iyrf+1965,'to',iyrl+1965,'\n') cpusort = system.time( o <- order(datain[,yrs+1])); o3=o[c(imed, ilow, iupp)] datain = datain[o,] tbl=matrix(c(round(100*(datain[c(imed, ilow, iupp),yrs+1]-1),2)),nrow=1) rownames(tbl)='region'; colnames(tbl)=c('Trend_Estimate','2.5%CI','97.5%CI') print(tbl) cat('\nAnnual Indices for Years:',minyr+1965,' to ',yrs+1965,'\nYear Annual Index 2.5% CI 97.5% CI\n') t3=Sys.time() tbl=NULL for (k in minyr:yrs) { cpusort = cpusort + system.time( o <- order(datain[,k])); o3=o[c(imed, ilow, iupp)] tbl=rbind(tbl,round(datain[o3,k],2)) } rownames(tbl)=1965+minyr:yrs; colnames(tbl)=c('Trend_Estimate','2.5%CI','97.5%CI') print(tbl) cat('time for annual indices:',Sys.time()-t3,'\n') if(codacnt == 0) cat('
No Results On File For This Species/Region#
\n') t2=Sys.time() cat('\nUsage Statistics:
\n
sort cpu time=',cpusort[3],'\n') cat('
cpu time reading coda files=',cpucoda[3],'\n
coda groups read=',codagrp,' lines=',codacnt,'\n') cat('
total cpu time',t2-t1,'
\n') } else cat('No species name found for ',aou,'\n