Dynamic_Imitation/f_estimate.cpp at master · niladridas/Dynamic_Imitation · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
/*
 * f_estimate.cpp
 *
 *  Created on: 29-Nov-2014
 *      Author: niladriisl
 */

#include </usr/include/eigen3/Eigen/LU>
#include </usr/include/eigen3/Eigen/SVD>
#include </usr/include/eigen3/Eigen/Core>
#include <cmath>
#include <boost/math/constants/constants.hpp>
#include <iostream>
#include "gaussPdf.h"
#include <math.h>

using namespace Eigen;

long double* f_estimate(const VectorXf input, const MatrixXf Mean,
		const MatrixXf Sigma, const VectorXf Priors, const VectorXf Xi_star) {

	// First calculate dimension of priors
	int dim_priors = Priors.rows(); // priors comes as a rows
	int dim_input = input.rows();

	MatrixXf input_mean_matrix(dim_input, dim_priors); //   input  position(xi) mean vector for 6 gaussians
	MatrixXf out_mean_matrix(dim_input, dim_priors); //
	MatrixXf Sigma_input(dim_input * dim_priors, dim_input); //  extracted  sigma for input
	MatrixXf Sigma_output_input(dim_input * dim_priors, dim_input); //   extracted sigma sigmadot for input
	MatrixXf Inverse_Sigma_input(dim_input * dim_priors, dim_input); // inverse of the sigma matrix taking 6, 2x2 matrix
	VectorXf p(dim_priors); //probability of input for 6 components
	VectorXf h(dim_priors); // h value in the linear regression equation
	MatrixXf A_matrix(dim_input * dim_priors, dim_input); // 'A' matrix in the linear regression equation
	MatrixXf B_matrix(dim_input, dim_priors); // 'B' matrix in the linear regression equation
	VectorXf output(dim_input); // output 'xi dot' of regression equation

	int i = 0;
	int k = 0;

	input_mean_matrix = Mean.block(0, 0, dim_input, dim_priors);
	out_mean_matrix = Mean.block(dim_input, 0, dim_input, dim_priors);

	//

	for (i = 0; i < dim_priors; i++) {
		Sigma_input.block(dim_input * i, 0, dim_input, dim_input) = Sigma.block(
				2 * dim_input * i, 0, dim_input, dim_input);
	}

	//
	for (i = 0; i < dim_priors; i++) {
		Sigma_output_input.block(dim_input * i, 0, dim_input, dim_input) =
				Sigma.block(2 * dim_input * i + dim_input, 0, dim_input,
						dim_input);
	}

	//

	for (i = 0; i < dim_input * dim_priors - dim_input; i = i + dim_input) //inverse
			{
		Inverse_Sigma_input.block(i, 0, dim_input, dim_input) =
				Sigma_input.block(i, 0, dim_input, dim_input).inverse();
	}

	//

	long double prob_input = 0;

	for (k = 0; k < dim_priors; k++) {
		double prob_input_tmp = gausspdf(input, input_mean_matrix.col(k),
				Sigma_input.block(dim_input * k, 0, dim_input, dim_input));
		double tmp_Pk = Priors(k) * prob_input_tmp;
		p(k) = tmp_Pk;
		prob_input = prob_input + tmp_Pk;

	}

	for (k = 0; k < dim_priors; k++) {
		double h_k = p(k) / prob_input;
		h(k) = h_k;
	}

	for (k = 0; k < dim_input * dim_priors - dim_input; k = k + dim_input) {
		A_matrix.block(k, 0, dim_input, dim_input) = Sigma_output_input.block(k,
				0, dim_input, dim_input)
				* Inverse_Sigma_input.block(k, 0, dim_input, dim_input);
	}

	for (k = 0; k < dim_priors; k++) {
		B_matrix.col(k) = out_mean_matrix.col(k)
				- A_matrix.block(dim_input * k, 0, dim_input, dim_input)
						* input_mean_matrix.col(k);
	}

	output.Zero(dim_input, 1);

	for (k = 0; k < dim_priors; k++) {
		output = output
				+ h(k)
						* (A_matrix.block(dim_input * k, 0, dim_input,
								dim_input) * input + B_matrix.col(k));

	}

	long double *output_array;
	output_array = new long double[dim_input];

	for (i = 0; i < dim_input; i++) {
		output_array[i] = output(i);
	}

	return output_array;
}