/[escript]/branches/clazy/escriptcore/src/DataLazy.cpp
ViewVC logotype

Diff of /branches/clazy/escriptcore/src/DataLazy.cpp

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

branches/schroedinger/escript/src/DataLazy.cpp revision 1899 by jfenwick, Mon Oct 20 05:13:24 2008 UTC branches/schroedinger_upto1946/escript/src/DataLazy.cpp revision 1993 by phornby, Fri Nov 7 04:52:15 2008 UTC
# Line 26  Line 26 
26  #include "DataTypes.h"  #include "DataTypes.h"
27  #include "Data.h"  #include "Data.h"
28  #include "UnaryFuncs.h"     // for escript::fsign  #include "UnaryFuncs.h"     // for escript::fsign
29    #include "Utils.h"
30    
31  /*  /*
32  How does DataLazy work?  How does DataLazy work?
# Line 78  using namespace boost; Line 79  using namespace boost;
79  namespace escript  namespace escript
80  {  {
81    
 const std::string&  
 opToString(ES_optype op);  
   
82  namespace  namespace
83  {  {
84    
# Line 95  enum ES_opgroup Line 93  enum ES_opgroup
93    
94    
95    
96  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","sin","cos","tan",  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
97                "sin","cos","tan",
98              "asin","acos","atan","sinh","cosh","tanh","erf",              "asin","acos","atan","sinh","cosh","tanh","erf",
99              "asinh","acosh","atanh",              "asinh","acosh","atanh",
100              "log10","log","sign","abs","neg","pos","exp","sqrt",              "log10","log","sign","abs","neg","pos","exp","sqrt",
101              "1/","where>0","where<0","where>=0","where<=0"};              "1/","where>0","where<0","where>=0","where<=0"};
102  int ES_opcount=32;  int ES_opcount=33;
103  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY,G_UNARY,G_UNARY,G_UNARY, //9  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
104              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 16              G_UNARY,G_UNARY,G_UNARY, //10
105              G_UNARY,G_UNARY,G_UNARY,                    // 19              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17
106              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 27              G_UNARY,G_UNARY,G_UNARY,                    // 20
107                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 28
108              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};
109  inline  inline
110  ES_opgroup  ES_opgroup
# Line 122  resultFS(DataAbstract_ptr left, DataAbst Line 122  resultFS(DataAbstract_ptr left, DataAbst
122      // that way, if interpolate is required in any other op we can just throw a      // that way, if interpolate is required in any other op we can just throw a
123      // programming error exception.      // programming error exception.
124    
125      FunctionSpace l=left->getFunctionSpace();
126      if (left->getFunctionSpace()!=right->getFunctionSpace())    FunctionSpace r=right->getFunctionSpace();
127      {    if (l!=r)
128          throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");    {
129      }      if (r.probeInterpolation(l))
130      return left->getFunctionSpace();      {
131        return l;
132        }
133        if (l.probeInterpolation(r))
134        {
135        return r;
136        }
137        throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
138      }
139      return l;
140  }  }
141    
142  // return the shape of the result of "left op right"  // return the shape of the result of "left op right"
# Line 136  resultShape(DataAbstract_ptr left, DataA Line 145  resultShape(DataAbstract_ptr left, DataA
145  {  {
146      if (left->getShape()!=right->getShape())      if (left->getShape()!=right->getShape())
147      {      {
148          throw DataException("Shapes not the same - shapes must match for lazy data.");        if (getOpgroup(op)!=G_BINARY)
149          {
150            throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
151          }
152          if (left->getRank()==0)   // we need to allow scalar * anything
153          {
154            return right->getShape();
155          }
156          if (right->getRank()==0)
157          {
158            return left->getShape();
159          }
160          throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
161      }      }
162      return left->getShape();      return left->getShape();
163  }  }
# Line 212  cout << "(1)Lazy created with " << m_sam Line 233  cout << "(1)Lazy created with " << m_sam
233    
234    
235    
236    
237  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
238      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
239      m_op(op)      m_op(op)
# Line 237  DataLazy::DataLazy(DataAbstract_ptr left Line 259  DataLazy::DataLazy(DataAbstract_ptr left
259  }  }
260    
261    
262    // In this constructor we need to consider interpolation
263  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
264      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
265      m_op(op)      m_op(op)
# Line 245  DataLazy::DataLazy(DataAbstract_ptr left Line 268  DataLazy::DataLazy(DataAbstract_ptr left
268     {     {
269      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
270     }     }
271    
272       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
273       {
274        FunctionSpace fs=getFunctionSpace();
275        Data ltemp(left);
276        Data tmp(ltemp,fs);
277        left=tmp.borrowDataPtr();
278       }
279       if (getFunctionSpace()!=right->getFunctionSpace())   // left needs to be interpolated
280       {
281        Data tmp(Data(right),getFunctionSpace());
282        right=tmp.borrowDataPtr();
283       }
284       left->operandCheck(*right);
285    
286     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
287     {     {
288      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
# Line 468  DataLazy::resolveUnary(ValueType& v, siz Line 506  DataLazy::resolveUnary(ValueType& v, siz
506    switch (m_op)    switch (m_op)
507    {    {
508      case SIN:        case SIN:  
509      tensor_unary_operation(m_samplesize, left, result, ::sin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
510      break;      break;
511      case COS:      case COS:
512      tensor_unary_operation(m_samplesize, left, result, ::cos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
513      break;      break;
514      case TAN:      case TAN:
515      tensor_unary_operation(m_samplesize, left, result, ::tan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
516      break;      break;
517      case ASIN:      case ASIN:
518      tensor_unary_operation(m_samplesize, left, result, ::asin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
519      break;      break;
520      case ACOS:      case ACOS:
521      tensor_unary_operation(m_samplesize, left, result, ::acos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
522      break;      break;
523      case ATAN:      case ATAN:
524      tensor_unary_operation(m_samplesize, left, result, ::atan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
525      break;      break;
526      case SINH:      case SINH:
527      tensor_unary_operation(m_samplesize, left, result, ::sinh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
528      break;      break;
529      case COSH:      case COSH:
530      tensor_unary_operation(m_samplesize, left, result, ::cosh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
531      break;      break;
532      case TANH:      case TANH:
533      tensor_unary_operation(m_samplesize, left, result, ::tanh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
534      break;      break;
535      case ERF:      case ERF:
536  #ifdef _WIN32  #ifdef _WIN32
# Line 523  DataLazy::resolveUnary(ValueType& v, siz Line 561  DataLazy::resolveUnary(ValueType& v, siz
561  #endif    #endif  
562      break;      break;
563      case LOG10:      case LOG10:
564      tensor_unary_operation(m_samplesize, left, result, ::log10);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
565      break;      break;
566      case LOG:      case LOG:
567      tensor_unary_operation(m_samplesize, left, result, ::log);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
568      break;      break;
569      case SIGN:      case SIGN:
570      tensor_unary_operation(m_samplesize, left, result, escript::fsign);      tensor_unary_operation(m_samplesize, left, result, escript::fsign);
571      break;      break;
572      case ABS:      case ABS:
573      tensor_unary_operation(m_samplesize, left, result, ::fabs);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
574      break;      break;
575      case NEG:      case NEG:
576      tensor_unary_operation(m_samplesize, left, result, negate<double>());      tensor_unary_operation(m_samplesize, left, result, negate<double>());
# Line 543  DataLazy::resolveUnary(ValueType& v, siz Line 581  DataLazy::resolveUnary(ValueType& v, siz
581      throw DataException("Programmer error - POS not supported for lazy data.");      throw DataException("Programmer error - POS not supported for lazy data.");
582      break;      break;
583      case EXP:      case EXP:
584      tensor_unary_operation(m_samplesize, left, result, ::exp);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
585      break;      break;
586      case SQRT:      case SQRT:
587      tensor_unary_operation(m_samplesize, left, result, ::sqrt);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
588      break;      break;
589      case RECIP:      case RECIP:
590      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
# Line 574  DataLazy::resolveUnary(ValueType& v, siz Line 612  DataLazy::resolveUnary(ValueType& v, siz
612    
613    
614    
615  #define PROC_OP(X) \  #define PROC_OP(TYPE,X)                               \
616      for (int i=0;i<steps;++i,resultp+=getNoValues()) \      for (int i=0;i<steps;++i,resultp+=resultStep) \
617      { \      { \
618         tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \         tensor_binary_operation##TYPE(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
619         lroffset+=leftStep; \         lroffset+=leftStep; \
620         rroffset+=rightStep; \         rroffset+=rightStep; \
621      }      }
# Line 600  DataLazy::resolveUnary(ValueType& v, siz Line 638  DataLazy::resolveUnary(ValueType& v, siz
638  // the whole sample as one big datapoint.  // the whole sample as one big datapoint.
639  // If one of the children is not expanded, then we need to treat each point in the sample  // If one of the children is not expanded, then we need to treat each point in the sample
640  // individually.  // individually.
641    // There is an additional complication when scalar operations are considered.
642    // For example, 2+Vector.
643    // In this case each double within the point is treated individually
644  DataTypes::ValueType*  DataTypes::ValueType*
645  DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
646  {  {
# Line 610  cout << "Resolve binary: " << toString() Line 651  cout << "Resolve binary: " << toString()
651    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
652    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
653    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
654    int steps=(bigloops?1:getNumDPPSample());        int steps=(bigloops?1:getNumDPPSample());
655    size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint    size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint
656    int leftStep=((leftExp && !rightExp)? getNoValues() : 0);    if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops
657    int rightStep=((rightExp && !leftExp)? getNoValues() : 0);    {
658        EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");
659        steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());
660        chunksize=1;    // for scalar
661      }    
662      int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
663      int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
664      int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0
665      // Get the values of sub-expressions      // Get the values of sub-expressions
666    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
667    const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note    const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
# Line 622  cout << "Resolve binary: " << toString() Line 670  cout << "Resolve binary: " << toString()
670    switch(m_op)    switch(m_op)
671    {    {
672      case ADD:      case ADD:
673      PROC_OP(plus<double>());          PROC_OP(/**/,plus<double>());
674      break;      break;
675      case SUB:      case SUB:
676      PROC_OP(minus<double>());      PROC_OP(/**/,minus<double>());
677      break;      break;
678      case MUL:      case MUL:
679      PROC_OP(multiplies<double>());      PROC_OP(/**/,multiplies<double>());
680      break;      break;
681      case DIV:      case DIV:
682      PROC_OP(divides<double>());      PROC_OP(/**/,divides<double>());
683        break;
684        case POW:
685           PROC_OP(<double (double,double)>,::pow);
686      break;      break;
687      default:      default:
688      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
# Line 799  DataLazy::intoString(ostringstream& oss) Line 850  DataLazy::intoString(ostringstream& oss)
850    }    }
851  }  }
852    
 // Note that in this case, deepCopy does not make copies of the leaves.  
 // Hopefully copy on write (or whatever we end up using) will take care of this.  
