text:
function model = buildMNXmodel(model_type)
%BUILDMNXMODEL  Construct a model structure from MNX database files.
%
%   MODEL = BUILDMNXMODEL(MODEL_TYPE)
%
%--------------------------------- INPUTS ---------------------------------
%
%   model_type      (OPTIONAL) Specify the type of model that will be
%                   built, which determines what information is included:
%
%                   'model' (DEFAULT) All metabolite, reaction, and
%                           compartment information is included in model
%                           construction. However, note that some mets and
%                           rxns will be excluded due to issues with lack
%                           of uniqueness and metabolites that don't exist
%                           in any reactions. Therefore, the model will not
%                           contain every rxn and met in the MXN database.
%
%                   'met'   Include only metabolite-related fields. This
%                           will include ALL metabolites in the MNX
%                           database, and will therefore contain more met
%                           information than if the 'model' option is
%                           specified.
%
%                   'rxn'   Included only reaction-related fields. No
%                           stoichiometry matrix will be constructed, but
%                           the model will include ALL reactions from the
%                           MNX database.
%
%                  'both'   Include fields related to either reactions or
%                           metabolites, but do not associate them with
%                           each other (i.e., no stoichimetry matrix). This
%                           will contain all the fields that would be
%                           obtained by running both the 'met' and 'rxn'
%                           options. More information will be present in
%                           this format than if the 'model' option is
%                           specified, because reactions and/or metabolites
%                           will not be removed to construct a functioning
%                           stoichiometry matrix.
%

% Specify directory containing MNX database files.
% * Note that the original files downloaded from the database were changed
%   from .tsv to .txt, and all header lines (those beginning with #) were
%   removed, except the last one, which was kept, but the leading '#' was
%   removed.
MNXdir = '';
if isempty(MNXdir)
    error('The MNX directory needs to be specified within the buildMNXmodel function.')
end

% handle inputs
if nargin < 1
    model_type = 'model';
end

%% Initialize model
% use similar field order as those from SBML
model=[];
model.id='MNXdatabase';
model.description='';
model.rxns={};
model.mets={};
model.S=[];
model.lb=[];
model.ub=[];
model.rev=[];
model.c=[];
model.b=[];
model.comps={};
model.compNames={};
model.rxnNames={};
model.grRules={};
model.rxnGeneMat=[];
model.eccodes={};
model.genes={};
model.metNames={};
model.metComps=[];
model.inchis={};
model.metFormulas={};
% model.metMiriams={};
model.metCharges=[];

% remove unnecessary fields for met-only or rxn-only structures
if strcmpi(model_type,'model')
    fprintf('*** Building MODEL structure from MNX Database (reactions AND metabolites, with stoichiometry matrix) ***\n\n')
    model.description = 'MetaNetX Database Model';
elseif strcmpi(model_type,'met')
    fprintf('*** Building information structure for MNX METABOLITES ***\n\n')
    model.description = 'MetaNetX Metabolite Data';
    model = rmfield(model,{'rxns','S','lb','ub','rev','c','b','comps','compNames','rxnNames','grRules','rxnGeneMat','eccodes','genes','metComps'});
elseif strcmpi(model_type,'rxn')
    fprintf('*** Building information structure for MNX REACTIONS ***\n\n')
    model.description = 'MetaNetX Reaction Data';
    model = rmfield(model,{'mets','S','lb','ub','rev','c','b','comps','compNames','rxnNames','grRules','rxnGeneMat','genes','metNames','metComps','inchis','metFormulas','metCharges'});
elseif strcmpi(model_type,'both')
    fprintf('*** Building information structure for MNX REACTIONS AND METABOLITES ***\n\n')
    model.description = 'MetaNetX Reaction and Metabolite Data';
    model = rmfield(model,{'S','lb','ub','rev','c','b','comps','compNames','grRules','rxnGeneMat','genes','metComps'});
else
    error('Invalid MODEL_TYPE argument.');
end

%% Load compartment properties

if strcmpi(model_type,'model')
    fprintf('Loading model compartment data... ');
    mnx = readtable([MNXdir,'comp_prop.txt'],'Delimiter','\t');

    % these compartment abbreviations will probably cause problems for some
    % functions, as they are not a single letter, but a full string
    model.comps = mnx.MNX_ID;
    model.compNames = mnx.Description;
    fprintf('Done.\n');
end

%% Load reaction properties

if ismember(model_type,{'model','rxn','both'})

    fprintf('Loading reaction data... ');
    mnx = readtable([MNXdir,'reac_prop.txt'],'Delimiter','\t');
    fprintf('Done.\n');

    if strcmp(model_type,'model')

        % remove reactions with ambiguous stoich coeff; e.g., (n), (N), (n-1), (2n), etc.
        fprintf('Processing reaction data... ');
        del_rxns = contains(mnx.Equation,'(');
        mnx(del_rxns,:) = [];

        % also remove some reactions that are duplicates, or imbalanced
        % (these were identified through a separate analysis)
        del_rxns = ismember(mnx.MNX_ID,{'MNXR109831','MNXR126812','MNXR117558','MNXR107762','MNXR135187','MNXR127223'});
        mnx(del_rxns,:) = [];

        % split each rxn equations into reactants and products
        eqns = split(mnx.Equation,' = ');
        reactants = eqns(:,1);
        products = eqns(:,2);

        % Remove reactions that have the same met as both a reactant and a product,
        % and with the same stoich coeff. This happens in cases where an enzyme has
        % been included in the reaction equation, or there is a mistake in the
        % reaction, or MNX is assuming that the two forms of a compound are
        % identical, when maybe it should not. Many are the result of the database
        % selecting the form of the compound that is most abundant at pH 7.3;
        % for example, the rxn: NH3 + H -> NH4 is  NH4 -> NH4 in the database.
        smet_reac = cellfun(@(r) regexp(r,'[0-9.]+ \w+@\w+','match'),reactants,'UniformOutput',false);
        smet_prod = cellfun(@(r) regexp(r,'[0-9.]+ \w+@\w+','match'),products,'UniformOutput',false);
        del_rxns = cellfun(@(reac,prod) ~isempty(intersect(reac,prod)),smet_reac,smet_prod);
        mnx(del_rxns,:) = [];
        reactants(del_rxns) = [];
        products(del_rxns) = [];

    else

        % add a field that contains a list of mets that participate in each reaction
        rxnMets = cellfun(@(r) regexp(r,'(\w+)@\w+','tokens'),mnx.Equation,'UniformOutput',false);
        model.rxnMets = cellfun(@(r) unique([r{:}]),rxnMets,'UniformOutput',false);

    end

    % extract reaction MNX IDs and other information
    model.rxns = mnx.MNX_ID;
    model.rxnBalanced = mnx.Balance;
    model.eccodes = mnx.EC;
    model.rxnNames = repmat({''},size(model.rxns));  % no rxnNames present
    model.rxnEqnMNX = mnx.Equation;
    model.rxnEqnNames = mnx.Description;

    if strcmp(model_type,'model')

        % associate every reaction with each of its constituent metabolites
        reactant_mets = cellfun(@(r) regexp(r,'\w+@\w+','match'),reactants,'UniformOutput',false);
        reactant_rxn_names = cellfun(@(rname,r) repmat({rname},1,numel(r)),mnx.MNX_ID,reactant_mets,'UniformOutput',false);
        reactant_rxn_met_pairs = [[reactant_mets{:}]',[reactant_rxn_names{:}]'];

        product_mets = cellfun(@(r) regexp(r,'\w+@\w+','match'),products,'UniformOutput',false);
        product_rxn_names = cellfun(@(rname,r) repmat({rname},1,numel(r)),mnx.MNX_ID,product_mets,'UniformOutput',false);
        product_rxn_met_pairs = [[product_mets{:}]',[product_rxn_names{:}]'];

        rxn_met_pairs = [reactant_rxn_met_pairs; product_rxn_met_pairs];

        % get a unique list of all metabolite-compartment combinations
        model.mets = unique([[reactant_mets{:}]';[product_mets{:}]']);
        model.metMNXID = regexprep(model.mets,'@.+$','');  % get list of metabolites without compartment names
        uniq_mets_nocomp = unique(model.metMNXID);  % unique list without compartments

        % extract metabolite compartments
        metComps = regexprep(model.mets,'^.+@','');
        [~,model.metComps] = ismember(metComps,model.comps);

        % obtain stoich coeffs
        reactant_stoichs = cellfun(@(r) -str2double(regexp(r,'^[0-9.]+ | [0-9.]+ ','match')), reactants,'UniformOutput',false);
        product_stoichs = cellfun(@(r) str2double(regexp(r,'^[0-9.]+ | [0-9.]+ ','match')), products,'UniformOutput',false);
        stoich_coeffs = [[reactant_stoichs{:}]';[product_stoichs{:}]'];

        % construct stoichiometry matrix
        [~,met_ind] = ismember(rxn_met_pairs(:,1),model.mets);
        [~,rxn_ind] = ismember(rxn_met_pairs(:,2),model.rxns);
        model.S = sparse(met_ind,rxn_ind,stoich_coeffs);

    end

    % get other rxn IDs
    rxnSources = mnx.Source;
    rxnSources(~contains(rxnSources,':')) = {''};  % only keep entries with a colon, which indicates they have a source beyond MNX
    rxnSourceNames = regexprep(rxnSources,':.+$','');  % retrieve source name
    rxnSourceIDs = regexprep(rxnSources,'^.+:','');  % retrieve source ID

    % add model fields corresponding to available source names
    model.rxnMNXID = model.rxns;
    model.rxnBiGGID = repmat({''},size(model.rxns));
    model.rxnBiGGID(ismember(rxnSourceNames,'bigg')) = rxnSourceIDs(ismember(rxnSourceNames,'bigg'));
    model.rxnKEGGID = repmat({''},size(model.rxns));
    model.rxnKEGGID(ismember(rxnSourceNames,'kegg')) = rxnSourceIDs(ismember(rxnSourceNames,'kegg'));
    model.rxnMetaCycID = repmat({''},size(model.rxns));
    model.rxnMetaCycID(ismember(rxnSourceNames,'metacyc')) = rxnSourceIDs(ismember(rxnSourceNames,'metacyc'));
    model.rxnREACTOMEID = repmat({''},size(model.rxns));
    model.rxnREACTOMEID(ismember(rxnSourceNames,'reactome')) = rxnSourceIDs(ismember(rxnSourceNames,'reactome'));
    model.rxnRheaID = repmat({''},size(model.rxns));
    model.rxnRheaID(ismember(rxnSourceNames,'rhea')) = rxnSourceIDs(ismember(rxnSourceNames,'rhea'));
    model.rxnSABIORKID = repmat({''},size(model.rxns));
    model.rxnSABIORKID(ismember(rxnSourceNames,'sabiork')) = rxnSourceIDs(ismember(rxnSourceNames,'sabiork'));
    model.rxnSEEDID = repmat({''},size(model.rxns));
    model.rxnSEEDID(ismember(rxnSourceNames,'seed')) = rxnSourceIDs(ismember(rxnSourceNames,'seed'));

    fprintf('Done.\n');

end

%% Load metabolite properties

if ismember(model_type,{'model','met','both'})

    fprintf('Loading metabolite data... ');
    mnx = readtable([MNXdir,'chem_prop.txt'],'Delimiter','\t');  % ~45 sec load time
    mnx_xref = readtable([MNXdir,'chem_xref.txt'],'Delimiter','\t');
    fprintf('Done.\n');

    fprintf('Processing metabolite data... ');
    if strcmp(model_type,'model')
        % retrieve row indices corresponding to model met list
        [~,met_ind] = ismember(model.metMNXID,mnx.MNX_ID);
        mnx = mnx(met_ind,:);
    else
        model.mets = mnx.MNX_ID;
        model.metMNXID = mnx.MNX_ID;
    end

    % extract information and add to model
    model.metNames = mnx.Description;
    model.metFormulas = mnx.Formula;
    model.metCharges = mnx.Charge;
    model.inchis = mnx.InChI;
    model.metSMILES = mnx.SMILES;

    % get other IDs
    metSources = mnx.Source;
    metSources(~contains(metSources,':')) = {''};  % only keep entries with a colon, which indicates they have a source beyond MNX
    metSourceNames = regexprep(metSources,':.+$','');  % retrieve source name
    metSourceIDs = regexprep(metSources,'^.+:','');  % retrieve source ID

    % remove rows of mnx_xref that don't contain external IDs
    mnx_xref(~contains(mnx_xref.XREF,':') | contains(mnx_xref.XREF,'deprecated:'),:) = [];

    % extract source and ID information from MNX xref data
    metSourceNamesX = regexprep(mnx_xref.XREF,':.+$','');  % retrieve source name
    metSourceIDsX = regexprep(mnx_xref.XREF,'^.+:','');  % retrieve source ID
    mnxID2extID = [mnx_xref.MNX_ID,metSourceNamesX,metSourceIDsX];

    % Split met name field by delimiters '|' and '; '. This takes a very
    % long time, and should only be done if necessary.
    splitnames = true;  % CHANGE TO FALSE FOR FASTER PROCESSING
    if ( splitnames )
        % change semicolon+space delimiters to vertical bar "|"
        regexprep(mnx_xref.Description,'; ','|');

        % break processing into 10 chunks
        chunk_ind = 1:round(length(mnx_xref.Description)/10):length(mnx_xref.Description);
        chunk_ind(end) = length(mnx_xref.Description)+1;

        % split met names
        mnxDescr = {};
        fprintf('\nSplitting MNX metabolite names: ');
        for i = 1:length(chunk_ind)-1
            fprintf('%u%% ',round(chunk_ind(i)/length(mnx_xref.Description)*100));
            mnxDescr = [mnxDescr; cellfun(@(x) strsplit(x,'|'),mnx_xref.Description(chunk_ind(i):chunk_ind(i+1)-1),'UniformOutput',false)];
        end
        fprintf('Done.\n');

        % flatten nested cells to obtain an Nx2 matrix of MNXID-name pairs
        fprintf('Organizing MNXID-name pairs... ');
        mnxIDs = cellfun(@(id,descr) repmat({id},1,numel(descr)),mnx_xref.MNX_ID,mnxDescr,'UniformOutput',false);
        model.mnxID2name = [[mnxIDs{:}]',[mnxDescr{:}]'];

        % remove repeated ID-name pairs
        [~,numericPair] = ismember(lower(model.mnxID2name),unique(lower(model.mnxID2name)));  % make numeric for faster processing
        [~,uniq_ind] = unique(numericPair,'rows');
        model.mnxID2name = model.mnxID2name(uniq_ind,:);
        fprintf('Done.\n');

    else
        % add unsplit name field as additional column
        mnxID2extID = [mnxID2extID,mnx_xref.Description];
    end

    % add model fields corresponding to available source names
    sNames = {'bigg','chebi','envipath','hmdb','kegg','lipidmaps','metacyc','reactome','sabiork','seed','slm'};
    fNames = {'metBiGGID','metChEBIID','metEnviPathID','metHMDBID','metKEGGID','metLIPIDMAPSID','metMetaCycID','metREACTOMEID','metSABIORKID','metSEEDID','metSLMID'};

    % This portion of code would be nice, but it takes forever to run.
%     for i = 1:length(sNames)
%         fprintf('Retrieving %s met IDs... ',upper(sNames{i}));
%         id2id = mnxID2extID(ismember(metSourceNamesX,sNames{i}),:);  % subset ID data for speed
%         ind = ismember(model.mets,id2id(:,1));  % subset met list for speed
%         model.(fNames{i}) = repmat({''},size(model.mets));
%         model.(fNames{i})(ind) = cellfun(@(id) id2id(ismember(id2id(:,1),id),3),model.mets(ind),'UniformOutput',false);
%         fprintf('Done.\n');
%     end

    % instead, just save the unprocessed information
    [~,ind] = ismember(mnxID2extID(:,2),sNames);
    mnxID2extID(:,2) = fNames(ind);  % rename to match field names
    model.mnxID2extID = mnxID2extID;

    % fill in other fields
    model.metBiGGID = repmat({''},size(model.mets));
    model.metBiGGID(ismember(metSourceNames,'bigg')) = metSourceIDs(ismember(metSourceNames,'bigg'));
    model.metChEBIID = repmat({''},size(model.mets));
    model.metChEBIID(ismember(metSourceNames,'chebi')) = metSourceIDs(ismember(metSourceNames,'chebi'));
    model.metEnviPathID = repmat({''},size(model.mets));
    model.metEnviPathID(ismember(metSourceNames,'envipath')) = metSourceIDs(ismember(metSourceNames,'envipath'));
    model.metHMDBID = repmat({''},size(model.mets));
    model.metHMDBID(ismember(metSourceNames,'hmdb')) = metSourceIDs(ismember(metSourceNames,'hmdb'));
    model.metKEGGID = repmat({''},size(model.mets));
    model.metKEGGID(ismember(metSourceNames,'kegg')) = metSourceIDs(ismember(metSourceNames,'kegg'));
    model.metLIPIDMAPSID = repmat({''},size(model.mets));
    model.metLIPIDMAPSID(ismember(metSourceNames,'lipidmaps')) = metSourceIDs(ismember(metSourceNames,'lipidmaps'));
    model.metMetaCycID = repmat({''},size(model.mets));
    model.metMetaCycID(ismember(metSourceNames,'metacyc')) = metSourceIDs(ismember(metSourceNames,'metacyc'));
    model.metREACTOMEID = repmat({''},size(model.mets));
    model.metREACTOMEID(ismember(metSourceNames,'reactome')) = metSourceIDs(ismember(metSourceNames,'reactome'));
    model.metSABIORKID = repmat({''},size(model.mets));
    model.metSABIORKID(ismember(metSourceNames,'sabiork')) = metSourceIDs(ismember(metSourceNames,'sabiork'));
    model.metSEEDID = repmat({''},size(model.mets));
    model.metSEEDID(ismember(metSourceNames,'seed')) = metSourceIDs(ismember(metSourceNames,'seed'));
    model.metSLMID = repmat({''},size(model.mets));  % SwissLipids database
    model.metSLMID(ismember(metSourceNames,'slm')) = metSourceIDs(ismember(metSourceNames,'slm'));

    fprintf('Done.\n');

    % %% Merge reactions with identical stoichiometry
    %
    % % group reactions with identical stoich columns
    % [~,~,rxn_groups] = unique(model.S','rows');
    %
    % % determine which groups have more than one member
    % rep_inds = ismember(rxn_groups,find(histcounts(rxn_groups,'BinMethod','integers') > 1));
    % rep_groups = unique(rxn_groups(rep_inds));

end

%% Final model adjustments

if strcmp(model_type,'model')
    model.c = zeros(size(model.rxns));
    model.b = zeros(size(model.mets));
    model.lb = -1000*ones(size(model.rxns));
    model.ub = 1000*ones(size(model.rxns));
    model.rev = true(size(model.rxns));

    model.genes = {'dummy'};
    model.grRules = repmat({''},size(model.rxns));
    model.rxnGeneMat = zeros(size(model.rxns));
end

if ismember(model_type,{'model','met','both'})
    % remove NAs from metFormulas
    model.metFormulas(ismember(model.metFormulas,'NA')) = {''};
    if strcmp(model_type,'model')
        % change format of mets from "met@comp" to "met[comp]"
        model.mets = strcat(regexprep(model.mets,'@','['),']');
    end
end

author: SysBioChalmers

repo: Human-GEM


text:
function []=ps_baseline_plot(linecolor)
%PS_BASELINE_PLOT plot single master baselines
%
%   Andy Hooper, March 2011
%
%   ======================================================================
%   ======================================================================
if nargin<1
    linecolor=[1 0 0];
end

load psver
psname=['ps',num2str(psver)];
ps=load(psname);

small_baseline_flag=getparm('small_baseline_flag');
if ~strcmpi(small_baseline_flag,'n')
   error('You are not in a PS directory')
end

m=ps.master_ix;

%clf

for i=1:length(ps.day)
    l=line([ps.day(i),ps.day(m)],[ps.bperp(i),ps.bperp(m)]);
    set(l,'color',linecolor,'linewidth',2)
end
hold on
set(gca,'fontsize',14)
p=plot(ps.day,ps.bperp,'k+');
set(p,'markersize',8,'linewidth',1)
hold off
datetick('x',10)
xlabel('Acquisition Date')
ylabel('Perpendicular Baseline (m)')

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function [xyQ,stateindex] = ExtractQ_XY_2D_POTTS(state,x,y,Q)

display('Extracting Q_[X, Y] values')
xyQ=cell(1,Q);
stateindex=xyQ;
for q = 1:Q
    stateindex{1,q}=find(state(:,:)==q);
    xyQ{1,q}=[];
    for count=1:prod(size(stateindex{1,q}))
        ElementRnCn(count,:)=FindRC_2D_QPOTTS(stateindex{1,q}(count),state);
        xyQ{1,q}=[xyQ{1,q};...
            x(ElementRnCn(count,1),ElementRnCn(count,2)) y(ElementRnCn(count,1),ElementRnCn(count,2))];
    end
end

author: Sable

repo: mcbench-benchmarks

path: github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/34784-monte-carlo-simulation-of-two-dimensional-grain-growth-code-version-no-1-basic/Monte Carlo Simulation Q-state Potts model 2D square-lattice - microstructure/ExtractQ_XY_2D_POTTS.m


text:
function logical_mask = build_logical_mask2(nt, maxNt, Ns)

%   Copyright (C) 2009 Cesare Magri
%   Version: 1.0.4

% -------
% LICENSE
% -------
% This software is distributed free under the condition that:
%
% 1. it shall not be incorporated in software that is subsequently sold;
%
% 2. the authorship of the software shall be acknowledged and the following
%    article shall be properly cited in any publication that uses results
%    generated by the software:
%
%      Magri C, Whittingstall K, Singh V, Logothetis NK, Panzeri S: A
%      toolbox for the fast information analysis of multiple-site LFP, EEG
%      and spike train recordings. BMC Neuroscience 2009 10(1):81;
%
% 3.  this notice shall remain in place in each source file.

logical_mask = (1:maxNt)'*ones(1,Ns)<=ones(maxNt,1)*nt(:)';

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [pvec, pstruct] = tapas_ph_binary_transp(r, ptrans)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2015 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pvec    = NaN(1,length(ptrans));
pstruct = struct;

pvec(1)       = tapas_sgm(ptrans(1),1); % v_0
pstruct.v_0   = pvec(1);
pvec(2)       = tapas_sgm(ptrans(2),1); % al_0
pstruct.al_0  = pvec(2);
pvec(3)       = exp(ptrans(3));         % S
pstruct.S     = pvec(3);

return;

author: translationalneuromodeling

repo: tapas


text:
function [G, ind] = subGrid(G,q,epsilon,varargin)
% sub-SO3Grid as epsilon neigborhood of a node
%
% Syntax
%   G = subGrid(G,midpoint,radius)
%
% Input
%  G        - @SO3Grid
%  midpoint - @quaternion
%  radius   - double
%
% Output
%  G - SO3Grid
%
% See also
%  SO3Grid/find S2Grid/subGrid

if nargin >= 3
  ind = find(G,q,epsilon,varargin{:});
  ind = any(ind,2);
elseif islogical(q)
  ind = q;
else
  ind = false(length(G),1);
  ind(q) = true;
end

G.a = G.a(ind);
G.b = G.b(ind);
G.c = G.c(ind);
G.d = G.d(ind);
G.i = G.i(ind);

G.gamma = subGrid(G.gamma,ind);
G.alphabeta  = subGrid(G.alphabeta,GridLength(G.gamma)>0);
G.gamma(GridLength(G.gamma)==0) = [];

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [data] = ft_checkdata(data, varargin)

% FT_CHECKDATA checks the input data of the main FieldTrip functions, e.g. whether
% the type of data strucure corresponds with the required data. If neccessary
% and possible, this function will adjust the data structure to the input
% requirements (e.g. change dimord, average over trials, convert inside from
% index into logical).
%
% If the input data does NOT correspond to the requirements, this function
% is supposed to give a elaborate warning message and if applicable point
% the user to external documentation (link to website).
%
% Use as
%   [data] = ft_checkdata(data, ...)
%
% Optional input arguments should be specified as key-value pairs and can include
%   feedback           = yes, no
%   datatype           = raw, freq, timelock, comp, spike, source, volume, dip
%   dimord             = any combination of time, freq, chan, refchan, rpt, subj, chancmb, rpttap, pos
%   senstype           = ctf151, ctf275, ctf151_planar, ctf275_planar, neuromag122, neuromag306, bti148, bti248, bti248_planar, magnetometer, electrode
%   inside             = logical, index
%   ismeg              = yes, no
%   hastrials          = yes, no
%   hasunits           = yes, no
%   hastrialdef        = yes, no
%   hasoffset          = yes, no (only applies to raw data)
%   hascumtapcnt       = yes, no (only applies to freq data)
%   hasdof             = yes, no
%   cmbrepresentation  = sparse, full (applies to covariance and cross-spectral density)
%
% For some options you can specify multiple values, e.g.
%   [data] = ft_checkdata(data, 'senstype', {'ctf151', 'ctf275'}), e.g. in megrealign
%   [data] = ft_checkdata(data, 'datatype', {'timelock', 'freq'}), e.g. in sourceanalysis

% Copyright (C) 2007-2009, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: ft_checkdata.m 2886 2011-02-16 09:42:35Z roboos $

% in case of an error this function could use dbstack for more detailled
% user feedback
%
% this function should replace/encapsulate
%   fixdimord
%   fixinside
%   fixprecision
%   fixvolume
%   data2raw
%   raw2data
%   grid2transform
%   transform2grid
%   fourier2crsspctrm
%   freq2cumtapcnt
%   sensortype
%   time2offset
%   offset2time
%
% other potential uses for this function:
%   time -> offset in freqanalysis
%   average over trials
%   csd as matrix

% get the optional input arguments
feedback      = keyval('feedback',      varargin); if isempty(feedback), feedback = 'no'; end
dtype         = keyval('datatype',      varargin); % should not conflict with the ft_datatype function
dimord        = keyval('dimord',        varargin);
stype         = keyval('senstype',      varargin); % senstype is a function name which should not be masked
ismeg         = keyval('ismeg',         varargin);
inside        = keyval('inside',        varargin); % can be logical or index
hastrials     = keyval('hastrials',     varargin);
hasunits      = keyval('hasunits',      varargin);
hastrialdef   = keyval('hastrialdef',   varargin); if isempty(hastrialdef), hastrialdef = 'no'; end
hasoffset     = keyval('hasoffset',     varargin); if isempty(hasoffset), hasoffset = 'no'; end
hasdimord     = keyval('hasdimord',     varargin); if isempty(hasdimord), hasdimord = 'no'; end
hascumtapcnt  = keyval('hascumtapcnt',  varargin);
hasdof        = keyval('hasdof',        varargin); if isempty(hasdof), hasdof = 'no'; end
haspow        = keyval('haspow',        varargin); if isempty(haspow), haspow = 'no'; end
cmbrepresentation = keyval('cmbrepresentation',  varargin);
channelcmb    = keyval('channelcmb',   varargin);
sourcedimord  = keyval('sourcedimord', varargin);
sourcerepresentation = keyval('sourcerepresentation', varargin);

% determine the type of input data
% this can be raw, freq, timelock, comp, spike, source, volume, dip
israw      = ft_datatype(data, 'raw');
isfreq     = ft_datatype(data, 'freq');
istimelock = ft_datatype(data, 'timelock');
iscomp     = ft_datatype(data, 'comp');
isspike    = ft_datatype(data, 'spike');
isvolume   = ft_datatype(data, 'volume');
issource   = ft_datatype(data, 'source');
isdip      = ft_datatype(data, 'dip');
ismvar     = ft_datatype(data, 'mvar');
isfreqmvar = ft_datatype(data, 'freqmvar');

% FIXME use the istrue function on ismeg and hasxxx options

if ~isequal(feedback, 'no')
  if israw
    nchan = length(data.label);
    ntrial = length(data.trial);
    fprintf('the input is raw data with %d channels and %d trials\n', nchan, ntrial);
  elseif isfreq
    nchan = length(data.label);
    nfreq = length(data.freq);
    if isfield(data, 'time'), ntime = num2str(length(data.time)); else ntime = 'no'; end
    fprintf('the input is freq data with %d channels, %d frequencybins and %s timebins\n', nchan, nfreq, ntime);
  elseif istimelock
    nchan = length(data.label);
    ntime = length(data.time);
    fprintf('the input is timelock data with %d channels and %d timebins\n', nchan, ntime);
  elseif iscomp
    ncomp = length(data.label);
    nchan = length(data.topolabel);
    fprintf('the input is component data with %d components and %d original channels\n', ncomp, nchan);
  elseif isspike
    nchan = length(data.label);
    fprintf('the input is spike data\n');
  elseif isvolume
    fprintf('the input is volume data with dimensions [%d %d %d]\n', data.dim(1), data.dim(2), data.dim(3));
  elseif issource
    nsource = size(data.pos, 1);
    fprintf('the input is source data with %d positions\n', nsource);
  elseif isdip
    fprintf('the input is dipole data\n');
  elseif ismvar
    fprintf('the input is mvar data\n');
  elseif isfreqmvar
    fprintf('the input is freqmvar data\n');
  end
end % give feedback

if isfreq || istimelock || iscomp || issource || isvolume
  % ensure consistency between the dimord string and the axes that describe the data dimensions
  %data = fixdimord(data, strcmp(sourcerepresentation, 'new'));
  data = fixdimord(data);
end

if istimelock
  % remove the unwanted fields
  if isfield(data, 'numsamples'),       data = rmfield(data, 'numsamples');       end
  if isfield(data, 'numcovsamples'),    data = rmfield(data, 'numcovsamples');    end
  if isfield(data, 'numblcovsamples'),  data = rmfield(data, 'numblcovsamples');  end
end

if issource && isvolume
  % it should be either one or the other
  % the choice here is to represent it as volume description since that is simpler to handle
  % remove the unwanted fields
  if isfield(data, 'pos'),    data = rmfield(data, 'pos');    end
  if isfield(data, 'xgrid'),  data = rmfield(data, 'xgrid');  end
  if isfield(data, 'ygrid'),  data = rmfield(data, 'ygrid');  end
  if isfield(data, 'zgrid'),  data = rmfield(data, 'zgrid');  end
  issource = false;
end

if issource || isvolume
  % these don't contain a dimord in the old representation
  % but in the frequency domain case they could
  % contain a .frequency field rather than a .freq field
  if isfield(data, 'frequency'),
    data.freq = data.frequency;
    data      = rmfield(data, 'frequency');
  end
end

if ~isempty(dtype)
  if ~isa(dtype, 'cell')
    dtype = {dtype};
  end

  okflag = 0;
  for i=1:length(dtype)
    % check that the data matches with one or more of the required ft_datatypes
    switch dtype{i}
      case 'raw'
        okflag = okflag + israw;
      case 'freq'
        okflag = okflag + isfreq;
      case 'timelock'
        okflag = okflag + istimelock;
      case 'comp'
        okflag = okflag + iscomp;
      case 'spike'
        okflag = okflag + isspike;
      case 'volume'
        okflag = okflag + isvolume;
      case 'source'
        okflag = okflag + issource;
      case 'dip'
        okflag = okflag + isdip;
      case 'mvar'
        okflag = okflag + ismvar;
      case 'freqmvar'
        okflag = okflag + isfreqmvar;
    end % switch dtype
  end % for dtype

  if ~okflag
    % try to convert the data
    for iCell = 1:length(dtype)
      if isequal(dtype(iCell), {'source'}) && isvolume
        data = volume2source(data);
        isvolume = 0;
        issource = 1;
        okflag = 1;
      elseif isequal(dtype(iCell), {'volume'}) && issource
        data = source2volume(data);
        isvolume = 1;
        issource = 0;
        okflag = 1;
      elseif isequal(dtype(iCell), {'raw'}) && issource
        data = data2raw(data);
        issource = 0;
        israw = 1;
        okflag = 1;
      elseif isequal(dtype(iCell), {'raw'}) && istimelock
        data = timelock2raw(data);
        istimelock = 0;
        israw = 1;
        okflag = 1;
      elseif isequal(dtype(iCell), {'timelock'}) && israw
        data = raw2timelock(data);
        israw = 0;
        istimelock = 1;
        okflag = 1;
      elseif isequal(dtype(iCell), {'raw'}) && isfreq
        data = freq2raw(data);
        isfreq = 0;
        israw = 1;
        okflag = 1;
      elseif isequal(dtype(iCell), {'raw'}) && iscomp
        data = comp2raw(data);
        iscomp = 0;
        israw = 1;
        okflag = 1;
      end
    end % for iCell
  end % if okflag

  if ~okflag
    % construct an error message
    if length(dtype)>1
      str = sprintf('%s, ', dtype{1:(end-2)});
      str = sprintf('%s%s or %s', str, dtype{end-1}, dtype{end});
    else
      str = dtype{1};
    end
    str = sprintf('This function requires %s data as input.', str);
    error(str);
  end % if okflag
end

if ~isempty(dimord)
  if ~isa(dimord, 'cell')
    dimord = {dimord};
  end

  if isfield(data, 'dimord')
    okflag = any(strcmp(data.dimord, dimord));
  else
    okflag = 0;
  end

  if ~okflag
    % construct an error message
    if length(dimord)>1
      str = sprintf('%s, ', dimord{1:(end-2)});
      str = sprintf('%s%s or %s', str, dimord{end-1}, dimord{end});
    else
      str = dimord{1};
    end
    str = sprintf('This function requires data with a dimord of %s.', str);
    error(str);
  end % if okflag
end

if ~isempty(stype)
  if ~isa(stype, 'cell')
    stype = {stype};
  end

  if isfield(data, 'grad') || isfield(data, 'elec')
    if any(strcmp(ft_senstype(data), stype));
      okflag = 1;
    else
      okflag = 0;
    end
  else
    okflag = 0;
  end

  if ~okflag
    % construct an error message
    if length(stype)>1
      str = sprintf('%s, ', stype{1:(end-2)});
      str = sprintf('%s%s or %s', str, stype{end-1}, stype{end});
    else
      str = stype{1};
    end
    str = sprintf('This function requires %s data as input, but you are giving %s data.', str, ft_senstype(data));
    error(str);
  end % if okflag
end

if ~isempty(ismeg)
  if isequal(ismeg, 'yes')
    okflag = isfield(data, 'grad');
  elseif isequal(ismeg, 'no')
    okflag = ~isfield(data, 'grad');
  end

  if ~okflag && isequal(ismeg, 'yes')
    error('This function requires MEG data with a ''grad'' field');
  elseif ~okflag && isequal(ismeg, 'no')
    error('This function should not be given MEG data with a ''grad'' field');
  end % if okflag
end

if ~isempty(inside)
  % TODO absorb the fixinside function into this code
  data   = fixinside(data, inside);
  okflag = isfield(data, 'inside');

  if ~okflag
    % construct an error message
    error('This function requires data with an ''inside'' field.');
  end % if okflag
end

%if isvolume
%  % ensure consistent dimensions of the volumetric data
%  % reshape each of the volumes that is found into a 3D array
%  param = parameterselection('all', data);
%  dim   = data.dim;
%  for i=1:length(param)
%    tmp  = getsubfield(data, param{i});
%    tmp  = reshape(tmp, dim);
%    data = setsubfield(data, param{i}, tmp);
%  end
%end

if issource || isvolume,
  % these are not used any more
  if isfield(data, 'xgrid'),  data = rmfield(data, 'xgrid');  end
  if isfield(data, 'ygrid'),  data = rmfield(data, 'ygrid');  end
  if isfield(data, 'zgrid'),  data = rmfield(data, 'zgrid');  end

  if isequal(hasunits, 'yes') && ~isfield(data, 'units')
    % calling convert_units with only the input data adds the units without converting
    data = ft_convert_units(data);
  end

  % the following section is to make a dimord-consistent representation of
  % volume and source data, taking trials, time and frequency into account
  if isequal(hasdimord, 'yes') && (~isfield(data, 'dimord') || ~strcmp(data.dimord,sourcedimord))

    % determine the size of the data
    if isfield(data, 'dimord'),
      dimtok = tokenize(data.dimord, '_');
      if ~isempty(strmatch('time', dimtok)), Ntime = length(data.time); else Ntime = 1; end
      if ~isempty(strmatch('freq', dimtok)), Nfreq = length(data.freq); else Nfreq = 1; end
    else
      Nfreq = 1;
      Ntime = 1;
    end

    %convert old style source representation into new style
    if isfield(data, 'avg') && isfield(data.avg, 'mom') && (isfield(data, 'freq') || isfield(data, 'frequency')) && strcmp(sourcedimord, 'rpt_pos'),
      %frequency domain source representation convert to single trial power
      Npos   = size(data.pos,1);
      Nrpt   = length(data.cumtapcnt);
      tmpmom = zeros(Npos, size(data.avg.mom{data.inside(1)},2));
      tmpmom(data.inside,:) = cat(1,data.avg.mom{data.inside});
      tmppow = zeros(Npos, Nrpt);
      tapcnt = [0;cumsum(data.cumtapcnt)];
      for k = 1:Nrpt
        Ntap = tapcnt(k+1)-tapcnt(k);
        tmppow(data.inside,k) = sum(abs(tmpmom(data.inside,(tapcnt(k)+1):tapcnt(k+1))).^2,2)./Ntap;
      end
      data.pow = tmppow';
      data     = rmfield(data, 'avg');
      if strcmp(inside, 'logical'),
        data     = fixinside(data, 'logical');
        data.inside = repmat(data.inside(:)',[Nrpt 1]);
      end
    elseif isfield(data, 'avg') && isfield(data.avg, 'mom') && (isfield(data, 'freq') || isfield(data, 'frequency')) && strcmp(sourcedimord, 'rpttap_pos'),
      %frequency domain source representation convert to single taper fourier coefficients
      Npos   = size(data.pos,1);
      Nrpt   = sum(data.cumtapcnt);
      data.fourierspctrm = complex(zeros(Nrpt, Npos), zeros(Nrpt, Npos));
      data.fourierspctrm(:, data.inside) = transpose(cat(1, data.avg.mom{data.inside}));
      data   = rmfield(data, 'avg');
    elseif isfield(data, 'avg') && isfield(data.avg, 'mom') && isfield(data, 'time') && strcmp(sourcedimord, 'pos_time'),
      Npos   = size(data.pos,1);
      Nrpt   = 1;
      tmpmom = zeros(Npos, size(data.avg.mom{data.inside(1)},2));
      tmpmom(data.inside,:) = cat(1,data.avg.mom{data.inside});
      data.mom = tmpmom;
      if isfield(data.avg, 'noise'),
        tmpnoise = data.avg.noise(:);
        data.noise = tmpnoise(:,ones(1,size(tmpmom,2)));
      end
      data = rmfield(data, 'avg');
      Ntime = length(data.time);
    elseif isfield(data, 'trial') && isfield(data.trial(1), 'mom') && isfield(data, 'time') && strcmp(sourcedimord, 'rpt_pos_time'),
      Npos   = size(data.pos,1);
      Nrpt   = length(data.trial);
      Ntime  = length(data.time);
      tmpmom = zeros(Nrpt, Npos, Ntime);
      for k = 1:Nrpt
        tmpmom(k,data.inside,:) = cat(1,data.trial(k).mom{data.inside});
      end
      data     = rmfield(data, 'trial');
      data.mom = tmpmom;
    elseif isfield(data, 'trial') && isstruct(data.trial)
      Nrpt = length(data.trial);
    else
      Nrpt = 1;
    end

    % start with an initial specification of the dimord and dim
    if (~isfield(data, 'dim') || ~isfield(data, 'dimord'))
      if issource
        % at least it should have a Nx3 pos
        data.dim    = size(data.pos, 1);
        data.dimord = 'pos';
      elseif isvolume
        % at least it should have a 1x3 dim
        data.dim    = data.dim;
        data.dimord = 'dim1_dim2_dim3';
      end
    end

    % add the additional dimensions
    if Nfreq>1
      data.dimord = [data.dimord '_freq'];
      data.dim    = [data.dim     Nfreq];
    end
    if Ntime>1
      data.dimord = [data.dimord '_time'];
      data.dim    = [data.dim     Ntime];
    end
    if Nrpt>1 && strcmp(sourcedimord, 'rpt_pos'),
      data.dimord = ['rpt_' data.dimord];
      data.dim    = [Nrpt   data.dim ];
    elseif Nrpt>1 && strcmp(sourcedimord, 'rpttap_pos'),
      data.dimord = ['rpttap_' data.dimord];
      data.dim    = [Nrpt   data.dim ];
    end

    % the nested trial structure is not compatible with dimord
    if isfield(data, 'trial') && isstruct(data.trial)
      param = fieldnames(data.trial);
      for i=1:length(param)
        if isa(data.trial(1).(param{i}), 'cell')
          concat = cell(data.dim(1), prod(data.dim(2:end)));
        else
          concat = zeros(data.dim(1), prod(data.dim(2:end)));
        end
        for j=1:length(data.trial)
          tmp = data.trial(j).(param{i});
          concat(j,:) = tmp(:);
        end % for each trial
        data.trial = rmfield(data.trial, param{i});
        data.(param{i}) = reshape(concat, data.dim);
      end % for each param
      data = rmfield(data, 'trial');
    end
  end

  % ensure consistent dimensions of the source reconstructed data
  % reshape each of the source reconstructed parameters
  if issource && prod(data.dim)==size(data.pos,1)
    dim = [prod(data.dim) 1];
  elseif issource && any(~cellfun('isempty',strfind(fieldnames(data), 'dimord')))
    dim = [size(data.pos,1) 1]; %sparsely represented source structure new style
  elseif isfield(data, 'dim'),
    dim = [data.dim 1];
  elseif isfield(data, 'dimord'),
    %HACK
    dimtok = tokenize(data.dimord, '_');
    for i=1:length(dimtok)
      if strcmp(dimtok(i), 'pos')
        dim(1,i) = size(getsubfield(data,dimtok{i}),1);
      elseif strcmp(dimtok(i), 'rpt')
        dim(1,i) = nan;
      else
        dim(1,i) = length(getsubfield(data,dimtok{i}));
      end
    end
    i = find(isnan(dim));
    if ~isempty(i)
      n = fieldnames(data);
      for ii=1:length(n)
        numels(1,ii) = numel(getfield(data,n{ii}));
      end
      nrpt = numels./prod(dim(setdiff(1:length(dim),i)));
      nrpt = nrpt(nrpt==round(nrpt));
      dim(i) = max(nrpt);
    end
    if numel(dim)==1, dim(1,2) = 1; end;
  end

  % these fields should not be reshaped
  exclude = {'cfg' 'fwhm' 'leadfield' 'q' 'rough'};
  if ~strcmp(inside, 'logical')
    % also exclude the inside/outside from being reshaped
    exclude = cat(2, exclude, {'inside' 'outside'});
  end

  param = setdiff(parameterselection('all', data), exclude);
  for i=1:length(param)
    if any(param{i}=='.')
      % the parameter is nested in a substructure, which can have multiple elements (e.g. source.trial(1).pow, source.trial(2).pow, ...)
      % loop over the substructure array and reshape for every element
      tok  = tokenize(param{i}, '.');
      sub1 = tok{1};  % i.e. this would be 'trial'
      sub2 = tok{2};  % i.e. this would be 'pow'
      tmp1 = getfield(data, sub1);
      for j=1:numel(tmp1)
        tmp2 = getfield(tmp1(j), sub2);
        tmp2 = reshape(tmp2, dim);
        tmp1(j) = setfield(tmp1(j), sub2, tmp2);
      end
      data = setfield(data, sub1, tmp1);
    else
      tmp  = getfield(data, param{i});
      tmp  = reshape(tmp, dim);
      data = setfield(data, param{i}, tmp);
    end
  end

end

if isequal(hastrials, 'yes')
  okflag = isfield(data, 'trial');
  if ~okflag
    error('This function requires data with a ''trial'' field');
  end % if okflag
end

if isequal(hastrialdef, 'yes')
  data = fixtrialdef(data);
end

if isequal(hasoffset, 'yes')
  okflag = isfield(data, 'offset');

  if ~okflag && isfield(data, 'time') && isa(data.time, 'cell')
    if ~isfield(data, 'fsample')
      data.fsample = 1/(data.time{1}(2)-data.time{1}(1));
    end
    for i=1:length(data.time);
      data.offset(i) = time2offset(data.time{i}, data.fsample);
    end
    okflag = 1;
  elseif ~okflag && ft_datatype(data, 'mvar')
    data.offset = 0;
    okflag = 1;
  end

  if ~okflag
    error('This function requires data with an ''offset'' field');
  end % if okflag

elseif isequal(hasoffset, 'no') && isfield(data, 'offset')
  data = rmfield(data, 'offset');
end % if hasoffset

if isequal(hascumtapcnt, 'yes') && ~isfield(data, 'cumtapcnt')
  error('This function requires data with a ''cumtapcnt'' field');
elseif isequal(hascumtapcnt, 'no') && isfield(data, 'cumtapcnt')
  data = rmfield(data, 'cumtapcnt');
end % if hascumtapcnt

if isequal(hasdof, 'yes') && ~isfield(data, 'hasdof')
  error('This function requires data with a ''dof'' field');
elseif isequal(hasdof, 'no') && isfield(data, 'hasdof')
  data = rmfield(data, 'cumtapcnt');
end % if hasdof

if ~isempty(cmbrepresentation)
  if istimelock
    data = fixcov(data, cmbrepresentation);
  elseif isfreq
    data = fixcsd(data, cmbrepresentation, channelcmb);
  elseif isfreqmvar
    data = fixcsd(data, cmbrepresentation, channelcmb);
  else
    error('This function requires data with a covariance, coherence or cross-spectrum');
  end
end % cmbrepresentation

if issource && ~isempty(sourcerepresentation)
  data = fixsource(data, 'type', sourcerepresentation);
end

if issource && ~strcmp(haspow, 'no')
 data = fixsource(data, 'type', sourcerepresentation, 'haspow', haspow);
end

if isfield(data, 'grad')
  % ensure that the gradiometer balancing is specified
  if ~isfield(data.grad, 'balance') || ~isfield(data.grad.balance, 'current')
    data.grad.balance.current = 'none';
  end
end

% This is necessary to ensure backward compatibility for data processed and stored in July 2010.
% We decided to change the field trialdef into sampleinfo.
if isfield(data, 'trialdef')
  data.sampleinfo = data.trialdef;
  data = rmfield(data, 'trialdef');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% represent the covariance matrix in a particular manner
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = fixcov(data, desired)
if isfield(data, 'cov')     && ~isfield(data, 'labelcmb')
  current = 'full';
elseif isfield(data, 'cov') &&  isfield(data, 'labelcmb')
  current = 'sparse';
else
  error('Could not determine the current representation of the covariance matrix');
end
if isequal(current, desired)
  % nothing to do
elseif strcmp(current, 'full') && strcmp(desired, 'sparse')
  % FIXME should be implemented
  error('not yet implemented');
elseif strcmp(current, 'sparse') && strcmp(desired, 'full')
  % FIXME should be implemented
  error('not yet implemented');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% represent the cross-spectral density matrix in a particular manner
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [data] = fixcsd(data, desired, channelcmb)

% FIXCSD converts univariate frequency domain data (fourierspctrm) into a bivariate
% representation (crsspctrm), or changes the representation of bivariate frequency
% domain data (sparse/full/sparsewithpow, sparsewithpow only works for crsspctrm or
% fourierspctrm)

% Copyright (C) 2010, Jan-Mathijs Schoffelen, Robert Oostenveld

if isfield(data, 'crsspctrm') && isfield(data, 'powspctrm')
  current = 'sparsewithpow';
elseif isfield(data, 'powspctrm')
  current = 'sparsewithpow';
elseif isfield(data, 'fourierspctrm') && ~isfield(data, 'labelcmb')
  current = 'fourier';
elseif ~isfield(data, 'labelcmb')
  current = 'full';
elseif isfield(data, 'labelcmb')
  current = 'sparse';
else
  error('Could not determine the current representation of the %s matrix', param);
end

% first go from univariate fourier to the required bivariate representation
if strcmp(current, 'fourier') && strcmp(desired, 'fourier')
  % nothing to do
elseif strcmp(current, 'fourier') && strcmp(desired, 'sparsewithpow')
  dimtok = tokenize(data.dimord, '_');
  if ~isempty(strmatch('rpttap',   dimtok)),
    nrpt = length(data.cumtapcnt);
    flag = 0;
  else
    nrpt = 1;
    flag = 1;
  end
  if ~isempty(strmatch('freq',  dimtok)), nfrq=length(data.freq);      else nfrq = 1; end
  if ~isempty(strmatch('time',  dimtok)), ntim=length(data.time);      else ntim = 1; end

  fastflag = all(data.cumtapcnt(:)==data.cumtapcnt(1));

  %create auto-spectra
  nchan     = length(data.label);
  if fastflag
    % all trials have the same amount of tapers
    powspctrm = zeros(nrpt,nchan,nfrq,ntim);
    ntap      = data.cumtapcnt(1);
    for p = 1:ntap
      powspctrm = powspctrm + abs(data.fourierspctrm(p:ntap:end,:,:,:,:)).^2;
    end
    powspctrm = powspctrm./ntap;
  else
    % different amount of tapers
    powspctrm = zeros(nrpt,nchan,nfrq,ntim)+i.*zeros(nrpt,nchan,nfrq,ntim);
    sumtapcnt = [0;cumsum(data.cumtapcnt(:))];
    for p = 1:nrpt
      indx   = (sumtapcnt(p)+1):sumtapcnt(p+1);
      tmpdat = data.fourierspctrm(indx,:,:,:);
      powspctrm(p,:,:,:) = (sum(tmpdat.*conj(tmpdat),1))./data.cumtapcnt(p);
    end
  end

  %create cross-spectra
  if ~isempty(channelcmb),
    ncmb      = size(channelcmb,1);
    cmbindx   = zeros(ncmb,2);
    labelcmb  = cell(ncmb,2);
    for k = 1:ncmb
      ch1 = find(strcmp(data.label, channelcmb(k,1)));
      ch2 = find(strcmp(data.label, channelcmb(k,2)));
      if ~isempty(ch1) && ~isempty(ch2),
        cmbindx(k,:)  = [ch1 ch2];
        labelcmb(k,:) = data.label([ch1 ch2])';
      end
    end

    crsspctrm = zeros(nrpt,ncmb,nfrq,ntim)+i.*zeros(nrpt,ncmb,nfrq,ntim);
    if fastflag
      for p = 1:ntap
        tmpdat1   = data.fourierspctrm(p:ntap:end,cmbindx(:,1),:,:,:);
        tmpdat2   = data.fourierspctrm(p:ntap:end,cmbindx(:,2),:,:,:);
        crsspctrm = crsspctrm + tmpdat1.*conj(tmpdat2);
      end
      crsspctrm = crsspctrm./ntap;
    else
      for p = 1:nrpt
        indx    = (sumtapcnt(p)+1):sumtapcnt(p+1);
        tmpdat1 = data.fourierspctrm(indx,cmbindx(:,1),:,:);
        tmpdat2 = data.fourierspctrm(indx,cmbindx(:,2),:,:);
        crsspctrm(p,:,:,:) = (sum(tmpdat1.*conj(tmpdat2),1))./data.cumtapcnt(p);
      end
    end
    data.crsspctrm = crsspctrm;
    data.labelcmb  = labelcmb;
  end
  data.powspctrm = powspctrm;
  data           = rmfield(data, 'fourierspctrm');
  if ntim>1,
    data.dimord = 'chan_freq_time';
  else
    data.dimord = 'chan_freq';
  end

  if nrpt>1,
    data.dimord = ['rpt_',data.dimord];
  end

  if flag, siz = size(data.crsspctrm); data.crsspctrm = reshape(data.crsspctrm, siz(2:end)); end
elseif strcmp(current, 'fourier') && strcmp(desired, 'sparse')

  if isempty(channelcmb), error('no channel combinations are specified'); end
  dimtok = tokenize(data.dimord, '_');
  if ~isempty(strmatch('rpttap',   dimtok)),
    nrpt = length(data.cumtapcnt);
    flag = 0;
  else
    nrpt = 1;
    flag = 1;
  end
  if ~isempty(strmatch('freq',  dimtok)), nfrq=length(data.freq); else nfrq = 1; end
  if ~isempty(strmatch('time',  dimtok)), ntim=length(data.time); else ntim = 1; end

  ncmb      = size(channelcmb,1);
  cmbindx   = zeros(ncmb,2);
  labelcmb  = cell(ncmb,2);
  for k = 1:ncmb
    ch1 = find(strcmp(data.label, channelcmb(k,1)));
    ch2 = find(strcmp(data.label, channelcmb(k,2)));
    if ~isempty(ch1) && ~isempty(ch2),
      cmbindx(k,:)  = [ch1 ch2];
      labelcmb(k,:) = data.label([ch1 ch2])';
    end
  end

  sumtapcnt = [0;cumsum(data.cumtapcnt(:))];
  fastflag  = all(data.cumtapcnt(:)==data.cumtapcnt(1));

  if fastflag && nrpt>1
    ntap = data.cumtapcnt(1);

    % compute running sum across tapers
    siz = [size(data.fourierspctrm) 1];

    for p = 1:ntap
      indx      = p:ntap:nrpt*ntap;

      if p==1.

        tmpc = zeros(numel(indx), size(cmbindx,1), siz(3), siz(4)) + ...
           1i.*zeros(numel(indx), size(cmbindx,1), siz(3), siz(4));
      end

      for k = 1:size(cmbindx,1)
        tmpc(:,k,:,:) = data.fourierspctrm(indx,cmbindx(k,1),:,:).*  ...
                   conj(data.fourierspctrm(indx,cmbindx(k,2),:,:));
      end

      if p==1
        crsspctrm = tmpc;
      else
        crsspctrm = tmpc + crsspctrm;
      end
    end
    crsspctrm = crsspctrm./ntap;
  else
    crsspctrm = zeros(nrpt, ncmb, nfrq, ntim);
    for p = 1:nrpt
      indx    = (sumtapcnt(p)+1):sumtapcnt(p+1);
      tmpdat1 = data.fourierspctrm(indx,cmbindx(:,1),:,:);
      tmpdat2 = data.fourierspctrm(indx,cmbindx(:,2),:,:);
      crsspctrm(p,:,:,:) = (sum(tmpdat1.*conj(tmpdat2),1))./data.cumtapcnt(p);
    end
  end
  data.crsspctrm = crsspctrm;
  data.labelcmb  = labelcmb;
  data           = rmfield(data, 'fourierspctrm');
  data           = rmfield(data, 'label');
  if ntim>1,
    data.dimord = 'chan_freq_time';
  else
    data.dimord = 'chan_freq';
  end

  if nrpt>1,
    data.dimord = ['rpt_',data.dimord];
  end

  if flag, siz = size(data.crsspctrm); data.crsspctrm = reshape(data.crsspctrm, siz(2:end)); end
elseif strcmp(current, 'fourier') && strcmp(desired, 'full')

  % this is how it is currently and the desired functionality of prepare_freq_matrices
  dimtok = tokenize(data.dimord, '_');
  if ~isempty(strmatch('rpttap',   dimtok)),
    nrpt = length(data.cumtapcnt);
    flag = 0;
  else
    nrpt = 1;
    flag = 1;
  end
  if ~isempty(strmatch('rpttap',dimtok)), nrpt=length(data.cumtapcnt); else nrpt = 1; end
  if ~isempty(strmatch('freq',  dimtok)), nfrq=length(data.freq);      else nfrq = 1; end
  if ~isempty(strmatch('time',  dimtok)), ntim=length(data.time);      else ntim = 1; end
  nchan     = length(data.label);
  crsspctrm = zeros(nrpt,nchan,nchan,nfrq,ntim);
  sumtapcnt = [0;cumsum(data.cumtapcnt(:))];
  for k = 1:ntim
    for m = 1:nfrq
      for p = 1:nrpt
        %FIXME speed this up in the case that all trials have equal number of tapers
        indx   = (sumtapcnt(p)+1):sumtapcnt(p+1);
        tmpdat = transpose(data.fourierspctrm(indx,:,m,k));
        crsspctrm(p,:,:,m,k) = (tmpdat*tmpdat')./data.cumtapcnt(p);
        clear tmpdat;
      end
    end
  end
  data.crsspctrm = crsspctrm;
  data           = rmfield(data, 'fourierspctrm');

  if ntim>1,
    data.dimord = 'chan_chan_freq_time';
  else
    data.dimord = 'chan_chan_freq';
  end

  if nrpt>1,
    data.dimord = ['rpt_',data.dimord];
  end

  % remove first singleton dimension
  if flag, siz = size(data.crsspctrm); data.crsspctrm = reshape(data.crsspctrm, siz(2:end)); end

elseif strcmp(current, 'fourier') && strcmp(desired, 'fullfast'),

  dimtok = tokenize(data.dimord, '_');
  nrpt = size(data.fourierspctrm, 1);
  nchn = numel(data.label);
  nfrq = numel(data.freq);
  if ~isempty(strmatch('time',  dimtok)), ntim=numel(data.time); else ntim = 1; end

  data.fourierspctrm = reshape(data.fourierspctrm, [nrpt nchn nfrq*ntim]);
  data.fourierspctrm(~isfinite(data.fourierspctrm)) = 0;
  crsspctrm = complex(zeros(nchn,nchn,nfrq*ntim));
  for k = 1:nfrq*ntim
    tmp = transpose(data.fourierspctrm(:,:,k));
    n   = sum(tmp~=0,2);
    crsspctrm(:,:,k) = tmp*tmp'./n(1);
  end
  data           = rmfield(data, 'fourierspctrm');
  data.crsspctrm = reshape(crsspctrm, [nchn nchn nfrq ntim]);
  if isfield(data, 'time'),
    data.dimord = 'chan_chan_freq_time';
  else
    data.dimord = 'chan_chan_freq';
  end

end % convert to the requested bivariate representation

% from one bivariate representation to another
if isequal(current, desired)
  % nothing to do

elseif (strcmp(current, 'full')       && strcmp(desired, 'fourier')) || ...
    (strcmp(current, 'sparse')        && strcmp(desired, 'fourier')) || ...
    (strcmp(current, 'sparsewithpow') && strcmp(desired, 'fourier'))
  % this is not possible
  error('converting the cross-spectrum into a Fourier representation is not possible');

elseif strcmp(current, 'full') && strcmp(desired, 'sparsewithpow')
  error('not yet implemented');
elseif strcmp(current, 'sparse') && strcmp(desired, 'sparsewithpow')
  % convert back to crsspctrm/powspctrm representation: useful for plotting functions etc
  indx     = labelcmb2indx(data.labelcmb);
  autoindx = indx(indx(:,1)==indx(:,2), 1);
  cmbindx  = setdiff([1:size(indx,1)]', autoindx);

  if strcmp(data.dimord(1:3), 'rpt')
    data.powspctrm = data.crsspctrm(:, autoindx, :, :);
    data.crsspctrm = data.crsspctrm(:, cmbindx,  :, :);
  else
    data.powspctrm = data.crsspctrm(autoindx, :, :);
    data.crsspctrm = data.crsspctrm(cmbindx,  :, :);
  end
  data.label    = data.labelcmb(autoindx,1);
  data.labelcmb = data.labelcmb(cmbindx, :);

  if isempty(cmbindx)
    data = rmfield(data, 'crsspctrm');
    data = rmfield(data, 'labelcmb');
  end

elseif strcmp(current, 'full') && strcmp(desired, 'sparse')
  dimtok = tokenize(data.dimord, '_');
  if ~isempty(strmatch('rpt',   dimtok)), nrpt=numel(data.cumtapcnt); else nrpt = 1; end
  if ~isempty(strmatch('freq',  dimtok)), nfrq=numel(data.freq);      else nfrq = 1; end
  if ~isempty(strmatch('time',  dimtok)), ntim=numel(data.time);      else ntim = 1; end
  nchan    = length(data.label);
  ncmb     = nchan*nchan;
  labelcmb = cell(ncmb, 2);
  cmbindx  = zeros(nchan, nchan);
  k = 1;
  for j=1:nchan
    for m=1:nchan
      labelcmb{k, 1} = data.label{m};
      labelcmb{k, 2} = data.label{j};
      cmbindx(m,j)   = k;
      k = k+1;
    end
  end

  % reshape all possible fields
  fn = fieldnames(data);
  for ii=1:numel(fn)
    if numel(data.(fn{ii})) == nrpt*ncmb*nfrq*ntim;
      if nrpt>1,
        data.(fn{ii}) = reshape(data.(fn{ii}), nrpt, ncmb, nfrq, ntim);
      else
        data.(fn{ii}) = reshape(data.(fn{ii}), ncmb, nfrq, ntim);
      end
    end
  end
  % remove obsolete fields
  data           = rmfield(data, 'label');
  try, data      = rmfield(data, 'dof'); end
  % replace updated fields
  data.labelcmb  = labelcmb;
  if ntim>1,
    data.dimord = 'chancmb_freq_time';
  else
    data.dimord = 'chancmb_freq';
  end

  if nrpt>1,
    data.dimord = ['rpt_',data.dimord];
  end

elseif strcmp(current, 'sparsewithpow') && strcmp(desired, 'sparse')

  % this representation for sparse data contains autospectra
  % as e.g. {'A' 'A'} in labelcmb
  if isfield(data, 'crsspctrm'),
    dimtok         = tokenize(data.dimord, '_');
    catdim         = match_str(dimtok, {'chan' 'chancmb'});
    data.crsspctrm = cat(catdim, data.powspctrm, data.crsspctrm);
    data.labelcmb  = [data.label(:) data.label(:); data.labelcmb];
    data           = rmfield(data, 'powspctrm');
  else
    data.crsspctrm = data.powspctrm;
    data.labelcmb  = [data.label(:) data.label(:)];
    data           = rmfield(data, 'powspctrm');
  end
  data = rmfield(data, 'label');

elseif strcmp(current, 'sparse') && strcmp(desired, 'full')
  dimtok = tokenize(data.dimord, '_');
  if ~isempty(strmatch('rpt',   dimtok)), nrpt=numel(data.cumtapcnt); else nrpt = 1; end
  if ~isempty(strmatch('freq',  dimtok)), nfrq=numel(data.freq);      else nfrq = 1; end
  if ~isempty(strmatch('time',  dimtok)), ntim=numel(data.time);      else ntim = 1; end

  if ~isfield(data, 'label')
    data.label = unique(data.labelcmb(:));
  end

  nchan     = length(data.label);
  ncmb      = size(data.labelcmb,1);
  cmbindx   = zeros(nchan,nchan);

  for k = 1:size(data.labelcmb,1)
    ch1 = find(strcmp(data.label, data.labelcmb(k,1)));
    ch2 = find(strcmp(data.label, data.labelcmb(k,2)));
    if ~isempty(ch1) && ~isempty(ch2),
      cmbindx(ch1,ch2) = k;
    end
  end

  complete = all(cmbindx(:)~=0);

  fn = fieldnames(data);
  for ii=1:numel(fn)
    if numel(data.(fn{ii})) == nrpt*ncmb*nfrq*ntim;
      if nrpt==1,
        data.(fn{ii}) = reshape(data.(fn{ii}), [nrpt ncmb nfrq ntim]);
      end

      tmpall = nan(nrpt,nchan,nchan,nfrq,ntim);

      for j = 1:nrpt
        for k = 1:ntim
          for m = 1:nfrq
            tmpdat = nan(nchan,nchan);
            indx   = find(cmbindx);
            if ~complete
              % this realizes the missing combinations to be represented as the
              % conjugate of the corresponding combination across the diagonal
              tmpdat(indx) = reshape(data.(fn{ii})(j,cmbindx(indx),m,k),[numel(indx) 1]);
              tmpdat       = ctranspose(tmpdat);
            end
            tmpdat(indx)    = reshape(data.(fn{ii})(j,cmbindx(indx),m,k),[numel(indx) 1]);
            tmpall(j,:,:,m,k) = tmpdat;
          end % for m
        end % for k
      end % for j

      % replace the data in the old representation with the new representation
      if nrpt>1,
        data.(fn{ii}) = tmpall;
      else
        data.(fn{ii}) = reshape(tmpall, [nchan nchan nfrq ntim]);
      end
    end % if numel
  end % for ii

  % remove obsolete fields
  try, data      = rmfield(data, 'powspctrm');  end
  try, data      = rmfield(data, 'labelcmb');   end
  try, data      = rmfield(data, 'dof');        end

  if ntim>1,
    data.dimord = 'chan_chan_freq_time';
  else
    data.dimord = 'chan_chan_freq';
  end

  if nrpt>1,
    data.dimord = ['rpt_',data.dimord];
  end

elseif strcmp(current, 'sparse') && strcmp(desired, 'fullfast')
  dimtok = tokenize(data.dimord, '_');
  if ~isempty(strmatch('rpt',   dimtok)), nrpt=numel(data.cumtapcnt); else nrpt = 1; end
  if ~isempty(strmatch('freq',  dimtok)), nfrq=numel(data.freq);      else nfrq = 1; end
  if ~isempty(strmatch('time',  dimtok)), ntim=numel(data.time);      else ntim = 1; end

  if ~isfield(data, 'label')
    data.label = unique(data.labelcmb(:));
  end

  nchan     = length(data.label);
  ncmb      = size(data.labelcmb,1);
  cmbindx   = zeros(nchan,nchan);

  for k = 1:size(data.labelcmb,1)
    ch1 = find(strcmp(data.label, data.labelcmb(k,1)));
    ch2 = find(strcmp(data.label, data.labelcmb(k,2)));
    if ~isempty(ch1) && ~isempty(ch2),
      cmbindx(ch1,ch2) = k;
    end
  end

  complete = all(cmbindx(:)~=0);

  fn = fieldnames(data);
  for ii=1:numel(fn)
    if numel(data.(fn{ii})) == nrpt*ncmb*nfrq*ntim;
      if nrpt==1,
        data.(fn{ii}) = reshape(data.(fn{ii}), [nrpt ncmb nfrq ntim]);
      end

      tmpall = nan(nchan,nchan,nfrq,ntim);

      for k = 1:ntim
        for m = 1:nfrq
          tmpdat = nan(nchan,nchan);
          indx   = find(cmbindx);
          if ~complete
            % this realizes the missing combinations to be represented as the
            % conjugate of the corresponding combination across the diagonal
            tmpdat(indx) = reshape(nanmean(data.(fn{ii})(:,cmbindx(indx),m,k)),[numel(indx) 1]);
            tmpdat       = ctranspose(tmpdat);
          end
          tmpdat(indx)    = reshape(nanmean(data.(fn{ii})(:,cmbindx(indx),m,k)),[numel(indx) 1]);
          tmpall(:,:,m,k) = tmpdat;
        end % for m
      end % for k

      % replace the data in the old representation with the new representation
      if nrpt>1,
        data.(fn{ii}) = tmpall;
      else
        data.(fn{ii}) = reshape(tmpall, [nchan nchan nfrq ntim]);
      end
    end % if numel
  end % for ii

  % remove obsolete fields
  try, data      = rmfield(data, 'powspctrm');  end
  try, data      = rmfield(data, 'labelcmb');   end
  try, data      = rmfield(data, 'dof');        end

  if ntim>1,
    data.dimord = 'chan_chan_freq_time';
  else
    data.dimord = 'chan_chan_freq';
  end

elseif strcmp(current, 'sparsewithpow') && strcmp(desired, 'full')
  % this is how is currently done in prepare_freq_matrices
  data = ft_checkdata(data, 'cmbrepresentation', 'sparse');
  data = ft_checkdata(data, 'cmbrepresentation', 'full');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert to new source representation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [output] = fixsource(input, varargin)

% FIXSOURCE converts old style source structures into new style source structures and the
% other way around
%
% Use as:
%   fixsource(input, type)
%    where input is a source structure,
%
% Typically, old style source structures contain
%   avg.XXX or trial.XXX fields
%
% The new style source structure contains:
%   source.pos
%   source.dim (optional, if the list of positions describes a 3D volume
%   source.XXX the old style subfields in avg/trial
%   source.XXXdimord string how to interpret the respective XXX field:
%     e.g. source.leadfield = cell(1,Npos), source.leadfielddimord = '{pos}_chan_ori'
%          source.mom       = cell(1,Npos), source.momdimord       = '{pos}_ori_rpttap'

type   = keyval('type',   varargin{:});
haspow = keyval('haspow', varargin{:});

if isempty(type),   type   = 'old'; end
if isempty(haspow), haspow = 'no';  end

fnames = fieldnames(input);
tmp    = cell2mat(strfind(fnames, 'dimord')); %get dimord like fields
if any(tmp>1),
  current = 'new';
elseif any(tmp==1),
  %don't know what to do yet data is JM's own invention
  current = 'old';
else
  current = 'old';
end

if strcmp(current, type),
  %do nothing
  output = input;

  %return
elseif strcmp(current, 'old') && strcmp(type, 'new'),
  %go from old to new

  if isfield(input, 'avg'),
    stuff  = getfield(input, 'avg');
    output = rmfield(input,  'avg');
  elseif isfield(input, 'trial'),
    stuff  = getfield(input, 'trial');
    output = rmfield(input,  'trial');
  else
    %this could occur later in the pipeline, e.g. when doing group statistics using individual subject
    %descriptive statistics
    error('the input does not contain an avg or trial field');
  end

  %-------------------------------------------------
  %remove and rename the specified fields if present
  removefields = {'xgrid';'ygrid';'zgrid';'method'};
  renamefields = {'frequency' 'freq'; 'csdlabel' 'orilabel'};
  fnames       = fieldnames(output);
  for k = 1:numel(fnames)
    ix = strmatch(fnames{k}, removefields);
    if ~isempty(ix),
      output = rmfield(output, fnames{k});
    end
    ix = strmatch(fnames{k}, renamefields(:,1), 'exact');
    if ~isempty(ix),
      output = setfield(output, renamefields{ix,2}, ...
                        getfield(output, renamefields{ix,1}));
      output = rmfield(output, fnames{k});
    end
  end

  %----------------------------------------------------------------------
  %put the stuff originally in avg or trial one level up in the structure
  fnames       = fieldnames(stuff(1));
  npos         = size(input.pos,1);
  nrpt         = numel(stuff);
  for k = 1:numel(fnames)
    if nrpt>1,
      %multiple trials
      %(or subjects FIXME not yet implemented, nor tested)
      tmp  = getfield(stuff(1), fnames{k});
      siz  = size(tmp);
      if isfield(input, 'cumtapcnt') && strcmp(fnames{k}, 'mom')
        %pcc based mom is orixrpttap
        %tranpose to keep manageable
        for kk = 1:numel(input.inside)
          indx = input.inside(kk);
          tmp{indx} = permute(tmp{indx}, [2 1 3]);
        end
        nrpttap = sum(input.cumtapcnt);
        sizvox  = [size(tmp{input.inside(1)}) 1];
        sizvox  = [nrpttap sizvox(2:end)];
      elseif strcmp(fnames{k}, 'mom'),
        %this is then probably not a frequency based mom
        nrpttap = numel(stuff);
        sizvox  = [size(tmp{input.inside(1)}) 1];
        sizvox  = [nrpttap sizvox];
      elseif iscell(tmp)
        nrpttap = numel(stuff);
        sizvox  = [size(tmp{input.inside(1)}) 1];
        sizvox  = [nrpttap sizvox];
      end

      if siz(1) ~= npos && siz(2) ==npos,
        tmp = transpose(tmp);
      end

      if iscell(tmp)
        %allocate memory for cell-array
        tmpall = cell(npos,1);
        for n = 1:numel(input.inside)
          tmpall{input.inside(n)} = zeros(sizvox);
        end
      else
        %allocate memory for matrix
        tmpall = zeros([npos nrpt siz(2:end)]);
      end

      cnt = 0;
      for m = 1:nrpt
        tmp = getfield(stuff(m), fnames{k});
        siz = size(tmp);
        if siz(1) ~= npos && siz(2) ==npos,
          tmp = transpose(tmp);
        end

        if ~iscell(tmp),
          tmpall(:,m,:,:,:) = tmp;
        else
          for n = 1:numel(input.inside)
            indx   = input.inside(n);
            tmpdat = tmp{indx};
            if isfield(input, 'cumtapcnt') && strcmp(fnames{k}, 'mom'),
              if n==1, siz1 = size(tmpdat,2); end
            else
              if n==1, siz1 = 1; end
            end
            tmpall{indx}(cnt+[1:siz1],:,:,:,:) = tmpdat;
            if n==numel(input.inside), cnt  = cnt + siz1;     end
          end
        end
      end
      output    = setfield(output, fnames{k}, tmpall);
      newdimord = createdimord(output, fnames{k}, 1);
      if ~isempty(newdimord)
        output    = setfield(output, [fnames{k},'dimord'], newdimord);
      end

    else
      tmp = getfield(stuff, fnames{k});
      siz = size(tmp);
      if isfield(input, 'cumtapcnt') && strcmp(fnames{k}, 'mom')
        %pcc based mom is orixrpttap
        %tranpose to keep manageable
        for kk = 1:numel(input.inside)
          indx = input.inside(kk);
          tmp{indx} = permute(tmp{indx}, [2 1 3]);
        end
      end
      if siz(1) ~= npos && siz(2) ==npos,
        tmp = transpose(tmp);
      end
      output    = setfield(output, fnames{k}, tmp);
      newdimord = createdimord(output, fnames{k});
      if ~isempty(newdimord)
        output    = setfield(output, [fnames{k},'dimord'], newdimord);
      end
    end
  end

  if isfield(output, 'csdlabel')
    output = setfield(output, 'orilabel', getfield(output, 'csdlabel'));
    output = rmfield(output,  'csdlabel');
  end

  if isfield(output, 'leadfield')
    % add dimord to leadfield as well. since the leadfield is not in
    % the original .avg or .trial field it has not yet been taken care of
    output.leadfielddimord = createdimord(output, 'leadfield');
  end

  if isfield(output, 'ori')
    % convert cell-array ori into matrix
    ori = zeros(3,npos) + nan;
    try,
      ori(:,output.inside) = cat(2, output.ori{output.inside});
    catch
      %when oris are in wrong orientation (row rather than column)
      for k = 1:numel(output.inside)
        ori(:,output.inside(k)) = output.ori{output.inside(k)}';
      end
     end
    output.ori = ori;
  end
  current = 'new';

elseif strcmp(current, 'new') && strcmp(type, 'old')
  %go from new to old
  error('not implemented yet');
end

if strcmp(current, 'new') && strcmp(haspow, 'yes'),

  %----------------------------------------------
  %convert mom into pow if requested and possible
  convert = 0;
  if isfield(output, 'mom') && size(output.mom{output.inside(1)},2)==1,
    convert = 1;
  else
    warning('conversion from mom to pow is not possible, either because there is no mom in the data, or because the dimension of mom>1. in that case call ft_sourcedescriptives first with cfg.projectmom');
  end

  if isfield(output, 'cumtapcnt')
    convert = 1 & convert;
  else
    warning('conversion from mom to pow will not be done, because cumtapcnt is missing');
  end

  if convert,
    npos = size(output.pos,1);
    nrpt = numel(output.cumtapcnt);
    tmpmom = cat(2,output.mom{output.inside});
    tmppow = zeros(npos, nrpt);
    tapcnt = [0;cumsum(output.cumtapcnt(:))];
    for k = 1:nrpt
      ntap = tapcnt(k+1)-tapcnt(k);
      tmppow(output.inside,k) = sum(abs(tmpmom((tapcnt(k)+1):tapcnt(k+1),:)).^2,1)./ntap;
    end
    output.pow       = tmppow;
    output.powdimord = ['pos_rpt_freq'];
  end

elseif strcmp(current, 'old') && strcmp(haspow, 'yes')
  warning('construction of single trial power estimates is not implemented here using old style source representation');

end

%--------------------------------------------------------
function [dimord] = createdimord(output, fname, rptflag);

if nargin==2, rptflag = 0; end

tmp = getfield(output, fname);

dimord = '';
dimnum = 1;
hasori = isfield(output, 'ori'); %if not, this is probably singleton and not relevant at the end

if iscell(tmp) && (size(output.pos,1)==size(tmp,dimnum) || size(output.pos,1)==size(tmp,2))
  dimord = [dimord,'{pos}'];
  dimnum = dimnum + 1;
elseif ~iscell(tmp) && size(output.pos,1)==size(tmp,dimnum)
  dimord = [dimord,'pos'];
  dimnum = dimnum + 1;
end

switch fname
  case 'cov'
    if hasori, dimord = [dimord,'_ori_ori']; end;
  case 'csd'
    if hasori, dimord = [dimord,'_ori_ori']; end;
  case 'csdlabel'
    dimord = dimord;
  case 'filter'
    dimord = [dimord,'_ori_chan'];
  case 'leadfield'
    %if hasori,
      dimord = [dimord,'_chan_ori'];
    %else
    %  dimord = [dimord,'_chan'];
    %end
  case 'mom'
    if isfield(output, 'cumtapcnt') && sum(output.cumtapcnt)==size(tmp{output.inside(1)},1)
      if hasori,
        dimord = [dimord,'_rpttap_ori'];
      else
        dimord = [dimord,'_rpttap'];
      end
    elseif isfield(output, 'time')
      if rptflag,
        dimord = [dimord,'_rpt'];
        dimnum = dimnum + 1;
      end
      if numel(output.time)==size(tmp{output.inside(1)},dimnum)
        dimord = [dimord,'_ori_time'];
      end
    end

    if isfield(output, 'freq') && numel(output.freq)>1,
      dimord = [dimord,'_freq'];
    end
  case 'nai'
    if isfield(output, 'freq') && numel(output.freq)==size(tmp,dimnum)
      dimord = [dimord,'_freq'];
    end
  case 'noise'
    if isfield(output, 'freq') && numel(output.freq)==size(tmp,dimnum)
      dimord = [dimord,'_freq'];
    end
  case 'noisecsd'
    if hasori, dimord = [dimord,'_ori_ori']; end
  case 'ori'
    dimord = '';
  case 'pow'
    if isfield(output, 'cumtapcnt') && numel(output.cumtapcnt)==size(tmp,dimnum)
      dimord = [dimord,'_rpt'];
      dimnum = dimnum + 1;
    end

    if isfield(output, 'freq') && numel(output.freq)>1 && numel(output.freq)==size(tmp,dimnum)
      dimord = [dimord,'_freq'];
      dimnum = dimnum+1;
    end

    if isfield(output, 'time') && numel(output.time)>1 && numel(output.time)==size(tmp,dimnum)
      dimord = [dimord,'_time'];
      dimnum = dimnum+1;
    end

    otherwise
      warning('skipping unknown fieldname %s', fname);
      %error(sprintf('unknown fieldname %s', fname));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = comp2raw(data)
% just remove the component topographies
data = rmfield(data, 'topo');
data = rmfield(data, 'topolabel');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = volume2source(data)
if isfield(data, 'dimord')
  % it is a modern source description
else
  % it is an old-fashioned source description
  xgrid = 1:data.dim(1);
  ygrid = 1:data.dim(2);
  zgrid = 1:data.dim(3);
  [x y z] = ndgrid(xgrid, ygrid, zgrid);
  data.pos = warp_apply(data.transform, [x(:) y(:) z(:)]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = source2volume(data)

if isfield(data, 'dimord')
  % it is a modern source description

  %this part depends on the assumption that the list of positions is describing a full 3D volume in
  %an ordered way which allows for the extraction of a transformation matrix
  %i.e. slice by slice
  try,
    if isfield(data, 'dim'),
      data.dim = pos2dim3d(data.pos, data.dim);
    else
      data.dim = pos2dim3d(data);
    end
  catch
  end
end

if isfield(data, 'dim') && length(data.dim)>=3,
  % it is an old-fashioned source description, or the source describes a regular 3D volume in pos
  xgrid = 1:data.dim(1);
  ygrid = 1:data.dim(2);
  zgrid = 1:data.dim(3);
  [x y z] = ndgrid(xgrid, ygrid, zgrid);
  ind = [x(:) y(:) z(:)];    % these are the positions expressed in voxel indices along each of the three axes
  pos = data.pos;            % these are the positions expressed in head coordinates
  % represent the positions in a manner that is compatible with the homogeneous matrix multiplication,
  % i.e. pos = H * ind
  ind = ind'; ind(4,:) = 1;
  pos = pos'; pos(4,:) = 1;
  % recompute the homogeneous transformation matrix
  data.transform = pos / ind;
end

% remove the unwanted fields
if isfield(data, 'pos'),    data = rmfield(data, 'pos');    end
if isfield(data, 'xgrid'),  data = rmfield(data, 'xgrid');  end
if isfield(data, 'ygrid'),  data = rmfield(data, 'ygrid');  end
if isfield(data, 'zgrid'),  data = rmfield(data, 'zgrid');  end

% make inside a volume
data = fixinside(data, 'logical');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = freq2raw(freq)

if strcmp(freq.dimord, 'rpt_chan_freq_time')
  dat = freq.powspctrm;
elseif strcmp(freq.dimord, 'rpttap_chan_freq_time')
  warning('converting fourier representation into raw data format. this is experimental code');
  dat = freq.fourierspctrm;
else
  error('this only works for dimord=''rpt_chan_freq_time''');
end

nrpt  = size(dat,1);
nchan = size(dat,2);
nfreq = size(dat,3);
ntime = size(dat,4);
data = [];
% create the channel labels like "MLP11@12Hz""
k = 0;
for i=1:nfreq
  for j=1:nchan
    k = k+1;
    data.label{k} = sprintf('%s@%dHz', freq.label{j}, freq.freq(i));
  end
end
% reshape and copy the data as if it were timecourses only
for i=1:nrpt
  data.time{i}  = freq.time;
  data.trial{i} = reshape(dat(i,:,:,:), nchan*nfreq, ntime);
  if any(isnan(data.trial{i}(1,:))),
    tmp = data.trial{i}(1,:);
    begsmp = find(isfinite(tmp),1, 'first');
    endsmp = find(isfinite(tmp),1, 'last' );
    data.trial{i} = data.trial{i}(:, begsmp:endsmp);
    data.time{i}  = data.time{i}(begsmp:endsmp);
  end
end
nsmp = cellfun('size',data.time,2);
seln = find(nsmp>1,1, 'first');
data.fsample = 1/(data.time{seln}(2)-data.time{seln}(1));

if isfield(freq, 'trialinfo'), data.trialinfo = freq.trialinfo; end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [data] = raw2timelock(data)

nsmp = cellfun('size',data.time,2);
data   = ft_checkdata(data, 'hastrialdef', 'yes');
ntrial = numel(data.trial);
nchan  = numel(data.label);
if ntrial==1
  data.time   = data.time{1};
  data.avg    = data.trial{1};
  data        = rmfield(data, 'trial');
  data.dimord = 'chan_time';
else
  % determine the location of the trials relative to the resulting combined time axis
  begtime  = cellfun(@min,  data.time);
  endtime  = cellfun(@max,  data.time);
  begsmp   = round((begtime - min(begtime)) * data.fsample) + 1;
  endsmp   = round((endtime - min(begtime)) * data.fsample) + 1;

  % create a combined time axis and concatenate all trials
  tmptime  = min(begtime):(1/data.fsample):max(endtime);
  tmptrial = zeros(ntrial, nchan, length(tmptime)) + nan;
  for i=1:ntrial
    tmptrial(i,:,begsmp(i):endsmp(i)) = data.trial{i};
  end

  % update the sampleinfo
  begpad = begsmp - min(begsmp);
  endpad = max(endsmp) - endsmp;
  if isfield(data, 'sampleinfo')
    data.sampleinfo = data.sampleinfo + [-begpad(:) endpad(:)];
  end

  % construct the output timelocked data
  % data.avg     = reshape(nanmean(tmptrial,     1), nchan, length(tmptime));
  % data.var     = reshape(nanvar (tmptrial, [], 1), nchan, length(tmptime))
  % data.dof     = reshape(sum(~isnan(tmptrial), 1), nchan, length(tmptime));
  data.trial   = tmptrial;
  data.time    = tmptime;
  data.dimord = 'rpt_chan_time';
  data = rmfield(data, 'fsample'); % fsample in timelock data is obsolete
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [data] = timelock2raw(data)
try
  nsmp = cellfun('size',data.time,2);
catch
  nsmp = size(data.time,2);
end
switch data.dimord
  case 'chan_time'
    data.trial{1} = data.avg;
    data.time     = {data.time};
    data          = rmfield(data, 'avg');
    seln = find(nsmp>1,1, 'first');
    data.fsample = 1/(data.time{seln}(2)-data.time{seln}(1));
  case 'rpt_chan_time'
    tmptrial = {};
    tmptime  = {};
    ntrial = size(data.trial,1);
    nchan  = size(data.trial,2);
    ntime  = size(data.trial,3);
    for i=1:ntrial
      tmptrial{i} = reshape(data.trial(i,:,:), [nchan, ntime]);
      tmptime{i}  = data.time;
    end
    data       = rmfield(data, 'trial');
    data.trial = tmptrial;
    data.time  = tmptime;
    seln = find(nsmp>1,1, 'first');
    data.fsample = 1/(data.time{seln}(2)-data.time{seln}(1));
  case 'subj_chan_time'
    tmptrial = {};
    tmptime  = {};
    ntrial = size(data.individual,1);
    nchan  = size(data.individual,2);
    ntime  = size(data.individual,3);
    for i=1:ntrial
      tmptrial{i} = reshape(data.individual(i,:,:), [nchan, ntime]);
      tmptime{i}  = data.time;
    end
    data       = rmfield(data, 'individual');
    data.trial = tmptrial;
    data.time  = tmptime;
    seln = find(nsmp>1,1, 'first');
    data.fsample = 1/(data.time{seln}(2)-data.time{seln}(1));
  otherwise
    error('unsupported dimord');
end
% remove the unwanted fields
if isfield(data, 'avg'),        data = rmfield(data, 'avg'); end
if isfield(data, 'var'),        data = rmfield(data, 'var'); end
if isfield(data, 'cov'),        data = rmfield(data, 'cov'); end
if isfield(data, 'dimord'),     data = rmfield(data, 'dimord'); end
if isfield(data, 'numsamples'), data = rmfield(data, 'numsamples'); end
if isfield(data, 'dof'),        data = rmfield(data, 'dof'); end

author: PatternRecognition

repo: OpenBMI


text:
function [assignment, cost] = assignmentsuboptimal2(distMatrix)
%ASSIGNMENTSUBOPTIMAL2    Compute suboptimal assignment
%   ASSIGNMENTSUBOPTIMAL2(DISTMATRIX) computes a suboptimal assignment
%   (minimum overall costs) for the given rectangular distance or cost
%   matrix, for example the assignment of tracks (in rows) to observations
%   (in columns). The result is a column vector containing the assigned
%   column number in each row (or 0 if no assignment could be done).
%
%   [ASSIGNMENT, COST] = ASSIGNMENTSUBOPTIMAL2(DISTMATRIX) returns the
%   assignment vector and the overall cost.
%
%   The algorithm searches the matrix for the minimum element and makes the
%   corresponding row-column assignment. After setting all elements in the
%   given row and column to infinity (i.e. forbidden assignment), the
%   search procedure is repeated until all assignments are done or only
%   infinite values are found.
%
%   This function and the corresponding mex-function can further be
%   improved by first sorting all elements instead of searching for the
%   minimum of all elements many times.
%
%   assignment.html  File Exchange  Donate via PayPal
%
%   Markus Buehren
%   Last modified 05.07.2011

% initialize
nOfRows    = size(distMatrix, 1);
assignment = zeros(nOfRows,1);
cost       = 0;

for n=1:nOfRows

  % find minimum distance observation-to-track pair
  [minDist, index1] = min(distMatrix, [], 1);
  [minDist, index2] = min(minDist);
  row = index1(index2);
  col = index2;

  if isfinite(minDist)

    % make the assignment
    assignment(row) = col;
    cost = cost + minDist;

    % delete observation-to-track pair
    distMatrix(row, :) = inf;
    distMatrix(:, col) = inf;

  else
    return
  end

end

author: GenLouvain

repo: GenLouvain


text:
% --------------------------------------------------------
% MDP Tracking
% Copyright (c) 2015 CVGL Stanford
% Licensed under The MIT License [see LICENSE for details]
% Written by Yu Xiang
% --------------------------------------------------------
%
% compute velocity
function v = compute_velocity(tracker)

fr = double(unique(tracker.frame_ids));
num = numel(fr);

% only use the past 3 frames
if num > 3
    fr = fr(num-2:num);
    num = 3;
end

% compute centers
centers = zeros(num, 2);
for i = 1:num
    index = find(tracker.frame_ids == fr(i));
    for j = 1:numel(index)
        ind = index(j);
        c = [(tracker.x1(ind)+tracker.x2(ind))/2 (tracker.y1(ind)+tracker.y2(ind))/2];
        centers(i,:) = centers(i,:) + c;
    end
    if numel(index)
        centers(i,:) = centers(i,:) / numel(index);
    end
end

count = 0;
vx = 0;
vy = 0;
cx = centers(:,1);
cy = centers(:,2);
for j = 2:num
    vx = vx + (cx(j)-cx(j-1)) / (fr(j) - fr(j-1));
    vy = vy + (cy(j)-cy(j-1)) / (fr(j) - fr(j-1));
    count = count + 1;
end
if count
    vx = vx / count;
    vy = vy / count;
end
v = [vx, vy];

author: yuxng


text:
function [AtA,A] = corrMatrix3D(obj)
% calucate 3D correlation matrix
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

nCha = size(obj.kCalib,4);

if(isreal(obj.kCalib))
    A = zeros(prod(obj.calibSize - obj.kernelSize + 1), prod(obj.kernelSize)*nCha,obj.measPara.precision);
else
    A = complex(zeros(prod(obj.calibSize - obj.kernelSize + 1), prod(obj.kernelSize)*nCha,obj.measPara.precision),zeros(prod(obj.calibSize - obj.kernelSize + 1), prod(obj.kernelSize)*nCha,obj.measPara.precision));
end
% A = [];
counter = 1;
for n=1:nCha
    if(isreal(obj.kCalib(:,:,:,n)))
        if(strcmp(obj.measPara.precision,'single'))
            A(:,counter:counter+prod(obj.kernelSize)-1) = im3colRSingle(obj.kCalib(:,:,:,n),obj.kernelSize).';
        else
            A(:,counter:counter+prod(obj.kernelSize)-1) = im3colR(obj.kCalib(:,:,:,n),obj.kernelSize).'; % before: tmp =
        end
    else
        if(strcmp(obj.measPara.precision,'single'))
            A(:,counter:counter+prod(obj.kernelSize)-1) = im3colCSingle(obj.kCalib(:,:,:,n),obj.kernelSize).';
        else
            A(:,counter:counter+prod(obj.kernelSize)-1) = im3colC(obj.kCalib(:,:,:,n),obj.kernelSize).';
        end
    end
    counter = counter + prod(obj.kernelSize);
% 	A = [A, tmp];
end

AtA = A'*A;

end

author: thomaskuestner


text:
 function [xs, info] = ir_pwls_os_rlalm(x, Ab, yi, R, varargin)
%function [xs, info] = ir_pwls_os_rlalm(x, Ab, yi, R, [options])
%|
%| penalized weighted least squares estimation / image reconstruction
%| using relaxed linearized augmented lagrangian method with
%| (optionally relaxed) ordered subsets. (relaxed OS-LALM)
%|
%| See ?. 201? IEEE T-MI by Hung Nien & J A Fessler
%| "Relaxed linearized algorithms for faster X-ray CT image reconstruction"
%|	http://dx.doi.org/10.1109/TMI.201?
%|
%|
%| cost(x) = (y-Ax)' W (y-Ax) / 2 + R(x)
%|
%| in
%|	x	[np 1]		initial estimate
%|	Ab	[nd np]		Gblock object (needs abs(Ab) method)
%|				or sparse matrix (implies nsubset=1)
%|	yi	[nb na]		measurements (noisy sinogram data)
%|	R	penalty		object (see Reg1.m), can be []
%|
%| option
%|	niter			# of iterations (default: 1)
%|	wi	[nb na]		weighting sinogram (default: [] for uniform)
%|	pixmax	[1] or [2]	max pixel value, or [min max] (default [0 inf])
%|	denom	[np 1]		precomputed denominator
%|	aai	[nb na]		precomputed row sums of |Ab|
%|	relax0	[1] or [2]	relax0 or (relax0, relax_rate) (default 1)
%|	rho			AL penalty parameter (default: [] for decreasing rho)
%|	alpha			over-relaxation parameter (default: 1.999)
%|	userfun	@		user-defined function handle (see default below)
%|					taking arguments (x, userarg{:})
%|	userarg	{}		user arguments to userfun (default {})
%|	chat
%|
%| out
%|	xs	[np niter]	iterates
%|	info	[niter 1]	time
%|
%| 2015-11-30, Hung Nien, based on ir_pwls_os_lalm
%| 2015-11-30, tweaks by Jeff Fessler

if nargin == 1 && streq(x, 'test'), ir_pwls_os_rlalm_test, return, end
if nargin < 4, ir_usage, end

% defaults
arg.niter = 1;
arg.isave = [];
arg.userfun = @userfun_default;
arg.userarg = {};
arg.pixmax = inf;
arg.chat = false;
arg.wi = [];
arg.aai = [];
arg.rho = []; % default: decreasing rho
arg.alpha = 1.999;
arg.relax0 = 1;
arg.denom = [];
arg.scale_nblock = true; % traditional scaling
arg.update_even_if_denom_0 = true;
arg = vararg_pair(arg, varargin);

arg.isave = iter_saver(arg.isave, arg.niter);

Ab = block_op(Ab, 'ensure'); % make it a block object (if not already)
nblock = block_op(Ab, 'n');
starts = subset_start(nblock);

cpu etic

wi = arg.wi;
if isempty(wi)
	wi = ones(size(yi), class(yi));
end
if isempty(arg.aai) && isempty(arg.denom)
	arg.aai = reshape(sum(abs(Ab)'), size(yi)); % a_i = sum_j |a_ij|
end

% check input sinogram sizes for OS
if (ndims(yi) ~= 2) || (size(yi,2) == 1 && nblock > 1)
	fail 'bad yi size'
end
if (ndims(wi) ~= 2) || (size(wi,2) == 1 && nblock > 1)
	fail 'bad wi size'
end

relax0 = arg.relax0(1);
if length(arg.relax0) == 1
	relax_rate = 0;
elseif length(arg.relax0) == 2
	relax_rate = arg.relax0(2);
else
	error relax
end

if length(arg.pixmax) == 2
	pixmin = arg.pixmax(1);
	pixmax = arg.pixmax(2);
elseif length(arg.pixmax) == 1
	pixmin = 0;
	pixmax = arg.pixmax;
else
	error pixmax
end

% likelihood denom, if not provided
denom = arg.denom;
if isempty(denom)
	denom = abs(Ab)' * col(arg.aai .* wi); % needs abs(Ab)
        arg = rmfield(arg, 'aai'); % clear to save memory - not needed below
end
if ~arg.update_even_if_denom_0
	% todo: this may not work for LALM because "denom" appears in numerator!
	denom(denom == 0) = inf; % trick: prevents pixels where denom=0 being updated
end

if isempty(R)
	pgrad = 0; % unregularized default
	Rdenom = 0;
end

alpha = arg.alpha;
if alpha<1 || alpha>2
	fail 'alpha should be between 1 and 2'
end

rho = arg.rho;
if isempty(rho)
	rho = @(k) pi/(alpha*k) * sqrt(1 - (pi/(2*(alpha*k)))^2) * (k>1) + (k==1);
else
	rho = @(k) rho; % constant user-specified value
end

[nb na] = size(yi);

x = x(:);
np = length(x);
xs = zeros(np, length(arg.isave), 'single');
if any(arg.isave == 0)
	xs(:, arg.isave == 0) = single(x);
end

%info = zeros(niter,?); % do not initialize since size may change

% initialization
iblock = nblock;
ia = iblock:nblock:na;
li = Ab{iblock} * x;
li = reshape(li, nb, length(ia));
resid = wi(:,ia) .* (li - yi(:,ia));
if arg.scale_nblock
	scale = nblock; % traditional way
else
	scale = na / numel(ia); % alternative - untested
end
zeta = scale * Ab{iblock}' * resid(:);

g = rho(1) * zeta;
h = denom .* x - zeta;

% iterate
for iter = 1:arg.niter
	ticker(mfilename, iter, arg.niter)

	relax = relax0 / (1 + relax_rate * (iter-1));

	% loop over subsets
	for iset = 1:nblock
		k = nblock*(iter-1)+iset;

		num = rho(k) * (denom .* x - h) + (1-rho(k)) * g;
		den = rho(k) * denom;
		if ~isempty(R)
			num = num + R.cgrad(R, x);
			den = den + R.denom(R, x);
		end
		x = x - relax * num ./ den;
		x = max(x, pixmin);
		x = min(x, pixmax);

		iblock = starts(iset);
		ia = iblock:nblock:na;

		li = Ab{iblock} * x;
		li = reshape(li, nb, length(ia));
		resid = wi(:,ia) .* (li - yi(:,ia));

		if arg.scale_nblock
			scale = nblock; % traditional way
		else
			scale = na / numel(ia); % alternative - untested
		end

		zeta = scale * Ab{iblock}' * resid(:); % A' * W * (y - A*x)
		g = (rho(k) * (alpha * zeta + (1-alpha) * g) + g) / (rho(k)+1);
		h = alpha * (denom .* x - zeta) + (1-alpha) * h;
	end

	if any(arg.isave == iter)
		xs(:, arg.isave == iter) = single(x);
	end
	info(iter,:) = arg.userfun(x, arg.userarg{:});
end

% default user function.
% using this evalin('caller', ...) trick, one can compute anything of interest
function out = userfun_default(x, varargin)
chat = evalin('caller', 'arg.chat');
if chat
%	x = evalin('caller', 'x');
	printm('minmax(x) = %g %g', min(x), max(x))
end
out = cpu('etoc');

% ir_pwls_os_rlalm_test()
% test with small quadratic case 1/2 |y - A x|_2^2 + 1/2 |C x|_2^2
% see also ir_ct_fan_beam_sqs_vs_lalm.m
function ir_pwls_os_rlalm_test
rng(0)
nd = 99; np = 38;
A = randn(nd,np);
A = abs(A);
R = Reg1(ones(np,1), 'type_penal', 'mat', 'type_diff', 'sparse', 'beta', 2^0);
C = R.C;
C = full(C);
if 0 % try to make it so unit step size is reasonable?
	D = abs(A)' * (abs(A) * ones(np,1)) + abs(C)' * (abs(C) * ones(np,1));
	A = A * diag(sqrt(1 ./ D));
	C = C * diag(sqrt(1 ./ D));
end
H = A' * A + C' * C;
pr cond(H)

xtrue = ones(np,1); xtrue(10:20) = 2;
yi = A * xtrue; yi = yi + 0.01 * max(yi(:)) * randn(nd,1);
xhat = H \ (A' * yi);
% nrms(xhat, xtrue)

x0 = zeros(np,1);
niter = 100;

printm 'Nesterov' % todo: replace with OGM
gradfun = @(x) A' * (A * x - yi) + C' * (C * x);
step = abs(A') * (abs(A) * ones(np,1)) + abs(C)' * (abs(C) * ones(np,1));
step = 1 / max(step); % todo: use D
xf = ir_fista(x0, 'gradfun', gradfun, 'step', step, 'niter', niter, ...
	'isave', 'all');

printm 'Relaxed LALM'
xl = ir_pwls_os_rlalm(x0, Gmatrix(A), yi, R, 'niter', niter, 'isave', 'all');

printm 'CG'
xc = qpwls_pcg1(x0, A, 1, yi, R.C, 'niter', niter, 'isave', 'all');

if im
	im plc 1 2
	im subplot 1
	plot([xtrue xhat xl(:,end) xc(:,end) xf(:,end)])
	im subplot 2
	xlabel 'iteration', ylabelf 'NRMSD $x^n$ vs $\hat{x}$'
	plot(	0:niter, nrms(xl, xhat), 'r-', ...
		0:niter, nrms(xf, xhat), 'g-', ...
		0:niter, nrms(xc, xhat), 'c-')
	legend('Relaxed LALM', 'FGM', 'CG')
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function str = get_structure( name )
%STRUCTURE Summary of this function goes here
%
% Structure used in this toolbox.
%
% individual structure:
% parameter: the parameter space point of the individual. it's a column-wise
% vector.
% objective: the objective space point of the individual. it's column-wise
% vector. It only have value after evaluate function is called upon the
% individual.
% estimation: Also a structure array of the individual. It's not used in
% MOEA/D but used in MOEA/D/GP. For every objective, the field contains the
% estimation from the GP model.
%
% estimation structure:
% obj: the estimated mean.
% std: the estimated standard deviation for the mean.
%
% subproblem structure:
% weight: the decomposition weight for the subproblem.
% optimal: the current optimal value of the current structure.
% curpoiont: the current individual of the subproblem.
% optpoint: the point that gain the optimal on the subproblem.
%

switch name
    case 'individual'
        str = struct('parameter',[],'objective'[],'estimation'[]);
    case 'subproblem'
        str = struct('weight',[],'optimal',[],'curpoint',[],'optpoint',[]);
    case 'estimation'
        str = struct();
    otherwise
end

author: sfvsfv

repo: Mathematical-modeling


text:
function area = triangleArea(pt1, pt2, pt3)
%TRIANGLEAREA Signed area of a triangle
%
%   AREA = triangleArea(P1, P2, P3)
%   Computes area of the triangle whose vertices are given by P1, P2 and
%   P3. Each vertex is a 1-by-2 row vector.
%
%   AREA = triangleArea(PTS)
%   Concatenates vertex coordinates in a 3-by-2 array. Each row of the
%   array contains coordinates of one vertex.
%
%
%   Example
%   % Compute area of a Counter-Clockwise (CCW) oriented triangle
%     triangleArea([10 10], [30 10], [10 40])
%     ans =
%         300
%
%   % Compute area of a Clockwise (CW) oriented triangle
%     triangleArea([10 40], [30 10], [10 10])
%     ans =
%         -300
%
%   See also
%   polygonArea, triangleArea3d
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-08-23,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% if data is given as one array, split vertices
if nargin == 1
    pt2 = pt1(2,:);
    pt3 = pt1(3,:);
    pt1 = pt1(1,:);
end

% compute individual vectors
v12 = bsxfun(@minus, pt2, pt1);
v13 = bsxfun(@minus, pt3, pt1);

% compute area from cross product
area = (v13(:,2) .* v12(:,1) - v13(:,1) .* v12(:,2)) / 2;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
classdef nme_load3p < mp.nm_element & mp.form_acp

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end     %% properties

    methods
        function name = name(obj)
            name = 'load3p';
        end

        function np = np(obj)
            np = 3;     %% this is a 3 port element
        end

        function obj = build_params(obj, nm, dm)
            build_params@mp.nm_element(obj, nm, dm);    %% call parent
            dme = obj.data_model_element(dm);

            %% constant complex power demand
            pd = [dme.pd1 dme.pd2 dme.pd3];
            qd = pd .* tan(acos( [dme.pf1 dme.pf2 dme.pf3] ));

            obj.s = pd(:) + 1j * qd(:);
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower