# The stats package for Statistics

package.checkclib('stats');  # has C dynamic link library been found ?

stats.initstring := 'stats v1.8.0 as of August 04, 2014\n';

create table stats.aux;

# mean of all the values in structure o; 6x faster than a C API implementation; 0.9.1, Jan 13, 2008;
# extended 1.3.3, January 30, 2011; extended 1.6.0, April 31, 2012
stats.mean := proc(o :: {table, sequence}) is
   local n := size(o);
   if n < 2 then return fail fi;  # 2.2.3 change
   return sadd o/n
end;


# quadratic mean of all the values in structure o; as fast as a C API implementation; 0.9.1, Jan 13, 2008;
# extended 1.3.3, January 30, 2011; extended 1.6.0, April 31, 2012
stats.qmean := proc(o :: {table, sequence}) is
   local n := size(o);
   if n < 2 then return fail fi;  # 2.2.3 change
   return sqrt(qsadd o/n)
end;


# The function computes the geometric mean of the values in the table or sequence o. It is a measure of
# central tendency. 14.03.2012; extended 1.6.0, April 31, 2012; extended 2.2.0, May 09, 2014
stats.gmean := proc(o :: {table, sequence}, option) is
   local r, n := 1, size o;
   if n < 2 then return fail fi;  # 2.2.3 change
   if option then
      r := smul(o)
   else
      for i in o do
         if i :: number then
            if i < 0 then  # switch into `complex mode` like Maple V Release 4 does
               mul r, i!0
            else
               mul r, i
            fi
         else
            argerror(i, 'stats.gmean', 'number expected in structure')
         fi
      od
   fi;
   return r^(1/n)
end;


# The function computes the harmonic mean of the values in the table or sequence o. It is useful with rates
# and ratios, as it provides the truest average. 1.6.0, March 31, 2012
stats.hmean := proc(o :: {table, sequence}) is
   local r, n := 0, size o;
   if n < 2 then return fail fi;  # 2.2.3 change
   for i in o do
      if i :: number then
         r := r + 1/i
      else
         argerror(i, 'stats.hmean', 'number expected in structure')
      fi
   od;
   return n/r
end;


# Variance of all the values in table t; as fast as a C API implementation; 0.9.1, Jan 13, 2008,
# tuned January 21, 2008; extended 1.3.3, January 30, 2011; extended 1.6.0, March 31, 2012
stats.var := proc(o :: {table, sequence}) is
   local m, n;
   n := size o;  # 2.2.3 change
   if n < 2 then return fail fi;
   m := sadd o/n;
   return (qsadd o - n * m*m) / n
end;


# Standard deviation of all the values in table t; 0.9.1, Jan 13, 2008
# tuned January 21, 2008; extended 1.3.3, January 30, 2011; extended 1.6.0, March 31, 2012; extended
# 2.2.0, May 10, 2014
stats.sd := proc(o :: {table, sequence}, option) is
   local m, n, r;
   n := size o;
   if n < 2 then return fail fi;  # 2.2.3 change
   m := sadd o/n;
   r := sqrt((qsadd o - n * m*m) / n);
   if option then
      return r / abs(m)
   else
      return r
   fi
end;


# Sample standard deviation of all the values in table t; 1.11.2
stats.ssd := proc(o :: {table, sequence}) is
   local m, n;
   n := size o;
   if n < 2 then return fail fi;
   m := sadd o/n;
   return sqrt((qsadd o - n * m*m) / (n-1))
end;


# Returns all values in the given table or sequence with the largest number of occurrence. Agena 1.3.3, January 30, 2011
stats.mode := proc(o) is
   local t, max, result;
   create dict t(size o);
   if o :: table then
      create table result
   elif o :: sequence then
      create seq result
   else
      argerror(o, 'stats.mode', 'table or sequence expected')
   fi;
   max := -infinity;
   for i in o do
      if t[i] then
         inc t[i], 1
      else
         t[i] := 1
      fi;
      if t[i] > max then
         max := t[i]
      fi
   od;
   for i, j in t do
      if max = j then
         insert i into result
      fi
   od;
   return result
end;


