Algorithme de Levenberg Marquardt

**TopCao** · 06/01/2010, 22h26

Bonjour à tous,

Je tiens tout d'abord à préciser que je suis débutant et totalement autodidact en C++.

Alors je dispose d'une série de données (48 valeurs) dont je sais qu'elles correspondent à une fonction du type :

S(b) = Sa.exp(-b.ADCa) + Sb.exp(-b.ADCb) + Sc.exp(-b.ADCc) + bruit

S étant le signal et b ma variable.

Je voudrais donc à l'aide d'un algorithme de Levenberg Marquardt retrouver les coefficients Sa, Sb, Sc, ADCa, ADCb et ADCc.

Je voulais mettre mes valeurs dans un tableau de pointeurs et générer mes 6 coeff mais je n'y arrive pas.

Je suis donc preneur de toute aide et je peux envoyer les bouts de code que j'ai commencé à taper.

Précisions : je prog en c++ avec Xcode (sous Mac) pour faire des progs en ligne de commande.

Merci à vous

alex

**Klaim** · 07/01/2010, 02h13

Il faudrait que tu montres le code que tu as déjà pour que tout ça soit bien clair stp. Comme ça on peut te guider.

**TopCao** · 07/01/2010, 18h17

Bonjour,

Voici mon code
3 fichiers :

itkCompex.h : mon header
itkComputeCompex.cxx : mon prog
itkCompex.txx : le template qui contient l'algo

Pour bien comprendre :
J'écris ce programme dans le but de l'exécuter sur une série d'images (d'ou la structure ITK rencontrée partout) mais je souhaite tester la fiabilité de l'algo de Levenberg sur une série de chiffres (en fait la valeur d'un voxel au cours du temps)

J'espère que vous n'allez pas être pris de convulsions en regardant mon code qui je l'admets est vraiment cracra et pas compilable en l'état !!

Merci pour vos conseils

Alex

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
/*
 *  itkCompex.h
 *  compex
 *
 *  Created by Alexandre Ladoux on 07/01/10.
 *  Copyright 2010 __LIMEC__. All rights reserved.
 *
 */
 
 
#ifndef _itk_Compex_h_
#define _itk_Compex_h_
 
#include <itkImageToImageFilter.h>
 
namespace itk
{
 
