% Compute minimum norm solutions on differenct shells of 3D source space, output for different shells separately, output only for specified locations possible

% lfdmat1: Leadfield matrix for given source space and sensor/electrode configuration
% dim: Number of dipole components, e.g. 3 for EEG (x,y,z), 2 for MEG (tangential)
% regpar: Regularization parameter for Tikhonov-Regularization (e.g. chosen by adjusting residual variance)
%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.                                                                      
% paths, names: list of paths and names for datafiles to be analyzed, if paths,names=='': File list can be created via a menu (choose files, then press 'Abbrechen')
% format: input data format, 'avr' for BESA-*.avr-files, 'mfxtxt' for BTI-MFX-ASCII-output (assuming first row elements to be latency values!)
% pathout: Pathname for output of results
% elplist: Channel list with channels to be excluded from calculations. If lfdmat1 has as many rows as there are remaining channels, it is used as it is, otherwise
% the rows specified by elplist are also excluded from lfdmat1 (and lfdmat1 is average referenced for EEG) and output as lfdmat
% locations: 3*nr_points-matrix with locations for which minimum norm values shall be output (e.g. electrode positions) for each shell
% from, to: Range of columns of data matrices to be analyzed, if from,to=='': All columns are included
% out_flag: 'abs': for every location, the modulo of dipole strength is output (only positive values),  'rad': only radial dipole component (EEG) is output (pos. and neg. values) 

function [diploc, lfdmat, data, G, paths, names] = invmap_reduce_find_regpar(lfdmat1, dim, diploc, regpar, paths, names, format, pathout, elplist, locations, from, to, out_flag, extension);

if nargin==0,
   disp(' [diploc, lfdmat, data, G, paths, names] = invmap_reduce_find_regpar(lfdmat1, dim, diploc, regpar, paths, names, format, pathout, elplist, locations, from, to, out_flag, (extension)); ');
   return;
end;
elplist=elplist
locations
from
to

if nargin==13,
   extension = '';
end;


chars = '';  	% Number of characters from input filename to be used as output filename, default ('') = all without extension
   
if isempty(names),
   if strcmp(format, 'avr')|strcmp(format, 'avrnoref'), ending = '*.avr'; end;
	if strcmp(format, 'mfxtxt'), ending = '*.asc'; end;
   [namelist,names,paths] = read_filelist(ending);
   nr_names = length(namelist(:,1));
else
   nr_names = length(names(:,1));
   for i=1:nr_names,
      namelist(i,1:length(deblank(names(i,:)))+length(deblank(paths(i,:)))) = sprintf('%s%s', deblank(paths(i,:)), deblank(names(i,:)));
   end;
end;
namelist

[data, m, n, latencies] = read_data(namelist, elplist, format, from, to);
nr_names = length(namelist(:,1));
nr_electrodes = m(1);
text = sprintf('Number of electrodes to be included in the following procedure: %d', nr_electrodes);
disp(text);

if isempty(diploc),			% If necessary, read dipole locations
	disp('');
	disp(' Reading dipole locations (diploc) ');
	disp('');
	diploc = read_matrix(655,3,'dipoles_symmetric2.dat')
	%diploc = read_matrix(655, 4,'loc_08_02_15.sph');
	diploc = diploc';
   size(diploc)
end;
if isempty(lfdmat1),
   disp('');
   disp(' Reading leadfield matrix (lfdmat) ');
   disp('');
	lfdmat1 = read_matrix(1965,148,'lfd_mat_148_1965');
	lfdmat1 = lfdmat1';
	disp(size(lfdmat1));
end
if nr_electrodes~=length(lfdmat1(:,1)),
   disp('');
   disp(' Resizing leadfield matrix ');
   disp('');
   lfdmat = reduce_columns(elplist, lfdmat1')';
   if strcmp(format, 'avr'), 
      disp('');
      disp(' Reaveragereferencing leadfield matrix! ');
      disp('');
      lfdmat = avg_ref(lfdmat); 
   end;
else
   lfdmat = lfdmat1;
end;


if length(regpar)>1,		% If interval is given, reg.curve is plotted and regpar is entered by hand
   disp('Plotting regularisation curve!');
   regularisation_curve(lfdmat, regpar, data, n);
   disp('Choose regularisation parameter!');
   answer = inputdlg('Choose regularisation parameter')
   regpar = str2num(char(answer));
   disp(regpar);
end;



disp('');
size(lfdmat)
disp(' Computing pseudoinverse (G) ');
disp('');
G = pinv_tikh(lfdmat,regpar);
disp(size(G));

disp('');
disp(' Computing residual variances (0-1) ');
disp('');
nsum = 0;
sumvar = 0;
for i=1:nr_names,
   inv = G*data(1:m(i),1:n(i),i);
   variances = res_var(lfdmat, inv, data(1:m(i),1:n(i),i))';
   minvar(i) = min(variances);
   maxvar(i) = max(variances);
   sumvar(i) = sum(variances);
   nsum = nsum+n(i);
   text = sprintf('Min: %f   Max: %f   Mean: %f', minvar(i), maxvar(i), sumvar(i)/n(i));
   tmp_name = deblank(names(i,:));
   tmp_name = [pathout,tmp_name(1:length(tmp_name)-4) extension '.var'];
   write_matrix([latencies(i,1:n(i))', variances(1:n(i))], tmp_name);
   disp(tmp_name);
   clear inv;
end;
text = sprintf('Over all:    Min: %f   Max: %f   Mean: %f', min(minvar), max(maxvar), sum(sumvar)/nsum); 
disp(text);


[q, r] = when_changes_radius(diploc(1:3,:), 0.001);

% Finde Quellenorte nahe den angegebenen Positionen 
if ~isempty(locations),
	for i=1:length(q),
   	for j=1:length(locations(1,:)),
     	if i==1,		% oberste Schicht
            diff(1,:) = diploc(1,1:q(1)) - locations(1,j);   
        	diff(2,:) = diploc(2,1:q(1)) - locations(2,j);
	        diff(3,:) = diploc(3,1:q(1)) - locations(3,j);
   	        offset = 0;			% Anzahl der Punkte in oberen Schichten
     	else        
 	  		diff(1,:) = diploc(1,q(i-1)+1:q(i)) - locations(1,j);   
        	diff(2,:) = diploc(2,q(i-1)+1:q(i)) - locations(2,j);
	        diff(3,:) = diploc(3,q(i-1)+1:q(i)) - locations(3,j);
   	        offset = q(i-1);  
	    end;       
	 	diff = norm_col(diff);
		[Y,I] = min(diff);
	   	elem(i,j) = offset+I;				% Indizes der nahesten Quellenraumelemente
	   	if strcmp(out_flag, 'abs')|strcmp(out_flag, 'ori')|strcmp(out_flag, 'oriabs')|strcmp(out_flag, 'absori'),
	        for k=1:dim,
	        	elem_G(i,(j-1)*dim+k) = dim*elem(i,j)-dim+k;
	        end
	    end
	    if strcmp(out_flag, 'rad'),
	      	elem_G(i,j) = elem(i,j)*dim;
	    end;
	  	koor(1:3,i,j) = diploc(1:3,elem(i,j));
	 	koor(4,i,j) = Y;
	    clear diff;
	  end;
  end;
end;   

disp('Computing and writing  inverse solutions (invsol)');

	offset = 0;
	for i=1:nr_names,
      for j=1:length(q),		
         if ~isempty(locations),
            rowelem = elem_G(j,:);
         else,
            if j>1, rowelem = dim*q(j-1)+1:dim*q(j);
            else, rowelem = 1:dim*q(1); end;
         end; 
         if strcmp(out_flag, 'abs') | strcmp(out_flag, 'absori')|strcmp(out_flag, 'oriabs'),
            invsolabs = inv_recon(G(rowelem,:), data(1:nr_electrodes,1:n(i),i), dim);
         end;
         if strcmp(out_flag, 'ori')|strcmp(out_flag, 'absori')|strcmp(out_flag, 'oriabs'),
            tmp = inv_recon(G(rowelem,:), data(1:nr_electrodes,1:n(i),i), 1);
            for k=1:dim,
               invsolori(1:length(rowelem)/dim,k:dim:dim*n(i)) = tmp(k:dim:length(rowelem), 1:n(i));
            end;
         end;   
         if (nargin==13)|isempty(chars),
            name = strtok(names(i,:), '.');
         else,
            name = names(i,1:chars);
         end;
         if strcmp(out_flag, 'abs')|strcmp(out_flag, 'absori')|strcmp(out_flag, 'oriabs'),
            filename = sprintf('%s%s_mn%d%s.avr', pathout, name, j, extension);
            disp(filename);
	         fid = fopen(filename, 'w');
   	      if length(latencies(1,:))>1, interval=latencies(1,2)-latencies(1,1); else, interval = 1; end;
      	   write_avr(invsolabs, filename, latencies(i,1), interval);
         end;
      	if strcmp(out_flag, 'ori')|strcmp(out_flag, 'absori')|strcmp(out_flag, 'oriabs'),
            filename = sprintf('%s/%s_mnori%d%s.avr', pathout, name, j, extension);
            disp(filename);
	         fid = fopen(filename, 'w');
   	      if length(latencies(1,:))>1, interval=(latencies(1,2)-latencies(1,1))/dim; else, interval = 1/dim; end;
            write_avr(invsolori, filename, latencies(i,1), interval);
            clear invsolori;
      	end;
		end;			% j
   end;		% i
   
   
