Dear Ole,
as you suggested I have changed the frequency grid (Q._GRID) to a larger
number and I have confirmed that also the jacobian matrix changed the
dimensions. However, after activating the backend
(Q.SENSOR_RESPONSE.BACKEND_DO = true) I did not achieve to solve the
problem and the next error message came out:
L2 = qpack2( Q, O, Y);
terminate called after throwing an instance of 'std::runtime_error'
what(): Run-time error in method: jacobianClose
Method jacobianClose needs input variable: sensor_response
/opt/arts/arts2.2/build/src/arts -r00 -o
/tmp/atmlab-tempera-tpa6df5bc9_d962_48a4_9b91_5f0713f5df5e
/tmp/atmlab-tempera-tpa6df5bc9_d962_48a4_9b91_5f0713f5df5e/cfile.arts:
Aborted
Regarding the f_backend8f.xml and f_backend.xml the number of elements
are 8 and 12 respectively. I attach again the function that set the
q_structure and also the backend files just in case that you want to
check is there is something wrong.
Once again, thank you very much for your help!
Best regards,
Fran
On 07/20/2015 06:46 PM, Ole Martin Christensen wrote:
I notice that your frequency grid has 14 elements, so another explanation could
be that the backend channels have not been activated in ARTS. You can check
this by changing your frequency grid (Q.F_GRID) to a larger number of
frequencies and see if the size of the jacobian matrix (J) you put into oem
changes.
if this is the case I think you might have to set Q.SENSOR_RESPONSE.BACKEND_DO
= true
to active the backend.
another possiblity is that you have to many elements in in your f_backend8f.xml
and f_backend.xml file? from your error I would almost guess that there are 14
elements in these files.
I hope this helps, or at least diagnoses the problem.
Ole Martin
________________________________________
From: Fran Navas [[email protected]]
Sent: Monday, July 20, 2015 6:09 PM
To: Ole Martin Christensen
Subject: Re: [arts-users] arts 2.2
Dear Ole,
I attach the functions for the setup of the Q and the Y structures and
also the code of the inversion. I hope that it is clear for you and if
you have any question just let me know.
After your comment I am thinking that the problem could come from the
Y.F variable. With ARTS version 2.0 this variable was empty in my
setting and it worked, but now with the last version of ARTS the next
error message appeared: ("Row size of Y(1).TNOISE must be 1 or equal to
the length of Y.F"). For this reason I set Y.F with the frequencies of
our radiometer and with the same dimension than Y.TNOISE, it means that
the frequencies are repeated following the block measurements (12
frequencies*9 elevations angles = 108 values). After this change I did
not get anymore the error message but the one that I explained in my
previous email.
Well, I hope that you can help me with this problem and I finally can
run the retrieval in ARTS 2.2 :) Thank you so much in advance for your help.
Best regards,
Fran
On 07/19/2015 02:58 PM, Ole Martin Christensen wrote:
Hi, Fran
(I'm adding the qpack mailing list here as well)
could you attach your Q file and perhaps the script you use to execute the
retrievals? I assume this problem could be due to some mismatch between the
number of channels specified in you Y.F (qpack Y struckture) and
Q.SENSOR_RESPONSE.
Ole Martin
________________________________________
From: [email protected] [[email protected]] on
behalf of Fran Navas [[email protected]]
Sent: Wednesday, July 15, 2015 10:15 AM
To: ARTS Users List
Subject: [arts-users] arts 2.2
Dear Patrick, Jana and all,
I am having some problems running our tropospheric temperature
retrievals (from filter bank measurements) in ARTS 2.2. These retrievals
were running in arts 2.0 until now. I have adapted the scripts in order
to run it in the last version and now everything look fine until the
different iterations start in the retrieval.
Below you can find the error message. I do not why but the dimensions of
the 'yf' that is calculated in the function: [R,yf,J] = calc_y_and_j(
comfun, Q, R, O, xnorm, x, iter ); are different to the ones from 'y'.
In my retrievals (y: [108*1]; yf: [126*1]). This difference in the
dimensions is the cause of the error. I guess that something is wrong in
the Q-structure but I don't find what can be. I would appreciate a lot
if someone could help me.
Thank you very much in advance!!!
Fran Navas
Error in ==> oem>calc_cost at 945
dd = y - yf;
Error in ==> oem>calc_cost_for_X at 956
[cost,cost_y,cost_x] = calc_cost( y, yf, Seinv, xa, x, Sxinv, O,
xnorm );
Error in ==> oem at 558
X = calc_cost_for_X( y, yf, Seinv, xa, x, Sxinv, O, xnorm, X,
iter );
Error in ==> qpack2 at 116
[X,R] = oem( O, Qm, R, @arts_oem, Sx, Se, Sxinv, [], xa, Y(m).Y );
--
*****************************************************************************
Francisco Navas Guzmán, PhD
Microwave Physics Group
Division of Atmospheric radiometry
Institute of Applied Physics
University of Bern
Address: | http://www.iap.unibe.ch
Sidlerstrasse 5 | Ph: +41-31 631 50 46
CH-3012 Bern (Switzerland) | Fax: +41-31 631 37 65
*****************************************************************************
_______________________________________________
arts-users mailing list
[email protected]
https://www.sat.ltu.se/mailman/listinfo/arts-users
--
*****************************************************************************
Francisco Navas Guzmán, PhD
Microwave Physics Group
Division of Atmospheric radiometry
Institute of Applied Physics
University of Bern
Address: | http://www.iap.unibe.ch
Sidlerstrasse 5 | Ph: +41-31 631 50 46
CH-3012 Bern (Switzerland) | Fax: +41-31 631 37 65
*****************************************************************************
--
*****************************************************************************
Francisco Navas Guzmán, PhD
Microwave Physics Group
Division of Atmospheric radiometry
Institute of Applied Physics
University of Bern
Address: | http://www.iap.unibe.ch
Sidlerstrasse 5 | Ph: +41-31 631 50 46
CH-3012 Bern (Switzerland) | Fax: +41-31 631 37 65
*****************************************************************************
<?xml version="1.0"?>
<arts format="ascii" version="1">
<Vector nelem="12">
5.1250000e+10
5.1750000e+10
5.2250000e+10
5.2850000e+10
5.3350000e+10
5.3850000e+10
5.4400000e+10
5.4900000e+10
5.5400000e+10
5.6000000e+10
5.6500000e+10
5.7000000e+10
</Vector>
</arts>
<?xml version="1.0"?>
<arts format="ascii" version="1">
<Vector nelem="8">
5.3350000e+10
5.3850000e+10
5.4400000e+10
5.4900000e+10
5.5400000e+10
5.6000000e+10
5.6500000e+10
5.7000000e+10
</Vector>
</arts>
%Q for doing temperatiure retrievals using a ground based radiometer.
%Created for Oliver Stähli
%Modified by F. Navas (2014.02.03)
function [Q,O] = setup_Q_tr(Y,TBtimestart,meteo, ILW,payP,paymeanT,paystdT);
cd('/home/tempera/data_processing_v1/retrieval/Inversion_FB_v2_arts22/matlab/aux_func/');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Atmlab settings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
arts_xmldata_path = atmlab( 'ARTS_XMLDATA_PATH' );
arts_includes = atmlab( 'ARTS_INCLUDES' );
if isnan( arts_xmldata_path )
error( 'You need to set ARTS_XMLDATA_PATH to run this file.' );
end
if isnan( arts_includes )
error( 'You need to ARTS_INCLUDES to run this file.' );
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q = qarts;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% General
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.INCLUDES = { fullfile( arts_includes, 'general.arts' ), ...
fullfile( arts_includes, 'continua.arts' ) };
Q.ATMOSPHERE_DIM = 1;
Q.STOKES_DIM = 1;
Q.J_DO = true;
Q.CLOUDBOX_DO = false;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Radiative transfer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Q.Y_UNIT = 'PlanckBT'; % Radiance units (1, RJBT or PlanckBT)
Q.Y_UNIT = 'RJBT';
Q.YCALC_WSMS = { 'yCalc' };
Q.PPATH_LMAX = 2000;
Q.PPATH_STEP_AGENDA = { 'ppath_stepGeometric' };
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Surface
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.R_GEOID = constants( 'EARTH_RADIUS' );
Q.Z_SURFACE = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Absorption
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.ABS_LINES = [];
Q.ABS_LINES_FORMAT = 'None';
%
% The option 'OnTheFly' is more accurate, and here also faster!
if 1
Q.ABSORPTION = 'OnTheFly';
else
Q.ABSORPTION = 'LoadTable';
Q.ABS_LOOKUP = fullfile( pwd, 'abs_lookup.xml' );
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pressure grid - needs to extend below surface, at the surface we use the
% surface pressure from meteo_Payerne
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
z_toa = 100e3;
%
a=p2z_simple(round(meteo.pressure*100))+10;
for ii=1:70
h(ii)=a+ii*20;
a=h(ii);
end
clear a
%Q.P_GRID = z2p_simple( [-250 0 h 26200:1000:z_toa] )';
Q.P_GRID = z2p_simple( [-250 0 h (h(end)+200):1000:z_toa] )';
Q.P_GRID(3,1)=round(meteo.pressure*100);
clear h
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Frequency grid of forward model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.F_GRID = [51.00e9,51.25e9,51.75e9,52.25e9,52.85e9,53.35e9,53.85e9,54.40e9, ...
54.90e9,55.40e9,56.00e9,56.50e9,57.00e9,57.50e9]';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Latitude and longitude
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.LAT_GRID = Y.LATITUDE;
Q.LON_GRID = Y.LONGITUDE;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sensor variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H = qartsSensor;
Q.ANTENNA_DIM = 1;
if ILW > 0.025
H.F_BACKEND = 'f_backend8f.xml';
H.BACKEND_CHANNEL_RESPONSE = 'backend_channel_response8f.xml';
else
H.F_BACKEND = 'f_backend.xml';
H.BACKEND_CHANNEL_RESPONSE = 'backend_channel_response.xml';
end
H.SENSOR_NORM = true;
Q.SENSOR_DO = true;
Q.SENSOR_RESPONSE = H;
% Defining the antenna pattern
antenna_pattern = dlmread('antenna.aa');
Q.SENSOR_RESPONSE.ANTENNA_RESPONSE.name = 'Antenna response';
Q.SENSOR_RESPONSE.ANTENNA_RESPONSE.gridnames = {'Polarisation','Frequency','Zenith angle','Azimut angle'};%
Q.SENSOR_RESPONSE.ANTENNA_RESPONSE.grids = {{'1'},5.385e10,antenna_pattern(:,1),0};
Q.SENSOR_RESPONSE.ANTENNA_RESPONSE.dataname = 'Response';
Q.SENSOR_RESPONSE.ANTENNA_RESPONSE.data(1,1,:,1) = antenna_pattern(:,2);
Q.SENSOR_RESPONSE.ANTENNA_DO = 1;
Q.SENSOR_RESPONSE.ANTENNA_LOS = [30:5:70]';
Q.MBLOCK_ZA_GRID = [20:2:80]';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set the zenith angle grid for each measurement block.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if TBtimestart>datenum(datestr('2012-06-05 08:00:00',31))&TBtimestart<...
% datenum(datestr('2012-06-29 12:10:00',31))
% Q.MBLOCK_ZA_GRID = [34.5:5:74.5]';
% else
% Q.MBLOCK_ZA_GRID = [30:5:70]';
% end
%- Thermal noise
f = qarts_get(H.F_BACKEND);
%Q.TNOISE_C = covmat1d_from_cfun(repmat(f,length(Q.MBLOCK_ZA_GRID),1),[],'drc');
Q.TNOISE_C = covmat1d_from_cfun(repmat(f,length(Q.SENSOR_RESPONSE.ANTENNA_LOS),1),[],'drc');
clear H f
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define L2 structure (beside retrieval quantities)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.L2_EXTRA = { 'cost', 'dx', 'xa', 'y', 'yf', 'bl', 'ptz', 'mresp', ...
'e', 'eo', 'es', 'date','A','tnoise' };
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Temperature Retrieval and HSE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.T.RETRIEVE = true;
Q.T.L2 = true;
Q.T.GRIDS = { Q.P_GRID(3:1:60), [], [] };
Q.T.ATMDATA = gf_artsxml( fullfile( arts_xmldata_path, 'climatology', ...
'cira86', 'cira86.t.xml' ), 'Temperature', 't_field' );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% apriori covariance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sizeTapriori=size(Q.T.GRIDS{1});
if ILW > 0.025
interpTstdpaytest=interp1([1 4 8 10 12 16 18 19 21 32 sizeTapriori(1,1)],...
[2 1.8 1.5 1.45 1.35 1.25 1.23 1.22 1.25 1.45 1.5],[1:1:sizeTapriori]);%%Stand: 2012.10.15
Q.T.SX = covmat1d_from_cfun( Q.T.GRIDS{1}, [Q.T.GRIDS{1} interpTstdpaytest'], ...
'exp', 0.14, .03, @log10 );
else
interpTstdpaytest=interp1([1 16 sizeTapriori(1,1)],[2 1.7 1.5],[1:1:sizeTapriori]);%%Stand: 2012.10.15 (apriori errors)
Q.T.SX = covmat1d_from_cfun( Q.T.GRIDS{1}, [Q.T.GRIDS{1} interpTstdpaytest'], ...
'exp', 0.2, .08, @log10 );
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% new a-priori covariance matrix from climatology (radiosondes)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%delete the rows with repeated values of pressure
% payP2(1,1)=payP(1);
% paystdT2(1,1)=paystdT(1);
% j=2;
% for i=2:size(payP,2)
% if payP(i)~=payP(i-1)
% payP2(j,1)=payP(i);
% paystdT2(j,1)=paystdT(i);
% j=j+1;
% end
% end
%
% Q.T.SX = covmat1d_from_cfun(Q.T.GRIDS{1}, [payP2 paystdT2],'exp', 0.2, .08, @log10 );
%
% Q.T.HSE = 'on';
% Q.T.METHOD = 'analytical';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine altitudes through HSE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.HSE.ON = true;
Q.HSE.ACCURACY = 0.1;
Q.HSE.P = 1013;
if Q.HSE.ON == false
Q.Z_ATMDATA = gf_artsxml( fullfile( arts_xmldata_path, 'climatology', ...
'cira86', 'cira86.z.xml' ), 'Altitudes', 'z_field' );
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Absorption in Troposphere and Stratosphere
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Oxygen
Q.ABS_SPECIES(1).TAG = {'O2-PWR93'};
Q.ABS_SPECIES(1).RETRIEVE = false;
% Water
Q.ABS_SPECIES(2).TAG = {'H2O-PWR98'};
Q.ABS_SPECIES(2).RETRIEVE = false;
% Nitrogen
Q.ABS_SPECIES(3).TAG = {'N2-SelfContStandardType'};
Q.ABS_SPECIES(3).RETRIEVE = false;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use Standardatmosphere separate for summer and winter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if Y.MONTH>4 && Y.MONTH<10
season = 'summer';
else
season = 'winter';
end
Q.ABS_SPECIES(1).ATMDATA = gf_artsxml(fullfile( arts_xmldata_path, ...
'atmosphere/fascod', ...
sprintf('midlatitude-%s.O2.xml',season)),'O2','vmr_field');
Q.ABS_SPECIES(2).ATMDATA = gf_artsxml(fullfile( arts_xmldata_path, ...
'atmosphere/fascod', ...
sprintf('midlatitude-%s.H2O.xml',season)),'H2O','vmr_field');
Q.ABS_SPECIES(3).ATMDATA = gf_artsxml(fullfile( arts_xmldata_path, ...
'atmosphere/fascod', ...
sprintf('midlatitude-%s.N2.xml',season)),'N2','vmr_field');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OEM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
O = qp2_l2( Q );
%
O.linear = false;
%
if ~O.linear
O.itermethod = 'GN';
%O.itermethod = 'ML';
%O.stop_dx = 0.001;
O.stop_dx = 0.01;
%O.maxiter = 10;
O.maxiter = 10;
O.ga_start = 100;
O.ga_max = 1e7;
O.ga_factor_ok = 10;
O.ga_factor_not_ok = 10;
end
_______________________________________________
qpack mailing list
[email protected]
https://www.sat.ltu.se/mailman/listinfo/qpack