	template <class TInputImage, class TOutputImage>
    class ITK_EXPORT Compex :
	public ImageToImageFilter<TInputImage, TOutputImage>
	{
 
	public:
 
		typedef Compex			Self;
		typedef ImageToImageFilter<TInputImage, TOutputImage> Superclass;
		typedef SmartPointer<Self> Pointer;
		typedef SmartPointer<const Self> ConstPointer;
 
		itkNewMacro  (Self);
		itkTypeMacro (Compex, ImageToImagFilter);
 
		typedef TInputImage  InputImageType;
		typedef TOutputImage OutputImageType;
 
		typedef typename InputImageType::PixelType   InputPixelType;
		typedef typename OutputImageType::PixelType  OutputPixelType;
		typedef typename OutputImageType::RegionType OutputImageRegionType;
		typedef OutputPixelType                      ScalarType;
 
		// Dé-commenter lors de l'application du process à des images
		// itkSetObjectMacro (DataSetImages, InputImageType);
		// itkGetObjectMacro (DataSetImages, InputImageType);
 
 
		itkSetMacro (v01, double);
		itkGetMacro (v01, double);
		itkSetMacro (v02, double);
		itkGetMacro (v02, double);
		itkSetMacro (v03, double);
		itkGetMacro (v03, double);
		itkSetMacro (v04, double);
		itkGetMacro (v04, double);
		itkSetMacro (v05, double);
		itkGetMacro (v05, double);
		itkSetMacro (v06, double);
		itkGetMacro (v06, double);
		itkSetMacro (v07, double);
		itkGetMacro (v07, double);
		itkSetMacro (v08, double);
		itkGetMacro (v08, double);
		itkSetMacro (v09, double);
		itkGetMacro (v09, double);
		itkSetMacro (v10, double);
		itkGetMacro (v10, double);
		itkSetMacro (v11, double);
		itkGetMacro (v11, double);
		itkSetMacro (v12, double);
		itkGetMacro (v12, double);
		itkSetMacro (v13, double);
		itkGetMacro (v13, double);
		itkSetMacro (v14, double);
		itkGetMacro (v14, double);		
		itkSetMacro (v15, double);
		itkGetMacro (v15, double);
		itkSetMacro (v16, double);
		itkGetMacro (v16, double);
		itkSetMacro (v17, double);
		itkGetMacro (v17, double);
		itkSetMacro (v18, double);
		itkGetMacro (v18, double);
		itkSetMacro (v19, double);
		itkGetMacro (v19, double);
		itkSetMacro (v20, double);
		itkGetMacro (v20, double);
		itkSetMacro (v21, double);
		itkGetMacro (v21, double);
		itkSetMacro (v22, double);
		itkGetMacro (v22, double);
		itkSetMacro (v23, double);
		itkGetMacro (v23, double);
		itkSetMacro (v24, double);
		itkGetMacro (v24, double);
		itkSetMacro (v25, double);
		itkGetMacro (v25, double);
		itkSetMacro (v26, double);
		itkGetMacro (v26, double);
		itkSetMacro (v27, double);
		itkGetMacro (v27, double);
		itkSetMacro (v28, double);
		itkGetMacro (v28, double);				
		itkSetMacro (v29, double);
		itkGetMacro (v29, double);
		itkSetMacro (v30, double);
		itkGetMacro (v30, double);
		itkSetMacro (v31, double);
		itkGetMacro (v31, double);
		itkSetMacro (v32, double);
		itkGetMacro (v32, double);
		itkSetMacro (v33, double);
		itkGetMacro (v33, double);
		itkSetMacro (v34, double);
		itkGetMacro (v34, double);
		itkSetMacro (v35, double);
		itkGetMacro (v35, double);
		itkSetMacro (v36, double);
		itkGetMacro (v36, double);
		itkSetMacro (v37, double);
		itkGetMacro (v37, double);
		itkSetMacro (v38, double);
		itkGetMacro (v38, double);
		itkSetMacro (v39, double);
		itkGetMacro (v39, double);
		itkSetMacro (v40, double);
		itkGetMacro (v40, double);
		itkSetMacro (v41, double);
		itkGetMacro (v41, double);
		itkSetMacro (v42, double);
		itkGetMacro (v42, double);		
		itkSetMacro (v43, double);
		itkGetMacro (v43, double);
		itkSetMacro (v44, double);
		itkGetMacro (v44, double);
		itkSetMacro (v45, double);
		itkGetMacro (v45, double);
		itkSetMacro (v46, double);
		itkGetMacro (v46, double);
		itkSetMacro (v47, double);
		itkGetMacro (v47, double);
		itkSetMacro (v48, double);
		itkGetMacro (v48, double);
 
 
	protected:
		Compex();
		~Compex(){};
 
 
		void BeforeThreadedGenerateData(void);
		void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, int threadId);
		void PrintSelf(std::ostream& os, Indent indent) const
		{
			Superclass::PrintSelf(os,indent);
		}
 
 
 
	private:
		CompexFilter (const Self&);
		void operator=(const Self&);
 
		double                  m_v01;
		double					m_v02;
		double                  m_v03;
		double					m_v04;
		double                  m_v05;
		double					m_v06;
		double                  m_v07;
		double					m_v08;
		double                  m_v09;
		double					m_v10;
		double                  m_v11;
		double					m_v12;
		double                  m_v13;
		double					m_v14;
		double                  m_v15;
		double					m_v16;
		double                  m_v17;
		double					m_v18;
		double                  m_v19;
		double					m_v20;
		double                  m_v21;
		double					m_v22;
		double                  m_v23;
		double					m_v24;
		double                  m_v25;
		double					m_v26;
		double                  m_v27;
		double					m_v28;
		double                  m_v29;
		double					m_v30;
		double                  m_v31;
		double					m_v32;
		double                  m_v33;
		double					m_v34;
		double                  m_v35;
		double					m_v36;
		double                  m_v37;
		double					m_v38;
		double                  m_v39;
		double					m_v40;
		double                  m_v41;
		double					m_v42;
		double                  m_v43;
		double					m_v44;
		double                  m_v45;
		double					m_v46;
		double                  m_v47;
		double					m_v48;
 
 
		// Dé-commenter lors de l'application du process à des images
		//typename InputImageType::Pointer m_B0Image;
 
	};
 
} // end of namespace
 
#include "itkCompex.txx"
 
#endif

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
// ATTENTION : COMMENTAIRE POUR MOI MEME
//
// Copie du template itkCompex.txx créé uniquement pour avoir l'indentation xcode
// non compilé par Xcode
// Copier le corps du fichier dans itkCompex.txx
// après toute modification et avant toute compilation
 
 
#ifndef _itk_Compex_txx_
#define _itk_Compex_txx_
 
