/[escript]/trunk/escript/py_src/generatediff
ViewVC logotype

Diff of /trunk/escript/py_src/generatediff

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 313 by gross, Mon Dec 5 07:01:36 2005 UTC revision 443 by gross, Fri Jan 20 06:22:38 2006 UTC
# Line 173  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 616  def testReduce(arg0,init_val,test_expr,p Line 623  def testReduce(arg0,init_val,test_expr,p
623               for i3 in range(arg0.shape[3]):               for i3 in range(arg0.shape[3]):
624                 out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2,i3]"))                           out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2,i3]"))          
625       return eval(post_expr)       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        
700    def unrollLoops(a,b,o,arg,tap="",x="x"):
701        out=""
702        if a.rank==1:
703                 z=""
704                 for i99 in range(a.shape[0]):
705                   if not z=="": z+="+"
706                   if o=="1":
707                      z+="(%s)*%s[%s]"%(a[i99]+b[i99],x,i99)
708                   else:
709                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i99],x,i99,b[i99],x,i99)
710    
711                 out+=tap+"%s=%s\n"%(arg,z)
712    
713        elif a.rank==2:
714            for i0 in range(a.shape[0]):
715                 z=""
716                 for i99 in range(a.shape[1]):
717                   if not z=="": z+="+"
718                   if o=="1":
719                      z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)
720                   else:
721                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i99],x,i99,b[i0,i99],x,i99)
722    
723                 out+=tap+"%s[%s]=%s\n"%(arg,i0,z)
724        elif a.rank==3:
725            for i0 in range(a.shape[0]):
726             for i1 in range(a.shape[1]):
727                 z=""
728                 for i99 in range(a.shape[2]):
729                   if not z=="": z+="+"
730                   if o=="1":
731                      z+="(%s)*%s[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],x,i99)
732                   else:
733                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i99],x,i99,b[i0,i1,i99],x,i99)
734    
735                 out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)
736        elif a.rank==4:
737            for i0 in range(a.shape[0]):
738             for i1 in range(a.shape[1]):
739               for i2 in range(a.shape[2]):
740                 z=""
741                 for i99 in range(a.shape[3]):
742                   if not z=="": z+="+"
743                   if o=="1":
744                      z+="(%s)*%s[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],x,i99)
745                   else:
746                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i99],x,i99,b[i0,i1,i2,i99],x,i99)
747    
748                 out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)
749        elif a.rank==5:
750            for i0 in range(a.shape[0]):
751             for i1 in range(a.shape[1]):
752               for i2 in range(a.shape[2]):
753                for i3 in range(a.shape[3]):
754                 z=""
755                 for i99 in range(a.shape[4]):
756                   if not z=="": z+="+"
757                   if o=="1":
758                      z+="(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],x,i99)
759                   else:
760                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99],x,i99,b[i0,i1,i2,i3,i99],x,i99)
761    
762                 out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
763        return out
764    
765    def unrollLoopsOfGrad(a,b,o,arg,tap=""):
766        out=""
767        if a.rank==1:
768                 for i99 in range(a.shape[0]):
769                   if o=="1":
770                      out+=tap+"%s[%s]=(%s)\n"%(arg,i99,a[i99]+b[i99])
771                   else:
772                      out+=tap+"%s[%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i99,a[i99],i99,b[i99])
773    
774        elif a.rank==2:
775            for i0 in range(a.shape[0]):
776                 for i99 in range(a.shape[1]):
777                   if o=="1":
778                      out+=tap+"%s[%s,%s]=(%s)\n"%(arg,i0,i99,a[i0,i99]+b[i0,i99])
779                   else:
780                      out+=tap+"%s[%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i99,a[i0,i99],i99,b[i0,i99])
781    
782        elif a.rank==3:
783            for i0 in range(a.shape[0]):
784             for i1 in range(a.shape[1]):
785                 for i99 in range(a.shape[2]):
786                   if o=="1":
787                      out+=tap+"%s[%s,%s,%s]=(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99]+b[i0,i1,i99])
788                   else:
789                      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])
790    
791        elif a.rank==4:
792            for i0 in range(a.shape[0]):
793             for i1 in range(a.shape[1]):
794               for i2 in range(a.shape[2]):
795                 for i99 in range(a.shape[3]):
796                   if o=="1":
797                     out+=tap+"%s[%s,%s,%s,%s]=(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99]+b[i0,i1,i2,i99])
798                   else:
799                     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])
800        return out
801    def unrollLoopsOfDiv(a,b,o,arg,tap=""):
802        out=tap+arg+"="
803        if o=="1":
804           z=0.
805           for i99 in range(a.shape[0]):
806                z+=b[i99,i99]+a[i99,i99]
807           out+="(%s)"%z    
808        else:
809           z=0.
810           for i99 in range(a.shape[0]):
811                z+=b[i99,i99]
812                if i99>0: out+="+"
813                out+="o*(%s)*x_ref[%s]**(o-1)"%(a[i99,i99],i99)
814           out+="+(%s)"%z    
815        return out
816    
817    def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
818        if where=="Function":
819           xfac_o=1.
820           xfac_op=0.
821           z_fac=1./2.
822           z_fac_s=""
823           zo_fac_s="/(o+1.)"
824        elif where=="FunctionOnBoundary":
825           xfac_o=1.
826           xfac_op=0.
827           z_fac=1.
828           z_fac_s="*dim"
829           zo_fac_s="*(1+2.*(dim-1.)/(o+1.))"
830        elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
831           xfac_o=0.
832           xfac_op=1.
833           z_fac=1./2.
834           z_fac_s=""
835           zo_fac_s="/(o+1.)"
836        out=""
837        if a.rank==1:
838                 zo=0.
839                 zop=0.
840                 z=0.
841                 for i99 in range(a.shape[0]):
842                      if i99==0:
843                        zo+=       xfac_o*a[i99]
844                        zop+=       xfac_op*a[i99]
845                      else:
846                        zo+=a[i99]
847                      z+=b[i99]
848    
849                 out+=tap+"%s=(%s)%s+(%s)%s"%(arg,zo,zo_fac_s,z*z_fac,z_fac_s)
850                 if zop==0.:
851                   out+="\n"
852                 else:
853                   out+="+(%s)*0.5**o\n"%zop
854        elif a.rank==2:
855            for i0 in range(a.shape[0]):
856                 zo=0.
857                 zop=0.
858                 z=0.
859                 for i99 in range(a.shape[1]):
860                      if i99==0:
861                        zo+=       xfac_o*a[i0,i99]
862                        zop+=       xfac_op*a[i0,i99]
863                      else:
864                        zo+=a[i0,i99]
865                      z+=b[i0,i99]
866    
867                 out+=tap+"%s[%s]=(%s)%s+(%s)%s"%(arg,i0,zo,zo_fac_s,z*z_fac,z_fac_s)
868                 if zop==0.:
869                   out+="\n"
870                 else:
871                   out+="+(%s)*0.5**o\n"%zop
872        elif a.rank==3:
873            for i0 in range(a.shape[0]):
874             for i1 in range(a.shape[1]):
875                 zo=0.
876                 zop=0.
877                 z=0.
878                 for i99 in range(a.shape[2]):
879                      if i99==0:
880                        zo+=       xfac_o*a[i0,i1,i99]
881                        zop+=       xfac_op*a[i0,i1,i99]
882                      else:
883                        zo+=a[i0,i1,i99]
884                      z+=b[i0,i1,i99]
885    
886                 out+=tap+"%s[%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,zo,zo_fac_s,z*z_fac,z_fac_s)
887                 if zop==0.:
888                   out+="\n"
889                 else:
890                   out+="+(%s)*0.5**o\n"%zop
891        elif a.rank==4:
892            for i0 in range(a.shape[0]):
893             for i1 in range(a.shape[1]):
894               for i2 in range(a.shape[2]):
895                 zo=0.
896                 zop=0.
897                 z=0.
898                 for i99 in range(a.shape[3]):
899                      if i99==0:
900                        zo+=       xfac_o*a[i0,i1,i2,i99]
901                        zop+=       xfac_op*a[i0,i1,i2,i99]
902    
903                      else:
904                        zo+=a[i0,i1,i2,i99]
905                      z+=b[i0,i1,i2,i99]
906    
907                 out+=tap+"%s[%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,zo,zo_fac_s,z*z_fac,z_fac_s)
908                 if zop==0.:
909                   out+="\n"
910                 else:
911                   out+="+(%s)*0.5**o\n"%zop
912    
913        elif a.rank==5:
914            for i0 in range(a.shape[0]):
915             for i1 in range(a.shape[1]):
916               for i2 in range(a.shape[2]):
917                for i3 in range(a.shape[3]):
918                 zo=0.
919                 zop=0.
920                 z=0.
921                 for i99 in range(a.shape[4]):
922                      if i99==0:
923                        zo+=       xfac_o*a[i0,i1,i2,i3,i99]
924                        zop+=       xfac_op*a[i0,i1,i2,i3,i99]
925    
926                      else:
927                        zo+=a[i0,i1,i2,i3,i99]
928                      z+=b[i0,i1,i2,i3,i99]
929                 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)
930                 if zop==0.:
931                   out+="\n"
932                 else:
933                   out+="+(%s)*0.5**o\n"%zop
934    
935        return out
936    def unrollLoopsSimplified(b,arg,tap=""):
937        out=""
938        if isinstance(b,float) or b.rank==0:
939                 out+=tap+"%s=(%s)*x[0]\n"%(arg,str(b))
940    
941        elif b.rank==1:
942            for i0 in range(b.shape[0]):
943                 out+=tap+"%s[%s]=(%s)*x[%s]\n"%(arg,i0,b[i0],i0)
944        elif b.rank==2:
945            for i0 in range(b.shape[0]):
946             for i1 in range(b.shape[1]):
947                 out+=tap+"%s[%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,b[i0,i1],i1)
948        elif b.rank==3:
949            for i0 in range(b.shape[0]):
950             for i1 in range(b.shape[1]):
951               for i2 in range(b.shape[2]):
952                 out+=tap+"%s[%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,b[i0,i1,i2],i2)
953        elif b.rank==4:
954            for i0 in range(b.shape[0]):
955             for i1 in range(b.shape[1]):
956               for i2 in range(b.shape[2]):
957                for i3 in range(b.shape[3]):
958                 out+=tap+"%s[%s,%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,i3,b[i0,i1,i2,i3],i3)
959        return out
960    
961    def unrollLoopsOfL2(b,where,arg,tap=""):
962        out=""
963        z=[]
964        if isinstance(b,float) or b.rank==0:
965           z.append(b**2)
966        elif b.rank==1:
967            for i0 in range(b.shape[0]):
968                 z.append(b[i0]**2)
969        elif b.rank==2:
970            for i1 in range(b.shape[1]):
971               s=0
972               for i0 in range(b.shape[0]):
973                  s+=b[i0,i1]**2
974               z.append(s)
975        elif b.rank==3:
976            for i2 in range(b.shape[2]):
977              s=0
978              for i0 in range(b.shape[0]):
979                 for i1 in range(b.shape[1]):
980                    s+=b[i0,i1,i2]**2
981              z.append(s)
982    
983        elif b.rank==4:
984          for i3 in range(b.shape[3]):
985             s=0
986             for i0 in range(b.shape[0]):
987               for i1 in range(b.shape[1]):
988                  for i2 in range(b.shape[2]):
989                     s+=b[i0,i1,i2,i3]**2
990             z.append(s)        
991        if where=="Function":
992           xfac_o=1.
993           xfac_op=0.
994           z_fac_s=""
995           zo_fac_s=""
996           zo_fac=1./3.
997        elif where=="FunctionOnBoundary":
998           xfac_o=1.
999           xfac_op=0.
1000           z_fac_s="*dim"
1001           zo_fac_s="*(2.*dim+1.)/3."
1002           zo_fac=1.
1003        elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1004           xfac_o=0.
1005           xfac_op=1.
1006           z_fac_s=""
1007           zo_fac_s=""    
1008           zo_fac=1./3.    
1009        zo=0.
1010        zop=0.
1011        for i99 in range(len(z)):
1012               if i99==0:
1013                   zo+=xfac_o*z[i99]
1014                   zop+=xfac_op*z[i99]
1015               else:
1016                   zo+=z[i99]
1017        out+=tap+"%s=sqrt((%s)%s"%(arg,zo*zo_fac,zo_fac_s)
1018        if zop==0.:
1019           out+=")\n"
1020        else:
1021           out+="+(%s))\n"%(zop*0.5**2)
1022        return out
1023    
1024  #=======================================================================================================  #=======================================================================================================
1025  # local reduction  # L2
1026  #=======================================================================================================  #=======================================================================================================
1027  for oper in [["length",0.,"out+%a1%**2","math.sqrt(out)"],  for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1028               ["maxval",-1.e99,"max(out,%a1%)","out"],    for data in ["Data","Symbol"]:
1029               ["minval",1.e99,"min(out,%a1%)","out"] ]:      for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1030    for case in case_set:           text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1031       for sh in shape_set:           tname="test_L2_on%s_from%s_rank%s"%(where,data,len(sh))
1032         if not case=="float" or len(sh)==0:           text+="   def %s(self):\n"%tname
1033           text=""           text+="      \"\"\"\n"
1034           text+="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"           text+="      tests L2-norm of %s on the %s\n\n"%(data,where)
1035           tname="def test_%s_%s_rank%s"%(oper[0],case,len(sh))           text+="      assumptions: self.domain supports integration on %s\n"%where
1036           text+="   %s(self):\n"%tname           text+="      \"\"\"\n"
1037           a=makeArray(sh,[-1.,1.])                       text+="      dim=self.domain.getDim()\n"
1038           a1=makeArray(sh,[-1.,1.])           text+="      w=%s(self.domain)\n"%where
1039           r1=testReduce(a1,oper[1],oper[2],oper[3])           text+="      x=w.getX()\n"
1040           r=testReduce(a,oper[1],oper[2],oper[3])           o="1"
1041             if len(sh)>0:
1042                sh_2=sh[:len(sh)-1]+(2,)
1043                sh_3=sh[:len(sh)-1]+(3,)            
1044                b_2=makeArray(sh[:len(sh)-1]+(2,),[-1.,1])
1045                b_3=makeArray(sh[:len(sh)-1]+(3,),[-1.,1])
1046             else:
1047                sh_2=()
1048                sh_3=()
1049                b_2=makeArray(sh,[-1.,1])
1050                b_3=makeArray(sh,[-1.,1])
1051    
1052             if data=="Symbol":
1053                val="s"
1054                res="sub"
1055             else:
1056                val="arg"
1057                res="res"
1058             text+="      if dim==2:\n"
1059             if data=="Symbol":
1060                   text+="        arg=Symbol(shape=%s,dim=dim)\n"%str(sh_2)
1061    
1062             text+="        %s=Data(0,%s,w)\n"%(val,sh_2)
1063             text+=unrollLoopsSimplified(b_2,val,tap="        ")
1064             text+=unrollLoopsOfL2(b_2,where,"ref",tap="        ")
1065             text+="\n      else:\n"
1066             if data=="Symbol":
1067                   text+="        arg=Symbol(shape=%s,dim=dim)\n"%str(sh_3)
1068             text+="        %s=Data(0,%s,w)\n"%(val,sh_3)        
1069             text+=unrollLoopsSimplified(b_3,val,tap="        ")
1070             text+=unrollLoopsOfL2(b_3,where,"ref",tap="        ")
1071             text+="\n      res=L2(arg)\n"
1072             if data=="Symbol":
1073                text+="      sub=res.substitute({arg:s})\n"
1074                text+="      self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1075                text+="      self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1076             else:
1077                text+="      self.failUnless(isinstance(res,float),\"wrong type of result.\")\n"
1078             text+="      self.failUnlessAlmostEqual(%s,ref,int(-log10(self.RES_TOL)),\"wrong result\")\n"%res
1079             t_prog+=text
1080    print t_prog
1081    1/0
1082    
1083    #=======================================================================================================
1084    # div
1085    #=======================================================================================================
1086    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1087      for data in ["Data","Symbol"]:
1088         for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1089             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1090             tname="test_div_on%s_from%s_%s"%(where,data,case)
1091             text+="   def %s(self):\n"%tname
1092             text+="      \"\"\"\n"
1093             text+="      tests divergence of %s on the %s\n\n"%(data,where)
1094             text+="      assumptions: %s(self.domain) exists\n"%case
1095             text+="                   self.domain supports div on %s\n"%where
1096             text+="      \"\"\"\n"
1097             if case=="ReducedSolution":
1098                text+="      o=1\n"
1099                o="1"
1100             else:
1101                text+="      o=self.order\n"
1102                o="o"
1103             text+="      dim=self.domain.getDim()\n"
1104             text+="      w_ref=%s(self.domain)\n"%where
1105             text+="      x_ref=w_ref.getX()\n"
1106             text+="      w=%s(self.domain)\n"%case
1107             text+="      x=w.getX()\n"
1108             a_2=makeArray((2,2),[-1.,1])
1109             b_2=makeArray((2,2),[-1.,1])
1110             a_3=makeArray((3,3),[-1.,1])
1111             b_3=makeArray((3,3),[-1.,1])
1112             if data=="Symbol":
1113                text+="      arg=Symbol(shape=(dim,),dim=dim)\n"
1114                val="s"
1115                res="sub"
1116             else:
1117                val="arg"
1118                res="res"
1119             text+="      %s=Vector(0,w)\n"%val
1120             text+="      if dim==2:\n"
1121             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1122             text+=unrollLoopsOfDiv(a_2,b_2,o,"ref",tap="        ")
1123             text+="\n      else:\n"
1124                    
1125           text+=mkText(case,"arg",a,a1)           text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1126           text+="      res=%s(arg)\n"%oper[0]           text+=unrollLoopsOfDiv(a_3,b_3,o,"ref",tap="        ")
1127           if case=="Symbol":                   text+="\n      res=div(arg,where=w_ref)\n"
1128               text+=mkText("array","s",a,a1)           if data=="Symbol":
1129               text+="      sub=res.substitute({arg:s})\n"                      text+="      sub=res.substitute({arg:s})\n"
1130               text+=mkText("array","ref",r,r1)           text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1131               res="sub"           text+="      self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1132             text+="      self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong function space of result.\")\n"%res
1133             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1134    
1135    
1136             t_prog+=text
1137    print t_prog
1138    1/0
1139    
1140    #=======================================================================================================
1141    # interpolation
1142    #=======================================================================================================
1143    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne","Solution","ReducedSolution"]:
1144      for data in ["Data","Symbol"]:
1145         for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1146          for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1147            if  where==case or \
1148                ( case in ["ContinuousFunction","Solution","ReducedSolution"] and where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"] ) or \
1149                ( case in ["FunctionOnContactZero","FunctionOnContactOne"] and where in ["FunctionOnContactZero","FunctionOnContactOne"] ) or \
1150                (case=="ContinuousFunction" and  where in ["Solution","ReducedSolution"]) or \
1151                (case=="Solution" and  where=="ReducedSolution") :
1152                
1153    
1154             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1155             tname="test_interpolation_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1156             text+="   def %s(self):\n"%tname
1157             text+="      \"\"\"\n"
1158             text+="      tests interpolation for rank %s %s onto the %s\n\n"%(len(sh),data,where)
1159             text+="      assumptions: self.domain supports inpterpolation from %s onto %s\n"%(case,where)
1160             text+="      \"\"\"\n"
1161             if case=="ReducedSolution" or where=="ReducedSolution":
1162                text+="      o=1\n"
1163                o="1"
1164             else:
1165                text+="      o=self.order\n"
1166                o="o"
1167             text+="      dim=self.domain.getDim()\n"
1168             text+="      w_ref=%s(self.domain)\n"%where
1169             text+="      x_ref=w_ref.getX()\n"
1170             text+="      w=%s(self.domain)\n"%case
1171             text+="      x=w.getX()\n"
1172             a_2=makeArray(sh+(2,),[-1.,1])
1173             b_2=makeArray(sh+(2,),[-1.,1])
1174             a_3=makeArray(sh+(3,),[-1.,1])
1175             b_3=makeArray(sh+(3,),[-1.,1])
1176             if data=="Symbol":
1177                text+="      arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1178                val="s"
1179                res="sub"
1180           else:           else:
1181               text+=mkText(case,"ref",r,r1)              val="arg"
1182               res="res"              res="res"
1183           if oper[0]=="length":                         text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
1184                 text+=mkTypeAndShapeTest(case,(),"res")           text+="      ref=Data(0,%s,w_ref)\n"%str(sh)
1185           else:                       text+="      if dim==2:\n"
1186              if case=="float" or case=="array":                   text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1187                 text+=mkTypeAndShapeTest("float",(),"res")           text+=unrollLoops(a_2,b_2,o,"ref",tap="        ",x="x_ref")
1188              else:                     text+="      else:\n"
1189                 text+=mkTypeAndShapeTest(case,(),"res")          
1190             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1191             text+=unrollLoops(a_3,b_3,o,"ref",tap="        ",x="x_ref")
1192             text+="      res=interpolate(arg,where=w_ref)\n"
1193             if data=="Symbol":
1194                text+="      sub=res.substitute({arg:s})\n"
1195             text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1196             text+="      self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong functionspace of result.\")\n"%res
1197             text+="      self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1198           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
1199           if case == "taggedData":           t_prog+=text
1200             t_prog_with_tags+=text  print test_header
1201    print t_prog
1202    print test_tail          
1203    1/0
1204    
1205    #=======================================================================================================
1206    # grad
1207    #=======================================================================================================
1208    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1209      for data in ["Data","Symbol"]:
1210         for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1211           for sh in [ (),(2,), (4,5), (6,2,2)]:
1212             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1213             tname="test_grad_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1214             text+="   def %s(self):\n"%tname
1215             text+="      \"\"\"\n"
1216             text+="      tests gradient for rank %s %s on the %s\n\n"%(len(sh),data,where)
1217             text+="      assumptions: %s(self.domain) exists\n"%case
1218             text+="                   self.domain supports gardient on %s\n"%where
1219             text+="      \"\"\"\n"
1220             if case=="ReducedSolution":
1221                text+="      o=1\n"
1222                o="1"
1223             else:
1224                text+="      o=self.order\n"
1225                o="o"
1226             text+="      dim=self.domain.getDim()\n"
1227             text+="      w_ref=%s(self.domain)\n"%where
1228             text+="      x_ref=w_ref.getX()\n"
1229             text+="      w=%s(self.domain)\n"%case
1230             text+="      x=w.getX()\n"
1231             a_2=makeArray(sh+(2,),[-1.,1])
1232             b_2=makeArray(sh+(2,),[-1.,1])
1233             a_3=makeArray(sh+(3,),[-1.,1])
1234             b_3=makeArray(sh+(3,),[-1.,1])
1235             if data=="Symbol":
1236                text+="      arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1237                val="s"
1238                res="sub"
1239           else:           else:
1240             t_prog+=text              val="arg"
1241                res="res"
1242             text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
1243             text+="      ref=Data(0,%s+(dim,),w_ref)\n"%str(sh)
1244             text+="      if dim==2:\n"
1245             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1246             text+=unrollLoopsOfGrad(a_2,b_2,o,"ref",tap="        ")
1247             text+="      else:\n"
1248            
1249             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1250             text+=unrollLoopsOfGrad(a_3,b_3,o,"ref",tap="        ")
1251             text+="      res=grad(arg,where=w_ref)\n"
1252             if data=="Symbol":
1253                text+="      sub=res.substitute({arg:s})\n"
1254             text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1255             text+="      self.failUnlessEqual(res.getShape(),%s+(dim,),\"wrong shape of result.\")\n"%str(sh)
1256             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1257    
1258    
1259             t_prog+=text
1260    print test_header
1261    print t_prog
1262    print test_tail          
1263    1/0
1264    
1265    
1266    #=======================================================================================================
1267    # integrate
1268    #=======================================================================================================
1269    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1270      for data in ["Data","Symbol"]:
1271        for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1272          for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1273            if (not case in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]) or where==case:  
1274             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1275             tname="test_integrate_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1276             text+="   def %s(self):\n"%tname
1277             text+="      \"\"\"\n"
1278             text+="      tests integral of rank %s %s on the %s\n\n"%(len(sh),data,where)
1279             text+="      assumptions: %s(self.domain) exists\n"%case
1280             text+="                   self.domain supports integral on %s\n"%where
1281    
1282             text+="      \"\"\"\n"
1283             if case=="ReducedSolution":
1284                text+="      o=1\n"
1285                o="1"
1286             else:
1287                text+="      o=self.order\n"
1288                o="o"
1289             text+="      dim=self.domain.getDim()\n"
1290             text+="      w_ref=%s(self.domain)\n"%where
1291             text+="      w=%s(self.domain)\n"%case
1292             text+="      x=w.getX()\n"
1293             a_2=makeArray(sh+(2,),[-1.,1])
1294             b_2=makeArray(sh+(2,),[-1.,1])
1295             a_3=makeArray(sh+(3,),[-1.,1])
1296             b_3=makeArray(sh+(3,),[-1.,1])
1297             if data=="Symbol":
1298                text+="      arg=Symbol(shape=%s)\n"%str(sh)
1299                val="s"
1300                res="sub"
1301             else:
1302                val="arg"
1303                res="res"
1304                
1305             text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
1306             if not len(sh)==0:
1307                text+="      ref=numarray.zeros(%s,numarray.Float)\n"%str(sh)
1308             text+="      if dim==2:\n"
1309             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1310             text+=unrollLoopsOfInteriorIntegral(a_2,b_2,where,"ref",tap="        ")
1311             text+="      else:\n"
1312            
1313             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1314             text+=unrollLoopsOfInteriorIntegral(a_3,b_3,where,"ref",tap="        ")
1315             if case in ["ContinuousFunction","Solution","ReducedSolution"]:
1316                 text+="      res=integrate(arg,where=w_ref)\n"
1317             else:
1318                 text+="      res=integrate(arg)\n"
1319    
1320             if data=="Symbol":
1321                text+="      sub=res.substitute({arg:s})\n"
1322             if len(sh)==0 and data=="Data":
1323                text+="      self.failUnless(isinstance(%s,float),\"wrong type of result.\")\n"%res
1324             else:
1325                if data=="Symbol":
1326                   text+="      self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1327                   text+="      self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1328                else:
1329                   text+="      self.failUnless(isinstance(res,numarray.NumArray),\"wrong type of result.\")\n"
1330                   text+="      self.failUnlessEqual(res.shape,%s,\"wrong shape of result.\")\n"%str(sh)
1331             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1332    
1333    
1334             t_prog+=text
1335    print test_header
1336    print t_prog
1337    print test_tail          
1338    1/0
1339    #=======================================================================================================
1340    # inverse
1341    #=======================================================================================================
1342    name="inverse"
1343    for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1344      for sh0 in [ (1,1), (2,2), (3,3)]:
1345                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1346                  tname="test_%s_%s_dim%s"%(name,case0,sh0[0])
1347                  text+="   def %s(self):\n"%tname
1348                  a_0=makeArray(sh0,[-1.,1])
1349                  for i in range(sh0[0]): a_0[i,i]+=2
1350                  if case0 in ["taggedData", "expandedData"]:
1351                      a1_0=makeArray(sh0,[-1.,1])
1352                      for i in range(sh0[0]): a1_0[i,i]+=3
1353                  else:
1354                      a1_0=a_0
1355                      
1356                  text+=mkText(case0,"arg",a_0,a1_0)
1357                  text+="      res=%s(arg)\n"%name
1358                  if case0=="Symbol":
1359                     text+=mkText("array","s",a_0,a1_0)
1360                     text+="      sub=res.substitute({arg:s})\n"
1361                     res="sub"
1362                     ref="s"
1363                  else:
1364                     ref="arg"
1365                     res="res"
1366                  text+=mkTypeAndShapeTest(case0,sh0,"res")
1367                  text+="      self.failUnless(Lsup(matrixmult(%s,%s)-kronecker(%s))<=self.RES_TOL,\"wrong result\")\n"%(res,ref,sh0[0])
1368                  
1369                  if case0 == "taggedData" :
1370                      t_prog_with_tags+=text
1371                  else:              
1372                      t_prog+=text
1373    
1374  print test_header  print test_header
1375  # print t_prog  # print t_prog
1376  print t_prog_with_tags  print t_prog_with_tags
1377  print test_tail            print test_tail          
1378  1/0  1/0
1379              
1380  #=======================================================================================================  #=======================================================================================================
1381  # tensor multiply  # trace
1382  #=======================================================================================================  #=======================================================================================================
1383  # oper=["generalTensorProduct",tensorProductTest]  def traceTest(r,offset):
1384  # oper=["matrixmult",testMatrixMult]      sh=r.shape
1385  oper=["tensormult",testTensorMult]      r1=1
1386        for i in range(offset): r1*=sh[i]
1387        r2=1
1388        for i in range(offset+2,len(sh)): r2*=sh[i]
1389        r_s=numarray.reshape(r,(r1,sh[offset],sh[offset],r2))
1390        s=numarray.zeros([r1,r2],numarray.Float)
1391        for i1 in range(r1):
1392            for i2 in range(r2):
1393                for j in range(sh[offset]): s[i1,i2]+=r_s[i1,j,j,i2]
1394        return s.resize(sh[:offset]+sh[offset+2:])
1395    name,tt="trace",traceTest
1396    for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1397      for sh0 in [ (4,5), (6,2,2),(3,2,3,4)]:
1398        for offset in range(len(sh0)-1):
1399                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1400                  tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1401                  text+="   def %s(self):\n"%tname
1402                  sh_t=list(sh0)
1403                  sh_t[offset+1]=sh_t[offset]
1404                  sh_t=tuple(sh_t)
1405                  sh_r=[]
1406                  for i in range(offset): sh_r.append(sh0[i])
1407                  for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])              
1408                  sh_r=tuple(sh_r)
1409                  a_0=makeArray(sh_t,[-1.,1])
1410                  if case0 in ["taggedData", "expandedData"]:
1411                      a1_0=makeArray(sh_t,[-1.,1])
1412                  else:
1413                      a1_0=a_0
1414                  r=tt(a_0,offset)
1415                  r1=tt(a1_0,offset)
1416                  text+=mkText(case0,"arg",a_0,a1_0)
1417                  text+="      res=%s(arg,%s)\n"%(name,offset)
1418                  if case0=="Symbol":
1419                     text+=mkText("array","s",a_0,a1_0)
1420                     text+="      sub=res.substitute({arg:s})\n"
1421                     res="sub"
1422                     text+=mkText("array","ref",r,r1)
1423                  else:
1424                     res="res"
1425                     text+=mkText(case0,"ref",r,r1)
1426                  text+=mkTypeAndShapeTest(case0,sh_r,"res")
1427                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1428                  
1429                  if case0 == "taggedData" :
1430                      t_prog_with_tags+=text
1431                  else:              
1432                      t_prog+=text
1433    
1434  for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:  print test_header
1435    # print t_prog
1436    print t_prog_with_tags
1437    print test_tail          
1438    1/0
1439    
1440    #=======================================================================================================
1441    # clip
1442    #=======================================================================================================
1443    oper_L=[["clip",clipTEST]]
1444    for oper in oper_L:
1445     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1446    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)]:
1447            if len(sh0)==0 or not case0=="float":
1448                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1449                  tname="test_%s_%s_rank%s"%(oper[0],case0,len(sh0))
1450                  text+="   def %s(self):\n"%tname
1451                  a_0=makeArray(sh0,[-1.,1])
1452                  if case0 in ["taggedData", "expandedData"]:
1453                      a1_0=makeArray(sh0,[-1.,1])
1454                  else:
1455                      a1_0=a_0
1456    
1457                  r=oper[1](a_0,-0.3,0.5)
1458                  r1=oper[1](a1_0,-0.3,0.5)
1459                  text+=mkText(case0,"arg",a_0,a1_0)
1460                  text+="      res=%s(arg,-0.3,0.5)\n"%oper[0]
1461                  if case0=="Symbol":
1462                     text+=mkText("array","s",a_0,a1_0)
1463                     text+="      sub=res.substitute({arg:s})\n"
1464                     res="sub"
1465                     text+=mkText("array","ref",r,r1)
1466                  else:
1467                     res="res"
1468                     text+=mkText(case0,"ref",r,r1)
1469                  text+=mkTypeAndShapeTest(case0,sh0,"res")
1470                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1471                  
1472                  if case0 == "taggedData" :
1473                      t_prog_with_tags+=text
1474                  else:              
1475                      t_prog+=text
1476    
1477    print test_header
1478    # print t_prog
1479    print t_prog_with_tags
1480    print test_tail          
1481    1/0
1482    
1483    #=======================================================================================================
1484    # maximum, minimum, clipping
1485    #=======================================================================================================
1486    oper_L=[ ["maximum",maximumTEST],
1487             ["minimum",minimumTEST]]
1488    for oper in oper_L:
1489     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1490      for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1491     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1492       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)]:
1493         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") \
1494            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)    
1495                use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"                use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1496    
1497                text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"                text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1498                # 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  
1499                text+="   def %s(self):\n"%tname                text+="   def %s(self):\n"%tname
1500                a_0=makeArray(sh0+sh_s,[-1.,1])                a_0=makeArray(sh0,[-1.,1])
1501                if case0 in ["taggedData", "expandedData"]:                if case0 in ["taggedData", "expandedData"]:
1502                    a1_0=makeArray(sh0+sh_s,[-1.,1])                    a1_0=makeArray(sh0,[-1.,1])
1503                else:                else:
1504                    a1_0=a_0                    a1_0=a_0
1505    
1506                a_1=makeArray(sh_s+sh1,[-1.,1])                a_1=makeArray(sh1,[-1.,1])
1507                if case1 in ["taggedData", "expandedData"]:                if case1 in ["taggedData", "expandedData"]:
1508                    a1_1=makeArray(sh_s+sh1,[-1.,1])                    a1_1=makeArray(sh1,[-1.,1])
1509                else:                else:
1510                    a1_1=a_1                    a1_1=a_1
1511                r=oper[1](a_0,a_1,sh_s)                r=oper[1](a_0,a_1)
1512                r1=oper[1](a1_0,a1_1,sh_s)                r1=oper[1](a1_0,a1_1)
1513                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)
1514                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)
1515                #text+="      res=matrixmult(arg0,arg1)\n"                text+="      res=%s(arg0,arg1)\n"%oper[0]
1516                text+="      res=tensormult(arg0,arg1)\n"                case=getResultCaseForBin(case0,case1)              
               #text+="      res=generalTensorProduct(arg0,arg1,offset=%s)\n"%(len(sh_s))  
1517                if case=="Symbol":                if case=="Symbol":
1518                   c0_res,c1_res=case0,case1                   c0_res,c1_res=case0,case1
1519                   subs="{"                   subs="{"
# Line 724  for case0 in ["float","array","Symbol"," Line 1533  for case0 in ["float","array","Symbol","
1533                else:                else:
1534                   res="res"                   res="res"
1535                   text+=mkText(case,"ref",r,r1)                   text+=mkText(case,"ref",r,r1)
1536                text+=mkTypeAndShapeTest(case,sh0+sh1,"res")                if len(sh0)>len(sh1):
1537                      text+=mkTypeAndShapeTest(case,sh0,"res")
1538                  else:
1539                      text+=mkTypeAndShapeTest(case,sh1,"res")
1540                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
1541                  
1542                if case0 == "taggedData" or case1 == "taggedData":                if case0 == "taggedData" or case1 == "taggedData":
1543                    t_prog_with_tags+=text                    t_prog_with_tags+=text
1544                else:                              else:              
1545                    t_prog+=text                    t_prog+=text
1546    
1547  print test_header  print test_header
1548  # print t_prog  # print t_prog
1549  print t_prog_with_tags  print t_prog_with_tags
1550  print test_tail            print test_tail          
1551  1/0  1/0
1552    
1553    
1554  #=======================================================================================================  #=======================================================================================================
1555  # outer/inner  # outer inner
1556  #=======================================================================================================  #=======================================================================================================
1557  oper=["inner",innerTEST]  oper=["outer",outerTEST]
1558  # oper=["outer",outerTEST]  # oper=["inner",innerTEST]
1559  for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:  for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1560    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)]:
1561     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
# Line 798  for case0 in ["float","array","Symbol"," Line 1614  for case0 in ["float","array","Symbol","
1614  print test_header  print test_header
1615  # print t_prog  # print t_prog
1616  print t_prog_with_tags  print t_prog_with_tags
1617    print test_tail          
1618    1/0
1619    
1620    #=======================================================================================================
1621    # local reduction
1622    #=======================================================================================================
1623    for oper in [["length",0.,"out+%a1%**2","math.sqrt(out)"],
1624                 ["maxval",-1.e99,"max(out,%a1%)","out"],
1625                 ["minval",1.e99,"min(out,%a1%)","out"] ]:
1626      for case in case_set:
1627         for sh in shape_set:
1628           if not case=="float" or len(sh)==0:
1629             text=""
1630             text+="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1631             tname="def test_%s_%s_rank%s"%(oper[0],case,len(sh))
1632             text+="   %s(self):\n"%tname
1633             a=makeArray(sh,[-1.,1.])            
1634             a1=makeArray(sh,[-1.,1.])
1635             r1=testReduce(a1,oper[1],oper[2],oper[3])
1636             r=testReduce(a,oper[1],oper[2],oper[3])
1637            
1638             text+=mkText(case,"arg",a,a1)
1639             text+="      res=%s(arg)\n"%oper[0]
1640             if case=="Symbol":        
1641                 text+=mkText("array","s",a,a1)
1642                 text+="      sub=res.substitute({arg:s})\n"        
1643                 text+=mkText("array","ref",r,r1)
1644                 res="sub"
1645             else:
1646                 text+=mkText(case,"ref",r,r1)
1647                 res="res"
1648             if oper[0]=="length":              
1649                   text+=mkTypeAndShapeTest(case,(),"res")
1650             else:            
1651                if case=="float" or case=="array":        
1652                   text+=mkTypeAndShapeTest("float",(),"res")
1653                else:          
1654                   text+=mkTypeAndShapeTest(case,(),"res")
1655             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1656             if case == "taggedData":
1657               t_prog_with_tags+=text
1658             else:
1659               t_prog+=text
1660    print test_header
1661    # print t_prog
1662    print t_prog_with_tags
1663    print test_tail          
1664    1/0
1665              
1666    #=======================================================================================================
1667    # tensor multiply
1668    #=======================================================================================================
1669    # oper=["generalTensorProduct",tensorProductTest]
1670    # oper=["matrixmult",testMatrixMult]
1671    oper=["tensormult",testTensorMult]
1672    
1673    for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1674      for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1675       for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1676         for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1677           for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:
1678              if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \
1679                   and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
1680                # 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
1681                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
1682                  case=getResultCaseForBin(case0,case1)  
1683                  use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1684                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1685                  # tname="test_generalTensorProduct_%s_rank%s_%s_rank%s_offset%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1),len(sh_s))
1686                  #tname="test_matrixmult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
1687                  tname="test_tensormult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
1688                  # if tname=="test_generalTensorProduct_array_rank1_array_rank2_offset1":
1689                  # print tnametest_generalTensorProduct_Symbol_rank1_Symbol_rank3_offset1
1690                  text+="   def %s(self):\n"%tname
1691                  a_0=makeArray(sh0+sh_s,[-1.,1])
1692                  if case0 in ["taggedData", "expandedData"]:
1693                      a1_0=makeArray(sh0+sh_s,[-1.,1])
1694                  else:
1695                      a1_0=a_0
1696    
1697                  a_1=makeArray(sh_s+sh1,[-1.,1])
1698                  if case1 in ["taggedData", "expandedData"]:
1699                      a1_1=makeArray(sh_s+sh1,[-1.,1])
1700                  else:
1701                      a1_1=a_1
1702                  r=oper[1](a_0,a_1,sh_s)
1703                  r1=oper[1](a1_0,a1_1,sh_s)
1704                  text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1705                  text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1706                  #text+="      res=matrixmult(arg0,arg1)\n"
1707                  text+="      res=tensormult(arg0,arg1)\n"
1708                  #text+="      res=generalTensorProduct(arg0,arg1,offset=%s)\n"%(len(sh_s))
1709                  if case=="Symbol":
1710                     c0_res,c1_res=case0,case1
1711                     subs="{"
1712                     if case0=="Symbol":        
1713                        text+=mkText("array","s0",a_0,a1_0)
1714                        subs+="arg0:s0"
1715                        c0_res="array"
1716                     if case1=="Symbol":        
1717                        text+=mkText("array","s1",a_1,a1_1)
1718                        if not subs.endswith("{"): subs+=","
1719                        subs+="arg1:s1"
1720                        c1_res="array"
1721                     subs+="}"  
1722                     text+="      sub=res.substitute(%s)\n"%subs
1723                     res="sub"
1724                     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1725                  else:
1726                     res="res"
1727                     text+=mkText(case,"ref",r,r1)
1728                  text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
1729                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1730                  if case0 == "taggedData" or case1 == "taggedData":
1731                      t_prog_with_tags+=text
1732                  else:              
1733                      t_prog+=text
1734    print test_header
1735    # print t_prog
1736    print t_prog_with_tags
1737  print test_tail            print test_tail          
1738  1/0  1/0
1739  #=======================================================================================================  #=======================================================================================================

Legend:
Removed from v.313  
changed lines
  Added in v.443

  ViewVC Help
Powered by ViewVC 1.1.26