# Returns the first, second, and third quartile of a _sorted_ table or sequence. Agena 1.3.3, January 30, 2011.
# See also: http://www.vias.org/tmdatanaleng/cc_quartile.html.
# stats.quartiles([-20, -4, -1, -1, 0, 1, 2, 3, 4, 4, 7, 7, 8, 11, 11, 12, 12, 15, 18, 22]) -> [0, 5.5, 12]
stats.quartiles := proc(o, p) is
   local r, n, l1, u1, iqr, a, b, m, pos25, pos75;
   if o :: table then
      if unassigned(p) then create table r(3) fi
   elif o :: sequence then
      if unassigned(p) then create sequence r(3) fi
   else
      argerror(o, 'stats.quartiles', 'table or sequence expected')
   fi;
   n := size o;
   if n < 2 then
      error('Error in `stats.quartiles`: structure must contain at least two elements.')  # 1.7.8 u 1
   fi;
   pos25, pos75 := roundf(0.25*(n + 1)), roundf(0.75*(n + 1));  # first and third quartile determined according to NIST
   if p :: number then
      case p
         of 1 then return o[pos25]
         #of 2 then return stats.median(o)
         of 2 then return if not even(n) then o[n\2 + 1] else (o[n\2] + o[n\2 + 1])/2 fi
         of 3 then return o[pos75]
         else
            error('Error in `stats.quartiles`: numbers 1, 2, or 3 expected as 2nd argument, got ' & p & '.')
      esac
   elif p :- 'null' then
      argerror(p, 'stats.quartiles', 'second argument is not a number')
   fi;
   r[1] := o[pos25];
   r[2] := if not even(n) then o[n\2 + 1] else (o[n\2] + o[n\2 + 1])/2 fi;  #stats.median(o);
   r[3] := o[pos75];
   if n > 4 then
      # determine interquartile range
      iqr := r[3] - r[1];
      # compute lower outlier limit L1
      l1 := r[1] - 1.5*iqr;
      a := 1; b := pos25 - 1;
      while a < b do
         m := (a + b) \ 2;
         #print('>>', a, o[a], b, o[b], m, o[m], o[m] > l1)
         if l1 > o[m] then
            a := m + 1
         else
            b := m
         fi
      od;
      try
         while o[m] <= l1 do inc m od;  # evade round-off errors
         r[4] := o[m]  # use index `m` instead of a or b !
      catch
         return r
      yrt;
      # compute upper outlier limit U1
      u1 := r[3] + 1.5*iqr;
      a := pos75 + 1; b := n;
      while a < b do
         m := (a + b) \ 2;
         #print('<<', a, o[a], b, o[b], m, o[m], o[m] < u1)
         if u1 > o[m] then
            a := m + 1
         else
            b := m
         fi
      od;
      try
         while o[m] >= u1 do dec m od;  # evade round-off errors
         r[5] := o[m];  # use index `m` instead of a or b !
      catch
         return r
      yrt;
      r[6] := iqr
   fi;
   return r
end;


# Returns the value below which a certain percent p of the elements in o fall, in other words: the
# p-th percentile is the score below which p percent of observations may be found.
stats.aux.percentile := proc(o, rank, n) is
   case rank
      of 1 then return o[1]
      of n then return o[n]
      else
         local k, d := entier(rank), frac(rank);
         if k + 1 > n then
            error('Error in `stats.percentile`: rank is greater than number of elements.')
         fi;
         return o[k] + d*(o[k+1] - o[k])
   esac
end;

stats.percentile := proc(o :: {table, sequence}, p :: number, option) is
   # o the list, p (an integer) the p-th percentile with 0 <= p < 100; 17.02.2012
   local rank, n;
   if float(p) then
      error('Error in `stats.percentile`: second argument must be an integer, got a float.')
   elif p < 0 or p > 99 then
      error('Error in `stats.percentile`: second argument must be in the range [0, 100).')
   elif option <> null and option :- string then
      error('Error in `stats.percentile`: option must be either \'nist\' or \'excel\'.')
   fi;
   n := size(o);
   if n = 0 then
      error('Error in `stats.percentile`: first argument must contain at least one element.')
   fi;
   case option
      of null then  # no option given: default definition as decribed in the English Wikipedia article `Percentile`
         rank := p/100*n + 1/2;
         return o[roundf(rank)]  # round rank to the nearest integer
      of 'nist' then   # Engineering Statistics Handbook: Percentile, NIST websites
         rank := p/100*(n+1);
         return stats.aux.percentile(o, rank, n)
      of 'excel' then  # as described in Wikipedia, "Pottel, Hans", Statistical flaws in Excel"
         rank := p/100*(n-1) + 1;
         return stats.aux.percentile(o, rank, n)
      else
         error('Error in `stats.percentile`: unknown third option \'' & option & '\'.');
   esac
