source: issm/trunk-jpl/src/m/contrib/gravity/vfsa.cpp@ 18434

Last change on this file since 18434 was 18434, checked in by abuzzi, 11 years ago

NEW: added vfsa code from atsu
File size: 13.6 KB
Line 
1#include <iostream>
2#include <cmath>
3#include <fstream>
4#include <string>
5#include <cstdio>
6#include <cstdlib>
7#include <time.h>
8#include <cassert>
9#include <gsl/gsl_multifit.h>
10
11using namespace std;
12
13class Matrix{/*{{{*/
14 private:
15 int M; /*Number of lines */
16 int N; /*Number if Columns */
17 double *values;
18 public:
19 Matrix(int m_in,int n_in){/*{{{*/
20 this->M = m_in;
21 this->N = n_in;
22 this->values = new double[M*N];
23 }/*}}}*/
24 ~Matrix(){/*{{{*/
25 delete [] this->values;
26 }/*}}}*/
27 void Echo(void){/*{{{*/
28 for(int i=0;i<M;i++){
29 for(int j=0;j<N;j++){
30 cout << " " << this->values[i*N+j];
31 }
32 cout << endl;
33 }
34 }/*}}}*/
35 void SetValue(int i,int j,double value){/*{{{*/
36 this->values[i*N+j] = value;
37 }/*}}}*/
38 double GetValue(int i,int j){/*{{{*/
39 return this->values[i*N+j];
40 }/*}}}*/
41 void GetSize(int* pM,int* pN){/*{{{*/
42 *pM = this->M;
43 *pN = this->N;
44 }/*}}}*/
45 double* GetPointer(void){/*{{{*/
46 return this->values;
47 }/*}}}*/
48 void MatrixSum(Matrix* A,Matrix* B){/*{{{*/
49 /*Check that sizes are compatible*/
50 int M_B,N_B,M_A,N_A;
51 B->GetSize(&M_B,&N_B);
52 A->GetSize(&M_A,&N_A);
53 assert(this->M==M_B && this->N==N_B);
54 assert(this->M==M_A && this->N==N_A);
55 for(int i=0;i<this->M;i++){
56 for(int j=0;j<this->N;j++){
57 this->SetValue(i,j,A->GetValue(i,j) + B->GetValue(i,j));
58 }
59 }
60 }/*}}}*/
61 void MatrixDiff(Matrix* A,Matrix* B){/*{{{*/
62 /*Check that sizes are compatible*/
63 int M_B,N_B,M_A,N_A;
64 B->GetSize(&M_B,&N_B);
65 A->GetSize(&M_A,&N_A);
66 assert(this->M==M_B && this->N==N_B);
67 assert(this->M==M_A && this->N==N_A);
68 for(int i=0;i<this->M;i++){
69 for(int j=0;j<this->N;j++){
70 this->SetValue(i,j,A->GetValue(i,j) - B->GetValue(i,j));
71 }
72 }
73 }/*}}}*/
74 void MatrixAbs(Matrix* A){/*{{{*/
75 for(int i=0;i<this->M;i++){
76 for(int j=0;j<this->N;j++){
77 this->SetValue(i,j,fabs(A->GetValue(i,j)));
78 }
79 }
80 }/*}}}*/
81 void MatrixEqual(Matrix* A){/*{{{*/
82 for(int i=0;i<this->M;i++){
83 for(int j=0;j<this->N;j++){
84 this->SetValue(i,j,A->GetValue(i,j));
85 }
86 }
87 }/*}}}*/
88 double MatrixInternSum(){/*{{{*/
89 double sum=0;
90 for(int i=0;i<this->M;i++){
91 for(int j=0;j<this->N;j++){
92 sum+=this->GetValue(i,j);
93 }
94 }
95 return sum;
96 }/*}}}*/
97 void ExtractLine(Matrix* A,int i){/*{{{*/
98 /* Check that the size of A is compatible */
99 int M_A,N_A;
100 A->GetSize(&M_A,&N_A);
101 assert(M_A==1 && this->N==N_A);
102 for(int j=0;j<this->N;j++){
103 A->SetValue(0,j,this->GetValue(i,j));
104 }
105 }/*}}}*/
106 void ExtractColumn(Matrix* A,int j){/*{{{*/
107 /* Check that the size of A is compatible */
108 int M_A,N_A;
109 A->GetSize(&M_A,&N_A);
110 assert(N_A==1 && this->M==M_A);
111 for(int i=0;i<this->M;i++){
112 A->SetValue(i,0,this->GetValue(i,j));
113 }
114 }/*}}}*/
115 void AddNumber(double a){/*{{{*/
116 for(int i=0;i<this->M;i++){
117 for(int j=0;j<this->N;j++){
118 this->SetValue(i,j,this->GetValue(i,j) + a);
119 }
120 }
121 }/*}}}*/
122};/*}}}*/
123
124/*Local prototypes{{{*/
125void makep(Matrix *Pobs,int nx,int ny, int dx, int dy);
126void vec2grid(Matrix *V,Matrix *V1,Matrix *V2,int nx, int ny);
127void msplit( Matrix *m, Matrix *m1,Matrix *m2,double dlevel);
128void plouff(Matrix *g,Matrix *Pobs,Matrix *Pp,Matrix * mesh,Matrix *rho,int dx,int dy, int dn,int m,int n,int l);
129void vec2gridsimple(Matrix *V,Matrix *V1,int nx, int ny);
130void reshape(Matrix* V,Matrix* V1,int nx,int ny);
131double misfit(Matrix* m0,Matrix* evalid,Matrix* gobs,int dlevel,Matrix* Pobs,Matrix* xobs,Matrix* yobs,Matrix* Pp,Matrix* rho1, Matrix* rho2,int dx,int dy,int dn,int nx,int ny, int mx,int my);
132void GSLsquarefit(Matrix** pX,Matrix* A,Matrix* B);
133/*}}}*/
134
135int main(){/*{{{*/
136
137 /* Seed the random number generator {{{*/
138 srand (time(NULL)); /*}}}*/
139 /* Define the variables {{{*/
140
141 int dx = 500; /* prism dimension in x-direction */
142 int dy = 500; /* prism dimension in y-direction */
143 int mx = 39; /* number of prisms in x-direction */
144 int my = 60; /* number of prisms in y-direction */
145 int nx = 39; /* number of data points in x-direction */
146 int ny = 60; /* number of data points in y-direction */
147 int dn = 15000; /* distance for neighbouting prisms for gravity calculation */
148 double ptval = 100.; /* max. amount to perturb model */
149 double ptval2 = 100.;
150
151 Matrix *Pobs=new Matrix(nx*ny,2); /* data positions */
152 makep(Pobs,nx,ny,dx,dy);
153// Pobs->Echo();
154
155
156 Matrix *Pp=new Matrix(mx*my,2); /* data positions */
157 makep(Pp,mx,my,dx,dy);
158// Pp->Echo();
159
160 double rhoi = 890; /* ice density */
161 double rhow = 1030; /* water density */
162// double rhos = 2013; /* sediment density */
163 double rhoc = 2670; /* bedrock density */
164
165 Matrix *Rho = new Matrix(1,2);
166 Rho->SetValue(0,0,rhoi);
167 Rho->SetValue(0,1,rhow);
168 Matrix *rho1 = new Matrix(1,3);
169 rho1->SetValue(0,0,rhoi);
170 rho1->SetValue(0,1,rhow);
171 rho1->SetValue(0,2,rhoc);
172 Matrix *rho2 = new Matrix(1,2);
173 rho2->SetValue(0,0,rhoi-rhoc);
174 rho2->SetValue(0,1,rhow-rhoc);
175
176 double ctr=1; /* parameter for filtering */
177 double sd=0.1;
178
179 Matrix *xobs= new Matrix(ny,nx);
180 Matrix *yobs= new Matrix(ny,nx);
181
182 vec2grid(Pobs,xobs,yobs,nx,ny);
183// xobs->Echo();
184// yobs->Echo();
185
186 double dlevel= 860; /* mean level of data acquisition */
187
188 double mmax = 1000; /* max value for layer interfaces */
189 double mmax2 = 1000;
190 double mmax3 = 1000;
191
192 /* control parameter for temperature schedule */
193
194 double ca=0.9; /* for acceptance */
195 double cm=0.5; /* for model perturbation */
196 /*}}}*/
197 /* load the data {{{*/
198
199 /* Observed gravity anomaly */
200
201 ifstream file1("store_fa500_36s.txt");
202 Matrix * gobs= new Matrix(nx*ny,1);
203 double inputnumber;
204 for(int i=0;i<ny*nx; i++){
205 file1 >> inputnumber;
206 gobs->SetValue(i,0, inputnumber*1e-5);
207 }
208 file1.close();
209// gobs->Echo();
210
211 /* load data about the ice thickness */
212
213 ifstream file2("store_flag_icethick500.txt");
214 Matrix * icethick= new Matrix(mx*my,1);
215 for(int s=0;s<mx*my; s++){
216 file2 >> inputnumber;
217 icethick->SetValue(s,0,inputnumber);
218 }
219 file2.close();
220// icethick->Echo();
221
222 /* load the batimethry data */
223
224 ifstream file3("store_flag_bathy500.txt");
225 Matrix * bathy= new Matrix(mx*my,1);
226 for(int s=0;s<mx*my; s++){
227 file3 >> inputnumber;
228 bathy->SetValue(s,0,inputnumber);
229 }
230 file3.close();
231// bathy->Echo();
232
233 /* id of grid to evaluate misfit */
234
235
236 ifstream file4("store_flag_eval500.txt");
237 Matrix * evalid= new Matrix(mx*my,1);
238 for(int s=0;s<mx*my; s++){
239 file4 >> inputnumber;
240 evalid->SetValue(s,0,inputnumber);
241 }
242 file4.close();
243// evalid->Echo();
244
245 /* initial guess of the model */
246
247 ifstream file5("m0_140114b.txt");
248 Matrix * mesh_ini= new Matrix(mx*my,3);
249 for(int s=0;s<mx*my; s++){
250 for(int j=0;j<3;j++){
251 file5 >> inputnumber;
252 mesh_ini->SetValue(s,j,inputnumber);
253 }
254 }
255 file5.close();
256// mesh_ini->Echo();
257 /*}}}*/
258 double q=0;
259 q=misfit(mesh_ini,evalid,gobs,dlevel,Pobs,xobs,yobs,Pp,rho1,rho2,dx,dy,dn,nx,ny,mx,my);
260 cout<< q;
261/* Matrix* G=new Matrix(4,3);
262 Matrix* df=new Matrix(4,1);
263 Matrix* o=new Matrix(4,1);
264 for(int i= 0;i<4;i++){
265 for(int j=0;j<3;j++){
266 if(j==2){G->SetValue(i,j,1);}
267 else{G->SetValue(i,j,i+j);}
268 }
269 df->SetValue(i,0,i*1e-3);
270 }
271 G->SetValue(0,0,1);
272 G->SetValue(0,1,2);
273
274 Matrix* M = NULL;
275 GSLsquarefit(&M,G,df);
276 for(int i= 0;i<4;i++){
277 o->SetValue(i,0,G->GetValue(i,0)*M->GetValue(0,0)+G->GetValue(i,1)*M->GetValue(1,0)+G->GetValue(i,2)*M->GetValue(2,0));}
278
279 G->Echo();
280 df->Echo();
281 M->Echo();
282 o->Echo();*/
283
284 return 0;
285}/*}}}*/
286
287void GSLsquarefit(Matrix** pX,Matrix* A,Matrix* B){/*{{{*/
288
289 /*GSL Matrices and vectors: */
290 int M,N;
291 double chisq;
292
293 /*Get system size*/
294 A->GetSize(&M,&N);
295
296 /*Initialize gsl matrices and vectors: */
297 gsl_matrix* a = gsl_matrix_alloc(M,N);
298 for(int i=0;i<M;i++){
299 for(int j=0;j<N;j++){
300 gsl_matrix_set (a,i,j,A->GetValue(i,j));
301 }
302 }
303 gsl_vector* b = gsl_vector_alloc(M);
304 for(int i=0;i<M;i++){
305 gsl_vector_set(b,i,B->GetValue(i,0));
306 }
307
308 gsl_vector* x = gsl_vector_alloc(N);
309 gsl_matrix* cov = gsl_matrix_alloc(N,N);
310
311 /*Least square fit: */
312 gsl_multifit_linear_workspace* work = gsl_multifit_linear_alloc(M,N);
313 gsl_multifit_linear (a, b, x, cov, &chisq, work);
314 gsl_multifit_linear_free (work);
315
316 /*Clean up and assign output pointer*/
317 Matrix* X = new Matrix(N,1);
318 for(int j=0;j<N;j++){
319 X->SetValue(j,0,gsl_vector_get(x,j));
320 }
321 *pX = X;
322
323 gsl_matrix_free(a);
324 gsl_vector_free(b);
325 gsl_matrix_free(cov);
326
327}/*}}}*/
328void makep(Matrix *Pobs,int nx,int ny, int dx, int dy){/*{{{*/
329 for(int i=0;i<ny;i++){
330 for(int j=0;j<nx;j++){
331 Pobs->SetValue(j+nx*i,0,j*dx);
332 Pobs->SetValue(j+nx*i,1,i*dy);
333 }
334 }
335}/*}}}*/
336void vec2grid(Matrix *V,Matrix *V1,Matrix *V2,int nx, int ny){/*{{{*/
337 for(int i=0;i<ny;i++){
338 for (int j=0;j<nx;j++){
339 V1->SetValue(i,j, V->GetValue(j+nx*i,0));
340 V2->SetValue(i,j, V->GetValue(j+nx*i,1));
341 }
342 }
343}/*}}}*/
344void msplit( Matrix *m, Matrix *m1,Matrix *m2,double dlevel){/*{{{*/
345 int sizem1,sizem2;
346 m->GetSize(&sizem1,&sizem2);
347 for(int i=0;i<sizem1;i++){
348 for(int j=0;j<sizem2+1;j++){
349 if(j<sizem2){
350 m1->SetValue(i,j,1e-10*(sizem2+1-j));
351 m2->SetValue(i,j,m->GetValue(i,j));
352 if(m->GetValue(i,j)<0){
353 m1->SetValue(i,j,m->GetValue(i,j));
354 m2->SetValue(i,j,i*1e-10);
355 }
356 }
357 else{
358 m1->SetValue(i,j,1e-10);
359 }
360 }
361 }
362 m1->AddNumber(dlevel);
363 m2->AddNumber(dlevel);
364}/*}}}*/
365void plouff(Matrix *g,Matrix *Pobs,Matrix *Pp,Matrix * mesh,Matrix *rho,int dx,int dy, int dn,int m,int n,int l){/*{{{*/
366 double gg=6.673e-11;
367 int si,sj,id,s;
368 double R,Q,P;
369 Matrix *xp= new Matrix(1,2);
370 Matrix *yp= new Matrix(1,2);
371 Matrix *xpp= new Matrix(1,2);
372 Matrix *ypp= new Matrix(1,2);
373 Matrix *U= new Matrix(l,4);
374 Matrix *U1=new Matrix(1,4);
375 Matrix *U2=new Matrix(1,4);
376 Matrix *gl= new Matrix(1,l-1);
377 bool test=true;
378 for(int c=0;c<n;c++){
379 g->SetValue(c,0,0);
380 for(int a=0;a<m;a++){
381 test=true;
382 xp->SetValue(0,0,Pp->GetValue(a,0)-Pobs->GetValue(c,0));
383 xp->SetValue(0,1,Pp->GetValue(a,0)-Pobs->GetValue(c,0)+dx);
384 if(xp->GetValue(0,0)<0 && xp->GetValue(0,0)<xp->GetValue(0,1) && xp->GetValue(0,0)*xp->GetValue(0,1)>=0){
385 xpp->SetValue(0,0,xp->GetValue(0,1));
386 xpp->SetValue(0,1,xp->GetValue(0,0));
387 xp->MatrixAbs(xpp);
388 }
389 yp->SetValue(0,0,Pp->GetValue(a,1)-Pobs->GetValue(c,1));
390 yp->SetValue(0,1,Pp->GetValue(a,1)-Pobs->GetValue(c,1)+dy);
391 if(yp->GetValue(0,0)<0 && yp->GetValue(0,0)<yp->GetValue(0,1) && yp->GetValue(0,0)*yp->GetValue(0,1)>=0){
392 ypp->SetValue(0,0,yp->GetValue(0,1));
393 ypp->SetValue(0,1,yp->GetValue(0,0));
394 yp->MatrixAbs(ypp);
395 }
396 P=sqrt(xp->GetValue(0,0)*xp->GetValue(0,0)+yp->GetValue(0,0)*yp->GetValue(0,0));
397 if(P>dn){
398 test=false;
399 for(int i=0;i<l-1;i++){
400 gl->SetValue(0,i,0);
401 }
402 }
403 if(test==true){
404 si=1;
405 sj=1;
406 id=0;
407 for(int i=0;i<2;i++){
408 si*=-1;
409 for(int j=0;j<2;j++){
410 si*=-1;
411 s=si*sj;
412 for(int k=0;k<l;k++){
413 R=sqrt(xp->GetValue(0,i)*xp->GetValue(0,i)+yp->GetValue(0,j)*yp->GetValue(0,j)+mesh->GetValue(a,k)*mesh->GetValue(a,k));
414 Q=atan(xp->GetValue(0,i)*yp->GetValue(0,j)/(mesh->GetValue(a,k)*R));
415 U->SetValue(k,id,s*(mesh->GetValue(a,k)*Q-xp->GetValue(0,i)*log(R+yp->GetValue(0,j))-yp->GetValue(0,j)*log(R+xp->GetValue(0,i))));
416 }
417 id++;
418 }
419 }
420 for(int b=0;b<l-1;b++){
421 U->ExtractLine(U1,b);
422 U->ExtractLine(U2,b+1);
423 gl->SetValue(0,b,rho->GetValue(0,b)*(U1->MatrixInternSum()*(-1)+U2->MatrixInternSum()));
424 }
425 }
426 g->SetValue(c,0,g->GetValue(c,0)+gg*gl->MatrixInternSum());
427 }
428 }
429 delete xp;
430 delete yp;
431 delete xpp;
432 delete ypp;
433 delete gl;
434 delete U;
435 delete U1;
436 delete U2;
437 }/*}}}*/
438void vec2gridsimple(Matrix *V,Matrix *V1,int nx, int ny){/*{{{*/
439 for(int i=0;i<ny;i++){
440 for (int j=0;j<nx;j++){
441 V1->SetValue(i,j, V->GetValue(j+nx*i,0));
442 }
443 }
444}/*}}}*/
445void reshape(Matrix* V,Matrix* V1,int nx,int ny){/*{{{*/
446 for (int i=0;i<ny;i++){
447 for(int j=0;j<nx;j++){
448 V1->SetValue(j+nx*i,0,V->GetValue(i,j));
449 }
450 }
451}/*}}}*/
452double misfit(Matrix* m0,Matrix* evalid,Matrix* gobs,int dlevel,Matrix* Pobs,Matrix* xobs,Matrix* yobs,Matrix* Pp,Matrix* rho1, Matrix* rho2,int dx,int dy,int dn,int nx,int ny, int mx,int my){/*{{{*/
453 Matrix* m1=new Matrix(mx*my,4);
454 Matrix* m2=new Matrix(mx*my,3);
455 Matrix* g1=new Matrix(nx*ny,1);
456 Matrix* g2=new Matrix(nx*ny,1);
457 Matrix* g=new Matrix(nx*ny,1);
458 Matrix* gcalgr=new Matrix(ny,nx);
459 Matrix* gcalvec=new Matrix(mx*my,1);
460 Matrix* df=new Matrix(mx*my,1);
461 Matrix* G=new Matrix(mx*my,3);
462 double a=0;
463 double b=0;
464 double e=0;
465 msplit(m0,m1,m2,dlevel);
466 plouff(g1,Pobs,Pp,m1,rho1,dx,dy,dn,mx*my,nx*ny,4);
467 plouff(g2,Pobs,Pp,m2,rho2,dx,dy,dn,mx*my,nx*ny,3);
468 g->MatrixSum(g1,g2);
469 vec2gridsimple(g,gcalgr,nx,ny);
470 reshape(gcalgr,gcalvec,mx,my);
471 for (int i=0;i<mx*my;i++){
472 df->SetValue(i,0,evalid->GetValue(i,0)*(gobs->GetValue(i,0)-gcalvec->GetValue(i,0)));
473 G->SetValue(i,0,evalid->GetValue(i,0)*Pobs->GetValue(i,0));
474 G->SetValue(i,1,evalid->GetValue(i,0)*Pobs->GetValue(i,1));
475 G->SetValue(i,2,evalid->GetValue(i,0));
476 }
477 Matrix* M = NULL;
478 GSLsquarefit(&M,G,df);
479
480 for (int i=0;i<my;i++){
481 for(int j=0;j<mx;j++){
482 gcalgr->SetValue(i,j,gcalgr->GetValue(i,j)+xobs->GetValue(i,j)*M->GetValue(0,0)+yobs->GetValue(i,j)*M->GetValue(1,0)+M->GetValue(2,0));
483 }
484 }
485 reshape(gcalgr,g,mx,my);
486 for (int i=0;i<mx*my;i++){
487 a=a+fabs(evalid->GetValue(i,0)*(gobs->GetValue(i,0)-g->GetValue(i,0)));
488 b=b+fabs(evalid->GetValue(i,0)*(gobs->GetValue(i,0)+g->GetValue(i,0)));
489 }
490 e=2*a/(a+b);
491 return e;
492}/*}}}*/
Note: See TracBrowser for help on using the repository browser.