
      program main
      implicit none
      double precision x(1000),y(1000),dx(1000),dy(1000),
     *  slope,fintcpt,rtmnsquare,xoutliers(100),youtliers(1000),
     *  x1(1000),y1(1000),k,b,sum1,sum2
      integer nsamp,i,numoutliers,nsamp1

      open(unit=1,file='testdata.txt')
      i=1
10    read(1,*,end=100)x(i),y(i)
      i=i+1
      goto 10
100   nsamp=i-1
                  
      goto 200
      
      
      nsamp=13
      x(1)=-1.0d0
      y(1)=3.0d0
      x(2)=-3.0d0
      y(2)=4.0d0
      x(3)=-5.0d0
      y(3)=5.0d0
      x(4)=-7.0d0
      y(4)=6.0d0
      x(5)=-9.0d0
      y(5)=7.0d0
      x(6)=-11.0d0
      y(6)=8.0d0
      x(7)=-3.0d0
      y(7)=1.0d0
      x(8)=-5.0d0
      y(8)=2.0d0
      x(9)=-7.0d0
      y(9)=3.0d0
      x(10)=-9.0d0
      y(10)=4.0d0
      x(11)=-11.0d0
      y(11)=5.0d0
      x(12)=-4.0d0
      y(12)=7.0d0
      x(13)=-12.0d0
      y(13)=1.0d0

200   slope=-2.0d0
      fintcpt=0.0d0 
     
      call OrthSoilRespRegres(nsamp,x,y,slope,fintcpt)
      write(*,*)slope,fintcpt
      pause
      
      slope=-2.0d0
      fintcpt=0.0d0 
      
      call orthlinreg_outlier(nsamp,x,y,slope,
     &  fintcpt,dx,dy,rtmnsquare,xoutliers,youtliers,
     &  numoutliers)
       write(*,*)slope/2.0d0,fintcpt,numoutliers
       do i=1,numoutliers
         write(*,*)xoutliers(i),youtliers(i)
       enddo
       end
     
      subroutine orthlinreg_outlier(nsamp0,x0,y0,slope,
     &  fintcpt,dx,dy,rtmnsquare,xoutliers,youtliers,
     &  numoutliers)
      implicit none
      integer nsamp0,numoutliers
      double precision x0(nsamp0),y0(nsamp0),slope,
     &  fintcpt,dx(nsamp0),dy(nsamp0),rtmnsquare,
     &  xoutliers(nsamp0),youtliers(nsamp0),xtest(nsamp0),
     &  ytest(nsamp0),slopetest,fintcpttest,dxtest(nsamp0),
     &  dytest(nsamp0),rtmnsquaretest,testmeasure(nsamp0),
     &  x(nsamp0),y(nsamp0)
      
      integer iwhichside,nsamptest,isitoutlier,
     &   isoutlier_1side,i,j,nsamp
      parameter (iwhichside=1)
      
      numoutliers=0
      nsamp=nsamp0
      do i=1,nsamp
        x(i)=x0(i)
        y(i)=y0(i)
      enddo
      
50    call orthlinreg(nsamp,x,y,slope,fintcpt,
     &      dx,dy,rtmnsquare)
     
      write(*,*)slope,fintcpt,rtmnsquare
      stop
      nsamptest=nsamp-1
      do i=1,nsamp
        do j=1,nsamp
          xtest(j)=x(j)
          ytest(j)=y(j)
        enddo
        xtest(i)=x(nsamp)
        ytest(i)=y(nsamp)
        call orthlinreg(nsamptest,xtest,ytest,slopetest,
     &    fintcpttest,dxtest,dytest,rtmnsquaretest)
!        write(*,*)i,slopetest,fintcpttest
     
!        testmeasure(i)=(slopetest-slope)**2+
!     &    (fintcpttest-fintcpt)**2
        
        testmeasure(i)=100.0d0*dabs(rtmnsquaretest-rtmnsquare)/
     &   rtmnsquare