end;


stats.aux.check_range := proc(a, b, procname :: string) is
   if assigned(a) then
      if a :- number then
         argerror(a, procname, 'second argument must be an integer')
      elif float(a) then
         error('Error in `stats.' & procname & '`: second argument must be an integer, got a float.')
      elif a < 0 or a > 99 then
         error('Error in `stats.' & procname & '`: second argument must be in the range [0, 100).')
      fi
   else
      a := 25
   fi;
   if assigned(b) then
      if b :- number then
         argerror(b, procname, 'third argument must be an integer')
      elif float(b) then
         error('Error in `stats.' & procname & '`: third argument must be an integer, got a float.')
      elif b < 0 or b > 99 then
         error('Error in `stats.' & procname & '`: third argument must be in the range [0, 100).')
      fi
   else
      b := 100 - a
   fi;
  return a, b
end;


# number of permutations; 14.03.2012, extended 18.04.2013
stats.numbperm := proc(x, y :: number) is
   if x :: set then
      return stats.numbperm(size x, y)
   elif x :- number then
      argerror(x, 'stats.numbperm', 'expected an integer or a set')
   fi;
   if float(x) or float(y) or x < 0 or y < 0 then  # faster than calling not(isnonnegint)
      return undefined
   elif x < y then
      return 0
   else
      return exp(lngamma(x+1) - lngamma(x-y+1))
   fi
end;


# Counts the combinations of the elements of n taken k at a time and returns a number. If n or m are
# negative or floats, the function returns `undefined`. 1.11.2 u1
stats.numbcomb := proc(x, y :: number) is
   if x :: set then
      return stats.numbcomb(size x, y)
   elif x :- number then
      argerror(x, 'stats.numbcomb', 'expected an integer or a set')
   fi;
   return binomial(x, y)
end;


# Sample skewness, a measure of the asymmetry of the probability distribution; returns 0 if
# a distribution is symmetric; a negative value if the left tail is longer; a positive value
# the right tail is longer. It computes the third moment about the mean and divides it by
# the third power of the standard deviation. See article on `Skewness` in the English Wikipedia.
# Multiply the result by sqrt(n*(n-1))/(n-2) [Wikipedia], or sqrt(n*(n-1))/n [Maple V R4], where
# n is the number of elements in the distribution, to get the population skewness. The latter
# however may return wrong results.
# 11.03.2012
stats.skewness := proc(o :: {table, sequence}) is
   local m, n, sd;
   n := size o;
   if n < 2 then return fail fi;
   sd := stats.sd(o)**3;
   if sd = 0 then return fail fi;
   m := sadd o/n;
   return stats.moment(o, 3, m) / sd
end;


# Interquartile Range
# determines the difference of the third and first quartile and is useful for determining
# the variability in a distribution o.
stats.iqr := proc(o :: {table, sequence}, a, b) is
   local s;
   #/ if not(o :: table or o :: sequence) then
      argerror(o, 'stats.iqr', 'first argument must be a table or sequence')
   fi; /#
   a, b := stats.aux.check_range(a, b, 'stats.iqr');
   if not stats.issorted(o) then
      s := copy(o);
      sort(s)
   else
      s := o
   fi;
   return stats.percentile(s, b) - stats.percentile(s, a)
end;


# counts the number of occurrences of each entry in a table or structure and returns a dictionary
# with its respective key the entry and its value the number of occurrences. You might optionally
# pass a procedure to be mapped on the structure before counting begins on the thus modified
# structure. Extended Agena 1.8.3, 08.10.2012.
stats.countentries := proc(o, f, ?) is
   create local dict result;
   if not(o :: table or o :: sequence) then
      argerror(o, 'stats.countentries', 'first argument must be a table or sequence')
   fi;
   if assigned(f) then
      if f :: procedure then
         o := map(f, o, unpack(varargs))
      else
         argerror(f, 'stats.countentries', 'second argument must be a procedure')
      fi
   fi;
   for i in toset(o) do
      result[i] := countitems(i, o)
   od;
   return result
end;

# Divides a numeric range defined by the pair p and its step size n into its subintervals, sorts all
# occurrences in the observation s (a sequence) into these subranges and finally counts all elements in
# these subranges.

# The function returns a table with the keys the respective left borders of the subranges and the values
# the number of counts in the respective subranges. It always also returns a second table which may include
# all those elements in s which are not part of the overall range defined by p.