853  DataAbstract*  DataAbstract*
854  DataLazy::deepCopy()  DataLazy::deepCopy()
855  {  {
856    if (m_op==IDENTITY)    switch (getOpgroup(m_op))
857    {    {
858      return new DataLazy(m_left);    // we don't need to copy the child here    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
859      case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
860      case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
861      default:
862        throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
863    }    }
   return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);  
864  }  }
865    
866    
# Line 825  DataLazy::getSlice(const DataTypes::Regi Line 877  DataLazy::getSlice(const DataTypes::Regi
877    throw DataException("getSlice - not implemented for Lazy objects.");    throw DataException("getSlice - not implemented for Lazy objects.");
878  }  }
879    
880    
881    // To do this we need to rely on our child nodes
882    DataTypes::ValueType::size_type
883    DataLazy::getPointOffset(int sampleNo,
884                     int dataPointNo)
885    {
886      if (m_op==IDENTITY)
887      {
888        return m_id->getPointOffset(sampleNo,dataPointNo);
889      }
890      if (m_readytype!='E')
891      {
892        collapse();
893        return m_id->getPointOffset(sampleNo,dataPointNo);
894      }
895      // at this point we do not have an identity node and the expression will be Expanded
896      // so we only need to know which child to ask
897      if (m_left->m_readytype=='E')
898      {
899        return m_left->getPointOffset(sampleNo,dataPointNo);
900      }
901      else
902      {
903        return m_right->getPointOffset(sampleNo,dataPointNo);
904      }
905    }
906    
907    // To do this we need to rely on our child nodes
908  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
909  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
910                   int dataPointNo) const                   int dataPointNo) const
911  {  {
912    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");    if (m_op==IDENTITY)
913      {
914        return m_id->getPointOffset(sampleNo,dataPointNo);
915      }
916      if (m_readytype=='E')
917      {
918        // at this point we do not have an identity node and the expression will be Expanded
919        // so we only need to know which child to ask
920        if (m_left->m_readytype=='E')
921        {
922        return m_left->getPointOffset(sampleNo,dataPointNo);
923        }
924        else
925        {
926        return m_right->getPointOffset(sampleNo,dataPointNo);
927        }
928      }
929      if (m_readytype=='C')
930      {
931        return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
932      }
933      throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
934    }
935    
936    // It would seem that DataTagged will need to be treated differently since even after setting all tags
937    // to zero, all the tags from all the DataTags would be in the result.
938    // However since they all have the same value (0) whether they are there or not should not matter.
939    // So I have decided that for all types this method will create a constant 0.
940    // It can be promoted up as required.
941    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
942    // but we can deal with that if it arrises.
943    void
944    DataLazy::setToZero()
945    {
946      DataTypes::ValueType v(getNoValues(),0);
947      m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
948      m_op=IDENTITY;
949      m_right.reset();  
950      m_left.reset();
951      m_readytype='C';
952      m_buffsRequired=1;
953  }  }
954    
955  }   // end namespace  }   // end namespace

Legend:
Removed from v.1899  
changed lines
  Added in v.1993

  ViewVC Help
Powered by ViewVC 1.1.26