text:
function obj = updateLegendMultipleAxes(obj, legIndex)

% x: ...[DONE]
% y: ...[DONE]
% traceorder: ...[DONE]
% font: ...[DONE]
% bgcolor: ...[DONE]
% bordercolor: ...[DONE]
% borderwidth: ...[DONE]
% xref: ...[DONE]
% yref: ...[DONE]
% xanchor: ...[DONE]
% yanchor: ...[DONE]

%=========================================================================%
%
%-GET NECESSARY INFO FOR MULTIPLE LEGENDS-%
%
%=========================================================================%

for traceIndex = 1:obj.State.Figure.NumPlots
    allNames{traceIndex} = obj.data{traceIndex}.name;
    allShowLegens(traceIndex) = obj.data{traceIndex}.showlegend;
    obj.data{traceIndex}.showlegend = false;
    obj.data{traceIndex}.legendgroup = obj.data{traceIndex}.name;

    axIndex = obj.getAxisIndex(obj.State.Plot(traceIndex).AssociatedAxis);
    [xSource, ySource] = findSourceAxis(obj, axIndex);
    xAxis = eval(sprintf('obj.layout.xaxis%d', xSource));
    yAxis = eval(sprintf('obj.layout.yaxis%d', xSource));

    allDomain(traceIndex, 1) = max(xAxis.domain);
    allDomain(traceIndex, 2) = max(yAxis.domain);
end

[~, groupIndex] = unique(allNames);

for traceIndex = groupIndex'
   obj.data{traceIndex}.showlegend = allShowLegens(traceIndex);
end

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

%-STANDARDIZE UNITS-%
legendUnits = get(obj.State.Legend(legIndex).Handle,'Units');
fontUnits = get(obj.State.Legend(legIndex).Handle,'FontUnits');
set(obj.State.Legend(legIndex).Handle,'Units','normalized');
set(obj.State.Legend(legIndex).Handle,'FontUnits','points');

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

%-LEGEND DATA STRUCTURE-%
legendData = get(obj.State.Legend(legIndex).Handle);

% only displays last legend as global Plotly legend
obj.layout.legend = struct();

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

%-layout showlegend-%
obj.layout.showlegend = strcmpi(legendData.Visible,'on');

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

%-legend (x,y) coordenates-%
obj.layout.legend.x = 1.005 * max(allDomain(:,1));
obj.layout.legend.y = 1.001 * max(allDomain(:,2));;

%-legend (x,y) refs-%
obj.layout.legend.xref = 'paper';
obj.layout.legend.yref = 'paper';

%-legend (x,y) anchors-%
obj.layout.legend.xanchor = 'left';
obj.layout.legend.yanchor = 'top';

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

if (strcmp(legendData.Box,'on') && strcmp(legendData.Visible, 'on'))

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

    %-legend traceorder-%
    obj.layout.legend.traceorder = 'normal';

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

    %-legend borderwidth-%
    obj.layout.legend.borderwidth = legendData.LineWidth;

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

    %-legend bordercolor-%
    col = 255*legendData.EdgeColor;
    obj.layout.legend.bordercolor = sprintf('rgb(%f,%f,%f)', col);

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

    %-legend bgcolor-%
    col = 255*legendData.Color;
    obj.layout.legend.bgcolor = sprintf('rgb(%f,%f,%f)', col);

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

    %-legend font size-%
    obj.layout.legend.font.size = legendData.FontSize;

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

    %-legend font family-%
    obj.layout.legend.font.family = matlab2plotlyfont(legendData.FontName);

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

    %-legend font colour-%
    col = 255*legendData.TextColor;
    obj.layout.legend.font.color = sprintf('rgb(%f,%f,%f)', col);

end

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

%-REVERT UNITS-%
set(obj.State.Legend(legIndex).Handle,'Units',legendUnits);
set(obj.State.Legend(legIndex).Handle,'FontUnits',fontUnits);

end

author: plotly


text:
function snowfall_display ( )

%*****************************************************************************80
%
%% SNOWFALL_DISPLAY plots yearly snowfall data side by side.
%
%  Discussion:
%
%    The file "snowfall.txt" contains snowfall records.
%    Each line contains the year, the snowfall in inches for 8 months, and
%    the total.
%
%    We want to plot the yearly snowfall curves for 2000 through 2010,
%    side by side in a 3D plot.
%
%    We can use PLOT3(X,Y,Z) to plot a single 3D curve.
%
%    We can use HOLD ON / HOLD OFF to force multiple 3D curves to appear
%    together.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 November 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SNOWFALL_DISPLAY\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Demonstrate how to display several curves side by side.\n' );
  fprintf ( 1, '  Here, each curve represents snowfall over a given year.\n' );
%
%  Load all the snowfall data.
%
  data = load ( 'snowfall.txt' );
%
%  Column 1 contains the years.
%
  years = data ( :, 1 );
  offset = years(1) - 1;
%
%  Columns 2 through 9 are snowfall amounts.
%
  snow = data ( :, 2:9);
%
%  Column 10 is the snowfall total.
%
  total = data ( :, 10 );
%
%  Create figure #1 and clear it.
%
  figure ( 1 )
  clf
%
%  Call HOLD ON so multiple curves share one image.
%
%  Call PLOT3 to draw each snowfall curve.
%
  hold on
  for year = 2000 : 2010
    x(1:8) = year;
    y(1:8) = ( 1 : 8 );
    z(1:8) = snow(year-offset,1:8);
    plot3 ( x, y, z, 'Linewidth', 3 );
  end
%
%  Annotate the plot.
%
  grid on
  xlabel ( 'Year' );
  ylabel ( 'Month' );
  set ( gca, 'yTick', 1 : 8 );
  set ( gca, 'yTickLabel', { 'Oct', 'Nov', 'Dec', 'Jan', 'Feb', 'Mar', 'Apr', 'May' } );
  zlabel ( 'Inches of Snow' );
  title ( 'Michigan Tech Snowfall Records', 'Fontsize', 16 );
%
%  Set the 3D viewing angles in degrees.
%
  azimuth = 45.0;
  elevation = 35.0;
  view ( azimuth, elevation );
%
%  Call HOLD OFF because we are done concatenating plots.
%
  hold off
%
%  Save a copy of the plot as a PNG file.
%
  filename = 'snowfall.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  PNG image saved as "%s"\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Va, success] = dcpf(B, Pbus, Va0, ref, pv, pq)
%DCPF  Solves a DC power flow.
%   [VA, SUCCESS] = DCPF(B, PBUS, VA0, REF, PV, PQ) solves for the bus
%   voltage angles at all but the reference bus, given the full system
%   B matrix and the vector of bus real power injections, the initial
%   vector of bus voltage angles (in radians), and column vectors with
%   the lists of bus indices for the swing bus, PV buses, and PQ buses,
%   respectively. Returns a vector of bus voltage angles in radians.
%
%   See also RUNDCPF, RUNPF.

%   MATPOWER
%   Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC)
%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
%   and Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% constant
Va_threshold = 1e5;     %% arbitrary threshold on |Va| for declaring failure

%% initialize result vector
Va = Va0;
success = 1;    %% successful by default

%% set up to trap non-singular matrix warnings
[lastmsg, lastid] = lastwarn;
lastwarn('');

%% update angles for non-reference buses
Va([pv; pq]) = B([pv; pq], [pv; pq]) \ ...
                (Pbus([pv; pq]) - B([pv; pq], ref) * Va0(ref));

[msg, id] = lastwarn;
%% Octave is not consistent in assigning proper warning id, so we'll just
%% check for presence of *any* warning
if ~isempty(msg) || max(abs(Va)) > Va_threshold
    success = 0;
end

%% restore warning state
lastwarn(lastmsg, lastid);

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function x = gen_form(L_p,x_s,A,b,K,M)
%GEN_FORM Transform a standard-form problem back to the general-form setting.
%
% x = gen_form(L_p,x_s,A,b,K,M)    (method 1)
% x = gen_form(L_p,x_s,x_0)        (method 2)
%
% Transforms the standard-form solution x_s back to the required
% solution to the general-form problem:
%    x = L_p*x_s + d ,
% where L_p and d depend on the method as follows:
%    method = 1: L_p = pseudoinverse of L, d  = K*(b - A*L_p*x_s)
%    method = 2: L_p = A-weighted pseudoinverse of L, d = x_0.
%
% Usually, the standard-form problem is generated by means of
% function std_form.
%
% Note that x_s may have more that one column.

% References: L. Elden, "Algorithms for regularization of ill-
% conditioned least-squares problems", BIT 17 (1977), 134-145.
% L. Elden, "A weighted pseudoinverse, generalized singular values,
% and constrained lest squares problems", BIT 22 (1982), 487-502.
% M. Hanke, "Regularization with differential operators.  An itera-
% tive approach", J. Numer. Funct. Anal. Optim. 13 (1992), 523-540.

% Per Christian Hansen, IMM, 06/12/93.

% Nargin determines which method.
if (nargin==6)
  [p,q] = size(x_s); Km = size(K,1);
  if (Km==0)
    x = L_p*x_s;
  else
    x = L_p*x_s + K*(M*(b*ones(1,q) - A*(L_p*x_s)));
  end
else
  x_0 = A; [p,q] = size(x_s);
  x = L_p*x_s + x_0*ones(1,q);
end

repo: REST


text:
% --------- DeepMIMO: A Generic Dataset for mmWave and massive MIMO ------%
% Authors: Ahmed Alkhateeb, Umut Demirhan, Abdelrahman Taha
% Date: March 17, 2022
% Goal: Encouraging research on ML/DL for MIMO applications and
% providing a benchmarking tool for the developed algorithms
% ---------------------------------------------------------------------- %

function [DeepMIMO_dataset, params] = DeepMIMO_generator(params)

    % -------------------------- DeepMIMO Dataset Generation -----------------%
    fprintf(' DeepMIMO Dataset Generation started')

    [params, params_inner] = validate_parameters(params);

    if params_inner.dynamic_scenario
        for f = 1:length(params_inner.list_of_folders)
            fprintf('\nGenerating Scene %i/%i', f, length(params_inner.list_of_folders))
            params.scenario_files = fullfile(params_inner.list_of_folders{f}, params.scenario); % The initial of all the scenario files
            DeepMIMO_scene{f} = generate_data(params, params_inner);
            param{f} = params;
        end

        DeepMIMO_dataset = DeepMIMO_scene;
        params = param;
        saveDataset = params{1}.saveDataset;
    else
        DeepMIMO_dataset = generate_data(params, params_inner);
        saveDataset = params.saveDataset;
    end

    % Saving the data
    if saveDataset
        fprintf('\n Saving the DeepMIMO Dataset ...')

        fileidx = 1;
        while isfile(sprintf('DeepMIMO_dataset/dataset_%i.mat', fileidx))
            fileidx = fileidx + 1;
        end
        sfile_DeepMIMO = sprintf('DeepMIMO_dataset/dataset_%i.mat', fileidx);
        dataset_params = params;
        save(sfile_DeepMIMO,'DeepMIMO_dataset', 'dataset_params', '-v7.3');

        fprintf('\n The generated DeepMIMO dataset is saved into %s file.', sfile_DeepMIMO);

    end

    fprintf('\n DeepMIMO Dataset Generation completed \n')

end

