Changeset 16


Ignore:
Timestamp:
03/22/11 10:35:04 (13 years ago)
Author:
jbrlod
Message:

make v1.1

Location:
trunk
Files:
4 added
47 deleted
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/data_out/plot_results.m

    r15 r16  
     1f=netcdf('exp_T.nc'); 
     2 
     3addpath('../scripts/matlab_toolbox'); 
    14sparam=cell(0,5); 
    25sparam(end+1,:)={'t',[12:1:24],[0:2:24],[-1:0.2:1],[-1:0.2:1]}; 
    … …  
    47sparam(end+1,:)={'u',[-0.2:0.02:0.2],[-0.2:0.015:0.2],[-0.14:0.02:0.14],[-0.14:0.02:0.14]}; 
    58sparam(end+1,:)={'v',[-0.2:0.02:0.2],[-0.2:0.015:0.2],[-0.14:0.02:0.14],[-0.14:0.02:0.14]}; 
    6 f=netcdf('exp_T.nc'); 
     9sparam(end+1,:)={'ssh',[-0.8:0.05:0.2],[-0.2:0.015:0.2],[-0.14:0.02:0.14],[-0.14:0.02:0.14]}; 
    710 
    811suff=cell(0,2); 
    … …  
    1316suff(end+1,:)={'a_c_true','parametre vrai (TP)'}; 
    1417suff(end+1,:)={'b_obs_48','forward du TP'}; 
     18 
     19suff2=cell(0,2); 
     20suff2(end+1,:)={'n_c_init','first guess (FG)'}; 
     21suff2(end+1,:)={'b_forw0_48','forward du FG'}; 
     22suff2(end+1,:)={'n_c_fin','param control final (AP)'}; 
     23suff2(end+1,:)={'b_forwfin_48','forward du PA'}; 
     24suff2(end+1,:)={'n_c_true','parametre vrai (TP)'}; 
     25suff2(end+1,:)={'b_obs_48','forward du TP'}; 
     26 
    1527 
    1628delta=cell(0,2); 
    … …  
    3850  for k=1:length(suff) 
    3951     
     52   if strcmp(sparam{j,1},'ssh') 
     53     param=f{[sparam{j,1} suff2{k,1}]}(:,:,:); 
     54      
     55      param=squeeze(param); 
     56     paramy{k}=mat2yao(param); 
     57     paramy{k}=mat2yao(param); 
     58     paramy{k}(:,1:2)=fliplr(paramy{k}(:,1:2)); 
     59     paramy{k}=sortrows(paramy{k},[1 2]); 
     60     figure(j); 
     61    subplot(3,2,k) 
     62    % contourf(longi(2:end-1,2:end-1),lati(2:end-1,2:end-1),param(2:end-1,2:end-1),sparam{j,2}); 
     63    % caxis([sparam{j,2}(1) sparam{j,2}(end)]) 
     64     %xlabel('longitude','Fontsize',9) 
     65     %ylabel('latitude','Fontsize',9) 
     66      plot_nemo_3D(paramy{k},longi,lati,zlev,1,sparam{j,2},1); 
     67 
     68     
     69    title(suff2{k,2}); 
     70     if (k==1)  
     71      colorbar 
     72     end 
    4073    
    41     param=f{[sparam{j,1} suff{k,1}]}(:,:,:,:); 
    42     param=squeeze(param); 
    43     paramy{k}=mat2yao(param); 
    44     paramy{k}(:,1:3)=fliplr(paramy{k}(:,1:3)); 
    45     paramy{k}=sortrows(paramy{k},[1 2 3]); 
    46     
     74   else 
     75     param=f{[sparam{j,1} suff{k,1}]}(:,:,:,:); 
     76     param=squeeze(param); 
     77     paramy{k}=mat2yao(param); 
     78     paramy{k}(:,1:3)=fliplr(paramy{k}(:,1:3)); 
     79     paramy{k}=sortrows(paramy{k},[1 2 3]); 
    4780    figure(j); 
    4881    subplot(3,2,k) 
    4982    plot_nemo_3D(paramy{k},longi,lati,zlev,1,sparam{j,2},1); 
    50     title(suff{k,2}) 
     83    title(suff{k,2})  
     84   
     85    
    5186    if (k==1)  
    5287      colorbar 
    … …  
    5994      colorbar 
    6095    end 
    61     
     96   end %if ~ssh 
    6297     
    63     end     
     98  end %for k    
     99   
    64100     %%Graphe des différences 
    65101    for k=1:size(delta,1) 
    66102       
    67     param=paramy{delta{k,1}(1)}; 
    68     param(:,end)=param(:,end)-paramy{delta{k,1}(2)}(:,end); 
    69      
    70     figure(j+2*size(sparam,1)) 
    71     subplot(2,2,k) 
    72     plot_nemo_3D(param,longi,lati,zlev,1,sparam{j,4},1); 
    73     colormap(cmapdif); 
    74     title(delta{k,2}); 
    75     if (k==2) 
    76       colorbar 
    77     end 
    78      
    79       figure(j+3*size(sparam,1)) 
    80     subplot(2,2,k) 
    81      plot_nemo_3D(param,longi,lati,zlev,lat0,sparam{j,5},2); 
    82     colormap(cmapdif); 
    83     title(delta{k,2}); 
    84     if (k==2) 
    85       colorbar 
    86     end 
    87     fprintf(1,[sparam{j,1} ' ' delta{k,2} '\n']); 
    88     [N,RPD,URPD,MAD,RMS,slope,intercept,r2]=make_stats(paramy{delta{k,1}(1)}(:,end),paramy{delta{k,1}(2)}(:,end)) 
     103       
     104      param=paramy{delta{k,1}(1)}; 
     105      param(:,end)=param(:,end)-paramy{delta{k,1}(2)}(:,end); 
     106       
     107      figure(j+2*size(sparam,1)) 
     108      subplot(2,2,k) 
     109      plot_nemo_3D(param,longi,lati,zlev,1,sparam{j,4},1); 
     110      colormap(cmapdif); 
     111      title(delta{k,2}); 
     112      if (k==2) 
     113        colorbar 
     114      end 
     115      if ~strcmp(sparam{j,1},'ssh') 
     116        figure(j+3*size(sparam,1)) 
     117        subplot(2,2,k) 
     118        plot_nemo_3D(param,longi,lati,zlev,lat0,sparam{j,5},2); 
     119        colormap(cmapdif); 
     120        title(delta{k,2}); 
     121        if (k==2) 
     122          colorbar 
     123        end 
     124      end % if ~ssh 
     125       
     126      if strcmp(sparam{j,1},'ssh') 
     127        n1=max(paramy{delta{k}(1)}(:,1)); 
     128        n2=max(paramy{delta{k}(1)}(:,2)); 
     129        iok=find(paramy{delta{k}(1)}(:,1)>1 & paramy{delta{k}(1)}(:,1)<n1 & paramy{delta{k}(1)}(:,2)>1 & ... 
     130                 paramy{delta{k}(1)}(:,2)<n2); 
     131      else 
     132        n1=max(paramy{delta{k}(1)}(:,1)); 
     133        n2=max(paramy{delta{k}(1)}(:,2)); 
     134        n3=max(paramy{delta{k}(1)}(:,3)); 
     135        iok=find(paramy{delta{k}(1)}(:,1)>1 & paramy{delta{k}(1)}(:,1)<n1 & paramy{delta{k}(1)}(:,2)>1 & ... 
     136                 paramy{delta{k}(1)}(:,2)<n2 & paramy{delta{k}(1)}(:,3)<n3); 
     137      end 
     138       
     139       
     140         %  fprintf(1,[int2str(j) sparam{j,1} ' ' delta{k,2} '\n']); 
     141      [N(j,k),RPD(j,k),URPD(j,k),MAD(j,k),RMS(j,k),slope(j,k),intercept(j,k),r2(j,k) ]=make_stats(paramy{delta{k,1}(1)}(iok,end),paramy{delta{k,1}(2)} (iok,end)); 
    89142 
    90      
    91   end 
     143    end %for k 
    92144   
    93145   
    … …  
    102154end 
    103155close(f) 
     156 
     157%sauvegarde du tableau 
     158 
     159fid=fopen('tab.txt','w'); 
     160for j=1:size(sparam,1) 
     161 
     162  fprintf(fid,'\n * ParamÚtre %s \n',sparam{j,1}); 
     163  fprintf(fid,'||Name  ||bias|| RMS || rel. err || r2    ||\n'); 
     164    for k=1:size(delta,1) 
     165       
     166 
     167fprintf(fid,'||%s %s || %4.3e ||  %4.3e ||   %4.3f  || %4.3f ||\n',sparam{j,1},delta{k,2}, MAD(j,k),RMS(j,k), ... 
     168        URPD(j,k),r2(j,k)); 
     169 
     170    end 
     171 
     172    fprintf(fid,'\n'); 
     173 
     174end 
     175fclose(fid); 

  • trunk/scripts/sinobad.i

    r14 r16  
    1616 
    1717xwriteout 0 true ../data_out/exp_T.nc 
    18 savestate sshn_c  1   ij   0    A       1       ../data_out/sshn_c_true.dat 
    19 savestate ta_c  1   ijk   0    A       1       ../data_out/ta_c_true.dat 
    20 savestate sa_c  1   ijk   0    A       1       ../data_out/sa_c_true.dat 
    21 savestate ua_c  1   ijk   0    A       1       ../data_out/ua_c_true.dat 
    22 savestate va_c  1   ijk   0    A       1       ../data_out/va_c_true.dat 
    2318 
    2419print_time OFF 
    … …  
    2924##OBSERVATION 
    3025xwriteout 50 obs_48 ../data_out/exp_T.nc 
    31 savestate sshb  1   ij   50    A       1       ../data_out/sshb_obs_48.dat    
    32 savestate tb  1   ijk   50    A       1       ../data_out/tb_obs_48.dat    
    33 savestate sb  1   ijk   50    A       1       ../data_out/sb_obs_48.dat    
    34 savestate ub  1   ijk   50    A       1       ../data_out/ub_obs_48.dat    
    35 savestate vb  1   ijk   50    A       1       ../data_out/vb_obs_48.dat    
    3626 
    3727loadobs sshb  1   ij   50    A       1       ../data_out/sshb_obs_48.dat D    
    … …  
    4939print_time OFF 
    5040FORWARD 
     41 
    5142xwriteout 0 init ../data_out/exp_T.nc 
    52 savestate sshn_c  1   ij   0    A       1       ../data_out/sshn_c_init.dat 
    53 savestate ta_c  1   ijk   0    A       1       ../data_out/ta_c_init.dat 
    54 savestate sa_c  1   ijk   0    A       1       ../data_out/sa_c_init.dat 
    55 savestate ua_c  1   ijk   0    A       1       ../data_out/ua_c_init.dat 
    56 savestate va_c  1   ijk   0    A       1       ../data_out/va_c_init.dat 
    57  
    5843xwriteout 50 forw0_48 ../data_out/exp_T.nc 
    59 savestate sshb  1   ij   50    A       1       ../data_out/sshb_forw0_48.dat    
    60 savestate tb  1   ijk   50    A       1       ../data_out/tb_forw0_48.dat    
    61 savestate sb  1   ijk   50    A       1       ../data_out/sb_forw0_48.dat    
    62 savestate ub  1   ijk   50    A       1       ../data_out/ub_forw0_48.dat    
    63 savestate vb  1   ijk   50    A       1       ../data_out/vb_forw0_48.dat    
    6444 
    6545##TEST DE LA FONCTION OBJECTIVE 
    … …  
    8363##SAUVEGARDE 
    8464xwriteout 0 fin ../data_out/exp_T.nc 
    85 savestate sshn_c  1   ij   0    A       1       ../data_out/sshn_c_fin.dat 
    86 savestate ta_c  1   ijk   0    A       1       ../data_out/ta_c_fin.dat 
    87 savestate sa_c  1   ijk   0    A       1       ../data_out/sa_c_fin.dat 
    88 savestate ua_c  1   ijk   0    A       1       ../data_out/ua_c_fin.dat 
    89 savestate va_c  1   ijk   0    A       1       ../data_out/va_c_fin.dat 
     65xwriteout 50 forwfin_48 ../data_out/exp_T.nc 
    9066 
    91 xwriteout 50 forwfin_48 ../data_out/exp_T.nc 
    92 savestate sshb  1   ij   50    A       1       ../data_out/sshb_forwfin_48.dat    
    93 savestate tb  1   ijk   50    A       1       ../data_out/tb_forwfin_48.dat    
    94 savestate sb  1   ijk   50    A       1       ../data_out/sb_forwfin_48.dat    
    95 savestate ub  1   ijk   50    A       1       ../data_out/ub_forwfin_48.dat    
    96 savestate vb  1   ijk   50    A       1       ../data_out/vb_forwfin_48.dat    
    97  
    98 goto FINRUN 
    99  
    100 ## analytiquement (2)  
    101 xistate_init 1 ../data_in/file_rest/GYRE_00000424_restart.nc 
    102 #xistate_init 1 ../data_in/file_rest/GYRE_00000400_restart.nc 
     67xrst_save GYRE_00000400_restart_CONTROL_T.nc 
    10368 
    10469 
    105  
    106  
    107 true_target_in_tab ta_c 
    108 savestate ta_c  1   ijk   0    A       1       ../data_out/ta_c_true.dat 
    109 print_time OFF 
    110 set_modeltime 0 
    111 FORWARD 
    112 savestate tb  1   ijk   0    A       3       ../data_out/tb_forw_init_424.dat    
    113 GOTO FINRUN 
    114  
    115 #load_shape_func ta_c ../data_in/shfs/tb400-450/shfs15_N00_tb.dat 
    116 #load_mean 1 ta_c ../data_in/shfs/tb400-450/mean_N00_tb.dat 
    117 #load_stdev_pca ta_c ../data_in/shfs/tb400-450/stdevEOF_N00_tb.dat 
    118  
    119 #loadstate ta_c 1 ijk 0 a 1 ../data_in/start_in_tb/ta_start400.dat D 
    120 #loadstate ta_c  1   ijk 0    A       1   ../data_out/tb_st10.dat D    
    121 #loadstate ta_c 1 ijk 0 a 1    ../data_out/ta_c_yao.dat D 
    122  
    123 #load_mean 0 ta_c | creating the observations  
    124  
    125 #setstate pca_ta 0 
    126 GOTO NOF  
    127 ## Faire une simulation direct  
    128 print_time OFF 
    129 set_modeltime 0 
    130 FORWARD 
    131 #savestate mod_name       ns  ijk   pdt  coding  format  file_name 
    132  
    133 #####savestate tb  1   ijk   0    A       3       ../data_out/tb_50.dat    
    134 #savestate tb  1   ijk   23    A       1       ../data_out/tb_obs20_2.dat    
    135 #savestate sshb  1   ij   23    A       1       ../data_out/sshb_obs20_0_1.dat    
    136  
    137 savestate tb  1   ijk   23    A       1       ../data_out/tb_obs_std_424.dat    
    138 savestate sshb  1   ij   23    A       1       ../data_out/sshb_obs_std_424.dat    
    139  
    140 #savestate tb  1   ijk   10    A       1       ../data_out/tb_st10.dat    
    141 #savestate tb  1   ijk   3    A       3       ../data_out/tb_sbd100.dat    
    142 #savestate sb  1   ijk   3    A       3       ../data_out/sb_sbd100.dat    
    143 #savestate ub  1   ijk   3    A       3       ../data_out/ub_sbd100.dat    
    144 #savestate vb  1   ijk   3    A       3       ../data_out/vb_sbd100.dat    
    145 #savestate sshb  1   ij   3    A       3       ../data_out/sshb_sbd100.dat    
    146 ###savestate ta  1   ijk   3    A       1       ../data_in/start_in_tb/ta_start400.dat    
    147                                                   | sauvegarder la 
    148                                                   | sortie du module tb 
    149                                                   | au pas du temps 100 
    150 #savestate sshb  1  ij   102    A       3       ../data_out/sshb100.dat 
    151                                                   | sauvegarder la 
    152                                                   | sortie du module sshb 
    153                                                   | au pas du temps 100 
    154  
    155 GOTO FINRUN 
    156 NOF 
    157 #loadobs tb  1   ijk   23    A       1       ../data_out/tb_obs20_2.dat D    
    158 loadobs tb  1   ijk   23    A       1       ../data_out/tb_obs_std_16_12_1_424.dat D    
    159 #loadobs sshb  1   ij   23    A       1       ../data_out/sshb_obs20.dat D    
    160 #loadstate ta_c 1 ijk 0 a 1 ../data_in/start_in_tb/GYRE_sortie_mohamed_ptb_198.dat D 
    161 #exit 
    162 print_cost ON 
    163 GOTO 3 
    164  
    165  
    166 GOTO TESTAD          |TESTLT |TESTAD |va vers TESTAD 
    167  
    168 outoobs sshb 1 12 
    169 outoobs tb 1 12 
    170  
    171 true_target_in_tab ta_c 
    172  
    173 GOTO TESTEND | TESTEND |TESTOF 
    174 FINTESTDF 
    175 #TEST DE LA PROGRAMMATION DES DERIVEES 
    176 TESTDF 
    177 testdf 3  3  1  4    f 1.e-10  0.001  10 
    178 stop 
    179 FINTESTDF 
    180  
    181 GOTO TESTEND 
    182 GOTO TESTEND 
    183 #----------------- 
    184 # TEST DE LINEAIRE TANGENT 
    185 TESTLT 
    186  testlt 1.e+00 0.5e-01 2. 15 
    187 stop 
    188 #----------------- 
    189 #TEST DE L'ADJOINT 
    190 TESTAD 
    191 testad 0.001 
    192 stop 
    193 testad 0.001 1.0e-6  1000   1.0e-16 
    194 stop 
    195 #------------------ 
    196 #TEST DE LA FONCTION OBJECTIVE 
    197 TESTOF 
    198 print_cost OFF 
    199 loadstate ta_c 1 ijk 0 a 1 ../data_in/start_in_tb/GYRE_sortie_mohamed_ptb_200.dat D 
    200 testof 1.e+00 0.5e-01 2. 10 1 
    201 TESTEND 
    202  
    203 #--------------------------------------------------------------------- 
    204 #----------------------> choice of the run  <------------------------- 
    205 #--------------------------------------------------------------------- 
    206 #-------> 1. RUNYAO 
    207 #-------> 2. RUNM2QN1 
    208 #-------> 3. RUNM1QN3 
    209  
    210 #charger l'état initial  
    211 loadstate ta_c 1 ijk 0 a 1 ../data_in/start_in_tb/GYRE_sortie_mohamed_ptb_200.dat D 
    212  
    213 GOTO 3 
    214 #--------------------------------------------------------------------- 
    215 1 
    216 !echo "RUN OPTIMIZATION WITH YRUN ......." 
    217 set_nbiter 20 
    218 setepsi ta_c 0.1 
    219 run 
    220 3 
    221 print_time OFF 
    222 set_modeltime 0 
    223 !echo "RUN OPTIMIZATION WITH M1QN3 ......." 
    224 setm_impres  5 
    225 setm_io      6 
    226 setm_mode    1 
    227 set_nbiter   100 
    228 setm_nsim    100 
    229 setm_dxmin   1.0e-4 
    230 setm_epsg    1.0e-10 
    231 setm_ddf1    1 
    232 runm 
    233  
    234 savestate ta_c  1   ijk   0    A       1       ../data_out/ta_c_yao.dat 
    235 savestate tb  1   ijk   23    A       1       ../data_out/tb_ret_424.dat    
    23670 
    23771goto FINRUN 

  • trunk/src/sinobad.h

    r14 r16  
    352352  int itest=atoi(argv[1]); 
    353353  if(itest==1){// from restart file 
    354     printf("start from rest_file "); 
     354    printf("start from rest_file\n "); 
    355355    neuler=0;//=1;// Restart from a file // Set time-step indicator at nit000 (leap-frog) 
    356356    // ifdef YE_*a_c alors *a=*b=*a_c 
Note: See TracChangeset for help on using the changeset viewer.