text:
function [D,P] = dijk(A,s,t)
%DIJK Shortest paths from nodes 's' to nodes 't' using Dijkstra algorithm.
% D = dijk(A,s,t)
%     A = n x n node-node weighted adjacency matrix of arc lengths
%         (Note: A(i,j) = 0   => Arc (i,j) does not exist;
%                A(i,j) = NaN => Arc (i,j) exists with 0 weight)
%     s = FROM node indices
%       = [] (default), paths from all nodes
%     t = TO node indices
%       = [] (default), paths to all nodes
%     D = |s| x |t| matrix of shortest path distances from 's' to 't'
%       = [D(i,j)], where D(i,j) = distance from node 'i' to node 'j'
%
%	(If A is a triangular matrix, then computationally intensive node
%   selection step not needed since graph is acyclic (triangularity is a
%   sufficient, but not a necessary, condition for a graph to be acyclic)
%   and A can have non-negative elements)
%
%	(If |s| >> |t|, then DIJK is faster if DIJK(A',t,s) used, where D is now
%   transposed and P now represents successor indices)
%
%  (Based on Fig. 4.6 in Ahuja, Magnanti, and Orlin, Network Flows,
%   Prentice-Hall, 1993, p. 109.)

% Copyright (c) 1998-2000 by Michael G. Kay
% Matlog Version 1.3 29-Aug-2000
%
%  Modified by JBT, Dec 2000, to delete paths

% Input Error Checking ******************************************************
error(nargchk(1,3,nargin));

[n,cA] = size(A);

if nargin < 2 | isempty(s), s = (1:n)'; else s = s(:); end
if nargin < 3 | isempty(t), t = (1:n)'; else t = t(:); end

if ~any(any(tril(A) ~= 0))			% A is upper triangular
   isAcyclic = 1;
elseif ~any(any(triu(A) ~= 0))	% A is lower triangular
   isAcyclic = 2;
else										% Graph may not be acyclic
   isAcyclic = 0;
end

if n ~= cA
   error('A must be a square matrix');
elseif ~isAcyclic & any(any(A < 0))
   error('A must be non-negative');
elseif any(s < 1 | s > n)
   error(['''s'' must be an integer between 1 and ',num2str(n)]);
elseif any(t < 1 | t > n)
   error(['''t'' must be an integer between 1 and ',num2str(n)]);
end
% End (Input Error Checking) ************************************************

A = A';		% Use transpose to speed-up FIND for sparse A

D = zeros(length(s),length(t));
P = zeros(length(s),n);

for i = 1:length(s)
   j = s(i);

   Di = Inf*ones(n,1); Di(j) = 0;

   isLab = logical(zeros(length(t),1));
   if isAcyclic ==  1
      nLab = j - 1;
   elseif isAcyclic == 2
      nLab = n - j;
   else
      nLab = 0;
      UnLab = 1:n;
      isUnLab = logical(ones(n,1));
   end

   while nLab < n & ~all(isLab)
      if isAcyclic
         Dj = Di(j);
      else	% Node selection
         [Dj,jj] = min(Di(isUnLab));
         j = UnLab(jj);
         UnLab(jj) = [];
         isUnLab(j) = 0;
      end

      nLab = nLab + 1;
      if length(t) < n, isLab = isLab | (j == t); end

      [jA,kA,Aj] = find(A(:,j));
      Aj(isnan(Aj)) = 0;

      if isempty(Aj), Dk = Inf; else Dk = Dj + Aj; end

      P(i,jA(Dk < Di(jA))) = j;
      Di(jA) = min(Di(jA),Dk);

      if isAcyclic == 1			% Increment node index for upper triangular A
         j = j + 1;
      elseif isAcyclic == 2	% Decrement node index for lower triangular A
         j = j - 1;
      end

      %disp( num2str( nLab ));
   end
   D(i,:) = Di(t)';
end


text:
function filename_final=mesh2vrml(X, Y, Z)

% Transforming data into VRML. Remember Z axis is pointing out of the
% screen towards you. X axis to the right and Y Axis is upwards.
xDimension=length(X(1,:));
xSpacing=abs(X(1,2)-X(1,1));
zDimension=length(Y(1,:));
zSpacing=abs(Y(2,1)-Y(1,1));
filename='myshape.wrl';
filename_final='elevation.wrl';
fid = fopen(filename, 'w+','n','UTF-8');
fprintf(fid,'#VRML V2.0 utf8');
fprintf(fid, '\n Transform { \n');
fprintf(fid,' \n translation %8.4f %8.4f %8.4f\n', -(max(max(X))-min(min(X)))/2, 0, -(max(max(Y))-min(min(Y)))/2);
fprintf(fid,' \n  rotation 0 0 1 0\n');
fprintf(fid,' \n  center %8.4f %8.4f %8.4f\n', -(max(max(X))+min(min(X)))/2, -(max(max(Z))+min(min(Z)))/2, 0);
fprintf(fid,' \n  scale 1 1 1\n');
fprintf(fid, '\n children Shape {\n');
fprintf(fid,'\n appearance	Appearance {\n');
fprintf(fid,'\n material	Material {}\n');
fprintf(fid,'\n	}\n');
fprintf(fid, '\n geometry	ElevationGrid { \n');
fprintf(fid, '\n xDimension %d\n', xDimension);
fprintf(fid, '\n xSpacing %d\n', xSpacing);
fprintf(fid, '\n zDimension %d\n', zDimension);
fprintf(fid, '\n zSpacing %d\n', zSpacing);
fprintf(fid,' height [');
for j=1:1:zDimension
    for i=1:1:xDimension
        fprintf(fid,'%f', Z(i,j));
        if (i==xDimension && j==zDimension)
        else
            fprintf(fid,',');
        end
    end
    fprintf(fid,'\n');
end
fprintf(fid,']\n');
fprintf(fid,'\n color NULL\n');
fprintf(fid,'\n colorPerVertex TRUE\n');
fprintf(fid,'\n normal NULL\n');
fprintf(fid,'\n normalPerVertex TRUE\n');
fprintf(fid, '\n texCoord NULL \n');
fprintf(fid, '\n ccw TRUE\n');
fprintf(fid, '\n solid FALSE\n');
fprintf(fid, '\n creaseAngle 0.0\n');
fprintf(fid,'		}\n');
fprintf(fid,'		}\n');
fprintf(fid,'		}\n');
fclose(fid);
copyfile(filename,filename_final, 'f');
delete(filename);

% Visualization code
figure_height=400;
view_width=350;
corner_x=20;
corner_y=10;
clearance=20;
fig=figure('Name', 'Preview', 'Position',[100 450 2*view_width+3*clearance figure_height+clearance]);
% Create the MATLAB GUI with two views of the aircraft
% First create the vrcanvas object for the first view. Specify the location
% and size of the view in pixels
% Create the plot of the translational coordinates of the aircraft
h = axes('Units', 'pixels', 'OuterPosition',[corner_x+view_width+3*clearance corner_y view_width figure_height]);
set(h, 'Position', [corner_x+clearance corner_y view_width figure_height]);
set(h, 'FontSize', 8);
surf(X,Y,Z);
hold on;
grid on;
title('Surface in MATLAB');
xlabel(h, 'x')
ylabel(h, 'y')
zlabel(h, 'z')
set(h, 'Units', 'normalized');
view(h, [90 0]);
world=vrworld(filename_final);
open(world);
c1 = vr.canvas(world, fig,[corner_x++view_width+3*clearance corner_y+clearance view_width-2*clearance  figure_height-2*clearance]);
set(c1, 'Units', 'normalized');

author: Sable

repo: mcbench-benchmarks


text:
function unique_num = r8col_tol_unique_count ( m, n, a, tol )

%*****************************************************************************80
%
%% R8COL_TOL_UNIQUE_COUNT counts tolerably unique entries in an R8COL.
%
%  Discussion:
%
%    Because the array is unsorted, this algorithm is O(N^2).
%
%    If the tolerance is large enough, then the concept of uniqueness
%    can become ambiguous.  If we have a tolerance of 1.5, then in the
%    list ( 1, 2, 3, 4, 5, 6, 7, 8, 9 ) is it fair to say we have only
%    one unique entry?  That would be because 1 may be regarded as unique,
%    and then 2 is too close to 1 to be unique, and 3 is too close to 2 to
%    be unique and so on.
%
%    This seems wrongheaded.  So I prefer the idea that an item is not
%    unique under a tolerance only if it is close to something that IS unique.
%    Thus, the unique items are guaranteed to cover the space if we include
%    a disk of radius TOL around each one.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 July 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows.
%
%    Input, integer N, the number of columns.
%
%    Input, real A(M,N), the unsorted array.
%
%    Input, real TOL, a tolerance for equality.
%
%    Output, integer UNIQUE_NUM, the number of unique columns of A.
%
  undx = zeros ( n, 1 );
%
%  Implicitly sort the array.
%
  indx = r8col_sort_heap_index_a ( m, n, a );
%
%  Consider entry I = 1.
%  It is unique, so set the number of unique items to K.
%  Set the K-th unique item to I.
%  Set the representative of item I to the K-th unique item.
%
  i = 1;
  k = 1;
  undx(k) = indx(i);
%
%  Consider entry I.
%
%  If it is unique, increase the unique count K, set the
%  K-th unique item to I, and set the representative of I to K.
%
%  If it is not unique, set the representative of item I to a
%  previously determined unique item that is close to it.
%
  for i = 2 : n

    unique = 1;

    for j = 1 : k
      diff = max ( abs ( a(1:m,indx(i)) - a(1:m,undx(j)) ) );
      if ( diff <= tol )
        unique = 0;
        break
      end
    end

    if ( unique )
      k = k + 1;
      undx(k) = indx(i);
    end

  end

  unique_num = k;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function featureReducedS = reduceDirFeatures(featuresAllDirS, reduceType)
% function featureReducedS = reduceDirFeatures(featuresAllDirS, reduceType)
%
% This function combines the directional feature values into a single
% value. The reduceType specifies how the different directions should be
% combined. The supported reduceTypes are 'avg', 'max', 'min', 'std',
% 'mad', 'iqr'
%
% APA, 4/11/2017

filedNamC = fieldnames(featuresAllDirS);

for iField = 1:length(filedNamC)
    fieldName = filedNamC{iField};
    switch lower(reduceType)
        case 'avg'
            featureReducedS.(fieldName) = mean([featuresAllDirS.(fieldName)]);
        case 'min'
            featureReducedS.(fieldName) = min([featuresAllDirS.(fieldName)]);
        case 'max'
            featureReducedS.(fieldName) = max([featuresAllDirS.(fieldName)]);
        case 'std'
            featureReducedS.(fieldName) = std([featuresAllDirS.(fieldName)]);
        case 'mad'
            featureReducedS.(fieldName) = mad([featuresAllDirS.(fieldName)]);
        case 'iqr'
            featureReducedS.(fieldName) = mad([featuresAllDirS.(fieldName)]);
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%% Copyright (C) 2014-2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym asin (@var{x})
%% Symbolic asin function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = asin (x)
%%   @result{} y = (sym) asin(x)
%% @end group
%% @end example
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @end defmethod

function y = asin(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('asin', x);
end

%!error asin (sym(1), 2)
%!assert (isequaln (asin (sym(nan)), sym(nan)))

%!shared x, d
%! d = 1;
%! x = sym('1');

%!test
%! f1 = asin(x);
%! f2 = asin(d);
%! assert( abs(double(f1) - f2) < 1e-15 )

%!test
%! D = [d d; d d];
%! A = [x x; x x];
%! f1 = asin(A);
%! f2 = asin(D);
%! assert( all(all( abs(double(f1) - f2) < 1e-15 )))

%!test
%! % round trip
%! y = sym('y');
%! A = asin (d);
%! f = asin (y);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B, -eps)

repo: octsympy


text:
function intrec(Z)
%INTREC Performs morphological character recognition
%INTREC(X) recognizes the characters (in integer form) present in the
%image X and displays it as output. INTREC uses morphological operations(like
%Dilation and Hit-or-Miss transform) to recognize characters.The charecter
%size of the integers present in the image should be exactly 26, otherwise
%it may not recognize it.
%
%Example:
%X=imread('testimage1.bmp');
%imshow(X)
%intrec(X)
%The digits found in the image are:
%0
%3
%5
%--------------------------------------------------------------------------
%Authors: Jahanzeb Rajput and Mohammad Fahad
%Department of Electrical Engineering
%University of Engineering and Technology, Lahore, Pakistan.
%--------------------------------------------------------------------------

A= imread ('zero.bmp');
B= imread ('one.bmp');
C= imread ('two.bmp');
D= imread ('three.bmp');
E= imread ('four.bmp');
F= imread ('five.bmp');
G= imread ('six.bmp');
H= imread ('seven.bmp');
I= imread ('eight.bmp');
J= imread ('nine.bmp');
%SE 2
SE = strel('square',3);

K=imdilate(A,SE);A2=K-A;

L=imdilate(B,SE);B2=L-B;

M=imdilate(C,SE);C2=M-C;

N=imdilate(D,SE);D2=N-D;

O=imdilate(E,SE);E2=O-E;

P=imdilate(F,SE);F2=P-F;

Q=imdilate(G,SE);G2=Q-G;

R=imdilate(H,SE);H2=R-H;

S=imdilate(I,SE);I2=S-I;

T=imdilate(J,SE);J2=T-J;
%-------------------------
%Hit or Miss
%-------------------------
disp('The digits found in the image are:');
if ~isempty(nonzeros(bwhitmiss(Z,A,A2)))
  disp('0');
end
%imshow('num1.bmp')
%U=bwhitmiss(Z,F,F2);
%figure
%imshow(U)
if ~isempty(nonzeros(bwhitmiss(Z,B,B2)))
    disp('1');
end

if ~isempty(nonzeros(bwhitmiss(Z,C,C2)))
   disp('2');
end

if ~isempty(nonzeros(bwhitmiss(Z,D,D2)))
    disp('3');
end

if ~isempty(nonzeros(bwhitmiss(Z,E,E2)))
   disp('4');
end

if ~isempty(nonzeros(bwhitmiss(Z,F,F2)))
   disp('5');
end

if ~isempty(nonzeros(bwhitmiss(Z,G,G2)))
   disp('6');
end

if ~isempty(nonzeros(bwhitmiss(Z,H,H2)))
   disp('7');
end

if ~isempty(nonzeros(bwhitmiss(Z,I,I2)))
   disp('8');
end

if ~isempty(nonzeros(bwhitmiss(Z,J,J2)))
   disp('9');
end


text:
function [u,P] = spm_fx_tfm_P(u,P)
% returns exogenous input and input dependent parameters
% FORMAT [u,P] = spm_fx_tfm_P(u,P)
%
% arguments:
% u  - inputs
% P  - parameters
%
% returns:
% u  - exogenous (conductance) inputs driving states
% P  - input dependent parameters
%
% tthis is a help are routine for the, microcircuit models equations of
% motion - it simply separates inputs into those affecting (driving) his
% neuronal states and those modulating parameters. It returns the exogenous
% (conductance) inputs and input dependent parameters.
%___________________________________________________________________________
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_fx_tfm_P.m 7679 2019-10-24 15:54:07Z spm $

% input dependent (intrinsic connection) parameters
%==========================================================================
j     = [4 3 2 1];
for i = 2:size(P.C,2)
    P.G(:,j(i - 1)) = P.G(:,j(i - 1)) + P.C(:,i) + u(i);
end

% exogenous inputs
%--------------------------------------------------------------------------
u     = exp(P.C(:,1))*u(1);

author: spm


text:
function neqn = p16_neqn ( )

%*****************************************************************************80
%
%% P16_NEQN returns the number of equations for problem P16.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NEQN, the number of equations.
%
  neqn = 4;

  return
end


text:
% alg_art1_test.m
% test ART etc.
% Copyright 2006-4-2, Jeff Fessler, University of Michigan

% generate data
if ~isvar('yi'), printm 'yi'
	ig = image_geom('nx', 16, 'ny', 14, 'dx', 1);
	ig.mask = ig.circ > 0;
	sg = sino_geom('par', 'nb', ig.nx, 'na', round(0.6 * ig.nx / 2)*2, ...
		'dr', 1);
	A = Gtomo2_strip(sg, ig);
	xtrue = ig.circ(5);
	yi = A * xtrue;
end

% ART1
if ~isvar('xart1'), printm 'art1'
	f.niter = 34;
	xinit = ig.zeros;

	xart1 = alg_art1(xinit(ig.mask), A', yi, 'niter', f.niter, 'isave', 'all');
	xart1 = ig.embed(xart1);
	im clf, im(xart1, 'ART1 iterates')
prompt
end

if ~isvar('yp')
	yp = A * xart1;
	resid = repmat(yi, [1 1 f.niter+1]) - yp;
end

	tmp = reshape(resid, [], f.niter+1);

if im
	im plc 3 3

	im(1, xtrue, 'xtrue')
	im(2, xart1, 'ART1')
	im(3, ig.mask, 'mask')
	im(4, yi, 'yi'), cbar
	im(5, yp(:,:,end), 'yp'), cbar
	im(6, yp(:,:,end)-yi, 'yp-yi'), cbar
	im(7, xart1(:,:,end), 'last iter'), cbar
	im(8, xart1(:,:,end)-xtrue, 'err'), cbar
	im subplot 9
	semilogy(0:f.niter, sqrt(mean(tmp.^2)))
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function identifiersDBs = getRegisteredDatabases()
% Get all databases registered with identifiers.org along with their
% pattern and and prefix data.
%
% USAGE:
%    identifiersDBs = getRegisteredDatabases()
%
% OUPUTS:
%    identifiersDBs:        A Struct with data on each database registered
%                           with identifiers.org (needs web access).
%                           The struct will have the following fields:
%                            * .name:       The name of the database
%                            * .pattern:    A regexp pattern allowed with the db.
%                            * .prefix:     The prefix of the db on identifiers.org

persistent databases

if isempty(databases)
    % load from external site
    try
        dbs = webread('https://registry.api.identifiers.org/resolutionApi/getResolverDataset');
    catch
        error('Could not load the databases registered with identifiers.org.\nThis is likely due to a missing internet connection.\nPlease try this again later');
    end
    dbs = dbs.payload.namespaces;
    % extract relevant information
    dbnames = {dbs.name};
    dbpatterns = {dbs.pattern};
    dbprefix = {dbs.prefix};
    % collate as struct
    databases = struct('name',dbnames,'pattern',dbpatterns,'prefix',dbprefix);
end

identifiersDBs = databases;

author: opencobra

repo: cobratoolbox


text:
function aboxes_out = post_process_candidate_detections_all_imgs(...
    aboxes_in, varargin)
% post_process_candidate_detections_all_imgs performs the post-processing
% step of non-maximum-suppression and optionally of box voting on the
% candidate bounding box detections of a set of images.
%
% INPUTS:
% 1) aboxes_in: contains the candidate bounding box detections of each image.
% It can be of two forms:
%   a) (the flag is_per_class must being set to false) a NI x 1 cell array
%   where NI is the number of images. The aboxes_in{i} element contains
%   the candidate bounding box detection of the i-th image in the form of a
%   NB_i x (4 + C) array where C is the number of categories and NB_i
%   is the number of candidate detection boxes of the i-th image. The first
%   4 columns of this array contain the bounding box coordinates in the
%   form of [x0,y0,x1,y1] (where the (x0,y0) and (x1,y1) are the top-left
%   and bottom-right corners) and the rest C columns contain the confidence
%   score of the bounding boxes for each of the C categories.
%   b) (the flag is_per_class must being set to true) a C x 1 cell array
%   where C is the number of categories. The aboxes_in{j} element in this case
%   is NI x 1 cell array, where NI is the number of images, with the candidate
%   bounding box detections of the j-th category for each image. The
%   element aboxes_in{j}{i} is a NB_{i,j} x 5 array, where NB_{i,j} is the
%   number candidate bounding box detections of the i-th image for the j-th
%   category. The first 4 columns of this array are the bounding box
%   coordinates in the form of [x0,y0,x1,y1] (where the (x0,y0) and (x1,y1)
%   are the top-left and bottom-right corners) and the 5-th column contains
%   the confidence score of each bounding box with respect to the j-th category.
% 2) The rest input arguments are given in the form of Name,Value pair
% arguments and are:
% 'threshold': is a C x 1 array, where C is the number of categories.
% It must contains the threshold per category that will be used for
% removing candidate boxes with low confidence prior to applying the
% non-max-suppression step.
% 'nms_iou_thrs': scalar value with the IoU threshold that will be used
% during the non-max-suppression step.
% 'is_per_class': boolean value that if set to false, then the 1.a) form
% of the aboxes_in input parameter will be expected; otherwise, if set to
% true then the 1.b) form of the aboxes_in input parameter will be expected
% 'max_per_image': scalar value with the maximum number of detection per
% image and per category. Default value: 200.
% 'do_bbox_voting': boolean value that if is set to True then the box voting
% step is applied. Default value: false
% 'box_ave_iou_thresh': scalar value with the minimum IoU threshold that
% is used in order to define the neighbors of bounding box during the box
% voting step. Default value: 0.5
% 'add_val': scalar value that is added to the confidence score of bounding
% boxes in order to compute the box weight during the box voting step.
% Default value: 1.5
% 'use_not_static_thr': boolean value that if set to true, then the threshold
% (param threshold) that is applied on the candidate bounding box
% detections will be modified (decreased from its input value) such that
% the average number of detections per image and per category to be
% ave_per_image. Default value: true
% 'ave_per_image': scalar value with the desired average number of
% detections per image and per category. Used only when the use_not_static_thr
% input parameter is set to true. default value: 5
%
% OUTPUT:
% 1) aboxes_out: is a C x 1 array where C is the number of
% categories. The aboxes_out{j} element is a NI x 1 cell array, where NI is
% the number of images, with the bounding box detections of the j-th
% category for each image. The element aboxes_out{j}{i} is a ND_{i,j} x 5 array,
% where ND_{i,j} is the number bounding box detections for the i-th image
% and the j-th category. The first 4 columns of this array contain the
% bounding box coordinates in the form of [x0,y0,x1,y1] (where the (x0,y0)
% and (x1,y1) are the top-left and bottom-right corners) and the 5-th
% column contains the confidence score of each bounding box with respect to
% the j-th category.
%
% This file is part of the code that implements the following ICCV2015 accepted paper:
% title: "Object detection via a multi-region & semantic segmentation-aware CNN model"
% authors: Spyros Gidaris, Nikos Komodakis
% institution: Universite Paris Est, Ecole des Ponts ParisTech
% Technical report: http://arxiv.org/abs/1505.01749
% code: https://github.com/gidariss/mrcnn-object-detection
%
% Part of the code in this file comes from the R-CNN code:
% https://github.com/rbgirshick/rcnn
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2015 Spyros Gidaris
%
% "Object detection via a multi-region & semantic segmentation-aware CNN model"
% Technical report: http://arxiv.org/abs/1505.01749
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

