# Diff of /trunk/escript/py_src/generatediff

revision 439 by gross, Fri Jan 20 01:46:22 2006 UTC revision 530 by gross, Wed Feb 15 07:11:00 2006 UTC
# Line 184  def makeArray(shape,rng): Line 184  def makeArray(shape,rng):
184         raise SystemError,"rank is restricted to 5"         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 524  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
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
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 696  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"):  def unrollLoops(a,b,o,arg,tap="",x="x"):
906      out=""      out=""
907      if a.rank==1:      if a.rank==1:
1003                 else:                 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])                   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      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=""):  def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
1023      if where=="Function":      if where=="Function":
1024         xfac_o=1.         xfac_o=1.
# Line 916  def unrollLoopsOfInteriorIntegral(a,b,wh Line 1138  def unrollLoopsOfInteriorIntegral(a,b,wh
1138                 out+="+(%s)*0.5**o\n"%zop                 out+="+(%s)*0.5**o\n"%zop
1139
1140      return out      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
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  # interpolation

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