text:
function [ fea, out ] = ex_diffusion1( varargin )
%EX_DIFFUSION1 2D Diffusion equation example on a unit square.
%
%   [ FEA, OUT ] = EX_DIFFUSION1( VARARGIN ) Diffusion equation on a unit
%   square with exact solutions. Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       isol        scalar {1}             Exact solution
%                           1                x*y
%                           2                x^2-y^2
%                           3                2*y/((1+x)^2+y^2)
%                           4                2*y/((1+x)^2+y^2)
%                           5                (sinh(pi*x)*sin(pi*y)+sinh(pi*y)*
%                                            sin(pi*x))/sinh(pi)
%       cd          scalar {0.1}           Diffusion coefficient
%       igrid       scalar 1/{0}           Cell type (0=quadrilaterals, 1=triangles)
%       hmax        scalar {1/40}          Max grid cell size
%       sfun        string {sflag1}        Shape function
%       iphys       scalar 0/{1}           Use physics mode to define problem    (=1)
%                                          or directly define fea.eqn/bdr fields (=0)
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { ...
  'isol',     1; ...
  'cd',       0.1; ...
  'igrid',    0; ...
  'hmax',     1/40; ...
  'sfun',     'sflag1'; ...
  'iphys',    1; ...
  'iplot',    1; ...
  'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;
switch opt.isol
  case 1
    refsol = 'x*y';
  case 2
    refsol = 'x^2-y^2';
  case 3
    refsol = '2*y/((1+x)^2+y^2)';
  case 4
    refsol = 'exp(pi*x)*sin(pi*y)';   % pi can be substituted for any constant.
  case 5
    refsol = '(sinh(pi*x)*sin(pi*y)+sinh(pi*y)*sin(pi*x))/sinh(pi)';
end

% Geometry definition.
gobj = gobj_rectangle();
fea.geom.objects = { gobj };

% Grid generation.
switch opt.igrid
  case -1
    fea.grid = rectgrid(round(1/opt.hmax));
    fea.grid = quad2tri(fea.grid);
  case 0
    fea.grid = rectgrid(round(1/opt.hmax));
  case 1
    fea.grid = gridgen(fea,'hmax',opt.hmax,'fid',fid);
end
n_bdr = max(fea.grid.b(3,:));           % Number of boundaries.

% Problem definition.
fea.sdim  = { 'x' 'y' };                % Coordinate names.
if ( opt.iphys==1 )

  fea = addphys(fea,@convectiondiffusion);   % Add convection and diffusion physics mode.
  fea.phys.cd.sfun = { opt.sfun };           % Set shape function.
  fea.phys.cd.eqn.coef{2,4}   = { opt.cd };  % Set diffusion coefficient.
  fea.phys.cd.bdr.sel         = [1 1 1 1];
  fea.phys.cd.bdr.coef{1,end} = repmat({refsol},1,n_bdr);   % Set Dirichlet boundary coefficient to reference solution.
  fea = parsephys(fea);                 % Check and parse physics modes.

else

  fea.dvar  = { 'c' };                  % Dependent variable name.
  fea.sfun  = { opt.sfun  };            % Shape function.

  % Define equation system.
  fea.eqn.a.form = { [2 3;2 3] };       % First row indicates test function space   (2=x-derivative + 3=y-derivative),
                                        % second row indicates trial function space (2=x-derivative + 3=y-derivative).
  fea.eqn.a.coef = { opt.cd };          % Diffusion coefficient used in assembling stiffness matrix.

  fea.eqn.f.form = { 1 };               % Test function space to evaluate in right hand side (1=function values).
  fea.eqn.f.coef = { 0 };               % Coefficient used in right hand side.

  % Define boundary conditions.
  fea.bdr.d     = cell(1,n_bdr);
 [fea.bdr.d{:}] = deal(refsol);         % Assign reference solution to all boundaries (Dirichlet).

  fea.bdr.n     = cell(1,n_bdr);        % No Neumann boundaries ('fea.bdr.n' empty).

end

% Parse and solve problem.
fea       = parseprob(fea);             % Check and parse problem struct.
fea.sol.u = solvestat(fea,'fid',fid);   % Call to stationary solver.

% Postprocessing.
s_err = ['abs(',refsol,'-c)'];
if ( opt.iplot>0 )
  figure
  subplot(2,1,1)
  postplot(fea,'surfexpr','c')
  title('Solution c')
  subplot(2,1,2)
  postplot(fea,'surfexpr',s_err,'evalstyle','exact')
  title('Error')
end

% Error checking.
if ( size(fea.grid.c,1)==4 )
  xi = [0;0];
else
  xi = [1/3;1/3;1/3];
end
err = evalexpr0(s_err,xi,1,1:size(fea.grid.c,2),[],fea);
ref = evalexpr0('c',xi,1,1:size(fea.grid.c,2),[],fea);
err = sqrt(sum(err.^2)/sum(ref.^2));

if( ~isempty(fid) )
  fprintf(fid,'\nL2 Error: %f\n',err)
  fprintf(fid,'\n\n')
end

out.err  = err;
out.pass = out.err<0.01;
if ( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function matRad_compileOmpMCInterface(dest,omcFolder)
% Compiles the ompMC interface (integrated as submodule)
%
% call
%   matRad_compileOmpMCInterface()
%   matRad_compileOmpMCInterface(dest)
%   matRad_compileOmpMCInterface(dest,sourceFolder)
%
% input:
%   dest:           (optional) destination for mex file. Default: location
%                   of this file
%   sourceFolder:   (optional) path to ompMC . Default assumes its checked
%                   out in the submodules folder of matRad
%
% References
%
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2020 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

matRad_cfg = MatRad_Config.instance();

env = matRad_getEnvironment();

if nargin < 1
    dest = fileparts(mfilename('fullpath'));
end

if nargin < 2
    omcFolder = [matRad_cfg.matRadRoot filesep 'submodules' filesep 'ompMC'];
end

sourceFolder = [omcFolder filesep 'src'];
interfaceFolder = [omcFolder filesep 'ucodes' filesep 'omc_matrad'];

mainFile = [interfaceFolder filesep 'omc_matrad.c'];

addFiles = {'ompmc.c','omc_utilities.c','omc_random.c'};
addFiles = cellfun(@(f) fullfile(sourceFolder,f),addFiles,'UniformOutput',false);

addFiles = strjoin(addFiles,' ');

if exist ('OCTAVE_VERSION','builtin')
    ccName = evalc('mkoctfile -p CC');
else
    myCCompiler = mex.getCompilerConfigurations('C','Selected');
    ccName = myCCompiler.ShortName;
end

%These settings have only been tested for MSVC and g++. You may need to adapt for other compilers
if ~isempty(strfind(ccName,'MSVC')) %Not use contains(...) because of octave
    flags{1,1} = 'COMPFLAGS';
    flags{1,2} = '/openmp';
    flags{2,1} = 'OPTIMFLAGS';
    flags{2,2} = '/O2';
else
    flags{1,1} = 'CFLAGS';
    flags{1,2} = '-std=gnu99 -fopenmp -O3';
    flags{2,1} = 'LDFLAGS';
    flags{2,2} = '-fopenmp';

end

includestring =  ['-I' sourceFolder];

flagstring = '';

%For Octave, the flags will be set in the environment, while they
%will be parsed as string arguments in MATLAB
for flag = 1:size(flags,1)
    if strcmp(env,'OCTAVE')
        preFlagContent = evalc(['mkoctfile -p ' flags{flag,1}]);
        if ~isempty(preFlagContent)
            preFlagContent = preFlagContent(1:end-1); %Strip newline
        end
        newContent = [preFlagContent ' ' flags{flag,2}];
        setenv(flags{flag,1},newContent);
        matRad_cfg.dispDebug('Set compiler flag %s to %s\n',flags{flag,1},newContent);
    else
        flagstring = [flagstring flags{flag,1} '="$' flags{flag,1} ' ' flags{flag,2} '" '];
    end
end

mexCall = ['mex -largeArrayDims ' flagstring ' ' includestring ' ' mainFile ' ' addFiles];
matRad_cfg.dispDebug('Compiler call: %s\n',mexCall);

currDir = pwd;
cd(dest);
eval(mexCall);
cd(currDir);
end

repo: matRad


text:
% Program to form Admittance And Impedance Bus Formation....
% with Transformer Tap setting..

function ybus = ybusppg()  % Returns ybus

linedata = linedata30(); % Calling "linedata3.m" for Line Data...
fb = linedata(:,1);     % From bus number...
tb = linedata(:,2);     % To bus number...
r = linedata(:,3);      % Resistance, R...
x = linedata(:,4);      % Reactance, X...
b = linedata(:,5);    % Ground Admittance, B/2...
a = linedata(:,6);      % Tap setting value..
z = r + i*x;            % Z matrix...
y = 1./z;               % To get inverse of each element...
b = i*b;                % Make B imaginary...

nbus = max(max(fb),max(tb));    % no. of buses...
nbranch = length(fb);           % no. of branches...
ybus = zeros(nbus,nbus);        % Initialise YBus...

 % Formation of the Off Diagonal Elements...
 for k = 1:nbranch
     ybus(fb(k),tb(k)) = ybus(fb(k),tb(k))-y(k)/a(k);
     ybus(tb(k),fb(k)) = ybus(fb(k),tb(k));
 end

 % Formation of Diagonal Elements....
 for m = 1:nbus
     for n = 1:nbranch
         if fb(n) == m
             ybus(m,m) = ybus(m,m) + y(n)/(a(n)^2) + b(n);
         elseif tb(n) == m
             ybus(m,m) = ybus(m,m) + y(n) + b(n);
         end
     end
 end

author: Sable

repo: mcbench-benchmarks


text:
function odf = rotate(odf,rot,varargin)
% rotate ODF
%
% Input
%  odf - @SO3Fun
%  q   - @rotation
%
% Output
%  rotated odf - @SO3Fun

ss = odf.SS.Laue;
if length(ss)>2 && ~any(rot == ss(:))
  warning('Rotating an ODF with specimen symmetry will remove the specimen symmetry')
  odf.SS = specimenSymmetry;
end

for i = 1:length(odf.components)
  odf.components{i} = odf.components{i}.rotate(rot,varargin{:});
end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function figure = plotlygetfile(file_owner, file_id)

    [un, key, domain] = signin;

    headers = struct(...
                    'name',...
                        {...
                            'plotly-username',...
                            'plotly-apikey',...
                            'plotly-version',...
                            'plotly-platform',...
                            'user-agent'
                        },...
                    'value',...
                        {...
                            un,...
                            key,...
                            plotly_version,...
                            'MATLAB',...
                            'MATLAB'
                        });

    url = [domain, '/apigetfile/', file_owner, '/', num2str(file_id)];

    [response_string, extras] = urlread2(url, 'Post', '', headers);
    response_handler(response_string, extras);
    response_object = loadjson(response_string);
    figure = response_object.payload.figure;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function area = area_quad ( quad_xy )

%*****************************************************************************80
%
%% AREA_QUAD returns the area of a quadrilateral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real QUAD_XY(2,4), the coordinates of the nodes.
%
%    Output, real AREA, the area of the quadrilateral.
%
  t1(1:2,1) = quad_xy(1:2,1);
  t1(1:2,2) = quad_xy(1:2,2);
  t1(1:2,3) = quad_xy(1:2,3);

  area1 = triangle_area ( t1 );

  t2(1:2,1) = quad_xy(1:2,3);
  t2(1:2,2) = quad_xy(1:2,4);
  t2(1:2,3) = quad_xy(1:2,1);

  area2 = triangle_area ( t2 );

  if ( area1 < 0.0 | area2 < 0.0 )

    t1(1:2,1) = quad_xy(1:2,2);
    t1(1:2,2) = quad_xy(1:2,3);
    t1(1:2,3) = quad_xy(1:2,4);

    area1 = triangle_area ( t1 );

    t2(1:2,1) = quad_xy(1:2,4);
    t2(1:2,2) = quad_xy(1:2,1);
    t2(1:2,3) = quad_xy(1:2,2);

    area2 = triangle_area ( t2 );

    if ( area1 < 0.0 | area2 < 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'AREA_QUAD - Fatal error!\n' );
      fprintf ( 1, '  The quadrilateral nodes seem to be listed in\n' );
      fprintf ( 1, '  the wrong order, or the quadrilateral is\n' );
      fprintf ( 1, '  degenerate.\n' );
      error ( 'AREA_QUAD - Fatal error!' );
    end

  end

  area = area1 + area2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef MOSD < ALGORITHM
% <multi> <real> <large/none> <constrained/none>
% Multiobjective steepest descent
% step --- 0.1 --- Step size

%------------------------------- Reference --------------------------------
% X. Liu and A. C. Reynolds, A multiobjective steepest descent method with
% applications to optimal well control, Computational Geosciences, 2016,
% 20: 355-374.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB Platform
% for Evolutionary Multi-Objective Optimization [Educational Forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            step = Algorithm.ParameterSet(0.1);

            %% Generate random population
            Population = Problem.Initialization(Problem.N);
            Archive    = UpdateArchive(Population,Problem.N);
            step0      = step;

            %% Optimization
            while Algorithm.NotTerminated(Archive)
                for i = 1 : Problem.N
                    gk  = FiniteDifference(Problem,Population(i));
                    gk1 = norm(gk(:,1));
                    gk2 = norm(gk(:,2));
                    if gk1 <= gk2
                        gs = gk(:,1);
                        gl = gk(:,2);
                    else
                        gs = gk(:,2);
                        gl = gk(:,1);
                    end
                    if gs'*(gs-gl) <= 0
                        d = -gs/norm(gs);
                    else
                        D = (((norm(gl).^2-gs'*gl)/(norm(gs).^2-gs'*gl))*(-gs))-gl;
                        d = D/norm(D);
                    end
                    Offspringdec = Population(i).dec+step*d';
                    Offspring    = Problem.Evaluation(Offspringdec);
                    Archive      = UpdateArchive([Archive,Offspring],Problem.N);
                    if ~any(Offspring.obj<Population(i).obj)
                        Population(i) = Archive(randi(end));
                        step = step/2;
                    else
                        Population(i) = Offspring;
                        step = min([2*step,step0]);
                    end
                end
            end
        end
    end
end

function df = FiniteDifference(Problem,X)
    if any(X.con>0)
        df = Problem.CalConGrad(X.dec)';
    else
        df = Problem.CalObjGrad(X.dec)';
    end
end

function P = UpdateArchive(P,N)
    P = P(NDSort(P.objs,P.cons,1)==1);
    if length(P) > N
        Choose = true(1,length(P));
        Dis    = pdist2(P.objs,P.objs);
        Dis(logical(eye(length(Dis)))) = inf;
        while sum(Choose) > N
            Remain   = find(Choose);
            Temp     = sort(Dis(Remain,Remain),2);
            [~,Rank] = sortrows(Temp);
            Choose(Remain(Rank(1))) = false;
        end
        P = P(Choose);
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% Fig. 9.33   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clear all;
close all;
%saturation nonlinearity (Figure 9.33)

N=0.1;
k= 1;
n=10
j=1;
for i=0.1:0.01:n
    Keq(j) = (2/pi)*(k*asin(N/(k*i))+(N/i)*sqrt(1-(N/(k*i))^2));
    j=j+1;
end;
plot([0,0.1:0.01:n],[1 Keq]);
axis([0 10 0 1.1])
title('Describing function for saturation nonlinearity')
xlabel('a')
ylabel('K_{eq}');
nicegrid;
hold off


text:
function determ = line_adj_determinant ( n )

%*****************************************************************************80
%
%% LINE_ADJ_DETERM returns the determinant of the LINE_ADJ matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM, the determinant.
%
       if ( mod ( n, 4 ) == 1 )
    determ =   0.0;
  elseif ( mod ( n, 4 ) == 2 )
    determ = - 1.0;
  elseif ( mod ( n, 4 ) == 3 )
    determ =   0.0;
  elseif ( mod ( n, 4 ) == 0 )
    determ = + 1.0;
  end

  return
end


text:
function hermite_cubic_test10 ( )

%*****************************************************************************80
%
%% HERMITE_CUBIC_TEST10 tests HERMITE_CUBIC_SPLINE_INTEGRAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 February 2011
%
%  Author:
%
%    John Burkardt
%
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_CUBIC_TEST10:\n' );
  fprintf ( 1, '  HERMITE_CUBIC_SPLINE_INTEGRAL integrates a Hermite\n' );
  fprintf ( 1, '  cubic spline over the definition interval [X1,XNN].\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  If the subintervals are equally spaced, the derivative\n' );
  fprintf ( 1, '  information has no effect on the result, except for\n' );
  fprintf ( 1, '  the first and last values, DN(1) and DN(NN).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '                            Exact       Computed\n' );
  fprintf ( 1, '     X1          XNN        Integral    Integral  Comment\n' );
  fprintf ( 1, '\n' );

  for test = 1 : 5
%
%  Equal spacing.
%
    if ( test == 1 )
      nn = 11;
      xn = linspace ( 0.0, pi, nn );
      fn = sin ( xn );
      dn = cos ( xn );
      integral_exact = - cos ( xn(nn) ) + cos ( xn(1) );
      comment = 'Equal spacing, correct DN';
%
%  Equal spacing, reset DN(2:NN-1) to random numbers.
%
    elseif ( test == 2 )
      nn = 11;
      xn = linspace ( 0.0, pi, nn );
      fn = sin ( xn );
      dn = cos ( xn );
      [ dn(2:nn-1), seed ] = r8vec_uniform_01 ( nn - 2, seed );
      dn(2:nn-1) = 1000.0 * dn(2:nn-1);
      integral_exact = - cos ( xn(nn) ) + cos ( xn(1) );
      comment = 'Equal spacing, DN(2:N-1) random';
%
%  Equal spacing, now reset all of DN to random numbers.
%
    elseif ( test == 3 )

      nn = 11;
      xn = linspace ( 0.0, pi, nn );
      fn = sin ( xn );
      [ dn, seed ] = r8vec_uniform_01 ( nn, seed );
      dn(1:nn) = 1000.0 * dn(1:nn);
      integral_exact = - cos ( xn(nn) ) + cos ( xn(1) );
      comment = 'Equal spacing, DN(1:N) random';
%
%  Variable spacing, correct data.
%
    elseif ( test == 4 )
      nn = 11;
      xn = linspace ( 0.0, pi^2, nn );
      xn = sqrt ( xn );
      fn = sin ( xn );
      dn = cos ( xn );
      integral_exact = - cos ( xn(nn) ) + cos ( xn(1) );
      comment = 'Variable spacing, correct DN';
%
%  Variable spacing, change one entry in DN.
%
    elseif ( test == 5 )
      nn = 11;
      xn = linspace ( 0.0, pi^2, nn );
      xn = sqrt ( xn );
      fn = sin ( xn );
      dn = cos ( xn );
      [ r, seed ] = r8_uniform_01 ( seed );
      dn( floor ( nn + 1 ) / 2 ) = 1000.0 * r;
      integral_exact = - cos ( xn(nn) ) + cos ( xn(1) );
      comment = 'Variable spacing, a single internal DN randomized.';
    end

    integral_computed = hermite_cubic_spline_integral ( nn, xn, fn, dn );

    fprintf ( 1, '  %10f  %10f  %10.6g  %10.6g  %s\n', ...
      xn(1), xn(nn), integral_exact, integral_computed, comment );

  end

  return
end


text:
% std_figtitle() - Generate plotting figure titles in a cell array
%
% Usage:
%   >> celltxt = std_figtitle(key1, val1, key2, val2);
%
% Inputs:
%  'subject'     - [string] Subject name
%  'datatype'    - [string] data type (For example 'erp')
%  'chanlabels'  - [cell array or string] channel names
%  'compnames'   - [cell array or string] component names
%  'condnames'   - [cell array] names of conditions
%  'cond2names'  - [cell array] names of conditions
%  'vals'        - [cell array or real array] value for plot panel (for
%                  example, 0 ms)
%  'valsunit'    - [cell array or string] unit value (i.e. ms)
%
% Optional statistical titles:
%  'condstat'    - ['on'|'off'] default is 'off'
%  'cond2stat'   - ['on'|'off'] default is 'off'
%  'mcorrect'    - [string] correction for multiple comparisons. Default is
%                  empty.
%  'statistics'  - [string] type of statictics
%  'threshold'   - [real] treshold value
%
% Optional grouping:
%  'condgroup'   - ['on'|'off'] group conditions together default is 'off'
%  'cond2group'  - ['on'|'off'] default is 'off'
%
% Outputs:
%   celltxt      - cell array of text string, one for each set of condition
%
% Example:
%  std_figtitle('subject', 'toto', 'condnames', { 'test1' 'test2' }, ...
%               'vals' , { [100] [4] }, 'valsunit', { 'ms' 'Hz' })
%
% Authors: Arnaud Delorme, SCCN/UCSD, Feb 2010

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

function [all_titles alllegends ] = std_figtitle(varargin)

if nargin < 1
    help std_figtitle;
    return;
end;

alllegends = {};
opt = finputcheck( varargin, { 'chanlabels'  {'cell' 'string'}   []     {};
                               'condnames'   {'cell' 'string'}   []     '';
                               'cond2names'  {'cell' 'string'}   []     '';
                               'condstat'    'string'            {'on' 'off'}     'off';
                               'cond2stat'   'string'            {'on' 'off'}     'off';
                               'condgroup'   'string'            {'on' 'off' 'together' 'apart'}     'off';
                               'cond2group'  'string'            {'on' 'off' 'together' 'apart'}     'off';
                               'plotmode'    'string'            {'normal' 'condensed'}              'normal';
                               'plotsubjects' 'string'           {'on' 'off'}                        'off';
                               'threshold'   'real'              []     NaN;
                               'statistics'  'string'            []     '';
                               'mcorrect'    'string'            []     '';
                               'datatype'    'string'            []     '';
                               'clustname'   'string'            []     '';
                               'compnames'   {'cell' 'string'}   []     {};
                               'vals'        {'cell' 'real'}     []     {}; % just for titles
                               'valsunit'    {'cell' 'string'}   []     {}; % just for titles
                               'subject'     'string'            []              '' }, 'std_figtitle'); %, 'ignore');
if isstr(opt), error(opt); end;
if ~iscell(opt.vals),       opt.vals       = { opt.vals }; end;
ncori  = length(opt.condnames);
nc2ori = length(opt.cond2names);
if ~isempty(opt.vals) && ~isempty(opt.vals{1}) && ~isnan(opt.vals{1}(1)), opt.condgroup = 'off'; opt.cond2group = 'off'; end;
if strcmpi(opt.plotmode, 'condensed'), opt.condgroup = 'on';  opt.cond2group = 'on';  end;
if strcmpi(opt.plotsubjects, 'on'),    opt.condgroup = 'off'; opt.cond2group = 'off'; end;
if strcmpi(opt.plotsubjects, 'on')
    opt.condgroup  = 'apart';
    opt.cond2group = 'apart';
end;
if strcmpi(opt.condgroup,  'on'), opt.condgroup =  'together'; end;
if strcmpi(opt.cond2group, 'on'), opt.cond2group = 'together'; end;

if strcmpi(opt.condgroup,  'together') &&  strcmpi(opt.cond2stat, 'on'),  opt.condgroup  = 'apart'; end;
if strcmpi(opt.cond2group, 'together') &&  strcmpi(opt.condstat , 'on') , opt.cond2group = 'apart'; end;
if ~( strcmpi(opt.condgroup,  'together') && strcmpi(opt.cond2group, 'together') )
    if strcmpi(opt.condgroup,  'together'), alllegends = opt.condnames ; opt.condnames  = ''; end;
    if strcmpi(opt.cond2group, 'together'), alllegends = opt.cond2names; opt.cond2names = ''; end;
end;
if ~iscell(opt.valsunit),   opt.valsunit   = { opt.valsunit }; end;
if ~iscell(opt.chanlabels), opt.chanlabels = { opt.chanlabels }; end;
if ~iscell(opt.condnames),  opt.condnames  = { opt.condnames }; end;
if ~iscell(opt.cond2names), opt.cond2names = { opt.cond2names }; end;
if isempty(opt.condnames),  opt.condnames{1}  = ''; end;
if isempty(opt.cond2names), opt.cond2names{1} = ''; end;

for c1 = 1:length(opt.condnames)
    for c2 = 1:length(opt.cond2names)

        % value (ms or Hz)
        % ----------------
        fig_title1 = '';
        for i = 1:length(opt.vals)
            if ~isempty(opt.vals{i}) && ~isnan(opt.vals{i}(1))
                if opt.vals{i}(1) == opt.vals{i}(end), fig_title1 = [ num2str(opt.vals{i}(1)) '' opt.valsunit{i} fig_title1];
                else                                   fig_title1 = [ num2str(opt.vals{i}(1)) '-' num2str(opt.vals{i}(2)) '' opt.valsunit{i} fig_title1 ];
                end;
                if length(opt.vals) > i, fig_title1 = [ ' & ' fig_title1 ]; end;
            end;
        end;

        % conditions
        % ----------
        if ~isempty(opt.condnames{c1})
            fig_title1 = [ value2str(opt.condnames{c1}) ', ' fig_title1];
        end;
        if ~isempty(opt.cond2names{c2})
            fig_title1 = [ value2str(opt.cond2names{c2}) ', ' fig_title1];
        end;

        % channel labels, component name, subject name and datatype
        % ---------------------------------------------------------
        fig_title2 = '';
        if length( opt.chanlabels ) == 1 && ~isempty( opt.chanlabels{1} )
            fig_title2 = [ opt.chanlabels{1} ', ' fig_title2 ];
        end;

        % cluster and component name
        % --------------------------
        if ~isempty( opt.clustname )
            if ~isempty( opt.compnames )
                if iscell( opt.compnames )
                     compstr = [ 'C' num2str(opt.compnames{min(c1,size(opt.compnames,1)),min(c2,size(opt.compnames,2))}) ];
                else compstr = [ opt.compnames ];
                end;
            else compstr = '';
            end;
            if ~isempty( opt.subject )
                if ~isempty( compstr )
                     fig_title2 = [ opt.subject '/' compstr ', ' fig_title2 ];
                else fig_title2 = [ opt.subject ', ' fig_title2 ];
                end;
            elseif ~isempty( compstr )
                fig_title2 = [ compstr ', ' fig_title2 ];
            end;

            if ~isempty( opt.datatype )
                 fig_title2 = [ opt.clustname ' ' opt.datatype ', ' fig_title2 ];
            else fig_title2 = [ opt.clustname ', ' fig_title2 ];
            end;

        else
            if ~isempty( opt.compnames )
                if iscell( opt.compnames )
                else fig_title2 = [ 'C' num2str(opt.compnames{c1,c2}) ', ' fig_title2 ];
                     fig_title2 = [ opt.compnames        ', ' fig_title2 ];
                end;
            end;

            % subject and data type
            % ---------------------
            if ~isempty( opt.subject )
                if ~isempty( opt.datatype )
                     fig_title2 = [ opt.subject ' ' opt.datatype ', ' fig_title2 ];
                else fig_title2 = [ opt.subject ', ' fig_title2 ];
                end;
            elseif ~isempty( opt.datatype )
                fig_title2 = [ opt.datatype ' - ' fig_title2 ];
            end;
        end;

        if strcmpi(opt.cond2group, 'together') && strcmpi(opt.condgroup, 'together')
            fig_title = fig_title2;
            if ~isempty(fig_title1) && strcmpi(fig_title1(end-1:end), ', '), fig_title1(end-1:end) = []; end;
            if ~isempty(fig_title1) && strcmpi(fig_title1(end-1:end), '- '), fig_title1(end-1:end) = []; end;
            alllegends{c1, c2} = fig_title1;
        else
            fig_title = [ fig_title2 fig_title1 ];
        end;
        if ~isempty(fig_title) && strcmpi(fig_title(end-1:end), ', '), fig_title(end-1:end) = []; end;
        if ~isempty(fig_title) && strcmpi(fig_title(end-1:end), '- '), fig_title(end-1:end) = []; end;
        all_titles{c1,c2}  = fig_title;

    end;
end;
if ~isempty(alllegends)
    alllegends = alllegends';
    alllegends = alllegends(:)';

    % convert legends to string if necessary
    % --------------------------------------
    for ileg = 1:length(alllegends), alllegends{ileg} = num2str(alllegends{ileg}); end;
end;

% statistic titles
% ----------------
if isnan(opt.threshold),
    basicstat = '(p-value)';
else
    if length(opt.threshold) >= 1
        basicstat = sprintf([ '(p<%.' thresh_pres(opt.threshold(1)) 'f)'], opt.threshold(1));
    end;
    if length(opt.threshold) >= 2
        basicstat = [ basicstat(1:end-1) sprintf([' red;p<%.' thresh_pres(opt.threshold(2))  'f brown)'], opt.threshold(2)) ];
    end;
    if length(opt.threshold) >= 3
        basicstat = [ basicstat(1:end-1) sprintf([';\np<%.' thresh_pres(opt.threshold(3)) 'f black)' ], opt.threshold(3)) ];
    end;
end;
if ~isempty(opt.statistics), basicstat = [ basicstat ' ' opt.statistics    ]; end;
if ~isempty(opt.mcorrect) && ~strcmpi(opt.mcorrect, 'none'),   basicstat = [ basicstat ' with ' opt.mcorrect ]; end;
if strcmpi(opt.condstat, 'on')
    rown = size(all_titles,1)+1;
    for c2 = 1:length(opt.cond2names)
        all_titles{rown, c2} = [ num2str(opt.cond2names{c2}) ' ' basicstat ];
    end;
end;
if strcmpi(opt.cond2stat, 'on')
    coln = size(all_titles,2)+1;
    for c1 = 1:length(opt.condnames)
        all_titles{c1, coln} = [ num2str(opt.condnames{c1}) ' ' basicstat ];
    end;
end;
if strcmpi(opt.condstat, 'on') && strcmpi(opt.cond2stat, 'on')
   all_titles{rown, coln} = [ 'Interaction ' basicstat ];
end;

function pres = thresh_pres(thresh_pres);
    if (round(thresh_pres*100)-thresh_pres*100) == 0
        pres = 2;
    elseif (round(thresh_pres*1000)-thresh_pres*1000) == 0
        pres = 3;
    else
        pres = -round(log10(thresh_pres))+1;
    end;
    pres = num2str(pres);

% convert
% -------
function str = value2str(value)

    if isstr(value)
        str = value;
    elseif isnumeric(value)
        if length(value) == 1
            str = num2str(value);
        else
            str = num2str(value(1));
            if length(value) <= 5
                for ind = 2:length(value)
                    str = [ str ' & ' num2str(value(ind)) ];
                end;
            else
                str = [ str ' & ' num2str(value(2)) ' & ...' ];
            end;
        end;
    else % cell array
        str = value{1};
        for ind = 2:length(value)
            str = [ str ' & ' value{ind} ];
        end;
    end;

author: PatternRecognition

repo: OpenBMI


text:
function square_grid_test01 ( )

%*****************************************************************************80
%
%% SQUARE_GRID_TEST01 tests SQUARE_GRID using the same parameters for all dimensions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 August 2014
%
%  Author:
%
%    John Burkardt
%
  a = [ -1.0, -1.0 ];
  b = [ +1.0, +1.0 ];
  c = [ 1, 1 ];
  ns = [ 3, 3 ];

  n = ns(1) * ns(2);

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SQUARE_GRID_TEST01\n' );
  fprintf ( 1, '  Create a grid using SQUARE_GRID.\n' );
  fprintf ( 1, '  Use the same parameters in every dimension.\n' );
  fprintf ( 1, '  Number of grid points N = %d\n', n );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I    NS     C      A         B\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 2
    fprintf ( 1, '  %4d  %4d  %4d  %8.4f  %8.4f\n', i, ns(i), c(i), a(i), b(i) );
  end

  x = square_grid ( n, ns, a, b, c );
  r8mat_transpose_print ( 2, n, x, '  Grid points:' );

  return
end


text:
function test_model = test_nrsfm(test,model_3d,test_em_iters)
% Use a trained shape model to predict keypoints on new instances
% given the bounding box and visible keypoints
% NOTE: The keypoints are already scaled to make all bboxes the same size
    params = get_params();
    P = test.points;
    MD = isnan(P(1:end/2,:));
    P(isnan(P)) = 0;
    tol = params.nrsfm.tol;
    max_em_iter = test_em_iters;
    try
        seg.poly_x = test.poly_x;
        seg.poly_y = test.poly_y;
    catch
        seg = test.mask;
    end
    % Do EM for inference with missing data
    [P3, RO, Tr, basis_coeffs,c] = em_sfm_known_shape(P, MD, model_3d.S_bar,...
        model_3d.defBasis, model_3d.sigma_sq, seg, test.imsize, test.bbox, tol, max_em_iter);

    %Set up test_model

    test_model.class = model_3d.class;
    test_model.S_bar = model_3d.S_bar;
    test_model.defBasis = model_3d.defBasis;
    test_model.part_names = model_3d.part_names;
    test_model.sigma_sq = model_3d.sigma_sq;
    test_model.points3 = P3;
    test_model.rots = RO;
    test_model.trs = Tr;
    test_model.def_weights = basis_coeffs;
    test_model.voc_image_id = test.voc_image_id;
    if(isfield(test,'voc_rec_id'))
        test_model.voc_rec_id = test.voc_rec_id;
    end
    if(isfield(test,'flip'))
        test_model.flip = test.flip;
    end
    test_model.bbox = test.bbox;
    test_model.points2 = test.points;
    test_model.params = params;
    test_model.c = c;
    test_model.seg = seg;
end

repo: CategoryShapes


text:
function mpc = case6ww
%CASE6WW  Power flow data for 6 bus, 3 gen case from Wood & Wollenberg.
%   Please see CASEFORMAT for details on the case file format.
%
%   This is the 6 bus example from pp. 104, 112, 119, 123-124, 549 of
%   "Power Generation, Operation, and Control, 2nd Edition",
%   by Allen. J. Wood and Bruce F. Wollenberg, John Wiley & Sons, NY, Jan 1996.

%   MATPOWER

%% MATPOWER Case Format : Version 2
mpc.version = '2';

%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 100;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [
	1	3	0	0	0	0	1	1.05	0	230	1	1.05	1.05;
	2	2	0	0	0	0	1	1.05	0	230	1	1.05	1.05;
	3	2	0	0	0	0	1	1.07	0	230	1	1.07	1.07;
	4	1	70	70	0	0	1	1	0	230	1	1.05	0.95;
	5	1	70	70	0	0	1	1	0	230	1	1.05	0.95;
	6	1	70	70	0	0	1	1	0	230	1	1.05	0.95;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
mpc.gen = [
	1	0	0	100	-100	1.05	100	1	200	50	0	0	0	0	0	0	0	0	0	0	0;
	2	50	0	100	-100	1.05	100	1	150	37.5	0	0	0	0	0	0	0	0	0	0	0;
	3	60	0	100	-100	1.07	100	1	180	45	0	0	0	0	0	0	0	0	0	0	0;
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [
	1	2	0.1	0.2	0.04	40	40	40	0	0	1	-360	360;
	1	4	0.05	0.2	0.04	60	60	60	0	0	1	-360	360;
	1	5	0.08	0.3	0.06	40	40	40	0	0	1	-360	360;
	2	3	0.05	0.25	0.06	40	40	40	0	0	1	-360	360;
	2	4	0.05	0.1	0.02	60	60	60	0	0	1	-360	360;
	2	5	0.1	0.3	0.04	30	30	30	0	0	1	-360	360;
	2	6	0.07	0.2	0.05	90	90	90	0	0	1	-360	360;
	3	5	0.12	0.26	0.05	70	70	70	0	0	1	-360	360;
	3	6	0.02	0.1	0.02	80	80	80	0	0	1	-360	360;
	4	5	0.2	0.4	0.08	20	20	20	0	0	1	-360	360;
	5	6	0.1	0.3	0.06	40	40	40	0	0	1	-360	360;
];

%%-----  OPF Data  -----%%
%% generator cost data
%	1	startup	shutdown	n	x1	y1	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
	2	0	0	3	0.00533	11.669	213.1;
	2	0	0	3	0.00889	10.333	200;
	2	0	0	3	0.00741	10.833	240;
];

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [row_sum,col_sum]=calc_sum(Image,rowsize,z)
col_sum=sum(Image)'/z;
for i=1:rowsize,
row_sum(i)=sum(Image(i,:))/z;
end


text:
function [block] = spm_vb_lambda(Y,block)
% Variational Bayes for GLM-AR models - Update lambda
% FORMAT [block] = spm_vb_lambda(Y,block)
%
% Y      - [T x N] time series
% block  - data structure (see spm_vb_glmar)
%__________________________________________________________________________
% Copyright (C) 2005-2014 Wellcome Trust Centre for Neuroimaging

% Will Penny and Nelson Trujillo-Barreto
% $Id: spm_vb_lambda.m 6079 2014-06-30 18:25:37Z spm $

if block.verbose
    disp('Updating lambda');
end

p = block.p;
k = block.k;
N = block.N;

for n=1:N
    if p > 0
        % Equation 77 in paper VB1
        Gn          = spm_vb_get_Gn (Y,block,n);
    else
        subblock_n  = [(n-1)*k+1:n*k];
        en          = Y(:,n) - block.X*block.w_mean(subblock_n,1);
        Gn          = trace(block.w_cov{n}*block.XTX) + en'*en;
    end
    % Equation 75 in paper VB1
    block.b_lambda(n,1)    = 1./(Gn./2 + 1./block.b_lambda_prior(n));
    block.mean_lambda(n,1) = block.c_lambda(n)*block.b_lambda(n);
end

author: spm


text:
%
%   colour transfer algorithm based on linear Monge-Kantorovitch solution
%
%   IR = colour_transfer_MKL(I_original, I_target, nbiterations);
%
%  (c) F. Pitie 2007
%
%  see reference:
%
%
function IR = colour_transfer_MKL(I0, I1)

if (ndims(I0)~=3)
    error('pictures must have 3 dimensions');
end

X0 = reshape(I0, [], size(I0,3));
X1 = reshape(I1, [], size(I1,3));

A = cov(X0);
B = cov(X1);

T = MKL(A, B);

mX0 = repmat(mean(X0), [size(X0,1) 1]);
mX1 = repmat(mean(X1), [size(X0,1) 1]);

XR = (X0-mX0)*T + mX1;

IR = reshape(XR, size(I0));

function [T] = MKL(A, B)
N = size(A,1);
[Ua,Da2] = eig(A);
Da2 = diag(Da2);
Da2(Da2<0) = 0;
Da = diag(sqrt(Da2 + eps));
C = Da*Ua'*B*Ua*Da;
[Uc,Dc2] = eig(C);
Dc2 = diag(Dc2);
Dc2(Dc2<0) = 0;
Dc = diag(sqrt(Dc2 + eps));
Da_inv = diag(1./(diag(Da)));
T = Ua*Da_inv*Uc*Dc*Uc'*Da_inv*Ua';

author: xjqicuhk

repo: SIMS

save_path: github-repos/MATLAB/xjqicuhk-SIMS


text:
% repeat a matrix
%
function [grad] = B_transpose(future_layers, curr_layer)
future_grad = GetFutureGrad(future_layers, curr_layer);

grad = future_grad';

end

author: singaxiong

repo: SignalGraph


text:
%% AMG TEST III: Robustness to time discretization
%
% We consider the linear finite element discretization of heat equation on
% the unstructured mesh with Neumann boundary conditions. We test the
% implicit time discretization with various time stepsizes dt = 1/h^4,
% 1/h^2, 1/h, 1.

%%
load lakemesh
%% Time step: h^4. Mass matrix will dominate.
close all;
dt = 'h^4';
[N,itStep,time,err] = amgtest(node,elem,2,[],dt);
%%
colHeaders = {'Unknowns','Iterations','Time (sec)','Error'};
makeHtmlTable([N itStep time err],[],[],colHeaders,[],6);
%%
r = showrate(N,time,3);
xlabel('N'); ylabel('Time');
title(['Complexity is N^{' num2str(r,3) '}'] ,'Fontsize', 14);

%% Time step: h^2
close all;
dt = 'h^2';
[N,itStep,time,err] = amgtest(node,elem,2,[],dt);
%%
colHeaders = {'Unknowns','Iterations','Time (sec)','Error'};
makeHtmlTable([N itStep time err],[],[],colHeaders,[],6);
%%
r = showrate(N,time,3);
xlabel('N'); ylabel('Time');
title(['Complexity is N^{' num2str(r,3) '}'] ,'Fontsize', 14);

%% Time step:  h
close all;
dt = 'h';
[N,itStep,time,err] = amgtest(node,elem,2,[],dt);
%%
colHeaders = {'Unknowns','Iterations','Time (sec)','Error'};
makeHtmlTable([N itStep time err],[],[],colHeaders,[],6);
%%
r = showrate(N,time,3);
xlabel('N'); ylabel('Time');
title(['Complexity is N^{' num2str(r,3) '}'] ,'Fontsize', 14);

%% Time step:  1. Like regularized Neumann problem.
close all;
dt = '1';
[N,itStep,time,err] = amgtest(node,elem,2,[],dt);
%%
colHeaders = {'Unknowns','Iterations','Time (sec)','Error'};
makeHtmlTable([N itStep time err],[],[],colHeaders,[],6);
%%
r = showrate(N,time,3);
xlabel('N'); ylabel('Time');
title(['Complexity is N^{' num2str(r,3) '}'] ,'Fontsize', 14);

repo: Aero-matlab


text:
function [ns_RESULT, nsLibraryInfo] = ns_GetLibraryInfo();

%ns_GetLibraryInfo   Get library version information
%
%   Usage:
%       [ns_RESULT, nsLibraryInfo] = ns_GetLibraryInfo()
%
%   Description:
%       Obtains information about the API library.
%
%   Return Values:
%       nsLibraryInfo   Structure to receive library version information.
%       ns_RESULT   This function returns ns_OK if the file is successfully
%                   opened. Otherwise one of the following error codes is
%                   generated:
%
%                       ns_LIBERROR	Library Error
%
%   Copyright (C) 2003 Neuroshare Project
%   Author: Almut Branner
%   Last modification: 8/11/2003

[ns_RESULT, nsLibraryInfo] = mexprog(2);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [poly,poly_idx,opendat,boudat] = extract_boundary(v_start,v_end,bnde,pts,order,opendat,boudat,type,type2)
% DESCRIPTION: Given a set of boundary edges and a starting and ending index
%              of a singly- or multi-polygonal region, organize them in a
%              winding order and/or add them to an existing opendat/boudat
%              structure. The program breaks up the boundary so that it
%              will never exceed 100,000 nodes to avoid memory allocation
%              issues adcprep.
%
% INPUTS:
%      v_start: the starting index of the boundary you want to trace
%        v_end: the ending index of the boundary you want to trace.
%         bnde: the indices of each boundary edge as a nbnde x 2 matrix
%         pts:  the x,y locations of all the points in the region
%               stored as an np x 2 matrix.
%         order:the order in which the traversal takes place
%               counter-clockwise (0) or clockwise (1).
%         opendat: open boundary information from a pre-existing grid
%          boudat: land boundary information from a pre-exist
% OUTPUTS:
%          poly: the boundary of each enclosing polygon sorted in winding-order
%                poly is returned as a cell-array of length number of polys.
%      poly_idx: indices of the polygon coordinates in the same format as
%                poly
%       opendat: open boundary information with appened open bou
%       boudat: land boundary information with appended land bou
%         type: flux(1) or elevation (2) type bc
%        type2: if flux, zero flux (20) or non-zero flux river (22) bc
%
% kjr,UND,CHL,2017
%
%                                           TRAVERSAL METHOD
% Pick any unvisited edge segment [v_start,v_next] and add these vertices to the polygon loop.
% Find the unvisited edge segment [v_i,v_j] that has either v_i = v_next or v_j = v_next and add the other vertex (the one not equal to v_next) to the polygon loop.
% Reset v_next as this newly added vertex, mark the edge as visited and continue from 2.
% Traversal is done when we get back to v_start.
% NOTE: that the signed area will be positive if the vertices are
% oriented counterclockwise, and will be negative if it is oriented clockwise
bnde= unique(bnde,'rows');
active = true(size(bnde,1),1);
p = 0;
exceed = 50e3; %100e3;

[rt,~] = find(v_start==bnde);
if isempty(rt), disp('v_start does not exist on boundary, check numbering'); return; end
r  = rt(order+1); % change this only here to from 1 to 2 or 2 to 1 to go left or right
tsel = bnde(r,:);
sel  = tsel(tsel~=v_start);
v_next = sel;
active(r) = 0;
% cut up boundary so land boundary never exceeds 100k nodes.
cut = true;
while cut
    p = p + 1;
    if(p > 1 )
        [rt,~] = find(v_next==bnde & active);
        r  = rt(1);
        tsel = bnde(r,:);
        sel  = tsel(tsel~=v_next);
        v_next = sel;
        active(r) = 0;
    end

    temp = [];
    temp2= [];

    temp  = pts(bnde(r,:)',:);
    temp2 = bnde(r,:)';
    if v_next ~= temp2(2)
      % swap
      temp = flipud(temp);
      temp2 = flipud(temp2);
    end

    k = 2;
    while v_next~=v_end % terminates when we reach v_end
        rt= (v_next==bnde(:,1) | v_next==bnde(:,2)) &  active;
        r = find(rt,1);
        tsel = bnde(r,:);
        sel=tsel(tsel~=v_next);
        k = k + 1;
        temp(k,:)= pts(sel,:);
        temp2(k,:)= sel;
        active(r) = 0;
        v_next = sel;
        % exceeded max xize, break
        if (k > exceed),disp('exceed'); break, end
        % reached ending vertex, break
        if(v_next==v_end), cut=false; disp('reached ending'); end
        % exhausted all edges and couldn't connect
        if(~any(active)), cut=false; disp('coudln''t conntect'); break, end
    end
    if length(temp) == 2
        cut = false;
    end
    poly{p}     = temp;
    poly_idx{p} = temp2;
    [area] = parea(poly{p}(:,1),poly{p}(:,2));
    if order == 0 % ccw
        if sign(area)<0
            poly{p} = flipud(poly{p});
            poly_idx{p} = flipud(poly_idx{p});
        end
    else % cw
        if sign(area)>0
            poly{p} = flipud(poly{p});
            poly_idx{p} = flipud(poly_idx{p});
        end
    end
end
for ii = 1 : p
    hold on; plot(poly{ii}(:,1),poly{ii}(:,2),'r-','linewi',2);
end

if ~exist('type','var')
   type = input('What kind of boundary is this, 1 (flux) or 2 (elevation)?');
end

% if populated
if ~isempty(boudat)
    if(type==1)
        nbou = boudat.nbou;
        nvel = boudat.nvel;
        nvell= boudat.nvell;
        nbvv = boudat.nbvv;
        ibtype = boudat.ibtype;
        if ~exist('type2','var')
           type2 = input('What kind of flux boundary is it, 20 (Mainland), 21 (Island), 22 (River)?');
        end
        for ii = 1 : length(poly)
            nbou = nbou + 1;
            nvell(nbou) = length(poly{ii}(:,1));
            nvel = nvel + nvell(nbou);
            nbvv(1:nvell(nbou),nbou) = int32(poly_idx{ii}(:));
            ibtype(nbou) = type2;
        end
        boudat.nbou = nbou ;
        boudat.nvel = nvel ;
        boudat.nvell = nvell ;
        boudat.ibtype = ibtype ;
        boudat.nbvv = nbvv ;
    end
end

% if populated
if ~isempty(opendat)
    if type==2
        nope = opendat.nope;
        nvdll= opendat.nvdll;
        neta = opendat.neta;
        ibtype=opendat.ibtype;
        nbdv  = opendat.nbdv;

        for ii = 1 : length(poly)
            nope = nope + 1;
            nvdll(nope) = length(poly_idx{ii}(:,1));
            neta = neta + nvdll(nope);
            ibtype(nope) = 0;
            nbdv(1:nvdll(nope),nope) = poly_idx{ii}(:,1);
        end
        % ocean boundary
        opendat.nope = nope ;
        opendat.neta = neta ;
        opendat.nvdll = nvdll ;
        opendat.ibtype = ibtype ;
        opendat.nbdv = nbdv ;

    end
end

% if empty
if isempty(opendat)
    if(type==2)
        % nothing populated
        nope = 0 ;
        nvdll = []  ;
        neta = 0 ;
        nbdv = [] ;
        for ii = 1 : length(poly)
            nope = nope + 1;
            nvdll(nope) = length(poly_idx{ii}(:,1));
            neta = neta + nvdll(nope);
            ibtype(nope) = 0;
            nbdv(1:nvdll(nope),nope) = poly_idx{ii}(:,1);
        end
        % ocean boundary
        opendat.nope = nope ;
        opendat.neta = neta ;
        opendat.nvdll = nvdll ;
        opendat.ibtype = ibtype ;
        opendat.nbdv = nbdv ;
    end
end

% if empty
if isempty(boudat)
    if type==1
        nbou = 0; nvel = 0; nvell = []; nbbv = [] ;
        if ~exist('type2','var')
           type2 = input('What kind of flux boundary is it, 20(island),22(River)?');
        end
        for ii = 1 : length(poly)
            nbou = nbou + 1;
            nvell(nbou) = length(poly{ii}(:,1));
            nvel = nvel + nvell(nbou);
            nbvv(1:nvell(nbou),nbou) = poly_idx{ii}(:);
            ibtype(nbou) = type2;
        end
        boudat.nbou = nbou ;
        boudat.nvel = nvel ;
        boudat.nvell = nvell ;
        boudat.ibtype = ibtype ;
        boudat.nbvv = nbvv ;
    end
end

end
% helper function, computes area of polygon
function [area]=parea(x,y)
n    = length(x);
xp   = [x; x(1)];
yp   = [y; y(1)];
area = 0;
for i = 1:n
    area = area + det([xp(i), xp(i+1); yp(i), yp(i+1)]);
end
area = 1/2*area;
end

author: CHLNDDEV


text: %ISNAN Dataset overload

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function model = gpReadFromFID(FID, varargin)

% GPREADFROMFID Load from a FID produced by the C++ implementation.
% FORMAT
% DESC loads in from a file stream the data format produced by the
% C++ GP implementation.
% ARG FID : the file ID from where the data is loaded.
% RETURN model : the model loaded in from the file.
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006, 2008, 2009
%
% SEEALSO : gpReadFromFile

% GP

try
  version = readVersionFromFID(FID);
catch
  if strcmp(lasterr, 'Incorrect file format')
    version = readDoubleFromFID(FID, 'gplvmVersion');
  else
    error(lasterr)
  end
end

if version > 0.2
  error('Incorrect file version.')
end

if version>0.11
  baseType = readStringFromFID(FID, 'baseType')
    if ~strcmp(baseType, 'dataModel') && ~strcmp(baseType, 'mapmodel')
    error('Incorrect base type in file.')
  end
  type = readStringFromFID(FID, 'type')
  if ~strcmp(type, 'gp')
    error('Incorrect type in file.')
  end
end
numData = readIntFromFID(FID, 'numData');
dataDim = readIntFromFID(FID, 'outputDim');
inputDim = readIntFromFID(FID, 'inputDim');
sparseApprox = readIntFromFID(FID, 'sparseApproximation');
numActive = readIntFromFID(FID, 'numActive');
if sparseApprox
  beta = modelReadFromFID(FID);
  beta = beta(1, 1);
end

learnScale = readBoolFromFID(FID, 'learnScale');
learnBias = readBoolFromFID(FID, 'learnBias');
scale = modelReadFromFID(FID);
bias = modelReadFromFID(FID);

kern = modelReadFromFID(FID);
noise = modelReadFromFID(FID);

if sparseApprox
  fixInducing = readBoolFromFID(FID, 'fixInducing');
  X_u = modelReadFromFID(FID);
end

X = varargin{1};
y = varargin{2};
% X = zeros(numData, inputDim);
% y = zeros(numData, dataDim);

warning('Noise model is ignored');

switch sparseApprox
 case 0
  approxType = 'ftc';
 case 1
  approxType = 'dtc';
 case 2
  approxType = 'fitc';
 case 3
  approxType = 'pitc';
 case 4
  approxType = 'dtcvar';
end
options = gpOptions(approxType);
options.numActive = numActive;
options.kern = kern;
if sparseApprox
  options.fixInducing = fixInducing;
end
model = gpCreate(inputDim, size(y, 2), X, y, options);
model.X = X;
if sparseApprox
  model.X_u = X_u;
  model.beta = beta;
end
model.scale = scale;
model.bias = bias;
model.m = gpComputeM(model);

% This forces kernel computation.
initParams = gpExtractParam(model);
model = gpExpandParam(model, initParams);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function res = spm_eeg_specest_morlet(S, data, time)
% Plugin for spm_eeg_tf implementing Morlet wavelet transform
% FORMAT res = spm_eeg_specest_morlet(S, data, time)
%
% S                     - input structure
% fields of S:
%    S.subsample   - factor by which to subsample the time axis (default - 1)
%  either
%    S.ncycles     - Morlet wavelet factor (default - 7)
%  or
%    S.timeres     - Fixed time window length in ms
%
%    S.frequencies - vector of frequencies (default - 0-48) at optimal frequency bins
%
% Output:
%  res -
%   If no input is provided the plugin returns a cfg branch for itself
%
%   If input is provided:
%      res.fourier - the complex output of wavelet transform
%      res.time    - time axis
%      res.freq    - frequency axis
%__________________________________________________________________________
% Copyright (C) 2010-2017 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_specest_morlet.m 7129 2017-07-04 16:24:53Z guillaume $

%-This part if for creating a config branch that plugs into spm_cfg_eeg_tf
% Any parameters can be specified and they are then passed to the plugin
% when it's called.
%--------------------------------------------------------------------------
if nargin == 0
    subsample         = cfg_entry;
    subsample.tag     = 'subsample';
    subsample.name    = 'Subsample';
    subsample.strtype = 'n';
    subsample.num     = [1 1];
    subsample.val     = {1};
    subsample.help    = {'Set to N to subsample the time axis to every Nth sample (to reduce the dataset size).'};

    ncycles         = cfg_entry;
    ncycles.tag     = 'ncycles';
    ncycles.name    = 'Number of wavelet cycles';
    ncycles.strtype = 'n';
    ncycles.num     = [1 1];
    ncycles.val     = {7};
    ncycles.help    = {'Number of wavelet cycles (a.k.a. Morlet wavelet factor).',...
        'This parameter controls the time-frequency trade-off',...
        'Increasing it increases the frequency resolution at the expense of time resolution.'};

    timeres         = cfg_entry;
    timeres.tag     = 'timeres';
    timeres.name    = 'Fixed time window length';
    timeres.strtype = 'r';
    timeres.num     = [1 1];
    timeres.val     = {0};
    timeres.help    = {'Fixed time window for all frequencies.',...
        'Specify time window length in ms.',...
        'Default valued of 0 specifies variable time window length.'};

    morlet      = cfg_branch;
    morlet.tag  = 'morlet';
    morlet.name = 'Morlet wavelet transform';
    morlet.val  = {ncycles, timeres, subsample};

    res = morlet;

    return;
end

%-Defaults
%--------------------------------------------------------------------------
try, S.subsample; catch, S.subsample = 1; end
try, S.ncycles;   catch, S.ncycles = 7;   end

dt = time(end) - time(1);
if ~isfield(S, 'frequencies') || isempty(S.frequencies)
    S.frequencies = (1/dt):max(1/dt, floor(dt)/dt):48;
end

%-Generate wavelets
%--------------------------------------------------------------------------
if ~isfield(S, 'timeres') || S.timeres == 0
    M = spm_eeg_morlet(S.ncycles, 1000*diff(time(1:2)), S.frequencies);
else
    M = spm_eeg_morlet(S.ncycles, 1000*diff(time(1:2)), S.frequencies, 1000./S.timeres);
end

%-Data dimensions
%--------------------------------------------------------------------------
Nchannels    = size(data, 1);
Nsamples     = size(data, 2);
Nfrequencies = numel(M);

%-Initialize output struct
%--------------------------------------------------------------------------
res.freq    = S.frequencies;
res.time    = time(1:S.subsample:end);
res.fourier = zeros(Nchannels, Nfrequencies, length(res.time));

%-Compute wavelet transform
%--------------------------------------------------------------------------
for j = 1:Nchannels
    for i = 1:Nfrequencies
        tmp = conv(data(j, :), M{i});

        % time shift to remove delay
        tmp = tmp([1:Nsamples] + (length(M{i})-1)/2);

        tmp = tmp(1:S.subsample:end);

        res.fourier(j, i, :) = tmp;
    end
end

%% The following is often faster but requires more memory
%for i = 1:Nfrequencies
%    H = spm_convmtx(M{i}',Nsamples); % faster than conv for large Nchannels
%    tmp = data * H';
%    % subsample and time shift to remove delay
%    tmp = tmp(:,(1:S.subsample:Nsamples) + (length(M{i})-1)/2);
%    res.fourier(:,i,:) = tmp;
%end


text:
function [rc,rcb,Pf,Pb] = pc2rcv(pc,R0)

%function [rc,rcb] = pc2rcv(pc,R0)
% Transforms normalized correlation matrices pc into
% forward and backward reflection coefficients rc and rcb.

% S. de Waele, March 2003.

if ~isstatv(pc), error('Partial correlations non-stationairy!'), end

s = kingsize(pc);
order = s(3)-1;
dim = s(1); I = eye(dim);
%if nargin == 1, R0 = I; end

rc = zeros(s);  rc(:,:,1)  = I;
Pf   = zeros(s); Pf(:,:,1)   = R0;

rcb = zeros(s);  rcb(:,:,1) = I;
Pb   = zeros(s); Pb(:,:,1)  = R0;

for p = 1:order,
   TsqrtPf = Tsqrt( Pf(:,:,p)); %square root M defined by: M=Tsqrt(M)*Tsqrt(M)'
   TsqrtPb= Tsqrt(Pb(:,:,p));
   %reflection coefficients
   rc(:,:,p+1) = -TsqrtPf *pc(:,:,p+1) *inv(TsqrtPb);
   rcb(:,:,p+1)= -TsqrtPb*pc(:,:,p+1)'*inv(TsqrtPf );
   %residual matrices
   Pf(:,:,p+1)  = (I-TsqrtPf *pc(:,:,p+1) *pc(:,:,p+1)'*inv(TsqrtPf ))*Pf(:,:,p);
   Pb(:,:,p+1) = (I-TsqrtPb*pc(:,:,p+1)'*pc(:,:,p+1) *inv(TsqrtPb))*Pb(:,:,p);
end %for p = 2:order,


text:
if 1
  % This is the setup used in the paper
  a=40;
  M=60;
  L=a*M;
  W=4;
  nrep=50;
else
  a=50;
  M=200;
  L=a*M;
  W=4;
  nrep=50;

end;

system('rm crossover.log');
%for gl=M:M:20*M
for gl=M:M:16*M
  s=sprintf('./time_dgt_fb %i %i %i %i %i %i >> crossover.log\n',a,M,L,W,gl,nrep);

  disp(s);
  system(s);
end;

s=sprintf('./time_dgt_fac %i %i %i %i %i > crossover.ref\n',a,M,L,W,nrep);

disp(s);
system(s);

system('rm crossover_real.log');
%for gl=M:M:20*M
for gl=M:M:16*M
  s=sprintf('./time_dgtreal_fb %i %i %i %i %i %i >> crossover_real.log\n',a,M,L,W,gl,nrep);

  disp(s);
  system(s);
end;

s=sprintf('./time_dgtreal_fac %i %i %i %i %i > crossover_real.ref\n',a,M,L,W,nrep);

disp(s);
system(s);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% prtBrvMc - Bayesian Random Variable Monte Carlo

classdef prtBrvMcmc

    methods (Abstract)
        [self, training, samples] = mcmc(self, x);
    end

    properties
        mcmcTotalIterations = 2e3; % Total number of iterations
        mcmcDiscardIterations = 1e3; % The first mcmcDiscardIterations will be removed
        mcmcRetainEveryOtherItation = 5; % Draw 1:mcmcRetainEveryOtherItation:end will be retained
        mcmcVerboseText = false; % display text
        mcmcVerbosePlot = false; % plot after each iteration
        mcmcVerboseMovie = false; % make a movie with each frame of plotting
        mcmcVerboseMovieFrames = []; % Where we store the frames
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [vec,rotMat]=ITRS2TEME(x,Jul1,Jul2,deltaTTUT1,xpyp,LOD)
%%ITRS2TEME Convert from the International Terrestrial Reference
%          System (ITRS) into the True Equator Mean Equinox (TEME) of date
%          coordinate system. The TEME system is non-standard and is
%          generally only used in the Specialized General Perturbations 4
%          (SGP4) orbital propagation algorithm. Note that the velocity
%          conversion does not include the (small) centrifugal effect of
%          polar motion.
%
%INPUTS: x The NXnumVec collection of vectors to convert. N can be 3, or
%          6. If the vectors are 3D, then they are position. 6D vectors
%          are assumed to be position and velocity, whereby the angular
%          velocity of the Earth's rotation is taken into account using a
%          non-relativistic formula.
% Jul1, Jul2 Two parts of a Julian date given in terrestrial time (TT).
%          The units of the date are days. The full date is the sum of
%          both terms. The date is broken into two parts to provide more
%          bits of precision. It does not matter how the date is split.
%  deltaTTUT1 An optional parameter specifying the difference between TT
%          and UT1 in seconds. This information can be obtained from
% http://www.iers.org/nn_11474/IERS/EN/DataProducts/EarthOrientationData/eop.html?__nnn=true
%          or
% http://www.usno.navy.mil/USNO/earth-orientation/eo-products
%          If this parameter is omitted or if an empty matrix is passed,
%          then the value provided by the function getEOP will be used
%          instead.
%     xpyp xpyp=[xp;yp] are the polar motion coordinates in radians
%          including the effects of tides and librations. If this
%          parameter is omitted or an empty matrix is passed the value
%          from the function getEOP will be used.
%      LOD The difference between the length of the day using terrestrial
%          time, international atomic time, or UTC without leap seconds
%          and the length of the day in UT1. This is an instantaneous
%          parameter (in seconds) proportional to the rotation rate of the
%          Earth. This is only needed if more than just position
%          components are being converted.
%
%The conversion from the TEME to the pseudo-Earth-Fixed (PEF) coordinate
%system is described in [1] and the relationship between the ITRS and the
%PEF is described in [2].
%
%The velocity transformation deals with the instantaneous rotational
%velocity of the Earth using a simple Newtonian velocity addition.
%Basically, the axis of rotation in the Pseudo-Ears-Fixed (PEF) frame is
%the z-axis (The PEF is akin to a less-accurate version of the TIRS).
%The rotation rate in that system is Constants.IERSMeanEarthRotationRate
%adjusted using the Length-of-Day (LOD) Earth Orientation Parameter (EOP).
%Thus, in the PEF, the angular velocity vector is [0;0;omega], where omega
%is the angular velocity accounting for the LOD EOP. Consequently, one
%account for rotation by transforming from the ITRS to the PEF,
%adding the cross product of Omega with the position in the PEF, and
%then converting to the TEME. This is a simple Newtonian conversion.
%
%The algorithm can be compiled for use in Matlab  using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%[vec,rotMat]=ITRS2TEME(x,Jul1,Jul2);
%or if more parameters are known,
%[vec,rotMat]=ITRS2TEME(x,Jul1,Jul2,deltaTTUT1,xpyp,LOD);
%
%REFERENCES:
%[1] D. A. Vallado, P. Crawford, R. Hujsak, and T. Kelso, "Implementing
%    the revised SGP4 in STK," in Proceedings of the AGI User Exchange,
%    Washington, DC, 17?18 Oct. 2006, slides. [Online].
%    Available: http://www.agi.com/downloads/events/2006-agi-user-exchange/8_revised_sgp4_vallado2.pdf
%[2] D. A. Vallado, J. H. Seago, and P. K. Seidelmann, "Implementation
%    issues surrounding the new IAU reference systems for astrodynamics,"
%    in Proceedings of the 16th AAS/AIAA Space Flight Mechanics
%    Conference, Tampa, FL, 22?26 Jan. 2006. [Online].
%    Available: http://www.centerforspace.com/downloads/files/pubs/AAS-06-134.pdf
%
%December 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% clc; clear all; close all;

%% hydro data
hydro = struct();
hydro = readAQWA(hydro, 'cylinder.AH1', 'cylinder.LIS');
hydro = radiationIRF(hydro,5,[],[],[],[]);
hydro = radiationIRFSS(hydro,[],[]);
hydro = excitationIRF(hydro,5,[],[],[],[]);
writeBEMIOH5(hydro)

%% Plot hydro data
% plotBEMIO(hydro)

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function [R, Z] = infermfa(X, LX, MX, PX)
%INFERMFA Infer MoFA using information from EM algorithm in MPPCA
%
%   [R, Z] = infermfa(X, LX, MX, PX)
%
% Computes local data representations and responsibilities of the datapoints
% to the clusters specified by LX, MX, and PX. Basically, this method performs
% an additional E-step of the EM-algorithm executed in MPPCA.
%
%

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

    % Initialize some variables
    [D N] = size(X);
    [D no_dims no_analyzers] = size(LX);
    const = -D / 2 * log(2 * pi);
    R  = zeros(no_analyzers, N);
    Z  = zeros(no_dims, N, no_analyzers);

    % Estimate the cluster centers based on input (= E-step)
    pii = 1 ./ PX;
    for kk=1:no_analyzers
        l        = LX(:,:,kk);
        ltpi     = (repmat(pii, [1 no_dims]) .* l)';
        ltpil    = ltpi * l;
        iltpil   = eye(no_dims) + ltpil;
        cc       = chol(iltpil);
        cci      = inv(cc);
        covz     = cci * cci';
        delta    = X - MX(:,kk * ones(1, N));
        meanz    = ((eye(no_dims) - ltpil * covz) * ltpi) * delta;

        Z(:,:,kk)  = meanz;
        R(kk,:)    = -.5 * (pii' * (delta .* delta) - sum(meanz .* (iltpil * meanz), 1)) - ...
                      sum(log(diag(cc)));
    end

    % Compute responsibilities of clusters to points
    R = R + const + .5 * sum(log(pii));
    R = exp(R - repmat(max(R, [], 1), [no_analyzers 1]));
    R = R ./ repmat(sum(R, 1), [no_analyzers 1]);

repo: cluster


text:
function [A,B,C]=spherTriAngles4Sides(a,b,c,r)
%%SPHERTRIANGLES4SIDES Given the lengths of the sides of a spherical
%       triangle, find the angles of the triangle. Note that the angles of
%       a spherical triangle do not sum to pi. This is just an application
%       of the spherical law of cosines. Note that unlike with a EUclidean
%       triangle, one cannot just enter arbitrary distances for a, b, and c
%       (even if much smaller than the circumference of the sphere),
%       because certain distance relations are not possible and will result
%       in non-real solutions.
%
%INPUTS: a, b, c The lengths of the sides of the spherical triangle
%                in any order. If one wishes to convert multiple triangles
%                at once, these can be matrices.
%              r The scalar radius of the sphere on which the distances are
%                measured. If omitted or an empty matrix is passed, the
%                default of 1 is used.
%
%OUTPUTS: A, B, C The angles of the spherical triangle in radians. A is the
%                 angle opposite side a, B the angle opposite side b and C
%                 the angle opposite side C. The dimensions correspond to
%                 those of the inputs.
%
%The spherical law of cosines for the sides of a spherical triangle is
%stated in [1].
%
%EXAMPLE:
%First, we choose 3 points on the sphere, find the distances between them
%and then show that we can get the angles. Then, we just set some arbitrary
%distances that do not form a possible triangle on the spehre and we see
%that the angular results are complex.
%Choose 3 points in term of latitude-longitude on a sphere.
% llPoints=deg2rad([19,  21, 27;
%                 -158,-147,-152]);
% r=1;
% a=greatCircleDistance(llPoints(:,1),llPoints(:,2),r);
% b=greatCircleDistance(llPoints(:,2),llPoints(:,3),r);
% c=greatCircleDistance(llPoints(:,3),llPoints(:,1),r);
% [A,B,C]=spherTriAngles4Sides(a,b,c,r);
% ADeg=rad2deg(A)
% BDeg=rad2deg(B)
% CDeg=rad2deg(C)
% a=1;
% b=2;
% c=3;
% [A,B,C]=spherTriAngles4Sides(a,b,c,r)
%
%REFERENCES:
%[1] Weisstein, Eric W. "Law of Cosines." From MathWorld--A Wolfram Web
%    Resource. https://mathworld.wolfram.com/LawofCosines.html
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(r))
    r=1;
end

a=a./r;
b=b./r;
c=c./r;

cosa=cos(a);
sina=sin(a);
cosb=cos(b);
sinb=sin(b);
cosc=cos(c);
sinc=sin(c);

A=acos((cosa-cosb.*cosc)./(sinb.*sinc));
B=acos((cosb-cosa.*cosc)./(sina.*sinc));
C=acos((cosc-cosa.*cosb)./(sina.*sinb));

end

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


text:
%NOFRICTION	return link object with zero friction
%
%	LINK = NOFRICTION(LINK)
%
%

% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function  l2 = nofriction(l, only)
%Link.nofriction Remove friction
%
% LN = L.nofriction() is a link object with the same parameters as L except
% nonlinear (Coulomb) friction parameter is zero.
%
% LN = L.nofriction('all') is a link object with the same parameters as L
% except all friction parameters are zero.

	l2 = Link(l);
	if (nargin == 2) && strcmpi(only(1:3), 'all')
		l2.B = 0;
	end
	l2.Tc = [0 0];
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function value = r8_shi ( x )

%*****************************************************************************80
%
%% R8_SHI evaluates the hyperbolic sine integral Shi of an R8 argument.
%
%  Discussion:
%
%    Shi ( x ) = Integral ( 0 <= t <= x ) sinh ( t ) dt / t
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the hyperbolic sine integral
%    Shi evaluated at X.
%
  persistent nshi
  persistent shics
  persistent xsml

  if ( isempty ( nshi ) )

    shics = [ ...
      0.0078372685688900950695200984317332, ...
      0.0039227664934234563972697574427225, ...
      0.0000041346787887617266746747908275, ...
      0.0000000024707480372882742135145302, ...
      0.0000000000009379295590763630457157, ...
      0.0000000000000002451817019520867353, ...
      0.0000000000000000000467416155257592, ...
      0.0000000000000000000000067803072389, ...
      0.0000000000000000000000000007731289, ...
      0.0000000000000000000000000000000711 ]';

    nshi = r8_inits ( shics, 10, 0.1 * r8_mach ( 3 ) );
    xsml = sqrt ( r8_mach ( 3 ) );

  end

  absx = abs ( x );

  if ( absx <= xsml )
    value = x;
  elseif ( absx <= 0.375 )
    value = x * ( 1.0 + r8_csevl ( 128.0 * x * x / 9.0 - 1.0, shics, nshi ) );
  else
    value = 0.5 * ( r8_ei ( x ) + r8_e1 ( x ) );
  end

  return
end

repo: jburkardt-m


text:
function suborder = lyness_suborder ( rule, suborder_num )

%*****************************************************************************80
%
%% LYNESS_SUBORDER returns the suborders for a Lyness rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    James Lyness, Dennis Jespersen,
%    Moderate Degree Symmetric Quadrature Rules for the Triangle,
%    Journal of the Institute of Mathematics and its Applications,
%    Volume 15, Number 1, February 1975, pages 19-32.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%
%    Input, integer SUBORDER_NUM, the number of suborders
%    of the rule.
%
%    Output, integer SUBORDER(SUBORDER_NUM), the suborders
%    of the rule.
%
  if ( rule == 0 )
    suborder = [ 1 ]';
  elseif ( rule == 1 )
    suborder = [ 3 ]';
  elseif ( rule == 2 )
    suborder = [ 1, 3 ]';
  elseif ( rule == 3 )
    suborder = [ 1, 3 ]';
  elseif ( rule == 4 )
    suborder = [ 1, 3, 3 ]';
  elseif ( rule == 5 )
    suborder = [ 3, 3  ]';
  elseif ( rule == 6 )
    suborder = [ 1, 3, 6  ]';
  elseif ( rule == 7 )
    suborder = [ 3, 3, 3 ]';
  elseif ( rule == 8 )
    suborder = [ 1, 3, 3 ]';
  elseif ( rule == 9 )
    suborder = [ 1, 3, 3, 3 ]';
  elseif ( rule == 10 )
    suborder = [ 3, 3, 6 ]';
  elseif ( rule == 11 )
    suborder = [ 1, 3, 3, 3, 6 ]';
  elseif ( rule == 12 )
    suborder = [ 1, 3, 3, 6 ]';
  elseif ( rule == 13 )
    suborder = [ 1, 3, 3, 6 ]';
  elseif ( rule == 14 )
    suborder = [ 1, 3, 3, 3, 6 ]';
  elseif ( rule == 15 )
    suborder = [ 1, 3, 3, 3, 6 ]';
  elseif ( rule == 16 )
    suborder = [ 3, 3, 3, 3, 3, 6 ]';
  elseif ( rule == 17 )
    suborder = [ 1, 3, 3, 3, 6 ]';
  elseif ( rule == 18 )
    suborder = [ 1, 3, 3, 3, 3, 6 ]';
  elseif ( rule == 19 )
    suborder = [ 1, 3, 3, 3, 3, 3, 6 ]';
  elseif ( rule == 20 )
    suborder = [ 3, 3, 3, 3, 3, 6, 6 ]';
  elseif ( rule == 21 )
    suborder = [ 1, 3, 3, 3, 3, 3, 6, 6 ]';
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LYNESS_SUBORDER - Fatal error!\n' );
    fprintf ( 1, '  Illegal RULE = %d\n', rule );
    error ( 'LYNESS_SUBORDER - Fatal error!' );
  end

  return
end


text:
function [ imat, mat, ifault ] = invmod ( mat, rmod, cmod, nrow, ifault )

%*****************************************************************************80
%
%% INVMOD inverts a matrix using modulo arithmetic.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 May 2013
%
%  Author:
%
%    Original FORTRAN77 version by Roger Payne.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Roger Payne,
%    Inversion of matrices with contents subject to modulo arithmetic,
%    Applied Statistics,
%    Volume 46, Number 2, 1997, pages 295-298.
%
%  Parameters:
%
%    Input/output, integer MAT(NROW*NROW).
%    On input, the matrix to be inverted.
%    On output, the product of the input matrix and IMAT.
%
%    Output, integer IMAT(NROW*NROW), the inverse matrix.
%    If IFAULT = -1 on output, then IMAT is only a left inverse.
%
%    Input, integer RMOD(NROW), the modulus for values in each row.
%
%    Input, integer CMOD(NROW), the modulus for values
%    in each column.
%
%    Input, integer NROW, the order of the matrix.
%
%    Output, integer IFAULT, an error flag.
%    0, no error was detected.
%    -1, only a left inverse could be formed.
%    1, the matrix contains elements that are negative, or too large.
%    2, the matrix contains nonzero elements in mixed modulus positions.
%    3, the matrix cannot be inverted.
%
  imat = [];
%
%  Check that elements in 'mixed-moduli' positions are all zero.
%
  n = 0;
  for i = 1 : nrow
    for j = 1 : nrow

      n = n + 1;

      if ( ( rmod(i) ~= cmod(j) ) && ( 0 < mat(n) ) )
        ifault = 2;
        return
      end

      if ( ( rmod(i) < mat(n) ) || ( mat(n) < 0 ) )
        ifault = 1;
        return
      end

    end
  end

  n = 0;
  for i = 1 : nrow
    for j = 1 : nrow
      n = n + 1;
      imat(n) = 0;
    end
  end
%
%  Sort rows and columns into ascending order of moduli
%
  [ mat, rmod, cmod, rsort, csort ] = msort ( mat, imat, rmod, cmod, nrow );
%
%  Complete initialization of inverse matrix
%
  for n = 1 : nrow + 1 : nrow * nrow
    imat(n) = 1;
  end
%
%  Invert the matrix.
%
  for ir = 1 : nrow

    kir = ( ir - 1 ) * nrow;

    if ( mat(kir+ir) == 0 )
%
%  Find a row JR below IR such that K(JR,IR)>0
%
      all_zero = 1;

      for kjr = kir + nrow + ir : nrow : nrow * nrow
        if ( 0 < mat(kjr) )
          all_zero = 0;
          break;
        end
      end
%
%  Column IR contains all zeros in rows IR or below:
%  look for a row above with zeros to left of column IR
%  and K(JR,IR)>0
%
      if ( all_zero )
        for kjr = ir : nrow : kir
          if ( 0 < mat(kjr) )
            for i = kjr - ir + 1 : kjr - 1
              if ( 0 < mat(i) )
                ifault = 3;
                return
              end
            end
            all_zero = 0;
            break
          end
        end
      end
%
%  Column IR contains all zeros
%
      if ( all_zero )
        continue;
      end
%
%  Switch row JR with row IR
%
      kjr = kjr - ir;

      for i = 1 : nrow

        k = mat(kir+i);
        mat(kir+i) = mat(kjr+i);
        mat(kjr+i) = k;

        k = imat(kir+i);
        imat(kir+i) = imat(kjr+i);
        imat(kjr+i) = k;

      end

    end
%
%  Find a multiplier N such that N*MAT(IR,IR)=1 mod(P{IR})
%
    k = mat(kir+ir);
    for n = 1 : rmod(ir) - 1
      if ( mod ( n * k, rmod(ir) ) == 1 )
        break;
      end
    end
%
%  Multiply row IR by N.
%
    if ( 1 < n )
      for i = kir + 1 : ir * nrow
        mat(i) = mat(i) * n;
        imat(i) = imat(i) * n;
      end
    end
%
%  Subtract MAT(JR,IR) * row IR from each row JR
%
    for kjr = 0 : nrow : nrow * nrow - 1
      n = rmod(ir) - mat(kjr+ir);
      if ( ( kjr ~= kir ) && ( n ~= 0 ) )
        for i = 1 : nrow
          mat(kjr+i)  = mod (  mat(kjr+i) + n *  mat(kir+i), cmod(i) );
          imat(kjr+i) = mod ( imat(kjr+i) + n * imat(kir+i), cmod(i) );
        end
      end
    end

  end
%
%  Check inversion was possible - that result has
%  non-zero elements only on diagonal.
%
  ifault = 0;
%
%  If we encounter a zero diagonal element, then only a left inverse
%  will be formed.
%
  for n = 1 : nrow + 1 : nrow * nrow
    if ( mat(n) == 0 )
      ifault = -1;
    end
    mat(n) = - mat(n);
  end

  for n = 1 : nrow * nrow
    if ( 0 < mat(n) )
      ifault = 3;
      return
    end
  end

  for n = 1 : nrow + 1 : nrow * nrow
    mat(n) = - mat(n);
  end
%
%  Unsort the rows and columns back into their original order.
%
  [ mat, imat, rmod, cmod, rsort, csort ] = musort ( mat, imat, rmod, ...
    cmod, rsort, csort, nrow );

  return
end


text:
function BC = createBC3D(meshvar)
% Creates a boundary condition structure from a mesh structure
% for a 3D structured mesh. The default boundary conditions on
% all boundaries are Neumann;
% The values of each boundary condition are defined as:
% BC.<position>. :
%               a, b, c, where
% a*grad(phi).e + b*phi = c
%
% SYNOPSIS:
%   BC = createBC3D(meshvar)
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Extract number of cells from the mesh structure
Nxyz = meshvar.dims;
Nx = Nxyz(1); Ny = Nxyz(2); Nz = Nxyz(3);

% Define the top, bottom, right, and left boundary conditions
% (default = Neumann, i.e., a = 1, b = 0, c = 0)
top.a = ones(Nx,Nz);
top.b = zeros(Nx,Nz);
top.c = zeros(Nx,Nz);
top.periodic = 0;

bottom.a = ones(Nx,Nz);
bottom.b = zeros(Nx,Nz);
bottom.c = zeros(Nx,Nz);
bottom.periodic = 0;

right.a = ones(Ny,Nz);
right.b = zeros(Ny,Nz);
right.c = zeros(Ny,Nz);
right.periodic = 0;

left.a = ones(Ny,Nz);
left.b = zeros(Ny,Nz);
left.c = zeros(Ny,Nz);
left.periodic = 0;

front.a = ones(Nx,Ny);
front.b = zeros(Nx,Ny);
front.c = zeros(Nx,Ny);
front.periodic = 0;

back.a = ones(Nx,Ny);
back.b = zeros(Nx,Ny);
back.c = zeros(Nx,Ny);
back.periodic = 0;

BC= BoundaryCondition(meshvar, left, right, bottom, top, back, front);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
%% ------ monetary policy variance and coefs switch independently ------ %%
%
% Both coefficients and variances in monetary policy equation
% change with two independent Markov processes

%% housekeeping
clear
close all
clc()

%% Create dataset
clc

do_plot=true;

scale=100;

[db,varlist0]=create_dataset(scale,do_plot);

varlist=fieldnames(varlist0);

%% set up Markov chains

markov_chains=struct('name','mpcoef',...
    'number_of_states',2,...
    'controlled_parameters',{{'c(1)','a0(1)','a1(1)','a2(1)'}},...
    'endogenous_probabilities',[],...
    'probability_parameters',[]);

markov_chains(2)=struct('name','mpvol',...
    'number_of_states',3,...
    'controlled_parameters',{{'s(1)'}},...
    'endogenous_probabilities',[],...
    'probability_parameters',[]);

%% Create the VAR

clc

nlags=2;

exog={};

constant=true;

panel=[];

sv0=svar(varlist,exog,nlags,constant,panel,markov_chains);

%% set up restrictions

% syntax is alag(eqtn,vname,chain_name,state)
%------------------------------------------------
lin_restr=cell(0,1);

nonlin_restr={'a0(3,FFR)>=0'};

which_chain=strcmp('mpcoef',{markov_chains.name});
% first equation or "FFR" equation:
%----------------------------------
for istate=1:markov_chains(which_chain).number_of_states

    mystate=int2str(istate);

    lin_restr=[lin_restr
        {
        ['a1(1,pi,mpcoef,',mystate,')=0']
        ['a2(1,pi,mpcoef,',mystate,')=0']
        ['a1(1,ygap,mpcoef,',mystate,')=0']
        ['a2(1,ygap,mpcoef,',mystate,')=0']
        ['a2(1,FFR,mpcoef,',mystate,')=0']
        }
        ]; %#ok<AGROW>

    nonlin_restr=[nonlin_restr
        {
        ['a1(1,FFR,mpcoef,',mystate,')>=0']
        ['a1(1,FFR,mpcoef,',mystate,')<=1']
        }
        ]; %#ok<AGROW>
end

lin_restr=[lin_restr
    {
    % second equation or "pi" equation
    %----------------------------------
    'a0(2,FFR)=0'
    'a1(2,FFR)=0'
    'a2(2,FFR)=0'
    'a1(2,ygap)=0'
    'a2(2,ygap)=0'
    % third equation or "ygap" equation
    %-----------------------------------
    'a1(3,FFR)=0'
    'a2(3,FFR)=0'
    'a1(3,pi)=0'
    'a2(3,pi)=0'
    'a0(3,pi)+a0(3,FFR)=0'
    }
    ];

restrictions=[lin_restr;nonlin_restr];

%% set priors

% priors for the VAR coefficients
%--------------------------------
var_prior=svar.prior_template();
var_prior.type='sz';
% priors for the mpcoef transition probabilities
%------------------------------------------------
switch_prior=struct();
switch_prior.mpcoef_tp_1_2={0.5,0.1,0.3,'beta'};
switch_prior.mpcoef_tp_2_1={0.5,0.1,0.3,'beta'};
% priors for the mpvol transition probabilities
%------------------------------------------------
switch_prior.dirichlet_1={0.1,'mpvol_tp_1_2',0.2,'mpvol_tp_1_3',0.2};
switch_prior.dirichlet_2={0.1,'mpvol_tp_2_1',0.2,'mpvol_tp_2_3',0.2};
switch_prior.dirichlet_3={0.1,'mpvol_tp_3_1',0.2,'mpvol_tp_3_2',0.2};

prior=struct();

prior.var=var_prior;

prior.nonvar=switch_prior;

%% Find posterior mode
clc

sv=sv0;

sv=estimate(sv,db,{'1960Q1','2015Q2'},prior,restrictions);

%% estimates
clc

pmode=posterior_mode(sv)

%% Printing estimates
clc

print_structural_form(sv)

%% Printing solution
clc

print_solution(sv)

%% plot smoothed state and regime probabilities
clc

close all

plot_probabilities(sv)

%% plots probabilities against data
close all

plot_data_against_probabilities(sv,'regime')

%% Impulse responses

myirfs=irf(sv);

%% Posterior sampling

%% Marginal data density

%% Out-of sample forecasts

%% Conditional forecast

author: jmaih


text:
%% Ex. 11 Another example for the usage of index for a matrix

A = [3 5; 2 4];
norm1 = 0;
for m = 1:2
for n = 1:2
 norm1 = norm1+A(m,n)^2;
end
end
norm1 = sqrt(norm1)

%Output:
%    norm1 = 7.348
% Remark: This program calculates the Euclidean norm of the A matrix

author: TheAlgorithms

repo: MATLAB-Octave


text:
clear all; close all;
I=imread('cameraman.tif');
I=im2double(I);
LEN=21;
THETA=11;
PSF=fspecial('motion', LEN, THETA);
J=imfilter(I, PSF, 'conv', 'circular');
noise_mean=0;
noise_var=0.0001;
K=imnoise(J, 'gaussian', noise_mean, noise_var);
figure;
subplot(121);  imshow(I);
subplot(122);  imshow(K);
NSR1=0;
L1=deconvwnr(K, PSF, NSR1);
NSR2=noise_var/var(I(:));
L2=deconvwnr(K, PSF, NSR2);
figure;
subplot(121);  imshow(L1);
subplot(122);  imshow(L2);

author: HuangCongQing


text:
function boolvec = artifact2boolvec(artifact, varargin)

% ARTIFACT2BOOLVEC converts between two representations of events or trials.
%
% FieldTrip uses a number of representations for events that are conceptually very similar
%   event    = structure with type, value, sample, duration and offset
%   trl      = Nx3 numerical array with begsample, endsample, offset
%   trl      = table with 3 columns for begsample, endsample, offset
%   artifact = Nx2 numerical array with begsample, endsample
%   artifact = table with 2 columns for begsample, endsample
%   boolvec  = 1xNsamples boolean vector with a thresholded TTL/trigger sequence
%   boolvec  = MxNsamples matrix vector with a thresholded TTL/trigger sequence
%
% If trl or artifact are represented as a MATLAB table, they can have additional
% columns. These additional columns have to be named and are not restricted to
% numerical values.
%
% This function makes a Boolean vector (or matrix when artifact is a cell-array of
% multiple artifact definitions) with 0 for artifact free sample and 1 for sample
% containing an artifact according to artifact specification. The length of the
% vector matches the last sample in the artifact definition, or endsample when
% specified.
%
% See also ARTIFACT2BOOLVEC, ARTIFACT2EVENT, ARTIFACT2TRL, BOOLVEC2ARTIFACT, BOOLVEC2EVENT, BOOLVEC2TRL, EVENT2ARTIFACT, EVENT2BOOLVEC, EVENT2TRL, TRL2ARTIFACT, TRL2BOOLVEC, TRL2EVENT

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

% trl and artifact are similar, except for the offset column
trl = artifact2trl(artifact);
boolvec = trl2boolvec(trl, varargin{:});


text:
%Calling function of the visualization functions
close all
clear all
colors={'r.' 'gx' 'b+' 'ys' 'm.' 'c.' 'k.' 'r*' 'g*' 'b*' 'y*' 'm*' 'c*' 'k*' };
path(path,'..\..\FUZZCLUST')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Selecting the data set
wine=1;
iris=0;
wisc=0;

if wine
    load winedat.txt
    data=winedat(:,1:end-1);
    C=winedat(:,end);
end

if iris
    load iris
    data=iris(:,1:4);
    C=zeros(length(data),1);
    for i=1:3
        C(find(iris(:,4+i)==1))=i;
    end
end
if wisc
   %Data preprocessing
    wisc=wk1read('wisconsin.wk1');
    NI=9;
    NT=length(wisc);
    data.X=[wisc(:,11) wisc(:,2:10)];
    data.X=sortrows(data.X,1);
    [I,J]=find(data.X(:,7)~=0);
    data.X=data.X(I,:);
    [I,J]=find(data.X(:,1)==2);
    data.X(I,1)=1;
    [I,J]=find(data.X(:,1)==4);
    data.X(I,1)=2;
    C=data.X(:,1);
    data=data.X(:,2:end);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %normalization of the data
data.X=data;
data=clust_normalize(data,'range');
%fuzzy c-means clustering
param.m=2;
param.c=2;
param.val=1;
param.vis=0;
result=Kmeans(data,param);
result=validity(result,data,param);
%Assignment for classification
[d1,d2]=max(result.data.f');
Cc=[];
for i=1:param.c
    Ci=C(find(d2==i));
    dum1=hist(Ci,1:param.c);
    [dd1,dd2]=max(dum1);
    Cc(i)=dd2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Principal Component Projection of the data and the cluster centers
param.q=2;
result = PCA(data,param,result);
%visualization
figure(1)
clf
for i=1:max(C)
    index=find(C==i);
    err=(Cc(d2(index))~=i);
    eindex=find(err);
    misclass(i)=sum(err);
    plot(result.proj.P(index,1),result.proj.P(index,2),[colors{i}])
    hold on
    plot(result.proj.P(index(eindex),1),result.proj.P(index(eindex),2),'o')
    hold on
end
    xlabel('y_1')
    ylabel('y_2')
    title('PCA projection')
%The error value of classification
perfclass=sum(misclass)/length(C)*100
%
plot(result.proj.vp(:,1),result.proj.vp(:,2),'r*');
%calculating realtion-indexes
result = samstr(data,result);
perf = [projeval(result,param) result.proj.e];
%
disp('Press any key.')
pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SAMMON mapping

proj.P=result.proj.P;   %Sammon uses the output of PCA for initializing
param.alpha = 0.4;
param.max=100;

figure(2)
result = Sammon(proj,data,result,param)
%visualization
clf
for i=1:max(C)
    index=find(C==i);
    err=(Cc(d2(index))~=i);
    eindex=find(err);
    misclass(i)=sum(err);
    plot(result.proj.P(index,1),result.proj.P(index,2),[colors{i}] )
    hold on
    plot(result.proj.P(index(eindex),1),result.proj.P(index(eindex),2),'o')
    hold on
end
    xlabel('y_1')
    ylabel('y_2')
    title('Conventional Sammon mapping')
%
plot(result.proj.vp(:,1),result.proj.vp(:,2),'r*');
%calculating realtion-indexes
result = samstr(data,result);
perfs = [projeval(result,param) result.proj.e];
%
disp('Press any key.')
pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Modified fuzzy SAMMON mapping

proj.P=result.proj.P; %FuzSam uses the output of Sammon for initializing
param.alpha = 0.4;
param.max=100;

figure(3)
result=FuzSam(proj,result,param);

clf
for i=1:max(C)
    index=find(C==i);
    err=(Cc(d2(index))~=i);
    eindex=find(err);
    misclass(i)=sum(err);
    plot(result.proj.P(index,1),result.proj.P(index,2),[colors{i}] )
    hold on
    plot(result.proj.P(index(eindex),1),result.proj.P(index(eindex),2),'o')
    hold on
end
    xlabel('y_1')
    ylabel('y_2')
    title('Fuzzy Sammon mapping')

plot(result.proj.vp(:,1),result.proj.vp(:,2),'r*')
%calculating realtion-indexes
result = samstr(data,result);
perff = [projeval(result,param) result.proj.e];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
perf
perfs
perff

author: Sable

repo: mcbench-benchmarks


text:
function h=hermfun(t,j)
%HERMFUN  Orthonormal Hermite functions. [with F. Rekibi]
%
%   H=HERMFUN(T,N) generates the fisrt N+1 orthonormal Hermite functions
%   [H0,...HN] on a time axis specfied by the column vector T.
%
%   HERMFUN uses the expression of Simons et al. 2003.
%
%   Note that H(:,1) is the zeroth-order Hermite function, which is equal
%   to a Gaussian.  H(:,2) is the first-order function, and so forth.
%
%   See also HERMPOLY.
%
%   'hermfun --f' generates a sample figure; compare with the Hermite
%   function figure at
%
%     http://en.wikipedia.org/wiki/Hermite_polynomials#Definition
%
%   Usage:  h=hermfun(t,n);
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2004--2015 F. Rekibi and J. M. Lilly
%                         --- type 'help jlab_license' for details

% 05.08.07  JML fixed bug to include N+1 columns

%   'hermfun --f' generates a sample figure; compare with Figure 2
%   of Simons, van der Hilst, and Zuber (2003), JGR 108 B5.

if strcmpi(t,'--f')
  type makefigs_hermfun
  makefigs_hermfun;
  return
end

if size(t,1)==1
    t=t';
end

H=hermpoly(t,j);

E=exp(-t.^2/2)*ones(1,j+1);
HE=H.*E;

h=zeros(length(t),j);
for k=1:j+1
	h(:,k)=HE(:,k)*frac(1,pi^(1/4))*frac(1,sqrt(2^(k-1)*factorial(k-1)));
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
%ISHOMOG Test if SE(3) homogeneous transformation
%
% ISHOMOG(T) is true (1) if the argument T is of dimension 4x4 or 4x4xN, else
% false (0).
%
% ISHOMOG(T, 'valid') as above, but also checks the validity of the rotation
% sub-matrix.
%
% Notes::
% - The first form is a fast, but incomplete, test for a transform is SE(3).
% - Does not work for the SE(2) case.
%
% See also ISROT, ISVEC.

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

function h = ishomog(tr, rtest)
    d = size(tr);
    if ndims(tr) >= 2
        h =  all(d(1:2) == [4 4]);

        if h && nargin > 1
            h = abs(det(tr(1:3,1:3)) - 1) < eps;
        end

    else
        h = false;
    end

author: Allopart

repo: rbpf-gmapping


text:
classdef BankAngleConstraint < AbstractConstraint
    %BankAngleConstraint Summary of this class goes here
    %   Detailed explanation goes here

    properties
        normFact = 1;
        event LaunchVehicleEvent
        eventNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        evalType(1,1) ConstraintEvalTypeEnum = ConstraintEvalTypeEnum.FixedBounds;
        stateCompType(1,1) ConstraintStateComparisonTypeEnum = ConstraintStateComparisonTypeEnum.Equals;
        stateCompEvent LaunchVehicleEvent
        stateCompNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;
    end

    methods
        function obj = BankAngleConstraint(event, lb, ub)
            obj.event = event;
            obj.lb = lb;
            obj.ub = ub;

             obj.id = rand();
        end

        function [lb, ub] = getBounds(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function [c, ceq, value, lwrBnd, uprBnd, type, eventNum, valueStateComp] = evalConstraint(obj, stateLog, celBodyData)
            type = obj.getConstraintType();

            switch obj.eventNode
                case ConstraintStateComparisonNodeEnum.FinalState
                    stateLogEntry = stateLog.getLastStateLogForEvent(obj.event);

                case ConstraintStateComparisonNodeEnum.InitialState
                    stateLogEntry = stateLog.getFirstStateLogForEvent(obj.event);

                otherwise
                    error('Unknown event node.');
            end

            ut = stateLogEntry.time;
            rVect = stateLogEntry.position;
            vVect = stateLogEntry.velocity;
            bodyInfo = stateLogEntry.centralBody;

            [bankAng,~,~] = stateLogEntry.attitude.getAeroAngles(ut, rVect, vVect, bodyInfo);
            value = rad2deg(bankAng);

            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                switch obj.stateCompNode
                    case ConstraintStateComparisonNodeEnum.FinalState
                        stateLogEntryStateComp = stateLog.getLastStateLogForEvent(obj.stateCompEvent).deepCopy();

                    case ConstraintStateComparisonNodeEnum.InitialState
                        stateLogEntryStateComp = stateLog.getFirstStateLogForEvent(obj.stateCompEvent).deepCopy();

                    otherwise
                        error('Unknown event node.');
                end

                cartElem = stateLogEntryStateComp.getCartesianElementSetRepresentation();
                cartElem = cartElem.convertToFrame(stateLogEntry.centralBody.getBodyCenteredInertialFrame());
                stateLogEntryStateComp.setCartesianElementSet(cartElem);

                ut = stateLogEntryStateComp.time;
                rVect = stateLogEntryStateComp.position;
                vVect = stateLogEntryStateComp.velocity;
                bodyInfo = stateLogEntryStateComp.centralBody;

                [bankAng,~,~] = stateLogEntryStateComp.attitude.getAeroAngles(ut, rVect, vVect, bodyInfo);
                valueStateComp = rad2deg(bankAng);
            else
                valueStateComp = NaN;
            end

            [c, ceq] = obj.computeCAndCeqValues(value, valueStateComp);

            lwrBnd = obj.lb;
            uprBnd = obj.ub;

            eventNum = obj.event.getEventNum();
        end

        function sF = getScaleFactor(obj)
            sF = obj.normFact;
        end

        function setScaleFactor(obj, sF)
            obj.normFact = sF;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesEvent(obj, event)
            tf = obj.event == event;
            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                tf = tf || obj.stateCompEvent == event;
            end
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end

        function tf = usesExtremum(obj, extremum)
            tf = false;
        end

        function tf = canUseSparseOutput(obj)
            tf = true;
        end

        function event = getConstraintEvent(obj)
            event = obj.event;
        end

        function type = getConstraintType(obj)
            type = 'Bank Angle';
        end

%         function name = getName(obj)
%             name = sprintf('%s - Event %i', obj.getConstraintType(), obj.event.getEventNum());
%         end

        function [unit, lbLim, ubLim, usesLbUb, usesCelBody, usesRefSc] = getConstraintStaticDetails(obj)
            unit = 'deg';
            lbLim = -360;
            ubLim = 360;
            usesLbUb = true;
            usesCelBody = false;
            usesRefSc = false;
        end

        function addConstraintTf = openEditConstraintUI(obj, lvdData)
%             addConstraintTf = lvd_EditGenericMAConstraintGUI(obj, lvdData);

            output = AppDesignerGUIOutput({false});
            lvd_EditGenericMAConstraintGUI_App(obj, lvdData, output);
            addConstraintTf = output.output{1};
        end
    end

    methods(Static)
        function constraint = getDefaultConstraint(~, ~)
            constraint = BankAngleConstraint(LaunchVehicleEvent.empty(1,0),0,0);
        end
    end
end

author: Arrowstar

repo: ksptot