text:
classdef gen80211aTestWaveform
    properties (Constant)
        %%
        fs  = 20e6 % sample rate
    end

    properties (Access = public)
        %%
        % carrier frequency
        fc
        % number of frames
        numPackets = 50
        % IEEE 802.11a Modulation and Coding Scheme [0-7]
        mcs = 2
        % Fix the WLAN frame length to a constant value? (true/false)
        constPktLen = true
        % Simulation seed
        % 0 - generate a random seed
        % xyz - seed is set to xyz
        seed = 0
    end

    properties (Access = public) % derived properties
        %%
        SampleTime_rx
        random_delay
        txPSDU
        txMPDU
        packetLen
        txWaveform
        MSDULenOctets
        MSDULenBits
        nht
    end

    properties (Access = public)
        %%
        noiseVar = 1e-6
        gain_per_packet
        rxWaveform
    end

    methods
        function obj = gen80211aTestWaveform(sim_settings, wlan_settings)
            mustBeMember(wlan_settings.fc,5180e6:20e6:5320e6);
            obj.fc = wlan_settings.fc;
            mustBeMember(wlan_settings.mcs,0:7);
            obj.mcs = wlan_settings.mcs;
            obj.constPktLen = wlan_settings.constPktLen;
            obj.seed = wlan_settings.seed;
            obj.numPackets = wlan_settings.numPackets;
            obj.SampleTime_rx = 1/obj.fs;
            obj.snr = sim_settings.snr;

            % Generate Non-HT Waveform
            obj = obj.genTestWaveform();
            if (sim_settings.SIM_STUDY == true)
                % apply channel impairments
                obj = obj.applyImpairments(sim_settings);
            end
        end
    end

    methods (Access = private)
        function obj = genTestWaveform(obj)
            % Generate 20 MHz Non-HT Waveform
            if (obj.seed == 0)
                obj.seed = int32(randi(intmax));
            end

            % If the data payload is of random length,
            % the minimum MSDU length is set to 250.
            % Otherwise, it is set to 2304.
            rng(obj.seed);
            if (obj.constPktLen)
                obj.MSDULenOctets = 2304*ones(obj.numPackets, 1);
            else
                obj.MSDULenOctets = randi([250 2304], obj.numPackets, 1);
            end
            rng(obj.seed);
            payload = randi([0 255], sum(obj.MSDULenOctets), 1);

            macConfig = wlanMACFrameConfig;
            macConfig.FrameType = 'Data';
            prev = 0;
            for ii = 1:obj.numPackets
                obj.nht{ii} = wlanNonHTConfig('MCS', obj.mcs, ...
                    'PSDULength', obj.MSDULenOctets(ii));

                macConfig.SequenceNumber = ii;
                [obj.txMPDU{ii}, lengthMPDU] = wlanMACFrame(payload(prev+1:prev+obj.MSDULenOctets(ii)), macConfig);
                prev = prev+obj.MSDULenOctets(ii);
                obj.nht{ii}.PSDULength = lengthMPDU;
                obj.txPSDU{ii} = reshape(de2bi(hex2dec(obj.txMPDU{ii}), 8)', [], 1);

                % WLAN waveform
                obj.txWaveform{ii} = wlanWaveformGenerator(obj.txPSDU{ii}, obj.nht{ii}, ...
                    'NumPackets',1,'IdleTime',20e-6).';
                obj.packetLen(ii) = length(obj.txWaveform{ii});
            end
            obj.MSDULenBits = (obj.MSDULenOctets+28)*8; % take overhead into account
            obj.txWaveform = cell2mat(obj.txWaveform).';
            obj.txWaveform = obj.txWaveform(:);

            % Start of the transmission is random
            obj.random_delay = round(2^17*(0.15+0.1*rand()));
            prepend_sig = zeros(obj.random_delay, 1);
            obj.txWaveform = [prepend_sig; obj.txWaveform];
        end

        function obj = applyImpairments(obj, sim_settings)
            % AWGN
            rng(obj.seed);
            awgnChan = comm.AWGNChannel('SignalPower',var(obj.txWaveform));
            awgnChan.NoiseMethod = 'Signal to noise ratio (SNR)';
            awgnChan.SNR = obj.snr;
            obj.rxWaveform = awgnChan(obj.txWaveform);
            % frequency/phase/timing offsets
            pfo = comm.PhaseFrequencyOffset('SampleRate', obj.fs, 'FrequencyOffset', 500, 'PhaseOffset', 30);
            obj.rxWaveform = pfo(obj.rxWaveform);

            % add channel impairments based on SIM_MODE
            if ((sim_settings.SIM_MODE == 2) || (sim_settings.SIM_MODE == 3))
                % Pass Non-HT waveform through a SISO channel.
                % The free-space path loss simulated is over a
                % (a) transmitter-to-receiver separation distance of 3 m,
                % (b) maximum Doppler shift of 3 Hz and
                % (c) an RMS path delay equal to 2x the sample time.
                dist = 30;
                pathLoss = 10^(-log10(4*pi*dist*(obj.fc/physconst('LightSpeed'))))*1e3;
                trms = 2/obj.fs;
                maxDoppShift = 3;
                rng(obj.seed);
                ch802 = comm.RayleighChannel('SampleRate',obj.fs,'MaximumDopplerShift',maxDoppShift,'PathDelays',trms);
                obj.rxWaveform = ch802(obj.rxWaveform)*pathLoss;
            end

            % arbitrary gain per packet
            if (sim_settings.GAIN_MODE == 0)
                obj.gain_per_packet = ones(obj.numPackets, 1);
                obj.rxWaveform = obj.txWaveform/150; % an arbitrary value to scale the entire waveform
            elseif (sim_settings.GAIN_MODE == 1)
                % Interpolated Gaussian random values as a model of gain
                % evolution. Values are scaled based on default AD9361 gain
                % table mappings.
                tmp_len = ceil(obj.numPackets/50)*50;
                tmp_g = randn(1,tmp_len/5);
                tmp_g = interp(tmp_g, 5);
                tmp_g = tmp_g-min(tmp_g);
                tmp_g = tmp_g/max(tmp_g);
                tmp_g = round(tmp_g*65);
                g = tmp_g(1:obj.numPackets);

                % Gain per packet in linear scale.
                obj.gain_per_packet = 10.^(g./20);
            end
            prev = obj.random_delay;
            for ii = 1:obj.numPackets
                obj.rxWaveform(prev+(1:obj.packetLen(ii)), :) = ...
                    obj.gain_per_packet(ii)*obj.rxWaveform(prev+(1:obj.packetLen(ii)), :);
                prev = prev+obj.packetLen(ii);
            end
        end
    end
end

author: analogdevicesinc


text:
function edgesEvalPlot( algs, nms, cols )
% Plot edge precision/recall results for directory of edge images.
%
% Enhanced replacement for plot_eval() from BSDS500 code:
%  http://www.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/
% Uses same format and is fully compatible with plot_eval. Use this
% function to plot the edge results created using edgesEvalDir.
%
% USAGE
%  edgesEvalPlot( algs, [nms], [cols] )
%
% INPUTS
%  algs       - {nx1} algorithm result directories
%  nms        - [{nx1}] algorithm names (for legend)
%  cols       - [{nx1}] algorithm colors
%
% OUTPUTS
%
% EXAMPLE
%
% See also edgesEvalDir
%
% Structured Edge Detection Toolbox      Version 3.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-microsoft.com]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the MSR-LA Full Rights License [see license.txt]

% parse inputs
if(nargin<2||isempty(nms)), nms={}; end; if(~iscell(nms)), nms={nms}; end
if(nargin<3||isempty(cols)), cols=repmat({'r','g','b','k','m'},1,100); end
if(~iscell(algs)), algs={algs}; end; if(~iscell(cols)), cols={cols}; end

% setup basic plot (isometric contour lines and human performance)
clf; box on; grid on; hold on;
line([0 1],[.5 .5],'LineWidth',2,'Color',.7*[1 1 1]);
for f=0.1:0.1:0.9, r=f:0.01:1; p=f.*r./(2.*r-f); %f=2./(1./p+1./r)
  plot(r,p,'Color',[0 1 0]); plot(p,r,'Color',[0 1 0]); end
if(1), h=plot(0.7235,0.9014,'o','MarkerSize',8,'Color',[0 .5 0],...
    'MarkerFaceColor',[0 .5 0],'MarkerEdgeColor',[0 .5 0]); end
set(gca,'XTick',0:0.1:1,'YTick',0:0.1:1);
grid on; xlabel('Recall'); ylabel('Precision');
axis equal; axis([0 1 0 1]);

