{

#code to do Streamflow disaggregatio..
# The monthly flows (XX) are transformed into an orthonormal space YY)
#  where they are all perpendicular to the plane X1+X2....+Xd = Z (the annual flow)
#  The rotation is obtained by the matrix qqv YY = qqv %*% XX
#   the result of this transformation is that the last column of YY = Z' = Z/sqrt(nd)
#   The rest of the vector is defined as U

#  Then the disaggregation is just resampling from the conditional pdf f(U | Z')


# Use Tarboton et al. (1997) paper for notations. This replaces their approach by
#  the K-NN bootstrap.



data=matrix(scan("Leesferry-mon-data.txt"),ncol=13,byrow=T)

#WYmonflow = matrix(scan("0725RebuiltNatFlow.txt"),ncol=12,byrow=T)
#converts ac-ft/month to cms
WYmonflow=data[,2:13]
WYmonflow = WYmonflow*.0004690502        

X=WYmonflow

#X=X[,1:2]              #example..


# test it for 3 variables...

n=length(X[,1])
nd=12
nd1=nd-1


XX=X[,1:nd]

dd=diag(1,nd,nd)

qqv=matrix(0,ncol=nd,nrow=nd)

qqv[,nd]=1/sqrt(nd)


for(i in (nd-1):1){

j1=i+1
yv=rep(0,nd)

for (j in j1:nd){
yv=yv+(qqv[,j] %*% dd[,i])*qqv[,j]}
xv=dd[,i]-yv
xv=xv/sqrt(sum(xv*xv))
qqv[,i]=xv

}

Rmat=t(qqv)

#another way of getting the R matrix..
qqv=diag(1,12)
qqv[,1]=1/sqrt(nd)
RR=qr.Q(qr(qqv))
xx=as.array(-RR)
 xy=rev(xx)
 xz=matrix(xy,ncol=12)
Rmatalt=t(xz)


YY=t(Rmat %*% t(XX))	#the last column of YY will be Z' = Z/sqrt(nd)

annflow=apply(XX,1,sum)	#the annual flow..

Zprime = annflow/sqrt(nd)


#generate an  aggregate flow = annflow, say xp

#get K-nearest neighbors of xp in the annflow vector

#Pick a neighbor with the probabilit metric W (i.e. 1/K weight function or distance..)


nsim=n
xsim=matrix(0,ncol=nd,nrow=nsim)

# for demo - just do a bootstrap of annflows

kk=round(sqrt(n))
W = 1:kk
W = 1/W
W = W/sum(W)
W = cumsum(W)

zz=round(runif(nsim,1,n))
annsim=annflow[zz]

for(i in 1:nsim){
xp=annsim[i]
xd=abs(xp-annflow)
xz=order(xd)

#select first K neighbors..
xz=xz[1:kk]

#pick a neighbor out of the K neighbros based on the weight
#metric W
xd = runif(1,0,1)
xy = c(xd,W)
xd = rank(xy)
i1 = xz[xd[1]]

zpsim=xp/sqrt(nd)
uu1=c(YY[i1,1:nd1],zpsim)


xpred=t(Rmat) %*% uu1
xsim[i,]=xpred

}

#check to see if the sum of simulated monthly flows are equal to the annual flows
# that went into the disagg.

zx=apply(xsim,1,sum)
  zy=cbind(zx,annsim)
  write(t(zy),file="test.txt",ncol=2)


}