/author of this script here
I’m not really a very good pythonist, so there may be silly ways to solve simple problems in my script.
You have all described my method correctly. But I think you looked at the wrong file, because for 2d my function is called general_potential
You can also start with a 30 * 30 or 40 * 40 well matrix, this will reduce the computation time, but it will already show some interesting quantum patterns
@weightt_an
Correct. I inadvertently picked up the 3D script. The 2D has 5 diagonals populated, not 7, and the data quantities aren’t as intense. Sorry to be such a putz. Welcome to the board!
Also diagmN and diagpN arrays lays not on ±3 diagonals, but on ±N diagonals. And diagmNN diagpNN on ± N^2 . With 0 … N^2 or 0 … N^3 numeration
OK! Took the tour bus to Schrodinger’s Well and shot some pictures! Thank you, @weightt_an for the ride!
eigs() was throwing a lot of exceptions about being unable to solve the Hamiltonian. I had no clear idea on how to set up seedmap.png - missed the memo on that, methinks - and just set a few pixels here and there.
eigs() solves this for the matrix’s Largest eigenvalues and vectors - Elevels sets the size of the ask - and don’t ask for all of them or you’ll be watching a spinning cursor for days.
As I noted before, I found eigs() to be unstable - it was not happy with the Hamiltonian being fed to it and crashed a lot. I’m pretty sure it’s because my seed map is not intelligently set up, but other parameters may have to be tweeked - i. e., we may have to go looking for parametric sweet spots. I’m not confident how G’MIC’s eig() is going to take to this, given how bad-tempered its scipy cousin is. There’s going to be the need for tuning, if it works at all. Not sure how to do that.
Finally, and - please! - I don’t mean this pejoratively! - this is “fun, Let’s explore this! Let’s explore that!” code. It is not efficient or “production quality.” It’s what springs up when lots of ideas are explored and some lines of investigation are dropped - but we don’t quite clean up the code. For example using np.linspace() just to get 1/(N-1)**2 is kind of a classic entry for the Rube Goldberg Machine Contest. But we’ve all written code like that and most of my G’MIC stuff springs from code like that. I’m not calling anybody out for having fun. I just raise this flag for the code translation effort: don’t take everything literally here.
Sorry I won’t be able to play in this sandbox anymore for a week or two - other projects call. Hope my notes are helpful. Have fun everybody!
Usually, the trick M + \lambda\;I_d where \lambda is some small positive value and I_d is the identity matrix is a good way to make it more stable 
(maybe not working here, I don’t know).
Yes… could be… but on the Python side, my problems may also stem from not using np.sparse.linalg.eigh(), which is designed for Hamiltonian matrices.
@David_Tschumperle: what’s the practical limit to size of matrices fed to eig()? 300,300,1,1? The vast, teeming G’MIC documentation team wants to know… 
As they say, don’t put your eigs in one basket ! 
…mutter… …mutter… I could also import @myselfhimself 'sgmic module into a Schrodinger’s Well IPython debug session and answer for myself how much/badly G’MIC’s eig() reacts to a Hamiltonian generated in @weightt_an code… Maybe there will be a little time tonight before bedtime… fingers crossed…
@grosgood There is also this command:
C:\Windows\System32>gmic h eigen
eigen (+):
Compute the eigenvalues and eigenvectors of selected symmetric matrices or
matrix fields.
If one selected image has 3 or 6 channels, it is regarded as a field of
2x2 or 3x3 symmetric matrices,
whose eigen elements are computed at each point of the field.
Example:
[#1] (1,0,0;0,2,0;0,0,3) +eigen
[#2] image.jpg structuretensors blur 2 eigen split[0] c
Tutorial: https://gmic.eu/oldtutorial/_eigen.shtml
EDIT: There seem to be this in the oldtutorial/eigen page that makes me wary:
@Reptorian
Yep. I know that guy from the Do Your Own Diffusion Tensor Fields Cookbook extravaganza. It may not have the capacity. @weightt_an test case is N=300, which sets the diagonal size of the Hamiltonian to N x N = 300 X 300 = 90,000 - a G’MIC image of 90000,90000,1,1 But @weightt_an also suggest in a post upstream from this (for trial sizes of N):
30 X 30 = 900 may bring us into the realm of -eigen, which tended to be unhappy after 1024x1024,1,1. Then again, that tutorial was written in 2013 on a 1.6 G’MIC release. Times have changed and @David_Tschumperle has bought the math processor and -eigen forward.
Anyhow, you were right in flagging diag() and eigs() - they sit at the computational heart of the graphic. Get past those, and I think you’ll have the makings of a very interesting filter. Hopefully, it won’t need a bazillion stabilization parameters…
@David_Tschumperle Do we have a easy way to get the diag vector from a set of vectors? It’s easy to do single vector though:
rep_diag:
{narg($*)},{narg($*)},1,1,"begin(vec_vals=[$*];);x==y?vec_vals[x];"
EDIT: Saw diag(V), but that doesn’t do from multiple vectors.
How about use diag(V) to take the five three 1x* vectors to square *x* images, then add them (arithmetically) each with offsets off the main diagonal, x+1, y+1, x+ 2N , y+2N, and for the other side of the main diagonal, x-1, y-1 and x-2N, y-2N. This feels doable in G’MIC. You’ll work with five three images of nearly the same size; each image having one of the former 1 X* vectors as its main diagonal, the image holding the main diagonal being biggest. Make a stack of them offset along the opposite diagonal, then add arithmetically, effectively flattening them into the main Hamitonian. Maybe do this with N=30 (900,900,1,1) images, or there abouts. Picture it?
Edit: Align with @weightt_an diagram.
Edit five three, because the two upper off-diagonals are duplicates of the two lower off-diagonals, being that the matrix as a whole is a Hamiltonian.
Nope still wrong
Looks like this
Great picture! Glad to see it!
Not what I was seeing last night in the debugger, which is what my image is based on, but that’s OK! Just a bug somewhere, that’s all. Maybe why I was seeing eigs() instability. Align things as you picture them here is a step forward. Thanks!
This is why in every discussion on the forum I ask for pretty pictures and numbers. Without them, it is easy to assume that we are talking about the same thing when we aren’t
(Every time I read Hamiltonian I keep on re-experiencing the adverts for and mentions of Hamilton the musical. Sorry, American friends, I am not remotedly interested.)
No sir. This is not Alexander Hamilton, the fellow on our (American) ten dollar bill, but Sir William Rowan Hamilton, a highly original and remarkable mind. Our great-grandchildren will still be using quaternions long after the musical is a forgotten bit of twenty-first century frivolity - that I’ll wager (Alas, I’m digressing off the thread again… Poor @Reptorian!)
Now, I think I know how to convert multiple set of vectors into diag.
So far, I got this.
rep_hamiltonian:
if s>1 to_gray. fi
n={$1}
npts={sqr($n)}
npts_n={$npts-$n}
npts_1={$npts-1}
#x_intervals image#
200,1,1,1,":begin(const ww=w-1);lerp(0,1,x/ww)"
increment={sqr(i(#-1,1))}
incrementValue={-1/$increment}
zeroV={4/$increment}
diag0={vector(#$npts,$zeroV)}
rm.
repeat $! l[$>]
$npts_n,1,1,1,"begin(
const inc_v="$incrementValue";
);
inc_v*i(#-1,x);"
$npts_1,1,1,1,"begin(
const inc_v="$incrementValue";
);
inc_v*i(#-2,x);"
diagmN=[{crop(#-2)}]
diagm1=[{crop(#-1)}]
rm[-2,-1]
size_diag={narg($diag0)}
$size_diag,$size_diag,1,1
5,1,1,1,"begin(
diagmN="$diagmN";
diagm1="$diagm1";
diag0=["$diag0"];
d_size="$size_diag";
initial_pos_x=[0,0,0,1,"$n"];
initial_pos_y=["$n",1,0,0,0];
);
xpos=initial_pos_x[x];
ypos=initial_pos_y[x];
for(p=0,p<d_size,p++,
x==0||x==4?(i(#-1,xpos,ypos)=diagmN[p];):
x==1||x==3?(i(#-1,xpos,ypos)=diagm1[p];):
(i(#-1,xpos,ypos)=diag0[p];);
xpos++;
ypos++;
if(xpos==d_size||ypos==d_size,break());
);
"
rm.
f[0] eig(crop(#-1))
endl done
Didn’t turn out how I wanted.
The 5,1,1,1 create a diag vector on the second last image.
Without the eig code. Here’s what it looks like:
From the chat @weightt_an and I had in discord, the second image needs a fix. But, it’s on the right track. On, +N/-N, there should be more. In other words, there should be more thickness. At this point, all there needs to be is adjusting to the code I have made, and see how to fix it.
OK - I believe the diagonals are off on displacement from the main diagonal. for N=30, the main diagonal should be 900 pixels long for a 900,900,1,1 Hamiltonian matrix image. You got that. There should be inner off-diagonals one pixel to the left and right of the main. You got that. The outer off-main diagonals should be offset N = 30. I’m running @weightt_an script in IPython debugger now and have created the Hamiltonian in his code - just before the eigs call. Snapshot of the array:
Okay, new code:
rep_hamiltonian:
if s>1 to_gray. fi
ge. {avg(im,iM)}
hyp_img={sqrt(sqr(w)+sqr(h))}
n={int(sqrt(($hyp_img^2)/2))}
npts={sqr($n)}
npts_n={$npts-$n}
npts_1={$npts-1}
#x_intervals image#
200,1,1,1,":begin(const ww=w-1);lerp(0,1,x/ww)"
increment={sqr(i(#-1,1))}
incrementValue={-1/$increment}
zeroV={4/$increment}
diag0={vector(#$npts,$zeroV)}
rm.
repeat $! l[$>]
$npts_n,1,1,1,"begin(
const inc_v="$incrementValue";
v=[crop(#-1)];
);
inc_v*v[x];"
$npts_1,1,1,1,"begin(
const inc_v="$incrementValue";
v=[crop(#-2)];
);
inc_v*v[x];"
diagmN=[{crop(#-2)}]
diagm1=[{crop(#-1)}]
rm[-2,-1]
size_diag={narg($diag0)}
$size_diag,$size_diag,1,1
5,1,1,1,"begin(
diagmN="$diagmN";
diagm1="$diagm1";
diag0=["$diag0"];
d_size="$size_diag";
initial_pos_x=[0,0,0,1,"$n"];
initial_pos_y=["$n",1,0,0,0];
);
xpos=initial_pos_x[x];
ypos=initial_pos_y[x];
for(p=0,p<d_size,p++,
x==0||x==4?(i(#-1,xpos,ypos)=diagmN[p];):
x==1||x==3?(i(#-1,xpos,ypos)=diagm1[p];):
(i(#-1,xpos,ypos)=diag0[p];);
xpos++;
ypos++;
if(xpos==d_size||ypos==d_size,break());
);
"
rm.
endl done
The command will automatically find the size of the Hamiltonian image. At this point, it’s close to being done.
Closeup:
Compared to Garry’s output, your line segments aren’t completely members of the 5 diagonals. It is probably a requirement for the filter to work as advertised. I am sure you are working on it. 
