text:
function vl_structuredNetwork_pairwiseModel_test_derivative
%vl_structuredNetwork_pairwiseModel_test_derivative tests vl_structuredNetwork_pairwiseModel

% create the following file from the variables of cnn_train_pairwiseModel.m right before the call of vl_structuredNetwork_pairwiseModel.m
load( 'vl_structuredNetwork_pairwiseModel_test_derivative.mat', 'net', 'im', 'gradients', 'labels', 'one' );
testEps = 1e-3;

fprintf('Computing the gradient ... ');
tStart = tic;
[lossValue, gradients, predictions] = vl_structuredNetwork_pairwiseModel(net, im, gradients, labels, one, ...
    'conserveMemory', true, ...
    'sync', true, ...
    'disableDropout', true ) ;
fprintf( '%f\n', toc(tStart) );

maxGroupTests = 100;
rng(1);

% test derivatives
for iLayer = length( net.layers ) : -1 : 1
    if ~isequal( net.layers{iLayer}.type, 'conv' )
        continue;
    end
    fprintf('Layer %d: %s\n', iLayer, net.layers{iLayer}.name);

    empiricalDerivative = zeros( numel(gradients{ iLayer }.dzdw{2}), 1, 'like', gradients{ iLayer }.dzdw{2});
    fprintf('Number of bias derivatives: %d\n', numel(empiricalDerivative));

    randOrder = randperm(numel(empiricalDerivative));
    numTests = min( numel(empiricalDerivative), maxGroupTests);
    for iValueIndex = 1 : numTests
        if mod(iValueIndex, 1000) == 0
            fprintf('Derivative #%d\n', iValueIndex);
        end
        iValue = randOrder(iValueIndex);

        initValue = net.layers{iLayer}.weights{2}(iValue);
        net.layers{iLayer}.weights{2}(iValue) = initValue + testEps;

        [lossValue_test, ~, predictions_test] = vl_structuredNetwork_pairwiseModel(net, im, gradients, labels, [], ...
            'conserveMemory', true, ...
            'sync', true, ...
            'disableDropout', true) ;

        empiricalDerivative(iValue) = sum(lossValue_test - lossValue) / testEps;

        net.layers{iLayer}.weights{2}(iValue) = initValue;
    end

    testIndices = randOrder(1 : numTests);
    emphiricalGradient = empiricalDerivative(testIndices);
    computedGradient = gradients{ iLayer }.dzdw{2}(testIndices);
    derivativeError = gather( norm(emphiricalGradient(:) - computedGradient(:)) / norm(computedGradient(:)) );

    fprintf('Relative error of bias derivatives: %f\n', derivativeError );
    fprintf('Norm of tested derivatives: %f\n', norm(computedGradient(:)) );

    empiricalDerivative = zeros( numel(gradients{ iLayer }.dzdw{1}), 1, 'like', gradients{ iLayer }.dzdw{1});
    fprintf('Number of filter derivatives: %d\n', numel(empiricalDerivative));

    randOrder = randperm(numel(empiricalDerivative));
    numTests = min( numel(empiricalDerivative), maxGroupTests);
    for iValueIndex = 1 : numTests
        if mod(iValueIndex, 1000) == 0
            fprintf('Derivative #%d\n', iValueIndex);
        end

        iValue = randOrder(iValueIndex);

        initValue = net.layers{iLayer}.weights{1}(iValue);
        net.layers{iLayer}.weights{1}(iValue) = initValue + testEps;

        [lossValue_test, ~, predictions_test] = vl_structuredNetwork_pairwiseModel(net, im, gradients, labels, [], ...
            'conserveMemory', true, ...
            'sync', true, ...
            'disableDropout', true) ;

        empiricalDerivative(iValue) = sum(lossValue_test - lossValue) / testEps;

        net.layers{iLayer}.weights{1}(iValue) = initValue;
    end

    testIndices = randOrder(1 : numTests);
    emphiricalGradient = empiricalDerivative(testIndices);
    computedGradient = gradients{ iLayer }.dzdw{1}(testIndices);
    derivativeError = gather( norm(emphiricalGradient(:) - computedGradient(:)) / norm(computedGradient(:)) );

    fprintf('Relative error of filter derivatives: %f\n', derivativeError );
    fprintf('Norm of tested derivatives: %f\n', norm(computedGradient(:)) );

end

end

author: aosokin


text:
function check = chi_square_noncentral_check ( a, b )

%*****************************************************************************80
%
%% CHI_SQUARE_NONCENTRAL_CHECK checks the parameters of the noncentral Chi Squared PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer A, the parameter of the PDF.
%    1.0 <= A.
%
%    Input, real B, the noncentrality parameter of the PDF.
%    0.0 <= B.
%
  if ( a < 1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CHI_SQUARE_NONCENTRAL_CHECK - Fatal error!\n' );
    fprintf ( 1, '  A < 1.\n' );
    check = 0;
    error ( 'CHI_SQUARE_NONCENTRAL_CHECK - Fatal error!' );
  end

  if ( b < 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CHI_SQUARE_NONCENTRAL_CHECK - Fatal error!\n' );
    fprintf ( 1, '  B < 0.\n' );
    check = 0;
    error ( 'CHI_SQUARE_NONCENTRAL_CHECK - Fatal error!' );
  end

  check = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef VehicleArticulatedNonlinear < VehicleDynamicsLateral.VehicleArticulated
    % VehicleArticulatedNonlinear Nonlinear articulated vehicle model.
    %
    % It inherits properties from VehicleArticulated.

    methods
        % Constructor
        function self = VehicleArticulatedNonlinear()
            self.mF0    = 5200;    % Mass at front axle of the tractor without semi-trailer [kg]
            self.mR0    = 2400;    % Mass at rear axle of the tractor without semi-trailer [kg]
            self.mF     = 6000;    % Mass at front axle of the tractor with semi-trailer [kg]
            self.mR     = 10000;   % Mass at rear axle of the tractor with semi-trailer [kg]
            self.mM     = 17000;   % Mass at the axle of the semi-trailer [kg]
            self.IT     = 46000;   % Moment of inertia of the tractor [kg.m2]
            self.IS     = 450000;  % Moment of inertia of the semi-trailer [kg.m2]
            self.lT     = 3.5;     % Wheelbase of the tractor [m]
            self.lS     = 7.7;     % Distance from articulation (fifth wheel) to semi-trailer axle [m]
            self.c      = -0.3;    % Position of articulation (fifth wheel) to semi-trailer axle [m]
            self.nF     = 2;       % Number of wheels at front axle
            self.nR     = 4;       % Number of wheels at rear axle
            self.nM     = 8;       % Number of wheels at semi-trailer axle
            self.wT     = 2.6;     % Tractor width [m]
            self.wS     = 2.4;     % Semi-trailer width [m]
            self.muy    = 0.3;     % Friction
            self.deltaf = 0;       % Steering angle front [rad]
            self.deltar = 0;       % Steering angle rear [rad]
            self.deltam = 0;       % Steering angle semi-trailer [rad]
            self.Fxf    = 0;       % Longitudinal force at front axle [N]
            self.Fxr    = 0;       % Longitudinal force at rear axle [N]
            self.Fxm    = 0;       % Longitudinal force at semi-trailer axle [N]
        end

        function dx = Model(self,t, states,tspan)
            % dx = VehicleModel.Model(t,states,tspan)
            % dx = VehicleModel.MassMatrix(t,states,tspan)
            %
            % Sintax
            % |dx = _VehicleModel_.Model(t,states,tspan)|
            %
            % |dx = _VehicleModel_.MassMatrix(t,states,tspan)|
            %
            % Arguments
            % The following table describes the input arguments:
            %
            % t - Time
            % states - Model state variables: [XT YT PSI PHI VT ALPHAT dPSI dPHI]
            % tspan - Time span

            % Vehicle parameters
            mT  = self.mT;
            mS  = self.mS;
            % IT = self.IT;
            % IS = self.IS;
            a   = self.a;
            b   = self.b;
            c   = self.c;
            d   = self.d;
            e   = self.e;
            nF  = self.nF;
            nR  = self.nR;
            nM  = self.nM;

            g = 9.81;

            FzF = self.mF * g;
            FzR = self.mR * g;
            FzM = self.mM * g;
            muy = self.muy;

            % States
            X       = states(1,1);
            Y       = states(2,1);
            PSI     = states(3,1);
            PHI     = states(4,1);
            VT      = states(5,1);
            ALPHAT  = states(6,1);
            dPSI    = states(7,1);
            dPHI    = states(8,1);

            if isa(self.deltaf,'function_handle')
                % Closed loop steering
                deltaf = self.deltaf([X;Y;PSI;PHI;VT;ALPHAT;dPSI;dPHI],t);
            elseif length(self.deltaf)>1
                % Open loop variable steering
                deltaf = interp1(tspan,self.deltaf,t);
            else
                % Open loop single value steering
                deltaf = self.deltaf;
            end

            if isa(self.deltar,'function_handle')
                % Closed loop steering
                deltar = self.deltar([X;Y;PSI;PHI;VT;ALPHAT;dPSI;dPHI],t);
            elseif length(self.deltar)>1
                % Open loop variable steering
                deltar = interp1(tspan,self.deltar,t);
            else
                % Open loop single value steering
                deltar = self.deltar;
            end

            if isa(self.deltam,'function_handle')
                % Closed loop steering
                deltam = self.deltam([X;Y;PSI;PHI;VT;ALPHAT;dPSI;dPHI],t);
            elseif length(self.deltam)>1
                % Open loop variable steering
                deltam = interp1(tspan,self.deltam,t);
            else
                % Open loop single value steering
                deltam = self.deltam;
            end

            % Slip angles
            ALPHAF = - deltaf + atan2(a*dPSI + VT*sin(ALPHAT), VT*cos(ALPHAT));
            ALPHAR = - deltar + atan2(VT*sin(ALPHAT) - b*dPSI, VT*cos(ALPHAT));
%           ALPHAM = atan2(((d + e)*(dPHI - dPSI) + VT * sin(ALPHAT + PHI) - b * dPSI * cos(PHI) - c * dPSI * cos(PHI)),(VT * cos(ALPHAT + PHI) + b * dPSI * sin(PHI) + c * dPSI * sin(PHI)));
            ALPHAM = - deltam + atan2((d + e)*(dPHI - dPSI) + VT*sin(ALPHAT + PHI) - dPSI*cos(PHI)*(b + c), VT*cos(ALPHAT + PHI) + dPSI*sin(PHI)*(b + c));
            % Longitudinal forces
            % Front
            if isa(self.Fxf,'function_handle')
                FxF = self.Fxf([X;Y;PSI;PHI;VT;ALPHAT;dPSI;dPHI],t);
            elseif length(self.Fxf)>1
                FxF = interp1(tspan,self.Fxf,t);
            else
                FxF = self.Fxf;
            end
            % Rear
            if isa(self.Fxr,'function_handle')
                FxR = self.Fxr([X;Y;PSI;PHI;VT;ALPHAT;dPSI;dPHI],t);
            elseif length(self.Fxr)>1
                FxR = interp1(tspan,self.Fxr,t);
            else
                FxR = self.Fxr;
            end
            % Semitrailer
            if isa(self.Fxm,'function_handle')
                FxM = self.Fxm([X;Y;PSI;PHI;VT;ALPHAT;dPSI;dPHI],t);
            elseif length(self.Fxm)>1
                FxM = interp1(tspan,self.Fxm,t);
            else
                FxM = self.Fxm;
            end

            % Lateral forces
            FyF = nF * self.tire.Characteristic(ALPHAF, FzF/nF, muy);
            FyR = nR * self.tire.Characteristic(ALPHAR, FzR/nR, muy);
            FyM = nM * self.tire.Characteristic(ALPHAM, FzM/nM, muy);

            t2 = cos(ALPHAT);
            t3 = cos(PSI);
            t4 = sin(PHI);
            t5 = sin(PSI);
            t6 = cos(deltam);
            t7 = sin(deltam);
            t8 = PSI+deltaf;
            t9 = PSI+deltar;
            t10 = dPHI.^2;
            t11 = dPSI.^2;
            t12 = ALPHAT+PHI;
            t13 = ALPHAT+PSI;
            t21 = -PHI;
            t22 = -PSI;
            t23 = -deltam;
            t14 = cos(t12);
            t15 = cos(t13);
            t16 = sin(t13);
            t17 = cos(t8);
            t18 = cos(t9);
            t19 = sin(t8);
            t20 = sin(t9);
            t24 = FyM.*d.*t6;
            t25 = FyM.*e.*t6;
            t26 = FxM.*d.*t7;
            t27 = FxM.*e.*t7;
            t28 = PHI+t22;
            t29 = PHI+t23;
            t31 = PSI+deltam+t21;
            t30 = cos(t28);
            t32 = sin(t28);
            t33 = cos(t29);
            t34 = sin(t29);
            t35 = cos(t31);
            t36 = sin(t31);
            t37 = VT.*d.*dPSI.*mS.*t14;
            f = [VT.*t15;VT.*t16;dPSI;dPHI;FxF.*t17+FxM.*t35+FxR.*t18-FyF.*t19-FyM.*t36-FyR.*t20-b.*mS.*t3.*t11-c.*mS.*t3.*t11-d.*mS.*t10.*t30-d.*mS.*t11.*t30+VT.*dPSI.*mS.*t16+VT.*dPSI.*mT.*t16+d.*dPHI.*dPSI.*mS.*t30.*2.0;FxF.*t19+FxM.*t36+FxR.*t20+FyF.*t17+FyM.*t35+FyR.*t18-b.*mS.*t5.*t11-c.*mS.*t5.*t11+d.*mS.*t10.*t32+d.*mS.*t11.*t32-VT.*dPSI.*mS.*t15-VT.*dPSI.*mT.*t15-d.*dPHI.*dPSI.*mS.*t32.*2.0;-t24-t25-t26-t27+t37+FxM.*b.*t34-FyM.*b.*t33+FxM.*c.*t34-FyM.*c.*t33+FyF.*a.*cos(deltaf)-FyR.*b.*cos(deltar)+FxF.*a.*sin(deltaf)-FxR.*b.*sin(deltar)+VT.*b.*dPSI.*mS.*t2+VT.*c.*dPSI.*mS.*t2-b.*d.*mS.*t4.*t10-c.*d.*mS.*t4.*t10+b.*d.*dPHI.*dPSI.*mS.*t4.*2.0+c.*d.*dPHI.*dPSI.*mS.*t4.*2.0;t24+t25+t26+t27-t37-b.*d.*mS.*t4.*t11-c.*d.*mS.*t4.*t11];

            dx = f;
        end

        function [value,isterminal,direction] = velocity(~,~,states)
            % If the velocity is less than 0.1m/s the integrator stops.
            % The MassMatrix is singular when the velocity is 0 m/s.
            value = states(5,1) - 0.1;
            isterminal = 1;
            direction = -1;
        end

        function M = MassMatrix(self,~,states)
            % Vehicle Parameters
            mT = self.mT;
            mS = self.mS;
            IT = self.IT;
            IS = self.IS;
            % a = self.a;
            b = self.b;
            c = self.c;
            d = self.d;
            % e = self.e;
            % deltaf = self.deltaf;
            % nF = self.nF;
            % nR = self.nR;
            % nM = self.nM;

            % g = 9.81;

            % FzF = self.mF * g;
            % FzR = self.mR * g;
            % FzM = self.mM * g;
            % muy = self.muy;

            % States
            PSI = states(3,1);
            PHI = states(4,1);
            VT = states(5,1);
            ALPHAT = states(6,1);
            % dPSI = states(7,1);
            % dPHI = states(8,1);

            t2 = cos(ALPHAT);
            t3 = cos(PHI);
            t4 = sin(ALPHAT);
            t5 = b+c;
            t6 = mS+mT;
            t7 = d.^2;
            t8 = ALPHAT+PHI;
            t9 = ALPHAT+PSI;
            t14 = -IS;
            t15 = -PSI;
            t10 = cos(t8);
            t11 = cos(t9);
            t12 = sin(t8);
            t13 = sin(t9);
            t16 = mS.*t7;
            t17 = b.*d.*mS.*t3;
            t18 = c.*d.*mS.*t3;
            t19 = PHI+t15;
            t20 = -t16;
            t21 = cos(t19);
            t22 = sin(t19);
            t23 = -t17;
            t24 = -t18;
            t25 = t14+t20+t23+t24;
            M = reshape([1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t6.*t11,t6.*t13,-mS.*(b.*t4+c.*t4+d.*t12),d.*mS.*t12,0.0,0.0,0.0,0.0,-VT.*t6.*t13,VT.*t6.*t11,-VT.*mS.*(b.*t2+c.*t2+d.*t10),VT.*d.*mS.*t10,0.0,0.0,0.0,0.0,mS.*(d.*t22.*2.0-t5.*sin(PSI).*2.0).*(-1.0./2.0),mS.*(d.*t21.*2.0+t5.*cos(PSI).*2.0).*(-1.0./2.0),IS+IT+t16+t17.*2.0+t18.*2.0+b.^2.*mS+c.^2.*mS+b.*c.*mS.*2.0,t25,0.0,0.0,0.0,0.0,d.*mS.*t22,d.*mS.*t21,t25,IS+t16],[8,8]);
        end
    end

end

author: andresmendes

repo: Vehicle-Dynamics-Lateral


text:
function p3 = rdivide(p1,p2)
%PREAL/RDIVIDE Overloaded rdivide (./) operator for class preal.

global useUnitsFlag

if ~(useUnitsFlag) % If physunits is disabled...
    p3=rdivide(double(p1),double(p2)); % ... treat as double.
    return
end

p1=preal(p1);
p2=preal(p2);
if isscalar(p1)
    p1=preal(p1.value*ones(size(p2)),p1.units);
elseif isscalar(p2)
    p2=preal(p2.value*ones(size(p1)),p2.units);
elseif ndims(p1)==ndims(p2)&size(p1)==size(p2)
    % do nothing
else
    error('Array dimensions must agree')
end
p3=preal(ones(size(p1)));
for k=1:numel(p3)
    p3(k).value=p1(k).value/p2(k).value;
    p3(k).units=p1(k).units-p2(k).units;
end

author: Sable

repo: mcbench-benchmarks


text:
function binname=mcpath(fname, ext)
%
% binname=mcpath(fname)
%
% get full executable path by prepending a command directory path
% parameters:
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%    fname: input, a file name string
%
% output:
%    binname: output, full file name located in the bin directory
%
%    if global variable ISO2MESH_BIN is set in 'base', it will
%    use [ISO2MESH_BIN filesep cmdname] as the command full path,
%    otherwise, let matlab pass the cmdname to the shell, which
%    will search command in the directories listed in system
%    $PATH variable.
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

p=getvarfrom({'caller','base'},'ISO2MESH_BIN');
binname=[];
if(isempty(p))
	% the bin folder under iso2mesh is searched first
	tempname=[fileparts(which(mfilename)) filesep 'bin' filesep fname];
	if(exist([fileparts(which(mfilename)) filesep 'bin'])==7)
        if(nargin>=2 && ~isempty(ext))
            if(exist([tempname ext],'file'))
                binname=[tempname ext];
            else
                binname=fname; % use binary without suffix on system PATH
            end
        else
		    binname=tempname;
        end
	else
		binname=fname;
	end
else
	binname=[p filesep fname];
end
% on 64bit windows machine, try 'exename_x86-64.exe' first
if(ispc && ~isempty(regexp(computer,'64', 'once')) && isempty(regexp(fname,'_x86-64$', 'once')))
    w64bin=regexprep(binname,'(\.[eE][xX][eE])*$','_x86-64.exe','emptymatch');
    if(exist(w64bin,'file'))
        binname=w64bin;
    end
end
% if no such executable exist in iso2mesh/bin, find it in PATH env variable
if(nargin>=2 && ~exist(binname,'file'))
    binname=fname;
end

author: fangq


% icp030001_new = poseSparsification( icp030001, '/home/yh/mapModel/2018/08.03/keep_0301.txt', 0.20);
% getDistribution(icp030001);
% hold on;
% getDistribution(icp030001_new);

% [pose10_new, pose10_keep] = ForKITTIPoses( '/home/yh/mapModel/2018/08.02/kitti10_pose.txt', '/home/yh/mapModel/2018/08.02/kitti10_keep.txt', pose10, 0.70);
%
% getDistribution(pose10_new);
% hold on;
% getDistribution(pose10_keep);

%% test

% icp0901_new = poseSparsification( icp0901, '/home/yh/mapModel/2018/08.03/test/keep_0901.txt', 0.20);
% getDistribution(icp0901);
% hold on;
% getDistribution(icp0901_new);

%% Chery test'

% chery_pose = poseSparsification( globalpose, '/home/yh/mapModel/2018/09.22/keep_chery.txt', 0.20);
% getDistribution(globalpose);
% hold on;
% getDistribution(chery_pose);

%% New YQ data pepare, the former are too dense\
%
% icp0301_new = poseSparsification( icp0301, '/home/yh/mapModel/2018/10.18/keep_0301.txt', 1.0);
% getDistribution(icp0301);
% hold on;
% getDistribution(icp0301_new);

%% Re-index for exp , 2019.07

% icp0301_new = poseSparsification_yq( '/home/yh/mapModel/2019/07.11/keep_0301.txt', '/home/yh/mapModel/2019/07.11/index_0301.txt', icp0301, 0.2);

% icp0301_new = poseSparsification_yq( '/home/yh/mapModel/2019/07.11/keep_0901.txt', '/home/yh/mapModel/2019/07.11/index_0901.txt', icp0901, 0.2);

% park_new = poseSparsification_yq( '/home/yh/mapModel/2019/07.23/keep_park.txt', '/home/yh/mapModel/2019/07.23/index_park.txt', parkpose, 0.2);

% kitti_new = poseSparsification_kitti( '/home/yh/mapModel/2019/07.21/keep_kitti.txt', '/home/yh/mapModel/2019/07.21/index_kitti.txt', Untitled, 0.5);

% kitti_new = poseSparsification_kitti( '/home/yh/mapModel/2019/07.27/keep_kitti.txt', '/home/yh/mapModel/2019/07.27/index_kitti.txt', X07, 0.5);


text:
function [poly, keepInds] = simplifyPolygon(poly, varargin)
%SIMPLIFYPOLYGON  Douglas-Peucker simplification of a polygon
%
%   POLY2 = simplifyPolygon(POLY, TOL)
%   Simplifies the input polygon using the Douglas-Peucker algorithm.
%
%   Example
%     elli = [20 30 40 20 30];
%     poly = ellipseToPolygon(elli, 500);
%     poly2 = simplifyPolygon(poly, 1); % use a tolerance equal to 1.
%     figure; hold on;
%     drawEllipse(elli);
%     drawPoint(poly2, 'mo');
%
%   See also
%   polygons2d, smoothPolygon, simplifyPolyline, resamplePolygon
%
%   References
%   http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
%

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

% call the simplifyPolyline function by ensuring the last vertex is present
poly = poly([1:end 1], :);
[poly, keepInds] = simplifyPolyline(poly, varargin{:});

% remove last vertex
poly(end, :) = [];
keepInds(end) = [];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
%compare max, min a mean manipulability for vanilla
%load('/Users/arturogilaparicio/Desktop/arte/robots/RETHINK/sgo_v0.5/experiment1/mat/experiment1A.mat')

[correct_manips]=eval_experiment(Gout, random_manips);

[y,i]=max(sum(correct_manips'));
max_manip = correct_manips(i,:);

[y,i]=min(sum(correct_manips'));
min_manip = correct_manips(i,:);

mean_manip = mean(correct_manips);

figure, plot(max_manip), hold
plot(min_manip)
plot(mean_manip)

repo: ARTE


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = svol_2(a, b, r, n, f, k, t)
% computes the SABR implied BS volatility due to Hagan's formula

	z = n./a.*(f*k).^((1-b)/2).*log(f./k);
	x = log((sqrt(1 - 2*r*z + z.^2) + z - r)/(1-r));
	Term1 = a ./ (f*k).^((1-b)/2) ./ (1 + (1-b)^2/24*log(f./k).^2 ...
        + (1-b)^4/1920*log(f./k).^4);
	Term2 = z ./ x;
	Term2(abs(x-z) < 1e-008) = 1;           % account for ATM
	Term3 = 1 + ((1-b)^2/24*a^2./(f*k).^(1-b) ...
        + r*b*n*a/4./(f*k).^((1-b)/2) + (2-3*r^2)/24*n.^2)*t;
	y = Term1.*Term2.*Term3;
end


text:
function [boolval, faxis] = checkfreq(varargin)

% last input is always the required string
tol      = varargin{end};
required = varargin{end-1};
varargin = varargin(1:end-2);

Ndata = numel(varargin);
Nfreq = zeros(1,Ndata);
faxis = zeros(1,0);
for i=1:Ndata
  Nfreq(i) = numel(varargin{i}.freq);
  faxis    = [faxis;varargin{i}.freq(:)];
end

if strcmp(required, 'unique')
  boolval = numel(unique(faxis))==numel(faxis) && ~all(isnan(faxis));
  % the second condition is included when the freq is set to dummy nan
elseif strcmp(required, 'identical')
  % the number of freq bins needs at least to be the same across
  % inputs
  boolval = all(Nfreq==Nfreq(1));
  if boolval
    % then check whether the axes are equal
    faxis   = reshape(faxis, Nfreq(1), []);
    boolval = all(all(abs(faxis - repmat(faxis(:,1), 1, Ndata))<tol)==1);
    faxis   = faxis(:,1);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [U,S,V] = svdsecon(X,k)
% Input:
% X : m x n matrix
% k : extracts the first k singular values
%
% Output:
% X = U*S*V' approximately (up to k)
%
% Description:
% Does equivalent to svds(X,k) but faster
% Requires that k < min(m,n) where [m,n] = size(X)
% This function is useful if k is much smaller than m and n
% or if X is sparse (see doc eigs)
%
% Vipin Vijayan (2014)

%X = bsxfun(@minus,X,mean(X,2));
[m,n] = size(X);
assert(k <= m && k <= n, 'k needs to be smaller than size(X,1) and size(X,2)');

if  m <= n
    C = double(X*X');
    [U,D] = eigs(C,k);
    clear C;
    if nargout > 2
        V = X'*U;
        s = sqrt(abs(diag(D)));
        V = bsxfun(@(x,c)x./c, V, s');
        S = diag(s);
    end
else
    C = double(X'*X);
    [V,D] = eigs(C,k);
    clear C;
    U = X*V; % convert evecs from X'*X to X*X'. the evals are the same.
    %s = sqrt(sum(U.^2,1))';
    s = sqrt(abs(diag(D)));
    U = bsxfun(@(x,c)x./c, U, s');
    S = diag(s);
end

author: zhoupc


text:
% (this is just an interface)
function varargout = mydatestr_inv (varargin)
    [varargout{1:nargout}] = mydatestri (varargin{:});
end

%!test
%! % just check whether it runs:
%! mydatestr_inv('2000', 'YYYY');


text:
function obj_val = test_binary_classifier(objective_function,classf, ...
				  prior,system,input_data)
% Returns the result of the objective function evaluated on the
% scores.
%
% Inputs:
%   objective_function: a function handle to the objective function
%                       to feed the scores into
%   classf: length T vector where T is the number of trials with entries +1 for target scores; -1
%           for non-target scores
%   prior: the prior (given to the system that produced the scores)
%   system: a function handle to the system to be run
%   input_data: the data to run the system on (to produce scores)
%
% Outputs
%   obj_val: the value returned by the objective function

if nargin==0
    test_this();
    return;
end

scores = system(input_data);
obj_val = evaluate_objective(objective_function,scores,classf,prior);

end

function test_this()

num_trials = 100;
input_data = randn(20,num_trials);

prior = 0.5;
maxiters = 1000;
classf = [ones(1,num_trials/2),-ones(1,num_trials/2)];
tar = input_data(:,1:num_trials/2);
non = input_data(:,num_trials/2+1:end);
[sys,run_sys,w0] = linear_fusion_factory(tar,non);

w = train_binary_classifier(@cllr_obj,classf,sys,[],w0,[],maxiters,[],prior,[],true);

system = @(data) run_sys(w,data);
test_binary_classifier(@cllr_obj,classf,prior,system,input_data)

end

author: nesl


text:
function main_cartpole_closedloop
close all;

T = 3;
N = 60;

solver = ocl.acados.Solver( ...
  T, ...
  'vars', @ocl.examples.cartpole.vars, ...
  'dae', @ocl.examples.cartpole.dae, ...
  'pathcosts', @ocl.examples.cartpole.pathcosts, ...
  'terminalcost', @ocl.examples.cartpole.terminalcost, ...
  'N', N, 'verbose', false);

solver.setInitialState('p', 0);
solver.setInitialState('v', 0);
solver.setInitialState('theta', 0);
solver.setInitialState('omega', 0);

solver.initialize('theta', [0 1], [0 0]);

solver.solve();

x0 = [0;0;0;0];

sim = ocl.Simulator(@ocl.examples.cartpole.vars, @ocl.examples.cartpole.dae);
sim.reset(x0);

draw_handles = ocl.examples.cartpole.draw_prepare(x0(1), x0(2), 0.8, 6);

% log window
log_fig = figure('menubar', 'none');
log_window = uicontrol(log_fig, 'Style', 'listbox', ...
  'Units',    'normalized', ...
  'Position', [0,0,1,1]);

data = struct;
data.T = T;
data.dt = T/N;
data.draw_handles = draw_handles;
data.t = 0;
data.current_state = x0;
data.force = {};

% control loop
control_timer = timer('TimerFcn', @(t,d) controller(t, d, solver, sim, log_window), ...
  'ExecutionMode', 'fixedRate', 'Period', data.dt, 'UserData', data);

start(control_timer);
cli(control_timer);
stop(control_timer);

end

function controller(t, ~, solver, sim, log_window)

dt = t.UserData.dt;
draw_handles = t.UserData.draw_handles;
time = t.UserData.t;
current_state = t.UserData.current_state;

solver.setInitialState('p', current_state(1));
solver.setInitialState('theta', current_state(2));
solver.setInitialState('v', current_state(3));
solver.setInitialState('omega', current_state(4));

[sol,~] = solver.solve();

u = sol.controls.F.value;

sim.current_state = current_state;

force = 0;
if ~isempty(t.UserData.force)
  force = t.UserData.force{end};
  t.UserData.force(end) = [];
end

x = sim.step(u(1) + force, dt);
x = x.value;

if abs(x(1)) > 6
  sim.current_state = [0;0;0;0];
  solver.initialize('p', [0 1], [0 0]);
  solver.initialize('theta', [0 1], [0 0]);
  solver.initialize('v', [0 1], [0 0]);
  solver.initialize('omega', [0 1], [0 0]);
  solver.initialize('F', [0 1], [0 0]);
end

% draw
ocl.examples.cartpole.draw(draw_handles, time, x, 0.8);

lines = splitlines(solver.stats());
set(log_window, 'String', lines);
drawnow

t.UserData.sol = sol;
t.UserData.t = time + dt;
t.UserData.current_state = sim.current_state;

end

function cli(t)

  cli_info = [newline, 'You are now in the command line interface. ', newline, ...
    'You can make use of ', ...
    'the following commands:', newline, newline, ...
    'q: quits the program', newline, ...
    'f: applies a force of 30N', newline, ...
    '<number>:  applies a force of <number> Newton', newline];

  disp(cli_info);

  terminated = false;
  while ~terminated
    m = input('  cli >>', 's');

    if strcmp(m, 'q') || strcmp(m, 'x') || strcmp(m, 'c')
      disp('exiting...')
      stop(t);
      terminated = true;
    elseif strcmp(m, 'f')
      disp('force!!')
      t.UserData.force{end+1} = 30*sign(rand-0.5);
    elseif ~isnan(str2double(m))
      F = str2double(m);
      disp(['force ', m, '!!!']);
      t.UserData.force{end+1} = F*sign(rand-0.5);
    else
      disp('Command not recognized!')
      disp(cli_info);
    end
  end

end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function [UV,F, res, edge_norms] = ...
  create_irregular_grid_with_min_angle( ...
    xRes, yRes, n, xWrap, yWrap, min_angle, ...
    no_interior, dart_threshold)
  % Deprecated.
  %
  % See also: create_irregular_grid

  warning('Deprecated. Call create_irregular_grid directly');
  [UV,F, res, edge_norms] = create_irregular_grid( ...
    xRes,yRes, ...
    'NumDarts',n,'xWrap',xWrap,'yWrap',yWrap, ...
    'MinAngle',min_angle,'NoInterior',no_interior, ...
    'DartThreshold',dart_threshold);
end

author: alecjacobson

repo: gptoolbox


text:
function [Hn,Xn] = myHist(Data)

nBins = round(length(Data) / 15);

[H,X] = hist(Data, nBins);

[Max,IMax] = max(H);

sumH = sum(H);
curSum = 0.0;

T = 0.990;

i1 = IMax;
i2 = IMax;
I1 = i1;
I2 = i2;
stop1 = 0;
stop2 = 0;
while ((i1>=1) || (i2<=length(X)))
    if (i1~=i2)
        if (stop1~=1) curSum = curSum + H(i1); end
        if (stop2~=1) curSum = curSum + H(i2); end
    else
        curSum = curSum + H(i1);
    end
    if (curSum <= sumH * T)
        I1 = i1;
        I2 = i2;
        if (i1>1)
            i1 = i1 - 1;
        else
            stop1 = 1;
        end
        if (i2<length(X))
            i2 = i2 + 1;
        else
            stop2 = 1;
        end
    else
        break;
    end
end

Xn = X(I1:I2);
Hn = H(I1:I2);

%plot(Xn,Hn + max(H));
%hold on;
%plot(X,H,'r')


text:
function devInfo = FindImaqDevices(desiredFormat)
%FINDIMAQDEVICES returns list of available image acquisition devices.
%   DEVINFO = FINDIMAQDEVICES returns a structure array with named fields.
%     .NAME is good for populating a GUI for user selection.
%     .CONSTRUCTOR aids programmatic connection to each device.
%     .ADAPTOR and .DEVICE help with nondefault object constructors.
%   The length of DEVINFO indicates how many devices detected.
%   If no devices detected, DEVINFO is empty.
%
%   DEVINFO = FINDIMAQDEVICES(DESIREDFORMAT) allows specification of the
%   desired video format. This string is inserted into the constructor
%   field for each device. DESIREDFORMAT is not checked, so if you specify
%   an invalid format then evaluating the constructor(s) will error.
%
%   NOTE: This function calls IMAQRESET, which closes and deletes any
%   previously active videoinput objects.
%
%   See also FINDCAM, IMAQHWINFO, IMAQRESET.

% Copyright 2003-2010 The MathWorks, Inc.

error(nargchk(0,1,nargin))
error(nargoutchk(1,1,nargout))

%make sure recently added plug-n-play devices also detectable
imaqreset

%installed adaptors
hwInfo = imaqhwinfo;
adaptors = hwInfo.InstalledAdaptors;

%available devices
numDevices = 0;  %start with empty list
for i=1:length(adaptors)  %check each adaptor
  adaptorInfo = imaqhwinfo(adaptors{i});
  devices = adaptorInfo.DeviceIDs;
  %add each device to list
  for j=1:length(devices)
    numDevices = numDevices+1;
    deviceInfo = imaqhwinfo(adaptors{i},devices{j});
    devInfo(numDevices).name = deviceInfo.DeviceName;
    constructor = deviceInfo.ObjectConstructor;
    if nargin>0  %add format to object constructor
      constructor = [constructor(1:end-1) ', ''' desiredFormat ''')'];
    end
    devInfo(numDevices).constructor = constructor;
    devInfo(numDevices).adaptor = adaptors{i};
    devInfo(numDevices).device = devices{j};
  end  %for j = devices
end  %for i = adaptors

%return empty idicator if no devices detected
if numDevices==0
  devInfo=[];
end


text:
function  [ons_secs, VOLLOCS, LOCS, verbose] = tapas_physio_create_scan_timing(...
    log_files, scan_timing, ons_secs, verbose)
% Extracts slice and volume scan onsets (triggers) from different vendor formats
%
%   [ons_secs, VOLLOCS, LOCS, verbose] = tapas_physio_create_scan_timing(...
%            log_files, scan_timing, ons_secs, verbose);
%
%
% IN
%   NOTE: The detailed description of all input structures can be found as
%   comments in tapas_physio_new
%
%   log_files    - file names (physiology and scan timing) and sampling rates
%
%   ons_secs     -  structure for time-dependent variables, i.e. onsets,
%                   specified in seconds, in particular
%                   .t          - time vector of phys time course
%
%   scan_timing         -  Parameters for sequence timing & synchronization
%   scan_tming.sqpar    -  slice and volume acquisition starts, TR,
%                          number of scans etc.
%   scan_timing.sync    -  synchronization options
%                          (e.g. from gradients, trigger, tics phys
%                           logfile to scan acquisition)
%
%       sqpar           - sequence timing parameters, used for computation
%                         of scan events from 'nominal' timing
%           .Nslices        - number of slices per volume in fMRI scan
%           .TR             - repetition time in seconds
%           .Ndummies       - number of dummy volumes
%           .Nscans         - number of full volumes saved (volumes in nifti file,
%                             usually rows in your design matrix)
%           .Nprep          - number of non-dummy, volume like preparation pulses
%                             before 1st dummy scan. If set, logfile is read from beginning,
%                             otherwise volumes are counted from last detected volume in the logfile
%           .time_slice_to_slice - time between the acquisition of 2 subsequent
%                             slices; typically TR/Nslices or
%                             minTR/Nslices, if minimal temporal slice
%                             spacing was chosen
%
%   verbose                 - defines output level (which graphics to plot
%                             and whether to save them)
%
% OUT
%   ons_secs    -  structure for time-dependent variables, i.e. onsets,
%                  specified in seconds, updated fields
%                   .spulse     - scan slice trigger events
%                   .svolpulse  - scan volume trigger events
%                   .spulse_per_vol
%                               - cell(nVolumes,1) of slice triggers per
%                                 volume
%                   .acq_codes  - acquisition codes (e.g. triggers) within
%                                 phys log files (e.g. Philips, Biopac)
%
%   VOLLOCS     - index locations in time vector (of physiological recordings),
%                             when volume scan events started
%   LOCS        - locations in time vector, when slice or volume scan
%                             events started
%
%   See also tapas_physio_new tapas_physio_main_create_regresssors

% Author: Lars Kasper
% Created: 2013-08-23
% Copyright (C) 2016 Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

sqpar   = scan_timing.sqpar;

% TODO: introduce auto that takes time stamps from default locations
% for different vendors
switch lower(scan_timing.sync.method)
    case 'nominal'
        [VOLLOCS, LOCS] = ...
            tapas_physio_create_scan_timing_nominal(ons_secs.t + ...
                    ons_secs.t_start, sqpar, log_files.align_scan);
    case {'gradient', 'gradient_log'}
        [VOLLOCS, LOCS, verbose] = ...
            tapas_physio_create_scan_timing_from_gradients_philips( ...
            log_files, scan_timing, verbose);
    case {'gradient_auto', 'gradient_log_auto'}
        [VOLLOCS, LOCS, verbose] = ...
            tapas_physio_create_scan_timing_from_gradients_auto_philips( ...
            log_files, scan_timing, verbose);
    case 'scan_timing_log'
        switch lower(log_files.vendor)
            case 'siemens'
                % for alignScan = 'last', in case logfile lasts longer than end of last scan
                % assuming t = 0 is already start of first scan
                durationPhyslogAfterEndOfLastScan = ons_secs.t(end) + ...
                    ons_secs.t_start - sqpar.Nscans*sqpar.TR;
                [VOLLOCS, LOCS] = ...
                    tapas_physio_create_scan_timing_nominal(ons_secs.t + ...
                    ons_secs.t_start, sqpar, log_files.align_scan, ...
                    durationPhyslogAfterEndOfLastScan);
            case 'siemens_tics'
                [VOLLOCS, LOCS, verbose] = ...
                    tapas_physio_create_scan_timing_from_tics_siemens( ...
                    ons_secs.t, ons_secs.t_start, log_files, verbose);
            case {'biopac_mat', 'biopac_txt', 'bids'}
                [VOLLOCS, LOCS, verbose] = ...
                    tapas_physio_create_scan_timing_from_acq_codes( ...
                    ons_secs.t + ons_secs.t_start, ons_secs.acq_codes, ...
                    sqpar, log_files.align_scan, verbose);
        end
    otherwise
        verbose = tapas_physio_log(...
            sprintf('unknown scan_timing.sync.method: %s', ...
            scan_timing.sync.method), verbose, 2);
end

% remove arbitrary offset in time vector now, since all timings have now
% been aligned to ons_secs.t
% ons_secs.t = ons_secs.t - ons_secs.t(1);

[ons_secs.svolpulse, ons_secs.spulse, ons_secs.spulse_per_vol, verbose] = ...
    tapas_physio_get_onsets_from_locs(...
    ons_secs.t, VOLLOCS, LOCS, sqpar, verbose);

author: translationalneuromodeling

repo: tapas


text:
% For transforming adjacency matrices to random orders,
% need dataset saved as "name_random.mat".

dataset = strvcat('MUTAG_random');
for ith_data = 1: size(dataset, 1)
    clear save_struct
    data = strcat(dataset(ith_data, :))
    dataname = strcat('../data/raw_data/', data, '.mat');
    load(dataname);
    X = eval(data(1:end-7));
    for i = 1: size(X, 2)
        am = full(X(i).am);
        l = X(i).nl.values;
        %[~, ~, am, l] = palette_wl(am, l); % transform to wl vertex order
        perm = randperm(length(l));
        X(i).am = am(perm, perm);
        X(i).nl.values = l(perm);
    end
    save_struct.(data(1:end-7)) = X;
    save(dataname, '-struct', 'save_struct', '-append');
end

author: muhanzhang

repo: DGCNN

save_path: github-repos/MATLAB/muhanzhang-DGCNN


text:
function varargout = circleArcAsCurve(arc, N)
%CIRCLEARCASCURVE Convert a circle arc into a series of points.
%
%   deprecated: replaced by circleArcToPolyline
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2006-05-22
% Copyright 2006 INRA - Cepia Software Platform

warning('matGeom:deprecated', ...
    'function "circleArcAsCurve" is deprecated, use "circleArcToPolygon" instead');

% format output
if nargout <= 1
    varargout = {circleArcToPolyline(arc, N)};
else
    [x, y] = circleArcToPolyline(arc, N);
    varargout = {x, y};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%--- help for dsge/bvar_dsge ---
%
%  INTERNAL FUNCTION: intermediary file for computing key elements for dsge-var
%
%  ::
%
%    [obj,retcode] = bvar_dsge(obj,varargin)
%
%  Args:
%
%     obj (rise | dsge): model object
%     varargin : ususal dsge options. The most important of which are:
%
%       - **dsgevar_lag** [integer|{4}]: number of lags in the VAR
%       - **dsgevar_constant** [false|{true}]: flag for having a constant in
%         the VAR
%       - **dsgevar_var_regime** [false|{true}]: use the VAR in simulations.
%         Otherwise use the DSGE
%       - **dsgevar_inner_param_uncertainty** [true|{false}]: make random draws
%         for the parameters around the mode for each simulation
%
%  Returns:
%     :
%
%     - **obj** [rise|dsge]: model object
%
%  Note:
%
%     - Because the BVAR-DSGE will have fewer variables than the DSGE, in
%       simulations, the missing variables will be stored as 0+1i in time series.
%       This is true for forecast, irf and simulate
%
%

author: jmaih


text:
function N = uplus(N)
%UPLUS   Unitary plus for CHEBOP2 objects.

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

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

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

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

% Last Modified by GUIDE v2.5 20-Sep-2013 16:08:24

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

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

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

Figure_handles.MU_display2=hObject;
filepath=mfilename('fullpath');   % find current running m-file fullpath
sep=filesep;
k=strfind(filepath, sep);
path=filepath(1:k(end)-1);

MU_icon(handles.Sag2Main_pushbutton,[path sep '..' sep '..' sep 'Resource' sep 'Icon' sep 'Exchange.gif']);
MU_icon(handles.Cor2Main_pushbutton,[path sep '..' sep '..' sep 'Resource' sep 'Icon' sep 'Exchange.gif']);

% Choose default command line output for MU_Matrix_Display_3Plane
handles.Parent = Figure_handles.MU_display;
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes MU_Matrix_Display_3Plane wait for user response (see UIRESUME)
% uiwait(MU_display_handles.MU_matrix_display_3Plane);

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

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

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

MU_display_handles=guidata(handles.Parent);
if MU_display_handles.V.Localizer.Local_switch==0;
    MU_display_handles.V.Localizer.Local_switch=1;
    MU_display_handles=MU_update_image(MU_display_handles.Matrix_display_axes,{MU_display_handles.TMatrix,MU_display_handles.Mask},MU_display_handles,0);
else
    MU_display_handles.V.Localizer.Local_switch=0;
    delete(MU_display_handles.V.Localizer.Local_h1);
    delete(MU_display_handles.V.Localizer.Local_h2);
end
guidata(handles.Parent, MU_display_handles);

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

MU_display_handles=guidata(handles.Parent);
if MU_display_handles.V.Localizer.Local_switch==0;
    MU_display_handles.V.Localizer.Local_switch=1;
    MU_display_handles=MU_update_image(MU_display_handles.Matrix_display_axes,{MU_display_handles.TMatrix,MU_display_handles.Mask},MU_display_handles,0);
else
    MU_display_handles.V.Localizer.Local_switch=0;
    delete(MU_display_handles.V.Localizer.Local_h1);
    delete(MU_display_handles.V.Localizer.Local_h2);
end
guidata(handles.Parent, MU_display_handles);

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

MU_display_handles=guidata(handles.Parent);

% initialize MU_display & more
MergeM=get(MU_display_handles.Matrix_name_edit,'String');
MU_display_new_handles=MU_display_initialize(MU_display_handles.MU_matrix_display, ...
                                             MU_display_handles,[MergeM '_s2m'], ...
                                             MU_display_handles.SMatrix);
if isempty(MU_display_new_handles)
    return;
end

MU_display_new_handles.Mask=permute(MU_display_handles.Mask,[3 1 2]);
MU_display_new_handles.SMatrix=permute(MU_display_handles.SMatrix,[3 1 2]); % Sagittal Matrix
MU_display_new_handles.CMatrix=permute(MU_display_handles.SMatrix,[3 2 1]); % Coronal Matrix
MU_display_new_handles.Slicer=1;
MU_display_new_handles.V.Localizer=MU_display_handles.V.Localizer;
MU_display_new_handles.V.Localizer.Local_point=[ceil(MU_display_new_handles.V.DimSize(2)/2) ...
                                                ceil(MU_display_new_handles.V.DimSize(1)/2)];

MU_display_new_handles=MU_update_image(MU_display_new_handles.Matrix_display_axes, ...
                                      {MU_display_new_handles.TMatrix,MU_display_new_handles.Mask},...
                                       MU_display_new_handles,0);
cla(handles.Matrix_display_axes2);
cla(handles.Matrix_display_axes3);
MU_update_ass_image({handles.Matrix_display_axes2,handles.Matrix_display_axes3}, ...
                    {MU_display_new_handles.SMatrix,MU_display_new_handles.CMatrix}, ...
                    MU_display_new_handles);

guidata(handles.Parent, MU_display_new_handles);

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

MU_display_handles=guidata(handles.Parent);

% initialize MU_display & more
MergeM=get(MU_display_handles.Matrix_name_edit,'String');
MU_display_new_handles=MU_display_initialize(MU_display_handles.MU_matrix_display, ...
                                             MU_display_handles,[MergeM '_c2m'], ...
                                             MU_display_handles.CMatrix);
if isempty(MU_display_new_handles)
    return;
end

MU_display_new_handles.Mask=permute(MU_display_handles.Mask,[3 2 1]);
MU_display_new_handles.SMatrix=permute(MU_display_handles.CMatrix,[3 1 2]); % Sagittal Matrix
MU_display_new_handles.CMatrix=permute(MU_display_handles.CMatrix,[3 2 1]); % Coronal Matrix
MU_display_new_handles.Slicer=1;
MU_display_new_handles.V.Localizer=MU_display_handles.V.Localizer;
MU_display_new_handles.V.Localizer.Local_point=[ceil(MU_display_new_handles.V.DimSize(2)/2) ...
                                                ceil(MU_display_new_handles.V.DimSize(1)/2)];

MU_display_new_handles=MU_update_image(MU_display_new_handles.Matrix_display_axes, ...
                                      {MU_display_new_handles.TMatrix,MU_display_new_handles.Mask},...
                                       MU_display_new_handles,0);
cla(handles.Matrix_display_axes2);
cla(handles.Matrix_display_axes3);
MU_update_ass_image({handles.Matrix_display_axes2,handles.Matrix_display_axes3}, ...
                    {MU_display_new_handles.SMatrix,MU_display_new_handles.CMatrix}, ...
                    MU_display_new_handles);

guidata(handles.Parent, MU_display_new_handles);

% --- Executes when user attempts to close MU_matrix_display_3Plane.
function MU_matrix_display_3Plane_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to MU_matrix_display_3Plane (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

global Figure_handles
Figure_handles=rmfield(Figure_handles,'MU_display2');

if ~ishandle(handles.Parent)
    delete(hObject);
    return;
end
MU_display_handles=guidata(handles.Parent);
MU_display_handles.V.Localizer.Local_switch=0;
MU_display_handles.Slicer=0;
MU_display_handles=rmfield(MU_display_handles,'SMatrix');
MU_display_handles=rmfield(MU_display_handles,'CMatrix');

try
    delete(MU_display_handles.V.Localizer.Local_h1);
    delete(MU_display_handles.V.Localizer.Local_h2);
catch me
end
set(MU_display_handles.MatrixCalc_pushbutton,'Enable','on'); % enable matrix calculation when slicer if off
guidata(handles.Parent, MU_display_handles);

% Hint: delete(hObject) closes the figure
delete(hObject);

% --- Executes on mouse motion over figure - except title and menu.
function MU_matrix_display_3Plane_WindowButtonMotionFcn(hObject, eventdata, handles)
% hObject    handle to MU_matrix_display_3Plane (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if ~ishandle(handles.Parent)
    delete(hObject);
end

% --- Executes during object deletion, before destroying properties.
function MU_matrix_display_3Plane_DeleteFcn(hObject, eventdata, handles)
% hObject    handle to MU_matrix_display_3Plane (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [panel, Vol] = DeterminePanelGeometry(inputgeo, figs)
% find coordinates for horseshoe vortices and control points and plot
% Aircraft-Fixed coordinate system:
%       x: forward along fuselage axis
%       y: out starboard wing
%       z: down
% Local Profile coordinate system:
%       x: along chord from leading edge toward trailing edge
%       z: up

plot_on = 1;

%% constants
halfSpan=inputgeo.b/2;     % half wingspan
sweep=inputgeo.sweep;       % wing sweep angle in radians
dihedral=inputgeo.dih;    %dihedral angle in radians
nc=inputgeo.nc;  %number of panels on camber line, upper and lower surfaces
ns = inputgeo.ns; % number of span segments per side
M = inputgeo.M;    %Freestream mach number
beta = sqrt(1-M^2); %Prandtl-Glauert correction

%% Root Airfoil data
y=0;            % span station
chord=inputgeo.c_r;        % chord length
alphaRoot=inputgeo.i_r;    % geometric angle of incidence at root in radians

% call function
[Croot,Troot,Aroot]=NacaCoord(inputgeo.root,y,nc); % Croot = nondimensional camber line, Troot = nondimensional surface;
%Croot coordinates = [x location in fraction chord, y location in fraction
%                       chord, z location in fraction chord]
%Troot coordinates = [x location in fraction chord, y location in fraction
%                       chord, z location in fraction chord]

%% Scale to chord length
Croot(:,1)=chord*Croot(:,1);Croot(:,3)=chord*Croot(:,3);
Troot(:,1)=chord*Troot(:,1);Troot(:,3)=chord*Troot(:,3);
Aroot=Aroot*chord^2;

%Calculate slope of camber line for each panel at the root
%Calculate the chord along each elemental panel
dzdxRoot = zeros(1,nc);
ccRoot = zeros(1,nc);
if nc == 1 %If there is only one chordwise element
    i=1;
    dzdxRoot(i)=-(Croot(i+1,3)-Croot(i,3))/(Croot(i+1,1)-Croot(i,1));
    ccRoot(i)=0.75*(Croot(i+1,1)-Croot(i,1));
else        %if there are more than one chordwise elements
    for i=1:nc
        dzdxRoot(i)=(Croot(i+1,3)-Croot(i,3))/(Croot(i+1,1)-Croot(i,1));
        if i==nc
            ccRoot(i)=0.75*(Croot(i+1,1)-Croot(i,1));
        else
            ccRoot(i)=0.75*(Croot(i+1,1)-Croot(i,1))+0.25*(Croot(i+2,1)-Croot(i+1,1));
        end
    end
end

if plot_on
    axes(figs.root); cla;
    plot(Croot(:,1)*-3.2808,Croot(:,3)*-3.2808,'r') %Convert to feet; plot with z-axis positive up and x-axis to the left
    hold on
    plot(Troot(:,1)*-3.2808,Troot(:,3)*-3.2808,'b')  %Convert to feet; plot with z-axis positive up and x-axis to the left
    axis equal
    title('Root airfoil')
end

%% Tip Airfoil data
y=-halfSpan;     % span station; left wing
chord=inputgeo.taper*inputgeo.c_r;        % chord length
alphaTip=inputgeo.i_r+inputgeo.twist;    % geometric angle of incidence at tip in radians

%call function
[Ctip,Ttip,Atip]=NacaCoord(inputgeo.tip,y,nc);

%% scale to chord length
Ctip(:,1)=chord*Ctip(:,1);Ctip(:,3)=chord*Ctip(:,3);
Ttip(:,1)=chord*Ttip(:,1);Ttip(:,3)=chord*Ttip(:,3);
Atip=Atip*chord^2;

%Calculate slope of camber line for each panel at the tip
%Calculate the chord along each elemental panel
dzdxTip = zeros(1,nc);
ccTip = zeros(1,nc);
if nc == 1 %If there is only one chordwise element
    dzdxTip(i)=-(Ctip(i+1,3)-Ctip(i,3))/(Ctip(i+1,1)-Ctip(i,1));
    ccTip(i)=0.75*(Ctip(i+1,1)-Ctip(i,1));
else        %if there are more than one chordwise elements
    for i=1:nc
        dzdxTip(i)=(Ctip(i+1,3)-Ctip(i,3))/(Ctip(i+1,1)-Ctip(i,1));
        if i==nc
            ccTip(i)=0.75*(Ctip(i+1,1)-Ctip(i,1));
        else
            ccTip(i)=0.75*(Ctip(i+1,1)-Ctip(i,1))+0.25*(Ctip(i+2,1)-Ctip(i+1,1));
        end
    end
end

if plot_on
    axes(figs.tip); cla;
    plot(Ctip(:,1)*-3.2808,Ctip(:,3)*-3.2808,'r') %Convert to feet; plot with z-axis positive up and x-axis to the left
    hold on
    plot(Ttip(:,1)*-3.2808,Ttip(:,3)*-3.2808,'b') %Convert to feet; plot with z-axis positive up and x-axis to the left
    axis equal
    title('Tip airfoil')
end

%% Root: Transform to Aircraft-Fixed Coordinate System:
% skin:
%Skin(number of spanwise sections, number of chordwise sections, coordinates in 3 dimensions (x,y,z))
Skin(1,1:2*nc+1,1:3) = Troot*[cos(alphaRoot) 0 sin(alphaRoot); 0 1 0; -sin(alphaRoot) 0 cos(alphaRoot)];

% Mean camber line:
Camber(1,1:nc+1,1:3) = Croot*[cos(alphaRoot) 0 sin(alphaRoot); 0 1 0; -sin(alphaRoot) 0 cos(alphaRoot)];

%% Tip: Transform to Aircraft-Fixed Coordinate System:
% sweep and dihedral
% Negative sign occurs because local axes are opposite to global
xLETip = -halfSpan*(tan(sweep));% x coordinate of Leading Edge at Tip
zLETip = -halfSpan*(tan(dihedral));% z coordinate of Leading Edge at Tip
% Skin:
LETip=ones(2*nc+1,1)*[xLETip 0 zLETip];
Skin(ns+1,1:2*nc+1,1:3) = Ttip*[cos(alphaTip) 0 sin(alphaTip); 0 1 0; -sin(alphaTip) 0 cos(alphaTip)]+LETip;

% Mean camber line
LETip=ones(nc+1,1)*[xLETip 0 zLETip];
Camber(ns+1,1:nc+1,1:3) = Ctip*[cos(alphaTip) 0 sin(alphaTip); 0 1 0; -sin(alphaTip) 0 cos(alphaTip)]+LETip;

%% Intermediate stations
for k=1:ns-1 % k = 0 corresponds to root, k = ns corresponds to tip
    eta=k/ns;
    for i=1:nc+1 % along camber
        Camber(k+1,i,:) = Camber(1,i,:)+eta*(Camber(ns+1,i,:)-Camber(1,i,:));
    end
    for i=1:2*nc+1 %along surface
        Skin(k+1,i,:) = Skin(1,i,:)+eta*(Skin(ns+1,i,:)-Skin(1,i,:));
    end
end

if plot_on
    %Plot wing surface
    axes(figs.fig); cla; hold on; axis equal;
    for k=1:ns+1
        plot3(Skin(k,:,1)*-3.2808,Skin(k,:,2)*3.2808,Skin(k,:,3)*-3.2808);  %Convert to feet; plot with z-axis positive up and x-axis to the left
        plot3(Skin(k,:,1)*-3.2808,-Skin(k,:,2)*3.2808,Skin(k,:,3)*-3.2808);  %Convert to feet; plot with z-axis positive up and x-axis to the left
        plot3(Camber(k,:,1)*-3.2808,Camber(k,:,2)*3.2808,Camber(k,:,3)*-3.2808,'r')  %Convert to feet; plot with z-axis positive up and x-axis to the left
    end
    for i=1:size(Skin,2)
    plot3(Skin(:,i,1)*-3.2808,Skin(:,i,2)*3.2808,Skin(:,i,3)*-3.2808);
    plot3(Skin(:,i,1)*-3.2808,-Skin(:,i,2)*3.2808,Skin(:,i,3)*-3.2808);
    end
end

%Determine panel properties
twist = zeros(1,ns);
dzdx = zeros(ns,nc);
cc = zeros(ns,nc);
CP = zeros(ns,nc,3);
BV = zeros(ns,nc,3);
BV1 = zeros(ns,nc,3);
BV2 = zeros(ns,nc,3);
BVhalfspan = zeros(ns,nc);
sweep_c4 = zeros(ns,nc);

for k = 1:ns
    %Determine local twist at each panel, dz/dx at control points, and
    %chord along left trailing leg of elemental panel
    eta=k/(ns+1);
    twist(k)=alphaRoot+eta*(alphaTip-alphaRoot);
    dzdx(k,:)=dzdxRoot+eta*(dzdxTip-dzdxRoot);
    cc(k,:)=ccRoot+eta*(ccTip-ccRoot);  %Validated
    %Determine control point coordinates, coordinates to center of bound
    %vortex, half span of bound vortex, and quarter-chord sweep angle
    for i = 1:nc
        CP(k,i,:)=0.5*(Camber(k,i,:)+Camber(k+1,i,:))+0.75*(0.5*(Camber(k,i+1,:)+Camber(k+1,i+1,:))-0.5*(Camber(k,i,:)+Camber(k+1,i,:)));
        BV(k,i,:)=0.5*(Camber(k,i,:)+Camber(k+1,i,:))+0.25*(0.5*(Camber(k,i+1,:)+Camber(k+1,i+1,:))-0.5*(Camber(k,i,:)+Camber(k+1,i,:)));
        BV1(k,i,:)=0.75*Camber(k+1,i,:)+0.25*Camber(k+1,i+1,:);  %Left bound vortex coordinate
        BV2(k,i,:)=0.75*Camber(k,i,:)+0.25*Camber(k,i+1,:);  %Right bound vortex coordinate
        BVhalfspan(k,i)= 0.5*(0.75*norm(shiftdim(Camber(k,i,2:3)-Camber(k+1,i,2:3)))+0.25*norm(shiftdim(Camber(k,i+1,2:3)-Camber(k+1,i+1,2:3))));
        %Only consider [Y,Z] distances due to definition of halfspan in
        %NASA paper; considering X distance as well makes halfspan increase
        %greatly as you increase sweep and invalidates calculations
        %Validated BVhalfspan calculation through (1:ns-1,1:nc)
        sweep_c4(k,i)=atan((Camber(k,i,1)-Camber(k+1,i,1))/(Camber(k,i,2)-Camber(k+1,i,2))); %Validated
    end
end

%Prandtl-Glauert corrections
CPprime = CP(:,:,1)/beta;
BVprime = BV(:,:,1)/beta;
BV1prime = BV1(:,:,1)/beta;
BV2prime = BV2(:,:,1)/beta;
sweep_c4_prime = atan(tan(sweep_c4)/beta);

for k = 1:ns
    for i  = 1:nc
        panel(k,i).CP=shiftdim(CP(k,i,:));  %Removing singleton dimensions
        panel(k,i).BV=shiftdim(BV(k,i,:)); %Removing singleton dimensions
        panel(k,i).BV1=shiftdim(BV1(k,i,:)); %Removing singleton dimensions
        panel(k,i).BV2=shiftdim(BV2(k,i,:)); %Removing singleton dimensions
        panel(k,i).dzdx=dzdx(k,i);
        panel(k,i).twist=twist(k);
        panel(k,i).s=BVhalfspan(k,i);
        panel(k,i).sweep=sweep_c4(k,i);
        panel(k,i).sweepprime=sweep_c4_prime(k,i);
        panel(k,i).CPprime=CPprime(k,i);
        panel(k,i).BVprime=BVprime(k,i);
        panel(k,i).BV1prime=BV1prime(k,i);
        panel(k,i).BV2prime=BV2prime(k,i);
        panel(k,i).cc=cc(k,i);
    end
end

%Determine wing volume
Vol = 1/3*(Aroot +sqrt(Aroot*Atip) + Atip)*inputgeo.b;


text:
function [fieldNames, fieldValues] = get_nonempty_fields(obj)
% Returns field names and values of all empty fields of MrCopyData-object
%
%   Y = MrCopyData()
%   Y.get_nonempty_fields(inputs)
%
% This is a method of class MrCopyData.
%
% IN
%
% OUT
%
% EXAMPLE
%   Y.get_nonempty_fields
%
%   See also MrCopyData

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-12-09
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.
%

fieldNames = cell(0,1);
fieldValues = cell(0,1);

%% Find all properties of MrCopyData to be searched
[sel, mobj] = get_properties_to_update(obj);

%% Loop over all properties distinguishing between variables and MrCopyData
% -objects which are treated recursively

for k = sel(:)'
    pname = mobj.Properties{k}.Name;
    currProp = obj.(pname);

    if isa(currProp, 'MrCopyData')
    %% Recursive operation on MrCopyData-property

        [fieldNamesTmp, fieldValuesTmp] = ...
            currProp.get_nonempty_fields();
        fieldNames = [fieldNames; fieldNamesTmp(:)];
        fieldValues = [fieldValues; fieldValuesTmp(:)];

    else
        %% Cell of MrCopyData also treated recursively for each cell element

        if iscell(currProp) ...
                && length(currProp) ...
                && isa(currProp{1}, 'MrCopyData')

            nCellProp = numel(currProp);
            currPropArray = currProp;
            for c = 1:nCellProp
                currProp = currPropArray{c};
                [fieldNamesTmp, fieldValuesTmp] = ...
                    currProp.get_nonempty_fields();
                fieldNames = [fieldNames; fieldNamesTmp(:)];
                fieldValues = [fieldValues; fieldValuesTmp(:)];
            end

        else
        %% No cell of MrCopyData, no MrCopyData...any other property, therefore
        % treat differently, and append output

            if ~isempty(currProp)
                fieldNamesTmp = {pname};
                fieldValuesTmp = {currProp};
                fieldNames = [fieldNames; fieldNamesTmp(:)];
                fieldValues = [fieldValues; fieldValuesTmp(:)];
            end
        end
    end
end

end


text:
function t = tvec_even_bracket3 ( nt, theta1, theta2 )

%*****************************************************************************80
%
%% TVEC_EVEN_BRACKET3 computes an evenly spaced set of angles between THETA1 and THETA2.
%
%  Example:
%
%    NT = 3
%    THETA1 = 30
%    THETA2 = 90
%
%    T = ( 40, 60, 80 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NT, the number of values to compute.
%
%    Input, real THETA1, THETA2, the limiting angles.
%
%    Output, real TVEC(NT), the evenly spaced angles.
%
  for i = 1 : nt
    t(i) = ( ( 2 * nt - 2 * i + 1 ) * theta1   ...
           + (          2 * i - 1 ) * theta2 ) ...
           / ( 2 * nt             );
  end

  return
end


text:
function res = exp_rule(match,repl)
% Create a substitution rule, which to replace a Pattern in an expression.
% Rule = exp_rule(Pattern,Replacement)
%
% In:
%   Pattern     : An expression; may contain blanks (see exp_blank()), named sub-expressions (see
%                 exp_pattern()) or generally any structure (see exp_symbol()).
%   Replacement : a replacement expression; may contain symbols and/or named patterns that reference
%                 parts inside Pattern.
%
% Out:
%   Rule        : a replacement rule
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-23

res = struct('head',{@Rule},'parts',{{match,repl}});

repo: REST


text:
% GET THE UNIT VECTOR
function [v_unit] = unit(v)
    v_unit = v/norm(v);
end

repo: OpenMAS


text:
function [c, ceq] = ma_latConstraint(stateLog, eventID, lbLat, ubLat, bodyIDApply, celBodyData, maData)
%ma_semiMajorAxisConstraint Summary of this function goes here
%   Detailed explanation goes here

    normFact = pi;

    if(ischar(eventID) && strcmpi(eventID,'final'))
        eventNum = max(stateLog(:,13));
    else
        [~, eventNum] = getEventByID(eventID, maData.script);
    end

    eventLog = stateLog(stateLog(:,13)==eventNum,:);
    finalEntry = eventLog(end,:);

    bodyID = finalEntry(8);

    if(bodyID == bodyIDApply || bodyIDApply==-1)
        bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);
        rVect = finalEntry(2:4)';
        vVect = finalEntry(5:7)';
        utSec = finalEntry(1);

        [lat, ~, ~] = getLatLongAltFromInertialVect(utSec, rVect, bodyInfo, vVect);

        lat = rad2deg(lat);
        if(lbLat == ubLat)
            c = [0 0];
            ceq(1) = lat - ubLat;
        else
            c(1) = lbLat - lat;
            c(2) = lat - ubLat;
            ceq = [0];
        end
        c = c/normFact;
        ceq = ceq/normFact;
    else
        c = [0 0];
        ceq = [0];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [hd_record]=tvhdecomp(beta_gibbs,D_record,strshocks_record,It,Bu,Y,n,m,p,k,T)

% preliminary tasks
% first create the hd_record and temp cells
hd_record=cell(n,n+1);
temp=cell(n,2);

% then initiate the Gibbs algorithm
for ii=1:It-Bu

% recover beta for the current iteration (one column for each period)
beta_iter=[];
for jj=1:T
beta_iter=[beta_iter beta_gibbs{jj,1}(:,ii)];
end

% recover D for the current iteration (one page for each period)
D_iter=repmat(reshape(D_record(:,ii),n,n),[1,1,T]);

% recover the series of period-specific orthogonal IRFs
[IRFcell]=bear.tvhdsim(beta_iter,D_iter,n,m,p,T,k);

% recover the structural disturbances
ETA=[];
for jj=1:n
ETA=[ETA;strshocks_record{jj,1}(ii,:)];
end
ETA=ETA';

% then compute the historical decomposition
   % loop over variables
   for jj=1:n
      % loop over shocks
      for kk=1:n
      % loop over shocks
      vshocks=ETA(:,kk);
         % loop over time periods
         for ll=1:T
         % initiate the vectors of IRFs and shocks
         virf=[];
            % then loop over IRF periods
            for mm=1:ll
            % create the vector of IRF coefficients
            virf=[virf IRFcell{mm,ll}(jj,kk)];
            end
         % compute then the contribution of shock kk for variable jj at period ll
         hd_record{jj,kk}(ii,ll)=virf*flipud(vshocks(1:ll,1));
         end
      end
   end

% then go for next Gibbs iteration
end

% step 6: compute the contributions of deterministic variables
% loop over rows of temp/hd_record
for ii=1:n
% fill the Ytot matrix in temp
% initial condition
temp{ii,1}=hd_record{ii,1};
   % sum over the remaining columns of hd_record
   for jj=2:n
   temp{ii,1}=temp{ii,1}+hd_record{ii,jj};
   end
% fill the Y matrix in temp
temp{ii,2}=repmat(Y(:,ii)',It-Bu,1);
% fill the Yd matrix in hd_record
hd_record{ii,n+1}=temp{ii,2}-temp{ii,1};
% go for next variable
end

author: european-central-bank

repo: BEAR-toolbox


text:
#!/usr/bin/env octave
%% Machine Learning Online Class - Exercise 1: Linear Regression

%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  linear exercise. You will need to complete the following functions
%  in this exericse:
%
%     warmUpExercise.m
%     plotData.m
%     gradientDescent.m
%     computeCost.m
%     gradientDescentMulti.m
%     computeCostMulti.m
%     featureNormalize.m
%     normalEqn.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%
% x refers to the population size in 10,000s
% y refers to the profit in $10,000s
%

%% Initialization
clear all; close all; clc

%% ==================== Part 1: Basic Function ====================
% Complete warmUpExercise.m
fprintf('Running warmUpExercise ... \n');
fprintf('5x5 Identity Matrix: \n');
warmUpExercise()

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ======================= Part 2: Plotting =======================
fprintf('Plotting Data ...\n')
data = csvread('ex1data1.txt');
X = data(:, 1); y = data(:, 2);
m = length(y); % number of training examples

% Plot Data
% Note: You have to complete the code in plotData.m
plotData(X, y);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =================== Part 3: Gradient descent ===================
fprintf('Running Gradient Descent ...\n')

X = [ones(m, 1), data(:,1)]; % Add a column of ones to x
theta = zeros(2, 1); % initialize fitting parameters

% Some gradient descent settings
iterations = 1500;
alpha = 0.01;

% compute and display initial cost
computeCost(X, y, theta)

% run gradient descent
theta = gradientDescent(X, y, theta, alpha, iterations);

% print theta to screen
fprintf('Theta found by gradient descent: ');
fprintf('%f %f \n', theta(1), theta(2));

% Plot the linear fit
hold on; % keep previous plot visible
plot(X(:,2), X*theta, '-')
legend('Training data', 'Linear regression')
hold off % don't overlay any more plots on this figure

% Predict values for population sizes of 35,000 and 70,000
predict1 = [1, 3.5] *theta;
fprintf('For population = 35,000, we predict a profit of %f\n',...
    predict1*10000);
predict2 = [1, 7] * theta;
fprintf('For population = 70,000, we predict a profit of %f\n',...
    predict2*10000);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ============= Part 4: Visualizing J(theta_0, theta_1) =============
fprintf('Visualizing J(theta_0, theta_1) ...\n')

% Grid over which we will calculate J
theta0_vals = linspace(-10, 10, 100);
theta1_vals = linspace(-1, 4, 100);

% initialize J_vals to a matrix of 0's
J_vals = zeros(length(theta0_vals), length(theta1_vals));

% Fill out J_vals
for i = 1:length(theta0_vals)
    for j = 1:length(theta1_vals)
	  t = [theta0_vals(i); theta1_vals(j)];
	  J_vals(i,j) = computeCost(X, y, t);
    end
end

% Because of the way meshgrids work in the surf command, we need to
% transpose J_vals before calling surf, or else the axes will be flipped
J_vals = J_vals';
% Surface plot
figure;
surf(theta0_vals, theta1_vals, J_vals)
xlabel('\theta_0'); ylabel('\theta_1');

% Contour plot
figure;
% Plot J_vals as 15 contours spaced logarithmically between 0.01 and 100
contour(theta0_vals, theta1_vals, J_vals, logspace(-2, 3, 20))
xlabel('\theta_0'); ylabel('\theta_1');
hold on;
plot(theta(1), theta(2), 'rx', 'MarkerSize', 10, 'LineWidth', 2);

fprintf('Program paused. Press enter to continue.\n');
pause;

author: SaveTheRbtz

repo: ml-class


text:
function [ xout, yout ] = quaerotate ( xin, yin )

%*****************************************************************************80
%
%% QUAEROTATE applies a rotation.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    26 June 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, real XIN, YIN, the coordinates of the point.
%
%    Output, real XOUT, YOUT, the coordinates of the point
%    after rotation.
%

%
%  Initialize the matrix of rotation.
%
  theta = 2.0 * pi / 3.0;
  a11 = cos ( theta );
  a22 = cos ( theta );
  a12 = - sin ( theta );
  a21 = -a12;
%
%  Apply the rotation matrix to the input vector.
%
  xout = a11 * xin + a12 * yin;
  yout = a21 * xin + a22 * yin;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function h = r8mat_house_form ( n, v )

%*****************************************************************************80
%
%% R8MAT_HOUSE_FORM constructs a Householder matrix from its compact form.
%
%  Discussion:
%
%    H(v) = I - 2 * v * v' / ( v' * v )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real V(N), the vector defining the Householder matrix.
%
%    Output, real H(N,N), the Householder matrix.
%

%
%  Compute the L2 norm of V.
%
  beta = sum ( v(1:n).^2 );
%
%  Form the matrix H.
%
  h = r8mat_identity ( n );

  for i = 1 : n
    for j = 1 : n
      h(i,j) = h(i,j) - 2.0 * v(i) * v(j) / beta;
    end
  end

  return
end


text:
function qr=v_rotmr2qr(mr)
%V_ROTMR2QR converts a matrix of real quaternion matrices to quaternion vectors
% Inputs:
%
%     MR(4m,4n,...)   mxn matrix of real quaternion matrices (each 4x4)
%
% Outputs:
%
%     QR(4m,n,...)   mxn matrix of real quaternion vectors (each 4x1)
%
% In matrix form, quaternions can be multiplied and added using normal matrix
% arithmetic. Each element of an mxn matrix of quaternions is itself a 4x4 block
% so the total dimension of MR is 4m x 4n.

%
%      Copyright (C) Mike Brookes 2000-2018
%      Version: $Id: v_rotmr2qr.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s=size(mr);
s(2)=s(2)/4;
mr=reshape(mr,s(1),[]);
qr=reshape(mr(:,1:4:end),s);
if ~nargout
    qr=qr(1:4); % select the first element
    v_rotqr2ro(qr(:)); % plot a rotated cube
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [alt, lat] = geodet1 (rmag, dec)

% geodetic latitude and altitude

% series solution

% input

%  rmag = geocentric radius (kilometers)
%  dec  = geocentric declination (radians)
%         (+north, -south; -pi/2 <= dec <= +pi/2)

% output

%  alt = geodetic altitude (kilometers)
%  lat = geodetic latitude (radians)
%        (+north, -south; -pi/2 <= lat <= +pi/2)

% global constants

%  req  = equatorial radius (kilometers)
%  flat = flattening factor (non-dimensional)

% reference: NASA TN D-7522

% Orbital Mechanics with Matlab

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

global req flat

n = req / rmag;
o = flat * flat;

a = 2 * dec;
p = sin(a);
q = cos(a);

a = 4 * dec;
r = sin(a);
s = cos(a);

% geodetic latitude (radians)

lat = dec + flat * n * p + o * n * r * (n - 0.25);

% geodetic altitude (kilometers)

alt = rmag + req * (flat * 0.5 * (1.0 - q) + o * (0.25 * n - 0.0625) * (1.0 - s) - 1.0);

author: Sable

repo: mcbench-benchmarks


text:
%--- help for solve_first_order ---
%
%  INTERNAL FUNCTION
%
%
%    Other functions named solve_first_order
%
%       dsge/solve_first_order
%


text:
function hermite_exactness ( n, x, w, p_max )

%*****************************************************************************80
%
%% HERMITE_EXACTNESS investigates exactness of Hermite quadrature.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points in the rule.
%
%    Input, real X(N), the quadrature points.
%
%    Input, real W(N), the quadrature weights.
%
%    Input, integer P_MAX, the maximum exponent.
%    0 <= P_MAX.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Quadrature rule for Hermite integral\n' );
  fprintf ( 1, '  Rule of order N = %d\n', n );
  fprintf ( 1, '  Degree          Relative Error\n' );
  fprintf ( 1, '\n' );

  for p = 0 : p_max

    s = hermite_integral ( p );

    v(1:n,1) = x(1:n) .^ p;

    q = w' * v;

    if ( s == 0.0 )
      e = abs ( q - s );
    else
      e = abs ( ( q - s ) / s );
    end

    fprintf ( 1, '  %6d  %24.16f\n', p, e );

  end

  return
end


text:
% [func, fx] = mcmc_init_hamiltonian(x_init, get_logpdf, get_dlogpdf, ...
% epsilon, L, func_x)
%
% epsilon is the parameter to the exponential distribution for sampling
% the momentum
%
% L is the number of simulation steps
function [func, fx] = mcmc_init_hamiltonian(x_init, get_logpdf, get_dlogpdf, ...
                                            epsilon, L, func_x)

if nargin < 6
  func_x = @func_x_default;
end

x_current = x_init;
[fx_current, dfx_current] = func_x(x_current);
logpdf_current = get_logpdf(fx_current);
dlogpdf_current = get_dlogpdf(dfx_current);

func = @hamiltonian;
fx = fx_current;

  function [x, fx] = hamiltonian(varargin)

  x = x_current;

  p = normrnd(0,1,size(x_current));
  p_current = p;

  lp_current = logpdf_current(varargin{:});
  dlp_current = dlogpdf_current(varargin{:});

% $$$   mycheckgrad(@chkgrad, x, 1e-6)
% $$$
% $$$     function [y,dy] = chkgrad(x)
% $$$     [fx, dfx] = func_x(x);
% $$$     f = get_logpdf(fx);
% $$$     df = get_dlogpdf(dfx);
% $$$     y = f(varargin{:});
% $$$     dy = df(varargin{:});
% $$$     end

  % Random step size
  e = exprnd(epsilon);

  % Make a half step for momentum at the beginning
  p = p + 0.5*e*dlp_current;

  for l=1:L
    x = x + e * p;

    if any(isnan(x))
      break;
    end

    [fx, dfx] = func_x(x);
    dlogpdf = get_dlogpdf(dfx);
    dlp = dlogpdf(varargin{:});
    if l < L
      p = p + e * dlp;
    else
      logpdf = get_logpdf(fx);
      lp = logpdf(varargin{:});
    end
  end

  if all(~isnan(x))

    % Make a half step for momentum at the end
    p = p + 0.5*e*dlp;

    % Negate momentum to make the proposal symmetric
    p = -p;

    if log(rand()) < ( (lp - 0.5*(p'*p)) - ...
                       (lp_current - 0.5*(p_current'*p_current)) )
      % Accept
      disp('Accept hamiltonian')
      x_current = x;
      fx_current = fx;
      dfx_current = dfx;
      logpdf_current = logpdf;
      dlogpdf_current = dlogpdf;
    else
      disp('Reject hamiltonian')
    end

  else
      disp('Reject hamiltonian')
  end

  x = x_current;
  fx = fx_current;
  end

end

function [fx, dfx] = func_x_default(x)
fx = x;
dfx = x;
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% designsim_gui_script
% Last modified 5/20/01 by Tor Wager		add optional trans2switch from GUI.  still need block.
% also added power calculation.

% sets defaults, reads GUI values from optimizeGUI, and runs simulation

disp('Designsim_gui_script')
disp('===============================')
disp('  ...Setting up variables from GUI.')
% returns vector of t values as t in workspace.

if ~exist('M') | ~isstruct(M),disp('	...No GA results (M struct) detected - using only random designs.'),end

clear SIM
clear H;

% ----------------------------------------------------------------
% * setup defauls
% ----------------------------------------------------------------
% fixed parameters
% --------------------
powerTvoxels = 30000;
powerTdf = 19;
format compact

% parameters from GUI
% --------------------
nonlinthreshold = 4;
noise_var = 0.3;
beta = [.5 .5 0 0;1 1 0 0];	% 20
c = [1 1 -1 -1];					% 13
niterations = 500;
conditions = [1 2 3 4];
ISI = 1.2;
restevery = 9;						% 26
restlength = 2;					% 25
HPlength = 60;						% 23
xc = [1 0];							% 24
trans2switch = 0;

% ----------------------------------------------------------------
% * setup input from GUI
% ----------------------------------------------------------------
% read values from GUI for GA
H = findobj('Tag', 'E1');                                                                                if ~isempty(H)
conditions = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E4');
ISI = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E12');
noise_var = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E13');
c = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E14');
niterations = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E15');
respMean = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E16');
respSD = str2num(get(H(1), 'String'));

H = findobj('Tag', 'C1');
addResponseVariability = get(H(1), 'Value');

% additional values required for random sim.

H = findobj('Tag', 'E2');
freqConditions = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E3');
scanLength = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E5');
TR = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E6');
cbalColinPowerWeights = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E7');
numGenerations = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E8');
sizeGenerations = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E9');
lowerLimit = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E10');
maxOrder = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E11');
plotFlag = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E20');
beta = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E23');
HPlength = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E24');
eval(['load ' get(H(1), 'String')]);

H = findobj('Tag', 'E25');
restlength = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E26');
restevery = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E27');
trans2switch = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E28');
trans2block = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E30');
NumStimthresh = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E31');
maxCbalDevthresh = str2num(get(H(1), 'String')); 	% counterbalancing deviation

H = findobj('Tag', 'E32');
maxFreqDevthresh = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E33');
nonlinthreshold = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E34');
dofirst = str2num(get(H(1), 'String'));

end % ~isempty(H) : means found the GUI

% ----------------------------------------------------------------
% * setup lists, HRF, and other parameters
% ----------------------------------------------------------------

HRF = spm_hrf(.1);
HRF = HRF/ max(HRF);

SIM.beta = beta;
SIM.contrasts = c;
SIM.noise_var = noise_var;
SIM.xc = xc;
SIM.ISI = ISI;
SIM.HPlength = HPlength;
SIM.restevery = restevery;
SIM.restlen = restlength;
SIM.niterations = niterations;
SIM.noise_var = noise_var;
SIM.freqConditions = freqConditions;

SIM.NumStimthresh = NumStimthresh;
SIM.maxCbalDevthresh = maxCbalDevthresh;
SIM.maxFreqDevthresh = maxFreqDevthresh;
SIM.nonlinthreshold = nonlinthreshold;

disp('  ...Generating design list.')
% for building RANDOM designs to test against GA design
fitnessMatrix = zeros(4,sizeGenerations);
numStim = ceil(scanLength / (ISI));
numRestStim = (ceil(numStim/(mean(restevery)+restlength)) - 1) * restlength;
numStimEachCond = ceil((numStim-numRestStim) * freqConditions);			            % row vector of stim in each cond
numsamps = ceil(numStim*ISI/TR);

SIM.numsamps = numsamps;
if exist('M') & isstruct(M), SIM.numsamps = size(M.modelatTR,1);,end
SIM.numStimEachCond = numStimEachCond;

% this is overwritten if trans2switch is used.
unsortedList = [];
for i = 1:size(conditions,2)									% unsorted list of all stims in proper frequencies
   startIndex = size(unsortedList,1)+1;
   unsortedList(startIndex:(startIndex + 2 * numStimEachCond(i)-1),1) = conditions(i);
end

% ----------------------------------------------------------------
% * jitter setup
% ----------------------------------------------------------------

dojitter = 0;
if size(freqConditions,2) == size(conditions,2) + 1
    conditions = [conditions max(conditions)+1];   % insert 0 as last condition if jitter
    dojitter = 1;
elseif size(freqConditions,2) == size(conditions,2)
else error('Frequency and condition vectors must be same length, or freq must be 1 longer for jitter')
end

% ----------------------------------------------------------------
% * criterion measures setup
% ----------------------------------------------------------------
if ~isempty(NumStimthresh) | ~isempty(maxCbalDevthresh) | ~isempty(maxFreqDevthresh)
	docriterion = 1;
	if isempty(NumStimthresh),NumStimthresh = 10000;,end
 	if isempty(maxCbalDevthresh),maxCbalDevthresh = 10000;,end
	if isempty(maxFreqDevthresh),maxFreqDevthresh = 10000;,end
else
	docriterion = 0;
end

maxrestthresh = 2;					% max number of rests in a row.

% ----------------------------------------------------------------
% * stochastic variability in response times - now not used
% ----------------------------------------------------------------
if addResponseVariability == 1

else

% ----------------------------------------------------------------
% * simulation - get smoothing matrix
% ----------------------------------------------------------------
	figure;
    disp(['  ...loading smoothing matrix with hrf LP and HP filter length of ' num2str(HPlength)])
    [hpS] = use_spm_filter(TR,SIM.numsamps,'none','specify',HPlength);
	[fullS] = use_spm_filter(TR,SIM.numsamps,'hrf','specify',HPlength);

% ----------------------------------------------------------------
% * simulation - setup ga simulation with different noise vectors
% ----------------------------------------------------------------
	if exist('M') & isstruct(M),
	    if ISI > TR,stimList = M.stimlist(1:numsamps+1,1);,end   	% cut off to save time if it's too long
                %stimList = insert_rests(stimList,M.restlist,restlength,numStim);not necessary;done in GA
                stimList = sampleInSeconds(M.stimlist,ISI,.1);				% hi-res sample at .1 s, matches HRF
		dmodel = getPredictors(stimList, HRF);
		dmodel = resample(dmodel,1,TR*10);
		dmodel = modelSaturation(dmodel,nonlinthreshold);

		if size(dmodel,1) > SIM.numsamps
			disp('	...truncating dmodel to number of samples')
			dmodel = dmodel(1:SIM.numsamps,:);
		elseif size(dmodel,1) < SIM.numsamps
		        warning('	...stimlist for GA is too short for this ISI!  Results may not be accurate for this ISI!')
	                dmodel = [dmodel; zeros(SIM.numsamps - size(dmodel,1),size(dmodel,2))];
	        end

% ----------------------------------------------------------------
% * simulation - run ga simulation with different noise vectors
% ----------------------------------------------------------------

  		disp(['  ...simulating GA results with noisy data at variance ' num2str(noise_var)])
	        % Test idealized GA - no response variability
		fprintf('		...original	model')
    	[dummy,SIM.ga.se,SIM.ga.t] = designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,0,xc,niterations);
 		fprintf('...high-pass filtered')
    	[dummy,SIM.ga.HPse,SIM.ga.HPt] = designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,hpS,xc,niterations);
		fprintf('...smoothed\n')
		[dummy,SIM.ga.HLse,SIM.ga.HLt] = designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,fullS,xc,niterations);
	end

    disp(['  ...simulating random designs with ' num2str(niterations) ' iterations...']);

SIM

% ----------------------------------------------------------------
% * simulation - setup random simulation - trans2switch
% ----------------------------------------------------------------

% This is part of the switching hack. uses 2 conditions, then figures out 4 columns from that.
% 5 is the rest condition.
if trans2switch
    if size(conditions,2) > 4
	NEWfreqConditions = [.5 .5 freqConditions(5)]; NEWconditions = [1 2 5];
    else NEWfreqConditions = [.5 .5];,NEWconditions = [1 2];
    end
	NEWnumStimEachCond = ceil((numStim-numRestStim) * NEWfreqConditions);			            % row vector of stim in each cond

	disp(['	...USING TRANS2SWITCH > num stim each cond = ' num2str(NEWnumStimEachCond)])

	for i = 1:size(NEWfreqConditions,2)									% unsorted list of all stims in proper frequencies
   		startIndex = size(unsortedList,1)+1;
   		unsortedList(startIndex:(startIndex + 2 * NEWnumStimEachCond(i)-1),1) = NEWconditions(i);
	end
end

% ----------------------------------------------------------------
% * simulation - setup random simulation - trans2block
% ----------------------------------------------------------------
if trans2block
	disp(['	...USING TRANS2BLOCK >'])
        freqConditions = [freqConditions./2 freqConditions./2];
        disp(['     freqConditions = ' num2str(freqConditions)])

        if dojitter
           maxcond = max(conditions)-1;
           conditions = [conditions(1:maxcond) maxcond+1:2*length(conditions)-1];
           disp(['     CONDITION ' num2str(conditions(end)) ' is blank jitter.'])
	else
           maxcond = max(conditions);
           conditions = [conditions maxcond+1:2*length(conditions)];
        end
        disp(['     conditions = ' num2str(conditions)])

	% ----------------------------------------------------------------
	% * setup special first trial
	% ----------------------------------------------------------------
	if isempty(dofirst)
		dofirst = 0;
	else
		disp(['	...block design: adding special predictor for first trial in block.'])
		try
			% this is really kludgy, but should give a rough estimate of how many first trials to expect.
			numfirst = numStim ./ mean(restevery);
			proportionfirst = numfirst./numStim;
		catch
			error('restevery not found...must use rest periods with transform2block!')
		end
		freqConditions = freqConditions - proportionfirst./size(freqConditions,2);
		freqConditions = [freqConditions proportionfirst];
		disp('	...first trial special condition - adjusted expected frequencies:')
		freqConditions
	end
end

% ----------------------------------------------------------------
% * simulation - setup random simulation - list and rest matrix
% ----------------------------------------------------------------
		disp('	...Setting up lists.')
                clear restMatrix, clear listMatrix
		% make rest matrix
		atleast = min(restevery);
		restPeriods = numStim / atleast + 5;	% make this many rows in restMatrix
		unsortedrestList = [];
		for i = 1:size(restevery,2):restPeriods
			unsortedrestList = [unsortedrestList;restevery'];
		end

		for i = 1:niterations
			restMatrix(:,i) = getRandom(unsortedrestList);
			listMatrix(:,i) = getRandom(unsortedList);
		end

		listMatrix = insert_rests(listMatrix,restMatrix,restlength,numStim);

		% insert rests before transforming
		for i = 1:niterations
			if trans2switch,
				listMatrix(:,i) = transform2switches(listMatrix(:,i));
		    end
			if trans2block,
             			listMatrix(:,i) = transform2block(listMatrix(:,i),restMatrix(:,i),restlength,dofirst);
      			end
                end

% ----------------------------------------------------------------
% * simulation - setup random simulation - criterion measures
% ----------------------------------------------------------------

              if docriterion
                disp('   ...Making sure all lists meet criteria...')

                % criterion measures
	        % ====================================================================
	        for i = 1:niterations
                	stimList = listMatrix(:,i);
			restList = restMatrix(:,i);

                	go = 0;tryindex = 1;

                	while ~go
                  		[maxNumStim,maxrest] = getMaxInARow(stimList);

                  		% ====== Frequency Deviation =======================================
                  		for j = 1:size(conditions,2)

         				freqMat(j) = sum(stimList == conditions(j)) / size(stimList,1);
      		  		end
                  		freqMat = freqMat ./ sum(freqMat);   % adjust for probe ''0''s
		  		maxFreqDev = max(abs(freqConditions(1:size(conditions,2)) - freqMat));

                  		% ====== Counterbalancing Deviation =======================================
                  		[cBal,dummy,maxDev] = getCounterBal(stimList, maxOrder,conditions,freqConditions);

                  		% ====== Test list for criteria satisfaction ==============================
		  		if maxNumStim < NumStimthresh & maxrest < maxrestthresh & maxDev < maxCbalDevthresh & maxFreqDev < maxFreqDevthresh
		     		go = 1;
		  		else
                     		restList        = getRandom(unsortedrestList);
			 		stimList		 = getRandom(unsortedList);
			 		stimList        = insert_rests(stimList,restList,restlength,numStim);
                     		if trans2switch,
                        		stimList = transform2switches(stimList);
                     		end
		     		if trans2block,
             					stimList = transform2block(stimList,restList,restlength,dofirst);
      			 		end

                     		tryindex = tryindex+1;
                  		end

		     	if mod(tryindex,100) == 0,disp(['Tried 100 times on model ' num2str(i)])
                          disp(['Avg. max stim in a row = ' num2str(maxNumStim)])
                          disp(['Avg. deviation from freqs = ' num2str(maxFreqDev)])
                          disp(['Avg. deviation from counterbal = ' num2str(maxDev)])
                          disp(['Avg. frequencies are = ' num2str(freqMat)])
                          disp(['Max restcond is ' num2str(maxrest)])
					conditions
					freqConditions
					freqMat
                     	end
				if mod(tryindex,500) == 0
					disp('Cannot meet criteria after 500 tries.')
					conditions
					freqConditions
					freqMat
					cBal
					error('Exiting now.')
				end
                	end  % end while

                	restMatrix(:,i) = restList;
			listMatrix(:,i) = stimList;
                	if mod(i,50) == 0, fprintf('%d ',i),end
                end  % end for ... iterations
	     end  % if docriterion

% ----------------------------------------------------------------
% * simulation - setup random simulation - make models
% ----------------------------------------------------------------

	for i = 1:niterations
      	% Test Random models

	  	% Prepare Stimulus List
      	        stimList = listMatrix(:,i);
		SIM.stimList = stimList;
		stimList = sampleInSeconds(stimList,ISI,.1);
		dmodel = getPredictors(stimList, HRF);
		dmodel = resample(dmodel,1,TR*10);
		dmodel = modelSaturation(dmodel,nonlinthreshold);

		if size(dmodel,1) > SIM.numsamps
			dmodel = dmodel(1:SIM.numsamps,:);
		elseif size(dmodel,1) < SIM.numsamps
			dmodel = [dmodel; zeros(SIM.numsamps - size(dmodel,1),size(dmodel,2))];
		end

% ----------------------------------------------------------------
% * simulation - run random simulation - no smoothing
% ----------------------------------------------------------------
% this gets the average over [numnoise] simulations with different noise vectors
   		if i == 1
	  		[dummy,SIM.rnd.se(:,i),SIM.rnd.t(:,i),GLM] = ...
			designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,0,xc);	%added 4th output to plot
			set(gcf,'Position',[88         514        1185         420])	;drawnow
			subplot(4,1,4); plot(dmodel(:,1));title('First predictor of model')
			drawnow
		else
			[dummy,SIM.rnd.se(:,i),SIM.rnd.t(:,i)] = ...
			designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,0,xc);
		end

% ----------------------------------------------------------------
% * simulation - run random simulation - HP filter
% ----------------------------------------------------------------
	% high-pass only
	[dummy, SIM.rnd.HPse(:,i),SIM.rnd.HPt(:,i)] = ...
	designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,hpS,xc);

% ----------------------------------------------------------------
% * simulation - run random simulation - high and low filter
% ----------------------------------------------------------------

      	% now with smoothing
      	[dummy, SIM.rnd.HLse(:,i),SIM.rnd.HLt(:,i)] = designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,fullS,xc);
      	fprintf('.')
	  if mod(i,30) == 0, fprintf(' %d ',i),end
	  if mod(i,90) == 0, fprintf('\n'),end

    end
	% for niterations

end
		% simulation with fixed response times.

% ----------------------------------------------------------------
% * power calculation for all sims
% ----------------------------------------------------------------

powerTcrit = tinv(1-.05/powerTvoxels,powerTdf);
for i = {'SIM.ga.t' 'SIM.ga.HPt' 'SIM.ga.HLt' 'SIM.rnd.t' 'SIM.rnd.HPt' 'SIM.rnd.HLt'}
	eval([i{1} 'power = sum(abs(' i{1} ') > powerTcrit) / length(' i{1} ');'])
end

	try
		plotSim(SIM)
	catch
		disp('Error plotting. command: plotSim(SIM)')
	end

save latestSIM SIM

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore