Array question

New to FreeBASIC? Post your questions here.
Post Reply
Topito216
Posts: 18
Joined: Dec 11, 2009 13:13
Location: Spain

Array question

Post by Topito216 »

Hello to all.
I doing a program of the adjust of minimum least squares with matrix of size 1700x7.
I am employing this #Include "fbmath.bi" in order to get de inverse of matrix and determinant.

As is a surveyor measuring program there are some mesures that out of tolerance and that have to remove of all equations systems.
The proper is to elimiinate the line or colum of matrix an recalculate. But it's dificult I have choosed to change all the elements of the line of the matrix by all zeros.
If the system are 1400 lines I haved to remove for example 3 lines and recalculate to get more acuracy.
This way works ok and in every iteration the calculate is more acuracy.
The problem that I dont understand it's that when the program finish crash.
Only when it finish, no while it's running. I cant understand.
I know it has do to with the changing by zeros but I cant understand why. I dont know if it's a problem of the math lib or of freebasic or a mathematic concept.
Some idea.
Thanks.
vdecampo
Posts: 2992
Joined: Aug 07, 2007 23:20
Location: Maryland, USA
Contact:

Post by vdecampo »

Most likely you are writing outside of the array bounds. I can't say for sure without seeing your source code. Compile with the -exx option and run your program at the command prompt and it should show you an error message.

-Vince
dodicat
Posts: 8262
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Array question

Post by dodicat »

Hi Topito216

You can try my MATRIXINVERT function.
It is self contained, and also gives the determinant.

I've got it set up here to generate a random matrix of any size you like.
But it is easily modified to accept your own input matrix if required.

I also have similar for complex components posted on the forum
Bottom of:

http://www.freebasic.net/forum/viewtopi ... or&start=0

Here's the one for reals:

Code: Select all



Function MATRIXINVERT(matrix() As Double,inverse() As Double) As Double
 '_____________________________________________________________   
    #macro setup(the_arrays_and_variables)
    Dim As Long x2,y2  'counters for invert
    Dim As Long x,y,z  'counters for lu
    Dim det As Double=1  'determinant
    Dim mm As Double     'system dimension
    Dim m As Long=Ubound(matrix)
    Dim As Integer x1,k1,g,sign,dimcount 'some counters
    sign=1                'initial determinate sign
Redim h(1 To m) As Double         'column of a unit matrix
Redim l(1 To m,1 To m) As Double  'lower triangular
Redim b(1 To m,1 To m) As Double  'upper triangular
Redim iv(1 To m) As Double        'intermediate vector
Redim inverse(1 To m,1 To m)      'the inverse
Redim h2(1 To m,1 To m) As Double 'unit matrix
Redim solution(1 To m) As Double  'a solution column
For i As Long=1 To m
    For j As Long=1 To m
        b(i,j)=matrix(i,j)    'take a copy of the matrix
        If i=j Then 
        l(i,j)=1.0            'make two unit matrices
        h2(i,j)=1.0        
    Else
        l(i,j)=0.0
        h2(i,j)=0.0
        Endif
    Next j
Next i
#endmacro
'Everything setup, ready to go!
'____________________________________________________________
 
     #macro cleanup(all_redundant_arrays)     
Erase l,b,h,h2,iv,solution
#endmacro
'You arrays will all be annihilated at the end
'_____________________________________________________________

#macro lu(make_triangular_matrices)
For y=1 To m-1
        If b(y,y)=0 Then          
        pivot("We don't want division by zero")
        If b(y,y)=0 Then          '
            cleanup("Sorry, the matrix is singular, must go now")
        Return 0
        End If
        Endif
    For x=y To m-1
        l(x+1,y)=b(x+1,y)/b(y,y)         'l() is lower triangular matrix
        mm=l(x+1,y)
        b(x+1,y)=0
        For z=y+1 To m
            b(x+1,z)=b(x+1,z)-mm*b(y,z)  'b() is upper triangular matrix
        Next z
    Next x
Next y
#endmacro
'Upper and Lower triangular matrices all made up
'___________________________________________________________

#macro pivot(the_matrix)
For dimcount=0 To Ubound(matrix)-1
For x1=1 To Ubound(matrix)
For k1=1 To Ubound(matrix)-1
    If k1+1+dimcount >ubound(matrix) Then Exit For
    If x1+dimcount>UBound(matrix) Then Exit For
If Abs(b(k1+1+dimcount,1+dimcount))>Abs(b(x1+dimcount,1+dimcount)) Then
        sign=-1*sign   'sign changes with each swap
For g As Long=1 To Ubound(matrix)
      Swap b(k1+1+dimcount,g),b(x1+dimcount,g)  'swap lower triangle row
      Swap h2(g,k1+1+dimcount),h2(g,x1+dimcount) 'swap unit row
Next g
End If
Next k1
Next x1
Next dimcount
#endmacro
'Biggest to top, smallest to the bottom, single rank sized
'__________________________________________________________
#macro ivector(intermediate_vector)
       For n As Long=1 To m
        iv(n)=h(n)
For j As Long=1 To n-1        
iv(n)=iv(n)-l(n,j)*iv(j)        
Next j
        Next n
#endmacro
'That's the intermediate vector made up
'___________________________________________________________
#macro fvector(solution_vector)
For n As Long=m To 1 Step -1
solution(n)=iv(n)/b(n,n)
For j As Long = m To n+1 Step -1
solution(n)=solution(n)-(b(n,j)*solution(j)/b(n,n))
Next j
Next n
#endmacro
'That's the final solution vector made up
'___________________________________________________________
#macro dtm(the_determinant)
For aa As Long=1 To m
    For bb As Long=1 To m
        If aa=bb Then
   det=det *b(aa,bb)           'multiply through the trace of b()
        End If
Next bb
Next aa
#endmacro
'The deteminant is found, and we trust it is not ZERO
'___________________________________________________________
#macro invert(the_matrix)
        for x2=1 to m
        For y2  =1 To m
          h(y2)=h2(x2,y2)  
        Next y2
        ivector("Get me the intermediate vector")
        fvector("Get me the solution vector")
        For z As Long=1 To m
                inverse(z,x2)=solution(z)
        Next z
        next x2
    
    dtm("Now get me the determinant")
    cleanup("Goodbye redundant arrays")
Return sign*det
#endmacro
'ALL DONE, ready to go, just waiting on MAIN

'**************************** MAIN ***************************
    setup("dimension the arrays and variables")
    pivot("the incoming matrix")    
    lu("make lower and upper triangular matrices")   
    invert("the matrix and return the determinant")           
'************************** END MAIN *************************

End Function



'THE ABOVE FUNCTION (MATRIXINVERT) IS SELF CONTAINED,
'RETURNING THE DETERMINANT AND BUILDING THE INVERSE
'_________________________________________________________


'THE FOLLOWING TWO ROUTINES ARE FOR THE EXAMPLE:
'MATRIXMULT and DECROUND
'ALTHOUGH
'MATRIXMULT IS A SELF CONTAINED MATRIX MULTIPLIER
'AND
'DECROUND IS A ROUNDING FUNCTION

   Sub matrixmult(m1() As Double,m2() As Double,ans() As Double)
   Dim rows As Long=Ubound(m1,1)
   Dim columns As Long=Ubound(m2,2)
   If Ubound(m1,2)<>UBound(m2,1) Then 
           Print "Can't do"
           Exit Sub
   Endif
   Redim ans(1 To rows,1 To columns)
   Dim rxc As Double
    For r As Long=1 To rows
        For c As Long=1 To columns
        rxc=0
        For k As Long = 1 To Ubound(m1,2)
            rxc=rxc+m1(r,k)*m2(k,c)'kj
        Next k
        ans(r,c)=rxc
        Next c
        Next r
    End Sub
    Function decround ( a As Double, b As Integer) As Double 
Dim y As Double
Dim i As Double
Dim r As Double
y = (Abs(a) - Int(Abs(a))) * (10 ^ b)
i = Int(y)
y = y - i
If y >= .5 Then i = i + 1
i = i / (10 ^ b)
r = Int(Abs(a)) + i
If a < 0 Then r = -r
decround = r
End Function

'INSTRUCTIONS:
'1)  MATRIXINVERT
'   Redim an empty array dimension 2 as a double to hold the inverse
'   e.g. redim myinverse(0,0) as double
'   dim or redim an array dimension 2 to accept the matrix you want to enter
'  e.g. redim mymatrix(1 to 4,1 to 4) as double
'  Call the function e.g. print MATRIXINVERT(mymatrix(),myinverse())
'   The function returns the determinant and fills myinverse() with the inverse

'2)  MATRIXMULT
'   Redim an empty array dimension 2 as a double to hold the product
'   e.g. redim answer(0,0) as double
'   as above, prepare two matrices to multiply
'   e.g. redim m1(1 to 5,1 to 3) as double
'        redim m2(1 to 3,1 to 4) as double
'   call the sub thus: matrixmult m1(),m2(),answer()
'   the product will be held in answer in this case answer(1 to 5,1 to 4)
'   NOTE:
'   To multyply say a single vector dim 4 by a matrix then 
'   the maths format applies i.e 
'   redim mat(1 to 4,1 to 4) as double (4 rows,4 columns)
'   redim vector(1 to 4,1 to 1) as double (4 rows,1 column)
'   redim answer(0,0) as double
'   matrixmult mat(),vector(),answer()  to fill the answer
'3) DECROUND Self explanitory e.g. decround(7.879965544,3)
    
    '************ EXAMPLE********************
    Function rr(first As Double, last As Double) As Double 'random generator
    Function = Rnd * (last - first) + first
End Function
    Dim m As Integer
start:
'ENTER A MATRIX
Input "ENTER DIMENSION OF MATRIX or (0 TO END)  ";m
If m=0 Then End
Print
'print "Please enter the components one by one"
Print
Redim b(1 To m,1 To m) As Double  'main matrix
Redim myinverse(0,0) As Double    'THE FUNCTION WILL FILL THIS

'MATRIX INPUT_________________________________
For aa As Integer=1 To m
    For bb As Integer=1 To m
       ' Print "Enter (";aa;",";bb;")"; 'OPTIONAL MANUAL INPUTS, uncomment the three bits
        'Input b(aa,bb)
        b(aa,bb)=rr(-1,1)
    Next bb
    'Print
Next aa
'_____________________________________________
Print
'DETERMINANT OUTPUT
Dim db As Double
db=MATRIXINVERT (b(),myinverse())
Print "Determinant =  ";db

'INVERSE OUTPUT__________________________________
'Print "INVERSE"
For aa As Integer=1 To m
For bb As Integer=1 To m
       ' Print "(";aa;",";bb;") = ";myinverse(aa,bb)  'optional printout of inverse    
Next bb
'Print
Next aa
'___________________________________________________

'NOW MULTIPLY MATRIX BY IT'S INVERSE TO CHECK
Redim product(0,0) As Double   
matrixmult b(),myinverse(),product()
Print "CHECK, MATRIX X INVERSE - SHOULD BE A UNiT MATRIX"
Print "If not then the matrix is singular,press a key"
sleep
Print
    For aa As Integer=1 To m
For bb As Integer=1 To m
        Print  decround(product(aa,bb),8);  'round to 8 places     
Next bb
Print
Next aa
Print
Goto start
   
integer
Posts: 410
Joined: Feb 01, 2007 16:54
Location: usa

Post by integer »

Are the lines that you "zero out" included in the least square calculation?
They should not be, but your statement seems to indicate that it is.

I have worked with many 2-dim least squares, and some 3-dim, however, I have never seen a least square using 7 dimensions.
Are ALL seven dimensions equally weighed?

Does the program work with only four, five or six lines of data?

Then at what point (before you use 1400 lines) does it crash?
Richard
Posts: 3096
Joined: Jan 15, 2007 20:44
Location: Australia

Post by Richard »

Maybe to eliminate an observation from the matrix, swap it with the good observation at the end and reduce the count. Repeat until all bad data is eliminated. Then solve the smaller system without the bad data.
TJF
Posts: 3809
Joined: Dec 06, 2009 22:27
Location: N47°, E15°
Contact:

Post by TJF »

@Topito216

You know about libFBla (FreeBasic Linear Algebra Library). Why don't you use it?
Richard wrote:Maybe to eliminate an observation from the matrix, swap it with the good observation at the end and reduce the count. Repeat until all bad data is eliminated. Then solve the smaller system without the bad data.
Using libFBla the swap method can shift the rows to the end and the redim method will reduce the matrix size. Then recalculate.
fxm
Moderator
Posts: 12564
Joined: Apr 22, 2009 12:46
Location: Paris suburbs, FRANCE

Post by fxm »

TJF wrote:... and the redim method will reduce the matrix size.
WARNING (extract of documentation):
With multiple dimensions, only the upper bound of only the first dimension may be safely increased. If the first dimension is reduced, the existing mappable data may be lost.
TJF
Posts: 3809
Joined: Dec 06, 2009 22:27
Location: N47°, E15°
Contact:

Post by TJF »

fxm wrote:
TJF wrote:... and the redim method will reduce the matrix size.
WARNING (extract of documentation):
With multiple dimensions, only the upper bound of only the first dimension may be safely increased. If the first dimension is reduced, the existing mappable data may be lost.
I didn't refer to the FB command REDIM. I spoke about the Redim_ method of the LA_M type. The said limitations are not valid for this method. Anyway, a matrix of that size shouldn't be used at all.


@Topito216

You shouldn't use a 1700x7 matrix. It's better to build the Gradient equations direct from the data (coefficient matrix and solution vector).

Here's an example

Code: Select all

#DEFINE LA_ScalarDouble
#INCLUDE "libFBla.bas"


' this function computes the value of polynom P at position X
FUNCTION Polynom(BYVAL X AS LA_S, BYVAL P AS LA_V PTR) AS LA_S
  WITH *P
    VAR r = .Val_(0)
    FOR i AS INTEGER = 1 TO .Ubound_()
      r += .Val_(i) * X ^ i
    NEXT
  END WITH : RETURN r
END FUNCTION

' this function computes the polynom regression
FUNCTION GetPolynom(BYVAL Order AS UINTEGER, _
                    BYVAL X AS LA_V PTR, BYVAL Y AS LA_V PTR) AS LA_V
  ' prepare variables for the gradient equations
  ' coefficient matrix m and solution vector n
  VAR m = LA_M(Order + 1), n = LA_V(order + 1)

  ' compute matrix (one triangle) and vector elements
  FOR k AS INTEGER = 0 TO X->Ubound_()
    FOR i AS INTEGER = 0 TO Order
      FOR j AS INTEGER = 0 TO i
        m.set_(j, i, m.val_(j, i) + X->val_(k) ^ (i + j))
      NEXT
      n.val_(i) = n.val_(i) + Y->val_(k) * X->val_(k) ^ i
    NEXT
  NEXT
  ' mirror triangle matrix
  FOR i AS INTEGER = 1 TO Order
    FOR j AS INTEGER = 0 TO i - 1
      m.set_(i, j, m.val_(j, i))
    NEXT
  NEXT

  ?"Gradient equations:"
  ?m
  ?
  ?"Gradient vector:"
  ?n
  ' compute the polynom parameters Beta
  RETURN n / m
