% This function provides a wrapper for the computeHGMex.  Precomputed data
% are computed and cached automatically.  The user can use text commands to
% query and control the caching process.
% 
%   If first argument is a string, it is treated as a command:
% 
%       computeHG('getProperties') - displays the current value of maxMem
%           and the amount of memory in use
%       computeHG('getCurMemInUse') - displays the amount of memory
%           currently in use
%       computeHG('getMaxMem') - displays current value of maxMem
%       computeHG('getData') - returns cached data tables
%       computeHG('setMaxMem', argument) - sets maxMem to argument
%       computeHG('clearData') - clears cache of data tables
%     
%   Otherwise, the function evalues a hypergeometric function:
% 
%       computeHG(mode, N, a, b, x, y) - computes the hypergeometric
%           function of one or two matrix arguments (computes a truncated
%           sum over partitions of size up to N).  Omit y to compute the 
%           hypergeometric function of single argument x.  mode indicates 
%           how much precomputation to do to speed up this and future calls 
%           to this function (acceptable values: 0, 1, or 2).
% 
% Written by Cy Chan, July 2007
% 
function [result, brokenDown] = computeHG(varargin)

    sizeOfInt = 4;

    persistent levelIndexTable;
    persistent dataArray;
    persistent maxMem;
    
    if (nargin == 0)
        fprintf('Usage: computeHG(mode, N, a, b, x, y)\n');
        fprintf('  Omit y for HG of single argument x.\n');
        fprintf('Usage: computeHG(command [, parameter])\n');
        fprintf('  Commands: ''getProperties'', ''getCurMemInUse'', ''getMaxMem'', ''getBackRefs'', ''setMaxMem'', ''clearBackRefs''.\n');
        return;
    end
        
    if (~isequal(varargin{1}, 'setMaxMem') && isempty(maxMem))
        disp('Default value of 500MB used for maxMem.  Set using command ''setMaxMem''');
        maxMem = 524288000;
    end

    % PROCESS TEXT COMMANDS

    if (ischar(varargin{1}))

        switch (varargin{1})
            case 'getProperties'
                temp1 = whos('dataArray');
                temp2 = whos('levelIndexTable');
                result.curMemInUse = temp1.bytes + temp2.bytes;
                result.maxMem = maxMem;
            case 'getCurMemInUse'
                temp1 = whos('dataArray');
                temp2 = whos('levelIndexTable');
                result = temp1.bytes + temp2.bytes;
            case 'getMaxMem'
                result = maxMem;
            case 'getData'
                result.levelIndexTable = levelIndexTable;
                result.dataArray = dataArray;
            case 'setMaxMem'
                if (nargin ~= 2 || ~isnumeric(varargin{2}) || ~isreal(varargin{2}) || length(varargin{2}) ~= 1)
                    error('Two arguments required for ''setMaxMem'' command.  Second argument must be a real scalar numeric value');
                end
                temp1 = whos('dataArray');
                temp2 = whos('levelIndexTable');
                if (varargin{2} < temp1.bytes + temp2.bytes)
                    error('Cannot set maxMem lower than current memory in use.  Use ''clearData'' to clear memory.');
                end
                maxMem = varargin{2};
                result = maxMem;
            case 'clearData'
                dataArray = [];
                levelIndexTable = [];
                result = dataArray;
            otherwise
                error('Invalid command.  Choose from: ''getProperties'', ''getCurMemInUse'', ''getMaxMem'', ''getData'', ''setMaxMem'', ''clearData''');
        end
        
        return;

    end
    
    % GET DATA INPUTS
    
    if (nargin ~= 5 && nargin ~= 6)
        fprintf('Usage: computeHG(mode, N, a, b, x, y)\n');
        fprintf('  Omit y for HG of single argument x.\n');
        fprintf('Usage: computeHG(command [, parameter])\n');
        fprintf('  Commands: ''getProperties'', ''getCurMemInUse'', ''getMaxMem'', ''getBackRefs'', ''setMaxMem'', ''clearBackRefs''.\n');
        return;
    end

    mode = varargin{1};
    N = varargin{2};
    a = varargin{3};
    b = varargin{4};
    x = varargin{5};
    if (nargin == 6)
        y = varargin{6};
    else
        y = NaN;
    end

    n = length(x);
    
    if (mode ~= 0 && mode ~= 1 && mode ~= 2)
        error('mode must be set to 0, 1, or 2');
    end

    % MODE = 0: USE LEVEL INDEX TABLE
    
    if (mode == 0)

        % increase the size of persistent levelIndexTable if necessary
        oldN = size(levelIndexTable, 1);
        oldn = size(levelIndexTable, 2);
        if (oldN < N || oldn < n)
            newN = max(N, size(levelIndexTable, 1));
            newn = max(n, size(levelIndexTable, 2));
            fprintf('Increasing size of levelIndexTable from (%d, %d, %d) to (%d, %d, %d) ...\n', oldN, oldn, oldN, newN, newn, newN);
            levelIndexTable = buildLevelIndexTableMex(newN, newn);
        end

        % run computation and return
        if (nargout == 2)
            [result, brokenDown] = computeHGMex(x, y, a, b, levelIndexTable, N);
        else
            result = computeHGMex(x, y, a, b, levelIndexTable, N);
        end
        
        return;

    end
    
    % MODE = 1 (or 2): USE BACK REFERENCE TABLE (and precomputed
    % coefficient data if MODE = 2)
    
    % try to find the correct backRefs structure in the persistent data array
    foundFlag = false;
    for (m = 1:length(dataArray))
        if (dataArray{m}.N == N)
            foundFlag = true;
            break;
        end
    end

    % check to see if we need to compute a new back reference table
    if (~foundFlag || length(dataArray{m}.outputLengths) < n)

        newSize = computeMemoryRequired(mode, N, n);
        if (newSize > maxMem)
            fprintf('maxMem (%d bytes) too low to precompute data (need %d bytes).\n', maxMem, newSize);
            error('Set higher maxMem using ''setMaxMem'' command.');
        end

        % clear old data for N out of memory
        if (foundFlag)
            fprintf('Evicting data for N = %d, n = %d ...\n', N, length(dataArray{m}.outputLengths));
            dataArray{m} = [];
        end

        % if not enough remaining memory, flush back reference cache.
        % TODO: Do something more intelligent than wiping the cache 
        % completely.  Need an eviction strategy.
        temp1 = whos('dataArray');
        temp2 = whos('levelIndexTable');
        curMemInUse = temp1.bytes + temp2.bytes;
        if (curMemInUse + newSize > maxMem)
            fprintf('Flushing all data to make room for new entry ...\n');
            dataArray = [];
            m = 1;
        end

        if (~foundFlag)
            m = length(dataArray) + 1;
        end
        tic; dataArray{m} = buildBackReferenceTableMex(N, n);
        fprintf('Time to compute back reference table: %g seconds\n', toc);

    end

    % check to see if we need to precompute coefficient data
    if (mode == 2)

        % if one matrix argument, figure out if we need to recompute
        % coeffData1 because n is too large for the current entry
        % 
        % if two matrix arguments, figure out if coeffData2 (a cell array)
        % contains an entry for n, and if so, set p to the correct index.
        % if not, set p to the index of a new entry.
        foundFlag = false;
        if (nargin == 5)
            if (isfield(dataArray{m}, 'coeffData1'))
                if (dataArray{m}.coeffData1.n >= n)
                    foundFlag = true;
                else
                    fprintf('Evicting coefficient data for N = %d, n = %d ...\n', N, dataArray{m}.coeffData1.n);
                    dataArray{m} = rmfield(dataArray{m}, 'coeffData1');
                end
            end
        else
            if (isfield(dataArray{m}, 'coeffData2'))
                for (p = 1:length(dataArray{m}.coeffData2))
                    if (dataArray{m}.coeffData2{p}.n == n)
                        foundFlag = true;
                        break;
                    end
                end
                if (foundFlag == false)
                    p = length(dataArray{m}.coeffData2) + 1;
                end
            else
                p = 1;
            end
        end
    
        % precompute coefficient data if needed
        if (~foundFlag)

            [totalSize, trash, newSize] = computeMemoryRequired(mode, N, n);
            if (totalSize > maxMem)
                fprintf('maxMem (%d bytes) too low to precompute data (need %d bytes).\n', maxMem, totalSize);
                error('Set higher maxMem using ''setMaxMem'' command.');
            end

            % flush data entries for other values of N if needed
            dataArray{m}.partitionSizes = []; % we're going to recompute this anyway
            temp1 = whos('dataArray');
            temp2 = whos('levelIndexTable');
            curMemInUse = temp1.bytes + temp2.bytes;
            if (curMemInUse + newSize > maxMem)
                fprintf('Flushing other data to make room for new entry ...\n');
                dataArray = dataArray(m);
                m = 1;
            end

            % flush coeffData2 entries for other values of n if needed
            temp1 = whos('dataArray');
            curMemInUse = temp1.bytes + temp2.bytes;
            if (curMemInUse + newSize > maxMem)
                fprintf('Flushing coeffData2 entries to make room for new entry ...\n');
                dataArray{m} = rmfield(dataArray{m}, 'coeffData2');
                p = 1;
            end
            
            if (nargin == 5)
                tic; [dataArray{m}.coeffData1, dataArray{m}.partitionSizes] = precomputeCoeffDataMex(N, n, dataArray{m}.outputLengths(n), 1);
            else
                tic; [dataArray{m}.coeffData2{p}, dataArray{m}.partitionSizes] = precomputeCoeffDataMex(N, n, dataArray{m}.outputLengths(n), 2);
            end
            fprintf('Time to precompute coefficient data: %g seconds\n', toc);
            
        end
        
    end
    
    % make the appropriate call to computeHGMex()
    if (nargout == 2)
        [result, brokenDown] = computeHGMex(x, y, a, b, dataArray{m});
    else
        result = computeHGMex(x, y, a, b, dataArray{m});
    end

return

function [totalSize, baseSize, coeffDataSize] = computeMemoryRequired(mode, N, n)

    sizeOfInt = 4;
    sizeOfDouble = 8;
    overheadPerField = 124;

    if (mode ~= 1 && mode ~= 2)
        error('Invalid mode selection');
    end
    
    [trash, outputLengths] = buildLevelIndexTableMex(N, n);
    
    baseSize = (1 ...                                       % N
              + n ...                                       % outputLengths
              + outputLengths(n - 1) * (n - 1) ...          % table
              + outputLengths(n) - outputLengths(n - 1) ... % extraReferences
              + outputLengths(n) ...                        % array
              + N * n) * sizeOfInt;                         % lastElementIndexTable
	baseSize = baseSize + 6 * overheadPerField;

    if (mode == 1)
        coeffDataSize = 0;
    else
        coeffDataSize = 1 * sizeOfInt ...                    % numMatrixArgs
                      + outputLengths(n) * sizeOfInt ...     % add
                      + outputLengths(n) * sizeOfDouble ...  % mult 
                      + outputLengths(n) * sizeOfInt;        % partitionSizes
        coeffDataSize = coeffDataSize + 5 * overheadPerField;
    end
    
    totalSize = baseSize + coeffDataSize;
    
return
