CRM-Pr. added.R

aalatt3 — Aug 15, 2013, 6:55 PM

#Bismilla-Reading the files and preparing the matrices------------CRM Model
ptm <- proc.time()
ms=4
PRD=read.csv(file="PRDCRM.csv",)
INJ=read.csv(file="INJCRM.csv")
PRD=do.call(cbind, PRD[ms])
PR=do.call(cbind, read.csv(file="PR.csv"))
INJ=do.call(cbind, INJ)
taoi.p=.2
taop=do.call(cbind, read.csv(file="taop.csv"))
#Preparing the matrices
pp.matrix=matrix(ncol=ncol(PRD), nrow=(nrow(PRD))) #pp matrix
pp.var=matrix(ncol=ncol(PRD), 1) #pp variance matrix
cpp.I=matrix(ncol=ncol(PRD), nrow=ncol(INJ)) #cpp.I matrix
inj.inj=matrix(ncol=ncol(INJ), nrow=ncol(INJ)) #inj-inj matrix
pp.q=matrix(ncol=ncol(PRD))  #pp.q matrix
ci.qj=matrix(ncol=ncol(PRD), nrow=ncol(INJ))   #ci.qj matrix
lhs=matrix(nrow=ncol(INJ)+2,1)   #lhs matrix
lambda=matrix(ncol=1, nrow=(1+ncol(INJ)))    #lambda matrix
master=matrix(ncol=ncol(INJ)+2, nrow=(ncol(INJ)+2))   #master matrix
i.avg=matrix(1,ncol(INJ))
i.dash=matrix(ncol=ncol(INJ), nrow=(nrow(INJ)))

#######################preparing the BHP matrix########################
##Preparing the 1st term (Pwf(n0)e)
bhp.1term=matrix(ncol=ncol(PR), nrow=nrow(PR))
for ( k in 1: ncol(PR))
  for (i in 1: nrow(PR)){
    bhp.1term[i,k]=PR[1,k]*exp((1-i)/taop[k,1])
  }
##Preparing the 2nd  term (Pwfkj) 
bhp.2term=matrix(ncol=ncol(PR), nrow=nrow(PR))
bhp.2term=PR
## Preparing the 3rd term(P'wfkj)





#############Preparing the i.dash values####################################
# preparing the i' values
taoi=do.call(cbind, read.csv(file="taoi.csv",header=FALSE))
for (i in 1:ncol(INJ)){
  for (j in 1:nrow (INJ)){
    temp=0
    for (k in 1:j){
      temp=(1/taoi[i,1])*exp((k-j)/taoi[i,1])*INJ[k,i]+temp
    }
    i.dash[j,i]=temp
  }
}

########################filling PP MATRIX FOR ALL PRODUCERS###################
for (j in 1 : ncol(PRD))
  for (i in 1: nrow(PRD))  {
    pp.matrix[i,1]=PRD[1,1]*exp((-i+1)/taoi.p)
  }
#######################filling THE COV. MATRIX FOR PP#########################
pp.var= var(pp.matrix)

#######################filling THE COV OF Cpp-I and its transpose#############
for (j in 1:ncol(cpp.I))
  for (i in 1: nrow(cpp.I)){
    cpp.I[i,j]=cov(pp.matrix[,1],i.dash[,i] )
  }
cpp.IT=t(cpp.I)
#######################filling THE INJ-INJ COV MATRIX##########################
for (j in 1: ncol(inj.inj))
  for (i in 1: ncol(inj.inj)){
    inj.inj[i,j]=cov(i.dash[,j], i.dash[,i])
  }
#######################filling SIGMA.PP-QQ MATRIX#############################
for ( j in 1: ncol (pp.q)){
  pp.q[1,j]=cov(pp.matrix[,j], PRD[,j])
}
######################filling CI-QJ MATRIX####################################
for ( j in 1: ncol(ci.qj))
  for (i in 1: nrow(ci.qj)){
    ci.qj[i,j]=cov(i.dash[,i], PRD[,j])
  }

################Preparing the avg. matrix for i.dash###########################
i.avg=colMeans(i.dash)
###################ASSEMBLING THE MASTER MATRIX################################
master[1,1]=pp.var[1,1]
master[1,2:(ncol(master)-1)]=cpp.IT[1,]
master[1,ncol(master)]=mean(pp.matrix)
master[2:(nrow(master)-1), 1]=cpp.I[,1]
master[nrow(master),1]=mean(pp.matrix)
master[2:(nrow(master)-1),2:(ncol(master)-1)]=inj.inj
master[nrow(master),2:(ncol
                       (master)-1)]=t(i.avg)
master[nrow(master),ncol(master)]=0
master[2:(nrow(master)-1), ncol(master)]=i.avg
##################################filling THE LHS Matrix#####################
lhs[1,1]=pp.q
lhs[2:(nrow(lhs)-1),1]=ci.qj[,1]
lhs[nrow(lhs),1]=mean(PRD)
##############lambda###################Resulsts##############################
lambda=solve(master)%*% lhs 
#for (i in 1 : nrow(lambda)){
if(lambda[i,1]<0) {lambda[i,1]=0}
#}
write.csv(lambda,file="lambda.csv")
#############################plotting########################################
coff=read.csv(file="BETAS.CSV")

png(file=paste("CRM Lambdas",".png"),width=900,height=900)
plot (lambda[2:nrow(lambda),1], type="l",xlab="Injector No.", 
      ylab="Lambda", col="green",axes=FALSE ,
      ylim=c(min(lambda,coff[,2] ),c(max(lambda, coff[,5]))))
axis(2)
axis(1,at=1:(ncol(INJ)), cex.axis=.7)
par(new=TRUE)
plot (coff[,ms+1], axes=FALSE, type="l", col="red",xlab=FALSE,ylab=FALSE)
#axis(4)
box()
title(main = (paste("Lambdas", lambda[1,1])))
dev.off()

pdf 
  2 

###############checking out#############################

# Preparing the lambda X i' matrix
lambda.i=matrix(ncol=ncol(INJ),nrow=nrow(INJ))
for (i in 1: ncol(INJ)){
  lambda.i[,i]=lambda[i+1, 1]*i.dash[,i]
}

#Wrapping up to calculate Qhat
#preparing the on col qhat

q.hat=matrix(nrow=nrow(INJ),1 )
for (i in 1:nrow(INJ)){
  q.hat[i,1]=sum(lambda.i[i,1:ncol(INJ)])+lambda[1,1]*pp.matrix[i,1]

}
diff=matrix(nrow(INJ), 1)
diff=q.hat-PRD
png(file=paste("Qobs vs Qhat",".png"),width=2000,height=900)
plot(PRD, type="l", col="red", axes=TRUE,ylim=c(min(diff,PRD,q.hat),max(PRD,q.hat)))

legend ('topright', c("Q Calc." ,"Q Obs","Difference"),lty=c(1,1,2),
        col=c("blue", "red", "brown"))

par(new=T)
plot(q.hat[,1], type="l", col="blue", axes=F)
par(new=T)
plot(diff, type="b", col="brown", axes=F,
     ylim=c(min(diff,PRD,q.hat),max(PRD,q.hat)))
dev.off()

pdf 
  2 

Difference=colSums(q.hat)-colSums(PRD)
proc.time() - ptm

   user  system elapsed 
   0.49    0.05    0.56