#### ***********************************************************************
#### These routines are for the Generalised Visual Fast Count Technique
#### for video data. The functions are available in the package emon -
#### which can be downloaded from
#### https://r-forge.r-project.org/projects/emon/
#### ***********************************************************************

library(emon)
ls(name="package:emon")

#************************************************
# Code to create Figure 1
#************************************************
#### Set appropriate working directory

setwd("C:\\Documents and Settings\\jtb00\\My Documents\\Fish Stuff\\Consultancy 2014\\Generalised VFC\\Data analysis")
s=5
x = rep(0, s)
m.vec = seq(0.5,20,0.5); lm.vec = length(m.vec)

#### nsims = 100,000 in paper, but takes a while to compute
nsims = 1000
est.vfc.2pos = rep(0, nsims); est.vfc.1pos = rep(0, nsims); est.vfc.3pos = rep(0, nsims)
est.mom.1pos = rep(0, nsims); est.mom.2pos = rep(0, nsims); est.mom.3pos = rep(0, nsims)
mean.mom.1pos = rep(0, lm.vec); mean.mom.2pos = rep(0, lm.vec); mean.mom.3pos = rep(0, lm.vec)
theory.1pos = rep(0, lm.vec); theory.2pos = rep(0, lm.vec); theory.3pos = rep(0, lm.vec)
#
for (k in 1:lm.vec) {
m = m.vec[k]

# Calculate mean of the MOM estimator
for (j in 1:nsims) {
x = rpois(5, m/s)
est.mom.1pos[j] = GVFCMOM(x, s, d=1, method='pois')
est.mom.2pos[j] = GVFCMOM(x, s, d=2, method='pois')
est.mom.3pos[j] = GVFCMOM(x, s, d=3, method='pois')
}
mean.mom.1pos[k] = mean(est.mom.1pos)
mean.mom.2pos[k] = mean(est.mom.2pos)
mean.mom.3pos[k] = mean(est.mom.3pos)
}

# Calculate the bias of the VFC estimator
for (k in 1:lm.vec) {
m = m.vec[k]
theory.3pos[k] = 100*(expected.pois(m, 5, 3) - m) / m
theory.2pos[k] = 100*(expected.pois(m, 5, 2) - m) / m
theory.1pos[k] = 100*(expected.pois(m, 5, 1) - m) / m
}

bias.mom.1pos = 100*(mean.mom.1pos - m.vec) / m.vec
bias.mom.2pos = 100*(mean.mom.2pos - m.vec) / m.vec
bias.mom.3pos = 100*(mean.mom.3pos - m.vec) / m.vec

#tiff("Figure1.tif", height=12, width=12, res=1000, pointsize=8, units="cm", compression='lzw')
plot(m.vec, theory.1pos, ylab='% bias', xlab=expression(paste('Mean per transect ', lambda)), type='l', ylim=c(-10,110))
lines(m.vec, bias.mom.1pos, col=gray(0.6))
lines(m.vec, theory.2pos, lty=2)
lines(m.vec, bias.mom.2pos, lty=2, col=gray(0.6))
lines(m.vec, theory.3pos, lty=3)
lines(m.vec, bias.mom.3pos, lty=3, col=gray(0.6))
legend('topright',legend=c("GVFC (d=1)","GVFCMOM (d=1)","GVFC (d=2)","GVFCMOM (d=2)","GVFC (d=3)","GVFCMOM (d=3)"),lty=c(1,1,2,2,3,3), col=c("black",gray(0.6),"black",gray(0.6),"black",gray(0.6)))
#dev.off()

#*************************************************
# Data analysis for the Folkestone Pomerania data
#*************************************************
S08 = c(1,1)
C11 = c(1,0,0,0,0,0,0,0,0)
S01 = c(0,0,0,0,0,0,0,0,0,0)
R06 = c(0,0,0,0,1,1)
C17 = c(1,0,0,0,0,0,0,0,0,0)
R04 = c(0,0,1,0,2)
C02 = c(1,0,0,1)
C09 = c(0,0,0,0,0,0,3,0,1)
C21 = c(1,1)
C26 = c(0,1,2)  
C05 = c(1,3)  
C07 = c(0,1,1)
C14 = c(2,0,2)
R02 = c(1,0,8)
R07 = c(0,0,0,1,3)
R08 = c(0,0,0,5,1)
S01 = c(0,0,0,0,0,0,0,0,0,0)
S06 = c(0,0,0,0,0,0,0,0,0,0)
S13 = c(0,0,0,0,0,0,0,0,0,0)
S15 = c(0,0,0,1,1)
S17 = c(0,0,1,1)

