From 173238768cd8a3a7b4468a3ea1f3722c79147800 Mon Sep 17 00:00:00 2001 From: kabdelhak <38032125+kabdelhak@users.noreply.github.com> Date: Tue, 8 Mar 2022 10:13:03 +0100 Subject: [PATCH] [NB] add function and algorithm differentiation (#8644) * [NB] add function and algorithm differentiation - algorithms and statements are differentiated properly - functions and function bodys are differentiated - recursion is handled - multiple occurences are handled - ToDo: strip interface variables such that only continuous are differentiated * [NB] improve function differentiation 1. if the function is builtin -> use hardcoded logic 2. if the function is not builtin -> check if there is a 'fitting' derivative defined. - 'fitting' means that all the zeroDerivative annotations have to hold 2.1 fitting function found -> use it 2.2 fitting function not found -> differentiate the body of the function ToDo: respect the 'order' of the derivative when differentiating! --- .../Compiler/NBackEnd/Classes/NBVariable.mo | 58 ++- .../NBackEnd/Modules/3_Post/NBJacobian.mo | 2 +- .../NBackEnd/Modules/3_Post/NBSolve.mo | 1 - .../Compiler/NBackEnd/Util/NBBackendUtil.mo | 19 + .../Compiler/NBackEnd/Util/NBDifferentiate.mo | 356 +++++++++++++++--- OMCompiler/Compiler/NFFrontEnd/NFComponent.mo | 9 +- OMCompiler/Compiler/NFFrontEnd/NFFunction.mo | 30 +- .../NFFrontEnd/NFFunctionDerivative.mo | 33 +- OMCompiler/Compiler/NFFrontEnd/NFInstNode.mo | 32 +- .../modelica/NBackend/functions/Makefile | 3 +- .../functions/function_annotation_der.mos | 110 +++--- .../NBackend/functions/function_diff.mos | 86 +++++ 12 files changed, 577 insertions(+), 162 deletions(-) create mode 100644 testsuite/simulation/modelica/NBackend/functions/function_diff.mos diff --git a/OMCompiler/Compiler/NBackEnd/Classes/NBVariable.mo b/OMCompiler/Compiler/NBackEnd/Classes/NBVariable.mo index c2b7e488096..350366c1d56 100644 --- a/OMCompiler/Compiler/NBackEnd/Classes/NBVariable.mo +++ b/OMCompiler/Compiler/NBackEnd/Classes/NBVariable.mo @@ -305,11 +305,11 @@ public output Boolean b; algorithm b := match Pointer.access(var) - case Variable.VARIABLE(backendinfo = BackendExtension.BACKEND_INFO(varKind = BackendExtension.DISCRETE_STATE())) then false; - case Variable.VARIABLE(backendinfo = BackendExtension.BACKEND_INFO(varKind = BackendExtension.DISCRETE())) then false; - case Variable.VARIABLE(backendinfo = BackendExtension.BACKEND_INFO(varKind = BackendExtension.PREVIOUS())) then false; - case Variable.VARIABLE(backendinfo = BackendExtension.BACKEND_INFO(varKind = BackendExtension.PARAMETER())) then false; - case Variable.VARIABLE(backendinfo = BackendExtension.BACKEND_INFO(varKind = BackendExtension.CONSTANT())) then false; + case Variable.VARIABLE(backendinfo = BackendExtension.BACKEND_INFO(varKind = BackendExtension.DISCRETE_STATE())) then false; + case Variable.VARIABLE(backendinfo = BackendExtension.BACKEND_INFO(varKind = BackendExtension.DISCRETE())) then false; + case Variable.VARIABLE(backendinfo = BackendExtension.BACKEND_INFO(varKind = BackendExtension.PREVIOUS())) then false; + case Variable.VARIABLE(backendinfo = BackendExtension.BACKEND_INFO(varKind = BackendExtension.PARAMETER())) then false; + case Variable.VARIABLE(backendinfo = BackendExtension.BACKEND_INFO(varKind = BackendExtension.CONSTANT())) then false; else true; end match; end isContinuous; @@ -818,21 +818,19 @@ public _ := match ComponentRef.node(cref) local InstNode qual; - Pointer old_var_ptr; Variable var; - case qual as InstNode.VAR_NODE() - algorithm - // get the variable pointer from the old cref to later on link back to it - old_var_ptr := BVariable.getVarPointer(cref); - // prepend the seed str and the matrix name and create the new cref_DIFF_DIFF - qual.name := PARTIAL_DERIVATIVE_STR + "_" + name; - cref := ComponentRef.append(cref, ComponentRef.fromNode(qual, ComponentRef.scalarType(cref))); - var := fromCref(cref); - // update the variable to be a jac var and pass the pointer to the original variable - // ToDo: tmps will get JAC_DIFF_VAR ! - var.backendinfo := BackendExtension.BackendInfo.setVarKind(var.backendinfo, BackendExtension.JAC_VAR()); - // create the new variable pointer and safe it to the component reference - (var_ptr, cref) := makeVarPtrCyclic(var, cref); + + // regular case for jacobians + case qual as InstNode.VAR_NODE() algorithm + // prepend the seed str and the matrix name and create the new cref_DIFF_DIFF + qual.name := PARTIAL_DERIVATIVE_STR + "_" + name; + cref := ComponentRef.append(cref, ComponentRef.fromNode(qual, ComponentRef.scalarType(cref))); + var := fromCref(cref); + // update the variable to be a jac var and pass the pointer to the original variable + // ToDo: tmps will get JAC_DIFF_VAR ! + var.backendinfo := BackendExtension.BackendInfo.setVarKind(var.backendinfo, BackendExtension.JAC_VAR()); + // create the new variable pointer and safe it to the component reference + (var_ptr, cref) := makeVarPtrCyclic(var, cref); then (); else algorithm @@ -841,6 +839,28 @@ public end match; end makePDerVar; + function makeFDerVar + "Creates a function derivative cref. Used in NBDifferentiation + for differentiating body vars of a function." + input output ComponentRef cref "old component reference to new component reference"; + algorithm + cref := match ComponentRef.node(cref) + local + InstNode qual; + + // for function differentiation (crefs are not lowered and only known locally) + case qual as InstNode.COMPONENT_NODE() algorithm + // prepend the seed str, matrix name locally not needed + qual.name := FUNCTION_DERIVATIVE_STR + "_" + qual.name; + cref := ComponentRef.fromNode(qual, ComponentRef.scalarType(cref)); + then cref; + + else algorithm + Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed for " + ComponentRef.toString(cref)}); + then fail(); + end match; + end makeFDerVar; + function makeStartVar "Creates a start variable pointer from a cref. Used in NBInitialization. e.g: angle -> $START.angle" diff --git a/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBJacobian.mo b/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBJacobian.mo index d284a0c587e..bc193adcd17 100644 --- a/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBJacobian.mo +++ b/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBJacobian.mo @@ -537,7 +537,6 @@ protected jacobianHT = optHT, // seed and temporary cref hashtable diffType = NBDifferentiate.DifferentiationType.JACOBIAN, funcTree = funcTree, - diffedFunctions = AvlSetPath.new(), scalarized = seedCandidates.scalarized ); @@ -546,6 +545,7 @@ protected (diffed_eqn_lst, diffArguments) := Differentiate.differentiateEquationPointerList(listReverse(eqn_lst), diffArguments, idx, name, getInstanceName()); diffedEquations := EquationPointers.fromList(diffed_eqn_lst); + funcTree := diffArguments.funcTree; // create equation data for jacobian // ToDo: split temporary and auxiliares once tearing is applied diff --git a/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBSolve.mo b/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBSolve.mo index a5ac3ede025..dbe29fbd2d5 100644 --- a/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBSolve.mo +++ b/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBSolve.mo @@ -385,7 +385,6 @@ public jacobianHT = NONE(), diffType = NBDifferentiate.DifferentiationType.SIMPLE, funcTree = funcTree, - diffedFunctions = AvlSetPath.new(), scalarized = false ); (derivative, diffArgs) := Differentiate.differentiateExpressionDump(residual, diffArgs, getInstanceName()); diff --git a/OMCompiler/Compiler/NBackEnd/Util/NBBackendUtil.mo b/OMCompiler/Compiler/NBackEnd/Util/NBBackendUtil.mo index dca6bf592b5..226965eb5fc 100644 --- a/OMCompiler/Compiler/NBackEnd/Util/NBBackendUtil.mo +++ b/OMCompiler/Compiler/NBackEnd/Util/NBBackendUtil.mo @@ -638,5 +638,24 @@ public end if; end isOnlyTimeDependentFold; + function isContinuous + input Expression exp; + output Boolean b; + algorithm + b := Expression.fold(exp, isContinuousFold, true); + end isContinuous; + + function isContinuousFold + input Expression exp; + input output Boolean b; + algorithm + if b then + b := match exp + case Expression.CREF() then BVariable.checkCref(exp.cref, BVariable.isContinuous); + else true; + end match; + end if; + end isContinuousFold; + annotation(__OpenModelica_Interface="backend"); end NBBackendUtil; diff --git a/OMCompiler/Compiler/NBackEnd/Util/NBDifferentiate.mo b/OMCompiler/Compiler/NBackEnd/Util/NBDifferentiate.mo index 2cf431ec24e..d87c6da9fc9 100644 --- a/OMCompiler/Compiler/NBackEnd/Util/NBDifferentiate.mo +++ b/OMCompiler/Compiler/NBackEnd/Util/NBDifferentiate.mo @@ -36,36 +36,39 @@ encapsulated package NBDifferentiate " public // OF imports + import Absyn; + import AbsynUtil; import BaseAvlTree; import AvlSetPath; import DAE; // NF imports + import Algorithm = NFAlgorithm; import BuiltinFuncs = NFBuiltinFuncs; import Call = NFCall; + import Class = NFClass; import ComponentRef = NFComponentRef; import Expression = NFExpression; - import NFInstNode.InstNode; - import NFFunction.Function; + import NFInstNode.{InstNode, CachedData}; import NFFlatten.{FunctionTree, FunctionTreeImpl}; + import NFFunction.Function; + import FunctionDerivative = NFFunctionDerivative; import Operator = NFOperator; import Prefixes = NFPrefixes; + import Sections = NFSections; import SimplifyExp = NFSimplifyExp; + import Statement = NFStatement; import Type = NFType; import NFPrefixes.Variability; import Variable = NFVariable; // Backend imports - import NBEquation.Equation; - import NBEquation.EquationAttributes; - import NBEquation.EquationPointers; - import NBEquation.IfEquationBody; - import NBEquation.WhenEquationBody; - import NBEquation.WhenStatement; + import NBEquation.{Equation, EquationAttributes, EquationPointers, IfEquationBody, WhenEquationBody, WhenStatement}; import BVariable = NBVariable; // Util imports import Array; + import BackendUtil = NBBackendUtil; import Error; import UnorderedMap; @@ -81,7 +84,6 @@ public Option> jacobianHT "seed and temporary cref hashtable x --> $SEED.MATRIX.x, y --> $pDer.MATRIX.y"; DifferentiationType diffType "Differentiation use case (time, simple, function, jacobian)"; FunctionTree funcTree "Function tree containing all functions and their known derivatives"; - AvlSetPath.Tree diffedFunctions "current functions, to prevent recursive differentiation"; Boolean scalarized "true if the variables are scalarized"; end DIFFERENTIATION_ARGUMENTS; @@ -94,7 +96,6 @@ public jacobianHT = NONE(), diffType = ty, funcTree = funcTree, - diffedFunctions = AvlSetPath.new(), scalarized = false ); end default; @@ -107,7 +108,7 @@ public function differentiateEquationPointerList "author: kabdelhak - Differentiates an array of equations wrapped in pointers." + Differentiates a list of equations wrapped in pointers." input output list> equations; input output DifferentiationArguments diffArguments; input Pointer idx; @@ -116,7 +117,6 @@ public protected Pointer diffArguments_ptr = Pointer.create(diffArguments); algorithm - // don't use EquationPointers.map because that would manipulate original eqn pointers equations := List.map(equations, function differentiateEquationPointer(diffArguments_ptr = diffArguments_ptr, name = name)); for eqn in equations loop Equation.createName(eqn, idx, context); @@ -177,6 +177,7 @@ public WhenEquationBody whenBody; Pointer diffArguments_ptr; EquationAttributes attr; + Algorithm alg; // ToDo: Element source stuff (see old backend) case Equation.SCALAR_EQUATION() algorithm @@ -231,21 +232,18 @@ public attr := differentiateEquationAttributes(eq.attr, diffArguments); then (Equation.WHEN_EQUATION(eq.size, whenBody, eq.source, attr), diffArguments); + case Equation.ALGORITHM() algorithm + (alg, diffArguments) := differentiateAlgorithm(eq.alg, diffArguments); + then (Equation.ALGORITHM(eq.size, alg, eq.source, eq.expand, eq.attr), diffArguments); + else algorithm // maybe add failtrace here and allow failing Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed for: " + Equation.toString(eq)}); then fail(); end match; - /* ToDo: - record ALGORITHM - Integer size "output size"; - Algorithm alg "Algorithm statements"; - DAE.ElementSource source "origin of algorithm"; - DAE.Expand expand "this algorithm was translated from an equation. we should not expand array crefs!"; - EquationAttributes attr "Additional Attributes"; - end ALGORITHM; +/* ToDo record AUX_EQUATION "Auxiliary equations are generated when auxiliary variables are generated that are known to always be solved in this specific equation. E.G. $CSE @@ -480,6 +478,8 @@ public algorithm // extract var pointer first to have following code more readable var_ptr := match exp + // function body expressions are not lowered (maybe do it?) + case _ guard(diffArguments.diffType == DifferentiationType.FUNCTION) then Pointer.create(NBVariable.DUMMY_VARIABLE); case Expression.CREF() then BVariable.getVarPointer(exp.cref); else algorithm Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed for: " + Expression.toString(exp)}); @@ -491,6 +491,30 @@ public Expression res; UnorderedMap jacobianHT; + + // ------------------------------------- + // Special rules for Type: FUNCTION + // (needs to be first because var_ptr is DUMMY) + // ------------------------------------- + + // Types: (FUNCTION) + // Any variable that is in the HT will be differentiated accordingly. 0 otherwise + case (Expression.CREF(), DifferentiationType.FUNCTION, SOME(jacobianHT)) algorithm + strippedCref := ComponentRef.stripSubscriptsAll(exp.cref); + if UnorderedMap.contains(strippedCref, jacobianHT) then + // get the derivative an reapply subscripts + derCref := UnorderedMap.getOrFail(strippedCref, jacobianHT); + derCref := ComponentRef.setSubscriptsList(listReverse(ComponentRef.subscriptsAll(exp.cref)), derCref); + res := Expression.fromCref(derCref); + else + res := Expression.makeZero(exp.ty); + end if; + then (res, diffArguments); + + // ------------------------------------- + // Generic Rules + // ------------------------------------- + // Types: (TIME) // differentiate time cref => 1 case (Expression.CREF(), DifferentiationType.TIME, _) @@ -567,12 +591,6 @@ public var_ptr := BVariable.makeStateVar(var_ptr, der_ptr); then (Expression.fromCref(derCref), diffArguments); - // ------------------------------------- - // Special rules for Type: FUNCTION - // ------------------------------------- - - // ToDo: Types (FUNCTION) all of this! - // ------------------------------------- // Special rules for Type: JACOBIAN // ------------------------------------- @@ -614,10 +632,14 @@ public end match; end differentiateComponentRef; - - // TODO: Copy Differentiate.mo function differentiateCalls function differentiateCall - "Differentiate builtin function calls" + "Differentiate builtin function calls + 1. if the function is builtin -> use hardcoded logic + 2. if the function is not builtin -> check if there is a 'fitting' derivative defined. + - 'fitting' means that all the zeroDerivative annotations have to hold + 2.1 fitting function found -> use it + 2.2 fitting function not found -> differentiate the body of the function + ToDo: respect the 'order' of the derivative when differentiating!" input output Expression exp "Has to be Expression.CALL()"; input output DifferentiationArguments diffArguments; protected @@ -630,13 +652,20 @@ public (exp, diffArguments) := match exp local Expression ret; - Boolean has_derviative_annotation = false; Call call, der_call; - Option func_opt; + Option func_opt, der_func_opt; list derivatives; Function func, der_func; list arguments = {}; Operator addOp, mulOp; + list> arguments_inputs; + Expression arg; + InstNode inp; + Boolean isCont, isReal; + // interface map. If the map contains a variable it has a zero derivative + // if the value is "true" it has to be stripped from the interface + // (it is possible that a variable has a zero derivative, but still appears in the interface) + UnorderedMap interface_map = UnorderedMap.new(stringHashDjb2Mod, stringEqual); // builtin functions case Expression.CALL(call = call as Call.TYPED_CALL()) guard(Function.isBuiltin(call.fn)) algorithm @@ -648,25 +677,44 @@ public func_opt := FunctionTreeImpl.getOpt(diffArguments.funcTree, call.fn.path); if Util.isSome(func_opt) then // The function is in the function tree - derivatives := Function.getDerivatives(Util.getOption(func_opt)); - ret := match derivatives - // no derivative of order 1 defined -> differentiate and add to the function tree - case {} algorithm - Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed because function differentiation is not yet implemented."}); - then fail(); - - // exactly one derivative of order 1 -> use it - case {der_func} algorithm - (arguments, diffArguments) := List.mapFold(call.arguments, differentiateExpression, diffArguments); - then Expression.CALL(Call.makeTypedCall(der_func, listAppend(call.arguments, arguments), call.var, call.purity)); - - // ERROR - more than one derivative of order 1 defined - // TODO pick first one according to MLS 3.5 section 12.7.1 - else algorithm - Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed because there were " + intString(listLength(derivatives)) - + " derivatives of order 1 (expected is exactly one).\n Derivatives:" + List.toString(derivatives, function Function.signatureString(printTypes = true), "", "", "\n", "")}); - then fail(); - end match; + SOME(func) := func_opt; + + // build interface map to check if a function fits + // save all inputs that would end up in a zero derivative in a map + arguments_inputs := List.zip(call.arguments, func.inputs); + for tpl in arguments_inputs loop + (arg, inp) := tpl; + // do not check for continuous if it is for functions (differentiating a function inside a function) + // crefs are not lowered there! assume it is continuous + isCont := (diffArguments.diffType == DifferentiationType.FUNCTION) or BackendUtil.isContinuous(arg); + isReal := Type.isRealRecursive(Expression.typeOf(arg)); // ToDo also records + if not (isCont and isReal) then + // add to map; if it is not continuous also already set to true (always removed from interface) + UnorderedMap.add(InstNode.name(inp), not isCont, interface_map); + end if; + end for; + + // try to get a fitting function from derivatives -> if none is found, differentiate + der_func_opt := Function.getDerivative(func, interface_map); + if Util.isSome(der_func_opt) then + SOME(der_func) := der_func_opt; + else + (der_func, diffArguments) := differentiateFunction(func, interface_map, diffArguments); + end if; + + for tpl in listReverse(arguments_inputs) loop + (arg, inp) := tpl; + // only keep the arguments which are not in the map or have value false + if not(UnorderedMap.contains(InstNode.name(inp), interface_map) and UnorderedMap.getSafe(InstNode.name(inp), interface_map)) then + arguments := arg :: arguments; + end if; + end for; + + // differentiate type arguments and append to original ones + (arguments, diffArguments) := List.mapFold(arguments, differentiateExpression, diffArguments); + arguments := listAppend(call.arguments, arguments); + + ret := Expression.CALL(Call.makeTypedCall(der_func, arguments, call.var, call.purity)); else // The function is not in the function tree and not builtin -> error Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() @@ -949,6 +997,214 @@ public end match; end differentiateBuiltinCall2Arg; + function differentiateFunction + input Function func; + output Function der_func; + input UnorderedMap interface_map; + input output DifferentiationArguments diffArguments; + algorithm + der_func := match func + local + InstNode node; + Pointer cls; + Class new_cls; + DifferentiationArguments funcDiffArgs; + UnorderedMap diff_map = UnorderedMap.new(ComponentRef.hash, ComponentRef.isEqual); + Sections sections; + list algorithms; + Absyn.Path new_path; + FunctionDerivative funcDer; + Function dummy_func; + CachedData cachedData; + + case der_func as Function.FUNCTION(node = node as InstNode.CLASS_NODE(cls = cls)) algorithm + new_cls := match Pointer.access(cls) + case new_cls as Class.INSTANCED_CLASS(sections = sections as Sections.SECTIONS()) algorithm + + // differentiate interface arguments + der_func.inputs := differentiateFunctionInterfaceNodes(der_func.inputs, interface_map, diff_map, true); + der_func.locals := differentiateFunctionInterfaceNodes(der_func.locals, interface_map, diff_map, true); + der_func.locals := listAppend(der_func.locals, list(InstNode.setComponentDirection(NFPrefixes.Direction.NONE, node) for node in der_func.outputs)); + der_func.outputs := differentiateFunctionInterfaceNodes(der_func.outputs, interface_map, diff_map, false); + + // prepare differentiation arguments + funcDiffArgs := DifferentiationArguments.default(); + funcDiffArgs.diffType := DifferentiationType.FUNCTION; + funcDiffArgs.funcTree := diffArguments.funcTree; + funcDiffArgs.jacobianHT := SOME(diff_map); + + // create "fake" function with correct interface to have the interface + // in the case of recursive differentiation (e.g. function calls itself) + dummy_func := func; + node.cls := Pointer.create(new_cls); + node.name := NBVariable.FUNCTION_DERIVATIVE_STR + "." + node.name; + // create "fake" function from new node (update cache to get correct derivative name) + der_func.path := AbsynUtil.prefixPath(NBVariable.FUNCTION_DERIVATIVE_STR, der_func.path); + cachedData := CachedData.FUNCTION({der_func}, true, false); + der_func.node := InstNode.setFuncCache(node, cachedData); + + // create fake derivative + funcDer := FunctionDerivative.FUNCTION_DER( + derivativeFn = der_func.node, + derivedFn = dummy_func.node, + order = Expression.INTEGER(1), + conditions = {}, // possibly needs updating + lowerOrderDerivatives = {} // possibly needs updating + ); + + // add fake derivative to function tree + dummy_func.derivatives := funcDer :: dummy_func.derivatives; + funcDiffArgs.funcTree := FunctionTreeImpl.add(funcDiffArgs.funcTree, dummy_func.path, dummy_func, FunctionTreeImpl.addConflictReplace); + + // differentiate function statements + (algorithms, funcDiffArgs) := List.mapFold(sections.algorithms, differentiateAlgorithm, funcDiffArgs); + + // add them to new node + sections.algorithms := algorithms; + new_cls.sections := sections; + node.cls := Pointer.create(new_cls); + cachedData := CachedData.FUNCTION({der_func}, true, false); + der_func.node := InstNode.setFuncCache(node, cachedData); + der_func.derivatives := {}; + then new_cls; + + else algorithm + Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed for class " + Class.toFlatString(Pointer.access(cls), func.node) + "."}); + then fail(); + end match; + + // add function to function tree + diffArguments.funcTree := FunctionTreeImpl.add(diffArguments.funcTree, der_func.path, der_func); + // add new function as derivative to original function + funcDer := FunctionDerivative.FUNCTION_DER( + derivativeFn = der_func.node, + derivedFn = func.node, + order = Expression.INTEGER(1), + conditions = {}, // possibly needs updating + lowerOrderDerivatives = {} // possibly needs updating + ); + func.derivatives := funcDer :: func.derivatives; + diffArguments.funcTree := FunctionTreeImpl.add(diffArguments.funcTree, func.path, func, FunctionTreeImpl.addConflictReplace); + then der_func; + + else algorithm + Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed for uninstanced function " + Function.signatureString(func) + "."}); + then fail(); + end match; + if Flags.isSet(Flags.DEBUG_DIFFERENTIATION) then + print("\n[BEFORE] " + Function.toFlatString(func) + "\n"); + print("\n[AFTER] " + Function.toFlatString(der_func) + "\n\n"); + end if; + end differentiateFunction; + + function differentiateFunctionInterfaceNodes + "differentiates function interface nodes (inputs, outputs, locals) and + adds them to the diff_map used for differentiation. Also returns the new + interface node lists for the differentiated function. + (outputs only have the differentiated and not the original interface nodes)" + input output list interface_nodes; + input UnorderedMap interface_map; + input UnorderedMap diff_map; + input Boolean keepOld; + protected + list new_nodes; + ComponentRef cref, diff_cref; + algorithm + new_nodes := if keepOld then listReverse(interface_nodes) else {}; + interface_nodes := list(node for node guard(not UnorderedMap.contains(InstNode.name(node), interface_map)) in interface_nodes); + for node in interface_nodes loop + cref := ComponentRef.fromNode(node, InstNode.getType(node)); + diff_cref := BVariable.makeFDerVar(cref); + UnorderedMap.add(cref, diff_cref, diff_map); + new_nodes := ComponentRef.node(diff_cref) :: new_nodes; + end for; + interface_nodes := listReverse(new_nodes); + end differentiateFunctionInterfaceNodes; + + function differentiateAlgorithm + input output Algorithm alg; + input output DifferentiationArguments diffArguments; + protected + list> statements; + list statements_flat; + list inputs, outputs; + algorithm + (statements, diffArguments) := List.mapFold(alg.statements, differentiateStatement, diffArguments); + statements_flat := List.flatten(statements); + (inputs, outputs) := Algorithm.getInputsOutputs(statements_flat); + alg := Algorithm.ALGORITHM(statements_flat, inputs, outputs, alg.source); + end differentiateAlgorithm; + + function differentiateStatement + input Statement stmt; + output list diff_stmts "two statements for 'Real' assignments (diff; original) and else one"; + input output DifferentiationArguments diffArguments; + algorithm + diff_stmts := match stmt + local + Statement diff_stmt; + Expression exp, lhs, rhs; + list branch_stmts_flat; + list> branch_stmts; + list>> branches = {}; + + // I. differentiate 'Real' assignment and return differentiated and original statement + case diff_stmt as Statement.ASSIGNMENT() guard(Type.isRealRecursive(Expression.typeOf(diff_stmt.lhs))) algorithm + (lhs, diffArguments) := differentiateExpression(diff_stmt.lhs, diffArguments); + (rhs, diffArguments) := differentiateExpression(diff_stmt.rhs, diffArguments); + diff_stmt.lhs := lhs; + diff_stmt.rhs := SimplifyExp.simplify(rhs); + then {diff_stmt, stmt}; + + // II. delegate differentiation to body and only return differentiated statement + case diff_stmt as Statement.FOR() algorithm + (branch_stmts, diffArguments) := List.mapFold(diff_stmt.body, differentiateStatement, diffArguments); + diff_stmt.body := List.flatten(branch_stmts); + then {diff_stmt}; + + case diff_stmt as Statement.WHILE() algorithm + (branch_stmts, diffArguments) := List.mapFold(diff_stmt.body, differentiateStatement, diffArguments); + diff_stmt.body := List.flatten(branch_stmts); + then {diff_stmt}; + + case diff_stmt as Statement.FAILURE() algorithm + (branch_stmts, diffArguments) := List.mapFold(diff_stmt.body, differentiateStatement, diffArguments); + diff_stmt.body := List.flatten(branch_stmts); + then {diff_stmt}; + + case diff_stmt as Statement.IF() algorithm + for branch in diff_stmt.branches loop + (exp, branch_stmts_flat) := branch; + (branch_stmts, diffArguments) := List.mapFold(branch_stmts_flat, differentiateStatement, diffArguments); + branches := (exp, List.flatten(branch_stmts)) :: branches; + end for; + diff_stmt.branches := listReverse(branches); + then {diff_stmt}; + + case diff_stmt as Statement.WHEN() algorithm + for branch in diff_stmt.branches loop + (exp, branch_stmts_flat) := branch; + (branch_stmts, diffArguments) := List.mapFold(branch_stmts_flat, differentiateStatement, diffArguments); + branches := (exp, List.flatten(branch_stmts)) :: branches; + end for; + diff_stmt.branches := listReverse(branches); + then {diff_stmt}; + + // III. assignments of non-Real are not differentiated, as well as empty statements + case Statement.ASSIGNMENT() then {stmt}; + case Statement.FUNCTION_ARRAY_INIT() then {stmt}; + case Statement.ASSERT() then {stmt}; + case Statement.TERMINATE() then {stmt}; + case Statement.NORETCALL() then {stmt}; + case Statement.RETURN() then {stmt}; + case Statement.BREAK() then {stmt}; + + else algorithm + Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed for:" + Statement.toString(stmt)}); + then fail(); + end match; + end differentiateStatement; + function differentiateBinary "Some of this is depcreated because of Expression.MULTARY(). Will always try to convert to MULTARY whenever possible. (commutativity)" diff --git a/OMCompiler/Compiler/NFFrontEnd/NFComponent.mo b/OMCompiler/Compiler/NFFrontEnd/NFComponent.mo index 584c4ed5b04..ce7481d4ae7 100644 --- a/OMCompiler/Compiler/NFFrontEnd/NFComponent.mo +++ b/OMCompiler/Compiler/NFFrontEnd/NFComponent.mo @@ -595,29 +595,30 @@ public output Boolean isInput = direction(component) == Direction.INPUT; end isInput; - function makeInput + function setDirection input output Component component; + input Direction direction; protected Attributes attr; algorithm () := match component case UNTYPED_COMPONENT(attributes = attr) algorithm - attr.direction := Direction.INPUT; + attr.direction := direction; component.attributes := attr; then (); case TYPED_COMPONENT(attributes = attr) algorithm - attr.direction := Direction.INPUT; + attr.direction := direction; component.attributes := attr; then (); else (); end match; - end makeInput; + end setDirection; function isOutput input Component component; diff --git a/OMCompiler/Compiler/NFFrontEnd/NFFunction.mo b/OMCompiler/Compiler/NFFrontEnd/NFFunction.mo index e964449d813..10cce72cb59 100644 --- a/OMCompiler/Compiler/NFFrontEnd/NFFunction.mo +++ b/OMCompiler/Compiler/NFFrontEnd/NFFunction.mo @@ -2570,17 +2570,31 @@ protected function getLocalDependenciesExp(locals = locals), deps); end getLocalDependenciesDim; - function getDerivatives - "returns all derivatives of a function of a certain order." + function getDerivative + "returns the first derivative that fits the interface_map + and returns NONE() if none of them fit." input Function original; - input Integer order = 1; - output list derivatives; + input UnorderedMap interface_map; + output Option derivative = NONE(); protected - list order_derivatives; + list derivatives; + Boolean perfect_fit; algorithm - order_derivatives := list(func for func guard(FunctionDerivative.getOrder(func) == order) in original.derivatives); - derivatives := List.flatten(list(getCachedFuncs(func.derivativeFn) for func in order_derivatives)); - end getDerivatives; + for func in original.derivatives loop + if FunctionDerivative.perfectFit(func, interface_map) then + derivative := SOME(List.first(getCachedFuncs(func.derivativeFn))); + return; + end if; + end for; + + // if no derivative could be found, set whole map to true + // so that the self-generated function removes as much as possible + if Util.isSome(derivative) then + for key in UnorderedMap.keyList(interface_map) loop + UnorderedMap.add(key, true, interface_map); + end for; + end if; + end getDerivative; end Function; annotation(__OpenModelica_Interface="frontend"); diff --git a/OMCompiler/Compiler/NFFrontEnd/NFFunctionDerivative.mo b/OMCompiler/Compiler/NFFrontEnd/NFFunctionDerivative.mo index 7f6ae3cb7ac..d8d124d7922 100644 --- a/OMCompiler/Compiler/NFFrontEnd/NFFunctionDerivative.mo +++ b/OMCompiler/Compiler/NFFrontEnd/NFFunctionDerivative.mo @@ -179,18 +179,29 @@ public subMod := SCode.NAMEMOD("derivative", mod); end toSubMod; - function getOrder - "returns the order of the given function derivative" - input FunctionDerivative funcDer; - output Integer order; + function perfectFit + "checks if the derivative is a perfect fit for specified interface map" + input FunctionDerivative fnDer; + input UnorderedMap interface_map; + output Boolean b = true; + protected + String name; + Condition cond; algorithm - order := match funcDer.order - case Expression.INTEGER(value = order) then order; - else algorithm - Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed because the order was not evaluated to be a constant: " + Expression.toString(funcDer.order)}); - then fail(); - end match; - end getOrder; + for condition in fnDer.conditions loop + (_, name, cond) := condition; + // if a zero derivative is required but the argument is not in the map + // this function derivative cannot be used + if cond == Condition.ZERO_DERIVATIVE and not UnorderedMap.contains(name, interface_map) then + b := false; return; + end if; + end for; + // if the function derivative is a perfect fit, add all conditions to the interface + for condition in fnDer.conditions loop + (_, name, cond) := condition; + UnorderedMap.add(name, true, interface_map); + end for; + end perfectFit; protected diff --git a/OMCompiler/Compiler/NFFrontEnd/NFInstNode.mo b/OMCompiler/Compiler/NFFrontEnd/NFInstNode.mo index 7950293c0e3..bb08e250d88 100644 --- a/OMCompiler/Compiler/NFFrontEnd/NFInstNode.mo +++ b/OMCompiler/Compiler/NFFrontEnd/NFInstNode.mo @@ -533,14 +533,14 @@ uniontype InstNode output String name; algorithm name := match node - case CLASS_NODE() then "class"; - case COMPONENT_NODE() then "component"; - case INNER_OUTER_NODE() then typeName(node.innerNode); - case REF_NODE() then "ref node"; - case NAME_NODE() then "name node"; - case IMPLICIT_SCOPE() then "implicit scope"; - case EMPTY_NODE() then "empty node"; - case VAR_NODE() then "var node"; + case CLASS_NODE() then "class"; + case COMPONENT_NODE() then "component"; + case INNER_OUTER_NODE() then typeName(node.innerNode); + case REF_NODE() then "ref node"; + case NAME_NODE() then "name node"; + case IMPLICIT_SCOPE() then "implicit scope"; + case EMPTY_NODE() then "empty node"; + case VAR_NODE() then "var node"; end match; end typeName; @@ -917,6 +917,22 @@ uniontype InstNode end match; end setDefinition; + function setComponentDirection + "creates new component!" + input Prefixes.Direction direction; + input output InstNode node; + algorithm + node := match node + case COMPONENT_NODE() algorithm + node.component := Pointer.create(Component.setDirection(Pointer.access(node.component), direction)); + then node; + + else algorithm + Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed for non component node: " + toString(node)}); + then fail(); + end match; + end setComponentDirection; + function info input InstNode node; output SourceInfo info; diff --git a/testsuite/simulation/modelica/NBackend/functions/Makefile b/testsuite/simulation/modelica/NBackend/functions/Makefile index 07c88a15e3a..469085a15cc 100644 --- a/testsuite/simulation/modelica/NBackend/functions/Makefile +++ b/testsuite/simulation/modelica/NBackend/functions/Makefile @@ -1,8 +1,9 @@ TEST = ../../../../rtest -v TESTFILES = \ -function_annotation_der.mos\ builtin_functions.mos\ +function_annotation_der.mos\ +function_diff.mos\ # test that currently fail. Move up when fixed. # Run make failingtest diff --git a/testsuite/simulation/modelica/NBackend/functions/function_annotation_der.mos b/testsuite/simulation/modelica/NBackend/functions/function_annotation_der.mos index 04245609927..44265429ecd 100644 --- a/testsuite/simulation/modelica/NBackend/functions/function_annotation_der.mos +++ b/testsuite/simulation/modelica/NBackend/functions/function_annotation_der.mos @@ -3,94 +3,86 @@ // status: correct loadString(" -model annotate_test - Real x; - Real y; +model function_annotation_der + Real a; + Real b(start=0.0, fixed=true); + Real c(start=0.0, fixed=true); + discrete Real k(start=1.0, fixed=true); equation - y = sin(x); - der(x) = f(y); -end annotate_test; + when a > 0 then + k = pre(k) + 1; + end when; + a = sin(b); + der(b) = f(a, k); + der(c) = f(a, b); +end function_annotation_der; function f input Real x; + input Real k; output Real y; algorithm - y := x^2; - annotation(derivative=df,Inline=false); + y := k*x^2; + annotation(derivative(zeroDerivative=k)=df,Inline=false); end f; function df input Real x; + input Real k; input Real der_x; output Real der_y; algorithm - der_y := 2*x*der_x; + der_y := k*2*x*der_x; annotation(Inline=false); end df; "); getErrorString(); -setCommandLineOptions("--newBackend -d=symjacdump"); getErrorString(); -simulate(annotate_test); getErrorString(); +setCommandLineOptions("--newBackend -d=debugDifferentiation"); getErrorString(); +simulate(function_annotation_der); getErrorString(); // Result: // true // "" // true // "" -// ################################################################################ +// ### debugDifferentiation | NBJacobian.jacobianSymbolic ### +// [BEFORE] [SCAL] (1) a = sin(b) ($RES_SIM_2) +// [AFTER ] [SCAL] (1) $pDER_ODE_JAC.a = cos(b) * $SEED_ODE_JAC.b ($RES_SIM_2) // -// [symjacdump] Creating symbolic Jacobians: +// ### debugDifferentiation | NBJacobian.jacobianSymbolic ### +// [BEFORE] [SCAL] (1) $DER.b = f(a, k) ($RES_SIM_1) +// [AFTER ] [SCAL] (1) $pDER_ODE_JAC.$DER.b = df(a, k, $pDER_ODE_JAC.a) ($RES_SIM_1) // -// ################################################################################ +// ### debugDifferentiation | NBJacobian.jacobianSymbolic ### +// [BEFORE] [SCAL] (1) $DER.c = f(a, b) ($RES_SIM_0) // -// ======================================== -// CONTINUOUS 1 ODE System -// ======================================== +// [BEFORE] function 'f' +// input Real 'x'; +// input Real 'k'; +// output Real 'y'; +// algorithm +// 'y' := 'k' * 'x' ^ 2.0; +// annotation(derivative(order = 1, zeroDerivative = 'k') = 'df', derivative(order = 1) = '$fDER.f', Inline = false); +// end 'f' // -// --- Alias of INI[1 | 2] --- -// BLOCK 1: Single Equation (status = Solve.EXPLICIT) -// ---------------------------------------- -// ### Variable: -// Real y -// ### Equation: -// [SCAL] (1) y = sin(x) ($RES_SIM_1) +// [AFTER] function '$fDER.f' +// input Real 'x'; +// input Real 'k'; +// input Real '$fDER_x'; +// input Real '$fDER_k'; +// output Real '$fDER_y'; +// protected +// Real 'y'; +// algorithm +// '$fDER_y' := k * (2.0 * x * $fDER_x) + $fDER_k * x ^ 2.0; +// 'y' := 'k' * 'x' ^ 2.0; +// annotation(Inline = false); +// end '$fDER.f' // -// --- Alias of INI[1 | 3] --- -// BLOCK 2: Single Equation (status = Solve.EXPLICIT) -// ---------------------------------------- -// ### Variable: -// Real $DER.x -// ### Equation: -// [SCAL] (1) $DER.x = f(y) ($RES_SIM_0) -// -// ======================================== -// [SIM] Jacobian ODE_JAC: -// ======================================== -// -// Result Equations (2/2) -// **************************************** -// (1) [SCAL] (1) $pDER_ODE_JAC.y = cos(x) * $SEED_ODE_JAC.x ($RES_ODE_JAC_0) -// (2) [SCAL] (1) $pDER_ODE_JAC.$DER.x = df(y, $pDER_ODE_JAC.y) ($RES_ODE_JAC_1) -// -// -// ======================================== -// Sparsity Pattern (nnz: 0) -// ======================================== -// -// ### Columns ### -// ---------------------------------------- -// (x) affects: {} -// -// ##### Rows ##### -// ---------------------------------------- -// ($DER.x) depends on: {} -// -// Sparsity Coloring Groups -// ---------------------------------------- -// {x} +// [AFTER ] [SCAL] (1) $pDER_ODE_JAC.$DER.c = $fDER.f(a, b, $pDER_ODE_JAC.a, $SEED_ODE_JAC.b) ($RES_SIM_0) // // record SimulationResult -// resultFile = "annotate_test_res.mat", -// simulationOptions = "startTime = 0.0, stopTime = 1.0, numberOfIntervals = 500, tolerance = 1e-06, method = 'dassl', fileNamePrefix = 'annotate_test', options = '', outputFormat = 'mat', variableFilter = '.*', cflags = '', simflags = ''", +// resultFile = "function_annotation_der_res.mat", +// simulationOptions = "startTime = 0.0, stopTime = 1.0, numberOfIntervals = 500, tolerance = 1e-06, method = 'dassl', fileNamePrefix = 'function_annotation_der', options = '', outputFormat = 'mat', variableFilter = '.*', cflags = '', simflags = ''", // messages = "LOG_SUCCESS | info | The initialization finished successfully without homotopy method. // LOG_SUCCESS | info | The simulation finished successfully. // " diff --git a/testsuite/simulation/modelica/NBackend/functions/function_diff.mos b/testsuite/simulation/modelica/NBackend/functions/function_diff.mos new file mode 100644 index 00000000000..50abdc15678 --- /dev/null +++ b/testsuite/simulation/modelica/NBackend/functions/function_diff.mos @@ -0,0 +1,86 @@ +// name: function_diff +// keywords: NewBackend +// status: correct + +loadString(" +model function_diff + Real x; + Real y; + Real z; +equation + y = sin(x); + der(x) = f(y); + der(z) = f(x); +end function_diff; + +function f + input Real a; + output Real b; +algorithm + if a >= 0 then + b := a^2; + else + b := f(-a); + end if; + annotation(Inline=false); +end f; +"); getErrorString(); + +setCommandLineOptions("--newBackend -d=debugDifferentiation"); getErrorString(); +simulate(function_diff); getErrorString(); +// Result: +// true +// "" +// true +// "" +// ### debugDifferentiation | NBJacobian.jacobianSymbolic ### +// [BEFORE] [SCAL] (1) y = sin(x) ($RES_SIM_2) +// [AFTER ] [SCAL] (1) $pDER_ODE_JAC.y = cos(x) * $SEED_ODE_JAC.x ($RES_SIM_2) +// +// ### debugDifferentiation | NBJacobian.jacobianSymbolic ### +// [BEFORE] [SCAL] (1) $DER.x = f(y) ($RES_SIM_1) +// +// [BEFORE] function 'f' +// input Real 'a'; +// output Real 'b'; +// algorithm +// if 'a' >= 0.0 then +// 'b' := 'a' ^ 2.0; +// elseif true then +// 'b' := 'f'(-'a'); +// end if; +// annotation(derivative(order = 1) = '$fDER.f', Inline = false); +// end 'f' +// +// [AFTER] function '$fDER.f' +// input Real 'a'; +// input Real '$fDER_a'; +// output Real '$fDER_b'; +// protected +// Real 'b'; +// algorithm +// if 'a' >= 0.0 then +// '$fDER_b' := 2.0 * a * $fDER_a; +// 'b' := 'a' ^ 2.0; +// elseif true then +// '$fDER_b' := '$fDER.f'(-'a', -'$fDER_a'); +// 'b' := 'f'(-'a'); +// end if; +// annotation(Inline = false); +// end '$fDER.f' +// +// [AFTER ] [SCAL] (1) $pDER_ODE_JAC.$DER.x = $fDER.f(y, $pDER_ODE_JAC.y) ($RES_SIM_1) +// +// ### debugDifferentiation | NBJacobian.jacobianSymbolic ### +// [BEFORE] [SCAL] (1) $DER.z = f(x) ($RES_SIM_0) +// [AFTER ] [SCAL] (1) $pDER_ODE_JAC.$DER.z = $fDER.f(x, $SEED_ODE_JAC.x) ($RES_SIM_0) +// +// record SimulationResult +// resultFile = "function_diff_res.mat", +// simulationOptions = "startTime = 0.0, stopTime = 1.0, numberOfIntervals = 500, tolerance = 1e-06, method = 'dassl', fileNamePrefix = 'function_diff', options = '', outputFormat = 'mat', variableFilter = '.*', cflags = '', simflags = ''", +// messages = "LOG_SUCCESS | info | The initialization finished successfully without homotopy method. +// LOG_SUCCESS | info | The simulation finished successfully. +// " +// end SimulationResult; +// "" +// endResult