! Calculate, for each dataset in XSCALE.HKL file, its CC against all others.
! Based on that information, analyse differences in multi-dimensions.
! Analyze the result and cluster, for further XSCALE processing.
!
! version 0.1 Kay Diederichs 25.9.2015
! version 0.2 Kay Diederichs 10.7.2017
! 13.4.2022 avoid floating overflow in sigma**2 calculation (has actually been seen) - untested code!!!!!!!!!
!
! ifort -ipo -qopenmp-link=static -traceback -warn -init=snan,arrays -fpe0 -static-intel -sox xscale_isocluster.f90 multidimensional-scaling_gelsy_subs.f90 resolution.f90 symops.f90 cc_mod.f90 myqsort.f90 -lmkl_lapack95_lp64 -mkl -o xscale_isocluster.fast
! ifort -qopenmp-link=static -check all,noarg_temp_created -traceback -warn all,noexternal -init=snan,arrays -fpe0 -static-intel -sox xscale_isocluster.f90 multidimensional-scaling_gelsy_subs.f90 resolution.o symops.o cc_mod.o myqsort.o -lmkl_lapack95_lp64 -mkl -o xscale_isocluster
!
! with previous compilation like: 
! ifort -c symops.f90 ../getrc.f90 ../resolution.f90 ../hash.f90 ../cc_mod.f90  myqsort.f90
! 
!
MODULE cross_mod
      REAL, ALLOCATABLE :: cc(:,:),ccano(:,:)
END MODULE cross_mod

MODULE xscaleitems_mod
      INTEGER nset,miniset
      REAL resmin,resmax
      CHARACTER*132, ALLOCATABLE :: inputfile(:)
      LOGICAL :: friedlaw=.TRUE.
      INTEGER :: spgnum=0
contains
      SUBROUTINE get_nset_reso(name,nbin)
      USE symops_mod, ONLY: getops,asuput, op,ncenop
      USE resolution_mod, ONLY: getres,a,b,c,alf,bet,gam
      IMPLICIT NONE
      CHARACTER*(*), INTENT(IN) :: name
      INTEGER, INTENT(INOUT) :: nbin
      CHARACTER line*80
      INTEGER :: ier,unit,iset,ih,ik,il,jh,jk,jl,iop,ihkl,oldhkl,ntot,nuniq
      REAL :: s2,resol,dum(6)
      REAL, ALLOCATABLE :: resminset(:),resmaxset(:)
 1000 FORMAT(a)
      
      miniset=HUGE(miniset)
      nset=0
      ntot=0
      nuniq=0
      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=',line(:LEN_TRIM(line))
          STOP 'read error in xscale_file'
        END IF
        IF (line(1:8)=='!FORMAT=') then
          IF (INDEX(line,'LAW=FALSE')>0) friedlaw=.FALSE.
        ELSE IF (line(1:18)=='!NUMBER_OF_IMAGES=') THEN
          READ(line(19:),*) nset 
          miniset=1
          WRITE(*,'(a)') TRIM(line)
        ELSE IF (line(1:7)=='! ISET=') then
          READ(line(8:),*) iset
          nset=MAX(nset,iset)
! to cope with XSCALE runs which result in output files with lowest ISET > 1 :         
          miniset=MIN(miniset,iset)    
        ELSE IF (line(1:20).eq.'!SPACE_GROUP_NUMBER=') then 
          READ(line(22:),*) spgnum
          WRITE(*,'(a)') TRIM(line)
          CALL getops(spgnum)
        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:14) == '!END_OF_HEADER') THEN
          EXIT
        END IF
      END DO
      IF (spgnum==0) STOP 'error: did not find !SPACE_GROUP_NUMBER='
      IF (a==0.) STOP 'error: did not find !UNIT_CELL_CONSTANTS='
      nset=nset-miniset+1
      ALLOCATE(inputfile(nset),resminset(nset),resmaxset(nset))
      resmin=HUGE(resmin)
      resmax=-1
      resminset=HUGE(resminset)
      resmaxset=-resminset
      DO
        READ(unit,1000,iostat=ier) line
        IF (ier/=0) THEN
          PRINT*,'ier=',ier
          PRINT*,'line=',TRIM(line)
          STOP 'read error in xds_file'
        END IF
        IF (line(1:12) == '!END_OF_DATA') EXIT
        READ(line,*,iostat=ier)ih,ik,il,dum,iset
        IF (ier/=0) THEN
          PRINT*,'ier=',ier
          PRINT*,'line=',TRIM(line)
          STOP 'read error 2 in xds_file'
        ELSE IF (dum(2)==0.) THEN
          PRINT*,'sigma=0'
          PRINT*,'line=',TRIM(line)
          STOP 'data error in xds_file'
        END IF
!          IF (ncheck==5) THEN
!            iset=1     ! reading XDS_ASCII.HKL, not XSCALE.HKL
!          ELSE
            iset=iset-miniset+1
!          END IF
! find min and max resolution
        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
