Subroutine GenericRegres(npoints,ny,y,nx,x,weity0,
     &weitx0,ndim0,beta_in_out,betamin0,betamax0,xmin0,xmax0,
     &iderivative,iregrestype0,shorty0,shortx0,fatbeta)
      implicit none
!iregrestype0=0, ordinary regression
!iregrestype0=1, orthogonal distance regression. Direct search methods
!  determine the shortest distance within the iteration
!iregrestype0=2, orthogonal distance regression. Direct search methods
!  expand the parameter vector to include x positions.
!iregrestype0=-1, implicit regression
!iderivative=0, no derivative provided
!iderivative=1, derivative provided
      include 'forgenericregres.h'
      integer npoints,ny,nx,iderivative,ndim0,iregrestype0
      double precision y(npoints,ny),x(npoints,nx),weity0(npoints,ny),
     &weitx0(npoints,nx),xmin0(npoints,nx),xmax0(npoints,nx),
     &beta_in_out(ndim0),betamin0(ndim0),betamax0(ndim0),
     &shorty0(npoints,ny),shortx0(npoints,nx),fatbeta
!
      integer i,j,INFO,ndim,k,n,i2,icompete,isitnaninf,nave
      double precision xtol,beta(ndim0+nx*npoints),
     &betacp(ndim0+nx*npoints),fatbetacp,beta0(ndim0+nx*npoints),
     &fatbeta0,ftol,gacontrol(12),ran2,ftol_relax,term1,term2,discount,
     &history(2000,ndim0+3),upper,lower,f1dim_generic,generic_pikaia
      parameter(xtol=1.0d-7,ftol=1.0d-7)
      external funkmin_generic,FCN_generic,f1dim_generic,generic_pikaia
!-----------------------------------------------------
      ndim=ndim0
      nxvars=nx
      nyvars=ny
      if((nx*npoints+ndim0).gt.1000)iregrestype0=0
      iregrestype=iregrestype0
      iknowder=iderivative
      nobs=npoints
      do i=1,npoints
        do j=1,nxvars
          xvars(i,j)=x(i,j)
          xmin(i,j)=xmin0(i,j)
          xmax(i,j)=xmax0(i,j)
          weitx0(i,j)=1.0d0
          weitx(i,j)=weitx0(i,j)
        enddo
        do j=1,nyvars
          yobs(i,j)=y(i,j)
          weity(i,j)=weity0(i,j)
        enddo
      enddo
      do i=1,ndim
        betamin(i)=betamin0(i)
        betamax(i)=betamax0(i)
        beta(i)=beta_in_out(i)
      enddo
      if(iregrestype.eq.2)iregrestype=1
c     gacontrol( 1) - number of individuals in a population (default
c                      is 100)
c     gacontrol( 2) - number of generations over which solution is
c                      to evolve (default is 500)
c     gacontrol( 3) - number of significant digits (i.e., number of
c                      genes) retained in chromosomal encoding (default
c                      is 6)  (Note: This number is limited by the
c                      machine floating point precision.  Most 32-bit
c                      floating point representations have only 6 full
c                      digits of precision.  To achieve greater preci-
c                      sion this routine could be converted to double
c                      precision, but note that this would also require
c                      a double precision random number generator, which
c                      likely would not have more than 9 digits of
c                      precision if it used 4-byte integers internally.)
c     gacontrol( 4) - crossover probability; must be  <= 1.0 (default
c                      is 0.85). If crossover takes place, either one
c                      or two splicing points are used, with equal
c                      probabilities
c     gacontrol( 5) - mutation mode; 1/2/3/4/5 (default is 2)
c                      1=one-point mutation, fixed rate
c                      2=one-point, adjustable rate based on fitness 
c                      3=one-point, adjustable rate based on distance
c                      4=one-point+creep, fixed rate
c                      5=one-point+creep, adjustable rate based on fitness
c                      6=one-point+creep, adjustable rate based on distance
c     gacontrol( 6) - initial mutation rate; should be small (default
c                      is 0.005) (Note: the mutation rate is the proba-
c                      bility that any one gene locus will mutate in
c                      any one generation.)
c     gacontrol( 7) - minimum mutation rate; must be >= 0.0 (default
c                      is 0.0005)
c     gacontrol( 8) - maximum mutation rate; must be <= 1.0 (default
c                      is 0.25)
c     gacontrol( 9) - relative fitness differential; range from 0
c                      (none) to 1 (maximum).  (default is 1.)
c     gacontrol(10) - reproduction plan; 1/2/3=Full generational
c                      replacement/Steady-state-replace-random/Steady-
c                      state-replace-worst (default is 3)
c     gacontrol(11) - elitism flag; 0/1=off/on (default is 0)
c                      (Applies only to reproduction plans 1 and 2)
c     gacontrol(12) - printed output 0/1/2=None/Minimal/Verbose
c                      (default is 0)
      idobounded=1