% load results for every algorithm (pr=[T,R,P,F])
n=length(algs); hs=zeros(1,n); res=zeros(n,9); prs=cell(1,n);
for i=1:n, a=[algs{i} '-eval'];
  pr=dlmread(fullfile(a,'eval_bdry_thr.txt')); pr=pr(pr(:,2)>=1e-3,:);
  [~,o]=unique(pr(:,3)); R50=interp1(pr(o,3),pr(o,2),max(pr(o(1),3),.5));
  res(i,1:8)=dlmread(fullfile(a,'eval_bdry.txt')); res(i,9)=R50; prs{i}=pr;
end;

% sort algorithms by ODS score
[~,o]=sort(res(:,4),'descend'); res=res(o,:); prs=prs(o);
cols=cols(o); if(~isempty(nms)), nms=nms(o); end

% plot results for every algorithm (plot best last)
for i=n:-1:1
  hs(i)=plot(prs{i}(:,2),prs{i}(:,3),'-','LineWidth',3,'Color',cols{i});
  fprintf('ODS=%.3f OIS=%.3f AP=%.3f R50=%.3f',res(i,[4 7:9]));
  if(~isempty(nms)), fprintf(' - %s',nms{i}); end; fprintf('\n');
end

% show legend if nms provided (report best first)
hold off; if(isempty(nms)), return; end
for i=1:n, nms{i}=sprintf('[F=.%i] %s',round(res(i,4)*100),nms{i}); end
if(1), hs=[h hs]; nms=['[F=.80] Human'; nms(:)]; end
legend(hs,nms,'Location','sw');

end

author: Cloud-CV

repo: object-proposals


text:
function h=similaritygraph(coord,S,lim,lcolor)
%function h=similaritygraph(coord,[S],[lim],[lcolor])
%
%PURPOSE
%
%To visualize a weighted graph
%
%INPUT
%
%[An argument in brackets is optional. If it isn't  given or it's
% an empty matrix/string, the function will use a default value.]
%
% coord    (Mx2 matrix) coordinates
% [S]      (MxM matrix) similarity matrix: if this not given, no edges
%            are drawn, only vertices (black dots).
% [lim]     (vector) 1xK of monotonically increasing values between
%            0 and 1 sets limits for line colors: default [0.5 0.7 0.9]
% [lcolor] (Kx3 matrix) of RGB colors that define line colors
%            for each interval ]lim(1),lim(2)], ]lim(2),lim(3)],
%            ..., ]lim(end-1),lim(end)]  default: shades of red
%
%OUTPUT
%
% h         (struct) various graphic handles and other info
%             (see details)
%
%DETAILS
%
%Draws a weighted graph where the black points are vertices and
%the red shaded lines between them are edges. The level is
%determined by dividing values of S (the weights) into bins between the
% <= lim(1)                 not drawn
% ]lim(1),lim(2)]           colormap(1,:)    default: light red
%  ...                        ...               ...
% ]lim(end-1),lim(end)]     colormap(end,:)  default: bright red
% >lim(2)                   not drawn
%
%NOTE The function always adds to a plot (turns 'hold on' temporarily).
%
%Output is a structure of handles to graphic objects etc
%
% h.graph   (vector) graphic handles to all edges (lines)
% h.vertex  (vector) graphic handles to all vertices (markers)
% h.example (vector) contains one example (graphic handle) of each
%             different line type (gray level); needed for legend
% h.text    (cell array of strings) contains labels for lines in h.example
%             legend(h.example,h.text) sets a proper legend for the
%             graph line colors
%
%If there are no values in ith bin, the shade is left out, and
%h.graph(i)=[], h.example(i)=NaN, and h.text(i)='' will be set.
%
%USED IN
% icassoGraph

%COPYRIGHT NOTICE
%This function is a part of Icasso software library
%Copyright (C) 2003-2005 Johan Himberg
%
%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 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.

isHold=ishold;
hold on;

if nargin<3|isempty(lim),
  lim=[.5 .75 .9 1];
end

if nargin<2,
  S=zeros(size(coord,1));
end

lim=lim(:)';

if any(diff(lim)<=0),
  error('Limits must be monotonically increasing');
end

h.edge=[];
h.vertex=[];
h.example=[];
h.text=[];

% Change this to alter vertex color and line attrib. and line color
% scale
MARKERSIZE=3;
MARKERCOLOR=[0 0 0];
LINEWIDTH=1;

Nclass=length(lim)-1;
N=size(S,1);

if all(S(:)<lim(1)),
  isLines=0;
else
  isLines=1;
end

if isLines,
  if nargin<4|isempty(lcolor)
    lcolor=redscale(Nclass+1);
    lcolor=redscale(Nclass+1,0.9);
    lcolor(1:end-1,:)=lcolor(2:end,:);
  elseif size(lcolor,1)~=Nclass,
    error('Wrong number of colors.');
  end

  for i=1:Nclass,
    S_=S;
    S_(S_<=lim(i) | S_>lim(i+1))=0;
    links(:,:,i)=S_;
    hold on;
    [dummy,dummy,hd]=som_grid(S_,[N 1],'coord',coord,'marker','none', ...
			      'linecolor',lcolor(i,:),'linewidth',i.*LINEWIDTH);
    if ~isempty(hd),
      h.example(i)=hd(1); h.edge{i}=hd;
      h.text{i}=[sprintf('%0.2f',lim(i)) '<s_{ij}\leq ' sprintf('%0.2f',lim(i+1))];
    else
      % No lines in this slot
      h.example(i)=NaN; h.edge{i}=[];
      h.text{i}='';
    end
  end

  [dummy,h.vertex]=som_grid(S_,[N 1],'coord',coord,'marker','o','line','none',...
			    'markercolor',MARKERCOLOR,'markersize',MARKERSIZE);
else

  [dummy,h.vertex]=som_grid('rect',[N 1],'coord',coord,'marker','o','line','none',...
			    'markercolor',MARKERCOLOR,'markersize', ...
			    MARKERSIZE);
  h.edge=[];
  h.example=h.vertex(1);
  lim=[];
  h.text{1}='all s_{ij} below threshold: no lines shown';
end

if length(MARKERSIZE)==N,
  set(h.node,'markerfacecolor','none');
end

if ~isHold,
  hold off;
end

if nargout==0,
  clear h;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = sphereMesh(sphere, varargin)
%SPHEREMESH Create a 3D mesh representing a sphere.
%
%   [V, F] = sphereMesh(S)
%   Creates a 3D mesh representing the sphere S given by [xc yc zy r].
%
%   [V, F] = sphereMesh();
%   Assumes sphere is the unit sphere centered at the origin.
%
%   [V, F] = sphereMesh(S, 'nTheta', NT, 'nPhi', NP);
%   Specifies the number of discretisation steps for the meridians and the
%   parallels. Default values are nTheta = 16 and nPhi = 32.
%
%
%   Example
%     s = [10 20 30 40];
%     [v, f] = sphereMesh(s);
%     drawMesh(v, f);
%     view(3); axis equal; light; lighting gouraud;
%
%   See also
%     meshes3d, drawSphere, ellipsoidMesh, cylinderMesh, surfToMesh
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2012-10-25, using Matlab 7.9.0.529 (R2009b)
% Copyright 2012-2022 INRA - Cepia Software Platform

if nargin == 0
    sphere = [0 0 0 1];
end

% number of meridians
nPhi    = 32;

% number of parallels
nTheta  = 16;

% process input arguments
while length(varargin) > 1
    paramName = varargin{1};
    switch lower(paramName)
        case 'ntheta', nTheta = varargin{2};
        case 'nphi', nPhi = varargin{2};
        otherwise
            error(['Could not recognise parameter: ' paramName]);
    end
    varargin(1:2) = [];
end

% extract sphere data
xc = sphere(:,1);
yc = sphere(:,2);
zc = sphere(:,3);
r  = sphere(:,4);

% compute spherical coordinates
theta   = linspace(0, pi, nTheta+1);
phi     = linspace(0, 2*pi, nPhi+1);

% convert to cartesian coordinates
sintheta = sin(theta);
x = xc + cos(phi') * sintheta * r;
y = yc + sin(phi') * sintheta * r;
z = zc + ones(length(phi),1) * cos(theta) * r;

% convert to FV mesh
[vertices, faces] = surfToMesh(x, y, z, 'yperiodic', true);

% format output
varargout = formatMeshOutput(nargout, vertices, faces);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [epo_X, epo_Y] = feature_extraction_FFT(epo)

chan = {'PO3','POz','PO4','O1','Oz','O2'};

epo = proc_selectChannels(epo, chan);
epo = proc_selectIval(epo, [0 4000]);  % [0 4000]
dataset = permute(epo.x, [3,1,2]);

[tr, dp, ch] = size(dataset); % tr: trial, dp: time, ch: channel

nominal = [];
for i=1:size(epo.y,2)
    nominal(i) = find(epo.y(:,i),1)-1;
end

%% Fast Fourier Transform (FFT)
X_arr=[]; % make empty array
for k=1:tr % trials
    x=squeeze(dataset(k, :,:)); % data
    N=length(x);    % get the number of points
    kv=0:N-1;        % create a vector from 0 to N-1
    T=N/epo.fs;         % get the frequency interval
    freq=kv/T;       % create the frequency range
    X=fft(x)/N*2;   % normalize the data
    cutOff = ceil(N/2); % get the only positive frequency

    % take only the first half of the spectrum
    X=abs(X(1:cutOff,:)); % absolute values to cut off
    freq = freq(1:cutOff); % frequency to cut off
    XX = permute(X,[3 1 2]);
    X_arr=[X_arr; XX]; % save in array
end

%% frequency band
% f_gt = [11 7 5];  % 5.45, 8.75, 12
f_last = find( freq > 30, 1); % freq < 30Hz
X_arr = X_arr(:,1:f_last,:); %

%% get features
epo_X = permute(X_arr,[2,3,4,1]);
epo_Y = nominal;

author: PatternRecognition

repo: OpenBMI


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = inversekinematic_Viper_s1300(robot, T)
%   Solves the inverse kinematic problem for the ADEPT Viper_S1300 robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=inversekinematic_Viper_s1300 returns 8 possible solutions, thus,
%   Q is a 6x8 matrix where each column stores 6 feasible joint values.
%
%
%   Example code:
%
%   robot=load_robot('ADEPT', 'Viper_s1300');
%   q = [0 0 0 0 0 0];
%   T = directkinematic(robot, q);
%   %Call the inversekinematic for this robot
%   qinv = inversekinematic(robot, T);
%   check that all of them are feasible solutions!
%   and every Ti equals T
%   for i=1:8,
%        Ti = directkinematic(robot, qinv(:,i))
%   end
%	See also DIRECTKINEMATIC.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function [q] = inversekinematic_Viper_s1300(robot, T)

%initialize q,
%eight possible solutions are generally feasible
q=zeros(6,8);

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry at the reference for this robot, distance from the wrist to
%the end effector
L6=d(6);

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Compute the position of the wrist, being W the Z component of the end effector's system
W = T(1:3,3);

% Pm: wrist position
Pm = [Px Py Pz]' - L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution, obtain theta1
% by geometric methods
q1 = atan2(Pm(2), Pm(1));

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);

q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);

q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%Arrange solutions, there are 8 possible solutions so far.
% if q1 is a solution, q1* = q1 + pi is also a solution.
% For each (q1, q1*) there are two possible solutions
% for q2 and q3 (namely, elbow up and elbow up solutions)
% So far, we have 4 possible solutions. Howefer, for each triplet (theta1, theta2, theta3),
% there exist two more possible solutions for the last three joints, generally
% called wrist up and wrist down solutions. For this reason,
%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];

%leave only the real part of the solutions
q=real(q);

%Note that in this robot, the joint q3 has a non-simmetrical range. In this
%case, the joint ranges from 60 deg to -219 deg, thus, the typical normalizing
%step is avoided in this angle (the next line is commented). When solving
%for the orientation, the solutions are normalized to the [-pi, pi] range
%only for the theta4, theta5 and theta6 joints.

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));

