Thanks for the explanation. I moved the doit()s to later in the script
where everything boils down to scalars and that works.
I had tried jacobean but sometimes that gives me AttributeError. Here
is the real thing I was working on:
def lpu_nonlin(Y, X, C):
"""
LPU for non-linear measurement models according to
Giovanni Mana and Francesca Pennecchi
Metrologia 44 (2007) 246-251
Inputs:
Y measurement model e.g. Y = (X[0] / X[1]) - X[2]
X vector of input quantities
C covariance matrix
Returns symbolic equations for:
expectation using LPU with 3rd degree Taylor series
variance using LPU with 3rd degree Taylor series
expectation using LPU with 1st degree Taylor series
variance using LPU with 1st degree Taylor series
"""
#Jacobean vector
J = sympy.Matrix([Y]).jacobian(X).T
#Hessian matrix
H = sympy.hessian(Y, X)
#Components of Mana & Pennacchi equations 2 and 3
JTC = J.T * C
JTCJ = (J.T * C * J)[0, 0]
HC = H * C
TrHC = HC.trace()
HC2 = HC**2
TrHC2 = HC2.trace()
dTrHC = sympy.Matrix([sympy.diff(TrHC, X)]).T
#dTrHC = sympy.Matrix([TrHC]).jacobean(X).T
#Expectation, Mana & Pennacchi eqn 2
exp3_y = (Y + TrHC / 2).doit()
#variance, Mana & Pennacchi eqn 3
var3_y = (JTCJ + TrHC2/2 + (JTC * dTrHC)[0, 0]).doit()
#according to GUM framework LPU
exp1_y = Y
var1_y = JTCJ.doit()
return exp3_y, var3_y, exp1_y, var1_y
For getting J I already changed from diff to jacobean as you suggested
but doing the same a few lines further on with TrHC fails (the
commented version, TrHC is a scalar function of X and C).
Why can that be? And why is the syntax of hessian and jacobean so
different?
Thanks again, Janwillem
On May 4, 11:55 am, Bastian Weber <bastian.we...@gmx-topmail.de>
wrote:
> Hi Jan,
>
>
>
> janwillem wrote:
> > I need some explanation on the workings of SymPy. As an example the
> > following script:
> > import sympy
> > X, F, B = sympy.symbols('XFB')
> > Y = X / F - B #eqn 1
> > DY = sympy.Matrix(sympy.diff(Y, (X, F, B))).T
> > print DY
>
> > I had expected (eqn 2): [1/F, -X/F**2, -1]
> > But get: [D(-B + X/F, X), D(-B + X/F, F), D(-B + X/F, B)]
>
> > Not after doing a trick of which I cannot remember why I tried it, I
> > get the desired result
> > DY = DY.subs({X:X, F:F, B:B})
> >>From the doc I had thought that DY.doin() would work but that gives
> > "raise AttributeError()".
> > So there is obviously something I do not understand, please some help
> > Janwillem
>
> .doit is a method of derivative-objects ( like D(-B + X/F, X) ) and not
> of Matrix-objects (like DY)
>
> What you could do is:
>
> DY.applyfunc(lambda arg: arg.doit())
>
> But maybe thats not very intuitve (at least as long as you are not
> familiar with pythons lambda functions).
>
> An alternative way to calculate your desired result would be
>
> DY = Matrix([Y]).jacobian([X,F,B])
>
> The reason why your subs-trick works is out of my knowledge.
>
> Regards,
> Bastian.
>
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to
> sympy+unsubscr...@googlegroups.com.
> For more options, visit this group
> athttp://groups.google.com/group/sympy?hl=en.
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to sy...@googlegroups.com.
To unsubscribe from this group, send email to
sympy+unsubscr...@googlegroups.com.
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
