Census Code

{* ***** BEGIN LICENSE BLOCK *****                                            *}
{* Version: MPL 1.1                                                           *}
{*                                                                            *}
{* The contents of this file are subject to the Mozilla Public License        *}
{* Version 1.1 (the "License"); you may not use this file except in           *}
{* compliance with the License. You may obtain a copy of the License at       *}
{* http://www.mozilla.org/MPL/                                                *}
{*                                                                            *}
{* Software distributed under the License is distributed on an "AS IS" basis, *}
{* WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License   *}
{* for the specific language governing rights and limitations under the       *}
{* License.                                                                   *}
{*                                                                            *}
{* The Original Code is Census                                                *}
{*                                                                            *}
{* The Initial Developer of the Original Code is Justin Wilkins               *}
{*                                                                            *}
{* Portions created by Justin Wilkins are Copyright (C)2001-2008              *}
{* Justin Wilkins. All Rights Reserved.                                       *}
{*                                                                            *}
{* Contributor(s):                                                            *}
{*                                                                            *}
{* ***** END LICENSE BLOCK *****                                              *}

unit loess_unit;

interface

uses
  Math;

procedure Loess(n, nsteps: Integer; delta, f: Double);
procedure Loest(x, y, w, rw: array of Double; xs, ys: Double;
  n, nleft, nright: Integer; userw, ok: Boolean);
procedure LoessSort(var List: array of Double;
  max, min: Integer);
function Find_Max(Arr: array of Double): Double;
function Find_Min(Arr: array of Double): Double;

implementation

procedure Loess(n, nsteps: Integer; delta, f: Double);

  {subroutine lowess(x, y, n, f, nsteps, delta, ys, rw, res)
      integer n
      integer nsteps
      real x(n), y(n), f, delta, ys(n), rw(n)
      real res(n)
      integer nright, min0, max0, i, j, ifix
      integer iter, last, m1, m2, ns, nleft
      real abs, cut, cmad, r, d1, d2
      real c1, c9, alpha, denom, float
      logical ok }
var
  x, y, ys, rw, res: array of Double;
  ok, userw: Boolean;
  m1, m2, ns, iter, nleft, nright, last, i, j: Integer;
  d1, d2, cmad, r, denom, alpha, cut, c1, c9: Double;
begin
  SetLength(x, n);
  SetLength(y, n);
  SetLength(ys, n);
  SetLength(rw, n);
  SetLength(res, n);

{      if (n .ge. 2) goto 1
         ys(1) = y(1)
         return
c at least two, at most n points      }

  if n < 2 then
    Exit;
 {  1  ns = max0(min0(ifix(f*float(n)), n), 2)
      iter = 1
         goto  3
   2     iter = iter+1
   3     if (iter .gt. nsteps+1) goto  22    }
  ns := Max(Min(Floor(f * n), n), 2);
  { for(iter=1; iter<=nsteps+1; iter=iter+1){      # robustness iterations
         nleft = 1; nright = ns
         last = 0        # index of prev estimated point
         i = 1   # index of current point     }
  for iter := 1 to nsteps + 1 do
  begin
    nleft := 1;
    nright := ns;
    last := 0;
    i := 1;
      {   repeat{
                 while(nright<n){
*  # move nleft, nright to right if radius decreases
                         d1 = x(i)-x(nleft)
                         d2 = x(nright+1)-x(i)
*  # if d1<=d2 with x(nright+1)==x(nright), lowest fixes
                         if (d1<=d2) break
*  # radius will not decrease by move right
                         nleft = nleft+1
                         nright = nright+1
                         }
    repeat
      while nright < n do
      begin
        d1 := x[i] - x[nleft];
        d2 := x[nright + 1] - x[i];
        if d1 <= d2 then
          Break;
        nleft := nleft + 1;
        nright := nright + 1;
      end;
            {     call lowest(x,y,n,x(i),ys(i),nleft,nright,res,iter>1,rw,ok)
*  # fitted value at x(i)
                 if (!ok) ys(i) = y(i)
*  # all weights zero - copy over value (all rw==0)
                 if (last<i-1) { # skipped points -- interpolate
                         denom = x(i)-x(last)    # non-zero - proof?
                         for(j=last+1; j<i; j=j+1){
                                 alpha = (x(j)-x(last))/denom
                                 ys(j) = alpha*ys(i)+(1.0-alpha)*ys(last)
                                 }
                       //  }
      if iter > 1 then
        userw := True
      else
        userw := False;
      Loest(x, y, res, rw, x[i], ys[i], n, nleft, nright, userw, ok);
      if ok = False then
        ys[i] := y[i];
      if last < i - 1 then
      begin
        denom := x[i] - x[last];
        for j := last + 1 to i - 1 do
        begin
          alpha := (x[j] - x[last]) / denom;
          ys[j] := (alpha * ys[i]) + ((1 - alpha) * ys[last]);
        end;
      end;
             {    last = i        # last point actually estimated
                 cut = x(last)+delta     # x coord of close points
                 for(i=last+1; i<=n; i=i+1){     # find close points
                         if (x(i)>cut) break     # i one beyond last pt within cut
                         if(x(i)==x(last)){      # exact match in x
                                 ys(i) = ys(last)
                                 last = i
                                 }
                       //  }
              //   i=max0(last+1,i-1)
   //   *  # back 1 point so interpolation within delta, but always go forward
   //             } until(last>=n)
      last := i;
      cut := x[last] + delta;
      for i := last + 1 to n do
      begin
        if x[i] > cut then
          Break;
        if x[i] = x[last] then
        begin
          ys[i] := ys[last];
          last := i;
        end;
      end;
      i := Max(last + 1, i - 1);
    until last >= n;
    {     do i = 1,n      # residuals
                 res(i) = y(i)-ys(i)
         if (iter>nsteps) break  # compute robustness weights except last time
         do i = 1,n
                 rw(i) = abs(res(i))
         call sort(rw,n)
         m1 = 1+n/2; m2 = n-m1+1
         cmad = 3.0*(rw(m1)+rw(m2))      # 6 median abs resid
         c9 = .999*cmad; c1 = .001*cmad
         do i = 1,n {
                 r = abs(res(i))
                 if(r<=c1) rw(i)=1.      # near 0, avoid underflow
                 else if(r>c9) rw(i)=0.  # near 1, avoid underflow
                 else rw(i) = (1.0-(r/cmad)**2)**2
                 }
       //  }
 // return
 // end
    for i := 1 to n do
      res[i] := y[i] - ys[i];
    if iter > nsteps then
      Break;
    for i := 1 to n do
      rw[i] := Abs(res[i]);
    LoessSort(rw, n, 0);
    m1 := Floor(1 + (n / 2));
    m2 := n - m1 + 1;
    cmad := 3 * (rw[m1] + rw[m2]);
    c9 := 0.999 * cmad;
    c1 := 0.001 * cmad;
    for i := 1 to n do
    begin
      r := Abs(res[i]);
      if r <= c1 then
        rw[i] := 1
      else
        if r > c9 then
          rw[i] := 0
        else
          rw[i] := Sqr(1 - Sqr(r / cmad));
    end;
  end;
end;

procedure Loest(x, y, w, rw: array of Double; xs, ys: Double;
  n, nleft, nright: Integer; userw, ok: Boolean);
 { subroutine lowest(x,y,n,xs,ys,nleft,nright,w,userw,rw,ok)
       integer n
      integer nleft, nright
      real x(n), y(n), xs, ys, w(n), rw(n)
      logical userw, ok
      integer nrt, j
      real abs, a, b, c, h, r
      real h1, sqrt, h9, amax1, range}
var
  nrt, j, k: Integer;
  a, b, c, h, r: Double;
  h1, h9, range: Double;
begin
 { range = x(n)-x(1)
  h = amax1(xs-x(nleft),x(nright)-xs)
  h9 = .999*h
  h1 = .001*h
  a = 0.0        # sum of weights         }

  range := x[n] - x[1];
  h := Max(xs - x[nleft], x[nright] - xs);
  h9 := 0.999 * h;
  h1 := 0.001 * h;
  a := 0;

{  for(j=nleft; j<=n; j=j+1){     # compute weights (pick up all ties on right)
         w(j)=0.
         r = abs(x(j)-xs)
         if (r<=h9) {    # small enough for non-zero weight
                 if (r>h1) w(j) = (1.0-(r/h)**3)**3
                 else      w(j) = 1.
                 if (userw) w(j) = rw(j)*w(j)
                 a = a+w(j)
                 }
        { else if(x(j)>xs)break   # get out at first zero wt on right
         }

  for j := nleft to n do
  begin
    w[j] := 0;
    r := Abs(x[j] - xs);
    if r <= h9 then
    begin
      if r > h1 then
        w[j] := IntPower(1 - IntPower((r / h), 3), 3)
      else
        w[j] := 1;
      if userw then
        w[j] := rw[j] * w[j];
      a := a + w[j];
    end
    else
      if x[j] > xs then
        Break;
  end;

{  nrt=j-1        # rightmost pt (may be greater than nright because of ties)
  if (a<=0.0) ok = FALSE
  else { # weighted least squares
         ok = TRUE
         do j = nleft,nrt
                 w(j) = w(j)/a   # make sum of w(j) == 1  }
  nrt := j - 1;
  if a <= 0 then
    ok := False
  else
    ok := True;
  for j := nleft to nrt do
    w[j] := w[j] / a;
    {     if (h>0.) {     # use linear fit
                 a = 0.0
                 do j = nleft,nrt
                         a = a+w(j)*x(j) # weighted center of x values
                 b = xs-a
                 c = 0.0
                 do j = nleft,nrt
                         c = c+w(j)*(x(j)-a)**2
                 if(sqrt(c)>.001*range) {       }
{*  # points are spread out enough to compute slope
                         b = b/c
                         do j = nleft,nrt
                                w(j) = w(j)*(1.0+b*(x(j)-a))
                         }
            //     }
  if h > 0 then
  begin
    a := 0;
    for j := nleft to nrt do
      a := a + w[j] * x[j];
    b := xs - a;
    c := 0;
    for j := nleft to nrt do
      c := c + w[j] * Sqr(x[j] - a);
    if Sqrt(c) > 0.001 * range then
    begin
      b := b / c;
      for j := nleft to nrt do
        w[j] := w[j] * (1 + b * (x[j] - a));
    end;
  end;
      {   ys = 0.0
         do j = nleft,nrt
                 ys = ys+w(j)*y(j)
         }

  ys := 0;
  for j := nleft to nrt do
    ys := ys + w[j] * y[j];
end;

procedure LoessSort(var List: array of Double;
  max, min: Integer);
var
  med_value: Double;
  hi, lo, i: Integer;
begin

    // If the list has <= 1 element, it's sorted.
  if (min >= max) then Exit;

    // Pick a dividing item randomly.
  i := min + Trunc(Random(max - min + 1));
  med_value := List[i];

    // Swap it to the front so we can find it easily.
  List[i] := List[min];

    // Move the items smaller than this into the left
    // half of the list. Move the others into the right.
  lo := min;
  hi := max;
  while (True) do
  begin
        // Look down from hi for a value < med_value.
    while (List[hi] >= med_value) do
    begin
      hi := hi - 1;
      if (hi <= lo) then Break;
    end;
    if (hi <= lo) then
    begin
            // We're done separating the items.
      List[lo] := med_value;
      Break;
    end;

        // Swap the lo and hi values.
    List[lo] := List[hi];

        // Look up from lo for a value >= med_value.
    lo := lo + 1;
    while (List[lo] < med_value) do
    begin
      lo := lo + 1;
      if (lo >= hi) then Break;
    end;
    if (lo >= hi) then
    begin
            // We're done separating the items.
      lo := hi;
      List[hi] := med_value;
      Break;
    end;

        // Swap the lo and hi values.
    List[hi] := List[lo];
  end; // while (True) do

    // Sort the two sublists.
  LoessSort(List, min, lo - 1);
  LoessSort(List, lo + 1, max);
end;

function Find_Max(Arr: array of Double): Double;
var
  i: Integer;
  M: Double;
begin
  M := Arr[Low(Arr)];
  for i := 1 to High(Arr) do
    if Arr[i] > M then M := Arr[i];
  Result := M;
end;

function Find_Min(Arr: array of Double): Double;
var
  i: Integer;
  M: Double;
begin
  M := Arr[Low(Arr)];
  for i := 1 to High(Arr) do
    if Arr[i] < M then M := Arr[i];
  Result := M;
end;

end.