!= silent.f90 12Nov2017 Silent sigma from sums of unequal size samples
! ... Keywords: silent bags unequal size samples loads Gaussian sigma Monte Carlo
! ... Silent from loads of bags to sigma ...............................
 IMPLICIT DOUBLE PRECISION (a-h, o-z)
 POINTER lotsize_vec(:)
 COMMON /params/ verbose
 LOGICAL verbose
 ALLOCATABLE x_vec(:), freqd(:), acfreq(:)
 INTERFACE
    SUBROUTINE readat (lots, lotsize_vec, zmu, sigma, criter, &
         ntr, iseed, ishow)
    IMPLICIT DOUBLE PRECISION (a-h, o-z)
    POINTER lotsize_vec(:)
    END
 END INTERFACE
 DATA one/ 1 /
 verbose = .TRUE.
 CALL readat (lots, lotsize_vec, zmu, sigma, criter, &
      ntrials, iseed, ishow)
 lot_tot = SUM(lotsize_vec); runs = ntrials * lot_tot
 IF (verbose) WRITE (*, 96730) lot_tot, LOG10(runs)
96730 FORMAT (' Total lots size,', T28, I10, T42, '| (sum of sizes)' / &
           ' (Trials).(total size) =', T30, '10 ^', F4.1, T42, &
           '| (run length)')
 IF (verbose) WRITE (*, 96740) SUM(lotsize_vec)*zmu/lots, SUM(lotsize_vec)*one/lots
96740 FORMAT (' Typical load (Av.mu),', T30, 1P, G12.4, T42, &
           '| Av = average(lotsize_vec),', 1P, G9.3)
 CALL randomset (iseed)
 CALL tempfile ('U', ifile, 'on')
! ... One criterion only is embedded in the simulation
 CALL simulate (ifile, criter, zmu, sigma, ntrials, lots, lotsize_vec)
 kl0 = 0; kls = 200
 ALLOCATE (x_vec(kl0:kls), freqd(kl0:kls), acfreq(kl0:kls))
 CALL classify (ifile, ntrials, kl0, kls, x_vec, freqd, acfreq, &
      xmin, xmax, xmode, average, stdev)
 INQUIRE (ifile, SIZE=ifilesize)
 IF (verbose) WRITE (*, 99090) one*ifilesize, 1.D-6*ifilesize
99090 FORMAT (' Temporary file size (bytes),', T30, EN12.2, T42, &
           '| MB, ', 1P, G8.3)
!! CALL tempfile ('U', ifile, 'off')
 IF (verbose) WRITE (*, 99110) xmin, xmax
99110 FORMAT (3X, 'min, max,', T18, 1P, 2G12.4, T42, '|')
 IF (verbose) WRITE (*, 99550) 'mode', xmode, 'average', average, 'stdev', stdev
99550 FORMAT (3X, A, ',', T30, 1P, G12.4, T42, '|')
 IF (verbose) WRITE (*, 99570) average-3*stdev, average+3*stdev
99570 FORMAT (' ave -/+ 3 std,', T18, 1P, 2G12.4, T42, '|')
 factor = SQRT(one*SUM(lotsize_vec))
 IF (verbose) WRITE (*, 99590) stdev*factor, stdev*factor/sigma
99590 FORMAT (T42, '|' / ' Estimated sigma,', T30, 1P, G12.4, T42, &
           '| est./sigma, ', 1P, G10.5, T67, '(~1 expectable)')
 STOP
 END !END!

 SUBROUTINE simulate (ifile, criter, zmu, sigma, ntrials, &
      lots, lotsize_vec)
! ... ..................................................................
 IMPLICIT DOUBLE PRECISION (a-h, o-z)
 DIMENSION lotsize_vec(lots)
 ALLOCATABLE sGauss_vec(:), w_vec(:), ave_vec(:)
 maxlot = MAXVAL(lotsize_vec)
 ALLOCATE (sGauss_vec(maxlot), w_vec(lots), ave_vec(lots))
 w_vec = lotsize_vec**criter / SUM(lotsize_vec**criter)
 REWIND (ifile)
 DO n=1, ntrials ! Simulate the sums (load weights)
    DO lot=1, lots; nss = lotsize_vec(lot)
       CALL rand_sGauss (nss, sGauss_vec(:nss)) ! standard Gaussian
       ave_vec(lot) = zmu + sigma * SUM(sGauss_vec(:nss)) / nss
    END DO
    equivweight = SUM(w_vec*ave_vec)
    WRITE (ifile) equivweight ! auxfile
 END DO
 DEALLOCATE (sGauss_vec, w_vec, ave_vec)
 RETURN
 END !END!

 FUNCTION standard_deviation (n, x_vec)
! ... ..................................................................
 IMPLICIT DOUBLE PRECISION (a-h, o-z)
 DIMENSION x_vec(n)
 ave = SUM(x_vec) / n
 var = SUM((x_vec-ave)**2) / (n - 1)
 standard_deviation = SQRT(var)
 RETURN
 END !END!

! not used
 SUBROUTINE combine_x (kriter, lots, lotsize_vec, ave_vec, stdev)
! ... Weighting criteria ...............................................
 IMPLICIT DOUBLE PRECISION (a-h, o-z)
 DIMENSION lotsize_vec(lots), ave_vec(lots)
 ALLOCATABLE w_vec(:)
 INTENT (IN)  :: lots, lotsize_vec, ave_vec
 INTENT (OUT) :: stdev
 DATA one/ 1 /
 ALLOCATE (w_vec(lots))
! ... Below: always no zero weight
 SELECT CASE (kriter)
 CASE (0) ! equal importance
 n = lots ! for simplicity
 ave_W = SUM(ave_vec) / lots
 variance = n * SUM((ave_vec-ave_w)**2) / (n - 1)
 iequiv = SUM(lotsize_vec)
 stdev = SQRT(variance / iequiv)
 CASE (1) ! alternative
 pow = 2 !! 3 * one / 2
 w_vec = lotsize_vec**pow / SUM(lotsize_vec**pow) ! normalize
 ave_W = SUM(w_vec*x_vec)
 variance = n * SUM(w_vec*(x_vec-ave_w)**2) / (n - 1)
 ave_lotsize = SUM(lotsize_vec) * one / lots
 stdev = SQRT(variance) * SQRT(ave_lotsize)
 END SELECT
 DEALLOCATE (w_vec)
 RETURN
 END !END!

 SUBROUTINE classify (ifile, ntr, kl0, kls, x_vec, freqd, acfreq, &
      xmin, xmax, xmode, average, stdev)
! ... ..................................................................
 IMPLICIT DOUBLE PRECISION (a-h, o-z)
 DIMENSION x_vec(kl0:kls), freqd(kl0:kls), acfreq(kl0:kls)
 ALLOCATABLE ifreq(:)
 xmin = HUGE(xmin); xmax = -xmin
 REWIND (ifile)
 DO i=1, ntr; READ (ifile) xi ! auxfile
    IF (xi < xmin) xmin = xi; IF (xi > xmax) xmax = xi
 END DO
 clwidth = (xmax - xmin) / (kls - kl0 + 1)
 x_vec = (/ (xmin + k * clwidth, k=kl0, kls) /)
 ALLOCATE (ifreq(kl0:kls+1))
 REWIND (ifile); ifreq = 0; sumx = 0
 DO i=1, ntr; READ (ifile) xi ! auxfile
    sumx = sumx + xi
    kl = kl0 + (xi - xmin) / clwidth; ifreq(kl) = ifreq(kl) + 1
 END DO
 ifreq(kls) = ifreq(kls) + ifreq(kls+1) ! absorb last
 freqd = ifreq / (ntr * clwidth)
 acfreq(kl0) = ifreq(kl0)
 DO i=kl0+1, kls; acfreq(i) = acfreq(i-1) + ifreq(i)
 END DO
! ... Mode (beware: MAXLOC always starts in 1)
 mode = LBOUND(ifreq, DIM=1) - 1 + MAXLOC(ifreq, DIM=1)
 DEALLOCATE (ifreq)
 acfreq = acfreq / ntr
 xmode = x_vec(mode); average = sumx / ntr
 REWIND (ifile); sumx2 = 0
 DO i=1, ntr; READ (ifile) xi ! auxfile
    sumx2 = sumx2 + (xi - average) * (xi - average)
 END DO
 stdev = SQRT(sumx2 / (ntr - 1))
! ... Welford method for comparison
! CALL stdev_Welford (ifile, ntr, aver_W, std_W)
 RETURN
 END SUBROUTINE !END!

 SUBROUTINE rand_sGauss (nn, sGauss_vec)
! ... ..................................................................
 IMPLICIT DOUBLE PRECISION (a-h, o-z)
 DIMENSION sGauss_vec(nn)
 ALLOCATABLE u_vec(:)
 DATA pi/ 3.141592653589793d0 /
 tini = TINY(tini)
! ... Box-Muller
 ALLOCATE (u_vec(-1:nn))
 CALL RANDOM_NUMBER (u_vec); WHERE (u_vec < tini) u_vec = tini
 DO i=2, nn, 2
    sqr = SQRT(-2*LOG(u_vec(i-1))); arg = 2 * pi *u_vec(i)
    sGauss_vec(i-1) = sqr * COS(arg)
    sGauss_vec(i)   = sqr * SIN(arg)
 END DO
 IF (nn/2*2 /= nn) THEN
    sqr = SQRT(-2*LOG(u_vec(-1))); arg = 2 * pi *u_vec(0)
    sGauss_vec(nn) = sqr * COS(arg) ! or SIN
 END IF
 DEALLOCATE (u_vec)
 RETURN
 END !END!

 SUBROUTINE randomset (iseed)
! ... Set random .......................................................
 ALLOCATABLE kseed(:)
 CALL RANDOM_SEED (SIZE=isize)
 ALLOCATE (kseed(isize))
95600 FORMAT (T42, '| (Random "size",', I3, ')')
 CALL RANDOM_SEED () ! Random set by processor or 'put' user seed
 IF (iseed /= 0) THEN; kseed = iseed; CALL RANDOM_SEED (PUT=kseed)
 END IF
 RETURN
 END !END!

 SUBROUTINE readat (lots, lotsize_vec, zmu, sigma, criter, &
      ntrials, iseed, ishow)
! ... ..................................................................
 IMPLICIT DOUBLE PRECISION (a-h, o-z)
 CHARACTER (LEN=3) :: hnoyes(0:1)=(/ ' No', 'Yes' /), &
      hbuffer*1024, i0*16
 LOGICAL verbose
 COMMON /params/ verbose
 POINTER temp_vec(:), lotsize_vec(:)
 INTERFACE
    SUBROUTINE vec_read (h, n, vector)
    IMPLICIT DOUBLE PRECISION (a-h, o-z)
    CHARACTER h*(*); POINTER vector(:)
    END
 END INTERFACE
 DATA one/ 1 /
! ... Get the sample ("lot") sizes (as reals: no available routine for int.s)
 CALL vec_read (hbuffer, lots, temp_vec)
! ... Sizes are integers
 ALLOCATE (lotsize_vec(lots)); lotsize_vec = NINT(temp_vec)
 DEALLOCATE (temp_vec)
 lot_tot = SUM(lotsize_vec)
 IF (verbose) WRITE (*, 95100) lots, lotsize_vec
95100 FORMAT (' Lots (N. of),', T32, I6, T42, '|' / &
           ' Lot (as sample) sizes:', T42, '|' / &
           (T42, '|', T8, 10I6))
 READ (*, *) zmu, sigma
 IF (verbose) WRITE (*, 95600) zmu, sigma
95600 FORMAT (' .mu, sigma,', T18, 1P, 2G12.4, T42, '|')
 READ (*, *, IOSTAT=nan) kriter, criter
 IF (nan /= 0) THEN; criter = kriter
    IF (verbose) WRITE (*, 95050) kriter
95050 FORMAT (' --- Criterion,', T30, I8, T42, '|')
 ELSE; IF (verbose) WRITE (*, 95060) criter
95060 FORMAT (' --- Criterion,', T30, 1P, G12.4, T42, '|')
 END IF
 READ (*, *) trials_lg, iseed
 zlg = LOG10(one*lot_tot)
 ntrials = NINT(10**trials_lg)
 IF (verbose) WRITE (*, 96500) trials_lg, iseed, TRIM(i0(ntrials)), &
      hnoyes(MIN(iseed,1))
96500 FORMAT (' Trials, seed,', T18, '10^', F5.3, I12, T42, &
           '| (', A, ')', T64, 'Repeatability: ', A3)
 READ (*, *) ishow; IF (verbose) WRITE (*, 98900) hnoyes(MIN(ishow,1))
98900 FORMAT (' Show coord.s ?', T35, A3, T42, '|' / &
   1X, 20('--'), '+-', 19('--'))
 RETURN
 END !END!

 FUNCTION i0 (n)
! ... ..................................................................
 CHARACTER i0*(*), hbuffer*32
 WRITE (hbuffer, *) n; i0 = ADJUSTL(hbuffer)
 RETURN
 END !END!

 SUBROUTINE from_excel (string)
! ... ..................................................................
 CHARACTER string*(*), tab, h_vec
 INTENT (INOUT) :: string
 LOGICAL :: ltest=.FALSE.
 ALLOCATABLE h_vec(:)
 string = ADJUSTL(string); leng = LEN_TRIM(string)
IF (ltest) print*,'_from_excel IN  string:_', trim(string) // '_'
 ipoint = 0; IF (INDEX(string(:leng), '.') > 0) ipoint = 1
 isemic = 0; IF (INDEX(string(:leng), ';') > 0) isemic = 1
 iamer = 1; IF (ipoint == 0 .OR. isemic > 0) iamer = 0
IF (ltest) print*,'_from_excel iAmer,', iamer
! ... Replacements
 ALLOCATE (h_vec(leng))
 h_vec = (/ (string(k:k), k=1, leng) /); tab = CHAR(9)
 SELECT CASE (iamer)
 CASE (0) ! ... De-tab and convert to Amer
 WHERE (h_vec == ',') h_vec = '.'
 WHERE (h_vec == tab) h_vec = ','
 WHERE (h_vec == ';') h_vec = ','
 CASE (1) ! ... De-tab
 WHERE (h_vec == tab) h_vec = ','
 END SELECT
 DO k=1, leng; string(k:k) = h_vec(k)
 END DO
 DEALLOCATE (h_vec)
IF (ltest) print*,'_from_excel OUT string:_', trim(string) // '_'
 RETURN
 END SUBROUTINE !END!

 SUBROUTINE vec_read (hwork, ndata, vector)
! ... Reads a vector and "converts" it to Amer .........................
 IMPLICIT DOUBLE PRECISION (a-h, o-z)
 CHARACTER hwork*(*)
 LOGICAL :: ltest=.FALSE.
 POINTER vector(:)
 CALL tempfile ('F', itempfile, 'on')
 REWIND (itempfile)
 ndata = 1
 DO; READ (*, "(A)") hwork; IF (hwork(:3) == 'EOD') EXIT
    hwork = ADJUSTL(hwork); IF (LEN_TRIM(hwork) == 0) CYCLE
IF (ltest) print*,'_vec_read 11 hwork:_', trim(hwork) // '_'
    CALL from_excel (hwork)
    CALL suppress_blanks (hwork)
IF (ltest) print*,'_vec_read 22 hwork:_', trim(hwork) // '_'
    WRITE (itempfile, "(A)") TRIM(hwork)
    leng = LEN_TRIM(hwork)
    DO i=2, leng; IF (hwork(i:i) == ',') ndata = ndata + 1
    END DO
 END DO
 ALLOCATE (vector(ndata))
 REWIND (itempfile)
 READ (itempfile, *) vector
 CALL tempfile ('F', itempfile, 'off')
 RETURN
 END SUBROUTINE !END!

 SUBROUTINE suppress_blanks (string)
! ... ..................................................................
 CHARACTER string*(*)
 LOGICAL :: ltest=.FALSE.
 DO; ib = INDEX(TRIM(string), ' '); IF (ib == 0) EXIT
IF (ltest) print*,CHAR(13), '_suppress 22 Leng, string:', &
         len_trim(string), ' _'//trim(string)//'_'
    IF (string(ib-1:ib-1) /= ',') string(ib:ib) = ','
IF (ltest) print*,'_suppress 33 ib,', ib, CHAR(13), &
         '_suppress string(ib+1:ib+1):_', string(ib+1:ib+1)//'_'
    IF (string(ib:ib) == ' ') THEN; leng = LEN_TRIM(string)
IF (ltest) print*,'_suppress 44 leng, string:', leng, ' _'//trim(string)//'_'
       string(ib:leng-1) = string(ib+1:leng)
       string(leng:leng) = ' '
IF (ltest) print*,'_suppress 55 Lg, string:', len_trim(string), &
            ' _'//trim(string)//'_'
    END IF
IF (ltest) print*,'_suppress 66 Leng, string:', len_trim(string), &
         ' _'//trim(string)//'_'
 END DO
 RETURN
 END SUBROUTINE !END!

 SUBROUTINE tempfile (hform, ifile, honoff)
! ... Manages a temporary file .........................................
 CHARACTER hform, honoff*(*), hpid*64, htime*10, form*11
 CHARACTER (LEN=255) :: hbuffer, scratch='/tmp/1038_'
 INTEGER GETPID
 IF (honoff == 'on') THEN
    WRITE (hpid, *) GETPID() ! till 2^15 = 32768 (<=five figures)
    hpid = ADJUSTL(hpid)
    CALL DATE_AND_TIME (TIME=htime)
! ... hhMMss.mmm --> smmm
    hbuffer = TRIM(hpid) // htime(6:6)//htime(8:10)
! ... 'ifile' is 'getpid' & 'smmm'
    READ (hbuffer, *) ifile
! ... 'File name is '...getpid' & 'ifile' & 'smmm'
    scratch = TRIM(scratch) // TRIM(hpid) // htime(8:10)
    SELECT CASE (hform)
    CASE ('U'); form = 'UNFORMATTED'
    CASE ('F'); form = 'FORMATTED'
    CASE DEFAULT; CALL error ('Wrong FORM for temporary file.')
    END SELECT
    OPEN (ifile, FILE=TRIM(scratch), FORM=TRIM(form), STATUS='UNKNOWN')
 ELSE IF (honoff == 'off') THEN; CLOSE (ifile, STATUS='DELETE')
 ELSE; CALL error ('Wrong file action.')
 END IF
 RETURN
 END SUBROUTINE !END!

 SUBROUTINE error (htext)
! ... ..................................................................
 CHARACTER htext*(*)
 WRITE (*, "(' %Err ', A, ' Stop.')") htext; STOP
 ENTRY warning (htext)
 WRITE (*, "(' %Wrn ', A)") htext
 RETURN
 END SUBROUTINE !END!