#include <itkCompex.h>
 
#include <itkImageRegionIterator.h>
#include <itkImageRegionConstIterator.h>
 
namespace itk
{
 
	template<class TInputImage, class TOutputImage>
	CompexFilter<TInputImage,TOutputImage>
	::CompexFilter()
	{
		m_v01  = 0.0;
		m_v02  = 0.0;
		m_v03  = 0.0;
		m_v04  = 0.0;
		m_v05  = 0.0;
		m_v06  = 0.0;
		m_v07  = 0.0;
		m_v08  = 0.0;
		m_v09  = 0.0;
		m_v10  = 0.0;
		m_v11  = 0.0;
		m_v12  = 0.0;
		m_v13  = 0.0;
		m_v14  = 0.0;
		m_v15  = 0.0;
		m_v16  = 0.0;
		m_v17  = 0.0;
		m_v18  = 0.0;
		m_v19  = 0.0;
		m_v20  = 0.0;
		m_v21  = 0.0;
		m_v22  = 0.0;
		m_v23  = 0.0;
		m_v24  = 0.0;
		m_v25  = 0.0;
		m_v26  = 0.0;
		m_v27  = 0.0;
		m_v28  = 0.0;
		m_v29  = 0.0;
		m_v30  = 0.0;
		m_v31  = 0.0;
		m_v32  = 0.0;
		m_v33  = 0.0;
		m_v34  = 0.0;
		m_v35  = 0.0;
		m_v36  = 0.0;
		m_v37  = 0.0;
		m_v38  = 0.0;
		m_v39  = 0.0;
		m_v40  = 0.0;
		m_v41  = 0.0;
		m_v42  = 0.0;
		m_v43  = 0.0;
		m_v44  = 0.0;
		m_v45  = 0.0;
		m_v46  = 0.0;
		m_v47  = 0.0;
		m_v48  = 0.0;
 
		// Dé-commenter lors de l'application du process à des images
		// m_B0Image = 0;
	}
 
 
	template<class TInputImage, class TOutputImage>
	void
	CompexFilter<TInputImage,TOutputImage>
	::BeforeThreadedGenerateData()
	{
		if( this->GetNumberOfInputs()==0 )
		{
			throw itk::ExceptionObject(__FILE__,__LINE__,"Error: Input is not set.");
		}
 
		if( m_B0Image.IsNull() )
		{
			throw itk::ExceptionObject(__FILE__,__LINE__,"Error: B0 image is not set.");
		}    
	}
 
 
	template<class TInputImage, class TOutputImage>
	void
	CompexFilter<TInputImage,TOutputImage>
	::ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, int threadId)
	{
 
		typedef ImageRegionIterator<OutputImageType>      IteratorOutputType;
		typedef ImageRegionConstIterator<InputImageType>  IteratorInputType;
 
 
		IteratorOutputType  itOut (this->GetOutput(),  outputRegionForThread);
		IteratorInputType   itIn  (this->GetInput(),   outputRegionForThread);
		//IteratorInputType   itB0  (this->GetB0Image(), outputRegionForThread);
 
 
		double v01 = this->GetBValue();
		double v02 = this->GetBValue();
		double v03 = this->GetBValue();
		double v04 = this->GetBValue();
		double v05 = this->GetBValue();
		double v06 = this->GetBValue();
		double v07 = this->GetBValue();
		double v08 = this->GetBValue();
		double v09 = this->GetBValue();
		double v10 = this->GetBValue();
		double v11 = this->GetBValue();
		double v12 = this->GetBValue();
		double v13 = this->GetBValue();
		double v14 = this->GetBValue();
		double v15 = this->GetBValue();
		double v16 = this->GetBValue();
		double v17 = this->GetBValue();
		double v18 = this->GetBValue();
		double v19 = this->GetBValue();
		double v20 = this->GetBValue();
		double v21 = this->GetBValue();
		double v22 = this->GetBValue();
		double v23 = this->GetBValue();
		double v24 = this->GetBValue();
		double v25 = this->GetBValue();
		double v26 = this->GetBValue();
		double v27 = this->GetBValue();
		double v28 = this->GetBValue();
		double v29 = this->GetBValue();
		double v30 = this->GetBValue();
		double v31 = this->GetBValue();
		double v32 = this->GetBValue();
		double v33 = this->GetBValue();
		double v34 = this->GetBValue();
		double v35 = this->GetBValue();
		double v36 = this->GetBValue();
		double v37 = this->GetBValue();
		double v38 = this->GetBValue();
		double v39 = this->GetBValue();
		double v40 = this->GetBValue();
		double v41 = this->GetBValue();
		double v42 = this->GetBValue();
		double v43 = this->GetBValue();
		double v44 = this->GetBValue();
		double v45 = this->GetBValue();
		double v46 = this->GetBValue();
		double v47 = this->GetBValue();
		double v48 = this->GetBValue();
 
 
 
// BOUCLE DE L'ALGORITHME DE LEVENBERG MARQUARDT		
 
 
		// Paramètres des données
		// int xRange = 128     ;	//Résolution en x quand image
		// int yRange = 128     ;	//Résolution en y quand image
		// int nSlices = 8      ;	//Résolution en z quand imageS
		int nbbValues = 48   ;  //Résolution en z
		int tableau[bValues] ;  //Tableau contenant les valeurs de b pour un voxel donné
 
 
		// Paramètres d'initialisation
		double Sa0  = 125	  ;
		double Sb0  = 99,001  ;
		double Sc0  = 77,28	  ;
		double ADCa = 0,0001  ;
		double ADCb = 0,00002 ;
		double ADCc = 0,00003 ;
 
		// Paramètres estimés
		double estSa0  = Sa0  ;
		double estSb0  = Sb0  ;
		double estSc0  = Sc0  ;
		double estADCa = ADCa ;
		double estADCb = ADCb ;
		double estADCc = ADCc ;
 
		// Fonction d'initialisation
		double estSignal = Sa0 * exp(-ADCa * bValue) + Sb0 * exp(-ADCb * bValue) +Sc0 * exp(-ADCc * bValue);
		int nbIter = 100;				// Nombre d'irérations
		double lambda = 0,001;          // Valeur initiale du facteur d'amortissement
		int nbParams = 6;				// Nombre de paramètres à estimer
 
 
		// Initialisation
		for (iter = 1; iter = nbIter; iter++) 
		{
 
			{
 
				tableau[0] = v01 ;				// tableau contenant à la suite les n valeurs de b
				...
				tableau[47] = v48 ;
 
 
				// ????? Evaluer la matrice jacobienne pour les paramètres estimés.
				double J = ???;
 
				// Faire l'erreur de distance pour les paramètres estimés
				estSignal = Sa0 * exp(-ADCa * bValue) + Sb0 * exp(-ADCb * bValue) +Sc0 * exp(-ADCc * bValue);
 
				// Calculer la matrice d'Hessian approximée avec J' la trasposée de J
				H = J' * J ;
 
				// appliquer le factor d'amortissement à la matrice hessienne
				double matriceLM = H + (lambda * nbParams);
 
				//incrémenter tous les paramètres
				double inc = - inv(matriceLM) * (J'*k);
 
 
 
			}		
 
		}
 
 
#endif

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
#include "GetPot.h"
 
#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkCompex.h>
 
//#include "itkGISImageIOFactory.h"
 
void PrintHelp (const char* exec)
{
	std::cout << "Usage: " << std::endl;
	std::cout << exec << "	<-v01 value01> 
	<-v02 value02> 
	<-v03 value03> 
	<-v04 value04> 
	<-v05 value06> 
	<-v06 value06> 
	<-v07 value07>
	<-v08 value08> 
	<-v09 value09> 
	<-v10 value10> 
	<-v11 value11> 
	<-v12 value12> 
	<-v13 value13>
	<-v14 value14> 
	<-v15 value15> 
	<-v16 value16> 
	<-v17 value17> 
	<-v18 value18>	
	<-v19 value19>
	<-v20 value20> 
	<-v21 value21> 
	<-v22 value22> 
	<-v23 value23> 
	<-v24 value24> 
	<-v25 value25>	
	<-v26 value26> 
	<-v27 value27> 
	<-v28 value28>
	<-v29 value29> 
	<-v30 value30> 
	<-v31 value31> 
	<-v32 value32> 
	<-v33 value33>	
	<-v34 value34>
	<-v35 value35> 
	<-v36 value36> 
	<-v37 value37> 
	<-v38 value38> 
	<-v39 value39> 
	<-v40 value40>	
	<-v41 value41> 
	<-v42 value42> 
	<-v43 value43> 
	<-v44 value44>
	<-v45 value45> 
	<-v46 value46> 
	<-v47 value47> 
	<-v48 value48> 
 
 
 
	" << std::endl;
}
 
 
int main (int narg, char *arg[])
{
 
	// parse arguments
	GetPot cl (narg, arg);
	if( cl.size() == 1 || cl.search (2,"--help","-h") ) {
		PrintHelp (cl[0]);
		return -1;
	}
 
 
	// scanner les parametres
	double bvalue = cl.follow (1000.0, 2, "-b", "-B");
 
	const char* input   = cl.follow ("", 2, "-i", "-I");
	const char* b0image = cl.follow ("", 2, "-b0", "-B0");
	const char* output  = cl.follow ("", 2, "-o", "-O");
 
	//itk::GISImageIOFactory::RegisterOneFactory();
 
 
	// lire l'image donnee en parametre
	typedef itk::Image<double, 3>           ImageType;
	typedef itk::ImageFileReader<ImageType> ReaderType;
 
	ImageType::Pointer InputImage = 0;
	{
		ReaderType::Pointer reader = ReaderType::New();
		reader->SetFileName ( input );
		try
		{
			reader->Update();
		}
		catch (itk::ExceptionObject &e)
		{
			std::cerr << e;
			return -1;
		}
 
		InputImage = reader->GetOutput(); 
	}
 
 
 
	ImageType::Pointer B0Image = 0;
	{
		ReaderType::Pointer reader = ReaderType::New();
		reader->SetFileName ( b0image );
		try
		{
			reader->Update();
		}
		catch (itk::ExceptionObject &e)
		{
			std::cerr << e;
			return -1;
		}
 
		B0Image = reader->GetOutput();
	}
 
 
 
	typedef itk::DiffusionWeightedImageToADCImageFilter<ImageType, ImageType> ADCCalculatorType;
 
	ADCCalculatorType::Pointer myFilter =  ADCCalculatorType::New();
	myFilter->SetBValue  ( bvalue );
	myFilter->SetB0Image ( B0Image );
	myFilter->SetInput   ( InputImage );
 
	// execute le filtre
	try
	{
		myFilter->Update();
	}
	catch (itk::ExceptionObject &e)
	{
		std::cerr << e;
		return -1;
	}
 
	InputImage = myFilter->GetOutput();
 
 
 
	// ecrire resultat
	typedef itk::ImageFileWriter<ImageType> WriterType;
	WriterType::Pointer writer = WriterType::New();
	writer->SetInput ( InputImage );
	writer->SetFileName ( output );
	try
	{
		writer->Update();
	}
	catch (itk::ExceptionObject &e)
	{
		std::cerr << e;
		return -1;
	}
 
	return 0;
 
}

**Davidbrcz** · 07/01/2010, 19h43

OMG, tu ne connais ni les boucles ni les tableaux ? Car déclarer à coup de copier/coller 48 variables, il faut être courageux (et un peu fou sur les bords)

**TopCao** · 07/01/2010, 19h47

Je ne sais pas récupérer autrement les 48 valeurs passées en arguments de mon prog dans le terminal

**Davidbrcz** · 07/01/2010, 20h25

Bah un tableau, une boucle et c'est reglé. Je laisse koala détailler tout ca, il a une bien meilleure prose que moi (et j'ai pas trop le temps aussi

)

Sinon, tu es sûr de vouloir attaquer des sujets comme les classes templates alors que tu ne maitrise même pas les boucles ?

**TopCao** · 07/01/2010, 21h59

Je n'ai pas le choix !

Mais je suis là pour apprendre

PS : 1. je sais faire des boucles 2. je sais à peu près manipuler des tableaux
Mais je ne sais pas faire de boucle sur la lectures des paramètres de mon programme.

Invité · 07/01/2010, 22h42

Comme tu as six paramètres, tu as intérêt à les mettre dans un tableau de taille fixe (ou un vector), plutot que de les définir comme six variables de noms différents.

Si tu tiens absolument aux noms, tu peux les redéfinir via une enum

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
 
enum {sa,sb,sc,adca,adcb,adcc,nbparams};
double params[nbparams];

comme ça tu peux accéder à tes paramètres par leur nom: params[sa], params[adcc]...

Il va ensuite te falloir définir trois fonctions :

- la fonction que tu minimise : quelque chose comme la somme sur tes 48 points des carrés des écarts entre ta valeur observée et celle prédite par un modèle ayant param[] pour valeurs des parametres... Donc un truc comme

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
 
double fom(double param[],vector<double> &datax, vector<double> &datay)
{
double res=0.0;
for(int i=0;i<datax.size();i++) {
  double v=datay[i]-f(param,datax[i]);
  res+=v*v;
}
return res;
}

(f() ta fonction modèle, datax et datay tes données)

- le gradient en une valeur de param[]: qui renvoie un tableau de six paramètres, sur le même principe que la précédente
- le hessien pour cette même valeur, qui renvoie une matrice, ou si tu préfères, six tableaux de 6 paramètres, idem

Une fois que tu as ces trois éléments, tout devient simple... Vu que ta fonction n'est pas très compliquée, tu as intérêt à tout calculer à la main, et à programmer les formules obtenues... (les dérivées de fonctions linéaires et d'exponentielles, ce n'est pas difficile.

Puis tu peux appliquer l'algorithme de mise à jour des pas à ces trois fonctions. Tu pars d'une valeur initiale param0[], tu calcules tes trois fonctions, et le pas, tu mets à jour, et tu recommences...

Francois

**TopCao** · 08/01/2010, 12h33

Envoyé par fcharton

Comme tu as six paramètres, tu as intérêt à les mettre dans un tableau de taille fixe (ou un vector), plutot que de les définir comme six variables de noms différents.

Si tu tiens absolument aux noms, tu peux les redéfinir via une enum

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
 
enum {sa,sb,sc,adca,adcb,adcc,nbparams};
double params[nbparams];

Merci pour ton aide
Alors si je comprends bien, nbparams a la valeur 6 (sa a la valeur 0) et on crée un tableau nommé params dont les 6 entrées seront toutes en double précision. Maintenant, ce que je ne comprends pas c'est par quel mécanisme celles-ci se verront attribué les noms définis dans enum même s'il me semble que cela fonctionne parfaitement.

Voilà comment je l'ai écrit du coup :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
enum {Sa = 10,Sb = 20,Sc = 30,ADCa = 0.001,ADCb = 0.002,ADCc = 0.003,nbParams};
double vector[nbParams];

J'avance doucement mais surement

**TopCao** · 08/01/2010, 12h47

Envoyé par TopCao

Merci pour ton aide

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
enum {Sa = 10,Sb = 20,Sc = 30,ADCa = 0.001,ADCb = 0.002,ADCc = 0.003,nbParams};
double vector[nbParams];

Je m'auto cite pour préciser qu'on ne peut initialiser les paramètres dans enum qu'avec des entiers....

**Goten** · 08/01/2010, 13h21

Attention les membres de l'enum doivent être des indices, pas contenir les valeurs.... là si tu fais ::
vector[sa] == vector[10] <== kaboum ton vecteur est plus petit.
Ce que te proposait fcharton c'est d'avoir un accés indexé du vecteur avec des les membres de l'enu :

vecteur[sa] == vecteur[0];
vecteur[sa] = 10;

**TopCao** · 08/01/2010, 13h27

Envoyé par Goten

Attention les membres de l'enum doivent être des indices, pas contenir les valeurs.... là si tu fais ::
vector[sa] == vector[10] <== kaboum ton vecteur est plus petit.
Ce que te proposait fcharton c'est d'avoir un accés indexé du vecteur avec des les membres de l'enu :

vecteur[sa] == vecteur[0];
vecteur[sa] = 10;

Ah d'accord, je n'avais pas saisi
Je me corrige

Merci

**TopCao** · 08/01/2010, 13h49

Bon alors la simple somme des moindres carrés ne compile pas...
J'ai encore besoin d'un peu de votre aide

merci !!

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
// ComputeLM
// Calcule les coefficients de la fonction signalDecay
// de type Sa * exp(-b * ADCa) + Sb * exp(-b * ADCb) + Sc * exp(-b * ADCc) + noise
// à partir d'un algorithme de Levenberg Macquardt
// pour n points de coordonnées (x,y)
 
#include <cstdio>
#include <iostream>
#include <cstdlib>
#include <math.h>
using namespace std;
 
// Déclare le nombre de valeurs de b puis sa valeur initiale
int nbbValue = 48;
int b = 0;
 
// Déclare un tableau contenant les valeurs x des points
// devrait prendre nbbValue comme paramètre mais ne marche pas
double dataX[48];
// Déclare un tableau contenant les valeurs y des points
// devrait prendre nbbValue comme paramètre mais ne marche pas
double dataY[48];
 
// Déclare le bruit comme étant une constante
// (ce qui est faux comme chacun sait !)
double noise = 0;
 
// Déclare tous les paramètres de la fonction mathématique
// signalDecay dans un tableau nommé vector de taille fixe (6)
// et initialise les valeurs comme suit (ne marche pô !)
enum {Sa,Sb,Sc,ADCa,ADCb,ADCc,nbParams};
double paramsVector[nbParams];
paramsVector[Sa] = 10;
paramsVector[Sb] = 20;
paramsVector[Sc] = 30;
paramsVector[ADCa] = 0.001;
paramsVector[ADCb] = 0;002;
paramsVector[ADCc] = 0.003;
 
//
// Fonction squareValue - Retourne le carré d'une valeur entrée en argument
//
double squareValue (double valueToSquare)
{
	return valueToSquare * valueToSquare;
}
 
// Fonction de comparaison de la valeur observée et de la valeur prédite
// retourne le carré de la différence
double compareValues (double paramsVector[], vector<double> &datax, vector<double> &datay)
{
	double sum = 0;
	for ( int i = 0 ; i < dataX.size() ; i = dataX[i++])
	{
		double realSignal = dataY[i];
		b = dataX[i];
		double expectedSignal = Sa * exp(-b * ADCa) + Sb * exp(-b * ADCb) + Sc * exp(-b * ADCc) + noise ;
		sum = squareValue(realSignal - expectedSignal);
	}
	return sum;
}
 
 
 
// Corps du programme
int main (int nNumberofArgs, char * pszArgs[])
{
	cout	<< "\nPour le moment, ce programme calcule la somme des différences au carré des valeurs espérées passées dans l'initialisation\n"
			<< "Il retourne la valeur obtenue\n"
			<< "Faire un CC des 48 points sous la forme x1 ... x48 suivis de y1 ... y48 et terminer par un nombre negatif\n"
			<< endl;
 
	// accumulation des suites de nombres
	// les range dans 2 jolis tableaux
 
	for (int i = 0 ; i = 47 ; i++)
	{
		double value = 0;
		cin >> value ;
		if (value < 0)
		{
			break ;
		}
		dataX[i] = value;
	}
	for (int j = 0 ; j = 47 ; j++)
		{
			double value = 0;
			cin >> value ;
			if (value < 0)
			{
				break ;
			}
			dataY[j] = value;
		}
 
 
	// affichage du résultat
		cout << "\n La somme des moindres carrés est : " << compareValues(paramsVector[], vector<double> &datax, vector<double> &datay) << "\n" << endl;
 
 
	cout << "A la prochaine - Kill prog dans 2 secondes..." << endl;
	// Attend pour terminer le programme que l'utilisateur lise
	// le contenu de la fenêtre
	system("SLEEP 2");
	return 0;
}

**Lavock** · 08/01/2010, 14h03

Pourrait tu précisais quelle erreur et la ligne s'il te plait ?

**TopCao** · 08/01/2010, 14h23

Envoyé par Lavock

Pourrait tu précisais quelle erreur et la ligne s'il te plait ?

Bon avec Xcode je n'arrive pas à comprendre le debug mode (le num des lignes ne s'affiche pas et je ne peux pas exporter le journal des erreurs)
Alors voici des snapshots thanks to Mac OsX

