text:
function bnet  = mk_incinerator_bnet(ns)
% MK_INCINERATOR_BNET The waste incinerator emissions example from Cowell et al p145
% function bnet  = mk_incinerator_bnet(ns)
%
% If ns is omitted, we use the scalars and binary nodes and the original params.
% Otherwise, we use random params of the desired size.
%
% Lauritzen, "Propogation of Probabilities, Means and Variances in Mixed Graphical Association Models",
% JASA 87(420): 1098--1108
% This example is reprinted on p145 of "Probabilistic Networks and Expert Systems",
% Cowell, Dawid, Lauritzen and Spiegelhalter, 1999, Springer.
% For a picture, see http://www.cs.berkeley.edu/~murphyk/Bayes/usage.html#cg_model

% node numbers
F = 1; W = 2; E = 3; B = 4; C = 5; D = 6; Min = 7; Mout = 8; L = 9;
names = {'F', 'W', 'E', 'B', 'C', 'D', 'Min', 'Mout', 'L'};
n = 9;
dnodes = [F W B];
cnodes = mysetdiff(1:n, dnodes);

% node sizes - all cts nodes are scalar, all discrete nodes are binary
if nargin < 1
  ns = ones(1, n);
  ns(dnodes) = 2;
  rnd = 0;
else
  rnd = 1;
end

% topology (p 1099, fig 1)
dag = zeros(n);
dag(F,E)=1;
dag(W,[E Min D]) = 1;
dag(E,D)=1;
dag(B,[C D])=1;
dag(D,[L Mout])=1;
dag(Min,Mout)=1;

% params (p 1102)
bnet = mk_bnet(dag, ns, 'discrete', dnodes, 'names', names);

if rnd
  for i=dnodes(:)'
    bnet.CPD{i} = tabular_CPD(bnet, i);
  end
  for i=cnodes(:)'
    bnet.CPD{i} = gaussian_CPD(bnet, i);
  end
else
  bnet.CPD{B} = tabular_CPD(bnet, B, 'CPT', [0.85 0.15]); % 1=stable, 2=unstable
  bnet.CPD{F} = tabular_CPD(bnet, F, 'CPT', [0.95 0.05]); % 1=intact, 2=defect
  bnet.CPD{W} = tabular_CPD(bnet, W, 'CPT', [2/7 5/7]); % 1=industrial, 2=household
  bnet.CPD{E} = gaussian_CPD(bnet, E, 'mean', [-3.9 -0.4 -3.2 -0.5], ...
			     'cov', [0.00002 0.0001 0.00002 0.0001]);
  bnet.CPD{D} = gaussian_CPD(bnet, D, 'mean', [6.5 6.0 7.5 7.0], ...
			     'cov', [0.03 0.04 0.1 0.1], 'weights', [1 1 1 1]);
  bnet.CPD{C} = gaussian_CPD(bnet, C, 'mean', [-2 -1], 'cov', [0.1 0.3]);
  bnet.CPD{L} = gaussian_CPD(bnet, L, 'mean', 3, 'cov', 0.25, 'weights', -0.5);
  bnet.CPD{Min} = gaussian_CPD(bnet, Min, 'mean', [0.5 -0.5], 'cov', [0.01 0.005]);
  bnet.CPD{Mout} = gaussian_CPD(bnet, Mout, 'mean', 0, 'cov', 0.002, 'weights', [1 1]);
end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
 function argpair = fatrix2_construct_args(ob, varargin)
%function argpair = fatrix2_construct_args(ob, options)
%|
%| Given a fatrix2_block object, return a cell array of name, value pairs
%| such that the object could be re-constructed by
%|	ob = fatrix2(argpair{:})
%| This is useful for constructing variations of the object.
%|
%| out
%|	argpair {2n}	{name1, value1, ...}
%|
%| Copyright 2010-12-02, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end

% user-definable properties
name_user = {'arg', 'dim', 'caller', 'odim', 'mask', 'meth'};

% private/internal properties
name_priv = {'dim', 'is_transpose', 'nblock', 'iblock'}

% handles
name_hand = {'back', 'forw', 'power', 'abs', 'free', 'gram', ...
	'mtimes_block', 'block_setup', 'blockify_data'};

argpair = {};

for ii=1:numel(name_user)
	name = name_user{ii};
	argpair{end+[1:2]} = {name, ob.(name)};
end

for ii=1:numel(name_priv)
	name = name_priv{ii};
	argpair{end+[1:2]} = {name, ob.(name)};
end

for ii=1:numel(name_hand)
	name = name_hand{ii};
	argpair{end+[1:2]} = {name, ob.(['handle_' name])};
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


function plotMask(mask)
colormap gray
imagesc(abs(mask));
end


text:
function patches = sampleIMAGES
% sampleIMAGES
% Returns 10000 patches for training

load '../data/IMAGES';

patchsize = 8;
numpatches = 10000;

% Initialize patches with zeros
patches = zeros(patchsize*patchsize, numpatches);
%% ---------- YOUR CODE HERE --------------------------------------
%  Instructions: Fill patches using data from IMAGES
%  IMAGES is a 3D array containing 10 images
%  For instance, IMAGES(21:30,21:30,1) is an image patch corresponding
%  to the pixels in the block (21,21) to (30,30) of Image 1

[ydim, xdim, num_images] = size(IMAGES);

for i = 1:numpatches
  img = randi(num_images);
  y_start = randi(ydim - patchsize + 1);
  x_start = randi(xdim - patchsize + 1);
  patch = IMAGES(y_start:y_start+patchsize-1, x_start:x_start+patchsize-1, img);
  patches(:,i) = patch(:);
end

%% ---------------------------------------------------------------
% For the autoencoder to work well we need to normalize the data
% Specifically, since the output of the network is bounded between [0,1]
% (due to the sigmoid activation function)
% It is important to make sure the input is also bounded between [0,1]
patches = normalizeData(patches);

end

%% ---------------------------------------------------------------
function patches = normalizeData(patches)

% Squash data to [0.1, 0.9] since we use sigmoid as the activation
% function in the output layer

% Remove DC
patches = bsxfun(@minus, patches, mean(patches));

% Truncate to +/-3 STD and scale to -1 to 1
pstd = 3 * std(patches(:));
patches = max(min(patches, pstd), -pstd) / pstd;

% Rescale to 0.1 to 0.9
patches = (patches + 1) * 0.4 + 0.1;

end

author: zellyn


text:
function test_stepsize_alg_demo()
% demonstration file for original stepsize algorithm.
%
% This file illustrates how to set user's own stepsize algorithm in case of linear
% regression problem. This demonstrates SGD and SVRG algorithms.
%
% This file is part of SGDLibrary.
%
% Created by H.Kasai on Sep. 25, 2017

    clc;
    clear;
    close all;

    %% generate synthetic data
    % set number of dimensions
    d = 10;
    % set number of samples
    n = 1000;
    % generate data
    data = logistic_regression_data_generator(n, d);

    %% define problem definitions
    problem = logistic_regression(data.x_train, data.y_train, data.x_test, data.y_test);

    %% perform algorithms SGD and SVRG
    options.w_init = data.w_init;
    options.step_init = 0.01;
    options.verbose = 2;

    options.step_alg = 'fix';
    [w_sgd_fix, info_sgd_fix] = sgd(problem, options);

    options.step_alg = 'decay';
    [w_sgd_decay, info_sgd_decay] = sgd(problem, options);

    options.step_alg = 'decay-2';
    [w_sgd_decay2, info_sgd_decay2] = sgd(problem, options);

    options.stepsizefun = @my_stepalg;  % set my_stepalg (user-defined stepsize algorithm)
    [w_sgd_my, info_sgd_my] = sgd(problem, options);

    %% display cost/optimality gap vs number of gradient evaluations
    display_graph('grad_calc_count','cost', {'SGD (fix)','SGD (decay)', 'SGD (decay-2)', 'SGD (My stepsize algorithm)'}, ...
            {w_sgd_fix, w_sgd_decay w_sgd_decay2, w_sgd_my}, {info_sgd_fix, info_sgd_decay, info_sgd_decay2, info_sgd_my});

end

%% define user-defined stepsize algorithm
function step = my_stepalg(iter, options)
    step = options.step_init / (10 + iter*0.5);
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function [eeg] = read_ns_eeg(filename, epoch)

% READ_NS_EEG read a NeuroScan 3.x or 4.x EEG File
%
% [eeg] = read_ns_eeg(filename, epoch)
%
%   filename     input Neuroscan .eeg file (version 3.x)
%   epoch        which epoch to read (default is all)
%
% The output data structure eeg has the fields:
%   eeg.data(..)    - epoch signal in uV (size: Nepoch x Nchan x Npnt)
% and
%   eeg.label       - electrode labels
%   eeg.nchan       - number of channels
%   eeg.npnt        - number of samplepoints in ERP waveform
%   eeg.time        - time for each sample
%   eeg.rate        - sample rate (Hz)
%   eeg.xmin        - prestimulus epoch start (e.g., -100 msec)
%   eeg.xmax        - poststimulus epoch end (e.g., 900 msec)
%   eeg.nsweeps     - number of accepted trials/sweeps

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

% read the neuroscan header
eeg = read_ns_hdr(filename);

% clear the variance part which is empty anyway
eeg = rmfield(eeg, 'variance');

% create a time axis
eeg.time = linspace(eeg.xmin, eeg.xmax, eeg.npnt);

% open the file and seek towards the place where the raw data is
fid = fopen_or_error(filename,'r','ieee-le');

% the default is to read all epochs
if nargin<2
  epoch = 1:eeg.nsweeps;
end

% determine whether it is 16 or 32 bit data
fseek(fid, 0, 'eof');
header_size = 900 + 75*eeg.nchan;
file_size   = ftell(fid);
sample_size = (file_size-header_size)/(eeg.nchan*eeg.npnt*eeg.nsweeps);
% note that the V4 format can have some extra information at the end of
% the file, causing the sample size to be slightly larger than 2 or 4
if floor(sample_size)==2
  epoch_size = eeg.nchan*eeg.npnt*2 + 13;
  datatype ='int16';
elseif floor(sample_size)==4
  datatype ='int32';
  epoch_size = eeg.nchan*eeg.npnt*4 + 13;
elseif floor(sample_size)>4
  % although this is not to be expected, it might be due to a very large "footer" compared to the data size
  % Olga Sysoeva reported that this extention to the datatype detection fixed it for her (see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=547)
  datatype ='int32';
  epoch_size = eeg.nchan*eeg.npnt*4 + 13;
end

% create empty storage for the data
data = zeros(length(epoch), eeg.nchan, eeg.npnt);
sweep = struct();

for i=1:length(epoch)
  fseek(fid, 900, 'bof');               % skip general header
  fseek(fid, 75*eeg.nchan, 'cof');          % skip channel headers
  status = fseek(fid, (epoch(i)-1)*epoch_size, 'cof');  % skip first epochs
  if status~=0
    ft_error('seek error while reading epoch data');
  end

  % fprintf('reading epoch %d at offset %d\n', epoch(i), ftell(fid));

  % read sweep header
  sweep(i).accept   = fread(fid, 1, 'uchar');
  sweep(i).type     = fread(fid, 1, 'ushort');
  sweep(i).correct  = fread(fid, 1, 'ushort');
  sweep(i).rt       = fread(fid, 1, 'float32');
  sweep(i).response = fread(fid, 1, 'ushort');
  sweep(i).reserved = fread(fid, 1, 'ushort');

  % read raw signal
  raw = fread(fid, eeg.nchan*eeg.npnt, datatype);
  if length(raw)~=eeg.nchan*eeg.npnt
    ft_error('fread error while reading epoch data');
  end
  data(i,:,:) = reshape(raw, [eeg.nchan, eeg.npnt]);

end

% convert raw signal to uV
for chan=1:eeg.nchan
  data(:,chan,:) = (data(:,chan,:) - eeg.baseline(chan)) * eeg.factor(chan);
end