function DeepMIMO_dataset = generate_data(params, params_inner)
    % Reading ray tracing data
    fprintf('\n Reading the channel parameters of the ray-tracing scenario %s', params.scenario)
    for t=1:params.num_active_BS
        bs_ID = params.active_BS(t);
        fprintf('\n Basestation %i', bs_ID);
        [TX{t}.channel_params, TX{t}.channel_params_BSBS, TX{t}.loc] = read_raytracing(bs_ID, params, params_inner.scenario_files);
    end

    % Constructing the channel matrices from ray-tracing
    for t = 1:params.num_active_BS
        fprintf('\n Constructing the DeepMIMO Dataset for BS %d', params.active_BS(t))
        c = progress_counter(params.num_user+params.enable_BS2BSchannels*params.num_active_BS);

        % BS transmitter location & rotation
        DeepMIMO_dataset{t}.loc = TX{t}.loc;
        DeepMIMO_dataset{t}.rotation = params_inner.array_rotation_BS(t,:);

        %----- BS-User Channels
        for user=1:params.num_user
            % Channel Construction

            if params.generate_OFDM_channels
                [DeepMIMO_dataset{t}.user{user}.channel]=construct_DeepMIMO_channel(params_inner.num_ant_BS(t, :), params_inner.array_rotation_BS(t,:), params_inner.ant_spacing_BS(t), params.num_ant_UE, params_inner.array_rotation_UE(user, :), params.ant_spacing_UE, TX{t}.channel_params(user), params);
            else
                [DeepMIMO_dataset{t}.user{user}.channel]=construct_DeepMIMO_channel_TD(params_inner.num_ant_BS(t, :), params_inner.array_rotation_BS(t,:), params_inner.ant_spacing_BS(t), params.num_ant_UE, params_inner.array_rotation_UE(user, :), params.ant_spacing_UE, TX{t}.channel_params(user), params);
                DeepMIMO_dataset{t}.user{user}.ToA = TX{t}.channel_params(user).ToA; %Time of Arrival/Flight of each channel path (seconds)
            end
            DeepMIMO_dataset{t}.user{user}.rotation = params_inner.array_rotation_UE(user, :);

            % Location, LOS status, distance, pathloss, and channel path parameters
            DeepMIMO_dataset{t}.user{user}.loc=TX{t}.channel_params(user).loc;
            DeepMIMO_dataset{t}.user{user}.LoS_status=TX{t}.channel_params(user).LoS_status;
            DeepMIMO_dataset{t}.user{user}.distance=TX{t}.channel_params(user).distance;
            DeepMIMO_dataset{t}.user{user}.pathloss=TX{t}.channel_params(user).pathloss;
            DeepMIMO_dataset{t}.user{user}.path_params=rmfield(TX{t}.channel_params(user),{'loc','distance','pathloss'});

            c.increment();
        end

        %----- BS-BS Channels
        if params.enable_BS2BSchannels
            for BSreceiver=1:params.num_active_BS
                % Channel Construction
                if params.generate_OFDM_channels
                    [DeepMIMO_dataset{t}.basestation{BSreceiver}.channel]=construct_DeepMIMO_channel(params_inner.num_ant_BS(t, :), params_inner.array_rotation_BS(t,:), params_inner.ant_spacing_BS(t), params_inner.num_ant_BS(BSreceiver, :), params_inner.array_rotation_BS(BSreceiver,:), params_inner.ant_spacing_BS(BSreceiver), TX{t}.channel_params_BSBS(BSreceiver), params);
                else
                    [DeepMIMO_dataset{t}.basestation{BSreceiver}.channel]=construct_DeepMIMO_channel_TD(params_inner.num_ant_BS(t, :), params_inner.array_rotation_BS(t,:), params_inner.ant_spacing_BS(t), params_inner.num_ant_BS(BSreceiver, :), params_inner.array_rotation_BS(BSreceiver,:), params_inner.ant_spacing_BS(BSreceiver), TX{t}.channel_params_BSBS(BSreceiver), params);
                    DeepMIMO_dataset{t}.basestation{BSreceiver}.ToA=TX{t}.channel_params_BSBS(BSreceiver).ToA; %Time of Arrival/Flight of each channel path (seconds)
                end
                DeepMIMO_dataset{t}.basestation{BSreceiver}.rotation = params_inner.array_rotation_BS(BSreceiver, :);

                % Location, LOS status, distance, pathloss, and channel path parameters
                DeepMIMO_dataset{t}.basestation{BSreceiver}.loc=TX{t}.channel_params_BSBS(BSreceiver).loc;
                DeepMIMO_dataset{t}.basestation{BSreceiver}.LoS_status=TX{t}.channel_params_BSBS(BSreceiver).LoS_status;
                DeepMIMO_dataset{t}.basestation{BSreceiver}.distance=TX{t}.channel_params_BSBS(BSreceiver).distance;
                DeepMIMO_dataset{t}.basestation{BSreceiver}.pathloss=TX{t}.channel_params_BSBS(BSreceiver).pathloss;
                DeepMIMO_dataset{t}.basestation{BSreceiver}.path_params=rmfield(TX{t}.channel_params_BSBS(BSreceiver),{'loc','distance','pathloss'});

                c.increment();
            end
        end
    end

end

author: DeepMIMO

repo: DeepMIMO-matlab


text:
classdef BoundingBoxGaterX <GaterX
% BOUNDINGBOXGATERX Class
%
% Summary of BoundingBoxGaterX:
% This is a class implementation of a bounding box gater, which gates
% measurements, given in the form of bounding boxes (x,y,w,h), depending on
% their proportional overlap with the estimated target bounding boxes.
%
% BoundingBoxGaterX Properties:
%   + OverlapThresh - The minimum overlap, proportional to the each
%                     target's bounding box.
%
% BoundingBoxGaterX Methods:
%   + BoundingBoxGaterX  - Constructor method
%   + gate - Perform gating and generate the validation matrix
%
% See also SystematicResamplerX
    properties
        OverlapThresh
    end

    properties (SetAccess = immutable)
        NumObsDims
    end

    methods
        function this = BoundingBoxGaterX(varargin)
        % BoundingBoxGaterX Constructor method
        %
        % Parameters
        % ----------
        % OverlapThresh: normalised scalar, optional
        %   The normalised ([0,1]) minimum percentage overlap between the estimated
        %   target bounding box and any given measurement, such that a measurement
        %   should be gated.
        %
        % Usage
        % -----
        % * eg = BoundingBoxGaterX(OverlapThresh) returns an object configured
        %   to exclude any measurements that cover less than OverlapThresh %
        %   of a given target state.
        % * eg = BoundingBoxGaterX('OverlapThresh', OverlapThresh) is also valid
        %
        % See also EllipsoidalGaterX/validate

            parser = inputParser;
            parser.addOptional('OverlapThresh',[]);
            parser.KeepUnmatched = true;
            parser.parse(varargin{:});
            this.OverlapThresh = parser.Results.OverlapThresh;

        end

        function [ValidationMatrix, GateVolumes] = gate(this,varargin)
        % GATE Perform bbox gating to generate the validation
        % matrix indicating the valid associations between Tracks and
        % measurements
        %
        % Parameters
        % ----------
        % TrackList: cell vector, optional
        %   A (1 x NumTracks) cell vector, where each cell contains a TrackX
        %   object. It is required that each TrackX object should have a FilterX
        %   object as a property.
        % TrackPredMeasMeans: matrix, optional
        %   A (NumObsDims x NumTracks) matrix, where each ith column represents
        %   the predicted measurement mean for the ith track. Can be used,
        %   in combination with TrackInnovErrCovars (see below), instead of
        %   a TrackList (see above).
        % TrackInnovErrCovars: matrix, optional
        %   A (NumObsDims x NumObsDims x NumTracks) matrix, where each
        %   (:,:,i) page represents the innovation covariance matrix for
        %   the ith track. Can be used, in combination with TrackInnovErrCovars
        %   (see above), instead of a TrackList (see above).
        % MeasurementList: matrix
        %   A (NumObsDims x NumMeasurements) matrix, where each jth column
        %   corresponds to the jth measurement.
        %
        % Usage
        % -----
        % * [ValidationMatrix] = gate(this,TrackList,MeasurementList)
        %   returns the (NumTracks x NumMeasurements) validation matrix
        %   ValidationMatrix indicating the valid associations between tracks
        %   in TrackList (1 x NumTracks array of TrackX objects) and measurements
        %   in MeasurmentList (NumObsDims x NumMeasurements matrix)
        % * [ValidationMatrix] = gate(this,TrackPredMeasMeans,TrackInnovErrCovars,MeasurementList)
        %   returns the (NumTracks x NumMeasurements) validation matrix
        %   ValidationMatrix indicating the valid associations between tracks
        %   in with predicted measurement means given by the (NumObsDims x NumTracks)
        %   TrackPredMeasMeans matrix, innovation covariances given by the
        %   (NumObsDims x NumObsDims x NumTracks) TrackInnovErrCovars matrix, and
        %   measurements in MeasurmentList (NumObsDims x NumMeasurements matrix), as well as a
        %   (1 x NumTracks) vector GateVolumes, containing the volume of the computed
        %   gate for each track.
        %
        % See also EllipsoidalGaterX/EllipsoidalGaterX

            if(nargin == 3)
                TrackList = varargin{1};
                MeasurementList = varargin{2};
                NumTracks = numel(TrackList);
            else
                TrackPredMeasMeans = varargin{1};
                MeasurementList = varargin{3};
                NumTracks = size(TrackPredMeasMeans,2);
            end

            [~, NumMeasurements] = size(MeasurementList);

            % Pre-allocate memory
            GateVolumes = zeros(1,NumTracks);
            ValidationMatrix = zeros(NumTracks,NumMeasurements);

            % Construct Validation and GateVolumes Matrices
            for trackInd = 1:NumTracks

                if(nargin==3)
                    % Store pointer to track
                    Track = TrackList{trackInd};

                    % Extract predicted measurement and innovation covariance from filter
                    PredMeasMean = Track.Filter.PredMeasMean;
                else
                    % Extract predicted measurement and innovation covariance from filter
                    PredMeasMean = TrackPredMeasMeans(:,trackInd);
                end
                if(NumMeasurements>0)
                    intersect_1 = rectint(PredMeasMean', MeasurementList');
                    intersect_1 = intersect_1./(PredMeasMean(3)*PredMeasMean(4));
                    intersect_2 = rectint(MeasurementList', PredMeasMean')';
                    intersect_2 = intersect_2./(MeasurementList(3,:).*MeasurementList(4,:));
                    ValidationMatrix(trackInd,:) = (intersect_1>this.OverlapThresh | intersect_2>this.OverlapThresh);
                end
                GateVolumes(trackInd) = PredMeasMean(3)*PredMeasMean(4);
            end
            this.ValidationMatrix = ValidationMatrix;
            this.GateVolumes = GateVolumes;
        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
%% hexahedral_hexagon_beam
% Below is a demonstration of the features of the |hexahedral_hexagon_beam| function

%%
clear; close all; clc;

%% Syntax
% |[Et,Vt,Ft]=hexahedral_hexagon_beam(r,ne,hz,nz,XYZ_centre);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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.
%
% This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function varargout = process_headpoints_remove( varargin )
% PROCESS_HEADPOINTS_REMOVE: Remove head points from the channel file below a certain level.
% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2015-2018

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Remove head points';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Import', 'Channel file'};
    sProcess.Index       = 63;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/ChannelFile#Automatic_registration';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'raw'};
    sProcess.OutputTypes = {'data', 'raw'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    sProcess.isSeparator = 1;
    % Title
    sProcess.options.title.Comment = ['Remove the head points below a certain Z-thresold.
' ...
                                      ' - Z=0 : Nasion, left ear, right ear
 - Z>0 : Towards the top of the head
 - Z<0 : Towards the neck.

'];
    sProcess.options.title.Type    = 'label';
    % Options: tolerance
    sProcess.options.zlimit.Comment = 'Limit Z-coordinate: ';
    sProcess.options.zlimit.Type    = 'value';
    sProcess.options.zlimit.Value   = {0, 'mm', 0};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Get options
    zValue = sProcess.options.zlimit.Value{1} / 1000;
    % Get all the channel files
    uniqueChan = unique({sInputs.ChannelFile});
    % Loop on all the channel files
    for i = 1:length(uniqueChan)
        % Get first input file for this subject
        strMsg = RemoveHeadpoints(uniqueChan{i}, zValue);
        % Report message
        if ~isempty(strMsg)
            bst_report('Info', sProcess, sInputs, strMsg);
        end
    end
    % Return all the files in input
    OutputFiles = {sInputs.FileName};
end

%% ===== REMOVE HEAD POINTS =====
function strMsg = RemoveHeadpoints(ChannelFile, zLimit)
    % Parse inputs
    if (nargin < 2) || isempty(zLimit)
        zLimit = [];
    end
    % Load channel file
    ChannelMat = in_bst_channel(ChannelFile);
    % Display warning: no head points
    if ~isfield(ChannelMat, 'HeadPoints') || isempty(ChannelMat.HeadPoints) || isempty(ChannelMat.HeadPoints.Label)
        strMsg = 'No head points in the file.';
        return;
    end
    % Initial head points number
    nPoints = length(ChannelMat.HeadPoints.Label);
    % Remove selected head points
    if isempty(zLimit)
        ChannelMat.HeadPoints = [];
        iDelete = 1:nPoints;
    else
        % Get EXTRA points
        iExtra = find(strcmpi(ChannelMat.HeadPoints.Type, 'EXTRA'));
        % Find the points below the z-threshold
        iDelete = find(ChannelMat.HeadPoints.Loc(3,iExtra) <= zLimit);
        % Remove the points
        if (length(iDelete) == nPoints)
            ChannelMat.HeadPoints = [];
        elseif ~isempty(iDelete)
            ChannelMat.HeadPoints.Loc(:,iExtra(iDelete)) = [];
            ChannelMat.HeadPoints.Label(iExtra(iDelete)) = [];
            ChannelMat.HeadPoints.Type(iExtra(iDelete))  = [];
        end
    end
    % Message: head points removed
    strMsg = sprintf('%d head points removed, %d points left.', length(iDelete), nPoints - length(iDelete));
    % History: Reamove all head points
    ChannelMat = bst_history('add', ChannelMat, 'headpoints', strMsg);
    % Save file back
    bst_save(file_fullpath(ChannelFile), ChannelMat, 'v7');
end

author: brainstorm-tools


text:
% SB2_DIAGNOSTIC  Helper function to output diagnostics in SPARSEBAYES
%
% USAGE (1):
%
% OPTIONS = SB2_DIAGNOSTIC(OPTIONS, ACTION)
%
% OUTPUT ARGUMENTS:
%
%	OPTIONS		Returned OPTIONS settings, possibly updated
%				(e.g. with file handle after 'open' action)
%
% INPUT ARGUMENTS:
%
%	OPTIONS		Options settings (SB2_USEROPTIONS)
%
%	ACTION		One of:
%				'OPEN' or 'START'			Begin diagnostics
%				'CLOSE', 'END' or 'FINISH'	End diagnostics
%
%
% USAGE (2):
%
%	SB2_DIAGNOSTIC(OPTIONS, LEVEL, MESSAGE, VARARGIN ...)
%
% OUTPUT ARGUMENTS: None
%
% INPUT ARGUMENTS:
%
%	OPTIONS		Options settings (SB2_USEROPTIONS)
%
%	LEVEL		Importance level of message (0 to 4)
%				NB: LOWER is MORE IMPORTANT
%
%	MESSAGE		Message string
%
%	VARARGIN	Optional arguments to fill in string placeholders
%
% EXAMPLE:
%
%	OPTIONS = SB2_Diagnostic(OPTIONS, 'start');
%	SB2_Diagnostic(OPTIONS,2,'Total of %d likelihood updates\n', 25);
%	OPTIONS = SB2_Diagnostic(OPTIONS, 'end');
%
% NOTES:
%
% This function offers a convenient way to output diagnostic information,
% either to screen or to a file.
%
% It will "filter" the incoming messages and only display/write them if
% their LEVEL is less than or equal to the DIAGNOSTICLEVEL as set by
% SB2_USEROPTIONS.
%
% e.g. after
%
%	OPTIONS	= SB2_UserOptions('diagnosticLevel',2)
%
% This message won't display:
%
%	SB2_Diagnostic(OPTIONS, 3, 'Pointless message')
%
% But this will:
%
%	SB2_Diagnostic(OPTIONS, 1, 'Really important message')
%
% Messages of level 0 are therefore most "important", and those of level
% 4 least important. Level-0 messages are always displayed, since the
% DIAGNOSTICLEVEL will always be at least zero.
%

%
% Copyright 2009, Vector Anomaly Ltd
%
% This file is part of the SPARSEBAYES library for Matlab (V2.0).
%
% SPARSEBAYES 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 2 of the License, or (at your option)
% any later version.
%
% SPARSEBAYES 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 SPARSEBAYES in the accompanying file "licence.txt"; if not, write to
% the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
% MA 02110-1301 USA
%
% Contact the author: m a i l [at] m i k e t i p p i n g . c o m
%
function OPTIONS = SB2_Diagnostic(OPTIONS, level, message_, varargin)

%
% Check if "level" is actually a string (this is a control call)
%
if ischar(level)
  switch upper(level)
   case {'OPEN','START'},
    if OPTIONS.diagnosticFID~=1
	  %
	  % If diagnosticFID is not 1, we must have previously set a
      % diagnostics file via SB2_USEROPTIONS
	  %
      OPTIONS.diagnosticFID	= fopen(OPTIONS.diagnosticFile_,'w');
      % Check for file open error
      if OPTIONS.diagnosticFID==-1
		error('Could not open diagnostic file ''%s''', value);
      end
    end
    %
   case {'CLOSE','END','FINISH'},
    if OPTIONS.diagnosticFID~=1
	  % diagnosticFID refers to a file opened above
      fclose(OPTIONS.diagnosticFID);
    end
  end
else
  %
  % Its a message call.
  %
  % Only output diagnostic messages if the "diagnostic level",
  % set via SB2_USEROPTIONS, is equal to or greater than the level of
  % the message to be displayed.
  %
  if level<=OPTIONS.diagnosticLevel
    fprintf(OPTIONS.diagnosticFID, repmat('    ', [1 level]));
    fprintf(OPTIONS.diagnosticFID, message_, varargin{:});
  end
end


text:
%% Default template
%% Visualize the Data

plot(ebsd)

%% Calculate an ODF

odf = calcODF(ebsd)

%% Detect grains

%segmentation angle
segAngle = 10*degree;

grains = calcGrains(ebsd,'threshold',segAngle);

%% Orientation of Grains

plot(grains)

author: mtex-toolbox

repo: mtex

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


text:
function newFaces = minConvexHull(points, varargin)
%MINCONVEXHULL Return the unique minimal convex hull of a set of 3D points
%
%   FACES = minConvexHull(PTS)
%   NODES is a set of 3D points  (as a Nx3 array). The function computes
%   the convex hull, and merge contiguous coplanar faces. The result is a
%   set of polygonal faces, such that there are no coplanar faces.
%   FACES is a cell array, each cell containing the vector of indices of
%   nodes given in NODES for the corresponding face.
%
%   FACES = minConvexHull(PTS, PRECISION)
%   Adjust the threshold for deciding if two faces are coplanar or
%   parallel. Default value is 1e-14.
%
%   Example
%     % extract square faces from a cube
%     [n, e, f] = createCube;
%     f2 = minConvexHull(n);
%     drawMesh(n, f2);
%
%     % Subdivides and smooths a mesh rpresenting a cube
%     [n, e, f] = createCube;
%     [n2, f2] = subdivideMesh(n, triangulateFaces(f), 4);
%     [n3, f3] = smoothMesh(n2, f2);
%     figure; drawMesh(n3, f3);
%     axis equal; view(3);
%     % merge coplanar faces, making apparent the faces of the original cube
%     f4 = minConvexHull(n3);
%     figure; drawMesh(n3, f4);
%     axis equal; view(3);
%
%
%   See also
%   meshes3d, mergeCoplanarFaces, drawMesh, convhull, convhulln
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2006-07-05
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

% HISTORY
%   20/07/2006 add tolerance for coplanarity test
%   21/08/2006 fix small bug due to difference of methods to test
%       coplanarity, sometimes resulting in 3 points of a face being not
%       coplanar! Also add control on precision
%   18/09/2007 ensure faces are given as horizontal vectors

% set up precision
acc = 1e-14;
if ~isempty(varargin)
    acc = varargin{1};
end

% triangulated convex hull. It is not uniquely defined.
faces = convhulln(points);

% compute centroid of the nodes
pointsCentroid = centroid(points);

% number of base triangular faces
N = size(faces, 1);

% compute normals of given faces
normals = planeNormal(createPlane(...
    points(faces(:,1),:), points(faces(:,2),:), points(faces(:,3),:)));

% initialize empty faces
newFaces = {};

% Processing flag for each triangle
% 1 : triangle to process, 0 : already processed
% in the beginning, every triangle face need to be processed
flag = ones(N, 1);

% iterate on each triangular face of the convex hull
for iFace = 1:N

    % check if face was already performed
    if ~flag(iFace)
        continue;
    end

    % indices of faces with same normal
    ind = find(abs(vectorNorm3d(cross(repmat(normals(iFace, :), [N 1]), normals)))<acc);
    ind = ind(ind~=iFace);

    % keep only coplanar faces (test coplanarity of points in both face)
    ind2 = iFace;
    for j = 1:length(ind)
        if isCoplanar(points([faces(iFace,:) faces(ind(j),:)], :), acc)
            ind2 = [ind2 ind(j)]; %#ok<AGROW>
        end
    end

    % compute order of the vertices in current face
    faceVertices = unique(faces(ind2, :));
    [tmp, I]  = angleSort3d(points(faceVertices, :)); %#ok<ASGLU>

    % create the new face, ensuring it is a row vector
    face = faceVertices(I);
    face = face(:)';

    % ensure face has normal pointing outwards
    outerNormal = meshFaceCentroids(points, face) - pointsCentroid;
    if dot(meshFaceNormals(points, face), outerNormal, 2) < 0
        face = face([1 end:-1:2]);
    end

    % add a new face to the list
    newFaces = [newFaces {face}]; %#ok<AGROW>

    % mark processed faces
    flag(ind2) = 0;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function SDP = chordal_relax_category_registration_v2(problem,varargin)
%% Apply a sparse, chordal, SECOND-order relaxation to category registration
%% Depending on multivariate polynomial package in SPOT
%% residual of v1: b(i) - R * (sum_k c_k ak(i)) - t
%% residual of v2: R*b(i) + t - sum_k c_k ak(i)
%% Heng Yang
%% July 06, 2021

params = inputParser;
params.CaseSensitive = false;

params.addParameter('checkMonomials',true, @(x) islogical(x));
params.addParameter('lambda',0.1, @(x) isscalar(x));

params.parse(varargin{:});

checkMonomials = params.Results.checkMonomials;
lambda         = params.Results.lambda;

fprintf('\n===================================================================')
fprintf('\nApplying Chordal SDP relaxation to category registration problem')
fprintf('\n===================================================================\n')
t0              = tic;

N               = problem.N;
K               = problem.K;
scene           = problem.scene;
shapes          = problem.shapes;
noiseBoundSq    = problem.noiseBoundSq;
tBound          = problem.translationBound;
tBoundSq        = tBound^2; % t'*t <= tBoundSq
cBoundSq        = problem.cBound^2; % should just be 1
barc2           = 1.0;

%% define POP variables
nrPrimalVars    = 9+3+K+N; % rotation: 9, translation: 3, binary: N, shape: K
p               = msspoly('p',nrPrimalVars);
r               = p(1:9);
R               = reshape(r,3,3);
col1 = R(:,1); col2 = R(:,2); col3 = R(:,3);
t               = p(10:12);
c               = p(12+1:12+K);
theta           = p(12+K+1:nrPrimalVars);
x               = [r;t;c];

%% define cost function
shape           = combine_shapes(shapes,c);
residuals       = [];
for i = 1:N
    distance            = R*scene(:,i) + t - shape(:,i);
    residuals           = [residuals; (distance' * distance) / noiseBoundSq];
end
f_cost = lambda * (c'*c); % regularization term
for i = 1:N
    f_cost = [f_cost;(1+theta(i))/2 * residuals(i) + (1-theta(i))/2 * barc2];
end

%% define constraints
h_r  = [1.0-col1'*col1;...
        1.0-col2'*col2;...
        1.0-col3'*col3;... % column unit length
        col1'*col2;...
        col2'*col3;...
        col3'*col1;... % colums orthogonal
        cross(col1,col2) - col3;...
        cross(col2,col3) - col1;...
        cross(col3,col1) - col2]; % columns righthandedness

h_theta = [];
for i = 1:N
    h_theta =[h_theta; 1-theta(i)^2];
end

g_t = tBoundSq - t'*t; % Translation bounded
g_c = [cBoundSq - c'*c;c]; % nonnegative and bounded shape parameters
g_x = [g_t;g_c];

%% Formulate the chordal sparse second-order relaxation
%% the 0-th block [1;x] * [1;x]'
basis0          = [1;x];
n0              = length(basis0);
basis_x0        = get_multiplier_basis(x,basis0,h_r(1));
pop0            = [mykron(basis_x0,h_r);...
                   mykron(basis0,basis0);...
                   f_cost(1)];
[~,degmat,coef_all] = decomp(pop0);
coef_all            = coef_all';
dim_loc0        = length(basis_x0) * length(h_r);
n0delta         = triangle_number(n0);
nterms          = size(degmat,1);
m_mom0          = n0delta - nterms;

assert(m_mom0==0,'The zero-th blk should have 0 moment constraints.')

coef_mom    = coef_all(:,dim_loc0+1:dim_loc0+n0^2);
coef_mom    = coef_mom';
B           = {};
B_normalize = {};

for i = 1:nterms
    [row,~,~]   = find(coef_mom(:,i));
    SDP_coli    = floor((row-1)./n0) + 1;
    SDP_rowi    = mod(row-1,n0) + 1;
    nnz         = length(SDP_rowi);

    Bi          = sparse(SDP_rowi,SDP_coli,ones(nnz,1),n0,n0);
    B{end+1}    = Bi;
    B_normalize{end+1} = Bi/nnz;
end

coef_loc0       = coef_all(:,1:dim_loc0);
A0_local        = {};

for i = 1:dim_loc0
    [rowi,~,vi] = find(coef_loc0(:,i));
    Ai      = sparse(n0,n0);
    for j   = 1:length(rowi)
        Ai  = Ai + vi(j) * B_normalize{rowi(j)};
    end
    A0_local = [A0_local;{Ai}];
end
A0_0    = sparse([1],[1],[1],n0,n0);
% The first block satisfies A0(X0) = b0;
A0      = [{A0_0};A0_local];
b0      = sparse(1,1,1,length(A0),1);

% Now build the cost matrix
coef_cost   = coef_all(:,dim_loc0+n0^2+1);
[row,~,v]   = find(coef_cost);
C           = sparse(n0,n0);
for i = 1:length(row)
    C       = C + v(i) * B_normalize{row(i)};
end

%% the 1-N blocks [1;x;theta(i);theta(i)*x] * [1;x;theta(i);theta(i)*x]'
%% Since there are (2+K) inequality constraints, it will generate (K+3)*N blocks
nrineq          = length(g_x);
Aall            = {};
Asuball         = {};
ball            = [];
Call            = {C};
A0append        = {};
for blkidx = 1:N
    basis       = [1;x;theta(blkidx);theta(blkidx)*x];
    n           = length(basis);
    basis_x     = get_multiplier_basis([x;theta(blkidx)],basis,h_r(1),0);
    basis_theta = get_multiplier_basis([x;theta(blkidx)],basis,h_theta(blkidx),0);
    basis_g     = [1;theta(blkidx)];

    out         = gen_chordal_subblk_catreg_v2(...
        basis,basis_x,h_r,basis_theta,h_theta(blkidx),g_x,basis_g,f_cost(1+blkidx));

    Acell       = out.A;
    Asub        = out.Asub;
    b           = out.b;
    C           = out.C;

    % add constraint that the top-left [1;x]*[1;x]' block is the same as
    % the 0-th block
    A0blk = {};
    for i = 1:n0
        for j = i:n0
            if i == j
                A0i  = sparse(i,j,-1,n0,n0);
                Ai   = sparse(i,j,1,n,n);
            else
                A0i  = sparse([i,j],[j,i],[-0.5,-0.5],n0,n0);
                Ai   = sparse([i,j],[j,i],[0.5,0.5],n,n);
            end
            A0blk    = [A0blk;{A0i}];
            Acell    = [Acell;{Ai}];
        end
    end

    ball             = [ball;b;sparse(n0delta,1)];
    A0append{end+1}  = A0blk;
    Aall{end+1}      = Acell;
    Call             = [Call;C];
    Asuball{end+1}   = Asub;
end

%% Convert to standard SDPT3 format
b           = [b0;ball];
blk         = cell( (nrineq+1)*N+1,2);
blk{1,1}    = 's'; blk{1,2} = n0;
n1          = out.blk{2,2};
n1delta     = triangle_number(n1);
step        = nrineq + 1;
for i = 1:N
    blk{step*i-nrineq+1,1} = 's';
    blk{step*i-nrineq+1,2} = n;
    for j = 1:nrineq
        blk{step*i-nrineq+j+1,1} = 's';
        blk{step*i-nrineq+j+1,2} = n1;
    end
end

A0t     = sparsesvec(blk(1,:),A0);
for i = 1:N
    A0t = [A0t,...
           sparse(n0delta,out.m),...
           sparsesvec(blk(1,:),A0append{i})];
end

ndelta  = triangle_number(n);
At      = {A0t};
for i = 1:N
    Ait = [sparse(ndelta,length(b0)),...
           sparse(ndelta,(i-1)*length(Aall{i})),...
           sparsesvec(blk(step*i-nrineq+1,:),Aall{i}),... % the moment matrix block
           sparse(ndelta,(N-i)*length(Aall{i}))];
    At  = [At;{Ait}];
    for j = 1:nrineq
        Aijt    = [sparse(n1delta,length(b0)),...
                    sparse(n1delta,(i-1)*length(Aall{i})),...
                    sparse(n1delta,out.m_mom+out.m_loc),... % the moment constraint
                    sparse(n1delta,(j-1)*n1delta),...
                    sparsesvec(blk(step*i-nrineq+j+1,:),Asuball{i}{j}),...
                    sparse(n1delta,(nrineq-j)*n1delta),...
                    sparse(n1delta,n0delta),...
                    sparse(n1delta,(N-i)*length(Aall{i}))];

        At  = [At;{Aijt}];
    end
end

SDP.blk = blk;
SDP.At  = At;
SDP.m   = length(b);
SDP.C   = Call;
SDP.b   = b;

tf    = toc(t0);
fprintf('\nDone in %g seconds.\n',tf);
fprintf('===================================================================\n')

%% Convert to Sedumi format
fprintf('Convert to Sedumi format ...\n')
t0    = tic;
sK.s  = [n0];
for i = 1:N
    sK.s        = [sK.s,n,n1*ones(1,nrineq)];
end

A0t     = sparsevec(blk(1,:),A0);
n0sq    = n0^2;
for i = 1:N
    A0t = [A0t,...
           sparse(n0sq,out.m),...
           sparsevec(blk(1,:),A0append{i})];
end

nsq     = n^2;
n1sq    = n1^2;
At      = {A0t};
for i = 1:N
    Ait = [sparse(nsq,length(b0)),...
           sparse(nsq,(i-1)*length(Aall{i})),...
           sparsevec(blk(step*i-nrineq+1,:),Aall{i}),... % the moment matrix block
           sparse(nsq,(N-i)*length(Aall{i}))];
    At  = [At;{Ait}];
    for j = 1:nrineq
        Aijt    = [sparse(n1sq,length(b0)),...
                    sparse(n1sq,(i-1)*length(Aall{i})),...
                    sparse(n1sq,out.m_mom+out.m_loc),... % the moment constraint
                    sparse(n1sq,(j-1)*n1delta),...
                    sparsevec(blk(step*i-nrineq+j+1,:),Asuball{i}{j}),...
                    sparse(n1sq,(nrineq-j)*n1delta),...
                    sparse(n1sq,n0delta),...
                    sparse(n1sq,(N-i)*length(Aall{i}))];

        At  = [At;{Aijt}];
    end
end

sdata.K     = sK;
sdata.At    = cat(1,At{:});
sdata.b     = b;

sc          = [];
for i = 1:length(Call)
    sc      = [sc;sparsevec(blk(i,:),Call(i))];
end
sdata.c     = sc;

SDP.sedumi   = sdata;

tf    = toc(t0);
fprintf('Done in %g seconds.\n',tf);
fprintf('===================================================================\n')

end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function pass = test_minandmax3est(pref)
% Test MINANDMAX3EST().

if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 1000*pref.cheb3Prefs.chebfun3eps;

f = chebfun3v();
pass(1) = isempty(minandmax3est(f));

% CHEBFUN3V (1 component)
f = chebfun3v(@(x,y,z) x, [ -2, 4, -1, 1, -1, 1 ]);
mM = minandmax3est(f);
pass(2) = ( length(mM) == 2 );
pass(3) = ( norm(mM - [ -2, 4 ]) < tol );

% CHEBFUN3V (3 components)
f = chebfun3v(@(x,y,z) x, @(x,y,z) y, @(x,y,z) z, [ -2, 4, 3, 17, -1, 42 ]);
mM = minandmax3est(f);
pass(4) = ( length(mM) == 6 );
pass(5) = ( norm(mM - [ -2, 4, 3, 17, -1, 42 ]) < tol );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function paramV = mtrFilename2Param(fileName,paramName)

% Will retrieve paramnames from files of the form *paramName_#_*.dat

sI = strfind(fileName,paramName)+length(paramName)+1;
eIVec = strfind(fileName,'_');
eI = min( eIVec(eIVec>sI) )-1;
if isempty(eI)
    eI = max(strfind(fileName,'.'))-1;
end
paramV = str2num(fileName(sI:eI));
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function extent = imPhysicalExtent(img, varargin)
%IMPHYSICALEXTENT Compute the physical extent of an image
%
%   BOX = imPhysicalExtent(IMG)
%   Returns the physical extent of the image IMG, such that the display
%   will be contained within the extent. IMG can be either 2D or 3D image.
%
%   BOX = imPhysicalExtent(IMG, SPACING)
%   Computes the extent by taking into account the resolution of the image.
%   The resolution is given in X-Y-Z order (ie., SPACING(1) corresponds to
%   the resolution of the second dimension of IMG).
%
%   BOX = imPhysicalExtent(IMG, SPACING, ORIGIN)
%   Also takes into account the stack origin, i.e. the coordinates of the
%   first voxel expressed in user coordinates.
%   The origin is given in X-Y-Z order.
%
%   BOX = imPhysicalExtent(IMG, 'spacing', SPACING, 'origin', ORIGIN)
%   Uses a syntax based on parameter name-value pairs.
%
%   BOX = imPhysicalExtent(DIM, SPACING, ORIGIN)
%   Computes the physical extent based on the size of the stack. This
%   syntax can be used to avoid passing an array as parameter. DIM is given
%   in X-Y-Z order.
%
%
%   Example
%   % Compute physical extent of MRI Human head
%     metadata = analyze75info('brainMRI.hdr');
%     I = analyze75read(metadata);
%     box = imPhysicalExtent(I, [1 1 2.5])
%     box =
%         0.5000  128.5000    0.5000  128.5000    -0.2500   67.2500
%
%   See also
%     imSize
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-03-03,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% size of image in each physical direction
sz = imSize(img);
% if length(dim) > 2
%     % compute size of the stack
%     sz = stackSize(img);
% else
%     % the size of the stack is given as input
%     sz = img;
% end

% default spacing and origin (matlab display convention)
nd = length(sz);
sp = ones(1, nd);
or = ones(1, nd);

% parse origin and spacing of the stack
if ~isempty(varargin)
    var = varargin{1};

    if isnumeric(var)
        % extract voxel spacing
        sp = var;

        if length(varargin) > 1
            % also extract voxel origin
            or = varargin{2};
        end

    elseif ischar(var)
        while length(varargin) > 2
            paramName = varargin{1};
            if strcmp(paramName, 'spacing')
                sp = varargin{2};

            elseif strcmp(paramName, 'origin')
                or = varargin{2};

            else
                error(['Unknown parameter: ' paramName]);
            end
        end
    end
end

% put extent in array
extent = (([zeros(nd, 1) sz'] - .5).* [sp' sp'] + [or' or'])';

% change array shape to get a single row
extent = extent(:)';

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [on,rt,ac] = spm_ADEM_cue_rt(DEM)
% returns reaction times and accuracy for ADEM_cued_response demo
% FORMAT [on,rt,ac] = spm_ADEM_cue_rt(DEM)
%
% DEM - DEM structure from ADEM_cued_response.m
%
% on  - cue onset
% ac  - accuracy
% rt  - reaction time
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_ADEM_set_rt.m 4231 2011-03-07 21:00:02Z karl $

% distance from target and cue contrast
%--------------------------------------------------------------------------
P   = DEM.pP.P{1};                            % location of targets
L   = DEM.pU.v{1}((1:2) + 2,:);               % location of finger
n   = length(P);                              % number of targets
N   = length(L);                              % number of time bins
D   = sparse(n,N);
for j = 1:N
    for i = 1:n
        D(i,j) = D(i,j) + sum((L(:,j) - P(:,i)).^2);
    end
end
D   = sqrt(D);                                % distance from targets
C   = DEM.pU.v{1}((1:n) + 4,:);               % contrast of targets

r   = 1/32;                                   % radius of proximity
c   = diff(C > 1/2,1,2) > 0;                  % target onset
on  = {};                                     % cue onset
ac  = {};                                     % accuracy
rt  = {};                                     % reaction time

% get performance
%--------------------------------------------------------------------------
for i = 1:n

    on{i} = find(c(i,:));
    for j = 1:length(on{i})
        try

            % minimum distance
            %--------------------------------------------------------------
            d        = D(i,(1:8) + on{i}(j))';
            ac{i}(j) = sqrt(min(d));

            % estimated reaction time
            %--------------------------------------------------------------
            X        = (1:length(d))' - 1;
            B        = pinv([X.^0 X.^1])*log(d);
            rt{i}(j) = (log(r) - B(1))/B(2);

        catch
            ac{i}(j) = NaN;
            rt{i}(j) = NaN;
        end
    end
end

% sort trials (over all targets
%--------------------------------------------------------------------------
on = spm_vec(on); [i j] = sort(on,1,'ascend'); on = on(j);
ac = spm_vec(ac); ac = ac(j);
rt = spm_vec(rt); rt = rt(j);

% remove first trial
%--------------------------------------------------------------------------
on(1) = [];
rt(1) = [];
ac(1) = [];

% convert spatial error to accuracy
%--------------------------------------------------------------------------
ac    = 1./ac;

% convert reaction time to ms
%--------------------------------------------------------------------------
dt    = 64/1000;
rt    = rt*1000*dt;
on    = on*dt;

author: spm


text:
%CA driver
%HPP-gas

clear all
clf

nx=52; %must be divisible by 4
ny=100;

z=zeros(nx,ny);
o=ones(nx,ny);
sand = z ;
sandNew = z;
gnd = z ;
diag1 = z;
diag2 = z;
and12 = z;
or12 = z;
sums = z;
orsum = z;

gnd(1:nx,ny-3)=1 ; % right ground line
gnd(1:nx,3)=1 ; % left ground line
gnd(nx/4:nx/2-2,ny/2)=1; %the hole line
gnd(nx/2+2:nx,ny/2)=1; %the hole line
gnd(nx/4, 1:ny) = 1; %top line
gnd(3*nx/4, 1:ny) = 1 ;%bottom line

%fill the left side
r = rand(nx,ny);
sand(nx/4+1:3*nx/4-1, 4:ny/2-1) = r(nx/4+1:3*nx/4-1, 4:ny/2-1)<0.3;
%sand(nx/4+1:3*nx/4-1, ny*.75:ny-4) = r(nx/4+1:3*nx/4-1, ny*.75:ny-4)<0.75;
%sand(nx/2,ny/2) = 1;
%sand(nx/2+1,ny/2+1) = 1;

imh = image(cat(3,z,sand,gnd));
set(imh, 'erasemode', 'none')
axis equal
axis tight

for i=1:1000
    p=mod(i,2); %margolis neighborhood

    %upper left cell update
    xind = [1+p:2:nx-2+p];
    yind = [1+p:2:ny-2+p];

    %See if exactly one diagonal is ones
    %only (at most) one of the following can be true!
    diag1(xind,yind) = (sand(xind,yind)==1) & (sand(xind+1,yind+1)==1) & ...
        (sand(xind+1,yind)==0) & (sand(xind,yind+1)==0);

    diag2(xind,yind) = (sand(xind+1,yind)==1) & (sand(xind,yind+1)==1) & ...
        (sand(xind,yind)==0) & (sand(xind+1,yind+1)==0);

    %The diagonals both not occupied by two particles
    and12(xind,yind) = (diag1(xind,yind)==0) & (diag2(xind,yind)==0);

    %One diagonal is occupied by two particles
    or12(xind,yind)  = diag1(xind,yind) | diag2(xind,yind);

    %for every gas particle see if it near the boundary
    sums(xind,yind) = gnd(xind,yind) | gnd(xind+1,yind) | ...
                        gnd(xind,yind+1) | gnd(xind+1,yind+1) ;

    % cell layout:
    % x,y    x+1,y
    % x,y+1  x+1,y+1
    %If (no walls) and (diagonals are both not occupied)
    %then there is no collision, so move opposite cell to current cell
    %If (no walls) and (only one diagonal is occupied)
    %then there is a collision so move ccw cell to the current cell
    %If (a wall)
    %then don't change the cell (causes a reflection)
    sandNew(xind,yind) = ...
        (and12(xind,yind)  & ~sums(xind,yind) & sand(xind+1,yind+1)) + ...
        (or12(xind,yind) & ~sums(xind,yind) & sand(xind,yind+1)) + ...
        (sums(xind,yind) & sand(xind,yind));

    sandNew(xind+1,yind) = ...
        (and12(xind,yind)  & ~sums(xind,yind) & sand(xind,yind+1)) + ...
        (or12(xind,yind) & ~sums(xind,yind) & sand(xind,yind))+ ...
        (sums(xind,yind) & sand(xind+1,yind));

    sandNew(xind,yind+1) = ...
        (and12(xind,yind)  & ~sums(xind,yind) & sand(xind+1,yind)) + ...
        (or12(xind,yind) & ~sums(xind,yind) & sand(xind+1,yind+1))+ ...
        (sums(xind,yind) & sand(xind,yind+1));

     sandNew(xind+1,yind+1) = ...
        (and12(xind,yind)  & ~sums(xind,yind) & sand(xind,yind)) + ...
        (or12(xind,yind) & ~sums(xind,yind) & sand(xind+1,yind))+ ...
        (sums(xind,yind) & sand(xind+1,yind+1));

    sand = sandNew;

    set(imh, 'cdata', cat(3,z,sand,gnd) )
    drawnow
end

author: Eurus-Holmes


text:
function r = rdivide(a,b)
%RDIVIDE      Hessian elementwise right division  a ./ b
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

  if prod(size(b))==1               % scalar denominator
    r = a / b;
    return
  end

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  scalar_a = ( prod(size(a))==1 );
  if ( ~scalar_a ) & ( ~isequal(size(a),size(b)) )
    error('dimensions do not match in hessian elementwise division')
  end

  N = getappdata(0,'INTLAB_HESSIAN_NUMVAR');

  % check for emptyness: cures Matlab bug
  % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
  % throughout this routine use e.g.  sparse(ib,jb,reshape(D(jb),size(sb)).*sb,N,m);
  % rather than  sparse(ib,jb,D(jb).*sb(:),N,m);

  if ~isa(a,'hessian')              % non-hessian (scalar or array) ./ hessian array
    if issparse(b.hx)
      m = prod(size(b.x));
      r.x = a ./ b.x;
      D = -r.x ./ b.x;
      [ib,jb,sb] = find(b.dx);
      r.dx = sparse(ib,jb,reshape(D(jb),size(sb)).*sb,N,m);
      bdx = sparse(ib,jb,reshape(1./b.x(jb),size(sb)).*sb,N,m);
      r.hx = adx2rhx(N,m,bdx,b.dx);
      [ib,jb,sb] = find(b.hx-r.hx);
      r.hx = sparse(ib,jb,reshape(D(jb),size(sb)).*sb,N^2,m);
    else
      r.x = a ./ b.x;
      D = r.x ./ b.x;
      D = repmat(-D(:).',N^2,1);
      r.dx = D(1:N,:) .* b.dx;
      bdx1 = b.dx./repmat(b.x(:).',N,1);
      index = repmat(1:N,N,1);
      r.hx = D .* ( b.hx - bdx1(index,:).*b.dx(index',:) );
    end
  elseif ~isa(b,'hessian')          % hessian ./ non-hessian array
    r.x = a.x ./ b;
    bb = 1 ./ (b(:).');
    if scalar_a                     % hessian scalar ./ non-hessian array
      r.dx = a.dx * bb;
      r.hx = a.hx * bb;
    else                            % hessian array ./ non-hessian array
      if issparse(a.hx)
        m = prod(size(a.x));
        [ia,ja,sa] = find(a.dx);
        r.dx = sparse(ia,ja,reshape(bb(ja),size(sa)).*sa,N,m);
        [ia,ja,sa] = find(a.hx);
        r.hx = sparse(ia,ja,reshape(bb(ja),size(sa)).*sa,N^2,m);
      else
        bb = repmat(bb,N^2,1);
        r.dx = a.dx .* bb(1:N,:);
        r.hx = a.hx .* bb;
      end
    end
  else                              % hessian ./ hessian array
    m = prod(size(b.x));
    r.x = a.x ./ b.x;
    if scalar_a                     % hessian scalar ./ hessian array
      if issparse(a.hx) | issparse(b.hx)
        bx = sparse(b.x(:));
        D = 1 ./ sqr(bx);
        Num = a.dx * bx.' - a.x * b.dx;
        [i,j,s] = find(Num);
        r.dx = sparse(i,j,s.*reshape(D(j),size(s)),N,m);
        Num = sparse(i,j,s./reshape(bx(j),size(s)),N,m);
        r.hx = a.hx * sparse(bx.') - adx2rhx(N,m,b.dx,Num) - a.x * b.hx;
        [i,j,s] = find(r.hx);
        r.hx = sparse(i,j,s.*reshape(D(j),size(s)),N^2,m);
      else
        Num = a.dx * b.x(:).' - a.x * b.dx;
        D = repmat( sqr((1./b.x(:)).') , N^2 , 1 );
        r.dx = Num .* D(1:N,:);
        index = repmat(1:N,N,1);
        Num = Num ./ repmat(b.x(:).',N,1);
        r.hx = ( a.hx * b.x(:).' - reshape(b.dx(index,:) .* Num(index',:),size(b.hx)) - a.x * b.hx ) .* D;
      end
    else                            % hessian array ./ hessian array
      if issparse(a.hx) | issparse(b.hx)
        [ia,ja,sa] = find(a.dx);
        [ib,jb,sb] = find(b.dx);
        ax = a.x(:);
        bx = b.x(:);
        Num = sparse(ia,ja,reshape(bx(ja),size(sa)).*sa,N,m) - sparse(ib,jb,reshape(ax(jb),size(sb)).*sb,N,m);
        D = sqr( 1 ./ bx );
        [ia,ja,sa] = find(Num);
        r.dx = sparse(ia,ja,reshape(D(ja),size(sa)).*sa,N,m);
        Num = sparse(ia,ja,sa./reshape(bx(ja),size(sa)),N,m);
        [ia,ja,sa] = find(a.hx);
        [ib,jb,sb] = find(b.hx);
        r.hx = sparse(ia,ja,reshape(bx(ja),size(sa)).*sa,N^2,m) - adx2rhx(N,m,b.dx,Num) - ...
               sparse(ib,jb,reshape(ax(jb),size(sb)).*sb,N^2,m);
        [ia,ja,sa] = find(r.hx);
        r.hx = sparse(ia,ja,sa.*reshape(D(ja),size(sa)),N^2,m);
      else
        ax = repmat(a.x(:).',N^2,1);
        bx = repmat(b.x(:).',N^2,1);
        Num = a.dx .* bx(1:N,:) - ax(1:N,:) .* b.dx;
        D = repmat( sqr((1./b.x(:)).') , N^2 , 1 );
        r.dx = Num .* D(1:N,:);
        index = repmat(1:N,N,1);
        Num = Num ./ bx(1:N,:);
        r.hx = ( a.hx .* bx - reshape(b.dx(index,:) .* Num(index',:),size(a.hx)) - ax .* b.hx ) .* D;
      end
    end
  end

  r = class(r,'hessian');

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
classdef LossTopBoxSmoothProb < dagnn.Loss
  % given top scoring box, it finds other boxes with at least overlap of
  % minOverlap and calculates the euclidean dist between top and other
  % boxes

  properties (Transient)
    gtIdx = []
    boxIdx = []
    probs = []
    minOverlap = 0.5
    nBoxes = 10
  end

  methods
    function outputs = forward(obj, inputs, params)
      if numel(inputs) ~= 4
        error('Number of inputs is not 2');
      end
      obj.gtIdx = [];
      obj.boxIdx = [];
      obj.probs = [];
      boxes  = double(gather(inputs{2})');
      scores = gather(squeeze(inputs{3}));
      labels = gather(squeeze(inputs{4}));

      if numel(boxes)<5
        return;
      end

      outputs{1} = zeros(1,'like',inputs{1});
      for c=1:numel(labels)
        if labels(c)<=0
          continue;
        end

        [so, si] = sort(scores(c,:),'descend');
        obj.gtIdx{c} = si(1);
        gtBox = boxes(:,obj.gtIdx{c});
        gtArea = (gtBox(3)-gtBox(1)+1) .* (gtBox(4)-gtBox(2)+1);

        bbs = boxes(:,si(2:min(obj.nBoxes,end)))';

        y1 = bbs(:,1);
        x1 = bbs(:,2);
        y2 = bbs(:,3);
        x2 = bbs(:,4);

        area = (x2-x1+1) .* (y2-y1+1);

        yy1 = max(gtBox(1), y1);
        xx1 = max(gtBox(2), x1);
        yy2 = min(gtBox(3), y2);
        xx2 = min(gtBox(4), x2);

        w = max(0.0, xx2-xx1+1);
        h = max(0.0, yy2-yy1+1);

        inter = w.*h;
        o = find((inter ./ (gtArea + area - inter))>obj.minOverlap);

        if isempty(o)
          continue;
        end

        obj.boxIdx{c} = si(o+1);
        obj.probs{c} = so(o+1);
        d = bsxfun(@minus,inputs{1}(:,:,:,obj.boxIdx{c}),inputs{1}(:,:,:,obj.gtIdx{c}));
        d = bsxfun(@times,d,obj.probs{c});
        outputs{1} = outputs{1} + 0.5 * sum(d(:).^2);
      end

      n = obj.numAveraged ;
      m = n + 1 ;
      obj.average = (n * obj.average + gather(outputs{1})) / m ;
      obj.numAveraged = m ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      derInputs = cell(1,4) ;
      derInputs{1} = zeros(size(inputs{1}),'like',inputs{1});
      for c=1:numel(obj.boxIdx)
        if isempty(obj.boxIdx{c}), continue; end
        derInputs{1}(:,:,:,obj.boxIdx{c}) = ...
          bsxfun(@minus,inputs{1}(:,:,:,obj.boxIdx{c}),inputs{1}(:,:,:,obj.gtIdx{c}));
        derInputs{1}(:,:,:,obj.boxIdx{c}) = bsxfun(@times,...
          reshape(obj.probs{c},[1 1 1 numel(obj.probs{c})]),derInputs{1}(:,:,:,obj.boxIdx{c}));
        derInputs{1}(:,:,:,obj.gtIdx{c}) = -sum(derInputs{1}(:,:,:,obj.boxIdx{c}),4);

      end
      derInputs{1} = derInputs{1} * derOutputs{1};
%       fprintf('LossTopBox l2 %f ',sqrt(sum(derInputs{1}(:).^2)));
      derParams = {} ;
    end

    function obj = LossTopBoxSmoothProb(varargin)
      obj.load(varargin) ;
      obj.loss = 'LossTopBoxSmoothProb';
    end

    function reset(obj)
      obj.gtIdx = [];
      obj.boxIdx = [];
      obj.probs = [];
      obj.average = 0 ;
      obj.numAveraged = 0 ;
    end

  end

end

author: hbilen

repo: WSDDN

save_path: github-repos/MATLAB/hbilen-WSDDN


text:
function [X,Y,Z] = idp3elli(idp,IDP,ns,NP)

%IDP3ELLI 3D ellipsoid from Gaussian IDP.
%   [X,Y,Z] = IDP3ELLI(x0,P,ns,NP) gives X, Y and Z coordinates of the
%   points corresponding to the 2 biggest semi-diametres of the ellipsoid
%   defined by the covariances matrix P and centered at x0:
%
%        (x-x0)'*(P^-1)*(x-x0) = ns^2.
%
%   where P and x0 are obtained by transforming the input inverse-depth
%   point idp and covariance IDP to an euclidean point:
%
%       x0 = idp2p(idp)
%       P = J*IDP*J'
%
%   being J the Jacobian of the conversion function. This conversion is
%   performed internally by the function PROPAGATEUNCERTAINTY.
%
%   The ellipsoid can be plotted in a 3D graphic by just creating a line
%   with line(X,Y,Z).
%
%   See also COV3ELLI, LINE, IDP2P, PROPAGATEUNCERTAINTY.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

[p,P] = propagateUncertainty(idp,IDP,@idp2p);

[X,Y,Z] = cov3elli(p,P,ns,NP);

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function graph_to_dot(adj, varargin)
%GRAPH_TO_DOT  Makes a GraphViz (AT&T) file representing an adjacency matrix
% graph_to_dot(adj, ...) writes to the specified filename.
%
% Optional arguments can be passed as name/value pairs: [default]
%
% 'filename' - if omitted, writes to 'tmp.dot'
% 'arc_label' - arc_label{i,j} is a string attached to the i-j arc [""]
% 'node_label' - node_label{i} is a string attached to the node i ["i"]
% 'width'     - width in inches [10]
% 'height'    - height in inches [10]
% 'leftright' - 1 means layout left-to-right, 0 means top-to-bottom [0]
% 'directed'  - 1 means use directed arcs, 0 means undirected [1]
%
% For details on graphviz, See http://www.research.att.com/sw/tools/graphviz
%
% See also dot_to_graph and draw_dot.

% First version written by Kevin Murphy 2002.
% Modified by Leon Peshkin, Jan 2004.
% Bugfix by Tom Minka, Mar 2004.

node_label = [];   arc_label = [];   % set default args
width = 10;        height = 10;
leftright = 0;     directed = 1;     filename = 'tmp.dot';

for i = 1:2:nargin-1                    % get optional args
  switch varargin{i}
    case 'filename', filename = varargin{i+1};
    case 'node_label', node_label = varargin{i+1};
    case 'arc_label', arc_label = varargin{i+1};
    case 'width', width = varargin{i+1};
    case 'height', height = varargin{i+1};
    case 'leftright', leftright = varargin{i+1};
    case 'directed', directed = varargin{i+1};
  end
end
% minka
if ~directed
  adj = triu(adj | adj');
end

fid = fopen(filename, 'w');
if directed
  fprintf(fid, 'digraph G {\n');
  arctxt = '->';
  if isempty(arc_label)
    labeltxt = '';
  else
    labeltxt = '[label="%s"]';
  end
else
  fprintf(fid, 'graph G {\n');
  arctxt = '--';
  if isempty(arc_label)
    labeltxt = '[dir=none]';
  else
    labeltext = '[label="%s",dir=none]';
  end
end
edgeformat = strcat(['%d ',arctxt,' %d ',labeltxt,';\n']);
fprintf(fid, 'center = 1;\n');
fprintf(fid, 'size=\"%d,%d\";\n', width, height);
if leftright
  fprintf(fid, 'rankdir=LR;\n');
end
Nnds = length(adj);
for node = 1:Nnds               %  process nodes
  if isempty(node_label)
    fprintf(fid, '%d;\n', node);
  else
    fprintf(fid, '%d [ label = "%s" ];\n', node, node_label{node});
  end
end
for node1 = 1:Nnds   % process edges
  arcs = find(adj(node1,:));         % children(adj, node);
  for node2 = arcs
    if  ~isempty(arc_label)
      fprintf(fid, edgeformat,node1,node2,arc_label{node1,node2});
    else
      fprintf(fid, edgeformat, node1, node2);
    end
  end
end
fprintf(fid, '}');
fclose(fid);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%HEATTRANSFER Heat transfer physics mode definition.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function w = getplot
f = figure;
plot(1:10);
w = webfigure(f);
close(f);

author: Sable

repo: mcbench-benchmarks


text: %MIN Dataset overload

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [solution, model] = solveCobraCPLEX(model, printLevel, basisReuse, conflictResolve, contFunctName, minNorm)
% Calls CPLEX to solve an LP or QP problem using the matlab API to cplex written by ILOG
%
% USAGE:
%
%    [solution, model] = solveCobraCPLEX(model, printLevel, basisReuse, conflictResolve, contFunctName, minNorm)
%
% INPUT:
%    model:              structure with mandatory and optional fields:
%
%                          * .A or .S: - `m` x `n` LHS matrix
%                          * .b - `m x 1` RHS vector
%                          * .c - `n x 1` Objective coeff vector
%                          * .lb - `n x 1` Lower bound vector
%                          * .ub - `n x 1` Upper bound vector
%                          * .osense - scalar Objective sense (-1 max, +1 min)
%                          * .rxns - (optional) cell array of reaction abbreviations (necessary for
%                            making a readable confilict resolution file).
%                          * .csense - (optional) Constraint senses, a string containting the constraint sense for
%                            each row in `A` ('E', equality, 'G' greater than, 'L' less than).
%                          * .LPBasis - (optional) Basis from previous solution of similar LP problem.
%                            See `basisReuse`
%
% OPTIONAL INPUTS:
%    printLevel:         Printing level in the CPLEX m-file and CPLEX C-interface.
%
%                          * 0 - Silent
%                          * 1 - Warnings and Errors
%                          * 2 - Summary information (Default)
%                          * 3 - More detailed information
%                          * > 10 - Pause statements, and maximal printing (debug mode)
%    basisReuse:         0 - Use this for one of solution of an LP (Default);
%                        1 - Returns a basis for reuse in the next LP i.e. outputs `model.LPBasis`
%    conflictResolve:    0 (Default);
%                        1 If LP problem is proven to be infeasible by CPLEX,
%                        it will print out a 'conflict resolution file',
%                        which indicates the irreducible infeasible set of
%                        equaltiy & inequality constraints that together,
%                        combine to make the problem infeasible. This is
%                        useful for debugging an LP problem if you want to
%                        try to resolve a constraint conflict
%    contFunctName:
%
%                        1. contFunctName = [] Use all default CLPEX control parameters, (Default);
%                        2.  contFunctName = someString e.g. 'someFunctionName'
%                        uses the user specified control parameters defined
%                        in `someFunctionName.m`;
%                        (see template function CPLEXParamSet for details).
%                        3. contFunctName = `cpxControl` structure (output from a file like `CPLEXParamSet.m`)
%    minNorm:            {(0), 1 , `n x 1` vector} If not zero then, minimise the Euclidean length
%                        of the solution to the LP problem. Gives the same objective,
%                        but minimises the square of flux. `minNorm` ~1e-6 should be
%                        high enough for regularisation yet keep the same objective
%
% OUTPUT:
%    solution:           Structure containing the following fields describing a LP solution:
%
%                          * .full:               Full LP solution vector
%                          * .obj:                Objective value
%                          * .rcost:              Lagrangian multipliers to the simple inequalties (Reduced costs)
%                          * .dual:               Lagrangian multipliers to the equalities
%                          * .nInfeas:            Number of infeasible constraints
%                          * .sumInfeas:          Sum of constraint violation
%                          * .stat:               COBRA Standardized solver status code:
%
%                            * 1 - Optimal solution
%                            * 2 - Unbounded solution
%                            * 0 - Infeasible
%                            * -1 - No solution reported (timelimit, numerical problem etc)
%                          * .origStat:           CPLEX status code. Use `cplexStatus(solution.origStat)` for
%                            more information from the CPLEX solver
%                          * .solver              solver used by `cplex`
%                          * .time                time taken to solve the optimization problem
%
% OPTIONAL OUTPUT:
%    model:              with field:
%
%                          * .LPBasis - When input basisReuse=1, we return a basis for reuse in the next LP
%
% CPLEX consists of 4 different LP solvers which can be used to solve sysbio optimization problems
% you can control which of the solvers, e.g. simplex vs interior point solver using the
% CPLEX control parameter cpxControl.LPMETHOD. At the moment, the solver is
% automatically chosen for you.
%
% .. Author: - Ronan Fleming 23 Oct  09  ILOG-CPLEX 12.1 via  matlab API

global CBT_LP_SOLVER
global CBT_MILP_SOLVER
global CBT_QP_SOLVER
global CBT_MIQP_SOLVER

% Note: Certain parts of the conflict resolution only work with tomlab cplex: http://tomwiki.com/CPLEX_Conflict_refiner
tomlabSelected = false;
if (strcmp(CBT_LP_SOLVER, 'tomlab_cplex') || strcmp(CBT_QP_SOLVER, 'tomlab_cplex') || strcmp(CBT_MILP_SOLVER, 'tomlab_cplex') || strcmp(CBT_MIQP_SOLVER, 'tomlab_cplex'))
    tomlabSelected = true;
end

if ~exist('printLevel','var')
    printLevel=0;
end
if ~exist('basisReuse','var')
    basisReuse=0;
end
if ~exist('conflictResolve','var')
    conflictResolve=0;
end

if ~exist('minNorm','var')
    minNorm=0;
end

if basisReuse
    if isfield(model,'LPBasis')
        basis=model.LPBasis;
        %use advanced starting information when optimization is initiated.
        %cpxControl.ADVIND=1;
    else
        basis=[];
    end
else
    basis=[];
    %do not use advanced starting information when optimization is initiated.
    cplex.Param.ADVIND=0;
end

if ~isfield(model,'A')
    if ~isfield(model,'S')
        error('Equality constraint matrix must either be a field denoted A or S.')
    end
    model.A=model.S;
end

if ~isfield(model,'csense')
    nMet=size(model.A);
    if printLevel>0
        fprintf('%s\n','Assuming equality constraints, i.e. S*v=b');
    end
    %assuming equality constraints
    model.csense(1:nMet,1)='E';
end

if ~isfield(model,'osense')
    %assuming maximisation
    model.osense=-1;
    if printLevel>0
        fprintf('%s\n','Assuming maximisation of objective');
    end
end

if size(model.A,2)~=length(model.c)
    error('dimensions of A & c are inconsistent');
end

if size(model.A,2)~=length(model.lb) || size(model.A,2)~=length(model.ub)
    error('dimensions of A & bounds are inconsistent');
end

%Conflict groups descriptor (cpxBuildConflict can be used to generate the input). Set this if
%conflict refinement is desired in the case that infeasibility is detected
%by CPLEX.
if conflictResolve
    if tomlabSelected
        [m_lin,n]=size(model.A);
        m_quad=0;
        m_sos=0;
        m_log=0;
        %determines how elaborate the output is
        mode='full';%'minimal';
        fprintf('%s\n%s\n','Building Structure for Conflict Resolution...','... this slows CPLEX down so should not be used for repeated LP');
        cpxBuildConflict(n,m_lin,m_quad,m_sos,m_log,mode);
    end
    prefix=pwd;
    suffix='CPLEX_conflict_file.txt';
end

% Initialize the CPLEX object
try
	cplex = Cplex('fba');
catch ME
	error('CPLEX not installed or licence server not up')
end

% Now populate the problem with the data
cplex.Model.sense = 'minimize';
if model.osense == 1
    %minimise linear objective
    cplex.Model.obj   = model.c;
else
    %maximise linear objective by reversing sign
    cplex.Model.obj   = - model.c;
end

cplex.Model.lb    = model.lb;
cplex.Model.ub    = model.ub;
cplex.Model.A     = model.A;

%cplex interface
if isfield(model,'csense')
    %set up constant vectors for CPLEX
    b_L(model.csense == 'E',1) = model.b(model.csense == 'E');
    b_U(model.csense == 'E',1) = model.b(model.csense == 'E');
    b_L(model.csense == 'G',1) = model.b(model.csense == 'G');
    b_U(model.csense == 'G',1) =  Inf;
    b_L(model.csense == 'L',1) = -Inf;
    b_U(model.csense == 'L',1) = model.b(model.csense == 'L');
    cplex.Model.lhs   = b_L;
    cplex.Model.rhs   = b_U;
else
    cplex.Model.lhs   = model.b;
    cplex.Model.rhs   = model.b;
end

%quadratic constraint matrix, size n x n
if sum(minNorm)~=0
    if length(minNorm)==1
        % same weighting of min norm for all variables
        cplex.Model.Q=model.osense*speye(length(model.c))*minNorm;
    else
        if length(minNorm)~=length(model.c)
            error('Either minNorm is a scalar, or is an n x 1 vector')
        else
            % individual weighting of min norm for all variables
            cplex.Model.Q=model.osense*spdiags(minNorm,0,length(model.c),length(model.c));
        end
    end
end

%set the solver parameters
if exist('contFunctName','var')
    if isstruct(contFunctName)
        cplex.Param=contFunctName;
    else
        if ~isempty(contFunctName)
            %calls a user specified function to create a CPLEX control structure
            %specific to the users problem. A TEMPLATE for one such function is
            %CPLEXParamSet
            %e.g. Param.lpmethod.Cur=0;
            cplex.Param=Param;
        end
    end
end

if printLevel==0
    cplex.DisplayFunc=[];
else
    %print level
    cplex.Param.barrier.display.Cur = printLevel;
    cplex.Param.simplex.display.Cur = printLevel;
    cplex.Param.sifting.display.Cur = printLevel;
end

%limit the processing to 3 threads
cplex.Param.threads.Cur = 3;

% Optimize the problem
cplex.solve();
solution.origStat   = cplex.Solution.status;

if printLevel>0 && solution.origStat~=1 && tomlabSelected
    %use tomlab code to print out exit meassage
    [ExitText,ExitFlag] = cplexStatus(solution.origStat);
    solution.ExitText=ExitText;
    solution.ExitFlag=ExitFlag;
    if any(model.c~=0)
        if isfield(cplex.Solution, 'objval')
            fprintf('\n%s%g\n',[ExitText ', Objective '], model.osense*cplex.Solution.objval);
        else
            fprintf('\n%s\n', ExitText);
        end
    end
end

if solution.origStat==1
    %extract the solution
    solution.obj        = model.osense*cplex.Solution.objval;
    solution.full       = cplex.Solution.x;
    solution.rcost      = cplex.Solution.reducedcost;
    solution.dual       = cplex.Solution.dual;
    solution.nInfeas    = NaN;
    solution.sumInfeas  = NaN;
    solution.solver     = cplex.Solution.method;
    solution.time       = cplex.Solution.time;
else
    solution.time=NaN;
    %conflict resolution
    if conflictResolve ==1
        cplex.refineConflict();
        cplex.writeConflict(suffix);
        if isfield(model,'mets') && isfield(model,'rxns')
            %this code reads the conflict resolution file and replaces the
            %arbitrary names with the abbreviations of metabolites and reactions
            [nMet,nRxn]=size(model.A);
            totAbbr=nMet+nRxn;
            conStrFind=cell(nMet+nRxn,1);
            conStrReplace=cell(nMet+nRxn,1);
            %only equality constraint rows
            for m=1:nMet
                conStrFind{m,1}=['c' int2str(m) ':'];
                conStrReplace{m,1}=[model.mets{m} ':  '];
            end
            %reactions
            for n=1:nRxn
                conStrFind{nMet+n,1}=['x' int2str(n) ' '];
                conStrReplace{nMet+n,1}=[model.rxns{n} ' '];
            end
            fid1 = fopen(suffix);
            fid2 = fopen(['COBRA_' suffix], 'w');
            while ~feof(fid1)
                tline{1}=fgetl(fid1);
                %replaces all occurrences of the string str2 within string str1
                %with the string str3.
                %str= strrep(str1, str2, str3)
                for t=1:totAbbr
                    tline= strrep(tline, conStrFind{t}, conStrReplace{t});
                end
                fprintf(fid2,'%s\n', tline{1});
            end
            fclose(fid1);
            fclose(fid2);
        else
            warning('Need reaction and metabolite abbreviations in order to make a readable conflict resolution file');
        end
        fprintf('%s\n',['Conflict resolution file written to: ' prefix filesep 'COBRA_' suffix]);
        fprintf('%s\n%s\n','The Conflict resolution file gives an irreducible infeasible subset ','of constraints which are making this LP Problem infeasible');
    else
        if printLevel>0
            fprintf('%s\n','No conflict resolution file. Perhaps set conflictResolve = 1 next time.');
        end
    end
end

% Try to give back COBRA Standardized solver status:
%           1   Optimal solution
%           2   Unbounded solution
%           0   Infeasible
%           -1  No solution reported (timelimit, numerical problem etc)
if solution.origStat==1
    solution.stat = 1;
else
    %use tomlab code to print out exit meassage
    if tomlabSelected
        [ExitText,ExitFlag] = cplexStatus(solution.origStat);
        solution.ExitText=ExitText;
        solution.ExitFlag=ExitFlag;
        if any(model.c~=0) && isfield(cplex,'Solution')
            if isfield(cplex.Solution,'x')
                fprintf('\n%s%g\n',[ExitText ', Objective '],  model.c'*cplex.Solution.x);
            end
        end
    end
    if solution.origStat==2
        solution.stat = 2;
    else
        if solution.origStat==3
            solution.stat = 0;
        else
            %this is a conservative view
            solution.stat = -1;
        end
    end
end

%return basis
if basisReuse
    model.LPBasis=basis;
end

if sum(minNorm)~=0 && printLevel>0
    fprintf('%s\n','This objective corresponds to a flux with minimum Euclidean norm.');
    if length(minNorm)==1
        fprintf('%s%d%s\n','The weighting for minimising the norm was ',minNorm,'.');
    else
        fprintf('%s%d%s\n','The sum of the weighting for minimising the norm was ',sum(minNorm),'.');
    end
    fprintf('%s\n','Check that the objective is the same without minimising the norm.');
end

author: opencobra

repo: cobratoolbox


text:
function [res,resL2,qualMeasOut]=SART(proj,geo,angles,niter,varargin)
% SART solves Cone Beam CT image reconstruction using Oriented Subsets
%              Simultaneous Algebraic Reconstruction Technique algorithm
%
%   SART(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem
%   using the projection data PROJ taken over ALPHA angles, corresponding
%   to the geometry described in GEO, using NITER iterations.
%
%   SART(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The
%   possible options in OPT are:
%
%
%   'lambda':      Sets the value of the hyperparameter. Default is 1
%
%   'lambda_red':  Reduction of lambda. Every iteration
%                  lambda=lambdared*lambda. Default is 0.99
%
%   'Init':        Describes different initialization techniques.
%                  'none'     : Initializes the image to zeros (default)
%                  'FDK'      : Initializes image to FDK reconstruction
%                  'multigrid': Initializes image by solving the problem in
%                               small scale and increasing it when relative
%                               convergence is reached.
%                  'image'    : Initialization using a user specified
%                               image. Not recommended unless you really
%                               know what you are doing.
%   'InitImg'      an image for the 'image' initialization. Avoid.
%
%   'Verbose'      1 or 0. Default is 1. Gives information about the
%                  progress of the algorithm.
%   'QualMeas'     Asks the algorithm for a set of quality measurement
%                  parameters. Input should contain a cell array of desired
%                  quality measurement names. Example: {'CC','RMSE','MSSIM'}
%                  These will be computed in each iteration.
% 'OrderStrategy'  Chooses the subset ordering strategy. Options are
%                  'ordered' : uses them in the input order, but divided
%                  'random'  : orders them randomly
%                  'angularDistance': chooses the next subset with the
%                                     biggest angular distance with the ones used.
% 'redundancy_weighting': true or false. Default is true. Applies data
%                         redundancy weighting to projections in the update step
%                         (relevant for offset detector geometry)
%  'groundTruth'  an image as grounf truth, to be used if quality measures
%                 are requested, to plot their change w.r.t. this known
%                 data.
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------

%% Deal with input parameters
blocksize=1;
[lambda,res,lambdared,verbose,QualMeasOpts,OrderStrategy,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin);

measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:)));
if ~any(isnan(gt(:)))
    QualMeasOpts{end+1}='error_norm';
    res_prev=gt;
    clear gt
end
if nargout<3 && measurequality
    warning("Image metrics requested but none catched as output. Call the algorithm with 3 outputs to store them")
    measurequality=false;
end
qualMeasOut=zeros(length(QualMeasOpts),niter);

resL2=zeros(1,niter);
if nargout>1
    computeL2=true;
else
    computeL2=false;
end

[alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy);
index_angles=cell2mat(orig_index);
angles_reorder=cell2mat(alphablocks);

% does detector rotation exists?
if ~isfield(geo,'rotDetector')
    geo.rotDetector=[0;0;0];
end
%% Create weighting matrices

% Projection weight, W

% Projection weight, W
W=computeW(geo,angles,gpuids);

% Back-Projection weight, V
V=computeV(geo,angles,alphablocks,orig_index,'gpuids',gpuids);

if redundancy_weights
    % Data redundancy weighting, W_r implemented using Wang weighting
    % reference: https://iopscience.iop.org/article/10.1088/1361-6560/ac16bc

    num_frames = size(proj,3);
    W_r = redundancy_weighting(geo);
    W_r = repmat(W_r,[1,1,num_frames]);
    % disp('Size of redundancy weighting matrix');
    % disp(size(W_r));
    W = W.*W_r; % include redundancy weighting in W
end

clear A x y dx dz;
%% hyperparameter stuff
nesterov=false;
if ischar(lambda)&&strcmp(lambda,'nesterov')
    nesterov=true;
    lambda=(1+sqrt(1+4))/2;
    gamma=0;
    ynesterov=zeros(size(res),'single');
    ynesterov_prev=ynesterov;
end
%% Iterate
offOrigin=geo.offOrigin;
offDetector=geo.offDetector;
rotDetector=geo.rotDetector;
DSD=geo.DSD;
DSO=geo.DSO;
% TODO : Add options for Stopping criteria
for ii=1:niter
    if (ii==1 && verbose==1);tic;end
    % If quality is going to be measured, then we need to save previous image
    % THIS TAKES MEMORY!
    if measurequality && ~strcmp(QualMeasOpts,'error_norm')
        res_prev = res; % only store if necesary
    end

    % reorder angles

    for jj=index_angles
        if size(offOrigin,2)==size(angles,2)
            geo.offOrigin=offOrigin(:,jj);
        end
        if size(offDetector,2)==size(angles,2)
            geo.offDetector=offDetector(:,jj);
        end
        if size(rotDetector,2)==size(angles,2)
            geo.rotDetector=rotDetector(:,jj);
        end
        if size(DSD,2)==size(angles,2)
            geo.DSD=DSD(jj);
        end
        if size(DSO,2)==size(angles,2)
            geo.DSO=DSO(jj);
        end
        % --------- Memory expensive-----------

        %         proj_err=proj(:,:,jj)-Ax(res,geo,angles(jj));       %                                 (b-Ax)
        %         weighted_err=W(:,:,jj).*proj_err;                   %                          W^-1 * (b-Ax)
        %         backprj=Atb(weighted_err,geo,angles(jj));           %                     At * W^-1 * (b-Ax)
        %         weigth_backprj=bsxfun(@times,1./V(:,:,jj),backprj); %                 V * At * W^-1 * (b-Ax)
        %         res=res+lambda*weigth_backprj;                      % x= x + lambda * V * At * W^-1 * (b-Ax)
        %------------------------------------
        %--------- Memory cheap(er)-----------
        if nesterov
            % The nesterov update is quite similar to the normal update, it
            % just uses this update, plus part of the last one.
            ynesterov=res+ bsxfun(@times,1./V(:,:,jj),Atb(W(:,:,index_angles(:,jj)).*(proj(:,:,index_angles(:,jj))-Ax(res,geo,angles_reorder(:,jj),'gpuids',gpuids)),geo,angles_reorder(:,jj),'gpuids',gpuids));
            res=(1-gamma)*ynesterov+gamma*ynesterov_prev;
        else
            res=res+lambda* bsxfun(@times,1./V(:,:,jj),Atb(W(:,:,index_angles(:,jj)).*(proj(:,:,index_angles(:,jj))-Ax(res,geo,angles_reorder(:,jj),'gpuids',gpuids)),geo,angles_reorder(:,jj),'gpuids',gpuids));
        end
        if nonneg
            res=max(res,0);
        end
    end

    % If quality is being measured
    if measurequality
        qualMeasOut(:,ii)=Measure_Quality(res,res_prev,QualMeasOpts);
    end

    if nesterov
        gamma=(1-lambda);
        lambda=(1+sqrt(1+4*lambda^2))/2;
        gamma=gamma/lambda;
    else
        lambda=lambda*lambdared;
    end
    if computeL2 || nesterov
        geo.offOrigin=offOrigin;
        geo.offDetector=offDetector;
        geo.DSD=DSD;
        geo.rotDetector=rotDetector;
        resL2(ii)=im3Dnorm(proj-Ax(res,geo,angles,'gpuids',gpuids),'L2'); % Compute error norm2 of b-Ax
        % If the error is not minimized.
        if  ii~=1 && resL2(ii)>resL2(ii-1)
            if verbose
                disp(['Convergence criteria met, exiting on iteration number:', num2str(ii)]);
            end
            return
        end
    end

    if (ii==1 && verbose==1)
        expected_time=toc*niter;
        disp('SART');
        disp(['Expected duration   :    ',secs2hms(expected_time)]);
        disp(['Expected finish time:    ',datestr(datetime('now')+seconds(expected_time))]);
        disp('');
    end
end

end

function initres=init_multigrid(proj,geo,alpha)

finalsize=geo.nVoxel;
% start with 64
geo.nVoxel=[64;64;64];
geo.dVoxel=geo.sVoxel./geo.nVoxel;
if any(finalsize<geo.nVoxel)
    initres=zeros(finalsize');
    return;
end
niter=100;
initres=zeros(geo.nVoxel');
while ~isequal(geo.nVoxel,finalsize)

    % solve subsampled grid
    initres=SART(proj,geo,alpha,niter,'Init','image','InitImg',initres,'Verbose',0,'gpuids',gpuids);

    % Get new dims.
    geo.nVoxel=geo.nVoxel*2;
    geo.nVoxel(geo.nVoxel>finalsize)=finalsize(geo.nVoxel>finalsize);
    geo.dVoxel=geo.sVoxel./geo.nVoxel;
    % Upsample!
    % (hopefully computer has enough memory............)
    [y, x, z]=ndgrid(linspace(1,size(initres,1),geo.nVoxel(1)),...
        linspace(1,size(initres,2),geo.nVoxel(2)),...
        linspace(1,size(initres,3),geo.nVoxel(3)));
    initres=interp3(initres,x,y,z);
    clear x y z
end
end

function [lambda,res,lambdared,verbose,QualMeasOpts,OrderStrategy,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,alpha,argin)
opts={'lambda','init','initimg','verbose','lambda_red','qualmeas','orderstrategy','nonneg','gpuids','redundancy_weighting','groundtruth'};
defaults=ones(length(opts),1);
% Check inputs
nVarargs = length(argin);
if mod(nVarargs,2)
    error('TIGRE:SART:InvalidInput','Invalid number of inputs')
end

% check if option has been passed as input
for ii=1:2:nVarargs
    ind=find(ismember(opts,lower(argin{ii})));
    if ~isempty(ind)
        defaults(ind)=0;
    else
        error('TIGRE:SART:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]);
    end
end

for ii=1:length(opts)
    opt=opts{ii};
    default=defaults(ii);
    % if one option is not default, then extract value from input
    if default==0
        ind=double.empty(0,1);jj=1;
        while isempty(ind)
            ind=find(isequal(opt,lower(argin{jj})));
            jj=jj+1;
        end
        if isempty(ind)
            error('TIGRE:SART:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]);
        end
        val=argin{jj};
    end

    switch opt
        % % % % % % % Verbose
        case 'verbose'
            if default
                verbose=1;
            else
                verbose=val;
            end
            if ~is2014bOrNewer
                warning('TIGRE: Verbose mode not available for older versions than MATLAB R2014b');
                verbose=false;
            end
            % % % % % % % hyperparameter, LAMBDA
        case 'lambda'
            if default
                lambda=1;
            elseif ischar(val)&&strcmpi(val,'nesterov')
                lambda='nesterov'; % just for lowercase/uppercase
            elseif length(val)>1 || ~isnumeric(val)
                error('TIGRE:SART:InvalidInput','Invalid lambda')
            else
                lambda=val;
            end
        case 'lambda_red'
            if default
                lambdared=1;
            else
                if length(val)>1 || ~isnumeric(val)
                    error('TIGRE:SART:InvalidInput','Invalid lambda')
                end
                lambdared=val;
            end
        case 'init'
            res=[];
            if default || strcmp(val,'none')
                res=zeros(geo.nVoxel','single');
                continue
            end
            if strcmp(val,'FDK')
                res=FDK(proj,geo,alpha);
                continue
            end
            if strcmp(val,'multigrid')
                res=init_multigrid(proj,geo,alpha);
                continue
            end
            if strcmp(val,'image')
                initwithimage=1; % it is used (10 lines below)
                continue
            end
            if isempty(res)
                error('TIGRE:SART:InvalidInput','Invalid Init option')
            end
            % % % % % % % ERROR
        case 'initimg'
            if default
                continue
            end
            if exist('initwithimage','var')
                if isequal(size(val),geo.nVoxel')
                    res=single(val);
                else
                    error('TIGRE:SART:InvalidInput','Invalid image for initialization');
                end
            end
        case 'qualmeas'
            if default
                QualMeasOpts={};
            else
                if iscellstr(val)
                    QualMeasOpts=val;
                else
                    error('TIGRE:SART:InvalidInput','Invalid quality measurement parameters');
                end
            end
        case 'orderstrategy'
            if default
                OrderStrategy='random';
            else
                OrderStrategy=val;
            end
        case 'nonneg'
            if default
                nonneg=true;
            else
                nonneg=val;
            end
        case 'gpuids'
            if default
                gpuids = GpuIds();
            else
                gpuids = val;
            end
        case 'redundancy_weighting'
            if default
                redundancy_weights = true;
            else
                redundancy_weights = val;
            end
        case 'groundtruth'
            if default
                gt=nan;
            else
                gt=val;
            end
        otherwise
            error('TIGRE:SART:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option']);
    end
end

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
%%*******************************************************************
%%  Convert problem from SeDuMi format to DSDP format
%%
%%  [AC,b] = readsedumi(At,b,C,K)
%%
%%  Input: At, b, C, K = Data in SeDuMi format.
%%
%%
%% DSDP: version 5.0
%% Copyright (c) 2004 by
%% S. Benson Y. Ye
%% Last modified: 2 Jan 04
%%******************************************************************

  function [AC,b]=readsedumi(At,bb,c,K);
%%
%%  First, load the matlab file containing At, c, b, and K
%%
  if (size(c,1) == 1), c = c'; end;
  if (size(bb,1) == 1), bb = bb'; end;
  [nn,mm] = size(At); if (max(size(c)) == 1); c = c*ones(nn,1); end;

  if ~isfield(K,'l'); K.l = 0; end
  if ~isfield(K,'q'); K.q = 0; end
  if ~isfield(K,'s'); K.s = 0; end
  if K.l == 0 | isempty(K.l); K.l = 0; end;
  if sum(K.q) == 0 | isempty(K.q); K.q = 0; end
  if sum(K.s) == 0 | isempty(K.s); K.s = 0; end
%%
%%
%%
  AC=cell(1,3);
block=1;
top=1;

   m = length(bb);
   b = -bb;

   if (K.l > 0)
      AC{block,1} = 'LP';
      AC{block,2} = K.l;
      A = At(top:top+n-1,:);
      CC = c(top:top+n-1,:);
      AC{block,3} = [-A(:,1:m) CC(:,1)];
      block = block+1;
      top = top+K.l;
   end
   if (K.q > 0)
      error(' Cannot accept SOCP cones\n');
   end

   if (K.s > 0)
      for i = 1:length(K.s)
         n = K.s(i);
         AC{block,1} = 'SDP';
         AC{block,2} = n;
         A = At(top:top+n^2-1,:);
         CC = c(top:top+n^2-1,:);
         indicies = triu(reshape(1:n^2,n,n));
         indicies = indicies(find(indicies));
         AC{block,3} = [-A(indicies,1:m) CC(indicies,1)];
         block = block+1;
         top = top+n*n;
       end
   end

%%

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus