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:         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:     
  
      
      subroutine nr_logistic_regression(X, Y, TX, B)
          implicit none
          real, intent(in) :: X(:,:)
          real, intent(in) :: Y(:)
  