text:
function [tasklist, taskfile] = LoadTaskFile_Reverb(taskFileRoot, dataset, datatype, distance, nCh, roomID)
switch dataset
    case 'train'
        datasetStr = 'tr';
    case 'dev'
        datasetStr = 'dt';
    case 'eval'
        datasetStr = 'et';
    otherwise
        fprintf('Error: unknown dataset: %s!\n', dataset);
        return;
end
if strcmpi(datatype, 'simu')
    datatypeStr = 'Sim';
else
    datatypeStr = 'Real';
end

chStr = {'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H'};

if strcmp(distance, 'cln')
    if strcmp(dataset, 'train')
        tasklist = {};
    else
        taskfile = [taskFileRoot '/' num2str(nCh) 'ch/' datatypeStr 'Data_' datasetStr '_for_' ...
            distance '_room' num2str(roomID)];
        tasklist = my_cat(taskfile)';
    end
else
    for i=1:nCh
        if strcmp(dataset, 'train')
            taskfile{i} = [taskFileRoot '/' num2str(nCh) 'ch/' datatypeStr 'Data_' datasetStr '_for_' ...
                num2str(nCh) 'ch_' chStr{i}];
        else
            taskfile{i} = [taskFileRoot '/' num2str(nCh) 'ch/' datatypeStr 'Data_' datasetStr '_for_' ...
                num2str(nCh) 'ch_' distance '_room' num2str(roomID) '_' chStr{i}];
        end
        tasklist(i,:) = my_cat(taskfile{i});
    end
end
end

author: singaxiong

repo: SignalGraph


text:
function n = norm(a)
n = sqrt(sum(a.^2));

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%SAD	Stereo matching using Zero Mean Sum of Absolute Differences metric
%
%	[DISP SCORES] = ZSAD(LEFT, RIGHT, LEFT_RIGHT, WIN_SIZE, DISP_RANGE)
%
%       DISP = Disparity map
%       SCORES = ZSAD scores at each disparity
%       LEFT, RIGHT = images
%       LEFT_RIGHT = 'l' if wrt left image, 'r' if wrt right image
%       WIN_SIZE = window size or [x_size y_size]
%       DISP_RANGE = [START_DISPARITY END_DISPARITY]
% Author: Jasmine E. Banks (jbanks@ieee.org)

% Copyright in this software is owned by CSIRO.  CSIRO grants permission to
% any individual or institution to use, copy, modify, and distribute this
% software, provided that:
%
% (a)     this copyright and permission notice appears in its entirety in or
% on (as the case may be) all copies of the software and supporting
% documentation;
%
% (b)     the authors of papers that describe software systems using this
% software package acknowledge such use by citing the paper as follows:
%
%     "Quantitative Evaluation of Matching Methods and Validity Measures for
%     Stereo Vision" by J. Banks and P. Corke, Int. J. Robotics Research,
%     Vol 20(7), 2001; and
%
% (c)     users of this software acknowledge and agree that:
%
%   (i) CSIRO makes no representations about the suitability of this software
%   for any purpose;
%
%   (ii) that the software is provided "as is" without express or implied
%   warranty; and
%
%   (iii) users of this software use the software entirely at their own risk.

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [result] = pf_bvalues(params)
% function [params] = pf_bvalues(params)
% --------------------------------------
% Calculation of b-values for the probabilistic forecast test.
% Input parameters:
%   params.mCatalog           Earthquake catalog
%   params.mPolygon           Polygon (defined by ex_selectgrid)
%   params.vX                 X-vector (defined by ex_selectgrid)
%   params.vY                 Y-vector (defined by ex_selectgrid)
%   params.vUsedNodes         Used nodes vX * vY defining the mPolygon (defined by ex_selectgrid)
%   params.bMap               Calculate a map (=1) or a cross-section (=0)
%   params.bNumber            Use constant number (=1) or constant radius (=0)
%   params.nNumberEvents      Number of earthquakes if bNumber == 1
%   params.fRadius            Radius of gridnode if bNumber == 0
%   params.nMinimumNumber     Minimum number of earthquakes per node
%   params.fSplitTime         Time at which the catalog will be split
%   params.bLearningPeriod    Fix duration of the learning period (1 = fix, 0 = use catalog from start)
%   params.fLearningPeriod    Duration of the learning period
%   params.bForecastPeriod    Fix duration of the forecasting period (1 = fix, 0 = use catalog till end)
%   params.fForecastPeriod    Duration of the forecasting period
%
% Output parameters:
%   Same as input parameters including
%   result.mValueGrid         Matrix of calculated b-values
%   result.vcsGridNames       Strings describing the values in mValueGrid
%
% Danijel Schorlemmer
% June 4, 2001

global bDebug;
if bDebug
  report_this_filefun(mfilename('fullpath'));
end

% Init result matrix
result = params;
result.mValueGrid = [];

for nNode_ = 1:length(params.mPolygon(:,1))
  % Create node catalog
  mNodeCatalog_ = params.mCatalog(params.caNodeIndices{nNode_}, :);
  % Calculate the b-values
  [fDeltaBValue, fBValueLearning, fBValueObserved, fStdDevLearning, fStdDevObserved] ...
    = pf_calcbvalues(mNodeCatalog_, params.fSplitTime, params.bLearningPeriod, params.fLearningPeriod, ...
    params.bForecastPeriod, params.fForecastPeriod, params.nMinimumNumber, params.nCalculateMC);
  % Store the results
  result.mValueGrid = [result.mValueGrid; fDeltaBValue fBValueLearning fBValueObserved fStdDevLearning fStdDevObserved];
end % for nNode
% Create the description-strings for the output-window
result.vcsGridNames = cellstr(char('b-value difference', 'b-value of learning period', 'b-value of observation period', ...
  'Std. dev. of learning period b-value ', 'Std. dev. of observation period b-value'));

author: CelsoReyes


text:
% checks the resolution limit of root-MUSIC
clear(); close all;

wavelength = 1; % normalized
d = wavelength / 2;
design = design_array_1d('ula', 10, d); % 10-element ULA
power_source = 1;
power_noise = 1;
snapshot_count = 100;

delta_theta = deg2rad(linspace(0.05, 5, 20)); % 0.1 deg - 10 deg
n_delta_theta = length(delta_theta);
n_repeat = 500;
res_prob = zeros(n_delta_theta, 1);
progressbar('reset', n_delta_theta);
for ii = 1:n_delta_theta
    n_resolved = 0;
    for rr = 1:n_repeat
        doas = [-delta_theta(ii)/2 delta_theta(ii)/2];
        [~, R] = snapshot_gen_sto(design, doas, wavelength, snapshot_count, power_noise, power_source);
        sp = rmusic_1d(R, 2, 2*pi*design.element_spacing/wavelength);
        % resolution condition:
        %   theta_1 \in (-delta_theta, 0)
        %   theta_2 \in (0, delta_theta)
        % This condition is pretty strict. It requires the estimated DOAs
        % are reasonably close to the correct ones.
        % You may change it to a milder condition.
        if sp.x_est(1) < 0 && sp.x_est(1) > -delta_theta(ii) ...
                && sp.x_est(2) > 0 && sp.x_est(2) < delta_theta(ii)
            n_resolved = n_resolved + 1;
        end
    end
    res_prob(ii) = n_resolved / n_repeat;
    progressbar('advance');
end
fprintf('\n');

plot(rad2deg(delta_theta), res_prob); grid on;
xlabel('Source separation (deg)'); ylabel('Resolution Probability');

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function pathname = homePageFullPath()
pathname = '';
currdir = pwd;
if ismac() || islinux()
    cd('~/');
    pathname = filesepStandard(pwd);
end
cd(currdir)

author: BUNPC


text:
%io_loadspec_bruk.m
%Chathura Kumaragamage, McGill University 2016.
%Jamie Near, McGill University 2016.
%
% USAGE:
% [out,ref]=io_loadspec_bruk(filename);
%
% DESCRIPTION:
% Reads in Bruker MRS data (fid.raw, fid.ref).
%
% op_loadspec_bruk outputs the data in structure format, with fields corresponding to time
% scale, fids, frequency scale, spectra, and header fields containing
% information about the acquisition.  The resulting matlab structure can be
% operated on by the other functions in this MRS toolbox.
%
% INPUTS:
% inDir   = Path to the scan number directory that contains the 'pdata' folder.
%
% OUTPUTS:
% out = Input dataset in FID-A structure format.
% ref = The Reference scan data (navigator echoes) in FID-A structure
%       format, if applicable.

function [out,ref]=io_loadspec_bruk(inDir)

%%%%%%%%chathu mod starts

