Merge some classes, use memory-mapped storage, vectorise more

[mohawk2]

No description provided.

[coveralls]

Pull Request Test Coverage Report for Build 1301452136

• 478 of 482 (99.17%) changed or added relevant lines in 40 files are covered.
• No unchanged relevant lines lost coverage.
• Overall coverage decreased (-0.3%) to 99.719%

---

Changes Missing Coverage	Covered Lines	Changed/Added Lines	%
lib/Photonic/Roles/Field.pm	3	7	42.86%

Totals
Change from base Build 1295774575:	-0.3%
Covered Lines:	1422
Relevant Lines:	1426

---

💛 - Coveralls

[mohawk2]

For EpsTensor classes, the above makes epsA etc be mandatory not provided to an evaluate method. This makes the objects properly immutable and means if one wants a different epsTensor for a different epsilon, just construct a new object. This is often considered good OO practice.

[wlmb]

For EpsTensor classes, the above makes epsA etc be mandatory not provided to an evaluate method. This makes the objects properly immutable and means if one wants a different epsTensor for a different epsilon, just construct a new object. This is often considered good OO practice.

Hi. I just started looking at your PR. There is a problem with the mandatory epsA and epsB and the inmutable object. A very big advantage of the Haydock method applied to the 2 components nonretarded case (the NR2 objects) is that the same Haydock coefficients may be used for arbitrary compositions, i.e., they don't have to be recalculated when eps[AB] are modified. That was the logic for using an evaluate method, the idea being that on the first call the Haydock coefficients would be calculated (this is the most expensive part of the calculation), but on subsequent calls, they would be recycled. How would this work if a new object is built for every pair of eps[AB]?

On the other hand, I'm still somewhat worried about the states being a pdl mapped to a file. Currently there are options to don't keep states except for the current and neighbor states (to conserve memory, when other states are not needed), to keep them in memory (for speed) or to keep them in a file. In the last two cases the may be accessed through an iterator. If they become a single pdl, will there still be three choices? I tested some time ago the file mapped pdl's, and I don't recall the details, but I remember that I was unable to create a large mapped pdl although it was smaller than the available memory. I will have to search my notes or experiment again, but if true, this could be an issue, as the size of the states pdl can become huge.

[mohawk2]

Hi. I just started looking at your PR. There is a problem with the mandatory epsA and epsB and the immutable object. A very big advantage of the Haydock method applied to the 2 components nonretarded case (the NR2 objects) is that the same Haydock coefficients may be used for arbitrary compositions, i.e., they don't have to be recalculated when eps[AB] are modified. That was the logic for using an evaluate method, the idea being that on the first call the Haydock coefficients would be calculated (this is the most expensive part of the calculation), but on subsequent calls, they would be recycled. How would this work if a new object is built for every pair of eps[AB]?

As I understand it, those coefficients live in the AllH object, which lives in the immutable object as nr. That doesn't go anywhere, and can be used to construct another object with different epsA and epsB without recalculation.

On the other hand, I'm still somewhat worried about the states being a pdl mapped to a file. Currently there are options to don't keep states except for the current and neighbor states (to conserve memory, when other states are not needed), to keep them in memory (for speed) or to keep them in a file. In the last two cases the may be accessed through an iterator. If they become a single pdl, will there still be three choices? I tested some time ago the file mapped pdl's, and I don't recall the details, but I remember that I was unable to create a large mapped pdl although it was smaller than the available memory. I will have to search my notes or experiment again, but if true, this could be an issue, as the size of the states pdl can become huge.

Please do so! The whole idea behind File::Map is that you can memory-map any file, no matter how large, immediately and without consuming (or being limited by) physical RAM. If you're saying you have found this not to be the case, that's a bug that can be fixed if you can show me how to reproduce it :-)

[mohawk2]

pdl> use PDL::IO::FastRaw
pdl> $PDL::BIGPDL = 1
pdl> $on_disk = mapfraw('fname', {Creat => 1, Dims => [500,1000,1000], Datatype=>double});
pdl> p $on_disk->info
PDL: Double D [500,1000,1000]

The above works instantly, and creates a 4GB file (taking up no physical RAM). I did try with 1000^3 (8GB of data) and it crashed with an "out of memory" which makes sense because the machine only has 2GB of physical RAM and 2GB of VM. I then did these sanity-check operations:

pdl> $on_disk(0:-1:2).=1 # takes about 20 secs
pdl> p $on_disk(0:10,0,0)
[
 [
  [1 0 1 0 1 0 1 0 1 0 1]
 ]
]
# closed the shell, and restarted it
pdl> use PDL::IO::FastRaw
pdl> $PDL::BIGPDL = 1
pdl> $on_disk = mapfraw('fname')
pdl> p $on_disk(0:10,0,0)
[
 [
  [1 0 1 0 1 0 1 0 1 0 1]
 ]
]

If using Linux and Bash, please make sure your ulimit are set right (here, unlimited for the various memory settings).

[wlmb]

As I understand it, those coefficients live in the AllH object, which lives in the immutable object as nr. That doesn't go anywhere, and can be used to construct another object with different epsA and epsB without recalculation.

You are right, so the user would recycle the already evaluated nr for succesive calculations (I have not tested yet recycling evaluated AllH's, but I guess they should work). Nevertheless, what would happen for EpsTensor, where the nr array is built from the geometry upon initialization? Do you contemplate the user providing an already evaluated nr?

On the other hand, I'm still somewhat worried about the states being a pdl mapped to a file...

Please do so! The whole idea behind File::Map is that you can memory-map any file, no matter how large, immediately and without consuming (or being limited by) physical RAM. If you're saying you have found this not to be the case, that's a bug that can be fixed if you can show me how to reproduce it :-)

I don't know if this is what you intend to use, but consider:

    pdl> $PDL::BIGPDL=1
    pdl> use PDL::IO::FastRaw
    pdl> $a=mapfraw("rem", {Creat=>1,Dims=>[3,200,200,200,100], Datatype=>9})  
    Out of memory!

So I couldn't create a pdl that corresponds to the 100 states as 3D complex vector fields on a modest 200x200x200 grid. Curiously, I can create the corresponding pdl without mapfraw for, say, 30 states

    pdl> $a= zeroes(cdouble,3,200,200,200,30)

and it uses about 56% of my memory (on a 16Gb laptop). If I create the pdl with mapfraw the program uses less than 0.1% and is much much faster. So my uneducated guess is that mapfraw reserves the memory it would require to store the whole ndarray and the liberates it, but if it can't reserve it, it fails.

[wlmb]

The above works instantly, and creates a 4GB file (taking up no physical RAM). I did try with 1000^3 (8GB of data) and it crashed with an "out of memory" which makes sense because the machine only has 2GB of physical RAM and 2GB of VM. I then did these sanity-check operations:

Ok. So you did tests similar to what I was doing. Going beyond the available memory would be quite slow, but is possible with the previous 'iterator' solution. But it seems not to be possible with the mapfraw solution. This worries me for 3D calculations, as they may become unfeasible.

[mohawk2]

