!-----------------------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------------------   
	PROGRAM Quantum_QSPR_Linear_at_Level_N_to_M_Scheme_S
!-----------------------------------------------------------------------------------------------------------
!
!	This program reads a quantum similarity matrix and a properties vector
!	then estimates the properties in a linear QQSPR level N to M calculation,
!	using a procedure Put One In, which is simialr to the Leave One Out.
!	Uses Scheme S to construct the QQSPR operator
!
!-----------------------------------------------------------------------------------------------------------
    Integer*4 MolMolMAX
    parameter(MOLMAX=70)  			!Maximum number of Molecules
    parameter(MolMolMAX=(MOLMAX*MOLMAX+MOLMAX)/2)
!
    CHARACTER (LEN=80):: GENERAL_TITLE, PROPERTY_TITLE
!-------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------
   	 double precision :: zero, one, Toler
     PARAMETER(zero=0.0d0, one=1.0d0, Toler=1.0D-8)
!
	integer*4 :: MoloM,MM,M,MxM					
    Integer*4 :: III,I,U					
    Integer*1 :: Prop_Standard
!
	Double precision :: Z, Prop, E_Prop, EE_Prop, Up_Dw, Mulum, Failures, RealComb, Meantrials
	allocatable:: Z(:), Prop(:), E_Prop(:), Up_Dw(:,:), EE_Prop(:)
!
	Double Precision :: Meanp,  Meanp2, Var, StD   
    Double Precision :: A_Error, Q_Error, Propes
    Double Precision :: P_Mean, EP_Mean, PEP_Cov, P_Var, EP_Var, CorrC
    Double Precision :: AA,BB,EError,AEError
!
!--------------------------------------------------------------------------------------------------
    Character*70 :: ID_Label_Pre,ID_Label_Suf
    Character*70 :: Property_IN
    CHARACTER*70 :: QSMatrix_IN
    Character*70 :: INW_QSMatrix_IN, INW_PRoperty_IN
    CHARACTER*70 :: AUXILY
    Character*200:: OUT_PRoperty
    Character*4  :: Cypher, CNUM
!--------------------------------------------------------------------------------------------------
	Logical*1	::Done
    allocatable	::Done(:)
!--------------------------------------------------------------------------------------------------
!
!		Reads a flag file in order to perform or not
!			optative operations into the data 
!
		open(403,File='Flagstaff_QQSPR.txt')
        read(403,*)Prop_Standard            !Standarizes data or not if /=0 or 0
        read(403,*)Inputo					!Mode of input 0 old  /=0 new
        read(403,*)M						!Number of elements taken from the Molom ones to estimate each Unknown property one
        close(403)
		MxM=(M*M+M)/2
!		
!		reads a label file, where it is found
!		information on the input and output file system
!		as used in the other MQSPS programs
!
		open(304,File='Labelstaff_QQSPR.txt')
!
		Read(304,*)ID_Label_Pre
        Read(304,*)ID_Label_Suf		
        if(Inputo==0)then
		Read(304,*)QSMatrix_IN
       	INW_QSMatrix_IN= TRIM(ID_Label_Pre)//TRIM(QSMatrix_IN)//TRIM(ID_Label_Suf)//'.txt'
        else
        read(304,'(a70)')INW_QSMatrix_IN
        INW_QSMatrix_IN=trim(INW_QSMatrix_IN)
        end if
        Read(304,*)Property_IN
!
        close(304)
!------------------------------------------------------------------------------------------------------------------  
!
!       Reads the optimal similarity matrix 
!
		write(*,'(a70)')INW_QSMatrix_IN
		Open(111,File=INW_QSMatrix_IN)
!
		if(Inputo==0)then
        READ(111,'(a80)')General_title
        READ(111,*)Molom,MM
!
		ALLOCATE(Z(MM))
!
   	 	Read(111,*)(iii,Z(iii),iii=1,MM)
!
		else
!
      	read(111,'(a70)')General_title
	    read(111,*)MoloM
      	read(111,'(a70)')AUXILY
      	read(111,'(a70)')AUXILY
      	read(111,'(a70)')AUXILY