# If an element in s equals the right border of a subinterval, then it is considered to be part of the next
# subinterval. But if an element in s equals the right border of the overall intervall p, it is considered
# part of the last subinterval.

# The function issues an error if it encounters a non-number in s, or if the left border in p is greater or
# equals to the right border in p.
#
# An example:
# > s := seq(0.1, 0.2, 0.3, 0.4, 1, 1.1, 2, 2.1);
#
# > stats.obcount(s, 0:2, 1):
# [0 ~ 4, 1 ~ 3]  [2.1]
#
# Algorithm for bisection taken from Niklaus Wirth, `Algorithmen und Datenstrukturen mit Modula-2`,
# 4th ed., 1986, p. 58
stats.obcount := proc(x :: sequence, p :: pair, n :: number) is  # Agena 1.7.8 u2, 03.09.2012
   local a, b, k, y, result, errors, ny;  # ny: 2.2.1 fix
   a, b := left(p), right(p);
   if a :- number or b :- number then
      error('Error in `stats.obcount`: pair must consist of numbers.')
   elif a >= b then
      error('Error in `stats.obcount`: left border is greater or equal to the right border.')
   fi;
   y := nseq(a, b, n);  # create all subintervals
   ny := size y;
   create table result(ny), table errors;
   # initialise with zero occurrences
   for i from a to b by n do
      result[i] := 0
   od;
   for i in x do
      if i :- number then
         error('Error in `stats.obcount`: ' & tostring(i) & ' is not a number, got ' & type(i))
      elif not(i in p) then  # ignore number
         insert i into errors
      else
         k := entier((i - a + n) / n) - abs(i = b);  # 1.12.5/6 improvement/patch
         inc result[y[k]]
      fi
   od;
   return result, errors
end;


# Sorts occurrences into subintervals, based on stats.obcount; Agena 1.12.4 u2, 07.08.2013
stats.obpart := proc(x :: sequence, p :: pair, n :: number, criterion, mapping) is
   local a, b, y, result, errors, f, mappinggiven, fni, k, ny;
   a, b := left(p), right(p);
   if a :- number or b :- number then
      error('Error in `stats.obpart`: pair must consist of numbers.')  # 1.12.9
   elif a >= b then
      error('Error in `stats.obpart`: left border is greater or equal to the right border.')
   fi;
   y := nseq(a, b, n);  # create all subintervals
   ny := size y;
   create table result(ny), table errors;
   # initialise sequences for the occurrences
   #for i from a to b - n by n do
   for i from a to b by n do
      result[i] := seq()
   od;
   if criterion then
      if criterion :- procedure then
         argerror(criterion, 'stats.obpart', 'expected a procedure for fourth argument')
      fi;
      f := criterion
   else
      f := << x -> x >>
   fi;
   if mapping then
      if mapping :- procedure then
         argerror(mapping, 'stats.obpart', 'expected a procedure for fifth argument')
      fi;
      mappinggiven := true
   fi;
   for i in x do
      fni := f(i);  # 1.12.4 u3, speed up computation
      if fni :- number then
         error('Error in `stats.obpart`: ' & tostring(fni) & ' is not a number.')  # 1.12.9
      elif not(fni in p) then  # ignore number
         insert i into errors
      else
         k := entier((fni - a + n) / n) - abs(fni = b);  # 1.12.5/6 improvement/patch
         insert if mapping then mapping(i) else i fi into result[y[k]]
      fi
   od;
   return result, errors
end;


# Computes the cumulative density function between the lower bound a and the upper bound
# b. If the mean mu is not given, it defaults to 0; if the standard deviations sigma is not given,
# it defaults to 1.
stats.cdf := proc(a :: number, b :: number, mu, sigma) is  # 1.11.0
   if unassigned(calc) then import calc fi;
   if unassigned(mu) then mu := 0 fi;
   if unassigned(sigma) then sigma := 1 fi;
   if sigma <= 0 then
      error('Error in `stats.cdf`: fourth argument must be positive.')
   fi;
   if a > b then
      error('Error in `stats.cdf`: second argument must be greater than fist argument.')
   fi;
   return stats.ndf(sigma) * calc.intde(<< x -> stats.nde(x, mu, sigma) >>, a, b)
end;