! count total and unique reflections:
        CALL asuput(ih,ik,il,jh,jk,jl,ihkl,iop,op,size(op,3)/ncenop)
        ntot=ntot+1
        IF (ntot==1) oldhkl=ihkl
        IF (oldhkl/=ihkl) THEN
          nuniq=nuniq+1
          oldhkl=ihkl
        END IF
      END DO   
      CLOSE(unit)
      nuniq=nuniq+1  ! don't forget the last unique reflection!
      WRITE(*,'(a,2f10.3)')'highest and lowest resol of any data set        :',resmin,resmax
      resmin=MAXVAL(resminset,MASK=resminset<HUGE(resminset))
      resmax=MINVAL(resmaxset,MASK=resmaxset>0.) 
      WRITE(*,'(a,2f10.3)')'highest and lowest COMMON resol of all data sets:',resmin,resmax
      WRITE(*,'(a,2i6)')  'data set numbers with highest low- and lowest high-resol:',&
                            MAXLOC(resminset)+miniset-1,MINLOC(resmaxset)+miniset-1
      WRITE(*,'(a,2i8)') 'total number of reflections, unique =',ntot,nuniq
      IF (nbin==0) THEN
        dum(1)=REAL(ntot)/(nuniq*nset)
        IF (dum(1)<1.) THEN
          nbin=1
        ELSE IF (dum(1)>2.) THEN
          nbin=10
        ELSE
          nbin=1 + NINT(9*dum(1)-9)
        END IF
        WRITE(*,'(a,f10.4,a,i0)') 'calculating nbin from x=(average multiplicity per data set): (1 if x<1, 10 if x>2) x=',&
                                  dum(1),' nbin=',nbin
      END IF
      END SUBROUTINE get_nset_reso
END MODULE xscaleitems_mod

SUBROUTINE xscale_file(unit,reslim1,reslim2,nbin,abcdef,nset)
      USE cross_mod, ONLY: cc,ccano
      USE symops_mod, ONLY: getops,asuput,is_centric, op,ncenop
      USE resolution_mod, ONLY: getres,getbin,getlim, a,b,c,alf,bet,gam
      USE cc_mod, ONLY: cc_type, get_cc,put_cc
      USE xscaleitems_mod, ONLY: inputfile,miniset
      IMPLICIT NONE
      
      INTEGER, INTENT(IN) :: unit,nbin
      REAL, INTENT(IN) :: reslim1,reslim2
      LOGICAL, INTENT(IN) :: abcdef(10)
      INTEGER, INTENT(IN) :: 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,ncommon(nbin)
      REAL :: vvr(nbin),s2,resol,fsq,sigma,dum(4),s2save,xx4,yy4
      CHARACTER :: line*200
      LOGICAL :: last=.FALSE.,centric
      INTEGER, ALLOCATABLE :: nap(:),nam(:),icontrib(:)            
      REAL, ALLOCATABLE :: xap(:),xam(:),sigap(:),sigam(:)
      TYPE(cc_type), ALLOCATABLE :: ccbin(:,:,:),ccanobin(:,:,:)

 1000 FORMAT(A)
      ntot=0
      ncheck=0
! Open input file; get spacegroup and cell constants
      DO
        READ(unit,1000,iostat=ier) line
        IF (ier/=0) THEN
          PRINT*,'ier=',ier
          PRINT*,'line=',line(:LEN_TRIM(line))
          STOP 'read error in xscale_file'
        END IF
        IF (line(1:20)=='!SPACE_GROUP_NUMBER=') then
          READ(line(21:),*) ispgr
          CALL getops(ispgr)
        ELSE IF (line(1:21).eq.'!UNIT_CELL_CONSTANTS=') then 
          READ(line(22:),*) a,b,c,alf,bet,gam
        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:13) == '!ITEM_ISET=10') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:12) == '!ITEM_IMG=10') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:7) == '! ISET=') THEN
          READ(line(8:),*) iset
          iset=iset-miniset+1
          i=INDEX(line,'INPUT_FILE=')
          IF (i>0) inputfile(iset)=line(i+11:)
        ELSE IF (line(1:14) == '!END_OF_HEADER') THEN
          EXIT
        END IF
      END DO
      IF (ncheck <6) STOP 'header items not found'  
      print*,'# of datasets=',nset
      
      ALLOCATE(icontrib(nset),nap(nset),nam(nset),xap(nset),xam(nset),cc(nset,nset),&
              sigap(nset),sigam(nset),ccbin(nbin,nset,nset),ccanobin(nbin,nset,nset),&
              ccano(nset,nset))
      icontrib=0
      nap=0
      nam=0
      xap=0
      xam=0
      sigap=0
      sigam=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
          IF (ncheck==5) THEN
            iset=1     ! reading XDS_ASCII.HKL, not XSCALE.HKL
          ELSE
            iset=iset-miniset+1
          END IF
          CALL getres(ih,ik,il,s2,resol)
!          IF (resol<resmin) resmin=resol
!          IF (resol>resmax) resmax=resol
          IF (resol>reslim1.OR.resol<reslim2) CYCLE
          IF (sigma<0.) THEN
            misfits = misfits + 1
            CYCLE
          END IF 
          IF (.NOT.abcdef(7)) THEN  ! switch off weighting
            sigma=1
          END IF
          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
          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
          IF (abcdef(1).AND.abcdef(10)) THEN   ! no -a or -A option
! isomorphous signal: store contributions of this unique reflection; I+ and I- averaged
            DO i=1,ncontrib-1
              kset=icontrib(i)   ! dataset number
              xx4=(xap(kset)+xam(kset))/(sigap(kset)+sigam(kset))
              DO j=i+1,ncontrib
                jset=icontrib(j)
                IF (jset<kset) STOP 'error jset<kset'
                yy4=(xap(jset)+xam(jset))/(sigap(jset)+sigam(jset))
                CALL put_cc(ccbin(ibin,kset,jset),xx4,yy4)
              END DO
            END DO
! isomorphous signal done
          ELSE IF (.NOT.abcdef(10)) THEN  ! -A option 
