Algorithme de Levenberg Marquardt

Version imprimable

06/01/2010, 22h26
TopCao

Algorithme de Levenberg Marquardt

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
07/01/2010, 02h13
Klaim

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.

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:

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:

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:

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;
 
}

07/01/2010, 19h43
Davidbrcz

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)
07/01/2010, 19h47
TopCao

Je ne sais pas récupérer autrement les 48 valeurs passées en arguments de mon prog dans le terminal
07/01/2010, 20h25
Davidbrcz

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 :aie:)

Sinon, tu es sûr de vouloir attaquer des sujets comme les classes templates alors que tu ne maitrise même pas les boucles ?
07/01/2010, 21h59
TopCao

Je n'ai pas le choix !

Mais je suis là pour apprendre :mouarf:

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.

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:

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:

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

08/01/2010, 12h33
TopCao
Citation:
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:

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:

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 :ccool:
08/01/2010, 12h47
TopCao
Citation:
Envoyé par TopCao

Merci pour ton aide

Code:

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....
08/01/2010, 13h21
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;
08/01/2010, 13h27
TopCao

Citation:

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

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

merci !!

Code:

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;
}

08/01/2010, 14h03
Lavock

Pourrait tu précisais quelle erreur et la ligne s'il te plait ?
08/01/2010, 14h23
TopCao

3 pièce(s) jointe(s)

Citation:

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
08/01/2010, 14h42
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:

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 [].

Citation:

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:

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:

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;
}

08/01/2010, 15h41
Lavock

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:

#include <vector>
08/01/2010, 15h54
Invité
Code:

1 2 for ( int i = 0 ; i < datax.size() ; i = dataX[i++])
Ca ne peut pas marcher...

Essaye avec...
Code:

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

Code:

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

Alors maintenant ma boucle donne ça :

Code:

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 !

Youpi !!

Mon programme compile enfin !!

Mais il n'en est pas pour autant moins buggé...
Code:

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 // // // ComputeLM - Compute Levenberg Macquardt // Le but final : calcule les coefficients de la fonction expectedFonction() // à partir d'un algorithme de Levenberg Macquardt // pour n points de coordonnées (x,y) // // #include <cstdio> #include <iostream> #include <cstdlib> #include <math.h> #include <vector> using namespace std; double noise = 0; // Déclare le bruit comme étant une constante, ce qui est faux comme chacun sait ! enum {a1,a2,a3,b1,b2,b3,nbParams}; // Déclare tous les paramètres de la fonction mathématique signalDecay double param[nbParams] = {1,2,3,0.1,0.2,0.3}; double expectedFonction (double params[] , double x) { return param[a1] * exp(-x * param[b1]) + param[a2] * exp(-x * params[b2]) + params[a3] * exp(-x * params[b3]) + noise ; } // Fonction de comparaison de la valeur observée et de la valeur prédite // retourne le carré de la différence double compareValues (double param[], vector<double> &datax, vector<double> &datay) { double res = 0.0; for ( int i = 0 ; i < datax.size() ; i++) { double sub = datay[i] - expectedFonction (param, datax[i]); res = sub * sub; } return res; } // Corps du programme int main (int nNumberofArgs, char * pszArgs[]) { vector<double> datax ; vector<double> datay ; 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 j = 0 ; j = 47 ; j++) { double value = 0; cin >> value ; if (value < 0) { break ; } datax[j] = value; } for (int k = 0 ; k = 47 ; k++) { double value = 0; cin >> value ; if (value < 0) { break ; } datay[k] = value; } // affichage du résultat cout << "\n La somme des moindres carrés est : "; cout << compareValues(param, datax, datay) ; cout << "\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; }
Citation:

pcnrad33:computeLM_bin alexandre$ ./ComputeLM

Pour le moment, ce programme calcule la somme des différences au carré des valeurs espérées passées dans l'initialisation
Il retourne la valeur obtenue
Faire un CC des 48 points sous la forme x1 ... x48 suivis de y1 ... y48 et terminer par un nombre negatif

0 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 6 4,764753586 3,805805753 3,058150472 2,47256061 2,011680023 1,647096725 1,357148516 1,12527586 0,938783974 0,787911213 0,665126903 0,564601286 0,481804684 0,41320382 0,356031286 0,308110167 0,267720304 0,233496076 0,204348059 0,179402818 0,157956496 0,139438942 0,123385871 0,109417205 0,097220146 0,086535912 0,077149292 0,068880392 0,061578087 0,055114802 0,049382336 0,044288505 0,039754428 0,035712331 0,032103757 0,028878093 0,025991369 0,023405261 0,021086263 0,019004997 0,017135638 0,015455428 0,013944265 0,012584358 0,011359929 0,010256961 0,009262982 0,008366877 -1
/Developer/SDKs/MacOSX10.6.sdk/usr/include/c++/4.2.1/debug/vector:198:
error: attempt to subscript container with out-of-bounds index 47, but
container only holds 0 elements.

Objects involved in the operation:
sequence "this" @ 0x0x7fff5fbffa10 {
type = NSt7__debug6vectorIdSaIdEEE;
}
Abort trap

Je dois mal parametrer mes fonctions...

08/01/2010, 21h02
Invité
Code:

for (int j = 0 ; j = 47 ; j++)

Ca va pas bien se passer, ca...

Ton vrai problème, c'est que tu déclares tes vecteur sans leur donner de taille,
Code:

1 2 3 vector<double> datax ; vector<double> datay ;
mais que tu essaies d'en accéder les éléments... Pour ajouter des éléments dans un vector<>, il faut soit lui donner une taille au départ
Code:

1 2 3 vector<double> datax(48) ; vector<double> datay(48) ;
soit effectuer des push_back() au fur et à mesure que tu ajoutes des données, par exemple, dans la boucle des datax
Code:

1 2 datax.push_back(val);
où val est la valeur que tu viens de lire (et que tu veux ajouter à ton tableau de données)

Francois
08/01/2010, 21h21
JolyLoic

Honnêtement, et sans méchanceté aucune, je pense que tu n'as pas actuellement le niveau nécessaire pour programmer. Tout ce à quoi tu pourras aboutir ici c'est à avoir un programme qui marche, mais que tu ne maîtrises pas vraiment, et avec des notions un peu floues voire erronées sur comment il marche.

Je pense qu'il pourrait être intéressant si tu dois faire ce genre d'activité de manière répétée (c'est à dire si ce programme n'est pas le seul que tu comptes faire au cours des 5 prochaines années) de prendre les choses plus progressivement en commençant par l'apprentissage de la programmation et d'un langage (éventuellement en parallèle du développement de ce code précis, pour lequel j'imagine que tu dois avoir des deadlines).

Des bons livres existent pour ça, les suivre te donnerait probablement une base bien plus solide.
08/01/2010, 22h17
koala01

Citation:

Envoyé par Davidbrcz

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 :aie:)

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

Pourquoi c'est toujours à moi qu'on pense dans ces cas là :question::ouin:

PS: modérateurs, n'hésitez pas à modérer :D
11/01/2010, 11h04
TopCao

Citation:

Envoyé par JolyLoic

Honnêtement, et sans méchanceté aucune, je pense que tu n'as pas actuellement le niveau nécessaire pour programmer. Tout ce à quoi tu pourras aboutir ici c'est à avoir un programme qui marche, mais que tu ne maîtrises pas vraiment, et avec des notions un peu floues voire erronées sur comment il marche.

Je pense qu'il pourrait être intéressant si tu dois faire ce genre d'activité de manière répétée (c'est à dire si ce programme n'est pas le seul que tu comptes faire au cours des 5 prochaines années) de prendre les choses plus progressivement en commençant par l'apprentissage de la programmation et d'un langage (éventuellement en parallèle du développement de ce code précis, pour lequel j'imagine que tu dois avoir des deadlines).

Des bons livres existent pour ça, les suivre te donnerait probablement une base bien plus solide.

Je suis totalement d'accord avec toi et je ne le prends pas mal du tout. Je dois "maitriser" un langage de programmation d'ici fin mai tout en avançant sur mes travaux de recherche en parallèle.
J'ai donc commencé à lire C++ pour les Nuls (une version de 2004) et je lis quelques tutos en fonctions de mes besoins.
Je veux bien des conseils sur les ouvrage essentiels qu'il faut avoir dans sa bibliothèque (et avoir lu surtout) pour se démerder seul en C++. Je rappelle
que je bosse sous Xcode (mais je peux télécharger un autre compilateur si on me le conseille), pour le moment en ligne de commande (je n'ai pas besoin d'interface visuelle) et que je dois intégrer ITK dans une seconde partie pour manipuler mes données (mais pour ça je peux directement me faire aider par un pro). De même si vous avez de bons tutos je suis preneur.
C'est un peu difficile de faire le tri parmi tous les livres et tutos dispos...

Merci pour vos précieux conseils. :ccool:
11/01/2010, 12h22
JolyLoic

Si l'anglais ne te pose pas trop de problèmes, un livre rapide (faut parfois un peu s'accrocher) mais bien expliqué et expliquant les bonnes choses est Accelerated C++.

Tu peux regarder le tutoriel d'Alp sur comment commencer en C++ :
http://alp.developpez.com/tutoriels/debuter-cpp/#LII
11/01/2010, 14h05
TopCao

Citation:

Envoyé par JolyLoic

Si l'anglais ne te pose pas trop de problèmes, un livre rapide (faut parfois un peu s'accrocher) mais bien expliqué et expliquant les bonnes choses est Accelerated C++.

Tu peux regarder le tutoriel d'Alp sur comment commencer en C++ :
http://alp.developpez.com/tutoriels/debuter-cpp/#LII

Ok je vais essayer de me procurer ce bouquin.
Merci
12/01/2010, 23h05
TopCao

Bon alors voilà
J'ai programmé plusieurs fonctions :

f(x) : retourne la valeur de la fonction en prenant les paramètres a1 etc.
smc : retourne la somme des moindre carré pour l'ensemble de mes valeurs de x en prenant les paramètres a1 etc.
grad_a1(x) à grad_b3(x) : retournent la valeur du gradient d'un paramètre pour un x donné
grad2_a1b1 etc : retourne la valeur de la dérivée seconde par paire de paramètres dans une matrice 6x6 (Hessien)

Comment désormais minimiser ma fonction smc avec tout ça ?? De plus si j'ai bien compris il me faut un paramètre lamba (?)

PS : j'attends mon livre avec impatience !
13/01/2010, 16h06
TopCao

up !

Si quelqu'un peut m'aider pour minimiser la fonction smc en tenant compte des gradients et du Hessien...