10    call funkmin_generic(ndim,beta,fatbeta)
      do i=1,ndim
        beta0(i)=beta(i)
        history(1,i)=beta(i)
      enddo
      fatbeta0=fatbeta
      history(1,ndim+1)=fatbeta
!entrance counter
      history(1,ndim+2)=1.0d0
!failure counter
      history(1,ndim+3)=0.0d0
!Is it a competition among different initial guesses?
      icompete=0
!j the total number of calls to nongradopt; k is the number of returns to the current best and reset
!to zero if a better minumum is found; n is the number of scouting points over the landscape of the cost function.
!The first initial guess provided by the user is always part of the set of scouting points.the rest consist of outcomes
!from calls to nongradopt if they are significantly different from the current best.
      j=0
      k=0
      n=1
      nave=1
      ftol_relax=ftol*1000.0d0
      discount=2.0d0
!relax the convergence criterion for scouting
30    fatbetacp=fatbeta
      do i=1,ndim
        betacp(i)=beta(i)
      enddo
      INFO=iregrestype
      call odr_leastsquare(ndim,FCN_generic,beta,nobs,
     &xvars(1:nobs,1:nxvars),nxvars,yobs(1:nobs,1:nyvars),
     &nyvars,weitx(1:nobs,1:nxvars),weity(1:nobs,1:nyvars),
     &iderivative,shortx(1:nobs,1:nxvars),
     &shorty(1:nobs,1:nyvars),fatbeta,INFO)
      call funkmin_generic(ndim,beta,fatbeta)
      if(isitnaninf(fatbeta).eq.1.or.fatbeta.gt.fatbetacp)then
        fatbeta=fatbetacp
        do i=1,ndim
          beta(i)=betacp(i)
        enddo
      else
        fatbetacp=fatbeta
        do i=1,ndim
          betacp(i)=beta(i)
        enddo
      endif
      call nongradopt(ndim,funkmin_generic,f1dim_generic,
     &beta,betamin,betamax,ftol_relax,fatbeta)
      if(isitnaninf(fatbeta).eq.1.or.fatbeta.gt.fatbetacp)then
        fatbeta=fatbetacp
        do i=1,ndim
          beta(i)=betacp(i)
        enddo
      endif
      if(fatbeta.gt.1.0d0)then
        term1=fatbeta*ftol_relax
      else
        term1=ftol_relax*10.0d0
      endif
      if(fatbeta.gt.fatbeta0)then
!failure
        if((fatbeta-fatbeta0).gt.term1)then
          if(icompete.eq.1)history(1,ndim+3)=history(1,ndim+3)+1.5d0
!even though fatbeta is much worse than fatbeta0, it is an output of optimization after all so
!include it in the set if it has not already been included in the set.
          i=1
          i2=1
40        if(dabs(history(i2,i)-beta(i)).gt.ftol_relax)then
            if(dabs(history(i2,ndim+1)-fatbeta).lt.term1)then
              history(i2,ndim+3)=history(i2,ndim+3)+1.0d0
              goto 60
            endif
            if(i2.ge.n)goto 50
            i2=i2+1
            i=1
            goto 40
          else
            if(i.ge.ndim)goto 60 
            i=i+1
            goto 40
          endif
50        n=n+1
          do i=1,ndim
            history(n,i)=beta(i)
          enddo
          history(n,ndim+1)=fatbeta
          history(n,ndim+2)=0.0d0
          history(n,ndim+3)=0.0d0
!use average only when there is imporvement
          nave=n
        else
!the difference is minimal even though fatbeta is larger than fatbeta0.
!Increment the counter for arriving at the same minimum.
          if(icompete.eq.1)history(1,ndim+3)=history(1,ndim+3)+1.0d0
          k=k+1
        endif
60      do i=1,ndim
          beta(i)=beta0(i)
        enddo
        fatbeta=fatbeta0
      else
!success
        if((fatbeta0-fatbeta).lt.term1)then
!negligible improvement. Increment the counter for arriving at the same minimum.
!no increment for the set of central initial guesses
          if(icompete.eq.1)history(1,ndim+3)=history(1,ndim+3)+0.5d0
          k=k+1
        else
