headers
Chris Lee | 1 Apr 14:41
Picon
Picon

cross correlation functions

Hi All,

I have two arrays of data and I want to perform a cross correlation on them. Will correlate(input, weights, output=None, mode='reflect', cval=0.0, origin=0) do this for me?

Would the structure be something like this (pseudo code):
arr1 = [bunch of data]
arr2 = [bunch of data]

input_arr = append(arr1, arr2)
weigths = zeros(input_arr.shape, float)
weigth[max(arr1.shape)::] = 1
out_mat = correlate(input_arr, weights)

Otherwise, I may have to loop... the prospect is not appealing as these arrays are reasonably large.

Cheers
Chris
***************************************************
Chris Lee
Laser Physics and Nonlinear Optics Group
MESA+ Research Institute for Nanotechnology
University of Twente
Phone: ++31 (0)53 489 3968
fax: ++31 (0)53 489 1102
***************************************************



_______________________________________________
SciPy-user mailing list
SciPy-user <at> scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
Permalink | Reply | 3 comments |
headers
Zachary Pincus | 1 Apr 17:30
Picon
Favicon

Re: How to fit a surface from a list of measured 3D points ?

Hello,

> I would like to estimate some geoemetric characteristic of this
> surface and view the variation of radius of curvature.
>
>> What's your eventual goal for the fit surface? There are a lot of
>> possible approaches possible with the tools in scipy (e.g. fit a
>> spline, as you have done, or fit a simpler parameteric surface in a
>> least-squares fashion, or go for a nonparametric approach and
>> interpolate between the points with radial basis functions).  The  
>> best
>> approach will likely depend on what you want to do with the  
>> surface...
>> but what you have seems like a reasonable start.
> I think a 2D surface spline could do the trick, but  I think I didn't
> fit it well :
> - the estimated z can go far below the range of my inputs ( [48-62] )
> - As you can see on small_data.jpg, the surface is really 3D, there is
> a kind of crease in my surface, and there are two z value for (x,y) ~
> (2,2). So I 'm not sure if this is the good approach

Hmm, good point. Can you rotate the data points in the 3D space so  
that the new z values do become a proper function in two dimensions?  
If not, then you'll have to:
a) fit a surface to all of the data in 3D (something done a lot by  
computer graphics and robotics people, who get point clouds as return  
data from LIDAR scanners and similar, and then try to fit the points  
to 3D surfaces for visualization / navigation)

b) Find locally-smooth patches and fit surfaces to these individually  
(the manifold-learning folks do this, e.g. "Hessian LLE"). Say you're  
interested in curvature around a given data point (x, y, z)... you  
could take the points within some neighborhood and then either fit  
them to a simple 3d surface (like some kind of paraboloid), or figure  
out (with e.g. PCA) the best projection of those data points to a  
plane, and then fit a surface to f(x, y) -> z for the transformed data.

or perhaps even c) just calculate what you need from the data points  
directly. If you just need very local curvature data, you could  
probably calculate that from a point and its nearest neighbors. (This  
is really just a degenerate case of b...)

Lots of tools for these tasks are in scipy, but nothing off-the-shelf  
that I know if.

Zach
Permalink | Reply | 3 comments |
headers
Jonathan Guyer | 1 Apr 19:36
Favicon

"smooth" plots

Lets say I want to plot y = f(x).

f() is non-analytical and expensive.

Are there any facilities in SciPy (or elsewhere in Python) for  
automating the selection of values of x that will produce a "smooth"  
plot of y vs. x? The only mentions I've seen for adaptive step size  
selection are for the ode integrators, which I'm not doing. Likewise,  
the splines and other smoothing functions don't seem appropriate  
because I don't have noisy data; I'm trying to figure out what data I  
even need to calculate (which will probably not be noisy at all).

Years ago I did things like this in Mathematica, which was really good  
at it (even though the actual plots were hideous (necessitating export  
of the plot points (which Mma was also really bad at))).
Permalink | Reply | 4 comments |
headers
LB | 1 Apr 19:40
Picon

Re: How to fit a surface from a list of measured 3D points ?

