/[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 517 by gross, Tue Feb 14 02:25:02 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    def makeNumberedArray(shape,s=1.):
188       out=numarray.zeros(shape,numarray.Float64)
189       if len(shape)==0:
190           out=s*1.
191       elif len(shape)==1:
192           for i0 in range(shape[0]):
193                       out[i0]=s*i0
194       elif len(shape)==2:
195           for i0 in range(shape[0]):
196              for i1 in range(shape[1]):
197                       out[i0,i1]=s*(i1+shape[1]*i0)
198       elif len(shape)==3:
199           for i0 in range(shape[0]):
200              for i1 in range(shape[1]):
201                 for i2 in range(shape[2]):
202                       out[i0,i1,i2]=s*(i2+shape[2]*i1+shape[2]*shape[1]*i0)
203       elif len(shape)==4:
204           for i0 in range(shape[0]):
205              for i1 in range(shape[1]):
206                 for i2 in range(shape[2]):
207                    for i3 in range(shape[3]):
208                       out[i0,i1,i2,i3]=s*(i3+shape[3]*i2+shape[3]*shape[2]*i1+shape[3]*shape[2]*shape[1]*i0)
209       else:
210           raise SystemError,"rank is restricted to 4"
211       return out        
212    
213  def makeResult(val,test_expr):  def makeResult(val,test_expr):
214     if isinstance(val,float):     if isinstance(val,float):
# Line 518  def mkCode(txt,args=[],intend=""): Line 550  def mkCode(txt,args=[],intend=""):
550        out=out.replace("%%a%s%%"%c,r)        out=out.replace("%%a%s%%"%c,r)
551      return out        return out  
552    
553    #=======================================================================================================
554    # slicing
555    #=======================================================================================================
556    for case0 in ["constData","taggedData","expandedData","Symbol"]:
557      for sh0 in [ (3,), (3,4), (3,4,3) ,(4,3,5,3)]:
558        # get perm:
559        if len(sh0)==2:
560            check=[[1,0]]
561        elif len(sh0)==3:
562            check=[[1,0,2],
563                   [1,2,0],
564                   [2,1,0],
565                   [2,0,2],
566                   [0,2,1]]
567        elif len(sh0)==4:
568            check=[[0,1,3,2],
569                   [0,2,1,3],
570                   [0,2,3,1],
571                   [0,3,2,1],
572                   [0,3,1,2] ,          
573                   [1,0,2,3],
574                   [1,0,3,2],
575                   [1,2,0,3],
576                   [1,2,3,0],
577                   [1,3,2,0],
578                   [1,3,0,2],
579                   [2,0,1,3],
580                   [2,0,3,1],
581                   [2,1,0,3],
582                   [2,1,3,0],
583                   [2,3,1,0],
584                   [2,3,0,1],
585                   [3,0,1,2],
586                   [3,0,2,1],
587                   [3,1,0,2],
588                   [3,1,2,0],
589                   [3,2,1,0],
590                   [3,2,0,1]]
591        else:
592             check=[]
593        
594        # create the test cases:
595        processed=[]
596        l=["R","U","L","P","C","N"]
597        c=[""]
598        for i in range(len(sh0)):
599           tmp=[]
600           for ci in c:
601              tmp+=[ci+li for li in l]
602           c=tmp
603        # SHUFFLE
604        c2=[]
605        while len(c)>0:
606            i=int(random.random()*len(c))
607            c2.append(c[i])
608            del c[i]
609        c=c2
610        for ci in c:
611          t=""
612          sh=()
613          for i in range(len(ci)):
614              if ci[i]=="R":
615                 s="%s:%s"%(1,sh0[i]-1)
616                 sh=sh+(sh0[i]-2,)            
617              if ci[i]=="U":
618                  s=":%s"%(sh0[i]-1)
619                  sh=sh+(sh0[i]-1,)            
620              if ci[i]=="L":
621                  s="2:"
622                  sh=sh+(sh0[i]-2,)            
623              if ci[i]=="P":
624                  s="%s"%(int(sh0[i]/2))
625              if ci[i]=="C":
626                  s=":"
627                  sh=sh+(sh0[i],)            
628              if ci[i]=="N":
629                  s=""
630                  sh=sh+(sh0[i],)
631              if len(s)>0:
632                 if not t=="": t+=","
633                 t+=s
634          N_found=False
635          noN_found=False
636          process=len(t)>0
637          for i in ci:
638             if i=="N":
639                if not noN_found and N_found: process=False
640                N_found=True
641             else:
642               if N_found: process=False
643               noNfound=True
644          # is there a similar one processed allready
645          if process and ci.find("N")==-1:
646             for ci2 in processed:
647               for chi in check:
648                   is_perm=True
649                   for i in range(len(chi)):
650                       if not ci[i]==ci2[chi[i]]: is_perm=False
651                   if is_perm: process=False
652          # if not process: print ci," rejected"
653          if process:
654           processed.append(ci)
655           text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
656           tname="test_getslice_%s_rank%s_%s"%(case0,len(sh0),ci)
657           text+="   def %s(self):\n"%tname
658           a_0=makeNumberedArray(sh0,s=1)
659           if case0 in ["taggedData", "expandedData"]:
660                a1_0=makeNumberedArray(sh0,s=-1.)
661           else:
662                a1_0=a_0
663           r=eval("a_0[%s]"%t)
664           r1=eval("a1_0[%s]"%t)
665           text+=mkText(case0,"arg",a_0,a1_0)
666           text+="      res=arg[%s]\n"%t
667           if case0=="Symbol":
668               text+=mkText("array","s",a_0,a1_0)
669               text+="      sub=res.substitute({arg:s})\n"
670               res="sub"
671               text+=mkText("array","ref",r,r1)
672           else:
673               res="res"
674               text+=mkText(case0,"ref",r,r1)
675           text+=mkTypeAndShapeTest(case0,sh,"res")
676           text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
677                  
678           if case0 == "taggedData" :
679                t_prog_with_tags+=text
680           else:              
681                t_prog+=text
682    
683    print test_header
684    # print t_prog
685    print t_prog_with_tags
686    print test_tail          
687    1/0
688    #============================================================================================
689  def innerTEST(arg0,arg1):  def innerTEST(arg0,arg1):
690      if isinstance(arg0,float):      if isinstance(arg0,float):
691         out=numarray.array(arg0*arg1)         out=numarray.array(arg0*arg1)
# Line 616  def testReduce(arg0,init_val,test_expr,p Line 784  def testReduce(arg0,init_val,test_expr,p
784               for i3 in range(arg0.shape[3]):               for i3 in range(arg0.shape[3]):
785                 out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2,i3]"))                           out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2,i3]"))          
786       return eval(post_expr)       return eval(post_expr)
787        
788    def clipTEST(arg0,mn,mx):
789         if isinstance(arg0,float):
790              return max(min(arg0,mx),mn)
791         out=numarray.zeros(arg0.shape,numarray.Float64)
792         if arg0.rank==1:
793            for i0 in range(arg0.shape[0]):
794                out[i0]=max(min(arg0[i0],mx),mn)
795         elif arg0.rank==2:
796            for i0 in range(arg0.shape[0]):
797             for i1 in range(arg0.shape[1]):
798                out[i0,i1]=max(min(arg0[i0,i1],mx),mn)
799         elif arg0.rank==3:
800            for i0 in range(arg0.shape[0]):
801             for i1 in range(arg0.shape[1]):
802               for i2 in range(arg0.shape[2]):
803                  out[i0,i1,i2]=max(min(arg0[i0,i1,i2],mx),mn)
804         elif arg0.rank==4:
805            for i0 in range(arg0.shape[0]):
806             for i1 in range(arg0.shape[1]):
807               for i2 in range(arg0.shape[2]):
808                 for i3 in range(arg0.shape[3]):
809                    out[i0,i1,i2,i3]=max(min(arg0[i0,i1,i2,i3],mx),mn)
810         return out
811    def minimumTEST(arg0,arg1):
812         if isinstance(arg0,float):
813           if isinstance(arg1,float):
814              if arg0>arg1:
815                  return arg1
816              else:
817                  return arg0
818           else:
819              arg0=numarray.ones(arg1.shape)*arg0
820         else:
821           if isinstance(arg1,float):
822              arg1=numarray.ones(arg0.shape)*arg1
823         out=numarray.zeros(arg0.shape,numarray.Float64)
824         if arg0.rank==0:
825              if arg0>arg1:
826                  out=arg1
827              else:
828                  out=arg0
829         elif arg0.rank==1:
830            for i0 in range(arg0.shape[0]):
831              if arg0[i0]>arg1[i0]:
832                  out[i0]=arg1[i0]
833              else:
834                  out[i0]=arg0[i0]
835         elif arg0.rank==2:
836            for i0 in range(arg0.shape[0]):
837             for i1 in range(arg0.shape[1]):
838              if arg0[i0,i1]>arg1[i0,i1]:
839                  out[i0,i1]=arg1[i0,i1]
840              else:
841                  out[i0,i1]=arg0[i0,i1]
842         elif arg0.rank==3:
843            for i0 in range(arg0.shape[0]):
844             for i1 in range(arg0.shape[1]):
845               for i2 in range(arg0.shape[2]):
846                 if arg0[i0,i1,i2]>arg1[i0,i1,i2]:
847                  out[i0,i1,i2]=arg1[i0,i1,i2]
848                 else:
849                  out[i0,i1,i2]=arg0[i0,i1,i2]
850         elif arg0.rank==4:
851            for i0 in range(arg0.shape[0]):
852             for i1 in range(arg0.shape[1]):
853               for i2 in range(arg0.shape[2]):
854                 for i3 in range(arg0.shape[3]):
855                  if arg0[i0,i1,i2,i3]>arg1[i0,i1,i2,i3]:
856                   out[i0,i1,i2,i3]=arg1[i0,i1,i2,i3]
857                  else:
858                   out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
859         return out
860        
861    def unrollLoops(a,b,o,arg,tap="",x="x"):
862        out=""
863        if a.rank==1:
864                 z=""
865                 for i99 in range(a.shape[0]):
866                   if not z=="": z+="+"
867                   if o=="1":
868                      z+="(%s)*%s[%s]"%(a[i99]+b[i99],x,i99)
869                   else:
870                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i99],x,i99,b[i99],x,i99)
871    
872                 out+=tap+"%s=%s\n"%(arg,z)
873    
874        elif a.rank==2:
875            for i0 in range(a.shape[0]):
876                 z=""
877                 for i99 in range(a.shape[1]):
878                   if not z=="": z+="+"
879                   if o=="1":
880                      z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)
881                   else:
882                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i99],x,i99,b[i0,i99],x,i99)
883    
884                 out+=tap+"%s[%s]=%s\n"%(arg,i0,z)
885        elif a.rank==3:
886            for i0 in range(a.shape[0]):
887             for i1 in range(a.shape[1]):
888                 z=""
889                 for i99 in range(a.shape[2]):
890                   if not z=="": z+="+"
891                   if o=="1":
892                      z+="(%s)*%s[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],x,i99)
893                   else:
894                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i99],x,i99,b[i0,i1,i99],x,i99)
895    
896                 out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)
897        elif a.rank==4:
898            for i0 in range(a.shape[0]):
899             for i1 in range(a.shape[1]):
900               for i2 in range(a.shape[2]):
901                 z=""
902                 for i99 in range(a.shape[3]):
903                   if not z=="": z+="+"
904                   if o=="1":
905                      z+="(%s)*%s[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],x,i99)
906                   else:
907                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i99],x,i99,b[i0,i1,i2,i99],x,i99)
908    
909                 out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)
910        elif a.rank==5:
911            for i0 in range(a.shape[0]):
912             for i1 in range(a.shape[1]):
913               for i2 in range(a.shape[2]):
914                for i3 in range(a.shape[3]):
915                 z=""
916                 for i99 in range(a.shape[4]):
917                   if not z=="": z+="+"
918                   if o=="1":
919                      z+="(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],x,i99)
920                   else:
921                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99],x,i99,b[i0,i1,i2,i3,i99],x,i99)
922    
923                 out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
924        return out
925    
926    def unrollLoopsOfGrad(a,b,o,arg,tap=""):
927        out=""
928        if a.rank==1:
929                 for i99 in range(a.shape[0]):
930                   if o=="1":
931                      out+=tap+"%s[%s]=(%s)\n"%(arg,i99,a[i99]+b[i99])
932                   else:
933                      out+=tap+"%s[%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i99,a[i99],i99,b[i99])
934    
935        elif a.rank==2:
936            for i0 in range(a.shape[0]):
937                 for i99 in range(a.shape[1]):
938                   if o=="1":
939                      out+=tap+"%s[%s,%s]=(%s)\n"%(arg,i0,i99,a[i0,i99]+b[i0,i99])
940                   else:
941                      out+=tap+"%s[%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i99,a[i0,i99],i99,b[i0,i99])
942    
943        elif a.rank==3:
944            for i0 in range(a.shape[0]):
945             for i1 in range(a.shape[1]):
946                 for i99 in range(a.shape[2]):
947                   if o=="1":
948                      out+=tap+"%s[%s,%s,%s]=(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99]+b[i0,i1,i99])
949                   else:
950                      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])
951    
952        elif a.rank==4:
953            for i0 in range(a.shape[0]):
954             for i1 in range(a.shape[1]):
955               for i2 in range(a.shape[2]):
956                 for i99 in range(a.shape[3]):
957                   if o=="1":
958                     out+=tap+"%s[%s,%s,%s,%s]=(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99]+b[i0,i1,i2,i99])
959                   else:
960                     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])
961        return out
962    def unrollLoopsOfDiv(a,b,o,arg,tap=""):
963        out=tap+arg+"="
964        if o=="1":
965           z=0.
966           for i99 in range(a.shape[0]):
967                z+=b[i99,i99]+a[i99,i99]
968           out+="(%s)"%z    
969        else:
970           z=0.
971           for i99 in range(a.shape[0]):
972                z+=b[i99,i99]
973                if i99>0: out+="+"
974                out+="o*(%s)*x_ref[%s]**(o-1)"%(a[i99,i99],i99)
975           out+="+(%s)"%z    
976        return out
977    
978    def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
979        if where=="Function":
980           xfac_o=1.
981           xfac_op=0.
982           z_fac=1./2.
983           z_fac_s=""
984           zo_fac_s="/(o+1.)"
985        elif where=="FunctionOnBoundary":
986           xfac_o=1.
987           xfac_op=0.
988           z_fac=1.
989           z_fac_s="*dim"
990           zo_fac_s="*(1+2.*(dim-1.)/(o+1.))"
991        elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
992           xfac_o=0.
993           xfac_op=1.
994           z_fac=1./2.
995           z_fac_s=""
996           zo_fac_s="/(o+1.)"
997        out=""
998        if a.rank==1:
999                 zo=0.
1000                 zop=0.
1001                 z=0.
1002                 for i99 in range(a.shape[0]):
1003                      if i99==0:
1004                        zo+=       xfac_o*a[i99]
1005                        zop+=       xfac_op*a[i99]
1006                      else:
1007                        zo+=a[i99]
1008                      z+=b[i99]
1009    
1010                 out+=tap+"%s=(%s)%s+(%s)%s"%(arg,zo,zo_fac_s,z*z_fac,z_fac_s)
1011                 if zop==0.:
1012                   out+="\n"
1013                 else:
1014                   out+="+(%s)*0.5**o\n"%zop
1015        elif a.rank==2:
1016            for i0 in range(a.shape[0]):
1017                 zo=0.
1018                 zop=0.
1019                 z=0.
1020                 for i99 in range(a.shape[1]):
1021                      if i99==0:
1022                        zo+=       xfac_o*a[i0,i99]
1023                        zop+=       xfac_op*a[i0,i99]
1024                      else:
1025                        zo+=a[i0,i99]
1026                      z+=b[i0,i99]
1027    
1028                 out+=tap+"%s[%s]=(%s)%s+(%s)%s"%(arg,i0,zo,zo_fac_s,z*z_fac,z_fac_s)
1029                 if zop==0.:
1030                   out+="\n"
1031                 else:
1032                   out+="+(%s)*0.5**o\n"%zop
1033        elif a.rank==3:
1034            for i0 in range(a.shape[0]):
1035             for i1 in range(a.shape[1]):
1036                 zo=0.
1037                 zop=0.
1038                 z=0.
1039                 for i99 in range(a.shape[2]):
1040                      if i99==0:
1041                        zo+=       xfac_o*a[i0,i1,i99]
1042                        zop+=       xfac_op*a[i0,i1,i99]
1043                      else:
1044                        zo+=a[i0,i1,i99]
1045                      z+=b[i0,i1,i99]
1046    
1047                 out+=tap+"%s[%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,zo,zo_fac_s,z*z_fac,z_fac_s)
1048                 if zop==0.:
1049                   out+="\n"
1050                 else:
1051                   out+="+(%s)*0.5**o\n"%zop
1052        elif a.rank==4:
1053            for i0 in range(a.shape[0]):
1054             for i1 in range(a.shape[1]):
1055               for i2 in range(a.shape[2]):
1056                 zo=0.
1057                 zop=0.
1058                 z=0.
1059                 for i99 in range(a.shape[3]):
1060                      if i99==0:
1061                        zo+=       xfac_o*a[i0,i1,i2,i99]
1062                        zop+=       xfac_op*a[i0,i1,i2,i99]
1063    
1064                      else:
1065                        zo+=a[i0,i1,i2,i99]
1066                      z+=b[i0,i1,i2,i99]
1067    
1068                 out+=tap+"%s[%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,zo,zo_fac_s,z*z_fac,z_fac_s)
1069                 if zop==0.:
1070                   out+="\n"
1071                 else:
1072                   out+="+(%s)*0.5**o\n"%zop
1073    
1074        elif a.rank==5:
1075            for i0 in range(a.shape[0]):
1076             for i1 in range(a.shape[1]):
1077               for i2 in range(a.shape[2]):
1078                for i3 in range(a.shape[3]):
1079                 zo=0.
1080                 zop=0.
1081                 z=0.
1082                 for i99 in range(a.shape[4]):
1083                      if i99==0:
1084                        zo+=       xfac_o*a[i0,i1,i2,i3,i99]
1085                        zop+=       xfac_op*a[i0,i1,i2,i3,i99]
1086    
1087                      else:
1088                        zo+=a[i0,i1,i2,i3,i99]
1089                      z+=b[i0,i1,i2,i3,i99]
1090                 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)
1091                 if zop==0.:
1092                   out+="\n"
1093                 else:
1094                   out+="+(%s)*0.5**o\n"%zop
1095    
1096        return out
1097    def unrollLoopsSimplified(b,arg,tap=""):
1098        out=""
1099        if isinstance(b,float) or b.rank==0:
1100                 out+=tap+"%s=(%s)*x[0]\n"%(arg,str(b))
1101    
1102        elif b.rank==1:
1103            for i0 in range(b.shape[0]):
1104                 out+=tap+"%s[%s]=(%s)*x[%s]\n"%(arg,i0,b[i0],i0)
1105        elif b.rank==2:
1106            for i0 in range(b.shape[0]):
1107             for i1 in range(b.shape[1]):
1108                 out+=tap+"%s[%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,b[i0,i1],i1)
1109        elif b.rank==3:
1110            for i0 in range(b.shape[0]):
1111             for i1 in range(b.shape[1]):
1112               for i2 in range(b.shape[2]):
1113                 out+=tap+"%s[%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,b[i0,i1,i2],i2)
1114        elif b.rank==4:
1115            for i0 in range(b.shape[0]):
1116             for i1 in range(b.shape[1]):
1117               for i2 in range(b.shape[2]):
1118                for i3 in range(b.shape[3]):
1119                 out+=tap+"%s[%s,%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,i3,b[i0,i1,i2,i3],i3)
1120        return out
1121    
1122    def unrollLoopsOfL2(b,where,arg,tap=""):
1123        out=""
1124        z=[]
1125        if isinstance(b,float) or b.rank==0:
1126           z.append(b**2)
1127        elif b.rank==1:
1128            for i0 in range(b.shape[0]):
1129                 z.append(b[i0]**2)
1130        elif b.rank==2:
1131            for i1 in range(b.shape[1]):
1132               s=0
1133               for i0 in range(b.shape[0]):
1134                  s+=b[i0,i1]**2
1135               z.append(s)
1136        elif b.rank==3:
1137            for i2 in range(b.shape[2]):
1138              s=0
1139              for i0 in range(b.shape[0]):
1140                 for i1 in range(b.shape[1]):
1141                    s+=b[i0,i1,i2]**2
1142              z.append(s)
1143    
1144        elif b.rank==4:
1145          for i3 in range(b.shape[3]):
1146             s=0
1147             for i0 in range(b.shape[0]):
1148               for i1 in range(b.shape[1]):
1149                  for i2 in range(b.shape[2]):
1150                     s+=b[i0,i1,i2,i3]**2
1151             z.append(s)        
1152        if where=="Function":
1153           xfac_o=1.
1154           xfac_op=0.
1155           z_fac_s=""
1156           zo_fac_s=""
1157           zo_fac=1./3.
1158        elif where=="FunctionOnBoundary":
1159           xfac_o=1.
1160           xfac_op=0.
1161           z_fac_s="*dim"
1162           zo_fac_s="*(2.*dim+1.)/3."
1163           zo_fac=1.
1164        elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1165           xfac_o=0.
1166           xfac_op=1.
1167           z_fac_s=""
1168           zo_fac_s=""    
1169           zo_fac=1./3.    
1170        zo=0.
1171        zop=0.
1172        for i99 in range(len(z)):
1173               if i99==0:
1174                   zo+=xfac_o*z[i99]
1175                   zop+=xfac_op*z[i99]
1176               else:
1177                   zo+=z[i99]
1178        out+=tap+"%s=sqrt((%s)%s"%(arg,zo*zo_fac,zo_fac_s)
1179        if zop==0.:
1180           out+=")\n"
1181        else:
1182           out+="+(%s))\n"%(zop*0.5**2)
1183        return out
1184  #=======================================================================================================  #=======================================================================================================
1185  # local reduction  # transpose
1186  #=======================================================================================================  #=======================================================================================================
1187  for oper in [["length",0.,"out+%a1%**2","math.sqrt(out)"],  def transposeTest(r,offset):
1188               ["maxval",-1.e99,"max(out,%a1%)","out"],      if isinstance(r,float): return r
1189               ["minval",1.e99,"min(out,%a1%)","out"] ]:      s=r.shape
1190    for case in case_set:      s1=1
1191       for sh in shape_set:      for i in s[:offset]: s1*=i
1192         if not case=="float" or len(sh)==0:      s2=1
1193           text=""      for i in s[offset:]: s2*=i
1194           text+="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"      out=numarray.reshape(r,(s1,s2))
1195           tname="def test_%s_%s_rank%s"%(oper[0],case,len(sh))      out.transpose()
1196           text+="   %s(self):\n"%tname      return numarray.resize(out,s[offset:]+s[:offset])
1197           a=makeArray(sh,[-1.,1.])              
1198           a1=makeArray(sh,[-1.,1.])  name,tt="transpose",transposeTest
1199           r1=testReduce(a1,oper[1],oper[2],oper[3])  for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1200           r=testReduce(a,oper[1],oper[2],oper[3])    for sh0 in [ (), (3,), (4,5), (6,2,2),(3,2,3,4)]:
1201        for offset in range(len(sh0)+1):
1202                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1203                  tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1204                  text+="   def %s(self):\n"%tname
1205                  sh_t=sh0[offset:]+sh0[:offset]
1206    
1207    #              sh_t=list(sh0)
1208    #              sh_t[offset+1]=sh_t[offset]
1209    #              sh_t=tuple(sh_t)
1210    #              sh_r=[]
1211    #              for i in range(offset): sh_r.append(sh0[i])
1212    #              for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])              
1213    #              sh_r=tuple(sh_r)
1214    
1215                  a_0=makeArray(sh0,[-1.,1])
1216                  if case0 in ["taggedData", "expandedData"]:
1217                      a1_0=makeArray(sh0,[-1.,1])
1218                  else:
1219                      a1_0=a_0
1220                  r=tt(a_0,offset)
1221                  r1=tt(a1_0,offset)
1222                  text+=mkText(case0,"arg",a_0,a1_0)
1223                  text+="      res=%s(arg,%s)\n"%(name,offset)
1224                  if case0=="Symbol":
1225                     text+=mkText("array","s",a_0,a1_0)
1226                     text+="      sub=res.substitute({arg:s})\n"
1227                     res="sub"
1228                     text+=mkText("array","ref",r,r1)
1229                  else:
1230                     res="res"
1231                     text+=mkText(case0,"ref",r,r1)
1232                  text+=mkTypeAndShapeTest(case0,sh_t,"res")
1233                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1234                  
1235                  if case0 == "taggedData" :
1236                      t_prog_with_tags+=text
1237                  else:              
1238                      t_prog+=text
1239    
1240    print test_header
1241    # print t_prog
1242    print t_prog_with_tags
1243    print test_tail          
1244    1/0
1245    #=======================================================================================================
1246    # L2
1247    #=======================================================================================================
1248    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1249      for data in ["Data","Symbol"]:
1250        for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1251             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1252             tname="test_L2_on%s_from%s_rank%s"%(where,data,len(sh))
1253             text+="   def %s(self):\n"%tname
1254             text+="      \"\"\"\n"
1255             text+="      tests L2-norm of %s on the %s\n\n"%(data,where)
1256             text+="      assumptions: self.domain supports integration on %s\n"%where
1257             text+="      \"\"\"\n"
1258             text+="      dim=self.domain.getDim()\n"
1259             text+="      w=%s(self.domain)\n"%where
1260             text+="      x=w.getX()\n"
1261             o="1"
1262             if len(sh)>0:
1263                sh_2=sh[:len(sh)-1]+(2,)
1264                sh_3=sh[:len(sh)-1]+(3,)            
1265                b_2=makeArray(sh[:len(sh)-1]+(2,),[-1.,1])
1266                b_3=makeArray(sh[:len(sh)-1]+(3,),[-1.,1])
1267             else:
1268                sh_2=()
1269                sh_3=()
1270                b_2=makeArray(sh,[-1.,1])
1271                b_3=makeArray(sh,[-1.,1])
1272    
1273             if data=="Symbol":
1274                val="s"
1275                res="sub"
1276             else:
1277                val="arg"
1278                res="res"
1279             text+="      if dim==2:\n"
1280             if data=="Symbol":
1281                   text+="        arg=Symbol(shape=%s,dim=dim)\n"%str(sh_2)
1282    
1283             text+="        %s=Data(0,%s,w)\n"%(val,sh_2)
1284             text+=unrollLoopsSimplified(b_2,val,tap="        ")
1285             text+=unrollLoopsOfL2(b_2,where,"ref",tap="        ")
1286             text+="\n      else:\n"
1287             if data=="Symbol":
1288                   text+="        arg=Symbol(shape=%s,dim=dim)\n"%str(sh_3)
1289             text+="        %s=Data(0,%s,w)\n"%(val,sh_3)        
1290             text+=unrollLoopsSimplified(b_3,val,tap="        ")
1291             text+=unrollLoopsOfL2(b_3,where,"ref",tap="        ")
1292             text+="\n      res=L2(arg)\n"
1293             if data=="Symbol":
1294                text+="      sub=res.substitute({arg:s})\n"
1295                text+="      self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1296                text+="      self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1297             else:
1298                text+="      self.failUnless(isinstance(res,float),\"wrong type of result.\")\n"
1299             text+="      self.failUnlessAlmostEqual(%s,ref,int(-log10(self.RES_TOL)),\"wrong result\")\n"%res
1300             t_prog+=text
1301    print t_prog
1302    1/0
1303    
1304    #=======================================================================================================
1305    # div
1306    #=======================================================================================================
1307    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1308      for data in ["Data","Symbol"]:
1309         for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1310             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1311             tname="test_div_on%s_from%s_%s"%(where,data,case)
1312             text+="   def %s(self):\n"%tname
1313             text+="      \"\"\"\n"
1314             text+="      tests divergence of %s on the %s\n\n"%(data,where)
1315             text+="      assumptions: %s(self.domain) exists\n"%case
1316             text+="                   self.domain supports div on %s\n"%where
1317             text+="      \"\"\"\n"
1318             if case=="ReducedSolution":
1319                text+="      o=1\n"
1320                o="1"
1321             else:
1322                text+="      o=self.order\n"
1323                o="o"
1324             text+="      dim=self.domain.getDim()\n"
1325             text+="      w_ref=%s(self.domain)\n"%where
1326             text+="      x_ref=w_ref.getX()\n"
1327             text+="      w=%s(self.domain)\n"%case
1328             text+="      x=w.getX()\n"
1329             a_2=makeArray((2,2),[-1.,1])
1330             b_2=makeArray((2,2),[-1.,1])
1331             a_3=makeArray((3,3),[-1.,1])
1332             b_3=makeArray((3,3),[-1.,1])
1333             if data=="Symbol":
1334                text+="      arg=Symbol(shape=(dim,),dim=dim)\n"
1335                val="s"
1336                res="sub"
1337             else:
1338                val="arg"
1339                res="res"
1340             text+="      %s=Vector(0,w)\n"%val
1341             text+="      if dim==2:\n"
1342             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1343             text+=unrollLoopsOfDiv(a_2,b_2,o,"ref",tap="        ")
1344             text+="\n      else:\n"
1345                    
1346           text+=mkText(case,"arg",a,a1)           text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1347           text+="      res=%s(arg)\n"%oper[0]           text+=unrollLoopsOfDiv(a_3,b_3,o,"ref",tap="        ")
1348           if case=="Symbol":                   text+="\n      res=div(arg,where=w_ref)\n"
1349               text+=mkText("array","s",a,a1)           if data=="Symbol":
1350               text+="      sub=res.substitute({arg:s})\n"                      text+="      sub=res.substitute({arg:s})\n"
1351               text+=mkText("array","ref",r,r1)           text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1352               res="sub"           text+="      self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1353             text+="      self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong function space of result.\")\n"%res
1354             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1355    
1356    
1357             t_prog+=text
1358    print t_prog
1359    1/0
1360    
1361    #=======================================================================================================
1362    # interpolation
1363    #=======================================================================================================
1364    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne","Solution","ReducedSolution"]:
1365      for data in ["Data","Symbol"]:
1366         for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1367          for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1368            if  where==case or \
1369                ( case in ["ContinuousFunction","Solution","ReducedSolution"] and where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"] ) or \
1370                ( case in ["FunctionOnContactZero","FunctionOnContactOne"] and where in ["FunctionOnContactZero","FunctionOnContactOne"] ) or \
1371                (case=="ContinuousFunction" and  where in ["Solution","ReducedSolution"]) or \
1372                (case=="Solution" and  where=="ReducedSolution") :
1373                
1374    
1375             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1376             tname="test_interpolation_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1377             text+="   def %s(self):\n"%tname
1378             text+="      \"\"\"\n"
1379             text+="      tests interpolation for rank %s %s onto the %s\n\n"%(len(sh),data,where)
1380             text+="      assumptions: self.domain supports inpterpolation from %s onto %s\n"%(case,where)
1381             text+="      \"\"\"\n"
1382             if case=="ReducedSolution" or where=="ReducedSolution":
1383                text+="      o=1\n"
1384                o="1"
1385             else:
1386                text+="      o=self.order\n"
1387                o="o"
1388             text+="      dim=self.domain.getDim()\n"
1389             text+="      w_ref=%s(self.domain)\n"%where
1390             text+="      x_ref=w_ref.getX()\n"
1391             text+="      w=%s(self.domain)\n"%case
1392             text+="      x=w.getX()\n"
1393             a_2=makeArray(sh+(2,),[-1.,1])
1394             b_2=makeArray(sh+(2,),[-1.,1])
1395             a_3=makeArray(sh+(3,),[-1.,1])
1396             b_3=makeArray(sh+(3,),[-1.,1])
1397             if data=="Symbol":
1398                text+="      arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1399                val="s"
1400                res="sub"
1401           else:           else:
1402               text+=mkText(case,"ref",r,r1)              val="arg"
1403               res="res"              res="res"
1404           if oper[0]=="length":                         text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
1405                 text+=mkTypeAndShapeTest(case,(),"res")           text+="      ref=Data(0,%s,w_ref)\n"%str(sh)
1406           else:                       text+="      if dim==2:\n"
1407              if case=="float" or case=="array":                   text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1408                 text+=mkTypeAndShapeTest("float",(),"res")           text+=unrollLoops(a_2,b_2,o,"ref",tap="        ",x="x_ref")
1409              else:                     text+="      else:\n"
1410                 text+=mkTypeAndShapeTest(case,(),"res")          
1411             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1412             text+=unrollLoops(a_3,b_3,o,"ref",tap="        ",x="x_ref")
1413             text+="      res=interpolate(arg,where=w_ref)\n"
1414             if data=="Symbol":
1415                text+="      sub=res.substitute({arg:s})\n"
1416             text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1417             text+="      self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong functionspace of result.\")\n"%res
1418             text+="      self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1419           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
1420           if case == "taggedData":           t_prog+=text
1421             t_prog_with_tags+=text  print test_header
1422    print t_prog
1423    print test_tail          
1424    1/0
1425    
1426    #=======================================================================================================
1427    # grad
1428    #=======================================================================================================
1429    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1430      for data in ["Data","Symbol"]:
1431         for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1432           for sh in [ (),(2,), (4,5), (6,2,2)]:
1433             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1434             tname="test_grad_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1435             text+="   def %s(self):\n"%tname
1436             text+="      \"\"\"\n"
1437             text+="      tests gradient for rank %s %s on the %s\n\n"%(len(sh),data,where)
1438             text+="      assumptions: %s(self.domain) exists\n"%case
1439             text+="                   self.domain supports gardient on %s\n"%where
1440             text+="      \"\"\"\n"
1441             if case=="ReducedSolution":
1442                text+="      o=1\n"
1443                o="1"
1444             else:
1445                text+="      o=self.order\n"
1446                o="o"
1447             text+="      dim=self.domain.getDim()\n"
1448             text+="      w_ref=%s(self.domain)\n"%where
1449             text+="      x_ref=w_ref.getX()\n"
1450             text+="      w=%s(self.domain)\n"%case
1451             text+="      x=w.getX()\n"
1452             a_2=makeArray(sh+(2,),[-1.,1])
1453             b_2=makeArray(sh+(2,),[-1.,1])
1454             a_3=makeArray(sh+(3,),[-1.,1])
1455             b_3=makeArray(sh+(3,),[-1.,1])
1456             if data=="Symbol":
1457                text+="      arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1458                val="s"
1459                res="sub"
1460           else:           else:
1461             t_prog+=text              val="arg"
1462                res="res"
1463             text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
1464             text+="      ref=Data(0,%s+(dim,),w_ref)\n"%str(sh)
1465             text+="      if dim==2:\n"
1466             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1467             text+=unrollLoopsOfGrad(a_2,b_2,o,"ref",tap="        ")
1468             text+="      else:\n"
1469            
1470             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1471             text+=unrollLoopsOfGrad(a_3,b_3,o,"ref",tap="        ")
1472             text+="      res=grad(arg,where=w_ref)\n"
1473             if data=="Symbol":
1474                text+="      sub=res.substitute({arg:s})\n"
1475             text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1476             text+="      self.failUnlessEqual(res.getShape(),%s+(dim,),\"wrong shape of result.\")\n"%str(sh)
1477             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1478    
1479    
1480             t_prog+=text
1481    print test_header
1482    print t_prog
1483    print test_tail          
1484    1/0
1485    
1486    
1487    #=======================================================================================================
1488    # integrate
1489    #=======================================================================================================
1490    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1491      for data in ["Data","Symbol"]:
1492        for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1493          for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1494            if (not case in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]) or where==case:  
1495             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1496             tname="test_integrate_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1497             text+="   def %s(self):\n"%tname
1498             text+="      \"\"\"\n"
1499             text+="      tests integral of rank %s %s on the %s\n\n"%(len(sh),data,where)
1500             text+="      assumptions: %s(self.domain) exists\n"%case
1501             text+="                   self.domain supports integral on %s\n"%where
1502    
1503             text+="      \"\"\"\n"
1504             if case=="ReducedSolution":
1505                text+="      o=1\n"
1506                o="1"
1507             else:
1508                text+="      o=self.order\n"
1509                o="o"
1510             text+="      dim=self.domain.getDim()\n"
1511             text+="      w_ref=%s(self.domain)\n"%where
1512             text+="      w=%s(self.domain)\n"%case
1513             text+="      x=w.getX()\n"
1514             a_2=makeArray(sh+(2,),[-1.,1])
1515             b_2=makeArray(sh+(2,),[-1.,1])
1516             a_3=makeArray(sh+(3,),[-1.,1])
1517             b_3=makeArray(sh+(3,),[-1.,1])
1518             if data=="Symbol":
1519                text+="      arg=Symbol(shape=%s)\n"%str(sh)
1520                val="s"
1521                res="sub"
1522             else:
1523                val="arg"
1524                res="res"
1525                
1526             text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
1527             if not len(sh)==0:
1528                text+="      ref=numarray.zeros(%s,numarray.Float)\n"%str(sh)
1529             text+="      if dim==2:\n"
1530             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1531             text+=unrollLoopsOfInteriorIntegral(a_2,b_2,where,"ref",tap="        ")
1532             text+="      else:\n"
1533            
1534             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1535             text+=unrollLoopsOfInteriorIntegral(a_3,b_3,where,"ref",tap="        ")
1536             if case in ["ContinuousFunction","Solution","ReducedSolution"]:
1537                 text+="      res=integrate(arg,where=w_ref)\n"
1538             else:
1539                 text+="      res=integrate(arg)\n"
1540    
1541             if data=="Symbol":
1542                text+="      sub=res.substitute({arg:s})\n"
1543             if len(sh)==0 and data=="Data":
1544                text+="      self.failUnless(isinstance(%s,float),\"wrong type of result.\")\n"%res
1545             else:
1546                if data=="Symbol":
1547                   text+="      self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1548                   text+="      self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1549                else:
1550                   text+="      self.failUnless(isinstance(res,numarray.NumArray),\"wrong type of result.\")\n"
1551                   text+="      self.failUnlessEqual(res.shape,%s,\"wrong shape of result.\")\n"%str(sh)
1552             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1553    
1554    
1555             t_prog+=text
1556    print test_header
1557    print t_prog
1558    print test_tail          
1559    1/0
1560    #=======================================================================================================
1561    # inverse
1562    #=======================================================================================================
1563    name="inverse"
1564    for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1565      for sh0 in [ (1,1), (2,2), (3,3)]:
1566                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1567                  tname="test_%s_%s_dim%s"%(name,case0,sh0[0])
1568                  text+="   def %s(self):\n"%tname
1569                  a_0=makeArray(sh0,[-1.,1])
1570                  for i in range(sh0[0]): a_0[i,i]+=2
1571                  if case0 in ["taggedData", "expandedData"]:
1572                      a1_0=makeArray(sh0,[-1.,1])
1573                      for i in range(sh0[0]): a1_0[i,i]+=3
1574                  else:
1575                      a1_0=a_0
1576                      
1577                  text+=mkText(case0,"arg",a_0,a1_0)
1578                  text+="      res=%s(arg)\n"%name
1579                  if case0=="Symbol":
1580                     text+=mkText("array","s",a_0,a1_0)
1581                     text+="      sub=res.substitute({arg:s})\n"
1582                     res="sub"
1583                     ref="s"
1584                  else:
1585                     ref="arg"
1586                     res="res"
1587                  text+=mkTypeAndShapeTest(case0,sh0,"res")
1588                  text+="      self.failUnless(Lsup(matrixmult(%s,%s)-kronecker(%s))<=self.RES_TOL,\"wrong result\")\n"%(res,ref,sh0[0])
1589                  
1590                  if case0 == "taggedData" :
1591                      t_prog_with_tags+=text
1592                  else:              
1593                      t_prog+=text
1594    
1595  print test_header  print test_header
1596  # print t_prog  # print t_prog
1597  print t_prog_with_tags  print t_prog_with_tags
1598  print test_tail            print test_tail          
1599  1/0  1/0
1600              
1601  #=======================================================================================================  #=======================================================================================================
1602  # tensor multiply  # trace
1603  #=======================================================================================================  #=======================================================================================================
1604  # oper=["generalTensorProduct",tensorProductTest]  def traceTest(r,offset):
1605  # oper=["matrixmult",testMatrixMult]      sh=r.shape
1606  oper=["tensormult",testTensorMult]      r1=1
1607        for i in range(offset): r1*=sh[i]
1608        r2=1
1609        for i in range(offset+2,len(sh)): r2*=sh[i]
1610        r_s=numarray.reshape(r,(r1,sh[offset],sh[offset],r2))
1611        s=numarray.zeros([r1,r2],numarray.Float)
1612        for i1 in range(r1):
1613            for i2 in range(r2):
1614                for j in range(sh[offset]): s[i1,i2]+=r_s[i1,j,j,i2]
1615        return s.resize(sh[:offset]+sh[offset+2:])
1616    name,tt="trace",traceTest
1617    for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1618      for sh0 in [ (4,5), (6,2,2),(3,2,3,4)]:
1619        for offset in range(len(sh0)-1):
1620                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1621                  tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1622                  text+="   def %s(self):\n"%tname
1623                  sh_t=list(sh0)
1624                  sh_t[offset+1]=sh_t[offset]
1625                  sh_t=tuple(sh_t)
1626                  sh_r=[]
1627                  for i in range(offset): sh_r.append(sh0[i])
1628                  for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])              
1629                  sh_r=tuple(sh_r)
1630                  a_0=makeArray(sh_t,[-1.,1])
1631                  if case0 in ["taggedData", "expandedData"]:
1632                      a1_0=makeArray(sh_t,[-1.,1])
1633                  else:
1634                      a1_0=a_0
1635                  r=tt(a_0,offset)
1636                  r1=tt(a1_0,offset)
1637                  text+=mkText(case0,"arg",a_0,a1_0)
1638                  text+="      res=%s(arg,%s)\n"%(name,offset)
1639                  if case0=="Symbol":
1640                     text+=mkText("array","s",a_0,a1_0)
1641                     text+="      sub=res.substitute({arg:s})\n"
1642                     res="sub"
1643                     text+=mkText("array","ref",r,r1)
1644                  else:
1645                     res="res"
1646                     text+=mkText(case0,"ref",r,r1)
1647                  text+=mkTypeAndShapeTest(case0,sh_r,"res")
1648                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1649                  
1650                  if case0 == "taggedData" :
1651                      t_prog_with_tags+=text
1652                  else:              
1653                      t_prog+=text
1654    
1655  for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:  print test_header
1656    # print t_prog
1657    print t_prog_with_tags
1658    print test_tail          
1659    1/0
1660    
1661    #=======================================================================================================
1662    # clip
1663    #=======================================================================================================
1664    oper_L=[["clip",clipTEST]]
1665    for oper in oper_L:
1666     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1667    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)]:
1668            if len(sh0)==0 or not case0=="float":
1669                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1670                  tname="test_%s_%s_rank%s"%(oper[0],case0,len(sh0))
1671                  text+="   def %s(self):\n"%tname
1672                  a_0=makeArray(sh0,[-1.,1])
1673                  if case0 in ["taggedData", "expandedData"]:
1674                      a1_0=makeArray(sh0,[-1.,1])
1675                  else:
1676                      a1_0=a_0
1677    
1678                  r=oper[1](a_0,-0.3,0.5)
1679                  r1=oper[1](a1_0,-0.3,0.5)
1680                  text+=mkText(case0,"arg",a_0,a1_0)
1681                  text+="      res=%s(arg,-0.3,0.5)\n"%oper[0]
1682                  if case0=="Symbol":
1683                     text+=mkText("array","s",a_0,a1_0)
1684                     text+="      sub=res.substitute({arg:s})\n"
1685                     res="sub"
1686                     text+=mkText("array","ref",r,r1)
1687                  else:
1688                     res="res"
1689                     text+=mkText(case0,"ref",r,r1)
1690                  text+=mkTypeAndShapeTest(case0,sh0,"res")
1691                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1692                  
1693                  if case0 == "taggedData" :
1694                      t_prog_with_tags+=text
1695                  else:              
1696                      t_prog+=text
1697    
1698    print test_header
1699    # print t_prog
1700    print t_prog_with_tags
1701    print test_tail          
1702    1/0
1703    
1704    #=======================================================================================================
1705    # maximum, minimum, clipping
1706    #=======================================================================================================
1707    oper_L=[ ["maximum",maximumTEST],
1708             ["minimum",minimumTEST]]
1709    for oper in oper_L:
1710     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1711      for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1712     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1713       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)]:
1714         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") \
1715            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)    
1716                use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"                use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1717    
1718                text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"                text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1719                # 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  
1720                text+="   def %s(self):\n"%tname                text+="   def %s(self):\n"%tname
1721                a_0=makeArray(sh0+sh_s,[-1.,1])                a_0=makeArray(sh0,[-1.,1])
1722                if case0 in ["taggedData", "expandedData"]:                if case0 in ["taggedData", "expandedData"]:
1723                    a1_0=makeArray(sh0+sh_s,[-1.,1])                    a1_0=makeArray(sh0,[-1.,1])
1724                else:                else:
1725                    a1_0=a_0                    a1_0=a_0
1726    
1727                a_1=makeArray(sh_s+sh1,[-1.,1])                a_1=makeArray(sh1,[-1.,1])
1728                if case1 in ["taggedData", "expandedData"]:                if case1 in ["taggedData", "expandedData"]:
1729                    a1_1=makeArray(sh_s+sh1,[-1.,1])                    a1_1=makeArray(sh1,[-1.,1])
1730                else:                else:
1731                    a1_1=a_1                    a1_1=a_1
1732                r=oper[1](a_0,a_1,sh_s)                r=oper[1](a_0,a_1)
1733                r1=oper[1](a1_0,a1_1,sh_s)                r1=oper[1](a1_0,a1_1)
1734                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)
1735                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)
1736                #text+="      res=matrixmult(arg0,arg1)\n"                text+="      res=%s(arg0,arg1)\n"%oper[0]
1737                text+="      res=tensormult(arg0,arg1)\n"                case=getResultCaseForBin(case0,case1)              
               #text+="      res=generalTensorProduct(arg0,arg1,offset=%s)\n"%(len(sh_s))  
1738                if case=="Symbol":                if case=="Symbol":
1739                   c0_res,c1_res=case0,case1                   c0_res,c1_res=case0,case1
1740                   subs="{"                   subs="{"
# Line 724  for case0 in ["float","array","Symbol"," Line 1754  for case0 in ["float","array","Symbol","
1754                else:                else:
1755                   res="res"                   res="res"
1756                   text+=mkText(case,"ref",r,r1)                   text+=mkText(case,"ref",r,r1)
1757                text+=mkTypeAndShapeTest(case,sh0+sh1,"res")                if len(sh0)>len(sh1):
1758                      text+=mkTypeAndShapeTest(case,sh0,"res")
1759                  else:
1760                      text+=mkTypeAndShapeTest(case,sh1,"res")
1761                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
1762                  
1763                if case0 == "taggedData" or case1 == "taggedData":                if case0 == "taggedData" or case1 == "taggedData":
1764                    t_prog_with_tags+=text                    t_prog_with_tags+=text
1765                else:                              else:              
1766                    t_prog+=text                    t_prog+=text
1767    
1768  print test_header  print test_header
1769  # print t_prog  # print t_prog
1770  print t_prog_with_tags  print t_prog_with_tags
1771  print test_tail            print test_tail          
1772  1/0  1/0
1773    
1774    
1775  #=======================================================================================================  #=======================================================================================================
1776  # outer/inner  # outer inner
1777  #=======================================================================================================  #=======================================================================================================
1778  oper=["inner",innerTEST]  oper=["outer",outerTEST]
1779  # oper=["outer",outerTEST]  # oper=["inner",innerTEST]
1780  for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:  for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1781    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)]:
1782     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 1835  for case0 in ["float","array","Symbol","
1835  print test_header  print test_header
1836  # print t_prog  # print t_prog
1837  print t_prog_with_tags  print t_prog_with_tags
1838    print test_tail          
1839    1/0
1840    
1841    #=======================================================================================================
1842    # local reduction
1843    #=======================================================================================================
1844    for oper in [["length",0.,"out+%a1%**2","math.sqrt(out)"],
1845                 ["maxval",-1.e99,"max(out,%a1%)","out"],
1846                 ["minval",1.e99,"min(out,%a1%)","out"] ]:
1847      for case in case_set:
1848         for sh in shape_set:
1849           if not case=="float" or len(sh)==0:
1850             text=""
1851             text+="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1852             tname="def test_%s_%s_rank%s"%(oper[0],case,len(sh))
1853             text+="   %s(self):\n"%tname
1854             a=makeArray(sh,[-1.,1.])            
1855             a1=makeArray(sh,[-1.,1.])
1856             r1=testReduce(a1,oper[1],oper[2],oper[3])
1857             r=testReduce(a,oper[1],oper[2],oper[3])
1858            
1859             text+=mkText(case,"arg",a,a1)
1860             text+="      res=%s(arg)\n"%oper[0]
1861             if case=="Symbol":        
1862                 text+=mkText("array","s",a,a1)
1863                 text+="      sub=res.substitute({arg:s})\n"        
1864                 text+=mkText("array","ref",r,r1)
1865                 res="sub"
1866             else:
1867                 text+=mkText(case,"ref",r,r1)
1868                 res="res"
1869             if oper[0]=="length":              
1870                   text+=mkTypeAndShapeTest(case,(),"res")
1871             else:            
1872                if case=="float" or case=="array":        
1873                   text+=mkTypeAndShapeTest("float",(),"res")
1874                else:          
1875                   text+=mkTypeAndShapeTest(case,(),"res")
1876             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1877             if case == "taggedData":
1878               t_prog_with_tags+=text
1879             else:
1880               t_prog+=text
1881    print test_header
1882    # print t_prog
1883    print t_prog_with_tags
1884    print test_tail          
1885    1/0
1886              
1887    #=======================================================================================================
1888    # tensor multiply
1889    #=======================================================================================================
1890    # oper=["generalTensorProduct",tensorProductTest]
1891    # oper=["matrixmult",testMatrixMult]
1892    oper=["tensormult",testTensorMult]
1893    
1894    for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1895      for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1896       for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1897         for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1898           for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:
1899              if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \
1900                   and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
1901                # 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
1902                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
1903                  case=getResultCaseForBin(case0,case1)  
1904                  use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1905                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1906                  # tname="test_generalTensorProduct_%s_rank%s_%s_rank%s_offset%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1),len(sh_s))
1907                  #tname="test_matrixmult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
1908                  tname="test_tensormult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
1909                  # if tname=="test_generalTensorProduct_array_rank1_array_rank2_offset1":
1910                  # print tnametest_generalTensorProduct_Symbol_rank1_Symbol_rank3_offset1
1911                  text+="   def %s(self):\n"%tname
1912                  a_0=makeArray(sh0+sh_s,[-1.,1])
1913                  if case0 in ["taggedData", "expandedData"]:
1914                      a1_0=makeArray(sh0+sh_s,[-1.,1])
1915                  else:
1916                      a1_0=a_0
1917    
1918                  a_1=makeArray(sh_s+sh1,[-1.,1])
1919                  if case1 in ["taggedData", "expandedData"]:
1920                      a1_1=makeArray(sh_s+sh1,[-1.,1])
1921                  else:
1922                      a1_1=a_1
1923                  r=oper[1](a_0,a_1,sh_s)
1924                  r1=oper[1](a1_0,a1_1,sh_s)
1925                  text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1926                  text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1927                  #text+="      res=matrixmult(arg0,arg1)\n"
1928                  text+="      res=tensormult(arg0,arg1)\n"
1929                  #text+="      res=generalTensorProduct(arg0,arg1,offset=%s)\n"%(len(sh_s))
1930                  if case=="Symbol":
1931                     c0_res,c1_res=case0,case1
1932                     subs="{"
1933                     if case0=="Symbol":        
1934                        text+=mkText("array","s0",a_0,a1_0)
1935                        subs+="arg0:s0"
1936                        c0_res="array"
1937                     if case1=="Symbol":        
1938                        text+=mkText("array","s1",a_1,a1_1)
1939                        if not subs.endswith("{"): subs+=","
1940                        subs+="arg1:s1"
1941                        c1_res="array"
1942                     subs+="}"  
1943                     text+="      sub=res.substitute(%s)\n"%subs
1944                     res="sub"
1945                     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1946                  else:
1947                     res="res"
1948                     text+=mkText(case,"ref",r,r1)
1949                  text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
1950                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1951                  if case0 == "taggedData" or case1 == "taggedData":
1952                      t_prog_with_tags+=text
1953                  else:              
1954                      t_prog+=text
1955    print test_header
1956    # print t_prog
1957    print t_prog_with_tags
1958  print test_tail            print test_tail          
1959  1/0  1/0
1960  #=======================================================================================================  #=======================================================================================================

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

  ViewVC Help
Powered by ViewVC 1.1.26