text:
function e = end(X,k,n)
%END Last index of indexing expression for tensor.
%
%   The expression X(end,:,:) will call END(X,1,3) to determine
%   the value of the first index.
%
%   See also TENSOR, TENSOR/SUBSREF, END.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if n > ndims(X)
  error('Subscript out of range.');
end
e = X.size(k);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [minidx, inputs, outputs] = utl_gridsearch(arg1, varargin)
% Exhaustive search over multiple possible arguments to a function.
% [MinIndex, Inputs, Outputs] = utl_gridsearch(Function, Domain...)
% [MinIndex, Inputs, Outputs] = utl_gridsearch(RangeFormat, Function, Domain...)
% [MinIndex, Inputs, Outputs] = utl_gridsearch(Options, Domain...)
%
% Grid search is a brute-force and assuption-free approach to optimizing parameters of a function.
% It amounts to testing all the possibilities (out of a predetermined set), obtaining the function
% value for each, and returning the parameter combination which gave the best result (here: smallest
% function value).
%
% The arguments to try for the given function are specified as a regular argument list, in which,
% depending on the chosen RangeFormat, either arrays among the arguments are interpreted as multiple
% possibilities to try, or arguments which contain a search() expression/data structure are
% interpreted as defining multiple possibilities. The most convenient way to specify multiple
% possibilities is by using the special function search(), which packs its multiple arguments into
% an annotated container dat structure which is interpreted by utl_gridsearch as a list of
% possibilities to try.
%
% Grid search can be run in parallel across multiple processes/machines, by passing the appropriate
% parallelization options.
%
%
% In:
%   RangeFormat : argument range format, optional; either 'direct' (search ranges are directly
%                 specified as arrays), or 'clauses' (search ranges are specified using the search()
%                 clause). default: 'direct'
%
%   Function    : the function to optimize; should always return a real number (or NaN) as first
%                 output, taken as the objective value of the function; may have additional output
%                 arguments, which are also captured and optionally returned by utl_gridsearch
%
%   Options     : cell array of name-value pairs to specify detailed options (alternative to
%                 RangeFormat/Function); possible names include:
%                 'func'    : the function to optimize (see Function, mandatory)
%
%                 'argform' : argument range format, either 'direct' or 'clauses' (see RangeFormat,
%                             default: direct)
%
%                 parallelization options (see par_beginschedule)
%                 'engine_gs' : the parallelization engine to be used for the grid search
%                               (default: 'local')
%                 'pool'    : node pool to be used for parallelization, when using the BLS scheduler
%                             (default: 'global')
%                 'policy'  : scheduling policy to be used, when using the BLS scheduler
%                             (default: 'global')
%
%   Domain...   : one-dimensional argument ranges to be searched;
%                 * in 'direct' mode, each array determines the possible values at that place in the
%                   function's argument list, all combinations of specified values are tried; for
%                   cell arrays, the search runs through all given cell contents at that position
%                 * in 'clauses' mode, search ranges for a given argument are specified using the
%                   search() expression in that place of the argument list
%
% Out:
%   MinIndex    : index of the minimum function return value (first return value if multiple ones
%                 are returned), out of all tried executions (see example); can be used as index
%                 into Inputs and/or Outputs
%
%   Inputs      : cell array of all tried input combinations (each a cell array of arguments)
%
%   Outputs     : cell array of all received function outputs, one for each input (each a cell array
%                 of outputs, since the function can have multiple outputs)
%
% Example:
%   % For the four equivalent calls,
%     utl_gridsearch(@f, 1:2, {'a',{'b'},5:7,{5:7}}, [], {}, 'xy', {[1 2 3]});
%     utl_gridsearch('direct', @f, 1:2, {'a',{'b'},5:7,{5:7}}, [], {}, 'xy', {[1 2 3]});
%     utl_gridsearch('clauses', @f, search(1,2), search('a',{'b'},5:7,{5:7}), [], {}, search('x','y'), [1 2 3]);
%     utl_gridsearch({'argform','clauses', 'func',@f}, search(1,2), search('a',{'b'},5:7,{5:7}), [], {}, search('x','y'), [1 2 3]);
%
%   % ... the following 16 = 2x4x1x1x2x1 executions of f are compared:
%       f(1,'a',[],{},'x', [1 2 3])
%       f(1,'a',[],{},'y', [1 2 3])
%       f(1,{'b'},[],{},'x', [1 2 3])
%       f(1,{'b'},[],{},'y', [1 2 3])
%       f(1,5:7,[],{},'x', [1 2 3])
%       f(1,5:7,[],{},'y', [1 2 3])
%       f(1,{5:7},[],{},'x', [1 2 3])
%       f(1,{5:7},[],{},'y', [1 2 3])
%       f(2,'a',[],{},'x', [1 2 3])
%       f(2,'a',[],{},'y', [1 2 3])
%       f(2,{'b'},[],{},'x', [1 2 3])
%       f(2,{'b'},[],{},'y', [1 2 3])
%       f(2,5:7,[],{},'x', [1 2 3])
%       f(2,5:7,[],{},'y', [1 2 3])
%       f(2,{5:7},[],{},'x', [1 2 3])
%       f(2,{5:7},[],{},'y', [1 2 3])
%
% See also:
%   search
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-06
dp;

% translate the possible formats for the control arguments into the cell-array-of-options form
if iscell(arg1)
    % utl_gridsearch({Options...}, Domain...)
    args = arg1;
elseif isa(arg1,'function_handle')
    % utl_gridsearch(Function, Domain...)
    args = {'func', arg1};
elseif ischar(arg1)
    % utl_gridsearch(RangeFormat, Function, Domain...)
    args = {'argform',arg1,'func',varargin{1}}; varargin = varargin(2:end);
else
    error('Unexpected control options format: %s',hlp_tostring(arg1));
end

% read options
opts = hlp_varargin2struct(args, 'func',mandatory, 'argform','direct', 'engine_gs','local', 'pool','global', 'policy','global');
if ~isa(opts.func,'function_handle')
    error('The Function argument must be a function handle.'); end

% rewrite arguments if given as clauses
if strcmp(opts.argform,'clauses')
    for i=1:length(varargin)
        varargin{i} = hlp_flattensearch(varargin{i},'cell'); end
elseif ~strcmp(opts.argform,'direct')
    error('Unsupported RangeFormat value: %s (allowed are ''clauses'' and ''direct'')',opts.argform);
end

% determine dimensions of the input grid
dims = cellfun('length',varargin);
pitches = [1 cumprod(max(1,dims))];
combos = prod(max(1,dims));

args = varargin;
for i=1:length(args)
    % pre-unpack cell args for singleton arguments
    if iscell(args{i}) && dims(i) == 1
        args{i} = args{i}{1}; end
end
for c=1:combos
    for i=find(dims>1)
        % select current value for non-singleton arguments
        args{i} = varargin{i}(1+floor(mod((c-1)/pitches(i),dims(i))));
        % and unpack cell args
        if iscell(varargin{i})
            args{i} = args{i}{1}; end
    end
    % generate tasks
    tasks{c} = [{@hlp_wrapresults,opts.func},args]; %#ok<AGROW>
    inputs{c} = args; %#ok<AGROW>
end

% schedule tasks
outputs = par_schedule(tasks,'engine',opts.engine_gs,'pool',opts.pool,'policy',opts.policy);

try
    minidx = argmin(cellfun(@(x)x{1},outputs));
catch
    minidx = NaN;
end

repo: REST


text:
% 5g ldpc encoding and decoding test
% reference: 3GPP TS 38.212 section 5.3.2
% author: Xiao, Shaoning 萧少宁
% license: MIT

function pass_flag = test_ldpc(varargin)

if nargin == 0
    base_graph_index = 1;
    Z_c = 384;
elseif nargin == 1
    base_graph_index = varargin{1};
    Z_c = 384;
elseif nargin == 2
    base_graph_index = varargin{1};
    Z_c = varargin{2};
end

SNR_list = 10;

trial_number = 1;

BLER = zeros(1, length(SNR_list));

TxRx.Decoder.LDPC.Iterations = 30;
TxRx.Decoder.LDPC.Type = 'OMS';
% TxRx.Decoder.LDPC.Type = 'MPA';

if base_graph_index == 1
    load base_graph_1_check_node_list
    base_graph_check_node_list = base_graph_1_check_node_list;
    LDPC.inf_bits = 22 * Z_c;
elseif base_graph_index == 2
    load base_graph_2_check_node_list
    base_graph_check_node_list = base_graph_2_check_node_list;
    LDPC.inf_bits = 10 * Z_c;
else
    error('wrong base graph index.');
end

