	program mimain
	implicit double precision (a-h,o-z)
	integer size,dim,grid1,grid2
	double precision noiselevel,zero
	include 'param.dat'
	parameter (nwcore = 5000000)
	double precision z(nwcore)
	
	double precision genXY

C       define the matrix
	double precision myX(size,dim), ftrue(size), Ker(size,size,dim)
C       y is double for Gaussian case; int for Bernoulli case
	double precision y(size)
	double precision theta(dim), Kth(size,size)
	double precision eye(size+1,size+1) 
	double precision lambda0, lambda1
C       arguments for tuning lambda0
	double precision KI(size,size), cbK(size+1,size+1), cby(size+1)
	double precision fest(size), sqrtmse(grid1)
	integer bestk
C       arguments for tuning lambda1
	common /c1/ myX
	common /c2/ ftrue
	common /c3/ Ker
C       arguments for LAPACK
	integer ipiv(size+1), info
C       tmp arguments
	double precision tmp


C       generate original myX, y, Ker, ftrue
	call simuXY(y)
	
C       initialize theta
	do 10, k = 1, dim
	   theta(k) = 1.0d+0
 10	continue
	
C       contruct BigKernel Kth: positive definite
	open(1,file = 'Kth', status = 'old',access = 'sequential', form = 'formatted')
	do 40, j = 1,size
	   do 30, i = 1,size
	      Kth(i,j) = zero
	      do 20, k = 1,dim
		 Kth(i,j) = Kth(i,j) + theta(k)*Ker(i,j,k)
 20	      continue
 30	   continue
	   write(1,1000) (Kth(i,j),i = 1,size)
 40	continue
	close(1)
	
	
C       generate eye
	do 60, j = 1, size+1
	   do 50, i = 1, size+1
	      if (i.ne.j) then
		 eye(i,j) = zero
	      else 
		 eye(i,j) = 1.0d+0
 	      endif
 50	   continue
 60	continue
	
	
C       generate cby, cbK and solve cbK*cb = cby
C       NOTE: cbK is NOT positive definite, thus cholesky is not proper
        do 160, k = 1,grid1
	   sqrtmse(k) = zero
	   write(*,*) "y"
C       Redefine cby each time, since it changed after LAPACK call
           do 70, i = 1,size+1
              if (i.le.size) then
		 cby(i) = y(i)
              else
		 cby(i) = zero
	      endif
 70	   continue
           write(*,*) 'cby', (cby(i),i=1,size+1)
           lambda0 = (2.0d+0)**(-k)
C       Define KI
	   do 90, j = 1, size
	      do 80, i = 1, size
		 if (i.eq.j) then
		    KI(i,j) = Kth(i,j) + lambda0*eye(i,j)
		 else 
		    KI(i,j) = Kth(i,j)
		 endif
 80		 continue
C               write(*,*) 'KI',j,(KI(i,j), i = 1,size)
 90	      continue
C       Define cbK
	    do 110, j = 1, size+1
	       do 100, i = 1, size+1
		  if (i.le.size .and. j.le.size) then
		     cbK(i,j) = KI(i,j)
		  else if (i.eq.size+1 .and. j.le.size) then
		     cbK(i,j) = 1.0d+0
		  else if (i.le.size .and. j.eq.size+1) then
		     cbK(i,j) = 1.0d+0
		  else
		     cbK(i,j) = zero
		  endif
 100	       continue
 110	    continue
C       ----- call LAPACK: 
C       call SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
	    call dgetrf(size+1,size+1,cbK,size+1,ipiv,info)
C       call SUBROUTINE DGETRS(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
C       after this call, cby changed
	    call dgetrs('N',size+1,1,cbK,size+1,ipiv,cby,size+1,info)
C       print out the solution cby
	    write(*,*) 'c',(cby(i),i =1,size)
	    write(*,*) 'b',cby(size+1)
C       compute fest
	    do 140, i = 1, size
	       fest(i) = zero
	       do 130, j = 1, size
		  fest(i) = fest(i) + Kth(j,i)*cby(j)
 130	       continue
	       fest(i) = fest(i) + cby(size+1)
 140	    continue
	    do 150, i = 1,size
	       sqrtmse(k) = sqrtmse(k) + (fest(i)-ftrue(i))**2
 150	    continue
	    sqrtmse(k) = sqrt(sqrtmse(k))
 160	 continue
	 
	
C	choose optimal lambda0
	 bestk = 1
	 tmp = sqrtmse(1)
	 do 170, k = 1, grid1
	    write(*,*) 'sqrtmse',k,sqrtmse(k)
	    if (k.lt.grid1 .and.tmp .gt. sqrtmse(k+1)) then
	       bestk = k+1
	       tmp = sqrtmse(k+1)
	    end if
 170	continue
	write(*,*) "bestk", bestk
	lambda0 = 2**(-bestk)
	
	
 1000	format(6(2X,D14.7))
	end
	








