library(stockassessment)

oldwd<-setwd("data")

catch.no<-read.ices('cn.dat')
catch.mean.weight<-read.ices('cw.dat')
dis.mean.weight<-read.ices('dw.dat')
land.mean.weight<-read.ices('lw.dat')
stock.mean.weight<-read.ices('sw.dat')
prop.mature<-read.ices('mo.dat')
natural.mortality<-read.ices('nm.dat')
surveys<-read.ices('survey.dat')
land.no<-read.ices('lf.dat')
dis.no<-catch.no-land.no 
prop.f<-read.ices('pf.dat')
prop.m<-read.ices('pm.dat')


# redirect output to avoid mgcv message generating warning icons
zz=file("dummy",open="wt")
sink(zz)
sink(zz,type="message")
library(mgcv)
sink(zz,type="message")
sink()
file.remove("dummy")


sdQ1 = read.table("q1sds.txt")
sdQ3 = read.table("q3sds-noage0.txt")
sdQ30 = read.table("q3sds-age0.txt")

attributes(surveys[[1]])$weight = (1/sdQ1^2)[,1:5]
attributes(surveys[[2]])$weight = (1/sdQ3^2)[,1:4]
attributes(surveys[[3]])$weight = (1/sdQ30^2)[,1]



setwd(oldwd)

# Modify to 6+ data
cutage<-6
low<-1
GE<-which(as.numeric(colnames(catch.no))>=cutage)
E<-which(as.numeric(colnames(catch.no))==cutage)
w<-catch.no[,GE]/rowSums(catch.no[,GE])
wex<-rbind(w,w[nrow(w),])
wD<-dis.no[,GE]/ifelse(rowSums(dis.no[,GE])>0,rowSums(dis.no[,GE]),1)
wL<-land.no[,GE]/rowSums(land.no[,GE])
catch.no[,E]<-rowSums(catch.no[,GE])
catch.no<-catch.no[,low:E]
prop.mature[,E]<-rowSums(prop.mature[,GE]*wex)
prop.mature<-prop.mature[,low:E]
stock.mean.weight[,E]<-rowSums(stock.mean.weight[,GE]*wex)
stock.mean.weight<-stock.mean.weight[,low:E]
catch.mean.weight[,E]<-rowSums(catch.mean.weight[,GE]*w)
catch.mean.weight<-catch.mean.weight[,low:E]
dis.mean.weight[,E]<-rowSums(dis.mean.weight[,GE]*wD)
dis.mean.weight<-dis.mean.weight[,low:E]
land.mean.weight[,E]<-rowSums(land.mean.weight[,GE]*wL)
land.mean.weight<-land.mean.weight[,low:E]
natural.mortality[,E]<-rowSums(natural.mortality[,GE]*wex)
natural.mortality<-natural.mortality[,low:E]
land.no[,E]<-rowSums(land.no[,GE])
land.no<-land.no[,low:E]
land.frac<-ifelse(catch.no>0,land.no/catch.no,1)
prop.f<-prop.f[,low:E]
prop.m<-prop.m[,low:E]




Ms <- natural.mortality
M <- Ms
M[(nrow(Ms)-9):nrow(Ms), 3:6] <- Ms[(nrow(Ms)-9):nrow(Ms), 3:6] + -log(1-0.15)




dat<-setup.sam.data(surveys=surveys,
                    residual.fleet=catch.no, 
                    prop.mature=prop.mature, 
                    stock.mean.weight=stock.mean.weight, 
                    catch.mean.weight=catch.mean.weight, 
                    dis.mean.weight=dis.mean.weight, 
                    land.mean.weight=land.mean.weight,
                    prop.f=prop.f, 
                    prop.m=prop.m, 
                    natural.mortality=M, 
                    land.frac=land.frac)
                    
save(dat, file="run/data.RData")