R/utilities.r

Defines functions get.trace geomean get.stat4 get.stat to_vpa_data read.vpa out.vpa make_summary_table convert_faa_perSPR get.SPR ref.F dyn.msy Generation.Time solv.Feq.Pope solv.Feq catch_equation caa.est.mat_wrong caa.est.mat calc.rel.abund

Documented in caa.est.mat caa.est.mat_wrong calc.rel.abund convert_faa_perSPR Generation.Time get.SPR out.vpa read.vpa ref.F to_vpa_data

#'
#' @import ggplot2
#' @import magrittr
#' @import dplyr
#' @import tidyr
#' @import tibble
#' @import readr
#' @import forcats
#' @import stringr
#' @import assertthat
#' @import patchwork
#' @importFrom magrittr %>%
#' @importFrom magrittr %T>%
#' @importFrom dplyr filter
#'
NULL

# 将来予測や管理基準値計算に使われる関数群 ----

# 個体群動態に関する基本的関数 ----

#'
#' 年齢別パラメータを与えて、1年分前進計算する
#'
#' @param min.age 最小の年齢(0歳がデフォルト)。数でなくて名前
#' @param max.age 最大の年齢。数でなくて名前なので、0からスタートして2歳までの場合には、min.age=0, max.age=2とするか、min.age=1, max.age=3とする
#'
#' @export
#' @encoding UTF-8

calc.rel.abund <- function(sel,Fr,na,M,waa,waa.catch=NULL,maa,min.age=0,max.age=Inf,Pope=TRUE,ssb.coef=0){

  if(is.null(waa.catch)) waa.catch <- waa
  rel.abund <- rep(NA, na)
  rel.abund[1] <- 1

  for (i in seq_len(na-1)) {
    rel.abund[i+1] <- rel.abund[i]*exp(-M[i]-sel[i]*Fr)
  }
  rel.abund[na] <- rel.abund[na-1]*exp(-M[na-1]-sel[na-1]*Fr)*(1-exp(-((max.age-min.age)-(na-2))*(M[na]+sel[na]*Fr)))/(1-exp(-M[na]-sel[na]*Fr))

  if(isTRUE(Pope)){
    ypr1 <- rel.abund*waa.catch[1:na]*(1-exp(-sel[1:na]*Fr))*exp(-M[1:na]/2)
  }
  else{
    # use Baranov catch equation
    ypr1 <- rel.abund*(1-exp(-sel[1:na]*Fr-M[1:na]))*sel[1:na]*Fr/
      (sel[1:na]*Fr+M[1:na])*waa.catch[1:na]
  }
  spr <- rel.abund*waa[1:na]*maa[1:na]*exp(-ssb.coef*(sel[1:na]*Fr+M[1:na]))
  return(list(rel.abund=rel.abund,ypr=ypr1,spr=spr))
}

#'
#' 年齢別生物パラメータとFと漁獲量を与えると与えた漁獲量と一致するFへの乗数を返す
#'
#' @param max_exploitation_rate 潜在的に漁獲できる漁獲量<入力した漁獲量の場合、潜在的に漁獲できる漁獲量の何%まで実際に漁獲するか
#' @param max_F F at ageの最大値となる値の上限をどこにおくか
#'
#' @export
#' @encoding UTF-8

caa.est.mat <- function(naa,saa,waa,M,catch.obs,Pope,set_max1=TRUE,max_exploitation_rate=0.99,max_F=exp(10)){

  tmpfunc <- function(logx,catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,M=M,out=FALSE,Pope=Pope){
    if(Pope==1 | Pope==TRUE) is.pope <- TRUE else is.pope <- FALSE
    x <- exp(logx)
    if(is.pope){
      caa <- naa*(1-exp(-saa*x))*exp(-M/2)
    }
    else{
      caa <- naa*(1-exp(-saa*x-M))*saa*x/(saa*x+M)
    }
    wcaa <- caa*waa
    if(out==FALSE){
        return(log((sum(wcaa,na.rm=T)-catch.obs)^2))
    }
    else{
      return(caa)
    }
  }

  C0 <- sum(tmpfunc(logx=100,catch.obs=catch.obs,naa=naa,saa=rep(1,length(saa)),waa=waa,M=M,Pope=Pope,out=TRUE) * waa)
  if(C0 < catch.obs){
    warning("The expected catch (", catch.obs, ") is over potential maximum catch (",round(C0,5),"). The expected catch is replaced by",round(C0,3),"x", max_exploitation_rate)
    catch.obs <- C0 * max_exploitation_rate

  }

  # caa.est.mat_wrongで初期値を設定するとはやそう
  tmp <- optimize(tmpfunc,c(-10,log(max_F)),catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,M=M,Pope=Pope,out=FALSE)#,tol=.Machine$double.eps)
  tmp2 <- tmpfunc(logx=tmp$minimum,catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,M=M,Pope=Pope,out=TRUE)
  realized.catch <- sum(tmp2*waa)
  if(!is.nan(realized.catch/catch.obs) & abs(realized.catch/catch.obs-1)>0.1) warning("expected catch:",catch.obs,", realized catch:",realized.catch)
  return(list(x=exp(tmp$minimum),caa=tmp2,realized.catch=realized.catch, expected.catch=catch.obs))
}

#'
#' 年齢別生物パラメータとFと漁獲量を与えると与えた漁獲量と一致するFへの乗数を返す (optimを使わないでやろうとしたけどだめだったやつ)
#'
#' これだと漁獲量をぴったりに返すFのvectorは得られるが、もとの選択率に一致しない
#'
#' @param naa numbers at age
#' @param saa selectivity at age
#' @param waa weight at age
#' @param M natural mortality at age
#' @param max_exploitation_rate 潜在的に漁獲できる漁獲量<入力した漁獲量の場合、潜在的に漁獲できる漁獲量の何%まで実際に漁獲するか
#' @param max_F F at ageの最大値となる値の上限をどこにおくか
#'
#' @export
#' @encoding UTF-8

caa.est.mat_wrong <- function(naa,saa,waa,M,catch.obs,Pope,max_exploitation_rate=0.99,max_F=exp(10)){

  C0 <- catch_equation(naa,exp(100),waa,M,Pope) %>% sum()
  if(C0 < catch.obs){
    warning("The expected catch (", catch.obs, ") is over potential maximum catch (",round(C0,5),"). The expected catch is replaced by",round(C0,3),"x", max_exploitation_rate)
    catch.obs <- C0 * max_exploitation_rate
  }

  caa_original <- catch_equation(naa,saa,1,M,Pope)
  catch_ratio <- catch.obs / sum(caa_original * waa)
  caa_expected <- caa_original * catch_ratio
  if(Pope==0){
    Fvector <- solv.Feq(caa_expected, naa, M)
  }
  else{
    Fvector <- solv.Feq.Pope(caa_expected, naa, M)
  }

  if( max_F < max(Fvector)) Fvector <- max_F * Fvector/max(Fvector)

  realized.catch <- catch_equation(naa,Fvector,waa,M,Pope)
  if(abs(sum(realized.catch)/catch.obs-1)>0.1) warning("expected catch:",catch.obs,", realized catch:",realized.catch)
  multi <- mean(Fvector/saa)
  return(list(x=multi,caa=realized.catch,realized.catch=sum(realized.catch), expected.catch=catch.obs, Fvector=Fvector))
}

# # 上の関数とどっちが使われているか不明,,,多分使われていないのでコメントアウトする
# caa.est <- function(naa,saa,waa,M,catch.obs,Pope){
#   saa <- saa/max(saa)
#   tmpfunc <- function(x,catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,M=M,out=FALSE,Pope=Pope){
#     if(isTRUE(Pope)){
#       caa <- naa*(1-exp(-saa*x))*exp(-M/2)
#     }
#     else{
#       caa <- naa*(1-exp(-saa*x-M))*saa*x/(saa*x+M)
#     }
#     wcaa <- caa*waa
#     if(out==FALSE){
#       return((sum(wcaa,na.rm=T)-catch.obs)^2)
#     }
#     else{
#       return(caa)
#     }
#   }
#   tmp <- optimize(tmpfunc,c(0,5),catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,M=M,Pope=Pope,out=FALSE)
#   tmp2 <- tmpfunc(x=tmp$minimum,catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,M=M,Pope=Pope,out=TRUE)
#   return(list(x=tmp$minimum,caa=tmp2))
# }

catch_equation <- function(naa,faa,waa,M,Pope=1){
  if(Pope==1 | Pope==TRUE) is.pope <- TRUE else is.pope <- FALSE
  if(is.pope){
    wcaa_mat <- naa*(1-exp(-faa))*exp(-M/2) * waa
  }
  else{
    wcaa_mat <- naa*(1-exp(-faa-M))*faa/(faa+M) * waa
  }
  return(wcaa_mat)
}

#' @export
#' @encoding UTF-8

solv.Feq <- function(cvec,nvec,mvec){
  Fres <- rep(0,length(cvec))
  # cat(nvec," ")
  for(i in 1:length(cvec)){
    F0 <- cvec[i]/nvec[i]
    F1 <- cvec[i]*(F0+mvec[i])/nvec[i]/(1-exp(-F0-mvec[i]))
    if(round(cvec[i],6)<round(nvec[i],6)){
      while(abs(F0-F1)>0.0001 ){
        F0 <- F1
        F1 <- cvec[i]*(F0+mvec[i])/nvec[i]/(1-exp(-F0-mvec[i]))
        if(F0-F1==-Inf) cat("\n",cvec[i]," ",nvec[i]," \n")
      }
      Fres[i] <- F1
    }
    else{
      Fres[i] <- 10
      cat("Warning: catch exceeded tot_num at: ",i," ",
          round(cvec[i],6)," ",round(nvec[i],6),"\n")
    }
  }
  Fres
}

#' @export
solv.Feq.Pope <- function(cvec,nvec,mvec){
  -log(1-cvec/nvec/exp(-mvec/2))
}

#' VPAの結果に格納されている生物パラメータから世代時間を計算
#'
#' @export
#' @encoding UTF-8

Generation.Time <- function(vpares,
                            maa.year=2014:2015,
                            maa=NULL,
                            M.year=2014:2015,
                            M=NULL,
                            age=NULL,
                            Plus = 19
){

  if(is.null(maa)){
    maa <- vpares$input$dat$maa
    maa <- rowMeans(maa[,colnames(maa) %in% maa.year,drop=F],na.rm=T)
    maa <- maa[!is.na(maa)]
  }
  if(is.null(M)){
    M <- vpares$input$dat$M
    M <- rowMeans(M[,colnames(M) %in% M.year,drop=F],na.rm=T)
    M <- M[!is.na(M)]
  }

  if(is.null(age)) age <- as.numeric(rownames(vpares$naa))
  if(Plus>0){
    maa <- c(maa, rep(1,Plus))
    M <- c(M, rep(M[length(M)],Plus))
    age <- c(age, max(age)+1:Plus)
  }
  A <- length(M)
  L <- c(1,exp(-cumsum(M[-A])))
  G <- sum(age*L*maa)/sum(L*maa)

  return(G)
}

### dynamics MSYを計算してみる
#' @export
#'

dyn.msy <- function(naa.past,naa.init=NULL,fmsy,a,b,resid,resid.year,waa,maa,M,assume_SR=TRUE, SR="HS"){

  if(assume_SR==TRUE){
    if (SR=="HS") SRF <- function(x,a,b) ifelse(x>b,b*a,x*a)
    if (SR=="BH") SRF <- function(x,a,b) a*x/(1+b*x)
    if (SR=="RI") SRF <- function(x,a,b) a*x*exp(-b*x)
    if (SR=="Mesnil") SRF <- function(x,a,b) 0.5*a*(x+sqrt(b^2+gamma^2/4)-sqrt((x-b)^2+gamma^2/4))
  }

  nyear <- length(resid)
  if(is.null(naa.init)) nage <- nrow(naa.past) else nage <- length(naa.init)
  naa <- matrix(0,nage,nyear)
  ssb <- numeric()
  if(is.null(naa.init)) naa[,1] <- naa.past[,colnames(naa.past)==min(resid.year)]
  else naa[,1] <- naa.init
  colnames(naa) <- resid.year
  if(is.null(naa.init)){
    waa <- waa[,colnames(naa.past)%in%resid.year]
    maa <- maa[,colnames(naa.past)%in%resid.year]
    M <- M[,colnames(naa.past)%in%resid.year]
  }
  for(i in 2:nyear){
    ssb[i-1] <- sum(naa[,i-1]*waa[,i-1]*maa[,i-1],na.rm=T)
    if(assume_SR==TRUE){
      naa[1,i] <- SRF(ssb[i-1],a,b)*exp(resid[i])
    }
    else{
      naa[1,i] <- naa.past[1,i]
    }
    for(j in 2:(nage-1)) naa[j,i] <- naa[j-1,i-1] * exp(-fmsy[j-1]-M[j-1,i-1])
    naa[nage,i] <- naa[nage-1,i-1] * exp(-fmsy[j-1]-M[j-1,i-1]) + naa[nage,i-1] * exp(-fmsy[nage]-M[nage,i-1])
  }
  i <- nyear ; ssb[i] <- sum(naa[,i]*waa[,i]*maa[,i], na.rm=TRUE)
  list(naa=naa,ssb=ssb)
}



# 再生産関係を仮定しない管理基準値計算 ----

#'
#' 再生産関係を仮定しない管理基準値計算(SPR,YPR,F0.1,Fmax)のための関数
#'
#' @param res VPAの出力結果(NULLも可)。ここがNULLの場合(VPAの出力結果を与えない場合)でも、Fcurrent, waa, maa, M, waa.catch, max.age, min.age, Popeを別途指定することによって、管理基準値計算ができるようになる
#' @param Fcurrent 仮定する選択率.NULLの場合,res$Fc.at.ageが使われる
#' @param waa 仮定する年齢別体重。直接の値を入れるか,waa.yearで年を指定するやり方のどちらでも動く。直接指定するほうが優先。
#' @param maa 仮定する年齢別成熟率。直接の値を入れるか,waa.yearで年を指定するやり方のどちらでも動く。直接指定するほうが優先。
#' @param M 仮定する年齢別死亡率。直接の値を入れるか,waa.yearで年を指定するやり方のどちらでも動く。直接指定するほうが優先。
#' @param waa.catch 仮定する年齢別体重(漁獲量計算用)。直接の値を入れるか,waa.yearで年を指定するやり方のどちらでも動く。直接指定するほうが優先。
#' @param M.year 年を指定して生物パラメータを仮定する場合.年の範囲の平均値が用いられる.NULLの場合,VPA最終年の値が使われる
#' @param waa.year 年を指定して生物パラメータを仮定する場合.年の範囲の平均値が用いられる.NULLの場合,VPA最終年の値が使われる
#' @param maa.year 年を指定して生物パラメータを仮定する場合.年の範囲の平均値が用いられる.NULLの場合,VPA最終年の値が使われる
#' @param rps.year Fmedの計算に使うRPSの年の範囲.NULLの場合,全範囲が用いられる
#' @param rps.vector Fmedの計算に使うRPSのベクトル。rps.yearよりもこちらが優先される。
#' @param max.age 加入年齢を0歳としたときに、SPR計算で考慮される最大の年齢の名前(何行目とかじゃないことに注意, デフォルトはInf)。min.ageも同様。VPA結果を与える場合にはVPA結果から自動的にもってくる。
#' @param min.age VPA結果を与える場合にはVPA結果から自動的にもってくるが、VPA結果を与えない場合、加入年齢を入力する
#' @param  pSPR = seq(10,90,by=10), # F\%SPRを計算するときの%SPR
#' @param d 0.001
#' @param  Fem.init 経験的管理基準値(Fmed, Fmean, Fhigh, Flow)の初期値 (default=0.5)
#' @param  Fmax.init Fmaxの初期値 (default=1.5)
#' @param  F0.1.init F0.1の初期値 (default=0.7)
#' @param  iterlim
#' @param  plot 結果のプロットを表示するかどうか
#' @param  Pope Popeの式を使うか
#' @param  F.range YPR, SPR曲線を書くときのFの範囲(Fの最大値のスケール)、かつ、F\%SPRを計算するときの初期値を決めるために利用される。F\%SPRの推定がうまくいかない場合はこの範囲を調整してください。
#' @encoding UTF-8
#'
#' @note F_SPRのF管理基準値の初期値は 与えられたFのもとでのSPR/目的のSPR を初期値とするように調整されるので不要。プラスグループを考慮するかどうかはVPAの結果のinput$plus.groupから自動判別される。
#'
#' @examples
#' \dontrun{
#' data(res_vpa_example)
#' # VPAデータを使う場合
#' res_refF1 <- ref.F(res=res_vpa_example,Fcurrent=frasyr::apply_year_colum(res_vpa_example$faa,2015:2017),
#'                 waa.year=2015:2017,maa.year=2015:2017,M.year=2015:2017)
#'
#' # 生物パラメータをデータとして与える場合
#' res_refF2 <- ref.F(res=NULL,Fcurrent=rep(0.1,5),
#'                    waa=rep(100,5),maa=c(0,0,1,1,1),M=rep(0.3,5),waa.catch=rep(100,5),
#'                    rps.vector=NULL, # Fmedを計算したりする場合のRPSのベクトル.NULLでもOK
#'                    Pope=TRUE,min.age=0,pSPR=c(30,40))
#' }
#'
#'
#' @export
#' @import tibble
#' @encoding UTF-8
#'

