# Diff of /trunk/escript/py_src/generatediff

revision 291 by gross, Fri Dec 2 03:10:06 2005 UTC revision 439 by gross, Fri Jan 20 01:46:22 2006 UTC
# Line 33  def wherepos(arg): Line 33  def wherepos(arg):
33     else:     else:
34        return 0.        return 0.
35
36
37  class OPERATOR:  class OPERATOR:
38      def __init__(self,nickname,rng=[-1000.,1000],test_expr="",math_expr=None,      def __init__(self,nickname,rng=[-1000.,1000],test_expr="",math_expr=None,
39                  numarray_expr="",symbol_expr=None,diff=None,name=""):                  numarray_expr="",symbol_expr=None,diff=None,name=""):
# Line 172  def makeArray(shape,rng): Line 173  def makeArray(shape,rng):
173               for i2 in range(shape[2]):               for i2 in range(shape[2]):
174                  for i3 in range(shape[3]):                  for i3 in range(shape[3]):
175                     out[i0,i1,i2,i3]=l*random.random()+rng[0]                     out[i0,i1,i2,i3]=l*random.random()+rng[0]
176       elif len(shape)==5:
177           for i0 in range(shape[0]):
178              for i1 in range(shape[1]):
179                 for i2 in range(shape[2]):
180                    for i3 in range(shape[3]):
181                      for i4 in range(shape[4]):
182                       out[i0,i1,i2,i3,i4]=l*random.random()+rng[0]
183     else:     else:
184         raise SystemError,"rank is restricted to 4"         raise SystemError,"rank is restricted to 5"
185     return out             return out
186
187
# Line 589  def testTensorMult(arg0,arg1,sh_s): Line 597  def testTensorMult(arg0,arg1,sh_s):
597               for i3 in range(arg0.shape[3]):               for i3 in range(arg0.shape[3]):
598                       out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]                       out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
599          return out          return out
600
601    def testReduce(arg0,init_val,test_expr,post_expr):
602         out=init_val
603         if isinstance(arg0,float):
604              out=eval(test_expr.replace("%a1%","arg0"))
605         elif arg0.rank==0:
606              out=eval(test_expr.replace("%a1%","arg0"))
607         elif arg0.rank==1:
608            for i0 in range(arg0.shape[0]):
609                   out=eval(test_expr.replace("%a1%","arg0[i0]"))
610         elif arg0.rank==2:
611            for i0 in range(arg0.shape[0]):
612             for i1 in range(arg0.shape[1]):
613                   out=eval(test_expr.replace("%a1%","arg0[i0,i1]"))
614         elif arg0.rank==3:
615            for i0 in range(arg0.shape[0]):
616             for i1 in range(arg0.shape[1]):
617               for i2 in range(arg0.shape[2]):
618                   out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2]"))
619         elif arg0.rank==4:
620            for i0 in range(arg0.shape[0]):
621             for i1 in range(arg0.shape[1]):
622               for i2 in range(arg0.shape[2]):
623                 for i3 in range(arg0.shape[3]):
624                   out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2,i3]"))
625         return eval(post_expr)
626
627    def clipTEST(arg0,mn,mx):
628         if isinstance(arg0,float):
629              return max(min(arg0,mx),mn)
630         out=numarray.zeros(arg0.shape,numarray.Float64)
631         if arg0.rank==1:
632            for i0 in range(arg0.shape[0]):
633                out[i0]=max(min(arg0[i0],mx),mn)
634         elif arg0.rank==2:
635            for i0 in range(arg0.shape[0]):
636             for i1 in range(arg0.shape[1]):
637                out[i0,i1]=max(min(arg0[i0,i1],mx),mn)
638         elif arg0.rank==3:
639            for i0 in range(arg0.shape[0]):
640             for i1 in range(arg0.shape[1]):
641               for i2 in range(arg0.shape[2]):
642                  out[i0,i1,i2]=max(min(arg0[i0,i1,i2],mx),mn)
643         elif arg0.rank==4:
644            for i0 in range(arg0.shape[0]):
645             for i1 in range(arg0.shape[1]):
646               for i2 in range(arg0.shape[2]):
647                 for i3 in range(arg0.shape[3]):
648                    out[i0,i1,i2,i3]=max(min(arg0[i0,i1,i2,i3],mx),mn)
649         return out
650    def minimumTEST(arg0,arg1):
651         if isinstance(arg0,float):
652           if isinstance(arg1,float):
653              if arg0>arg1:
654                  return arg1
655              else:
656                  return arg0
657           else:
658              arg0=numarray.ones(arg1.shape)*arg0
659         else:
660           if isinstance(arg1,float):
661              arg1=numarray.ones(arg0.shape)*arg1
662         out=numarray.zeros(arg0.shape,numarray.Float64)
663         if arg0.rank==0:
664              if arg0>arg1:
665                  out=arg1
666              else:
667                  out=arg0
668         elif arg0.rank==1:
669            for i0 in range(arg0.shape[0]):
670              if arg0[i0]>arg1[i0]:
671                  out[i0]=arg1[i0]
672              else:
673                  out[i0]=arg0[i0]
674         elif arg0.rank==2:
675            for i0 in range(arg0.shape[0]):
676             for i1 in range(arg0.shape[1]):
677              if arg0[i0,i1]>arg1[i0,i1]:
678                  out[i0,i1]=arg1[i0,i1]
679              else:
680                  out[i0,i1]=arg0[i0,i1]
681         elif arg0.rank==3:
682            for i0 in range(arg0.shape[0]):
683             for i1 in range(arg0.shape[1]):
684               for i2 in range(arg0.shape[2]):
685                 if arg0[i0,i1,i2]>arg1[i0,i1,i2]:
686                  out[i0,i1,i2]=arg1[i0,i1,i2]
687                 else:
688                  out[i0,i1,i2]=arg0[i0,i1,i2]
689         elif arg0.rank==4:
690            for i0 in range(arg0.shape[0]):
691             for i1 in range(arg0.shape[1]):
692               for i2 in range(arg0.shape[2]):
693                 for i3 in range(arg0.shape[3]):
694                  if arg0[i0,i1,i2,i3]>arg1[i0,i1,i2,i3]:
695                   out[i0,i1,i2,i3]=arg1[i0,i1,i2,i3]
696                  else:
697                   out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
698         return out
699    def unrollLoops(a,b,o,arg,tap="",x="x"):
700        out=""
701        if a.rank==1:
702                 z=""
703                 for i99 in range(a.shape[0]):
704                   if not z=="": z+="+"
705                   if o=="1":
706                      z+="(%s)*%s[%s]"%(a[i99]+b[i99],x,i99)
707                   else:
708                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i99],x,i99,b[i99],x,i99)
709
710                 out+=tap+"%s=%s\n"%(arg,z)
711
712        elif a.rank==2:
713            for i0 in range(a.shape[0]):
714                 z=""
715                 for i99 in range(a.shape[1]):
716                   if not z=="": z+="+"
717                   if o=="1":
718                      z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)
719                   else:
720                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i99],x,i99,b[i0,i99],x,i99)
721
722                 out+=tap+"%s[%s]=%s\n"%(arg,i0,z)
723        elif a.rank==3:
724            for i0 in range(a.shape[0]):
725             for i1 in range(a.shape[1]):
726                 z=""
727                 for i99 in range(a.shape[2]):
728                   if not z=="": z+="+"
729                   if o=="1":
730                      z+="(%s)*%s[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],x,i99)
731                   else:
732                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i99],x,i99,b[i0,i1,i99],x,i99)
733
734                 out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)
735        elif a.rank==4:
736            for i0 in range(a.shape[0]):
737             for i1 in range(a.shape[1]):
738               for i2 in range(a.shape[2]):
739                 z=""
740                 for i99 in range(a.shape[3]):
741                   if not z=="": z+="+"
742                   if o=="1":
743                      z+="(%s)*%s[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],x,i99)
744                   else:
745                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i99],x,i99,b[i0,i1,i2,i99],x,i99)
746
747                 out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)
748        elif a.rank==5:
749            for i0 in range(a.shape[0]):
750             for i1 in range(a.shape[1]):
751               for i2 in range(a.shape[2]):
752                for i3 in range(a.shape[3]):
753                 z=""
754                 for i99 in range(a.shape[4]):
755                   if not z=="": z+="+"
756                   if o=="1":
757                      z+="(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],x,i99)
758                   else:
759                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99],x,i99,b[i0,i1,i2,i3,i99],x,i99)
760
761                 out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
762        return out
763
765        out=""
766        if a.rank==1:
767                 for i99 in range(a.shape[0]):
768                   if o=="1":
769                      out+=tap+"%s[%s]=(%s)\n"%(arg,i99,a[i99]+b[i99])
770                   else:
771                      out+=tap+"%s[%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i99,a[i99],i99,b[i99])
772
773        elif a.rank==2:
774            for i0 in range(a.shape[0]):
775                 for i99 in range(a.shape[1]):
776                   if o=="1":
777                      out+=tap+"%s[%s,%s]=(%s)\n"%(arg,i0,i99,a[i0,i99]+b[i0,i99])
778                   else:
779                      out+=tap+"%s[%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i99,a[i0,i99],i99,b[i0,i99])
780
781        elif a.rank==3:
782            for i0 in range(a.shape[0]):
783             for i1 in range(a.shape[1]):
784                 for i99 in range(a.shape[2]):
785                   if o=="1":
786                      out+=tap+"%s[%s,%s,%s]=(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99]+b[i0,i1,i99])
787                   else:
788                      out+=tap+"%s[%s,%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99],i99,b[i0,i1,i99])
789
790        elif a.rank==4:
791            for i0 in range(a.shape[0]):
792             for i1 in range(a.shape[1]):
793               for i2 in range(a.shape[2]):
794                 for i99 in range(a.shape[3]):
795                   if o=="1":
796                     out+=tap+"%s[%s,%s,%s,%s]=(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99]+b[i0,i1,i2,i99])
797                   else:
798                     out+=tap+"%s[%s,%s,%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99],i99,b[i0,i1,i2,i99])
799        return out
800    def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
801        if where=="Function":
802           xfac_o=1.
803           xfac_op=0.
804           z_fac=1./2.
805           z_fac_s=""
806           zo_fac_s="/(o+1.)"
807        elif where=="FunctionOnBoundary":
808           xfac_o=1.
809           xfac_op=0.
810           z_fac=1.
811           z_fac_s="*dim"
812           zo_fac_s="*(1+2.*(dim-1.)/(o+1.))"
813        elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
814           xfac_o=0.
815           xfac_op=1.
816           z_fac=1./2.
817           z_fac_s=""
818           zo_fac_s="/(o+1.)"
819        out=""
820        if a.rank==1:
821                 zo=0.
822                 zop=0.
823                 z=0.
824                 for i99 in range(a.shape[0]):
825                      if i99==0:
826                        zo+=       xfac_o*a[i99]
827                        zop+=       xfac_op*a[i99]
828                      else:
829                        zo+=a[i99]
830                      z+=b[i99]
831
832                 out+=tap+"%s=(%s)%s+(%s)%s"%(arg,zo,zo_fac_s,z*z_fac,z_fac_s)
833                 if zop==0.:
834                   out+="\n"
835                 else:
836                   out+="+(%s)*0.5**o\n"%zop
837        elif a.rank==2:
838            for i0 in range(a.shape[0]):
839                 zo=0.
840                 zop=0.
841                 z=0.
842                 for i99 in range(a.shape[1]):
843                      if i99==0:
844                        zo+=       xfac_o*a[i0,i99]
845                        zop+=       xfac_op*a[i0,i99]
846                      else:
847                        zo+=a[i0,i99]
848                      z+=b[i0,i99]
849
850                 out+=tap+"%s[%s]=(%s)%s+(%s)%s"%(arg,i0,zo,zo_fac_s,z*z_fac,z_fac_s)
851                 if zop==0.:
852                   out+="\n"
853                 else:
854                   out+="+(%s)*0.5**o\n"%zop
855        elif a.rank==3:
856            for i0 in range(a.shape[0]):
857             for i1 in range(a.shape[1]):
858                 zo=0.
859                 zop=0.
860                 z=0.
861                 for i99 in range(a.shape[2]):
862                      if i99==0:
863                        zo+=       xfac_o*a[i0,i1,i99]
864                        zop+=       xfac_op*a[i0,i1,i99]
865                      else:
866                        zo+=a[i0,i1,i99]
867                      z+=b[i0,i1,i99]
868
869                 out+=tap+"%s[%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,zo,zo_fac_s,z*z_fac,z_fac_s)
870                 if zop==0.:
871                   out+="\n"
872                 else:
873                   out+="+(%s)*0.5**o\n"%zop
874        elif a.rank==4:
875            for i0 in range(a.shape[0]):
876             for i1 in range(a.shape[1]):
877               for i2 in range(a.shape[2]):
878                 zo=0.
879                 zop=0.
880                 z=0.
881                 for i99 in range(a.shape[3]):
882                      if i99==0:
883                        zo+=       xfac_o*a[i0,i1,i2,i99]
884                        zop+=       xfac_op*a[i0,i1,i2,i99]
885
886                      else:
887                        zo+=a[i0,i1,i2,i99]
888                      z+=b[i0,i1,i2,i99]
889
890                 out+=tap+"%s[%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,zo,zo_fac_s,z*z_fac,z_fac_s)
891                 if zop==0.:
892                   out+="\n"
893                 else:
894                   out+="+(%s)*0.5**o\n"%zop
895
896        elif a.rank==5:
897            for i0 in range(a.shape[0]):
898             for i1 in range(a.shape[1]):
899               for i2 in range(a.shape[2]):
900                for i3 in range(a.shape[3]):
901                 zo=0.
902                 zop=0.
903                 z=0.
904                 for i99 in range(a.shape[4]):
905                      if i99==0:
906                        zo+=       xfac_o*a[i0,i1,i2,i3,i99]
907                        zop+=       xfac_op*a[i0,i1,i2,i3,i99]
908
909                      else:
910                        zo+=a[i0,i1,i2,i3,i99]
911                      z+=b[i0,i1,i2,i3,i99]
912                 out+=tap+"%s[%s,%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,i3,zo,zo_fac_s,z*z_fac,z_fac_s)
913                 if zop==0.:
914                   out+="\n"
915                 else:
916                   out+="+(%s)*0.5**o\n"%zop
917
918        return out
919
920  #=======================================================================================================  #=======================================================================================================
921  # tensor multiply  # interpolation
922  #=======================================================================================================  #=======================================================================================================
923  # oper=["generalTensorProduct",tensorProductTest]  for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne","Solution","ReducedSolution"]:
924  # oper=["matrixmult",testMatrixMult]    for data in ["Data","Symbol"]:
925  oper=["tensormult",testTensorMult]       for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
926          for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
927            if  where==case or \
928                ( case in ["ContinuousFunction","Solution","ReducedSolution"] and where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"] ) or \
929                ( case in ["FunctionOnContactZero","FunctionOnContactOne"] and where in ["FunctionOnContactZero","FunctionOnContactOne"] ) or \
930                (case=="ContinuousFunction" and  where in ["Solution","ReducedSolution"]) or \
931                (case=="Solution" and  where=="ReducedSolution") :
932
933
934             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
935             tname="test_interpolation_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
936             text+="   def %s(self):\n"%tname
937             text+="      \"\"\"\n"
938             text+="      tests interpolation for rank %s %s onto the %s\n\n"%(len(sh),data,where)
939             text+="      assumptions: self.domain supports inpterpolation from %s onto %s\n"%(case,where)
940             text+="      \"\"\"\n"
941             if case=="ReducedSolution" or where=="ReducedSolution":
942                text+="      o=1\n"
943                o="1"
944             else:
945                text+="      o=self.order\n"
946                o="o"
947             text+="      dim=self.domain.getDim()\n"
948             text+="      w_ref=%s(self.domain)\n"%where
949             text+="      x_ref=w_ref.getX()\n"
950             text+="      w=%s(self.domain)\n"%case
951             text+="      x=w.getX()\n"
952             a_2=makeArray(sh+(2,),[-1.,1])
953             b_2=makeArray(sh+(2,),[-1.,1])
954             a_3=makeArray(sh+(3,),[-1.,1])
955             b_3=makeArray(sh+(3,),[-1.,1])
956             if data=="Symbol":
957                text+="      arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
958                val="s"
959                res="sub"
960             else:
961                val="arg"
962                res="res"
963             text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
964             text+="      ref=Data(0,%s,w_ref)\n"%str(sh)
965             text+="      if dim==2:\n"
966             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
967             text+=unrollLoops(a_2,b_2,o,"ref",tap="        ",x="x_ref")
968             text+="      else:\n"
969
970             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
971             text+=unrollLoops(a_3,b_3,o,"ref",tap="        ",x="x_ref")
972             text+="      res=interpolate(arg,where=w_ref)\n"
973             if data=="Symbol":
974                text+="      sub=res.substitute({arg:s})\n"
975             text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
976             text+="      self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong functionspace of result.\")\n"%res
977             text+="      self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
978             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
979             t_prog+=text
981    print t_prog
982    print test_tail
983    1/0
984
985  for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:  #=======================================================================================================
987    #=======================================================================================================
988    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
989      for data in ["Data","Symbol"]:
990         for case in ["ContinuousFunction","Solution","ReducedSolution"]:
991           for sh in [ (),(2,), (4,5), (6,2,2)]:
992             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
994             text+="   def %s(self):\n"%tname
995             text+="      \"\"\"\n"
996             text+="      tests gradient for rank %s %s on the %s\n\n"%(len(sh),data,where)
997             text+="      assumptions: %s(self.domain) exists\n"%case
998             text+="                   self.domain supports gardient on %s\n"%where
999             text+="      \"\"\"\n"
1000             if case=="ReducedSolution":
1001                text+="      o=1\n"
1002                o="1"
1003             else:
1004                text+="      o=self.order\n"
1005                o="o"
1006             text+="      dim=self.domain.getDim()\n"
1007             text+="      w_ref=%s(self.domain)\n"%where
1008             text+="      x_ref=w_ref.getX()\n"
1009             text+="      w=%s(self.domain)\n"%case
1010             text+="      x=w.getX()\n"
1011             a_2=makeArray(sh+(2,),[-1.,1])
1012             b_2=makeArray(sh+(2,),[-1.,1])
1013             a_3=makeArray(sh+(3,),[-1.,1])
1014             b_3=makeArray(sh+(3,),[-1.,1])
1015             if data=="Symbol":
1016                text+="      arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1017                val="s"
1018                res="sub"
1019             else:
1020                val="arg"
1021                res="res"
1022             text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
1023             text+="      ref=Data(0,%s+(dim,),w_ref)\n"%str(sh)
1024             text+="      if dim==2:\n"
1025             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1027             text+="      else:\n"
1028
1029             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1032             if data=="Symbol":
1033                text+="      sub=res.substitute({arg:s})\n"
1034             text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1035             text+="      self.failUnlessEqual(res.getShape(),%s+(dim,),\"wrong shape of result.\")\n"%str(sh)
1036             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1037
1038
1039             t_prog+=text
1041    print t_prog
1042    print test_tail
1043    1/0
1044
1045
1046    #=======================================================================================================
1047    # integrate
1048    #=======================================================================================================
1049    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1050      for data in ["Data","Symbol"]:
1051        for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1052          for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1053            if (not case in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]) or where==case:
1054             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1055             tname="test_integrate_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1056             text+="   def %s(self):\n"%tname
1057             text+="      \"\"\"\n"
1058             text+="      tests integral of rank %s %s on the %s\n\n"%(len(sh),data,where)
1059             text+="      assumptions: %s(self.domain) exists\n"%case
1060             text+="                   self.domain supports integral on %s\n"%where
1061
1062             text+="      \"\"\"\n"
1063             if case=="ReducedSolution":
1064                text+="      o=1\n"
1065                o="1"
1066             else:
1067                text+="      o=self.order\n"
1068                o="o"
1069             text+="      dim=self.domain.getDim()\n"
1070             text+="      w_ref=%s(self.domain)\n"%where
1071             text+="      w=%s(self.domain)\n"%case
1072             text+="      x=w.getX()\n"
1073             a_2=makeArray(sh+(2,),[-1.,1])
1074             b_2=makeArray(sh+(2,),[-1.,1])
1075             a_3=makeArray(sh+(3,),[-1.,1])
1076             b_3=makeArray(sh+(3,),[-1.,1])
1077             if data=="Symbol":
1078                text+="      arg=Symbol(shape=%s)\n"%str(sh)
1079                val="s"
1080                res="sub"
1081             else:
1082                val="arg"
1083                res="res"
1084
1085             text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
1086             if not len(sh)==0:
1087                text+="      ref=numarray.zeros(%s,numarray.Float)\n"%str(sh)
1088             text+="      if dim==2:\n"
1089             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1090             text+=unrollLoopsOfInteriorIntegral(a_2,b_2,where,"ref",tap="        ")
1091             text+="      else:\n"
1092
1093             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1094             text+=unrollLoopsOfInteriorIntegral(a_3,b_3,where,"ref",tap="        ")
1095             if case in ["ContinuousFunction","Solution","ReducedSolution"]:
1096                 text+="      res=integrate(arg,where=w_ref)\n"
1097             else:
1098                 text+="      res=integrate(arg)\n"
1099
1100             if data=="Symbol":
1101                text+="      sub=res.substitute({arg:s})\n"
1102             if len(sh)==0 and data=="Data":
1103                text+="      self.failUnless(isinstance(%s,float),\"wrong type of result.\")\n"%res
1104             else:
1105                if data=="Symbol":
1106                   text+="      self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1107                   text+="      self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1108                else:
1109                   text+="      self.failUnless(isinstance(res,numarray.NumArray),\"wrong type of result.\")\n"
1110                   text+="      self.failUnlessEqual(res.shape,%s,\"wrong shape of result.\")\n"%str(sh)
1111             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1112
1113
1114             t_prog+=text
1116    print t_prog
1117    print test_tail
1118    1/0
1119    #=======================================================================================================
1120    # inverse
1121    #=======================================================================================================
1122    name="inverse"
1123    for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1124      for sh0 in [ (1,1), (2,2), (3,3)]:
1125                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1126                  tname="test_%s_%s_dim%s"%(name,case0,sh0[0])
1127                  text+="   def %s(self):\n"%tname
1128                  a_0=makeArray(sh0,[-1.,1])
1129                  for i in range(sh0[0]): a_0[i,i]+=2
1130                  if case0 in ["taggedData", "expandedData"]:
1131                      a1_0=makeArray(sh0,[-1.,1])
1132                      for i in range(sh0[0]): a1_0[i,i]+=3
1133                  else:
1134                      a1_0=a_0
1135
1136                  text+=mkText(case0,"arg",a_0,a1_0)
1137                  text+="      res=%s(arg)\n"%name
1138                  if case0=="Symbol":
1139                     text+=mkText("array","s",a_0,a1_0)
1140                     text+="      sub=res.substitute({arg:s})\n"
1141                     res="sub"
1142                     ref="s"
1143                  else:
1144                     ref="arg"
1145                     res="res"
1146                  text+=mkTypeAndShapeTest(case0,sh0,"res")
1147                  text+="      self.failUnless(Lsup(matrixmult(%s,%s)-kronecker(%s))<=self.RES_TOL,\"wrong result\")\n"%(res,ref,sh0[0])
1148
1149                  if case0 == "taggedData" :
1150                      t_prog_with_tags+=text
1151                  else:
1152                      t_prog+=text
1153
1155    # print t_prog
1156    print t_prog_with_tags
1157    print test_tail
1158    1/0
1159
1160    #=======================================================================================================
1161    # trace
1162    #=======================================================================================================
1163    def traceTest(r,offset):
1164        sh=r.shape
1165        r1=1
1166        for i in range(offset): r1*=sh[i]
1167        r2=1
1168        for i in range(offset+2,len(sh)): r2*=sh[i]
1169        r_s=numarray.reshape(r,(r1,sh[offset],sh[offset],r2))
1170        s=numarray.zeros([r1,r2],numarray.Float)
1171        for i1 in range(r1):
1172            for i2 in range(r2):
1173                for j in range(sh[offset]): s[i1,i2]+=r_s[i1,j,j,i2]
1174        return s.resize(sh[:offset]+sh[offset+2:])
1175    name,tt="trace",traceTest
1176    for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1177      for sh0 in [ (4,5), (6,2,2),(3,2,3,4)]:
1178        for offset in range(len(sh0)-1):
1179                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1180                  tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1181                  text+="   def %s(self):\n"%tname
1182                  sh_t=list(sh0)
1183                  sh_t[offset+1]=sh_t[offset]
1184                  sh_t=tuple(sh_t)
1185                  sh_r=[]
1186                  for i in range(offset): sh_r.append(sh0[i])
1187                  for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])
1188                  sh_r=tuple(sh_r)
1189                  a_0=makeArray(sh_t,[-1.,1])
1190                  if case0 in ["taggedData", "expandedData"]:
1191                      a1_0=makeArray(sh_t,[-1.,1])
1192                  else:
1193                      a1_0=a_0
1194                  r=tt(a_0,offset)
1195                  r1=tt(a1_0,offset)
1196                  text+=mkText(case0,"arg",a_0,a1_0)
1197                  text+="      res=%s(arg,%s)\n"%(name,offset)
1198                  if case0=="Symbol":
1199                     text+=mkText("array","s",a_0,a1_0)
1200                     text+="      sub=res.substitute({arg:s})\n"
1201                     res="sub"
1202                     text+=mkText("array","ref",r,r1)
1203                  else:
1204                     res="res"
1205                     text+=mkText(case0,"ref",r,r1)
1206                  text+=mkTypeAndShapeTest(case0,sh_r,"res")
1207                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1208
1209                  if case0 == "taggedData" :
1210                      t_prog_with_tags+=text
1211                  else:
1212                      t_prog+=text
1213
1215    # print t_prog
1216    print t_prog_with_tags
1217    print test_tail
1218    1/0
1219
1220    #=======================================================================================================
1221    # clip
1222    #=======================================================================================================
1223    oper_L=[["clip",clipTEST]]
1224    for oper in oper_L:
1225     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1226    for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:    for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1227            if len(sh0)==0 or not case0=="float":
1228                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1229                  tname="test_%s_%s_rank%s"%(oper[0],case0,len(sh0))
1230                  text+="   def %s(self):\n"%tname
1231                  a_0=makeArray(sh0,[-1.,1])
1232                  if case0 in ["taggedData", "expandedData"]:
1233                      a1_0=makeArray(sh0,[-1.,1])
1234                  else:
1235                      a1_0=a_0
1236
1237                  r=oper[1](a_0,-0.3,0.5)
1238                  r1=oper[1](a1_0,-0.3,0.5)
1239                  text+=mkText(case0,"arg",a_0,a1_0)
1240                  text+="      res=%s(arg,-0.3,0.5)\n"%oper[0]
1241                  if case0=="Symbol":
1242                     text+=mkText("array","s",a_0,a1_0)
1243                     text+="      sub=res.substitute({arg:s})\n"
1244                     res="sub"
1245                     text+=mkText("array","ref",r,r1)
1246                  else:
1247                     res="res"
1248                     text+=mkText(case0,"ref",r,r1)
1249                  text+=mkTypeAndShapeTest(case0,sh0,"res")
1250                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1251
1252                  if case0 == "taggedData" :
1253                      t_prog_with_tags+=text
1254                  else:
1255                      t_prog+=text
1256
1258    # print t_prog
1259    print t_prog_with_tags
1260    print test_tail
1261    1/0
1262
1263    #=======================================================================================================
1264    # maximum, minimum, clipping
1265    #=======================================================================================================
1266    oper_L=[ ["maximum",maximumTEST],
1267             ["minimum",minimumTEST]]
1268    for oper in oper_L:
1269     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1270      for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1271     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1272       for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:       for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1273         for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:          if (len(sh0)==0 or not case0=="float") and (len(sh1)==0 or not case1=="float") \
1274            if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \             and (sh0==sh1 or len(sh0)==0 or len(sh1)==0) :
and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
# if len(sh_s)==1 and len(sh0+sh_s)==2 and (len(sh_s+sh1)==1 or len(sh_s+sh1)==2)): # test for matrixmult
if ( len(sh_s)==1 and len(sh0+sh_s)==2 and ( len(sh1+sh_s)==2 or len(sh1+sh_s)==1 )) or (len(sh_s)==2 and len(sh0+sh_s)==4 and (len(sh1+sh_s)==2 or len(sh1+sh_s)==3 or len(sh1+sh_s)==4)):  # test for tensormult
case=getResultCaseForBin(case0,case1)
1275                use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"                use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1276
1277                text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"                text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1278                # tname="test_generalTensorProduct_%s_rank%s_%s_rank%s_offset%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1),len(sh_s))                tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
#tname="test_matrixmult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
tname="test_tensormult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
# if tname=="test_generalTensorProduct_array_rank1_array_rank2_offset1":
# print tnametest_generalTensorProduct_Symbol_rank1_Symbol_rank3_offset1
1279                text+="   def %s(self):\n"%tname                text+="   def %s(self):\n"%tname
1280                a_0=makeArray(sh0+sh_s,[-1.,1])                a_0=makeArray(sh0,[-1.,1])
1281                if case0 in ["taggedData", "expandedData"]:                if case0 in ["taggedData", "expandedData"]:
1282                    a1_0=makeArray(sh0+sh_s,[-1.,1])                    a1_0=makeArray(sh0,[-1.,1])
1283                else:                else:
1284                    a1_0=a_0                    a1_0=a_0
1285
1286                a_1=makeArray(sh_s+sh1,[-1.,1])                a_1=makeArray(sh1,[-1.,1])
1287                if case1 in ["taggedData", "expandedData"]:                if case1 in ["taggedData", "expandedData"]:
1288                    a1_1=makeArray(sh_s+sh1,[-1.,1])                    a1_1=makeArray(sh1,[-1.,1])
1289                else:                else:
1290                    a1_1=a_1                    a1_1=a_1
1291                r=oper[1](a_0,a_1,sh_s)                r=oper[1](a_0,a_1)
1292                r1=oper[1](a1_0,a1_1,sh_s)                r1=oper[1](a1_0,a1_1)
1293                text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)                text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1294                text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)                text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1295                #text+="      res=matrixmult(arg0,arg1)\n"                text+="      res=%s(arg0,arg1)\n"%oper[0]
1296                text+="      res=tensormult(arg0,arg1)\n"                case=getResultCaseForBin(case0,case1)
#text+="      res=generalTensorProduct(arg0,arg1,offset=%s)\n"%(len(sh_s))
1297                if case=="Symbol":                if case=="Symbol":
1298                   c0_res,c1_res=case0,case1                   c0_res,c1_res=case0,case1
1299                   subs="{"                   subs="{"
# Line 651  for case0 in ["float","array","Symbol"," Line 1313  for case0 in ["float","array","Symbol","
1313                else:                else:
1314                   res="res"                   res="res"
1315                   text+=mkText(case,"ref",r,r1)                   text+=mkText(case,"ref",r,r1)
1316                text+=mkTypeAndShapeTest(case,sh0+sh1,"res")                if len(sh0)>len(sh1):
1317                      text+=mkTypeAndShapeTest(case,sh0,"res")
1318                  else:
1319                      text+=mkTypeAndShapeTest(case,sh1,"res")
1320                text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res                text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1321
1322                if case0 == "taggedData" or case1 == "taggedData":                if case0 == "taggedData" or case1 == "taggedData":
1323                    t_prog_with_tags+=text                    t_prog_with_tags+=text
1324                else:                              else:
1325                    t_prog+=text                    t_prog+=text
1326
1328  # print t_prog  # print t_prog
1329  print t_prog_with_tags  print t_prog_with_tags
1330  print test_tail            print test_tail
1331  1/0  1/0
1332
1333
1334  #=======================================================================================================  #=======================================================================================================
1335  # outer/inner  # outer inner
1336  #=======================================================================================================  #=======================================================================================================
1337  oper=["inner",innerTEST]  oper=["outer",outerTEST]
1338  # oper=["outer",outerTEST]  # oper=["inner",innerTEST]
1339  for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:  for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1340    for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:    for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1341     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
# Line 725  for case0 in ["float","array","Symbol"," Line 1394  for case0 in ["float","array","Symbol","
1395  # print t_prog  # print t_prog
1396  print t_prog_with_tags  print t_prog_with_tags
1397    print test_tail
1398    1/0
1399
1400    #=======================================================================================================
1401    # local reduction
1402    #=======================================================================================================
1403    for oper in [["length",0.,"out+%a1%**2","math.sqrt(out)"],
1404                 ["maxval",-1.e99,"max(out,%a1%)","out"],
1405                 ["minval",1.e99,"min(out,%a1%)","out"] ]:
1406      for case in case_set:
1407         for sh in shape_set:
1408           if not case=="float" or len(sh)==0:
1409             text=""
1410             text+="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1411             tname="def test_%s_%s_rank%s"%(oper[0],case,len(sh))
1412             text+="   %s(self):\n"%tname
1413             a=makeArray(sh,[-1.,1.])
1414             a1=makeArray(sh,[-1.,1.])
1415             r1=testReduce(a1,oper[1],oper[2],oper[3])
1416             r=testReduce(a,oper[1],oper[2],oper[3])
1417
1418             text+=mkText(case,"arg",a,a1)
1419             text+="      res=%s(arg)\n"%oper[0]
1420             if case=="Symbol":
1421                 text+=mkText("array","s",a,a1)
1422                 text+="      sub=res.substitute({arg:s})\n"
1423                 text+=mkText("array","ref",r,r1)
1424                 res="sub"
1425             else:
1426                 text+=mkText(case,"ref",r,r1)
1427                 res="res"
1428             if oper[0]=="length":
1429                   text+=mkTypeAndShapeTest(case,(),"res")
1430             else:
1431                if case=="float" or case=="array":
1432                   text+=mkTypeAndShapeTest("float",(),"res")
1433                else:
1434                   text+=mkTypeAndShapeTest(case,(),"res")
1435             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1436             if case == "taggedData":
1437               t_prog_with_tags+=text
1438             else:
1439               t_prog+=text
1441    # print t_prog
1442    print t_prog_with_tags
1443    print test_tail
1444    1/0
1445
1446    #=======================================================================================================
1447    # tensor multiply
1448    #=======================================================================================================
1449    # oper=["generalTensorProduct",tensorProductTest]
1450    # oper=["matrixmult",testMatrixMult]
1451    oper=["tensormult",testTensorMult]
1452
1453    for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1454      for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1455       for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1456         for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1457           for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:
1458              if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \
1459                   and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
1460                # if len(sh_s)==1 and len(sh0+sh_s)==2 and (len(sh_s+sh1)==1 or len(sh_s+sh1)==2)): # test for matrixmult
1461                if ( len(sh_s)==1 and len(sh0+sh_s)==2 and ( len(sh1+sh_s)==2 or len(sh1+sh_s)==1 )) or (len(sh_s)==2 and len(sh0+sh_s)==4 and (len(sh1+sh_s)==2 or len(sh1+sh_s)==3 or len(sh1+sh_s)==4)):  # test for tensormult
1462                  case=getResultCaseForBin(case0,case1)
1463                  use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1464                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1465                  # tname="test_generalTensorProduct_%s_rank%s_%s_rank%s_offset%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1),len(sh_s))
1466                  #tname="test_matrixmult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
1467                  tname="test_tensormult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
1468                  # if tname=="test_generalTensorProduct_array_rank1_array_rank2_offset1":
1469                  # print tnametest_generalTensorProduct_Symbol_rank1_Symbol_rank3_offset1
1470                  text+="   def %s(self):\n"%tname
1471                  a_0=makeArray(sh0+sh_s,[-1.,1])
1472                  if case0 in ["taggedData", "expandedData"]:
1473                      a1_0=makeArray(sh0+sh_s,[-1.,1])
1474                  else:
1475                      a1_0=a_0
1476
1477                  a_1=makeArray(sh_s+sh1,[-1.,1])
1478                  if case1 in ["taggedData", "expandedData"]:
1479                      a1_1=makeArray(sh_s+sh1,[-1.,1])
1480                  else:
1481                      a1_1=a_1
1482                  r=oper[1](a_0,a_1,sh_s)
1483                  r1=oper[1](a1_0,a1_1,sh_s)
1484                  text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1485                  text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1486                  #text+="      res=matrixmult(arg0,arg1)\n"
1487                  text+="      res=tensormult(arg0,arg1)\n"
1488                  #text+="      res=generalTensorProduct(arg0,arg1,offset=%s)\n"%(len(sh_s))
1489                  if case=="Symbol":
1490                     c0_res,c1_res=case0,case1
1491                     subs="{"
1492                     if case0=="Symbol":
1493                        text+=mkText("array","s0",a_0,a1_0)
1494                        subs+="arg0:s0"
1495                        c0_res="array"
1496                     if case1=="Symbol":
1497                        text+=mkText("array","s1",a_1,a1_1)
1498                        if not subs.endswith("{"): subs+=","
1499                        subs+="arg1:s1"
1500                        c1_res="array"
1501                     subs+="}"
1502                     text+="      sub=res.substitute(%s)\n"%subs
1503                     res="sub"
1504                     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1505                  else:
1506                     res="res"
1507                     text+=mkText(case,"ref",r,r1)
1508                  text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
1509                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1510                  if case0 == "taggedData" or case1 == "taggedData":
1511                      t_prog_with_tags+=text
1512                  else:
1513                      t_prog+=text