! treat as isomorphous signal but keep I+ and I- separate
            DO i=1,ncontrib-1
              kset=icontrib(i)   ! dataset number
              IF (sigap(kset)==0.) CYCLE
              xx4=xap(kset)/sigap(kset)
              DO j=i+1,ncontrib
                jset=icontrib(j)
                IF (jset<kset) STOP 'error jset<kset'
                IF (sigap(jset)==0.) CYCLE
                yy4=xap(jset)/sigap(jset)
                CALL put_cc(ccbin(ibin,kset,jset),xx4,yy4)
              END DO
            END DO
            DO i=1,ncontrib-1
              kset=icontrib(i)   ! dataset number
              IF (sigam(kset)==0.) CYCLE
              xx4=xam(kset)/sigam(kset)
              DO j=i+1,ncontrib
                jset=icontrib(j)
                IF (jset<kset) STOP 'error jset<kset'
                IF (sigam(jset)==0.) CYCLE
                yy4=xam(jset)/sigam(jset)
                CALL put_cc(ccbin(ibin,kset,jset),xx4,yy4)
              END DO
            END DO
! "isomorphous signal" with separate I+ and I- done
          ELSE IF (.NOT.abcdef(1)) THEN  ! -a option 
! anomalous signal: store contributions to this unique reflection
            DO i=1,ncontrib-1
              kset=icontrib(i)   ! dataset number
              IF (nap(kset)==0.OR.nam(kset)==0) CYCLE
              xx4=xap(kset)/sigap(kset) - xam(kset)/sigam(kset)
              DO j=i+1,ncontrib
                jset=icontrib(j)
                IF (nap(jset)==0.OR.nam(jset)==0) CYCLE
                yy4=xap(jset)/sigap(jset) - xam(jset)/sigam(jset)
                CALL put_cc(ccanobin(ibin,kset,jset),xx4,yy4)
              END DO
            END DO
! anomalous signal done
          END IF
          
          IF (last) EXIT  ! the big loop reading the XSCALE file    
          DO i=1,ncontrib   ! clear for next unique ref
            kset=icontrib(i)
            nap(kset)=0
            xap(kset)=0.
            nam(kset)=0
            xam(kset)=0.
            sigap(kset)=0.
            sigam(kset)=0.
          END DO
          ncontrib=0
          oldhkl=ihkl
          s2save=s2
        END IF
! store this observation
        IF (ncontrib==0) THEN   ! happens once for each unique ref
          centric=is_centric(jh,jk,jl,op,size(op,3)/ncenop)
        END IF
        IF (nap(iset)+nam(iset)==0) THEN  ! new iset
          ncontrib=ncontrib+1
          icontrib(ncontrib)=iset
        END IF
! in serial crystallography, sigma can become large in XSCALE
! avoid floating overflow by squaring the inverse of sigma 
        sigma=1./sigma
        sigma=sigma**2
        IF (iop>0.OR.centric) THEN  ! this is taken as I+
          xap(iset) = xap(iset) + fsq*sigma ! fsq/sigma**2
          nap(iset) = nap(iset) + 1
          sigap(iset) = sigap(iset) + sigma ! 1/sigma**2
        ELSE                        !  I-
          xam(iset) = xam(iset) + fsq*sigma ! fsq/sigma**2
          nam(iset) = nam(iset) + 1
          sigam(iset) = sigam(iset) + sigma ! 1/sigma**2
        END IF
      END DO   ! the big loop reading the XSCALE file 

! output of statistics

      WRITE(*,*) '# obs (w/o misfits), unique, misfits =',nobs,nref,misfits
!      WRITE(*,*) 'max and min resolution of data in file =',resmax,resmin
      WRITE(*,*) 'data between',reslim1,' and',reslim2,'A will be used'
      WRITE(*,*) 'resolution shell limits:'
      CALL getlim(vvr)
      ncommon=0
      PRINT '(40f7.3)',vvr
      IF (abcdef(1)) THEN
        vvr=0
        ntot=0
        cc=-2
        DO i=1,nset-1
          DO j=i+1,nset
            cc(j,i)=0
            n=0
            DO k=1,nbin
              IF (ccbin(k,i,j)%n > 5) THEN
                xx4 = get_cc(ccbin(k,i,j))
                cc(j,i) = cc(j,i) + ccbin(k,i,j)%n * xx4
                n       = n       + ccbin(k,i,j)%n
                ncommon(k) = ncommon(k) + ccbin(k,i,j)%n
                vvr(k) = vvr(k)  + xx4
                ntot(k)= ntot(k) + 1
              END IF
            END DO
            IF (n>0) THEN
              cc(j,i) = cc(j,i)/n
              cc(i,j) = cc(j,i)
            ELSE
              cc(j,i) = -2
            END IF
          END DO
        END DO
        PRINT*,'mean isomorphous correlation in resolution shells:'
        PRINT '(40f7.3)',vvr/MAX(1,ntot)
      ELSE
        PRINT*,'mean anomalous correlation in resolution shells:'
        vvr=0
        ntot=0
        ccano=-2
        DO i=1,nset-1
          DO j=i+1,nset
            ccano(j,i)=0
            n=0
            DO k=1,nbin
              IF (ccanobin(k,i,j)%n > 5) THEN
                xx4 = get_cc(ccanobin(k,i,j))
                ccano(j,i) = ccano(j,i) + ccanobin(k,i,j)%n * xx4
                n          = n          + ccanobin(k,i,j)%n
                ncommon(k) = ncommon(k) + ccbin(k,i,j)%n
                vvr(k) = vvr(k)  + xx4
                ntot(k)= ntot(k) + 1
              END IF
            END DO
            IF (n>0) THEN
              ccano(j,i) = ccano(j,i)/n
              ccano(i,j) = ccano(j,i)
            ELSE
              ccano(j,i) = -2
            END IF
          END DO
        END DO
        PRINT '(40f7.3)',vvr/MAX(1,ntot)
        cc=ccano
      END IF
      PRINT*,'average number of common reflections (used for CC) in resolution ranges:'
      PRINT '(40i6)',ncommon/MAX(1,ntot)