# ref.F
ref.F <- function(
  res=NULL, # VPAの結果のオブジェクト
  Fcurrent=NULL, # Fcurrentの仮定.NULLの場合,res$Fc.at.ageが使われる
  waa=NULL, # 仮定する生物パラメータ.直接の値を入れるか,年を指定するやり方のどちらでも動く。直接指定するほうが優先。
  maa=NULL,
  M=NULL,
  waa.catch=NULL,
  M.year=NULL,
  waa.year=NULL, # 年を指定して生物パラメータを仮定する場合.年の範囲の平均値が用いられる.NULLの場合,VPA最終年の値が使われる
  waa.catch.year=NULL,
  maa.year=NULL,
  rps.year = NULL, # Fmedの計算に使うRPSの年の範囲.NULLの場合,全範囲が用いられる
  rps.vector = NULL,
  max.age = Inf,
  min.age = 0,
  d = 0.001,
  Fem.init = 0.5,
  Fmax.init = 1.5, # Fmaxの初期値
  F0.1.init = 0.7, # F0.1の初期値
  pSPR = seq(10,90,by=10), # F%SPRを計算するときの%SPR
  iterlim=1000,
  plot=TRUE,
  Pope=NULL, # 2014.7.4追加
  F.range = seq(from=0,to=2,length=101)  # YPR, SPR曲線を書くときのFの範囲
){

  argname <- ls()
  arglist <- lapply(argname,function(x) eval(parse(text=x)))
  names(arglist) <- argname

  if(!is.null(res)){
    if(is.null(Pope)) Pope <- res$input$Pope

    naa <- res$naa
    ssb <- res$ssb
    ny <- ncol(naa)
    years <- dimnames(naa)[[2]]
    ages <- dimnames(naa)[[1]]

    if(is.null(Fcurrent)){
      Fcurrent <- res$Fc.at.age
    }
    sel <- Fcurrent/max(Fcurrent,na.rm=TRUE)

    if(is.null(waa.year)) waa.year <- rev(years)[1]
    if(is.null(maa.year)) maa.year <- rev(years)[1]
    if(is.null(M.year)) M.year <- rev(years)[1]
    if(is.null(waa.catch.year)) waa.catch.year <- rev(years)[1]

    if(is.null(waa))  waa <- apply_year_colum(res$input$dat$waa,waa.year)
    if(is.null(M))    M   <- apply_year_colum(res$input$dat$M,M.year)
    if(is.null(maa))  maa <- apply_year_colum(res$input$dat$maa,maa.year)
    if(is.null(waa.catch)){
        if(!is.null(res$input$dat$waa.catch)) waa.catch <- apply_year_colum(res$input$dat$waa.catch,waa.catch.year)
        else waa.catch <- waa
    }

## Remove duplicate code  
#    if(is.null(waa.catch)){
#      if(is.null(res$input$dat$waa.catch)){
#        waa.catch <- waa
#      }
#      else{
#       waa.catch <- apply_year_colum(res$input$dat$waa.catch,waa.year) # here, waa.year should be waa.catch.year
#      }
#    }

    na <- sum(!is.na(Fcurrent))
    ssb.coef <- ifelse(is.null(res$ssb.coef),0,res$ssb.coef)

    min.age <- min(as.numeric(rownames(res$naa)))
    if(min.age==0) slide.tmp <- TRUE else slide.tmp <- -1:-min.age
   
    if(!is.null(rps.year)){
      rps.data <- data.frame(year=as.numeric(names(colSums(ssb,na.rm=T))),
                             ssb=as.numeric(colSums(ssb,na.rm=T)),
                             recruit=as.numeric(c(naa[1,slide.tmp],rep(NA,min.age))))
      if (sum(is.na(rps.data$year))>0) rps.data <- rps.data[-which(is.na(rps.data$year)),]
      rps.data$rps <- rps <- rps.data$recruit/rps.data$ssb
      #  rps <- as.numeric(naa[1,]/colSums(ssb,na.rm=TRUE))

      #  if (is.null(rps.year)) rps.year <- years

      tmp <- rps.data$year %in% rps.year
      rps.q <- quantile(rps[tmp], na.rm=TRUE, probs=c(0.1,0.5,0.9))
      rps.q <- c(rps.q,mean(rps[tmp], na.rm=TRUE))
      names(rps.q)[4] <- "mean"
      spr.q <- 1/rps.q
    }
    else{
      rps <- NULL
      rps.q <- NULL
      spr.q <- NULL
      rps.data <- NULL
    }
  }
  if(is.null(res)){ # VPA結果を与えない場合

    if(is.vector(Fcurrent) && is.null(names(Fcurrent))){
        names(Fcurrent) <- 1:length(Fcurrent)
    }
    
    sel <- Fcurrent/max(Fcurrent,na.rm=TRUE)
    na <- length(Fcurrent)
    assertthat::assert_that(length(Fcurrent) == na)
    assertthat::assert_that(length(waa) == na)
    assertthat::assert_that(length(maa) == na)
    assertthat::assert_that(length(M)   == na)
    assertthat::assert_that(length(waa.catch) == na)
    assertthat::assert_that(!is.null(Pope))
    ssb.coef <- 0

    if(!is.null(rps.vector)){
      rps <- rps.data <- rps.vector
      rps.q <- quantile(rps, na.rm=TRUE, probs=c(0.1,0.5,0.9))
      rps.q <- c(rps.q,mean(as.numeric(rps), na.rm=TRUE))
      names(rps.q)[4] <- "mean"
      spr.q <- 1/rps.q
    }
    else{
      rps <- NULL
      rps.q <- NULL
      spr.q <- NULL
      rps.data <- NULL
    }
  }

  if(!is.null(res) && res$input$plus.group==FALSE){
    min.age <- min(as.numeric(rownames(res$naa)))
    max.age <- max(as.numeric(rownames(res$naa)))
  }

  calc.rel.abund2_ <- function(Fx,multi){
    calc.rel.abund(Fx,multi,na,M,waa,waa.catch,maa,
                   min.age=min.age,
                   max.age=max.age,
                   Pope=Pope,ssb.coef=ssb.coef)
  }

  original.spr <- calc.rel.abund2_(Fcurrent,1)
  original.spr0 <- calc.rel.abund2_(Fcurrent,0)
  original.perspr <- sum(original.spr$spr,na.rm=T)/sum(original.spr0$spr,na.rm=T)

  # Fcurrent
  Fcurrent_max_mean <- c(max(Fcurrent,na.rm=T), mean(Fcurrent,na.rm=T))

  # grid search
  Fcurrent_max <- Fcurrent_max_mean[1]
  F.range <- sort(c(F.range,  Fcurrent_max))
  spr0 <- sum(original.spr0$spr,na.rm=T)
  tmp <- lapply(F.range, function(x) calc.rel.abund2_(sel,x))
  ypr <- sapply(tmp,function(x) sum(x$ypr,na.rm=T))
  ypr_by_age <- purrr::map_dfr(tmp,function(x) x$ypr)
  colnames(ypr_by_age) <- str_c("TC-mean-A",colnames(ypr_by_age))
  pspr <- sapply(tmp,function(x) sum(x$spr,na.rm=T))/spr0*100
  ypr.spr <- data.frame(F.range=F.range,ypr=ypr,pspr=pspr)
  ypr.spr$Frange2Fcurrent  <- ypr.spr$F.range/Fcurrent_max
  ypr.spr <- cbind(ypr.spr, ypr_by_age)

  # F.spr

  spr.f.est <- function(log.p, out=FALSE, sub="med", spr0=NULL){
    Fr <- exp(log.p)

    tmp <- calc.rel.abund2_(sel,Fr)
    rel.abund <- tmp$rel.abund
    spr <- sum(tmp$spr,na.rm=T)
    if (isTRUE(out)) obj <- spr
    else{
      if(sub=="mean") obj <- (spr-spr.q[4])^2
      if(sub=="low") obj <- (spr-spr.q[3])^2
      if(sub=="med") obj <- (spr-spr.q[2])^2
      if(sub=="high") obj <- (spr-spr.q[1])^2
      if(is.numeric(sub)) obj <- (spr/spr0-sub/100)^2

    }

    return(obj)
  }

  spr0 <- spr.f.est(-Inf, out=TRUE)

  if(!is.null(rps)){
    Fmed.res <- nlm(spr.f.est, Fem.init, out=FALSE, sub="med", iterlim = iterlim)
    Fmean.res <- nlm(spr.f.est, Fem.init, out=FALSE, sub="mean", iterlim = iterlim)
    Flow.res <- nlm(spr.f.est, Fem.init, out=FALSE, sub="low", iterlim = iterlim)
    Fhigh.res <- nlm(spr.f.est, Fem.init, out=FALSE, sub="high", iterlim = iterlim)

    Fmean <- exp(Fmean.res$estimate)
    Fmed <- exp(Fmed.res$estimate)
    Flow <- exp(Flow.res$estimate)
    Fhigh <- exp(Fhigh.res$estimate)
  }
  else{
    Fmean <- Fmed <- Flow <- Fhigh <- NA
  }

  if (!is.null(pSPR)){
    FpSPR <- NULL

    for (i in pSPR){
      tmp <- which.min(abs(ypr.spr$pspr-i))+c(-1,1)
      tmp <- ifelse(tmp<=0,1,tmp)
      tmp <- ifelse(tmp>length(ypr.spr$F.range),length(ypr.spr$F.range),tmp)
      Fspr.init <- log(ypr.spr$F.range[tmp]) #original.perspr/i*100
      Fspr.init[Fspr.init==-Inf] <- -10
      FpSPR.res <- optimize(spr.f.est, Fspr.init,out=FALSE,sub=i,spr0=spr0)
      #nlm(spr.f.est, Fspr.init, out=FALSE, sub=i, spr0=spr0, iterlim=iterlim)
      #      cat("Estimate F%spr: initial value=", Fspr.init," : estimated value=",exp(FpSPR.res$estimate),"\n")
      FpSPR <- c(FpSPR, exp(FpSPR.res$minimum))
    }
    names(FpSPR) <- paste(pSPR,"%SPR",sep="")
  }

  # Fmax

  ypr.f.est <- function(log.p, out=FALSE){
    Fr <- exp(log.p)

    tmp <- calc.rel.abund2_(sel,Fr)
    rel.abund <- tmp$rel.abund
    ypr <- sum(tmp$ypr,na.rm=T)

    if (isTRUE(out)) obj <- ypr else obj <- -ypr

    return(obj)
  }

  Fmax.res <- nlm(ypr.f.est, log(Fmax.init), out=FALSE)

  Fmax <- exp(Fmax.res$estimate)

  # F0.1

  Fp <- function(log.p, out=FALSE){
    Fr <- exp(log.p)

    tmp <- calc.rel.abund2_(sel,Fr)
    rel.abund <- tmp$rel.abund
    ypr <- sum(tmp$ypr,na.rm=T)
    if (isTRUE(out)) obj <- ypr else obj <- -ypr

    return(obj)
  }

  F0.1.est <- function(log.p){
    p <- exp(log.p)
    ref.trend <- (Fp(log(d))-Fp(log(0)))/d
    trend <- (Fp(log(p+d)) - Fp(log(p)))/d

    obj <- (ref.trend/10 - trend)^2

    obj
  }

  F0.1.res <- nlm(F0.1.est,log(F0.1.init))

  F0.1 <- exp(F0.1.res$estimate)


  # output
  f.mean <- function(x) mean(x*sel, na.rm=T)

  Fmean <- c(Fmean, f.mean(Fmean))
  Fmed <- c(Fmed, f.mean(Fmed))
  Flow <- c(Flow, f.mean(Flow))
  Fhigh <- c(Fhigh, f.mean(Fhigh))
  Fmax <- c(Fmax, f.mean(Fmax))
  F0.1 <- c(F0.1, f.mean(F0.1))

  names(Fcurrent_max_mean) <- names(Fmed) <- names(Fmean) <- names(Flow) <- names(Fhigh) <- names(Fmax) <- names(F0.1) <- c("max","mean")
  Res <- list(min.age=min.age, max.age=max.age, rps.q=rps.q, spr.q=spr.q, Fcurrent=Fcurrent_max_mean,
              sel=sel, Fmed=Fmed, Flow=Flow, Fhigh=Fhigh, Fmax=Fmax, F0.1=F0.1, Fmean=Fmean,rps.data=rps.data)

  if (!is.null(pSPR)){
    FpSPR <- rbind(FpSPR, sapply(FpSPR, f.mean))
    rownames(FpSPR) <- c("max","mean")
    Res$FpSPR <- FpSPR
  }

  #---- make summary
  Res$summary <- as.data.frame(Res[substr(names(Res),1,1)=="F"])
  Res$summary <- rbind(Res$summary,Res$summary[1,]/Res$summary[1,1])
  dimnames(Res$summary)[[1]][3] <- "Fref/Fcur"

  Res$currentSPR <- list(SPR=sum(original.spr$spr,na.rm=T),
                         perSPR=original.perspr,
                         YPR=sum(original.spr$spr,na.rm=T),
                         Fcurrent=Fcurrent)

  Res$ypr.spr  <- ypr.spr #data.frame(F.range=F.range,ypr=ypr,spr=spr)
  Res$waa <- waa
  Res$waa.catch <- waa.catch
  Res$maa <- maa
  Res$arglist <- arglist
  Res$spr0 <- spr0

  class(Res) <- "ref"

  if(isTRUE(plot)){
    plot_Fref(Res)
  }
  return(Res)
}

#' 毎年のFの\%SPRやターゲットした\%SPRに相当するFの相対的な大きさを計算する
#'
#' VPA計算結果を使って毎年のF at ageがどのくらいのSPR, YPRに相当するかを計算する。また、各年のFが、目標としたSPR(target.SPR)を達成するためのF(Ftarget)の何倍(F/Ftarget)に相当するかも計算する。F/Ftargetは数値的に探索するが、そのときのF/Ftargetの上限をFmaxにて指定する。十分大きい値(デフォルトは10)を与えておけば大丈夫だが、Ftargetが非常に小さい数字になりそうな場合にはこの値をもっと大きくとるなどする。また、SPRの計算は、デフォルトでは等比級数の和の公式を使って、無限大の年齢までSPRを足しているが、max.ageを指定することで、有限の年齢までの和も計算できる。
#'
#' @param dres vpa関数の返り値
#' @param target.SPR 目標とするSPR。この値を入れると、結果の$ysdata$"F/Ftarget"で、その年のFが目標としたSPR(%)を達成するためのF(Ftarget)の何倍になっているかを返す。デフォルトは30が入っている。このとき、SPRを計算するための生物パラメータ(年齢別体重・成熟率・死亡率)はそれぞれの年で仮定されているものを用いる。
#' @param Fmax F/Ftargetを推定するときに探索するFの乗数の最大値
#' @param max.age SPRやYPRの計算をするときに最大何歳まで考慮するか(年齢のラベルではなく、ベクトルの長さ。デフォルトは無限大)。VPA計算でプラスグループを考慮していない(dres$input$plus.group==FALSE)場合には自動的に設定される。
#'
#'
#' @encoding UTF-8
#'
# #' @examples
# #' \dontrun{
# #'  data(res_vpa_example)
# #'  Fratio <- get.SPR(res_vpa_example,target.SPR=12)$ysdata$"F/Ftarget"
# #'  plot(colnames(res_vpa_example$naa),Fratio,type="b")
# #' }
#'
#' @export
#' @encoding UTF-8

get.SPR <- function(dres,target.SPR=30,Fmax=10,max.age=Inf){
  dres$ysdata <- matrix(0,ncol(dres$faa),5)
  dimnames(dres$ysdata) <- list(colnames(dres$faa),c("perSPR","YPR","SPR","SPR0","F/Ftarget"))
  for(i in 1:ncol(dres$faa)){
    dres$Fc.at.age <- dres$faa[,i] # Fc.at.ageに対象年のFAAを入れる
    if(!all(dres$Fc.at.age==0, na.rm=T)){
      byear <- colnames(dres$faa)[i] # 何年の生物パラメータを使うか

      a <- ref.F(dres,waa.year=byear,maa.year=byear,M.year=byear,rps.year=2000:2011,
                 pSPR=target.SPR,
                 F.range=c(seq(from=0,to=ceiling(max(dres$Fc.at.age,na.rm=T)*Fmax),
                               length=301),max(dres$Fc.at.age,na.rm=T)),plot=FALSE)
      # YPRと%SPR
      dres$ysdata[i,1:2] <- (as.numeric(rev(a$ypr.spr[which(a$ypr.spr$Frange2Fcurrent==1)[1],2:3])))
      # SPR
      dres$ysdata[i,3] <- a$spr0*dres$ysdata[i,1]/100
      # SPR0
      dres$ysdata[i,4] <- a$spr0
      # relative F
      dres$ysdata[i,5] <- 1/a$summary[3,grep("SPR",colnames(a$summary))][1]
    }
    else{
      break;
    }
  }
  dres$ysdata <- as.data.frame(dres$ysdata)
  dres$target.SPR <- target.SPR
  return(dres)
}

#' vpaデータ, 選択率の参照年1, 漁獲圧の参照年2を与えて、参照年2の漁獲圧のもとで参照年1の選択率となるF at ageを作成する関数。漁獲圧の変換は%SPR換算
#'
#' sel_year 選択率の参照年
#' faa_year 漁獲圧の参照年
#'
#' @encoding UTF-8
#' @export

convert_faa_perSPR <- function(res_vpa, sel_year, faa_year, Fcurrent_MSY=NULL, Fsel_MSY=NULL){
  if(is.null(Fcurrent_MSY)){
    Fcurrent_MSY <- apply_year_colum(res_vpa$faa,target_year=faa_year)
  }
  Fcurrent.per.spr <- ref.F(res_vpa,Fcurrent=Fcurrent_MSY,waa.year=faa_year,M.year=faa_year,
                            maa.year=faa_year,plot=FALSE,pSPR=NULL)$currentSPR$perSPR
  cat("%SPR in Fcurrent ",Fcurrent.per.spr,"\n")
  if(is.null(Fsel_MSY)){
    Fsel_MSY <- apply_year_colum(res_vpa$faa,target_year=sel_year)
  }

  Fmultiplier <- ref.F(res_vpa,Fcurrent=Fsel_MSY,waa.year=faa_year,M.year=faa_year,
                       maa.year=faa_year,pSPR=Fcurrent.per.spr*100,plot=FALSE)
  cat("----------------------------- \n ")
  cat("%SPR in Fsel_MSY",Fmultiplier$currentSPR$perSPR," = ")
  Fmultiplier <- rev(Fmultiplier$summary)[3,1] %>% unlist()
  cat("(F multiplier=",Fmultiplier,")\n")
  Fcurrent_MSY_new <- Fsel_MSY * Fmultiplier
  Fref_tmp <- ref.F(res_vpa,Fcurrent=Fcurrent_MSY_new,waa.year=faa_year,
                    M.year=faa_year,maa.year=faa_year,
                    pSPR=Fcurrent.per.spr*100,plot=FALSE)
  cat("%SPR in new-Fcurrent",Fref_tmp$currentSPR$perSPR,"\n")
  cat("----------------------------- \n ")
  matplot(cbind(Fcurrent_MSY_new,Fcurrent_MSY,Fsel_MSY),type="b",pch=1:3,ylab="F",xlab="Age (recruit age=1)")
  legend("bottomright",pch=1:3,col=1:3,c("Fcurrent","Fsel_MSY","new-Fcurrent"))
  return(Fcurrent_MSY_new)
}




#' @export
#' @encoding UTF-8
make_summary_table <- function(mat_data,side=1,probs=c(0.1,0.5,0.8)){
  res_mat <- cbind(apply(mat_data,side,mean),
                   t(apply(mat_data,side,quantile,probs=probs)))
  colnames(res_mat)[1] <- "average"

  as.data.frame(res_mat)
}


# 結果の入出力&サマリーの作成 ----

#'
#' VPA結果をcsvファイルに出力する
#'
#' @param res  VPAの結果
#' @param srres fit.SRの結果
#' @param msyres est.MSYの結果
#' @param fres_current future_vpaの結果(Fcurrent)
#' @param fres_HCR future_vpaの結果(F with HCR)
#' @param kobeII kobeII.matrixの結果
#' @param other_tables その他の表
#' @param filename csvファイルとpdfファイルの両方のファイル名を指定する場合(拡張子なしで指定)
#' @param csvname csvファイルのファイル名
#' @param pdfname pdfファイルのファイル名
#' @param
#'
#' @encoding UTF-8
#' @export

