gmane.comp.python.scientific.user

griddata not working after update to Python 2.7

Umut Yildiz — 2012-05-24T13:50:52

Dear All,

I used to use griddata in order to make my contourmaps. However, after I updated
my Python from 2.6 to 2.7 griddata is not working anymore.

I tried some workarounds but no success.

The countourmap that I produced before is here.
http://dl.dropbox.com/u/17983476/matplotlib/contour_dT_workingbefore.png

After the Python 2.7 update, it turns to the following.
http://dl.dropbox.com/u/17983476/matplotlib/contour_dT_broken.png

Here is the datafile.
http://dl.dropbox.com/u/17983476/matplotlib/contour_dT.dat

And the associated python script (which is also below).
http://dl.dropbox.com/u/17983476/matplotlib/contour_dT.py

The code that I was using before is here. I had to comment out #import griddata
line because this is the only way that it continues. Is this a bug in griddata,
or if there are new workarounds, I would be glad to know a new method to produce
my contourplots again.

Thanks a lot

----------------------------
#! /usr/bin python

import os
import sys
import math
from math import *
from numpy import *
#import griddata
from pylab import *
from matplotlib.ticker import FormatStrFormatter
params = {'axes.labelsize': 20,
'text.fontsize': 15,
'legend.fontsize': 14,
'xtick.labelsize': 20,
'ytick.labelsize': 20,
'text.usetex': True }
rcParams.update(params)

par1 = []
par2 = []
chis = []

rfile = file('contour_dT.dat','r')
line = rfile.readline()
data = line.split()

while len(data) >1:
par1.append(float(data[0]))
par2.append(float(data[1]))
chis.append(float(data[2]))
line = rfile.readline()
data = line.split()

par1 = array(par1)
par2 = array(par2)
chis = array(chis)

xi = linspace(3.2,7.8,50)
yi = linspace(15,300,50)
zi = griddata(par2,par1,chis,xi,yi)
levels = [0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.5,2,3,4,6,10,12,15,20,25,30,40,50]
CS = contourf(xi,yi,zi,levels,cmap=cm.jet)
CS2 = contour(CS, levels=CS.levels[::2],
                        colors = 'r',
                        hold='on')

cbar = colorbar(CS)
cbar.add_lines(CS2)

savefig("contour_dT.png")
show()

Some numpy funcs for PyPy

Dmitrey — 2012-05-24T11:32:29



hi all,
maybe you're aware of numpypy - numpy port for pypy (pypy.org) - Python
language implementation with dynamic compilation.

Unfortunately, numpypy developmnent is very slow due to strict quality
standards and some other issues, so for my purposes I have provided some
missing numpypy funcs, in particular

* atleast_1d, atleast_2d, hstack, vstack, cumsum, isscalar, asscalar,
asfarray, flatnonzero, tile, zeros_like, ones_like, empty_like,
where, searchsorted

* with "axis" parameter: nan(arg)min, nan(arg)max, all, any

and have got some OpenOpt / FuncDesigner functionality working faster
than in CPython.

File with this functions you can get here

Also you may be interested in some info at http://openopt.org/PyPy
Regards, Dmitrey.
_______________________________________________
SciPy-User mailing list
SciPy-User< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

efficiency of the simplex routine: R (optim) vsscipy.optimize.fmin

servant mathieu — 2012-05-24T08:15:51

Dear scipy users,

Again a question about optimization.

 I've just compared the efficiency of the simplex routine in R (optim) vs
scipy (fmin), when minimizing a chi-square. fmin is faster than optim,
but appears to be less efficient. In R, the value of the function is always
minimized step by step (there are of course some exceptions) while there is
lot of fluctuations in python. Given that the underlying simplex algorithm
is supposed to be the same, which mechanism is responsible for this
difference? Is it possible to constrain fmin so it could be more rigorous?

Cheers,
Mathieu
_______________________________________________
SciPy-User mailing list
SciPy-User< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

Linear algebra problem

federico vaggi — 2012-05-23T17:47:02

Hi,

I have a standard system, of type:

AX = B

Where:

X is a nxm matrix.

B is a mxn matrix.

A is a nxn matrix.

I have B and X, and I am trying to calculate A - however, B, is null matrix.

Using the default numpy solvers (
http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html#numpy.linalg.lstsq)
I obtain the trivial correct solution of A being a null matrix.

What's the most robust way to search for the solution which minimizes the
residuals, while still not returning a null matrix of A?  X, is usually
rank deficient, so I know I won't have an exact solution.

Thanks for all the help,

Federico
_______________________________________________
SciPy-User mailing list
SciPy-User< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

scipy.sparse.linalg.eigs is faster with k=8 than withk=1

Alejandro Weinstein — 2012-05-23T15:33:15

Hi:

I am using scipy.sparse.linalg.eigs with a 336x336 sparse matrix with
1144 nonzero entries. I only need the eigenvector corresponding to the
larger eigenvalue, so I was running the function with k=1. However, I
found that it is about 10 times faster to call the function with k=8.

I am testing this with the following code (available here:
https://gist.github.com/2775892):

#############################################################
import timeit

setup = """
import numpy as np
import scipy.sparse
import scipy.io
from scipy.sparse.linalg import eigs

P = scipy.io.mmread('P.mtx')
"""

n = 10
for k in range(1,20):
    code = 'eigs(P, k=%d)' % k
    t = timeit.timeit(stmt=code, setup=setup, number=n) / n
    print 'k: %2d, time: %5.1f ms' % (k, 1000*t)

#############################################################

The output is

k:  1, time: 301.7 ms
k:  2, time: 242.6 ms
k:  3, time: 352.0 ms
k:  4, time: 168.8 ms
k:  5, time: 148.1 ms
k:  6, time:  93.2 ms
k:  7, time:  70.0 ms
k:  8, time:  29.3 ms
k:  9, time:  45.2 ms
k: 10, time:  63.0 ms
k: 11, time: 209.1 ms
k: 12, time: 170.8 ms
k: 13, time: 120.2 ms
k: 14, time: 104.6 ms
k: 15, time: 115.0 ms
k: 16, time:  97.0 ms
k: 17, time:  94.0 ms
k: 18, time:  94.4 ms
k: 19, time:  74.3 ms

Is this behavior typical for eigs? In other words, should I always use
k set to values around 6 or 8, or is it matrix dependent?

I am also curious about this. I would expect that computing 1
eigenvector should never be slower than computing more eigenvectors.
Any about why is this happening?

Alejandro.

Array Selection Help -Part2-

Lothar Ulferts — 2012-05-23T14:06:43

Dear list,

based upon the thread "Array Selection Help"
http://thread.gmane.org/gmane.comp.python.scientific.user/19412 I would like to
modifiy the task a little bit:
The arr1 contains numbers (label) , arr2 floating points (values). But, here it
is different, zones shouldn't be 'Multiparts'. A zone should be split into its
untouched parts. 
For example:
label:
100 100 100 -99 -99 -99       100 100 100 -99 -99 -99
100 100 -99 -99 200 200       100 100 -99 -99 200 200
-99 -99 -99 -99 200 200  =>   -99 -99 -99 -99 200 200 
300 300 300 300 300 300       300 300 300 300 300 300
200 200 200 -99 100 100       300 400 400 -99 500 500
200 200 200 -99 100 100       400 400 400 -99 500 500

values:
1.5 1.9 1.8 0.3 0.1 0.1
1.5 1.7 0.6 0.3 2.5 2.9
0.6 0.6 0.8 0.4 2.1 2.1
3.1 3.2 3.3 3.4 3.5 3.6
4.7 4.7 4.0 0.1 1.0 1.4
4.3 4.0 4.9 0.3 1.2 1.1

Result of zonal min should be:
1.7 1.7 1.7 -99 -99 -99
1.7 1.7 -99 0.3 2.9 2.9
-99 -99 -99 -99 2.9 2.9
3.6 3.6 3.6 3.6 3.6 3.6
4.0 4.0 4.0 -99 1.0 1.0
4.0 4.0 4.0 -99 1.0 1.0


The question is, how do I tranform the left array to the right?

beste regards
Lothar

Wrong Step Response

otti — 2012-05-23T08:01:31

Hello everyone,

after reading some threads about step response problems in the mailing list
i couldnt come up with a proper solution for my problem.

I've got the following transfer function:

G = Vr/(sT1*(1+sT2))

Vr = 2/15
T1 = 0.2
T2 = 1

For this i like to plot the step response so i made the following script:

*---Begin Script*

# some imports
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
from scipy import signal

# Constants
T1 = 0.2
T2 = 1
Vr = float(2)/15
Tsim = 50

# Denumerator
N0 = np.array([T1*T2 ,T1 , 0])

# numerator
Z0 = Vr

# create tfcn
sys = sp.signal.lti(Z0, N0)

# create the step response
t = np.linspace(0, Tsim, 1000)
u = np.arange(len(t))
u = np.ones_like(u)
yout = sp.signal.lsim2(sys, T=t, U=u)[1]

plt.figure(1)
plt.plot(t, u, t, yout/yout.max())
plt.grid("on")
plt.xlabel("t")
plt.ylabel("h(t)")
plt.title("Sprungantwort")
plt.show()

*---End Script*

This gives me the following response, which i know is not the rigth one 
because i made this one
already in a lab at university with MATLAB but i don't have MATLAB at 
home because i prefer to use Numpy+Scipy+Matplotlib.

wrong step response

Here is the correct response:

correct step response

After several hours of reading an trying different approaches im kind of 
frustrated ^^.

I would be really thankfull for any help so that i can keep on going the 
Scipy track.


Greetings,
William
_______________________________________________
SciPy-User mailing list
SciPy-User< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

scipy.weave on windows vista compiling error

Tobias — 2012-05-21T11:36:13

Hello!

I'm trying to bring scipy.weave to work on Windows Vista with MinGw. 
I use Python 2.7.1, Scipy 0.10.1 and MinGw 4.6.2
My code looks like the following:

 from scipy import weave
 weave.inline("""print('Hello World!');""", [], compiler = 'mingw32-gcc')

I get the following error message:

 DistutilsPlatformError: don't know how to compile C/C++ code on platform 'nt' 
 with 'mingw32-gcc' compiler

With MinGw come a lot of different compiler. If I choose not 'mingw32-gcc' but 
for example 'gcc' weave complains 

 ValueError: invalid version number '4.'

So I also downloaded the olf MinGW 2.9.5 Version. Then I also get the Error 
"don't know how to compile C/C++ code on platform 'nt'"

The PATH-Variables are set and I also made shure, that there is write access to 
the MinGw and Python folders. Searching for this error message gave no 
progress, but maybe you could help me. Would be great!

All the best
Tobi

subclassing ndarray : i want slice to return ndarray,not subclass

pierre puiseux UPPA — 2012-05-21T09:56:10

Hello,
all my question is in the title. More precisely, like in scipy doc i try this :
======================================
import sys
import numpy as np

class ArrayChild(np.ndarray):

    def __new__(cls, array, info=None):
        obj = np.asarray(array).view(cls)
        obj.info = info
        return obj

    def __array_finalize__(self, obj):
        print>>sys.stderr, "__array_finalize__"
        if obj is None: return
        self.info = getattr(obj, 'info', None)

if __name__=='__main__':
    a = np.arange(6)
    a.shape=2,3
    a_child = ArrayChild(a)
    for x in a_child : 
        print>>sys.stderr, x, type(x)

=====================================
and i have the answer :
=====================================
__array_finalize__
__array_finalize__
[0 1 2] <class '__main__.ArrayChild'>
__array_finalize__
[3 4 5] <class '__main__.ArrayChild'>

=====================================
but i want
=====================================
__array_finalize__
__array_finalize__
[0 1 2] <class '__main__.np.ndarray'>
__array_finalize__
[3 4 5] <class '__main__.np.ndarray'>

=====================================
I've tried to redefine : __getslice__  like this, but it does not work.
=====================================
    def __getslice__(self, *args, **kwargs):
        return np.ndarray.__getslice__(self, *args, **kwargs)


=====================================
Some idea ?
Thanks to numpy/scipy team for that great job.
=====================================

Pierre Puiseux
Laboratoire de Mathématiques Appliquées
Université de Pau et des Pays de l'Adour
pierre.puiseux< at >univ-pau.fr
http://www.univ-pau.fr/~puiseux

_______________________________________________
SciPy-User mailing list
SciPy-User< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

SciPy installation troubles on CentOS 6.2

Magician — 2012-05-19T15:09:39

Hi All,


I'm trying to build SciPy from source code,
but I have some troubles.

My environment is below:
I and my colleagues (other users) want to use recent Python,
so I installed Python from sources, and I can't install SciPy
by using yum command.

Now I'm facing ATLAS compiling errors.
Configuration options are "--prefix=/usr/local/atlas-3.8.4 -Fa alg -fPIC".
I tried to build it for several times, and always I got errors as below:

It's very troublesome for me to build ATLAS by myself.
My purpose is just using SciPy on my Python.
Even if it's optimized not so good for my environment, it's OK.

Is there an easy or a sure way to build and install SciPy?


Magician

scipy test error: undefined symbol: ATL_buildinfo

Sergio Rojas — 2012-05-18T22:46:09

Hello guys,

 After struggling for a while, I just finished the installation
 of scipy using the intel mkl library.

 Ungracefully the tests ended with an error regarding
 "undefined symbol: ATL_buildinfo" (see below)

 Can this scipy installation live with this error? or
 Is it prone to throw out wrong computations?

 Also, why the number of ran test is not a fix number
 (running the tests using a scipy/Atlas installation
 finished succesfully indicating: Ran 5832 tests in 422.155s
 OK (KNOWNFAIL=14, SKIP=42), but the scipy/mkl tests ended with
 Ran 5833 tests in 448.386s FAILED (KNOWNFAIL=14, SKIP=34, errors=1))?

 what is the meaning of KNOWNFAIL=14?

 Sergio

 >$ python_mkl
 Python 2.7.2 (default, May 16 2012, 13:06:17)
 [GCC 4.6.1] on linux3
 Type "help", "copyright", "credits" or "license" for more information.
 >>> import scipy
 >>> scipy.show_config()
 umfpack_info:
 NOT AVAILABLE
 lapack_opt_info:
 libraries = ['mkl_lapack95_ilp64', 'mkl_lapack95_lp64', 'mkl_rt', 'mkl_intel
 _lp64', 'mkl_intel_thread', 'mkl_core', 'pthread']
 library_dirs = ['/home/srojas/myPROG/IntelC/composer_xe_2011_sp1.7.256/mkl/l
 ib/intel64']
 define_macros = [('SCIPY_MKL_H', None)]
 include_dirs = ['/home/srojas/myPROG/IntelC/composer_xe_2011_sp1.7.256/mkl/i
 nclude']
 blas_opt_info:
 libraries = ['mkl_rt', 'mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'pt
 hread']
 library_dirs = ['/home/srojas/myPROG/IntelC/composer_xe_2011_sp1.7.256/mkl/l
 ib/intel64']
 define_macros = [('SCIPY_MKL_H', None)]
 include_dirs = ['/home/srojas/myPROG/IntelC/composer_xe_2011_sp1.7.256/mkl/i
 nclude']
 lapack_mkl_info:
 libraries = ['mkl_lapack95_ilp64', 'mkl_lapack95_lp64', 'mkl_rt', 'mkl_intel
 _lp64', 'mkl_intel_thread', 'mkl_core', 'pthread']
 library_dirs = ['/home/srojas/myPROG/IntelC/composer_xe_2011_sp1.7.256/mkl/l
 ib/intel64']
 define_macros = [('SCIPY_MKL_H', None)]
 include_dirs = ['/home/srojas/myPROG/IntelC/composer_xe_2011_sp1.7.256/mkl/i
 nclude']
 blas_mkl_info:
 libraries = ['mkl_rt', 'mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'pt
 hread']
 library_dirs = ['/home/srojas/myPROG/IntelC/composer_xe_2011_sp1.7.256/mkl/l
 ib/intel64']
 define_macros = [('SCIPY_MKL_H', None)]
 include_dirs = ['/home/srojas/myPROG/IntelC/composer_xe_2011_sp1.7.256/mkl/i
 nclude']
 mkl_info:
 libraries = ['mkl_rt', 'mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'pt
 hread']
 library_dirs = ['/home/srojas/myPROG/IntelC/composer_xe_2011_sp1.7.256/mkl/l
 ib/intel64']
 define_macros = [('SCIPY_MKL_H', None)]
 include_dirs = ['/home/srojas/myPROG/IntelC/composer_xe_2011_sp1.7.256/mkl/i
 nclude']
 >>> scipy.test('full', verbose=2)
 Running unit tests for scipy
 NumPy version 1.6.1
 NumPy is installed in /home/srojas/myPROG/Python272GnuMKL/lib/python2.7/site-pac
 kages/numpy
 SciPy version 0.10.1
 SciPy is installed in /home/srojas/myPROG/Python272GnuMKL/lib/python2.7/site-pac
 kages/scipy
 Python version 2.7.2 (default, May 16 2012, 13:06:17) [GCC 4.6.1]
 nose version 1.1.2
 ...
 ... (DELETED OUTPUT)
 ...
 ======================================================================
 ERROR: Failure: ImportError (/home/srojas/myPROG/Python272GnuMKL/lib/python2.7/s
 ite-packages/scipy/linalg/atlas_version.so: undefined symbol: ATL_buildinfo)
 ----------------------------------------------------------------------
 Traceback (most recent call last):
 File "/home/srojas/myPROG/Python272GnuMKL/lib/python2.7/site-packages/nose/loa
 der.py", line 390, in loadTestsFromName
 addr.filename, addr.module)
 File "/home/srojas/myPROG/Python272GnuMKL/lib/python2.7/site-packages/nose/imp
 orter.py", line 39, in importFromPath
 return self.importFromDir(dir_path, fqname)
 File "/home/srojas/myPROG/Python272GnuMKL/lib/python2.7/site-packages/nose/imp
 orter.py", line 86, in importFromDir
 mod = load_module(part_fqname, fh, filename, desc)
 File "/home/srojas/myPROG/Python272GnuMKL/lib/python2.7/site-packages/scipy/li
 nalg/tests/test_atlas_version.py", line 6, in <module>
 import scipy.linalg.atlas_version
 ImportError: /home/srojas/myPROG/Python272GnuMKL/lib/python2.7/site-packages/sci
 py/linalg/atlas_version.so: undefined symbol: ATL_buildinfo

 ----------------------------------------------------------------------
 Ran 5833 tests in 448.386s

 FAILED (KNOWNFAIL=14, SKIP=34, errors=1)
 <nose.result.TextTestResult run=5833 errors=1 failures=0>
_______________________________________________
SciPy-User mailing list
SciPy-User< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

about scipy performance

Sergio Rojas — 2012-05-18T20:13:00

Running the performance test of scipy presented at
 [ http://software.intel.com/en-us/articles/numpy-scipy-with-mkl/ ] I found it strange a
 decrease in the time that takes the computation for a large K value (see for instance the value of
 tm corresponding to k=192 and k=200). Does this behaviour makes sense? Is there any test suite available
 to better check the performance of a scipy installation?

 Sergio
 PD. the output of np.show_config() is shown below the data.

 K TM GFLOPS
 64, 0.2182408, 34.91556
 80, 0.2414728, 39.50755
 96, 0.2986258, 38.37579
 104, 0.3231602, 38.43295
 112, 0.3429056, 39.01948
 120, 0.3645972, 39.33108
 128, 0.4074092, 37.55438
 144, 0.445568, 38.6473
 160, 0.4914636, 38.9449
 176, 0.5274322, 39.92931
 192, 0.6281674, 36.5826
 200, 0.6149654, 38.92902
 208, 0.6355928, 39.17603
 224, 0.673929, 39.79648
 240, 0.7244088, 39.67373
 256, 0.8612222, 35.60057
 384, 1.185954, 38.80421



 atlas_threads_info:
 libraries = ['lapack', 'ptf77blas', 'ptcblas', 'atlas']
 library_dirs = ['/home/srojas/myPROG/LapackLib_gfortran/Atlas64b/lib']
 define_macros = [('ATLAS_INFO', '"\\"3.9.72\\""')]
 language = f77
 include_dirs = ['/home/srojas/myPROG/LapackLib_gfortran/Atlas64b/include']
 blas_opt_info:
 libraries = ['ptf77blas', 'ptcblas', 'atlas']
 library_dirs = ['/home/srojas/myPROG/LapackLib_gfortran/Atlas64b/lib']
 define_macros = [('ATLAS_INFO', '"\\"3.9.72\\""')]
 language = c
 include_dirs = ['/home/srojas/myPROG/LapackLib_gfortran/Atlas64b/include']
 atlas_blas_threads_info:
 libraries = ['ptf77blas', 'ptcblas', 'atlas']
 library_dirs = ['/home/srojas/myPROG/LapackLib_gfortran/Atlas64b/lib']
 define_macros = [('ATLAS_INFO', '"\\"3.9.72\\""')]
 language = c
 include_dirs = ['/home/srojas/myPROG/LapackLib_gfortran/Atlas64b/include']
 lapack_opt_info:
 libraries = ['lapack', 'ptf77blas', 'ptcblas', 'atlas']
 library_dirs = ['/home/srojas/myPROG/LapackLib_gfortran/Atlas64b/lib']
 define_macros = [('ATLAS_INFO', '"\\"3.9.72\\""')]
 language = f77
 include_dirs = ['/home/srojas/myPROG/LapackLib_gfortran/Atlas64b/include']
 lapack_mkl_info:
 NOT AVAILABLE
 blas_mkl_info:
 NOT AVAILABLE
 mkl_info:
 NOT AVAILABLE
_______________________________________________
SciPy-User mailing list
SciPy-User< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

is it possible to constrain thescipy.optimize.curve_fit function?

servant mathieu — 2012-05-16T16:20:27

Dear scipy users,

I'm trying to fit to data a power law of the form :




def func (x, a,b, r):

      return r + a*np.power(x,-b)




I would like to constrain the curve_fit routine to only allow
positive parameter values. How is it possible to do so?



Kind regards,

Mathieu
_______________________________________________
SciPy-User mailing list
SciPy-User< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

ANN: IEP 3.0 - the Interactive Editor for Python

Almar Klein — 2012-05-15T07:55:01

Dear all,

We're pleased to announce version 3.0 of the Interactive Editor for Python.

IEP is a cross-platform Python IDE focused on interactivity and
introspection, which makes it very suitable for scientific computing. Its
practical design is aimed at simplicity and efficiency. IEP is written in
Python 3 and Qt. Binaries are available for Windows, Linux, and Mac.

Website: http://code.google.com/p/iep/
Discussion group: http://groups.google.com/group/iep_
Release notes:
http://code.google.com/p/iep/wiki/Release#version_3.0_(14-05-2012)

Regards,
  Rob, Ludo & Almar
_______________________________________________
SciPy-User mailing list
SciPy-User< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

Trajectory Integration via scipy.integrate or PyDSTool?

Sara Gallian — 2012-05-14T16:59:58

Hi all,
I'm trying to integrate particle trajectories in a finite domain, and I need to make them "reflect" at the boundaries (i.e. reverse the velocity perpendicular to the plane they collided with). I started by using  scipy.integrate.odeint , but since the integration steps are variable, simply checking the position and reversing the velocity won't work!
Can anyone suggest the quickest way to obtain this? would Vode be able to handle this, or should I try to learn to use  PyDSTool?
I'm running pretty late on a deadline, so any suggestion is more than appreciated :)
Thanks!
Sara

Using scipy.io.loadmat to read Matlab containers.Mapobject?

Gerrit Holl — 2012-05-14T14:51:54

Hi,

is it possible to use scipy.io.loadmat to read a Matlab containers.Map object?



In [15]: M = scipy.io.loadmat('/tmp/test.mat')

In [16]: M.keys()
Out[16]: ['__function_workspace__', 'None', '__version__',
'__header__', '__globals__']

In [17]: M["None"]
Out[17]:
MatlabOpaque([ ('cm', 'MCOS', 'containers.Map', [[3707764736L], [2L],
[1L], [1L], [1L], [1L]])],
      dtype=[('s0', '|O8'), ('s1', '|O8'), ('s2', '|O8'), ('arr', '|O8')])
In [15]: M = scipy.io.loadmat('/tmp/test.mat')

In [16]: M.keys()
Out[16]: ['__function_workspace__', 'None', '__version__',
'__header__', '__globals__']

In [17]: M["None"]
Out[17]:
MatlabOpaque([ ('cm', 'MCOS', 'containers.Map', [[3707764736L], [2L],
[1L], [1L], [1L], [1L]])],
      dtype=[('s0', '|O8'), ('s1', '|O8'), ('s2', '|O8'), ('arr', '|O8')])



regards,
Gerrit.

scipy.odr - Goodness of fit and parameter estimation for explicit orthogonal distance regression

Markus Baden — 2012-05-14T03:25:41

Hi list,

Currently, I am trying to fit a quadratic curve to a data set which has
much larger errors in the x than in the y direction. My errors are assumed
to be normally distributed and I want to estimate the confidence interval
of the fitted parameters. I have fitted the data two different ways. 1) I
neglect the x errors and fit the quadratic by minimizing the weighted
residuals (y -f) / sig_y via scipy.optimize.leastsq and 2) I use scipy.odr
to fit the parameters. Both result similar fitted parameters.

Now I am stuck with estimating the confidence intervals on these errors and
I have a couple of questions. For method 1) the reduced chi squared is bad
(much larger then 1), because when neglecting the x errors, none of the
points actually lie within a couple of standard deviations on the line.
However when including the x-errors they all fall nicely onto the line.

Is there a way for me to include the x errors into my minimization and
goodness of fit estimation via the reduced chi-squared? I came across a
remark in the CERN minuit documentation [1], that seemed to approximate the
function by a line over the point and then use this to convert x to y
errors. I also read something like this in numerical recipes [2], but not
for general least squares (only for the case of linear data). Does anyone
have any pointers into that direction?

My second question is related to method 2). Is there a way of accessing the
goodness of fit for ODR, similar to calculating the reduced chi-squared for
a fit that only has errors in the y? Another question along this line
concerns the scipy.odr.ODR.Output.sd_beta attribute. In the docstring it
says it is the standard error of the parameter, does that mean a 1 standard
deviation confidence interval? And how exactly are they calculated. Tried
to look at the source and in the odrpack guide, but unfortunately couldn't
figure that out.

I somehow have the feeling, that my problem should be quite standard
(goodness of fit with both x and y errors), but so far I could not find a
good explanation. Any pointers to textbooks or resources on the web would
be greatly appreciated.

Best regards,

Markus

[1] http://wwwasd.web.cern.ch/wwwasd/cgi-bin/listpawfaqs.pl/27
[2] Numerical Recipes in C, Second Edition, Chapter 15.3 "Straight-Line
Data with Errors in Both Coordinate"
_______________________________________________
SciPy-User mailing list
SciPy-User< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

Weighted KDE

Zachary Pincus — 2012-05-13T17:07:02

Hello all,

A while ago, someone asked on this list about whether it would be simple to modify scipy.stats.kde.gaussian_kde to deal with weighted data:
http://mail.scipy.org/pipermail/scipy-user/2008-November/018578.html

Anne and Robert assured the writer that this was pretty simple (modulo bandwidth selection), though I couldn't find any code that the original author may have generated based on that advice.

I've got a problem that could (perhaps) be solved neatly with weighed KDE, so I'd like to give this a go. I assume that at a minimum, to get basic gaussian_kde.evaluate() functionality:

(1) The covariance calculation would need to be replaced by a weighted-covariance calculation. (Simple enough.)

(2) In evaluate(), the critical part looks like this (and a similar stanza that loops over the points instead):
# if there are more points than data, so loop over data
for i in range(self.n):
    diff = self.dataset[:, i, newaxis] - points
    tdiff = dot(self.inv_cov, diff)
    energy = sum(diff*tdiff,axis=0) / 2.0
    result = result + exp(-energy)

I assume that, further, the 'diff' values ought to be scaled by the weights, too. Is this all that would need to be done? (For the integration and resampling, obviously, there would be a bit of other work...)

Thanks,
Zach

2D phase unwrapping

Gregor Thalhammer — 2012-05-11T19:29:51

Hi all,

I have beend searching for an implementation for phase unwrapping in 2D (and possibly also 3D). I found this old thread in the numpy mailing list

http://mail.scipy.org/pipermail/numpy-discussion/2008-November/038873.html

which mentions the C implementation from GERI: http://www.ljmu.ac.uk/GERI/90202.htm

While searching I found remarks by the authors that these algorithms have been incorporated into scipy, however I am unable to find them in current scipy or numpy. Am I missing something obvious? 
Instead I found this wrapper: https://github.com/pointtonull/pyunwrap
This seems to be based on a wrapper already mentioned in the above mentioned discussion, but the links mentioned there are dead. I added some setup.py, and with some small modifications I managed to compile the extension both on OS X and Windows. 

I would like to see these algorithms included in scipy, and I am willing to work on this. So now my questions: In the old numpy-discussion thread licensing issues are raised, can anybody tell more? On the GERI homepage they distribute their code under a non-commercial-use license, but the authors seem to agree on incorporating their code into scipy. 
Except from this, what else would be required?

Thanks for any advice
Gregor

Simple ndarray dim question?

Erik Kastman — 2012-05-10T21:15:43

Hi all,

Using SciPy and Matlab, I'm having trouble reconstructing an array to match what is given from a matlab cell array loaded using scipy.io.loadmat().

For example, say I create a cell containing a pair of double arrays in matlab and then load it using scipy.io (I'm using SPM to do imaging analyses in conjunction with pynifti and the like)

Matlab

    >> onsets{1} = [0 30 60 90]
    >> onsets{2} = [15 45 75 105]

Python

    >>> import scipy.io as scio
    >>> mat = scio.loadmat('onsets.mat')
    >>> mat['onsets'][0]
    array([[[ 0 30 60 90]], [[ 15  45  75 105]]], dtype=object)
    
    >>> mat['onsets'][0].shape
    
    (2,)

My question is this: **Why does this numpy array have the shape (2,) instead of (2,1,4)**? In real life I'm trying to use Python to parse a logfile and build these onsets cell arrays, so I'd like to be able to build them from scratch.

When I try to build the same array from the printed output, I get a different shape back:

    >>> new_onsets = array([[[ 0, 30, 60, 90]], [[ 15,  45,  75, 105]]], dtype=object)
    array([[[0, 30, 60, 90]],
    
           [[15, 45, 75, 105]]], dtype=object)
    
    >>> new_onsets.shape
    (2,1,4)

Unfortunately, the shape (vectors of doubles in a cell array) is coded in a spec upstream, so I need to be able to get this saved exactly in this format. Of course, it's not a big deal since I could just write the parser in matlab, but it would be nice to figure out what's going on and add a little to my [minuscule] knowledge of numpy.

Thanks in advance for any suggestions,
Erik

cross-posted to stack-overflow: http://stackoverflow.com/questions/10542263

How to assemble large sparse matrices effectively

cp3028 — 2012-05-09T18:36:39

Hello everyone,

I am working on an FEM project using Scipy. Now my problem is, that
the assembly of the sparse matrices is to slow. I compute the
contribution of every element in dense small matrices (one for each
element). For the assembly of the global matrices I loop over all
small dense matrices and set the matrice entries the following way:
...
[i,j] = someList[k][l]
Mglobal[i,j] = Mglobal[i,j] + Mlocal[k,l]
...

Mglobal is a lil_matrice of appropriate size, someList maps the
indexing variables.

Of course this is rather slow and consumes most of the matrice
assembly time. Is there a better way to assemble a large sparse matrix
from many small dense matrices? I tried scipy.weave but it doesn't
seem to work with sparse matrices

Yours,

cm