A linear curve / equation finder

Share your advanced PureBasic knowledge/code with the community.
infratec
Always Here
Always Here
Posts: 6817
Joined: Sun Sep 07, 2008 12:45 pm
Location: Germany

A linear curve / equation finder

Post by infratec »

I needed something to generate a linear equation from a few points to calc all other values.

Exactly I get a Zeitronix SDL file with Lambda values and in the Zeitronix software I can see calculated values and in the file were
only analog voltages wich dos not correspond to the voltage of the Lambda probe.

So I picked up a few points and needed to calculate a linear equation for this points.

First I searched in the net and found:
https://www.dcode.fr/function-equation-finder

But I thougt it would be nice to have this as PB code.

I used the calulation from:
https://blog.mbedded.ninja/mathematics/ ... e-fitting/

and borrowed the inverse matrix calculation from wilbert:
https://www.purebasic.fr/english/viewto ... 04#p523704

Combined, it results in LinearCurveFinder.pbi

Code: Select all

;
; General calulation of a linear curve fitting
; https://blog.mbedded.ninja/mathematics/curve-fitting/linear-curve-fitting/
;
; Calculation of an inverse matrix by wilbert
; https://www.purebasic.fr/english/viewtopic.php?p=523704#p523704
; 
; Topic
; https://www.purebasic.fr/english/viewtopic.php?p=586282#p586282
;

CompilerIf #PB_Compiler_IsMainFile
  EnableExplicit
CompilerEndIf



Procedure.i InverseMatrixDouble(*Matrix, *InvMatrix, n)              ; n x n matrix inverse
  
  Protected.l n1=n-1,p,k,t,v
  Protected.d r,m,*s.Double,*d.Double
  Protected Dim Mat.d(n1,n1)
  Protected Dim Piv.u(n1)
  
  
  CopyMemory(*Matrix, @Mat(), n*n*SizeOf(Double))               ; << copy *Matrix to Mat() >>
  
  For t=0 To n1 : Piv(t)=t : Next                               ; << init Piv() >>
  
  For p=0 To n1 : m=0
    For t=p To n1                                               ; << find highest absolute value >>
      r=Abs(Mat(p,Piv(t))) : If r>m : m=r : v=t : EndIf
    Next
    If m=0 : ProcedureReturn #False : EndIf                          ; << exit if all 0 >>
    Swap Piv(p),Piv(v) : k=Piv(p) : r=1/Mat(p,k) : Mat(p,k)=1        ; << get k and reciprocal for Mat(p,k) >>
    *d=@Mat(p,0) : For t=0 To n1 : *d\d*r : *d+8 : Next              ; << multiply pivot row by reciprocal >>
    *d=@Mat() 
    For t=0 To n1
      If t<>p
        *s=@Mat(p,0) : m=-Mat(t,k) : Mat(t,k)=0
        For v=0 To n1 : *d\d+*s\d*m : *s+8 : *d+8 : Next        ; << calculate new Mat() values >>
      Else
        *d+n<<3
      EndIf
    Next
  Next
  For t=0 To n1 : *d=*InvMatrix+Piv(t)*n<<3
    For v=0 To n1 : *d\d=Mat(t,Piv(v)) : *d+8 : Next            ; << set *InvMatrix >> 
  Next
  
  ProcedureReturn #True
  
EndProcedure


Procedure.i LinearCurveFinder(List PointList.Point(), *a.double, *b.Double)
  
  Protected.i Result
  Protected.d A11, A12, A21, A22, B11, B21
  Protected Dim DoubleArray.d(1, 1)
  
  
  ForEach PointList()
    A11 + PointList()\x * PointList()\x
    A12 + PointList()\x
    B11 + PointList()\x * PointList()\y
    B21 + PointList()\y
  Next
  A21 = A12
  A22 = ListSize(PointList())
  
;   Debug "A11: " + Str(A11)
;   Debug "A12: " + Str(A12)
;   Debug "B11: " + Str(B11)
;   Debug "B21: " + Str(B21)
;   Debug "A22: " + Str(A22)
  
  ; [ A11  A12 ] [ a ] = [ B11 ]
  ; [ A21  A22 ] [ b ] = [ B21 ]
  
  DoubleArray(0, 0) = A11
  DoubleArray(0, 1) = A21
  DoubleArray(1, 0) = A12
  DoubleArray(1, 1) = A22
  
  Result = InverseMatrixDouble(@DoubleArray(), @DoubleArray(), 2)
  If Result
    *a\d = DoubleArray(0, 0) * B11 + DoubleArray(0, 1) * B21
    *b\d = DoubleArray(1, 0) * B11 + DoubleArray(1, 1) * B21
  EndIf
  
  ProcedureReturn Result
  
EndProcedure


CompilerIf #PB_Compiler_IsMainFile
  
  #Width = 800
  #Height = 600
  
  Define.d a, b
  NewList PointList.point()
  Define.d MaxX, MaxY, MinX, MinY
  Define.d XFactor, YFactor, Factor
  Define.i FactorIsLargerThanOne, x, Event
  Define Equation$
  
  
  CompilerIf #True
    
    AddElement(PointList())
    PointList()\x = 2
    PointList()\y = 1
    
    AddElement(PointList())
    PointList()\x = 3
    PointList()\y = 7
    
    AddElement(PointList())
    PointList()\x = 6
    PointList()\y = 5
    
    AddElement(PointList())
    PointList()\x = 7
    PointList()\y = 9
    
  CompilerElse
    
    AddElement(PointList())
    PointList()\x = 420
    PointList()\y = 23
    
    AddElement(PointList())
    PointList()\x = 430
    PointList()\y = 25
    
    AddElement(PointList())
    PointList()\x = 450
    PointList()\y = 30
    
    AddElement(PointList())
    PointList()\x = 460
    PointList()\y = 32
    
    AddElement(PointList())
    PointList()\x = 520
    PointList()\y = 45
    
    AddElement(PointList())
    PointList()\x = 570
    PointList()\y = 55
    
  CompilerEndIf
  
  If LinearCurveFinder(PointList(), @a, @b)
    ;Debug "y = " + StrD(a, 6) + "x + " + StrD(b, 6)
    
    
    ForEach PointList()
      If MaxX < PointList()\x : MaxX = PointList()\x : EndIf
      If MaxY < PointList()\y : MaxY = PointList()\y : EndIf
      If MinX > PointList()\x : MinX = PointList()\x : EndIf
      If MinY > PointList()\y : MinY = PointList()\y : EndIf
    Next
    ;Debug StrD(MinX, 0) + " < X < " + StrD(MaxX, 0)
    ;Debug StrD(MinY, 0) + " < Y < " + StrD(MaxY, 0)
    
    Equation$ = "y = " + StrD(a, 6) + " x "
    If b >= 0
      Equation$ + "+"
    Else
      Equation$ + "-"
    EndIf
    Equation$ + " " + StrD(Abs(b), 6)
    
    OpenWindow(0, 0, 0, #Width, #Height, "Line Curve Finder:    " + Equation$, #PB_Window_MinimizeGadget|#PB_Window_ScreenCentered)
    CanvasGadget(0, 0, 0, #Width, #Height)
    
    XFactor = MaxX / #Width
    If XFactor > 1 : FactorIsLargerThanOne = #True : EndIf
    YFactor = MaxY / #Height
    If XFactor > 1 : FactorIsLargerThanOne = #True : EndIf
    
    If FactorIsLargerThanOne
      If XFactor > YFactor
        Factor = XFactor
      Else
        Factor = YFactor
      EndIf
    Else
      If XFactor > YFactor
        Factor = 1 / XFactor
      Else
        Factor = 1 / YFactor
      EndIf
    EndIf
    
    ;Debug "Factor: " + StrD(Factor)
    
    If StartVectorDrawing(CanvasVectorOutput(0))
      FlipCoordinatesY(#Height / 2)
      ForEach PointList()
        AddPathCircle(PointList()\x * Factor, PointList()\y * Factor, 4)
        ClosePath()
      Next
      VectorSourceColor(RGBA(255, 0, 0, 255))
      StrokePath(10)
      
      MovePathCursor(0, b * Factor)
      AddPathLine(MaxX * Factor, (a * MaxX + b) * Factor)
      ClosePath()
      VectorSourceColor(RGBA(0, 0, 255, 255))
      StrokePath(4)
      
      StopVectorDrawing()
    EndIf
    
    Repeat
      Event = WaitWindowEvent()
    Until Event = #PB_Event_CloseWindow
    
  Else
    Debug "Was not able to build the Inverse matrix"
  EndIf
  
CompilerEndIf
infratec
Always Here
Always Here
Posts: 6817
Joined: Sun Sep 07, 2008 12:45 pm
Location: Germany

Re: A linear curve / equation finder

Post by infratec »

Added a graphical output.
At least working with the 2 examples.
Post Reply