> Hmm, good point. Can you rotate the data points in the 3D space so  
> that the new z values do become a proper function in two dimensions?  
It may be possible, with some manual transformation of the data
points, but I would prefer a more generic approach if possible.

> If not, then you'll have to:
> a) fit a surface to all of the data in 3D (something done a lot by  
> computer graphics and robotics people, who get point clouds as return  
> data from LIDAR scanners and similar, and then try to fit the points  
> to 3D surfaces for visualization / navigation)
>
> b) Find locally-smooth patches and fit surfaces to these individually  
> (the manifold-learning folks do this, e.g. "Hessian LLE"). Say you're  
> interested in curvature around a given data point (x, y, z)... you  
> could take the points within some neighborhood and then either fit  
> them to a simple 3d surface (like some kind of paraboloid), or figure  
> out (with e.g. PCA) the best projection of those data points to a  
> plane, and then fit a surface to f(x, y) -> z for the transformed data.
>
> or perhaps even c) just calculate what you need from the data points  
> directly. If you just need very local curvature data, you could  
> probably calculate that from a point and its nearest neighbors. (This  
> is really just a degenerate case of b...)
>
> Lots of tools for these tasks are in scipy, but nothing off-the-shelf  
> that I know if.

The method c) seems the simplest at first sight but I see two issues
for this local approach :
  - the measured data are noisy. Using the nearest neighbor could give
a noisy result two, especially when looking at a radius of curvature
  - I don't see how to use this approach to plot the variation of
radius of curvature along the surface, It can give me an array of
radius of curvature, but as my data are not regularly spaced, it won't
be easy to handle.

Th method b) seems very fuzzy to me : I don't have any knowledge in
manifold-learning and I would have the second issue of the method c)
too.

