! Basic statistical analysis of "test snapshots" collected from different hosts
! with the BASH script produce_test.sh. Those data are stored in the files prec_host.o2.1 etc,
! where prec=32,64,80 and the last digit labels the run number. 
! We calculate average execution time by averaging over different runs and also evaluate dispersion
! and write results in prec_host.o2  
    program stata
    implicit none
    integer:: np(3),nr(3),nc
    real(16), parameter:: pi26=3.1415926535897932384626434_16 
    real(16), parameter:: e26= 2.7182818284590452353602875_16
    real(16):: x,y(0:1)
    real(16):: t(0:1,10,50),sigmaT(0:1,50),t_av(0:1,50),t2_av(0:1,50)
    real(16):: v(0:1,10,50),sigmaV(0:1,50),v_av(0:1,50),v2_av(0:1,50)
    real(8):: cnt
    integer(8):: two2p(50)
    integer:: p(50),ks,i,j,k,n,l,m
    character(128) :: h(10),dir,fn(10),fn1,cline
    character(4):: prec(3),tst(3)
    character(20):: ext,str

! Initialization
! no. of runs for different tests
    nr(1)=4 !8 !SUM
    nr(2)=4 !8 !E
    nr(3)=4 !8 !PI

! no. of lines in data files corresponding to different tests
!    np(1)=25
    np(1)=26
    np(2)=35 !36
    np(3)=35
    
    prec=(/"32_","64_","80_"/)
	
    tst(1)="_sum"
    tst(2)="_e"
    tst(3)="_pi"
! host names
      dir="./data2/"
      h(1)="cluster"
      h(2)="cluster13"
      h(3)="cluster31"
      h(4)="cluster42"
      h(5)="cluster51"
      h(6)="cluster52"
      h(7)="rut"
	  h(8)="nuke"
!!!!!! Make sure files prec_host_tst.o2.runnumber are present in dir
	do k=5,5 !hosts
	  do j=1,3 !precision 32,64 or 80
	    do i=1,3 !test type
	      nc=nr(i) 
	      do n=1,nc !average over runs (prec_prl_host.o2.runnumber)
! prepare the name fn of the data file
	         write(str,'(i2)') n
	         ext=".o2"
!	         ext=".o2native"
  fn(n)=trim(dir)//trim(prec(j))//trim(h(k))//trim(tst(i))//trim(ext)//"."//trim(adjustl(str))
	         print*, trim(fn(n))
! read in data from  fn(n)
             v(:,n,:)=0._16
             t(:,n,:)=0._16	         
	           open(1,file=trim(fn(n)),status='old')
	             read(1,'(a)') cline !read  non-numeric line
	               if(i>1) read(1,'(a)') cline !read  non-numeric line
!	               if(i>1) read(1,'(a)') cline !read  non-numeric line
	               do l=1,np(i) !read the rest of the file
	                 read(1,*,end=11) p(l),two2p(l),ks,v(0,n,l),t(0,n,l),v(1,n,l),t(1,n,l)
	               enddo ! end read from file
11                   continue
	           close(1)
	      enddo !n
	      
!!!!!! basic stat analysis: calculate mean values and times averaged over the runs and dispersion for each data point
	      t_av(:,:)=0._16
	      t2_av(:,:)=0._16
	      v_av(:,:)=0._16
	      v2_av(:,:)=0._16
          do m=0,1
	        do l=1,np(i)
	          cnt=0. !reset counter
	          do n=1,nc
	            v_av(m,l)=v_av(m,l) + v(m,n,l)
	            v2_av(m,l)=v2_av(m,l) + v(m,n,l)**2
	              x=t(m,n,l)
	              if(x > 1d-8) then 
	                t_av(m,l)=t_av(m,l)  + x
	                t2_av(m,l)=t2_av(m,l)+ x**2
	                cnt=cnt+1.
	              else !skip time=0 or <0 from calculation of average time 
	                print*,"***Skipping in "//trim(fn(n)),m, two2p(l), v(m,n,l), real(x)
!	                write(*,*) p(l),two2p(l),ks,v(m,n,l),x
	              end if
	          end do !n
! average time and dispersion
                  if(cnt > 1d-8) then
	            t_av(m,l)=t_av(m,l) /cnt !! average time
	            t2_av(m,l)=t2_av(m,l) /cnt
	          endif 
	            x=t2_av(m,l) - t_av(m,l)**2
	            sigmaT(m,l)=sqrt(x) !! dispersion, time
	            v_av(m,l)=v_av(m,l) /real(nc,16) ! average value
	            v2_av(m,l)=v2_av(m,l) /real(nc,16) 
	            x=v2_av(m,l) - v_av(m,l)**2
	            sigmaV(m,l)=sqrt(x) !! dispersion, value
	        end do !l
	      end do !m
20	format(i3,i15,i4,2e15.6)

!!!!!!!!!SUM 0 test
	  if(i==1) then 
	  
! write a new file with average values and dispersions
	      fn1=trim(dir)//trim(prec(j))//trim(h(k))//"_SUM"//trim(ext)
	      write(*,*) ">>>> "//trim(fn1)
	    open(2,file=trim(fn1))
        write(2,21)'# p','N=2**p','S(seq)','S_err','T(seq)','T_err','S(par)','S_err','T(par)','T_err'
!                     1   2        3        4       5        6       7        8       9        10    
	    do l=1,np(i)
        write(2,23) p(l),two2p(l),v_av(0,l),sigmaV(0,l),t_av(0,l),sigmaT(0,l),&
                                  v_av(1,l),sigmaV(1,l),t_av(1,l),sigmaT(1,l)
	    end do
	    close(2)
21  format(a3,a15,2(a18,3a15))
23  format(i3,i15,2(e18.8,3e15.6))

!!!!!!!!!!!!!!!!E test
	  elseif(i==2) then
	  
! write a new file with average values and dispersions
	      fn1=trim(dir)//trim(prec(j))//trim(h(k))//"_E"//trim(ext)
	      write(*,*) ">>>> "//trim(fn1)
	    open(2,file=trim(fn1))
!        write(2,31)'# p','N=2**p','E(seq)','E_err','T(seq)','T_err','E(par)','E_err','T(par)','T_err'
        write(2,31)'# p','N=2**p','E(seq)','(E-E26)/E26','T(seq)','T_err','E(par)','(E-E26)/E26','T(par)','T_err'
!                     1   2        3        4             5        6       7        8              9        10    
	    do l=1,np(i)
	      y(0)=(v_av(0,l)-e26)/e26
	      y(1)=(v_av(1,l)-e26)/e26
          write(2,33) p(l),two2p(l),v_av(0,l),y(0),t_av(0,l),sigmaT(0,l),&
                                    v_av(1,l),y(1),t_av(1,l),sigmaT(1,l)
	    end do
	    close(2)
	   
!!!!!!!!!!!!!!!!!PI test
	  elseif(i==3) then 

! write a new file with average values and dispersions
	      fn1=trim(dir)//trim(prec(j))//trim(h(k))//"_PI"//trim(ext)
	      write(*,*) ">>>> "//trim(fn1)
	    open(2,file=trim(fn1))
        write(2,31)'# p','N=2**p','PI(seq)','(PI-PI26)/PI26','T(seq)','T_err','E(par)','(PI-PI26)/PI26','T(par)','T_err'
!                     1   2        3         4                5        6       7         8               9        10    
	    do l=1,np(i)
	      y(0)=(v_av(0,l)-pi26)/pi26
	      y(1)=(v_av(1,l)-pi26)/pi26
          write(2,33) p(l),two2p(l),v_av(0,l),y(0),t_av(0,l),sigmaT(0,l),&
                                    v_av(1,l),y(1),t_av(1,l),sigmaT(1,l)
	    end do
	    close(2)
	  
	  endif    
!!!!!!!!!!!!!!!!!!!!!!	      
            enddo !i test
          enddo !j: precision 32,64,80
        enddo !k: host 

31  format(a3,a15,2(a26,3a15))
33  format(i3,i15,2(f26.21,3e15.6))
	end program