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

Diff of /trunk/escript/py_src/generateutil

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

revision 433 by gross, Tue Jan 17 23:54:38 2006 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 689  def minimumTEST(arg0,arg1): Line 901  def minimumTEST(arg0,arg1):
901                else:                else:
902                 out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]                 out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
903       return out       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    # transpose
1230    #=======================================================================================================
1231    def transposeTest(r,offset):
1232        if isinstance(r,float): return r
1233        s=r.shape
1234        s1=1
1235        for i in s[:offset]: s1*=i
1236        s2=1
1237        for i in s[offset:]: s2*=i
1238        out=numarray.reshape(r,(s1,s2))
1239        out.transpose()
1240        return numarray.resize(out,s[offset:]+s[:offset])
1241    
1242    name,tt="transpose",transposeTest
1243    for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1244      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+=unrollLoops(a_3,b_3,o,val,tap="        ")
1391             text+=unrollLoopsOfDiv(a_3,b_3,o,"ref",tap="        ")
1392             text+="\n      res=div(arg,where=w_ref)\n"
1393             if data=="Symbol":
1394                text+="      sub=res.substitute({arg:s})\n"
1395             text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1396             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:
1446                val="arg"
1447                res="res"
1448             text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
1449             text+="      ref=Data(0,%s,w_ref)\n"%str(sh)
1450             text+="      if dim==2:\n"
1451             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1452             text+=unrollLoops(a_2,b_2,o,"ref",tap="        ",x="x_ref")
1453             text+="      else:\n"
1454            
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
1464             t_prog+=text
1465    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:
1505                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  # inverse
1606  #=======================================================================================================  #=======================================================================================================

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

  ViewVC Help
Powered by ViewVC 1.1.26