function map= bktr_mapping(GSM,DM,map)

% in:
% GSM:    COBRA genome scale model
% DM:     COBRA dynamic model
% map:    optional, pratly filled map

% out:
% map:    structure containing the mapping of DM to GSM
% map.DRxns:   Dynamic Rxn Ids
% map.GSMRxns: Mappes IDs of the GSM (can be more than one for each dyn Rxn)
% map.DMets:   Dynamic Metabolite IDs
% map.GSMMets: Mapped Ids of the GSM (one to one)

%% Initialize
if nargin==2
    map.DRxns=DM.rxns;
    map.GSMRxns=cell(1,length(DM.rxns));
    map.DMets=DM.mets;
    map.GSMMets=cell(1,length(DM.mets));
    map.complete = false;
else
    disp('Mapped Metabolites:')
    for k=reshape(find(~cellfun(@isempty, map.GSMMets)),1,length(find(~cellfun(@isempty, map.GSMMets))))
        
        if 1==length(map.GSMMets(k))
            if strcmp(map.GSMMets{k},'rm')
                disp(['Metabolite ' DM.metNames{k} ' (' DM.mets{k} ') was removed.'])
                continue
            end
            
        end
        disp(['Metabolite ' DM.metNames{k} ' (' DM.mets{k} ') mapped to ' GSM.metNames{map.GSMMets{k}} ' (' GSM.mets{map.GSMMets{k}} ')'])
        
    end
    fprintf('\n')
end

%% Map Metabolites
idx_met_ChEBIID = 0;
idx_met_KEGGID = 0;
idx_met_Formula = 0;



