! compile with -Dstandalone for standalone version

! meanandvar -- Copyright 2015 Paolo Luchini
! http://CPLcode.net/Applications/Numerical/MeanandVar/

!( Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
  in the Software without restriction, including without limitation the rights
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the Software is
  furnished to do so, subject to the following conditions:
  
  The above copyright notice and this permission notice shall be included in all
  copies or substantial portions of the Software.
  
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE. !)

TEST=NO
BOOLEAN autots=YES

MEANANDVAR=STRUCTURE[REAL mean, bm(0..1), bt(0..1)
  INTEGER nt, bs, ts]

SUBROUTINE mvinit(MEANANDVAR mvp^) WITH mvp
  mean=0; bm=0; bt=0
  nt=0; bs=0; ts=1
END mvinit

SUBROUTINE meanandvar(MEANANDVAR mvp^;REAL xx) WITH mvp
  bm(0)=xx-mean+~
  INC bs
  IF bs>=ts THEN
    nt=~+bs
    mean=~+bm(0)/nt
    bt(1)=~+bm(0)∗bm(1)
    bm(1)=bm(0)∗(1-bs/nt)
    bt(0)=~+bm(0)∗bm(1)
    IF bt(1)>0.5∗bt(0) AND autots AND nt>bs THEN ts=CEILING(1.05∗ts)
    bs=0
    bm(0)=0
  END IF
END meanandvar

REAL FUNCTION var(MEANANDVAR mvp^) WITH mvp
    IF nt-3∗ts<=0 THEN RETURN 0
    RESULT=[bt(0)+2∗bt(1)]/nt/(nt-3∗ts)
    IF bt(1)>0 THEN RESULT=~∗(1+0.234∗[2∗bt(1)/bt(0)]^3.17)
END var

REAL FUNCTION rms(MEANANDVAR mvp^)=IF mvp.var>0 THEN SQRT(mvp.var) ELSE 0
REAL FUNCTION sum(MEANANDVAR mvp^)=mvp.mean∗mvp.nt
REAL FUNCTION sumvar(MEANANDVAR mvp^)=mvp.var∗mvp.nt^2
REAL FUNCTION sumrms(MEANANDVAR mvp^)=mvp.rms∗mvp.nt

#ifdef standalone

INTEGER startcolumn=1, numcolumn=1
IF COMMANDLINE.HI>=1 THEN
  IF COMMANDLINE(1)="-h" THEN
    WRITE <<helpends
    
Usage: meanandvar [-h] [col[-tocol]] [skip] [bsize]

Estimates the mean and the standard deviation of the estimate of the mean from
finite time series extracted from a stationary stochastic process. The process'
correlation function is assumed to be stationary and integrable but otherwise
unspecified. The estimate of the standard deviation is obtained by an unbiased
method of batch means with internally calculated adaptive batch size.

Input is accepted from stdin in a multicolumn format, possibly prefixed by
comment lines starting with "#". Each column is assumed to be a separate time
series. The first commandline parameter, if present, specifies the column, or
range of columns, to be operated on (default: 1). The second parameter, if
present, is a number of lines to be skipped at the beginning of the file
(default: 0). The third parameter, if present, is a fixed batch size to be used
(default: adaptive).

Output is its estimate of the mean and standard deviation for each selected
column, one per line. A final line contains the total number of samples received
and the automatic batch size (roughly proportional to correlation time).
helpends
    STOP
  END IF
  startcolumn=atoi(COMMANDLINE(1))
  mpos=strchr(COMMANDLINE(1),"-")
  IF mpos#NULL THEN
    numcolumn=atoi(mpos(1))-startcolumn+1
  END IF
END IF
skip=IF COMMANDLINE.HI>=2 THEN atoi(COMMANDLINE(2)) ELSE 0
MEANANDVAR mv(1..numcolumn)=0
DO mvinit(mv(i)) FOR ALL i
IF COMMANDLINE.HI>=3 THEN mv(1).ts=atoi(COMMANDLINE(3)); autots=NO
REAL xx
LOOP WHILE ungetc(getc(stdin),stdin)="#": READ
LOOP FOR i=1 TO skip: READ
LOOP WHILE READ xx
  LOOP FOR i=2 TO startcolumn: READ xx
  LOOP FOR col=1 TO numcolumn
    meanandvar[mv(col),xx]
    IF col<numcolumn THEN READ xx
  REPEAT
  READ
  IF mv(1).bs=0 THEN
    WITH mv(1)
    IF TEST AND nt-3∗ts>0 THEN
      sigma=var
      sigma10=[bt(0)+bt(1)]/nt/(nt-2∗ts)
      sigma0=bt(0)/nt/(nt-ts)
      WRITE nt,mean,sigma∗nt,sigma10∗nt,sigma0∗nt,ts, IF sigma>0 THEN SQRT(sigma) ELSE 0, IF sigma10>0 THEN SQRT(sigma10) ELSE 0, IF bt(0)>0 THEN bt(1)/bt(0) ELSE 0
    END IF
  END IF
REPEAT
IF NOT TEST THEN
  LOOP FOR col=LO TO HI
    sigma=mv(col).var
    WRITE mv(col).mean, IF sigma>0 THEN SQRT(sigma) ELSE sigma, ./.
  REPEAT; WRITE
  WRITE TO stderr mv(1).nt,mv(1).ts
END IF
#endif