As I understand it, those coefficients live in the AllH object, which lives in the immutable object as nr. That doesn't go anywhere, and can be used to construct another object with different epsA and epsB without recalculation.
You are right, so the user would recycle the already evaluated nr for successive calculations (I have not tested yet recycling evaluated AllH's, but I guess they should work). Nevertheless, what would happen for EpsTensor, where the nr array is built from the geometry upon initialization? Do you contemplate the user providing an already evaluated nr?

You could go:

my $eto1 = Photonic::LE::NR2::EpsTensor(geometry=>$g, epsA=>$ea1, epsB=>$eb1);
my $ett1 = $eto1->epsTensor;
my $eto2 = Photonic::LE::NR2::EpsTensor(geometry=>$g, epsA=>$ea2, epsB=>$eb2, nr=>$eto1->nr);
my $ett2 = $eto2->epsTensor;

This facility comes "for free" as part of Moo*, because nr is an attribute that can be either supplied, or will be lazy-built from the geometry. If you want, I can make a but method (which would copy all attributes "but" replacing some, and with init_arg=>undef for epsTensor, it would automatically not be copied) so you could do the above as:

# make $eto1
my $eto2 = $eto1->but(epsA=>$ea2, epsB=>$eb2);

I'll answer the mapfraw point separately.

[mohawk2]

I don't know if this is what you intend to use, but consider:

    pdl> $PDL::BIGPDL=1
    pdl> use PDL::IO::FastRaw
    pdl> $a=mapfraw("rem", {Creat=>1,Dims=>[3,200,200,200,100], Datatype=>9})  
    Out of memory!

So I couldn't create a pdl that corresponds to the 100 states as 3D complex vector fields on a modest 200x200x200 grid. Curiously, I can create the corresponding pdl without mapfraw for, say, 30 states

    pdl> $a= zeroes(cdouble,3,200,200,200,30)

and it uses about 56% of my memory (on a 16Gb laptop). If I create the pdl with mapfraw the program uses less than 0.1% and is much much faster. So my uneducated guess is that mapfraw reserves the memory it would require to store the whole ndarray and the liberates it, but if it can't reserve it, it fails.

That would also be my uneducated guess! I will take a look at this now.

[wlmb]

You are right, so the user would recycle the already evaluated nr for successive calculations (I have not tested yet recycling evaluated AllH's, but I guess they should work). Nevertheless, what would happen for EpsTensor, where the nr array is built from the geometry upon initialization? Do you contemplate the user providing an already evaluated nr?

You could go:

my $eto1 = Photonic::LE::NR2::EpsTensor(geometry=>$g, epsA=>$ea1, epsB=>$eb1);
my $ett1 = $eto1->epsTensor;
my $eto2 = Photonic::LE::NR2::EpsTensor(geometry=>$g, epsA=>$ea2, epsB=>$eb2, nr=>$eto1->nr);
my $ett2 = $eto2->epsTensor;

This facility comes "for free" as part of Moo*, because nr is an attribute that can be either supplied, or will be lazy-built from the geometry.

Ok. In this case, init_arg=>undef would have to go.

If you want, I can make a but method (which would copy all attributes "but" replacing some, and with init_arg=>undef for epsTensor, it would automatically not be copied) so you could do the above as:

# make $eto1
my $eto2 = $eto1->but(epsA=>$ea2, epsB=>$eb2);

The 'but' method seems a clean solution, so you don't have to carry all the other initialization parameters. Is this kind of solution in common usage? I didn't understand your remark above about the init_arg=>undef.

[mohawk2]

and it uses about 56% of my memory (on a 16Gb laptop). If I create the pdl with mapfraw the program uses less than 0.1% and is much much faster. So my uneducated guess is that mapfraw reserves the memory it would require to store the whole ndarray and the liberates it, but if it can't reserve it, it fails.

You've won a valuable prize! A new version of PDL with a fix will be on its way to CPAN as soon as it's passed CI. I knew there'd be at least one problem revealed by actually using it here! (There are also some other tweaks to the mmap-type stuff which I'd already made but not released yet)

[mohawk2]

This facility comes "for free" as part of Moo*, because nr is an attribute that can be either supplied, or will be lazy-built from the geometry.

Ok. In this case, init_arg=>undef would have to go.

Quite right.

If you want, I can make a but method (which would copy all attributes "but" replacing some, and with init_arg=>undef for epsTensor, it would automatically not be copied) so you could do the above as:

# make $eto1
my $eto2 = $eto1->but(epsA=>$ea2, epsB=>$eb2);

The 'but' method seems a clean solution, so you don't have to carry all the other initialization parameters. Is this kind of solution in common usage? I didn't understand your remark above about the init_arg=>undef.

I think you did in fact understand it, because you spotted the accidental contradiction, sorry.

The but method would need the init_arg to be undef so that a blind copy would leave out that attribute, whereas with the approach that uses the current code as-is, you'd just supply the nr attribute value. The two approaches are incompatible, and I'm inclined to not have this but method here. If we did, we'd have to have maybe an optional first param which would be an array-ref with attributes to leave out of the copy. It's workable, but not very elegant, and given that here the benefit would only be to not have to pass the geometry attribute, doesn't seem worthwhile.

[mohawk2]

and it uses about 56% of my memory (on a 16Gb laptop). If I create the pdl with mapfraw the program uses less than 0.1% and is much much faster. So my uneducated guess is that mapfraw reserves the memory it would require to store the whole ndarray and the liberates it, but if it can't reserve it, it fails.

You've won a valuable prize! A new version of PDL with a fix will be on its way to CPAN as soon as it's passed CI. I knew there'd be at least one problem revealed by actually using it here! (There are also some other tweaks to the mmap-type stuff which I'd already made but not released yet)

This has now been uploaded as 2.056 (which is quite a nice number, sort of), and on my machine I can now mapfraw a 16GB ndarray with no problems. Please try it :-)

[wlmb]

You've won a valuable prize! A new version of PDL with a fix will be on its way to CPAN as soon as it's passed CI. I knew there'd be at least one problem revealed by actually using it here! (There are also some other tweaks to the mmap-type stuff which I'd already made but not released yet)

This has now been uploaded as 2.056 (which is quite a nice number, sort of), and on my machine I can now mapfraw a 16GB ndarray with no problems. Please try it :-)

Ha! You're fast! Thanks for the prize! I appreciate it

[mohawk2]

Do you find it works right with mapfraw, and if so are you happy with the approach of memory-mapping the whole states? Are you content with the nr reuse? I'm putting a note in the Photonic::Roles::EpsTensor I'm just making now.

[wlmb]

Do you find it works right with mapfraw, and if so are you happy with the approach of memory-mapping the whole states? Are you content with the nr reuse? I'm putting a note in the Photonic::Roles::EpsTensor I'm just making now.

I'm just installing 2.056, so I haven't tested the new mapfraw, but after removing the limitation on the size of mapped ndarrays I guess memory mapping will be a good solution. Also, I'm content with the nr reuse, though my colleagues will have to change their codes. Maybe a temporal 'evaluate' that dies with a useful message would facilitate the transition.

[mohawk2]

As I'm re-reading the code again after making this EpsTensor role, I think there are more changes to be made that will make for much less code for your colleagues to need to write! For instance, I'm seeing the array of EpsL, and now starting to imagine how to vectorise the calculations of things.

[wlmb]

I just tested the new mapfraw. Seems to work nicely and fast.
In my implementation, epsA and epsB were assumed to be scalars. I guess they could be ndarrays themselves, so that a single epsTensor object could calculate a whole spectrum instead of a single number. Also, an array of characteristic functions B instead of a single function could be used to calculate in one step different geometries. The current code would not work as it assumes all dimensions of B are geometric dimensions, and uses functions such as 'sum' which consume all dimensions.

[mohawk2]

Also, more immediately: I don't see a need for a dummy evaluate, because you can't instantiate one of these objects without supplying the epsA (if needed), so at least that will be very clear.

[mohawk2]

I just tested the new mapfraw. Seems to work nicely and fast.

Great!

In my implementation, epsA and epsB were assumed to be scalars. I guess they could be ndarrays themselves, so that a single epsTensor object could calculate a whole spectrum instead of a single number. Also, an array of characteristic functions B instead of a single function could be used to calculate in one step different geometries. The current code would not work as it assumes all dimensions of B are geometric dimensions, and uses functions such as 'sum' which consume all dimensions.

Yes, there might be some adjustments needed. My feeling however is it will be surprisingly small, and then... threading! Multi-processing! ;-)

[wlmb]

Also, more immediately: I don't see a need for a dummy evaluate, because you can't instantiate one of these objects without supplying the epsA (if needed), so at least that will be very clear.

True...

[wlmb]

Yes, there might be some adjustments needed. My feeling however is it will be surprisingly small, and then... threading! Multi-processing! ;-)

That'll be nice.

[wlmb]

I just tested the new mapfraw. Seems to work nicely and fast.
Now I don't have to set $PDL::BIGPDL=1 to use big but mapped pdl's

[mohawk2]

I just tested the new mapfraw. Seems to work nicely and fast.
Now I don't have to set $PDL::BIGPDL=1 to use big but mapped pdl's

Ah, of course - the BIGPDL check is not done unless it's memory!

[wlmb]

Ah, of course - the BIGPDL check is not done unless it's memory!
The previous behavior (requiring BIGPDL) was misterious.

[mohawk2]

Ah, of course - the BIGPDL check is not done unless it's memory!

The previous behavior (requiring BIGPDL) was mysterious.

It's a sign of the brief allocation of memory that was being done. Not so mysterious, but still bad ;-)

[mohawk2]

The most recent commit renames the nr attribute to haydock since that's much less confusing for a newbie like me (and arrays of AllH objects are called that in eg the Green* classes).

[mohawk2]

The latest commit, switching from Moose to Moo (a lighter-weight system), makes the tests go from taking 12-13secs, to 9-10secs, which aids development.

[mohawk2]

The most recent commit makes a Roles::Field.

[mohawk2]

@wlmb There are a couple more changes I believe would be valuable, but I'd prefer to get your feedback on these changes so far first?

[wlmb]

Dear Ed,

On Sun, Dec 26, 2021 at 03:06:15PM -0800, mohawk2 wrote: @wlmb There are a couple more changes I believe would be valuable, but I'd prefer to get your feedback on these changes so far first?

Thanks! I'm sorry I have not reviewed all your previous changes. It was a busy end of semester. I expect to look at the previous changes carefully these coming days. Best regards and my best wishes. Luis

…

-- Reply to this email directly or view it on GitHub: #23 (comment) You are receiving this because you were mentioned. Message ID: ***@***.***>

-- o W. Luis Mochán, | tel:(52)(777)329-1734 /<(*) Instituto de Ciencias Físicas, UNAM | fax:(52)(777)317-5388 `>/ /\ Av. Universidad s/n CP 62210 | (*)/\/ \ Cuernavaca, Morelos, México | ***@***.*** /\_/\__/ GPG: 791EB9EB, C949 3F81 6D9B 1191 9A16 C2DF 5F0A C52B 791E B9EB

[mohawk2]

No problem! Glad to hear back from you already and I look forward to your thoughts when you have time. My hope is that these changes so far are acceptable, though naturally you may have questions and/or spot things I missed. Once this change-set is acceptable, then I'll look at the dimension-reordering stuff (basically, the aim is to avoid needing to move dimensions around for threading, which shows the dimension layout wasn't right yet).

[mohawk2]

@wlmb Could you take a look at these changes? :-)