**Lavock** · 08/01/2010, 14h42

dataX.size n'a pas de sens ! Tu voulais peut-être écrire datax.size ?

Tu as mis un ";" dans une declaration de paramVector. De plus, un tableau s'initialise comme ça :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

 type table[taille] = {v1,v2,v3,v4,v5,v6};

La tu assigne avant le main, est pas dans une fonction... c'est impossible.

Pour passer paramVector, pas besoin de mettre les [].

**TopCao** · 08/01/2010, 15h36

Envoyé par Lavock

dataX.size n'a pas de sens ! Tu voulais peut-être écrire datax.size ?

Tu as mis un ";" dans une declaration de paramVector. De plus, un tableau s'initialise comme ça :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

 type table[taille] = {v1,v2,v3,v4,v5,v6};

La tu assigne avant le main, est pas dans une fonction... c'est impossible.

Pour passer paramVector, pas besoin de mettre les [].

Alors j'ai corrigé comme suit mais il met dit que datax et vector ne sont pas déclarés

Pour moi dataX.size() retournait le nombre de valeurs dans dataX[]

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
// ComputeLM
// Calcule les coefficients de la fonction signalDecay
// de type Sa * exp(-b * ADCa) + Sb * exp(-b * ADCb) + Sc * exp(-b * ADCc) + noise
// à partir d'un algorithme de Levenberg Macquardt
// pour n points de coordonnées (x,y)
 
