# Data cleaning
dataclean_TW <- function(dat1,rmssp=F) {
 # hist(dat1$hbf,nclass=400)
  dat1 <- dat1[!is.na(dat1$hooks),] #
#  hist(dat1$hooks,nclass=250)
#  dat1 <- dat1[dat1$hooks<10000,] # clean up outliers
  dat1 <- dat1[dat1$hooks<4500,] # clean up outliers
  dat1 <- dat1[dat1$hooks>200,]
#  dat1 <- dat1[dat1$yft<750,]
#  dat1 <- dat1[dat1$bet<250,]
#  dat1 <- dat1[dat1$alb<250,]
#  dat1 <- dat1[is.na(dat1$hbf)==F,]
#  dat1 <- dat1[dat1$op_yr > 1976,]
#  dat1 <- dat1[dat1$yrqtr < 2017,]
#  dat1 <- dat1[dat1$EW==1,]
#  dat1 <- dat1[dat1$hooks >= 1000,]
  dat1[dat1$hbf %in% c(35,155,20000),"hbf"] <- 15
  dat1[dat1$hbf %in% c(26,30),"hbf"] <- 20
  dat1[dat1$hbf %in% c(25),"hbf"] <- 24
  lenzero <- function(x) sum(x > 0)
  if(rmssp) {
    ssp <- apply(dat1[,c("alb","bet","yft","ott","swo","mls","bum","blm","otb","skj","sha","oth","sbt")],1,lenzero)
    dat1 <- dat1[ssp > 1,]
  }
#  dat1$hbf <- as.factor(as.character(dat1$hbf))
  return(dat1)
  }

dataclean_TW_std <- function(dat1, doHBF=F) {
  dat1 <- dat1[!is.na(dat1$hooks),]
  dat1 <- dat1[dat1$hooks<10000,]
  dat1 <- dat1[dat1$hooks>1000,]
  dat1 <- dat1[dat1$yft + dat1$bet + dat1$alb > 0,]
  if(doHBF) dat1 <- dat1[dat1$hbf <= 25,]
  lenzero <- function(x) sum(x > 0)
  ssp <- apply(dat1[,c("alb","bet","yft","ott","swo","mls","bum","blm","otb","skj","sha","oth","sbt")],1,lenzero)
  dat1 <- dat1[ssp > 1,]
  remvec <- c("G")
  dat1 <- dat1[-grep(remvec,dat1$Remark),]
  return(dat1)
  }

dataclean_KR <- function(dat) {
 # hist(dat$hbf,nclass=400)
  if(sum(is.na(dat$alb))>0) dat[is.na(dat$alb),]$alb <- 0
  if(sum(is.na(dat$bet))>0) dat[is.na(dat$bet),]$bet <- 0
  if(sum(is.na(dat$blm))>0) dat[is.na(dat$blm),]$blm <- 0
  if(sum(is.na(dat$bum))>0) dat[is.na(dat$bum),]$bum <- 0
  if(sum(is.na(dat$mls))>0) dat[is.na(dat$mls),]$mls <- 0
  if(sum(is.na(dat$oth))>0) dat[is.na(dat$oth),]$oth <- 0
  if(sum(is.na(dat$pbf))>0) dat[is.na(dat$pbf),]$pbf <- 0
  if(sum(is.na(dat$sbt))>0) dat[is.na(dat$sbt),]$sbt <- 0
  if(sum(is.na(dat$sfa))>0) dat[is.na(dat$sfa),]$sfa <- 0
  if(sum(is.na(dat$sha))>0) dat[is.na(dat$sha),]$sha <- 0
  if(sum(is.na(dat$skj))>0) dat[is.na(dat$skj),]$skj <- 0
  if(sum(is.na(dat$swo))>0) dat[is.na(dat$swo),]$swo <- 0
  if(sum(is.na(dat$yft))>0) dat[is.na(dat$yft),]$yft <- 0
  dat <- dat[!is.na(dat$hooks),] # Clean up 294 NAs
  dat <- dat[dat$hooks<5000,] # clean up outliers
  dat <- dat[dat$hooks>200,]
 # dat <- dat[dat$yft<750,]
#  dat <- dat[dat$bet<250,]
#  dat <- dat[dat$alb<250,]
  dat <- dat[is.na(dat$hbf)==F,]
  dat <- dat[dat$op_yr > 1976,]
  dat <- dat[dat$yrqtr < 2017,]
  dat <- dat[dat$EW==1,]
  dat <- dat[dat$hooks >= 1000,]
  dat <- dat[dat$hbf >= 5,]
  return(dat)
  }

dataclean_JPIO <- function(dat,checktg=F,allHBF=F) {
  dat$op_yr <- as.numeric(dat$op_yr)
  dat$op_mon <- as.numeric(dat$op_mon)
  dat$op_day <- as.numeric(dat$op_day)
  dat <- dat[dat$op_day < 32,]
  dat <- dat[!is.na(dat$op_day),]
  dmy <- as.Date(paste(dat$op_yr,dat$op_mon,dat$op_day,sep="-"))
  dat <- dat[!is.na(dmy),]
  dat$lat <- as.numeric(dat$lat)
  dat <- dat[!is.na(dat$lat),]
  dat <- dat[!is.na(dat$lon),]
  dat$latcode <- as.numeric(dat$latcode)
  dat$lon <- as.numeric(dat$lon)
  dat$loncode <- as.numeric(dat$loncode)
  dat <- dat[dat$loncode %in% c(1,2),]
  dat$hbf <- as.numeric(dat$hbf)
#  dat$tonnage <- as.numeric(dat$tonnage)
  dat$hooks <- as.numeric(dat$hooks)
  dat$alb <- as.numeric(dat$alb)
  dat$bet <- as.numeric(dat$bet)
  dat$yft <- as.numeric(dat$yft)
  dat$swo <- as.numeric(dat$swo)
  dat$sbt <- as.numeric(dat$sbt)
  dat$blm <- as.numeric(dat$blm)
  dat$bum <- as.numeric(dat$bum)
  dat$mls <- as.numeric(dat$mls)
  if(sum(is.na(dat$alb))>0) dat[is.na(dat$alb),]$alb <- 0
  if(sum(is.na(dat$bet))>0) dat[is.na(dat$bet),]$bet <- 0
  if(sum(is.na(dat$yft))>0) dat[is.na(dat$yft),]$yft <- 0
  if(sum(is.na(dat$swo))>0) dat[is.na(dat$swo),]$swo <- 0
  if(sum(is.na(dat$sbt))>0) dat[is.na(dat$sbt),]$sbt <- 0
  if(sum(is.na(dat$blm))>0) dat[is.na(dat$blm),]$blm <- 0
  if(sum(is.na(dat$bum))>0) dat[is.na(dat$bum),]$bum <- 0
  if(sum(is.na(dat$mls))>0) dat[is.na(dat$mls),]$mls <- 0
  dat <- dat[!is.na(dat$hooks),]
  dat <- dat[dat$hooks<5000,] # clean up outliers
  dat <- dat[dat$hooks>200,]
#  dat <- dat[dat$yft<250,]
#  dat <- dat[dat$bet<250,]
#  dat <- dat[dat$alb<250,]
#  dat <- dat[dat$tonnage<30000 | is.na(dat$tonnage),]
#  dat[dat$fishingcat =="0",]
#  dat <- dat[dat$fishingcat !=".",]
#  dat <- dat[dat$fishingcat !="0",]
#  dat <- dat[dat$hbf != "  .",]
#  dat <- dat[is.na(dat$hbf)==FALSE | dat$op_yr < 1976,]
#  dat <- dat[dat$hbf < 26 | is.na(dat$hbf)==TRUE,]
  dat <- dat[(!is.na(dat$hbf) & dat$hbf < 26) | (dat$op_yr < 1976 & is.na(dat$hbf)),]
#  if(allHBF==F) {
#    dat[dat$hbf>22,]$hbf <- 22     # pool hbf > 22 into 22
#    dat <- dat[dat$hbf > 4,]   # remove swordfish targeting in R1 and R2
#    }
#  dat$ncrew <- as.numeric(dat$ncrew)
  if(checktg) dat <- dat[dat$target == 3 | is.na(dat$target),] # tuna target  (remove to avoid a change in 1994 - but recent trend is more important)
  return(dat)
  }

dataclean <- function(dat,checktg=F,allHBF=F) {
  dat$op_yr <- as.numeric(dat$op_yr)
  dat$op_mon <- as.numeric(dat$op_mon)
  dat$op_day <- as.numeric(dat$op_day)
  dat <- dat[dat$op_day < 32,]
  dat <- dat[!is.na(dat$op_day),]
  dat$lat <- as.numeric(dat$lat)
  dat$latcode <- as.numeric(dat$latcode)
  dat$lon <- as.numeric(dat$lon)
  dat$loncode <- as.numeric(dat$loncode)
  dat <- dat[dat$loncode %in% c(1,2),]
  dat$hbf <- as.numeric(dat$hbf)
#  dat$tonnage <- as.numeric(dat$tonnage)
  dat$hooks <- as.numeric(dat$hooks)
  dat$alb <- as.numeric(dat$alb)
  dat$bet <- as.numeric(dat$bet)
  dat$yft <- as.numeric(dat$yft)
  dat$swo <- as.numeric(dat$swo)
  dat$sbt <- as.numeric(dat$sbt)
  dat$blm <- as.numeric(dat$blm)
  dat$bum <- as.numeric(dat$bum)
  dat$mls <- as.numeric(dat$mls)
  if(sum(is.na(dat$alb))>0) dat[is.na(dat$alb),]$alb <- 0
  if(sum(is.na(dat$bet))>0) dat[is.na(dat$bet),]$bet <- 0
  if(sum(is.na(dat$yft))>0) dat[is.na(dat$yft),]$yft <- 0
  if(sum(is.na(dat$swo))>0) dat[is.na(dat$swo),]$swo <- 0
  if(sum(is.na(dat$sbt))>0) dat[is.na(dat$sbt),]$sbt <- 0
  if(sum(is.na(dat$blm))>0) dat[is.na(dat$blm),]$blm <- 0
  if(sum(is.na(dat$bum))>0) dat[is.na(dat$bum),]$bum <- 0
  if(sum(is.na(dat$mls))>0) dat[is.na(dat$mls),]$mls <- 0
  dat <- dat[!is.na(dat$hooks),]
  dat <- dat[dat$hooks<10000,] # clean up outliers
  dat <- dat[dat$hooks>200,]
  dat <- dat[dat$yft<250,]
  dat <- dat[dat$bet<250,]
  dat <- dat[dat$alb<250,]
  dat <- dat[dat$tonnage<30000 | is.na(dat$tonnage),]
#  dat[dat$fishingcat =="0",]
#  dat <- dat[dat$fishingcat !=".",]
  dat <- dat[dat$fishingcat !="0",]
#  dat <- dat[dat$hbf != "  .",]
  dat <- dat[is.na(dat$hbf)==F | dat$op_yr < 1976,]
  dat <- dat[dat$hbf < 26 | is.na(dat$hbf)==T,]
#  if(allHBF==F) {
#    dat[dat$hbf>22,]$hbf <- 22     # pool hbf > 22 into 22
#    dat <- dat[dat$hbf > 4,]   # remove swordfish targeting in R1 and R2
#    }
#  dat$ncrew <- as.numeric(dat$ncrew)
  if(checktg) dat <- dat[dat$target == 3 | is.na(dat$target),] # tuna target  (remove to avoid a change in 1994 - but recent trend is more important)
  return(dat)
  }

