! Calculate, for each dataset in nXDS_ASCII.HKL file, its CC against all others. Write output file.
! Version 0.0 Kay Diederichs 12.2.2018
!
! ifort -check all -traceback -warn -init=snan,arrays -fpe0 -static-intel -sox resolution.o symops.o cc_mod.o nxds_getcc.f90 -o nxds_getcc
! after previous ifort -c symops.f90 resolution.f90 cc_mod.f90
!
! gfortran -O2 -fsanitize=undefined -static-libgfortran -static-libubsan symops.f90 resolution.f90 cc_mod.f90 nxds_getcc.f90 -o nxds_getcc.gfortran
! 
!
! 29.1.2020 drop requirement that data set numbers are sorted

MODULE nxdsitems_mod
      INTEGER nset
      REAL resmin,resmax
      INTEGER spgnum
contains
      SUBROUTINE get_nset_reso(name)
      USE resolution_mod, ONLY: getres,a,b,c,alf,bet,gam
      IMPLICIT NONE
      CHARACTER*(*), INTENT(IN) :: name
      CHARACTER line*80
      INTEGER :: ier,unit,iset,ih,ik,il,item_img=0
      REAL :: s2,resol,dum(4),fsq,sigma
      REAL, ALLOCATABLE :: resminset(:),resmaxset(:)
 1000 FORMAT(a)
      
      OPEN(NEWUNIT=unit,FILE=name,STATUS='old',ACTION='READ',iostat=ier)
      IF (ier/=0) THEN
        PRINT*,'ier=',ier
        STOP 'could not open file'
      END IF
! Open input file; find number of datasets, and allocate array inputfile
      DO
        READ(unit,1000,iostat=ier) line
        IF (ier/=0) THEN
          PRINT*,'ier=',ier,' line=',TRIM(line)
          STOP 'read error'
        END IF
        IF (line(1:18)=='!NUMBER_OF_IMAGES=') then
          READ(line(19:),*) nset
          WRITE(*,'(a)') TRIM(line)
        ELSE IF (line(1:20).eq.'!SPACE_GROUP_NUMBER=') then 
          READ(line(22:),*) spgnum
          WRITE(*,'(a)') TRIM(line)
        ELSE IF (line(1:21).eq.'!UNIT_CELL_CONSTANTS=') then 
          READ(line(22:),*) a,b,c,alf,bet,gam
          WRITE(*,'(a)') TRIM(line)
        ELSE IF (line(1:10).eq.'!ITEM_IMG=') then 
          READ(line(11:),*) item_img
          WRITE(*,'(a)') TRIM(line)
        ELSE IF (line(1:14) == '!END_OF_HEADER') THEN
          EXIT
        END IF
      END DO
      IF (item_img==0) STOP 'ITEM_IMG not found. Use MERGE=FALSE in nXDS.INP!'
      IF (item_img/=10) STOP 'ITEM_IMG=10 expected'
      ALLOCATE(resminset(nset),resmaxset(nset))
      resmin=1000.
      resmax=-1.
      resminset=1000.
      resmaxset=-1.
! find min and max resolution
      DO
        READ(unit,1000,iostat=ier) line
        IF (ier/=0) THEN
          PRINT*,'ier=',ier
          PRINT*,'line=',TRIM(line)
          STOP 'read error'
        END IF
        IF (line(1:12) == '!END_OF_DATA') EXIT
        READ(line,*)ih,ik,il,fsq,sigma,dum,iset
        IF (sigma<0.) CYCLE
        CALL getres(ih,ik,il,s2,resol)
        IF (resol<resmin) resmin=resol
        IF (resol>resmax) resmax=resol      
        IF (resol<resminset(iset)) resminset(iset)=resol
        IF (resol>resmaxset(iset)) resmaxset(iset)=resol
      END DO     
      WRITE(*,'(a,2f10.3)')'highest and lowest resol of any image:',resmin,resmax
      resmin=MAXVAL(resminset,MASK=resminset<1000.)
      resmax=MINVAL(resmaxset,MASK=resmaxset>-1.) 
      WRITE(*,'(a,2f10.3)')'highest and lowest COMMON resol of all images:',resmin,resmax
      ier=MAXLOC(resminset,MASK=resminset<1000.,DIM=1)
      WRITE(*,'(a,i6,1x,f10.3)') 'index of image with lowest high-resol:',&
                                ier,resminset(ier)
      ier=MINLOC(resmaxset,MASK=resmaxset>-1.,DIM=1)
      WRITE(*,'(a,i6,1x,f10.3)') 'index of image with highest low-resol:',&
                                ier,resmaxset(ier)
      WRITE(*,'(a,i0)') 'count of images without reflections: ',COUNT(resmaxset==-1.)
      CLOSE(unit)
      END SUBROUTINE get_nset_reso

