subroutine samplingscheme(sampfunc,nparams,msamp,
     &  bestguess,guessconfid0,bmax,bmin,params,ysamp)
!
! This subroutine conducts samples from a given function using both random and
! systematic strategies. In the random strategy, sampled points randomly scatter 
! around the bestguess point. The closeness to the bestguess depends on the confidence
! on the bestguess. For any given parameter, the width of the last interval is 
! guessconfid0 X the width of the first interval. In the systematic strategy, two values 
! for each parameter are determined with the new bestguess after
! the random sampling in the middle. The two determined values of all parameters
! are systematically combined to form new additional sampling points. 
!
! After the two sampling strategies, the latest best guess is the first point in params.

      
      implicit none
      integer nparams,msamp
      double precision bestguess(nparams),guessconfid0,
     &    bmax(nparams),bmin(nparams),params(nparams,msamp),
     &    ysamp(msamp),guessconfid
      integer i,nright,nleft,ibest,msamptemp
      double precision accum,x1,delta,j,ybest,temp
      
      external sampfunc
      
      guessconfid=dmax1(1.0d0,guessconfid0)
      
      do i=1,nparams
        if(bestguess(i).lt.bmin(i).or.bestguess(i).gt.
     &    bmax(i))then
          write(*,*)'best guess out of bounds, sampling stops'
          stop
        endif
      enddo
      msamptemp=msamp-1-2*nparams
      
      if(msamptemp.lt.0)then
        write(*,*)'sampling number must be larger than',
     &      (1+2*nparams)
        stop
      endif
      if(msamptemp.gt.0)then
        if(msamptemp.ge.3)then
          if(mod(msamptemp,2).eq.0)then
            nright=msamptemp/2
            nleft=msamptemp/2
          else
            nright=msamptemp/2+1
            nleft=msamptemp/2
          endif
          do i=1,nparams
!first divide the right        
            x1=2.0d0*(bmax(i)-bestguess(i))/
     &        (dble(nright)*(guessconfid+1.0d0))
            delta=x1*(guessconfid-1.0d0)/dble(nright-1)
            accum=0.0d0
            do j=1,nright
              accum=accum+x1+dble(j-1)*delta
              params(i,j)=accum+bestguess(i)
            enddo
!
!then divide the left
            x1=2.0d0*(bestguess(i)-bmin(i))/
     &        (dble(nleft)*(guessconfid+1.0d0))
            delta=x1*(guessconfid-1.0d0)/dble(nleft-1)
            accum=0.0d0
            do j=1,nleft
              accum=accum+x1+dble(j-1)*delta
              params(i,j+nright)=bestguess(i)-accum
            enddo      
          enddo
        else
          if(msamptemp.eq.1)then
            do i=1,nparams
! arbitrarily take one value
              params(i,1)=bestguess(i)+(bmax(i)-bestguess(i))
     &                    *0.339354235d0
            enddo
          endif
          if(msamptemp.eq.2)then
            do i=1,nparams
! arbitrarily take two values              
              params(i,1)=bestguess(i)+(bmax(i)-bestguess(i))
     &                    *0.339354235d0
              params(i,2)=bestguess(i)-(bestguess(i)-bmin(i))
     &                    *0.339354235d0
            enddo
          endif
        endif
        call randpermut_dim_samp(msamptemp,nparams,
     &        params(1:nparams,1:msamptemp))
      endif
      msamptemp=msamptemp+1
      do i=1,nparams
        params(i,msamptemp)=bestguess(i)
      enddo
      do i=1,msamptemp
        call sampfunc(nparams,params(1:nparams,i:i),ysamp(i))
      enddo
      
      ibest=1
      ybest=ysamp(ibest)
      do i=2,msamptemp
        if(ysamp(i).lt.ybest)then
          ibest=i
          ybest=ysamp(i)
        endif
      enddo
      
      do i=1,nparams
        msamptemp=msamptemp+1
        params(i,msamptemp)=params(i,ibest)+
     &       (bmax(i)-params(i,ibest))/(1.0d0+guessconfid)
        msamptemp=msamptemp+1
        params(i,msamptemp)=params(i,ibest)-
     &       (params(i,ibest)-bmin(i))/(1.0d0+guessconfid)
        do j=1,i-1
          params(j,msamptemp-1)=params(j,ibest)
          params(j,msamptemp)=params(j,ibest)
        enddo
        do j=i+1,nparams
          params(j,msamptemp-1)=params(j,ibest)
          params(j,msamptemp)=params(j,ibest)
        enddo
      enddo
      do i=msamptemp-2*nparams+1,msamptemp
        call sampfunc(nparams,params(1:nparams,i:i),ysamp(i))
        if(ysamp(i).lt.ybest)then
          ibest=i
          ybest=ysamp(i)
        endif
      enddo
      do i=1,nparams
        temp=params(i,1)
        params(i,1)=params(i,ibest)
        params(i,ibest)=temp
      enddo
      temp=ysamp(1)
      ysamp(1)=ysamp(ibest)
      ysamp(ibest)=temp
      
      return
      end