-4

This piece of code:

K = 3
N = 3
E = [np.reshape(np.array(i), (K, N)) for i in itertools.product([0, 1, -1], repeat = K*N)]

print 'E = ', E

Generates all possible E matrices (dimensions 3x3) formed by 2 integers: 0 and 2, e.g.:

...
array([[0, 2, 2],
       [0, 0, 0],
       [2, 0, 0]]), array([[0, 2, 2],
       [0, 0, 0],
       [2, 0, 2]]), array([[0, 2, 2],
       [0, 0, 0],
       [2, 2, 0]])
...

Given this matrix equation:

A_SC = E * A     # Eqn. 1

where:

1) * stands for the standard matrix multiplication (rows, columns)

2) A_SC, E and A are 3x3 matrices,

3) E are all the possible integer matrices generated by the above code.

4) A is a known matrix:

A =np.array([[   0.288155519353E+01,   0.000000000000E+00,   0.568733333333E+01],
             [  -0.144077759676E+01,   0.249550000000E+01,   0.568733333333E+01],
             [  -0.144077759676E+01,  -0.249550000000E+01,   0.568733333333E+01]])

The A_SC matrix can be represented as 3 rows vectors: a1_SC, a2_SC and a3_SC:

       |a1_SC|
A_SC = |a2_SC|
       |a3_SC|

For a given E matrix, there is a A_SC matrix.

The following code:

1) loops over all possible E matrices,

2) calculates the A_SC matrix,

3) calculates the norm of a1_SC, a2_SC and a3_SC,

4) and calculates the determinant of the E matrix in that iteration:

for indx_E in E:
      A_SC = np.dot(indx_E,A)
      a1_SC = np.linalg.norm(A_SC[0])
      a2_SC = np.linalg.norm(A_SC[1])
      a3_SC = np.linalg.norm(A_SC[2])

      det_indx_E = np.linalg.det(indx_E)
      print 'a1_SC = ', a1_SC
      print 'a2_SC = ', a2_SC
      print 'a3_SC = ', a3_SC
      print 'det_indx_E = ', det_indx_E

The goal would be to obtain all those A_SC and E matrices (Eqn. 1) for which the norm of these 3 rows vectors is the same and greater than 10,

norm(a1_SC) = norm(a2_SC) = norm(a3_SC) > 10

And at the same time, the determinant of E has to be greater than 0.0. This condition can be expressed in this way: Just after this for loop, we can write an if loop:

tol_1 = 10
tol_2 = 0

for indx_E in E:

      A_SC = np.dot(indx_E,A)
      a1_SC = np.linalg.norm(A_SC[0])
      a2_SC = np.linalg.norm(A_SC[1])
      a3_SC = np.linalg.norm(A_SC[2])

      det_indx_E = np.linalg.det(indx_E)
      print 'a1_SC = ', a1_SC
      print 'a2_SC = ', a2_SC
      print 'a3_SC = ', a3_SC
      print 'det_indx_E = ', det_indx_E

      if  a1_SC > tol_1\
          and a2_SC > tol_1\
          and a3_SC > tol_1\
          and abs(a1_SC - a2_SC) == tol_2\
          and abs(a1_SC - a3_SC) == tol_2\
          and abs(a2_SC - a3_SC) == tol_2\
          and det_indx_E > 0.0:
             print 'A_SC = ', A_SC

             print 'a1_SC = ', a1_SC
             print 'a2_SC = ', a2_SC
             print 'a3_SC = ', a3_SC
             print 'det_indx_E = ', det_indx_E

             # Now, which is the `E` matrix for this `A_SC` ?
             #      A_SC = E * A     # Eqn. 1
             #      A_SC * inv(A) = E * A * inv(A)  # Eqn. 2
             #
             #      ------------------------------
             #     | A_SC * inv(A) = E  # Eqn. 3  |
             #      ------------------------------
             E_sol = np.dot(A_SC, np.linalg.inv(A))
             print 'E_sol = ', E_sol

