/**************************/
/*  Load R-libraries used */
/**************************/

library("boot")
library("class")
library("foreign")
library("MASS")
library("sp")
library("maptools")
library("msm")
library("nnet")
library("lattice")
library("Matrix")
library("spatial")
library("SparseM")
library("spam")
library("fields")
library("tripack")
library("tools")
library("deldir")
library("nlme")
library("coda")
library("spdep")

/********************************************/
/*   Simulate noisy image on a grid         */
/********************************************/

/* define Moore's Basis functions: 
s=range of angular distance: q/2pi
n=number of basis functions
d=Dmax
*/

b<-function(s,n,d){
        B<-(d/2)*(((sin(pi*n-2*pi*d*s)/(pi*n-2*pi*d*s))-(sin(pi*n+2*pi*d*s)/(pi*n+2*pi*d*s))))
        return(B)
}

/*generate noiseless image intensities on a grid  

s=angular distance:q/2pi
grid size: 100x100
beam-center: (50,50)
s range: (0.001,0.075) or q range: (0.006,0.47)
a<-c(5.281465e-05, 4.178242e-05, 1.686120e-05 ,5.454271e-06, 3.379092e-06, 3.007528e-06, 2.055125e-06, 8.427587e-07)      # 8 Moore's coefficients
d<-58.5     # Dmax

*/
    
x1<-seq(0.001,sqrt(2)*0.075,length=100)     
y1<-seq(0.001,sqrt(2)*0.075,length=100)             # x1, y1 defines the grid corresponding to the s range
S<-matrix(0,ncol=length(x1), nrow=length(y1))       # matrix to store the s distance from beam-center
K<-matrix(0,ncol=length(x1), nrow=length(y1))  
I_data<-matrix(0,ncol=length(x1), nrow=length(y1))  # Noiseless intensity matrix


for (i in 1:length(x1))
      { 
        for(j in 1:length(y1))
         { 
           S[i,j]<-sqrt((x1[i]-round(x1[50],3))^2+(y1[j]-round(x1[50],3))^2)
           K[i,j]<-(a[1]*b(S[i,j],1,d)+a[2]*b(S[i,j],2,d)+a[3]*b(S[i,j],3,d)+a[4]*b(S[i,j],4,d)+a[5]*b(S[i,j],5,d)+a[6]*b(S[i,j],6,d)+a[7]*b(S[i,j],7,d)+a[8]*b(S[i,j],8,d))
           I_data[i,j]<-((K[i,j])/S[i,j])*100000
        }
       }


/*generate the noisy image intensities */

nbd<-cell2nb(100,100,type="queen")   # setup spatial neighorhood structure of image pixels: queen type nbd chosen
N<-nb2listw(nbd,style="W")
NN=similar.listw(N)
W=as_dsTMatrix_listw(NN)
A=as_dsCMatrix_IrW(W,rho=0.5)
AA=as(A,"CsparseMatrix")
Dinv=as(Diagonal(x=1/I_data),"CsparseMatrix")
X<-(Dinv%*%AA)
XpX<-t(X)%*%X
error=rnorm(10000,mean=0,sd=1)
Xpe<-t(X)%*%error
sol<-as.vector(solve(XpX,Xpe))
E<-matrix(sol,nrow=100)          #error matrix

data=I_data+E    # final noisy image


/***************************************/

/*     Analysis of a pixilated image   */

/***************************************/

/* likelihood function for the spatial independent case */

lpostf<-function(a1,a2,a3,a4,a5,a6,a7,a8,y=data,cxx,cyy,d){
              cx=floor(cxx)
              cy=floor(cyy)
              del_cx=cxx-cx
              del_cy=cyy-cy
                           
              I<-rep(0,10000)
              res<-rep(0,10000)              
              S<-rep(0,10000)
              Kfunc<-rep (0,10000)
              x1<-seq(0.001,sqrt(2)*0.075,length=100)
              y1<-seq(0.001,sqrt(2)*0.075,length=100)
              
              g=expand.grid(x=x1,y=y1)
              S<-sqrt((g$x-x1[cx]+0.000001-del_cx*0.001061273)^2+(g$y-y1[cy]+0.000001-del_cy*0.001061273)^2)

              a<-c(a1,a2,a3,a4,a5,a6,a7,a8)/100000
              Kfunc<-(a[1]*b(S,1,d)+a[2]*b(S,2,d)+a[3]*b(S,3,d)+a[4]*b(S,4,d)+a[5]*b(S,5,d)+a[6]*b(S,6,d)+a[7]*b(S,7,d)+a[8]*b(S,8,d))
              U<-ifelse((Kfunc > 0)==T,Kfunc,0.00001)
              I<-(U*100000)/S
              res<-(y-I)
              Vhat=as.vector(I)           

             llike=(sum(((res^2)/(2*Vhat)),na.rm=T))+0.5*sum(log(Vhat),na.rm=T)
             
              prior_part<-sum(log(a))+sum((0.22*(log(a[1:8])-1)^2))+2*d-109*log(d)+(((cxx-50)^2)/8)+(((cyy-50)^2)/8)
              
               logpost<-(-(prior_part+llike)) 
              return(logpost)
             }

/* MCMC updates of all parameters except the autocorrelation 

NUMIT: number of iterations of MCMC
s,s1,s2,s3: std errors used to increment parameter values for MCMC updates
*/



MHsampler = function(NUMIT,s,s1,s2,s3)
  {
    n = length(as.vector(data))
    cat("n=",n,"\n")

    ## NUMIT   matrix to store Markov chain values
    ## each row corresponds to one of 11 parameters in order: a's,cx,cy,d
    ## each column corresponds to a single state of the Markov chain
    mchain = matrix(NA,15,NUMIT)
    acc_a12d=0         # acceptance rates for parameter updates
    acc_a3=0
    acc_a4=0
    acc_a5=0
    acc_a6=0
    acc_a7=0
    acc_a8=0
    acc_cx=0
    acc_cy=0
    
    ## starting values for Markov chain
        
    log_curr<-lpostf(5.26,4.18,1.70,0.55,0.34,0.30,0.21,0.09,data,50.00,49.99,58.71) # initial call to evaluate logl


    rho1d=-0.99    ### pairwise correlations for joint updates of a1,a2 and Dmax
    rho2d= 0.93
    rho12=-0.93
    
    
    mchain[,1] =c(5.26,4.18,1.70,0.55,0.34,0.30,0.21,0.09,50.00,49.99,58.71,log_curr,rho1d,rho2d,rho12) #initial parameter values
    
     mu12d=c(mchain[1,1],mchain[2,1],mchain[11,1])
     sa1=0.03
     sa2=0.007
     Sigma12d=matrix(c(sa1^2,rho12*sa1*sa2,rho1d*sa1*s2,
       rho12*sa1*sa2,sa2^2,rho2d*sa2*s2,rho1d*sa1*s2,rho2d*sa2*s2,s2^2),3,3)   ### 3 D correlation matrix for joint updates
     
    for (i in 2:NUMIT)
      {
       
        ## most upto date state for each parameter
        curra1 = mchain[1,i-1]
        curra2 = mchain[2,i-1]
        curra3 = mchain[3,i-1]
        curra4 = mchain[4,i-1]
        curra5 = mchain[5,i-1]
        curra6 = mchain[6,i-1]
        curra7 = mchain[7,i-1]
        curra8 = mchain[8,i-1]        
        currcx = mchain[9,i-1]
        currcy = mchain[10,i-1]
        currd = mchain[11,i-1]
        
        
        cat(" Iteration:",i)


         ## sample from full conditional distribution of theta  ( joint Metropolis-Hastings update)
        propa12d = mvrnorm(1,mu12d,Sigma12d) # draw one sample at random )
        propa1=propa12d[1]
        propa2=propa12d[2]
        propd=propa12d[3]

        ## Metropolis accept-reject step (in log scale)
        log_prop<-lpostf(propa1,propa2,curra3,curra4,curra5,curra6,curra7,curra8,data,currcx,currcy,propd)
        logMHratio1=log_prop-log_curr
 
        logalpha1 = min(0,logMHratio1) # alpha = min(1,MHratio)
        if ((log(runif(1))<(log(1)+logalpha1)) ) # accept if unif(0,1)<alpha, i.e. accept with probability alpha, else stay at current state
          {
            acc_a12d = acc_a12d + 1 # increment count of accepted proposals

            curra1 = propa1
            curra2 = propa2
            currd = propd
            log_curr<-log_prop
           }

        ## sample from full conditional distribution of theta (Metropolis-Hastings update)
        propa3 = rtnorm(1,mean=curra3,sd=s,0,Inf) # draw one sample at random from exp )

        ## Metropolis accept-reject step (in log scale)
        log_prop<-lpostf(curra1,curra2,propa3,curra4,curra5,curra6,curra7,curra8,data,currcx,currcy,currd)
        logMHratio3=log_prop-log_curr+log(1-pnorm(-curra3/s))-log(1-pnorm(-propa3/s))
                     

        logalpha3 = min(0,logMHratio3) # alpha = min(1,MHratio)
        if ((log(runif(1))<(log(1)+logalpha3))) # accept if unif(0,1)<alpha, i.e. accept with probability alpha, else stay at current state
          {
            acc_a3 = acc_a3 + 1 # increment count of accepted proposals
            curra3 = propa3
            log_curr<-log_prop
           }

        ## sample from full conditional distribution of theta (Metropolis-Hastings update)
        propa4 = rtnorm(1,mean=curra4,sd=s,0,Inf) # draw one sample at random)

        ## Metropolis accept-reject step (in log scale)
        log_prop<-lpostf(curra1,curra2,curra3,propa4,curra5,curra6,curra7,curra8,data,currcx,currcy,currd)
        logMHratio4=log_prop-log_curr+log(1-pnorm(-curra4/s))-log(1-pnorm(-propa4/s))
                     
        logalpha4 = min(0,logMHratio4) # alpha = min(1,MHratio)
        if ((log(runif(1))<(log(1)+logalpha4))) # accept if unif(0,1)<alpha, i.e. accept with probability alpha, else stay at current state
          {
            acc_a4 = acc_a4 + 1 # increment count of accepted proposals
            curra4 = propa4
            log_curr<-log_prop
           }
        
       ## sample from full conditional distribution of theta (Metropolis-Hastings update)
      propa5 = rtnorm(1,mean=curra5,sd=s1,0,Inf) # draw one sample at random)

        ## Metropolis accept-reject step (in log scale)
        log_prop<- lpostf(curra1,curra2,curra3,curra4,propa5,curra6,curra7,curra8,data,currcx,currcy,currd)
        logMHratio5 = log_prop - log_curr+log(1-pnorm(-curra5/s1))-log(1-pnorm(-propa5/s1))
                             
        logalpha5 = min(0,logMHratio5) # alpha = min(1,MHratio)
        if ((log(runif(1))<(log(1)+logalpha5))) # accept if unif(0,1)<alpha, i.e. accept with probability alpha, else stay at current state
          {
            acc_a5 = acc_a5 + 1 # increment count of accepted proposals
            curra5 = propa5
            log_curr<-log_prop
           }
        
         ## sample from full conditional distribution of theta (Metropolis-Hastings update)
      propa6 = rtnorm(1,mean=curra6,sd=s1,0,Inf) # draw one sample at random from exp )

        ## Metropolis accept-reject step (in log scale)
        log_prop<- lpostf(curra1,curra2,curra3,curra4,curra5,propa6,curra7,curra8,data,currcx,currcy,currd)
        logMHratio6 = log_prop - log_curr+log(1-pnorm(-curra6/s1))-log(1-pnorm(-propa6/s1))
                     
        logalpha6 = min(0,logMHratio6) # alpha = min(1,MHratio)
        if ((log(runif(1))<(log(1)+logalpha6))) # accept if unif(0,1)<alpha, i.e. accept with probability alpha, else stay at current state
          {
            acc_a6 = acc_a6 + 1 # increment count of accepted proposals
            curra6 = propa6
            log_curr<-log_prop
           }

        ## sample from full conditional distribution of theta (Metropolis-Hastings update)
       propa7 = rtnorm(1,mean=curra7,sd=s1,0,Inf) # draw one sample at random from exp 

        ## Metropolis accept-reject step (in log scale)
        log_prop<- lpostf(curra1,curra2,curra3,curra4,curra5,curra6,propa7,curra8,data,currcx,currcy,currd)
        logMHratio7 = log_prop - log_curr+log(1-pnorm(-curra7/s1))-log(1-pnorm(-propa7/s1))
                            
        logalpha7 = min(0,logMHratio7) # alpha = min(1,MHratio)
        if ((log(runif(1))<(log(1)+logalpha7))) # accept if unif(0,1)<alpha, i.e. accept with probability alpha, else stay at current state
          {
            acc_a7 = acc_a7 + 1 # increment count of accepted proposals
            curra7 = propa7
            log_curr<-log_prop
           }
       
       ## sample from full conditional distribution of theta (Metropolis-Hastings update)
       propa8 = rtnorm(1,mean=curra8,sd=s1,0,Inf) # draw one sample at random from exp )

        ## Metropolis accept-reject step (in log scale)
        log_prop<- lpostf(curra1,curra2,curra3,curra4,curra5,curra6,curra7,propa8,data,currcx,currcy,currd)
        logMHratio8 = log_prop - log_curr+log(1-pnorm(-curra8/s1))-log(1-pnorm(-propa8/s1))
                           
        logalpha8 = min(0,logMHratio8) # alpha = min(1,MHratio)
        if ((log(runif(1))<(log(1)+logalpha8))) # accept if unif(0,1)<alpha, i.e. accept with probability alpha, else stay at current state
          {
            acc_a8 = acc_a8 + 1 # increment count of accepted proposals
            curra8 = propa8
            log_curr<-log_prop
           }
       
       ## sample from full conditional distribution of theta (Metropolis-Hastings update)
       propcx = rtnorm(1,mean=currcx,sd=s3,0,Inf) # draw one sample at random from exp )

        ## Metropolis accept-reject step (in log scale)
        log_prop<- lpostf(curra1,curra2,curra3,curra4,curra5,curra6,curra7,curra8,data,propcx,currcy,currd)
        logMHratio9 = log_prop - log_curr+log(1-pnorm(-currcx/s3))-log(1-pnorm(-propcx/s3))
                     
       logalpha9 = min(0,logMHratio9) # alpha = min(1,MHratio)
        if ((log(runif(1))<(log(1)+logalpha9))) # accept if unif(0,1)<alpha, i.e. accept with probability alpha, else stay at current state
          {
            acc_cx = acc_cx + 1 # increment count of accepted proposals
            currcx = propcx
            log_curr<-log_prop
           }
      
      ## sample from full conditional distribution of theta (Metropolis-Hastings update)
       propcy = rtnorm(1,mean=currcy,sd=s3,0,Inf) # draw one sample at random from exp )

        ## Metropolis accept-reject step (in log scale)
        log_prop<- lpostf(curra1,curra2,curra3,curra4,curra5,curra6,curra7,curra8,data,currcx,propcy,currd)
        logMHratio10 = log_prop - log_curr+log(1-pnorm(-currcy/s3))-log(1-pnorm(-propcy/s3))
                     
        logalpha10 = min(0,logMHratio10) # alpha = min(1,MHratio)
        if ((log(runif(1))<(log(1)+logalpha10))) # accept if unif(0,1)<alpha, i.e. accept with probability alpha, else stay at current state
          {
            acc_cy = acc_cy + 1 # increment count of accepted proposals
            currcy = propcy
            log_curr<-log_prop
           }
       
        cat("\n")

          ## update chain with new values
        
     mchain[,i] = c(curra1,curra2,curra3,curra4,curra5,curra6,curra7,curra8,currcx,currcy,currd,log_curr,rho1d,rho2d,rho12)
   batch=c(seq(3000,25000, length=11)) 
   batch_index=match(i,batch,nomatch=0)
   if (batch_index > 0){ 
   rho1d=cor(mchain[1,1:batch[batch_index]],mchain[11,1:batch[batch_index]])     ### update the correltion matrix of the parameters a11, a12 and Dmax at regular interval
   rho2d=cor(mchain[2,1:batch[batch_index]],mchain[11,1:batch[batch_index]])
   rho12=cor(mchain[1,1:batch[batch_index]],mchain[2,1:batch[batch_index]])
   sa1=sqrt(var(mchain[1,1:batch[batch_index]]))
   sa2=sqrt(var(mchain[2,1:batch[batch_index]]))
   s2=sqrt(var(mchain[11,1:batch[batch_index]]))
   mua1=mean(mchain[1,1:batch[batch_index]])
   mua2=mean(mchain[2,1:batch[batch_index]])
   mud=mean(mchain[11,1:batch[batch_index]])

   mu12d=c(mua1,mua2,mud)
   Sigma12d=matrix(c(sa1^2,rho12*sa1*sa2,rho1d*sa1*s2,rho12*sa1*sa2,sa2^2,rho2d*sa2*s2,rho1d*sa1*s2,rho2d*sa2*s2,s2^2),3,3)
      }
     }

   cat("Markov chain algorithm ran for ",NUMIT,"\n", "acc.rate for a12d=",acc_a12d/(NUMIT-1),"\n","acc.rate for a3=",acc_a3/(NUMIT-1),"\n",
       "acc.rate for a4=",acc_a4/(NUMIT-1),"\n","acc.rate for a5=",acc_a5/(NUMIT-1),"\n","acc.rate for a6=",acc_a6/(NUMIT-1),"\n","acc.rate for a7=",acc_a7/(NUMIT-1),"\n",
        "acc.rate for cx=",acc_cx/(NUMIT-1),"\n","acc.rate for cy=",acc_cy/(NUMIT-1),"\n")
    return(mchain)  
          }

/*Function call for MCMC updates*/


a1_hat1<-MHsampler(25000,s=0.002,s1=0.001,s2=0.27,s3=0.003)


/* Estimation of autocorrelation parameter(s)

K: size of the sub-regions in the image
residual:  residual image obtained after removing the trend

*/

bb=seq(1,1024,K)
a_entry=bb[-length(bb)]+1
aa=c(1,a_entry)
rhob=matrix(0,length(aa),length(bb))
for (i in 1:length(aa)){
     for (j in 1:length(bb)){
           block=as.vector(residual[aa[i]:bb[i],aa[j]:bb[j]])
           if (any(is.na(block))==FALSE)
           {
           rhob[i,j]=sp.correlogram(nbd,block,order=5,method="corr",style="W",zero.policy=T)$res[1]
            }
           else
           {
           rhob[i,j]=NA                          
             }
  }
}
/***********************************************/

/*    Spatial analysis of simulated image data */

/***********************************************/

/*define function to remove NA's */
new_sum=function(x) {
                      l=sum(x) 
                       is.na(l)
                    }

/*define nbd of each pixel: excluding the edge pixels*/

x2=2:99
g1=expand.grid(x=x2,y=x2)
nbd=matrix(c(g1$x+100*(g1$y-1),ifelse(g1$y>1,100*(g1$y-2)+g1$x,0),ifelse(g1$x>1,100*(g1$y-1)+g1$x-1,0),ifelse(g1$x+1<=100,g1$x+1+100*(g1$y-1),0),ifelse(g1$y+1<=100,g1$x+100*g1$y,0)),9604,5)


### Likelihood function for spatial dependent case ###

lpostf<-function(a1,a2,a3,a4,a5,a6,a7,a8,y=data,cxx,cyy,d,nd=nbd){
              cx=floor(cxx)
              cy=floor(cyy)
              del_cx=cxx-cx
              del_cy=cyy-cy
                           
              I<-rep(0,10000)
              res<-rep(0,10000)              
              S<-rep(0,10000)
              Kfunc<-rep (0,10000)
              x1<-seq(0.001,sqrt(2)*0.075,length=100)
              y1<-seq(0.001,sqrt(2)*0.075,length=100)
              
              g=expand.grid(x=x1,y=y1)
              S<-sqrt((g$x-x1[cx]+0.000001-del_cx*0.001061273)^2+(g$y-y1[cy]+0.000001-del_cy*0.001061273)^2)

              a<-c(a1,a2,a3,a4,a5,a6,a7,a8)/100000
              Kfunc<-(a[1]*b(S,1,d)+a[2]*b(S,2,d)+a[3]*b(S,3,d)+a[4]*b(S,4,d)+a[5]*b(S,5,d)+a[6]*b(S,6,d)+a[7]*b(S,7,d)+a[8]*b(S,8,d))
              U<-ifelse((Kfunc > 0)==T,Kfunc,0.00001)
              I<-(U*100000)/S
              res<-(y-I)
              Vinv=1/as.vector(I)
             
              suml=0
                                    
             Res1=res[nd[,1]]
             Res2=res[nd[,2]]
             Res3=res[nd[,3]]
             Res4=res[nd[,4]]
             Res5=res[nd[,5]]
            
             Var1=Vinv[nd[,1]]
             Var2=Vinv[nd[,2]]
             Var3=Vinv[nd[,3]]
             Var4=Vinv[nd[,4]]
             Var5=Vinv[nd[,5]]
 
             A1=(Var1+0.25*(Var2+Var3+Var4+Var5))*(Res1^2)
             A2=(Var2+0.015625*Var1)*(Res2^2)
             A3=(Var3+0.015625*Var1)*(Res3^2)
             A4=(Var4+0.015625*Var1)*(Res4^2)
             A5=(Var5+0.015625*Var1)*(Res5^2)

             A12=2*(0.5*Var2+0.125*Var1)*(Res2*Res1)
             A13=2*(0.5*Var3+0.125*Var1)*(Res3*Res1)
             A14=2*(0.5*Var4+0.125*Var1)*(Res4*Res1)
             A15=2*(0.5*Var5+0.125*Var1)*(Res5*Res1)

             A23=2*(0.015625*Var1)*(Res2*Res3)
             A24=2*(0.015625*Var1)*(Res2*Res4)
             A25=2*(0.015625*Var1)*(Res2*Res5)

             A34=2*(0.015625*Var1)*(Res3*Res4)
             A35=2*(0.015625*Var1)*(Res3*Res5)
             A45=2*(0.015625*Var1)*(Res4*Res5)

             detV=log(Var1)+log(Var2)+log(Var3)+log(Var4)+log(Var5)
             
             
             llike=0.125*((A1+A2+A3+A4+A5+A12+A13+A14+A15+A23+A24+A25+A34+A35+A45)-0.5*detV)
            
              suml=sum(llike,na.rm=T)
             
              prior_part<-sum(log(a))+sum((0.22*(log(a[1:8])-1)^2))+2*d-109*log(d)+(((cxx-50)^2)/8)+(((cyy-50)^2)/8)
              
               logpost<-(-(prior_part+suml)) 
              return(logpost)
             }