The method a) is what I had initially in mind, but I didn't see how to
do this in scipy :-(
I believed at first that I could make a sort of parametric bispline
fit with the functions available in scipy.interpolate, but I didn't
succeed in.
Do you have any example or hint for doing this kind of treatment in
scipy ?

LB
Permalink | Reply | 2 comments |
headers
josef.pktd | 1 Apr 19:51
Picon

Re: How to fit a surface from a list of measured 3D points ?

On Wed, Apr 1, 2009 at 1:40 PM, LB <berthe.loic <at> gmail.com> wrote:
>> Hmm, good point. Can you rotate the data points in the 3D space so
>> that the new z values do become a proper function in two dimensions?
> It may be possible, with some manual transformation of the data
> points, but I would prefer a more generic approach if possible.
>
>> If not, then you'll have to:
>> a) fit a surface to all of the data in 3D (something done a lot by
>> computer graphics and robotics people, who get point clouds as return
>> data from LIDAR scanners and similar, and then try to fit the points
>> to 3D surfaces for visualization / navigation)
>>
>> b) Find locally-smooth patches and fit surfaces to these individually
>> (the manifold-learning folks do this, e.g. "Hessian LLE"). Say you're
>> interested in curvature around a given data point (x, y, z)... you
>> could take the points within some neighborhood and then either fit
>> them to a simple 3d surface (like some kind of paraboloid), or figure
>> out (with e.g. PCA) the best projection of those data points to a
>> plane, and then fit a surface to f(x, y) -> z for the transformed data.
>>
>> or perhaps even c) just calculate what you need from the data points
>> directly. If you just need very local curvature data, you could
>> probably calculate that from a point and its nearest neighbors. (This
>> is really just a degenerate case of b...)
>>
>> Lots of tools for these tasks are in scipy, but nothing off-the-shelf
>> that I know if.
>
> The method c) seems the simplest at first sight but I see two issues
> for this local approach :
>  - the measured data are noisy. Using the nearest neighbor could give
> a noisy result two, especially when looking at a radius of curvature
>  - I don't see how to use this approach to plot the variation of
> radius of curvature along the surface, It can give me an array of
> radius of curvature, but as my data are not regularly spaced, it won't
> be easy to handle.
>
> Th method b) seems very fuzzy to me : I don't have any knowledge in
> manifold-learning and I would have the second issue of the method c)
> too.
>
> The method a) is what I had initially in mind, but I didn't see how to
> do this in scipy :-(
> I believed at first that I could make a sort of parametric bispline
> fit with the functions available in scipy.interpolate, but I didn't
> succeed in.
> Do you have any example or hint for doing this kind of treatment in
> scipy ?
>
> LB

If you have noisy data, then a kernel ridge regression or fitting a
gaussian process might be more appropriate than just interpolation.

I posted a simple example code for it a while ago, and pymvpa has a
more complete implementation. Around 900 points should still be ok,
since it builds the square distance matrix (900,900)

I don't know anything about your curvature measure, but it should be
worth a try.

Josef
Permalink | Reply | 0 comments |
headers
Robert Kern | 1 Apr 21:22
Picon
Gravatar

Re: "smooth" plots

On Wed, Apr 1, 2009 at 12:36, Jonathan Guyer <guyer <at> nist.gov> wrote:
> Lets say I want to plot y = f(x).
>
> f() is non-analytical and expensive.
>
> Are there any facilities in SciPy (or elsewhere in Python) for
> automating the selection of values of x that will produce a "smooth"
> plot of y vs. x? The only mentions I've seen for adaptive step size
> selection are for the ode integrators, which I'm not doing. Likewise,
> the splines and other smoothing functions don't seem appropriate
> because I don't have noisy data; I'm trying to figure out what data I
> even need to calculate (which will probably not be noisy at all).

Most of the spline fitting in scipy.interpolate is for exact
interpolation and not noisy data. If your function is smooth (C0 and
C1 continuous), and there are no singularities, then the most
straightforward thing to do would be to do a coarse sampling of your
function, spline fit it, and then evaluate the spline finely. This is,
of course, approximate, but probably not too bad.

There is no adaptive sampling in scipy. That would require much more
information about the screen geometry and line thicknesses. We leave
that to the plotting frameworks. You might want to check SAGE for this
functionality; they would often have need for it.

--

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
  -- Umberto Eco
Permalink | Reply | 1 comment |
headers
Jonathan Guyer | 1 Apr 21:29
Favicon

Re: "smooth" plots


On Apr 1, 2009, at 3:22 PM, Robert Kern wrote:

> straightforward thing to do would be to do a coarse sampling of your
> function, spline fit it, and then evaluate the spline finely. This is,
> of course, approximate, but probably not too bad.

That's probably a reasonable approach for my immediate need.

> There is no adaptive sampling in scipy. That would require much more
> information about the screen geometry and line thicknesses.

I fully expected to have to tell it what I thought was "good" or  
"smooth", but I'd hoped to find something that would automate the  
guesses. I've got some papers on it around here someplace; I just was  
hoping to exploit somebody else's labor.

> You might want to check SAGE for this
> functionality; they would often have need for it.

Good idea, thanks.
Permalink | Reply | 0 comments |
headers
josef.pktd | 1 Apr 21:41
Picon

Re: cross correlation functions

On Wed, Apr 1, 2009 at 8:41 AM, Chris Lee <c.j.lee <at> tnw.utwente.nl> wrote:
> Hi All,
> I have two arrays of data and I want to perform a cross correlation on them.
> Will correlate(input, weights, output=None, mode='reflect', cval=0.0,
> origin=0) do this for me?
> Would the structure be something like this (pseudo code):
> arr1 = [bunch of data]
> arr2 = [bunch of data]
> input_arr = append(arr1, arr2)
> weigths = zeros(input_arr.shape, float)
> weigth[max(arr1.shape)::] = 1
> out_mat = correlate(input_arr, weights)
> Otherwise, I may have to loop... the prospect is not appealing as these
> arrays are reasonably large.
> Cheers
> Chris

What's the dimension of your arr1, arr2? there are 3 different
correlate in numpy/scipy, for 1d numpy.correlate and for 2d
scipy.signal.correlate2d might do more directly what you want.

Josef
Permalink | Reply | 2 comments |
headers
Gideon Simpson | 2 Apr 01:32
Picon
Favicon

blitz error

I tried this little bit of code:

def linear_p(qp_hat, qm_hat, vg, kappa, n1, il, ik):

	linp = np.zeros(qp_hat.shape) + 0.j
	expr = "linp = -vg * qp_hat * il + kappa * ik * n1 * qm_hat"
	weave.blitz(expr, check_size=1)
	return linp

and got the following error:

Traceback (most recent call last):
   File "./soliton1_blitz.py", line 185, in <module>
     k1p = filter(k1p) + linop_p(Ep_hat, Em_hat)
   File "./soliton1_blitz.py", line 79, in <lambda>
     vg, kappa, n1, il, ik)
   File "/Users/gideon/Code/CME/CMEops_blitz.py", line 10, in linear_p
     weave.blitz(expr, check_size=1)
   File "/opt/lib/python2.5/site-packages/scipy/weave/blitz_tools.py",  
line 62, in blitz
     **kw)
   File "/opt/lib/python2.5/site-packages/scipy/weave/ 
inline_tools.py", line 462, in compile_function
     verbose=verbose, **kw)
   File "/opt/lib/python2.5/site-packages/scipy/weave/ext_tools.py",  
line 367, in compile
     verbose = verbose, **kw)
   File "/opt/lib/python2.5/site-packages/scipy/weave/build_tools.py",  
line 272, in build_extension
     setup(name = module_name, ext_modules = [ext],verbose=verb)
   File "/opt/lib/python2.5/site-packages/numpy/distutils/core.py",  
line 184, in setup
     return old_setup(**new_attr)
   File "/sw/lib/python2.5/distutils/core.py", line 168, in setup
     raise SystemExit, "error: " + str(msg)
scipy.weave.build_tools.CompileError: error: Command "g++ -fno-strict- 
aliasing -mno-fused-madd -DNDEBUG -g -fwrapv -O3 -Wall -I/opt/lib/ 
python2.5/site-packages/scipy/weave -I/opt/lib/python2.5/site-packages/ 
scipy/weave/scxx -I/opt/lib/python2.5/site-packages/scipy/weave/blitz - 
I/opt/lib/python2.5/site-packages/numpy/core/include -I/sw/include/ 
python2.5 -c /Users/gideon/.python25_compiled/ 
sc_24bab5a447df6354e489345507f46f0d2.cpp -o /var/folders/yL/ 
yLD5tRJiGWa7oyM6mWJUN++++TI/-Tmp-/gideon/python25_intermediate/ 
compiler_eded75b75a17d5bc9afe97cc30cf08c0/Users/ 
gideon/.python25_compiled/sc_24bab5a447df6354e489345507f46f0d2.o"  
failed with exit status 1

Could it be because I'm using complex numbers?

-gideon
Permalink | Reply | 1 comment |
headers
Chris Lee | 2 Apr 08:52
Picon
Picon

Re: cross correlation functions

Thanks,

I think the 2d function is what I want. I have two 1-d arrays and I  
would like a 2-d cross correlation of the two.

Cheers
Chris

On Apr 1, 2009, at 9:41 PM, josef.pktd <at> gmail.com wrote:

> On Wed, Apr 1, 2009 at 8:41 AM, Chris Lee <c.j.lee <at> tnw.utwente.nl>  
> wrote:
>> Hi All,
>> I have two arrays of data and I want to perform a cross correlation  
>> on them.
>> Will correlate(input, weights, output=None, mode='reflect', cval=0.0,
>> origin=0) do this for me?
>> Would the structure be something like this (pseudo code):
>> arr1 = [bunch of data]
>> arr2 = [bunch of data]
>> input_arr = append(arr1, arr2)
>> weigths = zeros(input_arr.shape, float)
>> weigth[max(arr1.shape)::] = 1
>> out_mat = correlate(input_arr, weights)
>> Otherwise, I may have to loop... the prospect is not appealing as  
>> these
>> arrays are reasonably large.
>> Cheers
>> Chris
>
> What's the dimension of your arr1, arr2? there are 3 different
> correlate in numpy/scipy, for 1d numpy.correlate and for 2d
> scipy.signal.correlate2d might do more directly what you want.
>
> Josef
> _______________________________________________
> SciPy-user mailing list
> SciPy-user <at> scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user

***************************************************
Chris Lee
Laser Physics and Nonlinear Optics Group
MESA+ Research Institute for Nanotechnology
University of Twente
Phone: ++31 (0)53 489 3968
fax: ++31 (0)53 489 1102
***************************************************
Permalink | Reply | 1 comment |

Gmane