anchor.diff:              it = it+1
           end do
   
  +    end subroutine nr_logistic_regression
  +
  +    subroutine logistic_regression(X, Y, B)
  +        implicit none
  +        real, intent(in) :: X(:,:)
  +        real, intent(in) :: Y(:)
  +        real(8), intent(out) :: B(:)
  +
  +        real(8), allocatable :: Yp(:), P(:), BN(:), V(:,:)
  +        real, allocatable :: YN(:), XV(:,:)
  +        integer :: nvars, ntimes, i, t, f, it
  +        real :: d
  +
  +        nvars = size(X,2) -1
  +        ntimes = size(Y)
  +
  +        allocate(BN(nvars+1))
  +        allocate(Yp(ntimes))
  +        allocate(YN(ntimes))
  +        allocate(P(ntimes))
  +        allocate(V(ntimes,ntimes))
  +        allocate(XV(ntimes, nvars+1))
  +
  +        do t = 1, ntimes, 1
  +            if(Y(t) .ne. 0.0) then
  +                Yp(t) = 1.0d0
  +            else
  +                Yp(t) = 0.0d0
  +            end if
  +        end do
  +
  +        B = 0.0d0
  +        f = 0
  +        i = 0
  +        it = 0
  +        do while (f /= 1)
  +            P = 1.0d0 / (1.0d0 + exp(-matmul(X, B))) 
  +            if (ANY(P > 0.97)) then
  +                !print *, "WARNING: logistic regression diverging"
  +                f = 1
  +            else
  +
  +                YN = Yp - P
  +                V = 0.0d0
  +                do t = 1, ntimes, 1
  +                    V(t,t) = P(t)*(1.0d0-P(t))
  +                enddo
  +                XV = matmul(V,X)
  +                call lapack_least_squares(XV, YN, BN)
  +
  +                f = 1
  +                do i = 1, nvars+1, 1
  +                    if(BN(i) .GT. 1.0E-04 .OR. BN(i) .LT. -1.0E-04) then
  +                        f = 0
  +                    end if
  +                end do
  +                if(it > 8) then
  +                    !print *, "WARNING: logistic regression failed to converge"
  +                    f = 1
  +                endif
  +
  +                B = B + BN
  +
  +            endif
  +            it = it+1
  +        end do
  +
       end subroutine logistic_regression
   
   
  
positive.code:         LDY = m
          
          nrhs = 1
          
          ! First find the optimal work size (LWORK)
          LWORK = -1
          CALL SGELS( 'N', M, N, NRHS, X, LDX, Y, LDY, WORK, LWORK, INFO )
          LWORK = MIN( 1000, INT( WORK( 1 ) ) )
          
          ! Now solve the equations X*B = Y
          CALL SGELS( 'N', M, N, NRHS, X, LDX, Y, LDY, WORK, LWORK, INFO )
          
          ! print*, "=============================="
          ! print*, X
          
          print*, "=============================="
          print*, Y(1:n)
          B(1:n) = Y(1:n)
  
      end subroutine lapack_least_squares
  
  
      !
  
negative.code:             y = y + coefficients(i) * x(i)
          end do
          
      end function nr_compute_regression
      
      function compute_regression(x, training_x, training_y, coefficients) result(y)
          implicit none
          real,    intent(in),    dimension(:)   :: x
          real,    intent(in),    dimension(:,:) :: training_x
          real,    intent(in),    dimension(:)   :: training_y
          real(8), intent(inout), dimension(:)   :: coefficients
          real,    allocatable,   dimension(:,:) :: training_x_lp
          real,    allocatable,   dimension(:)   :: training_y_lp
          
          real :: y
          
          integer :: nvars, ntimes, i
          
          nvars  = size(x)
          ntimes = size(training_x,1)
          
          allocate(training_x_lp(ntimes, nvars))
          training_x_lp = training_x
          
          allocate(training_y_lp(ntimes))
          training_y_lp = training_y
          
          if (maxval(training_y) == 0) then
              y=0
              return
          endif
          
          if (any(maxval(abs(training_x),dim=1) == 0.0)) then
              y=0
              return
          endif
          
          call lapack_least_squares(training_x_lp, training_y_lp, coefficients)
          
          y = 0
          if (maxval(abs(coefficients)) == 0) then
              y = sum(training_y) / size(training_y)
          endif
          do i=1,nvars
              y = y + coefficients(i) * x(i)
          end do
          
      end function compute_regression
  
      
      subroutine lapack_least_squares(X, Y, B, working_space)
          implicit none
          real, intent(inout), dimension(:,:) :: X
          real, intent(inout), dimension(:)   :: Y
          real(8), intent(inout), dimension(:)   :: B
          real, intent(inout), dimension(:), optional :: working_space ! temporary working space
          ! it could be more efficient to keep one space allocated and pass it in every time, 
          ! in practice, having a modest "WORK" allocated on the stack seems to be at least as fast. 
          real, dimension(1000) :: WORK
          integer :: n, m, nrhs, LDX, LDY, LWORK, INFO
          
  