Re: [Numpy-discussion] f2py with allocatable arrays
On 4. juli 2012, at 02:23, Sturla Molden wrote: > Den 03.07.2012 20:38, skrev Casey W. Stark: >> >> Sturla, this is valid Fortran, but I agree it might just be a bad >> idea. The Fortran 90/95 Explained book mentions this in the >> allocatable dummy arguments section and has an example using an array >> with allocatable, intent(out) in a subrountine. You can also see this >> in the PDF linked from >> http://fortranwiki.org/fortran/show/Allocatable+enhancements. > > Ok, so it's valid Fortran 2003. I never came any longer than to Fortran > 95 :-) Make sure any Fortran code using this have the extension .f03 -- > not .f95 or .f90 -- or it might crash horribly. > To be pedantic: to my knowledge, the common convention is .f for fixed and .f for free form source code. As is stated in the link, "..the Fortran standard itself does not define any extension..." http://fortranwiki.org/fortran/show/File+extensions As one example, ifort doesn't even want to read files with the .f95 suffix. You'll have to pass it a flag stating that "yep, that's a fortran file all right". I use the .f90 suffix everywhere, but maybe that's just me. Paul ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] f2py with allocatable arrays
Den 04.07.2012 01:59, skrev Sturla Molden: > But neither was the case here. The allocatable was a dummy variable in > a subroutine's interface, declared with intent(out). That is an error > the compiler should trap, because it is doomed to segfault. Ok, so the answer here seems to be: In Fortran 90 is compiler dependent. In Fortran 95 it is an error. In extensions to Fortran 95 it is legal. In Fortran 2003 and 2008 it is legal. Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] f2py with allocatable arrays
Den 03.07.2012 20:38, skrev Casey W. Stark: > > Sturla, this is valid Fortran, but I agree it might just be a bad > idea. The Fortran 90/95 Explained book mentions this in the > allocatable dummy arguments section and has an example using an array > with allocatable, intent(out) in a subrountine. You can also see this > in the PDF linked from > http://fortranwiki.org/fortran/show/Allocatable+enhancements. Ok, so it's valid Fortran 2003. I never came any longer than to Fortran 95 :-) Make sure any Fortran code using this have the extension .f03 -- not .f95 or .f90 -- or it might crash horribly. > I have never used array pointers in Fortran, but I might give this a > shot. Are there any problems returning pointers to arrays back to > python though? In the Fortran 2003 ISO C binding you can find the methods F_C_POINTER and C_F_POINTER that will convert between C and Fortran pointers. A Fortran pointer is an array struct, much like a NumPy view array, and has dimensions and strides. It can reference contiguous and discontiguous blocks of memory. A C pointer is just a memory address. You must therefore use special methods to convert between C and Fortran pointers. There is an extension to Cython called fwrap that will generate this kind of boiler-plate code for conversion between NumPy and Fortran using the ISO C bindings in Fortran 2003. It is an alternative to f2py, though less mature. At the binary level, a Fortran pointer and an allocatable array usually have the same representation. But there are semantic differences between them. In Fortran 2003, pointers are less useful than in Fortran 95, except when interfacing with C through the ISO C bindings. That is because you can put allocatable arrays as memers in derived types, and (as I learned today) use them as dummy variables. In Fortran 95 those would be common cases for using pointers. By the way: The Fortran counter part to a C pointer is a Cray pointer, which is a common non-standard extension to the language. A Cray pointer, when supported, is a pair of variables: one storing the address and the other dereferencing the address. Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] f2py with allocatable arrays
Den 03.07.2012 19:24, skrev Pearu Peterson: > > One can have allocatable arrays in module data block, for instance, where > they a global In Fortran 2003 one can also have allocatable arrays as members in derived types. But neither was the case here. The allocatable was a dummy variable in a subroutine's interface, declared with intent(out). That is an error the compiler should trap, because it is doomed to segfault. Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] "import numpy" performance
On Mon, Jul 2, 2012 at 2:59 PM, Nathaniel Smith wrote:
> On Mon, Jul 2, 2012 at 10:06 PM, Robert Kern wrote:
>> On Mon, Jul 2, 2012 at 9:43 PM, Benjamin Root wrote:
>>>
>>> On Mon, Jul 2, 2012 at 4:34 PM, Nathaniel Smith wrote:
>>
I think this ship has sailed, but it'd be worth looking into lazy
importing, where 'numpy.fft' isn't actually imported until someone
starts using it. There are a bunch of libraries that do this, and one
would have to fiddle to get compatibility with all the different
python versions and make sure you're not killing performance (might
have to be in C) but something along the lines of
class _FFTModule(object):
def __getattribute__(self, name):
mod = importlib.import_module("numpy.fft")
_FFTModule.__getattribute__ = mod.__getattribute__
return getattr(mod, name)
fft = _FFTModule()
>>>
>>> Not sure how this would impact projects like ipython that does
>>> tab-completion support, but I know that that would drive me nuts in my basic
>>> tab-completion setup I have for my regular python terminal. Of course, in
>>> the grand scheme of things, that really isn't all that important, I don't
>>> think.
>>
>> We used to do it for scipy. It did interfere with tab completion. It
>> did drive many people nuts.
>
> Sounds like a bug in your old code, or else the REPLs have gotten
> better? I just pasted the above code into both ipython and python
> prompts, and typing 'fft.' worked fine in both cases. dir(fft)
> works first try as well.
>
> (If you try this, don't forget to 'import importlib' first, and note
> importlib is 2.7+ only. Obviously importlib is not necessary but it
> makes the minimal example less tedious.)
For anyone interested, I worked out a small lazy-loading class that
we use in nitime [1], which does not need importlib and thus works
on python versions before 2.7 and also has a bit of repr pretty
printing.
I wrote about this to Scipy-Dev [2], and in the original nitime
PR [3] tested that it works in python 2.5, 2.6, 2.7, 3.0, 3.1 and
3.2.
Since that time, we've only changed how we deal with the one
known limitation: reloading a lazily-loaded module was a noop in
that PR, but now throws an error (there's one line commented out
if the noop behavior is preferred).
Here's a link to the rendered docs [4], but if you just grab the
LazyImport class from [1], you can do
fft = LazyImport('numpy.fft')
1. https://github.com/nipy/nitime/blob/master/nitime/lazyimports.py
2. http://mail.scipy.org/pipermail/scipy-dev/2011-September/016606.html
3. https://github.com/nipy/nitime/pull/88
4.
http://nipy.sourceforge.net/nitime/api/generated/nitime.lazyimports.html#module-nitime.lazyimports
best,
--
Paul Ivanov
314 address only used for lists, off-list direct email at:
http://pirsquared.org | GPG/PGP key id: 0x0F3E28F7
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Buildbot status
On Mon, Jul 2, 2012 at 6:31 PM, Travis Oliphant wrote: > Ondrej should have time to work on this full time in the coming days. That's great; having Ondrej on this full time will help a great deal. > NumFocus can provide some funding needed for maintaining servers, etc, but > keeping build bots active requires the efforts of multiple volunteers. > If anyone has build machines to offer, please let Ondrej know so he can > coordinate getting Jenkins slaves onto them and hooking them up to the master. I'd be glad if we could discuss it mainly on list, just to keep everyone in the loop. For now, I think we need to answer the two questions mentioned above: 1) What happens to the current installation on buildbot.scipy.org? 2) If we're not keeping buildbot, or if we want additional systems, which ones should we use? Jenkins? and then also 3) Which build slaves should we employ? We have the current build slaves, the nipy ones have been volunteered, and then there's the GCC build farm mentioned by Fernando. Ondrej, perhaps you can comment on what you had in mind? If we have a clear plan of action before you start off, we can all help out in putting the pieces together. Stéfan ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Would a patch with a function for incrementing an array with advanced indexing be accepted?
Hi, Here is code example that work only with different index: import numpy x=numpy.zeros((5,5)) x[[0,2,4]]+=numpy.random.rand(3,5) print x This won't work if in the list [0,2,4], there is index duplication, but with your new code, it will. I think it is the most used case of advanced indexing. At least, for our lab:) Fred On Mon, Jul 2, 2012 at 7:48 PM, John Salvatier wrote: > Hi Fred, > > That's an excellent idea, but I am not too familiar with this use case. What > do you mean by list in 'matrix[list]'? Is the use case, just incrementing > in place a sub matrix of a numpy matrix? > > John > > > On Fri, Jun 29, 2012 at 11:43 AM, Frédéric Bastien wrote: >> >> Hi, >> >> I personnaly can't review this as this is too much in NumPy internal. >> >> My only comments is that you could add a test and an example in the >> doc for matrix[list]. I think it will be the most used case. >> >> Fred >> >> On Wed, Jun 27, 2012 at 7:47 PM, John Salvatier >> wrote: >> > I've submitted a pull request ( https://github.com/numpy/numpy/pull/326 >> > ). >> > I'm new to the numpy and python internals, so feedback is greatly >> > appreciated. >> > >> > >> > On Tue, Jun 26, 2012 at 12:10 PM, Travis Oliphant >> > wrote: >> >> >> >> >> >> On Jun 26, 2012, at 1:34 PM, Frédéric Bastien wrote: >> >> >> >> > Hi, >> >> > >> >> > I think he was referring that making NUMPY_ARRAY_OBJECT[...] syntax >> >> > support the operation that you said is hard. But having a separate >> >> > function do it is less complicated as you said. >> >> >> >> Yes. That's precisely what I meant. Thank you for clarifying. >> >> >> >> -Travis >> >> >> >> > >> >> > Fred >> >> > >> >> > On Tue, Jun 26, 2012 at 1:27 PM, John Salvatier >> >> > wrote: >> >> >> Can you clarify why it would be super hard? I just reused the code >> >> >> for >> >> >> advanced indexing (a modification of PyArray_SetMap). Am I missing >> >> >> something >> >> >> crucial? >> >> >> >> >> >> >> >> >> >> >> >> On Tue, Jun 26, 2012 at 9:57 AM, Travis Oliphant >> >> >> >> >> >> wrote: >> >> >>> >> >> >>> >> >> >>> On Jun 26, 2012, at 11:46 AM, John Salvatier wrote: >> >> >>> >> >> >>> Hello, >> >> >>> >> >> >>> If you increment an array using advanced indexing and have repeated >> >> >>> indexes, the array doesn't get repeatedly >> >> >>> incremented, >> >> >>> http://comments.gmane.org/gmane.comp.python.numeric.general/50291. >> >> >>> I wrote a C function that does incrementing with repeated indexes >> >> >>> correctly. >> >> >>> The branch is here (https://github.com/jsalvatier/numpy see the >> >> >>> last >> >> >>> two >> >> >>> commits). Would a patch with a cleaned up version of a function >> >> >>> like >> >> >>> this be >> >> >>> accepted into numpy? I'm not experienced writing numpy C code so >> >> >>> I'm >> >> >>> sure it >> >> >>> still needs improvement. >> >> >>> >> >> >>> >> >> >>> This is great. It is an often-requested feature. It's *very >> >> >>> difficult* >> >> >>> to do without changing fundamentally what NumPy is. But, yes this >> >> >>> would be >> >> >>> a great pull request. >> >> >>> >> >> >>> Thanks, >> >> >>> >> >> >>> -Travis >> >> >>> >> >> >>> >> >> >>> >> >> >>> ___ >> >> >>> NumPy-Discussion mailing list >> >> >>> NumPy-Discussion@scipy.org >> >> >>> http://mail.scipy.org/mailman/listinfo/numpy-discussion >> >> >>> >> >> >> >> >> >> >> >> >> ___ >> >> >> NumPy-Discussion mailing list >> >> >> NumPy-Discussion@scipy.org >> >> >> http://mail.scipy.org/mailman/listinfo/numpy-discussion >> >> >> >> >> > ___ >> >> > NumPy-Discussion mailing list >> >> > NumPy-Discussion@scipy.org >> >> > http://mail.scipy.org/mailman/listinfo/numpy-discussion >> >> >> >> ___ >> >> NumPy-Discussion mailing list >> >> NumPy-Discussion@scipy.org >> >> http://mail.scipy.org/mailman/listinfo/numpy-discussion >> > >> > >> > >> > ___ >> > NumPy-Discussion mailing list >> > NumPy-Discussion@scipy.org >> > http://mail.scipy.org/mailman/listinfo/numpy-discussion >> > >> ___ >> NumPy-Discussion mailing list >> NumPy-Discussion@scipy.org >> http://mail.scipy.org/mailman/listinfo/numpy-discussion > > > > ___ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion > ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] f2py with allocatable arrays
Hi all. Thanks for the speedy responses! I'll try to respond to all... The first idea is to split up the routine into two -- one to compute the final size of the arrays, and the second to fill them in. I might end up doing this, because it is simplest, but it means creating the initial conditions twice, throwing them away the first time. Not a huge deal because this setup is not that big (a few minutes), but it will take twice as long. Sturla, this is valid Fortran, but I agree it might just be a bad idea. The Fortran 90/95 Explained book mentions this in the allocatable dummy arguments section and has an example using an array with allocatable, intent(out) in a subrountine. You can also see this in the PDF linked from http://fortranwiki.org/fortran/show/Allocatable+enhancements. I have never used array pointers in Fortran, but I might give this a shot. Are there any problems returning pointers to arrays back to python though? As for storing the arrays in the module data block, I would like to avoid that approach if possible. It doesn't make sense for these to be module level components when the size and values depend on the input to a subroutine in the module. Best, Casey On Tue, Jul 3, 2012 at 10:24 AM, Pearu Peterson wrote: > > > On Tue, Jul 3, 2012 at 5:20 PM, Sturla Molden wrote: > >> >> As for f2py: Allocatable arrays are local variables for internal use, >> and they are not a part of the subroutine's calling interface. f2py only >> needs to know about the interface, not the local variables. >> > > One can have allocatable arrays in module data block, for instance, where > they a global. f2py supports wrapping these allocatable arrays to python. > See, for example, > > > http://cens.ioc.ee/projects/f2py2e/usersguide/index.html#allocatable-arrays > > Pearu > > > ___ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion > > ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] f2py with allocatable arrays
On Tue, Jul 3, 2012 at 5:20 PM, Sturla Molden wrote: > > As for f2py: Allocatable arrays are local variables for internal use, > and they are not a part of the subroutine's calling interface. f2py only > needs to know about the interface, not the local variables. > One can have allocatable arrays in module data block, for instance, where they a global. f2py supports wrapping these allocatable arrays to python. See, for example, http://cens.ioc.ee/projects/f2py2e/usersguide/index.html#allocatable-arrays Pearu ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] "import numpy" performance
On Mon, Jul 2, 2012 at 12:17 PM, Andrew Dalke wrote: > In this email I propose a few changes which I think are minor > and which don't really affect the external NumPy API but which > I think could improve the "import numpy" performance by at > least 40%. +1 -- I think I remember that thread -- at the time, I was experiencing some really, really slow inport times myself -- it turned out to be something really wierd with my system (though I don't remember exactly what), but numpy still is too big an import. Another note -- I ship stuff with py2exe and friends a fair bit -- numpy's "Import a whole bunch of stuff you may well not be using" approach means I have to include all that stuff, or hack the heck out of numpy -- not ideal. > 1) remove "add_newdocs" and put the docstrings in the C code > 'add_newdocs' still needs to be there, > > The code says: > > # This is only meant to add docs to objects defined in C-extension modules. > # The purpose is to allow easier editing of the docstrings without > # requiring a re-compile. +1 -- isn't it better for the docs to be with the code, anyway? > 2) Don't optimistically assume that all submodules are > needed. For example, some current code uses > import numpy numpy.fft.ifft > +1 see above -- really, what fraction of code uses fft and polynomial, and ... "namespaces are one honking great idea" I appreciate the legacy, and the easy-of-use at the interpreter, but it sure would be nice to clean this up -- maybe keep the leegacy by having a new import: import just_numpy as np that would import the core stuff, and offer the "extra" packages as specific imports -- ideally, we'd dpreciate the old way, and reccoment the extra importing for the future, and some day have "numpy" and "numpy_plus". (Kind of like pylab, I suppose) lazy importing may work OK, too, though more awkward for py2exe and friends, and perhaps a bit "magic" for my taste. > 3) Especially: don't always import 'numpy.testing' +1 > I have not worried about numpy import performance for > 4 years. While I have been developing scientific software > for 20 years, and in Python for 15 years, it has been > in areas of biology and chemistry which don't use arrays. remarkable -- I use arrays for everything! most of which are not classic big arrays you process with lapack type stuff ;-) > yeah, it's just using the homogenous array most of the time. exactly -- I know Travis says: "if you're going to use numpy arrays, use numpy", but they really are pretty darn handy even if you just use them as containers. Ben root wrote: > Not sure how this would impact projects like ipython that does tab-completion > support, > but I know that that would drive me nuts in my basic tab-completion setup I > have for >my regular python terminal. Of course, in the grand scheme of things, that >really > isn't all that important, I don't think. I do think it's important to support easy interactive use, Ipyhton, etc -- with nice tab completion, easy access to doc string, etc. But it should alo be possible to not have all that where it isn't required -- hence my "import numpy_plus" type proposal. I never did get why the polynomial stuff was added to core numpy -Chris -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception chris.bar...@noaa.gov ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Change in memmap behaviour
On Tue, Jul 3, 2012 at 10:35 AM, Thouis (Ray) Jones wrote:
> On Mon, Jul 2, 2012 at 11:52 PM, Sveinung Gundersen
> wrote:
>>
>> On 2. juli 2012, at 22.40, Nathaniel Smith wrote:
>>
>>> On Mon, Jul 2, 2012 at 6:54 PM, Sveinung Gundersen
>>> wrote:
[snip]
Your actual memory usage may not have increased as much as you think,
since memmap objects don't necessarily take much memory -- it sounds
like you're leaking virtual memory, but your resident set size
shouldn't go up as much.
As I understand it, memmap objects retain the contents of the memmap in
memory after it has been read the first time (in a lazy manner). Thus, when
reading a slice of a 24GB file, only that part recides in memory. Our
system
reads a slice of a memmap, calculates something (say, the sum), and then
deletes the memmap. It then loops through this for consequitive slices,
retaining a low memory usage. Consider the following code:
import numpy as np
res = []
vecLen = 3095677412
for i in xrange(vecLen/10**8+1):
x = i * 10**8
y = min((i+1) * 10**8, vecLen)
res.append(np.memmap('val.float64', dtype='float64')[x:y].sum())
The memory usage of this code on a 24GB file (one value for each nucleotide
in the human DNA!) is 23g resident memory after the loop is finished (not
24g for some reason..).
Running the same code on 1.5.1rc1 gives a resident memory of 23m after the
loop.
>>>
>>> Your memory measurement tools are misleading you. The same memory is
>>> resident in both cases, just in one case your tools say it is
>>> operating system disk cache (and not attributed to your app), and in
>>> the other case that same memory, treated in the same way by the OS, is
>>> shown as part of your app's resident memory. Virtual memory is
>>> confusing...
>>
>> But the crucial difference is perhaps that the disk cache can be cleared by
>> the OS if needed, but not the application memory in the same way, which must
>> be swapped to disk? Or am I still confused?
>>
>> (snip)
>>
Great! Any idea on whether such a patch may be included in 1.7?
>>>
>>> Not really, if I or you or someone else gets inspired to take the time
>>> to write a patch soon then it will be, otherwise not...
>>>
>>> -N
>>
>> I have now tried to add a patch, in the way you proposed, but I may have
>> gotten it wrong..
>>
>> http://projects.scipy.org/numpy/ticket/2179
>
> I put this in a github repo, and added tests (author credit to Sveinung)
> https://github.com/thouis/numpy/tree/mmap_children
>
> I'm not sure which branch to issue a PR request against, though.
Looks good to me, thanks to both of you!
Obviously should be merged to master; beyond that I'm not sure. We
definitely want it in 1.7, but I'm not sure if that's been branched
yet or not. (Or rather, it has been branched, but then maybe it was
unbranched again? Travis?) Since it was a 1.6 regression it'd make
sense to cherrypick to the 1.6 branch too, just in case it gets
another release.
-n
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] f2py with allocatable arrays
Den 03.07.2012 11:54, skrev George Nurser: >> module zp >>implicit none >>contains >>subroutine ics(..., num_particles, particle_mass, positions, velocities) >> use data_types, only : dp >> implicit none >> ... inputs ... >> integer, intent(out) :: num_particles >> real (kind=dp), intent(out) :: particle_mass >> real (kind=dp), intent(out), dimension(:, :), allocatable :: positions, >> velocities >> ... >>end subroutine >> end module >> This looks like a Fortran error. deallocate is called automatically on allocatable arrays when a subroutine exit (guarranteed from Fortran 95, optional in Fortran 90). Allocatable arrays might even be put on the stack depending on the requested size. I am not even sure intent(out) on an allocatable array is legal or an error the compiler should trap, but it is obviously very bad Fortran. So "positions" and "velocities" should here be declared pointer, not allocatable. As for f2py: Allocatable arrays are local variables for internal use, and they are not a part of the subroutine's calling interface. f2py only needs to know about the interface, not the local variables. Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] f2py with allocatable arrays
On 7/2/2012 7:17 PM, Casey W. Stark wrote: Hi numpy. Does anyone know if f2py supports allocatable arrays, allocated inside fortran subroutines? The old f2py docs seem to indicate that the allocatable array must be created with numpy, and dropped in the module. Here's more background to explain... I have a fortran subroutine that returns allocatable positions and velocities arrays. I wish I could get rid of the allocatable part, but you don't know how many particles it will create until the subroutine does some work (it checks if each particle it perturbs ends up in the domain). module zp implicit none contains subroutine ics(..., num_particles, particle_mass, positions, velocities) use data_types, only : dp implicit none ... inputs ... integer, intent(out) :: num_particles real (kind=dp), intent(out) :: particle_mass real (kind=dp), intent(out), dimension(:, :), allocatable :: positions, velocities ... end subroutine end module I tested this with a fortran driver program and it looks good, but when I try with f2py, it cannot compile. It throws the error "Error: Actual argument for 'positions' must be ALLOCATABLE at (1)". I figure this has something to do with the auto-generated "*-f2pywrappers2.f90" file, but the build deletes the file. If anyone knows an f2py friendly way to rework this, I would be happy to try. I'm also fine with using ctypes if it can handle this case. Best, Casey ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion Hi, ... not sure what version of numpy you are using but it can be done with the f2py included with numpy 1.6.x. Here is a (hopefully intelligible) code fragment from a working application: module Grid_ ! a representation of a two-dimensional, rectangular, (image) sensor grid ! -- this component is used in the "smoothing" portion of the Thematic MaP solution algorithm (smoothness prior probabilities) use runtime_ implicit none integer, private, parameter :: begin_ = 1, end_ = 2, row_ = 1, column_ = 2 integer, private, dimension(begin_:end_) :: rows__ = (/0,0/), columns__ = (/0,0/) ! grid rows and columns extent integer, allocatable, dimension(:,:) :: pixel_neighbors ! the neighbors of a specified pixel -- ((row,column), ..., (row,column)) !+ this array is managed (allocated, deallocated, populated) by the neighbors_of() procedure !+ when allocated, the first dimension is always 2 -- pixel (row,column) !+ it would be preferable for this array to be the return result, of procedure neighbors_of(), ! but F2PY does not seem to support an allocatable array as a function result contains [snip] function neighbors_of (pixel) result(completed) ! determines all second-order (grid) neighbors of *pixel* !-- returns .true. if successful !-- use the error_diagnostic() procedure if .false. is returned !-- upon successful completion of this procedure, Grid_.pixel_neighbors contains the neighbors of *pixel* !-- *pixel* may not be on the grid; in this case, Grid_.pixel_neighbors is not allocated and .false. is returned ! integer, dimension(row_:column_), intent(in) :: pixel ! pixel to be considered ... f2py does not support this integer, dimension(2), intent(in) :: pixel ! pixel under consideration (row,column) logical :: completed ! integer, dimension(begin_:end_) :: neighbor_rows, neighbor_columns ... f2py does not support this ! integer, dimension(row_:column_) :: neighbor ... f2py does not support this integer, dimension(2) :: neighbor_rows, neighbor_columns ! each is: (begin,end) integer, dimension(2) :: neighbor ! (row,column) integer :: row, column, code, count character (len=100) :: diagnostic completed = .false. count = 0 ! *pixel* has no neighbors if (allocated (pixel_neighbors)) deallocate (pixel_neighbors) if (is_interior (pixel)) then count = 8 ! interior pixels have eight, second-order neighbors neighbor_rows(begin_) = pixel(row_)-1 neighbor_rows(end_) = pixel(row_)+1 neighbor_columns (begin_) = pixel(column_)-1 neighbor_columns (end_) = pixel(column_)+1 else if (is_border (pixel)) then count = 5 ! non-corner, border pixels have five, second-order neighbors, but ... if (is_corner (pixel)) count = 3 ! corner pixels have three, second-order neighbors neighbor_rows(begin_) = max (pixel(row_)-1, rows__(begin_)) neighbor_rows(end_) = min (p
Re: [Numpy-discussion] f2py with allocatable arrays
> Hi numpy. > > Does anyone know if f2py supports allocatable arrays, allocated inside > fortran subroutines? The old f2py docs seem to indicate that the > allocatable array must be created with numpy, and dropped in the module. > Here's more background to explain... > > I have a fortran subroutine that returns allocatable positions and > velocities arrays. I wish I could get rid of the allocatable part, but you > don't know how many particles it will create until the subroutine does > some > work (it checks if each particle it perturbs ends up in the domain). > > module zp > implicit none > contains > subroutine ics(..., num_particles, particle_mass, positions, velocities) > use data_types, only : dp > implicit none > ... inputs ... > integer, intent(out) :: num_particles > real (kind=dp), intent(out) :: particle_mass > real (kind=dp), intent(out), dimension(:, :), allocatable :: > positions, > velocities > ... > end subroutine > end module > > I tested this with a fortran driver program and it looks good, but when I > try with f2py, it cannot compile. It throws the error "Error: Actual > argument for 'positions' must be ALLOCATABLE at (1)". I figure this has > something to do with the auto-generated "*-f2pywrappers2.f90" file, but > the > build deletes the file. > > If anyone knows an f2py friendly way to rework this, I would be happy to > try. I'm also fine with using ctypes if it can handle this case. Can you split your code so that you have a subroutine which calculates the number of particles only, and call this subroutine from your 'original' routine? If yes, then you might be able to say somehting like real (kind=dp), intent(out), dimension(get_particle_number(WHATEVER) :: Not entirely sure, though ... Can someone with more f2py knowledge confirm this? ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] f2py with allocatable arrays
Can you interface your fortran program twice? First time return the number of particles, dimensions etc to python python then creates work array of right size Second interface pass work array as in/out array, dimension in fortran argument list, to fortran fortran copies allocatable arrays to argument arrays clumsy, I know. George Nurser On 3 July 2012 02:17, Casey W. Stark wrote: > Hi numpy. > > Does anyone know if f2py supports allocatable arrays, allocated inside > fortran subroutines? The old f2py docs seem to indicate that the allocatable > array must be created with numpy, and dropped in the module. Here's more > background to explain... > > I have a fortran subroutine that returns allocatable positions and > velocities arrays. I wish I could get rid of the allocatable part, but you > don't know how many particles it will create until the subroutine does some > work (it checks if each particle it perturbs ends up in the domain). > > module zp > implicit none > contains > subroutine ics(..., num_particles, particle_mass, positions, velocities) > use data_types, only : dp > implicit none > ... inputs ... > integer, intent(out) :: num_particles > real (kind=dp), intent(out) :: particle_mass > real (kind=dp), intent(out), dimension(:, :), allocatable :: positions, > velocities > ... > end subroutine > end module > > I tested this with a fortran driver program and it looks good, but when I > try with f2py, it cannot compile. It throws the error "Error: Actual > argument for 'positions' must be ALLOCATABLE at (1)". I figure this has > something to do with the auto-generated "*-f2pywrappers2.f90" file, but the > build deletes the file. > > If anyone knows an f2py friendly way to rework this, I would be happy to > try. I'm also fine with using ctypes if it can handle this case. > > Best, > Casey > > ___ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion > ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Change in memmap behaviour
On Mon, Jul 2, 2012 at 11:52 PM, Sveinung Gundersen wrote:
>
> On 2. juli 2012, at 22.40, Nathaniel Smith wrote:
>
>> On Mon, Jul 2, 2012 at 6:54 PM, Sveinung Gundersen
>> wrote:
>>> [snip]
>>>
>>>
>>>
>>> Your actual memory usage may not have increased as much as you think,
>>> since memmap objects don't necessarily take much memory -- it sounds
>>> like you're leaking virtual memory, but your resident set size
>>> shouldn't go up as much.
>>>
>>>
>>> As I understand it, memmap objects retain the contents of the memmap in
>>> memory after it has been read the first time (in a lazy manner). Thus, when
>>> reading a slice of a 24GB file, only that part recides in memory. Our system
>>> reads a slice of a memmap, calculates something (say, the sum), and then
>>> deletes the memmap. It then loops through this for consequitive slices,
>>> retaining a low memory usage. Consider the following code:
>>>
>>> import numpy as np
>>> res = []
>>> vecLen = 3095677412
>>> for i in xrange(vecLen/10**8+1):
>>> x = i * 10**8
>>> y = min((i+1) * 10**8, vecLen)
>>> res.append(np.memmap('val.float64', dtype='float64')[x:y].sum())
>>>
>>> The memory usage of this code on a 24GB file (one value for each nucleotide
>>> in the human DNA!) is 23g resident memory after the loop is finished (not
>>> 24g for some reason..).
>>>
>>> Running the same code on 1.5.1rc1 gives a resident memory of 23m after the
>>> loop.
>>
>> Your memory measurement tools are misleading you. The same memory is
>> resident in both cases, just in one case your tools say it is
>> operating system disk cache (and not attributed to your app), and in
>> the other case that same memory, treated in the same way by the OS, is
>> shown as part of your app's resident memory. Virtual memory is
>> confusing...
>
> But the crucial difference is perhaps that the disk cache can be cleared by
> the OS if needed, but not the application memory in the same way, which must
> be swapped to disk? Or am I still confused?
>
> (snip)
>
>>>
>>> Great! Any idea on whether such a patch may be included in 1.7?
>>
>> Not really, if I or you or someone else gets inspired to take the time
>> to write a patch soon then it will be, otherwise not...
>>
>> -N
>
> I have now tried to add a patch, in the way you proposed, but I may have
> gotten it wrong..
>
> http://projects.scipy.org/numpy/ticket/2179
I put this in a github repo, and added tests (author credit to Sveinung)
https://github.com/thouis/numpy/tree/mmap_children
I'm not sure which branch to issue a PR request against, though.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
