inmost_autodiff.h 31.4 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1
2
3
4
5
6
7
#pragma once
#ifndef INMOST_AUTODIFF_H_INCLUDED
#define INMOST_AUTODIFF_H_INCLUDED
#include "inmost_common.h"
#include "inmost_mesh.h"
#include <sstream> //for debug

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

#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
33
#define AD_INV   5 //inversion 1/x, depricated
Kirill Terekhov's avatar
Kirill Terekhov committed
34
#define AD_POW   6 //power operation x^y
Kirill Terekhov's avatar
Kirill Terekhov committed
35
36
#define AD_SQRT  7 //square root operation

Kirill Terekhov's avatar
Kirill Terekhov committed
37
38
39

#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
40
41
#define AD_VAL   19 // calculate only value of current expression

Kirill Terekhov's avatar
Kirill Terekhov committed
42
43
44
45
46
47
48
49
50
51
52
//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
53
54
55
56
57
58
59
#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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#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
76
#define AD_TYPE_VALUE     4
Kirill Terekhov's avatar
Kirill Terekhov committed
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110

namespace INMOST
{
	class Automatizator; //forward declaration
	typedef dynarray<INMOST_DATA_REAL_TYPE, 2048> precomp_values_t;
	class expr;
	
	

	class Automatizator
	{
	public:
		typedef INMOST_DATA_REAL_TYPE(*func_callback)(Storage * current_element, void * user_data);
		typedef std::pair<Storage *, INMOST_DATA_REAL_TYPE> stencil_pair;
		typedef dynarray<stencil_pair, 64> stencil_pairs;
		typedef void(*stencil_callback)(Storage * current_element, stencil_pairs & out_stencil, void * user_data);
	private:
		typedef struct{ Tag t; MIDType domain_mask; } tagdomain;
		typedef small_hash<INMOST_DATA_ENUM_TYPE, tagdomain, 128> const_tag_type;
		typedef struct{ tagdomain d; Tag indices; } tagpair;
		typedef small_hash<INMOST_DATA_ENUM_TYPE, tagpair, 128> tagpairs_type;
		typedef std::vector<tagpair> index_enum;
		typedef struct { Tag elements, coefs; } stencil_tag;
		typedef struct { std::string name; INMOST_DATA_ENUM_TYPE kind; MIDType domainmask; void * link; } stencil_kind_domain;
		typedef small_hash<INMOST_DATA_ENUM_TYPE, stencil_kind_domain, 128> stencil_type;
		typedef struct func_name_callback_t { std::string name; func_callback func; } func_name_callback;
		typedef small_hash<INMOST_DATA_ENUM_TYPE, func_name_callback, 128> func_type;
	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
111
112
113
114
			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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
		} table;
		typedef table * table_ptr;
	private:
		typedef dynarray<INMOST_DATA_REAL_TYPE, 128> values_container;
		typedef small_hash<INMOST_DATA_ENUM_TYPE, 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;
		INMOST_DATA_ENUM_TYPE first_num;
		INMOST_DATA_ENUM_TYPE last_num;
		Mesh * m;
Kirill Terekhov's avatar
Kirill Terekhov committed
129
130
131
132
#if defined(NEW_VERSION)
		void DerivativeFill(expr & var, INMOST_DATA_ENUM_TYPE element, INMOST_DATA_ENUM_TYPE parent, Solver::Row & entries, INMOST_DATA_REAL_TYPE multval, void * user_data);
		INMOST_DATA_REAL_TYPE EvaluateSub(expr & var,INMOST_DATA_ENUM_TYPE element,INMOST_DATA_ENUM_TYPE parent, void * user_data);
#else
Kirill Terekhov's avatar
Kirill Terekhov committed
133
134
		INMOST_DATA_REAL_TYPE DerivativePrecompute(const expr & var, Storage * e, precomp_values_t & values, void * user_data);
		void DerivativeFill(const expr & var, Storage * e, Solver::Row & entries, precomp_values_t & values, INMOST_DATA_REAL_TYPE multval, void * user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
135
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
	public:
		Automatizator(Mesh * m);
		~Automatizator();
		__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, MIDType domain_mask = 0);
		INMOST_DATA_ENUM_TYPE RegisterStencil(std::string name, stencil_callback func, MIDType 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, MIDType domain_mask = 0);
		INMOST_DATA_ENUM_TYPE RegisterStaticTag(Tag t, MIDType domain_mask = 0);
		void EnumerateDynamicTags();
		__INLINE Tag                 GetDynamicValueTag(INMOST_DATA_ENUM_TYPE ind) { return reg_tags[ind].d.t; }
		__INLINE Tag                 GetDynamicIndexTag(INMOST_DATA_ENUM_TYPE ind) { return reg_tags[ind].indices; }
		__INLINE MIDType             GetDynamicMask(INMOST_DATA_ENUM_TYPE ind) { return reg_tags[ind].d.domain_mask; }
		__INLINE Tag                 GetStaticValueTag(INMOST_DATA_ENUM_TYPE ind) { return reg_ctags[ind].t; }
		__INLINE MIDType             GetStaticMask(INMOST_DATA_ENUM_TYPE ind) { return reg_ctags[ind].domain_mask; }
		__INLINE INMOST_DATA_REAL_TYPE GetDynamicValue(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(Storage * e, INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_ENUM_TYPE comp = 0) { return e->IntegerArray(GetDynamicIndexTag(ind))[comp]; }
		__INLINE bool                isDynamicValid(Storage * e, INMOST_DATA_ENUM_TYPE ind) { MIDType mask = GetDynamicMask(ind); return mask == 0 || e->GetMarker(mask); }
		__INLINE INMOST_DATA_REAL_TYPE GetStaticValue(Storage * e, INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_ENUM_TYPE comp = 0) { return e->RealArray(GetStaticValueTag(ind))[comp]; }
		__INLINE bool                isStaticValid(Storage * e, INMOST_DATA_ENUM_TYPE ind) { MIDType mask = GetStaticMask(ind); return mask == 0 || e->GetMarker(mask); }
Kirill Terekhov's avatar
Kirill Terekhov committed
158
159
160
161
#if defined(NEW_VERSION)
		INMOST_DATA_REAL_TYPE Evaluate(expr & var, Storage * e, void * user_data);
		INMOST_DATA_REAL_TYPE Derivative(expr & var, Storage * e, Solver::Row & out, Storage::real multiply, void * user_data);
#else
Kirill Terekhov's avatar
Kirill Terekhov committed
162
		INMOST_DATA_REAL_TYPE Evaluate(const expr & var, Storage * e, void * user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
163
164
165
166
167
168
169
170
171
172
		INMOST_DATA_REAL_TYPE Derivative(const expr & var, Storage * e, Solver::Row & out, Storage::real multiply,  void * user_data);
#endif
		__INLINE INMOST_DATA_REAL_TYPE                                   GetIndex(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(Storage *e, INMOST_DATA_ENUM_TYPE tagind) { return static_cast<INMOST_DATA_ENUM_TYPE>(e->IntegerArray(GetDynamicIndexTag(tagind)).size()); }
		__INLINE Mesh *                                                  GetMesh() { return m; }
		__INLINE INMOST_DATA_REAL_TYPE *                                 GetTableArguments(INMOST_DATA_ENUM_TYPE tableind) {return reg_tables[tableind]->args;}
		__INLINE INMOST_DATA_REAL_TYPE *                                 GetTableValues(INMOST_DATA_ENUM_TYPE tableind) {return reg_tables[tableind]->vals;}
		__INLINE INMOST_DATA_ENUM_TYPE                                   GetTableSize(INMOST_DATA_ENUM_TYPE tableind) {return reg_tables[tableind]->size;}
		__INLINE INMOST_DATA_REAL_TYPE                                   GetTableValue(INMOST_DATA_ENUM_TYPE tableind, INMOST_DATA_REAL_TYPE arg) { return reg_tables[tableind]->get_value(arg); }
		__INLINE INMOST_DATA_REAL_TYPE                                   GetTableDerivative(INMOST_DATA_ENUM_TYPE tableind, INMOST_DATA_REAL_TYPE arg) { return reg_tables[tableind]->get_derivative(arg); }
Kirill Terekhov's avatar
Kirill Terekhov committed
173
		__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]->get_both(arg); }
Kirill Terekhov's avatar
Kirill Terekhov committed
174
		INMOST_DATA_ENUM_TYPE                                            GetStencil(INMOST_DATA_ENUM_TYPE stnclind, Storage * elem, void * user_data, stencil_pairs & ret);
Kirill Terekhov's avatar
Kirill Terekhov committed
175
176
177
		__INLINE INMOST_DATA_REAL_TYPE GetFunction(INMOST_DATA_ENUM_TYPE funcid, Storage * elem, void * user_data) { return reg_funcs[funcid].func(elem, user_data); }
	};

Kirill Terekhov's avatar
Kirill Terekhov committed
178
#if defined(NEW_VERSION)
Kirill Terekhov's avatar
Kirill Terekhov committed
179
180
181
182
183
184
185
186
	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
187
			INMOST_DATA_ENUM_TYPE op, priority;
Kirill Terekhov's avatar
Kirill Terekhov committed
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
			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
235
236
			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
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
			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
280
				case AD_SQRT:
Kirill Terekhov's avatar
Kirill Terekhov committed
281
					return left.i == other.left.i;
Kirill Terekhov's avatar
Kirill Terekhov committed
282
283
				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
284
285
286
287
288
				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
289
290
291
				case AD_COND_TYPE:
				case AD_COND_MARK:
					return left.i == other.left.i;
Kirill Terekhov's avatar
Kirill Terekhov committed
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
				}
				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