!        write(*,*)i,testmeasure(i)
      enddo
      
      isitoutlier=isoutlier_1side(nsamp,testmeasure,iwhichside)
      if(isitoutlier.lt.1.or.isitoutlier.gt.nsamp)return
! outlier detected
      numoutliers=numoutliers+1
      xoutliers(numoutliers)=x(isitoutlier)
      youtliers(numoutliers)=y(isitoutlier)
      x(isitoutlier)=x(nsamp)
      y(isitoutlier)=y(nsamp)
      nsamp=nsamp-1
      if(nsamp.le.2)then
        write(*,*)'No enough good data left'
        stop
      endif
      goto 50
      return
      end 
          
! orthogonal linear regression
      subroutine orthlinreg(nsamp,x,y,slope0,fintcpt0,
     &      dx,dy,rtmnsquare)
      implicit none
      integer nsamp
      double precision x(nsamp),y(nsamp),dx1(nsamp),
     &    dy1(nsamp),slope(2),fintcpt(2),dx2(nsamp),
     &    dy2(nsamp),slope0,fintcpt0,dx(nsamp),dy(nsamp)
      integer i,j 
      double precision w,u,v,xbar,ybar,root1,root2,
     & a,b,c,rtmnsquare1,rtmnsquare2,rtmnsquare
      

      
      xbar=0.0d0
      ybar=0.0d0
      w=0.0d0
      u=0.0d0
      v=0.0d0
      do i=1,nsamp
        xbar=xbar+x(i)
        ybar=ybar+y(i)
        w=w+x(i)*x(i)
        u=u+y(i)*y(i)
        v=v+x(i)*y(i)
      enddo
      xbar=xbar/dble(nsamp)
      ybar=ybar/dble(nsamp)
      w=w/dble(nsamp)
      u=u/dble(nsamp)
      v=v/dble(nsamp)
      a=v-xbar*ybar
      b=w-u-xbar*xbar+ybar*ybar
      c=xbar*ybar-v
      call quadraticroots(a,b,c,root1,root2)
      slope(1)=root1
      slope(2)=root2
      fintcpt(1)=ybar-slope(1)*xbar
      fintcpt(2)=ybar-slope(2)*xbar
      rtmnsquare1=0.0d0
      rtmnsquare2=0.0d0
      do i=1,nsamp
        dx1(i)=(y(i)-fintcpt(1)-x(i)*slope(1))*
     &    slope(1)/(1.0d0+slope(1)*slope(1))
        dy1(i)=-(y(i)-fintcpt(1)-x(i)*slope(1))/
     &          (1.0d0+slope(1)*slope(1))     
        rtmnsquare1=rtmnsquare1+dx1(i)**2+dy1(i)**2
                
        dx2(i)=(y(i)-fintcpt(2)-x(i)*slope(2))*
     &    slope(2)/(1.0d0+slope(2)*slope(2))
        dy2(i)=-(y(i)-fintcpt(2)-x(i)*slope(2))/
     &          (1.0d0+slope(2)*slope(2))
        rtmnsquare2=rtmnsquare2+dx2(i)**2+dy2(i)**2
      enddo
      rtmnsquare1=dsqrt(rtmnsquare1/dble(nsamp))
      rtmnsquare2=dsqrt(rtmnsquare2/dble(nsamp))
      if(rtmnsquare1.gt.rtmnsquare2)then
        rtmnsquare=rtmnsquare2
        slope0=slope(2)
        fintcpt0=fintcpt(2)
        do i=1,nsamp
          dx(i)=dx2(i)
          dy(i)=dy2(i)
        enddo
      else
        rtmnsquare=rtmnsquare1
        slope0=slope(1)
        fintcpt0=fintcpt(1)
        do i=1,nsamp
          dx(i)=dx1(i)
          dy(i)=dy1(i)
        enddo
      endif
      return
      end
      
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!
      subroutine OrthSoilRespRegres(npoints,x0,y0,slope,fintcpt)
      implicit none
