@Alejandro Jakubi Not if it helps provide an immediate workaround to resolve the issue at hand for the given individual.

@Alejandro Jakubi Not if it helps provide an immediate workaround to resolve the issue at hand for the given individual.

@Alec Mihailovs This is very interesting. I did a quick check on 32bit Linux, with comparable results. I did all tests on local hardisk, under /tmp.

I also tested fread'ind into a 2^26 length Array in the case that the file only had 5 bytes previously written to it in total. The `fread1` command did the right thing, and did not introduce garbage into the entries past the 5th position.

At length 2^26 the timings differed between fread1 and readbytes was about a factor of one hundred. Very impressive. It's not even clear that the result isn't even better than that, since the fread1 timing varies a bit and is so small that it maybe mostly noise.

It was mostly in the range of sizes 2^22 through 2^24 that readbytes started to slow down. It performed fine at, say 2^18.

One thing that I noticed (Linux 32bit) was that at size 2^26 the readbytes command took this amount of time in various Maple releases:

Maple 14:  5.240 sec

Maple 13.01: 5.644 sec

Maple 12.02: 5.112 sec

Maple 11.02: 1.832 sec

Maple 10.06: 1.552 sec

But 0.05 sec is still much better than 1.552 sec.

I used this code below, in a file. The file "aaa" must first be written out to currentdir(), and contain the data, naturally. Eg, writebytes("aaa",[3,4,5,6,7]) for the "short test".

### cut here ### snip 8< ### check.mpl ###

restart:
kernelopts(printbytes=false):

str,st,ba,bu:=
time[real](),time(),kernelopts(bytesalloc),kernelopts(bytesused):
   B:=Array(1..2^26,98,datatype=integer[1]):
time[real]()-str,time()-st,
kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;

str,st,ba,bu:=
time[real](),time(),kernelopts(bytesalloc),kernelopts(bytesused):
   readbytes("aaa",B):
time[real]()-str,time()-st,
kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;

Vector[row](B[1..10]),Vector[row](B[-10..-1]);

restart:
kernelopts(printbytes=false):

fopen1:=define_external('fopen',
    'filename'::string,
    'mode'::string,
    'RETURN'::integer[4],
    'LIB'="libc.so.6"):
fread1:=define_external('fread',
    'ptr'::REF(ARRAY(datatype=integer[1])),
    'size'::integer[4],
    'count'::integer[4],
    'stream'::integer[4],
    'RETURN'::integer[4],
    'LIB'="libc.so.6"):
fclose1:=define_external('fclose',
    'stream'::integer[4],
    'RETURN'::integer[4],
    'LIB'="libc.so.6"):

B:=Array(1..2^26,98,datatype=integer[1]):

str,st,ba,bu:=
time[real](),time(),kernelopts(bytesalloc),kernelopts(bytesused):
   p:=fopen1("aaa","rb"):
   fread1(B,1,2^26,p):
   fclose1(p):
time[real]()-str,time()-st,
kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;

Vector[row](B[1..10]),Vector[row](B[-10..-1]);

### cut here ### snip 8< ### end of file ###

ps. It can be handy to use integer[kernelopts('wordsize')/8] although Windows 64 sometimes makes that problematic.

It's probably worth mentioning that calling readbytes("aaa",B,N) for N the size didn't seem to speed it up.

So, now I wonder how fast endianness can be accomodated (in a C wrapper, obviously) for wider formats. And I wonder about ImageTools:-Read and ImportMatrix:-ReadBinaryFile. It's hard to tell, since as I am sure you know iolib(4,..) is a kernel built-in.

acer

reality := realities[reality];

@herclau You don't have to explain to me the use of SVD for Least Squares solving.

Look:

> y.U.Ss.V;

         [0.00520590831446885, 1.00159598094724, 0.999421478032667]

