function varargout = Arith07(xC)
% Arith07     Arithmetic encoder or decoder 
% Vectors of integers are arithmetic encoded, 
% these vectors are collected in a cell array, xC.
% If first argument is a cell array the function do encoding,
% else decoding is done.
%
% [y, Res] = Arith07(xC);                    % encoding
% y = Arith07(xC);                           % encoding
% xC = Arith07(y);                           % decoding
% ------------------------------------------------------------------% Arguments:
%  y        a column vector of non-negative integers (bytes) representing 
%           the code, 0 <= y(i) <= 255. 
%  Res      a matrix that sum up the results, size is (NumOfX+1)x4
%           one line for each of the input sequences, the columns are
%           Res(:,1) - number of elements in the sequence
%           Res(:,2) - unused (=0) 
%           Res(:,3) - bits needed to code the sequence
%           Res(:,4) - bit rate for the sequence, Res(:,3)/Res(:,1)
%           Then the last line is total (which include bits needed to store NumOfX)
%  xC       a cell array of column vectors of integers representing the
%           symbol sequences. (should not be to large integers)
%           If only one sequence is to be coded, we must make the cell array
%           like: xC=cell(2,1); xC{1}=x; % where x is the sequence
% ------------------------------------------------------------------
% Note: this routine is extremely slow since it is all Matlab code
% SOME NOTES ON THE FUNCTION  
% This function is almost like Arith06, but some important changes have
% been done. Arith06 is buildt almost like Huff06, but this close connection
% is removed in Arith07. This imply that to understand the way Arith06
% works you should read the dokumentation for Huff06 and especially the
% article on Recursive Huffman Coding. To understand how Arith07 works it is 
% only confusing to read about the recursive Huffman coder, Huff06.
%
% Arith07 code each of the input sequences in xC in sequence, the
% method used for each sequence depends on what kind (xType) the 
% sequence is. We have
%  xType  Explanation
%   0     Empty sequence (L=0)
%   1     Only one symbol (L=1) and x>1
%   9     Only one symbol (L=1) and x=0/1
%  11     Only one symbol (L=1) and x<0
%   2     all symbols are equal, L>1, x(1)=x(2)=...=x(L)>1
%   8     all symbols are equal, L>1, x(1)=x(2)=...=x(L)=0/1
%  10     all symbols are equal, L>1, x(1)=x(2)=...=x(L)<0
%   3     non-negative integers, 0<=x(l)<=1000
%   4     not to large integers, -1000<=x(l)<=1000
%   5     large non-negative integers possible, 0<=x(l)
%   6     also possible with large negative numbers
%   7     A binary sequence, x(l)=0/1. 
% Many functions are defined in this m-file:
%  PutBit, GetBit     read and write bit from/to bitstream
%  PutABit, GetABit   Arithmetic encoding of a bit, P{0}=P{1}=0.5
%  PutVLIC, GetVLIC   Variable Length Integer Code
%  PutS(x,S,C), GetS(S,C)  A symbol x, which is in S, is aritmetic coded
%                     according to the counts given by C.
%                     Ex: A binary variable with the givene probabilities
%                     P{0}=0.8 and P{1}=0.2, PutS(x,[0,1],[4,1]).
%  PutN(x,N), GetN(N) Arithmetic coding of a number x in range 0<=x<=N where
%                     we assume equal probability, P{0}=P{1}=...=P{N}=1/(N+1)

%----------------------------------------------------------------------
% Copyright (c) 1999-2001.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  karl.skretting@tn.his.no   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:
% Ver. 1.0  19.05.2001  KS: Function made (based on Arith06 and Huff06)
%----------------------------------------------------------------------

% these global variables are used to read from or write to the compressed sequence
global y Byte BitPos      
% and these are used by the subfunctions for arithmetic coding
global high low range ub hc lc sc K code

Mfile='Arith07';
K=24;  % number of bits to use in integers  (4 <= K <= 24)
Display=1;      % display progress and/or results

% check input and output arguments, and assign values to arguments
if (nargin < 1); 
   error([Mfile,': function must have input arguments, see help.']); 
end
if (nargout < 1); 
   error([Mfile,': function must have output arguments, see help.']); 
end

if (~iscell(xC))
   Encode=0;Decode=1;
   y=xC(:);            % first argument is y
   y=[y;0;0;0;0];      % add some zeros to always have bits available
else
   Encode=1;Decode=0;
   NumOfX = length(xC);
end