c
C  ODRPACK ARGUMENT DEFINITIONS
C      ==> FCN      NAME OF THE USER SUPPLIED FUNCTION SUBROUTINE
C      ==> N        NUMBER OF OBSERVATIONS 
C      ==> M        COLUMNS OF DATA IN THE EXPLANATORY VARIABLE
C      ==> NP       NUMBER OF PARAMETERS
C      ==> NQ       NUMBER OF RESPONSES PER OBSERVATION
C     <==> BETA     FUNCTION PARAMETERS
C      ==> Y        RESPONSE VARIABLE
C      ==> LDY      LEADING DIMENSION OF ARRAY Y
C      ==> X        EXPLANATORY VARIABLE
C      ==> LDX      LEADING DIMENSION OF ARRAY X
C      ==> WE       "EPSILON" WEIGHTS
C      ==> LDWE     LEADING DIMENSION OF ARRAY WE
C      ==> LD2WE    SECOND DIMENSION OF ARRAY WE
C      ==> WD       "DELTA" WEIGHTS
C      ==> LDWD     LEADING DIMENSION OF ARRAY WD
C      ==> LD2WD    SECOND DIMENSION OF ARRAY WD
C      ==> IFIXB    INDICATORS FOR "FIXING" PARAMETERS (BETA)
C      ==> IFIXX    INDICATORS FOR "FIXING" EXPLANATORY VARIABLE (X)
C      ==> LDIFX    LEADING DIMENSION OF ARRAY IFIXX
C      ==> JOB      TASK TO BE PERFORMED 
C      ==> NDIGIT   GOOD DIGITS IN SUBROUTINE FUNCTION RESULTS
C      ==> TAUFAC   TRUST REGION INITIALIZATION FACTOR
C      ==> SSTOL    SUM OF SQUARES CONVERGENCE CRITERION
C      ==> PARTOL   PARAMETER CONVERGENCE CRITERION
C      ==> MAXIT    MAXIMUM NUMBER OF ITERATIONS
C      ==> IPRINT   PRINT CONTROL 
C      ==> LUNERR   LOGICAL UNIT FOR ERROR REPORTS 
C      ==> LUNRPT   LOGICAL UNIT FOR COMPUTATION REPORTS 
C      ==> STPB     STEP SIZES FOR FINITE DIFFERENCE DERIVATIVES WRT BETA
C      ==> STPD     STEP SIZES FOR FINITE DIFFERENCE DERIVATIVES WRT DELTA
C      ==> LDSTPD   LEADING DIMENSION OF ARRAY STPD
C      ==> SCLB     SCALE VALUES FOR PARAMETERS BETA
C      ==> SCLD     SCALE VALUES FOR ERRORS DELTA IN EXPLANATORY VARIABLE 
C      ==> LDSCLD   LEADING DIMENSION OF ARRAY SCLD
C     <==> WORK     DOUBLE PRECISION WORK VECTOR
C      ==> LWORK    DIMENSION OF VECTOR WORK
C     <==  IWORK    INTEGER WORK VECTOR
C      ==> LIWORK   DIMENSION OF VECTOR IWORK
C     <==  INFO     STOPPING CONDITION 
 
C  PARAMETERS SPECIFYING MAXIMUM PROBLEM SIZES HANDLED BY THIS DRIVER
C     MAXN          MAXIMUM NUMBER OF OBSERVATIONS 
C     MAXM          MAXIMUM NUMBER OF COLUMNS IN EXPLANATORY VARIABLE
C     MAXNP         MAXIMUM NUMBER OF FUNCTION PARAMETERS
C     MAXNQ         MAXIMUM NUMBER OF RESPONSES PER OBSERVATION

