text:
function [ candidates, scores ] = sample_bing_windows( im, num_samples)
%SAMPLE_BING_WINDOWS Will generate equaly distributed windows in space,
%following Bing sizes
%   Bing uses 29 specific sizes, this method spread this sizes homogenously
%   inside the image
scores = [];
im_wh = [size(im, 2), size(im, 1)];

original_bing_window_sizes = [[512 512 ]; [256, 512 ]; [128, 512 ]; [64, 512 ]; ...
    [512, 256 ]; [256, 256 ]; [128, 256 ]; [64, 256 ]; [32, 256 ]; ...
    [512, 128 ]; [256, 128 ]; [128, 128 ]; [64, 128 ]; [32, 128 ]; [16, 128 ]; ...
    [512, 64 ]; [256, 64 ]; [128, 64 ]; [64, 64 ]; [32, 64 ]; [16, 64 ]; ...
    [128, 32 ]; [64, 32 ]; [32, 32 ]; [16, 32 ];  ...
    [64, 16 ]; [32, 16 ]; [16, 16 ]];

%original_bing_window_sizes = [[64, 128]; [32, 128 ]];

original_num_window_sizes = size(original_bing_window_sizes, 1);

bing_window_sizes = [];

% we filter window sizes to fit inside the image
for i=1:size(original_bing_window_sizes, 1),
    window_wh = original_bing_window_sizes(i, :);

    if(sum(window_wh < im_wh) == 2),
        bing_window_sizes = [bing_window_sizes; window_wh];
    else
        % the window is disregarded
    end

end

num_window_sizes = size(bing_window_sizes, 1);

assert(num_samples > num_window_sizes);

candidates = [];
if num_window_sizes ~= original_num_window_sizes,
    % we add one candidate that covers the whole image size
   candidates =  [0, 0, im_wh];
end

use_v0 = false;

if use_v0,
    % will have as many splits in x and y axis (non-square strides)

    sqrt_num_samples_size = floor(sqrt(num_samples / num_window_sizes));
    samples_per_size = sqrt_num_samples_size * sqrt_num_samples_size;
    assert(samples_per_size > 3);
    extra_samples = num_samples - (samples_per_size * num_window_sizes) - size(candidates, 1);
    delta_samples = ((sqrt_num_samples_size + 1) * (sqrt_num_samples_size + 1)) - samples_per_size;

    divisions_per_size = repmat(sqrt_num_samples_size, 1, num_window_sizes);

    start_extra_size_index = 6; % because we like 256x256 (and smaller)
    for i = start_extra_size_index:num_window_sizes,
        if extra_samples > delta_samples,
            divisions_per_size(i) = divisions_per_size(i) + 1;
            extra_samples = extra_samples - delta_samples;
        else
            fprintf('Added %i extra sizes\n', i - 1);
            break;
        end
    end

    for i = 1:num_window_sizes,

        num_divisions = divisions_per_size(i);
        window_wh = bing_window_sizes(i, :);
        x1 = 1;
        y1 = 1;
        x2 = im_wh(1) - window_wh(1);
        y2 = im_wh(2) - window_wh(2);

        assert(x2 > 1);
        assert(y2 > 1);

        [xx, yy] = meshgrid(linspace(x1, x2, num_divisions), ...
            linspace(y1, y2, num_divisions));
        top_left_xy = [xx(:) yy(:)];
        for j = 1:size(top_left_xy, 1),
            xy = top_left_xy(j, :);
            window = [xy, xy + window_wh];
            candidates = [candidates; window];
        end

    end

else
    % v1
    % will use square strides

    num_samples_per_size = floor(num_samples / num_window_sizes);

    stride_per_size = zeros(1, num_window_sizes);
    total_placed_samples = size(candidates, 1);
    for i = 1:num_window_sizes,

        window_wh = bing_window_sizes(i, :);
        x2 = im_wh(1) - window_wh(1);
        y2 = im_wh(2) - window_wh(2);

        assert(x2 > 0);
        assert(y2 > 0);

        block_area = (x2 * y2) / num_samples_per_size;
        stride = sqrt(block_area);
        assert(stride > 0);

        num_samples_placed = compute_num_samples_placed(im_wh, window_wh, stride);

        while(num_samples_placed > num_samples_per_size)
            stride = stride + 1; % larger stride, less placements
            num_samples_placed = compute_num_samples_placed(im_wh, window_wh, stride);
        end

        assert(num_samples_placed <= num_samples_per_size);

        total_placed_samples = total_placed_samples + num_samples_placed;
        assert(total_placed_samples <= num_samples);

        stride_per_size(i) = stride;
    end

    assert(total_placed_samples <= num_samples);

    sqrt_num_samples_size = floor(sqrt(num_samples / num_window_sizes));
    samples_per_size = sqrt_num_samples_size * sqrt_num_samples_size;
%     assert(samples_per_size > 3);
    assert(samples_per_size > 0);
    extra_samples = num_samples - total_placed_samples;
    assert(extra_samples >= 0);

    start_extra_size_index = 6; % because we like 256x256 (and smaller)
    for i = start_extra_size_index:num_window_sizes,

        window_wh = bing_window_sizes(i, :);
        stride = stride_per_size(i);

        num_samples_placed = compute_num_samples_placed(im_wh, window_wh, stride);

        new_stride = stride * 0.75;
        %new_stride = stride - 1;
        new_num_samples_placed = compute_num_samples_placed(im_wh, window_wh, new_stride);

        delta_samples = new_num_samples_placed - num_samples_placed;

        if extra_samples > delta_samples,
            stride_per_size(i) = new_stride;
            total_placed_samples = total_placed_samples + delta_samples;
            extra_samples = extra_samples - delta_samples;
        else
            fprintf('Added %i extra sizes\n', i - start_extra_size_index);
            break;
        end
    end

    if extra_samples > 0,
        fprintf('%i extra_samples remaining\n', extra_samples);
    end

    for i = 1:num_window_sizes,

        window_wh = bing_window_sizes(i, :);
        stride = stride_per_size(i);
        x1 = 1;
        y1 = 1;
        x2 = im_wh(1) - window_wh(1);
        y2 = im_wh(2) - window_wh(2);

        x_dots = x1 + (mod(x2 - x1, stride) / 2): stride : x2;
        y_dots = y1 + (mod(y2 - y1, stride) / 2): stride : y2;

        assert(x2 > 0);
        assert(y2 > 0);

        [xx, yy] = meshgrid(x_dots, y_dots);
        top_left_xy = [xx(:) yy(:)];
        for j = 1:size(top_left_xy, 1),
            xy = top_left_xy(j, :);
            window = [xy, xy + window_wh];
            candidates = [candidates; window];
        end

    end

    assert(size(candidates, 1) == total_placed_samples);
end

assert(size(candidates, 1) <= num_samples);
end

function [ret] = compute_num_samples_placed(im_wh, window_wh, stride)
x1 = 1;
y1 = 1;
x2 = im_wh(1) - window_wh(1);
y2 = im_wh(2) - window_wh(2);

x_dots = x1 + (mod(x2 - x1, stride) / 2): stride : x2;
y_dots = y1 + (mod(y2 - y1, stride) / 2): stride : y2;
num_samples_placed = size(x_dots, 2) * size(y_dots, 2);
ret = num_samples_placed;
end

author: hosang

repo: detection-proposals


text:
% ----- Global evaluation parameters -----
% Range of binning factors used for this evaluation.
binningFactors = [2, 3, 4];
% Range of compression levels (index 1 means uncoded; indices 2-5 mean
% coded with the corresponding QP setting using H.265/HEVC coder.
compressions   = [NaN, 10, 20, 30, 40];
% Range of number of frames for SR at the different binning factors.
numberOfFrames = [5, ...    % Number of frames for binning factor 2
                  11, ...   % Number of frames for binning factor 3
                  17];      % Number of frames for binning factor 4
% Index of the first reference frame for sliding window processing.
startReferenceFrame = 9;

% ----- Settings considered for this evaluation -----
% This is the index of the SR method that is evaluated.
sr_method           = 0:length(SRMethods);
% This is the index of the binning factor that is evaluated.
binning_val         = 1:3;
% This is the index of the compression setting.
compress_val        = 1:5;
% This is the index that indicates the number of input frames.
numberOfFrames_val  = 1;
% This is the index of the scene according to datasets.mat.
scenes_val          = 1:14;
% This is the index of the sliding window that is evaluated.
sliding_val         = 1;

author: thomas-koehler

repo: SupER


text:
function bst_which(filename, exploreMode)
% BST_WHICH:  Locates a file on disk, using available system tools
%
% USAGE:  bst_which(filename, 'explorer')
%         bst_which(filename, 'terminal')

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

% Parse inputs
isTerminal = strcmpi(exploreMode, 'terminal');

% Check file type
fileType = file_gettype(filename);
switch (fileType)
    case 'link'
        filename = file_resolve_link(filename);
    case 'data'
        isRaw = (length(filename) > 9) && ~isempty(strfind(filename, 'data_0raw'));
        if isRaw
            FileMat = in_bst_data(filename, 'F');
            filename = FileMat.F.filename;
        end
end
% Folder / file
if isdir(filename)
    filepath = filename;
else
    filepath = bst_fileparts(filename);
end

% === WINDOWS ===
if strncmp(computer,'PC',2)
    % Open Windows CMD
    if isTerminal
        dos(['cmd.exe /k cd ' filepath]);
    % Open Windows Explorer
    else
        dos(['explorer /e,/select,' filename]);
    end

% === MAC OS ===
elseif strncmp(computer,'MAC',3)
    % Open Terminal
    if isTerminal
        system(['open -a Terminal.app']);
    % Open Finder
    else
        system(['open -a Finder.app "' filepath '"']);
    end

% === OTHER SYSTEMS ===
else
    % === TERMINAL ===
    if isTerminal
        xdg_deskptop = getenv('XDG_CURRENT_DESKTOP');
        % KDE default terminal
        if strcmp('KDE', xdg_deskptop) && (system('which konsole > /dev/null') == 0)
            [status, cmdout] = system(['export LD_LIBRARY_PATH=/usr/lib/x86_64-linux-gnu && konsole --workdir "' filepath '" &']);
            if status == 0, return; end
        end
        % MATE default terminal
        if strcmp('MATE', xdg_deskptop) && (system('which mate-terminal > /dev/null') == 0)
            [status, cmdout] = system(['mate-terminal --working-directory="' filepath '" &']);
            if status == 0, return; end
        end
        % GNOME terminal
        if system('which gnome-terminal > /dev/null') == 0
            [status, cmdout] = system(['gnome-terminal --working-directory="' filepath '" &']);
            if status == 0, return; end
        end
        % XTERM
        if system('which xterm > /dev/null')
            [status, cmdout] = system(['xterm -e ''cd "' filepath '" && /bin/bash'' &']);
            if status == 0, return; end
        end
        % Error
        disp(cmdout);
        error('No terminal emulator found for your operating system.');

    % === FILE EXPLORERS ===
    else
        % Any X Desktop Group (XDG) compliant
        ixXdg = system('which xdg-open > /dev/null');
        if (ixXdg == 0)
            [status, cmdout] = system(['xdg-open "' filepath '"']);
            if status == 0, return; end
        end
        % DOLPHIN (KDE)
        isKde = system('which dolphin > /dev/null');
        if (isKde == 0)
            [status, cmdout] = system(['export LD_LIBRARY_PATH=/usr/lib/x86_64-linux-gnu && dolphin "' filepath '" &']);
            if status == 0, return; end
        end
        % NAUTILUS (GNOME)
        isGnome = system('which nautilus > /dev/null');
        if (isGnome == 0)
            [status, cmdout] = system(['nautilus "' filepath '" &']);
            if status == 0, return; end
        end
        % Error
        disp(cmdout);
        error('No file manager found for your operating system.');
    end
end

author: brainstorm-tools


text:
function g = sigmoid(z)
%SIGMOID Compute sigmoid functoon
%   J = SIGMOID(z) computes the sigmoid of z.

% You need to return the following variables correctly
g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the sigmoid of each value of z (z can be a matrix,
%               vector or scalar).

g = 1 ./ (1 + e .^ (-z));

% =============================================================

end

author: gopaczewski

repo: coursera-ml


text:
function sequentialNames(varargin)
%SEQUENTIALNAMES Sets names of all unnamed layers based on type and order
% OBJ.SEQUENTIALNAMES() replaces empty layer names with new names
% sequentially, based on their function handle and execution order. OBJ
% is the network output layer.
% Example layer names are 'conv1', 'conv2', 'pool1'.
%
% OBJ1.SEQUENTIALNAMES(OBJ2, ...) or SEQUENTIALNAMES(OBJ1, OBJ2, ...) is
% the same, but considering a network with several outputs OBJ1, OBJ2, ...
%
% SEQUENTIALNAMES(..., MODIFIER) also specifies a function handle to be
% evaluated on each name, possibly modifying it (e.g. to append a prefix or
% suffix).
%
% See also Layer.workspaceNames.

% Copyright (C) 2016 Joao F. Henriques.
% 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).

  if ~isempty(varargin) && isa(varargin{end}, 'function_handle')
    modifier = varargin{end} ;
    varargin(end) = [] ;
  else
    modifier = @deal ;
  end

  % list all layer objects, in forward order
  assert(~isempty(varargin), 'Not enough input arguments.') ;
  objs = cell(1, numel(varargin)) ;
  for i = 1:numel(varargin)
    objs{i} = varargin{i}.find() ;
  end
  objs = [objs{:}] ;  % flatten array of results

  % automatically set missing names, according to the execution order
  for k = 1:numel(objs)
    if isa(objs{k}, 'Input')  % input1, input2, ...
      if isempty(objs{k}.name)
        n = nnz(cellfun(@(o) isa(o, 'Input'), objs(1:k))) ;
        objs{k}.name = modifier(sprintf('input%i', n)) ;  %#ok<*AGROW>
      end
    elseif ~isempty(objs{k}.func)  % conv1, conv2...
      name = objs{k}.name ;
      if isempty(name)
        n = nnz(cellfun(@(o) isequal(o.func, objs{k}.func), objs(1:k))) ;
        name = func2str(objs{k}.func) ;  %#ok<*PROP>

        if strncmp(name, 'vl_nn', 5), name(1:5) = [] ; end  % shorten names
        if strcmp(name, 'bnorm_wrapper'), name = 'bnorm' ; end
        name = sprintf('%s%i', name, n) ;
        objs{k}.name = modifier(name) ;
      end

      % set dependant Param names for vl_nnconv and vl_nnbnorm
      suffix = [] ;
      if isequal(objs{k}.func, @vl_nnconv)
        suffix = {[], 'filters', 'biases'} ;
      elseif isequal(objs{k}.func, @vl_nnbnorm_wrapper)
        suffix = {[], 'g', 'b', 'moments'} ;
      end
      inputs = objs{k}.inputs ;
      if ~isempty(suffix)
        for i = 2 : min(numel(suffix), numel(inputs))
          if isa(inputs{i}, 'Param') && isempty(inputs{i}.name)
            inputs{i}.name = modifier(sprintf('%s_%s', name, suffix{i})) ;
          end
        end
      end

      % set dependant Param names for general layers: name_p1, name_p2, ...
      ps = find(cellfun(@(o) isa(o, 'Param'), inputs)) ;
      for i = 1:numel(ps)
        if isempty(inputs{ps(i)}.name)
          inputs{ps(i)}.name = modifier(sprintf('%s_p%i', name, i)) ;
        end
      end
    end
  end

  % name any remaining objects by class
  for k = 1:numel(objs)
    if isempty(objs{k}.name)
      n = nnz(cellfun(@(o) isa(o, class(objs{k})), objs(1:k))) ;
      objs{k}.name = modifier(sprintf('%s%i', lower(class(objs{k})), n)) ;
    end
  end
