Your current code has fragments like this:
weighted_load = load_obj(f'./sym.sobj')
print(f'weighted == weighted_load {weighted == weighted_load}')
save(weighted, f'./sym.sobj')
In particular, "sym.sobj" is read before it is written, so 'weighted_load'
will be None. If I put the "save" lines earlier, I see this:
weighted == weighted_load True
realm == realm_load: True
realm != realm_load: False
determinant
-7.116356115499656e-14
-7.116356115499656e-14
-7.116356115499656e-14
Having said all of that, I think that determinants tend to be numerically
unstable, and it would not be surprising for a number so close to zero to
display that instability. A small change in one of the entries of the
matrix, due to how the decimals are represented, could easily change the
determinant.
On Wednesday, January 4, 2023 at 9:12:55 AM UTC-8 ivana...@gmail.com wrote:
> I found I got different determinant value every time I ran my program. I
> tried to save the current matrix and load the matrix in the next run, and I
> compared two matrices. The compare results are not stable either. Sometime
> they are the same and not the same at the same time (A == B is False and A
> != B is False either)!
>
> Here is the minimal case I could figure out. I tried to replace the
> testing matrix to a simpler one, but the issue could not be reproduced
> after replacing. I am sorry about giving you some many lines of code. I
> guess it's not easy to read code on Google Groups, so the file is attached.
>
> from functools import cache
>
> from sage.all_cmdline import *
>
>
> l = 1
>
>
> def load_obj(filename):
> try:
> return load(filename)
> except FileNotFoundError:
> return None
>
>
> @cache
> def expected_distance(n, l):
> if n <= -2:
> return 0
> if n == -1:
> return -2 * x
> if n == 0:
> return 1
> return ((
> - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
> - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n
> - 2, l)
> ) / (8 * (n + 1) * x)).full_simplify()
>
> # Generate the matrix for testing. I tried to figure out a minimal
> reproducible
> # code, but I failed.
> depth = 8
> functions = [expected_distance(i, l) for i in range(depth * 2 - 1)]
> hankel = matrix.hankel(functions[:depth], functions[depth:depth * 2 - 1],
> SR)
> weighted = matrix(depth, depth, [(hankel[i][j] / hankel[0][j] / hankel[i][
> 0])
> for i in range(depth)
> for j in range(depth)])
>
> # Only with this useless line, you can see the first and the second cases
> # listed in the comment (*) below.
> weighted = weighted.submatrix(0, 0, depth, depth)
>
> weighted_load = load_obj(f'./sym.sobj')
> print(f'weighted == weighted_load {weighted == weighted_load}')
> save(weighted, f'./sym.sobj')
>
> realm = weighted.substitute({x: -0.0070218})
> realm2 = weighted.substitute({x: -0.0070218})
> realm_load = load_obj(f'./real.sobj')
> save(realm, f'./real.sobj')
>
> # (*) There are three possible output of this two lines
> # 1.
> # realm == realm_load: False
> # realm != realm_load: False
> # 2.
> # realm == realm_load: True
> # realm != realm_load: False
> # 3.
> # realm == realm_load: False
> # realm != realm_load: True
> #
> # The possibilities of each case are almost the same.
> print(f'realm == realm_load: {realm == realm_load}')
> print(f'realm != realm_load: {realm != realm_load}')
>
> # The determinant values are different in defferent runs
> print('determinant')
> print(realm.det())
> print(realm2.det())
> print(realm_load.det())
>
> # Only print something when the third case happened.
> for i in range(realm.nrows()):
> for j in range(realm.ncols()):
> if realm[i][j] != realm_load[i][j]:
> print(f'{(i, j)}')
> print(realm[i][j])
> print(realm_load[i][j])
>
>
--
You received this message because you are subscribed to the Google Groups
"sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sage-devel+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sage-devel/a7e473a2-7dfa-4e69-9f15-651ef6e77c09n%40googlegroups.com.
