text:
% EJEMPLO SENCILLO DE INTERFAZ ENTRE ARTE Y COPPELIA.

function irb140_abb_pick_and_place()
    % init basic data. Robots and number of collision objects.
    coppelia = [];
    coppelia.n_collision_objects = 0;
    coppelia.robot.name = 'R1';
    coppelia.robot.n_joints = 6;
    coppelia.robot.end_effector.n_joints = 2;

    coppelia = coppelia_start(coppelia);
    % escribid una funcion que este entre coppelia start y coppelia_stop
    pick_and_place(coppelia)

    coppelia_stop(coppelia);
end

function pick_and_place(coppelia)
%load in arte
global robot
robot = load_robot('ABB', 'IRB140');

    %T1, pre pick point
    T1 = [0  1  0 0.52;
         1  0  0  -0.1;
          0  0 -1 0.45;
          0  0  0 1];
    % pick point
    T2 = [0  1 0 0.52;
          1 0 0  -0.1;
         0 0 -1 0.33;
         0 0 0 1];
    % Release point
    T3 = [-1 0 0 0.1;
         0 1 0  -0.5;
         0 0 -1 0.6;
          0 0 0 1];
    T4 = [-1 0 0 0.1;
         0 1 0  -0.5;
        0 0 -1 0.5;
        0 0 0 1];
    % Podeis ir a una posicion q arbitraria
    q = [0.1 0.1 0.1 0.8 0.8 0.8]';
    set_joint_target_trajectory(coppelia, q)

    %podeis ir a una posicion que sea la cinematica inversa deseada
    qinv = inversekinematic(robot, T1);
    set_joint_target_trajectory(coppelia, qinv(:,1))
    % Esto abre la pinza
    open_gripper(coppelia);
    % podéis esparar un múltiplo del dt de Coppelia (step size de simulación de
    % coppelia, tip. 50 ms. Esperar a que se abra la pinza
    coppelia_wait(coppelia, 10)
    %podeis ir a una posicion que sea la cinematica inversa deseada
    qinv = inversekinematic(robot, T2);
    set_joint_target_trajectory(coppelia, qinv(:,1))

    % cerrar pinza
    close_gripper(coppelia);
    % Esperar a que se cierre
    coppelia_wait(coppelia, 5)
    % subir
    qinv = inversekinematic(robot, T1);
    set_joint_target_trajectory(coppelia, qinv(:,1))

    qinv = inversekinematic(robot, T3);
    set_joint_target_trajectory(coppelia, qinv(:,1))

    qinv = inversekinematic(robot, T4);
    set_joint_target_trajectory(coppelia, qinv(:,1))

    % Esto abre la pinza
    open_gripper(coppelia);
    coppelia_wait(coppelia, 5)

    % VOLVEMOS A la posicion inicial
    set_joint_target_trajectory(coppelia, q)
end

repo: ARTE


text:
function [r]=RMSE(err)
% r = RMSE(err)
% author: Liang Xiong (lxiong@cs.cmu.edu)

r = sqrt(mean(err(:).^2));

author: andrewssobral

repo: lrslibrary


text:
function [z,ze]=doubleLengthSum(x,xe,y,ye)
%%DOUBLELENGTHSUM Given two doublelength values, which could, for example,
%       be returned by exactPairMult or exactPairSum, find their sum as
%       another doublelength value.
%
%INPUTS: x,xe A real doublelength floating point number. The exact number
%              is x+xe (added with infinite precision) and
%              abs(xe)<=abs(x+xe)*2^(-t)/(1+2^(-t)), where t is the number
%              of bits in the mantissa, which is 53, since it is assumed
%              that these values are floating point doubles.
%        y, ye A second real doublelength floating point number that should
%              be added to x,xe.
%
%OUTPUTS: z, ze The doublelength floating point number that is the sum of
%               (x,xe) and (y,ye). If the inputs satisfy the proper
%               splitting of the number into two parts, (the relation in
%               size between abs(xe) and abs(x+xe) and similarly for y and
%               ye), then this will also satisfy the splitting of the
%               numbers. This does not necessarily hold if denormalized
%               numbers are encountered.
%
%This function implements the Add2 algorithm of [1].
%
%Note that this assumes that the processor rounding mode has been set to
%round to "nearest," which is the default in Matlab, but which can be
%changed using the function setProcRoundingMode.
%
%EXAMPLE:
% [z,ze]=doubleLengthSum(1,1e-30,1e20,1e-1)
%One gets z=1e20 and ze=1.1. The 1e-30 component is truncated away, because
%it is too small to represent with two doubles.
%
%REFERENCES:
%[1] T. J. Dekker, "A Floating Point Technique for Extending the Available
%    Precision," Numerische Mathematik, vol. 18, no. 3, Jun. 1971, pp.
%    224-242.
%
%November 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

r=x+y;
if(abs(x)>abs(y))
    s=x-r+y+ye+xe;
else
    s=y-r+x+xe+ye;
end

z=r+s;
ze=r-z+s;

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function d = dot(obj1,obj2)
% inner product between embeddings
%
% Syntax
%   d = dot(e1,e2)
%
% Input
%  e1,e2 - @embedding
%
% Output
%  d - double
%

d = 0;
for i = 1:length(obj1.u)

  id = -(1:obj1.rank(i));
  d = d + EinsteinSum(obj1.u{i},id,obj2.u{i},id);

end

end

author: mtex-toolbox

repo: mtex

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


text:
function [status, results] = mrtrix_tensor_generate (in_file, out_file, b_file, verbose)

%
% Calculate diffusion tensors.
%
% Parameters
% ----------
% in_file: The name of a diffusion file in .mif format
% out_file: The name of the resulting dti file in .mif format
% b_file: The name of a mrtrix format gradient vector file (default: 'encoding.b')
%
% Returns
% -------
% status: whether (0) or not (1) the operation succeeded
% results: the results of the operation in the terminal
%
% Notes
% -----
% http://www.brain.org.au/software/mrtrix/tractography/preprocess.html

if notDefined('verbose')
   verbose = true;
end
if notDefined('b_file')
    b_file = 'encoding.b';
end

% This command generates  tensors:
cmd_str = sprintf('dwi2tensor %s %s -grad %s',in_file, outfile, b_file);

% Send it to mrtrix:
[status,results] = mrtrix_cmd(cmd_str,verbose);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function drawAcrobot(t,z,p)

clf; hold on;

length = p.l1+p.l2;
axis equal; axis(length*[-1,1,-1,1]); axis off;

[p1,p2] = acrobotKinematics(z,p);
pos = [[0;0],p1,p2];

plot(0,0,'ks','MarkerSize',25,'LineWidth',4)
plot(pos(1,:),pos(2,:),'Color',[0.1, 0.8, 0.1],'LineWidth',4)
plot(pos(1,:),pos(2,:),'k.','MarkerSize',50)

title(sprintf('Acrobot Animation,  t = %6.4f', t));

drawnow; pause(0.001);

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function pass = test_eigs()
% TAD, 10 Jan 2014

tol = 1e-7;

dom = [-pi/2, pi/2];
D2 = operatorBlock.diff(dom, 2);
E = functionalBlock.eval(dom);
El = E(dom(1));
Er = E(dom(end));
L = linop(D2);
L = addbc(L, El, 0);
L = addbc(L, Er, 0);

e_true = flipud(-(1:6).'.^2);

%%
prefs = cheboppref;
prefs.discretization = @chebcolloc2;
[V, D] = eigs(L, 6, prefs);
e = diag(D);
err(1) = norm(e - e_true, inf);
% Check that we actually computed eigenfunctions
err(2) = norm(L*V-V*D);
%%
prefs.discretization = @ultraS;
[V, D] = eigs(L, 6, 0, prefs);
e = diag(D);
err(3) = norm(e - e_true, inf);
err(4) = norm(L*V-V*D);
%%
prefs.discretization = @chebcolloc1;
[V, D] = eigs(L, 6, prefs);
e = diag(D);
err(5) = norm(e - e_true, inf);
err(6) = norm(L*V-V*D);
%%
pass = err < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [nKeq, nGC, nNone] = standardGibbsFormationEnergyStats(modelT, figures)
% Generate the stats +/- pie chart on the provinence of the metabolite standard Gibbs energies.
%
% USAGE:
%
%    [nKeq, nGC, nNone] = standardGibbsFormationEnergyStats(modelT, figures)
%
% INPUTS:
%    modelT:
%    figures:
%
% OUTPUTS:
%    nKeq:
%    nGC:
%    nNone:
%
% .. Author: - Ronan M.T. Fleming

[nMet,nRxn]=size(modelT.S);

nKeq=eps;
nGC=eps;
nNone=eps;

noneBool=false(nMet,1);
p=1;
noneAbbr=[];
for m=1:nMet
    if strcmp(modelT.met(m).dGft0Source,'Keq')
        nKeq=nKeq+1;
    else
        if isnan(modelT.met(m).dGft0Source)
            nNone=nNone+1;
            noneBool(m,1)=1;
            abbr=modelT.mets{m};
            abbr=abbr(1:end-3);
            noneAbbr{p,1}=abbr;
            p=p+1;
        else
            nGC=nGC+1;
        end
    end
end

if figures
    figure1=figure;
    h=pie([nKeq,nGC,nNone],{{'K_{eq}',int2str(nKeq)},{'Group Contribution',int2str(nGC)},{'None',int2str(nNone)}});
    textObjs = findobj(h,'Type','text');
    set(textObjs,'FontSize',16);
    title('Provinence of reactant \Delta_{f}G^{\o} data','FontSize',16);
    saveas(figure1,'dfG0provinencePieChart','fig');

end

fprintf('%s\n',['Number of unique reactants without dGf0: ' int2str(length(unique(noneAbbr))) ]);

author: opencobra

repo: cobratoolbox


text:
function test_suite = test_angle3Points
% One-line description here, please.
%   output = testAngle3Points(input)
%
%   Example
%   testAngle3Points
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-04-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testSimple(testCase) %#ok<*DEFNU>
% all points inside window, possibly touching edges

p1 = [10 0];
p2 = [0 0];
p3 = [0 10];
angle = angle3Points(p1, p2, p3);
testCase.assertEqual(pi/2, angle, 'AbsTol', .01);

function testBundledInput(testCase)
% all points inside window, possibly touching edges

p1 = [10 0];
p2 = [0 0];
p3 = [0 10];
angle = angle3Points([p1; p2; p3]);
testCase.assertEqual(pi/2, angle, 'AbsTol', .01);

function testArray(testCase)
% all points inside window, possibly touching edges

p1 = [10 0; 20 0];
p2 = [0 0;0 0];
p3 = [0 10; 0 20];
angle = angle3Points(p1, p2, p3);

testCase.assertEqual(2, size(angle, 1));
testCase.assertEqual([pi/2;pi/2], angle, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%  Figure 10.70      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% fig10_87.m is a script to generate Fig. 10.70 the linear RTP response
% LQG method with internal model
% RTP chamber demo Example
clf;
a=[-0.068208813989939 0.014929776357245 0.000000065782442;...
   0.045809672136430 -0.118134773528570 0.021802006306129;...
   0.000000433637498 0.046839194867347 -0.100884399149872];
b3=[0.37873508304055 0.110575403544586 0.022912893467038;...
   0.000000002974222 0.449046982554739 0.073572900808161;...
   0.000000027324116 0.000702342292121 0.417797228783230];
c3=eye(3,3);
d=0*eye(3,3);
sysp=ss(a,b3,c3,d);
[eol]=eig(a);
[zol]=tzero(sysp);
[zol22]=tzero(a,b3(:,2),c3(2,:),d(2,2));
% Combine 3 lamps into a single actuator
b=b3(:,1)+b3(:,2)+b3(:,3);
% Select center temperature
c=c3(2,:)
aa=[zeros(1,1) c;zeros(3,1) a];
bb=[zeros(1,1);b];
%
% weight temperature differences
q1hat=[eye(1,1) zeros(1,3);0 2 -1 -1; 0 -1 2 -1; 0 -1 -1 2]+1e-6*eye(4,4);
q1hat=diag([1 10 10 10])*q1hat;
q2hat=eye(1,1);
[kk,sric,ee]=lqr(aa,bb,q1hat,q2hat);
k1=kk(1,1);
ko=kk(:,2:4);
acl=[a-b*ko b;-k1*c zeros(1,1)];
bcl=[zeros(3,1);k1];
ccl=[c zeros(1,1)];
dcl=zeros(1,1);
[ecl]=eig(acl);
[zcl]=tzero(acl,bcl,ccl,dcl);
%CL DC gain
[cldcgain]=dcl-ccl*inv(acl)*bcl;
%CL Step Response
%[y,t]=step(acl,bcl,ccl,dcl);
%s1y=plot(y(:,:,1));
%grid;
%pause;
%Control effort
%cclu=[-ko eye(1,1)];
%dclu=[0];
%[uu,t]=step(acl,bcl,cclu,dclu);
%su=plot(uu(:,:,1));
%grid;
%pause
%Step in all channels
%t=0:.1:100;
%u=[25*ones(1,251) 25*ones(1,500)0*ones(1,250)];
%u10=[u'];
%sysc1=ss(acl,bcl,ccl,dcl);
%[yy1,t]=lsim(sysc1,u10,t);
%plot(t,yy1)
%grid
%pause;
%cclu=[-ko eye(1,1)];
%dclu=[0];
%syscl2=ss(acl,bcl,cclu,dclu);
%[uu1,t]=lsim(syscl2,u10,t);
%su=plot(uu1(:,:,1));
%grid;
%pause
% Estimator design
qe=eye(1,1);
re=0.001*eye(1,1);
[ll,pp,el]=lqe(a,b,c,qe,re);
ac=zeros(1,1);
bc=-k1;
cc=eye(1,1);
dc=-ko;
acle=[a, b*cc, -b*ko;
   bc*c, ac, zeros(1,3);
   ll*c, b*cc, a-ll*c-b*ko];
bcle=[zeros(3,1);-bc;zeros(3,1)];
ccle=[c, zeros(1,4)];
dcle=zeros(1,1);
[ecle]=eig(acle);
[zcle]=tzero(acle,bcle,ccle,dcle);
dcgain=dcle-ccle*inv(acle)*bcle;
t=0:.1:100;
R=[0:.1:25, 25*ones(1,500), 0*ones(1,250)];
R11=[R'];
syscl=ss(acle,bcle,ccle,dcle);
[yy,t]=lsim(syscl,R11,t);
plot(t,yy,'--');
grid on;
hold on;
plot(t,R11,'-');
xlabel('Time (sec)');
ylabel('Temperature (K)');
title('Fig. 10.70 (a) Internal model controller: temperature tracking response');
%legend('y')
pause;
hold off;
cclu=[zeros(1,3), eye(1,1), -ko];
dclu=zeros(1,1);
syscu=ss(acle,bcle,cclu,dclu);
[uuu,t]=lsim(syscu,R11,t);
plot(t(1:753,:),uuu(1:753,:));
hold on;
plot(t(752:818,:),0*ones(67,3));
hold on;
plot(t(753:818,:),uuu(753:818,:),'--');
hold on;
plot(t(819:1001,:),uuu(819:1001,:));
xlabel('Time (sec)');
ylabel('Lamp voltage');
title('Fig. 10.70 (b) Internal model controller: control effort');
legend('u')
grid;
hold off;

author: Sable

repo: mcbench-benchmarks


text:
function [ grid_weight, grid_point ] = sparse_grid_cfn ( dim_num, level_max, ...
  rule, point_num )

%*****************************************************************************80
%
%% SPARSE_GRID_CFN computes a sparse grid based on a CFN 1D rule.
%
%  Discussion:
%
%    The 1D quadrature rule is assumed to be Closed Fully Nested.
%
%    Closed Fully Nested rules include Clenshaw Curtis rules.
%
%    A Smolyak construction is used to create a multidimensional sparse grid.
%
%    The user specifies:
%    * the spatial dimension of the quadrature region,
%    * the level that defines the Smolyak grid.
%    * the quadrature rule.
%    * the number of points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, controls the size of the final
%    sparse grid.
%
%    Input, integer RULE, the index of the rule.
%    1, "CC", Clenshaw Curtis Closed Fully Nested rule.
%    2, "F1", Fejer 1 Open Fully Nested rule.
%    3, "F2", Fejer 2 Open Fully Nested rule.
%    4, "GP", Gauss Patterson Open Fully Nested rule.
%    5, "GL", Gauss Legendre Open Weakly Nested rule.
%    6, "GH", Gauss Hermite Open Weakly Nested rule.
%    7, "LG", Gauss Laguerre Open Non Nested rule.
%
%    Input, integer POINT_NUM, the number of points in the grid,
%    as determined by LEVELS_INDEX_SIZE.
%
%    Output, real GRID_WEIGHT(POINT_NUM), the weights.
%
%    Output, real GRID_POINT(DIM_NUM,POINT_NUM), the points.
%
  if ( rule ~= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SPARSE_GRID_CFN - Fatal error!\n' );
    fprintf ( 1, '  Illegal input rule index = %d\n', rule );
    error ( 'SPARSE_GRID_CFN - Fatal error!' );
  end
%
%  Determine the index vector, relative to the full product grid,
%  that identifies the points in the sparse grid.
%
  [ grid_index, grid_base ] = levels_index_cfn ( dim_num, level_max, point_num );
%
%  Compute the physical coordinates of the abscissas.
%
  if ( 0 == level_max )
    order_max = 1;
  else
    order_max = 2^level_max + 1;
  end

  grid_point = zeros ( dim_num, point_num );

  for point = 1 : point_num
    for dim = 1 : dim_num

      if ( rule == 1 )
        grid_point(dim,point) = ...
          cc_abscissa ( order_max, grid_index(dim,point) + 1 );
      end

    end
  end
%
%  Gather the weights.
%
  grid_weight = sparse_grid_weights_cfn ( dim_num, level_max, rule, ...
    point_num, grid_index );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function hRD = plugin_rd_rx
% Reference design definition

% Call the common reference design definition function
hRD = AnalogDevices.adrv9361z7035.common.plugin_rd('ccbob_lvds','Rx');
AnalogDevices.adrv9361z7035.common.add_io(hRD, 'Rx', 'ccbob_lvds');

author: analogdevicesinc


text:
function [z, J, H] = fun(x)

[x1,x2,x3] = split(x);

z = [...
    x1*x2+sin(x3);
    x2/x1-cos(x3)];

if nargout > 1
    J = [...
              x2   x1 cos(x3)
        -x2/x1^2 1/x1 sin(x3)];

    if nargout > 2
        H(:,:,1) = [...
                    0       1 0
            2*x2/x1^3 -1/x1^2 0];

        H(:,:,2) = [...
                  1 0 0
            -1/x1^2 0 0];

        H(:,:,3) = [...
            0 0 -sin(x3)
            0 0  cos(x3)];
    end
end
return

%% jac
syms x1 x2 x3 real
x = [x1;x2;x3];

z = fun(x);

J = jacobian(z,x)
H1 = diff(J,'x1')
H2 = diff(J,'x2')
H3 = diff(J,'x3')

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
% eeg_mergelocs() - merge channel structure while preserving channel
%                   order
%
%      >> mergedlocs = eeg_mergelocs(loc1, loc2, loc3, ...);
%
% Inputs:
%     loc1     - EEGLAB channel location structure
%     loc2     - second EEGLAB channel location structure
%
% Output:
%     mergedlocs - merged channel location structure
%
% Author: Arnaud Delorme, August 2006

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

function alllocs = eeg_mergelocs(varargin)

persistent warning_shown;

% sort by length
% --------------
len = cellfun(@length, varargin);
[tmp so] = sort(len, 2, 'descend');
varargin = varargin(so);

alllocs = varargin{1};
for index = 2:length(varargin)

    % fuse while preserving order (assumes the same channel order)
    % ------------------------------------------------------------
    tmplocs = varargin{index};
    newlocs = myunion(alllocs, tmplocs);

    if length(newlocs) > length(union({ alllocs.labels }, { tmplocs.labels }))

        if isempty(warning_shown)
            disp('Warning: different channel montage order for the different datasets');
            warning_shown = 1;
        end;

        % trying to preserve order of the longest array
        %----------------------------------------------
        if length(alllocs) < length(tmplocs)
            tmp     = alllocs;
            alllocs = tmplocs;
            tmplocs = tmp;
        end;
        allchans = { alllocs.labels tmplocs.labels };
        [uniquechan ord1 ord2 ]  = unique( allchans );

        [tmp rminds] = intersect( uniquechan, { alllocs.labels });
        ord1(rminds) = [];
        tmplocsind = ord1-length(alllocs);

        newlocs = concatlocs(alllocs, tmplocs(tmplocsind));

    end;
    alllocs = newlocs;
end;

% union of two channel location structure
% without loosing the order information
% ---------------------------------------
function alllocs = myunion(locs1, locs2)

   labs1 = { locs1.labels };
   labs2 = { locs2.labels };

   count1 = 1;
   count2 = 1;
   count3 = 1;
   alllocs = locs1; alllocs(:) = [];
   while count1 <= length(locs1) | count2 <= length(locs2)

       if count1 > length(locs1)
           alllocs = copyfields(alllocs, count3, locs2(count2));
           %alllocs(count3) = locs2(count2);
           count2 = count2 + 1;
           count3 = count3 + 1;
       elseif count2 > length(locs2)
           alllocs = copyfields(alllocs, count3, locs1(count1));
           %alllocs(count3) = locs1(count1);
           count1 = count1 + 1;
           count3 = count3 + 1;
       elseif strcmpi(labs1{count1}, labs2{count2})
           alllocs = copyfields(alllocs, count3, locs1(count1));
           %alllocs(count3) = locs1(count1);
           count1 = count1 + 1;
           count2 = count2 + 1;
           count3 = count3 + 1;
       elseif isempty(strmatch(labs1{count1}, labs2, 'exact'))
           alllocs = copyfields(alllocs, count3, locs1(count1));
           %alllocs(count3) = locs1(count1);
           count1 = count1 + 1;
           count3 = count3 + 1;
       else
           alllocs = copyfields(alllocs, count3, locs2(count2));
           %alllocs(count3) = locs2(count2);
           count2 = count2 + 1;
           count3 = count3 + 1;
       end;

   end;

% concatenate channel structures with different fields
function loc3 = concatlocs(loc1, loc2);

    fields1 = fieldnames(loc1);
    fields2 = fieldnames(loc2);
    if isequal(fields1, fields2)
        loc3 = [ loc1 loc2 ];
    else
        loc3 = loc1;
        for index = 1:length(loc2)
            loc3 = copyfields(loc3, length(loc1)+index, loc2(index));
        end;
    end;

% copy fields of a structure to another one
function [struct1] = copyfields(struct1, index1, struct2)

    fields1 = fieldnames(struct1);
    fields2 = fieldnames(struct2);

    if isequal(fields1, fields2)
        struct1(index1) = struct2;
    else
        for index = 1:length(fields2)
            struct1 = setfield(struct1, {index1}, fields2{index}, getfield(struct2, fields2{index}));
        end;
    end;

author: PatternRecognition

repo: OpenBMI


text:
function TorF = hasPQcap(gen, hilo)
%HASPQCAP  Checks for P-Q capability curve constraints.
%   TORF = HASPQCAP(GEN, HILO) returns a column vector of 1's and 0's. The 1's
%   correspond to rows of the GEN matrix which correspond to generators which
%   have defined a capability curve (with sloped upper and/or lower bound on
%   Q) and require that additional linear constraints be added to the OPF.
%
%   The GEN matrix in version 2 of the MATPOWER case format includes columns
%   for specifying a P-Q capability curve for a generator defined as the
%   intersection of two half-planes and the box constraints on P and Q. The
%   two half planes are defined respectively as the area below the line
%   connecting (Pc1, Qc1max) and (Pc2, Qc2max) and the area above the line
%   connecting (Pc1, Qc1min) and (Pc2, Qc2min).
%
%   If the optional 2nd argument is 'U' this function returns true only for
%   rows corresponding to generators that require the upper constraint on Q.
%   If it is 'L', only for those requiring the lower constraint. If the 2nd
%   argument is not specified or has any other value it returns true for rows
%   corresponding to gens that require either or both of the constraints.
%
%   It is smart enough to return true only if the corresponding linear
%   constraint is not redundant w.r.t the box constraints.

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

[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;

%% default value
if nargin < 2
    hilo = 'B';     %% look at both top and bottom by default
end

%% for which gens is it specified
k = find( gen(:, PC1) | gen(:, PC2) );
ng = size(gen, 1);

if isempty(k)
    TorF = zeros(ng, 1);
else
    %% eliminate cases where QMIN = QMAX = QC
    kk = find( gen(k, QMIN) == gen(k, QMAX) & ...
                gen(k, QMIN) == gen(k, QC1MAX) & ...
                gen(k, QMIN) == gen(k, QC1MIN) & ...
                gen(k, QMIN) == gen(k, QC2MAX) & ...
                gen(k, QMIN) == gen(k, QC2MIN) );
    k(kk) = [];

    %% check for errors in capability curve data
    if any( gen(k, PC1) >= gen(k, PC2) )
        error('hasPQcap: must have Pc1 < Pc2');
    end
    if any( gen(k, QC2MAX) <= gen(k, QC2MIN) & gen(k, QC1MAX) <= gen(k, QC1MIN) )
        error('hasPQcap: capability curve defines an empty set');
    end

    %% for which gens is it specified
    k = find( gen(:, PC1) ~= gen(:, PC2) );
    L = zeros(ng, 1);
    U = zeros(ng, 1);
    dPc = gen(k, PC2) - gen(k, PC1);

    if ~strcmp(hilo, 'U')       %% include lower constraint
        dQc = gen(k, QC2MIN) - gen(k, QC1MIN);
        Qmin_at_Pmin = gen(k, QC1MIN) + (gen(k, PMIN) - gen(k, PC1)) .* ...
            dQc ./ dPc;
        Qmin_at_Pmax = gen(k, QC1MIN) + (gen(k, PMAX) - gen(k, PC1)) .* ...
            dQc ./ dPc;
        L(k) = Qmin_at_Pmin > gen(k, QMIN) | Qmin_at_Pmax > gen(k, QMIN);
    end

    if ~strcmp(hilo, 'L')       %% include upper constraint
        dQc = gen(k, QC2MAX) - gen(k, QC1MAX);
        Qmax_at_Pmin = gen(k, QC1MAX) + (gen(k, PMIN) - gen(k, PC1)) .* ...
            dQc ./ dPc;
        Qmax_at_Pmax = gen(k, QC1MAX) + (gen(k, PMAX) - gen(k, PC1)) .* ...
            dQc ./ dPc;
        U(k) = Qmax_at_Pmin < gen(k, QMAX) | Qmax_at_Pmax < gen(k, QMAX);
    end

    TorF = L | U;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [v] = spm_mar_gen (w,A,C,n,ndisc)
% Generate data from MAR model
% FORMAT [v] = spm_mar_gen (w,A,C,n,ndisc)
%
% Generates n time steps of the MAR(p) process
%
%     v(k,:)' = w' + A1*v(k-1,:)' +...+ Ap*v(k-p,:)' + eta(k,:)',
%
%  where A=[A1 ... Ap] is the coefficient matrix, and w is a vector of
%  intercept terms that is included to allow for a nonzero mean of the
%  process. The vectors eta(k,:) are independent Gaussian noise
%  vectors with mean zero and covariance matrix C.
%
%  This function is adapted from the ARFIT toolbox by Neumaier and
%  Schneider
%___________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mar_gen.m 1143 2008-02-07 19:33:33Z spm $

m       = size(C,1);                  % dimension of state vectors
p       = size(A,2)/m;                % order of process

if (p ~= round(p))
    error('Bad arguments.');
end

if (length(w) ~= m | min(size(w)) ~= 1)
    error('Dimensions of arguments are mutually incompatible.')
end
w       = w(:)';                      % force w to be row vector

% Discard the first ndisc time steps; if ndisc is not given as input
% argument, use default
if (nargin < 5)
    ndisc = 10^3;
end

% Compute Cholesky factor of covariance matrix C
R       = chol(C);                    % R is upper triangular

% Get ndisc+n independent Gaussian pseudo-random vectors with
% covariance matrix C=R'*R
randvec = randn([ndisc+n,m])*R;

% Add intercept vector to random vectors
randvec = randvec + ones(ndisc+n,1)*w;

% Get transpose of system matrix A (use transpose in simulation because
% we want to obtain the states as row vectors)
AT      = A';

% Take the p initial values of the simulation to equal the process mean,
% which is calculated from the parameters A and w
if any(w)
    %  Process has nonzero mean    mval = inv(B)*w'    where
    %             B = eye(m) - A1 -... - Ap;
    %  Assemble B
    B    = eye(m);
    for j=1:p
        B = B - A(:, (j-1)*m+1:j*m);
    end
    %  Get mean value of process
    mval = w / B';

    %  The optimal forecast of the next state given the p previous
    %  states is stored in the vector x. The vector x is initialized
    %  with the process mean.
    x    = ones(p,1)*mval;
else
    %  Process has zero mean
    x    = zeros(p,m);
end

% Initialize state vectors
u      = [x; zeros(ndisc+n,m)];

% Simulate n+ndisc observations. In order to make use of Matlab's
% vectorization capabilities, the cases p=1 and p>1 must be treated
% separately.
if p==1
    for k=2:ndisc+n+1;
        x(1,:) = u(k-1,:)*AT;
        u(k,:) = x + randvec(k-1,:);
    end
else
    for k=p+1:ndisc+n+p;
        for j=1:p;
            x(j,:) = u(k-j,:)*AT((j-1)*m+1:j*m,:);
        end
        u(k,:) = sum(x)+randvec(k-p,:);
    end
end

% return only the last n simulated state vectors
v = u(ndisc+p+1:ndisc+n+p,:);

author: spm


text:
function [ modulated ] = WGFT_gmod(g,k,V)

% Modulate a signal g on a graph whose eigenvectors are the columns of V
% kernel g is input in the vertex domain
% k is in 1:N

N=size(V,1);
dc=V(:,1);
factor=1./dc;
modulated=factor.*V(:,k).*g;

% modulated=sqrt(N)*V(:,k).*g; OLD CODE - ADDED FACTOR TO ADJUST FOR
% NORMALIZED
end

repo: gspbox


text:
function val=erfComplex(z)
%%ERFCOMPLEX Evaluate the error function in a manner valid for complex
%            arguemnts. Matlab's built-in function erf cannot handle
%            complex arguments. The error function is defined as
%            (2/pi)*integral_0^z exp(-t^2) dt.
%
%INPUTS: z A scalar, vector, or matrix of values at which one wishes to
%          evaluate the error function.
%
%OUTPUTS: val A matrix having the same dimensions as z in which the values
%             of the scaled error function are computed for the values in
%             z.
%
%For values of abs(z)<0.08, six terms of the maclaurin series in [1] are
%used. These should be sufficient for convergence. For larger values, the
%Faddeeva function is used through the relation
%erf(z)=1-exp(-z^2)*Faddeeva(1i*z);
%
%REFERENCES:
%[1] Weisstein, Eric W. "Erf." From MathWorld--A Wolfram Web Resource
%    http://mathworld.wolfram.com/Erf.html
%
%November 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numVals=numel(z);

zList=z;
val=zeros(size(z));

for curVal=1:numVals
    z=zList(curVal);
    y=imag(z);
    %If the argument is real, the built-in erf function is faster.
    if(y==0)
        x=real(z);
        val(curVal)=erf(x);
        continue;
    end

    if(abs(z)<0.08)
        %Use the Maclaurin series. For abs(z)<0.08, convergence occurs
        %within 6 terms.
        NMax=6;
        z2=z^2;

        %Initialize to the n=0 term
        sumVal=z;
        prodVal=z;
        for n=1:NMax
            prodVal=prodVal*z2*(1-2*n)/(n*(1+2*n));
            sumVal=sumVal+prodVal;
        end

        val(curVal)=2/sqrt(pi)*sumVal;
    else%Compute using the Faddeeva function.
        val(curVal)=1-exp(-z^2)*Faddeeva(1i*z);
    end
end
end

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


text:
function fem2d_poisson_rectangle_linear ( nx, ny )

%*****************************************************************************80
%
%% MAIN is the main routine for FEM2D_POISSON_RECTANGLE_LINEAR.
%
%  Discussion:
%
%    This program solves
%
%      - d2U(X,Y)/dx2 - d2U(X,Y)/dy2 = F(X,Y)
%
%    in a rectangular region in the plane.
%
%    Along the boundary of the region, Dirichlet conditions
%    are imposed:
%
%      U(X,Y) = G(X,Y)
%
%    The code uses continuous piecewise linear basis functions on
%    triangles determined by a uniform grid of NX by NY points.
%
%    u    =      sin ( pi * x ) * sin ( pi * y ) + x
%
%    dudx = pi * cos ( pi * x ) * sin ( pi * y ) + 1
%    dudy = pi * sin ( pi * x ) * cos ( pi * y )
%
%    d2udx2 = - pi * pi * sin ( pi * x ) * sin ( pi * y )
%    d2udy2 = - pi * pi * sin ( pi * x ) * sin ( pi * y )
%
%    rhs  = 2 * pi * pi * sin ( pi * x ) * sin ( pi * y )
%
%  THINGS YOU CAN EASILY CHANGE:
%
%    1) Change NX or NY, the number of nodes in the X and Y directions.
%    2) Change XL, XR, YB, YT, the left, right, bottom and top limits
%       of the rectangle.
%    3) Change the exact solution in the EXACT routine, but make sure you also
%       modify the formula for RHS in the assembly portion of the program%
%
%  HARDER TO CHANGE:
%
%    4) Change from "linear" to "quadratic" triangles;
%    5) Change the region from a rectangle to a general triangulated region;
%    6) Handle Neumann boundary conditions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NX, NY, the number of nodes in the X and Y directions.
%
%  Local Parameters:
%
%    Local, sparse real A(NODE_NUM,NODE_NUM), the finite element system matrix.
%
%    Local, real B(NODE_NUM), the finite element right hand side.
%
%    Local, real C(NODE_NUM), the finite element coefficient vector.
%
%    Local, integer ELEMENT_NODE(3,ELEMENT_NUM), the indices of the nodes
%    that form each element.
%
%    Local, integer ELEMENT_NUM, the number of elements.
%
%    Local, integer NODE_NUM, the number of nodes.
%
%    Local, real NODE_XY(2,NODE_NUM), the X and Y coordinates of each node.
%
%    Local, real XL, the X coordinate of the left boundary.
%
%    Local, real XR, the X coordinate of the right boundary.
%
%    Local, real YB, the Y coordindate of the bottom boundary.
%
%    Local, real YT, the Y coordinate of the top boundary.
%
  timestamp ( );

  element_order = 3;

  xl = 0.0;
  xr = 1.0;
  yb = 0.0;
  yt = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_POISSON_RECTANGLE_LINEAR\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solution of the Poisson equation:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  - Uxx - Uyy = F(x,y) inside the region,\n' );
  fprintf ( 1, '       U(x,y) = G(x,y) on the boundary of the region.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The region is a rectangle, defined by:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ' %f = XL<= X <= XR = %f\n', xl, xr );
  fprintf ( 1, ' %f = YB<= Y <= YT = %f\n', yb, yt );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The finite element method is used, with piecewise\n' );
  fprintf ( 1, '  linear basis functions on 3 node triangular\n' );
  fprintf ( 1, '  elements.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The corner nodes of the triangles are generated by an\n' );
  fprintf ( 1, '  underlying grid whose dimensions are\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  NX =                       %d\n', nx );
  fprintf ( 1, '  NY =                       %d\n', ny );
%
%  NODE COORDINATES
%
%  Numbering of nodes is suggested by the following 5x10 example:
%
%    J=5 | K=41  K=42 ... K=50
%    ... |
%    J=2 | K=11  K=12 ... K=20
%    J=1 | K= 1  K= 2     K=10
%        +--------------------
%          I= 1  I= 2 ... I=10
%
  node_num = nx * ny;

  fprintf ( 1, '  Number of nodes =          %d\n', node_num );

  node_xy = zeros(2,node_num);

  k = 0;
  for j = 1 : ny
    for i = 1 : nx

      k = k + 1;

      node_xy(1,k) = ( ( nx - i     ) * xl   ...
                     + (      i - 1 ) * xr ) ...
                     / ( nx     - 1 );

      node_xy(2,k) = ( ( ny - j     ) * yb   ...
                     + (      j - 1 ) * yt ) ...
                     / ( ny     - 1 );

    end
  end
%
%  ELEMENT array
%
%  Organize the nodes into a grid of 3-node triangles.
%  Here is part of the diagram for a 5x10 example:
%
%    |  \ |  \ |  \ |
%    |   \|   \|   \|
%   21---22---23---24--
%    |\ 8 |\10 |\12 |
%    | \  | \  | \  |
%    |  \ |  \ |  \ |  \ |
%    |  7\|  9\| 11\|   \|
%   11---12---13---14---15---16---17---18---19---20
%    |\ 2 |\ 4 |\ 6 |\  8|                   |\ 18|
%    | \  | \  | \  | \  |                   | \  |
%    |  \ |  \ |  \ |  \ |      ...          |  \ |
%    |  1\|  3\|  5\| 7 \|                   |17 \|
%    1----2----3----4----5----6----7----8----9---10
%
  element_num = 2 * ( nx - 1 ) * ( ny - 1 );

  fprintf ( 1, '  Number of elements =       %d\n', element_num );

  element_node = zeros ( element_order, element_num );

  k = 0;

  for j = 1 : ny - 1
    for i = 1 : nx - 1
%
%     (I,J+1)-
%      |  \
%      |   \
%      |    \  |
%    (I,J)---(I+1,J)
%
      k = k + 1;
      element_node(1,k) = i     + ( j - 1 ) * nx;
      element_node(2,k) = i + 1 + ( j - 1 ) * nx;
      element_node(3,k) = i     +   j       * nx;
%
%    (I,J+1)--(I+1,J+1)
%      |  \    |
%          \   |
%           \  |
%         -(I+1,J)
%
      k = k + 1;
      element_node(1,k) = i + 1 +   j       * nx;
      element_node(2,k) = i     +   j       * nx;
      element_node(3,k) = i + 1 + ( j - 1 ) * nx;

    end
  end
%
%  ASSEMBLE THE SYSTEM
%
%  Assemble the coefficient matrix A and the right-hand side B of the
%  finite element equations, ignoring boundary conditions.
%
  b = zeros(node_num,1);
  a = sparse ( [], [], [], node_num, node_num );

  for e = 1 : element_num

    i1 = element_node(1,e);
    i2 = element_node(2,e);
    i3 = element_node(3,e);

    area = 0.5 * ...
      ( node_xy(1,i1) * ( node_xy(2,i2) - node_xy(2,i3) ) ...
      + node_xy(1,i2) * ( node_xy(2,i3) - node_xy(2,i1) ) ...
      + node_xy(1,i3) * ( node_xy(2,i1) - node_xy(2,i2) ) );
%
%  Consider each quadrature point.
%  Here, we use the midside nodes as quadrature points.
%
    for q1 = 1 : 3

      q2 = mod ( q1, 3 ) + 1;

      nq1 = element_node(q1,e);
      nq2 = element_node(q2,e);

      xq = 0.5 * ( node_xy(1,nq1) + node_xy(1,nq2) );
      yq = 0.5 * ( node_xy(2,nq1) + node_xy(2,nq2) );
      wq = 1.0 / 3.0;
%
%  Consider each test function in the element.
%
      for ti1 = 1 : element_order

        ti2 = mod ( ti1,     3 ) + 1;
        ti3 = mod ( ti1 + 1, 3 ) + 1;

        nti1 = element_node(ti1,e);
        nti2 = element_node(ti2,e);
        nti3 = element_node(ti3,e);

        qi = 0.5 * ( ...
            ( node_xy(1,nti3) - node_xy(1,nti2) ) * ( yq - node_xy(2,nti2) ) ...
          - ( node_xy(2,nti3) - node_xy(2,nti2) ) * ( xq - node_xy(1,nti2) ) ) ...
          / area;
        dqidx = - 0.5 * ( node_xy(2,nti3) - node_xy(2,nti2) ) / area;
        dqidy =   0.5 * ( node_xy(1,nti3) - node_xy(1,nti2) ) / area;

        rhs = 2.0 * pi * pi * sin ( pi * xq ) * sin ( pi * yq );

        b(nti1) = b(nti1) + area * wq * rhs * qi;
%
%  Consider each basis function in the element.
%
        for tj1 = 1 : element_order

          tj2 = mod ( tj1,     3 ) + 1;
          tj3 = mod ( tj1 + 1, 3 ) + 1;

          ntj1 = element_node(tj1,e);
          ntj2 = element_node(tj2,e);
          ntj3 = element_node(tj3,e);

          qj = 0.5 * ( ...
              ( node_xy(1,ntj3) - node_xy(1,ntj2) ) * ( yq - node_xy(2,ntj2) ) ...
            - ( node_xy(2,ntj3) - node_xy(2,ntj2) ) * ( xq - node_xy(1,ntj2) ) ) ...
              / area;
          dqjdx = - 0.5 * ( node_xy(2,ntj3) - node_xy(2,ntj2) ) / area;
          dqjdy =   0.5 * ( node_xy(1,ntj3) - node_xy(1,ntj2) ) / area;

          a(nti1,ntj1) = a(nti1,ntj1) ...
            + area * wq * ( dqidx * dqjdx + dqidy * dqjdy );

        end

      end

    end

  end
%
%  BOUNDARY CONDITIONS
%
%  If the K-th variable is at a boundary node, replace the K-th finite
%  element equation by a boundary condition that sets the variable to U(K).
%
  k = 0;

  for j = 1 : ny

    for i = 1 : nx

      k = k + 1;

      if ( i == 1 | i == nx | j == 1 | j == ny )

        [ u, dudx, dudy ] = exact ( node_xy(1,k), node_xy(2,k) );

        a(k,1:node_num) = 0.0;
        a(k,k)          = 1.0;
        b(k)            = u;

      end
    end
  end
%
%  SOLVE the linear system A * C = B.
%
  c = a \ b;
%
%  COMPARE computed and exact solutions at the nodes.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '     K     I     J          X           Y        U               U\n' );
  fprintf ( 1, '                                                 exact           computed \n' );

  k = 0;

  for j = 1 : ny
    fprintf ( 1, '\n' );
    for i = 1 : nx

      k = k + 1;

      [ u, dudx, dudy ] = exact ( node_xy(1,k), node_xy(2,k) );

      fprintf ( 1, '  %4d  %4d  %4d  %10f  %10f  %14e  %14e  %14e\n', ...
        k, i, j, node_xy(1,k), node_xy(2,k), u, c(k), abs ( u - c(k) ) );

    end

  end
%
%  Compute error integrals.
%
  if ( 1 )

  el2 = 0.0;
  eh1 = 0.0;

  for e = 1 : element_num

    i1 = element_node(1,e);
    i2 = element_node(2,e);
    i3 = element_node(3,e);

    area = 0.5 * ...
      ( node_xy(1,i1) * ( node_xy(2,i2) - node_xy(2,i3) ) ...
      + node_xy(1,i2) * ( node_xy(2,i3) - node_xy(2,i1) ) ...
      + node_xy(1,i3) * ( node_xy(2,i1) - node_xy(2,i2) ) );
%
%  Consider each quadrature point.
%  Here, we use the midside nodes as quadrature points.
%
    for q1 = 1 : 3

      q2 = mod ( q1, 3 ) + 1;

      nq1 = element_node(q1,e);
      nq2 = element_node(q2,e);

      xq = 0.5 * ( node_xy(1,nq1) + node_xy(1,nq2) );
      yq = 0.5 * ( node_xy(2,nq1) + node_xy(2,nq2) );
      wq = 1.0 / 3.0;

      uh = 0.0;
      dudxh = 0.0;
      dudyh = 0.0;

      for tj1 = 1 : element_order

        tj2 = mod ( tj1,     3 ) + 1;
        tj3 = mod ( tj1 + 1, 3 ) + 1;

        ntj1 = element_node(tj1,e);
        ntj2 = element_node(tj2,e);
        ntj3 = element_node(tj3,e);

        qj = 0.5 * ( ...
            ( node_xy(1,ntj3) - node_xy(1,ntj2) ) * ( yq - node_xy(2,ntj2) ) ...
          - ( node_xy(2,ntj3) - node_xy(2,ntj2) ) * ( xq - node_xy(1,ntj2) ) ) ...
              / area;
        dqjdx = - 0.5 * ( node_xy(2,ntj3) - node_xy(2,ntj2) ) / area;
        dqjdy =   0.5 * ( node_xy(1,ntj3) - node_xy(1,ntj2) ) / area;

        uh    = uh    + c(ntj1) * qj;
        dudxh = dudxh + c(ntj1) * dqjdx;
        dudyh = dudyh + c(ntj1) * dqjdy;

      end

      [ u, dudx, dudy ] = exact ( xq, yq );

      el2 = el2 + ( uh - u )^2 * area;
      eh1 = eh1 + ( ( dudxh - dudx )^2 + ( dudyh - dudy )^2 ) * area;

    end

  end

  el2 = sqrt ( el2 );
  eh1 = sqrt ( eh1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '*********************************************\n' );
  fprintf ( 1, '*                                           *\n' );
  fprintf ( 1, '*  ERRORS:                                  *\n' );
  fprintf ( 1, '*    L2 error =          %14f     *\n', el2 );
  fprintf ( 1, '*    H1-seminorm error = %14f     *\n', eh1 );
  fprintf ( 1, '*                                           *\n' );
  fprintf ( 1, '*********************************************\n' );

  end
%
%  WRITE the data to files.
%
  node_filename = 'rectangle_nodes.txt';

  r8mat_write ( node_filename, 2, node_num, node_xy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wrote the node file "%s"\n', node_filename );

  element_filename = 'rectangle_elements.txt';

  i4mat_write ( element_filename, element_order, element_num, element_node );

  fprintf ( 1, '  Wrote the element file "%s"\n', element_filename );

  value_filename = 'rectangle_solution.txt';

  r8mat_write ( value_filename, 1, node_num, c' );

  fprintf ( 1, '  Wrote the solution value file "%s"\n', value_filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_POISSON_RECTANGLE_LINEAR:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function [ u, dudx, dudy ] = exact ( x, y )

%*****************************************************************************80
%
%% EXACT calculates the exact solution and its first derivatives.
%
%  Discussion:
%
%    The function specified here depends on the problem being
%    solved.  The user must be sure to change both EXACT and RHS
%    or the program will have inconsistent data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the coordinates of a point
%    in the region, at which the exact solution is to be evaluated.
%
%    Output, real U, DUDX, DUDY, the value of
%    the exact solution U and its derivatives dUdX
%    and dUdY at the point (X,Y).
%
  u    =      sin ( pi * x ) * sin ( pi * y ) + x;
  dudx = pi * cos ( pi * x ) * sin ( pi * y ) + 1.0;
  dudy = pi * sin ( pi * x ) * cos ( pi * y );

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

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%

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

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

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

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

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

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function timestamp ( )

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

  return
end


text:
% demo_AAM

% This demo shows how to use the AAM and MESMA code, with a quick example
% on Gaussian artificial data sets.

% Initialization
d=40;         % number of spectral bands
num=1000;     % number of pixels in input data set
p=3;          % number of spectral libraries
N=[10,20,30]; % sizes of spectral libraries

% Generate data sets
L{1}=randn(d,N(1));
L{2}=randn(d,N(2));
L{3}=randn(d,N(3));
x=randn(d,num);

% Call AAM
[index1, abundances1, reconstruction1, error1]=AAM(x,L,3);

% Call MESMA
[index2, abundances2, reconstruction2, error2]=MESMA_bruteforce (x,L);

author: zhouyuanzxcv

repo: Hyperspectral


text:
% mvm_demo

% OVERVIEW: This script provides a demo on how to execute the MVM function,
% part of the MVToolbox in the Physionet-Cardiovascular-Signal-Toolbox.
% This demo script does the following,
% 1. Loads a sample ECG with annotations from the ./MVtoolbox/Demos/testdata subfolder
% 2. Remove power line interference at 60 Hz and baseline wander signal
% from the ecg
% 3. Perform arrhythmia detection on each of the segment_size min analysis windows,
% avoid performing MVM analysis on windows with positive arrythmia
% detection
% 4. Compute MVM on segment_size min windows not affected by arrhythmia
% 5. Store the output in the following variables,
%
%       variable                description (Units)
%       energyinband_array      QRS variability energy measured for each analysis window in the beatquency
%                               domain (beats^(-2))
%       sqi_array               The signal quality measured for each analysis
%                               window (unitless with value in the closed interval [0,1]. higher values mean cleaner signals)
%       heart_rate_estimate     The estimate of the median heart rate for each
%                               analysis window is given in this array (beats per minute)
%
% The expected values for each of the output variables is as follows,
%       energyinband_array = [8.6362e-08,NaN,1.1978e-07]
%       sqi_array = [0.9988,NaN,1]
%       heart_rate_estimate = [80,NaN,79.893]
%
% REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Ismail Sadiq
%	COPYRIGHT (C) 2019
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information. The license may be found in
%       the Documents folder of the Physionet-Cardiovascular-Signal-Toolbox.

% Load test ecg data
addpath(genpath('../../../../../PhysioNet-Cardiovascular-Signal-Toolbox-master/')); % add all dependencies from cardiovascular signal toolbox
sig = load('./testdata/mvm/sample_ecg_datam');
siginfo = readheader('./testdata/mvm/sample_ecg_datam.hea');
Fs = siginfo.freq;
ecg = (sig.val - siginfo.adczero)./siginfo.gain; % Adjust signal according to gain and dc offset
clear siginfo sig;

% Read QRSon, Q, R, S, QRSoff fiducial point annotations. Multiple
ann.QRSon = read_ann('./testdata/mvm/sample_ecgm','qrson'); ann.Q = read_ann('./testdata/mvm/sample_ecgm','q'); ann.R = read_ann('./testdata/mvm/sample_ecgm','r'); ann.S = read_ann('./testdata/mvm/sample_ecgm', 's'); ann.QRSoff = read_ann('./testdata/mvm/sample_ecgm', 'qrsoff');
%figure(1); plot(ecg); hold on; % Check annotation if needed
%stem(ann.QRSon, ecg(ann.QRSon)); stem(ann.Q, ecg(ann.Q)); stem(ann.R, ecg(ann.R)); stem(ann.S, ecg(ann.S)); stem(ann.QRSoff, ecg(ann.QRSoff)); hold off

% Read filter coefficients for power line interference (pli) filter
b = csvread('../Tools/MVM/pli_filter_coef_fs1000.csv');

% Add path to helper functions
addpath('../Tools/MVM/')

% Remove baseline wander
baseline = medianfilter_is(ecg, Fs);
ecg_baselinefree = ecg - baseline;

% Remove pli
% Filter each signal for pli
ecg_filtered = filter(b,1,ecg_baselinefree); clear ecg_baselinefree baseline
% Remove delay
delay = mean(grpdelay(b));
ecg_filtered(1:delay) = [];

% Detect arrhythmia

% Compute feature
minute_duration = 1*60*Fs;   % Duration of 1 minute in samples
segment_size = 5;   % Size of MVM analysis window in minutes
ii = 1; % Index to keep track of each non overlapping segment_size minute window

energyinband_array = NaN(1,ceil(length(ecg_filtered)/(segment_size*minute_duration))); % Array for storing the MVM energy evalauted for each segment_size minute analysis window, initialze to length of max possible non overlapping segment_size minute windows
sqi_array = energyinband_array; % Array for storing the signal quality for each analysis window, same dimensions as energyinband_array
heart_rate_est = energyinband_array; % Array for storing the heart rate estimate for each analysis window, same dimensions as energyinband_array
while (ii <= length(ecg_filtered))

    % Use atrial fibirillation arrhythmia detection in 60 beat non overlapping windows for each segment_size min analysis window
    if (ii+segment_size*minute_duration-1 > length(ecg_filtered))
        % Length of last analysis window is less than segment_size minutes, perform analysis till end of record
        endi = length(ecg_filtered);
    else
        endi = ii+segment_size*minute_duration-1;
    end

    % Get current 60 beat window for arrhythmia detection
    cur_win_r = ann.R(ann.R > ii); cur_win_r = cur_win_r(cur_win_r < endi);
    Index_test = NaN(1,ceil(length(cur_win_r)/60)); % Variable for storing arrhythmia detection for each 60 beat window
    for jj = 1:60:length(cur_win_r)
        % Get the current 60 beat window for arrhythmia detection
        if (jj+59 <= length(cur_win_r))
            cur_r = cur_win_r(jj:jj+59);
        else
            cur_r = cur_win_r(jj:end);
        end

        % If fewer than 12 beats or greater than 60 beats in the current
        % window do not perform arrhythmia detection
        if (length(diff(cur_r)) >= 12 && length(diff(cur_r)) <= 60)
            features = AF_features(diff(cur_r),Fs);
            Index_test(ceil(jj/60)) = SVM_AFdetection_withoutTrainingModel(features,1);
        else
            disp('Please input a RR interval time series with number of beats between 12 and 60')
        end

    end

    if (nansum(Index_test) == 0)
        % If no arrhythmia detection in current segment_size min window, compute MVM

        % Get qrson, q, r, s, qrsoff annotations for current window
        ann_win.R = cur_win_r-ii+1;
        ann_win.QRSon = ann.QRSon(ann.QRSon > ii); ann_win.QRSon = ann_win.QRSon(ann_win.QRSon < endi)-ii+1;
        ann_win.Q = ann.Q(ann.Q > ii); ann_win.Q = ann_win.Q(ann_win.Q < endi)-ii+1;
        ann_win.S = ann.S(ann.S > ii); ann_win.S = ann_win.S(ann_win.S < endi)-ii+1;
        ann_win.QRSoff = ann.QRSoff(ann.QRSoff > ii); ann_win.QRSoff = ann_win.QRSoff(ann_win.QRSoff < endi)-ii+1;

        % Get ecg in current window
        ecg_filtered_win = ecg_filtered(ii:endi);
        % Morph Var QRS
        disp(['Evaluating QRS MVM for analysis window ' num2str(ceil(ii/(segment_size*minute_duration))) ' starting at time (s) ' num2str(ii/Fs)])
        debug = 0;
        normalize = 1;
        [energyinband_array(ceil(ii/(segment_size*minute_duration))),sqi_array(ceil(ii/(segment_size*minute_duration))),heart_rate_est(ceil(ii/(segment_size*minute_duration)))] = ComputeMVM(ecg_filtered_win,ann_win,Fs,segment_size,normalize);
        % Analysis of current window successful, shift to next segment_size minute window
        ii = ii + segment_size*minute_duration;
    else
        % If arrhythmia detected in any consecutive 60 beats, shift analysis window by 1 minute.
        disp('Detected arrhythmia. shifting analysis window by 1 min.');
        ii = ii + minute_duration;
    end

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function  [ p_xy, p_type ] = state_initialize ( nb, ni, np, r );

%*****************************************************************************80
%
%% STATE_INITIALIZE initializes the state.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 December 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NB, the number of boundary points.
%
%    Input, integer NI, the number of interior points.
%
%    Input, integer NP, the total number of points.
%
%    Input, real R, the radius of the circle.
%
%    Output, real P_XY(2,NP), the coordinates of the points.
%
%    Output, integter P_TYPE(1,NP), the type of each point.
%    1, the point is constrained to the boundary.
%    2, the point is constrained to the interior.
%
  if ( 0 )
    p_xy = disk_sample_uniform ( np, r );
  else
    p_xy = disk_sample_nonuniform ( np, r );
  end

  p_type = zeros ( 1, np );
  p_type(1,1:nb)    = 1;
  p_type(1,nb+1:np) = 2;

  p_norm(1,1:np) = sqrt ( sum ( p_xy.^2, 1 ) );

  p_xy(1,1:nb) = r * p_xy(1,1:nb) ./ p_norm(1:nb);
  p_xy(2,1:nb) = r * p_xy(2,1:nb) ./ p_norm(1:nb);

  return
end


text:
%% A script to generate copy'able functions for quad-curl problem
% used in "Error analysis of a decoupled finite element method for quad-curl problems"
% https://arxiv.org/abs/2102.03396
%
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

clear; clc;
syms x y z
% U = [sin(y), -sin(z), sin(x)];
% U = [0, 0, (sin(x))^2*(sin(y))^2*sin(z)];
U = [0, 0, (sin(x)).^3.*(sin(y)).^3.*(sin(z)).^2];
% U = [0, 0, x^2*(1-x)^2*y^2*(1-y)^2*z^2*(1-z)^2];

X = [x y z];

%%
curlu = curl(U,X);
curlcurlu = curl(curlu,X);
tricurlu = curl(curlcurlu,X);
quadcurlu = curl(tricurlu,X);
pentacurlu = curl(quadcurlu,X);
hexacurlu = curl(pentacurlu,X);

vectorsStr = {'curlu', 'curlcurlu', 'tricurlu', ...
              'quadcurlu', 'pentacurlu', 'hexacurlu'};

%%
for idx = 1:length(vectorsStr)
    vector = eval(vectorsStr{idx});
    fprintf("function s = %s(p)\n",vectorsStr{idx})
    fprintf("x = p(:,1); y = p(:,2); z = p(:,3);\n")
    for j = 1:3
        vStr = char(vector(j));
        vStr = strrep(vStr, '*', '.*');
        vStr = strrep(vStr, '^', '.^');
        fprintf("s(:,%d) = %s;\n", j, vStr)
    end
    fprintf('end\n\n')
end

repo: ifem


text:
function out = innerProduct(f, g)
%INNERPRODUCT   Inner product of two ONEFUN objects.
%   INNERPRODUCT(F, G) returns the L2 inner product (on [-1,1]) of the two
%   ONEFUN objects F and G (conjugate linear in F).
%
% See also SUM.

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

% Deal with empty case:
if ( isempty(f) || isempty(g) )
    out = [];
    return
end

if ( ~isa(f, 'onefun') || ~isa(g, 'onefun') )
    error('CHEBFUN:SINGFUN:innerProduct:input', ...
        'innerProduct() only operates on two ONEFUN objects.');
end

m = size(f, 2);
n = size(g, 2);
cf = conj(f);

% Loop over columns of f and g:
out = zeros(m, n);
for j = 1:m
    fj = extractColumns(cf, j);    % jth column of f.
    for k = 1:n
        gk = extractColumns(g, k); % kth column of g.
        % Call SUM:
        out(j,k) = sum(fj.*gk);
    end
end

end


text:
function matrix = pseudoRotation(skel, channels, ind2Compute)

% PSEUDOGRADIENT ths position function of a node is composed with 2
% terms. The first one is the position of the parent node, the second is an
% additional term to get the real position. This function computes the
% derivative of the second term for the node 'ind2Compute' in respect with
% the channels 'rotInd2' given a set of channels.
%
%	Description:
%
%	[gradient matrix] = PSEUDOGRADIENT(skel, channels,gradientDirection, ind,rotInd2, ind2Compute)
%
%	 Returns:
%	  GRADIENT - the coponents of the gradient (x,y,z)
%     MATRIX - the rotation matrix
%	 Arguments:
%	  SKEL - a skeleton structure
%	  CHANNELS - the channels where the gradient has to be computes
%     GRADIENTDIRECTION - the direction of the gradient
%     IND - a recursive paremeter... start with 0 in general
%     ROTIND2 - the gradient has to be computed in respect with this
%     channel
%     IND2COMPUTE - the node ID where the gradient has to be computed

skelTmp = skel.tree(ind2Compute);
parent = skelTmp.parent;
skelParent = skel.tree(parent);
matrix = eye(3);

while (parent~=0)

    skelParent = skel.tree(parent);

rotVal = zeros(1, 3);
    for j = 1:length(skelParent.rotInd)
        rind = skelParent.rotInd(j);
        if rind
            rotVal(j) = channels(rind);
        else
            rotVal(j) = 0;
        end
    end

    tdof = rotationMatrix(deg2rad(rotVal(1)), ...
        deg2rad(rotVal(2)), ...
        deg2rad(rotVal(3)), ...
        skelParent.order);
    torient = rotationMatrix(deg2rad(skelParent.axis(1)), ...
        deg2rad(skelParent.axis(2)), ...
        deg2rad(skelParent.axis(3)), ...
        skelParent.axisOrder);
    torientInv = rotationMatrix(deg2rad(-skelParent.axis(1)), ...
        deg2rad(-skelParent.axis(2)), ...
        deg2rad(-skelParent.axis(3)), ...
        skelParent.axisOrder(end:-1:1));

    matrix = matrix*torientInv*tdof*torient;

    parent = skelParent.parent;

end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function datapt = lvd_KinematicStateTasks(stateLogEntry, subTask, inFrame)
%lvd_KinematicStateTasks Summary of this function goes here
%   Detailed explanation goes here

    cartElem = stateLogEntry.getCartesianElementSetRepresentation();
    cartElem = cartElem.convertToFrame(inFrame);

    datapt = -1;
    switch subTask
        % Cartesian Elements
        case 'rVectX'
            rVect = cartElem.rVect;
            datapt = rVect(1);
        case 'rVectY'
            rVect = cartElem.rVect;
            datapt = rVect(2);
        case 'rVectZ'
            rVect = cartElem.rVect;
            datapt = rVect(3);
        case 'vVectX'
            vVect = cartElem.vVect;
            datapt = vVect(1);
        case 'vVectY'
            vVect = cartElem.vVect;
            datapt = vVect(2);
        case 'vVectZ'
            vVect = cartElem.vVect;
            datapt = vVect(3);

        % Keplerian Elements
        case 'sma'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            datapt = kepElemSet.sma;
        case 'ecc'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            datapt = kepElemSet.ecc;
        case 'inc'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            datapt = rad2deg(kepElemSet.inc);
        case 'raan'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            datapt = rad2deg(kepElemSet.raan);
        case 'arg'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            datapt = rad2deg(kepElemSet.arg);
        case 'tru'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            datapt = rad2deg(kepElemSet.tru);
        case 'mean'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            datapt = rad2deg(kepElemSet.getMeanAnomaly());

        % Geographic Elements
        case 'lat'
            geoElemSet = cartElem.convertToGeographicElementSet();
            datapt = rad2deg(geoElemSet.lat);
        case 'long'
            geoElemSet = cartElem.convertToGeographicElementSet();
            datapt = rad2deg(geoElemSet.long);
        case 'alt'
            geoElemSet = cartElem.convertToGeographicElementSet();
            datapt = geoElemSet.alt;
        case 'velAz'
            geoElemSet = cartElem.convertToGeographicElementSet();
            datapt = rad2deg(geoElemSet.velAz);
        case 'velEl'
            geoElemSet = cartElem.convertToGeographicElementSet();
            datapt = rad2deg(geoElemSet.velEl);
        case 'velMag'
            geoElemSet = cartElem.convertToGeographicElementSet();
            datapt = geoElemSet.velMag;

        % Universal Elements - don't need to do Kep repeated elements
        case 'c3'
            univElemSet = cartElem.convertToUniversalElementSet();
            datapt = univElemSet.c3;
        case 'tau'
            univElemSet = cartElem.convertToUniversalElementSet();
            datapt = univElemSet.tau;

        % Misc.
        case 'period'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            datapt = kepElemSet.getPeriod();
        case 'rPe'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            datapt = kepElemSet.getRadiusPeriapsis();
        case 'rApo'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            datapt = kepElemSet.getRadiusApoapsis();
        case 'altPeri'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            datapt = kepElemSet.getAltitudePeriapsis();
        case 'altApo'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            datapt = kepElemSet.getAltitudeApoapsis();
        case 'H1'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            [datapt, ~, ~, ~] = kepElemSet.getEquinoctialElements();
        case 'K1'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            [~, datapt, ~, ~] = kepElemSet.getEquinoctialElements();
        case 'H2'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            [~, ~, datapt, ~] = kepElemSet.getEquinoctialElements();
        case 'K2'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            [~, ~, ~, datapt] = kepElemSet.getEquinoctialElements();
        case 'FPA'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            datapt = rad2deg(kepElemSet.getFlightPathAngle());
        case 'hyperVelUnitVectX'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            [datapt, ~, ~, ~, ~, ~] = kepElemSet.getOutboundHyperbolicVelocityElements();
        case 'hyperVelUnitVectY'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            [~, datapt, ~, ~, ~, ~] = kepElemSet.getOutboundHyperbolicVelocityElements();
        case 'hyperVelUnitVectZ'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            [~, ~, datapt, ~, ~, ~] = kepElemSet.getOutboundHyperbolicVelocityElements();
        case 'hyperVelUnitVectRA'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            [~, ~, ~, datapt, ~, ~] = kepElemSet.getOutboundHyperbolicVelocityElements();
        case 'hyperVelUnitVectDec'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            [~, ~, ~, ~, datapt, ~] = kepElemSet.getOutboundHyperbolicVelocityElements();
        case 'hyperVelMag'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            [~, ~, ~, ~, ~, datapt] = kepElemSet.getOutboundHyperbolicVelocityElements();

        case 'radius'
            datapt = cartElem.getRadiusMagnitude();
        case 'velocity'
            datapt = cartElem.getVelocityMagnitude();
        case 'horzVel'
            [datapt, ~] = cartElem.getHorzVertVelocities();
        case 'vertVel'
            [~, datapt] = cartElem.getHorzVertVelocities();

        case 'longDriftRate'
            kepElemSet = cartElem.convertToKeplerianElementSet();
            datapt = kepElemSet.getLongDriftRate() * (3600*180/pi); %convert rad/sec to deg/hr

        case 'centralBodyId'
            datapt = stateLogEntry.centralBody.id;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [gly_crop, gly_inrange] = crop_candidates(img_frame, curr_samples, template_size)
%create gly_crop, gly_inrange

nsamples = size(curr_samples,1);
c = prod([template_size 3]);
gly_inrange = zeros(nsamples,1);
gly_crop = zeros(c,nsamples);

for n = 1:nsamples
    curr_afnv = curr_samples(n, :);

    %    [img_cut, gly_inrange(n)] = IMGaffine_r(img_frame, curr_afnv, template_size);
    [img_cut, gly_inrange(n)] = IMGaffine_c(img_frame, curr_afnv, template_size);

    gly_crop(:,n) = reshape(img_cut, c , 1);
end

author: yuxng


text:
clear all
close all
path(path,'..\..\FUZZCLUST')
colors={'r.' 'gx' 'b+' 'ys' 'md' 'cv' 'k.' 'r*' 'g*' 'b*' 'y*' 'm*' 'c*' 'k*' };

%the data
data.X=nDexample(5,250,2,1);
%normalization
data=clust_normalize(data,'range');

%parameters
param.c=3;
param.m=2;
param.val=1;
param.vis=1;
%Kmeans clustering
result=kmedoid(data,param);

%validation
result = validity(result,data,param);
%
result.validity


text:
function [dbBool, noDocMets] = checkFormulaValidty(model)
% Assesses whether metabolites in model are likely to be documented in databases.
%
% USAGE:
%
%    [dbBool, noDocMets] = checkFormulaValidty(model)
%
% INPUT:
%    model:        Model structure array
%
% OUTPUTS:
%    dbBool:       A boolean vector where the number of rows is equal to the number of
%                  metabolites in model. Contains a logical 1 in rows for metabolites that
%                  are likely to be documented in databases and a logical 0 elsewhere.
%    noDocMets:    A cell array containing a list of all metabolites in model that
%                  are not likely to be documented in any database.
%
% .. Author: - Ronan Fleming

[nMets, nRxns] = size(model.S);

dbBool = ones(nMets,1);

for n = 1:nMets;

    if ~isempty(strfind(model.mets{n},'_hs'))
        dbBool(n)=0;
    end
    if ~isempty(strfind(model.mets{n},'hs_'))
        dbBool(n)=0;
    end
    if ~isempty(strfind(model.mets{n},'retn'))
        dbBool(n)=0;
    end
    if ~isempty(strfind(model.mets{n},'vitd'))
        dbBool(n)=0;
    end
    if ~isempty(strfind(model.mets{n},'rtotal'))
        dbBool(n)=0;
    end
    if ~isempty(strfind(model.mets{n},'deg'))
        dbBool(n)=0;
    end
    if ~isempty(strfind(model.metFormulas{n},'R'))
        dbBool(n)=0;
    end
    if ~isempty(strfind(model.metFormulas{n},'FULLR'))
        dbBool(n)=0;
    end
    if ~isempty(strfind(model.metFormulas{n},'X'))
        dbBool(n)=0;
    end

end

noDocMets = model.mets(~dbBool, 1);


text:
function plota(a,peval,savethis)
% plota(a,peval,savethis)

% figure
% plot(a','.:')
% xlabel('time');
% ylabel('a');
% grid on
% setfontsizefigure(12)
% if savethis
%     SaveImageFULL([peval.res_dir '/a'])
% end

figure
bar(mean(a,2))
hold on
errorbar(1:size(a,1), mean(a,2),std(a,[],2),'+r')
xlabel('component')
ylabel('mean(intensity)')
grid on
setfontsizefigure(12)
if savethis
    SaveImageFULL([peval.res_dir '/a_mean'])
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function dtiExportFibers(handles,varargin)
%
%  dtiExportFibers(handles,options);
%
% Exports non-mrDiffusion fiber/path file formats. Currently supports:
%
% 1. MetroTrac (*.pdb) format
%
% Options may include:
% * 'current' | 'selected' | 'all': save only the current fiber group, allow
% user to select some groups to save, or save all the fibers (default is
% current)
%
%
% HISTORY:
% 2007.06.15 AJS: wrote it.
%

fn = handles.defaultPath;
if(exist(fullfile(fn, 'fibers'),'dir')); fn = fullfile(fn, 'fibers', filesep); end
if isempty(handles.fiberGroups), warndlg('No fiber groups'); return; end

% Parse options
if(~isempty(varargin))
    if(~isempty(strmatch('all', lower(varargin)))); saveMode = 3;
    elseif(~isempty(strmatch('sel', lower(varargin)))); saveMode = 2;
    else saveMode = 1;
    end
    normFlag = ~isempty(strmatch('norm', lower(varargin)));
end

if(normFlag)
    if(length(handles.t1NormParams)>1)
        [id,OK] = selector(1:length(handles.t1NormParams), ...
        {handles.t1NormParams.name}, 'Select the coord space...');
        if(~OK); disp('user canceled.'); return; end
    else
        id = 1;
    end
    if(~isfield(handles.t1NormParams(id), 'name'))
        % in the old days, we didn't name the coord space.
        warning('no name associated with this coord space; using MNI.'); %#ok<WNTAG>
        handles.t1NormParams(id).name = 'MNI';
    end
    %coordSpace = handles.t1NormParams(id).name;
    if(isfield(handles.t1NormParams(id),'sn'))
         xform = handles.t1NormParams(id);
        if(~isfield(xform,'deformX')||isempty(xform.deformX))
        	disp('Computing inverse deformation...');
            [xform.deformX, xform.deformY, xform.deformZ] = mrAnatInvertSn(xform.sn);
            handles.t1NormParams(id).deformX = xform.deformX;
            handles.t1NormParams(id).deformY = xform.deformY;
            handles.t1NormParams(id).deformZ = xform.deformZ;
        end
        % xform that goes from acpc space to the deformation field space.
        xform.inMat = inv(xform.sn.VF.mat);
    else
        disp(['There is no xform specified- saving untransformed coords as "' coordSpace '".']);
        xform = [];
    end
else
    %coordSpace = 'acpc';
    xform = [];
end

% Save 'em.
if(saveMode==1)
    fg = handles.fiberGroups(handles.curFiberGroup);
    fn = fullfile(fn, [strrep(fg.name, ' ', '_') '.pdb']);
    [f, p] = uiputfile({'*.pdb';'*.*'}, ['Export Fibers "' fg.name '"...'], fn);
    if(isnumeric(f)), disp('Export current fibers canceled.'); return; end
    % For pdb format lets save out in ACPC space
    %mtrExportFibers(fg, fullfile(p, f), handles.xformToAcpc, xform);
    % Assume fiber groups are in AcPc space
    mtrExportFibers(fg, fullfile(p, f));
    disp('Current fibers exported.');
elseif(saveMode==2 || saveMode==3)
    p = uigetdir(fn, 'Fiber Group directory');
    if isnumeric(p), disp('Export FGs ... canceled.'), return; end
    if(saveMode==2)
        sList = dtiSelectFGs(handles,'Export FGs');
        if isempty(sList), disp('Export Fiber Groups ... canceled.'); return; end
    else
        sList = 1:length(handles.fiberGroups);
    end
    FGs = handles.fiberGroups(sList);
    for ii=1:length(FGs)
        fg = FGs(ii);
        fName = [strrep(fullfile(p,fg.name), ' ', '_') '.pdb'];
        mtrExportFibers(fg, fName, handles.xformToAcpc, xform);
    end
    disp(sprintf('Saved %.0f Fiber Groups.\n', length(FGs)))
end

return;

author: vistalab

repo: vistasoft