SparkFun Forums 

Where electronics enthusiasts find answers.

Your source for all things Atmel.
By josheeg
#156341
Linear Discriminent Analisis Arduino Duo code... pattern matching...
Can someone test it with a duo or run the pc example agents larger sets of data or confurm the numbers are slightly different than the excel spread sheet and I am using doubles does excel use floats?

Binary sketch size: 183,028 bytes (of a 524,288 byte maximum) - 34% used

So I know it fits...

Eigen library is needed and there are arduino tutorials for the arduino duo you need to select that as the board type.
Eigen website shows how to setup the pc enviroment I am using code blocks ide.
Code: Select all
//Joshua Wojnas
//lda example reference http://people.revoledu.com/kardi/tutorial/LDA/Numerical%20Example.html
//reference http://eigen.tuxfamily.org/dox/TutorialMatrixArithmetic.html
// Example By: RandomVibe
// Eigen Doc: http://eigen.tuxfamily.org/dox/
// Quick Reference: http://eigen.tuxfamily.org/dox/QuickRefPage.html

#include <iostream>
#include <math.h>
#include <Eigen312.h>     // Calls main Eigen3.1.2 matrix class library
#include <Dense>
#include <LU>             // Calls inverse, determinant, LU decomp., etc.
using namespace Eigen;    // Eigen related statement; simplifies syntax for declaration of matrices
using Eigen::MatrixXd;
using namespace std;

void print_mtxf(const Eigen::MatrixXf& K);


void setup() {

    Serial.begin(9600);
   
    // DECLARE MATRICES
    //--------------------
  //data matrix x row col
MatrixXd x(4,2);
//cur
x(0,0) = 2.95;//g0
x(1,0) = 2.53;
x(2,0) = 3.57;
x(3,0) = 3.16;

//dia
x(0,1) = 6.63;//g0
x(1,1) = 7.79;
x(2,1) = 5.65;
x(3,1) = 5.47;

MatrixXd x1(3,2);

x1(0,0) = 2.58;//g1
x1(1,0) = 2.16;
x1(2,0) = 3.27;

x1(0,1) = 4.46;//g1
x1(1,1) = 6.22;
x1(2,1) = 3.52;

MatrixXd ui(1,2);//group/feature

ui(0,0)=(x(0,0) + x(1,0) + x(2,0) + x(3,0))/4.0;
ui(0,1)=(x(0,1) + x(1,1) + x(2,1) + x(3,1))/4.0;

//transpose ui to uit
MatrixXd uit(1,2);//group/feature
uit=ui.transpose();


MatrixXd ui1(1,2);//group/feature

ui1(0,0)=(x1(0,0) + x1(1,0) + x1(2,0))/3.0;
ui1(0,1)=(x1(0,1) + x1(1,1) + x1(2,1))/3.0;

MatrixXd ui1t(1,2);//group/feature
ui1t=ui1.transpose();

MatrixXd u(1,2);//all group/feature
u(0,0)=(x(0,0) + x(1,0) + x(2,0) + x(3,0) + x1(0,0) + x1(1,0) + x1(2,0))/7.0;

u(0,1)=(x(0,1) + x(1,1) + x(2,1) + x(3,1) + x1(0,1) + x1(1,1) + x1(2,1))/7.0;

MatrixXd ximinu(4,2);
//cur
ximinu(0,0) = x(0,0) - u(0,0);//f0
ximinu(1,0) = x(1,0) - u(0,0);
ximinu(2,0) = x(2,0) - u(0,0);
ximinu(3,0) = x(3,0) - u(0,0);
//dia
ximinu(0,1) = x(0,1) - u(0,1);//f1
ximinu(1,1) = x(1,1) - u(0,1);
ximinu(2,1) = x(2,1) - u(0,1);
ximinu(3,1) = x(3,1) - u(0,1);
std::cout << ximinu << std::endl;
cout << "" << endl;

MatrixXd ximinu1(3,2);
//cur
ximinu1(0,0) = x1(0,0) - u(0,0);
ximinu1(1,0) = x1(1,0) - u(0,0);
ximinu1(2,0) = x1(2,0) - u(0,0);
//dia
ximinu1(0,1) = x1(0,1) - u(0,1);//g1
ximinu1(1,1) = x1(1,1) - u(0,1);
ximinu1(2,1) = x1(2,1) - u(0,1);

MatrixXd ximinut(4,2);
ximinut= ximinu.transpose();

MatrixXd ximinu1t(3,2);
ximinu1t= ximinu1.transpose();

MatrixXd ci(2,2);
ci = ( ximinut * ximinu ) /4.0;

MatrixXd ci1(2,2);
ci1 = ( ximinu1t * ximinu1 ) /3.0;

MatrixXd c(2,2);
c(0,0) = 4.0/7.0 * ci(0,0) + 3.0/7.0 * ci1(0,0);
c(1,0) = 4.0/7.0 * ci(1,0) + 3.0/7.0 * ci1(1,0);
c(0,1) = 4.0/7.0 * ci(0,1) + 3.0/7.0 * ci1(0,1);
c(1,1) = 4.0/7.0 * ci(1,1) + 3.0/7.0 * ci1(1,1);

MatrixXd cinverse(2,2);
cinverse=c.inverse();
 
MatrixXd xk(1,2);//new data
xk(0,0) = 2.81;//f1
xk(0,1) = 5.46;

MatrixXd xkt(1,2);//new data
xkt=xk.transpose();

MatrixXd lnp1(1,1);//p1 4/7
lnp1(0,0) = 0.0;
lnp1(0,0)=log(4.0/7.0);

MatrixXd lnp2(1,1);//p1 3/7
lnp2(0,0) = 0.0;
lnp2(0,0)=log(4.0/7.0);

MatrixXd f1(1,1);//ui   * cinverse * xkt  -1.0/2.0 * ui  * cinverse * uit + lnp1
MatrixXd f2(1,1);//ui1  * cinverse * xkt -1.0/2.0 * ui1 * cinverse * ui1t + lnp2

f1 = ui   * cinverse * xkt  -1.0/2.0 * ui  * cinverse * uit + lnp1;
f2 = ui1  * cinverse * xkt -1.0/2.0 * ui1 * cinverse * ui1t + lnp2;
 
if (f1(0,0)>f2(0,0)) cout << "the data is in group 1" << endl;
if (f1(0,0)<f2(0,0)) cout << "the data is in group 2" << endl;
    // Print Result
    //----------------------------
     //print_mtxf(K);      // Print Matrix Result (passed by reference)
   
}




void loop() {
  // put your main code here, to run repeatedly:
 
}




// PRINT MATRIX (float type)
//-----------------------------
void print_mtxf(const Eigen::MatrixXf& X) 
{
    int i, j, nrow, ncol;
   
    nrow = X.rows();
    ncol = X.cols();

    Serial.print("nrow: "); Serial.println(nrow);
    Serial.print("ncol: "); Serial.println(ncol);       
    Serial.println();
   
    for (i=0; i<nrow; i++)
    {
        for (j=0; j<ncol; j++)
        {
            Serial.print(X(i,j), 6);   // print 6 decimal places
            Serial.print(", ");
        }
        Serial.println();
    }
    Serial.println();
}
Here is the pc version...
Code: Select all
//Josh W
//lda example reference http://people.revoledu.com/kardi/tutorial/LDA/Numerical%20Example.html
//reference http://eigen.tuxfamily.org/dox/TutorialMatrixArithmetic.html