END SUBROUTINE xscale_file

SUBROUTINE call_cluster(result,nset,jdim,cl,numclus)
      USE qsort_m, ONLY: myqsort
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: nset,jdim,numclus ! numclus=desired # of clusters
      REAL, INTENT(IN) :: result(jdim,nset)
      INTEGER, INTENT(OUT) :: cl(nset)  ! return cluster number; negative means "halo"
! local variables
      REAL, ALLOCATABLE :: dist(:,:),rho(:),delta(:)
      INTEGER, ALLOCATABLE :: icl(:),isort(:),nsort(:)
      INTEGER i,j,k,nclust,nc,nh
      REAL :: xx,yy !zz
      
      ALLOCATE(dist(nset,nset),rho(nset),delta(nset),icl(nset))	
      DO k=1,nset
	dist(k,k)=0.
        xx=SQRT(DOT_PRODUCT(result(:,k),result(:,k)))  ! length of k
        IF (xx<EPSILON(xx)) THEN  ! result(:,i) set to zero in mds_scalarproduct_iter
          PRINT*,'sorry, must stop because 1 or more datasets (see numbers given above)'
          PRINT*,'does not have enough reflections in common with other datasets'
          STOP ' remove these from XSCALE.INP and re-run XSCALE!'
        END IF
	DO i=1,k-1
          yy=SQRT(DOT_PRODUCT(result(:,i),result(:,i)))  ! length of i
          yy=DOT_PRODUCT(result(:,i),result(:,k))/(xx*yy) ! cos(angle(k,i))
          IF (ABS(yy)>1.) yy=yy/ABS(yy)  ! avoid acos problem
          dist(i,k)=acos(yy)   ! angle(k,i)
!          yy=ATAN2(result(2,i),result(1,i))
!          zz=ATAN2(result(2,k),result(1,k))
!write(101,'(2i4,4f12.4)') i,k,dist(i,k),ABS(yy-zz),yy,zz
          dist(k,i)=dist(i,k)
	END DO
      END DO
IF (numclus==1) THEN
  cl=1
  RETURN
END IF

      CALL findcluster(nset,dist,rho,delta,cl,nclust,icl,"gaussian",numclus) !"gaussian" o

      ALLOCATE(isort(nclust),nsort(nclust))  ! for sorting by size
      DO i=1,nclust
        nc=0
        nh=0
        DO j=1,nset
          IF (ABS(cl(j))==i) THEN
            nc=nc+1
            IF (cl(j)>0) nh=nh+1
          END IF
        END DO
        isort(i)=i
        nsort(i)=-nc
        WRITE(*,'(''cluster:'',i3,'' center:'',i6,'' elements:'',i6,'' core:'',i6,'' halo:'',i6)')&
                  i,icl(i),nc,nh,nc-nh
      END DO
! sort clusters by size, first=biggest
      CALL myqsort(nclust,nsort,isort)
! reverse lookup table
      DO i=1,nclust
        nsort(isort(i))=i
      END DO
! now use isort as copy of icl
      isort=icl(:nclust)
      DO i=1,nclust
        icl(nsort(i))=isort(i)  ! cluster center
      END DO
     
! change cl array accordingly
      DO i=1,nset
        j=SIGN(1,cl(i))         ! j is -1 or 1
        cl(i)=j*nsort(ABS(cl(i)))
      END DO
! repeat output but now sorted
      PRINT*,' '
      PRINT*,'sorted by cluster size:'
      PRINT*,' '
      DO i=1,nclust
        nc=0
        nh=0
        DO j=1,nset
          IF (ABS(cl(j))==i) THEN
            nc=nc+1
            IF (cl(j)>0) nh=nh+1
          END IF
        END DO
        WRITE(*,'(''cluster:'',i3,'' center:'',i6,'' elements:'',i6,'' core:'',i6,'' halo:'',i6)')&
                  i,icl(i),nc,nh,nc-nh
      END DO
      
      
! overwrite the cl array by assigning each point to the closest cluster center
!      DO i=1,nset
!        xx=HUGE(xx)
!        DO j=1,nclust
!          IF (dist(i,icl(j)) < xx) THEN
!            xx=dist(i,icl(j))
!            cl(i)=j
!          END IF
!        END DO
!      END DO
      
END SUBROUTINE call_cluster

! KD 9.9.14 Matlab-clustering algorithm from Science 344 (2014), 1492-1496 
! additionally, find best number of clusters based on the following cutoff decision:
! if a cluster has more halo points than core points, this is no longer a cluster
! compilation: ifort -C -traceback -fpe0 -w findcluster.f90 myqsort.o
! also see http://cran.r-project.org/web/packages/densityClust/
! and http://eric-yuan.me/clustering-fast-search-find-density-peaks/
      SUBROUTINE findcluster(nd,dist,rho,delta,cl,nclust,icl,kernel,numclus)
      USE qsort_m    ! for myqsort and mymedian
      IMPLICIT none
      INTEGER, INTENT(IN) :: nd,numclus
      REAL, INTENT(IN) :: dist(nd,nd)
      REAL, INTENT(OUT) :: rho(nd),delta(nd)
      INTEGER, INTENT(OUT) :: nclust,cl(nd),icl(nd)
      CHARACTER(*), INTENT(IN) :: kernel
      
      REAL, ALLOCATABLE :: xx(:),bord_rho(:),gamma(:)
      INTEGER, ALLOCATABLE :: ordrho(:),nneigh(:),ind(:)
      INTEGER :: i,j,n,ii,jj,nc,nh,itry,nhalo,olddf
      REAL :: dc
      REAL :: maxd,rhomin,deltamin,rho_aver,psi,oldpsi
      LOGICAL :: finish
      
! find typical distance, and max distance
      ALLOCATE(xx((nd*(nd-1))/2))	       
      maxd=0.
      n=0
      oldpsi=0.
      DO i=1,nd-1
        DO j=i+1,nd
	  n=n+1
	  IF (dist(i,j)/=dist(j,i)) STOP 'asymmetric distance matrix?!'
	  xx(n)=dist(i,j)
	  maxd=max(dist(i,j),maxd)
	  oldpsi=oldpsi+dist(i,j)**2
	END DO
      END DO
      oldpsi=oldpsi/((nd-1)/2.)
      WRITE(*,*)'average percentage of neighbours (hard coded): 2'
      i=max(1,nint(n*0.02))
      CALL mymedian(n,i,xx)  
      dc=xx(i)
IF (dc==0.) dc=MINVAL(xx,MASK=xx>0.)  ! otherwise floating invalid below. KD 4.1.18
      DEALLOCATE(xx)
      ALLOCATE(ordrho(nd),nneigh(nd),bord_rho(nd),gamma(nd),ind(nd))
      
      rho=0.
      IF (kernel=='gaussian') THEN  ! change to switch from gaussian to cut-off 
        WRITE(*,*)'Computing Rho with gaussian kernel of radius:',dc  
        DO i=1,nd-1
          DO j=i+1,nd
            rho(i)=rho(i)+EXP(-(dist(i,j)/dc)**2)
            rho(j)=rho(j)+EXP(-(dist(i,j)/dc)**2)
          END DO
        END DO      
      ELSE
        WRITE(*,*)'Computing Rho with cut-off kernel of radius:',dc  
        DO i=1,nd-1
          DO j=i+1,nd
            IF (dist(i,j)<dc) THEN
              rho(i)=rho(i)+1.
              rho(j)=rho(j)+1.
            END IF
          END DO
        END DO
      END IF

      delta=-rho   ! use delta as temporary array, to sort rho descendingly
      ordrho=(/(i,i=1,nd)/)
      CALL myqsort(nd,delta,ordrho) ! only ordrho is used afterwards
      
      delta(ordrho(1))=-1.
      nneigh(ordrho(1))=0
      DO ii=2,nd
        delta(ordrho(ii))=maxd
	nneigh(ordrho(ii))=-1 ! since not all nneigh are usually spec'd, set -1
        DO jj=1,ii-1
          IF (dist(ordrho(ii),ordrho(jj))<delta(ordrho(ii))) THEN
            delta(ordrho(ii))=dist(ordrho(ii),ordrho(jj))
            nneigh(ordrho(ii))=ordrho(jj)   ! not all nneigh are specified here
          END IF
        END DO
      END DO
      delta(ordrho(1))=maxval(delta)
      DO i=1,nd
        ind(i)=i
        gamma(i)=rho(i)*delta(i)
      END DO
      CALL myqsort(nd,gamma,ind)
      WRITE(*,*) 'highest 10 values of rho*delta, and their rho, delta:'
      DO i=nd,MAX(1,nd-9),-1
        WRITE(*,*) nd-i+1,ind(i),gamma(i),rho(ind(i)),delta(ind(i))
      END DO

      IF (numclus>0) THEN
        itry=numclus
        CALL findcluster1
        WRITE(*,*)'psi=',psi,' nhalo=',nhalo
        RETURN
      END IF
! automatic determination of # of clusters:
      finish=.FALSE.
      olddf=0
      DO itry=1,100   ! itry goes through the number of clusters
        CALL findcluster1  
	IF (nclust==olddf) CYCLE
        WRITE(*,*)'psi=',psi,' nhalo=',nhalo,' F=',(oldpsi-psi)/(nclust-olddf)/(psi/(nd-nclust))