% solve for the last three joints
% for any of the possible combinations (theta1, theta2, theta3)
for i=1:2:size(q,2),
    % use solve_spherical_wrist2 for the particular orientation
    % of the systems in this ABB robot
    % use either the geometric or algebraic method.
    % the function solve_spherical_wrist2 is used due to the relative
    % orientation of the last three DH reference systems.

    %use either one algebraic method or the geometric
    %qtemp = solve_spherical_wrist(robot, q(:,i), T, 1, 'geometric'); %wrist up
    qtemp = solve_spherical_wrist(robot, q(:,i), T, 1,'algebraic'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    %qtemp = solve_spherical_wrist(robot, q(:,i), T, -1, 'geometric'); %wrist down
    qtemp = solve_spherical_wrist(robot, q(:,i), T, -1, 'algebraic'); %wrist down
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%DH table parameters with which we calculate theta2, using
%geometric methods
L2=a(2);
L3=d(4);

%Offset distance between the centers of the reference systems of the links
%2 and 3
A1 = a(3);

%See geometry of the robot. Considering L4 calculate the offset between the
%centers of the reference systems of the links 2 and 3
L4 = sqrt(A1^2 + L3^2);

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

%Distance between the system 1 to the wrist
r = sqrt(p1(1)^2 + p1(2)^2);

beta = atan2(-p1(2), p1(1));
gamma = real(acos((L2^2+r^2-L4^2)/(2*r*L2)));   %Theorem of the cosine

if ~isreal(gamma)
    disp('WARNING:inversekinematic_Viper_s1300: the point is not reachable for this configuration, imaginary solutions');
    %gamma = real(gamma);
end

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta2. We add pi/2 to offset lags in our reference systems
q2(1) = pi/2 - beta - gamma; %elbow up
q2(2) = pi/2 - beta + gamma; %elbow down

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%DH table parameters with which we calculate theta3, using
%geometric methods
L2=a(2);
L3=d(4);

A1 = a(3);

%See geometry of the robot, like in the function q2
L4 = sqrt(A1^2 + L3^2);

%The delta angle is fixed because they are the lines that make up the gap
%between the links 2, 3 and 4
delta = real(acos((A1^2+L4^2-L3^2)/(2*A1*L4)));   %Theorem of the cosine

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

%Same as the function q2
r = sqrt(p1(1)^2 + p1(2)^2);

%Real angle between the links 2 and 3
ro = real(acos((L2^2 + L4^2 - r^2)/(2*L2*L4)));   %Theorem of the cosine

if ~isreal(ro)
   disp('WARNING:inversekinematic_Viper_s1300: the point is not reachable for this configuration, imaginary solutions');
   %ro = real(ro);
end

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta3. We add pi to offset lags in our reference systems
q3(1) = pi - ro - delta; %elbow up
q3(2) = pi + ro - delta; %elbow down

%remove complex solutions for q for the q1+pi solutions
function  qreal = arrange_solutions(q)
qreal=q(:,1:4);

%sum along rows if any angle is complex, for any possible solutions, then v(i) is complex
v = sum(q, 1);

for i=5:8,
    if isreal(v(i))
        qreal=[qreal q(:,i)]; %store the real solutions
    end
end

qreal = real(qreal);

repo: ARTE


text:
% [tInc,tIncCh] = hmrR_MotionArtifactByChannel(data, probe, mlActMan, mlActAuto, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh)
%
% UI NAME:
% Motion_Artifacts_By_Channel
%
%
% Identifies motion artifacts in an input data matrix d. If any active
% data channel exhibits a signal change greater than std_thresh or
% amp_thresh, then a segment of data around that time point is marked as a
% motion artifact. The channel-wise motion artifacts are recorded in the
% output matrix tIncCh. If any channel has a motion artifact, it is
% indicated by the vector tInc.
%
%
% INPUTS:
% data: SNIRF data structure data, containing time course data
% probe:   SNIRF data structure probe, containing probe source/detector geometry
% tIncMan: Cell array of vectors corresponding to the number of time bases in data.
%          tIncMan has been manually excluded. 0-excluded. 1-included. Vector same length as d.
% mlActMan: Cell array of vectors, one for each time base in data, specifying
%        active/inactive channels with 1 meaning active, 0 meaning inactive
% mlActAuto: Cell array of vectors, one for each time base in data, specifying
%        active/inactive channels with 1 meaning active, 0 meaning inactive
% tMotion: Check for signal change indicative of a motion artifact over
%     time range tMotion. Units of seconds.
% tMask: Mark data over +/- tMask seconds around the identified motion
%     artifact as a motion artifact. Units of seconds.
% STDEVthresh: If the signal d for any given active channel changes by more
%     that stdev_thresh * stdev(d) over the time interval tMotion, then
%     this time point is marked as a motion artifact. The standard deviation is
%     determined for each channel independently.
% AMPthresh: If the signal d for any given active channel changes by more
%     that amp_thresh over the time interval tMotion, then this time point
%     is marked as a motion artifact.
%
%
% OUTPUTS:
% tInc: a vector of length time points with 1's indicating data included
%       and 0's indicate motion artifact
% tIncCh: a matrix with #time points x #channels, with 1's indicating data
%       included and 0's indicating motion artifacts on a channel by
%       channel basis
%
% USAGE OPTIONS:
% Motion_Artifacts_By_Channel:  [tIncAuto, tIncAutoCh] = hmrR_MotionArtifactByChannel(dod, probe, mlActMan, mlActAuto, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh)
%
% PARAMETERS:
% tMotion: 0.5
% tMask: 1.0
% STDEVthresh: 50.0
% AMPthresh: 5.00
%
% PREREQUISITES:
% Intensity_to_Delta_OD: dod = hmrR_Intensity2OD( intensity )
%
% LOG:
% K. Perdue
% kperdue@nmr.mgh.harvard.edu
% Sept. 23, 2010
% Modified by DAB 3/17/2011 to not act on a single channel by default
% DAB 4/20/2011 added tMotion and tMask and hard coded stdev option.
% DAB 8/4/2011 added amp_thresh to work at same time as std_thresh
%
% TO DO:
% Consider tIncMan

function [tInc, tIncCh] = hmrR_MotionArtifactByChannel(data, probe, mlActMan, mlActAuto, tIncMan, tMotion, tMask, std_thresh, amp_thresh)

tInc   = cell(length(data), 1);
tIncCh = cell(length(data), 1);

% Input processing.  Check required inputs, throw errors if necessary.
if nargin<3
    error('First and second inputs must be data and probe SNIRF objects of type DataClass and ProbeClass.')
end
if ~isa(data, 'DataClass')
    error('First input must be data SNIRF objects of type DataClass.')
end
if ~isa(probe, 'ProbeClass')
    error('Second input must be probe SNIRF objects of type ProbeClass.')
end
if isempty(tIncMan)
    tIncMan = cell(length(data),1);
end
if isempty(mlActMan)
    mlActMan = cell(length(data),1);
end
if isempty(mlActAuto)
    mlActAuto = cell(length(data),1);
end

for iBlk = 1:length(data)

    d           = data(iBlk).GetDataTimeSeries();
    fs          = data(iBlk).GetTime();
    MeasList    = data(iBlk).GetMeasList();

    if length(fs)~=1
        fs = 1/(fs(2)-fs(1));
    end

    if isempty(tIncMan{iBlk})
        tIncMan{iBlk} = ones(size(d,1),1);
    end
    tInc{iBlk}   = ones(size(d,1),1);
    tIncCh{iBlk} = tIncCh_Initialize(tIncCh{iBlk}, d, MeasList);

    % get list of active channels
    mlActMan{iBlk} = mlAct_Initialize(mlActMan{iBlk}, MeasList);
    mlActAuto{iBlk} = mlAct_Initialize(mlActAuto{iBlk}, MeasList);
    lstAct1 = mlAct_Matrix2IndexList(mlActAuto{iBlk}, MeasList);
    lstAct2 = mlAct_Matrix2IndexList(mlActMan{iBlk}, MeasList);
    lstAct = unique([lstAct1(:)', lstAct2(:)']);

    % set artifact buffer for tMask seconds on each side of spike
    art_buffer = round(tMask*fs); % time in seconds times sample rate

    % LOOP OVER CHANNELS
    for iCh = 1:length(lstAct)

        lstActTmp = find(MeasList(:,1)==MeasList(lstAct(iCh),1) & ...
                         MeasList(:,2)==MeasList(lstAct(iCh),2) );

        % calculate std_diff for each channel
        std_diff = std( d(2:end, lstActTmp) - d(1:end-1, lstActTmp), 0, 1);

        % calculate max_diff across channels for different time delays
        max_diff = zeros(size(d,1)-1, length(lstActTmp));
        for ii = 1:round(tMotion*fs)
            max_diff=max([abs(d((ii+1):end,lstActTmp)-d(1:(end-ii),lstActTmp)); zeros(ii-1,length(lstActTmp))], max_diff);
        end

        % find indices with motion artifacts based on std_thresh or amp_thresh
        bad_inds = zeros(size(max_diff));
        mc_thresh = std_diff*std_thresh;
        for ii = 1:length(lstActTmp)
            bad_inds(:,ii) = max( [ max_diff(:,ii)>mc_thresh(ii), max_diff(:,ii)>amp_thresh ], [], 2);
        end
        bad_inds = find(max(bad_inds,[],2)==1);

        % Eliminate time points before or after motion artifacts
        if ~isempty(bad_inds)
            bad_inds = repmat(bad_inds, 1, 2*art_buffer+1) + repmat(-art_buffer:art_buffer, length(bad_inds), 1);
            bad_inds = bad_inds((bad_inds>0) & (bad_inds<=(size(d, 1)-1)));

            % exclude points that were manually excluded
            bad_inds(find(tIncMan{iBlk}(bad_inds)==0)) = [];

            % Set t and diff of data to 0 at the bad inds
            tInc{iBlk}(1+bad_inds) = 0; % bad inds calculated on diff so add 1
            tIncCh{iBlk}(1+bad_inds, lstAct(iCh)) = 0;
        end

    end % loop over channels

end

author: BUNPC


text:
function to_extract = find_node_indices(obj, varargin)
% Identify nodes to extract from brainpathway object given an integer
% vector of which nodes to select, cell array of strings for node names,
% or combination of the two. Returns a logical vector of nodes in set
%
% obj : a brainpathway object
%
% Optional inputs:
% a cell array of strings containing node names (matched to obj.node_labels)
% an integer vector of which nodes
%
% 'flatten' (not used here - would collapse atlas regions if used in select_atlas_regions)
%
% Outputs:
% to_extract: a logical

% -------------------------------------------------------------------------
% DEFAULTS AND INPUTS
% -------------------------------------------------------------------------

strings_to_find = [];
integers_to_find = [];
doflatten = false;

% optional inputs with default values
for i = 1:length(varargin)

    if iscell(varargin{i})
        strings_to_find = varargin{i};

    elseif isnumeric(varargin{i})
        integers_to_find = varargin{i};

    elseif ischar(varargin{i})
        switch varargin{i}

            case 'flatten', doflatten = true;

                %             case 'xxx', xxx = varargin{i+1}; varargin{i+1} = [];

            otherwise
                warning(['Unknown input string option:' varargin{i} '. Assuming it might be an atlas label. Place atlas labels in a cell array']);
                strings_to_find{1} = varargin{i};
        end
    end
end

% -------------------------------------------------------------------------
% INIT
% -------------------------------------------------------------------------

k = length(obj.region_indx_for_nodes); % size(obj.node_dat, 2);
to_extract = false(1, k);

% -------------------------------------------------------------------------
% FIND BY STRING
% -------------------------------------------------------------------------

for i = 1:length(strings_to_find)

    % Find which names match
    wh = ~cellfun(@isempty, strfind(obj.node_labels, strings_to_find{i}));

    to_extract = to_extract | wh;

end

% -------------------------------------------------------------------------
% FIND BY NUMBERS
% -------------------------------------------------------------------------

to_extract(integers_to_find) = true;

if ~any(to_extract)
    error('No nodes identified to extract.');
end

end % subfunction

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef math_model_cpf_accs < mp.math_model_cpf_acc & mp.mm_shared_pfcpf_accs
%MP.MATH_MODEL_CPF_ACCS  MATPOWER mathematical model for continuation power flow (CPF) problem.
%   ?
%
%   MP.MATH_MODEL_CPF_ACCS ... power flow ...
%
%   Properties
%       ? - ?
%
%   Methods
%       ?

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by 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.

%     properties
%     end

    methods
        function tag = form_tag(obj)
            tag = 'accs';
        end

        function name = form_name(obj)
            name = 'AC-cartesian-power';
        end

        function obj = add_node_balance_constraints(obj, nm, dm, mpopt)
            %% power balance constraints
            ad = obj.aux_data;
            fcn = @(x)node_balance_equations_cpf(obj, x, nm);
            obj.add_nln_constraint({'Pmis', 'Qmis', 'Vmis'}, [ad.npv+ad.npq;ad.npq;ad.npv], 1, fcn, []);
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function x = tsvread(f,v,hdr)
% Load text and numeric data from file
% FORMAT x = tsvread(f,v,hdr)
% f   - filename (can be gzipped) {txt,mat,csv,tsv,json}
% v   - name of field to return if data stored in a structure [default: '']
%       or index of column if data stored as an array
% hdr - detect the presence of a header row for csv/tsv [default: true]
%
% x   - corresponding data array or structure
%__________________________________________________________________________
%
% Based on spm_load.m from SPM12.
%__________________________________________________________________________

% Copyright (C) 2018, Guillaume Flandin, Wellcome Centre for Human Neuroimaging
% Copyright (C) 2018--, BIDS-MATLAB developers

%-Check input arguments
%--------------------------------------------------------------------------
if ~exist(f,'file')
    error('Unable to read file ''%s''',f);
end

if nargin < 2, v = ''; end
if nargin < 3, hdr = true; end % Detect

%-Load the data file
%--------------------------------------------------------------------------
[~,~,ext] = fileparts(f);
switch ext(2:end)
    case 'txt'
        x = load(f,'-ascii');
    case 'mat'
        x  = load(f,'-mat');
    case 'csv'
        % x = csvread(f); % numeric data only
        x = dsvread(f,',',hdr);
    case 'tsv'
        % x = dlmread(f,'\t'); % numeric data only
        x = dsvread(f,'\t',hdr);
    case 'json'
        x = bids.util.jsondecode(f);
    case 'gz'
        fz  = gunzip(f,tempname);
        sts = true;
        try
            x   = tsvread(fz{1});
        catch
            sts = false;
        end
        delete(fz{1});
        rmdir(fileparts(fz{1}));
        if ~sts, error('Cannot load ''%s''.',f); end
    otherwise
        try
            x = load(f);
        catch
            error('Unknown file format.');
        end
end

%-Return relevant subset of the data if required
%--------------------------------------------------------------------------
if isstruct(x)
    if isempty(v)
        fn = fieldnames(x);
        if numel(fn) == 1 && isnumeric(x.(fn{1}))
            x = x.(fn{1});
        end
    else
        if ischar(v)
            try
                x = x.(v);
            catch
                error('Data do not contain array ''%s''.',v);
            end
        else
            fn = fieldnames(x);
            try
                x = x.(fn{v});
            catch
                error('Invalid data index.');
            end
        end
    end
elseif isnumeric(x)
    if isnumeric(v)
        try
            x = x(:,v);
        catch
            error('Invalid data index.');
        end
    elseif ~isempty(v)
        error('Invalid data index.');
    end
end

%==========================================================================
% function x = dsvread(f,delim)
%==========================================================================
function x = dsvread(f,delim,header)
% Read delimiter-separated values file into a structure array
%  * header line of column names will be used if detected
%  * 'n/a' fields are replaced with NaN

%-Input arguments
%--------------------------------------------------------------------------
if nargin < 2, delim = '\t'; end
if nargin < 3, header = true; end % true: detect, false: no
delim = sprintf(delim);
eol   = sprintf('\n');

%-Read file
%--------------------------------------------------------------------------
S   = fileread(f);
if isempty(S), x = []; return; end
if S(end) ~= eol, S = [S eol]; end
S   = regexprep(S,{'\r\n','\r','(\n)\1+'},{'\n','\n','$1'});

%-Get column names from header line (non-numeric first line)
%--------------------------------------------------------------------------
h   = find(S == eol,1);
hdr = S(1:h-1);
var = regexp(hdr,delim,'split');
N   = numel(var);
n1  = isnan(cellfun(@str2double,var));
n2  = cellfun(@(x) strcmpi(x,'NaN'),var);
if header && any(n1 & ~n2)
    hdr     = true;
    try
        var = genvarname(var);
    catch
        var = matlab.lang.makeValidName(var,'ReplacementStyle','hex');
        var = matlab.lang.makeUniqueStrings(var);
    end
    S       = S(h+1:end);
else
    hdr     = false;
    fmt     = ['Var%0' num2str(floor(log10(N))+1) 'd'];
    var     = arrayfun(@(x) sprintf(fmt,x),(1:N)','UniformOutput',false);
end

% remove double delim
try
    while ~isempty(strfind(S,[delim delim]))
        S = strrep(S,[delim delim],delim);
    end
end

%-Parse file
%--------------------------------------------------------------------------
if exist('OCTAVE_VERSION','builtin') % bug #51093
    S = strrep(S,delim,'#');
    delim = '#';
end
if ~isempty(S)
    d = textscan(S,'%s','Delimiter',delim);
else
    d = {[]};
end
if rem(numel(d{1}),N), error('Varying number of delimiters per line.'); end
d = reshape(d{1},N,[])';
allnum = true;
for i=1:numel(var)
    sts = true;
    dd = zeros(size(d,1),1);
    for j=1:size(d,1)
        if strcmp(d{j,i},'n/a')
            dd(j) = NaN;
        else
            dd(j) = str2double(d{j,i}); % i,j considered as complex
            if isnan(dd(j)), sts = false; break; end
        end
    end
    if sts
        x.(var{i}) = dd;
    else
        x.(var{i}) = d(:,i);
        allnum     = false;
    end
end

if ~hdr && allnum
    x = struct2cell(x);
    x = [x{:}];
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
classdef magrcros < ZmapVGridFunction
    % MAGRCROS calcualtes Z-values in cross section
    %
    % TODO make this the first x-section function

    properties
        cutoff          datetime                % time of cut
        use_fixed_start logical     = false;
        periodA_start   datetime                % windows start time [where pre-window data starts]
        use_fixed_end   logical     = false;
        periodB_end     datetime                % windows end time [where post-window data ends]
        window_duration duration    = ZmapGlobal.Data.compare_window_dur;
        bin_dur         duration
    end

    properties(Constant)
        PlotTag         = 'zsection';
        ReturnDetails   = cell2table({ ... TODO update this. it hasn't been really done.
            ... VariableNames, VariableDescriptions, VariableUnits
            'AST',                  'Z-value comparing rate before to rate after cutoff','';...
            'LTA',                  'Z-value comparing rate outside window to inside window','';...
            'RUB',                  'Z-value comparing rate before cutoff [before window] to rate inside window','';...
            'PCT',                  'Compare ','';...
            'nBeforeCutoff',        '','';...
            'nAfterCutoff',         '','';...
            'nInWindow',            '','';...
            'nNotInWindow',         '','';...
            'DistAlongStrike',    'Distance along strike','km'...
            }, 'VariableNames', {'Names','Descriptions','Units'});

        CalcFields      = {...
            'AST','LTA','RUB','PCT',...
            'nBeforeCutoff','nAfterCutoff','nInWindow','nNotInWindow'}

        ParameterableProperties = ["cutoff", "use_fixed_start",...
            "periodA_start", "periodB_end", "use_fixed_end", "window_duration", "bin_dur"]

        References="";
        %Negative z-values indicate an increase in the seismicity rate, positive values a decrease.
            unit_options = {'seconds','hours','days','years'};
            unit_functions = {@seconds, @hours, @days, @years};
    end
    methods
        function obj=magrcros(zap,varargin)
            obj@ZmapVGridFunction(zap,'AST');
            report_this_filefun();

            % set the dependent variables here
            obj.periodA_start = min(obj.RawCatalog.Date);
            obj.periodB_end = max(obj.RawCatalog.Date);

            obj.parseParameters(varargin);
            warning('ZMAP:unimplemented','apparently still broken');

            obj.StartProcess();

            % consider this for future: uimenu(op1,'Label','Show Circles ','MenuSelectedFcn',@(~,~)plotcirc)
        end

        function InteractiveSetup(obj)
            % get the grid parameter

            % these are provided by the eventselection
            %ni = 100;
            %ra = ZG.ra;

            if isempty(obj.cutoff)
                % set the default cutoff to the time of the biggest event in catalog.
                % selects first event if multiple events are same size
                biggest=find(obj.RawCatalog.Magnitude==max(obj.RawCatalog.Magnitude) , 1);
                obj.cutoff = obj.RawCatalog.Date(biggest);
            end
            % make the interface
            %
            zdlg = ZmapDialog();
            zdlg.AddEventSelector('evsel', obj.EventSelector);
            if ~isempty(obj.Shape)
                zdlg.AddCheckbox('useGridFromShape', 'Limit grid to polygon', true,[],...
                    'Only evaluate for gridpoints within the polygon region. Does not restrict the catalog');
                zdlg.AddCheckbox('useCatFromShape', 'Limit catalog to polygon', false,[],...
                    'Only evaluate for events within the shape region. Does not restrict the grid');

            end

            % these provided by the grid
            %dd = 1.00;
            %dx = 1.00 ;
            %zdlg.AddEdit('dx_km','Spacing along strike [km]',dx,'spacing in horizontal plane');
            %zdlg.AddEdit('dd_km','Spacing in depth [km]',dd,'spacing in vertical plane');

            %zdlg.AddEdit('bin_dur','Time steps in days',ZG.bin_dur,'time steps in days');

            default_unit = find(obj.unit_options == "days");
            unitizer = obj.unit_functions{default_unit};

            zdlg.AddCheckbox('use_fixed_start', 'Fix StartTime', obj.use_fixed_end, 'fixed_start',...
                'Otherwise, the StartTime will depend on the catalog');

            zdlg.AddEdit('fixed_start',         'Start time',    obj.periodA_start,...
                'window size in specified units');

            zdlg.AddCheckbox('use_fixed_end',   'Fix EndTime',  obj.use_fixed_start, 'fixed_end',...
                'Otherwise, the StartTime will depend on the catalog');

            zdlg.AddEdit('fixed_end',           'End time',     obj.periodB_end,...
                'end time');

            zdlg.AddEdit('cutoff','Please enter date & time of cut:', obj.cutoff, 'Cutoff Date as yyyy-mm-dd hh:MM:ss');

            zdlg.AddEdit('win_dur',             'Window Size',  unitizer(obj.window_duration),...
                'window size in specified units');

            zdlg.AddPopup('win_dur_unit', 'Window Size Units:', obj.unit_options,   default_unit,...
                'Chooose units for window duration');
            zdlg.AddEdit('n_bins_in_window',    'Number of bins within window',...
                round(obj.window_duration/obj.ZG.bin_dur),...
                'Number of windows used to divide up the window [INTEGER]. this determines the bin size');

            [zparam,okPressed]=zdlg.Create('Name', 'Z-value xsection input parameters');
            if ~okPressed
                return
            end

            %dd=zparam.dd_km;
            %dx=zparam.dx_km;

            if isfield(zparam,'useCatFromShape') && zparam.useCatFromShape
                errordlg('not yet implemented: use catalog from shape. First limit catalog by shape, THEN call this function');
                return
            end

            SetValuesFromDialog(obj,zparam);
            obj.doIt();

        end  %  InteractiveSetup

        function SetValuesFromDialog(obj, res)
            if res.use_fixed_start
                obj.periodA_start = min(obj.RawCatalog.Date);
            else
                obj.periodA_start = res.fixed_start;
            end

            if res.use_fixed_end
                obj.periodB_end = max(obj.RawCatalog.Date);
            else
                obj.periodB_end = res.fixed_end;
            end
            unitizer = obj.unit_functions{res.win_dur_unit};
            obj.window_duration = unitizer(res.win_dur);
            obj.bin_dur = obj.window_duration/res.n_bins_in_window;
            obj.cutoff = res.cutoff;
            obj.EventSelector=res.evsel;

            if isfield(res,'useGridFromShape') && res.useGridFromShape
                obj.Grid = obj.Grid.MaskWithShape(obj.Shape);
            end
        end

        function Calculate(obj)

            % this is how to get a xsection grid from ZmapMainWindow (aa is handle to window)
            % xsz  = aa.xsec_zap('A - A''')
            % gr = xsz.Grid.MaskWithShape(aa.shape)
            % and the above should already be pulled into the object.

            assert(obj.cutoff > obj.periodA_start && obj.cutoff < obj.periodB_end,...
                'Cutoff date should be some time after %s and before %s', char(obj.periodA_start), char(obj.periodB_end));
            assert(obj.periodA_start < obj.periodB_end,'Invalid dates: Start date is after End date');
            assert (obj.cutoff + obj.window_duration > obj.periodB_end,...
                'window extends past end date. Manually set end date or change window');

            edges_for_cov = unique([obj.cutoff : - obj.bin_dur : obj.periodA_start, obj.cutoff : obj.bin_dur : obj.periodB_end]);

            obj.gridCalculations(@calculation_function);
            obj.Result.periodA_start = obj.periodA_start;
            obj.Result.periodB_end = obj.periodB_end;
            obj.Result.cutoff = obj.cutoff;
            obj.Result.bin_dur = obj.bin_dur;

            function out = calc_probability(old)
                %calculate probabliity, where old is one of the zmaps.
                % salvaged from vi_cucro
                valueMap = old;
                l = valueMap < 2.57;
                valueMap(l) = 2.65;
                pr = 0.0024 + 0.03*(valueMap - 2.57).^2;
                pr = (1-1./(exp(pr)));
                out = pr;
            end

            % catsave3('magrcros');

            % Plot the results
            %det = 'nop'
            %in2 = 'nocal'
            %menucros() -> which was at one point hooked up to incube, but not while I've been manipulating it. CGR

            function out=calculation_function(catalog)
                if catalog.Count <= 4
                    out = nan(size(obj.CalcFields));
                    return
                end

                %% do prelim calculations, so that they are only done ONCE per grid point

                % this had all been done with histogram and bins, but makes much more sense
                % to let the time periods dictate duration

                nPerBin = histcounts(catalog.Date,edges_for_cov);
                idxBeforeCutoff = edges_for_cov(2:end) <= obj.cutoff;
                idxAfterCutoff = ~idxBeforeCutoff;
                idxInWindow = idxAfterCutoff  & edges_for_cov(2:end) <= obj.cutoff + obj.window_duration;
                idxNotInWindow = ~idxInWindow;

                nBeforeCutoff = nPerBin(idxBeforeCutoff);
                nInWindow = nPerBin(idxInWindow); % NALL(2) : # IN window
                nAfterCutoff = nPerBin(idxAfterCutoff);
                nNotInWindow = nPerBin(idxNotInWindow);

                covBeforeCutoff = cov(nPerBin(idxBeforeCutoff));
                covAfterCutoff = cov(nPerBin(idxAfterCutoff));
                covInWindow = cov(nPerBin(idxInWindow));
                covNotInWindow = cov(nPerBin(idxNotInWindow));

                durBeforeCutoff = sum(idxBeforeCutoff) * obj.bin_dur;
                durAfterCutoff = sum(idxAfterCutoff) * obj.bin_dur;
                durInWindow = sum(idxInWindow) * obj.bin_dur;
                durNotInWindow = sum(idxNotInWindow) * obj.bin_dur;

                meanBeforeCutoff = mean(nBeforeCutoff);
                meanAfterCutoff  = mean(nAfterCutoff);
                meanInWindow = mean(nInWindow);
                meanNotInWindow = mean(nNotInWindow);

                out = [calc_ast() calc_lta() calc_rubberband() calc_percent(),...
                    sum(nBeforeCutoff), sum(nAfterCutoff), sum(nInWindow), sum(nNotInWindow)];

                % now comes the actual calculations for
                function out=calc_percent()
                    out =  -((meanBeforeCutoff-meanAfterCutoff)./meanBeforeCutoff)*100;
                    assert(numel(out)==1)
                end

                function out=calc_rubberband()
                    % loop over all point for rubber band : compare window to before cutoff
                    term1 = covBeforeCutoff / days(durBeforeCutoff);
                    term2 = covInWindow / days(durInWindow);
                    out =  (meanBeforeCutoff - meanInWindow) ./ (sqrt(term1 + term2));
                    assert(numel(out)==1)
                end

                function out=calc_ast()
                    % make the AST function map : compare before and after cutoff
                    term1 =covBeforeCutoff / days(durBeforeCutoff);
                    term2 = covAfterCutoff / days(durAfterCutoff);
                    out = (meanBeforeCutoff - meanAfterCutoff) ./ (sqrt(term1 + term2));
                    assert(numel(out)==1)
                end

                function out=calc_lta()
                    % Calculate LTA: compare window to everything else
                    term1 = covNotInWindow / days(durNotInWindow);
                    term2 = covInWindow  / days(durInWindow);
                    out =  (meanNotInWindow - meanInWindow)./(sqrt(term1+term2));
                    assert(numel(out)==1)
                end

                % removed something that calculated 'maz', since it was never referenced elsewhere, and
                % the original code was seriously convoluted. -CGR

            end %  calculation_function
        end % Calculate

    end % methods
    methods(Static)
        function h = AddMenuItem(parent, zapFcn, varargin) %xsec_zap
            % create a menu item
            label = 'Z-value section map';
            h = uimenu(parent, 'Label', label,...
                'MenuSelectedFcn', @(~,~)XZfun.magrcros(zapFcn()),...
                varargin{:});
        end

        %{
        %% UNMODIFIED/UNIMPLEMENTED BY CGR
        function obj=my_load()
            % Load exist z-grid
            [file1,path1] = uigetfile(['*.mat'],'z-value gridfile');
            if length(path1) > 1

                load([path1 file1])
                det = 'nop'
                in2 = 'nocal'
                menucros
            else
                return
            end
        end
        %}

    end % static methods

end

author: CelsoReyes


text:
%AOK TFR 2.0
%Deming Zhang
%deming.zhang@gmail.com
%

function varargout = tfagui(varargin)

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

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

% --- Executes just before TFA is made visible.
function TFA_OpeningFcn(hObject, eventdata, handles, varargin)
% Choose default command line output for TFA
vars = evalin('base','who');
set(handles.listbox_var,'String',vars);

handles.fs=1;
handles.lag=32;
handles.fftlen=512;
handles.tstep=1;
handles.vol=2;
handles.sig=[];
handles.tfr=[];
handles.t=[];
handles.f=[];
handles.begin=1;
handles.end=1;
handles.res=1;
handles.color=JET;
handles.times=6;
handles.vmanner=1;%linear
handles.realf=1;
%Ridges
handles.ci=[];
handles.ri=[];

handles.output = hObject;
% Update handles structure
guidata(hObject, handles);
set(gcf,'Color',[1 1 1]);

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

% --- Outputs from this function are returned to the command line.
function varargout = TFA_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;

% --- Executes on slider movement.
function slider_vol_Callback(hObject, eventdata, handles)
h=get(hObject);
handles.vol=h.Value;
handles.output = hObject;
guidata(hObject, handles);

% --- Executes on selection change in listbox_color.
function listbox_color_Callback(hObject, eventdata, handles)
h=get(hObject);
list_entries = get(hObject,'String');
index_selected = get(hObject,'Value');
str = list_entries{index_selected};
handles.color=str2num(str);
handles.output = hObject;
guidata(hObject, handles);
colormap(handles.color);

% --- Executes on button press in pushbutton_apply.
function pushbutton_apply_Callback(hObject, eventdata, handles)
s=handles.sig;
if imag(s)==0,      s=[s imag(hilbert(s))];
else                s=[real(s) imag(s)];end
s=s(handles.begin:handles.end,:);
res=handles.res;
if res>1,    s=resample(s,10,round(res*10));
elseif res>0 & res~=1,    s=resample(s,10,round(10*res));end
%slen=handles.slen;
slen=length(s);
fs=handles.fs;
lag=handles.lag;
fftlen=handles.fftlen;
tstep=handles.tstep;
vol=handles.vol;

fid = fopen('signal','w');
fprintf(fid,'%d %d %d %d %f\n',slen,lag,fftlen,tstep,vol);
fprintf(fid,'%f %f\n',s');
fclose(fid);
tic
!ATFR.exe signal >NUL
toc
Nt=floor((slen+lag-1)/tstep+1);
fid=fopen('TFR','rb');
TFR=fread(fid,[fftlen,Nt],'double');
fclose(fid);
!del TFR signal >NUL

t=linspace(0,(slen-1)/fs,slen);
f=linspace(-fs*(0.5-1/fftlen),fs*(0.5-1/fftlen),fftlen);

suffix=round(lag/(2*tstep));
TFR=TFR(:,suffix+1:Nt-suffix);
TFR=[TFR(1+fftlen/2:end,:);TFR(1:fftlen/2,:)];

handles.t=t;
handles.f=f;
%save 'd:\tfr.mat' TFR;
handles.tfr=TFR;

updateAOK(handles);

axes(handles.ax_pow);
p=( abs(fftshift(fft(s(:,1)+i*s(:,2)))) );

plot(p,linspace(f(1),f(end),length(p)));
if handles.realf==0,
axis([min(p)-1,max(p)+1,f(1),f(end)]);ylabel('Frequency(Hz)');title('Spectrum');
else axis([min(p)-1,max(p)+1,0,f(end)]);ylabel('Frequency(Hz)');title('Spectrum');
end
axes(handles.ax_wav);
ts=0:slen/(fs*(slen-1)):slen/fs;
if handles.realf==0,
plot(ts,s);axis tight;ylabel('WaveForm');
else plot(ts,s(:,1));axis tight;ylabel('WaveForm');
end
handles.output = hObject;
guidata(hObject, handles);

function edit_lag_Callback(hObject, eventdata, handles)
h=get(hObject);
handles.lag=str2num(h.String);
handles.output = hObject;
guidata(hObject, handles);

function edit_step_Callback(hObject, eventdata, handles)
h=get(hObject);
handles.tstep=str2num(h.String);
handles.output = hObject;
guidata(hObject, handles);

function edit_fft_Callback(hObject, eventdata, handles)
h=get(hObject);
handles.fftlen=str2num(h.String);
handles.output = hObject;
guidata(hObject, handles);

% --- Executes on selection change in listbox_manner.
function listbox_manner_Callback(hObject, eventdata, handles)
h=get(hObject);
handles.vmanner=h.Value;
handles.output = hObject;
guidata(hObject, handles);
updateAOK(handles);
% --------------------------------------------------------------------

% --- Executes on selection change in listbox_var.
function listbox_var_Callback(hObject, eventdata, handles)
sig=get_var_names(handles);
sig=sig(:);
handles.sig=sig;
handles.slen=length(sig);
handles.begin=1;
handles.end=handles.slen;
set(handles.edit_b,'String',num2str(handles.begin));
set(handles.edit_e,'String',num2str(handles.end));

if imag(sig)==0,
    axes(handles.ax_wav);
    plot([0:length(sig)-1]/handles.fs,sig);axis tight;ylabel('WaveForm');
    axes(handles.ax_pow);
    p=( abs(fftshift(fft(sig))) );
    plot(p(length(p)/2+1:end),linspace(0,0.5*handles.fs,length(p)/2));
else
    axes(handles.ax_wav);
    plot([0:length(sig)-1]/handles.fs,[real(sig) imag(sig)]);axis tight;ylabel('WaveForm');
    axes(handles.ax_pow);
    p=( abs(fftshift(fft(sig))) );
    plot(p,linspace(-0.5*handles.fs,0.5*handles.fs,length(p)));
end

handles.output = hObject;
guidata(hObject, handles);

function sig = get_var_names(handles)
list_entries = get(handles.listbox_var,'String');
index_selected = get(handles.listbox_var,'Value');
sig = list_entries{index_selected(1)};
sig=evalin('base',sig);

% --- Executes on button press in pushbutton_update.
function pushbutton_update_Callback(hObject, eventdata, handles)
vars = evalin('base','who');
set(handles.listbox_var,'String',vars);

% --- Executes during object creation, after setting all properties.

function edit_fs_Callback(hObject, eventdata, handles)
h=get(hObject);
handles.fs=str2num(h.String);
handles.output = hObject;
guidata(hObject, handles);

function updateAOK(handles)
t=handles.t;
f=handles.f;
if handles.realf==1,
    len=ceil(length(f)/2);
    tfr=handles.tfr(len:end,:);
    f=f(len:end);
else
    tfr=handles.tfr;
end
pt=0.005*2^handles.times;
if handles.vmanner==1,
    tfr=tfr.^pt;
else
    tfr=10*log10(tfr.^pt);
end
axes(handles.ax_aok);imagesc(t,f,tfr);xlabel('Time(s)');axis xy;ylabel('Frequency(Hz)');title('Time Frequency Representation');
%axes(handles.ax_aok);imagesc(t,f,tfr);xlabel('Time(s)');axis xy;ylabel('Frequency(Hz)');title('AOK Time-Frequency Representation');

% --- Executes on slider movement.
function slider_times_Callback(hObject, eventdata, handles)
h=get(hObject);
handles.times=h.Value;
updateAOK(handles);
handles.output = hObject;
guidata(hObject, handles);

% --- Executes on slider movement.
function edit_res_Callback(hObject, eventdata, handles)
h=get(hObject);
handles.res=str2num(h.String);
handles.output = hObject;
guidata(hObject, handles);

% --- Executes on slider movement.
function edit_b_Callback(hObject, eventdata, handles)
h=get(hObject);
handles.begin=str2num(h.String);
handles.output = hObject;
guidata(hObject, handles);
% --- Executes on slider movement.
function edit_e_Callback(hObject, eventdata, handles)
h=get(hObject);
handles.end=str2num(h.String);
handles.output = hObject;
guidata(hObject, handles);

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

% Hints: get(hObject,'String') returns contents of edit_cmd as text
%        str2double(get(hObject,'String')) returns contents of edit_cmd as a double
evalin('base',get(hObject,'String'));
set(hObject,'String','');
vars = evalin('base','who');
set(handles.listbox_var,'String',vars);

function axis_wav_Callback(hObject, eventdata, handles)
set(handles.edit_b,'String',num2str(handles.begin));
set(handles.edit_e,'String',num2str(handles.end));
handles.output = hObject;
guidata(hObject, handles);

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[ri,ci]=ridges(handles.tfr,12);
ri=handles.f(ri);
ci=handles.t(ci);
axes(handles.ax_aok);hold on;plot(ci,ri,'r.');hold off;
handles.ri=ri;handles.ci=ci;
handles.output = hObject;
guidata(hObject, handles);

function [ri,ci] = ridges(m,par)
 if nargin < 2,    par=3;  end;
 m = abs(m);
 [nrows,ncols] = size(m);
 ridges   = zeros(size(m));
 t      = 1:nrows;
 tplus  =[nrows 1:nrows-1 ];% rshift(t);
 tminus =[2:nrows 1];% lshift(t);
for i = 1:ncols,
 	x = m(:,i)';
    for j = 1:par,
		x=max([x(t); x(tplus); x(tminus)]);
    end;
	x = x(:);
	thresh = trimmean(x,5);
    ridges(:,i) = (~(m(:,i)<x)).*(m(:,i)>0.01*thresh);
 end;
 [ri,ci]=find(ridges>0);

% --- Executes on button press in pushbutton_load.
function pushbutton_load_Callback(hObject, eventdata, handles)
[fn, pn]=uigetfile('*.*', 'Mat or Ascii File to Load:');
if fn==0, return; end;
evalin('base', ['load(''' pn fn ''');']);
vars = evalin('base','who');
set(handles.listbox_var,'String',vars);

% --- Executes on button press in checkbox_real.
function checkbox_real_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox_real (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
h=get(hObject);
handles.realf=h.Value;
handles.output = hObject;
guidata(hObject, handles);

% Hint: get(hObject,'Value') returns toggle state of checkbox_real

author: Sable

repo: mcbench-benchmarks


text:
DO_TRAIN = 0;
DO_TEST = 1;

%% Training

if DO_TRAIN

datadir = '~/data/eccv08/msrc21/';

if 0
    load(fullfile(datadir, 'msSegmentFeatures_tu_tmp.mat'));
    load(fullfile(datadir, 'trainTestFn_tu.mat'));
    %train = train(randperm(numel(train)));
end

maxdata = 1000000;
maxtree = 1000;
bagsize = 25000;

% good vs bad segment classifier
if 1

[data, lab, w] = msFormatClassifierData(...
    segfeatures(train, :), seggood(train, :), trainw(train, :), maxdata);
nfs = [1 2 4 8 16 32 48 64 94];
valset = (ceil(0.1*(numel(lab))):ceil(0.4*numel(lab)));
ntrees = 100;
bagfract = bagsize / (numel(lab)-numel(valset));
[segnf, segerr] = rfSelectPoolSize(data, (lab+1)/2+1, bagfract, ntrees, ...
    maxtree, nfs, valset, w);
end
ntrees = 200;
[dt_seg, segoob] = rfTrainForest(data, (lab+1)/2+1, ...
    bagfract, ntrees, maxtree, segnf, [], w);

% label classifier
[data, lab, w] = msFormatClassifierData(...
    segfeatures(train, :), seglabel(train, :), trainw(train, :), maxdata);
%valset = (1:ceil(0.2*numel(lab)));
ntrees = 100;
bagfract = bagsize / (numel(lab)-numel(valset));
[labnf, laberr] = rfSelectPoolSize(data, lab, bagfract, ntrees, ...
    maxtree, nfs, valset, w);

ntrees = 500;
[dt_lab, laboob] = rfTrainForest(data, lab, ...
    bagfract, ntrees, maxtree, labnf, [], w);

save('~/data/rfexperiments/msrcRandForestJointClassifiers.mat', ...
    'maxdata', 'maxtree', 'nfs', 'segnf', 'segerr', 'labnf', 'laberr', ...
    'dt_seg', 'segoob', 'dt_lab', 'laboob');

end

%% Testing

if DO_TEST

if 0
    load(fullfile(datadir, 'msrc_imsegs.mat'));
    load(fullfile(datadir, 'testLabels.mat'));
    load(fullfile(datadir, 'msMultipleSegmentations_tu_tmp.mat'));
end

pg = msTestRf(imsegs(test), segfeatures(test, :), smaps(test), dt_lab, dt_seg, 0);
ignore = [5 8];
DO_SP = 1;
[acc, cm, classcount] = msAnalyzeResult(imsegs(test), labels(test), pg, DO_SP, ignore);
cm(ignore, :) = []; cm(:, ignore) = [];

save('~/data/rfexperiments/msrcResult_rf.mat', 'pg', 'acc', 'cm');

end

author: zouchuhang

repo: LayoutNet


text:
function gee_its_simple_test
%GEE_ITS_SIMPLE_TEST tests the "Gee! It's Simple!" package
% Exhaustive test of the "Gee! It's Simple!" package.  Returns the largest
% relative residual for any solution to A*x=b (using the inf norm).  This test
% exercises all statements in the package.  Note that the rand state is
% modified.
%
% Example:
%   gee_its_simple_test ;
%
% See also: gee_its_simple, gee_its_short, rand, mldivide, gee_its_simple_resid

% Copyright 2007, Timothy A. Davis.
% http://www.cise.ufl.edu/research/sparse

%-------------------------------------------------------------------------------
% error-handling tests
%-------------------------------------------------------------------------------

fprintf ('\nTesting error handling (expect error and warning messages):\n\n');

gunk = 0 ;
ok = 0 ;

lasterr ('') ;

try
    % too many inputs
    gee_its_simple_factorize (A,gunk) ;                                     %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % too many outputs
    [LU,p,rcnd,gunk] = gee_its_simple_factorize (A) ;                       %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % too few inputs
    gee_its_simple_factorize ;                                              %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % too many inputs
    x = gee_its_simple (A,b,gunk) ;                                         %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % too many outputs
    [x,rcnd,gunk] = gee_its_simple (A,b) ;                                  %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % too few inputs
    x = gee_its_simple ;                                                    %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % too many inputs
    x = gee_its_simple_forwardsolve (A,b,gunk) ;                            %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % too many outputs
    [x,gunk] = gee_its_simple_forwardsolve (A,b) ;                          %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % too few inputs
    x = gee_its_simple_forwardsolve ;                                       %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % too many inputs
    x = gee_its_simple_backsolve (A,b,gunk) ;                               %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % too many outputs
    [x,gunk] = gee_its_simple_backsolve (A,b) ;                             %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % too few inputs
    x = gee_its_simple_backsolve ;                                          %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % rectangular A
    x = gee_its_simple (eye (4,3), ones (4,1)) ;                            %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % A is 3D
    x = gee_its_simple (ones (9,3,3), ones (9,1)) ;                         %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % b is 3D
    x = gee_its_simple (eye (3,3), ones (3,3,3)) ;                          %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % dimensions of A and b do not matrix
    x = gee_its_simple (eye (3,3), ones (4,1)) ;                            %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % dimensions of L and b do not matrix
    x = gee_its_simple_forwardsolve (eye (3,3), ones (4,1)) ;               %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

try
    % dimensions of U and b do not matrix
    x = gee_its_simple_backsolve (eye (3,3), ones (4,1)) ;                  %#ok
catch
    ok = ok + 1 ;
    disp (lasterr) ;
end
fprintf ('\n') ;

% singular matrix
lastwarn ('') ;
x = gee_its_simple (0, 1) ;                                                 %#ok
[msg, id] = lastwarn ;
if (~isempty (msg) & ~isempty (id))                                         %#ok
    ok = ok + 1 ;
end
fprintf ('\n') ;

% ill-conditioned matrix
lastwarn ('') ;
x = gee_its_simple ([1e30 2e30 ; 1 1], [1 ; 1]) ;                           %#ok
[msg, id] = lastwarn ;
if (~isempty (msg) & ~isempty (id))                                         %#ok
    ok = ok + 1 ;
end

if (ok ~= 20)
    error ('test failed') ;
end

fprintf ('\n\nError-handing tests complete (all error messages and warnings\n');
fprintf ('shown above were expected).  Now testing for accuracy:\n\n') ;

%-------------------------------------------------------------------------------
% compare accuracy vs. backslash
%-------------------------------------------------------------------------------

maxerr1 = 0 ;   % largest residual for A\b (gee_its_sweet)
maxerr2 = 0 ;   % largest residual for gee_its_simple (A,b)
maxerr3 = 0 ;   % largest residual for gee_its_short (A,b)
maxerr4 = 0 ;   % largest residual for gee_its_too_short (A,b)
rmax = 0 ;      % largest relative difference in rcond

nmax = 50 ;     % largest dimension of A to test
cmax = 10 ;     % largest number of columns of b to test
ntrials = 2 ;   % number of trials for each x=A\b

rand ('state', 0) ;

for n = 0:nmax
    fprintf ('.') ;
    for c = 0:cmax
        for trial = 1:ntrials

            % set up the system
            A = rand (n) ;
            b = rand (n,c) ;

            % solve it four different ways
            x1 = gee_its_sweet (A,b) ;      % this is just a one-liner: x=A\b
            x2 = gee_its_simple (A,b) ;
            x3 = gee_its_short (A,b) ;
            x4 = gee_its_too_short (A,b) ;

            % get the relative residuals
            err1 = gee_its_simple_resid (A, x1, b) ;
            err2 = gee_its_simple_resid (A, x2, b) ;
            err3 = gee_its_simple_resid (A, x3, b) ;
            err4 = gee_its_simple_resid (A, x4, b) ;

            maxerr1 = max (maxerr1, err1) ;
            maxerr2 = max (maxerr2, err2) ;
            maxerr3 = max (maxerr3, err3) ;
            maxerr4 = max (maxerr4, err4) ;

            if (max ([err1 err2 err3]) > 1e-14)
                error ('test failed') ;
            end

            % test rcond
            if (n > 0)
                [L,U,p] = lu (A) ;                                          %#ok
                r1 = min (abs (diag (U))) / max (abs (diag (U))) ;
                [LU,p,r2] = gee_its_simple_factorize (A) ;
                if (r1 ~= 0)
                    r = abs (r1 - r2) / r1 ;
                    rmax = max (rmax, r) ;
                end
                if (r > 1e-10)
                    error ('test failed') ;
                end
            end
        end
    end
end
fprintf ('\n') ;

fprintf ('max residual for backslash:         %g\n', maxerr1) ;
fprintf ('max residual for gee_its_simple:    %g\n', maxerr2) ;
fprintf ('max residual for gee_its_short:     %g\n', maxerr3) ;
fprintf ('max residual for gee_its_too_short: %g (no pivoting!)\n', maxerr4) ;

fprintf ('\n\nAll tests passed OK\n') ;


text:
function planC = clearCache(planC, doseNum)
    indexS = planC{end};
    planC{indexS.dose}(doseNum).cachedColor = [];
    planC{indexS.dose}(doseNum).cachedTime = [];
    planC{indexS.dose}(doseNum).numCachedSlices = 0;
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR