###########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

##===================================================================
#This group of Maple procedures is supposed to be the general routine
#that does not depend on the dimensions. The routines here are basic
#to all mapel program and therefore are called by other programs

###########Header files###########################################
#FreeParameter(solt,newpara,choice) 
#        -> determine the free parameter in a set.
# ConvertToPoly(EQ)
#        -> if the equations in EQ are rational functions
#           the replace it with its numerator. 
#ConvertToSingular(EQ, filename)
#        -> convert EQ into a format used by Singular.
#           note: the parameters must in the form c(j), j=0,1,2,...
#SimplifySet(solt)
#        -> simplify the elements in a set
#ConvertSetIntoSum(solt)
#        -> ouput the sum of the square of the elements of the set
###################################################################

#################################################################
##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 
####################################################################

with(linalg):
with(LinearAlgebra):

##Find out all the free parameter in a solution set
##The criterion used here is to find left=right, that is the left side
##and the right side are the same
##When choice=0, do not do substitution
##otherwise, replace the free parameters with the ordered parameters:
##c1,c2,c3,c4,c5,...
##So in solt, the parameter should NOT use c.number

FreeParameter:=proc(solt,newpara,choice)
local i,ind,para,var:
var:={}: ind:=0:
if(solt={}) then
  printf("It is the empty set!\n"):
else
  for i from 1 to nops(solt) do
      para:=op(1,solt[i]):
      if( para=op(2,solt[i])) then
         ind:=ind+1:
         if (choice<>0) then
            var:=var union { para=evaln(newpara||ind) }:
         else
            var:=var union { para }:
         fi:
      fi:
  od:
fi:
print(var):
save var, VARIABLE:
printf("\nTotal number of free parameters is %d\n", ind):
var:=var:
end;

ConvertToPoly:=proc(solt)
local i,poly, SOLT, DENOM;
SOLT:={}:
DENOM:={}:
if(solt={}) then
  printf(`It is the empty set!\n`):
else
  for i from 1 to nops(solt) do
      poly:=expand(solt[i]):
      poly:=convert(poly, fraction):
      SOLT:=SOLT union {numer(poly)}:
      DENOM:=DENOM union {denom(poly)}:
  od:
fi:
print("The denominators are:", DENOM):
SOLT:=SOLT:
end;


##Find the multiplicity of a matrix mask
##Input: the matrix of polynomials of the matrix mask.
##Output: result -> the multiplicity of the matrix mask.
##If poly is not of type Matrix, then output 0.
Multiplicity:=proc(poly)
    local tmp,result:
    tmp:=op(1,poly): 
    if ( tmp[1]=tmp[2] ) then
         result:=tmp[1]:
    else
         printf(`Elements of a mask must be square matrices!\n`):
         printf(`Program is terminated!\n`):
         quit:
    fi:
    result:=result:
end proc:

##fresh the matrix mask
FreshMask:=proc(poly)
   local supp,i,j,result:
   supp:=Multiplicity(poly):
   result:=Array(supp[1]..supp[2],supp[3]..supp[4],0):
   for i from supp[1] to supp[2] do
   for j from supp[3] to supp[4] do
       result[i,j]:=1*poly[i,j]:
   od:od:
   result:=result:
end proc:

SimplifySet:=proc(solt)
local i,poly, SOLT, DENOM;
SOLT:={}:
DENOM:={}:
if(solt={}) then
  printf(`It is the empty set!\n`):
else
  for i from 1 to nops(solt) do
      poly:=simplify(convert(solt[i], fraction)):
      SOLT:=SOLT union {poly}:
  od:
fi:
SOLT:=SOLT:
end;


##output the sum of the squares of the elements in a set
ConvertSetIntoSum:=proc(solt)
local i,poly;
poly:=0:
if(solt={}) then
  printf(`It is the empty set!\n`):
else
  for i from 1 to nops(solt) do
      poly:=poly+solt[i]*solt[i]:
  od:
fi:
poly:=expand(poly):
end;



PolyOrder:=(A,B)->evalb( (nops(A) < nops(B))):

ConvertToSingular:=proc(EQ, filename)
local i, GEQ,LEQ;
GEQ:=EQ minus {0}:
LEQ:=convert(GEQ, list):
LEQ:=sort(LEQ, PolyOrder):

interface(quiet=true);
writeto(filename):
for i from 1 to nops(LEQ) do
    printf("%s,\n",convert(LEQ[i],string)):
od:
printf("0;\n"):
writeto(terminal):
end;

OutputInMatlab:=proc(EQ, filename, whichform)
local i,len,poly,total, GEQ,LEQ;
GEQ:=EQ minus {0}:
LEQ:=convert(GEQ, list):
LEQ:=sort(LEQ, PolyOrder):
len:=nops(LEQ):
total:=array(1..len):
interface(quiet=true);
writeto(filename):
if (whichform=0) then
   printf("function val=myfunc(c)\n"):
   for i from 1 to len do
      poly:=expand(evalf(LEQ[i],32)):
       printf("val(%d)=%s;\n",i,convert(poly,string)):
   od:
fi:
if (whichform=1) then
   printf("function val=myfunc(c)\n"):
   for i from 1 to len do
       total[i]:=expand(evalf(LEQ[i],32)):
       printf("total(%d)=%s;\n",i,convert(total[i],string)):
   od:
   printf("val="):
   for i from 1 to len-1 do
       printf("total(%d)*total(%d)+",i,i):
   od:
   printf("total(%d)*total(%d);\n",len,len):
fi:
writeto(terminal):
end;

AbsmaxInSet:=proc(solt)
local Max, i, val;
Max:=-1.0:
for i from 1 to nops(solt) do
    val:=abs(evalf(op(i,solt), 32)):
    if(val>Max) then
       Max:=val:
    fi:
od:
printf(`The absolution max in a set is %16.12f\n`, Max):
Max:=Max:
end:

ConvertArrayIntoSet:=proc(A)
      local len,i,result:
      len:=nops(A):
      result:={}:
      for i from 1 to len do
          result:=result union {convert(A[i], set)}:
      od:
      print(len,result):
      result:=result:
end proc: