# Diff of /trunk/escript/py_src/generatediff

revision 283 by gross, Wed Nov 30 23:35:29 2005 UTC revision 439 by gross, Fri Jan 20 01:46:22 2006 UTC
# Line 33  def wherepos(arg): Line 33  def wherepos(arg):
33     else:     else:
34        return 0.        return 0.
35
36
37  class OPERATOR:  class OPERATOR:
38      def __init__(self,nickname,rng=[-1000.,1000],test_expr="",math_expr=None,      def __init__(self,nickname,rng=[-1000.,1000],test_expr="",math_expr=None,
39                  numarray_expr="",symbol_expr=None,diff=None,name=""):                  numarray_expr="",symbol_expr=None,diff=None,name=""):
# Line 172  def makeArray(shape,rng): Line 173  def makeArray(shape,rng):
173               for i2 in range(shape[2]):               for i2 in range(shape[2]):
174                  for i3 in range(shape[3]):                  for i3 in range(shape[3]):
175                     out[i0,i1,i2,i3]=l*random.random()+rng[0]                     out[i0,i1,i2,i3]=l*random.random()+rng[0]
176       elif len(shape)==5:
177           for i0 in range(shape[0]):
178              for i1 in range(shape[1]):
179                 for i2 in range(shape[2]):
180                    for i3 in range(shape[3]):
181                      for i4 in range(shape[4]):
182                       out[i0,i1,i2,i3,i4]=l*random.random()+rng[0]
183     else:     else:
184         raise SystemError,"rank is restricted to 4"         raise SystemError,"rank is restricted to 5"
185     return out             return out
186
187
# Line 431  def mkText(case,name,a,a1=None,use_taggi Line 439  def mkText(case,name,a,a1=None,use_taggi
439           if case=="float":           if case=="float":
440             if isinstance(a,float):             if isinstance(a,float):
441                  t_out+="      %s=%s\n"%(name,a)                  t_out+="      %s=%s\n"%(name,a)
442             elif len(a)==1:             elif a.rank==0:
443                  t_out+="      %s=%s\n"%(name,a)                  t_out+="      %s=%s\n"%(name,a)
444             else:             else:
445                  t_out+="      %s=numarray.array(%s)\n"%(name,a.tolist())                  t_out+="      %s=numarray.array(%s)\n"%(name,a.tolist())
446           elif case=="array":           elif case=="array":
447             if isinstance(a,float):             if isinstance(a,float):
448                  t_out+="      %s=numarray.array(%s)\n"%(name,a)                  t_out+="      %s=numarray.array(%s)\n"%(name,a)
449             elif len(a)==1:             elif a.rank==0:
450                  t_out+="      %s=numarray.array(%s)\n"%(name,a)                  t_out+="      %s=numarray.array(%s)\n"%(name,a)
451             else:             else:
452                  t_out+="      %s=numarray.array(%s)\n"%(name,a.tolist())                  t_out+="      %s=numarray.array(%s)\n"%(name,a.tolist())
453           elif case=="constData":           elif case=="constData":
454             if isinstance(a,float):             if isinstance(a,float):
455                t_out+="      %s=Data(%s,self.functionspace)\n"%(name,a)                t_out+="      %s=Data(%s,self.functionspace)\n"%(name,a)
456             elif len(a)==1:             elif a.rank==0:
457                t_out+="      %s=Data(%s,self.functionspace)\n"%(name,a)                t_out+="      %s=Data(%s,self.functionspace)\n"%(name,a)
458             else:             else:
459                t_out+="      %s=Data(numarray.array(%s),self.functionspace)\n"%(name,a.tolist())                t_out+="      %s=Data(numarray.array(%s),self.functionspace)\n"%(name,a.tolist())
# Line 453  def mkText(case,name,a,a1=None,use_taggi Line 461  def mkText(case,name,a,a1=None,use_taggi
461             if isinstance(a,float):             if isinstance(a,float):
462                t_out+="      %s=Data(%s,self.functionspace)\n"%(name,a)                t_out+="      %s=Data(%s,self.functionspace)\n"%(name,a)
463                t_out+="      %s.setTaggedValue(1,%s)\n"%(name,a1)                t_out+="      %s.setTaggedValue(1,%s)\n"%(name,a1)
464             elif len(a)==1:             elif a.rank==0:
465                t_out+="      %s=Data(%s,self.functionspace)\n"%(name,a)                t_out+="      %s=Data(%s,self.functionspace)\n"%(name,a)
466                t_out+="      %s.setTaggedValue(1,%s)\n"%(name,a1)                t_out+="      %s.setTaggedValue(1,%s)\n"%(name,a1)
467             else:             else:
# Line 464  def mkText(case,name,a,a1=None,use_taggi Line 472  def mkText(case,name,a,a1=None,use_taggi
472                if isinstance(a,float):                if isinstance(a,float):
473                   t_out+="      %s=Data(%s,self.functionspace)\n"%(name,a)                   t_out+="      %s=Data(%s,self.functionspace)\n"%(name,a)
474                   t_out+="      %s.setTaggedValue(1,%s)\n"%(name,a1)                   t_out+="      %s.setTaggedValue(1,%s)\n"%(name,a1)
475                elif len(a)==1:                elif a.rank==0:
476                   t_out+="      %s=Data(%s,self.functionspace)\n"%(name,a)                   t_out+="      %s=Data(%s,self.functionspace)\n"%(name,a)
477                   t_out+="      %s.setTaggedValue(1,%s)\n"%(name,a1)                   t_out+="      %s.setTaggedValue(1,%s)\n"%(name,a1)
478                else:                else:
# Line 475  def mkText(case,name,a,a1=None,use_taggi Line 483  def mkText(case,name,a,a1=None,use_taggi
483                t_out+="      msk_%s=whereNegative(self.functionspace.getX()[0]-0.5)\n"%name                t_out+="      msk_%s=whereNegative(self.functionspace.getX()[0]-0.5)\n"%name
484                if isinstance(a,float):                if isinstance(a,float):
485                     t_out+="      %s=msk_%s*(%s)+(1.-msk_%s)*(%s)\n"%(name,name,a,name,a1)                     t_out+="      %s=msk_%s*(%s)+(1.-msk_%s)*(%s)\n"%(name,name,a,name,a1)
486                elif len(a)==1:                elif a.rank==0:
487                     t_out+="      %s=msk_%s*numarray.array(%s)+(1.-msk_%s)*numarray.array(%s)\n"%(name,name,a,name,a1)                     t_out+="      %s=msk_%s*numarray.array(%s)+(1.-msk_%s)*numarray.array(%s)\n"%(name,name,a,name,a1)
488                else:                else:
489                     t_out+="      %s=msk_%s*numarray.array(%s)+(1.-msk_%s)*numarray.array(%s)\n"%(name,name,a.tolist(),name,a1.tolist())                     t_out+="      %s=msk_%s*numarray.array(%s)+(1.-msk_%s)*numarray.array(%s)\n"%(name,name,a.tolist(),name,a1.tolist())
490           elif case=="Symbol":           elif case=="Symbol":
491             if isinstance(a,float):             if isinstance(a,float):
492                t_out+="      %s=Symbol(shape=())\n"%(name)                t_out+="      %s=Symbol(shape=())\n"%(name)
493             elif len(a)==1:             elif a.rank==0:
494                t_out+="      %s=Symbol(shape=())\n"%(name)                t_out+="      %s=Symbol(shape=())\n"%(name)
495             else:             else:
496                t_out+="      %s=Symbol(shape=%s)\n"%(name,str(a.shape))                t_out+="      %s=Symbol(shape=%s)\n"%(name,str(a.shape))
# Line 517  def mkCode(txt,args=[],intend=""): Line 525  def mkCode(txt,args=[],intend=""):
525        out=out.replace("%%a%s%%"%c,r)        out=out.replace("%%a%s%%"%c,r)
526      return out        return out
527
528    def innerTEST(arg0,arg1):
529        if isinstance(arg0,float):
530           out=numarray.array(arg0*arg1)
531        else:
532           out=(arg0*arg1).sum()
533        return out
534
535    def outerTEST(arg0,arg1):
536        if isinstance(arg0,float):
537           out=numarray.array(arg0*arg1)
538        elif isinstance(arg1,float):
539           out=numarray.array(arg0*arg1)
540        else:
541           out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
542        return out
543
544    def tensorProductTest(arg0,arg1,sh_s):
545        if isinstance(arg0,float):
546           out=numarray.array(arg0*arg1)
547        elif isinstance(arg1,float):
548           out=numarray.array(arg0*arg1)
549        elif len(sh_s)==0:
550           out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
551        else:
552           l=len(sh_s)
553           sh0=arg0.shape[:arg0.rank-l]
554           sh1=arg1.shape[l:]
555           ls,l0,l1=1,1,1
556           for i in sh_s: ls*=i
557           for i in sh0: l0*=i
558           for i in sh1: l1*=i
559           out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,ls,l1))
560           out2=numarray.zeros((l0,l1),numarray.Float)
561           for i0 in range(l0):
562              for i1 in range(l1):
563                  for i in range(ls): out2[i0,i1]+=out1[i0,i,i,i1]
564           out=out2.resize(sh0+sh1)
565        return out
566
567    def testMatrixMult(arg0,arg1,sh_s):
568         return numarray.matrixmultiply(arg0,arg1)
569
570
571    def testTensorMult(arg0,arg1,sh_s):
572         if len(arg0)==2:
573            return numarray.matrixmultiply(arg0,arg1)
574         else:
575            if arg1.rank==4:
576              out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2],arg1.shape[3]),numarray.Float)
577              for i0 in range(arg0.shape[0]):
578               for i1 in range(arg0.shape[1]):
579                for i2 in range(arg0.shape[2]):
580                 for i3 in range(arg0.shape[3]):
581                  for j2 in range(arg1.shape[2]):
582                   for j3 in range(arg1.shape[3]):
583                         out[i0,i1,j2,j3]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2,j3]
584            elif arg1.rank==3:
585              out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2]),numarray.Float)
586              for i0 in range(arg0.shape[0]):
587               for i1 in range(arg0.shape[1]):
588                for i2 in range(arg0.shape[2]):
589                 for i3 in range(arg0.shape[3]):
590                  for j2 in range(arg1.shape[2]):
591                         out[i0,i1,j2]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2]
592            elif arg1.rank==2:
593              out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
594              for i0 in range(arg0.shape[0]):
595               for i1 in range(arg0.shape[1]):
596                for i2 in range(arg0.shape[2]):
597                 for i3 in range(arg0.shape[3]):
598                         out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
599            return out
600
601    def testReduce(arg0,init_val,test_expr,post_expr):
602         out=init_val
603         if isinstance(arg0,float):
604              out=eval(test_expr.replace("%a1%","arg0"))
605         elif arg0.rank==0:
606              out=eval(test_expr.replace("%a1%","arg0"))
607         elif arg0.rank==1:
608            for i0 in range(arg0.shape[0]):
609                   out=eval(test_expr.replace("%a1%","arg0[i0]"))
610         elif arg0.rank==2:
611            for i0 in range(arg0.shape[0]):
612             for i1 in range(arg0.shape[1]):
613                   out=eval(test_expr.replace("%a1%","arg0[i0,i1]"))
614         elif arg0.rank==3:
615            for i0 in range(arg0.shape[0]):
616             for i1 in range(arg0.shape[1]):
617               for i2 in range(arg0.shape[2]):
618                   out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2]"))
619         elif arg0.rank==4:
620            for i0 in range(arg0.shape[0]):
621             for i1 in range(arg0.shape[1]):
622               for i2 in range(arg0.shape[2]):
623                 for i3 in range(arg0.shape[3]):
624                   out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2,i3]"))
625         return eval(post_expr)
626
627    def clipTEST(arg0,mn,mx):
628         if isinstance(arg0,float):
629              return max(min(arg0,mx),mn)
630         out=numarray.zeros(arg0.shape,numarray.Float64)
631         if arg0.rank==1:
632            for i0 in range(arg0.shape[0]):
633                out[i0]=max(min(arg0[i0],mx),mn)
634         elif arg0.rank==2:
635            for i0 in range(arg0.shape[0]):
636             for i1 in range(arg0.shape[1]):
637                out[i0,i1]=max(min(arg0[i0,i1],mx),mn)
638         elif arg0.rank==3:
639            for i0 in range(arg0.shape[0]):
640             for i1 in range(arg0.shape[1]):
641               for i2 in range(arg0.shape[2]):
642                  out[i0,i1,i2]=max(min(arg0[i0,i1,i2],mx),mn)
643         elif arg0.rank==4:
644            for i0 in range(arg0.shape[0]):
645             for i1 in range(arg0.shape[1]):
646               for i2 in range(arg0.shape[2]):
647                 for i3 in range(arg0.shape[3]):
648                    out[i0,i1,i2,i3]=max(min(arg0[i0,i1,i2,i3],mx),mn)
649         return out
650    def minimumTEST(arg0,arg1):
651         if isinstance(arg0,float):
652           if isinstance(arg1,float):
653              if arg0>arg1:
654                  return arg1
655              else:
656                  return arg0
657           else:
658              arg0=numarray.ones(arg1.shape)*arg0
659         else:
660           if isinstance(arg1,float):
661              arg1=numarray.ones(arg0.shape)*arg1
662         out=numarray.zeros(arg0.shape,numarray.Float64)
663         if arg0.rank==0:
664              if arg0>arg1:
665                  out=arg1
666              else:
667                  out=arg0
668         elif arg0.rank==1:
669            for i0 in range(arg0.shape[0]):
670              if arg0[i0]>arg1[i0]:
671                  out[i0]=arg1[i0]
672              else:
673                  out[i0]=arg0[i0]
674         elif arg0.rank==2:
675            for i0 in range(arg0.shape[0]):
676             for i1 in range(arg0.shape[1]):
677              if arg0[i0,i1]>arg1[i0,i1]:
678                  out[i0,i1]=arg1[i0,i1]
679              else:
680                  out[i0,i1]=arg0[i0,i1]
681         elif arg0.rank==3:
682            for i0 in range(arg0.shape[0]):
683             for i1 in range(arg0.shape[1]):
684               for i2 in range(arg0.shape[2]):
685                 if arg0[i0,i1,i2]>arg1[i0,i1,i2]:
686                  out[i0,i1,i2]=arg1[i0,i1,i2]
687                 else:
688                  out[i0,i1,i2]=arg0[i0,i1,i2]
689         elif arg0.rank==4:
690            for i0 in range(arg0.shape[0]):
691             for i1 in range(arg0.shape[1]):
692               for i2 in range(arg0.shape[2]):
693                 for i3 in range(arg0.shape[3]):
694                  if arg0[i0,i1,i2,i3]>arg1[i0,i1,i2,i3]:
695                   out[i0,i1,i2,i3]=arg1[i0,i1,i2,i3]
696                  else:
697                   out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
698         return out
699    def unrollLoops(a,b,o,arg,tap="",x="x"):
700        out=""
701        if a.rank==1:
702                 z=""
703                 for i99 in range(a.shape[0]):
704                   if not z=="": z+="+"
705                   if o=="1":
706                      z+="(%s)*%s[%s]"%(a[i99]+b[i99],x,i99)
707                   else:
708                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i99],x,i99,b[i99],x,i99)
709
710                 out+=tap+"%s=%s\n"%(arg,z)
711
712        elif a.rank==2:
713            for i0 in range(a.shape[0]):
714                 z=""
715                 for i99 in range(a.shape[1]):
716                   if not z=="": z+="+"
717                   if o=="1":
718                      z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)
719                   else:
720                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i99],x,i99,b[i0,i99],x,i99)
721
722                 out+=tap+"%s[%s]=%s\n"%(arg,i0,z)
723        elif a.rank==3:
724            for i0 in range(a.shape[0]):
725             for i1 in range(a.shape[1]):
726                 z=""
727                 for i99 in range(a.shape[2]):
728                   if not z=="": z+="+"
729                   if o=="1":
730                      z+="(%s)*%s[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],x,i99)
731                   else:
732                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i99],x,i99,b[i0,i1,i99],x,i99)
733
734                 out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)
735        elif a.rank==4:
736            for i0 in range(a.shape[0]):
737             for i1 in range(a.shape[1]):
738               for i2 in range(a.shape[2]):
739                 z=""
740                 for i99 in range(a.shape[3]):
741                   if not z=="": z+="+"
742                   if o=="1":
743                      z+="(%s)*%s[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],x,i99)
744                   else:
745                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i99],x,i99,b[i0,i1,i2,i99],x,i99)
746
747                 out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)
748        elif a.rank==5:
749            for i0 in range(a.shape[0]):
750             for i1 in range(a.shape[1]):
751               for i2 in range(a.shape[2]):
752                for i3 in range(a.shape[3]):
753                 z=""
754                 for i99 in range(a.shape[4]):
755                   if not z=="": z+="+"
756                   if o=="1":
757                      z+="(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],x,i99)
758                   else:
759                      z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99],x,i99,b[i0,i1,i2,i3,i99],x,i99)
760
761                 out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
762        return out
763
765        out=""
766        if a.rank==1:
767                 for i99 in range(a.shape[0]):
768                   if o=="1":
769                      out+=tap+"%s[%s]=(%s)\n"%(arg,i99,a[i99]+b[i99])
770                   else:
771                      out+=tap+"%s[%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i99,a[i99],i99,b[i99])
772
773        elif a.rank==2:
774            for i0 in range(a.shape[0]):
775                 for i99 in range(a.shape[1]):
776                   if o=="1":
777                      out+=tap+"%s[%s,%s]=(%s)\n"%(arg,i0,i99,a[i0,i99]+b[i0,i99])
778                   else:
779                      out+=tap+"%s[%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i99,a[i0,i99],i99,b[i0,i99])
780
781        elif a.rank==3:
782            for i0 in range(a.shape[0]):
783             for i1 in range(a.shape[1]):
784                 for i99 in range(a.shape[2]):
785                   if o=="1":
786                      out+=tap+"%s[%s,%s,%s]=(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99]+b[i0,i1,i99])
787                   else:
788                      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])
789
790        elif a.rank==4:
791            for i0 in range(a.shape[0]):
792             for i1 in range(a.shape[1]):
793               for i2 in range(a.shape[2]):
794                 for i99 in range(a.shape[3]):
795                   if o=="1":
796                     out+=tap+"%s[%s,%s,%s,%s]=(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99]+b[i0,i1,i2,i99])
797                   else:
798                     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])
799        return out
800    def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
801        if where=="Function":
802           xfac_o=1.
803           xfac_op=0.
804           z_fac=1./2.
805           z_fac_s=""
806           zo_fac_s="/(o+1.)"
807        elif where=="FunctionOnBoundary":
808           xfac_o=1.
809           xfac_op=0.
810           z_fac=1.
811           z_fac_s="*dim"
812           zo_fac_s="*(1+2.*(dim-1.)/(o+1.))"
813        elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
814           xfac_o=0.
815           xfac_op=1.
816           z_fac=1./2.
817           z_fac_s=""
818           zo_fac_s="/(o+1.)"
819        out=""
820        if a.rank==1:
821                 zo=0.
822                 zop=0.
823                 z=0.
824                 for i99 in range(a.shape[0]):
825                      if i99==0:
826                        zo+=       xfac_o*a[i99]
827                        zop+=       xfac_op*a[i99]
828                      else:
829                        zo+=a[i99]
830                      z+=b[i99]
831
832                 out+=tap+"%s=(%s)%s+(%s)%s"%(arg,zo,zo_fac_s,z*z_fac,z_fac_s)
833                 if zop==0.:
834                   out+="\n"
835                 else:
836                   out+="+(%s)*0.5**o\n"%zop
837        elif a.rank==2:
838            for i0 in range(a.shape[0]):
839                 zo=0.
840                 zop=0.
841                 z=0.
842                 for i99 in range(a.shape[1]):
843                      if i99==0:
844                        zo+=       xfac_o*a[i0,i99]
845                        zop+=       xfac_op*a[i0,i99]
846                      else:
847                        zo+=a[i0,i99]
848                      z+=b[i0,i99]
849
850                 out+=tap+"%s[%s]=(%s)%s+(%s)%s"%(arg,i0,zo,zo_fac_s,z*z_fac,z_fac_s)
851                 if zop==0.:
852                   out+="\n"
853                 else:
854                   out+="+(%s)*0.5**o\n"%zop
855        elif a.rank==3:
856            for i0 in range(a.shape[0]):
857             for i1 in range(a.shape[1]):
858                 zo=0.
859                 zop=0.
860                 z=0.
861                 for i99 in range(a.shape[2]):
862                      if i99==0:
863                        zo+=       xfac_o*a[i0,i1,i99]
864                        zop+=       xfac_op*a[i0,i1,i99]
865                      else:
866                        zo+=a[i0,i1,i99]
867                      z+=b[i0,i1,i99]
868
869                 out+=tap+"%s[%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,zo,zo_fac_s,z*z_fac,z_fac_s)
870                 if zop==0.:
871                   out+="\n"
872                 else:
873                   out+="+(%s)*0.5**o\n"%zop
874        elif a.rank==4:
875            for i0 in range(a.shape[0]):
876             for i1 in range(a.shape[1]):
877               for i2 in range(a.shape[2]):
878                 zo=0.
879                 zop=0.
880                 z=0.
881                 for i99 in range(a.shape[3]):
882                      if i99==0:
883                        zo+=       xfac_o*a[i0,i1,i2,i99]
884                        zop+=       xfac_op*a[i0,i1,i2,i99]
885
886                      else:
887                        zo+=a[i0,i1,i2,i99]
888                      z+=b[i0,i1,i2,i99]
889
890                 out+=tap+"%s[%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,zo,zo_fac_s,z*z_fac,z_fac_s)
891                 if zop==0.:
892                   out+="\n"
893                 else:
894                   out+="+(%s)*0.5**o\n"%zop
895
896        elif a.rank==5:
897            for i0 in range(a.shape[0]):
898             for i1 in range(a.shape[1]):
899               for i2 in range(a.shape[2]):
900                for i3 in range(a.shape[3]):
901                 zo=0.
902                 zop=0.
903                 z=0.
904                 for i99 in range(a.shape[4]):
905                      if i99==0:
906                        zo+=       xfac_o*a[i0,i1,i2,i3,i99]
907                        zop+=       xfac_op*a[i0,i1,i2,i3,i99]
908
909                      else:
910                        zo+=a[i0,i1,i2,i3,i99]
911                      z+=b[i0,i1,i2,i3,i99]
912                 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)
913                 if zop==0.:
914                   out+="\n"
915                 else:
916                   out+="+(%s)*0.5**o\n"%zop
917
918        return out
919
920  #=======================================================================================================  #=======================================================================================================
921  # basic binary operations (tests only!)  # interpolation
922  #=======================================================================================================  #=======================================================================================================
923  oper_range=[-5.,5.]  for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne","Solution","ReducedSolution"]:
924  for oper in [["add" ,"+",[-5.,5.]],    for data in ["Data","Symbol"]:
925               ["mult","*",[-5.,5.]],       for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
926               ["quotient" ,"/",[-5.,5.]],        for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
927               ["power" ,"**",[0.01,5.]]]:          if  where==case or \
928     for case0 in case_set:              ( case in ["ContinuousFunction","Solution","ReducedSolution"] and where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"] ) or \
929       for case1 in case_set:              ( case in ["FunctionOnContactZero","FunctionOnContactOne"] and where in ["FunctionOnContactZero","FunctionOnContactOne"] ) or \
930         for sh in shape_set:              (case=="ContinuousFunction" and  where in ["Solution","ReducedSolution"]) or \
931           for sh_p in shape_set:              (case=="Solution" and  where=="ReducedSolution") :
932             if len(sh_p)>0:
933                resource=[-1,1]
934             else:           text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
935                resource=[1]           tname="test_interpolation_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
936             for sh_d in resource:           text+="   def %s(self):\n"%tname
937              if sh_d>0:           text+="      \"\"\"\n"
938                 sh0=sh           text+="      tests interpolation for rank %s %s onto the %s\n\n"%(len(sh),data,where)
939                 sh1=sh+sh_p           text+="      assumptions: self.domain supports inpterpolation from %s onto %s\n"%(case,where)
940             text+="      \"\"\"\n"
941             if case=="ReducedSolution" or where=="ReducedSolution":
942                text+="      o=1\n"
943                o="1"
944             else:
945                text+="      o=self.order\n"
946                o="o"
947             text+="      dim=self.domain.getDim()\n"
948             text+="      w_ref=%s(self.domain)\n"%where
949             text+="      x_ref=w_ref.getX()\n"
950             text+="      w=%s(self.domain)\n"%case
951             text+="      x=w.getX()\n"
952             a_2=makeArray(sh+(2,),[-1.,1])
953             b_2=makeArray(sh+(2,),[-1.,1])
954             a_3=makeArray(sh+(3,),[-1.,1])
955             b_3=makeArray(sh+(3,),[-1.,1])
956             if data=="Symbol":
957                text+="      arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
958                val="s"
959                res="sub"
960             else:
961                val="arg"
962                res="res"
963             text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
964             text+="      ref=Data(0,%s,w_ref)\n"%str(sh)
965             text+="      if dim==2:\n"
966             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
967             text+=unrollLoops(a_2,b_2,o,"ref",tap="        ",x="x_ref")
968             text+="      else:\n"
969
970             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
971             text+=unrollLoops(a_3,b_3,o,"ref",tap="        ",x="x_ref")
972             text+="      res=interpolate(arg,where=w_ref)\n"
973             if data=="Symbol":
974                text+="      sub=res.substitute({arg:s})\n"
975             text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
976             text+="      self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong functionspace of result.\")\n"%res
977             text+="      self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
978             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
979             t_prog+=text
981    print t_prog
982    print test_tail
983    1/0
984
985    #=======================================================================================================
987    #=======================================================================================================
988    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
989      for data in ["Data","Symbol"]:
990         for case in ["ContinuousFunction","Solution","ReducedSolution"]:
991           for sh in [ (),(2,), (4,5), (6,2,2)]:
992             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
994             text+="   def %s(self):\n"%tname
995             text+="      \"\"\"\n"
996             text+="      tests gradient for rank %s %s on the %s\n\n"%(len(sh),data,where)
997             text+="      assumptions: %s(self.domain) exists\n"%case
998             text+="                   self.domain supports gardient on %s\n"%where
999             text+="      \"\"\"\n"
1000             if case=="ReducedSolution":
1001                text+="      o=1\n"
1002                o="1"
1003             else:
1004                text+="      o=self.order\n"
1005                o="o"
1006             text+="      dim=self.domain.getDim()\n"
1007             text+="      w_ref=%s(self.domain)\n"%where
1008             text+="      x_ref=w_ref.getX()\n"
1009             text+="      w=%s(self.domain)\n"%case
1010             text+="      x=w.getX()\n"
1011             a_2=makeArray(sh+(2,),[-1.,1])
1012             b_2=makeArray(sh+(2,),[-1.,1])
1013             a_3=makeArray(sh+(3,),[-1.,1])
1014             b_3=makeArray(sh+(3,),[-1.,1])
1015             if data=="Symbol":
1016                text+="      arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1017                val="s"
1018                res="sub"
1019             else:
1020                val="arg"
1021                res="res"
1022             text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
1023             text+="      ref=Data(0,%s+(dim,),w_ref)\n"%str(sh)
1024             text+="      if dim==2:\n"
1025             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1027             text+="      else:\n"
1028
1029             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1032             if data=="Symbol":
1033                text+="      sub=res.substitute({arg:s})\n"
1034             text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1035             text+="      self.failUnlessEqual(res.getShape(),%s+(dim,),\"wrong shape of result.\")\n"%str(sh)
1036             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1037
1038
1039             t_prog+=text
1041    print t_prog
1042    print test_tail
1043    1/0
1044
1045
1046    #=======================================================================================================
1047    # integrate
1048    #=======================================================================================================
1049    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1050      for data in ["Data","Symbol"]:
1051        for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1052          for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1053            if (not case in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]) or where==case:
1054             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1055             tname="test_integrate_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1056             text+="   def %s(self):\n"%tname
1057             text+="      \"\"\"\n"
1058             text+="      tests integral of rank %s %s on the %s\n\n"%(len(sh),data,where)
1059             text+="      assumptions: %s(self.domain) exists\n"%case
1060             text+="                   self.domain supports integral on %s\n"%where
1061
1062             text+="      \"\"\"\n"
1063             if case=="ReducedSolution":
1064                text+="      o=1\n"
1065                o="1"
1066             else:
1067                text+="      o=self.order\n"
1068                o="o"
1069             text+="      dim=self.domain.getDim()\n"
1070             text+="      w_ref=%s(self.domain)\n"%where
1071             text+="      w=%s(self.domain)\n"%case
1072             text+="      x=w.getX()\n"
1073             a_2=makeArray(sh+(2,),[-1.,1])
1074             b_2=makeArray(sh+(2,),[-1.,1])
1075             a_3=makeArray(sh+(3,),[-1.,1])
1076             b_3=makeArray(sh+(3,),[-1.,1])
1077             if data=="Symbol":
1078                text+="      arg=Symbol(shape=%s)\n"%str(sh)
1079                val="s"
1080                res="sub"
1081             else:
1082                val="arg"
1083                res="res"
1084
1085             text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
1086             if not len(sh)==0:
1087                text+="      ref=numarray.zeros(%s,numarray.Float)\n"%str(sh)
1088             text+="      if dim==2:\n"
1089             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1090             text+=unrollLoopsOfInteriorIntegral(a_2,b_2,where,"ref",tap="        ")
1091             text+="      else:\n"
1092
1093             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1094             text+=unrollLoopsOfInteriorIntegral(a_3,b_3,where,"ref",tap="        ")
1095             if case in ["ContinuousFunction","Solution","ReducedSolution"]:
1096                 text+="      res=integrate(arg,where=w_ref)\n"
1097             else:
1098                 text+="      res=integrate(arg)\n"
1099
1100             if data=="Symbol":
1101                text+="      sub=res.substitute({arg:s})\n"
1102             if len(sh)==0 and data=="Data":
1103                text+="      self.failUnless(isinstance(%s,float),\"wrong type of result.\")\n"%res
1104             else:
1105                if data=="Symbol":
1106                   text+="      self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1107                   text+="      self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1108              else:              else:
1109                 sh1=sh                 text+="      self.failUnless(isinstance(res,numarray.NumArray),\"wrong type of result.\")\n"
1110                 sh0=sh+sh_p                 text+="      self.failUnlessEqual(res.shape,%s,\"wrong shape of result.\")\n"%str(sh)
1111                         text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1112              if (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \
1113                 len(sh0)<5 and len(sh1)<5:
1114             t_prog+=text
1116    print t_prog
1117    print test_tail
1118    1/0
1119    #=======================================================================================================
1120    # inverse
1121    #=======================================================================================================
1122    name="inverse"
1123    for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1124      for sh0 in [ (1,1), (2,2), (3,3)]:
1125                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1126                  tname="test_%s_%s_dim%s"%(name,case0,sh0[0])
1127                  text+="   def %s(self):\n"%tname
1128                  a_0=makeArray(sh0,[-1.,1])
1129                  for i in range(sh0[0]): a_0[i,i]+=2
1130                  if case0 in ["taggedData", "expandedData"]:
1131                      a1_0=makeArray(sh0,[-1.,1])
1132                      for i in range(sh0[0]): a1_0[i,i]+=3
1133                  else:
1134                      a1_0=a_0
1135
1136                  text+=mkText(case0,"arg",a_0,a1_0)
1137                  text+="      res=%s(arg)\n"%name
1138                  if case0=="Symbol":
1139                     text+=mkText("array","s",a_0,a1_0)
1140                     text+="      sub=res.substitute({arg:s})\n"
1141                     res="sub"
1142                     ref="s"
1143                  else:
1144                     ref="arg"
1145                     res="res"
1146                  text+=mkTypeAndShapeTest(case0,sh0,"res")
1147                  text+="      self.failUnless(Lsup(matrixmult(%s,%s)-kronecker(%s))<=self.RES_TOL,\"wrong result\")\n"%(res,ref,sh0[0])
1148
1149                  if case0 == "taggedData" :
1150                      t_prog_with_tags+=text
1151                  else:
1152                      t_prog+=text
1153
1155    # print t_prog
1156    print t_prog_with_tags
1157    print test_tail
1158    1/0
1159
1160    #=======================================================================================================
1161    # trace
1162    #=======================================================================================================
1163    def traceTest(r,offset):
1164        sh=r.shape
1165        r1=1
1166        for i in range(offset): r1*=sh[i]
1167        r2=1
1168        for i in range(offset+2,len(sh)): r2*=sh[i]
1169        r_s=numarray.reshape(r,(r1,sh[offset],sh[offset],r2))
1170        s=numarray.zeros([r1,r2],numarray.Float)
1171        for i1 in range(r1):
1172            for i2 in range(r2):
1173                for j in range(sh[offset]): s[i1,i2]+=r_s[i1,j,j,i2]
1174        return s.resize(sh[:offset]+sh[offset+2:])
1175    name,tt="trace",traceTest
1176    for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1177      for sh0 in [ (4,5), (6,2,2),(3,2,3,4)]:
1178        for offset in range(len(sh0)-1):
1179                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1180                  tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1181                  text+="   def %s(self):\n"%tname
1182                  sh_t=list(sh0)
1183                  sh_t[offset+1]=sh_t[offset]
1184                  sh_t=tuple(sh_t)
1185                  sh_r=[]
1186                  for i in range(offset): sh_r.append(sh0[i])
1187                  for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])
1188                  sh_r=tuple(sh_r)
1189                  a_0=makeArray(sh_t,[-1.,1])
1190                  if case0 in ["taggedData", "expandedData"]:
1191                      a1_0=makeArray(sh_t,[-1.,1])
1192                  else:
1193                      a1_0=a_0
1194                  r=tt(a_0,offset)
1195                  r1=tt(a1_0,offset)
1196                  text+=mkText(case0,"arg",a_0,a1_0)
1197                  text+="      res=%s(arg,%s)\n"%(name,offset)
1198                  if case0=="Symbol":
1199                     text+=mkText("array","s",a_0,a1_0)
1200                     text+="      sub=res.substitute({arg:s})\n"
1201                     res="sub"
1202                     text+=mkText("array","ref",r,r1)
1203                  else:
1204                     res="res"
1205                     text+=mkText(case0,"ref",r,r1)
1206                  text+=mkTypeAndShapeTest(case0,sh_r,"res")
1207                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1208
1209                  if case0 == "taggedData" :
1210                      t_prog_with_tags+=text
1211                  else:
1212                      t_prog+=text
1213
1215    # print t_prog
1216    print t_prog_with_tags
1217    print test_tail
1218    1/0
1219
1220    #=======================================================================================================
1221    # clip
1222    #=======================================================================================================
1223    oper_L=[["clip",clipTEST]]
1224    for oper in oper_L:
1225     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1226      for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1227            if len(sh0)==0 or not case0=="float":
1228                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1229                  tname="test_%s_%s_rank%s"%(oper[0],case0,len(sh0))
1230                  text+="   def %s(self):\n"%tname
1231                  a_0=makeArray(sh0,[-1.,1])
1232                  if case0 in ["taggedData", "expandedData"]:
1233                      a1_0=makeArray(sh0,[-1.,1])
1234                  else:
1235                      a1_0=a_0
1236
1237                  r=oper[1](a_0,-0.3,0.5)
1238                  r1=oper[1](a1_0,-0.3,0.5)
1239                  text+=mkText(case0,"arg",a_0,a1_0)
1240                  text+="      res=%s(arg,-0.3,0.5)\n"%oper[0]
1241                  if case0=="Symbol":
1242                     text+=mkText("array","s",a_0,a1_0)
1243                     text+="      sub=res.substitute({arg:s})\n"
1244                     res="sub"
1245                     text+=mkText("array","ref",r,r1)
1246                  else:
1247                     res="res"
1248                     text+=mkText(case0,"ref",r,r1)
1249                  text+=mkTypeAndShapeTest(case0,sh0,"res")
1250                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1251
1252                  if case0 == "taggedData" :
1253                      t_prog_with_tags+=text
1254                  else:
1255                      t_prog+=text
1256
1258    # print t_prog
1259    print t_prog_with_tags
1260    print test_tail
1261    1/0
1262
1263    #=======================================================================================================
1264    # maximum, minimum, clipping
1265    #=======================================================================================================
1266    oper_L=[ ["maximum",maximumTEST],
1267             ["minimum",minimumTEST]]
1268    for oper in oper_L:
1269     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1270      for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1271       for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1272         for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1273            if (len(sh0)==0 or not case0=="float") and (len(sh1)==0 or not case1=="float") \
1274               and (sh0==sh1 or len(sh0)==0 or len(sh1)==0) :
1275                use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"                use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1276
1277                text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"                text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1278                tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))                tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1279                text+="   def %s(self):\n"%tname                text+="   def %s(self):\n"%tname
1280                a_0=makeArray(sh0,oper[2])                a_0=makeArray(sh0,[-1.,1])
1281                if case0 in ["taggedData", "expandedData"]:                if case0 in ["taggedData", "expandedData"]:
1282                    a1_0=makeArray(sh0,oper[2])                    a1_0=makeArray(sh0,[-1.,1])
1283                else:                else:
1284                    a1_0=a_0                    a1_0=a_0
1285
1286                a_1=makeArray(sh1,oper[2])                a_1=makeArray(sh1,[-1.,1])
1287                if case1 in ["taggedData", "expandedData"]:                if case1 in ["taggedData", "expandedData"]:
1288                    a1_1=makeArray(sh1,oper[2])                    a1_1=makeArray(sh1,[-1.,1])
1289                else:                else:
1290                    a1_1=a_1                    a1_1=a_1
1291                r1=makeResult2(a1_0,a1_1,"%a1%"+oper[1]+"%a2%")                r=oper[1](a_0,a_1)
1292                r=makeResult2(a_0,a_1,"%a1%"+oper[1]+"%a2%")                r1=oper[1](a1_0,a1_1)
1293                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)
1294                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)
1295                text+="      res=%s(arg0,arg1)\n"%oper[0]                text+="      res=%s(arg0,arg1)\n"%oper[0]
1296                  case=getResultCaseForBin(case0,case1)
1297                  if case=="Symbol":
1298                     c0_res,c1_res=case0,case1
1299                     subs="{"
1300                     if case0=="Symbol":
1301                        text+=mkText("array","s0",a_0,a1_0)
1302                        subs+="arg0:s0"
1303                        c0_res="array"
1304                     if case1=="Symbol":
1305                        text+=mkText("array","s1",a_1,a1_1)
1306                        if not subs.endswith("{"): subs+=","
1307                        subs+="arg1:s1"
1308                        c1_res="array"
1309                     subs+="}"
1310                     text+="      sub=res.substitute(%s)\n"%subs
1311                     res="sub"
1312                     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1313                  else:
1314                     res="res"
1315                     text+=mkText(case,"ref",r,r1)
1316                  if len(sh0)>len(sh1):
1317                      text+=mkTypeAndShapeTest(case,sh0,"res")
1318                  else:
1319                      text+=mkTypeAndShapeTest(case,sh1,"res")
1320                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1321
1322                  if case0 == "taggedData" or case1 == "taggedData":
1323                      t_prog_with_tags+=text
1324                  else:
1325                      t_prog+=text
1326
1328    # print t_prog
1329    print t_prog_with_tags
1330    print test_tail
1331    1/0
1332
1333
1334    #=======================================================================================================
1335    # outer inner
1336    #=======================================================================================================
1337    oper=["outer",outerTEST]
1338    # oper=["inner",innerTEST]
1339    for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1340      for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1341       for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1342         for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1343            if (len(sh0)==0 or not case0=="float") and (len(sh1)==0 or not case1=="float") \
1344               and len(sh0+sh1)<5:
1345                  use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1346
1347                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1348                  tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1349                  text+="   def %s(self):\n"%tname
1350                  a_0=makeArray(sh0,[-1.,1])
1351                  if case0 in ["taggedData", "expandedData"]:
1352                      a1_0=makeArray(sh0,[-1.,1])
1353                  else:
1354                      a1_0=a_0
1355
1356                  a_1=makeArray(sh1,[-1.,1])
1357                  if case1 in ["taggedData", "expandedData"]:
1358                      a1_1=makeArray(sh1,[-1.,1])
1359                  else:
1360                      a1_1=a_1
1361                  r=oper[1](a_0,a_1)
1362                  r1=oper[1](a1_0,a1_1)
1363                  text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1364                  text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1365                  text+="      res=%s(arg0,arg1)\n"%oper[0]
1366                case=getResultCaseForBin(case0,case1)                              case=getResultCaseForBin(case0,case1)
1367                if case=="Symbol":                if case=="Symbol":
1368                   c0_res,c1_res=case0,case1                   c0_res,c1_res=case0,case1
# Line 584  for oper in [["add" ,"+",[-5.,5.]], Line 1383  for oper in [["add" ,"+",[-5.,5.]],
1383                else:                else:
1384                   res="res"                   res="res"
1385                   text+=mkText(case,"ref",r,r1)                   text+=mkText(case,"ref",r,r1)
1386                if isinstance(r,float):                              text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
text+=mkTypeAndShapeTest(case,(),"res")
else:
text+=mkTypeAndShapeTest(case,r.shape,"res")
1387                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
1388
1389                if case0 == "taggedData" or case1 == "taggedData":                if case0 == "taggedData" or case1 == "taggedData":
1390                    t_prog_with_tags+=text                    t_prog_with_tags+=text
1391                else:                              else:
1392                    t_prog+=text                    t_prog+=text
1393
1395  # print t_prog  # print t_prog
1396  print t_prog_with_tags  print t_prog_with_tags
1397  print test_tail  print test_tail
1398  1/0  1/0
1399
1400  #=======================================================================================================  #=======================================================================================================
1401    # local reduction
1402    #=======================================================================================================
1403    for oper in [["length",0.,"out+%a1%**2","math.sqrt(out)"],
1404                 ["maxval",-1.e99,"max(out,%a1%)","out"],
1405                 ["minval",1.e99,"min(out,%a1%)","out"] ]:
1406      for case in case_set:
1407         for sh in shape_set:
1408           if not case=="float" or len(sh)==0:
1409             text=""
1410             text+="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1411             tname="def test_%s_%s_rank%s"%(oper[0],case,len(sh))
1412             text+="   %s(self):\n"%tname
1413             a=makeArray(sh,[-1.,1.])
1414             a1=makeArray(sh,[-1.,1.])
1415             r1=testReduce(a1,oper[1],oper[2],oper[3])
1416             r=testReduce(a,oper[1],oper[2],oper[3])
1417
1418             text+=mkText(case,"arg",a,a1)
1419             text+="      res=%s(arg)\n"%oper[0]
1420             if case=="Symbol":
1421                 text+=mkText("array","s",a,a1)
1422                 text+="      sub=res.substitute({arg:s})\n"
1423                 text+=mkText("array","ref",r,r1)
1424                 res="sub"
1425             else:
1426                 text+=mkText(case,"ref",r,r1)
1427                 res="res"
1428             if oper[0]=="length":
1429                   text+=mkTypeAndShapeTest(case,(),"res")
1430             else:
1431                if case=="float" or case=="array":
1432                   text+=mkTypeAndShapeTest("float",(),"res")
1433                else:
1434                   text+=mkTypeAndShapeTest(case,(),"res")
1435             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1436             if case == "taggedData":
1437               t_prog_with_tags+=text
1438             else:
1439               t_prog+=text
1441    # print t_prog
1442    print t_prog_with_tags
1443    print test_tail
1444    1/0
1445
1446    #=======================================================================================================
1447    # tensor multiply
1448    #=======================================================================================================
1449    # oper=["generalTensorProduct",tensorProductTest]
1450    # oper=["matrixmult",testMatrixMult]
1451    oper=["tensormult",testTensorMult]
1452
1453    for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1454      for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1455       for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1456         for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1457           for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:
1458              if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \
1459                   and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
1460                # 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
1461                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
1462                  case=getResultCaseForBin(case0,case1)
1463                  use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1464                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1465                  # tname="test_generalTensorProduct_%s_rank%s_%s_rank%s_offset%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1),len(sh_s))
1466                  #tname="test_matrixmult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
1467                  tname="test_tensormult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
1468                  # if tname=="test_generalTensorProduct_array_rank1_array_rank2_offset1":
1469                  # print tnametest_generalTensorProduct_Symbol_rank1_Symbol_rank3_offset1
1470                  text+="   def %s(self):\n"%tname
1471                  a_0=makeArray(sh0+sh_s,[-1.,1])
1472                  if case0 in ["taggedData", "expandedData"]:
1473                      a1_0=makeArray(sh0+sh_s,[-1.,1])
1474                  else:
1475                      a1_0=a_0
1476
1477                  a_1=makeArray(sh_s+sh1,[-1.,1])
1478                  if case1 in ["taggedData", "expandedData"]:
1479                      a1_1=makeArray(sh_s+sh1,[-1.,1])
1480                  else:
1481                      a1_1=a_1
1482                  r=oper[1](a_0,a_1,sh_s)
1483                  r1=oper[1](a1_0,a1_1,sh_s)
1484                  text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1485                  text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1486                  #text+="      res=matrixmult(arg0,arg1)\n"
1487                  text+="      res=tensormult(arg0,arg1)\n"
1488                  #text+="      res=generalTensorProduct(arg0,arg1,offset=%s)\n"%(len(sh_s))
1489                  if case=="Symbol":
1490                     c0_res,c1_res=case0,case1
1491                     subs="{"
1492                     if case0=="Symbol":
1493                        text+=mkText("array","s0",a_0,a1_0)
1494                        subs+="arg0:s0"
1495                        c0_res="array"
1496                     if case1=="Symbol":
1497                        text+=mkText("array","s1",a_1,a1_1)
1498                        if not subs.endswith("{"): subs+=","
1499                        subs+="arg1:s1"
1500                        c1_res="array"
1501                     subs+="}"
1502                     text+="      sub=res.substitute(%s)\n"%subs
1503                     res="sub"
1504                     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1505                  else:
1506                     res="res"
1507                     text+=mkText(case,"ref",r,r1)
1508                  text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
1509                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1510                  if case0 == "taggedData" or case1 == "taggedData":
1511                      t_prog_with_tags+=text
1512                  else:
1513                      t_prog+=text
1515    # print t_prog
1516    print t_prog_with_tags
1517    print test_tail
1518    1/0
1519    #=======================================================================================================
1520  # basic binary operation overloading (tests only!)  # basic binary operation overloading (tests only!)
1521  #=======================================================================================================  #=======================================================================================================
1522  oper_range=[-5.,5.]  oper_range=[-5.,5.]
# Line 613  for oper in [["add" ,"+",[-5.,5.]], Line 1529  for oper in [["add" ,"+",[-5.,5.]],
1529       for sh0 in shape_set:       for sh0 in shape_set:
1530         for case1 in case_set:         for case1 in case_set:
1531           for sh1 in shape_set:           for sh1 in shape_set:
1532             if (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \             if not case0=="array" and \
1533                   (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \
1534                 (sh0==() or sh1==() or sh1==sh0) and \                 (sh0==() or sh1==() or sh1==sh0) and \
1535                 not (case0 in ["float","array"] and  case1 in ["float","array"]):                 not (case0 in ["float","array"] and  case1 in ["float","array"]):
1536                  use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1537                text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"                text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1539                text+="   def %s(self):\n"%tname                text+="   def %s(self):\n"%tname
# Line 671  for oper in [["add" ,"+",[-5.,5.]], Line 1589  for oper in [["add" ,"+",[-5.,5.]],
1589                      t_prog+=text                      t_prog+=text
1590
1591
# print u_prog
# 1/0
1593  print t_prog  # print t_prog
1594    # print t_prog_with_tags
1595    print t_prog_failing
1596    print test_tail
1597    1/0
1598    #=======================================================================================================
1599    # basic binary operations (tests only!)
1600    #=======================================================================================================
1601    oper_range=[-5.,5.]
1602    for oper in [["add" ,"+",[-5.,5.]],
1603                 ["mult","*",[-5.,5.]],
1604                 ["quotient" ,"/",[-5.,5.]],
1605                 ["power" ,"**",[0.01,5.]]]:
1606       for case0 in case_set:
1607         for case1 in case_set:
1608           for sh in shape_set:
1609             for sh_p in shape_set:
1610               if len(sh_p)>0:
1611                  resource=[-1,1]
1612               else:
1613                  resource=[1]
1614               for sh_d in resource:
1615                if sh_d>0:
1616                   sh0=sh
1617                   sh1=sh+sh_p
1618                else:
1619                   sh1=sh
1620                   sh0=sh+sh_p
1621
1622                if (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \
1623                   len(sh0)<5 and len(sh1)<5:
1624                  use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1625                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1626                  tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1627                  text+="   def %s(self):\n"%tname
1628                  a_0=makeArray(sh0,oper[2])
1629                  if case0 in ["taggedData", "expandedData"]:
1630                      a1_0=makeArray(sh0,oper[2])
1631                  else:
1632                      a1_0=a_0
1633
1634                  a_1=makeArray(sh1,oper[2])
1635                  if case1 in ["taggedData", "expandedData"]:
1636                      a1_1=makeArray(sh1,oper[2])
1637                  else:
1638                      a1_1=a_1
1639                  r1=makeResult2(a1_0,a1_1,"%a1%"+oper[1]+"%a2%")
1640                  r=makeResult2(a_0,a_1,"%a1%"+oper[1]+"%a2%")
1641                  text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1642                  text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1643                  text+="      res=%s(arg0,arg1)\n"%oper[0]
1644
1645                  case=getResultCaseForBin(case0,case1)
1646                  if case=="Symbol":
1647                     c0_res,c1_res=case0,case1
1648                     subs="{"
1649                     if case0=="Symbol":
1650                        text+=mkText("array","s0",a_0,a1_0)
1651                        subs+="arg0:s0"
1652                        c0_res="array"
1653                     if case1=="Symbol":
1654                        text+=mkText("array","s1",a_1,a1_1)
1655                        if not subs.endswith("{"): subs+=","
1656                        subs+="arg1:s1"
1657                        c1_res="array"
1658                     subs+="}"
1659                     text+="      sub=res.substitute(%s)\n"%subs
1660                     res="sub"
1661                     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1662                  else:
1663                     res="res"
1664                     text+=mkText(case,"ref",r,r1)
1665                  if isinstance(r,float):
1666                     text+=mkTypeAndShapeTest(case,(),"res")
1667                  else:
1668                     text+=mkTypeAndShapeTest(case,r.shape,"res")
1669                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1670
1671                  if case0 == "taggedData" or case1 == "taggedData":
1672                      t_prog_with_tags+=text
1673                  else:
1674                      t_prog+=text
1676    # print t_prog
1677    print t_prog_with_tags
1678    print test_tail
1679  1/0  1/0
1680
1681  # print t_prog_with_tagsoper_range=[-5.,5.]  # print t_prog_with_tagsoper_range=[-5.,5.]
1682  for oper in [["add" ,"+",[-5.,5.]],  for oper in [["add" ,"+",[-5.,5.]],
1683               ["sub" ,"-",[-5.,5.]],               ["sub" ,"-",[-5.,5.]],
# Line 1375  def X(): Line 2377  def X():
2377                                ref_diff=(makeResult(trafo[j0,j1,j2,j3]*a_in+finc,f)-makeResult(trafo[j0,j1,j2,j3]*a_in,f))/finc                                ref_diff=(makeResult(trafo[j0,j1,j2,j3]*a_in+finc,f)-makeResult(trafo[j0,j1,j2,j3]*a_in,f))/finc
2378                                t_prog+="      self.failUnlessAlmostEqual(dvdin[%s,%s,%s,%s],%s,self.places,\"%s-derivative: wrong derivative of %s\")\n"%(j0,j1,j2,j3,ref_diff,str(sh_in),str((j0,j1,j2,j3)))                                t_prog+="      self.failUnlessAlmostEqual(dvdin[%s,%s,%s,%s],%s,self.places,\"%s-derivative: wrong derivative of %s\")\n"%(j0,j1,j2,j3,ref_diff,str(sh_in),str((j0,j1,j2,j3)))
2379
2380  #==================  #
case="inner"
for arg0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
for arg1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
sh1=sh0
if (len(sh0)==0 or not arg0=="float") and (len(sh1)==0 or not arg1=="float"):
tname="test_%s_%s_rank%s_%s_rank%s"%(case,arg0,len(sh0),arg1,len(sh1))
t_prog+="   def %s(self):\n"%tname
a0=makeArray(sh0,[-1,1])
a0_1=makeArray(sh0,[-1,1])
a1=makeArray(sh1,[-1,1])
a1_1=makeArray(sh1,[-1,1])
t_prog+=mkText(arg0,"arg0",a0,a0_1)
t_prog+=mkText(arg1,"arg1",a1,a1_1)
t_prog+="      res=%s(arg0,arg1)\n"%case

print t_prog
1/0
2381
2382  #==================  #==================
2383  cases=["Scalar","Vector","Tensor", "Tensor3","Tensor4"]  cases=["Scalar","Vector","Tensor", "Tensor3","Tensor4"]

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