!reset the counter for arriving at a better minimum.
!Increment the set of central initial guesses
          if(dabs(discount-2.0d0).lt.ftol)then
            discount=dmax1(0.001d0,(fatbeta0-fatbeta)/1000.0d0)
          endif
          k=0
          n=n+1
          do i=1,ndim
            history(n,i)=beta(i)
          enddo
          history(n,ndim+1)=fatbeta
          history(n,ndim+2)=0.0d0
          history(n,ndim+3)=0.0d0
        endif
        do i=1,ndim
          beta0(i)=beta(i)
        enddo
        fatbeta0=fatbeta
      endif
      j=j+1
      if(j.lt.20.and.k.lt.2)then
        if(j.lt.10)then
          term1=0.01d0+dmin1(history(1,ndim+3)*0.025d0,0.9d0)
          history(1,ndim+2)=history(1,ndim+2)+1.0d0
          do i=1,ndim
            lower=history(1,i)-term1*(history(1,i)-betamin(i))
            upper=history(1,i)+term1*(betamax(i)-history(1,i))
            beta(i)=lower+ran2()*(upper-lower)
          enddo
          icompete=1
          goto 70
        endif
!try average
        if(n.gt.nave)then
          term1=1.0d0/(history(1,ndim+1)+1.0d-5)
          do i=2,n
            term1=term1+1.0d0/(history(i,ndim+1)+1.0d-5)
          enddo
          do i=1,ndim
            beta(i)=history(1,i)/(term1*(history(1,ndim+1)+1.0d-5))
            do icompete=2,n
              beta(i)=beta(i)+history(icompete,i)/
     &(term1*(history(icompete,ndim+1)+1.0d-5))
            enddo
          enddo
          nave=n
          icompete=0
          goto 70
        endif
!try different initial guesses
        if(ran2().gt.0.2d0)then
!guess around the best
          icompete=1
          term1=history(1,ndim+1)+
     &discount*history(1,ndim+2)*history(1,ndim+3)
          do i=2,n
            term2=history(i,ndim+1)+
     &discount*history(i,ndim+2)*history(i,ndim+3)
            if(term2.le.term1)then
              term1=term2
              do i2=1,ndim+3
                history(n+1,i2)=history(i,i2)
                history(i,i2)=history(1,i2)
                history(1,i2)=history(n+1,i2)
              enddo
            endif
          enddo
          term1=0.01d0+dmin1(history(1,ndim+2)*history(1,ndim+3)*
     &0.015d0,0.9d0)
          history(1,ndim+2)=history(1,ndim+2)+1.0d0
          do i=1,ndim
            lower=history(1,i)-term1*(history(1,i)-betamin(i))
            upper=history(1,i)+term1*(betamax(i)-history(1,i))
            beta(i)=lower+ran2()*(upper-lower)
          enddo
        else
!completely random guess
          do i=1,ndim
            beta(i)=betamin(i)+ran2()*(betamax(i)-betamin(i))
          enddo
          icompete=0
        endif
70      call funkmin_generic(ndim,beta,fatbeta)
        goto 30
      else
        if((ftol_relax-ftol).gt.ftol)then
          if(k.le.1)then
            n=n+1
            do i=1,ndim+3
              history(n,i)=history(1,i)
            enddo
            do i=1,ndim
              history(1,i)=beta(i)
            enddo
            history(1,ndim+1)=fatbeta
            history(1,ndim+2)=0.0d0
            history(1,ndim+3)=0.0d0
            do i=1,n
              do icompete=1,ndim
                betacp(icompete)=history(i,icompete)
              enddo
              fatbetacp=history(i,ndim+1)
              call RepeatCompassSearch(ndim,betacp,fatbetacp,
     &betamin,betamax,funkmin_generic,f1dim_generic,ftol_relax)
              call funkmin_generic(ndim,betacp,fatbetacp)
              if(isitnaninf(fatbetacp).eq.0.and.fatbetacp.lt.
     &fatbeta)then
                do icompete=1,ndim
                  beta(icompete)=betacp(icompete)
                enddo
                fatbeta=fatbetacp
              endif
            enddo
            do i=1,ndim
              beta0(i)=beta(i)
            enddo
            fatbeta0=fatbeta
            icompete=1
            j=0
          else
            icompete=0
          endif
          ftol_relax=ftol
          goto 30
        endif
      endif

      goto 110

      call RepeatCompassSearch(ndim,beta,fatbeta,
     &betamin,betamax,funkmin_generic,f1dim_generic,xtol)
      call funkmin_generic(ndim,beta,fatbeta)
      if(isitnaninf(fatbeta).eq.1.or.fatbeta.ge.fatbeta0)then
!if RepeatCompassSearch cannot improve, we end the search
        do i=1,ndim
          beta(i)=beta0(i)
        enddo
        fatbeta=fatbeta0
        goto 110
      endif
      do i=1,12
        gacontrol(i)=-1.0d0
      enddo
      gacontrol(1)=250.0d0
      gacontrol(2)=5000.0d0
      gacontrol(3)=8.0d0
      do i=1,ndim
        beta0(i)=(beta(i)-betamin(i))/(betamax(i)-betamin(i))
      enddo
      fatbeta0=fatbeta
      call pikaia(generic_pikaia,ndim,gacontrol,beta0,fatbeta0,j)
      fatbeta0=1.0d+100
      if(j.eq.0)then
        do i=1,ndim
          beta0(i)=betamin(i)+beta0(i)*(betamax(i)-betamin(i))
        enddo
        call funkmin_generic(ndim,beta0,fatbeta0)
      endif
80    if(isitnaninf(fatbeta0).eq.1.or.fatbeta0.gt.fatbeta)then
        fatbeta0=fatbeta
        do i=1,ndim
          beta0(i)=beta(i)
        enddo
      else
        do i=1,ndim
          beta(i)=beta0(i)
        enddo
        fatbeta=fatbeta0
      endif
!
      INFO=iregrestype
      call odr_leastsquare(ndim,FCN_generic,beta,nobs,
     &xvars(1:nobs,1:nxvars),nxvars,yobs(1:nobs,1:nyvars),
     &nyvars,weitx(1:nobs,1:nxvars),weity(1:nobs,1:nyvars),
     &iderivative,shortx(1:nobs,1:nxvars),
     &shorty(1:nobs,1:nyvars),fatbeta,INFO)
      call funkmin_generic(ndim,beta,fatbeta)
      if(isitnaninf(fatbeta).eq.1.or.fatbeta.gt.fatbeta0)then
        do i=1,ndim
          beta(i)=beta0(i)
        enddo
        fatbeta=fatbeta0
      endif
      iregrestype=iregrestype0
      if(iregrestype.eq.2)then
        do i=1,npoints
          do j=1,nx
            ndim=ndim+1
            beta(ndim)=shortx(i,j)
            betamin(ndim)=xmin0(i,j)
            betamax(ndim)=xmax0(i,j)
            if(beta(ndim).lt.betamin(ndim).or.
     &beta(ndim).gt.betamax(ndim))then
              beta(ndim)=x(i,j)
            endif
          enddo
        enddo
        call funkmin_generic(ndim,beta,fatbeta)
      endif
      j=0
100   j=j+1
      fatbeta0=fatbeta
      do i=1,ndim
        beta0(i)=beta(i)
      enddo
      call nongradopt(ndim,funkmin_generic,
     &f1dim_generic,beta,betamin,betamax,ftol,fatbeta)
      call funkmin_generic(ndim,beta,fatbeta)
      if(isitnaninf(fatbeta).eq.1.or.fatbeta.ge.fatbeta0)then
        fatbeta=fatbeta0
        do i=1,ndim
          beta(i)=beta0(i)
        enddo
        goto 110
      endif
      fatbetacp=fatbeta
      do i=1,ndim
        betacp(i)=beta(i)
      enddo
      call RepeatCompassSearch(ndim,betacp,fatbetacp,
     &betamin,betamax,funkmin_generic,f1dim_generic,xtol)
      call funkmin_generic(ndim,betacp,fatbetacp)
      if(isitnaninf(fatbetacp).eq.1.or.fatbetacp.ge.fatbeta)then
        goto 110
      else
        fatbeta=fatbetacp
        do i=1,ndim
          beta(i)=betacp(i)
        enddo
      endif
      if(j.ge.2.or.fatbeta.eq.fatbeta0)goto 110
      if(dabs(fatbeta0-fatbeta).gt.ftol)then
        do i=1,ndim
          betacp(i)=beta(i)-beta0(i)
          beta0(i)=beta(i)
        enddo
        fatbeta0=fatbeta
        call linmin(beta,betamin,betamax,betacp,ndim,
     &f1dim_generic,fatbeta)
        call funkmin_generic(ndim,beta,fatbeta)
        if(isitnaninf(fatbeta).eq.0.and.fatbeta.lt.fatbeta0)goto 100
        fatbeta=fatbeta0
        do i=1,ndim
          beta(i)=beta0(i)
        enddo
      endif
110   call funkmin_generic(ndim,beta,fatbeta)
      do i=1,ndim0
        beta_in_out(i)=beta(i)
      enddo
      do i=1,npoints
        do j=1,nyvars
          shorty0(i,j)=shorty(i,j)
        enddo
        do j=1,nxvars
          shortx0(i,j)=shortx(i,j)
        enddo
      enddo
      return
      end subroutine GenericRegres
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$