# Data preparation
dataprep_TW <- function(dat,alldat=F) {
  dat$dmy <- as.Date(paste(dat$op_yr,dat$op_mon,dat$op_day,sep="-"))
  a <- dat$embark_date
  a <- gsub(" ","0",a)
  dat$embark_dmy <- as.Date(a,format="%Y%m%d")
  a <- dat$debark_date
  a <- gsub(" ","0",a)
  dat$debark_dmy <- as.Date(a,format="%Y%m%d")
  a <- dat$op_start_date
  a <- gsub(" ","0",a)
  dat$op_start_dmy <- as.Date(a,format="%Y%m%d")
  a <- dat$op_end_date
  a <- gsub(" ","0",a)
  dat$op_end_dmy <- as.Date(a,format="%Y%m%d")
#  hist(dat$op_start_dmy,breaks="months")
  dat$lat <- dat$op_lat
  dat$lon <- dat$op_lon

  dat$tonnage <- as.factor(substring(dat$callsign,1,1))
  a <- levels(dat$tonnage)
  levs <- cbind(c("0","1","2","3","4","5","6","7","8"),c("< 5","5-10","10-20","20-50","50-100","100-200","200-500","500-1000",">= 1000"))
  levels(dat$tonnage) <- levs[match(a,levs[,1]),2]
#  table(dat$tonnage,substring(dat$callsign,1,1))

#  table(dat$NS,useNA="always")
#  table(dat$EW,useNA="always")

  dat$lat[is.na(dat$NS)==F & dat$NS%in% c(3,7)] <- (dat$lat[is.na(dat$NS)==F & dat$NS%in% c(3,7)]+1) * -1
  dat$lon[is.na(dat$EW)==F & dat$EW==2] <- 360 - (dat$lon[is.na(dat$EW)==F & dat$EW==2] + 1)
  dat <- dat[is.na(dat$lon) | dat$lon >= 10,]
  dat <- dat[is.na(dat$lon) | dat$lon < 130,]
  dat <- dat[is.na(dat$lat) | dat$lat < 29,]

  la <- as.integer(substring(dat$op_area,1,2))
  lo <- as.integer(substring(dat$op_area,3,4))
  dat$lat5 <- 5*(la - 73)/2+2.5
  dat$lat5[la%%2 == 0] <- -(5 * (la[la%%2 == 0] - 74))/2-2.5   #

  dat$lon5 <- -(5*(lo - 1)/2+2.5)
  dat$lon5[lo%%2 == 0] <- 5*(lo[lo%%2 == 0] - 2)/2 + 2.5
  # Indian ocean regions for YFT and BET?
  #regY   1 is N of 10N and W of 80,
  dat$regY <- 0
  dat[dat$lat5 >=  10 & dat$lon5 < 80,]$regY <- 1
  dat[dat$lat5 <  10 & dat$lat5 >=  -10 & dat$lon5 >= 40 & dat$lon5 < 75,]$regY <- 2
  dat[dat$lat5 <  -10 & dat$lat5 >=  -15 & dat$lon5 >= 60 & dat$lon5 < 75,]$regY <- 2
  dat[dat$lat5 <  -10 & dat$lat5 >= -30 & dat$lon5 >= 20 & dat$lon5 < 60,]$regY <- 3
  dat[dat$lat5 <  -30 & dat$lat5 >= -40 & dat$lon5 >= 20 & dat$lon5 < 40,]$regY <- 3
  dat[dat$lat5 <  -15 & dat$lat5 >= -40 & dat$lon5 >= 60 & dat$lon5 < 120,]$regY <- 4
  dat[dat$lat5 <  -30 & dat$lat5 >= -40 & dat$lon5 >= 40 & dat$lon5 < 60,]$regY <- 4
  dat[dat$lat5 < 10 & dat$lat5 >= -15 & dat$lon5 >= 75 & dat$lon5 < 100,]$regY <- 5
  dat[dat$lat5 < -5 & dat$lat5 >= -15 & dat$lon5 >= 100 & dat$lon5 < 110,]$regY <- 5
  dat[dat$lat5 < -10 & dat$lat5 >= -15 & dat$lon5 >= 110 & dat$lon5 < 130,]$regY <- 5
  dat[dat$lat5 < 30 & dat$lat5 >= 10 & dat$lon5 >= 80 & dat$lon5 < 100,]$regY <- 6
  #regB   North of 15S and west of 80 is R1, or north of  20 and west of 45; north of 15S and east of 80 is R2; north of 35S is R3
  dat$regB <- 0
  dat[dat$lat5 <  10 & dat$lat5 >=  -15 & dat$lon5 >= 20 & dat$lon5 < 80,]$regB <- 1
  dat[dat$lat5 <  10 & dat$lat5 >=  -20 & dat$lon5 >= 20 & dat$lon5 < 46,]$regB <- 1
  dat[dat$lat5 < 10 & dat$lat5 >= -15 & dat$lon5 >= 80 & dat$lon5 < 100,]$regB <- 2
  dat[dat$lat5 < -3 & dat$lat5 >= -15 & dat$lon5 >= 100 & dat$lon5 < 110,]$regB <- 2
  dat[dat$lat5 < -7 & dat$lat5 >= -15 & dat$lon5 >= 110 & dat$lon5 < 130,]$regB <- 2
  dat[dat$lat5 <  -20 & dat$lat5 >= -35 & dat$lon5 >= 20 & dat$lon5 < 120,]$regB <- 3
  dat[dat$lat5 <  -15 & dat$lat5 >= -35 & dat$lon5 >= 46 & dat$lon5 < 120,]$regB <- 3

  dat$bt1 <- !(is.na(dat$bait1) | dat$bait1==0)
  dat$bt2 <- !(is.na(dat$bait2) | dat$bait2==0)
  dat$bt3 <- !(is.na(dat$bait3) | dat$bait3==0)
  dat$bt4 <- !(is.na(dat$bait4) | dat$bait4==0)
  dat$bt5 <- !(is.na(dat$bait5) | dat$bait5==0)

  dat$vessid <- as.factor(as.numeric(as.factor(dat$callsign)))
  dat$tripid <- as.factor(paste(dat$vessid,dat$op_start_date))
  dat$tripidmon <- as.factor(paste(dat$vessid,dat$op_yr,dat$op_mon))
  dat$moon <- lunar.illumination(dat$dmy)

  dat$yrqtr <- dat$op_yr + floor((dat$op_mon-1)/3)/4 + 0.125
  dat$latlong <- paste(dat$lat5,dat$lon5,sep="_")
  dat$sbt <- dat$sbf + dat$pbf
  dat$sbt_w <- dat$sbf_w + dat$pbf_w
  dat$Total <- with(dat,alb + bet + yft + pbf + sbf + ott + swo + mls + bum + blm + otb + skj + sha + oth)
  dat$Total2 <- apply(dat[,c("bet","yft","alb")],1,sum)
  noms <- c("vessid","yrqtr","latlong","op_yr","op_mon","hbf","hooks","tonnage","tripid","tripidmon","moon",
  "alb","bet","yft","ott","swo","mls","bum","blm","otb","skj","sha","oth","sbt","Total","Total2",
  "alb_w","bet_w","yft_w","ott_w","swo_w","mls_w","bum_w","blm_w","otb_w","skj_w","sha_w","oth_w","sbt_w",
  "sst","bt1","bt2","bt3","bt4","bt5","hookdp","target","rem",
  "op_end_date","dmy","embark_dmy","debark_dmy","op_start_dmy","op_end_dmy","lat","lon","lat5","lon5","regY","regB","oil","foc")
  dat <- dat[,noms]
  return(dat)
  }

dataprep_KR <- function(dat,alldat=F) {
  dat$dmy <- as.Date(dat$DATE)
  dat$hbf <- round(dat$hooks / dat$floats,0)
  dat$moon <- lunar.illumination(dat$dmy)

  dat$lat <- dat$Lat01
  dat$lon <- dat$Long01
  dat$lat[dat$NS==2] <- (dat$Lat01[dat$NS==2]+1) * -1
  dat$lon[dat$EW==2] <- 360 - (dat$Long01[dat$EW==2] + 1)
  dat <- dat[dat$lon >= 0,]
  dat <- dat[dat$lat < 29,]

  dat$lat5 <- 5 * floor(dat$lat/5) + 2.5
  dat$lon5 <- 5 * floor(dat$lon/5) + 2.5

  # Indian ocean regions for YFT and BET?
  #regY   1 is N of 10N and W of 80,
  dat$regY <- 0
  dat[dat$lat >=  10 & dat$lon < 80,]$regY <- 1
  dat[dat$lat <  10 & dat$lat >=  -10 & dat$lon >= 40 & dat$lon < 75,]$regY <- 2
  dat[dat$lat <  -10 & dat$lat >=  -15 & dat$lon >= 60 & dat$lon < 75,]$regY <- 2
  dat[dat$lat <  -10 & dat$lat >= -30 & dat$lon >= 20 & dat$lon < 60,]$regY <- 3
  dat[dat$lat <  -30 & dat$lat >= -40 & dat$lon >= 20 & dat$lon < 40,]$regY <- 3
  dat[dat$lat <  -15 & dat$lat >= -40 & dat$lon >= 60 & dat$lon < 120,]$regY <- 4
  dat[dat$lat <  -30 & dat$lat >= -40 & dat$lon >= 40 & dat$lon < 60,]$regY <- 4
  dat[dat$lat < 10 & dat$lat >= -15 & dat$lon >= 75 & dat$lon < 100,]$regY <- 5
  dat[dat$lat < -5 & dat$lat >= -15 & dat$lon >= 100 & dat$lon < 110,]$regY <- 5
  dat[dat$lat < -10 & dat$lat >= -15 & dat$lon >= 110 & dat$lon < 130,]$regY <- 5
  dat[dat$lat < 30 & dat$lat >= 10 & dat$lon >= 80 & dat$lon < 100,]$regY <- 6
  #regB   North of 15S and west of 80 is R1, or north of  20 and west of 45; north of 15S and east of 80 is R2; north of 35S is R3
  dat$regB <- 0
  dat[dat$lat <  10 & dat$lat >=  -15 & dat$lon >= 20 & dat$lon < 80,]$regB <- 1
  dat[dat$lat <  10 & dat$lat >=  -20 & dat$lon >= 20 & dat$lon < 46,]$regB <- 1
  dat[dat$lat < 10 & dat$lat >= -15 & dat$lon >= 80 & dat$lon < 100,]$regB <- 2
  dat[dat$lat < -3 & dat$lat >= -15 & dat$lon >= 100 & dat$lon < 110,]$regB <- 2
  dat[dat$lat < -7 & dat$lat >= -15 & dat$lon >= 110 & dat$lon < 130,]$regB <- 2
  dat[dat$lat <  -20 & dat$lat >= -35 & dat$lon >= 20 & dat$lon < 120,]$regB <- 3
  dat[dat$lat <  -15 & dat$lat >= -35 & dat$lon >= 46 & dat$lon < 120,]$regB <- 3

  dat$yrqtr <- dat$op_yr + floor((dat$op_mon-1)/3)/4 + 0.125
  dat$latlong <- paste(dat$lat5,dat$lon5,sep="_")
  dat$vessid <- as.factor(as.numeric(dat$VESSEL_NAME_rev))
  #dat$vessid <- as.factor(as.numeric(dat$VESSEL_CD))
  dat$tripidmon <- paste(dat$vessid,dat$op_yr,dat$op_mon)
  dat$Total <- with(dat,alb + bet + yft + sbt + pbf + sfa + mls + bum + blm + swo + sha + skj + oth)
  dat$Total2 <- apply(dat[,c("bet","yft","alb")],1,sum)
  return(dat)
  }

dataprep_KRSBT <- function(dat,alldat=F) {
  dat$dmy <- as.Date(dat$DATE)
  dat$hbf <- round(dat$hooks / dat$floats,0)
  dat$moon <- lunar.illumination(dat$dmy)
  dat$lat <- dat$Lat01
  dat$lon <- dat$lonx <- dat$Long01
  dat$lat[dat$NS==2] <- (dat$Lat01[dat$NS==2]+1) * -1
  dat$lon[dat$EW==2] <- 360 - (dat$Long01[dat$EW==2] + 1)
  dat$lonx[dat$EW==2] <- -(dat$Long01[dat$EW==2] + 1)
  dat <- dat[dat$lon >= 0,]
  dat <- dat[dat$lat < 29,]

  dat$lat5 <- 5 * floor(dat$lat/5) + 2.5
  dat$lon5 <- 5 * floor(dat$lon/5) + 2.5
  dat$lonx5 <- 5 * floor(dat$lonx/5) + 2.5

  #regSBT   preliminary boudnaries - update later
  dat$regSBT <- 0
  dat[dat$lat <  -30 & dat$lon < 60,]$regSBT <- 1
  dat[dat$lat <  -30 & dat$lon >= 60,]$regSBT <- 2
  dat$lonfac <- factor(dat$lonx,levels=seq(-25,190,1))
  dat$latfac <- factor(dat$lat,levels=seq(-60,10,1))
  dat$lon5fac <- factor(dat$lonx,levels=seq(-25,190,1))
  dat$lat5fac <- factor(dat$lat,levels=seq(-60,10,1))

  dat$qtr <- floor((dat$op_mon-1)/3)+1
  dat$yrqtr <- dat$op_yr + floor((dat$op_mon-1)/3)/4 + 0.125
  dat$latlong <- paste(dat$lat5,dat$lon5,sep="_")
  dat$vessid <- as.factor(as.numeric(dat$VESSEL_NAME_rev))
  dat$tripidmon <- paste(dat$vessid,dat$op_yr,dat$op_mon)
  return(dat)
  }

dataprep_JPIO <- function(dat,alldat=T) {
  dat$dmy <- as.Date(paste(dat$op_yr,dat$op_mon,dat$op_day,sep="-"))
  dat$tripidmon <- as.factor(paste(dat$vessid,dat$op_yr,dat$op_mon))
  dat$moon <- lunar.illumination(dat$dmy)

  dat$lat_raw <- dat$lat
  dat$lon_raw <- dat$lon
  dat$lat[dat$latcode==2] <- (dat$lat_raw[dat$latcode==2]+1) * -1
  dat$lon[dat$loncode==2] <- 360 - (dat$lon_raw[dat$loncode==2] + 1)
  dat <- dat[dat$lon >= 0,]

  dat$lat5 <- 5 * floor(dat$lat/5) + 2.5
  dat$lon5 <- 5 * floor(dat$lon/5) + 2.5

  dat$regY <- 0
  dat[dat$lat5 >=  10 & dat$lon5 < 80 & !is.na(dat$lat5),]$regY <- 1
  dat[dat$lat5 <  10 & dat$lat5 >=  -10 & dat$lon5 >= 40 & dat$lon5 < 75 & !is.na(dat$lat5),]$regY <- 2
  dat[dat$lat5 <  -10 & dat$lat5 >=  -15 & dat$lon5 >= 60 & dat$lon5 < 75 & !is.na(dat$lat5),]$regY <- 2
  dat[dat$lat5 <  -10 & dat$lat5 >= -30 & dat$lon5 >= 20 & dat$lon5 < 60 & !is.na(dat$lat5),]$regY <- 3
  dat[dat$lat5 <  -30 & dat$lat5 >= -40 & dat$lon5 >= 20 & dat$lon5 < 40 & !is.na(dat$lat5),]$regY <- 3
  dat[dat$lat5 <  -15 & dat$lat5 >= -40 & dat$lon5 >= 60 & dat$lon5 < 120 & !is.na(dat$lat5),]$regY <- 4
  dat[dat$lat5 <  -30 & dat$lat5 >= -40 & dat$lon5 >= 40 & dat$lon5 < 60 & !is.na(dat$lat5),]$regY <- 4
  dat[dat$lat5 < 10 & dat$lat5 >= -15 & dat$lon5 >= 75 & dat$lon5 < 100 & !is.na(dat$lat5),]$regY <- 5
  dat[dat$lat5 < -5 & dat$lat5 >= -15 & dat$lon5 >= 100 & dat$lon5 < 110 & !is.na(dat$lat5),]$regY <- 5
  dat[dat$lat5 < -10 & dat$lat5 >= -15 & dat$lon5 >= 110 & dat$lon5 < 130 & !is.na(dat$lat5),]$regY <- 5
  dat[dat$lat5 < 30 & dat$lat5 >= 10 & dat$lon5 >= 80 & dat$lon5 < 100 & !is.na(dat$lat5),]$regY <- 6
  #regB   North of 15S and west of 80 is R1, or north of  20 and west of 45; north of 15S and east of 80 is R2; north of 35S is R3
  dat$regB <- 0
  dat[dat$lat5 <  10 & dat$lat5 >=  -15 & dat$lon5 >= 20 & dat$lon5 < 80 & !is.na(dat$lat5),]$regB <- 1
  dat[dat$lat5 <  10 & dat$lat5 >=  -20 & dat$lon5 >= 20 & dat$lon5 < 46 & !is.na(dat$lat5),]$regB <- 1
  dat[dat$lat5 < 10 & dat$lat5 >= -15 & dat$lon5 >= 80 & dat$lon5 < 100 & !is.na(dat$lat5),]$regB <- 2
  dat[dat$lat5 < -3 & dat$lat5 >= -15 & dat$lon5 >= 100 & dat$lon5 < 110 & !is.na(dat$lat5),]$regB <- 2
  dat[dat$lat5 < -7 & dat$lat5 >= -15 & dat$lon5 >= 110 & dat$lon5 < 130 & !is.na(dat$lat5),]$regB <- 2
  dat[dat$lat5 <  -20 & dat$lat5 >= -35 & dat$lon5 >= 20 & dat$lon5 < 120 & !is.na(dat$lat5),]$regB <- 3
  dat[dat$lat5 <  -15 & dat$lat5 >= -35 & dat$lon5 >= 46 & dat$lon5 < 120 & !is.na(dat$lat5),]$regB <- 3

  dat$callsign[dat$callsign == "      "] <- ".     "
  dat$vessid <- as.numeric(as.factor(paste(dat$callsign)))
  if (alldat==F) dat <- dat[dat$vessid != 1,]
  dat$vessid <- as.numeric(as.factor(dat$vessid))

  dat$yrqtr <- dat$op_yr + floor((dat$op_mon-1)/3)/4 + 0.125
  dat$latlong <- paste(dat$lat5,dat$lon5,sep="_")
  dat$Total <- with(dat,alb + bet + yft + sbt + swo + mls + bum + blm)
  dat$Total2 <- apply(dat[,c("bet","yft","alb")],1,sum)

  dat$trip_yr <- as.numeric(substr(as.character(dat$trip_st),1,4))
#  dat <- dat[dat$trip_yr > 1945 | is.na(dat$trip_yr),]

  dat$tripid <- paste(dat$vessid,dat$trip_st,sep="_")
  dat$tripid[dat$vessid == 1] <- NA
  dat$tripid[dat$trip_st == 0] <- NA

  dat$hbf[dat$op_yr < 1976 & is.na(dat$hbf)] <- 5
#  a <- table(dat$tripid)
#  dat$sets_per_trip <- NA
#  dat$sets_per_trip <- a[match(dat$tripid,names(a))]
  noms <- c("vessid","yrqtr","latlong","op_yr","op_mon","hbf","hooks","tripid","tripidmon","trip_yr","moon",
  "alb","bet","yft","swo","mls","bum","blm","sbt","Total","Total2","dmy","lat","lon","lat5","lon5","regY","regB")
  dat <- dat[,noms]

  return(dat)
  }

dataprep <- function(dat,alldat=F) {
  dat$lat_raw <- dat$lat
  dat$lon_raw <- dat$lon
  dat$lat[dat$latcode==2] <- (dat$lat_raw[dat$latcode==2]+1) * -1
  dat$lon[dat$loncode==2] <- 360 - (dat$lon_raw[dat$loncode==2] + 1)
  dat <- dat[dat$lon >= 0,]

  dat$lat5 <- 5 * floor(dat$lat/5)
  dat$lon5 <- 5 * floor(dat$lon/5)

  dat$reg <- 0
  dat[dat$lat <  40 & dat$lat >=  20 & dat$lon >= 110 & dat$lon < 170,]$reg <- 1
  dat[dat$lat <  40 & dat$lat >=  20 & dat$lon >= 170 & dat$lon < 210,]$reg <- 2
  dat[dat$lat <  20 & dat$lat >= -10 & dat$lon >= 110 & dat$lon < 170,]$reg <- 3
  dat[dat$lat <  20 & dat$lat >= -10 & dat$lon >= 170 & dat$lon < 210,]$reg <- 4
  dat[dat$lat < -10 & dat$lat >= -40 & dat$lon >= 140 & dat$lon < 170,]$reg <- 5
  dat[dat$lat < -10 & dat$lat >= -40 & dat$lon >= 170 & dat$lon < 210,]$reg <- 6
  dat[dat$lat <  40 & dat$lat >=  20 & dat$lon >= 210,]$reg <- 7
  dat[dat$lat <  20 & dat$lat >= -40 & dat$lon >= 210,]$reg <- 8

  dat$subreg <- 0
  dat[dat$lat <  40 & dat$lat >=  20 & dat$lon >= 110 & dat$lon < 170,]$subreg <- 1
  dat[dat$lat <  40 & dat$lat >=  20 & dat$lon >= 170 & dat$lon < 210,]$subreg <- 2
  dat[dat$lat <  20 & dat$lat >=   0 & dat$lon >= 110 & dat$lon < 150,]$subreg <- 3.1
  dat[dat$lat <  20 & dat$lat >=   0 & dat$lon >= 150 & dat$lon < 170,]$subreg <- 3.2
  dat[dat$lat <   0 & dat$lat >= -10 & dat$lon >= 110 & dat$lon < 150,]$subreg <- 3.3
  dat[dat$lat <   0 & dat$lat >= -10 & dat$lon >= 150 & dat$lon < 170,]$subreg <- 3.4
  dat[dat$lat <  20 & dat$lat >= -10 & dat$lon >= 170 & dat$lon < 180,]$subreg <- 4.1
  dat[dat$lat <  20 & dat$lat >= -10 & dat$lon >= 180 & dat$lon < 210,]$subreg <- 4.2
  dat[dat$lat < -10 & dat$lat >= -40 & dat$lon >= 140 & dat$lon < 170,]$subreg <- 5
  dat[dat$lat < -10 & dat$lat >= -40 & dat$lon >= 170 & dat$lon < 210,]$subreg <- 6
  dat[dat$lat <  40 & dat$lat >=  20 & dat$lon >= 210,]$subreg <- 7
  dat[dat$lat <  20 & dat$lat >= -40 & dat$lon >= 210,]$subreg <- 8


  dat$callsign[dat$callsign == "      "] <- ".     "
  dat$vessid <- as.numeric(as.factor(paste(dat$callsign)))
  if (alldat==F) dat <- dat[dat$vessid != 1,]
  dat$vessid <- as.numeric(as.factor(dat$vessid))

  dat$yrqtr <- dat$op_yr + floor((dat$op_mon)/3)/4 + 0.125
  dat$latlong <- paste(dat$lat5,dat$lon5,sep=".")
  dat <- dat[dat$yrqtr < 2012,]
  #dat <- dat[dat$reg %in% 1:6,]

  dat$newfishingcat <- NA
  dat <- dat[dat$fishingcat<=3,]
  dat <- dat[dat$op_yr < 1967 | dat$op_yr > 1970 | dat$fishingcat < 3,]
  dat <- dat[dat$op_yr <= 1957 | dat$fishingcat != ".",]
  dat[dat$op_yr <=1957 & dat$reg %in% c(1),]$newfishingcat <- 1
  dat[dat$op_yr <=1957 & dat$reg %in% c(0,2:8),]$newfishingcat <- 2
  dat[dat$op_yr >1957 & dat$op_yr<=1993 & dat$fishingcat==1,]$newfishingcat <- 1
  dat[dat$op_yr >1993 & dat$fishingcat==3,]$newfishingcat <- 1
  dat[dat$op_yr >1957 & dat$op_yr<=1966 & dat$fishingcat %in% c(2,3),]$newfishingcat <- 2
  dat[dat$op_yr >1966 & dat$op_yr<=1970 & dat$fishingcat %in% c(2),]$newfishingcat <- 2
  dat[dat$op_yr >=1971 & dat$op_yr<=1993 & dat$fishingcat %in% c(2,3),]$newfishingcat <- 2
  dat[dat$op_yr >1993 & dat$fishingcat %in% c(1,2),]$newfishingcat <- 2

  dat <- dat[dat$yrqtr > 1945,]

  dat$trip_yr <- as.numeric(substr(as.character(dat$trip_st),1,4))
  dat <- dat[dat$trip_yr > 1945 | is.na(dat$trip_yr),]

  dat$tripid <- paste(dat$vessid,dat$trip_st,sep="_")
  dat$tripid[dat$vessid == 1] <- NA
  dat$tripid[dat$trip_st == "       0"] <- NA

  a <- table(dat$tripid)
  dat$sets_per_trip <- NA
  dat$sets_per_trip <- a[match(dat$tripid,names(a))]
  return(dat)
  }


get.coefs <- function(model,nyrs) {
  coefs <- exp(c(0,summary(model)$coefficients[2:nyrs,1]))
  coefs <- coefs/mean(coefs)
  return(coefs)
  }

get.bin.coefs2 <- function(model,nyrs,dat) {
  mn <- logit(mean(dat[,3]!=0))
  a <- c(0,summary(model)$coefficients[2:nyrs,1])
  a <- a - mean(a) + mn
  coefs <- inv.logit(a)
#  coefs <- coefs/mean(coefs)
  return(coefs)
  }

get.bin.coefs <- function(model,nyrs,dat) {
  mn <- logit(mean(model$y))
  a <- c(0,summary(model)$coefficients[2:nyrs,1])
  a <- a - mean(a) + mn
  coefs <- inv.logit(a)
  return(coefs)
  }

get.coefs.summ <- function(modelsumm,nyrs) {
  coefs <- exp(c(0,modelsumm$coefficients[2:nyrs,1]))
  coefs <- coefs/mean(coefs)
  return(coefs)
  }

get.summ.coefs <- function(model,nyrs) {
  coefs <- exp(c(0,model$coefficients[2:nyrs,1]))
  coefs <- coefs/mean(coefs)
  }

get.bin.coefs.summ <- function(modelsumm,nyrs,dat) {
  mn <- logit(mean(modelsumm$deviance.resid > -0.001))
  a <- c(0,modelsumm$coefficients[2:nyrs,1])
  a <- a - mean(a) + mn
  coefs <- inv.logit(a)
  return(coefs)
  }


#get.coefs <- function(model,op_yr) {
#  a <- grep("yrqtr",names(model$coefficients))
#  estyrs <- as.numeric(gsub("as.factor(yrqtr)","",names(model$coefficients)[a],fixed=T))
#  coefs <- exp(c(0,summary(model)$coefficients[2:nyrs,1]))
#  coefs <- coefs/mean(coefs)
#  }

plot.slope.ratio <- function(coefs1,coefs2,yr,titl) {
  # base goes first, boat goes second
  windows(height=14,width=12)
  par(mfrow=c(2,1),mar=c(4,4,3,1))
  plot(yr,coefs1/coefs2,xlim=c(1976,2010),ylim=c(0,2),xlab="Year",ylab="Ratio of coefficients")
  title(main=titl,cex.main=1.2)
  lin <-lm(coefs1/coefs2 ~ yr)
  logl <-lm(log(coefs1/coefs2) ~ yr)
  lines(yr,exp(logl$coefficients[1] + logl$coefficients[2]*yr),lty=3)
  tt <- paste(prettyNum(100*logl$coefficients[2],digits=2,format="f"),
  "% \261 ",prettyNum(100*summary(logl)$coefficients[2,2],digits=2,format="f"),", p = ",
  prettyNum(summary(logl)$coefficients[2,4],
  digits=2,format="f"),sep="")
  text(min(yr)+5, 1.9, tt,font=2,col="red",cex=1.1)
  par(mar=c(5,4,1,1))
  plot(yr,coefs1,type="l",ylab="Relative abundance estimate",xlab="Year",ylim=c(0,2.5))
  lines(yr,coefs2,col="red")
  }

plot.slope.ratio2 <- function(coefs1,coefs2,yr1,yr2,titl) {
  # base goes first, boat goes second
  windows(height=14,width=12)
  par(mfrow=c(2,1),mar=c(4,4,3,1))
  yy <- c(yr1,yr2)
  yr <- seq(min(yy),max(yy),0.25)
  coefs1 <- coefs1[match(yr,yr1)]
  coefs2 <- coefs2[match(yr,yr2)]
  plot(yr,coefs1/coefs2,ylim=c(0,2),xlab="Year",ylab="Ratio of coefficients")
  title(main=titl,cex.main=1.2)
  lin <-lm(coefs1/coefs2 ~ yr)
  logl <-lm(log(coefs1/coefs2) ~ yr)
  lines(yr,exp(logl$coefficients[1] + logl$coefficients[2]*yr),lty=3)
  tt <- paste(prettyNum(100*logl$coefficients[2],digits=2,format="f"),
  "% \261 ",prettyNum(100*summary(logl)$coefficients[2,2],digits=2,format="f"),", p = ",
  prettyNum(summary(logl)$coefficients[2,4],digits=2,format="f"),sep="")
  text(min(yr)+10, 1.9, tt,font=2,col="red",cex=1.1)
  par(mar=c(5,4,1,1))
  plot(yr,coefs1,type="l",ylab="Relative abundance estimate",xlab="Year",ylim=c(0,3.5))
  lines(yr,coefs2,col="red")
  }

plot.agg.slope.ratio <- function(coefs1,coefs2,aggyr,opyr,titl,lab1="coefs1",lab2="coefs2",fname=NULL) {
  # agg goes first, op goes second
  windows(height=14,width=12)
  par(mfrow=c(2,1),mar=c(4,4,3,1))
  myr <- aggyr[match(opyr,aggyr)]
  coefs1 <- coefs1[match(opyr,aggyr)]
  plot(myr,coefs1/coefs2,ylim=c(0,2),xlab="Year",ylab="Ratio of coefficients")
  title(main=titl,cex.main=1.2)
  legend("bottomright",legend=paste(lab1,"/",lab2),col=1,lty=3)
  lin <-lm(coefs1/coefs2 ~ myr)
  logl <-lm(log(coefs1/coefs2) ~ myr)
  lines(myr,exp(logl$coefficients[1] + logl$coefficients[2]*myr),lty=3)
  tt <- paste(prettyNum(100*logl$coefficients[2],digits=2,format="f"),
  "% \261 ",prettyNum(100*summary(logl)$coefficients[2,2],digits=2,format="f"),", p = ",
  prettyNum(summary(logl)$coefficients[2,4],digits=2,format="f"),sep="")
  text(min(as.numeric(as.character(myr)))+10, 1.9, tt,font=2,col="red",cex=1.1)
  par(mar=c(5,4,1,1))
  plot(myr,coefs1,type="l",ylab="Relative abundance estimate",xlab="Year",ylim=c(0,3.5))
  lines(myr,coefs2,col="red")
  legend("topright",legend=c(lab1,lab2),col=1:2,lty=c(1,1))
  if(is.null(fname)==F) savePlot(paste(fname,lab1,lab2),type="png")
  }

plot.agg.slope.ratio_2000 <- function(coefs1,coefs2,aggyr,opyr,titl,lab1="coefs1",lab2="coefs2") {
  # agg goes first, op goes second
  windows(height=14,width=12)
  par(mfrow=c(2,1),mar=c(4,4,3,1))
  myr <- aggyr[match(opyr,aggyr)]
  coefs1 <- coefs1[match(opyr,aggyr)]
  plot(myr,coefs1/coefs2,xlim=c(1976,2010),ylim=c(0,2),xlab="Year",ylab="Ratio of coefficients")
  title(main=titl,cex.main=1.2)
  legend("bottomright",legend=paste(lab1,"/",lab2),col=1,lty=3)

  lin <-lm(coefs1/coefs2 ~ myr)
  logl <-lm(log(coefs1/coefs2) ~ myr)
  logl2000 <-lm(log(coefs1[myr<=2000]/coefs2[myr<=2000]) ~ myr[myr<=2000])
  lines(myr[myr<=2000],exp(logl2000$coefficients[1] + logl2000$coefficients[2]*myr[myr<=2000]),lty=3)
  tt <- paste(prettyNum(100*logl2000$coefficients[2],digits=2,format="f"),
  "% \261 ",prettyNum(100*summary(logl2000)$coefficients[2,2],digits=2,format="f"),", p = ",
  prettyNum(summary(logl2000)$coefficients[2,4],
  digits=2,format="f"),sep="")
  text(min(myr)+5, 1.9, tt,font=2,col="red",cex=1.1)
  par(mar=c(5,4,1,1))
  plot(myr,coefs1,type="l",ylab="Relative abundance estimate",xlab="Year",ylim=c(0,2.5))
  lines(myr,coefs2,col="red")
  legend("topleft",legend=c(lab1,lab2),col=1:2,lty=c(1,1))
  }

plot.res <- function(summ1,summ2,name1,name2) {
  nyrs <- length(grep("yrqtr",rownames(summ1$coefficients)))+1
  coefs1 <- get.summ.coefs(summ1,nyrs)
  coefs2 <- get.summ.coefs(summ2,nyrs)
  yrs <- rownames(summ1$coefficients)[grep("yrqtr",rownames(summ1$coefficients))]
  yrs <- c(1980.25,as.numeric(substring(yrs,17)))
  windows(height=13,width=11)
  par(mfrow=c(2,1),mar=c(5,4,1,2))
  plot.slope.ratio(coefs1,coefs2,yrs,"")
  plot(yrs,coefs1,type="l",ylab="Relative abundance estimate",xlab="Year",ylim=c(0,2.5))
  lines(yrs,coefs2,col="red")
  legend("topleft",legend=c(name1, name2),lty=c(1,1),col=c("black","red"))
}

make_formula <- function(runsp,modtype,addboat,splitboat=F,addmain=F,addbranch=F,addother=F,addalb=F) {
  fmla <- "~ as.factor(yrqtr) + as.factor(latlong) + poly(hbf,6)"
  modhead <- switch(modtype,logn=paste("log(",runsp,"/hooks+0.01)"),deltabin=paste(runsp,"!= 0"),deltapos=paste("log((",runsp,")/hooks)"),qp=runsp,propn=runsp)
  fmla <- paste(modhead,fmla)
  if(modtype %in% c("deltabin","qp")) fmla <- paste(fmla,"+ ns(hooks,6)")
  if(addboat & !splitboat) fmla <- paste(fmla,"+ as.factor(vessid)")
  if(addboat & splitboat) fmla <- paste(fmla,"+ as.factor(splitvess)")
  if(addmain) fmla <- paste(fmla,"+ as.factor(mainline) + as.factor(mainline):ns(hbf,6)")
  if(addbranch) fmla <- paste(fmla,"+ as.factor(branchline)")
  if(addother) fmla <- paste(fmla,"+ as.factor(other)")
  if(addalb) fmla <- paste(fmla,"+ as.factor(alb_cat)")
  return(fmla)
  }

make_formula_IO <- function(runsp,modtype,addboat,dohbf=T,splitboat=F,addcl=NA,addpca=NA) {
  fmla <- "~ yrqtr + latlong"
  if(dohbf) fmla <- paste(fmla,"+ ns(hbf,10)")
  modhead <- switch(modtype,logn=paste("log(",runsp,"/hooks+mn)"),deltabin=paste(runsp,"!= 0"),deltapos=paste("log((",runsp,")/hooks)"),qp=runsp,propn=runsp)
  fmla <- paste(modhead,fmla)
#  if(modtype %in% c("deltabin","qp")) fmla <- paste(fmla,"+ ns(hooks,10)")
  fmla <- paste(fmla,"+ ns(hooks,10)")
  if(addboat & !splitboat) fmla <- paste(fmla,"+ vessid")
  if(addboat & splitboat) fmla <- paste(fmla,"+ splitvess")
  if(!is.na(addcl)) fmla <- paste(fmla,"+",addcl)
  if(!is.na(addpca)) fmla <- paste(fmla,"+",addpca)
  return(fmla)
  }

make_formula_SBT <- function(runsp,modtype,addboat,dohbf=T,splitboat=F,addcl=NA,addpca=NA,addhooks=T,addmoon=F,addmonth=F) {
  fmla <- "~ op_yr + latlong"
  if(dohbf) fmla <- paste(fmla,"+ ns(hbf,10)")
  modhead <- switch(modtype,logn=paste0("log(",runsp,"/hooks+mn)"),deltabin=paste(runsp,"!= 0"),negbin=runsp,deltapos=paste("log((",runsp,")/hooks)"),qp=runsp,propn=runsp)
  fmla <- paste(modhead,fmla)
  if(addhooks) fmla <- paste(fmla,"+ ns(hooks,10)")
  if(addboat & !splitboat) fmla <- paste(fmla,"+ vessid")
  if(addboat & splitboat) fmla <- paste(fmla,"+ splitvess")
  if(addmonth) fmla <- paste(fmla,"+ ns(op_mon,df=4)")
  if(addmoon) fmla <- paste(fmla,"+ ns(moon,df=4)")
  if(!is.na(addcl)) fmla <- paste(fmla,"+",addcl)
  if(!is.na(addpca)) fmla <- paste(fmla,"+",addpca)
  return(fmla)
  }

aggregate_data <- function(dat,sp) {
  tab1 <- aggregate(dat[,3], list(dat$yrqtr, dat$latlong, dat$hbf), sum)
  tab2 <- aggregate(dat$hooks, list(dat$yrqtr, dat$latlong, dat$hbf), sum)
  taball <- cbind(tab1[,1:4],tab2[,4])
  names(taball) <-  c("yrqtr", "latlong","hbf",sp,"hooks")
  taball$yrqtr <- as.numeric(as.character(taball$yrqtr))
  taball$latlong <- as.character(taball$latlong)
  taball$hbf <- as.numeric(as.character(taball$hbf))
  taball[,4] <- as.numeric(as.character(taball[,4]))
  taball$hooks <- as.numeric(as.character(taball$hooks))
  return(taball)
  }

mk_wts <- function(dat,wttype,catch=NULL,sp=NULL) {
  if(wttype=="equal") wts <- NULL
  if(wttype=="propn") wts <- catch
  if(wttype=="area") {
    a <- tapply(dat$latlong,list(dat$latlong,dat$yrqtr),length)
    i <- match(dat$latlong,rownames(a))
    j <- match(dat$yrqtr,colnames(a))
    n <- mapply("[", list(a), i, j)
    wts <- 1/n
    }
  if(wttype=="catch") {
    if(is.null(catch)) catch <- tapply(dat[,sp],list(dat$latlong),sum)
    a <- tapply(dat$latlong,list(dat$latlong,dat$yrqtr),length)
    i <- match(dat$latlong,rownames(a))
    j <- match(dat$yrqtr,colnames(a))
    n <- mapply("[", list(a), i, j)
    cwts <- mapply("[", list(catch), i)/sum(catch)
    wts <- cwts/n
    }
  return(wts)
  }

mk_wts_op_yr <- function(dat,wttype,catch=NULL,sp=NULL) {
  if(wttype=="equal") wts <- NULL
  if(wttype=="propn") wts <- catch
  if(wttype=="area") {
    a <- tapply(dat$latlong,list(dat$latlong,dat$op_yr),length)
    i <- match(dat$latlong,rownames(a))
    j <- match(dat$op_yr,colnames(a))
    n <- mapply("[", list(a), i, j)
    wts <- 1/n
    }
  if(wttype=="catch") {
    if(is.null(catch)) catch <- tapply(dat[,sp],list(dat$latlong),sum)
    a <- tapply(dat$latlong,list(dat$latlong,dat$op_yr),length)
    i <- match(dat$latlong,rownames(a))
    j <- match(dat$op_yr,colnames(a))
    n <- mapply("[", list(a), i, j)
    cwts <- mapply("[", list(catch), i)/sum(catch)
    wts <- cwts/n
    }
  return(wts)
  }

mk_wts_integer <- function(dat,wttype,catch=NULL) {
  if(wttype=="equal") wts <- NULL
  if(wttype=="area") {
    a <- tapply(dat$latlong,list(dat$latlong,dat$yrqtr),length)
    i <- match(dat$latlong,rownames(a))
    j <- match(dat$yrqtr,colnames(a))
    n <- mapply("[", list(a), i, j)
    wts <- 1/n
    wts <- floor(10 * max(n) * wts)
    }
  if(wttype=="catch") {
    if(is.null(catch)) catch <- tapply(dat$bet,list(dat$latlong),sum)
    a <- tapply(dat$latlong,list(dat$latlong,dat$yrqtr),length)
    i <- match(dat$latlong,rownames(a))
    j <- match(dat$yrqtr,colnames(a))
    n <- mapply("[", list(a), i, j)
    cwts <- mapply("[", list(catch), i)/sum(catch)
    wts <- cwts/n
    }
  return(wts)
  }

make_strat <- function(dat) {
  a <- tapply(dat$latlong,list(dat$latlong,dat$yrqtr),length)
  i <- match(dat$latlong,rownames(a))
  j <- match(dat$yrqtr,colnames(a))
  n <- mapply("[", list(a), i, j)
  return(n)
  }

samp_data <- function(dat,n,nsamp) {
  p <- nsamp / n
  r <- runif(length(n))
  d2 <- dat[r<p,]
  return(d2)
  }

samp_data2 <- function(dat,p) {
  nn <- dim(dat)[1]
  nsamp <- floor(nn * p)
  r <- runif(nn)
  d1 <- dat[order(r),]
  d2 <- d1[1:nsamp,]
  return(d2)
  }

# p(bet !=0) ~ yrqtr + latlong + poly(hbf,6) + ns(hooks,6) + vessel id
# log(bet) ~ yrqtr + latlong + poly(hbf,6) + vessel id

# select the dataset
select_data <- function(indat,runreg,runsp,mt,minqtrs=2,maxqtrs=500,addmain=F,addother=F,addalb=F,fc="both",bait="no23",llstrat=5,doboth=F,addcl=F) {
  gdat <- indat[indat$reg==runreg,]
  if(runreg==9) gdat <- indat[indat$reg==3 & indat$lat >= -5 & indat$lat <10,]
  if(llstrat!=5) gdat$latlong <- paste(llstrat*floor(gdat$lat/llstrat),llstrat*floor(gdat$lon/llstrat),sep="_")
  if(mt=="deltapos") gdat <- gdat[gdat[,runsp] > 0,]
  if(bait=="no23") gdat <- gdat[(gdat$bait!=2 & gdat$bait !=3) | is.na(gdat$bait),]
  a <- table(gdat$vessid,gdat$yrqtr)
  a <- apply(a>0,1,sum)
  table(a)
  a <- a[a >= minqtrs & a <= maxqtrs]
  gdat <- gdat[gdat$vessid %in% names(a),]
  a <- table(gdat$yrqtr);a
  a <- a[a>=100]
  gdat <- gdat[gdat$yrqtr %in% names(a),]
  a <- table(gdat$vessid);a
  a <- a[a>=100]
  gdat <- gdat[gdat$vessid %in% names(a),]
  if(sum(is.na(gdat$hbf)) > 0) gdat[is.na(gdat$hbf),]$hbf <- 5
  gdat <- gdat[gdat$hbf >= 5,]
  fcnum <- switch(fc,"both"=0,"OS"=1,"DW"=2)
  if (fc!="both") gdat <- gdat[gdat$newfishingcat==fcnum,]
  if(addmain) {
    gdat <- gdat[gdat$target==3,c("vessid","hooks","yft", "bet", "alb", "hbf", "yrqtr","latlong","mainline","branchline")]
    gdat <- gdat[gdat$mainline %in% c(1,2),]
    gdat <- gdat[gdat$branchline %in% c(1,2),]
    } else {
    gdat <- gdat[,c("vessid","hooks","yft", "bet", "alb", "hbf", "yrqtr","latlong")]
    }
  if(addother) {
      a <- (gdat[,switch(runsp,"bet"="yft","yft"="bet")]+1)/gdat$hooks
      divs <- c(0,0.1, 0.5, 0.9,1)
      b <- quantile(a,divs)
      gdat$other <- findInterval(a,b)
      }
  if(addalb) {
      a <- (gdat[,"alb"]+1)/gdat$hooks
      divs <- c(0,0.1, 0.5, 0.9,1)
      b <- quantile(a,divs)
      gdat$alb_cat <- findInterval(a,b)
      }
  a <- grep(switch(runsp,"bet"="yft","yft"="bet"),names(gdat))
  if(!doboth) gdat <- gdat[,-a]
  a <- grep("alb",names(gdat),fixed=T)[1]
  gdat <- gdat[,-a]
  return(gdat)
  }

select_data_JPIO <- function(indat,runreg,runsp,mt,minqtrs=2,maxqtrs=500,llstrat=5,addcl=NA,addpca=NA,samp=1) {
  gdat <- indat[indat$reg==runreg,]
  if(sum(is.na(gdat$hbf)) > 0) gdat[is.na(gdat$hbf),]$hbf <- 5
  gdat <- gdat[gdat$hbf >= 5,]
  if(llstrat!=5) gdat$latlong <- paste(llstrat*floor(gdat$lat/llstrat),llstrat*floor(gdat$lon/llstrat),sep="_")
  if(mt=="deltapos") gdat <- gdat[gdat[,runsp] > 0,]
  a <- table(gdat$vessid,gdat$yrqtr)
  a <- apply(a>0,1,sum)
  table(a)
  a <- a[a >= minqtrs & a <= maxqtrs]
  gdat <- gdat[gdat$vessid %in% names(a),]
  a <- table(gdat$yrqtr);a
  a <- a[a>=100]
  gdat <- gdat[gdat$yrqtr %in% names(a),]
  a <- table(gdat$vessid);a
  a <- a[a>=100]
  gdat <- gdat[gdat$vessid %in% names(a),]
  vars <- c("vessid","hooks","yft", "bet", "hbf", "yrqtr","latlong", "moon")
  if(!is.na(addcl)) vars <- c(vars,addcl)
  if(!is.na(addpca)) vars <- c(vars,addpca)
  gdat <- gdat[,vars]
  a <- grep(switch(runsp,"bet"="yft","yft"="bet"),names(gdat))
  gdat <- gdat[,-a]
  if(!is.na(samp)) gdat <- samp_data2(gdat,samp)
  gdat$vessid <- as.factor(gdat$vessid)
  gdat$latlong <- as.factor(gdat$latlong)
  gdat$yrqtr <- as.factor(gdat$yrqtr)
  return(gdat)
  }

select_data_KRSBT <- function(indat,runreg,runsp,mt,minqtrs=2,maxqtrs=500,minperyr=100,minpervess=100,minperll=100,llstrat=5,addcl=NA,addpca=NA,samp=1) {
  gdat <- indat
  if(runreg < 15) gdat <- indat[indat$sarea==runreg,]
  if(runreg==789)  gdat <- indat[indat$sarea %in% c(7,8,9),]
  if(sum(is.na(gdat$hbf)) > 0) gdat[is.na(gdat$hbf),]$hbf <- 5
  #gdat <- gdat[gdat$hbf >= 5,]
  if(llstrat!=5) gdat$latlong <- paste(llstrat*floor(gdat$lat/llstrat),llstrat*floor(gdat$lon/llstrat),sep="_")
  if(mt=="deltapos") gdat <- gdat[gdat[,runsp] > 0,]
  a <- table(gdat$vessid,gdat$yrqtr)
  a <- apply(a>0,1,sum)
  table(a)
  a <- a[a >= minqtrs & a <= maxqtrs]
  gdat <- gdat[gdat$vessid %in% names(a),]
  a <- table(gdat$yrqtr);a
  a <- a[a>=minperyr]
  gdat <- gdat[gdat$yrqtr %in% names(a),]
  a <- table(gdat$vessid);a
  a <- a[a>=minpervess]
  gdat <- gdat[gdat$vessid %in% names(a),]
  a <- table(gdat$latlong);a
  a <- a[a>=minperll]
  gdat <- gdat[gdat$latlong %in% names(a),]
  vars <- c("vessid","hooks","sbt", "hbf", "yrqtr", "op_yr","latlong","moon","op_mon")
  if(!is.na(addcl)) vars <- c(vars,addcl)
  if(!is.na(addpca)) vars <- c(vars,addpca)
  gdat <- gdat[,vars]
  if(!is.na(samp)) gdat <- samp_data2(gdat,samp)
  gdat$vessid <- factor(gdat$vessid)
  gdat$latlong <- factor(gdat$latlong)
  gdat$yrqtr <- factor(gdat$yrqtr,levels=seq(min(gdat$yrqtr),max(gdat$yrqtr),0.25))
  gdat$op_yr <- factor(gdat$op_yr,levels=seq(min(gdat$op_yr),max(gdat$op_yr),1))
  return(gdat)
  }

select_data_TWIO <- function(indat,runreg,runsp,mt,minqtrs=2,maxqtrs=500,llstrat=5,addcl=NA,addpca=NA,samp=1) {
  gdat <- indat[indat$reg==runreg,]
  if(sum(is.na(gdat$hbf)) > 0) gdat[is.na(gdat$hbf),]$hbf <- 5
  gdat <- gdat[gdat$hbf >= 5,]
  if(llstrat!=5) gdat$latlong <- paste(llstrat*floor(gdat$lat/llstrat),llstrat*floor(gdat$lon/llstrat),sep="_")
  if(mt=="deltapos") gdat <- gdat[gdat[,runsp] > 0,]
  a <- table(gdat$vessid,gdat$yrqtr)
  a <- apply(a>0,1,sum)
  table(a)
  a <- a[a >= minqtrs & a <= maxqtrs]
  gdat <- gdat[gdat$vessid %in% names(a),]
  a <- table(gdat$yrqtr);a
  a <- a[a>=100]
  gdat <- gdat[gdat$yrqtr %in% names(a),]
  a <- table(gdat$vessid);a
  a <- a[a>=100]
  gdat <- gdat[gdat$vessid %in% names(a),]
  vars <- c("vessid","hooks","yft", "bet", "yrqtr","latlong", "moon","bt1","bt2","bt3","bt4","bt5")
  if(!is.na(addcl)) vars <- c(vars,addcl)
  if(!is.na(addpca)) vars <- c(vars,addpca)
  gdat <- gdat[,vars]
  a <- grep(switch(runsp,"bet"="yft","yft"="bet"),names(gdat))
  gdat <- gdat[,-a]
  if(!is.na(samp)) gdat <- samp_data2(gdat,samp)
  gdat$vessid <- as.factor(gdat$vessid)
  gdat$latlong <- as.factor(gdat$latlong)
  gdat$yrqtr <- as.factor(gdat$yrqtr)
  return(gdat)
  }

select_data_JointIO <- function(indat,runreg,runsp,mt,minqtrs=2,maxqtrs=500,llstrat=5,addcl=NA,addpca=NA,samp=1) {
  gdat <- indat[indat$reg==runreg,]
  if(sum(is.na(gdat$hbf)) > 0) gdat[is.na(gdat$hbf),]$hbf <- 5
  gdat <- gdat[gdat$hbf >= 5,]
  if(llstrat!=5) gdat$latlong <- paste(llstrat*floor(gdat$lat/llstrat),llstrat*floor(gdat$lon/llstrat),sep="_")
  if(mt=="deltapos") gdat <- gdat[gdat[,runsp] > 0,]
  a <- table(gdat$vessid,gdat$yrqtr)
  a <- apply(a>0,1,sum)
  table(a)
  a <- a[a >= minqtrs & a <= maxqtrs]
  gdat <- gdat[gdat$vessid %in% names(a),]
  a <- table(gdat$yrqtr);a
  a <- a[a>=100]
  gdat <- gdat[gdat$yrqtr %in% names(a),]
  a <- table(gdat$vessid);a
  a <- a[a>=100]
  gdat <- gdat[gdat$vessid %in% names(a),]
  vars <- c("vessid","hooks","yft", "bet", "yrqtr","latlong","hbf")
  if(!is.na(addcl)) vars <- c(vars,addcl)
  if(!is.na(addpca)) vars <- c(vars,addpca)
  gdat <- gdat[,vars]
  a <- grep(switch(runsp,"bet"="yft","yft"="bet"),names(gdat))
  gdat <- gdat[,-a]
  if(!is.na(samp)) gdat <- samp_data2(gdat,samp)
  gdat$vessid <- as.factor(gdat$vessid)
  gdat$latlong <- as.factor(gdat$latlong)
  gdat$yrqtr <- as.factor(gdat$yrqtr)
  return(gdat)
  }



select_data_TW <- function(indat,runreg,runsp,mt,minqtrs=2,maxqtrs=500,addmain=F,addother=F,addalb=F,fc="both",bait="no23",llstrat=5,doboth=F) {
  gdat <- indat[indat$reg==runreg,]
  if(llstrat!=5) gdat$latlong <- paste(llstrat*floor(gdat$lat/llstrat),llstrat*floor(gdat$lon/llstrat),sep="_")
  if(mt=="deltapos") gdat <- gdat[gdat[,runsp] > 0,]
  if(bait=="no23") gdat <- gdat[(gdat$bait!=2 & gdat$bait !=3) | is.na(gdat$bait),]
  a <- table(gdat$vessid,gdat$yrqtr)
  a <- apply(a>0,1,sum)
  table(a)
  a <- a[a >= minqtrs & a <= maxqtrs]
  gdat <- gdat[gdat$vessid %in% names(a),]
  a <- table(gdat$yrqtr);a
  a <- a[a>=100]
  gdat <- gdat[gdat$yrqtr %in% names(a),]
  a <- table(gdat$vessid);a
  a <- a[a>=100]
  gdat <- gdat[gdat$vessid %in% names(a),]
  if(sum(is.na(gdat$hbf)) > 0) gdat[is.na(gdat$hbf),]$hbf <- 5
  gdat <- gdat[gdat$hbf >= 5,]
  fcnum <- switch(fc,"both"=0,"OS"=1,"DW"=2)
  if (fc!="both") gdat <- gdat[gdat$newfishingcat==fcnum,]
  if(addmain) {
    gdat <- gdat[gdat$target==3,c("vessid","hooks","yft", "bet", "alb", "hbf", "yrqtr","latlong","mainline","branchline")]
    gdat <- gdat[gdat$mainline %in% c(1,2),]
    gdat <- gdat[gdat$branchline %in% c(1,2),]
    } else {
    gdat <- gdat[,c("vessid","hooks","yft", "bet", "alb", "hbf", "yrqtr","latlong")]
    }
  if(addother) {
      a <- (gdat[,switch(runsp,"bet"="yft","yft"="bet")]+1)/gdat$hooks
      divs <- c(0,0.1, 0.5, 0.9,1)
      b <- quantile(a,divs)
      gdat$other <- findInterval(a,b)
      }
  if(addalb) {
      a <- (gdat[,"alb"]+1)/gdat$hooks
      divs <- c(0,0.1, 0.5, 0.9,1)
      b <- quantile(a,divs)
      gdat$alb_cat <- findInterval(a,b)
      }
  a <- grep(switch(runsp,"bet"="yft","yft"="bet"),names(gdat))
  if(!doboth) gdat <- gdat[,-a]
  a <- grep("alb",names(gdat),fixed=T)[1]
  gdat <- gdat[,-a]
  return(gdat)
  }

qqDist <- function (x, standardise = F, add.median = F, ...)
{
    n <- length(x)
    seq.length <- min(1000, n)
    if (standardise) {
        SEQ <- seq(1, 2 * n + 1, length = seq.length)/2
        U <- qnorm(qbeta(0.975, SEQ, rev(SEQ)))
        L <- qnorm(qbeta(0.025, SEQ, rev(SEQ)))
        if (add.median)
            M <- qnorm(qbeta(0.5, SEQ, rev(SEQ)))
    }
    else {
        SD <- sqrt(var(x) * (n + 1)/n)
        SEQ <- seq(1, 2 * n + 1, length = seq.length)/2
        U <- mean(x) + SD * qt(qbeta(0.975, SEQ, rev(SEQ)), n -
            1)
        L <- mean(x) + SD * qt(qbeta(0.025, SEQ, rev(SEQ)), n -
            1)
        if (add.median)
            M <- mean(x) + SD * qt(qbeta(0.5, SEQ, rev(SEQ)),
                n - 1)
    }
    X <- qnorm((SEQ - 0.25)/(n + 0.5))
    qqnorm(x, main = "", ...)
    lines(X, U, type = "l",col=2)
    lines(X, L, type = "l",col=2)
    if (add.median)
        lines(X, M, type = "l",col=2)
    invisible()
}

plotdiags <- function(res,ti="") {
  hist(res,nclass=200,freq=F,xlab="Residuals",main=ti)
  lines((-30:30)/10,dnorm((-30:30)/10,sd=sd(res)),col=2)
  sdres <- res/sd(res)
  qqDist(sdres,add.median=T)
  }

splitvessels <- function(indat,period) {
  vess <- unique(indat$vessid)
  indat$oldvess <- indat$vessid
  indat$vessid <- ""
  minyr <- maxyr <- vess
  for (i in 1:length(vess)) {
    a <- grep(vess[i],indat$oldvess)
    minyr <- min(indat[a,]$yrqtr)
    indat[a,]$vessid <- paste(indat[a,]$oldvess,floor((indat[a,]$yrqtr - minyr)/period))
    }
  return(indat)
  }

mt <- "deltabin"
mt <- "deltapos"

plot_effects <- function(model,indat,addmain=F,addbranch=F,addalb=F,addother=F,ti="") {
  cf <- model$coefficients
  pred <- predict(model,data=indat,type="terms",se.fit=T)
  fishlab <- switch(runsp,yft="Yellowfin",bet="Bigeye"); methlab <- switch(mt,deltabin="Delta-binomial",deltapos="Delta-positive",logl="Lognormal(+0.01)",propn="Proportion Bigeye")
  nfigs <- 6 + addmain + addbranch + addalb + addother
  mf <- switch(nfigs-5,c(2,3),c(3,3),c(3,3),c(3,3),c(3,4))
  hw <- c(14,19)
  windows(height=hw[1],width=hw[2])
  par(mfrow=mf,mar=c(5,4,2,1),oma=c(0,0,3,0))
  pr <- pred$fit ; prse <- pred$se.fit
  termlist <- dimnames(pred$fit)[[2]]
  llpos <- grep("latlong",termlist)
  hbfpos  <- grep("hbf",termlist)
  mainpos <- grep("mainline",termlist)
  vesspos <- grep("vessid",termlist)
  branchpos <- grep("branchline",termlist)
  if(length(grep("hooks",termlist))>0) db <- T else db <- F

  index <- sort(unique(indat$yrqtr))
  b <- match(index,indat$yrqtr)
  out <- exp(pr[b,1])
  se1 <- exp(pr[b,1] - 1.96 * prse[b,1])
  se2 <- exp(pr[b,1] + 1.96 * prse[b,1])
  plot(as.numeric(as.character(index)),out,ylim=c(0,3),xlab="Year",ylab="Effect size",main="Year effects")
  segments(as.numeric(as.character(index)), se1,  as.numeric(as.character(index)), se2, lty=1, col="slate grey")
  points(as.numeric(as.character(index)), out, pch=16)

  index <- sort(unique(indat$latlong))
  b <- match(index,indat$latlong)
  out <- exp(pr[b,llpos])
  se1 <- exp(pr[b,llpos] - 1.96 * prse[b,llpos])
  se2 <- exp(pr[b,llpos] + 1.96 * prse[b,llpos])
  plot(as.numeric(as.character(index)),out,ylim=c(0,3),xlab="Latlong",ylab="Effect size",main="Spatial effects")
  segments(as.numeric(as.character(index)), se1,  as.numeric(as.character(index)), se2, lty=1, col="slate grey")
  points(as.numeric(as.character(index)), out, pch=16)
#  plot(indat$latlong,pr[,2],ylim=c(-.75,.75),xlab="Latitude",ylab="Effect size")

  ##plot coefficients
  library(maps)
  ll <- as.numeric(indat$latlong)
  index <- sort(unique(ll))
  lats <- trunc(ll)
  lons <- 1000*abs((ll - trunc(ll)))
  alat <- round(lats[match(index, ll)]*2)/2 + 2.5
  alon <- round(lons[match(index, ll)]*2)/2 + 2.5
# lat <- trunc((dat$newlat[match(index, dat$latlong)] + 50)/5) * 5 + 2.5 - 50
# long <- trunc((dat$newlong[match(index, dat$latlong)])/5) * 5 + 2.5
  coefs <- exp(pr[match(index, ll),llpos])
  coefs2 <- tapply(coefs, list(alon, alat), mean)
  image(sort(as.numeric(unique(alon))), sort(unique(alat)), coefs2, zlim=c(0.5, 2.5), ylab="Lat", xlab="Long",col=rev(heat.colors(12)))
  contour(sort(unique(alon)), sort(unique(alat)), coefs2, levels= c(0,1,2,2.5), add=TRUE,col=4)
  maps::map(add=TRUE)

  index <- sort(unique(indat$vessid))
  b <- match(index,indat$vessid)
  out <- exp(pr[b,vesspos])
  se1 <- exp(pr[b,vesspos] - 1.96 * prse[b,vesspos])
  se2 <- exp(pr[b,vesspos] + 1.96 * prse[b,vesspos])
  plot(as.factor(index),out,type="n",ylim=c(0,3),xlab="Vessel ID",ylab="Effect size",main="")
  segments(match(index,index), se1,  match(index,index), se2, lty=1, col="slate grey")
  points(as.factor(index), out, pch=16)

  vts <- tapply(pr[,vesspos],list(indat$yrqtr,indat$vessid),mean)
  y<-exp(as.vector(vts))
  x<-rep(as.numeric(dimnames(vts)[[1]]),dim(vts)[[2]])
  plot(x,y,type="p",ylim=c(0,3),xlab="Vessel ID",ylab="Effect size",main="Vessel effects",pch=16,cex=0.9)
  mn <- tapply(exp(pr[,vesspos]),indat$yrqtr,mean)
  lines(as.numeric(dimnames(vts)[[1]]), mn, col=2,lwd=2)
#  plot(as.factor(indat$vessid),pr[,vesspos],ylim=c(-.75,.75),xlab="Vessel",ylab="Effect size")

  rsp <- pr[,hbfpos[1]]
  rsp.se <- prse[,hbfpos[1]]
  if(addmain) {
    a <- indat$mainline==1
    rsp2 <- pr[,hbfpos[1]] + pr[,hbfpos[2]] + pr[,mainpos[1]]
    rsp2.se <- exp(prse[,hbfpos[1]] + prse[,hbfpos[2]] + prse[,mainpos[1]])
    b <- match(sort(unique(indat[a,]$hbf)),indat[a,]$hbf)
    plot(indat[a,]$hbf[b],exp(rsp[a][b]),ylim=c(0,10),xlab="Nylon mainline HBF",ylab="Effect size",main=paste("Nylon HBF"))
    lines(indat[a,]$hbf[b],exp(rsp[a][b]+1.96*rsp.se[a][b]),lty=2)
    lines(indat[a,]$hbf[b],exp(rsp[a][b]-1.96*rsp.se[a][b]),lty=2)
    b <- match(sort(unique(indat[a==F,]$hbf)),indat[a==F,]$hbf)
    plot(indat[a==F,]$hbf[b],exp(rsp2[a==F][b]),ylim=c(0,3),xlab="Other mainline HBF",ylab="Effect size",main = "Other HBF")
    lines(indat[a==F,]$hbf[b],exp(rsp2[a==F][b] + 1.96*rsp2.se[a==F][b]),lty=2)
    lines(indat[a==F,]$hbf[b],exp(rsp2[a==F][b] - 1.96*rsp2.se[a==F][b]),lty=2)
    } else {
    b <- match(sort(unique(indat$hbf)),indat$hbf)
    plot(indat$hbf[b],exp(pr[,hbfpos][b]),ylim=c(0,3),xlab="HBF",ylab="Effect size",main="HBF effects")
    lines(indat$hbf[b],exp(rsp[b]+2*rsp.se[b]),lty=2)
    lines(indat$hbf[b],exp(rsp[b]-2*rsp.se[b]),lty=2)
    }
  }

plot_effects_IO <- function(model,indat,addcl=F,addpca=F,ti="",dohbf=T) {
  cf <- model$coefficients
  pred <- predict(model,data=indat,type="terms",se.fit=T)
  fishlab <- switch(runsp,yft="Yellowfin",bet="Bigeye"); methlab <- switch(mt,deltabin="Delta-binomial",deltapos="Delta-positive",logl="Lognormal(+0.01)",propn="Proportion Bigeye")
  nfigs <- 6 + addcl + addpca
  mf <- switch(nfigs-5,c(2,3),c(3,3),c(3,3),c(3,3),c(3,4))
  hw <- c(14,19)
  windows(height=hw[1],width=hw[2])
  par(mfrow=mf,mar=c(5,4,2,1),oma=c(0,0,3,0))
  pr <- pred$fit ; prse <- pred$se.fit
  termlist <- dimnames(pred$fit)[[2]]
  llpos <- grep("latlong",termlist)
  hbfpos  <- grep("hbf",termlist)
  vesspos <- grep("vessid",termlist)
  branchpos <- grep("branchline",termlist)
  if(length(grep("hooks",termlist))>0) db <- T else db <- F

  index <- sort(unique(indat$yrqtr))
  b <- match(index,indat$yrqtr)
  out <- exp(pr[b,1])
  se1 <- exp(pr[b,1] - 1.96 * prse[b,1])
  se2 <- exp(pr[b,1] + 1.96 * prse[b,1])
  plot(as.numeric(as.character(index)),out,ylim=c(0,3),xlab="Year",ylab="Effect size",main="Year effects")
  segments(as.numeric(as.character(index)), se1,  as.numeric(as.character(index)), se2, lty=1, col="slate grey")
  points(as.numeric(as.character(index)), out, pch=16)

  index <- sort(unique(indat$latlong))
  b <- match(index,indat$latlong)
  out <- exp(pr[b,llpos])
  se1 <- exp(pr[b,llpos] - 1.96 * prse[b,llpos])
  se2 <- exp(pr[b,llpos] + 1.96 * prse[b,llpos])
  plot(as.numeric(index),out,ylim=c(0,3),xlab="Latlong",ylab="Effect size",main="Spatial effects")
  segments(as.numeric(index), se1,  as.numeric(index), se2, lty=1, col="slate grey")
  points(as.numeric(index), out, pch=16)
#  plot(indat$latlong,pr[,2],ylim=c(-.75,.75),xlab="Latitude",ylab="Effect size")

  ##plot coefficients
  library(maps)
  ll <- as.character(indat$latlong)
  index <- sort(unique(ll))
  index2 <- unlist(strsplit(index,"\\_"))
  pos1 <- 2*(1:length(index)) -1
  pos2 <- pos1 + 1
  lats <- as.numeric(index2[pos1])
  lons <- as.numeric(index2[pos2])
  alat <- round(lats[match(ll,index)]*2)/2 + 2.5
  alon <- round(lons[match(ll,index)]*2)/2 + 2.5
# lat <- trunc((dat$newlat[match(index, dat$latlong)] + 50)/5) * 5 + 2.5 - 50
# long <- trunc((dat$newlong[match(index, dat$latlong)])/5) * 5 + 2.5
  coefs <- exp(pr[match(ll,index),llpos])
  coefs2 <- tapply(coefs, list(alon, alat), mean)
  image(sort(as.numeric(unique(alon))), sort(unique(alat)), coefs2, zlim=c(0.5, 2.5), ylab="Lat", xlab="Long",col=rev(heat.colors(12)))
  contour(sort(unique(alon)), sort(unique(alat)), coefs2, levels= c(0,1,2,2.5), add=TRUE,col=4)
  maps::map(add=TRUE)

  if("vessid" %in% termlist) {
    index <- sort(unique(indat$vessid))
    b <- match(index,indat$vessid)
    out <- exp(pr[b,vesspos])
    se1 <- exp(pr[b,vesspos] - 1.96 * prse[b,vesspos])
    se2 <- exp(pr[b,vesspos] + 1.96 * prse[b,vesspos])
    ind <- as.factor(index)
    plot(match(ind,ind),out,type="n",ylim=c(0,3),xlab="Vessel ID",ylab="Effect size",main="")
    segments(match(ind,ind), se1,  match(ind,ind), se2, lty=1, col="slate grey")
    points(match(ind,ind), out, pch=16)
  
    vts <- tapply(pr[,vesspos],list(indat$yrqtr,indat$vessid),mean)
    y<-exp(as.vector(vts))
    x<-rep(as.numeric(dimnames(vts)[[1]]),dim(vts)[[2]])
    plot(x,y,type="p",ylim=c(0,3),xlab="Vessel ID",ylab="Effect size",main="Vessel effects",pch=16,cex=0.9)
    mn <- tapply(exp(pr[,vesspos]),indat$yrqtr,mean)
    lines(as.numeric(dimnames(vts)[[1]]), mn, col=2,lwd=2)
  #  plot(as.factor(indat$vessid),pr[,vesspos],ylim=c(-.75,.75),xlab="Vessel",ylab="Effect size")
  }

  if(dohbf) {
    rsp <- pr[,hbfpos[1]]
    rsp.se <- prse[,hbfpos[1]]
    b <- match(sort(unique(indat$hbf)),indat$hbf)
    plot(indat$hbf[b],exp(pr[,hbfpos][b]),ylim=c(0,3),xlab="HBF",ylab="Effect size",main="HBF effects")
    lines(indat$hbf[b],exp(rsp[b]+2*rsp.se[b]),lty=2)
    lines(indat$hbf[b],exp(rsp[b]-2*rsp.se[b]),lty=2)
    }

  if(addcl) {
    index <- sort(unique(indat$cl))
    b <- match(index,indat$cl)
    out <- exp(pr[b,1])
    se1 <- exp(pr[b,1] - 1.96 * prse[b,1])
    se2 <- exp(pr[b,1] + 1.96 * prse[b,1])
    plot(as.numeric(as.character(index)),out,ylim=c(0,3),xlab="Cluster category",ylab="Effect size",main="Cluster effects")
    segments(as.numeric(as.character(index)), se1,  as.numeric(as.character(index)), se2, lty=1, col="slate grey")
    points(as.numeric(as.character(index)), out, pch=16)
    }

  if(addpca) {
    ind <- indat[order(indat$pca),]
    out <- exp(pr[b,1])
    se1 <- exp(pr[b,1] - 1.96 * prse[b,1])
    se2 <- exp(pr[b,1] + 1.96 * prse[b,1])
    plot(as.numeric(as.character(index)),out,ylim=c(0,3),xlab="Cluster category",ylab="Effect size",main="Cluster effects")
    segments(as.numeric(as.character(index)), se1,  as.numeric(as.character(index)), se2, lty=1, col="slate grey")
    points(as.numeric(as.character(index)), out, pch=16)
    }
  title(main=ti,cex.main=1.5,outer=T)
  }

plot_effects_SBT <- function(model,indat,addcl=F,addpca=F,ti="",dohbf=T,domonth=domonth) {
  cf <- model$coefficients
  pred <- predict(model,data=indat,type="terms",se.fit=T)
  fishlab <- "SBT"; methlab <- switch(mt,deltabin="Delta-binomial",deltapos="Delta-positive",logl="Lognormal(+0.01)",propn="Proportion Bigeye")
  nfigs <- 5 + addcl + addpca + dohbf
  mf <- switch(nfigs-4,c(2,3),c(2,3),c(3,3),c(3,3),c(3,3),c(3,4))
  hw <- c(14,19)
  windows(height=hw[1],width=hw[2])
  par(mfrow=mf,mar=c(5,4,2,1),oma=c(0,0,3,0))
  pr <- pred$fit ; prse <- pred$se.fit
  termlist <- dimnames(pred$fit)[[2]]
  llpos <- grep("latlong",termlist)
  hbfpos  <- grep("hbf",termlist)
  vesspos <- grep("vessid",termlist)
  branchpos <- grep("branchline",termlist)
  monthpos <- grep("op_mon",termlist)
  if(length(grep("hooks",termlist))>0) db <- T else db <- F

  index <- sort(unique(indat$yrqtr))
  b <- match(index,indat$yrqtr)
  out <- exp(pr[b,1])
  se1 <- exp(pr[b,1] - 1.96 * prse[b,1])
  se2 <- exp(pr[b,1] + 1.96 * prse[b,1])
  plot(as.numeric(as.character(index)),out,ylim=c(0,3),xlab="Year",ylab="Effect size",main="Year effects")
  segments(as.numeric(as.character(index)), se1,  as.numeric(as.character(index)), se2, lty=1, col="slate grey")
  points(as.numeric(as.character(index)), out, pch=16)

  index <- sort(unique(indat$latlong))
  b <- match(index,indat$latlong)
  out <- exp(pr[b,llpos])
  se1 <- exp(pr[b,llpos] - 1.96 * prse[b,llpos])
  se2 <- exp(pr[b,llpos] + 1.96 * prse[b,llpos])
  plot(as.numeric(index),out,ylim=c(0,3),xlab="Latlong",ylab="Effect size",main="Spatial effects")
  segments(as.numeric(index), se1,  as.numeric(index), se2, lty=1, col="slate grey")
  points(as.numeric(index), out, pch=16)

  ##plot coefficients
  library(maps)
  ll <- as.character(indat$latlong)
  index <- sort(unique(ll))
  index2 <- unlist(strsplit(index,"\\_"))
  pos1 <- 2*(1:length(index)) -1
  pos2 <- pos1 + 1
  lats <- as.numeric(index2[pos1])
  lons <- as.numeric(index2[pos2])
  lons[lons > 190] <- lons[lons > 190] - 360
  alat <- round(lats[match(ll,index)]*2)/2 + 2.5
  alon <- round(lons[match(ll,index)]*2)/2 + 2.5
  alat <- factoralong(alat,5)
  alon <- factoralong(alon,5)
  coefs <- exp(pr[match(ll,index),llpos])
  coefs2 <- tapply(coefs, list(alon, alat), mean)
  edge1 <- as.numeric(dimnames(coefs2)[[1]])
  edge2 <- as.numeric(dimnames(coefs2)[[2]])
  image(edge1,edge2,coefs2, zlim=c(0.5, 2.5), ylab="Lat", xlab="Long",col=rev(heat.colors(12)))
  contour(edge1,edge2,coefs2, levels= c(0,1,2,2.5), add=TRUE,col=4)
  maps::map(add=TRUE,fill=T)

  index <- sort(unique(indat$vessid))
  b <- match(index,indat$vessid)
  out <- exp(pr[b,vesspos])
  se1 <- exp(pr[b,vesspos] - 1.96 * prse[b,vesspos])
  se2 <- exp(pr[b,vesspos] + 1.96 * prse[b,vesspos])
  ind <- as.factor(index)
  plot(match(ind,ind),out,type="n",ylim=c(0,3),xlab="Vessel ID",ylab="Effect size",main="")
  segments(match(ind,ind), se1,  match(ind,ind), se2, lty=1, col="slate grey")
  points(match(ind,ind), out, pch=16)

  vts <- tapply(pr[,vesspos],list(indat$yrqtr,indat$vessid),mean)
  y<-exp(as.vector(vts))
  x<-rep(as.numeric(dimnames(vts)[[1]]),dim(vts)[[2]])
  plot(x,y,type="p",ylim=c(0,3),xlab="Vessel ID",ylab="Effect size",main="Vessel effects",pch=16,cex=0.9)
  mn <- tapply(exp(pr[,vesspos]),indat$yrqtr,mean)
  lines(as.numeric(dimnames(vts)[[1]]), mn, col=2,lwd=2)

  if(dohbf) {
    rsp <- pr[,hbfpos[1]]
    rsp.se <- prse[,hbfpos[1]]
    b <- match(sort(unique(indat$hbf)),indat$hbf)
    plot(indat$hbf[b],exp(pr[,hbfpos][b]),ylim=c(0,3),xlab="HBF",ylab="Effect size",main="HBF effects")
    lines(indat$hbf[b],exp(rsp[b]+2*rsp.se[b]),lty=2)
    lines(indat$hbf[b],exp(rsp[b]-2*rsp.se[b]),lty=2)
    }

  if(addcl) {
    index <- sort(unique(indat$cl))
    b <- match(index,indat$cl)
    out <- exp(pr[b,1])
    se1 <- exp(pr[b,1] - 1.96 * prse[b,1])
    se2 <- exp(pr[b,1] + 1.96 * prse[b,1])
    plot(as.numeric(as.character(index)),out,ylim=c(0,3),xlab="Cluster category",ylab="Effect size",main="Cluster effects")
    segments(as.numeric(as.character(index)), se1,  as.numeric(as.character(index)), se2, lty=1, col="slate grey")
    points(as.numeric(as.character(index)), out, pch=16)
    }

  if(addpca) {
    ind <- indat[order(indat$pca),]
    out <- exp(pr[b,1])
    se1 <- exp(pr[b,1] - 1.96 * prse[b,1])
    se2 <- exp(pr[b,1] + 1.96 * prse[b,1])
    plot(as.numeric(as.character(index)),out,ylim=c(0,3),xlab="Cluster category",ylab="Effect size",main="Cluster effects")
    segments(as.numeric(as.character(index)), se1,  as.numeric(as.character(index)), se2, lty=1, col="slate grey")
    points(as.numeric(as.character(index)), out, pch=16)
    }

  if(domonth) {
    index <- sort(unique(indat$op_mon))
    b <- match(index,indat$op_mon)
    out <- exp(pr[b,monthpos])
    se1 <- exp(pr[b,monthpos] - 1.96 * prse[b,monthpos])
    se2 <- exp(pr[b,monthpos] + 1.96 * prse[b,monthpos])
    ind <- as.factor(index)
    plot(match(ind,ind),out,type="n",ylim=c(0,3),xlab="Month",ylab="Effect size",main="")
    segments(match(ind,ind), se1,  match(ind,ind), se2, lty=1, col="slate grey")
    points(match(ind,ind), out, pch=16)
    }
  title(main=ti,cex.main=1.5,outer=T)
  }

# plot_effects_SBT_yr <- function(model,indat,addcl=F,addpca=F,ti="",dohbf=T,domonth=domonth,domoon=F) {
#   cf <- model$coefficients
#   pred <- predict(model,data=indat,type="terms",se.fit=T)
#   fishlab <- "SBT"; methlab <- switch(mt,deltabin="Delta-binomial",deltapos="Delta-positive",logl="Lognormal(+0.01)",propn="Proportion Bigeye")
#   nfigs <- 5 + addcl + addpca + dohbf + domonth + domoon
#   mf <- switch(nfigs-4,c(2,3),c(2,3),c(3,3),c(3,3),c(3,3),c(3,4))
#   hw <- c(14,19)
#   windows(height=hw[1],width=hw[2])
#   par(mfrow=mf,mar=c(5,4,2,1),oma=c(0,0,3,0))
#   pr <- pred$fit ; prse <- pred$se.fit
#   termlist <- dimnames(pred$fit)[[2]]
#   llpos <- grep("latlong",termlist)
#   hbfpos  <- grep("hbf",termlist)
#   vesspos <- grep("vessid",termlist)
#   branchpos <- grep("branchline",termlist)
#   monthpos <- grep("op_mon",termlist)
#   moonpos <- grep("moon",termlist)
#   if(length(grep("hooks",termlist))>0) db <- T else db <- F
# 
#   index <- sort(unique(indat$op_yr))
#   b <- match(index,indat$op_yr)
#   out <- exp(pr[b,1])
#   se1 <- exp(pr[b,1] - 1.96 * prse[b,1])
#   se2 <- exp(pr[b,1] + 1.96 * prse[b,1])
#   plot(as.numeric(as.character(index)),out,ylim=c(0,4),xlab="Year",ylab="Effect size",main="Year effects")
#   segments(as.numeric(as.character(index)), se1,  as.numeric(as.character(index)), se2, lty=1, col="slate grey")
#   points(as.numeric(as.character(index)), out, pch=16)
# 
#   index <- sort(unique(indat$latlong))
#   b <- match(index,indat$latlong)
#   out <- exp(pr[b,llpos])
#   se1 <- exp(pr[b,llpos] - 1.96 * prse[b,llpos])
#   se2 <- exp(pr[b,llpos] + 1.96 * prse[b,llpos])
#   plot(as.numeric(index),out,ylim=c(0,4),xlab="Latlong",ylab="Effect size",main="Spatial effects")
#   segments(as.numeric(index), se1,  as.numeric(index), se2, lty=1, col="slate grey")
#   points(as.numeric(index), out, pch=16)
# 
#   ##plot coefficients
#   library(maps)
#   ll <- as.character(indat$latlong)
#   index <- sort(unique(ll))
#   index2 <- unlist(strsplit(index,"\\_"))
#   pos1 <- 2*(1:length(index)) -1
#   pos2 <- pos1 + 1
#   lats <- as.numeric(index2[pos1])
#   lons <- as.numeric(index2[pos2])
#   lons[lons > 190] <- lons[lons > 190] - 360
#   alat <- round(lats[match(ll,index)]*2)/2 + 2.5
#   alon <- round(lons[match(ll,index)]*2)/2 + 2.5
#   alat <- factoralong(alat,5)
#   alon <- factoralong(alon,5)
#   coefs <- exp(pr[match(ll,index),llpos])
#   coefs2 <- tapply(coefs, list(alon, alat), mean)
#   edge1 <- as.numeric(dimnames(coefs2)[[1]])
#   edge2 <- as.numeric(dimnames(coefs2)[[2]])
#   image(edge1,edge2,coefs2, zlim=c(0.5, 2.5), ylab="Lat", xlab="Long",col=rev(heat.colors(12)))
#   contour(edge1,edge2,coefs2, levels= c(0,1,2,2.5), add=TRUE,col=4)
#   maps::map(add=TRUE,fill=T)
# 
#   index <- sort(unique(indat$vessid))
#   b <- match(index,indat$vessid)
#   out <- exp(pr[b,vesspos])
#   se1 <- exp(pr[b,vesspos] - 1.96 * prse[b,vesspos])
#   se2 <- exp(pr[b,vesspos] + 1.96 * prse[b,vesspos])
#   ind <- as.factor(index)
#   plot(match(ind,ind),out,type="n",ylim=c(0,4),xlab="Vessel ID",ylab="Effect size",main="")
#   segments(match(ind,ind), se1,  match(ind,ind), se2, lty=1, col="slate grey")
#   points(match(ind,ind), out, pch=16)
# 
#   vts <- tapply(pr[,vesspos],list(indat$yrqtr,indat$vessid),mean)
#   y<-exp(as.vector(vts))
#   x<-rep(as.numeric(dimnames(vts)[[1]]),dim(vts)[[2]])
#   plot(x,y,type="p",ylim=c(0,3),xlab="Vessel ID",ylab="Effect size",main="Vessel effects",pch=16,cex=0.9)
#   mn <- tapply(exp(pr[,vesspos]),indat$yrqtr,mean)
#   lines(as.numeric(dimnames(vts)[[1]]), mn, col=2,lwd=2)
# 
#   if(dohbf) {
#     rsp <- pr[,hbfpos[1]]
#     rsp.se <- prse[,hbfpos[1]]
#     b <- match(sort(unique(indat$hbf)),indat$hbf)
#     plot(indat$hbf[b],exp(pr[,hbfpos][b]),ylim=c(0,4),xlab="HBF",ylab="Effect size",main="HBF effects")
#     lines(indat$hbf[b],exp(rsp[b]+2*rsp.se[b]),lty=2)
#     lines(indat$hbf[b],exp(rsp[b]-2*rsp.se[b]),lty=2)
#     }
# 
#   if(addcl) {
#     index <- sort(unique(indat$cl))
#     b <- match(index,indat$cl)
#     out <- exp(pr[b,1])
#     se1 <- exp(pr[b,1] - 1.96 * prse[b,1])
#     se2 <- exp(pr[b,1] + 1.96 * prse[b,1])
#     plot(as.numeric(as.character(index)),out,ylim=c(0,4),xlab="Cluster category",ylab="Effect size",main="Cluster effects")
#     segments(as.numeric(as.character(index)), se1,  as.numeric(as.character(index)), se2, lty=1, col="slate grey")
#     points(as.numeric(as.character(index)), out, pch=16)
#     }
# 
#   if(addpca) {
#     ind <- indat[order(indat$pca),]
#     out <- exp(pr[b,1])
#     se1 <- exp(pr[b,1] - 1.96 * prse[b,1])
#     se2 <- exp(pr[b,1] + 1.96 * prse[b,1])
#     plot(as.numeric(as.character(index)),out,ylim=c(0,4),xlab="Cluster category",ylab="Effect size",main="Cluster effects")
#     segments(as.numeric(as.character(index)), se1,  as.numeric(as.character(index)), se2, lty=1, col="slate grey")
#     points(as.numeric(as.character(index)), out, pch=16)
#     }
# 
#   if(domonth) {
#     index <- sort(unique(indat$op_mon))
#     b <- match(index,indat$op_mon)
#     out <- exp(pr[b,monthpos])
#     se1 <- exp(pr[b,monthpos] - 1.96 * prse[b,monthpos])
#     se2 <- exp(pr[b,monthpos] + 1.96 * prse[b,monthpos])
#     ind <- index
#     plot(ind,out,type="n",ylim=c(0,3),xlab="Month",ylab="Effect size",main="")
#     segments(ind, se1,  ind, se2, lty=1, col="slate grey")
#     points(ind, out, pch=16)
#     }
#   if(domoon) {
#     index <- sort(unique(indat$moon))
#     b <- match(index,indat$moon)
#     out <- exp(pr[b,moonpos])
#     se1 <- exp(pr[b,moonpos] - 1.96 * prse[b,moonpos])
#     se2 <- exp(pr[b,moonpos] + 1.96 * prse[b,moonpos])
#     ind <- index
#     plot(ind,out,type="n",ylim=c(0,3),xlab="Lunar illumination",ylab="Effect size",main="")
#     lines(ind, se1,  lty=2, col="slate grey")
#     lines(ind, se2, lty=2, col="slate grey")
#     points(ind, out, pch=16)
#     }
#   title(main=ti,cex.main=1.5,outer=T)
#   }
# 
plot.pacific <- function(plot_title="",lims=c(100,260,-45,45)) {
  plot(1,1, yaxt="n", xaxt="n", type="n", xlim=c(lims[1]+10,lims[2]-10), ylim=c(lims[3]+5,lims[4]-5), ylab="", xlab="", bg="lightblue")
  polygon(c(lims[1],lims[2],lims[2],lims[1]), c(lims[3],lims[3],lims[4],lims[4]), col="lightblue")
#  polygon(eez[,1], eez[,2], lwd=1, col="white")
#  lines(eez[,1], eez[,2], lwd=1, col="slate grey")
  maps::map('world2',  yaxt="n", xaxt="n", add=T, resolution=1)
  maps::map('world2',  region = c("USA","Hawaii","Mexico","Japan","China","South Korea","North Korea","Philippines","Vietnam","Laos","Taiwan","Fiji", "Vanuatu", "Malaysia", "Australia", "New Zealand", "Indonesia", "New Caledonia", "Papua New Guinea", "Solomon Islands"), fill=T, add=T, yaxt="n", xaxt="n", col="black", density=50)
  box(lwd=3)
  lines(c(210, 210, 230, 230), c(45, -2.5, -2.5, -45), lwd=2, lty=2)
  lines(c(170, 170), c(-35, 40), lwd=2, lty=1)
  lines(c(120, 210), c(20, 20), lwd=2, lty=1)
  lines(c(120, 230), c(-10, -10), lwd=2, lty=1)
  lines(c(120, 230), c(-35, -35), lwd=2, lty=1)
  axis(1, at=seq(lims[1],lims[2],by=10), labels=F)
  axis(2, at=seq(lims[3],lims[4],by=5), labels=F)
  latseq <- seq(lims[3]+10,lims[4]-10,by=10) ;latseq2 <- as.character(latseq)
  lonseq <- seq(lims[1]+20,lims[2]-20,by=20) ;lonseq2 <- as.character(lonseq)
  latseq2[latseq < 0] <- paste(abs(latseq[latseq < 0]),"S",sep="")
  latseq2[latseq > 0] <- paste(latseq[latseq > 0],"N",sep="")
  lonseq2[lonseq < 180] <- paste(lonseq2[lonseq < 180],"E",sep="")
  lonseq2[lonseq > 180] <- paste(360-lonseq[lonseq > 180],"W",sep="")
  axis(2, at=latseq, labels=latseq2, cex.axis=0.75)
  axis(1, at=lonseq, labels=lonseq2, cex.axis=0.75)
  mtext(side=3, line=0.5, plot_title)
}

plot.pacific <- function(plot_title="",lims=c(100,260,-45,45)) {
  plot(1,1, yaxt="n", xaxt="n", type="n", xlim=c(lims[1]+10,lims[2]-10), ylim=c(lims[3]+5,lims[4]-5), ylab="", xlab="", bg="lightblue")
  polygon(c(lims[1],lims[2],lims[2],lims[1]), c(lims[3],lims[3],lims[4],lims[4]), col="lightblue")
#  polygon(eez[,1], eez[,2], lwd=1, col="white")
#  lines(eez[,1], eez[,2], lwd=1, col="slate grey")
  maps::map('world2',  yaxt="n", xaxt="n", add=T, resolution=1)
  maps::map('world2',  region = c("USA","Hawaii","Mexico","Japan","China","South Korea","North Korea","Philippines","Vietnam","Laos","Taiwan","Fiji", "Vanuatu", "Malaysia", "Australia", "New Zealand", "Indonesia", "New Caledonia", "Papua New Guinea", "Solomon Islands"), fill=T, add=T, yaxt="n", xaxt="n", col="black", density=50)
  box(lwd=3)
  lines(c(210, 210, 230, 230), c(45, -2.5, -2.5, -45), lwd=2, lty=2)
  lines(c(170, 170), c(-35, 40), lwd=2, lty=1)
  lines(c(120, 210), c(20, 20), lwd=2, lty=1)
  lines(c(120, 230), c(-10, -10), lwd=2, lty=1)
  lines(c(120, 230), c(-35, -35), lwd=2, lty=1)
  axis(1, at=seq(lims[1],lims[2],by=10), labels=F)
  axis(2, at=seq(lims[3],lims[4],by=5), labels=F)
  latseq <- seq(lims[3]+10,lims[4]-10,by=10) ;latseq2 <- as.character(latseq)
  lonseq <- seq(lims[1]+20,lims[2]-20,by=20) ;lonseq2 <- as.character(lonseq)
  latseq2[latseq < 0] <- paste(abs(latseq[latseq < 0]),"S",sep="")
  latseq2[latseq > 0] <- paste(latseq[latseq > 0],"N",sep="")
  lonseq2[lonseq < 180] <- paste(lonseq2[lonseq < 180],"E",sep="")
  lonseq2[lonseq > 180] <- paste(360-lonseq[lonseq > 180],"W",sep="")
  axis(2, at=latseq, labels=latseq2, cex.axis=0.75)
  axis(1, at=lonseq, labels=lonseq2, cex.axis=0.75)
  mtext(side=3, line=0.5, plot_title)
}

plot_IO <- function (plot_title = "", uselims = c(20, 130, -50, 25), sp = "YFT")
{
    lims <- uselims
    plot(1, 1, yaxt = "n", xaxt = "n", type = "n", xlim = c(lims[1], lims[2]), ylim = c(lims[3], lims[4]),
        ylab = "", xlab = "", bg = "lightblue")
    polygon(c(lims[1]-5, lims[2]+5, lims[2]+5, lims[1]-5), c(lims[3]-5, lims[3]-5, lims[4]+5, lims[4]+5), col = "lightblue")
    if (sp == "ALB") {
        lines(c(34.5, 44.2), c(-20, -20), lwd = 3, col = "slate grey",lty = 1)
        lines(c(49, 120), c(-20, -20), lwd = 3, col = "slate grey",lty = 1)
        xoffset <- 5
        yoffset <- 2.5
        text(115 + xoffset, -11 - yoffset, "N", col = "red",  cex = 1.5)
        text(115 + xoffset, -40 - yoffset, "S", col = "red",  cex = 1.5)
    }
    if (sp %in% c("YFT")) {
        lines(c(20, 120), c(-40, -40), lwd = 3, col = "slate grey",lty = 1)
        lines(c(20, 20), c(-40, -35), lwd = 3, col = "slate grey",lty = 1)
        lines(c(40, 40), c(-40, -30), lwd = 3, col = "slate grey",lty = 1)
        lines(c(40, 40), c(-10, 10), lwd = 3, col = "slate grey",lty = 1)
        lines(c(60, 60), c(-30, -10), lwd = 3, col = "slate grey",lty = 1)
        lines(c(75, 75), c(-15, 10), lwd = 3, col = "slate grey",lty = 1)
        lines(c(100, 100), c(-5, 10), lwd = 3, col = "slate grey",lty = 1)
        lines(c(110, 110), c(-10, -5), lwd = 3, col = "slate grey",lty = 1)
        lines(c(120, 120), c(-40, -15), lwd = 3, col = "slate grey",lty = 1)
        lines(c(130, 130), c(-15, -10), lwd = 3, col = "slate grey",lty = 1)
        lines(c(40, 60), c(-30, -30), lwd = 3, col = "slate grey",lty = 1)
        lines(c(60, 130), c(-15, -15), lwd = 3, col = "slate grey",lty = 1)
        lines(c(40, 60), c(-10, -10), lwd = 3, col = "slate grey",lty = 1)
        lines(c(110, 130), c(-10, -10), lwd = 3, col = "slate grey",lty = 1)
        lines(c(100, 110), c(-5, -5), lwd = 3, col = "slate grey",lty = 1)
        lines(c(40, 100), c(10, 10), lwd = 3, col = "slate grey",lty = 1)
        text(67.5, 15, "1", col = "red",  cex = 1.5)
        text(60, -2.5, "2", col = "red",  cex = 1.5)
        text(35, -35, "3", col = "red",  cex = 1.5)
        text(80, -27.5, "4", col = "red",  cex = 1.5)
        text(85, -2.5, "5", col = "red",  cex = 1.5)
        text(90, 15, "6", col = "red",  cex = 1.5)
     }
    if (sp %in% c("BETcore","YFTcore")) {
        lines(c(20, 140), c(-35, -35), lwd = 3, col = "slate grey",lty = 1)
        lines(c(20, 125.5), c(-15, -15), lwd = 3, col = "slate grey",lty = 1)
        lines(c(20, 100), c(10, 10), lwd = 3, col = "slate grey",lty = 1)
        lines(c(80, 80), c(10, -15), lwd = 3, col = "slate grey",lty = 1)
        xoffset <- 5
        yoffset <- 2.5
        text(115 + xoffset, -11 - yoffset, "N", col = "red",  cex = 1.5)
        text(115 + xoffset, -40 - yoffset, "S", col = "red",  cex = 1.5)
    }
    maps::map("world", yaxt = "n", xaxt = "n", add = T, resolution = 1, interior=F,fill=T)
    box(lwd = 3)
    axis(1, at = seq(lims[1], lims[2], by = 10), labels = F)
    axis(2, at = seq(lims[3], lims[4], by = 5), labels = F)
    latseq <- seq(lims[3] + 10, lims[4] - 10, by = 10)
    latseq2 <- as.character(latseq)
    lonseq <- seq(lims[1] + 20, lims[2] - 10, by = 20)
    lonseq2 <- as.character(lonseq)
    latseq2[latseq < 0] <- paste(abs(latseq[latseq < 0]), "S",
        sep = "")
    latseq2[latseq > 0] <- paste(latseq[latseq > 0], "N", sep = "")
    lonseq2[lonseq < 180] <- paste(lonseq2[lonseq < 180], "E",
        sep = "")
    lonseq2[lonseq > 180] <- paste(360 - lonseq[lonseq > 180],
        "W", sep = "")
    axis(2, at = latseq, labels = latseq2, cex.axis = 0.75)
    axis(1, at = lonseq, labels = lonseq2, cex.axis = 0.75)
    mtext(side = 3, line = 0.5, plot_title, font = 2, cex = 1.1)
}

storesumm <- function(mod,fname) {
  summry <- summary(mod)
  save(summry,file=paste(fname,".RData"))
  a <- capture.output(summry);
  cat(a,file=paste(fname,".txt"),sep="\n",append=F)
  }

combine_delta_xl <- function(fnamedelta,fnamepos) {
    xl_delta <- read.csv(paste(fnamedelta,".csv",sep=""))
    xl_pos <- read.csv(paste(fnamepos,".csv",sep=""))
    pos <- match(xl_pos[,2],xl_delta[,2])
    coefs.boat <- xl_delta[pos,5] * xl_pos[,5]
    coefs.base <- xl_delta[pos,3] * xl_pos[,3]
    yrpos <- xl_delta[pos,2]
    fishlab <- switch(runsp,yft="Yellowfin",bet="Bigeye");
    plot.slope.ratio(coefs.base,coefs.boat,yrpos,titl=paste("Region",runreg,fishlab,"Delta lognormal combined"))
#    par(mar=c(5,4,1,1))
#    plot(yrpos,coefs.base,type="l",ylab="Relative abundance estimate",xlab="Year",ylim=c(0,2.5))
#    lines(yrpos,coefs.boat,col="red")
    fname2 <- gsub("deltabin","deltacomb",fnamedelta)
    savePlot(paste(fname2,".png",sep=""),type="png")
    write.csv(cbind(yrpos,coefs.base,coefs.boat),file=paste(fname2,".csv",sep=""))
    graphics.off()
    }

combine_delta_xl_indices <- function(fnamedelta,fnamepos) {
    xl_delta <- read.csv(paste(fnamedelta,".csv",sep=""))
    xl_pos <- read.csv(paste(fnamedelta,".csv",sep=""))
    coefs.boat <- xl_delta[,3] * xl_pos[,3]
    fishlab <- switch(runsp,yft="Yellowfin",bet="Bigeye");
    yrpos <- xl_delta[,2]
    windows();par(mar=c(5,4,1,1))
    plot(yrpos,coefs.boat,type="l",ylab="Relative abundance estimate",xlab="Year",ylim=c(0,2.5))
    fname2 <- gsub("deltabin","deltacomb",fnamedelta)
    savePlot(paste(fname2,".png",sep=""),type="png")
    write.csv(cbind(yrpos,coefs.boat),file=paste(fname2,".csv",sep=""))
    graphics.off()
    }

combine_delta <- function(fnamedelta,fnamepos) {
    load(paste("model.",fnamedelta,".base.RData",sep=""))
    load(paste("model.",fnamedelta,".boat.RData",sep=""))
    yrbin <- as.numeric(model.base$xlevels[[1]])
    a <- length(yrbin)
    coefsdeltabin.base <- get.coefs(model.base,a)
    coefsdeltabin.boat <- get.coefs(model.boat,a)
    rm(model.base,model.boat); gc()
    load(paste("model.",fnamepos,".base.RData",sep=""))
    load(paste("model.",fnamepos,".boat.RData",sep=""))
    yrpos <- as.numeric(model.base$xlevels[[1]])
    a <- length(yrpos)
    coefsdeltapos.base <- get.coefs(model.base,a)
    coefsdeltapos.boat <- get.coefs(model.boat,a)
    rm(model.base,model.boat); gc()
    a <- match(names(coefsdeltapos.base),names(coefsdeltabin.base))
    coefs.base <- coefsdeltabin.base[a] * coefsdeltapos.base
    coefs.boat <- coefsdeltabin.boat[a] * coefsdeltapos.boat
    fishlab <- switch(runsp,yft="Yellowfin",bet="Bigeye");
    plot.slope.ratio(coefs.base,coefs.boat,yrpos,titl=paste("Region",runreg,fishlab,"Delta lognormal combined"))
    par(mar=c(5,4,1,1))
    plot(yrpos,coefs.base,type="l",ylab="Relative abundance estimate",xlab="Year",ylim=c(0,2.5))
    lines(yrpos,coefs.boat,col="red")
    fname2 <- paste(fname," deltacomb",sep="")
    savePlot(paste(fname2,".png",sep=""),type="png")
    write.csv(cbind(yr,coefs.base,coefs.boat),file=paste(fname,".csv",sep=""))
    graphics.off()
    }

plot_catchmap <- function(indat=a,vbl,dcd,latlim=c(-40,20),lonlim=c(20,130),brk=seq(0,1,.05),brk2=seq(0,1,.1),ti="",ti2=dcd) {
  plot(1:5,1:5,ylim=latlim,xlim=lonlim,type="n",xlab="Longitude",ylab="Latitude")
  image(lonlim,latlim,matrix(1,nrow=2,ncol=2),add=T,col="light blue")
  indat <- cbind(indat,vbl)
  indat <- indat[indat$lon <= 140,]
  div <- min(diff(sort(unique(indat$lon))))
  test <- indat[indat$decade==dcd,]
  if(dim(test)[1] > 0) {
    a1 <- with(indat[indat$decade==dcd,],tapply(vbl,list(factor(lon,levels=seq(min(lon),max(lon),div)),factor(lat,levels=seq(min(lat),max(lat),div))),sum))
    image(as.numeric(rownames(a1)),as.numeric(colnames(a1)),a1,add=T,col=rev(heat.colors(length(brk)-1)),breaks=brk)
    contour(as.numeric(rownames(a1)),as.numeric(colnames(a1)),a1,add=T,levels=brk2)
    maps::map(database="world",add=T, interior=F,fill=T)
  }
  title(paste(ti2,ti))
  }


plot_cpuemap <- function(indat=a,vb1,vb2,dcd,latlim=c(-40,40),lonlim=c(120,210),brk=seq(0,1,.05),brk2=seq(0,1,.1),ti="") {
  plot(1:5,1:5,ylim=latlim,xlim=lonlim,type="n",xlab="Longitude",ylab="Latitude")
  indat <- cbind(indat,vb1,vb2)
  indat <- indat[indat$lon <= 210,]
  a1 <- with(indat[indat$decade==dcd,],tapply(vb1,list(lon,lat),sum)/tapply(vb2,list(lon,lat),sum))
  image(as.numeric(rownames(a1)),as.numeric(colnames(a1)),a1,add=T,col=rev(heat.colors(length(brk)-1)),breaks=brk)
  contour(as.numeric(rownames(a1)),as.numeric(colnames(a1)),a1,add=T,levels=brk2)
  maps::map(database="world",add=T, interior=F,fill=T)
  title(paste(dcd,ti))
  }

boxplots_spCL <- function(dat=dataset,cl="kmeans",ti="",outL=T,nsp=13) {
  if(nsp==8) windows(30,20);par(mfrow=c(2,4),mar=c(3,3,3,3))
  if(nsp==12) windows(30,20);par(mfrow=c(3,4),mar=c(3,3,3,3))
  if(nsp==13) windows(30,20);par(mfrow=c(3,5),mar=c(3,3,3,3))
  wd=table(dat[,cl])
  for(sp in allsp) {
    boxplot(dat[,sp]/dat$hooks ~ dat[,cl],width=wd,main=sp,outline=outL,pars = list(boxwex = 1))
    }
  title(paste0("R",r," ",ti," ",cl),outer=T,line=-1,cex.main=1.5)
  savePlot(paste0("boxplots_spCL",cl,ti,"R",r,".png"),type="png")
  }

aggregate_by_trip <- function(dat=newdat,flds) {
  dt1 <- data.table(dat)
  setkey(dt1,TRIP_NUM)
  df2 <- dt1[, lapply(.SD, mean, na.rm=T), by=list(TRIP_NUM),.SDcols=flds]
  df2 <- data.frame(df2)
  df3 <- cbind(TRIP_NUM=df2[,1],df2[,flds]/apply(df2[,flds],1,sum))
  return(df3)
  }

# PCA by set
Mode <- function(x) {
# from http://stackoverflow.com/questions/2547402/standard-library-function-in-r-for-finding-the-mode
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

make_newdat <- function(model,datx) {
  vessidx <- Mode(datx$vessid)
  latlongx <- Mode(datx$latlong)
  newdat <- expand.grid(yrqtr=sort(unique(datx$yrqtr)),vessid=vessidx,latlong=latlongx,hooks=median(datx$hooks),hbf=median(datx$hbf))
  predresp <- predict(model,newdata=newdat,type="response",se.fit=T)
  predterms <- predict(model,newdata=newdat,type="terms",se.fit=T)
  return(list(newdat=newdat,predresp=predresp,predterms=predterms))
  }

make_newdat2 <- function(model,datx) {
  vessidx <- Mode(datx$vessid)
  latlongx <- Mode(datx$latlong)
  newdat <- expand.grid(yrqtr=sort(unique(datx$yrqtr)),vessid=vessidx,latlong=latlongx,hooks=median(datx$hooks),hbf=median(datx$hbf),moon=0.5)
  predresp <- predict(model,newdata=newdat,type="response",se.fit=T)
  predterms <- predict(model,newdata=newdat,type="terms",se.fit=T)
  return(list(newdat=newdat,predresp=predresp,predterms=predterms))
  }

make_newdat_TW <- function(model,datx) {
  vessidx <- Mode(datx$vessid)
  latlongx <- Mode(datx$latlong)
  newdat <- expand.grid(yrqtr=sort(unique(datx$yrqtr)),vessid=vessidx,latlong=latlongx,hooks=median(datx$hooks),
    bt1=sort(unique(datx$bt1))[1],bt2=sort(unique(datx$bt2))[1],bt3=sort(unique(datx$bt3))[1],
    bt4=sort(unique(datx$bt4))[1],bt5=sort(unique(datx$bt5))[1],moon=0.5)
  predresp <- predict(model,newdata=newdat,type="response",se.fit=T)
  predterms <- predict(model,newdata=newdat,type="terms",se.fit=T)
  return(list(newdat=newdat,predresp=predresp,predterms=predterms))
  }

make_newdat_JP <- function(model,datx) {
  vessidx <- Mode(datx$vessid)
  latlongx <- Mode(datx$latlong)
  newdat <- expand.grid(yrqtr=sort(unique(datx$yrqtr)),vessid=vessidx,latlong=latlongx,hooks=median(datx$hooks),
            moon=0.5,hbf=median(datx$hbf,na.rm=T))
  predresp <- predict(model,newdata=newdat,type="response",se.fit=T)
  predterms <- predict(model,newdata=newdat,type="terms",se.fit=T)
  return(list(newdat=newdat,predresp=predresp,predterms=predterms))
  }

ls.objects <- function (pos = 1, pattern, order.by, decreasing = FALSE, head = FALSE, n = 5)
{
    napply <- function(names, fn) sapply(names, function(x) fn(get(x,
        pos = pos)))
    names <- ls(pos = pos, pattern = pattern)
    obj.class <- napply(names, function(x) as.character(class(x))[1])
    obj.mode <- napply(names, mode)
    obj.type <- ifelse(is.na(obj.class), obj.mode, obj.class)
    obj.size <- napply(names, object.size)
    obj.dim <- t(napply(names, function(x) as.numeric(dim(x))[1:2]))
    vec <- is.na(obj.dim)[, 1] & (obj.type != "function")
    obj.dim[vec, 1] <- napply(names, length)[vec]
    out <- data.frame(obj.type, obj.size, obj.dim)
    names(out) <- c("Type", "Size", "Rows", "Columns")
    if (!missing(order.by))
        out <- out[order(out[[order.by]], decreasing = decreasing),
            ]
    if (head)
        out <- head(out, n)
    out
}

lsos <- function(...,n=10) {
  ls.objects(...,order.by="Size",decreasing=TRUE,head=TRUE,n=n)
  }

limit_vessels <- function(gdat,minqtrs,maxqtrs=10000) {
  a <- table(gdat$vessid,gdat$yrqtr)
  a <- apply(a>0,1,sum)
  table(a)
  a <- a[a >= minqtrs & a <= maxqtrs]
  gdat <- gdat[gdat$vessid %in% names(a),]
  a <- table(gdat$yrqtr);a
  a <- a[a>=100]
  gdat <- gdat[gdat$yrqtr %in% names(a),]
  a <- table(gdat$vessid);a
  a <- a[a>=100]
  gdat <- gdat[gdat$vessid %in% names(a),]
  return(gdat)
  }

factoralong <- function(x,diffs) {
  stt <- min(x,na.rm=T)
  stp <- max(x,na.rm=T)
  x <- factor(x,levels=seq(stt,stp,diffs))
  return(x)
  }