# stats.trimmean
#/ Returns the arithmetic mean of the interior of a distribution o (of type table or sequence), where the number f el [0, 1)
  determines the fraction of the data that is to be excluded from the margins. The number p of data to be excluded is always
  rounded down to the nearest even number. The function then does not take into account p/2 points from the left margin and
  p/2 points from the right margin when calculating the average using Kahan/Ozawa round-off error prevention. The function
  does not sort the observation.

  The return is a number. If the structure is empty, `fail` is returned. /#

stats.trimmean := proc(o :: {table, sequence}, f :: number) is
   local n, a;
   n := size o;
   if n < 2 then
      return fail
   elif f < 0 or f >= 1 then
      error('Error in `stats.trimmean`: factor must be in the range [0, 1).')
   fi;
   a := int(f * n);
   if not even(a) then dec a fi;
   div a, 2;
   return stats.amean(o[1 + a to n - a])
end;


# Returns the normalised Herfindahl–Hirschman index, an indicator of the amount of competition in economy.
# A value of 0 means that there is absolute competition, i.e. that all companies have the same share, and
# 1 means that there is a monopoly.
#
# It is also a good measure to determine the stability of a distribution, with a value tending to zero indicating
# that the number of outliers is low, and a value tending to 1 that there are outliers.

stats.herfindahl := proc(o :: {table, sequence}) is
   local m, rs, h;
   for i in o do
      if not isnonneg(i) then
         error('Error in `stats.herfindahl`: all elements in structure must be non-negative.')
      fi
   od;
   m := sadd(o);
   rs := recip size o;
   h := sadd(<< x -> (x/m)**2 >> @ o);  # compute Herfindahl-Hirschmann index
   return (h - rs)/(1 - rs)  # normalise it
end;


#/ stats.chauvenet

Receives a table or sequence obj of numbers and checks them for outliers using the formula:

   p := n * erfc( (abs(x - stats.amean(obj)) / stats.sd(obj)) )

where n is the number of observations in a distribution, and x a sample of it.

In the first form, it checks whether the number x is an outlier and returns true or false.

In the second form, it iterates obj again and again as long as it does not find an outlier,
and returns the outliers in a structure defined by the type of obj.

By default, if p < 0.5, where 0.5 is the magical Chauvenet number, an outlier is detected.
If you pass the option bailout=c, then c, a nonnegative number, will be the threshold.

If you pass the option jump=true, as soon as an outlier is detected, it is removed from the
distribution and then the whole evaluation process is restarted immediately with a reduced
distribution along with a recomputed mean and deviation.

If you do not, all remaining items are also checked according to the current criteria - after
the last item has been checked, only then the outliers are removed from the distribution, the
mean and deviation are recomputed and another iteration begins.

Further information: `Cleaning Data the Chauvenet Way`, by Lily Lin and Paul D. Sherman,
published at the South East SAS Users Group's website http://www.sesug.org. /#

stats.chauvenet := proc(s :: {table, sequence}, ?) is  # initiated Agena 2.2.0 RC 4
   local r, probab, n, dev, x, m, bailout, jump, meanfn, devfn, outlier;
   bailout := 0.5;  # Chauvenet's magical number
   meanfn := stats.amean; devfn := stats.sd;
   for i in varargs do
      if i :: pair then
         local lhs, rhs := left(i), right(i);
         if lhs = 'bailout' and rhs :: number and rhs >= 0 then  # use a threshold different from 0.5
            bailout := rhs
         elif lhs = 'jump' and rhs :: boolean then
            jump := rhs   # true = immediately stop current iteration, delete outlier, update all data and
                          # recommence evaluation
         elif lhs = 'mean' and rhs :: procedure then  # 2.2.3, experimental
            meanfn := rhs
         elif lhs = 'dev' and rhs :: procedure then  # 2.2.3, experimental
            devfn := rhs
         elif lhs = 'outlier' and rhs in {'lower', 'upper'} then  # 2.2.6, return only lower or upper outliers
            outlier := rhs
         else
            error('Error in `stats.chauvenet`: unknown option `' & lhs & '`.')
         fi
      elif i :: number then  # check a specific sample
         x := i
      else
         error('Error in `stats.chauvenet`: invalid argument.')
      fi
   od;
   n := size s;  # 2.2.3 fix
   if n < 2 then return fail fi;
   if unassigned(x) then  # no specific sample to be checked but entire distribution ?
      local o, m, sd, flag, probab;
      if s :: sequence then
         create sequence r
      else
         create table r
      fi;
      o := copy(s);  # copy the distribution so that the function acts non-destructively
      do
         m := meanfn(o);
         sd := devfn(o);
         for i in copy(o) do  # copy, for removing outliers from a structure during an iteration may cause undefined behaviour
            flag := true;
            probab := n * erfc( (abs(i - m) / sd) );
            if probab < bailout then  # outlier ?
               skip when assigned(outlier) and (outlier = 'lower' and i >= m or outlier = 'upper' and i < m);
               insert i into r;
               delete i from o;  # delete all occurrences of observation i from distribution
               flag := false;
               break when jump   # immediately restart estimation and recalculate size n, amean, and standard deviation
            fi
         od;
         n := size o
      until flag or n < 2  # traversed whole observation but did not find any further outlier or observation too small ?
   else  # check specific observation only
      m := meanfn(s);
      probab := n * erfc( (abs(x - m) / devfn(s)) );
      r :=  probab < bailout and
         (unassigned(outlier) or (outlier = 'lower' and x < m or outlier = 'upper' and x > m))
   fi;
   return r  # return outliers or the result of the test of a specific observation
