###########Copyright and Notice==============
#Copyright since 2001, Bin Han, University of Alberta,
#All Rights Reserved 

#The use of all the following software is permitted for
#non-commercial, educational, and research use only. The software
#and/or related materials are provided "as-is" without warranty of any
#kind including any warranties of performance or merchantability or
#fitness for a particular use or purpose or for any purpose
#whatsoever, for the licensed product, however used. In no event shall
#University of Alberta and/or Bin Han be liable for any damages and/or
#costs, including but not limited to incidental or consequential
#damages of any kind, including economic damage or injury to property
#and lost profits, regardless of whether University of Alberta shall
#be advised, have reason to know, or in fact shall know of the
#possibility. User bears all risk relating to quality and performance
#of the software and/or related materials. 

#Any use other than non-commercial, educational, or research, or any
#redistribution in original or modified form requires prior written
#authorization from the copyright holder.
 
#Report errors, mistakes, and bugs to Bin Han at bhan@ualberta.ca.
#Send comments/suggestions to Bin Han at bhan@ualberta.ca

##===================================================================
## MAPLE 6.0 routines for CBC algorithms on constructing 1-dimensional
## biorthogonal multiwavelets with a general dilation factor.

#This file consists of a set of useful MAPLE subroutines to handle 1D
#matrix masks with a general dilation factor.
###  Representation format of matrix mask:
#      poly(x)=sum(Mask[i]*x^i; i \in Z):
#So poly is a matrix of polynomials with the variables (x) reserved

#Usage:
#  the program uses Matrix and Vector from LinearAlgebra.
#      Do not use matrix and vector from linalg!!!

##D1PrintMask:=proc(poly,choice) 
##  Print a 1D matrix mask using various formats.
##  Input: the matrix of polynomials of the matrix mask.
##  Output: Display the matrix mask on screen in various formats.
##        and output the matrix mask in array format
##  choice<-1, Do not display, output matrix mask in numerical array
##           format (6 decimal places)
##  choice=-1, Do not display, output matrix mask in exact value array
##  choice=0,  Display the mask and ouput the exact value form,
##  otherwise, Display and output numerical form (6 Digits)

##D1PrintElement:=proc(poly,xp,choice) 
##  Print the element of a 1D matrix mask at poistion [xp].
##  Input: the matrix of polynomials of the matrix mask, xp,yp
##  Output: Display the element of the matrix mask on screen.
##  choice>=0, exact format, otherwise, numerical to 6 digits

##D1BiorthEqs:=proc(primal_poly,dual_poly,dilation,orth,choice) 
##  Find the equations of the (bi)orthogonal relation:
##      Primalpoly*\bar{Dualpoly} is an interpolatory mask w.r.t
##      Dilation*Z.
##  Input: primal_poly, dual_poly, dilation, orth, choice
##  orth=0 means dual filter: the primal-dual relation
##  orth=1 means orthogonal filter: lowpass-highpass relation
##  choice=0, equations in set for exact dual relation
##  otherwise, equations in set for dual relation without
##  zero-position 

##D1SumRuleEqs:=proc(poly,dilation,sr,knownvalue,choice)    
##  Find the equations for sum rules of order sr w.r.t.
##     Dilation*Z.
##  Input: poly,dilation,sr,choice,knownvalue
##  Output: result -> the equations for sum rules of order sr
##  choice:=0, then use standard definition of sum rules for
##      scalar mask. That is, sum(Mask[i+Mj](i+Mj)^\mu, j\in Z^s)
##      is independent of the integer i.
##  choice:=1, then each coset equals the knownvalue. That is,
##      sum(Mask[i+Mj](M^(-1)i+j)^\mu, j)=knownvalue[\mu].
##  choice=2, the y vector values are known for matrix mask.
##      that is, yvect[\mu]=knownvalue[\mu]
##  otherwise, we output the equations of sum rules with
##      parameterized y vectors (in yvect)
##  Do not assume that the sum of the mask is = |det(dilation)|
##  The definition of sum rules is defined in Eqs. (1.6) in
##    Bin Han, Approximation properties and construction of
##    Hermite Interpolants and Biorthogonal Multiwavelets,
##    Journal of Approximation Theory, Vol. 110, (2001) 18-53

##D1CBCAlgorithm:=proc(poly,dilation,sr,choice)
##  Given a matrix mask, find the \tilde y vector values by
##  the CBC algorithm
##  Input: poly,dilation,sr,choice
##  Output: result -> the \tilde y vectors in the sum rules
##        in Theorem 3.1 in the following paper
##        Bin Han, Approximation properties and construction of
##        Hermite Interpolants and Biorthogonal Multiwavelets,
##        Journal of Approximation Theory, Vol. 110, (2001) 18-53
##  Choice=0, output the moment for scalar mask:
##        h^a(\mu):=sum( \tilde a(i+Mj)(M^(-1)i+j)^\mu/\mu!, j)
##  Choice=otherwise, output the \tilde y vectors in Theorem 3.1

##D1MaskOfSym:=proc(sz,dilation,hermiteorder,interpolating,cc) 
##  Parametrize the mask using prescribed symmetry group {-1}.
##  Input: sz,dilation,hermiteorder,interpolating
##  Output: a matrix of polynomial for the parameterized matrix mask
##  When interpolating=0, it is not an interpolanting mask
##  w.r.t. dilation*Z, otherwise it is an interpolating matrix mask.
##cc is the name of the parameters: cc1, cc2, cc3,...

##D1YVectForHermiteMask:=proc(hmord,sr)
##  Given hermiteorder and sum rule order, generate the y vector for
##  refinable Hermite mask

##There are also many other useful maple routines in this file.

#################################################################
##Some directly related references:

#For symmetry of a refinable function and mask with a dilation matrix
#    see the papers

#    Bin Han, Symmetry property and construction of wavelets with a
#    general dilation matrix, Linear Algebra and Its Applications,
#    (2001), to appear.

#    Bin Han, Computing the smoothness exponent of a symmetric
#    multivariate refinable function, (2001), preprint

#    Bin Han, Thomas P.-Y. Yu, Bruce Piper,
##   Multivariate Refinable Hermite Interpolants, (2002), preprint.

#For CBC algorithms on constructing biorthogonal multiwavelets.
#  CBC stands for Coset By Coset
#  These program is mainly based on the following paper

#     Bin Han, Hermite interpolants and biorthogonal multiwavelets
#     with arbitrary order of vanishing moments (1999) SPIE
#     Proc. Vol. 3813, pp. 147--161. 
#     The CBC algorithm for the multivariate multiwavelets with
#     a general dilation matrix was established.

# Also, see the following papers for the CBC algorithm

#     Bin Han, Analysis and Construction of Optimal Multivariate
#     Biorthogonal Wavelets With Compact Support, SIAM Journal on
#     Mathematical Analysis, Vol. 31, No.2 (1999/2000), 274--304. 
#     CBC algorithm for interpolatory masks with dilation
#     matrix 2I_s was first introduced in this paper.

#     Di-Rong Chen, Bin Han and Sherman D. Riemenschneider,
#     Construction of Multivariate Biorthogonal Wavelets With
#     Arbitrary Vanishing Moments, Advances in Computational
#     Mathematics, Vol. 13 No. 2 (2000), 131-165. 
#     CBC alogrithm for the scalar case with a general dilation
#     matrix was established.

#     Bin Han and Rong-Qing Jia, Quincunx Fundamental Refinable
#     Functions and Quincunx Biorthogonal Wavelets, Mathematics of
#     Computation, Vol. 71, No. 237, (2002), 165--196.
#     CBC algorithm for quincunx wavelets was discussed

#      Bin Han, Construction of multivariate biorthogonal wavelets by
#      CBC algorithm. Wavelet analysis and multiresolution methods
#      (Urbana-Champaign, IL, 1999), 105--143, Lecture Notes in Pure
#      and Appl. Math., 212, Dekker, New York, 2000.
#      Survey paper on CBC algorithm.

