PROGRAM MAIN

      implicit none
C     allocare xs,jcs,irs
C      nlam = 10
C      nnz = 10

      INTEGER*4 :: no,ni,nlam,maxit,ne,nx,nnz
      INTEGER*4 :: lmu,nlp,jerr
      integer*4 :: i, j
      REAL*4 :: flmin,thr,jd,parm,isd,intr,ka,ulam

      integer*4, dimension (:), allocatable :: irs, jcs, pcs;
      real*4, dimension (:), allocatable :: xs, y, w, vp;
      real*4, dimension (:,:), allocatable :: cl;


C      real xs(nnz)
C      real y(no)
C      real w(no)
C     real vp(ni)
C      real cl(2,ni)
C      real irs(nnz)
C      real jcs(nnz)

      integer*4, dimension (:), allocatable :: ia, nin;
      real*4, dimension (:), allocatable :: a0, rsq, alm;
      real*4, dimension (:,:), allocatable :: ca;
                                      
      
      EXTERNAL spelnet

C      PRINT*, 'Enter number of repeats'
C      READ*, nnz, no, ni, nlam


C     leggere file esterni

C      CALL read_idl_files("nnz.dat",nnz,1,1)
C      CALL read_idl_files("no.dat",no,1,1)
C      CALL read_idl_files("ni.dat",ni,1,1)
C      CALL read_idl_files("nlam.dat",nlam,1,1)

      OPEN(1, FILE="nnz.dat")
      READ(1, "(I16)") nnz
      CLOSE(1)

      OPEN(1, FILE="no.dat")
      READ(1, "(I16)") no
      CLOSE(1)

      OPEN(1, FILE="ni.dat")
      READ(1, "(I16)") ni
      CLOSE(1)

      OPEN(1, FILE="nlam.dat")
      READ(1, "(I16)") nlam
      CLOSE(1)

C     TO_DISPLAY
      PRINT *, 'nnz =' , nnz
      PRINT *, 'no =' , no
      PRINT *, 'ni =' , ni
      PRINT *, 'nlam =' , nlam

C     INPUT ARG
      ! allocate memory
      ! INTEGER INDEXES
      allocate(irs(nnz),stat=jerr)
      allocate(jcs(nnz),stat=jerr)
      allocate(pcs(ni+1),stat=jerr)

      ! REAL VECTORS / MATRICES
      allocate(xs(nnz),stat=jerr)
      allocate(y(no),stat=jerr)
      allocate(w(no),stat=jerr)
      allocate(vp(ni),stat=jerr)
      allocate(cl(2,ni),stat=jerr)
      
C     OUTPUT ARG
      allocate(ia(ni),stat=jerr)
      allocate(nin(nlam),stat=jerr)
      allocate(a0(nlam),stat=jerr)
      allocate(rsq(nlam),stat=jerr)
      allocate(alm(nlam),stat=jerr)
      allocate(ca(ni,nlam),stat=jerr)


      OPEN(1, FILE="irs.dat")
      DO i = 1, nnz
          READ(1, "(I16)") irs(i)
      END DO
      CLOSE(1)

      OPEN(1, FILE="jcs.dat")
      DO i = 1, nnz
          READ(1, "(I16)") jcs(i)
      END DO
      CLOSE(1)

      OPEN(1, FILE="pcs.dat")
      DO i = 1, ni+1
          READ(1, "(I16)") pcs(i)
C          PRINT *, 'pcs = ', pcs(i)
      END DO
      CLOSE(1)



      OPEN(1, FILE="xs.dat")
      DO i = 1, nnz
          READ(1, "(E16.8)") xs(i)
C          READ(1, "(F16.9)") xs(i)
      END DO
      CLOSE(1)

      OPEN(1, FILE="y.dat")
      DO i = 1, no
          READ(1, "(E16.8)") y(i)
      END DO
      CLOSE(1)

      OPEN(1, FILE="w.dat")
      DO i = 1, no
          READ(1, "(E16.8)") w(i)
      END DO
      CLOSE(1)

      OPEN(1, FILE="vp.dat")
      DO i = 1, ni
          READ(1, "(E16.8)") vp(i)
      END DO
      CLOSE(1)

      OPEN(1, FILE="cl.dat")
      DO i = 1, 2
          DO j = 1, ni
              READ(1, "(E16.8)") cl(i,j)
C              PRINT *, 'cl = ', i, j, cl(i,j)
          END DO
      END DO
      CLOSE(1)


      
C      CALL read_idl_files("xs.dat",xs,nnz,1)
C      CALL read_idl_files("y.dat",y,no,1)
C      CALL read_idl_files("w.dat",w,no,1)
C      CALL read_idl_files("vp.dat",vp,ni,1)
C      CALL read_idl_files("irs.dat",irs,nnz,1)
C      CALL read_idl_files("jcs.dat",jcs,nnz,1)
C      CALL read_idl_files("cl.dat",cl,2,ni)
      



C     ulam = 0 default 
      ulam = 0.0      
C     flmin
      flmin = 0.0001
C     thr
      thr = 1E-7
C     isd standardization = True
C     options.standardize Logical flag for x variable standardization, prior to
C                     fitting the model sequence. The coefficients are
C                     always returned on the original scale. Default is
C                     standardize = true. If variables are in the same
C                     units already, you might not wish to standardize. See
C                     details below for y standardization with
C                     family='gaussian'.
      isd = 1.0
C     options.intr    Should intercept(s) be fitted (default=true) or set
C                     to zero (false).
      intr = 1.0
         
C     with parm = 1.0 it does lasso, no elastic net
      parm = 1.0
C     with jd = 0 it considers all variables
      jd = 0.0
C     maxit (max number of iteration)
      maxit = 1E+5
      
C     gtype = 'naive';
      ka = 2.0
      if (ni < 500) then
C         gtype = 'covariance';
          ka = 1.0
      end if
      PRINT *, 'ka = ', ka
      
C     limit the numbers of predictors;
      ne = ni + 1;
      PRINT *, 'ne = ', ne

C     options.pmax        Limit the maximum number of variables ever to be
C                         nonzero. Default is min(dfmax * 2 + 20, nvars).
      nx = ni
      PRINT *, 'nx = ', nx
      PRINT *, 'isd = ', isd
      PRINT *, 'intr = ', intr
      PRINT *, 'parm = ', parm
      PRINT *, 'jd = ', jd
      PRINT *, 'maxit = ', maxit
      PRINT *, 'thr = ', thr
      call spelnet(ka,parm,no,ni,xs,pcs,irs,y,w,jd,vp,cl,ne,nx,
     $           nlam,flmin,ulam,thr,isd,intr,maxit,lmu,a0,ca,ia,nin,
     $           rsq,alm,nlp,jerr)

      

C     function new_lam = fix_lam(lam)
      alm(1)=exp( 2*log(alm(2)) - log(alm(3)) )
      
      OPEN(1, FILE="lmu.dat", STATUS='new')
          WRITE(1, "(I16)") lmu
      CLOSE(1)
      
      
      OPEN(1, FILE="ca.dat", STATUS='new')
      DO i = 1, nlam !lmu
          DO j = 1, ni
              WRITE(1, "(E16.8)") ca(j,i)
          END DO
      END DO
      CLOSE(1)
      
      
      OPEN(1, FILE="a0.dat", STATUS='new')
      DO i = 1, nlam !nlam
              WRITE(1, "(E16.8)") a0(i)
      END DO
      CLOSE(1)
      
      
      OPEN(1, FILE="ia.dat", STATUS='new')
      DO i = 1, ni
              WRITE(1, "(I16)") ia(i)
      END DO
      CLOSE(1)
      
      
      OPEN(1, FILE="nin.dat", STATUS='new')
      DO i = 1, nlam !lmu
              WRITE(1, "(I16)") nin(i)
      END DO
      CLOSE(1)
      
      OPEN(1, FILE="rsq.dat", STATUS='new')
      DO i = 1, nlam !lmu
              WRITE(1, "(E16.8)") rsq(i)
      END DO
      CLOSE(1)
      
      
      OPEN(1, FILE="alm.dat", STATUS='new')
      DO i = 1, nlam !lmu
              WRITE(1, "(E16.8)") alm(i)
C              PRINT * , 'alm = ', i , alm(i)
      END DO
      CLOSE(1)
      
      
      OPEN(1, FILE="nlp.dat", STATUS='new')
          WRITE(1, "(I16)") nlp
      CLOSE(1)
      
      
      OPEN(1, FILE="jerr.dat", STATUS='new')
          WRITE(1, "(I16)") jerr
C         PRINT *, 'jerr = ', jerr
      CLOSE(1)
      


C     scrivere
C     [lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr] 

C      CALL write_idl_files("lmu.dat",lmu,ni,1)
C      CALL write_idl_files("a0.dat",a0,nlam,1)
C      CALL write_idl_files("ca.dat",ca,ni,nlam)
C      CALL write_idl_files("ia.dat",ia,ni,1)
C      CALL write_idl_files("nin.dat",nin,nlam,1)
C     CALL write_idl_files("rsq.dat",rsq,nlam,1)
C      CALL write_idl_files("alm.dat",alm,nlam,1)
C      CALL write_idl_files("nlp.dat",nlp,1,1)
C     CALL write_idl_files("jerr.dat",jerr,1,1)

      END