C
Code by Michael Kupferschmid
C
C     This program generates linear programming problems having
C     random constant columns, solves them, and collects the values
C     of the resulting regression coefficients.
C
C     I/O unit  use
C     --------  -----------------------------
C         1     input data
C         2     L2 regression BETA(2) values
C         3     LAV regression BETA(2) values
C         5     R read from keyboard
C         6     results printed on screen
C
C     variable  meaning
C     --------  -------
C     ARG       command-line specification of R
C     BETA      coefficients in linear fit y=BETA(1)+BETA(2)*x
C     BETAS1    L1 slopes
C     BETAS2    L2 slopes
C     DFLOAT    Fortran function returns REAL*8 for INTEGER*4
C     DR250     routine returns next random vector uniform on [0,1]
C     DY        disturbance added to Y(I), uniform on [-R,R]
C     EPS       vector of pseudorandom numbers uniform on [0,1]
C     GETARG    system routine gets one command-line argument
C     GETFIL    routine attaches a file for a given purpose
C     HI        upper boundary of a histogram bin
C     I         index on data pairs read in
C     IROUND    function returns REAL*8 rounded to nearest INTEGER*4
C     J         index on histogram bins
C     K         index on experiments
C     L         a histogram bin count
C     L1REGW    routine does weighted least-abs-values regression
C     L2REGW    routine does weighted linear regression
C     LO        lower boundary of a histogram bin
C     M         number of dots to print for a histogram bar
C     N         number of data pairs read in
C     NDOTS     total number of dots in the histogram
C     NEXP      number of experiments to run
C     NMAX      maximum number of data points allowed
C     R         amplitude of the random disturbances added to Y
C     RC        return codes
C     W         weights (all 1.D0)
C     X         x values
C     XIN       an x value
C     Y         y values
C     YIN       a y value
C
      PARAMETER(NMAX=1000)
      REAL*8 BETA(2),W(NMAX),X(NMAX),Y(NMAX),EPS(NMAX)
      REAL*8 DY,HI,LO,R,XIN,YIN
      PARAMETER(NDOTS=30,NEXP=30*NDOTS)
      REAL*8 BETAS1(NEXP),BETAS2(NEXP)
      CHARACTER*13 ARG
      INTEGER*4 RC
C
C ------------------------------------------------------------------
C
C     read and count (x,y) pairs of input data
      OPEN(UNIT=1,FILE='westwood.data')
      DO 1 I=1,NMAX
           READ(1,*,END=2,ERR=3) XIN,YIN
           N=I
           X(N)=XIN
           Y(N)=YIN
           W(N)=1.D0
    1 CONTINUE
      WRITE(0,901) NMAX
  901 FORMAT('more than',I5,' data pairs')
      STOP
    3 WRITE(0,902) N+1
  902 FORMAT('read error at line',I5)
      STOP
    2 IF(N.LT.2) THEN
         WRITE(0,903) N
  903    FORMAT('N=',I1,' does not permit a regression')
         STOP
      ENDIF
C
C     find out how much noise to add to the Y vector
      CALL GETARG(1,ARG)
      READ(ARG,*) R
C
C     run experiments
      DO 4 K=1,NEXP
C          add uniformly-distributed noise to the Y vector
           CALL DR250(N,EPS)
           DO 5 I=1,N
                DY=2.D0*R*(EPS(I)-0.5D0)
                Y(I)=Y(I)+DY
    5      CONTINUE
C
C          compute the L2 regression and remember the results
           CALL L2REGW(X,Y,W,N, BETA,RC)
           IF(RC.NE.0) THEN
              WRITE(0,904) RC
  904         FORMAT('least-squares regression failed, RC=',I2)
              STOP
           ELSE
              BETAS2(K)=BETA(2)
           ENDIF
C
C          compute the L1 regression and remember the results
           CALL L1REGW(X,Y,W,N, BETA,RC)
           IF(RC.NE.0) THEN
              WRITE(0,905) RC
  905         FORMAT('LAV regression failed, RC=',I2)
              STOP
           ELSE
              BETAS1(K)=BETA(2)
           ENDIF
    4 CONTINUE
C
C     histogram the L2 results
      CALL GETFIL('L2 histogram',12,' ',0,2,3)
      WRITE(2,906) R
  906 FORMAT('L2, r=',F4.1)
      DO 6 J=1,8
           LO=-1.5D0+DFLOAT(J-1)
           HI=LO+1.D0
           L=0
           DO 7 K=1,NEXP
                IF(BETAS2(K).GE.LO .AND. BETAS2(K).LT.HI) L=L+1
    7      CONTINUE
           M=IROUND(DFLOAT(L)/DFLOAT(NEXP/NDOTS))
           WRITE(2,907) LO,('.',K=1,M)
  907      FORMAT(F9.2,1X,30A1)
    6 CONTINUE
C
C     histogram the L1 results
      CALL GETFIL('L1 histogram',12,' ',0,3,3)
      WRITE(3,908) R
  908 FORMAT('L1, r=',F4.1)
      DO 8 J=1,8
           LO=-1.5D0+DFLOAT(J-1)
           HI=LO+1.D0
           L=0
           DO 9 K=1,NEXP
                IF(BETAS1(K).GE.LO .AND. BETAS1(K).LT.HI) L=L+1
    9      CONTINUE
           M=IROUND(DFLOAT(L)/DFLOAT(NEXP/NDOTS))
           WRITE(3,907) LO,('.',K=1,M)
    8 CONTINUE
      STOP
      END