324
				case AD_SQRT:
Kirill Terekhov's avatar
Kirill Terekhov committed
325
					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
326
327
				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
328
329
330
331
332
				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
333
334
335
				case AD_COND_TYPE:
				case AD_COND_MARK:
					return left.i == other.left.i;
Kirill Terekhov's avatar
Kirill Terekhov committed
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
				}
				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
366
				case AD_SQRT:
Kirill Terekhov's avatar
Kirill Terekhov committed
367
					return left.i == r[other.left.i];
Kirill Terekhov's avatar
Kirill Terekhov committed
368
369
				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
370
371
372
373
374
				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
375
376
377
				case AD_COND_TYPE:
				case AD_COND_MARK:
					return left.i == other.left.i;
Kirill Terekhov's avatar
Kirill Terekhov committed
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
				}
				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
393
394
				case AD_COND_MARK:
				case AD_COND_TYPE:
Kirill Terekhov's avatar
Kirill Terekhov committed
395
396
397
398
399
400
401
402
403
404
405
406
				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
407
				case AD_SQRT:
Kirill Terekhov's avatar
Kirill Terekhov committed
408
409
				case AD_SIN: return AD_TYPE_UNARY;
				case AD_COND: return AD_TYPE_TERNARY;
Kirill Terekhov's avatar
Kirill Terekhov committed
410
				case AD_VAL: return AD_TYPE_VALUE;