! F calculated from http://en.wikipedia.org/wiki/F-test#Regression_problems 16.9.14
! the one-sided critical value of F for 
! p(1,10)=0.05 is 241.9       p(1,10)=0.10 is ~60
! p(1,infinity)=0.05 is 254.3  (the second parameter in the bracket is # points)
! so, (roughly) if F > 100 , the improvement is statistically significant
        oldpsi=psi
	olddf=nclust
	IF (finish) EXIT  ! with this itry, we found #(halo) > #(core)
      END DO
      itry=itry-1
      WRITE(*,*) ''
      WRITE(*,*) 'best number of clusters:',itry
      CALL findcluster1       
      
CONTAINS     
      
      SUBROUTINE findcluster1

      psi=0.
      nhalo=0
      deltamin=huge(deltamin)
      rhomin=huge(rhomin)
      DO i=nd+1-itry,nd
        deltamin=min(deltamin,delta(ind(i)))
        rhomin=min(rhomin,rho(ind(i)))
      END DO
      WRITE(*,*)''
      WRITE(*,*)'deltamin, rhomin=',deltamin,rhomin      
      nclust=0
      DO i=1,nd
        cl(i)=-1
      END DO
      DO i=1,nd
        IF ( (rho(i)>=rhomin) .AND. (delta(i)>=deltamin)) THEN
           nclust=nclust+1
           cl(i)=nclust
           icl(nclust)=i
        END IF
      END DO
      IF (nclust==olddf) RETURN  ! don't repeat old calculation
      WRITE(*,*)'NUMBER OF CLUSTERS:', nclust,' :'
      
! assignment
      DO i=1,nd
        IF (cl(ordrho(i))==-1) THEN
	  IF (nneigh(ordrho(i))==-1) STOP 'assert error nneigh -1'
	  IF (nneigh(ordrho(i))== 0) STOP 'assert error nneigh  0'
          cl(ordrho(i))=cl(nneigh(ordrho(i)))
        END IF
      END DO
     
! find halo, and set the cluster number to negative
      IF (nclust>1) THEN
        DO i=1,nclust
          bord_rho(i)=0.
        END DO
        DO i=1,nd-1
          DO j=i+1,nd
            IF ((cl(i)/=cl(j)) .AND. (dist(i,j)<=dc)) THEN
     	      rho_aver=(rho(i)+rho(j))/2.
     	      IF (rho_aver>bord_rho(cl(i))) THEN
     	        bord_rho(cl(i))=rho_aver
     	      END IF
     	      IF (rho_aver>bord_rho(cl(j))) THEN
     	        bord_rho(cl(j))=rho_aver
     	      END IF
            END IF
          END DO
        END DO
        DO i=1,nd
          IF (rho(i)<bord_rho(cl(i))) THEN
            cl(i)=-cl(i)
          END IF
        END DO
      END IF
      DO i=1,nclust
        nc=0
        nh=0
        DO j=1,nd
          IF (ABS(cl(j))==i) THEN
            nc=nc+1                ! nc counts all cluster points
	    psi=psi+dist(j,icl(i))**2   ! for the F calculation
	    IF (cl(j)>0) nh=nh+1   ! nh counts core points of cluster
          END IF
        END DO
!        WRITE(*,'(''cluster:'',i3,'' center:'',i6,'' elements:'',i6,'' core:'',i6,'' halo:'',i6)')&
!                  i,icl(i),nc,nh,nc-nh
	nhalo=nhalo+nc-nh
        IF (nh < nc-nh) finish=.TRUE.
      END DO
      END SUBROUTINE findcluster1
      END SUBROUTINE findcluster


PROGRAM xscalecluster
      USE cross_mod, ONLY: cc
      USE resolution_mod, ONLY: setbin, a,b,c,alf,bet,gam
      USE xscaleitems_mod, ONLY:get_nset_reso,nset,resmin,resmax, &
                                inputfile,friedlaw,spgnum
      USE qsort_m
      IMPLICIT NONE
      REAL :: reslim1=0,reslim2=0,phi,temp,wscale=1.
      CHARACTER (LEN=132) :: arg='',name=''
      INTEGER :: unit,ier,nbin=0,mode=2,i,j,k,narg,jdim=3,ncluster,numclus=1,icen(9)=0
      INTEGER, ALLOCATABLE :: purpose(:,:),cluster(:),inx(:),iunit(:)
      REAL, ALLOCATABLE :: result(:,:),dist(:),clustcen(:),tempresult(:,:)
      LOGICAL :: abcdef(10)=.TRUE.,skiparg=.FALSE.
      INTEGER,DIMENSION(8) :: itime1

 
      narg=COMMAND_ARGUMENT_COUNT()
      WRITE(*,'(a)') 'xscale_isocluster KD 2020-07-27. -h option shows options'
      WRITE(*,'(a)') 'Licensed for academic use only; no redistribution.' ! Expires 2020-07-31'
      WRITE(*,'(a)') 'Please cite Assmann,Wang and Diederichs (2020) Acta Cryst D76, 636.'

      CALL DATE_AND_TIME(VALUES=itime1)
!      IF (itime1(1)*10000 + itime1(2)*100 + itime1(3)>20200731) 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: xscale_isocluster -dmin <lowres> -dmax <highres> -nbin <nbin> -mode <1 or 2> -dim <dim>',&
                ' -clu <#> -cen <#,#,#,...> -<aAisw> XSCALE_FILE_NAME'
        WRITE(*,'(a)') 'highres, lowres and nbin (default=1 for low, and 10 for high multiplicity) have the usual meanings.'
        WRITE(*,'(a)') 'The default of highres and lowres is the COMMON resolution range of all datasets in XSCALE_FILE.'  
        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)') 'clu (by default 1) is number of clusters.'
        WRITE(*,'(a)') 'cen (by default automatically) is a set of cluster centers (up to 9) specified by their ISET.'
        WRITE(*,'(a)') 'cen must be specified after clu, and the number of values given must match clu.' 
        WRITE(*,'(a)') ' -a: use anomalous (instead of isomorphous) signal only (unsuitable for partial datasets)'
        WRITE(*,'(a)') ' -A: account for anomalous signal by separating I+ and I- (suitable for partial datasets)'
        WRITE(*,'(a)') ' -i: write individual pseudo-PDB files to visualize clusters'
        WRITE(*,'(a)') ' -s: scale (default 1) for WEIGHT=1/cos(scale*angle) values written to XSCALE.*.INP file(s)'
        WRITE(*,'(a)') ' -w: no weighting of intensities with their sigmas'
        WRITE(*,'(a)') 'The XSCALE.?.INP output files have angles [deg] w.r.t. cluster centers'
        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
          WRITE(*,'(a,f5.2)')' dmin is',reslim2
        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
          WRITE(*,'(a,f5.2)')' dmax is',reslim1
        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
          WRITE(*,'(a,i2)')' nbin is',nbin
        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
          WRITE(*,'(a,i2)')' mode for resolution shells is',mode
        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 
            skiparg=.TRUE.
          END IF
          WRITE(*,'(a,i2)')' dimension is',jdim
        ELSE IF (arg(1:4)=='-CLU'.OR.arg(1:4)=='-clu') THEN
          IF (arg(5:5)=='=') THEN
            READ(arg(6:),*) numclus
          ELSE IF (i<narg) THEN 
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) numclus 
            skiparg=.TRUE.
          END IF
          WRITE(*,'(a,i2)')' number of clusters is',numclus
        ELSE IF (arg(1:4)=='-CEN'.OR.arg(1:4)=='-cen') THEN
          IF (arg(5:5)=='=') THEN
            READ(arg(6:),*) icen(:numclus) 
          ELSE IF (i<narg) THEN 
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) icen(:numclus) 
            skiparg=.TRUE.
          END IF
          WRITE(*,'(a,9i6)')' cluster centers are',icen(:numclus)
        ELSE IF (arg(1:2)=='-s') THEN
          IF (arg(3:3)=='=') THEN
            READ(arg(4:),*) wscale 
          ELSE IF (i<narg) THEN 
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) wscale 
            skiparg=.TRUE.
          END IF
          WRITE(*,'(a,f5.2)')' scale for WEIGHT=1/cos(scale*angle) is',wscale
        ELSE IF (arg(1:1)=='-') THEN
          IF (INDEX(arg,'w')>0) THEN
            abcdef(7)=.FALSE.
            WRITE(*,*)'option -w found: sigma weighting switched off'
          END IF
          IF (INDEX(arg,'a')>0) THEN
            abcdef(1)=.FALSE.
            WRITE(*,*)'option -a found: base calculations on anomalous (instead of isomorphous) signal'
          END IF 
          IF (INDEX(arg,'A')>0) THEN
            abcdef(10)=.FALSE.
            WRITE(*,*)'option -A found: consider I+ and I- as separate reflections'
          END IF 
          IF (INDEX(arg,'i')>0) THEN
            abcdef(2)=.FALSE.
            WRITE(*,*)'option -i found: write out pseudo-PDB files to visualize clusters'
          END IF 
        ELSE  ! file name
          name=arg
        END IF
      END DO 
      IF (abcdef(7)) WRITE(*,*) 'statistics weighted by sigma values'
      IF (.NOT.abcdef(1).AND..NOT.abcdef(10)) STOP '-a and -A cannot be used together'
      WRITE(*,*)
      WRITE(*,'(a,a)') 'XSCALE file is ',TRIM(name)
      CALL get_nset_reso(name,nbin)
      IF (reslim1==0.) reslim1=resmin
      IF (reslim2==0.) reslim2=resmax
      CALL setbin(nbin,max(reslim1,reslim2),min(reslim1,reslim2),mode)
      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
      CALL xscale_file(unit,max(reslim1,reslim2),min(reslim1,reslim2),nbin,&
           abcdef,nset)
      CLOSE(unit)
      ALLOCATE(purpose(nset,nset),result(jdim,nset),cluster(nset), &
               inx(nset),dist(nset),tempresult(jdim,nset))
      tempresult=0
! declare some matrix elements as unknown
      purpose=0      ! at first, all matrix elements declared known and useful
      DO i=1,nset
        purpose(i,i)=-1     ! diagonal is always unknown
        DO j=i+1,nset
          IF (cc(j,i)>=-1) CYCLE
          purpose(j,i)=-1   ! unknown matrix elements
          purpose(i,j)=purpose(j,i)
        END DO
      END DO
      phi=0.
      result=0.
      WRITE(*,'(a,i0,1x,i0,a,i0,a)') 'finding ',nset,jdim,'-vectors from ', &
                                  COUNT(purpose==0)/2,' unique correlation coefficients'
      IF (nset*jdim>COUNT(purpose==0)/2) THEN
        WRITE(*,*) 'under-determined. Decrease dim!'
        STOP 'cannot continue'
      END IF
! check if each set has at least dim ccs
      j=0
      DO i=1,nset
        IF (COUNT(cc(:,i)>=-1)<jdim) THEN
          WRITE(*,*)'data set ',i,' has only ', COUNT(cc(:,i)>=-1),' data sets it overlaps with'
          j=j+1
        END IF
      END DO
      IF (j>0) STOP 'cannot continue'
      CALL solve_mds_neu(nset,jdim,cc,purpose,result)
! align result with x,y,z axes
!      CALL inertia_align(nset,jdim,result)
! make the first result vector have positive components in all directions
      DO i=1,jdim
        IF (result(i,1)<0) result(i,:)=-result(i,:) ! this conserves scalar products
      END DO
      CALL wpdb(100*result,nset,jdim,'iso.pdb','A')
      PRINT*,'in COOT, use Draw/Cell & Symmetry/Show unit cells - the 0 is the origin'

! cluster:    
      CALL call_cluster(result,nset,jdim,cluster,numclus)
      ncluster=maxval(ABS(cluster))
      IF (numclus>0) ncluster=MIN(ncluster,numclus)
      IF (ncluster==1) THEN
        PRINT*,'Looking for a single cluster (-clu 1). No halo (outlier) info available.'
      ELSE IF (ncluster>9) THEN
        PRINT*,'there are',ncluster,' clusters. Only 9 will be processed further.'
      ELSE
        PRINT*,'there are',ncluster,' clusters.'
      END IF
      IF (ncluster>9) ncluster=9
      ALLOCATE(clustcen(jdim),iunit(ncluster))
! prepare XSCALE.INP.# files
      DO i=1,ncluster