Byte=0;BitPos=1;         % ready to read/write into first position
low=0;high=2^K-1;ub=0;   % initialize the coder
code=0;
% we just select some probabilities which we belive will work quite well   
TypeS=[  3,  4,  5,  7,  6,  1,  9, 11,  2,  8, 10,  0];
TypeC=[100, 50, 50, 50, 20, 10, 10, 10,  4,  4,  2,  1];

if Encode
   Res=zeros(NumOfX,4);
   % initalize the global variables
   y=zeros(10,1);    % put some zeros into y initially
   % start encoding, first write VLIC to give number of sequences
   PutVLIC(NumOfX);
   % now encode each sequence continuously
   Ltot=0;
   for num=1:NumOfX
      x=xC{num};
      x=full(x(:));        % make sure x is a non-sparse column vector
      L=length(x);
      Ltot=Ltot+L;
      y=[y(1:Byte);zeros(50+2*L,1)];  % make more space available in y
      StartPos=Byte*8-BitPos+ub;      % used for counting bits
      % find what kind of sequenqe we have, save this in xType
      if (L==0)
         xType=0;
      elseif L==1
         if (x==0) | (x==1)       % and it is 0/1
            xType=9; 
         elseif (x>1)             % and it is not 0/1 (and positive)
            xType=1;
         else                     % and it is not 0/1 (and negative) 
            xType=11;
         end
      else
         % now find some more info about x 
         maxx=max(x);
         minx=min(x);
         rangex=maxx-minx+1;
         if (rangex==1)                    % only one symbol
            if (maxx==0) | (maxx==1)       % and it is 0/1
               xType=8; 
            elseif (maxx>1)                % and it is not 0/1 (and positive)
               xType=2;
            else                           % and it is not 0/1 (and negative) 
               xType=10;
            end
         elseif (minx == 0) & (maxx == 1)   % a binary sequence
            xType=7;
         elseif (minx >= 0) & (maxx <= 1000)
            xType=3;
         elseif (minx >= 0) 
            xType=5;
         elseif (minx >= -1000) & (maxx <= 1000)
            xType=4;
         else
            xType=6;
         end
      end        % if (L==0)
      if Display >= 2
         disp([Mfile,': sequence ',int2str(num),' has xType=',int2str(xType)]);
      end
      % 
      if sum(xType==[4,6])        % negative values are present          
         I=find(x);               % non-zero entries in x
         Sg=(sign(x(I))+1)/2;     % the signs will be needed later, 0/1
         x=abs(x);   
      end
      if sum(xType==[5,6])        % we take the logarithms of the values          
         I=find(x);               % non-zero entries in x
         xa=x;                    % additional bits
         x(I)=floor(log2(x(I)));
         xa(I)=xa(I)-2.^x(I);
         x(I)=x(I)+1;
      end
      %  now do the coding of this sequence
      PutS(xType,TypeS,TypeC);
      if (xType==1)       % one symbol and x=x(1)>1
         PutVLIC(x-2);
      elseif (xType==2)   % L>1 but only one symbol, x(1)=x(2)=...=x(L)>1
         PutVLIC(L-2);
         PutVLIC(x(1)-2);
      elseif sum(xType==[3,4,5,6])  % now 'normalized' sequences: 0 <= x(i) <= 1000
         PutVLIC(L-2);     
         M=max(x);
         PutN(M,1000);         % some bits for M
         % initialize model
         T=[ones(M+2,1);0];
         Tu=flipud((-1:(M+1))');   % (-1) since ESC never is used in Tu context
         % and code the symbols in the sequence x
         for l=1:L
            sc=T(1);
            m=x(l); 
            hc=T(m+1);lc=T(m+2);
            if hc==lc      % unused symbol, code ESC symbol first
               hc=T(M+2);lc=T(M+3);
               EncodeSymbol;      % code escape with T table
               sc=Tu(1);hc=Tu(m+1);lc=Tu(m+2);  % symbol with Tu table
               Tu(1:(m+1))=Tu(1:(m+1))-1;       % update Tu table
            end
            EncodeSymbol;  % code actual symbol with T table (or Tu table)
            % update T table, MUST be identical in Encode and Decode
            % this avoid very large values  in T even if sequence is very long
            T(1:(m+1))=T(1:(m+1))+1; 
            if (rem(l,5000)==0)  
               dT=T(1:(M+2))-T(2:(M+3));
               dT=floor(dT*7/8+1/8);
               for m=(M+2):(-1):1; T(m)=T(m+1)+dT(m); end;
            end
         end          % for l=1:L
         % this end the "elseif sum(xType==[3,4,5,6])"-clause
      elseif (xType==7)   % L>1 and   0 <= x(i) <= 1
         PutVLIC(L-2);
         EncodeBin(x,L);    % code this sequence a special way
      elseif (xType==8)   % L>1 and   0 <= x(1)=x(2)=...=x(L) <= 1
         PutVLIC(L-2);
         PutABit(x(1));
      elseif (xType==9)   % L=1 and   0 <= x(1) <= 1
         PutABit(x(1));
      elseif (xType==10)       % L>1 and  x(1)=x(2)=...=x(L) <= -1
         PutVLIC(L-2);
         PutVLIC(-1-x(1));
      elseif (xType==11)       % L=1 and   x(1) <= -1
         PutVLIC(-1-x);
      end         % if (xType==1)   
      % additional information should be coded as well
      if 0         % first the way it is not done any more
         if sum(xType==[4,6])           % sign must be stored          
            for i=1:length(Sg); PutABit(Sg(i)); end;   
         end
         if sum(xType==[5,6])        % additional bits must be stored          
            for i=1:L
               for ii=(x(i)-1):(-1):1
                  PutABit(bitget(xa(i),ii));
               end
            end
         end
      else         % this is how we do it
         if sum(xType==[4,6])           % sign must be stored          
            EncodeBin(Sg,length(I));    % since length(I)=length(Sg)
         end
         if sum(xType==[5,6])        % additional bits must be stored          
            b=zeros(sum(x)-length(I),1);  % number of additional bits
            bi=0;
            for i=1:L
               for ii=(x(i)-1):(-1):1
                  bi=bi+1;
                  b(bi)=bitget(xa(i),ii);
               end
            end
            if (bi~=(sum(x)-length(I)))
               error([Mfile,': logical error, bi~=(sum(x)-length(I)).']); 
            end
            EncodeBin(b,bi);    % since bi=(sum(x)-length(I))
         end
      end
      %
      EndPos=Byte*8-BitPos+ub;    % used for counting bits
      bits=EndPos-StartPos;
      Res(num,1)=L;
      Res(num,2)=0;
      Res(num,3)=bits;
      if L>0; Res(num,4)=bits/L; else Res(num,4)=bits; end;
      if Display
         disp([Mfile,': Sequence ',int2str(num),' of ',int2str(L),' symbols ',...
               'encoded using ',int2str(bits),' bits.']);
      end
   end         % for num=1:NumOfX
   % flush the arithmetic coder
   PutBit(bitget(low,K-1));     
   ub=ub+1;
   while ub>0
      PutBit(~bitget(low,K-1));
      ub=ub-1;
   end
   % flush is finished
   y=y(1:Byte);   
   varargout(1) = {y};
   if (nargout >= 2) 
      % now calculate results for the total
      Res(NumOfX+1,1)=Ltot;
      Res(NumOfX+1,2)=0;
      Res(NumOfX+1,3)=Byte*8;
      if (Ltot>0); Res(NumOfX+1,4)=Byte*8/Ltot; else Res(NumOfX+1,4)=Byte*8; end;
      varargout(2) = {Res}; 
   end
end         % if Encode

if Decode
   for k=1:K
      code=code*2;
      code=code+GetBit;   % read bits into code
   end
   NumOfX=GetVLIC;   % first read number of sequences
   xC=cell(NumOfX,1);
   for num=1:NumOfX
      % find what kind of sequenqe we have, xType, stored first in sequence
      xType=GetS(TypeS,TypeC);
      % now decode the different kind of sequences, each the way it was stored 
      if (xType==0)       % empty sequence, no more symbols coded
         x=[];
      elseif (xType==1)       % one symbol and x=x(1)>1
         x=GetVLIC+2;
      elseif (xType==2)   % L>1 but only one symbol, x(1)=x(2)=...=x(L)>1
         L=GetVLIC+2;
         x=ones(L,1)*(GetVLIC+2);
      elseif sum(xType==[3,4,5,6])  % now 'normalized' sequences: 0 <= x(i) <= 1000
         L=GetVLIC+2;
         x=zeros(L,1);
         M=GetN(1000);      % M is max(x)
         % initialize model
         T=[ones(M+2,1);0];
         Tu=flipud((-1:(M+1))');   % (-1) since ESC never is used in Tu context
         % and decode the symbols in the sequence x
         for l=1:L
            sc=T(1);
            range=high-low+1;
            count=floor(( (code-low+1)*sc-1 )/range);
            m=2; while (T(m)>count); m=m+1; end; 
            hc=T(m-1);lc=T(m);m=m-2;
            RemoveSymbol;
            if (m>M)     % decoded ESC symbol, find symbol from Tu table
               sc=Tu(1);range=high-low+1;
               count=floor(( (code-low+1)*sc-1 )/range);
               m=2; while (Tu(m)>count); m=m+1; end; 
               hc=Tu(m-1);lc=Tu(m);m=m-2;
               RemoveSymbol;
               Tu(1:(m+1))=Tu(1:(m+1))-1;   % update Tu table
            end
            x(l)=m;
            % update T table, MUST be identical in Encode and Decode
            % this avoid very large values  in T even if sequence is very long
            T(1:(m+1))=T(1:(m+1))+1; 
            if (rem(l,5000)==0)  
               dT=T(1:(M+2))-T(2:(M+3));
               dT=floor(dT*7/8+1/8);
               for m=(M+2):(-1):1; T(m)=T(m+1)+dT(m); end;
            end
         end          % for l=1:L
         % this end the "elseif sum(xType==[3,4,5,6])"-clause
      elseif (xType==7)   % L>1 and   0 <= x(i) <= 1
         L=GetVLIC+2;
         x=DecodeBin(L);    % decode this sequence a special way
      elseif (xType==8)   % L>1 and   0 <= x(1)=x(2)=...=x(L) <= 1
         L=GetVLIC+2;
         x=ones(L,1)*GetABit;
      elseif (xType==9)   % L=1 and   0 <= x(1) <= 1
         x=GetABit;
      elseif (xType==10)       % L>1 and  x(1)=x(2)=...=x(L) <= -1
         L=GetVLIC+2;
         x=ones(L,1)*(-1-GetVLIC);
      elseif (xType==11)       % L=1 and   x(1) <= -1
         x=(-1-GetVLIC);
      end         % if (xType==0)   
      % additional information should be decoded as well
      L=length(x);
      I=find(x);
      if 0
         if sum(xType==[4,6])           % sign must be retrieved          
            Sg=zeros(size(I));
            for i=1:length(I); Sg(i)=GetABit; end;   % and the signs   (0/1)
            Sg=Sg*2-1;                               % (-1/1)
         end
         if sum(xType==[5,6])        % additional bits must be retrieved
            xa=zeros(L,1);
            for i=1:L
               for ii=2:x(i)
                  xa(i)=2*xa(i)+GetABit;
               end
            end
            x(I)=2.^(x(I)-1);
            x=x+xa;
         end
      else
         if sum(xType==[4,6])           % sign must be retrieved          
            Sg=DecodeBin(length(I));    % since length(I)=length(Sg)
            Sg=Sg*2-1;                  % (-1/1)
         end
         if sum(xType==[5,6])        % additional bits must be retrieved
            bi=sum(x)-length(I);     % number of additional bits
            b=DecodeBin(bi);  
            bi=0;
            xa=zeros(L,1);
            for i=1:L
               for ii=2:x(i)
                  bi=bi+1;
                  xa(i)=2*xa(i)+b(bi);
               end
            end
            x(I)=2.^(x(I)-1);
            x=x+xa;
         end
      end
      if sum(xType==[4,6])           % sign must be used          
         x(I)=x(I).*Sg;
      end
      % now x is the retrieved sequence      
      xC{num}=x;
   end         % for num=1:NumOfX
   varargout(1) = {xC}; 
end

return     % end of main function, Arith07
%----------------------------------------------------------------------
%----------------------------------------------------------------------

% --- The functions for binary sequences: EncodeBin and DecodeBin -------
% These function may call themselves recursively
function EncodeBin(x,L)
global y Byte BitPos
global high low range ub hc lc sc K code

Display=0;
x=x(:);
if (length(x)~=L); error('EncodeBin: length(x) not equal L.'); end;

% first we try some different coding methods to find out which one
% that might do best. Many more methods could have been tried, for
% example methods to check if x is a 'byte' sequence, or if the bits
% are grouped in other ways. The calling application is the best place
% to detect such dependencies, since it will (might) know the process and
% possible also its statistical (or deterministic) properties.
% Here we just check some few coding methods: direct, split, diff, diff+split
% The main variables used for the different methods are:
%  direct: x, I, J, L11, b0
%  split: x is split into x1 and x2, L1, L2, b1 is bits needed
%  diff: x3 is generated from x, I3, J3 L31, b2
%  diff+split: x3 is split into x4 and x5, L4, L5, b3
%