anchor.diff:          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)
  +
  +        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
   
   
       !
  
positive.code:     use nr
  contains
      
      function compute_regression(x, training_x, training_y, coefficients) result(y)
          implicit none
          real,    intent(in),    dimension(:)   :: x
          real,    intent(in),    dimension(:,:) :: training_x
  
negative.code:     
  
      
      subroutine nr_logistic_regression(X, Y, TX, B)
          implicit none
          real, intent(in) :: X(:,:)
          real, intent(in) :: Y(:)
  