out.vpa <- function(res=NULL,    # VPA result
                    srres=NULL,  # fit.SR result
                    msyres=NULL, # est.MSY result
                    fres_current=NULL,   # future projection result
                    fres_HCR=NULL,   # future projection result
                    kobeII=NULL, # kobeII result
                    kobe.ratio=NULL, # kobe.ratio
                    other_tables=NULL, # other table
                    filename="vpa", # filename without extension
                    csvname=NULL,
                    pdfname=NULL,
                    ci.future=c(0.1,0.5,0.9)
){
  old.par <- par()
  exit.func <- function(){
    dev.off()
    options(warn=0)
  }
  on.exit(exit.func())

  if(!is.null(filename)){
    csvname <- paste(filename,".csv",sep="")
    pdfname <- paste(filename,".pdf",sep="")
  }

  pdf(pdfname)
  par(mfrow=c(3,2),mar=c(3,3,2,1))
  options(warn=-1)

  write.table2 <- function(x,title.tmp="",is.plot=TRUE,...){
    if(is.plot){
      if(!is.null(dim(x))){
        matplot(colnames(x),t(x),type="b",ylim=c(0,max(x,na.rm=T)),pch=substr(rownames(x),1,1))
      }
      else{
        barplot(x)
      }
      title(title.tmp)
    }
    if(!is.null(dim(x))){
      tmp <- matrix("",nrow(x)+1,ncol(x)+1)
      tmp[-1,-1] <- as.character(unlist(x))
      tmp[-1,1] <- rownames(x)
      tmp[1,-1] <- colnames(x)
    }
    else{
      tmp <- x
    }
    write.table(tmp,append=T,sep=",",quote=FALSE,file=csvname,col.names=F,row.names=F,...)
  }

  pd <- packageDescription("frasyr")
  if(is.null(pd$GithubSHA1)) pd$GithubSHA1 <- "local" # fraysrをload_allした場合コミット番号が記録されないため
  write(paste("# frasyr(@",pd$GithubSHA1,") outputs at ",date()," & ",getwd(), sep=""),file=csvname)

  if(!is.null(res)){
    write("# VPA results",file=csvname, append=T)

    write("\n# catch at age", file=csvname,append=T)
    write.table2(res$input$dat$caa,title.tmp="Catch at age")

    write("\n# maturity at age", file=csvname,append=T)
    write.table2(res$input$dat$maa,title.tmp="Maturity at age")

    write("\n# weight at age for biomass calculation", file=csvname,append=T)
    write.table2(res$input$dat$waa,title.tmp="Weight at age (for biomass)")

    if(!is.null(res$input$dat$waa.catch)){
      write("\n# weight at age for catch calculation", file=csvname,append=T)
      write.table2(res$input$dat$waa.catch,title.tmp="Weight at age (for catch)")
    }

    write("\n# M at age",file=csvname,append=T)
    write.table2(res$input$dat$M,title.tmp="M at age")

    write("\n# fishing mortality at age",file=csvname,append=T)
    write.table2(res$faa,title.tmp="F at age")

    write("\n# numbers at age",file=csvname,append=T)
    write.table2(res$naa,title.tmp="Numbers at age")

    write("\n# total and spawning biomass ",file=csvname,append=T)
    x <- rbind(colSums(res$ssb, na.rm=T),colSums(res$baa, na.rm=T),colSums(res$wcaa, na.rm=T))
    rownames(x) <- c("Spawning biomass","Total biomass","Catch biomass")
    write.table2(x,title.tmp="Total and spawning biomass")

    write("\n# YPR & SPR history ",file=csvname,append=T)
    get.SPR(res)$ysdata %>% rownames_to_column(var="year") %>%
                   as_tibble() %>% select(-"F/Ftarget") %>%
                   write_csv(path=csvname,append=T, col_names=TRUE)
  }

  if(!is.null(srres)){

    get_summary_ <- function(srres){
        as_tibble(srres$pars) %>% mutate(AICc   = srres$AICc,
                                         AIC    = srres$AIC,
                                         method = srres$input$method,
                                         type   = srres$input$SR,
                                         AR     = srres$input$AR,
                                         out.AR = srres$input$out.AR)
    }

    if("fit.SR" %in% class(srres)){
      write("\n# SR fit data",file=csvname,append=T)
      srres$input$SRdata %>% as_tibble() %>%  mutate(weight=srres$input$w) %>%
        write_csv(path=csvname,append=T,col_names=TRUE)
      write("\n# SR fit resutls",file=csvname,append=T)
      sr_summary <- get_summary_(srres)
      write_csv(sr_summary,path=csvname,append=T,
                col_names=TRUE)
    }
    if("fit.SRregime" %in% class(srres)){
      write("\n# SR fit data",file=csvname,append=T)
      srres$input$SRdata %>% as_tibble() %>%  mutate(weight=srres$input$w) %>%
        write_csv(path=csvname,append=T,col_names=TRUE)

      write("\n# SR fit resutls",file=csvname,append=T)
      tibble(AICc   =srres$AICc,
             AIC    =srres$AIC,
             method=srres$input$method,
             type  =srres$input$SR) %>%
        write_csv(path=csvname,append=T,col_names=TRUE)

      partable <- srres$regime_pars
      if(!is.null(srres$steepness)) partable <- partable %>% left_join(srres$steepness)
      # tentative
      write_csv(partable, path=csvname,append=T,col_names=TRUE)
    }
    if("SRfit.average" %in% class(srres)){
      write("\n# SR fit data",file=csvname,append=T)
      srres[[1]]$input$SRdata %>% as_tibble() %>%  mutate(weight=srres$input$w) %>%
        write_csv(path=csvname, append=T, col_names=TRUE)

      write("\n# SR fit resutls",file=csvname,append=T)
      sr_summary <- purrr::map_dfr(srres, function(x) get_summary_(x), .id="id")
      write_csv(sr_summary,path=csvname,append=T,
                col_names=TRUE)
    }
  }

  if(!is.null(msyres)){
    write("\n# MSY Reference points",file=csvname,append=T)
    write_csv(msyres$summary,path=csvname,append=T,
              col_names=TRUE)
  }

  tmpfunc <- function(fres, label=""){
    if(class(fres)=="future_new"){
      fres <- format_to_old_future(fres)
    }

    write(str_c("\n# future F at age",label), file=csvname,append=T)
    write.table2(apply(fres$faa,c(1,2),mean),title.tmp="Average future F at age")

    write(str_c("\n# future numbers at age",label), file=csvname,append=T)
    write.table2(apply(fres$naa,c(1,2),mean),title.tmp="Average future numbers at age")

    write(str_c("\n# future maturity at age",label), file=csvname,append=T)
    write.table2(apply(fres$maa,c(1,2),mean),title.tmp="Average maturity numbers at age")

    write(str_c("\n# future weight (for biomass) at age",label), file=csvname,append=T)
    write.table2(apply(fres$waa,c(1,2),mean),title.tmp="Average weight numbers at age")

    write(str_c("\n# future weight (for catch) at age",label), file=csvname,append=T)
    write.table2(apply(fres$waa.catch,c(1,2),mean),title.tmp="Average weight numbers at age")

    write(str_c("\n# future total biomass",label), file=csvname,append=T)
    make_summary_table(fres$vbiom,1,probs=ci.future) %>%
      rownames_to_column(var="year") %>%
      write_csv(path=csvname,append=TRUE, col_names = TRUE)

    write(str_c("\n# future total catch",label), file=csvname,append=T)
    make_summary_table(fres$vwcaa,1,probs=ci.future) %>%
      rownames_to_column(var="year") %>%
      write_csv(path=csvname,append=TRUE, col_names = TRUE)
  }

  if(!is.null(fres_current)){
    write("\n# future projection under F current",file=csvname,append=T)
    tmpfunc(fres_current, label="- Fcurrent")
  }

  if(!is.null(fres_HCR)){
    write("\n# future projection under HCR",file=csvname,append=T)
    tmpfunc(fres_HCR, label="- HCR")
  }

  if(!is.null(kobeII)){
    write("\n# Kobe II table",file=csvname,append=T)
    kobeII.table_name <- names(kobeII)
    for(i in 1:length(kobeII.table_name)){
      tmptable <- kobeII[kobeII.table_name[i]][[1]]
      if(nrow(tmptable)>0){
        write(str_c("\n# ",kobeII.table_name[i]),file=csvname,append=T)
        write_csv(tmptable,path=csvname,append=TRUE,
                  col_names = TRUE)
      }
    }
  }

  if(!is.null(kobe.ratio)){
    write("\n# Kobe ratio",file=csvname,append=T)
    kobe.ratio %>%
        write_csv(path=csvname,append=T, col_names=TRUE)
  }

  if(!is.null(other_tables)){
    for(i in seq_len(length(other_tables))){
      write(str_c("\n# ", names(other_tables)[i]), file=csvname,append=T)
      other_tables[[i]] %>%
        write_csv(path=csvname,append=T, col_names=TRUE)
    }
  }
}

#' csvファイルとしてまとめられた資源計算結果を読み込んでRのオブジェクトにする
#'
#' @param release.label 放流データを読み込みたい場合。ここで指定したラベルをつけて、年齢を行、年を列とするデータを入れる(catch at ageと同じ構造だが、年齢は必要な年齢だけ取り出した形で良い)
#' @param tfile 資源計算結果がまとめられたcsvファイルの名前
#' @param Pope  VPA計算時にどっちを使っているかここで設定する(TRUE or FALSE)。デフォルトはNULLで、その場合にはcaa,faa,naaの関係から自動判別するが、自動判別の結果が出力されるので、それをみて正しく判断されているか確認してください。
#' @param plus.group プラスグループを考慮するかどうか。こちらについても、NULLの場合にはfaaとnaaの関係から自動判別するが、結果を一応確認すること。
#' @param release_alive.label 放流魚のうち生き残って加入したものの尾数。年齢✕年の行列。データのない年齢は省略可。(release_all.labelとrelease_aliverate.labelが与えられる場合、release_alive.labelは与えないこと。加入年齢よりも高齢のデータは入力できるが、現時点では再生産関係や将来予測では考慮されない)
#' @param release_all.label 全放流尾数。年齢✕年の行列。データのない年齢は省略可。
#' @param release_ratealive.label 全放流尾数が生き残って加入したする比率。年齢✕年の行列。データのない年齢は省略可。
#'
#' @encoding UTF-8
#'
#'
#'
#' @export

read.vpa <- function(tfile,
                     caa.label="catch at age",
                     maa.label="maturity at age",
                     waa.label="weight at age",
                     waa.biomass.label="weight at age for biomass calculation",
                     waa.catch.label="weight at age for catch calculation",
                     M.label="M at age",
                     faa.label="fishing mortality at age",
                     Fc.label="Current F",
                     naa.label="numbers at age",
                     release_alive.label="release alive dat", # for older version "release dat"
                     release_all.label="release all dat",
                     release_ratealive.label="release alive rate dat",
                     Blimit=NULL,
                     Pope=NULL,
                     plus.group=NULL,
                     fc.year=NULL){

  tmpdata <- read.csv(tfile,header=F,as.is=F,colClasses="character")

  tmpfunc <- function(tmpdata,label,type=NULL){
    flags <- which(substr(tmpdata[,1],1,1)=="#")
    flags <- c(flags, nrow(tmpdata)+1)
    flag.name <- tmpdata[flags,1]
    flag.name <- gsub("#","",flag.name)
    flag.name <- gsub(" ","",flag.name)
    get.flag <- which(flag.name==gsub(" ","",label))
    {if(length(get.flag)>0){
      tdat <- tmpdata[(flags[get.flag]+1):(flags[min(which(flags[get.flag]<flags))]-1),]
      if(!is.null(type)){
        tdat <- tdat[,!apply(tdat=="",2,all)]
        tdat <- as.numeric(tdat)
      }
      else{
        tmp.names <- list(tdat[-1,1],tdat[1,-1])
        tdat <- tdat[,!apply(tdat=="",2,all)]
        tdat <- tdat[!apply(tdat=="",1,all),]
        tdat <- tdat[,!apply(is.na(tdat),2,all)]
        tdat <- tdat[!apply(is.na(tdat),1,all),]
        tdat <- sapply((tdat[-1,-1]),as.numeric)
        tmp.names <- lapply(tmp.names,function(x) x[x!=""])
        tmp.names <- lapply(tmp.names,function(x) x[!is.na(x)])
        if(is.null(dim(tdat))) dim(tdat) <- sapply(tmp.names,length)
        dimnames(tdat) <- tmp.names
        tdat <- as.data.frame(tdat)
      }
    }
      else{
        tdat <- NULL
      }}
    tdat
  }

  dres <- list()
  dres$naa <- tmpfunc(tmpdata,naa.label)
  dres$faa <- tmpfunc(tmpdata,faa.label)
  dres$Fc.at.age <- tmpfunc(tmpdata,Fc.label,type="Fc")
  #  dres$Fc.at.age <- dres$Fc.at.age[!is.na(dres$Fc.at.age)]
  dres$Fc.at.age <- dres$Fc.at.age[1:nrow(dres$naa)]
  #  if(length(dres$Fc.at.age)!=nrow(dres$naa)) stop("Dimension of Fc.at.age and numbers at age is differerent.")

  dres$input <- list()
  dres$input$dat <- list()
  dres$input$dat$maa <- tmpfunc(tmpdata,maa.label)
  dres$input$dat$caa <- tmpfunc(tmpdata,caa.label)
  dres$input$dat$M <- tmpfunc(tmpdata,M.label)
  dres$input$dat$waa <- tmpfunc(tmpdata,waa.label)
  if(is.null(dres$input$dat$waa)) dres$input$dat$waa <- tmpfunc(tmpdata,waa.biomass.label)
  dres$input$dat$waa.catch <- tmpfunc(tmpdata,waa.catch.label)
  if(is.null(dres$input$dat$waa.catch)) dres$input$dat$waa.catch <- dres$input$dat$waa

  # for release data (only release alive)
  release.old <- tmpfunc(tmpdata,"release dat")
  dres$input$dat$release.alive <- tmpfunc(tmpdata,release_alive.label)
  if(!is.null(release.old) && is.null(dres$input$dat$release.alive)) dres$input$dat$release.alive <- release.old
  if(!is.null(dres$input$dat$release.alive)) assertthat::assert_that(is.null(dres$input$dat$release.all) && is.null(dres$input$dat$release.aliverate))

  # for release data (with release all and rate)
  dres$input$dat$release.all <- tmpfunc(tmpdata,release_all.label)
  dres$input$dat$release.ratealive <- tmpfunc(tmpdata,release_ratealive.label)
  if(!is.null(dres$input$dat$release.all)){
    assertthat::assert_that(!is.null(dres$input$dat$release.ratealive),
                            is.null(dres$input$dat$release.alive),
                            all(dim(dres$input$dat$release.ratealive) == dim(dres$input$dat$release.all)))
    dres$input$dat$release.alive <- dres$input$dat$release.ratealive * dres$input$dat$release.all
  }

  # create ssb & baa data
  dres$ssb <- dres$input$dat$waa * dres$input$dat$maa * dres$naa
  dres$ssb <- as.data.frame(dres$ssb)

  dres$baa <- dres$input$dat$waa * dres$naa
  dres$baa <- as.data.frame(dres$baa)

  # create total catch
  dres$wcaa <- dres$input$dat$waa.catch * dres$input$dat$caa
  dres$wcaa <- as.data.frame(dres$wcaa)

  dres$Blimit <- Blimit

  ## catch at ageの計算時にpopeの近似式を使っているかどうか、通常は外から情報として与えてほしいところだが、与えられない場合、入力されたcaa,faa,naaの関係を見て、Popeで計算されているのかそうでないのかを判断してdres$input$Popeに入れる
  if(is.null(Pope)){
    caa.pope  <- dres$naa*(1-exp(-dres$faa))*exp(-dres$input$dat$M/2)
    diff.pope <- mean(unlist(dres$input$dat$caa/caa.pope),na.rm=T)

    faa <- dres$faa
    M <- dres$input$dat$M
    caa.bara <- dres$naa*faa/(faa+M)*(1-exp(-faa-M))
    diff.bara <- mean(unlist(dres$input$dat$caa/caa.bara),na.rm=T)

    if(abs(1-mean(diff.bara))>abs(1-mean(diff.pope))){
      dres$input$Pope <- TRUE

      cat("Pope is TRUE... OK? (mean difference=", 1-mean(diff.pope),")\n")
    }
    else{
      dres$input$Pope <- FALSE
      cat("Pope is FALSE... OK? (mean difference=", 1-mean(diff.bara),")\n")
    }
  }
  else{
    dres$input$Pope <- Pope
  }

  ## プラスグループを考慮しているかどうかを判別する
  if(is.null(plus.group)){
    plus.group <- detect_plus_group(dres)
    cat("Plus group is TRUE... OK?\n")
  }
  dres$input$plus.group <- plus.group

  if(is.null(dres$Fc.at.age) && !is.null(fc.year)) dres$Fc.at.age <- apply(dres$faa[,colnames(dres$faa)%in%fc.year],1,mean)

  # その他、他関数で必要になるVPAへのインプット
  dres$input$last.catch.zero <- FALSE
  class(dres) <- "vpa"

  return(dres)
}

#' tidy形式のVPAへのインプットデータをdata_dandlerへ渡す形式に変換する(暫定版)
#'
#' @export
#' @encoding UTF-8
#'

to_vpa_data <- function(x, label_name){
  vpa_data <- x %>% dplyr::filter(label==label_name) %>% as.data.frame()
  if(label_name!="abund"){
    rownames(vpa_data) <- str_c("X",vpa_data$year)
    vpa_data <- vpa_data %>% select(-year, -label)
    vpa_data$value <- vpa_data$fishery <- NULL
    vpa_data <- as.data.frame(t(vpa_data))
    rownames(vpa_data) <- str_sub(rownames(vpa_data),5,5)
  }
  else{
    vpa_data <- vpa_data %>%
      select(-starts_with("age_")) %>%
      pivot_wider(names_from=year, values_from=value)
    abund_name <- vpa_data$fishery
    vpa_data <- as.data.frame(vpa_data)
    vpa_data$label <- vpa_data$fishery <- NULL
    rownames(vpa_data) <- abund_name
    colnames(vpa_data) <- str_c("X",colnames(vpa_data))
  }
  vpa_data
}


#' @export
#' @encoding UTF-8
#'

get.stat <- function(fout,eyear=0,tmp.year=NULL, use_new_output=FALSE){

  if(isTRUE(use_new_output)){
    fout <- format_to_old_future(fout)
    col.target <- TRUE
  }
  else{
    col.target <- ifelse(fout$input$N==0,1,-1)
  }


  tmp <- as.numeric(fout$vssb[(nrow(fout$vssb)-eyear):nrow(fout$vssb),col.target])
  if(is.null(tmp.year)) tmp.year <- (nrow(fout$vwcaa)-eyear):nrow(fout$vwcaa)
  a <- data.frame("catch.mean"=mean(fout$vwcaa[tmp.year,col.target]),
                  "catch.sd"=sd(fout$vwcaa[tmp.year,col.target]),
                  "catch.geomean"=geomean(fout$vwcaa[tmp.year,col.target]),
                  "catch.median"=median(fout$vwcaa[tmp.year,col.target],na.rm=T),
                  "catch.L10"=quantile(fout$vwcaa[tmp.year,col.target],na.rm=T,probs=0.1),
                  "catch.H10"=quantile(fout$vwcaa[tmp.year,col.target],na.rm=T,probs=0.9),
                  "ssb.mean"=mean(fout$vssb[tmp.year,col.target]),
                  "ssb.sd"=sd(fout$vssb[tmp.year,col.target]),
                  "ssb.geomean"=geomean(fout$vssb[tmp.year,col.target]),
                  "ssb.median"=median(fout$vssb[tmp.year,col.target],na.rm=T),
                  "ssb.L10"=quantile(fout$vssb[tmp.year,col.target],na.rm=T,probs=0.1),
                  "ssb.H10"=quantile(fout$vssb[tmp.year,col.target],na.rm=T,probs=0.9),
                  "biom.mean"=mean(fout$vbiom[tmp.year,col.target]),
                  "biom.sd"=sd(fout$vbiom[tmp.year,col.target]),
                  "biom.geomean"=geomean(fout$vbiom[tmp.year,col.target]),
                  "biom.median"=median(fout$vbiom[tmp.year,col.target],na.rm=T),
                  "biom.L10"=quantile(fout$vbiom[tmp.year,col.target],na.rm=T,probs=0.1),
                  "biom.H10"=quantile(fout$vbiom[tmp.year,col.target],na.rm=T,probs=0.9),
                  "cbiom.mean"  = mean   (fout$vbiom_catch[tmp.year,col.target]),
                  "cbiom.sd"     =sd     (fout$vbiom_catch[tmp.year,col.target]),
                  "cbiom.geomean"=geomean(fout$vbiom_catch[tmp.year,col.target]),
                  "cbiom.median" =median (fout$vbiom_catch[tmp.year,col.target],na.rm=T),
                  "cbiom.L10"   =quantile(fout$vbiom_catch[tmp.year,col.target],na.rm=T,probs=0.1),
                  "cbiom.H10"   =quantile(fout$vbiom_catch[tmp.year,col.target],na.rm=T,probs=0.9),
                  "rec.mean"=mean(unlist(fout$naa[1,,])[tmp.year,col.target]),
                  "rec.sd"=sd(unlist(fout$naa[1,,])[tmp.year,col.target]),
                  "rec.geomean"=geomean(unlist(fout$naa[1,,])[tmp.year,col.target]),
                  "rec.median"=median(unlist(fout$naa[1,,])[tmp.year,col.target],na.rm=T),
                  "rec.L10"=quantile(unlist(fout$naa[1,,])[tmp.year,col.target],na.rm=T,probs=0.1),
                  "rec.H10"=quantile(unlist(fout$naa[1,,])[tmp.year,col.target],na.rm=T,probs=0.9),
                  #                    "lower.HSpoint"=lhs,
                  "Fref2Fcurrent"=fout$multi,
                  fmulti=fout$multi
  )
  a$U.mean <- a$catch.mean/a$cbiom.mean
  a$U.median <- a$catch.median/a$cbiom.median
  a$U.geomean <- a$catch.geomean/a$cbiom.geomean

  a$catch.CV <- a$catch.sd/a$catch.mean
  a$ssb.CV <- a$ssb.sd/a$ssb.mean
  a$biom.CV <- a$biom.sd/a$biom.mean
  a$rec.CV <- a$rec.sd/a$rec.mean

  #    Faa <- as.data.frame(t(fout$multi * fout$input$res0$Fc.at.age))
  # Faa <- as.data.frame(t(fout$multi * fout$currentF))
  if("finalmeanF" %in% names(fout)) Faa <- as.data.frame(t(fout$finalmeanF)) else Faa <- as.data.frame(t(fout$multi * fout$futureF))
  colnames(Faa) <- paste("F",dimnames(fout$naa)[[1]],sep="")
  res.stat1 <- cbind(a,Faa) # ここまで、get.stat

  tmpfunc_ <- function(x, nage, agename, label){
    x.mat <- matrix(0,nage,5)
    for(i in 1:nage){
      x.mat[i,1] <- mean   (x[i,tmp.year,col.target])
      x.mat[i,2] <- median (x[i,tmp.year,col.target])
      x.mat[i,3] <- geomean(x[i,tmp.year,col.target])
      x.mat[i,4:5] <- quantile(x[i,tmp.year,col.target],probs=c(0.1,0.9),na.rm=T)
    }
    x.mat <- as.numeric(x.mat)
    names(x.mat) <- c(paste(label,"-mean-A",agename,sep=""),
                       paste(label,"-median-A",agename,sep=""),
                       paste(label,"-geomean-A",agename,sep=""),
                       paste(label,"-L10-A",agename,sep=""),
                       paste(label,"-H10-A",agename,sep=""))
    x.mat
  }

  nage <- dim(fout$naa)[[1]]
  agename <- dimnames(fout$naa)[[1]]

  tb  <- fout$naa * fout$waa
  ctb <- fout$naa * fout$waa.catch
  ssb <- fout$naa * fout$waa *fout$maa
  tc  <- fout$wcaa

  tb.mat  <- tmpfunc_(tb,  nage, agename, "TB")
  ctb.mat <- tmpfunc_(ctb, nage, agename, "CTB")
  ssb.mat <- tmpfunc_(ssb, nage, agename, "SSB")
  tc.mat  <- tmpfunc_(tc,  nage, agename, "TC")

  res.stat2 <- as.data.frame(t(c(tb.mat,ctb.mat,tc.mat,ssb.mat)))
  res.stat  <- cbind(res.stat1,res.stat2) %>% as_tibble()
  return(res.stat)
}

