/**********************************************************
*******  PAV - Pool-Adjacent-Violators  *******************

Author: Niels H. Veldhuijzen  mailto:nhveldh@universal.nl
Date:   July, 1999

This code is written and submitted for public, non-commercial use.
Although I used the procedure on many datasets successfully, 
I cannot guarantee perfect performance in all circumstances. 
Please send comments and error reports to:
        nhveldh@universal.nl
If this code is used in any project, 
appropriate attribution will be appreciated.

The PAV algorithm has been described in many places.
My reference is:
 R.E. Barlow, D.J. Bartholomew, J.M. Bremner, H.D. Brunk (1972):
 Statistical Inference under Order Restrictions.
 Wiley, Chichester. ISBN 0 471 04970 0.

The PAV procedure takes one input parameter:
a row or column vector.
It produces a column vector of the same
length as the input vector,
such that the elements in the output vector
are monotonic nondecreasing
in their rank numbers.
I.e., a vector X is mapped onto a vector X* such
that if (I gt J), (X*[I] ge X*[J]) and X* approximates X
in a least-squares sense.

The procedure PAVSUB performs one pass of the PAV-algoritm.
Let S be a vector in which each element contains
the sum of a block of order violating elements
of the original vector.
With each element of S there corresponds
an element of the vector N
which holds the number of elements in the block.
The vector X contains S./N,
the mean values of all the blocks.
Before the first pass, S is equal to the input vector X,
and N contains all ones.
The local vector D is filled with
the number of elements of X which
form an order violating block.
If all entries of D are equal to one,
there are no more offending blocks.
The number of elements in D unequal to one are counted in C; this variable is the stop criterion.

Use:             xstar=pav(x)
Input:           x, a row or column vector with k components
Output:          xstar, a column vector with k components
************************************************************
***********************************************************/
proc(4)=pavsub(s,n,x);
local d,c,t,p,i;
d=((x-shiftr(x',-1,x[rows(x)]+1)') .lt 0).*seqa(1,1,rows(x));
d=selif(d,d .ne 0);
d=d-shiftr(d',1,0)'; @ sizes of order violating blocks @
c=sumc(d-ones(rows(d),1));
if c ne 0; @ There are violating blocks @
    t={};
    p={};
    i=1;
    do while i le rows(d)-1;
        t=t|sumc(s[seqa(1,1,d[i])]); @ Form new vector of sums @
        p=p|sumc(n[seqa(1,1,d[i])]); @ Form new vector of counts @
        s=trimr(s,d[i],0);
        n=trimr(n,d[i],0);
        i=i+1;
    endo;
    s=t|sumc(s); @ Update s @
    n=p|sumc(n); @ Update n @
    x=s./n;      @ Update x, the block means @
endif;
retp(s,n,x,c);
endp;
/**********************************************************/
proc(1)=pav(x);
local s,y,n,c,i,j;
/*
Initialise s, y, n and c.
The vector s is the vector of current block sums,
n is the vector of current block sizes,
and y is the vector of current block means.
*/
x=vecr(x); @ make input a column vector @
s=x;
y=x;
n=ones(rows(x),1);
c=1; @ criterion: number of order violating blocks with more than one
element @
do while c ne 0;
    {s,n,y,c}=pavsub(s,n,y);
endo;
/*
Ready, so expand the vector y into s
by repeating each block mean y[i]
precisely n[i] times.
*/
s={};
i=1;
do while i le rows(n);
    s=s|ones(n[i],1)*y[i];
    i=i+1;
endo;
retp(s);
endp;