| Re: How to speed up KRIG2D by 30x [message #86129 is a reply to message #86125] |
Thu, 10 October 2013 04:44   |
![Mike[5] is currently offline](theme/default/images/xoffline.png.pagespeed.ic.XRkd1fkXye.png "Mike[5] is currently offline")
Mike[5]
Messages: 5
Registered: October 2013
|
Junior Member |
|
|
Hi Chris,
Thank you for your prompt response.
You are right, I did notice the covariances were not right. I rewrote them as variograms instead, like below:
FUNCTION Krig_sphere, rad, c ;Return Spherical variogram Fcn
r = rad/c[0] < 1.0
return, c[1] + c[2] * (r*(1.5 - 0.5*r^2))
end
FUNCTION Krig_expon, rad, c ;Return Exponential variogram Fcn
return, c[1] + c[2] * (1.0 - exp(-3.0*rad/c[0]))
end
Another small change I made was to use the two lines below instead of the double loop to fill the array coefficients
a[0,0] = distance_measure( transpose([[x],[y]]), /MATRIX )
a = call_function(fname, a, t) ; Get coefficient matrix
On Wednesday, October 9, 2013 6:32:53 PM UTC+1, Chris Torrence wrote:
> Hi Mike,
>
>
>
> This is fantastic. I'm working on adding in your change. However, I just found a different problem. If you look at the "Krig_Sphere" function within krig2d.pro, the code doesn't match the docs. The documentation states that for spherical covariance:
>
>
>
> C(d) = C1 - 1.5 C1 (d/A) + 0.5 C1 (d/A)^3 if d < A
>
> = C0 + C1 if d = 0
>
> = 0 if d > A
>
>
>
> Note that I threw in a factor of C1 on the first line (the docs are wrong).
>
>
>
> Here is the code:
>
> FUNCTION Krig_sphere, d, t
>
> r = d/t[0]
>
> v = t[1] + t[2] * (r * (1.5 - 0.5 * r *r) > 0)
>
> z = where(d eq 0, count)
>
> if count ne 0 then v[z] = 0
>
> return, (t[1] + t[2]) - v
>
> end
>
>
>
> We are not clipping to 0 for d > A. In fact, for d > A, the function actually starts to go back up and levels off at the constant C1. You can see this with the following plot:
>
> p = plot(krig_sphere(findgen(40), [10, 0.5, 1]))
>
>
>
> I think the code should be something more like:
>
> FUNCTION Krig_sphere, d, t
>
> r = d/t[0]
>
> v = t[2] * (1 - r * (1.5 - 0.5*r*r))
>
> v[WHERE(d eq 0, /NULL)] = t[1] + t[2]
>
> v[WHERE(d gt t[0], /NULL)] = 0
>
> return, v
>
> end
>
>
>
> Thoughts?
>
>
>
> -Chris
>
> IDL Product Lead
>
> ExelisVIS
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
