text:
% loadeeg() - load a binary data file in Neuroscan .eeg file format.
%
% Usage:
% >> signal = loadeeg(filename);
% >> [signal, accept, typeeeg, rt, response, chan_names, pnts, ...
%      ntrials, srate, xmin, xmax] = loadeeg( filename, chanlist, ...
%        triallist, typerange, accepttype, rtrange, responsetype);
%
% Inputs:
%   filename     - [string] Input Neuroscan .eeg file
%   chanlist     - [integer array] Only import selected channels
%                  Ex: [3,4:10] {default: 'all'}
%   triallist    - [integer array] Only import selected trials {default: import all}
%   typerange    - [integer array] Only import trials of selected type
%                  {default: import all}
%   accepttype   - [integer array] Only import trials with the selected
%                  'accept' field values {default: import all}
%   rtrange      - [float array] [min max] (ms) Only import trials with subject
%                  reaction times in this range {default: import all}
%   responsetype - [integer array] Only import trials with selected
%                  response type values {default: import all}
%   format       - ['short'|'int32'] data format. Neuroscan v4.3+ assume 32-bit data
%                  while older versions assume 16-bit. {default: 'short' = 16-bit}
%
% Outputs:
%   signal       - Output signal of size (trials, points)
%   accept       - [1/0] vector of values for the accept field (one per trial)
%   typeeeg      - [Integer] values for the accept type (size 1,trials)
%   rt           - [float] values for the accept rt (size trials)
%   response     - [Integer] values for the accept response (size 1,trials)
%   chan_names   - ['string' array] channel names
%   pnts         - Number of time points per trial
%   ntrials      - Number of trials
%   srate        - Sampling rate (Hz)
%   xmin         - Trial start latency (ms)
%   xmax         - Trial end latency (ms)
%
% Example:
%   % Load .eeg data into an array named 'signal'
%   >> [signal]=loadeeg( 'test.eeg' );
%   % Plot the signal in the first channel, first trial
%   >> plot( signal(1,:) );
%
% Author: Arnaud Delorme, CNL, Salk Institute, 2001
%
% See also: pop_loadeeg(), eeglab()

% .eeg binary file format
%   data are organised into an array of Number_of_electrode x (Number_of_points_per_trial*Number_of_sweeps)
%   for a file with 32 electrodes, 700 points per trial and 300 sweeps, the resulting array is
%   of 32 collumn and 700*300 rows (300 consecutive blocs of 700 points)

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

function [signal, accept, typeeeg, rt, response, chan_names, pnts, ...
	nsweeps, rate, xmin, xmax]=loadeeg( FILENAME, chanlist, ...
	TrialList, typerange, acceptype, rtrange, responsetype, format)
if nargin<1
	fprintf('Not enought arguments\n');
	help loadeeg
	return;
end;
if nargin<2 chanlist='all'; end;
if nargin<3 TrialList='all'; end;
if nargin<4 typerange='all'; end;
if nargin<5 acceptype='all'; end;
if nargin<6 rtrange  ='all'; end;
if nargin<7 responsetype='all'; end;
if nargin<8 format='short'; end;
format = lower(format);

if ~strcmpi(format, 'short') && ~strcmpi(format, 'int32'), error('loadeeg: format error'); end;

% open file for reading
% ---------------------
fid=fopen(FILENAME,'r','ieee-le');
if fid<0
	fprintf('Error LOADEEG: File %s not found\n', FILENAME);
	return;
end;

% read general part of the erp header and set variables
% -----------------------------------------------------
rev = char(fread(fid, 20, 'uchar'))'; % revision number
erp = fread(fid, 342, 'uchar');	% skip the firsts 362 from BOF (368 bytes in orginal comment?)
nsweeps = fread(fid, 1, 'ushort');	% number of sweeps
erp = fread(fid, 4, 'uchar'); 	% skip 4 bytes
pnts= fread(fid, 1, 'ushort');	% number of point per waveform
chan= fread(fid, 1, 'ushort');  % number of channels
erp = fread(fid, 4, 'uchar'); 	% skip 4 bytes
rate= fread(fid, 1, 'ushort');  % sample rate (Hz)
erp = fread(fid, 127, 'uchar');	% skip 127 bytes
xmin= fread(fid, 1, 'float32'); % in s
xmax= fread(fid, 1, 'float32'); % in s
erp = fread(fid, 387, 'uchar');	% skip 387 bytes

fprintf('number of channels         : %d\n', chan);
fprintf('number of points per trial : %d\n', pnts);
fprintf('sampling rate (Hz)         : %f\n', rate);
fprintf('xmin (s)                   : %f\n', xmin);
fprintf('xmax (s)                   : %f\n', xmax);

% read electrode configuration
% ----------------------------
fprintf('Electrode configuration\n');
for elec = 1:chan
   	channel_label_tmp = fread(fid, 10, 'uchar');
	chan_names(elec,:) = channel_label_tmp';
	for index = 2:9 if chan_names(elec,index) == 0 chan_names(elec,index)=' '; end; end;
	erp = fread(fid, 47-10, 'uchar');
	baseline(elec) = fread(fid, 1, 'ushort');
	erp = fread(fid, 10, 'uchar');
	sensitivity(elec) = fread(fid, 1, 'float32');
	erp = fread(fid, 8, 'uchar');
	calib(elec) = fread(fid, 1, 'float32');
	fprintf('%s: baseline: %d\tsensitivity: %f\tcalibration: %f\n', ...
            char(chan_names(elec,1:4)), baseline(elec), sensitivity(elec), calib(elec));
	factor(elec) = calib(elec) * sensitivity(elec) / 204.8;
end;
%fprintf('Electrode configuration\n');
%for elec = 1:chan
%	erp = fread(fid, 47, 'uchar');
%	baseline(elec) = fread(fid, 1, 'ushort');
%	erp = fread(fid, 10, 'uchar');
%	sensitivity(elec) = fread(fid, 1, 'float32');
%	erp = fread(fid, 8, 'uchar');
%	calib(elec) = fread(fid, 1, 'float32');
%	fprintf('baseline: %d\tsensitivity: %f\tcalibration: %f\n', baseline(elec), sensitivity(elec), calib(elec));
%	factor(elec) = calib(elec) * sensitivity(elec) / 204.8;
%end;

xsize    = chan * pnts;
buf_size = chan * pnts ;			% size in shorts

% set tags for conditions
% -----------------------
if isstr(chanlist) && strcmpi(chanlist, 'all'), chanlist = [1:chan]; end;
if isstr(TrialList) && strcmpi(TrialList, 'all'),	trialtagI     = 1; else trialtagI     = 0; end;
if isstr(acceptype) && strcmpi(acceptype, 'all'),	acceptagI     = 1; else acceptagI     = 0; end;
if isstr(typerange) && strcmpi(typerange, 'all'), 	typetagI      = 1; else typetagI      = 0; end;
if isstr(responsetype) && strcmpi(responsetype, 'all'),	responsetagI  = 1; else responsetagI  = 0; end;
if isstr(rtrange) && strcmpi(rtrange, 'all'), 	rttagI        = 1; else rttagI        = 0; end;

count_selected = 1;
fprintf('Reserving array (can take some time)\n');
if isstr(TrialList)
    signal = zeros( chan, pnts*nsweeps);
else
    signal = zeros( chan, pnts*length(TrialList));
end;
fprintf('Array reserved, scanning file\n');

for sweep = 1:nsweeps

	% read sweeps header
	% ------------------
	s_accept   = fread(fid, 1, 'uchar');
	s_type     = fread(fid, 1, 'ushort');
	s_correct  = fread(fid, 1, 'ushort');
	s_rt       = fread(fid, 1, 'float32');
	s_response = fread(fid, 1, 'ushort');
	s_reserved = fread(fid, 1, 'ushort');

	unreaded_buf = 1;

	% store the sweep or reject the sweep
	% -----------------------------------
	if trialtagI trialtag = 1;        else trialtag = ismember(sweep, TrialList); end;
	if acceptagI acceptag = 1;        else acceptag =  ismember(s_accept, acceptype); end;
	if typetagI  typetag  = 1; 	  else typetag  =  ismember(s_type, typerange); end;
	if responsetagI responsetag  = 1; else responsetag  = ismember(s_response, responsetype); end;
	if rttagI       rttag  = 1; 	  else rttag  =  ismember(s_rt, rtrange); end;

	if typetag
		if trialtag
			if acceptag
				if responsetag
					if rttag

						buf = fread(fid, [chan pnts], format);
						unreaded_buf = 0;

						% copy information to array
						% -------------------------
						accept(count_selected)   = s_accept;
						typeeeg(count_selected)  = s_type;
						rt(count_selected)       = s_rt;
						response(count_selected) = s_response;

						% demultiplex the data buffer and convert to microvolts
						% -----------------------------------------------------
						for elec = 1:chan
							buf(elec, :) = (buf(elec, :)-baseline(elec)-0.0)*factor(elec);
						end;
                        try
                            signal(:,[((count_selected-1)*pnts+1):count_selected*pnts]) = buf;
                            count_selected = count_selected + 1;
                            if not(mod(count_selected,10)) fprintf('%d sweeps selected out of %d\n', count_selected-1, sweep); end;
                        catch,
                            disp('Warning: File truncated, aborting file read.');
                            count_selected = count_selected + 1;
                            break;
                        end;
					end;
				end;
			end;
		end;
	end;

	if unreaded_buf
        if strcmpi(format, 'short'), fseek(fid, buf_size*2, 'cof');
        else                         fseek(fid, buf_size*4, 'cof');
        end;
    end;
end;
nsweeps = count_selected-1;
fclose(fid);

% restrincting array
% ---------------------------------------
fprintf('rereservation of variables\n');
signal = signal(chanlist, 1:(count_selected-1)*pnts);
chan_names = chan_names(chanlist,:);

return;

author: PatternRecognition

repo: OpenBMI


text:
function [Fc,Vc]=patchDetach(varargin)

% function [Fc,Vc]=patchDetach(F,V,scaleFactor)
%-------------------------------------------------------------------------
%
%
% Change log:
% 2010
% 2019/04/23 Updated function header text
%-------------------------------------------------------------------------

%%

[Fc,Vc]=scalePatch(varargin{:});

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function e=PerecptronTst(x,y,w,b);
%==========================================
% Testing phase
%==========================================
tic
[l,p]=size(x);
e=0; % number of test errors
for i=1:l
    xx=x(i,:); % take one row
    ey=xx*w+b; % apply the perceptron classification rule
    if ey>=0.5
       ey=1;
    else
       ey=0;
    end
    if y(i)~=ey;
       e=e+1;
    end;
end
toc

author: Sable

repo: mcbench-benchmarks


text:
function J = computeCost(traj_obj, gamma)

tau_vec = traj_obj.tau_vec;

A_yaw = augmentA_yaw(tau_vec);
% Permutation matrix to rearrange b
C_yaw = permutMat_yaw(tau_vec);
% Cost matrix
Q_yaw = augmentQ(tau_vec,2);
% C*A^(-T)*Q*A^(-1)*C^(T) = R
R_yaw = comp_R_appA(A_yaw,C_yaw,Q_yaw);
% Known derivatives including position.
bFyaw = bF_yaw(0, tau_vec);
b_srtdYaw = b_srtd_yaw(R_yaw, bFyaw);
% P_yaw = (C_yaw*A_yaw)\b_srtdYaw;
J_yaw = b_srtdYaw'*R_yaw*b_srtdYaw;

% Compute cost
J = traj_obj.J + J_yaw + gamma*sum(tau_vec);

end

author: yorgoon

repo: minimum-snap-geometric-control


text:
function sudoku_test ( )

%*****************************************************************************80
%
%% SUDOKU_TEST tests the SUDOKU library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2007
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SUDOKU_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the SUDOKU library.\n' );

  sudoku_test01 ( );
  sudoku_test02 ( );
  sudoku_test03 ( );
  sudoku_test04 ( );
  sudoku_test05 ( );
  sudoku_test06 ( );
  sudoku_test07 ( );
  sudoku_test08 ( );
  sudoku_test09 ( );
  sudoku_test10 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SUDOKU_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% DEMUSPSIVM1 Try the IVM on the USPS digits data with RBF kernel.

% IVM

dataSetName = 'usps';
experimentNo = 1;

randn('seed', 1e5)
rand('seed', 1e5)

[X, y, XTest, yTest] = mapLoadData(dataSetName);

capitalName = dataSetName;
capitalName(1) = upper(capitalName(1));

options = ivmOptions;
options.kern = {'rbf', 'lin', 'bias', 'white'};
options.numActive = 500;

mu = zeros(size(yTest));
varSigma = zeros(size(yTest));