END FUNCTION


' here is the data input (both vectors must have the same size)
' the positions X
VAR x = LA_V("-4, -3, -2, -1, 0, 1, 2, 3, 4")
' and the measurement values Y
VAR y = LA_V("0.03, 0.094, 0.213, 0.388, 0.643, 0.987, 1.426, 1.969, 2.632")


' you may consider to scale your data into the interval [-1, 1] here to
' avoid numerical issues.


' the order of the polynom
VAR order = 3

' compute the polynom parameters Beta
VAR beta = GetPolynom(order, @x, @y)
?
?"Polynom parameters:"
?beta

' compute standard deviation
VAR stabw = 0.0
FOR i AS INTEGER = 0 TO x.Ubound_()
  stabw += (Polynom(x.Val_(i), @beta) - y.Val_(i)) ^ 2
NEXT
?
?"Standard deviation:"
?stabw

?
?"Press a key to reduce the set of data."

SLEEP

' the record no 3 will be removed (last record is swapped to pos 3)
VAR i = 3
x.Val_(i) = x.Val_(x.Ubound_())
x.Redim_(x.Ubound_() - 1)
y.Val_(i) = y.Val_(y.Ubound_())
y.Redim_(y.Ubound_() - 1)

' compute the polynom parameters Beta again
VAR beta1 = GetPolynom(order, @x, @y)
?
?"Polynom parameters:"
?beta1

' compute standard deviation again
stabw = 0.0
FOR i AS INTEGER = 0 TO x.Ubound_()
  stabw += (Polynom(x.Val_(i), @beta1) - y.Val_(i)) ^ 2
NEXT
?
?"Standard deviation:"
?stabw

SLEEP
Topito216
Posts: 18
Joined: Dec 11, 2009 13:13
Location: Spain

Least squares adjust

Post by Topito216 »

First of all I have to thank all of you for answer to me so fast. It's fantastic to know there are people ready to help. I'll try to answer all your post in this with post with my badl english lenguage. I can write a little better but I'll need double time to do it.

The courious is that the crash happens when the program exit and It did the job ok. No while is running , jus to exit.

Vdecampo:
You are allright. The problem was that I try to write outside the range of a unidimensional matrix. Your help "-exx" was very usefull because it detetec 4 fails that are not the origin of the crash program but show wrong information.
like this:

Code: Select all

Dim As Double X(1 To 7,1 To 1)
Print "betta= ";x(2,1),,x(1,2)*180/PI;" grados sexa"
The second number for to see the angle in degres mustbe x(2,1).
But this was not the crash origin.
The origin of the probglem was to write at jacobiana(0, 1 to 7) when elimina(n)=0:

Code: Select all

Dim elimina(1 To numerodelineas*3) As Integer
Dim As Double jacobiana (1 To numerodelineas*3,1 To 7)
For n=1 To UBound (elimina,1)
jacobiana(elimina(n),1)=0
jacobiana(elimina(n),2)=0
jacobiana(elimina(n),3)=0
jacobiana(elimina(n),4)=0
jacobiana(elimina(n),5)=0
jacobiana(elimina(n),6)=0
jacobiana(elimina(n),7)=0
independiente(elimina(n),1)=0
End If
Next n
This extension for to compile -exx is very usefull if you know more tricky like that or the site where to learn please tell me. Thanks.

Dodicat
Your program works ok and I'll have in mind because it's seem easy to understand but my knoledge about FB is low and there are commands like #macro that I dont understand but with the time and writing and writing i'll learn more.


integer

To change a entire line or column to zeros works ok and it's the same efect that remove the entire line. I have confirm it doing test.
I was wrong the problem was not to operate matrix with entires lines of zeros. The problem was to write outside the limits definition of matrix.

The matrix are bidimensional. It's a system model of the movement of a radiotelescope with horizontal and vertical movementm, even torsion, displacement of parabolic,intersection of elevation an horizontal axis, etc.
Are 7 incognites and thousands of mesures from a taquimeter that coordinates of retroreflector cubes at radio telescope.

As you can imagine the model is:
AX = T
where:
A = is the Jacobian matrix of partial derivatives with respect to the parameters. Dimension (k , 7). k maybe 10 or 1000000
X = x − x0 is the vector of estimated corrections to the a priori values of the parameters, Dimension (7,1). Incognites to solve
T = L − F(x0) is the vector of residuals, that is, the observed minus computed observations. (k,1)
Extracting the parameters to estimate from equation 6 and weighting the input observations
results in:
X = (A'WA)^-1*(A'WT)
This is a tipical MinimumLeastSquares.
As you see I have to multiply, add, transpose and invert matrix.

I have calculate the system with (k) 1000000 measures and works perfect.

Richard

It's a good idea and tried to do it but finally change entires by zeros works ok and I did not to have to swaping lines that it's a little more dificult.

TJF

I know your library But I AM SORRY , my low level of knowledge of functions, librarys and methods make imposible to understand your library. I have only the file libFBla.bas but it's not enough for me.

Your example:

Code: Select all

#INCLUDE "libFBla.bas"

#DEFINE english

' used to print out headers and values
#MACRO HEADER(_G_, _E_, _T_)
 #IFDEF english
   ?:?_E_
 #ELSE
   ?:?_G_
 #ENDIF
 ?_T_
#ENDMACRO


HEADER("Loesung eines linearen Gleichungssystemes", _
       "Solution of a linear equation system", "")
?"   3x - 3y + z = 0"
?"      + 4y - z = 5"
?"   2x - 2y + z = 1"

VAR A = LA_M("3, -3,  1" NL _
             "0,  4, -1" NL _
             "2, -2,  1")
HEADER("Koeffizientenmatrix A:", _
       "Coefficient matrix A:", A)

VAR B = LA_V("0, 5, 1")
HEADER("Loesungsvektor B:", _
       "Solution vector B:", B)

VAR r = B / A
HEADER("Ergebnisvektor r = B / A:", _
       "Result vector r = B / A:", r)

HEADER("Ergebnisse:", _
       "Results:", "")
?"x = ";r.Val_(0)
?"y = ";r.Val_(1)
?"z = ";r.Val_(2)



HEADER("Matrizenrechnung", _
       "Matrix algebra", "")
VAR m1 = LA_M( _
  ".11, .12, .13" NL _
  ".21, .22, .23")
?"m1:" : ?m1

VAR m2 = LA_M( _
  "11, 12" NL _
  "21, 22" NL _
  "31, 32")
?"m2:" : ?m2

VAR m3 = m1 * m2
?"m3 = m1 * m2:" : ?m3

#IFNDEF __FB_UNIX__
SLEEP
#ENDIF
I't very easy to understant but:
How pick to a VAR position like as matrix (example (m(5,4)) in order to change the value at that position
How to inver the variable m2 or m1
How calculate determinant
How to transpose,
etc.

I know it's not too dificult but I could not calculate for example the determinant of a matrix. Without examples it's very dificult for me to understand how to use your library. I know it's my problem and I apreciate your help but I'll have to wait untill my knowledge upgrade.

That's the reason because I use th FBmath. It has a pdf file with explanation and examples to operate with matrix.
In an other hand in order to add,rest,multiply, I have create the following subroutines. Are very simply but they works:

Code: Select all

'subrutinas matrices

Sub TRASPONER(MA()As Double,MT()As Double)
	' TRASPONE UNA MATRIZ DE ENTRADA EN UNA DE SALIDA
	Dim As Integer I,J
		For I=1 To UBound (MA,2)
		For J=1 To UBound (MA,1)
			 MT(I,J)=MA(J,I)
		Next
	Next
End Sub

Sub SUMA(MA() As Double,MB() As double)
	' SUMA DOS MATRICES Y EL RESULTADO SE CARGA EN LA PRIMERA
	If UBound (MA,1)<>UBound (MB,1) Then Print "NO SE PUEDEN SUMAR, LA MATRIZ ES A": Exit Sub
	If UBound (MA,2)<>UBound (MB,2) Then Print "NO SE PUEDEN SUMAR, LA MATRIZ ES A": Exit Sub
	Dim As Integer I,J
	For I=1 To UBound (MA,1)
		For J=1 To UBound (MA,2)
			 MA(I,J)=MA(I,J)+MB(I,J)
		Next
	Next
End Sub
Sub RESTA(MA() As Double,MB() As Double) 
	' RESTA DOS MATRICES Y EL RESULTADO SE CARGA EN LA PRIMERA
	If UBound (MA,1)<>UBound (MB,1) Then Print "NO SE PUEDEN RESTAR, LA MATRIZ ES A": Exit Sub
	If UBound (MA,2)<>UBound (MB,2) Then Print "NO SE PUEDEN RESTAR, LA MATRIZ ES A": Exit Sub
	Dim As Integer I,J
	For I=1 To UBound (MA,1)
		For J=1 To UBound (MA,2)
			 MA(I,J)=MA(I,J)-MB(I,J)
		Next
	Next
End Sub
Sub PINTAMATRIZ(MA() As Double)
	Dim As Integer I,J
	Print "Rango: ";UBound(ma,1),UBound(ma,2)
	
	If ubound(MA,2)>1 Then	
		For I=1 To UBound (MA,1)
		For J=1 To UBound (MA,2)
			Print MA(I,J),
		Next
		Print ""
	Next
	End If
	
	If ubound(MA,2)=1 Then	'ATENCION A LAS UNIDADES 10000
		For J=1 To UBound (MA,1)
			Print MA(J,1)
		
		Next
		Print ""
		End If
	
End Sub
Sub MULTIPLICA(MA() As Double,MB() As Double,MC()As Double) 
	'MULTIPLICA LAS DOS MATRICES Y EL RESULTADO SALE POR LA TERCERA
Dim As Integer I,J,K
	For I=1 To UBound(MA,1)
		For J=1 To UBound(MB,2)
			MC(I,J)=0
			For K=1 To UBound(MA,2)
				 MC(I,J)=MC(I,J)+MA(I,K)*MB(K,J)
			Next K
		Next J
	Next I
End Sub
Sub MULTIPLICAV(MA() As Double,V() As Double,MC()As Double) 
	'MULTIPLICA LAS DOS MATRICES Y EL RESULTADO SALE POR LA TERCERA
Dim As Integer I,J,K
	For I=1 To UBound(MA,1)
	
			MC(I,1)=0
			For K=1 To UBound(MA,2)
				 MC(I,1)=MC(I)+MA(I,K)*V(K,1)
			Next K
		
	Next I
End Sub

For inverse and determinant and all the operations I have choose the fbmath routine and my own subroutines:

Code: Select all



TRASPONER (jacobiana(),AT())
MULTIPLICA (AT(),jacobiana(),Normal())

GaussJordan Normal(), Det ' CALCULA EL DETERMINANTE DE LA MATRIZ Y LA INVERSA SE CARGA EN LA MATRIZ DE ENTRADA
SELECT CASE MathErr
  CASE MatOk
    'PINTAMATRIZ (Normal())
   

    'PRINT : PRINT "Determinant =", Det
    'Print ""
    Sleep retraso
    'GaussJordan NORMAL(), Det
    IF MathErr = MatOk THEN
      'PINTAMATRIZ (NORMAL())
    END IF
  CASE MatErrDim
    PRINT "Non-compatible dimensions!"
  CASE MatSing
    PRINT "Quasi-singular matrix!"
END Select

MULTIPLICA (AT(),independiente(),T())
MULTIPLICA (Normal(),T(),X())
MULTIPLICA (JACOBIANA(),X(),V())
RESTA (V(),INDEPENDIENTE())
The program works very well finally and I wonder when it shows results from operations of complicated formulas.



The next step is to use Robust Estimators for residuals in order to have a good criterion to remove bad observations

Thanks to all again.

Bye[/code]
Richard
Posts: 3096
Joined: Jan 15, 2007 20:44
Location: Australia

Post by Richard »

@ Topito216.
I have some questions because I have an interest in antennas.

What is the diameter and focal length of the parabola you are measuring?
What is the maximum frequency ( ? GHz) you will observe with the antenna?

You have many reflectors on the surface of the parabola, do you place the tachymeter;
(1) fixed on the ground nearby, (2) at the vertex, or (3) at the focus ?

Can you produce a phase error map (= +/-mm) of the parabola surface by using coherent detection of satellite signals ?
Topito216
Posts: 18
Joined: Dec 11, 2009 13:13
Location: Spain

Radio Telescope

Post by Topito216 »

I answer:

Richard wrote:@ Topito216.
I have some questions because I have an interest in antennas.

What is the diameter and focal length of the parabola you are measuring?
I dont know the focal length but,
Nasmyth-Cassegrain
Diameter principal: 40m
Diámeter secundario: 3.2m




What is the maximum frequency ( ? GHz) you will observe with the antenna? 22 Gz

You have many reflectors on the surface of the parabola, do you place the tachymeter;
(1) fixed on the ground nearby, (2) at the vertex, or (3) at the focus ? The goal is not mesure the paraboloid surface, it's to mesure the dimensions of antenna and the invariant reference point.

Can you produce a phase error map (= +/-mm) of the parabola surface by using coherent detection of satellite signals ?
The holography technique is used to get a model of the surface of the antenna with acuracy of 1e-5 meters. Now people is working on it but there is not results yet.
TJF
Posts: 3809
Joined: Dec 06, 2009 22:27
Location: N47°, E15°
Contact:

Post by TJF »

Topito216 wrote:How pick to a VAR position like as matrix (example (m(5,4)) ...
libFBla.bas, line 580 ff wrote:' Get value at Row and Col
FUNCTION LA_M.Val_(BYVAL R AS UINTEGER, BYVAL C AS UINTEGER) AS LA_S
...
Topito216 wrote:... in order to change the value at that position
libFBla.bas, line 595 ff wrote:' Set a new Value at Row and Col
SUB LA_M.Set_(BYVAL R AS UINTEGER, BYVAL C AS UINTEGER, _
BYVAL V AS LA_S = .0)
...
Topito216 wrote:How to inver the variable m2 or m1
libFBla.bas, line 842 ff wrote:' Creates the inverse of a matrix
FUNCTION Inv(BYREF M AS LA_M) AS LA_M
...
m2 cannot be inverted (3 x 2 matrix).
Topito216 wrote:How calculate determinant
libFBla.bas, line 611 ff wrote:' Compute the determinat (returns sub-determinant if Ro < Co)
FUNCTION LA_M.Det_() AS LA_S
...
Topito216 wrote:How to transpose,
libFBla.bas, line 751 ff wrote:' Creates a transposed matrix
FUNCTION Transp OVERLOAD(BYREF M AS LA_M) AS LA_M
...
etc ... Be honest, did you look into the file at all?


Here's a libFBla version of your code:

Code: Select all

#INCLUDE ONCE "libFBla.bas"

' No need for
' DIM AS DOUBLE jacobiana(...),AT(...),Normal(...),independiente(...),T(...),X(...),V(...),Det

' ...

' TRASPONER (jacobiana(),AT())
VAR AT = Transp(jacobiana)

' MULTIPLICA (AT(),jacobiana(),Normal())
VAR Normal = AT * jacobiana

'GaussJordan Normal(), Det ' CALCULA EL DETERMINANTE DE LA MATRIZ Y LA INVERSA SE CARGA EN LA MATRIZ DE ENTRADA
'SELECT CASE MathErr
  'CASE MatOk
    ''PINTAMATRIZ (Normal())


    ''PRINT : PRINT "Determinant =", Det
    ''Print ""
    'Sleep retraso
    ''GaussJordan NORMAL(), Det
    'IF MathErr = MatOk THEN
      ''PINTAMATRIZ (NORMAL())
    'END IF
  'CASE MatErrDim
    'PRINT "Non-compatible dimensions!"
  'CASE MatSing
    'PRINT "Quasi-singular matrix!"
'END Select
Normal = Inv(Normal) ' on error output by libFBla

' MULTIPLICA (AT(),independiente(),T())
VAR T = AT * independiente

' MULTIPLICA (Normal(),T(),X())
VAR X = Normal * T

' MULTIPLICA (jacobiana(),X(),V())
VAR V = jacobiana * X

' RESTA (V(),independiente())
V -= independiente
Sorry, no parentheses here. And no indices to get out of range.

Instead of

Code: Select all

Normal = Inv(Normal)
' ...
VAR X = Normal * T
u can use

Code: Select all

VAR X = T / Normal ' (as shown in the example)
But OK, if you think your version is easier to understand ... and to maintain ...


Soon, I'll upload a new version with smallish bugfixes regarding formated output.
Topito216
Posts: 18
Joined: Dec 11, 2009 13:13
Location: Spain

Post by Topito216 »

TFJ.

Thanks for your help.

Believe me if I tell you that I tried during 2 hours run your library but it was imposible for me. I tried the first lines of the bas file as:

Code: Select all

#DEFINE LA_ScalarDouble
#INCLUDE "libFBla.bas"
VAR B = LA_V("0, 5, 1")
Var C = LA_V("4,5,3")
Var D= LL_Che(B,C)
Print LL_Che(B,C)
Cause I could not to make work the code I stoped and decided to make my own subroutines beacuse ther are easy and for the moment I don need more complicate routines.


With your news examples I realized that it's easy to use your library and I can use sub and functions for the momento but I dont understanf the use of macros....

For example that was easy:

Code: Select all

#DEFINE LA_ScalarDouble
#INCLUDE "libFBla.bas"
VAR B = LA_V("0, 5, 1")
Var C = LA_V("4,5,3")
Print Dist(B,C)
Print (Angle(b,c))

sleep
I found imcompatibility between libfbla.bas and fbmath.bi.
I use fbmath to get the value of constant pi to convert the angle in radians to degres multipliying by 180/pi

Code: Select all

#DEFINE LA_ScalarDouble
#INCLUDE "libFBla.bas"
#Include "fbmath.bi"
This error:
Program Files (x86)\FreeBASIC\inc\fbmath.bi(387) error 4: Duplicated definition in 'BYREF Det AS DOUBLE)'

I'd like to use both libraries because I need also some functions from fbmath.bi.

I can't this works? Operators and properties.

Code: Select all

#DEFINE LA_ScalarDouble
#INCLUDE "libFBla.bas"
'#Include "fbmath.bi"

Dim s As String


VAR B = LA_V("0, 1, 1")
Var C = LA_M("0,2,1")

s=C.CAST()
print B.Len_() 

sleep
An other question is that the matrix that I use are dimensioned from 1 to x. Then the index 1 is the first position at row or column but your definition of matrix is with the first position at index 0 like r.val(0) expresion. It's posible to change without problems?

One of the bigger problem I found, was I dont know how to create an matrix or vector empty but with dimensions determined.

I tried something like tht but it does not worked.

Code: Select all

#DEFINE LA_ScalarDouble
#INCLUDE "libFBla.bas"

Dim As LA_V P
Dim As New LA_V P
P.SET_(1)=8
Print P.VAL_(1)
Print P.LEN_()

Sleep


I think that with a little bit more help I'll get to take advantage of that complete library.

Best regards.
TJF
Posts: 3809
Joined: Dec 06, 2009 22:27
Location: N47°, E15°
Contact:

Post by TJF »

Topito216 wrote:I think that with a little bit more help I'll get to take advantage of that complete library.
It's not only you who needs help. I as a developer need some feedback to find weak points. Your last post is helpful for me.
Topito216 wrote:... it was imposible for me. I tried the first lines of the bas file as:

Code: Select all

#DEFINE LA_ScalarDouble
#INCLUDE "libFBla.bas"
VAR B = LA_V("0, 5, 1")
Var C = LA_V("4,5,3")
Var D= LL_Che(B,C)
Print LL_Che(B,C)
Cause I could not to make work the code I stoped and decided to make my own subroutines beacuse ther are easy and for the moment I don need more complicate routines.


With your news examples I realized that it's easy to use your library and I can use sub and functions for the momento but I dont understanf the use of macros....
LL_Che is a (multi-line) MACRO for internal use (and it was badly documented). At your skill-level you shouldn't use it.

To get the length of a vector use the property Len_ (like you did in some code: VectorName.Len_() ).
Topito216 wrote:I found imcompatibility between libfbla.bas and fbmath.bi.
I use fbmath to get the value of constant pi to convert the angle in radians to degres multipliying by 180/pi

Code: Select all

#DEFINE LA_ScalarDouble
#INCLUDE "libFBla.bas"
#Include "fbmath.bi"
This error:
Program Files (x86)\FreeBASIC\inc\fbmath.bi(387) error 4: Duplicated definition in 'BYREF Det AS DOUBLE)'

I'd like to use both libraries because I need also some functions from fbmath.bi.
Thanks for the report! That's because both libraries have a FUNCTION named GaussJordan. I fixed this in version 0.2. It works when you include libFBla.bas before fbmath.bi (as it is in your code).

To fix it with your version you can cover one of the libraries by a NAMESPACE.
Topito216 wrote:I can't this works? Operators and properties.

Code: Select all

#DEFINE LA_ScalarDouble
#INCLUDE "libFBla.bas"
'#Include "fbmath.bi"

Dim s As String


VAR B = LA_V("0, 1, 1")
Var C = LA_M("0,2,1")

s=C.CAST()
print B.Len_() 

sleep
Sorry, Var C = LA_M("0,2,1") doesn't work due to a bug. I didn't test to declare a matrix with only one row yet. Thanks for your report! It's fixed in version 0.2. (Workaround: VAR C = Transp(LA_V("0,2,1")))

Instead of s=C.CAST() use s=C.
Topito216 wrote:An other question is that the matrix that I use are dimensioned from 1 to x. Then the index 1 is the first position at row or column but your definition of matrix is with the first position at index 0 like r.val(0) expresion. It's posible to change without problems?
I don't see any reason for starting an index at 1 (, exept your personal preference). But I know lots of technical reasons why an index should start at zero (, see also other libraries like ATLAS).

You can change it without problems by using a macro to shift the index like

Code: Select all

VAR B = LA_V("1, 2, 3")

#DEFINE Index(N) (N - 1) ' or use a short name like I_
?B.VAL_(Index(2)) ' should output 2
Topito216 wrote:One of the bigger problem I found, was I dont know how to create an matrix or vector empty but with dimensions determined.

I tried something like tht but it does not worked.

Code: Select all

#DEFINE LA_ScalarDouble
#INCLUDE "libFBla.bas"

#DEFINE LA_ScalarDouble
#INCLUDE "libFBla.bas"

Dim As LA_V P
Dim As New LA_V P
P.SET_(1)=8
Print P.VAL_(1)
Print P.LEN_()

Sleep#DEFINE LA_ScalarDouble
#INCLUDE "libFBla.bas"

Dim As LA_V P
Dim As New LA_V P
P.SET_(1)=8
Print P.VAL_(1)
Print P.LEN_()

Sleep
Dim As New LA_V P
P.SET_(1)=8
Print P.VAL_(1)
Print P.LEN_()

Sleep
  • When you declare a variable of type LA_V it has 1 element by default ( Dim As LA_V P ).
  • When you need a vector with more than one element just declare like VAR P = LA_V(8) and get 8 elements in this case (index 0 to 7, all elements are set to 0.0 by default).
  • When you like to set initial values use VAR P = LA_V(8, -1) to set all 8 elements to -1 (version 0.2).
  • Find more about this in the comments before the CONSTRUCTORS of LA_V and LA_M.
  • As in your code: to get a value in a vector use ? P.Val_(0) (PROPERTY in LA_V type).
  • To set a value in a vector use P.Val_(0) = 8 (PROPERTY in LA_V type).
  • Getting a matrix value looks similar: ? MatrixName.Val_(Row, Column). But instead of a PROPERTY here a FUNCTION is used. This is neccessary because in FreeBasic a PROPERTY can have only one index (not useable for a matrix).
  • The limitation of a FreeBasic-PROPERTY is the reason why setting a matrix value is different: MatrixName.Set_(Row, Column, Value).
I think I can upload version 0.2 today at the german webside.
Post Reply