Commit d5434456 authored by Kirill Terekhov's avatar Kirill Terekhov

Cleanup deprecated code from autodiff

parent de11b0be
......@@ -13,633 +13,43 @@
namespace INMOST
{
Automatizator * Automatizator::CurrentAutomatizator = NULL;
bool print_ad_ctor = false;
bool GetAutodiffPrint() {return print_ad_ctor;}
void SetAutodiffPrint(bool set) {print_ad_ctor = set;}
bool CheckCurrentAutomatizator() {return Automatizator::HaveCurrent();}
Automatizator * Automatizator::CurrentAutomatizator = NULL;
bool print_ad_ctor = false;
bool GetAutodiffPrint() {return print_ad_ctor;}
void SetAutodiffPrint(bool set) {print_ad_ctor = set;}
bool CheckCurrentAutomatizator() {return Automatizator::HaveCurrent();}
void FromBasicExpression(Sparse::Row & entries, const basic_expression & expr)
{
Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
expr.GetJacobian(1.0,merger);
merger.RetriveRow(entries);
merger.Clear();
}
void AddBasicExpression(Sparse::Row & entries, INMOST_DATA_REAL_TYPE multme, INMOST_DATA_REAL_TYPE multit, const basic_expression & expr)
{
Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
merger.PushRow(multme,entries);
expr.GetJacobian(multit,merger);
merger.RetriveRow(entries);
merger.Clear();
}
void FromGetJacobian(const basic_expression & expr, INMOST_DATA_REAL_TYPE mult, Sparse::Row & r)
{
Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
expr.GetJacobian(mult,merger);
merger.AddRow(1.0,r);
merger.RetriveRow(r);
merger.Clear();
}
#if defined(NEW_VERSION)
//! returns offset from the end of precomputed z
void Automatizator::DerivativeFill(expr & var, INMOST_DATA_ENUM_TYPE element, INMOST_DATA_ENUM_TYPE parent, Sparse::RowMerger & entries, INMOST_DATA_REAL_TYPE multval, void * user_data)
void FromBasicExpression(Sparse::Row & entries, const basic_expression & expr)
{
INMOST_DATA_ENUM_TYPE voffset = var.values_offset(element), doffset = var.derivatives_offset(element);
INMOST_DATA_ENUM_TYPE k = static_cast<INMOST_DATA_ENUM_TYPE>(var.data.size()-1);
Storage e = Storage(m,var.current_stencil[element].first);
INMOST_DATA_REAL_TYPE lval, rval, ret;
var.values[doffset+k] = multval;
expr::expr_data * arr = &var.data[0], *it;
//for (expr::data_type::reverse_iterator it = var.data.rbegin(); it != var.data.rend(); ++it)
do
{
it = arr+k;
switch (it->op)
{
case AD_EXT:
assert(parent != ENUMUNDEF);
it->left.e->values[it->left.e->derivatives_offset(parent)+it->right.i] += var.values[doffset+k];
break;
case AD_COND:
{
expr & next = var.values[voffset+it->left.i] > 0.0 ? *it->right.q->left.e : *it->right.q->right.e;
DerivativeFill(next, 0, element, entries,var.values[doffset+k], user_data);
}
break;
case AD_PLUS:
var.values[doffset+it->left.i] += var.values[doffset+k];
var.values[doffset+it->right.i] += var.values[doffset+k];
break;
case AD_MINUS:
var.values[doffset+it->left.i] += var.values[doffset+k];
var.values[doffset+it->right.i] -= var.values[doffset+k];
break;
case AD_MULT:
rval = var.values[voffset+it->right.i];
lval = var.values[voffset+it->left.i];
var.values[doffset+it->left.i] += var.values[doffset+k]*rval;
var.values[doffset+it->right.i] += var.values[doffset+k]*lval;
break;
case AD_DIV:
rval = var.values[voffset+it->right.i];
lval = var.values[voffset+it->left.i];
var.values[doffset+it->left.i] += var.values[doffset+k]/rval;
var.values[doffset+it->right.i] -= var.values[doffset+k]*lval/(rval*rval);
break;
case AD_POW:
ret = var.values[voffset+k];
rval = var.values[voffset+it->right.i];
lval = var.values[voffset+it->left.i];
var.values[doffset+it->right.i] += var.values[doffset+k] * ::log(lval);
var.values[doffset+it->left.i] += var.values[doffset+k] * ret * rval / lval;
break;
case AD_SQRT:
ret = var.values[voffset+k];
var.values[doffset+it->left.i] += 0.5 * var.values[doffset+k] / ret;
break;
case AD_ABS:
lval = var.values[voffset+it->left.i];
var.values[doffset+it->left.i] += var.values[doffset+k] * (lval > 0.0 ? 1.0 : -1.0);
break;
case AD_EXP:
ret = var.values[voffset+k];
var.values[doffset+it->left.i] += var.values[doffset+k] * ret;
break;
case AD_LOG:
lval = var.values[voffset+it->left.i];
var.values[doffset+it->left.i] += var.values[doffset+k] / lval;
break;
case AD_SIN:
lval = var.values[voffset+it->left.i];
var.values[doffset+it->left.i] += var.values[doffset+k] * ::cos(lval);
break;
case AD_COS:
lval = var.values[voffset+it->left.i];
var.values[doffset+it->left.i] -= var.values[doffset+k] * ::sin(lval);
break;
case AD_COND_MARK: break;
case AD_COND_TYPE: break;
case AD_CONST: break;
case AD_MES: break;
case AD_VAL:
{
expr & next = *it->right.e;
next.current_stencil.resize(1);
next.current_stencil[0] = stencil_pair(e->GetHandle(),1.0);
next.resize_for_stencil();
var.values[doffset+it->left.i] += var.values[doffset+k] * EvaluateSub(next,0,element,user_data);
break;
}
default:
if (it->op >= AD_FUNC) {}
else if (it->op >= AD_TABLE)
{
lval = var.values[voffset+it->left.i];
var.values[doffset+it->left.i] += var.values[doffset+k] * reg_tables[it->op-AD_TABLE]->get_derivative(lval);
}
else if (it->op >= AD_STNCL)
{
for (INMOST_DATA_ENUM_TYPE j = 0; j < it->left.e->current_stencil.size(); ++j) if( it->left.e->current_stencil[j].first != NULL )
{
DerivativeFill(*it->left.e, j, ENUMUNDEF, entries, var.values[doffset+k] * it->left.e->current_stencil[j].second, user_data);
}
}
else if (it->op >= AD_CTAG) {}
else if (it->op >= AD_TAG)
{
if (isDynamicValid(e, it->op))
{
entries[GetDynamicIndex(e,it->op,it->left.i)] += var.values[doffset+k];
}
}
else assert(false);
}
} while(k-- != 0);
}
INMOST_DATA_REAL_TYPE Automatizator::EvaluateSub(expr & var, INMOST_DATA_ENUM_TYPE element, INMOST_DATA_ENUM_TYPE parent, void * user_data)
{
INMOST_DATA_ENUM_TYPE k = 0, offset = var.values_offset(element);
Storage e = Storage(m,var.current_stencil[element].first);
expr::expr_data * arr = &var.data[0], *it;
//for (expr::data_type::iterator it = var.data.begin(); it != var.data.end(); ++it)
do
{
it = arr+k;
switch (it->op)
{
case AD_EXT:
assert(parent != ENUMUNDEF);
var.values[offset+k] = it->left.e->values[it->left.e->values_offset(parent)+it->right.i];
break;
case AD_COND:
{
expr & next = var.values[offset+it->left.i] > 0.0 ? *it->right.q->left.e : *it->right.q->right.e;
next.current_stencil.resize(1);
next.current_stencil[0] = stencil_pair(e->GetHandle(),1.0);
next.resize_for_stencil();
var.values[offset+k] = EvaluateSub(next, 0,element, user_data);
}
break;
case AD_PLUS: var.values[offset+k] = var.values[offset+it->left.i] + var.values[offset+it->right.i]; break;
case AD_MINUS: var.values[offset+k] = var.values[offset+it->left.i] - var.values[offset+it->right.i]; break;
case AD_MULT: var.values[offset+k] = var.values[offset+it->left.i] * var.values[offset+it->right.i]; break;
case AD_DIV: var.values[offset+k] = var.values[offset+it->left.i] / var.values[offset+it->right.i]; break;
case AD_POW: var.values[offset+k] = ::pow(var.values[offset+it->left.i], var.values[offset+it->right.i]); break;
case AD_SQRT: var.values[offset+k] = ::sqrt(var.values[offset+it->left.i]); break;
case AD_ABS: var.values[offset+k] = ::fabs(var.values[offset+it->left.i]); break;
case AD_EXP: var.values[offset+k] = ::exp(var.values[offset+it->left.i]); break;
case AD_LOG: var.values[offset+k] = ::log(var.values[offset+it->left.i]); break;
case AD_SIN: var.values[offset+k] = ::sin(var.values[offset+it->left.i]); break;
case AD_COS: var.values[offset+k] = ::cos(var.values[offset+it->left.i]); break;
case AD_CONST: var.values[offset+k] = it->left.r; break;
case AD_COND_TYPE: var.values[offset+k] = ((e->GetElementType() & it->left.i)? 1.0 : -1.0); break;
case AD_COND_MARK: var.values[offset+k] = e->GetMarker(it->left.i) ? 1.0 : -1.0; break;
case AD_MES: assert(!(e->GetElementType() & (ESET | MESH))); m->GetGeometricData(e->GetHandle(), MEASURE, &var.values[offset+k]); break;
case AD_VAL: var.values[offset+k] = var.values[offset+it->left.i]; break;
default:
if (it->op >= AD_FUNC) var.values[offset+k] = reg_funcs[it->op-AD_FUNC].func(e, user_data);
else if (it->op >= AD_TABLE) var.values[offset+k] = reg_tables[it->op-AD_TABLE]->get_value(var.values[offset+it->left.i]);
else if (it->op >= AD_STNCL)
{
it->left.e->current_stencil.clear();
GetStencil(it->op,e,user_data,it->left.e->current_stencil);
it->left.e->resize_for_stencil();
var.values[offset+k] = 0.0;
for (INMOST_DATA_ENUM_TYPE j = 0; j < it->left.e->current_stencil.size(); ++j) if( it->left.e->current_stencil[j].first != NULL )
var.values[offset+k] += EvaluateSub(*it->left.e,j,ENUMUNDEF,user_data) * it->left.e->current_stencil[j].second;
}
else if (it->op >= AD_CTAG) var.values[offset+k] = GetStaticValue(e, it->op, it->left.i);
else if (it->op >= AD_TAG) var.values[offset+k] = GetDynamicValue(e, it->op, it->left.i);
else assert(false);
}
//k++;
} while(++k != var.data.size());
return var.values[offset+var.data.size()-1];
Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
expr.GetJacobian(1.0,merger);
merger.RetriveRow(entries);
merger.Clear();
}
INMOST_DATA_REAL_TYPE Automatizator::Evaluate(expr & var, const Storage & e, void * user_data)
void AddBasicExpression(Sparse::Row & entries, INMOST_DATA_REAL_TYPE multme, INMOST_DATA_REAL_TYPE multit, const basic_expression & expr)
{
var.current_stencil.resize(1);
var.current_stencil[0] = stencil_pair(e->GetHandle(),1.0);
var.resize_for_stencil();
return EvaluateSub(var,0,ENUMUNDEF,user_data);
Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
merger.PushRow(multme,entries);
expr.GetJacobian(multit,merger);
merger.RetriveRow(entries);
merger.Clear();
}
INMOST_DATA_REAL_TYPE Automatizator::Derivative(expr & var, const Storage & e, Sparse::Row & out, Storage::real multiply, void * user_data)
{
Sparse::RowMerger & m = GetMerger();
INMOST_DATA_REAL_TYPE ret;
var.current_stencil.resize(1);
var.current_stencil[0] = stencil_pair(e->GetHandle(),1.0);
var.resize_for_stencil();
ret = EvaluateSub(var,0,ENUMUNDEF,user_data);
m.PushRow(1.0,out);
DerivativeFill(var, 0, ENUMUNDEF, m, multiply, user_data);
m.RetriveRow(out);
m.Clear();
return ret*multiply;
}
#else
INMOST_DATA_REAL_TYPE Automatizator::DerivativePrecompute(const expr & var, const Storage & e, precomp_values_t & values, void * user_data)
void FromGetJacobian(const basic_expression & expr, INMOST_DATA_REAL_TYPE mult, Sparse::Row & r)
{
assert(var.op != AD_NONE);
INMOST_DATA_REAL_TYPE lval, rval, ret = 0.0;
switch (var.op)
{
case AD_COND_TYPE:
lval = (e->GetElementType() & (*(ElementType*)&var.left))? 1.0 : -1.0;
return lval;
case AD_COND_MARK:
lval = e->GetMarker(*(MarkerType *)&var.left)? 1.0 : -1.0;
return lval;
case AD_COND:
lval = Evaluate(*var.left, e, user_data);
rval = DerivativePrecompute(*(lval > 0.0 ? var.right->left : var.right->right), e, values, user_data);
values.push_back(lval);
return rval*var.coef;
case AD_PLUS:
lval = DerivativePrecompute(*var.left, e, values, user_data);
rval = DerivativePrecompute(*var.right, e, values, user_data);
return (lval + rval)*var.coef;
case AD_MINUS:
lval = DerivativePrecompute(*var.left, e, values, user_data);
rval = DerivativePrecompute(*var.right, e, values, user_data);
return (lval - rval)*var.coef;
case AD_MULT:
lval = DerivativePrecompute(*var.left, e, values, user_data);
rval = DerivativePrecompute(*var.right, e, values, user_data);
values.push_back(lval);
values.push_back(rval);
return (lval * rval)*var.coef;
case AD_DIV:
lval = DerivativePrecompute(*var.left, e, values, user_data);
rval = DerivativePrecompute(*var.right, e, values, user_data);
values.push_back(lval);
values.push_back(rval);
return (lval / rval)*var.coef;
case AD_POW:
lval = DerivativePrecompute(*var.left, e, values, user_data);
rval = DerivativePrecompute(*var.right, e, values, user_data);
values.push_back(lval);
values.push_back(rval);
ret = ::pow(lval, rval);
values.push_back(ret);
return ret*var.coef;
case AD_SQRT:
lval = DerivativePrecompute(*var.left, e, values, user_data);
ret = ::sqrt(lval);
values.push_back(ret);
return ret*var.coef;
case AD_INV:
lval = DerivativePrecompute(*var.left, e, values, user_data);
values.push_back(lval);
return var.coef / lval;
case AD_ABS:
lval = DerivativePrecompute(*var.left, e, values, user_data);
values.push_back(lval);
return ::fabs(lval)*var.coef;
case AD_EXP:
lval = DerivativePrecompute(*var.left, e, values, user_data);
ret = ::exp(lval);
values.push_back(ret);
return ret*var.coef;
case AD_LOG:
lval = DerivativePrecompute(*var.left, e, values, user_data);
values.push_back(lval);
return ::log(lval)*var.coef;
case AD_SIN:
lval = DerivativePrecompute(*var.left, e, values, user_data);
values.push_back(lval);
return ::sin(lval)*var.coef;
case AD_COS:
lval = DerivativePrecompute(*var.left, e, values, user_data);
values.push_back(lval);
return ::cos(lval)*var.coef;
case AD_CONST:
return var.coef;
case AD_MES:
assert(!(e->GetElementType() & (ESET | MESH)));
m->GetGeometricData(e->GetHandle(), MEASURE, &ret);
return ret*var.coef;
case AD_VAL:
lval = DerivativePrecompute(*var.left, e, values, user_data);
rval = Evaluate(*var.right,e,user_data);
values.push_back(rval);
return lval*var.coef;
}
if (var.op >= AD_FUNC)
{
ret = reg_funcs[var.op-AD_FUNC].func(e, user_data);
return ret*var.coef;
}
else if (var.op >= AD_TABLE)
{
lval = DerivativePrecompute(*var.left, e, values, user_data);
values.push_back(lval);
ret = reg_tables[var.op-AD_TABLE]->get_value(lval);
return ret*var.coef;
}
else if (var.op >= AD_STNCL)
{
stencil_kind_domain & st = reg_stencils[var.op-AD_STNCL];
assert(st.domainmask == 0 || e->GetMarker(st.domainmask));
if (st.kind == 0)
{
Storage::reference_array elems = e->ReferenceArray(static_cast<stencil_tag *>(st.link)->elements);
Storage::real_array coefs = e->RealArray(static_cast<stencil_tag *>(st.link)->coefs);
assert(elems.size() == coefs.size());
for (INMOST_DATA_ENUM_TYPE k = 0; k < elems.size(); ++k) if( elems.at(k) != InvalidHandle() )
{
lval = DerivativePrecompute(*var.left, elems[k], values, user_data);
ret += lval * coefs[k];
}
}
else if (st.kind == 1)
{
stencil_pairs get_st;
reinterpret_cast<stencil_callback>(st.link)(e, get_st, user_data);
for (INMOST_DATA_ENUM_TYPE k = 0; k < get_st.size(); ++k) if( get_st[k].first != InvalidHandle() )
{
lval = DerivativePrecompute(*var.left, Storage(m,get_st[k].first), values, user_data);
ret += lval * get_st[k].second;
}
}
return ret*var.coef;
}
else if (var.op >= AD_CTAG)
{
INMOST_DATA_ENUM_TYPE comp = *(INMOST_DATA_ENUM_TYPE *)(&var.left);
ret = GetStaticValue(e, var.op, comp);
return ret * var.coef;
}
else if (var.op >= AD_TAG)
{
INMOST_DATA_ENUM_TYPE comp = *(INMOST_DATA_ENUM_TYPE *)(&var.left);
ret = GetDynamicValue(e, var.op, comp);
return ret*var.coef;
}
assert(false);
return 0.0;
Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
expr.GetJacobian(mult,merger);
merger.AddRow(1.0,r);
merger.RetriveRow(r);
merger.Clear();
}
//! returns offset from the end of precomputed values
void Automatizator::DerivativeFill(const expr & var, const Storage & e, Sparse::RowMerger & entries, precomp_values_t & values, INMOST_DATA_REAL_TYPE multval, void * user_data)
{
assert(var.op != AD_NONE);
INMOST_DATA_REAL_TYPE lval, rval, ret;
switch (var.op)
{
case AD_COND_MARK:
case AD_COND_TYPE:
return;
case AD_COND:
lval = values.back(); values.pop_back();
DerivativeFill(*(lval > 0.0 ? var.right->left : var.right->right), e, entries, values, multval*var.coef, user_data);
return;
case AD_PLUS:
DerivativeFill(*var.right, e, entries, values, multval*var.coef, user_data);
DerivativeFill(*var.left, e, entries, values, multval*var.coef, user_data);
return;
case AD_MINUS:
DerivativeFill(*var.right, e, entries, values, -multval*var.coef, user_data);
DerivativeFill(*var.left, e, entries, values, multval*var.coef, user_data);
return;
case AD_MULT:
rval = values.back(); values.pop_back();
lval = values.back(); values.pop_back();
DerivativeFill(*var.right, e, entries, values, lval*multval*var.coef, user_data);
DerivativeFill(*var.left, e, entries, values, rval*multval*var.coef, user_data);
return;
case AD_DIV:
rval = values.back(); values.pop_back();
lval = values.back(); values.pop_back();
DerivativeFill(*var.right, e, entries, values, -multval * lval / (rval*rval) * var.coef, user_data);
DerivativeFill(*var.left, e, entries, values, multval / rval*var.coef, user_data);
return;
case AD_POW:
ret = values.back(); values.pop_back();
rval = values.back(); values.pop_back();
lval = values.back(); values.pop_back();
DerivativeFill(*var.right, e, entries, values, multval * ret * ::log(lval) * var.coef, user_data);
DerivativeFill(*var.left, e, entries, values, multval * ret * rval / lval * var.coef, user_data);
return;
case AD_SQRT:
ret = values.back(); values.pop_back();
DerivativeFill(*var.left, e, entries, values, 0.5 * multval / ret * var.coef, user_data);
return;
case AD_INV:
lval = values.back(); values.pop_back();
DerivativeFill(*var.left, e, entries, values, -multval / (lval * lval) * var.coef, user_data);
return;
case AD_ABS:
lval = values.back(); values.pop_back();
DerivativeFill(*var.left, e, entries, values, multval * var.coef * (lval > 0 ? 1.0 : -1.0), user_data);
return;
case AD_EXP:
ret = values.back(); values.pop_back();
DerivativeFill(*var.left, e, entries, values, multval * ret * var.coef, user_data);
return;
case AD_LOG:
lval = values.back(); values.pop_back();
DerivativeFill(*var.left, e, entries, values, multval * var.coef / lval, user_data);
return;
case AD_SIN:
lval = values.back(); values.pop_back();
DerivativeFill(*var.left, e, entries, values, multval * var.coef * ::cos(lval), user_data);
return;
case AD_COS:
lval = values.back(); values.pop_back();
DerivativeFill(*var.left, e, entries, values, -multval * var.coef * ::sin(lval), user_data);
return;
case AD_CONST: return;
case AD_MES: return;
case AD_VAL:
rval = values.back(); values.pop_back();
DerivativeFill(*var.left, e, entries, values, multval * var.coef * rval, user_data);
return;
}
if (var.op >= AD_FUNC)
{
return;
}
else if (var.op >= AD_TABLE)
{
lval = values.back(); values.pop_back();
DerivativeFill(*var.left, e, entries, values, multval * var.coef * reg_tables[var.op-AD_TABLE]->get_derivative(lval), user_data);
return;
}
else if (var.op >= AD_STNCL)
{
stencil_kind_domain & st = reg_stencils[var.op-AD_STNCL];
assert(st.domainmask == 0 || e->GetMarker(st.domainmask));
if (st.kind == 0)
{
Storage::reference_array elems = e->ReferenceArray(static_cast<stencil_tag *>(st.link)->elements);
Storage::real_array coefs = e->RealArray(static_cast<stencil_tag *>(st.link)->coefs);
assert(elems.size() == coefs.size());
for (INMOST_DATA_ENUM_TYPE k = elems.size(); k > 0; --k) if( elems.at(k-1) != InvalidHandle() )
DerivativeFill(*var.left, elems[k - 1], entries, values, var.coef * coefs[k - 1] * multval, user_data);
}
else if (st.kind == 1)
{
stencil_pairs get_st;
reinterpret_cast<stencil_callback>(st.link)(e, get_st, user_data);
for (INMOST_DATA_ENUM_TYPE k = static_cast<INMOST_DATA_ENUM_TYPE>(get_st.size()); k > 0; --k) if( get_st[k-1].first != InvalidHandle() )
DerivativeFill(*var.left, Storage(m,get_st[k - 1].first), entries, values, var.coef * get_st[k - 1].second*multval, user_data);
}
return;
}
else if (var.op >= AD_CTAG) return;
else if (var.op >= AD_TAG)
{
if (isDynamicValid(e, var.op))
{
INMOST_DATA_ENUM_TYPE ind = GetDynamicIndex(e, var.op, *(INMOST_DATA_ENUM_TYPE *)(&var.left));
entries[ind] += multval * var.coef;
}
return;
}
assert(false);
return;
}
INMOST_DATA_REAL_TYPE Automatizator::Evaluate(const expr & var, const Storage & e, void * user_data)
{
assert(var.op != AD_NONE);
switch (var.op)
{
case AD_COND_MARK: return e->GetMarker(*(MarkerType *)&var.left) ? 1.0 : -1.0;
case AD_COND_TYPE: return (e->GetElementType() & (*(ElementType *)&var.left)) ? 1.0 : -1.0;
case AD_COND: return Evaluate(*(Evaluate(*var.left, e, user_data) > 0.0 ? var.right->left : var.right->right), e, user_data)*var.coef;
case AD_PLUS: return (Evaluate(*var.left, e, user_data) + Evaluate(*var.right, e, user_data))*var.coef;
case AD_MINUS: return (Evaluate(*var.left, e, user_data) - Evaluate(*var.right, e, user_data))*var.coef;
case AD_MULT: return (Evaluate(*var.left, e, user_data) * Evaluate(*var.right, e, user_data))*var.coef;
case AD_DIV: return (Evaluate(*var.left, e, user_data) / Evaluate(*var.right, e, user_data))*var.coef;
case AD_INV: return var.coef / Evaluate(*var.left, e, user_data);
case AD_POW: return ::pow(Evaluate(*var.left, e, user_data), Evaluate(*var.right, e, user_data))*var.coef;
case AD_SQRT: return ::sqrt(Evaluate(*var.left, e, user_data))*var.coef;
case AD_ABS: return ::fabs(Evaluate(*var.left, e, user_data))*var.coef;
case AD_EXP: return ::exp(Evaluate(*var.left, e, user_data))*var.coef;
case AD_LOG: return ::log(Evaluate(*var.left, e, user_data))*var.coef;
case AD_SIN: return ::sin(Evaluate(*var.left, e, user_data))*var.coef;
case AD_COS: return ::cos(Evaluate(*var.left, e, user_data))*var.coef;
case AD_CONST: return var.coef;
case AD_MES: assert(!(e->GetElementType() & (ESET | MESH))); Storage::real ret; m->GetGeometricData(e->GetHandle(), MEASURE, &ret); return ret*var.coef;
case AD_VAL: return Evaluate(*var.left,e,user_data)*var.coef;
}
if (var.op >= AD_FUNC) return reg_funcs[var.op-AD_FUNC].func(e, user_data);
if (var.op >= AD_TABLE) return reg_tables[var.op-AD_TABLE]->get_value(Evaluate(*var.left, e, user_data))*var.coef;
if (var.op >= AD_STNCL)
{
INMOST_DATA_REAL_TYPE ret = 0.0;
stencil_kind_domain & st = reg_stencils[var.op-AD_STNCL];
assert(st.domainmask == 0 || e->GetMarker(st.domainmask));
if (st.kind == 0)
{
Storage::reference_array elems = e->ReferenceArray(static_cast<stencil_tag *>(st.link)->elements);
Storage::real_array coefs = e->RealArray(static_cast<stencil_tag *>(st.link)->coefs);
assert(elems.size() == coefs.size());
for (INMOST_DATA_ENUM_TYPE k = 0; k < elems.size(); ++k) if( elems.at(k) != InvalidHandle() )
ret += var.coef * Evaluate(*var.left, elems[k], user_data) * coefs[k];
}
else if (st.kind == 1)
{
stencil_pairs get_st;
reinterpret_cast<stencil_callback>(st.link)(e, get_st, user_data);
for (INMOST_DATA_ENUM_TYPE k = 0; k < get_st.size(); ++k) if ( get_st[k].first != InvalidHandle() )
ret += var.coef * Evaluate(*var.left, Storage(m,get_st[k].first), user_data) * get_st[k].second;
}
return ret;
}
if (var.op >= AD_CTAG) return GetStaticValue(e, var.op, *(INMOST_DATA_ENUM_TYPE *)(&var.left))*var.coef;
if (var.op >= AD_TAG) return GetDynamicValue(e, var.op, *(INMOST_DATA_ENUM_TYPE *)(&var.left))*var.coef;
assert(false);
return 0.0;
}
INMOST_DATA_REAL_TYPE Automatizator::Derivative(const expr & var, const Storage & e, Sparse::Row & out, Storage::real multiply, void * user_data)
{
Sparse::RowMerger & m = GetMerger();
INMOST_DATA_REAL_TYPE ret;
precomp_values_t values;
ret = DerivativePrecompute(var, e, values, user_data);
m.PushRow(1.0,out);
DerivativeFill(var, e, m, values, multiply, user_data);
m.RetriveRow(out);
m.Clear();
return ret*multiply;
}
#endif
Automatizator::Automatizator(Mesh * m) :first_num(0), last_num(0), m(m) {}
Automatizator::~Automatizator()
{
for (unsigned k = 0; k < index_tags.size(); k++)
index_tags[k].indices = m->DeleteTag(index_tags[k].indices);
for (table_type::iterator it = reg_tables.begin(); it != reg_tables.end(); ++it)
{
delete[] (*it)->args;
delete[] (*it)->vals;
delete (*it);
}
for (stencil_type::iterator it = reg_stencils.begin(); it != reg_stencils.end(); ++it)
if (it->kind == 0)
delete static_cast<stencil_tag *>(it->link);
}
INMOST_DATA_ENUM_TYPE Automatizator::RegisterFunc(std::string name, func_callback func)
{
INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_funcs.size()) + AD_FUNC;
func_name_callback v;
v.name = name;
v.func = func;
reg_funcs.push_back(v);
return ret;
}
//register stencil that can be got from tags
INMOST_DATA_ENUM_TYPE Automatizator::RegisterStencil(std::string name, Tag elements_tag, Tag coefs_tag, MarkerType domain_mask)
{
INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_stencils.size()) + AD_STNCL;
stencil_kind_domain st;
stencil_tag * save = new stencil_tag;
st.name = name;
save->coefs = coefs_tag;
save->elements = elements_tag;
st.kind = 0;
st.link = static_cast<void *>(save);
st.domainmask = domain_mask;
reg_stencils.push_back(st);
return ret;
}
//register stencil that can be got from function
INMOST_DATA_ENUM_TYPE Automatizator::RegisterStencil(std::string name, stencil_callback func, MarkerType domain_mask)
{
INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_stencils.size()) + AD_STNCL;
stencil_kind_domain st;
st.name = name;
st.kind = 1;
st.link = reinterpret_cast<void *>(func);
st.domainmask = domain_mask;
reg_stencils.push_back(st);
return ret;
}
INMOST_DATA_ENUM_TYPE Automatizator::RegisterTable(std::string name, INMOST_DATA_REAL_TYPE * Arguments, INMOST_DATA_REAL_TYPE * Values, INMOST_DATA_ENUM_TYPE size)
{
INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_tables.size()) + AD_TABLE;
table_ptr t = new table;
t->name = name;