SUBROUTINE nxds_file(name,reslim1,reslim2,nbin,weighted,nset,ccarray,nccarray)
      USE symops_mod, ONLY: getops,asuput,op,ncenop
      USE resolution_mod, ONLY: getres,getbin,getlim
      USE cc_mod, ONLY: cc_type, get_cc,put_cc
      IMPLICIT NONE
      CHARACTER*(*), INTENT(IN) :: name
      INTEGER, INTENT(IN) :: nbin
      REAL, INTENT(IN) :: reslim1,reslim2
      LOGICAL, INTENT(IN) :: weighted
      INTEGER, INTENT(IN) :: nset
      REAL, INTENT(OUT) :: ccarray(nset,nset)
      INTEGER, INTENT(OUT) :: nccarray(nset,nset)
      INTEGER :: ih,ik,il,oldhkl=0,jh,jk,jl,ihkl=0,iop,i,j,k,n,ispgr,ncheck,ier, &
                 iset,jset,nref=0,kset,nobs=0,ibin,misfits=0,ntot(nbin),ncontrib=0, &
                 modprint=1000,unit,narray(40)=0,kset1,jset1
!      INTEGER, PARAMETER :: i8kind=SELECTED_INT_KIND(15)  
!      INTEGER(i8kind) :: ipointer,newsize
      INTEGER(8) :: ipointer,newsize
      REAL :: vvr(nbin),s2,resol,fsq,sigma,dum(4),s2save,xx4,yy4
      CHARACTER :: line*200
      LOGICAL :: last=.FALSE.
      INTEGER, ALLOCATABLE :: na(:),icontrib(:),ifound(:,:)
      REAL, ALLOCATABLE :: xa(:),siga(:)
      TYPE(cc_type), ALLOCATABLE :: ccbin(:,:),temp_ccbin(:,:)

 1000 FORMAT(A)
 
      ntot=0
      ncheck=0
! Open input file; get spacegroup and cell constants
      OPEN(NEWUNIT=unit,FILE=name,STATUS='old',ACTION='READ',iostat=ier)
      IF (ier/=0) THEN
        PRINT*,'ier=',ier
        STOP 'could not open file'
      END IF
      DO
        READ(unit,1000,iostat=ier) line
        IF (ier/=0) THEN
          PRINT*,'ier=',ier
          PRINT*,'line=',TRIM(line)
          STOP 'read error'
        END IF
        IF (line(1:20)=='!SPACE_GROUP_NUMBER=') then
          READ(line(21:),*) ispgr
          CALL getops(ispgr)
        ELSE IF (line(1:9) == '!ITEM_H=1') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:9) == '!ITEM_K=2') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:9) == '!ITEM_L=3') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:12) == '!ITEM_IOBS=4') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:19) == '!ITEM_SIGMA(IOBS)=5') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:12) == '!ITEM_IMG=10') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:14) == '!END_OF_HEADER') THEN
          EXIT
        END IF
      END DO
      IF (ncheck <6) STOP 'not all items H,K,L,IOBS,SIGMA(IOBS),IMG found'  
! ccbin needs most memory - 44 bytes per entry if REAL(8), 24 if REAL(4)
      WRITE(*,'(a,i0,a)') 'memory allocated is at most nbin*nset*nset= ',NINT((24*(nbin/1.e6)*nset)*nset),' MB'
      newsize=MAX(1,nset/10)
      newsize=newsize*nset   ! prevent INTEGER(4) overflow
      ALLOCATE(icontrib(nset),na(nset),xa(nset),ccbin(nbin,newsize),siga(nset),ifound(nset,nset))
      icontrib=0
      na=0
      xa=0
      siga=0
      ifound=0
      ipointer=0
! read reflections
      DO
        READ(unit,1000) line
        IF (line(1:12) == '!END_OF_DATA') last=.TRUE. 
        IF (.NOT.last) THEN
          READ(line,*)ih,ik,il,fsq,sigma,dum,iset
          CALL getres(ih,ik,il,s2,resol)
          IF (resol>reslim1.OR.resol<reslim2) CYCLE
          IF (sigma<0.) THEN
            misfits = misfits + 1
            CYCLE
          END IF 
          IF (.NOT.weighted) sigma=1  ! switch off weighting
          CALL asuput(ih,ik,il,jh,jk,jl,ihkl,iop,op,SIZE(op,3)/ncenop)
          nobs=nobs+1
          IF (nobs==1) THEN   ! initialize
            oldhkl=ihkl
            s2save=s2
          END IF
        ELSE IF (nobs==0) THEN  ! does this ever happen?
          EXIT ! the reflection loop
        END IF
        IF (ihkl /= oldhkl.OR.last) THEN  ! process unique reflection
! find resolution bin
	  CALL getbin(s2save,ibin)
          nref=nref+1                  ! update # of unique refs
          IF (MOD(nref,modprint)==0) THEN
            PRINT*,nref,' unique reflections processed'
            modprint=modprint*2
          END IF
! isomorphous signal: store contributions of this unique reflection
          DO i=1,ncontrib-1
            kset=icontrib(i)   ! dataset number
            xx4=xa(kset)/siga(kset)
            DO j=i+1,ncontrib
              jset=icontrib(j)
              IF (jset<kset) THEN   ! revert order if necessary. KD 29.1.2020
                jset1=kset
                kset1=jset
              ELSE
                jset1=jset
                kset1=kset
              END IF
              yy4=xa(jset)/siga(jset)
              IF (ifound(kset1,jset1)==0) THEN
                ipointer=ipointer+1
                IF (ipointer>SIZE(ccbin,DIM=2,KIND=8)) THEN
                  newsize=nset
                  newsize=newsize*nset  ! this allows for INTEGER(8) values
                  newsize=MIN(newsize,2*SIZE(ccbin,DIM=2,KIND=8))
                  ALLOCATE(temp_ccbin(nbin,newsize))
                  temp_ccbin(:,1:SIZE(ccbin,DIM=2,KIND=8))=ccbin
                  CALL MOVE_ALLOC(FROM=temp_ccbin,TO=ccbin)
                  PRINT*,'allocating',((nbin*newsize)*24)/(1024*1024),' MB'
                END IF
                ifound(kset1,jset1)=ipointer
              END IF
              CALL put_cc(ccbin(ibin,ifound(kset1,jset1)),xx4,yy4)
            END DO
          END DO
! isomorphous signal done
          IF (last) EXIT  ! the big loop reading the XSCALE file    
          DO i=1,ncontrib   ! clear for next unique ref
            kset=icontrib(i)
            na(kset)=0
            xa(kset)=0.
            siga(kset)=0.
          END DO
          ncontrib=0
          oldhkl=ihkl
          s2save=s2
        END IF
! store this observation
        IF (na(iset)==0) THEN  ! new iset
          ncontrib=ncontrib+1
          icontrib(ncontrib)=iset
        END IF
        xa(iset) = xa(iset) + fsq/sigma**2
        na(iset) = na(iset) + 1
        siga(iset) = siga(iset) + 1./sigma**2
      END DO   ! the big loop reading the XSCALE file 
      CLOSE(unit)
      
! output

      WRITE(*,*) '# obs (w/o misfits), unique, misfits =',nobs,nref,misfits
      WRITE(*,*) 'data between',reslim1,' and',reslim2,'A will be used'
      WRITE(*,*) 'resolution shell limits:'
      CALL getlim(vvr)
      WRITE(*,'(40f7.3)') vvr
      WRITE(*,*) 'mean correlation coefficient in resolution shells:'
      vvr=0
      ntot=0
      ccarray=-2.
      DO i=1,nset-1
        DO j=i+1,nset
          IF (ifound(i,j)==0) CYCLE
          yy4=0.
          n=0
          DO k=1,nbin
            IF (ccbin(k,ifound(i,j))%n > 3) THEN
              xx4 = get_cc(ccbin(k,ifound(i,j)))
              yy4 = yy4 + ccbin(k,ifound(i,j))%n * xx4
              n  = n  + ccbin(k,ifound(i,j))%n
              vvr(k) = vvr(k)  + xx4
              ntot(k)= ntot(k) + 1
            END IF
          END DO
          IF (n>0) THEN
            nccarray(i,j)=n
            nccarray(j,i)=n
            ccarray(i,j)=yy4/n
            ccarray(j,i)=ccarray(i,j)
            IF (n<41) narray(n)=narray(n)+1
          END IF
        END DO
      END DO
      WRITE(*,'(40f7.3)') vvr/MAX(1,ntot)   
      WRITE(*,*) 
      WRITE(*,'(a)') ' Frequency of number of reflections n leading to pairwise correlation coefficients:'
!      WRITE(*,'(a,i6,9i7,20i6)') 'n=',(i,i=4,21)
!      WRITE(*,'(  i8,9i7,20i6)',iostat=ier) narray(4:21)
      WRITE(*,'(a ,18i7)'   ) 'n=',(i,i=4,21)
      WRITE(*,'(2x,18es7.0)') REAL(narray(4:21))

END SUBROUTINE nxds_file
END MODULE nxdsitems_mod

SUBROUTINE outcc(nset,jdim,ccarray,nccarray)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: nset,jdim
      REAL, INTENT(INOUT) :: ccarray(nset,nset)
      INTEGER, INTENT(INOUT) :: nccarray(nset,nset)
      INTEGER, ALLOCATABLE :: mapping(:)
      INTEGER :: i,j,iset,unit,n,isetold=0,notminus2(jdim-1)
      INTEGER(8) :: n8=0
      
      ALLOCATE(mapping(nset))
      
! Iteratively remove images that correlate with less than jdim other images
! It may still happen that there are two or more groups of images that
! only correlate within their own group. CC_ANALYSIS checks this. 
      WRITE(*,*) 
      WRITE(*,'(a,i0,a)') 'Iteratively remove images that correlate with less than <dim>=',jdim,' other images -'
      DO
        mapping=0
        iset=0
        DO i=1,nset
          n=0
          DO j=1,nset
            IF (ccarray(j,i)>-2.) THEN
              n=n+1
              IF (n<jdim) notminus2(n)=j
            END IF
          END DO
          IF (n>=jdim) THEN
            iset=iset+1
            mapping(i)=iset
          ELSE
            DO j=1,n
              ccarray(notminus2(j),i)=-2.
              ccarray(i,notminus2(j))=-2.
            END DO
          END IF
        END DO
        IF (isetold==iset) EXIT
        WRITE(*,*)'Images with enough information for CC_ANALYSIS:',iset
        isetold=iset
      END DO
      
      OPEN(NEWUNIT=unit,FILE='nxds_getcc.dat',STATUS='unknown')
      DO i=1,nset-1
        IF (mapping(i)==0) CYCLE
        DO j=i+1,nset
          IF (mapping(j)==0.OR.ccarray(i,j)<-1.5) CYCLE
          IF (ccarray(i,j)>=0.) THEN
            WRITE(unit,'(2(i0,1x),f6.4,1x,i0)')mapping(i),mapping(j),ccarray(i,j),nccarray(i,j)
          ELSE
            WRITE(unit,'(2(i0,1x),f6.3,1x,i0)')mapping(i),mapping(j),ccarray(i,j),nccarray(i,j)
          END IF
          n8=n8+1
        END DO
      END DO
      CLOSE(unit)
      OPEN(NEWUNIT=unit,FILE='nxds_getcc.map',STATUS='unknown')
      DO i=1,nset
        IF (mapping(i)>0) WRITE(unit,'(i0,1x,i0)') i,mapping(i)
      END DO
      
      WRITE(*,*) 
      WRITE(*,*)'Frequency of images that correlate with n others, after iterative removal:'
      mapping=0
      DO i=1,nset
        j=COUNT(ccarray(:,i)>-2.)
        IF (j>0) mapping(j)=mapping(j)+1
      END DO
      WRITE(*,'(a,i2,29i4)')        'n=',(i,i=1,30)
      WRITE(*,'(     30i4)',iostat=j) mapping(1:30)

      WRITE(*,*) 
      WRITE(*,'(a,i0,a)')'nxds_getcc.dat written with ',n8,' lines'
      WRITE(*,'(a,i0,a)')'nxds_getcc.map written with ',iset,' lines'

END SUBROUTINE outcc

