mm-409
===
Subject: Re: 2D FFT for fast polynomial multiplication?
> Another question about this: how is the fourier-transformed function
> represented? Is it unnecessary to expand the range of sampling points
> (which would require an iFFT )? After all, more and more complex
> functions (polynomials of higher degree) are supposed to require more
> and more frequency sampling points, aren't they? So if I multiply two
> fourier-transformed Polynomials of the same m sampling points each,
> wouldn't I need to evaluate new sampling points, which might require an
> iFFT?
To multiply (convolve) two polynomials (arrays) of degree N and M, the
usual method is to zero-pad both arrays to degree at least N+M-1, do two
FFTs of the same order, multiply the Fourier amplitudes, then IFFT with
the same order to get your resulting polynomial of degree N+M-1.
(If N >> M or vice-versa, you can save some time by using overlap-add or
similar tricks to express as a sequence of convolutions of order
2*min(M,N).)
(Polynomial multiplication is just a linear convolution, for which there
are zillions of references online and elsewhere. Try Numerical Recipes,
to pick one at random.)
===
Subject: Re: 2D FFT for fast polynomial multiplication?
> To multiply (convolve) two polynomials (arrays) of degree N and M, the
> usual method is to zero-pad both arrays to degree at least N+M-1, do two
> FFTs of the same order, multiply the Fourier amplitudes, then IFFT with
> the same order to get your resulting polynomial of degree N+M-1.
But don't you need longer zero padding parts if you multiply more than two
factors?
===
Subject: Re: 2D FFT for fast polynomial multiplication?
>To multiply (convolve) two polynomials (arrays) of degree N and M, the
>usual method is to zero-pad both arrays to degree at least N+M-1, do two
>FFTs of the same order, multiply the Fourier amplitudes, then IFFT with
>the same order to get your resulting polynomial of degree N+M-1.
> But don't you need longer zero padding parts if you multiply more than
two
> factors?
Yes, of course. That's what he said in
Always pad to n, the size of the resulting polynomial.
===
Subject: Re: 2D FFT for fast polynomial multiplication?
> Another question about this: how is the fourier-transformed function
> represented? Is it unnecessary to expand the range of sampling points
> (which would require an iFFT )? After all, more and more complex
> functions (polynomials of higher degree) are supposed to require more
> and more frequency sampling points, aren't they? So if I multiply two
> fourier-transformed Polynomials of the same m sampling points each,
> wouldn't I need to evaluate new sampling points, which might require
> an iFFT?
> To multiply (convolve) two polynomials (arrays) of degree N and M, the
> usual method is to zero-pad both arrays to degree at least N+M-1, do two
> FFTs of the same order, multiply the Fourier amplitudes, then IFFT with
> the same order to get your resulting polynomial of degree N+M-1.
Ok. So this rules out the O(n) multiplication of fourier-transformed
representations because for padding, the polynomial must be in
power-coefficient-form (or however you'd call it), and the transform costs
O(n log n)
> (If N >> M or vice-versa, you can save some time by using overlap-add or
> similar tricks to express as a sequence of convolutions of order
> 2*min(M,N).)
But the runtime is still O(size of result * log(size of result)), isn't it?
> (Polynomial multiplication is just a linear convolution, for which there
> are zillions of references online and elsewhere. Try Numerical Recipes,
> to pick one at random.)
Great, I'll have a look at Numerical Recipes.
Best regards,
===
Subject: Re: 2D FFT for fast polynomial multiplication?
> (If N >> M or vice-versa, you can save some time by using overlap-add
> or similar tricks to express as a sequence of convolutions of order
> 2*min(M,N).)
> But the runtime is still O(size of result * log(size of result)), isn't
it?
Not quite; it's O((M+N) * log min(M,N)) in that case. So, there are
some savings, but only in the log.
===
Subject: Re: 2D FFT for fast polynomial multiplication?
>> c) If I have not only two, but k polynomials to be multiplied, what
>> is the time complexity of determining the product?
O(k n log n).
> Are you sure about the factor k? Keep in mind that n is the number of
> coefficients in the *result* - and if there are more factors to start
> with, the result is also bigger so the k is already contained in n, so
> to speak.
Right, sorry; I was thinking of cyclic convolutions (i.e. a fixed degree
n, mod z^n - 1). This also affects the answer to question (b), which is
basically subsumed by (or at least similar to) your question (c).
Suppose you have k polynomials of degree m, with n = km.
If you transform them all at once to multiply together in the Fourier
domain, you need to zero-pad them all to size n, so it is O(k n log n).
If you multiply them pairwise in the obvious way, you get degree 2m,
then degree 3m, then degree 4m, and so on until degree n. This has
complexity:
O(2m log 2m + 3m log 3m + ... + km log km)
= O(m log m * (2 + 3 + ... + k))
+ O(m * (2 log 2 + 3 log 3 + ... + k log k))
= O(k n log (n/k)) + O(k n log k)
I didn't bother working out the second term in detail, so that term is
only an upper bound (as O(..) technically means anyway). Besides, you
can use overlap-add tricks to do the convolution of a degree-m
polynomial with a degree-Nm polynomial in O(Nm log m) time, and this
gets rid of the second term entirely.
On the other hand, if you multiply them in equal-degree pairs in a tree
with depth lg k, the complexity is:
O(k/2 * 2m log 2m + k/4 * 4m log 4m + ... + 1 * km log km)
= O(km * (log 2m + log 4m + ... + log km))
= O(km * ([log m] * [lg k] + [log 2] * [1 + 2 + 3 + ... lg k])
= O(km * ([log m] * [lg k] + [log 2] * [lg k]^2)
= O(km * (log m * k) * (log k))
= O(n (log n) (log k))
which is better than the previous approach. I'm just working this out
on the fly, so hopefully I'm not making any algebra errors, but you
should be able to work it out yourself
> Or, which would suffice:
> which literature could I use as reference to justify a statement of the
> product of multiple multivariate Polynomials can be determined in O(n
> log n) time, where n is the size of the result?
theory behind the fast-multiplication-via-FFT. There are lots of
references on the web for this, so you should read up on it. Once you
understand it, the multivariate case is the same: it's just a
multivariate convolution, and the proof of the convolution theorem is
the same (just change the indices into vectors of indices).
People use the FFT for 2d convolutions all the time in image processing,
I'm sure you can find some links there.
> Err..the above thing *was* the third question, wasn't it? Why should it
> be redundant?
(You had three question marks in question (c). The third question mark
was redundant with the first.)
===
Subject: Re: learning how to use Matlab
> Hello matlab users. I am trying to figure out how to write a Matlab
> program to evaluate the function y=ln(1/1-x)
> for any user specified value of x, where x is a number <1(note that ln
> is the natural logarithm to the base e) Use an if structure to
> varify that the value passed to the program is legal. If the value of
> x is legal, calculate y(x). If not wright a suitable error message and
> quit.
> Please give me some pointers on how to write this program or even
> better if you can solve it email me, it would be better.
Try something like this (I didn't check it):
function y = f(x)
if (x >= 1.0)
disp('bad argument: x >= 1.0');
return;
end
y = log(1.0/(1.0 - x);
end
===
Subject: How to calculate the total coverage area of a few circles?
In my simulation, N circles with the same radius r are randomly
placed. Let P_i denote the center of circle i. For any i, p_i lies
within the coverage range of at least one other clicle, i.e. at least
one other circle contains p_i. How to calculate the total coverage
area of the N overlaped circle? The method should be easy to be
implemented by programming for simulation.
Any comments is welcome.
Thank you very much.
Leng Supeng
===
Subject: Pointers to Solve Modified Transport Problem using Matlab
I was taking a look at Vivekanand's problem (message copied below)
and came across the answer provided by Erwin (copied below). Was
curious to know whether we could solve the problem simply by
considering a bounded transportation problem.
If instead of having 1 source with tiered costs for transportation,
i.e. X monetary units for the first 1000 units, Y monetary units for
the next 1000 units, couldn't we formulate it as following:
We split the source into 2 seperate sources A and B. Consider the
cost of source A as X monetary units with a 1000 units upperbound. And
the cost of source B as Y monetary units with a 1000 units lowerbound
and a 2000units upperbound, and so forth.
This would translate the original transportation problem into
another bounded transportation problem, and we could solve it using a
regular bounded transportation problem algorithm.
Wouldn't this work? Erwin mentioned non-convexities, can anyone plz
explain. And wouldn't the method just described take that into
consideration and simplify the problem?
-Auro
--------------- Original Message ---------------------
I believe this will lead to non-convexities. You could
formulate this as a mixed-integer programming problem
and it could even be formulated using so-called
SOS2 variables. For a small example see:
http://www.gams.com/~erwin/interpolation.pdf
----------------------------------------------------------------
erwin@gams.com, http://www.gams.com/~erwin
----------------------------------------------------------------
> I am trying to solve a typical transportation problem with a small
> twist. The cost of transportation is not a fixed value but has a
> tiered values. Which means that instead of a fixed cost of
> transporting, the cost of transporting first 1000 units is some value
> and the cost of transporting next 1000 is another value less than the
> previous value (due to discounts and others).
> Help me understand if there is any standard algorithm to solve this in
> one pass. I can think of an algorithm that takes multiple passes for
> each tier and solves this.
===
Subject: Counter Intuitive Results: Optimising an FFT routine for Speed
me some good leads.
I present justification for a lot of the comments that drew
(constructive) criticism below. Firstly, let me summarise the feedback
that has pointed me in the right direction :-
> Steven
> 20.13: What's the best way of making my program efficient?
> A: By picking good algorithms and implementing them carefully.
Count me in for that, there was no intention or suggestion
of micro-optimisation or use of inline assembler in my post. What I am
trying to find out is what algorithm is effective and what C syntax is
more efficient, albeit only on the small Turbo-C development system
that I have.
The biggest speed gains in an FFT (and most algorithms) come from much
higher-
> level transformations. For example, our FFTW (www.fftw.org) ...
I see that I previously downloaded
FFTW 3.0.1 is the latest official version of FFTW
a while back, I am overwhelmed with volume :
458 *.c files (2.9MB)
42 *.h files
Could take me a while to assimilate all of this.
> you should use a different routine than clock()
I thought I did quite well on the Athlon /Win98, synchronising on the
20ms ticks and counting in between I could improve resolution to
600ns. Still dependent on the system time sharing on the CPU though.
> The best resolution comes from reading the CPU cycle counter directly (see
> www.fftw.org/download.html for code to do this on many CPUs).
Yes, I knew there was something called the software stopwatch,
hopefully I will come right with CYCLE.H.
I can't go and peep at this URL at the mo' as my C-FAQ download is
still stepping through each answer page. CYCLE.H has quite a bit of
code to it, I suppose I just chuck out most of the CPU specific code
and keep
_( Pentium cycle counter)
typedef unsigned long long ticks;
static __inline__ ticks getticks(void) {
ticks ret;
__asm__ __volatile__(rdtsc: =A (ret));
return ret; }
> Tristan
I followed my usual method of
QUERY GOOGLE -> READ -> FILTER -> POST
and using this I selected the four newsgroups that had hits on
C Speed Optimisation FFT
> comp.lang.oberon
===
Subject: Re: Optimizing array bounds checking
> comp.lang.functional
===
Subject: Re: Academia isn't interested...
....
But what you *can* do is generate low-level code from high-level
specifications. This has already been well-demonstrated in e.g. the
FFTW
(written in OCaml), which generates extremely efficient C-algorithms
for the FFT from specifications and specialises them for the given
architecture. This, however, is high-level programming again, not
low-level programming. There is no escape from imperativeness if you
work with von-Neumann architectures on a low level.
> Eric
> You'd have to examine the actual generated code to have a hope
> of figuring out what's really happening in any particular case.
Don't want to do that, anyway my solution has to be fairly portable,
from me (design engineer) it goes to software engineer, rebuilt under
Visual C++, and then targetted to TMS320 (or equiv) DSP.
> Steven
> You realize, of course, that your first mistake was starting with the
> radix-2 Numerical Recipes in C code and spending your time attempting
> micro-optimizations. The biggest speed gains in an FFT (and most
> algorithms) come from much higher-level transformations.
Actually I started with someone else's Fortran->C implementation of an
algorithm that looks like it came from the Cooley-Tukey paper. Then I
stepped to the Numerical Recipes (C) algorithm (Danielson - Lanczos?)
(which by the way still carries the legacy of Fortran array indexing)
This took -72% off the time.
I have refined this in about five major steps, so that at the point of
posting I was
-23% off the time of the NumRep routine as published (only high
level tweaks).
I will perform my last trial of converting the original code to double
without changing the indexing, just to see what I would have got,
before ...
Now, however, I will trawl through all the source in FFTW and see what
I can find.
> If you just want a fast FFT and are not doing this as a learning
experience, you'll
> be much better off downloading an existing optimized code.
Even though I will probably do better with an FFTW routine, the
lessons learned will be valuable for other pieces of time critical
code.
> A decent compiler should transform array references a[j] in a loop
> into separate incremented pointers for you, if it's advantageous to do
> so. By being clever you can easily just confuse the compiler's
> optimizer unless you know what you're doing.
The two sets of indices do not change sequentially, which is why I
hoped that me forcing the use of pointers would be faster.
> Kurt
> I've done the above many a time. The only way I've been able to cope is
> by learning the instruction set and reading the (equivalent) assembly
code.
> It's pretty easy to tell when you've confused the optimizer.
This would be interesting given unlimited time, I do know that the DSP
code was definitely hand tweaked in the past as portability and
maintainability were far lower priority than performance.
Paul
===
Subject: Re: Counter Intuitive Results: Optimising an FFT routine for Speed
postmaster@pwi.mailshell.com says...
> The best resolution comes from reading the CPU cycle counter directly
(see
> www.fftw.org/download.html for code to do this on many CPUs).
> Yes, I knew there was something called the software stopwatch,
> hopefully I will come right with CYCLE.H.
> I can't go and peep at this URL at the mo' as my C-FAQ download is
> still stepping through each answer page. CYCLE.H has quite a bit of
> code to it, I suppose I just chuck out most of the CPU specific code
> and keep
> _( Pentium cycle counter)
> typedef unsigned long long ticks;
> static __inline__ ticks getticks(void) {
> ticks ret;
> __asm__ __volatile__(rdtsc: =A (ret));
> return ret; }
This, and most of this thread is OT here, but since you've already
posted this, I don't want others to be bit by this in the future.
Danger: The above is TOTALLY UNRELIABLE on SMP systems. Getting
dependent upon this will bite you, as even notebook computers
are now available with dual CPUs. There is no synchronization
between the TSCs on the individual CPUs. So, you are as likely
as not to get bogus time values if your process gets bounced
between processors.
===
Subject: Re: Counter Intuitive Results: Optimising an FFT routine for Speed
[...]
> http://www.eskimo.com/~scs/C-faq/top.html
> through all the URLs at this site.
The other versions link,
,
points you to a compressed plain-text copy of the whole thing at
. You can uncompress it
with uncompress or gunzip. Apparently it's the most up-to-date
version available.
===
Subject: Re: Counter Intuitive Results: Optimising an FFT routine for Speed
> Don't want to do that, anyway my solution has to be fairly portable,
> from me (design engineer) it goes to software engineer, rebuilt under
> Visual C++, and then targetted to TMS320 (or equiv) DSP.
Whoa, whoa, hold on, you want this to be fast on a DSP chip? But you're
planning to benchmark and optimize it on a general-purpose CPU?
Sigh...
===
Subject: Re: Counter Intuitive Results: Optimising an FFT routine for Speed
NNTP-Proxy-Relay: library2.airnews.net
> tempr=wr*data[ixj]-wi*data[ixj+1];
> tempi=wr*data[ixj+1]+wi*data[ixj];
> data[ixj]=data[ixData]-tempr;
> data[ixj+1]=data[ixData+1]-tempi;
> data[ixData] += tempr;
> data[ixData+1] += tempi;
It has been quite a few years since I last looked at FFT code, and that was
a casual look during some self-education time. This is the fundamental
radix-2 butterfly, isn't it?
===
Subject: Re: Counter Intuitive Results: Optimising an FFT routine for Speed
> tempr=wr*data[ixj]-wi*data[ixj+1];
> tempi=wr*data[ixj+1]+wi*data[ixj];
> data[ixj]=data[ixData]-tempr;
> data[ixj+1]=data[ixData+1]-tempi;
> data[ixData] += tempr;
> data[ixData+1] += tempi;
> It has been quite a few years since I last looked at FFT code, and that
was
> a casual look during some self-education time. This is the fundamental
> radix-2 butterfly, isn't it?
Looks like it, and programmed in a way that is guaranteed to run slow by
preventing optimisations.
===
Subject: Re: Counter Intuitive Results: Optimising an FFT routine for Speed
> Why should accessing vectorised data using pointer references be so
> much slower than indexing - surely Turbo-C cannot have pointer
> de-referencing so badly wrong?
x86 has a shortage of registers. By forcing the compiler to store
intermediate results in pointers, you might be running it out of registers
and causing thrashing of the register file.
Just a guess - you'd need to examine the assembly to see where it is really
going wrong.
If you need to optimise code, try not to use an old compiler for a new
target. It can't take advantage of any changes to the architecture since
when it was written.
===
Subject: Re: Counter Intuitive Results: Optimising an FFT routine for Speed
<3f82fb30$0$574$b45e6eb0@senator-bedfellow.mit.eduDon't want to do
that, anyway my solution has to be fairly portable,
> from me (design engineer) it goes to software engineer, rebuilt
> under Visual C++, and then targetted to TMS320 (or equiv) DSP.
> Whoa, whoa, hold on, you want this to be fast on a DSP chip? But
> you're planning to benchmark and optimize it on a general-purpose CPU?
> Sigh...
Optimising for the TMS320 processors (they vary a *lot* BTW) is *way* OT
for this group since on the ones I've used there are all sorts of trick
you can play in assembler that you can't (and the compiler I used did
not) use.
The OP should find a group or mailing list specific to the processor he
wants to use and probably obtain a library written in assembler to make
full use of the facilities the processor has. Also read the processor
manual, ISTR the TMS320C1x/2x/5x books had example FFT code...
===
Mark Gordon
Paid to be a Geek & a Senior Software Developer
Although my email address says spamtrap, it is real and I read it.
===
Subject: Counter Intuitive Results: Optimising an FFT routine for Speed
OK - I feel somewhat exonerated.
I tried various combinations of options with the high iteration code
segment. Indeed the pointerized version DID RUN SIGNIFICANTLY FASTER
(as in 21%) when I revert to 32 bit floating point here instead of
double (for those entering the thread here I repeat the salient code).
This is strange but credible, I understood that all internal
arithmetic was carried out with doubles, and casts to and from float
performed as necessary to match the variables data types. This does
not appear to be the case with Turbo-C
Now I can put my present trials to rest and start hunting through
FFTW.
BTW - I could not get CYCLE.h to work, Turbo-C does not support
__INLINE__ or __ASM__. Good thing I spent the time to get my
****PORTABLE**** timer to work.
/*
_1 curReal = wr *data[ ixj ] -wi *data[ ixj
+1];
_1 curImag = wr *data[ ixj+1 ] +wi *data[ ixj
];
_1 data[ ixj] = data[ ixData ] -curReal;
_1 data[ ixj+1] = data[ ixData+1] -curImag;
_1 data[ ixData] += curReal;
_1 data[ ixData+1] += curImag;
*/
#if defined(FFT_FLOAT)
# define srcPtr0_ floatPtr_1__
# define srcPtr1_ floatPtr_2__
# define tgtPtr0_ floatPtr_3__
# define tgtPtr1_ floatPtr_4__
#endif
srcPtr1_ = srcPtr0_ = &( data[ixj]); ++srcPtr1_;
tgtPtr1_ = tgtPtr0_ = &( data[ixData]); ++tgtPtr1_;
curReal = wr **srcPtr0_ -wi **srcPtr1_;
curImag= wr **srcPtr1_ +wi **srcPtr0_;
*srcPtr0_ = *tgtPtr0_ -curReal;
*srcPtr1_ = *tgtPtr1_ -curImag;
*tgtPtr0_ += curReal;
*tgtPtr1_ += curImag;
# undef srcPtr0_
# undef srcPtr1_
# undef tgtPtr0_
# undef tgtPtr1_
This is a typical execution time with full double variables
FFT_NRC starting to execute 5000 loops
FFT_NRC time = 923647 ±88(3.5s) ns
Here is the comparison between 32 bit floats, FFT_NRC_prev is the _1
code commented out above and FFT_NRC is the pointerized code :
FFT_NRC starting to execute 5000 loops
FFT_NRC time = 519175 ±88(3.5s) ns
FFT_NRC_prev starting to execute 5000 loops
FFT_NRC_prev time = 658008 ±88(3.5s) ns
FFT_NRC time change :FFT_NRC_Prev = -138.833 ± 0.124(3.5s)
[Micro]s
(-21.10%)
FFT_NRC_prev starting to execute 1000 loops
FFT_NRC_prev time = 654509 ±451(3.5s) ns
FFT_NRC starting to execute 1000 loops
FFT_NRC time = 515313 ±451(3.5s) ns
===
Subject: Re: Counter Intuitive Results: Optimising an FFT routine for Speed
> This is strange but credible, I understood that all internal
> arithmetic was carried out with doubles, and casts to and from float
> performed as necessary to match the variables data types. This does
> not appear to be the case with Turbo-C
The x86 registers are extended precision (or at least double, depending
upon the mode), so your compiler doesn't have much choice in the matter.
Nor does this affect the speed.
Single precision is almost always faster than double precision, though,
at least for significant transform sizes, simply because of cache usage.
(Double precision variables also need to be 8-byte aligned on x86 or the
speed dies horribly and unpredictably; not all compilers do this.)
> BTW - I could not get CYCLE.h to work, Turbo-C does not support
> __INLINE__ or __ASM__. Good thing I spent the time to get my
> ****PORTABLE**** timer to work.
You asked for something with more resolution than clock(), so I told you
the best resolution you can get; if you're happy with clock(), why did
you ask? Unfortunately, accessing the cycle counter is necessarily
compiler/CPU-specific. (You also inexplicably grabbed the code from
cycle.h labeled as gcc-specific.)
BTW, clock() may be portable, but any assumptions about its resolution
are not.
...
Not that any of this is relevant to the DSP chip you eventually plan to
run on.
===
Subject: Re: 3-d integration routine is needed
Thank you for all your suggestions, I am going to work on these suggestions
and try to come up with a solution. By the way my integration limits are
all
scalar which makes the problem even simpler.
> google but I could not find anything. Do you have any suggestions? Any
> reliable routine which involves c/c++/fortran is good enough for me.
===
Subject: Re: 3-d integration routine is needed
> google but I could not find anything. Do you have any suggestions? Any
> reliable routine which involves c/c++/fortran is good enough for me.
http://www.mat.univie.ac.at/~neum/software.html#integral
offers a number of choices. [~ is a tilde]
===
Subject: Re: 3-d integration routine is needed
> google but I could not find anything. Do you have any suggestions? Any
> reliable routine which involves c/c++/fortran is good enough for me.
in case your integral is a convolution, you may apply FFTs threetimes
for fast integrations.
Axel
===
Subject: Re: 3-d integration routine is needed
Quick and dirty solution:
Use any one dimensional integration routine recursively. Look at netlib
for a reliable one.
Mike
> google but I could not find anything. Do you have any suggestions? Any
> reliable routine which involves c/c++/fortran is good enough for me.
===
Subject: LSODAR
Has anyone had any success trying to solve two independent sets of
ODE's within a single program using LSODAR?
Reading the comments in the LSODAR I can see that with this kind of
problem the user must call the DSRCAR routine to restore and save
LSODAR's common block before and after the LSODAR call for each ODE
set (ie each problem requires the preservation of its own LSODAR
common block)
I have implemented the calls to DSRCAR in my program but could not get
sensible results. So as a test case I am trying to solve the same
simple set of ODE's twice, and should therefore get two sets of
identical results. At the moment only the first set gives me the
correct answer. If I only solve for one or the other ODE sets, I get
sensible results. It just wont work with two sets.
The only way I have managed to get round this problem is to set istate
=1 and t =0 on every call and force lsodar to re-initialise itself
each time it is called. ie its now solving many problems for time t =
0 to tout, but using updated the y array values from the provious
solution. I am then incrementing time externally from lsodar. This is
obviously not very efficient.
Does anyone have any ideas as to why I might not be able to solve the
problem by using the DSRCAR subroutine?
I am using the 8th of August 2001 release of DLSODAR.
I am currently NOT providing a jacobian routine. This might be
important as the two sets of results only start to differ once the
stiff method is used (I think!).
===
Subject: Re: LSODAR
> Has anyone had any success trying to solve two independent sets of
> ODE's within a single program using LSODAR?
> Reading the comments in the LSODAR I can see that with this kind of
> problem the user must call the DSRCAR routine to restore and save
> LSODAR's common block before and after the LSODAR call for each ODE
> set (ie each problem requires the preservation of its own LSODAR
> common block)
> I have implemented the calls to DSRCAR in my program but could not
get
> sensible results. So as a test case I am trying to solve the same
> simple set of ODE's twice, and should therefore get two sets of
> identical results. At the moment only the first set gives me the
> correct answer. If I only solve for one or the other ODE sets, I get
> sensible results. It just wont work with two sets.
> The only way I have managed to get round this problem is to set
istate
> =1 and t =0 on every call and force lsodar to re-initialise itself
> each time it is called. ie its now solving many problems for time t
=
> 0 to tout, but using updated the y array values from the provious
> solution. I am then incrementing time externally from lsodar. This
is
> obviously not very efficient.
> Does anyone have any ideas as to why I might not be able to solve
the
> problem by using the DSRCAR subroutine?
> I am using the 8th of August 2001 release of DLSODAR.
> I am currently NOT providing a jacobian routine. This might be
> important as the two sets of results only start to differ once the
> stiff method is used (I think!).
>
Try
RSAV = 0.d0
ISAV = 0
just before calling SRCAR.
===
===
Subject: Re: LSODAR
> Has anyone had any success trying to solve two independent sets of
> ODE's within a single program using LSODAR?
Reading the comments in the LSODAR I can see that with this kind of
> problem the user must call the DSRCAR routine to restore and save
> LSODAR's common block before and after the LSODAR call for each ODE
> set (ie each problem requires the preservation of its own LSODAR
> common block)
I have implemented the calls to DSRCAR in my program but could not
> get
> sensible results. So as a test case I am trying to solve the same
> simple set of ODE's twice, and should therefore get two sets of
> identical results. At the moment only the first set gives me the
> correct answer. If I only solve for one or the other ODE sets, I get
> sensible results. It just wont work with two sets.
The only way I have managed to get round this problem is to set
> istate
> =1 and t =0 on every call and force lsodar to re-initialise itself
> each time it is called. ie its now solving many problems for time t
> =
> 0 to tout, but using updated the y array values from the provious
> solution. I am then incrementing time externally from lsodar. This
> is
> obviously not very efficient.
Does anyone have any ideas as to why I might not be able to solve
> the
> problem by using the DSRCAR subroutine?
I am using the 8th of August 2001 release of DLSODAR.
> I am currently NOT providing a jacobian routine. This might be
> important as the two sets of results only start to differ once the
> stiff method is used (I think!).
Try
> RSAV = 0.d0
> ISAV = 0
> just before calling SRCAR.
Gerry,
Could you please expand on what you mean by this. RSAV and ISAV are
arrays. Do you mean set their contents to zero before the first call
to SRCAR? If so why ? - In my first call of SRCAR, JOB = 1 .Anything
in ISAV and RSAV will be overwritten. Setting them to zero at any
other time will surely result in the loss of the common block data.
Have you successfully used the SRCAR routine before? I was assuming
the simulation wasnt working because there is some additional common
data which this routine does not backup.
===
Subject: Re: LSODAR
> Has anyone had any success trying to solve two independent sets
of
> ODE's within a single program using LSODAR?
Reading the comments in the LSODAR I can see that with this kind
of
> problem the user must call the DSRCAR routine to restore and
save
> LSODAR's common block before and after the LSODAR call for each
ODE
> set (ie each problem requires the preservation of its own LSODAR
> common block)
I have implemented the calls to DSRCAR in my program but could
not
> get
> sensible results. So as a test case I am trying to solve the
same
> simple set of ODE's twice, and should therefore get two sets of
> identical results. At the moment only the first set gives me the
> correct answer. If I only solve for one or the other ODE sets, I
get
> sensible results. It just wont work with two sets.
The only way I have managed to get round this problem is to set
> istate
> =1 and t =0 on every call and force lsodar to re-initialise
itself
> each time it is called. ie its now solving many problems for
time t
> =
> 0 to tout, but using updated the y array values from the
provious
> solution. I am then incrementing time externally from lsodar.
This
> is
> obviously not very efficient.
Does anyone have any ideas as to why I might not be able to
solve
> the
> problem by using the DSRCAR subroutine?
I am using the 8th of August 2001 release of DLSODAR.
> I am currently NOT providing a jacobian routine. This might be
> important as the two sets of results only start to differ once
the
> stiff method is used (I think!).
>
Try
> RSAV = 0.d0
> ISAV = 0
> just before calling SRCAR.
> Gerry,
> Could you please expand on what you mean by this. RSAV and ISAV are
> arrays. Do you mean set their contents to zero before the first call
> to SRCAR? If so why ? - In my first call of SRCAR, JOB = 1 .Anything
> in ISAV and RSAV will be overwritten. Setting them to zero at any
> other time will surely result in the loss of the common block data.
> Have you successfully used the SRCAR routine before? I was assuming
> the simulation wasnt working because there is some additional common
> data which this routine does not backup.
,
My previous post was confusing as it related to ODEPACK prior to the
June/2001 revision which had no SRCAR and JOB and one had to 'RESTORE,
FROM RSAV AND ISAV, THE CONTENTS OF THE INTERNAL COMMON BLOCKS USED BY
LSODA. ' Now it's easier.
The only advise I can offer is that:
if IOPT<>0 , reinitialize rwork and iwork, ie,
rwork = 0.d0 !F95 syntax, actually only 5-7 need to be zeroed.
iwork = 0 !F95 syntax, actually only 5-9 need to be zeroed.
If rwork/iwork are not reinitialized for each ivp you have no
assurance that they are zero and, consequently, you might not get
LSODA(R)'s default values for inputs other than the ones you elect to
override (eg, iwork(5)=1 for verbose messages or iwork(6) to change
the default number of steps the solver can do per step). If IOPT=1, a
value of zero for any of these optional inputs will cause the default
value to be used;
and for successive ivp's, reinitialize t and istate, ie,
t = 0.d0
istate = 1
Let me know if this was of any help.
===
===
Subject: A simple question?
Content-transfer-encoding: 8bit
While studying Carmichael numbers I came across this curiosity:
Why is it that for any prime P, P > 29 that P^2 = K (mod (29^2-1))
where K results always in one of (1, 11^2, 13^2, 17^2, 19^2, 23^2).
This is true for all tested primes, but not for all composites.
It appears that aprox. 13% of carmichael numbers fail on this simple
test up to 10^12.
P^
===
#############
Imagination is more important than knowledge - A. Einstein
===
Subject: Re: A simple question?
> While studying Carmichael numbers I came across this curiosity:
> Why is it that for any prime P, P > 29 that P^2 = K (mod (29^2-1))
> where K results always in one of (1, 11^2, 13^2, 17^2, 19^2, 23^2).
> This is true for all tested primes, but not for all composites.
> It appears that aprox. 13% of carmichael numbers fail on this simple
> test up to 10^12.
I'm not sure why this is in sci.math.num-analysis, but:
The same thing happens for any number that has no common
factor with 29^2-1. You can verify this by checking just
the numbers from 0 to 29^2-2. Of course any prime > 29
satisfies that condition :-). So it appears that approximately
13% of Carmichael numbers up to 10^12 have 2, 3, 5, or 7
as factors.
So why are there so few squares mod 29^2-1? Well, it's
2^3.3.5.7; mod those factors there are, respectively,
1,1,2,3 nonzero squares. So there are only 6 squares
mod 29^2-1 (other than the non-coprime ones), and it's
not all that surprising that they're represented by
1 and 5 squares of small primes.
===
Gareth McCaughan
.sig under construc
===
Subject: who to invert 3*3 singular matrix
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision:
1.9 primary) id h973a5P15817;
I will be very appreciated if you could invert the following 3*3
matrix;
fist row entries are (2.25,.433,0)
2nd row entries are (.433,.75,0)
3rd entries are (0,0,0)
thanks alot
===
Subject: Re: who to invert 3*3 singular matrix
X-Enigmail-Version: 0.76.1.0
X-Enigmail-Supports: pgp-inline, pgp-mime
> I will be very appreciated if you could invert the following 3*3
> matrix;
> fist row entries are (2.25,.433,0)
> 2nd row entries are (.433,.75,0)
> 3rd entries are (0,0,0)
alot
This one is a possibility, the Moore-Penrose pseudo-inverse:
(0.5, -0.2889, 0)
(-0.2889, 1.5, 0)
(0, 0, 0)
Dan
===
Subject: Re: who to invert 3*3 singular matrix
> I will be very appreciated if you could invert the following 3*3
> matrix;
> fist row entries are (2.25,.433,0)
> 2nd row entries are (.433,.75,0)
> 3rd entries are (0,0,0)
alot
Why do you think it is called singular?
===
^^^^^^^^^^^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/
===
Subject: Re: who to invert 3*3 singular matrix
I will be very appreciated if you could invert the following 3*3
>matrix;
>fist row entries are (2.25,.433,0)
> 2nd row entries are (.433,.75,0)
> 3rd entries are (0,0,0)
>thanks alot
>
never try the impossible.
look in your textbook.
===
Subject: Need Your Expertise Please!
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision:
1.9 primary) id h97LQCi23500;
Hi Everyone!
I have an odd request, here goes...
I purchased an Expedition with the keypad entry - no one knows what
the combination is and I don't have time or money to take it to the
dealer...
I'd like to try combinations at my leisure until I find the one that
works. What I need is to know every possible combination using the
numbers 1,3,5,7,9. Numbers may repeat, so that needs to be taken into
consideration as well...
My email is bblocksmith@earthlink.net
Sincerely,
Jeanne
===
Subject: Re: Need Your Expertise Please!
> I have an odd request, here goes...
> I purchased an Expedition with the keypad entry - no one knows what
> the combination is and I don't have time or money to take it to the
> dealer...
> I'd like to try combinations at my leisure until I find the one that
> works. What I need is to know every possible combination using the
> numbers 1,3,5,7,9. Numbers may repeat, so that needs to be taken into
> consideration as well...
> My email is bblocksmith@earthlink.net
> Sincerely,
> Jeanne
How many numbers long is the combination? If you have a set
of 5 numbers allowing repeats, the number of possible combos
is 5^n where n is the length of the combination sequence.
However it is hard to believe someone who can afford an Expedition
can't afford to have a dealer read out the combination.
===
^^^^^^^^^^^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/
===
Subject: Re: Open source Conjugate Gradient
This solves the normal equation AT.A.x = AT.b.
/*-----------------------------------------------------
! Conjugate Gradient Method for a Sparse Linear System
!------------------------------------------------------
! Ref.: Numerical Recipes, Cambridge University Press
! 1986 (chapter 2.10).
!
! C++ demo1 version by J-P Moreau, Paris.
!-------------------------------------------------------
!
! Solve sparse linear system (size=10):
!
! 2x1 -x10 = 0
! 2x2 -x9 -x10 = 0
! 2x3 -x8 -x9 = 0
! 2x4 -x7 -x6 = 0
! 2x5 -x6 -x7 = 0
! -x5 +2x6 = 11
! -x4 -x5 +2x7 = 0
! -x3 -x4 +2x8 = 0
! -x2 -x3 +2x9 = 0
! -x1 -x2 +2x10 = 0
!
!
! Structure of sparse matrix A:
!
! X X
! X XX
! X XX
! X XX
! XXX
! XX
! XX X
! XX X
! XX X
! XX X
!
! Any key to continue...
!
! Number of iterations: 10
!
! Solution vector:
!
! 1.0000 3.0000 5.0000 7.0000 9.0000
! 10.0000 8.0000 6.0000 4.0000 2.0000
!
! RSQ = 1.386462e-018
!
!
! Product A.X:
!
! 0.0000 0.0000 0.0000 0.0000 0.0000
! 11.0000 0.0000 0.0000 0.0000 0.0000
!
!-----------------------------------------------------*/
#include
#include
#define SIZE 11
#define EPS 1e-6
//global variables
double A[SIZE][SIZE];
double B[SIZE], B1[SIZE], X[SIZE];
double RSQ;
int i,j,N;
void ASUB(double *,double *);
void ATSUB(double *,double *);
void SPARSE(double *B,int M,double *X,double *RSQ) {
/*------------------------------------------------------------------
!Solves the linear system A.x = b for the vector X of length N,
!given the right-hand vector B, and given two subroutines,
!ASUB(XIN,XOUT) and ATSUB(XIN,XOUT), which respectively calculate
!A.x and AT.x (AT for Transpose of A) for x given as first argument,
!returning the result in their second argument.
!These subroutines should take every advantage of the sparseness
!of the matrix A. On input, X should be set to a first guess of the
!desired solution (all zero components is fine). On output, X is
!the solution vector, and RSQ is the sum of the squares of the
!components of the residual vector A.x - b. If this is not small,
!then the matrix is numerically singular and the solution represents
!a least-squares approximation.
!-----------------------------------------------------------------*/
//LABEL: L1
double G[SIZE],H[SIZE],XI[SIZE],XJ[SIZE];
int IRST,ITER,J;
EPS2=N*EPS*EPS;
IRST=0;
L1:IRST++;
ASUB(X,XI);
RP=0.0;
BSQ=0.0;
for (J=1; J<=M; J++) {
BSQ=BSQ+B[J]*B[J];
XI[J]=XI[J]-B[J];
RP=RP+XI[J]*XI[J];
}
ATSUB(XI,G);
for (J=1; J<=M; J++) {
G[J]=-G[J];
H[J]=G[J];
}
for (ITER=1; ITER<=10*M; ITER++) { //main loop
ASUB(H,XI);
ANUM=0.0;
ADEN=0.0;
for (J=1; J<=M; J++) {
ANUM=ANUM+G[J]*H[J];
ADEN=ADEN+XI[J]*XI[J];
}
if (ADEN==0.0) printf( Very singular matrix.n);
ANUM=ANUM/ADEN;
for (J=1; J<=M; J++) {
XI[J]=X[J];
X[J]+=ANUM*H[J];
}
ASUB(X,XJ);
*RSQ=0.0;
for (J=1; J<=M; J++) {
XJ[J]=XJ[J]-B[J];
*RSQ=*RSQ+XJ[J]*XJ[J];
}
if (*RSQ==RP || *RSQ<=BSQ*EPS2) {
printf( Number of iterations: %dn,ITER);
return; //normal return
}
if (*RSQ > RP) {
for (J=1; J<=M; J++) X[J]=XI[J];
if (IRST>=3) return; //return if too many restarts
goto L1;
}
RP=*RSQ;
ATSUB(XJ,XI); //compute gradient for next iteration
GG=0.0;
DGG=0.0;
for (J=1; J<=M; J++) {
GG=GG+G[J]*G[J];
DGG=DGG+(XI[J]+G[J])*XI[J];
}
if (GG==0.0) {
printf( Number of iterations: %dn, ITER);
return; //rare but normal return
}
for (J=1; J<=M; J++) {
G[J]=-XI[J];
}
} //main loop
printf( Too many iterations.nn);
}
/*In these versions of ASUB1 and ATSUB1,
we do not take any advantage of
sparseness! On the contrary,
ASUB and ATSUB take into account the
sparseness of the matrix A.
N.B. For big values of N, the difference
in computation time may be important. */
void ASUB1(double *X,double *V) {
int I,J;
for (I=1; I<=N; I++) {
V[I]=0.0;
for (J=1; J<=N; J++)
V[I]=V[I]+A[I][J]*X[J];
}
}
void ASUB(double *X,double *V) {
int I;
V[1]=A[1][1]*X[1]+A[1][10]*X[10];
for (I=2; I<6; I++)
V[I]=A[I][I]*X[I]+A[I][N-I+1]*X[N-I+1]+A[I][N-I+2]*X[N-I+2];
V[6]=A[6][6]*X[6]+A[6][5]*X[5];
for (I=7; I<11; I++)
V[I]=A[I][I]*X[I]+A[I][N-I+1]*X[N-I+1]+A[I][N-I+2]*X[N-I+2];
}
void ATSUB1(double *X,double *V) {
int I,J;
for (I=1; I<=N; I++) {
V[I]=0.0;
for (J=1; J<=N; J++)
V[I]=V[I]+A[J][I]*X[J];
}
}
void ATSUB(double *X,double *V) {
int I;
V[1]=A[1][1]*X[1]+A[10][1]*X[10];
for (I=2; I<6; I++)
V[I]=A[I][I]*X[I]+A[N-I+1][I]*X[N-I+1]+A[N-I+2][I]*X[N-I+2];
V[6]=A[6][6]*X[6]+A[5][6]*X[5];
for (I=7; I<11; I++)
V[I]=A[I][I]*X[I]+A[N-I+1][I]*X[N-I+1]+A[N-I+2][I]*X[N-I+2];
}
// show structure of sparse matrix A
void ShowMat() {
int i,j;
printf(n);
for (i=1; i<=N; i++) {
printf( );
for (j=1; j<=N; j++)
if (A[i][j]!=0.0) printf(X); else printf( );
printf(n);
}
printf(n);
}
void main() {
// matrix A
N=10;
for (i=1; i<=N; i++)
for (j=1; j<=N; j++)
A[i][j]=0.0;
A[1][1]=2.0; A[1][10]=-1.0;
A[2][2]=2.0; A[2][9]=-1.0; A[2][10]=-1.0;
A[3][3]=2.0; A[3][8]=-1.0; A[3][9]=-1.0;
A[4][4]=2.0; A[4][7]=-1.0; A[4][8]=-1.0;
A[5][5]=2.0; A[5][6]=-1.0; A[5][7]=-1.0;
A[6][5]=-1.0; A[6][6]=2.0;
A[7][4]=-1.0; A[7][5]=-1.0; A[7][7]=2.0;
A[8][3]=-1.0; A[8][4]=-1.0; A[8][8]=2.0;
A[9][2]=-1.0; A[9][3]=-1.0; A[9][9]=2.0;
A[10][1]=-1.0; A[10][2]=-1.0; A[10][10]=2.0;
// vector B
for (i=1; i<=N; i++) B[i]=0.0; B[6]=11.0;
// initial guess for X
for (i=1; i<=N; i++) X[i]=0.0;
printf(n Structure of sparse matrix A:n);
ShowMat();
printf( Any key to continue... );
getchar();
printf(n);
SPARSE(B,N,X,&RSQ);
printf(n Solution vector:nn);
for (i=1; i<6; i++) printf(%8.4f,X[i]); printf(n);
for (i=6; i<=N; i++) printf(%8.4f,X[i]); printf(n);
printf(n RSQ = %enn, RSQ);
// verification A.x = B
ASUB(X,B1);
printf(n Product A.X:nn);
for (i=1; i<6; i++) printf(%8.4f,B1[i]); printf(n);
for (i=6; i<=N; i++) printf(%8.4f,B1[i]); printf(nn);
}
// end of file tsparse1.cpp
===
Subject: MATLAB code to solve MCFPs
Have been looking for the past couple of days for some robust and
efficient MATLAB routines to solve Minimum Cost Flow Problems (of the
order of 1000+ nodes). Would really appreciate it if somebody could
help me find the same.
-Auro
===
Subject: need help with a basic problem
Need to maximize NPV(net present value) given budget constaint.
Amt available for period 1 is $50 & for period 2 is $20
Cashflows and NPV required each period given below :
Project NPV Cash Outflow Cash Outflow
Period1 Period2
1 $14.00 $12.00 $3.00
2 17 54 7
3 17 6 6
4 15 6 2
5 40 30 35
6 12 6 6
7 14 48 4
8 10 36 3
9 12 18 3
Set this up as a LP problem and solve.
===
Subject: Counter Intuitive Results: Optimising an FFT routine for Speed
ratio is now approaching the noise level asymptote, though there was
at least one good point :
> compressed plain-text copy of the whole thing at
> X-Last-Modified: February 7, 1999
Oh and by the way, since I have been cross-examined on cross-posting :
> alt.comp.lang.learn.c-c++,comp.answers,alt.answers,news.answers
is the general reader of
alt.answers,news.answers
going to be interested in the finer points of C?
> Whoa, whoa, hold on, you want this to be fast on a DSP chip? But you're
> planning to benchmark and optimize it on a general-purpose CPU?
Absolutely right - I am doing high level system engineeting here, I
will pass on the algorithm implementation and source code (probably
FFTW derived, but may still be Numerical Recipes version enhanced on
several counts) to the S/W engineer, he implements a version in Visual
C++ for the client (PC) workstation, and the embedded version gets
ported to whatever DSP is chosen for the new hardware platform. The
DSP speed is the most critical, but speed on the PC client is also
important. The lessons learned probably remain with me.
> The OP should find a group or mailing list specific to the processor
.....
The reality of the situation is that this won't happen.
> The x86 registers are extended precision (or at least double, depending
> upon the mode), so your compiler doesn't have much choice in the matter.
> Nor does this affect the speed.
Surely the implication of this is that float variable arithmetic
should be slower than double arithmetic, because of the extra casts?
> Single precision is almost always faster than double precision, though,
> at least for significant transform sizes, simply because of cache usage.
should be the same speed as for float (there should be cast overheads
though to slow the float version down).
> You asked for something with more resolution than clock(), so I told you
> the best resolution you can get; if you're happy with clock(), why did
> you ask?
> You also inexplicably grabbed the code from cycle.h labeled as
gcc-specific.
toolbox, when I have something apart from Turbo-C I am sure I will use
it.
2. What resolution is this, are we talking of about 1ns, 10ns or 100ns
(and I do remember that you or someone else said it is for relative
comparisons anyway, but there will be a precise answer for a fixed CPU
and clock speed)?
3. I am not happy with clock, I am distinctly unhappy with same. Not
only is 18.2ms an age for most fast applications, it only actually
can get a resolution of
666ns (900MHz Athlon, Win-98)
and worst case
1145[Micro]s (998MHz Celeron, Win-2k).
4. GCC specific? - I grabbed the code for
/*
* Pentium cycle counter
*/
#if (defined(__GNUC__) || defined(__ICC)) && defined(__i386__) ...
__asm__ __volatile__(rdtsc: =A (ret));
I have also tried
/*
* X86-64 cycle counter
*/
#if defined(__GNUC__) && defined(__x86_64__) &&
!defined(HAVE_TICK_COUNTER)
asm volatile(rdtsc : =a (a), =d (d));
This at least generates the compiler error that in-line assembly is
not allowed, the __ASM__ keyword not being recognised at all in the
Pentium example.
Any chance of being able to get the assembled hex codes for these
assembler statements and patching these into the executable. I
imagine inserting a unique sequence of straight C statements that a
search and replace executable could replace with the opcodes for
rdtsc: =A (ret)
Seems bizarre I know, but pragmatism was never thus handicapped.
____________________________________________________________________
Another great Dane has made free
With a question of Be or Not be.
Now might Schr.9adinger's puss,
In descending by Schuss,
Leave one track on each side of a tree?
____________________________________________________________________
Quantum mechanics, hmmm. You put a cat in a box, along with a hammer
and
some poison and a radioactive isotope ... I forget exactly how this
goes.
Anyway, keep some bandages on hand, because I guarantee the cat won't
be
happy.
____________________________________________________________________
===
Subject: Re: Counter Intuitive Results: Optimising an FFT routine for Speed
> ratio is now approaching the noise level asymptote, though there was
> at least one good point :
> compressed plain-text copy of the whole thing at
> Well its plain text, and it is contiguous, which is a lot better than
> my 331 html files totalling 982kB, but up to date it is not :
> X-Last-Modified: February 7, 1999
I said it seems to be the most up-to-date-version available; as far as
I can tell, it is.
===
Keith Thompson (The_Other_Keith) kst@cts.com
San Diego Supercomputer Center <*>
Schroedinger does Shakespeare: To be *and* not to be
===
Subject: Re: Counter Intuitive Results: Optimising an FFT routine for Speed
>Whoa, whoa, hold on, you want this to be fast on a DSP chip? But you're
>planning to benchmark and optimize it on a general-purpose CPU?
> Absolutely right - I am doing high level system engineeting here
The point is that optimizing for a DSP chip (with all sorts of
special-purpose instructions etc.) is completely different from
optimizing for a general-purpose CPU, and the same code, or even the
same *style* of code, is unlikely to be fast on both kinds of architecture.
> Surely the implication of this is that float variable arithmetic
> should be slower than double arithmetic, because of the extra casts?
There isn't a one-to-one mapping from C operators to hardware
operations. The cast from float to double (or actually, to 80-bit
extended precision) is done costlessly in hardware as part of the
load/store.
> 4. GCC specific? - I grabbed the code for
> #if (defined(__GNUC__) || defined(__ICC)) && defined(__i386__) ...
Nota bene: __GNUC__. (If you scroll down in the file, you'll see
alternate code for Visual C++.) Anyway, I think that this is the least
of your worries.
Good luck to you. I hereby declare this thread dead.
===
Subject: Make movie out of sequential print screens?
I have a complicated display that I can save as an image using
ctrl+printscrn. The problem is that the display updates over a long
period of time and I would like to save each frame without having to
personally perform ctrl+printscrn. After each frame is saved, I would
like to stitch them together to make a AVI or MPEG file.
Is there a way do this?
===
Subject: Re: Make movie out of sequential print screens?
This question is way off topic.
The easiest format to get to might be Motion JPEG.
>I have a complicated display that I can save as an image using
>ctrl+printscrn. The problem is that the display updates over a long
>period of time and I would like to save each frame without having to
>personally perform ctrl+printscrn. After each frame is saved, I would
>like to stitch them together to make a AVI or MPEG file.
>Is there a way do this?
===
Subject: Numerical optimization problem: smallest hypersphere?
Hi all,
given a euclidian vector space V (dimension < 10) and a (possibly rather
large) set of vectors x_i within V, I am looking for the center y and
radius r of the smallest hypersphere containing these points.
Therefore, I need to calculate
r = min_y max_i | x_i - y |
efficiently, as I need to do many of these optimizations.
As far as I can tell, the objective function f(y) = max_i | x_i - y | is
convex and Lipschitz continuous (if I recall the definition correctly),
thus it should possess a unique minimum. However, because of the max
it is not differentiable at the optimum. So my question is: which
algorithm for unconstrained optimization should I consider?
Playing around with a Nelder-Mead type optimizer gave somewhat
discouraging results. The tested algorithm proved to be highly
sensitive to the initial guess and converged to non-optimal points
where it appeared to stick. This is definitely not good. Since it is
also sometimes said that there is no firm mathematical basis for this
method, I am looking for something else.
Can anyone point to algorithms better suited to the above problem?
Freely available descriptions/write-ups would be great, as well as a
reference implementation (any higher-level language).
===
Subject: Re: Numerical optimization problem: smallest hypersphere?
>Hi all,
given a euclidian vector space V (dimension < 10) and a (possibly rather
>large) set of vectors x_i within V, I am looking for the center y and
>radius r of the smallest hypersphere containing these points.
Therefore, I need to calculate
r = min_y max_i | x_i - y |
efficiently, as I need to do many of these optimizations.
As far as I can tell, the objective function f(y) = max_i | x_i - y | is
>convex and Lipschitz continuous (if I recall the definition correctly),
>thus it should possess a unique minimum. However, because of the max
>it is not differentiable at the optimum. So my question is: which
>algorithm for unconstrained optimization should I consider?
Playing around with a Nelder-Mead type optimizer gave somewhat
>discouraging results. The tested algorithm proved to be highly
>sensitive to the initial guess and converged to non-optimal points
>where it appeared to stick. This is definitely not good. Since it
is
>also sometimes said that there is no firm mathematical basis for this
>method, I am looking for something else.
Can anyone point to algorithms better suited to the above problem?
>Freely available descriptions/write-ups would be great, as well as a
>reference implementation (any higher-level language).
>--
-ha
minimizing the max is the same as minimizing the max squared.
name the value of the max M.
then
norm(x_i - y )^2 <= M for all i and M is the smallest such number
hence minimize the function
f(y,M) def = M
subject to the constraints
(x_i-y)'(x_i-y) -M <= 0 for all i
happily enough this is a convex optimization problem with a linear
objective function and quadratic constraints , the number of unkowns
being quite modest (dim(y)+M ).
since you can provide also a reasonable initial and feasible guess
y0=center of mass M=the maximum radius
any reasonable nonlinear optimization code should do.
for problems of this type second order cone programming codes should
be good.
see
http://plato.la.asu.edu/topics/problems/nlores.html
for codes
if for some reason you do not want to use one of these packages then
here even a simple augmented Lagrangian code should do quite well
since the problem has such a nice structure
===
Subject: SVD of a large matrix or different solution.
X-ID: SAj69yZ1oeg5jGI1U3aYlua3yoMv+tEQP-8U+0I8LnTwpMJZz1m5wT
I'm trying to fit a 4-6 dimensional spline to scattered
datapoints.
So, essentially, I have to solve Ax=b.
The corresponding matrix A gets quite big, say, 10kx70k , with about
1k*70k of non zero entries.
Since the matrix is most probably rank deficient, I think a SVD is the
right way to go. Am I correct?
The matrix is quite large, so a out-of-core Algorithm seems in order.
Does anyone have any pointers to a library capable of decomposing such
a matrix?
Or some references/suggestions?
----
Jan C. Bernauer
===
Subject: Re: SVD of a large matrix or different solution.
I'm trying to fit a 4-6 dimensional spline to scattered
datapoints.
>So, essentially, I have to solve Ax=b.
>The corresponding matrix A gets quite big, say, 10kx70k , with about
>1k*70k of non zero entries.
Since the matrix is most probably rank deficient, I think a SVD is the
>right way to go. Am I correct?
The matrix is quite large, so a out-of-core Algorithm seems in order.
Does anyone have any pointers to a library capable of decomposing such
>a matrix?
Or some references/suggestions?
>----
> Jan C. Bernauer
for sure you have a 70kx10k problem , i.e. a linear least squares problem
with 10k unknowns and a sparse matrix. or, so to say
f(x)=norm(Ax-b,2)^2
a quadratic function of 10k variables to be minimized, with cheap
operations to compute the gradient and Hessian times a vector.
no problem.
see
http://plato.la.asu.edu/topics/problems/nlolsq.html
for ready to use codes. you may also have a look at the
approximation section of our guide which has pointers to scattered data
fits.
===
Subject: Re: SVD of a large matrix or different solution.
X-ID: Ek1jf-ZTQe8ktlQhUjyAHK6qNOtyuUj9wZ3VGGCDwS-gIA++OHPQYO
nospamspellucci@fb04373.mathematik.tu-darmstadt.de ( Spellucci)
>for sure you have a 70kx10k problem , i.e. a linear least squares problem
>with 10k unknowns and a sparse matrix. or, so to say
> f(x)=norm(Ax-b,2)^2
Exactly!
> a quadratic function of 10k variables to be minimized, with cheap
> operations to compute the gradient and Hessian times a vector.
> no problem.
> see
> http://plato.la.asu.edu/topics/problems/nlolsq.html
> for ready to use codes. you may also have a look at the
> approximation section of our guide which has pointers to scattered data
> fits.
----
Jan C. Bernauer
===
Subject: Re: SVD of a large matrix or different solution.
> I'm trying to fit a 4-6 dimensional spline to scattered
datapoints.
> So, essentially, I have to solve Ax=b.
> The corresponding matrix A gets quite big, say, 10kx70k , with about
> 1k*70k of non zero entries.
> Since the matrix is most probably rank deficient, I think a SVD is the
> right way to go. Am I correct?
Probably not. Since you want to fit 70000 coefficients to
10000 data points, you may well get a completely spurious solution.
(Try fitting 7 coefficients to 1 data point!)
Therefore you need to check whatever approach you use by fitting
the model to 95% of your data and determine the quality of the fit
using the remaining 5%. (This is called cross validation.)
A better way to proceed is to add equations kJx=0 where J is an
approximation to second derivative operators and k a suitable
scaling constant. Then solve (A^TA+k^2 J^TJ)x=A^Tb by conjugate gradients
(CG).
For background see, e.g.
http://www.mat.univie.ac.at/~neum/papers.html#regtutorial
[~ is a tilde]
For CG see any book on iterative methods.
> The matrix is quite large, so a out-of-core Algorithm seems in order.
With CG, you can generate the multiplication maps x to Ax, A^Tx, Jx, J^Tx on
the
fly and need not much storage.
> Or some references/suggestions?
http://www.mat.univie.ac.at/~neum/stat.html#fit
has pointers to scattered data fit routines.
===
Subject: Re: SVD of a large matrix or different solution.
I'm trying to fit a 4-6 dimensional spline to scattered
datapoints.
> So, essentially, I have to solve Ax=b.
> The corresponding matrix A gets quite big, say, 10kx70k , with about
> 1k*70k of non zero entries.
Since the matrix is most probably rank deficient, I think a SVD is the
> right way to go. Am I correct?
>Probably not. Since you want to fit 70000 coefficients to
>10000 data points, you may well get a completely spurious solution.
>(Try fitting 7 coefficients to 1 data point!)
Sorry, got that one in the wrong order. 70000 data points to 10000
coefficients. Somehow I'm used to the screen dimension
ordering ;)
>Therefore you need to check whatever approach you use by fitting
>the model to 95% of your data and determine the quality of the fit
>using the remaining 5%. (This is called cross validation.)
>A better way to proceed is to add equations kJx=0 where J is an
>approximation to second derivative operators and k a suitable
>scaling constant. Then solve (A^TA+k^2 J^TJ)x=A^Tb by conjugate gradients
(CG).
>For background see, e.g.
>http://www.mat.univie.ac.at/~neum/papers.html#regtutorial
>[~ is a tilde]
>For CG see any book on iterative methods.
I will look into that, thanks.
> The matrix is quite large, so a out-of-core Algorithm seems in order.
>With CG, you can generate the multiplication maps x to Ax, A^Tx, Jx, J^Tx
on the
>fly and need not much storage.
That would make things a lot easier.
> Or some references/suggestions?
>http://www.mat.univie.ac.at/~neum/stat.html#fit
>has pointers to scattered data fit routines.
I'll check that, too!
----
Jan C. Bernauer
===
Subject: the definition of frequency in fftw
Hi all,
I am using fftw as my fourier transform tools. I am going to take
numerical derivative with FFT. i.e.
D[u,t] = ifft(i*w*fft(u))
where w is the frequency and i is for the sqrt(-1).
According to numerical recipe in C/C++, I got the angular frequency
is of the form
w = 2*pi*n/(N * dt ), n = -N/2, ... , N/2
where N is the length of the signal and dt is the sampling interval.
However, from the document of fftw, the k-th output corresponds to the
frequency k/T
here T is the total sampling period.
According to the definition from fftw, I take
w = n/N
It is quite different from the definition in the previous one (dT and
2*pi has gone !!!). A testing program show that the result is correct
only when I take
D[u, t] = ifftw( i*w*fftw(u) ) and w=n/N;
===
Subject: Re: the definition of frequency in fftw
>According to numerical recipe in C/C++, I got the angular frequency
>is of the form
> w = 2*pi*n/(N * dt ), n = -N/2, ... , N/2
It says nothing of the kind. See the paragraph after Eq. (12.1.8).
> However, from the document of fftw, the k-th output corresponds to the
> frequency k/T
> here T is the total sampling period.
This is the *frequency*, not the angular frequency (which is what you
need for differentiation); for the latter you multiply by 2*pi.
> According to the definition from fftw, I take
> w = n/N
Besides the 2*pi factor, you are forgetting about aliasing. You have a
choice of whether to call something a frequency of n/N or (n-N)/N or ...
The best approach for numerical differentiation (e.g. for solving
differential equations) is usually to take the angular frequency with
the smallest possible magnitude, i.e. for FFTW:
w = 2*pi*n/(N*dT), n = 0, ..., [N/2]-1, -[N/2], ..., -1
where [..] denotes the greatest integer (floor) function.
(This is the same as what the Numerical Recipes authors give, but to
actually use it for differentiation in their case you would multiply by
-i*w instead of i*w, because they use the opposite sign convention from
FFTW.)
All of this follows simply by inspection of the DFT (or rather the IDFT)
definition, by the way.
===
Subject: Re: the definition of frequency in fftw
> However, from the document of fftw, the k-th output corresponds to the
> frequency k/T
> here T is the total sampling period.
> This is the *frequency*, not the angular frequency (which is what you
> need for differentiation); for the latter you multiply by 2*pi.
Yes. In matlab (it also apply fftw, right?), I take the an angular
frequency as the differentitator.
> According to the definition from fftw, I take
> w = n/N
> Besides the 2*pi factor, you are forgetting about aliasing. You have a
> choice of whether to call something a frequency of n/N or (n-N)/N or ...
> The best approach for numerical differentiation (e.g. for solving
> differential equations) is usually to take the angular frequency with
> the smallest possible magnitude, i.e. for FFTW:
> w = 2*pi*n/(N*dT), n = 0, ..., [N/2]-1, -[N/2], ..., -1
> where [..] denotes the greatest integer (floor) function.
Actually, I do define w like your way. I take N = 2^k, k is an
positive integer.
and n = 0, 1, ... N/2, -N/2+1, ..., -1
w = 2*pi*n/(N*dT), (*)
> (This is the same as what the Numerical Recipes authors give, but to
> actually use it for differentiation in their case you would multiply by
> -i*w instead of i*w, because they use the opposite sign convention from
> FFTW.)
It is not the same. In fftw, w = n/N. In which , 2*pi and dT has gone
!!! If I take -i*w as the differentitator and define w as (*), the
result is quite different from what we expect unless we forget 2*pi
and dT. ??????
===
Subject: Re: the definition of frequency in fftw
> Yes. In matlab (it also apply fftw, right?), I take the an angular
> frequency as the differentitator.
This follows trivially from the differentiation rule for exponentials;
it's not specific to matlab or fftw or numerical recipes. The fact that
you can't derive this for yourself indicates something important.
> It is not the same. In fftw, w = n/N. In which , 2*pi and dT has gone
> !!!
Where in FFTW does it say that? Nowhere. The FFTW manual says that the
*frequency*, not the angular frequency (= 2*pi * frequency), is n/N. As
for the factor of dT, that is just a question of units, and the FFTW
manual explicitly gives the alternative formula n/T, where T is the
total sampling time (i.e. T = N*dT). The FFTW manual also explicitly
mentions the aliasing of positive and negative frequencies.
You need to make some effort (on your own) to have some understanding of
what you are doing; I get the sense that you just want to blindly copy
formulas.
===
Subject: Re: Solve Benedict Webb Rubin Equation of state
> Can any body assist me in solving the BWR equation for all its roots
> analytically.
>
What is the equation? If it is a transcendental equation or a polynomial
of degree larger than 4, there is no analytic (i.e. closed-form) solution.
You then have to do it numerically. Getting ALL the roots may be hard.
===
^^^^^^^^^^^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/
===
Subject: Searching for introductory statistical modeling text
I am looking for an introductory (graduate level physicist) text that
describes modern statistical modeling. In particular, I want to
understand blurbs such as:
The first technique uses bootstrapping and is thought to be as accurate
using less data than the second method which employs only brute force.
For the first method, a data set of size uniformly distributed between 25
and 50 was generated and analyzed. The second method has a data set size
between 100 and 200, which was subsequently generated and analyze it.
This process was repeated 1000 times.
For variance reduction, we want the random numbers used in the two methods
to be the same for each of the 1000 comparisons...
Statistics and I have never been the best of friends... words like
bootstrapping, variance reduction, and interference are all lost
on
me at the moment. I know what a mean, covariances, kurtosis, skew,
pdf, etc is...
Jamie
===
Subject: correlation metric
I'm looking for some good metrics which may be derived through
a correlation matrix. I.e. I'm only interested in one single
output number as a function of the correlation matrix.
If the correlation matrix is of size MxM and there is full
correlation between all elements then I would like to get say, 1.
With correlation between only two elements M-1 is a good number while
under no correlation M would be a nice number.
Any tips ? References would be nice. Thank you.
===
Subject: Re: correlation metric
> I'm looking for some good metrics which may be derived through
> a correlation matrix. I.e. I'm only interested in one single
> output number as a function of the correlation matrix.
> If the correlation matrix is of size MxM and there is full
> correlation between all elements then I would like to get say, 1.
> With correlation between only two elements M-1 is a good number while
> under no correlation M would be a nice number.
> Any tips ? References would be nice. Thank you.
Try (trace[R])^2 / trace[R^2], where R is the correlation matrix.
===
Subject: Re: correlation metric
> Try (trace[R])^2 / trace[R^2], where R is the correlation matrix.
Is there a name for this metric ? Is it used in some applications?
===
Subject: Re: correlation metric
> Is there a name for this metric ?
Not that I know of. In this case I would call the result the
effective dimensionality of your set of M variables.
> ... Is it used in some applications?
(trace[C])^2/trace[C^2], where C is a double-centered covariance
matrix, is Box's adjusted degrees of freedom (to correct for
non-sphericity) in repeated-measures analysis of variance.
For a frequency distribution {f1,f2,...} over unordered categories,
(Sum[f])^2/Sum[f^2] gives the effective number of categories over
which the observations are spread.
For data sets composed of i.i.d. observations in which each case has
associated with it an arbitrary weight w that is used in calculating
summary statistics (e.g., mean[x] = Sum[w*x]/Sum[w]), the effective
sample size is (Sum[w])^2/Sum[w^2].
There are almost certainly others, but none come to mind just now.
===
Subject: Re: Percentile Estimation in a Monte Carlo Simulation
> I am trying to estimate a fairly extreme percentile (99.7) for a
> distribution obtained from Monte Carlo simulation. I am facing a
> problem because of the instability of data at the extremes of the
> distribution, as in, if the percentile values are calculated for
> different simulation runs, they seem to jump around a lot. I am
> currently using the MATLAB prctile function to estimate the percentile
> and am limited to estimate it with 1000 iterations.
> I would really appreciate if anybody can suggest me a better algorithm
> than what is being used in the prctile function, which can provide a
> closer estimate to the true percentile.
> Uttam
David Jones has offered sound advice. If you have Excel (it has a
percentile function ... subject to comments users here may like to
make!), you may wish to try it as an alternative as you will not be
limited to 1000 iterations (rather strange that Matlab restricts you
to this rather useless limitation).
===
Subject: Linux + MS Matlab + Wine = ?
X-ZLExpiry: -1
X-ZLReceiptConfirm: N
Mail-To-News-Contact: postmaster@nym.alias.net
Is it possible to run the windows version of matlab under Linux
with WINdows Emulator? I know there is a Linux Matlab, but I only
have the Windows version.
BTW doesn't WINE implement the complete Windows API? How
come some software does run on it? Is it because of bugs
in such software or because it relies on undocumented
features?
===
Subject: Re: Linux + MS Matlab + Wine = ?
> Is it possible to run the windows version of matlab under Linux
> with WINdows Emulator? I know there is a Linux Matlab, but I only
> have the Windows version.
> BTW doesn't WINE implement the complete Windows API? How
> come some software does run on it? Is it because of bugs
> in such software or because it relies on undocumented
> features?
Have you considered using GNU Octave, a Matlab replacement written by
Eaton and others?
http://www.octave.org/
It is free software, i.e. available under the GNU General Public License.
Written for Linux and ported to run under Windows using cygwin, it
might
serve your purpose nicely.
regards, chip
===
Subject: Re: Linux + MS Matlab + Wine = ?
X-Enigmail-Version: 0.76.1.0
X-Enigmail-Supports: pgp-inline, pgp-mime
> Is it possible to run the windows version of matlab under Linux
> with WINdows Emulator? I know there is a Linux Matlab, but I only
> have the Windows version.
> BTW doesn't WINE implement the complete Windows API? How
> come some software does run on it? Is it because of bugs
> in such software or because it relies on undocumented
> features?
I have tried to install it under various Wine versions in a Mandrake
installation and it is the only program which has been able to freeze
the computer completely - so no I don't think it is worth trying :-)
But with windows 98 running under win4lin there is no problem using the
windows version of matlab.
.
===
Subject: Re: Linux + MS Matlab + Wine = ?
> Is it possible to run the windows version of matlab under Linux
> with WINdows Emulator? I know there is a Linux Matlab, but I only
> have the Windows version.
If you encounter problems using Wine, try VMware.
Evaulation versions should be available from their web site.
It's pretty sure that you wouldn't have trouble with Matlab using
that emulator. But it doesn't come for free...
Bernhard
===
Subject: Re: Linux + MS Matlab + Wine = ?
: Is it possible to run the windows version of matlab under Linux
: with WINdows Emulator? I know there is a Linux Matlab, but I only
: have the Windows version.
: BTW doesn't WINE implement the complete Windows API? How
: come some software does run on it? Is it because of bugs
: in such software or because it relies on undocumented
: features?
Not every software uses everything of the Windows API. In fact, Software
not
written by Microsoft normally only uses only a small part of the Api...
Bye
===
Uwe Bonnes bon@elektron.ikp.physik.tu-darmstadt.de
Institut fuer Kernphysik Schlossgartenstrasse 9 64289 Darmstadt
--------- Tel. 06151 162516 -------- Fax. 06151 164321 ----------
===
Subject: Re: Linux + MS Matlab + Wine = ?
> BTW doesn't WINE implement the complete Windows API? How
> come some software does run on it? Is it because of bugs
> in such software or because it relies on undocumented
> features?
Wine is alpha software, it's under development, so it simply isn't
finished.
===
Subject: 2 axis matrix question!
first off im only in grade 9 so i wont know much, if u can provide me
with a straight forward formula that would be excellent.
Ok here's my drill,
I have an object on a 2d matrix and I only know the position and
facing angle of that object. I want to know which point is exactly or
approximately 5 units away from the object going along the angle the
object is facing. the object can be anywhere on the matrix and can
face any angle, but still it's destination must be 5 units along its
facing angle.
Please someone help me. I need this for reprogramming a game.
===
Subject: Supreme
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision:
1.9 primary) id h99Hjt410223;
I'm in my second week of my first year of university and hoping for a
little help with my analysis homework as I am finding it quite a
struggle. The question is:
Let A and B be bounded, nonempty sets of real numbers.
Is it always true that sup(AB) = supAsupB, where AB = {xy: x belongs
to A, y belongs to B}?
I know that it is not true that sup(AB) = supAsupB but I am having
difficulty in finding a counterexample. Not quite sure where to start
looking.
I'd appreciate any help or hints on this.
Also if anyone can recommend some good books for beginners trying to
learn Analysis that would be great.
Katherine
===
Subject: Re: Supreme
> I'm in my second week of my first year of university and hoping for a
> little help with my analysis homework as I am finding it quite a
> struggle. The question is:
> Let A and B be bounded, nonempty sets of real numbers.
> Is it always true that sup(AB) = supAsupB, where AB = {xy: x belongs
> to A, y belongs to B}?
> I know that it is not true that sup(AB) = supAsupB but I am having
> difficulty in finding a counterexample. Not quite sure where to start
> looking.
> I'd appreciate any help or hints on this.
I'll give one hint. sup() is all about the >= relation.
Is it always true that b>=c ==> ab>=ac?
> Also if anyone can recommend some good books for beginners trying to
> learn Analysis that would be great.
Clapham's Introduction to real analysis is very short
and pretty clearly written, but it seems to be out of
print.
A book I remember fondly is Ralph Boas's A primer of real
functions. I read this as a slightly precocious school
pupil a year or two before going to university, and loved
it, which at least suggests that it's clearly written and
accessible to beginners. There's been a new edition since
then, and in any case I haven't a copy to hand, so I can't
make any guarantees :-).
===
Gareth McCaughan
.sig under construc
===
Subject: Help on 3D peak fit
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision:
1.9 primary) id h9A1krV13760;
I have ploted data from a matrix as a 3D graph or a contour map (using
Origin). These plots are showing multiple peaks and I would like to
fit these peaks with 3D Gaussians to get numerical parameters of each
peak.
Does anyone know a way and a software to fit these peak and get the
parameters?
Thank You
===
Subject: Re: Help on 3D peak fit
>I have ploted data from a matrix as a 3D graph or a contour map (using
>Origin). These plots are showing multiple peaks and I would like to
>fit these peaks with 3D Gaussians to get numerical parameters of each
>peak.
>Does anyone know a way and a software to fit these peak and get the
>parameters?
>Thank You
>
the way to go depends....
if your model is
f(x,y) = sum a_i exp(-((x-x_i)^2+(y_y_i)^2)/s_i^2)
then you have a sum of radial functions and can use codes for fitting by
radial
basis functions. there are some in
http://www.netlib.org/toms
if your model is more general, namely
f(x,y)= sum a_i exp(-[x-x_i,y-y_i]A_i[x_x_i;y-y_i]) (in matlab
notation)
with positive definite 2 by 2 matrices A_i, then you have a nonlinear least
squares problem in a quite general form.. you can avoid the side condition
on A_i
by replacing it beforehand by
A_i =L_i L_i'
with triangular matrices L_i as the unknowns, but nevertheless there
remains
a condition that the diagonal elements of the L_i should by >= eps>0 for
some
eps. solving this as a least squares problem is in principle routine
and for codes see
http://plato.la.asu.edu/topics/problems/nlolsq.html
but I should warn you: with more than 3 peaks say (i.i. i=1,2,3)
you will get a very hard problem and even good codes may fail or use a
tremendous effort to get some reasonable result. you also must be
able to provide reasonable initial guesses for all parameters.
the a_i appear linear hence you have a so called separable problem and
could in principle express the a_i as functions of the other parameters.
but this makes the problem even stronger nonlinear hence i discourage this
===
Subject: Re: Help on 3D peak fit
>I have plotted data from a matrix as a 3D graph or a contour map (using
>Origin). These plots are showing multiple peaks and I would like to
>fit these peaks with 3D Gaussians to get numerical parameters of each
>peak.
>Does anyone know a way and a software to fit these peak and get the
>parameters?
> but I should warn you: with more than 3 peaks say (i.i. i=1,2,3)
> you will get a very hard problem and even good codes may fail or use a
> tremendous effort to get some reasonable result. you also must be
> able to provide reasonable initial guesses for all parameters.
This is done by discarding temporarily all data note 'belonging' to each
peak
(in turn), and fitting a Gaussian to them, Then solving
a linear least squares problem to weight the resulting Gaussians.
(Requires that you have enough data.)
===
Subject: Re: Help on 3D peak fit
>I have plotted data from a matrix as a 3D graph or a contour map (using
>Origin). These plots are showing multiple peaks and I would like to
>fit these peaks with 3D Gaussians to get numerical parameters of each
>peak.
You should look for clustering or mixture modeling. Something like the
(*)EM
algorithm could solve your problem (if you have enough data).
> This is done by discarding temporarily all data note 'belonging' to each
peak
> (in turn), and fitting a Gaussian to them, Then solving
> a linear least squares problem to weight the resulting Gaussians.
> (Requires that you have enough data.)
Yes. It is faster (in my experience) to start by assigning uniform
probabilities that each data belong to the peak i, and by iterating
something
called the SEM algorithm. It will give you a good start for more classical
algorithm (like EM).
===
===
Subject: Newton-Raphson for complex argument
Given the classic Newton-Raphson solution:
w_(n+1) = w_n - f(w_n)/f'(w-n)
where the process is repeated until the difference between
w_(n+1) and w_n is sufficiently small the operation is
obvious enough when w is a real number.
But... what is the test when w is complex (u+iv)?
Is there any in print reference material for numerical
applications for complex numbers?
===
Subject: Re: Newton-Raphson for complex argument
> Given the classic Newton-Raphson solution:
> w_(n+1) = w_n - f(w_n)/f'(w-n)
> where the process is repeated until the difference between
> w_(n+1) and w_n is sufficiently small the operation is
> obvious enough when w is a real number.
> But... what is the test when w is complex (u+iv)?
> Is there any in print reference material for numerical
> applications for complex numbers?
The test is still the same you just use complex arithmetic and
continue until the differences between w_(n+1) and w_n; lets call it
u+iv; norm (u^2+v^2) or abs sqrt( norm )is sufficient small.
I don't have any print reference material. But one application is
finding zeros of polynomial with either real or complex coefficients.
If you are into that you can download a windows baszed polynomial
solver using the most common method like Newton, Laguerre, Halley and
Graeffe to solve these problems from
www.hvks.com/numercial.htm
===
Subject: Re: Newton-Raphson for complex argument
> Given the classic Newton-Raphson solution:
> w_(n+1) = w_n - f(w_n)/f'(w-n)
> where the process is repeated until the difference between
> w_(n+1) and w_n is sufficiently small the operation is
> obvious enough when w is a real number.
> But... what is the test when w is complex (u+iv)?
|w_(n+1)-w_n| is sufficiently small
> Is there any in print reference material for numerical
> applications for complex numbers?
You may wish to look at Section 5.6 Complex zeros
of my book
Introduction to Numerical Analysis
Cambridge Univ. Press, Cambridge 2001.
For a review of the book, see
===
Subject: Re: Newton-Raphson for complex argument
> Given the classic Newton-Raphson solution:
> w_(n+1) = w_n - f(w_n)/f'(w-n)
^
w_(n+1) = w_n - f(w_n)/f'(w_n)
> where the process is repeated until the difference between
> w_(n+1) and w_n is sufficiently small the operation is
> obvious enough when w is a real number.
> But... what is the test when w is complex (u+iv)?
> Is there any in print reference material for numerical
> applications for complex numbers?
You do it the same way. However, you can run into the same instabilities
as for the real case, esp. near a saddle point where f(z) -> 0.
The basis of the Laguerre method for finding the roots of polynomials
is complex Newton-Raphson. Have a look at Acton, Numerical Methods
that Work (American Methematical Society).
===
^^^^^^^^^^^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/
===
Subject: Looking for a software for an estimation with confidence intervals
A function a(n) is presumed to be asymptotically equivalent to
(n!)^s*A^n*n^beta*M when n->+infinity
I have data (ni,a(ni)) i=1..N and I would like to identify s, A, beta and
M.
I need error estimates on these values (estimation of these values with
confidence intervals)
I presume that several softwares can be used for that, but I don't know
which one to choose (Matlab ? SAS ? SPSS ? others ?). I have no experience
with statistical softwares.
I've used Maple LeastSquares function, applied on ln(a(n)), it works fine
for having values of s, A, beta and M; but it does not give me confidence
intervals.
All suggestions are welcome.
===
Subject: Re: Looking for a software for an estimation with confidence
intervals
I've forgotten a precision: the values of a(n) in the data are very large,
and known with hundreds of digits. Probably high precision will be needed
to
perform the estimation.
===
Subject: Questions on Function Approximation (using FPGAs)
,
I am working on (compound) function approximation with one input
variable using piecewise polynomial approximation with non-linear
joints. These approximations are implemented in hardware using Xilinx
FPGAs.
Example of such functions include: f(x)=sqrt(-ln(x)) or
f(x)=x*ln(x) where x = [0,1), which are used for Gaussian noise
generation (Box-Muller method) and Entropy calculation
respectively.
Does anyone know any other real-life applications where compound
functions need to be approximated?
My second question is on the function f(x)=sqrt(-ln(x)) over x =
[0,1). This function is highly non-linear and approaches infinity
as x gets close to zero. This requires floating point
implementation (due to the large polynomial coefficients, which I
want to avoid). Are there any transformations I am apply to the
function to decompose it 2 or more functions that are more linear?
(Note that ln(x) is also highly non-linear over x = [0,1))
Dong-U Lee
===
Subject: Re: Questions on Function Approximation (using FPGAs)
> ,
> I am working on (compound) function approximation with one input
> variable using piecewise polynomial approximation with non-linear
> joints. These approximations are implemented in hardware using Xilinx
FPGAs.
> Example of such functions include: f(x)=sqrt(-ln(x)) or
> f(x)=x*ln(x) where x = [0,1), which are used for Gaussian noise
> generation (Box-Muller method) and Entropy calculation
> respectively.
> Does anyone know any other real-life applications where compound
> functions need to be approximated?
sqrt(1/x) is used so often that more than one modern cpu contains
built-in lookup tables to generate a good starting point for a NR
iteration.
> My second question is on the function f(x)=sqrt(-ln(x)) over x =
> [0,1). This function is highly non-linear and approaches infinity
> as x gets close to zero. This requires floating point
> implementation (due to the large polynomial coefficients, which I
> want to avoid). Are there any transformations I am apply to the
> function to decompose it 2 or more functions that are more linear?
> (Note that ln(x) is also highly non-linear over x = [0,1))
My first guess would be to look for some kind of rational approximation,
even if this does require a final division.
Terje
===
-
almost all programming can be viewed as an exercise in caching
===
Subject: Questions on Function Approximation (using FPGAs)
,
I am working on (compound) function approximation with one input
variable using piecewise polynomial approximation with non-linear
joints. These approximations are implemented in hardware using Xilinx
FPGAs.
Example of such functions include: f(x)=sqrt(-ln(x)) or
f(x)=x*ln(x) where x = [0,1), which are used for Gaussian noise
generation (Box-Muller method) and Entropy calculation
respectively.
Does anyone know any other real-life applications where compound
functions need to be approximated?
My second question is on the function f(x)=sqrt(-ln(x)) over x =
[0,1). This function is highly non-linear and approaches infinity
as x gets close to zero. This requires floating point
implementation (due to the large polynomial coefficients, which I
want to avoid). Are there any transformations I am apply to the
function to decompose it 2 or more functions that are more linear?
(Note that ln(x) is also highly non-linear over x = [0,1))
Dong-U Lee
===
Subject: Re skyline matrix solver
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision:
1.9 primary) id h9ADf1530198;
You can use Spooles, among others:
http://www.netlib.org/linalg/spooles/spooles.2.2.html
===
Subject: Maybe it is a new number factorization algorithm.
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision:
1.9 primary) id h9ADf3S30202;
Hi everyone,
First I want to sorry for my language mistakes, because it is not my
nation language.
It is possible that I had just found a new number factorization
algorithm. It is grounded by equation system solving and nothing more.
You create an equation system, solve it and get results. For example:
Your number: N:= 120;
You create equation system: Eqs:= create_my_equation_system(N);
Solve it: solve(Eqs);
after my algorithm solves the equation system it returns this results
(similar):
1 * 120, 2 * 60, 3 * 40, 4 * 30, 5 * 24, 6 * 20, 8 * 15, 10 * 12.
It works fine with any length of number. But solving time is not
small, jet :-)
Before I publish this my algorithm I want to check out if there is no
any other similar algorithm to my. In this case I want to register it.
Well my question would be to ask if you know any existing similar
algorithm to my (or you think it is similar) please send it, or write
where I can find it.
Mindaugas Rukas
===
Subject: roots of quadratic matrix equation?
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision:
1.9 primary) id h9ABpO622561;
Is there a formula for the roots of the following quadratic matrix
equation?
QAQ - BQ + C = 0
where A, B, C and Q are all symmetric nxn matrices, n>1?
thanks
===
Subject: Re: roots of quadratic matrix equation?
> Is there a formula for the roots of the following quadratic matrix
> equation?
> QAQ - BQ + C = 0
> where A, B, C and Q are all symmetric nxn matrices, n>1?
> thanks
Do you really need a closed formula ? you could solve the equation
above by numerically means: if F(A) = QAQ - BQ + C then
dF/dA = 2*AQ - B. Now you can use Newton iteration:
delta_A = - inv (dF/dA) * F(A)
A = A + delta_A
Uwe.
===
Dr. rer. nat. Uwe Schmitt http://www.procoders.net
schmitt@procoders.net A service to open source is a service to
mankind.
===
Subject: golf shaft flex equation
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision:
1.9 primary) id h9B3OEc22958;
The reference numbers you are referring to are part of the Brunswick
Frequency Coefficient Measurement system for measuring/anotating shaft
stiffness in a golf club.
They are based on the resonant frequency of a 43.5 length of shaft,
clamped at one end (5) with a weight of 200 grams at the tip. The
frequency 'slope' of these shafts is 8.6 cpm per inch. ie as the shaft
gets shorter the frequency increases at that rate. (or vice versa) The
FCM numbers, such as 2.0, 4.5, 6.0, 7.0 etc are merely related to the
measured cycles per minute and represent (cpm-200)/10, eg. FCM 4.5 is
the 'standard' of a 43.5 shaft/club system which resonates at 245
cpm. You can arrange a simple chart with 'frequency' on the vertical
axis, and shaft length on the horizontal axis. Spot the relevant
frequency standards at the 43.5 length (ie the 2.0, 4.5, 6.0 's etc),
and draw straight lines from these points back down to the minimum
lengths, sloping down (from right to left) at the rate of 8.6 cpm per
inch of reducing length. Thus you can then read off the appropriate
frequencies at the different shaft lengths, given the 'stiffness'
nomenclature of the actual club, or if you have a particular club of
known frequency you can read off from the chart the relevant stiffness
'standard' of that club.
===
===
Subject: Bode Plots from Data. - Matlab
by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision:
1.9 primary) id h9BC2QM22328;
I have tons of data of an oscillation of some control laws. I am
having a difficult time finding a way to create a bode plot(frequency
response) of the data. I need to do this in Matlab
If any of you could lend me a hand, you can e-mail me at
fkatzenb@usermail.com or katzenbergfl@navair.navy.mil
===
Subject: Discrete Chebyshev: best algorithm?
I'm interested in discrete chebyshev approximations. I have
to use linear functions, in the meaning that I have to use
y=mx+q
and I am interested in using the best available algorithm.
I know that there is an algorithm from Borrodale and
lips (for the general problem) that is implemented in
the nag library (by the way, how can I get it in C?.. how
algorithms by Megiddo and Seidel that work in
O(n) for linear programming in small dimension. I am
currently using M. Hohmeyer implementation of Seidel
randomized algorithm.
But what is the best one? Is Borrodale' and lips' still
preferable?
===
Subject: Re: Discrete Chebyshev: best algorithm?
> I'm interested in discrete chebyshev approximations. I have
> to use linear functions, in the meaning that I have to use
> y=mx+q
> and I am interested in using the best available algorithm.
> I know that there is an algorithm from Borrodale and
> lips (for the general problem) that is implemented in
> the nag library (by the way, how can I get it in C?.. how
> much money?)
nag_linf_fit is part of the NAG C Library - you can buy it.
(see http://www.nag.co.uk/numeric/CL/manual/pdf/E02/e02gcc.pdf )
It seems to use a modified version of ACM algo 495, for which the
Fortran source is available at
http://www.netlib.org/toms/495
If neither efficiency nor readability is an issue for your application
you could try to throw the Fortran source into f2c and see what
comes out ...
> algorithms by Megiddo and Seidel that work in
> O(n) for linear programming in small dimension. I am
> currently using M. Hohmeyer implementation of Seidel
> randomized algorithm.
> But what is the best one? Is Borrodale' and lips' still
> preferable?
> Thank you very much
> m
===
Subject: 4th order Bspline
I have a stupid question about the 4th order Bspline. In the paper(
see link), they use summation to denote a B-spline. But why they use a
vector (4 j)~T.
Can anyone help me to understand it?
thanks.
http://www.ri.cmu.edu/pub_files/pub2/wu_yu_te_1998_1/wu_yu_te_1998_1.pdf
===
Subject: Re: 4th order Bspline
the 4th order Bspline is denoted in equation 1 in the paper below.thanks.
> I have a stupid question about the 4th order Bspline. In the paper(
> see link), they use summation to denote a B-spline. But why they use a
> vector (4 j)~T.
> Can anyone help me to understand it?
> http://www.ri.cmu.edu/pub_files/pub2/wu_yu_te_1998_1/wu_yu_te_1998_1.pdf
===
Subject: Re: mathematical extrapolation of measurement values
> I need to find an algorithm to extrapolate measurement values gathered
> by a server I am writing. However, I don't know any rules these values
> adhere to, since I might use that algorithm on several sets of
> measurement values with different maths lying beneath.
> Best idea I have come up with so far is to put a curve of nth potence
> through the measurement values of mine and hope that the resulting
> curve allows some extrapolation into the future. So far I have
> realized this for a 2nd degree curve like f(x)=ax*x+b*x+c using three
> points of my measurements. Doesn't work too well, though.
> Requirements and other stuff:
> I know that my measurements won't show any abrupt changes. The values
> I'm watching tend to change rather slow. Furthermore, implementation
> time and calculation time are both restricted. I can't spend much time
> for either - meaning I can spend about 1-2 days for the implementation
> of the algorithm and I'll need to be able to calculate (roughly
> guessed) 1000+ extrapolations per second on a publically available PC
> (making me refrain from iterations). Additionally, the calculation
> should need only as few measurement values as reasonably possible,
> since I don't have a history of my values and need to store my
> calculation-values separately. Extrapolation should start with only a
> hand full of values gathered.
Here's an approach that worked well in retrofitting the (classical BDF)
derivative of the PID controller fcn. Create a sliding interval with
three fcn values, fit it with a cubic spline and then use the endpoint
derivative for extrapolation.
You can adapt spline from netlib/sfmm dir for maximum efficiency
by return(ing) right after the spline derivative (b(n), for n=3) is
calculated.
--
------------------------------------------------------
Applied Algorithms http://sdynamix.com
===
Subject: Runge Kutta ABC's inner working ?
Do you kwow where I can found a smooth introduction
to Runge Kutta methods (internet site).
I know there is google but I want your opinion
about a site for the non mathematician layman.
Thank you for sharing
===
Subject: Orthogonal polynomials , weight w(x)=exp(-k*x) , (a,b)=(0,1).
Content-Length: 941
Originator: rusin@vesuvius
Suppose that (P_n(x))_{n>=0} , (P_n(1)=1 ,n=0,1,...) ,
is the polynomial system orthogonal on [0,1] with
respect to weight function w(x)=e^{-x} .
In other words, for all non-negative integers n,j
(P_n,P_j):=Integral_{x=0 to x=1} P_n(x)P_j(x)w(x) dx = 0 when n=/=j
and (P_n,P_n) =/=0 .
I am interested to find coefficients (A_n),(B_n) and (C_n) from
the three-term recurrence
P_{n+1}(x) = (A_n*x+B_n)*P_n(x) - C_n*P_{n-1}(x) , n=1,2,...,
with P_0(x)=1 , P_1(x)=(e-1)*x+(2-e).
The same question in a more general case when w(x)=e^{-k*x},k in R, k=/=0 .
Integrals of form Integral_{x=0 to x=1}e^{-k*x}f(x) dx
arise in quantum mechanical computations (for instance see R. Mach ,
J.Math.Phys 25 (1984)).
Of interest is also the symmetric weight w(x)=e^{-x(1-x)} .
I need the coefficients in order to find some Gaussian quadratures
for above weights.
If possible, please some references. Thanking in advance, Alex.
=========
===
Subject: Re: Orthogonal polynomials , weight w(x)=exp(-k*x) , (a,b)=(0,1).
Content-Length: 1471
Originator: rusin@vesuvius
>Suppose that (P_n(x))_{n>=0} , (P_n(1)=1 ,n=0,1,...) ,
>is the polynomial system orthogonal on [0,1] with
>respect to weight function w(x)=e^{-x} .
>In other words, for all non-negative integers n,j
> (P_n,P_j):=Integral_{x=0 to x=1} P_n(x)P_j(x)w(x) dx = 0 when n=/=j
>and (P_n,P_n) =/=0 .
I am interested to find coefficients (A_n),(B_n) and (C_n) from
>the three-term recurrence
>P_{n+1}(x) = (A_n*x+B_n)*P_n(x) - C_n*P_{n-1}(x) , n=1,2,...,
> with P_0(x)=1 , P_1(x)=(e-1)*x+(2-e).
>The same question in a more general case when w(x)=e^{-k*x},k in R, k=/=0
.
> Integrals of form Integral_{x=0 to x=1}e^{-k*x}f(x) dx
>arise in quantum mechanical computations (for instance see R. Mach ,
> J.Math.Phys 25 (1984)).
>Of interest is also the symmetric weight w(x)=e^{-x(1-x)} .
>I need the coefficients in order to find some Gaussian quadratures
>for above weights.
>If possible, please some references. Thanking in advance, Alex.
>=========
>
you could compute the coefficients using the scalar product itself, see e.g.
the
formulae in stoer-bulirsch. but it is not a good idea to use them for
computing
weights and nodes of quadrature formulae because of problems with roundoff
walter gautschi did a lot of work on this and his code orthpol in
http://www.netlib.org/toms will serve your needs perfectly.
===
Subject: Simplifying logical proposition
I would like to know if anyone could help me understanding how Maple did
the
following simplification (using bsimp):
Proposition: (not p or q) and ((q and not r) or (r and not q))
Simplified proposition: r and not p and not q
Thank you
===
Subject: Re: Simplifying logical proposition
> I would like to know if anyone could help me understanding how Maple did
the
> following simplification (using bsimp):
If it does, good question since what you post below are not equivalent. So
those trying to prove it can stop.
> Proposition: (not p or q) and ((q and not r) or (r and not q))
> Simplified proposition: r and not p and not q
p q r ~p or q q and ~r r and ~q (q and ~r) (~p or q)
or (r and ~q)
and ((q and ~r) or (r and ~q))
F F F T F F F
F
F F T T F T T
T
F T F T T F T
T
F T T T F F F
F
T F F F F F F
F
T F T F F T T
F
T T F T T F T
T
T T T T F F F
F
> Simplified proposition: r and not p and not q
p q r r and ~p (r and ~p) and ~q (parenthesis not really
needed)
F F F F F
F F T T T
F T F F F
F T T T F
T F F F F
T F T F F
T T F F F
T T T F F
===
Subject: FFT Spectral Analysis - SpectraLAB, SpectraPLUS,
FFT Spectral Analysis - SpectraLAB, SpectraPLUS,
SpectraPRO, SpectraRTA, acoustic calculation, simulation
and sound system design
SpectraLAB v4.32.16c
SpectraLAB v4.32.17
SpectraPLUS v2.32.03c
SpectraPLUS v2.32.04
SpectraPRO v3.32.16c
SpectraPRO332 v3.32.17
SpectraRTA v1.32.12c
SpectraRTA132 v1.32.13
for more info, please send e-mail
===
Subject: calculating log of the ratio of determinants
I have to calculate the log ratio of two symmetric p.d 500x500
determinants S1 and S2, the determinants are too large and will not
fit into a double variable.
Since I cannot calculate the determinants exactly I can think of the
following ways to do it
1]
calculate S2^{-1}S1, find out its eigenvalues and add up the sum of
their logs
2]
calculate sum of the logs of the ev's of S1, find out sum of the ev's
of S2 and subtract
Is there any reason to prefer one over the another?
===
Subject: Re: calculating log of the ratio of determinants
> I have to calculate the log ratio of two symmetric p.d 500x500
> determinants S1 and S2, the determinants are too large and will not
> fit into a double variable.
> Since I cannot calculate the determinants exactly I can think of the
> following ways to do it
> 1]
> calculate S2^{-1}S1, find out its eigenvalues and add up the sum of
> their logs
> 2]
> calculate sum of the logs of the ev's of S1, find out sum of the ev's
> of S2 and subtract
> Is there any reason to prefer one over the another?
The second is usually much cheaper and is often used in practice.
you can also add the binary exponents and multiply the mantissae
(but if your matrix would be much larger you would occasionally
need to rescale the mantissa product by some power of 2.
===
Subject: sign function
===
Subject: Re: sign function
Please do not post the same question _separately_ to different newsgroups.
If you feel that you _must_ post the same question to more than one group,
crossposting is much better.
I already responded to your question in k12.ed.math . Interested readers
may see the response there, once its moderator posts my message.
David Cantrell
===
Subject: How can I get the frequency of a periodical function
I have n datapoints {(t1,x1),(t2,x2),...,(tn,xn)} of a periodical
Signal t->x(t). Now I search the Fourier-Approximation for this
Signal. When I search the approximation with fft, so I got the
magnitudes ak and bk to the harmonical functions sin(k*w0*t) and
cos(k*w0*t) (where w0=2*pi/(tn-t1) and T=tn-t1).
My problem is, that I don't know the frequency of the signal. The
datapoints are dristributed over a range tn-t1=k*T (k>1,real). How can
I find T?
===
Subject: FFTW and polynomail Multiplication
I am trying to use FFTW to square a polynomial of large degree N. I
have tried the following
----------------------------------------------------------
fftw_complex *in, *out;
in = fftw_malloc(N* sizeof (fftw_complex));
out = fftw_malloc(N* sizeof (fftw_complex));
for(int i=0;i < N/2; i++)
/* in is my row-vector of Array coefficients */
in[i][0]= MY ARRAY COEFFICIENTS;
/* all imaginary parts are 0 */
in[i][1]= 0;
}
p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_MEASURE);
fftw_execute(p);
for(int i=0;i < N; i++)
out[i][0]= out[i][0]*out[i][0];
}
p = fftw_plan_dft_1d(N,out,in, FFTW_BACKWARD, FFTW_MEASURE);
fftw_execute(p);
print(in);
----------------------------------------------------------
now when I print... I get all zeroes in the real part of 'in'.
what am I doing wrong? Is there any hope for me?
===
Subject: Re: FFTW and polynomail Multiplication
> I am trying to use FFTW to square a polynomial of large degree N. I
> have tried the following [....................................]
> now when I print... I get all zeroes in the real part of 'in'.
See question 3.16 in the FFTW FAQ.
===
Subject: Numerical Puzzle
I have been working on a puzzle sponsored by my university but are
unable to solve it. I am asking if someone here can solve it. If you
solve it in time you are eligible for a PocketPC or Microsoft
software. The hint is Fortune Teller Origami and the title is
2914849 steps to go? Better save all my powers! The puzzle is here:
http://www.acm.uiuc.edu/puzzlehack/ph2/puzzles2/1/puzzle1.jpg