end

repo: MFT


text:
function [SegLabel,NcutDiscrete,NcutEigenvectors,NcutEigenvalues,W,imageEdges]= NcutImage(I,nbSegments);
%  [SegLabel,NcutDiscrete,NcutEigenvectors,NcutEigenvalues,W,imageEdges]= NcutImage(I);
%  Input: I = brightness image
%         nbSegments = number of segmentation desired
%  Output: SegLable = label map of the segmented image
%          NcutDiscrete = Discretized Ncut vectors
%
% Timothee Cour, Stella Yu, Jianbo Shi, 2004.

if nargin <2,
   nbSegments = 10;
end

[W,imageEdges] = ICgraph(I);

[NcutDiscrete,NcutEigenvectors,NcutEigenvalues] = ncutW(W,nbSegments);

%% generate segmentation label map
[nr,nc,nb] = size(I);

SegLabel = zeros(nr,nc);
for j=1:size(NcutDiscrete,2),
    SegLabel = SegLabel + j*reshape(NcutDiscrete(:,j),nr,nc);
end

author: jwyang


text:
%  INTERNAL FUNCTION: creates impact matrix for contemporaneous and future shocks
%
%  ::
%
%    M=stochastic_impact_cumulator(model,y0,nsteps)
%    M=stochastic_impact_cumulator(model,y0,nsteps,y_pos)
%    M=stochastic_impact_cumulator(model,y0,nsteps,y_pos,e_pos)
%    M=stochastic_impact_cumulator(model,y0,nsteps,y_pos,e_pos,states)
%    [M,ufkst,states,PAI,TT,Q]=stochastic_impact_cumulator(...)
%
%  Args:
%
%     - **model** [struct]:
%
%       - **T** [1 x h cell]: solution of the model,
%         T{regime}=[Ty,Tsig,Te_0,...Te_k]
%       - **sstate** [1 x h cell]: steady state in each regime
%       - **state_cols** [vector|{1:ny}]: location of state endogenous
%         variables in y0 (see below).
%       - **k** [scalar]: anticipation horizon (Beyond the current period)
%       - **Qfunc** [empty|function_handle]: endogenous transition matrix
%
%     - **y0** [ny x 1 vector]: initial conditions
%     - **nsteps** [integer]: number of forecast steps
%     - **y_pos** [vector]: location of restricted endogenous variables.
%     - **e_pos** [vector]: location of restricted shocks
%     - **states** [empty|nsteps x 1 vector]: states visited for each forecast
%       step.
%
%  Returns:
%     :
%
%     - **M** [struct]:
%
%       - **R** [matrix]: convoluted restrictions of shocks stemming from the
%         restrictions on endogenous variables
%       - **ufkst** [matrix]: unconditional forecasts for the restricted
%         endogenous variables (Excluding the initial conditions!!!).
%       - **const** [vector]: impact of the constant (steady state + risk) for
%         the restricted endogenous variables.
%       - **S** [matrix]: direct restrictions on shocks
%       - **nshocks** [integer]: number of shocks
%       - **ny** [integer]: number of endogenous variables
%
%     - **ufkst** [ny x (nsteps+1) matrix]: Unconditional forecasts mean (with
%       initial condition at the beginning!!!)
%     - **states** [nsteps x 1 vector]: regimes visited for each forecast step.
%     - **PAI** [h x nsteps matrix]: choice probabilities of regimes for each
%       step
%     - **TT** [matrix]: convolution of autoregressive terms for the
%       restrictions
%     - **Q** [h x h x nsteps array]: Time series of transition matrices
%
%

author: jmaih


text:
function G = filterGauss( dims, mu, C, show )
% n-dimensional Gaussian filter.
%
% Creates an image of a Gaussian with arbitrary covariance matrix. The
% dimensionality and size of the filter is determined by dims (eg dims=[10
% 10] creates a 2D filter of size 10x10). If mu==[], it is calculated to be
% the center of the n-dim image.  C can be a full nxn covariance matrix, or
% an nx1 vector of variance.  In the latter case C is calculated as
% C=diag(C). If C=[]; then C=(dims/6).^2, ie it is transformed into a
% vector of variances such that along each dimension the variance is equal
% to (siz/6)^2.
%
% USAGE
%  G = filterGauss( dims, [mu], [C], [show] )
%
% INPUTS
%  dims    - n element vector of dimensions of final Gaussian
%  mu      - [] n element vector specifying the mean
%  C       - [] nxn cov matrix, nx1 set of vars, or variance
%  show    - [0] figure to use for optional display
%
% OUTPUTS
%  G       - image of the created Gaussian
%
% EXAMPLE
%  g = filterGauss( 21, [], 4, 1); %1D
%  sig=3; G = filterGauss( 4*[sig sig] + 1, [], [sig sig].^2, 2 ); %2D
%  R = rotationMatrix( [1,1,0], pi/4 );
%  C = R'*[10^2 0 0; 0 5^2 0; 0 0 16^2]*R;
%  G3 = filterGauss( [51,51,51], [], C, 3 ); %3D
%
% See also NORMPDF2
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

nd = length( dims );
if( nargin<2 || isempty(mu)); mu=(dims+1)/2; end
if( nargin<3 || isempty(C)); C=(dims/6).^2; end
if( nargin<4 || isempty(show) ); show=0; end

if( length(mu)~=nd ); error('invalid mu'); end

if( nd==1 ) % fast special case
  xs = 1:dims(1);
  G = exp(-(xs-mu).*(xs-mu)/(2*C))';

else
  % make C have correct dimensions
  if( numel(C)==1 ); C=repmat(C,[1 nd]); end
  if( size(C,1)==1 || size(C,2)==1 ); C=diag(C); end
  if( any(size(C)~=nd)); error( 'invalid C'); end

  % get vector of grid locations
  temp = cell(1,nd);
  for d=1:nd; temp{d} = 1:dims(d); end
  [ temp{:}] = ndgrid( temp{:} );
  xs = zeros( nd, prod(dims) );
  for d=1:nd; xs( d, : ) = temp{d}(:)'; end

  % evaluate the Gaussian at those points
  G = normpdf2( xs, mu, C );
  if( nd>1 ); G = reshape( G, dims ); end
end

% suppress very small values and normalize
G(G<eps*max(G(:))*10) = 0;
G = G/sum(G(:));

% display
if( show );
  if( nd==3 )
    filterVisualize( G, show, .2 );
  else
    filterVisualize( G, show );
  end
end

author: zouchuhang

repo: LayoutNet


text:
function a = any(f, dim)
%ANY   True if any element of a SINGFUN is a nonzero number. ANY ignores
%      entries that are NaN (Not a Number).
%   ANY(F, DIM), where DIM is the dimension along which ANY is taken. Note
%   that array-valued SINGFUN is not supported. If DIM is 1, then ANY
%   returns a logical value which is TRUE if any element of F is nonzero.
%   If DIM is 2, ANY returns a SMOOTHFUN which takes the value 1 if F is
%   nonzero.  In this case, F must either be identically zero or have no
%   roots in its domain.  Otherwise, garbage is returned without warning.
%
%   ANY(F) is shorthand for ANY(F, 1).

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information

% Parse inputs:
if ( nargin < 2 )
    dim = 1;
end

a = any(f.smoothPart, dim);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function status = fseek( fid, offset, origin )
%FSEEK Set file position indicator.
% STATUS = FSEEK(FID, OFFSET, ORIGIN) mimics the behavior of MATLAB's FSEEK
% function, but allows for seeking past the end of a file (just appends
% zeros to the end of the file until it reaches desired location).
%
%   OFFSET values are interpreted as follows:
%       >= 0   Move position indicator OFFSET bytes after ORIGIN.
%       < 0    Move position indicator OFFSET bytes before ORIGIN.
%
%   ORIGIN values are interpreted as follows:
%       'bof' or -1   Beginning of file
%       'cof' or  0   Current position in file
%       'eof' or  1   End of file
%
%   STATUS is 0 on success and -1 on failure.  If an error occurs, use
%   FERROR to get more information.
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

status = fseek(fid, offset, origin); % Try MATLAB's fseek first
if status < 0 % Didn't work
    [message, errnum] = ferror(fid);
    if errnum == -27 % "Offset is bad - after end-of-file or last character written."
        % This is not an error in our case.  We will do special handling to
        % allow this fseek to skip to position past current end of file.
        switch origin
            case {-1, 'bof'}
                offset = offset - ftell(fid); % Convert offset to be from current position, not beginning of file
            case {0, 'cof'} % No change required
            case {1, 'eof'}
                fseek(fid,0,'eof');
        end
        if offset>0
            % Use fwrite skip parameter to append to end of file an
            % arbitrary amount
            status = -(fwrite(fid, uint8(0), 'uint8', offset-1) ~= 1);
        end
    end
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function [ params ] = tapas_BLT_runIntroTrain( params )
% tapas_BLT_runIntroTrain Presents written instructions and logs the event
%
% About:
%   Coded by F.Petzschner 19. April 2017
%	Last change: F.Petzschner 19. April 2017

for i = 1:params.nIntroTrain

    Screen('DrawTexture', params.screen.window, params.text.IntroTrain{i}, [], params.screen.fit, 0);
    Screen('Flip', params.screen.window);
    wait2_escapeOption(params.dur.showIntroScreen, params);
    detecting = 1;

    if i == params.nIntroTrain
        wait2(3000);
        break
    else
        while detecting
            keyCode = 0;
            keyCode = detectkey(params);
            if keyCode == params.keys.two
                    Screen('DrawTexture', params.screen.window, params.text.imgbreak, [], params.screen.fit, 0);
                    Screen('Flip', params.screen.window);
                    wait2(1000);
                    keyCode2 = detectkey(params);
                    if keyCode2 == params.keys.one
                       detecting = 0;
                    end
            elseif keyCode == params.keys.one
                detecting = 0;
            end
        end
    end

    fprintf(['Intro Training Slide...',num2str(i)])
end

author: translationalneuromodeling

repo: tapas


text:
  function [phantom, params] = rect_im(ig, params, varargin)
%|function [phantom, params] = rect_im(ig, params, options)
%|
%| generate rectangle phantom image from parameters:
%|	[x_center y_center x_width y_width angle_degrees amplitude]
%|
%| in
%|	ig			image_geom() object
%|	params			rect parameters, if empty, use the default
%|				params = rect_im_default_parameters(xfov, yfov)
%|
%| options
%|	'oversample'	int	oversampling factor, for grayscale boundaries
%|	'hu_scale'	float	use 1000 to scale shepp-logan to HU (default: 1)
%|	'fov'		float	default ig.fov
%|
%| out
%|	phantom		[nx ny]	image
%|
%| 2008-08-04, Yong Long, adapted from ellipse_im()
%| Copyright 2006-2-2, Jeff Fessler, University of Michigan

if nargin == 1 && streq(ig, 'test'), rect_im_test, return, end
if nargin < 1, ir_usage, end

arg.oversample = 1;
arg.replace = 0;
arg.hu_scale = 1;
arg.fov = [];
arg.chat = false;
arg = vararg_pair(arg, varargin);
if isempty(arg.fov), arg.fov = ig.fov; end

if ~isvar('params') || isempty(params)
	params = rect_im_default_parameters(arg.fov, arg.fov);
end

params(:,6) = params(:,6) * arg.hu_scale;

do_fast = params(:,5) == 0; % non-rotated ones are faster

if arg.chat
	printm('%d of %d fast', sum(do_fast), numel(do_fast))
end

if any(do_fast)
	phantom = rect_im_fast(params(do_fast,:), ...
		ig.nx, ig.ny, ig.dx, ig.dy, ig.offset_x, ig.offset_y, ...
		arg.replace);
else
	phantom = 0;
end

if any(~do_fast)
	phantom = phantom + rect_im_slow(params(~do_fast,:), ...
		ig.nx, ig.ny, ig.dx, ig.dy, ig.offset_x, ig.offset_y, ...
		arg.oversample, arg.replace);
end

% rect_im_fast()
% fast version for non-rotated rectangles
% using exact integration over each pixel so over-sampling is irrelevant
function phantom = rect_im_fast(params, ...
	nx, ny, dx, dy, offset_x, offset_y, replace)

if size(params,2) ~= 6
	error 'bad rect parameter vector size'
end

phantom = zeros(nx, ny);

wx = (nx-1)/2 + offset_x;
wy = (ny-1)/2 + offset_y;
x1 = ([0:nx-1]' - wx) * dx; % col
y1 = ([0:ny-1] - wy) * dy; % row

fun = @(x1, x2, wx) ... % integrated rect(x/wx) function from x1 to x2
	max(min(x2, wx/2) - max(x1, -wx/2), 0);

%ticker reset
ne = nrow(params);
for ie = 1:ne
%	ticker(mfilename, ie, ne)

	rect = params(ie, :);
	cx = rect(1);	wx = rect(3);
	cy = rect(2);	wy = rect(4);
	theta = deg2rad(rect(5));
	if theta ~= 0, fail 'theta=0 required', end
	x = x1 - cx;
	y = y1 - cy;
	tx = fun(x-abs(dx)/2, x+abs(dx)/2, wx) / abs(dx);
	ty = fun(y-abs(dy)/2, y+abs(dy)/2, wy) / abs(dy);
	tmp = single(tx) * single(ty); % outer product (separable)
	if replace
		phantom(tmp > 0) = rect(6);
	else
		phantom = phantom + rect(6) * tmp;
	end
end

% rect_im_slow()
% slower version that handles rotation too and uses over-sampling
function phantom = rect_im_slow(params, ...
	nx, ny, dx, dy, offset_x, offset_y, over, replace)

if size(params,2) ~= 6
	error 'bad rect parameter vector size'
end

phantom = zeros(nx*over, ny*over);

wx = (nx*over-1)/2 + offset_x*over;
wy = (ny*over-1)/2 + offset_y*over;
xx = ([0:nx*over-1] - wx) * dx / over;
yy = ([0:ny*over-1] - wy) * dy / over;
[xx yy] = ndgrid(xx, yy);

ticker reset
ne = nrow(params);
for ie = 1:ne
	ticker(mfilename, ie, ne)

	rect = params(ie, :);
	cx = rect(1);	wx = rect(3);
	cy = rect(2);	wy = rect(4);
	theta = deg2rad(rect(5));
	x = cos(theta) * (xx-cx) + sin(theta) * (yy-cy);
	y = -sin(theta) * (xx-cx) + cos(theta) * (yy-cy);
	tmp = abs(x / wx) < 1/2 & abs(y / wy) < 1/2;
	if replace
		phantom(tmp > 0) = rect(6);
	else
		phantom = phantom + rect(6) * tmp;
	end
end

phantom = downsample2(phantom, over);

%
% default params for rectangles
% the first four columns are unitless for fov=64
%
function params = rect_im_default_parameters(xfov, yfov)
f = 1/64;
params = [ ...
	0	0	50	50	0	1
	10	-16	25	16	0	-0.5
	-13	15	13	13	1*45	1
...
	-18	0	1	1	0	1
	-12	0	1	1	0	1
	-6	0	1	1	0	1
	0	0	1	1	0	1
	6	0	1	1	0	1
	12	0	1	1	0	1
	18	0	1	1	0	1
];

params(:,[1 3]) = params(:,[1 3]) * xfov / 64;
params(:,[2 4]) = params(:,[2 4]) * yfov / 64;

% rect_im_test()
function rect_im_test
ig = image_geom('nx', 2^8, 'ny', 2^8, 'fov', 100);
pic = rect_im(ig, [[0.5 0 3 20]*ig.dx 0 1], 'oversample', 3);
jf_equal(unique(pic), [0;1])
im(pic)
pic = rect_im(ig, [], 'oversample', 3, 'chat', im);
%unique(pic)
im(ig.x, ig.y, pic, 'default rects'), cbar

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function pass = test_mixed_tech( pref )
% This tests the chebfun2 constructor for mixed technology.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e2 * pref.cheb2Prefs.chebfun2eps;

pass = [];

%% 'trig' in one or both dimensions

% Construct from function handle
u = @(x,y) sin(2*pi*x).*cos(2*pi*y);

f1 = chebfun2(u, 'trigx');
f2 = chebfun2(u, 'periodicy');
f3 = chebfun2(u, 'trig');
f4 = chebfun2(u);

pass(end+1) = ( norm(f1 - f2) < tol );
pass(end+1) = ( norm(f2 - f3) < tol );
pass(end+1) = ( norm(f3 - f4) < tol );

% Construct from data
m = 31;
n = 32;
u = @(x,y) sin(2*pi*x).*cos(2*pi*y);
[xc, yc] = chebpts2(n, m);
[xt, yt] = meshgrid(trigpts(n), trigpts(m));

f1 = chebfun2(u(xt,yc), 'trigx');
f2 = chebfun2(u(xc,yt), 'periodicy');
f3 = chebfun2(u(xt,yt), 'trig');
f4 = chebfun2(u(xc,yc));

pass(end+1) = ( norm(f1 - f2) < tol );
pass(end+1) = ( norm(f2 - f3) < tol );
pass(end+1) = ( norm(f3 - f4) < tol );

%% 'equi' in one or both dimensions
m = 31;
n = 32;
u = @(x,y) sin(x.*y);
[xc, yc] = chebpts2(n, m);
[xe, ye] = meshgrid(linspace(-1,1,n), linspace(-1,1,m));

f1 = chebfun2(u(xe,yc), 'equix');
f2 = chebfun2(u(xc,ye), 'equiy');
f3 = chebfun2(u(xe,ye), 'equi');
f4 = chebfun2(u(xc,yc));

pass(end+1) = ( norm(f1 - f2) < tol );
pass(end+1) = ( norm(f2 - f3) < tol );
pass(end+1) = ( norm(f3 - f4) < tol );

%% 'coeffs' in one or both dimensions
m = 31;
n = 32;
u = @(x,y) sin(x.*y);
[xc, yc] = chebpts2(n, m);

vals_vals     = u(xc,yc);
vals_coeffs   = chebtech2.vals2coeffs( vals_vals );
coeffs_vals   = chebtech2.vals2coeffs( vals_vals.' ).';
coeffs_coeffs = chebtech2.vals2coeffs( vals_coeffs.' ).';

f1 = chebfun2(u);
f2 = chebfun2(coeffs_vals, 'coeffsx');
f3 = chebfun2(vals_coeffs, 'coeffsy');
f4 = chebfun2(coeffs_coeffs, 'coeffs');

pass(end+1) = ( norm(f1 - f2) < tol );
pass(end+1) = ( norm(f2 - f3) < tol );
pass(end+1) = ( norm(f3 - f4) < tol );

%% 'trig' and 'coeffs' in one or both dimensions
m = 31;
n = 32;
u = @(x,y) sin(2*pi*x).*cos(2*pi*y);
[xc, yc] = chebpts2(n, m);
[xt, yt] = meshgrid(trigpts(n), trigpts(m));

trigvals_chebvals     = u(xt,yc);
chebvals_trigvals     = u(xc,yt);
trigvals_trigvals     = u(xt,yt);
trigvals_chebcoeffs   = chebtech2.vals2coeffs( trigvals_chebvals     );
trigcoeffs_trigvals   =  trigtech.vals2coeffs( trigvals_trigvals.'   ).';
trigcoeffs_chebvals   =  trigtech.vals2coeffs( trigvals_chebvals.'   ).';
trigcoeffs_chebcoeffs =  trigtech.vals2coeffs( trigvals_chebcoeffs.' ).';
chebvals_trigcoeffs   =  trigtech.vals2coeffs( chebvals_trigvals     );

f1 = chebfun2(u);
f2 = chebfun2(trigcoeffs_chebvals, 'coeffsx', 'trigx');
f3 = chebfun2(trigvals_chebcoeffs, 'coeffsy', 'trigx');
f4 = chebfun2(trigcoeffs_chebcoeffs, 'coeffs', 'trigx');
f5 = chebfun2(trigcoeffs_trigvals, 'trig', 'coeffsx');
f6 = chebfun2(chebvals_trigcoeffs, 'trigy', 'coeffsy');

pass(end+1) = ( norm(f1 - f2) < tol );
pass(end+1) = ( norm(f2 - f3) < tol );
pass(end+1) = ( norm(f3 - f4) < tol );
pass(end+1) = ( norm(f4 - f5) < tol );
pass(end+1) = ( norm(f5 - f6) < tol );

%% Preferences in one or both dimensions
p = pref;
p.tech = @trigtech;

f = chebfun2(@(x,y) sin(2*pi*x).*y, {p, []});
g = chebfun2(@(x,y) sin(2*pi*x).*y);
pass(end+1) = ( norm(f - g) < tol );

f = chebfun2(@(x,y) sin(2*pi*y).*x, {[], p});
g = chebfun2(@(x,y) sin(2*pi*y).*x);
pass(end+1) = ( norm(f - g) < tol );

f = chebfun2(@(x,y) sin(2*pi*x).*sin(2*pi*y), {p, p});
g = chebfun2(@(x,y) sin(2*pi*x).*sin(2*pi*y));
pass(end+1) = ( norm(f - g) < tol );

%% The last argument takes precedence
f = chebfun2(@(x,y) sin(2*pi*x).*y, 'trig', 'trigx');
g = chebfun2(@(x,y) sin(2*pi*x).*y, 'trigx');
pass(end+1) = ( norm(f - g) < tol );

m = 31;
n = 32;
u = @(x,y) sin(x.*y);
[xx, yy] = meshgrid(linspace(-1,1,n), linspace(-1,1,m));
uu = u(xx,yy);
f = chebfun2(uu, 'equix', 'equi');
g = chebfun2(uu, 'equi');
pass(end+1) = ( norm(f - g) < tol );

prefx = pref;
prefx.tech = @trigtech;
f = chebfun2(@(x,y) x, {prefx, []}, pref);
g = chebfun2(@(x,y) x, pref);
pass(end+1) = ( norm(f - g) < tol );

end


text:
function v_leg = chebcoeffs2legvals(c_cheb)
%CHEBCOEFFS2LEGVALS  Convert Chebyshev coefficients to Legendre values.
%   V_LEG = CHEBCOEFFS2LEGVALS(C_CHEB) converts the vector C_CHEB of Chebyshev
%   coefficients to a vector V_LEG of Legendre values such that
%       C_CHEB(1)*T0(X_k) + ... + C_CHEB(N)*T{N-1} = V_LEG_k, k = 0, ... N-1,
%   where X_k are the Legendre nodes returned by LEGPTS().
%
% See also LEGPTS.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% This command is a wrapper for chebfun.ndct.
v_leg = chebfun.ndct(c_cheb);

end


text:
function exportIVH(structNum,scanSet,Opt)
% exportDVH
% This function exports DVH in an EXCEL format. Opt give the option to
% choose if you want Normalised or Absolute DVH
% written DK
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.
%
% Usage
% exportDVH(structNum,Opt)

% for command line help document
if ~exist('structNum')& ~exist('scanSet')& ~exist('Opt')
    prompt = {'Enter the structure Number';'Enter the scan number'; 'Enter "abs" for absolute OR "nor" for normalized DVH'};
            dlg_title = 'Export IVH';
            num_lines = 1;
            def = {'';'';''};
            outPutQst = inputdlg(prompt,dlg_title,num_lines,def);
            if isempty(outPutQst{1}) | isempty(outPutQst{2})| isempty(outPutQst{3})
                warning('Need to enter all the inputs');
                return
            else
                structNum = str2num(outPutQst{1});
                scanSet = str2num(outPutQst{2});
                Opt = outPutQst{3};
            end
end

path = uigetdir( cd,'Select destination Directory for DVH export');
global planC
indexS = planC{end};
structureCell = planC{indexS.structures};
optS = opts4Exe([getCERRPath,'CERROptions.json']);
%loop over all the structures that need to be exported
for i = 1:length(structNum)
    name = structureCell(structNum(i)).structureName;
    [scansV, volsV] = getIVH(structNum(i), scanSet(i), planC);
    [scanBinsV, volsHistV] = doseHist(scansV, volsV, optS.DVHBinWidth);
    cumVolsV = cumsum(volsHistV);
    cumVols2V  = cumVolsV(end) - cumVolsV;  %cumVolsV is the cumulative volume lt that corresponding scan
    switch upper(Opt)
        case 'ABS'
%             if abs flag is set just export the values as it is
            fVol = cumVols2V;
        case 'NOR'
            %Normalizing the volume
            fVol = cumVols2V/cumVolsV(end);
    end
    M = [scanBinsV; fVol];
    %Export only NumPts points (for Dr. Myerson)
    NumPts = 200;
    if size(M,2)>NumPts
        indAll = round(linspace(1,size(M,2),NumPts));
        M = M(:,indAll);
    end
    xlswrite(fullfile(path,name), M');
end
clear fVol scansV volsV cumVolsV cumVols2V scanBinsV volsHistV optS

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [path, totalCost, farthestPreviousHop, farthestNextHop] = dijkstra(n, netCostMatrix, s, d, farthestPreviousHop, farthestNextHop)
% path: the list of nodes in the path from source to destination;
% totalCost: the total cost of the path;
% farthestNode: the farthest node to reach for each node after performing
% the routing;
% n: the number of nodes in the network;
% s: source node index;
% d: destination node index;

% clear;
% noOfNodes  = 50;
% rand('state', 0);
% figure(1);
% clf;
% hold on;
% L = 1000;
% R = 200; % maximum range;
% netXloc = rand(1,noOfNodes)*L;
% netYloc = rand(1,noOfNodes)*L;
% for i = 1:noOfNodes
%     plot(netXloc(i), netYloc(i), '.');
%     text(netXloc(i), netYloc(i), num2str(i));
%     for j = 1:noOfNodes
%         distance = sqrt((netXloc(i) - netXloc(j))^2 + (netYloc(i) - netYloc(j))^2);
%         if distance <= R
%             matrix(i, j) = 1;   % there is a link;
%             line([netXloc(i) netXloc(j)], [netYloc(i) netYloc(j)], 'LineStyle', ':');
%         else
%             matrix(i, j) = inf;
%         end;
%     end;
% end;
%
%
% activeNodes = [];
% for i = 1:noOfNodes,
%     % initialize the farthest node to be itself;
%     farthestPreviousHop(i) = i;     % used to compute the RTS/CTS range;
%     farthestNextHop(i) = i;
% end;
%
% [path, totalCost, farthestPreviousHop, farthestNextHop] = dijkstra(noOfNodes, matrix, 1, 15, farthestPreviousHop, farthestNextHop);
% path
% totalCost
% if length(path) ~= 0
%     for i = 1:(length(path)-1)
%         line([netXloc(path(i)) netXloc(path(i+1))], [netYloc(path(i)) netYloc(path(i+1))], 'Color','r','LineWidth', 0.50, 'LineStyle', '-.');
%     end;
% end;
% hold off;
% return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% all the nodes are un-visited;
visited(1:n) = 0;

distance(1:n) = inf;    % it stores the shortest distance between each node and the source node;
parent(1:n) = 0;

distance(s) = 0;
for i = 1:(n-1),
    temp = [];
    for h = 1:n,
         if visited(h) == 0   % in the tree;
             temp=[temp distance(h)];
         else
             temp=[temp inf];
         end
     end;
     [t, u] = min(temp);    % it starts from node with the shortest distance to the source;
     visited(u) = 1;       % mark it as visited;
     for v = 1:n,           % for each neighbors of node u;
         if ( ( netCostMatrix(u, v) + distance(u)) < distance(v) )
             distance(v) = distance(u) + netCostMatrix(u, v);   % update the shortest distance when a shorter path is found;
             parent(v) = u;                                     % update its parent;
         end;
     end;
end;

path = [];
if parent(d) ~= 0   % if there is a path!
    t = d;
    path = [d];
    while t ~= s
        p = parent(t);
        path = [p path];

        if netCostMatrix(t, farthestPreviousHop(t)) < netCostMatrix(t, p)
            farthestPreviousHop(t) = p;
        end;
        if netCostMatrix(p, farthestNextHop(p)) < netCostMatrix(p, t)
            farthestNextHop(p) = t;
        end;

        t = p;
    end;
end;

totalCost = distance(d);

return;

author: Sable

repo: mcbench-benchmarks


text:
% Check if an input file is in BNT format
%
% Checks format of the input file. Returns TRUE if the file is in BNT
% format or FALSE otherwise. If file doesn't exist returns FALSE.
%
function result = isBntFormat(filename)
    result = false;
    if exist(filename, 'file') == 0
        return;
    end

    fid = data.safefopen(filename, 'r');
    str = fgets(fid);

    result = ischar(str) && ~isempty(strfind(lower(str), 'name')) && ~isempty(strfind(lower(str), 'version'));
end

author: kavli-ntnu


text:
function weight = olmar1_expert(data, weight_o, epsilon, W)
% This program generates portfolio for a specified parameter setting.
% OLMAR-1 expert
%
% function weight = olmar1_expert(data, weight_o, epsilon, W)
%
% weight: experts portfolio, used for next rebalance/combination
%
% data: market sequence vectors
% weight_o: last portfolio
% epsilon: mean reversion threshold
% W: window size for calculating moving average
%
% Example: weight = olmar1_expert(data, weight_o, epsilon, W)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[T, N] = size(data);

if (T < W+1)
    data_phi = data(T, :);
else
    data_phi = zeros(1, N);
    tmp_x = ones(1, N);
    for i = 1:W,
        data_phi = data_phi + 1./tmp_x;
        tmp_x = tmp_x.*data(T-i+1, :);
    end

    data_phi = data_phi*(1/W);
end

% Step 3: Suffer loss
ell = max([0 epsilon - data_phi*weight_o]);

% Step 4: Set parameter
x_bar = mean(data_phi);
denominator = (data_phi - x_bar)*(data_phi - x_bar)';
if (~eq(denominator, 0.0)),
    lambda = ell / denominator;
else  % Zero volatility
    lambda = 0;
end

% Step 5: Update portfolio
weight = weight_o + lambda*(data_phi' - x_bar);

% Step 6: Normalize portfolio
weight = simplex_projection(weight, 1);

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function trans = localToGlobal3d(varargin)
%LOCALTOGLOBAL3D Transformation matrix from local to global coordinate system
%
%   TRANS = localToGlobal3d(CENTER, THETA, PHI, PSI)
%   Compute the transformation matrix from a local (or modelling)
%   coordinate system to the global (or world) coordinate system.
%   This is a low-level function, used by several drawing functions.
%
%   The transform is defined by:
%   - CENTER: the position of the local origin into the World coordinate
%       system
%   - THETA: colatitude, defined as the angle with the Oz axis (between 0
%       and 180 degrees), positive in the direction of the of Oy axis.
%   - PHI: azimut, defined as the angle of the normal with the Ox axis,
%       between 0 and 360 degrees
%   - PSI: intrinsic rotation, corresponding to the rotation of the object
%       around the direction vector, between 0 and 360 degrees
%
%   The resulting transform is obtained by applying (in that order):
%   - Rotation by PSI   around he Z-axis
%   - Rotation by THETA around the Y-axis
%   - Rotation by PHI   around the Z-axis
%   - Translation by vector CENTER
%   This corresponds to Euler ZYZ rotation, using angles PHI, THETA and
%   PSI.
%
%   The 'eulerAnglesToRotation3d' function may better suit your needs as
%   it is more 'natural'.
%
%   Example
%   localToGlobal3d
%
%   See also
%   transforms3d, createEulerAnglesRotation
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-06-19,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

%   HISTORY
%   19/08/2009 fix bug in parsing center with 4 args
%   2011-06-21 use degrees

% extract the components of the transform
if nargin == 1
    % all components are bundled in  the first argument
    var     = varargin{1};
    center  = var(1:3);
    theta   = var(4);
    phi     = var(5);
    psi     = 0;
    if length(var) > 5
        psi = var(6);
    end

elseif nargin == 4
    % arguments = center, then the 3 angles
    center  = varargin{1};
    theta   = varargin{2};
    phi     = varargin{3};
    psi     = varargin{4};

elseif nargin > 4
    % center is given in 3 arguments, then 3 angles
    center  = [varargin{1} varargin{2} varargin{3}];
    theta   = varargin{4};
    phi     = varargin{5};
    psi     = 0;
    if nargin > 5
        psi = varargin{6};
    end
end

% conversion from degrees to radians
k = pi / 180;

% rotation around normal vector axis
rot1    = createRotationOz(psi * k);

% colatitude
rot2    = createRotationOy(theta * k);

% longitude
rot3    = createRotationOz(phi * k);

% shift center
tr      = createTranslation3d(center);

% create final transform by concatenating transforms
trans   = tr * rot3 * rot2 * rot1;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function hdr = mrReadEfile(EfileName)
% header = mrReadEfile(EfileName)
%
% Returns header info from Efile produced by G. Glover's auto-recon program.
% Output is structure having field names corresponding to Glover's parameters.
%
% Ress 9/01
% $Author: bob $
% $Date: 2008/02/14 19:50:01 $
% ras, 07/05: imported into mrVista 2.0.
% ras, 10/08: reads alignment points for each slice.
if notDefined('EfileName')
	EfileName = mrvSelectFile('r', {'E*.7' '*.*'}, 'Select an E file ...');
end

hdr = [];

fid = fopen(EfileName);
if fid == -1
	Alert(['Could not open file: ' EfileName])
	return
end

% initialize empty fields for alignment info
% * ras 10/08: after some empirical testing, I believe gw_point1 represents
% the upper left-hand corner of each slice in scanner space, gw_point2
% represents the lower left-hand corner, and gw_point3 represents the
% uppper right-hand corner. In the output header structure, I will store
% each set of points as a 3xnSlices array. The rows are the L->R, P->A, and
% I->S position of each point, respectively; and the columns are different
% slices.
hdr.gw_point1 = [];
hdr.gw_point2 = [];
hdr.gw_point3 = [];

nLine = 0;
while 1
	line = fgetl(fid);
	if ~ischar(line), break, end
	ieq = findstr(line, '=');
	lhs = deblank(line(1:ieq-1));
	rhs = line(ieq+2:end);
	switch lhs
		case 'rev'
			% Revision
			hdr.rev = rhs;
		case 'date of scan'
			% Date
			hdr.date = rhs;
			% Make this more reasonable
			hdr.date=[hdr.date(1:6),hdr.date(8:end)];
		case 'time of scan'
			% Time
			hdr.time = rhs;
		case 'patient name'
			% Name
			hdr.name = rhs;
		case 'psd'
			% PSD
			hdr.psd = rhs;
		case 'coil'
			% Coil
			hdr.coil = rhs;
		case 'examnum/seriesnum'
			vals = explode('/', rhs);
			hdr.examNum = str2num(vals{1});
			hdr.seriesNum = str2num(vals{2});
		case 'series description'
			hdr.seriesDescription = rhs;
		case 'patient id'
			hdr.patientID = rhs;
		case 'exam description'
			hdr.examDescription = rhs;
		case 'series description'
			hdr.seriesDescription = rhs;
		case 'num echoes'
			hdr.nechos = getNumberFromString(rhs);
		case 'slquant'
			% Number of slices
			hdr.slquant = getNumberFromString(rhs);
		case 'num time frames'
			% Number of frames
			hdr.nframes = getNumberFromString(rhs);
		case {'numextra discards' 'nextra discards'}
			% Number of discarded shots
			hdr.nextra = getNumberFromString(rhs);
		case 'nshot'
			% Number of interleaves
			hdr.nshots = getNumberFromString(rhs);
		case 'FOV'
			% Field of view (mm)
			hdr.FOV = getNumberFromString(rhs);
		case 'slice thick'
			% Slice thickness (mm)
			hdr.sliceThickness = getNumberFromString(rhs);
		case 'skip'
			% Slice spacing (skip)
			hdr.skip = getNumberFromString(rhs);
		case 'TR'
			% Echo time, TR (ms)
			hdr.TR = getNumberFromString(rhs);
		case 'TE'
			% Repetition time, TE (ms)
			hdr.TE = getNumberFromString(rhs);
		case 'time/frame'
			% Acquisition time (ms)
			hdr.tAcq = getNumberFromString(rhs);
		case 'equiv matrix size'
			% Equivalent matrix size
			hdr.equivMatSize = getNumberFromString(rhs);
		case 'imgsize'
			% Image size
			hdr.imgsize = getNumberFromString(rhs);
		case 'pixel size'
			% Pixel size (mm)
			hdr.pixel = getNumberFromString(rhs);
		case 'freq'
			% Frequency (MHz)
			hdr.freq = getNumberFromString(rhs)/1.e6;
		case 'flip angle'
			% Flip Angle (deg)
			hdr.flipAngle = str2num(rhs);
		case 'R1, R2, TG (mps)'
			% Gains [R1, R2, TG]
			hdr.R1 = getNumberFromString(rhs, 1);
			hdr.R2 = getNumberFromString(rhs, 2);
			hdr.TG = getNumberFromString(rhs, 3);
		case 'BigEndian'
			hdr.BigEndian = str2num(rhs);
		case 'start_loc'
			hdr.start_loc = [sscanf(rhs, '%f %f %f', 3)]';
		case 'gw_point1'
			hdr.gw_point1(:,end+1) = sscanf(rhs, '%f %f %f', 3);
		case 'gw_point2'
			hdr.gw_point2(:,end+1) = sscanf(rhs, '%f %f %f', 3);
		case 'gw_point3'
			hdr.gw_point3(:,end+1) = sscanf(rhs, '%f %f %f', 3);
		otherwise
			% try to record the remaining fields, if the descriptions are
			% valid field names
			try
				hdr.(lhs) = rhs;
			catch
				% oh well...
% 				fprintf('Missed: %s, %s\n', lhs, rhs);  % debug
			end
	end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function computeModelCoefficients_batchHN_Old(pathExperiments,nExp,imbalance,nBatch,matlabPATH,seed)
% -------------------------------------------------------------------------
% function computeModelCoefficients_batchHN(pathExperiments,nExp,fSetNames,imbalance,matlabPATH)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes the final logistic regression coefficients and
% bootstrap confidence intervals of the final models obtained for all
% outcomes analyzed in the HN study. See ref.[1] for more details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathExperiments: Full path to the directory containing all experiments.
%                     --> Ex: '/myProject/WORKSPACE/CV-BASED_RESULTS'
% 2. nExp: Numerical value specifying the number of experiments to analyze.
%          --> Ex: 10
% 3. imbalance: String specifying the type of imbalance-adjustement strategy
%               employed. Either 'IABR' for imbalance-adjusted bootstrap
%               resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%               logistic regression (see ref.[2]).
%               --> Ex: 'IALR'
% 4. nBatch: Number of parallel batch.
%            --> Ex: 8
% 5. matlabPATH: Full path to the MATLAB excutable on the system.
%                --> 'matlab' if a symbolic link to the matlab executable
%                     was previously created.
% -------------------------------------------------------------------------
% OUTPUTS: Final coefficients, model response and confidence intervals
%          saved in the corresponding ../experiment/outcome/fSetName
%          folder.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
time = 30; % Number of seconds to wait before checking if parallel computations are done

for exp = 1:nExp
    cd(fullfile(pathExperiments,['Experiment',num2str(exp)])), load('training')
    cd('FINAL_MODELS'), pathFinalModels = pwd;
    mkdir('batchLog_Coeff'), cd('batchLog_Coeff'), pathBatch = pwd;
    nameOutcomes = fieldnames(training); nOutcomes = numel(nameOutcomes);
    for o = 1:nOutcomes
        outcomes.(nameOutcomes{o}) = training.(nameOutcomes{o}).outcome;
    end
    setNames = fieldnames(training.(nameOutcomes{1}).text);
    [param] = batchExperiments(setNames,outcomes,nBatch); nBatch = length(param);

    % PRODUCE BATCH COMPUTATIONS
    save('workspace','pathFinalModels','training','param','pathBatch','imbalance','seed'), pause(3)
    for i = 1:nBatch
        nameScript = ['batch',num2str(i),'_script.m'];
        fid = fopen(nameScript,'w');
        fprintf(fid,'tic\n');
        fprintf(fid,'load(''workspace'')\n');
        for j = 1:numel(param{i})
            fprintf(fid,['cd(fullfile(pathFinalModels,param{',num2str(i),'}{',num2str(j),'}{2},param{',num2str(i),'}{',num2str(j),'}{1}))\n']);
            fprintf(fid,'load(''finalModel'')\n');
            fprintf(fid,['fprintf(''COMPUTING THE LOGISTIC REGRESSION COEFFICIENTS OF THE FINAL MODEL OF "',param{i}{j}{2},'" OUTCOME, "',param{i}{j}{1},'" FEATURE SET ... '')']);
            fprintf(fid,'\n');
            fprintf(fid,['[coeff,response,modelCI] = computeModelCoefficients_HN(finalModel.Data,training.',param{i}{j}{2},'.outcome,imbalance,seed);\n']);
            fprintf(fid,['fprintf(''DONE!\\n'')']);
            fprintf(fid,'\n');
            fprintf(fid,'save(''coeff'',''coeff''), save(''response'',''response''), save(''modelCI'',''modelCI'')\n');
        end
        fprintf(fid,'cd(pathBatch)\n');
        fprintf(fid,['system(''touch batch',num2str(i),'_end'');\n']);
        fprintf(fid,'clear all\n');
        fprintf(fid,'toc');
        fclose(fid);
        system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
    end
    waitBatch(pathBatch,time,nBatch)
    delete('workspace.mat')
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [ precision, recall ] = ComputePrecisionRecall( tau, sigma, tp,tr,k,fsc_k )
%COMPUTEPRECISIONRECALL Summary of this function goes here
%   Detailed explanation goes here

if nargin == 2
    tr = 0.7;   % recall threshold
    tp = 0.6;   % precision threshold
    k = 2;      % min number of matches, used in penalizing split & merge
    fsc_k = 0.8;% penalize value of split or merge
end

tot_gt = 0;
recall_accum = 0;

tot_detected = 0;
precision_accum = 0;

num_images = numel(tau);
assert(num_images == numel(sigma));

flag_gt = cell(1, num_images);
flag_det = cell(1, num_images);

for ifile=1:num_images

    [num_gt, num_detected] = size( sigma{ifile} );
    tot_gt = tot_gt + num_gt;
    tot_detected = tot_detected + num_detected;

    % ----- mark unprocessed
    flag_gt{ifile} = zeros(num_gt, 1);
    flag_det{ifile} = zeros(num_detected, 1);

    % ---------------------------------------
    % check one-to-one match
    % ---------------------------------------
    for i_gt=1:num_gt

        num_detected_in_sigma = numel( find( sigma{ifile}(i_gt,:)>tr) );
        num_detected_in_tau = numel( find( tau{ifile}(i_gt,:)>tp) );

        if num_detected_in_sigma == 1 && num_detected_in_tau == 1
            recall_accum = recall_accum + 1.0;
            precision_accum = precision_accum + 1.0;

            % Mark the ground truth and detection, do not process twice
            flag_gt{ifile}(i_gt) = 1;
            idx_det = sigma{ifile}(i_gt,:)>tr;
            flag_det{ifile}(idx_det) = 1;
        end
    end

    % ---------------------------------------
    % check one-to-many match (split)
    % one gt with many detected rectangles
    % ---------------------------------------
    for i_gt=1:num_gt

        if flag_gt{ifile}(i_gt) > 0
            continue;
        end

        num_nonzero_in_sigma = sum( sigma{ifile}(i_gt,:)>0 );
        if num_nonzero_in_sigma >= k

            % -------------------------------------------------------------
            % Search the possible "many" partners for this	"one" rectangle
            % -------------------------------------------------------------

            % ----- satisfy 1st condition
            % only select unprocessed data
            idx_detected_in_tau = find( (tau{ifile}(i_gt,:)'>=tp) & (flag_det{ifile}==0) );
            num_detected_in_tau = numel( idx_detected_in_tau );

            if num_detected_in_tau == 1
                % Only one of the many-rectangles qualified ->
                % This match degraded to a one-to-one match
                if ( (tau{ifile}(i_gt, idx_detected_in_tau) >= tp) && ...
                        (sigma{ifile}(i_gt, idx_detected_in_tau) >= tr) )
                    recall_accum = recall_accum + 1.0;
                    precision_accum = precision_accum + 1.0;
                end
            else
                % satisfy 2nd condition
                if sum( sigma{ifile}(i_gt,idx_detected_in_tau) ) >= tr

                    % Mark the "one" rectangle
                    flag_gt{ifile}(i_gt) = 1;

                    % Mark all the "many" rectangles
                    flag_det{ifile}(idx_detected_in_tau) = 1;

                    recall_accum = recall_accum + fsc_k;
                    precision_accum = precision_accum + num_detected_in_tau * fsc_k;

                end
            end

        end

        % No match
        recall_accum = recall_accum + 0;
        precision_accum = precision_accum + 0;

    end

    % ---------------------------------------
    % check many-to-one match (merge)
    % one detected rectangle with many gt
    % ---------------------------------------
    for i_test=1:num_detected

        if flag_det{ifile}(i_test) > 0
            continue;
        end

        num_nonzero_in_tau = sum( tau{ifile}(:,i_test)>0 );
        if num_nonzero_in_tau >= k

            % satisfy 1st condition
            % only select unprocessed data
            idx_detected_in_sigma = find( (sigma{ifile}(:,i_test)>=tr) & (flag_gt{ifile}==0) );
            num_detected_in_sigma = numel( idx_detected_in_sigma );

            if num_detected_in_sigma == 1
                % Only one of the many-rectangles qualified ->
                % This match degraded to a one-to-one match
                if ( (tau{ifile}(idx_detected_in_sigma, i_test) >= tp) && ...
                        (sigma{ifile}(idx_detected_in_sigma, i_test) >= tr) )
                    recall_accum = recall_accum + 1.0;
                    precision_accum = precision_accum + 1.0;
                end
            else
                % satisfy 2nd condition
                if sum( tau{ifile}(idx_detected_in_sigma,i_test) ) >= tp
                    % Mark the "one" rectangle
                    flag_det{ifile}(i_test) = 1;

                    % Mark all the "many" rectangles
                    flag_gt{ifile}(idx_detected_in_sigma) = 1;

                    recall_accum = recall_accum + num_detected_in_sigma*fsc_k;
                    precision_accum = precision_accum + fsc_k;
%                     recall_accum = recall_accum + num_detected_in_sigma;
%                     precision_accum = precision_accum + 1.0;
                end
            end

        end

        % No match
        recall_accum = recall_accum + 0;
        precision_accum = precision_accum + 0;
    end

end
recall = recall_accum / tot_gt;
precision = precision_accum / tot_detected;

end

author: cs-chan

repo: Total-Text-Dataset


text:
function varargout = deadhub(varargin)

switch class(varargin{1})

    case 'double'
        error('Overloaded SDPVAR/SIN CALLED WITH DOUBLE. Report error')

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = struct('convexity','convex','monotonicity','none','definiteness','positive','model','callback');
        operator.bounds     = @bounds;
        operator.convexhull = @convexhull;

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error('SDPVAR/SIN called with CHAR argument?');
end

function [L,U] = bounds(xL,xU,lambda)

fL = deadhub(xL,lambda);
fU = deadhub(xU,lambda);
U = max(fL,fU);
L = min(fL,fU);
if xL<0 & xU>0
    L = 0;
end

function [Ax, Ay, b, K] = convexhull(xL,xU,lambda)

K.l = 0;
K.f = 0;
fL = deadhub(xL,lambda);
fU = deadhub(xU,lambda);
if xL>=-lambda & xU<=lambda
     Ax = 0;Ay = 1;b = 0;K.f = 1;
elseif xU < -3*lambda
    Ax = 1;Ay = 1;b = 2*lambda^2;K.f = 1;
elseif xL > 3*lambda
    Ax = -1;Ay = 1;b = 2*lambda^2;K.f = 1;
else
    dfL = derivative(xL,lambda);
    dfU = derivative(xU,lambda);
    [Ax,Ay,b,K] = convexhullConvex(xL,xU,fL,fU,dfL,dfU);
end

function df=derivative(x,lambda)
ax = abs(x);
if ax<lambda
    df=0;
elseif ax>3*lambda
    df = lambda;
elseif ax<=3*lambda
    df = 0.25*(2*ax-6*lambda)+lambda;
end
if x<0
    df=-df;
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function varargout = cumprod(varargin)
%CUMPROD  Indefinite product integral of a CHEBFUN2.
%   G = CUMPROD(F) returns the CHEBFUN2 G = exp( cumsum(log(F)) )
%
%   G = CUMPROD(F, DIM) returns the CHEBFUN2 G = exp( cumsum(log(F), DIM) )
%
% See also CUMSUM, SUM, PROD.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

[varargout{1:nargout}] = cumprod@separableApprox(varargin{:});

end


text:
function imgOut = TumblinTMO(img, L_da, Ld_Max, C_Max, L_wa)
%
%        imgOut = TumblinTMO(img, L_da, Ld_Max, C_Max, L_wa)
%
%
%        Input:
%           -img: an HDR image
%           -L_da: adaptation display luminance in [10,30] cd/m^2
%           -Ld_Max: maximum display luminance in [80, 180] cd/m^2
%           -C_Max: maximum LDR monitor contrast typically between 30 to 100
%
%        Output:
%           -imgOut: a tone mapped image in [0,1]
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%     The paper describing this technique is:
%     "Two Methods for Display of High Contrast Images"
% 	  by JACK TUMBLIN, JESSICA K. HODGINS, and BRIAN K. GUENTER
%     in ACM Transactions on Graphics, Vol. 18, No. 1, January 1999, Pages 56?94.
%

%is it a three color channels image?
check13Color(img);

checkNegative(img);

%check parameters
if(~exist('L_da', 'var'))
    L_da = 20;
end

if(~exist('Ld_Max', 'var'))
    Ld_Max = 100;
end

if(~exist('C_Max', 'var'))
    C_Max = 100;
end

if(~exist('L_wa', 'var'))
    L_wa = logMean(lum(img)); %luminance world adaptation
end

%compute luminance channel
L = lum(img);

%range compression
gamma_w = StevensCSF(L_wa);
gamma_d = StevensCSF(L_da);
gamma_wd = gamma_w / (1.855 + 0.4 * log(L_da));
m = C_Max.^((gamma_wd - 1) / 2.0);
Ld = L_da * m .* ((L ./ L_wa).^(gamma_w / gamma_d));
Ld = Ld / Ld_Max;

%change luminance
imgOut = ChangeLuminance(img, L, Ld);

end

author: banterle


text:
%% Demo of *plot_littlewood_paley_1d*

%% Usage
% littlewood = *plot_littlewood_paley_1d*(filters) (see
% <matlab:doc('plot_littlewood_paley_1d') plot_littlewood_paley_1d>).
%
%% Description
% *plot_littlewood_paley* computes, at every frequency, the
% Littlewood-Paley sum of a filter bank, i.e. the total power spectral
% density
% \sum_{j, \theta} |\hat{\psi_j} (\omega)|^2 + |\hat{\phi_J}(\omega)|^2
% If this sum is between $(1-epsilon)$ and $1$ for small $epsilon,
% the associated wavelet transform is proved to be contractive and
% almost unitary.

% In this demo, we display the Littlewood-Paley sum of a dyadic Morlet
% wavelet filter bank, with a very low averaging size of 8 samples. The
% Littlewood-Paley sum is shown in red, while the lowpass filter phi and
% the bandpass filters psi are respectively shown in green and blue.

figure;
T = 2^3;
interpolation = 2^10;
filt_opt.Q = 1;
filt_opt.J = T_to_J(T,filt_opt);
dyadic_filters = morlet_filter_bank_1d(T*interpolation,filt_opt);

plot_littlewood_paley_1d(dyadic_filters);
title('Q = 1 ; T = 8 samples (interpolated)');

% A more realistic example is constructed with an averaging size of 4096
% samples and a quality factor of 8. These values are typical in audio
% signal processing. The lowpass filter has such a narrow bandwidth that it
% is almost not visible in this second plot.

figure;
T = 2^12;
filt_opt.Q = 8;
filt_opt.J = T_to_J(T,filt_opt);
audio_filters = morlet_filter_bank_1d(T,filt_opt);

plot_littlewood_paley_1d(audio_filters);
title('Q = 8 ; T = 4096 samples');

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [ m, d ] = thanksgiving_canada ( y )

%*****************************************************************************80
%
%% THANKSGIVING_CANADA computes Canadian Thanksgiving for a Common year.
%
%  Discussion:
%
%    Canadian Thanksgiving occurs on the second Monday in October.
%
%  Example:
%
%    Input:
%
%      Y = 2002
%
%    Output:
%
%      M = 11
%      D = 28
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year.
%
%    Output, integer M, D, the month and day of Thanksgiving.
%

%
%  Determine the day of the week for 8 October, the earliest day
%  that Thanksgiving can occur.
%
  m = 10;
  d = 8;
  f = 0.0;

  w = ymdf_to_weekday_common ( y, m, d, f );
%
%  If W = 2 means this day is Monday, and day D is Thanksgiving.
%  Otherwise, figure out how to increment W to 2;
%  The same increment makes D the correct day number.
%
  if ( w < 2 )
    d = d + 2 - w;
  elseif ( 2 < w )
    d = d + 2 + 7 - w;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%ControlCenter
mex SearchCycleNode.c

n=30; % The number of nodes in the tree

CostMatric = rand( n );  % randomly generating a cost matrix
for p=1:n
    CostMatric( p,p ) = 0;
end

% 1. search maximum directed spanning tree by unsymmetrical cost matrix "CostMatric", with specified root node "Root".

Root = 3;    % Root of the tree is predefined ahead .
[MaxTree1,MaxCost1] =  DirectedMaximumSpanningTree( CostMatric,Root )
h = view(biograph( MaxTree1 ));

% 2. % 1. search maximum directed spanning tree by unsymmetrical cost matrix "CostMatric", with no specified root node "Root".

[MaxTree2,MaxCost2] =  MaximalDirectedMSF( CostMatric )
h = view(biograph( MaxTree2 ));

% 3. search minimum directed spanning tree by unsymmetrical cost matrix "CostMatric", with specified root node "Root".

Root = 3;    % Root of the tree is predefined ahead .
[MaxTree3,MaxCost3] =  DirectedMinimalSpanningTree( CostMatric,Root )
h = view(biograph( MaxTree3 ));

% 2. % 1. search minimum directed spanning tree by unsymmetrical cost matrix "CostMatric", with no specified root node "Root".

[MaxTree4,MaxCost4] =  MinimalDirectedMSF( CostMatric )
h = view(biograph( MaxTree4 ));


text:
function A = resize(B,R,C)
% A = resize(B,R,C)   Crop or zero-pad B to have R rows and C columns.
%	I'm sure this must already be provided, but how to know?
% dpwe 1995jan21

% Copyright (c) 1995 Dan Ellis <dpwe@ee.columbia.edu>
% released under GPL - see file COPYRIGHT

A = zeros(R,C);
[r,c] = size(B);

mr = min(r,R);
mc = min(c,C);

A(1:mr,1:mc) = B(1:mr, 1:mc);

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code was written by Marcello Torchio, University of Pavia.
% Please send comments or questions to
% marcello.torchio01@ateneopv.it
%
% Copyright 2017: 	Marcello Torchio, Lalo Magni, and Davide M. Raimondo, University of Pavia
%					Bhushan Gopaluni, University of British Columbia
%                 	Richard D. Braatz, MIT.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RETREIVEDATA retreives the data which is returned in the results structure.

function [ce_t, cs_barrato_t, T_t, jflux_t, Phis_t, Phie_t, cs_star_t, SOC_t, film_t, js_t, Q_t,t_tot] = retreiveData(ce_t, cs_barrato_t, T_t, jflux_t, Phis_t, Phie_t,cs_star_t, SOC_t, film_t, js_t, Q_t,t_tot, y, t, param)
	% Extract differential variables after the integration process
	ce_t            = [ce_t;y(param.ce_indices)];
	cs_barrato_t    = [cs_barrato_t;y(param.cs_average_indices)];
	T_t             = [T_t;y(param.T_indices)];
	film_t          = [film_t;y(param.film_indices)];
    Q_t             = [Q_t;y(param.Q_indices)];
	% Extract the algebraic variables after the integration process
	jflux_t         = [jflux_t;y(param.jflux_indices)];
	Phis_t          = [Phis_t;y(param.Phis_indices)];
	Phie_t          = [Phie_t;y(param.Phie_indices)];
	js_t            = [js_t;y(param.js_indices)];

	% Check if Fick's law of diffusion is used. This is required to define the
    % correct way how to evaluate the SOC.
    if(param.SolidPhaseDiffusion~=3)
        cs_average = cs_barrato_t(end,param.Np+1:end);
    else
        start_index = param.Nr_p*param.Np+1;
        end_index   = start_index+param.Nr_n-1;
        cs_average  = zeros(param.Nn,1);
        for n=1:param.Nn
            cs_average(n)   = 1/param.Rp_n*(param.Rp_n/param.Nr_n)*sum(cs_barrato_t(end,start_index:end_index));
            start_index     = end_index + 1;
            end_index       = end_index + param.Nr_n;
        end
    end
    Csout           = sum(cs_average);
    Sout            = 100*(1/param.len_n*(param.len_n/(param.Nn))*Csout/param.cs_maxn);
    cs_star_t       = [cs_star_t;surfaceConcentration(cs_barrato_t(end,:)',jflux_t(end,:)',Q_t(end,:)',T_t(end,:)',param)'];
    SOC_t           = [SOC_t;Sout];
    t_tot           = [t_tot;t];
end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function gmm_12 = product_gmm(gmm_1,gmm_2,mrf)
%PRODUCT_GMM takes the product between two GMMs

K1 = size(gmm_1.Priors,2);

k = 1;

for i=1:K1

    I = find(mrf.A(i,:));

    if ~isempty(I)

        jl = [mrf.Aj_index(i,I,1)' mrf.Aj_index(i,I,2)'];

        for j=1:size(jl,1)

                [Mu,Sigma] = product_gauss(gmm_1.Mu(:,i),gmm_1.Sigma(:,:,i),gmm_2.Mu(:,jl(2)),gmm_2.Sigma(:,:,jl(2)));
                gmm_12.Mu(:,k) = Mu;
                gmm_12.Sigma(:,:,k) = Sigma;
                k = k + 1;
        end

    end
end

K = size(gmm_12.Mu,2);
gmm_12.Priors = ones(1,K)./K;

end

author: epfl-lasa


text:
%
% Author:  Stilian Stoev (C), sstoev@math.bu.edu
%
% This function computes the discrete convolution of the vectors
% a and b, efficiently, by using the FFT algorithm.
%
% input:
%  a <- a vector
%  b <- a vector with the same dimensions as a
%  force <- if force == 1, then the FFT's are "forced" to be
%           of dyadic complexity.
%
% output:
%  c <- the discrete convolution of a and b, that is,
%       (a(1)*b(1), a(1)*b(2)+a(2)*b(1), ...)
%
% usage:
%   c = fftconv(a,b,force);
%
% Written by sstoev@math.bu.edu
%
function c = fftconv(a,b,force);

na = length(a);
nb = length(b);

if force,
 n = 2^(fix(log2(na+nb))+1);
else
 n=na+nb;
end;
A = fft(a,n);
B = fft(b,n);

c = ifft(A.*B,n);
c = real(c(1:na+nb-1));


text:
function [err_p,elerr_p] = cdpost_bc(viscosity,aez,fez,elerror,xy,ev,ebound)
%cdpost_bc postprocesses local Poisson error estimator
%   [err_p,elerr_p] = cdpost_bc(viscosity,aez,fez,elerror,xy,ev,ebound);
%   input
%          viscosity   viscosity parameter
%          aez         elementwise Poisson problem matrices
%          fez         elementwise rhs vectors
%          elerror     elementwise error estimate (without BC imposition)
%          xy          vertex coordinate vector
%          ev          element mapping matrix
%          ebound      element edge boundary matrix
%   output
%          err_p       global error estimate
%          elerr_p     elementwise error estimate
%   IFISS function: DJS; 5 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
      x=xy(:,1); y=xy(:,2);
      nel=length(ev(:,1));
      lev=[ev,ev(:,1)]; elerr_p=elerror;
%
% recompute contributions from elements with Dirichlet boundaries
      nbde=length(ebound(:,1));
      ebdy = zeros(nel,1);
      edge = zeros(nel,1);
% isolate boundary elements
      for el = 1:nbde
      ee = ebound(el,1);
      ebdy(ee) = ebdy(ee)+1; edge(ee)=ebound(el,2);
      end
%
% two edge elements
      k2=find(ebdy==2);
      nel2b=length(k2);
% loop over two edge elements
      for el = 1:nel2b
      el2e=k2(el);
      kk=find(ebound(:,1) == el2e);
      edges=ebound(kk,2);
% set up original matrix and RHS vector
	  ae=squeeze(aez(el2e,1:5,1:5));
      fe=fez(el2e,:)';
% set up local coordinates and impose interpolated error as Dirichlet bc
      xl=x(lev(el2e,:)); yl=y(lev(el2e,:));
      [bae,fe] = localbc_pcd(ae,fe,edges,xl,yl);
% solve local problem
        err=(viscosity*bae)\fe;
        elerr_p(el2e,1) = err'*fe/viscosity;
      end
% end of element loop
%
% one edge elements
      k1=find(ebdy==1);
      nel1b=length(k1);
% loop over one edge elements
      for el = 1:nel1b
      el1e=k1(el);
      kk=find(ebound(:,1) == el1e);
      edges=ebound(kk,2);
% set up original matrix and RHS vector
      fe=fez(el1e,:)';
	  ae=squeeze(aez(el1e,1:5,1:5));
% set up local coordinates and impose interpolated error as Dirichlet bc
      xl=x(lev(el1e,:)); yl=y(lev(el1e,:));
      [bae,fe] = localbc_pcd(ae,fe,edges,xl,yl);
% solve for local estimate
        err=(viscosity*bae)\fe;
        elerr_p(el1e,1) = err'*fe/viscosity;
      end
% end of element loop
%
      err_p = sqrt(sum(elerr_p));
      elerr_p = sqrt(elerr_p);
      fprintf('done\n')
      fprintf('estimated global error (in energy):  %10.6e\n',err_p)
 return


text:
function [Cp1,Cp2] = Cpre_q1q1(xy,ev)
%Cpre_q1q1  generate stabilizations for least sqrs commutator for Q1-Q1
%   [Cp1,Cp2] = Cpre_q1q1(xy,ev);
%   input
%          xy         Q2 nodal coordinate vector
%          ev         element mapping matrix
%   output
%          Cp1        pressure stabilization 1 for preconditioner
%          Cp2        pressure stabilization 2 for preconditioner
%   IFISS function: HCE; 7 July 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

nngpt=4;
x=xy(:,1); y=xy(:,2);
xp=xy(:,1); yp=xy(:,2);
nvtx=length(x); nu=2*nvtx; np=length(xp);
nel=length(ev(:,1));
%
% initialise global matrices
Cp1 = sparse(np,np);
Cp2 = sparse(np,np);
%
% Gauss point integration rules
if (nngpt==4)        % 2x2 Gauss points
   gpt=1.0e0/sqrt(3.0e0);
   s(1) = -gpt; t(1) = -gpt; wt(1)=1;
   s(2) =  gpt; t(2) = -gpt; wt(2)=1;
   s(3) =  gpt; t(3) =  gpt; wt(3)=1;
   s(4) = -gpt; t(4) =  gpt; wt(4)=1;
elseif (nngpt==1)   % 1x1 Gauss point
   s(1) =    0; t(1) =    0; wt(1)=4;
else
   error('Check Gauss point integration specification')
end
%
% inner loop over elements
for ivtx = 1:4
   xl_v(:,ivtx) = x(ev(:,ivtx));
   yl_v(:,ivtx) = y(ev(:,ivtx));
end

%
% element assembly into global matrices
cpe1 = zeros(nel,4,4);
cpe2 = zeros(nel,4,4);

% loop over Gauss points
for igpt = 1:nngpt
   sigpt=s(igpt);
   tigpt=t(igpt);
   wght=wt(igpt);
   [jac,invjac,phi,dphidx,dphidy] = deriv(sigpt,tigpt,xl_v,yl_v);
   for j = 1:4
      for i = 1:4
          cpe1(:,i,j) = cpe1(:,i,j) + .25 * wght*(phi(:,i)-0.25) .*(phi(:,j)-0.25) ;
          cpe2(:,i,j) = cpe2(:,i,j) +  (16*jac(:)) .\ (wght*(phi(:,i)-0.25) .*(phi(:,j)-0.25));
      end
   end
% end of Gauss point loop
end

%%  element assembly into global matrices
for krow=1:4
   nrow=ev(:,krow);
   for kcol=1:4
      ncol=ev(:,kcol);
      Cp1 = Cp1 + sparse(nrow,ncol,cpe1(:,krow,kcol),nvtx,nvtx);
      Cp2 = Cp2 + sparse(nrow,ncol,cpe2(:,krow,kcol),nvtx,nvtx);
   end
end


text:
function g = p26_g ( n, x )

%*****************************************************************************80
%
%% P26_G evaluates the gradient for problem 26.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real G(N), the gradient of the objective function.
%
  g = zeros ( n, 1 );

  jroot = 5;
  k = 500;
  fi = k;
  dfidx1 = 0.0;
  dfidx2 = 0.0;

  for j = 1 : jroot * jroot

    j1 = mod ( j-1, jroot ) + 1;
    a1 = -32 + j1 * 16;

    j2 = ( j - 1 ) / jroot;
    a2 = -32 + j2 * 16;

    fj = ( j + ( x(1) - a1 )^6 + ( x(2) - a2 )^6 );

    fi = fi + 1.0 / fj;

    dfidx1 = dfidx1 - ( 1.0 / fj^2 ) * 6.0 * ( x(1) - a1 )^5;
    dfidx2 = dfidx2 - ( 1.0 / fj^2 ) * 6.0 * ( x(2) - a2 )^5;

  end

  g(1) = - ( 1.0 / fi^2 ) * dfidx1;
  g(2) = - ( 1.0 / fi^2 ) * dfidx2;

  return
end


text:
function region = p31_region ( )

%*****************************************************************************80
%
%% P31_REGION returns the name of the integration region for problem 31.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string REGION, the name of the integration region.
%
  region = 'box';

  return
end


text:
% fastranc.m -  adaptive noise cancellation for fmrib_fastr.m
%   [out,y]=fmribanc(refs,d,N,mu)
%
%   Inputs:
%       refs: Reference data
%       d:   Input data
%       N:   filter order (length)
%       mu:  Step size
%
%   Outputs:
%       out: error signal (signal after noise cancellation)
%       y: filtered out noise
%
%   If a binary version of this program is present in same directory
%   it will be run instead.  The binary version is MUCH fastr.  However,
%   if it is not present or unavailable for your platform, you can use
%   the matlab compiler to compile the C Mex-File  fastranc.c by typing
%   >> mex fastranc.c
%   from your command line after you CD into the fmrib1.0 plugin directory.
%   This will create the appropriate binary file for your platform
%   (if you have the MATLAB compiler). You can check what binaries are
%   available by looking at the extention of the fastranc.* files:
%   fastranc.mexa64   -  Linux 64 bit (compiled on Xeon EM64T but might
%                          be the same for athalon 64)
%   fastranc.mexaxp   -  Alpha 64 processor,
%   fastranc.mexglx   -  Linux 32 bit processor.
%   fastranc.dll      -  Windows.
%
%   Check your matlab help for more info on compiling MEX files.
%
%
%   Author:  Rami K. Niazy
%
%   Copyright (c) 2004 University of Oxford

% Copyright (C) 2004 University of Oxford
% Author:   Rami K. Niazy, FMRIB Centre
%           rami@fmrib.ox.ac.uk
%
% 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

function [out,y]=fastranc(refs,d,N,mu)

nargchk(4,4,nargin);
refs=refs(:);
d=d(:);
mANC=length(d);
if length(d)~=length(refs)
error('Reference and Input data must be of the same length','fmribanc() error!');
end

W=zeros(N+1,1);
r=flipud([0;refs(1:N)]);
out=zeros(mANC,1);
y=zeros(mANC,1);

for E=N+1:mANC

    %-----------------calc------------------------------------------

    r=[refs(E);r(1:end-1)];
    y(E)=sum(W.*r);
    out(E)=d(E)-y(E);
    W=W+2*mu*out(E)*r;

    %---------------------------------------------------------------
end
return;

author: PatternRecognition

repo: OpenBMI


text:
function [p1,p2]=rbline(varargin)
%function to draw a rubberband line and return the start and end points
%Usage: [p1,p2]=rbline;     uses current axes
% or    [p1,p2]=rbline(h);  uses axes refered to by handle, h

% Created by Sandra Martinka March 2002

switch nargin
case 0
  h=gca;
case 1
  h=varargin{1};
  axes(h);
otherwise
  disp('Too many input arguments.');
end
cudata=get(gcf,'UserData'); %current UserData
hold on;
k=waitforbuttonpress;
p1=get(h,'CurrentPoint');       %get starting point
p1=p1(1,1:2);                   %extract x and y
lh=plot(p1(1),p1(2),'+:');      %plot starting point
udata.p1=p1;
udata.h=h;
udata.lh=lh;
set(gcf,'UserData',udata,'WindowButtonMotionFcn','wbmf','DoubleBuffer','on');
k=waitforbuttonpress;
p2=get(h,'Currentpoint');       %get end point
p2=p2(1,1:2);                   %extract x and y
set(gcf,'UserData',cudata,'WindowButtonMotionFcn','','DoubleBuffer','off'); %reset UserData, etc..
delete(lh);

function wbmf
%window motion callback function

utemp=get(gcf,'UserData');
ptemp=get(utemp.h,'CurrentPoint');
ptemp=ptemp(1,1:2);
set(utemp.lh,'XData',[utemp.p1(1),ptemp(1)],'YData',[utemp.p1(2),ptemp(2)]);


text:
% This M file creates the "LossSpec_Diode_Library.mat " file
% containing loss and thermal diode models provided in the
% "LossModelib.mdl' libray.
% Diode specifications are saved in the "LossSpec_Diode" structure.
%
% Data provided below for three different commercial
% IGBT/Diode half-bridge packs are extracted from manufacturer data sheets
% (see pdf files provided in the same directory).

% Pierre Giroux, Hydro-Quebec (IREQ), March 2012

%% How to to add your own diode specifications
% -------------------------------------------
%   1)  Open the "LossModelib.mdl' library
%   2)  Duplicate a cell containing specifications of one existing IGBT
%       and insert it at the end of this file
%       (see "Add your own specifications here" cell)
%   3)  Modify this new set of data according to your specifications
%   4)  Run this script file in order to update the "LossSpec_Diode_Library.mat " file
%   5)  In the "LossModelib.mdl' library
%            - select the IGBT or half-bridge model to update
%            - open the Mask Editor
%            - select "Diode_Type" parameter
%            - In the "Popups" window, enter a new description
%              after existing diode descriptions
%   6) Save the "LossModelib.mdl' library.
% Note: Future release of SPS will provide an easy-to-use GUI to add or
% edit device characteristics.

% Pierre Giroux, Hydro-Quebec (IREQ), March 2012

% If=  Forward current (A)
% Vf= Forward voltage (V)
% Erec= Reverse recovery energy (mJ)
% Tj= Junction temperature (degrees Celsius)
% Vcc= Supply voltage (V)
% Note1: Reverse recovery energy (Erec)is function of Vcc, If and Tj
% Note2: Diode on-state losses = Vf (function of If and Tj) * If
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
clear all

%% Fuji Electric - Diode, IGBT Module, 2 in one-package, 600V, 150A

Type         = 'Diode';
Manufacturer = 'Fuji Electric';
PartNo       = '2MBI150U2A 060';
Description  = 'Diode, IGBT Module, 2 in one-package, 600V, 150A';

% Typical reverse recovery energy (Erec) vs If
Vcc_Erec = 300;
Tj_Erec  = [ 25 125 ];
If_Erec  = [ 0  25    50    75    100   125   150   175   200    225  ];
Erec     = [ 0  0.27  0.48  0.62  0.69  0.73  0.75  0.77  0.795  0.81
             0  0.48  0.82  1.07  1.23  1.37  1.47  1.58  1.67   1.79 ];
%
% Diode - Typical on-state characteristics
Tj_OnState=  [ 25 125 ];
If_OnState = [ 0  0.1  3     10    30    70    111   167   240   300  ];
Vf_OnState = [ 0  0.43 0.78  0.92  1.11  1.31  1.47  1.66  1.89  2.07
               0  0.42 0.55  0.70  0.96  1.23  1.47  1.72  2.02  2.25 ];
%-------------------------------------------------------------------------
%
% Thermal Impedance:
%
Rth_jc=0.46; % Junction-to-Case thermal resistance(K/W)
% Based on transient thermal impedance vs. time curve
Cth_j=0.13;    % Junction thermal capacitance (Joule/Kelvin)
%
%-------------------------------------------------------------------------
LossSpec_Diode(1)= struct('Type',Type, 'Manufacturer',Manufacturer, 'PartNo',PartNo, ...
    'Description',Description, 'Vcc_Erec',Vcc_Erec, ...
    'Tj_Erec',Tj_Erec, 'If_Erec',If_Erec, 'Erec',Erec, ...
    'Tj_OnState',Tj_OnState, 'If_OnState',If_OnState, 'Vf_OnState',Vf_OnState, ...
    'Rth_jc',Rth_jc, 'Cth_j',Cth_j);
%
clear Vcc_Erec Tj_Erec If_Erec Erec Tj_OnState ...
      If_OnState Vf_OnState Manufacturer PartNo Description Type ...
      Rth_jc Cth_j
%=========================end of specifications===========================

%% ABB - Diode, IGBT Module, ABB HiPak, 1700V, 800A
Type='Diode';
Manufacturer= 'ABB';
PartNo= '5SNE 0800M170100';
Description= 'Diode, IGBT Module, ABB HiPak, 1700V, 800A';

% For Tj=25 deg. C, manufacturer specifications provide Erec at only one current (800 A).
% The Erec values for Tj=25 deg. C and Tj=125 deg. C, 800 A specified below
% are used to deduce the Erec(Ic) curve @ Tj=25 deg. C from the Erec(Ic) curve @ Tj=125 deg. C
% assuming that Erec stays proportionnal to its value  @ Tj= 125 deg. C

Tj_800A=   [ 25   125 ]; % Junction temperature (C)
Erec_800A= [ 150  270];  % Reverse recovery energy  at Tj_800A (mJ) for Ic=800A
%
% Typical reverse recovery energy (Erec) vs If
Vcc_Erec= 900;
Tj_Erec=  [ 25 125 ];
If_Erec=   [ 0 101  196  301  500  703  800 902  1100 1300 1500 1600 ];
% Erec @ 125 C (per data sheet curves)
Erec(2,:)= [ 0 82.0 115  145  201  245  270 282  309  330  340  342  ];
% Eon @ 25 C (calculated using data sheet spec at 25 C)
Erec(1,:)=Erec_800A(1)/Erec_800A(2) * Erec(2,:);
%
% Diode - Typical on-state characteristics
Tj_OnState=  [ 25 125 ];
If_OnState = [ 0  0.01  96.4 142  195  259  365  500  669  900  1100 1400 1500];
Vf_OnState = [ 0  0.80  1.09 1.17 1.23 1.29 1.38 1.48 1.58 1.71 1.83 2.00 2.06
               0  0.50  0.86 0.97 1.06 1.15 1.28 1.43 1.58 1.78 1.94 2.18 2.26 ];
%-------------------------------------------------------------------------
%
% Thermal Impedance:
%
Rth_jc=0.036; % Junction-to-Case thermal resistance(K/W)
% Based on transient thermal impedance vs. time curve
Cth_j=3.729;    % Junction thermal capacitance (Joule/Kelvin)
%
%-------------------------------------------------------------------------
LossSpec_Diode(2)= struct('Type',Type, 'Manufacturer',Manufacturer, 'PartNo',PartNo, ...
    'Description',Description, 'Vcc_Erec',Vcc_Erec, ...
    'Tj_Erec',Tj_Erec, 'If_Erec',If_Erec, 'Erec',Erec, ...
    'Tj_OnState',Tj_OnState, 'If_OnState',If_OnState, 'Vf_OnState',Vf_OnState, ...
    'Rth_jc',Rth_jc, 'Cth_j',Cth_j);
%
clear Vcc_Erec Tj_Erec If_Erec Erec Tj_OnState ...
      If_OnState Vf_OnState Manufacturer PartNo Description Type ...
      Rth_jc Cth_j Tj_800A Erec_800A
%=========================end of specifications===========================

%% ABB - Diode, Half-bridge IGBT Module, 3300V, 250A;
Type='Diode';
Manufacturer= 'ABB';
PartNo= '5SNG 0250P330300';
Description= 'Half-bridge IGBT Module, 3300V, 250A';

% For Tj=25 deg. C, manufacturer specifications provide Erec at only one current (250 A).
% The Erec values for Tj=25 deg. C and Tj=125 deg. C, 250 A specified below
% are used to deduce the Erec(Ic) curve @ Tj=25 deg. C from the Erec(Ic) curve @ Tj=125 deg. C
% assuming that Erec stays proportionnal to its value  @ Tj= 125 deg. C

Tj_250A=   [ 25   125 ];  % Junction temperature (C)
Erec_250A= [ 165  280 ];  % Reverse recovery energy at Tj_250A (mJ) for If=250A
%
% Typical reverse recovery energy (Erec) vs If
Vcc_Erec= 1800;
Tj_Erec=  [ 25 125 ];
If_Erec=   [ 0  26.2  48.9  99   150  199  250  299  350  400  450  ];
% Erec @ 125 C (per data sheet curves)
Erec(2,:)= [ 0  85.9  113   165  211  248  280  306  325  336  340  ];
% Eon @ 25 C (calculated using data sheet spec at 25 C)
Erec(1,:)=Erec_250A(1)/Erec_250A(2) * Erec(2,:);
%
% Diode - Typical on-state characteristics
Tj_OnState=  [ 25 125 ];
If_OnState = [ 0  0.01  14.0  43    66.5  101  162   250   401   500  ];
Vf_OnState = [ 0  0.75  1.13  1.40  1.53  1.65 1.81  2.00  2.25  2.39
               0  0.50  0.90  1.20  1.37  1.55 1.81  2.11  2.50  2.71 ];
%-------------------------------------------------------------------------
%
% Thermal Impedance:
%
Rth_jc=0.102; % Junction-to-Case thermal resistance(K/W)
% Based on transient thermal impedance vs. time curve
Cth_j=1.3235;    % Junction thermal capacitance (Joule/Kelvin)
%
%-------------------------------------------------------------------------
LossSpec_Diode(3)= struct('Type',Type, 'Manufacturer',Manufacturer, 'PartNo',PartNo, ...
    'Description',Description, 'Vcc_Erec',Vcc_Erec, ...
    'Tj_Erec',Tj_Erec, 'If_Erec',If_Erec, 'Erec',Erec, ...
    'Tj_OnState',Tj_OnState, 'If_OnState',If_OnState, 'Vf_OnState',Vf_OnState, ...
    'Rth_jc',Rth_jc, 'Cth_j',Cth_j);
%
clear Vcc_Erec Tj_Erec If_Erec Erec Tj_OnState ...
      If_OnState Vf_OnState Manufacturer PartNo Description Type ...
      Rth_jc Cth_j Tj_250A Erec_250A
%=========================end of specifications===========================

%% Add your own specifications here
% Type='Diode';
% Manufacturer= 'My Manufacturer';
% PartNo= 'XXX';
% Description= 'My Module, YYY V, ZZZ A;
% ....

%% Save diode thermal data structure "LossSpec_Diode" in "LossSpec_Diode_Library.mat" file
save LossSpec_Diode_Library LossSpec_Diode
%


text:
function Offspring = Operator(Problem,Parent1,Parent2)
% Crossover and mutation used in MOEA/D-M2M
% proM --- 1 --- The expectation of number of bits doing mutation

%------------------------------- 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".
%--------------------------------------------------------------------------

    %% Parameter setting
    Parent1 = Parent1.decs;
    Parent2 = Parent2.decs;
    [N,D]   = size(Parent1);

    %% Crossover
    rc = (2*rand(N,1)-1).*(1-rand(N,1).^(-(1-Problem.FE/Problem.maxFE).^0.7));
    OffDec = Parent1 + repmat(rc,1,D).*(Parent1-Parent2);

    %% Mutation
    rm    = 0.25*(2*rand(N,D)-1).*(1-rand(N,D).^(-(1-Problem.FE/Problem.maxFE).^0.7));
    Site  = rand(N,D) < 1/D;
    Lower = repmat(Problem.lower,N,1);
    Upper = repmat(Problem.upper,N,1);
    OffDec(Site) = OffDec(Site) + rm(Site).*(Upper(Site)-Lower(Site));

	%% Set the infeasible decision variables to feasible values
    temp1 = OffDec < Lower;
    temp2 = OffDec > Upper;
    rnd   = rand(N,D);
    OffDec(temp1) = Lower(temp1) + 0.5*rnd(temp1).*(Parent1(temp1)-Lower(temp1));
    OffDec(temp2) = Upper(temp2) - 0.5*rnd(temp2).*(Upper(temp2)-Parent1(temp2));
    Offspring     = Problem.Evaluation(OffDec);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function T = compute_tensor_field_random(n,options)

% compute_tensor_field_random - create 2D TF
%
%   T = compute_tensor_field_random(n,options);
%
%   Copyright (c) 2008 Gabriel Peyre

options.null = 0;
sigma_flow = getoptions(options, 'sigma_tensor', 50*n/256);
niter_tensor = getoptions(options, 'niter_tensor', 8);
verb = getoptions(options, 'verb', 1);

T = randn(n,n,3);
en = linspace(.5,1,n^2);
an = linspace(.1,1,n^2);
or = linspace(-pi/2,pi/2,n^2);
for i=1:niter_tensor
    if verb
        progressbar(i,niter_tensor);
    end
    T = perform_blurring(T,sigma_flow);
    U = perform_tensor_mapping(T,+1);
%    clf;
%   for k=1:3
%        subplot(1,3,k); a = U(:,:,k);
%        hist(a(:),100); axis tight;
%    end
%    drawnow;
    U(:,:,1) = perform_histogram_equalization(U(:,:,1), en);
    U(:,:,2) = perform_histogram_equalization(U(:,:,2), an);
    U(:,:,3) = perform_histogram_equalization(U(:,:,3), or);
    T = perform_tensor_mapping(U,-1);
end

author: gpeyre

repo: matlab-toolboxes


text:
function [mappedX, mapping] = charting(X, no_dims, no_analyzers, max_iterations, eig_impl)
%CHARTING Performs manifold charting on dataset X
%
%   [mappedX, mapping] = charting(X, no_dims, no_analyzers, max_iterations, eig_impl)
%
% Performs manifold charting on dataset X to reduce its dimensionality to
% no_dims dimensions. The variable no_analyzers determines the number of
% local factor analyzers that is used in the mixture of factor analyzers
% (default = 40). The variable max_iterations sets the maximum number of
% iterations that is employed in the EM algorithm (default = 200).
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    if ~exist('no_dims', 'var')
        no_dims = 2;
    end
    if ~exist('no_analyzers', 'var')
        no_analyzers = 40;
    end
    if ~exist('max_iterations', 'var')
        max_iterations = 200;
    end

    % Initialize some parameters
    tol = 1e-10;                        % regularization parameter
    min_std = 1e-3;                     % minimum STD of Gaussians
    kf = no_analyzers * (no_dims + 1);

    % Construct MFA model on the data
    disp('Running EM algorithm and compute local factor analyzers...');
    [LX, MX, PX] = mppca(X', no_dims, no_analyzers, tol, max_iterations, min_std);
    [R, Z] = infermfa(X', LX, MX, PX);

    % Adds last entry = 1 in posterior mean to handle means of factor analyzers
    Z(no_dims + 1,:,:) = 1;
    Z = permute(Z, [1 3 2]);

    % Construct blockdiagonal matrix D
    disp('Performing manifold charting...');
    D = zeros((no_dims + 1) * no_analyzers, (no_dims + 1) * no_analyzers);
    for i=1:no_analyzers
        Ds = zeros(no_dims + 1, no_dims + 1);
        for j=1:size(X, 1)
            Ds = Ds + R(i, j) .* (Z(:,i,j) * Z(:,i,j)');
        end
        D((i - 1) * (no_dims + 1) + 1:i * (no_dims + 1), (i - 1) * (no_dims + 1) + 1:i * (no_dims + 1)) = Ds;
    end

    % Construct responsibility weighted local representation matrix U
    R = reshape(R, [1 no_analyzers size(X, 1)]);
    U = reshape(bsxfun(@times, R, Z), [kf size(X, 1)])';

    % Solve generalized eigenproblem
    if strcmp(eig_impl, 'Matlab')
        options.disp = 0;
        options.isreal = 1;
        options.issym = 1;
        [V, lambda] = eigs(D - U' * U, U' * U, no_dims + 1, 'SM', options);
    else
        options.Disp = 0;
        options.LSolver = 'bicgstab';
        [V, lambda] = jdqz(D - U' * U, U' * U, no_dims + 1, 'SM', options);
    end
    [lambda, ind] = sort(diag(lambda));
    V = V(:,ind(2:end));

    % Compute final lowdimensional data representation
    mappedX = U * V;

    % Store mapping data for out-of-sample extension
    if nargout > 1
        mapping.LX = LX;
        mapping.MX = MX;
        mapping.PX = PX;
        mapping.V = V;
    end

repo: cluster