tic
% Learn an IVM for each digit
for trainData = 0:9
  index = trainData+1;

  % Train the IVM.
  model = ivmRun(X, y(:, index), options);

  % Make prediction for this digit.
  [mu(:, index), varSigma(:, index)] = ivmPosteriorMeanVar(model, XTest);
  mu(:, index) = mu(:, index) + model.noise.bias;
  yPred = sign(mu(:, index));
  testError(index) = 1-sum(yPred==yTest(:, index))/size(yTest, 1);
  fprintf('Digit %d, test error %2.4f\n', trainData, testError(index));

  % Deconstruct IVM for saving.
  [kernStore{index}, noiseStore{index}, ...
   ivmInfoStore{index}] = ivmDeconstruct(model);
  save(['dem' capitalName num2str(experimentNo)], 'testError', ...
       'ivmInfoStore', 'kernStore', 'noiseStore')
end
overallTime = toc;

% Make prediction for all digits.
[void, yPred] = max(mu, [], 2);
[void, yTest] = max(yTest, [], 2);
yPred = yPred - 1;
yTest = yTest - 1;
overallError = 1 - sum(yPred == yTest)/size(yTest, 1);

confusMat = zeros(10);
for i = 1:length(yPred)
  confusMat(yPred(i)+1, yTest(i)+1) = confusMat(yPred(i)+1, yTest(i)+1) + 1;
end
save(['dem' capitalName num2str(experimentNo)], 'testError', ...
     'ivmInfoStore', 'kernStore', ...
     'noiseStore', 'overallError', ...
     'confusMat', 'overallTime');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function anat = nifti2mrVistaAnat(ni)
% Take a nifti and re-orient to mrVista anatomy convention
%
%   anat = nifti2mrVistaAnat(ni)
%
%   ni: can be a matrix, a nifti struct, or a path to a nifti file.
%
% Our preferred NIFTI format is [sagittal(L:R), coronal(P:A), axial(I:S)] format.
% mrLoadRet coords are in [axial(S:I), coronal(A:P), sagittal(L:R)] format.
% This function permutes our preferred NIFTI format into our mrLoadRet format.
%
% April, 2009: JW

