Derivative of conditional

[Marco Morandini]

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.

[Martin Sandve Alnæs]

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.

[Marco Morandini]

@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 return f_4 * (1 * f_4 + 1 * f_4) + 1 * f_4 * f_4 + 1 as dtt and v_1 as at ?

Please bear with me if I'm still misunderstanding how ufl/ffc works.

Thanks

[Martin Sandve Alnæs]

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 return conditional(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.

[Marco Morandini]

No, returning conditional(c, dt, df) does not work. After adding a suitable check_arity for conditionals and removing the checks

        ffc_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.

[Martin Sandve Alnæs]

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"

[Marco Morandini]

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.

[Martin Sandve Alnæs]

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.

[Marco Morandini]

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 ?

[Martin Sandve Alnæs]

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).

[Marco Morandini]

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 (with parameters["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

[Martin Sandve Alnæs]

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.

[Martin Sandve Alnæs]

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.

[Marco Morandini]

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.

[Marco Morandini]

File test_conditional.py is another reduced testcase.

I think that uflacs gets confused in graph_rebuild.py because ws is a conditional returning a vector.

The workaround I'm using right now is simple (i.e. UnWrap2 instead of UnWrap 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.

[Martin Sandve Alnæs]

Thanks for the test, I'll take a look. I think you're right I didn't consider vector valued conditionals.

[Martin Sandve Alnæs]

Removing milestone: 1.7 (automated comment)