Just to be clear, this is the entire code:

A =np.array([[   0.288155519353E+01,   0.000000000000E+00,   0.568733333333E+01],
                 [  -0.144077759676E+01,   0.249550000000E+01,   0.568733333333E+01],
                 [  -0.144077759676E+01,  -0.249550000000E+01,   0.568733333333E+01]])

K = 3
N = 3
E = [np.reshape(np.array(i), (K, N)) for i in itertools.product([0, 1, -1], repeat = K*N)]

print 'type(E) = ', type(E)
print 'E = ', E
print 'len(E) = ', len(E)

tol_1 = 10
tol_2 = 0

for indx_E in E:

      A_SC = np.dot(indx_E,A)
      a1_SC = np.linalg.norm(A_SC[0])
      a2_SC = np.linalg.norm(A_SC[1])
      a3_SC = np.linalg.norm(A_SC[2])

      det_indx_E = np.linalg.det(indx_E)
      print 'a1_SC = ', a1_SC
      print 'a2_SC = ', a2_SC
      print 'a3_SC = ', a3_SC
      print 'det_indx_E = ', det_indx_E

      if  a1_SC > tol_1\
          and a2_SC > tol_1\
          and a3_SC > tol_1\
          and abs(a1_SC - a2_SC) == tol_2\
          and abs(a1_SC - a3_SC) == tol_2\
          and abs(a2_SC - a3_SC) == tol_2\
          and det_indx_E > 0.0:
             print 'A_SC = ', A_SC

             print 'a1_SC = ', a1_SC
             print 'a2_SC = ', a2_SC
             print 'a3_SC = ', a3_SC
             print 'det_indx_E = ', det_indx_E

             # Now, which is the `E` matrix for this `A_SC` ?
             #      A_SC = E * A     # Eqn. 1
             #      A_SC * inv(A) = E * A * inv(A)  # Eqn. 2
             #
             #      ------------------------------
             #     | A_SC * inv(A) = E  # Eqn. 3  |
             #      ------------------------------
             E_sol = np.dot(A_SC, np.linalg.inv(A))
             print 'E_sol = ', E_sol

The problem is that no A_SC (and therefore no E_sol) are printed. If you run this code, all norms and determinants are printed at each iteration, for example the following:

a1_SC =  12.7513326014
a2_SC =  12.7513326014
a3_SC =  12.7513326014
det_indx_E =  8.0

This would be a perfect candidate, because it satisfies

a1_SC = a2_SC = a3_SC = 12.7513326014 > 10.0

and

determinant > 0.0

However, no A_SC (and therefore no E_sol) are printed... why is this happening?

For example, this E matrix:

        2 0 0
  E =   0 2 0
        0 0 2

has det = 8.0, and is a candidate, because it has:

a1_SC = a2_SC = a3_SC = 12.7513326014 > 10.0
DavidC.
  • 669
  • 8
  • 26
  • @MohitC When I say `Just to be clear, this is the entire code:` there I have pasted the entire code, as a minimal working example – DavidC. Oct 17 '17 at 09:56
  • 1
    If working with floats, it is not guaranteed that if `(str(a) == str(b))` (the rendered values look equal), then `(a - b) == 0`. Maybe you should try to set `tol_2` to something like 1e-6 and test if `abs(a1_SC - a2_SC) <= tol_2` – Ronald Oct 17 '17 at 09:59
  • 2
    It's not at all minimal, but your problem is using `==` on floats. If you had bothered to error trace you'd see all your `abs(. . . . ) == tol2` checks are likely `False` because of floating point errors. Use `np.isclose` – Daniel F Oct 17 '17 at 10:00
  • @Ronald Thank you very much for pointing that out. It is true that the rendered values look equal, and I agree with you that this does not mean that it is going to be true that `abs(a1_SC - a2_SC) == tol_2`. I have set `tol_2` to `1E-6` and performed `abs(a1_SC - a2_SC) <= tol_2` instead of `abs(a1_SC - a2_SC) == tol_2`. This is prints out all `A_SC` matrices for which the 7 conditions in the `if` loop are satisfied. Therefore, this solves the problem. Again, thank you so much. I quite do not understand why this post deserves 4 negative votes, though. – DavidC. Oct 17 '17 at 17:09
  • @DanielF I apologise if the `entire code` seems not to be minimal. Why would you say it is not a MWE? If you copy and paste that `entire code`, you will get exactly all the errors I mention in the post. If this post seems long is because I wanted to explain in detail the problem – DavidC. Oct 17 '17 at 17:09
  • @DanielF Anyway, thank you so much for pointing out the idea of checking `False` or `True` in the `abs` conditions. I have tried `print abs(a1_SC - a2_SC) == tol_2` before the `if` loop, and all the messages are `False`, so this confirms where was the error, thank you so much. Also, I don't quite understand the tone in "If you had bothered to error trace"... I simply did not have that idea. If this is the reason for 4 negative votes, I personally think it is not fair. I dedicated quite a lot of time to explain the problem in detail. Thank you very much for pointing out where was the problem. – DavidC. Oct 17 '17 at 17:10
  • [mcve] starts with "minimal" for a reason. The "entire code" is never minimal unless you haven't done any work yourself - for which you will be rightfully downvoted. The minimal code to reproduce the problem in this case is: " `abs(a - b) == 0` is always `False` even when `a` and `b` seem to be equal." - which is the question you'd ask if you'd traced your error back through the code and found the point that didn't make sense. – Daniel F Oct 18 '17 at 09:30
  • @DanielF Thanks for this. `The "entire code" is never minimal unless you haven't done any work yourself - for which you will be rightfully downvoted.` I just wanted to say that all that "entire code" was the result of my effort and work done all by my myself. My problem is that perhaps I did not summarized the "entire code" into the key point. Also, I did not have that brilliant idea to trace back the error by checking the `True` or `False` answer by asking for `print abs(a1_SC - a2_SC) == tol_2` before the `if` loop. Anyway, thanks a lot for all your more than useful help :) – DavidC. Oct 18 '17 at 10:09
  • By "work" I meant debugging work, not coding. Programming isn't just about writing code. Basically the philosophy of SO on troubleshooting questions is: if your problem isn't condensed down to a few lines of code and a clear "Here's what I expected this code to do" and "here's what happened," then your question 1) isn't valuable to later readers and 2) implies your time is more valuable than the answerer's - both of which attract downvotes. – Daniel F Oct 19 '17 at 06:10
  • Check out [Rubber Duck Problem Solving](https://blog.codinghorror.com/rubber-duck-problem-solving/) - the process of asking a good question will often yield the answer before you even finish asking it. Take it from someone who's deleted ten times as many half-finished questions then they've hit "submit" on - it works! – Daniel F Oct 19 '17 at 06:12
  • @DanielF Thanks for this – DavidC. Oct 19 '17 at 15:42

1 Answers1

0

The simple answer would be not to mistake double output as string with real underlying precision. Simplest change:

tol_2 = 1e-8

and change the tol_2 related conditions to limitation:

        and abs(a1_SC - a2_SC) <= tol_2\
        and abs(a1_SC - a3_SC) <= tol_2\
        and abs(a2_SC - a3_SC) <= tol_2\

That should solve your issue.

Remember that when you do not explicitly use symbolic computations on a computer, you should always be ready for numerical errors even in the simplest of examples

And if you need to check for STRICT, mathematical equality of some things - you have to use symbolic math packages and related machinery.

In case if the required equality is more of a 'physical' sense (like, what force is enough for pulling that box) - then the approach I described would be OK since some 'error' is always present in a physical world, you just have to state required tolerance (which we do with tol_2 in this case)

E P
  • 116
  • 8
  • Probably more info on tolerances and python floating-point comparisons in [this question](https://stackoverflow.com/questions/5595425/what-is-the-best-way-to-compare-floats-for-almost-equality-in-python) – E P Oct 17 '17 at 10:27