> V^%T.Ss.U^%T.y;

                            [0.00520590831446883]
                            [                   ]
                            [   1.00159598094724]
                            [                   ]
                            [  0.999421478032667]

Those two above are just transpositions of each other. Earlier, I suggested you use y.U.Ss.V. But of course you can use V^%T.Ss.U^%T.y as well, because they are just transpositions of each other! I simply chose the version which required no additional cost of actually doing transpositon.

You constructed Ss as the inverse of DiagonalMatrix(S[1 .. 3]) so of course Ss is symmetric. And so when multiplying y by V^%T.Ss.U^%T you are (least squares) solving M.x=y.

In your P2 you actually used V.U^%T.y.Ss which is no good. I don't know where you get that formula from.

Remember, SingularValues returns Vt which is the transpose of V. Forgetting that is not enough to generate your strange formula used in your P2, I think.

For fun (but perhaps not optimally using the 'thin' option, I didn't check), since Maple uses SVD to compute the non-square float pseudo-inverse,

> MatrixInverse(M).y;

                            [0.00520590831446895]
                            [                   ]
                            [   1.00159598094724]
                            [                   ]
                            [  0.999421478032667]

@herclau You don't have to explain to me the use of SVD for Least Squares solving.

Look:

> y.U.Ss.V;

         [0.00520590831446885, 1.00159598094724, 0.999421478032667]

> V^%T.Ss.U^%T.y;

                            [0.00520590831446883]
                            [                   ]
                            [   1.00159598094724]
                            [                   ]
                            [  0.999421478032667]

Those two above are just transpositions of each other. Earlier, I suggested you use y.U.Ss.V. But of course you can use V^%T.Ss.U^%T.y as well, because they are just transpositions of each other! I simply chose the version which required no additional cost of actually doing transpositon.

You constructed Ss as the inverse of DiagonalMatrix(S[1 .. 3]) so of course Ss is symmetric. And so when multiplying y by V^%T.Ss.U^%T you are (least squares) solving M.x=y.

In your P2 you actually used V.U^%T.y.Ss which is no good. I don't know where you get that formula from.

Remember, SingularValues returns Vt which is the transpose of V. Forgetting that is not enough to generate your strange formula used in your P2, I think.

For fun (but perhaps not optimally using the 'thin' option, I didn't check), since Maple uses SVD to compute the non-square float pseudo-inverse,

> MatrixInverse(M).y;

                            [0.00520590831446895]
                            [                   ]
                            [   1.00159598094724]
                            [                   ]
                            [  0.999421478032667]

I wonder whether this is due to some inefficiency related to (lack of decent) buffering. I guess that I mean: it doesn't sound so good to read from disk a byte at a time (even if not importing integer[1] and even if endianness is a concern). I wonder what the performance would be from just using something like Linux fread (buffering according to cache size, say)?

acer

I believe that this has been mentioned before. But a reminder serves its purpose.

p.s.! Welcome back, Alec, whether here for long or not.

acer

@herclau So first one must correct the cos term your P2 and P3, by multiplying the call to cos() by a factor of 0.5 as mentioned. After doing so then one can get a plot that matches your followup plot on the restricted range 4.0 to 5.7, as above. OK. And then you are unsatisfied with these results of your by-hand use of SVD to do the least-squares fitting, demonstrated by how the red curve is not matching well in that restricted plot.

Did you want,

> y.U.Ss.V;

         [0.00520590828213008, 1.00159598094467, 0.999421478049389]

instead of your,

> V.U^%T.y.Ss;

          [0.0103501902515756, 0.967599719533608, 1.05732003209726]

?

Using y.U.Ss.V gets the overlaid plot below, where the red and green curves coincide so much that only the green curve (on top) is easily visible,

acer

@herclau So first one must correct the cos term your P2 and P3, by multiplying the call to cos() by a factor of 0.5 as mentioned. After doing so then one can get a plot that matches your followup plot on the restricted range 4.0 to 5.7, as above. OK. And then you are unsatisfied with these results of your by-hand use of SVD to do the least-squares fitting, demonstrated by how the red curve is not matching well in that restricted plot.

Did you want,

> y.U.Ss.V;

         [0.00520590828213008, 1.00159598094467, 0.999421478049389]

instead of your,

> V.U^%T.y.Ss;

          [0.0103501902515756, 0.967599719533608, 1.05732003209726]

?

Using y.U.Ss.V gets the overlaid plot below, where the red and green curves coincide so much that only the green curve (on top) is easily visible,

acer

@hirnyk ...that doesn't really explain how the smooth overlaid curve in the parent post was obtained in Mathematica. Which is presumably why the Answer above already linked to that help page you snapped. What value of Mma's ListPlot's InterpolationOrder would fill in such a good a match for sin+cos/2 by mere splines, in the tight corners where there are few data points?

While on topic, the original post had another funny thing, at least to my eye. The 'errors' introduced by the addition of that final random array, with points between 0.0 and 0.1, look a little small on the plot. it doesn't look like the points vary from the sin-cos/2 function by that much (ie, not by 0.05 in y direction on average). Maybe my eyes are going. (Note, the Answer above uses 0.0..0.01 for those 'errors'. But it is the original Mma's 0.0..0.1 range for them that seems under-represented in the plotted data.)

Or maybe the original plot was made with some (other) command variant(s).

@hirnyk ...that doesn't really explain how the smooth overlaid curve in the parent post was obtained in Mathematica. Which is presumably why the Answer above already linked to that help page you snapped. What value of Mma's ListPlot's InterpolationOrder would fill in such a good a match for sin+cos/2 by mere splines, in the tight corners where there are few data points?

While on topic, the original post had another funny thing, at least to my eye. The 'errors' introduced by the addition of that final random array, with points between 0.0 and 0.1, look a little small on the plot. it doesn't look like the points vary from the sin-cos/2 function by that much (ie, not by 0.05 in y direction on average). Maybe my eyes are going. (Note, the Answer above uses 0.0..0.01 for those 'errors'. But it is the original Mma's 0.0..0.1 range for them that seems under-represented in the plotted data.)

Or maybe the original plot was made with some (other) command variant(s).

@Maple_Maple I am not sure what you mean by "related at T". The solutions I showed each specified what T (along with all the other variables) had to be.

Do you mean that you think that the only valid answers are ones in which the other variables are defined in terms of T? If that's so, then run,

solve({fff},AllSolutions);

yourself and see what from the whole solution set might satisfy your need.

You could try your luck with something like,

eliminate({fff}, indets({fff},name) minus {constants} minus {T});

but I recall having a problem with that -- some eventual murky error message issued from `eliminate/recursive`.

The command solve({fff},AllSolutions) eventually returns, after about an hour for me in Maple 14, giving even more solutions (involving RootOf, some of them long) but also including those I showed above.

note to self: there are some duplicate entries in the solutions I showed above (differing only in the particular name of the assumed _Z parameters).

@Maple_Maple I am not sure what you mean by "related at T". The solutions I showed each specified what T (along with all the other variables) had to be.

Do you mean that you think that the only valid answers are ones in which the other variables are defined in terms of T? If that's so, then run,

solve({fff},AllSolutions);

yourself and see what from the whole solution set might satisfy your need.

You could try your luck with something like,

eliminate({fff}, indets({fff},name) minus {constants} minus {T});

but I recall having a problem with that -- some eventual murky error message issued from `eliminate/recursive`.

The command solve({fff},AllSolutions) eventually returns, after about an hour for me in Maple 14, giving even more solutions (involving RootOf, some of them long) but also including those I showed above.

note to self: there are some duplicate entries in the solutions I showed above (differing only in the particular name of the assumed _Z parameters).

@epostma Erik, you are a gentleman and a scholar.

An entertaining post!

That financial summary at the end is the best part.

acer