#include <cstdio>
#include <iostream>
#include <cstdlib>
#include <math.h>
using namespace std;
 
// Déclare le nombre de valeurs de b puis sa valeur initiale
int nbbValue = 48;
int b = 0;
 
// Déclare un tableau contenant les valeurs x des points
// devrait prendre nbbValue comme paramètre mais ne marche pas
double dataX[48];
// Déclare un tableau contenant les valeurs y des points
// devrait prendre nbbValue comme paramètre mais ne marche pas
double dataY[48];
 
// Déclare le bruit comme étant une constante
// (ce qui est faux comme chacun sait !)
double noise = 0;
 
// Déclare tous les paramètres de la fonction mathématique
// signalDecay dans un tableau nommé vector de taille fixe (6)
// et initialise les valeurs comme suit (ne marche pô !)
enum {Sa,Sb,Sc,ADCa,ADCb,ADCc,nbParams};
double paramsVector[nbParams] = {10,20,30,0.001,0.002,0.003};
 
 
//
// Fonction squareValue - Retourne le carré d'une valeur entrée en argument
//
double squareValue (double valueToSquare)
{
	return valueToSquare * valueToSquare;
}
 
// Fonction de comparaison de la valeur observée et de la valeur prédite
// retourne le carré de la différence
double compareValues (double paramsVector, vector<double> &datax, vector<double> &datay)
{
	double sum = 0;
	for ( int i = 0 ; i < datax.size() ; i = dataX[i++])
	{
		double realSignal = dataY[i];
		b = dataX[i];
		double expectedSignal = Sa * exp(-b * ADCa) + Sb * exp(-b * ADCb) + Sc * exp(-b * ADCc) + noise ;
		sum = squareValue(realSignal - expectedSignal);
	}
	return sum;
}
 
 
 