for SNR_list_index = 1:length(SNR_list)

    sigma_square = 10^(-SNR_list(SNR_list_index)/10);

    for trial = 1:trial_number

        tx_bits = randi([0, 1], LDPC.inf_bits, 1);

        %         [encoded_bits, LDPC.H, Z_c] = ldpc_encode(tx_bits, base_graph_index);
        [encoded_bits, LDPC.H, Z_c] = ldpc_encode_optimized(tx_bits, base_graph_index);

        [LDPC.par_bits, LDPC.tot_bits] = size(LDPC.H);

        symbols = 1 - 2 * encoded_bits;

        noise = randn(size(symbols));

        waveform = symbols + noise * sqrt(sigma_square);

        LLR_received = 2 * waveform / sigma_square;

        rx_bits = decLDPC_layered(TxRx, LDPC, LLR_received.', base_graph_check_node_list, Z_c);

%         f0=1./(1+exp(-2*LLR_received));
%         f1 = 1 - f0;
%         f = [f0, f1].';
%         rx_bits_prime = ldpc_decode(f0,f1, full(LDPC.H), TxRx.Decoder.LDPC.Iterations);
%
%         ldpc_dec = comm.LDPCDecoder(LDPC.H);
%         LLR_received_extended = [zeros(2*Z_c, 1); LLR_received];
%         decoded_bits = step(ldpc_dec, LLR_received_extended);
%
%         [rx_bits_prime, success, k] = ldpc_decode_1(f, full(LDPC.H), 2);
        rx_bits = rx_bits(:);

        if ~isequal(tx_bits, rx_bits)
            BLER(SNR_list_index) = BLER(SNR_list_index) + 1;
        end

%         if ~isequal(rx_bits_prime(1:LDPC.inf_bits), tx_bits)
%             BLER(SNR_list_index) = BLER(SNR_list_index) + 1;
%         end
%
%         if ~isequal(double(decoded_bits), tx_bits)
%             BLER(SNR_list_index) = BLER(SNR_list_index) + 1;
%         end

    end

end

fprintf('BLER is %f\n', BLER/trial_number);

if BLER < 1e-10
    pass_flag = true;
else
    pass_flag = false;
end

end

author: xiaoshaoning


text:
function varargout = lvd_EditVehiclePointGUI(varargin)
% LVD_EDITVEHICLEPOINTGUI MATLAB code for lvd_EditVehiclePointGUI.fig
%      LVD_EDITVEHICLEPOINTGUI, by itself, creates a new LVD_EDITVEHICLEPOINTGUI or raises the existing
%      singleton*.
%
%      H = LVD_EDITVEHICLEPOINTGUI returns the handle to a new LVD_EDITVEHICLEPOINTGUI or the handle to
%      the existing singleton*.
%
%      LVD_EDITVEHICLEPOINTGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in LVD_EDITVEHICLEPOINTGUI.M with the given input arguments.
%
%      LVD_EDITVEHICLEPOINTGUI('Property','Value',...) creates a new LVD_EDITVEHICLEPOINTGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before lvd_EditVehiclePointGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to lvd_EditVehiclePointGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help lvd_EditVehiclePointGUI

% Last Modified by GUIDE v2.5 12-Jan-2021 17:14:15

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lvd_EditVehiclePointGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @lvd_EditVehiclePointGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before lvd_EditVehiclePointGUI is made visible.
function lvd_EditVehiclePointGUI_OpeningFcn(hObject, eventdata, handles, varargin)
    % This function has no output args, see OutputFcn.
    % hObject    handle to figure
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    % varargin   command line arguments to lvd_EditVehiclePointGUI (see VARARGIN)

    % Choose default command line output for lvd_EditVehiclePointGUI
    handles.output = hObject;

    centerUIFigure(hObject);

    point = varargin{1};
    setappdata(hObject, 'point', point);

	populateGUI(handles, point);

    % Update handles structure
    guidata(hObject, handles);

    % UIWAIT makes lvd_EditVehiclePointGUI wait for user response (see UIRESUME)
    uiwait(handles.lvd_EditVehiclePointGUI);

function populateGUI(handles, point)
    set(handles.pointNameText,'String',point.getName());

% --- Outputs from this function are returned to the command line.
function varargout = lvd_EditVehiclePointGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
    if(isempty(handles))
        varargout{1} = false;
    else
        point = getappdata(hObject, 'point');

        point.setName(handles.pointNameText.String);

        varargout{1} = true;
        close(handles.lvd_EditVehiclePointGUI);
    end

% --- Executes on button press in saveAndCloseButton.
function saveAndCloseButton_Callback(hObject, eventdata, handles)
% hObject    handle to saveAndCloseButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    errMsg = validateInputs(handles);

    if(isempty(errMsg))
        uiresume(handles.lvd_EditVehiclePointGUI);
    else
        msgbox(errMsg,'Invalid Point Inputs','error');
    end

function errMsg = validateInputs(handles)
    errMsg = {};

    if(isempty(strtrim(handles.pointNameText.String)))
        errMsg = {'Point name must contain more than white space and must not be empty.'};
    end

% --- Executes on button press in cancelButton.
function cancelButton_Callback(hObject, eventdata, handles)
% hObject    handle to cancelButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    close(handles.lvd_EditVehiclePointGUI);

function pointNameText_Callback(hObject, eventdata, handles)
% hObject    handle to pointNameText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of pointNameText as text
%        str2double(get(hObject,'String')) returns contents of pointNameText as a double

% --- Executes during object creation, after setting all properties.
function pointNameText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to pointNameText see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on key press with focus on lvd_EditVehiclePointGUI or any of its controls.
function lvd_EditVehiclePointGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditVehiclePointGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'enter'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'escape'
            close(handles.lvd_EditVehiclePointGUI);
    end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef SpectrumAnalyzer
    properties
        fs = 16000;
        nCh = 1;
        nFFT = 512;
        nBin = 257;
        frameSize = 400;
        frameShift = 160;
        frameOverlap = 240;
        frameRate = 100;
        window = 'hamming';
        doDithering = 1;
        doDCRemoval = 0;
        useGPU = 0;
    end

    methods
        function obj = SpectrumAnalyzer(fs, nCh)
            obj.fs = fs;
            if nargin<2
                obj.nCh = 1;
            else
                obj.nCh = nCh;
            end
        end

        function spec = Analyze(obj, data)
            spec = sfft_multi(data, obj.frameSize, obj.frameShift, obj.nFFT, obj.window, obj.doDCRemoval, obj.useGPU, obj.doDithering);
            spec = spec(1:obj.nBin,:,:);
        end

        function wav = Synthesize(obj, complexSpectrum)
            [nBin,D2,D3] = size(complexSpectrum);
            if D2>1 && D3>1     % if both D2 and D3 > 1, assume D2 is the number of channels, while D3 is the number of frames (same format as sfft_multi).
                complexSpectrum = permute(complexSpectrum, [1 3 2]);
            end
            for i=1:size(complexSpectrum,3)
                data = squeeze(complexSpectrum(:,:,i));
                data(obj.nFFT/2+2:obj.nFFT,:) = conj(data(obj.nFFT/2:-1:2,:));
                xi = real(ifft(data));
                wav(:,i) = obj.ola(xi);
            end
        end

        function output = DisplaySpectrogramFromWav(obj,wav, newFigure, doCMN, channelSelector)
            if nargin<3; newFigure = 1; end
            if newFigure; figure; end
            if nargin<4; doCMN = 0; end
            if nargin>=5
                data = obj.Analyze(wav(:,channelSelector));
            else
                data = obj.Analyze(wav);
            end
            data = data(1:obj.nBin,:,:);

            [nBin, nCh, nFr] = size(data);
            if doCMN==0
                displayMatrix = log(abs(reshape(data, nBin*nCh, nFr)));
            elseif doCMN == 1
                displayMatrix = CMN(log(abs(reshape(data, nBin*nCh, nFr)))')';
            end
            imagesc(displayMatrix);
            if nargout>0
                output = displayMatrix;
            end
        end

        function output = DisplaySpectrogram(obj,data, newFigure, doCMN)
            if nargin<3; newFigure = 1; end
            if newFigure; figure; end
            if nargin<4; doCMN = 0; end
            [nBin, nCh, nFr] = size(data);
            if doCMN==0
                displayMatrix = log(abs(reshape(data, nBin*nCh, nFr)));
            elseif doCMN == 1
                displayMatrix = CMN(log(abs(reshape(data, nBin*nCh, nFr)))')';
            end
            imagesc(displayMatrix);
            if nargout>0
                output = displayMatrix;
            end
        end
    end

    methods (Access = protected)
        function xfinal = ola(obj, xi)
            nFr = size(xi,2);
            k=1;
            xfinal=zeros(nFr*obj.frameShift,1);
            x_old = zeros(obj.frameOverlap,1);
            for j=1:nFr
                xfinal(k:k+obj.frameOverlap-1) = x_old + xi(1:obj.frameOverlap,j);
                x_old = xi(1+obj.frameShift:obj.frameSize,j);
                k = k + obj.frameShift;
            end
        end
    end
end

author: singaxiong

repo: SignalGraph


text:
function [found] = findUndirectedEdge(edge, edges_id)

id1 = edge(1);
id2 = edge(2);

foundForward = find(edges_id(:,1)==id1 & edges_id(:,2)==id2);
if length(foundForward) == 1
  found = 1;
  return;
end

foundBack = find(edges_id(:,2)==id1 & edges_id(:,1)==id2);
if length(foundBack) == 1
  found = 1;
  return;
end

if length(foundForward) > 1 || length(foundBack) > 1
  error('Multiple edge!')
end

found = 0;

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function command = generate_python_command(script, paths)
% generate_python_command Generate command line for a Python tracker
%
% This function generates the appropritate command string that will
% run the python executable and execute the given script that includes your
% tracker implementation.
%
% Input:
% - script (string): Name of the tracker script to be executed.
% - paths (cell): An array of strings that denote directories to be added to Matlab path.
%
% Output:
% - command (string): Generated command string.
%

trax_python = get_global_variable('trax_python');
python_exec = get_global_variable('python');

% If path to python trax implementation is set then we attempt to export it to tracker
if ~isempty(trax_python)
    paths = cat(1, {trax_python}, paths(:));
end

path_string = strjoin(cellfun(@(p) sprintf('sys.path.append(''%s'');', p), paths, 'UniformOutput', false), '');

% Attempt to locate python interpreter
if isempty(python_exec)
    if ispc()
	    path_separator = ';';
	    exec_name = 'python.exe';
    else
	    path_separator = ':';
	    exec_name = 'python';
    end

    system_paths = strsplit(getenv('PATH'), path_separator);

    for i = 1:numel(system_paths)
        if exist(fullfile(system_paths{i}, exec_name), 'file') == 2
            python_exec = fullfile(system_paths{i}, exec_name);
	    break;
        end
    end

end

if isempty(python_exec)
    error('Unable to locate Python interpreter, please set the "python" global variable manually.');
end

if ispc()
	python_executable = ['"', python_exec, '"'];
	python_flags = {};
else
	python_executable = python_exec;
	python_flags = {};
end

argument_string = strjoin(python_flags, ' ');

python_script = sprintf('import sys; %s import %s', path_string, script);

command = sprintf('%s %s -c "%s"', python_executable, argument_string, python_script);

author: votchallenge

repo: toolkit-legacy


text:
function [Ds, D] = spm_eeg_inv_extract(D)
% Exports source activity using the MAP projector
% FORMAT [Ds] = spm_eeg_inv_extract(D)
% Requires:
%
%     D.inv{i}.source.XYZ   - (n x 3) matrix of MNI coordinates
%
% Optional:
%
%     D.inv{i}.source.rad   - radius (mm) of VOIs (default 5 mm)
%     D.inv{i}.source.label - label(s) for sources (cell array)
%     D.inv{i}.source.fname - output file name
%     D.inv{i}.source.type  - output type ('evoked'/'trials')
%__________________________________________________________________________
% Copyright (C) 2011 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak, Laurence Hunt, Karl Friston
% $Id: spm_eeg_inv_extract.m 7009 2017-02-07 18:14:16Z vladimir $

% SPM data structure
%==========================================================================
try
    inv = D.inv{D.val};
catch
    inv = D.inv{end};
    D.val = numel(D.inv);
end

Ds = [];

% defaults
%--------------------------------------------------------------------------
try, XYZ   = inv.source.XYZ;   catch,  return;                          end
try, rad   = inv.source.rad;   catch,  rad   = 5;                       end
try, label = inv.source.label; catch,  label = {};                      end
try
    fname  = inv.source.fname;
catch
    fname  = fullfile(D.path, ['i' D.fname]);
end
try
    type   =  inv.source.type;
catch
    if isequal(D.type, 'evoked')
        type = 'evoked';
    else
        type = 'trials';
    end
end

if numel(label)<size(XYZ, 1)
    for i = (numel(label)+1):size(XYZ, 1)
        label{i} = ['Source ' num2str(i)];
    end
end

% find relevant vertices
%==========================================================================
vert  = inv.mesh.tess_mni.vert(inv.inverse.Is, :);        % vertices
Ns    = size(XYZ, 1);                        % number of sources
svert = {};
for i = 1:Ns
    dist = sqrt(sum([vert(:,1) - XYZ(i,1), ...
                     vert(:,2) - XYZ(i,2), ...
                     vert(:,3) - XYZ(i,3)].^2, 2));
    if rad > 0
        svert{i} = find(dist < rad);
    else
        [junk,svert{i}] = min(dist);
        XYZ(i, :) = vert(svert{i}, :);
    end
end

% report
%--------------------------------------------------------------------------
Iy = find(cellfun('isempty', svert));
if ~isempty(Iy)
    disp(['No proximal vertices found for source(s) ' num2str(Iy)]);
    disp('These sources will be discarded');

    svert(Iy) = [];
    label(Iy) = [];
    Ns        = numel(label);

    if Ns == 0
        disp('No valid sources - Exiting');
        return;
    end
end

js      = spm_vec(svert);

% get info from inv
%--------------------------------------------------------------------------
scale   = inv.inverse.scale;
J       = inv.inverse.J;
U       = inv.inverse.U; % spatial projector (contains multiple cells if multimodal)
T       = inv.inverse.T; % S in Friston 2008 - temporal projector
TT      = T*T';
M       = inv.inverse.M(js, :);
Ic      = inv.inverse.Ic;
It      = inv.inverse.It;
Np      = length(It);

try
    trial = inv.inverse.trials;
catch
    trial = D.condlist;
end

% get source data
%==========================================================================
MYi    = 1;
switch(type)
    case 'evoked'
        Ne     = length(trial);
        MY     = zeros(size(M,1), Ne*Np);

        for i = 1:Ne
            MY(:, MYi:(MYi+Np-1)) = J{i}(js,:)*T';
            MYi                   = MYi + Np;
        end

        clabel = trial;
    case 'trials'
        Ne     = length(D.indtrial(trial, 'GOOD'));
        MY     = zeros(size(M,1), Ne*Np);
        clabel = {};

        for i = 1:numel(trial)

            c      = D.indtrial(trial{i}, 'GOOD');
            clabel = [clabel D.conditions(c)];

            % conditional expectation of contrast (J*W) and its energy
            %--------------------------------------------------------------
            Nt    = length(c);
            spm_progress_bar('Init',Nt,sprintf('extracting data condition %d',i),'trials');

            for j = 1:Nt

                if ~strcmp(D.modality(1,1), 'Multimodal')

                    % unimodal data
                    %------------------------------------------------------
                    Y     = D(Ic{1},It,c(j))*TT;
                    Y     = U{1}*Y*scale;

                else

                    % multimodal data
                    %------------------------------------------------------
                    for k = 1:length(U)
                        Y       = D(Ic{k},It,c(j))*TT;
                        UY{k,1} = U{k}*Y*scale(k);
                    end
                    Y = spm_cat(UY);
                end

                MY(:, MYi:(MYi+Np-1)) = M*Y;
                MYi                   = MYi+Np;

                spm_progress_bar('Set',j)
            end
            spm_progress_bar('Clear')
        end
end

% compute regional response in terms of first eigenvariate
%==========================================================================
iS    = [0 cumsum(cellfun('length', svert))];
Y     = zeros(Ns, size(MY, 2));

spm_progress_bar('Init', Ns, 'extracting eigenvariates', 'sources');

for i = 1:Ns
    y     = MY((iS(i)+1):iS(i+1),:)';

    [m,n]   = size(y);
    if n>1
        if m > n
            [v,s,v] = svd(y'*y);
            s       = diag(s);
            v       = v(:,1);
            u       = y*v/sqrt(s(1));
        elseif m>1
            [u,s,u] = svd(y*y');
            s       = diag(s);
            u       = u(:,1);
            v       = y'*u/sqrt(s(1));
        end
        d       = sign(sum(v));
        u       = u*d;
        v       = v*d;
        Y(i, :) = u'*sqrt(s(1)/n);
    else
        Y(i, :) = y';
    end

    spm_progress_bar('Set',i);
end

spm_progress_bar('Clear')

% create source dataset
%-----------------------------------------------------------------------
Ds = clone(D, fname, [Ns Np Ne], 3);
Ds = chanlabels(Ds, 1:Ns, label);
Ds = timeonset(Ds, D.time(It(1)));
Ds = chantype(Ds, 1:Ns, 'LFP');
Ds = conditions(Ds, 1:Ne, clabel);
Ds = sensors(Ds, 'MEG', []);
Ds = sensors(Ds, 'EEG', []);
Ds = fiducials(Ds, []);
Ds = rmfield(Ds, fieldnames(Ds));
Ds = condlist(Ds, trial);

spm_progress_bar('Init', Ne, 'writing out data', 'trials');
for i = 1:Ne
    Ds(:,:,i) = Y(:,((i - 1)*Np + 1):i*Np);
    spm_progress_bar('Set',i);
end
spm_progress_bar('Clear');

% store results in the original inv
%-----------------------------------------------------------------------
D.inv{D.val}.source.type  = type;
D.inv{D.val}.source.XYZ   = XYZ;
D.inv{D.val}.source.label = label;
D.inv{D.val}.source.rad   = rad;
D.inv{D.val}.source.fname = fname;

% save headers
%-----------------------------------------------------------------------
save(Ds);

author: spm


text:
function result=lvq_predict(P,Tc,Num_Compet,w1,w2)
n=size(P,2);
result=zeros(2,n);
result(1,:)=Tc;
for i=1:n
    d=zeros(Num_Compet,1);
    for j=1:Num_Compet
        d(j)=sqrt(sse(w1(j,:)'-P(:,i)));
    end
    n1=compet(-1*d);
    n2=purelin(w2*n1);
    result(2,i)=vec2ind(n2);
end
Num_Correct=length(find(result(2,:)==Tc));
accuracy=Num_Correct/n;
disp(['accuracy=' num2str(accuracy*100) '%(' num2str(Num_Correct) '/' num2str(n) ')']);

author: HuangCongQing


text:
%kputpixel 'Display Pixel Values Associated with Region of Image'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros putpixel.pane file
%
% Parameters:
% InputFile: i 'Input File', required: 'File containing pixels to display'
% Integer: valuewidth 'Width of Pixel Value Display ', default: 7: 'Specifies the number of pixel values along width'
% Integer: valueheight 'Height of Pixel Value Display', default: 7: 'Specifies the number of pixel values along height'
% Integer: x 'X Placement', default: -1: 'X location for GUI autoplacement'
% Integer: width 'Window Width ', default: 256: 'Window width'
% Integer: y 'Y Placement', default: -1: 'Y location for GUI autoplacement'
% Integer: height 'Window Height', default: 256: 'Window height'
% Double: update 'Update time', default: 2: 'How often to check input file for modification'
%
% Example: kputpixel(i, {'i','';'valuewidth',7;'valueheight',7;'x',-1;'width',256;'y',-1;'height',256;'update',2})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% putpixel - Display Pixel Values Associated with Region of Image
%
%  DESCRIPTION
% Putpixel is a visualization program that displays value data as
% a grid of values.  These value data points may be pixels in an image.
%
% The number of map pixel values that are displayed along width and height
% can be specified using the [-valuewidth] and [-valueheight] arguments.  These
% options control the number of pixel values that are displayed in
% the horizontal and vertical directions of the display.
%
% It is possible to specify the map pixel background color using the
% [-showcolormap] argument.  When set to "TRUE", the pixel under the
% pointer will be displayed as the background color for the label object in
% which the map values are displayed.  When set to "FALSE", the background
% color is displayed as black.
%
% By default, the visualization display will share its colormap will all
% other applications running at the same time; that is, it does not not use a
% .I "private colormap",
% but rather makes use of the
% .I "default colormap".
% It is possible to specify that the visualization display allocate its
% own private colormap or "grab" all available colors for its own use.
% This is done using the [-priv] argument.  When [-priv] is specified,
% moving the mouse pointer into the display window will cause the display
% to have its private colormap installed; moving the pointer out of the
% display window will cause the private colormap to be de-installed.  This
% results in the "technoflashing" phenomenon characteristic of private
% colormap installation.
%
% Color allocation can also be controlled and set to either
% .I "read-only"
% or
% .I "read/write"
% using the [-alloc] argument.  When set to
% .I "read-only",
% once a color cell has been allocated, it can have its color set only once;
% from then on, the color cell can be shared by multiple applications, but
% not changed.  If the visualization display requires the color displayed
% to change, it must re-allocate the color cell, forcing a re-display of
% data.  This can be an expensive procedure.  In contrast, after a
% .I "read/write"
% color cell is allocated, it can have its color changed at any time
% without re-allocation; the data being displayed does not need to be
% redisplayed, and the color update process is much more efficient.
% However, the colors used in the visualization display cannot be shared by
% other applications.
%
% The input file containing the data values to be displayed is monitored
% by default; the file is checked every 2 seconds for change, and if it
% has been modified, the displayed values are updated accordingly.  The
% interval at which the input file is checked for modification may be
% specified using the [-update] argument.
%
% On creation, the window which displays the value data may be placed
% manually (the default method), or placed automatically.  For automatic
% placement of the window which displays the value data, specify the
% desired location in device (screen) coordinates using the [-x] and
% [-y] arguments.
%
% While the display window should be created with a default size that is
% appropriate to display the data, a width and height for the window can
% be specified explicitly using the [-width] and [-height] arguments.
% Alternatively, you may allow the display window to size itself initially,
% and then use the window manager to interactively modify its size as desired.
%
%
%
%  EXAMPLES
% % putpixel -i image:mandril
% % putpixel -i image:ball
%
%  "SEE ALSO"
% putdata(1)
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kputpixel(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kputpixel(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'valuewidth', 7;'valueheight', 7;'x', -1;'width', 256;'y', -1;'height', 256;'update', 2};
maxval={0,1,1,1000,1000,1000,1000,1};
minval={0,1,1,-1,-1,-1,-1,1};
istoggle=[0,1,1,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','Integer','Integer','Integer','Integer','Integer','Integer','Double'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'putdata"  -pixel'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% test code for tt_qconv()
%
% April 26, 2011
% Vladimir Kazeev
% vladimir.kazeev@gmail.com
% INM RAS
% Moscow, Russia
%

% a,b,eps are parameters
d=10;
a=-1;
b=1;
eps=1E-10;
%

n=2^d;
h=(b-a)/n;
mesh=((1-n:n)'-1/2)*h;
x=exp(-1000*(mesh-0.1).^2);

mesh=((1-n/2:n/2)'-1/2)*h;
y=zeros(n,1);
y(n/2-n/8:n/2+n/8-1)=ones(n/4,1);
y=exp(-1000*(mesh-0.5).^2);

tt_x=full_to_tt(reshape(x,2*ones(1,d+1)),eps);
tt_y=full_to_tt(reshape(y,2*ones(1,d)),eps);

tt_z=tt_qconv(tt_x,tt_y,d);
z=h*tt_qtofull(tt_z,1);

f1=figure; hold on;
set(0,'CurrentFigure',f1)
plot(mesh,x(n/2:n*3/2-1),'r');
f2=figure; hold on;
set(0,'CurrentFigure',f2)
plot(mesh,y,'b');
f3=figure; hold on;
set(0,'CurrentFigure',f3)
plot(mesh,z,'g');

author: oseledets

repo: TT-Toolbox


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------

% This function defines all parameters to tune in the training of MCG
% params = nyud_params(); %sf_mUCM_multi_3sc_u_4r_12k_params();

%% Compute all points from combining these base hierarchies
% The order of combination is hard-coded, see pareto_learning to change it
% The result will be loaded if already computed
function pareto_chose_point(params, num_regions)
  pareto = pareto_learning(params);

  %% Plot pareto front results

  % Play with the following two number to get the desired point

  % Which measure used to get the Pareto
  pareto_meas = 2;  % Jaccard at class level

  % Which particular point in the curve we choose
  pareto_id = -1;   % Around 14k candidates

  if(pareto_id == -1)
    pareto_id = find( pareto.st{pareto_meas}{params.n_r_cand}.mean_n_masks > num_regions, 1, 'first');
  end

  figure;
  x_limit = 1e6;
  titles{1} = {'Maximum achievable quality',['(Pascal segvoc ' params.gt_set_pareto ')'], 'Jaccard at instance level (Ji)'};
  titles{2} = {'Maximum achievable quality',['(Pascal segvoc ' params.gt_set_pareto ')'], 'Jaccard at class level (Jc)'};
  measures  = {'jaccard_object', 'jaccard_class'};
  colors = {'r','g','b','k','m','c'};
  for jj=1:length(measures)
      subplot(1,2,jj); hold on;

      % Singletons, pairs, and triplets
      % Pareto on the same measure that evaluated
      for ii=1:params.n_r_cand
          plot(pareto.st{jj}{ii}.mean_n_masks,pareto.st{jj}{ii}.(measures{jj}),[colors{ii} '-'])
      end

      for ii=1:params.n_r_cand
          % Pareto on one measure, evaluated on the other
          kk = setdiff([1 2],jj);
          plot(pareto.st{kk}{ii}.mean_n_masks,pareto.st{kk}{ii}.(measures{jj}),[colors{ii} '--'])
      end

      % Plot working point
      tmp = pareto.st{pareto_meas}{params.n_r_cand}.(measures{jj});
      plot(pareto.st{pareto_meas}{params.n_r_cand}.mean_n_masks(pareto_id),tmp(pareto_id),'k*')
      text(70,0.85,{[num2str(pareto.st{pareto_meas}{params.n_r_cand}.mean_n_masks(pareto_id)) ' c/i'],...
                    ['J=' num2str(tmp(pareto_id))]})

      % Make plot nicer
      title(strrep(titles{jj},'_','\_'))
      if (params.n_r_cand==3)
          legend({'Singletons same','Pairs same','Triplets same','Singletons diff','Pairs diff','Triplets diff'},4)
      elseif (params.n_r_cand==4)
          legend({'Singletons same','Pairs same','Triplets same','4-tuples','Singletons diff','Pairs diff','Triplets diff','4-tuples diff'},4)
      end
      grid minor
      grid on
      axis([10,x_limit,0.3,0.9])
      set(gca,'XScale','log')
  end

  %% Once you have found the desired point, save the parameters to file, i.e.,
  % the number of candidates from each hierarchy and for singletons, pairs, etc.
  sel_parameters  = pareto.pars{pareto_meas}{params.n_r_cand}(:,pareto_id);
  sampled_cands = pareto.base_stats{pareto_meas}.n_cands;
  total_cands = 0;
  clear n_cands;
  for ii=1:length(sel_parameters)
      if sel_parameters(ii)==0
          n_cands(ii) = 0; %#ok<*SAGROW>
      else
          n_cands(ii) = sampled_cands(sel_parameters(ii));
          total_cands = total_cands + pareto.base_stats{ii}.mean_n_masks(sel_parameters(ii));
      end
  end
  assert(abs(total_cands-pareto.st{pareto_meas}{params.n_r_cand}.mean_n_masks(pareto_id))<1e-5)
  n_cands = reshape(n_cands,length(params.hiers),params.n_r_cand)';
  n_cands

  save(params.files.pareto_point,'n_cands','pareto_id');
  disp(['Saved: ' params.files.pareto_point])
end
%% Write the curves to file (to put them in the paper)
% out_dir = '/Users/jpont/Publications/2014_CVPR/LaTeX/data/obj_cands/';
%
% % Write up to singletons, pairs, and triplets (w.r.t pareto_meas)
% write_jaccard_to_file(pareto.st{pareto_meas}{1},fullfile(out_dir,[params.gt_set_pareto '_' params.hiers_id '_pareto_singletons.txt']));
% write_jaccard_to_file(pareto.st{pareto_meas}{2},fullfile(out_dir,[params.gt_set_pareto '_' params.hiers_id '_pareto_pairs.txt']));
% write_jaccard_to_file(pareto.st{pareto_meas}{3},fullfile(out_dir,[params.gt_set_pareto '_' params.hiers_id '_pareto_triplets.txt']));
% write_jaccard_to_file(pareto.st{pareto_meas}{4},fullfile(out_dir,[params.gt_set_pareto '_' params.hiers_id '_pareto_4tuples.txt']));
%
% % Write selected pareto point
% % params = sf_mUCM_sub_params('sub1');
% load(params.files.pareto_point)
% fid = fopen(fullfile(out_dir,[params.gt_set_pareto '_' params.hiers_id '_pareto_selected_point.txt']),'w');
% fprintf(fid, 'ncands\tjac_class\tjac_instance\n');
% fprintf(fid, '%d\t%f\t%f\n', [pareto.st{pareto_meas}{4}.mean_n_masks(pareto_id); pareto.st{pareto_meas}{4}.jaccard_class(pareto_id); pareto.st{pareto_meas}{4}.jaccard_object(pareto_id)]);
% fclose(fid);

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function sparse_grid_cc_dataset ( dim_num, level_max )

%*****************************************************************************80
%
%% SPARSE_GRID_CC_DATASET is the main program for SPARSE_GRID_CC_DATASET.
%
%  Discussion:
%
%    This program computes a sparse grid quadrature rule based on 1D
%    Clenshaw-Curtis rules and writes it to a file..
%
%    The user specifies:
%    * the spatial dimension of the quadrature region,
%    * the level that defines the Smolyak grid.
%
%  License:
%
%    This software is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_CC_DATASET\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute the abscissas and weights of a quadrature rule\n' );
  fprintf ( 1, '  associated with a sparse grid derived from a Smolyak\n' );
  fprintf ( 1, '  construction based on 1D Clenshaw-Curtis rules.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Inputs to the program include:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    DIM_NUM, the spatial dimension.\n' );
  fprintf ( 1, '    (typically in the range of 2 to 10)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    LEVEL_MAX, the "level" of the sparse grid.\n' );
  fprintf ( 1, '    (typically in the range of 0, 1, 2, 3, ...\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Output from the program includes:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    A printed table of the abscissas and weights.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    * A set of 3 files that define the quadrature rule.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    "cc_d?_level?_r.txt", a file of the ranges.\n' );
  fprintf ( 1, '    "cc_d?_level?_w.txt", a file of the weights;\n' );
  fprintf ( 1, '    "cc_d?_level?_x.txt", a file of the abscissas;\n' );
%
%  Get the spatial dimension.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    dim_num = input ( '  Enter the value of DIM_NUM (1 or greater): ' );
  elseif ( ischar ( dim_num ) )
    dim_num = str2num ( dim_num );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension requested is = %d\n', dim_num );
%
%  Get the level.
%
  if ( nargin < 2 )
    fprintf ( 1, '\n' );
    level_max = input ( '  Enter the value of LEVEL_MAX (0 or greater): ' );
  elseif ( ischar ( level_max ) )
    level_max = str2num ( level_max );
  end

  level_min = max ( 0, level_max + 1 - dim_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEVEL_MIN = %d\n', level_min );
  fprintf ( 1, '  LEVEL_MAX = %d\n', level_max );
%
%  How many distinct points will there be?
%
  point_num = sparse_grid_cfn_size ( dim_num, level_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of distinct abscissas in the\n' );
  fprintf ( 1, '  quadrature rule is determined from the spatial\n' );
  fprintf ( 1, '  dimension DIM_NUM and the level LEVEL_MAX.\n' );
  fprintf ( 1, '  For the given input, this value will be = %d\n', point_num );

  r = zeros ( dim_num, 2 );
%
%  Compute the weights and points.
%
  r(1:dim_num,1) = -1.0;
  r(1:dim_num,2) = +1.0;

  [ w, x ] = sparse_grid_cc ( dim_num, level_max, point_num );

  r8mat_transpose_print_some ( dim_num, point_num, x, 1, 1, ...
    dim_num, 10, '  First 10 grid points:' );

  r8vec_print_some ( point_num, w, 1, 10, '  First 10 weights:' );

  weight_sum = sum ( w(1:point_num) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Weights sum to   %24.16f\n', weight_sum );
  fprintf ( 1, '  Correct value is %24.16f\n', 2.0^dim_num );
%
%  Write the rule to files.
%
  r_filename = sprintf ( 'cc_d%d_level%d_r.txt', dim_num, level_max );
  w_filename = sprintf ( 'cc_d%d_level%d_w.txt', dim_num, level_max );
  x_filename = sprintf ( 'cc_d%d_level%d_x.txt', dim_num, level_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Creating R file = "%s".\n', r_filename );

  r8mat_write ( r_filename, dim_num, 2, r );

  fprintf ( 1, '  Creating W file = "%s".\n', w_filename );

  r8mat_write ( w_filename, 1, point_num, w );

  fprintf ( 1, '  Creating X file = "%s".\n', x_filename );

  r8mat_write ( x_filename, dim_num, point_num, x );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_CC_DATASET:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function test_level = abscissa_level_closed_nd ( level_max, dim_num, ...
  test_num, test_val )

%*******************************************************************************
%
%% ABSCISSA_LEVEL_CLOSED_ND: first level at which given abscissa is generated.
%
%  Discussion:
%
%    We assume an underlying product grid.  In each dimension, this product
%    grid has order 2^LEVEL_MAX + 1.
%
%    We will say a sparse grid has total level LEVEL if each point in the
%    grid has a total level of LEVEL or less.
%
%    The "level" of a point is determined as the sum of the levels of the
%    point in each spatial dimension.
%
%    The level of a point in a single spatial dimension I is determined as
%    the level, between 0 and LEVEL_MAX, at which the point's I'th index
%    would have been generated.
%
%
%    This description is terse and perhaps unenlightening.  Keep in mind
%    that the product grid is the product of 1D grids,
%    that the 1D grids are built up by levels, having
%    orders (total number of points ) 1, 3, 5, 9, 17, 33 and so on,
%    and that these 1D grids are nested, so that each point in a 1D grid
%    has a first level at which it appears.
%
%    Our procedure for generating the points of a sparse grid, then, is
%    to choose a value LEVEL_MAX, to generate the full product grid,
%    but then only to keep those points on the full product grid whose
%    LEVEL is less than or equal to LEVEL_MAX.
%
%
%    Note that this routine is really just testing out the idea of
%    determining the level.  Our true desire is to be able to start
%    with a value LEVEL, and determine, in a straightforward manner,
%    all the points that are generated exactly at that level, or
%    all the points that are generated up to and including that level.
%
%    This allows us to generate the new points to be added to one sparse
%    grid to get the next, or to generate a particular sparse grid at once.
%
%  License:
%
%    This software is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer LEVEL_MAX, controls the size of the final sparse grid.
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer TEST_NUM, the number of points to be tested.
%
%    Input, integer TEST_VAL(DIM_NUM,TEST_NUM), the indices of the points
%    to be tested.  Normally, each index would be between 0 and 2^LEVEL_MAX.
%
%    Output, integer TEST_LEVEL(TEST_NUM), the value of LEVEL at which the
%    point would first be generated, assuming that a standard sequence of
%    nested grids is used.
%
  test_level = zeros ( test_num, 1 );

  if ( level_max == 0 )
    test_level(1:test_num) = 0;
    return
  end

  order = 2^level_max + 1;

  for j = 1 : test_num

    test_level(j) = index_to_level_closed ( dim_num, test_val(1:dim_num,j), ...
      order, level_max );

  end

  return
end
function value = cc_abscissa ( n, i )

%*******************************************************************************
%
%% CC_ABSCISSA returns the I-th abscissa of the Clenshaw Curtis rule.
%
%  Discussion:
%
%    The abscissas are numbered from left to right.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the rule.
%
%    Input, integer I, the index of the desired abscissa.  1 <= I <= N.
%
%    Output, real VALUE, the value of the I-th abscissa in the
%    rule of order N.
%
  if ( n < 1 )
    value = - Inf;
  elseif ( i < 1 || n < i )
    value = - Inf;
  elseif ( n == 1 )
    value = 0.0;
  elseif ( 2 * ( n - i ) == n - 1 )
    value = 0.0;
  else
    value = cos ( ( n - i ) * pi / ( n - 1 ) );
  end

  return
end
function w = cc_weights ( n )

%*******************************************************************************
%
%% CC_WEIGHTS computes Clenshaw Curtis weights.
%
%  License:
%
%    This software is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Charles Clenshaw, Alan Curtis,
%    A Method for Numerical Integration on an Automatic Computer,
%    Numerische Mathematik,
%    Volume 2, Number 1, December 1960, pages 197-205.
%
%  Parameters:
%
%    Input, integer N, the order of the rule.
%
%    Output, real W(N), the weights of the rule.
%
  w = zeros ( n, 1 );

  if ( n == 1 )
    w(1) = 2.0;
    return
  end

  theta = zeros ( n, 1 );

  for i = 1 : n
    theta(i) = ( i - 1 ) * pi / ( n - 1 );
  end

  for i = 1 : n

    w(i) = 1.0;

    for j = 1 : floor ( ( n - 1 ) / 2 )

      if ( 2 * j == ( n - 1 ) )
        b = 1.0;
      else
        b = 2.0;
      end

      w(i) = w(i) - b * cos ( 2.0 * j * theta(i) ) / ( 4 * j * j - 1 );

    end

  end

  w(1)     =       w(1)     / ( n - 1 );
  w(2:n-1) = 2.0 * w(2:n-1) / ( n - 1 );
  w(n)     =       w(n)     / ( n - 1 );

  return
end
function value = choose ( n, k )

%*******************************************************************************
%
%% CHOOSE computes the binomial coefficient C(N,K).
%
%  Discussion:
%
%    The value is calculated in such a way as to avoid overflow and
%    roundoff.  The calculation is done in integer arithmetic.
%
%    The formula used is:
%
%      C(N,K) = N! / ( K! * (N-K)! )
%
%  License:
%
%    This software is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    ML Wolfson, HV Wright,
%    Algorithm 160:
%    Combinatorial of M Things Taken N at a Time,
%    Communications of the ACM,
%    Volume 6, Number 4, April 1963, page 161.
%
%  Parameters:
%
%    Input, integer N, K, are the values of N and K.
%
%    Output, integer VALUE, the number of combinations of N
%    things taken K at a time.
%
  mn = min ( k, n - k );

  if ( mn < 0 )

    value = 0;

  elseif ( mn == 0 )

    value = 1;

  else

    mx = max ( k, n - k );
    value = mx + 1;

    for i = 2 : mn
      value = ( value * ( mx + i ) ) / i;
    end

  end

  return
end
function [ a, more, h, t ] = comp_next ( n, k, a, more, h, t )

%*******************************************************************************
%
%% COMP_NEXT computes the compositions of the integer N into K parts.
%
%  Discussion:
%
%    A composition of the integer N into K parts is an ordered sequence
%    of K nonnegative integers which sum to N.  The compositions (1,2,1)
%    and (1,1,2) are considered to be distinct.
%
%    The routine computes one composition on each call until there are no more.
%    For instance, one composition of 6 into 3 parts is
%    3+2+1, another would be 6+0+0.
%
%    On the first call to this routine, set MORE = FALSE.  The routine
%    will compute the first element in the sequence of compositions, and
%    return it, as well as setting MORE = TRUE.  If more compositions
%    are desired, call again, and again.  Each time, the routine will
%    return with a new composition.
%
%    However, when the LAST composition in the sequence is computed
%    and returned, the routine will reset MORE to FALSE, signaling that
%    the end of the sequence has been reached.
%
%    This routine originally used a SAVE statement to maintain the
%    variables H and T.  I have decided (based on an wasting an
%    entire morning trying to track down a problem) that it is safer
%    to pass these variables as arguments, even though the user should
%    never alter them.  This allows this routine to safely shuffle
%    between several ongoing calculations.
%
%    There are 28 compositions of 6 into three parts.  This routine will
%    produce those compositions in the following order:
%
%     I         A
%     -     ---------
%     1     6   0   0
%     2     5   1   0
%     3     4   2   0
%     4     3   3   0
%     5     2   4   0
%     6     1   5   0
%     7     0   6   0
%     8     5   0   1
%     9     4   1   1
%    10     3   2   1
%    11     2   3   1
%    12     1   4   1
%    13     0   5   1
%    14     4   0   2
%    15     3   1   2
%    16     2   2   2
%    17     1   3   2
%    18     0   4   2
%    19     3   0   3
%    20     2   1   3
%    21     1   2   3
%    22     0   3   3
%    23     2   0   4
%    24     1   1   4
%    25     0   2   4
%    26     1   0   5
%    27     0   1   5
%    28     0   0   6
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2008
%
%  Author:
%
%    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Albert Nijenhuis, Herbert Wilf,
%    Combinatorial Algorithms for Computers and Calculators,
%    Second Edition,
%    Academic Press, 1978,
%    ISBN: 0-12-519260-6,
%    LC: QA164.N54.
%
%  Parameters:
%
%    Input, integer N, the integer whose compositions are desired.
%
%    Input, integer K, the number of parts in the composition.
%
%    Input, integer A(K), the previous composition.  On the first call,
%    with MORE = FALSE, set A = [].  Thereafter, A should be the
%    value of A output from the previous call.
%
%    Input, logical MORE.  The input value of MORE on the first
%    call should be FALSE, which tells the program to initialize.
%    On subsequent calls, MORE should be TRUE, or simply the
%    output value of MORE from the previous call.
%
%    Input, integer H, T, two internal parameters needed for the
%    computation.  The user may need to initialize these before the
%    very first call, but these initial values are not important.
%    The user should not alter these parameters once the computation
%    begins.
%
%    Output, integer A(K), the next composition.
%
%    Output, logical MORE, will be TRUE unless the composition
%    that is being returned is the final one in the sequence.
%
%    Output, integer H, T, the updated values of the two internal
%    parameters.
%
  if ( ~more )

    t = n;
    h = 0;
    a(1) = n;
    a(2:k) = 0;

  else

    if ( 1 < t )
      h = 0;
    end

    h = h + 1;
    t = a(h);
    a(h) = 0;
    a(1) = t - 1;
    a(h+1) = a(h+1) + 1;

  end

  more = ( a(k) ~= n );

  return
end
function value = i4_modp ( i, j )

%*******************************************************************************
%
%% I4_MODP returns the nonnegative remainder of I4 division.
%
%  Discussion:
%
%    If
%      NREM = I4_MODP ( I, J )
%      NMULT = ( I - NREM ) / J
%    then
%      I = J * NMULT + NREM
%    where NREM is always nonnegative.
%
%    The MOD function computes a result with the same sign as the
%    quantity being divided.  Thus, suppose you had an angle A,
%    and you wanted to ensure that it was between 0 and 360.
%    Then mod(A,360) would do, if A was positive, but if A
%    was negative, your result would be between -360 and 0.
%
%    On the other hand, I4_MODP(A,360) is between 0 and 360, always.
%
%  Example:
%
%        I     J     MOD  I4_MODP    Factorization
%
%      107    50       7       7    107 =  2 *  50 + 7
%      107   -50       7       7    107 = -2 * -50 + 7
%     -107    50      -7      43   -107 = -3 *  50 + 43
%     -107   -50      -7      43   -107 =  3 * -50 + 43
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 1999
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the number to be divided.
%
%    Input, integer J, the number that divides I.
%
%    Output, integer VALUE, the nonnegative remainder when I is
%    divided by J.
%
  if ( j == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_MODP - Fatal error!\n' );
    fprintf ( 1, '  Illegal divisor J = %d\n', j );
    error ( 'I4_MODP - Fatal error!' );
  end

  value = mod ( i, j );

  if ( value < 0 )
    value = value + abs ( j );
  end

  return
end
function value = index_to_level_closed ( dim_num, t, order, level_max )

%*******************************************************************************
%
%% INDEX_TO_LEVEL_CLOSED determines the level of a point given its index.
%
%  License:
%
%    This software is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer T(DIM_NUM), the grid indices of a point in a 1D closed rule.
%    0 <= T(I) <= ORDER.
%
%    Input, integer ORDER, the order of the rule.
%
%    Input, integer LEVEL_MAX, the level with respect to which the
%    index applies.
%
%    Output, integer VALUE, the first level on which
%    the point associated with the given index will appear.
%
  value = 0;

  for dim = 1 : dim_num

    s = i4_modp ( t(dim), order );

    if ( s == 0 )

      level = 0;

    else

      level = level_max;

      while ( mod ( s, 2 ) == 0 )
        s = floor ( s / 2 );
        level = level - 1;
      end

    end

    if ( level == 0 )
      level = 1;
    elseif ( level == 1 )
      level = 0;
    end

    value = value + level;

  end

  return
end
function order = level_to_order_closed ( dim_num, level )

%*******************************************************************************
%
%% LEVEL_TO_ORDER_CLOSED converts a level to an order for closed rules.
%
%  Discussion:
%
%    Sparse grids can naturally be nested.  A natural scheme is to use
%    a series of one-dimensional rules arranged in a series of "levels"
%    whose order roughly doubles with each step.
%
%    The arrangement described here works naturally for the Clenshaw Curtis
%    and Newton Cotes closed rules.
%
%    The following table shows how the growth will occur:
%
%    Level    Order
%
%    0          1
%    1          3 =  2 + 1
%    2          5 =  4 + 1
%    3          9 =  8 + 1
%    4         17 = 16 + 1
%    5         33 = 32 + 1
%
%    For the Clenshaw Curtis and Newton Cotes Closed rules, the point growth
%    is nested.  If we have ORDER points on a particular LEVEL, the next
%    level includes all these old points, plus ORDER-1 new points, formed
%    in the gaps between successive pairs of old points.
%
%    Level    Order = New + Old
%
%    0          1   =  1  +  0
%    1          3   =  2  +  1
%    2          5   =  2  +  3
%    3          9   =  4  +  5
%    4         17   =  8  +  9
%    5         33   = 16  + 17
%
%    In this routine, we assume that a vector of levels is given,
%    and the corresponding orders are desired.
%
%  License:
%
%    This software is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL(DIM_NUM), the nesting level.
%
%    Output, integer ORDER(DIM_NUM), the order (number of points) of the rule.
%
  order = zeros ( dim_num, 1 );

  for dim = 1 : dim_num

    if ( level(dim) < 0 )
      order(dim) = -1;
    elseif ( level(dim) == 0 )
      order(dim) = 1;
    else
      order(dim) = ( 2^level(dim) ) + 1;
    end

  end

  return
end
function indx = multigrid_index0 ( dim_num, order_1d, order_nd )

%*******************************************************************************
%
%% MULTIGRID_INDEX0 returns an indexed multidimensional grid.
%
%  Discussion:
%
%    For dimension DIM, the second index of INDX may vary from
%    0 to ORDER_1D(DIM)-1.
%
%  License:
%
%    This software is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension of the points.
%
%    Input, integer ORDER_1D(DIM_NUM), the order of the
%    rule in each dimension.
%
%    Input, integer ORDER_ND, the product of the entries of ORDER_1D.
%
%    Output, integer INDX(DIM_NUM,ORDER_ND), the indices of the points in
%    the grid.  The second dimension of this array is equal to the
%    product of the entries of ORDER_1D.
%
  p = 0;
  indx = zeros(dim_num,order_nd);

  a = [];
  more = 0;

  while ( 1 )

    [ a, more ] = vec_colex_next2 ( dim_num, order_1d, a, more );

    if ( ~more )
      break
    end

    p = p + 1;

    indx(1:dim_num,p) = a(1:dim_num)';

  end

  return
end
function grid_index = multigrid_scale_closed ( dim_num, order_nd, level_max, ...
  level_1d, grid_index )

%*******************************************************************************
%
%% MULTIGRID_SCALE_CLOSED renumbers a grid as a subgrid on a higher level.
%
%  Discussion:
%
%    This routine takes a grid associated with a given value of
%    LEVEL, and multiplies all the indices by a power of 2, so that
%    the indices reflect the position of the same points, but in
%    a grid of level LEVEL_MAX.
%
%  License:
%
%    This software is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer ORDER_ND, the number of points in the grid.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Input, integer LEVEL_1D(DIM_NUM), the level in each dimension.
%
%    Input, integer GRID_INDEX(DIM_NUM,POINT_NUM), the index
%    values for each grid point, based in the level for which the grid
%    was generated.
%
%    Output, integer GRID_INDEX(DIM_NUM,POINT_NUM), the index
%    values for each grid point, appropriate for the grid as a subgrid
%    of a grid of level LEVEL_MAX.
%
  for dim = 1 : dim_num

    if ( level_1d(dim) == 0 )

      if ( 0 == level_max )
        order_max = 1;
      else
        order_max = 2^level_max + 1;
      end

      grid_index(dim,1:order_nd) = floor ( ( order_max - 1 ) / 2 );

    else

      factor = 2^( level_max - level_1d(dim) );

      grid_index(dim,1:order_nd) = grid_index(dim,1:order_nd) * factor;

    end

  end

  return
end
function w_nd = product_weights_cc ( dim_num, order_1d, order_nd )

%*******************************************************************************
%
%% PRODUCT_WEIGHTS_CC computes weights for a Clenshaw Curtis product rule.
%
%  Discussion:
%
%    This routine computes the weights for a quadrature rule which is
%    a product of 1D closed rules of varying order.
%
%  License:
%
%    This software is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer ORDER_1D(DIM_NUM), the order of the 1D rules.
%
%    Input, integer ORDER_ND, the order of the product rule.
%
%    Output, real W_ND(DIM_NUM,ORDER_ND), the product rule weights.
%
  w_nd(1:order_nd) = 1.0;

  for dim = 1 : dim_num

    w_1d = cc_weights ( order_1d(dim) );

    w_nd = r8vec_direct_product2 ( dim, order_1d(dim), w_1d, dim_num, ...
      order_nd, w_nd );

  end

  return
end
function r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*******************************************************************************
%
%% R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, an optional title.
%
  incx = 5;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d       ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%12f', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function w = r8vec_direct_product2 ( factor_index, factor_order, ...
  factor_value, factor_num, point_num, w )

%*******************************************************************************
%
%% R8VEC_DIRECT_PRODUCT2 creates a direct product of R8VEC's.
%
%  Discussion:
%
%    To explain what is going on here, suppose we had to construct
%    a multidimensional quadrature rule as the product of K rules
%    for 1D quadrature.
%
%    The product rule will be represented as a list of points and weights.
%
%    The J-th item in the product rule will be associated with
%      item J1 of 1D rule 1,
%      item J2 of 1D rule 2,
%      ...,
%      item JK of 1D rule K.
%
%    In particular,
%      X(J) = ( X(1,J1), X(2,J2), ..., X(K,JK))
%    and
%      W(J) = W(1,J1) * W(2,J2) * ... * W(K,JK)
%
%    So we can construct the quadrature rule if we can properly
%    distribute the information in the 1D quadrature rules.
%
%    This routine carries out that task for the weights W.
%
%    Another way to do this would be to compute, one by one, the
%    set of all possible indices (J1,J2,...,JK), and then index
%    the appropriate information.  An advantage of the method shown
%    here is that you can process the K-th set of information and
%    then discard it.
%
%  Example:
%
%    Rule 1:
%      Order = 4
%      W(1:4) = ( 2, 3, 5, 7 )
%
%    Rule 2:
%      Order = 3
%      W(1:3) = ( 11, 13, 17 )
%
%    Rule 3:
%      Order = 2
%      W(1:2) = ( 19, 23 )
%
%    Product Rule:
%      Order = 24
%      W(1:24) =
%        ( 2 * 11 * 19 )
%        ( 3 * 11 * 19 )
%        ( 4 * 11 * 19 )
%        ( 7 * 11 * 19 )
%        ( 2 * 13 * 19 )
%        ( 3 * 13 * 19 )
%        ( 5 * 13 * 19 )
%        ( 7 * 13 * 19 )
%        ( 2 * 17 * 19 )
%        ( 3 * 17 * 19 )
%        ( 5 * 17 * 19 )
%        ( 7 * 17 * 19 )
%        ( 2 * 11 * 23 )
%        ( 3 * 11 * 23 )
%        ( 5 * 11 * 23 )
%        ( 7 * 11 * 23 )
%        ( 2 * 13 * 23 )
%        ( 3 * 13 * 23 )
%        ( 5 * 13 * 23 )
%        ( 7 * 13 * 23 )
%        ( 2 * 17 * 23 )
%        ( 3 * 17 * 23 )
%        ( 5 * 17 * 23 )
%        ( 7 * 17 * 23 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer FACTOR_INDEX, the index of the factor being processed.
%    The first factor processed must be factor 1.
%
%    Input, integer FACTOR_ORDER, the order of the factor.
%
%    Input, real FACTOR_VALUE(FACTOR_ORDER), the factor values for
%    factor FACTOR_INDEX.
%
%    Input, integer FACTOR_NUM, the number of factors.
%
%    Input, integer POINT_NUM, the number of elements in the direct product.
%
%    Input/output, real W(POINT_NUM), the elements of the
%    direct product, updated by the latest factor.
%
%  Local Parameters:
%
%    Local, integer START, the first location of a block of values to set.
%
%    Local, integer CONTIG, the number of consecutive values to set.
%
%    Local, integer SKIP, the distance from the current value of START
%    to the next location of a block of values to set.
%
%    Local, integer REP, the number of blocks of values to set.
%
  persistent contig;
  persistent rep;
  persistent skip;

  if ( factor_index == 1 )
    contig = 1;
    skip = 1;
    rep = point_num;
    w(1:point_num) = 1.0;
  end

  rep = rep / factor_order;
  skip = skip * factor_order;

  for j = 1 : factor_order

    start = 1 + ( j - 1 ) * contig;

    for k = 1 : rep
      w(start:start+contig-1) = w(start:start+contig-1) * factor_value(j);
      start = start + skip;
    end

  end

  contig = contig * factor_order;

  return
end
function r8vec_print_some ( n, a, i_lo, i_hi, title )

%*******************************************************************************
%
%% R8VEC_PRINT_SOME prints "some" of an R8VEC.
%
%  Discussion:
%
%    An R8VEC is a vector of R8 values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%
%    Input, real A(N), the vector to be printed.
%
%    Input, integer MAX_PRINT, the maximum number of lines to print.
%
%    Input, string TITLE, an optional title.
%
  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  fprintf ( 1, '\n' );

  for i = max ( 1, i_lo ) : min ( n, i_hi )
    fprintf ( 1, '  %8d  %12f\n', i, a(i) );
  end

  return
end
function len = s_len_trim ( s )

%*******************************************************************************
%
%% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LEN, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
 function [ grid_weight, grid_point ] = sparse_grid_cc ( dim_num, ...
   level_max, point_num  )

%*******************************************************************************
%
%% SPARSE_GRID_CC computes a sparse grid of Clenshaw Curtis points.
%
%  Discussion:
%
%    This program computes a quadrature rule and writes it to a file.
%
%    The quadrature rule is associated with a sparse grid derived from
%    a Smolyak construction using a closed 1D quadrature rule.
%
%    The user specifies:
%    * the spatial dimension of the quadrature region,
%    * the level that defines the Smolyak grid.
%    * the closed 1D quadrature rule (Clenshaw-Curtis or Newton-Cotes Closed).
%
%  License:
%
%    This software is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, controls the size of the final sparse grid.
%
%    Input, integer POINT_NUM, the number of points in the grid, as determined
%    by SPARSE_GRID_CC_SIZE.
%
%    Output, real GRID_WEIGHTS(POINT_NUM), the weights.
%
%    Output, real GRID_POINTS(DIM_NUM,POINT_NUM), the points.
%

%
%  Determine the index vector, relative to the full product grid,
%  that identifies the points in the sparse grid.
%
  grid_index = sparse_grid_cc_index ( dim_num, level_max, point_num );
%
%  Compute the physical coordinates of the abscissas.
%
  if ( 0 == level_max )
    order_max = 1;
  else
    order_max = 2^level_max + 1;
  end

  grid_point = zeros ( dim_num, point_num );

  for point = 1 : point_num
    for dim = 1 : dim_num
      grid_point(dim,point) = ...
        cc_abscissa ( order_max, grid_index(dim,point) + 1 );
    end
  end
%
%  Gather the weights.
%
  grid_weight = sparse_grid_cc_weights ( dim_num, level_max, point_num, ...
    grid_index );

  return
end
function grid_index = sparse_grid_cc_index ( dim_num, level_max, point_num )

%*******************************************************************************
%
%% SPARSE_GRID_CC_INDEX indexes the points forming a sparse grid.
%
%  Discussion:
%
%    The points forming the sparse grid are guaranteed to be a subset
%    of a certain product grid.  The product grid is formed by DIM_NUM
%    copies of a 1D rule of fixed order.  The orders of the 1D rule,
%    (called ORDER_1D) and the order of the product grid, (called ORDER)
%    are determined from the value LEVEL_MAX.
%
%    Thus, any point in the product grid can be identified by its grid index,
%    a set of DIM_NUM indices, each between 1 and ORDER_1D.
%
%    This routine creates the GRID_INDEX array, listing (uniquely) the
%    points of the sparse grid.
%
%    An assumption has been made that the 1D rule is closed (includes
%    the interval endpoints) and nested (points that are part of a rule
%    of a given level will be part of every rule of higher level).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Input, integer POINT_NUM, the total number of points in the grids.
%
%    Output, integer GRID_INDEX(DIM_NUM,POINT_NUM), a list of point indices,
%    representing a subset of the product grid of level LEVEL_MAX,
%    representing (exactly once) each point that will show up in a
%    sparse grid of level LEVEL_MAX.
%
  grid_index = zeros ( dim_num, point_num );
%
%  The outer loop generates LEVELs from 0 to LEVEL_MAX.
%
  point_num2 = 0;

  for level = 0 : level_max
%
%  The middle loop generates the next partition LEVEL_1D(1:DIM_NUM)
%  that adds up to LEVEL.
%
    level_1d = [];
    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );
%
%  Transform each 1D level to a corresponding 1D order.
%
      order_1d = level_to_order_closed ( dim_num, level_1d );
%
%  The product of the 1D orders gives us the number of points in this grid.
%
      order_nd = prod ( order_1d(1:dim_num) );
%
%  The inner (hidden) loop generates all points corresponding to given grid.
%
      grid_index2 = multigrid_index0 ( dim_num, order_1d, order_nd );
%
%  Adjust these grid indices to reflect LEVEL_MAX.
%
      grid_index2 = multigrid_scale_closed ( dim_num, order_nd, level_max, level_1d, ...
        grid_index2 );
%
%  Determine the first level of appearance of each of the points.
%
      grid_level = abscissa_level_closed_nd ( level_max, dim_num, order_nd, ....
        grid_index2 );
%
%  Only keep those points which first appear on this level.
%
      for point = 1 : order_nd

        if ( grid_level(point) == level )

          point_num2 = point_num2 + 1;

          grid_index(1:dim_num,point_num2) = grid_index2(1:dim_num,point);

        end

      end

      if ( ~more )
        break
      end

    end

  end

  return
end
function point_num = sparse_grid_cfn_size ( dim_num, level_max )

%*****************************************************************************80
%
%% SPARSE_GRID_CFN_SIZE sizes a sparse grid using Closed Fully Nested rules.
%
%  Discussion:
%
%    The grid is defined as the sum of the product rules whose LEVEL
%    satisfies:
%
%      0 <= LEVEL <= LEVEL_MAX.
%
%    This calculation is much faster than a previous method.  It simply
%    computes the number of new points that are added at each level in the
%    1D rule, and then counts the new points at a given DIM_NUM dimensional
%    level vector as the product of the new points added in each dimension.
%
%    This approach will work for nested families, and may be extensible
%    to other families, and to mixed rules.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 December 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Output, integer POINT_NUM, the total number of unique
%    points in the grids.
%

%
%  Special case.
%
  if ( level_max < 0 )
    point_num = 0;
    return
  end

  if ( level_max == 0 )
    point_num = 1;
    return
  end
%
%  Construct the vector that counts the new points in the 1D rule.
%
  new_1d = zeros ( level_max+1, 1 );

  new_1d(0+1) = 1;
  new_1d(1+1) = 2;

  j = 1;
  for i = 2 : level_max
    j = j * 2;
    new_1d(i+1) = j;
  end
%
%  Count the number of points by counting the number of new points
%  associated with each level vector.
%
  level_1d = zeros ( dim_num, 1 );

  point_num = 0;

  for level = 0 : level_max

    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );

      point_num = point_num + prod ( new_1d(level_1d(1:dim_num)+1) );

      if ( ~more )
        break
      end

    end

  end

  return
end
function grid_weight = sparse_grid_cc_weights ( dim_num, level_max, point_num, ...
  grid_index )

%*******************************************************************************
%
%% SPARSE_GRID_CC_WEIGHTS gathers the weights.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Input, integer POINT_NUM, the total number of points in the grids.
%
%    Input, integer GRID_INDEX(DIM_NUM,POINT_NUM), a list of point indices,
%    representing a subset of the product grid of level LEVEL_MAX,
%    representing (exactly once) each point that will show up in a
%    sparse grid of level LEVEL_MAX.
%
%    Output, real GRID_WEIGHT(POINT_NUM), the weights
%    associated with the sparse grid points.
%
  if ( level_max == 0 )
    grid_weight(1:point_num) = 2.0^dim_num;
    return
  end

  grid_weight(1:point_num) = 0.0;

  level_min = max ( 0, level_max + 1 - dim_num );

  for level = level_min : level_max
%
%  The middle loop generates the next partition LEVEL_1D(1:DIM_NUM)
%  that adds up to LEVEL.
%
    level_1d = [];
    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );
%
%  Transform each 1D level to a corresponding 1D order.
%
      order_1d = level_to_order_closed ( dim_num, level_1d );
%
%  The product of the 1D orders gives us the number of points in this grid.
%
      order_nd = prod ( order_1d(1:dim_num) );
%
%  Generate the indices of the points corresponding to the grid.
%
      grid_index2 = multigrid_index0 ( dim_num, order_1d, order_nd );
%
%  Compute the weights for this grid.
%
      grid_weight2 = product_weights_cc ( dim_num, order_1d, order_nd );
%
%  Adjust the grid indices to reflect LEVEL_MAX.
%
      grid_index2 = multigrid_scale_closed ( dim_num, order_nd, level_max, ...
        level_1d, grid_index2 );
%
%  Now determine the coefficient.
%
      coeff = (-1)^( level_max - level ) ...
        * choose ( dim_num - 1, level_max - level );

      for point2 = 1 : order_nd

        for point = 1 : point_num

          if ( all ( ...
            grid_index2(1:dim_num,point2) == grid_index(1:dim_num,point) ...
          ) )
            grid_weight(point) = grid_weight(point) ...
              + coeff * grid_weight2(point2);
            break
          end

        end

      end

      if ( ~more )
        break
      end

    end

  end

  return
end
function timestamp ( )

%*******************************************************************************
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
function [ a, more ] = vec_colex_next2 ( dim_num, base, a, more )

%*******************************************************************************
%
%% VEC_COLEX_NEXT2 generates vectors in colex order.
%
%  Discussion:
%
%    The vectors are produced in colexical order, starting with
%    (0,0,...,0),
%    (1,0,...,0),
%    ...
%    (BASE(1)-1,BASE(2)-1,...,BASE(DIM_NUM)-1).
%
%  Example:
%
%    DIM_NUM = 2,
%    BASE = [ 3, 3]
%
%    0   0
%    1   0
%    2   0
%    0   1
%    1   1
%    2   1
%    0   2
%    1   2
%    2   2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dennis Stanton, Dennis White,
%    Constructive Combinatorics,
%    Springer, 1986,
%    ISBN: 0387963472,
%    LC: QA164.S79.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer BASE(DIM_NUM), the base to be used in each dimension.
%
%    Input, integer A(DIM_NUM), except on the first call, this should
%    be the output value of A on the last call.
%
%    Input, logical MORE, should be FALSE on the first call, and
%    thereafter should be the output value of MORE from the previous call.
%
%    Output, integer A(DIM_NUM), the next vector.
%
%    Output, logical MORE, is TRUE if another vector was computed.
%    If MORE is FALSE on return, then ignore the output value A, and
%    stop calling the routine.
%
  if ( ~more )

    a(1:dim_num) = 0;
    more = 1;

  else

    for i = 1 : dim_num

      a(i) = a(i) + 1;

      if ( a(i) < base(i) )
        return
      end

      a(i) = 0;

    end

    more = 0;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [hAx, prim, varargin] = parseDrawInput(prim, valFun, type, defOpts, varargin)
%PARSEDRAWINPUT Parse input arguments for drawing functions: draw*.
%
%   INPUT:
%       PRIM: The primitive object: line, plane, ...
%       VALFUN: An anonymous Function to validate PRIM
%       TYPE: The drawing type of PRIM: 'line', 'patch', ...
%       DEFOPTS: The default drawing options for PRIM as struct
%
%   OUTPUT:
%       HAX: The current or specified axes for drawing
%       PRIM: validated PRIM

% ------
% Author: oqilipo
% E-mail: N/A
% Created: 2017-10-13, using R2017b
% Copyright 2017-2022

% Check if first input argument is an axes handle
if isAxisHandle(prim)
    hAx = prim;
    prim = varargin{1};
    varargin(1) = [];
else
    hAx = gca;
end

% Check if the primitive is valid
p=inputParser;
addRequired(p,'prim',valFun)
parse(p,prim)

% parse input arguments if there are any
if ~isempty(varargin)
    if length(varargin) == 1
        if isstruct(varargin{1})
            % if options are specified as struct, need to convert to
            % parameter name-value pairs
            varargin = [fieldnames(varargin{1}) struct2cell(varargin{1})]';
            varargin = varargin(:)';
        else
            % if option is a single argument, assume it is the color
            switch type
                case 'line'
                    varargin = {'Color', varargin{1}};
                case 'patch'
                    varargin = {'FaceColor', varargin{1}};
            end
        end
    end
else
    % If no arguments are given, use the default options
    varargin = [fieldnames(defOpts) struct2cell(defOpts)]';
    varargin = varargin(:)';
end

end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function coinWrite(prob,filename,type)
%COINWRITE  Write a Mathematical File From Matlab Values
%
% coinWrite(prob,filename,type) writes the problem prob to the file specified
% by filename and problem type specified by type. If you do not specify a
% file extension it will default to the type specified.
%
%   Available File Types:
%       - MPS  [.mps]
%       - QPS  [.qps]
%       - LP   [.lp]  (Appears to be CPLEX LP version - Simplified)
%
% You may specify a full path to the file, or if you specify a filename
% only, it will be written to the current MATLAB directory.
%
% The routines underneath use COIN-OR utilities for File IO. See
% attached EPL License.

%   Copyright (C) 2011 Jonathan Currie (I2C2)

%Quick Checks
if(~exist('type','var') || isempty(type))
    %See if contained within file name
    dots = strfind(filename,'.');
    if(isempty(dots))
        error('You must supply the file type to this function');
    else %hopefully got a good ext, will check below
        type = filename(dots(end)+1:end);
    end
end
if(~ischar(filename))
    error('Filename must be a char array!');
end

%Problem Checks
if(~isempty(prob.fun) || ~isempty(prob.nlcon))
    error('You cannot write Nonlinear problems using this interface');
end
if(isa(prob.f,'function_handle') || isa(prob.H,'function_handle'))
    error('f and H must be vectors / matrices!');
end
if(~isempty(prob.Q))
    error('You cannot write quadratic constraints using this interface');
end
if(prob.sense == -1)
    error('You cannot explicity write maximization problems using this interface'); %fix up later
end

%Determine file type
switch(lower(type))
    case 'mps'
        filetype = 'mps';
    case 'qps'
        filetype = 'qps';
    case 'lp'
        filetype = 'lp';
    otherwise
        error('Unknown file type %s',type);
end

%Check problem type vs file type
if(isfield(prob,'type') && ~isempty(prob.type))
    switch(lower(prob.type))
        case {'lp','bilp','milp'}
            if(strcmpi(filetype,'qps'))
                error('LP/BILP/MILPs should be written to MPS or LP files');
            end
        case {'qp','miqp'}
            if(~strcmpi(filetype,'qps'))
                error('QP/MIQPs should be written to QPS files');
            end
        otherwise
            error('Problem type %s is not supported by this interface',upper(prob.type));
    end
end

% Make sure we have file extension
if(isempty(strfind(filename,['.' filetype])))
    filename = [filename '.' filetype];
end

% If filename includes absolute path, we can skip which
if(~isempty(strfind(filename,':')))
    p = filename;
else %Locate the full path to the file
    p = [cd filesep filename];
end

%Convert to coin problem
if(~isempty(prob.b) || ~isempty(prob.beq))
    [prob.A,prob.rl,prob.ru] = gen2row(prob.A,prob.b,prob.Aeq,prob.beq);
end
if(isempty(prob.A))
    prob.A = spalloc(0,length(prob.f),0); %error thrown internally if not ok
end

%Ensure sparsity
prob.A = sparse(prob.A);
prob.H = sparse(prob.H);

%Ensure Bounds
if(isempty(prob.lb))
    prob.lb = -Inf(size(prob.f));
end
if(isempty(prob.ub))
    prob.ub = Inf(size(prob.f));
end

%Setup Integer Vars
if(isstruct(prob.int))
    prob.int = prob.int.str;
end
if(isempty(prob.int))
    prob.int = zeros(size(prob.f));
elseif(ischar(prob.int)) %char string
    prob.int = lower(prob.int);
    i = zeros(size(prob.f));
    int = prob.int == 'i';
    bin = prob.int == 'b';
    i(int) = 1;
    i(bin) = 1;
    prob.lb(bin) = 0;
    prob.ub(bin) = 1;
    prob.int = i;
elseif(length(prob.int) ~= length(prob.f) || any(prob.int > 1)) %find indicies (might be wrong though!)
    i = zeros(size(prob.f));
    i(prob.int) = 1;
    prob.int = i;
end
prob.int = int8(prob.int);

%Setup SOS
if(~isempty(prob.sos))
    [r,c] = size(prob.sos.type);
    if(r > c)
        prob.sos_type = str2num(prob.sos.type); %#ok<ST2NM>
    else
        prob.sos_type = str2num(prob.sos.type'); %#ok<ST2NM>
    end
    if(length(prob.sostype) == 1)
        prob.sos_index = {prob.sos.index};
        prob.sos_weight = {prob.sos.weight};
    end
else
    prob.sos_type = [];
    prob.sos_index = [];
    prob.sos_weight = [];
end
%Setup Objc
if(~isempty(prob.objbias) && prob.objbias == 0)
    prob.objbias = 0;
end

%Call coin-or mex
coinW(prob,p,filetype);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus