/[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 433 by gross, Tue Jan 17 23:54:38 2006 UTC revision 437 by gross, Fri Jan 20 00:16:58 2006 UTC
# Line 173  def makeArray(shape,rng): Line 173  def makeArray(shape,rng):
173               for i2 in range(shape[2]):               for i2 in range(shape[2]):
174                  for i3 in range(shape[3]):                  for i3 in range(shape[3]):
175                     out[i0,i1,i2,i3]=l*random.random()+rng[0]                     out[i0,i1,i2,i3]=l*random.random()+rng[0]
176       elif len(shape)==5:
177           for i0 in range(shape[0]):
178              for i1 in range(shape[1]):
179                 for i2 in range(shape[2]):
180                    for i3 in range(shape[3]):
181                      for i4 in range(shape[4]):
182                       out[i0,i1,i2,i3,i4]=l*random.random()+rng[0]
183     else:     else:
184         raise SystemError,"rank is restricted to 4"         raise SystemError,"rank is restricted to 5"
185     return out             return out        
186    
187    
# Line 689  def minimumTEST(arg0,arg1): Line 696  def minimumTEST(arg0,arg1):
696                else:                else:
697                 out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]                 out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
698       return out       return out
699    def unrollLoops(a,b,o,arg,tap=""):
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)*x[%s]"%(a[i99]+b[i99],i99)
707                   else:
708                      z+="(%s)*x[%s]**o+(%s)*x[%s]"%(a[i99],i99,b[i99],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)*x[%s]**o+(%s)*x[%s]"%(a[i0,i99],i99,b[i0,i99],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)*x[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],i99)
731                   else:
732                      z+="(%s)*x[%s]**o+(%s)*x[%s]"%(a[i0,i1,i99],i99,b[i0,i1,i99],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)*x[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],i99)
744                   else:
745                      z+="(%s)*x[%s]**o+(%s)*x[%s]"%(a[i0,i1,i2,i99],i99,b[i0,i1,i2,i99],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)*x[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],i99)
758                   else:
759                      z+="(%s)*x[%s]**o+(%s)*x[%s]"%(a[i0,i1,i2,i3,i99],i99,b[i0,i1,i2,i3,i99],i99)
760    
761                 out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
762        return out
763    
764    def unrollLoopsOfGrad(a,b,o,arg,tap=""):
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    
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    
922    #=======================================================================================================
923    # integrate
924    #=======================================================================================================
925    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
926      for data in ["Data","Symbol"]:
927        for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
928          for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
929            if (not case in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]) or where==case:  
930             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
931             tname="test_integrate_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
932             text+="   def %s(self):\n"%tname
933             text+="      \"\"\"\n"
934             text+="      tests integral of rank %s %s on the %s\n\n"%(len(sh),data,where)
935             text+="      assumptions: %s(self.domain) exists\n"%case
936             text+="                   self.domain supports integral on %s\n"%where
937    
938             text+="      \"\"\"\n"
939             if case=="ReducedSolution":
940                text+="      o=1\n"
941                o="1"
942             else:
943                text+="      o=self.order\n"
944                o="o"
945             text+="      dim=self.domain.getDim()\n"
946             text+="      w_ref=%s(self.domain)\n"%where
947             text+="      w=%s(self.domain)\n"%case
948             text+="      x=w.getX()\n"
949             a_2=makeArray(sh+(2,),[-1.,1])
950             b_2=makeArray(sh+(2,),[-1.,1])
951             a_3=makeArray(sh+(3,),[-1.,1])
952             b_3=makeArray(sh+(3,),[-1.,1])
953             if data=="Symbol":
954                text+="      arg=Symbol(shape=%s)\n"%str(sh)
955                val="s"
956                res="sub"
957             else:
958                val="arg"
959                res="res"
960                
961             text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
962             if not len(sh)==0:
963                text+="      ref=numarray.zeros(%s,numarray.Float)\n"%str(sh)
964             text+="      if dim==2:\n"
965             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
966             text+=unrollLoopsOfInteriorIntegral(a_2,b_2,where,"ref",tap="        ")
967             text+="      else:\n"
968            
969             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
970             text+=unrollLoopsOfInteriorIntegral(a_3,b_3,where,"ref",tap="        ")
971             if case in ["ContinuousFunction","Solution","ReducedSolution"]:
972                 text+="      res=integrate(arg,where=w_ref)\n"
973             else:
974                 text+="      res=integrate(arg)\n"
975    
976             if data=="Symbol":
977                text+="      sub=res.substitute({arg:s})\n"
978             if len(sh)==0 and data=="Data":
979                text+="      self.failUnless(isinstance(%s,float),\"wrong type of result.\")\n"%res
980             else:
981                if data=="Symbol":
982                   text+="      self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
983                   text+="      self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
984                else:
985                   text+="      self.failUnless(isinstance(res,numarray.NumArray),\"wrong type of result.\")\n"
986                   text+="      self.failUnlessEqual(res.shape,%s,\"wrong shape of result.\")\n"%str(sh)
987             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
988    
989    
990             t_prog+=text
991    print test_header
992    print t_prog
993    print test_tail          
994    1/0
995    #=======================================================================================================
996    # grad
997    #=======================================================================================================
998    for where in ["Function","FunctionOnBoundary"]:
999      for data in ["Data","Symbol"]:
1000         for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1001           for sh in [ (),(2,), (4,5), (6,2,2)]:
1002             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1003             tname="test_grad_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1004             text+="   def %s(self):\n"%tname
1005             text+="      \"\"\"\n"
1006             text+="      tests gradient for rank %s %s on the %s\n\n"%(len(sh),data,where)
1007             text+="      assumptions: %s(self.domain) exists\n"%where
1008             if where=="FunctionOnBoundary": text+="                   self.domain supports gradient on the boundary\n"
1009             text+="      \"\"\"\n"
1010             if case=="ReducedSolution":
1011                text+="      o=1\n"
1012                o="1"
1013             else:
1014                text+="      o=self.order\n"
1015                o="o"
1016             text+="      dim=self.domain.getDim()\n"
1017             text+="      w_ref=%s(self.domain)\n"%where
1018             text+="      x_ref=w_ref.getX()\n"
1019             text+="      w=%s(self.domain)\n"%case
1020             text+="      x=w.getX()\n"
1021             a_2=makeArray(sh+(2,),[-1.,1])
1022             b_2=makeArray(sh+(2,),[-1.,1])
1023             a_3=makeArray(sh+(3,),[-1.,1])
1024             b_3=makeArray(sh+(3,),[-1.,1])
1025             if data=="Symbol":
1026                text+="      arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1027                val="s"
1028                res="sub"
1029             else:
1030                val="arg"
1031                res="res"
1032             text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
1033             text+="      ref=Data(0,%s+(dim,),w_ref)\n"%str(sh)
1034             text+="      if dim==2:\n"
1035             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1036             text+=unrollLoopsOfGrad(a_2,b_2,o,"ref",tap="        ")
1037             text+="      else:\n"
1038            
1039             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1040             text+=unrollLoopsOfGrad(a_3,b_3,o,"ref",tap="        ")
1041             text+="      res=grad(arg,where=w_ref)\n"
1042             if data=="Symbol":
1043                text+="      sub=res.substitute({arg:s})\n"
1044             text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1045             text+="      self.failUnlessEqual(res.getShape(),%s+(dim,),\"wrong shape of result.\")\n"%str(sh)
1046             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1047    
1048    
1049             t_prog+=text
1050    print test_header
1051    print t_prog
1052    print test_tail          
1053    1/0
1054    
1055  #=======================================================================================================  #=======================================================================================================
1056  # inverse  # inverse
1057  #=======================================================================================================  #=======================================================================================================

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

  ViewVC Help
Powered by ViewVC 1.1.26