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('\n

Trend 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('\n


Usage 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\n')