At present, the Grid package function CoordInterpGrid will not accept a
linspace-generated vector as the x-coordinates, making precise,
coordinate-based interpolation a little problematic.
--Peter
On Thursday, April 24, 2014 6:51:40 AM UTC-7, J Luis wrote:
>
> I also agree that the Matlab behavior is more intuitive, but even it can
> fail for the same reason discussed here. After being bitten by one of such
> case that took me a while to debug, when I need the exact ends I now always
> use linspace()
>
> Quinta-feira, 24 de Abril de 2014 14:28:16 UTC+1, Peter Simon escreveu:
>>
>> *Matlab behavior:*
>>
>> >> format long
>> >> x = 0:pi/100:pi;
>> >> length(x)
>>
>> ans =
>>
>> 101
>>
>> >> x(1) - 0
>>
>> ans =
>>
>> 0
>>
>> >> x(end) - pi
>>
>> ans =
>>
>> 0
>>
>> >> diff([max(diff(x)), min(diff(x))])
>>
>> ans =
>>
>> -4.440892098500626e-16
>>
>>
>> *Julia behavior:*
>>
>> julia> x = 0:pi/100:pi;
>>
>> julia> length(x)
>> 100
>>
>> julia> x[1]-0
>> 0.0
>>
>> julia> x[end]-pi
>> -0.031415926535897754
>>
>> julia> diff([maximum(diff(x)), minimum(diff(x))])
>> 1-element Array{Float64,1}:
>> -4.44089e-16
>>
>>
>> I highlighted the important differences in red. IMO the Matlab behavior
>> is more intuitive. If you choose an increment that "very nearly"
>> (apparently within a few ulp, as stated by Stefan) evenly divides the
>> difference between the first and last entries in the colon notation, Matlab
>> assumes that you want the resulting vector to begin and end with exactly
>> these values. This makes sense in most situations, IMO. E.g. when using a
>> range to specify sample points in a numerical integration scheme, or when
>> using a range to specify x-coordinates in an interpolation scheme. In both
>> of these cases you don't want the last point to be adjusted away from what
>> you specified using the colon notation.
>>
>> --Pete
>>
>> --Peter
>>
>> On Thursday, April 24, 2014 5:28:19 AM UTC-7, andrew cooke wrote:
>>>
>>>
>>> how does matlab differ from julia here? (i though the original problem
>>> was related to "fundamental" issues of how you represent numbers on a
>>> computer, not language specific).
>>>
>>>
>>> On Thursday, 24 April 2014 02:24:59 UTC-3, Peter Simon wrote:
>>>>
>>>> Thanks, Simon, that construct works nicely to solve the problem I
>>>> posed.
>>>>
>>>> I have to say, though, that I find Matlab's colon range behavior more
>>>> intuitive and generally useful, even if it isn't as "exact" as Julia's.
>>>>
>>>> --Peter
>>>>
>>>> On Wednesday, April 23, 2014 7:17:23 PM UTC-7, Simon Kornblith wrote:
>>>>>
>>>>> pi*(0:0.01:1) or similar should work.
>>>>>
>>>>> On Wednesday, April 23, 2014 7:12:58 PM UTC-4, Peter Simon wrote:
>>>>>>
>>>>>> Thanks for the explanation--it makes sense now. This question arose
>>>>>> for me because of the example presented in
>>>>>> https://groups.google.com/d/msg/julia-users/CNYaDUYog8w/QH9L_Q9Su9YJ:
>>>>>>
>>>>>> x = [0:0.01:pi]
>>>>>>
>>>>>> used as the set of x-coordinates for sampling a function to be
>>>>>> integrated (ideally over the interval (0,pi)). But the range defined in
>>>>>> x
>>>>>> has a last entry of 3.14, which will contribute to the inaccuracy of the
>>>>>> integral being sought in that example. As an exercise, I was trying to
>>>>>> implement the interpolation solution described later in that thread by
>>>>>> Cameron McBride: "BTW, another possibility is to use a spline
>>>>>> interpolation on the original data and integrate the spline evaluation
>>>>>> with quadgk()". It seems that one cannot use e.g. linspace(0,pi,200)
>>>>>> for
>>>>>> the x values, because CoordInterpGrid will not accept an array as its
>>>>>> first
>>>>>> argument, so you have to use a range object. But the range object has a
>>>>>> built-in error for the last point because of the present issue. Any
>>>>>> suggestions?
>>>>>>
>>>>>> Thanks,
>>>>>>
>>>>>> --Peter
>>>>>>
>>>>>> On Wednesday, April 23, 2014 3:24:10 PM UTC-7, Stefan Karpinski wrote:
>>>>>>>
>>>>>>> The issue is that float(pi) < 100*(pi/100). The fact that pi is not
>>>>>>> rational – or rather, that float64(pi) cannot be expressed as the
>>>>>>> division
>>>>>>> of two 24-bit integers as a 64-bit float – prevents rational lifting of
>>>>>>> the
>>>>>>> range from kicking in. I worried about this kind of issue when I was
>>>>>>> working on FloatRanges, but I'm not sure what you can really do, aside
>>>>>>> from
>>>>>>> hacks where you just decide that things are "close enough" based on
>>>>>>> some ad
>>>>>>> hoc notion of close enough (Matlab uses 3 ulps). For example, you can't
>>>>>>> notice that pi/(pi/100) is an integer – because it isn't:
>>>>>>>
>>>>>>> julia> pi/(pi/100)
>>>>>>> 99.99999999999999
>>>>>>>
>>>>>>>
>>>>>>> One approach is to try to find a real value x such that
>>>>>>> float64(x/100) == float64(pi)/100 and float64(x) == float64(pi). If any
>>>>>>> such value exists, it makes sense to do a lifted FloatRange instead of
>>>>>>> the
>>>>>>> default naive stepping seen here. In this case there obviously exists
>>>>>>> such
>>>>>>> a real number – π itself is one such value. However, that doesn't quite
>>>>>>> solve the problem since many such values exist and they don't
>>>>>>> necessarily
>>>>>>> all produce the same range values – which one should be used? In this
>>>>>>> case,
>>>>>>> π is a good guess, but only because we know that's a special and
>>>>>>> important
>>>>>>> number. Adding in ad hoc special values isn't really satisfying or
>>>>>>> acceptable. It would be nice to give the right behavior in cases where
>>>>>>> there is only one possible range that could have been intended (despite
>>>>>>> there being many values of x), but I haven't figured out how determine
>>>>>>> if
>>>>>>> that is the case or not. The current code handles the relatively
>>>>>>> straightforward case where the start, step and stop values are all
>>>>>>> rational.
>>>>>>>
>>>>>>>
>>>>>>> On Wed, Apr 23, 2014 at 5:59 PM, Peter Simon <[email protected]>wrote:
>>>>>>>
>>>>>>>> The first three results below are what I expected. The fourth
>>>>>>>> result surprised me:
>>>>>>>>
>>>>>>>> julia> (0:pi:pi)[end]
>>>>>>>> 3.141592653589793
>>>>>>>>
>>>>>>>> julia> (0:pi/2:pi)[end]
>>>>>>>> 3.141592653589793
>>>>>>>>
>>>>>>>> julia> (0:pi/3:pi)[end]
>>>>>>>> 3.141592653589793
>>>>>>>>
>>>>>>>> julia> (0:pi/100:pi)[end]
>>>>>>>> 3.1101767270538954
>>>>>>>>
>>>>>>>> Is this behavior correct?
>>>>>>>>
>>>>>>>> Version info:
>>>>>>>> julia> versioninfo()
>>>>>>>> Julia Version 0.3.0-prerelease+2703
>>>>>>>> Commit 942ae42* (2014-04-22 18:57 UTC)
>>>>>>>> Platform Info:
>>>>>>>> System: Windows (x86_64-w64-mingw32)
>>>>>>>> CPU: Intel(R) Core(TM) i7 CPU 860 @ 2.80GHz
>>>>>>>> WORD_SIZE: 64
>>>>>>>> BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY)
>>>>>>>> LAPACK: libopenblas
>>>>>>>> LIBM: libopenlibm
>>>>>>>>
>>>>>>>>
>>>>>>>> --Peter
>>>>>>>>
>>>>>>>>
>>>>>>>
