Derivative of conditional
Consider
element = FiniteElement("Lagrange", interval, 1)
u = Coefficient(element)
du = TrialFunction(element)
v = TestFunction(element)
f = conditional(le(abs(u), 0.00015478), u, sin(u)/u)
fd = diff(f, u)
L = v * f * dx
A = inner(v, dot(fd, du)) * dx
# A = derivative(L, u, du)
forms = [L, A]
where the sin(u)/u is replaced by its Taylor expansion for small u in order to avoid the singularity for u=0. The code for derivative will expand this into
dc = conditional(c, 1, 0)
return dc*dt + (1.0 - dc)*df
This is problematic because both expressions are going to be computed in the code generated by ffc, and summed together, one multiplied by 1, the other by 0. The final result is a nan derivative for u=0.
Comments (18)
-
-
reporter @martinal , apologize for asking, but I'm still fighting this - and loosing. Is there an easy way to transform an expression like
f_4 * (v_1 * f_4 + v_1 * f_4) + v_1 * f_4 * f_4 + v_1
into something like
(f_4 * (1 * f_4 + 1 * f_4) + 1 * f_4 * f_4 + 1) * v_1
?
If yes, do you think it would be ok, at least for ffc (I've not looked into uflacs yet), to have
def conditional(self, o, unused_dc, dt, df): if isinstance(dt, Zero) and isinstance(df, Zero): # Assuming dt and df have the same indices here, which should be the case return dt else: # collect arguments c, t, f = o.ufl_operands (dtt, at) = the_sought_transformation(dt) (dft, af) = the_sought_transformation(df) #check that at == af, otherwise raise error dc = conditional(c, dtt, dft) * at return dc
where
the_sought_transformation(dt)
should returnf_4 * (1 * f_4 + 1 * f_4) + 1 * f_4 * f_4 + 1
asdtt
andv_1
asat
?Please bear with me if I'm still misunderstanding how ufl/ffc works.
Thanks
-
1) the_sought_transformation is only easy to do when every subexpression is scalar.
2) uflacs does something like that at an intermediate stage (no, you can't use that at the ufl level)
3) the_sought_transformation is not sufficient, it wouldn't deal with for example
grad(v)
4) the case that the
conditional(c, 1, 0)
fix was added for will still break with your code here.I suggest that you create a branch of ufl and change
conditional
to returnconditional(c, dt, df)
, which I believe should work in your case, and hope that you don't get problems with the case that doesn't work for. -
reporter No, returning
conditional(c, dt, df)
does not work. After adding a suitablecheck_arity
for conditionals and removing the checksffc_assert(len(true) == 1 and firstkey(true) == (),\ "True value of Condtional should only be one function: " + repr(true)) ffc_assert(len(false) == 1 and firstkey(false) == (),\ "False value of Condtional should only be one function: " + repr(false))
in ffc quadraturetransformer.py I end with this kind of generated code:
C[0] = ((F0 <= 0.7645)) ? FE0[ip][k]*0.123 : FE0[ip][k]; // Number of operations for primary indices: 16 for (unsigned int j = 0; j < 2; j++) { for (unsigned int k = 0; k < 2; k++) { // Number of operations to compute entry: 4 A[j*2 + k] += FE0[ip][j]*C[0]*W2[ip]*det; } // end loop over 'k' } // end loop over 'j'
because conditionals are special-handled, and evaluated early, in ffc quadraturegenerator.py. And I fear that moving them inside the loops (if possible at all) would lead to inefficient code. That's why I asked.
-
I think I have a fix in uflacs. Can you test it with a real example?
Install the latest next branch of ufl and uflacs, then set this at the top of your program:
from dolfin import * parameters["form_compiler"]["representation"] = "uflacs"
-
- changed milestone to 1.7
-
assigned issue to
- changed version to dev
-
reporter Thanks! The "real" example fails with
... File "/home/marco/local/Fenics/lib/python2.7/site-packages/uflacs/analysis/graph_rebuild.py", line 267, in rebuild_with_scalar_subexpressions ffc_assert(len(vs) == len(ws), "Expecting one symbol for each expression.") File "/home/marco/local/Fenics/lib64/python2.7/site-packages/ffc/log.py", line 44, in ffc_assert condition or error(*message) File "<string>", line 1, in <lambda> File "/home/marco/local/Fenics/lib/python2.7/site-packages/ufl/log.py", line 158, in error raise self._exception_type(self._format_raw(*message)) Exception: Expecting one symbol for each expression.
but I'm not 100% sure the test is correct: up to now I've never been able to run it! ffc with CONDITIONAL_WORKAROUND = False can however build it, although incorrectly.
I can send the test to you or attach it somewhere (three files, I'm dealing with finite rotations) should you think it could be useful. Do you need the complete backtrace and/or the real test?
I think this simple test
from dolfin import * parameters["form_compiler"]["representation"] = "uflacs" def coefa(phi): return conditional(le(phi, 0.7645), phi * phi * phi * 0.123, 0.987 * phi * phi * phi) mesh = UnitIntervalMesh(1) U_SPACE = FunctionSpace(mesh, "CG", 1) u = Function(U_SPACE) v = TestFunction(U_SPACE) ut = TrialFunction(U_SPACE) functional = inner(v, (backrot_u)) * dx f = assemble(functional) print as_backend_type(f).vec().view() J = assemble(derivative(functional, u, ut)) print as_backend_type(J).mat().view()
is still compiled incorrectly. Inside tabulate_tensor (for J) I see
sv3[0] = weights3[iq] * sp3[0]; sv3[1] = w0 <= 0.7645; sv3[2] = 0.987 * w0; sv3[3] = sv3[2] * w0; sv3[4] = sv3[2] + sv3[2]; sv3[5] = sv3[4] * w0; sv3[6] = sv3[5] + sv3[3]; sv3[7] = (sv3[1] ? 2 : sv3[6]);
where the ``true``` value seems to be wrong. In the meantime I'm trying to fix ffc; I hope to be on the right track.
-
If you can send me the real test that's very useful.
Also the simple test doesn't run: backrot_u is not defined, and it doesn't actually use the conditional. If you can post a working simple test, that's also very useful.
-
reporter - attached RotCoef.py
- attached beam.py
- attached RotVector.py
- attached beamtest.py
Apologize, I erased one extra line after pasting the simple test. I'm attaching the files this time.
The simple test, that concerns me more, is beamtest.py.
The complex testcase is beam.y, and requires RotCoef.py and RotVector.py. This is really work in progress: I would avoid wasting your time for some silly error on my part (of course, we can develop it together should you be interested).
In beam.py you'll find there two versions: version 1 builds, version 2 does not build. I think they should give the same result, but I may well be wrong.
P.S. shouldn't you change also the implementation of conditional in forward_ad.py ?
-
Considering the simple case first. It compiles just fine here with uflacs. The code looks like a perfect replicate of the input. And the result is of course zero when you assemble something*u and u=0.
(P.S. that code is old and unused and should instead be removed).
-
reporter - attached beamtest2.py
This is the output I get with beamtest.py
marco@pao:~/Projects/Chierichetti/AnbaDolfin/FiniteRotations/T> python beamtest.py Vec Object:() 1 MPI processes type: mpi Process [0] 0 0 None Calling FFC just-in-time (JIT) compiler, this may take some time. Mat Object:() 1 MPI processes type: seqaij row 0: (0, 0.666667) (1, 0.333333) row 1: (0, 0.333333) (1, 0.666667) None
The generated code for the Jacobian matrix is
double sv3[9]; sv3[0] = weights3[iq] * sp3[0]; sv3[1] = w0 <= 0.7645; sv3[2] = 0.987 * w0; sv3[3] = sv3[2] * w0; sv3[4] = sv3[2] + sv3[2]; sv3[5] = sv3[4] * w0; sv3[6] = sv3[5] + sv3[3]; sv3[7] = (sv3[1] ? 2 : sv3[6]); sv3[8] = sv3[0] * sv3[7]; for (int ia0 = 0; ia0 < 2; ++ia0) { for (int ia1 = 0; ia1 < 2; ++ia1) { A[2 * ia0 + ia1] += sv3[8] * FE1_C0_Q3[0][iq][ia0 - 0] * FE1_C0_Q3[0][iq][ia1 - 0]; } }
I'm not sure it is correct. I've double-checked removing and reinstalling ufl. ffc, and uflacs, and clearing $HOME/.instant .
Do you see a different output?
For the newly attached beamtest2.py I get
marco@pao:~/Projects/Chierichetti/AnbaDolfin/FiniteRotations/T> python beamtest2.py Vec Object:() 1 MPI processes type: mpi Process [0] 0 0 None Mat Object:() 1 MPI processes type: seqaij row 0: (0, 0) (1, 0) row 1: (0, 0) (1, 0) None
and the code generated for the Jacobian matrix is
double sv3[14]; sv3[0] = weights3[iq] * sp3[0]; sv3[1] = w0 <= 0.7645; sv3[2] = w0 + w0; sv3[3] = sv3[2] * w0; sv3[4] = w0 * w0; sv3[5] = sv3[3] + sv3[4]; sv3[6] = 0.123 * sv3[5]; sv3[7] = 0.987 * w0; sv3[8] = sv3[7] + sv3[7]; sv3[9] = sv3[8] * w0; sv3[10] = sv3[7] * w0; sv3[11] = sv3[9] + sv3[10]; sv3[12] = (sv3[1] ? sv3[6] : sv3[11]); sv3[13] = sv3[0] * sv3[12]; for (int ia0 = 0; ia0 < 2; ++ia0) { for (int ia1 = 0; ia1 < 2; ++ia1) { A[2 * ia0 + ia1] += sv3[13] * FE1_C0_Q3[0][iq][ia0 - 0] * FE1_C0_Q3[0][iq][ia1 - 0]; } }
that, apart for some missed optimizations, I think is correct.
For the FFC side, this https://bitbucket.org/MarcoMorandini/marcomorandini-ffc/branch/marcomorandini%2Ffix-ufl-issue70 is what I was able to do up to now. Does it make sense?
The file
optimisedquadraturetransformer.py
is still failing (withparameters["form_compiler"]["optimize"] = True
):File "/home/marco/local/Fenics/lib64/python2.7/site-packages/ffc/representation.py", line 229, in _compute_integral_ir parameters) File "/home/marco/local/Fenics/lib64/python2.7/site-packages/ffc/quadrature/quadraturerepresentation.py", line 93, in compute_integral_ir itg_data.integral_type, cell) File "/home/marco/local/Fenics/lib64/python2.7/site-packages/ffc/quadrature/quadraturerepresentation.py", line 159, in _transform_integrals_by_type terms = _transform_integrals(transformer, integrals_dict, integral_type) File "/home/marco/local/Fenics/lib64/python2.7/site-packages/ffc/quadrature/quadraturerepresentation.py", line 202, in _transform_integrals terms = transformer.generate_terms(integral.integrand(), integral_type) File "/home/marco/local/Fenics/lib64/python2.7/site-packages/ffc/quadrature/quadraturetransformerbase.py", line 770, in generate_terms value, ops, sets = self._create_entry_data(val, integral_type) File "/home/marco/local/Fenics/lib64/python2.7/site-packages/ffc/quadrature/optimisedquadraturetransformer.py", line 659, in _create_entry_data for b in value.get_unique_vars(BASIS)]) KeyError: Symbol('(bC[0] ? 0.123*(F0*(F0*FE0[ip][k] + F0*FE0[ip][k]) + F0*F0*FE0[ip][k]) : (0.987*F0*F0*FE0[ip][k] + F0*(0.987*F0*FE0[ip][k] + 0.987*F0*FE0[ip][k])))', BASIS)
Is this because I've introduced a new boolean type? At any rate, if you think the approach is sound I will try to finish the work.
Thanks
-
I did something wrong, I see the same output now. I need to revisit the fix in uflacs.
I'm not the author of that code in ffc and my intention is to replace it with uflacs so don't ask me.
-
There's a new fix in the martinal/fix-conditionals branch of uflacs. It needs more detailed testing before I trust this code, but give it a try will you?
NB! I've reset next to master because this wasn't well enough tested.
-
reporter - attached beam2.py
With martinal/fix-conditionals the "version 1" form of beam.py builds and runs fine. A simple test for a beam with a distributed transverse load gives reasonable results. I get the same results with my ffc branch. While this is not a correctness proof, it is at least reassuring. Uflacs generated code is much faster than ffc.
However, by complicating a little the form (beam2.py) I get an error. It is the same error message I get with the commented "version 2" of beam.py. It's not clear to me whether this is an error on my side or a different bug.
-
reporter - attached test_conditional.py
File
test_conditional.py
is another reduced testcase.I think that uflacs gets confused in
graph_rebuild.py
becausews
is a conditional returning a vector.The workaround I'm using right now is simple (i.e.
UnWrap2
instead ofUnWrap
in the test), and I could live with it. With the workaround I get almost correct results for a difficult beam test case.I'm still investigating whether the remaining issues are due to the handling of conditionals or to an error on my side.
-
Thanks for the test, I'll take a look. I think you're right I didn't consider vector valued conditionals.
-
- changed status to resolved
-
- removed milestone
Removing milestone: 1.7 (automated comment)
- Log in to comment
This is known. Just documenting a possible solution I just thought of (possibly incomprehensible for others):
Let d/du conditional(c, t, f) = conditional(c, d/du t, d/du f) again, like it was before (the current situation is a fix for a form compiler issue)
In uflacs, let factorization phase handle cases like these that currently confuse it:
conditional(c, f*u, 0) -> conditional(c, f, 0)*u
conditional(c, f*u, g*u) -> conditional(c, f, g)*u
conditional(c, f*u, g*u.dx(0)) -> conditional(c, f, 0)*u + conditional(c, 0, g)*u.dx(0)
However if I do the changes in ufl and uflacs, standard ffc will likely break in the case
conditional(c, 1, 0)
was fixing initially.I may not have time to do this anytime soon.