! write out iso.#.pdb 
        WRITE(name,'(a,i1,a)')'iso.',i,'.pdb'
        k=0
        DO j=1,nset
          IF (ABS(cluster(j))==i) THEN
            k=k+1
            tempresult(:,k)=result(:,j)
          END IF
        END DO
        IF (.NOT.abcdef(2)) CALL wpdb(100*tempresult,k,jdim,name,CHAR(MIN(90,64+i)))   ! segid=A...Z
        IF (jdim==4.AND.i==1) PRINT*,'the B-factor column has 100+4th_coordinate'
        IF (jdim>=5.AND.i==1) PRINT*,'the B-factor column has 100+4th_coordinate, and occ is 1+5th_coordinate/100' 
! write XSCALE.#.INP
        WRITE(name,'(a,i1,a)')'XSCALE.',i,'.INP'
        OPEN(NEWUNIT=iunit(i),FILE=name)
        WRITE(iunit(i),'(a,6f9.3)') 'UNIT_CELL_CONSTANTS=',a,b,c,alf,bet,gam
        WRITE(iunit(i),'(a,i4)') 'SPACE_GROUP_NUMBER=',spgnum
        WRITE(iunit(i),'(a,i1,a)') 'OUTPUT_FILE=XSCALE.',i,'.HKL'
        IF (.NOT.friedlaw) WRITE(iunit(i),'(a)') 'FRIEDEL''S_LAW=FALSE'
        WRITE(iunit(i),'(a)') 'SAVE_CORRECTION_IMAGES=FALSE'
        WRITE(iunit(i),'(a)') 'PRINT_CORRELATIONS=FALSE'
        WRITE(iunit(i),'(a)') 'WFAC1=1.25 ! XDS/XSCALE defaults are 1.0/1.5'
      END DO
! give the user information about each data set
!      WRITE(*,*)' set  strength isomorphism cluster'
!      WRITE(*,*)'    (1:strong) (0:center) (<0:halo)'
!      DO i=1,nset
!        temp=SQRT(DOT_PRODUCT(result(:,i),result(:,i)))
!        WRITE(*,'(i5,2f10.3,i7)') i,temp,result(2,i), cluster(i)

!!        WRITE(*,'(i5,2f10.3,i7)') i,MIN(1.,temp), ATAN2(result(2,i),result(1,i)), cluster(i)
!!        WRITE(*,'(i5,2f10.3,i7)') i,MIN(1.,temp), result(1,i)/temp, cluster(i)
!      END DO
      DO i=1,ncluster
! find centers of clusters:
        clustcen=0
        k=0    
        DO j=1,nset
          IF (DOT_PRODUCT(result(:,j),result(:,j))<0.0001) CYCLE  ! too short
          IF (ABS(cluster(j))==i) THEN    ! also use halo
            clustcen(:)=clustcen(:) + result(:,j)
            k = k + 1
          END IF
        END DO
        IF (k>0) THEN   
          clustcen=clustcen/k        !  center of cluster
          IF (icen(i)>0) clustcen=result(:,icen(i))  ! user-specified cluster center
          temp=SQRT(DOT_PRODUCT(clustcen,clustcen))
          clustcen=clustcen/temp     !  set length=1
          WRITE(*,'(a,i0,a,10f8.2)') 'direction of cluster ',i,' is ',clustcen 
! sort cluster members by ABS(distance) to center, and write to XSCALE.INP.#
          k=0
          DO j=1,nset
            temp=SQRT(DOT_PRODUCT(result(:,j),result(:,j)))
            IF (temp<0.1) cluster(j)=-ABS(cluster(j))   ! mark as halo
            IF (ABS(cluster(j))==i) THEN 
              temp=DOT_PRODUCT(result(:,j),clustcen) / temp ! cos(angle)
              IF (ABS(temp)>1.) temp=temp/ABS(temp)  ! avoid acos problems
              k=k+1
              dist(k)=acos(temp)   ! lowest angle should be sorted to pos 1
              inx(k)=j
            END IF
          END DO
          CALL myqsort(k,dist,inx)
          DO j=1,k
            WRITE(iunit(i),'(a,a)') 'INPUT_FILE=',TRIM(inputfile(inx(j)))
            temp=SQRT(DOT_PRODUCT(result(:,inx(j)),result(:,inx(j))))
            WRITE(iunit(i),'(a,2i7,a,f10.3,f8.1,i7)') &
                 '!new, old ISET=',j,inx(j),' length=CC*,angle,cluster=',temp,dist(j)*180/3.14159,cluster(inx(j))
!            temp=sumres(ibin)*SQRT(-3./(2.*3.14159**2)*log(MAX(0.01,dist(j))))
!            WRITE(iunit(i),'(a,f10.3,a)') 'RMSD=',temp,'!INCLUDE_RESOLUTION_RANGE=00 00'
            WRITE(iunit(i),'(a)') '!INCLUDE_RESOLUTION_RANGE=00 00'
            temp=1./COS(MIN(1.57,wscale*dist(j)))  ! limit WEIGHT if wscale*dist(j) > acos(0) 
            WRITE(iunit(i),'(a,F10.3)') 'WEIGHT=',temp 
          END DO        
        END IF
      END DO
      DO i=1,ncluster
        WRITE(name,'(a,i1,a)')'XSCALE.',i,'.INP'
        WRITE(*,*)TRIM(name),' was written'
      END DO
      WRITE(*,*)'comment out INPUT_FILEs with large angle, or low strength, or maybe cluster <0 (halo)!'
      STOP
END PROGRAM xscalecluster