end;


# stats.extrema
# Expects a sequence or table obj of points xk:yk and the number delta and determines the local minima and maxima.
# A value yk is considered an extrema if the difference to its surrounding is at least delta.
# The function returns two structures of pairs, i.e. points, the first one including the local minima, the second
# one the local maxima. The type of the structures is determined by the type of obj.
# This is the Agena port of the `peakdet` function written by Eli Billauer for MATLAB, see: http://billauer.co.il/peakdet.html.

stats.extrema := proc(v :: {sequence, table}, delta :: number) is  # 2.2.0 u1
   local maxtab, mintab, min, max, minpos, maxpos, this, x, y, lookformax;
   if v :: table then
      create table maxtab, table mintab
   else
      create sequence maxtab, sequence mintab
   fi;
   if delta <= 0 then
      error('Error in `stats.extrema`: second argument must be positive.')
   fi;
   min := infinity; max := -infinity;
   lookformax := true;
   for i in v do
      if i :- pair then
         error('Error in `stats.extrema`: first argument must include pairs only.')
      fi;
      x := left(i); y := right(i);
      if x :- number or y :- number then
         error('Error in `stats.extrema`: pairs must consist of numbers.')
      fi;
      if y > max then
         max := y; maxpos := x
      fi;
      if y < min then
         min := y; minpos := x
      fi;
      if lookformax then
         if y < max - delta then
            insert maxpos:max into maxtab;
            min := y; minpos := x;
            lookformax := false
         fi
      else
         if y > min + delta then
            insert minpos:min into mintab;
            max := y; maxpos := x;
            lookformax := true
         fi
      fi
   od;
   return mintab, maxtab
end;


# Returns a function computing the z-score (standard score) of a sample x in the table or sequence obj - the number of standard deviations x is above or below the mean according to the formula: z = (x - mu)/sd.

# The resulting function returns a positive number z if x is above the mean and a negative number z if it is below. It does, however, not check whether x is part of obj. The result is computed using Kahan-Ozawa round-off error prevention for mu and sd.

stats.zscore := proc(o:: {table, sequence}) is  # 2.2.7
   local mean, sd;
   mean := stats.amean(o);
   sd := sqrt(stats.kosumdata(o, 2, mean) / size o);
   return proc(x) is
      return (x - mean)/sd
   end
end;


# A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise
# See: http://www.csse.uwa.edu.au/~pk/research/matlabfns/Misc/dbscan.m and
# http://www.dbs.ifi.lmu.de/Publikationen/Papers/KDD-96.final.frame.pdf

#/
% DBSCAN DBSCAN clustering algorithm, Matlab implementation
%
% Usage:  [C, ptsC, centres] = dbscan(P, E, minPts)
%
% Arguments:
%         P - dim x Npts array of points.
%         E - Distance threshold.
%    minPts - Minimum number of points required to form a cluster.
%
% Returns:
%         C - Cell array of length Nc listing indices of points associated with
%             each cluster.
%      ptsC - Array of length Npts listing the cluster number associated with
%             each point.  If a point is denoted as noise (not enough nearby
%             elements to form a cluster) its cluster number is 0.
%   centres - dim x Nc array of the average centre of each cluster.

% Reference:
% Martin Ester, Hans-Peter Kriegel, Jörg Sander, Xiaowei Xu (1996). "A
% density-based algorithm for discovering clusters in large spatial databases
% with noise".  Proceedings of the Second International Conference on Knowledge
% Discovery and Data Mining (KDD-96). AAAI Press. pp. 226-231.
% Also see: http://en.wikipedia.org/wiki/DBSCAN

% Copyright (c) 2013 Peter Kovesi
% Centre for Exploration Targeting
% The University of Western Australia
% peter.kovesi at uwa edu au
%
% 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, 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.

% PK January 2013

function [C, ptsC, centres] = dbscan(P, E, minPts)

    [dim, Npts] = size(P);

    ptsC  = zeros(Npts,1);
    C     = {};
    Nc    = 0;               % Cluster counter.
    Pvisit = zeros(Npts,1);  % Array to keep track of points that have been visited.

    for n = 1:Npts
       if ~Pvisit(n)                            % If this point not visited yet
           Pvisit(n) = 1;                       % mark as visited
           neighbourPts = regionQuery(P, n, E); % and find its neighbours

           if length(neighbourPts) < minPts-1  % Not enough points to form a cluster
               ptsC(n) = 0;                    % Mark point n as noise.

           else                % Form a cluster...
               Nc = Nc + 1;    % Increment number of clusters and process
                               % neighbourhood.

               C{Nc} = [n];    % Initialise cluster Nc with point n
               ptsC(n) = Nc;   % and mark point n as being a member of cluster Nc.

               ind = 1;        % Initialise index into neighbourPts array.

               % For each point P' in neighbourPts ...
               while ind <= length(neighbourPts)

                   nb = neighbourPts(ind);

                   if ~Pvisit(nb)        % If this neighbour has not been visited
                       Pvisit(nb) = 1;   % mark it as visited.

                       % Find the neighbours of this neighbour and if it has
                       % enough neighbours add them to the neighbourPts list
                       neighbourPtsP = regionQuery(P, nb, E);
                       if length(neighbourPtsP) >= minPts
                           neighbourPts = [neighbourPts  neighbourPtsP];
                       end
                   end

                   % If this neighbour nb not yet a member of any cluster add it
                   % to this cluster.
                   if ~ptsC(nb)
                       C{Nc} = [C{Nc} nb];
                       ptsC(nb) = Nc;
                   end

                   ind = ind + 1;  % Increment neighbour point index and process
                                   % next neighbour
               end
           end
       end
    end

    % Find centres of each cluster
    centres = zeros(dim,length(C));
    for n = 1:length(C)
        for k = 1:length(C{n})
            centres(:,n) = centres(:,n) + P(:,C{n}(k));
        end
        centres(:,n) = centres(:,n)/length(C{n});
    end

end % of dbscan

#/
%------------------------------------------------------------------------
% Find indices of all points within distance E of point with index n
% This function could make use of a precomputed distance table to avoid
% repeated distance calculations, however this would require N^2 storage.
% Not a big problem either way if the number of points being clustered is
% small.   For large datasets this function will need to be optimised.

% Arguments:
%              P - the dim x Npts array of data points
%              n - Index of point of interest
%              E - Distance threshold

function neighbours = regionQuery(P, n, E)

    E2 = E^2;
    [dim, Npts] = size(P);
    neighbours = [];

    for i = 1:Npts
        if i ~= n
            % Test if distance^2 < E^2
            v = P(:,i)-P(:,n);
            dist2 = v'*v;
            if dist2 < E2
               neighbours = [neighbours i];
            end
        end
    end

end % of regionQuery /#


# The functions finds clusters in a sequence of n-dimensional points represented as pairs x:y (2D) or sequences
# (n-dimensional) and returns a table with the individual clusters along with their respective points.

# It also returns a register of the size of the whole distribution listing the cluster number associated with
# each point, where the point in this case is represented by its integral position in the sequence.

# eps is the maximum allowed distance between two points that shall belong to the same neighbourhood. minPts is
# the minimum number of points that shall constitute a neighbourhood.

# By specifying the 'select' option along with a function returning a boolean, e.g. 'select':<< x -> right x < 1 >>,
# only points satisfying the given criterion are examined. By specifying the 'method' option, you can control how
# the function determines clusters: 'method':'original' uses the classic one, 'method':'modified' uses a much faster
# and memory-saving implementation that contrary to the original method immediately flags neighbours of neighbours
# as being visited and thus does not examine them again in further passes. The default is 'original'.

stats.dbscan := proc(P :: sequence, eps :: number, minPts :: number, ?) is
   local npts, ptsC, C, nc, Pvisit, neighbourPts, neighbourPtsP, centres, nb, ind, sneighbourPts,
      method, selector;
   method := 'original';
   selector := << x -> true >>;  # consider all points by default
   npts := size P;
   return when npts = 0;
   for i in varargs do
      if i :: pair then
         case left i
            of 'method' then
               if right(i) in {'original', 'modified'} then
                  method := right i
               else
                  error('Error in `stats.dbscan`: invalid value `' & tostring(right i) & '` for method option.')
               fi
            of 'select' then
               if (right i) :: procedure then
                  selector := right i
               else
                  error('Error in `stats.dbscan`: procedure expected for select option, got ' & type(right i) & '.')
               fi
            else
               error('Error in `stats.dbscan`: unknown option `' & tostring(left i) & '`.')
         esac
      fi
   od;
   ptsC := nreg(<< x -> 0 >>, 1, npts);    # we create a zero register to track noise and cluster elements
   create table C;
   nc := 0;  # Cluster counter.
   Pvisit := nreg(<< x -> 0 >>, 1, npts);  # register to keep track of points that have been visited
   for n to npts do
      if Pvisit[n] = 0 and selector(P[n]) then  # If this point not visited yet and satisfies selection criterium,
         Pvisit[n] := 1;                   # mark as visited
         neighbourPts := stats.neighbours(P, n, eps, 2, true);  # and find its neighbours
         if size neighbourPts >= minPts then  # Enough points to form a cluster, then form a cluster, otherwise ignore point as noise
            inc nc;                        # Increment number of clusters and process neighbourhood.
            C[nc] := seq(P[n]);            # Initialise cluster nc with point n
            ptsC[n] := nc;                 # and mark point n as being a member of cluster nc.
            ind := 1;                      # Initialise index into neighbourPts array.
            while ind <= size neighbourPts do  # For each point P' in neighbourPts ...
               nb := neighbourPts[ind];
               if Pvisit[nb] = 0 then      # If this neighbour has not yet been visited,
                  Pvisit[nb] := 1;         # mark it as visited.
                  # Find the neighbours of this neighbour and if it has enough neighbours add them to the neighbourPts list
                  skip when not selector(P[nb]);  # ignore point if it does not satisfy criterion
                  neighbourPtsP := stats.neighbours(P, nb, eps, 2, true);  # get the indices
                  if size neighbourPtsP >= minPts then
                     if method = 'modified' then
                        for i in neighbourPtsP do
                           if Pvisit[i] = 0 then  # not yet visited ?
                              Pvisit[i] := 1;     # mark as visited
                              skip when not selector(P[i]);  # ignore point if it does not satisfy criterion
                              insert i into neighbourPts  # insert neighbour of neigbour in P-neighbour list and do not evaluate again
                           fi
                        od;
                        # avoid out-of-memory errors, duplicates have been found in real observations
                        neighbourPts := unique(neighbourPts)
                     else
                        if size neighbourPtsP >= minPts then  # avoid out-of-memory errors
                           neighbourPts := unique(neighbourPts union neighbourPtsP)
                        fi
                     fi
                  fi
               fi;
               # If this neighbour nb not yet a member of any cluster add it to this cluster if neighbour satisfies criterion
               if ptsC[nb] = 0 and selector(P[nb]) then
                  insert P[nb] into C[nc];
                  ptsC[nb] := nc
               fi;
               inc ind  # Increment neighbour point index and process next neighbour.
            od
         fi
      fi
   od;
   for i to size C do
      if C[i] then
         C[i] := unique(C[i]);  # sometimes, there are multiple occurences of a point in a cluster
         sort(C[i], << x, y -> left x < left y>>)
      fi
   od;
   # Find centres of each cluster. Contrary to the Mathlab version, we assume we are in 2D space only
   #centres := nreg(<< x -> 0 >>, 1, size C);
   #for n to size C do
   #   for k to size C[n] do
   #      #centres(:,n) = centres(:,n) + P(:,C{n}(k));
   #      centres[n] := centres[n] + P(:,C{n}(k));
   #   od;
   #   # centres(:,n) = centres(:,n)/length(C{n});
   #   centres(:,n) = centres(:,n)/length(C{n});
   #od;
   # C - Cell array of length nc listing indices of points associated with each cluster.
   # ptsC - Array of length npts listing the cluster number associated with each point.  If a point is denoted as noise
   #     (not enough nearby elements to form a cluster) its cluster number is 0.
   # centres - dim 2 nc array of the average centre of each cluster.
   return C, ptsC; # centres
end;