%get the original number of Repetitions (all repatitions from scanner is generated as
%a 1D vector. Need to split before further processing
averages=1; %Because Bruker does the averaging online.
ADC_OFFSET=70;      %(offset between ADC on and ADC acquire, empirically determined)
method_fid=fopen([inDir '/method']);
line=fgets(method_fid);
index=findstr(line,'$PVM_DigNp=');
while isempty(index)
    line=fgets(method_fid);
    index=findstr(line,'$PVM_DigNp=');
end
equals_index=findstr(line,'=');
rawDataPoints=line(equals_index+1:end);
rawDataPoints=str2double(rawDataPoints);
fclose(method_fid);

%NOW LOAD IN THE RAW DATA.  FIRST TRY USING THE FID.RAW FILE.  IF THAT DOES
%NOT WORK, THEN USE THE REGULAR FID FID.
try
    fid_data=fread(fopen([inDir '/fid.raw']),'int');
    real_fid = fid_data(2:2:length(fid_data));
    imag_fid = fid_data(1:2:length(fid_data));
    fids_raw=real_fid-1i*imag_fid;

    method_fid=fopen([inDir '/method']);
    line=fgets(method_fid);
    index=findstr(line,'$PVM_NAverages=');
    while isempty(index)
        line=fgets(method_fid);
        index=findstr(line,'$PVM_NAverages=');
    end
    equals_index=findstr(line,'=');
    rawAverages=line(equals_index+1:end);
    rawAverages=str2double(rawAverages);
    fclose(method_fid);
    fids_raw=reshape(fids_raw,[],rawAverages);
    fids_trunc=fids_raw(ADC_OFFSET:end,:);
    %         fids_trunc=fids_trunc';
    fids=padarray(fids_trunc, [ADC_OFFSET-1,0],'post');
    specs=fftshift(ifft(fids,[],1),1);
    sz=size(specs);
    out.flags.averaged=0;
    %specify the dims
    dims.t=1;
    dims.coils=0;
    dims.averages=2;
    dims.subSpecs=0;
    dims.extras=0;
catch
    %NOW TRY USING FID FILE
    display('WARNING: /fid.raw not found. Using /fid ....')
    %specs=fread(fopen([inDir '/pdata/1/2dseq']),'int');
    fid_data=fread(fopen([inDir '/fid']),'int');
    real_fid = fid_data(2:2:length(fid_data));
    imag_fid = fid_data(1:2:length(fid_data));
    fids_raw=real_fid-1i*imag_fid;

    rawAverages=size(real_fid,1)./rawDataPoints;

    if mod(size(real_fid,1),rawDataPoints) ~=0
        display 'number of repetitions cannot be accurately found';
    end
    fids_raw=reshape(fids_raw,rawDataPoints,[]);

    fids_trunc=fids_raw(ADC_OFFSET:end,:);
    fids=padarray(fids_trunc, [ADC_OFFSET-1,0],'post');
    %convert back to time domain
    %if the length of Fids is odd, then you have to do a circshift of one to
    %make sure that you don't introduce a small frequency shift into the fids
    %vector.
    specs=fftshift(ifft(fids,[],1),1);

    %calculate the size;
    sz=size(specs);
    out.flags.averaged=1;
    %specify the dims
    dims.t=1;
    dims.coils=0;
    dims.averages=0;
    dims.subSpecs=0;
    dims.extras=0;
end

%NOW TRY LOADING IN THE REFERENCE SCAN DATA (IF IT EXISTS)
isRef=false;
if exist([inDir '/fid.ref'])
    isRef=true;
    ref_data=fread(fopen([inDir '/fid.ref']),'int');
    real_ref = ref_data(2:2:length(ref_data));
    imag_ref = ref_data(1:2:length(ref_data));
    fids_ref=real_ref-1i*imag_ref;

    fids_ref=reshape(fids_ref,[],rawAverages);
    fids_ref_trunc=fids_ref(ADC_OFFSET:end,:);
    %         fids_trunc=fids_trunc';
    reffids=padarray(fids_ref_trunc, [ADC_OFFSET-1,0],'post');
    refspecs=fftshift(ifft(reffids,[],1),1);
    sz=size(refspecs);
    out.flags.averaged=0;
    %specify the dims
    refdims.t=1;
    refdims.coils=0;
    refdims.averages=2;
    refdims.subSpecs=0;
    refdims.extras=0;
end

%%%%%Chathu mod ends

%Now get some of the relevent spectral parameters

%First get the spectral width
method_fid=fopen([inDir '/method']);
line=fgets(method_fid);
index=findstr(line,'$PVM_DigSw=');
while isempty(index)
    line=fgets(method_fid);
    index=findstr(line,'$PVM_DigSw=');
end
equals_index=findstr(line,'=');
spectralwidth=line(equals_index+1:end);
spectralwidth=str2double(spectralwidth);
fclose(method_fid);

%Now get the transmitter frequency
acqp_fid=fopen([inDir '/acqp']);
line=fgets(acqp_fid);
index=findstr(line,'$BF1=');
while isempty(index)
    line=fgets(acqp_fid);
    index=findstr(line,'$BF1=');
end
equals_index=findstr(line,'=');
txfrq=line(equals_index+1:end);
txfrq=str2double(txfrq);
txfrq=txfrq*1e6;
fclose(acqp_fid);

%B0
Bo=txfrq/42577000;

%Spectral width in PPM
spectralwidthppm=spectralwidth/(txfrq/1e6);

%Now get the TE
method_fid=fopen([inDir '/method']);
line=fgets(method_fid);
index=findstr(line,'$PVM_EchoTime=');
%%%%%CHATHU MOD for FIDS with Pulse acquire
loop_count=0;
TE_present=true;
while isempty(index) && TE_present
    line=fgets(method_fid);
    index=findstr(line,'$PVM_EchoTime=');
    loop_count=loop_count+1;
    if loop_count>10000
        TE_present=false;
    end
end

if TE_present
    equals_index=findstr(line,'=');
    te=line(equals_index+1:end);
    te=str2double(te);
else
    te=0;
end

%%%%%CHATHU MOD for FIDS with Pulse acquire END
fclose(method_fid);

%Now get the TR
method_fid=fopen([inDir '/method']);
line=fgets(method_fid);
index=findstr(line,'$PVM_RepetitionTime=');
while isempty(index)
    line=fgets(method_fid);
    index=findstr(line,'$PVM_RepetitionTime=');
end
equals_index=findstr(line,'=');
tr=line(equals_index+1:end);
tr=str2double(tr);
fclose(method_fid);

%Now get the sequence
method_fid=fopen([inDir '/method']);
line=fgets(method_fid);
index=findstr(line,'$Method=');
while isempty(index)
    line=fgets(method_fid);
    index=findstr(line,'$Method=');
end
equals_index=findstr(line,'=');
sequence=line(equals_index+1:end);
sequence=strtrim(sequence);
fclose(method_fid);

%Specify the number of subspecs.  For now, this will always be one.
subspecs=1;
rawSubspecs=1;

%calculate the ppm scale
ppm=[4.65+(spectralwidthppm/2):-spectralwidthppm/(length(specs)-1):4.65-(spectralwidthppm/2)];

%calculate the dwelltime:
dwelltime=1/spectralwidth;

%calculate the time scale
t=[0:dwelltime:(sz(1)-1)*dwelltime];

%FILLING IN DATA STRUCTURE FOR THE FID.RAW DATA
out.fids=fids;
out.specs=specs;
out.sz=sz;
out.ppm=ppm;
out.t=t;
out.spectralwidth=spectralwidth;
out.dwelltime=dwelltime;
out.txfrq=txfrq;
out.date=date;
out.dims=dims;
out.Bo=Bo;
out.averages=rawAverages;
out.rawAverages=rawAverages;
out.subspecs=subspecs;
out.rawSubspecs=rawSubspecs;
out.seq=sequence;
out.te=te;
out.tr=tr;
out.pointsToLeftshift=0;

%FILLING IN THE FLAGS FOR THE FID.RAW DATA
out.flags.writtentostruct=1;
out.flags.gotparams=1;
out.flags.leftshifted=0;
out.flags.filtered=0;
out.flags.zeropadded=0;
out.flags.freqcorrected=0;
out.flags.phasecorrected=0;

out.flags.addedrcvrs=1;
out.flags.subtracted=0;
out.flags.writtentotext=0;
out.flags.downsampled=0;
out.flags.avgNormalized=0;
if out.dims.subSpecs==0
    out.flags.isFourSteps=0;
else
    out.flags.isFourSteps=(out.sz(out.dims.subSpecs)==4);
end

if isRef
    %FILLING IN DATA STRUCTURE FOR THE FID.REF DATA
    ref.fids=reffids;
    ref.specs=refspecs;
    ref.sz=sz;
    ref.ppm=ppm;
    ref.t=t;
    ref.spectralwidth=spectralwidth;
    ref.dwelltime=dwelltime;
    ref.txfrq=txfrq;
    ref.date=date;
    ref.dims=dims;
    ref.Bo=Bo;
    ref.averages=rawAverages;
    ref.rawAverages=rawAverages;
    ref.subspecs=subspecs;
    ref.rawSubspecs=rawSubspecs;
    ref.seq=sequence;
    ref.te=te;
    ref.tr=tr;
    ref.pointsToLeftshift=0;

    %FILLING IN THE FLAGS FOR THE FID.REF DATA
    ref.flags.writtentostruct=1;
    ref.flags.gotparams=1;
    ref.flags.leftshifted=0;
    ref.flags.filtered=0;
    ref.flags.zeropadded=0;
    ref.flags.freqcorrected=0;
    ref.flags.phasecorrected=0;

    ref.flags.addedrcvrs=1;
    ref.flags.subtracted=0;
    ref.flags.writtentotext=0;
    ref.flags.downsampled=0;
    ref.flags.avgNormalized=0;
    if ref.dims.subSpecs==0
        ref.flags.isFourSteps=0;
    else
        ref.flags.isFourSteps=(ref.sz(ref.dims.subSpecs)==4);
    end
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function y = num2str(varargin)
% Avoids '-0' in result by the simple procedure of adding 0 to the input
% if it's numeric. ('-0' arises in Octave because of "negative 0". It turns
% out that adding "(positive) 0" to "negative 0" gives "(positive) 0".
if isnumeric(varargin{1})
    varargin{1} = varargin{1} + 0;
end
% Fixes alignment in certain cases
% http://stackoverflow.com/q/34483961/2586922
% https://savannah.gnu.org/bugs/?46770
x = varargin{1}; x = reshape(x, size(x,1),[]);
if nargin==1 || ischar(varargin{1}) || isnumeric(varargin{2}) || any(imag(x(:)))
    y = builtin('num2str', varargin{:});
else
    fmt = varargin{2}; y = sprintf([fmt '\n'], x.'); y = regexp(y, '\n', 'split'); y = y(1:end-1).';
    y = cellfun(@fliplr, y, 'uniformoutput', false); y = char(y); y = fliplr(y);
    y = reshape(y.',[],size(x,1)).'; y = strtrim(y);
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function [ T ] = computeTrans( imgs )
Thresh = 10;
confidence = 0.99;
inlierRatio = 0.3;
epsilon = 1.5;

nImgs = size(imgs, 4);

T = zeros(3, 3, nImgs);
T(:, :, 1) = eye(3);
[f2, d2] = getSIFTFeatures(imgs(:, :, :, 1), Thresh);
for i = 2 : nImgs
    f1 = f2;
    d1 = d2;
    [f2, d2] = getSIFTFeatures(imgs(:, :, :, i), Thresh);
    [matches, ~] = getMatches(f1, d1, f2, d2);
    [T(:, :, i),~] = RANSAC(confidence, inlierRatio, 1, matches, epsilon);
end
end

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [ r, code, w, m ] = s_to_format ( s )

%*****************************************************************************80
%
%% S_TO_FORMAT reads a FORTRAN format from a string.
%
%  Discussion:
%
%    This routine will read as many characters as possible until it reaches
%    the end of the string, or encounters a character which cannot be
%    part of the format.  This routine is limited in its ability to
%    recognize FORTRAN formats.  In particular, we are only expecting
%    a single format specification, and cannot handle extra features
%    such as 'ES' and 'EN' codes, '5X' spacing, and so on.
%
%    Legal input is:
%
%       0 nothing
%       1 blanks
%       2 optional '('
%       3 blanks
%       4 optional repeat factor R
%       5 blanks
%       6 CODE ( 'A', 'B', 'E', 'F', 'G', 'I', 'L', 'O', 'Z', '*' )
%       7 blanks
%       8 width W
%       9 optional decimal point
%      10 optional mantissa M
%      11 blanks
%      12 optional ')'
%      13 blanks
%
%  Example:
%
%    S                 R   CODE   W    M
%
%    'I12              1   I      12   0
%    'E8.0'            1   E       8   0
%    'F10.5'           1   F      10   5
%    '2G14.6'          2   G      14   6
%    '*'               1   *      -1  -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string containing the
%    data to be read.  Reading will begin at position 1 and
%    terminate at the end of the string, or when no more
%    characters can be read.
%
%    Output, integer R, the repetition factor, which defaults to 1.
%
%    Output, character CODE, the format code.
%
%    Output, integer W, the field width.
%
%    Output, integer M, the mantissa width.
%
  LEFT = 1;
  RIGHT = -1;

  state = 0;
  paren_sum = 0;
  pos = 0;
  s_length = s_len_trim ( s );

  r = 0;
  w = 0;
  code = '?';
  m = 0;

  while ( pos < s_length )

    pos = pos + 1;
    c = s(pos);
%
%  BLANK character:
%
    if ( c == ' ' )

      if ( state == 4 )
        state = 5;
      elseif ( state == 6 )
        state = 7;
      elseif ( state == 10 )
        state = 11;
      elseif ( state == 12 )
        state = 13;
      end
%
%  LEFT PAREN
%
    elseif ( c == '(' )

      if ( state < 2 )
        paren_sum = paren_sum + LEFT;
      else
        state = -1;
        break;
      end
%
%  DIGIT (R, F, or W)
%
    elseif ( ch_is_digit ( c ) )

      if ( state <= 3 )
        state = 4;
        r = ch_to_digit ( c );
      elseif ( state == 4 )
        d = ch_to_digit ( c );
        r = 10 * r + d;
      elseif ( state == 6 || state == 7 )
        if ( code == '*' )
          state = -1;
          break;
        end
        state = 8;
        w = ch_to_digit ( c );
      elseif ( state == 8 )
        d = ch_to_digit ( c );
        w = 10 * w + d;
      elseif ( state == 9 )
        state = 10;
        m = ch_to_digit ( c );
      elseif ( state == 10 )
        d = ch_to_digit ( c );
        m = 10 * m + d;
      else
        state = -1;
        break;
      end
%
%  DECIMAL POINT
%
    elseif ( c == '.' )

      if ( state == 8 )
        state = 9;
      else
        state = -1;
        break;
      end
%
%  RIGHT PAREN
%
    elseif ( c == ')' )

      paren_sum = paren_sum + RIGHT;

      if ( paren_sum ~= 0 )
        state = -1;
        break;
      end

      if ( state == 6 && code == '*' )
        state = 12;
      elseif ( 6 <= state )
        state = 12;
      else
        state = -1;
        break;
      end
%
%  Code
%
    elseif ( ch_is_format_code ( c ) )

      if ( state < 6 )
        state = 6;
        code = c;
      else
        state = -1;
        break;
      end
%
%  Unexpected character
%
    else

      state = -1;
      break;

    end

  end

  if ( paren_sum ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'S_TO_FORMAT - Fatal error!\n' );
    fprintf ( 1, '  Parentheses mismatch.\n' );
    error ( 'S_TO_FORMAT - Fatal error!' );
  end

  if ( state < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'S_TO_FORMAT - Fatal error!\n' );
    fprintf ( 1, '  Parsing error.\n' );
    error ( 'S_TO_FORMAT - Fatal error!' );
  end

  if ( r == 0 )
    r = 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pdf = weibull_discrete_pdf ( x, a, b )

%*****************************************************************************80
%
%% WEIBULL_DISCRETE_PDF evaluates the discrete Weibull PDF.
%
%  Discussion:
%
%    PDF(X)(A,B) = ( 1 - A )**X**B - ( 1 - A )**(X+1)**B.
%
%    WEIBULL_DISCRETE_PDF(X)(A,1) = GEOMETRIC_PDF(X)(A)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X, the argument of the PDF.
%    0 <= X
%
%    Input, real A, B, the parameters that define the PDF.
%    0 <= A <= 1,
%    0 < B.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x < 0 )
    pdf = 0.0;
  else
    pdf = ( 1.0 - a )^(x^b) - ( 1.0 - a )^((x+1)^b);
  end

  return
end


text:
% lsinit;
% find first two points and establish correct scale
%

if size(alist,2)==0,
  % evaluate at absolutely smallest point
  alp=0;
  if amin>0, alp=amin; end;
  if amax<0, alp=amax; end;
  % new function value
  falp=feval(func,data,x+alp*p);
  alist=alp;flist=falp;
elseif size(alist,2)==1,
  % evaluate at absolutely smallest point
  alp=0;
  if amin>0, alp=amin; end;
  if amax<0, alp=amax; end;
  if alist~=alp,
    % new function value
    falp=feval(func,data,x+alp*p);
    alist=[alist,alp];flist=[flist,falp];
  end;
end;
aamin=min(alist);aamax=max(alist);
if amin>aamin | amax<aamax,
  alist,amin,amax
  error('non-admissible step in alist');
end;

% establish correct scale
if aamax-aamin<=scale,
  alp1=max(amin,min(-scale,amax));
  alp2=max(amin,min(+scale,amax));
  alp=inf;
  if aamin-alp1>=alp2-aamax, alp=alp1; end;
  if alp2-aamax>=aamin-alp1, alp=alp2; end;
  if alp<aamin | alp>aamax,
    % new function value
    falp=feval(func,data,x+alp*p);
    alist=[alist,alp];flist=[flist,falp];
  end;
end;

if size(alist,2)==1,
  scale,aamin,aamax,alp1,alp2
  error('lsinit bug: no second point found');
end;

lssort;

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function tnet = vl_simplenn_tidy(net)
%VL_SIMPLENN_TIDY  Fix an incomplete or outdated SimpleNN network.
%   NET = VL_SIMPLENN_TIDY(NET) takes the NET object and upgrades
%   it to the current version of MatConvNet. This is necessary in
%   order to allow MatConvNet to evolve, while maintaining the NET
%   objects clean. This function ignores custom layers.
%
%   The function is also generally useful to fill in missing default
%   values in NET.
%
%   See also: VL_SIMPLENN().

% Copyright (C) 2014-15 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

tnet = struct('layers', {{}}, 'meta', struct()) ;

% copy meta information in net.meta subfield
if isfield(net, 'meta')
  tnet.meta = net.meta ;
end

if isfield(net, 'classes')
  tnet.meta.classes = net.classes ;
end

if isfield(net, 'normalization')
  tnet.meta.normalization = net.normalization ;
end

% copy layers
for l = 1:numel(net.layers)
  defaults = {'name', sprintf('layer%d', l), 'precious', false};
  layer = net.layers{l} ;

  % check weights format
  switch layer.type
    case {'conv', 'convt', 'bnorm'}
      if ~isfield(layer, 'weights')
        layer.weights = {...
          layer.filters, ...
          layer.biases} ;
        layer = rmfield(layer, 'filters') ;
        layer = rmfield(layer, 'biases') ;
      end
  end
  if ~isfield(layer, 'weights')
    layer.weights = {} ;
  end

  % check that weights inlcude moments in batch normalization
  if strcmp(layer.type, 'bnorm')
    if numel(layer.weights) < 3
      layer.weights{3} = ....
        zeros(numel(layer.weights{1}),2,'single') ;
    end
  end

  % fill in missing values
  switch layer.type
    case {'conv', 'pool'}
      defaults = [ defaults {...
        'pad', 0, ...
        'stride', 1, ...
        'opts', {}}] ;

    case 'convt'
      defaults = [ defaults {...
        'crop', 0, ...
        'upsample', 1, ...
        'numGroups', 1, ...
        'opts', {}}] ;

    case {'pool'}
      defaults = [ defaults {...
        'method', 'max', ...
        'pad', 0, ...
        'stride', 1, ...
        'opts', {}}] ;

    case 'relu'
      defaults = [ defaults {...
        'leak', 0}] ;

    case 'dropout'
      defaults = [ defaults {...
        'rate', 0.5}] ;

    case {'normalize', 'lrn'}
      defaults = [ defaults {...
        'param', [5 1 0.0001/5 0.75]}] ;

    case {'pdist'}
      defaults = [ defaults {...
        'noRoot', false, ...
        'aggregate', false, ...
        'p', 2, ...
        'epsilon', 1e-3, ...
        'instanceWeights', []} ];
  end

  for i = 1:2:numel(defaults)
    if ~isfield(layer, defaults{i})
      layer.(defaults{i}) = defaults{i+1} ;
    end
  end

  % save back
  tnet.layers{l} = layer ;
end

repo: CREST-Release


text:
%CONVEXHULL  Finds the convex hull of a point set
%
%     hull = cv.convexHull(points)
%     hull = cv.convexHull(points, 'OptionName', optionValue, ...)
%
% ## Input
% * __points__ Input 2D point set, stored in numeric array
%   (Nx2/Nx1x2/1xNx2) or cell array of 2-element vectors (`{[x,y], ...}`).
%
% ## Output
% * __hull__ Output convex hull. It is either an integer vector of indices or
%   vector of points. In the first case, the hull elements are 0-based indices
%   of the convex hull points in the original array (since the set of convex
%   hull points is a subset of the original point set). In the second case,
%   hull elements are the convex hull points themselves. In case output is the
%   hull points, it has the same type as the input.
%
% ## Options
% * __ReturnPoints__ Operation flag. In case of a matrix, when the flag is
%   true, the function returns convex hull points (Mx2 matrix). Otherwise, it
%   returns indices of the convex hull points (vector of length M). In case
%   the input is a cell-array, when the flag is true, the function return
%   convex hull points (as a cell-array of points). Otherwise it returns
%   indices of points (vector of length M). default true
% * __Clockwise__ Orientation flag. If it is true, the output convex hull is
%   oriented clockwise. Otherwise, it is oriented counter-clockwise. The usual
%   screen coordinate system is assumed so that the origin is at the top-left
%   corner, x axis is oriented to the right, and y axis is oriented downwards.
%   default false
%
% The function cv.convexHull finds the convex hull of a 2D point set using the
% Sklansky's algorithm [Sklansky82] that has `O(N logN)` complexity in the
% current implementation. See the sample `convexhull_demo.m` that demonstrates
% the usage of the function.
%
% ## References
% [Sklansky82]:
% > Jack Sklansky. "Finding the convex hull of a simple polygon".
% > Pattern Recognition Letters, 1(2):79-83, 1982.
%
% See also: convhull, convhulln, bwconvhull, boundary
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function x = sample_vector(p,col)
%SAMPLE_VECTOR    Sample from multiple categorical distributions.
% X = SAMPLE_VECTOR(P) returns a row vector of cols(P) integers, where
% X(I) = SAMPLE(P(:,I)).
%
% X = SAMPLE_VECTOR(P,COL) returns a row vector of length(COL) integers, where
% X(I) = SAMPLE(P(:,COL(I))).
% This is equivalent to SAMPLE_VECTOR(P(:,COL)), but faster.

% Written by Tom Minka
% (c) Microsoft Corporation. All rights reserved.

cdf = cumsum(p);
if nargin >= 2
  cdf = cdf(:,col);
end
if any(cdf(end,:) <= 0)
  error('distribution is all zeros');
end
u = rand(1,cols(cdf)).*cdf(end,:);
x = col_sum(cdf < repmat(u,rows(cdf),1)) + 1;

author: andrewssobral

repo: lrslibrary


text:
function result = arg_report(type,func,args)
% Report information of a certain Type about the given Function.
% Result = arg_report(Type,Function,Arguments)
%
% Functions that declare their arguments via arg_define() make their parameter specification
% accessible to outside functions. This can be used to display auto-generated settings dialogs, to
% record function calls, and so on. Some other functions can enable reports, as well (like
% declare_properties() and expose_handles()).
%
% In:
%   Type : Type of information to report. The currently defined report types are as follows:
%
%          'vals' : Report the values of the function's arguments as a struct, possibly with
%                   sub-structs, after assignment of Arguments. (provided by arg_define)
%          'nvps' : Report the values of the function's arguments as cell array of name-value pairs.
%                   (provided by arg_define)
%          'lean' : Report a lean specification of the function's arguments as a struct array, with
%                   fields as in arg_specifier, including help text, etc. This is after assignment
%                   of Arguments. (provided by arg_define)
%          'rich' : Report a rich specification of the function's arguments as a struct array, with
%                   fields as in arg_specifier. In addition to lean this includes information about
%                   alternative (non-default) settings, for use in GUI generation. (provided by
%                   arg_define)
%
%          'properties' : Report declared properties of the function, if any. Arguments must be empty.
%                         (provided by declare_properties)
%
%          'handle': Report function handles to scoped and nested functions within the Function. The
%                    names of those functions are given as a cell array in Arguments, and the result
%                    is a cell array of associated function handles (or one handle if only one
%                    function was requested).
%
%   Function : Handle to a function that supports the given report type.
%
%   Arguments : Cell array of inputs to be passed to the function.
%
% Out:
%   Result : the reported data.
%
% Examples:
%   % for a function call with some arguments assigned, obtain a struct with all parameter
%   % names and values, including defaults
%   params = arg_report('vals',@myfunction,{4,10,true,'option1','xxx','option5',10})
%
%   % obtain a specification of all function arguments, with defaults, help text, type, shape, and other
%   % meta-data (with a subset of settings customized according to arguments)
%   spec = arg_report('rich',@myfunction,myarguments)
%
%   % obtain a report of properties of the function (declared via declared_properties() within the
%   % function)
%   props = arg_report('properties',@myfunction)
%
% See also:
%   arg_define, declare_properties, expose_handles
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-09-24

% Copyright (C) Christian Kothe, SCCN, 2010, christian@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify it under the terms of the GNU
% General Public License as published by the Free Software Foundation; either version 2 of the
% License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
% even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with this program; if not,
% write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA

% check inputs
if ischar(func)
    func = str2func(func); end
if ~isa(func,'function_handle')
    error('The given Function argument must be a function handle.'); end
if ~ischar(type)
    error('The given Type argument must be a string.'); end
if nargin < 3
    if strcmp(type,'properties')
        % for the properties report we implicitly pad the arguments with blanks to allow functions
        % to define up to this number of traditional non-varargin arguments prior to the varargin
        args(1:15) = {[]};
    else
        args = {};
    end
elseif isequal(args,[])
    args = {};
elseif ~iscell(args)
    error('Arguments must be a cell array, if given.');
end

% disable the expression features of functions called by arg_report
persistent handle_expressions;
if isempty(handle_expressions)
    handle_expressions = exist('exp_eval','file'); end

result = {};
try
    % evaluate the funtion with two special arguments appended
    if handle_expressions
        hlp_scope({'disable_expressions',1},func,args{:},'__arg_report__',type);
    else
        func(args{:},'__arg_report__',type);
    end
catch report
    if strcmp(report.identifier,'BCILAB:arg:report_args')
        global tracking; %#ok<TLEV>
        % get the ticket of the report
        ticket = sscanf(report.message(end-4:end),'%u');
        % read out the payload and return the ticket
        result = tracking.arg_sys.reports{ticket};
        tracking.arg_sys.tickets.addLast(ticket);
    elseif ~strcmp(type,'properties')
        % got a genuine error
        settings = dbstatus;
        if any(strcmp({settings.cond},'error')) && ~hlp_resolve('disable_dbstop_if_error_msg',false)
            % rerun the code to propagate it properly (we're not using rethrow here to not confuse dbstop if error)
            disp_once('arg_report: caught error while in "dbstop if error" mode; re-running offending code to get break point...\n  error traceback was: %s\n',hlp_handleerror(report));
            if handle_expressions
                hlp_scope({'disable_expressions',1},func,args{:},'__arg_report__',type);
            else
                func(args{:},'__arg_report__',type);
            end
        else
            rethrow(report);
        end
    end
end

repo: REST


text:
function tmp_str=cat_str(cellstrs)

tmp_str=[];
for c_idx=1:length(cellstrs)
    if ~isempty(tmp_str)
        tmp_str=[tmp_str '_'];
    end
    tmp_str=[tmp_str cellstrs{c_idx}];
end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
%% initializeWecSim
% Initialize WEC-Sim run, called from wecSim.m
%
%% Following Classes are required to be defined in the WEC-Sim input file:
%%
%% simu = simulationClass();                                               - To Create the Simulation Variable
%%
%% waves = waveClass('<wave type');                                       - To create the Wave Variable and Specify Type
%%
%% body(<body number>) = bodyClass('<hydrodynamics data file name>.h5');   - To initialize bodyClass:
%%
%% constraint(<constraint number>) = constraintClass('<Constraint name>'); - To initialize constraintClass
%%
%% pto(<PTO number>) = ptoClass('<PTO name>');                             - To initialize ptoClass
%%
%% mooring(<mooring number>) = mooringClass('<Mooring name>');             - To initialize mooringClass (only needed when mooring blocks are used)
%%
%%

%% Start WEC-Sim log
% Clear old input, plots, log file and start new log file.
diary off
clear body waves simu output pto constraint ptoSim mooring values names InParam
try delete('*.log'); end
diary('simulation.log')

%% Add 'temp' directory
% Store root directory of this *.m file
projectRootDir = pwd;
if exist('pctDir')
    projectRootDir = [projectRootDir filesep pctDir];
end

% Create 'temp' directory if it doesn't exist and add to 'temp' path
warning('off','MATLAB:MKDIR:DirectoryExists')
if exist('pctDir')
    mkdir ([pctDir filesep 'temp'])
else
    if mkdir('temp') == 0
        mkdir 'temp'
    end
end
addpath(fullfile(projectRootDir,'temp'),'-end');

% Save Simulink-generated helper files to 'temp' directory
Simulink.fileGenControl('set',...
    'CacheFolder',fullfile(projectRootDir,'temp'))

%% Read input file
tic

% Set input parameters based on how the simulation is called
if exist('runWecSimCML','var') && runWecSimCML==1
    % wecSim input from wecSimInputFile.m of case directory in the standard manner
    fprintf('\nWEC-Sim Input From Standard wecSimInputFile.m Of Case Directory... \n');
    bdclose('all');
    run('wecSimInputFile');
else
    % Get global reference frame parameters
    values = get_param([bdroot,'/Global Reference Frame'],'MaskValues');    % Cell array containing all Masked Parameter values
    names = get_param([bdroot,'/Global Reference Frame'],'MaskNames');      % Cell array containing all Masked Parameter names
    j = find(strcmp(names,'InputMethod'));
    if strcmp(values{j},'Input File')
        % wecSim input from input file selected in Simulink block
        fprintf('\nWEC-Sim Input From File Selected In Simulink... \n');
        i = find(strcmp(names,'InputFile'));
        run(values{i});
    elseif strcmp(values{j},'Custom Parameters')
        % wecSim input from custom parameters in Simulink block
        fprintf('\nWEC-Sim Input From Custom Parameters In Simulink... \n');
        inputFile = 'wecSimInputFile_simulinkCustomParameters';
        writeInputFromBlocks(inputFile);
        run(inputFile);
    end
end
clear values names i j;

% Read Inputs for Multiple Conditions Run
try fprintf('wecSimMCR Case %g\n',imcr); end

if exist('mcr','var') == 1
    for n=1:length(mcr.cases(1,:))
        if strcmp(mcr.header{n},'LoadFile')                                 % added this flag for system ID efforts!
            load(mcr.cases{imcr,n});
        elseif iscell(mcr.cases)
            eval([mcr.header{n} '= mcr.cases{imcr,n};']);
        else
            eval([mcr.header{n} '= mcr.cases(imcr,n);']);
        end
    end; clear n combine;
    try
        waves.spectrumFile = ['..' filesep pctDir filesep '..' filesep waves.spectrumFile];
        waves.elevationFile = ['..' filesep pctDir filesep '..' filesep waves.elevationFile];
    end
end

% Waves and Simu: check inputs
waves.checkInputs();
simu.checkInputs();

% Constraints: count & set orientation
if exist('constraint','var') == 1
    simu.numConstraints = length(constraint(1,:));
    for ii = 1:simu.numConstraints
        constraint(ii).setNumber(ii);
        constraint(ii).setOrientation();
    end; clear ii
end

% PTOs: count & set orientation & set pretension
if exist('pto','var') == 1
    simu.numPtos = length(pto(1,:));
    for ii = 1:simu.numPtos
        pto(ii).setNumber(ii);
        pto(ii).setOrientation();
        pto(ii).setPretension();
    end; clear ii
end

% Mooring Configuration: count
if exist('mooring','var') == 1
    simu.numMoorings = length(mooring(1,:));
    for ii = 1:simu.numMoorings
        mooring(ii).setLoc();
        mooring(ii).setNumber(ii);
    end; clear ii
end

% Bodies: count, check inputs, read hdf5 file, and check inputs
numHydroBodies = 0;
numNonHydroBodies = 0;
numDragBodies = 0;
hydroBodLogic = zeros(length(body(1,:)),1);
nonHydroBodLogic = zeros(length(body(1,:)),1);
dragBodLogic = zeros(length(body(1,:)),1);
for ii = 1:length(body(1,:))
    body(ii).setNumber(ii);
    body(ii).checkInputs(simu.explorer);
    if body(ii).nonHydro==0
        if numNonHydroBodies > 0 || numDragBodies > 0
            error('All hydro bodies must be specified before any drag or non-hydro bodies.')
        end
        numHydroBodies = numHydroBodies + 1;
        hydroBodLogic(ii) = 1;
    elseif body(ii).nonHydro==1
        numNonHydroBodies = numNonHydroBodies + 1;
        nonHydroBodLogic(ii) = 1;
    elseif body(ii).nonHydro==2
        numDragBodies = numDragBodies + 1;
        dragBodLogic(ii) = 1;
    end
end
simu.numHydroBodies = numHydroBodies; clear numHydroBodies
simu.numDragBodies = numDragBodies; clear numDragBodies
for ii = 1:simu.numHydroBodies
    body(ii).setDOF(simu.numHydroBodies,simu.b2b);

    % Determine if hydro data needs to be reloaded from h5 file, or if hydroData
    % was stored in memory from a previous run.
    if exist('totalNumOfWorkers','var') ==0 && exist('mcr','var') == 1 && simu.reloadH5Data == 0 && imcr > 1
        body(ii).loadHydroData(hydroData(ii));
    else
        % check for correct h5 file
        h5Info = dir(body(ii).h5File);
        h5Info.bytes;
        if h5Info.bytes < 1000
            error(['This is not the correct *.h5 file. Please install git-lfs to access the correct *.h5 file, or run \hydroData\bemio.m to generate a new *.h5 file'])
        end
        clearvars h5Info
        % Read hydro data from BEMIO and load into the bodyClass
        tmp_hydroData = readBEMIOH5(body(ii).h5File, body(ii).number, body(ii).meanDrift);
        body(ii).loadHydroData(tmp_hydroData);
        clear tmp_hydroData
    end
end; clear ii

% Cable Configuration: count, set Cg/Cb, PTO loc, L0 and initialize bodies
if exist('cable','var')==1
    simu.numCables = length(cable(1,:));
    for ii = 1:simu.numCables
        cable(ii).setNumber(ii);
        cable(ii).setCg();
        cable(ii).setCb();
        cable(ii).setTransPTOLoc();
        cable(ii).setLength();
        cable(ii).dragForcePre(simu.rho);
        cable(ii).setVolume(simu.rho);
        cable(ii).setOrientation();
        cable(ii).linearDampingMatrix();
    end
end

% PTO-Sim: read input, count
if exist('ptoSim','var') == 1
    simu.numPtoSim = length(ptoSim(1,:));
    for ii = 1:simu.numPtoSim
        ptoSim(ii).number = ii;
    end; clear ii
end

toc

%% Pre-processing start
tic
fprintf('\nWEC-Sim Pre-processing ...   \n');
if exist('pctDir')
    cd(pctDir);
end

%% HydroForce Pre-Processing: Wave Setup & HydroForcePre.
% simulation setup
simu.setup();

% wave setup
if any(hydroBodLogic == 1)
    % When hydro bodies (and an .h5) are present, define the wave using those
    % parameters.
    waves.setup(body(1).hydroData.simulation_parameters.w, body(1).hydroData.simulation_parameters.waterDepth, simu.rampTime, simu.dt, simu.maxIt, simu.time, simu.gravity, simu.rho);
    % Check that direction and freq are within range of hydro data
    if  min(waves.direction) <  min(body(1).hydroData.simulation_parameters.direction) || max(waves.direction) >  max(body(1).hydroData.simulation_parameters.direction)
        error('waves.direction outside of range of available hydro data')
    end
    if strcmp(waves.type,'elevationImport')~=1 && strcmp(waves.type,'noWave')~=1 && strcmp(waves.type,'noWaveCIC')~=1
        if  min(waves.omega) <  min(body(1).hydroData.simulation_parameters.w) || max(waves.omega) >  max(body(1).hydroData.simulation_parameters.w)
            error('waves.omega outside of range of available hydro data')
        end
    end
else
    % When no hydro bodies (and no .h5) are present, define the wave using
    % input file parameters
    waves.setup([], [], simu.rampTime, simu.dt, simu.maxIt, simu.time, simu.gravity, simu.rho);
end

% Nonlinear hydro
for kk = 1:length(body(1,:))
    if (body(kk).nonlinearHydro > 0) || (simu.paraview.option == 1)
        body(kk).importBodyGeometry(simu.domainSize)
    end
end; clear kk

% hydroForcePre
idx = find(hydroBodLogic==1);
if ~isempty(idx)
    for kk = 1:length(idx)
        it = idx(kk);
        body(it).hydroForcePre(waves.omega,waves.direction,simu.cicTime,waves.bem.count,simu.dt,...
            simu.rho,simu.gravity,waves.type,waves.waveAmpTime,simu.stateSpace,simu.b2b);
    end; clear kk idx
end

% nonHydroPre
idx = find(nonHydroBodLogic==1);
if ~isempty(idx)
    for kk = 1:length(idx)
        it = idx(kk);
        body(it).nonHydroForcePre(simu.rho);
    end; clear kk idx
end

% dragBodyPre
idx = find(dragBodLogic == 1);
if ~isempty(idx)
    for kk = 1:length(idx)
        it = idx(kk);
        body(it).dragForcePre(simu.rho);
    end; clear kk idx
end

% Check cicEndTime
if waves.typeNum~=0 && waves.typeNum~=10
    for iBod = 1:simu.numHydroBodies
        if simu.cicEndTime > max(body(iBod).hydroData.hydro_coeffs.radiation_damping.impulse_response_fun.t)
            error('simu.cicEndTime is larger than the length of the IRF');
        end
    end
end

% Check that the hydro data for each body is given for the same frequencies
for ii = 1:simu.numHydroBodies
    if length(body(1).hydroData.simulation_parameters.w) ~= length(body(ii).hydroData.simulation_parameters.w)
        error('BEM simulations for each body must have the same number of frequencies')
    else
        for jj = 1:length(body(1).hydroData.simulation_parameters.w)
            if body(1).hydroData.simulation_parameters.w(jj) ~= body(ii).hydroData.simulation_parameters.w(jj)
                error('BEM simulations must be run with the same frequencies.')
            end; clear jj;
        end
    end
end; clear ii;

% Check for elevationImport with nonlinearHydro
for ii = 1:simu.numHydroBodies
    if strcmp(waves.type,'elevationImport') && body(ii).nonlinearHydro == 1
        error(['Cannot run WEC-Sim with nonlinear hydro (body(ii).nonlinearHydro) and "elevationImport" wave type'])
    end
end

% Check for elevationImport with morisonElement
for ii = 1:simu.numHydroBodies
    if strcmp(waves.type,'elevationImport') && body(ii).morisonElement.option ~= 0
        error(['Cannot run WEC-Sim with Morison Element (body(ii).morisonElement.option>0) and "elevationImport" wave type'])
    end
end

% Check for morisonElement inputs for body(ii).morisonElement.option == 1 || body(ii).morisonElement.option == 2
for ii = 1:length(body(1,:))
    if body(ii).morisonElement.option == 1
        if body(ii).nonHydro ~=1
            [rgME,~] = size(body(ii).morisonElement.rgME);
            for jj = 1:rgME
                if true(isfinite(body(ii).morisonElement.z(jj,:))) == true
                    warning(['"body.morisonElement.z" is not used for "body.morisonElement.option = 1". Check body ',num2str(ii),' element ',num2str(jj)])
                end
                if length(body(ii).morisonElement.cd(jj,:)) ~= 3 || length(body(ii).morisonElement.ca(jj,:)) ~= 3 || length(body(ii).morisonElement.area(jj,:)) ~= 3
                    error(['cd, ca, and area coefficients for each elelement for "body.morisonElement.option = 1" must be of size [1x3] and all columns of data must be real and finite. Check body ',num2str(ii),' element ',num2str(jj),' coefficients'])
                end
            end; clear jj
        else
            if body(ii).morisonElement.option == 1 || body(ii).morisonElement.option == 2
                warning(['Morison elements are not available for non-hydro bodies. Please check body ',num2str(ii),' inputs.'])
            end
        end
    elseif body(ii).morisonElement.option == 2
        if body(ii).nonHydro ~=1
            [rgME,~] = size(body(ii).morisonElement.rgME);
            for jj = 1:rgME
                if body(ii).morisonElement.cd(jj,3) ~= 0 || body(ii).morisonElement.ca(jj,3) ~= 0 || body(ii).morisonElement.area(jj,3) ~= 0
                    warning(['cd, ca, and area coefficients for "body.morisonElement.option == 2" must be of size [1x2], third column of data is not used. Check body ',num2str(ii),' element ',num2str(jj),' coefficients'])
                end
            end; clear jj
        else
            if body(ii).morisonElement.option ==1 || body(ii).morisonElement.option ==2
                warning(['Morison elements are not available for non-hydro bodies. Please check body ',num2str(ii),' inputs.'])
            end
        end
    end; clear ii
end

%% Set variant subsystems options
for ii=1:length(body(1,:))
    if body(ii).nonHydro==0
        % Nonlinear FK Force Variant Subsystem
        eval(['nonLinearHydro_' num2str(ii) ' = body(ii).nonlinearHydro;']);
        eval(['sv_b' num2str(ii) '_linearHydro = Simulink.Variant(''nonLinearHydro_', num2str(ii), '==0'');']);
        eval(['sv_b' num2str(ii) '_nonlinearHydro=Simulink.Variant(''nonLinearHydro_', num2str(ii), '>0'');']);
        % Nonlinear Wave Elevation Variant Subsystem
        eval(['sv_b' num2str(ii) '_meanFS=Simulink.Variant(''nonLinearHydro_', num2str(ii), '<2'');']);
        eval(['sv_b' num2str(ii) '_instFS=Simulink.Variant(''nonLinearHydro_', num2str(ii), '==2'');']);
    end
end; clear ii;
% Morison Element
for ii=1:length(body(1,:))
    if body(ii).nonHydro ~=1
        eval(['morisonElement_' num2str(ii) ' = body(ii).morisonElement.option;'])
        eval(['sv_b' num2str(ii) '_MEOff = Simulink.Variant(''morisonElement_' num2str(ii) '==0'');'])
        eval(['sv_b' num2str(ii) '_MEOn = Simulink.Variant(''morisonElement_' num2str(ii) '==1 || morisonElement_' num2str(ii) '==2'');'])
    end
end; clear ii;
% Radiation Damping
if waves.typeNum==0 || waves.typeNum==10 %'noWave' & 'regular'
    radiation_option = 1;
elseif simu.stateSpace == 1
    radiation_option = 3;
else
    radiation_option = 2;
end
sv_constantCoeff=Simulink.Variant('radiation_option==1');
sv_convolution=Simulink.Variant('radiation_option==2');
sv_stateSpace=Simulink.Variant('radiation_option==3');
% Wave type
typeNum = waves.typeNum;
sv_noWave=Simulink.Variant('typeNum<10');

% Passive Yaw
for ii=1:length(body(1,:))
    eval(['yaw_' num2str(ii) ' = body(ii).yaw.option;'])
    eval(['sv_regularWaves_b' num2str(ii) '= Simulink.Variant(''typeNum>=10 && typeNum<20 && yaw_', num2str(ii), '==0'');'])
    eval(['sv_regularWavesYaw_b' num2str(ii) '= Simulink.Variant(''typeNum>=10 && typeNum<20 && yaw_' num2str(ii) '==1'');'])
    eval(['sv_irregularWaves_b' num2str(ii) '= Simulink.Variant(''typeNum>=20 && typeNum<30 && yaw_' num2str(ii) '==0'');'])
    eval(['sv_irregularWavesYaw_b' num2str(ii) '= Simulink.Variant(''typeNum>=20 && typeNum<30 && yaw_' num2str(ii) '==1'');'])
end; clear ii

sv_udfWaves=Simulink.Variant('typeNum>=30');
% Body2Body
B2B = simu.b2b;
sv_noB2B=Simulink.Variant('B2B==0');
sv_B2B=Simulink.Variant('B2B==1');
numBody=simu.numHydroBodies;
% nonHydro
for ii=1:length(body(1,:))
    eval(['nhbody_' num2str(ii) ' = body(ii).nonHydro;'])
    eval(['sv_b' num2str(ii) '_hydroBody = Simulink.Variant(''nhbody_' num2str(ii) '==0'');'])
    eval(['sv_b' num2str(ii) '_nonHydroBody = Simulink.Variant(''nhbody_' num2str(ii) '==1'');'])
    eval(['sv_b' num2str(ii) '_dragBody = Simulink.Variant(''nhbody_' num2str(ii) '==2'');'])
end; clear ii

%Efficiency model for hydraulic motor PTO-Sim block
try
    for ii=1:length(ptoSim(1,:))
        eval(['EffModel_' num2str(ii) ' = ptoSim(ii).hydraulicMotor.effModel;']);
        eval(['sv_b' num2str(ii) '_AnalyticalEfficiency = Simulink.Variant(''EffModel_', num2str(ii), '==1'');']);
        eval(['sv_b' num2str(ii) '_TabulatedEfficiency = Simulink.Variant(''EffModel_', num2str(ii), '==2'');']);
    end; clear ii;
end

% Visualization Blocks
if ~isempty(waves.marker.location) && typeNum < 30
    visON = 1;
else
    visON = 0;
end
sv_visualizationON  = Simulink.Variant('visON==1');
sv_visualizationOFF = Simulink.Variant('visON==0');

%% End Pre-Processing and Output All the Simulation and Model Setting
toc
simu.listInfo(waves.typeNum);
waves.listInfo
fprintf('\nList of Body: ');
fprintf('Number of Bodies = %u \n',simu.numHydroBodies)
for i = 1:simu.numHydroBodies
    body(i).listInfo
end;  clear i
fprintf('\nList of PTO(s): ');
if (exist('pto','var') == 0)
    fprintf('No PTO in the system\n')
else
    fprintf('Number of PTOs = %G \n',length(pto(1,:)))
    for i=1:simu.numPtos
        pto(i).listInfo
    end; clear i
end
fprintf('\nList of Constraint(s): ');
if (exist('constraint','var') == 0)
    fprintf('No Constraint in the system\n')
else
    fprintf('Number of Constraints = %G \n',length(constraint(1,:)))
    for i=1:simu.numConstraints
        constraint(i).listInfo
    end; clear i
end
fprintf('\n')

%% Load simMechanics file & Run Simulation
tic
fprintf('\nSimulating the WEC device defined in the SimMechanics model %s...   \n',simu.simMechanicsFile)
% Modify some stuff for simulation
for iBod = 1:simu.numHydroBodies
    body(iBod).adjustMassMatrix(simu.adjMassFactor,simu.b2b);
end; clear iBod
warning('off','Simulink:blocks:TDelayTimeTooSmall');
warning('off','Simulink:blocks:BusSelDupBusCreatorSigNames');
warning('off','MATLAB:loadlibrary:FunctionNotFound');
warning('off','MATLAB:loadlibrary:parsewarnings');
warning('off','MATLAB:printf:BadEscapeSequenceInFormat');
warning('off','Simulink:blocks:DivideByZero');
warning('off','sm:sli:setup:compile:SteadyStateStartNotSupported')
set_param(0, 'ErrorIfLoadNewModel', 'off')
% Load parameters to Simulink model
simu.loadSimMechModel(simu.simMechanicsFile);

toc

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


function [y,varargout] = bootstrap_(x,Nboot,replacement)

if nargin<3
    replacement = 0;
else
    if replacement>1
        error('replacement is logical, 0 or 1.')
    end
end

[T,N] = size(x);

y = nan(T,N,Nboot);

for nboot =1 :Nboot
    if replacement ==1
        % replacement
        order_(:,nboot)  = randi(T,[T,1]);
    else
        % no replacement
        order_(:,nboot)  = randperm(T);
    end
    y(:,:,nboot) = x(order_(:,nboot),:);
end

varargout = {order_};


text:
function y = vec(x)
    y = x(:);
end

author: stephenbeckr

repo: fastRPCA


text:
% -------------------------------------------------------------------------
%
%     Demo software for V-BM4D volumetric data denoising
%               Release ver. 1.0  (11 December 2014)
%
% -------------------------------------------------------------------------
%
% The software implements the V-BM4D denoising algorithm described in:
%
%    M. Maggioni, G. Boracchi, A. Foi, K. Egiazarian, "Video Denoising
%      Using Separable 4D Nonlocal Spatiotemporal Transforms",
%      Proc. SPIE Electronic Imaging 2011, San Francisco, CA, USA.
%
%    M. Maggioni, G. Boracchi, A. Foi, K. Egiazarian, "Video Denoising,
%      Deblocking and Enhancement Through Separable 4-D Nonlocal
%      Spatiotemporal Transforms", IEEE Trans. on Image Proc.,
%      Vol. 21, No. 9, Sep. 2012. doi:10.1109/TIP.2012.2199324
%
% -------------------------------------------------------------------------
%
% authors:               Matteo Maggioni
%                        Alessandro Foi
%
% web page:              http://www.cs.tut.fi/~foi/GCF-BM3D
%
% contact:               firstname.lastname@tut.fi
%
% -------------------------------------------------------------------------
% Copyright (c) 2010-2014 Tampere University of Technology.
% All rights reserved.
% This work should be used for nonprofit purposes only.
% -------------------------------------------------------------------------
%
% Disclaimer
% ----------
%
% Any unauthorized use of these routines for industrial or profit-oriented
% activities is expressively prohibited. By downloading and/or using any of
% these files, you implicitly agree to all the terms of the TUT limited
% license (included in the file Legal_Notice.txt).
% -------------------------------------------------------------------------

clear all;close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Modifiable parameters
file_name = '';      % Video file name. If empty, a prompt dialog will appear
sigma     = 25;      % Noise standard deviation. it should be in the same
                     % intensity range of the video
profile = 'lc';      % V-BM4D parameter profile
                     %  'lc' --> low complexity
                     %  'np' --> normal profile
do_wiener = 1;       % Wiener filtering
                     %   1 --> enable Wiener filtering
                     %   0 --> disable Wiener filtering
sharpen = 1;         % Sharpening
                     %   1 --> disable sharpening
                     %  >1 --> enable sharpening
deflicker = 1;       % Deflickering
                     %   1 --> disable deflickering
                     %  <1 --> enable deflickering
verbose = 1;         % Verbose mode

est_noise = 0;       % Enable noise estimation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       MODIFY BELOW THIS POINT ONLY IF YOU KNOW WHAT YOU ARE DOING       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% If no file_name is given, let the user select a file from local disk
if isempty(file_name)
    [ file_name, folder ] = uigetfile({'*.avi'}, 'Select Video');
    if isequal(file_name,0)
        error('No file selected.')
    end
    file_name = [folder,file_name];
end
if strcmpi(file_name(end-2:end),'mat')
    % if a matlab file is given, the noise-free video is assumed to be
    % saved as a 3-D matrix called "y"
    load(file_name)
else
    % video formats mast be readable by VideoReader ('help VideoReader')
    y = read_video(file_name);
end

% Scaling data
S = 255;
I_MAX = 255/S;
% Create synthetic Gaussian noise
randn('seed',0);
y = cast(y, 'single')/S;
sigma = sigma/S;
z = y + sigma*randn(size(y));
% Set sigma to -1 to enable noise estimation
if est_noise
    sigma = -1;
end

% V-BM4D filtering
disp('Denoising started')
tic;
y_est = vbm4d( z, sigma, profile, do_wiener, sharpen, deflicker, verbose );
time = toc;
PSNR = 10*log10(I_MAX^2/mean((y(:)-y_est(:)).^2));
fprintf('Denoising completed (%.1fs / %.1ffps): PSNR %.2fdB \n', time, size(z,ndims(z))/time, PSNR);

% Show results
f = round(size(y,3)/2);
figure,imshow([y(:,:,f) z(:,:,f) y_est(:,:,f)],[])

repo: DeSCI


text:
function anim = computeSMFromW( isProj, W, varargin )
% Recover rigid structure and motion from 2D calibrated measurements W
%
% Essential matrix decomposition
% Reference: HZ2, p259, Result 9.19, and p. 294
% isProj && nFrame==2 && isCalibrated
% fast
%
% Normalized 8-point algorithm
% Reference: HZ2, p279 and alg 11.1 p. 282
% isProj && nFrame==2 && ~isCalibrated
% fast
%
% Projective Factorization
% Reference: HZ2, p445, Algorithm 18.2
% Sturm and Triggs ECCV 96
% A Factorization Based Algorithm for Multi-Image Projective Structure and
% Motion
% isProj && nFrame>2 && method==0
%
% Projective Factorization
% Reference: Iterative Extensions of the Sturm/Triggs Algorithm:
% Convergence and Nonconvergence, from Oliensis, Hartley, PAMI 07
% isProj && nFrame>2 && method==Inf
% slower
%
% Gold Standard for Affine camera matrix
% Reference: HZ2, p351, Algorithm 14.1
% ~isProj && nFrame==2 && onlyErrorFlag
% fast
%
% Tomasi Kanade without the metric constraint
% Affine camera matrix, MLE estimation (Tomasi Kanade)
% Reference: HZ2, p437, Algorithm 18.1
% ~isProj && nFrame>2
% fast
%
% Tomasi Kanade with the metric constraint
% ~isProj && nFrame>2 && isCalibrated && method=0
% fast
%
% Tomasi Kanade with the metric constraint (more accurate computation)
% ~isProj && nFrame>2 && isCalibrated && method=inf
% slow
%
% Tomasi Kanade and autocalibration
% ~isProj && nFrame>2 && ~isCalibrated
% slow
%
% if doAffineUpgrade and/or doMetricUpgrade are true, affine and metric
% upgrades from Chandraker IJCV 2009 is applied
% slow
%
% If there are any missing entries:
% Low-Rank Matrix Fitting Based on Subspace Perturbation Analysis
% with Applications to Structure from Motion
% Hongjun Jia, Aleix M. Martinez, PAMI 08
%
% Returns the structure and camera parameters in an Animation object
%
% USAGE
%   anim = computeSMFromW( isProj, W, 'method',method )
%   anim = computeSMFromW( isProj, W,'method',Inf,'onlyErrorFlag',true )
%
% INPUTS
%  isProj     - flag indicating if the camera is projective
%  W          - [] [ 2 x nPoint x nFrame ] 2D projected features (NaN if
%               missing entry
%  varargin   - list of paramaters in quotes alternating with their values
%       - 'isCalibrated' [false] flag indicating if the cameras are
%                        calibrated (intrinsic parameters are identity)
%       - 'doAffineUpgrade' [false] flag indicating if the affine upgrade
%                        is computed
%       - 'doMetricUpgrade' [false] flag indicating if the metric upgrade
%                        is computed
%       - 'K',[] [3 x 1 ], [ 3 x nFrame ] calibration parameters
%                       (or [5 x 1 ], [ 5 x nFrame ] when projective)
%                       If given, they won't be optimized upon, you will
%                       need to run bundle adjustment after that
%       - 'KFull',[] [3 x 3 ] or [ 3 x 3 x nFrame ] calibration matrices
%       - 'method', [inf] method for performing SFM (see above for details)
%       - 'onlyErrorFlag', [false] flag indicating if only the error is
%         needed when nFrame==2 && method==Inf && ~isProj
%       - 'nItrSba', [100] number of bundle adjustment iterations
%       - 'tolAffine', [1e-5] tolerance for affine upgrade as explained in
%                          affineUpgrade
%       - 'tolMetric', [1e-5] tolerance for affine upgrade as explained in
%                          metricUprgade
%
% OUTPUTS 1
%  anim      - Animation object
%
% OUTPUTS 2
%  err       - error, acoording to Gold Standard for Affine camera matrix
%
% EXAMPLE
%
% See also
%
% Vincent's Structure From Motion Toolbox      Version 3.1.1
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

[ isCalibrated doAffineUpgrade doMetricUpgrade K KFull method ...
  onlyErrorFlag nItrSBA tolAffine tolMetric ] = getPrmDflt( varargin, ...
  { 'isCalibrated' false 'doAffineUpgrade' false 'doMetricUpgrade' false ...
  'K' [] 'KFull' [] 'method' inf 'onlyErrorFlag' false 'nItrSBA' 100 ...
  'tolAffine' 1e-5 'tolMetric' 1e-5}, 1);

% Remove points that do not have two views with no NaN's
good_point_mask = sum(~any(isnan(W),1),3) >=2;
W = W(:,good_point_mask,:);

nFrame = size(W,3); nPoint = size(W,2);

if doMetricUpgrade; doAffineUpgrade=true; end
if nFrame==2
  doAffineUpgrade=false;
  if ~isCalibrated; doMetricUpgrade=false; end
end

if size(W,1)==3; W = normalizePoint(W,3); end

% create an animation object that wil contain the output
anim=Animation(); anim.isProj=isProj; WOri=W;

% If calibrated, apply inv(K)
if ~isempty(K); anim.K=K; end
if ~isempty(KFull); anim.KFull=KFull; end
KFull=anim.KFull;
if ~isempty(KFull)
  isCalibrated=true;
  doMetricUpgrade=true;
  anim.KFull=[];
  % unapply the internal parameter matrix to the measurements
  if size(KFull,3)==1
    W = normalizePoint(multiTimes(inv(KFull),normalizePoint(W,-3),1.2),3);
  else
    invKFull = zeros(3,3,nFrame);
    for i=1:nFrame; invKFull(:,:,i) = inv(KFull(:,:,i)); end
    W = normalizePoint(multiTimes(inv(KFull),normalizePoint(W,-3),2),3);
  end
end

if isProj
  [P,S] = computeSMFromWProjective( W, method, isCalibrated );
else
  [P,S] = computeSMFromWAffine( W, method, onlyErrorFlag );
end
if onlyErrorFlag; anim=P; return; end

% fill the Animation object with the results
anim.S=S; anim.P=P; anim.W=W;

% perform an affine upgrade if requested
H=[];
if anim.isProj && ~isCalibrated
  % do bundle adjustment
  if nItrSBA > 0
    anim = bundleAdjustment( anim, 'nItr', nItrSBA );
  else
    anim = bundleAdjustment( anim, 'nItr', 20 );
  end

  % modify P so that P(:,:,1)==eye(3,4)
  anim=anim.setFirstPRtToId();

  % if we want to do an affine upgrade
  if doAffineUpgrade
    % only apply the quasi affine upgrade if we are under octave
    % (as Yalmip does not work there)
    if exist('OCTAVE_VERSION','builtin')==5
      warning('Only performing a quasi-affine upgrade under Octave');
      [ HEye, Hqa ] = affineUpgrade(anim);
      H=HEye*Hqa;
    else
      [ HEye, Hqa, pInf ] = affineUpgrade(anim, 'tol', tolAffine);
      % if the camera is calibrated or if we do not do a metric upgrade
      % stop here
      if ~doMetricUpgrade
        % apply H
        H=eye(4); H(4,1:3)=-pInf;
        H=HEye*H;
      end
    end
    if doMetricUpgrade
      % perform a metric upgrade if requested
      if exist('OCTAVE_VERSION','builtin')==5
        H=metricUpgrade(anim, 'isCalibrated', isCalibrated);
      else
        [ H, anim.KFull ]=metricUpgrade(anim, 'pInf', pInf, ...
          'isCalibrated', isCalibrated, 'tol', tolMetric);
      end
    end
  end
end
if anim.isProj && ~isCalibrated
  if doMetricUpgrade
    H=metricUpgrade(anim, 'isCalibrated', isCalibrated, 'method', method);
  end
end

% apply the homography H
if ~isempty(H)
  anim.P=multiTimes(anim.P,H,1);
  anim.S=normalizePoint(H\normalizePoint(anim.S,-4),4);
end

% recover rotations and translations
if doMetricUpgrade && (exist('OCTAVE_VERSION','builtin')~=5 || isCalibrated)
  P=anim.P;
  if ~isCalibrated && ~isempty(anim.KFull)
    if size(anim.KFull,3)==1; P=multiTimes(inv(anim.KFull),P,1.2);
    else P=multiDiv(anim.KFull,P,2);
    end
  end
  R=zeros(3,3,nFrame);
  if anim.isProj
    for i=1:nFrame
      P(:,:,i)=P(:,:,i)/nthroot(det(P(:,1:3,i)),3);
      R(:,:,i)=rotationMatrix(P(:,1:3,i));
    end
    anim.t=reshape(P(:,4,:),3,nFrame);
  else
    for i=1:nFrame
      R(:,:,i)=rotationMatrix(rotationMatrix(P(1:2,1:3,i)));
    end
    anim.t=reshape(P(1:2,4,:),2,nFrame); anim.t(3,:)=0;
  end
  anim.R=R;

  anim=anim.setFirstPRtToId();
end

% update the mask if any
mask=reshape(any(~isnan(W),1),nPoint,nFrame);
if any(mask==0); anim.mask=mask; end

% do a final bundle adjustment
if nItrSBA > 0; anim = bundleAdjustment( anim, 'nItr', nItrSBA ); end

W = anim.W; S = anim.S;
anim.W = zeros(2, length(good_point_mask), nFrame);
anim.S = zeros(3, length(good_point_mask), nFrame);
anim.W = W; anim.S = S;

% Re-assign the KFull from the input arguments if any
if ~isempty(KFull); anim.KFull=KFull; anim.W=WOri; end
end

author: vrabaud


text:
function [Zi, h] = ft_plot_topo(chanX, chanY, dat, varargin)

% FT_PLOT_TOPO interpolates and plots the 2-D spatial topography of the
% potential or field distribution over the head
%
% Use as
%   ft_plot_topo(x, y, val, ...)
%
% Optional arguments should come in key-value pairs and can include
%   'gridscale'     = scalar, number of points along both directions for interpolation (default = 67)
%   'datmask'       = vector of same dimensions as val
%   'mask'          = cell-array with line segments that forms the mask (see FT_PREPARE_LAYOUT)
%   'outline'       = cell-array with line segments that for the outline (see  FT_PREPARE_LAYOUT)
%   'isolines'      = vector with values for isocontour lines (default = [])
%   'interplim'     = string, 'electrodes' or 'mask' (default = 'electrodes')
%   'interpmethod'  = string, 'nearest', 'linear', 'natural', 'cubic' or 'v4' (default = 'v4')
%   'style'         = can be 'surf', 'iso', 'isofill', 'surfiso', 'imsat', 'imsatiso', 'colormix'
%   'clim'          = [min max], limits for color scaling
%   'shading'       = string, 'none', 'flat', 'interp' (default = 'flat')
%   'parent'        = handle which is set as the parent for all plots
%   'tag'           = string, the name assigned to the object. All tags with the same name can be deleted in a figure, without deleting other parts of the figure.
%   'box'           = draw a box around the local axes, can be 'yes' or 'no'
%
% It is possible to plot the object in a local pseudo-axis (c.f. subplot), which is specfied as follows
%   'hpos'          = horizontal position of the lower left corner of the local axes
%   'vpos'          = vertical position of the lower left corner of the local axes
%   'width'         = width of the local axes
%   'height'        = height of the local axes
%   'hlim'          = horizontal scaling limits within the local axes
%   'vlim'          = vertical scaling limits within the local axes
%
% See also FT_PLOT_TOPO3D, FT_PLOT_LAYOUT, FT_TOPOPLOTER, FT_TOPOPLOTTFR

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

% these are for speeding up the plotting on subsequent calls
persistent previous_argin previous_maskimage

% get the optional input arguments
hpos          = ft_getopt(varargin, 'hpos',         0);
vpos          = ft_getopt(varargin, 'vpos',         0);
width         = ft_getopt(varargin, 'width',        []);
height        = ft_getopt(varargin, 'height',       []);
gridscale     = ft_getopt(varargin, 'gridscale',    67); % 67 in original
shading       = ft_getopt(varargin, 'shading',      'flat');
interplim     = ft_getopt(varargin, 'interplim',    'electrodes');
interpmethod  = ft_getopt(varargin, 'interpmethod', 'v4');
style         = ft_getopt(varargin, 'style',        'surfiso'); % can be 'surf', 'iso', 'isofill', 'surfiso', 'imsat', 'imsatiso', 'colormix'
tag           = ft_getopt(varargin, 'tag',          '');
isolines      = ft_getopt(varargin, 'isolines');
datmask       = ft_getopt(varargin, 'datmask');
mask          = ft_getopt(varargin, 'mask');
outline       = ft_getopt(varargin, 'outline');
clim          = ft_getopt(varargin, 'clim', []);
parent        = ft_getopt(varargin, 'parent', []);
box           = ft_getopt(varargin, 'box', false);

% convert the yes/no strings into boolean values
box = istrue(box);

% check for nans in the data, they can be still left incase people want to mask non channels.
if any(isnan(dat))
  ft_warning('the data passed to ft_plot_topo contains NaNs, these channels will be removed from the data to prevent interpolation errors, but will remain in the mask');
  flagNaN = true;
else
  flagNaN = false;
end
NaNind = isnan(dat);

% everything is added to the current figure
holdflag = ishold;
if ~holdflag
  hold on
end

% layout units can be arbitrary (e.g. pixels for .mat files)
% so we need to compute the right scaling factor
% create a matrix with all coordinates
% from positions, mask, and outline
allCoords = [chanX(:) chanY(:)];
if ~isempty(mask)
  for k = 1:numel(mask)
    allCoords = [allCoords; mask{k}];
  end
end
if ~isempty(outline)
  for k = 1:numel(outline)
    allCoords = [allCoords; outline{k}];
  end
end

naturalWidth = (max(allCoords(:, 1))-min(allCoords(:, 1)));
naturalHeight = (max(allCoords(:, 2))-min(allCoords(:, 2)));

if isempty(width) && isempty(height)
  xScaling = 1;
  yScaling = 1;
elseif isempty(width) && ~isempty(height)
  % height specified, auto-compute width while maintaining aspect ratio
  yScaling = height/naturalHeight;
  xScaling = yScaling;
elseif ~isempty(width) && isempty(height)
  % width specified, auto-compute height while maintaining aspect ratio
  xScaling = width/naturalWidth;
  yScaling = xScaling;
else
  % both width and height specified
  xScaling = width/naturalWidth;
  yScaling = height/naturalHeight;
end

% keep original channel positions
chanXorg = chanX;
chanYorg = chanY;
chanX = chanX(:) * xScaling + hpos;
chanY = chanY(:) * yScaling + vpos;

if strcmp(interplim, 'electrodes')
  hlim = [min(chanX) max(chanX)];
  vlim = [min(chanY) max(chanY)];
elseif (strcmp(interplim, 'mask') || strcmp(interplim, 'mask_individual')) && ~isempty(mask)
  hlim = [inf -inf];
  vlim = [inf -inf];
  for i=1:length(mask)
    hlim = [min([hlim(1); mask{i}(:, 1)*xScaling+hpos]) max([hlim(2); mask{i}(:, 1)*xScaling+hpos])];
    vlim = [min([vlim(1); mask{i}(:, 2)*yScaling+vpos]) max([vlim(2); mask{i}(:, 2)*yScaling+vpos])];
  end
else
  hlim = [min(chanX) max(chanX)];
  vlim = [min(chanY) max(chanY)];
end

% check if all mask point are inside the limits otherwise redefine mask
newpoints = [];
if length(mask)==1
  % which channels are outside
  inside  = inside_contour([chanX chanY], mask{1});
  outside = ~inside;
  newpoints = [chanX(outside) chanY(outside)];
end

% try to speed up the preparation of the mask on subsequent calls
%current_argin = {chanX, chanY, gridscale, mask, datmask}; % old: to be plotted channel positions must be the same
current_argin = {chanXorg, chanYorg, gridscale, mask, datmask, interplim}; % new: unscaled channel positions must be the same over calls
if isequal(current_argin, previous_argin)
  % don't construct the binary image, but reuse it from the previous call
  maskimage = previous_maskimage;
elseif ~isempty(mask)
  % convert the mask into an indexed image
  maskimage = zeros(gridscale); %false(gridscale);
  %hlim      = [min(chanX) max(chanX)];
  %vlim      = [min(chanY) max(chanY)];
  xi        = linspace(hlim(1), hlim(2), gridscale);   % x-axis for interpolation (row vector)
  yi        = linspace(vlim(1), vlim(2), gridscale);   % y-axis for interpolation (row vector)
  [Xi, Yi]   = meshgrid(xi', yi);
  if ~isempty(newpoints) && (hpos == 0 || vpos == 0)
    ft_warning('Some points fall outside the outline, please consider using another layout')
    % FIXME: I am not sure about it, to be tested!
    %     tmp = [mask{1};newpoints];
    %     indx = convhull(tmp(:, 1), tmp(:, 2));
    %     mask{1} = tmp(indx, :);
    % NOTE: if you set hpos and/or vpos, newpoints is not empty, but nothing
    % needs to be fixed (this fixme screws up things, then)
  end
  for i=1:length(mask)
    mask{i}(:, 1) = mask{i}(:, 1)*xScaling+hpos;
    mask{i}(:, 2) = mask{i}(:, 2)*yScaling+vpos;
    mask{i}(end+1, :) = mask{i}(1, :);                   % force them to be closed
    maskimage(inside_contour([Xi(:) Yi(:)], mask{i})) = i; %true;
  end

else
  maskimage = [];
end

% adjust maskimage to also mask channels as specified in maskdat
if ~isempty(datmask)
  xi           = linspace(hlim(1), hlim(2), gridscale);   % x-axis for interpolation (row vector)
  yi           = linspace(vlim(1), vlim(2), gridscale);   % y-axis for interpolation (row vector)
  maskimagetmp = griddata(chanX', chanY, double(datmask), xi', yi, 'nearest'); % interpolate the mask data
  if isempty(maskimage)
    maskimage = maskimagetmp;
  else
    maskimage = (maskimage + maskimagetmp) > 1.01;
  end
end

% take out NaN channels if interpmethod does not work with NaNs
if flagNaN && strcmp(interpmethod, 'v4')
  dat(NaNind) = [];
  chanX(NaNind) = [];
  chanY(NaNind) = [];
end

% convert chanX, chanY and dat to double
dat = double(dat);
chanX = double(chanX);
chanY = double(chanY);

%interpolate data
xi         = linspace(hlim(1), hlim(2), gridscale);       % x-axis for interpolation (row vector)
yi         = linspace(vlim(1), vlim(2), gridscale);       % y-axis for interpolation (row vector)

if ~ft_platform_supports('griddata-vector-input')
  % in GNU Octave, griddata does not support vector
  % positions; make a grid to get the locations in vector form
  [xi,yi]=meshgrid(xi,yi);
  xi=xi';
end

if ~isempty(maskimage) && strcmp(interplim, 'mask_individual')
  % do the interpolation for each set of electrodes within a mask, useful
  % for ECoG data with multiple grids, to avoid cross talk
  Zi = zeros(size(maskimage));
  for i=1:max(maskimage(:))
    chansel = inside_contour([chanX chanY], mask{i});
    [Xi, Yi, tmpZi]  = griddata(chanX(chansel)', chanY(chansel), dat(chansel), xi', yi, interpmethod);
    Zi(maskimage==i) = tmpZi(maskimage==i);
  end
else
  [Xi, Yi, Zi] = griddata(chanX', chanY, dat, xi', yi, interpmethod); % interpolate the topographic data
end

if ~isempty(maskimage)
  % apply mask to the data to hide parts of the interpolated data (outside the circle) and channels that were specified to be masked
  % this combines the input options mask and maskdat
  Zi(maskimage==0) = NaN;
end

% The topography should be plotted prior to the isolines to ensure that it is exported correctly, see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2496
if strcmp(style, 'surf') || strcmp(style, 'surfiso')

  deltax = xi(2)-xi(1); % length of grid entry
  deltay = yi(2)-yi(1); % length of grid entry
  h = surf(Xi-deltax/2, Yi-deltay/2, zeros(size(Zi)), Zi, 'EdgeColor', 'none', 'FaceColor', shading);
  set(h, 'tag', tag);
  if ~isempty(parent)
    set(h, 'Parent', parent);
  end

  % this is an alternative to using NaN values to hide parts of the interpolated data (outside the circle), but fails in eps files
  % the NaN solution is more robust and works for all renderers and file formats
  %   if ~isempty(maskimage)
  %     set(h, 'facealpha', 'flat');
  %     set(h, 'alphadatamapping', 'scaled');
  %     set(h, 'alphadata', maskimage);
  %   end

elseif strcmp(style, 'imsat') || strcmp(style, 'imsatiso')
  % Plot the surface in an alternate style (using imagesc and saturation masking) so that it can be nicely saved to a vectorized format

  % set mask and cdat, and check for clim
  if isempty(clim)
    ft_error('clim is required for style = ''imsat'' or style = ''imsatiso''')
  end

  cmap    = get(gcf, 'colormap');
  rgbcdat = cdat2rgb(Zi, cmap, clim, maskimage~=0);

  h = imagesc(xi, yi, rgbcdat, clim);
  set(h, 'tag', tag);

elseif strcmp(style, 'colormix')
  % Plot the surface in an alternate style (using imagesc and saturation masking) so that it can be nicely saved to a vectorized format

  % set mask and cdat, and check for clim
  if isempty(clim)
    error('clim is required for style = ''colormix''')
  end

  if ~isempty(datmask)
    % use maskimagetmp in combination with maskimage, maskimagetmp is
    % scaled between 0 and 1
    maskimagetmp = maskimagetmp./max(maskimagetmp(:));
    Zmask        = double(maskimage~=0);
    Zmask(Zmask>0) = maskimagetmp(Zmask>0);
  else
    Zmask        = double(maskimage~=0);
  end

  cmap    = get(gcf, 'colormap');
  rgbcdat = bg_rgba2rgb([1 1 1], Zi, cmap, clim, Zmask, 'rampup', [0 1]);

  h = imagesc(xi, yi, rgbcdat, clim);
  set(h,'tag',tag);

end

% Plot the outline of the head, ears and nose
for i=1:length(outline)
  xval = outline{i}(:, 1) * xScaling  + hpos;
  yval = outline{i}(:, 2) * yScaling + vpos;
  ft_plot_vector(xval, yval, 'Color', 'k', 'LineWidth', 2, 'tag', tag, 'parent', parent);
end

% Create isolines
if strcmp(style, 'iso') || strcmp(style, 'surfiso') ||  strcmp(style, 'imsatiso')
  if ~isempty(isolines)
    [cont, h] = contour(Xi, Yi, Zi, isolines, 'k');
    set(h, 'tag', tag);
    if ~isempty(parent)
      set(h, 'Parent', parent);
    end
  end
end

% Plot filled contours
if strcmp(style, 'isofill') && ~isempty(isolines)
  [cont, h] = contourf(Xi, Yi, Zi, isolines, 'k');
  set(h, 'tag', tag);
  if ~isempty(parent)
    set(h, 'Parent', parent);
  end
end

if box
  % this plots a box around the original hpos/vpos with appropriate width/height
  boxposition = zeros(1,4);
  boxposition(1) = hpos - width/2;
  boxposition(2) = hpos + width/2;
  boxposition(3) = vpos - height/2;
  boxposition(4) = vpos + height/2;
  ft_plot_box(boxposition);
end

% apply clim if it was given
if ~isempty(clim)
  caxis(clim)
end

% remember the current input arguments, so that they can be
% reused on a subsequent call in case the same input argument is given
previous_argin     = current_argin;
previous_maskimage = maskimage;

if ~holdflag
  hold off
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function plot_savePC(hFig, obj, idxPC)

% Copy object
hFig.UserData.PC{idxPC} = obj.copy;

% Keep only active points to save memory
if sum(hFig.UserData.PC{idxPC}.act) < hFig.UserData.PC{idxPC}.noPoints

    % Deactivate logging
    originalLogLevel = msg('O', 'GetLogLevel');
    msg('O', 'SetLogLevel', 'off');

    % Reconstruct
    hFig.UserData.PC{idxPC}.reconstruct;

    % Reactivate logging
    msg('O', 'SetLogLevel', originalLogLevel);

end

% Fields saved in user property of point cloud (hFig.UserData.PC{idxPC}.U):
% PC.U.hPlot                   handles of plotted points
% PC.U.hMenu                   handle of menu
% PC.U.Visible                 logical value that defines if PC is visible or not
% PC.U.p                       plot parameters
% PC.U.percentOfVisiblePoints  percent of points visible within area defined by the parameter 'Limits'
% PC.U.plotNormals             plot normals? (true or false)
% PC.U.hPlotNormals            handle of normals

% Initialize plot handle
hFig.UserData.PC{idxPC}.U.hPlot = [];

% Initialize menu handle
hFig.UserData.PC{idxPC}.U.hMenu = [];

% Visible
hFig.UserData.PC{idxPC}.U.Visible = true;

% Default parameters
hFig.UserData.PC{idxPC}.U.p.Color        = 'A.z';
hFig.UserData.PC{idxPC}.U.p.MarkerSize   = 1;
hFig.UserData.PC{idxPC}.U.p.MaxPoints    = 1e6;
hFig.UserData.PC{idxPC}.U.p.ColormapName = 'parula';
hFig.UserData.PC{idxPC}.U.p.CAxisLim     = [];
hFig.UserData.PC{idxPC}.U.p.Limits       = [-Inf Inf; -Inf Inf; -Inf Inf];

% Plot normals
hFig.UserData.PC{idxPC}.U.plotNormals = false;

% Initialize plot handle for normals
hFig.UserData.PC{idxPC}.U.hPlotNormals = [];

end

author: pglira


text:
% Stress recovery in plane stress analysis of plates
% Plane stress analysis of a thin plate under tension at its extremes
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Warning : On running this the workspace memory will be deleted. Save any
% if any data present before running the code !!
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%--------------------------------------------------------------------------
% Code written by : Siva Srinivas Kolukula                                |
%                   Senior Research Fellow                                |
%                   Structural Mechanics Laboratory                       |
%                   Indira Gandhi Center for Atomic Research              |
%                   India                                                 |
% E-mail : allwayzitzme@gmail.com                                         |
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%
% Variable descriptions
%   k = element matrix for stiffness
%   f = element vector
%   stiffness = system matrix
%   force = system vector
%   displacement = system nodal displacement vector
%   coordinates = coordinate values of each node
%   nodes = nodal connectivity of each element
%   index = a vector containing system dofs associated with each element
%   gausspoint = matrix containing sampling points for bending term
%   gaussweight = matrix containing weighting coefficients for bending term
%   bcdof = a vector containing dofs associated with boundary conditions
%   bcval = a vector containing boundary condition values associated with
%           the dofs in 'bcdof'
%   B = matrix for kinematic equation for plane stress
%   D = matrix for material property for plane stress
%   stressGP = stress at the Gauss points
%----------------------------------------------------------------------------

%--------------------------------------------------------------------------
%  input data
%--------------------------------------------------------------------------
clear
clc
disp('Please wait Programme is under Run')
%--------------------------------------------------------------------------
% Input data for nodal coordinate values
%--------------------------------------------------------------------------
load coordinates.dat ;
%--------------------------------------------------------------------------
% Input data for nodal connectivity for each element
%--------------------------------------------------------------------------
load nodes.dat ;

nel = length(nodes) ;                  % number of elements
nnel=4;                                % number of nodes per element
ndof=2;                                % number of dofs per node (UX,UY)
nnode = length(coordinates) ;          % total number of nodes in system
sdof=nnode*ndof;                       % total system dofs
edof=nnel*ndof;                        % degrees of freedom per element

% Units are in SI system
a = 1 ;                           % Length of the plate (along X-axes)
b = 1 ;                           % Length of the plate (along Y-axes)
elementsalongX = 10 ;             % Number of elements along X-axes
elementsalongY = 10 ;

E = 2.1*10^11;                      % Youngs modulus
nu = 0.3;                           % Poisson's ratio
t = 0.0254;                         % plate thickness
rho = 7840. ;                       % Density of the plate

nglx = 2; ngly = 2;         % 2x2 Gauss-Legendre quadrature
nglxy=nglx*ngly;            % number of sampling points per element

%--------------------------------------------------------------------------
% Input data for boundary conditions
%--------------------------------------------------------------------------
% (0,0) and (1,0) are fixed

  bcdof = [ 1 2 241 242] ;
  bcval = zeros(1,length(bcdof)) ;
%--------------------------------------------------------------------------
%  initialization of matrices and vectors
%--------------------------------------------------------------------------

force = zeros(sdof,1);                % system force vector
stiffness = zeros(sdof,sdof);         % system stiffness matrix
displacement = zeros(sdof,1);         % system displacement vector
eldepl = zeros(edof,1) ;              % element displacement vector
index = zeros(edof,1);                % index vector
B = zeros(3,edof);              % kinematic matrix for bending
D = zeros(3,3);                 % constitutive matrix for bending

%--------------------------------------------------------------------------
% force vector
%--------------------------------------------------------------------------
P = 1e5 ;       % Load
rightedge = find(coordinates(:,1)==a);
rightdof = 2*rightedge-ones(length(rightedge),1);
force(rightdof) = -P*b/(elementsalongY+1) ;
leftedge = find(coordinates(:,1)==0);
leftdof = 2*leftedge-ones(length(leftedge),1) ;
force(leftdof) = P*b/(elementsalongY+1) ;

%--------------------------------------------------------------------------
%  computation of element matrices and vectors and their assembly
%--------------------------------------------------------------------------
[Gausspoint,Gaussweight]=GaussQuadrature(nglx);     % sampling points & weights
D = E/(1-nu^2)*[1 nu 0 ; nu 1 0; 0 0 (1-nu)/2] ;    % Constituent Matrix for Plane stress

for iel=1:nel           % loop for the total number of elements

for i=1:nnel
nd(i)=nodes(iel,i);         % extract connected node for (iel)-th element
xx(i)=coordinates(nd(i),1);    % extract x value of the node
yy(i)=coordinates(nd(i),2);    % extract y value of the node
end

k = zeros(edof,edof);        % initialization of stiffness matrix

%--------------------------------------------------------------------------
%  numerical integration for stiffness matrix
%--------------------------------------------------------------------------

for intx=1:nglx
xi = Gausspoint(intx,1);                  % sampling point in x-axis
wtx = Gaussweight(intx,1);               % weight in x-axis
for inty=1:ngly
eta = Gausspoint(inty,1);                  % sampling point in y-axis
wty = Gaussweight(inty,1) ;              % weight in y-axis

[shape,dhdr,dhds] = shapefunctions(xi,eta);     % compute shape functions and
                                    % derivatives at sampling point
jacobian = Jacobian(nnel,dhdr,dhds,xx,yy);  % compute Jacobian

detjacob=det(jacobian);                 % determinant of Jacobian
invjacob=inv(jacobian);                 % inverse of Jacobian matrix

[dhdx,dhdy]=shapefunctionderivatives(nnel,dhdr,dhds,invjacob); % derivatives w.r.t.
                                               % physical coordinate

B=fekineps(nnel,dhdx,dhdy);          % kinematic matrix for stiffness

%--------------------------------------------------------------------------
%  compute element stiffness matrix
%--------------------------------------------------------------------------

k = k+B'*D*B*wtx*wty*detjacob;

end
end                      % end of numerical integration loop for bending term

index = elementdof(nd,nnel,ndof); % extract system dofs associated with element

stiffness = assemble(stiffness,k,index);    % assemble element stiffness matrices

end

%--------------------------------------------------------------------------
%   apply boundary conditions
%--------------------------------------------------------------------------

[stiffness,force] = constraints(stiffness,force,bcdof,bcval);

%--------------------------------------------------------------------------
% Solve the matrix equation
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Solve the matrix equation
%--------------------------------------------------------------------------
displacement = stiffness\force ;

%--------------------------------------------------------------------------
% Stress Recovery
%--------------------------------------------------------------------------
% Calculaing the stress at the Gauss points
[vonmises,stressGP,stressEXP,stressPR,stressG] = StressRecovery(D,displacement);

% Calculating the Element nodal stresses

% Nodal Averaging
type = 'average' ;
% type = 'sum' ;
stressEXPavg = NodalAveraging(stressEXP,type);
stressPRavg = NodalAveraging(stressPR,type) ;

%#### OUTPUTS ####

display('=======SYSTEM MAIN PROPERTIES (isotropic material)=======');
fprintf('       Elasticity  constant  % .2e N/m**2 \n',E);
fprintf('       Poission ratio        % .5f         \n',nu);
fprintf('       thickness of plate    % .5f m       \n',t);
fprintf('       length of plate       % .5f m       \n',a);
fprintf('       breadth of plate      % .5f m       \n',b);
fprintf(' \n ')

display('======SYSTEM GEOMETRICALY PROPERTIES======');
fprintf('        System total element no   % .f \n  ',nel);
fprintf('      System total nodes        % .f \n  ',nnode);
fprintf('      System degree of freedom  % .f \n  ',sdof);

fprintf(' \n ')
% Output of displacements
mytable(displacement,1);

% Output of Stresses at Gauss points
display('======STRESS VALUES AT THE GAUSS POINTS=======')
mytable(stressGP,3);

% Output Von Mises stress
mytable(vonmises,2) ;

% Output of element nodal stresses
display('====ELEMENT NODAL STRESS USING EXTRAPOLATION====')
mytable(stressEXP,3) ;
display('====ELEMENT NODAL STRESS USING PATCH RECOVERY====')
mytable(stressPR,3) ;

% Output of Nodal averaged stresses
display('====NODAL AVERAGED STRESSES OF EXTRAPOLATION====')
mytable(stressEXPavg,4);
display('====NODAL AVERAGED STRESSES OF PATCH RECOVERY====')
mytable(stressPRavg,4);

author: Sable

repo: mcbench-benchmarks


text:
classdef MOEADDE_F2 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing MOEA/D-DE

%------------------------------- Reference --------------------------------
% H. Li and Q. Zhang, Multiobjective optimization problems with complicated
% Pareto sets, MOEA/D and NSGA-II, IEEE Transactions on Evolutionary
% Computation, 2009, 13(2): 284-302.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = [0,-ones(1,obj.D-1)];
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            J1 = 3 : 2 : obj.D;
            J2 = 2 : 2 : obj.D;
            PopObj(:,1) = X(:,1)         + 2*mean((X(:,J1)-sin(repmat(6*pi*X(:,1),1,length(J1))+repmat(J1*pi/obj.D,size(X,1),1))).^2,2);
            PopObj(:,2) = 1-sqrt(X(:,1)) + 2*mean((X(:,J2)-sin(repmat(6*pi*X(:,1),1,length(J2))+repmat(J2*pi/obj.D,size(X,1),1))).^2,2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^0.5;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% Compute error statistics on groundtruth vs NovAtel
distance = ((ground_truth(:,1)-gps_x).^2 + (ground_truth(:,2)-gps_y).^2).^0.5;
raw_rms = rms(distance);
raw_max = max(distance);
fprintf('\nRMS 2D Raw: %0.4fm\n',raw_rms)
fprintf('Max 2D Raw: %0.4fm\n',raw_max)

%% Compute error statistics on groundtruth vs filter output
xyz = out_profile(:,2:4);
llh = ecef2lla(xyz);
[x,y] = deg2utm(llh(:,1),llh(:,2));
x = x-min_x;
y = y-min_y;
h = -llh(:,3);

distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2).^0.5;
filter_rms = rms(distance);
filter_max = max(distance);
fprintf('RMS 2D Filter: %0.4fm, %0.2f%% change\n',filter_rms,(raw_rms-filter_rms)/raw_rms*100);
fprintf('Max 2D Filter: %0.4fm, %0.2f%% change\n',filter_max,(raw_max-filter_max)/raw_max*100);

%% 3D Stats
distance = ((ground_truth(:,1)-gps_x).^2 + (ground_truth(:,2)-gps_y).^2 + (ground_truth(:,3)-gps_h).^2).^0.5;
raw_rms = rms(distance);
raw_max = max(distance);
fprintf('\nRMS 3D Raw: %0.4fm\n',raw_rms)
fprintf('Max 3D Raw: %0.4fm\n',raw_max)
distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2 + (ground_truth_full(:,3)-h).^2).^0.5;
filter_rms = rms(distance);
filter_max = max(distance);
fprintf('RMS 3D Filter: %0.4fm, %0.2f%% change\n',filter_rms,(raw_rms-filter_rms)/raw_rms*100);
fprintf('Max 3D Filter: %0.4fm, %0.2f%% change\n',filter_max,(raw_max-filter_max)/raw_max*100);

%% Rotation stats

author: awerries

repo: kalman-localization


text:
function demoami
% DEMOAMI demos function ami for bivariate time series data
%
% USAGE:
%       demoami
%
% INPUT:
%
% OUTPUT:
%      amis:   vector of average mutual information for time lags of 0:nLags
%     corrs:   vector of correlation (or autocorrelation for univariate) for
%              time lags of 0:nLags
%
% EXAMPLES:
%
% See also AMI, PROBXY, RHIST

% Copyright 2004-2005 by Durga Lal Shrestha.
% eMail: durgals@hotmail.com
% $Date: 2005/06/27
% $Revision: 1.1.0 $ $Date: 2005/07/01 $

% ***********************************************************************
clear all
close all
clc
mydata = load('data.txt')';
lag = 25;
nBins = [15 15];
[iy ry] = ami(mydata,nBins,lag);


text:
% Description:
%
%     Return a list of unique values.  A must be a cell array of vectors of
%     numeric data.  B will be sorted in the order of the first occurrance of
%     each unique value in A.
%
% Syntax:
%
%     [ B, I, J ] = mUnique(A)
%
% Inputs:
%
%     A - [ 1 x N ] (cell)
%
% Outputs:
%
%     B - [ 1 x U ] (cell)
%     I - [ 1 x U ] (int)
%     J - [ 1 x N ] (int)
%
% Details:
%
% Examples:
%
% Notes:
%
% Author(s):
%
%     William Gruner (williamgruner@gmail.com)
%
% References:
%
% Acknowledgements:
%
%     Many thanks to Dr. Erik Erhardt and Dr. Elena Allen of the Mind Research
%     Network (www.mrn.org) for their continued collaboration.
%
% Version:
%
%     $Author: williamgruner $
%     $Date: 2010-04-01 11:39:09 -0600 (Thu, 01 Apr 2010) $
%     $Revision: 482 $

function [ B, I, J ] = mUnique(A)

    temp = cell(0);

    for i = 1 : length(A)
        temp{i} = num2str(A{i}, '%d ');
    end

    [ B, I, J ] = unique(temp);

    I = sort(I);
    B = A(I);


text:
function d=jensen_shannon_divergence(XI,XJ)
  % Implementation of the Jensen-Shannon Divergence to use with pdist
  % (cf. "The Earth Movers' Distance as a Metric for Image Retrieval",
  %      Y. Rubner, C. Tomasi, L.J. Guibas, 2000)
  %
  % @author: B. Schauerte
  % @date:   2009
  % @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/

  % Copyright 2009 B. Schauerte. All rights reserved.
  %
  % Redistribution and use in source and binary forms, with or without
  % modification, are permitted provided that the following conditions are
  % met:
  %
  %    1. Redistributions of source code must retain the above copyright
  %       notice, this list of conditions and the following disclaimer.
  %
  %    2. Redistributions in binary form must reproduce the above copyright
  %       notice, this list of conditions and the following disclaimer in
  %       the documentation and/or other materials provided with the
  %       distribution.
  %
  % THIS SOFTWARE IS PROVIDED BY B. SCHAUERTE ''AS IS'' AND ANY EXPRESS OR
  % IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  % WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  % DISCLAIMED. IN NO EVENT SHALL B. SCHAUERTE OR CONTRIBUTORS BE LIABLE
  % FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  % CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  % SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
  % BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
  % WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
  % OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
  % ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  %
  % The views and conclusions contained in the software and documentation
  % are those of the authors and should not be interpreted as representing
  % official policies, either expressed or implied, of B. Schauerte.

  d=jeffrey_divergence(XI,XJ);
  d=d / 2;

author: opencobra

repo: cobratoolbox