text:
function [ grid_order, grid_point ] = cc_grids_minmax ( dim_num, q_min, ...
  q_max, grid_num, point_num )

%*****************************************************************************80
%
%% CC_GRIDS_MINMAX computes CC orders and grids with Q_MIN <= Q <= Q_MAX.
%
%  Discussion:
%
%    The necessary dimensions of GRID_ORDER and GRID_POINT can be
%    determined by calling CC_GRIDS_MINMAX first.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 October 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer Q_MIN, Q_MAX, the minimum and maximum values of
%    Q, the sum of the orders in each spatial coordinate.
%
%    Input, integer GRID_NUM, the number of Clenshaw Curtis
%    grids whose Q value is between Q_MIN and Q_MAX.
%
%    Input, integer POINT_NUM, the total number of points in the grids.
%
%    Output, integer GRID_ORDER(DIM_NUM,GRID_NUM), contains, for each
%    grid, the order of the Clenshaw-Curtis rule in each dimension.
%
%    Output, real GRID_POINT(DIM_NUM,POINT_NUM), contains
%    a list of all the abscissas of all the rules, listed one grid at
%    a time.  If a point occurs in several grids, it will be listed
%    several times.
%

%
%  Outer loop generates Q's from Q_MIN to Q_MAX.
%
  point_num = 0;
  grid_num = 0;

  for q = q_min : q_max
%
%  Middle loop generates next partition that adds up to Q.
%
    more = 0;
    order_1d = [];

    while ( 1 )

      [ order_1d, more ] = compnz_next ( q, dim_num, order_1d, more );
%
%  Inner (hidden) loop generates all CC points corresponding to given grid.
%
      order_nd = prod ( order_1d(1:dim_num) );

      grid_point(1:dim_num,point_num+1:point_num+order_nd) = cc_grid ( ...
        dim_num, order_1d, order_nd );

      point_num = point_num + order_nd;

      grid_num = grid_num + 1;
      grid_order(1:dim_num,grid_num) = order_1d(1:dim_num);

      if ( ~more )
        break
      end

    end

  end

  return
end


text:
[data]=xlsread('LRP.xls');
datevector=[1972:0.25:2014.75]';
subplot(2,2,1)
plot(datevector,data(:,4),'r--','LineWidth',0.5);
hold on
plot(datevector,data(:,5),'r','LineWidth',2.5);
hold on
plot(datevector,data(:,6),'r--','LineWidth',0.5);
hold on
plot(datevector,data(:,13),'b--','LineWidth',0.5);
hold on
plot(datevector,data(:,14),'b','LineWidth',2.5);
hold on
plot(datevector,data(:,15),'b--','LineWidth',0.5);
hold on
plot(datevector,data(:,1),'k','LineWidth',2.8);
hold off
% label the endogenous variables
axis tight
title('log GDP')
box off
subplot(2,2,2)
plot(datevector,data(:,7),'r--','LineWidth',0.5);
hold on
plot(datevector,data(:,8),'r','LineWidth',2.5);
hold on
plot(datevector,data(:,9),'r--','LineWidth',0.5);
hold on
plot(datevector,data(:,16),'b--','LineWidth',0.5);
hold on
plot(datevector,data(:,17),'b','LineWidth',2.5);
hold on
plot(datevector,data(:,18),'b--','LineWidth',0.5);
hold on
plot(datevector,data(:,2),'k','LineWidth',2.8);
hold off
% label the endogenous variables
axis tight
title('Inflation')
%set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
box off
subplot(2,2,3)
plot(datevector,data(:,10),'r--','LineWidth',0.5);
hold on
h1=plot(datevector,data(:,11),'r','LineWidth',2.5);
hold on
plot(datevector,data(:,12),'r--','LineWidth',0.5);
hold on
plot(datevector,data(:,19),'b--','LineWidth',0.5);
hold on
h2=plot(datevector,data(:,20),'b','LineWidth',2.5);
hold on
plot(datevector,data(:,21),'b--','LineWidth',0.5);
hold on
h3=plot(datevector,data(:,3),'k','LineWidth',2.8);
hold off
% label the endogenous variables
axis tight
title('Effective Federal Funds Rate')
legend([h1 h2 h3],{'Prior for the long run','Normal Wishart prior','Actual data'});
box off

ncolumns=ceil(n^0.5);
nrows=ceil(n/ncolumns);
dcNW=figure;
for i=1:n;
set(dcNW,'name','Deterministic component');
subplot(nrows,ncolumns,i)
plot(decimaldates1,median(hd_record{n^2+i}),'r','LineWidth',2.8);
hold on
plot(decimaldates1,quantile(hd_record{n^2+i},0.16),'b--','LineWidth',1.8);
hold on
plot(decimaldates1,quantile(hd_record{n^2+i},0.84),'b--','LineWidth',1.8);
hold on
plot(decimaldates1,data_endo(5:end,i),'k','LineWidth',1.8);
hold off
% label the endogenous variables
axis tight
title(endo{i,1})
set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
box off
end

author: european-central-bank

repo: BEAR-toolbox


text:
% Copyright 2019 Jonas Koenemann, Moritz Diehl, University of Freiburg
% Redistribution is permitted under the 3-Clause BSD License terms. Please
% ensure the above copyright notice is visible in any derived work.
%
function stairs(x, y, varargin)
  ocl.utils.checkStartup()

  x = ocl.Variable.getValue(x);
  y = ocl.Variable.getValue(y);

  stairs(x,y,'LineWidth', 3, varargin{:})
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function model_wrl = GB_read_vrml_fast_BA(filename)
% BA: I created this for single output as structure

MIMICS_10_KEY = '#coordinates written in 1mm / 10000';
MIMICS_13_KEY = '#coordinates written in 1mm / 0';

fid = fopen(filename, 'r');
if (fid == -1)
    fprintf(1, '\aError: unable to open input IV file %s.\n',filename);
    return,
end;

input = fscanf(fid,'%c');
fclose(fid);

% so now lets find the middle points section
if ~isempty(regexpi(input,'(.*)texCoord(.*)')) % Bardiya Added this, in case the model has texture in it
    tokens = regexpi(input,'(.*)point\s+(.*)point\s+\[([^\]]*)\](.*)coordIndex\s+\[([\d\s\,\n\r\.-]*)\](.*)','tokens');
    %now lets get the points
    last_ind = strfind(tokens{1}{2},']');
    pts = sscanf(tokens{1}{2}(2:last_ind-1),'%f %f %f,',[3 inf])';
    connects = sscanf(tokens{1}{5},'%d, %d, %d, %d,',[4 inf])';
else
    tokens = regexpi(input,'(.*)point\s+\[([^\]]*)\](.*)coordIndex\s+\[([\d\s\,\n\r\.-]*)\](.*)','tokens');
    %now lets get the points
    pts = sscanf(tokens{1}{2},'%f %f %f,',[3 inf])';
    connects = sscanf(tokens{1}{4},'%d, %d, %d, %d,',[4 inf])';
end

% tokens = regexpi(input,'(.*)texCoord(.*)');

%check for the way mimics 10, outputs pts
if (isempty(pts)),
    tok = regexp(tokens{1}{2},'^[\s]*#[^\r\n]*(.*)','tokens');
    pts = sscanf(tok{1}{1},'%f %f %f,',[3 inf])';
end;

%check for alternate format, with no commas between the numbers, just at
%the end of the row
if (size(connects,1)==1 && size(connects,2)==1),
    connects = sscanf(tokens{1}{4},'%d %d %d %d',[4 inf])';
end;

%% Look for scale
scale = [1 1 1];

%check in headers
for i=[1 3 5],
    scaleToken = regexpi(tokens{1}{i} ,'scale\s+([\d\.]+)\s+([\d\.]+)\s+([\d\.]+)','tokens');
    if (size(scaleToken,1)>0),
        scale = str2double(scaleToken{1,1});
    end;
end;

%if our scale is not exactly [1 1 1]
if (sum(scale~=1)>0),
    for i=1:3,
        pts(:,i) = pts(:,i)*scale(i);
    end;
end;

% check for the existance of the mimics 10 or 13 key, if so, the output was
% in m, and must be converted to mm
if (~isempty(strfind(input,MIMICS_10_KEY)) || ~isempty(strfind(input,MIMICS_13_KEY)))
    pts = pts .* 1000;
end;

%% BA Output
if abs(mean(mean(abs(pts)))) < 0.5
    disp('Converting Units: m to mm');
    pts = pts.*1e3;
end

model_wrl.pts = pts;
model_wrl.conn = connects;

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%
% Bayesian Optimization of Combinatorial Structures
%
% Copyright (C) 2018 R. Baptista & M. Poloczek
%
% BOCS 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.
%
% BOCS 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 BOCS.  If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) 2018 MIT & University of Arizona
% Authors: Ricardo Baptista & Matthias Poloczek
% E-mails: rsb@mit.edu & poloczek@email.arizona.edu
%

clear; close all; clc
addpath(genpath('../algorithms'))
addpath(genpath('../stat_model'))
addpath(genpath('../test_problems/IsingModel'))
addpath(genpath('../tools'))

% Save inputs in struct
inputs = struct;
inputs.n_vars     = 12;
inputs.evalBudget = 100;
inputs.n_init     = 20;
inputs.lambda     = 1e-4;
inputs.estimator  = 'horseshoe';

% Generate random 3x3 graphical model
Theta   = rand_ising_grid(9);
Moments = ising_model_moments(Theta);

% Save objective function and regularization term
inputs.model    = @(x) KL_divergence_ising(Theta, Moments, x);
inputs.penalty  = @(x) inputs.lambda*sum(x,2);

% Generate initial samples for statistical models
inputs.x_vals   = sample_models(inputs.n_init, inputs.n_vars);
inputs.y_vals   = inputs.model(inputs.x_vals);

% Run BOCS-SDP and BOCS-SA (order 2)
B_SA  = BOCS(inputs.model, inputs.penalty, inputs, 2, 'SA');
B_SDP = BOCS(inputs.model, inputs.penalty, inputs, 2, 'sdp');

% Plot results
figure
hold on
plot(1:length(B_SA.objVals),  cummin(B_SA.objVals), '-r');
plot(1:length(B_SDP.objVals), cummin(B_SDP.objVals), '-.b');
set(gca,'YScale','log')
xlabel('$t$','interpreter','latex')
ylabel('Best $f(x)$','interpreter','latex')
legend({'BOCS - SA','BOCS - SDP'},'location','northeast')
hold off

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
function isMonge=isMongeMatrix(A)
%%ISMONGEMATRIX Determine whether a matrix is a Monge matrix. An n1Xn2
%          matrix is a Monge matrix if A(i1,j1)+A(i2,j2)<=A(i1,j2)+A(i2,j1)
%          for all 1<=i1<i2<=n1 and 1<=j1<j2<=n2. Monge matrices have
%          special properties in some discrete optimization problems.
%
%INPUTS: A An n1Xn2 matrix.
%
%OUTPUTS: isMonge This is true if A is a Monge matrix and false otherwise.
%
%EXAMPLE:
% A=[35, 25, 42, 88;
%    84, 63, 75, 96;
%    51,  9,  4, 25;
%    93, 16,  6, 18];
% isMonge=isMongeMatrix(A)
%A is a Monge matrix so isMonge will be 1.
%
%February 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n1=size(A,1);
n2=size(A,2);

isMonge=true;
for i2=2:n1
    for i1=1:(i2-1)
        for j2=2:n2
            for j1=1:(j2-1)
                if(A(i1,j1)+A(i2,j2)>A(i1,j2)+A(i2,j1))
                   isMonge=false;
                   return;
                end
            end
        end
    end
end

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


function [ report ] = codeReport( varargin )
%CODEREPORT Builds a report of the code health
%
% This function generates a Markdown report on the code health. At the moment
% this is limited to the McCabe (cyclomatic) complexity of a function and its
% subfunctions.
%
% This makes use of |checkcode| in MATLAB.
%
% Usage:
%
%   CODEREPORT('function', functionName) to determine which function is
%   analyzed. (default: matlab2tikz)
%
%   CODEREPORT('complexityThreshold', integer ) to set above which complexity, a
%   function is added to the report (default: 10)
%
%   CODEREPORT('stream', stream) to set to which stream/file to output the report
%   (default: 1, i.e. stdout). The stream is used only when no output argument
%   for `codeReport` is specified!.
%
% See also: checkcode, mlint

    SM = StreamMaker();
    %% input options
    ipp = m2tInputParser();
    ipp = ipp.addParamValue(ipp, 'function', 'matlab2tikz', @ischar);
    ipp = ipp.addParamValue(ipp, 'complexityThreshold', 10, @isnumeric);
    ipp = ipp.addParamValue(ipp, 'stream', 1, SM.isStream);
    ipp = ipp.parse(ipp, varargin{:});

    stream = SM.make(ipp.Results.stream, 'w');

    %% generate report data
    data = checkcode(ipp.Results.function,'-cyc','-struct');
    [complexityAll, mlintMessages] = splitCycloComplexity(data);

    %% analyze cyclomatic complexity
    categorizeComplexity = @(x) categoryOfComplexity(x, ...
                                 ipp.Results.complexityThreshold, ...
                                 ipp.Results.function);

    complexityAll = arrayfun(@parseCycloComplexity, complexityAll);
    complexityAll = arrayfun(categorizeComplexity, complexityAll);

    complexity = filter(complexityAll, @(x) strcmpi(x.category, 'Bad'));
    complexity = sortBy(complexity, 'line', 'ascend');
    complexity = sortBy(complexity, 'complexity', 'descend');

    [complexityStats] = complexityStatistics(complexityAll);

    %% analyze other messages
    %TODO: handle all mlint messages and/or other metrics of the code

    %% format report
    dataStr = complexity;
    dataStr = arrayfun(@(d) mapField(d, 'function',  @markdownInlineCode), dataStr);
    if ~isempty(dataStr)
        dataStr = addFooterRow(dataStr, 'complexity', @sum, {'line',0, 'function',bold('Total')});
    end
    dataStr = arrayfun(@(d) mapField(d, 'line',         @integerToString), dataStr);
    dataStr = arrayfun(@(d) mapField(d, 'complexity',   @integerToString), dataStr);

    report = makeTable(dataStr, {'function', 'complexity'}, ...
                                {'Function', 'Complexity'});

    %% command line usage
    if nargout == 0
        if ismember(stream.name, {'stdout','stderr'})
            stream.print('%s\n', codelinks(report, ipp.Results.function));
        else
            stream.print('%s\n', report);
        end

        figure('name',sprintf('Complexity statistics of %s', ipp.Results.function));
        h = statisticsPlot(complexityStats, 'Complexity', 'Number of functions');
        for hh = h
            plot(hh, [1 1]*ipp.Results.complexityThreshold, ylim(hh), ...
                 'k--','DisplayName','Threshold');
        end
        legend(h(1),'show','Location','NorthEast');

        clear report
    end

end
%% CATEGORIZATION ==============================================================
function [complexity, others] = splitCycloComplexity(list)
% splits codereport into McCabe complexity and others
    filter = @(l) ~isempty(strfind(l.message, 'McCabe complexity'));
    idxComplexity = arrayfun(filter, list);
    complexity = list( idxComplexity);
    others     = list(~idxComplexity);
end
function [data] = categoryOfComplexity(data, threshold, mainFunc)
% categorizes the complexity as "Good", "Bad" or "Accepted"
  TOKEN = '#COMPLEX'; % token to signal allowed complexity

  try %#ok
    helpStr = help(sprintf('%s>%s', mainFunc, data.function));
    if ~isempty(strfind(helpStr, TOKEN))
        data.category = 'Accepted';
        return;
    end
  end
  if data.complexity > threshold
      data.category = 'Bad';
  else
      data.category = 'Good';
  end
end

%% PARSING =====================================================================
function [out] = parseCycloComplexity(in)
% converts McCabe complexity report strings into a better format
    out = regexp(in.message, ...
                 'The McCabe complexity of ''(?<function>[A-Za-z0-9_]+)'' is (?<complexity>[0-9]+).', ...
                 'names');
    out.complexity = str2double(out.complexity);
    out.line = in.line;
end

%% DATA PROCESSING =============================================================
function selected = filter(list, filterFunc)
% filters an array according to a binary function
    idx = logical(arrayfun(filterFunc, list));
    selected = list(idx);
end
function [data] = mapField(data, field, mapping)
    data.(field) = mapping(data.(field));
end
function sorted = sortBy(list, fieldName, mode)
% sorts a struct array by a single field
% extra arguments are as for |sort|
    values = arrayfun(@(m)m.(fieldName), list);
    [dummy, idxSorted] = sort(values(:), 1, mode); %#ok
    sorted = list(idxSorted);
end

function [stat] = complexityStatistics(list)
% calculate some basic statistics of the complexities

    stat.values     = arrayfun(@(c)(c.complexity), list);
    stat.binCenter  = sort(unique(stat.values));

    categoryPerElem = {list.category};
    stat.categories = unique(categoryPerElem);
    nCategories = numel(stat.categories);

    groupedHist = zeros(numel(stat.binCenter), nCategories);
    for iCat = 1:nCategories
        category = stat.categories{iCat};
        idxCat = ismember(categoryPerElem, category);
        groupedHist(:,iCat) = hist(stat.values(idxCat), stat.binCenter);
    end

    stat.histogram  = groupedHist;
    stat.median     = median(stat.values);
end
function [data] = addFooterRow(data, column, func, otherFields)
% adds a footer row to data table based on calculations of a single column
footer = data(end);
for iField = 1:2:numel(otherFields)
    field = otherFields{iField};
    value = otherFields{iField+1};
    footer.(field) = value;
end
footer.(column) = func([data(:).(column)]);
data(end+1) = footer;
end

%% FORMATTING ==================================================================
function str = integerToString(value)
% convert integer to string
    str = sprintf('%d',value);
end
function str = markdownInlineCode(str)
% format as inline code for markdown
    str = sprintf('`%s`', str);
end
function str = makeTable(data, fields, header)
% make a markdown table from struct array
    nData = numel(data);
    str = '';
    if nData == 0
        return; % empty input
    end
    %TODO: use gfmTable from makeTravisReport instead to do the formatting

    % determine column sizes
    nFields = numel(fields);
    table = cell(nFields, nData);
    columnWidth = zeros(1,nFields);
    for iField = 1:nFields
        field = fields{iField};
        table(iField, :) = {data(:).(field)};
        columnWidth(iField) = max(cellfun(@numel, table(iField, :)));
    end
    columnWidth = max(columnWidth, cellfun(@numel, header));
    columnWidth = columnWidth + 2; % empty space left and right
    columnWidth([1,end]) = columnWidth([1,end]) - 1; % except at the edges

    % format table inside cell array
    table = [header; table'];
    for iField = 1:nFields
        FORMAT = ['%' int2str(columnWidth(iField)) 's'];

        for jData = 1:size(table, 1)
            table{jData, iField} = strjust(sprintf(FORMAT, ...
                                           table{jData, iField}), 'center');
        end
    end

    % insert separator
    table = [table(1,:)
             arrayfun(@(n) repmat('-',1,n), columnWidth, 'UniformOutput',false)
             table(2:end,:)]';

    % convert cell array to string
    FORMAT = ['%s' repmat('|%s', 1,nFields-1) '\n'];
    str = sprintf(FORMAT, table{:});

end

function str = codelinks(str, functionName)
% replaces inline functions with clickable links in MATLAB
str = regexprep(str, '`([A-Za-z0-9_]+)`', ...
                ['`<a href="matlab:edit ' functionName '>$1">$1</a>`']);
%NOTE: editing function>subfunction will focus on that particular subfunction
% in the editor (this also works for the main function)
end
function str = bold(str)
str = ['**' str '**'];
end

%% PLOTTING ====================================================================
function h = statisticsPlot(stat, xLabel, yLabel)
% plot a histogram and box plot
    nCategories = numel(stat.categories);
    colors = colorscheme;

    h(1) = subplot(5,1,1:4);
    hold all;
    hb = bar(stat.binCenter, stat.histogram, 'stacked');

    for iCat = 1:nCategories
        category = stat.categories{iCat};

        set(hb(iCat), 'DisplayName', category, 'FaceColor', colors.(category), ...
                   'LineStyle','none');
    end

    %xlabel(xLabel);
    ylabel(yLabel);

    h(2) = subplot(5,1,5);
    hold all;

    boxplot(stat.values,'orientation','horizontal',...
                        'boxstyle',   'outline', ...
                        'symbol',     'o', ...
                        'colors',  colors.All);
    xlabel(xLabel);

    xlims = [min(stat.binCenter)-1 max(stat.binCenter)+1];
    c     = 1;
    ylims = (ylim(h(2)) - c)/3 + c;

    set(h,'XTickMode','manual','XTick',stat.binCenter,'XLim',xlims);
    set(h(1),'XTickLabel','');
    set(h(2),'YTickLabel','','YLim',ylims);
    linkaxes(h, 'x');
end
function colors = colorscheme()
% recognizable color scheme for the categories
 colors.All      = [  0 113 188]/255;
 colors.Good     = [118 171  47]/255;
 colors.Bad      = [161  19  46]/255;
 colors.Accepted = [236 176  31]/255;
end


text:
function ExchangeRxnMatrix = createXMatrix(compoundsIn, transport, compartment)
% Creates a matrix full of exchange reactions based
% on the input list (creates an exchange reaction for each of the
% metabolites present in the model)
%
% USAGE:
%
%    ExchangeRxnMatrix = createXMatrix(compoundsIn, transport, compartment)
%
% INPUTS:
%    compoundsIn:         SU matrix
%    transport:           if 1, transport reactions will be defined as well for
%                         every compound (default: 0, which corresponds to only
%                         exchange reactions (from the specified compartment)
%    compartment:         (default = '[c]') --> transport from cytoplasm [c] to
%                         extracellular space [e] (or sink)
%                         any other compartment (except '[e]' creates transport from [c] to
%                         the compartment, from the compartment to the external compartment
%                         and an exchanger in the external compartment.
% OUTPUT:
%   ExchangeRxnMatrix:    model containing all exchange reactions for all
%                         compounds in compoundsIn
%
% .. Author: - IT 11-10-07
%    Modifications TP October 2017

if ~exist('transport','var') || isempty(transport)
    transport = 0;
end
if ~exist('compartment','var') || isempty(compartment)
    compartment = '[c]';
else
    if strcmp(compartment,'[e]')
        error('Cannot use the external compartment as target compartment');
    end
end

%Bring the compounds in the right dimension.
compoundsIn = columnVector(compoundsIn);

showprogress(0,'Exchange reaction list ...');

ExchangeRxnMatrix = createModel;

cnt=1;
HTABLE = java.util.Hashtable;

%removes the compartment from compound name
compoundsIn=regexprep(compoundsIn,'(\w*)]','');
compoundsIn=regexprep(compoundsIn,'[','');
compoundsIn=regexprep(compoundsIn,']','');
compoundsIn=regexprep(compoundsIn,' ','');

compounds=unique(compoundsIn);

%creates sparse matrix with corresponding dimensions
if  transport==0
    ExchangeRxnMatrix.S=spalloc(length(compounds),length(compounds),length(compounds));
    %Create one sink reaction per input metabolit = 1 S entry, one model
    %metabolite and one reaction per input metabolite
elseif transport == 1
    if (strcmp(compartment,'[c]')==1)
        ExchangeRxnMatrix.S=spalloc(2 * length(compounds),2*length(compounds),3*length(compounds));
        %This has 1 transporter and 1 Exchanger = 3 non zero entries per input metabolite.
        %2 model metabolites per input metabolites and 2 reactions per
        %input metabolite.
    else
        ExchangeRxnMatrix.S=spalloc(3 * length(compounds),3 * length(compounds),5*length(compounds));
        %Generate 2 transporters and 1 Exchanger = 5 non zero entries per input metabolite.
        %3 model metabolites (c,e,comp) per input metabolites and 3 reactions per
        %input metabolite.
    end
end

ExchangeRxnMatrix.mets=compounds;
for i=1:length(compounds)
    HTABLE.put(compounds{i}, i);
end

for i=1:length(compounds)
    if ~isempty(compounds(i))
        if transport == 0
            tmp = ['sink_' compounds(i) compartment];
            ExchangeRxnMatrix.rxns(cnt,1) = strcat(tmp(1),tmp(2),tmp(3));
            tmp = ['sink for ' compounds(i)];
            ExchangeRxnMatrix.rxnsNames(cnt,1) =   strcat(tmp(1),tmp(2));
            tmp = ['1 ' compounds(i) compartment ' <==>'];
            ExchangeRxnMatrix.rxnFormulas(cnt,1) =  strcat(tmp(1),tmp(2),tmp(3));
            tmp = [compounds(i) compartment];
            ExchangeRxnMatrix.mets(i) = strcat(tmp(1),tmp(2));
            %   ExchangeRxnMatrix.grRules{cnt}='';
            [ExchangeRxnMatrix] = addReactionGEM(ExchangeRxnMatrix,ExchangeRxnMatrix.rxns(cnt,1),ExchangeRxnMatrix.rxnsNames(cnt,1),ExchangeRxnMatrix.rxnFormulas(cnt,1),1,-10000,10000);
            cnt = cnt + 1;

        elseif transport == 1 %currently only this branch is taken.
            tmp = ['Ex_' compounds(i) '[e]'];
            ExchangeRxnMatrix.rxns(cnt,1) = strcat(tmp(1),tmp(2),tmp(3));
            tmp = ['Exchange of ' compounds(i)];
            ExchangeRxnMatrix.rxnsNames(cnt,1) =   strcat(tmp(1),tmp(2));
            tmp = ['1 ' compounds(i) '[e] <==>'];
            ExchangeRxnMatrix.rxnFormulas(cnt,1) =  strcat(tmp(1),tmp(2),tmp(3));
            tmp = [compounds(i) '[e]'];
            ExchangeRxnMatrix.mets(i) = strcat(tmp(1),tmp(2));
            HTABLE.put(ExchangeRxnMatrix.mets{i}, i);
            %   ExchangeRxnMatrix.grRules{cnt}='';
            [ExchangeRxnMatrix, HTABLE] = addReactionGEM(ExchangeRxnMatrix,ExchangeRxnMatrix.rxns(cnt,1),ExchangeRxnMatrix.rxnsNames(cnt,1),ExchangeRxnMatrix.rxnFormulas(cnt,1),1,-10000,10000,[],[],[],[],[], HTABLE);
            cnt = cnt + 1;

            if (strcmp(compartment,'[c]')==1)
                % creates transport reaction from [c] to [e]
                tmp = [compounds{i} 'tr'];
                ExchangeRxnMatrix.rxns{cnt,1} = tmp;
                tmp = ['Transport of ' compounds{i}];
                ExchangeRxnMatrix.rxnsNames{cnt,1} =   tmp;
                tmp = ['1 ' compounds{i} '[e] <==> 1 ' compounds{i} '[c]'];
                ExchangeRxnMatrix.rxnFormulas{cnt,1} =  tmp;
                tmp = [compounds{i} '[c]'];
                ExchangeRxnMatrix.mets{length(ExchangeRxnMatrix.mets)+1,1} = tmp;
                HTABLE.put(ExchangeRxnMatrix.mets{end}, length(ExchangeRxnMatrix.mets));
                %  ExchangeRxnMatrix.grRules{cnt}='';
                [ExchangeRxnMatrix, HTABLE] = addReactionGEM(ExchangeRxnMatrix,ExchangeRxnMatrix.rxns(cnt,1),ExchangeRxnMatrix.rxnsNames(cnt,1),ExchangeRxnMatrix.rxnFormulas(cnt,1),1,-10000,10000,[],[],[],[],[], HTABLE);
                cnt = cnt + 1;

            else % keep this branch the same for now.
                % creates transport reaction from [c] to [p]
                tmp = [compounds{i} 't' compartment(2:end-1) 'r'];
                ExchangeRxnMatrix.rxns{cnt,1} = tmp;
                tmp = ['[c] to ' compartment ' Transport of ' compounds{i}];
                ExchangeRxnMatrix.rxnsNames{cnt,1} =   tmp;
                tmp = ['1 ' compounds{i} compartment ' <==> 1 ' compounds{i} '[c]'];
                ExchangeRxnMatrix.rxnFormulas{cnt,1} =  tmp;
                tmp = [compounds{i} '[c]'];
                ExchangeRxnMatrix.mets{length(ExchangeRxnMatrix.mets)+1,1} = tmp;
                tmp = [compounds{i} compartment];
                ExchangeRxnMatrix.mets{length(ExchangeRxnMatrix.mets)+1,1} = tmp;
                %  ExchangeRxnMatrix.grRules{cnt}='';
                [ExchangeRxnMatrix] = addReactionGEM(ExchangeRxnMatrix,ExchangeRxnMatrix.rxns(cnt,1),ExchangeRxnMatrix.rxnsNames(cnt,1),ExchangeRxnMatrix.rxnFormulas(cnt,1),1,-10000,10000);
                cnt = cnt + 1;

                % creates transport reaction from [p] to [e]
                tmp = [compounds{i} 'tr'];
                ExchangeRxnMatrix.rxns{cnt,1} = tmp;
                tmp = [compartment ' to [e] Transport of ' compounds{i}];
                ExchangeRxnMatrix.rxnsNames{cnt,1} =   tmp;
                tmp = ['1 ' compounds{i} '[e] <==> 1 ' compounds{i} compartment];
                ExchangeRxnMatrix.rxnFormulas{cnt,1} =  tmp;
                %    ExchangeRxnMatrix.grRules{cnt}='';
                [ExchangeRxnMatrix] = addReactionGEM(ExchangeRxnMatrix,ExchangeRxnMatrix.rxns(cnt,1),ExchangeRxnMatrix.rxnsNames(cnt,1),ExchangeRxnMatrix.rxnFormulas(cnt,1),1,-10000,10000);
                cnt = cnt + 1;
            end
        end
    end

    showprogress(i/length(compounds));
end

author: opencobra

repo: cobratoolbox


text:
function outstr = GE_dumpImageHeader(img)
%
% outstr = GE_dumpImageHeader(img)
%
% Writes the image header to the string outstr
%
% Souheil J. Inati
% Dartmouth College
% May 2000
% souheil.inati@dartmouth.edu
%

outstr = sprintf('im_suid: %s\n',         char(img.im_suid));
outstr  = strcat(outstr, sprintf('im_uniq: %d\n',         img.im_uniq));
outstr  = strcat(outstr, sprintf('im_diskid: %s\n',       char(img.im_diskid)));
outstr  = strcat(outstr, sprintf('im_exno: %d\n',         img.im_exno));
outstr  = strcat(outstr, sprintf('im_seno: %d\n',         img.im_seno));
outstr  = strcat(outstr, sprintf('im_no: %d\n',           img.im_no));
outstr  = strcat(outstr, sprintf('im_datetime: %d\n',     img.im_datetime));
outstr  = strcat(outstr, sprintf('im_actual_dt: %d\n',    img.im_actual_dt));
outstr  = strcat(outstr, sprintf('sctime: %f\n',          img.sctime));
outstr  = strcat(outstr, sprintf('slthick: %f\n',         img.slthick));
outstr  = strcat(outstr, sprintf('imatrix_X: %d\n',       img.imatrix_X));
outstr  = strcat(outstr, sprintf('imatrix_Y: %d\n',       img.imatrix_Y));
outstr  = strcat(outstr, sprintf('dfov: %f\n',            img.dfov));
outstr  = strcat(outstr, sprintf('dfov_rect: %f\n',       img.dfov_rect));
outstr  = strcat(outstr, sprintf('dim_X: %f\n',           img.dim_X));
outstr  = strcat(outstr, sprintf('dim_Y: %f\n',           img.dim_Y));
outstr  = strcat(outstr, sprintf('pixsize_X: %f\n',       img.pixsize_X));
outstr  = strcat(outstr, sprintf('pixsize_Y: %f\n',       img.pixsize_Y));
outstr  = strcat(outstr, sprintf('pdid: %s\n',            sprintf('%s',img.pdid)));
outstr  = strcat(outstr, sprintf('contrastIV: %s\n',      cstring(img.contrastIV)));
outstr  = strcat(outstr, sprintf('contrastOral: %s\n',    cstring(img.contrastOral)));
outstr  = strcat(outstr, sprintf('contmode: %d\n',        img.contmode));
outstr  = strcat(outstr, sprintf('serrx: %d\n',           img.serrx));
outstr  = strcat(outstr, sprintf('imgrx: %d\n',           img.imgrx));
outstr  = strcat(outstr, sprintf('screenformat: %d\n',    img.screenformat));
outstr  = strcat(outstr, sprintf('plane: %d\n',           img.plane));
outstr  = strcat(outstr, sprintf('scanspacing: %f\n',     img.scanspacing));
outstr  = strcat(outstr, sprintf('im_compress: %d\n',     img.im_compress));
outstr  = strcat(outstr, sprintf('im_scouttype: %d\n',    img.im_scouttype));
outstr  = strcat(outstr, sprintf('loc_ras: %s\n',         char(img.loc_ras)));
outstr  = strcat(outstr, sprintf('loc: %f\n',          img.loc));
outstr  = strcat(outstr, sprintf('ctr_R: %f\n',        img.ctr_R));
outstr  = strcat(outstr, sprintf('ctr_A: %f\n',        img.ctr_A));
outstr  = strcat(outstr, sprintf('ctr_S: %f\n',        img.ctr_S));
outstr  = strcat(outstr, sprintf('norm_R: %f\n',       img.norm_R));
outstr  = strcat(outstr, sprintf('norm_A: %f\n',       img.norm_A));
outstr  = strcat(outstr, sprintf('norm_S: %f\n',       img.norm_S));
outstr  = strcat(outstr, sprintf('tlhc_R: %f\n',       img.tlhc_R));
outstr  = strcat(outstr, sprintf('tlhc_A: %f\n',       img.tlhc_A));
outstr  = strcat(outstr, sprintf('tlhc_S: %f\n',       img.tlhc_S));
outstr  = strcat(outstr, sprintf('trhc_R: %f\n',       img.trhc_R));
outstr  = strcat(outstr, sprintf('trhc_A: %f\n',       img.trhc_A));
outstr  = strcat(outstr, sprintf('trhc_S: %f\n',       img.trhc_S));
outstr  = strcat(outstr, sprintf('brhc_R: %f\n',       img.brhc_R));
outstr  = strcat(outstr, sprintf('brhc_A: %f\n',       img.brhc_A));
outstr  = strcat(outstr, sprintf('brhc_S: %f\n',       img.brhc_S));
outstr  = strcat(outstr, sprintf('forimgrev: %s\n',    char(img.forimgrev)));
outstr  = strcat(outstr, sprintf('tr: %d\n',           img.tr));
outstr  = strcat(outstr, sprintf('ti: %d\n',           img.ti));
outstr  = strcat(outstr, sprintf('te: %d\n',           img.te));
outstr  = strcat(outstr, sprintf('te2: %d\n',          img.te2));
outstr  = strcat(outstr, sprintf('numecho: %d\n',      img.numecho));
outstr  = strcat(outstr, sprintf('echonum: %d\n',      img.echonum));
outstr  = strcat(outstr, sprintf('tbldlta: %f\n',      img.tbldlta));
outstr  = strcat(outstr, sprintf('nex: %f\n',          img.nex));
outstr  = strcat(outstr, sprintf('contig: %d\n',       img.contig));
outstr  = strcat(outstr, sprintf('hrtrate: %d\n',      img.hrtrate));
outstr  = strcat(outstr, sprintf('tdel: %d\n',         img.tdel));
outstr  = strcat(outstr, sprintf('saravg: %f\n',       img.saravg));
outstr  = strcat(outstr, sprintf('sarpeak: %f\n',      img.sarpeak));
outstr  = strcat(outstr, sprintf('monsar: %d\n',          img.monsar));
outstr  = strcat(outstr, sprintf('trgwindow: %d\n',       img.trgwindow));
outstr  = strcat(outstr, sprintf('reptime: %f\n',         img.reptime));
outstr  = strcat(outstr, sprintf('imgpcyc: %d\n',         img.imgpcyc));
outstr  = strcat(outstr, sprintf('xmtgain: %d\n',         img.xmtgain));
outstr  = strcat(outstr, sprintf('rcvgain1: %d\n',        img.rcvgain1));
outstr  = strcat(outstr, sprintf('rcvgain2: %d\n',        img.rcvgain2));
outstr  = strcat(outstr, sprintf('mr_flip: %d\n',         img.mr_flip));
outstr  = strcat(outstr, sprintf('mindat: %d\n',          img.mindat));
outstr  = strcat(outstr, sprintf('cphase: %d\n',          img.cphase));
outstr  = strcat(outstr, sprintf('swappf: %d\n',          img.swappf));
outstr  = strcat(outstr, sprintf('pauseint: %d\n',        img.pauseint));
outstr  = strcat(outstr, sprintf('pausetime: %f\n',       img.pausetime));
outstr  = strcat(outstr, sprintf('obplane: %d\n',         img.obplane));
outstr  = strcat(outstr, sprintf('slocfov: %d\n',         img.slocfov));
outstr  = strcat(outstr, sprintf('xmtfreq: %d\n',         img.xmtfreq));
outstr  = strcat(outstr, sprintf('autoxmtfreq: %d\n',     img.autoxmtfreq));
outstr  = strcat(outstr, sprintf('autoxmtgain: %d\n',     img.autoxmtgain));
outstr  = strcat(outstr, sprintf('prescan_r1: %d\n',      img.prescan_r1));
outstr  = strcat(outstr, sprintf('prescan_r2: %d\n',      img.prescan_r2));
outstr  = strcat(outstr, sprintf('user_bitmap: %d\n',     img.user_bitmap));
outstr  = strcat(outstr, sprintf('cenfreq: %d\n',         img.cenfreq));
outstr  = strcat(outstr, sprintf('imode: %d\n',           img.imode));
outstr  = strcat(outstr, sprintf('iopt: %d\n',            img.iopt));
outstr  = strcat(outstr, sprintf('pseq: %d\n',            img.pseq));
outstr  = strcat(outstr, sprintf('pseqmode: %d\n',        img.pseqmode));
outstr  = strcat(outstr, sprintf('psdname: %s\n',         cstring(img.psdname)));
outstr  = strcat(outstr, sprintf('psd_datetime: %d\n',    img.psd_datetime));
outstr  = strcat(outstr, sprintf('psd_iname: %s\n',       cstring(img.psd_iname)));
outstr  = strcat(outstr, sprintf('ctyp: %d\n',            img.ctyp));
outstr  = strcat(outstr, sprintf('cname: %s\n',         cstring(img.cname)));
outstr  = strcat(outstr, sprintf('surfctyp: %d\n',      img.surfctyp));
outstr  = strcat(outstr, sprintf('surfcext: %d\n',      img.surfcext));
outstr  = strcat(outstr, sprintf('rawrunnum: %d\n',     img.rawrunnum));
outstr  = strcat(outstr, sprintf('cal_fldstr: %d\n',    img.cal_fldstr));
outstr  = strcat(outstr, sprintf('supp_tech: %d\n',     img.supp_tech));
outstr  = strcat(outstr, sprintf('vbw: %f\n',           img.vbw));
outstr  = strcat(outstr, sprintf('slquant: %d\n',       img.slquant));
outstr  = strcat(outstr, sprintf('gpre: %d\n',          img.gpre));
outstr  = strcat(outstr, sprintf('intr_del: %d\n',      img.intr_del));
outstr  = strcat(outstr, sprintf('user0: %f\n',         img.user0));
outstr  = strcat(outstr, sprintf('user1: %f\n',         img.user1));
outstr  = strcat(outstr, sprintf('user2: %f\n',         img.user2));
outstr  = strcat(outstr, sprintf('user3: %f\n',         img.user3));
outstr  = strcat(outstr, sprintf('user4: %f\n',         img.user4));
outstr  = strcat(outstr, sprintf('user5: %f\n',         img.user5));
outstr  = strcat(outstr, sprintf('user6: %f\n',         img.user6));
outstr  = strcat(outstr, sprintf('user7: %f\n',         img.user7));
outstr  = strcat(outstr, sprintf('user8: %f\n',         img.user8));
outstr  = strcat(outstr, sprintf('user9: %f\n',         img.user9));
outstr  = strcat(outstr, sprintf('user10: %f\n',        img.user10));
outstr  = strcat(outstr, sprintf('user11: %f\n',        img.user11));
outstr  = strcat(outstr, sprintf('user12: %f\n',        img.user12));
outstr  = strcat(outstr, sprintf('user13: %f\n',        img.user13));
outstr  = strcat(outstr, sprintf('user14: %f\n',        img.user14));
outstr  = strcat(outstr, sprintf('user15: %f\n',        img.user15));
outstr  = strcat(outstr, sprintf('user16: %f\n',        img.user16));
outstr  = strcat(outstr, sprintf('user17: %f\n',        img.user17));
outstr  = strcat(outstr, sprintf('user18: %f\n',        img.user18));
outstr  = strcat(outstr, sprintf('user19: %f\n',        img.user19));
outstr  = strcat(outstr, sprintf('user20: %f\n',          img.user20));
outstr  = strcat(outstr, sprintf('user21: %f\n',          img.user21));
outstr  = strcat(outstr, sprintf('user22: %f\n',          img.user22));
outstr  = strcat(outstr, sprintf('user23: %f\n',          img.user23));
outstr  = strcat(outstr, sprintf('user24: %f\n',          img.user24));
outstr  = strcat(outstr, sprintf('im_alloc_key: %s\n',    sprintf('%s',img.im_alloc_key)));
outstr  = strcat(outstr, sprintf('im_lastmod: %d\n',      img.im_lastmod));
outstr  = strcat(outstr, sprintf('im_verscre: %s\n',      char(img.im_verscre)));
outstr  = strcat(outstr, sprintf('im_verscur: %s\n',      char(img.im_verscur)));
outstr  = strcat(outstr, sprintf('im_pds_a: %d\n',        img.im_pds_a));
outstr  = strcat(outstr, sprintf('im_pds_c: %d\n',        img.im_pds_c));
outstr  = strcat(outstr, sprintf('im_pds_u: %d\n',        img.im_pds_u));
outstr  = strcat(outstr, sprintf('im_checksum: %d\n',     img.im_checksum));
outstr  = strcat(outstr, sprintf('im_archived: %d\n',     img.im_archived));
outstr  = strcat(outstr, sprintf('im_complete: %d\n',     img.im_complete));
outstr  = strcat(outstr, sprintf('satbits: %d\n',         img.satbits));
outstr  = strcat(outstr, sprintf('scic: %d\n',            img.scic));
outstr  = strcat(outstr, sprintf('satxloc1: %d\n',        img.satxloc1));
outstr  = strcat(outstr, sprintf('satxloc2: %d\n',        img.satxloc2));
outstr  = strcat(outstr, sprintf('satyloc1: %d\n',        img.satyloc1));
outstr  = strcat(outstr, sprintf('satyloc2: %d\n',        img.satyloc2));
outstr  = strcat(outstr, sprintf('satzloc1: %d\n',        img.satzloc1));
outstr  = strcat(outstr, sprintf('satzloc2: %d\n',        img.satzloc2));
outstr  = strcat(outstr, sprintf('satxthick: %d\n',       img.satxthick));
outstr  = strcat(outstr, sprintf('satythick: %d\n',       img.satythick));
outstr  = strcat(outstr, sprintf('satzthick: %d\n',       img.satzthick));
outstr  = strcat(outstr, sprintf('flax: %d\n',            img.flax));
outstr  = strcat(outstr, sprintf('venc: %d\n',            img.venc));
outstr  = strcat(outstr, sprintf('thk_disclmr: %d\n',     img.thk_disclmr));
outstr  = strcat(outstr, sprintf('ps_flag: %d\n',         img.ps_flag));
outstr  = strcat(outstr, sprintf('ps_status: %d\n',       img.ps_status));
outstr  = strcat(outstr, sprintf('image_type: %d\n',      img.image_type));
outstr  = strcat(outstr, sprintf('vas_collapse: %d\n',    img.vas_collapse));
outstr  = strcat(outstr, sprintf('user23n: %f\n',         img.user23n));
outstr  = strcat(outstr, sprintf('user24n: %f\n',         img.user24n));
outstr  = strcat(outstr, sprintf('proj_alg: %d\n',        img.proj_alg));
outstr  = strcat(outstr, sprintf('proj_name: %s\n',       cstring(img.proj_name)));
outstr  = strcat(outstr, sprintf('x_axis_rot: %f\n',      img.x_axis_rot));
outstr  = strcat(outstr, sprintf('y_axis_rot: %f\n',      img.y_axis_rot));
outstr  = strcat(outstr, sprintf('z_axis_rot: %f\n',      img.z_axis_rot));
outstr  = strcat(outstr, sprintf('thresh_min1: %d\n',     img.thresh_min1));
outstr  = strcat(outstr, sprintf('thresh_max1: %d\n',     img.thresh_max1));
outstr  = strcat(outstr, sprintf('thresh_min2: %d\n',     img.thresh_min2));
outstr  = strcat(outstr, sprintf('thresh_max2: %d\n',     img.thresh_max2));
outstr  = strcat(outstr, sprintf('echo_trn_len: %d\n',    img.echo_trn_len));
outstr  = strcat(outstr, sprintf('frac_echo: %d\n',       img.frac_echo));
outstr  = strcat(outstr, sprintf('prep_pulse: %d\n',      img.prep_pulse));
outstr  = strcat(outstr, sprintf('cphasenum: %d\n',       img.cphasenum));
outstr  = strcat(outstr, sprintf('var_echo: %d\n',        img.var_echo));
outstr  = strcat(outstr, sprintf('ref_img: %s\n',         char(img.ref_img)));
outstr  = strcat(outstr, sprintf('sum_img: %s\n',         char(img.sum_img)));
outstr  = strcat(outstr, sprintf('img_window: %d\n',      img.img_window));
outstr  = strcat(outstr, sprintf('img_level: %d\n',       img.img_level));
outstr  = strcat(outstr, sprintf('slop_int_1: %d\n',      img.slop_int_1));
outstr  = strcat(outstr, sprintf('slop_int_2: %d\n',      img.slop_int_2));
outstr  = strcat(outstr, sprintf('slop_int_3: %d\n',      img.slop_int_3));
outstr  = strcat(outstr, sprintf('slop_int_4: %d\n',      img.slop_int_4));
outstr  = strcat(outstr, sprintf('slop_int_5: %d\n',      img.slop_int_5));
outstr  = strcat(outstr, sprintf('slop_float_1: %f\n',    img.slop_float_1));
outstr  = strcat(outstr, sprintf('slop_float_2: %f\n',    img.slop_float_2));
outstr  = strcat(outstr, sprintf('slop_float_3: %f\n',    img.slop_float_3));
outstr  = strcat(outstr, sprintf('slop_float_4: %f\n',    img.slop_float_4));
outstr  = strcat(outstr, sprintf('slop_float_5: %f\n',    img.slop_float_5));
outstr  = strcat(outstr, sprintf('slop_str_1: %s\n',      cstring(img.slop_str_1)));
outstr  = strcat(outstr, sprintf('slop_str_2: %s\n',      cstring(img.slop_str_2)));
outstr  = strcat(outstr, sprintf('scanactno: %d\n',       img.scanactno));
outstr  = strcat(outstr, sprintf('vasflags: %d\n',        img.vasflags));
outstr  = strcat(outstr, sprintf('vencscale: %f\n',       img.vencscale));
outstr  = strcat(outstr, sprintf('integrity: %d\n',       img.integrity));
outstr  = strcat(outstr, sprintf('fphase: %d\n',          img.fphase));
outstr  = strcat(outstr, sprintf('freq_dir: %d\n',        img.freq_dir));
outstr  = strcat(outstr, sprintf('vas_mode: %d\n',        img.vas_mode));
outstr  = strcat(outstr, sprintf('image_uid: %s\n',       sprintf('%s',img.image_uid)));
outstr  = strcat(outstr, sprintf('sop_uid: %s\n',         sprintf('%s',img.sop_uid)));
outstr  = strcat(outstr, sprintf('dont_use_1: %d\n',      img.dont_use_1));
outstr  = strcat(outstr, sprintf('dont_use_2: %d\n',      img.dont_use_2));
outstr  = strcat(outstr, sprintf('dont_use_3: %d\n',      img.dont_use_3));
outstr  = strcat(outstr, sprintf('pscopts: %d\n',         img.pscopts));
outstr  = strcat(outstr, sprintf('asoffsetx: %d\n',       img.asoffsetx));
outstr  = strcat(outstr, sprintf('asoffsety: %d\n',       img.asoffsety));
outstr  = strcat(outstr, sprintf('asoffsetz: %d\n',           img.asoffsetz));
outstr  = strcat(outstr, sprintf('unoriginal: %d\n',          img.unoriginal));
outstr  = strcat(outstr, sprintf('interleaves: %d\n',         img.interleaves));
outstr  = strcat(outstr, sprintf('effechospace: %d\n',        img.effechospace));
outstr  = strcat(outstr, sprintf('viewsperseg: %d\n',         img.viewsperseg));
outstr  = strcat(outstr, sprintf('rbpm: %d\n',                img.rbpm));
outstr  = strcat(outstr, sprintf('rtpoint: %d\n',             img.rtpoint));
outstr  = strcat(outstr, sprintf('rcvrtype: %d\n',            img.rcvrtype));
outstr  = strcat(outstr, sprintf('dbdt: %f\n',                img.dbdt));
outstr  = strcat(outstr, sprintf('dbdtper: %f\n',             img.dbdtper));
outstr  = strcat(outstr, sprintf('estdbdtper: %f\n',          img.estdbdtper));
outstr  = strcat(outstr, sprintf('estdbdtts: %f\n',           img.estdbdtts));
outstr  = strcat(outstr, sprintf('saravghead: %f\n',          img.saravghead));
outstr  = strcat(outstr, sprintf('neg_scanspacing: %f\n',     img.neg_scanspacing));
outstr  = strcat(outstr, sprintf('offsetfreq: %d\n',          img.offsetfreq));
outstr  = strcat(outstr, sprintf('user_usage_tag: %d\n',      img.user_usage_tag));
outstr  = strcat(outstr, sprintf('user_fill_mapMSW: %d\n',    img.user_fill_mapMSW));
outstr  = strcat(outstr, sprintf('user_fill_mapLSW: %d\n',    img.user_fill_mapLSW));
outstr  = strcat(outstr, sprintf('user25: %f\n',              img.user25));
outstr  = strcat(outstr, sprintf('user26: %f\n',              img.user26));
outstr  = strcat(outstr, sprintf('user27: %f\n',              img.user27));
outstr  = strcat(outstr, sprintf('user28: %f\n',              img.user28));
outstr  = strcat(outstr, sprintf('user29: %f\n',              img.user29));
outstr  = strcat(outstr, sprintf('user30: %f\n',              img.user30));
outstr  = strcat(outstr, sprintf('user31: %f\n',              img.user31));
outstr  = strcat(outstr, sprintf('user32: %f\n',              img.user32));
outstr  = strcat(outstr, sprintf('user33: %f\n',              img.user33));
outstr  = strcat(outstr, sprintf('user34: %f\n',              img.user34));
outstr  = strcat(outstr, sprintf('user35: %f\n',              img.user35));
outstr  = strcat(outstr, sprintf('user36: %f\n',              img.user36));
outstr  = strcat(outstr, sprintf('user37: %f\n',              img.user37));
outstr  = strcat(outstr, sprintf('user38: %f\n',              img.user38));
outstr  = strcat(outstr, sprintf('user39: %f\n',              img.user39));
outstr  = strcat(outstr, sprintf('user40: %f\n',              img.user40));
outstr  = strcat(outstr, sprintf('user41: %f\n',              img.user41));
outstr  = strcat(outstr, sprintf('user42: %f\n',              img.user42));
outstr  = strcat(outstr, sprintf('user43: %f\n',              img.user43));
outstr  = strcat(outstr, sprintf('user44: %f\n',              img.user44));
outstr  = strcat(outstr, sprintf('user45: %f\n',              img.user45));
outstr  = strcat(outstr, sprintf('user46: %f\n',              img.user46));
outstr  = strcat(outstr, sprintf('user47: %f\n',              img.user47));
outstr  = strcat(outstr, sprintf('user48: %f\n',              img.user48));
outstr  = strcat(outstr, sprintf('slop_int_6: %d\n',          img.slop_int_6));
outstr  = strcat(outstr, sprintf('slop_int_7: %d\n',          img.slop_int_7));
outstr  = strcat(outstr, sprintf('slop_int_8: %d\n',          img.slop_int_8));
outstr  = strcat(outstr, sprintf('slop_int_9: %d\n',          img.slop_int_9));
outstr  = strcat(outstr, sprintf('mr_padding: %s\n',          cstring(img.mr_padding)));

return

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [y,K] = spm_rand_power_law(csd,Hz,dt,N)
% generates random variates with a power law spectral density
% FORMAT [y,K] = spm_rand_power_law(csd,Hz,dt,N)
% csd - spectral densities (one per column)
% Hz  - frequencies
% dt  - sampling interval
% N   - number of time bins
%
% y   - random variate
% K   - convolution (kernel) operator: y(:,i) = K*randn(N,1)
%
% see also: spm_rand_mar; spm_Q
%__________________________________________________________________________
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_rand_power_law.m 6759 2016-03-27 19:45:17Z karl $

% create random process
%==========================================================================

% create AR representation and associated convolution kernels
%--------------------------------------------------------------------------
p     = fix(length(Hz) - 1);
for i = 1:size(csd,2);
    ccf    = spm_csd2ccf(csd(:,i),Hz,dt);
    mar    = spm_ccf2mar(ccf,p);
    A      = [1 -mar.a'];
    P      = spdiags(ones(N,1)*A,-(0:mar.p),N,N);
    y(:,i) = P\randn(N,1)*sqrt(max(ccf));
end

% convolution kernel
%--------------------------------------------------------------------------
if nargout > 1
    K = inv(P)*sqrt(max(ccf));
end

return

% NB: alternative scheme - create random process
%==========================================================================
[m,n] = size(csd);
w     = (0:(N - 1))/dt/N;
dHz   = Hz(2) - Hz(1);
g     = zeros(N,n);
for i = 1:m
    j      = find(w > Hz(i),1);
    s      = sqrt(csd(i,:)).*(randn(1,n) + 1j*randn(1,n));
    g(j,:) = s;
    j      = N - j + 2;
    g(j,:) = conj(s);

end
y     = real(ifft(g));
y     = y*sqrt(mean(sum(G)*dHz)/mean(var(y)));

author: spm


text:
function [slug] = g2slug(g)
% Convert mass from grams to slugs.
% Chad Greene 2012
slug = g*(0.002204622621849)/32.17405;

author: Sable

repo: mcbench-benchmarks


text:
function lambda = idem_random_eigenvalues ( n, rank, key )

%*****************************************************************************80
%
%% IDEM_RANDOM_EIGENVALUES returns the eigenvalues of the IDEM_RANDOM matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, integer RANK, the rank of A.
%
%    Input, integer KEY, a positive value that selects the data.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  lambda(1:rank,1) = 1.0;
  lambda(rank+1:n,1) = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ripples,sd,bad] = FindRipples(filtered,varargin)

%FindRipples - Find hippocampal ripples (100~200Hz oscillations).
%
%  USAGE
%
%    [ripples,stdev,noise] = FindRipples(filtered,<options>)
%
%    Ripples are detected using the normalized squared signal (NSS) by
%    thresholding the baseline, merging neighboring events, thresholding
%    the peaks, and discarding events with excessive duration.
%    Thresholds are computed as multiples of the standard deviation of
%    the NSS. Alternatively, one can use explicit values, typically obtained
%    from a previous call.
%
%    filtered       ripple-band filtered LFP (one channel).
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'thresholds'  thresholds for ripple beginning/end and peak, in multiples
%                   of the stdev (default = [2 5])
%     'durations'   minimum inter-ripple interval, and minimum and maximum
%                   ripple durations, in ms (default = [30 20 100])
%     'baseline'    interval used to compute normalization (default = all)
%     'restrict'    same as 'baseline' (for backwards compatibility)
%     'frequency'   sampling rate (in Hz) (default = 1250Hz)
%     'stdev'       reuse previously computed stdev
%     'show'        plot results (default = 'off')
%     'noise'       noisy ripple-band filtered channel used to exclude ripple-
%                   like noise (events also present on this channel are
%                   discarded)
%    =========================================================================
%
%  OUTPUT
%
%    ripples        for each ripple, [start_t peak_t end_t peakNormalizedPower]
%    stdev          standard deviation of the NSS (can be reused subsequently)
%    noise          ripple-like activity recorded simultaneously on the noise
%                   channel (for debugging info)
%
%  SEE
%
%    See also FilterLFP, RippleStats, SaveRippleEvents, PlotRippleStats.

% Copyright (C) 2004-2011 by Michaël Zugaro, initial algorithm by Hajime Hirase
%
% This program 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.

% Default values
frequency = 1250;
show = 'off';
restrict = [];
sd = [];
lowThresholdFactor = 2; % Ripple envoloppe must exceed lowThresholdFactor*stdev
highThresholdFactor = 5; % Ripple peak must exceed highThresholdFactor*stdev
minInterRippleInterval = 30; % in ms
minRippleDuration = 20; % in ms
maxRippleDuration = 100; % in ms
noise = [];

% Check number of parameters
if nargin < 1 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help FindRipples'' for details).');
end

% Check parameter sizes
if ~isdmatrix(filtered) | size(filtered,2) ~= 2,
	error('Parameter ''filtered'' is not a Nx2 matrix (type ''help FindRipples'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property (type ''help FindRipples'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'thresholds',
			thresholds = varargin{i+1};
			if ~isdvector(thresholds,'#2','>0'),
				error('Incorrect value for property ''thresholds'' (type ''help FindRipples'' for details).');
			end
			lowThresholdFactor = thresholds(1);
			highThresholdFactor = thresholds(2);
		case 'durations',
			durations = varargin{i+1};
			if ~isdvector(durations,'#2','>0') && ~isdvector(durations,'#3','>0'),
				error('Incorrect value for property ''durations'' (type ''help FindRipples'' for details).');
			end
			if length(durations) == 2,
				minInterRippleInterval = durations(1);
				maxRippleDuration = durations(2);
			else
				minInterRippleInterval = durations(1);
				minRippleDuration = durations(2);
				maxRippleDuration = durations(3);
			end
		case 'frequency',
			frequency = varargin{i+1};
			if ~isdscalar(frequency,'>0'),
				error('Incorrect value for property ''frequency'' (type ''help FindRipples'' for details).');
			end
		case 'show',
			show = varargin{i+1};
			if ~isstring(show,'on','off'),
				error('Incorrect value for property ''show'' (type ''help FindRipples'' for details).');
			end
		case {'baseline','restrict'},
			restrict = varargin{i+1};
			if ~isempty(restrict) & ~isdvector(restrict,'#2','<'),
				error('Incorrect value for property ''restrict'' (type ''help FindRipples'' for details).');
			end
		case 'stdev',
			sd = varargin{i+1};
			if ~isdscalar(sd,'>0'),
				error('Incorrect value for property ''stdev'' (type ''help FindRipples'' for details).');
			end
		case 'noise',
			noise = varargin{i+1};
			if ~isdmatrix(noise) | size(noise,1) ~= size(filtered,1) | size(noise,2) ~= 2,
				error('Incorrect value for property ''noise'' (type ''help FindRipples'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help FindRipples'' for details).']);
	end
end

% Parameters
windowLength = round(frequency/1250*11);

% Square and normalize signal
signal = filtered(:,2);
squaredSignal = signal.^2;
window = ones(windowLength,1)/windowLength;
keep = [];
if ~isempty(restrict),
	keep = filtered(:,1)>=restrict(1)&filtered(:,1)<=restrict(2);
end

[normalizedSquaredSignal,sd] = unity(Filter0(window,sum(squaredSignal,2)),sd,keep);

% Detect ripple periods by thresholding normalized squared signal
thresholded = normalizedSquaredSignal > lowThresholdFactor;
start = find(diff(thresholded)>0);
stop = find(diff(thresholded)<0);
% Exclude last ripple if it is incomplete
if length(stop) == length(start)-1,
	start = start(1:end-1);
end
% Exclude first ripple if it is incomplete
if length(stop)-1 == length(start),
    stop = stop(2:end);
end
% Correct special case when both first and last ripples are incomplete
if start(1) > stop(1),
	stop(1) = [];
	start(end) = [];
end
firstPass = [start,stop];
if isempty(firstPass),
	disp('Detection by thresholding failed');
	return
else
	disp(['After detection by thresholding: ' num2str(length(firstPass)) ' events.']);
end

% Merge ripples if inter-ripple period is too short
minInterRippleSamples = minInterRippleInterval/1000*frequency;
secondPass = [];
ripple = firstPass(1,:);
for i = 2:size(firstPass,1)
	if firstPass(i,1) - ripple(2) < minInterRippleSamples,
		% Merge
		ripple = [ripple(1) firstPass(i,2)];
	else
		secondPass = [secondPass ; ripple];
		ripple = firstPass(i,:);
	end
end
secondPass = [secondPass ; ripple];
if isempty(secondPass),
	disp('Ripple merge failed');
	return
else
	disp(['After ripple merge: ' num2str(length(secondPass)) ' events.']);
end

% Discard ripples with a peak power < highThresholdFactor
thirdPass = [];
peakNormalizedPower = [];
for i = 1:size(secondPass,1)
	[maxValue,maxIndex] = max(normalizedSquaredSignal([secondPass(i,1):secondPass(i,2)]));
	if maxValue > highThresholdFactor,
		thirdPass = [thirdPass ; secondPass(i,:)];
		peakNormalizedPower = [peakNormalizedPower ; maxValue];
	end
end
if isempty(thirdPass),
	disp('Peak thresholding failed.');
	return
else
	disp(['After peak thresholding: ' num2str(length(thirdPass)) ' events.']);
end

% Detect negative peak position for each ripple
peakPosition = zeros(size(thirdPass,1),1);
for i=1:size(thirdPass,1),
	[minValue,minIndex] = min(signal(thirdPass(i,1):thirdPass(i,2)));
	peakPosition(i) = minIndex + thirdPass(i,1) - 1;
end

% Discard ripples that are way too short
time = filtered(:,1);
ripples = [time(thirdPass(:,1)) time(peakPosition) time(thirdPass(:,2)) peakNormalizedPower];
duration = ripples(:,3)-ripples(:,1);
ripples(duration<minRippleDuration/1000,:) = [];
disp(['After min duration test: ' num2str(size(ripples,1)) ' events.']);

% Discard ripples that are way too long
duration = ripples(:,3)-ripples(:,1);
ripples(duration>maxRippleDuration/1000,:) = [];
disp(['After max duration test: ' num2str(size(ripples,1)) ' events.']);

% If a noisy channel was provided, find ripple-like events and exclude them
bad = [];
if ~isempty(noise),
	% Square and pseudo-normalize (divide by signal stdev) noise
	squaredNoise = noise(:,2).^2;
	window = ones(windowLength,1)/windowLength;
	normalizedSquaredNoise = unity(Filter0(window,sum(squaredNoise,2)),sd,[]);
	excluded = logical(zeros(size(ripples,1),1));
	% Exclude ripples when concomittent noise crosses high detection threshold
	previous = 1;
	for i = 1:size(ripples,1),
		j = FindInInterval(noise,[ripples(i,1),ripples(i,3)],previous);
		previous = j(2);
		if any(normalizedSquaredNoise(j(1):j(2))>highThresholdFactor),
			excluded(i) = 1;
		end
	end
	bad = ripples(excluded,:);
	ripples = ripples(~excluded,:);
	disp(['After noise removal: ' num2str(size(ripples,1)) ' events.']);
end

% Optionally, plot results
if strcmp(show,'on'),
	figure;
	if ~isempty(noise),
		MultiPlotXY([time signal],[time squaredSignal],[time normalizedSquaredSignal],[time noise(:,2)],[time squaredNoise],[time normalizedSquaredNoise]);
		nPlots = 6;
		subplot(nPlots,1,3);
 		ylim([0 highThresholdFactor*1.1]);
		subplot(nPlots,1,6);
  		ylim([0 highThresholdFactor*1.1]);
	else
		MultiPlotXY([time signal],[time squaredSignal],[time normalizedSquaredSignal]);
%  		MultiPlotXY(time,signal,time,squaredSignal,time,normalizedSquaredSignal);
		nPlots = 3;
		subplot(nPlots,1,3);
  		ylim([0 highThresholdFactor*1.1]);
	end
	for i = 1:nPlots,
		subplot(nPlots,1,i);
		hold on;
  		yLim = ylim;
		for j=1:size(ripples,1),
			plot([ripples(j,1) ripples(j,1)],yLim,'g-');
			plot([ripples(j,2) ripples(j,2)],yLim,'k-');
			plot([ripples(j,3) ripples(j,3)],yLim,'r-');
			if i == 3,
				plot([ripples(j,1) ripples(j,3)],[ripples(j,4) ripples(j,4)],'k-');
			end
		end
		for j=1:size(bad,1),
			plot([bad(j,1) bad(j,1)],yLim,'k-');
			plot([bad(j,2) bad(j,2)],yLim,'k-');
			plot([bad(j,3) bad(j,3)],yLim,'k-');
			if i == 3,
				plot([bad(j,1) bad(j,3)],[bad(j,4) bad(j,4)],'k-');
			end
		end
		if mod(i,3) == 0,
			plot(xlim,[lowThresholdFactor lowThresholdFactor],'k','linestyle','--');
			plot(xlim,[highThresholdFactor highThresholdFactor],'k-');
		end
	end
end

function y = Filter0(b,x)

if size(x,1) == 1,
	x = x(:);
end

if mod(length(b),2) ~= 1,
	error('filter order should be odd');
end

shift = (length(b)-1)/2;

[y0 z] = filter(b,1,x);

y = [y0(shift+1:end,:) ; z(1:shift,:)];

function [U,stdA] = unity(A,sd,restrict)

if ~isempty(restrict),
	meanA = mean(A(restrict));
	stdA = std(A(restrict));
else
	meanA = mean(A);
	stdA = std(A);
end
if ~isempty(sd),
	stdA = sd;
end

U = (A - meanA)/stdA;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function res = norm1(X)
	%% ================== File info ==========================
	% Author		: Tiep Vu (http://www.personal.psu.edu/thv102/)
	% Time created	: Tue Jan 26 22:32:25 2016
	% Last modified	: Tue Jan 26 22:32:27 2016
	% Description	:
	%	return norm 1 of the input matrix X
	%% ================== end File info ==========================
	res = sum(abs(vec(X)));
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
load('ROIs/CC')

h = guidata(1);
coords = h.rois(h.curRoi).coords;
% Re-index starting at 1
coords = coords-repmat(min(coords),size(coords,1),1)+1;
coords = round(coords);
voxels = zeros(max(coords));
ind = sub2ind(size(voxels), coords(:,1), coords(:,2), coords(:,3));
voxels(ind) = 1;
voxels = smooth3(voxels,'gaussian',5);
voxels = voxels>.2;

voxels = dtiCleanImageMask(voxels, 5);
mm = [1 1 1];
id = -1;
[id,wasOpen] = mrmCheckMeshServer(id, 'localhost');
[mesh,lights,tenseMesh] = mrmBuildMesh(uint8(voxels), mm, 'localhost', id);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [struct_irf_record D_record gamma_record]=mairfres(beta_gibbs,sigma_gibbs,It,Bu,IRFperiods,n,p,k1,signrestable,signresperiods)

% function [struct_irf_record D_record gamma_record Qdraw Qsuccess]=mairfres(beta_gibbs,sigma_gibbs,It,Bu,IRFperiods,n,p,k1,signrestable,signresperiods)
% runs the gibbs sampler to obtain draws from the posterior distribution of IRFs, orthogonalised with a sign restriction setting
% inputs:  - matrix 'beta_gibbs': record of the gibbs sampler draws for the beta vector
%          - matrix 'sigma_gibbs': record of the gibbs sampler draws for the sigma matrix (vectorised)
%          - integer 'It': total number of iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'Bu': number of burn-in iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'IRFperiods': number of periods for IRFs
%          - integer 'n': the number of endogenous variables in the model
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k1': the number of coefficients related to the endogenous variables for each equation in the model
%          - cell 'signrestable': table recording the sign restriction input from the user
%          - cell 'signresperiods': table containing the periods corresponding to each restriction
% outputs: - cell 'struct_irf_record': record of the gibbs sampler draws for the orthogonalised IRFs
%          - matrix 'D_record': record of the gibbs sampler draws for the structural matrix D
%          - matrix 'gamma_record': record of the gibbs sampler draws for the structural disturbances variance-covariance matrix gamma
%          - integer 'Qdraw': total number of draws of the Q matrix
%          - integer 'Qsuccess': number of successful draws of the Q matrix

% preliminary tasks
% create first the cell that will store the results from the simulations
struct_irf_record=cell(n,n);
% storage cell
storage1=cell(It-Bu,1);
storage2=cell(It-Bu,1);

% now identify all the periods concerned with restrictions
% first expand the non-empty entries in signresperiods since they are only expressed in intervals: transform into list
% for instance, translate [1 4] into [1 2 3 4]; I don't think this can done without a loop
temp=cell2mat(signresperiods(~cellfun(@isempty,signresperiods)));
periods=[];
for ii=1:size(temp,1)
periods=[periods temp(ii,1):temp(ii,2)];
end
% suppress duplicates and sort
periods=sort(unique(periods))';
% count the total number of restriction periods (required for IRF matrix)
nperiods=size(periods,1);

% Identify the restriction matrices
% create five cells, corresponding to the three possible restrictions:
% one cell for sign restrictions, three cells for magnitude restrictions, one cell for zero restrictions
Scell=cell(1,n);
Mcell=cell(1,n);
Mlcell=cell(1,n);
Mucell=cell(1,n);
Zcell=cell(1,n);

% Check if value and periods restrictions correspond to each other
if sum(sum(~cellfun(@isempty,signresperiods) == ~cellfun(@isempty,signrestable))) == n^2
    % All cells with sign restrictions also specify the horizon over which
    % these are applied
else
    disp('Warning: Value restrictions do not correspond to period restrictions one to one')
    pause(1)
end

% loop over rows and columns of the period matrix
for ii=1:n
   for jj=1:n
      % if entry (ii,jj) of the period matrix is not empty...
      if ~isempty(signresperiods{ii,jj}) && ~isempty(signrestable{ii,jj})
      % ... then there is a restriction over one (or several) periods
      % loop overt those periods
         for kk=signresperiods{ii,jj}(1,1):signresperiods{ii,jj}(1,2)
         % identify the position of the considered period within the list of all periods (required to build the matrix)
         position=find(periods==kk);
         % now create the restriction matrix: this will depend on the type of restriction
            % if it is a positive sign restriction...
            if strcmp(signrestable{ii,jj},'+')
            % ... then input a 1 entry in the corresponding S matrix
            Scell{1,jj}=[Scell{1,jj};zeros(1,n*nperiods)];
            Scell{1,jj}(end,(position-1)*n+ii)=1;
            % if it is a negative sign restriction...
            elseif strcmp(signrestable{ii,jj},'-')
            % ... then input a -1 entry in the corresponding S matrix
            Scell{1,jj}=[Scell{1,jj};zeros(1,n*nperiods)];
            Scell{1,jj}(end,(position-1)*n+ii)=-1;
            % if it is a zero restriction...
            elseif strcmp(signrestable{ii,jj},'0')
            % ... then input a 1 entry in the corresponding Z matrix
            Zcell{1,jj}=[Zcell{1,jj};zeros(1,n*nperiods)];
            Zcell{1,jj}(end,(position-1)*n+ii)=1;
            % else, a non-empty entry being neither a sign nor a zero restriction has to be a magnitude restriction
            else
            % fill the corresponding M matrices:
            % input a 1 in M
            Mcell{1,jj}=[Mcell{1,jj};zeros(1,n*nperiods)];
            Mcell{1,jj}(end,(position-1)*n+ii)=1;
            % input the lower value of the interval in Ml
            temp=str2num(signrestable{ii,jj});
            Mlcell{1,jj}=[Mlcell{1,jj};temp(1,1)];
            % input the upper value of the interval in Mu
            Mucell{1,jj}=[Mucell{1,jj};temp(1,2)];
            end
         end
      end
   end
end

% now check what kind of restrictions apply among sign, zero and magnitude restrictions
% check for sign restrictions: if there are any, at least one entry in the cell Scell is non-empty
if sum(~cellfun(@isempty,Scell))~=0
signres=1;
else
signres=0;
end
% similarly check for zero restrictions
if sum(~cellfun(@isempty,Zcell))~=0
zerores=1;
else
zerores=0;
end
% and finally, check for magnitude restrictions
if sum(~cellfun(@isempty,Mcell))~=0
magnres=1;
else
magnres=0;
end

not_successful = 0;
hbar = bear.parfor_progressbar(It-Bu,'Progress of Sign Restriction Draws');  %create the progress bar

% initiate Gibbs algorithm
parfor ii=1:It-Bu
% initiate the variable 'success'; this variable will be used to check whether the restrictions are satisfied
success=0;

   while success==0
   not_successful = not_successful+1;
   % switch 'success' to 1; it will be turned back to zero if at any time Q is detected as a candidate not satisfying the restrictions
   success=1;

   % draw the vector of VAR coefficients
   % select first a draw index randomly
   index=floor(rand*(It-Bu))+1;
   % then draw a random set of beta and sigma corresponding to this index (this is done to make it possible to draw, if required, an infinite number of values from the gibbs sampler record, with equal probability on each value)
   B=reshape(beta_gibbs(:,ii),k1,n);
   sigma=reshape(sigma_gibbs(:,ii),n,n);
   hsigma=chol(bear.nspd(sigma),'lower');
   % obtain orthogonalised IRFs
   [irfmatrix ortirfmatrix]=bear.mairfsim(B,hsigma,p,n,max(IRFperiods,max(periods)));

   % generate the stacked IRF matrix
   stackedirfmat=[];
      for jj=1:numel(periods)
      stackedirfmat=[stackedirfmat;ortirfmatrix(:,:,periods(jj,1)+1)];
      end
   % draw an entire random matrix Q satisfying the zero restrictions
   [Q]=bear.qzerores(n,Zcell,stackedirfmat);

   % generate the candidate matrix of structural IRFs
   candidate=stackedirfmat*Q;

   % verify the restrictions in turn (sign, magnitude, or both)
   % the zero restrictions don't have to be verified, there are satisfied by construction
      % consider first the case of pure sign restrictions
      if signres==1 && magnres==0
         % loop over structural shocks; stop as soon as the draw is detected as a fail
         jj=1;
         while success==1 && jj<=n
            % if the corresponding entry in Scell is not empty, it contains a restriction to be checked
            if ~isempty(Scell{1,jj})
               % check if the restrictions hold
               if all(Scell{1,jj}*candidate(:,jj)>=0)
               % if the restrictions do not hold, there may still be a possibility by switching the sign of the Q column
               elseif all(Scell{1,jj}*((-1)*candidate(:,jj))>=0)
               Q(:,jj)=-Q(:,jj);
               % else, if there is no way to have Q succesful, count it as a fail and switch the variable success to 0
               else
               success=0;
               end
            end
         jj=jj+1;
         end
      % consider now the case of pure magnitude restrictions
      elseif signres==0 && magnres==1
         % loop over structural shocks
         jj=1;
         while success==1 && jj<=n
            % if the corresponding entry is not empty, it contains a restriction to be checked
            if ~isempty(Mcell{1,jj})
               % check if the restriction holds
               if all((Mcell{1,jj}*candidate(:,jj)-Mlcell{1,jj}).*(Mucell{1,jj}-Mcell{1,jj}*candidate(:,jj))>=0)
               % if they do not hold, there may still be a possibility by switching the sign of the Q column
               elseif all((Mcell{1,jj}*(-1)*candidate(:,jj)-Mlcell{1,jj}).*(Mucell{1,jj}-Mcell{1,jj}*(-1)*candidate(:,jj))>=0)
               Q(:,jj)=-Q(:,jj);
               % else, if there is no way to have Q succesful, count it as a fail and switch the variable success to 0
               else
               success=0;
               end
            end
         jj=jj+1;
         end
      % consider then the case of mixed sign and magnitude restrictions
      elseif signres==1 && magnres==1
         % loop over structural shocks
         jj=1;
         while success==1 && jj<=n
         % for a given structural shock, there may be no restrictions, only one type, or both types
            % if both types
            if ~isempty(Scell{1,jj}) && ~isempty(Mcell{1,jj})
               % check both restrictions
               if all(Scell{1,jj}*candidate(:,jj)>=0) && all((Mcell{1,jj}*candidate(:,jj)-Mlcell{1,jj}).*(Mucell{1,jj}-Mcell{1,jj}*candidate(:,jj))>=0)
               % if not, try to switch the sign of the corresponding column in Q
               elseif all(Scell{1,jj}*((-1)*candidate(:,jj))>=0) && all((Mcell{1,jj}*(-1)*candidate(:,jj)-Mlcell{1,jj}).*(Mucell{1,jj}-Mcell{1,jj}*(-1)*candidate(:,jj))>=0)
               Q(:,jj)=-Q(:,jj);
               % else, if there is no way to have Q succesful, count it as a fail and switch the variable success to 0
               else
               success=0;
               end
            % if only sign restrictions
            elseif ~isempty(Scell{1,jj}) && isempty(Mcell{1,jj})
               % check if the restrictions hold
               if all(Scell{1,jj}*candidate(:,jj)>=0)
               % if the restrictions do not hold, there may still be a possibility by switching the sign of the Q column
               elseif all(Scell{1,jj}*((-1)*candidate(:,jj))>=0)
               Q(:,jj)=-Q(:,jj);
               % else, if there is no way to have Q succesful, count it as a fail and switch the variable success to 0
               else
               success=0;
               end
            % if only magnitude restrictions
            elseif isempty(Scell{1,jj}) && ~isempty(Mcell{1,jj})
               % check if the restriction holds
               if all((Mcell{1,jj}*candidate(:,jj)-Mlcell{1,jj}).*(Mucell{1,jj}-Mcell{1,jj}*candidate(:,jj))>=0)
               % if they do not hold, there may still be a possibility by switching the sign of the Q column
               elseif all((Mcell{1,jj}*(-1)*candidate(:,jj)-Mlcell{1,jj}).*(Mucell{1,jj}-Mcell{1,jj}*(-1)*candidate(:,jj))>=0)
               Q(:,jj)=-Q(:,jj);
               % else, if there is no way to have Q succesful, count it as a fail and switch the variable success to 0
               else
               success=0;
               end
            end
         jj=jj+1;
         end
      % finally, the only possible remaining case is that of pure zero restrictions, satisfied by construction: no need to check
      else
      end

   % now, if sucess is still equal to 1, it means that the draw was successful for all the restrictions: keep it
   % otherwise, if success has been switched to 0, there was at least one failure: discard the draw and try with a new draw
   end

   % store
   for jj=1:IRFperiods
   storage1{ii,1}(:,:,jj)=ortirfmatrix(:,:,jj)*Q;
   end
   storage2{ii,1}=hsigma*Q;

   hbar.iterate(1);   % update progress by one iteration

end

close(hbar);   %close progress bar

% reorganise storage
% loop over iterations
for ii=1:It-Bu
   % loop over IRF periods
   for jj=1:IRFperiods
      % loop over variables
      for kk=1:n
         % loop over shocks
         for ll=1:n
         struct_irf_record{kk,ll}(ii,jj)=storage1{ii,1}(kk,ll,jj);
         end
      end
   end
D_record(:,ii)=storage2{ii,1}(:);
gamma_record(:,ii)=bear.vec(eye(n));
end

fprintf('Accepted Draws in Percent of Total Number of Draws: %f', 100*(It-Bu)/(not_successful + It-Bu))


text:
function [xs,phid,wlsden,timeiter,niter] = wls_cd(A, W, yy, xhat, niter, tol)
%	weighted least squares minimized by coordinate descent algorithm
%	Input
%		A	[nn,np]		system matrix
%		W	[nn,nn]		inverse data covariance
%		yy	[nn,1]		noisy data
%		xhat	[np,1] or [nx,ny]	initial estimate
%		tol			stopping criterion for convergence
%		niter			max # of iterations
%
%	Output
%		xs	[np,niter+1]	estimates each iteration
%		phid	[niter+1,1]	objective decrease
%		wlsden	[np,1]		curvature of surrogate
%		timeiter [niter+1,1]	CPU time each iteration
%		niter			# total iterations
%
% Saowapak Sotthivirat
% 8/21/00

[nx,ny]=size(xhat);
np = ncol(A);

yy = yy(:);
xhat = xhat(:);
xs(:,1) = xhat;

% Precompute the curvature
wlsden = diag(A'*W*A);

wii = diag(W);
res = A * xhat - yy;
phi0 = res'*W*res/2;
ii = 1;
phid(ii) = 0;
t0 = cputime;

while phid(ii) < tol & ii <= niter

	for j=1:np
		xjold = xhat(j);
		xjnew = xjold;
		phijdot = sum(A(:,j).*wii.*res);
		diff = phijdot/wlsden(j);
		xjnew = xjnew - diff;
		xhat(j) = max(0,xjnew);
		res = res - A(:,j) * diff;
	end

	ii = ii+1;
	xs(:,ii) = xhat;
	phid(ii) = phi0 - res'*W*res/2;
	timeiter(ii) = cputime - t0;
end
niter = ii;

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function f_lim = sphArrayNoiseThreshold(R, Nmic, maxG_db, maxN, arrayType, dirCoeff)
%SPHARRAYNOISETHRESHOLD Returns freq.limits for noise amplification in an SMA
%
%   SPHARRAYNOISETHRESHOLD returns the frequencies that the noise in the
%   output channels of a SMA, after performing the SHT and equalization of
%   the output signals, reaches a certain user-defined threshold maxG_db.
%   The frequencies are computed only at the lower range of each order,
%   where its response decays rapidly, ignoring for example the nulls of an
%   open array at the higher frequencies. The estimation of the limits are
%   based on a linear approximation of the log-log response found e.g. in
%
%       Sector-based Parametric Sound Field Reproduction in the Spherical Harmonic Domain
%       A Politis, J Vilkamo, V Pulkki
%       IEEE Journal of Selected Topics in Signal Processing 9 (5), 852 - 866
%
%   Inputs:
%         R:          array radius
%         Nmic:       number of microphones
%         maxG_db:    max allowed amplification for the noise level,
%               maxG_db = 20*log10(maxG)
%         maxN:       maximum order that the array supports
%         arrayType:  'open', 'rigid' or 'directional' for an open or rigid spherical
%               array of omnidirectional microphones, or for an array of
%               first-order directional microphones
%         dirCoeff: (optional) if arrayType=='directional', then dirCoeff
%               expresses the first-order directivity of microphones
%               from 0-1 (1:omni, 0.5:cardioid, 0:dipole)
%
%         f_lim:    frequency points that the threhsold is reached,
%               for orders 1:maxN
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SPHARRAYNOISETHRESHOLD.M - 11/7/2013
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<6
    dirCoeff = [];
end

c = 343;
f_lim = zeros(1,maxN);
for n=1:maxN
    bn = sphModalCoeffs(n, 1, arrayType, dirCoeff)/(4*pi);
    bn = bn(end);
    maxG = 10^(maxG_db/10);
    kR_lim = (maxG*Nmic*abs(bn).^2)^(-10*log10(2)/(6*n));
    f_lim(n) = kR_lim*c/(2*pi*R);
end

end

author: polarch

repo: Spherical-Array-Processing


text:
function [c, mass] = masscenter(F)
%------------------------------------------------------------------------------
%
% This function computes the center of mass of gridfunction F (seen as a
% density distribution) .
%
% Beware if F assumes not mere positive values, does it make sense?
%
% See also: m00, m10, m01, mu10, mu01
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: December 11, 2000.
%  2000 Stichting CWI, Amsterdam
%------------------------------------------------------------------------------
denomina = m00(F);
if denomina == 0
  error(' masscenter - demoninator vanishes ')
else
  cx = m10(F)/denomina;
  cy = m01(F)/denomina;
end
% whether cx, cy < 0 etc. could be checked here (see above warning).
if nargout == 1
  c = [cx cy];
elseif nargout == 2
  c = [cx cy];
  mass = denomina;
else
  error(' masscenter - wrong number of output arguments ')
end
%------------------------------------------------------------------------------


text:
function obj = updateStemseries(obj,dataIndex)

%-store original stem handle-%
stem_group = obj.State.Plot(dataIndex).Handle;

%-get children-%
stem_child = get(stem_group ,'Children');

%------------------------------------------------------------------------%

%-update line-%
obj.State.Plot(dataIndex).Handle = stem_child(1);
updateLineseries(obj,dataIndex);
stem_temp_data = obj.data{dataIndex};

%------------------------------------------------------------------------%

%-scatter mode-%
stem_temp_data.mode = 'lines+markers';

%------------------------------------------------------------------------%

%-update marker-%
obj.State.Plot(dataIndex).Handle = stem_child(2);
updateLineseries(obj,dataIndex);

%------------------------------------------------------------------------%

stem_temp_data.marker = obj.data{dataIndex}.marker;

%-hide every other marker-%
color_temp = cell(1,length(stem_temp_data.x));
line_color_temp = cell(1,length(stem_temp_data.x));

for n = 1:3:length(stem_temp_data.x)
    color_temp{n} = 'rgba(0,0,0,0)';
    color_temp{n+1} = stem_temp_data.marker.color;
    color_temp{n+2} = 'rgba(0,0,0,0)';
    line_color_temp{n} = 'rgba(0,0,0,0)';
    line_color_temp{n+1} = stem_temp_data.marker.line.color;
    line_color_temp{n+2} = 'rgba(0,0,0,0)';
end

% add new marker/line colors
stem_temp_data.marker.color = color_temp;
stem_temp_data.marker.line.color = line_color_temp;

%------------------------------------------------------------------------%

%-revert handle-%
obj.State.Plot(dataIndex).Handle = stem_group;
obj.data{dataIndex} = stem_temp_data;

end

author: plotly


text:
function [ t, det, info ] = ctrdi ( t, ldt, n, job )

%*****************************************************************************80
%
%% CTRDI computes the determinant and inverse of a complex triangular matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complex T(LDT,N), the triangular matrix.  The zero
%    elements of the matrix are not referenced, and the corresponding
%    elements of the array can be used to store other information.
%
%    Input, integer LDT, the leading dimension of T.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer JOB.
%    010, no determinant,    inverse, matrix is lower triangular.
%    011, no determinant,    inverse, matrix is upper triangular.
%    100,    determinant, no inverse.
%    110,    determinant,    inverse, matrix is lower triangular.
%    111,    determinant,    inverse, matrix is upper triangular.
%
%    Output, complex T(LDT,N), if an inverse was requested, then T has
%    been overwritten by its inverse.
%
%    Output, complex DET(2), the determinant of the original matrix,
%    if requested.  Otherwise not referenced.
%    Determinant = DET(1) * 10.0**DET(2) with 1.0 <= cabs1 ( DET(1) ) < 10.0
%    or DET(1) == 0.0.  Also, DET(2) is strictly real.
%
%    Output, integer INFO.
%    0, an inverse was requested and the matrix is nonsingular.
%    K, an inverse was requested, but the K-th diagonal element
%    of T is zero.
%
  det = [];
  info = 0;

  if ( floor ( job / 100 ) ~= 0 )

    det(1) = 1.0;
    det(2) = 0.0;

    for i = 1 : n

      det(1) = det(1) * t(i,i);

      if ( cabs1 ( det(1) ) == 0.0 )
        break
      end

      while ( cabs1 ( det(1) ) < 1.0 )
        det(1) = det(1) * 10.0;
        det(2) = det(2) - 1.0;
      end

      while ( 10.0 <= cabs1 ( det(1) ) )
        det(1) = det(1) / 10.0;
        det(2) = det(2) + 1.0;
      end

    end

  end
%
%  Compute inverse of upper triangular matrix.
%
  if ( mod ( floor ( job / 10 ), 10 ) ~= 0 )

    if ( mod ( job, 10 ) ~= 0 )

      info = 0;

      for k = 1 : n

        if ( cabs1 ( t(k,k) ) == 0.0 )
          info = k;
          break
        end

        t(k,k) = 1.0 / t(k,k);
        temp = -t(k,k);
        t(1:k-1,k) = t(1:k-1,k) * temp;

        for j = k+1 : n
          temp = t(k,j);
          t(k,j) = 0.0;
          t(1:k,j) = t(1:k,j) + temp * t(1:k,k);
        end

      end
%
%  Compute inverse of lower triangular matrix.
%
    else

      info = 0;

      for k = n : -1: 1

        if ( cabs1 ( t(k,k) ) == 0.0 )
          info = k;
          break
        end

        t(k,k) = 1.0 / t(k,k);

        if ( k ~= n )
          temp = -t(k,k);
          t(k+1:n,k) = t(k+1:n,k) * temp;
        end

        for j = 1 : k-1
          temp = t(k,j);
          t(k,j) = 0.0;
          t(k:n,j) = t(k:n,j) + temp * t(k:n,k);
        end

      end

    end

  end

  return
end


text:
% a basic 3D ZTE/PETRA sequence
% achieves "TE" about 70 us and possibly below

%% high-level sequence parameters

fov=256e-3; dx=2e-3;            % Define FOV and resolution
alpha=4;                        % flip angle
Nr=300;                         % number of readout points (with some oversampling)
R= 8;                           % acceleration/undersampling (for the outer shell)
R_inner= 1;                     % acceleration/undersampling (angular direction of the inner area)
xSpoil=3;%0.6;  %0.6 for 2mm       % the amount of spoiling after the end of the readout (used to ramp to the next point)

Kmax=1/2/dx;
dK=1/fov;

%% more detailed and derived params

rf_duration=10e-6;              % duration of the excitation pulse
ro_duration=300e-6;             % read-out time: controls RO bandwidth and T2-blurring
minRF_to_ADC_time=50e-6;        % the parameter wich defines TE together with ro_discard
%ro_discard=4;                   % how many ADC samples are contaminated by RF switching artifacts and alike
rfSpoilingInc=117;              % RF spoiling increment

% system limits
sys = mr.opts('MaxGrad', 36, 'GradUnit', 'mT/m', ...
    'MaxSlew', 180, 'SlewUnit', 'T/m/s', 'rfRingdownTime', 10e-6, ...
    'rfDeadTime', 100e-6, 'adcDeadTime', 10e-6, 'gamma',42.576e6); % 1H: 42.576e6  23Na: 11.262e6

seq=mr.Sequence(sys);           % Create a new sequence object
seq_sar=mr.Sequence(sys);       % Create an auxillary sequence object for SAR testing

%% create main sequence elements

% %create alpha-degree block pulse
% rf = mr.makeBlockPulse(alpha*pi/180,'Duration',rf_duration,'system',sys);
%or create alpha-degree gaussian pulse
rf = mr.makeGaussPulse(alpha*pi/180,'Duration',rf_duration,'timeBwProduct',3,'system',sys);
Tenc=rf_duration/2+minRF_to_ADC_time+ro_duration; %encoding time
Tg=sys.rfDeadTime+rf_duration/2+Tenc+sys.adcDeadTime; % constant gradient time
Tt=ceil((Tenc*(1+xSpoil)-Tg)/sys.gradRasterTime)*sys.gradRasterTime; % transition time
Ag=Kmax/Tenc; % read gradient grdient amplitude
% % Gr=mr.makeTrapezoid('z','Amplitude',Ag,'flatTime',Tg,'riseTime',0); % this graient has no ramps, I wonder if the Matlab mr library would like it...
% Gr=mr.makeExtendedTrapezoid('z','times',[0 Tg],'amplitudes',[Ag Ag]); % this is a constant graient with no ramps
TR=Tg+Tt;

adc=mr.makeAdc(Nr,'Duration', ro_duration, 'Delay', sys.rfDeadTime+rf_duration+minRF_to_ADC_time);

SamplesBookkeeping=[];

%rfbw=1/rf_duration;
rfbw=mr.calcRfBandwidth(rf);
fprintf('Read gradient amplitude %g mT/m, effective "slice thinckess" %f mm\n', 1e3*Ag/sys.gamma, 1/Ag*rfbw*1e3);
FO=50e-3*Ag;
NF=0; % TR is broken now for NF~=0...

%% generate the sampling set on a surface of a sphere
[phi, theta, im]=spherical_samples(Kmax,dK,R);
Ns=length(phi);
SamplesBookkeeping=[SamplesBookkeeping Ns];

%% main ZTE loop
fprintf('Populating ZTE loop (%d TRs)\n', Ns);
tic;
populate_subsequence(sys,seq,rf,adc,phi,theta,im,Ns,Ag,TR,Tt,NF,FO);
toc

%% SPI loop
KstartZTE=Ag*(rf_duration/2+minRF_to_ADC_time+adc.dwell);
nKspi=floor(KstartZTE/dK);
dKspi=KstartZTE/(nKspi+1);
Tenc=rf_duration/2+minRF_to_ADC_time+adc.dwell/2;
fprintf('Populating SPI loop (%d spheres)\n', nKspi);
tic;
for s=nKspi:-1:1
    % generate the sampling set on a surface of a sphere
    [phi, theta, im]=spherical_samples(dKspi*s,dKspi,R_inner); % normally no acceleration
    Ns=length(phi);
    SamplesBookkeeping=[SamplesBookkeeping Ns];
    % update gradients
    Ag=dKspi*s/Tenc; % read gradient grdient amplitude
    % actually create the next sampling sphere
    fprintf('Populating sphere %d (%d TRs)\n',s,Ns);
    fprintf('Effective "slice thinckess" %f mm\n',  1/Ag*rfbw*1e3);
    populate_subsequence(sys,seq,rf,adc,phi,theta,im,Ns,Ag,TR,Tt,NF,FO);
end
% sample the centere of k-space
seq.addBlock(mr.makeDelay(Tt));
seq.addBlock(rf, adc, mr.makeDelay(TR-Tt));
SamplesBookkeeping=[SamplesBookkeeping 1];
toc
fprintf('Total number of SPI samples: %d; a Cartesian patch would require %d\n', sum(SamplesBookkeeping(2:end)), ceil(nKspi^3*pi*4/3));

% %% OLD!!! SPI loop
% Kspi=ceil(Ag*(rf_duration/2+minRF_to_ADC_time+adc.dwell)/dK);
% rf.delay=rf.delay+Tt;
% Tspi1=mr.calcDuration(rf)+sys.rfRingdownTime;
% adc.delay=adc.delay+Tt-Tspi1;
% Tspi2=TR-Tspi1;
%
% g=mr.makeTrapezoid('x','system',sys,'area',dK*Kspi,'duration',minRF_to_ADC_time);
% gx=g;gx.channel='x';
% gy=g;gy.channel='y';
% gz=g;gz.channel='z';
%
% gxr=g;gxr.channel='x';
% gyr=g;gyr.channel='y';
% gzSpoil=mr.makeTrapezoid('x','system',sys,'area',Kmax*(1+xSpoil),'duration',minRF_to_ADC_time);
%
% fprintf('Populating SPI loop (~%d TRs)\n', round(pi/6*((2*Kspi+1)^3)));
% tic;
%
% % SPI loop itself
% gxr.amplitude=0;
% gyr.amplitude=0;
% for k=-Kspi:Kspi
%     for j=-Kspi:Kspi
%         for i=-Kspi:Kspi
%             if i^2+j^2+k^2>Kspi^2
%                 continue;
%             end
%             % refocusing and spoiling for the previous TR
%             seq.addBlock( rf, gxr, gyr, gzSpoil, mr.makeDelay(Tspi1));
%             % new TR
%             gx.amplitude=g.amplitude/Kspi*i;
%             gy.amplitude=g.amplitude/Kspi*j;
%             gy.amplitude=g.amplitude/Kspi*k;
%             seq.addBlock( rf, adc, gx, gy, gz, mr.makeDelay(Tspi2));
%             gxr.amplitude=-gx.amplitude;
%             gyr.amplitude=-gy.amplitude;
%         end
%     end
% end
% toc

%% check whether the timing of the sequence is correct
[ok, error_report]=seq.checkTiming;

if (ok)
    fprintf('Timing check passed successfully\n');
else
    fprintf('Timing check failed! Error listing follows:\n');
    fprintf([error_report{:}]);
    fprintf('\n');
end

%%
%seq.plot('TimeRange',[0 1]);

seq.setDefinition('FOV', [fov fov fov]);
seq.setDefinition('Name', 'petra');
seq.setDefinition('SamplesPerShell', SamplesBookkeeping);

seq.write('zte_petra.seq');
return

%% create an RF-only version of the sequence (e.g. for the SAR or signal evolution testing)

tic;
[total_duration, total_numBlocks]=seq.duration();
for iB=1:total_numBlocks
    b=seq.getBlock(iB);
    bd=seq.blockDurations(iB);
    bs={mr.makeDelay(bd)};
    if ~isempty(b.rf)
        bs{end+1}=b.rf;
    end
    if ~isempty(b.adc)
        bs{end+1}=b.adc;
    end
    seq_sar.addBlock(bs);
end
toc
seq_sar.write('zte_petra_sar.seq');
return

% %% test binary storing
%
% seq.writeBinary('zte_petra.bin');
% seq_bin=mr.Sequence();
% seq_bin.readBinary('zte_petra.bin');
% seq_bin.write('zte_petra_bin.seq');
% return

%% visualize the 3D k-space
tic;
[kfa,~,kf]=seq.calculateKspacePP();
toc

%
figure;plot3(kf(1,:),kf(2,:),kf(3,:));
hold on;plot3(kfa(1,:),kfa(2,:),kfa(3,:),'r.');

%% local functions
function [phi, theta, im]=spherical_samples(Kr,dK,R)
% the number of samples equals the ceil of the area of the sphere divided
% by the area around every sample
Ns=ceil(4*pi*((Kr/dK)^2)/R);
np=0:(Ns-1);
alpha_gold=pi*(3-sqrt(5));
phi=np*alpha_gold;
%theta=0.5*pi*sqrt(np/Ns); % from Davide Piccini  https://doi.org/10.1002/mrm.22898
theta=acos(1-2*np/(Ns-1));  % from  Anton Semechko (2020). Suite of functions to perform uniform sampling of a sphere (https://github.com/AntonSemechko/S2-Sampling-Toolbox), GitHub. Retrieved October 3, 2020.

xp=sin(theta).*cos(phi);
yp=sin(theta).*sin(phi);
zp=cos(theta);
%figure; sphere; colormap gray;
%hold on; plot3(xp,yp,zp,'.');

% looking for the optimal interleaving factor
nm=round(Ns/2); % middle of the trajetory
sr=round(sqrt(Ns));
v0=[xp(nm); yp(nm); zp(nm)];
v=[xp(nm+(1:sr)); yp(nm+(1:sr)); zp(nm+(1:sr))];
%figure;plot(vecnorm((v-v0(:,ones(size(v,2),1)))));
[dKm,im]=min(vecnorm((v-v0(:,ones(size(v,2),1)))));

% fprintf('requested dK=%g achieved minimal dK=%g, min acceleration: %g, Ns=%d\n', dK/sqrt(R), dKm*Kr, (dKm*Kr/dK)^2,Ns);

% v=[xp; yp; zp]*Kr/dK/sqrt(R)*100;
% md=zeros(1,Ns);
% d=zeros(1,Ns);
% for i=1:Ns
%     d=vecnorm(v(:,i*ones(1,Ns))-v);
%     d(i)=NaN;
%     md(i)=min(d);
% end
% fprintf('dKmin=%g%%, dKmed=%g%% dKmax=%g%%\n', min(md),median(md),max(md));

end

function populate_subsequence(sys,seq,rf,adc,phi,theta,im,Ns,Ag,TR,Tt,FN, FO)
Azc=Ag*(TR-Tt)/(TR+Tt);  %*0.35;
%Azc=0; % for MoCo we need "no-gradient" event blocks to be able to apply updates

Gr=mr.makeExtendedTrapezoid('z','times',[0 TR-Tt],'amplitudes',[Ag Ag]); % this is a constant graient with no ramps

% pre-ramp the gradient to Azc
if abs(Azc)>eps
    Tpr=max(2,ceil(Azc/sys.maxSlew/sys.gradRasterTime))*sys.gradRasterTime;
    assert(Tpr<=TR);
    % this "dummy TR" doe not have an RF pulse, which is not good when it is called for the inner shells... but otherwise there would be no enough spoiling...
    seq.addBlock(mr.align('right',mr.makeDelay(TR),mr.makeExtendedTrapezoid('z','system',sys,'times',[0,Tpr],'amplitudes',[0,Azc])));
end

% the loop itself
for j=1:im
    Glast=struct('x',0,'y',0,'z',Azc);
    for i=j:im:Ns
        Gcr=mr.rotate('z',phi(i),mr.rotate('y',theta(i),Gr));
        Gcurr=struct('x',0,'y',0,'z',0);
        for g=1:length(Gcr)
            Gcurr.(Gcr{g}.channel)=Gcr{g}.waveform(1);
        end
        seq.addBlock( ...
            mr.makeExtendedTrapezoid('x','system',sys,'times',[0,Tt],'amplitudes',[Glast.x,Gcurr.x]), ...
            mr.makeExtendedTrapezoid('y','system',sys,'times',[0,Tt],'amplitudes',[Glast.y,Gcurr.y]), ...
            mr.makeExtendedTrapezoid('z','system',sys,'times',[0,Tt],'amplitudes',[Glast.z,Gcurr.z]));
        for f=-FN:FN % 'slices' loop (frequency offsets)
            rf.freqOffset=f*FO;
            rf.phaseOffset=-2*pi*f*FO*rf.t(end)/2;
            seq.addBlock( [{rf, adc}, Gcr] );
        end
        Glast=Gcurr;
    end
    rf_aux=rf;
    rf_aux.delay=rf.delay+Tt;
    if (j==im)
        Azc=0; % on the last interleave we ramp down to 0
    end
    seq.addBlock( rf_aux, ...
        mr.makeExtendedTrapezoid('x','system',sys,'times',[0,TR],'amplitudes',[Glast.x,0]), ...
        mr.makeExtendedTrapezoid('y','system',sys,'times',[0,TR],'amplitudes',[Glast.y,0]), ...
        mr.makeExtendedTrapezoid('z','system',sys,'times',[0,TR],'amplitudes',[Glast.z,Azc]));
end
end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function sVF1 = minus(sVF1, sVF2)
%
% Syntax
%   sVF = sVF1 - sVF2
%   sVF = sVF - v
%
% Input
%  sVF1 - @S2VectorFieldTri
%  sVF2 - @S2VectorField
%  v    - @vector3d
%
% Output
%  sF - @S2FunTri
%

sVF1 = sVF1 + (-sVF2);

author: mtex-toolbox

repo: mtex

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


text:
function inside = p13_inside ( m, n, point )

%*****************************************************************************80
%
%% P13_INSIDE reports if a point is inside the region in problem 13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real POINT(M,N), the coordinates of the points.
%
%    Output, logical INSIDE(N), is TRUE if the point is in the region.
%
  inside(1:n) =                                           ...
       ( 45.0 <= point(1,1:n) & point(1,1:n) <= 55.0 &    ...
         30.0 <= point(2,1:n) & point(2,1:n) <= 90.0    ) ...
    |                                                     ...
       ( 900.0 <=                                         ...
         ( point(1,1:n) - 50.0 ).^2                       ...
       + ( point(2,1:n) -  0.0 ).^2 &                     ...
         ( point(1,1:n) - 50.0 ).^2                       ...
       + ( point(2,1:n) -  0.0 ).^2 <= 1600.0 );

  return
end