Mailing List Archive: 49091 messages
  • Home
  • Script library
  • AltME Archive
  • Mailing list
  • Articles Index
  • Site search
 

GPS: minimizing sum of squares algorithm

 [1/14] from: jjmmes:ya:hoo:es at: 8-Nov-2002 19:41


Has anybody written this algorithm in Rebol ? I need this to georeference gps data, i.e. find a, b, c, d, e, f below given three calibration points x = a + b lat + c lon y = d + e lat + f lon It's available as a math library in most languages but maybe some rebol guru has it too ! Regards

 [2/14] from: jan:skibinski:sympatico:ca at: 8-Nov-2002 14:22


Jose, I could certainly help you with this but I would need to know more about the shape of the functions b c e f in order to model them appriopriately. For example, depending on what those functions represent a polynomial approximation might be of no use or it might be quite good actually. A reference to existing georeference algorithms could be a good starting point. Regards, Jan jose wrote:

 [3/14] from: tomc:darkwing:uoregon at: 8-Nov-2002 12:12


Hi Jose I implemented some linear alrerbra functions a few years back. I will see if I can dig them up over the weekend On Fri, 8 Nov 2002, [iso-8859-1] jose wrote:

 [4/14] from: jjmmes::yahoo:es at: 9-Nov-2002 16:23


a,b,c,d,e,f are constants. I basically want to solve a set of 6 equations with six unknowns, using one of the existing linear algebra algorithms. The algorithm implementation I'm looking for is also known as RMS, root mean squares but couldn't find a good algorithm description in the web Regards Jose --- Jan Skibinski <[jan--skibinski--sympatico--ca]> escribió: > Jose,

 [5/14] from: jan:skibinski:sympatico:ca at: 9-Nov-2002 11:01


Jose,
> a,b,c,d,e,f are constants. > I basically want to solve a set of 6 equations with
<<quoted lines omitted: 4>>
> Regards > Jose
Ah, then this is quite simple. I thought you used Rebol syntax with equal signs put there by mistake instead of ( :), while in fact you wrote mathematical equations with implied * between 'b and 'lat, etc. Sorry for that. I'll post a solution later today. Jan

 [6/14] from: jan:skibinski:sympatico:ca at: 10-Nov-2002 3:40


Jose, As promised, here is a solution to your GPS minimization problem, or rather a part of it. The only function I am attaching here is the main GPS function. I have not tested it, so it might be buggy, therefore threat it as the outline only. I developed and tested today various pieces that the gps function is calling, including a solver for a system of linear equations. If you want to use them just ask for it and I will send it to you. I did not plan to develop any linear algebra package at this time and whatever I developed today is not optimized and quite rough around the edges. For this reason I am not publicly posting it to the library. But you can use any other solver, or any other package with the gps function below. All the best, Jan comment { Georeference GPS data modeling: ------------------------------ Find six constants relating [x y] coordinates and [long lat] coordinates, which we abbreviate here as [u v]: x = a1 + a2 u + a3 v y = b1 + b2 u + b3 v To find a1, a2, a3, b1, b2, b3 we need at least 3 sets of data for x, y, u and v - the more the merrier: [x1 x2 x3 ... xn] [y1 y2 y3 ... yn] [u1 u2 u3 ... un] [v1 v2 v3 ... vn] Minimization of the least square error F: F = Sum (xi - a1 - a2 ui - a3 vi)**2 + (yi - b1 - b2 ui - b3 vi) ** 2 Six stationary conditions for extremum of F: dF/da1 = 0, dF/da2=0, dF/a3=0, dF/db1=0, dF/db2=0, dF/db3=0 Simplification: ------------- I assumed below the same average variance for all the samples,which is good enough if your data is reasonably correct, otherwise the model would need to deal with individual variances, and the description would not be as simple as this one. And because the variance is constant it will not appear anywhere in the equations that follow. Since variables x and y are independent the system of six equation separates itself into two system of three equation each. We have to separately solve two matrix equations, M a = p M b = q The matrix M is: | n | Sum ui | Sum vi | | | | | | Sum ui | Sum ui ui | Sum ui vi | | | | | | Sum vi | Sum ui vi | Sum vi vi | Vector p-transposed = [Sum xi, Sum xi ui, Sum xi vi] Vector q-transposed = [Sum yi, Sum yi ui, Sum yi vi] Internally, we will represent matrix M as block of three block-columns. } gps: func [ {A block of two blocks [[a1 a2 a3][b1 b2 b3]] containing estimated model parameters, computed from 4 sets of input data 'xs 'ys 'us 'vs } xs [block!] {[x1 x2 .... xn]} ys [block!] {[y1 y2 ..... yn]} us [block!] {[u1 u2 ..... un]} vs [block!] {[v1 v2 ..... vn]} /local n m p q result [block!] ][ simple-sum: func[x][foldl :+ 0 x] n: length? xs m: copy [] append/only m reduce [n simple-sum us simple-sum vs] append/only m reduce [simple-sum us product-vv us us product-vv us vs] append/only m reduce [simple-sum vs product-vv us vs product-vv vs vs] p: reduce [simple-sum xs product-vv xs us product-vv xs vs] q: reduce [simple-sum ys product-vv ys us product-vv ys vs] result: copy [] append/only result solve m p append/only result solve m q result ] comment { 'product-vv stands for scalar inner product of two vectors. 'Solve uses my home grown algorithm based on Gram-Schmidt orthogonalization procedure. A typical approach here would be to use the LU decomposition technique with pivoting. }

 [7/14] from: jjmmes::yahoo::es at: 10-Nov-2002 18:11


Jan, Many thanks for your help. Can you pls send me your functions so that I can test your gps functions I find that there are very few math related scripts available in rebol, so your scripts are most valuable! Thanks --- Jan Skibinski <[jan--skibinski--sympatico--ca]> escribió: > Jose,

 [8/14] from: jan:skibinski:sympatico:ca at: 10-Nov-2002 15:14


Jose and all,
> Many thanks for your help. Can you pls send me your > functions so that I can test your gps functions > > I find that there are very few math related scripts > available in rebol, so your scripts are most valuable! >
Jose, I just sent you the script you requested. Here are just few general comments about linear algebra scripts. I am sure several people here developed their own Linear Algebra packages in Rebol. Tom Conlin for one. I bet Joel Neely has one too. I have developed several Haskell modules that deal with some aspects of Linear Algebra, and go way beyond the basic stuff, such as LU decomposition and a solver of systems of linear equations. I needed support for nonsymmetric complex matrices, eigenvalues, eigenvectors, and all that stuff. I have also developed a Quantum Vector module which supports the <bra| |ket> Dirac's notation, and several application modules from Quantum Mechanics, one of which is a prototype of "Haskell Simulator of Quantum Computer". I did no plan however to port this stuff to Rebol, knowing well that this is not in a general interest of Rebolites. But your GPS problem forced me to port some of the stuff to Rebol. Here are few of my observations about my yesterday's porting process: ------------------------------------------------------ The process was rather easy, because of the clarity of the Haskell syntax. I did not spend much time on optimization and I used the same concepts in Rebol as in Haskell. So, with few exceptions, the recursive stuff in Haskell, remains recursive in Rebol. The original functions were by defnition "indexless" - they came from my web article "Indexless Linear Algebra Algorithms." The resulting translation is almost indexless too. The list comprehension notation , such as: p = [bra_ket u ak | ak <- a] were easy to translate to p: copy [] foreach ak a [insert tail result (bra-ket u :ak)] However, since Haskell is lazy the expressions are not evaluated until they are needed. So if during recursion, some items such as 'u cannot be evaluated because say u: first a but 'a is already empty, nothing wrong happens in Haskell because of its laziness. Not so in Rebol. The solution is to use functions, even if they represent constants. So the above 'p would be wrapped inside a function instead. p: func [ /local result ][ result: copy [] foreach ak a [insert tail result (bra-ket u :ak)] result ] Haskell has several standard words that are identical as in rebol but with completely different semantics: 'head means 'first, 'tail means 'next, etc. But one can easily get used to such things. Few things have bitten me though. I am still not used to certain "in-place" behaviour of certain functions, such as 'replace, or worst of all 'reverse, which returns empty block. This bites badly in pipelines of expressions. Regards, Jan

 [9/14] from: jjmmes:ya:hoo:es at: 10-Nov-2002 21:36


Can you pls send a link to your web article Indexless Linear Algebra Algorithms ? The link at numeric-quest seems broken ! Thanks again --- Jan Skibinski <[jan--skibinski--sympatico--ca]> escribió: >

 [10/14] from: tomc:darkwing:uoregon at: 10-Nov-2002 13:55


Hi Jose, Hmmm I recall "root mean square" from Calculus not Linear and it was the sqrt of the area under f(x^2) between a couple of points if this is the correct interpretation I am not seeing how it will help you find your constants. I think you are closer to your mark in the subject line and may be looking for a "linear-least-squares" fit. (or if you are less lucky a "non-linear-least-squares" fit) disclaimer: I have only been through standard undergraduate classes and have not used this stuff in several years (winter 99, and never in the real world). if I were you I would try to view what I needed as 2 problems each with 3 unknowns since your x and y equations seem independent and the amount of work will normally square with the number of unknowns, that is 2*3^2 < 6^2 looking at just the x half, a few gps samples you could have in something similar to your notation 1*a + lat1*b + lon1*c = x1 1*a + lat2*b + lon2*c = x2 1*a + lat3*b + lon3*c = x3 or in matrix form | 1 lat1 lon1 | |a| |x1| | 1 lat2 lon2 | * |b| = |x2| | 1 lat3 lon3 | |c| |x3| in a more standard notation this is typically represented as A x = b the linear-least-squares solution 'x' to Ax = b is x = (A_t * A)^-1 * A_t * b where A_t is the transpose of A (transpose means the indices for a row and column are swapped or that a matrix is reflected about its diagonal) and the (A_t A)^-1 means the inverse of the matrix that results from multiplying the transpose of A by A (order matters) more specifically it is the matrix you would have to multiply (A_t * A) by to get the "Identity matrix" where the identity matrix is a square matrix with 1's from the top-left to lower-right and 0's elsewhere a matrix must be square to have an inverse There is no guarantee a square matrix has an inverse. I am not going into inverting a matrix because there are alternatives. all in all it is pretty expensive to go this way but at least it is straight forward if you want to plow through it. If you are going to do allot of these there is a shortcut which stems from the rearrangement A_t A x = A_t b where the two sides of the equation can be expressed as Augmented Matrices going back to something that resembles your notation translated to matrix form we get || 1 1 1| |1 lat1 lon1| : a| | 1 1 1 : x1| ||lat1 lat2 lat3| * |1 lat2 lon2| : b| = |lat1 lat2 lat3 : x2| ||lon1 lon2 lon3| |1 lat3 lon3| : c| |lon1 lon2 lon3 : x3| where the [a b c] & [x1 x2 x3] vectors are augmenting the square matrices these augmented matrices (representing your simultaneous equations) may be solved by adding multiples of the rows in a matrix to one another. to end up with 1s, 0s, and solutions (if there are any) there are methods by Gauss, and an improvement called "Gauss-Jordan" which have been around for hundreds of years. in the end you end up with identity matrices augmented with the solutions |1 0 0 : i*a| |1 0 0 : p| |0 1 0 : j*b| = |0 1 0 : q| |0 0 1 : k*c] |0 0 1 : r] so a = p/i b = q/j c = r/k there is a function "naive-gauss" in code I wrote to double check the longhand results required in my homework. (the professor would not even allow a four function calc in the class) and was not finished beyond the point where it served my immediate purpose http://darkwing.uoregon.edu/~tomc/core/math/matrix.r that should work to reduce your augmented equations. there is another "gauss-jordan" which would be better if it did not have a bug I never got around to tracking down note: the "comp" function in the file was derived from Eric Long's compare best of luck On Sat, 9 Nov 2002, [iso-8859-1] jose wrote:

 [11/14] from: jan:skibinski:sympatico:ca at: 10-Nov-2002 16:48


Hi Jose,
> Can you pls send a link to your web article > "Indexless Linear Algebra Algorithms" ? > > The link at numeric-quest seems broken ! > > Thanks again >
I know, I broke it. The site is still defunct. But you can retrieve its past contents from the Wayback Machine. Go to www.archive.org, type in a name of a site you are interested in, and voila! Isn't that neat? While in Numeric Quest, select page with Haskell modules. The recommended way is to openit in a separtate window, so you will break out of the stupid frames. All so-called "literate modules", which look like HTML pages, used to be directly executable in Hugs interpreter. Just rename a file giving it a literate extension "lhs" rather than html and you should be able to run it - provided that Wayback Machine did not do too much damage while reformatting the pages - like substituting ">" by >, or by adding java scripts to them. But the Orthogonals.html (from the last archive) seems to be accepted by Hugs in spite of all that reformatting. Enjoy, Jan

 [12/14] from: g:santilli:tiscalinet:it at: 11-Nov-2002 11:31


Hi Jan, On Sunday, November 10, 2002, 9:14:54 PM, you wrote: JS> Few things have bitten me though. I am still not used JS> to certain "in-place" behaviour of certain functions, such JS> as 'replace, or worst of all 'reverse, which returns empty JS> block. This bites badly in pipelines of expressions. I don't know why REVERSE returns the tail of the block. You can use HEAD REVERSE in expressions, or patch it like this: use [_reverse] [ _reverse: :reverse reverse: func [ "Reverses a series." value [series! tuple! pair!] /part "Limits to a given length or position." range [integer! series!] ] [ head either part [_reverse/part value range] [_reverse value] ] ] Regards, Gabriele. -- Gabriele Santilli <[g--santilli--tiscalinet--it]> -- REBOL Programmer Amigan -- AGI L'Aquila -- REB: http://web.tiscali.it/rebol/index.r

 [13/14] from: jan:skibinski:sympatico:ca at: 11-Nov-2002 11:33


Hi Gabriele,
> I don't know why REVERSE returns the tail of the block. You can > use HEAD REVERSE in expressions, or patch it like this:
<<quoted lines omitted: 9>>
> ] > ]
Thanks, that's interesting application of 'use! But I'd rather stay away from patching and just get used to the standard behaviour of 'reverse. Jan

 [14/14] from: jjmmes::yahoo:es at: 11-Nov-2002 21:21


Hi Tom, Thanks for helping me out. I think the general approach to solving this is with Jan's code, specially when you want to acount for 4+ data points. It turns out that with 3 data points I only have 2 sets of independent equations, u: f(x,y) and v: f'(x,y), so I only have to solve one set of 3 equations with 3 unknowns and can get the unknows expressed as a function of u,v . Regards, Jose --- Tom Conlin <[tomc--darkwing--uoregon--edu]> escribió:
> Hi Jose, > Hmmm
<<quoted lines omitted: 153>>
> > > > c, d, e, f below given three calibration > points
=== message truncated ===

Notes
  • Quoted lines have been omitted from some messages.
    View the message alone to see the lines that have been omitted