#include <iostream>
#include <math.h>
#include <Eigen/Dense>
#include <Eigen/LU>
using Eigen::MatrixXd;
using namespace std;
using namespace Eigen;
int main()
{
cout << "Hello Pattern matching Linear Discriminent Analisis!" << endl;
cout << "" << endl;
cout << "x data" << endl;

//data matrix x row col
MatrixXd x(4,2);
//cur
x(0,0) = 2.95;//g0
x(1,0) = 2.53;
x(2,0) = 3.57;
x(3,0) = 3.16;

//dia
x(0,1) = 6.63;//g0
x(1,1) = 7.79;
x(2,1) = 5.65;
x(3,1) = 5.47;

std::cout << x << std::endl;
cout << "End of x data" << endl;
cout << "" << endl;

cout << "x1 data" << endl;

MatrixXd x1(3,2);

x1(0,0) = 2.58;//g1
x1(1,0) = 2.16;
x1(2,0) = 3.27;

x1(0,1) = 4.46;//g1
x1(1,1) = 6.22;
x1(2,1) = 3.52;

std::cout << x1 << std::endl;
cout << "End of x1 data" << endl;
cout << "" << endl;

//group adv
cout << "x data adverage ui" << endl;

MatrixXd ui(1,2);//group/feature

ui(0,0)=(x(0,0) + x(1,0) + x(2,0) + x(3,0))/4.0;
ui(0,1)=(x(0,1) + x(1,1) + x(2,1) + x(3,1))/4.0;
std::cout << ui << std::endl;
cout << "" << endl;

cout << "transpose ui to uit" << endl;
//transpose ui to uit
MatrixXd uit(1,2);//group/feature
uit=ui.transpose();
std::cout << uit << std::endl;
cout << "" << endl;

cout << "x1 data adverage ui1" << endl;

MatrixXd ui1(1,2);//group/feature

ui1(0,0)=(x1(0,0) + x1(1,0) + x1(2,0))/3.0;
ui1(0,1)=(x1(0,1) + x1(1,1) + x1(2,1))/3.0;
std::cout << ui1 << std::endl;
cout << "" << endl;

//transpose ui to uit
std::cout << "transpose ui1 to ui1t" << std::endl;
MatrixXd ui1t(1,2);//group/feature
ui1t=ui1.transpose();
std::cout << ui1t << std::endl;
cout << "" << endl;


cout << "x & x1 data adverage u" << endl;
MatrixXd u(1,2);//all group/feature
u(0,0)=(x(0,0) + x(1,0) + x(2,0) + x(3,0) + x1(0,0) + x1(1,0) + x1(2,0))/7.0;

u(0,1)=(x(0,1) + x(1,1) + x(2,1) + x(3,1) + x1(0,1) + x1(1,1) + x1(2,1))/7.0;

std::cout << u << std::endl;
cout << "" << endl;


cout << "mean corrected data xig - u" << endl;
MatrixXd ximinu(4,2);
//cur
ximinu(0,0) = x(0,0) - u(0,0);//f0
ximinu(1,0) = x(1,0) - u(0,0);
ximinu(2,0) = x(2,0) - u(0,0);
ximinu(3,0) = x(3,0) - u(0,0);
//dia
ximinu(0,1) = x(0,1) - u(0,1);//f1
ximinu(1,1) = x(1,1) - u(0,1);
ximinu(2,1) = x(2,1) - u(0,1);
ximinu(3,1) = x(3,1) - u(0,1);
std::cout << ximinu << std::endl;
cout << "" << endl;

cout << "mean corrected data xi1 - u" << endl;


MatrixXd ximinu1(3,2);
//cur
ximinu1(0,0) = x1(0,0) - u(0,0);
ximinu1(1,0) = x1(1,0) - u(0,0);
ximinu1(2,0) = x1(2,0) - u(0,0);
//dia
ximinu1(0,1) = x1(0,1) - u(0,1);//g1
ximinu1(1,1) = x1(1,1) - u(0,1);
ximinu1(2,1) = x1(2,1) - u(0,1);

std::cout << ximinu1 << std::endl;

cout << " " << endl;

cout << "Transpose matricies" << endl;
cout << "xi - u T" << endl;

MatrixXd ximinut(4,2);
ximinut= ximinu.transpose();
std::cout << ximinut<< std::endl;
cout << " " << endl;

cout << "xi1 - u T" << endl;
MatrixXd ximinu1t(3,2);
ximinu1t= ximinu1.transpose();
std::cout << ximinu1t << std::endl;
cout << " " << endl;

cout << "Covariance matrix of group ci" << endl;
MatrixXd ci(2,2);
ci = ( ximinut * ximinu ) /4.0;
std::cout << ci << std::endl;
cout << "" << endl;

cout << "Covariance matrix of group ci1" << endl;
MatrixXd ci1(2,2);
ci1 = ( ximinu1t * ximinu1 ) /3.0;
std::cout << ci1 << std::endl;
cout << "" << endl;

cout << "Pooled within group Covariance matrix c" << endl;
MatrixXd c(2,2);
c(0,0) = 4.0/7.0 * ci(0,0) + 3.0/7.0 * ci1(0,0);
c(1,0) = 4.0/7.0 * ci(1,0) + 3.0/7.0 * ci1(1,0);
c(0,1) = 4.0/7.0 * ci(0,1) + 3.0/7.0 * ci1(0,1);
c(1,1) = 4.0/7.0 * ci(1,1) + 3.0/7.0 * ci1(1,1);
std::cout << c << std::endl;
cout << "" << endl;

cout << "inverse of Pooled within group Covariance matrix cinverse" << endl;
MatrixXd cinverse(2,2);
 cinverse=c.inverse();
std::cout << cinverse << std::endl;
cout << "" << endl;

cout << "Probability of a group" << endl;
cout << "x = 4/7 x1 = 3/7 p(4/7,3/7)" << endl;
Vector2d p(4.0/7.0,3.0/7.0);
std::cout << p << std::endl;
cout << "" << endl;

//discriminent function
//MatrixXd fi(2,2);
//double xk = 1;
//formula for calculatin likelyhood of data in a group...
//fi = uig  * cinverse xkt - 1/2uig cinverse uitg + ln probability
//fi  = ui  * cinverse * x.transpose() - 1/2 * ui * cinverse;// * uit;// +.3;
//fi1 = ui1 * cinverse * x1.transpose() - 1/2 * ui1 * cinverse ui1.transpose + ln(3/7);
//std::cout << fi << std::endl;
MatrixXd xk(1,2);//new data
xk(0,0) = 2.81;//f1
xk(0,1) = 5.46;

MatrixXd xkt(1,2);//new data
xkt=xk.transpose();
//ui  * cinverse * xkt
//1.0/2.0 * ui * cinverse * uit
cout << "ln(p1) =" << endl;

MatrixXd lnp1(1,1);//p1 4/7
lnp1(0,0) = 0.0;
lnp1(0,0)=log(4.0/7.0);

MatrixXd lnp2(1,1);//p1 3/7
lnp2(0,0) = 0.0;
lnp2(0,0)=log(4.0/7.0);

std::cout << lnp1 << std::endl;
cout << "" << endl;
cout << "The LDA of two features of data f1 f2" << endl;

//std::cout << ui   * cinverse * xkt  -1.0/2.0 * ui  * cinverse * uit + lnp1  << std::endl;
//std::cout << ui1  * cinverse * xkt -1.0/2.0 * ui1 * cinverse * ui1t + lnp2 << std::endl;
MatrixXd f1(1,1);//ui   * cinverse * xkt  -1.0/2.0 * ui  * cinverse * uit + lnp1
MatrixXd f2(1,1);//ui1  * cinverse * xkt -1.0/2.0 * ui1 * cinverse * ui1t + lnp2

f1 = ui   * cinverse * xkt  -1.0/2.0 * ui  * cinverse * uit + lnp1;
f2 = ui1  * cinverse * xkt -1.0/2.0 * ui1 * cinverse * ui1t + lnp2;
std::cout << f1  << std::endl;
std::cout << f2  << std::endl;

cout << "The larger of the discriminent functions results is what group the data is in" << endl;
if (f1(0,0)>f2(0,0)) cout << "the data is in group 1" << endl;
if (f1(0,0)<f2(0,0)) cout << "the data is in group 2" << endl;

cout << "" << endl;


cout << "End of program!" << endl;

    return 0;
}