anchor.diff:              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)
  +
  +    
  +    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
           
  
positive.code:             it = it+1
          end do
  
      end subroutine logistic_regression
  
  
  
negative.code:         LDY = m
          
          nrhs = 1
  
          if (present(working_space)) then
              ! First find the optimal work size (LWORK)
              LWORK = -1
              CALL SGELS( 'N', M, N, NRHS, X, LDX, Y, LDY, working_space, LWORK, INFO )
              LWORK = MIN( size(working_space), INT( working_space( 1 ) ) )
              
              ! Now solve the equations X*B = Y
              CALL SGELS( 'N', M, N, NRHS, X, LDX, Y, LDY, working_space, LWORK, INFO )
          else
              ! 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 )
          endif
          if (INFO/=0) then
              write(*,*) "ERROR in SGELS with argument:", 0-INFO
              Y = 0
          endif
  
          B(1:n) = Y(1:n)
  
      end subroutine lapack_least_squares
  
      subroutine weighted_least_squares(X, Y, B, W)
          implicit none
          real, intent(inout), dimension(:,:) :: X
          real, intent(inout), dimension(:)   :: Y, W
          real(8), intent(inout), dimension(:)   :: B
          real, dimension(:,:), allocatable :: X_weighted, W_full
          real, dimension(:),   allocatable :: Y_weighted, Bs, output_y
          
          real, dimension(10000) :: WORK
          integer :: i, n, m, nrhs, LDX, LDY, LDW, LWORK, INFO, p
          
          ! ntimes
          m = size(X,1)
          LDX = m
          LDY = m
          LDW = m
          p   = m
          
          ! nvars
          n = size(X,2)
          
          nrhs = 1
  
          allocate(output_y(p))
          allocate(Bs(n))
          allocate(W_full(m,m))
          W_full = 0
          do i=1,m
              W_full(i,i) = 1/W(i)
          enddo
          
          
          ! print*, "entering sggglm"
          LWORK = 10000
          call sggglm(m, n, p, X, LDX, W_full, LDW, Y, Bs, output_y, WORK, LWORK, INFO)
          
          B(1:n) = Bs(1:n)
  
      end subroutine weighted_least_squares
  
  
      !
  