---
title: "04 Create IC format from estimation results"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{04 Create IC format from estimation results}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

## Introduction
This document shows how to take the output from doEstimationForAllStrata() and format it into the InterCatch format. This is done by using the function exportEstimationResultsToInterCatchFormat() in the RDBEScore package.

## Load the package

```{r setup}
library(RDBEScore)
```

## Step 1) load and prepare some test data 

```{r}
# load an example dataset

# H1 directory
dirH1 <- "../tests/testthat/h1_v_20250211/"

fixProbs <- function(myH1RawObject){

  # Edit our data so that we have SRSWOR on each level and calculate the probs
  myH1RawObject[["VS"]]$VSselectMeth <- "SRSWOR"
  myH1RawObject[["VS"]]$VSincProb <- myH1RawObject[["VS"]]$VSnumSamp / myH1RawObject[["VS"]]$VSnumTotal
  myH1RawObject[["VS"]]$VSselProb <- 1/myH1RawObject[["VS"]]$VSnumTotal
  myH1RawObject[["FT"]]$FTselectMeth <- "SRSWOR"
  myH1RawObject[["FT"]]$FTincProb <- myH1RawObject[["FT"]]$FTnumSamp / myH1RawObject[["FT"]]$FTnumTotal
  myH1RawObject[["FT"]]$FTselProb <- 1/myH1RawObject[["FT"]]$FTnumTotal
  myH1RawObject[["FO"]]$FOselectMeth <- "SRSWOR"
  myH1RawObject[["FO"]]$FOincProb <- myH1RawObject[["FO"]]$FOnumSamp / myH1RawObject[["FO"]]$FOnumTotal
  myH1RawObject[["FO"]]$FOselProb <- 1/myH1RawObject[["FO"]]$FOnumTotal
  myH1RawObject[["SS"]]$SSselectMeth <- "SRSWOR"
  myH1RawObject[["SS"]]$SSincProb <- myH1RawObject[["SS"]]$SSnumSamp / myH1RawObject[["SS"]]$SSnumTotal
  myH1RawObject[["SS"]]$SSselProb <- 1/myH1RawObject[["SS"]]$SSnumTotal
  myH1RawObject[["SA"]]$SAselectMeth <- "SRSWOR"
  myH1RawObject[["SA"]]$SAincProb <- myH1RawObject[["SA"]]$SAnumSamp / myH1RawObject[["SA"]]$SAnumTotal
  myH1RawObject[["SA"]]$SAselProb <- 1/myH1RawObject[["SA"]]$SAnumTotal

  # Update our test data with some random sample measurements (it didn't include these)
  # set the random seed
  set.seed(1234)
  myH1RawObject[['SA']]$SAsampWtLive <- round(runif(n=nrow(myH1RawObject[['SA']]),min = 1, max = 100))
  myH1RawObject[['SA']]$SAsampWtMes <- round(runif(n=nrow(myH1RawObject[['SA']]),min = 1, max = 100))

  myH1RawObject

}


## Step 1) load and prepare some test data
myH1RawObject <- createRDBESDataObject(input  = dirH1)
#Filter our data for WGRDBES-EST TEST 1, 1965, H1
myCatchFraction <- "Lan"
mySpecies <- 1019159

myFields <- c("DEyear","DEhierarchy","DEsampScheme","DEstratumName","SAspeCode","SScatchFra")
myValues <- c(1965,1,"National Routine","DE_stratum1_H1",mySpecies,myCatchFraction)
myH1RawObject <- filterRDBESDataObject(myH1RawObject,
                                       fieldsToFilter = myFields,
                                       valuesToFilter = myValues )
myH1RawObject <- findAndKillOrphans(myH1RawObject)
myH1 <- fixProbs(myH1RawObject)

```

## Step 2) Estimate for live weight for the data


```{r}
myTestData <- createRDBESEstObject(myH1, 1, stopTable = "SA")
# Get rid of rows that don't have an SA row
myTestData <- myTestData[!is.na(myTestData$SAid),]
# Estimate using the data
myStrataResults <- doEstimationForAllStrata(myTestData, "SAsampWtLive")

```

## Step 3) Export the results to InterCatch format (just HI and SI for this data)

```{r}

# Get our estimated values for the PSU
psuEstimates <- myStrataResults[myStrataResults$recType == "VS",]

# This is the data we will export to IC format
dataToOutput <- data.frame(Country = myH1[["SD"]]$SDctry,
                           Year = myH1[["DE"]]$DEyear,
                           SeasonType = NA,
                           Season = NA,
                           Fleet= NA,
                           AreaType = "Stratum",
                           FishingArea = psuEstimates$stratumName,
                           DepthRange = "NA",
                           Species = mySpecies ,
                           Stock = mySpecies,
                           CatchCategory = substr(myCatchFraction, 1, 1),
                           ReportingCategory = "A",
                           Usage = "H",
                           SamplesOrigin = "O",
                           UnitCATON = "kg",
                           CATON = psuEstimates$est.total/1000.0,
                           varCATON = psuEstimates$var.total/1000000.0,
                           Sex = "NA",
                           PlusGroup = "NA",
                           MeanWeight = -9,
                           unitMeanWeight = "kg",
                           unitCANUM = "n",
                           UnitAgeOrLength = "year",
                           NumberCaught = -9,
                           MeanWeight = -9,
                           CANUMtype = "Age",
                           AgeLength = c(1,1,1,2,2,2,3,3,3),
                           NumSamplesAge = c(300,300,300,200,200,200,100,100,100))

tempIC <- exportEstimationResultsToInterCatchFormat(dataToOutput)
print(tempIC)
```

See also other package vignettes:

```{r,echo=F, results='asis', eval=T}
# note: vignettes are first compiled and only then html, R and Rmd are moved to folder doc
vignettes <- data.frame(Rmd=dir()[grepl(dir(), pat=".Rmd")])
# determines index of vignette being compiled
thisVignetteRmd<-paste0(gsub(":","",gsub(" ", "-",as.character(rmarkdown::metadata$title))),".Rmd")
# prints list with links
for (i in seq_along(vignettes$Rmd)) {
  if(tolower(vignettes$Rmd[i])!=tolower(thisVignetteRmd)){
  cat(sprintf("- [%s](%s)\n", gsub("-"," ",gsub(".Rmd","",vignettes$Rmd[i])), gsub(".Rmd",".html",vignettes$Rmd[i])))
}
}
```

#END