% Check format of input argument (Ideally, we shoudl require NIFTI and not
% allow a data array with no header. But we don't want to break things.)
if isnumeric(ni) || islogical(ni), data = ni; end
if ischar(ni), ni = niftiRead(ni); end

% Reorient to RAS
if isstruct(ni),
    ni = niftiApplyCannonicalXform(ni);
    data = niftiGet(ni, 'data');
end

% Permute
anat = permute(data, [3 2 1]);

% then flip dims 1 and 2
anat = flip(flip(anat, 1),2);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function pendulum()
% Matlab code for simulating a simple pendulum

P.g = 9.81;     %(m/s^2) gravity acceleration
P.l = 1.0;      %(m) length of pendulum
P.c = 1.0;     %(1/s) viscous damping constant
P.m = 1.0;     %(kg) pendulum mass

th0 = 0.5;  %(rad) initial angle of pendulum, from -j axis
w0 = 0.0;   %(rad/s) initial angular rate of pendulum
z0 = [th0;w0];   %Initial state vector

tSpan = [0,4]; %(s) [start, end] times for simulation
dtMax = 0.01;  %(s) maximum allowable time step

nStep = ceil(diff(tSpan)/dtMax);  %Number of simulation steps
t = linspace(tSpan(1),tSpan(2),nStep);  %(s) time vector
z = zeros(2,nStep); %  Initialize the state matrix
z(:,1) = z0;

%Run the simulation, using Euler integration
for i=2:nStep
    dt = t(i)-t(i-1);
    zNow = z(:,i-1);
    z(:,i) = zNow + dt*dynamics(zNow,P);
end

%Generate a plot of the result:
figure(1); clf;
subplot(2,1,1); plot(t,z(1,:));
xlabel('Time (s)'); ylabel('Angle (rad)');
title('Simple Damped Pendulum');
subplot(2,1,2); plot(t,z(2,:));
xlabel('Time (s)'); ylabel('Rate (rad/s)');

end

function dz = dynamics(z,P)
% Compute the dynamics of a simple damped pendulum:

th = z(1,:); w = z(2,:); %Unpack the state vector

%Dynamics:
dth = w;
dw = -(P.g/P.l)*sin(th) - (P.c/(P.m*P.l*P.l))*w;

dz = [dth; dw];  %Pack up state derivative
end

author: MatthewPeterKelly

repo: dscTutorials


text:
function plot_DR30_miss(plot_obj,plot_args,legend_string)
% Plots the point on the minDCF curve at which the number of misses
% falls below 30 (see Doddington's Rule of 30).
% Inputs:
%   plot_args: A cell array of arguments to be passed to 'plot' that control
%     the appearance of the curve. See Matlab's help on 'plot' for information.
%   legend_string: Optional. A string to describe this curve in the legend.
%     This string will be prepended with the system name.

if ischar(plot_args)
    plot_args = {plot_args};
end

figure(plot_obj.fh);
assert(iscell(plot_args))
lh = plot(plot_obj.plo(plot_obj.dr30Miss),plot_obj.minDCF(plot_obj.dr30Miss),plot_args{:});

if exist('legend_string','var') && ~isempty(legend_string)
    assert(ischar(legend_string))
    plot_obj.add_legend_entry(lh,legend_string,true);
end

end

author: nesl


text:
function err = getL2error(node,elem,uexact,uh,quadOrder,varargin)
%% GETL2ERROR L2 norm of the approximation error.
%
%  err = getL2error(node,elem,@uexact,uh) computes the L2 norm of the error
%  between the exact solution uexact and a finite element approximation uh
%  on a mesh described by node and elem.
%
%  The input parameter uexact is a function handle and uh is a column array
%  which could be:
%    - P0 element i.e. discontinuous and piecewise constant
%    - P1 element i.e. continuous and piecewise linear
%    - CR element i.e. piecewise linear and continuous at mid pts of edges
%    - P2 element i.e. continuous and piecewise quadratic
%    - P1+P0 element which could happen in the fluid application
%
%  err = getL2error(node,elem,@uexact,uh,quadOrder) computes error
%  using the quadrature rule with order quadOrder (up to 5). The default
%  order is 3.
%
%  Example: compute L2 error of piecewise linear interpolation
%
%     [node,elem] = squaremesh([0,1,0,1],0.25);
%     for k = 1:4
%         exactu = inline('sin(pi*pxy(:,1)).*sin(pi*pxy(:,2))','pxy');
%         uI = exactu(node);
%         N(k) = size(node,1);
%         err(k) = getL2error(node,elem,exactu,uI);
%         [node,elem] = uniformrefine(node,elem);
%     end
%     showrate(N,err);
%
% The cubic element is added by Jie Zhou.
%
% See also getH1error, getL2error3, getH1error3, quadpts.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

Nu = length(uh);    N = size(node,1);   NT = size(elem,1);
% Euler formula N-NE+NT = c % rough estimateus using Euler formula
NE = N + NT;    NP2 = N + NE;   NP3 = N + 2*NE + NT;
%% Default quadrature orders for different elements
if Nu > N+NT-5
    elem2dof = dofP2(elem);
    NP2 = max(elem2dof(:));
    NE = NP2 - N;
    NP3 = N+2*NE+NT;
end

%% Default quadrature orders for different elements
if ~exist('quadOrder','var')
    switch Nu
        case NT     % piecewise constant function P0
            quadOrder = 2;
        case N      % piecewise linear function P1 element
            quadOrder = 3;
        case NE     % piecewise linear function CR element
            quadOrder = 3;
        case N+NT   % piecewise linear function + constant function
            quadOrder = 3;
        case NP2    % piecewise quadratic function
            quadOrder = 4;
        case NE+NT  % weak Galerkin element
            quadOrder = 3;
        case NP3    % P3 element
            quadOrder = 5;
    end
end

%% compute L2 error element-wise using quadrature rule with order quadOrder
err = zeros(NT,1);
[lambda,weight] = quadpts(quadOrder);
% basis function at quadrature points
switch Nu
    case N    % P1 piecewise linear function
        phi = lambda; % linear bases
    case N+NT % P1+P0
        phi = lambda; % linear bases
    case NE  % CR nonconforming P1 element
        phi = 1-2*lambda;
        elem2edge = elem2dof(:,4:6) - N;
    case NP2 % P2 piecewise quadratic elements
        phi(:,1) = lambda(:,1).*(2*lambda(:,1)-1);
        phi(:,2) = lambda(:,2).*(2*lambda(:,2)-1);
        phi(:,3) = lambda(:,3).*(2*lambda(:,3)-1);
        phi(:,4) = 4*lambda(:,2).*lambda(:,3);
        phi(:,5) = 4*lambda(:,1).*lambda(:,3);
        phi(:,6) = 4*lambda(:,2).*lambda(:,1);
    case NE+NT  % weak Galerkin element
%             uhp = uh(1:NT); % only count the interior part
        phi = 1-2*lambda;
        elem2edge = elem2dof(:,4:6) - N + NT;
    case 2*NE+NT+N % P3 piecewise cubic elements
        phi(:,1) = 0.5*(3*lambda(:,1)-1).*(3*lambda(:,1)-2).*lambda(:,1);
        phi(:,2) = 0.5*(3*lambda(:,2)-1).*(3*lambda(:,2)-2).*lambda(:,2);
        phi(:,3) = 0.5*(3*lambda(:,3)-1).*(3*lambda(:,3)-2).*lambda(:,3);
        phi(:,4) = 9/2*lambda(:,3).*lambda(:,2).*(3*lambda(:,2)-1);
        phi(:,5) = 9/2*lambda(:,3).*lambda(:,2).*(3*lambda(:,3)-1);
        phi(:,6) = 9/2*lambda(:,1).*lambda(:,3).*(3*lambda(:,3)-1);
        phi(:,7) = 9/2*lambda(:,1).*lambda(:,3).*(3*lambda(:,1)-1);
        phi(:,8) = 9/2*lambda(:,1).*lambda(:,2).*(3*lambda(:,1)-1);
        phi(:,9) = 9/2*lambda(:,1).*lambda(:,2).*(3*lambda(:,2)-1);
        phi(:,10) = 27*lambda(:,1).*lambda(:,2).*lambda(:,3);
        elem2dof = dofP3(elem);
end
nQuad = size(lambda,1);
for p = 1:nQuad
    % evaluate uh at quadrature point
    switch Nu
        case NT   % P0 piecewise constant function
            uhp = uh;
        case N    % P1 piecewise linear function
            uhp = uh(elem(:,1))*phi(p,1) + ...
                  uh(elem(:,2))*phi(p,2) + ...
                  uh(elem(:,3))*phi(p,3);
        case N+NT % P1+P0
            uhp = uh(elem(:,1))*phi(p,1) + ...
                  uh(elem(:,2))*phi(p,2) + ...
                  uh(elem(:,3))*phi(p,3);
            uhp = uhp + uh(N+1:end);
        case NE  % CR nonconforming P1 element
            uhp = uh(elem2edge(:,1))*phi(p,1) + ...
                  uh(elem2edge(:,2))*phi(p,2) + ...
                  uh(elem2edge(:,3))*phi(p,3);
        case NP2 % P2 piecewise quadratic function
            uhp = uh(elem2dof(:,1)).*phi(p,1) + ...
                  uh(elem2dof(:,2)).*phi(p,2) + ...
                  uh(elem2dof(:,3)).*phi(p,3) + ...
                  uh(elem2dof(:,4)).*phi(p,4) + ...
                  uh(elem2dof(:,5)).*phi(p,5) + ...
                  uh(elem2dof(:,6)).*phi(p,6);
        case NP3
            uhp = uh(elem2dof(:,1)).*phi(p,1) + ...
                  uh(elem2dof(:,2)).*phi(p,2) + ...
                  uh(elem2dof(:,3)).*phi(p,3) + ...
                  uh(elem2dof(:,4)).*phi(p,4) + ...
                  uh(elem2dof(:,5)).*phi(p,5) + ...
                  uh(elem2dof(:,6)).*phi(p,6) + ...
                  uh(elem2dof(:,7)).*phi(p,7) + ...
                  uh(elem2dof(:,8)).*phi(p,8) + ...
                  uh(elem2dof(:,9)).*phi(p,9) + ...
                  uh(elem2dof(:,10)).*phi(p,10);
        case NE+NT  % weak Galerkin element
%             uhp = uh(1:NT); % only count the interior part
            uhp = uh(elem2edge(:,1))*phi(p,1) + ...
                  uh(elem2edge(:,2))*phi(p,2) + ...
                  uh(elem2edge(:,3))*phi(p,3);
    end
    % quadrature points in the x-y coordinate
    pxy = lambda(p,1)*node(elem(:,1),:) ...
        + lambda(p,2)*node(elem(:,2),:) ...
        + lambda(p,3)*node(elem(:,3),:);
    err = err + weight(p)*(uexact(pxy) - uhp).^2;
end
%% Modification
% area of triangles
ve2 = node(elem(:,1),:)-node(elem(:,3),:);
ve3 = node(elem(:,2),:)-node(elem(:,1),:);
area = 0.5*abs(-ve3(:,1).*ve2(:,2)+ve3(:,2).*ve2(:,1));
err = area.*err;
err(isnan(err)) = 0; % singular values, i.e. uexact(p) = infty, are excluded
err = sqrt(sum(err));

repo: Aero-matlab


text:
function varargout = cumsum2(varargin)
%CUMSUM2   Double indefinite integral of a DISKFUN.
%   This is not allowed and returns an error.  This function exists so that the
%   error message is meaningful to a DISKFUN user.

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

error('CHEBFUN:DISKFUN:CUMSUM2:fail',...
      'The indefinite integral is not supported for a diskfun.');
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [] = printTextSummaryForHMMData( data_struct, model )

numObj = length(data_struct);
Ts = [ data_struct(:).T];

switch model.obsModel.type
    case 'Multinomial'
        data_struct(1).numVocab = length(model.obsModel.params.lambda);
        nV = data_struct(1).numVocab;
        obsDescr = sprintf( 'with %d distinct symbols\n', nV );
        if isfield( data_struct, 'nEmissions' )
            Es = [ data_struct(:).nEmissions ] ;
            obsDescr = [obsDescr sprintf('\t   between %d-%d emissions per timestep ( median = %.1f ) \n', min(Es), max(Es), median(Es)) ];
        else
            symbolInfo = sprintf('\t   all timesteps emit exactly %d symbols\n',  size( data_struct(1).obs,1 ) );
            obsDescr = [obsDescr symbolInfo];
        end
    case 'Bernoulli'
        D = size( data_struct(1).obs, 2 );
        obsDescr = sprintf( 'across %d dimensions', D );
    case 'Gaussian'
        D = size( data_struct(1).obs, 2 );
        obsDescr = sprintf( 'with %d dimensions', D );
    case 'AR-Gaussian'
        D = size( data_struct(1).obs, 2 );
        obsDescr = sprintf( 'with %d dimensions', D );
end

fprintf( 'Dataset Summary: \n\t %d time series objects \n', numObj );
fprintf('\t   between %d-%d timesteps per object ( median = %.1f ) \n', min(Ts), max(Ts), median(Ts));
fprintf('\t   %s emissions %s\n', model.obsModel.type, obsDescr );

author: michaelchughes

repo: NPBayesHMM


text:
function [filt] = notchfilter(dat,Fs,Fl,N)

% NOTCHFILTER line noise reduction filter for EEG/MEG data
%
% [filt] = notchfilter(dat, Fsample, Fline)
%
% where
%   dat        data matrix (Nchans X Ntime)
%   Fsample    sampling frequency in Hz
%   Fline      line noise frequency (would normally be 50Hz)
%   N          optional filter order, default is 4
%
% if Fline is specified as 50, a band of 48-52 is filtered out
% if Fline is specified as [low high], that band is filtered out

% original      (c) 2003, Pascal Fries
% modifications (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$

if nargin<4
  % set the default filter order
  N = 4;
end

Nchans   = size(dat,1);
Nsamples = size(dat,2);

% use a digital FIR filter
Fn = Fs/2;           % Nyquist frequency
if length(Fl)==1
  % default use a notch-width of 2Hz in both directions
  % otherwise use the specified band
  Fl = [Fl-2 Fl+2];
end
[B, A] = butter(N, [min(Fl)/Fn max(Fl)/Fn], 'stop');
filt = filtfilt(B, A, dat')';

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%channel_estimation.m
% for LS/DFT Channel Estimation with linear/spline interpolation
clear all;
close all;
clf;
% figure(1), clf
% figure(2), clf
Nfft=512;
Ng=Nfft/8;
Nofdm=Nfft+Ng;
Nsym=100;
Nps=32;
Np=Nfft/Nps;
Nd=Nfft-Np; % Pilot spacing, Numbers of pilots and data per OFDM symbol
Nbps=4;
M=2^Nbps; % Number of bits per (modulated) symbol
mod_object = modem.qammod('M',M, 'SymbolOrder','gray');
demod_object = modem.qamdemod('M',M, 'SymbolOrder','gray');
Es=1;
A=sqrt(3/2/(M-1)*Es); % Signal energy and QAM normalization factor
%fs = 10e6;  ts = 1/fs;  % Sampling frequency and Sampling period
SNRs = [30];
sq2=sqrt(2);
for i=1:length(SNRs)
   SNR = SNRs(i);
   rand('seed',1);
   randn('seed',1);
   MSE = zeros(1,6);
   nose = 0;
   for nsym=1:Nsym
      Xp = 2*(randn(1,Np)>0)-1;    % Pilot sequence generation
      %Data = ((2*(randn(1,Nd)>0)-1) + j*(2*(randn(1,Nd)>0)-1))/sq2; % QPSK modulation
      msgint=randint(1,Nfft-Np,M);    % bit generation
      Data = modulate(mod_object,msgint)*A;
      %Data = modulate(mod_object, msgint);
      %Data = modnorm(Data,'avpow',1)*Data;   % normalization
      ip = 0;
      pilot_loc = [];
      for k=1:Nfft
         if mod(k,Nps)==1
            X(k) = Xp(floor(k/Nps)+1);
            pilot_loc = [pilot_loc k];
            ip = ip+1;
         else
             X(k) = Data(k-ip);
         end
      end
      x = ifft(X,Nfft);                            % IFFT
      xt = [x(Nfft-Ng+1:Nfft) x];                  % Add CP
      h = [(randn+j*randn) (randn+j*randn)/2];     % generates a (2-tap) channel
      H = fft(h,Nfft);
      channel_length = length(h); % True channel and its time-domain length
      H_power_dB = 10*log10(abs(H.*conj(H)));      % True channel power in dB
      y_channel = conv(xt, h);                     % Channel path (convolution)
      sig_pow = mean(y_channel.*conj(y_channel));
      %y_aw(1,1:Nofdm) = y(1,1:Nofdm) + ...
      %   sqrt((10.^(-SNR/10))*sig_pow/2)*(randn(1,Nofdm)+j*randn(1,Nofdm)); % Add noise(AWGN)
      yt = awgn(y_channel,SNR,'measured');
      y = yt(Ng+1:Nofdm);                   % Remove CP
      Y = fft(y);                           % FFT
      for m=1:3
         if m==1
             H_est = LS_CE(Y,Xp,pilot_loc,Nfft,Nps,'linear');
             method='LS-linear'; % LS estimation with linear interpolation
         elseif m==2
             H_est = LS_CE(Y,Xp,pilot_loc,Nfft,Nps,'spline');
             method='LS-spline'; % LS estimation with spline interpolation
         else
             H_est = MMSE_CE(Y,Xp,pilot_loc,Nfft,Nps,h,SNR);
             method='MMSE'; % MMSE estimation
         end
         H_est_power_dB = 10*log10(abs(H_est.*conj(H_est)));
         h_est = ifft(H_est);
         h_DFT = h_est(1:channel_length);
         H_DFT = fft(h_DFT,Nfft); % DFT-based channel estimation
         H_DFT_power_dB = 10*log10(abs(H_DFT.*conj(H_DFT)));
         if nsym==1
           figure(1)
           subplot(319+2*m)
           plot(H_power_dB,'b','linewidth',1);
           grid on;
           hold on;
           plot(H_est_power_dB,'r:+','Markersize',4,'linewidth',1);
           %axis([0 32 -6 10])
           title(method);
           xlabel('Subcarrier Index');
           ylabel('Power[dB]');
           legend('True Channel',method,4);
           set(gca,'fontsize',9)
           subplot(320+2*m)
           plot(H_power_dB,'b','linewidth',1);
           grid on;
           hold on;
           plot(H_DFT_power_dB,'r:+','Markersize',4,'linewidth',1);
           %axis([0 32 -6 10])
           title([method ' with DFT']);
           xlabel('Subcarrier Index');
           ylabel('Power[dB]');
           legend('True Channel',[method ' with DFT'],4);
           set(gca,'fontsize',9)
         end
         MSE(m) = MSE(m) + (H-H_est)*(H-H_est)';
         MSE(m+3) = MSE(m+3) + (H-H_DFT)*(H-H_DFT)';
      end
      Y_eq = Y./H_est;
      if nsym>=Nsym-10
          figure(2)
          subplot(221)
          plot(Y,'.','Markersize',5)
          axis([-1.5 1.5 -1.5 1.5])
          axis('equal')
          set(gca,'fontsize',9)
          hold on,
          subplot(222)
          plot(Y_eq,'.','Markersize',5)
          axis([-1.5 1.5 -1.5 1.5])
          axis('equal')
          set(gca,'fontsize',9)
          hold on,
      end
      ip = 0;
      for k=1:Nfft
         if mod(k,Nps)==1
             ip=ip+1;
         else
             Data_extracted(k-ip)=Y_eq(k);
         end
      end
      msg_detected = demodulate(demod_object,Data_extracted/A);
      nose = nose + sum(msg_detected~=msgint);
   end
   MSEs(i,:) = MSE/(Nfft*Nsym);
end
Number_of_symbol_errors=nose
figure(3)%, clf
semilogy(SNRs',MSEs(:,1),'-x', SNRs',MSEs(:,3),'-o')
legend('LS-linear','MMSE')
fprintf('MSE of LS-linear/LS-spline/MMSE Channel Estimation = %6.4e/%6.4e/%6.4e\n',MSEs(end,1:3));
fprintf('MSE of LS-linear/LS-spline/MMSE Channel Estimation with DFT = %6.4e/%6.4e/%6.4e\n',MSEs(end,4:6));

author: LyricYang


text:
function x = error_f_inverse ( y )

%*****************************************************************************80
%
%% ERROR_F_INVERSE inverts the error function ERF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real Y, the value of the error function.
%
%    Output, real X, the value such that ERROR_F(X) = Y.
%
  z = ( y + 1.0 ) / 2.0;

  x = normal_01_cdf_inv ( z );

  x = x / sqrt ( 2.0 );

  return
end


text:
function [S,SE] = sample_edges(V,E,samples_per_edge)
  % SAMPLE_EDGES Compute samples_per_edge extra points along each edge in E
  % defined over vertices of V.
  %
  % S = sample_edges(V,E,samples_per_edge)
  %
  % Inputs:
  %   V  vertices over which edges are defined, # vertices by dim
  %   E  edge list, # edges by 2
  %   samples_per_edge  number of extra samples to be computed along edge not
  %     including start and end points
  % Outputs:
  %   S  sampled vertices, size less than # edges * (2+samples_per_edge) by dim,
  %   always begins with V so that E is also defined over S
  %
  %

  dim = size(V,2);

  % trivial case
  if(isempty(E))
    S = [];
    SE = [];
  elseif(samples_per_edge < 0)
    S = V;
    SE = [];
  else
    % fraction parameter
    t = linspace(0.0,1.0,samples_per_edge+2);
    % get rid of start and end points
    t = t(2:(end-1));
    % repeat for each coordinate
    t = reshape(repmat(t,dim,1),1,size(t,2)*dim);
    % repeat for each edge
    t = repmat(t,size(E,1),1);
    % repeat start coords

    sp = repmat(V(E(:,1),:),1,samples_per_edge);
    % repeat end coords
    ep = repmat(V(E(:,2),:),1,samples_per_edge);
    % lerp from start point to end point for each coordinate
    S = sp.*(1-t) + ep.*t;
    % reshape to list coordinates
    S = S';
    S = reshape(S,dim,prod(size(S))/dim)';
    S = [V;S];
    % Determine edges between samples
    E1 = repmat((1:(samples_per_edge-1))',size(E,1),1);
    E2 = repmat((2:(samples_per_edge))',size(E,1),1);
    off = size(V,1) + ...
      reshape( ...
        repmat( ...
          samples_per_edge*((1:size(E,1))-1), ...
          samples_per_edge-1, ...
          1), ...
        (samples_per_edge-1)*size(E,1),...
        1);
    SE = repmat(off,1,2) + [E1 E2];
    % end point edges
    SE = [ ...
      SE; ...
      E(:,1) size(V,1) + (1+samples_per_edge*((1:size(E,1))-1))'; ...
      E(:,2) size(V,1) + (samples_per_edge*((1:size(E,1))))'];
  end
end

author: alecjacobson

repo: gptoolbox


text:
% Driver script for solving the 1D Euler equations
Globals1D;

% Polynomial order used for approximation
N = 4;

% Generate simple mesh
[Nv, VX, K, EToV] = MeshGen1D(0, 2*pi, 80);

% Initialize solver and construct grid and metric
StartUp1D;

% Set up initial conditions -- Sod's problem
MassMatrix = inv(V')/V;
%cx = ones(Np,1)*sum(MassMatrix*x,1)/2;

%u = ones(Np,K).*( sin(2*pi*cx) );
u0 = 0.5 + sin(x);
FinalTime = 1.5;

% Solve Problem
[u] = iBurgers1D(u0,FinalTime);
snapnow;

% Plot Figure
%plotrange=[0,2*pi,min(min(u0))-0.2,max(max(u0))+0.2];
%plot(x,u,x,u0,'-+'); axis(plotrange); ylabel('u(x)'); xlabel('x');


text:
function A = metric_05 ( p )

%*****************************************************************************80
%
%% METRIC_05 evaluates metric #5 at any point.
%
%  Discussion:
%
%    This routine evaluates the matrix that determines the metric
%    at a point.
%
%    It is not diagonal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 May 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P(2), the point at which the metric matrix is to
%    be evaluated.
%
%    Output, real A[2,2], the metric matrix.
%
  A = [ 2.0, 3.0; 3.0, 5.0 ];

  return
end


text:
function [image1,image2,img3]=mosaic_map(img1,img2,d)

[m1,n1,p1] = size(img1);
m11 = ceil(m1/d);
n11 = ceil(n1/d);
for i=1:2:m11
    for j=2:2:n11
        img1((i-1)*d+1:i*d,(j-1)*d+1:j*d,:)=0;
    end
end
for i=2:2:m11
    for j=1:2:n11
        img1((i-1)*d+1:i*d,(j-1)*d+1:j*d,:)=0;
    end
end
image1=img1(1:m1,1:n1,:);

%%
[m2,n2,p2] = size(img2);
m22 = ceil(m2/d);
n22 = ceil(n2/d);
for i=1:2:m22
    for j=1:2:n22
       img2((i-1)*d+1:i*d,(j-1)*d+1:j*d,:)=0;
    end
end
for i=2:2:m22
    for j=2:2:n22
       img2((i-1)*d+1:i*d,(j-1)*d+1:j*d,:)=0;
    end
end
image2=img2(1:m2,1:n2,:);
%%
img3=image1+image2;

end

author: LJY-RS

repo: RIFT-multimodal-image-matching


text:
function conv_layer_idx = get_conv_layer_idx_from_layer_idx(layer_idx)
    global config;
    c = 0;
    for m = 1:layer_idx
        if strfind(config.forward_pass_scheme{m}, 'conv')
            c = c + 1;
        end
    end
    conv_layer_idx = c;
end

author: jimmy-ren


text:
% Handles experimental focus swapping

function focus_swap(hObject, eventdata, handles)
    x_size     = str2double(get(handles.x_size_edit,'String'));
    y_size     = str2double(get(handles.y_size_edit,'String'));

    if x_size < y_size
        positionf = [.18,.065,.87,.87];
    else
        positionf = [.32,.3,.6,.65];
    end

    focused = getappdata(0,'focused');

    if focused == hObject %Return when clicking an object already in focus
        return;
    end

    current = getappdata(hObject,'pos');

    set(focused, 'Position', current);      %Swap the axes positions
    axis image
    setappdata(focused,'pos',current);

    set(hObject, 'Position', positionf);
    axis image
    setappdata(0,'focused', hObject);

    if get(handles.show_colorbar,'Value')
        pos = get(hObject, 'Position');
        colorbar('Position',[.947, pos(2), .015, pos(4)],'DrawMode','fast')
    end

    axes(focused)
    colorbar('delete');
    drawnow

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% M=2 scattering + freq scatt, mult Q1, cv parameters

run_name = 'DSS_Table2_TIMIT_m2_freq_multQ1';

src = phone_src('/path/to/timit');

[train_set,test_set,valid_set] = phone_partition(src);

N = 2^13;
T_s = 2560;

filt1_opt.filter_type = {'gabor_1d','morlet_1d'};
filt1_opt.Q = [8 1];
filt1_opt.J = T_to_J(512,filt1_opt);

sc1_opt.M = 2;

ffilt1_opt.filter_type = 'morlet_1d';
ffilt1_opt.J = 6;

fsc1_opt.M = 1;

Wop1 = wavelet_factory_1d(N, filt1_opt, sc1_opt);
fWop1 = wavelet_factory_1d(64, ffilt1_opt, fsc1_opt);

scatt_fun1 = @(x)(log_scat(renorm_scat(scat(x,Wop1))));
fscatt_fun1 = @(x)(func_output(@scat_freq,2,scatt_fun1(x),fWop1));
format_fun1 = @(x)(format_scat(fscatt_fun1(x)));

filt2_opt = filt1_opt;
filt2_opt.Q = [1 1];
filt2_opt.J = T_to_J(512,filt2_opt);

sc2_opt = sc1_opt;

ffilt2_opt = ffilt1_opt;
ffilt2_opt.J = 4;

fsc2_opt = fsc1_opt;

Wop2 = wavelet_factory_1d(N, filt2_opt, sc2_opt);
fWop2 = wavelet_factory_1d(16, ffilt2_opt, fsc2_opt);

scatt_fun2 = @(x)(log_scat(renorm_scat(scat(x,Wop2))));
fscatt_fun2 = @(x)(func_output(@scat_freq,2,scatt_fun2(x),fWop2));
format_fun2 = @(x)(format_scat(fscatt_fun2(x)));

duration_fun = @(x,obj)(32*duration_feature(x,obj));

features = {format_fun1, format_fun2, duration_fun};

for k = 1:length(features)
	fprintf('testing feature #%d...',k);
	tic;
	if nargin(features{k}) == 1
		sz = size(features{k}(randn(N,1)));
	else
		sz = size(features{k}(randn(N,1),struct('u1',1,'u2',N)));
	end
	aa = toc;
	fprintf('OK (%.2fs) (size [%d,%d])\n',aa,sz(1),sz(2));
end

database_opt.input_sz = N;
database_opt.output_sz = T_s;
database_opt.obj_normalize = 2;
database_opt.collapse = 1;

db = prepare_database(src,features,database_opt);
db.features = single(db.features);
db = svm_calc_kernel(db,'gaussian','triangle',[db.indices{train_set}]);

optt.kernel_type = 'gaussian';
optt.gamma = 2.^[-14:2:-10];
optt.C = 2.^[2:2:6];
optt.search_depth = 2;
optt.full_test_kernel = 1;

[dev_err_grid,C_grid,gamma_grid] = ...
	svm_adaptive_param_search(db,train_set,valid_set,optt);

[dev_err,ind] = min(dev_err_grid{end});
C = C_grid{end}(ind);
gamma = gamma_grid{end}(ind);

optt1 = optt;
optt1.C = C;
optt1.gamma = gamma;

model = svm_train(db,train_set,optt1);
labels = svm_test(db,model,test_set);
err = classif_err(labels,test_set,db.src);

save([run_name '.mat'],'labels','err','C','gamma');

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


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

%%
clear; close all; clc;

%% Syntax
% |setFEBioPath(FEBioPathSpec);|

%% Description
% This function sets the FEBio path defined by FEBioPathSpec. Setting the
% path means it is stored in the GIBBON configuration file:
% /GIBBON/config/FEBioPath.txt.
%
% See also |getFEBioPath|

%% Examples
% This example shows how the febio path can be specified and set for use in
% GIBBON

%%
% The FEBio path should be the full path to the FEBio executable file (and
% should include the extendion for Windows).

%%
% These are examples of paths defined for Windows and Linux (Ubuntu):
% Windows:
% C:\Program Files\febio-2.9.1\bin\FEBio2.exe
%
% Linux:
% /home/kevin/FEBio-2.9.1/bin/febio2

%%
% To avoid this demo from undoing the path settings for a user the current
% path is here instead uploaded with |getFEBioPath|

%Define FEBio path
febioPath=getFEBioPath; % Example: febioPath='/home/kevin/FEBio-2.9.1/bin/febio2'

%%
% Setting the FEBio path
setFEBioPath(febioPath);

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode


text:
function stla_to_tri_surface_fast ( stla_file_name, node_file_name, ...
  triangle_file_name )

%*****************************************************************************80
%
%% STLA_TO_TRI_SURFACE_FAST converts an ASCII STL file to a TRI_SURFACE file.
%
%  Discussion:
%
%    This function makes use of the "fast" routines STLA_SIZE_FAST and
%    STLA_READ_FAST, which may be much faster when reading large files.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string STLA_FILE_NAME, the name of the ASCII STL file to be read.
%
%    Input, string NODE_FILE_NAME, the name of the TRI_SURFACE node file to
%    be written.
%
%    Input, string TRIANGLE_FILE_NAME, the name of the TRI_SURFACE triangle
%    file to be written.
%

%
%  Get the size of the STLA data.
%
  [ solid_num, node_num, triangle_num, text_num ] = stla_size_fast ( stla_file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  STLA file "%s".\n', stla_file_name );
  fprintf ( 1, '  Number of solids   = %d\n', solid_num );
  fprintf ( 1, '  Number of nodes    = %d\n', node_num );
  fprintf ( 1, '  Number of faces    = %d\n', triangle_num );
  fprintf ( 1, '  Number of textures = %d\n', text_num );
%
%  Read the STLA data.
%
  [ node_xyz, triangle_node, triangle_normal ] = stla_read_fast ( stla_file_name, ...
    node_num, triangle_num );
%
%  Write the data to the TRI_SURFACE files.
%
  dim_num = size ( node_xyz, 1 );
  order_num = 3;

  tri_surface_write ( node_file_name, triangle_file_name, dim_num, node_num, ...
    order_num, triangle_num, node_xyz, triangle_node );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  TRI_SURFACE data written to:\n' );
  fprintf ( 1, '  node file     "%s".\n', node_file_name );
  fprintf ( 1, '  triangle_file "%s".\n', triangle_file_name );

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

%*****************************************************************************80
%
%% I4MAT_WRITE0 writes an I4MAT file with no header.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%
%    Input, logical HEADER, is TRUE if the header is to be included.
%

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

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE0 - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE0 - Error!' );
    return;
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

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

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

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

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE0 - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE0 - Error!' );
    return;
  end
%
%  Write the data.
%
%  For greater precision, try:
%
%     fprintf ( output_unit, '  %24,16f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %14f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function [ node_xyz, face_node, face_normal ] = stla_read_fast ( ...
  input_file_name, node_num, face_num )

%*****************************************************************************80
%
%% STLA_READ_FAST reads graphics information from an ASCII StereoLithography file.
%
%  Discussion:
%
%    This routine should be faster than STLA_READ, since it does little
%    error checking, and uses native MATLAB commands like SSCANF to
%    read the data directly.
%
%  Modified:
%
%    16 February 2007
%
%  Author:
%
%    Bill Seufzer,
%    NASA Langley
%
%  Reference:
%
%    3D Systems, Inc,
%    Stereolithography Interface Specification,
%    October 1989.
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Input, integer NODE_NUM, the number of vertices defined.
%
%    Input, integer FACE_NUM, the number of faces defined.
%
%    Output, real NODE_XYZ(3,NODE_NUM), the coordinates of points.
%
%    Output, integer FACE_NODE(3,FACE_NUM), the nodes that make up each face.
%
%    Output, real FACE_NORMAL(3,FACE_NUM), the normal vector
%    at each face.
%
  node_xyz = [];
  face_node = [];
  face_normal = [];
%
%  If no file input, try to get one from the user.
%
  if ( nargin < 1 )
    input_file_name = input ( 'Enter the name of the ASCII STL file.' );
    if ( isempty ( input_file_name ) )
      return
    end
  end
%
%  If no size input, scan the file.
%
  if ( nargin < 3 )
    [ solid_num, node_num, face_num, text_num ] = stla_size_fast ( ...
      input_file_name );
  end
%
%  Open the file.
%
  input_file_unit = fopen ( input_file_name, 'r' );

  if ( input_file_unit < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'STLA_READ_FAST - Fatal error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'STLA_READ_FAST - Fatal error!' )
    return
  end
%
%  Allocate space.
%
  solid = 0;
  node = 0;
  face = 0;

  node_xyz = zeros(3,node_num);
  face_node = zeros(3,face_num);
  face_normal = zeros(3,face_num);
%
%  Read the text, one line at a time.
%
  while ( 1 )

    text = fgetl ( input_file_unit );
%
%  On end of input, FGETL returns a -1 instead of a character array.
%
    if ( ~ischar ( text ) )
      break
    end
%
%  Left justify the string (no initial blanks).
%  Lowercase the string.
%
    text = lower ( strjust ( text, 'left' ) );
%
%  Locate the first word by finding the indices of the whitespace characters.
%
    blanks = find ( isspace ( text ) == 1 );

    if ( isempty ( blanks ) )
      word1 = text;
    else
      word1 = text(1:blanks(1)-1);
    end
%
%  Action depends on initial keyword.
%
    switch word1

       case 'solid'

         solid = solid + 1;

       case 'facet'

         face = face + 1;
         face_normal(:,face) = sscanf ( text, '%*s %*s %f %f %f' );

       case 'vertex'

         node = node + 1;
         node_xyz(:,node) = sscanf ( text, '%*s %f %f %f' );

       case 'endloop'

         face_node(:,face) = [ node-2; node-1; node ];

       otherwise

     end

  end
%
%  Close the file.
%
  fclose ( input_file_unit );

  return
end
function [ solid_num, node_num, face_num, text_num ] = stla_size_fast ( ...
  input_file_name )

%*****************************************************************************80
%
%% STLA_SIZE_FAST determines sizes associated with an STLA file.
%
%  Discussion:
%
%    This routine assumes that the file is a legal STLA file.
%
%    To perform checks on the file, call STLA_CHECK first.
%
%    Note that the counts for the number of nodes and edges are
%    overestimates, since presumably, most nodes will be defined several
%    times, once for each face they are part of, and most edges will
%    be defined twice.
%
%  Modified:
%
%    16 February 2007
%
%  Author:
%
%    Bill Seufzer,
%    NASA, Langley
%
%  Reference:
%
%    3D Systems, Inc,
%    Stereolithography Interface Specification,
%    October 1989.
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer SOLID_NUM, the number of solids defined.
%    Presumably, this is 1.
%
%    Output, integer NODE_NUM, the number of vertices defined.
%
%    Output, integer FACE_NUM, the number of faces defined.
%
%    Output, integer TEXT_NUM, the number of lines of text.
%
  solid_num = -1;
  node_num = -1;
  face_num = -1;
  text_num = -1;
%
%  If no file input, try to get one from the user.
%
  if ( nargin < 1 )
    input_file_name = input ( 'Enter the name of the ASCII STL file.' );
    if ( isempty ( input_file_name ) )
      return
    end
  end

  input_file_unit = fopen ( input_file_name, 'r' );

  if ( input_file_unit < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'STLA_SIZE_FAST - Fatal error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'STLA_SIZE_FAST - Fatal error!' )
    return
  end

  solid_num = 0;
  node_num = 0;
  face_num = 0;
  text_num = 0;
%
%  Read the text, one line at a time.
%
  while ( 1 )

    text = fgetl ( input_file_unit );
%
%  On end of input, FGETL returns a -1 instead of a character array.
%
    if ( ~ischar ( text ) )
      break
    end

    text_num = text_num + 1;
%
%  Left justify the string (no initial blanks).
%  Lowercase the string.
%
    text = lower ( strjust ( text, 'left' ) );
%
%  Locate the first word by finding the indices of the whitespace characters.
%
    blanks = find ( isspace ( text ) == 1 );
%
%  If there are no blanks, the first word is the whole string.
%  Otherwise, the first word is the string before the first blank.
%
    if ( isempty ( blanks ) )
      word1 = text;
    else
      word1 = text(1:blanks(1)-1);
    end
%
%  Count occurrences of "SOLID", "FACET" or "VERTEX".
%
    switch word1

      case 'solid'
        solid_num = solid_num + 1;
      case 'facet'
        face_num = face_num + 1;
      case 'vertex'
        node_num = node_num + 1;
      otherwise

    end

  end
%
%  Close the file.
%
  fclose ( input_file_unit );

  return
end
function timestamp ( )

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

  return
end
function tri_surface_write ( node_file_name, triangle_file_name, dim_num, ...
  node_num, order_num, triangle_num, node_xyz, triangle_node )

%*****************************************************************************80
%
%% TRI_SURFACE_WRITE writes graphics information to a pair of TRI_SURFACE files.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string NODE_FILE_NAME, the name of the node file.
%
%    Input, string TRIANGLE_FILE_NAME, the name of the triangle file.
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer NODE_NUM, the number of points.
%
%    Input, integer ORDER_NUM, the order of the triangles.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, real NODE_XYZ(DIM_NUM,NODE_NUM), the node coordinates.
%
%    Input, integer TRIANGLE_NODE(ORDER_NUM,TRIANGLE_NUM),
%    the nodes that form the triangles.
%
  r8mat_write0 ( node_file_name, dim_num, node_num, node_xyz );

  i4mat_write0 ( triangle_file_name, order_num, triangle_num, ...
    triangle_node );

  return
end


text:
classdef TestMultiply
    %TestMultiply

    methods (Static)
        function test_rbg_image
            img = cv.imread(fullfile(mexopencv.root(),'test','img001.jpg'), 'ReduceScale',2);
            img = img/2;

            out = cv.multiply(img, img);
            validateattributes(out, {class(img)}, {'size',size(img)});
            if mexopencv.require('images')
                expected = my_multiply(img, img);
                assert(isequal(out, expected));
            end

            out = cv.multiply(img, img, 'DType','uint16');
            validateattributes(out, {'uint16'}, {'size',size(img)});
            if mexopencv.require('images')
                expected = my_multiply(img, img, 'uint16');
                assert(isequal(out, expected));
            end
        end

        function test_2d_matrix
            A = uint8([10 35; 20 25]);
            B = uint8([10 10; 20 10]);

            out = cv.multiply(A, B);
            validateattributes(out, {'uint8'}, {'size',size(A)});
            expected = uint8([100 255; 255 250]);  % A .* B
            assert(isequal(out, expected));

            out = cv.multiply(A, 2);
            validateattributes(out, {'uint8'}, {'size',size(A)});
            expected = A .* 2;
            assert(isequal(out, expected));

            out = cv.multiply(2, A);
            validateattributes(out, {'uint8'}, {'size',size(A)});
            expected = 2 .* A;
            assert(isequal(out, expected));

            out = cv.multiply(2, 2);
            validateattributes(out, {'double'}, {'scalar'});
            expected = 2 * 2;
            assert(isequal(out, expected));

            out = cv.multiply(A, B, 'Scale',2);
            validateattributes(out, {'uint8'}, {'size',size(A)});
            expected = uint8([200 255; 255 255]);  % 2*A .* B
            assert(isequal(out, expected));
        end

        function test_output_depth
            out = cv.multiply(int16(20), uint8(20), 'DType','single');
            validateattributes(out, {'single'}, {'scalar'});
            assert(isequal(out, 400));
        end

        function test_saturation
            out = cv.multiply(uint8(20), uint8(30), 'DType','uint8');
            validateattributes(out, {'uint8'}, {'scalar'});
            assert(isequal(out, 255));

            out = cv.multiply(uint8(20), uint8(30), 'DType','uint16');
            validateattributes(out, {'uint16'}, {'scalar'});
            assert(isequal(out, 600));

            % saturation not applied for int32 (result overflows)
            out = cv.multiply(int32(50000), int32(50000), 'DType','int32');
            %validateattributes(out, {'int32'}, {'scalar', '<',0});

            % unlike MATLAB where result saturates
            %out = int32(50000) * int32(50000);
            %assert(isequal(out, intmax('int32')));
        end

        function test_error_argnum
            try
                cv.multiply();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

function out = my_multiply(src1, src2, dtype)
    %MY_MULTIPLY  Similar to cv.multiply using immultiply from IPT

    if nargin < 3, dtype = class(src1); end

    % multiply two images with specified output class
    out = immultiply(cast(src1, dtype), cast(src2, dtype));
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%% This file will be used to plot the different simulation data
% Coded By: K
% Last Updated: 2019/06/23
%%
clc;clear all; close all;

%% Add path
addpath('./Datas')

%% Load data
load('SimulationData_V_1.mat')

%% Plot
close all
% figure(1)
% waterfall(x,tspan,real(u))
% colormap([0 0 0]);
% view(42,55)
% set(gca,'FontSize',18);
% set(gcf,'Position',[100 100 600 400]);
% set(gcf,'PaperPositionMode','auto');
% grid on

figure(2)
Spacer1=10;
Spacer2=15;
h1=surf(x,tspan,real(u))
set(h1,'LineStyle','none')
view(42,55)
set(gca,'FontSize',40);
set(gcf,'Position',[100 100 600 400]);
set(gcf,'PaperPositionMode','auto');
grid on

hold on
h2=surf(x(1:Spacer1:end),tspan(1:Spacer2:end),real(u(1:Spacer2:length(tspan),1:Spacer1:length(x))))
set(h2,'LineStyle','-')

% figure(3)
% pcolor(x,tspan,real(u))
% shading interp
% set(gca,'FontSize',18);
% set(gcf,'Position',[100 100 600 400]);
% set(gcf,'PaperPositionMode','auto');
% grid on

author: dynamicslab

repo: SINDy-PI


text:
function title = p12_title ( option )

%*****************************************************************************80
%
%% P12_TITLE sets the title for problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, string TITLE, the title of the problem.
%
  if ( option == 1 )
    title = 'Materially nonlinear problem, NPOLYS = 2, NDERIV = 1.';
  elseif ( option == 2 )
    title = 'Materially nonlinear problem, NPOLYS = 4, NDERIV = 1.';
  elseif ( option == 3 )
    title = 'Materially nonlinear problem, NPOLYS = 4, NDERIV = 2.';
  elseif ( option == 4 )
    title = 'Materially nonlinear problem, NPOLYS = 6, NDERIV = 1.';
  elseif ( option == 5 )
    title = 'Materially nonlinear problem, NPOLYS = 6, NDERIV = 2.';
  elseif ( option == 6 )
    title = 'Materially nonlinear problem, NPOLYS = 6, NDERIV = 3.';
  end

  return
end

author: johannesgerer


text:
function [c,s,px,py,pz,rho,u,lambda] = unbalance3d_dose_2( source,target,gamma,H,reverse)

layer1=0.4;
layer2=1.0;

%profile on;
%% Input
%Source/Target

%load('deformedDoseOutputForLung_490pts_new/doseGrid.mat');
% Parameters
[nx,ny,nz]=size(source);
Hx=H.Hx;
Hy=H.Hy;
Hz=H.Hz;
% Hx=abs(xDoseV(length(xDoseV))-xDoseV(1))/10;
% Hy=abs(yDoseV(length(yDoseV))-yDoseV(1))/10;
% Hz=abs(zDoseV(length(zDoseV))-zDoseV(1))/10;

n=2;
nt    = 10;
hx    = Hx/nx;
hy    = Hy/ny;
hz    = Hz/nz;
ht    = 1/nt;

%% Input Normalization
source=double(source);
target=double(target);
if reverse=='r'
    source=255-source;
    target=255-target;
end
tr0  = sum(source(:));
tr1  = sum(target(:));
source = source+tr0/nx/ny/nz/15;
target = target+tr0/nx/ny/nz/15;

g=0;
g0=g+tr1/nx/ny/nz;
g1=g+tr0/nx/ny/nz;

rho0=g0*ones([size(source),2]);
rho0(:,:,:,1) = source;
rho1=g1*ones([size(target),2]);
rho1(:,:,:,1) = target;

%% Parameters Initialization
dim   = [n,nx,ny,nz,nt];
dim1=[n,round(nx*layer1),round(ny*layer1),round(nz*layer1),nt];
h     = [hx,hy,hz,ht];
h1=[Hx/dim1(2),Hy/dim1(3),Hz/dim1(4),ht];
param = paraminit(rho0(:),rho1(:),dim,h,gamma);
rho0_0 = fastResample3d(rho0,[dim1(2),dim1(3),dim1(4)]);
rho1_0 = fastResample3d(rho1,[dim1(2),dim1(3),dim1(4)]);
param1 = paraminit(rho0_0,rho1_0, dim1, h1, gamma);

m=param.m;
%% Set initial values for p rho u
px1     = zeros(n,dim1(2)-1,dim1(3),dim1(4),nt);
py1     = zeros(n,dim1(2),dim1(3)-1,dim1(4),nt);
pz1     = zeros(n,dim1(2),dim1(3),dim1(4)-1,nt);
u1      = zeros(m,dim1(2),dim1(3),dim1(4),nt);

tc=linspace(ht*3/2,1-ht/2,nt-1);
rho0_1 = fastResample3d(rho0,[dim1(2),dim1(3),dim1(4)]);
rho1_1 = fastResample3d(rho1,[dim1(2),dim1(3),dim1(4)]);
rho_1=ones(1,nt-1).*rho0_1(:)+tc.*(rho1_1(:)-rho0_1(:));

lambda1=zeros(dim1);
% initial values
init1.px     = px1;
init1.py     = py1;
init1.pz     = pz1;
init1.rho    = rho_1;
init1.u      = u1;
init1.lambda=lambda1;
%% iterations
max_iter=30;
tic
[px1,py1,pz1,rho_1,u1,lambda1,errorhist1] = SQPOPT(@(px,py,pz,rho,u,param) costgradf(px,py,pz,rho,u,param),init1,max_iter,param1);
toc
optcost=errorhist1(end,1);
iter_num=errorhist1(end,2);
fprintf('layer 1 Computed cost is: %f, iteration number=%d\n',optcost,iter_num);

%% results1
px1     = reshape(px1,[n,dim1(2)-1,dim1(3),dim1(4),nt]);
py1     = reshape(py1,[n,dim1(2),dim1(3)-1,dim1(4),nt]);
pz1     = reshape(pz1,[n,dim1(2),dim1(3),dim1(4)-1,nt]);
rho_1    = reshape(rho_1,[n,dim1(2),dim1(3),dim1(4),(nt-1)]);
u1      = reshape(u1,[m,dim1(2),dim1(3),dim1(4),nt]);
lambda1 = reshape(lambda1,[n,dim1(2),dim1(3),dim1(4),nt]);

px2     = zeros(n,dim(2)-1,dim(3),dim(4),nt);
py2     = zeros(n,dim(2),dim(3)-1,dim(4),nt);
pz2     = zeros(n,dim(2),dim(3),dim(4)-1,nt);
u2      = zeros(m,dim(2),dim(3),dim(4),nt);
rho_2    = zeros(n,dim(2),dim(3),dim(4),(nt-1));
lambda2=zeros(dim);
for k=1:nt
    temp=px1(:,:,:,:,k);
    temp=permute(temp,[2,3,4,1]);
    interp4M = fastResample3d(temp,[dim(2)-1,dim(3),dim(4)]);
    px2(:,:,:,:,k) = permute(interp4M,[4,1,2,3]);

    temp=py1(:,:,:,:,k);
    temp=permute(temp,[2,3,4,1]);
    interp4M = fastResample3d(temp,[dim(2),dim(3)-1,dim(4)]);
    py2(:,:,:,:,k) = permute(interp4M,[4,1,2,3]);

    temp=pz1(:,:,:,:,k);
    temp=permute(temp,[2,3,4,1]);
    interp4M = fastResample3d(temp,[dim(2),dim(3),dim(4)-1]);
    pz2(:,:,:,:,k) = permute(interp4M,[4,1,2,3]);

    if(k~=nt)
        temp=rho_1(:,:,:,:,k);
        temp=permute(temp,[2,3,4,1]);
        interp4M = fastResample3d(temp,[dim(2),dim(3),dim(4)]);
        rho_2(:,:,:,:,k) = permute(interp4M,[4,1,2,3]);
    end
    temp=u1(:,:,:,:,k);
    temp=permute(temp,[2,3,4,1]);
    interp4M = fastResample3d(temp,[dim(2),dim(3),dim(4)]);
    u2(:,:,:,:,k) = permute(interp4M,[4,1,2,3]);

    temp=lambda1(:,:,:,:,k);
    temp=permute(temp,[2,3,4,1]);
    interp4M = fastResample3d(temp,[dim(2),dim(3),dim(4)]);
    lambda2(:,:,:,:,k) = permute(interp4M,[4,1,2,3]);

end

% initial values
init.px     = px2;
init.py     = py2;
init.pz     = pz2;
init.rho    = rho_2;
init.u      = u2;
init.lambda=lambda2;
%% iterations
max_iter=20;
tic
[px,py,pz,rho,u,lambda,errorhist] = SQPOPT(@(px,py,pz,rho,u,param) costgradf(px,py,pz,rho,u,param),init,max_iter,param);
toc
optcost=errorhist(end,1);
iter_num=errorhist(end,2);
fprintf('layer 2 Computed cost is: %f, iteration number=%d\n',optcost,iter_num);

[c,s]=cs_results(px,py,pz,rho,u,param);
%% results
px     = reshape(px,[n,(nx-1),ny,nz,nt]);
py     = reshape(py,[n,nx,(ny-1),nz,nt]);
pz     = reshape(pz,[n,nx,ny,(nz-1),nt]);
rho    = reshape(rho,[n,nx,ny,nz,(nt-1)]);
u      = reshape(u,[m,nx*ny*nz*nt]);
lambda = reshape(lambda,[n,nx*ny*nz*nt]);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% CLASSIFY
% See also
%
% Clustering:
%   demoCluster        - Clustering demo.
%   demoGenData        - Generate data drawn form a mixture of Gaussians.
%   kmeans2            - Fast version of kmeans clustering.
%   meanShift          - meanShift clustering algorithm.
%   meanShiftIm        - Applies the meanShift algorithm to a joint spatial/range image.
%   meanShiftImExplore - Visualization to help choose sigmas for meanShiftIm.
%
% Calculating distances efficiently:
%   distMatrixShow     - Useful visualization of a distance matrix of clustered points.
%   pdist2             - Calculates the distance between sets of vectors.
%   softMin            - Calculates the softMin of a vector.
%
% Principal components analysis:
%   pca                - Principal components analysis (alternative to princomp).
%   pcaApply           - Companion function to pca.
%   pcaRandVec         - Generate random vectors in PCA subspace.
%   pcaVisualize       - Visualization of quality of approximation of X given principal comp.
%   visualizeData      - Project high dim. data unto principal components (PCA) for visualization.
%
% Confusion matrix display:
%   confMatrix         - Generates a confusion matrix according to true and predicted data labels.
%   confMatrixShow     - Used to display a confusion matrix.
%
% Radial Basis Functions (RBFs):
%   rbfComputeBasis    - Get locations and sizes of radial basis functions for use in rbf network.
%   rbfComputeFtrs     - Evaluate features of X given a set of radial basis functions.
%   rbfDemo            - Demonstration of rbf networks for regression.
%
% Fast random fern/forest classification/regression code:
%   fernsClfApply      - Apply learned fern classifier.
%   fernsClfTrain      - Train random fern classifier.
%   fernsInds          - Compute indices for each input by each fern.
%   fernsRegApply      - Apply learned fern regressor.
%   fernsRegTrain      - Train boosted fern regressor.
%   forestApply        - Apply learned forest classifier.
%   forestTrain        - Train random forest classifier.
%
% Fast boosted decision tree code:
%   adaBoostTrain      - Train boosted decision tree classifier.
%   adaBoostApply      - Apply learned boosted decision tree classifier.
%   binaryTreeTrain    - Train binary decision tree classifier.
%   binaryTreeApply    - Apply learned binary decision tree classifier.

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function fname_out = savecase(fname, varargin)
%SAVECASE  Saves a MATPOWER case file, given a filename and the data.
%   SAVECASE(FNAME, CASESTRUCT)
%   SAVECASE(FNAME, CASESTRUCT, VERSION)
%   SAVECASE(FNAME, BASEMVA, BUS, GEN, BRANCH)
%   SAVECASE(FNAME, BASEMVA, BUS, GEN, BRANCH, GENCOST)
%   SAVECASE(FNAME, BASEMVA, BUS, GEN, BRANCH, AREAS, GENCOST)
%   SAVECASE(FNAME, COMMENT, CASESTRUCT)
%   SAVECASE(FNAME, COMMENT, CASESTRUCT, VERSION)
%   SAVECASE(FNAME, COMMENT, BASEMVA, BUS, GEN, BRANCH)
%   SAVECASE(FNAME, COMMENT, BASEMVA, BUS, GEN, BRANCH, GENCOST)
%   SAVECASE(FNAME, COMMENT, BASEMVA, BUS, GEN, BRANCH, AREAS, GENCOST)
%
%   FNAME = SAVECASE(FNAME, ...)
%
%   Writes a MATPOWER case file, given a filename and data struct or list of
%   data matrices. The FNAME parameter is the name of the file to be created or
%   overwritten. If FNAME ends with '.mat' it saves the case as a MAT-file
%   otherwise it saves it as an M-file. Optionally returns the filename,
%   with extension added if necessary. The optional COMMENT argument is
%   either string (single line comment) or a cell array of strings which
%   are inserted as comments. When using a MATPOWER case struct, if the
%   optional VERSION argument is '1' it will modify the data matrices to
%   version 1 format before saving.

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

%% define named indices into bus, gen, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;

%% default arguments
if ischar(varargin{1}) || iscell(varargin{1})
    if ischar(varargin{1})
        comment = {varargin{1}};    %% convert char to single element cell array
    else
        comment = varargin{1};
    end
    [args{1:(length(varargin)-1)}] = deal(varargin{2:end});
else
    comment = {''};
    args = varargin;
end
mpc_ver = '2';              %% default MATPOWER case file version
if isstruct(args{1})        %% 1st real argument is a struct
    mpc = args{1};
    if length(args) > 1
        mpc.version = args{2};
        mpc_ver = mpc.version;
    end
    baseMVA = mpc.baseMVA;
    bus     = mpc.bus;
    gen     = mpc.gen;
    branch  = mpc.branch;
    if isfield(mpc, 'areas')
        areas   = mpc.areas;
    end
    if isfield(mpc, 'gencost')
        gencost = mpc.gencost;
    end
else                        %% 1st real argument is NOT a struct
    baseMVA = args{1};
    bus     = args{2};
    gen     = args{3};
    branch  = args{4};
    mpc.baseMVA = baseMVA;
    mpc.bus     = bus;
    mpc.gen     = gen;
    mpc.branch  = branch;
    if length(args) == 5
        gencost = args{5};
        mpc.gencost = gencost;
    end
    if length(args) == 6
        areas   = args{5};
        gencost = args{6};
        mpc.areas   = areas;
        mpc.gencost = gencost;
    end
end

%% modifications for version 1 format
if strcmp(mpc_ver, '1')
    %% remove extra columns of gen
    if size(gen, 2) >= MU_QMIN
        gen = gen(:, [1:PMIN, MU_PMAX:MU_QMIN]);
    else
        gen = gen(:, 1:PMIN);
    end
    %% use the version 1 values for column names
    shift = MU_PMAX - PMIN - 1;
    tmp = num2cell([MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] - shift);
    [MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = deal(tmp{:});

    %% remove extra columns of branch
    if size(branch, 2) >= MU_ST
        branch = branch(:, [1:BR_STATUS, PF:MU_ST]);
    elseif size(branch, 2) >= QT
        branch = branch(:, [1:BR_STATUS, PF:QT]);
    else
        branch = branch(:, 1:BR_STATUS);
    end
    %% use the version 1 values for column names
    shift = PF - BR_STATUS - 1;
    tmp = num2cell([PF, QF, PT, QT, MU_SF, MU_ST] - shift);
    [PF, QF, PT, QT, MU_SF, MU_ST] = deal(tmp{:});
end

%% verify valid filename
[pathstr, fcn_name, extension] = fileparts(fname);
if isempty(extension)
    extension = '.m';
end
if regexp(fcn_name, '\W')
    old_fcn_name = fcn_name;
    fcn_name = regexprep(fcn_name, '\W', '_');
    fprintf('WARNING: ''%s'' is not a valid function name, changed to ''%s''\n', old_fcn_name, fcn_name);
end
fname = fullfile(pathstr, [fcn_name extension]);

%% open and write the file
if strcmp(upper(extension), '.MAT')     %% MAT-file
    if strcmp(mpc_ver, '1')
        if exist('areas', 'var') && exist('gencost', 'var')
            save(fname, 'baseMVA', 'bus', 'gen', 'branch', 'areas', 'gencost');
        else
            save(fname, 'baseMVA', 'bus', 'gen', 'branch');
        end
    else
        save(fname, 'baseMVA', 'mpc');
    end
else                                %% M-file
    %% open file
    [fd, msg] = fopen(fname, 'wt');     %% print it to an M-file
    if fd == -1
        error(['savecase: ', msg]);
    end

    %% function header, etc.
    if strcmp(mpc_ver, '1')
        if exist('areas', 'var') && exist('gencost', 'var') && ~isempty(gencost)
            fprintf(fd, 'function [baseMVA, bus, gen, branch, areas, gencost] = %s\n', fcn_name);
        else
            fprintf(fd, 'function [baseMVA, bus, gen, branch] = %s\n', fcn_name);
        end
        prefix = '';
    else
        fprintf(fd, 'function mpc = %s\n', fcn_name);
        prefix = 'mpc.';
    end
    if isempty(comment{1})
        comment{1} = sprintf('%s', upper(fcn_name));
    else
        comment{1} = sprintf('%s  %s', upper(fcn_name), comment{1});
    end
    for k = 1:length(comment)
        fprintf(fd, '%%%s\n', comment{k});
    end
    fprintf(fd, '\n%%%% MATPOWER Case Format : Version %s\n', mpc_ver);
    if ~strcmp(mpc_ver, '1')
        fprintf(fd, 'mpc.version = ''%s'';\n', mpc_ver);
    end
    fprintf(fd, '\n%%%%-----  Power Flow Data  -----%%%%\n');
    fprintf(fd, '%%%% system MVA base\n');
    fprintf(fd, '%sbaseMVA = %.9g;\n', prefix, baseMVA);

    %% bus data
    ncols = size(bus, 2);
    fprintf(fd, '\n%%%% bus data\n');
    fprintf(fd, '%%\tbus_i\ttype\tPd\tQd\tGs\tBs\tarea\tVm\tVa\tbaseKV\tzone\tVmax\tVmin');
    if ncols >= MU_VMIN             %% opf SOLVED, save with lambda's & mu's
        fprintf(fd, '\tlam_P\tlam_Q\tmu_Vmax\tmu_Vmin');
    end
    fprintf(fd, '\n%sbus = [\n', prefix);
    if ncols < MU_VMIN              %% opf NOT SOLVED, save without lambda's & mu's
        fprintf(fd, '\t%d\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%d\t%.9g\t%.9g\t%.9g\t%d\t%.9g\t%.9g;\n', bus(:, 1:VMIN).');
    else                            %% opf SOLVED, save with lambda's & mu's
        fprintf(fd, '\t%d\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%d\t%.9g\t%.9g\t%.9g\t%d\t%.9g\t%.9g\t%.4f\t%.4f\t%.4f\t%.4f;\n', bus(:, 1:MU_VMIN).');
    end
    fprintf(fd, '];\n');

    %% generator data
    ncols = size(gen, 2);
    fprintf(fd, '\n%%%% generator data\n');
    fprintf(fd, '%%\tbus\tPg\tQg\tQmax\tQmin\tVg\tmBase\tstatus\tPmax\tPmin');
    if ~strcmp(mpc_ver, '1')
        fprintf(fd, '\tPc1\tPc2\tQc1min\tQc1max\tQc2min\tQc2max\tramp_agc\tramp_10\tramp_30\tramp_q\tapf');
    end
    if ncols >= MU_QMIN             %% opf SOLVED, save with mu's
        fprintf(fd, '\tmu_Pmax\tmu_Pmin\tmu_Qmax\tmu_Qmin');
    end
    fprintf(fd, '\n%sgen = [\n', prefix);
    if ncols < MU_QMIN              %% opf NOT SOLVED, save without mu's
        if strcmp(mpc_ver, '1')
            fprintf(fd, '\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%d\t%.9g\t%.9g;\n', gen(:, 1:PMIN).');
        else
            fprintf(fd, '\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g;\n', gen(:, 1:APF).');
        end
    else
        if strcmp(mpc_ver, '1')
            fprintf(fd, '\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%d\t%.9g\t%.9g\t%.4f\t%.4f\t%.4f\t%.4f;\n', gen(:, 1:MU_QMIN).');
        else
            fprintf(fd, '\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.4f\t%.4f\t%.4f\t%.4f;\n', gen(:, 1:MU_QMIN).');
        end
    end
    fprintf(fd, '];\n');

    %% branch data
    ncols = size(branch, 2);
    fprintf(fd, '\n%%%% branch data\n');
    fprintf(fd, '%%\tfbus\ttbus\tr\tx\tb\trateA\trateB\trateC\tratio\tangle\tstatus');
    if ~strcmp(mpc_ver, '1')
        fprintf(fd, '\tangmin\tangmax');
    end
    if ncols >= QT                  %% power flow SOLVED, save with line flows
        fprintf(fd, '\tPf\tQf\tPt\tQt');
    end
    if ncols >= MU_ST               %% opf SOLVED, save with mu's
        fprintf(fd, '\tmu_Sf\tmu_St');
        if ~strcmp(mpc_ver, '1')
            fprintf(fd, '\tmu_angmin\tmu_angmax');
        end
    end
    fprintf(fd, '\n%sbranch = [\n', prefix);
    if ncols < QT                   %% power flow NOT SOLVED, save without line flows or mu's
        if strcmp(mpc_ver, '1')
            fprintf(fd, '\t%d\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%d;\n', branch(:, 1:BR_STATUS).');
        else
            fprintf(fd, '\t%d\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%d\t%.9g\t%.9g;\n', branch(:, 1:ANGMAX).');
        end
    elseif ncols < MU_ST            %% power flow SOLVED, save with line flows but without mu's
        if strcmp(mpc_ver, '1')
            fprintf(fd, '\t%d\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%d\t%.4f\t%.4f\t%.4f\t%.4f;\n', branch(:, 1:QT).');
        else
            fprintf(fd, '\t%d\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%d\t%.9g\t%.9g\t%.4f\t%.4f\t%.4f\t%.4f;\n', branch(:, 1:QT).');
       end
    else                            %% opf SOLVED, save with lineflows & mu's
        if strcmp(mpc_ver, '1')
            fprintf(fd, '\t%d\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%d\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f;\n', branch(:, 1:MU_ST).');
        else
            fprintf(fd, '\t%d\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%d\t%.9g\t%.9g\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f;\n', branch(:, 1:MU_ANGMAX).');
        end
    end
    fprintf(fd, '];\n');

    %% OPF data
    if (exist('areas', 'var') && ~isempty(areas)) || ...
        (exist('gencost', 'var') && ~isempty(gencost))
        fprintf(fd, '\n%%%%-----  OPF Data  -----%%%%');
    end
    if exist('areas', 'var') && ~isempty(areas)
        %% area data
        fprintf(fd, '\n%%%% area data\n');
        fprintf(fd, '%%\tarea\trefbus\n');
        fprintf(fd, '%sareas = [\n', prefix);
        if ~isempty(areas)
            fprintf(fd, '\t%d\t%d;\n', areas(:, 1:2).');
        end
        fprintf(fd, '];\n');
    end
    if exist('gencost', 'var') && ~isempty(gencost)
        %% generator cost data
        fprintf(fd, '\n%%%% generator cost data\n');
        fprintf(fd, '%%\t1\tstartup\tshutdown\tn\tx1\ty1\t...\txn\tyn\n');
        fprintf(fd, '%%\t2\tstartup\tshutdown\tn\tc(n-1)\t...\tc0\n');
        fprintf(fd, '%sgencost = [\n', prefix);
        if ~isempty(gencost)
            n1 = 2 * max(gencost(gencost(:, MODEL) == PW_LINEAR,  NCOST));
            n2 =     max(gencost(gencost(:, MODEL) == POLYNOMIAL, NCOST));
            n = max([n1; n2]);
            if size(gencost, 2) < n + 4
                error('savecase: gencost data claims it has more columns than it does');
            end
            template = '\t%d\t%.9g\t%.9g\t%d';
            for i = 1:n
                template = [template, '\t%.9g'];
            end
            template = [template, ';\n'];
            fprintf(fd, template, gencost(:, 1:n+4).');
        end
        fprintf(fd, '];\n');
    end

    if ~strcmp(mpc_ver, '1')
        %% generalized OPF user data
        if (isfield(mpc, 'A') && ~isempty(mpc.A)) || ...
                (isfield(mpc, 'N') && ~isempty(mpc.N))
            fprintf(fd, '\n%%%%-----  Generalized OPF User Data  -----%%%%');
        end

        %% user constraints
        if isfield(mpc, 'A') && ~isempty(mpc.A)
            %% A
            fprintf(fd, '\n%%%% user constraints\n');
            print_sparse(fd, sprintf('%sA', prefix), mpc.A);
            if isfield(mpc, 'l') && ~isempty(mpc.l) && ...
                    isfield(mpc, 'u') && ~isempty(mpc.u)
                fprintf(fd, 'lbub = [\n');
                fprintf(fd, '\t%.9g\t%.9g;\n', [mpc.l mpc.u].');
                fprintf(fd, '];\n');
                fprintf(fd, '%sl = lbub(:, 1);\n', prefix);
                fprintf(fd, '%su = lbub(:, 2);\n\n', prefix);
            elseif isfield(mpc, 'l') && ~isempty(mpc.l)
                fprintf(fd, '%sl = [\n', prefix);
                fprintf(fd, '\t%.9g;\n', mpc.l);
                fprintf(fd, '];\n\n');
            elseif isfield(mpc, 'u') && ~isempty(mpc.u)
                fprintf(fd, '%su = [\n', prefix);
                fprintf(fd, '\t%.9g;\n', mpc.u);
                fprintf(fd, '];\n');
            end
        end

        %% user costs
        if isfield(mpc, 'N') && ~isempty(mpc.N)
            fprintf(fd, '\n%%%% user costs\n');
            print_sparse(fd, sprintf('%sN', prefix), mpc.N);
            if isfield(mpc, 'H') && ~isempty(mpc.H)
                print_sparse(fd, sprintf('%sH', prefix), mpc.H);
            end
            if isfield(mpc, 'fparm') && ~isempty(mpc.fparm)
                fprintf(fd, 'Cw_fparm = [\n');
                fprintf(fd, '\t%.9g\t%d\t%.9g\t%.9g\t%.9g;\n', [mpc.Cw mpc.fparm].');
                fprintf(fd, '];\n');
                fprintf(fd, '%sCw    = Cw_fparm(:, 1);\n', prefix);
                fprintf(fd, '%sfparm = Cw_fparm(:, 2:5);\n', prefix);
            else
                fprintf(fd, '%sCw = [\n', prefix);
                fprintf(fd, '\t%.9g;\n', mpc.Cw);
                fprintf(fd, '];\n');
            end
        end

        %% user vars
        if isfield(mpc, 'z0') || isfield(mpc, 'zl') || isfield(mpc, 'zu')
            fprintf(fd, '\n%%%% user vars\n');
        end
        if isfield(mpc, 'z0') && ~isempty(mpc.z0)
            fprintf(fd, '%sz0 = [\n', prefix);
            fprintf(fd, '\t%.9g;\n', mpc.z0);
            fprintf(fd, '];\n');
        end
        if isfield(mpc, 'zl') && ~isempty(mpc.zl)
            fprintf(fd, '%szl = [\n', prefix);
            fprintf(fd, '\t%.9g;\n', mpc.zl);
            fprintf(fd, '];\n');
        end
        if isfield(mpc, 'zu') && ~isempty(mpc.zu)
            fprintf(fd, '%szu = [\n', prefix);
            fprintf(fd, '\t%.9g;\n', mpc.zu);
            fprintf(fd, '];\n');
        end

        %% (optional) generator unit types
        if isfield(mpc, 'gentype') && iscell(mpc.gentype)
            ng = length(mpc.gentype);
            if size(mpc.gen, 1) ~= ng
                warning('savecase: gentype field does not have the expected dimensions (length = %d, expected %d)', ng, size(mpc.gen, 1));
            end

            fprintf(fd, '\n%%%% generator unit type (see GENTYPES)\n');
            fprintf(fd, '%sgentype = {\n', prefix);
            for k = 1:ng
                fprintf(fd, '\t''%s'';\n', mpc.gentype{k});
            end
            fprintf(fd, '};\n');
        end

        %% (optional) generator fuel types
        if isfield(mpc, 'genfuel') && iscell(mpc.genfuel)
            ng = length(mpc.genfuel);
            if size(mpc.gen, 1) ~= ng
                warning('savecase: genfuel field does not have the expected dimensions (length = %d, expected %d)', ng, size(mpc.gen, 1));
            end

            fprintf(fd, '\n%%%% generator fuel type (see GENFUELS)\n');
            fprintf(fd, '%sgenfuel = {\n', prefix);
            for k = 1:ng
                fprintf(fd, '\t''%s'';\n', mpc.genfuel{k});
            end
            fprintf(fd, '};\n');
        end

        %% (optional) bus names
        if isfield(mpc, 'bus_name') && iscell(mpc.bus_name)
            nb = length(mpc.bus_name);
            if size(mpc.bus, 1) ~= nb
                warning('savecase: bus_name field does not have the expected dimensions (length = %d, expected %d)', nb, size(mpc.bus, 1));
            end

            fprintf(fd, '\n%%%% bus names\n');
            fprintf(fd, '%sbus_name = {\n', prefix);
            for k = 1:nb
                fprintf(fd, '\t''%s'';\n', strrep(mpc.bus_name{k}, '''', ''''''));
            end
            fprintf(fd, '};\n');
        end

        %% execute userfcn callbacks for 'savecase' stage
        if isfield(mpc, 'userfcn')
            run_userfcn(mpc.userfcn, 'savecase', mpc, fd, prefix);
        end
    end

    %% close file
    if fd ~= 1
        fclose(fd);
    end
end

if nargout > 0
    fname_out = fname;
end

function print_sparse(fd, varname, A)

[i, j, s] = find(A);
[m, n] = size(A);

if isempty(s)
    fprintf(fd, '%s = sparse(%d, %d);\n', varname, m, n);
else
    fprintf(fd, 'ijs = [\n');
    if m == 1           %% i, j, s are row vectors
        fprintf(fd, '\t%d\t%d\t%.9g;\n', [i; j; s]);
    else                %% i, j, s are column vectors
        fprintf(fd, '\t%d\t%d\t%.9g;\n', [i j s].');
    end
    fprintf(fd, '];\n');
    fprintf(fd, '%s = sparse(ijs(:, 1), ijs(:, 2), ijs(:, 3), %d, %d);\n', varname, m, n);
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function savewav(y,filename,fs)
%SAVEWAV saves floating point audio data in a *.wav with 32bit precision
%
%   Usage: savewav(y,filename,fs)
%
%   Input parameters:
%       y         - audio data in floating point [NFrames x NChannels]
%       filename  - name of *.wav file (file extension will NOT be added
%                   automatically)
%       fs        - sample rate / Hz
%
%   Since wavwrite has been removed in MATLAB 2015b and audiowrite does not
%   support audio data with a channel number higher than 256, this is an
%   alternative to save massive multichannel. This function does only support
%   floating point data and saves it 32bit floating point precision.
%
%   See also: audiowrite, audioread, fwrite

% Copyright (c) 2015 Fiete Winter
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.

%% Input Check
narginchk(3,3);
if ~isfloat(y)
  error('%s: only floating point supported', upper(mfilename));
end
if ~ndims(y) > 2
  error('%s: only audio data with [NFrames x NChannels] supported', ...
    upper(mfilename));
end

%%
[NFrames, NChannels] = size(y);  % Number of Frames and Channels
NSamples = NFrames*NChannels;  % total Samples
Nbits = 32;  % 32 floating point
Nbyte = 4;  % bytes per Samples

fileID = fopen(filename,'w', 'ieee-le');  % open wav-file (little Endian)

fwrite(fileID, 'RIFF', 	'char*1');  % RIFF header
% number of bytes following (chunk size)
Noverallbytes = 4 + 4; % for 'WAVE' and 'fmt '
Noverallbytes = Noverallbytes + 4;  % for the number of bytes in 'fmt '
Noverallbytes = Noverallbytes + 16; % bytes in 'fmt '
Noverallbytes = Noverallbytes + 12; % for 'fact'
Noverallbytes = Noverallbytes + 8 + Nbyte*NSamples;  % 'data' chunk
fwrite(fileID, Noverallbytes, 'uint32');

fwrite(fileID, 'WAVE', 	'char*1');  % WAVE header

%% fmt chunk
fwrite(fileID, 'fmt ', 	'char*1');  % mind the whitespace!
fwrite(fileID, 16, 	'uint32');  % number of bytes following (chunk size)
fwrite(fileID, 3, 'uint16');  % data type (3 for float)
fwrite(fileID, NChannels, 'uint16');  % number of channels
fwrite(fileID, fs, 'uint32');  % samplerate
fwrite(fileID, fs * NChannels * Nbyte, 'uint32');  % bytes per second
fwrite(fileID, NChannels * Nbyte, 'uint16');  % frame size
fwrite(fileID, Nbits, 'uint16');  % bits per sample

%% fact chunk (floating point only)
fwrite(fileID, 'fact', 	'char*1');
fwrite(fileID, 4, 	'uint32');  % number of bytes following (chunk size)
fwrite(fileID, NFrames, 	'uint32');  % number of frames

%% data chunk
fwrite(fileID, 'data', 	'char*1');
fwrite(fileID, Nbyte*NSamples, 'uint32');  % number of bytes following (chunk size)
y = y.';  % prepare data for interlacing
fwrite(fileID, y(:), 'float32');  % interlaced data

%%
fclose(fileID);

end

author: sfstoolbox

repo: sfs-matlab


text:
function Subs = som_ind2sub(msize,inds)

%SOM_IND2SUB Map grid subscripts from linear index.
%
% Subs = som_ind2sub(msize,inds)
%
%  sub = som_ind2sub([10 15],44);
%  sub = som_ind2sub(sMap,44);
%  sub = som_ind2sub(sMap.msize,44);
%  Subs = som_ind2sub([10 15],[44 13 91]');
%
%  Input and output arguments:
%   msize  (struct) map or topology struct
%          (vector) size 1 x m, specifies the map grid size
%   inds   (vector) size n x 1, linear indeces of n map units
%
%   Subs   (matrix) size n x m, the subscripts
%
% See also SOM_SUB2IND.

% Contributed to SOM Toolbox vs2, February 2nd, 2000 by Juha Vesanto
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta juuso 300798

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isstruct(msize),
  if strcmp(msize.type,'som_map'), msize = msize.topol.msize;
  elseif strcmp(msize.type,'som_topol'), msize = msize.msize;
  else error('Invalid first argument.'); end
end

n = length(msize);
k = [1 cumprod(msize(1:end-1))];
inds = inds - 1;
for i = n:-1:1,
  Subs(:,i) = floor(inds/k(i))+1;
  inds = rem(inds,k(i));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: ilarinieminen

repo: SOM-Toolbox


text:
function [F] = mci_nmm_r2p2_dfdp (x,u,P,M)
% State Jacobian for two region, two parameter NMM
% FORMAT [F] = mci_nmm_r2p2_dfdp (x,u,P,M)
%
% x         State
% u         Inputs
% P         Parameters
% M         Model structure
%
% F         F(i,j) = df(x)_i/dtheta_j
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny and Biswa Sengupta
% $Id: mci_nmm_r2p2_dfdx.m 6548 2015-09-11 12:39:47Z will $

% 18 state variables
F = zeros(18,18);

% f         Flow, dx/dt

curr_P=M.can_P; % Canonical parameter set

% 2 free parameters
curr_P.A{1}(2,1)=P(1); % Forward connection, w_21
curr_P.A{2}(1,2)=P(2); % Backward connection, w_12

f = mci_nmm_fx_delay(x,u,curr_P,M);

author: spm


text:
function [km2] = cm22km2(cm2)
% Convert area from square centimeters to square kilometers.
% Chad A. Greene 2012
km2 = cm2*1E-10;

author: Sable

repo: mcbench-benchmarks


text:
function ell = equivalentEllipsoid(points)
% Equivalent ellipsoid of a set of 3D points.
%
%   ELL = equivalentEllipsoid(PTS)
%   Compute the equivalent ellipsoid of the set of points PTS. The result
%   is an ellipsoid defined by:
%   ELL = [XC YC ZC  A B C  PHI THETA PSI]
%   where [XC YC ZY] is the center, [A B C] are the lengths of the
%   semi-axes (in decreasing order), and [PHI THETA PSI] are Euler angles
%   representing  the ellipsoid orientation, in degrees.
%
%   Example
%     pts = randn(300, 3);
%     pts = transformPoint3d(pts, createScaling3d([6 4 2]));
%     pts = transformPoint3d(pts, createRotationOx(pi/6));
%     pts = transformPoint3d(pts, createRotationOy(pi/4));
%     pts = transformPoint3d(pts, createRotationOz(pi/3));
%     pts = transformPoint3d(pts, createTranslation3d([5 4 3]));
%     elli = equivalentEllipsoid(pts);
%     figure; drawPoint3d(pts); axis equal;
%     hold on; drawEllipsoid(elli, ...
%         'drawEllipses', true, 'EllipseColor', 'b', 'EllipseWidth', 3);
%
%   See also
%     spheres, drawEllipsoid, equivalentEllipse, principalAxes
%     rotation3dToEulerAngles

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2011-03-12,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% number of points
n = size(points, 1);

% compute centroid
center = mean(points);

% compute the covariance matrix
covPts = cov(points)/n;

% perform a principal component analysis with 2 variables,
% to extract equivalent axes
[U, S] = svd(covPts);

% extract length of each semi axis
radii = sqrt(5) * sqrt(diag(S)*n)';

% sort axes from greater to lower
[radii, ind] = sort(radii, 'descend');

% format U to ensure first axis points to positive x direction
U = U(ind, :);
if U(1,1) < 0
    U = -U;
    % keep matrix determinant positive
    U(:,3) = -U(:,3);
end

% convert axes rotation matrix to Euler angles
angles = rotation3dToEulerAngles(U);

% concatenate result to form an ellipsoid object
ell = [center, radii, angles];

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%%********************************************************************
%% infeaspt: generate an initial point for sdp.m
%%
%%  [X0,y0,Z0] = infeaspt(blk,At,C,b,options,scalefac);
%%
%%  options = 1  if want X0,Z0 to be scaled identity matrices
%%          = 2  if want X0,Z0 to be scalefac*(identity matrices).
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%********************************************************************

   function [X0,y0,Z0] = infeaspt(blk,At,C,b,options,scalefac);
%%
   if (nargin < 5); options = 1; end;
   if (options == 1); scalefac = []; end;
   if (options == 2) & (nargin < 6); scalefac = 1000; end;
   if (scalefac <= 0); error('scalefac must a positive number'); end;
   state = rand('state');
   rand('state',0);
%%
   if ~iscell(At); At = {At}; end;
   if ~iscell(C);  C = {C}; end;
   m = length(b);
   if all(size(At) == [size(blk,1) m]);
      convertyes = zeros(size(blk,1),1);
      for p = 1:size(blk,1)
         if strcmp(blk{p,1},'s') & all(size(At{p,1}) == sum(blk{p,2}))
            convertyes(p) = 1;
         end
      end
      if any(convertyes)
         At = svec(blk,At,ones(size(blk,1),1));
      end
   end;
%%
   %%[blk,At,C,b] = validate(blk,At,C,b);
%%
   X0 = cell(size(C)); Z0 = cell(size(C));
   m = length(b);
   for p = 1:size(blk,1);
      pblk = blk(p,:);
      blktmp = pblk{2};
      n = length(C{p});
      y0 = zeros(m,1);
      b2 = 1 + abs(b');
      if (options == 1);
         if strcmp(pblk{1},'s');
            normAni = [];
            X0{p} = sparse(n,n); Z0{p} = sparse(n,n);
            ss = [0, cumsum(blktmp)];
            tt = [0, cumsum(blktmp.*(blktmp+1)/2)];
            for i = 1:length(pblk{2})
               if ~isempty(At{p,1})
                  pos = [tt(i)+1 : tt(i+1)];
                  Ai = At{p,1}(pos,:);
                  normAni = 1+sqrt(sum(Ai.*Ai));
               end
               if (length(At(p,:)) >= 2) %% for low rank constraints
                  dd = At{p,3};
                  qq = [0, cumsum(pblk{3})]; normtmp = ones(1,length(pblk{3}));
                  idxD = [0; find(diff(dd(:,1))); size(dd,1)];
                  for k = 1:length(pblk{3})
                      idx = [qq(k)+1 : qq(k+1)];
                      idx2 = [idxD(k)+1: idxD(k+1)];
                      Ak = At{p,2}(:,idx);
                      ii = dd(idx2,2)-qq(k); %% undo cumulative indexing
                      jj = dd(idx2,3)-qq(k);
                      len = pblk{3}(k);
                      Dk = spconvert([ii,jj,dd(idx2,4); len,len,0]);
                      tmp = Ak'*Ak*Dk;
                      normtmp(1,k) = 1+sqrt(sum(sum(tmp.*tmp')));
                  end
                  normAni = [normAni, normtmp];
               end
               pos = [ss(i)+1 : ss(i+1)];  ni = length(pos);
               tmp = C{p}(pos,pos);
               normCni = 1+sqrt(sum(sum(tmp.*tmp)));
               const  = 10; %%--- old: const = 1;
               constX = max([const,sqrt(ni),ni*(b2./normAni)]);
               constZ = max([const,sqrt(ni),normAni,normCni]);
               X0{p}(pos,pos) = constX*spdiags(1+1e-10*rand(ni,1),0,ni,ni);
               Z0{p}(pos,pos) = constZ*spdiags(1+1e-10*rand(ni,1),0,ni,ni);
            end
         elseif strcmp(pblk{1},'q');
            s = 1+[0, cumsum(blktmp)];
            len = length(blktmp);
            normC = 1+norm(C{p});
            normA = 1+sqrt(sum(At{p,1}.*At{p,1}));
            idenqX = zeros(sum(blktmp),1);
            idenqZ = zeros(sum(blktmp),1);
            idenqX(s(1:len)) = max([1,b2./normA])*sqrt(blktmp') ;
            idenqZ(s(1:len)) = max([sqrt(blktmp); max([normA,normC])*ones(1,len)])';
            idenqX(s(1:len)) = idenqX(s(1:len)).*(1+1e-10*rand(len,1));
            idenqZ(s(1:len)) = idenqZ(s(1:len)).*(1+1e-10*rand(len,1));
            X0{p} = idenqX;
            Z0{p} = idenqZ;
         elseif strcmp(pblk{1},'l');
            normC = 1+norm(C{p});
            normA = 1+sqrt(sum(At{p,1}.*At{p,1}));
            const = 10; %%--- old: const =1;
            constX = max([const,sqrt(n),sqrt(n)*b2./normA]);
            constZ = max([const,sqrt(n),normA,normC]);
            X0{p} = constX*(1+1e-10*rand(n,1));
            Z0{p} = constZ*(1+1e-10*rand(n,1));
         elseif strcmp(pblk{1},'u');
            X0{p} = sparse(n,1);
            Z0{p} = sparse(n,1);
         else
            error(' blk: some fields not specified correctly');
         end;
      elseif (options == 2);
         if strcmp(pblk{1},'s');
            n = sum(blktmp);
            X0{p} = scalefac*spdiags(1+1e-10*rand(n,1),0,n,n);
            Z0{p} = scalefac*spdiags(1+1e-10*rand(n,1),0,n,n);
         elseif strcmp(pblk{1},'q');
            s = 1+[0, cumsum(blktmp)];
            len = length(blktmp);
            idenq = zeros(sum(blktmp),1);
            idenq(s(1:len)) = 1+1e-10*rand(len,1);
            X0{p} = scalefac*idenq;
            Z0{p} = scalefac*idenq;
         elseif strcmp(pblk{1},'l');
            X0{p} = scalefac*(1+1e-10*rand(n,1));
            Z0{p} = scalefac*(1+1e-10*rand(n,1));
         elseif strcmp(pblk{1},'u');
            X0{p} = sparse(n,1);
            Z0{p} = sparse(n,1);
         else
            error(' blk: some fields not specified correctly');
         end
      end
   end
   rand('state',state);
%%********************************************************************

repo: REST


text:
function up = predictor(xp,x,u)

up = interp1(x,u,xp,'spline');