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

Diff of /trunk/escript/py_src/generatediff

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

revision 439 by gross, Fri Jan 20 01:46:22 2006 UTC revision 587 by gross, Fri Mar 10 02:26:50 2006 UTC
# Line 179  def makeArray(shape,rng): Line 179  def makeArray(shape,rng):
179               for i2 in range(shape[2]):               for i2 in range(shape[2]):
180                  for i3 in range(shape[3]):                  for i3 in range(shape[3]):
181                    for i4 in range(shape[4]):                    for i4 in range(shape[4]):
182                     out[i0,i1,i2,i3,i4]=l*random.random()+rng[0]                     out[i0,i1,i2,i3,i4]=l*ranm.random()+rng[0]
183     else:     else:
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*int(8*random.random()+1)
194       elif len(shape)==2:
195           for i0 in range(shape[0]):
196              for i1 in range(shape[1]):
197                       out[i0,i1]=s*int(8*random.random()+1)
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*int(8*random.random()+1)
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*int(8*random.random()+1)
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 480  def mkText(case,name,a,a1=None,use_taggi Line 505  def mkText(case,name,a,a1=None,use_taggi
505                   t_out+="      %s.setTaggedValue(1,numarray.array(%s))\n"%(name,a1.tolist())                   t_out+="      %s.setTaggedValue(1,numarray.array(%s))\n"%(name,a1.tolist())
506                t_out+="      %s.expand()\n"%name                          t_out+="      %s.expand()\n"%name          
507             else:             else:
508                t_out+="      msk_%s=whereNegative(self.functionspace.getX()[0]-0.5)\n"%name                t_out+="      msk_%s=whereZero(self.functionspace.getX()[0],1.e-8)\n"%name
509                if isinstance(a,float):                if isinstance(a,float):
510                     t_out+="      %s=msk_%s*(%s)+(1.-msk_%s)*(%s)\n"%(name,name,a,name,a1)                     t_out+="      %s=(1.-msk_%s)*(%s)+msk_%s*(%s)\n"%(name,name,a,name,a1)
511                elif a.rank==0:                elif a.rank==0:
512                     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)
513                else:                else:
# Line 497  def mkText(case,name,a,a1=None,use_taggi Line 522  def mkText(case,name,a,a1=None,use_taggi
522    
523           return t_out           return t_out
524    
525  def mkTypeAndShapeTest(case,sh,argstr):  def mkTypeAndShapeTest(case,sh,argstr,name=""):
526      text=""      text=""
527      if case=="float":      if case=="float":
528           text+="      self.failUnless(isinstance(%s,float),\"wrong type of result.\")\n"%argstr           text+="      self.failUnless(isinstance(%s,float),\"wrong type of result%s.\")\n"%(argstr,name)
529      elif case=="array":      elif case=="array":
530           text+="      self.failUnless(isinstance(%s,numarray.NumArray),\"wrong type of result.\")\n"%argstr           text+="      self.failUnless(isinstance(%s,numarray.NumArray),\"wrong type of result%s.\")\n"%(argstr,name)
531           text+="      self.failUnlessEqual(%s.shape,%s,\"wrong shape of result.\")\n"%(argstr,str(sh))           text+="      self.failUnlessEqual(%s.shape,%s,\"wrong shape of result%s.\")\n"%(argstr,str(sh),name)
532      elif case in ["constData","taggedData","expandedData"]:          elif case in ["constData","taggedData","expandedData"]:    
533           text+="      self.failUnless(isinstance(%s,Data),\"wrong type of result.\")\n"%argstr           text+="      self.failUnless(isinstance(%s,Data),\"wrong type of result%s.\")\n"%(argstr,name)
534           text+="      self.failUnlessEqual(%s.getShape(),%s,\"wrong shape of result.\")\n"%(argstr,str(sh))           text+="      self.failUnlessEqual(%s.getShape(),%s,\"wrong shape of result%s.\")\n"%(argstr,str(sh),name)
535      elif case=="Symbol":      elif case=="Symbol":
536           text+="      self.failUnless(isinstance(%s,Symbol),\"wrong type of result.\")\n"%argstr           text+="      self.failUnless(isinstance(%s,Symbol),\"wrong type of result%s.\")\n"%(argstr,name)
537           text+="      self.failUnlessEqual(%s.getShape(),%s,\"wrong shape of result.\")\n"%(argstr,str(sh))           text+="      self.failUnlessEqual(%s.getShape(),%s,\"wrong shape of result%s.\")\n"%(argstr,str(sh),name)
538      return text      return text
539    
540  def mkCode(txt,args=[],intend=""):  def mkCode(txt,args=[],intend=""):
# Line 525  def mkCode(txt,args=[],intend=""): Line 550  def mkCode(txt,args=[],intend=""):
550        out=out.replace("%%a%s%%"%c,r)        out=out.replace("%%a%s%%"%c,r)
551      return out        return out  
552    
553    
554    #=======================================================================================================
555    # eigenvalues and eigen vectors 2D:
556    #=======================================================================================================
557    alltests= \
558      [ ("case0",[[0.0, 0.0], [0.0, 0.0]],[0.0, 0.0]) \
559       , ("case3",[[-1.0, 0.0], [0.0, -1.0]],[-1.0, -1.0]) \
560       , ("case5",[[-0.99999999999999967, -6.4606252205695602e-16], [-6.4606252205695602e-16, -0.99999999999999967]],[-1.0, -1.0]) \
561       , ("case6",[[0.0, 0.0], [0.0, 0.0001]],[0.0, 0.0001]) \
562       , ("case7",[[0.0001, 0.0], [0.0, 0.0]],[0.0, 0.0001]) \
563       , ("case8",[[6.0598371831785722e-06, 2.3859213977648625e-05], [2.3859213977648629e-05, 9.3940162816821425e-05]],[0.0, 0.0001]) \
564       , ("case9",[[1.0, 0.0], [0.0, 2.0]],[1.0, 2.0]) \
565       , ("case10",[[2.0, 0.0], [0.0, 1.0]],[1.0, 2.0]) \
566       , ("case11",[[1.0605983718317855, 0.23859213977648688], [0.23859213977648688, 1.9394016281682138]],[1.0, 2.0]) \
567       , ("case12",[[1.0, 0.0], [0.0, 1000000.0]],[1.0, 1000000.0]) \
568       , ("case13",[[1000000.0, 0.0], [0.0, 1.0]],[1.0, 1000000.0]) \
569       , ("case14",[[60599.311233413886, 238591.90118434647], [238591.90118434647, 939401.68876658613]],[1.0, 1000000.0]) \
570       ]
571    dim=2
572    for case in ["constData","taggedData","expandedData"]:
573       n=0
574       while n<len(alltests):
575          text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
576          tname="test_eigenvalues_and_eigenvectors_%s_dim%s_%s"%(case,dim,alltests[n][0])
577          text+="   def %s(self):\n"%tname
578          a_0=numarray.array(alltests[n][1],numarray.Float64)
579          ev_0=numarray.array(alltests[n][2],numarray.Float64)
580          if case in ["taggedData", "expandedData"]:
581             a1_0=numarray.array(alltests[n+1][1],numarray.Float64)
582             ev1_0=numarray.array(alltests[n+1][2],numarray.Float64)
583             n+=2
584          else:
585             a1_0=a_0                  
586             ev1_0=ev_0
587             n+=1
588          text+=mkText(case,"arg",a_0,a1_0)
589          text+="      res=eigenvalues_and_eigenvectors(arg)\n"
590          text+=mkText(case,"ref_ev",ev_0,ev1_0)
591          text+=mkTypeAndShapeTest(case,(dim,),"res[0]"," for eigenvalues")
592          text+=mkTypeAndShapeTest(case,(dim,dim),"res[1]"," for eigenvectors")
593          text+="      self.failUnless(Lsup(res[0]-ref_ev)<=self.RES_TOL*Lsup(ref_ev),\"wrong eigenvalues\")\n"
594          for i in range(dim):
595              text+="      self.failUnless(Lsup(matrixmult(arg,res[1][:,%s])-res[0][%s]*res[1][:,%s])<=self.RES_TOL*Lsup(res[0]),\"wrong eigenvector %s\")\n"%(i,i,i,i)
596          text+="      self.failUnless(Lsup(matrixmult(transpose(res[1]),res[1])-kronecker(%s))<=self.RES_TOL,\"eigenvectors are not orthonormal\")\n"%dim
597          if case == "taggedData" :
598               t_prog_with_tags+=text
599          else:              
600               t_prog+=text
601    print test_header
602    print t_prog
603    print t_prog_with_tags
604    print test_tail          
605    1/0
606    #=======================================================================================================
607    # get slices
608    #=======================================================================================================
609    from esys.escript import *
610    for case0 in ["constData","taggedData","expandedData"]:
611       for sh0 in [ (3,), (3,4), (3,4,3) ,(4,3,3,3)]:
612        # get perm:
613        if len(sh0)==2:
614            check=[[1,0]]
615        elif len(sh0)==3:
616            check=[[1,0,2],
617                   [1,2,0],
618                   [2,1,0],
619                   [2,0,2],
620                   [0,2,1]]
621        elif len(sh0)==4:
622            check=[[0,1,3,2],
623                   [0,2,1,3],
624                   [0,2,3,1],
625                   [0,3,2,1],
626                   [0,3,1,2] ,          
627                   [1,0,2,3],
628                   [1,0,3,2],
629                   [1,2,0,3],
630                   [1,2,3,0],
631                   [1,3,2,0],
632                   [1,3,0,2],
633                   [2,0,1,3],
634                   [2,0,3,1],
635                   [2,1,0,3],
636                   [2,1,3,0],
637                   [2,3,1,0],
638                   [2,3,0,1],
639                   [3,0,1,2],
640                   [3,0,2,1],
641                   [3,1,0,2],
642                   [3,1,2,0],
643                   [3,2,1,0],
644                   [3,2,0,1]]
645        else:
646             check=[]
647        
648        # create the test cases:
649        processed=[]
650        l=["R","U","L","P","C","N"]
651        c=[""]
652        for i in range(len(sh0)):
653           tmp=[]
654           for ci in c:
655              tmp+=[ci+li for li in l]
656           c=tmp
657        # SHUFFLE
658        c2=[]
659        while len(c)>0:
660            i=int(random.random()*len(c))
661            c2.append(c[i])
662            del c[i]
663        c=c2
664        for ci in c:
665          t=""
666          sh=()
667          sl=()
668          for i in range(len(ci)):
669              if ci[i]=="R":
670                 s="%s:%s"%(1,sh0[i]-1)
671                 sl=sl+(slice(1,sh0[i]-1),)
672                 sh=sh+(sh0[i]-2,)            
673              if ci[i]=="U":
674                  s=":%s"%(sh0[i]-1)
675                  sh=sh+(sh0[i]-1,)
676                  sl=sl+(slice(0,sh0[i]-1),)            
677              if ci[i]=="L":
678                  s="2:"
679                  sh=sh+(sh0[i]-2,)            
680                  sl=sl+(slice(2,sh0[i]),)            
681              if ci[i]=="P":
682                  s="%s"%(int(sh0[i]/2))
683                  sl=sl+(int(sh0[i]/2),)            
684              if ci[i]=="C":
685                  s=":"
686                  sh=sh+(sh0[i],)            
687                  sl=sl+(slice(0,sh0[i]),)            
688              if ci[i]=="N":
689                  s=""
690                  sh=sh+(sh0[i],)
691              if len(s)>0:
692                 if not t=="": t+=","
693                 t+=s
694          if len(sl)==1: sl=sl[0]
695          N_found=False
696          noN_found=False
697          process=len(t)>0
698          for i in ci:
699             if i=="N":
700                if not noN_found and N_found: process=False
701                N_found=True
702             else:
703               if N_found: process=False
704               noNfound=True
705          # is there a similar one processed allready
706          if process and ci.find("N")==-1:
707             for ci2 in processed:
708               for chi in check:
709                   is_perm=True
710                   for i in range(len(chi)):
711                       if not ci[i]==ci2[chi[i]]: is_perm=False
712                   if is_perm: process=False
713          # if not process: print ci," rejected"
714          if process:
715           processed.append(ci)
716           for case1 in ["array","constData","taggedData","expandedData"]:
717              text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
718              tname="test_setslice_%s_rank%s_%s_%s"%(case0,len(sh0),case1,ci)
719              text+="   def %s(self):\n"%tname
720              a_0=makeNumberedArray(sh0)
721              if case0 in ["taggedData", "expandedData"]:
722                  a1_0=makeNumberedArray(sh0)
723              else:
724                  a1_0=a_0*1.
725    
726              a_1=makeNumberedArray(sh)
727              if case1 in ["taggedData", "expandedData"]:
728                  a1_1=makeNumberedArray(sh)
729              else:
730                  a1_1=a_1*1.
731    
732              text+=mkText(case0,"arg",a_0,a1_0)                                  
733              text+=mkText(case1,"val",a_1,a1_1)                                  
734              text+="      arg[%s]=val\n"%t
735              a_0.__setitem__(sl,a_1)
736              a1_0.__setitem__(sl,a1_1)
737              if Lsup(a_0-a1_0)==0.:
738                 text+=mkText("constData","ref",a_0,a1_0)
739              else:
740                 text+=mkText("expandedData","ref",a_0,a1_0)
741              text+="      self.failUnless(Lsup(arg-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"
742                  
743              if case0 == "taggedData" or case1 == "taggedData":
744                t_prog_with_tags+=text
745              else:              
746                t_prog+=text
747    
748    print test_header
749    # print t_prog
750    print t_prog_with_tags
751    print test_tail          
752    1/0
753    
754    #=======================================================================================================
755    # (non)symmetric part
756    #=======================================================================================================
757    from esys.escript import *
758    for name in ["symmetric", "nonsymmetric"]:
759     f=1.
760     if name=="nonsymmetric": f=-1
761     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
762      for sh0 in [ (3,3), (2,3,2,3)]:
763                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
764                  tname="test_%s_%s_rank%s"%(name,case0,len(sh0))
765                  text+="   def %s(self):\n"%tname
766                  a_0=makeNumberedArray(sh0,s=1.)
767                  r_0=(a_0+f*transpose(a_0))/2.
768                  if case0 in ["taggedData", "expandedData"]:
769                     a1_0=makeNumberedArray(sh0,s=-1.)
770                     r1_0=(a1_0+f*transpose(a1_0))/2.
771                  else:
772                      a1_0=a_0                  
773                      r1_0=r_0
774                  text+=mkText(case0,"arg",a_0,a1_0)
775                  text+="      res=%s(arg)\n"%name
776                  if case0=="Symbol":
777                     text+=mkText("array","s",a_0,a1_0)
778                     text+="      sub=res.substitute({arg:s})\n"
779                     res="sub"
780                     text+=mkText("array","ref",r_0,r1_0)
781                  else:
782                     res="res"
783                     text+=mkText(case0,"ref",r_0,r1_0)  
784                  text+=mkTypeAndShapeTest(case0,sh0,"res")
785                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
786                  
787                  if case0 == "taggedData" :
788                      t_prog_with_tags+=text
789                  else:              
790                      t_prog+=text
791    print test_header
792    print t_prog
793    # print t_prog_with_tags
794    print test_tail          
795    1/0
796    
797    #=======================================================================================================
798    # eigenvalues
799    #=======================================================================================================
800    import numarray.linear_algebra
801    name="eigenvalues"
802    for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
803      for sh0 in [ (1,1), (2,2), (3,3)]:
804                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
805                  tname="test_%s_%s_dim%s"%(name,case0,sh0[0])
806                  text+="   def %s(self):\n"%tname
807                  a_0=makeArray(sh0,[-1.,1])
808                  a_0=(a_0+numarray.transpose(a_0))/2.
809                  ev=numarray.linear_algebra.eigenvalues(a_0)
810                  ev.sort()
811                  if case0 in ["taggedData", "expandedData"]:
812                      a1_0=makeArray(sh0,[-1.,1])
813                      a1_0=(a1_0+numarray.transpose(a1_0))/2.
814                      ev1=numarray.linear_algebra.eigenvalues(a1_0)
815                      ev1.sort()
816                  else:
817                      a1_0=a_0                  
818                      ev1=ev
819                  text+=mkText(case0,"arg",a_0,a1_0)
820                  text+="      res=%s(arg)\n"%name
821                  if case0=="Symbol":
822                     text+=mkText("array","s",a_0,a1_0)
823                     text+="      sub=res.substitute({arg:s})\n"
824                     res="sub"
825                     text+=mkText("array","ref",ev,ev1)
826                  else:
827                     res="res"
828                     text+=mkText(case0,"ref",ev,ev1)  
829                  text+=mkTypeAndShapeTest(case0,(sh0[0],),"res")
830                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
831                  
832                  if case0 == "taggedData" :
833                      t_prog_with_tags+=text
834                  else:              
835                      t_prog+=text
836    print test_header
837    # print t_prog
838    print t_prog_with_tags
839    print test_tail          
840    1/0
841    
842    #=======================================================================================================
843    # get slices
844    #=======================================================================================================
845    for case0 in ["constData","taggedData","expandedData","Symbol"]:
846      for sh0 in [ (3,), (3,4), (3,4,3) ,(4,3,5,3)]:
847        # get perm:
848        if len(sh0)==2:
849            check=[[1,0]]
850        elif len(sh0)==3:
851            check=[[1,0,2],
852                   [1,2,0],
853                   [2,1,0],
854                   [2,0,2],
855                   [0,2,1]]
856        elif len(sh0)==4:
857            check=[[0,1,3,2],
858                   [0,2,1,3],
859                   [0,2,3,1],
860                   [0,3,2,1],
861                   [0,3,1,2] ,          
862                   [1,0,2,3],
863                   [1,0,3,2],
864                   [1,2,0,3],
865                   [1,2,3,0],
866                   [1,3,2,0],
867                   [1,3,0,2],
868                   [2,0,1,3],
869                   [2,0,3,1],
870                   [2,1,0,3],
871                   [2,1,3,0],
872                   [2,3,1,0],
873                   [2,3,0,1],
874                   [3,0,1,2],
875                   [3,0,2,1],
876                   [3,1,0,2],
877                   [3,1,2,0],
878                   [3,2,1,0],
879                   [3,2,0,1]]
880        else:
881             check=[]
882        
883        # create the test cases:
884        processed=[]
885        l=["R","U","L","P","C","N"]
886        c=[""]
887        for i in range(len(sh0)):
888           tmp=[]
889           for ci in c:
890              tmp+=[ci+li for li in l]
891           c=tmp
892        # SHUFFLE
893        c2=[]
894        while len(c)>0:
895            i=int(random.random()*len(c))
896            c2.append(c[i])
897            del c[i]
898        c=c2
899        for ci in c:
900          t=""
901          sh=()
902          for i in range(len(ci)):
903              if ci[i]=="R":
904                 s="%s:%s"%(1,sh0[i]-1)
905                 sh=sh+(sh0[i]-2,)            
906              if ci[i]=="U":
907                  s=":%s"%(sh0[i]-1)
908                  sh=sh+(sh0[i]-1,)            
909              if ci[i]=="L":
910                  s="2:"
911                  sh=sh+(sh0[i]-2,)            
912              if ci[i]=="P":
913                  s="%s"%(int(sh0[i]/2))
914              if ci[i]=="C":
915                  s=":"
916                  sh=sh+(sh0[i],)            
917              if ci[i]=="N":
918                  s=""
919                  sh=sh+(sh0[i],)
920              if len(s)>0:
921                 if not t=="": t+=","
922                 t+=s
923          N_found=False
924          noN_found=False
925          process=len(t)>0
926          for i in ci:
927             if i=="N":
928                if not noN_found and N_found: process=False
929                N_found=True
930             else:
931               if N_found: process=False
932               noNfound=True
933          # is there a similar one processed allready
934          if process and ci.find("N")==-1:
935             for ci2 in processed:
936               for chi in check:
937                   is_perm=True
938                   for i in range(len(chi)):
939                       if not ci[i]==ci2[chi[i]]: is_perm=False
940                   if is_perm: process=False
941          # if not process: print ci," rejected"
942          if process:
943           processed.append(ci)
944           text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
945           tname="test_getslice_%s_rank%s_%s"%(case0,len(sh0),ci)
946           text+="   def %s(self):\n"%tname
947           a_0=makeNumberedArray(sh0,s=1)
948           if case0 in ["taggedData", "expandedData"]:
949                a1_0=makeNumberedArray(sh0,s=-1.)
950           else:
951                a1_0=a_0
952           r=eval("a_0[%s]"%t)
953           r1=eval("a1_0[%s]"%t)
954           text+=mkText(case0,"arg",a_0,a1_0)
955           text+="      res=arg[%s]\n"%t
956           if case0=="Symbol":
957               text+=mkText("array","s",a_0,a1_0)
958               text+="      sub=res.substitute({arg:s})\n"
959               res="sub"
960               text+=mkText("array","ref",r,r1)
961           else:
962               res="res"
963               text+=mkText(case0,"ref",r,r1)
964           text+=mkTypeAndShapeTest(case0,sh,"res")
965           text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
966                  
967           if case0 == "taggedData" :
968                t_prog_with_tags+=text
969           else:              
970                t_prog+=text
971    
972    print test_header
973    # print t_prog
974    print t_prog_with_tags
975    print test_tail          
976    1/0
977    #============================================================================================
978  def innerTEST(arg0,arg1):  def innerTEST(arg0,arg1):
979      if isinstance(arg0,float):      if isinstance(arg0,float):
980         out=numarray.array(arg0*arg1)         out=numarray.array(arg0*arg1)
# Line 696  def minimumTEST(arg0,arg1): Line 1146  def minimumTEST(arg0,arg1):
1146                else:                else:
1147                 out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]                 out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
1148       return out       return out
1149        
1150  def unrollLoops(a,b,o,arg,tap="",x="x"):  def unrollLoops(a,b,o,arg,tap="",x="x"):
1151      out=""      out=""
1152      if a.rank==1:      if a.rank==1:
# Line 797  def unrollLoopsOfGrad(a,b,o,arg,tap=""): Line 1248  def unrollLoopsOfGrad(a,b,o,arg,tap=""):
1248                 else:                 else:
1249                   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])
1250      return out      return out
1251    def unrollLoopsOfDiv(a,b,o,arg,tap=""):
1252        out=tap+arg+"="
1253        if o=="1":
1254           z=0.
1255           for i99 in range(a.shape[0]):
1256                z+=b[i99,i99]+a[i99,i99]
1257           out+="(%s)"%z    
1258        else:
1259           z=0.
1260           for i99 in range(a.shape[0]):
1261                z+=b[i99,i99]
1262                if i99>0: out+="+"
1263                out+="o*(%s)*x_ref[%s]**(o-1)"%(a[i99,i99],i99)
1264           out+="+(%s)"%z    
1265        return out
1266    
1267  def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):  def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
1268      if where=="Function":      if where=="Function":
1269         xfac_o=1.         xfac_o=1.
# Line 916  def unrollLoopsOfInteriorIntegral(a,b,wh Line 1383  def unrollLoopsOfInteriorIntegral(a,b,wh
1383                 out+="+(%s)*0.5**o\n"%zop                 out+="+(%s)*0.5**o\n"%zop
1384    
1385      return out      return out
1386    def unrollLoopsSimplified(b,arg,tap=""):
1387        out=""
1388        if isinstance(b,float) or b.rank==0:
1389                 out+=tap+"%s=(%s)*x[0]\n"%(arg,str(b))
1390    
1391        elif b.rank==1:
1392            for i0 in range(b.shape[0]):
1393                 out+=tap+"%s[%s]=(%s)*x[%s]\n"%(arg,i0,b[i0],i0)
1394        elif b.rank==2:
1395            for i0 in range(b.shape[0]):
1396             for i1 in range(b.shape[1]):
1397                 out+=tap+"%s[%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,b[i0,i1],i1)
1398        elif b.rank==3:
1399            for i0 in range(b.shape[0]):
1400             for i1 in range(b.shape[1]):
1401               for i2 in range(b.shape[2]):
1402                 out+=tap+"%s[%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,b[i0,i1,i2],i2)
1403        elif b.rank==4:
1404            for i0 in range(b.shape[0]):
1405             for i1 in range(b.shape[1]):
1406               for i2 in range(b.shape[2]):
1407                for i3 in range(b.shape[3]):
1408                 out+=tap+"%s[%s,%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,i3,b[i0,i1,i2,i3],i3)
1409        return out
1410    
1411    def unrollLoopsOfL2(b,where,arg,tap=""):
1412        out=""
1413        z=[]
1414        if isinstance(b,float) or b.rank==0:
1415           z.append(b**2)
1416        elif b.rank==1:
1417            for i0 in range(b.shape[0]):
1418                 z.append(b[i0]**2)
1419        elif b.rank==2:
1420            for i1 in range(b.shape[1]):
1421               s=0
1422               for i0 in range(b.shape[0]):
1423                  s+=b[i0,i1]**2
1424               z.append(s)
1425        elif b.rank==3:
1426            for i2 in range(b.shape[2]):
1427              s=0
1428              for i0 in range(b.shape[0]):
1429                 for i1 in range(b.shape[1]):
1430                    s+=b[i0,i1,i2]**2
1431              z.append(s)
1432    
1433        elif b.rank==4:
1434          for i3 in range(b.shape[3]):
1435             s=0
1436             for i0 in range(b.shape[0]):
1437               for i1 in range(b.shape[1]):
1438                  for i2 in range(b.shape[2]):
1439                     s+=b[i0,i1,i2,i3]**2
1440             z.append(s)        
1441        if where=="Function":
1442           xfac_o=1.
1443           xfac_op=0.
1444           z_fac_s=""
1445           zo_fac_s=""
1446           zo_fac=1./3.
1447        elif where=="FunctionOnBoundary":
1448           xfac_o=1.
1449           xfac_op=0.
1450           z_fac_s="*dim"
1451           zo_fac_s="*(2.*dim+1.)/3."
1452           zo_fac=1.
1453        elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1454           xfac_o=0.
1455           xfac_op=1.
1456           z_fac_s=""
1457           zo_fac_s=""    
1458           zo_fac=1./3.    
1459        zo=0.
1460        zop=0.
1461        for i99 in range(len(z)):
1462               if i99==0:
1463                   zo+=xfac_o*z[i99]
1464                   zop+=xfac_op*z[i99]
1465               else:
1466                   zo+=z[i99]
1467        out+=tap+"%s=sqrt((%s)%s"%(arg,zo*zo_fac,zo_fac_s)
1468        if zop==0.:
1469           out+=")\n"
1470        else:
1471           out+="+(%s))\n"%(zop*0.5**2)
1472        return out
1473    #=======================================================================================================
1474    # transpose
1475    #=======================================================================================================
1476    def transposeTest(r,offset):
1477        if isinstance(r,float): return r
1478        s=r.shape
1479        s1=1
1480        for i in s[:offset]: s1*=i
1481        s2=1
1482        for i in s[offset:]: s2*=i
1483        out=numarray.reshape(r,(s1,s2))
1484        out.transpose()
1485        return numarray.resize(out,s[offset:]+s[:offset])
1486    
1487    name,tt="transpose",transposeTest
1488    for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1489      for sh0 in [ (), (3,), (4,5), (6,2,2),(3,2,3,4)]:
1490        for offset in range(len(sh0)+1):
1491                  text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1492                  tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1493                  text+="   def %s(self):\n"%tname
1494                  sh_t=sh0[offset:]+sh0[:offset]
1495    
1496    #              sh_t=list(sh0)
1497    #              sh_t[offset+1]=sh_t[offset]
1498    #              sh_t=tuple(sh_t)
1499    #              sh_r=[]
1500    #              for i in range(offset): sh_r.append(sh0[i])
1501    #              for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])              
1502    #              sh_r=tuple(sh_r)
1503    
1504                  a_0=makeArray(sh0,[-1.,1])
1505                  if case0 in ["taggedData", "expandedData"]:
1506                      a1_0=makeArray(sh0,[-1.,1])
1507                  else:
1508                      a1_0=a_0
1509                  r=tt(a_0,offset)
1510                  r1=tt(a1_0,offset)
1511                  text+=mkText(case0,"arg",a_0,a1_0)
1512                  text+="      res=%s(arg,%s)\n"%(name,offset)
1513                  if case0=="Symbol":
1514                     text+=mkText("array","s",a_0,a1_0)
1515                     text+="      sub=res.substitute({arg:s})\n"
1516                     res="sub"
1517                     text+=mkText("array","ref",r,r1)
1518                  else:
1519                     res="res"
1520                     text+=mkText(case0,"ref",r,r1)
1521                  text+=mkTypeAndShapeTest(case0,sh_t,"res")
1522                  text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1523                  
1524                  if case0 == "taggedData" :
1525                      t_prog_with_tags+=text
1526                  else:              
1527                      t_prog+=text
1528    
1529    print test_header
1530    # print t_prog
1531    print t_prog_with_tags
1532    print test_tail          
1533    1/0
1534    #=======================================================================================================
1535    # L2
1536    #=======================================================================================================
1537    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1538      for data in ["Data","Symbol"]:
1539        for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1540             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1541             tname="test_L2_on%s_from%s_rank%s"%(where,data,len(sh))
1542             text+="   def %s(self):\n"%tname
1543             text+="      \"\"\"\n"
1544             text+="      tests L2-norm of %s on the %s\n\n"%(data,where)
1545             text+="      assumptions: self.domain supports integration on %s\n"%where
1546             text+="      \"\"\"\n"
1547             text+="      dim=self.domain.getDim()\n"
1548             text+="      w=%s(self.domain)\n"%where
1549             text+="      x=w.getX()\n"
1550             o="1"
1551             if len(sh)>0:
1552                sh_2=sh[:len(sh)-1]+(2,)
1553                sh_3=sh[:len(sh)-1]+(3,)            
1554                b_2=makeArray(sh[:len(sh)-1]+(2,),[-1.,1])
1555                b_3=makeArray(sh[:len(sh)-1]+(3,),[-1.,1])
1556             else:
1557                sh_2=()
1558                sh_3=()
1559                b_2=makeArray(sh,[-1.,1])
1560                b_3=makeArray(sh,[-1.,1])
1561    
1562             if data=="Symbol":
1563                val="s"
1564                res="sub"
1565             else:
1566                val="arg"
1567                res="res"
1568             text+="      if dim==2:\n"
1569             if data=="Symbol":
1570                   text+="        arg=Symbol(shape=%s,dim=dim)\n"%str(sh_2)
1571    
1572             text+="        %s=Data(0,%s,w)\n"%(val,sh_2)
1573             text+=unrollLoopsSimplified(b_2,val,tap="        ")
1574             text+=unrollLoopsOfL2(b_2,where,"ref",tap="        ")
1575             text+="\n      else:\n"
1576             if data=="Symbol":
1577                   text+="        arg=Symbol(shape=%s,dim=dim)\n"%str(sh_3)
1578             text+="        %s=Data(0,%s,w)\n"%(val,sh_3)        
1579             text+=unrollLoopsSimplified(b_3,val,tap="        ")
1580             text+=unrollLoopsOfL2(b_3,where,"ref",tap="        ")
1581             text+="\n      res=L2(arg)\n"
1582             if data=="Symbol":
1583                text+="      sub=res.substitute({arg:s})\n"
1584                text+="      self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1585                text+="      self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1586             else:
1587                text+="      self.failUnless(isinstance(res,float),\"wrong type of result.\")\n"
1588             text+="      self.failUnlessAlmostEqual(%s,ref,int(-log10(self.RES_TOL)),\"wrong result\")\n"%res
1589             t_prog+=text
1590    print t_prog
1591    1/0
1592    
1593    #=======================================================================================================
1594    # div
1595    #=======================================================================================================
1596    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1597      for data in ["Data","Symbol"]:
1598         for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1599             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1600             tname="test_div_on%s_from%s_%s"%(where,data,case)
1601             text+="   def %s(self):\n"%tname
1602             text+="      \"\"\"\n"
1603             text+="      tests divergence of %s on the %s\n\n"%(data,where)
1604             text+="      assumptions: %s(self.domain) exists\n"%case
1605             text+="                   self.domain supports div on %s\n"%where
1606             text+="      \"\"\"\n"
1607             if case=="ReducedSolution":
1608                text+="      o=1\n"
1609                o="1"
1610             else:
1611                text+="      o=self.order\n"
1612                o="o"
1613             text+="      dim=self.domain.getDim()\n"
1614             text+="      w_ref=%s(self.domain)\n"%where
1615             text+="      x_ref=w_ref.getX()\n"
1616             text+="      w=%s(self.domain)\n"%case
1617             text+="      x=w.getX()\n"
1618             a_2=makeArray((2,2),[-1.,1])
1619             b_2=makeArray((2,2),[-1.,1])
1620             a_3=makeArray((3,3),[-1.,1])
1621             b_3=makeArray((3,3),[-1.,1])
1622             if data=="Symbol":
1623                text+="      arg=Symbol(shape=(dim,),dim=dim)\n"
1624                val="s"
1625                res="sub"
1626             else:
1627                val="arg"
1628                res="res"
1629             text+="      %s=Vector(0,w)\n"%val
1630             text+="      if dim==2:\n"
1631             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1632             text+=unrollLoopsOfDiv(a_2,b_2,o,"ref",tap="        ")
1633             text+="\n      else:\n"
1634            
1635             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1636             text+=unrollLoopsOfDiv(a_3,b_3,o,"ref",tap="        ")
1637             text+="\n      res=div(arg,where=w_ref)\n"
1638             if data=="Symbol":
1639                text+="      sub=res.substitute({arg:s})\n"
1640             text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1641             text+="      self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1642             text+="      self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong function space of result.\")\n"%res
1643             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1644    
1645    
1646             t_prog+=text
1647    print t_prog
1648    1/0
1649    
1650  #=======================================================================================================  #=======================================================================================================
1651  # interpolation  # interpolation
Legend:
Removed from v.439  
changed lines
  Added in v.587
  ViewVC Help
Powered by ViewVC 1.1.26