Xmipp
v3.23.11-Nereus
xmipp
libraries
dimred
gplvm.cpp
Go to the documentation of this file.
1
/***************************************************************************
2
*
3
* Authors: Francisco Sanz Encinas franciscosanz89@gmail.com (2013)
4
*
5
* Unidad de Bioinformatica of Centro Nacional de Biotecnologia , CSIC
6
*
7
* This program is free software; you can redistribute it and/or modify
8
* it under the terms of the GNU General Public License as published by
9
* the Free Software Foundation; either version 2 of the License, or
10
* (at your option) any later version.
11
*
12
* This program is distributed in the hope that it will be useful,
13
* but WITHOUT ANY WARRANTY; without even the implied warranty of
14
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15
* GNU General Public License for more details.
16
*
17
* You should have received a copy of the GNU General Public License
18
* along with this program; if not, write to the Free Software
19
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
20
* 02111-1307 USA
21
*
22
* All comments concerning this program package may be sent to the
23
* e-mail address 'xmipp@cnb.csic.es'
24
***************************************************************************/
25
26
#include "
gplvm.h
"
27
#include <
data/numerical_tools.h
>
28
29
void
GPLVM::setSpecificParameters
(
double
sigma)
30
{
31
this->sigma=
sigma
;
32
}
33
34
double
GPLVM::objectiveFunction
()
35
{
36
sumY2
.
initZeros
(
MAT_YSIZE
(
Y
));
37
FOR_ALL_ELEMENTS_IN_MATRIX2D
(
Y
)
38
VEC_ELEM
(
sumY2
,
i
)+=
MAT_ELEM
(
Y
,
i
,
j
)*
MAT_ELEM
(
Y
,
i
,
j
);
39
40
double
aux=1.0/(2*
sigma
*
sigma
);
41
K
.
resizeNoCopy
(
MAT_YSIZE
(
Y
),
MAT_YSIZE
(
Y
));
42
FOR_ALL_ELEMENTS_IN_MATRIX2D
(
K
)
43
{
44
double
aux2=0;
45
for
(
int
k
=0;
k
<
MAT_XSIZE
(
Y
);
k
++)
46
aux2+=
MAT_ELEM
(
Y
,
i
,
k
)*
MAT_ELEM
(
Y
,
j
,
k
);
47
MAT_ELEM
(
K
,
i
,
j
)=exp((2*aux2-
VEC_ELEM
(
sumY2
,
i
)-
VEC_ELEM
(
sumY2
,
j
))*aux);
48
}
49
50
matrixOperation_AAt
(*
X
,
tmp
);
51
tmp
=
K
.
inv
() *
tmp
;
52
53
double
d
=
MAT_YSIZE
(*
X
);
54
double
n
=
MAT_XSIZE
(*
X
);
55
return
(-(d*n/2) *
log
( 2 * M_PI) - (n/2) *
log
(
K
.
det
()+2.2e-308) - 0.5 *
tmp
.
trace
());
56
}
57
58
double
gplvmObjectiveFuntion
(
double
*p,
void
*
prm
)
59
{
60
auto
*gpvlm=(
GPLVM
*)prm;
61
Matrix2D<double>
&
Y
=gpvlm->Y;
62
memcpy(&
MAT_ELEM
(Y,0,0),&(p[1]),
MAT_XSIZE
(Y)*
MAT_YSIZE
(Y)*
sizeof
(
double
));
63
double
c
=gpvlm->objectiveFunction();
64
return
c
;
65
}
66
67
void
GPLVM::reduceDimensionality
()
68
{
69
// Compute the first guess
70
PCA::reduceDimensionality
();
71
72
// Refine
73
Matrix1D<double>
pY(
MAT_XSIZE
(
Y
)*
MAT_YSIZE
(
Y
)),
steps
(
MAT_XSIZE
(
Y
)*
MAT_YSIZE
(
Y
));
74
steps.
initConstant
(1);
75
memcpy(&
VEC_ELEM
(pY,0),&
MAT_ELEM
(
Y
,0,0),
VEC_XSIZE
(pY)*
sizeof
(
double
));
76
double
goal;
77
int
iter
;
78
powellOptimizer
(pY,1,
VEC_XSIZE
(pY),&
gplvmObjectiveFuntion
,
this
,0.01,goal,iter,steps,
true
);
79
memcpy(&
MAT_ELEM
(
Y
,0,0),&
VEC_ELEM
(pY,0),
VEC_XSIZE
(pY)*
sizeof
(
double
));
80
}
GPLVM
Definition:
gplvm.h:39
FOR_ALL_ELEMENTS_IN_MATRIX2D
#define FOR_ALL_ELEMENTS_IN_MATRIX2D(m)
Definition:
matrix2d.h:104
gplvm.h
MAT_YSIZE
#define MAT_YSIZE(m)
Definition:
matrix2d.h:124
numerical_tools.h
VEC_ELEM
#define VEC_ELEM(v, i)
Definition:
matrix1d.h:245
DimRedAlgorithm::Y
Matrix2D< double > Y
Output data.
Definition:
dimred_tools.h:147
VEC_XSIZE
#define VEC_XSIZE(m)
Definition:
matrix1d.h:77
c
doublereal * c
Definition:
numerical_recipes.cpp:2230
DimRedAlgorithm::X
Matrix2D< double > * X
Pointer to input data.
Definition:
dimred_tools.h:141
Matrix2D::trace
T trace() const
Definition:
matrix2d.cpp:1304
Matrix2D::inv
void inv(Matrix2D< T > &result) const
Definition:
matrix2d.cpp:663
powellOptimizer
void powellOptimizer(Matrix1D< double > &p, int i0, int n, double(*f)(double *x, void *), void *prm, double ftol, double &fret, int &iter, const Matrix1D< double > &steps, bool show)
Definition:
numerical_tools.cpp:44
iter
glob_prnt iter
Definition:
numerical_recipes.cpp:4698
GPLVM::reduceDimensionality
void reduceDimensionality()
Reduce dimensionality.
Definition:
gplvm.cpp:67
i
#define i
Definition:
numerical_recipes.cpp:2493
k
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
d
doublereal * d
Definition:
numerical_recipes.cpp:2230
Matrix2D::resizeNoCopy
void resizeNoCopy(int Ydim, int Xdim)
Definition:
matrix2d.h:534
matrixOperation_AAt
void matrixOperation_AAt(const Matrix2D< double > &A, Matrix2D< double > &C)
Definition:
matrix2d.cpp:449
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition:
matrix2d.h:116
Matrix2D< double >
GPLVM::setSpecificParameters
void setSpecificParameters(double sigma=1.0)
Set specific parameters.
Definition:
gplvm.cpp:29
log
void log(Image< double > &op)
Definition:
image_operate.cpp:228
GPLVM::sigma
double sigma
Definition:
gplvm.h:42
PCA::reduceDimensionality
void reduceDimensionality()
Reduce dimensionality.
Definition:
pca.cpp:28
GPLVM::tmp
Matrix2D< double > tmp
Definition:
gplvm.h:53
Matrix2D::det
T det() const
Definition:
matrix2d.cpp:38
Matrix1D::initZeros
void initZeros()
Definition:
matrix1d.h:592
j
#define j
Definition:
numerical_recipes.cpp:2493
steps
double steps
Definition:
numerical_recipes.cpp:5166
GPLVM::objectiveFunction
double objectiveFunction()
Objective function.
Definition:
gplvm.cpp:34
GPLVM::K
Matrix2D< double > K
Definition:
gplvm.h:53
Matrix1D< double >
MAT_XSIZE
#define MAT_XSIZE(m)
Definition:
matrix2d.h:120
prm
ProgClassifyCL2D * prm
Definition:
mpi_classify_CL2D.cpp:35
GPLVM::sumY2
Matrix1D< double > sumY2
Definition:
gplvm.h:54
Matrix1D::initConstant
void initConstant(T val)
Definition:
matrix1d.cpp:83
n
int * n
Definition:
numerical_recipes.cpp:2229
gplvmObjectiveFuntion
double gplvmObjectiveFuntion(double *p, void *prm)
Definition:
gplvm.cpp:58
Generated by
1.8.13