Kirill Terekhov's avatar
Kirill Terekhov committed
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
				}
				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
432
433
434
435
436
437
438
439
440
441
		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
442
443
444
		void relink_data()
		{
			for (data_type::iterator it = data.begin(); it != data.end(); ++it)
Kirill Terekhov's avatar
Kirill Terekhov committed
445
				if (it->op == AD_COND )
Kirill Terekhov's avatar
Kirill Terekhov committed
446
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
447
448
449
450
451
452
453
454
455
456
					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
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
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
				}
		}
		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
510
511
512
513
			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
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
			}
		}
		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
534
					from_to[i] = static_cast<INMOST_DATA_ENUM_TYPE>(data.size());
Kirill Terekhov's avatar
Kirill Terekhov committed
535
536
537
538
539
					data.push_back(*it);
					replace_operands(data.back(),from_to);
				}
				i++;
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
540
			return from_to[static_cast<int>(other.size()) - 1];
Kirill Terekhov's avatar
Kirill Terekhov committed
541
542
543
544
545
546
547
548
549
550
		}
		
	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)); }
		expr(INMOST_DATA_ENUM_TYPE op, const expr & operand) : data(operand.data) { relink_data();  data.push_back(expr_data(op, data.size() - 1, ENUMUNDEF)); }
Kirill Terekhov's avatar
Kirill Terekhov committed
551
552
553
554
555
556
		expr(const expr & operand, const expr & multiplyer) : data(operand.data) { relink_data();  data.push_back(expr_data(AD_VAL, data.size() - 1, multiplyer)); }
		expr(const expr & cond, const expr & if_true, const expr & if_false) :data(cond.data)
		{ 
			relink_data();  
			data.push_back(expr_data(data.size() - 1, expr_data(&if_true, &if_false))); 
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
		expr(INMOST_DATA_ENUM_TYPE op, const expr & l, const expr & r) :data(l.data) 
		{
			relink_data();
			INMOST_DATA_ENUM_TYPE lp = data.size() - 1;
			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;
	};
	

	__INLINE expr operator +(INMOST_DATA_REAL_TYPE left, const expr & right) { return expr(left) + right; }
	__INLINE expr operator -(INMOST_DATA_REAL_TYPE left, const expr & right) { return expr(left) - right; }
	__INLINE expr operator *(INMOST_DATA_REAL_TYPE left, const expr & right) { return expr(left) * right; }
	__INLINE expr operator /(INMOST_DATA_REAL_TYPE left, const expr & right) { return expr(left) / right; }
Kirill Terekhov's avatar
Kirill Terekhov committed
597
598
599
600
601
602
603
604
605
606
	__INLINE expr ad_pow(const expr & v, const expr & n) { return expr(AD_POW, v, n); }
	//__INLINE expr ad_pow(const expr & v, INMOST_DATA_REAL_TYPE n) { return expr(AD_POW, v, expr(n)); }
	//__INLINE expr ad_pow(const expr & v, INMOST_DATA_ENUM_TYPE n) { return expr(AD_POW, v, expr(static_cast<INMOST_DATA_REAL_TYPE>(n))); }
	__INLINE expr ad_abs(const expr & v) { return expr(AD_ABS, v); }
	__INLINE expr ad_exp(const expr & v) { return expr(AD_EXP, v); }
	__INLINE expr ad_log(const expr & v) { return expr(AD_LOG, v); }
	__INLINE expr ad_sin(const expr & v) { return expr(AD_SIN, v); }
	__INLINE expr ad_cos(const expr & v) { return expr(AD_COS, v); }
	__INLINE expr ad_sqrt(const expr & v) {return expr(AD_SQRT, v);}
	__INLINE expr ad_val(const expr & v, const expr & multiplyer = expr(0.0)) {return expr(v, multiplyer);}
Kirill Terekhov's avatar
Kirill Terekhov committed
607
	__INLINE expr measure() { return expr(AD_MES,ENUMUNDEF); }
Kirill Terekhov's avatar
Kirill Terekhov committed
608
609
610
	__INLINE expr condition(const expr & cond, const expr & if_gt, const expr & if_le) { return expr(cond, if_gt, if_le); }
	__INLINE expr condition_etype(ElementType etypes, const expr & if_true, const expr & if_false) { return expr(expr(AD_COND_TYPE, etypes), if_true, if_false); }
	__INLINE expr condition_marker(MIDType marker, const expr & if_true, const expr & if_false) { return expr(expr(AD_COND_MARK, marker), if_true, if_false); }
Kirill Terekhov's avatar
Kirill Terekhov committed
611
612
613
614
	__INLINE expr stencil(INMOST_DATA_ENUM_TYPE stncl, const expr & v) { assert(stncl >= AD_STNCL && stncl < AD_TABLE); return expr(stncl, &v); }
	__INLINE expr tabval(INMOST_DATA_ENUM_TYPE tabl, const expr & v) { assert(tabl >= AD_TABLE && tabl < AD_FUNC); return expr(tabl, v); }
	__INLINE expr tagval(INMOST_DATA_ENUM_TYPE reg_tag, INMOST_DATA_ENUM_TYPE comp = 0) { assert(reg_tag >= AD_TAG && reg_tag < AD_STNCL); return expr(reg_tag, comp); }
	__INLINE expr funcval(INMOST_DATA_ENUM_TYPE reg_func) { assert(reg_func >= AD_FUNC); return expr(reg_func, ENUMUNDEF); }
Kirill Terekhov's avatar
Kirill Terekhov committed
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
#else

	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) {}
		expr(INMOST_DATA_ENUM_TYPE tag_op, INMOST_DATA_ENUM_TYPE comp) :op(tag_op), coef(1), left(reinterpret_cast<expr *>(comp)), right(NULL) {}
		~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; }
		__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);
		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;
	};
	

	__INLINE expr operator +(const INMOST_DATA_REAL_TYPE & left, const expr & right) { return expr(AD_PLUS, new expr(left), new expr(right)); }
	__INLINE expr operator -(const INMOST_DATA_REAL_TYPE & left, const expr & right) { return expr(AD_MINUS, new expr(left), new expr(right)); }
	__INLINE expr operator *(const INMOST_DATA_REAL_TYPE & left, const expr & right) { expr var(right); var.coef *= left; return var; }
	__INLINE expr operator /(const INMOST_DATA_REAL_TYPE & left, const expr & right) { expr var; var.op = AD_INV; var.coef = left; var.left = new expr(right); var.right = NULL; return var; }
	__INLINE expr ad_pow(const expr & v, const expr n) { return expr(AD_POW, new expr(v), new expr(n)); }
	__INLINE expr ad_abs(const expr & v) { return expr(AD_ABS, new expr(v), NULL); }
	__INLINE expr ad_exp(const expr & v) { return expr(AD_EXP, new expr(v), NULL); }
	__INLINE expr ad_log(const expr & v) { return expr(AD_LOG, new expr(v), NULL); }
	__INLINE expr ad_sin(const expr & v) { return expr(AD_SIN, new expr(v), NULL); }
	__INLINE expr ad_cos(const expr & v) { return expr(AD_COS, new expr(v), NULL); }
	__INLINE expr ad_sqrt(const expr & v) {return expr(AD_SQRT, new expr(v), NULL);}
	__INLINE expr ad_val(const expr & v, const expr & multiplyer = expr(0.0)) {return expr(AD_VAL,new expr(v), new expr(multiplyer));}
	__INLINE expr measure() { return expr(AD_MES, NULL, NULL); }
	__INLINE expr condition_etype(ElementType etype, const expr & if_true, const expr & if_false) { return expr(AD_COND, new expr(AD_COND_TYPE,etype), new expr(AD_ALTR, new expr(if_true), new expr(if_false))); }
	__INLINE expr condition_marker(MIDType marker, const expr & if_true, const expr & if_false) { return expr(AD_COND, new expr(AD_COND_MARK,marker), new expr(AD_ALTR, new expr(if_true), new expr(if_false))); }
	__INLINE expr condition(const expr & cond, const expr & if_true, const expr & if_false) { return expr(AD_COND, new expr(cond), new expr(AD_ALTR, new expr(if_true), new expr(if_false))); }
	__INLINE expr stencil(INMOST_DATA_ENUM_TYPE stncl, const expr & v) { assert(stncl >= AD_STNCL && stncl < AD_TABLE); return expr(stncl, new expr(v), NULL); }
	__INLINE expr tabval(INMOST_DATA_ENUM_TYPE tabl, const expr & v) { assert(tabl >= AD_TABLE && tabl < AD_FUNC); return expr(tabl, new expr(v), NULL); }
	__INLINE expr tagval(INMOST_DATA_ENUM_TYPE reg_tag, INMOST_DATA_ENUM_TYPE comp = 0) { assert(reg_tag >= AD_TAG && reg_tag < AD_STNCL); return expr(reg_tag, comp); }
	__INLINE expr funcval(INMOST_DATA_ENUM_TYPE reg_func) { assert(reg_func >= AD_FUNC); return expr(reg_func, NULL, NULL); }
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
700
701
702
703
}

#endif
#endif