inmost_autodiff.h 38.6 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1

Kirill Terekhov's avatar
Kirill Terekhov committed
2 3 4 5
#ifndef INMOST_AUTODIFF_H_INCLUDED
#define INMOST_AUTODIFF_H_INCLUDED
#include "inmost_common.h"
#include "inmost_mesh.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
6
#include "inmost_solver.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
7
#include "inmost_variable.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
8 9
#include <sstream> //for debug

Kirill Terekhov's avatar
Kirill Terekhov committed
10
//#define NEW_VERSION
Kirill Terekhov's avatar
Kirill Terekhov committed
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

#if defined(USE_AUTODIFF) && (!defined(USE_MESH))
#warning "USE_AUTODIFF require USE_MESH"
#undef USE_AUTODIFF
#endif

//#define DPRNT



#if defined(USE_AUTODIFF)
#include <math.h>


#define FILTER_EPS 1e-12

#define AD_SPACE 16384 //typical number of registered entries

#define AD_NONE  0 //this should generate an error
//binary operations below
#define AD_PLUS  1 //sum x+y
#define AD_MINUS 2 //substraction x-y
#define AD_MULT  3 //multiplication x*y
#define AD_DIV   4 //devision x/y
Kirill Terekhov's avatar
Kirill Terekhov committed
35
#define AD_INV   5 //inversion 1/x, depricated
Kirill Terekhov's avatar
Kirill Terekhov committed
36
#define AD_POW   6 //power operation x^y
Kirill Terekhov's avatar
Kirill Terekhov committed
37 38
#define AD_SQRT  7 //square root operation

Kirill Terekhov's avatar
Kirill Terekhov committed
39 40 41

#define AD_PRECOMP 15 // this expression points to precomputed replacement, expr * left is position in array of precomputed values, expr * right points to original expression

Kirill Terekhov's avatar
Kirill Terekhov committed
42 43
#define AD_VAL   19 // calculate only value of current expression

Kirill Terekhov's avatar
Kirill Terekhov committed
44 45 46 47 48 49 50 51 52 53 54
//unary operations below
#define AD_COS   20 //cos(x)

#define AD_ABS   22 // |x|
#define AD_EXP   23 // exp(x)
#define AD_LOG   24 // log(x)
#define AD_SIN   25 // sin(x)

#define AD_CONST 50 // expr * left represents const value
#define AD_MES   51 // volume(e) for cell, area(e) for face, length(e) for edge

Kirill Terekhov's avatar
Kirill Terekhov committed
55 56 57 58 59 60 61
#define AD_COND  60 //condition
#define AD_COND_TYPE 61 //condition on element type
#define AD_COND_MARK 62 //condition on marker
#define AD_ALTR  65 //alternatives for condition



Kirill Terekhov's avatar
Kirill Terekhov committed
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
#define AD_EXT   99 // pointer to external expression data

#define AD_TAG   100               //tag with dynamic data
#define AD_CTAG  (AD_TAG+AD_SPACE)   //tag with constant data
#define AD_STNCL (AD_CTAG+AD_SPACE)  //stencil that sets current elements and their coefs
#define AD_TABLE (AD_STNCL+AD_SPACE) //table of values
#define AD_FUNC  (AD_TABLE+AD_SPACE) //register function that calculates some value from element


//types of operands

#define AD_TYPE_INVALID   ENUMUNDEF
#define AD_TYPE_ENDPOINT  0
#define AD_TYPE_UNARY     1
#define AD_TYPE_BINARY    2
#define AD_TYPE_TERNARY   3
Kirill Terekhov's avatar
Kirill Terekhov committed
78
#define AD_TYPE_VALUE     4
Kirill Terekhov's avatar
Kirill Terekhov committed
79 80 81 82 83 84 85 86

namespace INMOST
{
	class Automatizator; //forward declaration
	typedef dynarray<INMOST_DATA_REAL_TYPE, 2048> precomp_values_t;
	class expr;
	
	
Kirill Terekhov's avatar
Kirill Terekhov committed
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
  class Residual
  {
    Sparse::Matrix jacobian;
    Sparse::Vector residual;
  public:
    Residual(std::string name = "", INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0, INMOST_MPI_Comm _comm = INMOST_MPI_COMM_WORLD)
      : jacobian(name,start,end,_comm),residual(name,start,end,_comm) {}
    Residual(const Residual & other)
      : jacobian(other.jacobian), residual(other.residual)
    {}
    Residual & operator =(Residual const & other)
    {
      jacobian = other.jacobian;
      residual = other.residual;
      return *this;
    }
    Sparse::Matrix & GetJacobian() {return jacobian;}
    const Sparse::Matrix & GetJacobian() const {return jacobian;}
    Sparse::Vector & GetResidual() {return residual;}
    const Sparse::Vector & GetResidual() const {return residual;}
    INMOST_DATA_ENUM_TYPE GetFirstIndex() const {return residual.GetFirstIndex();}
    INMOST_DATA_ENUM_TYPE GetLastIndex() const {return residual.GetLastIndex();}
    void GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const 
    {
      start = residual.GetFirstIndex(); 
      end = residual.GetLastIndex();
    }
    void SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end)
    {
      jacobian.SetInterval(beg,end);
      residual.SetInterval(beg,end);
    }
    multivar_expression_reference operator [](INMOST_DATA_ENUM_TYPE row)
    {
      return multivar_expression_reference(residual[row],&jacobian[row]);
    }
    void ClearResidual()
    {
      for(Sparse::Vector::iterator it = residual.Begin(); it != residual.End(); ++it) (*it) = 0.0;
    }
    void ClearJacobian()
    {
      for(Sparse::Matrix::iterator it = jacobian.Begin(); it != jacobian.End(); ++it)
        it->Clear();
    }
    void Clear()
    {
#if defined(USE_OMP)
#pragma omp for
#endif
      for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k) 
      {
        residual[k] = 0.0;
        jacobian[k].Clear();
      }
    }
    INMOST_DATA_REAL_TYPE Norm()
    {
      INMOST_DATA_REAL_TYPE ret = 0;
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
147
#pragma omp for reduction(+:ret)
Kirill Terekhov's avatar
Kirill Terekhov committed
148 149 150
#endif
      for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k) 
        ret += residual[k]*residual[k];
Kirill Terekhov's avatar
Kirill Terekhov committed
151 152 153 154
#if defined(USE_MPI)
        INMOST_DATA_REAL_TYPE tmp = ret;
        MPI_Allreduce(&tmp, &ret, 1, INMOST_MPI_DATA_REAL_TYPE, MPI_SUM, jacobian.GetCommunicator());
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178
      return sqrt(ret);
    }
    /// Normalize entries in jacobian and right hand side
    void Rescale()
    {
#if defined(USE_OMP)
#pragma omp for
#endif
      for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k)
      {
        INMOST_DATA_REAL_TYPE norm = 0.0;
        for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q)
          norm += jacobian[k].GetValue(q)*jacobian[k].GetValue(q);
        norm = sqrt(norm);
        if( norm )
        {
          norm = 1.0/norm;
          residual[k] *= norm;
          for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q)
            jacobian[k].GetValue(q) *= norm;
        }
      }
    }
  };
Kirill Terekhov's avatar
Kirill Terekhov committed
179 180 181 182

	class Automatizator
	{
	public:
Kirill Terekhov's avatar
Kirill Terekhov committed
183 184
		typedef INMOST_DATA_REAL_TYPE(*func_callback)(const Storage & current_element, void * user_data);
		typedef std::pair<HandleType, INMOST_DATA_REAL_TYPE> stencil_pair;
Kirill Terekhov's avatar
Kirill Terekhov committed
185
		typedef dynarray<stencil_pair, 64> stencil_pairs;
Kirill Terekhov's avatar
Kirill Terekhov committed
186
		typedef void(*stencil_callback)(const Storage & current_element, stencil_pairs & out_stencil, void * user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
187
	private:
Kirill Terekhov's avatar
Kirill Terekhov committed
188
    static Automatizator * CurrentAutomatizator;
Kirill Terekhov's avatar
Kirill Terekhov committed
189 190 191 192 193
#if defined(USE_OMP)
    std::vector<Sparse::RowMerger> merger;
#else
    Sparse::RowMerger merger;
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
194
		typedef struct{ Tag t; MarkerType domain_mask; } tagdomain;
Kirill Terekhov's avatar
Kirill Terekhov committed
195
		typedef dynarray<tagdomain, 128> const_tag_type;
Kirill Terekhov's avatar
Kirill Terekhov committed
196
		typedef struct{ tagdomain d; Tag indices; } tagpair;
Kirill Terekhov's avatar
Kirill Terekhov committed
197
		typedef dynarray<tagpair, 128> tagpairs_type;
Kirill Terekhov's avatar
Kirill Terekhov committed
198 199
		typedef std::vector<tagpair> index_enum;
		typedef struct { Tag elements, coefs; } stencil_tag;
Kirill Terekhov's avatar
Kirill Terekhov committed
200
		typedef struct { std::string name; INMOST_DATA_ENUM_TYPE kind; MarkerType domainmask; void * link; } stencil_kind_domain;
Kirill Terekhov's avatar
Kirill Terekhov committed
201
		typedef dynarray<stencil_kind_domain, 128> stencil_type;
Kirill Terekhov's avatar
Kirill Terekhov committed
202
		typedef struct func_name_callback_t { std::string name; func_callback func; } func_name_callback;
Kirill Terekhov's avatar
Kirill Terekhov committed
203
		typedef dynarray<func_name_callback, 128> func_type;
Kirill Terekhov's avatar
Kirill Terekhov committed
204 205 206 207 208 209 210
	public:
		typedef struct
		{
			std::string name;
			INMOST_DATA_REAL_TYPE * args;
			INMOST_DATA_REAL_TYPE * vals;
			INMOST_DATA_ENUM_TYPE size;
Kirill Terekhov's avatar
Kirill Terekhov committed
211 212 213 214
			INMOST_DATA_ENUM_TYPE binary_search(INMOST_DATA_REAL_TYPE arg);
			INMOST_DATA_REAL_TYPE get_value(INMOST_DATA_REAL_TYPE arg);
			INMOST_DATA_REAL_TYPE get_derivative(INMOST_DATA_REAL_TYPE arg);
			std::pair<INMOST_DATA_REAL_TYPE, INMOST_DATA_REAL_TYPE> get_both(INMOST_DATA_REAL_TYPE arg);
Kirill Terekhov's avatar
Kirill Terekhov committed
215 216 217 218
		} table;
		typedef table * table_ptr;
	private:
		typedef dynarray<INMOST_DATA_REAL_TYPE, 128> values_container;
Kirill Terekhov's avatar
Kirill Terekhov committed
219 220 221 222 223 224 225
		typedef dynarray<table_ptr, 32>              table_type;
		const_tag_type        reg_ctags;
		index_enum            index_tags;
		tagpairs_type         reg_tags;
		table_type            reg_tables;
		stencil_type          reg_stencils;
		func_type             reg_funcs;
Kirill Terekhov's avatar
Kirill Terekhov committed
226 227 228
		INMOST_DATA_ENUM_TYPE first_num;
		INMOST_DATA_ENUM_TYPE last_num;
		Mesh * m;
Kirill Terekhov's avatar
Kirill Terekhov committed
229
#if defined(NEW_VERSION)
Kirill Terekhov's avatar
Kirill Terekhov committed
230
		INMOST_DATA_REAL_TYPE                                            EvaluateSub(expr & var,INMOST_DATA_ENUM_TYPE element,INMOST_DATA_ENUM_TYPE parent, void * user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
231
		void                                                             DerivativeFill(expr & var, INMOST_DATA_ENUM_TYPE element, INMOST_DATA_ENUM_TYPE parent, Sparse::RowMerger & entries, INMOST_DATA_REAL_TYPE multval, void * user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
232
#else
Kirill Terekhov's avatar
Kirill Terekhov committed
233
		INMOST_DATA_REAL_TYPE                                            DerivativePrecompute(const expr & var, const Storage & e, precomp_values_t & values, void * user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
234
		void                                                             DerivativeFill(const expr & var, const Storage & e, Sparse::RowMerger & entries, precomp_values_t & values, INMOST_DATA_REAL_TYPE multval, void * user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
235
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
236 237 238
	public:
		Automatizator(Mesh * m);
		~Automatizator();
Kirill Terekhov's avatar
Kirill Terekhov committed
239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257
		__INLINE INMOST_DATA_ENUM_TYPE                                   GetFirstIndex() { return first_num; }
		__INLINE INMOST_DATA_ENUM_TYPE                                   GetLastIndex() { return last_num; }
		INMOST_DATA_ENUM_TYPE                                            RegisterFunc(std::string name, func_callback func);
		INMOST_DATA_ENUM_TYPE                                            RegisterStencil(std::string name, Tag elements_tag, Tag coefs_tag, MarkerType domain_mask = 0);
		INMOST_DATA_ENUM_TYPE                                            RegisterStencil(std::string name, stencil_callback func, MarkerType domain_mask = 0);
		INMOST_DATA_ENUM_TYPE                                            RegisterTable(std::string name, INMOST_DATA_REAL_TYPE * Arguments, INMOST_DATA_REAL_TYPE * Values, INMOST_DATA_ENUM_TYPE size);
		INMOST_DATA_ENUM_TYPE                                            RegisterDynamicTag(Tag t, ElementType typemask, MarkerType domain_mask = 0);
		INMOST_DATA_ENUM_TYPE                                            RegisterStaticTag(Tag t, MarkerType domain_mask = 0);
		void                                                             EnumerateDynamicTags();
		__INLINE Tag                                                     GetDynamicValueTag(INMOST_DATA_ENUM_TYPE ind) { return reg_tags[ind-AD_TAG].d.t; }
		__INLINE Tag                                                     GetDynamicIndexTag(INMOST_DATA_ENUM_TYPE ind) { return reg_tags[ind-AD_TAG].indices; }
		__INLINE MarkerType                                              GetDynamicMask(INMOST_DATA_ENUM_TYPE ind) { return reg_tags[ind-AD_TAG].d.domain_mask; }
		__INLINE Tag                                                     GetStaticValueTag(INMOST_DATA_ENUM_TYPE ind) { return reg_ctags[ind-AD_CTAG].t; }
		__INLINE MarkerType                                              GetStaticMask(INMOST_DATA_ENUM_TYPE ind) { return reg_ctags[ind-AD_CTAG].domain_mask; }
		__INLINE INMOST_DATA_REAL_TYPE                                   GetDynamicValue(const Storage & e, INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_ENUM_TYPE comp = 0) { return e->RealArray(GetDynamicValueTag(ind))[comp]; }
		__INLINE INMOST_DATA_ENUM_TYPE                                   GetDynamicIndex(const Storage & e, INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_ENUM_TYPE comp = 0) { return e->IntegerArray(GetDynamicIndexTag(ind))[comp]; }
		__INLINE bool                                                    isDynamicValid(const Storage & e, INMOST_DATA_ENUM_TYPE ind) { MarkerType mask = GetDynamicMask(ind); return mask == 0 || e->GetMarker(mask); }
		__INLINE INMOST_DATA_REAL_TYPE                                   GetStaticValue(const Storage & e, INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_ENUM_TYPE comp = 0) { return e->RealArray(GetStaticValueTag(ind))[comp]; }
		__INLINE bool                                                    isStaticValid(const Storage & e, INMOST_DATA_ENUM_TYPE ind) { MarkerType mask = GetStaticMask(ind); return mask == 0 || e->GetMarker(mask); }
Kirill Terekhov's avatar
Kirill Terekhov committed
258
#if defined(NEW_VERSION)
Kirill Terekhov's avatar
Kirill Terekhov committed
259
		INMOST_DATA_REAL_TYPE                                            Evaluate(expr & var, const Storage & e, void * user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
260
		INMOST_DATA_REAL_TYPE                                            Derivative(expr & var, const Storage & e, Sparse::Row & out, Storage::real multiply, void * user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
261
#else
Kirill Terekhov's avatar
Kirill Terekhov committed
262
		INMOST_DATA_REAL_TYPE                                            Evaluate(const expr & var, const Storage & e, void * user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
263
		INMOST_DATA_REAL_TYPE                                            Derivative(const expr & var, const Storage & e, Sparse::Row & out, Storage::real multiply,  void * user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
264
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
265 266
		__INLINE INMOST_DATA_REAL_TYPE                                   GetIndex(const Storage & e, INMOST_DATA_ENUM_TYPE tagind, INMOST_DATA_ENUM_TYPE comp = 0) { return e->IntegerArray(GetDynamicIndexTag(tagind))[comp]; }
		__INLINE INMOST_DATA_ENUM_TYPE                                   GetComponents(const Storage & e, INMOST_DATA_ENUM_TYPE tagind) { return static_cast<INMOST_DATA_ENUM_TYPE>(e->IntegerArray(GetDynamicIndexTag(tagind)).size()); }
Kirill Terekhov's avatar
Kirill Terekhov committed
267
		__INLINE Mesh *                                                  GetMesh() { return m; }
Kirill Terekhov's avatar
Kirill Terekhov committed
268 269 270 271 272 273 274 275
		__INLINE INMOST_DATA_REAL_TYPE *                                 GetTableArguments(INMOST_DATA_ENUM_TYPE tableind) {return reg_tables[tableind-AD_TABLE]->args;}
		__INLINE INMOST_DATA_REAL_TYPE *                                 GetTableValues(INMOST_DATA_ENUM_TYPE tableind) {return reg_tables[tableind-AD_TABLE]->vals;}
		__INLINE INMOST_DATA_ENUM_TYPE                                   GetTableSize(INMOST_DATA_ENUM_TYPE tableind) {return reg_tables[tableind-AD_TABLE]->size;}
		__INLINE INMOST_DATA_REAL_TYPE                                   GetTableValue(INMOST_DATA_ENUM_TYPE tableind, INMOST_DATA_REAL_TYPE arg) { return reg_tables[tableind-AD_TABLE]->get_value(arg); }
		__INLINE INMOST_DATA_REAL_TYPE                                   GetTableDerivative(INMOST_DATA_ENUM_TYPE tableind, INMOST_DATA_REAL_TYPE arg) { return reg_tables[tableind-AD_TABLE]->get_derivative(arg); }
		__INLINE std::pair<INMOST_DATA_REAL_TYPE, INMOST_DATA_REAL_TYPE> GetTableBoth(INMOST_DATA_ENUM_TYPE tableind, INMOST_DATA_REAL_TYPE arg) { return reg_tables[tableind-AD_TABLE]->get_both(arg); }
		INMOST_DATA_ENUM_TYPE                                            GetStencil(INMOST_DATA_ENUM_TYPE stnclind, const Storage & elem, void * user_data, stencil_pairs & ret);
		__INLINE INMOST_DATA_REAL_TYPE                                   GetFunction(INMOST_DATA_ENUM_TYPE funcid, const Storage & elem, void * user_data) { return reg_funcs[funcid-AD_FUNC].func(elem, user_data); }
Kirill Terekhov's avatar
Kirill Terekhov committed
276 277 278 279 280 281 282 283
    Sparse::RowMerger &                                              GetMerger() 
    {
#if defined(USE_OMP)
      return merger[omp_get_thread_num()];
#else
      return merger;
#endif
    }
Kirill Terekhov's avatar
Kirill Terekhov committed
284 285 286 287 288 289 290 291
    /// Remove global current automatizator.
    static void RemoveCurrent() {CurrentAutomatizator = NULL;}
    /// Set current global automatizator, so that variable will be optimized with row merger.
    static void MakeCurrent(Automatizator * aut) {CurrentAutomatizator = aut;}
    /// Check that there is an automatizator.
    static bool HaveCurrent() {return CurrentAutomatizator != NULL;}
    /// Retrive the automatizator.
    static Automatizator * GetCurrent() {return CurrentAutomatizator;}
Kirill Terekhov's avatar
Kirill Terekhov committed
292 293
	};

Kirill Terekhov's avatar
Kirill Terekhov committed
294
#if defined(NEW_VERSION)
Kirill Terekhov's avatar
Kirill Terekhov committed
295 296 297 298 299 300 301 302
	class expr
	{
		
		typedef small_hash<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE, 1024> replace_type;
		
		class expr_data
		{
			typedef union { INMOST_DATA_REAL_TYPE r; INMOST_DATA_ENUM_TYPE i; expr * e; expr_data * q; } operand;
Kirill Terekhov's avatar
Kirill Terekhov committed
303
			INMOST_DATA_ENUM_TYPE op, priority;
Kirill Terekhov's avatar
Kirill Terekhov committed
304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350
			operand left, right;
		public:
			expr_data() : op(AD_NONE) { memset(&left, 0, sizeof(operand)); memset(&right, 0, sizeof(operand)); }
			expr_data(const expr_data & other) : op(other.op) 
			{
				if (other.op == AD_STNCL)
					left.e = new expr(*other.left.e);
				else if (other.op == AD_ALTR)
				{
					left.e = new expr(*other.left.e);
					right.e = new expr(*other.right.e);
				}
				else if (other.op == AD_COND)
				{
					left = other.left;
					right.q = new expr_data(*other.right.q);
				}
				else
				{
					left = other.left; 
					right = other.right;
				}
			}
			expr_data & operator =(expr_data const & other) 
			{ 
				op = other.op; 
				if (other.op == AD_STNCL)
					left.e = new expr(*other.left.e);
				else if (other.op == AD_ALTR)
				{
					left.e = new expr(*other.left.e);
					right.e = new expr(*other.right.e);
				}
				else if (other.op == AD_COND)
				{
					left = other.left;
					right.q = new expr_data(*other.right.q);
				}
				else
				{
					left = other.left;
					right = other.right;
				}
				return *this;
			}
			expr_data(INMOST_DATA_REAL_TYPE val) : op(AD_CONST) { left.r = val; }
			expr_data(INMOST_DATA_ENUM_TYPE new_op, INMOST_DATA_ENUM_TYPE l, INMOST_DATA_ENUM_TYPE r) : op(new_op) { left.i = l; right.i = r; }
Kirill Terekhov's avatar
Kirill Terekhov committed
351 352
			expr_data(INMOST_DATA_ENUM_TYPE new_op, INMOST_DATA_ENUM_TYPE l, INMOST_DATA_REAL_TYPE r) : op(new_op) { left.i = l; right.r = r; }
			expr_data(INMOST_DATA_ENUM_TYPE new_op, INMOST_DATA_ENUM_TYPE l, const expr & e) : op(new_op) { left.i = l; right.e = new expr(e); }
Kirill Terekhov's avatar
Kirill Terekhov committed
353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395
			expr_data(INMOST_DATA_ENUM_TYPE l, const expr_data & e) : op(AD_COND) { left.i = l; right.q = new expr_data(e); }
			expr_data(INMOST_DATA_ENUM_TYPE op, const expr * e) : op(op) { left.e = new expr(*e); }
			expr_data(const expr * a, const expr * b) : op(AD_ALTR) { left.e = new expr(*a); right.e = new expr(*b); }
			~expr_data()
			{
				if (op == AD_STNCL)
				{
					delete left.e;
					op = AD_NONE;
				}
				else if (op == AD_COND)
				{
					delete right.q;
					op = AD_NONE;
				}
				else if (op == AD_ALTR)
				{
					delete left.e;
					delete right.e;
					op = AD_NONE;
				}
			}
			bool operator ==(const expr_data & other) const
			{
				assert(op != AD_NONE);
				if (op == AD_EXT) return left.e->data[right.i] == other;
				if (other.op == AD_EXT) return *this == other.left.e->data[other.right.i];
				if (op != other.op) return false;
				switch (op)
				{
				case AD_MES: return true;
				case AD_CONST: return fabs(left.r - other.left.r) < 1e-9;
				case AD_PLUS:
				case AD_MINUS:
				case AD_MULT:
				case AD_DIV:
				case AD_POW:
					return left.i == other.left.i && right.i == other.right.i;
				case AD_COS:
				case AD_ABS:
				case AD_EXP:
				case AD_LOG:
				case AD_SIN:
Kirill Terekhov's avatar
Kirill Terekhov committed
396
				case AD_SQRT:
Kirill Terekhov's avatar
Kirill Terekhov committed
397
					return left.i == other.left.i;
Kirill Terekhov's avatar
Kirill Terekhov committed
398 399
				case AD_VAL:
					return left.i == other.right.i && (right.e == other.right.e || *right.e == *other.right.e);
Kirill Terekhov's avatar
Kirill Terekhov committed
400 401 402 403 404
				case AD_ALTR:
					if (left.e == other.left.e && right.e == other.right.e) return true;
					return *left.e == *other.left.e && *right.e == *other.right.e;
				case AD_COND:
					return left.i == other.left.i && (right.q == other.right.q || *right.q == *other.right.q);
Kirill Terekhov's avatar
Kirill Terekhov committed
405 406 407
				case AD_COND_TYPE:
				case AD_COND_MARK:
					return left.i == other.left.i;
Kirill Terekhov's avatar
Kirill Terekhov committed
408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439
				}
				if (op >= AD_FUNC) return true;
				if (op >= AD_TABLE) return left.i == other.left.i;
				if (op >= AD_STNCL) return (left.e == other.left.e) || (*left.e == *other.left.e);
				if (op >= AD_CTAG) return left.i == other.left.i;
				if (op >= AD_TAG) return left.i == other.left.i;
				assert(false); //cannot reach here
				return false;
			}
			bool deep_cmp(const expr_data & other, const expr & my_set, const expr & other_set ) const
			{
				assert(op != AD_NONE);
				
				if (op == AD_EXT) return left.e->data[right.i].deep_cmp(other,*left.e,other_set);
				if (other.op == AD_EXT) return deep_cmp(other.left.e->data[other.right.i],my_set, *other.left.e);
				if (op != other.op) return false;
				switch (op)
				{
				case AD_MES: return true;
				case AD_CONST: return fabs(left.r - other.left.r) < 1e-9;
				case AD_PLUS:
				case AD_MINUS:
				case AD_MULT:
				case AD_DIV:
				case AD_POW:
					return my_set.data[left.i].deep_cmp(other_set.data[other.left.i],my_set,other_set) 
						&& my_set.data[right.i].deep_cmp(other_set.data[other.right.i],my_set,other_set);
				case AD_COS:
				case AD_ABS:
				case AD_EXP:
				case AD_LOG:
				case AD_SIN:
Kirill Terekhov's avatar
Kirill Terekhov committed
440
				case AD_SQRT:
Kirill Terekhov's avatar
Kirill Terekhov committed
441
					return my_set.data[left.i].deep_cmp(other_set.data[other.left.i], my_set, other_set);
Kirill Terekhov's avatar
Kirill Terekhov committed
442 443
				case AD_VAL:
					return (right.e == other.right.e || *right.e == *other.right.e) && my_set.data[left.i].deep_cmp(other_set.data[other.left.i], my_set, other_set);
Kirill Terekhov's avatar
Kirill Terekhov committed
444 445 446 447 448
				case AD_ALTR:
					if (left.e == other.left.e && right.e == other.right.e) return true;
					return *left.e == *other.left.e && *right.e == *other.right.e;
				case AD_COND:
					return my_set.data[left.i].deep_cmp(other_set.data[other.left.i], my_set, other_set) && (right.q == other.right.q || *right.q == *other.right.q);
Kirill Terekhov's avatar
Kirill Terekhov committed
449 450 451
				case AD_COND_TYPE:
				case AD_COND_MARK:
					return left.i == other.left.i;
Kirill Terekhov's avatar
Kirill Terekhov committed
452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481
				}
				if (op >= AD_FUNC) return true;
				if (op >= AD_TABLE) return my_set.data[left.i].deep_cmp(other_set.data[other.left.i], my_set, other_set);
				if (op >= AD_STNCL) return (left.e == other.left.e) || (*left.e == *other.left.e);
				if (op >= AD_CTAG) return left.i == other.left.i;
				if (op >= AD_TAG) return left.i == other.left.i;
				assert(false); //cannot reach here
				return false;
			}
			bool cmp(const expr_data & other,replace_type & r) const
			{
				assert(op != AD_NONE);
				if (op == AD_EXT) return left.e->data[right.i] == other;
				if (other.op == AD_EXT) return *this == other.left.e->data[other.right.i];
				if (op != other.op) return false;
				switch (op)
				{
				case AD_MES: return true;
				case AD_CONST: return fabs(left.r - other.left.r) < 1e-9;
				case AD_PLUS:
				case AD_MINUS:
				case AD_MULT:
				case AD_DIV:
				case AD_POW:
					return left.i == r[other.left.i] && right.i == r[other.right.i];
				case AD_COS:
				case AD_ABS:
				case AD_EXP:
				case AD_LOG:
				case AD_SIN:
Kirill Terekhov's avatar
Kirill Terekhov committed
482
				case AD_SQRT:
Kirill Terekhov's avatar
Kirill Terekhov committed
483
					return left.i == r[other.left.i];
Kirill Terekhov's avatar
Kirill Terekhov committed
484 485
				case AD_VAL:
					return left.i == r[other.left.i] && (right.e == other.right.e || *right.e == *other.right.e);
Kirill Terekhov's avatar
Kirill Terekhov committed
486 487 488 489 490
				case AD_ALTR:
					if (left.e == other.left.e && right.e == other.right.e) return true;
					return *left.e == *other.left.e && *right.e == *other.right.e;
				case AD_COND:
					return left.i == r[other.left.i] && (right.q == other.right.q || *right.q == *other.right.q);
Kirill Terekhov's avatar
Kirill Terekhov committed
491 492 493
				case AD_COND_TYPE:
				case AD_COND_MARK:
					return left.i == other.left.i;
Kirill Terekhov's avatar
Kirill Terekhov committed
494 495 496 497 498 499 500 501 502 503 504 505 506 507 508
				}
				if (op >= AD_FUNC) return true;
				if (op >= AD_TABLE) return left.i == r[other.left.i];
				if (op >= AD_STNCL) return (left.e == other.left.e) || (*left.e == *other.left.e);
				if (op >= AD_CTAG) return left.i == other.left.i;
				if (op >= AD_TAG) return left.i == other.left.i;
				assert(false); //cannot reach here
				return false;
			}
			INMOST_DATA_ENUM_TYPE type() const
			{ 
				assert(op != AD_NONE);
				if (op == AD_EXT) return left.e->data[right.i].type();
				switch (op)
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
509 510
				case AD_COND_MARK:
				case AD_COND_TYPE:
Kirill Terekhov's avatar
Kirill Terekhov committed
511 512 513 514 515 516 517 518 519 520 521 522
				case AD_CONST:
				case AD_MES: return AD_TYPE_ENDPOINT;
				case AD_PLUS:
				case AD_MINUS:
				case AD_MULT:
				case AD_DIV:
				case AD_POW:
				case AD_ALTR: return AD_TYPE_BINARY;
				case AD_COS:
				case AD_ABS:
				case AD_EXP:
				case AD_LOG:
Kirill Terekhov's avatar
Kirill Terekhov committed
523
				case AD_SQRT:
Kirill Terekhov's avatar
Kirill Terekhov committed
524 525
				case AD_SIN: return AD_TYPE_UNARY;
				case AD_COND: return AD_TYPE_TERNARY;
Kirill Terekhov's avatar
Kirill Terekhov committed
526
				case AD_VAL: return AD_TYPE_VALUE;
Kirill Terekhov's avatar
Kirill Terekhov committed
527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547
				}
				if (op >= AD_FUNC) return AD_TYPE_ENDPOINT;
				if (op >= AD_TABLE) return AD_TYPE_UNARY;
				if (op >= AD_STNCL) return AD_TYPE_UNARY;
				if (op >= AD_CTAG) return AD_TYPE_ENDPOINT;
				if (op >= AD_TAG) return AD_TYPE_ENDPOINT;
				assert(false); //cannot reach here
				return AD_TYPE_INVALID;
			}
			bool is_func() { return op >= AD_FUNC; }
			bool is_table() { return op >= AD_TABLE && op < AD_FUNC; }
			bool is_stncl() { return op >= AD_STNCL && op < AD_TABLE; }
			bool is_ctag() { return op >= AD_CTAG && op < AD_STNCL; }
			bool is_tag() { return op >= AD_TAG && op < AD_CTAG; }
			friend class expr;
			friend class Automatizator;
		};
		friend class expr_data;
		//typedef std::vector<expr_data> data_type;
		typedef std::vector<expr_data> data_type;
		data_type data;
Kirill Terekhov's avatar
Kirill Terekhov committed
548 549 550 551 552 553 554 555 556 557
		std::vector<INMOST_DATA_REAL_TYPE> values;
		Automatizator::stencil_pairs current_stencil;
		__INLINE INMOST_DATA_ENUM_TYPE values_offset(INMOST_DATA_ENUM_TYPE element) {return element * static_cast<INMOST_DATA_ENUM_TYPE>(data.size());}
		__INLINE INMOST_DATA_ENUM_TYPE derivatives_offset(INMOST_DATA_ENUM_TYPE element) {return element * static_cast<INMOST_DATA_ENUM_TYPE>(data.size()) + static_cast<INMOST_DATA_ENUM_TYPE>(data.size()*current_stencil.size());}
		void resize_for_stencil()
		{
			values.resize(current_stencil.size()*data.size()*2);
			memset(&values[current_stencil.size()*data.size()],0,sizeof(INMOST_DATA_REAL_TYPE)*current_stencil.size()*data.size());
		}

Kirill Terekhov's avatar
Kirill Terekhov committed
558 559 560
		void relink_data()
		{
			for (data_type::iterator it = data.begin(); it != data.end(); ++it)
Kirill Terekhov's avatar
Kirill Terekhov committed
561
				if (it->op == AD_COND )
Kirill Terekhov's avatar
Kirill Terekhov committed
562
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
563 564 565 566 567 568 569 570 571 572
					for (data_type::iterator jt = it->right.q->left.e->data.begin(); jt != it->right.q->left.e->data.end(); ++jt)
					{
						if (jt->op == AD_EXT) jt->left.e = this;
					}
					it->right.q->left.e->relink_data();
					for (data_type::iterator jt = it->right.q->right.e->data.begin(); jt != it->right.q->right.e->data.end(); ++jt) 
					{
						if (jt->op == AD_EXT) jt->left.e = this;
					}
					it->right.q->left.e->relink_data();
Kirill Terekhov's avatar
Kirill Terekhov committed
573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625
				}
		}
		void link_data(data_type & other, replace_type & r)
		{
			INMOST_DATA_ENUM_TYPE i = 0, j;
			for (data_type::iterator it = other.begin(); it != other.end(); ++it)
			{
				j = 0;
				if (it->op == AD_EXT)
				{
					it->left.e = this;
					it->right.i = r[it->right.i];
				}
				else for (data_type::iterator jt = data.begin(); jt != data.end(); ++jt)
				{
					if (*it == *jt)
					{
						it->op = AD_EXT;
						it->left.e = this;
						it->right.i = j;
						break;
					}
					j++;
				}
				i++;
			}
		}
		void replace_operands(expr_data & e, replace_type & r)
		{
			if (e.op == AD_EXT) replace_operands(e.left.e->data[e.right.i],r);
			else if (e.op == AD_ALTR)
			{
				link_data(e.left.e->data,r);
				link_data(e.right.e->data,r);
			}
			else if (e.is_stncl()) {}
			else switch (e.type())
			{
			case AD_TYPE_UNARY:
				assert(r.is_present(e.left.i));
				e.left.i = r[e.left.i];
				break;
			case AD_TYPE_BINARY:
				assert(r.is_present(e.left.i));
				assert(r.is_present(e.right.i));
				e.left.i = r[e.left.i];
				e.right.i = r[e.right.i];
				break;
			case AD_TYPE_TERNARY:
				assert(r.is_present(e.left.i));
				e.left.i = r[e.left.i];
				replace_operands(*e.right.q, r);
				break;
Kirill Terekhov's avatar
Kirill Terekhov committed
626 627 628 629
			case AD_TYPE_VALUE:
				e.left.i = r[e.left.i];
				link_data(e.right.e->data,r);
				break;
Kirill Terekhov's avatar
Kirill Terekhov committed
630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649
			}
		}
		INMOST_DATA_ENUM_TYPE merge_data(const data_type & other)
		{
			replace_type from_to;
			INMOST_DATA_ENUM_TYPE i = 0, j;
			for (data_type::const_iterator it = other.begin(); it != other.end(); ++it)
			{
				j = 0;
				for (data_type::iterator jt = data.begin(); jt != data.end(); ++jt)
				{
					if (jt->cmp(*it,from_to))
					{
						from_to[i] = j;
						break;
					}
					j++;
				}
				if (!from_to.is_present(i))
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
650
					from_to[i] = static_cast<INMOST_DATA_ENUM_TYPE>(data.size());
Kirill Terekhov's avatar
Kirill Terekhov committed
651 652 653 654 655
					data.push_back(*it);
					replace_operands(data.back(),from_to);
				}
				i++;
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
656
			return from_to[static_cast<int>(other.size()) - 1];
Kirill Terekhov's avatar
Kirill Terekhov committed
657 658 659 660 661 662 663 664 665
		}
		
	public:
		expr() : data() { }
		expr(const expr & other) : data(other.data) { relink_data(); }
		expr(INMOST_DATA_ENUM_TYPE op, const expr * sub) :data() { data.push_back(expr_data(op,sub)); }
		expr(const expr * l, const expr * r) :data() { data.push_back(expr_data(l,r)); }
		expr(INMOST_DATA_REAL_TYPE val) : data() { data.push_back(expr_data(val)); }
		expr(INMOST_DATA_ENUM_TYPE op, INMOST_DATA_ENUM_TYPE comp) : data() { data.push_back(expr_data(op, comp, ENUMUNDEF)); }
Kirill Terekhov's avatar
Kirill Terekhov committed
666 667
		expr(INMOST_DATA_ENUM_TYPE op, const expr & operand) : data(operand.data) { relink_data();  data.push_back(expr_data(op, static_cast<INMOST_DATA_ENUM_TYPE>(data.size() - 1), ENUMUNDEF)); }
		expr(const expr & operand, const expr & multiplyer) : data(operand.data) { relink_data();  data.push_back(expr_data(AD_VAL, static_cast<INMOST_DATA_ENUM_TYPE>(data.size() - 1), multiplyer)); }
Kirill Terekhov's avatar
Kirill Terekhov committed
668 669 670
		expr(const expr & cond, const expr & if_true, const expr & if_false) :data(cond.data)
		{ 
			relink_data();  
Kirill Terekhov's avatar
Kirill Terekhov committed
671
			data.push_back(expr_data(static_cast<INMOST_DATA_ENUM_TYPE>(data.size() - 1), expr_data(&if_true, &if_false))); 
Kirill Terekhov's avatar
Kirill Terekhov committed
672
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
673 674 675
		expr(INMOST_DATA_ENUM_TYPE op, const expr & l, const expr & r) :data(l.data) 
		{
			relink_data();
Kirill Terekhov's avatar
Kirill Terekhov committed
676
			INMOST_DATA_ENUM_TYPE lp = static_cast<INMOST_DATA_ENUM_TYPE>(data.size() - 1);
Kirill Terekhov's avatar
Kirill Terekhov committed
677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706
			INMOST_DATA_ENUM_TYPE rp = merge_data(r.data); 
			data.push_back(expr_data(op, lp, rp));
		}
		~expr() {}
		expr & operator =(expr const & other) {data = other.data; return *this;}
		expr operator +() const { return expr(*this); }
		expr operator -() const { expr zero(0.0); return zero - *this; }
		expr operator +(const expr & other) const {return expr(AD_PLUS,*this,other);}
		expr operator -(const expr & other) const { return expr(AD_MINUS, *this, other); }
		expr operator *(const expr & other) const { return expr(AD_MULT, *this, other); }
		expr operator /(const expr & other) const { return expr(AD_DIV, *this, other); }
		expr operator +(INMOST_DATA_REAL_TYPE other) const { return expr(AD_PLUS, *this, expr(other)); }
		expr operator -(INMOST_DATA_REAL_TYPE other) const { return expr(AD_MINUS, *this, expr(other)); }
		expr operator *(INMOST_DATA_REAL_TYPE other) const { return expr(AD_MULT, *this, expr(other)); }
		expr operator /(INMOST_DATA_REAL_TYPE other) const { return expr(AD_DIV, *this, expr(other)); }
		bool operator ==(const expr & other) //this should account for possible reorder! here x*y is not the same as y*x
		{
			if (data.size() != other.data.size()) return false;
			data_type::iterator it = data.begin();
			data_type::const_iterator jt = other.data.begin();
			while (it != data.end() && jt != other.data.end())
			{
				if (!it->deep_cmp(*jt, *this, other)) return false;
				++it; ++jt;
			}
			return true;
		}
		
		friend class Automatizator;
	};
707
}
Kirill Terekhov's avatar
Kirill Terekhov committed
708

709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732
__INLINE INMOST::expr operator +(INMOST_DATA_REAL_TYPE left, const INMOST::expr & right) { return INMOST::expr(left) + right; }
__INLINE INMOST::expr operator -(INMOST_DATA_REAL_TYPE left, const INMOST::expr & right) { return INMOST::expr(left) - right; }
__INLINE INMOST::expr operator *(INMOST_DATA_REAL_TYPE left, const INMOST::expr & right) { return INMOST::expr(left) * right; }
__INLINE INMOST::expr operator /(INMOST_DATA_REAL_TYPE left, const INMOST::expr & right) { return INMOST::expr(left) / right; }
__INLINE INMOST::expr ad_pow(const INMOST::expr & v, const INMOST::expr & n) { return INMOST::expr(AD_POW, v, n); }
//__INLINE INMOST::expr ad_pow(const INMOST::expr & v, INMOST_DATA_REAL_TYPE n) { return INMOST::expr(AD_POW, v, INMOST::expr(n)); }
//__INLINE INMOST::expr ad_pow(const INMOST::expr & v, INMOST_DATA_ENUM_TYPE n) { return INMOST::expr(AD_POW, v, INMOST::expr(static_cast<INMOST_DATA_REAL_TYPE>(n))); }
__INLINE INMOST::expr ad_abs(const INMOST::expr & v) { return INMOST::expr(AD_ABS, v); }
__INLINE INMOST::expr ad_exp(const INMOST::expr & v) { return INMOST::expr(AD_EXP, v); }
__INLINE INMOST::expr ad_log(const INMOST::expr & v) { return INMOST::expr(AD_LOG, v); }
__INLINE INMOST::expr ad_sin(const INMOST::expr & v) { return INMOST::expr(AD_SIN, v); }
__INLINE INMOST::expr ad_cos(const INMOST::expr & v) { return INMOST::expr(AD_COS, v); }
__INLINE INMOST::expr ad_sqrt(const INMOST::expr & v) {return INMOST::expr(AD_SQRT, v);}
__INLINE INMOST::expr ad_val(const INMOST::expr & v, const INMOST::expr & multiplyer = INMOST::expr(0.0)) {return INMOST::expr(v, multiplyer);}
__INLINE INMOST::expr measure() { return INMOST::expr(AD_MES,ENUMUNDEF); }
__INLINE INMOST::expr condition(const INMOST::expr & cond, const INMOST::expr & if_gt_zero, const INMOST::expr & if_le_zero) { return INMOST::expr(cond, if_gt_zero, if_le_zero); }
__INLINE INMOST::expr condition_etype(INMOST::ElementType etypes, const INMOST::expr & if_true, const INMOST::expr & if_false) { return INMOST::expr(INMOST::expr(AD_COND_TYPE, etypes), if_true, if_false); }
__INLINE INMOST::expr condition_marker(INMOST::MarkerType marker, const INMOST::expr & if_true, const INMOST::expr & if_false) { return INMOST::expr(INMOST::expr(AD_COND_MARK, marker), if_true, if_false); }
__INLINE INMOST::expr stencil(INMOST_DATA_ENUM_TYPE stncl, const INMOST::expr & v) { assert(stncl >= AD_STNCL && stncl < AD_TABLE); return INMOST::expr(stncl, &v); }
__INLINE INMOST::expr tabval(INMOST_DATA_ENUM_TYPE tabl, const INMOST::expr & v) { assert(tabl >= AD_TABLE && tabl < AD_FUNC); return INMOST::expr(tabl, v); }
__INLINE INMOST::expr tagval(INMOST_DATA_ENUM_TYPE reg_tag, INMOST_DATA_ENUM_TYPE comp = 0) { assert(reg_tag >= AD_TAG && reg_tag < AD_STNCL); return INMOST::expr(reg_tag, comp); }
__INLINE INMOST::expr funcval(INMOST_DATA_ENUM_TYPE reg_func) { assert(reg_func >= AD_FUNC); return INMOST::expr(reg_func, ENUMUNDEF); }
namespace INMOST
{
Kirill Terekhov's avatar
Kirill Terekhov committed
733
#else
734
}
Kirill Terekhov's avatar
Kirill Terekhov committed
735

736 737 738 739 740 741 742
__INLINE INMOST::expr operator +(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right);
__INLINE INMOST::expr operator -(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right);
__INLINE INMOST::expr operator *(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right);
__INLINE INMOST::expr operator /(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right);
  
namespace INMOST
{
Kirill Terekhov's avatar
Kirill Terekhov committed
743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759
	class expr
	{
		INMOST_DATA_ENUM_TYPE op;
		INMOST_DATA_REAL_TYPE coef;
		expr * left;
		expr * right;
	public:
		expr() : op(AD_NONE), coef(1), left(NULL), right(NULL) { }
		expr(const expr & other) : op(other.op), coef(other.coef)
		{
			if ((other.op >= AD_TAG && other.op < AD_STNCL) || other.op == AD_COND_TYPE || other.op == AD_COND_MARK) left = other.left; //copy component number
			else if (other.left != NULL) left = new expr(*other.left); else left = NULL;
			if (other.right != NULL) right = new expr(*other.right); else right = NULL;
		}
		expr(INMOST_DATA_REAL_TYPE val) : op(AD_CONST), coef(val), left(NULL), right(NULL) {}
		//expr(INMOST_DATA_ENUM_TYPE val) : op(AD_CONST), coef(val), left(NULL), right(NULL) {}
		expr(INMOST_DATA_ENUM_TYPE new_op, expr * l, expr * r) : op(new_op), coef(1), left(l), right(r) {}
Kirill Terekhov's avatar
Kirill Terekhov committed
760 761 762 763
		expr(INMOST_DATA_ENUM_TYPE tag_op, INMOST_DATA_ENUM_TYPE comp) :op(tag_op), coef(1), right(NULL) 
		{
			*(INMOST_DATA_ENUM_TYPE *)(&left) = comp;
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797
		~expr()
		{
			if (op < AD_COS)
			{
				delete left;
				delete right;
			}
			else if (op < AD_CONST)
				delete left;
			else if (op >= AD_STNCL && op < AD_TABLE)
				delete left;
		}
		expr & operator =(expr const & other)
		{
			op = other.op; coef = other.coef;
			if (other.op >= AD_TAG && other.op < AD_STNCL)
			{
				left = other.left; //copy component number
				right = other.right;
			}
			else if (other.left != NULL) left = new expr(*other.left); else left = NULL;
			if (other.right != NULL) right = new expr(*other.right); else right = NULL;
			return *this;
		}
		expr operator +() { return expr(*this); }
		expr operator -() { expr var(*this); var.coef *= -1.0; return var; }
		expr operator +(const expr & other) const { return expr(AD_PLUS, new expr(*this), new expr(other)); }
		expr operator -(const expr & other) const { return expr(AD_MINUS, new expr(*this), new expr(other)); }
		expr operator *(const expr & other) const { return expr(AD_MULT, new expr(*this), new expr(other)); }
		expr operator /(const expr & other) const { return expr(AD_DIV, new expr(*this), new expr(other)); }
		expr operator +(const INMOST_DATA_REAL_TYPE & other) const { return expr(AD_PLUS, new expr(*this), new expr(other)); }
		expr operator -(const INMOST_DATA_REAL_TYPE & other) const { return expr(AD_MINUS, new expr(*this), new expr(other)); }
		expr operator *(const INMOST_DATA_REAL_TYPE & other) const { expr var(*this); var.coef *= other; return var; }
		expr operator /(const INMOST_DATA_REAL_TYPE & other) const { expr var(*this); var.coef /= other; return var; }
798 799 800 801
		__INLINE friend expr (::operator +)(const INMOST_DATA_REAL_TYPE & left, const expr & right);
		__INLINE friend expr (::operator -)(const INMOST_DATA_REAL_TYPE & left, const expr & right);
		__INLINE friend expr (::operator *)(const INMOST_DATA_REAL_TYPE & left, const expr & right);
		__INLINE friend expr (::operator /)(const INMOST_DATA_REAL_TYPE & left, const expr & right);
Kirill Terekhov's avatar
Kirill Terekhov committed
802 803 804 805 806
		bool operator ==(const expr & other) {return this == &other || (op == other.op && left == other.left && right == other.right);}
		bool is_endpoint() { return op >= AD_TAG && op < AD_STNCL; }
		friend class Automatizator;
	};
	
807
}
Kirill Terekhov's avatar
Kirill Terekhov committed
808

809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830
__INLINE INMOST::expr operator +(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right) { return INMOST::expr(AD_PLUS, new INMOST::expr(left), new INMOST::expr(right)); }
__INLINE INMOST::expr operator -(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right) { return INMOST::expr(AD_MINUS, new INMOST::expr(left), new INMOST::expr(right)); }
__INLINE INMOST::expr operator *(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right) { INMOST::expr var(right); var.coef *= left; return var; }
__INLINE INMOST::expr operator /(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right) { INMOST::expr var; var.op = AD_INV; var.coef = left; var.left = new INMOST::expr(right); var.right = NULL; return var; }
__INLINE INMOST::expr ad_pow(const INMOST::expr & v, const INMOST::expr n) { return INMOST::expr(AD_POW, new INMOST::expr(v), new INMOST::expr(n)); }
__INLINE INMOST::expr ad_abs(const INMOST::expr & v) { return INMOST::expr(AD_ABS, new INMOST::expr(v), NULL); }
__INLINE INMOST::expr ad_exp(const INMOST::expr & v) { return INMOST::expr(AD_EXP, new INMOST::expr(v), NULL); }
__INLINE INMOST::expr ad_log(const INMOST::expr & v) { return INMOST::expr(AD_LOG, new INMOST::expr(v), NULL); }
__INLINE INMOST::expr ad_sin(const INMOST::expr & v) { return INMOST::expr(AD_SIN, new INMOST::expr(v), NULL); }
__INLINE INMOST::expr ad_cos(const INMOST::expr & v) { return INMOST::expr(AD_COS, new INMOST::expr(v), NULL); }
__INLINE INMOST::expr ad_sqrt(const INMOST::expr & v) {return INMOST::expr(AD_SQRT, new INMOST::expr(v), NULL);}
__INLINE INMOST::expr ad_val(const INMOST::expr & v, const INMOST::expr & multiplyer = INMOST::expr(0.0)) {return INMOST::expr(AD_VAL,new INMOST::expr(v), new INMOST::expr(multiplyer));}
__INLINE INMOST::expr measure() { return INMOST::expr(AD_MES, NULL, NULL); }
__INLINE INMOST::expr condition_etype(INMOST::ElementType etype, const INMOST::expr & if_true, const INMOST::expr & if_false) { return INMOST::expr(AD_COND, new INMOST::expr(AD_COND_TYPE,etype), new INMOST::expr(AD_ALTR, new INMOST::expr(if_true), new INMOST::expr(if_false))); }
__INLINE INMOST::expr condition_marker(INMOST::MarkerType marker, const INMOST::expr & if_true, const INMOST::expr & if_false) { return INMOST::expr(AD_COND, new INMOST::expr(AD_COND_MARK,marker), new INMOST::expr(AD_ALTR, new INMOST::expr(if_true), new INMOST::expr(if_false))); }
__INLINE INMOST::expr condition(const INMOST::expr & cond, const INMOST::expr & if_gt_zero, const INMOST::expr & if_le_zero) { return INMOST::expr(AD_COND, new INMOST::expr(cond), new INMOST::expr(AD_ALTR, new INMOST::expr(if_gt_zero), new INMOST::expr(if_le_zero))); }
__INLINE INMOST::expr stencil(INMOST_DATA_ENUM_TYPE stncl, const INMOST::expr & v) { assert(stncl >= AD_STNCL && stncl < AD_TABLE); return INMOST::expr(stncl, new INMOST::expr(v), NULL); }
__INLINE INMOST::expr tabval(INMOST_DATA_ENUM_TYPE tabl, const INMOST::expr & v) { assert(tabl >= AD_TABLE && tabl < AD_FUNC); return INMOST::expr(tabl, new INMOST::expr(v), NULL); }
__INLINE INMOST::expr tagval(INMOST_DATA_ENUM_TYPE reg_tag, INMOST_DATA_ENUM_TYPE comp = 0) { assert(reg_tag >= AD_TAG && reg_tag < AD_STNCL); return INMOST::expr(reg_tag, comp); }
__INLINE INMOST::expr funcval(INMOST_DATA_ENUM_TYPE reg_func) { assert(reg_func >= AD_FUNC); return INMOST::expr(reg_func, NULL, NULL); }
namespace INMOST
{
Kirill Terekhov's avatar
Kirill Terekhov committed
831
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
832 833 834 835
}

#endif
#endif