%       program hd1dsp.m
%
%       Single precision version of hd1dmp13.m
%       The subintervals are determined by iterating the mapping Nitp times
%       rather than the full Nit times. We do so to keep the minimum
%       interval size from becoming too small and to insure the
%       necessary separation between subintervals.
%       Current version:  August 13, 2020
%
%       Computes Hausdorff dimension of the Cantor set
%       E_{a_1, a_2,...,a_R}
%       consisting of all points x in the interval (0,1) with
%       the continued fraction expansion of
%       x = [a_1, a_2, a_3, ...] =                  1
%                                        ----------------------
%                                          a_1 +     1
%                                                --------
%                                                a_2 +    1
%                                                     ---------
%                                                     a_3 + ...
%
%      where each a_n is one of the integers a_1, ... a_R.
%
%  Uses the routines thetabxitpsp.m, thetabxitsp.m, and eval1dPnitsp.m
%
%     We approximate by continuous piecewise polynomials of degree <= deg, 
%     deg >= 1, after first iterating the mappings theta_m.
%
%  deg = degree of piecewise polynomila
%  xlen = length of eigenvector
%  hlim = maximum allowable length of a subinterval
%  R = number of mappings b used in the computation
%  i.e., the number of a_i in the set E.
%  Nit = number of times the mappings theta_m are iterated
%  Nit = 0 means there is no iteration of the mappings
%  Nitp = number of times mappings theta_m are iterated
%  for the purpose of defining subintervals.
%  Here, Nitp=-1 means we use the basic interval [cin,din]
%  Nitp = 0 means we only use the original mappings
%  omega = (J(1), \ldots, J(Nit + 1)) \in B_{mu},  J(1) = b
%  theta_{omega} = (theta_{J(1)} \circ theta_{J(2)} \circ \ldots
%     \circ theta_{J((Nit)} \circ theta_{J(Nit +1)}
%  

format long
format longG  % Toolbox shows all digits in case of 'long' formats.
% Input data:
       tol = mp('10^(-16)')
%
% deg = the degree of the piecewise polynomial
      prompt = 'Input degree of piecewise polynomial ';  
      deg = input(prompt)

% hlim = maximum allowable length of a subinterval
      prompt = 'Input maximum allowable length of a subinterval ';  
      hlim = input(prompt)
%
      mu = 10; 
% mu = hmax/hmin; hmax = max interval size; hmin = min interval size

%   The value of R must be input for each run
%   To do E[1,3], input R = 2.
      prompt = 'Input R ';  
      R = input(prompt)
      blen = R

% Input the values of the vector b in increasing order 
% one at a time at the prompt. To do E[1,3], input 1 at the first
% prompt and 3 at the second prompt.
% Note that vector b is stored in bin in decreasing order
% so the largest value of b is in bin(blen) and the 
% smallest value of b is in bin(1).
        bin = zeros(blen,1);
        for j =1: blen
            prompt = 'Input values of vector b in increasing order one at a time after each prompt.'; 
            bin(blen+1-j) = input(prompt);
        end
        bin

% We take the initial interval [cin,din] to depend on bin(1) and b(blen).
		   quot = bin(blen)/bin(1);
		   cin =  - bin(blen)/2 + sqrt((bin(blen)/2)^2 + quot)
                   din =  cin/quot

% Input the number of times the theta_b are to be iterated
       prompt = 'Input Nit '; 
       Nit =  input(prompt)
% The mappings theta_m are iterated to produce a new operator
% whose Lipschitz constant has been reduced from the Lipschitz constant
% of the original operator.  This is required to satisfy the hypotheses
% of the theory.
% Compute the matrix of permutations of the vector b, 
% allowing Nit repetitions.
       Pbs = permn([bin], Nit+1);

% Compute the number of rows in this matrix, which gives the
% number of terms in the sum of the function whose eigenvalues we seek.
       Pblen = size(Pbs,1)

% Although we can define a mesh on the domain [cin,din], itself, 
% we can reduce the number of mesh points by computing
% a new domain contained in the interval [cin,din] by
% applying iterates of the mappings theta_m to the interval [cin,din].
% We first find an unsorted version, [cend(i), dend(i)]
      cend = zeros(Pblen,1);
      dend = zeros(Pblen,1);
      parfor i = 1:Pblen
     	   cend(i) = thetabxitsp(din,Pbs(i,:),Nit);
	   dend(i) = thetabxitsp(cin,Pbs(i,:),Nit);
	   if (cend(i) > dend(i))
             [cend(i), dend(i)] = deal(dend(i), cend(i)); % reverses arguments
           end
      end
% We next obtain a sorted (by left endpoint) version of these intervals 
                   [co,cindxo] = sort(cend);
                   [do,dindxo] = sort(dend);
%    Sort the elements of Pbs in the same order as the co's.
                Pb = Pbs(cindxo,:);

% Although it is possible to use the domain obtained by
% Nit iterations of the mappings, this can lead to very 
% small subinterval sizes. Instead, we determine a new domain 
% by iterating only Nitp number of times, while still using 
% the mappings obtained by Nit iterations to calculate the mapping L.
% Nitp is determined so that the length of the smallest interval
% will be >= hlim/mu.
% The interval with smallest length is theta_{omega}([cin,din])
% where omega consists of all entries equal to bin(1).

% In addition to the criteria that the minimum subinterval length
% be >= hlim/mu, we also require that the gap bewteen two subintervals,
% be sufficiently large.
% slength = length of shortest subinterval
slength =0;
mingapc = -1;
Nitp = Nit+1;
while ((slength < hlim/mu) || (mingapc < 0) && (Nitp > 0))
   Nitp = Nitp - 1;
   if Nitp == 0
     slength = din - cin;
   else
     omegaG = permn(bin(1), Nitp+1);
     slength = abs(thetabxitsp(cin,omegaG,Nitp) - thetabxitsp(din,omegaG,Nitp));
   end

if (Nitp >=0)
   Pbsp = permn([bin], Nitp+1);
% Compute the number of rows in this matrix.
   Pblenp = size(Pbsp,1);
else
   Pblenp = 1;
end
Pblenp

% The new domain will be the union of
% the intervals [c(i), d(i)], i = 1, ..., Pblenp.
% We first find an unsorted version, [cendp(i), dendp(i)]
% We assume that the intervals produced by this procedure are disjoint.
%
   cendp = zeros(Pblenp,1);
   dendp = zeros(Pblenp,1);
   if Nitp == -1
	   cendp(1) = cin;
           dendp(1) = din;
   else
      parfor i = 1:Pblenp
        cendp(i) = thetabxitsp(din,Pbsp(i,:),Nitp);
        dendp(i) = thetabxitmp(cin,Pbsp(i,:),Nitp);
        if (cendp(i) > dendp(i))
         [cendp(i), dendp(i)] = deal(dendp(i), cendp(i)); % reverses arguments
        end
      end
   end

% We next obtain a sorted (by left endpoint) version of these intervals 
		   [c,cindx] = sort(cendp);
		   [d,dindx] = sort(dendp);
% Sort the elements of Pbsp in the same order as the c's.
                   Pbp = Pbsp(cindx,:);

%      We approximate a function on each of the intervals [c_i,d_i]
%      by a piecewise polynomial of degree <= deg.
%      We next break up each of the intervals [c_i,d_i],
%      i = 1, ... Pblenp, into subintervals with endpoints [a_{i,j}, b_{i,j}].
%      To determine the number of such subintervals on each interval,
%      first calculate the length of the intervals, the maximum of these 
%      lengths, and then the total length of all subintervals.
                   dc = d-c;
                   dcmax = max(d-c);
		   totlen = sum(dc)
%                  ratio = totlen/Pblenp;
%      Since we want no subinterval to have length >= hlim
%      we further divide each of the intervals [c(i),d(i)] 
%      into Ns(i) subintervals, with Ns(i) >= 1.
		   Ns = ceil(dc/hlim);
                   Ns = max(Ns,1);
%
%      Then calculate the length, h(j), of each of the
%       Ns(j) subintervals of [c(j), d(j)].
		   h = gdivide(dc,Ns);
%                        
% In addition to the criteria that the minimum subinterval length
% be >= hlim/mu, we also require that the gap bewteen two subintervals,
% c(i+1) - d(i), i = 1,... Pblenp-1,
% be >= h(i)*(sin(pi/(2*deg+2)))^2
		   for i = 1:Pblenp-1
            	      gap(i) = c(i+1) - d(i);
   		      gapcomp(i) = gap(i)- h(i)*(sin(pi/(2*deg+2)))^2;
                   end
                   mingapc = min(gapcomp)
end
%      Compute total number of subintervals
                   Nssum = sum(Ns)
		   hmax = max(h)
         	   hexpr = hmax^(deg)

%      The subintervals will have endpoints [a_{i,j}, b_{i,j}] given by:
        for i = 1:Pblenp
               for j= 1:Ns(i)
		   a(i,j) = c(i) + (j-1)*h(i);
               end
               for j= 1: Ns(i)-1
                   b(i,j) = a(i,j+1);
               end
		   b(i,Ns(i)) = d(i);
         end
%
%     We compute the index number, aix(i,j), of the points a(i,j)
      parfor j= 1:Ns(1)
		   aix(1,j) = 1 + (j-1)*deg;
      end
     for i = 2:Pblenp
        for j= 1:Ns(i)
		   aix(i,j) = aix(i-1,Ns(i-1)) + 1 + j*deg;
        end
     end
%	aix
%    
%  On each of these subintervals, we interpolate a function
%  at the deg + 1 extended Chebyshev points.  These are the two end points
%  of the subinterval plus deg -1 interior points given by the 
%  following formulas.

% The extended Chebyshev points in [-1,1] are given by
%
    frac = pi/(2*deg + 2);
    parfor k =1:deg+1
  	chvhat(k) =  - cos(mp(2*k-1)*frac)/cos(frac);
    end
%
% The extended Chebyshev points in [a(i,j),b(i,j)] are given by	chv(i,j,k)  
% These are calculated below along with a single index of all mesh points.
% in increasing order, which we denote by x
    index =0;
    for i=1:Pblenp
      for j= 1:Ns(i)
           chv(i,j,1) = a(i,j);
           x(index + 1) = chv(i,j,1);
          for k=2:deg
       	   chv(i,j,k) = a(i,j) + h(i)*(1 + chvhat(k))/2;
		   x(index + k) = chv(i,j,k);
          end
		   index = index + deg;
                   chv(i,j,deg+1) = b(i,j);
      end

		   index = index +1;
		   x(index) = d(i);
    end
% Note that the total number of interpolation points xlen will be  = index
        xlen = index

%
% For each of the mesh points x and values of J in Pbs, compute the functions
% aphib = thetabxitp
%
%  Note that elements in Pb are the sorted version of elements in Pbs.
      aphib = zeros(xlen,Pblen);
      parfor k = 1:xlen
         for l = 1:Pblen 
           aphib(k,l) = thetabxitpsp(x(k), Pb(l,:), Nit);
         end
      end
% disp ('Computed aphib = thetabxitpsp.')
%
%   We next compute the Lagrange basis functions Lx corresponding 
%   to each of these degrees of freedom at the point thetabxkl.
%   The indices k,l in L(k,l,i) determine the point x = thetabx(k,l)
%   and the index i is the ith Lagrange basis function on the subinterval
%   containing the point x = thetabx(k,l), and whose Lagrange 
%   interpolation points are x(ix(k,l)), ..., x(ix(k,l) + deg).
%     ixkl = ix(k,l);
denom = ones(deg+1,1);
parfor  i = 1: deg+1
        for j = 1:deg +1
	   if (i~=j)
	     denom(i) = (chvhat(i) - chvhat(j))*denom(i);
           end
        end
end
% display('Computed denom')
numer = ones(xlen,Pblen);
%  To compute the numerator, we first determine
%  for each of the points thetabxkl (depending on x(k), Pb(l,:), and Nit),
%  the index number ix(k,l) of the left mesh point 
%  of the subinterval to which it belongs. There are Nssum subintervals 
%  with each one having deg-1 interior mesh points.
%  Note that we must decide how to choose the subinterval 
%  when the points lie on the boundary of two subintervals.
%  The index number refers to the sorted intervals.  
%
		   ix = zeros(xlen,Pblen);
		   Lx = ones(xlen,Pblen,deg+1);
% Compute ix and Lx
%

parfor k = 1:xlen
   for l = 1:Pblen
      thetabxkl = thetabxitsp(x(k),Pb(l,:),Nit);
%  We first determine which interval [c(jl), d(jl)] thetabxkl belongs to.
          Pbshort = Pb(l,1:Nitp+1);
	  [cout,jl,ib] = intersect(Pbp,Pbshort,'rows');
%  Even though theoretically, thetabxkl belongs to the interval
%  [c(jl), d(jl)], because of roundoff error, it
%  could happen that thetabxkl < c(jl). We fix this as follows.
          thetabxkl = max(thetabxkl,c(jl));
          jg = 1;
       if (Ns(jl) > 1)
          testv =  max(thetabxkl -a(jl,2: Ns(jl)),0);
          if (max(testv >0))
            jg = find(testv,1,'last')+ 1;
          end
       end
       ix(k,l) = aix(jl,jg);
       ixkl = ix(k,l);

%
% If the index is equal to xlen, then reset to deg less
      if (ixkl == xlen)
         ixkl = ixkl - deg;
      end

% Compute (2/h(jl))^r times augmented numerator of the Lagrange polynomial.
           fact = thetabxkl - x(ixkl: ixkl +deg);
           numer(k,l) = prod(2*fact/h(jl));
% Compute the extra factor in numer(k,l) depending on i 
           for i = 1:deg+1
% Combine the numerator with the denominator and the extra factor
               if (abs(fact(i)) ~= 0)
                 Lx(k,l,i) = numer(k,l)*h(jl)/(2*denom(i)*fact(i));
               end
           end
   end
end

% disp ('Computed Lagrange basis functions.')
%
disp ('Compute value of Hausdorff dimension s.')
%
% Compute the root of log(max evalue(s)) = 0 using the Matlab routine fzero.
%  
% tolerance must be changed with fewer digits
options = optimset('Display','iter', 'TolX', 10^(-15));
fun = @(s) log(eval1dPnitsp(s, ix, aphib, xlen, Pblen, Lx, deg));
disp ('Value of Hausdorff dimension s.')
% The initial interval must be changed depending on the problem
% The initial interval must be changed depending on the problem
      prompt = 'Input lendpt of interval containing Hausdorff dimension s';  
      lendpt = input(prompt)
      prompt = 'Input rendpt of interval containing Hausdorff dimension s';  
      rendpt = input(prompt)
% The general case
[s fvAL, exitflag output] = fzero(fun, [lendpt, rendpt],options)
% The case E[1,2]
% [s fvAL, exitflag output] = fzero(fun, [0.53, 0.54],options)
%
% The case E[1,3]
%[s fvAL, exitflag output] = fzero(fun, [0.45,0.46],options)
%
% The case E[2,3]
%[s fvAL, exitflag output] = fzero(fun, [0.33 , 0.34],options)
%
% The case E[10,11]
%[s fvAL, exitflag output] = fzero(fun, [0.14, 0.15],options) 
%
% The case E[100,1000]
%[s fvAL, exitflag output] = fzero(fun, [0.06, 0.07],options)
%
% The case E[2,4,6,8,10]
%[s fvAL, exitflag output] = fzero(fun, [0.51, 0.52],options)
%
% The case E[1, ..., 10]
%[s fvAL, exitflag output] =  fzero(fun, [0.92, 0.93],options)
%
mua = eval1dPnitsp(s, ix, aphib, xlen, Pblen, Lx, deg)
disp ('The eigenvalue of the matrix minus one.')
disp ('The value should be approximately 0.')
		  mua1 = mua - 1
% Compute hmax and hmin
%hmax = max(h)
%hmin = min(h)

%       end program hd1dsp.m

