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

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

  ViewVC Help
Powered by ViewVC 1.1.26