get.stat3 <- get.stat


get.stat4 <- function(fout,Brefs,
                      refyear=c(2019:2023,2028,2038)){
  col.target <- ifelse(fout$input$N==0,1,-1)
  years <- as.numeric(rownames(fout$vwcaa))

  if(is.null(refyear)){
    refyear <- c(seq(from=min(years),to=min(years)+5),
                 c(min(years)+seq(from=10,to=20,by=5)))
  }

  catch.mean <- rowMeans(fout$vwcaa[years%in%refyear,col.target])
  names(catch.mean) <- str_c("Catch",names(catch.mean))
  catch.mean <- as_tibble(t(catch.mean))

  Btarget.prob <- rowMeans(fout$vssb[years%in%refyear,col.target]>Brefs$Btarget) %>%
    t() %>% as_tibble()
  names(Btarget.prob) <- str_c("Btarget_prob",names(Btarget.prob))

  #    Blow.prob <- rowMeans(fout$vssb[years%in%refyear,col.target]>Brefs$Blow) %>%
  #        t() %>% as_tibble()
  #    names(Blow.prob) <- str_c("Blow_prob",names(Blow.prob))

  Blimit.prob <- rowMeans(fout$vssb[years%in%refyear,col.target]<Brefs$Blimit) %>%
    t() %>% as_tibble()
  names(Blimit.prob) <- str_c("Blimit_prob",names(Blimit.prob))

  Bban.prob <- rowMeans(fout$vssb[years%in%refyear,col.target]<Brefs$Bban) %>%
    t() %>% as_tibble()
  names(Bban.prob) <- str_c("Bban_prob",names(Bban.prob))

  return(bind_cols(catch.mean,Btarget.prob,Blimit.prob,Bban.prob))
}

geomean <- function(x)
{
  ifelse(all(x > 0), exp(mean(log(x))), NA)
}

get.trace <- function(trace){
  trace <- trace  %>% as_tibble() %>%
    select(starts_with("TC-mean"),ssb.mean,fmulti,catch.CV) %>%
    mutate(label=as.character(1:nrow(.)))

  trace <- trace %>% gather(value=value,key=age,-label,-fmulti,-ssb.mean,-catch.CV) %>%
    mutate(age=str_extract(age, "(\\d)+")) %>%
    mutate(age=factor(age)) %>%
    mutate(age=fct_reorder(age,length(age):1))
  return(trace)
}


#' 列が年である行列に対して、年を指定するとその年のあいだの平均値(等)を返す関数
#'
#' @encoding UTF-8
#' @export
#'
#'

apply_year_colum <- function(mat,target_year,stat="mean"){
  if(target_year[1]<0){
    target_year <- rev(colnames(mat))[-target_year] %>%
      as.numeric() %>% sort()
  }
  mat <- as.data.frame(mat)
  apply(mat[as.character(target_year)],1,get(stat))
}


# VPAや将来予測の結果をtidyに ----
convert_df <- function(df,name){
  df %>%
    as_tibble %>%
    mutate(age = as.numeric(rownames(df))) %>%
    gather(key=year, value=value, -age, convert=TRUE) %>%
    group_by(year) %>%
    #        summarise(value=sum(value)) %>%
    mutate(type="VPA",sim="s0",stat=name)
}

#'
#' @export
#' @encoding UTF-8
convert_2d_future <- function(df, name, label="tmp"){
  df %>%
    as_tibble %>%
    mutate(year=rownames(df)) %>%
    gather(key=sim, value=value, -year, convert=TRUE) %>%
    mutate(year=as.numeric(year), stat=name, label=label)
}

#' future_vpaの結果オブジェクトのリストをtibble形式に変換する関数
#'
#' make_kobeII_tableに入れるオブジェクトを作る
#'
#' @param fout_list future_vpaの結果のオブジェクトのリスト
#'
#' @encoding UTF-8
#' @export
#'


convert_future_list_table <- function(fout_list,name_vector=NULL,beta_vector=NULL,label="tmp"){


  if(is.null(name_vector)){
    if(!is.null(names(fout_list))) name_vector <- names(fout_list)
    else name_vector <- 1:length(fout_list)
  }

  if(is.null(beta_vector)){
    beta_vector <- name_vector
  }

  res <- purrr::map_dfr(1:length(fout_list),
                 function(i){
                   convert_future_table(fout_list[[i]],label=name_vector[i]) %>%
                     rename(HCR_name=label) %>%
                     mutate(beta=beta_vector[i])
                 })
  return(res)
}

#' future_vpaの結果オブジェクトをtibble形式に変換する関数
#'
#' @param fout future_vpaの結果のオブジェクト
#'
#' @encoding UTF-8
#' @export
#'

convert_future_table <- function(fout,label="tmp"){

  if(class(fout)=="future_new") fout <- format_to_old_future(fout)

  U_table <- fout$vwcaa/fout$vbiom_catch
  if(is.null(fout$Fsakugen)) fout$Fsakugen <- -(1-fout$faa[1,,]/fout$currentF[1])
  if(is.null(fout$recruit))  fout$recruit <- fout$naa[1,,]

  ssb      <- convert_2d_future(df=fout$vssb,   name="SSB",     label=label)
  catch    <- convert_2d_future(df=fout$vwcaa,  name="catch",   label=label)
  cbiomass  <- convert_2d_future(df=fout$vbiom_catch,  name="cbiomass", label=label)
  biomass  <- convert_2d_future(df=fout$vbiom,  name="biomass", label=label)
  U_table  <- convert_2d_future(df=U_table,     name="U",       label=label)
  beta_gamma    <- convert_2d_future(df=fout$alpha,  name="beta_gamma",   label=label)
  Fsakugen <- convert_2d_future(df=fout$Fsakugen, name="Fsakugen",   label=label)
  recruit  <- convert_2d_future(df=fout$recruit, name="Recruitment",   label=label)
  if(!is.null(fout$Fratio)){
    Fratio <- convert_2d_future(df=fout$Fratio, name="Fratio",   label=label)
  }
  else{
    Fratio <- NULL
  }

  Fsakugen_ratio <- Fsakugen %>%
    mutate(value=value+1)
  Fsakugen_ratio$stat <- "Fsakugen_ratio"

  bind_rows(ssb,catch,biomass,cbiomass,beta_gamma,Fsakugen,Fsakugen_ratio,recruit, U_table, Fratio)
}


convert_vector <- function(vector,name){
  vector %>%
    as_tibble %>%
    mutate(year = as.integer(names(vector))) %>%
    mutate(type="VPA",sim="s0",stat=name,age=NA)
}

#' VPAの結果オブジェクトをtibble形式に変換する関数
#'
#' @param vpares vpaの結果のオブジェクト
#' @encoding UTF-8
#'
#'
#' @export

convert_vpa_tibble <- function(vpares,SPRtarget=NULL){

  if (is.null(vpares$input$dat$waa.catch)) {
    vpares$input$dat$waa.catch <- vpares$input$dat$waa
    if (class(vpares)=="sam") {
      total.catch <- colSums(vpares$caa*vpares$input$dat$waa,na.rm=T)
    } else {
      total.catch <- colSums(vpares$input$dat$caa*vpares$input$dat$waa,na.rm=T)
    }
  } else {
    total.catch <- colSums(vpares$input$dat$caa*vpares$input$dat$waa.catch,na.rm=T)
  }

  # ここでcbiomassを定義する(今後もbioamss, ssbを計算するときは極力ssb, biomを使わないようにする)
  ssb <- vpares$naa * vpares$input$dat$maa * vpares$input$dat$waa
  biomass <- vpares$naa * vpares$input$dat$waa
  cbiomass <- vpares$naa * vpares$input$dat$waa.catch
  U <- total.catch/colSums(cbiomass, na.rm=T)
  SSB <- convert_vector(colSums(ssb,na.rm=T),"SSB") %>%
    dplyr::filter(value>0&!is.na(value))
  Biomass <- convert_vector(colSums(biomass,na.rm=T),"biomass") %>%
    dplyr::filter(value>0&!is.na(value))
  cBiomass <- convert_vector(colSums(cbiomass,na.rm=T),"cbiomass") %>%
    dplyr::filter(value>0&!is.na(value))
  FAA <- convert_df(vpares$faa,"fishing_mortality") %>%
    dplyr::filter(value>0&!is.na(value))
  Recruitment <- convert_vector(colSums(vpares$naa[1,,drop=F]),"Recruitment") %>%
    dplyr::filter(value>0&!is.na(value))

  if(!is.null(SPRtarget)){
    if(is.null(vpares$input$dat$waa.catch)) waa.catch <- vpares$input$dat$waa
    else waa.catch <- vpares$input$dat$waa.catch
    Fratio <- purrr::map_dfc(1:ncol(vpares$naa),
                             function(i){
                               tmp <- !is.na(vpares$faa[,i])
                               calc_Fratio(faa=vpares$faa[tmp,i],
                                           maa=vpares$input$dat$maa[tmp,i],
                                           waa=vpares$input$dat$waa[tmp,i],
                                           M  =vpares$input$dat$M[tmp,i],
                                           waa.catch=waa.catch[tmp],
                                           SPRtarget=SPRtarget,
                                           plus_group=vpares$input$plus.group)
                             })
    colnames(Fratio) <- colnames(vpares$naa)
    Fratio <- convert_df(Fratio,"Fratio")
  }
  else{
    Fratio <- NULL
  }

  all_table <- bind_rows(SSB,
                         Biomass,
                         cBiomass,
                         convert_vector(U[U>0],"U"),
                         convert_vector(total.catch[total.catch>0],"catch"),
                         convert_df(vpares$naa,"fish_number"),
                         FAA,
                         convert_df(vpares$input$dat$waa,"weight"),
                         convert_df(vpares$input$dat$maa,"maturity"),
                         convert_df(vpares$input$dat$caa,"catch_number"),
                         convert_df(vpares$input$dat$M,  "natural_mortality"),
                         Recruitment,
                         Fratio)
}

#' fit.SRの結果をtibble形式に治す
#'
#' @param SR_result fit.SRの結果のオブジェクト
#' @encoding UTF-8
#'
#' @export
#'

convert_SR_tibble <- function(res_SR){
  if(class(res_SR)=="fit.SR"){
    resSRtibble <-bind_rows(tibble(value=as.numeric(res_SR$pars),type="parameter",name=names(res_SR$pars)),
                            res_SR$pred %>% mutate(type="prediction",name="prediction"),
                            res_SR$input$SRdata %>% as_tibble() %>%
                              mutate(type="observed",name="observed",residual=res_SR$resid))
    if(!is.null(res_SR$steepness)) resSRtibble<-bind_rows(resSRtibble,tibble(value=as.numeric(res_SR$steepness),type="parameter",name=names(res_SR$steepness)))
  }
  if(class(res_SR)=="fit.SRregime"){ # regimeあり
    resSR1 <- pivot_longer(res_SR$regime_pars,col=-regime) %>% mutate(type="parameter")
    resSR2 <- res_SR$pred %>% mutate(type="prediction",name="prediction")
    resSR3 <- res_SR$input$SRdata %>% as_tibble() %>%
                              mutate(type="observed",name="observed",residual=res_SR$regime_resid$resid)

    resSRtibble<-bind_rows(resSR1,resSR2,resSR3)

    if(!is.null(res_SR$steepness)) {
      for(j in 1:nrow(res_SR$steepness)){
        res_steepness <- res_SR$steepness[j,]
        res_steepness_tibble <- pivot_longer(res_steepness,col=-regime) %>% mutate(type="parameter")
        resSRtibble<- bind_rows(resSRtibble,res_steepness_tibble)
      }
    }
  }
  return(resSRtibble)
}

#' 管理基準値の表を作成する
#'
#' @param refs_base est.MSYから得られる管理基準値の表
#' @encoding UTF-8
#'
#' @export
#'

make_RP_table <- function(refs_base){
  #    require(formattable)
  #    require(tidyverse,quietly=TRUE)
  table_output <- refs_base %>%
    select(-RP_name) %>% # どの列を表示させるか選択する
    # 各列の有効数字を指定
    mutate(SSB=round(SSB,-floor(log10(min(SSB)))),
           SSB2SSB0=round(SSB2SSB0,2),
           Catch=round(Catch,-floor(log10(min(Catch)))),
           Catch.CV=round(Catch.CV,2),
           U=round(U,2),
           Fref2Fcurrent=round(Fref2Fcurrent,2)) %>%
    rename("管理基準値"=RP.definition,"親魚資源量"=SSB,"B0に対する比"=SSB2SSB0,
           "漁獲量"=Catch,"漁獲量の変動係数"=Catch.CV,"漁獲率"=U,"努力量の乗数"=Fref2Fcurrent)

  table_output  %>%
    # 表をhtmlで出力
    formattable::formattable(list(`親魚資源量`=color_bar("olivedrab"),
                                       `漁獲量`=color_bar("steelblue"),
                                       `漁獲率`=color_bar("orange"),
                                       `努力量の乗数`=color_bar("tomato")))

  #    return(table_output)

}

#' 管理基準値表から目的の管理基準値を取り出す関数
#'
#' @param refs_base est.MSYから得られる管理基準値の表
#' @param RP_name 取り出したい管理基準値の名前
#' @encoding UTF-8
#'
#' @export
#'

derive_RP_value <- function(refs_base,RP_name){
  refs_base[refs_base$RP.definition%in%RP_name,]
}


# kobe II matrix など、パフォーマンスを計算する関数 ----

#'
#' beta.simluationの結果などを読んで、kobeII talbeに整形する関数
#'
#' @param kobeII_data beta.simulationまたはconvert_future_list_tableの返り値
#' @param res_vpa VPAの結果
#' @param year.catch 平均漁獲量の表を出力する期間。その他year.ssb(平均親魚量), year.ssbtarget(SBtargetを上回る確率)、、なども同様
#'
#' @details
#' tidy形式になっているkobeII_dataにおいて、HCR_name, betaの列のラベルの組み合わせを一つの管理方式として、その管理方式ごとに少尉予測の結果を集計する
#'
#' @export
#'
#' @encoding UTF-8