#      Bin Han, Hermite interpolants and biorthogonal multiwavelets
#      with arbitrary order of vanishing moments (1999) SPIE
#      Proc. Vol. 3813, pp. 147--161. 
#      Survey paper on CBC algorithm for 1-dimensional multiwavelets.

#      Papers can be downloaded at
##          http://www.ualberta.ca/~bhan/publ.htm

# Program was developed by Bin Han at University of Alberta
# Version 1 on January 6, 2002.
# Initial Tests of the program have been done on January 6, 2002.
# Report bugs, mistakes, errors, comments, suggestions etc. to
# Bin Han, bhan@ualberta.ca, htpp://www.ualbert.ca/~bhan 
####################################################################

################################################################
##Main Maple 6.0 routines.
#################################################################
read "basicall":
default_D1dilation:=2:

##Find the support of a 1D matrix mask.
##Input: the matrix of polynomials of the matrix mask.
##Output: result -> support of mask=[supp[1],supp[2]] 
D1Support:=proc(poly)
    local i,j,multi,POLY,result:
    result:=Vector(2):
    multi:=Multiplicity(poly):
    for i from 1 to multi do
    for j from 1 to multi do
        POLY:=expand(poly[i,j]):
        result[1]:=max(result[1],degree(subs({x=1/x},POLY),x)):
        result[2]:=max(result[2],degree(POLY,x)):
    od:od:
    result[1]:=-result[1]:
    result:=result:
end proc:

##Print a 1D matrix mask using various formats.
##Input: the matrix of polynomials of the matrix mask.
##Output: Display the matrix mask on screen in various formats.
##        and output the matrix mask in array format
##choice<-1, Do not display, output matrix mask in numerical array
##           format (6 decimal places)
##choice=-1, Do not display, output matrix mask in exact value array
##choice=0,  Display the mask and ouput the exact value form,
##otherwise, Display and output numerical form (6 Digits)
D1PrintMask:=proc(poly,choice)
    local i,j,ii,supp,multi,POLY,val,printmask,result:
    multi:=Multiplicity(poly):
    supp:=D1Support(poly):
    result:=Array(supp[1]..supp[2],1..multi,1..multi):
    if (choice>=0) then
      printmask:=Matrix(multi,(supp[2]-supp[1]+1)*multi):
    fi:
    for i from 1 to multi do
    for j from 1 to multi do
        POLY:=expand(poly[i,j]):
        for ii from supp[1] to supp[2] do
            val:=coeff(POLY,x,ii):
            if ((choice<(-1))or(choice>0)) then val:=evalf(val,6): fi:
            result[ii,i,j]:=val:
            if (choice>=0) then
               printmask[i,(ii-supp[1])*multi+j]:=val:
            fi:
        od:
    od:od:
    if (choice>=0) then print(printmask,supp): fi:
    result:=result:
end proc:

##Print the element of a 1D matrix mask at poistion [xp].
##Input: the matrix of polynomials of the matrix mask, xp,yp
##Output: Display the element of the matrix mask on screen.
##choice>=0, exact format, otherwise, numerical to 6 digits
D1PrintElement:=proc(poly,xp,choice)
    local i,j,multi,POLY,val,result:
    multi:=Multiplicity(poly):
    result:=Matrix(multi,multi):
    for i from 1 to multi do
    for j from 1 to multi do
        POLY:=expand(poly[i,j]):
        val:=coeff(POLY,x,xp):
        if (choice<0) then val:=evalf(val,6): fi:
        result[i,j]:=val:
    od:od:
    print(result):
    result:=result:
end proc:

##Find the equations of the (bi)orthogonal relation:
##      Primalpoly*\bar{Dualpoly} is an interpolatory mask w.r.t
##      Dilation*Z.
##Input: primal_poly, dual_poly, dilation, orth, choice
##orth=0 means dual filter: the primal-dual relation
##orth=1 means orthogonal filter: lowpass-highpass relation
##choice=0, equations in set for exact dual relation
##otherwise, equations in set for dual relation without zero-position