C  PARAMETER DECLARATIONS AND SPECIFICATIONS
      INTEGER    LDIFX,LDSCLD,LDSTPD,LDWD,LDWE,LDX,LDY,LD2WD,LD2WE,
     +           LIWORK,LWORK,MAXM,MAXN,MAXNP,MAXNQ
      PARAMETER (MAXM=25,MAXN=50000,MAXNP=30,MAXNQ=1,
     +           LDY=MAXN,LDX=MAXN,
     +           LDWE=1,LD2WE=1,LDWD=1,LD2WD=1,
     +           LDIFX=MAXN,LDSTPD=1,LDSCLD=1,
     +           LWORK=18 + 11*MAXNP + MAXNP**2 + MAXM + MAXM**2 + 
     +                 4*MAXN*MAXNQ + 6*MAXN*MAXM + 2*MAXN*MAXNQ*MAXNP +  
     +                 2*MAXN*MAXNQ*MAXM + MAXNQ**2 + 
     +                 5*MAXNQ + MAXNQ*(MAXNP+MAXM) + LDWE*LD2WE*MAXNQ,
     +           LIWORK=20+MAXNP+MAXNQ*(MAXNP+MAXM))
C  VARIABLE DECLARATIONS 
      INTEGER          I,INFO,IPRINT,J,JOB,L,LUNERR,LUNRPT,M,MAXIT,N,
     +                 NDIGIT,NP,NQ
      INTEGER          IFIXB(MAXNP),IFIXX(LDIFX,MAXM),IWORK(LIWORK)
      DOUBLE PRECISION PARTOL,SSTOL,TAUFAC
      DOUBLE PRECISION BETA(MAXNP),SCLB(MAXNP),SCLD(LDSCLD,MAXM),
     +                 STPB(MAXNP),STPD(LDSTPD,MAXM),
     +                 WD(LDWD,LD2WD,MAXM),WE(LDWE,LD2WE,MAXNQ),
     +                 WORK(LWORK),X(LDX,MAXM),Y(LDY,MAXNQ)
c
      integer npoints,i1
      double precision x0(npoints),y0(npoints),slope,fintcpt
      
    
      
      EXTERNAL         OrthRespFCN
c
C  SPECIFY DEFAULT VALUES FOR DODRC ARGUMENTS
      WE(1,1,1)  = -1.0D0
      WD(1,1,1)  = -1.0D0
      IFIXB(1)   = -1
!      IFIXX(1,1) = -1
!      JOB        = 00023
      JOB=20
      NDIGIT     = -1
      TAUFAC     = -1.0D0
      SSTOL      = -1.0D0
      PARTOL     = -1.0D0
      MAXIT      = -1
!      IPRINT     = -1
!      IPRINT=0
      IPRINT=-1
      LUNERR     = -1
      LUNRPT     = -1
      STPB(1)    = -1.0D0
      STPD(1,1)  = -1.0D0
      SCLB(1)    = -1.0D0
      SCLD(1,1)  = -1.0D0

      MAXIT      = 200000
C  SET UP ODRPACK REPORT FILES
      LUNERR  =   9
      LUNRPT  =   9     
c
      N=npoints
      M=1
      NP=2
      NQ=1

            
      do I=1,N
        do i1=1,M
          X(I,i1)=x0(I)
        enddo
        Y(I,1)=y0(I)
      enddo
      BETA(1)=slope
      BETA(2)=fintcpt
      
C  READ PROBLEM DATA, AND SET NONDEFAULT VALUE FOR ARGUMENT IFIXX
      DO 10 I=1,N
         DO 15 J=1, M
           IFIXX(I,J) = 1
15       CONTINUE
10    CONTINUE
60        CALL DODRC(OrthRespFCN,
     +           N,M,NP,NQ,
     +           BETA,
     +           Y,LDY,X,LDX,
     +           WE,LDWE,LD2WE,WD,LDWD,LD2WD,
     +           IFIXB,IFIXX,LDIFX,
     +           JOB,NDIGIT,TAUFAC,
     +           SSTOL,PARTOL,MAXIT,
     +           IPRINT,LUNERR,LUNRPT,
     +           STPB,STPD,LDSTPD,
     +           SCLB,SCLD,LDSCLD,
     +           WORK,LWORK,IWORK,LIWORK,
     +           INFO)

      slope=BETA(1)
      fintcpt=BETA(2)
      return
      END
c
      SUBROUTINE OrthRespFCN(N,M,NP,NQ,
     +               LDN,LDM,LDNP,
     +               BETA,XPLUSD,
     +               IFIXB,IFIXX,LDIFX,
     +               IDEVAL,F,FJACB,FJACD,
     +               ISTOP)
      implicit none      
C  SUBROUTINE ARGUMENTS
C      ==> N        NUMBER OF OBSERVATIONS
C      ==> M        NUMBER OF COLUMNS IN EXPLANATORY VARIABLE
C      ==> NP       NUMBER OF PARAMETERS
C      ==> NQ       NUMBER OF RESPONSES PER OBSERVATION
C      ==> LDN      LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING N
C      ==> LDM      LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING M
C      ==> LDNP     LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING NP
C      ==> BETA     CURRENT VALUES OF PARAMETERS
C      ==> XPLUSD   CURRENT VALUE OF EXPLANATORY VARIABLE, I.E., X + DELTA
C      ==> IFIXB    INDICATORS FOR "FIXING" PARAMETERS (BETA)
C      ==> IFIXX    INDICATORS FOR "FIXING" EXPLANATORY VARIABLE (X)
C      ==> LDIFX    LEADING DIMENSION OF ARRAY IFIXX
C      ==> IDEVAL   INDICATOR FOR SELECTING COMPUTATION TO BE PERFORMED
C     <==  F        PREDICTED FUNCTION VALUES
C     <==  FJACB    JACOBIAN WITH RESPECT TO BETA
C     <==  FJACD    JACOBIAN WITH RESPECT TO ERRORS DELTA
C     <==  ISTOP    STOPPING CONDITION, WHERE
C                     0 MEANS CURRENT BETA AND X+DELTA WERE
C                       ACCEPTABLE AND VALUES WERE COMPUTED SUCCESSFULLY
C                     1 MEANS CURRENT BETA AND X+DELTA ARE
C                       NOT ACCEPTABLE;  ODRPACK SHOULD SELECT VALUES 
C                       CLOSER TO MOST RECENTLY USED VALUES IF POSSIBLE
C                    -1 MEANS CURRENT BETA AND X+DELTA ARE
C                       NOT ACCEPTABLE;  ODRPACK SHOULD STOP

C  INPUT ARGUMENTS, NOT TO BE CHANGED BY THIS ROUTINE:
      INTEGER          I,IDEVAL,ISTOP,L,LDIFX,LDM,LDN,LDNP,M,N,NP,NQ
      DOUBLE PRECISION BETA(NP),XPLUSD(LDN,M)
      INTEGER          IFIXB(NP),IFIXX(LDIFX,M)
C  OUTPUT ARGUMENTS:
      DOUBLE PRECISION F(LDN,NQ),FJACB(LDN,LDNP,NQ),FJACD(LDN,LDM,NQ)
            
C  CHECK FOR UNACCEPTABLE VALUES FOR THIS PROBLEM
c
!      
      IF (MOD(IDEVAL,10).GE.1) THEN
        DO 110 L = 1,NQ
          DO 100 I = 1,N
            F(I,L)=BETA(2)+BETA(1)*XPLUSD(I,1)
  100     CONTINUE
  110   CONTINUE
      END IF

C  COMPUTE DERIVATIVES WITH RESPECT TO BETA
      IF (MOD(IDEVAL/10,10).GE.1) THEN
        DO 210 L = 1,NQ
          DO 200 I = 1,N
            FJACB(I,1,L)=XPLUSD(I,1)
            FJACB(I,2,L)=1.0d0
  200     CONTINUE
  210   CONTINUE
      ENDIF
      
C  COMPUTE DERIVATIVES WITH RESPECT TO DELTA
      IF (MOD(IDEVAL/100,10).GE.1) THEN
        DO 310 L = 1,NQ
          DO 300 I = 1,N
            FJACD(I,1,L)=BETA(1)
  300     CONTINUE
  310   CONTINUE
          
      END IF
      RETURN
      END
!
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