// Corps du programme
int main (int nNumberofArgs, char * pszArgs[])
{
	cout	<< "\nPour le moment, ce programme calcule la somme des différences au carré des valeurs espérées passées dans l'initialisation\n"
			<< "Il retourne la valeur obtenue\n"
			<< "Faire un CC des 48 points sous la forme x1 ... x48 suivis de y1 ... y48 et terminer par un nombre negatif\n"
			<< endl;
 
	// accumulation des suites de nombres
	// les range dans 2 jolis tableaux
 
	for (int i = 0 ; i = 47 ; i++)
	{
		double value = 0;
		cin >> value ;
		if (value < 0)
		{
			break ;
		}
		dataX[i] = value;
	}
	for (int j = 0 ; j = 47 ; j++)
		{
			double value = 0;
			cin >> value ;
			if (value < 0)
			{
				break ;
			}
			dataY[j] = value;
		}
 
 
	// affichage du résultat
		cout << "\n La somme des moindres carrés est : ";
		cout << compareValues(paramsVector, vector<double> &datax, vector<double> &datay) << "\n" << endl;
 
 
	cout << "A la prochaine - Kill prog dans 2 secondes..." << endl;
	// Attend pour terminer le programme que l'utilisateur lise
	// le contenu de la fenêtre puis appuie sur une touche
	system("SLEEP 2");
	return 0;
}