lat = c(51.00624439, 51.02346979, 50.99115206, 51.0173252, 51.02944899,
   51.02654, 51.00969395, 51.02204645, 51.03532721, 51.04490968,
   51.01566, 51.02015696, 51.02585864, 51.01347, 51.02003, 51.0233,
   50.99115206, 51.0007123, 51.0147602, 51.04182475, 51.03110881)

long = c(1.27062894, 1.305595503, 1.251393586, 1.3330022, 1.282221582,
   1.327234, 1.238538388, 1.279034619, 1.258740908, 1.2914384, 1.31293,
   1.244711054, 1.239145274, 1.30674, 1.3207, 1.3246, 1.251393586, 1.276897481,
   1.299778301, 1.266223243, 1.247370643)

#####
##### Maximum likelihood estimates
#####

library(MASS)

ast.all = c(S08, C11, S01, R06, C17, R04, C02, C09, C21, C26, C05, C07,
            C14, R02, R07, R08, S01, S06, S13, S15, S17)

#### Multiply by 6.967 as this converts to per 100m

ast.nb = fitdistr(ast.all, "Negative Binomial")

#### Check formulae in the paper

ast.po = fitdistr(ast.all, "Poisson")
paper.po.se = sqrt(sum(ast.all))/length(ast.all)
paper.po.se; ast.po                #### correct

paper.nb.se = sqrt(ast.nb$est[2]*(ast.nb$est[1]+ast.nb$est[2]) /
       (length(ast.all)*ast.nb$est[1]))
paper.nb.se; ast.nb                #### correct


##### Maximum likelihood estimate
ast.nb$est[2] * 6.967

##### 95% confidence intervals
6.967*(ast.nb$e[2] - 1.96*ast.nb$s[2])  # mu - 1.96*se(mu)
6.967*(ast.nb$e[2] + 1.96*ast.nb$s[2])  # mu + 1.98*se(mu)
##### (1.89,4.27)


#####*******************************************************************
##### GVFCMOM Point estimates
#####*******************************************************************

# Need to calculate a contant C so that estimates are per 100m^2
# This needs to be a function of the number of minute-long trawl segments
# because the video trawls are not all the same length.
# 1 knot=0.514444 metres per second
# speed of video trwal is 0.5 knots
# width of video window approximately 0.93m

const = function(S) {
#**************************************************************
# Calculates C to convert mean figures to mean per 100m^2 as a
# function of the number of segments S
#**************************************************************
distperminute = 0.514444*60*0.5   # 15.4333
areaperminute = 0.93*distperminute
C = 100 / (areaperminute*S)
C}

# S=7 for C26
# S=9 for C11 and R06
# S=10 for the rest of the stations

con7 = const(7)
con9 = const(9)
con10 = const(10)

S08.est = con10 * GVFCMOM(S08, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
C11.est = con9 * GVFCMOM(C11, s=9, d=2, method='nb', k=0.42, lowint=0, highint=100)
S01.est = con10 * GVFCMOM(S01, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
R06.est = con9 * GVFCMOM(R06, s=9, d=2, method='nb', k=0.42, lowint=0, highint=100)
C17.est = con10 * GVFCMOM(C17, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
R04.est = con10 * GVFCMOM(R04, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
C02.est = con10 * GVFCMOM(C02, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)

C09.est = con10 * GVFCMOM(C09, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
C21.est = con10 * GVFCMOM(C21, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
C26.est = con7 * GVFCMOM(C26, s=7, d=2, method='nb', k=0.42, lowint=0, highint=100)
C05.est = con10 * GVFCMOM(C05, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
C07.est = con10 * GVFCMOM(C07, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
C14.est = con10 * GVFCMOM(C14, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
R02.est = con10 * GVFCMOM(R02, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)

R07.est = con10 * GVFCMOM(R07, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
R08.est = con10 * GVFCMOM(R08, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
S01.est = con10 * GVFCMOM(S01, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
S06.est = con10 * GVFCMOM(S06, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
S13.est = con10 * GVFCMOM(S13, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
S15.est = con10 * GVFCMOM(S15, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)
S17.est = con10 * GVFCMOM(S17, s=10, d=2, method='nb', k=0.42, lowint=0, highint=100)

momests = c(S08.est, C11.est, S01.est, R06.est, C17.est, R04.est, C02.est,
            C09.est, C21.est, C26.est, C05.est, C07.est, C14.est, R02.est,
            R07.est, R08.est, S01.est, S06.est, S13.est, S15.est, S17.est)

par(pty='s')

####
#### "momests" contains the GVFCMOM estimates. Write them to a file for
#### plotting on the map
#### 

results = cbind(long, lat, momests)
write.csv(results, file = "folkestone results.csv", row.names=FALSE)
 
tran = sqrt(momests / pi)