[Numpy-discussion] matlab to numpy help
I found a matlab script that I want to translate into numpy, but have difficulties with understanding indexing in matlab. I haven't used matlab very much and I was hoping that someone could help with the following: It says: Uns = ones(1,m); ... and then later Xijm(:,Uns); m is a positive integer. The problem is that Xijm should be a 1-dimensional array (I think), so what is Uns doing in the second statement? Is this some way to expand Xijm into a second dimension? What would be a numpy equivalent? Gabriel ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] ising model: f2py vs cython comparison
On Sun, December 23, 2007 3:29 am, Ondrej Certik wrote:
> Hi,
>
> I need to write 2D Ising model simulation into my school, so I wrote
> it in Python, then rewrote it in Fortran + f2py, and also Cython:
>
> http://hg.sharesource.org/isingmodel/
>
> And Cython solution is 2x faster than f2py. I understand, that I am
> comparing many things - wrappers, my fortran coding skills
> vs Cython C code generation and gcc vs gfortran. But still any ideas
> how to speed up fortran+f2py solution is very welcomed.
Though the problem is 2D, your implementations are essentially
1D. If you would treat the array A as 2D array (and avoid calling
subroutine p) then you would gain some 7% speed up in Fortran.
When using -DF2PY_REPORT_ATEXIT for f2py then a summary
of timings will be printed out about how much time was spent
in Fortran code and how much in the interface. In the given case
I get (nsteps=5):
Overall time spent in ...
(a) wrapped (Fortran/C) functions : 1962 msec
(b) f2py interface, 60 calls :0 msec
(c) call-back (Python) functions:0 msec
(d) f2py call-back interface, 0 calls :0 msec
(e) wrapped (Fortran/C) functions (acctual) : 1962 msec
that is, most of the time is spent in Fortran function and no time
in wrapper. The conclusion is that the cause of the
difference in timings is not in f2py or cpython generated
interfaces but in Fortran and C codes and/or compilers.
Some idiom used in Fortran code is just slower than in C..
For example, in C code you are doing calculations using
float precision but in Fortran you are forcing double precision.
HTH,
Pearu
PS: Here follows a setup.py file that I used to build the
extension modules instead of the Makefile:
#file: setup.py
def configuration(parent_package='',top_path=None):
from numpy.distutils.misc_util import Configuration
config = Configuration('',parent_package,top_path)
config.add_extension('mcising', sources=['mcising.f'],
define_macros = [('F2PY_REPORT_ATEXIT',1)]
)
#config.add_extension('pyising', sources=['pyising.pyx'])
return config
from numpy.distutils.core import setup
setup(configuration = configuration)
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] ising model: f2py vs cython comparison
On Sun, December 23, 2007 3:29 am, Ondrej Certik wrote: > Hi, > > I need to write 2D Ising model simulation into my school, so I wrote > it in Python, then rewrote it in Fortran + f2py, and also Cython: > > http://hg.sharesource.org/isingmodel/ > > And Cython solution is 2x faster than f2py. I understand, that I am > comparing many things - wrappers, my fortran coding skills > vs Cython C code generation and gcc vs gfortran. But still any ideas > how to speed up fortran+f2py solution is very welcomed. When using g77 compiler instead of gfortran, I get a speed up 4.8 times. Btw, a line in a if statement of the fortran code should read `A(p(i,j,N)) = - A(p(i,j,N))`. Pearu ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] gfortran/g77+f2py vs gcc+Cython speed comparison
Hi,
I needed to write 2D Ising model simulation into my school and I
decided to compare the two possible solutions how to do it, so I of
course wrote
it in Python, then rewrote it in Fortran + f2py, and also Cython. What
is better? Read below. :) But for the impatient, I am going to use
Cython, reasons below.
CCing to Cython, numpy (f2py is discussed there), and sage-devel
(there are people there who could be interested in these kinds of
comparisons).
The code is available at:
http://hg.sharesource.org/isingmodel/
How to play with that - just do this (after installing Mercurial):
$ hg clone http://hg.sharesource.org/isingmodel/
[...]
$ cd isingmodel
$ hg up db7dd01cdc26 # just to be sure
that we are talking about the same revision / code
$ make
[...]
$ time python simulate.py
[...]
real0m2.026s
user0m1.988s
sys 0m0.020s
This runs Cython code. Then apply this patch to run fortran code instead:
$ hg di
diff -r db7dd01cdc26 simulate.py
--- a/simulate.py Sun Dec 23 02:23:30 2007 +0100
+++ b/simulate.py Sun Dec 23 02:24:33 2007 +0100
@@ -31,8 +31,8 @@ def MC(mu = 1, temp = 2, dim = 20, steps
J=1 #coupling constant
k=1 #Boltzman constant
-#from mcising import mc
-from pyising import mc
+from mcising import mc
+#from pyising import mc
B = D1(A)
mc(B, dim, steps, temp, H, mu, J, k)
return D2(B)
And then again (and apply the patch below, otherwise it might not work for you):
$ time python simulate.py
[...]
real0m3.600s
user0m3.528s
sys 0m0.052s
So it's a lot slower.
We are comparing many things here - wrappers, my fortran coding skills
vs Cython C code generation and gcc vs gfortran. So I wrote to numpy
mailinglist. First Travis (the author of numpy) suggested:
"
My experience with similar kinds of comparisons is that gnu fortran
compilers are not very good, especially on 2-d problems. Try using a
different fortran compiler to see if speeds improve.
"
Then Pearu (the author of f2py) suggested:
"
Though the problem is 2D, your implementations are essentially
1D. If you would treat the array A as 2D array (and avoid calling
subroutine p) then you would gain some 7% speed up in Fortran.
When using -DF2PY_REPORT_ATEXIT for f2py then a summary
of timings will be printed out about how much time was spent
in Fortran code and how much in the interface. In the given case
I get (nsteps=5):
Overall time spent in ...
(a) wrapped (Fortran/C) functions : 1962 msec
(b) f2py interface, 60 calls :0 msec
(c) call-back (Python) functions:0 msec
(d) f2py call-back interface, 0 calls :0 msec
(e) wrapped (Fortran/C) functions (acctual) : 1962 msec
that is, most of the time is spent in Fortran function and no time
in wrapper. The conclusion is that the cause of the
difference in timings is not in f2py or cpython generated
interfaces but in Fortran and C codes and/or compilers.
Some idiom used in Fortran code is just slower than in C..
For example, in C code you are doing calculations using
float precision but in Fortran you are forcing double precision.
HTH,
Pearu
PS: Here follows a setup.py file that I used to build the
extension modules instead of the Makefile:
#file: setup.py
def configuration(parent_package='',top_path=None):
from numpy.distutils.misc_util import Configuration
config = Configuration('',parent_package,top_path)
config.add_extension('mcising', sources=['mcising.f'],
define_macros = [('F2PY_REPORT_ATEXIT',1)]
)
#config.add_extension('pyising', sources=['pyising.pyx'])
return config
from numpy.distutils.core import setup
setup(configuration = configuration)
"
and then quickly added
"
When using g77 compiler instead of gfortran, I get a speed
up 4.8 times.
Btw, a line in a if statement of the fortran code
should read `A(p(i,j,N)) = - A(p(i,j,N))`.
"
(btw I have no idea how it could work for me without the A(p(i,j,N)) =
- A(p(i,j,N)) fix, quite embarassing)
So then we discussed it on #sympy IRC:
* Now talking on #sympy
hi pearu, thanks a lot for testing it!
4.8 speedup, jesus christ. so the gfortran sucks a lot
hey ondrej
* ondrej is trying g77
gortran has advantage of being Fortran 90 compiler
g77 is depricated in Debian and dead upstram (imho)
yes but I guess those guys who are maintaing g77/gfortran,
those do not crunch numbers much;)
g77 has a long development history and is well optimized
I fear the fortran is a bad investment for new projects
I know g77 is well optimized, but it's dead
gfortran is a write from scratch and is too young
do you think many people still write fortran? I use it just
because g77 was faster than gcc
g77 is certainly not dead, scientist use it a lot because of its speed
btw`A(p(i,j,N)) = - A(p(i,j,N))`.
means my fortran code was broken
doesn't it?
some think that fortran is a dead language, some use it a
Re: [Numpy-discussion] matlab to numpy help
Gabriel J.L. Beckers wrote: > I found a matlab script that I want to translate into numpy, but have > difficulties with understanding indexing in matlab. I haven't used > matlab very much and I was hoping that someone could help with the > following: > > It says: > > Uns = ones(1,m); > > ... and then later > > Xijm(:,Uns); > > > m is a positive integer. The problem is that Xijm should be a > 1-dimensional array (I think), so what is Uns doing in the second > statement? Is this some way to expand Xijm into a second dimension? What > would be a numpy equivalent? What this is doing depends on exactly what Xijm is. Matlab started out with 2-D arrays only--*everything* was a 2-D array (matrix) of double-precision numbers. (Even strings were represented as double-precision matrices.) Later, support for more dimensions was tacked on, but a vector in Matlab is still a 2-D array, and by default it is a single row: >> Xijm = [2,3,4] Xijm = 2 3 4 >> size(Xijm) ans = 1 3 >> Xijm(:,ones(1,2)) ans = 2 2 To make a row into a column, you can transpose it: >> XijmT = Xijm.' XijmT = 2 3 4 >> XijmT(:,ones(1,2)) ans = 2 2 3 3 4 4 Numpy is much more general, and supports arrays of any number of dimensions right from the start. My guess is that in your Matlab application, X is a column (XijmT), as in the second example above, and the indexing is repeating the columns as shown. Usually this is done to facilitate an array operation with another array that, in this example, has shape (3,2). In numpy it is not necessary to make a new array with the columns repeated because broadcasting achieves the same result more efficiently: In [1]:import numpy In [2]:XijmT = numpy.array([2,3,4]) In [3]:XijmT.shape Out[3]:(3,) In [4]:A = numpy.arange(6) In [5]:A.shape = 3,2 In [6]:A Out[6]: array([[0, 1], [2, 3], [4, 5]]) In [9]:XijmT[:,numpy.newaxis] + A Out[9]: array([[2, 3], [5, 6], [8, 9]]) The indexing with numpy.newaxis makes a 2-D view of the 1-D array, allowing the column dimension to be broadcast to match that of A. Eric > > Gabriel > > > ___ > Numpy-discussion mailing list > Numpy-discussion@scipy.org > http://projects.scipy.org/mailman/listinfo/numpy-discussion ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] setting items in a matrix
> To learn array basics: > http://pages.physics.cornell.edu/~myers/teaching/ComputationalMethods/python/arrays.html> > http://www.scipy.org/Numpy_Example_List#head-a8a8874581c2ebfc69a37ab513974a229ff3bfaa> > http://www.scipy.org/Numpy_Example_List#head-a261b8dd10bda6a5fc268fe4f4171acee2f83968> > http://homes.esat.kuleuven.be/~python/doc/numpy/array.html> > > To learn matrix basics: > http://www.scipy.org/NumPy_for_Matlab_Users> Thanks you so much! These links are really great and that's what I was really looking for. Nice! Kind regards and nice X-mas, Timmie ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] matlab to numpy help
if m = 5, then Uns = ones(1,m) means 1x5 matrix(vector) of 1: [1 1 1 1 1] Xijm(:,Uns) is called Tony's trick where you actually replicate the Xijm (vector i guess) 5 times (according to the Uns.) e.g. http://xtargets.com/snippets/posts/show/55 Matlab introduced repmat() function that does the job "explicitly", using Tony's trick mostly. Hope it helps, Alex On Dec 23, 7:52 pm, "Gabriel J.L. Beckers" <[EMAIL PROTECTED]> wrote: > I found a matlab script that I want to translate into numpy, but have > difficulties with understanding indexing in matlab. I haven't used > matlab very much and I was hoping that someone could help with the > following: > > It says: > > Uns = ones(1,m); > > ... and then later > > Xijm(:,Uns); > > m is a positive integer. The problem is that Xijm should be a > 1-dimensional array (I think), so what is Uns doing in the second > statement? Is this some way to expand Xijm into a second dimension? What > would be a numpy equivalent? > > Gabriel > > ___ > Numpy-discussion mailing list > [EMAIL PROTECTED]://projects.scipy.org/mailman/listinfo/numpy-discussion ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
