basho_stats

%% -------------------------------------------------------------------
%%
%% stats: Statistics Suite for Erlang
%%
%% Copyright (c) 2009 Dave Smith (dizzyd@dizzyd.com)
%%
%% This file is provided to you under the Apache License,
%% Version 2.0 (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.apache.org/licenses/LICENSE-2.0
%%
%% Unless required by applicable law or agreed to in writing,
%% software distributed under the License is distributed on an
%% "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
%% KIND, either express or implied.  See the License for the
%% specific language governing permissions and limitations
%% under the License.
%%
%% -------------------------------------------------------------------
-module(basho_stats_histogram).

-export([new/3,
         update/2, update_all/2,
         quantile/2,
         counts/1,
         observations/1,
         summary_stats/1]).

-ifdef(EQC).
-export([prop_count/0, prop_quantile/0]).
-endif.

-include("stats.hrl").

-record(hist, { n = 0,
                min,
                max,
                bin_scale,
                bin_step,
                bins,
                capacity,
                stats }).

%% ===================================================================
%% Public API
%% ===================================================================

new(MinVal, MaxVal, NumBins) ->
    #hist { min = MinVal,
            max = MaxVal,
            bin_scale = NumBins / (MaxVal - MinVal),
            bin_step = (MaxVal - MinVal) / NumBins,
            bins = gb_trees:empty(),
            capacity = NumBins,
            stats = basho_stats_sample:new() }.


%%
%% Update the histogram with a new observation.
%%
%% NOTE: update/2 caps values within #hist.min and #hist.max;
%% if you provide a value outside those boundaries the first or last
%% bin, respectively, get updated and the histogram is consequently
%% skewed.
%%
update(Value, Hist) ->
    Bin = which_bin(Value, Hist),
    case gb_trees:lookup(Bin, Hist#hist.bins) of
        {value, Counter} ->
            ok;
        none ->
            Counter = 0
    end,
    Hist#hist { n = Hist#hist.n + 1,
                bins = gb_trees:enter(Bin, Counter + 1, Hist#hist.bins),
                stats = basho_stats_sample:update(Value, Hist#hist.stats)}.


update_all(Values, Hist) ->
    lists:foldl(fun(Value, H) -> update(Value, H) end,
                Hist, Values).
    


%%
%% Estimate the quantile from the histogram. Quantile should be a value
%% between 0 and 1. Returns 'NaN' if the histogram is currently empty.
%%
quantile(_Quantile, #hist { n = 0 }) ->
    'NaN';
quantile(Quantile, Hist)
  when Quantile > 0; Quantile < 1 ->
    %% Sort out how many complete samples we need to satisfy the requested quantile
    MaxSamples = Quantile * Hist#hist.n,

    %% Now iterate over the bins, until we have gathered enough samples
    %% to satisfy the request. The resulting bin is an estimate.
    Itr = gb_trees:iterator(Hist#hist.bins),
    case quantile_itr(gb_trees:next(Itr), 0, MaxSamples) of
        max ->
            Hist#hist.max;
        EstBin ->
            %% We have an estimated bin -- determine the lower bound of said
            %% bin
            Hist#hist.min + (EstBin / Hist#hist.bin_scale)
    end.

%%
%% Get the counts for each bin in the histogram
%%
counts(Hist) ->
    [bin_count(I, Hist) || I <- lists:seq(0, Hist#hist.capacity-1)].


%%
%% Number of observations that are present in this histogram
%%
observations(Hist) ->
    Hist#hist.n.

%%
%% Return basic summary stats for this histogram
%%
summary_stats(Hist) ->
    basho_stats_sample:summary(Hist#hist.stats).


%% ===================================================================
%% Internal functions
%% ===================================================================

which_bin(Value, Hist) ->
    Bin = trunc((Value - Hist#hist.min) * Hist#hist.bin_scale),
    Lower = Hist#hist.min + (Bin * Hist#hist.bin_step),
    Upper = Hist#hist.min + ((Bin + 1) * Hist#hist.bin_step),

    if
        Value > Upper ->
            erlang:min(Bin + 1, Hist#hist.capacity - 1);
        Value =< Lower ->
            erlang:max(Bin - 1, 0);
        Value == Hist#hist.max ->
            Hist#hist.capacity-1;
        true ->
            Bin
    end.

            
quantile_itr(none, _Samples, _MaxSamples) ->
    max;
quantile_itr({Bin, Counter, Itr2}, Samples, MaxSamples) ->
    Samples2 = Samples + Counter,
    if
        Samples2 < MaxSamples ->
            %% Not done yet, move to next bin
            quantile_itr(gb_trees:next(Itr2), Samples2, MaxSamples);
        true ->
            %% We only need some of the samples in this bin; we make
            %% the assumption that values within the bin are uniformly
            %% distributed.
            Bin + ((MaxSamples - Samples) / Counter)
    end.


bin_count(Bin, Hist) ->
    case gb_trees:lookup(Bin, Hist#hist.bins) of
        {value, Count} ->
            Count;
         none ->
            0
    end.
            
%% ===================================================================
%% Unit Tests
%% ===================================================================

-ifdef(EUNIT).

simple_test() ->
    %% Pre-calculated tests
    [7,0] = counts(update_all([10,10,10,10,10,10,14], new(10,18,2))).

-ifdef(EQC).


qc_count_check(Min, Max, Bins, Xs) ->
    LCounts = counts(update_all(Xs, new(Min, Max, Bins))),
    RCounts = basho_stats_utils:r_run(Xs,
        ?FMT("hist(x, seq(~w,~w,length.out=~w), plot=FALSE)$counts",
                                         [Min, Max, Bins+1])),
    case LCounts == RCounts of
        true ->
            true;
        _ ->
            io:format("LCounts ~p, RCounts ~p~n", [LCounts, RCounts]),
            false
    end.


prop_count() ->
    ?FORALL({Min, Bins, Xlen}, {choose(0, 99), choose(2, 20), choose(2, 100)},
        ?LET(Max, choose(Min+1, 100),
            ?LET(Xs, vector(Xlen, choose(Min, Max)),
                ?WHENFAIL(
                    begin
                            io:format("Min ~p, Max ~p, Bins ~p, Xs ~w~n",
                                [Min, Max, Bins, Xs]),
                            Command = ?FMT("hist(x, seq(~w,~w,length.out=~w), plot=FALSE)$counts",
                                [Min, Max, Bins+1]),
                            InputStr = [integer_to_list(I) || I <- Xs],
                            io:format(?FMT("x <- c(~s)\n", [string:join(InputStr, ",")])),
                            io:format(?FMT("write(~s, ncolumns=1, file=stdout())\n", [Command]))
                    end,
                    qc_count_check(Min, Max, Bins, Xs))))).

qc_count_test() ->
    true = eqc:quickcheck(prop_count()).

qc_quantile_check(Q, Min, Max, Bins, Xs) ->
    Hist = new(Min, Max, Bins),
    LCounts = counts(update_all(Xs, Hist)),
    Lq = quantile(Q * 0.01, update_all(Xs, Hist)),
    [Rq] = basho_stats_utils:r_run(Xs, ?FMT("quantile(x, ~4.2f, type=4)", [Q * 0.01])),
    case abs(Lq - Rq) < 1 of
        true ->
            true;
        false ->
            ?debugMsg("----\n"),
            ?debugFmt("Q: ~p Min: ~p Max: ~p Bins: ~p\n", [Q, Min, Max, Bins]),
            ?debugFmt("Lq: ~p != Rq: ~p\n", [Lq, Rq]),
            ?debugFmt("Xs: ~w\n", [Xs]),
            false
    end.

prop_quantile() ->
    %% Loosey-goosey checking of the quantile estimation against R's more precise method.
    %%
    %% To ensure a minimal level of accuracy, we ensure that we have between 50-200 bins
    %% and between 100-500 data points.
    %%
    %% TODO: Need to nail down the exact error bounds
    %%
    %% XXX since we try to generate the quantile from the histogram, not the
    %% original data, our results and Rs don't always agree and this means the
    %% test will occasionally fail. There's not an easy way to fix this.
    ?FORALL({Min, Bins, Xlen, Q}, {choose(1, 99), choose(50, 200), choose(100, 500),
                                    choose(0,100)},
             ?LET(Max, choose(Min+1, 100),
                  ?LET(Xs, vector(Xlen, choose(Min, Max)),
                      ?WHENFAIL(
                          begin
                                  io:format("Min ~p, Max ~p, Bins ~p, Q ~p, Xs ~w~n",
                                      [Min, Max, Bins, Q, Xs]),
                                  Command = ?FMT("quantile(x, ~4.2f, type=4)", [Q * 0.01]),
                                  InputStr = [integer_to_list(I) || I <- Xs],
                                  io:format(?FMT("x <- c(~s)\n", [string:join(InputStr, ",")])),
                                  io:format(?FMT("write(~s, ncolumns=1, file=stdout())\n", [Command]))
                          end,
                          qc_quantile_check(Q, Min, Max, Bins, Xs))))).

qc_quantile_test() ->
    true = eqc:quickcheck(prop_quantile()).

-endif.
-endif.