**Lavock** · 08/01/2010, 15h41

Tu veux avoir quoi par cette instruction ?
Si c'est la taille de dataX, alors il n'y à pas de méthode .size... Tu connais déjà sa taille puisqu'elle est fixe.

Sinon, pense à :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

#include <vector>

Invité · 08/01/2010, 15h54

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
for ( int i = 0 ; i < datax.size() ; i = dataX[i++])

Ca ne peut pas marcher...

Essaye avec...

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
for ( int i = 0 ; i < datax.size() ; i++)

Ensuite, dans

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

double expectedSignal = Sa * exp(-b * ADCa) + Sb * exp(-b * ADCb) + Sc * exp(-b * ADCc) + noise ;

Il faut mettre les valeurs des parametres que tu passes dans Params...

Francois

**TopCao** · 08/01/2010, 17h32

Alors maintenant ma boucle donne ça :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
double compareValues (double paramsVector, vector<double> &datax, vector<double> &datay)
{
	double sum = 0;
	for ( int i = 0 ; i < datax.size() ; i++)
	{
		double realSignal = dataY[i];
		b = dataX[i];
		double expectedSignal = paramsVector[Sa] * exp(-b * paramsVector[ADCa]) + paramsVector[Sb] * exp(-b * paramsVector[ADCb]) + paramsVector[Sc] * exp(-b * paramsVector[ADCc]) + noise ;
		sum = realSignal - expectedSignal;
		sum = sum * sum ;
	}
	return sum;
}

Mais j'ai toujours du mal à capter ce qu'est mon vecteur dans l'histoire du coup.
Je vais y bosser et je reviens !

Algorithme de Levenberg Marquardt

C++

Discussions similaires

Partager

Partager