SAGA API  v9.5
mat_formula.cpp
Go to the documentation of this file.
1 
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Application Programming Interface //
9 // //
10 // Library: SAGA_API //
11 // //
12 //-------------------------------------------------------//
13 // //
14 // mat_formula.cpp //
15 // //
16 // Copyright (C) 2002 by Andre Ringeler //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. //
22 // //
23 // This library is free software; you can redistribute //
24 // it and/or modify it under the terms of the GNU Lesser //
25 // General Public License as published by the Free //
26 // Software Foundation, either version 2.1 of the //
27 // License, or (at your option) any later version. //
28 // //
29 // This library is distributed in the hope that it will //
30 // be useful, but WITHOUT ANY WARRANTY; without even the //
31 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
32 // PARTICULAR PURPOSE. See the GNU Lesser General Public //
33 // License for more details. //
34 // //
35 // You should have received a copy of the GNU Lesser //
36 // General Public License along with this program; if //
37 // not, see <http://www.gnu.org/licenses/>. //
38 // //
39 //-------------------------------------------------------//
40 // //
41 // e-mail: aringel@saga-gis.org //
42 // //
43 // contact: Andre Ringeler //
44 // Institute of Geography //
45 // University of Goettingen //
46 // Goldschmidtstr. 5 //
47 // 37077 Goettingen //
48 // Germany //
49 // //
50 //-------------------------------------------------------//
51 // //
52 // Based on //
53 // FORMULC.C 2.22 as of 2/19/98 //
54 // Copyright (c) 1993 - 98 by Harald Helfgott //
55 // //
56 // Modified for Grid Data by Andre Ringeler 2001 //
57 // Modified for Function-Fitting by Andre Ringeler 2002 //
58 // Converted to C++ by Andre Ringeler 2002 //
59 // //
60 // Modified to fit SAGA needs by Olaf Conrad 2002 //
61 // //
63 
64 //---------------------------------------------------------
65 #include <string.h>
66 #include <stdlib.h>
67 #include <stdio.h>
68 #include <stdarg.h>
69 #include <errno.h>
70 #include <ctype.h>
71 
72 #include "mat_tools.h"
73 #include "grid.h"
74 
75 
77 // //
78 // //
79 // //
81 
82 //---------------------------------------------------------
83 #define MAX_PARMS 3
84 #define MAX_CTABLE 255
85 
86 #define STD_FNC_NUM 19
87 
88 //---------------------------------------------------------
89 #define GET_VALUE_BUFSIZE 500
90 
91 //---------------------------------------------------------
92 #define EPSILON 1e-9
93 
94 
96 // //
97 // //
98 // //
100 
101 //---------------------------------------------------------
102 static double f_atan2(double x, double val)
103 {
104  return( atan2(x, val) );
105 }
106 
107 //---------------------------------------------------------
108 static double f_pow(double x, double val)
109 {
110  return( pow(x, val) );
111 }
112 
113 //---------------------------------------------------------
114 static double f_gt(double x, double val)
115 {
116  return( x > val ? 1.0 : 0.0 );
117 }
118 
119 //---------------------------------------------------------
120 static double f_lt(double x, double val)
121 {
122  return( x < val ? 1.0 : 0.0 );
123 }
124 
125 //---------------------------------------------------------
126 static double f_eq(double x, double val)
127 {
128  return( fabs(x - val) < EPSILON ? 1.0 : 0.0 );
129 }
130 
131 //---------------------------------------------------------
132 static double f_min(double a, double b)
133 {
134  return( a < b ? a : b );
135 }
136 
137 //---------------------------------------------------------
138 static double f_max(double a, double b)
139 {
140  return( a > b ? a : b );
141 }
142 
143 //---------------------------------------------------------
144 static double f_pi(void)
145 {
146  return( M_PI );
147 }
148 
149 //---------------------------------------------------------
150 static double f_int(double x)
151 {
152  return( (int)(x) );
153 }
154 
155 //---------------------------------------------------------
156 static double f_sqr(double x)
157 {
158  return( x*x );
159 }
160 
161 //---------------------------------------------------------
162 static double f_fmod(double x, double val)
163 {
164  return( fmod(x, val) );
165 }
166 
167 //---------------------------------------------------------
168 static double f_rand_u(double min, double max)
169 {
170  return( CSG_Random::Get_Uniform(min, max) );
171 }
172 
173 //---------------------------------------------------------
174 static double f_rand_g(double mean, double stdv)
175 {
176  return( CSG_Random::Get_Gaussian(mean, stdv) );
177 }
178 
179 //---------------------------------------------------------
180 static double f_and(double x, double y)
181 {
182  return( x != 0.0 && y != 0.0 ? 1.0 : 0.0 );
183 }
184 
185 //---------------------------------------------------------
186 static double f_or(double x, double y)
187 {
188  return( x != 0.0 || y != 0.0 ? 1.0 : 0.0 );
189 }
190 
191 //---------------------------------------------------------
192 static double f_ifelse(double condition, double x, double y)
193 {
194  return( condition ? x : y );
195 // return( fabs(condition) >= EPSILON ? x : y );
196 }
197 
198 
200 // //
201 // //
202 // //
204 
205 //---------------------------------------------------------
206 static CSG_Formula::TSG_Function gSG_Functions[MAX_CTABLE] =
207 {
208  { "exp" , (TSG_Formula_Function_1)exp , 1, false }, // 1
209  { "ln" , (TSG_Formula_Function_1)log , 1, false }, // 2
210  { "sin" , (TSG_Formula_Function_1)sin , 1, false }, // 3
211  { "cos" , (TSG_Formula_Function_1)cos , 1, false }, // 4
212  { "tan" , (TSG_Formula_Function_1)tan , 1, false }, // 5
213  { "asin" , (TSG_Formula_Function_1)asin , 1, false }, // 6
214  { "acos" , (TSG_Formula_Function_1)acos , 1, false }, // 7
215  { "atan" , (TSG_Formula_Function_1)atan , 1, false }, // 8
216  { "atan2" , (TSG_Formula_Function_1)f_atan2 , 2, false }, // 9
217  { "abs" , (TSG_Formula_Function_1)fabs , 1, false }, // 10
218  { "sqrt" , (TSG_Formula_Function_1)sqrt , 1, false }, // 11
219  { "gt" , (TSG_Formula_Function_1)f_gt , 2, false }, // 12
220  { "lt" , (TSG_Formula_Function_1)f_lt , 2, false }, // 13
221  { "eq" , (TSG_Formula_Function_1)f_eq , 2, false }, // 14
222  { "pi" , (TSG_Formula_Function_1)f_pi , 0, false }, // 15
223  { "int" , (TSG_Formula_Function_1)f_int , 1, false }, // 16
224  { "mod" , (TSG_Formula_Function_1)f_fmod , 2, false }, // 17
225  { "ifelse", (TSG_Formula_Function_1)f_ifelse, 3, false }, // 18
226  { "log" , (TSG_Formula_Function_1)log10 , 1, false }, // 19
227  { "pow" , (TSG_Formula_Function_1)f_pow , 2, false }, // 20
228  { "sqr" , (TSG_Formula_Function_1)f_sqr , 1, false }, // 21
229  { "rand_u", (TSG_Formula_Function_1)f_rand_u, 2, true }, // 22
230  { "rand_g", (TSG_Formula_Function_1)f_rand_g, 2, true }, // 23
231  { "and" , (TSG_Formula_Function_1)f_and , 2, false }, // 24
232  { "or" , (TSG_Formula_Function_1)f_or , 2, false }, // 25
233  { "min" , (TSG_Formula_Function_1)f_min , 2, false }, // 26
234  { "max" , (TSG_Formula_Function_1)f_max , 2, false }, // 27
235  { NULL , (TSG_Formula_Function_1) NULL , 0, false }
236 };
237 
238 
240 // //
241 // //
242 // //
244 
245 //---------------------------------------------------------
246 CSG_Formula::CSG_Formula(void)
247 {
248  m_Formula.code = NULL;
249  m_Formula.ctable = NULL;
250 
251  m_bError = false;
252 
253  m_ctable = NULL;
254  m_error = NULL;
255 
256  //-----------------------------------------------------
257  m_Functions = (TSG_Function *)SG_Calloc(MAX_CTABLE, sizeof(TSG_Function));
258 
259  for(int i=0; i<MAX_CTABLE; i++)
260  {
261  m_Functions[i].Name = gSG_Functions[i].Name;
262  m_Functions[i].Function = gSG_Functions[i].Function;
263  m_Functions[i].nParameters = gSG_Functions[i].nParameters;
264  m_Functions[i].bVarying = gSG_Functions[i].bVarying;
265  }
266 }
267 
268 //---------------------------------------------------------
269 CSG_Formula::~CSG_Formula(void)
270 {
271  Destroy();
272 
273  SG_Free(m_Functions);
274 }
275 
276 //---------------------------------------------------------
277 bool CSG_Formula::Destroy(void)
278 {
279  SG_FREE_SAFE(m_Formula.code);
280  SG_FREE_SAFE(m_Formula.ctable);
281 
282  m_bError = false;
283 
284  return( true );
285 }
286 
287 
289 // //
290 // //
291 // //
293 
294 //---------------------------------------------------------
295 CSG_String CSG_Formula::Get_Help_Operators(bool bHTML, const CSG_String Additional[][2])
296 {
297  const int nOperators = 35;
298 
299  CSG_String Operators[nOperators][2] =
300  {
301  { "+" , _TL("Addition") },
302  { "-" , _TL("Subtraction") },
303  { "*" , _TL("Multiplication") },
304  { "/" , _TL("Division") },
305  { "abs(x)" , _TL("Absolute Value") },
306  { "mod(x, y)" , _TL("Returns the floating point remainder of x/y") },
307  { "int(x)" , _TL("Returns the integer part of floating point value x") },
308  { "sqr(x)" , _TL("Square") },
309  { "sqrt(x)" , _TL("Square Root") },
310  { "exp(x)" , _TL("Exponential") },
311  { "pow(x, y)" , _TL("Returns x raised to the power of y") },
312  { "x ^ y" , _TL("Returns x raised to the power of y") },
313  { "ln(x)" , _TL("Natural Logarithm") },
314  { "log(x)" , _TL("Base 10 Logarithm") },
315  { "pi()" , _TL("Returns the value of Pi") },
316  { "sin(x)" , _TL("Sine, expects radians") },
317  { "cos(x)" , _TL("Cosine, expects radians") },
318  { "tan(x)" , _TL("Tangent, expects radians") },
319  { "asin(x)" , _TL("Arcsine, returns radians") },
320  { "acos(x)" , _TL("Arccosine, returns radians") },
321  { "atan(x)" , _TL("Arctangent, returns radians") },
322  { "atan2(x, y)" , _TL("Arctangent of x/y, returns radians") },
323  { "min(x, y)" , _TL("Returns the minimum of values x and y") },
324  { "max(x, y)" , _TL("Returns the maximum of values x and y") },
325  { "gt(x, y)" , _TL("Returns true (1), if x is greater than y, else false (0)") },
326  { "x > y" , _TL("Returns true (1), if x is greater than y, else false (0)") },
327  { "lt(x, y)" , _TL("Returns true (1), if x is less than y, else false (0)") },
328  { "x < y" , _TL("Returns true (1), if x is less than y, else false (0)") },
329  { "eq(x, y)" , _TL("Returns true (1), if x equals y, else false (0)") },
330  { "x = y" , _TL("Returns true (1), if x equals y, else false (0)") },
331  { "and(x, y)" , _TL("Returns true (1), if both x and y are true (i.e. not 0)") },
332  { "or(x, y)" , _TL("Returns true (1), if at least one of both x and y is true (i.e. not 0)") },
333  { "ifelse(c, x, y)", _TL("Returns x, if condition c is true (i.e. not 0), else y") },
334  { "rand_u(x, y)" , _TL("Random number, uniform distribution with minimum x and maximum y") },
335  { "rand_g(x, y)" , _TL("Random number, Gaussian distribution with mean x and standard deviation y") }
336  };
337 
338  //-----------------------------------------------------
339  int i;
340  CSG_String s;
341 
342  if( bHTML )
343  {
344  s += "<table border=\"0\">";
345 
346  for(i=0; i<nOperators; i++)
347  {
348  CSG_String op = Operators[i][0]; op.Replace("<", "&lt;");
349 
350  s += "<tr><td><b>" + op + "</b></td><td>" + Operators[i][1] + "</td></tr>";
351  }
352 
353  if( Additional )
354  {
355  for(i=0; !Additional[i][0].is_Empty(); i++)
356  {
357  CSG_String op = Additional[i][0]; op.Replace("<", "&lt;");
358 
359  s += "<tr><td><b>" + op + "</b></td><td>" + Additional[i][1] + "</td></tr>";
360  }
361  }
362 
363  s += "</table>";
364  }
365  else
366  {
367  for(i=0; i<nOperators; i++)
368  {
369  s += Operators[i][0] + " - " + Operators[i][1] + "\n";
370  }
371 
372  if( Additional )
373  {
374  for(i=0; !Additional[i][0].is_Empty(); i++)
375  {
376  s += Additional[i][0] + " - " + Additional[i][1] + "\n";
377  }
378  }
379  }
380 
381  return( s );
382 }
383 
384 //---------------------------------------------------------
385 bool CSG_Formula::Get_Error(CSG_String &Message)
386 {
387  if( m_bError )
388  {
389  Message = CSG_String::Format("%s %s %d\n", _TL("Error in formula"), _TL("at position"), m_Error_Position);
390 
391  if( m_Error_Position < 0 || m_Error_Position >= (int)m_sFormula.Length() )
392  {
393  Message += m_sFormula;
394  }
395  else
396  {
397  Message += m_sFormula.Left (m_Error_Position) + " ["
398  + m_sFormula [m_Error_Position] + "] "
399  + m_sFormula.Right(m_sFormula.Length() - (m_Error_Position + 1));
400  }
401 
402  Message += "\n";
403  Message += m_sError;
404  Message += "\n";
405 
406  return( true );
407  }
408 
409  return( false );
410 }
411 
412 //---------------------------------------------------------
413 void CSG_Formula::_Set_Error(const CSG_String &Error)
414 {
415  if( Error.is_Empty() )
416  {
417  m_bError = false;
418  m_sError .Clear();
419  }
420  else
421  {
422  m_bError = true;
423  m_sError = Error;
424  }
425 }
426 
427 
429 // //
431 
432 //---------------------------------------------------------
433 void CSG_Formula::Set_Variable(char Var, double Value)
434 {
435  m_Parameters[Var - 'a'] = Value;
436 }
437 
438 //---------------------------------------------------------
439 bool CSG_Formula::Set_Formula(const CSG_String &Formula)
440 {
441  if( Formula.Length() > 0 )
442  {
443  Destroy();
444 
445  m_sFormula = Formula;
446  m_Formula = _Translate(Formula, "abcdefghijklmnopqrstuvwxyz", &m_Length, &m_Error_Position);
447 
448  if( m_Formula.code != NULL )
449  {
450  return( true );
451  }
452  }
453 
454  Destroy();
455 
456  return( false );
457 }
458 
459 
461 // //
463 
464 //---------------------------------------------------------
465 double CSG_Formula::Get_Value(void) const
466 {
467  return( _Get_Value(m_Parameters, m_Formula) );
468 }
469 
470 //---------------------------------------------------------
471 double CSG_Formula::Get_Value(double x) const
472 {
473  double Parameters[32];
474 
475  memcpy(Parameters, m_Parameters, 32 * sizeof(double));
476 
477  Parameters['x'-'a'] = x;
478 
479  return( _Get_Value(Parameters, m_Formula) );
480 }
481 
482 //---------------------------------------------------------
483 double CSG_Formula::Get_Value(const CSG_Vector &Values) const
484 {
485  return( Get_Value(Values.Get_Data(), Values.Get_N()) );
486 }
487 
488 //---------------------------------------------------------
489 double CSG_Formula::Get_Value(double *Values, int nValues) const
490 {
491  double Parameters[32];
492 
493  for(int i=0; i<nValues; i++)
494  {
495  Parameters[i] = Values[i];
496  }
497 
498  return( _Get_Value(Parameters, m_Formula) );
499 }
500 
501 //---------------------------------------------------------
502 double CSG_Formula::Get_Value(const char *Args, ...) const
503 {
504  double Parameters[32];
505 
506  va_list ap;
507 
508  va_start(ap, Args);
509 
510  while( *Args )
511  {
512  Parameters[(*Args++) - 'a'] = va_arg(ap, double);
513  }
514 
515  va_end(ap);
516 
517  return( _Get_Value(Parameters, m_Formula) );
518 }
519 
520 //---------------------------------------------------------
521 double CSG_Formula::_Get_Value(const double *Parameters, TSG_Formula func) const
522 {
523  double x, y, z, buffer[GET_VALUE_BUFSIZE];
524 
525  double *bufp = buffer; // points to the first free space in the buffer
526  char *function = func.code;
527  double *ctable = func.ctable;
528  double result;
529 
530  if( !function )
531  {
532  return( 0 ); // _Set_Error(_TL("empty coded function"));
533  }
534 
535  for( ; ; )
536  {
537  switch( *function++ )
538  {
539  case '\0':
540  return( buffer[0] );
541 
542  case 'D':
543  *bufp++ = ctable[*function++];
544  break;
545 
546  case 'V':
547  *bufp++ = Parameters[(*function++) - 'a'];
548  break;
549 
550  case 'M':
551  result = -(*--bufp);
552  *bufp++ = result;
553  break;
554 
555  case '+':
556  y = *(--bufp);
557  result = y + *(--bufp);
558  *bufp++ = result;
559  break;
560 
561  case '-':
562  y = *--bufp;
563  result = *(--bufp) - y;
564  *bufp++ = result;
565  break;
566 
567  case '*':
568  y = *(--bufp);
569  result = *(--bufp) * y;
570  *bufp++ = result;
571  break;
572 
573  case '/':
574  y = *--bufp;
575  result = *(--bufp) / y;
576  *bufp++ = result;
577  break;
578 
579  case '^':
580  y = *--bufp;
581  result = pow(*(--bufp), y);
582  *bufp++ = result;
583  break;
584 
585  case '=':
586  y = *--bufp;
587  result = y == *(--bufp) ? 1.0 : 0.0;
588  *bufp++ = result;
589  break;
590 
591  case '>':
592  y = *--bufp;
593  result = y < *(--bufp) ? 1.0 : 0.0;
594  *bufp++ = result;
595  break;
596 
597  case '<':
598  y = *--bufp;
599  result = y > *(--bufp) ? 1.0 : 0.0;
600  *bufp++ = result;
601  break;
602 
603  case '&':
604  y = *--bufp;
605  result = y && *(--bufp) ? 1.0 : 0.0;
606  *bufp++ = result;
607  break;
608 
609  case '|':
610  y = *--bufp;
611  result = y || *(--bufp) ? 1.0 : 0.0;
612  *bufp++ = result;
613  break;
614 
615  case 'F':
616  switch (m_Functions[*function].nParameters)
617  {
618  case 0:
619  *bufp++ = ((TSG_Formula_Function_0)m_Functions[*function++].Function)();
620  break;
621 
622  case 1:
623  x = *--bufp;
624  *bufp++ = ((TSG_Formula_Function_1)m_Functions[*function++].Function)(x);
625  break;
626 
627  case 2:
628  y = *--bufp;
629  x = *--bufp;
630  *bufp++ = ((TSG_Formula_Function_2)m_Functions[*function++].Function)(x, y);
631  break;
632 
633  case 3:
634  z = *--bufp;
635  y = *--bufp;
636  x = *--bufp;
637  *bufp++ = ((TSG_Formula_Function_3)m_Functions[*function++].Function)(x, y, z);
638  break;
639 
640  default:
641  return( 0 ); // _Set_Error(_TL("I2: too many parameters"));
642  }
643  break;
644 
645  default:
646  return( 0 ); // _Set_Error(_TL("I1: unrecognizable operator"));
647  }
648  }
649 
650 // if( (bufp - buffer) != 1 ) // _Set_Error(_TL("I3: corrupted buffer"));
651 
652  return( buffer[0] );
653 }
654 
655 
657 // //
658 // //
659 // //
661 
662 //---------------------------------------------------------
663 const char * CSG_Formula::Get_Used_Variables(void)
664 {
665  static CSG_String ret;
666 
667  ret.Clear();
668 
669  for(int i=0; i<'z'-'a'; i++)
670  {
671  if( m_Vars_Used[i] == true )
672  {
673  ret.Append((char)(i + 'a'));
674  }
675  }
676 
677  return( ret );
678 }
679 
680 
682 // //
683 // //
684 // //
686 
687 //---------------------------------------------------------
688 // int varying; Does the result of the function vary
689 // even when the parameters stay the same?
690 // varying = 1 for e.g. random - number generators.
691 // Result: 0 is rendered if there is an error
692 // 1 is rendered otherwise
693 //
694 bool CSG_Formula::Add_Function(const char *Name, TSG_Formula_Function_1 Function, int nParameters, bool bVarying)
695 {
696  if( nParameters < 0 || nParameters > 3 )
697  {
698  _Set_Error(_TL("invalid number of parameters"));
699 
700  return( false );
701  }
702 
703  TSG_Function *pFunction;
704 
705  for(pFunction=m_Functions; pFunction->Function && strcmp(Name, pFunction->Name); pFunction++)
706  {}
707 
708  if( pFunction->Function != NULL ) // old function is superseded
709  {
710  pFunction->Function = Function;
711  pFunction->nParameters = nParameters;
712  pFunction->bVarying = bVarying;
713 
714  _Set_Error();
715 
716  return( true );
717  }
718 
719  if( (pFunction - m_Functions) >= MAX_CTABLE - 1 )
720  {
721  _Set_Error(_TL("function table full"));
722 
723  return( false );
724  }
725 
726  pFunction->Name = Name;
727  pFunction->Function = Function;
728  pFunction->nParameters = nParameters;
729  pFunction->bVarying = bVarying;
730 
731  _Set_Error();
732 
733  return( true );
734 }
735 
736 
738 // //
740 
741 //---------------------------------------------------------
742 int CSG_Formula::_Get_Function(int i, char *Name, int *nParameters, int *bVarying)
743 {
744  if( !m_Functions[i].Function )
745  {
746  _Set_Error(_TL("index out of bounds"));
747 
748  return( 0 );
749  }
750 
751  strcpy(Name , m_Functions[i].Name);
752  *nParameters = m_Functions[i].nParameters;
753  *bVarying = m_Functions[i].bVarying ? 1 : 0;
754 
755  _Set_Error();
756 
757  return( 1 );
758 }
759 
760 //---------------------------------------------------------
761 // If the function exists, _Get_Function() returns the index
762 // of its name in the table. Otherwise, it returns -1.
763 //
764 int CSG_Formula::_Get_Function(const char *Name)
765 {
766  TSG_Function *pFunction;
767 
768  for(pFunction=m_Functions; pFunction->Function && strcmp(Name, pFunction->Name); pFunction++)
769  {}
770 
771  if( pFunction->Function == NULL )
772  {
773  _Set_Error(_TL("function not found"));
774 
775  return( -1 );
776  }
777 
778  _Set_Error();
779 
780  return( (int)(pFunction - m_Functions) );
781 }
782 
783 
785 // //
786 // //
787 // //
789 
790 //---------------------------------------------------------
791 inline int CSG_Formula::_is_Operand(char c)
792 {
793  return( (c == '+')
794  || (c == '-')
795  || (c == '*')
796  || (c == '/')
797  || (c == '^')
798  || (c == '=')
799  || (c == '<')
800  || (c == '>')
801  || (c == '&')
802  || (c == '|')
803  );
804 }
805 
806 //---------------------------------------------------------
807 inline int CSG_Formula::_is_Operand_Code(char c)
808 {
809  return( (c == '+')
810  || (c == '-')
811  || (c == '*')
812  || (c == '/')
813  || (c == '^')
814  || (c == '=')
815  || (c == '<')
816  || (c == '>')
817  || (c == '&')
818  || (c == '|')
819  || (c == 'M')
820  );
821 }
822 
823 //---------------------------------------------------------
824 inline int CSG_Formula::_is_Number(char c)
825 {
826  return( isdigit(c) || c == '.' || c == 'E' );
827 }
828 
829 
831 // //
832 // //
833 // //
835 
837 // Interpreting functions //
839 
840 //---------------------------------------------------------
841 CSG_Formula::TSG_Formula CSG_Formula::_Translate(const char *sourc, const char *args, int *leng, int *error)
842 {
843  const char *scan, *scarg;
844  char *result, *source, *code, *nfunc;
845  size_t size_estim;
846  double *ctable;
847  TSG_Formula returned;
848 
849  //-----------------------------------------------------
850  *leng = 0;
851  *error = 0;
852  m_error = NULL;
853  returned.code = NULL;
854  returned.ctable = NULL;
855 
856  //-----------------------------------------------------
857  source = (char *)SG_Malloc((strlen(sourc) + 1) * sizeof(char));
858 
859  if( source == NULL )
860  {
861  _Set_Error(_TL("no memory"));
862 
863  return( returned );
864  }
865 
866  strcpy(source, sourc);
867 
868  for(scan=source; *scan!='\0'; scan++)
869  {
870  if( islower(*scan) && !isalpha(*(scan + 1)) && (scan == source || !isalpha(*(scan - 1))) )
871  {
872  for(scarg=args; *scarg!='\0' && *scarg != *scan; scarg++)
873  {}
874 
875  if( *scarg == '\0' )
876  {
877  _Set_Error(_TL("undeclared parameter"));
878 
879  m_error = scan;
880  *error = (int)(m_error - source);
881 
882  SG_Free(source);
883 
884  return (returned);
885  }
886  }
887  }
888 
889  //-----------------------------------------------------
890  size_estim = _max_size(source);
891 
892  if( !(code = (char *)SG_Malloc(size_estim)) )
893  {
894  _Set_Error(_TL("no memory"));
895 
896  *error = -1;
897 
898  SG_Free(source);
899 
900  return (returned);
901  }
902 
903 
904  //-----------------------------------------------------
905  m_pctable = 0;
906 
907  if( !(m_ctable = (double *)SG_Malloc(MAX_CTABLE * sizeof(double))) )
908  {
909  _Set_Error(_TL("no memory"));
910 
911  *error = -1;
912 
913  SG_Free(source);
914  SG_Free(code);
915 
916  return (returned);
917  }
918 
919  ctable = m_ctable;
920 
921  //-----------------------------------------------------
922  _Set_Error();
923 
924  result = _i_trans(code, (char *)source, (char *)source + strlen(source));
925 
926  if( !result || m_bError )
927  {
928  *error = (int)(m_error ? m_error - source : -1);
929 
930  SG_Free(source);
931  SG_Free(code);
932  SG_Free(m_ctable);
933 
934  return (returned);
935  }
936  else
937  {
938  *result = '\0';
939  *error = -1;
940  *leng = (int)(result - code);
941 
942  if( ((*leng) + 1) * sizeof(char) > size_estim )
943  {
944  _Set_Error(_TL("I4: size estimate too small"));
945 
946  SG_Free(source);
947 
948  return( returned );
949  }
950  else if( ((*leng) + 1) * sizeof(char) < size_estim )
951  {
952  nfunc = (char *)SG_Malloc(((*leng) + 1) * sizeof(char));
953 
954  if (nfunc)
955  {
956  memcpy(nfunc, code, ((*leng) + 1) * sizeof(char));
957  SG_Free(code);
958  code = nfunc;
959  }
960  }
961 
962  if( m_pctable < MAX_CTABLE )
963  {
964  ctable = (double *)SG_Malloc(m_pctable * sizeof(double));
965 
966  if( ctable )
967  {
968  memcpy(ctable, m_ctable, m_pctable * sizeof(double));
969 
970  SG_Free(m_ctable);
971  }
972  else
973  {
974  ctable = m_ctable;
975  }
976  }
977  else
978  {
979  ctable = m_ctable;
980  }
981 
982  returned.code = code;
983  returned.ctable = ctable;
984 
985  _Set_Error();
986 
987  SG_Free(source);
988 
989  return (returned);
990  }
991 }
992 
993 
995 // //
996 // //
997 // //
999 
1000 //---------------------------------------------------------
1001 char * CSG_Formula::_i_trans(char *function, char *begin, char *end)
1002 {
1003  char tempch, *scan, *endf, *tempu, *temp3, *temps = NULL, *par_buf, *paramstr[MAX_PARMS];
1004  int i, pars, space, n_function;
1005  double tempd;
1006 
1007  if( begin >= end )
1008  {
1009  _Set_Error(_TL("missing operand"));
1010  m_error = begin;
1011  return NULL;
1012  }
1013 
1014  for(pars = 0, scan = begin; scan < end && pars >= 0; scan++)
1015  {
1016  if( *scan == '(' ) { pars++; } else
1017  if( *scan == ')' ) { pars--; }
1018  }
1019 
1020  if( pars < 0 || pars > 0 )
1021  {
1022  _Set_Error(_TL("unmatched parentheses"));
1023  m_error = scan - 1;
1024  return NULL;
1025  }
1026 
1027  for(pars = 0, scan = end - 1; scan >= begin; scan--)
1028  {
1029  if( *scan == '(' ) pars++; else
1030  if( *scan == ')' ) pars--; else
1031  if( !pars && (*scan == '+' || ((*scan == '-') && scan != begin)) && (scan == begin || *(scan - 1) != 'E') )
1032  break;
1033  }
1034 
1035  if( scan >= begin )
1036  {
1037  if ((tempu = _i_trans(function, begin, scan)) &&
1038  (temp3 = _i_trans(tempu, scan + 1, end)))
1039  {
1040  *temp3++ = *scan;
1041  temp3 = _comp_time(function, temp3, 2);
1042  if( m_bError )
1043  return NULL;
1044  else
1045  return temp3;
1046  }
1047  else
1048  return NULL;
1049  }
1050 
1051  for (pars = 0, scan = end - 1; scan >= begin; scan--)
1052  {
1053  if (*scan == '(') pars++;
1054  else if (*scan == ')')
1055  pars--;
1056  else if (!pars &&(*scan == '*' || *scan == '/'))
1057  break;
1058  }
1059  if (scan >= begin)
1060  {
1061  if ((tempu = _i_trans(function, begin, scan)) &&
1062  (temp3 = _i_trans(tempu, scan + 1, end)))
1063  {
1064  *temp3++ = *scan;
1065  temp3 = _comp_time(function, temp3, 2);
1066  if (m_bError)
1067  return NULL;
1068  else
1069  return temp3;
1070  }
1071  else
1072  return NULL;
1073  }
1074 
1075  /* unary minus */
1076  if (*begin == '-')
1077  {
1078  tempu = _i_trans(function, begin + 1, end);
1079  if (tempu)
1080  {
1081  *tempu++ = 'M';
1082  tempu = _comp_time(function, tempu, 1);
1083  if (m_bError)
1084  return NULL;
1085  else
1086  return tempu;
1087  }
1088  else
1089  return NULL;
1090  }
1091 
1092  for (pars = 0, scan = end - 1; scan >= begin; scan--)
1093  {
1094  if (*scan == '(') pars++;
1095  else if (*scan == ')')
1096  pars--;
1097  else if (!pars &&(*scan == '^'))
1098  break;
1099  else if (!pars &&(*scan == '='))
1100  break;
1101  else if (!pars &&(*scan == '>'))
1102  break;
1103  else if (!pars &&(*scan == '<'))
1104  break;
1105  else if (!pars &&(*scan == '&'))
1106  break;
1107  else if (!pars &&(*scan == '|'))
1108  break;
1109  }
1110 
1111  if (scan >= begin)
1112  {
1113  if ((tempu = _i_trans(function, begin, scan)) &&
1114  (temp3 = _i_trans(tempu, scan + 1, end)))
1115  {
1116  *temp3++ = *scan;
1117  temp3 = _comp_time(function, temp3, 2);
1118  if (m_bError)
1119  return NULL;
1120  else
1121  return temp3;
1122  }
1123  else
1124  return NULL;
1125  }
1126 
1127  /* erase white space */
1128  while (isspace(*begin))
1129  begin++;
1130  while (isspace(*(end - 1)))
1131  end--;
1132 
1133  if (*begin == '(' && *(end - 1) == ')')
1134  return _i_trans(function, begin + 1, end - 1);
1135 
1136  if (end == begin + 1 && islower(*begin))
1137  {
1138  *function++ = 'V';
1139  *function++ = *begin;
1140  return function;
1141  }
1142 
1143  tempch = *end;
1144  *end = '\0';
1145  tempd = strtod(begin, (char **)&tempu);
1146  *end = tempch;
1147 
1148  if( (char *)tempu == end )
1149  {
1150  *function++ = 'D';
1151  if (m_pctable < MAX_CTABLE)
1152  {
1153  m_ctable[m_pctable] = tempd;
1154  *function++ = (char)m_pctable++;
1155  }
1156  else
1157  {
1158  _Set_Error(_TL("too many constants"));
1159  m_error = begin;
1160  return NULL;
1161  }
1162  return function;
1163  }
1164 
1165  /*function*/
1166  if (!isalpha(*begin) && *begin != '_')
1167  {
1168  _Set_Error(_TL("syntax error"));
1169  m_error = begin;
1170  return NULL;
1171  }
1172 
1173  for(endf = begin + 1; endf < end &&(isalnum(*endf) || *endf == '_'); endf++)
1174  {}
1175 
1176  tempch = *endf;
1177  *endf = '\0';
1178  if ((n_function = _Get_Function(begin)) == -1)
1179  {
1180  *endf = tempch;
1181  m_error = begin;
1182  return NULL;
1183  }
1184  *endf = tempch;
1185  if (*endf != '(' || *(end - 1) != ')')
1186  {
1187  _Set_Error(_TL("improper function syntax"));
1188  m_error = endf;
1189  return NULL;
1190  }
1191 
1192  if( m_Functions[n_function].nParameters == 0 )
1193  {
1194  /*function without parameters(e.g. pi()) */
1195  space = 1;
1196  for (scan = endf + 1; scan <(end - 1); scan++)
1197  if (!isspace(*scan))
1198  space = 0;
1199  if (space)
1200  {
1201  *function++ = 'F';
1202  *function++ = n_function;
1203  function = _comp_time(function - 2, function, 0);
1204  if (m_bError)
1205  return NULL; /* internal error in _comp_time */
1206  else
1207  return function;
1208  }
1209  else
1210  {
1211  m_error = endf + 1;
1212  _Set_Error(_TL("too many parameters"));
1213  return NULL;
1214  }
1215  }
1216  else
1217  { /*function with parameters*/
1218  tempch = *(end - 1);
1219  *(end - 1) = '\0';
1220  par_buf = (char *)SG_Malloc(sizeof(char) * (strlen(endf + 1) + 1));
1221 
1222  if (!par_buf)
1223  {
1224  _Set_Error(_TL("no memory"));
1225  m_error = NULL;
1226  return NULL;
1227  }
1228 
1229  strcpy(par_buf, endf + 1);
1230  *(end - 1) = tempch;
1231 
1232  for (i = 0; i < m_Functions[n_function].nParameters; i++)
1233  {
1234  if ((temps = _my_strtok((i == 0) ? par_buf : NULL)) == NULL)
1235  break;
1236  paramstr[i] = temps;
1237  }
1238 
1239  if (temps == NULL)
1240  {
1241  SG_Free(par_buf);
1242  m_error = end - 2;
1243  _Set_Error(_TL("too few parameters"));
1244  return NULL;
1245  }
1246 
1247  if ((temps = _my_strtok(NULL)) != NULL)
1248  {
1249  SG_Free(par_buf);
1250  m_error =(temps - par_buf) +(endf + 1);
1251  _Set_Error(_TL("too many parameters"));
1252  return NULL;
1253  }
1254 
1255  tempu = function;
1256  for (i = 0; i < m_Functions[n_function].nParameters; i++)
1257  if( !(tempu = _i_trans(tempu, paramstr[i], paramstr[i] + strlen(paramstr[i]))) )
1258  {
1259  m_error =(m_error - par_buf) +(endf + 1);
1260  SG_Free(par_buf);
1261 
1262  return NULL;
1263  }
1264 
1265  /* OK */
1266  SG_Free(par_buf);
1267  *tempu++ = 'F';
1268  *tempu++ = n_function;
1269  tempu = _comp_time(function, tempu, m_Functions[n_function].nParameters);
1270  if (m_bError)
1271  return NULL; /* internal error in _comp_time */
1272  else
1273  return tempu;
1274  }
1275 }
1276 
1277 //---------------------------------------------------------
1278 char * CSG_Formula::_comp_time(char *function, char *fend, int npars)
1279 {
1280  char *scan, temp;
1281  int i;
1282  double tempd;
1283  TSG_Formula trans;
1284 
1285  scan = function;
1286  for (i = 0; i < npars; i++)
1287  {
1288  if (*scan++ != 'D')
1289  return fend;
1290  scan++;
1291  }
1292 
1293  if (!((scan == fend -(sizeof((char) 'F') + sizeof(char))
1294  && *(fend - 2) == 'F' && m_Functions[*(fend - 1)].bVarying == 0) ||
1295  (scan == fend - sizeof(char)
1296  && _is_Operand_Code(*(fend - 1))))
1297  )
1298  return fend;
1299 
1300  temp = *fend;
1301  *fend = '\0';
1302 
1303  trans.code = function;
1304  trans.ctable = m_ctable;
1305  tempd = _Get_Value(m_Parameters, trans);
1306  *fend = temp;
1307  *function++ = 'D';
1308  m_pctable -= npars;
1309  *function++ =(char) m_pctable;
1310  m_ctable[m_pctable++] = tempd;
1311 
1312  return function;
1313 }
1314 
1315 //---------------------------------------------------------
1316 int CSG_Formula::_max_size(const char *source)
1317 {
1318  int numbers = 0;
1319  int functions = 0;
1320  int operators = 0;
1321  int variables = 0;
1322 
1323  const size_t var_size = 2 * sizeof(char);
1324  const size_t num_size = sizeof(char) + sizeof(double);
1325  const size_t op_size = sizeof(char);
1326  const size_t end_size = sizeof('\0');
1327 
1328  for(int i=0; i<'z'-'a'; i++)
1329  {
1330  m_Vars_Used[i] = false;
1331  }
1332 
1333  for(const char *scan=source; *scan; scan++)
1334  {
1335  if( isalpha(*scan) && (*scan != 'E') )
1336  {
1337  if( isalpha(*(scan + 1)) || isdigit(*(scan + 1)) )
1338  {
1339  // must be a function name (combination of letters and digits, e.g. sin(..), atan2(..))
1340  }
1341  else if( *(scan + 1) == '(' )
1342  {
1343  functions++;
1344  }
1345  else
1346  {
1347  variables++;
1348  m_Vars_Used[(int)(*scan - 'a')] = true;
1349  }
1350  }
1351  }
1352 
1353  if( _is_Operand(*source) )
1354  {
1355  operators++;
1356  }
1357 
1358  if( *source != '\0' )
1359  {
1360  for(const char *scan=source+1; *scan; scan++)
1361  {
1362  if( _is_Operand(*scan) && *(scan - 1) != 'E' )
1363  {
1364  operators++;
1365  }
1366  }
1367  }
1368 
1369  const char *scan = source;
1370 
1371  while( *scan )
1372  {
1373  if( _is_Number(*scan) || ((*scan == '+' || *scan == '-') && scan > source && *(scan - 1) == 'E') )
1374  {
1375  numbers++;
1376  scan++;
1377 
1378  while( _is_Number(*scan) ||((*scan == '+' || *scan == '-') && scan > source && *(scan - 1) == 'E') )
1379  scan++;
1380  }
1381  else
1382  {
1383  scan++;
1384  }
1385  }
1386 
1387  return( numbers*num_size + operators*op_size + functions*num_size + variables*var_size + end_size );
1388 }
1389 
1390 //---------------------------------------------------------
1391 char * CSG_Formula::_my_strtok(char *s)
1392 {
1393  static char *token = NULL;
1394 
1395  if( s != NULL )
1396  token = s;
1397  else if( token != NULL )
1398  s = token;
1399  else
1400  return( NULL );
1401 
1402  for(int pars=0; *s != '\0' && (*s != ',' || pars != 0); s++)
1403  {
1404  if (*s == '(') ++pars;
1405  if (*s == ')') --pars;
1406  }
1407 
1408  if( *s == '\0' )
1409  {
1410  s = token;
1411  token = NULL;
1412  }
1413  else
1414  {
1415  *s = '\0';
1416  char *next_token = s + 1;
1417  s = token;
1418  token = next_token;
1419  }
1420 
1421  return( s );
1422 }
1423 
1424 
1426 // //
1427 // //
1428 // //
1430 
1431 //---------------------------------------------------------
CSG_Formula::Add_Function
bool Add_Function(const char *Name, TSG_Formula_Function_1 Function, int nParameters, bool bVarying=false)
Definition: mat_formula.cpp:694
CSG_Formula::Get_Error
bool Get_Error(CSG_String &Message)
Definition: mat_formula.cpp:385
SG_FREE_SAFE
#define SG_FREE_SAFE(PTR)
Definition: api_core.h:205
CSG_Formula::SSG_Function::Function
TSG_Formula_Function_1 Function
Definition: mat_tools.h:1897
CSG_Formula::CSG_Formula
CSG_Formula(void)
Definition: mat_formula.cpp:246
gSG_Functions
static CSG_Formula::TSG_Function gSG_Functions[MAX_CTABLE]
Definition: mat_formula.cpp:206
CSG_String::Append
CSG_String & Append(const CSG_String &String)
Definition: api_string.cpp:355
_TL
#define _TL(s)
Definition: api_core.h:1489
f_fmod
static double f_fmod(double x, double val)
Definition: mat_formula.cpp:162
f_rand_g
static double f_rand_g(double mean, double stdv)
Definition: mat_formula.cpp:174
CSG_Vector::Get_Data
double * Get_Data(void) const
Definition: mat_tools.h:384
CSG_String::Length
size_t Length(void) const
Definition: api_string.cpp:172
f_pi
static double f_pi(void)
Definition: mat_formula.cpp:144
CSG_Random::Get_Uniform
static double Get_Uniform(void)
Definition: mat_tools.cpp:274
MAX_CTABLE
#define MAX_CTABLE
Definition: mat_formula.cpp:84
SG_Malloc
SAGA_API_DLL_EXPORT void * SG_Malloc(size_t size)
Definition: api_memory.cpp:65
CSG_Random::Get_Gaussian
static double Get_Gaussian(double mean, double stddev)
Definition: mat_tools.cpp:298
CSG_Formula::SSG_Function
Definition: mat_tools.h:1894
MAX_PARMS
#define MAX_PARMS
Definition: mat_formula.cpp:83
grid.h
CSG_Formula::Set_Formula
bool Set_Formula(const CSG_String &Formula)
Definition: mat_formula.cpp:439
SG_Free
SAGA_API_DLL_EXPORT void SG_Free(void *memblock)
Definition: api_memory.cpp:83
f_sqr
static double f_sqr(double x)
Definition: mat_formula.cpp:156
f_rand_u
static double f_rand_u(double min, double max)
Definition: mat_formula.cpp:168
f_min
static double f_min(double a, double b)
Definition: mat_formula.cpp:132
CSG_String::Replace
size_t Replace(const CSG_String &sOld, const CSG_String &sNew, bool bReplaceAll=true)
Definition: api_string.cpp:563
SG_Calloc
SAGA_API_DLL_EXPORT void * SG_Calloc(size_t num, size_t size)
Definition: api_memory.cpp:71
f_ifelse
static double f_ifelse(double condition, double x, double y)
Definition: mat_formula.cpp:192
mat_tools.h
f_pow
static double f_pow(double x, double val)
Definition: mat_formula.cpp:108
CSG_Vector
Definition: mat_tools.h:360
CSG_Formula::Get_Help_Operators
static CSG_String Get_Help_Operators(bool bHTML=true, const CSG_String Additional[][2]=NULL)
Definition: mat_formula.cpp:295
CSG_Formula::Set_Variable
void Set_Variable(char Variable, double Value)
Definition: mat_formula.cpp:433
CSG_Formula::Get_Used_Variables
const char * Get_Used_Variables(void)
Definition: mat_formula.cpp:663
TSG_Formula_Function_0
double(* TSG_Formula_Function_0)(void)
Definition: mat_tools.h:1858
CSG_Vector::Get_N
int Get_N(void) const
Definition: mat_tools.h:382
CSG_Formula::Destroy
bool Destroy(void)
Definition: mat_formula.cpp:277
CSG_Formula::SSG_Function::nParameters
int nParameters
Definition: mat_tools.h:1899
f_eq
static double f_eq(double x, double val)
Definition: mat_formula.cpp:126
f_int
static double f_int(double x)
Definition: mat_formula.cpp:150
CSG_String::Format
static CSG_String Format(const char *Format,...)
Definition: api_string.cpp:270
f_or
static double f_or(double x, double y)
Definition: mat_formula.cpp:186
CSG_String::Left
CSG_String Left(size_t count) const
Definition: api_string.cpp:705
CSG_Formula::Get_Value
double Get_Value(void) const
Definition: mat_formula.cpp:465
CSG_String::Clear
void Clear(void)
Definition: api_string.cpp:259
TSG_Formula_Function_2
double(* TSG_Formula_Function_2)(double, double)
Definition: mat_tools.h:1860
M_PI
#define M_PI
Definition: mat_tools.h:98
CSG_String
Definition: api_core.h:563
CSG_Formula::SSG_Function::Name
const char * Name
Definition: mat_tools.h:1895
TSG_Formula_Function_1
double(* TSG_Formula_Function_1)(double)
Definition: mat_tools.h:1859
CSG_String::is_Empty
bool is_Empty(void) const
Definition: api_string.cpp:178
f_and
static double f_and(double x, double y)
Definition: mat_formula.cpp:180
CSG_Formula::~CSG_Formula
virtual ~CSG_Formula(void)
Definition: mat_formula.cpp:269
CSG_Formula::TSG_Function
struct CSG_Formula::SSG_Function TSG_Function
CSG_Formula::SSG_Function::bVarying
bool bVarying
Definition: mat_tools.h:1901
TSG_Formula_Function_3
double(* TSG_Formula_Function_3)(double, double, double)
Definition: mat_tools.h:1861
GET_VALUE_BUFSIZE
#define GET_VALUE_BUFSIZE
Definition: mat_formula.cpp:89
f_lt
static double f_lt(double x, double val)
Definition: mat_formula.cpp:120
f_gt
static double f_gt(double x, double val)
Definition: mat_formula.cpp:114
f_max
static double f_max(double a, double b)
Definition: mat_formula.cpp:138
f_atan2
static double f_atan2(double x, double val)
Definition: mat_formula.cpp:102
EPSILON
#define EPSILON
Definition: mat_formula.cpp:92
CSG_String::Right
CSG_String Right(size_t count) const
Definition: api_string.cpp:693