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 <psimo...@gmail.com>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
>>>>>>>
>>>>>>>
>>>>>>