make_kobeII_table <- function(kobeII_data,
                              res_vpa,
                              year.catch=(max(as.numeric(colnames(res_vpa$naa)))+1:10),
                              year.ssb=(max(as.numeric(colnames(res_vpa$naa)))+1:10),
#                              year.Fsakugen=(max(as.numeric(colnames(res_vpa$naa)))+1:10),
                              year.ssbtarget=(max(as.numeric(colnames(res_vpa$naa)))+1:10),
                              year.ssblimit=(max(as.numeric(colnames(res_vpa$naa)))+1:10),
                              year.ssbban=(max(as.numeric(colnames(res_vpa$naa)))+1:10),
                              year.ssbmin=(max(as.numeric(colnames(res_vpa$naa)))+1:10),
#                              year.ssbmax=(max(as.numeric(colnames(res_vpa$naa)))+1:10),
#                              year.aav=(max(as.numeric(colnames(res_vpa$naa)))+1:10),
#                              year.risk=(max(as.numeric(colnames(res_vpa$naa)))+1:10),
#                              year.catchdiff=(max(as.numeric(colnames(res_vpa$naa)))+1:10),
                              Btarget=0,
                              Blimit=0,
                              Bban=0){
  # 平均漁獲量
  (catch.mean <- kobeII_data %>%
     dplyr::filter(year%in%year.catch,stat=="catch") %>% # 取り出す年とラベル("catch")を選ぶ
     group_by(HCR_name,beta,year) %>%
     summarise(catch.mean=mean(value)) %>%  # 値の計算方法を指定(漁獲量の平均ならmean(value))
     # "-3"とかの値で桁数を指定
     spread(key=year,value=catch.mean) %>% ungroup() %>%
     arrange(HCR_name,desc(beta)) %>% # HCR_nameとbetaの順に並び替え
     mutate(stat_name="catch.mean"))

  # 平均親魚
  (ssb.mean <- kobeII_data %>%
      dplyr::filter(year%in%year.ssb,stat=="SSB") %>%
      group_by(HCR_name,beta,year) %>%
      summarise(ssb.mean=mean(value)) %>%
      spread(key=year,value=ssb.mean) %>% ungroup() %>%
      arrange(HCR_name,desc(beta)) %>% # HCR_nameとbetaの順に並び替え
      mutate(stat_name="ssb.mean"))

  # 平均親魚
  (biomass.mean <- kobeII_data %>%
      dplyr::filter(year%in%year.ssb,stat=="biomass") %>%
      group_by(HCR_name,beta,year) %>%
      summarise(biomass.mean=mean(value)) %>%
      spread(key=year,value=biomass.mean) %>% ungroup() %>%
      arrange(HCR_name,desc(beta)) %>% # HCR_nameとbetaの順に並び替え
      mutate(stat_name="biomass.mean"))

  # 親魚, 下5%
  (ssb.ci05 <- kobeII_data %>%
      dplyr::filter(year%in%year.ssb,stat=="SSB") %>%
      group_by(HCR_name,beta,year) %>%
      summarise(ssb.ci05=quantile(value,probs=0.05)) %>%
      spread(key=year,value=ssb.ci05) %>% ungroup() %>%
      arrange(HCR_name,desc(beta)) %>% # HCR_nameとbetaの順に並び替え
      mutate(stat_name="ssb.ci05"))

  # 親魚, 上5%
  (ssb.ci95 <- kobeII_data %>%
      dplyr::filter(year%in%year.ssb,stat=="SSB") %>%
      group_by(HCR_name,beta,year) %>%
      summarise(ssb.ci95=quantile(value,probs=0.95)) %>%
      spread(key=year,value=ssb.ci95) %>% ungroup() %>%
      arrange(HCR_name,desc(beta)) %>% # HCR_nameとbetaの順に並び替え
      mutate(stat_name="ssb.ci95"))

  # 1-currentFに乗じる値=currentFからの努力量の削減率の平均値(実際には確率分布になっている)
  ## (Fsakugen.table <- kobeII_data %>%
  ##     dplyr::filter(year%in%year.Fsakugen,stat=="Fsakugen") %>% # 取り出す年とラベル("catch")を選ぶ
  ##     group_by(HCR_name,beta,year) %>%
  ##     summarise(Fsakugen=round(mean(value),2)) %>%
  ##     spread(key=year,value=Fsakugen) %>% ungroup() %>%
  ##     arrange(HCR_name,desc(beta)) %>% # HCR_nameとbetaの順に並び替え
  ##     mutate(stat_name="Fsakugen"))

  # SSB>SSBtargetとなる確率
  ssbtarget.table <- kobeII_data %>%
    dplyr::filter(year%in%year.ssbtarget,stat=="SSB") %>%
    group_by(HCR_name,beta,year) %>%
    summarise(ssb.over=round(100*mean(value>Btarget))) %>%
    spread(key=year,value=ssb.over) %>%
    ungroup() %>%
    arrange(HCR_name,desc(beta))%>%
    mutate(stat_name="Pr(SSB>SSBtarget)")

  # SSB>SSBlimとなる確率
  ssblimit.table <- kobeII_data %>%
    dplyr::filter(year%in%year.ssblimit,stat=="SSB") %>%
    group_by(HCR_name,beta,year) %>%
    summarise(ssb.over=round(100*mean(value>Blimit))) %>%
    spread(key=year,value=ssb.over)%>%
    ungroup() %>%
    arrange(HCR_name,desc(beta))%>%
    mutate(stat_name="Pr(SSB>SSBlim)")

  # SSB>SSBbanとなる確率
  ssbban.table <- kobeII_data %>%
    dplyr::filter(year%in%year.ssbban,stat=="SSB") %>%
    group_by(HCR_name,beta,year) %>%
    summarise(ssb.over=round(100*mean(value>Bban))) %>%
    spread(key=year,value=ssb.over)%>%
    ungroup() %>%
    arrange(HCR_name,desc(beta))%>%
    mutate(stat_name="Pr(SSB>SSBban)")

  # SSB>SSBmin(過去最低親魚量を上回る確率)
  ssb.min <- min(unlist(colSums(res_vpa$ssb, na.rm=T)))
  ssbmin.table <- kobeII_data %>%
    dplyr::filter(year%in%year.ssbmin,stat=="SSB") %>%
    group_by(HCR_name,beta,year) %>%
    summarise(ssb.over=round(100*mean(value>ssb.min))) %>%
    spread(key=year,value=ssb.over)%>%
    ungroup() %>%
    arrange(HCR_name,desc(beta))%>%
    mutate(stat_name="Pr(SSB>SSBmin)")

  # SSB>SSBmax(過去最低親魚量を上回る確率)
  ## ssb.max <- max(unlist(colSums(res_vpa$ssb, na.rm=T)))
  ## ssbmax.table <- kobeII_data %>%
  ##   dplyr::filter(year%in%year.ssbmax,stat=="SSB") %>%
  ##   group_by(HCR_name,beta,year) %>%
  ##   summarise(ssb.over=round(100*mean(value>ssb.max))) %>%
  ##   spread(key=year,value=ssb.over)%>%
  ##   ungroup() %>%
  ##   arrange(HCR_name,desc(beta))%>%
  ##   mutate(stat_name="Pr(SSB>SSBmax)")

  # オプション: Catch AAV mean
  ## calc.aav <- function(x)sum(abs(diff(x)))/sum(x[-1])
  ## catch.aav.table <- kobeII_data %>%
  ##   dplyr::filter(year%in%year.aav,stat=="catch") %>%
  ##   group_by(HCR_name,beta,sim) %>%
  ##   dplyr::summarise(catch.aav=(calc.aav(value))) %>%
  ##   group_by(HCR_name,beta) %>%
  ##   summarise(catch.aav.mean=mean(catch.aav)) %>%
  ##   arrange(HCR_name,desc(beta))%>%
  ##   mutate(stat_name="catch.aav")

  ## # risk
  ## calc.aav2 <- function(x){
  ##     xx <- sum(x[-1]/x[-length(x)]<0.5,na.rm=T) # 0が2つ続く場合にNAが発生するがそこは計算から除く
  ##     return(xx)
  ## }

  ##   #  calc.aav2 <- function(x){ xx <- sum(x[-1]/x[-length(x)]<0.5); if(is.na(xx)) browser() else return(xx)}
  ## catch.risk <- kobeII_data %>%
  ##   dplyr::filter(year%in%year.risk,stat=="catch") %>%
  ##   group_by(HCR_name,beta,sim) %>%
  ##   dplyr::summarise(catch.aav=calc.aav2(value)) %>%
  ##   group_by(HCR_name,beta) %>%
  ##   summarise(value=mean(catch.aav>0)) %>%
  ##   arrange(HCR_name,desc(beta))%>%
  ##   mutate(stat_name="catch.risk")

##   bban.risk <- kobeII_data %>%
##     dplyr::filter(year%in%year.risk & stat=="SSB") %>%
##     group_by(HCR_name,beta,sim) %>%
##     dplyr::summarise(Bban.fail=sum(value<Bban)) %>%
##     group_by(HCR_name,beta) %>%
##     summarise(value=mean(Bban.fail>0)) %>%
##     arrange(HCR_name,desc(beta))%>%
##     mutate(stat_name="bban.risk")

##   blimit.risk <- kobeII_data %>%
##     dplyr::filter(year%in%year.risk,stat=="SSB") %>%
##     group_by(HCR_name,beta,sim) %>%
##     dplyr::summarise(Blimit.fail=sum(value<Blimit)) %>%
##     group_by(HCR_name,beta) %>%
##     summarise(value=mean(Blimit.fail>0)) %>%
##     arrange(HCR_name,desc(beta))%>%
##     mutate(stat_name="blimit.risk")

##   overfishing.risk <- kobeII_data %>%
##     dplyr::filter(year%in%year.risk,stat=="Fratio") %>%
##     group_by(HCR_name,beta,sim) %>%
##     dplyr::summarise(overfishing=sum(value>1)) %>%
##     group_by(HCR_name,beta) %>%
##     summarise(value=mean(overfishing>0)) %>%
##     arrange(HCR_name,desc(beta))%>%
##     mutate(stat_name="overfishing.risk")

##   redzone.risk1 <- kobeII_data %>%
##     dplyr::filter(year%in%year.risk,stat=="Fratio")

##   redzone.risk2 <- kobeII_data %>%
##     dplyr::filter(year%in%year.risk,stat=="SSB") %>%
##     mutate(Bratio=value/Btarget) %>%
##     mutate(Fratio=redzone.risk1$value) %>%
##     mutate(is.redzone=(Bratio < 1 & round(Fratio,3) > 1)) %>%
##     arrange(HCR_name,beta,sim)

## #  browser()

##   redzone.risk <- redzone.risk2 %>%
##     group_by(HCR_name,beta,sim) %>%
##      dplyr::summarise(sum.redzone=sum(is.redzone==TRUE)) %>%
##     group_by(HCR_name,beta) %>%
##     dplyr::summarise(value=mean(sum.redzone>0)) %>%
##     arrange(HCR_name,desc(beta)) %>%
##     mutate(stat_name="redzone.risk")

  ## if(!is.null(Bspecific)){
  ##   bspecific.risk <- kobeII_data %>%
  ##     dplyr::filter(year%in%year.risk,stat=="SSB") %>%
  ##     group_by(HCR_name,beta,sim) %>%
  ##     dplyr::summarise(Bspecific.fail=sum(value<Bspecific)) %>%
  ##     group_by(HCR_name,beta) %>%
  ##     summarise(value=mean(Bspecific.fail>0)) %>%
  ##     arrange(HCR_name,desc(beta))%>%
  ##     mutate(stat_name="bspecific.risk")
  ## }else{
  ##   bspecific.risk <- NA
  ## }

  # kobe statistics
  ## overssbtar <- kobeII_data %>%
  ##   dplyr::filter(year%in%year.risk,stat=="SSB") %>%
  ##   mutate(is.over.ssbtar= value > Btarget)
  ## overFtar <- kobeII_data %>%
  ##   dplyr::filter(year%in%year.risk,stat=="Fratio") %>%
  ##   mutate(is.over.Ftar= round(value,2) > 1)
  ## overssbtar$is.over.Ftar <- overFtar$is.over.Ftar

  ## kobe.stat <- overssbtar %>%
  ##     mutate("red"   =(is.over.ssbtar==FALSE) & (is.over.Ftar==TRUE),
  ##            "green" =(is.over.ssbtar==TRUE ) & (is.over.Ftar==FALSE),
  ##            "yellow"=(is.over.ssbtar==FALSE) & (is.over.Ftar==FALSE),
  ##            "orange"=(is.over.ssbtar==TRUE ) & (is.over.Ftar==TRUE)) %>%
  ##   group_by(HCR_name,beta,year) %>%
  ##   summarise(red.prob=mean(red,na.rm=T),
  ##             green.prob=mean(green,na.rm=T),
  ##             yellow.prob=mean(yellow,na.rm=T),
  ##             orange.prob=mean(orange,na.rm=T))  %>%
  ##     mutate(stat_name="kobe.stat")


  res_list <- list(catch.mean   = catch.mean,
                   ssb.mean         = ssb.mean,
                   biomass.mean         = biomass.mean,
                   ssb.lower05percent            = ssb.ci05,
                   ssb.upper95percent            = ssb.ci95,
                   prob.over.ssbtarget  = ssbtarget.table,
                   prob.over.ssblimit   = ssblimit.table,
                   prob.over.ssbban     = ssbban.table,
                   prob.over.ssbmin     = ssbmin.table)
#                   prob.over.ssbmax     = ssbmax.table,
                   ## catch.aav       = catch.aav.table,
                   ## kobe.stat       = kobe.stat,
                   ## catch.risk = catch.risk,
                   ## overfishing.risk = overfishing.risk,
                   ## redzone.risk = redzone.risk,
                   ## bban.risk = bban.risk,
                   ## blimit.risk = blimit.risk,
#                   bspecific.risk = bspecific.risk)
  return(res_list)

}


#'
#' 様々なベータやinput設定をもとに複数の将来予測を繰り返して結果をtibbleで返す関数
#'
#' @param finput future_vpaの返り値の$input
#' @param beta_vector betaを単純に変える場合のbetaのベクトル
#' @param year_beta_change betaを変更する年の範囲。NULLの場合には全部の年を変える。
#' @param datainput_setting_original make_future_dataに与えて引数を作り直す場合の make_future_dataの返り値.
#' @param datainput_setting_extra betaに加えて変えたい追加の設定。beta_vectorと同じ長さ分必要で、beta_vectorも一緒にかかる
#' @param label_name HCR_nameのラベルの名前. NULLの場合、beta_vectorが使われる
#' @param ncore 並列計算する場合のコア数。動くかどうか不明。
#'
#' @description
#' その他、calculate_all_pmに渡す引数
#'
#' @encoding UTF-8
#' @export
#'
#'

beta.simulation <- function(finput,beta_vector,
                            year.lag=0,type="old",year_beta_change=NULL,
                            datainput_setting_extra   =NULL,
                            datainput_setting_original=NULL,
                            label_name = NULL,
                            output_type="tidy", ncore = 1,
                            save_detail=rep(0,length(beta_vector)), ...){

  if(!is.null(datainput_setting_extra)) assertthat::assert_that(length(datainput_setting_extra)==length(beta_vector))
  if(!is.null(label_name)){
    assertthat::assert_that(length(label_name)==length(beta_vector))
  }
  else{
    label_name <- beta_vector
  }
  assertthat::assert_that(length(beta_vector)==length(save_detail))

  tb <- tb2 <- NULL
  future_year <- dimnames(finput$tmb_data$HCR_mat)[[1]]
  if(!is.null(year_beta_change)){
    year_column_beta_change <- future_year %in% year_beta_change
  }
  else{
    year_column_beta_change <- TRUE
  }

  res_list <- purrr::map(rep(NA, length(beta_vector)), function(x) x)

  if(ncore==1){
    for(i in 1:length(beta_vector)){
      if(type=="old"){
        stop("old function of future.vpa is not supported now")
      }
      else{
        finput$tmb_data$HCR_mat[year_column_beta_change,,"beta"] <- beta_vector[i]
        if(!is.null(finput$MSE_input_data)) finput$MSE_input_data$input$HCR_beta <- beta_vector[i]
        if(!is.null(datainput_setting_extra)){
          finput$tmb_data <- redo_future(datainput_setting_original, datainput_setting_extra[[i]], only_data=TRUE)$data
        }
        fres_base <- do.call(future_vpa,finput)
        if(save_detail[i]==1) res_list[[i]] <- fres_base
#        fres_base <- format_to_old_future(fres_base)
      }
      tmp <- convert_future_table(fres_base,label=label_name[i]) %>%
        rename(HCR_name=label)  %>% mutate(beta=beta_vector[i])
      tb <- bind_rows(tb,tmp)

      tmp <- calculate_all_pm(fres_base,...) %>%
          mutate(HCR_name=label_name[i], beta=beta_vector[i])
      tb2 <- bind_rows(tb2,tmp)
    }
  }
  else{
    library(foreach)
    cl <- parallel::makeCluster(ncore, type="FORK")
    doParallel::registerDoParallel(cl)

    tb <- foreach::foreach(i=1:length(beta_vector))%dopar%{
      finput$tmb_data$HCR_mat[year_column_beta_change,,"beta"] <- beta_vector[i]
      if(!is.null(finput$MSE_input_data)) finput$MSE_input_data$input$HCR_beta <- beta_vector[i]
      if(!is.null(datainput_setting_extra)){
        finput$tmb_data <- redo_future(datainput_setting_original, datainput_setting_extra[[i]], only_data=TRUE)$data
      }
      do.call(future_vpa,finput)
    }
    parallel::stopCluster(cl)

    res_list <- tb
    res_list[save_detail==0] <- NULL

    tb2 <- purrr::map_dfr(seq_len(length(tb)), function(i){
      calculate_all_pm(tb[[i]],...) %>%
        rename(HCR_name=label_name[i],beta=beta_vector[i])
    })

    tb <- purrr::map_dfr(seq_len(length(tb)), function(i){
      format_to_old_future(tb[[i]]) %>%
        convert_future_table(label=label_name[i]) %>%
        rename(HCR_name=label)  %>% mutate(beta=beta_vector[i])
    })


  }

  if(sum(save_detail)==0) return(lst(tb,tb2))
  else return(list(tb=tb, tb2=tb2, res_list=res_list))
}


#' MSYを達成するときの\%SPRを計算する(calc_future_perSPRのwrapper)
#'
#' @encoding UTF-8
#' @export
#'

calc_perspr <- function(...){
  calc_future_perSPR(...)
}

#' 将来予測の結果オブジェクトから生物パラメータを取り出して、その生物パラメータをベースに、Fvectorで与えたF at ageに対応するSPRを計算して返す
#'
#' @param fout 将来予測のアウトプット(finputがない場合)。future_vpaの結果はformat_to_old_future関数をかまさないと動かない。
#' @param res_vpa Popeの式を使うかどうか、plus_groupの設定のために利用。実際、漁獲量は計算していないので、不要といえば不要。is_popleとplus_groupが設定されていてもこちらを優先する。
#' @param is_pope res_vpaがない場合、popeの式を使うか
#' @param plus_group es_vpaがない場合、プラスグループを考慮するか
#' @param Fvector Fのベクトル
#' @param target.col 将来予測の何列目の年を取り出すか(NULLの場合、最後の年)
#' @param target.year 将来予測の何年目を取り出すか(年の名前)(NULLの場合、最後の年)
#' @param SPRtarget これを与えると、Fvectorを何倍すればここで指定した\%SPRと一致するかを返すようになる
#'
#' @encoding UTF-8
#' @export
#'

calc_future_perSPR <- function(fout=NULL,
                               res_vpa=NULL,
                               biopar=NULL,
                               Fvector,
                               is_pope=NULL,
                               plus_group=NULL,
                        Fmax=10,
                        target.col=NULL,
                        target.year=NULL,
                        SPRtarget=NULL,
                        SPR_unit="digit" # or "%"
){

  if(!is.null(res_vpa)){
    info_source <- "vpa"
    is_pope <- res_vpa$input$Pope
    plus_group <- res_vpa$input$plus.group
  }
  if(!is.null(fout)){
    info_source  <- "future"
    if(class(fout)=="future_new") fout <- format_to_old_future(fout)
  }
  if(!is.null(biopar))  info_source  <- "bio" # bioが優先される

  fout.tmp <- fout

  SPR_multi <- ifelse(SPR_unit=="%", 100, 1)

  # 将来予測結果が与えられた場合
  if(info_source=="future"){
    # シミュレーションが複数回ある場合には、その平均値を用いる
    if(is.null(target.col) && is.null(target.year)){
      waa.tmp       <- fout.tmp$waa      [,dim(fout.tmp$waa)      [[2]],] %>% apply(1,mean)
      waa.catch.tmp <- fout.tmp$waa.catch[,dim(fout.tmp$waa.catch)[[2]],] %>% apply(1,mean)
      maa.tmp       <- fout.tmp$maa      [,dim(fout.tmp$maa)      [[2]],] %>% apply(1,mean)
      M.tmp         <- fout.tmp$M        [,dim(fout.tmp$M)        [[2]],] %>% apply(1,mean)
    }
    else{
      # 年の範囲を指定する場合、年で平均してから、シミュレーション回数で平均する
      if(!is.null(target.year)){
        if(!is.list(target.year)){
          target.year.char <- as.character(target.year)
          waa.tmp       <- fout.tmp$waa      [,target.year.char,,drop=FALSE] %>% apply(c(1,3),mean) %>% apply(1,mean)
          waa.catch.tmp <- fout.tmp$waa.catch[,target.year.char,,drop=FALSE] %>% apply(c(1,3),mean) %>% apply(1,mean)
          maa.tmp       <- fout.tmp$maa      [,target.year.char,,drop=FALSE] %>% apply(c(1,3),mean) %>% apply(1,mean)
          M.tmp         <- fout.tmp$M        [,target.year.char,,drop=FALSE] %>% apply(c(1,3),mean) %>% apply(1,mean)
      }
        else{
        waa.tmp       <- fout.tmp$waa      [,as.character(target.year$waa),,drop=FALSE] %>% apply(c(1,3),mean) %>% apply(1,mean)
        waa.catch.tmp <- fout.tmp$waa.catch[,as.character(target.year$waa.catch),,drop=FALSE] %>% apply(c(1,3),mean) %>% apply(1,mean)
        maa.tmp       <- fout.tmp$maa      [,as.character(target.year$maa),,drop=FALSE] %>% apply(c(1,3),mean) %>% apply(1,mean)
        M.tmp         <- fout.tmp$M        [,as.character(target.year$M),,drop=FALSE] %>% apply(c(1,3),mean) %>% apply(1,mean)
      }
    }
    if(!is.null(target.col)){
      waa.tmp       <- fout.tmp$waa[,target.col,]       %>% apply(1,mean)
      waa.catch.tmp <- fout.tmp$waa.catch[,target.col,] %>% apply(1,mean)
      maa.tmp       <- fout.tmp$maa[,target.col,]       %>% apply(1,mean)
      M.tmp         <- fout.tmp$M[,target.col,]         %>% apply(1,mean)
    }
    }}

  if(info_source=="vpa"){ # 将来予測結果が与えられない場合にはVPA結果からもってくる
    if(!is.list(target.year)){
      target.year.char <- as.character(target.year)
      waa.tmp       <- res_vpa$input$dat$waa      [target.year.char] %>% apply(1,mean)
      maa.tmp       <- res_vpa$input$dat$maa      [target.year.char] %>% apply(1,mean)
      M.tmp         <- res_vpa$input$dat$M        [target.year.char] %>% apply(1,mean)
      if(!is.null(res_vpa$input$dat$waa.catch)){
        waa.catch.tmp <- res_vpa$input$dat$waa.catch[target.year.char] %>% apply(1,mean)
      }
      else{
        waa.catch.tmp <- waa.tmp
      }
    }
    else{
      waa.tmp       <- res_vpa$input$dat$waa      [as.character(target.year$waa)      ] %>% apply(1,mean)
      maa.tmp       <- res_vpa$input$dat$maa      [as.character(target.year$maa)      ] %>% apply(1,mean)
      M.tmp         <- res_vpa$input$dat$M        [as.character(target.year$M)    ] %>% apply(1,mean)
      if(!is.null(res_vpa$input$dat$waa.catch)){
        waa.catch.tmp <- res_vpa$input$dat$waa.catch[as.character(target.year$waa.catch)] %>% apply(1,mean)
      }
      else{
        waa.catch.tmp <- waa.tmp
      }
    }}

  if(info_source=="bio"){
    waa.tmp <- biopar$waa
    maa.tmp <- biopar$maa
    M.tmp <- biopar$M
    waa.catch.tmp <- biopar$waa.catch
  }

  # 緊急措置。本来ならどこをプラスグループとして与えるかを引数として与えないといけない
  # 現状で、すべてのカラムがゼロ=資源計算では考慮されていないセルとして認識されている
  allsumpars <- waa.tmp+waa.catch.tmp+maa.tmp+M.tmp
  waa.tmp <- waa.tmp[allsumpars!=0]
  waa.catch.tmp <- waa.catch.tmp[allsumpars!=0]
  maa.tmp <- maa.tmp[allsumpars!=0]
  M.tmp <- M.tmp[ allsumpars!=0]
  Fvector <- Fvector %>%  as.numeric()
  Fvector <- Fvector[allsumpars!=0 & !is.na(allsumpars)]

  # SPRを計算
  if(!is.null(SPRtarget)) SPRtarget_tmp <- SPRtarget/SPR_multi*100 else SPRtarget_tmp <- NULL
  tmp <- calc_Fratio(Fvector,waa=waa.tmp,maa=maa.tmp,M=M.tmp,SPRtarget=SPRtarget_tmp,
                     waa.catch=waa.catch.tmp,Pope=is_pope,
                     return_SPR=TRUE,plus_group=plus_group)
  if(is.null(SPRtarget))  return(ifelse(length(tmp)==1,1*SPR_multi,tmp$SPR_original/100*SPR_multi))
  else{
    tmp$SPR_original <- tmp$SPR_original/100*SPR_multi
    tmp$SPR_est <- tmp$SPR_est/100*SPR_multi
    tmp$SPR_target <- tmp$SPR_target/100*SPR_multi
    tmp$waa <- waa.tmp
    tmp$waa.catch <- waa.catch.tmp
    tmp$maa <- maa.tmp
    tmp$M <- M.tmp
    return(tmp)
  }
}