PROGRAM nxds_getcc
      USE resolution_mod, ONLY: setbin
      USE nxdsitems_mod, ONLY:get_nset_reso,resmin,resmax,nset,nxds_file
      IMPLICIT NONE
      REAL :: reslim1=0,reslim2=0
      CHARACTER (LEN=132) :: arg='',name=''
      INTEGER :: nbin=1,mode=2,i,narg,jdim=3
      LOGICAL :: weighted=.TRUE.,skiparg=.FALSE.
      INTEGER,DIMENSION(8) :: itime1
      REAL, ALLOCATABLE :: ccarray(:,:)
      INTEGER, ALLOCATABLE :: nccarray(:,:)

      narg=COMMAND_ARGUMENT_COUNT()
      WRITE(*,'(a)') 'nxds_getcc KD 2020-01-29. -h option shows options'
      WRITE(*,'(a)') 'Academic use only; no redistribution. Expires 2021-12-31'
      CALL DATE_AND_TIME(VALUES=itime1)
      IF (itime1(1)*10000 + itime1(2)*100 + itime1(3)>20211231) STOP 'expired'
      IF (narg==1) THEN
        CALL GET_COMMAND_ARGUMENT(1, arg)
        IF (arg(1:2)=='-h') skiparg=.TRUE.
      END IF
      IF (narg==0.OR.skiparg) THEN
        WRITE(*,'(a,a)') 'usage: nxds_getcc -dmin <lowres> -dmax <highres> -nbin <nbin> -mode <1 or 2> -dim <dim>',&
                ' [-w] nXDS_ASCII_name'
        WRITE(*,'(a)') 'dmax, dmin and nbin (default 1) have the usual meanings. The default'
        WRITE(*,'(a)') '  of dmax, dmin is the common resolution range of all datasets in nXDS_ASCII_name.'  
        WRITE(*,'(a)') 'mode can be 1 (equal volumes of resolution shells) or 2 (increasing volumes; default).'
        WRITE(*,'(a)') 'dim is number of dimensions (default 3).'
        WRITE(*,'(a)') ' -w: no weighting of intensities with their sigmas'
        STOP
      END IF

      DO i=1,narg
        IF (skiparg) THEN
          skiparg=.FALSE.
          CYCLE
        END IF
        CALL GET_COMMAND_ARGUMENT(i, arg)
        IF      (arg(1:5)=='-DMIN'.OR.arg(1:5)=='-dmin') THEN
          IF (arg(6:6)=='=') THEN
            READ(arg(7:),*) reslim2
          ELSE IF (i<narg) THEN
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) reslim2
            skiparg=.TRUE.
          ENDIF
        ELSE IF (arg(1:5)=='-DMAX'.OR.arg(1:5)=='-dmax') THEN
          IF (arg(6:6)=='=') THEN
            READ(arg(7:),*) reslim1
          ELSE IF (i<narg) THEN
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) reslim1 
            skiparg=.TRUE.
          END IF
        ELSE IF (arg(1:5)=='-NBIN'.OR.arg(1:5)=='-nbin') THEN
          IF (arg(6:6)=='=') THEN
            READ(arg(7:),*) nbin
          ELSE IF (i<narg) THEN 
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) nbin 
            skiparg=.TRUE.
          END IF
        ELSE IF (arg(1:5)=='-MODE'.OR.arg(1:5)=='-mode') THEN
          IF (arg(6:6)=='=') THEN
            READ(arg(7:),*) mode
          ELSE IF (i<narg) THEN 
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) mode 
            skiparg=.TRUE.
          END IF
        ELSE IF (arg(1:4)=='-DIM'.OR.arg(1:4)=='-dim') THEN
          IF (arg(5:5)=='=') THEN
            READ(arg(6:),*) jdim
          ELSE IF (i<narg) THEN 
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) jdim 
            WRITE(*,'(a,i2)')' dimension is',jdim
            skiparg=.TRUE.
          END IF
        ELSE IF (arg(1:1)=='-') THEN
          IF (INDEX(arg,'w')>0) THEN
            weighted=.FALSE.
            WRITE(*,*)'option -w found: sigma weighting switched off'
          END IF
        ELSE  ! file name
          name=arg
        END IF
      END DO 
      IF (weighted) WRITE(*,*) 'statistics weighted by sigma values'
      WRITE(*,*)
      WRITE(*,'(a,a)') 'nXDS_ASCII file is ',TRIM(name)
      WRITE(*,'(a,i0)')'number of resolution shells: nbin= ',nbin
      WRITE(*,'(a,i0)')'number of dimensions: dim= ',jdim
      WRITE(*,'(a,i0)')'type of resolution shells: mode= ',mode
      CALL get_nset_reso(name)
      IF (reslim1==0.) reslim1=resmin
      IF (reslim2==0.) reslim2=resmax
      ALLOCATE(ccarray(nset,nset),nccarray(nset,nset))
      CALL setbin(nbin,max(reslim1,reslim2),min(reslim1,reslim2),mode)
      CALL nxds_file(name,max(reslim1,reslim2),min(reslim1,reslim2),nbin,&
           weighted,nset,ccarray,nccarray)
      CALL outcc(nset,jdim,ccarray,nccarray)
      WRITE(*,'(a,i0,a)') 'ready for e.g. cc_analysis -dim ',jdim,' nxds_getcc.dat out.dat'

END PROGRAM nxds_getcc
