function [ D, ind ] = GDictFast( name, pars )
% GDictFast generates dictionaries by translating prototype functions 
%
% Usage:    D = GDictFast( name, pars );
%
% Inputs:
%   name    type of dictionary (string format)
%   pars    parameters (cell format)
%           pars = {L,sp,j,a,b}
% Description of the parameters
%   L       number of points partitioning the interval sp=[sp(1) sp(2)]
%   sp      array sp=[c d] stands for interval [c,d]
%   j       array containing all the scales to be considered note that larger j implies 
%           finer scale
%   a       dilation factor (positive integer)
%   b       translation factor (b=1 for basis, keep 1/b being integer) for 'chui_cubic' 
%           and 'chui_cubicb' use b=1 only
%
% Outputs:
%   D       desired dictionary
%   ind     array of indices separating different scales in D 
%
% Comments:
% The dictionaries are constructed in the following way: a prototype function for every
% scale is computed for values x=linspace(sp(1),sp(2),L) (x is a dicretization of the 
% interval sp containing L nodes). Then this prototype function is shifted on x by an 
% appropriate number of nodes. For this end  (L-1)*(b*a^(-jmax))/(sp(2)-sp(1)) must be
% an integer where jmax stands for the finest scale parameter. In the case of 'chui_cubic'
% or 'chui_cubicb' dictionaries only the inner functions are constructed by the method 
% above. The boundary function are calculated using given analytical expressions.
%
% Only scales with at least one inner function are considered.
%
% References:
%   M. Andrle and L. Rebollo-Neira, "Spline wavelet dictionaries for non-linear signal 
%   approximation", preprint, 2005.
%   M. Andrle and L. Rebollo-Neira, "Cardinal B-spline dictionaries on a compact interval", 
%   Applied and Computational Harmonic Analysis, Vol(18,3) 336-346 (2005).

% See http://www.ncrg.aston.ac.uk/Projects/BiOrthog/ for more details
tic;
funname='GDictFast';
if nargin~=2; error('Need 2 input arguments');end
if ~ischar(name); error('NAME must be a string');end
name = lower(name);

L=pars{1}; c=pars{2}(1); d=pars{2}(2);
j=pars{3};b=pars{5};a=pars{4};j=sort(j);
if b<=0; error('I expect b>0');end;
if 1/b~=round(1/b)
    fprintf('\nIt is not good idea to use %s for b=%f',funname,b);
    error('\nIt is exact only if 1/b is an integer');
end
%x is a dicretization of the interval sp=[sp(1) sp(2)]
x=linspace(c,d,L);
%int is the length of the interval sp
int=d-c; 

if name(end)=='b'
    supp='supp=order;';
    namef=name(1:end-1);
else
    supp='supp=2*order-1;';
    namef=name;
end

order=str2num(namef);
if isempty(order)
    if strcmp(namef,'chui_cubic')
        order=4;
        fprintf('This family is implemented only on [0,K] for b=1, K integer\n');
    else
        error('Unknown name')
    end
else
    namef=order;
end
eval(supp); %evaluate the length of support 

zz=ceil(log(supp/int)/log(a));
%eliminate too coarse scales
%only scales with at least one inner function are considered
j=j(j>=zz); 
Lm=length(j);
if isempty(j)
    error('Please take j>%d, no inner function found',zz);
end
jmax=max(j);
%for the accuracy the parameters must satisfy the following conditions
pom=(L-1)*(b*a^(-jmax))/int;
if pom~=round(pom)
    fprintf('This method cannot be accurate with L=%d.\n',L);
    if int/(b*a^(-jmax))~=int/(b*a^(-jmax))
        fprintf('No integer L satisfying L=%g*J/(b*(a^(-%g)))+1 found\n',int,jmax);
        error('Try to change interval sp, or parameters a and b')
    else    
    fprintf('Suggested values of L are:\n');J=1:12;
    fprintf('%d, ',int*J/(b*(a^(-jmax)))+1);fprintf('\b\b.\n');
    error('Please choose L as L=%g*J/(b*(a^(-%g)))+1 for J being a positive integer',int,jmax);
    end
end
%end of the accuracy test

fprintf('\nInterval [%g,%g] divided into %d points, step=%g.\n',c,d,L,x(2)-x(1));
fprintf('\nDictionary ''%s'' for a=%g, b=%g, and j=[',name,a,b);
fprintf('%d ',j);fprintf('\b] is being constructed.\n');
fprintf('The smallest distance between two consecutive functions is %g.\n',b*a^(-jmax));

ind=zeros(Lm+1,1);
if name(end)=='b'
    %only B-splines, no wavelets
    for i=1:Lm
        ind(i+1)=ind(i)+NumFun(namef,'scal.f.',c,d,j(i),a,b);
    end
    D=zeros(L,ind(Lm+1));
    for mi=1:Lm,
        D(:,ind(mi)+1:ind(mi+1))=ScalLevel(namef,x,j(mi),a,b);
        fprintf('%2dth scale level %4d scaling f., indices %4d :%4d \n',j(mi),ind(mi+1)-ind(mi),ind(mi)+1,ind(mi+1));
    end
    ind(1)=[];
else
    %multiresolution-like dictionaries
    ind(1)=NumFun(namef,'scal.f.',c,d,j(1),a,b);
    for i=1:Lm
        ind(i+1)=ind(i)+NumFun(namef,'wavelet',c,d,j(i),a,b);
    end
    D=zeros(L,ind(Lm+1));
    D(:,1:ind(1))=ScalLevel(namef,x,j(1),a,b);
    fprintf('\n%2dth scale level %4d scaling f., indices %4d :%4d \n',j(1),ind(1),1,ind(1));
    for mi=1:Lm
        D(:,ind(mi)+1:ind(mi+1))=WavLevel(namef,x,j(mi),a,b);
        fprintf('%2dth scale level %4d wavelets,   indices %4d :%4d \n',j(mi),ind(mi+1)-ind(mi),ind(mi)+1,ind(mi+1));
    end
end
tt=toc;
fprintf('\n%s generated %d functions in %g seconds\n',funname,size(D,2),tt);
fprintf('   (%g functions per second)\n\n',size(D,2)/tt);


%Copyright (C) 2006 Miroslav ANDRLE and Laura REBOLLO-NEIRA
%
%This program is free software; you can redistribute it and/or modify it under the terms 
%of the GNU General Public License as published by the Free Software Foundation; either 
%version 2 of the License, or (at your option) any later version.
%
%This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%See the GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License along with this program;
%if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
%Boston, MA  02110-1301, USA.