/********************************************/

/*     P(r) from estimated coeffcients      */

/********************************************/

/* trapezoidal function: to calculate area under curve */

trap.rule <- function(x,y) {
                  sum(diff(x)*(y[-1]+y[-length(y)]))/2
                           }


/* Plot reference P(r): for simulated data */

/*
a<-c( 5.281465e-05, 4.178242e-05, 1.686120e-05 ,5.454271e-06, 3.379092e-06, 3.007528e-06, 2.055125e-06, 8.427587e-07)
d<-58.5       #Dmax
*/

r<-seq(0,d,length=500)   # evaluate P(r) on 500 discrete r values between (0,d)   
Pr_ref=8*pi*r*(a[1]*sin((pi*r*1)/d)+a[2]*sin((pi*r*2)/d)+a[3]*sin((pi*r*3)/d)+a[4]*sin((pi*r*4)/d)+a[5]*sin((pi*r*5)/d)+a[6]*sin((pi*r*6)/d)+a[7]*sin((pi*r*7)/d)+a[8]*sin((pi*r*8)/d))
Pr_ref_const=trap.rule(r,Pr_ref)  # normalizing constant
plot((Pr_ref/Pr_ref_const)~r, xlab="r", ylab="P(r)",col="gray",type="l",lwd=1,ylim=c(0,0.04))


/* Estimated P(r)'s  using esimated parameters from MCMC 

Numcurves: 10000  # of sampled posterior P(r)'s 
a1_hat1: MCMC samples of parameters

*/

Pr_est<-matrix(NA,nrow=500,ncol=Numcurves)
R<-matrix(NA,nrow=500,ncol=Numcurves)
Pr_std<-matrix(NA,nrow=500,ncol=Numcurves)  ##standardized P(r) such that the area under the curve sum to 1

for (k in 1:Numcurves){
    R[,k]=seq(0,a1_hat[11,k+15000],length=500)
    Pr_est[,k]=8*pi*R[,k]*(a1_hat1[1,k+15000]*sin((pi*R[,k]*1)/R[500,k])+a1_hat1[2,k+15000]*sin((pi*R[,k]*2)/R[500,k])+a1_hat1[3,k+15000]*sin((pi*R[,k]*3)/R[500,k])+a1_hat1[4,k+15000]*sin((pi*R[,k]*4)/R[500,k])+a1_hat1[5,k+15000]*sin((pi*R[,k]*5)/R[500,k])+a1_hat1[6,k+15000]*sin((pi*R[,k]*6)/R[500,k])+a1_hat1[7,k+15000]*sin((pi*R[,k]*7)/R[500,k])+a1_hat1[8,k+15000]*sin((pi*R[,k]*8)/R[500,k]))
    Pr_std[,k]=Pr_est[,k]/trap.rule(R[,k],Pr_est[,k])
      }



/*
 Integrated error with respect to reference P(r) for both 2-D and 1-D method

x=matrix of estimated parameter values from MCMC updates
a=coefficients for reference P(r)
am=coefficients from 1-D method

*/

integrated_error=function(x,a,am){
d=a[length(a)]
ll=dim(x)[2]
dm=am[length(am)] 

r=seq(0,d,length=500)
Porig=8*pi*r*(a[1]*sin((pi*r*1)/d)+a[2]*sin((pi*r*2)/d)+a[3]*sin((pi*r*3)/d)+a[4]*sin((pi*r*4)/d)+a[5]*sin((pi*r*5)/d)+a[6]*sin((pi*r*6)/d))
const_mcmc=trap.rule(r,Porig)
#cat("c=",const_mcmc)
rm=seq(0,dm,length=500)
Pmoore=8*pi*rm*(am[1]*sin((pi*rm*1)/dm)+am[2]*sin((pi*rm*2)/dm)+am[3]*sin((pi*rm*3)/dm)+am[4]*sin((pi*rm*4)/dm)+am[5]*sin((pi*rm*5)/dm)+am[6]*sin((pi*rm*6)/dm))    
const_moore=trap.rule(rm,Pmoore)
       
    for (k in 1:ll){
                    R[,k]=seq(0,x[9,k+n],length=500)
                    PR[,k]=8*pi*R[,k]*(x[1,k+n]*sin((pi*R[,k]*1)/R[500,k])+x[2,k+n]*sin((pi*R[,k]*2)/R[500,k])+x[3,k+n]*sin((pi*R[,k]*3)/R[500,k])+x[4,k+n]*sin((pi*R[,k]*4)/R[500,k])+x[5,k+n]*sin((pi*R[,k]*5)/R[500,k])+x[6,k+n]*sin((pi*R[,k]*6)/R[500,k]))
                    num[k]=trap.rule(R[,k],PR[,k]) 
                    abs_diff[k]=sum(abs((PR[,k]/num[k])-(Porig/const_mcmc)))
                    sum_moore=sum(abs((Pmoore/const_moore)-(Porig/const_mcmc)))
                          }
     return(list(abs_diff=abs_diff, sum_moore=sum_moore))
                         }

/* call integrated_error */

R=matrix(0,500,dim(a1_hat1)[2])   #a1_hat1: matrix of estimated parameter values from 2-D method
PR=matrix(0,500,dim(a1_hat1)[2])
PR1=matrix(0,500,dim(a1_hat1)[2])
abs_diff=rep(0,dim(a1_hat1)[2])
num=rep(0,dim(a1_hat1)[2])
amoore: #coefficients from 1D method   

error_2D=integrated_error(a1_hat1,a,amoore)$abs_diff
error_1D=integrated_error(a1_hat1,a,amoore)$sum_moore



/*
Generate crossvalidation dataset (mask) for experimental image
*/

num=133 # of cells sampled in each x- and y- direction

x1=sort(sample(1:1024,num))
y1=sort(sample(1:1024,num))
l1=numeric(1024)
l2=numeric(1024) 

l1[x1]=1
a1=l1
l2[y1]=1
a2=l2

g=expand.grid(x=a1,y=a2)
sumg=g$x+g$y
tt=ifelse(sumg >1,NA,0)     ##keep all cells except the sampled cells for test 
test_mask=matrix(tt,1024,1024)    
vv=ifelse(sumg >1,0,NA)     ##keep only the sampled cells for validation
valid_mask=matrix(vv,1024,1024)




valid_mask=matrix(vv,1024,1024)