D1BiorthEqs:=proc(primal_poly,dual_poly,dilation,orth,choice)
    local i,j,ii,m,supp,multi,POLY,result;
    multi:=Multiplicity(dual_poly):
    m:=abs(dilation): #absolute value of the dilation factor
    POLY:=Matrix(multi,multi):
    for i from 1 to multi do
    for j from 1 to multi do
	POLY[j,i]:=subs({x=1/x},dual_poly[i,j]):
    od:od:
    POLY:=`.`(primal_poly,POLY):
    supp:=D1Support(POLY):
    result:={}:
    if (choice=0) then
       for i from 1 to multi do POLY[i,i]:=POLY[i,i]-m*orth: od:
     else
       for i from 1 to multi do
       for j from 1 to multi do
           POLY[i,j]:=expand(POLY[i,j]-coeff(POLY[i,j],x,0)):
       od:od:
    fi:
    for i from 1 to multi do
    for j from 1 to multi do
    for ii from supp[1] to supp[2] do
        #if ii \in dilation*Z, then its coefficient must be 0
        if ((abs(ii) mod m)=0) then
           result:=result union { coeff(POLY[i,j],x,ii) }:
        fi:
    od:od:od:
    result:=result:
end proc:

##Find the moment array of a matrix mask.
##Input: poly,dilation,sr
##Output: result -> the array of the moments of a matrix mask:
##When choice=0, output the whole moments
## Ja[\mu]:=sum(Mask[j]*(M^(-1)j)^\mu, j)
## choice=otherwise, output the coset moments
## Ja[gep, \mu]:=sum(Mask[gep+Mj](M^(-1)gep+j)^\mu, j)

D1MomentOfMask:=proc(poly,dilation,sr,choice)
    local i0,i,j,ii,iii,m,multi,POLY,supp,val,mask,tmp,result:
    m:=abs(dilation):
    multi:=Multiplicity(poly):
    supp:=D1Support(poly):
    mask:=Array(supp[1]..supp[2],1..multi,1..multi):
    for i from 1 to multi do
    for j from 1 to multi do
        POLY:=expand(poly[i,j]):
        for ii from supp[1] to supp[2] do
            mask[ii,i,j]:=coeff(POLY,x,ii):
        od:
    od:od:
    #compute the moment array of a matrix mask
    if (choice=0) then
       result:=Array(0..(sr-1),1..multi,1..multi):
    else
       result:=Array(0..(m-1),0..(sr-1),1..multi,1..multi):
    fi:
    for iii from 0 to (sr-1) do #\mu=iii
    for ii from supp[1] to supp[2] do
        if(iii=0) then val:=1: else val:=(ii/dilation)^iii: fi:
        for i from 1 to multi do
        for j from 1 to multi do
            tmp:=val*mask[ii,i,j]/(iii!):
            if (choice=0) then
               result[iii,i,j]:=result[iii,i,j]+tmp:
            else
               i0:=(ii+abs(ii)*m) mod m:
               result[i0,iii,i,j]:=result[i0,iii,i,j]+tmp:
            fi:
        od:od:
    od:od:
    result:=result:
end proc:

##Find the equations for sum rules of order sr w.r.t.
##     Dilation*Z.
##Input: poly,dilation,sr,choice,knownvalue
##Output: result -> the equations for sum rules of order sr
##choice:=0, then use standard definition of sum rules for
##      scalar mask. That is, sum(Mask[i+Mj](i+Mj)^\mu, j\in Z^s)
##      is independent of the integer i. 
##choice:=1, then each coset equals the knownvalue. That is,
##      sum(Mask[i+Mj](M^(-1)i+j)^\mu, j)=knownvalue[\mu]. 
##choice=2, the y vector values are known for matrix mask.
##      that is, yvect[\mu]=knownvalue[\mu]
##otherwise, we output the equations of sum rules with
##      parameterized y vectors (in yvect)
##Do not assume that the sum of the mask is = |det(dilation)|
##The definition of sum rules is defined in Eqs. (1.6) in
##    Bin Han, Approximation properties and construction of
##    Hermite Interpolants and Biorthogonal Multiwavelets,
##    Journal of Approximation Theory, Vol. 110, (2001) 18-53

D1SumRuleEqs:=proc(poly,dilation,sr,knownvalue,choice)
    local i,j,k,ii,jj,m,gep,multi,POLY,supp,val,Ja,tmp,yvect,result;
    m:=abs(dilation):
    multi:=Multiplicity(poly):
    supp:=D1Support(poly):
    Ja:=D1MomentOfMask(poly,dilation,sr,1):
    result:=Array(0..sr):
    if ((multi=1)and((choice=0)or(choice=1))) then #multi=1
       for k from 0 to (sr-1) do #\mu=k
           result[k+1]:={}:
           for gep from 0 to (m-1) do
               if (choice=0) then
               result[k+1]:=result[k+1] union {Ja[gep,k,1,1]-Ja[0,k,1,1]}:
               else
               result[k+1]:=result[k+1] union {Ja[gep,k,1,1]-knownvalue[k]}:
               fi:
           od:
       od:
    else
       if (choice=2) then
          yvect:=knownvalue:
       else
          yvect:=Array(0..(sr-1),1..multi):
          jj:=1:
          for i from 0 to (sr-1) do #\mu=i
          for ii from 1 to multi do
              yvect[i,ii]:=evaln(yvect||jj):
              jj:=jj+1:
          od:od:
       fi:
       for k from 0 to (sr-1) do #\mu=k
           result[k+1]:={}:
       for gep from 0 to (m-1) do          
       for i from 1 to multi do
               tmp:=-((1/dilation)^k)*yvect[k,i]:
               for ii from 0 to k do #\nu=ii <= \mu
               for j from 1 to multi do
                   tmp:=tmp+(-1)^ii*Ja[gep,ii,j,i]*yvect[k-ii,j]:
               od:od:
               result[k+1]:=result[k+1] union {tmp}:
            od:
       od:od:
    fi:
    result[0]:={}:
    for i from 1 to (sr) do
        result[0]:=result[0] union convert(result[i],set):
    od:
    result:=result:
end proc:

##Given a matrix mask, find the \tilde y vector values by
##the CBC algorithm
##Input: poly,dilation,sr,choice
 ##Output: result -> the \tilde y vectors in the sum rules
##        in Theorem 3.1 in the following paper
##        Bin Han, Approximation properties and construction of
##        Hermite Interpolants and Biorthogonal Multiwavelets,
##        Journal of Approximation Theory, Vol. 110, (2001) 18-53
##Choice=0, output the moment for scalar mask:
##        h^a(\mu):=sum( \tilde a(i+Mj)(M^(-1)i+j)^\mu/\mu!, j)
##Choice=otherwise, output the \tilde y vectors in Theorem 3.1