ip = inputParser;
ip.addParamValue('is_per_class',      false, @islogical);
ip.addParamValue('do_bbox_voting',    false, @islogical);
ip.addParamValue('box_ave_iou_thresh',  0.5, @isnumeric);
ip.addParamValue('use_not_static_thr', true, @islogical);
ip.addParamValue('ave_per_image',         5, @isnumeric);
ip.addParamValue('max_per_image',       200, @isnumeric);
ip.addParamValue('add_val',             1.5, @isnumeric);
ip.addParamValue('use_gpu',           false, @islogical);
ip.addParamValue('nms_iou_thrs',        0.3, @isnumeric);
ip.addParamValue('threshold',            -4, @isnumeric);

ip.parse(varargin{:});
opts = ip.Results;

nms_iou_thrs       = opts.nms_iou_thrs;
threshold         = opts.threshold;
is_per_class       = opts.is_per_class;
do_bbox_voting     = opts.do_bbox_voting;
box_ave_iou_thresh = opts.box_ave_iou_thresh;
use_not_static_thr = opts.use_not_static_thr;
max_per_image      = opts.max_per_image;
ave_per_image      = opts.ave_per_image;
add_val            = opts.add_val;
use_gpu            = opts.use_gpu;

if ~is_per_class
    num_imgs    = length(aboxes_in);
    num_classes = size(aboxes_in{1},2) - 4;
else
    num_classes = length(aboxes_in);
    num_imgs    = length(aboxes_in{1});
end

max_per_set   = ceil(ave_per_image * num_imgs);
aboxes_out    = cell(num_classes, 1);
for i = 1:num_classes, aboxes_out{i} = cell(num_imgs, 1); end

thresh_val_all = ones(num_classes,1,'single') * threshold;

for i = 1:num_imgs
    % get the candidate bounding box detection of one image
    if ~is_per_class
        assert(size(aboxes_in{i},2) == (4 + num_classes));
        bbox_cand_dets = aboxes_in{i};
    else
        bbox_cand_dets = cell(num_classes,1);
        for j = 1:num_classes
            assert(size(aboxes_in{j}{i},2) == 5);
            bbox_cand_dets{j} = aboxes_in{j}{i};
        end
    end

    % post-process the candidate bounding box detection of one image
    bbox_detections_per_class = post_process_candidate_detections(bbox_cand_dets, ...
        'thresholds',thresh_val_all, 'nms_iou_thrs',nms_iou_thrs,...
        'max_per_image',max_per_image,'do_bbox_voting',do_bbox_voting,...
        'box_ave_iou_thresh',box_ave_iou_thresh,'add_val',add_val, ...
        'use_gpu',use_gpu);

    for j = 1:num_classes, aboxes_out{j}{i} = bbox_detections_per_class{j}; end

    if (mod(i, 1000) == 0) && use_not_static_thr
        for j = 1:num_classes
            [aboxes_out{j}, thresh_val_all(j)] = keep_top_k(aboxes_out{j}, i, max_per_set, thresh_val_all(j));
        end
    end
end

if use_not_static_thr
    disp(thresh_val_all(:)');
    for i = 1:num_classes
        % go back through and prune out detections below the found threshold
        aboxes_out{i} = prune_detections(aboxes_out{i}, thresh_val_all(i));
    end
end

total_num_bboxes = zeros(num_classes, 1);
for j = 1:num_classes
    total_num_bboxes(j) = sum(cellfun(@(x) size(x,1), aboxes_out{j}, 'UniformOutput', true));
end

fprintf('Average number of bounding boxes per class\n');
disp(total_num_bboxes(:)' / num_imgs);
fprintf('Average number of bounding boxes in total\n');
disp(sum(total_num_bboxes) / num_imgs);
end

function [boxes, thresh] = keep_top_k(boxes, end_at, top_k, thresh)
% keep top K detections
X = cat(1, boxes{1:end_at});
if isempty(X), return; end

scores = sort(X(:,end), 'descend');
thresh = scores(min(length(scores), top_k));
for image_index = 1:end_at
    bbox = boxes{image_index};
    keep = find(bbox(:,end) >= thresh);
    boxes{image_index} = bbox(keep,:);
end

end

function bbox_dets = prune_detections(bbox_dets, thresh)
for j = 1:length(bbox_dets)
    if ~isempty(bbox_dets{j})
        bbox_dets{j}(bbox_dets{j}(:,end) < thresh ,:) = [];
    end
end
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function[h]=hlines(y,sty)
%HLINES   Add horizontal lines to a plot.
%
%   HLINES(Y,STYLE) puts horizontal lines at locations Y having style
%   described by the string STYLE.  HLINES draws the lines to fit the
%   current axis.
%
%   The SYTLE strings follow the format specified in LINESTYLE.  Thus
%   STYLE='2b--' draws blue dotted lines of width 2.
%
%   HLINES(Y) uses the default STYLE='g--', a green dashed line.
%
%   H=HLINES(...) returns an array of handles H to the lines.
%
%   See also VLINES, DLINES.
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2000--2008 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(y,'--t')
    return
end

if nargin==1
	sty='g--';
end

hold on
y=y(:)';
ax=axis;
x1=ax(1)*ones(size(y));
x2=ax(2)*ones(size(y));
y=[y;y];
x=[x1;x2];
h=plot(x,y);
linestyle(h,sty);
axis(ax)

if nargout==0
   clear h
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function value = r8vec_distance ( dim_num, v1, v2 )

%*****************************************************************************80
%
%% R8VEC_DISTANCE returns the Euclidean distance between two R8VEC's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, real V1(DIM_NUM), V2(DIM_NUM), the vectors.
%
%    Output, real VALUE, the Euclidean distance
%    between the vectors.
%
  value = sqrt ( sum ( ( v1(1:dim_num) - v2(1:dim_num) ).^2 ) );

  return
end


text:
mask = 5*double(instances(:,:,1))+5*double(labels(:,:,1));
maskColor = zeros(size(image,1)*size(image,2),3);
uniquemask = unique(mask(:));
for i =1:length(uniquemask)

    sel = mask(:)==uniquemask(i);

    maskColor(sel, :) = repmat(ObjectColor(uniquemask(i)),sum(sel),1);

end
maskColor(find(mask(:)==0),:) = 1;
maskColor = reshape(maskColor,[size(mask,1),size(mask,2),3]);

figure

imshow(maskColor);
imwrit(maskColor,'NYUmask.png')
%%
load('/n/fs/modelnet/SUN3DV2/prepareGT/cls.mat')

addpath('/n/fs/modelnet/SUN3DV2/roomlayout/')
fullname = '/n/fs/sun3d/data/rgbd_voc/000414_2014-06-04_19-49-13_260595134347_rgbf000044-resize'
data = readframe(fullname);
groundTruthBbs  = data.groundtruth3DBB;

sequenceName = getSequenceName(fullname);
gtRoom3D = GroundTruthBox(sequenceName,0);
cameraXYZ = data.anno_extrinsics'*gtRoom3D;
cameraXYZ([2 3],:) = cameraXYZ([3 2],:);
cameraXYZ(3,:) = - cameraXYZ(3,:);
cameraXYZ = data.Rtilt * cameraXYZ;

my_mhCorner3D = cameraXYZ; %data.Rtilt*data.anno_extrinsics'*gtCorner3D;
visulize_wholeroom(groundTruthBbs,cls,fullname,my_mhCorner3D)

author: thusiyuan


text:
function kmph = ftps2kmph(ftps)
%FTPS2KMPH Convert speed from feet per second to kilometers per hour
%
%  kmph = FTPS2KMPH(ftps) converts speeds from feet per second to
%   kilometers per hour.
%
%  See also FTPS2KTS, FTPS2MPH, FTPS2MPS, KMPH2FTPS.

% Jonathan Sullivan
% Original: May 2011
% jonathan.sullivan@ll.mit.edu

kmph = ftps*1.09728;


text:
classdef GFMMOEA < ALGORITHM
% <multi/many> <real/integer/label/binary/permutation>
% Generic front modeling based multi-objective evolutionary algorithm
% theta --- 0.2 --- Penalty parameter
% fPFE  --- 0.1 --- Frequency of employing generic front modeling

%------------------------------- Reference --------------------------------
% Y. Tian, X. Zhang, R. Cheng, C. He, and Y. Jin, Guiding evolutionary
% multi-objective optimization with generic front modeling, IEEE
% Transactions on Cybernetics, 2020, 50(3): 1106-1119.
%------------------------------- 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
            [theta,fPFE] = Algorithm.ParameterSet(0.2,0.1);

            %% Generate random population
            Population = Problem.Initialization();
            FrontNo    = NDSort(Population.objs,inf);
            zmin       = min(Population.objs,[],1);
            % Calculate the fitness of each solution
            [P,A]     = deal(ones(1,Problem.M));
            [App,Dis] = CalFitness(Population.objs-repmat(zmin,length(Population),1),P,A);
            Dis       = sort(Dis,2);
            Crowd     = Dis(:,1) + 0.1*Dis(:,2);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = TournamentSelection(2,Problem.N,FrontNo,-theta*App-(1-theta)*Crowd);
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                zmin       = min([zmin;Offspring.objs],[],1);
                if ~mod(ceil(Problem.FE/Problem.N),ceil(fPFE*ceil(Problem.maxFE/Problem.N))) || fPFE == 0
                    [P,A] = GFM(Population(FrontNo==1).objs-repmat(zmin,sum(FrontNo==1),1));
                end
                [Population,FrontNo,App,Crowd] = EnvironmentalSelection([Population,Offspring],P,A,zmin,theta,Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%CODEGENERATOR.GENSLBLOCKCORIOLIS Generate Simulink block for Coriolis matrix
%
% cGen.genslblockcoriolis() generates a robot-specific Simulink block to compute
% Coriolis/centripetal matrix.
%
% Notes::
% - Is called by CodeGenerator.gencoriolis if cGen has active flag genslblock
% - The Coriolis matrix is stored row by row to avoid memory issues.
% - The Simulink block recombines the the individual blocks for each row.
% - Access to generated function is provided via subclass of SerialLink
%   whose class definition is stored in cGen.robjpath.
%
% Author::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.CodeGenerator, CodeGenerator.gencoriolis.

% Copyright (C) 2012-2014, by Joern Malzahn
%
% 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
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.
function [ ] = genslblockcoriolis( CGen )

%% Open or create block library
bdclose('all')                                                              % avoid problems with previously loaded libraries
load_system('simulink');
if ~(exist([CGen.slibpath,simulinkext]) == 2)                                  % Create new block library if none exists
 CGen.createnewblocklibrary;
end
open_system(CGen.slibpath);
set_param(CGen.slib,'lock','off');

[q,qd] = CGen.rob.gencoords;

%% Generate Coriolis Block
CGen.logmsg([datestr(now),'\tGenerating Simulink Block for the robot Coriolis matrix\n']);
nJoints = CGen.rob.n;
symname = 'coriolis';

CGen.logmsg([datestr(now),'\t\t... enclosing subsystem ']);
CoriolisBlock = [CGen.slib,'/',symname];
if ~isempty(find_system(CGen.slib,'SearchDepth',1,'Name',symname))                    % Delete previously generated block
    delete_block(CoriolisBlock)
    save_system;
end

% Subsystem in which individual rows are concatenated
add_block('built-in/SubSystem',CoriolisBlock);                               % Add new inertia matrix block
add_block('simulink/Math Operations/Matrix Concatenate'...
    , [CoriolisBlock,'/coriolis']...
    , 'NumInputs',num2str(nJoints)...
    , 'ConcatenateDimension','1');
add_block('simulink/Sinks/Out1',[CoriolisBlock,'/out']);
add_block('simulink/Sources/In1',[CoriolisBlock,'/q']);
add_block('simulink/Sources/In1',[CoriolisBlock,'/qd']);
add_line(CoriolisBlock,'coriolis/1','out/1');
CGen.logmsg('\t%s\n',' done!');

for kJoints = 1:nJoints
    CGen.logmsg([datestr(now),'\t\t... Embedded Matlab Function Block for joint ',num2str(kJoints),': ']);

    % Generate Embedded Matlab Function block for each row
    symname = ['coriolis_row_',num2str(kJoints)];
    fname = fullfile(CGen.sympath,[symname,'.mat']);

    if exist(fname,'file')
        tmpStruct = load(fname);
    else
        error ('genslblockcoriolis:SymbolicsNotFound','Save symbolic expressions to disk first!')
    end

    blockaddress = [CoriolisBlock,'/',symname];
    if doesblockexist(CGen.slib,symname)
        delete_block(blockaddress);
        save_system;
    end

    CGen.logmsg('%s',' block creation');
    symexpr2slblock(blockaddress,tmpStruct.(symname),'vars',{q,qd});

    % connect output
    CGen.logmsg('%s',', output wiring');
    if ( verLessThan('matlab','7.11.0.584') ) && ( isequal(tmpStruct.(symname),zeros(1,nJoints)) )
        % There is a bug in earlier Matlab versions. If the symbolic
        % vector is a zero vector, then the Simulink Embedded Matlab
        % Function block outputs a scalar zero. We need to concatenate
        % a row vector of zeros here, which we have to construct on our
        % own.
        add_block('simulink/Math Operations/Matrix Concatenate'...              % Use a matrix concatenation block ...
            , [CoriolisBlock,'/DimCorrection',num2str(kJoints)]...               % ... named with the current row number ...
            , 'NumInputs',num2str(nJoints),'ConcatenateDimension','2');         % ... intended to concatenate zero values for each joint ...
        % ... columnwise. This will circumvent the bug.

        for iJoints = 1:nJoints                                                 % Connect signal lines from the created block (which outputs
            add_line(CoriolisBlock...                                            % a scalar zero in this case) with the bugfix block.
                , [symname,'/1']...
                , ['DimCorrection',num2str(kJoints),'/', num2str(iJoints)]);
        end

        add_line(CoriolisBlock,['DimCorrection',num2str(kJoints)...              % Connect the fixed row with other rows.
            , '/1'],['coriolis/', num2str(kJoints)]);

    else
        add_line(CoriolisBlock,[symname,'/1']...          % In case that no bug occurs, we can just connect the rows.
            , ['coriolis/', num2str(kJoints)]);
    end

    % Vector generalized joint values
    add_line(CoriolisBlock,'q/1',[symname,'/1']);
    % Vector generalized joint velocities
    add_line(CoriolisBlock,'qd/1',[symname,'/2']);

    CGen.logmsg('\t%s\n','row complete!');
end
addterms(CoriolisBlock); % Add terminators where needed
distributeblocks(CoriolisBlock);
CGen.logmsg([datestr(now),'\tCoriolis matrix block complete\n']);

%% Cleanup
% Arrange blocks
distributeblocks(CGen.slib);

% Lock, save and close library
set_param(CGen.slib,'lock','on');
save_system(CGen.slib,CGen.slibpath);
close_system(CGen.slib);

end

author: petercorke

repo: robotics-toolbox-matlab


text:
function obj = t1_pnp_func_new(in1,in2,in3)
coefs_tq1_1 = in2(1);
coefs_tq1_2 = in2(4);
coefs_tq2_1 = in2(2);
coefs_tq1_3 = in2(7);
coefs_tq2_2 = in2(5);
coefs_tq3_1 = in2(3);
coefs_tq1_4 = in2(10);
coefs_tq2_3 = in2(8);
coefs_tq3_2 = in2(6);
coefs_tq1_5 = in2(13);
coefs_tq2_4 = in2(11);
coefs_tq3_3 = in2(9);
coefs_tq1_6 = in2(16);
coefs_tq2_5 = in2(14);
coefs_tq3_4 = in2(12);
coefs_tq1_7 = in2(19);
coefs_tq2_6 = in2(17);
coefs_tq3_5 = in2(15);
coefs_tq1_8 = in2(22);
coefs_tq2_7 = in2(20);
coefs_tq3_6 = in2(18);
coefs_tq1_9 = in2(25);
coefs_tq2_8 = in2(23);
coefs_tq3_7 = in2(21);
coefs_tq2_9 = in2(26);
coefs_tq3_8 = in2(24);
coefs_tq3_9 = in2(27);
coefs_tq1_10 = in2(28);
coefs_tq2_10 = in2(29);
coefs_tq3_10 = in2(30);
pinvG1_1 = in1(1);
pinvG1_2 = in1(4);
pinvG1_3 = in1(7);
q0 = in3(1,:);
q1 = in3(2,:);
q2 = in3(3,:);
q3 = in3(4,:);
obj = q0.^2.*(coefs_tq1_1.*pinvG1_1+coefs_tq2_1.*pinvG1_2+coefs_tq3_1.*pinvG1_3)+q1.^2.*(coefs_tq1_5.*pinvG1_1+coefs_tq2_5.*pinvG1_2+coefs_tq3_5.*pinvG1_3)+q2.^2.*(coefs_tq1_8.*pinvG1_1+coefs_tq2_8.*pinvG1_2+coefs_tq3_8.*pinvG1_3)+q3.^2.*(coefs_tq1_10.*pinvG1_1+coefs_tq2_10.*pinvG1_2+coefs_tq3_10.*pinvG1_3)+q0.*q1.*(coefs_tq1_2.*pinvG1_1+coefs_tq2_2.*pinvG1_2+coefs_tq3_2.*pinvG1_3)+q0.*q2.*(coefs_tq1_3.*pinvG1_1+coefs_tq2_3.*pinvG1_2+coefs_tq3_3.*pinvG1_3)+q0.*q3.*(coefs_tq1_4.*pinvG1_1+coefs_tq2_4.*pinvG1_2+coefs_tq3_4.*pinvG1_3)+q1.*q2.*(coefs_tq1_6.*pinvG1_1+coefs_tq2_6.*pinvG1_2+coefs_tq3_6.*pinvG1_3)+q1.*q3.*(coefs_tq1_7.*pinvG1_1+coefs_tq2_7.*pinvG1_2+coefs_tq3_7.*pinvG1_3)+q2.*q3.*(coefs_tq1_9.*pinvG1_1+coefs_tq2_9.*pinvG1_2+coefs_tq3_9.*pinvG1_3);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function varargout = gui_loadset(varargin)
% load a dataset via the GUI
% [Dataset,Action] = gui_loadset()
% [Dataset,Action] = gui_loadset('Initialpath',Initialpath)
%
% Out:
%   Dataset     : the dataset, an expression (or [], if the user pressed Cancel)
%   Action      : either 'OK' or 'Cancel'
%
%                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                           2010-10-27

% Last Modified by GUIDE v2.5 05-Jan-2012 07:14:18

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

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

% --- Executes just before gui_loadset is made visible.
function gui_loadset_OpeningFcn(hObject, eventdata, handles, varargin)
if length(varargin) < 2
    initialpath = env_translatepath('data:/');
else
    initialpath = env_translatepath(varargin{2});
end

% run file dialog...
[handles.filename,handles.filepath] = uigetfile({'*.*', 'any supported file'; '*.set', 'EEGLAB dataset'; '*.vhdr', 'BrainProducts Vision Recorder file'; ...
    '*.bdf', 'BioSEMI BDF file'; '*.cnt', 'Neuroscan CNT / ANT EEProbe CNT file'; '*.avr', 'ANT EEProbe average / Megis/BESA file'; '*.mat', 'MATLAB file'; ...
    '*.raw', 'Electrical Geodesics / ERPSS / Yokogawa raw data file '; '*.sna', 'Snapmaster SNA file'; '*.eeg', 'Neuroscan EEG / ANT EEProbe / BrainVision data file'; ...
    '*.edf', 'European data format file'; '*.ds', 'CTF folder'; '*.rdf', 'ERPSS data file'; '*.asc', 'INStep ASC file'; '*.m4d', '4D pdf file'; ...
    '*.smr', 'Cambridge Electronic Design file'; '*.egis', 'Electrical Geodesics file'; '*.ave', 'Yokogawa / Electrical Geodesics file'; ...
    '*.gave', 'Electrical Geodesics file'; '*.ses', 'Electrical Geodesics file '; '*.swf', 'Megis/BESA file'; '*.avg', 'NeuroScan file'; '*.nxe', 'Nexstim file'; ...
    '*.dat', 'BrainVision data file'; '*.seg', 'BrainVision data file'; '*.res4', 'VSM MedTech file'; '*.meg4', 'VSM MedTech file'; ...
    '*.m4d', 'Neuromag Elektra / 4D Neuroimaging file'; '*.pdf', 'Neuromag Elektra / 4D Neuroimaging file'; '*.xyz', 'Neuromag Elektra / 4D Neuroimaging file'},'Select dataset(s) to load',initialpath);

% Choose default command line output for gui_loadset
handles.output = {[],'Cancel'};
% Update handles structure
guidata(hObject, handles);

if ~isnumeric(handles.filename)
    % UIWAIT makes gui_loadset wait for user response (see UIRESUME)
    uiwait(handles.figure1);
else
    % otherwise we directly close the figure again...
    close(handles.figure1);
end

% --- Outputs from this function are returned to the command line.
function varargout = gui_loadset_OutputFcn(hObject, eventdata, handles)
try
    varargout = handles.output;
    close(handles.figure1);
catch
    varargout = {[],'Cancel'};
end

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% construct optional arguments
args = {};
channels = str2num(get(handles.edit2,'String'));
if ~isempty(channels)
    args = [args {'channels' channels}]; end
timerange = str2num(get(handles.edit5,'String'));
if ~isempty(timerange)
    args = [args {'timerange' timerange}]; end
samplerange = str2num(get(handles.edit4,'String'));
if ~isempty(samplerange)
    args = [args {'samplerange' samplerange}]; end
types = get(handles.edit8,'String');
if ~isempty(types)
    % neither quoted nor a cell array --> put in quotes
    if types(1) ~= '''' && types(1) ~= '{'
        types = ['''' types '''']; end
    % evaluate & make a cell if necessary
    types = eval(types);
    if ~iscell(types)
        types = {types}; end
    args = [args {'types' types}];
end
miscopts = get(handles.edit6,'String');
if ~isempty(miscopts)
    args = [args evalin('base',['{' miscopts '}'])]; end
if isfield(handles,'markerchannel')
    args = [args {'markerchannel', {handles.markerchannel}}]; end
if ~isfield(handles,'insert_params')
    handles.insert_params = {}; end

% construct the loader expression
if iscell(handles.filename)
    for c=1:length(handles.filename)
        dataset{c} = io_loadset([handles.filepath,handles.filename{c}],args{:});
        for p=1:length(handles.insert_params)
            dataset{c} = set_insert_markers(dataset{c},handles.insert_params{p}); end

    end
    dataset = set_concat(dataset{:});
else
    dataset = io_loadset([handles.filepath,handles.filename],args{:});
    for p=1:length(handles.insert_params)
        dataset = set_insert_markers(dataset,handles.insert_params{p}); end
end

% optionally check if the dataset is loadable
if get(handles.checkbox1,'Value') == get(handles.checkbox1,'Max')
    tmp = exp_eval(dataset); end

% assign in workspace
assignin('base',get(handles.edit1,'String'),dataset);
handles.output = {dataset,'OK'};
guidata(hObject,handles);
uiresume(handles.figure1);

function pushbutton2_Callback(hObject, eventdata, handles)
handles.output = {[],'Cancel'};
guidata(hObject,handles);
uiresume(handles.figure1);

function pushbutton3_Callback(hObject, eventdata, handles)
handles.insert_params = gui_eventinsertion();
guidata(hObject,handles);

function pushbutton6_Callback(hObject, eventdata, handles)
tmp = arg_guidialog(@set_infer_markers,'Invoke',false);
if ~isempty(tmp)
    handles.markerchannel = tmp; end
guidata(hObject,handles);

function edit1_Callback(hObject, eventdata, handles)
gui_check_varname(hObject,'lastdata');

function edit2_Callback(hObject, eventdata, handles)
gui_check_vector(hObject);

function edit3_Callback(hObject, eventdata, handles)
gui_check_scalar(hObject);

function edit4_Callback(hObject, eventdata, handles)
gui_check_range(hObject);

function edit5_Callback(hObject, eventdata, handles)
gui_check_range(hObject);

function edit8_Callback(hObject, eventdata, handles)

function pushbutton5_Callback(hObject, eventdata, handles)
% help button
env_doc io_loadset

function figure1_KeyPressFcn(hObject, eventdata, handles)
if strcmp(eventdata.Key,'return')
    pushbutton1_Callback(handles.pushbutton1, eventdata, handles); end
if strcmp(eventdata.Key,'escape')
    pushbutton2_Callback(handles.pushbutton2, eventdata, handles); end

% --- auto-generated trash ---

% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
function edit2_CreateFcn(hObject, eventdata, handles)
function edit3_CreateFcn(hObject, eventdata, handles)
function edit4_CreateFcn(hObject, eventdata, handles)
function edit5_CreateFcn(hObject, eventdata, handles)
function edit6_Callback(hObject, eventdata, handles)
function edit6_CreateFcn(hObject, eventdata, handles)
function checkbox1_Callback(hObject, eventdata, handles)
function edit8_CreateFcn(hObject, eventdata, handles)
function figure1_CreateFcn(hObject, eventdata, handles)

repo: REST


text:
% LORENZ system
% Evaluate ensemble of models (e.g. for varying training data length, noise
% level, etc.) in MPC

clear all, close all, clc
addpath('../utils');

WORKING = 0;
SystemModel = 'LORENZ';

%% Parameters
InputSignalType = 'sphs';
ModelName = 'NARX'; % DMDc, NARX, SINDYc
TrainAlg = 'trainbr'; %trainlm, trainbr % Choose: trainbr: eta_vec ~= 0, trainlm: eta_vec = 0
select_model = ModelName;
dep_trainlength = 1;
dep_noise = 0;

% Dependency on training length, noise-free training data
% Ntrain_vec = [5:15,20:5:95,100:100:1000];%,1500:500:3000];
% eta_vec = 0;
% Nr = 1;

% Dependency on training length, fixed noise level of data
Ntrain_vec = [5:15,20:5:95,100:100:1000,1250,1500,2000,3000];
eta_vec = 0.05;
Nr = 1;

% Some parameters
N_LENGTHS = length(Ntrain_vec);
Nmodels = N_LENGTHS;

%% Folders
figpath0 = ['../FIGURES/EX_',SystemModel,'_Dependencies/'];
datapath0 = ['../DATA/EX_',SystemModel,'_Dependencies/'];

if strcmp(ModelName,'NARX')
    figpath = [figpath0, ModelName, '/',TrainAlg,'/TL/'];
    datapath = [datapath0, ModelName, '/',TrainAlg,'/TL/'];
else
    figpath = [figpath0, ModelName, '/TL/'];
    datapath = [datapath0, ModelName, '/TL/'];
end

mkdir(figpath)
mkdir(datapath)
%% Load all models
switch ModelName
    case 'DMDc'
        Models(1:Nmodels) = struct('name',[],'sys',[],'Ndelay', [], 'Ttraining', [],'dt',[], 'Err', [], 'ErrM', []);
    case 'SINDYc'
        Models(1:Nmodels) = struct('name',[],'polyorder',[],'usesine',[],'Xi',[],'dt',[],'N',[], 'Ttraining', [], 'Err', [], 'ErrM', []);
    case 'NARX'
        Models(1:Nmodels) = struct('name',[],'net',[],'stateDelays',[],'inputDelays',[], 'hiddenSizes', [], 'dt',[], 'Ttraining', [], 'Err', [], 'ErrM', []);
end

if eta_vec == 0
    for iModel = 1:Nmodels
        iN = iModel;
        load(fullfile(datapath,['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)), ...
            '_Eta',sprintf('%03g',100*eta_vec),'_iR1','.mat']))
        Models(iModel) = Model;
    end
elseif eta_vec~=0
    for iModel = 1:Nmodels
        iN = iModel;
        load(fullfile(datapath,['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)), ...
            '_Eta',sprintf('%03g',100*eta_vec),'_BEST_MODEL.mat']))
        Models(iModel) = Model;
    end
end

%% Run MPC for all models
runMPC

%% Save results
save(fullfile(datapath,['EX_',SystemModel,'_MPC_',ModelName,'_',InputSignalType,'_TrainLength','.mat']),'Results')

%% Show performance
figure,hold on
for iM = 1:Nmodels
    Jend = cumsum(Results(iM).J);
    plot(Ntrain_vec(iM),Jend(end),'ok')
end
set(gca,'xscale','log','yscale','log')

author: eurika-kaiser

repo: SINDY-MPC


text:
% NOLM_Bifurcation_Diagram.
% Copyright Springer 2013 A.L. Steele and S. Lynch.
clear
N=49999;halfN=N/2;
lambda=1.55E-6;n2=3.2E-20;Aeff=30E-12;L=80;
E1(1)=0;Pmax=40;phi=0;Gauss(1)=0;
Pin(1:N)=0; Pout(1:N)=0;
phif=0*pi; E6p=0;
kappa1=0.25;kappa2=0.8;kappa3=0.8;
rootk1 = sqrt(kappa1); irootk1 = 1i*sqrt(1-kappa1);
rootk2 = sqrt(kappa2); irootk2 = 1i*sqrt(1-kappa2);
rootk3 = sqrt(kappa3); irootk3 = 1i*sqrt(1-kappa3);
G=sqrt((1-kappa2)*(1-kappa3));
% Ramp the power up and down
for n=1:N

    Ein = sqrt(Pmax*exp(-0.02*((n*Pmax/N-Pmax/2))^2));
    E1 = rootk3*Ein+irootk3*E6p;
    P1 = abs(E1)^2;

    E3 = rootk1*E1;
    E4 = irootk1*E1;

    phic=2*pi*n2*L*(2-kappa1)*P1/(lambda*Aeff);
    phicc=2*pi*n2*L*(1+kappa1)*P1/(lambda*Aeff);

    E3p = E3*exp(-1i*(phi+phic));
    E4p = E4*exp(-1i*(phi+phicc));

    E5 = rootk1*E3p+irootk1*E4p;

    Pout(n) = abs(rootk2*E5)^2;
    Pin(n) = abs(Ein)^2;
    E6 = irootk2*E5;

    E6p = E6*exp(1i*phif);

end

% Plot the bifurcation diagrams
figure(1)
clf
fsize=15;
subplot(2,1,1)
hold on
plot(Pout(1:N),'.','MarkerSize',1)
plot(Pin(1:N),'.','MarkerSize',1);
xlabel('Number of Ring Passes','FontSize',fsize);
ylabel('Output Power','FontSize',fsize);
hold off

subplot(2,1,2)
hold on
plot(Pin,Pout,'.','MarkerSize',1);
xlabel('Input Power','FontSize',fsize);
ylabel('Output Power','FontSize',fsize);
hold off


text:
 function [x, mask0] = mri_phase_denoise(yi, varargin)
%function [x, mask0] = mri_phase_denoise(yi, [options]) <- recommended usage
%|function [x] = mri_phase_denoise(yi, l2b, niter, chat, wthresh) <- old way
%|
%| in
%|	yi	[(N)]	noisy complex image: y = mag .* exp(1i * x)
%|			can be any size: 1D, 2D, 3D, ...
%|
%| options
%|	l2b		log_2(beta), regularization parameter
%|	order		regularization order (default: 1, for historical
%|				reasons, but 2 is probably preferable)
%|	niter		# of iterations
%|	chat		1 to show pictures
%|	wthresh		fraction of magnitude maximum to include in fitting
%|	init		initial image for iterations
%|	isave		which iterations to save.  (default: last)
%|	'pl'	1|0	1 for new PL method (recommend), 0 for old way (default)
%| out
%|	x	[(N)]	cleaned up phase estimate (radians)
%|	mask0	[(N)]	logical: 1 for high-magnitude pixels
%|
%| Example of weighted "denoising" of MRI phase images.
%| This is a "simple" way to estimate good field inhomogeneity maps
%| from the usual approach of two readouts with a short delay.
%| It also smoothly interpolates over regions with signal voids.
%|
%| Caution: the sign of the field map estimated here is the opposite (negative)
%| of the sign of the field map needed for input to the Gmri object.
%|
%| Copyright 1999, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if ischar(yi)
	[x, mask0] = mri_phase_denoise_test(yi, varargin{:});
	if ~nargout, clear x mask0, end
return
end

% defaults
arg.chat = 0;
arg.l2b = -5;
arg.order = 1;
arg.niter = 150;
arg.isave = [];
arg.wthresh = 0.4;
arg.fmax = 0.05; % fraction of max threshold for trimmed median for wi_ml
arg.init = [];
arg.clim = []; % limits for phase display
arg.pl = false; % PL
arg.wi_ml = false; % use wi based on ML instead of threshold

% backword compatible for old argument list: l2b, niter, chat, wthresh
if length(varargin) && isnumeric(varargin{1})
	arg.l2b = varargin{1};
	if length(varargin) >= 2, arg.niter = varargin{2}; end
	if length(varargin) >= 3, arg.chat = varargin{3}; end
	if length(varargin) >= 4, arg.wthresh = varargin{4}; end
else
	arg = vararg_pair(arg, varargin);
end
if isempty(arg.isave), arg.isave = arg.niter; end

mag = abs(yi);
yi = angle(yi);
if arg.chat
	im plc 2 3
	im(1, mag, 'magnitude'), cbar
	im(2, yi, 'raw phase map', arg.clim), cbar
end

dim_yi = size(yi);

%
% specify weights: this needs more work to be automatic!
%
mask0 = mag > arg.wthresh * max(mag(:)); % ignore pixels with "too small" magnitude
if arg.chat
	im(3, mask0, 'weights'), cbar
	im(4, mask0 .* yi, 'masked phase', arg.clim), cbar
end

%mean(mag(:))
%median(mag(:))
%clf, hist(mag(:), 100), pause
%median(mag(mag(:) > 0.05 * max(mag(:))))

if arg.wi_ml
	wi = mri_phase_wi_ml(mag, arg.fmax);
else
	wi = mask0;
end
%W = diag_sp(wi(:));
W = Gdiag(wi);

%
% initial phase image
%
if isempty(arg.init)
	arg.init = yi;
	arg.init(mask0 == 0) = mean(yi(mask0 == 0));
end
if arg.chat
	im(5, arg.init, 'Initial phase', arg.clim), cbar
end

%G = diag_sp(ones(prod(dim_yi),1));
G = Gdiag(ones([dim_yi 1]));

%
% regularizer
%
if arg.order ~= 2, warn('order=2 recommended'), end
mask1 = true(size(yi)); % estimate / extrapolate to *all* pixels
R = Reg1(mask1, 'beta', 2^arg.l2b, 'order', arg.order);
%R = Robject(mask1, 'beta', 2^arg.l2b, 'order', arg.order);
%	'type_denom', 'matlab', ...

if 0 % old way
	[C, wjk] = C2sparse('tight', mask1, 8); % todo: cut?  do not use!
	C = spdiag(sqrt(wjk), 'nowarn') * C; % caution: missing prior to 2005-11-28
	C = sqrt(2^arg.l2b) * C;
end

% report expected blur (at image center)
if 1
	qpwls_psf(G, R, 1, mask1, W);
end

%
% run qpwls algorithm for regularized fitting
%
xinit = arg.init(mask1);
if arg.pl
	med = median(mag(mag(:) > 0.05 * max(mag(:))));
	data = {yi(:), (mag(:)/med).^2}; handle = @phase_dercurv; % PL
%profile on % todo!
	x = pl_pcg_qs_ls(xinit, G, data, handle, R, ...
		'niter', arg.niter, 'isave', arg.isave);
%profile report
	x = embed(x, mask1);
else
	warn 'recommend using "pl" option.  use wls for historical only'
%	data = {yi(:), wi(:)}; handle = @wls_dercurv; % qpwls
	x = qpwls_pcg1(xinit, G, W, yi(:), R.C, ...
		'niter', arg.niter, 'isave', arg.isave);
	x = embed(x, mask1);
end

if 0 % old way
	x = qpwls_pcg(x, G, W, yi(:), 0, R.C, 1, arg.niter);
	x = reshape(x, [dim_yi arg.niter]);
	x = x(:,:,end);
end

if arg.chat
	if ndims(yi) == 2
		im(6, x(:,:,end), 'QPWLS-CG phase', arg.clim), cbar
	else % 3d
		im(6, x(:,:,:,end), 'QPWLS-CG phase', arg.clim), cbar
	end
end

%
% mri_phase_wi_ml()
% wi based on ML estimation
% trick: normalize by median of non-background so that beta is "universal"
%
function wi = mri_phase_wi_ml(mag, fmax)
med = median(mag(mag(:) > fmax * max(mag(:))));
wi = (mag / med).^2;

%
% phase_dercurv()
% wi * (1 - cos(yi - li))
%
function [deriv, curv] = phase_dercurv(data, li, varargin)
yi = data{1};
wi = data{2};
deriv = wi .* sin(li - yi);
curv = wi;

%
% built-in test/example
%
function [xq, mask0] = mri_phase_denoise_test(type, varargin)

% read data
yi = ir_get_data(fullfile('mri','2001-phase-data','phfit.mat'));
clim = [-0.5 1.5];

%if nargout
%	order = 1;
%else
%	order = 2;
%end
[xq mask0] = mri_phase_denoise(yi, ...
	'clim', clim, 'init', [], 'chat', 1, varargin{:});
%	'init', 5*randn(size(yi)));
%	'init', 5*ones(size(yi)));
xq = xq(:,:,end);
if im
	title 'QPWLS-CG phase (simple wi)'
end

cpu etic
xpl = mri_phase_denoise(yi, 'pl', 1, varargin{:});
cpu etoc 'PL time'
im(4, xpl, 'PL-CG phase', clim), cbar

cpu etic
xml_qpwls = mri_phase_denoise(yi, 'wi_ml', 1, varargin{:});
cpu etoc 'PWLS time'
im(5, xml_qpwls, 'QPWLS-CG (ML wj)', clim), cbar

max_percent_diff(xpl, xml_qpwls)
%max_percent_diff(xpl, xq)
nrms(xml_qpwls, xpl)
nrms(xq, xpl)
%im(4, xpl-xq), cbar

% ir_savefig fig_mr_phase_pl


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%

for i = 1:size(X, 1),
  min_d = inf;
  for k = 1:K,
    dis = X(i, :)' - centroids(k, :)';
    dis_sq = dis' * dis;
    if min_d > dis_sq,
      idx(i) = k;
      min_d = dis_sq;
    endif
  endfor
endfor

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

end

author: fewtime

repo: ML

save_path: github-repos/MATLAB/fewtime-ML


text:
function res = isCounterClockwise(p1, p2, p3, varargin)
%ISCOUNTERCLOCKWISE Compute the relative orientation of 3 points.
%
%   CCW = isCounterClockwise(P1, P2, P3);
%   Computes the orientation of the 3 points. The returns is:
%   +1 if the path P1->P2->P3 turns Counter-Clockwise (i.e., the point P3
%       is located "on the left" of the line P1-P2)
%   -1 if the path turns Clockwise (i.e., the point P3 lies "on the right"
%       of the line P1-P2)
%   0  if the point P3 is located on the line segment [P1 P2].
%
%   This function can be used in more complicated algorithms: detection of
%   line segment intersections, convex hulls, point in triangle...
%
%   CCW = isCounterClockwise(P1, P2, P3, EPS);
%   Specifies the threshold used for detecting colinearity of the 3 points.
%   Default value is 1e-12 (absolute).
%
%   Example
%   isCounterClockwise([0 0], [10 0], [10 10])
%   ans =
%       1
%   isCounterClockwise([0 0], [0 10], [10 10])
%   ans =
%       -1
%   isCounterClockwise([0 0], [10 0], [5 0])
%   ans =
%       0
%
%   See also
%   points2d, isPointOnLine, isPointInTriangle, polygonArea
%
%   References
%     Algorithm adapated from Sedgewick's book.
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2010-04-09
% Copyright 2010-2022 INRAE - Cepia Software Platform

% get threshold value
eps = 1e-12;
if ~isempty(varargin)
    eps = varargin{1};
end

% ensure all data have same size
np = max([size(p1, 1) size(p2, 1) size(p3,1)]);
if np > 1
    if size(p1,1) == 1
        p1 = repmat(p1, np, 1);
    end
    if size(p2,1) == 1
        p2 = repmat(p2, np, 1);
    end
    if size(p3,1) == 1
        p3 = repmat(p3, np, 1);
    end
end

% init with 0
res = zeros(np, 1);

% extract vector coordinates
x0  = p1(:, 1);
y0  = p1(:, 2);
dx1 = p2(:, 1) - x0;
dy1 = p2(:, 2) - y0;
dx2 = p3(:, 1) - x0;
dy2 = p3(:, 2) - y0;

% check non colinear cases
res(dx1 .* dy2 > dy1 .* dx2) =  1;
res(dx1 .* dy2 < dy1 .* dx2) = -1;

% case of colinear points
ind = abs(dx1 .* dy2 - dy1 .* dx2) < eps;
res(ind( (dx1(ind) .* dx2(ind) < 0) | (dy1(ind) .* dy2(ind) < 0) )) = -1;
res(ind(  hypot(dx1(ind), dy1(ind)) <  hypot(dx2(ind), dy2(ind)) )) =  1;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function geometry_test068 ( )

%*****************************************************************************80
%
%% GEOMETRY_TEST068 tests the SPHERE_DISTANCE routines.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  ntest = 6;

  name = [ ...
       'Atlanta, Georgia  '; ...
       'North Pole        '; ...
       'South Pole        '; ...
       'Timbuktu          '; ...
       'San Antonio, Texas'; ...
       'Savannah, Georgia ' ];
  lat_d =  [ 33, 90, -90, 16, 29, 32 ];
  lat_m =  [ 11,  0,   0, 49, 25,  5 ];
  long_d = [ 82,  0,   0,  3, 98, 81 ];
  long_m = [ 34,  0,   0,  0, 30,  6 ];
  radius = 3957.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'GEOMETRY_TEST068\n' );
  fprintf ( 1, '  POINTS_DIST_SPHERE_3D measures the distance between\n' );
  fprintf ( 1, '  two points on a sphere.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  All tests uses RADIUS = %f\n', radius );
  fprintf ( 1, '  which is the radius of the earth in miles.\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ntest-1

    lat1 = dms_to_radians ( lat_d(i), lat_m(i), 0.0 );
    long1 = dms_to_radians ( long_d(i), long_m(i), 0.0 );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Distance from %s\n', name(i,:) );

    for j = i+1 : ntest

      lat2 = dms_to_radians ( lat_d(j), lat_m(j), 0.0 );
      long2 = dms_to_radians ( long_d(j), long_m(j), 0.0 );

      dist1 = sphere_distance1 ( lat1, long1, lat2, long2, radius );
      dist2 = sphere_distance2 ( lat1, long1, lat2, long2, radius );
      dist3 = sphere_distance3 ( lat1, long1, lat2, long2, radius );

      fprintf ( 1, '             to %s is %16.8f  %16.8f  %16.8f\n', ...
        name(j,:), dist1, dist2, dist3 );

    end

  end

  return
end


text:
function calpak_test636 ( )

%*****************************************************************************80
%
%% CALPAK_TEST636 tests YEAR_TO_SCALIGER_COMMON.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST636\n' );
  fprintf ( 1, '  For a Common year,\n' );
  fprintf ( 1, '  YEAR_TO_SCALIGER_COMMON determines the Scaliger indices.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Year  Julian / Metonic / Indiction\n' );
  fprintf ( 1, '\n' );

  for y = -4713 : -4675

    y2 = y_astronomical_to_common ( y );
    sy = y_to_s_common ( y2 );
    [ c1, c2, c3, r1, r2, r3 ] = year_to_scaliger_common ( y2 );
    fprintf ( 1, '  %10s  %3d  %3d    %3d  %3d    %3d  %3d\n', ...
    sy, c1, r1, c2, r2, c3, r3 );

  end

  return
end


text:
function y = gam_lpdf(x,a,b)
%GAM_LPDF Log of Gamma probability density function (lpdf).
%
%   Y = GAM_LPDF(X,A,B) Returns the gamma lpdf with
%   shape A and inverse scale B, at the values in X.
%
%   The size of X is the common size of the input arguments. A
%   scalar input functions as a constant matrix of the same size as
%   the other inputs.
%
%   The parameterization is as in Gelman, Carlin, Stern, Dunson, Vehtari,
%   and Rubin (2013). Bayesian Data Analysis, third edition.
%
% Copyright (c) 1998-2004 Aki Vehtari

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

if nargin < 3,
  error('Requires three input arguments.');
end

%y = b.^a/gamma(a)*x^(a-1)*exp(-b*x);
y = a.*log(b)-gammaln(a)+(a-1).*log(x)-b.*x;

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function icadat = ica(fmridat_obj, varargin)
% Spatial ICA of an fmri_data object
%   - icadat = ica(fmridat_obj, [number of ICs to save])
%   - icadat is also an fmri_data object, with .dat field voxels x components
%
% :Notes:
%   - icasig = W * mixedsig
%   - icasig = icadat.dat' = W * fmridat_obj.dat'
%
% A is scaled version of fmridat_obj.dat' * icadat.dat
%
% A and W are stored in additional_info field of icadat

nic = 30;
if length(varargin) > 0, nic = varargin{1}; end

neig =  size(fmridat_obj.dat, 2);

[icasig, A, W] = icatb_fastICA(double(fmridat_obj.dat'), 'lastEig', neig, ...
    'numOfIC', nic, 'stabilization', 'on', 'verbose', 'on');

icadat = fmri_data;

icadat.dat = icasig';
icadat.mask = fmridat_obj.mask;
icadat.volInfo = fmridat_obj.volInfo;
icadat.removed_voxels = fmridat_obj.removed_voxels;
icadat.removed_images = 0;

icadat.additional_info{1} = A;
icadat.additional_info{2} = W;
icadat.history{1} = 'Spatial ICA of fmri_data object.';
icadat.history{2} = 'Stored A and W matrices in cells 1,2 of additional_info';
icadat.source_notes = 'fastICA algorithm on images (spatial)';

% Show on orthviews

tmpicadat = icadat;
for i = 1:size(tmpicadat.dat, 2)
    d = tmpicadat.dat(:, i);
    tmpicadat.dat(d < prctile(d, 85) & d > prctile(d, 15), i) = 0;

end

orthviews(tmpicadat);
spm_orthviews_white_background

end

% m1 = prctile(icadat.dat(:), 5);
% m2 = prctile(icadat.dat(:), 95);
%
% spm_orthviews('window', 1:min(24, size(icadat.dat, 2)), [m1 m2]);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% notes

% How lookback/horizon changes accuracy
% exploratory data analysis
% How does sector affect price movements
% Try a mixture of SVM models?

%%

SNPIndices = zeros(6356, 1);
otherIndices = [];

for i=1:505
    boolIdx = strcmp(Symbol(i), names);
    if (any(boolIdx))
        SNPIndices = (SNPIndices | boolIdx);
        otherIndices = [otherIndices i];
    end
end

%%
% have names and Sector, Names

newSectors = cell(494, 1);

for i=1:494

    linInd = find(strcmp(Symbol(i), NAMES));
    newSectors{linInd} = Sector{i};

end

%% data loading and plotting tests

today = 30;
lookback = 10;
horizon = 2;
[X, y, ~, ~] = loadData(today, lookback, horizon);

plot(today-lookback+1:today, X(1,:));
hold on;
plot(today+horizon, y(1), 'r*', 'MarkerSize', 12);

%%
%% try some clustering
k = 4;
[idx, means] = kmeans(PRICE_HISTORY, k);

%%
figure();
for i=1:k
    subplot(2, 3, i);
    hold on;
    plot(PRICE_HISTORY(idx == i, :)');
end

%%
figure();
for i=1:k
    subplot(2, 3, i);
    hold on;
    plot(means(i,:));
end

%%

[X, y] = loadData2(30, 10);

author: gudbrandtandberg


text:
function [FV] = mesh_grow(FV,origin,dist),

% mesh_grow - explode vertices of mesh by specific distance
%
% FV = mesh_grow(FV,origin,dist)
%
% FV is a struct with fields:
%
% FV.vertices   - Nx3 matrix of Cartesian vertex coordindates (X,Y,Z)
% FV.faces      - Mx3 matrix of triangulation of FV.vertices
%
% origin        - 1x3 row vector, usually (0,0,0)
%
% dist          - how far to explode the mesh away from the origin;
%                 this distance is relative to current distance from
%                 the origin, not the total distance from the origin.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:57 $

% Licence:  GNU GPL, no implied or express warranties
% History:  10/2002, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    xo = origin(1); yo = origin(2); zo = origin(3);

    Nvert = size(FV.vertices,1);

    fprintf('...mesh explosion...'); tic;

    for v = 1:Nvert,

        x = FV.vertices(v,1);
        y = FV.vertices(v,2);
        z = FV.vertices(v,3);

        % Find direction cosines for line from centre to vertex
        d = sqrt( (x-xo)^2 + (y-yo)^2 + (z-zo)^2 );

        l = (x-xo)/d; % cos alpha
        m = (y-yo)/d; % cos beta
        n = (z-zo)/d; % cos gamma

        % now decrease d by dist
        d = d + dist;

        % locate vertex at this new distance
        x = (l * d) + xo;
        y = (m * d) + yo;
        z = (n * d) + zo;

        FV.vertices(v,:) = [ x y z ];
    end

    t = toc; fprintf('...done (%5.2f sec)\n',t);

return

author: PatternRecognition

repo: OpenBMI


text:
function plotData(x, y)
%PLOTDATA Plots the data points x and y into a new figure
%   PLOTDATA(x,y) plots the data points and gives the figure axes labels of
%   population and profit.

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the training data into a figure using the
%               "figure" and "plot" commands. Set the axes labels using
%               the "xlabel" and "ylabel" commands. Assume the
%               population and revenue data have been passed in
%               as the x and y arguments of this function.
%
% Hint: You can use the 'rx' option with plot to have the markers
%       appear as red crosses. Furthermore, you can make the
%       markers larger by using plot(..., 'rx', 'MarkerSize', 10);

figure; % open a new figure window

plot(x, y, 'rx', 'MarkerSize', 10);
xlabel('Population of City in 10,000s');
ylabel('Profit in $10,000s');

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

end

author: sfvsfv

repo: Mathematical-modeling


text:
function varargout = iff_internal(varargin)
X = varargin{1};
Y = varargin{2};

if nargin == 2
    zero_tolerance = 1e-5;
else
    zero_tolerance = abs(varargin{3});
end

% Normalize data
if isa(Y,'constraint')
    Y=lmi(Y,[],[],1);
end
if isa(X,'constraint')
    X=lmi(X,[],[],1);
end
if isa(X,'lmi') & isa(Y,'sdpvar')
    temp = X;
    X = Y;
    Y = temp;
end

if isa(X,'sdpvar') & isa(Y,'sdpvar')
    % user presumably works with binary variables
    varargout{1} = [Y == X];

elseif (isa(X,'sdpvar') & isa(Y,'lmi'))
    % Binary variable iff a constraint holds
    switch settype(Y)
        case 'elementwise'  % binary X <--> Y(:)>=0
            varargout{1} = binary_iff_lp(X,-sdpvar(Y),zero_tolerance);
        case 'equality'
            varargout{1} = binary_iff_equality(X,sdpvar(Y),zero_tolerance);
        case 'multiple'
            Y1 = Y(find(is(Y,'equality')));
            Y2 = Y(find(is(Y,'elementwise')));
            % Glue using intermediate variables holding truth of each set
            % of constraints. X is true iff both di are true and v.v.
            % Hence, we make a recursive call with two new models
            binvar d1 d2
            C1 = iff_internal(d1,Y1,zero_tolerance);
            C2 = iff_internal(d2,Y2,zero_tolerance);
            C3 = [X <= d1,X <= d2,X >= d1+d2-1];
            varargout{1} = [C1, C2, C3];
          %  varargout{1} = [iff_internal(X,Y1),iff_internal(X,Y2)];
        otherwise
            error('IFF not implemented for this case');
    end
elseif isa(X,'lmi') & isa(Y,'lmi')
    % Constraint holds iff constraint holds.
    % Glue using an intermediate variable
    binvar d
    C1 =  iff_internal(d,X,zero_tolerance);
    C2 =  iff_internal(d,Y,zero_tolerance);
    varargout{1} = [C1,C2];
end

function F = binary_iff_lp(X,f,zero_tolerance)
% X == 1    <=>   f<=0
[M,m,infbound] = derivebounds(f);
if infbound
    warning('You have unbounded variables in IFF leading to a lousy big-M relaxation.');
end

[nf,mf]=size(f);
if nf*mf==1
    F = linearnegativeconstraint_iff_binary(f,X,M,m,zero_tolerance);
else
    f = reshape(f,nf*mf,1);
    di = binvar(nf*mf,1);
    F = linearnegativeconstraint_iff_binary(f,di,M,m,zero_tolerance);
    if length(X)==1
        % X is true if any di
        F = [F, X>=sum(di)-length(di)+1, X <= di];
    else
        % This must be a vectorized X(i) iff f(i)
        F = [F, di == X];
    end

    % di=0 means the ith hypeplane is violated
    % X=1 means we are in the polytope
    % F  = [f <= M*(1-X), f>=eps+(m-eps).*di, X>=sum(di)-length(di)+1, X <= di];

    % Add some cuts for c < a'x+b < d
    [bA] = getbase(f);
    b = bA(:,1);
    A = bA(:,2:end);
    S = zeros(0,length(di));
    for i = 1:length(b)
        j = findrows(abs(A),abs(A(i,:)));
        j = j(j > i);
        if length(j)==1
            S(end+1,[i j]) = 1;
        end
    end
    if size(S,1) > 0
        % Add cut cannot be outside both constraints
        F = F + (S*di >= 1);
    end
end

function F = binary_iff_equality(X,Y,zero_tolerance)

% Things like iff(x,y==1) is sent as X, 1-Y
if isLogicalVector(X) && isLogicalVector(Y)
    if isequal(getbase(Y),[0 -1]) | isequal(getbase(Y),[-1 1])
        Y = -Y;
    end
     F = [1-X == Y];
     return
 end
Y = Y(:);
d = binvar(length(Y),3);
% We have to model every single line of equality.
% di1 : Yi < -eps
% di2 : -eps < Yi < eps
% di3 : eps < Yi
C = sum(d,2)==1;
[M,m,infbounds] = derivebounds(Y);
for i = 1:length(Y)
    if all(ismember(depends(Y(i)),yalmip('binvariables'))) & all(getbase(Y(i))==fix(getbase(Y(i))))
        eps = 0.5;
    else
        eps = 0;
    end
    C1 =  binary_iff_lp(d(i,1),Y(i)+eps,abs(zero_tolerance));             % Y <-eps
    C2 =  binary_iff_lp(d(i,2),[Y(i)-eps;-eps-Y(i)],abs(zero_tolerance)); % -eps < Y < eps
    C3 =  binary_iff_lp(d(i,3),eps-Y(i),abs(zero_tolerance));             % eps < Y
    C = [C, C1,C2,C3];
    % Cut off some trivial cases
    if m(i)>=0
        C = [C,d(i,1)==0];
    elseif M(i)<=0
        C = [C,d(i,3)==0];
    end
end
% X is true if all di2 are true and v.v
F = [C,X <= d(:,2), X >= sum(d(:,2))-length(Y)+1];

function yes = isLogicalVector(X)
yes = 0;
X = sdpvar(X);
B = getbase(X);
if all(ismember(depends(X),yalmip('binvariables')))
    b = B(:,1);
    if all(b==0 | b==1)
        A = B(:,2:end);
        if all(A == 0 | A == 1 | A == -1)
            if all(sum(A | A,2)==1)
                if all(sum(B,2)<= 1)
                    yes = 1;
                end
            end
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%DELTA2TR Convert differential motion  to a homogeneous transform
%
% T = DELTA2TR(D) is a homogeneous transform (4x4) representing differential
% translation and rotation. The vector D=(dx, dy, dz, dRx, dRy, dRz)
% represents an infinitessimal motion, and is an approximation to the spatial
% velocity multiplied by time.
%
% See also TR2DELTA.

% Copyright (C) 1993-2015, 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 delta = delta2tr(d)
    d = d(:);
    delta = eye(4,4) + [skew(d(4:6)) d(1:3); 0 0 0 0];

author: Allopart

repo: rbpf-gmapping


text:
function t_most_suc(quiet, create_plots, create_pdfs, savedir)
%T_MOST_SUC  Tests of stochastic unit commitment optimizations.
%
%   T_MOST_SUC(QUIET, CREATE_PLOTS, CREATE_PDFS, SAVEDIR)
%   Can generate summary plots and save them as PDFs in a directory of
%   your choice.
%   E.g. t_most_suc(0, 1, 1, '~/Downloads/suc_plots')

%   MOST
%   Copyright (c) 2015-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MOST.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/most for more info.

if nargin < 4
    savedir = '.';              %% save in current working directory by default
    if nargin < 3
        create_pdfs = 0;        %% do NOT save plots to PDF files
        if nargin < 2
            create_plots = 0;   %% do NOT create summary plots of results
            if nargin < 1
                quiet = 0;      %% verbose by default
            end
        end
    end
end
if create_plots
    if create_pdfs
        fname = 'suc-ex';
    else
        fname = '';
    end
    pp = 0;     %% plot counter
end

solvers = {'CPLEX', 'GLPK', 'GUROBI', 'MOSEK', 'OT'};
fcn = {'cplex', 'glpk', 'gurobi', 'mosek', 'intlinprog'};
% solvers = {'CPLEX'};
% fcn = {'cplex'};
% solvers = {'OT'};
% fcn = {'intlinprog'};
% solvers = {'GUROBI'};
% fcn = {'gurobi'};
% solvers = {'DEFAULT'};
% fcn = {'most'};
ntests = 37;
t_begin(ntests*length(solvers), quiet);

if quiet
    verbose = 0;
else
    verbose = 0;
end
% verbose = 2;

if have_feature('octave')
    if have_feature('octave', 'vnum') >= 4
        file_in_path_warn_id = 'Octave:data-file-in-path';
    else
        file_in_path_warn_id = 'Octave:load-file-in-path';
    end
    s1 = warning('query', file_in_path_warn_id);
    warning('off', file_in_path_warn_id);
end

casefile = 'ex_case3b';
solnfile =  't_most_suc_soln';
soln = load(solnfile);
mpopt = mpoption;
mpopt = mpoption(mpopt, 'out.gen', 1);
mpopt = mpoption(mpopt, 'verbose', verbose);
% mpopt = mpoption(mpopt, 'opf.violation', 1e-6, 'mips.gradtol', 1e-8, ...
%         'mips.comptol', 1e-8, 'mips.costtol', 1e-8);
mpopt = mpoption(mpopt, 'model', 'DC');
mpopt = mpoption(mpopt, 'most.price_stage_warn_tol', 10);

%% solver options
if have_feature('cplex')
    %mpopt = mpoption(mpopt, 'cplex.lpmethod', 0);       %% automatic
    %mpopt = mpoption(mpopt, 'cplex.lpmethod', 1);       %% primal simplex
    mpopt = mpoption(mpopt, 'cplex.lpmethod', 2);       %% dual simplex
    %mpopt = mpoption(mpopt, 'cplex.lpmethod', 3);       %% network simplex
    %mpopt = mpoption(mpopt, 'cplex.lpmethod', 4);       %% barrier
    mpopt = mpoption(mpopt, 'cplex.opts.mip.tolerances.mipgap', 0);
    mpopt = mpoption(mpopt, 'cplex.opts.mip.tolerances.absmipgap', 0);
    mpopt = mpoption(mpopt, 'cplex.opts.threads', 2);
end
if have_feature('glpk')
    mpopt = mpoption(mpopt, 'glpk.opts.mipgap', 0);
    mpopt = mpoption(mpopt, 'glpk.opts.tolint', 1e-10);
    mpopt = mpoption(mpopt, 'glpk.opts.tolobj', 1e-10);
end
if have_feature('gurobi')
    %mpopt = mpoption(mpopt, 'gurobi.method', -1);       %% automatic
    %mpopt = mpoption(mpopt, 'gurobi.method', 0);        %% primal simplex
    mpopt = mpoption(mpopt, 'gurobi.method', 1);        %% dual simplex
    %mpopt = mpoption(mpopt, 'gurobi.method', 2);        %% barrier
    mpopt = mpoption(mpopt, 'gurobi.threads', 2);
    mpopt = mpoption(mpopt, 'gurobi.opts.MIPGap', 0);
    mpopt = mpoption(mpopt, 'gurobi.opts.MIPGapAbs', 0);
end
if have_feature('mosek')
    sc = mosek_symbcon;
    %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_FREE);            %% default
    %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_INTPNT);          %% interior point
    %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_PRIMAL_SIMPLEX);  %% primal simplex
    mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);     %% dual simplex
    %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_FREE_SIMPLEX);    %% automatic simplex
    %mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_X', 0);
    mpopt = mpoption(mpopt, 'mosek.opts.MSK_IPAR_MIO_NODE_OPTIMIZER', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);
    mpopt = mpoption(mpopt, 'mosek.opts.MSK_IPAR_MIO_ROOT_OPTIMIZER', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);
    mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_ABS_RELAX_INT', 1e-9);
    %mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_REL_RELAX_INT', 0);
    mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_REL_GAP', 0);
    mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_ABS_GAP', 0);
end
if have_feature('intlinprog')
    %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'interior-point');
    %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'active-set');
    %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'simplex');
    mpopt = mpoption(mpopt, 'linprog.Algorithm', 'dual-simplex');
    %mpopt = mpoption(mpopt, 'intlinprog.RootLPAlgorithm', 'primal-simplex');
    mpopt = mpoption(mpopt, 'intlinprog.RootLPAlgorithm', 'dual-simplex');
    mpopt = mpoption(mpopt, 'intlinprog.TolCon', 1e-9);
    mpopt = mpoption(mpopt, 'intlinprog.TolGapAbs', 0);
    mpopt = mpoption(mpopt, 'intlinprog.TolGapRel', 0);
    mpopt = mpoption(mpopt, 'intlinprog.TolInteger', 1e-6);
    %% next line is to work around a bug in intlinprog
    % (Technical Support Case #01841662)
    mpopt = mpoption(mpopt, 'intlinprog.LPPreprocess', 'none');
end
if ~verbose
    mpopt = mpoption(mpopt, 'out.all', 0);
end
% mpopt = mpoption(mpopt, 'out.all', -1);

%% define named indices into data matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
[CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
    CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
    CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
    CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
    CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
    CT_MODCOST_X] = idx_ct;

%% load base case file
mpc = loadcase(casefile);

nb = size(mpc.bus, 1);
nl = size(mpc.branch, 1);
ng = size(mpc.gen, 1);

xgd = loadxgendata('ex_xgd_uc', mpc);
[iwind, mpc, xgd] = addwind('ex_wind_uc', mpc, xgd);
profiles_d = getprofiles('ex_wind_profile_d', iwind);
profiles_d = getprofiles('ex_load_profile', profiles_d);
profiles_s = getprofiles('ex_wind_profile', iwind);
profiles_s = getprofiles('ex_load_profile', profiles_s);
nt = size(profiles_d(1).values, 1);

mpc0 = mpc;
xgd0 = xgd;

for s = 1:length(solvers)
    if ~have_feature(fcn{s})    %% check if we have the solver
        t_skip(ntests, sprintf('%s not installed', solvers{s}));
    else
        mpopt = mpoption(mpopt, 'opf.dc.solver', solvers{s});
        mpopt = mpoption(mpopt, 'most.solver', mpopt.opf.dc.solver);
        mpopt = mpoption(mpopt, 'most.storage.cyclic', 1);

        mpc = mpc0;
        xgd = xgd0;

        t = sprintf('%s : deterministic : ', solvers{s});
        mdi = loadmd(mpc, nt, xgd, [], [], profiles_d);
        mdo = most(mdi, mpopt);
        ms = most_summary(mdo);
        t_ok(mdo.QP.exitflag > 0, [t 'success']);
        ex = soln.determ;
        t_is(ms.f, ex.f, 8, [t 'f']);
        t_is(ms.Pg, ex.Pg, 8, [t 'Pg']);
        t_is(ms.Rup, ex.Rup, 8, [t 'Rup']);
        t_is(ms.Rdn, ex.Rdn, 8, [t 'Rdn']);
        t_is(ms.Pf, ex.Pf, 8, [t 'Pf']);
        t_is(ms.u, ex.u, 8, [t 'u']);
        t_is(ms.lamP, ex.lamP, 8, [t 'lamP']);
        t_is(ms.muF, ex.muF, 8, [t 'muF']);
        % determ = most_summary(mdo);
        if create_plots
            pp = pp + 1;
            plot_case('Base : Deterministic', mdo, ms, 500, 150, savedir, pp, fname);
        end
        % keyboard;

        t = sprintf('%s : individual trajectories : ', solvers{s});
        transmat_s = cell(1, nt);
        I = speye(3);
        [transmat_s{:}] = deal(I);
        transmat_s{1} = [ 0.158655253931457; 0.682689492137086; 0.158655253931457 ];
        mdi = loadmd(mpc, transmat_s, xgd, [], [], profiles_s);
        mdi = filter_ramp_transitions(mdi, 0.1);
        mdo = most(mdi, mpopt);
        ms = most_summary(mdo);
        t_ok(mdo.QP.exitflag > 0, [t 'success']);
        ex = soln.transprob1;
        t_is(ms.f, ex.f, 5, [t 'f']);
        t_is(ms.Pg, ex.Pg, 6, [t 'Pg']);
        t_is(ms.Rup, ex.Rup, 6, [t 'Rup']);
        t_is(ms.Rdn, ex.Rdn, 6, [t 'Rdn']);
        t_is(ms.Pf, ex.Pf, 6, [t 'Pf']);
        t_is(ms.u, ex.u, 8, [t 'u']);
        % t_is(ms.lamP, ex.lamP, 5, [t 'lamP']);
        % t_is(ms.muF, ex.muF, 5, [t 'muF']);
        % transprob1 = most_summary(mdo);
        if create_plots
            pp = pp + 1;
            plot_case('Individual Trajectories', mdo, ms, 500, 150, savedir, pp, fname);
        end
        % keyboard;

        t = sprintf('%s : full transition probabilities : ', solvers{s});
%        transmat_sf = transmat_s;
        transmat_sf = ex_transmat(nt);
%         transmat_sf = cell(1, nt);
%         T = [ 0.158655253931457; 0.682689492137086; 0.158655253931457 ];
%         [transmat_sf{:}] = deal(T * ones(1,3));
%         transmat_sf{1} = T;
        mdi = loadmd(mpc, transmat_sf, xgd, [], [], profiles_s);
%        mdi = filter_ramp_transitions(mdi, 0.9);
        mdo = most(mdi, mpopt);
        ms = most_summary(mdo);
        t_ok(mdo.QP.exitflag > 0, [t 'success']);
        ex = soln.transprobfull;
        t_is(ms.f, ex.f, 3, [t 'f']);
        t_is(ms.Pg, ex.Pg, 3, [t 'Pg']);
        t_is(ms.Rup, ex.Rup, 3, [t 'Rup']);
        t_is(ms.Rdn, ex.Rdn, 3, [t 'Rdn']);
        t_is(ms.Pf, ex.Pf, 8, [t 'Pf']);
        t_is(ms.u, ex.u, 8, [t 'u']);
        % t_is(ms.lamP, ex.lamP, 5, [t 'lamP']);
        % t_is(ms.muF, ex.muF, 5, [t 'muF']);
        % transprobfull = most_summary(mdo);
        if create_plots
            pp = pp + 1;
            plot_case('Full Transition Probabilities', mdo, ms, 500, 150, savedir, pp, fname);
        end
        % keyboard;

        t = sprintf('%s : full transition probabilities + cont : ', solvers{s});
        mdi = loadmd(mpc, transmat_sf, xgd, [], 'ex_contab', profiles_s);
%        mdi = filter_ramp_transitions(mdi, 0.9);
        mdo = most(mdi, mpopt);
        ms = most_summary(mdo);
        t_ok(mdo.QP.exitflag > 0, [t 'success']);
        ex = soln.transprobcont;
        t_is(ms.f, ex.f, 3.5, [t 'f']);
        t_is(ms.Pg, ex.Pg, 6, [t 'Pg']);
        t_is(ms.Rup, ex.Rup, 6, [t 'Rup']);
        t_is(ms.Rdn, ex.Rdn, 6, [t 'Rdn']);
        t_is(ms.Pf, ex.Pf, 6, [t 'Pf']);
        t_is(ms.u, ex.u, 8, [t 'u']);
        % t_is(ms.lamP, ex.lamP, 5, [t 'lamP']);
        % t_is(ms.muF, ex.muF, 5, [t 'muF']);
        % transprobcont = most_summary(mdo);
        if create_plots
            pp = pp + 1;
            plot_case('+ Contingencies', mdo, ms, 500, 150, savedir, pp, fname);
        end
        % keyboard;

        t = sprintf('%s : + storage : ', solvers{s});
        if mpopt.out.all
            fprintf('Add storage\n');
        end
        [iess, mpc, xgd, sd] = addstorage('ex_storage', mpc, xgd);
        mdi = loadmd(mpc, transmat_sf, xgd, sd, 'ex_contab', profiles_s);
        mdo = most(mdi, mpopt);
        ms = most_summary(mdo);
        t_ok(mdo.QP.exitflag > 0, [t 'success']);
        ex = soln.wstorage;
        t_is(ms.f, ex.f, 3, [t 'f']);
        t_is(ms.Pg, ex.Pg, 3, [t 'Pg']);
        t_is(ms.Rup, ex.Rup, 3, [t 'Rup']);
        t_is(ms.Rdn, ex.Rdn, 4, [t 'Rdn']);
        t_is(ms.Pf, ex.Pf, 3, [t 'Pf']);
        t_is(ms.u, ex.u, 8, [t 'u']);
        % t_is(ms.lamP, ex.lamP, 5, [t 'lamP']);
        % t_is(ms.muF, ex.muF, 5, [t 'muF']);
        % wstorage = most_summary(mdo);
        if create_plots
            pp = pp + 1;
            plot_case('+ Storage', mdo, ms, 500, 150, savedir, pp, fname);
            create_plots = 0;   %% don't do them again
        end
        % keyboard;
    end
end

if have_feature('octave')
    warning(s1.state, file_in_path_warn_id);
end

t_end;

% save t_most_suc_soln determ transprob1 transprobfull transprobcont wstorage
% determ.u
% transprob1.u
% transprobcont.u
% transprobfull.u
% wstorage.u

function h = plot_case(label, md, ms, maxq, maxp, mypath, pp, fname)

if nargin < 8
    fname = '';
end

%% colors:  blue     red               yellow           purple            green
cc = {[0 0.45 0.74], [0.85 0.33 0.1], [0.93 0.69 0.13], [0.49 0.18 0.56], [0.47 0.67 0.19]};

ig = (1:3)';
id = 4;
iw = 5;
is = 6;

subplot(3, 1, 1);
md.mpc = rmfield(md.mpc, 'genfuel');
plot_uc(md, [], 'title', label);
ylabel('Unit Commitment', 'FontSize', 16);
ah = gca;
ah.YAxisLocation = 'left';

subplot(3, 1, 2);
x = (1:ms.nt)';
Pg = md.results.ExpectedDispatch;
y1 = Pg(ig, :)';
if ms.ng == 6
    y1 = [y1 max(-Pg(is, :), 0)' max(Pg(is, :), 0)'];
end
y2 = -sum(Pg([id; iw], :), 1)';
[ah1, h1, h2] = plotyy(x, y1, x, y2);
axis(ah1(1), [0.5 12.5 0 maxq]);
axis(ah1(2), [0.5 12.5 0 maxq]);
% ah1(1).XLim = [0.5 12.5];
% ah1(2).XLim = [0.5 12.5];
% ah1(1).YLim = [0 300];
% ah1(2).YLim = [0 450];
ah1(1).YTickMode = 'auto';
ah1(2).YTickMode = 'auto';
ah1(1).XTick = 1:12;
nn = 3;
for j = 1:3
    h1(j).LineWidth = 2;
    h1(j).Color = cc{j};
end
if ms.ng == 6
    h1(4).LineWidth = 2;
    h1(4).Color = cc{5};
    h1(4).LineStyle = ':';
    h1(5).LineWidth = 2;
    h1(5).Color = cc{5};
end
h2.LineWidth = 2;
h2.Color = cc{4};
h2.LineStyle = ':';
ah1(2).YColor = cc{4};
%title('Generation & Net Load', 'FontSize', 16);
ylabel(ah1(1), 'Generation, MW', 'FontSize', 16);
ylabel(ah1(2), 'Net Load, MW', 'FontSize', 16);
xlabel('Period', 'FontSize', 16);
set(ah1(1), 'FontSize', 14);
set(ah1(2), 'FontSize', 14);
if ms.ng == 6
    legend('Gen 1', 'Gen 2', 'Gen 3', 'Storage Charge', 'Storage Discharge', 'Location', [0.7 0.6 0 0]);
else
    legend('Gen 1', 'Gen 2', 'Gen 3', 'Location', [0.7 0.58 0 0]);
end

subplot(3, 1, 3);
if length(size(ms.lamP)) == 4
    elamP = sum(sum(ms.lamP, 4), 3) ./ (ones(ms.nb,1) * md.StepProb);
%     elamP = sum(sum(ms.lamP, 4), 3);
elseif length(size(ms.lamP)) == 3
    elamP = sum(ms.lamP, 3);
else
    elamP = ms.lamP;
end

y1 = elamP';
plot(x, y1, 'LineWidth', 2);
% title('Nodal Price', 'FontSize', 16);
ylabel('Nodal Price, $/MWh', 'FontSize', 16);
xlabel('Period', 'FontSize', 16);
axis([0.5 12.5 0 maxp]);
ah = gca;
set(ah, 'FontSize', 14);
ah.XTick = 1:12;
legend('Bus 1', 'Bus 2', 'Bus 3', 'Location', [0.7 0.28 0 0]);

if nargin > 7 && ~isempty(fname)
    h = gcf;
    set(h, 'PaperSize', [11 8.5]);
    set(h, 'PaperPosition', [0.25 0.25 10.5 8]);
    print('-dpdf', fullfile(mypath, sprintf('%s-%d', fname, pp)));
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function test_bug1514

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY ft_spike_select ft_selectdata

% the following structure corresponds to the one explained in
% ft_datatype_spike

spike = [];
spike.label = {'unit1' 'unit2' 'unit3'};
spike.timestamp = {randn(1,504) randn(1,50) randn(1,101)};
spike.time = {randn(1,504) randn(1,50) randn(1,101)};
spike.trial = {ceil(rand(1,504)*100) ceil(rand(1,50)*100) ceil(randn(1,101)*100)};
spike.trialtime = 3*randn(100,2);
spike.waveform = {randn(1,32,504) randn(1,32,50) randn(1,32,101)};
spike.waveformdimord = '{chan}_lead_time_spike';
spike.fourierspctrm = {randn(504,2,20), randn(50,2,20), randn(101,2,20)};
spike.fourierspctrmdimord = '{chan}_spike_lfplabel_freq';
spike.lfplabel = {'lfpchan1', 'lfpchan2'};
spike.freq = 1:20;

cfg = [];
cfg.trials = 1:50;
output = ft_selectdata(cfg, spike);

cfg = [];
cfg.channel = 2;
output = ft_selectdata(cfg, spike);

cfg = [];
cfg.latency = 'prestim';
output = ft_selectdata(cfg, spike);

cfg = [];
cfg.trials = 1:50;
cfg.channel = 2;
cfg.latency = 'prestim';
output = ft_selectdata(cfg, spike);

assert(all(output.time{1}<=0));
assert(all(output.trial{1}>=1 & output.trial{1}<=50));
assert(isequal(output.label, {'unit2'}));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Example of using KST class for interfacing with KUKA iiwa robots

% Drawing a sequare with a pen using the KUKA iiwa robot

% An example script, it is used to show how to use the different
% functions of the KUKA Sunrise Toolbox

% First start the server on the KUKA iiwa controller
% Then run the following script in Matlab

% be careful when runnning the script, be sure that no human, nor obstacles
% are around the robot

% This example works with Sunrise application version KST_1.7  and higher.
% Copyright Mohammad SAFEEA, 26th-June-2018

close all;clear;clc;
warning('off')
%% Create the robot object
ip='172.31.1.147'; % The IP of the controller
arg1=KST.LBR7R800; % choose the robot iiwa7R800 or iiwa14R820
arg2=KST.Medien_Flansch_elektrisch; % choose the type of flange
Tef_flange=eye(4); % transofrm matrix of EEF with respect to flange
iiwa=KST(ip,arg1,arg2,Tef_flange); % create the object

%% Start a connection with the server
flag=iiwa.net_establishConnection();
if flag==0
    return;
end
pause(1);
disp('Doing some stuff')

%% move to initial position
jPos_init={0,pi*20/180,0,-pi*70/180,0,pi*90/180,0}; % initial confuguration
VEL=0.15; % relative velocity
iiwa.movePTPJointSpace(jPos_init, VEL); % point to point motion in joint space

%% Get position roientation of end effector
disp('Cartesian position')
Pos=iiwa.getEEFPos();
disp(Pos);
z_1=Pos{3}; % save initial hight level
z0=448+3; % go to writing position
VEL=50; % velocity of the end effector, mm/sec

%% Insert the Pen at box level
Pos{3}=z0; %% first point
iiwa.movePTPLineEEF(Pos, VEL);

disp=50*2; % length of the side of the sequare

%% move in x positive direction
Pos{1}=Pos{1}+disp;
iiwa.movePTPLineEEF(Pos, VEL)

%% move in y negative direction
Pos{2}=Pos{2}-disp;
iiwa.movePTPLineEEF(Pos, VEL)

%% move in x negative direction
Pos{1}=Pos{1}-disp;
iiwa.movePTPLineEEF(Pos, VEL)

%% move in y positive direction
Pos{2}=Pos{2}+disp;
iiwa.movePTPLineEEF(Pos, VEL)

%% Go back to initial position
Pos{3}=z_1;
iiwa.movePTPLineEEF(Pos, VEL)

%% turn off the server
iiwa.net_turnOffServer();

repo: KST-Kuka-Sunrise-Toolbox