% store the epoch information and data in the output structure
eeg.data = squeeze(data);
eeg.sweep = squeeze(sweep');

fclose(fid);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function im = vl_impattern(varargin)
% VL_IMPATTERN Generate an image from a stock pattern
%   IM=VLPATTERN(NAME) returns an instance of the specified
%   pattern. These stock patterns are useful for testing algoirthms.
%
%   All generated patterns are returned as an image of class
%   DOUBLE. Both gray-scale and colour images have range in [0,1].
%
%   VL_IMPATTERN() without arguments shows a gallery of the stock
%   patterns. The following patterns are supported:
%
%   Wedge::
%     The image of a wedge.
%
%   Cone::
%     The image of a cone.
%
%   SmoothChecker::
%     A checkerboard with Gaussian filtering on top. Use the
%     option-value pair 'sigma', SIGMA to specify the standard
%     deviation of the smoothing and the pair 'step', STEP to specfity
%     the checker size in pixels.
%
%   ThreeDotsSquare::
%     A pattern with three small dots and two squares.
%
%   UniformNoise::
%     Random i.i.d. noise.
%
%   Blobs:
%     Gaussian blobs of various sizes and anisotropies.
%
%   Blobs1:
%     Gaussian blobs of various orientations and anisotropies.
%
%   Blob:
%     One Gaussian blob. Use the option-value pairs 'sigma',
%     'orientation', and 'anisotropy' to specify the respective
%     parameters. 'sigma' is the scalar standard deviation of an
%     isotropic blob (the image domain is the rectangle
%     [-1,1]^2). 'orientation' is the clockwise rotation (as the Y
%     axis points downards). 'anisotropy' (>= 1) is the ratio of the
%     the largest over the smallest axis of the blob (the smallest
%     axis length is set by 'sigma'). Set 'cut' to TRUE to cut half
%     half of the blob.
%
%   A stock image::
%     Any of 'box', 'roofs1', 'roofs2', 'river1', 'river2', 'spotted'.
%
%   All pattern accept a SIZE parameter [WIDTH,HEIGHT]. For all but
%   the stock images, the default size is [128,128].

% Author: Andrea Vedaldi

% AUTORIGHTS

if nargin > 0
  pattern=varargin{1} ;
  varargin=varargin(2:end) ;
else
  pattern = 'gallery' ;
end

patterns = {'wedge','cone','smoothChecker','threeDotsSquare', ...
            'blob', 'blobs', 'blobs1', ...
            'box', 'roofs1', 'roofs2', 'river1', 'river2'} ;

% spooling
switch lower(pattern)
  case 'wedge', im = wedge(varargin) ;
  case 'cone', im = cone(varargin) ;
  case 'smoothchecker', im = smoothChecker(varargin) ;
  case 'threedotssquare', im = threeDotSquare(varargin) ;
  case 'uniformnoise', im = uniformNoise(varargin) ;
  case 'blob', im = blob(varargin) ;
  case 'blobs', im = blobs(varargin) ;
  case 'blobs1', im = blobs1(varargin) ;
  case {'box','roofs1','roofs2','river1','river2','spots'}
    im = stockImage(pattern, varargin) ;
  case 'gallery'
    clf ;
    num = numel(patterns) ;
    for p = 1:num
      vl_tightsubplot(num,p,'box','outer') ;
      imagesc(vl_impattern(patterns{p}),[0 1]) ;
      axis image off ;
      title(patterns{p}) ;
    end
    colormap gray ;
    return ;
  otherwise
    error('Unknown patter ''%s''.', pattern) ;
end

if nargout == 0
  clf ; imagesc(im) ; hold on ;
  colormap gray ; axis image off ;
  title(pattern) ;
  clear im ;
end

function [u,v,opts,args] = commonOpts(args)
opts.size = [128 128] ;
[opts,args] = vl_argparse(opts, args) ;
ur = linspace(-1,1,opts.size(2)) ;
vr = linspace(-1,1,opts.size(1)) ;
[u,v] = meshgrid(ur,vr);

function im = wedge(args)
[u,v,opts,args] = commonOpts(args) ;
im = abs(u) + abs(v) > (1/4) ;
im(v < 0) = 0 ;

function im = cone(args)
[u,v,opts,args] = commonOpts(args) ;
im = sqrt(u.^2+v.^2) ;
im = im / max(im(:)) ;

function im = smoothChecker(args)
opts.size = [128 128] ;
opts.step = 16 ;
opts.sigma = 2 ;
opts = vl_argparse(opts, args) ;
[u,v] = meshgrid(0:opts.size(1)-1, 0:opts.size(2)-1) ;
im = xor((mod(u,opts.step*2) < opts.step),...
         (mod(v,opts.step*2) < opts.step)) ;
im = double(im) ;
im = vl_imsmooth(im, opts.sigma) ;

function im = threeDotSquare(args)
[u,v,opts,args] = commonOpts(args) ;
im = ones(size(u)) ;
im(-2/3<u & u<2/3 & -2/3<v & v<2/3) = .75 ;
im(-1/3<u & u<1/3 & -1/3<v & v<1/3) = .50 ;
[drop,i] = min(abs(v(:,1))) ;
[drop,j1] = min(abs(u(1,:)-1/6)) ;
[drop,j2] = min(abs(u(1,:))) ;
[drop,j3] = min(abs(u(1,:)+1/6)) ;
im(i,j1) = 0 ;
im(i,j2) = 0 ;
im(i,j3) = 0 ;

function im = blobs(args)
[u,v,opts,args] = commonOpts(args) ;
im = zeros(size(u)) ;
num = 5 ;
square = 2 / num ;
sigma = square / 2 / 3 ;
scales = logspace(log10(0.5), log10(1), num) ;
skews = linspace(1,2,num) ;
for i=1:num
  for j=1:num
    cy = (i-1) * square + square/2 - 1;
    cx = (j-1) * square + square/2 - 1;
    A = sigma * diag([scales(i) scales(i)/skews(j)])  * [1 -1 ; 1 1] / sqrt(2)  ;
    C = inv(A'*A) ;
    x = u - cx ;
    y = v - cy ;
    im = im + exp(-0.5 *(x.*x*C(1,1) + y.*y*C(2,2) + 2*x.*y*C(1,2))) ;
  end
end
im = im / max(im(:)) ;

function im = blob(args)
[u,v,opts,args] = commonOpts(args) ;
opts.sigma = 0.15 ;
opts.anisotropy = .5 ;
opts.orientation = 2/3 * pi  ;
opts.cut = false ;
opts = vl_argparse(opts, args) ;
im = zeros(size(u)) ;
th = opts.orientation ;
R = [cos(th) -sin(th) ; sin(th) cos(th)] ;
A = opts.sigma * R * diag([opts.anisotropy 1]) ;
T = [0;0] ;
[x,y] = vl_waffine(inv(A),-inv(A)*T,u,v) ;
im = exp(-0.5 *(x.^2 + y.^2)) ;
if opts.cut
  im = im .* double(x > 0) ;
end

function im = blobs1(args)
[u,v,opts,args] = commonOpts(args) ;
opts.number = 5 ;
opts.sigma = [] ;
opts = vl_argparse(opts, args) ;
im = zeros(size(u)) ;
square = 2 / opts.number ;
num = opts.number ;
if isempty(opts.sigma)
  sigma = 1/6 * square ;
else
  sigma = opts.sigma * square ;
end
rotations = linspace(0,pi,num+1) ;
rotations(end) = [] ;
skews = linspace(1,2,num) ;
for i=1:num
  for j=1:num
    cy = (i-1) * square + square/2 - 1;
    cx = (j-1) * square + square/2 - 1;
    th = rotations(i) ;
    R = [cos(th) -sin(th); sin(th) cos(th)] ;
    A = sigma * R * diag([1 1/skews(j)]) ;
    C = inv(A*A') ;
    x = u - cx ;
    y = v - cy ;
    im = im + exp(-0.5 *(x.*x*C(1,1) + y.*y*C(2,2) + 2*x.*y*C(1,2))) ;
  end
end
im = im / max(im(:)) ;

function im = uniformNoise(args)
opts.size = [128 128] ;
opts.seed = 1 ;
opts = vl_argparse(opts, args) ;
state = vl_twister('state') ;
vl_twister('state',opts.seed) ;
im = vl_twister(opts.size([2 1])) ;
vl_twister('state',state) ;

function im = stockImage(pattern,args)
opts.size = [] ;
opts = vl_argparse(opts, args) ;
switch pattern
  case 'river1',  path='river1.jpg' ;
  case 'river2',  path='river2.jpg' ;
  case 'roofs1',  path='roofs1.jpg' ;
  case 'roofs2',  path='roofs2.jpg' ;
  case 'box',     path='box.pgm' ;
  case 'spots',   path='spots.jpg' ;
end
im = imread(fullfile(vl_root,'data',path)) ;
im = im2double(im) ;
if ~isempty(opts.size)
  im = imresize(im, opts.size) ;
  im = max(im,0) ;
  im = min(im,1) ;
end

author: Cloud-CV

repo: object-proposals


text:
function J = computeCostMulti(X, y, theta)
%COMPUTECOSTMULTI Compute cost for linear regression with multiple variables
%   J = COMPUTECOSTMULTI(X, y, theta) computes the cost of using theta as the
%   parameter for linear regression to fit the data points in X and y

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta
%               You should set J to the cost.

J = sum((X * theta .- y) .^ 2) / (2 * m);

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

end

author: sfvsfv

repo: Mathematical-modeling


text:
% OtQueryPrior Move to the previous record of the query.
%
% Usage:
%    OtQueryPrior(hquery)

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-04-11 Created (Johan Meijdam, VORtech)

function iret = OtQueryPrior(hquery)
   iret = Ot('query.prior', hquery);
end

author: Sable

repo: mcbench-benchmarks


text:
function systemsound(string)
% systemsound v 1.0.0
% Programmed by Drew Weymouth
%
% Plays a Microsoft Windows system sound specified by an input string.
%   string may be an exact filename of a .wav file in the C:\Windows\Media directory
%   or a string with the "Windows" prefix removed from the filename.
%   Input string is case insensitive, and need not include '.wav' extension.
% Examples:
%   systemsound('critical stop')
%   systemsound('exclamation')
% Please Note: 'notify' and 'windows notify' are different sounds
%   likewise, 'recycle' and 'windows recycle', 'ringin' and 'windows ringin',
%   and 'ding' and 'windows ding' are different.
% Tested with Windows 7, results may vary for older versions.

dir= 'C:\Windows\Media\';
string= lower(string);
try [y,fs]= wavread(sprintf('%s%s',dir,string));
catch
    try [y,fs]= wavread(sprintf('%swindows %s',dir,string));
    catch, error('Invalid selection.')
    end
end
sound(y,fs,16)

end


text:
function [repairedData, coefficients, grossErrors] = repair_corrupted_data(dictionary, corruptedData)
% repair_missing_data.m
%
%   Given a set of complete data vectors, and another set of vectors with
%   missing entries, use L1-minimization to repair the incomplete vectors.
%   This function required the CVX package for semidefinite programming.
%
%
% Inputs:
%   dictionary           - a matrix whose columns are data vectors that are
%                          used as overcomplete basis to represent other
%                          vectors. note that some or all of these columns
%                          may contain gross errors.
%   corruptedData        - a matrix whose columns are data vectors with
%                          some elements that may have gross errors. these
%                          vectors will be repaired using the dictionary.
%
% Outputs:
%   repairedData        - a matrix whose columns are repaired versions of
%                         the given incomplete data vectors
%   coefficients        - a matrix whose ith column is  a list of
%                         coefficients used to represent the ith vector in
%                         in corrupted data as a linear combination of
%                         vectors in the dictionary
%   grossErrors         - a matrix whose ij-th entry is true if the ij-th
%                         entry of corruptedData was corrupted by gross
%                         errors.
% Dependencies:
%   CVX package
%
% Nov. '07  Shankar Rao -- srrao@uiuc.edu

% Copyright 2007, University of Illinois. All rights reserved.

VERBOSE = false;
GROSS_ERROR_THRESHOLD = 0.5;

[dimensionCount, sampleCount] = size(dictionary);
corruptedSampleCount = size(corruptedData,2);

repairedData = corruptedData;
grossErrors = false(size(corruptedData));
coefficients = zeros(sampleCount+dimensionCount, corruptedSampleCount);
normalizedData = dictionary ./ repmat(sqrt(sum(dictionary.^2, 1)), dimensionCount, 1);
I = eye(dimensionCount);
X = [ normalizedData I];
for sampleIndex = 1:corruptedSampleCount
    y = corruptedData(:, sampleIndex);

    if VERBOSE,
        disp(sprintf('Repairing sample %d of %d', sampleIndex, corruptedSampleCount));
    end

    state = cvx_quiet(true);
    cvx_begin
        variable c(sampleCount+dimensionCount);
        minimize(norm(c,1));
            subject to
                X * c == y;
    cvx_end
    cvx_quiet(state);

    yhat = X(:, 1:sampleCount)*c(1:sampleCount);
    coefficients(:, sampleIndex) = c;
    errors = c(sampleCount+1:end);
    grossErrors(:, sampleIndex) = abs(errors) > GROSS_ERROR_THRESHOLD;
    repairedData(grossErrors(:, sampleIndex), sampleIndex) = yhat(grossErrors(:, sampleIndex));
end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
function all_XTr_pseudo = turn_training_simultaneous_data_into_pseudo_populations(all_XTr_simul, all_YTr, feature_to_siteID_mapping)

% A helper function takes simultaneously recorded training data and turns them into pseudo-populations.
%  This is useful for assessing whether training on shuffled data leads to similar performance (and decision boundaries)
%  as training with simultaneously recorded data (i.e., this allows one to compute I_diag as discussed by Averbeck, Latham and Pouget,
%  in Neural corelations, population coding, and computation, Nature Neuroscience May 2006).
%
%  Notes about the code:
%   1.  Care is taken so that if population vectors contain neuron's responses from multiple time periods,
%        when the data is shuffled, all data from the same trial (from a given site) will be permuted together.
%        This is important because mixing responses from the same trial (and site) into different population vectors could
%        lead to biased results since there are dependencies between time points in a trial.
%
%   2.  This code is slow.  It might be possible to write a faster version of the code (particularly
%        if one doesn't have totally different randomization for each label, etc).
%
%  Inputs:
%    all_XTr_simul{iTimePeriod}{iCV} = [num_features x num_points]  The simultaneously population training data
%    all_YTr The training labels  (can be a cell array for each time period or a vector having the same labels for all time periods)
%    feature_to_siteID_mapping{iTimePeriod} = a vector (length of num_features) that lists which site an a given feature
%      in all_XTr_simul came from (this is important so that all features from a given site are moved together).
%
%  Output:  all_XTr_pseudo, which is the simultaneous data converted into pseudo-populations

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

%     This code is part of the Neural Decoding Toolbox.
%     Copyright (C) 2011 by Ethan Meyers (emeyers@mit.edu)
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 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/>.

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

for iTimePeriod = 1:length(all_XTr_simul)
    for iCV = 1:length(all_XTr_simul{1})

        if iscell(all_YTr)
            curr_labels = all_YTr{iTimePeriod};
        else
            curr_labels = all_YTr;
        end
        unique_labels = unique(curr_labels);

        unique_site_numbers = unique(feature_to_siteID_mapping{iTimePeriod});

        all_XTr_pseudo{iTimePeriod}{iCV} = NaN .* ones(size(all_XTr_simul{iTimePeriod}{iCV}));

        for iLabel = unique_labels'

            curr_label_inds = find(curr_labels == iLabel);

            for iSite = unique_site_numbers

                curr_rand_label_inds = curr_label_inds(randperm(length(curr_label_inds)));

                site_feature_inds  = find(feature_to_siteID_mapping{iTimePeriod} == iSite);

                all_XTr_pseudo{iTimePeriod}{iCV}(site_feature_inds, curr_label_inds) = all_XTr_simul{iTimePeriod}{iCV}(site_feature_inds, curr_rand_label_inds);   %[num_neuron_features x num_trials]
            end
        end

    end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [PBarVals,dPBarValsdTheta,d2PBarValsdTheta2]=LegendreCos(theta,M,scalFactor)
%%LEGENDRECOS Evaluate the fully associated Legendre functions of
%           cos(theta) of degree n and order m for all n from 0 to M and
%           for each n, m goes from 0 to n. Such a function is typically
%           written as \bar{P}_{nm}(cos(theta)). Also evaluate
%           D{\bar{P}_{nm}(cos(theta))} and D2{\bar{P}_{nm}(cos(theta))},
%           where D{} is the first derivative with respect to theta and
%           D2{} is the second derivative with respect to theta. All of the
%           values can be scaled by a factor of scalFac, if desired, to
%           help prevent overflows with high degrees and orders.
%
%INPUTS: theta An angle in radians.
%       maxDeg The maximum degree and order of the output. This must be
%              >=0.
%   scalFactor A scale factor to help prevent overflow of the results. If
%              this parameter is omitted or an empty matrix is passed, the
%              default of 1 is used.
%
%OUTPUTS: PBarUVals An instance of the CountingClusterSet class such that
%                   PBarVals(n+1,m+1)=scalFac*\bar{P}_{nm}(cos(theta)).
%                   To extract all coefficients as a vector just call
%                   PBarUVals(:).
%   dPBarValsdTheta An instance of the CountingClusterSet class such that
%                   dPBarValsdTheta(n+1,m+1)=scalFac*D{\bar{P}_{nm}(cos(theta))}
% d2PBarValsdTheta2 An instance of the CountingClusterSet class such that
%                   d2PBarValsdTheta2(n+1,m+1)=scalFac*D2{\bar{P}_{nm}(cos(theta))}
%
%The definition of associated Legendre polynomial underlying this function
%is
% P_(nm)(x)=(1-x^2)^(m/2)*Dm{P(n,x)}
%where Dm{} represents the mth-order derivative and P(n,x) is the Legendre
%polynomial having the form
%P(n,x)=(1/2^n)*sum_{k=0}^n binomial(n,n)*(x-1)^(n-k)*(x+k)^k
%This definition lacks the Condon-Shortley phase. Adding the phase is the
%same as multiplying P_(nm)(x) by (-1)^m. The polynomials are fully
%normalized, meaning that each one is multiplied by
%N(n,m)=sqrt((2*n+1)*(2-KDelta(m))*factorial(n-m)/factorial(n+m))
%This is normalization 0 in the changeSpherHarmonicNorm function.
%
%Fully normalized associated Legendre polynomials often arise when
%computing spherical harmonic expansions. When simply evaluating an
%expansion, one will typically use fully normalized Legendre functions that
%have been divided by abs(sin(theta))^m as in [1]. However, when fitting
%coefficients for such an expansion, one will typically need the function
%values directly, which is what this function provides.
%
%This function just calls the function NALegendreCosRat and then multiplies
%out the denominators.
%
%EXAMPLE:
%One can see that this function returns the same values as the legendre
%function of cos(theta) when the proper normalization is selected. However,
%this function returns more values than legendre, as Matlab's legendre
%function just returns a vector for a single value of n.
% M=3;
% theta=-0.5;
% PBarVecMatLab=legendre(M,cos(theta));
% PBarVals=LegendreCos(theta,M);
% arePolyVals=true;
% origNormType=0;
% normType=8;
% PBarVals=changeSpherHarmonicNorm(PBarVals,origNormType,normType,arePolyVals);
% PBarVec=PBarVals(M+1,:);
% ratioCompare=PBarVecMatLab./PBarVec
%One will see that the ratio is essentially all ones.
%
%REFERENCES:
%[1] S. A. Holmes and W. E. Featherstone, "A unified approach to the
%    Clenshaw summation and the recursive computation of very high degree
%    and order normalised associated Legendre functions," Journal of
%    Geodesy, vol. 76, no. 5, pp. 279-299, May 2002.
%
%June 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(scalFactor))
    scalFactor=1;
end

u=abs(sin(theta));

if(nargin<2)
    PBarVals=NALegendreCosRat(theta,M,scalFactor);
    for n=1:M
        um=u;
        for m=1:n
            PBarVals(n+1,m+1)=PBarVals(n+1,m+1)*um;
            um=um*u;
        end
    end
elseif(nargout==2)
    [PBarVals,dPBarValsdTheta]=NALegendreCosRat(theta,M,scalFactor);
    for n=1:M
        um=u;
        for m=1:n
            PBarVals(n+1,m+1)=PBarVals(n+1,m+1)*um;
            dPBarValsdTheta(n+1,m+1)=dPBarValsdTheta(n+1,m+1)*um;
            um=um*u;
        end
    end
else%nargout==3
    [PBarVals,dPBarValsdTheta,d2PBarValsdTheta2]=NALegendreCosRat(theta,M,scalFactor);
    for n=1:M
        um=u;
        for m=1:n
            PBarVals(n+1,m+1)=PBarVals(n+1,m+1)*um;
            dPBarValsdTheta(n+1,m+1)=dPBarValsdTheta(n+1,m+1)*um;
            d2PBarValsdTheta2(n+1,m+1)=d2PBarValsdTheta2(n+1,m+1)*um;
            um=um*u;
        end
    end
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [ a, b ] = p53_lim ( )

%*****************************************************************************80
%
%% P53_LIM returns the integration limits for problem 53.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function r=geodistance(ci,cf,m)

%GEODISTANCE: Calculates the distance in meters between two points on earth surface.
%
% Usage:  r = geodistance( coordinates1 , coordinates2 , method ) ;
%
%	  Where coordinates1 = [longitude1,latitude1] defines the
%	  initial position and coordinates2 = [longitude2,latitude2]
%	  defines the final position.
%	  Coordinates values should be specified in decimal degrees.
%	  Method can be an integer between 1 and 23, default is m = 6.
%         Methods 1 and 2 are based on spherical trigonometry and a
%         spheroidal model for the earth, respectively.
%	  Methods 3 to 24 use Vincenty's formulae, based on ellipsoid
%         parameters.
%         Here it follows the correspondence between m and the type of
%         ellipsoid:
%
%         m =  3 -> ANS ,        m =  4 -> GRS80,    m = 5 -> WGS72,
%         m =  6 -> WGS84,       m =  7 -> NSWC-9Z2,
%         m =  8 -> Clarke 1866, m =  9 -> Clarke 1880,
%         m = 10 -> Airy 1830,
%         m = 11 -> Bessel 1841 (Ethiopia,Indonesia,Japan,Korea),
%         m = 12 -> Bessel 1841 (Namibia),
%         m = 13 -> Sabah and Sarawak (Everest,Brunei,E.Malaysia),
%         m = 14 -> India 1830, m = 15 -> India 1956,
%         m = 16 -> W. Malaysia and Singapore 1948,
%         m = 17 -> W. Malaysia 1969,
%         m = 18 -> Helmert 1906, m = 19 -> Helmert 1960,
%         m = 20 -> Hayford International 1924,
%         m = 21 -> Hough 1960, m = 22 -> Krassovsky 1940,
%         m = 23 -> Modified Fischer 1960,
%         m = 24 -> South American 1969.
%
%	  Important notes:
%
%	 1)South latitudes are negative.
%	 2)East longitudes are positive.
%	 3)Great circle distance is the shortest distance between two points
%          on a sphere. This coincides with the circumference of a circle which
%          passes through both points and the centre of the sphere.
%	 4)Geodesic distance is the shortest distance between two points on a spheroid.
%	 5)Normal section distance is formed by a plane on a spheroid containing a
%          point at one end of the line and the normal of the point at the other end.
%          For all practical purposes, the difference between a normal section and a
%          geodesic distance is insignificant.
%	 6)The method m=2 assumes a spheroidal model for the earth with an average
%          radius of 6364.963 km. It has been derived for use within Australia.
%          The formula is estimated to have an accuracy of about 200 metres over 50 km,
%          but may deteriorate with longer distances.
%          However, it is not symmetric when the points are exchanged.
%
%  Examples: A = [150 -30]; B = [150 -31]; L = [151 -80];
%            [geodistance(A,B,1) geodistance(A,B,2) geodistance(A,B,3)]
%            [geodistance(A,L,1) geodistance(A,L,2) geodistance(A,L,3)]
%            geodistance([0 0],[2 3])
%            geodistance([2 3],[0 0])
%            geodistance([0 0],[2 3],1)
%            geodistance([2 3],[0 0],1)
%            geodistance([0 0],[2 3],2)
%            geodistance([2 3],[0 0],2)
%            for m = 1:24
%            r(m) = geodistance([150 -30],[151 -80],m);
%            end
%            plot([1:m],r), box on, grid on
%
%***************************************************************************************
% Second version: 07/11/2007
% Third  version: 03/08/2010
%
% Contact: orodrig@ualg.pt
%
% Any suggestions to improve the performance of this
% code will be greatly appreciated.
%
% Reference: Geodetic Calculations Methods
%            Geoscience Australia
%            (http://www.ga.gov.au/geodesy/calcs/)
%
%***************************************************************************************
% Changed a little bit by Mao,
% So, it works in vector computaiton now.
%
% Contact: argansos@hotmail.com
% Ganquan Mao, 15/09/2011
%
%***************************************************************************************

if size(ci,2)~=2 || size(cf,2)~=2
    error('ci cf must be a n*2 matrix!')
end

if size(ci,1)~=size(cf,1)
    error('ci cf must have same size!')
end

r=[];%for convergence

if nargin==2
    m=6;
end

lambda1=ci(:,1)*pi/180;
   phi1=ci(:,2)*pi/180;

lambda2=cf(:,1)*pi/180;
   phi2=cf(:,2)*pi/180;

L=lambda2-lambda1;

if m==1%great circle distance, based on spherical trigonometry
    r=180*1.852*60*acos(sin(phi1).*sin(phi2)...
      +cos(phi1).*cos(phi2).*cos(lambda2-lambda1))./pi;
    r=1000*abs(r);

elseif m==2%spheroidal model for the earth
    term1=111.08956*(ci(:,2)-cf(:,2)+0.000001);
    term2=cos(phi1+((phi2-phi1)/2));
    term3=(cf(:,1)-ci(:,1)+0.000001)/(cf(:,2)-ci(:,2)+0.000001);
    r=1000*abs(term1./cos(atan(term2.*term3)));

else%apply Vincenty's formulae (as long as the points are not coincident)
    alla=[0,0,6378160,6378137,6378135,6378137,6378145,6378206.4,...
          6378249.145,6377563.396,6377397.155,6377483.865,6377298.556,...
          6377276.345,6377301.243,6377304.063,6377295.664,6378200,...
          6378270,6378388,6378270,6378245,6378155,6378160];

    allf=[0,0,1/298.25,1/298.257222101,1/298.26,1/298.257223563,...
          1/298.25,1/294.9786982,1/293.465,1/299.3249646,1/299.1528128,...
          1/299.1528128,1/300.8017,1/300.8017,1/300.8017,1/300.8017,...
          1/300.8017,1/298.3,1/297,1/297,1/297,1/298.3,1/298.3,1/298.25];

    a=alla(m);
    f=allf(m);

    b=a*(1-f);

    axa=a^2;
    bxb=b^2;

    U1=atan((1-f)*tan(phi1));
    U2=atan((1-f)*tan(phi2));

    lambda=L;
    lambda_old=sqrt(-1);%there is no way a complex number is going to coincide with a real number!

    ntrials=0;%just in case...

    while any(abs(lambda-lambda_old)>1e-9)

        ntrials=ntrials+1;

        lambda_old=lambda;
        sin_sigma=sqrt((cos(U2).*sin(lambda)).^2+(cos(U1).*sin(U2)...
                  -sin(U1).*cos(U2).*cos(lambda)).^2);
        cos_sigma=sin(U1).*sin(U2)+cos(U1).*cos(U2).*cos(lambda);
        sigma=atan2(sin_sigma,cos_sigma);
        sin_alpha=cos(U1).*cos(U2).*sin(lambda)./sin_sigma;
        cos2_alpha=1-sin_alpha.^2;
        cos_2sigmam=cos_sigma-2*sin(U1).*sin(U2)./cos2_alpha;

        C=(f/16)*cos2_alpha.*(4+f*(4-3*cos2_alpha));

        lambda=L+f*(1-C).*sin_alpha.*(sigma+C.*sin_sigma.*(cos_2sigmam...
               +C.*cos_sigma.*(-1+2*(cos_2sigmam).^2)));

       % stop the function if convergence is not achieved
        if ntrials>1000
            disp('Convergence failure...')
            return
        end

    end

    % convergence achieved? get the distance
    uxu=cos2_alpha*(axa-bxb)/bxb;
    A=1+(uxu/16384).*(4096+uxu.*(-768+uxu.*(320-175*uxu)));
    B=(uxu/1024).*(256+uxu.*(-128+uxu.*(74-47*uxu)));
    delta_sigma=B.*sin_sigma.*(cos_2sigmam+(B/4).*(cos_sigma.*(-1+...
                2*cos_2sigmam.^2)-(B/6).*cos_2sigmam.*(-3+...
                4*sin_sigma.^2).*(-3+4*cos_2sigmam.^2)));
    r=b*A.*(sigma-delta_sigma);

end

r(isnan(r))=0;

end

author: goGPS-Project


text:
function [Population,succ1,succ2,Fitness] = First_Stage_EnvironmentalSelection(Population,N,isOrigin)
% The environmental selection of SPEA2

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

% This function is written by Kangjia Qiao

    %% Calculate the fitness of each solution
    if isOrigin==1
        Fitness = CalFitness(Population.objs,Population.cons);
    else
        Fitness = CalFitness(Population.objs);
    end

    %% Environmental selection
    Next = Fitness < 1;
    if sum(Next) < N
        [~,Rank] = sort(Fitness);
        Next(Rank(1:N)) = true;
    elseif sum(Next) > N
        Del  = Truncation(Population(Next).objs,sum(Next)-N);
        Temp = find(Next);
        Next(Temp(Del)) = false;
    end

    % Population for next generation
    Population = Population(Next);
    Fitness    = Fitness(Next);
    % Sort the population
    [Fitness,rank] = sort(Fitness);
    Population = Population(rank);

    %% calculate success rate1
    off_index = Next(1+N:2*N);
    succ1 = zeros(1,length(off_index));
    for j = 1:length(off_index)
        if off_index(j) == 1
            succ1(j)=1;
        end
    end
    %% calculate success rate2
    off_index = Next(1+2*N:end);
    succ2 = zeros(1,length(off_index));
    for j = 1:length(off_index)
        if off_index(j) == 1
            succ2(j)=1;
        end
    end
end

function Del = Truncation(PopObj,K)
% Select part of the solutions by truncation

    %% Truncation
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Del = false(1,size(PopObj,1));
    while sum(Del) < K
        Remain   = find(~Del);
        Temp     = sort(Distance(Remain,Remain),2);
        [~,Rank] = sortrows(Temp);
        Del(Remain(Rank(1))) = true;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% Script to make a mrSESSION directory from analyze-format data
% You need to fill in many of the parameters below.
% Note that you need to have converted everything (Inplanes and TSeries) to analyze format before
% you run this script.
% Assume we run this in the destination dir.

clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set up directory info

anatomyInplaneDir=26; % This is the location of the IRRARE T1 weighted data

% Where are the source files for the data?
% sourceFileBaseDir='/raid/MRI/data/wade/Nikos_Data/Monkey/Spring2003/original/C01.j51/converted/'
% sourceFileBaseDir='//biac1/wandell/data/monkey/retinal-lesion/A01/Ky1/converted/'
 sourceFileBaseDir='G:/data/monkey/retinal-lesion/convertedScans/A01/ky1/converted/'

sourceSubDirArray=[16]; %locations of the functional scans, ex. [1:4];
nFunctionals=length(sourceSubDirArray);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Enter mrSESSION info - you will need to edit many of these parameters
% These steps replace mrInitRet.

MLR_ScanOffset=0;  % mrLoadRet needs scans numbered 1:n; we need to convert from scans numbered 16:25, for example
                   % the offsest lets us start the mrLoadRet numbering at a
                   % different number other than 1, when adding scans to a
                   % previously created directory.

M.mrLoadRetVersion=3.0100;
M.sessionCode='A01_Ky1_Ret'
M.description='MPI: Ret scans A01, session Ky1'
M.subject='A01'
M.examNum='A01_Ky1';

% INPLANES
I.FOV=128; % in mm
I.fullSize=[256 256]; % pixels
I.voxelSize=[0.5 0.5 2] % in mm
I.spacing=0
I.nSlices=17
I.examNum=M.examNum
I.crop=[0 0;255 255]; % top left, bottom right coordinates. In this case, no crop.
I.cropSize=[256 256];

M.inplanes=I;

% FUNCTIONALS
F.PfileName=0; % To be filled in later when we loop over the functionals
F.totalFrames=132;
F.junkFirstFrames=0;
F.nFrames=F.totalFrames-F.junkFirstFrames;
F.slices=1:I.nSlices;
F.fullSize=[256 256]; % Resolution, in pixels    %%128?
F.cropSize=[256 256]; % no need to crop image    %%128?
F.crop=[0 0;127 127];
F.voxelSize=[I.FOV./F.fullSize,  I.voxelSize(3)];
F.effectiveResolution=F.voxelSize;
F.framePeriod=6.5;  % The length of 1 TR (TR = 'effective frame duration' = tAcq*nshots) times the number of interleaves
                    % for 15 slices, TR = 6; for 17 slices, TR = 6.5

% RECONPARAMS
R.nframes=F.totalFrames;
R.nshots=8;         %check excel file for this
R.FOV=I.FOV;
R.sliceThickness=I.voxelSize(3);
R.skip=0
R.TR=F.framePeriod;

R.tAcq=806;        %tAcq is the actual TR. Multiply this by the number of
                    %shots to get the effective frame duration. This is 8
                    %for 17 slice monkeys - check excel file.

R.slquant=I.nSlices;

R.equivMatSize=256;  %%128?
R.imgsize=128;  %%128?
F.reconParams=R;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now we loop through all the n directories listed in sourceSubDirArray
% We will create a 1xn functionals structure in mrSESSION by copying F a lot. The only thing we'll change will
% be F.PfileName We also create dataTYPES at this point. dataTYPES will be a single structure with some
% 1xn sub-structures.
%            scanParams: [1x14 struct]
%            blockedAnalysisParams: [1x14 struct]
%            eventAnalysisParams: [1x14 struct]
% eventAnalysisParams are all 0
% Scan params and blockedAnalysisParams are to be set in the loop by copying in the structures we define next.....

%SCAN PARAMETERS
S.annotation='retinotopy';
S.nFrames=F.nFrames;
S.framePeriod=F.framePeriod
S.slices=F.slices;
S.cropSize=F.cropSize;

%BLOCKED ANALYSIS PARAMETERS
B.blockedAnalysis=1;
B.detrend=1;      % -1 = linear detrend; 1 = high pass filter; varied opinions on which is better.
B.inhomoCorrect=1;
B.nCycles=11; % For the color expts=4. For the retinotopies, this is 12??? For J02 this is 11;

B.temporalNormalization=0; % This is essential for the MPG scans that were converted brucker to analyze with the BUG.
                           % For all good/normal scans set this to 0.

E.eventAnalysis=0;

nScans=length(sourceSubDirArray);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make more Dirs. Be sure to run makemrLoadRet3Session in the destination
% directory.

disp('Generating sub directories');

a=mkdir('Inplane'); % Holds data in 'slice' format
a=mkdir('Gray'); % Holds data in 3D restricted to the cortical gray matter
a=mkdir('Volume'); % Barely used. Holds data in 3D
a=mkdir('Raw'); % Original FID files, IRRARE files etc.
a=mkdir('Raw/Anatomy');
a=mkdir('Raw/Anatomy/Inplane');
a=mkdir('Inplane/ROIs'); % Regions of interest
a=mkdir('Inplane/Original');
a=mkdir('Inplane/Original/TSeries');
a=mkdir('Gray/ROIs');

tSerDir=['Inplane/Original/TSeries/'];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert the inplanes and generate the anat.mat file.

inpIndex=sprintf('%02d',anatomyInplaneDir);

rotateInplanes=1; % Number of 90degree ccw rotations to apply to everything to make it look nice.
                    % # rotations = rotateInplanes * 90

% anatFile=fullfile(sourceFileBaseDir,int2str(anatomyInplaneDir),[M.subject,'_',inpIndex])
% anat=analyze2mrLoadRetInplanes(fullfile(sourceFileBaseDir,int2str(anatomyInplaneDir),[M.subject,'_',inpIndex]),'./');
% anatFile=fullfile(sourceFileBaseDir,int2str(anatomyInplaneDir))
% anat=analyze2mrLoadRetInplanes(anatFile,'./');
anat=analyze2mrLoadRetInplanes(fullfile(sourceFileBaseDir,int2str(anatomyInplaneDir),'2dseq_001'),'./Raw/Anatomy/Inplane/');

[x,y,nSlices]=size(anat);

% Run through the anat matrix and do rotations
for t=1:nSlices
    anat(:,:,t)=rot90(anat(:,:,t),rotateInplanes);
end

save('Inplane/anat.mat','anat');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create dataTYPES structure

dataTYPES.name='Original';

for t=1:nFunctionals
    F.PfileName=['P_MPI_',int2str(sourceSubDirArray(t))];
    M.functionals(t)=F;

    % Also fill in the dataTYPES structure.
    dataTYPES.blockedAnalysisParams(t)=B;
    dataTYPES.eventAnalysisParams(t)=E;
    dataTYPES.scanParams(t)=S;
    dataTYPES.name='Original';

end

global mrSESSION;
global dataTYPES;
global HOMEDIR;

mrSESSION=M;

saveSession;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now we're ready to import and convert the functional data

disp ('Saved mrSESSION: Now converting tSeries data...');

for thisScan=1:nFunctionals

    % make the directory and name of the tSeries files to be created here
    thisTSerDir=fullfile(tSerDir,['Scan',int2str(thisScan+MLR_ScanOffset)])
    a=mkdir(thisTSerDir);

    % input path and name of input analyze files (2dseq_)
    thisFunc=int2str(sourceSubDirArray(thisScan));
    thisFunc02Padded=sprintf('%02d',sourceSubDirArray(thisScan));
    % sourceRoot=fullfile(sourceFileBaseDir,thisFunc,[M.subject,'_',thisFunc02Padded,'_EPI_']);
    sourceRoot=fullfile(sourceFileBaseDir,thisFunc,'2dseq_');

    % imStack=analyze2mrLoadRetTSeries(inFileRoot,outFileRoot,nVols,firstVolIndex,doRotate,scaleFact,flipudFlag,fliplrFlag)
    % where doRotate can be any integer. The images are rotated anticlockwise by doRotate*90degrees.
    % scaleFact is a 2*1 scale factor for scaling the x and y dimensions independently.
    % if flipudFlag or fliprlFlag is set, the images are flipped up/down or left/right after rotation
    imStack=analyze2mrLoadRet3TSeries(sourceRoot,thisTSerDir,F.nFrames,F.junkFirstFrames+1,1,1,0,0);

    fprintf('\nDone scan %d (directory %d)\n',thisScan,sourceSubDirArray(thisScan));

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [hmg, HMG_u, HMG_s, HMG_k, HMG_c] = invPinHoleHmg(u,s,k,c)

% INVPINHOLEHMG Retro-project anchored homogeneous point AHP.
%   HMG = INVPINHOLEHMG(U,S) gives the retroprojected anchored homogeneous
%   point (HMG) of a pixel U at inverse-depth S, from a canonical pin-hole
%   camera, that is, with calibration parameters
%     u0 = 0, v0 = 0, au = 1, av = 1
%   It uses reference frames {RDF,RD} (right-down-front for the 3D world
%   points and right-down for the pixel), according to this scheme:
%
%         / z (forward)
%        /
%       +------- x                 +------- u
%       |                          |
%       |      3D : P=[x;y;z]      |     image : U=[u;v]
%       | y                        | v
%
%   HMG = INVPINHOLEHMG(U,S,K) allows the introduction of the camera's
%   calibration parameters:
%     K = [u0 v0 au av]'
%
%   HMG = INVPINHOLEHMG(U,S,K,C) allows the introduction of the camera's radial
%   distortion correction parameters:
%     C = [c2 c4 c6 ...]'
%   so that the new pixel is corrected following the distortion equation:
%     U = U_D * (1 + K2*R^2 + K4*R^4 + ...)
%   with R^2 = sum(U_D.^2), being U_D the distorted pixel in the image
%   plane for a camera with unit focal length.
%
%   If U is a pixels matrix, INVPINHOLEHMG(U,...) returns a HMGS matrix HMG,
%   with these matrices defined as
%     U   = [U1 ... Un];       Ui   = [ui;vi]
%     HMG = [HMG1 ... HMGn];   HMGi = [xi,yi,zi,ri]
%   where xi, yi, zi are the non-homogeneous parts of HMG, defining the
%   optical ray, and ri is the inverse of the distance (wrongly named
%   "inverse depth")
%
%   [HMG,HMG_u,HMG_s,HMG_k,HMG_c] returns the Jacobians of HMG wrt U, S, K and C. It
%   only works for single pixels U=[u;v], and for distortion correction
%   vectors C of up to 3 parameters C=[c2;c4;c6]. See UNDISTORT for
%   information on longer distortion vectors.
%
%   See also RETRO, UNDISTORT, DEPIXELLISE, PINHOLEHMG.

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

if nargout == 1 % only point

    switch nargin
        case 2
            v = retro(u,1);
        case 3
            v = retro(depixellise(u,k),1);
        case 4
            v = retro(undistort(depixellise(u,k),c),1);
    end
    n = normvec(v);
    hmg  = [n;s];

else % Jacobians

    if size(u,2) > 1
        error('Jacobians not available for multiple pixels')
    else

        switch nargin
            case 2
                [v, V_u] = invPinHole(u,1);
                [n,N_v]  = normvec(v,1);
                hmg      = [n;s];
                HMG_v    = [N_v;zeros(1,3)];
                HMG_s    = [0;0;0;1];
                HMG_u    = HMG_v*V_u;

            case 3
                [v, V_u, V_1, V_k] = invPinHole(u,1,k);
                [n,N_v]  = normvec(v,1);
                hmg      = [n;s];
                HMG_v    = [N_v;zeros(1,3)];
                HMG_s    = [0;0;0;1];
                HMG_u    = HMG_v*V_u;
                HMG_k    = HMG_v*V_k;

            case 4
                [v, V_u, V_1, V_k, V_c] = invPinHole(u,1,k,c);
                [n,N_v]  = normvec(v,1);
                hmg      = [n;s];
                HMG_v    = [N_v;zeros(1,3)];
                HMG_s    = [0;0;0;1];
                HMG_u    = HMG_v*V_u;
                HMG_k    = HMG_v*V_k;
                HMG_c    = HMG_v*V_c;

        end
    end

end

return

%% jacobians
syms u v s u0 v0 au av c2 c4 c6 real
U=[u;v];
k=[u0;v0;au;av];
c=[c2;c4;c6];

% [hmg,HMG_u,HMG_s] = invPinHoleHmg(U,s);
% [hmg,HMG_u,HMG_s,HMG_k] = invPinHoleHmg(U,s,k);
[hmg,HMG_u,HMG_s,HMG_k,HMG_c] = invPinHoleHmg(U,s,k,c);

simplify(HMG_u - jacobian(hmg,U))
simplify(HMG_s - jacobian(hmg,s))
simplify(HMG_k - jacobian(hmg,k))
% simplify(HMG_c - jacobian(hmg,c))

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


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

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

% Last Modified by GUIDE v2.5 20-Feb-2013 10:37:02

% Copyright 2013 The MathWorks, Inc.

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

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

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

% Initialize the GUI layout
sketchplot(gcf);
set(get(handles.plotstreamline,'XLabel'),'String','x (m)')
set(get(handles.plotstreamline,'YLabel'),'String','y (m)')
set(get(handles.plotstreamline,'Title'),'String','Normalized Stream Function \Psi')
set(get(handles.plotvector,'XLabel'),'String','x (m)')
set(get(handles.plotvector,'YLabel'),'String','y (m)')
set(get(handles.plotvector,'Title'),'String','Velocity Profile')
% Create the help button on toolbar
[X, map] = imread(fullfile(...
    matlabroot,'toolbox','matlab','icons','csh_icon.gif'));
icon = ind2rgb(X,map);
uipushtool(handles.uitoolbar1,'CData',icon,...
    'TooltipString','Help',...
    'ClickedCallback',@AppHelp);

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

% Update handles structure
guidata(hObject, handles);

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

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

% Get default command line output from handles structure
varargout{1} = handles.output;

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

function handles = getpara(handles)
% Get user input parameters from the GUI
handles.para.D1 = str2double(get(handles.D1,'String'));
handles.para.D2 = str2double(get(handles.D2,'String'));
handles.para.L1 = str2double(get(handles.L1,'String'));
handles.para.L2 = str2double(get(handles.L2,'String'));
handles.para.U1 = str2double(get(handles.U1,'String'));
handles.para.nx = str2double(get(handles.nx,'String'));
handles.para.ny = str2double(get(handles.ny,'String'));
handles.para.nmax = str2double(get(handles.nmax,'String'));
cla(handles.plotstreamline),cla(handles.plotvector),
axis([handles.plotstreamline handles.plotvector], [0 handles.para.L1+handles.para.L2 0 max(handles.para.D1,handles.para.D2)]),
% Draw the channel shape according to input
if handles.para.D1 < handles.para.D2 % converging channel
    cornerposition = [0 handles.para.D1 handles.para.L1 handles.para.D2-handles.para.D1];
elseif handles.para.D1 > handles.para.D2 % diverging channel
    cornerposition = [handles.para.L1 handles.para.D2 handles.para.L2 handles.para.D1-handles.para.D2];
else
    cornerposition = [0 handles.para.D1 1 1]; % straight channel
end
rectangle('Parent',handles.plotstreamline,'Position',cornerposition,'FaceColor',[0.94 0.94 0.94])
rectangle('Parent',handles.plotvector,'Position',cornerposition,'FaceColor',[0.94 0.94 0.94])

function inviscidchannel(handles)
%% Initialization
% get input parameters
handles = getpara(handles);
% store input parameters with shorter names
D1 = handles.para.D1;
D2 = handles.para.D2;
L1 = handles.para.L1;
L2 = handles.para.L2;
U1 = handles.para.U1;
nx = handles.para.nx;
ny = handles.para.ny;
nmax = handles.para.nmax;
% compute grid length
dx = (L1+L2)/(nx-1);
dy = max(D1,D2)/(ny-1);
% nondimensionalization
dX = dx/D1;
dY = dy/D1;
PSI = ones(ny,nx);
%% Boundary conditions
PSI(1,:) = 0; % bottom
PSI(floor(D1/dy+1),1:floor(L1/dx+1)) = 1; % left part of top
PSI(floor(D2/dy+1),floor(L1/dx+1):nx) = 1; % right part of top
for j = 1:floor(D1/dy+1)
    PSI(j,1) = (j-1)*dY; % inlet
end
for j = 1:floor(D2/dy+1)
    PSI(j,nx) = D1/D2*(j-1)*dY; % outlet
end
%% Iteration to solve for stream function PSI
err = 1e-6; % convergence criteria
n = 0;
while n < nmax
    n = n + 1;
    tempPSI = PSI;
    % left part of the channel
    for i = 2:floor(L1/dx+1)
        for j = 2:(floor(D1/dy+1)-1)
            PSI(j,i) = 1/(2/dX^2+2/dY^2)*((tempPSI(j,i+1)+PSI(j,i-1))/dX^2+(tempPSI(j+1,i)+PSI(j-1,i))/dY^2);
        end
    end
    % right part of the channel
    for i = floor(L1/dx+1):nx-1
        for j = 2:(floor(D2/dy+1)-1)
            PSI(j,i) = 1/(2/dX^2+2/dY^2)*((tempPSI(j,i+1)+PSI(j,i-1))/dX^2+(tempPSI(j+1,i)+PSI(j-1,i))/dY^2);
        end
    end
    % checking for convergence
    if max(max(abs(PSI-tempPSI))) <= err
        break
    end
end
%% Prompt convergence result in a message box
if n < nmax
    helpdlg({'Converged.';sprintf('Total number of iterations: %g\n',n)},...
        'Converged');
else
    warndlg({'NOT Converged!';sprintf('Change n_x, n_y, n_max or convergence criteria.\n')},...
        'NOT Converged');
end
%% Solve for velocity u and v from stream function
u = zeros(ny,nx);
v = zeros(ny,nx);
psi = PSI*U1*D1; % dimentionalization
for j = 2:ny-1
    u(j,:) = (psi(j+1,:)-psi(j-1,:))/2/dy;
end
for i = 2:nx-1
    v(:,i) = -(psi(:,i+1)-psi(:,i-1))/2/dx;
end
% velocity at boundaries (inviscid, so slip occurs)
u(1,:) = u(2,:);
u(ny,:) = u(ny-1,:);
v(:,1) = v(:,2);
v(:,nx) = v(:,nx-1);
%% Plot results
[X, Y] = meshgrid(0:dx:(L1+L2),0:dy:max(D1,D2));
% contour plot of stream function
contourf(handles.plotstreamline,X,Y,PSI),colormap(handles.plotstreamline,'jet')
% vector plot of velocity u and v
quiver(handles.plotvector,X,Y,u,v),
axis([handles.plotstreamline handles.plotvector],[0 L1+L2 0 max(D1,D2)]),
% plot the shape of the channel based on input parameters
if D1 < D2
    cornerposition = [0 D1 L1 D2-D1];
elseif D1 > D2
    cornerposition = [L1 D2 L2 D1-D2];
else
    cornerposition = [0 D1 1 1];
end
rectangle('Parent',handles.plotstreamline,'Position',cornerposition,'FaceColor',[0.94 0.94 0.94])
rectangle('Parent',handles.plotvector,'Position',cornerposition,'FaceColor',[0.94 0.94 0.94])
set(get(handles.plotstreamline,'XLabel'),'String','x (m)')
set(get(handles.plotstreamline,'YLabel'),'String','y (m)')
set(get(handles.plotstreamline,'Title'),'String','Normalized Stream Function \Psi')
set(get(handles.plotvector,'XLabel'),'String','x (m)')
set(get(handles.plotvector,'YLabel'),'String','y (m)')
set(get(handles.plotvector,'Title'),'String','Velocity Profile')

function sketchplot(h)
%Sketch the plot of the channel to show geometric parameters
% and governing equations. You can learn how to draw sketches
% programmatically and write LaTeX. Don't worry about it otherwise,
% since it has nothing to do with computation.

% Create doublearrow
d1 = annotation(h,'doublearrow',[0.06 0.06],...
    [0.6 0.8]);

% Create doublearrow
l1 = annotation(h,'doublearrow',[0.06 0.15],...
    [0.6 0.6]);

% Create doublearrow
l2 = annotation(h,'doublearrow',[0.15 0.3],...
    [0.6 0.6]);

% Create doublearrow
d2 = annotation(h,'doublearrow',[0.3 0.3],...
    [0.7 0.6]);

% Create arrow U
annotation(h,'arrow',[0.02 0.06],...
    [0.7 0.7],...
    'HeadStyle','Plain',...
    'HeadWidth',5);

% Create textbox U
annotation(h,'textbox',...
    [0.02 0.67 0.03 0.07],...
    'String',{'U'},...
    'FitBoxToText','off',...
    'LineStyle','none');

% Create textbox
annotation(h,'textbox',...
    [0.06 0.65 0.03 0.07],...
    'String',{'D1'},...
    'FitBoxToText','off',...
    'LineStyle','none');

% Create textbox
annotation(h,'textbox',...
    [0.1 0.58 0.03 0.07],...
    'String',{'L1'},...
    'FitBoxToText','off',...
    'LineStyle','none');

% Create textbox
annotation(h,'textbox',...
    [0.22 0.58 0.03 0.07],...
    'String',{'L2'},...
    'FitBoxToText','off',...
    'LineStyle','none');

% Create textbox
annotation(h,'textbox',...
    [0.25 0.61 0.03 0.07],...
    'String',{'D2'},...
    'FitBoxToText','off',...
    'LineStyle','none');

% Create channel
annotation(h,'line',...
    [0.05 0.15],[0.8 0.8],...
    'LineWidth',2);

% Create channel
annotation(h,'line',...
    [0.15 0.31],[0.7 0.7],...
    'LineWidth',2);

% Create channel
annotation(h,'line',...
    [0.15 0.15],[0.7 0.8],...
    'LineWidth',2);

% Create channel
annotation(h,'line',...
    [0.05 0.31],[0.58 0.58],...
    'LineWidth',2);

% Create textbox equations
annotation(h,'textbox',[0.32 0.73 0.03 0.07],...
    'String',{'\partial^2\psi/\partialx^2 + \partial^2\psi/\partialy^2 = 0','','u = \partial\psi/\partialy','', 'v = -\partial\psi/\partialx'},...
    'FitBoxToText','on',...
    'LineStyle','none');

function AppHelp(varargin)
% Create help message
dlgname = 'About Inviscid Channel App';
txt = {'Inviscid flow in a channel (Laplace equation of stream function)';
    'solved with explicit method';
    '';
    'Start - Start the simulation';
    'Solution Code - Open a pre-written script to show the solution code';
    '';
    '* Try D1 > D2, D1 < D2 or D1 = D2.';
    '* Use zoom, pan and data curser in the toolbar to interact with the';
    '  plots as usual';
    '* If result does not  converge, try:';
    '       a) increasing n_max (total number of iterations)';
    '       b) decreasing n_x, n_y (number of grids in each direction)';
    '';
    'Copyright 2013 The MathWorks, Inc.'};
helpdlg(txt,dlgname);

function D1_Callback(hObject, eventdata, handles)
% hObject    handle to D1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
getpara(handles);
% Hints: get(hObject,'String') returns contents of D1 as text
%        str2double(get(hObject,'String')) returns contents of D1 as a double

function D2_Callback(hObject, eventdata, handles)
% hObject    handle to D2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
getpara(handles);
% Hints: get(hObject,'String') returns contents of D2 as text
%        str2double(get(hObject,'String')) returns contents of D2 as a double

function L1_Callback(hObject, eventdata, handles)
% hObject    handle to L1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
getpara(handles);
% Hints: get(hObject,'String') returns contents of L1 as text
%        str2double(get(hObject,'String')) returns contents of L1 as a double

function L2_Callback(hObject, eventdata, handles)
% hObject    handle to L2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
getpara(handles);
% Hints: get(hObject,'String') returns contents of L2 as text
%        str2double(get(hObject,'String')) returns contents of L2 as a double

% --- Executes on button press in solution.
function solution_Callback(hObject, eventdata, handles)
% hObject    handle to solution (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
open('inviscidChannelScript.m');


text:
function [res] = mne_add_coil_defs(chs,accuracy,templates)
%
% [res] = mne_add_coil_defs(chs,accuracy,coil_def_templates)
%
% Add transformed coil definitions to the channel info
%
% chs        - original channel definitions
% accuracy   - desired accuracy (0, 1, or 2, defaults to 1)
% templates  - coil definition templates
%              (defaults to $MNE_ROOT/setup/mne/coil_def.dat or $MNE_ROOT/share/mne/coil_def.dat)
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%   Revision 1.8  2006/05/03 18:53:05  msh
%   Approaching Matlab 6.5 backward compatibility
%
%   Revision 1.7  2006/04/26 00:14:24  msh
%   Return empty structure from fiff_pick_channels evoked and fiff_pick_types_evoked if there is no match.
%   Accuracy checking was incorrect in mne_add_coil_defs
%
%   Revision 1.6  2006/04/24 00:07:24  msh
%   Fixed error in the computation of the surface-based coordinate system.
%
%   Revision 1.5  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.4  2006/04/21 17:31:07  msh
%   Added the examples.
%   Modified the formatting of some informative output.
%
%   Revision 1.3  2006/04/18 20:44:46  msh
%   Added reading of forward solution.
%   Use length instead of size when appropriate
%
%   Revision 1.2  2006/04/17 15:01:34  msh
%   More small improvements.
%
%   Revision 1.1  2006/04/17 11:52:15  msh
%   Added coil definition stuff
%
%
me='MNE:mne_add_coil_defs';

global FIFF;
if isempty(FIFF)
    FIFF = fiff_define_constants();
end

if nargin == 2 || nargin == 1 || isempty(templates)
    %
    %   We need the templates
    %
    try
        templates = mne_load_coil_def;
    catch
        error(me,'%s',mne_omit_first_line(lasterr));
    end
    if nargin == 1
        accuracy = 1;
    end
elseif nargin ~= 3
    error(me,'Incorrect number of arguments');
end

fprintf(1,'\t%d coil definition templates available\n',length(templates));

if accuracy ~= 0 && accuracy ~= 1 && accuracy ~= 2
    error(me,'Accuracy should attain one of the values 0, 1, 2');
end

res    = chs;
nchan  = length(chs);
ntemp  = length(templates);
nmeg   = 0;
neeg   = 0;
eeg_cm = zeros(1,3);

first = true;
for k = 2:nchan
    if chs(k).kind == FIFF.FIFFV_EEG_CH || chs(k).kind == FIFF.FIFFV_MEG_CH || ...
            chs(k).kind == FIFF.FIFFV_REF_MEG_CH
        if first
            coord_frame = chs(k).coord_frame;
            first = false;
        elseif chs(k).coord_frame ~= coord_frame
            error(me,'All coils and electrodes are not in the same coordinate system (%s)',chs(k).ch_name);
        end
    end
end

lower_half = hex2dec('FFFF');
for k = 1:nchan
    res(k).coil_def = [];
    if chs(k).kind == FIFF.FIFFV_EEG_CH && ~isempty(chs(k).eeg_loc)
        %
        %   Set up a "coil definition" for an EEG electrode
        %
        temp             = templates(1);             %   Just pick one
        temp.class       = 1000;                     %   EEG electrode class
        temp.id          = chs(k).coil_type;         %   Copy from the channel info
        temp.accuracy    = 1;                        %   Normal accuracy no matter what
        temp.num_points  = size(chs(k).eeg_loc,2);   %   Number of points is normally two
        temp.size        = 10e-3;                    %   Rather arbitrary
        temp.baseline    = 0;                        %   Baseline does not apply
        temp.description = 'EEG electrode';          %   Useful description
        temp.coildefs    = zeros(temp.num_points,7); %   Location data
        temp.coildefs(1,1)   = 1;                    %   The electrode location
        temp.coildefs(1,2:4) = chs(k).eeg_loc(:,1)';
        nn  = chs(k).eeg_loc(:,1);                   %   The normal definition is probably immaterial
        len = nn'*nn;
        if len > 0
            temp.coildefs(1,5:7) = nn';
        end
        if temp.num_points > 1                       %   Add the reference electrode
            temp.coildefs(2,1)   = -1;
            temp.coildefs(2,2:4) = chs(k).eeg_loc(:,2);
            nn  = chs(k).eeg_loc(:,2);
            len = nn'*nn;
            if len > 0
                temp.coildefs(2,5:7) = nn';
            end
        end
        temp.FV          = [];                            %   Empty for now
        eeg_cm           = eeg_cm + temp.coildefs(1,2:4); %   Accumulate center of mass
        res(k).coil_def  = temp;                          %   There it is
        neeg = neeg + 1;
    elseif chs(k).kind == FIFF.FIFFV_MEG_CH || chs(k).kind == FIFF.FIFFV_REF_MEG_CH
        %
        %   Set up an MEG coil definition
        %
        coil_type = bitand(double(chs(k).coil_type),lower_half);
        temp = [];
        for t = 1:ntemp
            if templates(t).id == coil_type && templates(t).accuracy == accuracy
                temp = templates(t);
                break;
            end
        end
        %
        %  Did we find the template?
        %
        if isempty(temp)
            error(me,'Could not find an MEG coil template (coil type = %d accuracy = %d) for channel %s', ...
                coil_type,accuracy,chs(k).ch_name);
        end
        %
        %  Transform the template using the coil transformation
        %
        np = size(temp.coildefs,1);
        %
        %  Transform the integration points and the normals
        %
        trans = chs(k).coil_trans(1:3,:);
        temp.coildefs(:,2:4) = (trans*[ temp.coildefs(:,2:4)' ; ones(1,np) ])';
        temp.coildefs(:,5:7) = (trans*[ temp.coildefs(:,5:7)' ; zeros(1,np) ])';
        %
        %  Transform the vertices
        %
        if ~isempty(temp.FV)
            nvert = size(temp.FV.vertices,1);
            temp.FV.vertices = (trans*[temp.FV.vertices' ; ones(1,nvert)])';
        end
        %
        %   Attach the coil definition to the result
        %
        res(k).coil_def = temp;
        nmeg = nmeg + 1;
    end
end
%
%     Create the EEG coils drawings now that we have the center of mass
%
if neeg > 0
    eeg_cm = eeg_cm/neeg;
    for k = 1:nchan
        if res(k).kind == FIFF.FIFFV_EEG_CH && ~isempty(res(k).eeg_loc)
            res(k).coil_def.FV = draw_electrode(res(k).coil_def,eeg_cm);
        end
    end
end

fprintf(1,'\t%d MEG coil definitions and %d EEG electrodes set up\n',nmeg,neeg);

return;

    function FV = draw_electrode(def,cm)

        % create a patch for an EEG electrode

        if def.class == 1000
            % round disk
            Radius  = def.size/2;           % radius
            Len_cir = 15;                   % number of points for circle
            r0      = def.coildefs(1,2:4);
            %
            %   Make a template circle
            %
            rr(1,:) = cos(2*pi*[0:Len_cir]/Len_cir);
            rr(2,:) = sin(2*pi*[0:Len_cir]/Len_cir);
            rr(3,:) = 0.0;
            rr(4,:) = ones(1,Len_cir+1);
            %
            %    Scaling
            %
            T1 = [ [ Radius 0 0 0 ] ; [ 0 Radius 0 0 ] ; [ 0 0 1 0 ...
                ] ; [  0 0 0 1 ] ];
            %
            %   Rotation (this is only approximately correct...)
            %
            %   The normal is parallel to the vector from CM to the electrode location
            %
            nn     = (def.coildefs(1,2:4) - cm)';
            sizenn = norm(nn);
            if sizenn > 0
                nn = nn/sizenn;
                [ u, s, v ]  = svd(eye(3,3) - nn*nn');
                if nn'*u(:,3) < 0
                    u = - u;
                end
                T2 = [ u [ 0 0 0 ]' ; [ 0 0 0 1 ]];
            else
                T2 = eye(4,4);
            end
            %
            %   Translation
            %
            T3 = [ [ 1 0 0 r0(1) ] ; [ 0 1 0 r0(2) ] ; [ 0 0 1 r0(3) ...
                ] ; [  0 0 0 1 ] ];
            %
            %    Transfrom the data and create the structure
            %
            rr = T3*T2*T1*rr;
            FV.vertices = [[rr(1,:)' rr(2,:)' rr(3,:)']];
            FV.faces    = [1:length(FV.vertices)];
        end
        return;
    end

end


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
% ##2
%==============================================================================
%
%   - data                 Hand, Omega=(0,20)x(0,25),
%   - viewer               viewImage2D
%   - interpolation        splineInter
%   - distance             SSD
%   - transformation       affine2D
%   - regularizer          mfElastic
%   - optimizer            Gauss-Newton
% ===============================================================================

close all, help(mfilename);

% load data, set viewer, interpolator, transformation, distance
setup2DHNSPData
imgModel('reset','imgModel','splineInter','regularizer','moments','theta',1e-2);
distance('reset','distance','SSD');
trafo('reset','trafo','affine2D');
regularizer('reset','regularizer','mfElastic','alpha',5e2,'mu',1,'lambda',0);

% run MLIR
[yc,wc,his] = MLIR(ML,'maxLevel',8);

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


text:
function p = prior_gaussian(varargin)
%PRIOR_GAUSSIAN  Gaussian prior structure
%
%  Description
%    P = PRIOR_GAUSSIAN('PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    creates Gaussian prior structure in which the named
%    parameters have the specified values. Any unspecified
%    parameters are set to default values.
%
%    P = PRIOR_GAUSSIAN(P, 'PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    modify a prior structure with the named parameters altered
%    with the specified values.
%
%    Parameters for Gaussian prior [default]
%      mu       - location [0]
%      s2       - scale squared (variance) [1]
%      mu_prior - prior for mu [prior_fixed]
%      s2_prior - prior for s2 [prior_fixed]
%
%  See also
%    PRIOR_*

% Copyright (c) 2000-2001,2010 Aki Vehtari
% Copyright (c) 2010 Jaakko Riihim�ki

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'PRIOR_GAUSSIAN';
  ip.addOptional('p', [], @isstruct);
  ip.addParamValue('mu',0, @(x) isscalar(x) && x>0);
  ip.addParamValue('mu_prior',[], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('s2',1, @(x) isscalar(x) && x>0);
  ip.addParamValue('s2_prior',[], @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  p=ip.Results.p;

  if isempty(p)
    init=true;
    p.type = 'Gaussian';
  else
    if ~isfield(p,'type') && ~isequal(p.type,'Gaussian')
      error('First argument does not seem to be a valid prior structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('mu',ip.UsingDefaults)
    p.mu = ip.Results.mu;
  end
  if init || ~ismember('s2',ip.UsingDefaults)
    p.s2 = ip.Results.s2;
  end
  % Initialize prior structure
  if init
    p.p=[];
  end
  if init || ~ismember('mu_prior',ip.UsingDefaults)
    p.p.mu=ip.Results.mu_prior;
  end
  if init || ~ismember('s2_prior',ip.UsingDefaults)
    p.p.s2=ip.Results.s2_prior;
  end

  if init
    % set functions
    p.fh.pak = @prior_gaussian_pak;
    p.fh.unpak = @prior_gaussian_unpak;
    p.fh.lp = @prior_gaussian_lp;
    p.fh.lpg = @prior_gaussian_lpg;
    p.fh.recappend = @prior_gaussian_recappend;
  end
end

function [w, s] = prior_gaussian_pak(p)

  w=[];
  s={};
  if ~isempty(p.p.mu)
    w = p.mu;
    s=[s; 'Gaussian.mu'];
  end
  if ~isempty(p.p.s2)
    w = [w log(p.s2)];
    s=[s; 'log(Gaussian.s2)'];
  end
end

function [p, w] = prior_gaussian_unpak(p, w)

  if ~isempty(p.p.mu)
    i1=1;
    p.mu = w(i1);
    w = w(i1+1:end);
  end
  if ~isempty(p.p.s2)
    i1=1;
    p.s2 = exp(w(i1));
    w = w(i1+1:end);
  end
end

function lp = prior_gaussian_lp(x, p)

  lp = 0.5*sum(-log(2*pi) -log(p.s2)- 1./p.s2 .* sum((x-p.mu).^2,1));

  if ~isempty(p.p.mu)
    lp = lp + p.p.mu.fh.lp(p.mu, p.p.mu);
  end
  if ~isempty(p.p.s2)
    lp = lp + p.p.s2.fh.lp(p.s2, p.p.s2) + log(p.s2);
  end
end

function lpg = prior_gaussian_lpg(x, p)

  lpg = (1./p.s2).*(p.mu-x);

  if ~isempty(p.p.mu)
    lpgmu = sum((1./p.s2).*(x-p.mu)) + p.p.mu.fh.lpg(p.mu, p.p.mu);
    lpg = [lpg lpgmu];
  end
  if ~isempty(p.p.s2)
    lpgs2 = (sum(-0.5*(1./p.s2-1./p.s2.^2.*(x-p.mu).^2 )) + p.p.s2.fh.lpg(p.s2, p.p.s2)).*p.s2 + 1;
    lpg = [lpg lpgs2];
  end
end

function rec = prior_gaussian_recappend(rec, ri, p)
% The parameters are not sampled in any case.
  rec = rec;
  if ~isempty(p.p.mu)
    rec.mu(ri,:) = p.mu;
  end
  if ~isempty(p.p.s2)
    rec.s2(ri,:) = p.s2;
  end
end


text:
%%********************************************************************
%% minEpts: minimum volume ellipsoid containing given
%%          points, V(:,1),...,V(:,m).
%%
%% max log(det(B))
%% s.t. || Bx + d || <= 1, for all x = V(:,1),...,V(:,m)
%%
%% [blk,At,C,b,OPTIONS,B,d] = minEpts(V,solve)
%%
%% For problem formulation, see
%% Vandenberghe, Boyd, Wu, Determinant maximization with linear
%% matrix inequalities, SIAM J. Matrix Analysis and Applications,
%% 19 (1998), pp. 499--533.
%%
%%********************************************************************

  function  [blk,At,C,b,OPTIONS,B,d] = minEpts(V,solve)

  if (nargin < 2); solve = 0; end;
%%
%% form data matrices
%%
  [p,m] =  size(V); N = (p+1)*m;
  blk{1,1} = 's'; blk{1,2} = (p+1)*ones(1,m);
  blk{2,1} = 's'; blk{2,2} = p;
%%
  count = 0;
  for j = 1:p
     s1 = V(j,:)'; i1 = [j:(p+1):(p+1)*(m-1)+j]; j1 = [p+1:p+1:(p+1)*m];
     tmp = sparse(i1,j1,s1,N,N);
     tmp = tmp + tmp';
     count = count + 1;
     F{1,count} = -tmp;
     F{2,count} = sparse(j,j,-1,p,p);
  end
  for j = 2:p
     for k = 1:j-1
        s1 = V(k,:)'; i1 = [j:(p+1):(p+1)*(m-1)+j]; j1 = [p+1:p+1:(p+1)*m];
        s2 = V(j,:)'; i2 = [k:(p+1):(p+1)*(m-1)+k]; j2 = [p+1:p+1:(p+1)*m];
        tmp  = sparse(i1,j1,s1,N,N);
        tmp  = tmp + sparse(i2,j2,s2,N,N);
        tmp = tmp + tmp';
        count = count + 1;
        F{1,count} = -tmp;
        F{2,count} = sparse([j k],[k j],[-1 -1],p,p);
     end
  end
  for j = 1:p
     s1 = ones(m,1);
     i1 = [j:(p+1):(p+1)*(m-1)+j]; j1 = [p+1:p+1:(p+1)*m];
     tmp = sparse(i1,j1,s1,N,N);
     tmp = tmp + tmp';
     count = count + 1;
     F{1,count} = -tmp;
     F{2,count} = sparse(p,p);
  end
  At = svec(blk,F,ones(2,1));
  C{1,1} = speye(N);
  C{2,1} = zeros(p);
  b = zeros(p*(p+1)/2+p,1);
  parbarrier{1,1} = 0;
  parbarrier{2,1} = 1;
  OPTIONS.parbarrier = parbarrier;
%%
%% || Bx + d || <= 1.
%% x'*(B'*B)*x + 2(B'*d)'*x  + d'*d <= 1.
%%
  if (solve)
     [obj,X,y,Z,info] = sqlp(blk,At,C,b,OPTIONS);
     if (length(y) ~= p*(p+3)/2)
        error('length of y not compatible with p');
     end
     B = diag(y(1:p));
     tmp = p;
     for k = 1:p-1
         B(k+1:p,k) = y(tmp + [1:p-k]);
         B(k,k+1:p) = B(k+1:p,k)';
         tmp = tmp + p-k;
     end;
     d = y(tmp+1:length(y));
  else
     B = []; d = [];
  end
%%********************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%  Figure 7.51      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% script for fig. 7.51
%
clf;
dp=[1 1 0];
np=[1];
nc=conv([1 1.001],[8.32 0.8]);
dc=conv([1 4.08],[1 0.0196]);
num=conv(np,nc);
den=conv(dp,dc);
sys=tf(num,den);
rlocus(sys);


text:
function [Handles] = ctf_channel_plot(ctf,surf,CHAN,Handles)

% ctf_channel_plot - plot ctf.sensor data
%
% [Handles] = ctf_channel_plot(ctf [,surf][,CHAN][,H])
%
% ctf  - a struct returned by ctf_read
% surf - a surface tesselation struct, surf.vertices and surf.faces;
%        the default is to create a convex hull of the sensor locations
% CHAN - see ctf_channel_select for options (default, CHAN = 'meg')
%
% Handles - handles struct, eg:
%           H.Hf - figure handle to plot into (default, H.Hf = figure;)
%           H.Hs - sensor handles
%           H.Hp - patch handles
%           H.Hl - light handles
%
% All inputs, except ctf, are optional.  An empty input [] will be replaced
% with the default.
%
%      <>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%      <                                                      >
%      <                    DISCLAIMER:                       >
%      <                                                      >
%      < THIS PROGRAM IS INTENDED FOR RESEARCH PURPOSES ONLY. >
%      < THIS PROGRAM IS IN NO WAY INTENDED FOR CLINICAL OR   >
%      <                    OFFICIAL USE.                     >
%      <                                                      >
%      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>
%

% $Revision: 1.1 $ $Date: 2009-01-30 03:49:26 $

% Copyright (C) 2004  Darren L. Weber
%
% 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.

% Modified: 05/2004, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('ctf','var'), error('no ctf input'); end
if isempty(ctf), error('no ctf input'); end

% check the figure handle for initialization
if ~exist('Handles','var'),
  Handles.Hf = figure;
  colordef(Handles.Hf,'white')
end
if isempty(Handles),
  Handles.Hf = figure;
  colordef(Handles.Hf,'white')
end
if ~isfield(Handles,'Hf'),
  Handles.Hf = figure;
  colordef(Handles.Hf,'white')
end
if isempty(Handles.Hf),
  Handles.Hf = figure;
  colordef(Handles.Hf,'white')
end

% check other fields of the Handles
if ~isfield(Handles,'Hp'), Handles.Hp = {}; end
if ~isfield(Handles,'Hs'), Handles.Hs = {}; end
if ~isfield(Handles,'Hl'), Handles.Hl = []; end

usesurf = 1; % assume input surf
if ~exist('surf','var'), usesurf = 0; end
if isempty(surf), usesurf = 0; end

if ~exist('CHAN','var'), CHAN = 'meg'; end
if CHAN,
  % Get the channel indices to plot
  [CHAN,type] = ctf_channel_select(ctf,CHAN);
  sensors.vertices = zeros(length(CHAN),3);
  for i = 1:length(CHAN),
    sensors.vertices(i,:) = ctf.sensor.info(CHAN(i)).location(:,1)';
  end
end

ver = '$Revision: 1.1 $ $Date: 2009-01-30 03:49:26 $';
fprintf('\nCTF_CHANNEL_PLOT [v %s]\n',ver(11:15));

figure(Handles.Hf); hold on
set(Handles.Hf,'Renderer','OpenGL')

if CHAN,
  % plot the sensor cloud
  Handles.Hs{end+1} = scatter3(sensors.vertices(:,1),...
                               sensors.vertices(:,2),...
                               sensors.vertices(:,3),...
                               50,'bo','filled');
end

if ~usesurf,
  % nothing provided, assume we want
  % to create a convex hull of the sensors
  surf.vertices = sensors.vertices;
  surf.faces = convhulln(surf.vertices);
end
Nvert = length(surf.vertices);

% assume surface is a skin color
skinRGB = [1, 0.85, 0.75];
vertexColors = ones(Nvert,1) * skinRGB;

% plot the surface patch (eg, sensors or scalp)
Handles.Hp{end+1} = patch('faces',surf.faces,'vertices',surf.vertices,...
                          'facecolor','interp',...
                          'FaceVertexCData', vertexColors,...
                          'edgecolor', 'none',...
                          'facealpha',0.5);

% view parameters
view([40,20])
daspect([1 1 1])
axis equal
axis tight
grid on

lighting phong
material dull

% only add a light if it doesn't exist already
if ~isfield(Handles,'Hl'),
  Handles.Hl = camlight;
elseif isempty(Handles.Hl),
  Handles.Hl = camlight;
end

% The CTF MEG Head Coordinate System is
% +X anterior, +Y left, +Z superior
xlabel('x (cm)')
ylabel('y (cm)')
zlabel('z (cm)')

rotate3d on

return

author: PatternRecognition

repo: OpenBMI


text:
function matched_points=it_features_matching(upl_1, lor_1, upl_2, lor_2, th, cal, matches, past_weigth)

%% ITERATIVE FEATURES MATCHING
%% As features are selected the best matching in any iteration.
%% Then it is recursively refined till reaching the threshold size
%% subwindows

%% Input
%% iml,imr       : templates to be matched
%% upl_1, upl_2  : image coordinates of the templates' upper left corner
%% lor_1, lor_2  : image coordinates of the templates' lower right corner
%% th            : matching threshold
%% cal           : (1=rectified images with epipolar horizontal lines; 0=uncalibrated images )
%% matches       : N x 4 matrix comprised of the already matched points
%% past weigth   : present match weigth

%% Output
%% matched_points : Correspondences matrix (N x 4)
%%                  [v1_1 u1_1 v2_1 u2_1
%%                   v1_2 u1_2 v2_2 u2_2
%%                  ...
%%                   v1_N u1_N v1_N u1_N]

global IM1 IM2 IM1_X IM1_Y IM1_X2 IM1_Y2 SIGMA_MAX;
offset=0;
imw1=IM1(upl_1(1):lor_1(1),upl_1(2):lor_1(2));
imw2=IM2(upl_2(1):lor_2(1),upl_2(2):lor_2(2));

[mt1,mt2,upl_mt1,lor_mt1,upl_mt2,lor_mt2,pres_weigth,err]=temp_temp_matching(imw1,imw2,upl_1,lor_1,upl_2,lor_2,th,cal);
pres_weigth=pres_weigth+past_weigth;
[tysize,txsize]=size(mt1);

if (~err)  %% returned valid match

    if (tysize<8 | txsize<8)  %% last step. The central position of the matched subwindows are added to the correspondences set.
        matched_points=[matches; (upl_mt1+lor_mt1)/2, (upl_mt2+lor_mt2)/2, pres_weigth];

    else
        %% matched subwindows tiling
        y1=floor(tysize/2)-offset;
        y2=floor(tysize/2)+offset;
        x1=floor(txsize/2)-offset;
        x2=floor(txsize/2)+offset;

        %hi-le templates
        t1_1=mt1( 1:y2 , 1:x2 );
        upl_1_1=upl_mt1 + [0,0];
        lor_1_1=upl_mt1 + [y2-1,x2-1];

        t2_1=mt2( 1:y2  , 1:x2 );
        upl_2_1=upl_mt2 + [0,0];
        lor_2_1=upl_mt2 + [y2-1,x2-1];

        %hi-ri templates
        t1_2=mt1( 1:y2 , x1:txsize );
        upl_1_2=upl_mt1 + [0,x1-1];
        lor_1_2=upl_mt1 + [y2-1,txsize-1];

        t2_2=mt2( 1:y2 , x1:txsize );
        upl_2_2=upl_mt2 + [0,x1-1];
        lor_2_2=upl_mt2 + [y2-1,txsize-1];

        %lo-le templates
        t1_3=mt1( y1:tysize , 1:x2 );
        upl_1_3=upl_mt1 + [y1-1,0];
        lor_1_3=upl_mt1 + [tysize-1,x2-1];

        t2_3=mt2( y1:tysize , 1:x2 );
        upl_2_3=upl_mt2 + [y1-1,0];
        lor_2_3=upl_mt2 + [tysize-1,x2-1];

        %lo-ri templates
        t1_4=mt1( y1:tysize , x1:txsize );
        upl_1_4=upl_mt1 + [y1-1,x1-1];
        lor_1_4=upl_mt1 + [tysize-1,txsize-1];

        t2_4=mt2( y1:tysize , x1:txsize );
        upl_2_4=upl_mt2 + [y1-1,x1-1];
        lor_2_4=upl_mt2 + [tysize-1,txsize-1];

        matched_points=[matches;
                             it_features_matching(upl_1_1,lor_1_1,upl_2_1,lor_2_1,th,cal,matches,pres_weigth);
                             it_features_matching(upl_1_2,lor_1_2,upl_2_2,lor_2_2,th,cal,matches,pres_weigth);
                             it_features_matching(upl_1_3,lor_1_3,upl_2_3,lor_2_3,th,cal,matches,pres_weigth);
                             it_features_matching(upl_1_4,lor_1_4,upl_2_4,lor_2_4,th,cal,matches,pres_weigth)];

     end

 else
     matched_points=matches;
 end


text:
function [expressionRxns, parsedGPR, gene_used, signifRxns] = mapExpressionToReactions(model, expressionData, minSum)
% Determines the expression data associated to each reaction present in
% the model
%
% USAGE:
%
%    [expressionRxns parsedGPR, gene_used] = mapExpressionToReactions(model, expressionData)
%    [expressionRxns, parsedGPR, gene_used, signifRxns] =  mapExpressionToReactions(model, expressionData, minSum)
%
% INPUTS:
%	model                   model strusture
%	expressionData          mRNA expression data structure
%       .gene               	cell array containing GeneIDs in the same
%                               format as model.genes
%       .value                  Vector containing corresponding expression
%                               value (FPKM/RPKM)
%       .sig:               [optional field] Vector containing significance values of
%                           expression corresponding to expression values in
%                           expressionData.value (ex. p-values)
%
% OPTIONAL INPUT:
%    minSum:         instead of using min and max, use min for AND and Sum
%                    for OR (default: false, i.e. use min)
%
% OUTPUTS:
%   expressionRxns:         n x 1 non-negative value for reaction expression, corresponding to model.rxns.
%                           expressionRxns(j) is NaN when there is no expression data for the genes corresponding to reaction j.
%   parsedGPR:              cell matrix containing parsed GPR rule
%   gene_used:              gene identifier, corresponding to model.rxns, from GPRs
%                           whose value (expression and/or significance) was chosen for that
%                           reaction
%
% OPTIONAL OUTPUTS:
%   signifRxns:              significance of reaction expression, corresponding to model.rxns.

%
% Authors:
%       - Anne Richelle, May 2017 - integration of new extraction methods
%       - Chaitra Sarathy, Oct 2019, add significance value as optional input

if ~exist('minSum','var')
    minSum = false;
end

if isfield(expressionData, 'sig')
    exprSigFlag = 1;
else
    exprSigFlag = 0;
end

% Extracting GPR data from model
parsedGPR = GPRparser(model,minSum);

if exprSigFlag == 0

    % Find wich genes in expression data are used in the model
    % Returns vectors of gene identifiers and corresponding gene expression
    % levels for each gene present in the model ('model.genes').
    [gene_id, gene_expr] = findUsedGenesLevels(model,expressionData);

    % Link the gene to the model reactions
    % Map gene expression to reaction expression using the GPR rules. An AND
    % will be replaced by MIN and an OR will be replaced by MAX.
    [expressionRxns,  gene_used] = selectGeneFromGPR(model, gene_id, gene_expr, parsedGPR, minSum);

else

    [gene_id, gene_expr, gene_sig] = findUsedGenesLevels(model, expressionData);
    [expressionRxns,  gene_used, signifRxns] = selectGeneFromGPR(model, gene_id, gene_expr, parsedGPR, minSum, gene_sig);

end

author: opencobra

repo: cobratoolbox


text:
function walshMatrix=walsh(N)
hadamardMatrix = hadamard(N);
HadIdx = 0:N-1;                          % Hadamard index
M = log2(N)+1;                           % Number of bits to represent the index

binHadIdx = fliplr(dec2bin(HadIdx,M))-'0'; % Bit reversing of the binary index
binSeqIdx = zeros(N,M-1);                  % Pre-allocate memory
for k = M:-1:2
    % Binary sequency index
    binSeqIdx(:,k) = xor(binHadIdx(:,k),binHadIdx(:,k-1));
end
SeqIdx = binSeqIdx*pow2((M-1:-1:0)');    % Binary to integer sequency index
walshMatrix = hadamardMatrix(SeqIdx+1,:); % 1-based indexing


text:
function y=Afun1(x);
if x<1
      y=x+1;
else
      y=1+1/x;
end

author: Eurus-Holmes