!      
      	MM=(MoloM*MoloM+MoloM)/2
!
		ALLOCATE(Z(MM))
!      
!		retrieve QSM from new input
!
      	Do iii=1,MM
        read(111,*,END=979)molmom,mol,mom,Sobmax
        if(molmom<=MM)then
		Z(iii)=Sobmax
        else
        exit
        end if
        end do
!
		end if        
!        
979        CLOSE(111)
!
		ALLOCATE(Prop(MoloM), E_Prop(MoloM), Up_Dw(MoloM,2), Done(MoloM))  
!
!		Reads the property vector
!             
    	INW_Property_IN=TRIM(ID_Label_Pre)//TRIM(Property_IN)//TRIM(ID_Label_Suf)//'.txt'
		Open(222,File=INW_Property_IN)
        Read(222,'(a80)')Property_title
        Read(222,*)MPM
        if(MPM/=MoloM)then
          write(*,*)'****** Property number and Molecule Number not Coincident : Stop and Correct!'
          pause
          Stop
          end if
!
		Read(222,*)(Prop(iii),iii=1,MoloM)
        close(222) 
!
		if(Prop_Standard/=0)then
        Write(*,*)' Experimental Property Standarized values used'
        P_Mean=Zero
        do I=1,MoloM
          P_Mean=P_mean+Prop(I)
          end do
          P_Mean=P_Mean/MoloM
          P_Var=Zero
          do I=1,MoloM
            P_Var=P_Var+(Prop(I)-P_Mean)**2
            end do
            P_Var=P_Var/MoloM
          do I=1,MoloM
            Prop(I)=(Prop(I)-P_Mean)/sqrt(P_Var)
            end do
!           
        write(*,*)(I,Prop(I),I=1,MoloM)
        end if                 
!
!--------------------------------------------------------------------------------------
!
!*      Start QQSPR procedure
!
        Cypher=CNUM(M)				!Converts the Approximation level into a character to be used as a file label further ahead
		
		if(Prop_Standard==0)then
            OUT_PRoperty = TRIM(ID_Label_Pre)//'Linear_QQSPR_Scheme_S_'//Trim(PROPERTY_IN)//'_No_Standard_'//&
        	TRIM(ID_Label_Suf)//'_Level_'//TRIM(Cypher)//'.txt'
          else
            OUT_PRoperty = TRIM(ID_Label_Pre)//'Linear_QQSPR_Scheme_S_'//Trim(PROPERTY_IN)//'_Standard_'//&
        	TRIM(ID_Label_Suf)//'_Level_'//TRIM(Cypher)//'.txt'
        end if
		Open(313,File=OUT_PRoperty)
!
		write(*,*)' Level of Estimation : ', M	
        Meantrials=zero	
        Takeonein: do U=1,MoloM
        write(*,*)' Estimating element U: ',U
!
		Done(U)=.true.
		Call Order_M( U, M, MxM, Molom, MM, Prop, Z, Meanp, Meanp2, MuluM, Failures, RealComb)
!
        if(RealComb<=zero)then
          write(313,*)'Sorry Estimation Property Failure for property U : ',U
          write(313,*)'Skipping to next item!'
          Done(U)=.false.
          Cycle Takeonein 
          end if
!
		write(*,*)' Mean : ',Meanp
        write(*,*)' # Combinations : ',MuluM
        write(*,*)' # Failures     : ',Failures
        Write(*,*)' # Used Trials  : ',RealComb
        write(*,*)' Mean sq : ',Meanp2
        Meantrials=Meantrials+RealComb
!        
		E_Prop(U)=Meanp/RealComb
        Var=(Meanp2/RealComb)-E_Prop(U)**2
!        
        if(Var>zero)then; StD=SQRT(Var); else; 
          StD=SQRT(Abs(VAR)); Write(*,*)' Variance small or eroneous : ',Var; 
          Write(313,*)' Variance small or erroneous : ',Var; end if
        Up_Dw(U,1)=E_Prop(U)-2*StD
        Up_Dw(U,2)=E_Prop(U)+2*StD
!        
        write(*,*)' Stand. Dv. : ',StD
!
		end do Takeonein !U 
!
!--------------------------------------------------------------------------------------
!               
		A_Error=zero; Q_Error=zero
        Propes=zero
		do I=1,MoloM
        if(Done(I)) then
        A_Error=A_Error+Abs(Prop(I)-E_Prop(I))
        Q_Error=Q_error+(Prop(I)-E_Prop(I))**2
        Propes=Propes+1
        end if
        end do
!
		A_Error=A_Error/Propes; Q_Error=Q_Error/Propes 
!        
        Write(313,*)' * Property estimation by a QQSPR Procedure * '
        Write(313,*)' * Linear Level: ', M,' * '
        Write(313,*)' * Scheme S Computation * '
        If(Prop_Standard/=0)then; write(313,*)' * Property Values Standarized *' 
          				    else; write(313,*)' * Property Values NOT Standarized *'; end if
        Write(313,'(a80)')Property_title
        Write(313,*)' * Initial Number of Molecules : ',MoloM
        Write(313,*)' * Number of Molecules taken IN: ',Propes
        do I=1,MoloM
          if(Done(I))write(313,fmt='(1X,I3,4(1X,E12.5))')I,Prop(I),E_Prop(I),Up_Dw(I,1),Up_Dw(I,2)
          end do
!          
          Write(313,*)' ............................................................................. '
          Write(313,*)' Mean Absolute Error  : ',A_Error
          Write(313,*)' Mean Quadratic Error : ',Q_Error
          Write(313,*)' Sqrt Quadratic Error : ',Sqrt(Q_Error)
!
!		Compute a correlation coefficient between estimates and experimental values
!
		P_Mean=zero; EP_Mean=zero; PEP_Cov=zero; PEP_spro=zero
!        
		do I=1,Molom
        	if(Done(I))then		
        	P_Mean=P_Mean+Prop(I); EP_Mean=EP_Mean+E_Prop(I)
            end if
        	end do
!
        P_Mean=P_Mean/Propes; EP_Mean=EP_Mean/Propes
        P_Var=Zero; EP_Var=zero
        do I=1,MoloM
          if(Done(I))then
          	P_Var=P_Var+(Prop(I)-P_Mean)**2
            EP_Var=EP_Var+(E_Prop(I)-EP_Mean)**2
        	PEP_spro=PEP_spro+(Prop(I)-P_mean)*(E_Prop(I)-EP_Mean)
            end if
            end do
!            
		P_Var=P_Var/Propes
        EP_Var=EP_Var/Propes
        PEP_Cov=PEP_spro/Propes
!
        CorrC=PEP_Cov/sqrt(P_Var*EP_Var)
        AA=PEP_Cov/EP_Var
        BB=P_Mean-EP_Mean*AA
!
		Meantrials=Meantrials/Propes
!
		Write(313,*)' ............................................................................... ' 
        Write(313,*)' * Mean Number of trials       : ', Meantrials      
        Write(313,*)' * Property Mean               : ', P_Mean 
        Write(313,*)' * Property Variance           : ', P_Var
        Write(313,*)' * Estimated Property Mean     : ', EP_Mean
        Write(313,*)' * Estimated Property Variance : ', EP_Var
        Write(313,*)' * Covariance                  : ', PEP_Cov 
		Write(313,*)' ............................................................................... '         
        Write(313,*)' * Correlation Coefficient     : ', CorrC 
        Write(313,*)' * Linear Coefficients: Prop=A*E_Prop+B *' 
        Write(313,*)' * A : ',AA, ' * B : ',BB  
		Write(313,*)' ............................................................................... ' 
!
!		Studies the linear relationship    
!
		Allocate(EE_Prop(MoloM))
		do I=1,MoloM
        EE_Prop(I)=AA*E_Prop(i)+BB
        end do
        EError=zero; AEError=zero
		Do I=1,MoloM 
        DPP=Prop(I)-EE_Prop(i)
        AEError=AEError+Abs(DPP)       
		EError=EError+(DPP)**2
        end do
        AEError=Aeerror/Molom
        EError=EError/MoloM
   		Write(313,*)' ............................................................................... '
        write(313,*)' Mean Absolute estimation error : ', AEError
        write(313,*)' Mean Quadratic estimation error: ', EError
		Write(313,*)' ............................................................................... ' 
        write(313,*)' Molecule   				Property '
        write(313,*)'		              Experimental			Estimated              Abs Error'	
        do I=1,MoloM
        write(313,*)I,Prop(I),EE_Prop(I),abs(Prop(i)-EE_Prop(i)) 
        end do         
		Write(313,*)' ............................................................................... ' 
!
        Close(313)               
!        
!--------------------------------------------------------------------------------------
!
		deallocate(Z,Prop,E_Prop,EE_Prop,Up_DW)
        WRITE(*,*)' Molecular Quantum Similarity Suite'
		write(*,*)' QQSPR Computation Success!'
    	pause
        Stop
!
!--------------------------------------------------------------------------------------------------------------------
	end  Quantum_QSPR_Linear_at_Level_N_to_M_Scheme_S             !program
!...............................................................................................................
!
		subroutine Order_M( U, M, MxM, Molom, MM, Prop, Z, Meanp, Meanp2, Comb, Failures, RealComb) 
!
		Double precision zero,Toler
        parameter (Zero=0.0d0,Toler=1.0D-8)
!        
        logical*1 fail,skip
!
		integer*4 U, M, Molom, MM, I, J, UU, UUU, MxM
        integer*4 INX (M), IU(M), IJ(MxM), IxI, IoI, IoJ
!
        double precision Meanp, Meanp2, Momo
        Double precision Z(MM),Prop(Molom)
        Double Precision Dzi(M,M)
        Double Precision PP(M),PPP(M),ZZ(MxM),ZZI(MxM),Lambda(M) 
        Double Precision Failures,Comb, RealComb
!        
		  Meanp=zero; Meanp2=zero 
          UUU=(U*U-U)/2; UU=U+UUU
!*         
!*----------------------------------------INITIALIZE NESTED DO-------
!*
      DO  I=1,M; Inx(I)=I;  END DO
      Comb=Zero; Failures=Zero; RealComb=Zero
!*
!*----------------------------------------------------- CODE --------
     UNUS: DO 
!*
!*  Use here the combinations of MoloM molecules taken from M by M index vector
!*
        Comb=Comb+1
!
		Skip=.false.
		Scoop: do I=1,M
        	if(Inx(I)==U)then; Skip=.True.
        	Comb=Comb-1
            Exit Scoop
        	end if
        end do Scoop
!                  
        If(.not.Skip) then
		do I=1,M; PP(I)=Prop(Inx(I)); end do
! 
        do I=1,M; IoI=(I*I-I)/2; IxI=(Inx(I)*(Inx(I)-1))/2 
          If(U<Inx(I))then; IU(I)=U+IxI; else; IU(I)=Inx(I)+UUU; end if
          Do J=1,I;  IoJ= J+IoI;
            if(Inx(J)<Inx(I))then; IJ(IoJ)= Inx(J)+ IxI; else; IJ(IoJ)=Inx(I)+(Inx(J)*(Inx(J)-1))/2; end if
            end do !J
          end do !I
!
!		Construct the shifted similarity matrix            
!
		do I=1,M; IoI=(I*I-I)/2; Do J=1,M; 
        If(I>=J) then; IoJ=J+IoI; Else; IoJ=I+(J*J-J)/2; end if
        Dzi(I,J)=Z(IJ(IoJ))-(Z(IU(J))+Z(IU(I))+Z(UU)); end do; end do
!
!	 multiply the Property vector by the transpose and construct a symmetric matrix Z(T)Z to be inverted
!
		do I=1,M
        PPP(I)=zero; Do K=1,M; PPP(I)=PPP(I)+Dzi(K,I)*PP(K); end do
        end do