if ~isempty(find(cellfun(@isempty, map.GSMMets)))
todos= reshape(find(cellfun(@isempty, map.GSMMets)),1,length(find(cellfun(@isempty, map.GSMMets)))); % match all metabolites that were not matched yet;

    for k=todos
        % by ChEBIID
        if all([isfield(DM,'metChEBIID') , isfield(GSM,'metChEBIID')])
            if ~isempty(DM.metChEBIID{k})
                a=ismember(GSM.metChEBIID,[DM.metChEBIID{k}]);
                idx_met_ChEBIID=find(a);
            end
        end
        % by KEGGID
        if all([isfield(DM,'metKEGGID') , isfield(GSM,'metKEGGID') ])
            if ~isempty(DM.metKEGGID{k})
                a=ismember(GSM.metKEGGID,[DM.metKEGGID{k}]);
                idx_met_KEGGID=find(a);
            end
        end
        % by Formula
        if all([isfield(DM,'metFormula') , isfield(GSM,'metFormula') ])
            if ~isempty(DM.metFormula{k})
                a=ismember(GSM.metFormula,[DM.metFormula{k}]);
                 idx_met_Formula=find(a);
            end
        end
        % by metID
        [a, idx_met_ID]=ismember([DM.mets{k }],GSM.mets);
        idx_met_ID(~a)=[];
        
        
        
        all_matches = unique([idx_met_ChEBIID(:); idx_met_KEGGID(:); idx_met_Formula(:); idx_met_ID(:)]);
        all_matches(all_matches==0)=[];
        if isempty(all_matches) % no match found
            disp(['No match found for Metabolite ' DM.metNames{k} ' (' DM.mets{k} ')'])
            selected_match=input('Input ID (field model.mets from the genome scale model) or remove (rm) cancel (N): ','s');
            if strcmp(selected_match,'N')
                disp('Mapping failed (no match)! Please revisit dynmaic model file!')
                return
            elseif ismember(selected_match,GSM.mets)
                selected_match=find(strcmp(selected_match,GSM.mets));
                disp(['Metabolite ' DM.metNames{k} ' (' DM.mets{k} ') mapped to ' GSM.metNames{selected_match} ' (' GSM.mets{selected_match} ')'])
                map.GSMMets{k} = selected_match; %store in the list
            elseif strcmp(selected_match, 'rm') % Tag fo removal (e.g. sums of  compounds)
                map.GSMRxns{k} = 'rm';
                disp(['Metabolite ' DM.metNames{k} ' (' DM.mets{k} ') was removed.'])
            else
                disp('Mapping failed (input metabolite not in model)! Please revisit dynmaic model file!')
                return
            end
            
        elseif length(all_matches)==1 % unique match found
            disp(['Unique match found for Metabolite ' DM.metNames{k} ' (' DM.mets{k} '): ' GSM.metNames{all_matches} ' (' GSM.mets{all_matches} ')'])
            selected_match=input('Agreed? Yes (Y) or choose other metabolite by ID (field model.mets from the genome scale model) or remove (rm) cancel (N): ','s');
            if strcmp(selected_match,'Y')
                map.GSMMets{k} = all_matches; %store in the list
                disp('Mappe unique match as above!')
            elseif strcmp(selected_match,'N')
                disp('Mapping failed (no match)! Please revisit dynmaic model file!')
                return
            elseif ismember(selected_match,GSM.mets)
                selected_match=find(strcmp(selected_match,GSM.mets));
                disp(['Metabolite ' DM.metNames{k} ' (' DM.mets{k} ') mapped to ' GSM.metNames{selected_match} ' (' GSM.mets{selected_match} ')'])
                map.GSMMets{k} = selected_match; %store in the list
            elseif strcmp(selected_match, 'rm') % Tag fo removal (e.g. sums of  compounds)
                map.GSMRxns{k} = 'rm';
                disp(['Metabolite ' DM.metNames{k} ' (' DM.mets{k} ') was removed.'])
            else
                disp('Mapping failed (input metabolite not in model)! Please revisit dynmaic model file!')
                return
            end
            
            
        else % more than one match found (selection required)
            disp(['Severeal matches found for Metabolite ' DM.metNames{k} ' (' DM.mets{k} '): '])
            if isfield(GSM,'metChEBIID'), matchChEBIID = [DM.metChEBIID(k); GSM.metChEBIID(all_matches)]; else matchChEBIID = cell(size(all_matches)+[0,1]);end
            if isfield(GSM,'metKEGGID'),   matchKEGGID = [DM.metKEGGID(k); GSM.metKEGGID(all_matches)]; else matchKEGGID = cell(size(all_matches)+[0,1]);end
            if isfield(GSM,'metFormula'), matchFormula = [DM.metFormula(k);GSM.metFormula(all_matches)]; else matchFormula = cell(size(all_matches)+[0,1]);end
            
            disp(table(['Dyn Met',num2cell(1:length(all_matches))]', ...
                [DM.metNames{k} ;GSM.metNames(all_matches)],...
                [DM.mets{k}; GSM.mets(all_matches)], ...
                matchChEBIID,matchKEGGID,matchFormula,...
                'VariableNames',{'MatchNumber' ,'MetName', 'ID', 'ChEBIID', 'KEGGID', 'Formula'}))
            
            selected_match=input('Please select (#, comma separated list possible): ','s');
            selected_match=str2double(strsplit(selected_match,','));
            
            while ~all(ismember(selected_match,1:length(all_matches)))
                selected_match=input('Wrong input, try again to select (#, comma separated list possible) or choose other metabolite by ID (field model.mets from the genome scale model)or remove (rm) or cancel (N): ','s');
                
                if strcmp(selected_match,'N')
                    disp('Mapping failed! Please revisit dynmaic model file!')
                    return
                elseif strcmp(selected_match,'rm')
                    map.GSMMets{k} = 'rm';
                    break
                elseif ismember(selected_match,GSM.mets)
                    break
                end
                %             selected_match=str2double(strsplit(selected_match,','));
            end
            
            if strcmp(selected_match,'rm')
                disp(['Metabolite ' DM.metNames{k} ' (' DM.mets{k} ') was removed.'])
            elseif ismember(selected_match,1:length(all_matches))
                disp(['Metabolite ' DM.metNames{k} ' (' DM.mets{k} ') mapped to ' GSM.metNames{all_matches(selected_match)} ' (' GSM.mets{all_matches(selected_match)} ')'])
                map.GSMMets{k} = all_matches(selected_match(:))'; %store in the list
            else ismember(selected_match,GSM.mets)
                selected_match=find(strcmp(selected_match,GSM.mets));
                disp(['Metabolite ' DM.metNames{k} ' (' DM.mets{k} ') mapped to ' GSM.metNames{selected_match} ' (' GSM.mets{selected_match} ')'])
                map.GSMMets{k} = selected_match; %store in the list
            end
        end
    end
end
fprintf('\n')
fprintf('Metabolite matching sucessfull! \n\n')

%% Map Reactions
if nargin > 2
    disp('Mapped Reactions:')
    
    for k=reshape(find(~cellfun(@isempty, map.GSMRxns)),1,length(find(~cellfun(@isempty, map.GSMRxns))))
        if 1==length(map.GSMRxns(k))
            if strcmp(map.GSMRxns{k},'rm')
                disp(['Reaction ' DM.rxnNames{k} ' (' DM.rxns{k} ') was removed.'])
                continue
            end
        end
        
        disp(['Reaction ' DM.rxnNames{k} ' (' DM.rxns{k} ') mapped to ' GSM.rxnNames{map.GSMRxns{k}} ' (' GSM.rxns{map.GSMRxns{k}} ')'])
        
    end
    fprintf('\n')
end


% pull all GSM rections from matched metabolites
pull_rxn_list = findRxnsFromMets(GSM, GSM.mets([map.GSMMets{:}]));

if ~isempty(find(cellfun(@isempty, map.GSMRxns)))
    todos = reshape(find(cellfun(@isempty, map.GSMRxns)),1, length(find(cellfun(@isempty, map.GSMRxns))));

    for k = todos
        % mapping by stoichiometry (exact match)
        part_mets = [map.GSMMets{find(DM.S(:,k))}];% find mapped substrates and products of current rxn
        idx_rxn_stoich = find(all(GSM.S(part_mets,:),1) ) ; % look for reaction in GSM with the same/comparable stoichiometry
        
        idx_rxn_genes=[];
        try % mapping by genes (if both models are corrrectly annotaded with genes)
            Kth_genes = findGenesFromRxns(DM,DM.rxns(k));
            [~, idx_rxn_genes]=findRxnsFromGenes(GSM,[Kth_genes{:}],0,1); % find all rxns in the GSm catalyzed by the genes of the kth rxn of DM
            idx_rxn_genes = unique(idx_rxn_genes(:,1));
            idx_rxn_genes=idx_rxn_genes(ismember(idx_rxn_genes, pull_rxn_list)); % check if they utilize the mapped metabolites
        catch
        end
        
             
        % mapping by Kegg IDs
        idx_rxn_KEGGID=0;
        if all([isfield(DM,'rxnKEGGID') , isfield(GSM,'rxnKEGGID')])
            if ~isempty(DM.rxnKEGGID{k}) 
                try
                a=ismember(GSM.rxnKEGGID,[DM.rxnKEGGID{k}]);
                idx_rxn_KEGGID=find(a);
                catch 
                end
            end
        end
        
        % Mapping by E.C. numbers
        idx_rxn_ECNumbers=0;
        if all([isfield(DM,'rxnECNumbers'), isfield(GSM,'rxnECNumbers')])
            if ~isempty(DM.rxnECNumbers{k}),
                try
                a=ismember(GSM.rxnECNumbers,[DM.rxnECNumbers{k}]);
                idx_rxn_ECNumbers=find(a);
                catch
                end
            end
        end
        
        
        %%% Select matches
        all_matches = unique([idx_rxn_stoich(:); idx_rxn_KEGGID(:); idx_rxn_genes(:); idx_rxn_ECNumbers(:)]);
        all_matches(all_matches==0)=[];
        
        if isempty(all_matches) % no match found
            disp(['No match found for Reaction ' DM.rxnNames{k} ' (' DM.rxns{k} ')'])
            selected_match=input('Input ID (field model.rxns from the genome scale model, comma separated list possible) remove reaction (rm) or cancel (N): ','s');
            selected_match= strsplit(selected_match,{',',' ','\t'});
            if strcmp(selected_match,'N')
                disp('Mapping failed (no match)! Please revisit dynmaic model file!')
                return
            elseif all(ismember(selected_match,GSM.rxns))
                for i=1:length(selected_match)
                    selected_match{i}=find(strcmp(selected_match(i),GSM.rxns));
                end
                
                disp(['Reaction ' DM.rxnNames{k} ' (' DM.rxns{k} ') mapped to ' GSM.rxnNames{[selected_match{:}]} ' (' GSM.rxns{[selected_match{:}]} ')'])
                map.GSMRxns{k} = cell2mat(selected_match); %store in the list
            elseif strcmp(selected_match, 'rm') % Tag fo removal (biomass reaction or consumer pseudo rxn)
                map.GSMRxns{k} = 'rm';
                disp(['Reaction ' DM.rxnNames{k} ' (' DM.rxns{k} ') was removed.'])
            else
                idx_not_found=~ismember(selected_match,GSM.rxns);
                disp(['Mapping failed (input reaction ' [selected_match{idx_not_found}] ' not in model)! Please revisit dynamic model file!'])
                return
            end
            
        elseif length(all_matches)==1 % unique match found
            disp(['Unique match found for Reaction ' DM.rxnNames{k} ' (' DM.rxns{k} '): ' GSM.rxnNames{all_matches} ' (' GSM.rxns{all_matches} ')'])
            selected_match=input('Agreed? Yes (Y) or choose other reaction by ID (field model.rxns from the genome scale model) or remove (rm) cancel (N): ','s');
            if strcmp(selected_match,'Y')
                map.GSMRxns{k} = all_matches; %store in the list
                disp('Mapped unique match as above!')
            elseif strcmp(selected_match,'N')
                disp('Mapping failed (no match)! Please revisit dynmaic model file!')
                return
            elseif strcmp(selected_match, 'rm') % Tag fo removal (e.g. sums of  compounds)
                map.GSMRxns{k} = 'rm';
                disp(['Reaction ' DM.rxnNames{k} ' (' DM.rxns{k} ') was removed.'])
            else

               selected_match= strsplit(selected_match,{',',' ','\t'});
                if all(ismember(selected_match,GSM.rxns))
                    
                    selected_match=ismember(selected_match,GSM.rxns);                    
                    disp(['Reaction ' DM.rxnNames{k} ' (' DM.rxns{k} ') mapped to ' GSM.rxnNames{[selected_match{:}]} ' (' GSM.rxns{[selected_match{:}]} ')'])
                    map.GSMRxns{k} = selected_match(:); %store in the list
                else
                    disp('Mapping failed (input metabolite not in model)! Please revisit dynmaic model file!')
                    return
                end
            end
            
            
            
        else % more than one match found (selection required)
            disp(['Severeal matches found for Reaction ' DM.rxnNames{k} ' (' DM.rxns{k} '): '])
            if isfield(GSM,'rxnKEGGID') && isfield(DM,'rxnKEGGID'),   matchKEGGID = [DM.rxnKEGGID(k); GSM.rxnKEGGID(all_matches)]; else matchKEGGID = cell(size(all_matches)+[1,0]);end
            if isfield(GSM,'rxnECNumber') && isfield(DM,'rxnECNumber'), matchECNumber = [DM.rxnECNumber(k);GSM.rxnECNumber(all_matches)]; else matchECNumber = cell(size(all_matches)+[1,0]);end
            if isfield(GSM,'grRules') && isfield(DM,'grRules'), matchGenes = [DM.grRules(k);GSM.grRules(all_matches)]; else matchGenes = cell(size(all_matches)+[1,0]);end
            
            
            
            table(['DynRxn',num2cell(1:length(all_matches))]', ...
                [DM.rxnNames{k} ;GSM.rxnNames(all_matches)],...
                [DM.rxns{k}; GSM.rxns(all_matches)], ...
                matchKEGGID,matchECNumber, matchGenes,...
                'VariableNames',{'MatchNumber' ,'RxnName', 'ID', 'KEGGID', 'ECNumber', 'Genes'}),
            selected_match=input('Please select (#, comma separated list possible): ','s');
            selected_match=str2double(strsplit(selected_match,','));
            while ~all(ismember(selected_match,1:length(all_matches)))
                selected_match=input('Wrong input, try again to select (#) or remove (rm) or cancel (N): ','s');
                if strcmp(selected_match,'N')
                    disp('Mapping failed! Please revisit dynmaic model file!')
                    return
                elseif strcmp(selected_match,'rm')
                    map.GSMRxns{k} ='rm';
                    break
                end
                selected_match=str2double(strsplit(selected_match,','));
            end
            
            if strcmp(selected_match,'rm')
                disp(['Reaction ' DM.rxnNames{k} ' (' DM.rxns{k} ') was removed.'])
            else
                disp(['Reaction ' DM.rxnNames{k} ' (' DM.rxns{k} ') mapped to ' GSM.rxnNames{all_matches(selected_match)} ' (' GSM.rxns{all_matches(selected_match)} ')'])
                map.GSMRxns{k} = all_matches(selected_match); %store in the list
            end
        end
        
        
    end
end
fprintf('\n')
fprintf('Reaction matching sucessfull! \n\n')


% map.complete=true;