# Analysis of Dispersion according to C. R. Rao (1965)  Linear Statistical 
#  Inference and Its Applications.  John Wiley & Sons, New York, 522pp.
# 1) function mlreg revised from brewScor.R subroutine which solved for 1 Y.
# 2) Input of a test file from Steele & Torrie Table 15.12
#    a. I can not extract the matrix data from the object YX (2/4/2003)
#    b. Now I can (2/5/2003)
# 3) Make it possible for fn mlreg to print out a particular matrix (e.g. the design matrix X'X)
# 4) Extract the names of variables from the attributes of yx, the input matrix.
# 5) Create a unique time based name to store files in using Sys.time fn.
# 6) Store some files X'Xi, X'Y, Y'Y for further use using:
 #    if(capabilities("fifo")) {
 #      zz <- fifo("foo", "w+")
 #      writeLines("abc", zz)
 #      print(readLines(zz))
 #      close(zz)
 #      unlink("foo")
 #    }
# 7) Start calculating the reduction due to a factor.
 #   Factor effects are calculated. Now their significance must be established.

# xcelnum<- read.csv('edperscore.csv', header = TRUE, sep = ",")
require(sm)
require(MASS)

# define linear regression solving routine for multiple Ys and multiple Xs
# YX = Data Y appended to a design matrix X
# NY identifies the number of Y columns in the data set then the rest are assumed to be X
# PX tells the function to print out the X'X matrix
# there are ways to temporarily store results so that it will be available to routines outside the fn
#
mlreg<-function(yxtyx,NY,PX){
 bx<- NY+1
 nc<- ncol(yxtyx)
# yxtyx<-  t(YX)%*% YX 
 Syy<-yxtyx[1:NY,1:NY] 
 Sxy<-yxtyx[bx:nc,1:NY] 
 Sxx<-yxtyx[bx:nc,bx:nc]
 if (PX==1) {cat("Design Matrix, X'X:\n"); print(Sxx)} 
 xtxi<-ginv(Sxx)  # G inverse fn from require(MASS)
# Define calculation of R-squared as alternate output
 B<- xtxi %*% Sxy
# Rsq<- 1 - (Syy - t(Sxy) %*% B)/(Syy-t(Sxy[,1])%*%Sxy[,1]/nrow(YX))
 syxtyx<-yxtyx
 syxtyx[bx:nc, 1:NY]<-B 
 syxtyx[bx:nc,bx:nc]<- xtxi
 syxtyx 
}

FLAM<- function(L, N, P, Q, M) {
X<-1
}

# -------------------------------------------------- 52
# End of Function deffinitions

# Begin Main Program
# Data input
 cat("Sys.time:", btime<- Sys.time(), '\n')
# yx<- read.csv('ST15d12.csv', header = TRUE, sep =",")
 yx<- read.csv('LuisSmoothed.csv', header = TRUE, sep =",")
 attach(yx)
 ny <- 2

#
# matrify the input data
 YX<- matrix(0, NR<-nrow(yx), NC<-ncol(yx), byrow=FALSE)
 for (i in 1:NC) { for (j in 1:NR) {YX[j,i]<- (yx[[i]][j])}}
 cat("Day Label:", dlabel<- format(Sys.time(), "%Y%b%d"),'\n')
 cat("DayTime Label:", tlabel<- format(Sys.time(), "%a%X"),'\n')
yxtyx<-  t(YX)%*% YX
tSyy<- yxtyx[1:ny, 1:ny] 
#
# Solution of Overall Matrix Equation and extraction of Solution Beta
 Beta<- (SYXtYX<-mlreg(yxtyx, ny, 1))[(ny+1):NC,1:ny]
 cat("Estimated Parameter Matrix, Beta\n")
 print(Beta)
 hisdat<-Beta[2:nrow(Beta),]
 YXnames<-attributes(yx)[1]
 Ynames<-YXnames[[1]][1:ny]
 cat("Y Variable names:\n")
 print(Ynames)
 Xnames<-YXnames[[1]][(ny+1):NC]
 noU<-Xnames[2:nrow(Beta)]
 cat("Design column names:\n")
 print(Xnames)
#
# barplot of Solution Beta (minus u estimation).
 barplot(t(hisdat),beside=TRUE, 
   xlab="FACTORS",
   col = c("red", "green", "blue", "orange", "purple", "brown", "black")[1:ny],
   names.arg=noU, main="Estimated Parameters, B", 
   legend = Ynames,
   sub="color identifies the trait affected by a factor")
# Normalize factors within a trait
 maxB<- 0*(1:ny)
 for (i in 1:ny) maxB[i]<- max(abs(hisdat[,i]))
 Nhisdat<- matrix(maxB, nrow(hisdat), ncol(hisdat), byrow=TRUE)
 Nhisdat<- hisdat/Nhisdat
pause() # between simple and normalized solution output.
 barplot(t(Nhisdat),beside=TRUE, 
   xlab="FACTORS",
   col = c("red", "green", "blue", "orange", "purple", "brown", "black")[1:ny],
   names.arg=noU, main="Normalized Estimated Parameters, B", 
   legend = Ynames,
   sub="color identifies the trait affected by a factor")
#
# Factor solutions: 105
#
iSyy<- matrix(0,ny,ny)
Ksol<- ask("How many direct solutions are to be calculated?")
DiSyy<- array(0,c(ny,ny,Ksol+2))
cat("Big Array defined:\n")

df<- 0*(1:(Ksol+2))
cat("df vector defined:\n")
df[1]<- NR

DiSyy[1:ny, 1:ny, 1]<- tSyy 
for (i in 1:Ksol) {
  namsol<- ask(paste("Enter name of factor", i, "to be solved"))
  isol<- ask(paste("Enter the indices for", namsol))
  df[i]<- length(isol)
  idim<- ny + df[i]
  iYXtYX<- matrix(0, idim, idim)
  index<- c(1:ny, ny+isol) 
  iYXtYX[1:idim, 1:idim]<- yxtyx[index, index]
   ir<- 1:ny
   ic<- (ny+1):ncol(iYXtYX)
   iBeta<- (iSYXtYX<-mlreg(iYXtYX, ny, 0))[ic, ir]
   iXtY<- iSYXtYX[ir, ic]
   iYtXB<- matrix(iXtY, length(ir), length(ic)) %*% iBeta
   DiSyy[1:ny, 1:ny, i+1 ]<- iYtXB
   iSyy<- iSyy + iYtXB
  }

DiSyy[1:ny, 1:ny, Ksol+2]<- tSyy-iSyy 
cat("Array SS:\n")
 print(DiSyy)
# after the loops: 123
cat("Total SS:\n")
 print(tSyy)  
cat("Design SS:\n")
 print(iSyy)  
cat("Residual SS:\n")
 print(tSyy-iSyy)  
cat("Wilkes Lambda:\n")
 print(Lambda<- det(tSyy-iSyy)/det(iSyy) ) 
cat("Chi Square:\n")
 P<- ny
 Q<- NC-ny
 N<- NR
 M<- N - 0.5*(P+Q+1)
 print(ChiSq<-  (-M*log(Lambda)))
 
 #P <- ny
 #Q <- NC-ny
 #N <- NR
 #M<- N - 0.5*(P+Q+1)
 #STATS<- c(Lambda, (-M*log(Lambda)), P*Q, N, P, Q, M, FLAM(Lambda,N,P,Q,M))
# print(STATS)
# Define Wiles Lambda function
# -----------------------------------------------------
# cat("Wilkes Lambda:\n")
# m<- NR - 
# print(-m*log(Lambda))
#
# J calculation of statistics over Wilkes lambda  
#NB. Computation of Lambda and appending of Chi-Sq and F stats with FLAM
# LAMB=: 4 :0
# R0=. x.
# RV=. y.
# lambda=: (det R0)% det R0 + RV
# N=: RDF + Q=: $,V2
# P=: 1}.$R0
# if. (lambda > 0) do. goto_L2. else. goto_L1. end.
# label_L1. MESS=: MESS, 'lambda negative, '
# lambda=: 1
# label_L2. lambda=: 3 3 $(_1 }. lambda,(-M*^.lambda),(P*Q),N,P,Q,M=: N - 0.5 * P+Q+1), FLAM lambda
# )

# End of program and timestamp:
 cat("Sys.time:",etime<- Sys.time(), '\n')
 cat("Run Time:",etime-btime, 'secs.\n')