D1CBCAlgorithm:=proc(poly,dilation,sr,choice)
    local i,ii,iii,j,k,m,multi,Ja,Ja0MinusmI,yvect,SE,tmp,tmp1,result:
    multi:=Multiplicity(poly):
    m:=abs(dilation):
    Ja:=D1MomentOfMask(poly,dilation,sr,0):
    if ((multi=1)and(choice=0)) then
       if (Ja[0,1,1]=0) then
          printf(`The sum of the primal mask must be nonzero!\n`):
          printf(`Program is terminated at subroutine D1CBCAlgorithm!\n`):
          quit:
       fi:          
       result:=Array(0..(sr-1)):
       result[0]:=1:
       for i from 1 to (sr-1) do #\mu=i
           tmp:=0:
           for ii from 0 to (i-1) do #\nu=ii < \mu
               tmp:=tmp+(-1)^(i-ii)*result[ii]*Ja[i-ii,1,1]:
           od:
           result[i]:=-tmp/Ja[0,1,1]:
       od:
    else
       result:=Array(0..(sr-1),1..multi):
       Ja0MinusmI:=matrix(multi,multi):
       for i from 1 to multi do
           for j from 1 to multi do
               Ja0MinusmI[i,j]:=-Ja[0,i,j]:
           od:
           Ja0MinusmI[i,i]:=Ja0MinusmI[i,i]+m:
       od:
       ##generate y0 vector here
       tmp1:=kernel(Ja0MinusmI, 'nulldimension');
       if (nulldimension <> 1) then
           printf(`There are something wrong with the mask!\n`):
           printf(`Program is terminated at subroutine D1CBCAlgorithm!\n`):
           quit;
       fi:
       tmp:=tmp1[1]:
       for i from 1 to multi do result[0,i]:=tmp[i]: od:

       #generate all other \tilde y vectors (they must be unique)
       for iii from 1 to (sr-1) do #\mu=iii
           for i from 1 to multi do
               Ja0MinusmI[i,i]:=dilation^(1+iii)-Ja[0,i,i]:
           od:
           SE:=inverse(Ja0MinusmI):
           for i from 1 to multi do
               tmp[i]:=0:
               for ii from 0 to (iii-1) do #\nu=ii < \mu
               for k from 1 to multi do
                   tmp[i]:=tmp[i]+dilation^(iii-ii)*Ja[iii-ii,i,k]*result[ii,k]:
               od:od:
           od:
           for i from 1 to multi do
               result[iii,i]:=0:
               for j from 1 to multi do
                   result[iii,i]:=result[iii,i]+SE[i,j]*tmp[j]:
               od:
           od:
       od:
    fi:
    result:=result:
end proc:

##Parametrize the mask using prescribed symmetry group {-1}.
##Input: sz,dilation,hermiteorder,interpolating
##Output: a matrix of polynomial for the parameterized matrix mask
##When interpolating=0, it is not an interpolanting mask
##w.r.t. dilation*Z, otherwise it is an interpolating matrix mask.
##cc is the name of the parameters: cc1, cc2, cc3,...
D1MaskOfSym:=proc(sz,dilation,hermiteorder,interpolating,cc)
    local i,j,ii,jj,iii,jjj,k,ind,invdilation,multi,m,mask,VAR,result:
    multi:=hermiteorder+1:
    mask:=Array((-sz)..sz,1..multi,1..multi,0):
    m:=abs(dilation):
    result:=Matrix(multi):  
    if (intepolating<>0) then
       for i from 1 to multi do
           mask[0,i,i]:=(1/dilation)^(i-1):
       od:
    fi:
    ind:=1: VAR:={}:
    for ii from 0 to sz do
        if((interpolating=0)or( ((ii+abs(ii)*m) mod m)<>0)) then
            for i from 1 to multi do
            for j from 1 to multi do
                mask[ii,i,j]:=evaln(cc||ind):
                mask[-ii,i,j]:=(-1)^(i+j)*mask[ii,i,j]:
		VAR:=VAR union { evaln(cc||ind) }:
                ind:=ind+1:
            od:od:
        fi:
     od:
     for i from 1 to multi do
     for j from 1 to multi do
         result[i,j]:=0:
         for ii from (-sz) to sz do
             result[i,j]:=result[i,j]+mask[ii,i,j]*(x^ii):
         od:
     od:od:
     save VAR, FreeParameter:
     result:=result:
end proc:

##Given hermiteorder and sum rule order, generate the y vector for
##refinable Hermite mask
D1YVectForHermiteMask:=proc(hmord,sr)
     local i,j,result:
     if(hmord=0) then
        result:=Array(0..(sr-1),0):
        result[0]:=1:
     else
        result:=Array(0..(sr-1),1..(hmord+1),0):
        for i from 0 to hmord do result[i,i+1]:=1: od:
     fi:
     result:=result:
end proc: