On Fri, Feb 28, 2014 at 4:51 PM, Ozzy Lash <[hidden email]> wrote:
> > > > On Fri, Feb 28, 2014 at 3:16 AM, Juan Pablo Carbajal <[hidden email]> > wrote: >> >> On Fri, Feb 28, 2014 at 3:11 AM, Dmitri A. Sergatskov >> <[hidden email]> wrote: >> > >> > Back to the original problem. >> > I think if we pre-scale argument such that X is about MAXINT, we get the >> > max >> > possible accuracy in this situation. In this particular case >> > >> > >> > >> > Gamma = 1.62e7; >> > duration = 10/Gamma; >> > dt = 0.0025/Gamma; >> > t = 0:dt:duration; >> > sc = 2^31 / duration ; >> > y = mod (sc*t, sc*0.2/Gamma)/sc; >> > >> > [a,b] = find(abs(y) < eps); >> > >> > sum(a) >> > ans = 51 >> > >> > b = >> > >> > Columns 1 through 11: >> > >> > 1 81 161 241 321 401 481 561 641 721 >> > 801 >> > >> > Columns 12 through 22: >> > >> > 881 961 1041 1121 1201 1281 1361 1441 1521 1601 >> > 1681 >> > >> > Columns 23 through 33: >> > >> > 1761 1841 1921 2001 2081 2161 2241 2321 2401 2481 >> > 2561 >> > >> > Columns 34 through 44: >> > >> > 2641 2721 2801 2881 2961 3041 3121 3201 3281 3361 >> > 3441 >> > >> > Columns 45 through 51: >> > >> > 3521 3601 3681 3761 3841 3921 4001 >> > >> > >> > Perhaps we should include such pre-scaling in the code of rem and mod >> > if we try for matlab compatibility... >> > >> > Dmitri. >> > -- >> > >> > >> > _______________________________________________ >> > Help-octave mailing list >> > [hidden email] >> > https://mailman.cae.wisc.edu/listinfo/help-octave >> > >> >> Hi I am opening a bug for mod, flagged as incompatible result. I get >> the following comaring betwene octave 3.8.1 and matlab2008b >> Gamma = 1.62e7; >> duration = 10/Gamma; >> dt = 0.0025/Gamma; >> t = 0:dt:duration; >> y = mod (t, 0.2/Gamma); >> find (y==0,3,'first') >> >> octave >> 1 241 401 >> >> r2008b >> 1 81 161 >> >> Reading the help of mod in matlab it says >> MOD(x,y) is x - n.*y where n = floor(x./y) if y ~= 0. If y is not an >> integer and the quotient x./y is within roundoff error of an integer, >> then n is that integer. >> >> So indeed matlab is giving a result considerig roundoff error, I >> assume they do something like >> function m = mod_ml(x,y) >> if fix(y) != y >> err = abs (x./y - round(x./y)) < sqrt (eps); >> m = mod (x,y); >> m(err) = 0; >> endif >> endfunction >> >> Shall I fix our mod function? >> >> @Renaud: you can use mod_ml for the moment and tell us if it fixes our >> issues >> _______________________________________________ >> Help-octave mailing list >> [hidden email] >> https://mailman.cae.wisc.edu/listinfo/help-octave > > > Juan, > > What is the line: > > m(err) = 0; > > in your mod_ml function trying to do. I tried it in an ancient version of > octave (3.0.5) and it didn't increment n at all, but if I changed the line > to: > > m=0; > > It gave 4001 increments. > > Bill The function mod_ml returns a zero when the relation when x is divisible y, i.e. x is congruent modulo y with zero. m(err) =0 replaces all answers of mod that weren't zero but should have been if one accepts the matlab behavior of defining a tolerance interval of length 2*eps around each integer. _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave |
On Fri, Feb 28, 2014 at 5:07 PM, Juan Pablo Carbajal
<[hidden email]> wrote: > On Fri, Feb 28, 2014 at 4:51 PM, Ozzy Lash <[hidden email]> wrote: >> >> >> >> On Fri, Feb 28, 2014 at 3:16 AM, Juan Pablo Carbajal <[hidden email]> >> wrote: >>> >>> On Fri, Feb 28, 2014 at 3:11 AM, Dmitri A. Sergatskov >>> <[hidden email]> wrote: >>> > >>> > Back to the original problem. >>> > I think if we pre-scale argument such that X is about MAXINT, we get the >>> > max >>> > possible accuracy in this situation. In this particular case >>> > >>> > >>> > >>> > Gamma = 1.62e7; >>> > duration = 10/Gamma; >>> > dt = 0.0025/Gamma; >>> > t = 0:dt:duration; >>> > sc = 2^31 / duration ; >>> > y = mod (sc*t, sc*0.2/Gamma)/sc; >>> > >>> > [a,b] = find(abs(y) < eps); >>> > >>> > sum(a) >>> > ans = 51 >>> > >>> > b = >>> > >>> > Columns 1 through 11: >>> > >>> > 1 81 161 241 321 401 481 561 641 721 >>> > 801 >>> > >>> > Columns 12 through 22: >>> > >>> > 881 961 1041 1121 1201 1281 1361 1441 1521 1601 >>> > 1681 >>> > >>> > Columns 23 through 33: >>> > >>> > 1761 1841 1921 2001 2081 2161 2241 2321 2401 2481 >>> > 2561 >>> > >>> > Columns 34 through 44: >>> > >>> > 2641 2721 2801 2881 2961 3041 3121 3201 3281 3361 >>> > 3441 >>> > >>> > Columns 45 through 51: >>> > >>> > 3521 3601 3681 3761 3841 3921 4001 >>> > >>> > >>> > Perhaps we should include such pre-scaling in the code of rem and mod >>> > if we try for matlab compatibility... >>> > >>> > Dmitri. >>> > -- >>> > >>> > >>> > _______________________________________________ >>> > Help-octave mailing list >>> > [hidden email] >>> > https://mailman.cae.wisc.edu/listinfo/help-octave >>> > >>> >>> Hi I am opening a bug for mod, flagged as incompatible result. I get >>> the following comaring betwene octave 3.8.1 and matlab2008b >>> Gamma = 1.62e7; >>> duration = 10/Gamma; >>> dt = 0.0025/Gamma; >>> t = 0:dt:duration; >>> y = mod (t, 0.2/Gamma); >>> find (y==0,3,'first') >>> >>> octave >>> 1 241 401 >>> >>> r2008b >>> 1 81 161 >>> >>> Reading the help of mod in matlab it says >>> MOD(x,y) is x - n.*y where n = floor(x./y) if y ~= 0. If y is not an >>> integer and the quotient x./y is within roundoff error of an integer, >>> then n is that integer. >>> >>> So indeed matlab is giving a result considerig roundoff error, I >>> assume they do something like >>> function m = mod_ml(x,y) >>> if fix(y) != y >>> err = abs (x./y - round(x./y)) < sqrt (eps); >>> m = mod (x,y); >>> m(err) = 0; >>> endif >>> endfunction >>> >>> Shall I fix our mod function? >>> >>> @Renaud: you can use mod_ml for the moment and tell us if it fixes our >>> issues >>> _______________________________________________ >>> Help-octave mailing list >>> [hidden email] >>> https://mailman.cae.wisc.edu/listinfo/help-octave >> >> >> Juan, >> >> What is the line: >> >> m(err) = 0; >> >> in your mod_ml function trying to do. I tried it in an ancient version of >> octave (3.0.5) and it didn't increment n at all, but if I changed the line >> to: >> >> m=0; >> >> It gave 4001 increments. >> >> Bill > > The function mod_ml returns a zero when the relation when x is > divisible y, i.e. x is congruent modulo y with zero. > m(err) =0 > replaces all answers of mod that weren't zero but should have been if > one accepts the matlab behavior of defining a tolerance interval of > length 2*eps around each integer. sorry, the interval has 2*sqrt(eps) _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave |
Regarding the m(err) =0; I must have cut and pasted wrong the first time. I see what you are trying to do. I was just thinking of m as a scalar (I was passing x and y as scalars as in the original script). It works now. Sorry about that. _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave |
This post was updated on .
In reply to this post by Gordon Haverland
Thank you all for your answers. I will try to answer every question that was asked to me.
@ Gordon : You are right, I am drawing pseudo-random numbers at a certain point in the code. I am not aware of Mersenne-Twister and /dev/random stuff, so that I have simply used the rand() function the simplest way possible, in order to get a pseudo-random number whose distribution is uniform in [0,1]. For the entropy of the system, here is what I got (I am running Linux, no problem of translation under other OS) : [renaud@localhost ~]$ cat /proc/sys/kernel/random/entropy_avail 3464 @ Juan-Pablo : Indeed, I have printed the value taken by "n", value which stops increasing at 27 under octave and goes as expected to 51 under matlab, which I have explained in a previous message. I should have mentioned it again, sorry for that. Concerning the outputs you asked for, they look to be the same (visual constatation). @ Bill : This is exactly what I expect from the code. Your suggestion about the modulo is interesting. However, since your example using a comparison to an epsilon yields different results, have you got another idea on how I could perform the modulo operation ? @ Dmitri and Briankaz : Why calling modulo on non-integer numbers is a bad idea ? Does it come from roundoff errors in the numerical representation of real numbers ? Have you got a suggestion for performing such an operation with no risk of failure ? Basically, I want to record the simulation results only at some multiple integers of the time step dt. @ Dmitri (again) : can you please develop your explanation in the message posted at 1:21 am @ Juan-Pablo : your "mod_ml" function seems to work perfectly (both on the much simplified and older version I posted and on the latest version of the code). I am still waiting before becoming too enthusiastic (I have to perform further checkings tonight) but I think your function fixed the problem :) Now, is the problem I had before related to a bad usage of modulo or related to a bad coded modulo function (particularly the example posted by Dmitri surprised me) ? Thanks to all of you, Renaud |
In reply to this post by Gordon Haverland
On 02/27/2014 05:57 PM, ghaverla wrote:
That could be a possible explanation, and if so, a simple solution would be to use /dev/urandom instead. /dev/urandom does not exhaust the entropy because it just runs a PRNG off the /dev/random seed.From what I read, Octave is supposed to make use of a Mersenne-Twister for its core 0-1 RNG. And I am guessing you are drawing random numbers from somewhere in your calculations. If somehow you are getting random numbers from /dev/random instead of the Mersenne-Twister, you could be depleting your system of entropy, and hence every random you generate ends up being 0. /proc/sys/kernel/random/entropy_avail can help. My description assumes Linux, I don't know how this might translate to other hardware/OS. But I have seen entirely too many project descriptions, such as a playing card game, where they are drawing random numbers from /dev/random (directly or indirectly). _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave |
In reply to this post by rchretien
On Fri, Feb 28, 2014 at 7:33 PM, rchretien <[hidden email]> wrote:
> @ Juan-Pablo : your "mod_ml" function seems to work perfectly (both on the > much simplified and older version I posted and on the latest version of the > code). I wait before becoming too enthusiastic (I have to perform further > checkings tonight) but I think your function fixed the problem :) Now, is > the problem I had before related to a bad usage of modulo or related to a > bad coded modulo function ? It is neither of them. It is a matter of convention. Matlab's version of mod assumes that when x/y is almost an integer then it should be an integer (and therefore mod returns zero). Octave's mod leaves this to the user. so when you call mod with a real valued modulus, you should be aware of round off issues. The advantage of Matlab's convetion is that most user wont worry and wont notice the problem. The advantage of Octave it that it gives you an opportunity to learn, as it might be your case now :D. The disadvantage of Matlab's function it that it might hide other problems that happen to fall under this criteria. The disadvantage of octaves is that user should be aware of what their are doing, which is not always the case (included me!), since one can't be an expert on everything. Some people though may argument in favor of one or another, I think it is almost a philosophical issue. It is probable that due to compatibility with matlab we will change mode to reproduce their criteria. _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave |
In reply to this post by Przemek Klosowski-7
On Fri, 28 Feb 2014 13:46:57 -0500
Przemek Klosowski <[hidden email]> wrote: > On 02/27/2014 05:57 PM, ghaverla wrote: > > From what I read, Octave is supposed to make use of a > > Mersenne-Twister for its core 0-1 RNG. And I am guessing you are > > drawing random numbers from somewhere in your calculations. If > > somehow you are getting random numbers from /dev/random instead of > > the Mersenne-Twister, you could be depleting your system of > > entropy, and hence every random you generate ends up being > > 0. /proc/sys/kernel/random/entropy_avail can help. My description > > assumes Linux, I don't know how this might translate to other > > hardware/OS. > > > > But I have seen entirely too many project descriptions, such as a > > playing card game, where they are drawing random numbers > > from /dev/random (directly or indirectly). > That could be a possible explanation, and if so, a simple solution > would be to use /dev/urandom instead. /dev/urandom does not exhaust > the entropy because it just runs a PRNG off the /dev/random seed. I think what Octave does by default, is what Perl can do optionally, and that is replace the rand() in libc, with a Mersenne-Twister. I would disagree with using /dev/urandom. Switching to a long period RNG is a better solution. But, the OP now knows how to query the system for how much entropy it has stored up. Gord _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave |
In reply to this post by Juan Pablo Carbajal-2
Hi,
Thanks for the explanation, I shall be cautious with those roundoff errors when dealing with a mod called on real valued numbers. The difference of behaviour between the two softwares surprised me, but as you mentioned, the good point is that I was given an opportunity to learn something that could have caused me more trouble in my future researchs. Anyway, I have run the code again on the clusters with your mod_ml function and everything is looking fine (through calculations are not over yet). Thanks to every participant of that discussion for helping me. Best regards, Renaud PS : Should I rename the topic with a more appropriate name ? |
Free forum by Nabble | Edit this page |