options(echo = FALSE);

compute_SBfactors <- function( B1seriesFilename = "",
                               SBfactors = c(),
                               SByears = c() ) {

########################################################################
##
## PURPOSE:
## --------
##     R function to rescale input seasonal break factors such that they
##     'perfectly balance' an input B1 series after application.
##
##     By Frank Masci; Time Series Analysis Section
##
## VERSION; LAST MODIFIED:
## -----------------------
##     Vsn 3.0; 24/9/2006 
##
## INPUT ARGUMENTS:
## ----------------
##     B1seriesFilename - filename of B1 series from SEASABS in format
##                        "YR  MNTH  VALUE" before SB corrections are applied.
##     SBfactors - vector  of approximate SB priors from S*I chart for each
##                 period (set equal to 1 if you don't want to apply a prior).
##     SByears - vector of corresponding "break years" for each period. For a
##               given period j, equals last year before which <S*I> level changes
##               abruptly. Note if SBfactors[j]=1, then SByears[j] is ignored. 
##
## OUTPUT:
## -------
##     Table of rescaled seasonal break factors for each 'break year' for
##     input into SEASABS.
##
## EXAMPLE:
## --------
##     The following function call reads in some approximate SB factors derived
##     apriori from the S*I chart for March and August at break years 2001/2002
##     and 1999/2000 respectively and rescales them such that the B1 total will
##     remain the same after application. Note that SB factors of unity are not
##     rescaled. Execute the following in R after "sourcing" this file, e.g:
##
##     source("S:\R\TSAlib\compute_SBfactors.r");
##     compute_SBfactors("B1_series_beforeSB.txt",
##                       c(1,1,1.26,1,1,1,1,0.76,1,1,1,1),
##                       c(2001,2001,2001,2001,2001,2001,2001,1999,2001,2001,2001,2001));
##
########################################################################

# Parse inputs and perform some sanity checks.

s <- scan(B1seriesFilename, list(0,0,0));
years <- s[[1]];
periods <- s[[2]];
values <- s[[3]];

freq <- 4;
for(j in 1:12) {
  if(periods[[j]] == 12) {
    freq <- 12;
  }
}

periodstring <- seq(1,freq,by=1);

if( length(SBfactors) != length(periodstring) ) {
  write(paste("\n\n*** Error: length of \"SBfactors\" vector wrong; quitting...\n"),file="");
  q();
}

if( length(SByears) != length(periodstring) ) {
  write(paste("\n\n*** Error: length of \"SByears\" vector wrong; quitting...\n"),file="");
  q();
}

#-------------------------------------------------------------------
# Determine start/end years for each period.

startyear <- 0;
Nsb <- 0;
for(j in 1:freq) {
  if( SBfactors[[j]] != 1 ) {
    Nsb <- Nsb + 1;
  }
  for(i in 1:length(values)) {
    if( periods[[i]] == j ) {
      startyear[[j]] = years[[i]];
      break;
    }
  }
}

revperiods <- periods[c(seq(length(values), 1, by=-1))];
revyears <- years[c(seq(length(values), 1, by=-1))];

endyear <- 0;
for(j in 1:freq) {
  for(i in 1:length(values)) {
    if( revperiods[[i]] == j ) {
      endyear[[j]] = revyears[[i]];
      break;
    }
  }
}

# sanity check that elements of input SByears vector are in range.

for(i in 1:freq) {
  if( (SByears[[i]] - startyear[[i]]) <= 0 ) {
    write(paste("\n\n*** Error: break year less than or equal to start year for period",i, "; quitting...\n"),file="");
    q();
  }
  if( (SByears[[i]] - endyear[[i]]) >= 0 ) {
    write(paste("\n\n*** Error: break year greater than or equal to end year for period",i,"; quitting...\n"),file="");
    q();
  }
}

if( Nsb == 0 ) {
  write(paste("\n\n*** Error: all SB factors are = 1. Are you sure you want this? Quitting...\n"),file="");
}

#-------------------------------------------------------------------
# Compute B1 totals before/after applying seasonal break factors.
# For B1totalAfter, break each period i into: for years <= SByears[i]
# and then for > SByears[i].

SBfactors <- (1 / SBfactors);
B1totalBefore <- 0;
B1totalAfter <- 0;
B1totalAfterLT_SByears <- 0;

for(j in 1:freq) {
  B1totalBefore[[j]] = 0;
  B1totalAfter[[j]] = 0;
  B1totalAfterLT_SByears[[j]] = 0;

  for(i in 1:length(values)) {
    if(periods[[i]] == j) {
      B1totalBefore[[j]] = B1totalBefore[[j]] + values[[i]];
      if( years[[i]] <= SByears[[j]] ) {        
        B1totalAfter[[j]] = B1totalAfter[[j]] + SBfactors[[j]]*values[[i]];
        B1totalAfterLT_SByears[[j]] = B1totalAfter[[j]];
      } else {
        B1totalAfter[[j]] = B1totalAfter[[j]] + values[[i]];
      }
    }
  }
}

#-------------------------------------------------------------------
# Apportion excess/deficit in B1 after SB correction equally amongst
# all periods with SB factor != 1. Compute correction factors to input
# SBfactors != 1.

DiffperPeriod <- (sum(B1totalBefore) - sum(B1totalAfter)) / Nsb;

g <- 0;
for(j in 1:freq) {
  if( SBfactors[[j]] != 1 ) {
    g[[j]] = (DiffperPeriod / B1totalAfterLT_SByears[[j]]) + 1;
  } else {
    g[[j]] = 1;
  }
}

finalSBfactors <- g * SBfactors;

#-------------------------------------------------------------------
# Check B1 is balanced with final factors.

B1totalAfterFinal <- 0;
for(j in 1:freq) {
  B1totalAfterFinal[[j]] = 0;
  for(i in 1:length(values)) {
    if(periods[[i]] == j) {
      if( years[[i]] <= SByears[[j]] ) {        
        B1totalAfterFinal[[j]] = B1totalAfterFinal[[j]] + finalSBfactors[[j]]*values[[i]];
      } else {
        B1totalAfterFinal[[j]] = B1totalAfterFinal[[j]] + values[[i]];
      }
    }
  }
}

#-------------------------------------------------------------------
# Outputs:

write(paste("\nTotal B1 BEFORE SB correction = ",sum(B1totalBefore)),file="");
write(paste("Total B1 AFTER SB correction = ",sum(B1totalAfterFinal)),file="");
write(paste("\nPeriod / Break years / SB Factors:"),file="");
write(paste("-----------------------------------------"),file="");
SByearPlus1 <- SByears + 1;
write(paste("  ",periodstring," ",SByears,"/",SByearPlus1,"  ",1/finalSBfactors),file="");
write(paste("-----------------------------------------\n"),file="");

}
#------------------------------END-----------------------------------

# Examples, execute as:
# "S:\\R\\rw2001\\bin\\Rterm.exe" -q --no-restore --no-save < compute_SBfactors_v3.r

compute_SBfactors("B1_Agriculture_beforeSB.txt",
                  c(1,1,1.26,1,1,1,1,0.76,1,1,1,1),
                  c(2001,2001,2001,2001,2001,2001,2001,1999,2001,2001,2001,2001));

compute_SBfactors("B1_withSB_noiseless.txt",
                  c(0.96,0.96,0.96,0.96,0.96,1.362,0.96,1,0.96,0.96,0.96,0.960),
                  c(1989,1989,1989,1989,1989,1989,1989,1989,1989,1989,1989,1989));