!
		do I=1,M; IoI=(I*I-I)/2; do J=1,I; IoJ=J+IoI; ZZ(IoJ)=zero
        do K=1,M; ZZ(IoJ)=ZZ(IoJ)+Dzi(K,I)*Dzi(K,J); end do
          end do; end do        
!        
        ZZI=ZZ !Transfer the shifted similarity matrix to the inverse site		
!
!		Invert the shifted similarity matrix via a Gauss-Jordan procedure
!*
		call GJdef2(M,MxM,ZZI,fail)
!*        
        if(.not.fail)then 						!If Inversion fails do not count this combination
        RealComb=RealComb+1
!
!	proceed to compute coefficients and estimate means
!
		Do I=1,M; Lambda(I)=zero; IoI=(I*I-I)/2
        do J=1,M; If(I>=J)then; IoJ=J+IoI; else; IoJ=I+(J*J-J)/2; end if 
          Lambda(I)=Lambda(I)+ZZI(IoJ)*PPP(J)
          end do; end do
!*
!		Estimate unknown value for molecule U
!
		Momo=zero
        do I=1,M
          Momo=Momo+Lambda(I)*(Z(IU(I))-Z(UU))
          end do
!*
        Meanp=Meanp+Momo; Meanp2=Meanp2+Momo**2
!
		else
!        
        Failures=Failures+1
!
        end if          	!.Not.fail
!
		end if 				!.Not.Skip 
!*
!* End of use combination index
!*
!*-----------------------------------NESTED DO LOOP------------------
	SUCUBUS: do  k=M,1,-1; k1=k-1
	TANATHOS: do  Inxk=Inx(k)+1,MoloM
      do i=1,k1; if (Inxk<=Inx(i)) CYCLE TANATHOS; END DO
      Inx(k)=Inxk
!*
      SIXTIUS: do i=k+1,M; i1=i-1
      	 KONDOR: do  Inxi=1,MoloM
         do j=1,i1; if (Inxi<=Inx(j)) CYCLE KONDOR; END DO
      	    EXIT KONDOR
            END DO KONDOR
        IF(Inxi>MoloM)CYCLE TANATHOS
          Inx(i)=Inxi
        END DO SIXTIUS
!*-------------------------------------------------------------------
      CYCLE UNUS
      END DO TANATHOS
      END DO SUCUBUS
      EXIT UNUS
      END DO UNUS
!
!*----------------------------------------------END DO LOOP---------
!*
		end subroutine Order_M
!..................................................................................................
!
!	Gauss-Jordan Inversion procedure Wilkinson-Reinsch 'Linear Algebra' p47
!
	Subroutine GJdef2(n,nn,A,fail)
!
	Implicit none
    integer*4 n,nn
    double precision toler, one
    parameter(toler=1.0d-9, one=1.0d0)
!
!	nn=n*(n+1)/2 is the matrix dimension
!    
    double precision  A(nn), h(n)
    Integer*4 i,ii,ij,k,m
    Double precision p,q
    Double precision r,quot
    Logical*1 fail
!
	fail=.false.
	m=0
    do k=n,1,-1     
      p=A(1); if(p<=toler)then; fail=.true.; return; end if
      ii=1
!
!	use p=a(s) and ii=s if the lower triangle is stored in a(s:s+nn-1)
!
		do i=2,n      
    	m=ii; ii=ii+i  ! m=i*(i-1)/2, ii=i*(i+1)/2
    	q=a(m+1); r=one/p; quot=q*r
    	if(i>k)then; h(i)=quot; else; h(i)=-quot; end if
      	do ij=m+2,ii;  a(ij-i)=a(ij)+q*h(ij-m); end do
        end do !i
!
	m=m-1 !m=n*(n-1)/2-1 , ii=n*(n+1)/2        
    A(ii)=r
    do i=2,n; A(m+i)=h(i); end do
!
	end do !k
!
	end Subroutine GJdef2          
!-------------------------------------------------------------------------------------