#' kobeIItable から任意の表を指名して取り出す
#'
#' @param kobeII_table \code{make_kobeII_table}の出力
#' @param name \code{kobeII_table}の要素名
#'
#' @encoding UTF-8
pull_var_from_kobeII_table <- function(kobeII_table, name) {
  table <- kobeII.table[[name]]
  table %>%
    dplyr::arrange(desc(beta)) %>%
    dplyr::select(-HCR_name, -stat_name)
}

#' kobeIItableから取り出した表を整形
#'
#' - 報告書に不要な列を除去する
#' - 単位を千トンに変換
#' @param beta_table \code{pull_var_from_kobeII_table}で取得した表
#' @param divide_by 表の値をこの値で除する.トンを千トンにする場合には1000
#' @param round TRUEなら値を丸める.漁獲量は現状整数表示なのでデフォルトはTRUE
format_beta_table <- function(beta_table, divide_by = 1, round = TRUE) {
  beta   <- beta_table %>%
    dplyr::select(beta) %>%
    magrittr::set_colnames("\u03B2") # greek beta in unicode
  values <- beta_table %>%
    dplyr::select(-beta) / divide_by
  if (round == TRUE) return(cbind(beta, round(values)))
  cbind(beta, values)
}

#' 値の大小に応じて表の背景にグラデーションをつける
#' @param beta_table \code{format_beta_table}で整形したβの表
#' @param color 表の背景となる任意の色
colorize_table <- function(beta_table, color) {
  beta_table %>%
    formattable::formattable(list(formattable::area(col = -1) ~
                                    formattable::color_tile("white", color)))
}

#' 表を画像として保存
#'
#' # @inheritParams \code{\link{formattable::as.htmlwidget}}
#' # @inheritParams \code{\link{htmltools::html_print}}
#' # @inheritParams \code{\link{webshot::webshot}}
#' @param table ファイルとして保存したい表
#' @examples
#' \dontrun{
#' your_table %>%
#'  export_formattable(file = "foo.png")
#' }
#' @export
export_formattable <- function(table, file, width = "100%", height = NULL,
                               background = "white", delay = 0.1) {
  widget <- formattable::as.htmlwidget(table, width = width, height = height)
  path   <- htmltools::html_print(widget, background = background, viewer = NULL)
  url    <- paste0("file:///", gsub("\\\\", "/", normalizePath(path)))
  webshot::webshot(url,
                   file = file,
                   selector = ".formattable_widget",
                   delay = delay)
}

#' kobeIItableから任意の表を取得し画像として保存
#'
#' @inheritParams \code{\link{pull_var_from_kobeII_table}}
#' @inheritParams \code{\link{format_beta_table}}
#' @inheritParams \code{\link{colorize_table}}
#' @inheritParams \code{\link{export_formattable}}
export_kobeII_table <- function(name, divide_by, color, fname, kobeII_table) {
  kobeII_table %>%
    pull_var_from_kobeII_table(name) %>%
    format_beta_table(divide_by = divide_by) %>%
    colorize_table(color) %>%
    export_formattable(fname)
}

#' β調整による管理効果を比較する表を画像として一括保存
#'
#' @inheritParams \code{\link{pull_var_from_kobeII_table}}
#' @param fname_ssb 「平均親魚量」の保存先ファイル名
#' @param fname_catch 「平均漁獲量」の保存先ファイル名
#' @param fname_ssb_above_target 「親魚量が目標管理基準値を上回る確率」の保存先ファイル名
#' @param fname_ssb_above_limit 「親魚量が限界管理基準値を上回る確率」の保存先ファイル名
#' @examples
#' \dontrun{
#' export_kobeII_tables(kobeII.table)
#' }
#' @export
export_kobeII_tables <- function(kobeII_table,
                                 fname_ssb = "tbl_ssb.png",
                                 fname_catch = "tbl_catch.png",
                                 fname_ssb_above_target = "tbl_ssb>target.png",
                                 fname_ssb_above_limit = "tbl_ssb>limit.png") {
  blue   <- "#96A9D8"
  green  <- "#B3CE94"
  yellow <- "#F1C040"

  purrr::pmap(list(name = c("ssb.mean", "catch.mean",
                            "prob.over.ssbtarget", "prob.over.ssblimit"),
                   divide_by = c(1000, 1000, 1, 1),
                   color = c(blue, green, yellow, yellow),
                   fname = c(fname_ssb, fname_catch,
                             fname_ssb_above_target, fname_ssb_above_limit)),
              .f = export_kobeII_table,
              kobeII_table = kobeII_table)
}




#'
#' 将来予測の結果のリストを入れると、代表的なパフォーマンス指標をピックアップする
#'
#' @param future_list future_vpaまたはfuture.vpaの返り値のリスト
#' @param res_vpa vpaの返り値
#' @param ABC_year 特にABCに注目したい年
#' @param is_MSE MSEの結果を使うかどうか。MSEの結果の場合には、ABCや加入尾数の真の値との誤差も出力する
#' @param indicator_year ABC_yearから相対的に何年後を指標として取り出すか
#' @param Btarget 目標管理基準値の値
#' @param Blimit 限界管理基準値の値
#' @param Bban 禁漁水準の値
#' @param type 出力の形式。"long"は縦長(ggplotに渡すときに便利)、"wide"は横長(数字を直接比較するときに便利)
#' @param biomass.unit 資源量の単位
#'
#' @export
#' @encoding UTF-8
#'

get_performance <- function(future_list,res_vpa,ABC_year=2021,
                            is_MSE=FALSE,
                            indicator_year=c(0,5,10),Btarget=0, Blimit=0, Bban=0,
                            type=c("long","wide")[1],biomass.unit=10000,...){

  future_list_original <- future_list
  future_list <- purrr::map(future_list,
                            function(x) if(class(x)=="future_new")
                              format_to_old_future(x) else x)

  if(is.null(names(future_list))) names(future_list) <- 1:length(future_list)

  future_tibble <- purrr::map_dfr(1:length(future_list),
                                  function(i) convert_future_table(future_list[[i]],
                                                                   label=names(future_list)[i]) %>%
                                    rename(HCR_name=label) %>% mutate(beta=NA))

  kobe_res <- make_kobeII_table(future_tibble,res_vpa,
                                year.ssb   = ABC_year+indicator_year,
                                year.catch = ABC_year+indicator_year,
                                year.ssbtarget = ABC_year+indicator_year,
                                year.ssblimit  = ABC_year+indicator_year,
                                year.ssbban=NULL, year.ssbmin=NULL, year.ssbmax=NULL,
                                year.aav = c(ABC_year,ABC_year-1),
                                Btarget= Btarget,
                                Blimit = Blimit,
                                Bban   = Bban)

  if(isTRUE(is_MSE)){
    error_table <- purrr::map_dfr(1:length(future_list_original), function(i){
      if("SR_MSE" %in% names(future_list_original[[i]])){
        plot_bias_in_MSE(future_list_original[[i]], out="stat") %>%
          dplyr::filter(year %in% (ABC_year + indicator_year)) %>%
          group_by(year, stat) %>%
          summarise(mean_error=mean(Relative_error_normal)) %>%
          mutate(HCR_name=names(future_list)[i])
      }
      else{
        NULL
      }
    })
    error_table <- error_table %>%
      gather(key=stat_name,value=value,-HCR_name,-year,-stat) %>%
      mutate(unit="",stat_category="推定バイアス") %>%
      mutate(stat_year_name=str_c(stat_category,year)) %>%
      ungroup(year) %>%
      mutate(year=as.character(year)) %>%
      mutate(stat_name=stat) %>%
      select(-stat)
  }
  else{
    error_table <- NULL
  }

  junit <- c("","十","百","千","万")[log10(biomass.unit)+1]

  stat_data <- tibble(stat_name=c("ssb.mean","catch.mean","Pr(SSB>SSBtarget)","Pr(SSB>SSBlim)",
                                  "catch.aav"),
                      stat_category=c("平均親魚量 ", "平均漁獲量 ", "目標上回る確率 ", "限界上回る確率 ",
                                      "漁獲量変動"))

  kobe_res <- purrr::map_dfr(kobe_res[c("ssb.mean", "catch.mean", "prob.over.ssbtarget",
                                        "prob.over.ssblimit", "catch.aav")],
                             function(x) x %>% select(-beta) %>%
                               gather(key=year,value=value,-HCR_name,-stat_name)) %>%
    mutate(value=ifelse(stat_name %in% c("ssb.mean", "catch.mean"), value/biomass.unit, value)) %>%
    mutate(unit =ifelse(stat_name %in% c("ssb.mean", "catch.mean"), str_c(junit, "トン"), "%")) %>%
    mutate(unit =ifelse(stat_name %in% c("catch.aav"), "", unit)) %>%
    left_join(stat_data) %>%
    mutate(stat_year_name=str_c(stat_category,year))

  kobe_res <- bind_rows(kobe_res,error_table)

  if(type=="wide"){
    kobe_res <- kobe_res  %>%
      select(-year, -stat_name, -stat_category) %>%
      spread(key=HCR_name,value=value) %>%
      select(2:ncol(.),1)
  }

  return(tibble::lst(kobe_res,error_table))
}

#'
#' 短期的将来予測における複数の管理方策のパフォーマンスを比較する表を出力する
#'
#' @param future_list 将来予測の結果のリスト
#' @param res_vpa VPAの結果
#' @param ... get_performanceで必要な引数
#'
#' @encoding UTF-8
#' @export
#'

compare_future_performance <- function(future_list,res_vpa,res_MSY,
                                       biomass.unit=1000,is_MSE=FALSE,...){
  perform_res <- get_performance(future_list=future_list, res_vpa=res_vpa,
                                 Btarget=derive_RP_value(res_MSY$summary,"Btarget0")$SSB,
                                 Blimit =derive_RP_value(res_MSY$summary,"Blimit0")$SSB,
                                 Bban   =derive_RP_value(res_MSY$summary,"Bban0")$SSB,
                                 type="long",is_MSE=is_MSE,biomass.unit=biomass.unit,...)

  g1_ssb0 <- perform_res$kobe_res %>% dplyr::filter(stat_name=="ssb.mean") %>%
    ggplot() +
    geom_bar(aes(x=HCR_name,y=value,fill=stat_category),stat="identity") +
    facet_wrap(stat_category~year,ncol=1)+coord_flip()+
    geom_label(aes(x=HCR_name,y=max(value)/2,label=str_c(round(value),unit)),
               alpha=0.5)+
    theme_SH() + theme(legend.position="top")+xlab("Senario") +
    guides(fill=guide_legend(title=""))+
    scale_fill_manual(values=c("lightblue"))

  g1_ssb <- g1_ssb0 +
    geom_hline(yintercept=derive_RP_value(res_MSY$summary,"Btarget0")$SSB/biomass.unit,
               col="#00533E")+
    geom_hline(yintercept=derive_RP_value(res_MSY$summary,"Blimit0")$SSB/biomass.unit,
               col="#edb918")

  g1_catch0 <- g1_ssb0 %+% dplyr::filter(perform_res$kobe_res, stat_name=="catch.mean")+
    scale_x_discrete(labels=rep("",4))+
    scale_fill_manual(values=c("lightgreen"))

  g1_catch <- g1_catch0 +
    geom_hline(yintercept=derive_RP_value(res_MSY$summary,"Btarget0")$Catch/10000,
               col="#00533E")+
    geom_hline(yintercept=rev(colSums(res_vpa$wcaa,na.rm=T))[1]/biomass.unit,
               col="gray",lty=2)

  g1_probtar <- g1_catch0 %+% dplyr::filter(perform_res$kobe_res, stat_name=="Pr(SSB>SSBtarget)")+
    scale_fill_manual(values=c("lightblue"))+
    geom_hline(yintercept=c(0,50,100),col="gray",lty=2)

  g1_problim <- g1_catch0 %+% dplyr::filter(perform_res$kobe_res, stat_name=="Pr(SSB>SSBlim)")+
    scale_fill_manual(values=c("gray"))+
    geom_hline(yintercept=c(0,50,100),col="gray",lty=2)

  g1_error_table <- perform_res$error_table  %>%
    ggplot() +
    geom_bar(aes(x=HCR_name,y=value,fill=stat_category),stat="identity") +
    facet_grid(year~stat_name)+coord_flip()+
    geom_label(aes(x=HCR_name,y=max(value)/2,label=str_c(round(value,2),unit)),
               alpha=0.5)+
    theme_SH() + theme(legend.position="top")+xlab("Senario") +
    guides(fill=guide_legend(title=""))+
    scale_fill_manual(values=c("lightpink"))

  g1_performance <- gridExtra::marrangeGrob(list(g1_ssb,g1_probtar,g1_catch,g1_problim),
                                            widths=c(1.3,1,1,1),nrow=1,ncol=4)

  list(g1_performance, g1_error_table, perform_res)
  #ggsave(g1_performance, filename="g1_performance.png" ,path=output_folder,
  #           width=20, height= 10)

  #    ggsave(g1_error_table, filename="g1_error_table.png" ,path=output_folder,
  #       width=15, height= 8)

}


#'
#' calculate F/Ftarget based on F_\%SPR multiplier
#'
#' @param faa F at age
#' @param waa weight at age
#' @param maa maturity at age
#' @param M natural morality at age
#' @param SPRtarget target SPR (NULLの場合には最適化しない)
#' @param return_SPR return SPR as well as Fratio
#' @param plus_group プラスグループを考慮するかどうか
#' @param max.age SPR計算を打ち切る最大の年。デフォルトはInf
#'
#' もともとのF at ageの最大がexp(-7)よりも小さい場合にはFratio=0となる。一方で、F at ageをすごく大きくしても指定されたSPRを実現できないような場合のFratioの上限値を50とする。
#'
#' @export
#' @encoding UTF-8
#'


calc_Fratio <- function(faa, waa, maa, M, SPRtarget=30, waa.catch=NULL,Pope=TRUE, return_SPR=FALSE, plus_group=TRUE, max.age=Inf){

  if(plus_group==FALSE) max.age  <- length(faa)

  calc.rel.abund2_ <- function(sel,Fr){
    calc.rel.abund(sel,Fr,na=length(faa),M=M, waa=waa, waa.catch=waa.catch,
                     min.age=1,max.age=max.age,Pope=Pope,ssb.coef=0,maa=maa)
  }

  tmpfunc <- function(x,SPR0=0,...){
    SPR_tmp <- calc.rel.abund2_(faa,exp(x))$spr %>% sum()
    sum(((SPR_tmp/SPR0*100)-SPRtarget)^2)
  }

  if(max(faa, na.rm=T)<exp(-7)){ return(0) }

  else{
    tmp <- !is.na(faa)
    SPR0 <- calc.rel.abund2_(faa,0)$spr %>% sum()
    SPR_original <- calc.rel.abund2_(faa,1)$spr %>% sum()
    SPR_original <- SPR_original/SPR0*100
    if(!is.null(SPRtarget)){
        opt_res <- optimize(tmpfunc,interval=c(-7,log(50)),SPR0=SPR0)
        SPR_est <- calc.rel.abund2_(faa,exp(opt_res$minimum))$spr %>% sum()
        SPR_est <- SPR_est/SPR0 * 100
#        if(abs(SPR_est-SPRtarget)>0.01) {return(NA)}
        Fratio <- 1/exp(opt_res$minimum)
    }
    else{
      SPR_est <- SPR_original
      Fratio <- 1
    }

    if(isTRUE(return_SPR)){
        list(Fratio=Fratio, SPR_est=SPR_est, SPR_target=SPRtarget, SPR_original=SPR_original)
    }
    else{
        return(Fratio)
    }
  }
}


#'
#' 使うフォルダ名を与えると一連の結果の関数を読み込む関数
#'
#' @export
#' @encoding UTF-8
#'

load_folder <- function(folder_name){

  tmpfunc <- function(folder_name, file_name){
    if(isTRUE(file.exists(str_c(folder_name,"/",file_name)))){
      if(isTRUE(str_detect(file_name, pattern=".rda"))){
        a <- load(str_c(folder_name,"/",file_name))
        a <- get(a)
      }
      if(isTRUE(str_detect(file_name, pattern=".csv"))){
        a <- read_csv(str_c(folder_name,"/",file_name))
      }
    }
    else{
      a <- NA
    }
    return(a)
  }

  file_name <- c("res_MSY.rda","res_SR.rda","res_future_0.8HCR.rda","kobeII.table.rda","model_selection.csv")

  res_all <- list()
  for(i in 1:length(file_name)){
    res_all[[i]] <- purrr::map(folder_name, function(x) tmpfunc(x,file_name[i]))
    res_all[[i]] <- res_all[[i]][!is.na(res_all[[i]])]
  }
  names(res_all) <- file_name

  res_all$res_vpa <- purrr::map(res_all$res_MSY.rda, function(x) if(!is.na(x)) x$res_vpa else NA)
  res_all$res_vpa <- res_all$res_vpa[!is.na(res_all$res_vpa)]
  invisible(res_all)
}


#' Make Kobe ratio (?) for frasyr_tool
#'
#' @param result_vpa object createb dy vpa()
#' @param result_msy object created by script 1do_MSYest.R of SC meeting
#' @return A tibble object
#' @export
make_kobe_ratio <- function(result_vpa, result_msy) {

  assertthat::assert_that(
    assertthat::has_name(result_vpa, c("ssb")),
    assertthat::has_name(result_msy, c("summary")),
    assertthat::has_name(result_msy$summary, c("perSPR"))
  )

  return_kobe_ratio <- function() {
    tibble::tibble(year   = get_year(),
                   Fratio = get_f_ratio(),
                   Bratio = calc_b_ratio()) %>%
      ad_hoc_filter()
  }

  get_year <- function() {
    colnames(result_vpa$ssb)
  }

  get_f_ratio <- function() {
    target_spr  <- derive_RP_value(result_msy$summary,"Btarget0")$perSPR * 100
    spr_history <- get.SPR(result_vpa,
                           target_spr, Fmax = 7)

    assertthat::assert_that(
      assertthat::validate_that(is.list(spr_history)),
      assertthat::has_name(spr_history, "ysdata"),
      assertthat::has_name(spr_history$ysdata, "F/Ftarget"),
      assertthat::are_equal(rownames(spr_history$ysdata),
                            colnames(result_vpa$ssb))
    )

    force(spr_history$ysdata$`F/Ftarget`)
  }


  calc_b_ratio <- function() {
    ssb        <- colSums(result_vpa$ssb, na.rm=T)
    target_ssb <- derive_RP_value(result_msy$summary,"Btarget0")$SSB

    force(ssb / target_ssb)
  }

  ad_hoc_filter <- function(koberatio) {
    dplyr::filter(koberatio, !is.na(Bratio)) # When does Bratio become 'NA'?
  }

  return_kobe_ratio()
}

