A18 ==== Here is a way with symbolic lists. (* Here are some sample lists*) lst1={A,B}; lst2={a,b,c}; (*Join into one list*) list1=Join[lst1,lst2]; (*Do the outer product to get all possible ordered pairs*) list2=Partition[Flatten[Outer[List,list1,list1]],2] {{A, A}, {A, B}, {A, a}, {A, b}, {A, c}, {B, A}, {B, B}, {B, a}, {B, b}, {B, c}, {a, A}, {a, B}, {a, a}, {a, b}, {a, c}, {b, A}, {b, B}, {b, a}, {b, b}, {b, c}, {c, A}, {c, B}, {c, a}, {c, b}, {c, c}} (*Turn the pairs into products*) list3=list2/.{x_,y_}¬{x y} {{A^2}, {A*B}, {a*A}, {A*b}, {A*c}, {A*B}, {B^2}, {a*B}, {b*B}, {B*c}, {a*A}, {a*B}, {a^2}, {a*b}, {a*c}, {A*b}, {b*B}, {a*b}, {b^2}, {b*c}, {A*c}, {B*c}, {a*c}, {b*c}, {c^2}} (*Union weeds out repeats*) list4=Union[list3,list3] {{a^2}, {a*A}, {A^2}, {a*b}, {A*b}, {b^2}, {a*B}, {A*B}, {b*B}, {B^2}, {a*c}, {A*c}, {b*c}, {B*c}, {c^2}} (*Now turn the products back into pairs*) (*This is the step the requires symbols*) list5=list4/.{x_^2}¬{x,x}/.{x_*y_}¬{x,y} {{a, a}, {a, A}, {A, A}, {a, b}, {A, b}, {b, b}, {a, B}, {A, B}, {b, B}, {B, B}, {a, c}, {A, c}, {b, c}, {B, c}, {c, c}} (*List it out*) (*If you have numbers, you can now use a Rule to replace*) Sort[lst]//TableForm -- Kevin J. McCann Joint Center for Earth Systems Technology (JCET) Department of Physics UMBC Baltimore MD 21250 > Dear group, I have the following question regarding a lengthy calculation > using Mathematica: For given w points in x direction and h points in y direction, I can > construct all the points using h=10; w=8; > points=Flatten[Transpose[Outer[List,Range[w],Range[h]]],1] Next, I need to find all the possible pairs of point including points > themselves, i.e., pair AA. I can use pairs=Outer[List,points,points,1] Then, I have to clear those pairs that repeat themselves, i.e., pair AB and > pair BA. Also, when w and h are of the order of 1000s, the computation > takes a very long time. Is there a better way of doing the second part of Sincerely Cheng > ==== ================================================ > Cheng Liu, Ph.D. > MST-8, Structure/Property Relations > Materials Science and Technology Division > Los Alamos National Laboratory > Los Alamos, New Mexico 87545 ==== ================================================ ==== Oops! I didn't read the question properly. I hope I have got it right this time. Union[Map[Sort, pairs]] does what you want. -- Steve Luttrell West Malvern, UK > Dear group, I have the following question regarding a lengthy calculation > using Mathematica: For given w points in x direction and h points in y direction, I can > construct all the points using h=10; w=8; > points=Flatten[Transpose[Outer[List,Range[w],Range[h]]],1] Next, I need to find all the possible pairs of point including points > themselves, i.e., pair AA. I can use pairs=Outer[List,points,points,1] Then, I have to clear those pairs that repeat themselves, i.e., pair AB and > pair BA. Also, when w and h are of the order of 1000s, the computation > takes a very long time. Is there a better way of doing the second part of Sincerely Cheng > ==== ================================================ > Cheng Liu, Ph.D. > MST-8, Structure/Property Relations > Materials Science and Technology Division > Los Alamos National Laboratory > Los Alamos, New Mexico 87545 ==== ================================================ ==== Cases[pairs, _?(#[[1]] != #[[2]] &)] does what you want. -- Steve Luttrell West Malvern, UK > Dear group, I have the following question regarding a lengthy calculation > using Mathematica: For given w points in x direction and h points in y direction, I can > construct all the points using h=10; w=8; > points=Flatten[Transpose[Outer[List,Range[w],Range[h]]],1] Next, I need to find all the possible pairs of point including points > themselves, i.e., pair AA. I can use pairs=Outer[List,points,points,1] Then, I have to clear those pairs that repeat themselves, i.e., pair AB and > pair BA. Also, when w and h are of the order of 1000s, the computation > takes a very long time. Is there a better way of doing the second part of Sincerely Cheng > ==== ================================================ > Cheng Liu, Ph.D. > MST-8, Structure/Property Relations > Materials Science and Technology Division > Los Alamos National Laboratory > Los Alamos, New Mexico 87545 ==== ================================================ Reply-To: ==== If we have inaccurately known parameters, I think Interval arithmetic does a far better job of assessing the situation. As for impossible demands on memory and time, the computation that took 992.3 seconds for you took 32.8 seconds for me. Anyway, it can be done faster AND more accurately without bignums: Timing[pts = With[{ss = ser}, Table[({#1, ss} & )[x], {x, 50, 70, 1/10}]]; ] ListPlot[pts, PlotJoined -> True, PlotRange -> All]; MaxMemoryUsed[] {10.640999999999998*Second, Null} In any case, we spent far more time writing code and evaluating results than waiting on execution. If anything, your examples suggest only that machine precision AND bignum computations are suspect. The results may or may not be worth the pixels they take up on my screen, and unless I compute in some alternative way instead -- or use progressively more digits in bignums until things settle down -- I can only guess at their reliability. For an application such as your example, I think the best solution is to use infinite precision for a limited number of points, and then Interpolation. It's safer than using SetPrecision because it doesn't involve guessing how many digits of precision to use, and it's far faster because it doesn't involve testing higher and higher levels of precision. The choice of points for exact computation may be tricky, but there are adaptive algorithms for that. Here's an interesting way to proceed, for instance: ser = Normal[Series[Cos[x], {x, 0, 200}]]; Timing[pts = Table[{x, ser}, {x, 50, 70, 1/2}];] f = Interpolation[pts]; Timing[plot1 = Plot[f[x], {x, 50, 70}, PlotPoints -> 30, PlotDivision -> 3];] Cases[plot1, Line[a__] -> a, Infinity][[1, All, 1]]; Timing[newPts = Union[pts, ({x, ser} /. x -> #) & /@ (Rationalize[#, 1/100] & /@ %)];] g = Interpolation[newPts, InterpolationOrder -> 5]; plot1 = Plot[Cos[x] - g[x], {x, 50, 70}, PlotRange -> All]; {1.703000000000003*Second, Null} {0.1560000000000059*Second, Null} {4.968999999999994*Second, Null} {0.546999999999997*Second, Null} Length[pts] Length[newPts] 41 124 I used only a few points for the first plot and it already looked good. Just to be sure, I used Plot to select more points, and used infinite precision computation again for those points. The final Plot shows error limited to about 10^-6. Increasing InterpolationOrder decreases errors significantly, too, at fairly small cost. Bobby -----Original Message----- >However, if the coefficients and powers of your example series were not > perfectly known, what then? We need to distinguish between having an exact value which is imperfectly known and not having an exact value - having a range of values. If I may speculate a little more on real uses: - If we have inaccurately known parameters that do have definite values we may still want to calculae accurately over possible ranges of the parameter; and if the definite values give distinctive outcomes then testing with high accuracy inputs is a way of getting a more accurate determination of the real value - rather like using an inverse function. - if parameters do not have a definite value then we are into statistics, however we might still need to know the outcomes of inputing accurate values to get an idea of the behaviour of the process. Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 ----- Original Message ----- > adding arbitrary digits serves no purpose. Yes, plots may get smoother > as more digits are added, but they would not converge to a correct > result -- merely to a precise one. (In the chemistry industry where my wife works, the difference between > accuracy and precision is well known. Precision means getting the same > answer over and over --- whether it's right or not. Accuracy means > getting the right answer --- whether it's precise or not. It's low > variance versus small bias.) Modify your example like this: ser = N@Normal[Series[Cos[#], {#, 0, 200}]]; > Timing[pts = With[{ss = > ser}, Table[SetPrecision[{#, ss}, 80] &@x, {x, 50., 70., .1}]];] > ListPlot[pts, PlotJoined -> True, PlotRange -> All]; > MaxMemoryUsed[] Once the series coefficients have lost precision, you can't get it back > again. Furthermore, in using SetPrecision, there's a danger that one > could THINK he has regained it. Bobby -----Original Message----- > Sent: Tuesday, October 15, 2002 3:18 AM > To: mathgroup@smc.vnet.net > Greetings, I have read with great interest this lively debate on numerical > prcesion > and > accuracy. As I work in the fields of finance and economics, where we > feel > ourselves blessed if we get three digits of accuracy, I'm curious as > to > what > scientific endeavors require 50+ digits of precision? As I recall > there > are > some areas, such as high energy physics and some elements of > astronomy, > that > might require so many digits in some circumstances. Are there others? > -Mark > Mark, There may be occasions when the outcome of a real process is so > sensitive > to changes in input that unless we know very precisely what the input is > then we can know very little about the outcome - chaotic processes are > of > this kind. The difficulty is real and no amount of computer power or > clever > progamming will do much about it. Another situation is when the the process is not so sensitive but > calculating with our formula or programme introduces accumulates > significant > errors. Here is a very artificial example of the latter (I time the computation > and > find the MaximumMemory used in the session as we go through the > example): ser=Normal[Series[Cos[#],{#,0,200}]]; MaxMemoryUsed[] 1714248 Calculating with machine number does not show much of a pattern ( I > have > deleted the graphics - please evaluate the code), > pts= With[{ss=ser},Table[ {#,ss}&[x], > {x,50.,70., .1}]];//Timing > ListPlot[pts, PlotJoined->True]; > MaxMemoryUsed[] {5.11 Second,Null} 1723840 Using bigfloat inputs with precision 20 shows some pattern: pts= With[{ss=ser},Table[ {#,ss}&[SetPrecision[x,20]], > {x,50.,70., .1}]];//Timing > ListPlot[pts, PlotJoined->True,PlotRange[Rule]All]; > MaxMemoryUsed[] {17.52 Second,Null} 1759664 > Precision 40 does very well: pts= With[{ss=ser},Table[ {#,ss}&[SetPrecision[x,40]], > {x,50.,70., .1}]];//Timing > ListPlot[pts, PlotJoined->True,PlotRange[Rule]All]; > MaxMemoryUsed[] {19.38 Second,Null} 1797072 Now we might think the correct outcomes are showing up, and use an > interpolating function for further , and faster, calculation. f=Interpolation[pts] InterpolatingFunction[{{50.000000,70.00000}},<>] pts= Table[ f[x],{x,50, 70, .1}];//Timing > ListPlot[pts, PlotJoined->True,PlotRange[Rule]All]; > MaxMemoryUsed[] {0.33 Second,Null} > As a matter of interest, this is what happens if we substitute exact > numbers > (rationals and integers) for reals-- > the computation takes an excessively long time and quite a bit more > memory. pts= With[{ss=ser},Table[ {#,ss}&[SetPrecision[x,Infinity]], > {x,50.,70., .1}]];//Timing > ListPlot[pts, PlotJoined->True,PlotRange[Rule]All]; > MaxMemoryUsed[] {992.28 Second,Null} 2413808 This also shows that we may in fact want to replace exact inputs with > bigfloats. > I should be interested to hear of other example, really real one in > particular. I imagine that there are many situations where trends and > shapes > are more important than specific values. -- > Allan --------------------- > Allan Hayes > Mathematica Training and Consulting > Leicester UK > www.haystack.demon.co.uk > hay@haystack.demon.co.uk > Voice: +44 (0)116 271 4198 > ==== Greetings, Is anyone aware of Mathematica code implementing statistical cluster analys= is, e.g., k-means method, etc.? -Mark ==== Bobby, You rightly point out that care should be exercised when using (high precision) bigfloats, but this should not obscure the proper use of them. I have suggested some uses that are valid subject to circumstances (raising precision) or essential (converting exact numbers to bigfloats to avoid impossible demands on memory and time) - Daniel Lichtblau gave others. >However, if the coefficients and powers of your example series were not > perfectly known, what then? We need to distinguish between having an exact value which is imperfectly known and not having an exact value - having a range of values. If I may speculate a little more on real uses: - If we have inaccurately known parameters that do have definite values we may still want to calculae accurately over possible ranges of the parameter; and if the definite values give distinctive outcomes then testing with high accuracy inputs is a way of getting a more accurate determination of the real value - rather like using an inverse function. - if parameters do not have a definite value then we are into statistics, however we might still need to know the outcomes of inputing accurate values to get an idea of the behaviour of the process. Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 ----- Original Message ----- > result -- merely to a precise one. (In the chemistry industry where my wife works, the difference between > accuracy and precision is well known. Precision means getting the same > answer over and over --- whether it's right or not. Accuracy means > getting the right answer --- whether it's precise or not. It's low > variance versus small bias.) Modify your example like this: ser = N@Normal[Series[Cos[#], {#, 0, 200}]]; > Timing[pts = With[{ss = > ser}, Table[SetPrecision[{#, ss}, 80] &@x, {x, 50., 70., .1}]];] > ListPlot[pts, PlotJoined -> True, PlotRange -> All]; > MaxMemoryUsed[] Once the series coefficients have lost precision, you can't get it back > again. Furthermore, in using SetPrecision, there's a danger that one > could THINK he has regained it. Bobby -----Original Message----- > Sent: Tuesday, October 15, 2002 3:18 AM > To: mathgroup@smc.vnet.net > Greetings, I have read with great interest this lively debate on numerical > prcesion > and > accuracy. As I work in the fields of finance and economics, where we > feel > ourselves blessed if we get three digits of accuracy, I'm curious as > to > what > scientific endeavors require 50+ digits of precision? As I recall > there > are > some areas, such as high energy physics and some elements of > astronomy, > that > might require so many digits in some circumstances. Are there others? > -Mark > Mark, There may be occasions when the outcome of a real process is so > sensitive > to changes in input that unless we know very precisely what the input is > then we can know very little about the outcome - chaotic processes are > of > this kind. The difficulty is real and no amount of computer power or > clever > progamming will do much about it. Another situation is when the the process is not so sensitive but > calculating with our formula or programme introduces accumulates > significant > errors. Here is a very artificial example of the latter (I time the computation > and > find the MaximumMemory used in the session as we go through the > example): ser=Normal[Series[Cos[#],{#,0,200}]]; MaxMemoryUsed[] 1714248 Calculating with machine number does not show much of a pattern ( I > have > deleted the graphics - please evaluate the code), > pts= With[{ss=ser},Table[ {#,ss}&[x], > {x,50.,70., .1}]];//Timing > ListPlot[pts, PlotJoined->True]; > MaxMemoryUsed[] {5.11 Second,Null} 1723840 Using bigfloat inputs with precision 20 shows some pattern: pts= With[{ss=ser},Table[ {#,ss}&[SetPrecision[x,20]], > {x,50.,70., .1}]];//Timing > ListPlot[pts, PlotJoined->True,PlotRange[Rule]All]; > MaxMemoryUsed[] {17.52 Second,Null} 1759664 > Precision 40 does very well: pts= With[{ss=ser},Table[ {#,ss}&[SetPrecision[x,40]], > {x,50.,70., .1}]];//Timing > ListPlot[pts, PlotJoined->True,PlotRange[Rule]All]; > MaxMemoryUsed[] {19.38 Second,Null} 1797072 Now we might think the correct outcomes are showing up, and use an > interpolating function for further , and faster, calculation. f=Interpolation[pts] InterpolatingFunction[{{50.000000,70.00000}},<>] pts= Table[ f[x],{x,50, 70, .1}];//Timing > ListPlot[pts, PlotJoined->True,PlotRange[Rule]All]; > MaxMemoryUsed[] {0.33 Second,Null} > As a matter of interest, this is what happens if we substitute exact > numbers > (rationals and integers) for reals-- > the computation takes an excessively long time and quite a bit more > memory. pts= With[{ss=ser},Table[ {#,ss}&[SetPrecision[x,Infinity]], > {x,50.,70., .1}]];//Timing > ListPlot[pts, PlotJoined->True,PlotRange[Rule]All]; > MaxMemoryUsed[] {992.28 Second,Null} 2413808 This also shows that we may in fact want to replace exact inputs with > bigfloats. > I should be interested to hear of other example, really real one in > particular. I imagine that there are many situations where trends and > shapes > are more important than specific values. -- > Allan --------------------- > Allan Hayes > Mathematica Training and Consulting > Leicester UK > www.haystack.demon.co.uk > hay@haystack.demon.co.uk > Voice: +44 (0)116 271 4198 > Reply-To: ==== You're using SetPrecision when infinite precision is a meaningful option -- when there's no doubt about the coefficients and powers in the series. Bignums clearly make the computation faster in that case. However, if the coefficients and powers of your example series were not perfectly known, what then? If they begin life as machine numbers, adding arbitrary digits serves no purpose. Yes, plots may get smoother as more digits are added, but they would not converge to a correct result -- merely to a precise one. (In the chemistry industry where my wife works, the difference between accuracy and precision is well known. Precision means getting the same answer over and over --- whether it's right or not. Accuracy means getting the right answer --- whether it's precise or not. It's low variance versus small bias.) Modify your example like this: ser = N@Normal[Series[Cos[#], {#, 0, 200}]]; Timing[pts = With[{ss = ser}, Table[SetPrecision[{#, ss}, 80] &@x, {x, 50., 70., .1}]];] ListPlot[pts, PlotJoined -> True, PlotRange -> All]; MaxMemoryUsed[] Once the series coefficients have lost precision, you can't get it back again. Furthermore, in using SetPrecision, there's a danger that one could THINK he has regained it. Bobby -----Original Message----- feel > ourselves blessed if we get three digits of accuracy, I'm curious as to what > scientific endeavors require 50+ digits of precision? As I recall there are > some areas, such as high energy physics and some elements of astronomy, that > might require so many digits in some circumstances. Are there others? > -Mark Mark, There may be occasions when the outcome of a real process is so sensitive to changes in input that unless we know very precisely what the input is then we can know very little about the outcome - chaotic processes are of this kind. The difficulty is real and no amount of computer power or clever progamming will do much about it. Another situation is when the the process is not so sensitive but calculating with our formula or programme introduces accumulates significant errors. Here is a very artificial example of the latter (I time the computation and find the MaximumMemory used in the session as we go through the example): ser=Normal[Series[Cos[#],{#,0,200}]]; MaxMemoryUsed[] 1714248 Calculating with machine number does not show much of a pattern ( I have deleted the graphics - please evaluate the code), pts= With[{ss=ser},Table[ {#,ss}&[x], {x,50.,70., .1}]];//Timing ListPlot[pts, PlotJoined->True]; MaxMemoryUsed[] {5.11 Second,Null} 1723840 Using bigfloat inputs with precision 20 shows some pattern: pts= With[{ss=ser},Table[ {#,ss}&[SetPrecision[x,20]], {x,50.,70., .1}]];//Timing ListPlot[pts, PlotJoined->True,PlotRange[Rule]All]; MaxMemoryUsed[] {17.52 Second,Null} 1759664 Precision 40 does very well: pts= With[{ss=ser},Table[ {#,ss}&[SetPrecision[x,40]], {x,50.,70., .1}]];//Timing ListPlot[pts, PlotJoined->True,PlotRange[Rule]All]; MaxMemoryUsed[] {19.38 Second,Null} 1797072 Now we might think the correct outcomes are showing up, and use an interpolating function for further , and faster, calculation. f=Interpolation[pts] InterpolatingFunction[{{50.000000,70.00000}},<>] pts= Table[ f[x],{x,50, 70, .1}];//Timing ListPlot[pts, PlotJoined->True,PlotRange[Rule]All]; MaxMemoryUsed[] {0.33 Second,Null} As a matter of interest, this is what happens if we substitute exact numbers (rationals and integers) for reals-- the computation takes an excessively long time and quite a bit more memory. pts= With[{ss=ser},Table[ {#,ss}&[SetPrecision[x,Infinity]], {x,50.,70., .1}]];//Timing ListPlot[pts, PlotJoined->True,PlotRange[Rule]All]; MaxMemoryUsed[] {992.28 Second,Null} 2413808 This also shows that we may in fact want to replace exact inputs with bigfloats. I should be interested to hear of other example, really real one in particular. I imagine that there are many situations where trends and shapes are more important than specific values. -- Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 > ==== Here's a very simple solution: v = {{100, 200}, {150, 250}, {120, 270}, {300, 400}}; Interval @@ v List @@ % Interval[{100, 270}, {300, 400}] {{100, 270}, {300, 400}} DrBob -----Original Message----- second list that contains the overall upper and lower edges of the overlapping sub-ranges. A simple example : {{100,200},{150,250},{120,270},{300,400}} would result in {{100,270},{300,400}}. In the real case, the input list has several hundred elements and the output list typically has five elements. I have a working solution based on loops, but there must be a more elegant one. I would be very grateful for any suggestions. John Leary ==== Dear Fellows in MathGroup, I have a list of 17,000+ {x,y} pairs of data each x value is a positive integer from 1 to 100+ each y value is a positive real number As a *short* example, let's consider: data = {{3,1},{4,3},{3,2},{1,10},{4,2},{1,6},{5,2},{2,5},{7,1}} I want to group the data by the x value and report the arithmetic average > of the y values in each group. For the example, i want to report: output = {{1,8},{2,5},{3,1.5},{4,2.5},{5,2},{6,0},{7,1}} In this example, x=6 does not occur so i report the average y[6] = 0. Can anyone suggest a way to do this efficiently?/ many thanks > dave +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ > David E. Burmaster, Ph.D. > Alceon Corporation > POBox 382069 (new Box number effective 1 Sep 2001) > Harvard Square Station > Cambridge, MA 02238-2069 (new ZIP code effective 1 Sep 2001) Voice 617-864-4300 Web http://www.Alceon.com > +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Probably most efficient would be to iterate over the list, bin, and then average the bins. averageByBin[data:{{_,_}..}] := Module[ {len, binsizes, averages}, len = Max[Map[First,data]]; binsizes = Table[0,{len}]; averages = Table[0,{len}]; Map [ (binsizes[[#[[1]]]]++; averages[[#[[1]]]] += #[[2]]) &, data]; Do [If[binsizes[[j]]==0, binsizes[[j]]++], {j,len}]; Transpose[{Range[len],N[averages]/binsizes}] ] In[30]:= averageByBin[data] Out[30]= {{1, 8.}, {2, 5.}, {3, 1.5}, {4, 2.5}, {5, 2.}, {6, 0.}, {7, 1.}} If you separate out integer first values from real second values in the pairs, you can enter two separate lists and take advantage of Compile to make it faster still. Daniel Lichtblau Wolfram Research ==== Dear Fellows in MathGroup, I have a list of 17,000+ {x,y} pairs of data each x value is a positive integer from 1 to 100+ each y value is a positive real number As a *short* example, let's consider: data = {{3,1},{4,3},{3,2},{1,10},{4,2},{1,6},{5,2},{2,5},{7,1}} I want to group the data by the x value and report the arithmetic average of the y values in each group. For the example, i want to report: output = {{1,8},{2,5},{3,1.5},{4,2.5},{5,2},{6,0},{7,1}} In this example, x=6 does not occur so i report the average y[6] = 0. Can anyone suggest a way to do this efficiently?/ many thanks dave +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ David E. Burmaster, Ph.D. Alceon Corporation POBox 382069 (new Box number effective 1 Sep 2001) Harvard Square Station Cambridge, MA 02238-2069 (new ZIP code effective 1 Sep 2001) Voice 617-864-4300 Web http://www.Alceon.com +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ==== One typical idea for your application is to sort the data and then use Split. For example, ({#1[[1,1]], Plus @@ #1[[All,2]]/Length[#1]} & ) /@ Split[Sort[data], #1[[1]] == #2[[1]] & ]; The only difference between the above function and your desired result is that when there is no data for a particular integer, the average for that integer does not appear in the answer. Carl Woll Physics Dept U of Washington > Dear Fellows in MathGroup, I have a list of 17,000+ {x,y} pairs of data each x value is a positive integer from 1 to 100+ each y value is a positive real number As a *short* example, let's consider: data = {{3,1},{4,3},{3,2},{1,10},{4,2},{1,6},{5,2},{2,5},{7,1}} I want to group the data by the x value and report the arithmetic average > of the y values in each group. For the example, i want to report: output = {{1,8},{2,5},{3,1.5},{4,2.5},{5,2},{6,0},{7,1}} In this example, x=6 does not occur so i report the average y[6] = 0. Can anyone suggest a way to do this efficiently?/ many thanks > dave +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ > David E. Burmaster, Ph.D. > Alceon Corporation > POBox 382069 (new Box number effective 1 Sep 2001) > Harvard Square Station > Cambridge, MA 02238-2069 (new ZIP code effective 1 Sep 2001) Voice 617-864-4300 Web http://www.Alceon.com > +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ==== >Dear Fellows in MathGroup, I have a list of 17,000+ {x,y} pairs of data each x value is a positive integer from 1 to 100+ each y value is a positive real number As a *short* example, let's consider: data = {{3,1},{4,3},{3,2},{1,10},{4,2},{1,6},{5,2},{2,5},{7,1}} I want to group the data by the x value and report the arithmetic average >of the y values in each group. For the example, i want to report: output = {{1,8},{2,5},{3,1.5},{4,2.5},{5,2},{6,0},{7,1}} In this example, x=6 does not occur so i report the average y[6] = 0. Can anyone suggest a way to do this efficiently?/ Block[{data = {{3,1},{4,3},{3,2},{1,10},{4,2},{1,6}, {5,2},{2,5},{7,1}}, sd = Table[{0,0.},{10}]}, Off[Infinity::indet]; Off[General::dbyz]; (sd[[#[[1]],2]]++;sd[[#[[1]],1]]+=#[[2]])&/@data; sd=MapIndexed[{#2[[1]],Divide@@#}&,sd]; On[Infinity::indet]; On[General::dbyz]; sd/.Indeterminate->0.] --> {{1, 8.}, {2, 5.}, {3, 1.5}, {4, 2.5}, {5, 2.}, {6, 0.}, {7, 1.}, {8, 0.}, {9, 0.}, {10, 0.}} This makes 3 linear passes overall on the data. For large data sets, that may be a problem. DH ==== > Dear Fellows in MathGroup, I have a list of 17,000+ {x,y} pairs of data each x value is a positive integer from 1 to 100+ each y value is a positive real number As a *short* example, let's consider: data = {{3,1},{4,3},{3,2},{1,10},{4,2},{1,6},{5,2},{2,5},{7,1}} I want to group the data by the x value and report the arithmetic average > of the y values in each group. For the example, i want to report: output = {{1,8},{2,5},{3,1.5},{4,2.5},{5,2},{6,0},{7,1}} In this example, x=6 does not occur so i report the average y[6] = 0. Can anyone suggest a way to do this efficiently?/ many thanks > dave Dave, One way: data={{3,1},{4,3},{3,2},{1,10},{4,2},{1,6},{5,2},{2,5},{7,1}}; f[x_] = 0; ((f[#1[[1,1]]] = Plus @@ #1[[All,2]]/Length[#1]) & ) /@ Split[Sort[data], #1[[1]] == #2[[1]] & ] {8, 5, 3/2, 5/2, 2, 1} Table[f[i], {i, 1, 8}] {8, 5, 3/2, 5/2, 2, 0, 1, 0} -- Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ > David E. Burmaster, Ph.D. > Alceon Corporation > POBox 382069 (new Box number effective 1 Sep 2001) > Harvard Square Station > Cambridge, MA 02238-2069 (new ZIP code effective 1 Sep 2001) Voice 617-864-4300 Web http://www.Alceon.com > +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ==== Perhaps Numerical Recipes in C++ ----- Original Message ----- constants > and addition for example. I currently do this all with Mathematica but for > large symbolic calculations involving many operations Mathematica becomes quite > slow. (at least with my coding ability ;) ) Does anyone have any > recommendations for a c++ book along these lines? Simply put, I basically > want to write c++ code that functions similarly to mathematica but that > is specialized to particlar types of algebraic problems and therefore > runs faster. Any suggestions? -chris ==== I'm looking for a book recommendation. I've been using mathematica and I love it. I'd like to speed up some of my Mathematica codes though. What I really want to do is code some basic algrebraic steps in a lower level language. I need to somehow set up c/c++ code for operators, how they multiply, and also deal with multiplicative constants and addition for example. I currently do this all with Mathematica but for large symbolic calculations involving many operations Mathematica becomes quite slow. (at least with my coding ability ;) ) Does anyone have any recommendations for a c++ book along these lines? Simply put, I basically want to write c++ code that functions similarly to mathematica but that is specialized to particlar types of algebraic problems and therefore runs faster. Any suggestions? -chris ==== > I've been trying to use PlotVectorField for the > following differential equation: > dy/dt = 0.08*y*(1-y/1000) > but I haven't been successful yet. > I tried to do the following: > f[t_, y_] := {1, 0.08*y*(1 - y/1000)} > < PlotVectorField[f[t, y], {t, 0, 80}, {y, 0, 1400}]; > but I'm getting a meaningless plot so I'd appreciate > if someone could tell me what is what I'm doing wrong. > Ruben Just a thought, but have you checked your function for typos? As you've written it, f is only a function of y .... Dave. ==== > Greetings This problem can be solved by conventional programming, but I wonder if > there is an elegant Mathematica solution ? A list contains pairs of values, with each pair representing the lower > and > upper edge of a sub-range. Some of the sub-ranges partially overlap, > some > fully overlap, others don't overlap at all. The problem is to produce > a > second list that contains the overall upper and lower edges of the > overlapping sub-ranges. A simple example : {{100,200},{150,250},{120,270},{300,400}} would > result > in {{100,270},{300,400}}. In the real case, the input list has several hundred elements and the > output list typically has five elements. I have a working solution based on loops, but there must be a more > elegant > one. I would be very grateful for any suggestions. I'm not sure about elegance, but I will tell you my approach to the problem. The general algorithm would seem to be: sort the ranges such that the first range has a left edge smaller than all other ranges. If two ranges have matching left edges, sort them according to their right edge. While there are two or more ranges in the list, operate on the rest of the list, compare the first range to the result of operating on the rest of the list. If the right edge of the first edge is larger than the left edge of the first element of the result return a list with the first element being a range with the left edge of the first range and the right being the larger of the right edge of the first element or the right edge of the first range in the result of operating on the list. That statement probably isn't very clear, it's a good thing I'm not employed as a teacher. Here is code which is probably easier to follow: (*Handle some degenerate cases here *) compress[{}] := {}; compress[lst : {{_?NumericQ, _?NumericQ}}] := lst (*This pattern is the stopping point of the recursion*) compress[rng : {_?NumericQ, _?NumericQ}, {}] := {rng} (*This function operates on the rest of the list then creates the first element appropriately*) compress[rng : {_? NumericQ, _?NumericQ}, lst : {{_?NumericQ, _?NumericQ} ..}] := With[{tl = compress[First[lst], Rest[lst]]}, If[rng[[2]] > Last[tl][[2]], {rng}, If[rng[[2]] > tl[[1, 1]], {{rng[[1]], If[tl[[1, 2]] > rng[[2]], tl[[1, 2]], rng[[2]]]}, Sequence @@ Rest[tl]}, {rng, Sequence @@ tl}]]] (*This function sorts the list properly then starts the recursion*) compress[lst : {{_?NumericQ, _?NumericQ} ..}] := With[{s = Sort[lst]}, compress[First[s], Rest[s]]] You will probably have to increase $RecursionLimit. This algorithm ran in 2 seconds on a list of 1000 elements on a 1GHz G4 PowerMac. There are probably optimizations that can be made though. Ssezi ==== Ruben, Your plot looks bad because PlotVectorField uses AspectRatio->Automatic, which means that the t and y axes are at the same scale. This makes the plot very high and narrow. Try the following... PlotVectorField[f[t, y], {t, 0, 80}, {y, 0, 1400}, AspectRatio -> 1, Frame -> True, ImageSize -> 450]; David Park djmp@earthlink.net http://home.earthlink.net/~djmp/ but I'm getting a meaningless plot so I'd appreciate if someone could tell me what is what I'm doing wrong. Ruben __________________________________________________ Do you Yahoo!? http://faith.yahoo.com Reply-To: ==== You'll get more efficient methods from others, but I think the following is instructive: data = {{3, 1}, {4, 3}, {3, 2}, {1, 10}, {4, 2}, {1, 6}, {5, 2}, {2, 5}, {7, 1}}; ClearAll[total, count] total[x_] := 0 count[x_] := 0 {total[#[[1]]] += #[[2]], count[#[[1]]]++} & /@ data; ?total ?count {#, total[#]/count[#]} & /@ Union[data[[All, 1]]] DrBob -----Original Message----- of the y values in each group. For the example, i want to report: output = {{1,8},{2,5},{3,1.5},{4,2.5},{5,2},{6,0},{7,1}} In this example, x=6 does not occur so i report the average y[6] = 0. Can anyone suggest a way to do this efficiently?/ many thanks dave +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ David E. Burmaster, Ph.D. Alceon Corporation POBox 382069 (new Box number effective 1 Sep 2001) Harvard Square Station Cambridge, MA 02238-2069 (new ZIP code effective 1 Sep 2001) Voice 617-864-4300 Web http://www.Alceon.com +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Reply-To: ==== It's not clear what you mean by pairs that repeat themselves, but does this do what you want? pairs = Flatten[Outer[List, points, points, 1], 1]; Map[Union, Map[Sort, pairs, 2], 1] DrBob -----Original Message----- themselves, i.e., pair AA. I can use pairs=Outer[List,points,points,1] Then, I have to clear those pairs that repeat themselves, i.e., pair AB and pair BA. Also, when w and h are of the order of 1000s, the computation takes a very long time. Is there a better way of doing the second part of Sincerely Cheng ==== ================================================ Cheng Liu, Ph.D. MST-8, Structure/Property Relations Materials Science and Technology Division Los Alamos National Laboratory Los Alamos, New Mexico 87545 ==== ================================================ ==== David, You will probably get a lot of answers for this. Here is my entry. data = {{3, 1}, {4, 3}, {3, 2}, {1, 10}, {4, 2}, {1, 6}, {5, 2}, {2, 5}, {7, 1}}; First I will show it step-by-step. nmax = 10; Union[Join[data, Table[{i, 0}, {i, 1, nmax}]]] Split[%, #1[[1]] == #2[[1]] & ] Map[Last, %, {2}] (Plus @@ #1/Length[#1] & ) /@ % Transpose[{Range[nmax], %}] giving {{1, 0}, {1, 6}, {1, 10}, {2, 0}, {2, 5}, {3, 0}, {3, 1}, {3, 2}, {4, 0}, {4, 2}, {4, 3}, {5, 0}, {5, 2}, {6, 0}, {7, 0}, {7, 1}, {8, 0}, {9, 0}, {10, 0}} {{{1, 0}, {1, 6}, {1, 10}}, {{2, 0}, {2, 5}}, {{3, 0}, {3, 1}, {3, 2}}, {{4, 0}, {4, 2}, {4, 3}}, {{5, 0}, {5, 2}}, {{6, 0}}, {{7, 0}, {7, 1}}, {{8, 0}}, {{9, 0}}, {{10, 0}}} {{0, 6, 10}, {0, 5}, {0, 1, 2}, {0, 2, 3}, {0, 2}, {0}, {0, 1}, {0}, {0}, {0}} {16/3, 5/2, 1, 5/3, 1, 0, 1/2, 0, 0, 0} {{1, 16/3}, {2, 5/2}, {3, 1}, {4, 5/3}, {5, 1}, {6, 0}, {7, 1/2}, {8, 0}, {9, 0}, {10, 0}} This wraps it into one statement. nmax = 10; Transpose[{Range[nmax], (Plus @@ #1/Length[#1] & ) /@ Map[Last, Split[Union[Join[data, Table[{i, 0}, {i, 1, nmax}]]], #1[[1]] == #2[[1]] & ], {2}]}] {{1, 16/3}, {2, 5/2}, {3, 1}, {4, 5/3}, {5, 1}, {6, 0}, {7, 1/2}, {8, 0}, {9, 0}, {10, 0}} This times a case of 20000 pairs on an 800MHz machine. data2 = Table[{Random[Integer, {1, 100}], Random[Real, {0, 5}]}, {20000}]; nmax = 100; data = data2; Timing[Transpose[{Range[nmax], (Plus @@ #1/Length[#1] & ) /@ Map[Last, Split[Union[Join[data, Table[{i, 0}, {i, 1, nmax}]]], #1[[1]] == #2[[1]] & ], {2}]}]; ] {0.55 Second, Null} David Park djmp@earthlink.net http://home.earthlink.net/~djmp/ output = {{1,8},{2,5},{3,1.5},{4,2.5},{5,2},{6,0},{7,1}} In this example, x=6 does not occur so i report the average y[6] = 0. Can anyone suggest a way to do this efficiently?/ many thanks dave +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ David E. Burmaster, Ph.D. Alceon Corporation POBox 382069 (new Box number effective 1 Sep 2001) Harvard Square Station Cambridge, MA 02238-2069 (new ZIP code effective 1 Sep 2001) Voice 617-864-4300 Web http://www.Alceon.com ==== I'm not 100% sure of this, but how about sortout[intervals_] := ({First[#1], Last[#1]} & ) /@ Sort /@ Flatten /@ Split[Sort[intervals], #1[[2]] >= #2[[1]] & ] ? Example: intrvls = Table[x = Random[Integer, {0, 30}]; {x, x + 2}, {10}] {{18, 20}, {2, 4}, {14, 16}, {1, 3}, {0, 2}, {16, 18}, {7, 9}, {19, 21}, {6, 8}, {6, 8}} sortout[intrvls] {{0, 4}, {6, 9}, {14, 21}} --- Selwyn Hollis > Greetings This problem can be solved by conventional programming, but I wonder if > there is an elegant Mathematica solution ? A list contains pairs of values, with each pair representing the lower and > upper edge of a sub-range. Some of the sub-ranges partially overlap, some > fully overlap, others don't overlap at all. The problem is to produce a > second list that contains the overall upper and lower edges of the > overlapping sub-ranges. A simple example : {{100,200},{150,250},{120,270},{300,400}} would result > in {{100,270},{300,400}}. In the real case, the input list has several hundred elements and the > output list typically has five elements. I have a working solution based on loops, but there must be a more elegant > one. I would be very grateful for any suggestions. John Leary > ==== John, Use the Interval routine. Interval @@ {{100, 200}, {150, 250}, {120, 270}, {300, 400}} List @@ % giving... Interval[{100, 270}, {300, 400}] {{100, 270}, {300, 400}} David Park djmp@earthlink.net http://home.earthlink.net/~djmp/ in {{100,270},{300,400}}. In the real case, the input list has several hundred elements and the output list typically has five elements. I have a working solution based on loops, but there must be a more elegant one. I would be very grateful for any suggestions. John Leary ==== >I have a list of 17,000+ {x,y} pairs of data each x value is a positive integer from 1 to 100+ each y value is a positive real number As a *short* example, let's consider: data = {{3,1},{4,3},{3,2},{1,10},{4,2},{1,6},{5,2},{2,5},{7,1}} I want to group the data by the x value and report the arithmetic average >of the y values in each group. For the example, i want to report: output = {{1,8},{2,5},{3,1.5},{4,2.5},{5,2},{6,0},{7,1}} In this example, x=6 does not occur so i report the average y[6] = 0. Can anyone suggest a way to do this efficiently?/ > The basic approach is (Plus @@ #)/Length[#] & /@ Split[Sort[data], #1[[1]] == #2[[1]] &] {{1, 8}, {2, 5}, {3, 3/2}, {4, 5/2}, {5, 2}, {7, 1}} To fill in the gaps: dataAvg[data_] := Module[{val = First /@ data, xData}, xData = Join[data, {#, 0} & /@ Complement[Range[Max[val]], val]]; (Plus @@ #)/Length[#] & /@ Split[Sort[xData], #1[[1]] == #2[[1]] &]]; dataAvg[data] {{1, 8}, {2, 5}, {3, 3/2}, {4, 5/2}, {5, 2}, {6, 0}, {7, 1}} Bob Hanlon ==== >IÇm new, so IÇm sorry if the question is so easy, I really think that >is easy, but I donÇt know how to do it. I have this: b=n-m a=x-b Y=3a + 4a^2 and the program show me this: >3(x-n+m) + 4(x-n+m)^2 or, something like that, the problem is that I want the program show >me Y in function of b, or sometimes in function of a, something like >this: Y=3(x-b) + 4(x-b)^2 or Y=3a + 4a^2 > Solve[{b == n - m, a == x - b, Y == 3a + 4a^2}, Y] {{Y -> a*(4*a + 3)}} Solve[{b == n - m, a == x - b, Y == 3a + 4a^2}, Y, a] // FullSimplify {{Y -> (4*b - 4*x - 3)*(b - x)}} Bob Hanlon ==== Boole (defined in the AddOn package Calculus`Integration`) is not written to deal with symbolic parameters. So you may do better to use better to use the built-in UnitStep function for which Mathematica knows more rules. In fact, if you do not mind getting an error message you can get something like your answer by mixing UnitStep and Boole: << Calculus`Integration` In[2]:= FullSimplify[Integrate[Boole[0 < x < y < 1]*UnitStep[z - y], {z, 0, 1}],y>0] Integrate::region: The region defined by 0ç.8czç.8c1&&00] Out[3]= -(-1+y) UnitStep[x,1-y,-x+y] Note that this expresses exactly the same condition as (1-y)Boole[0 Suppose f(x,y,z)=Boole[0 over, say, z, i.e. Integrate[f[x,y,z],{z,0,1}]. One would expect to see > the > output like this (1-y)Boole[0 as an argument in the Boole function). Instead, an error appears > (warning) that the integration cannot be performed. How to resolve this issue so it produces a desired answer? Janusz. > ==== >This problem can be solved by conventional programming, but I wonder if there is an elegant Mathematica solution ? A list contains pairs of values, with each pair representing the lower >and >upper edge of a sub-range. Some of the sub-ranges partially overlap, some fully overlap, others don't overlap at all. The problem is to produce >a >second list that contains the overall upper and lower edges of the >overlapping sub-ranges. A simple example : {{100,200},{150,250},{120,270},{300,400}} would result in {{100,270},{300,400}}. In the real case, the input list has several hundred elements and the >output list typically has five elements. I have a working solution based on loops, but there must be a more elegant one. I would be very grateful for any suggestions. > lst = {{100, 200}, {150, 250}, {120, 270}, {300, 400}}; List @@ Union @@ Interval /@ lst {{100, 270}, {300, 400}} List @@ Interval[Sequence @@ lst] {{100, 270}, {300, 400}} Bob Hanlon ==== In my opinion, the best way to learn C++ and code algorithms is to simply start using it. There is oodles of source on the web and at www.sf.net to help you along too. Matt -- http://mffm.darktech.org WSOLA TimeScale Audio Mod : http://mffmtimescale.sourceforge.net/ FFTw C++ : http://mffmfftwrapper.sourceforge.net/ Vector Bass : http://mffmvectorbass.sourceforge.net/ Multimedia Time Code : http://mffmtimecode.sourceforge.net/ ==== do any of you have a working qsum author indentification workbook you are willing to share? thanks ==== >-----Original Message----- >Sent: Wednesday, October 16, 2002 8:26 PM >To: mathgroup@smc.vnet.net >Dear group, I have the following question regarding a lengthy calculation >using Mathematica: For given w points in x direction and h points in y direction, I can >construct all the points using h=10; w=8; > points=Flatten[Transpose[Outer[List,Range[w],Range[h]]],1] Next, I need to find all the possible pairs of point including points >themselves, i.e., pair AA. I can use pairs=Outer[List,points,points,1] Then, I have to clear those pairs that repeat themselves, >i.e., pair AB and >pair BA. Also, when w and h are of the order of 1000s, the >computation >takes a very long time. Is there a better way of doing the >second part of Sincerely Cheng > ==== ================================================ >Cheng Liu, Ph.D. >MST-8, Structure/Property Relations >Materials Science and Technology Division >Los Alamos National Laboratory >Los Alamos, New Mexico 87545 ==== ================================================ Cheng, you didn't tell about a specific order of your resulting pairs, so you not need Transpose in your first line: points = Flatten[Outer[List, Range[w], Range[h]], 1] Now build the pairs (for sake of clarity , let me call the points points = {p1, p2, p3, p4, p5}, of course you don't do that): Flatten[ With[{l = Length[points]}, Array[If[#1>#2, Unevaluated[Sequence[]], points[[{#1, #2}]]]&, {l, l}]], 1] {{p1, p1}, {p1, p2}, {p1, p3}, {p1, p4}, {p1, p5}, {p2, p2}, {p2, p3}, {p2, p4}, {p2, p5}, {p3, p3}, {p3, p4}, {p3, p5}, {p4, p4}, {p4, p5}, {p5, p5}} Alternatively you might do pairs = Outer[List, points, points, 1] Flatten[MapIndexed[Drop[#1, First[#2] - 1] &, pairs], 1] {{p1, p1}, {p1, p2}, {p1, p3}, {p1, p4}, {p1, p5}, {p2, p2}, {p2, p3}, {p2, p4}, {p2, p5}, {p3, p3}, {p3, p4}, {p3, p5}, {p4, p4}, {p4, p5}, {p5, p5}} ...try out which is faster. -- Hartmut Wolf ==== When I first sent my answer I thought there was no interaction between Boole and UnitStep at all,and that one could safely use UnitStep after loading the Calculus`Integration` package, but there seems to be more to it than I had assumed and it is not necessarily for the best. Consider first the following: In[1]:= FullSimplify[Integrate[UnitStep[x,y-x,z-y,1-z],{z,0,1}],y>0] Out[1]= -(-1+y) UnitStep[x,1-y,-x+y] No surprises here. Now let's load the package: In[2]:= <0] Integrate::region: The region defined by 0.89.81óz.89.81ó1&&x.89 .81«0&&-x+y.89.81«0&&1-z.89.8126 40&&-y+z.89.81«0 could not be broken down into cylinders. Integrate::region: The region defined by y.89.81óz.89.81ó1&&x.89 .81«0&&-x+y.89.81«0&&1-z.89.8126 40 could not be broken down into cylinders. Integrate::region: The region defined by -1.89.81óz.89.81ó-y&&x21 1.81«0&&-x+y.89.81«0 could not be broken down into cylinders. General::stop: Further output of Integrate::region will be suppressed during this calculation. Out[3]= -(-1+y) UnitStep[x,1-y,-x+y] Clearly an attempt was made to decompose this into cylinders with respect to z (using CAD) which of course failed. Fortunately we still get the right answer. Secondly, the package actually contains some interesting functions which have been commented out and were apparently intended for future development. On of them is: removeUnitStep[expr_] := ReplaceRepeated[expr, UnitStep[e__] :> Boole[Apply[And, Map[(# .89.81« 0) &, {e}]]]]; Using it we get: removeUnitStep[FullSimplify[Integrate[UnitStep[x,y-x,z-y,1- z],{z,0,1}],y>0]] (error messages removed) -(-1+y) Boole[x.89.81«0&&x.89.81óy&&y.89[ CapitalARing]ó1] Mathematica can't convert this into (1-y)Boole[0<=x<=y<=1] but it can do the converse: Map[LogicalExpand,(1-y)Boole[0.89.81óx.89[Capita lARing]óy.89.81ó1],{2}] (1-y) Boole[0.89.81óx&&x.89.81ó y&&y.89.81ó1] Andrzej Kozlowski Yokohama, Japan http://www.mimuw.edu.pl/~akoz/ http://platon.c.u-tokyo.ac.jp/andrzej/ > Boole (defined in the AddOn package Calculus`Integration`) is not > written to deal with symbolic parameters. So you may do better to use > better to use the built-in UnitStep function for which Mathematica > knows more rules. In fact, if you do not mind getting an error message > you can get something like your answer by mixing UnitStep and Boole: << Calculus`Integration` In[2]:= > FullSimplify[Integrate[Boole[0 < x < y < 1]*UnitStep[z - y], > {z, 0, 1}],y>0] Integrate::region: The region defined by 0.8cÁ.8cåz.8cÁ .8cå1&&0 could > not be > broken down into cylinders. Out[2]= > -(-1+y) Boole[0 already included in Boole, but the Integration package does not have > rules for combining Bool and UnitStep (don't forget that Boole is not a > built in function!). If you dispense with Bool altogether you get: In[3]:= > FullSimplify[Integrate[UnitStep[x,y-x,z-y,1-z],{z,0,1}],y>0] Out[3]= > -(-1+y) UnitStep[x,1-y,-x+y] Note that this expresses exactly the same condition as > (1-y)Boole[0 Yokohama, Japan > http://www.mimuw.edu.pl/~akoz/ > http://platon.c.u-tokyo.ac.jp/andrzej/ > and also try to use proper Mathematica >> Suppose f(x,y,z)=Boole[0> over, say, z, i.e. Integrate[f[x,y,z],{z,0,1}]. One would expect to >> see >> the >> output like this (1-y)Boole[0> as an argument in the Boole function). Instead, an error appears >> (warning) that the integration cannot be performed. >> How to resolve this issue so it produces a desired answer? >> Janusz. > > Reply-To: ==== --- not as a poor a showing for my simple-minded method as I feared. data2 = Table[{Random[Integer, {1, 100}], Random[Real, {0, 5}]}, {20000}]; nmax = 100; data = data2; Timing[dave = Transpose[{Range[nmax], (Plus @@ #1/Length[#1] &) /@ Map[Last, Split[Union[Join[data, Table[{i, 0}, {i, 1, nmax}]]], #1[[1]] == #2[[1]] &], {2}]}];] ClearAll[total, count] total[x_] := 0 count[x_] := 0 Timing[{total@#[[1]] += #[[2]], count[#[[1]]]++} & /@ data; brt = {#, total[#]/count[#]} & /@ Union[data[[All, 1]]]; ] {0.20299999999999985*Second, Null} {0.5619999999999998*Second, Null} DrBob -----Original Message----- Split[%, #1[[1]] == #2[[1]] & ] Map[Last, %, {2}] (Plus @@ #1/Length[#1] & ) /@ % Transpose[{Range[nmax], %}] giving {{1, 0}, {1, 6}, {1, 10}, {2, 0}, {2, 5}, {3, 0}, {3, 1}, {3, 2}, {4, 0}, {4, 2}, {4, 3}, {5, 0}, {5, 2}, {6, 0}, {7, 0}, {7, 1}, {8, 0}, {9, 0}, {10, 0}} {{{1, 0}, {1, 6}, {1, 10}}, {{2, 0}, {2, 5}}, {{3, 0}, {3, 1}, {3, 2}}, {{4, 0}, {4, 2}, {4, 3}}, {{5, 0}, {5, 2}}, {{6, 0}}, {{7, 0}, {7, 1}}, {{8, 0}}, {{9, 0}}, {{10, 0}}} {{0, 6, 10}, {0, 5}, {0, 1, 2}, {0, 2, 3}, {0, 2}, {0}, {0, 1}, {0}, {0}, {0}} {16/3, 5/2, 1, 5/3, 1, 0, 1/2, 0, 0, 0} {{1, 16/3}, {2, 5/2}, {3, 1}, {4, 5/3}, {5, 1}, {6, 0}, {7, 1/2}, {8, 0}, {9, 0}, {10, 0}} This wraps it into one statement. nmax = 10; Transpose[{Range[nmax], (Plus @@ #1/Length[#1] & ) /@ Map[Last, Split[Union[Join[data, Table[{i, 0}, {i, 1, nmax}]]], #1[[1]] == #2[[1]] & ], {2}]}] {{1, 16/3}, {2, 5/2}, {3, 1}, {4, 5/3}, {5, 1}, {6, 0}, {7, 1/2}, {8, 0}, {9, 0}, {10, 0}} This times a case of 20000 pairs on an 800MHz machine. data2 = Table[{Random[Integer, {1, 100}], Random[Real, {0, 5}]}, {20000}]; nmax = 100; data = data2; Timing[Transpose[{Range[nmax], (Plus @@ #1/Length[#1] & ) /@ Map[Last, Split[Union[Join[data, Table[{i, 0}, {i, 1, nmax}]]], #1[[1]] == #2[[1]] & ], {2}]}]; ] {0.55 Second, Null} David Park djmp@earthlink.net http://home.earthlink.net/~djmp/ For the example, i want to report: output = {{1,8},{2,5},{3,1.5},{4,2.5},{5,2},{6,0},{7,1}} In this example, x=6 does not occur so i report the average y[6] = 0. Can anyone suggest a way to do this efficiently?/ many thanks dave +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ David E. Burmaster, Ph.D. Alceon Corporation POBox 382069 (new Box number effective 1 Sep 2001) Harvard Square Station Cambridge, MA 02238-2069 (new ZIP code effective 1 Sep 2001) Voice 617-864-4300 Web http://www.Alceon.com +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ==== >-----Original Message----- >Sent: Wednesday, October 16, 2002 8:26 PM >To: mathgroup@smc.vnet.net >Greetings This problem can be solved by conventional programming, but I >wonder if >there is an elegant Mathematica solution ? A list contains pairs of values, with each pair representing >the lower and >upper edge of a sub-range. Some of the sub-ranges partially >overlap, some >fully overlap, others don't overlap at all. The problem is to >produce a >second list that contains the overall upper and lower edges of the >overlapping sub-ranges. A simple example : {{100,200},{150,250},{120,270},{300,400}} >would result >in {{100,270},{300,400}}. In the real case, the input list has several hundred elements and the >output list typically has five elements. I have a working solution based on loops, but there must be a >more elegant >one. I would be very grateful for any suggestions. John Leary John, several proposals (without any attempt to moduralize): (1) use IntervalUnion: List @@ IntervalUnion @@ Interval /@ list (2) use Split (it's a little bit tricky to be correct): high = Sequence[]; {#[[1, 1]], Max[#[[All, -1]]]} & /@ Split[Sort[list], (high = Max[high,Last[#1]]) >= First[#2] || (high = Last[#1])&] (3) procedural programming: maxExtends[list_] := (sl = Sort[list]; length = Length[sl]; r = collect[]; i = 1; While[i <= length, {low, high} = sl[[i]]; If[++i <= length, {curlow, curhigh} = sl[[i]]; While[high >= curlow && (high = Max[high, curhigh]; ++i <= length), {curlow, curhigh} = sl[[i]] ]]; r = collect[r, {low, high}] ]; List @@ Flatten[r]) Let's do some benchmarks: 10,000 Intervals: list = {# - Random[], # + Random[]} & /@ NestList[# + Random[] &, 0, 10000]; List @@ IntervalUnion @@ Interval /@ list // Length // Timing {2.503 Second, 1181} high = Sequence[]; {#[[1, 1]], Max[#[[All, -1]]]} & /@ Split[Sort[ list], (high = Max[high, Last[#1]]) >= First[#2] || (high = Last[#1]) &] // Length // Timing {2.934 Second, 1181} maxExtends[list] // Length // Timing {3.926 Second, 1181} The corresponding results for 100,000 Intervals: {27.329 Second, 11266} {30.234 Second, 11266} {35.791 Second, 11266} and for 500,000 Intervals {144.058 Second, 56728} {154.782 Second, 56728} {181.111 Second, 56728} To look at scaling behaviour I just collected the prior results IntervalUnion: {%355, %345 , %350}[[All, 1, 1]] {2.503, 27.329, 144.058} % // {#[[2]]/(10*#[[1]]), #[[3]]/(5*#[[2]])} & {1.09185, 1.05425} Split: {%357, %347, %352}[[All, 1, 1]] {2.934, 30.234, 154.782} % // {#[[2]]/(10*#[[1]]), #[[3]]/(5*#[[2]])} & {1.03047, 1.02389} Procedural: {%358, %348, %353}[[All, 1, 1]] {3.926, 35.791, 181.111} % // {#[[2]]/(10*#[[1]]), #[[3]]/(5*#[[2]])} & {0.91164, 1.01205} Due to Sort, the Split and the Procedural versions should behave as O[n log n], I'm not shure whether IntervalUnion does (seems to be a little bit more progressive at costs). -- Hartmut Wolf ==== I use NDSolve to approximate a 6-dimensional and highly non-linear dynamic system numerically. QUESTION: What are usual (and simple) techniques to check the reliability of the resulting numerical solution? So far I tried to use another solution algorithm by switching the Option Method from Method->Automatic to Method->RungeKutta. Apart from numerical differences, the solution qualitatively remained unaffected (I am mainly interested in the qualitative characteristics of the solution). I would be grateful for any hint! ==== Here is a version which compresses the list going forward, it does not seem to be significantly faster than the other version but the memory use should be slightly smaller: (*Handle some degenerate cases here *) fcompress[{}]:={}; (*This pattern is the stopping point of the recursion*) fcompress[lst:{{_?NumericQ,_?NumericQ}}]:=lst (*This function needs to be called with a sorted list*) fcompress[lst:{{_?NumericQ,_?NumericQ}..}]:=With[{rng1=lst[[1]],rng2=lst [[2]]},If[rng1[[2]] lst = {{100, 200}, {150, 250}, {120, 270}, {300, 400}}; List @@ Union @@ Interval /@ lst {{100, 270}, {300, 400}} List @@ Interval[Sequence @@ lst] {{100, 270}, {300, 400}} > Forgot about the Interval function, but to continue the pattern List@@Interval@@lst also works. I personally prefer not to use Sequence objects if I can avoid them. Sseziwa PS - It is much faster than my previous posts. ==== Dear group, very clear. Suppose that I have a list of points {p1,p2, ..., pn}, I try to find all possible pairs of them. The pairs may include {pi,pi}, but {pi,pj} and {pj,pi} are considered the same and only one is counted. That said. After some try and error, I came to the following way: h=4;w=5; points=Flatten[Outer[List,Range[w],Range[h]],1]; pairs=Flatten[Map[Outer[List,{#},Drop[points,Position[points,#][[1,1]]-1],1] [ [1]]&,points],1]; The speed of the above calculation is reasonably fast. But I run into the memory problem. For example, for h=64 and w=64, the length of the list pairs will be w*h (w*h+1)/2 = 8,390,656. In my case, The numbers for h and w will be a lot larger than 64. How can I get around this memory problem or that's the dead end for my calculation (I do have 1 GB physical mem in Cheng >Dear group, I have the following question regarding a lengthy calculation >using Mathematica: For given w points in x direction and h points in y direction, I can >construct all the points using h=10; w=8; > points=Flatten[Transpose[Outer[List,Range[w],Range[h]]],1] Next, I need to find all the possible pairs of point including points >themselves, i.e., pair AA. I can use pairs=Outer[List,points,points,1] Then, I have to clear those pairs that repeat themselves, i.e., pair AB and >pair BA. Also, when w and h are of the order of 1000s, the computation >takes a very long time. Is there a better way of doing the second part of Sincerely Cheng > ==== ================================================ >Cheng Liu, Ph.D. >MST-8, Structure/Property Relations >Materials Science and Technology Division >Los Alamos National Laboratory >Los Alamos, New Mexico 87545 ==== ================================================ ==== ================================================ Cheng Liu, Ph.D. MST-8, Structure/Property Relations Materials Science and Technology Division Los Alamos National Laboratory Los Alamos, New Mexico 87545 ==== ================================================ ==== Just two small (maybe relevant) notes: ----- Original Message ----- even some more than just today's machine precision. > Sorry, but that's not as profound as it sounds. The speed of light is > indeed a very specific number, but that doesn't mean we can measure it > precisely. Instead, like E or Pi or 2 or Sqrt[7], it's a defined > constant and -- unlike E or Pi or Sqrt[7] -- the definition doesn't > allow us to compute it with arbitrary precision. Yes, it's defined now > so that it can pretend to unlimited precision -- but that only means > meters (or seconds, take your pick) aren't defined precisely. The speed of light is *postulated* (rather then defined) constant, and (since 1983.) in the international metric (aka SI) system of units one meter is *defined* to be the distance traveled by light during (exactely) 1 / 299792458 seconds (second is defined some other way). So it follows that the (exact) value for the speed of light is 299792458 m / s. However, the problem is not in the correctness of this particular value (it is actually just a convinient convention), but in the postulate that this value is the same for all observers. To be specific, if there are two observers moving relative to eachother and measuring the speed of light using two (in principle) _identical_ experimental devices, the question is whether they get the same result, ie. c' = c. Experimentally speaking, one wihses to know how much (if at all) the quantity (c'-c)/c differs from zero. If someone measures a nonzero value, that _would_ actually be for a Nobel... Marko > -----Original Message----- > To: mathgroup@smc.vnet.net > In the real world of physics there are several subatomic level > processes > which can only be distinguished by small changes in the n-th decimal > place. > But there is one example which is fairly easy to comprehend, and that is > the > constancy of the speed of light in a vacuum regardless of reference > frame, > as proposed in Einstein's special theory of relativity. If this were > true > only to the 9th or 10th decimal place, or, for that matter, to the > 50th > place, then whoever managed to show that it was not really a constant > would > certainly be in Nobel Prize territory, and much of modern physics would > need > a rewrite. Kevin Greetings, I have read with great interest this lively debate on numerical > prcesion > and > accuracy. As I work in the fields of finance and economics, where we > feel > ourselves blessed if we get three digits of accuracy, I'm curious as > to > what > scientific endeavors require 50+ digits of precision? As I recall > there > are > some areas, such as high energy physics and some elements of > astronomy, > that > might require so many digits in some circumstances. Are there > others? > -Mark > Reply-To: ==== Sorry, but that's not as profound as it sounds. The speed of light is indeed a very specific number, but that doesn't mean we can measure it precisely. Instead, like E or Pi or 2 or Sqrt[7], it's a defined constant and -- unlike E or Pi or Sqrt[7] -- the definition doesn't allow us to compute it with arbitrary precision. Yes, it's defined now so that it can pretend to unlimited precision -- but that only means meters (or seconds, take your pick) aren't defined precisely. For anything we can measure (or even COUNT, in the real world), I suspect 16-digit machine precision is more than enough. Bobby -----Original Message----- as proposed in Einstein's special theory of relativity. If this were true only to the 9th or 10th decimal place, or, for that matter, to the 50th place, then whoever managed to show that it was not really a constant would certainly be in Nobel Prize territory, and much of modern physics would need a rewrite. Kevin > Greetings, I have read with great interest this lively debate on numerical prcesion > and > accuracy. As I work in the fields of finance and economics, where we feel > ourselves blessed if we get three digits of accuracy, I'm curious as to > what > scientific endeavors require 50+ digits of precision? As I recall there > are > some areas, such as high energy physics and some elements of astronomy, > that > might require so many digits in some circumstances. Are there others? > -Mark ==== I have 2 or more separate Plots which have different y but the same x axes. Like: Plot[Sin[x],{x,0,10}]; Plot[1000Sin[x],{x,0,10}]; At the display and printout, the y axes are not aligned. Even using the PlotRegion and ImageSize Options doesn't help. The only way I found was to align it manually with the mouse. Is there a package that solves the problem? The problem is, I have many, many plots to align... I use Mathematica 4.0.2.0 on a Windows 2000 PC. Who can help? Max ==== I have 2 or more separate Plots which have different y but the same x axes. Like: Plot[Sin[x],{x,0,10}]; Plot[1000Sin[x],{x,0,10}]; At the display and printout, the y axes are not aligned. Even using the PlotRegion and ImageSize Options doesn't help. The only way I found was to align it manually with the mouse. Is there a package that solves the problem? The problem is, I have many, many plots to align... I use Mathematica 4.0.2.0 on a Windows 2000 PC. Who can help? Max Answers please to: ulbrich@biochem.mpg.de ==== > A simple example : {{100,200},{150,250},{120,270},{300,400}} would result > in {{100,270},{300,400}}. In the real case, the input list has several hundred elements and the > output list typically has five elements. I have a working solution based on loops, but there must be a more elegant > one. I would be very grateful for any suggestions. Until recently, this could have been tedious, but now, tada!: In[31] := data = {{100, 200}, {150, 250}, {120, 270}, {300, 400}}; List @@ IntervalUnion @@ Interval /@ data Out[31] = {{100, 270}, {300, 400}} Tom Burton ==== >-----Original Message----- >Sent: Wednesday, October 16, 2002 8:26 PM >To: mathgroup@smc.vnet.net >Dear Fellows in MathGroup, I have a list of 17,000+ {x,y} pairs of data each x value is a positive integer from 1 to 100+ each y value is a positive real number As a *short* example, let's consider: data = {{3,1},{4,3},{3,2},{1,10},{4,2},{1,6},{5,2},{2,5},{7,1}} I want to group the data by the x value and report the >arithmetic average >of the y values in each group. For the example, i want to report: output = {{1,8},{2,5},{3,1.5},{4,2.5},{5,2},{6,0},{7,1}} In this example, x=6 does not occur so i report the average y[6] = 0. Can anyone suggest a way to do this efficiently?/ many thanks >dave +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ >David E. Burmaster, Ph.D. >Alceon Corporation >POBox 382069 (new Box number effective 1 Sep 2001) >Harvard Square Station >Cambridge, MA 02238-2069 (new ZIP code effective 1 Sep 2001) Voice 617-864-4300 Web http://www.Alceon.com >+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Dave, my first attempt used the same reasoning as Daniel Lichtblau proposed (and came out only slightly faster than his). However the Sort/Split idea as brought forward by Bob Hanlon and Allan Hayes is much faster; Bobby Treat's version turns out to be slower (than Daniel's and mine). To reach comparable results, I slightly modified Allan's solution (which was fastest): data = Table[{Random[Integer, {1, 98}], Random[]}, {20000}]; (f[x_] = 0; ((f[#1[[1, 1]]] = Plus @@ #1[[All, 2]]/Length[#1]) &) /@ Split[Sort[data], #1[[1]] == #2[[1]] &]; r4 = {#, f[#]} & /@ Range[98];) // Timing {3.045 Second, Null} So I reconsidered that idea and found a solution which is nearly twice as fast: binnedAverage2[data_, max_] := Module[{v, i, ix, ixx, ixxx}, {i, v} = With[{rr = Range[max]}, Transpose[Sort[Join[data, Transpose[{rr, rr - rr}]]]]]; ix = Split[i]; ixx = FoldList[Plus[#1, Length[#2]] &, 0, ix]; ixxx = Transpose[Transpose[Partition[ixx, 2, 1]] + {1, 0}]; Transpose[{First /@ ix, Plus @@ #/Max[Length[#] - 1, 1] &[Take[v, #]] & /@ ixxx}]] (r7 = binnedAverage2[data, 98]); // Timing {1.612 Second, Null} r7 == r4 True -- Hartmut Wolf ==== The notebook code below contains ways of speeding up some uses of Table, Sort and Split and avoiding building up inefficienlty large intermediate expressions. I put them together while looking into the various solutions in the thread grouping and averaging {x,y} pairs of data -- they helped give what seems to be the fastest uncompiled solution yet. Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 To make a notebook from the following, copy between the two lines +++++++, past into a notebook, click Yes in the panel that pops up. +++++++ Notebook[{ Cell[CellGroupData[{ Cell[TextData[{ Speedups n, StyleBox[Table, Sort, Split, FontSize->18] }], Title], Cell[Allan Hayes October 2002, Subsubtitle], Cell[BoxData[ (Off[General::, General::])], Input, InitializationCell->True], Cell[CellGroupData[{ Cell[General Techniques, Section], Cell[CellGroupData[{ Cell[Table, Subsection], Cell[TextData[{ The following illustrates a generally applicable point: that inserting large expressions directly using a function or , StyleBox[With, FontFamily->Courier], can lead to inefficiently large expressions. This can be avoided by storing in a local variable. }], Text], Cell[BoxData[{ ((ranges = Table[{1, 200}, {10000}];)), n, ((data = Range[100000];))}], Input], Cell[This is quite slow., Text], Cell[BoxData[ (takesX[data_, ranges_] := [IndentingNewLine](Take[data, #] &) /@ ranges)], Input], Cell[CellGroupData[{ Cell[BoxData[ ((takesX[data, ranges];) // Timing)], Input], Cell[BoxData[ ({3.789999999999992` Second, Null})], Output] }, Closed]], Cell[TextData[{ The use of the local variable , StyleBox[dd, FontFamily->Courier], below gives much quicker code }], Text], Cell[BoxData[ (takes[data_, ranges_] := [IndentingNewLine]Module[{dd = data}, (Take[dd, #] &) /@ ranges])], Input, InitializationCell->True, FontColor->RGBColor[0, 0, 1]], Cell[CellGroupData[{ Cell[BoxData[ ((takes[data, ranges];) // Timing)], Input], Cell[BoxData[ ({1.4899999999999807` Second, Null})], Output] }, Closed]], Cell[TextData[{ We might have expected that the direct insertion of , StyleBox[data, FontFamily->Courier], into , StyleBox[Take[data,#]&, FontFamily->Courier New], would be quicker, but it looks as if the very big expression that is built up when using this technique slows things down.nIn the second form , StyleBox[dd, FontFamily->Courier], is evaluated once as each expressionntt, StyleBox[Take[dd,#]&[{1,200}], FontFamily->Courier New], nis evaluated. }], Text], Cell[TextData[{ As an extra check, use , StyleBox[With, FontFamily->Courier], to input directly: }], Text], Cell[BoxData[ (takesY[data_, ranges_] := [IndentingNewLine]With[{dd = data}, (Take[dd, #] &) /@ ranges])], Input], Cell[CellGroupData[{ Cell[BoxData[ ((takesY[data, ranges];) // Timing)], Input], Cell[BoxData[ ({3.240000000000002` Second, Null})], Output] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[Sort, Subsection], Cell[BoxData[{ ((k = 2;)), [IndentingNewLine], ((data = ((SeedRandom[0]; Table[Random[ Integer, {1, 100}], {100000}, {k}]));))}], Input], Cell[Let's sort the data:, Text], Cell[CellGroupData[{ Cell[BoxData[ ((Sort[data];) // Timing)], Input], Cell[BoxData[ ({1.9299999999999997` Second, Null})], Output] }, Closed]], Cell[< Sorting by the first coordinate might be expected to be quicker, but the following slows the calculation down markedly! >, Text], Cell[CellGroupData[{ Cell[BoxData[ ((Sort[data, #1[([1])] [LessEqual] #2[([1])] &];) // Timing)], Input], Cell[BoxData[ ({150.16` Second, Null})], Output] }, Closed]], Cell[TextData[{ But using , StyleBox[Ordering, FontFamily->Courier], , as below will usually speed things up (the advantage increases with the size of k). }], Text], Cell[BoxData[ (sortn[dat_, n_] := [IndentingNewLine]dat[([Ordering[ dat[([All, n])]]])])], Input, InitializationCell->True, FontColor->RGBColor[0, 0, 1]], Cell[CellGroupData[{ Cell[BoxData[ ((sortn[data, 1];) // Timing)], Input], Cell[BoxData[ ({1.4799999999999898` Second, Null})], Output] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[Split, Subsection], Cell[BoxData[{ ((data = ((SeedRandom[0]; Table[Random[ Integer, {1, 100}], {100000}, {2}]));)), n, ((sd = Sort[data];))}], Input], Cell[TextData[{ Split , StyleBox[data, FontFamily->Courier], by the first coordinate }], Text], Cell[CellGroupData[{ Cell[BoxData[ ((Split[sd, #1[([1])] [Equal] #2[([1])] &];) // Timing)], Input], Cell[BoxData[ ({10.549999999999999` Second, Null})], Output] }, Closed]], Cell[The following (following Hartmut Wolf) is quicker , Text], Cell[BoxData[ (splitnX[dat_, n_] := [IndentingNewLine]Module[{s, p, tks}, [IndentingNewLine]s = Split[dat[([All, n])]]; p = FoldList[ #1 + #2 &, 0, Length /@ s]; [IndentingNewLine](Take[ dat, #] &) /@ Transpose[{Drop[p, (-1)] + 1, Rest[p]}][IndentingNewLine]])], Input], Cell[CellGroupData[{ Cell[BoxData[ ((splitnX[sd, 1];) // Timing)], Input], Cell[BoxData[ ({4.939999999999998` Second, Null})], Output] }, Closed]], Cell[TextData[{ And using the trick we found above when looking at, StyleBox[ Table , FontFamily->Courier New], it gets still quicker }], Text], Cell[BoxData[ (splitn[dat_, n_] := [IndentingNewLine]Module[{s, p, tks, dd}, [IndentingNewLine]ss = (s = Split[dat[([All, n])]]); p = FoldList[ #1 + #2 &, 0, Length /@ s]; [IndentingNewLine]dd = dat; [IndentingNewLine](Take[dd, #] &) /@ Transpose[{Drop[p, (-1)] + 1, Rest[p]}][IndentingNewLine]])], Input, InitializationCell->True, FontColor->RGBColor[0, 0, 1]], Cell[CellGroupData[{ Cell[BoxData[ ((splitn[sd, 1];) // Timing)], Input], Cell[BoxData[ ({2.09` Second, Null})], Output] }, Closed]] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[Application to Binned Averages, Section], Cell[< Daniel Lichtblau's elegant compiled solution is quickest. The uncompiled version, 2., comes fairly close, but will probably fall further behind with very big inputs. >, Text], Cell[TextData[{ To make the timings more comparable, each section below starts with , StyleBox[Quit, FontFamily->Courier New], . This must be evaluated separately before the rest of the code. You will be asked if you want to evaluate the initialization cells when continuing with rest of the code - click Yes. }], Text], Cell[CellGroupData[{ Cell[1. Hartmut Wolf, Subsection, CellDingbat->None, CellMargins->{{10, Inherited}, {Inherited, Inherited}}], Cell[BoxData[ (Quit)], Input], Cell[BoxData[ ((data = ((SeedRandom[0]; Table[Random[ Integer, {1, 100}], {100000}, {2}]));))], Input], Cell[BoxData[ (((binnedAverage1[data_, max_] := Module[{v, i, ix, ixx, ixxx}, {i, v} = With[{rr = Range[max]}, Transpose[ Sort[Join[data, Transpose[{rr, rr - rr}]]]]]; [IndentingNewLine]ix = Split[i]; [IndentingNewLine]ixx = FoldList[Plus[#1, Length[#2]] &, 0, ix]; [IndentingNewLine]ixxx = Transpose[ Transpose[Partition[ixx, 2, 1]] + {1, 0}]; [IndentingNewLine]Transpose[{First /@ ix, ((((Plus @@ #))/Max[Length[#] - 1, 1] &)[ Take[v, #]] &) /@ ixxx}]])(n) ))], Input], Cell[CellGroupData[{ Cell[BoxData[ (Table[((((r1 = binnedAverage1[data, 100]));) // Timing) // First, {5}])], Input], Cell[BoxData[ ({3.790000000000001` Second, 4.609999999999999` Second, 4.17` Second, 5.` Second, 4.010000000000002` Second})], Output] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[2. Allan Hayes, Subsection, CellDingbat->None, CellMargins->{{10, Inherited}, {Inherited, Inherited}}], Cell[TextData[{ Note that , StyleBox[Tr[lst], FontFamily->Courier], StyleBox[ is used , FontFamily->Times New Roman], StyleBox[ins, FontFamily->Times New Roman], tead o, StyleBox[f , FontFamily->Times New Roman], StyleBox[Plus@@lst, FontFamily->Courier], StyleBox[ , FontFamily->Times New Roman], StyleBox[to sum, FontFamily->Times New Roman], StyleBox[ , FontFamily->Times New Roman], StyleBox[lst, FontFamily->Courier], . }], Text], Cell[BoxData[ (Quit)], Input], Cell[BoxData[ ((data = ((SeedRandom[0]; Table[Random[ Integer, {1, 100}], {100000}, {2}]));))], Input], Cell[BoxData[ (binnedAverage2[data_, max_: 0] := [IndentingNewLine]Module[{spl, avs, inds, dd, zs}, [IndentingNewLine]spl = splitn[sortn[data, 1], 1]; [IndentingNewLine]avs = (((Tr[#1[([All, 2])]]/ Length[#1])) &) /@ spl; [IndentingNewLine]inds = spl[([All, 1, 1])]; [IndentingNewLine]dd = Range[Max[Last[inds], max]]; [IndentingNewLine]zs = dd - dd; [IndentingNewLine]zs[([inds])] = avs; [IndentingNewLine]Transpose[{dd, zs}][IndentingNewLine]])], Input], Cell[CellGroupData[{ Cell[BoxData[ (Table[(binnedAverage2[data] // Timing) // First, {5}])], Input], Cell[BoxData[ ({2.8100000000000005` Second, 3.0700000000000003` Second, 2.799999999999999` Second, 2.799999999999999` Second, 2.75` Second})], Output] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[3. Daniel Lichtblau, Subsection, CellDingbat->None, CellMargins->{{10, Inherited}, {Inherited, Inherited}}], Cell[BoxData[ (Quit)], Input], Cell[BoxData[ ((data = ((SeedRandom[0]; Table[Random[ Integer, {1, 100}], {100000}, {2}]));))], Input], Cell[BoxData[ StyleBox[(averageByBin3[ data : {{_, _} .. }] := [IndentingNewLine]Module[n {res = Transpose[data]}, n res = averagebyBinC[res[([1])], res[([2])]]; n Transpose[{Range[Length[res]], res}]n]), FormatType->StandardForm]], Input], Cell[BoxData[ StyleBox[((n)((averagebyBinC = Compile[{{xvals, _Integer, 1}, {yvals, _Real, 1}}, n Module[{len = Max[xvals], binsizes, averages, len2 = Length[xvals], indx}, n binsizes = Table[0, {len}]; n averages = Table[0. , {len}]; n Do [nindx = xvals[([j])]; n(binsizes[([indx])]++); n averages[([indx])] += yvals[([j])], n{j, len2}n]; n binsizes = Map[If [# == 0, 1, #] &, binsizes]; n averages / binsizesn]n];)(n) )), FormatType->StandardForm]], Input], Cell[CellGroupData[{ Cell[BoxData[ (Table[((averageByBin3[data];) // Timing) // First, {5}])], Input], Cell[BoxData[ ({2.3099999999999996` Second, 2.1400000000000006` Second, 2.3599999999999994` Second, 2.3100000000000005` Second, 2.3599999999999994` Second})], Output] }, Closed]] }, Closed]] }, Closed]] }, Open ]] }, ScreenRectangle->{{0, 1024}, {0, 723}}, AutoGeneratedPackage->None, WindowToolbars->RulerBar, WindowSize->{459, 310}, WindowMargins->{{250, Automatic}, {Automatic, -20}}, ShowSelection->True ] +++++++ ==== greetings: can someone suggest a method for capturing what is printed via ?Global`* to a list? also i would like ?@ in a list. michael Reply-To: kuska@informatik.uni-leipzig.de ==== Last /@ Select[ Flatten[Messages[Evaluate[ToExpression[#]]] & /@ Names[Global`*]], MatchQ[#, RuleDelayed[HoldPattern[_], _]] &] Jens greetings: can someone suggest a method for capturing what is printed > via ?Global`* to a list? also i would like ?@ in a list. michael ==== check out the Names function. --- / FROM iMic AT 02.10.18 09:39 (Yesterday) / --- > greetings: can someone suggest a method for capturing what is printed > via ?Global`* to a list? also i would like ?@ in a list. michael -- Daniel Reeves -- http://ai.eecs.umich.edu/people/dreeves/ You know, it's at times like this when I'm trapped in a Vogon airlock with a man from Betelgeuse and about to die of asphyxiation in deep space that I really wish I'd listened to what my mother told me when I was young! Why, what did she tell you? I don't know, I didn't listen! ==== I need to color the surface of a regular polyhedron (an icosahedron, specifically) according to a relatively simple function of the spatial coordinates of the surface. I easily got a nice icosahedron of the appropriate size, but thus far I have been unable to resolve the coloring issue. Any advice is appreciated. Francis Reply-To: ng ==== Try to turn off lighting (in options: Lighting->False ). I presume you are using graphics primitives for your surfaces and have read the help file for SurfaceColor[ frontcolor, backcolor ] Then I think it should work. Nicholas > I need to color the surface of a regular polyhedron (an icosahedron, > specifically) according to a relatively simple function of the spatial > coordinates of the surface. I easily got a nice icosahedron of the > appropriate size, but thus far I have been unable to resolve the > coloring issue. Any advice is appreciated. Francis > ==== >can someone suggest a method for capturing what is printed >via ?Global`* to a list? also i would like ?@ in a list. > Names[Global`*] Names[@] Bob Hanlon ==== > Dear group, I have the following question regarding a lengthy calculation > using Mathematica: For given w points in x direction and h points in y direction, I can > construct all the points using h=10; w=8; > points=Flatten[Transpose[Outer[List,Range[w],Range[h]]],1] Next, I need to find all the possible pairs of point including points > themselves, i.e., pair AA. I can use pairs=Outer[List,points,points,1] Then, I have to clear those pairs that repeat themselves, i.e., pair > AB and > pair BA. Also, when w and h are of the order of 1000s, the computation > takes a very long time. Is there a better way of doing the second > part of Have you tried using the KSubsets function from DiscreteMath`Combinatorica`? What you want are all KSubsets of length 2 plus the original set of points duplicated for example: Needs[DiscreteMath`Combinatorica`] pairs=Join[KSubsets[points,2],Transpose[{points,points}]] Sseziwa ==== Dear All, I am currently beginning to use the 4.2 version of Mathematica on winNT. I find that the front end is very crash prone (a lot more then 4.1)... it often crashes on editing things. I am wondering if anyone experienced the same problems and if any solution had already been found. -- ________________________________________________________ Dr. Nicolas Fressengeas - - - http://www.ese-metz.fr/~fresseng Sup.8elec / Laboratoire Mat.8eriaux Optiques, Photonique et Syst.8fmes 2 rue E.Belin, 57070 METZ Cedex Plan d'acc.8fs: http://www.iti.fr/PlanPerso/23704/1 When everything else fails, read the instructions... ==== I would like to replace Dt[x_] in a complex expression. For example consider Dt[x]/.Dt[arg_]->f[arg] where x is a pure Symbol (has no value). However, Mathematica refuses apply the rule. Various other attempts, for example Unevaluted[Dt[arg_]]:>f[arg] and HoldForm[...] were unsuccessful. I managed to circumvent the problem by Dt[x]/.Dt->fun However, I were more happy if I understood why the first attempt was not successful. Can someone explain? Johannes Ludsteck <><><><><><><><><><><><> Johannes Ludsteck Economics Department University of Regensburg Universitaetsstrasse 31 93053 Regensburg ==== I'm not altogether sure I understood your problem, but it seems to me that if In[1]:= lay={{{ob11},{ob12},{ob13},{ob14}},{{ob21},{ob22},{ob23},{ob24}}}; then In[2]:= Flatten[lay,1] Out[2]= {{ob11},{ob12},{ob13},{ob14},{ob21},{ob22},{ob23},{ob24}} gets rid of the intermediate layer. Alternatively, List/@Flatten[lay]. Tomas Garza Mexico City ----- Original Message ----- f[gg] gives {{{ob11}, {ob12}, {ob13}},{{ob21}, {ob22}, {ob23}, {ob24}}, ... > ,{{obn1}, {obn2}, {obn3}}} where the number of objects in each layer can vary. The additional intermediate layer in the output of f prevents the feedback > to f, when use function like nest. I looked up and tried several method, and it seems to be easy to get rid of > all inner layers, or the innerest layer, or the outmost layer (use > Sequence). Is there a way to get rid of the middle layer as describe above? > Sincerely, > JT > _________________________________________________________________ > Surf the Web without missing calls! Get MSN Broadband. > http://resourcecenter.msn.com/access/plans/freeactivation.asp > ==== gg={{ob1},{ob2},{ob3}}; However, the output of f having a form like: f[gg] gives {{{ob11}, {ob12}, {ob13}},{{ob21}, {ob22}, {ob23}, {ob24}}, ... ,{{obn1}, {obn2}, {obn3}}} where the number of objects in each layer can vary. The additional intermediate layer in the output of f prevents the feedback to f, when use function like nest. I looked up and tried several method, and it seems to be easy to get rid of all inner layers, or the innerest layer, or the outmost layer (use Sequence). Is there a way to get rid of the middle layer as describe above? Sincerely, JT _________________________________________________________________ Surf the Web without missing calls!æGet MSN Broadband. http://resourcecenter.msn.com/access/plans/freeactivation.asp ==== > I use NDSolve to approximate a 6-dimensional and highly non-linear > dynamic system numerically. QUESTION: What are usual (and simple) techniques to check the > reliability of the resulting numerical solution? I would check the the boundary conditions and then plot the residual error from the differential equation. Here is an example with all output deleted. possibilities = NDSolve[{y[x]^4 + y'[x]^3 == 0, y[0] - 1 == 0}, y[x], {x, 0, 1}] yy[x_] = y[x] /. possibilities SHOW RESULTS Needs[Graphics`Colors`]; SetOptions[Plot, PlotStyle -> ({Thickness[0.01], #1} & ) /@ {Red,Green,Blue} ]; Plot[Evaluate[Re[yy[x]]], {x, 0, 1}] Plot[Evaluate[Im[yy[x]]], {x, 0, 1}] CHECK INITIAL CONDITIONS yy[0] - 1 PREPARE AND SHOW RESIDUAL res[x_] = yy[x]^4 + yy'[x]^3 Plot[Evaluate[Re[res[x]]], {x, 0, 1}] Plot[Evaluate[Im[res[x]]], {x, 0, 1}] Plot[Evaluate[Re[res[x]]], {x, 0, 0.001}, PlotRange -> All] Plot[Evaluate[Im[res[x]]], {x, 0, 0.001}, PlotRange -> All] ==== If there is a preserved quantity in the ODE system (and the residuals referred to a previous post on this thread certainly are) you can use it to perform an independent check. It must be noted though that the residuals are not appropriate for an entirely independent check in most cases, since most numerical solvers use them for exactly this kind of check. ODEs from physics usualy have one or more 'integrals'(preserved quantities); typicaly the energy of the system. Example: the ODEs {x'[t]==-y[t],y'[t]==x[t]} describe a rotation around the origin with fixed angular velocity. The preserved quantity is ofcourse the distance from the origin d[t]=Sqrt[x[t]^2+y[t]^2]. One could use Abs[d[0]-d[t]] as an independent check for the validity of the approximate NDSolve solution. My advice is study carefuly the ODEs before feeding them to NDSolve; analytical results (such as the existence of preserved quantities) can radicaly improve the usefulness of numerical results. Orestis ==== I apologize for the length of this post, but I don't see how else to be precise about my question. The short story is this: -- copy and paste the lines below into Mathematica -- execute -- the result is a really big expresssion, but I want the terms in this expression to be grouped in powers of my variable b. However, Collect[] doesn't appear to be working right. I call Collect[%, b] and I think I'm getting % back. At the very least, I can clearly see more than one term in the expression that includes b^9, for example. I'd be very grateful if someone a little more knowledgeable than I could execute these lines and see if they can get Collect[] to work. Of, if this is how Collect[] is supposed to work, what command should I be using? Troy. So, here's what I did. First, to get my equation: denom = Sqrt[(B^2 - r^2)^2 + 4*(r^2)*(b^2)] cnu = (2*b^2 - B^2 + r^2)/denom snu = -2*b*Sqrt[B^2 - b^2]/denom sif = 2*r*b/denom cif = (r^2 - B^2)/denom pdr = -Cos[ds]*Sin[q]*(snu*cif + cnu*sif) - Sin[ds]*(cnu*cif - snu*sif) HH = -(B^2 - b^2)*V^2/(r^2) + (((B*V)^2)/( r^2) - 2*w*b*V*Cos[q]*Cos[ds] + (w* r)^2 - (w*r*pdr)^2)*(Cos[qr])^2 Now, my equation is really HH == 0, but there's some manipulations I want to do first. I don't know Mathematica well, so all I could see to do was to perform operations on HH, then put the equation together. H2 = Expand[HH] H3 = Collect[HH, Sqrt[B^2 - b^2]] H4 = H3*( (4 * (b*r)^2 + (B^2 - r^2)^2)^2 ) H7 = H4*(r^2) H8 = Collect[ Cancel[H7], Sqrt[-b^2 + B^2] ] H9 = Equal[H8, 0] /. Equal[ aa_ + Sqrt[B^2 - b^2]*bb_, 0] -> Equal[ Sqrt[B^2 - b^2]*bb, -aa] H10 = Thread[#^2 &[H9], Equal] // ExpandAll H11 = H10 /. Equal[ mm_ , nn_] -> Equal[ mm - nn , 0] H12 = H10 /. Equal[ qq_ , 0] -> qq For H3 and H8, Collect[] seems to work. The command to get H10 I copied from a post by Andrzej Kozlowski. H13 = Collect[H11,b] H14 = Collect[H12,b] but the results don't seem to be 'collected' polynomials. ==== > I apologize for the length of this post, but I don't see how else to be > precise about my question. The short story is this: > -- copy and paste the lines below into Mathematica > -- execute > -- the result is a really big expresssion, but I want the terms in this > expression > to be grouped in powers of my variable b. However, Collect[] doesn't appear to be working right. I call Collect[%, > b] and I think I'm getting % back. At the very least, I can clearly see > more than one term in the expression that includes b^9, for example. I'd be very grateful if someone a little more knowledgeable than I could > execute these lines and see if they can get Collect[] to work. Of, if > this is how Collect[] is supposed to work, what command should I be > using? > Troy. So, here's what I did. First, to get my equation: denom = Sqrt[(B^2 - r^2)^2 + 4*(r^2)*(b^2)] > cnu = (2*b^2 - B^2 + r^2)/denom > snu = -2*b*Sqrt[B^2 - b^2]/denom > sif = 2*r*b/denom > cif = (r^2 - B^2)/denom pdr = -Cos[ds]*Sin[q]*(snu*cif + > cnu*sif) - Sin[ds]*(cnu*cif - snu*sif) HH = -(B^2 - b^2)*V^2/(r^2) + (((B*V)^2)/( > r^2) - 2*w*b*V*Cos[q]*Cos[ds] + (w* > r)^2 - (w*r*pdr)^2)*(Cos[qr])^2 Now, my equation is really HH == 0, but there's some manipulations I > want to do first. I don't know Mathematica well, so all I could see to > do was to perform operations on HH, then put the equation together. H2 = Expand[HH] > H3 = Collect[HH, Sqrt[B^2 - b^2]] > H4 = H3*( (4 * (b*r)^2 + (B^2 - r^2)^2)^2 ) > H7 = H4*(r^2) > H8 = Collect[ Cancel[H7], Sqrt[-b^2 + B^2] ] > H9 = Equal[H8, 0] /. Equal[ > aa_ + Sqrt[B^2 - b^2]*bb_, 0] -> Equal[ Sqrt[B^2 - b^2]*bb, -aa] > H10 = Thread[#^2 &[H9], Equal] // ExpandAll > H11 = H10 /. Equal[ mm_ , nn_] -> Equal[ mm - nn , 0] > H12 = H10 /. Equal[ qq_ , 0] -> qq For H3 and H8, Collect[] seems to work. > The command to get H10 I copied from a post by Andrzej Kozlowski. H13 = Collect[H11,b] > H14 = Collect[H12,b] but the results don't seem to be 'collected' polynomials. You have mistyped the definition of H12 - you want the last replacement applied to H11, not H10! As written H12 is equal to H10 (the rule does not apply). Sorry! S. Bonanos http://www.inp.demokritos.gr/~sbonano/SB.html ==== It's difficult to guess what you wanted to attain and what your problems were. Perhaps this example might help you: In[11]:= f[x_, _] /; Mod[x, 2Pi] == 0 := x/(2Pi) In[12]:= f[x_, y_] := x + y In[13]:= epsilon = $MachineEpsilon Out[13]= 2.220446049250313*^-16 In[14]:= f[6Pi(1 + epsilon Sin[#]), #] & /@ Range[0, 4Pi, Pi 10^-1] Out[14]= {3, 3., 19.4779, 19.792, 20.1062, 20.4204, 20.7345, 21.0487, 21.3628, 3., 3, 22.3053, 22.6195, 22.9336, 23.2478, 23.5619, 23.8761, 24.1903, 24.5044, 24.8186, 3, 3., 25.7611, 26.0752, 26.3894, 26.7035, 27.0177, 27.3319, 27.646, 3., 3, 28.5885, 28.9027, 29.2168, 29.531, 29.8451, 30.1593, 30.4734, 30.7876, 31.1018, 3} Look close at the values returned (compare with epsilon Sin[Range[0, 4Pi, Pi 10^-1]]), also to recognize the dangers of such a definition. -- Hartmut ==== I pretend to know how many times the function f has to be evaluated when using NDSolve to solve a differential equation. I've tried : cont=0; f[x_]:=(cont++;x^2); NDSolve[ {y''[x]+y[x]==f[x],y[0]==1,y'[0]==0},y,{x,0,6}] I obtain the solution of the ODE but cont returns 1 instead the number of funtion evaluations. How can I obtain it? Higinio Ramos ==== Hi Higinio, you need to make sure that the function f[x] is evaluated only when it gets a numeric argument. Otherwise, Mathematica expands f[x] symbolically to the expression x^2 before NDSolve is called. Writing your code as follows will solve your problem: cont=0; f[x_?NumericQ]:=(cont++;x^2); NDSolve[ {y''[x]+y[x]==f[x],y[0]==1,y'[0]==0},y,{x,0,6}] Eckhard Hennig ==== cont=0; f[x_?NumericQ]:=(cont++;x^2) NDSolve[ {y''[x]+y[x]==f[x],y[0]==1,y'[0]==0},y,{x,0,6}] -Jim ==== Higinio, Try the following: Clear[f] cont=0; f[x_?NumericQ]:=(cont++;x^2); NDSolve[ {y''[x]+y[x]==f[x],y[0]==1,y'[0]==0},y,{x,0,6}] The problem you are encountering occurs because NDSolve evaluates it's inputs first, and so NDSolve is trying to solve the differential equation y''[x]+y[x]==x^2 instead of y''[x]+y[x]==f[x]. By restricting the definition of f to evaluate only when it's input is numeric, NDSolve will keep the f[x] in the differential equation. I've included a Clear[f] statement just in case the old definition of f was still there. Carl Woll Physics Dept U of Washington Reply-To: kuska@informatik.uni-leipzig.de ==== Hi, and after cont = 0; f[x_?NumericQ] := (cont++; x^2); NDSolve[{y''[x] + y[x] == f[x], y[0] == 1, y'[0] == 0}, y, {x, 0, 6}] cont is 93. Thats because *you* gave only a blank pattern for f[x_] and NDSolve[] evaluate it to compile the right hand sides or translate it to an internal function. So *your* function is only once when Mathematica can include it into it's internal function for the right hand sides. Jens ==== I was looking for a function that can convert numbers between bases. For example, going from Base 3 to Base 9, or from Base 16 to Base 8. It seems like I have to a few steps and then use BaseForm to convert to the desired end base. Is there a single command or simple function to do such conversions? For example, ConvertBase[from_, to_, num_] would convert from base to base 9 using the number num as the number to convert (it would best be allowed to let users input in that specific base since they'd be going from base 3 to base 9, for example so maybe num is needed at all). ==== I never liked the function, but now I have need of it. I hope someone can help me see the way to implement it. For a class I am teaching I want to write a simple function that does a random search for a function's maxima. I feed the function N randomly generated inputs and look at the outputs. Provided N is large enough, the largest of the outputs I obtain should be near the function's maximum value. What I want to be able to do is to identify the input that gave me that largest output. Any ideas? My only idea is this: create a table of pairs {output, input{ and find the element of the table that has the largest first entry. This will give me the input I seek. Is there a slick way of doing this in mathematica? -- Jason Miller, Ph.D. Division of Mathematics and Computer Science Truman State University 100 East Normal St. Kirksville, MO 63501 http://vh216801.truman.edu 660.785.7430 ==== Here's one way: argMax::usage = argMax[f, domain] returns the element of domain for which f of that element is maximal -- breaks ties in favor of first occurrence.; SetAttributes[argMax, HoldFirst]; argMax[f_, dom_List] := And here's another way (defining it in terms of tupleMax)... tupleMax::usage = tupleMax[list] returns the tuple that is lexicographically maximal.; tupleMax[l_List] := Fold[If[OrderedQ[{#1, #2}], #2, #1]&, First[l], Rest[l]] argMax[f_, dom_List] := tupleMax[{f[#],#}& /@ dom][[2]] --- / FROM Jason Miller AT 02.09.18 06:18 (Today) / --- -- -- -- -- -- -- -- -- -- -- -- -- Daniel Reeves http://ai.eecs.umich.edu/people/dreeves/ In science it often happens that scientists say, 'You know that's a really good argument; my position is mistaken,' and then they would actually change their minds and you never hear that old view from them again. They really do it. It doesn't happen as often as it should, because scientists are human and change is sometimes painful. But it happens every day. I cannot recall the last time something like that happened in politics or religion. -- Carl Sagan, 1987 CSICOP Keynote Address Reply-To: kuska@informatik.uni-leipzig.de ==== Hi, fun[x_, y_] := x^2 + y^2 RandomMin[n_Integer] := First[Sort[ {#, fun @@ #} & /@ Table[{Random[], Random[]}, {n}], Last[#1] < Last[#2] &]] should do it. Jens Reply-To: kuska@informatik.uni-leipzig.de ==== Hi, does help you ? Jens ==== I've enhanced the code for PartialEvaluation (see my post, which I incorrectly attaced to someone else's reply). It is somewhat more flexible now and does more error checking. It finds the rule in DownValues that matches f[args] (allowing for more than one rule in DownValues); it allows the user to specify the position of the main CompoundExpression; and it allows the user to specify an expression that gets appended to the truncated CompoundExpression (and hence evaluated and returned). The package is just a bit too large to post here. It can be downloaded from my web site at http://www.markfisher.net/~mefisher/mma/mathematica.html Nevertheless, I can give an outline of the code here (absent the error checking stuff). PartialEvaluation[f[args], n, expr]: (* get the DownValues and turn them off *) dv = DownValues[f]; DownValues[f] = {}; (* find the rule that matches *) matches = Position[MatchQ[f[args], #]& /@ dv[[All, 1]], True]; match = dv[[ matches[[1, 1]] ]]; (* find the main CompoundExpression *) ppos = Position[match, HoldPattern[CompoundExpression[__]]]; pos = First[Sort[ppos]]; (* extract, truncate, append to, and reinsert it *) held = Extract[match, pos, Hold]; held = ReplacePart[held, Sequence, {1, 0}]; held = Take[held, n]; held = Join[held, Hold[expr]]; held = ReplacePart[Hold[Evaluate[held]], CompoundExpression, {1, 0}]; match = ReplacePart[match, held, pos, 1]; (* apply the modified rule and restore DownValues *) result = f[args] /. match; DownValues[f] = dv; result --Mark Reply-To: kuska@informatik.uni-leipzig.de ==== Hi, Clear[f] f[0, _] := q f[Pi, _] := p f[2Pi, _] := r f[x_, y_] := x^2*y^2 ??? Jens ==== There were some typos in my last message. It should have read: Let the known values of the function at x={0,2Pi,4Pi} be {f0,f2,f4}. Your conditions can be met by f(x,y)=x(x-2Pi)(x-4Pi)g(y)+h(x), where {h(0),h(2Pi),h(4Pi)}={f0,f2,f4}. An h(x) can be found by quadratic interpolation: h[x]= a + bx + cx^2; Solve[{a==f0,a+2Pi b +(2Pi)^2 c==f2, a+4Pi b+(4Pi)^2 c==f4},{a,b,c}] ==== Let the known values of the function at x={0,2Pi,4Pi} be {f0,f2,f4}. Your conditions can be met by f(x,y)=x(x-2Pi)(x-4Pi)g(y)+h(x), where {h(0),h(2Pi),h(4Pi)}={f0,f2,f4}. An h(x) can be found by quadratic interpolation: h[x]= a + bx + cx^2; Solve[{a==f0,a+2Pi b +(2Pi)^2+(2Pi)^2 c==f2, a+4Pi b+(4Pi)^2 c==f4] Maurice Reply-To: murray@math.umass.edu ==== The same principles that allow particular cases for defining a function of a single variable would apply here, because Mathematica applies particular rules before it applies general ones. For example: f[0, y_] := 0 f[2 Pi, y_] := y - 2 f[4 Pi, y_] := y - 4 f[x_, y_] := x^2 + y^3 This will do exactly what it looks like it does! If you have a general family of particular cases, say at all even integral multiples of Pi, then you could use something like the following in place of the first three lines above: f[k_ Pi, y_] := y - k /; IntegerQ[k] && EvenQ[k] There are variants as to where to place the condition IntegerQ[k] && EvenQ[k], for example: f[k_ Pi , y_] /; IntegerQ[k] && EvenQ[k] := y - k f[k_ Pi /; IntegerQ[k] && EvenQ[k], y_] := y - k -- Murray Eisenberg murray@math.umass.edu Mathematics & Statistics Dept. Lederle Graduate Research Tower phone 413 549-1020 (H) University of Massachusetts 413 545-2859 (W) 710 North Pleasant Street Amherst, MA 01375 ==== Fabio, f[x_ /; IntegerQ[Rationalize[x/(2Pi)]], y] := g[x] f[x_, y_] := h[x, y] f[# Pi, y] & /@ Range[10] {h[Pi, y], g[2*Pi], h[3*Pi, y], g[4*Pi], h[5*Pi, y], g[6*Pi], h[7*Pi, y], g[8*Pi], h[9*Pi, y], g[10*Pi]} f[2.0Pi, y] g[6.28319] f[3.5, y] h[3.5, y] f[4.01Pi, 5] h[12.5978, 5] You can use a second argument in Rationalize to control how small a rationalization error is allowed and there is always the question about what domain around 2nPi you want to include in the first definition. David Park djmp@earthlink.net http://home.earthlink.net/~djmp/ ==== Daniel, I'd say your algorithm (as an ALGORITHM) is precisely the same as Husain's. The implementation is different, but the logical steps are the same (except perhaps for <= vs. <), and so is the logical structure of the tree produced. Even the storage method is pretty much the same. Using {} instead of Null and List instead of 'node', and putting values in the middle position, doesn't add up to a change that could affect algorithmic complexity. Even the fact that Flatten gives a sorted List with your implementation is something you get for free, and the advantage comes only when it's time to output a sorted list. It has nothing to do with the work of building the tree. Algorithmic complexity (absent hidden factors under the kernel's control) is precisely the same for your solution and Husain's (mine too, I think). The speed differences we're seeing have to do with nuances of the kernel, pattern matching costs, etc. Our three nearly equivalent algorithms simply incur different hidden costs. That's the case with many of the problems that turn into speed wars among us. We compare implementation speed (including hidden costs), not pure algorithmic complexity. Compiled code or packed arrays often make all the difference, and it usually happens in the background where we can't see it. As you've hinted below, the vagaries of Update can make a linear algorithm quadratic, but may NOT so plague a slightly different implementation of the very same algorithm. So... I'm suggesting we might want to be careful when we leap to conclusions about algorithmic complexity based on a few timed experiments. Bobby Treat -----Original Message----- I had tried variations on this but they all seemed ever so slightly slower. Might be version an OS dependent. That Null that Goes by the name of tree sort. Using Flatten is basically a shortcut for walking the tree left-to-right. It seems I did not look hard enough at the original post or the follow-ups. In point of fact, regarding speed, there IS no obvious advantage of my code over that first version. The advantage is there, of course, as timing checks indicate. But it is quite far from obvious. The reason is related to the (rare) need for Update. It is pointed out in the manual that infinite evaluation is not in fact possible (no surprise thus far). A consequence is that Mathematica requires internal optimization features to keep reevaluation to a minimum. Update may be required in cases where such optimizations are overzealous. What we have in the case of this example is, roughly, the opposite. The code that determines need for reevaluation is simply not working sufficiently well. It does make the correct determination, but only after looking over the entire structure that is being evaluated. As we have a structure growing linearly in the number of iterations, and it gets mostly traversed each iteration, the algorithmic complexity becomes at least O(n^2). Not a disaster unless the appropriate complexity is smaller, which was the case for this example. This is not a bug insofar as it is a (regretably) necessary consequence of Mathematica evaluation semantics. The problems it can cause are all the same quite undesirable. We've made some inroads on this problem by enlarging substantially the class of expressions for which reevaluation may be quickly short-circuited, with the idea of fixing all but the most pathological of examples. The status of that work is unfortunately not known to me at present, though I'll try to find out about it. Daniel Lichtblau Wolfram Research ==== I'm curious, why the built-in boolean functions And and Or aren't commutative? Attributes/@{Or,And}//ColumnForm {Flat, HoldAll, OneIdentity, Protected} {Flat, HoldAll, OneIdentity, Protected} However, the arithmetical functions are Orderless: Attributes/@{Plus,Times}//ColumnForm {Flat, Listable, NumericFunction, OneIdentity, Orderless, Protected} {Flat, Listable, NumericFunction, OneIdentity, Orderless, Protected} ------------------- Evgeni Trifonov Vladivostok, Russia ==== Evgeni, Note that both And and Or can return a value without evaluating all of their arguments. For example, if the first argument to And is False, there is no reason to look at the other arguments. Suppose a user knows that one of the arguments is usually false, and so would like to look at that argument before any of the others, to avoid unneeded computations of the other arguments. He would put that argument first. If Mathematica turned around and sorted the arguments (which is what happens when a function is orderless), then that argument might end up being evaluated last. If the arguments take a significant amount of time to compute, then sorting first may cause the function to take much longer to evaluate. At any rate, if the above is not a concern for you, you may always change the attributes of And and Or to anything you want. Carl Woll Physics Dept U of Washington ----- Original Message ----- ==== I have defined two functions (Math 4.1) In[1]:= integrate[c_ y_,t_]:=c*integrate[y,t]/;FreeQ[c,t]; In[2]:= integrate[Exp[(a_.)*tau], t_]:=(E^(a*t)-1)/a/; FreeQ[a,t]&&FreeQ[a,tau]; when I applied these functions the solution is different if the coefficient is the letter d or lower (solution fine) or f or higher (solution wrong) . Here is an example: Using letter f In[3]:=integrate[f Exp[-0.8316*tau], t] Out[3]:=integrate[f, t]/E^(0.8316*tau) (*wrong*) Using letter a (the solution is rigth In[4]:=integrate[a Exp[-0.8316*tau], t] Out[4]:=-1.20250*a*(-1 + E^(-0.8316*t)) What it wrong? Guillermo Sanchez --------------------------------------------- This message was sent using Endymion MailMan. ==== Neither of your solutions are wrong according to your rules. Lets step through them In[3]:=integrate[f Exp[-0.8316*tau], t] lexographical (sp) order as In[3]:=integrate[Exp[-0.8316*tau] f, t] It then applies rule 1: integrate[c_ y_,t_]:=c*integrate[y,t]/;FreeQ[c,t] Well, The term Exp[-0.8316*tau] contains no t, so according to your rule it is equal to Exp[-0.8316*tau] integrate[f, t] At this point none of your rules apply, so the kernel stops. When you use a instead of f the ordering of the constant and the exponential are reversed, and you get the behaviour that you actually desired. A simple fix is to replace the condition in rule 1 by FreeQ[c,t]&&FreeQ[c,tau] You may however want to think carefully about exactly what you want your function integrate to do. For example, with the rules given, the variable t and the variable tau seem to be playing the same role. Erich ==== << Graphics`ImplicitPlot` does the job. Matthias Bode Sal. Oppenheim jr. & Cie. KGaA Koenigsberger Strasse 29 D-60487 Frankfurt am Main GERMANY Mobile: +49(0)172 6 74 95 77 Internet: http://www.oppenheim.de -----Ursprí.b9ngliche Nachricht----- Gesendet: Mittwoch, 18. September 2002 08:09 An: mathgroup@smc.vnet.net Betreff: Drawing an ellipse I am new with Mathematica, I have one question, I know that the sollution might be very easy, but I wasn't able to find it by now. I would like to draw an ellipse, the formula let's say is as follows: 0.09 x^2 +0.04 x y + 0.06 y^2 = 4 Maciej ==== Solve[eqn, var] and InverseFunction are designed to do what you want, in principle. Try Solve[(a*h+b*h^2) == f, h] For your eqn, however, it is impossible to solve for h since h appears as a factor in a sum *and* as an exponent. Try Solve[Exp[-h]*(a*h+b*h^2) == f, h] and study Mathematica's message. Matthias Bode Sal. Oppenheim jr. & Cie. KGaA Koenigsberger Strasse 29 D-60487 Frankfurt am Main GERMANY Mobile: +49(0)172 6 74 95 77 Internet: http://www.oppenheim.de -----Ursprí.b9ngliche Nachricht----- Gesendet: Mittwoch, 18. September 2002 08:10 An: mathgroup@smc.vnet.net Betreff: Invert a function Hi Everybody, How can I invert a funcion? For example In:=f[h_]:=Exp[-h](a*h+b*h^2); I want to invert it as h as a function of f How can I do that? Please suggest. Raj ==== This is the exact code you use? Maybe I am missing something here, but you seem to pass only the String obj to the setMathCommand method, not the content of the obj variable; er... this won't work; you are using the drawArea variable, without ever declaring or assigning something to it! pff... hm... you could try this: import com.wolfram.jlink.MathCanvas; public class My_DisplayGraphicsViaJava { public static void displayGraphics(String cmds[], MathCanvas drawArea){ for (int i = 0; i < cmds.length-1; i++){ drawArea.setMathCommand(cmds[i]); drawArea.repaintNow(); Thread.sleep(200); } } } save this to a file called My_DisplayGraphicsViaJava.java and compile it. Put the class file in a directory that is in your CLASSPATH (or look up how the CLASSPATH works in J/Link 2.0 in the J/Link documentation, its properly explained there); If all that has worked, load the class: LoadJavaClass[My_DisplayGraphicsViaJava]; then use the static method displayGraphics, by passing it the String array and the drawArea variablen (which is the MathCanvas object, you have instantiated in your Mathematica Notebook): My_DisplayGraphicsViaJava`displayGraphics[obj, drawArea]; obj being the String array (or string List in M_) that holds your Graphics Expressions (or what you want to show); Yes! ... use the Java code I proposed above, that should at least compile properly; Place the compiled class file in a directory on your CLASSPATH (consult J/Link docu if you don't know about CLASSPATH) and then try the Mathematica code again. I hope this was not too confusing; I don't have a working Mathematica here, so I can't verify, if my proposed code would really work; there migth be typos/errors in my code. murphee ==== I am trying to write an animation using JLink.....basically as follows...(using Mathematica 4.1, JLink 2.01 and java 1.4.1)... frm = createWindow[]; (* basically to create a MathFrame and a drawArea = MathCanvas *) Here is my drawing routine.... Map[(obj = drawMembrane[#]; drawArea@setMathCommand[obj]; drawArea@repaintNow[]; Pause[.0001];) &, t ] t = time list ={0,.25,etc}.... and drawMembrane is a function that uses the time value to make a graphics plot......unfortunately, this routine draws the next time plot very slowly, plus I also notice that without the Pause[] statement, Mathematica only draws the graph for the final time value... Following an example in the JLink documentation, I have written the following java routine....which I saved in the same dir as MathFrame, etc.... public class My_DisplayGraphicsViaJava { public void displayGraphics(String cmds[]) { // declare String Array for Graphics Objects for (int i = 0; i < cmds.length-1; i++) drawArea.setMathCommand(cmds[i]); drawArea.repaintNow(); Thread.sleep(200);} } and then in Mathematica I have LoadJavaClass[My_DisplayGraphicsViaJava]; My_DisplayGraphicsViaJava`displayGraphics[obj]; (* where obj is the array of Graphics plots created *) But needless to say this dont work....LoadJavaClass responds by saying Class not Found.....so what should I do to correct this? Also , on the java routine, was i supposed to add: Extends MathFrame, etc....how will my java routine know what drawArea.setMathCommand is, etc.....am I also supposed to compile my java routine? thanks....jerry blimbaum NSWC panama city, fl ==== Dear colleagues, I'm trying to plot the next exponential function Plot[Exp[1/(Abs[x-1]-Abs[x-2])],Range] If the Range is positive there arenÇt problems. However if the Range is negative Mathematica plots a strange beast... Have anyone found anything of this sort? I would apreciate to know how to resolve this problem. Juan Egea Garcia Bullas - Murcia - Espa.96a ==== The function you are plotting is constant for negative x. You might want to zoom in on the line and show you a picture of numerical noise. Erich ==== o = Options[Graphics3D]; (* save the original options *) SetOptions[Graphics3D, ...]; (* mess with them *) SetOptions[Graphics3D, Sequence@@o]; (* restore them *) --- / FROM David Sagan AT 02.09.18 06:29 (Today) / --- -- -- -- -- -- -- -- -- -- -- -- -- Daniel Reeves http://ai.eecs.umich.edu/people/dreeves/ The idea of programming in a low level language like C will seem as specialized and esoteric as programming in microcode or assembler seems today. -- Stephen Wolfram, creator of Mathematica Reply-To: kuska@informatik.uni-leipzig.de ==== Hmm, oldOptions=Options[Graphics3D]; (* do something with it *) SetOptions[Graphics3D, Sequence @@ oldOptions] works fine. Jens ==== I need to solve algebraic-differential-equations (some of index 1, maybe some of higher index as well) from within Mathematica 4.2. Is there an implementation of numerical (like DASSL) or symbolic algorithms capable of doing that around? I know that some numerical software systems which may have DASSL-routines that do what I want. Has anybody suggestions on which of these systems is easy to interface with Mathematica? Is there a C++ library that contains such routines and can be interfaced with Mathematica? Yours, Reinhard Oldenburg ==== Reinhard, have you tried the NDAESolve package that comes with the circuit analysis toolbox Analog Insydes? You can get a free evaluation version of Analog Insydes from www.analog-insydes.de. Eckhard Hennig ==== You can use All: In[1]:= mat = {{a,b,c},{d,e,f}}; In[2]:= mat[[All,2]] Out[2]= {b, e} -- Bhuvanesh, Wolfram Research. ==== Dear Colleagues, I intend to make an animation in which ball A rolls down on an inclined plane from the left whilst ball B - starting from the same height - rolls down Cosh[t]'s path from the right. x-axis is time t, y-axis is height h. Ball A is fine; ball B - which should arrive at h=0 before A - is beyond my means. Matthias Bode Sal. Oppenheim jr. & Cie. KGaA Koenigsberger Strasse 29 D-60487 Frankfurt am Main GERMANY Mobile: +49(0)172 6 74 95 77 Internet: http://www.oppenheim.de ==== Matthias, The simplest way to get the equation of motion is to set up the Lagrangian. Let's assume a 1 kg mass. Then the kinetic energy is KE = Simplify[(1/2)*(x'[t]^2 + D[Cosh[x[t]],t]^2)] and the potential energy is PE = 9.8*Cosh[x[t]] The Lagrangian is L = KE - PE and the equation of motion is diffeq = Simplify[ D[D[L, x'[t]], t] ] == Simplify[ D[L, x[t]] ] Now solve and animate ... xx[t_] = x[t]/. First[ NDSolve[{diffeq, x[0] == -1, x'[0] == 0}, x[t], {t, 0, 5}]] curve = Plot[Cosh[x], {x, -1, 1}] Do[ Show[curve, Graphics[Disk[{xx[t], Cosh[xx[t]]}, 0.025]], {t, 0, 5, 0.1}] ---- Selwyn Hollis ==== Having noticed your statement ... BEYOND MY MEANS I thing you aren't yet familiar with Lagrangian formalism. It's quite easy to derive a general equation of motion for a point mass, subjected to gravity and to moving on a curve f = y(x) (i.e. f = Cosh[#]&). 1) I'll leave re-deriving equation to you, here is what I've got (just copy paste it).: !(getEq[ f_] := [IndentingNewLine](x'')[ t] + (x')[t]^2 (2 (f')[x[t]] (f'')[x[t]])/(1 + (f' 1) 2) Next, you integrate it .: getSol[f_, x0_] := Module[{tStop}, NDSolve[{getEq[f], x[0] == x0, x'[0] == 0} , x , {t, 0, 10} ] // First ] 3) Here follows animation code, specially for linear versus cosh case, apply my initial conditions (below) makeDuo[f_, g_, x0_, fpt_] := Module[{ solf = getSol[f, x0], solg = getSol[g, x0], tf, tg, maxT, minT }, tf = solf[[1, 2, 1, 1, 2]]; tg = solg[[1, 2, 1, 1, 2]]; maxT = Max[{tf, tg}]; minT = Min[{tf, tg}]; Do[ Plot[{f@x, g@x} , {x, 0, x0} AbsolutePointSize[10], Hue[0], Point[{x[t], f@x[t]} /. solf], Hue[.6], Point[{x[t], g@x[t]} /. solg] } ] , {t, 0, minT, minT/fpt} ] ] 4) My initial conditions. In my opinion, you weren't true on this. Saying STARTING FROM THE SAME HEIGHT is not enough - you should specify x as well, thus specifing starting POINT and not just height y. makeDuo[(Cosh[1.] - 1) # &, Cosh[#] - 1 &, 1, 23] Feel free to modify the code, it should be easy. bye, Borut | Dear Colleagues, | | I intend to make an animation in which | | ball A rolls down on an inclined plane from the left whilst | | ball B - starting from the same height - rolls down Cosh[t]'s path from the | right. | | x-axis is time t, y-axis is height h. | | Ball A is fine; ball B - which should arrive at h=0 before A - is beyond my | means. | | | Matthias Bode | Sal. Oppenheim jr. & Cie. KGaA | Koenigsberger Strasse 29 | D-60487 Frankfurt am Main | GERMANY | Mobile: +49(0)172 6 74 95 77 | Internet: http://www.oppenheim.de | | ==== As I derived a generalization for a 3D parameterized curve yesterday, I'd noticed a mistake in my equation posted below, a factor '2' in expression involving x'[t]^2, should be '1'. Since this forum is of an alt. type, I've published the whole notebook at http://www2.arnes.si/~gljpoljane22/math/FallingCurve3D.nb Bye, Borut p.s. A 'fill-the-gap' riddle for those interested in physics lore. Richard Feynman once said: Science is like _ _ _, sometimes something useful comes out, but that is not the reason why we are doing it. | ... | 1) I'll leave re-deriving equation to you, here is what I've got (just copy | paste it).: | | !(getEq[ | f_] := [IndentingNewLine](x'')[ | t] + (x')[t]^2 (2 (f')[x[t]] (f'')[x[t]])/(1 + | (f' | 1) | ... ==== David, There is no command such as RestoreDefaults. But you can do the following: saveGraphics3DOptions = Options[Graphics3D]; Make your plots. SetOptions[Graphics3D, Sequence @@ saveGraphics3DOptions]; The example in the Help for SetOptions shows a similar example. My own personal preference would be to not change the Graphics3D options, but to add the option change to the plot itself. David Park djmp@earthlink.net http://home.earthlink.net/~djmp/ Sender: steve@smc.vnet.net ==== I have the following code and am trying to include a legend. The BarChart appears, but with no legend. I have followed the 'help files', but I think I'm missing something. In my code, SortTime and SessTime are value pairs. Photog is a string array containing labels. /* *Begin code here */ ShowLegend[BarChart[SortTime, SessTime, RGBColor[0.5, 0.5, 1]}, {{{RGBColor[1, 0, 0], Sort Time }, {RGBColor[0.5, 0.5, 1], Session Time }}, pete ==== How's this: maxima[f_, x_List] := Last@Sort@({f[#], #} & /@ x) maxima[Sin,Range[0,2Pi,0.1]] {0.999574,1.6} That gives you both f[x] and x at the maximum. Bobby Treat -----Original Message----- that largest output. Any ideas? My only idea is this: create a table of pairs {output, input{ and find the element of the table that has the largest first entry. This will give me the input I seek. Is there a slick way of doing this in mathematica? -- Jason Miller, Ph.D. Division of Mathematics and Computer Science Truman State University 100 East Normal St. Kirksville, MO 63501 http://vh216801.truman.edu 660.785.7430 ==== want base 12 numbers, for instance, stored in -- a List of digits? A character string? A list of digits, highest power to lowest, is easy to handle: ClearAll[convertBase] convertBase[(from_Integer)? n:{_Integer..}] /; Max[n] < from := from], to] convertBase[5, 3, {1, 2}] {2, 1} Bobby Treat -----Original Message----- example, ConvertBase[from_, to_, num_] would convert from base to base 9 using the number num as the number to convert (it would best be allowed to let users input in that specific base since they'd be going from base 3 to base 9, for example so maybe num is needed at all). ==== I believe that's not the adjacency matrix Thomas asked for. It doesn't have zeroes on the diagonal (it isn't square) and it doesn't have ones to indicate that two actors are associated with the same event. Instead, it shows connections between actors and events, which is actually more useful, as I'll demonstrate. (In addition -- though it doesn't really matter -- Jens-Peer switched the 'actors' and 'events' nomenclature within the function.) If you multiply it by its transpose, you get something else that's useful: lst = {{1, A}, {1, B}, {2, B}, {3, C}, {3, D}, {1, D}, {1, C}}; AdjacenceMatrix[lst : {{_, _} ..}] := Module[{actors, events, adj}, {actors, events} = Union /@ Transpose[lst]; adj = Table[0, {Length[actors]}, {Length[events]}]; Scan[(Part[adj, Sequence @@ #] = 1) &, lst /. MapIndexed[Rule[#1, First[#2]] &, events]]; adj] MatrixForm[a = AdjacenceMatrix[lst]] MatrixForm[b = a.Transpose[a]] Matrix 'b' records how many events two actors have in common. On the diagonal, it shows the total number of events each actor is connected to. It's easy to put zeroes on the diagonal: MatrixForm[c = b (1 - IdentityMatrix[Length[b]])] To get the originally intended incidence matrix, this works: However, I think matrices 'a' and 'b' are actually more useful, and 'a' easily leads to all the others. Bobby -----Original Message----- AdjacenceMatrix[lst : {{_, _} ..}] := Module[ {actors,events adj}, {events, actors} = Union /@ Transpose[lst]; adj = Table[0, {Length[events]}, {Length[actors]}]; Scan[(Part[adj, Sequence @@ #] = 1) &, lst /. MapIndexed[Rule[#1, First[#2]] &, actors]]; adj ] you get In[]:=AdjacenceMatrix[lst] Out[]={{1, 1, 0, 1}, {0, 1, 0, 0}, {0, 0, 1, 1}} Jens in to can't ==== Although it is not outstanding, I prefer to write your equation in the form: 9*x^2+4*x*y+6*y^2==400 In order to plot the ellipse you can use the package Graphics`ImplicitPlot` in this way: In[1]:= << Graphics`ImplicitPlot` In[2]:= ImplicitPlot[ 9*x^2 + 4*x*y + 6*y^2 == 400, {x, -7, 7}, Germ.87n Buitrago A. ----- Original Message ----- ==== first this in not an ellipse. the ellipse generic formula is : (x^2)/(a^2)+(y^2)/(b^2) = 0. I guest that you want report a graph of the solution y of the equation: 0.09 x^2 +0.04 x y + 0.06 y^2 = 4 for x varying into a given range. As you can veryfing using: Solve[0.09 x^2 + 0.04 x y + 0.06 y^2 == 4, y] the solutins are two: 0.1414213562373095*Sqrt[48.00000000000001 + 0.*x - 1.*x^2])}, 0.1414213562373095*Sqrt[48.00000000000001 + 0.*x - 1.*x^2])}} than you can construct a function: function1[x_] = 8.333333333333334*(-0.04*x - 0.1414213562373095*Sqrt[48.00000000000001 - *x^2])}, approximatively -6