Current reconstruction using Bayesian inverse filter.
[syntax]
[Zact,Jinfo,bayes_parm,vb_parm,MEGinfo,Jext,Pointlist] ...
= vbmeg_current_reconstruct_z2(proj_root, curr_parm)
[input]
proj_root: <<string>> VBMEG project root directory.
curr_parm: <<struct>> Parameters for current estimation.
--- fields of curr_parm
dbayesfile : <<string>> Dynamics estimation file (.dbayes.mat).
currfile : <<string>> Cortical current file (.curr.mat), created
by this function.
trial_average: <optional> <<bool>> If true,
= [ON] : average current over all sessions
= OFF : current for each session
ix_area : <optional> Vertex indices to calculate estimated
current. If 'ix_area' is empty or not given, cortical
currents in the active region are calculated.
tsubsmpl : <optional> <bosolete> Specify subsampled time
index. If 'tsubsmpl' is empty or not given, time
subsampling is not done.
dsampf : <optional> <<int>> Specify frequency of
downsampling. This value must be smaller than the
original sampling frequency of M/EEG data.
overlap_mode : <optional> <<int>>
= 0 : current is averaged over overlapped time window
= 1 : current is not averaged for overlapped window
current time series of each time windows
are concatenated sequentially for spectral analysis
verbose : <<bool>> Verbose flag.
---
[note] !! currently unconfirmed.
If following field is given, these values are used instead of
bayes_parm field in result file:
---
curr_parm.basisfile
curr_parm.megfile
curr_parm.twin_meg
curr_parm.Tperiod
curr_parm.Tnext
curr_parm.extra.basisfile (for extra dipole)
---
[output]
Zact : active current
Zact(n,t,:) is the current at the vertex 'ix_act(n)' & the time 't'
Zact(Nact,Nsample) for trial_average = ON
Zact(Nact,Nsample,Ntrials) for trial_average = OFF
Nact : # of active region,
Nsample : # of time sample,
Ntrials : # of trials in all session]
--- fields of Jinfo
version : <<string>> Version of cortical current file.
curr_type : <<string>> 'Z-current'. It can be 'J-current' for VBMEG
version 0.8 or older.
Wact : <<float matrix>> Smoothing Gaussian filter, mapping from
Z-current to J-current.
ix_act : <<int vector>>: Vertex indices of Z-current.
ix_act_ex : <<int vector>>: Vertex indices of J-current.
Lact : <<int>> Number of current direction at one vertex.
Tsample : <<int vector>> Time sample indices of the original MEG
data. length(Tsample) == size(Zact,2) == size(Jact,2).
Tmsec : <<float vector>> Time in msec.
SampleFreq: <<float>> Sample frequency of cortical current, not
original M/EEG signal [Hz].
Pretrigger: <<int>> Time points of the length of the pretrigger
period of cortical current data. It is neither actual time
nor time points of the original M/EEG signal.
Ntrial : <<int>> Number of trials of estimated current.
patch_norm: <<bool>> Cortical current is patch size normalized
(Am/m^2) or not (Am).
Tix : <<L x 1 cell>> Time sample indices of each time window.
Zact(:,Tix{n},:) is the set of Z-current within the n-th
time window.
---
* Inverse filter calculation is done in vb_invfilter_calc
[history]
2006-09-03 M. Sato
* Non-overlapped concatenation mode is added for spectral analysis
2008-08-19 Taku Yoshioka
Extra dipole support
2008-09-30 Taku Yoshioka
Minor change for variables in Jinfo
2008-10-23 Taku Yoshioka
Bug fix for current estimation without extra dipoles
2009-04-02 Taku Yoshioka
Parameter name changed within this code for readability
(just replacing 'resultfile' to bayesfile)
2010-03-01 M. Sato
Bug fix for Wact index and fieldname(tsubsamp->tsubsmpl)
2010-12-06 taku-y
[enhancement] curr_parm.dsampf supported.
[minor] Following fields of Jinfo set in this function:
SampleFreq
Pretrigger
Tmsec
[trivial] Jinfo.version = vb_version.
2010-12-07 taku-y
[trivial] Jinfo.version = vbmeg('version');
2011-05-11 takiu-y
[debug] Jinfo.Tmsec corrected.
2011-06-28 taku-y
[minor] Jinfo.Tix added.
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [Zact,Jinfo,bayes_parm,vb_parm,MEGinfo,Jext,Pointlist] ... 0002 = vb_current_reconstruct_z_dynamics(proj_root, curr_parm) 0003 % Current reconstruction using Bayesian inverse filter. 0004 % 0005 % [syntax] 0006 % [Zact,Jinfo,bayes_parm,vb_parm,MEGinfo,Jext,Pointlist] ... 0007 % = vbmeg_current_reconstruct_z2(proj_root, curr_parm) 0008 % 0009 % [input] 0010 % proj_root: <<string>> VBMEG project root directory. 0011 % curr_parm: <<struct>> Parameters for current estimation. 0012 % --- fields of curr_parm 0013 % dbayesfile : <<string>> Dynamics estimation file (.dbayes.mat). 0014 % currfile : <<string>> Cortical current file (.curr.mat), created 0015 % by this function. 0016 % trial_average: <optional> <<bool>> If true, 0017 % = [ON] : average current over all sessions 0018 % = OFF : current for each session 0019 % ix_area : <optional> Vertex indices to calculate estimated 0020 % current. If 'ix_area' is empty or not given, cortical 0021 % currents in the active region are calculated. 0022 % tsubsmpl : <optional> <bosolete> Specify subsampled time 0023 % index. If 'tsubsmpl' is empty or not given, time 0024 % subsampling is not done. 0025 % dsampf : <optional> <<int>> Specify frequency of 0026 % downsampling. This value must be smaller than the 0027 % original sampling frequency of M/EEG data. 0028 % overlap_mode : <optional> <<int>> 0029 % = 0 : current is averaged over overlapped time window 0030 % = 1 : current is not averaged for overlapped window 0031 % current time series of each time windows 0032 % are concatenated sequentially for spectral analysis 0033 % verbose : <<bool>> Verbose flag. 0034 % --- 0035 % 0036 % [note] !! currently unconfirmed. 0037 % If following field is given, these values are used instead of 0038 % bayes_parm field in result file: 0039 % --- 0040 % curr_parm.basisfile 0041 % curr_parm.megfile 0042 % curr_parm.twin_meg 0043 % curr_parm.Tperiod 0044 % curr_parm.Tnext 0045 % curr_parm.extra.basisfile (for extra dipole) 0046 % --- 0047 % 0048 % [output] 0049 % Zact : active current 0050 % 0051 % Zact(n,t,:) is the current at the vertex 'ix_act(n)' & the time 't' 0052 % Zact(Nact,Nsample) for trial_average = ON 0053 % Zact(Nact,Nsample,Ntrials) for trial_average = OFF 0054 % Nact : # of active region, 0055 % Nsample : # of time sample, 0056 % Ntrials : # of trials in all session] 0057 % --- fields of Jinfo 0058 % version : <<string>> Version of cortical current file. 0059 % curr_type : <<string>> 'Z-current'. It can be 'J-current' for VBMEG 0060 % version 0.8 or older. 0061 % Wact : <<float matrix>> Smoothing Gaussian filter, mapping from 0062 % Z-current to J-current. 0063 % ix_act : <<int vector>>: Vertex indices of Z-current. 0064 % ix_act_ex : <<int vector>>: Vertex indices of J-current. 0065 % Lact : <<int>> Number of current direction at one vertex. 0066 % Tsample : <<int vector>> Time sample indices of the original MEG 0067 % data. length(Tsample) == size(Zact,2) == size(Jact,2). 0068 % Tmsec : <<float vector>> Time in msec. 0069 % SampleFreq: <<float>> Sample frequency of cortical current, not 0070 % original M/EEG signal [Hz]. 0071 % Pretrigger: <<int>> Time points of the length of the pretrigger 0072 % period of cortical current data. It is neither actual time 0073 % nor time points of the original M/EEG signal. 0074 % Ntrial : <<int>> Number of trials of estimated current. 0075 % patch_norm: <<bool>> Cortical current is patch size normalized 0076 % (Am/m^2) or not (Am). 0077 % Tix : <<L x 1 cell>> Time sample indices of each time window. 0078 % Zact(:,Tix{n},:) is the set of Z-current within the n-th 0079 % time window. 0080 % --- 0081 % 0082 % * Inverse filter calculation is done in vb_invfilter_calc 0083 % 0084 % [history] 0085 % 2006-09-03 M. Sato 0086 % * Non-overlapped concatenation mode is added for spectral analysis 0087 % 2008-08-19 Taku Yoshioka 0088 % Extra dipole support 0089 % 2008-09-30 Taku Yoshioka 0090 % Minor change for variables in Jinfo 0091 % 2008-10-23 Taku Yoshioka 0092 % Bug fix for current estimation without extra dipoles 0093 % 2009-04-02 Taku Yoshioka 0094 % Parameter name changed within this code for readability 0095 % (just replacing 'resultfile' to bayesfile) 0096 % 2010-03-01 M. Sato 0097 % Bug fix for Wact index and fieldname(tsubsamp->tsubsmpl) 0098 % 2010-12-06 taku-y 0099 % [enhancement] curr_parm.dsampf supported. 0100 % [minor] Following fields of Jinfo set in this function: 0101 % SampleFreq 0102 % Pretrigger 0103 % Tmsec 0104 % [trivial] Jinfo.version = vb_version. 0105 % 2010-12-07 taku-y 0106 % [trivial] Jinfo.version = vbmeg('version'); 0107 % 2011-05-11 takiu-y 0108 % [debug] Jinfo.Tmsec corrected. 0109 % 2011-06-28 taku-y 0110 % [minor] Jinfo.Tix added. 0111 % 0112 % Copyright (C) 2011, ATR All Rights Reserved. 0113 % License : New BSD License(see VBMEG_LICENSE.txt) 0114 0115 if ~isempty(proj_root) 0116 dbayesfile = fullfile(proj_root, curr_parm.dbayesfile); 0117 else 0118 dbayesfile = curr_parm.bayesfile; 0119 end 0120 0121 0122 % 0123 % Verbose level setting 0124 % (note: 'verbose_level' is not related to input variable 0125 % 'curr_parm.verbose'. So two configurations relating to message display 0126 % are coexisting in this function.) 0127 % 0128 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0129 global vbmeg_inst 0130 verbose_const = vb_define_verbose; 0131 0132 if isempty(vbmeg_inst) | ~isfield(vbmeg_inst,'verbose_level'), 0133 verbose_level = verbose_const.VERBOSE_LEVEL_NOTICE; 0134 else 0135 verbose_level = vbmeg_inst.verbose_level; 0136 end 0137 0138 % 0139 % load VBMEG estimated result 0140 % 0141 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0142 load(dbayesfile, 'Model','vb_parm','dbayes_parm'); 0143 bayes_file = fullfile(proj_root, dbayes_parm.bayes_file); 0144 load(bayes_file, 'bayes_parm', 'Model_ext', 'Pointlist'); 0145 0146 % 0147 % check parameter of 'curr_parm' 0148 % 0149 % Values of 'curr_parm' fields dominates over 0150 % those of 'bayes_parm' in bayesfile 0151 % 0152 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0153 [bayes_parm,ix_area,trial_average,tsubsamp,overlap_mode,dsampf] ... 0154 = check_arg(bayes_parm, curr_parm); 0155 0156 if ~isempty(proj_root); 0157 bayes_parm_abs = vb_parm_absolute_path(proj_root, bayes_parm); 0158 else 0159 bayes_parm_abs = bayes_parm; 0160 end 0161 0162 % create temp area file and add diffusion parcellation area 0163 % to the copied area file. then use it instead of original area file. 0164 area_key = 'dmri_parcellation_area'; 0165 [bayes_parm_abs, dbayes_parm] = ... 0166 inner_replace_area_file_with_diffusion_area_file(... 0167 bayes_parm_abs, dbayes_parm, area_key); 0168 0169 % 0170 % MEG data preparation 0171 % B : MEG data 0172 % 0173 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0174 [B,Ntrials,Nch,Tsample,Twindow,Tmsec] ... 0175 = vb_megdata_preparation(bayes_parm_abs); 0176 MEGinfo = vb_load_meg_info(bayes_parm_abs.megfile{1}); 0177 0178 % 0179 % Preparation of lead fields 0180 % Gact : leadfield of focal window 0181 % ix_act : Vertex index corresponding to active current Zact 0182 % ix_act_ex : Vertex index corresponding to active current Jact 0183 % Wact : Spatial smoothing matrix of focal window 0184 % 0185 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0186 0187 % Focal window 0188 vb_disp('--- Lead field matrix of focal window'); 0189 0190 lf_parm.brainfile = bayes_parm_abs.brainfile; 0191 lf_parm.areafile = bayes_parm_abs.areafile; 0192 lf_parm.patch_norm = bayes_parm_abs.patch_norm; 0193 lf_parm.expand_spatial_filter = bayes_parm_abs.expand_spatial_filter; 0194 lf_parm.spatial_smoother = bayes_parm_abs.spatial_smoother; 0195 lf_parm.basisfile = bayes_parm_abs.basisfile; 0196 lf_parm.area_key = bayes_parm_abs.area_key; 0197 lf_parm.reduce = bayes_parm_abs.reduce; 0198 lf_parm.Rfilt = bayes_parm_abs.Rfilt; 0199 lf_parm.remove_area_key = []; 0200 0201 [Gact, ix_act, ix_act_ex, Wact, Lact] = ... 0202 vb_leadfield_preparation(lf_parm); 0203 0204 % Extra dipole 0205 if isfield(bayes_parm_abs,'extra') & ~isempty(bayes_parm_abs.extra), 0206 vb_struct2vars(bayes_parm_abs,{'extra'}); 0207 fprintf('--- Lead field matrix of extra dipoles \n'); 0208 for n=1:length(extra.basisfile) 0209 tmp = vb_load_basis(extra.basisfile{n}); 0210 Gext{n} = tmp'; 0211 end 0212 else 0213 Gext = []; 0214 end 0215 0216 % 0217 % Area index in which current is calculated 0218 % 0219 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0220 if ~isempty(ix_area), 0221 % Select vertex index 'ix_area' within the active current region 0222 [jx_area_ex, ix_area_ex] = vb_index2indexsub(ix_area, ix_act_ex); 0223 else 0224 jx_area_ex = 1:length(ix_act_ex); 0225 end 0226 0227 Wact = Wact(jx_area_ex,:); 0228 jx_act = find( sum(Wact, 1) > 0); 0229 Wact = Wact(:,jx_act); 0230 0231 % active index of Z-current 0232 ix_act = ix_act(jx_act); 0233 % active index of J-current 0234 ix_act_ex = ix_act_ex(jx_area_ex); 0235 0236 % # of active vertex 0237 Njact_area = length(jx_act); 0238 0239 % # of extra dipoles 0240 if ~isempty(Gext), Njext = size(Gext{1},2); 0241 else Njext = 0; end 0242 0243 % 0244 % Constant 0245 % 0246 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0247 Nsession = length(B); % Number of sessions 0248 Ntotal = sum(Ntrials); % Total number of trials in all sessions 0249 0250 % 0251 % Temporal subsampling index 0252 % 0253 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0254 sampf = MEGinfo.SampleFreq; 0255 0256 if ~isempty(dsampf), 0257 tsubsamp = ceil(1:sampf/dsampf:Tsample); 0258 Jinfo.SampleFreq = dsampf; 0259 else 0260 if isempty(tsubsamp), tsubsamp = 1:Tsample; end 0261 Jinfo.SampleFreq = sampf; 0262 end 0263 0264 %Jinfo.Tmsec = Tmsec(tsubsamp); 0265 [tmp,ix] = min(abs(Tmsec(tsubsamp))); 0266 Jinfo.Pretrigger = ix; 0267 0268 % 0269 % Temporal smoothing window weight 0270 % 0271 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0272 [Tweight ,Tindex, Nindex, Nsample] = ... 0273 vb_overlapped_timewindow_weight(Twindow, Tsample, tsubsamp, overlap_mode); 0274 0275 Nwindow = length(Nindex); % # of time window 0276 Jinfo.Tix = Nindex; 0277 0278 if overlap_mode == 1, 0279 vb_disp('Non-overlapped concatenation mode'); 0280 end 0281 0282 % 0283 % Initialization 0284 % 0285 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0286 if trial_average == ON 0287 % Current averaged over trials 0288 Zact = zeros(Njact_area,Nsample); 0289 Jext = zeros(Njext,Nsample); 0290 else 0291 % Current for each trial 0292 Zact = zeros(Njact_area,Nsample,Ntotal); 0293 Jext = zeros(Njext,Nsample,Ntotal); 0294 end 0295 0296 % 0297 % Estimated current variance 0298 % 0299 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0300 a_inv = Model.a; 0301 Cov = Model.Cov; 0302 if ~isempty(Model_ext), e_inv = Model_ext.a; 0303 else e_inv = []; end 0304 0305 % 0306 % Store current sources in bayes.mat 0307 % 0308 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0309 if trial_average == ON 0310 if Nsession == 1 0311 if size(Model.Z{1},3) == 1 0312 Zact = Model.Z{1}; 0313 Jext = []; 0314 else 0315 error('The number of trials must be one.'); 0316 end 0317 else 0318 error('The number of sessions must be one.'); 0319 end 0320 else 0321 error('Trials must be averaged.'); 0322 end 0323 0324 % 0325 % Current Info 0326 % 0327 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0328 Jinfo.version = '2.0-0.b.9'; %vbmeg('version'); 0329 Jinfo.curr_type = 'Z-current'; 0330 0331 Jinfo.Lact = Lact; 0332 Jinfo.Wact = Wact; 0333 Jinfo.ix_act = ix_act; 0334 Jinfo.ix_act_ex = ix_act_ex; 0335 Jinfo.NJact = Njact_area; 0336 Jinfo.Nsession = Nsession; 0337 Jinfo.jactdir = []; 0338 Jinfo.Ntrial = Ntrials; 0339 0340 % ix_act : Vertex index corresponding to active current Zact 0341 % ix_act_ex : Vertex index corresponding to active current Jact 0342 % Wact : Spatial smoothing matrix of focal window 0343 % Jact = Wact * Zact 0344 0345 % Actual time corresponding to columns of Zact, supporting overlap mode 0346 % and non-overlapped concatenation mode 0347 Tid_all = []; 0348 Nid_all = []; 0349 for j=1:Nwindow 0350 Tid_all = [Tid_all Tindex{j}]; 0351 Nid_all = [Nid_all Nindex{j}]; 0352 end 0353 0354 if overlap_mode==false, 0355 ix = unique(Tid_all); 0356 Jinfo.Tmsec = Tmsec(ix); 0357 else 0358 Jinfo.Tmsec = Tmsec(Tid_all); 0359 end 0360 0361 %Tstart = bayes_parm.twin_meg(1); 0362 %Tend = bayes_parm.twin_meg(2); 0363 %if isempty(tsubsamp) 0364 % Jinfo.Tsample = Tstart:Tend; 0365 %else 0366 % Jinfo.Tsample = tsubsamp + Tstart - 1; 0367 %end 0368 0369 return; 0370 0371 %%%% --------------- 0372 function [bayes_parm, ix_area, trial_average, tsubsamp, overlap_mode,dsampf] = ... 0373 check_arg(bayes_parm,curr_parm) 0374 0375 if isfield(curr_parm,'basisfile'), 0376 bayes_parm.basisfile = curr_parm.basisfile; 0377 end; 0378 if isfield(curr_parm,'megfile'), 0379 bayes_parm.megfile = curr_parm.megfile ; 0380 end; 0381 if isfield(curr_parm,'twin_meg'), 0382 bayes_parm.twin_meg = curr_parm.twin_meg ; 0383 end; 0384 if isfield(curr_parm,'Tperiod'), 0385 bayes_parm.Tperiod = curr_parm.Tperiod ; 0386 end; 0387 if isfield(curr_parm,'Tnext'), 0388 bayes_parm.Tnext = curr_parm.Tnext ; 0389 end; 0390 0391 if ~isfield(curr_parm,'trial_average'), 0392 trial_average = ON; 0393 else 0394 trial_average = curr_parm.trial_average; 0395 end; 0396 0397 bayes_parm.trial_average = trial_average; 0398 0399 if ~isfield(curr_parm,'ix_area'), 0400 ix_area = []; 0401 else 0402 ix_area = curr_parm.ix_area; 0403 end; 0404 if ~isfield(curr_parm,'tsubsmpl'), 0405 tsubsamp = []; 0406 else 0407 tsubsamp = curr_parm.tsubsmpl; 0408 end; 0409 if ~isfield(curr_parm,'overlap_mode'), 0410 overlap_mode = 0; 0411 else 0412 overlap_mode = curr_parm.overlap_mode; 0413 end; 0414 0415 if isfield(curr_parm,'extra'), 0416 if isfield(curr_parm.extra,'basisfile'), 0417 bayes_parm.extra.basisfile = curr_parm.extra.basisfile; 0418 end; 0419 end; 0420 0421 if ~isfield(curr_parm,'dsampf'), 0422 dsampf = []; 0423 else 0424 dsampf = curr_parm.dsampf; 0425 end; 0426 0427 return; 0428 0429 function [bayes_parm_abs, dbayes_parm] = ... 0430 inner_replace_area_file_with_diffusion_area_file(... 0431 bayes_parm_abs, ... 0432 dbayes_parm, area_key) 0433 % copy specified area file as a temporary file and dMRI are into the file, 0434 % then return the filename. 0435 % [Input] 0436 % bayes_parm_abs : initial value is calculated by this bayes_parm. 0437 % dnayes_parm : dynamics bayes estimation parameter. 0438 % area_key : save dmri area into a temporary file by this name. 0439 % [Output] 0440 % bayes_parm_abs.areafile : temporary area filename including dMRI area. 0441 0442 if ~isfield(dbayes_parm, 'area_key') || isempty(dbayes_parm.area_key) 0443 dbayes_parm.area_key = area_key; 0444 end 0445 if ~isfield(dbayes_parm, 'area_key_global') || isempty(dbayes_parm.area_key_global) 0446 dbayes_parm.area_key_global = area_key; 0447 end 0448 0449 % find area key 0450 found = false; 0451 load(bayes_parm_abs.areafile, 'Area'); 0452 for k=1:length(Area) 0453 if strcmp(Area{k}.key, area_key) 0454 found = true; 0455 end 0456 end 0457 if found, return; end 0458 0459 % bayes_parm_abs.areafile 0460 load(dbayes_parm.dmri_file, 'vbmeg_area_ix'); 0461 0462 % copy area file to temp area file and put dMRI area into the file. 0463 temp_area_file = [tempname(tempdir), '.area.mat']; 0464 copyfile(bayes_parm_abs.areafile, temp_area_file); 0465 area_new.Iextract = vbmeg_area_ix; 0466 area_new.key = area_key; 0467 vb_add_area(temp_area_file, area_new); 0468 0469 % use temp_area_file instead of original area file. 0470 bayes_parm_abs.areafile = temp_area_file;