#' Source specific lines in an R file
#'
#' @param file character string with the path to the file to source.
#' @param lines numeric vector of lines to source in \code{file}.
#'
#' @export

source_lines <- function(file, lines, encoding="UTF-8",...){
    source(textConnection(readLines(file, encoding=encoding)[lines]),...)
}

#' re-calculate projection with different arguments
#'
#' @param data_future make_future_dataの返り値。これ全体でなくて、$inputのみでもOK
#'
#' @export

redo_future <- function(data_future, input_data_list, SR_sd=NULL, SR_b=NULL, only_data=FALSE,is_regime=(class(data_future$input$res_SR)=="fit.SRregime"), ...){

  if("input" %in% names(data_future)) input_data <- data_future$input else input_data <- data_future

  if(! all(names(input_data_list) %in% names(input_data))) stop("names of input_data_list is invalid!")

  input_data[names(input_data_list)] <- input_data_list

  if(!is.null(SR_sd)){
    if(is_regime){
      cat("This is regime future\n")
      input_data$res_SR$regime_pars$sd[] <- SR_sd
      input_data$res_SR$pars$sd[] <- SR_sd
    }
    else{
      input_data$res_SR$pars$sd[] <- SR_sd
    }
    if(input_data$resid_type!="lognormal") warning("SR_sd=0 cannot be used in nonparametric recruitment")
  }

  if(!is.null(SR_b)){
    if(is_regime){
      cat("This is regime future\n")
      input_data$res_SR$regime_pars$b <- SR_b
    }
    else{
      input_data$res_SR$pars$b <- SR_b
    }
  }

  future_data <- safe_call(make_future_data,input_data)
  if(only_data==TRUE) return(future_data) else future_vpa(future_data$data,...)
}

#'
#' 将来予測においてFが非常に小さい場合には決定論的予測と確率論的予測の平均がほぼ一致するかを確認するための関数
#'
#' HSを仮定していて、折れ点の下側から将来予測するような場合には、将来予測の最初のほうは一致しない。マッチングする年代などを工夫する必要がある
#'
#' @export
#'

test_sd0_future <- function(data_future,nsim=NULL,nyear=10,future_range=NULL,...){

  is_regime <- !is.null(data_future$input$regime_shift_option)
  {if(is_regime){
    sd_org <- data_future$input$res_SR$regime_pars %>%
      dplyr::filter(regime==data_future$input$regime_shift_option$future_regime) %>%
      select(sd)
  }
  else{
    sd_org <- data_future$input$res_SR$pars$sd
  }}

  # determine sample size
  tol <- 0.007
  if(is.null(nsim)) nsim <- round((0.5/tol)^2) # 1%以下(0.7%)の誤差が期待されるnsim
  cat("nsim for checking sd=0:",nsim,"\n")

  # run 2 funture projections
  res1 <- redo_future(data_future, list(nsim=nsim, nyear=nyear), multi_init=0.01, ...)
  res2 <- redo_future(data_future, list(nsim=2   , nyear=nyear), multi_init=0.01, SR_sd=0, ...)

  a <- try(compare_future_res12(res1,res2,future_range=future_range))

  cat("* Fをなるべく小さくした場合の将来予測において、決定論的予測と確率的予測の平均値がほぼ同じになるか?(ここでnoが出る場合にはバグの可能性があるので管理者に連絡してください)(モデル平均・バックワードリサンプリングの場合にはnoになっちゃいます(今後改善)): ",ifelse(class(a)=="try-error", "not ","OK\n"))

  return(lst(res1,res2,a))
}

compare_future_res12 <- function(res1,res2,tol=0.01, future_range=NULL){
  nyear <- dim(res1$naa)[[2]]
  if(is.null(future_range)) future_range <- res1$input$tmb_data$start_random_rec_year:nyear
  if(dim(res1$naa)[[3]]==dim(res2$naa)[[3]]){
      mean_difference_in_naa <- mean(abs(1-res1$naa[,future_range,]/res2$naa[,future_range,]),na.rm=T)
      mean_difference_in_wcaa <- mean(abs(1-res1$wcaa[,future_range,]/res2$wcaa[,future_range,]),na.rm=T)
  }
  else{
      mean_difference_in_naa <- 1-mean(apply(res1$naa[,future_range,],c(1,2),mean,na.rm=TRUE)/
                                     apply(res2$naa[,future_range,],c(1,2),mean,na.rm=TRUE),na.rm=T)
      mean_difference_in_wcaa <- 1-mean(apply(res1$wcaa[,future_range,],c(1,2),mean,na.rm=TRUE)/
                                      apply(res2$wcaa[,future_range,],c(1,2),mean,na.rm=TRUE),na.rm=T)
  }

  cat("mean_difference in naa=", mean_difference_in_naa,"\n")
  expect_equal(mean_difference_in_naa,0,tol=tol)

  cat("mean_difference in wcaa=", mean_difference_in_wcaa,"\n")
  expect_equal(mean_difference_in_wcaa,0,tol=tol)
}

#'
#' MSEの計算が正しいかを確認するための関数
#'
#' 1. sd=0の場合の結果の比較
#' 2. MSE_nsim=2にしてMSE_sd=0にしても良いかどうか
#'
#' @export
#'

check_MSE_sd0 <- function(data_future, data_MSE=NULL, nsim_for_check=10000, tol=c(0.01,0.01,0.01)){

  data_future_sd0 <- redo_future(data_future,list(nyear=5,nsim=5),only_data=TRUE,SR_sd=0)
  data_future     <- redo_future(data_future,list(nyear=5,nsim=5),only_data=TRUE)
  data_future_10000 <- redo_future(data_future,list(nyear=5,nsim=nsim_for_check),only_data=TRUE)

  if(!is.null(data_MSE)) data_MSE <- redo_future(data_MSE,list(nyear=5,nsim=5),only_data=TRUE)
  else data_MSE <- data_future

    # check MSE program is correct?
  res1 <- future_vpa(tmb_data=data_future_sd0$data,
                         optim_method="none",multi_init = 1,SPRtarget=0.3,
                         do_MSE=TRUE, MSE_input_data=data_future_sd0,MSE_nsim=2)
  res2 <- future_vpa(tmb_data=data_future_sd0$data,
                     optim_method="none",multi_init = 1,SPRtarget=0.3,
                     do_MSE=FALSE, MSE_input_data=data_future_sd0,MSE_nsim=2)
  a1 <- try(compare_future_res12(res1,res2,tol=tol[1]))
  cat("* 真の個体群動態のSD=0のとき, MSEした結果と単純シミュレーションの結果が一致するか?(加入の残差にリサンプリングを使っている場合にはSD=0でも決定論的予測にならないので一致しなくても大丈夫。対数正規分布残差でここがOKにならない場合にはバグの可能性があるので管理者に連絡してください) ",ifelse(class(a1)=="try-error", "not ",""),"OK\n")

  # check sd in MSE=0 is OK?
  res1.time <- system.time(
    res1 <- future_vpa(tmb_data=data_future$data,
                     optim_method="none",
                     multi_init = 1,SPRtarget=0.3,
                     do_MSE=TRUE, MSE_input_data=data_MSE,MSE_nsim=nsim_for_check))

  res2.time <- system.time(
    res2 <- future_vpa(tmb_data=data_future$data,
                     optim_method="none",
                     multi_init = 1,SPRtarget=0.3,
                     do_MSE=TRUE, MSE_input_data=data_MSE,MSE_nsim=2,MSE_sd=0))
  a2 <- try(compare_future_res12(res1,res2,tol=tol[2]))
  cat("* MSEしたとき、ABC算定年までの加入のSDを0として、決定論的な漁獲量をABCとしても、nsim_for_check回分の確率計算の平均漁獲量をABCした場合と同じになるか?(OKであればMSE_sd=0にして計算時間を短縮できる、加入の残差にリサンプリングを使っている場合には同じにならない:MSE_sd=0は使えない) ",ifelse(class(a2)=="try-error", "not ",""),"OK\n")
  cat("res1.time=",res1.time,"res2.time=",res2.time,"\n")

  # first year catch
  res1 <- future_vpa(tmb_data=data_future_10000$data,
                     optim_method="none",multi_init = 1,SPRtarget=0.3,
                     do_MSE=FALSE,MSE_input_data=data_future_10000)

  res2 <- future_vpa(tmb_data=data_future_10000$data,
                     optim_method="none",multi_init = 1,SPRtarget=0.3,
                     do_MSE=TRUE, MSE_nsim=2, MSE_sd=0, MSE_input_data=data_future_10000)
  # ここのtorelanceはそんなに高くない
  a3 <- try(expect_equal(mean(get_wcatch(res1)["2019",])/
               mean(get_wcatch(res2)["2019",]),
               1,tol=tol[3]))
  cat("* HCRを導入する最初の年のABCは通常の将来予測の平均漁獲量と、MSEを十分回数実施したときの平均漁獲量と一致するはず。それぞれnsim_for_check回数分計算した場合、一致するか?(上の2つが通っている場合、ここもOKになるはず。OKにならなかったらバグの可能性があるので管理者に連絡してください): ",ifelse(class(a3)=="try-error", "not ",""),"OK\n")
  return(lst(a1,a2,a3,res1,res2))
}


#'
#'
#'  @export
#'

take_interval <- function(prob,target){
    x <- which(abs(diff(sign(prob-target)))>0)
    c(x,x+1)
}


#'
#' 将来予測やVPAの結果から生物パラメータをとりだす
#'
#' @param res_obj VPAか将来予測の結果のオブジェクト。どちらでも良い。
#' @param derive_year 生物パラメータとF at ageを取り出す期間(年の名前で指定)
#' @param stat 取り出した期間のパラメータをここで指定する関数で処理する。基本は平均する(mean)
#'
#' 将来予測結果を入れる場合には複数のシミュレーション、年の間の結果をすべてstatする
#'
#' @export
#'

derive_biopar <- function(res_obj=NULL, derive_year=NULL, stat=mean){

  derive_year <- as.character(derive_year)

  if(!is.null(res_obj$input$dat$caa)){
    res_obj$input$dat$faa <- res_obj$faa
    derive_char <- c("M","waa","maa","faa")
    if(!is.null(res_obj$input$dat$waa.catch)) derive_char <- c(derive_char,"waa.catch")
    bio_par <- purrr::map_dfc(res_obj$input$dat[derive_char],
                              function(x) apply(x[,derive_year,drop=F],1,stat))
    if(!is.null(res_obj$input$dat$waa.catch)) bio_par$waa.catch <- bio_par$waa
  }

  if(class(res_obj)[[1]]=="future"||class(res_obj)[[1]]=="future_new"){
    bio_list <- res_obj[c("waa","faa")]
    if(is.null(res_obj$maa)) bio_list$maa <- res_obj$input$tmb_data$maa else bio_list$maa <- res_obj$maa
    bio_list$M <- res_obj$input$tmb_data$M
    if(is.null(bio_list$M)) bio_list$M <- res_obj$M
    if(!is.null(res_obj$waa_catch_mat)) bio_list$waa.catch <- res_obj$waa_catch_mat else bio_list$waa.catch <- res_obj$waa.catch
    bio_par <- purrr::map_dfc(bio_list,
                   function(x) apply(x[,derive_year,,drop=F],1,stat))
  }

  bio_par <- bio_par[apply(bio_par,1,sum)!=0,]
  bio_par <- bio_par[!is.na(apply(bio_par,1,sum)),]
  return(bio_par)
}


#' 与えられた個体群動態でプラスグループが考慮されているかどうか
#' @param dres VPAの結果
#'
#' @export

detect_plus_group <- function(dres){
  naa2 <- dres$naa[,2]
  plus_age <- max(which(!is.na(naa2)))
  naa2_plus <- calc_forward(naa=dres$naa,faa=dres$faa,M=dres$input$dat$M,t=1,plus_age=plus_age,plus_group=TRUE)[,2]
  naa2_noplus <- calc_forward(naa=dres$naa,faa=dres$faa,M=dres$input$dat$M,t=1,plus_age=plus_age,plus_group=FALSE)[,2]
  if(sum((naa2-naa2_plus)^2,na.rm=T)<sum((naa2-naa2_noplus)^2,na.rm=T)) plus.group <- TRUE else plus.group <- FALSE
  return(plus.group)
}


#' future_vpaの返り値に要約統計量を加えるための内部関数
#'
#' @param res_future future_futureの返り値
#' @param target 指定するとtargetで指定した列の値のみが抽出される。NULLの場合は全シミュレーションの平均値
#'
#' @export
#'

derive_future_summary <- function(res_future, target=NULL){

  assertthat::assert_that(class(res_future) == "future_new")

  if(is.null(target)){
    tmpfunc <- function(x, fun=mean) apply(x,1,fun)
  }
  if(!is.null(target)){
    tmpfunc <- function(x) x[,target]
  }

  Fmean <- apply(res_future$faa,c(2,3),sum)

  tibble(
    year    = as.numeric(dimnames(res_future$SR_mat[,,"ssb"])[[1]]),
    SSB     = tmpfunc(res_future$SR_mat[,,"ssb"]),
    biomass = tmpfunc(res_future$SR_mat[,,"biomass"]),
    cbiomass = tmpfunc(res_future$SR_mat[,,"cbiomass"]),
    recruit = tmpfunc(res_future$SR_mat[,,"recruit"]),
    intercept = tmpfunc(res_future$SR_mat[,,"intercept"]),
    deviance = tmpfunc(res_future$SR_mat[,,"deviance"]),
    deviance_sd = tmpfunc(res_future$SR_mat[,,"deviance"],fun=sd),
    catch   = tmpfunc(res_future$HCR_realized[,,"wcatch"]),
    beta    = tmpfunc(res_future$HCR_mat[,,"beta"]),
    Blimit  = tmpfunc(res_future$HCR_mat[,,"Blimit"]),
    Bban    = tmpfunc(res_future$HCR_mat[,,"Bban"]),
    beta_gamma = tmpfunc(res_future$HCR_realized[,,"beta_gamma"]),
    Fmean      = tmpfunc(Fmean),
    Fratio     = tmpfunc(res_future$HCR_realized[,,"Fratio"]))
}


#'
#' F at ageをVPAの結果から\%SPRで変換したりするための関数
#'
#' @param res_vpa VPAの結果オブジェクト(Popeの設定やF at ageをこちらからとってくる)
#' @param data_future 将来予測のためのデータ(生物パラメータを将来予測期間から撮ってくる場合に必要)
#' @param faa_vector 漁獲圧を代表するベクトル
#' @param faa_vector_year 漁獲圧を取り出すときの年の範囲(VPA期間限定)
#' @param faa_bio_year 漁獲圧をSPRに換算するときに生物パラメータを取り出す年の範囲(data_futureがある場合将来予測年も指定可能。生物パラメータが密度によって変わる場合)。list(waa = 2014:2018, waa.catch = 2014:2018, maa = 2016:2018,M   = 2014:2018)とすると生物パラメータによって異なる期間の指定も可能。下のsaa_bio_yearも同様
#' @param faa_bio 漁獲圧をSPRに換算するとき用いる生物パラメータそのものtibble(waa=waa, maa=maa, M=M)として指定する
#' @param saa_vector 選択率を代表するベクトル
#' @param saa_vector_year 選択率を取り出すときの年の範囲(VPA期間限定)
#' @param saa_bio_year 選択率をSPRに換算するときに生物パラメータを取り出す年の範囲(data_futureがある場合将来予測年も指定可能。生物パラメータが密度によって変わる場合)
#' @param saa_bio 選択率をSPRに換算するとき用いる生物パラメータそのものtibble(waa=waa, maa=maa, M=M)として指定する
#'
#' @export
#'
#'

convert_Fvector <- function(res_vpa=NULL,
                            res_future = NULL,
                            faa_vector=NULL,
                            faa_vector_year=NULL,
                            faa_bio_year=NULL,
                            faa_bio=NULL,
                            saa_vector=NULL,
                            saa_vector_year=NULL,
                            saa_bio_year=NULL,
                            saa_bio=NULL){

  assert_that(
    (is.null(faa_vector)   | is.null(faa_vector_year)),
    (is.null(saa_vector)   | is.null(saa_vector_year)),
    (is.null(saa_bio_year) | is.null(saa_bio)),
    (is.null(faa_bio_year) | is.null(faa_bio))
  )

  if(is.null(faa_vector)) faa_vector <- apply_year_colum(res_vpa$faa,target_year=faa_vector_year)
  if(is.null(saa_vector)) saa_vector <- apply_year_colum(res_vpa$faa,target_year=saa_vector_year)

  # faa_vectorが何%のSPRにあたるか
  faa_perSPR <- calc_future_perSPR(fout    = res_future,
                                   res_vpa = res_vpa,
                                   Fvector = faa_vector,
                                   target.year = faa_bio_year,
                                   biopar = faa_bio)
  cat("%SPR in faa=", faa_perSPR,"\n")
  # saa_vectorがfaa_vectorに相当する%SPRになるためには何倍にしないといけないか
  saa_multiplier <- calc_future_perSPR(fout=res_future,
                                       res_vpa=res_vpa,
                                       Fvector=saa_vector,
                                       target.year=saa_bio_year,
                                       SPRtarget=faa_perSPR,
                                       biopar=saa_bio)
  # faaの漁獲圧の大きさに相当するsaaの選択率を持ったF at age
  Fvector <- saa_vector/saa_multiplier$Fratio
  return(lst(Fvector, faa_perSPR))
}


# 縦の行列の年によって足し算する; 関数の汎用版
#' @export

rowtapply2 <- function(a0,FUN.name){
    FUN <- get(FUN.name)
    yname <- floor(as.numeric(rownames(a0)))
    res <- matrix(0,length(unique(yname)),ncol(a0))
    for(i in 1:ncol(a0)){
        res[,i] <- tapply(a0[,i],yname,FUN)
    }
    dimnames(res) <- list(unique(yname),colnames(a0))
    res
}

#'
#' kobe.tableをさらにsummaryする
#'
#' @param target_threshold c(60, 50)みたいな2つの長さのベクトル。一番目はbeta=0.8のときのtargetを上回る確率、2番めは50%のときの。資源状態が良い場合には1番目の値は2番めの値よりも大きいが、資源状態が悪いと1番目の値は2番めよりも小さくなる。その場合には自動的にc(100,50)となるように置き換わる(つまりランク3は出現しない)
#' @param risk_threshold c(0.2,15) みたいな2つの長さのベクトル。一番目はbeta=0.8のときに10年間でずっとthresholdを上回る確率、2番めは50%のとき。資源状態が良い場合には1番目の値は2番めの値よりも小さくなる。
#' @param ssbpercent_summary_year 目標管理基準値を上回るかどうかを判断する年
#' @param ssb_summary_year パフォーマンス指標として取り出すSSBの年
#' @param catch_summary_year パフォーマンス指標として取り出すCatchの年
#' @param SBmsy ここで与えた相対値としてSBのパフォーマンス指標を示す
#' @param MSY ここで与えた相対値としてMSYのパフォーマンス指標を示す
#' @param Bthreshold_label ランクづけに利用するリスクをkobe.tableから持ってくるときのkobe.tableのリストの名前
#'
#' @export

summary_kobe_table <- function(kobeII.table,
                               target_threshold,
                               risk_threshold,
                               ssbpercent_summary_year=c(2031),
                               ssb_summary_year=c(2025,2031),
                               catch_summary_year=c(2021,2026,2031),
                               SBmsy=1, MSY=1,
                               Bthreshold_label="blimit.risk",
                               sort_result_table=FALSE){
  tmpfunc_ <- function(x,header)  str_c(header,"_",x)

  # riskのしきい値をどうするか?の設定が必要
  ssb_table <-  kobeII.table$ssb.mean %>% select(as.character(ssb_summary_year))/SBmsy
  ssb_table <-  ssb_table  %>% as_tibble() %>% rename_all(tmpfunc_, header="SSB")
  catch_table <- kobeII.table$catch.mean %>% select(as.character(catch_summary_year))/MSY
  catch_table <- catch_table %>% as_tibble() %>% rename_all(tmpfunc_, header="Catch")

  summary_HCR <-
      bind_cols(
          kobeII.table$ssb.mean %>% select("HCR_name","beta"),
          kobeII.table$prob.over.ssbtarget %>% select(as.character(ssbpercent_summary_year)) %>% rename_all(tmpfunc_, header="SSB_prob"),
         ssb_table,
         catch_table,
          kobeII.table[Bthreshold_label][[1]] %>% ungroup %>% select(value) %>% dplyr::rename(risk_Bthreshold=value),
          kobeII.table$bban.risk %>% ungroup %>% select(value) %>% dplyr::rename(risk_Bban=value),
          kobeII.table$catch.risk %>% ungroup %>% select(value) %>% dplyr::rename(risk_catch=value)
      )

  if(target_threshold[1]<target_threshold[2]){ target_threshold[1] <- 100; risk_threshold[1] <- 0}
  summary_HCR$SSB_prob_XXXX <- summary_HCR[str_c("SSB_prob_", ssbpercent_summary_year)]
  summary_HCR$Catch_XXXX <- summary_HCR[str_c("Catch_", catch_summary_year[1])]
  summary_HCR <- summary_HCR %>%
      mutate(category_target = case_when((SSB_prob_XXXX >= target_threshold[1]) ~ 2,
                                         (SSB_prob_XXXX <  target_threshold[1] & SSB_prob_XXXX >= (target_threshold[2]-0.5)) ~ 1,
                                         (SSB_prob_XXXX < (target_threshold[2]-0.5)) ~ 0),
             category_risk   = case_when((risk_Bthreshold <= risk_threshold[1]) ~ 2,
                                         (risk_Bthreshold >  risk_threshold[1] & risk_Bthreshold <= risk_threshold[2]) ~ 1,
                                         (risk_Bthreshold >  risk_threshold[2]) ~ 0)) %>%
      mutate(category_risk = ifelse(sum(risk_threshold)==0, 2, category_risk)) %>%
      mutate(category_combine = str_c(category_target,category_risk)) %>%
      mutate(category         = case_when(category_combine=="10" ~ 1,
                                          category_combine=="20" ~ 1.5,
                                          category_combine=="11" ~ 2,
                                          category_combine=="12" ~ 2.5,
                                          category_combine=="21" ~ 2.5,
                                          category_combine=="22" ~ 3,
                                          category_target==0     ~ 0)) %>%
      select(category, HCR_name:risk_catch, category_risk, category_target)


  if(sort_result_table==TRUE){
    summary_HCR <- summary_HCR %>%
        arrange(desc(category), desc(Catch_XXXX))
  }
  return(summary_HCR)
}

#'
#' VPAデータが1年分追加されたダミーデータを生成する
#'
#' @export
#'

create_dummy_vpa <- function(res_vpa){

  res_vpa_updated <- res_vpa

  add_1year <- function(naa){
    nyear <- ncol(naa)
    year_name <- colnames(naa) %>% as.numeric()
    year_name <- c(year_name,max(year_name)+1)
    nage  <- nrow(naa)
    empty_matrix <- matrix(0, nage, nyear+1)
    dimnames(empty_matrix) <- list(rownames(naa), year_name)

    empty_matrix[,-nyear] <- as.matrix(naa)
    empty_matrix[, nyear] <- naa[,nyear]
    as.data.frame(empty_matrix)
  }

  res_vpa_updated$naa           <- add_1year(res_vpa$naa)
  res_vpa_updated$faa           <- add_1year(res_vpa$faa)
  res_vpa_updated$input$dat$waa <- add_1year(res_vpa$input$dat$waa)
  res_vpa_updated$input$dat$maa <- add_1year(res_vpa$input$dat$maa)
  res_vpa_updated$input$dat$M   <- add_1year(res_vpa$input$dat$M  )
  res_vpa_updated$input$dat$caa <- add_1year(res_vpa$input$dat$caa)

  if(!is.null(res_vpa_updated$input$dat$release.all))  res_vpa_updated$input$dat$release.all <- res_vpa_updated$input$dat$release.all %>% add_1year()
  if(!is.null(res_vpa_updated$input$dat$release.alive))  res_vpa_updated$input$dat$release.alive <- res_vpa_updated$input$dat$release.alive %>% add_1year()
  if(!is.null(res_vpa_updated$input$dat$release.ratealive))  res_vpa_updated$input$dat$release.ratealive <- res_vpa_updated$input$dat$release.ratealive %>% add_1year()
  if(!is.null(res_vpa_updated$input$dat$release.dat))  res_vpa_updated$input$dat$release.all <- res_vpa_updated$input$dat$release.dat %>% add_1year()

  if(!is.null(res_vpa_updated$input$dat$waa.catch))
    res_vpa_updated$input$dat$waa.catch <- add_1year(res_vpa$input$dat$waa.catch)

  res_vpa_updated$baa <- res_vpa_updated$input$dat$waa * res_vpa_updated$naa
  res_vpa_updated$ssb <- res_vpa_updated$input$dat$maa * res_vpa_updated$baa
  res_vpa_updated$wcca <- res_vpa_updated$input$dat$caa * res_vpa_updated$input$dat$waa

  return(res_vpa_updated)
}


#'
#' res_futureから考えられるほぼすべてのパフォーマンス指標をとりだし、2行のtibbleで返す
#'
#' @param res_future future_vpaの返り値
#' @param SBtarget それを上回る・下回る確率を計算する。目標管理基準値
#' @param SBlimit それを上回る・下回る確率を計算する。限界管理基準値
#' @param SBban それを上回る・下回る確率を計算する。禁漁水準
#' @param SBmin それを上回る・下回る確率を計算する。過去最低親魚量
#' @param MSY それを上回る・下回る確率を計算する。MSY
#' @param is_scale TRUEの場合、親魚量はSBtargetで、漁獲量はMSYで割った相対値が出力される。それ以外はそのまま
#' @param unit 確率を出力するときの単位。100を入れると%単位で結果が返される
#' @param period_extra デフォルトではSSBminなどを一度でも下回るなど、期間を指定して計算する統計量はABC_year + 0:9, 0:4, 5:9, 1:10, 1:4, 6:10で決め打ちしているが、それ以外の期間を指定したいときにここの引数で与える
#' @param type "AS": age-structured from frasyr, "PM": production model from frapmr
#'
#'
# #' @examples
# #' \dontrun{
# #'   data(res_future_HSL2)
# #'   calculate_all_pm(res_future_HSL2, 0, 0, 0, 0)
# #' }
#'
#' @export
#'

calculate_all_pm <- function(res_future, SBtarget=-1, SBlimit=-1, SBban=-1, SBmin=-1, MSY=-1, is_scale=FALSE, unit=1, period_extra=NULL, type="AS", fun_period=mean){
    # by year performance (start_future_year:last_year)
    # mean, median, ci5%, ci10%, ci90%, ci95%, CV,
    # ssb, biomass, number by age, catch weight
    # probability > SBtarget, SBmin, SBlimit, SBban, SBmax

    get_annual_pm <- function(mat,fun,label){
        x <- apply(mat,1,fun)
        tibble(stat=str_c(label,"_",names(x)),value=x)
    }

    fun_list <- list(ci0.05=function(x) quantile(x,0.05, na.rm=TRUE),
                     ci0.10 = function(x) quantile(x,0.1, na.rm=TRUE),
                     ci0.95 = function(x) quantile(x,0.95, na.rm=TRUE),
                     ci0.90 = function(x) quantile(x,0.90, na.rm=TRUE),
                     cv = function(x) sqrt(var(x, na.rm=TRUE))/mean(x, na.rm=TRUE),
                     mean = function(x) mean(x, na.rm=TRUE),
                     median = function(x) median(x, na.rm=TRUE),
                     prob_target = function(x) mean(x>SBtarget, na.rm=TRUE)*unit,
                     prob_limit  = function(x) mean(x>SBlimit, na.rm=TRUE)*unit,
                     prob_ban    = function(x) mean(x>SBban, na.rm=TRUE)*unit,
                     prob_min    = function(x) mean(x>SBmin, na.rm=TRUE)*unit)

  if(is_scale){
    scale_ssb <- SBtarget
    scale_catch <- MSY
    SBlimit     <- SBlimit/SBtarget
    SBban       <- SBban/SBtarget
    SBmin       <- SBmin/SBtarget
    SBtarget    <- SBtarget/SBtarget
  }
  else{
    scale_ssb <- 1
    scale_catch <- 1
  }

  if(type=="AS"){
    year_future <- dimnames(res_future$naa)[[2]][res_future$input$tmb_data$future_initial_year:dim(res_future$naa)[[2]]]
    age_label <- str_c("A",dimnames(res_future$naa)[[1]])
    year_label <- dimnames(res_future$naa)[[2]]

    ssb_mat <- res_future$SR_mat[year_future,,"ssb"] / scale_ssb
    catch_mat <- res_future$HCR_realized[year_future,,"wcatch"] / scale_catch
    biom_mat <- apply(res_future$naa * res_future$waa,c(2,3),sum)
  }
  if(type=="PM"){
    year_future <- res_future$mat_year$year[!res_future$mat_year$is_est]
    age_label   <- NULL
    year_label  <- res_future$mat_year$year

    if(!is_scale){
      ssb_mat   <- res_future$mat_stat["B",,] 
      catch_mat <- res_future$mat_stat["C",,] 
      biom_mat  <- res_future$mat_stat["B",,] 
    }
    else{
      ssb_mat   <- res_future$mat_stat["Bratio",,] 
      catch_mat <- res_future$mat_stat["Cratio",,] 
      biom_mat  <- res_future$mat_stat["Bratio",,]       
    }
  }

  stat_data <- NULL
  for(i in 1:length(fun_list)){
    fun <- fun_list[[i]]
    funname <- names(fun_list)[i]

    if(type=="AS"){
      x1 <- purrr::map_dfr(seq_len(length(age_label)),
                       function(x) get_annual_pm(res_future$naa [x,year_future,],
                                                 fun,str_c(funname,"_naa_", age_label[x])))
      x2 <- purrr::map_dfr(seq_len(length(age_label)),
                     function(x) get_annual_pm(res_future$wcaa[x,year_future,],
                                               fun,str_c(funname,"_wcaa_",age_label[x])))
      x3 <- purrr::map_dfr(seq_len(length(age_label)),
                       function(x) get_annual_pm(res_future$faa [x,year_future,],
                                                 fun,str_c(funname,"_faa_", age_label[x])))      
    }
    else{
      x1 <- x2 <- x3 <- NULL
    }
    tmp <- bind_rows(
      x1,x2,x3,
      get_annual_pm(ssb_mat,  fun  ,str_c(funname,"_ssb")),
      get_annual_pm(catch_mat,fun  ,str_c(funname,"_catch")),
      get_annual_pm(biom_mat, fun  ,str_c(funname,"_biom"))
    )
    stat_data <- bind_rows(stat_data, tmp)
  }

  # temporal scale
  # by term performance
  # - management term1, ABC_year + 0:9
  # - management term2, ABC_year + 0:4
  # - management term3, ABC_year + 5:9
  # - management term4, ABC_year + 1:10
  # - management term5, ABC_year + 1:4
  # - management term6, ABC_year + 6:10
  #
  # mean(ssb), mean(biomass), mean(catch), mean(AAV), mean(CV), Pr(SBany>SBtarget), Pr(SBany>SBlimit),Pr(SBany>SBban),Pr(SBany>SBmin), AAV_min (査読コメントを参照), max_AAV_min, mean(min_catch)

  get_period_pm <- function(mat, fun_name, year_period, sum_fun_name, fun_name_char){
    tmp <- rownames(mat)%in%year_period
    if(sum(tmp)!=length(year_period)){
        cat("part of years does not match future projection year\n")
        return(NA)
    }
    else{
      if(type=="PM"){
        if(fun_name_char==c("prob_limit_any")){
          mat <- sweep(mat, 2, SBlimit, FUN="/")
        }
        if(fun_name_char==c("prob_ban_any")){
          mat <- sweep(mat, 2, SBban, FUN="/")
        }
        if(fun_name_char==c("prob_min_any")){
          mat <- sweep(mat, 2, SBmin, FUN="/")
        }                
      }
      mat1 <- mat[tmp,]
      res <- apply(mat1,2,fun_name) %>% sum_fun_name()
      return(res)
    }
  }

  if(type=="AS"){
    ABC_year     <- year_label[res_future$input$tmb_data$start_ABC_year] %>% as.numeric()
  }
  if(type=="PM"){
    ABC_year     <- res_future$mat_year$year[res_future$mat_year$is_manage] %>% min()
  }
  period_range <- list(0:9, 0:4, 5:9, 1:10, 1:4, 6:10)
  if(!is.null(period_extra)){
    for(i in 1:length(period_extra)){
      if(sum(purrr::map_lgl(period_range, function(x) all(x==period_extra[[i]])))==0){
        period_range <- c(period_range, period_extra[i])
      }
    }
  }
  period_list  <- purrr::map(period_range, function(x) ABC_year + x)
  names(period_list) <- purrr::map_chr(period_list, function(x) str_c(range(x),collapse="."))

  av <- function(x){
    av_value <- (x[-1]-x[-length(x)])/x[-length(x)]
    av_value <- av_value[!is.nan(av_value) & av_value<Inf]
    if(length(av_value)==0) return(NA) else av_value
  }

  calc_mdr_ <- function(x){
    if(!is.na(av(x)[1])){
#      if(min(av(x))>0) cat(min(av(x)),"\n")
      return(min(c(av(x),0), na.rm=TRUE))
    }
    else{
      return(NA)
    }
  }

  calc_mdr0_ <- function(x){
    if(!is.na(av(x)[1])){
      return(min(av(x), na.rm=TRUE))
    }
    else{
      return(NA)
    }
  }

  mean2 <- function(x) fun_period(x,na.rm=TRUE)
    
  fun_list2 <- list(cv     = function(x) sd(x, na.rm=TRUE)/mean(x,na.rm=TRUE),
                    mean   = function(x) mean(x,na.rm=TRUE),
                    median   = function(x) median(x,na.rm=TRUE),
                    aav   = function(x) mean(abs(av(x)),na.rm=TRUE),
                    mav   = function(x) median(abs(av(x)),na.rm=TRUE),                    
                    adr = function(x){ x0 <- av(x) ;
                                       x0[x0>0] <- NA ;
                                       mean(x0,na.rm=TRUE)    },
                    mdr = calc_mdr_,
                    mdr0 = calc_mdr0_,
                    min_value = function(x){
                        if(all(is.na(x))) NA else min(x, na.rm=TRUE)
                    },
                    max_value = function(x){
                        if(all(is.na(x))) NA else max(x, na.rm=TRUE)
                    },
                    # 以下、na.rm=FALSEに変更。PMで全てのデータがNAであった場合、na.rm=TRUEとすると0という数字が入ってしまうので、PMで使うならここはFALSEとすべき
                    # ただし、一部だけNAが入るような場合にはどうすべきか?そういう事例がでたときに要検討
                    prob_target_any  = function(x) ifelse(sum(x<SBtarget,na.rm=FALSE)>0,1,0),
                    prob_limit_any  = function(x){
                      if(type=="PM") SBlimit <- rep(1,length(x))
                      ifelse(sum(x<SBlimit,na.rm=FALSE)>0,1,0)
                    },
                    prob_half_any  = function(x) ifelse(sum(x[-1]<0.5*x[-length(x)])>0,1,0),
                    prob_ban_any    = function(x){
                      if(type=="PM") SBban <- rep(1,length(x))                      
                      ifelse(sum(x<SBban,na.rm=FALSE)>0,1,0)
                    },
                    prob_min_any    = function(x){
                      if(type=="PM") SBmin <- rep(1,length(x))                      
                      ifelse(sum(x<SBmin,na.rm=FALSE)>0,1,0)
                    })

  mat_list <- lst(ssb=ssb_mat, biom=biom_mat, catch=catch_mat)
  for(j in seq_len(length(fun_list2))){
    for(i in seq_len(length(period_list))){
      stat_data <- bind_rows(
        stat_data,
        purrr::map_dfr(1:length(mat_list),
                       function(x)
                           tibble(stat=str_c(names(fun_list2)[j],names(mat_list)[x],names(period_list)[i],sep="_"),
                                value=get_period_pm(mat_list[[x]], fun_list2[[j]], period_list[[i]], mean2, fun_name_char=names(fun_list2)[j])))
        )
    }}
  return(stat_data)
}

#'
#' @export
#'


derive_all_stat <- function(res, word_vector, exact=FALSE){
  if(exact==FALSE){
    for(i in 1:length(word_vector)){
      res <- res %>% dplyr::filter(str_detect(stat,word_vector[i]))
    }
    if(nrow(res)==0) stop("no match") else return(res)
  }
  else{
    res <- res %>% dplyr::filter(stat==word_vector)
    return(res)
  }
}


#' リスクテーブルを作る
#'
#' @param target_threshold カテゴリ分けに使うSSB>SSBtargetの確率。1番目がbeta=0.8のときの値、2番目が50%のときの値(なのでだいたい50%になるはず)
#' @param limit_threshold  カテゴリ分けに使うSSBall>SSBthresholdの確率。1番目がbeta=0.8のときの値、2番目が50%のときの値(なので1番目よりも2番めのほうが小さいはず)
#'
#' @export
#'

rank_HCR <- function(summary_HCR,
                     target_threshold,
                     risk_threshold,
                     target_prob_name="SBtar_prob",
                     risk_prob_name="SBlim_prob"){

  SSB_prob_vector <- summary_HCR[target_prob_name] %>% unlist %>% as.numeric
  risk_prob_vector <- summary_HCR[risk_prob_name]  %>% unlist %>% as.numeric
  target_threshold <- target_threshold
  risk_threshold   <- risk_threshold

  if(target_threshold[1]<target_threshold[2]) target_threshold[1] <- Inf

  summary_HCR <- summary_HCR %>%
      mutate(category_target = case_when((SSB_prob_vector >= target_threshold[1]) ~ 2,
                                         (SSB_prob_vector <  target_threshold[1] & SSB_prob_vector >= (target_threshold[2]*0.99)) ~ 1,
                                         (SSB_prob_vector < (target_threshold[2]*0.99)) ~ 0),
             category_risk   = case_when((risk_prob_vector <= risk_threshold[1]) ~ 2,
                                         (risk_prob_vector >  risk_threshold[1] & risk_prob_vector <= risk_threshold[2]) ~ 1,
                                         (risk_prob_vector >  risk_threshold[2]) ~ 0))

  if(sum(risk_threshold)==0) summary_HCR$category_risk[] <- 2

  summary_HCR <- summary_HCR %>%
      mutate(category_combine = str_c(category_target,category_risk)) %>%
      mutate(category         = case_when(category_combine=="10" ~ 1,
                                          category_combine=="20" ~ 1.5,
                                          category_combine=="11" ~ 2,
                                          category_combine=="12" ~ 2.5,
                                          category_combine=="21" ~ 2.5,
                                          category_combine=="22" ~ 3,
                                          category_target==0     ~ 0))
  return(summary_HCR)

}

#'
#' beverton-holtのh,R0とbioparsを与えるとa,bを返す関数
#'
#' @export

get.ab.bh <- function(h,R0,biopars){

    get.SPR0 <- function(M,maa,waa,output="simple"){
        nage <- length(M)
        S <- exp(-M)
        N <- numeric()
        N[1] <- 1
        for(i in 2:(nage-1)) N[i] <- N[i-1]*S[i-1]
        N[nage] <- N[nage-1] * S[nage]/(1-S[nage])
        SPR0 <- sum(N * maa * waa)
        if(output=="simple") return(SPR0) else return(listN2(N,SPR0))
    }
    SPR0 <- get.SPR0(biopars$M,biopars$maa,biopars$waa)
    S0 <- R0*SPR0
    beta <- (5*h-1)/(4*h*R0)
    alpha <- SPR0*(1-h)/(4*h)
    a <- 1/alpha
    b <- beta/alpha
    return(tibble::lst(SPR0,R0,h,S0,a,b))
}

#' 漁獲量の上限設定をしたときの設定がちゃんと生きているかどうかを確かめる
#'
#' @export
#'

check_fix_CVoption <- function(res_future){
  wcatch <- res_future$HCR_realized[,,"wcatch"]
  wcatch[-1,]/wcatch[-nrow(wcatch),]
}
        
#' @export      

format_type <- function(){
    tribble(~name, ~col, ~ lty,
            "target", "#00533E", "dashed",
            "limit",  "#EDB918", "dotdash",
            "ban",   "#C73C2E", "dotted")
}
ichimomo/frasyr documentation built on May 3, 2024, 1:30 a.m.