text:
function fem = update1320nbr(fem_struct,jj)
% update1320nbr takes a node with connectivity between 13 and 20
%  and adds six new nodes and uses delanunay trianglulation to
%  create the new elements. The region is then sprung. An nei is
%  not required but if one is not present, it will be generated.
%
%  Variables
%  fem_struct -- finite element structure from opnml.
%  jj -- the node number that has 12 elements connected to it.
%  fem -- the updated finite element structure.
%        Note:  Only fem.x,fem.y,fem.z,fem.e, and fem.nei
%               get updated.  fem.bnd does not need to be updated.
%
%  Usage -- fem = update1320nbr(fem_struct,jj)
%
%  FileName: update1320nbr.m
%  Written by: Ben Holladay (SEAP Student 2004)
%  Date: July 7,2004
%  Modified:  August 30, 2004 -- Chris Massey
%             June 14,2005 -- Ben Holladay
%                   Fix: Ensure that delaunay doesn't create elements
%                   outside the patch.
% Last Modifed:
%         Sept. 19, 2006 -- Chris Massey, NRL Code 7322, S.S.C. MS 39529
%                      % Add test to make sure the correct nodal neighbor
%                        connectivity existed for the requested update
%         July 14, 2008 -- Chris Massey, NRL Code 7322
%               moved test for 4 neighbors to after the test for nei.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Checks if the fem_struct has an nei. If not one is generated.
try
    nei = fem_struct.nei;
catch
    fem_struct.nei = ele2nei(fem_struct.e,fem_struct.x,fem_struct.y);
    nei = fem_struct.nei;
end

tempnc = sum(~ismember((fem_struct.nei(jj,:)),0));
if (tempnc < 13 | tempnc > 20)
    error(['There are ',num2str(tempnc),' nodal neighbors for node number ',...
           num2str(jj),'. This routine only works for 13-20 nodal neighbors.']);
end

%Sets the intial fem_struct variables.
x = fem_struct.x;
y = fem_struct.y;
z = fem_struct.z;
enodes = fem_struct.e;

%Adds new nodes,determines connectivity, and labels the bad node,
%neighboring nodesm and neighboring elements.
badnode = jj;
[nbrelem,~] = find(enodes == badnode);
[nnodes,trnc] = size(nei);
temp1 = sum(ismember(nei(jj,:),0));
nc = trnc - temp1;
nbrnode = nei(badnode,1:nc);
newnode = nnodes + [1:6];
vod = ([1:nc,1:nc]);
inrad = max(sqrt((x(nbrnode)-x(badnode)).^2 + (y(nbrnode)-y(badnode)).^2));

%Computes the center angle of each element.
b2 = (x(nbrnode(:))-x(badnode)).^2 + (y(nbrnode(:))-y(badnode)).^2;
a2 = (x(nbrnode(:))-x(nbrnode([2:nc,1]))).^2 + (y(nbrnode(:))-y(nbrnode([2:nc,1]))).^2;
A = (180/pi)*acos((b2+b2([2:nc,1])-a2)./(2*sqrt(b2).*sqrt(b2([2:nc,1]))));
ang = A;

%Sets the patterns for how to add the nodes.
addnode(1,:) = ([2,2,2,2,2,3]);
addnode(2,:) = ([2,2,3,2,2,3]);
addnode(3,:) = ([2,3,2,3,2,3]);
addnode(4,:) = ([3,3,2,3,3,2]);
addnode(5,:) = ([2,3,3,3,3,3]);
addnode(6,:) = ([3,3,3,3,3,3]);
addnode(7,:) = ([3,3,3,3,3,4]);
addnode(8,:) = ([3,3,4,3,3,4]);

%Determines the appropriate pattern for adding nodes
addnodepat = addnode((nc - 12),:);

%Labels reference nodes to be use to determine the new nodes locations.
spb(1) = 1;
spb(2) = spb(1) + addnodepat(1);
spb(3) = spb(2) + addnodepat(2);
spb(4) = spb(3) + addnodepat(3);
spb(5) = spb(4) + addnodepat(4);
spb(6) = spb(5) + addnodepat(5);

%Determines the total angle that each node will cover.
testang(1) = sum([ang([spb(1):(spb(2)-1)])]);
testang(2) = sum([ang([spb(2):(spb(3)-1)])]);
testang(3) = sum([ang([spb(3):(spb(4)-1)])]);
testang(4) = sum([ang([spb(4):(spb(5)-1)])]);
testang(5) = sum([ang([spb(5):(spb(6)-1)])]);
testang(6) = sum([ang([spb(6):nc])]);

%Test to make sure no node is trying to cover more than half of the
%circle.
test = find(testang > 180);
while ~isempty(test)
    hightemp = find(testang > 180);
    hhtempplus = hightemp + 1;
    if hhtempplus == 7; break; end
    spb(vod(hhtempplus)) = spb(vod(hhtempplus)) - 1;
    testang(1) = sum([ang([spb(1):(spb(2)-1)])]);
    testang(2) = sum([ang([spb(2):(spb(3)-1)])]);
    testang(3) = sum([ang([spb(3):(spb(4)-1)])]);
    testang(4) = sum([ang([spb(4):(spb(5)-1)])]);
    testang(5) = sum([ang([spb(5):(spb(6)-1)])]);
    testang(6) = sum([ang([spb(6):nc])]);
    test = find(testang > 180);
end

%First guess of the new coordinates and bathymetry.
M(1,1) = sum(x([(nbrnode([spb(1):spb(2)])),badnode,badnode]));
M(1,2) = sum(y([(nbrnode([spb(1):spb(2)])),badnode,badnode]));
M(1,3) = sum(z([(nbrnode([spb(1):spb(2)])),badnode,badnode]));

M(2,1) = sum(x([(nbrnode([spb(2):spb(3)])),badnode,badnode]));
M(2,2) = sum(y([(nbrnode([spb(2):spb(3)])),badnode,badnode]));
M(2,3) = sum(z([(nbrnode([spb(2):spb(3)])),badnode,badnode]));

M(3,1) = sum(x([(nbrnode([spb(3):spb(4)])),badnode,badnode]));
M(3,2) = sum(y([(nbrnode([spb(3):spb(4)])),badnode,badnode]));
M(3,3) = sum(z([(nbrnode([spb(3):spb(4)])),badnode,badnode]));

M(4,1) = sum(x([(nbrnode([spb(4):spb(5)])),badnode,badnode]));
M(4,2) = sum(y([(nbrnode([spb(4):spb(5)])),badnode,badnode]));
M(4,3) = sum(z([(nbrnode([spb(4):spb(5)])),badnode,badnode]));

M(5,1) = sum(x([(nbrnode([spb(5):spb(6)])),badnode,badnode]));
M(5,2) = sum(y([(nbrnode([spb(5):spb(6)])),badnode,badnode]));
M(5,3) = sum(z([(nbrnode([spb(5):spb(6)])),badnode,badnode]));

M(6,1) = sum(x([(nbrnode([spb(6):nc,1])),badnode,badnode]));
M(6,2) = sum(y([(nbrnode([spb(6):nc,1])),badnode,badnode]));
M(6,3) = sum(z([(nbrnode([spb(6):nc,1])),badnode,badnode]));

M(7,1) = x(badnode);
M(7,2) = y(badnode);
M(7,3) = z(badnode);

M(1,:) = M(1,:) ./ (addnodepat(1) + 3);
M(2,:) = M(2,:) ./ (addnodepat(2) + 3);
M(3,:) = M(3,:) ./ (addnodepat(3) + 3);
M(4,:) = M(4,:) ./ (addnodepat(4) + 3);
M(5,:) = M(5,:) ./ (addnodepat(5) + 3);
M(6,:) = M(6,:) ./ (addnodepat(6) + 3);

x([([newnode(1:6)]),badnode]) = ([M(:,1)]);
y([([newnode(1:6)]),badnode]) = ([M(:,2)]);
z([([newnode(1:6)]),badnode]) = ([M(:,3)]);

%Uses delaunay trianglulation to created a new and updated element list.
tempx = x([badnode,newnode(:)',nbrnode(:)']);
tempy = y([badnode,newnode(:)',nbrnode(:)']);
tri = delaunay(tempx,tempy);
tri2 = tri;

%Convertst the local element list from tri back to its global numberiing.
for i = 1:nc+7
    temp = find(i == tri);
    tmp = ([badnode,newnode(:)',nbrnode(:)']);
    tri2(temp) = tmp(i);
end

%Checks if the elements are in CCW order.
xnodes = x(tri2);
ynodes = y(tri2);
temparea = 0.5*(xnodes(:,1).*(ynodes(:,2)-ynodes(:,3))+xnodes(:,2).*...
    (ynodes(:,3)-ynodes(:,1))+xnodes(:,3).*(ynodes(:,1)-ynodes(:,2)));
Itri = find(temparea < 0);
tri3 = tri2;
tri3(Itri,2) = tri2(Itri,3);
tri3(Itri,3) = tri2(Itri,2);

%Removes elements that would be created outside the patch.
%!!MODIFICATION!!
oldbnd = detbndy(enodes(nbrelem,:));
newbnd = detbndy(tri3);
oldbnd = sort(oldbnd,2);
newbnd = sort(newbnd,2);
temp11 = find(ismember(newbnd,oldbnd,'rows') ~= 1);
while isempty(temp11) == 0
    newbnd = detbndy(tri3);
    newbnd = sort(newbnd,2);
    temp11 = find(ismember(newbnd,oldbnd,'rows') ~= 1);
    temp12 = newbnd(temp11,:);
    for it = 1:size(temp12,1)
        [tr,tc] = find(tri3 == temp12(it,1) | tri3 == temp12(it,2));
        [ta,tb,tc] = unique(tr);
        temp13 = setdiff(1:length(tr),tb);
        temp14 = tr(temp13);
        tri3(temp14,:) = [];
    end
end

%Finds and throws out bad elements and adds new elements to the global list.
tempa = setdiff(1:1:length(enodes),nbrelem);
enodes = enodes(tempa,1:3);
nelems = size(enodes,1);
tri3r = size(tri3,1);
enodes([(nelems + 1):(nelems + tri3r)],:) = tri3;

%Springs the region.
stoptol = 10e-8 * inrad;
tol = stoptol + 10;
i = 1;
imax = 20;
while i <= imax && tol > stoptol
    M2(7,1) = sum([M(2,1),M(1,1),M(3,1),M(5,1),M(4,1),M(6,1)]);
    M2(7,2) = sum([M(2,2),M(1,2),M(3,2),M(5,2),M(4,2),M(6,2)]);
    M2(7,3) = sum([M(2,3),M(1,3),M(3,3),M(5,3),M(4,3),M(6,3)]);

    [r1,c1] = find(enodes == newnode(1));
    nbr1 = unique(enodes(r1,:));
    temp1 = find(nbr1 == newnode(1));
    nbr1(temp1) = [];
	M2(1,1) = sum(x([nbr1]));
    M2(1,2) = sum(y([nbr1]));
    M2(1,3) = sum(z([nbr1]));

    [r2,c] = find(enodes == newnode(2));
    nbr2 = unique(enodes(r2,:));
    temp2 = find(nbr2 == newnode(2));
    nbr2(temp2) = [];
	M2(2,1) = sum(x([nbr2]));
    M2(2,2) = sum(y([nbr2]));
    M2(2,3) = sum(z([nbr2]));

    [r3,c3] = find(enodes == newnode(3));
    nbr3 = unique(enodes(r3,:));
    temp3 = find(nbr3 == newnode(3));
    nbr3(temp3) = [];
	M2(3,1) = sum(x([nbr3]));
    M2(3,2) = sum(y([nbr3]));
    M2(3,3) = sum(z([nbr3]));

    [r4,c4] = find(enodes == newnode(4));
    nbr4 = unique(enodes(r4,:));
    temp4 = find(nbr4 == newnode(4));
    nbr4(temp4) = [];
	M2(4,1) = sum(x([nbr4]));
    M2(4,2) = sum(y([nbr4]));
    M2(4,3) = sum(z([nbr4]));

    [r5,c5] = find(enodes == newnode(5));
    nbr5 = unique(enodes(r5,:));
    temp5 = find(nbr5 == newnode(5));
    nbr5(temp5) = [];
	M2(5,1) = sum(x([nbr5]));
    M2(5,2) = sum(y([nbr5]));
    M2(5,3) = sum(z([nbr5]));

	[r6,c6] = find(enodes == newnode(6));
    nbr6 = unique(enodes(r6,:));
    temp6 = find(nbr6 == newnode(6));
    nbr6(temp6) = [];
	M2(6,1) = sum(x([nbr6]));
    M2(6,2) = sum(y([nbr6]));
    M2(6,3) = sum(z([nbr6]));

    M2(1,:) = M2(1,:) ./ (size((nbr1),1));
    M2(2,:) = M2(2,:) ./ (size((nbr2),1));
    M2(3,:) = M2(3,:) ./ (size((nbr3),1));
    M2(4,:) = M2(4,:) ./ (size((nbr4),1));
    M2(5,:) = M2(5,:) ./ (size((nbr5),1));
    M2(6,:) = M2(6,:) ./ (size((nbr6),1));
    M2(7,:) = M2(7,:) ./ 6;

	tol = max(sqrt((M2(:,1)-M(:,1)).^2 + (M2(:,2)-M(:,2)).^2));

    M = M2;
    x([([newnode(1:6)]),badnode]) = ([M(:,1)]);
    y([([newnode(1:6)]),badnode]) = ([M(:,2)]);
    z([([newnode(1:6)]),badnode]) = ([M(:,3)]);
    i = i + 1;
end

%Updates the nei for the region.
tempnode = ([newnode,badnode,nbrnode]);
tempcount = size(tempnode,2);
for i = 1:tempcount
    [tr,tc] = find(tempnode(i) == enodes);
    tempnbr = unique(enodes(tr,:));
    tmp1 = find(tempnbr == tempnode(i));
    tempnbr(tmp1) = [];
    tempnbrx = x(tempnbr) - x(tempnode(i));
    tempnbry = y(tempnbr) - y(tempnode(i));
    angle = (atan2(tempnbry,tempnbrx))*(180/pi);
    tempang = find(angle < 0);
    angle(tempang) = angle(tempang) + 360;
    [tmp2,tempi] = sort(angle);
    tempnbr1 = tempnbr(tempi);
    tmpr = size(tempnbr,1);
    nei(tempnode(i),:) = ([tempnbr1(1:tmpr)',zeros(1,trnc-tmpr)]);
end

%Createst the output structure.
fem = fem_struct;
fem.x = x;
fem.y = y;
fem.z = z;
fem.e = enodes;
fem.nei = nei;
return

author: CHLNDDEV


text:
function [x,y,K]=postprocessSDP(newx,newy,prepinfo,newK)
% [x,y,K]=postprocessSDP(newx,newy,prepinfo,newK)
% postprocessSDP: Postprocesses an SDP solution using the info from
% preprocessing (preprocessSDP)
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi, preprocessSDP

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

if issparse(newx)
    x=sparse([]);
else
    x=[];
end
y=newy;
K=newK;

for j=length(prepinfo):-1:1
    op=prepinfo{j};
    switch op(1)
        case 0
            %Do nothing
            x=[newx(end-op(2)+1:end);x];
            newx=newx(1:end-op(2));
        case 1
            %Diagonal matrix
            x=[reshape(diag(newx(1:op(2))),op(2)^2,1);x];
            newx=newx(op(2)+1:end);
            K.l=K.l-op(2);
            K.s(end+1)=op(2);
            K.rsdpN=K.rsdpN+1;
    end
end

repo: REST


text:
function jac = p18_jac ( neqn, t, y )

%*****************************************************************************80
%
%% P18_JAC evaluates the jacobian for problem p18.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the jacobian.
%
%    Output, real JAC(NEQN,NEQN), the jacobian matrix.
%
  jac = zeros ( neqn, neqn );

  d = ( sqrt ( ( y(1).^2 + y(2).^2 ) ) ).^5;

  jac(1,3) = 1.0;
  jac(2,4) = 1.0;
  jac(3,1) = ( 2.0 * y(1).^2 - y(2).^2 ) / d;
  jac(3,2) = 3.0 * y(1) * y(2) / d;
  jac(4,1) = 3.0 * y(1) * y(2) / d;
  jac(4,2) = ( - y(1).^2 + 2.0 * y(2).^2 ) / d;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function linplus_test42 ( )

%*****************************************************************************80
%
%% TEST42 tests R8GE_NP_TRF, R8GE_NP_TRM, R8GE_NP_TRS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  m = 10;
  n = m;
  nrhs = 1;

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST42\n' );
  fprintf ( 1, '  For a matrix in general storage,\n' );
  fprintf ( 1, '  R8GE_NP_TRF factors without pivoting,\n' );
  fprintf ( 1, '  R8GE_NP_TRS solves factored systems.\n' );
  fprintf ( 1, '  R8GE_NP_TRM computes A*X for factored A.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M =    %d\n', m );
  fprintf ( 1, '  Matrix columns N = %d\n', n );
%
%  Set the matrix.
%
  [ a, seed ] = r8ge_random ( m, n, seed );
%
%  Set the desired solution.
%
  x(1:n) = 1.0E+00;
%
%  Compute the corresponding right hand side.
%
  b = r8ge_mxv ( m, n, a, x );
%
%  Factor the matrix.
%
  [ a_lu, info ] = r8ge_np_trf ( m, n, a );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST42 - Fatal error!\n' );
    fprintf ( 1, '  R8GE_NP_TRF declares the matrix is singular!\n' );
    fprintf ( 1, '  The value of INFO is %d\n', info );
    return
  end
%
%  Solve the linear system.
%
  b_mat = r8vec_to_r8ge ( n, nrhs, b );
  [ x_mat, info ] = r8ge_np_trs ( n, nrhs, 'N', a_lu, b_mat );
  x = r8ge_to_r8vec ( n, nrhs, x_mat );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST42 - Fatal error!\n' );
    fprintf ( 1, '  R8GE_TRS returned an error condition!\n' );
    fprintf ( 1, '  The value of INFO is %d\n', info );
    return
  end

  r8vec_print ( n, x, '  Solution:' );
%
%  Set the desired solution.
%
  x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
  job = 0;
  b = r8ge_np_trm ( m, n, a_lu, x, job );
%
%  Solve the system
%
  b_mat = r8vec_to_r8ge ( n, nrhs, b );
  [ x_mat, info ] = r8ge_np_trs ( n, nrhs, 'N', a_lu, b_mat );
  x = r8ge_to_r8vec ( n, nrhs, x_mat );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST42 - Fatal error!\n' );
    fprintf ( 1, '  R8GE_TRS returned an error condition!\n' );
    fprintf ( 1, '  The value of INFO is %d\n', info );
    return
  end

  r8vec_print ( n, x, '  Solution:' );
%
%  Set the desired solution.
%
  x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
  job = 1;
  b = r8ge_np_trm ( m, n, a_lu, x, job );
%
%  Solve the system.
%
  b_mat = r8vec_to_r8ge ( n, nrhs, b );
  [ x_mat, info ] = r8ge_np_trs ( n, nrhs, 'T', a_lu, b_mat );
  x = r8ge_to_r8vec ( n, nrhs, x_mat );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST42 - Fatal error!\n' );
    fprintf ( 1, '  R8GE_TRS returned an error condition!\n' );
    fprintf ( 1, '  The value of INFO is %d\n', info );
    return
  end

  r8vec_print ( n, x, '  Solution of transposed system:' );

  return
end


text:
function G = gsp_jtv_graph(G,T,fs,param)
% GSP_JTV_GRAPH Add time information to the graph structure
%   Usage:  G = gsp_jtv_graph(G);
%           G = gsp_jtv_graph(G,T);
%           G = gsp_jtv_graph(G,T,fs);
%           G = gsp_jtv_graph(G,T,fs,param);
%
%   Input parameters:
%         G          : Graph structure
%         T          : Time length
%         fs         : Sampling frequency (default 1)
%         param      : Structure of optional parameters
%
%   Output parameters:
%         G          : Time-Vertex Graph structure
%
%   This function adds the time domain information to the structure G.
%   The fields stored inside G.jtv are:
%   - T         : Time length
%   - fs        : Sampling frequency
%   - NFFT      : Number of frequency point (default [])
%   - Transform : Time Fourier basis: 'dft' or 'dct'
%   - Extension : Signal will be zero padded in time to take into account negative lag
%   - Lag       : Lag axis
%   - omega     : Frequency axis
%   - DiffT     : Time Gradient
%   - LT        : Time Laplacian
%
%   Additional parameters
%   ---------------------
%
%   * *param.transform*  : Fourier basis: 'dft' or 'dct'. (default 'dft')
%   * *param.approx*     : Stencil approx for the gradient: 'forward' or 'backward'. (default 'forward')
%   * *param.extension*  : Signal will be zero padded in time when needed. (default '0')
%   * *param.NFFT*       : Number of frequency point. (default [])
%

% Author :  Francesco Grassi
% Date   : September 2016

if nargin<4
    param = struct;
end

if ~isstruct(G)
    error('G is not a valid graph');
end

if nargin<2 || ~isnumeric(T) || T<1
    error('Time length T must be a numeric value greater than 0');
end

if nargin<3 || isempty(fs)
    fs=1;
end

if ~isnumeric(fs) || fs<0
    error('Sampling frequency fs must be a numeric value greater than 0');
end

if ~isfield(param,'NFFT'),      param.NFFT = []; end
if ~isfield(param,'transform'), param.transform = 'dft'; end
if ~isfield(param,'approx'),    param.approx = 'forward'; end
if ~isfield(param,'extension'), param.extension = 0; end

%% JTV PARAMETERS
G.jtv = struct;
G.jtv.T = T;
G.jtv.fs = fs;
G.jtv.NFFT = param.NFFT;
G.jtv.transform = param.transform;
G.jtv.extension = param.extension;
G.jtv.omega = gsp_jtv_fa(G);

if G.jtv.extension
    G.jtv.lag = 2*G.jtv.T-1;
else
    G.jtv.lag = T;
end

%% DIFFERENTIAL OPERATORS
z = ones(T,1);

% G.jtv.LT = spdiags([-z 2*z -z], [-1 0 1], T, T);

G.jtv.DiffT = spdiags([-z z], [0 1], T, T);

switch param.transform
    case 'dft'
        G.jtv.DiffT(T,1) = 1;
    case 'dct'
        G.jtv.DiffT(T,T-1) = 1;
    otherwise
        error('Unknown transform');
end

G.jtv.LT = G.jtv.DiffT'*G.jtv.DiffT;

if strcmpi(param.approx,'backward')
    G.jtv.DiffT = -G.jtv.DiffT.';
end

% toeplitz([-1 zeros(T-2,1) 1],[-1 1 zeros(T-2,1)]); %periodic forward
% toeplitz([1 -1 zeros(T-2,1)],[1 zeros(T-2,1) -1]); %periodic backward (= -transpose(periodic forward))

end

repo: gspbox


text:
% LIONSAMBA example script
% Continuous custom current profile: this script shows how to run a continuous custom current
% profile charge/discharge.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

% Clear the workspace
clear
close all

%% Parameters
% Define the integration times.
t0 = 0;
tf = 10^4;

% Define the initial state structure
initialState.Y  = [];
initialState.YP = [];

% Define the parameters structure.
param{1}               = Parameters_init;

param{1}.Np            = 10;
param{1}.Ns            = 10;
param{1}.Nn            = 10;

param{1}.hcell         = 1;
param{1}.Tref          = 298.15;

param{1}.AbsTol        = 1e-6;
param{1}.RelTol        = 1e-6;

param{1}.CutoffSOC     = 20;

param{1}.SolidPhaseDiffusion = 1;

if(param{1}.SolidPhaseDiffusion == 3)
    multp = param{1}.Nr_p;
    multn = param{1}.Nr_n;
else
    multp = 1;
    multn = 1;
end
I1C = 29.5;

C_rate = 1.5;
%% Discharge section
% Discharge the battery to the 20% of SOC
out = startSimulation(t0,tf,initialState,-25,param);

% Store the Jacobian matrix for future computations. This is possible
% because the different scenarios share the same model structure, and the
% only quantity which differs is the applied current density.
param{1}.JacobianFunction = out.JacobianFun;

% Update the initial states
initialState = out.initialState;

% Redefine the cutoff SOC in order to avoid simulation interruptions
param{1}.CutoffSOC     = 2;

% Rest the battery with no current applied
out2 = startSimulation(0,5000,initialState,0,param);

% Update the initial states
initialState = out2.initialState;

% Apply a custom current profile
param{1}.OperatingMode = 4;

% Define which external function is going to be used from the simulato to
% gather the value of the applied current density at each time step
param{1}.CurrentDensityFunction = @getInputCurrent;

% Run the simulation
out3 = startSimulation(0,5000,initialState,0,param);

% Update the initial states
initialState = out3.initialState;

% Go back to galvanostatic operating conditions
param{1}.OperatingMode= 1;

% Run the simulation
out4 = startSimulation(0,5000,initialState,0,param);
%% Plot the results

% Concatenate all the time results
time = [out.time{1};out2.time{1}+out.time{1}(end);out3.time{1}+out2.time{1}(end)+out.time{1}(end);out4.time{1}+out3.time{1}(end)+out2.time{1}(end)+out.time{1}(end)];

figure(1)
plot(time,[out.Voltage{1};out2.Voltage{1};out3.Voltage{1};out4.Voltage{1}],'LineWidth',6)
hold on
xlabel('Time [s]')
ylabel('Voltage [V]')
grid on
box on
title('Cell Voltage')

figure(2)
plot(time,[out.SOC{1};out2.SOC{1};out3.SOC{1};out4.SOC{1}],'LineWidth',6)
hold on
xlabel('Time [s]')
ylabel('SOC [%]')
grid on
box on
title('Cell SOC')

figure(3)
plot(time,[out.Temperature{1}(:,end);out2.Temperature{1}(:,end);out3.Temperature{1}(:,end);out4.Temperature{1}(:,end)],'LineWidth',6)
hold on
xlabel('Time [s]')
ylabel('Temperature[K]')
grid on
box on
title('Cell Temperature')

figure(4)
plot(time,[out.curr_density;out2.curr_density;out3.curr_density;out4.curr_density],'LineWidth',6)
hold on
xlabel('Time [s]')
ylabel('Applied Current [A/m^2]')
grid on
box on
title('Cell Input Current')

figure(5)
plot(time,[out.cs_average{1}(:,[1 param{1}.Np*multp]);out2.cs_average{1}(:,[1 param{1}.Np*multp]);out3.cs_average{1}(:,[1 param{1}.Np*multp]);out4.cs_average{1}(:,[1 param{1}.Np*multp])],'LineWidth',6)
hold on
line([time(1) time(end)],[param{1}.cs_maxp param{1}.cs_maxp],'LineWidth',6,'LineStyle','--')
xlabel('Time [s]')
ylabel('Positive electrode average concentration [mol/m^3]')
grid on
box on
legend('c^* x=0','c^* x=l_p','c^* max')

figure(6)
plot(time,[out.cs_average{1}(:,[(param{1}.Np+1)*multn (param{1}.Np+param{1}.Nn)*multn]);out2.cs_average{1}(:,[param{1}.Np+1 param{1}.Np+param{1}.Nn]*multn);out3.cs_average{1}(:,[param{1}.Np+1 param{1}.Np+param{1}.Nn]*multn);out4.cs_average{1}(:,[param{1}.Np+1 param{1}.Np+param{1}.Nn]*multn)],'LineWidth',6)
hold on
line([time(1) time(end)],[param{1}.cs_maxn param{1}.cs_maxn],'LineWidth',6,'LineStyle','--')
xlabel('Time [s]')
ylabel('Negative electrode average concentration [mol/m^3]')
grid on
box on
legend('c^* x=l_p+l_s','c^* x=l_p+l_s+l_n','c^* max')

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function [sampft, sampyt] = gp_rnd(gp, x, y, varargin)
%GP_RND  Random draws from the posterior Gaussian process
%
%  Description
%    [SAMPFT, SAMPYT] = GP_RND(GP, X, Y, XT, OPTIONS) takes a
%    Gaussian process structure, record structure (from gp_mc) or
%    array (from gp_ia) GP together with a matrix XT of input
%    vectors, matrix X of training inputs and vector Y of training
%    targets, and returns a random sample SAMPFT and SAMPYT from
%    the posterior distribution p(ft|x,y,xt) and the predictive
%    distribution p(yt|x,y,xt) at locations XT.
%
%    OPTIONS is optional parameter-value pair
%      nsamp     - determines the number of samples (default = 1).
%      predcf    - index vector telling which covariance functions are
%                  used for prediction. Default is all (1:gpcfn)
%      tstind    - a vector defining, which rows of X belong to which
%                  training block in *IC type sparse models. Default is [].
%                  See also GP_PRED.
%      z         - optional observed quantity in triplet (x_i,y_i,z_i)
%                  Some likelihoods may use this. For example, in case of
%                  Poisson likelihood we have z_i=E_i, that is, expected
%                  value
%                  for ith case.
%      fcorr     - Method used for latent marginal posterior corrections.
%                  Default is 'off'. Possible methods are 'fact' for EP
%                  and either 'fact', 'cm2' or 'lr' for Laplace. If method is
%                  'on', 'fact' is used for EP and 'cm2' for Laplace.
%                  See GP_PREDCM and Cseke & Heskes (2011) for more information.
%      splitnormal - determines if the samples are drawn from
%                  split-Normal approximation (Geweke, 1989) in the
%                  case of non-Gaussian likelihood.
%                  Possible values are 'on' (default) and 'off'.
%      n_scale   - the maximum number of the most significant principal
%                  component directories to scale in the split-Normal
%                  approximation (default 50).
%
%    If likelihood is non-Gaussian and gp.latent_method is Laplace the
%    samples are drawn from split-Normal approximation (see option
%    splitnormal), and if gp.latent_method is EP the samples are drawn
%    from the normal approximation.
%
%  Reference
%    Cseke & Heskes (2011). Approximate Marginals in Latent Gaussian
%    Models. Journal of Machine Learning Research 12 (2011), 417-454
%
%    Geweke, J. (1989).  Bayesian inference in econometric models using
%    Monte Carlo integration. Econometrica 57:1317-1339.
%
%  See also
%    GP_PRED, GP_PAK, GP_UNPAK

%  Internal comments
%    - sampling with FIC, PIC and CS+FIC forms full nxn matrix and
%       works only when sampling for the training inputs
%    - The code is not optimized
%
% Copyright (c) 2007-2010 Jarno Vanhatalo
% Copyright (c) 2008      Jouni Hartikainen
% Copyright (c) 2011      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.

ip=inputParser;
ip.FunctionName = 'GP_RND';
ip.addRequired('gp',@(x) isstruct(x) || iscell(x));
ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addOptional('xt', [], @(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))))
ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('zt', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('predcf', [], @(x) isempty(x) || ...
                 isvector(x) && isreal(x) && all(isfinite(x)&x>0))
ip.addParamValue('tstind', [], @(x) isempty(x) || iscell(x) ||...
                 (isvector(x) && isreal(x) && all(isfinite(x)&x>0)))
ip.addParamValue('nsamp', 1, @(x) isreal(x) && isscalar(x))
ip.addParamValue('fcorr', 'off', @(x) ismember(x, {'fact','cm2','off','on','lr'}))
ip.addParamValue('splitnormal', 'on', @(x) (islogical(x) && isscalar(x))|| ...
                 ismember(x,{'on' 'off' 'full'}))
ip.addParamValue('n_scale',50, @(x) isnumeric(x) && x>=0);
if numel(varargin)==0 || isnumeric(varargin{1})
  % inputParser should handle this, but it doesn't
  ip.parse(gp, x, y, varargin{:});
else
  ip.parse(gp, x, y, [], varargin{:});
end
xt=ip.Results.xt;
z=ip.Results.z;
zt=ip.Results.zt;
tn = size(x,1);
predcf=ip.Results.predcf;
tstind=ip.Results.tstind;
nsamp=ip.Results.nsamp;
fcorr=ip.Results.fcorr;
splitnormal = ip.Results.splitnormal;
n_scale = min(tn, floor(ip.Results.n_scale));
if isempty(xt)
  xt=x;
  if isempty(tstind)
    if iscell(gp)
      gptype=gp{1}.type;
    else
      gptype=gp.type;
    end
    switch gptype
      case {'FULL' 'VAR' 'DTC' 'SOR'}
        tstind = [];
      case {'FIC' 'CS+FIC'}
        tstind = 1:size(x,1);
      case 'PIC'
        if iscell(gp)
          tstind = gp{1}.tr_index;
        else
          tstind = gp.tr_index;
        end
    end
  end
end

if isfield(gp, 'latent_method')
  if iscell(gp)
      gplatmet=gp{1}.latent_method;
  else
      gplatmet=gp.latent_method;
  end
  if ~strcmp(gplatmet, 'Laplace') && strcmp(splitnormal,'on')
      % splitnormal is applicable only with Laplace
      splitnormal='off';
  end
end

sampyt=[];
if isstruct(gp) && numel(gp.jitterSigma2)==1
  % Single GP
  if isfield(gp, 'monotonic') && gp.monotonic
      [gp,x,y,z,xt,zt] = gp.fh.setUpDataForMonotonic(gp,x,y,z,xt,zt);
  end

  if (isfield(gp.lik.fh,'trcov') && ~isfield(gp,'lik_mono')) ...
      || isfield(gp, 'latentValues')
    % ====================================================
    %    Gaussian likelihood or MCMC with latent values
    % ====================================================

    if nargout > 1
      [Eft, Covft, ~, Eyt, Covyt] ...
        = gp_jpred(gp,x,y,xt,'z',z,'predcf',predcf,'tstind',tstind);
    else
      [Eft, Covft] = gp_jpred(gp,x,y,xt,'z',z, ...
                              'predcf',predcf,'tstind',tstind);
    end
    rr = randn(size(Eft,1),nsamp);
    sampft = bsxfun(@plus, Eft, chol(Covft,'lower')*rr);
    if nargout > 1
      sampyt = bsxfun(@plus, Eyt, chol(Covyt,'lower')*rr);
    end

  else
    % =============================
    %    non-Gaussian likelihood
    % =============================

    if nargout > 1
      warning('Sampling of yt is not implemented for non-Gaussian likelihoods');
    end
    if strcmp(splitnormal,'on') && isfield(gp,'meanf')
        warning('Split-Normal is not implemented for mean functions');
        splitnormal='off';
    end
    if strcmp(splitnormal,'on') && isfield(gp.lik, 'nondiagW')
        warning('Split-Normal is not implemented for likelihoods with nondiagonal Hessian')
        splitnormal='off';
    end
    if strcmp(splitnormal,'on') && ~isequal(fcorr, 'off')
        warning('Split-Normal is not used when latent marginal posterior corrections are used')
        splitnormal='off';
    end

    if strcmp(splitnormal,'on')
      % ------------------
      %    Splitnormal on
      % ------------------

      switch gp.type

        % ---------------------------
        case 'FULL'

          [e,~,~, param] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);
          if isnan(e)
            error('Laplace-algorithm returned NaN');
          end
          [Ef, L, W] = deal(param.f, param.L, param.La2);

          % Evaluate the covariance
          K = gp_trcov(gp,x);
          K_ss = gp_trcov(gp,xt,predcf);
          K_nf = gp_cov(gp,xt,x,predcf);
          if W >= 0
            % Likelihood is log concave
            if issparse(K) && issparse(L)
              if issparse(W)
                sqrtWK = sqrt(W)*K;
              else
                sqrtWK = sparse(1:tn,1:tn,sqrt(W),tn,tn)*K;
              end
              Covf = K - sqrtWK'*ldlsolve(L,sqrtWK);
            else
              V = linsolve(L,bsxfun(@times,sqrt(W),K),struct('LT',true));
              Covf = K - V'*V;
            end
          else
            % Likelihood is not log concave
            V = bsxfun(@times,L,W');
            Covf = K - K*(diag(W) - V'*V)*K;
          end

        % ---------------------------
        case 'FIC'

          if ~isempty(tstind) && length(tstind) ~= tn
            error('tstind (if provided) has to be of same length as x.')
          end

          % Ensure the direction of tstind
          tstind = tstind(:);

          [e,~,~, param] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);
          if isnan(e)
            error('Laplace-algorithm returned NaN');
          end
          Ef = param.f;

          u = gp.X_u;
          K_uu = gp_trcov(gp,u,predcf);
          K_uu = (K_uu+K_uu')./2;          % ensure the symmetry
          Luu = chol(K_uu, 'lower');

          B = Luu\gp_cov(gp,u,x,predcf);
          B2 = Luu\gp_cov(gp,u,xt,predcf);
          K = B'*B;
          K_ss = B2'*B2;
          K_nf = B2'*B;
          clear K_uu Luu B B2

          % If the prediction is made for training set, evaluate Lav
          % also for prediction points.
          if ~isempty(tstind)
            kss = gp_trvar(gp,xt,predcf);
            % Replace diagonals
            K(1:length(K)+1:numel(K)) = kss(tstind);
            K_ss(1:length(K_ss)+1:numel(K_ss)) = kss;
            % Replace common samples in K_nf
            K_nf(tstind+(0:size(xt,1):(tn-1)*size(xt,1))') = kss(tstind);
            clear kss
          else
            % Add lambda
            % Replace diagonals
            K(1:length(K)+1:numel(K)) = gp_trvar(gp,x,predcf);
            K_ss(1:length(K_ss)+1:numel(K_ss)) = gp_trvar(gp,xt,predcf);
          end

          Wsqrt = -gp.lik.fh.llg2(gp.lik, y, Ef, 'latent', z);
          if any(Wsqrt < 0)
            error('FIC not implemented for non-log-concave likelihoods')
          end
          Wsqrt = sqrt(Wsqrt);
          L = chol(eye(size(K)) + bsxfun(@times,bsxfun(@times,Wsqrt,K),Wsqrt'), 'lower');

          % Evaluate the covariance for the posterior of f
          V = linsolve(L,bsxfun(@times,Wsqrt,K),struct('LT',true));
          Covf = K - V'*V;

        % ---------------------------
        case {'PIC' 'PIC_BLOCK'}

          u = gp.X_u;
          K_fu = gp_cov(gp,x,u,predcf);         % f x u
          K_uu = gp_trcov(gp,u,predcf);         % u x u, noiseles covariance K_uu
          K_uu = (K_uu+K_uu')./2;               % ensure the symmetry of K_uu

          ind = gp.tr_index;
          Luu = chol(K_uu)';
          B=Luu\(K_fu');
          B2 = Luu\(gp_cov(gp,u,xt,predcf));
          clear Luu

          [e,~,~, p] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);
          if isnan(e)
            error('Laplace-algorithm returned NaN');
          end
          [Ef, La2] = deal(p.f, p.La2);

          % Evaluate the variance
          sqrtW = -gp.lik.fh.llg2(gp.lik, y, Ef, 'latent', z);
          if any(sqrtW < 0)
            error('PIC not implemented for non-log-concave likelihoods')
          end
          sqrtW = sqrt(sqrtW);
          % Components for (I + W^(1/2)*(Qff + La2)*W^(1/2))^(-1) = Lahat^(-1) - L2*L2'
          Lahat = cell(length(ind),1);
          for i=1:length(ind)
            Lahat{i} = eye(length(ind{i})) + bsxfun(@times,bsxfun(@times,sqrtW(ind{i}),La2{i}),sqrtW(ind{i})');
          end
          sKfu = bsxfun(@times, sqrtW, K_fu);
          iLasKfu = zeros(size(K_fu));
          for i=1:length(ind)
            iLasKfu(ind{i},:) = Lahat{i}\sKfu(ind{i},:);
          end
          A2 = K_uu + sKfu'*iLasKfu; A2=(A2+A2')./2;
          L2 = iLasKfu/chol(A2);

          % NOTE!
          % This is done with full matrices at the moment.
          % Needs to be rewritten.
          K_ss = B2'*B2;
          K_nf = B2'*B;
          K = B'*B;
          C = -L2*L2';
          for i=1:length(ind)
            % Originally
            % >> La = gp_trcov(gp, xt(tstind{i},:), predcf) - B2(:,tstind{i})'*B2(:,tstind{i});
            % >> K_ss(ind{i},ind{i}) =  K_ss(ind{i},ind{i}) + La;
            % changed into (works if xt is not the same as x)
            % >> La = gp_trcov(gp, xt(tstind{i},:), predcf) - B2(:,tstind{i})'*B2(:,tstind{i});
            % >> K_ss(tstind{i},tstind{i}) =  K_ss(tstind{i},tstind{i}) + La;
            % which is implemented in the line bellow
            K_ss(tstind{i},tstind{i}) =  gp_trcov(gp, xt(tstind{i},:), predcf);
            K_nf(tstind{i},ind{i}) = gp_cov(gp, xt(tstind{i},:), x(ind{i},:),predcf);
            K(ind{i},ind{i}) = K(ind{i},ind{i}) + La2{i};
            C(ind{i},ind{i}) =  C(ind{i},ind{i}) + inv(Lahat{i});
          end

          Covf = K - K * bsxfun(@times,bsxfun(@times,sqrtW,C),sqrtW') * K;

        % ---------------------------
        case 'CS+FIC'

          if ~isempty(tstind) && length(tstind) ~= tn
            error('tstind (if provided) has to be of same length as x.')
          end

          % Ensure the direction of tstind
          tstind = tstind(:);

          % Indexes to all non-compact support and compact support covariances.
          cscf = cellfun(@(x) isfield(x,'cs'), gp.cf);
          predcf1 = find(~cscf);
          predcf2 = find(cscf);
          if ~isempty(predcf)
            predcf1 = intersect(predcf1,predcf);
            predcf2 = intersect(predcf2,predcf);
          end

          % Laplace approximation
          [e,~,~, param] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);
          if isnan(e)
            error('Laplace-algorithm returned NaN');
          end
          Ef = param.f;

          u = gp.X_u;

          K_uu = gp_trcov(gp,u,predcf1);
          K_uu = (K_uu+K_uu')./2;
          Luu = chol(K_uu)';

          B=Luu\(gp_cov(gp,u,x,predcf1));
          B2=Luu\(gp_cov(gp,u,xt,predcf1));
          K = B'*B;
          K_ss = B2'*B2;
          K_nf = B2'*B;
          clear K_uu Luu B B2

          if ~isempty(tstind)
            % Add Lambda
            kss = gp_trvar(gp,xt,predcf1);
            % Replace diagonals
            K(1:length(K)+1:numel(K)) = kss(tstind);
            K_ss(1:length(K_ss)+1:numel(K_ss)) = kss;
            % Replace common samples in K_nf
            K_nf(tstind+(0:size(xt,1):(tn-1)*size(xt,1))') = kss(tstind);

            % Add CS covariance
            kss = gp_trcov(gp, xt, predcf2);
            K = K + kss(tstind,tstind);
            K_ss = K_ss + kss;
            clear kss

          else
            % Add lambda
            % Replace diagonals
            K(1:length(K)+1:numel(K)) = gp_trvar(gp,x,predcf1);
            K_ss(1:length(K_ss)+1:numel(K_ss)) = gp_trvar(gp,xt,predcf1);

            % Add CS covariance
            K = K + gp_trcov(gp, x, predcf2);
            K_ss = K_ss + gp_trcov(gp, xt, predcf2);

          end

          % Add CS covariance for K_nf
          K_nf = K_nf + gp_cov(gp, xt, x, predcf2);

          sqrtW = -gp.lik.fh.llg2(gp.lik, y, Ef, 'latent', z);
          if any(sqrtW < 0)
            error('CS+FIC not implemented for non-log-concave likelihoods')
          end
          sqrtW = sqrt(sqrtW);
          L = chol(eye(size(K)) + bsxfun(@times,bsxfun(@times,sqrtW,K),sqrtW'), 'lower');
          V = linsolve(L,bsxfun(@times,sqrtW,K),struct('LT',true));
          Covf = K - V'*V;

        otherwise
          error('Split-Normal not implemented for %s', gp.type)

      end

      % Scaling n_scale principal component directions
      % Split-normal approximation: Geweke (1989)

      % N.B. svd of sparse matrix is not generally sparse (unless Covf is
      % block diagonal or can be made such by reordering the dimensions).
      % Scaling using the Cholesky factorisation could possibly be done
      % sparsely but the directions would be different (not split-Normal).
      % Thus full matrices has to be used.
      [V, D, ~] = svd(full(Covf));
      T = real(V) * sqrt(real(D)); % Ensuring the real

      L = chol(K,'lower');

      % Optimise the skewness
      delta = [-6:0.5:-0.5, 0.5:0.5:6];
      D = bsxfun(@plus, kron(delta,T(:,1:n_scale)), Ef);
      ll = zeros(1,n_scale*length(delta)); % Here looping is faster than cellfun
      for i = 1:n_scale*length(delta)
        ll(i) = 2*gp.lik.fh.ll(gp.lik,y,D(:,i),z);
      end
      ll = ll - sum(D.*(L'\(L\D)));
      ll0 = 2*gp.lik.fh.ll(gp.lik, y, Ef, z) - Ef'*(L'\(L\Ef));
      ll(ll >= ll0) = NaN;
      f = bsxfun(@rdivide, abs(delta), reshape(sqrt(ll0-ll),n_scale,length(delta)));
      clear D V ll ll0
      q = max(f(:,delta>0),[],2);
      r = max(f(:,delta<0),[],2);
      q(isnan(q)) = 1;
      r(isnan(r)) = 1;
      % Safety limits [1/10, 10]
      q = min(max(q,0.1),10);
      r = min(max(r,0.1),10);
      if n_scale < tn
        q = [q;ones(tn-n_scale,1)];
        r = [r;ones(tn-n_scale,1)];
      end

      % Draw samples from split-Normal approximated posterior of f and
      % sample from the conditional distribution of ft for each sample of f
      u = rand(tn,nsamp);
      c = r./(r+q);
      Covft = K_ss - K_nf*(L'\(L\K_nf'));
      Covft_ind = find(diag(Covft)>gp.jitterSigma2);
      if length(Covft_ind) == size(xt,1)
        % The following evaluates samples of ft directly by
        % Eft + chol(Covft)*randn,
        % where Eft = K_nf*K\f and Covft = K_ss - K_nf*K\K_nf'
        sampft = K_nf*(L'\(L\ ...
          bsxfun(@plus,Ef,T*((bsxfun(@times,q,double(bsxfun(@ge,u,c))) ...
          +bsxfun(@times,-r,double(bsxfun(@lt,u,c)))).*abs(randn(tn,nsamp)))) ...
          )) + chol(Covft,'lower')*randn(size(xt,1),nsamp);
      else
        % Zero variance dimensions in p(ft|X,xt,f)
        % Evaluate first Eft for each f
        sampft = K_nf*(L'\(L\ ...
            bsxfun(@plus,Ef,T*((bsxfun(@times,q,double(bsxfun(@ge,u,c))) ...
            +bsxfun(@times,-r,double(bsxfun(@lt,u,c)))).*abs(randn(tn,nsamp)))) ));
        % Add variability only from non-zero variance dimensions
        sampft(Covft_ind,:) = sampft(Covft_ind,:) ...
          + chol(Covft(Covft_ind,Covft_ind),'lower')*randn(length(Covft_ind),nsamp);
      end

    else
      % -------------------
      %    split-Normal off
      % -------------------
      [Eft, Covft] = gp_jpred(gp,x,y,xt,'z',z, ...
        'predcf',predcf, 'tstind',tstind, 'fcorr','off');
      sampft = bsxfun(@plus, Eft, chol(Covft,'lower')*randn(size(xt,1),nsamp));

      if ~isequal(fcorr, 'off')
        % Do marginal corrections for samples
        [pc_predm, fvecm] = gp_predcm(gp, x, y, xt, 'z', z, 'ind', 1:size(xt,1), 'fcorr', fcorr);
        for i=1:size(xt,1)
          % Remove NaNs and zeros
          pc_pred=pc_predm(:,i);
          dii=isnan(pc_pred)|pc_pred==0;
          pc_pred(dii)=[];
          fvec=fvecm(:,i);
          fvec(dii)=[];
          % compute cdf
          cumsumpc = cumsum(pc_pred)/sum(pc_pred);
          % Remove non-unique values from grid vector & distribution
          [cumsumpc, inds] = unique(cumsumpc);
          fvec = fvec(inds);
          % use inverse cdf to make marginal transformation
          fsnc = normcdf(sampft(i,:), Eft(i,1), sqrt(diag(Covft(i,i))));
          sampft(i,:) = interp1(cumsumpc,fvec,fsnc);
        end
      end
    end

  end

elseif isstruct(gp) && numel(gp.jitterSigma2)>1
  % MCMC
  nmc=size(gp.jitterSigma2,1);
  % resample nsamp cases from nmc samples
  % deterministic resampling has low variance and small bias for
  % equal weights
  gi=resampdet(ones(nmc,1),nsamp,1);
  sampft = zeros(size(xt,1),nsamp);
  if nargout>=2
    sampyt = zeros(size(xt,1),nsamp);
  end
  for i1=1:nmc
    sampind = find(gi==i1);
    nsampi = length(sampind);
    if nsampi>0
      Gp = take_nth(gp,i1);
      if isfield(Gp,'latent_method') && isequal(Gp.latent_method,'MCMC')
        Gp = rmfield(Gp,'latent_method');
      end
      if isfield(gp, 'latentValues') && ~isempty(gp.latentValues)
        % Non-Gaussian likelihood. The latent variables should be used in
        % place of observations
        for ni=1:nsampi
          if ni==1
            % use stored latent values
            f = gp.latentValues(i1,:);
          else
            % need to resample latent values
            opt=scaled_mh();
            f=scaled_mh(f, opt, Gp, x, y, z);
          end
          if nargout<2
            sampft(:,sampind(ni)) = gp_rnd(Gp, x, f', xt, 'nsamp', 1, ...
                             'z', z, 'zt', zt, 'predcf', predcf, ...
                             'tstind', tstind);
          else
            [tsampft, tsampyt] = gp_rnd(Gp, x, f', xt, 'nsamp', ...
                                        1, 'z', z, 'zt', zt, ...
                                        'predcf', predcf, 'tstind', tstind);

            sampft(:,sampind(ni)) = tsampft;
            sampyt(:,sampind(ni)) = tsampyt;
          end

        end
      else
        % Gaussian likelihood
        if nargout<2
          sampft(:,sampind) = gp_rnd(Gp, x, y, xt, 'nsamp', nsampi, ...
                           'z', z, 'zt', zt, 'predcf', predcf, ...
                           'tstind', tstind);
        else
          [tsampft, tsampyt] = gp_rnd(Gp, x, y, xt, 'nsamp', ...
                                      nsampi, 'z', z, 'zt', zt, ...
                                      'predcf', predcf, 'tstind', tstind);
          sampft(:,sampind) = tsampft;
          sampyt(:,sampind) = tsampyt;
        end
      end
    end
  end
elseif iscell(gp)
  % gp_ia
  ngp=length(gp);
  if isfield(gp{1},'ia_weight')
    gw = zeros(ngp,1);
    for i1=1:ngp
      gw(i1)=gp{i1}.ia_weight;
    end
  else
    gw=ones(ngp,1);
  end
  % resample nsamp cases from nmc samples
  % stratified resampling has has higher variance than deterministic
  % resampling, but has a smaller bias, and thus it should be more
  % suitable for unequal weights
  gi=resampstr(gw,nsamp,1);
  sampft = zeros(size(xt,1),nsamp);
  if nargout>=2
    sampyt = zeros(size(xt,1),nsamp);
  end
  for i1 = 1:ngp
    nsampi=sum(gi==i1);
    if nsampi>0
      if nargout<2
        sampft(:,gi==i1) = gp_rnd(gp{i1}, x, y, xt, 'nsamp', nsampi, ...
                         'z', z, 'zt', zt, 'predcf', predcf, ...
                         'tstind', tstind, 'fcorr', fcorr);
      else
        [tsampft, tsampyt] = gp_rnd(gp{i1}, x, y, xt, 'nsamp', nsampi, ...
                                    'z', z, 'zt', zt, 'predcf', predcf, ...
                                    'tstind', tstind, 'fcorr', fcorr);
        sampft(:,gi==i1) = tsampft;
        sampyt(:,gi==i1) = tsampyt;
      end
    end
  end
end

author: gpstuff-dev

repo: gpstuff

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


text:
% fits quadratic to y = f(x) and outputs yq = f(xq)
% x = nsamps x 1
% y = nsamps x nt
% xq = nsampsq x 1
% yq = nsampsq x nt
function yq = fitQuad(x,y,xq)

X    = [ones(numel(x),1) x x.^2];
XT   = (X'*X)\X';

a    = XT * y;

XQ   = [ones(numel(xq),1) xq xq.^2];
yq   = XQ * a;

author: cortex-lab


text:
%Image descriptor based on Histogram of Orientated Gradients for gray-level images. This code
%was developed for the work: O. Ludwig, D. Delgado, V. Goncalves, and U. Nunes, 'Trainable
%Classifier-Fusion Schemes: An Application To Pedestrian Detection,' In: 12th International IEEE
%Conference On Intelligent Transportation Systems, 2009, St. Louis, 2009. V. 1. P. 432-437. In
%case of publication with this code, please cite the paper above.

function H=HOG(Im)
nwin_x=3;%set here the number of HOG windows per bound box
nwin_y=3;
B=9;%set here the number of histogram bins
[L,C]=size(Im); % L num of lines ; C num of columns
H=zeros(nwin_x*nwin_y*B,1); % column vector with zeros
m=sqrt(L/2);
if C==1 % if num of columns==1
    Im=im_recover(Im,m,2*m);%verify the size of image, e.g. 25x50
    L=2*m;
    C=m;
end
Im=double(Im);
step_x=floor(C/(nwin_x+1));
step_y=floor(L/(nwin_y+1));
cont=0;
hx = [-1,0,1];
hy = -hx';
grad_xr = imfilter(double(Im),hx);
grad_yu = imfilter(double(Im),hy);
angles=atan2(grad_yu,grad_xr);
magnit=((grad_yu.^2)+(grad_xr.^2)).^.5;
for n=0:nwin_y-1
    for m=0:nwin_x-1
        cont=cont+1;
        angles2=angles(n*step_y+1:(n+2)*step_y,m*step_x+1:(m+2)*step_x);
        magnit2=magnit(n*step_y+1:(n+2)*step_y,m*step_x+1:(m+2)*step_x);
        v_angles=angles2(:);
        v_magnit=magnit2(:);
        K=max(size(v_angles));
        %assembling the histogram with 9 bins (range of 20 degrees per bin)
        bin=0;
        H2=zeros(B,1);
        for ang_lim=-pi+2*pi/B:2*pi/B:pi
            bin=bin+1;
            for k=1:K
                if v_angles(k)<ang_lim
                    v_angles(k)=100;
                    H2(bin)=H2(bin)+v_magnit(k);
                end
            end
        end

        H2=H2/(norm(H2)+0.01);
        H((cont-1)*B+1:cont*B,1)=H2;
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function [x,fval,exitflag,info,Opt] = opti_fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,opts)
%OPTI_FMINCON Solve a NLP using an OPTI NLP Solver (Matlab Overload)
%
%   [x,fval,exitflag,info] = opti_fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)
%   solves the constrained nonlinear optimization min f(x) where fun is
%   the nonlinear function to be minimized [fun(x)], starting at x0. A,b
%   are linear inequality constraints, Aeq,beq are linear equality
%   constraints, lb,ub are decision variable bounds and nonlcon are the
%   nonlinear constraints in Matlab form.
%
%   [x,fval,exitflag,info] = opti_fmincon(fun,...,opts) allows the
%   user to specify optiset options. This includes specifying a solver via
%   the 'solver' field of optiset.
%
%   [x,...,info,Opt] = opti_fmincon(fun,...) returns the internally
%   built OPTI object.

%   Copyright (C) 2011 Jonathan Currie (IPL)

%Handle missing arguments
if nargin < 10, opts = optiset; end
if nargin < 9, nonlcon = []; end
if nargin < 8, ub = []; end
if nargin < 7, lb = []; end
if nargin < 6, beq = []; end
if nargin < 5, Aeq = []; end
if nargin < 4, b = []; end
if nargin < 3, A = []; end
if nargin < 2, error('You must supply at least 2 arguments to opti_fmincon'); end

%Opti enforces a column for x0, so make sure it is here too
if(size(x0,2) > 1), x0 = x0(:); end

%Sort out Fun + Grad
[f,g] = detGrad(fun,x0);

%Sort out NLCON
[nlcon,nlrhs,nle] = detNlcon(nonlcon,x0);

%Build OPTI Object
Opt = opti('fun',f,'grad',g,'nlmix',nlcon,nlrhs,nle,'ineq',A,b,'eq',Aeq,beq,'bounds',lb,ub,'x0',x0,'options',opts);

%Solve
[x,fval,exitflag,info] = solve(Opt);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function setmap(s, alternateMap)
%SETMAP   sets the default colormap used with spectrograms
%   USAGE: setmap(spectralobject, alternateMap)
%
%   EXAMPLE: setmap(s, bone(64));
%     this sets the map to the 'bone' colormap, with 64 shading  values.
%     If you don't specify the # of values, MatLab will assign the colormap
%     to your current figure.  If no figure exists, MatLab will create one.
%
%   GETMAP and SETMAP are different from GET and SET in that they affect
%   all spectrograms, regardless of the spectralobject used.
%
%   Please Type HELP GRAPH3D to see a number of useful colormaps.

% VERSION: 1.0 of spectralobject
% AUTHOR: Celso Reyes
% LASTUPDATE: 1/30/2007

global SPECTRAL_MAP
%size(SPECTRAL_MAP)
sizecheck = size(alternateMap);
if sizecheck(2) ~= 3,
    warning('Map ignored.  alternateMap has  incorrect # columns; should be 3');
    return
end

SPECTRAL_MAP = alternateMap;

author: geoscience-community-codes

repo: GISMO


text:
function [Population,RankSolution] = EnvironmentalSelection_noCon(Population,N,alpha,gamma,para)

%------------------------------- 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".
%--------------------------------------------------------------------------

% This function is written by Mengjun Ming

    %% Parameter
    popSize   = length(Population);
    NumSeq    = (1:popSize)';
    RankConvg = zeros(popSize,1);
    RankDivs  = zeros(popSize,1);

    %% Modify the infeasible solutions
    PopObj = Population.objs;
    PopCon = Population.cons;
    z      = min(PopObj,[],1);
    n      = max(PopObj,[],1);

    Infeasible_all = any(PopCon>0,2);
    phi_max = max(sum(max(0,PopCon(Infeasible_all,:)),2));

    M          = length(z);
    [W,~]      = UniformPoint(N,M);
    [~,Region] = min(pdist2(PopObj-z,W,'cosine'),[],2);
    PopObj_2   = PopObj;
    for i = 1:size(W,1)
        index = find(Region==i);
        if (~isempty(index))
            Objs_temp  = PopObj_2(index,:);
            Cons_temp  = PopCon(index,:);
            Infeasible = any(Cons_temp>0,2);
            if (sum(Infeasible)~=0)
                F_max = max(Objs_temp,[],1);
                PopObj_2(index(Infeasible),:) = Objs_temp(Infeasible,:)+(sum(max(0,Cons_temp(Infeasible,:)),2)/phi_max).^(exp(para)/max(gamma,0.000001)).*(F_max-Objs_temp(Infeasible,:));
            end
        end
    end

    %% Non-dominated sorting
    Dominate = false(popSize);
    for i = 1 : popSize-1
        for j = i+1 : popSize
            k = any(PopObj_2(i,:)<PopObj_2(j,:)) - any(PopObj_2(i,:)>PopObj_2(j,:));
            if k == 1
                Dominate(i,j) = true;
            elseif k == -1
                Dominate(j,i) = true;
            end
        end
    end

    %% Calculate S(i)
    S = sum(Dominate,2);

    %% Calculate R(i)
    R = zeros(1,popSize);
    for i = 1 : popSize
        R(i) = sum(S(Dominate(:,i)));
    end
    FrontNo = R + 1;

    %% Calculate the crowding distance of each solution
    PopObj   = Population.objs;
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Distance = sort(Distance,2);
    CrowdDis = Distance(:,floor(sqrt(popSize)));

    %% Add a middle column
    MiddleLevel = zeros(popSize,1);

    %% Environmental selection
    Next = FrontNo == 1;
    if sum(Next) <= N
        [~,indx_Convg] = sortrows([FrontNo',-CrowdDis]);
    elseif sum(Next) > N
        Del  = Truncation(Population(Next).objs,sum(Next)-N);
        Temp = find(Next);
        Next(Temp(Del)) = false;
        MiddleLevel(Temp(Del)) = 1;
        [~,indx_Convg] = sortrows([FrontNo',MiddleLevel,-CrowdDis]);
    end
    RankConvg(indx_Convg) = NumSeq;

    %% Environmental selection -- diversity
    FrontNo_D = ones(popSize,1);
    for i = 1:size(W,1)
        index = find(Region==i);
        if (~isempty(index))
            Objs_temp = PopObj_2(index,:);
            g_temp = sum((Objs_temp-z).*W(i,:),2);
            [~,index_FrontNo_D] = sort(g_temp);
            FrontNo_D(index(index_FrontNo_D)) = (1:length(g_temp))';
        end
    end
    [~,indx_divs] = sortrows([FrontNo_D,-CrowdDis]);
    RankDivs(indx_divs) = NumSeq;

    %% Population for next generation
    RankSolution = alpha*RankConvg+(1-alpha)*RankDivs;
    [~,Rank]     = sort(RankSolution);
    Population   = Population(Rank(1:N));
    RankSolution = 1 : N;
end

function Del = Truncation(PopObj,K)
% Select part of the solutions by truncation

    %% Truncation
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Del = false(1,size(PopObj,1));
    while sum(Del) < K
        Remain   = find(~Del);
        Temp     = sort(Distance(Remain,Remain),2);
        [~,Rank] = sortrows(Temp);
        Del(Remain(Rank(1))) = true;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function MW = perform_jp2_rescaling(MW,Jmin,dir)

n = size(MW,1);
Jmax = log2(n)-1;

for j=Jmax:-1:Jmin
    q_min = 1;
    if j==Jmin
        q_min = 0;
    end
    for q=q_min:3
        [selx,sely] = compute_quadrant_selection(j,q);
        MW(selx,sely) = 2^(-(Jmax-j+1)*dir) * MW(selx,sely);
    end
end

author: gpeyre

repo: matlab-toolboxes


text:
function p = fig2plotly(varargin)

%------------------------------FIG2PLOTLY---------------------------------%

% Convert a MATLAB figure to a Plotly Figure

% [CALL]:

% p = fig2plotly
% p = fig2plotly(fig_han)
% p = fig2plotly(fig_han, 'property', value, ...)

% [INPUTS]: [TYPE]{default} - description/'options'

% fig_han: [handle]{gcf} - figure handle
% fig_struct: [structure array]{get(gcf)} - figure handle structure array

% [VALID PROPERTIES / VALUES]:

% filename: [string]{'untitled'} - filename as appears on Plotly
% fileopt: [string]{'new'} - 'new, overwrite, extend, append'
% world_readable: [boolean]{true} - public(true) / private(false)
% link: [boolean]{true} - show hyperlink (true) / no hyperlink (false)
% open: [boolean]{true} - open plot in browser (true)

% [OUTPUT]:

% p - plotlyfig object

% [ADDITIONAL RESOURCES]:

% For full documentation and examples, see https://plot.ly/matlab

%-------------------------------------------------------------------------%

%--FIGURE INITIALIZATION--%
if nargin == 0
    varargin{1} = gcf;
end

%-------------------------------------------------------------------------%

%--CONSTRUCT PLOTLY FIGURE OBJECT--%
p = plotlyfig(varargin{:});

%-------------------------------------------------------------------------%

%--MAKE CALL TO PLOTLY--%
p.plotly;

%-------------------------------------------------------------------------%

end

author: plotly


text:
function quenchWriteFibersFromState(stateFile,fg,fgBins,outName,outType)
%
% quenchWriteFibersFromState([stateFile],[fg],[fgBins],[outName],[outType='pdb'])
%
% OVERVIEW: Take a quench state file and from that file read the unique
%           subgroup indices and save out a fiber group (pdb) for each
%           unique value.
%
% NOTE:
%           This function will be parsing out a text file. If the text file
%           has been altered in any way this will fail spectacularly.
%
% INPUTS:
%       stateFile - a Quench stateFile. -- if there is one qst file in the
%                   directory where this function is called then it
%                   will be selected and used. If there are multiple then
%                   it will prompt the user if it's not provided.
%
%       fg        - a pbd fiber group associated with the qst file. This
%                   will be read in from the qst file. If the fg is not
%                   within the same directory then the user will be
%                   prompted if not provided.
%
%       fgbins    - The bins in the qst for each of the fiber groups user
%                   wants pulled out of the state file. Defaults to all.
%
%       outName   - A [numel(fgBins) X 1] cell array containing the pdb
%                   filenames. If not provided the names will be taken from
%                   the qst file by default.
%
%       outType   - Type of file to save out: 'pdb' or 'mat'
%
% WEB RESOURCES:
%       mrvBrowseSVN('quenchWriteFibersFromState');
%
%
% EXAMPLE USAGE:
%
%   1. Default usage where the '.qst' and '.pdb' file are in the same
%      directory and only one '.qst' file exists.
%       >> quenchWriteFibersFromState;
%
%   2. Same as 1 above, but return only bins 1 4 8 from the 'qst' file and
%      retain the names used within the 'qst' file:
%       >> quenchWriteFibersFromState([],[],[1 4 8],[],[]);
%
%   3. Same as 2, but give each group a specific name:
%       >> quenchWriteFibersFromState([],[],[1 4 8],{'fg_1','fg_2','fg_3'},[]);
%
%   4. Provide all inputs:
%       >> sateFile = 'greatFileName.qst';
%       >> fg = 'greatFileName.pdb';
%       >> fgBins = [1 12 29];
%       >> outName = {'fg_1','fg_2','fg_3');
%       >> outType = 'pdb';
%       >> quenchWriteFibersFromState(stateFile,fg,fgBins,outName,outType);
%
%
% (C) Stanford VISTA, 2012 [lmp]
%

%% Check inputs

% QST file: If there is only one qst file in the CWD then we load that one
% If there is more than one prompt user for the file
if notDefined('stateFile')
    qstName = dir('*qst');
    s = size(qstName);
    if s(1) > 1
       disp('More that one Quench state file found in this directory... Opening dialog...');
        stateFile = mrvSelectFile('r','*.qst','Selcte QUENCH state file',pwd);
    else
       stateFile = fullfile(pwd,qstName.name);
    end
end

% Check that fg is a struct or a file that we can read in
% Perhaps we should read this in later, when we need it.
if notDefined('fg')
    fg = [];
end

% Set outname to empy if they didn't pass it in - we'll set it later
if notDefined('outName')
    outName = '';
end

% They didn't pass in fGroups so we set it to empty and use all unique
% values later
if notDefined('fgBins')
   fgBins = [];
end

% By default we save out pdb files
if notDefined('outType')
    outType = 'pdb';
end

%% Read the qst file

% Open the state file for reading
qst = fopen(stateFile,'r');

% Find out how many lines are in the qst text file
num = 1;
tline = fgetl(qst);
while ischar(tline)
   tline = fgetl(qst);
   num   = num+1;
end
fclose(qst);

% Open the file back up for reading
qst = fopen(stateFile,'r');

% Loop over the file and read each line individually into a cell array - Q
Q = cell(1,num);
for ii = 1:num
   tline = fgetl(qst);
   Q{ii} = tline;
end

fclose(qst);

% Q is now a cell array containing each line of the qst file
Q = Q';

%% Load the fiber group and begin to parse the Quench State .qst file

% If the user did not pass in a fiber group then 'fg' will be empty and we
% will try to load the fiber group from the path listed within the quench
% state file
if isempty(fg)
    disp('Loading fibers from path in the ''qst'' file...');
    fg = Q{14}; % The path in the qst file
    try
        fg = fgRead(fg);
    catch ME
        % That didn't work so we prompt the user to select the fiber group
        disp(ME);
        fg = mrvSelectFile('r','*.pdb','Select PDB file',pwd);
        fg = fgRead(fg);
    end

% It's not empty and it's a file - load it.
elseif ~isstruct(fg) && exist(fg,'file')
    disp('Loading fibers from file ...');
    fg = fgRead(fg);

% It was passed in as a structure - so we're done
elseif isfield(fg,'fibers')
    disp('Fibers already loaded.');
end

% Sanity check: Check the number of inds against the number of fibers in
% fg.fibers
if str2num(Q{18}(15:end)) ~= numel(fg.fibers)
    error('The number of fibers referenced in the ''qst'' file does not match that of the fiber group passed in!')
end

% The fiber indices from the qst file
fInds = str2num(Q{19});

% The unique fiber group indices
fGroups = unique(fInds);

% The total number of fiber groups
nGroups = numel(fGroups);

%% Get names for the fiber groups

% Initialize the fName cell array
fName = cell(nGroups,1);

% Set up the fName cell array for naming each fiber group
if notDefined('outName') || numel(outName) ~= numel(fgBins)
    % Starting number of lines before first name. Names occur every 5 lines.
    nind = 28;
    for nn = 1:nGroups
        % Remove any spaces in the file names and replace with '_'
        fName{nn} = regexprep(Q{nind}(7:end),'\W','_');
        % Increase the line ind by 5 each pass.
        nind = nind+5;
    end
else
    fName = outName;
end

% If outName was not passed in but fgBins are reuested then take those
% fgBins out of the fNames structure ans use only those for the names
if notDefined('outName') && ~notDefined('fgBins');
   outName = fName(fgBins);
   fName   = outName;
end

%% Create the fiber groups and write them out

fprintf('Creating %s groups from ''qst'' file...',num2str(numel(fGroups)));

% The ID of unique fiber group bins - if passed in then use those
if ~notDefined('fgBins')
   fGroups = fgBins;
end

% For each unique fg index in fgGroups assign the fibers to a new group and
% write them out
for ff = 1:numel(fGroups)

    % fgI are the indices for a given fibergroup (ff) - fInds is the long
    % array of indices from the state file
    fgI = find(fInds == fGroups(ff));

    % Keep only fgI from fg.fibers
    fgOut = fgExtract(fg,fgI,'keep');

    % Set the name from the fName cell array
    fgOut.name = fName{ff};

    % Show the user what's going on
    fprintf('\n  %s - %s',num2str(ff), fName{ff});

    % Save out the fiber group
    fgWrite(fgOut,fgOut.name,outType);

    % Clear for memory
    clear fgOut
end

fprintf('\nDone.\n');

return

%%
%#ok<*ST2NM>

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function sVF = uminus(sVF)
%
% Syntax
%   sVF = -sVF
%

sVF.sF = -sVF.sF;

end

author: mtex-toolbox

repo: mtex

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


text:
function toms097_test02 ( )

%*****************************************************************************80
%
%% TOMS097_TEST02 tests R8MAT_SHORTEST_PATH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 March 2014
%
%  Author:
%
%    John Burkardt
%
  n = 6;

  a = [ ...
     0.0, -1.0, -1.0, -1.0, -1.0, -1.0; ...
     2.0,  0.0, -1.0, -1.0, -1.0,  5.0; ...
     5.0,  7.0,  0.0, -1.0,  2.0, -1.0; ...
    -1.0,  1.0,  4.0,  0.0, -1.0,  2.0; ...
    -1.0, -1.0, -1.0,  3.0,  0.0,  4.0; ...
    -1.0,  8.0, -1.0, -1.0,  3.0,  0.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TOMS097_TEST02\n' );
  fprintf ( 1, '  R8MAT_SHORTEST_PATH uses Floyd''s algorithm to find the\n' );
  fprintf ( 1, '  shortest distance between all pairs of nodes\n' );
  fprintf ( 1, '  in a directed graph, starting from the initial array\n' );
  fprintf ( 1, '  of direct node-to-node distances.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  In the initial direct distance array, if\n' );
  fprintf ( 1, '    A(I,J) = -1,\n' );
  fprintf ( 1, '  this indicates there is NO directed link from\n' );
  fprintf ( 1, '  node I to node J.  In that case, the value of\n' );
  fprintf ( 1, '  of A(I,J) is essentially "infinity".\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial direct-link distance matrix:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '%10.4f', a(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  for j = 1 : n
    for i = 1 : n
      if ( a(i,j) == -1.0 )
        a(i,j) = Inf;
      end
    end
  end

  a = r8mat_shortest_path ( n, a );

  for j = 1 : n
    for i = 1 : n
      if ( a(i,j) == Inf )
        a(i,j) = -1.0;
      end
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  In the final shortest distance array, if\n' );
  fprintf ( 1, '    A(I,J) = -1,\n' );
  fprintf ( 1, '  this indicates there is NO directed path from\n' );
  fprintf ( 1, '  node I to node J.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Final distance matrix:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '%10.4f', a(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
%% Visualization with MATLAB Builder for .NET
% Copyright (c) 2009, The MathWorks, Inc.
%
% This demo / example shows how to implement an C# application that uses
% MATLAB built components as a visualization engine.
% The main window can be resized, and there is a movable slitter between
% the textbox and the picturebox.
% After the functions are called for the first time, visualization
% should be fast enough to enable a real-time experience when
% resizing the figure.

% There is a compiled exe and related DLLs in the folder bin\release

%%
% 1. Compile the dll with
mcc  -B 'dotnet:NETSquare,NETSquareclass,2.0,Private,local' makesquare.m ...
    makebitmap.m makefigure.m

% 2. Add generated DLL and MWArray DLL to project references

% 3. Compile and run the project


function combine_box_region_nets(box_net, region_net, pi_net)
if(~exist(box_net, 'file'))
    fprintf('Box net does not exist!\n'); return;
end
if(~exist(region_net, 'file'))
    fprintf('Region net does not exist!\n');return;
end
caffe('init', 'prototxts/rcnn_extract_fc7.prototxt', box_net);
box_layers=caffe('get_weights');
caffe('init', 'prototxts/rcnn_extract_fc7.prototxt', region_net);
region_layers=caffe('get_weights');

for i=1:numel(box_layers)
    box_layers(i).layer_names=[box_layers(i).layer_names '_1'];
end
for i=1:numel(box_layers)
    region_layers(i).layer_names=[region_layers(i).layer_names '_2'];
end
layers=[box_layers; region_layers];
caffe('read_prototxt', 'prototxts/pinetwork_extract_fc7.prototxt');
caffe('change_weights', layers);
caffe('save_net', pi_net);


text:
%-------------------------------------------------------------------------------------------------------------
% This is an implementation of the TWSC algorithm for additive white Gaussian noise
% noise removal.
%
% Author:  Jun Xu, csjunxu@comp.polyu.edu.hk / nankaimathxujun@gmail.com
%          The Hong Kong Polytechnic University
%
% Please refer to the following paper if you find this code helps:
%
% @article{TWSC_ECCV2018,
% 	author = {Jun Xu and Lei Zhang and David Zhang},
% 	title = {A Trilateral Weighted Sparse Coding Scheme for Real-World Image Denoising},
% 	journal = {ECCV},
% 	year = {2018}
% }
%
% Please see the file License.txt for the license governing this code.
%-------------------------------------------------------------------------------------------------------------
clear;
Original_image_dir  =    'cleanimages/'; % put clean images in this folder
Sdir = regexp(Original_image_dir, '\', 'split');
fpath = fullfile(Original_image_dir, '*.png');
im_dir  = dir(fpath);
im_num = length(im_dir);

method = 'TWSC';
write_MAT_dir = [method '/Results_AWGN/'];
write_sRGB_dir  = [write_MAT_dir method '/'];
if ~isdir(write_sRGB_dir)
    mkdir(write_sRGB_dir);
end
for nSig = [15 25 35 50 75]
    %% Parameters
    Par.innerIter = 2;
    Par.win = 30;
    Par.lambda1 = 0;
    Par.ps = 8;
    Par.outerIter = 10;
    Par.step = 3;
    Par.nlspini = 90;
    Par.nlspgap = 10;
    if 0 < nSig <= 20
        Par.outerIter = 8;
        Par.delta = .07;
        Par.nlspini = 70;
        Par.lambda2 = .9;
    elseif 20 < nSig <= 30
        Par.delta = .06;
        Par.lambda2 = .76;
    elseif 30 < nSig <= 40
        Par.delta = .07;
        Par.lambda2 = .78;
    elseif 40 < nSig <= 60
        Par.nlspini = 120;
        Par.nlspgap = 15;
        Par.delta = .05;
        Par.lambda2 = .72;
    elseif 60 < nSig <= 80
        Par.ps = 9;
        Par.outerIter = 14;
        Par.step = 4;
        Par.nlspini = 140;
        Par.delta = .05;
        Par.lambda2 = .68; % .66
    else
        disp('Please tune the above parameters by yourself, thanks!');
    end
    % record all the results in each iteration
    Par.PSNR = zeros(Par.outerIter, im_num, 'double');
    Par.SSIM = zeros(Par.outerIter, im_num, 'double');
    T512 = [];
    T256 = [];
    for i = 1:im_num
        Par.nlsp = Par.nlspini;  % number of non-local patches
        Par.image = i;
        Par.nSig = nSig/255;
        Par.I =  im2double( imread(fullfile(Original_image_dir, im_dir(i).name)) );
        S = regexp(im_dir(i).name, '\.', 'split');
        randn('seed',0);
        Par.nim =   Par.I + Par.nSig*randn(size(Par.I));
        fprintf('%s :\n',im_dir(i).name);
        PSNR =   csnr( Par.nim*255, Par.I*255, 0, 0 );
        SSIM      =  cal_ssim( Par.nim*255, Par.I*255, 0, 0 );
        fprintf('The initial value of PSNR = %2.4f, SSIM = %2.4f \n', PSNR,SSIM);
        time0 = clock;
        [im_out, Par]  =  TWSC_Sigma_AWGN(Par);
        if size(Par.I,1) == 512
            T512 = [T512 etime(clock,time0)];
            fprintf('Total elapsed time = %f s\n', (etime(clock,time0)) );
        elseif size(Par.I,1) ==256
            T256 = [T256 etime(clock,time0)];
            fprintf('Total elapsed time = %f s\n', (etime(clock,time0)) );
        end
        im_out(im_out>1)=1;
        im_out(im_out<0)=0;
        % calculate the PSNR
        Par.PSNR(Par.outerIter, Par.image)  =   csnr( im_out*255, Par.I*255, 0, 0 );
        Par.SSIM(Par.outerIter, Par.image)      =  cal_ssim( im_out*255, Par.I*255, 0, 0 );
        imname = sprintf([write_sRGB_dir method '_nSig' num2str(nSig) '_oIte' num2str(Par.outerIter) '_iIte' num2str(Par.innerIter) '_ps' num2str(Par.ps) '_step' num2str(Par.step) '_nlspini' num2str(Par.nlspini) '_nlspgap' num2str(Par.nlspgap) '_delta' num2str(Par.delta) '_l1' num2str(Par.lambda1) '_l2' num2str(Par.lambda2) '_' im_dir(i).name]);
        imwrite(im_out,imname);
        fprintf('%s : PSNR = %2.4f, SSIM = %2.4f \n',im_dir(i).name, Par.PSNR(Par.outerIter, Par.image),Par.SSIM(Par.outerIter, Par.image)     );
    end
    PSNR = Par.PSNR(end,:);
    mPSNR=mean(PSNR,2);
    SSIM = Par.SSIM(end,:);
    mSSIM=mean(SSIM,2);
    mT512 = mean(T512);
    sT512 = std(T512);
    mT256 = mean(T256);
    sT256 = std(T256);
    fprintf('The average PSNR = %2.4f, SSIM = %2.4f. \n', mPSNR,mSSIM);
    name = sprintf([write_MAT_dir method '_nSig' num2str(nSig) '_oIte' num2str(Par.outerIter) '_iIte' num2str(Par.innerIter) '_ps' num2str(Par.ps) '_step' num2str(Par.step) '_nlspini' num2str(Par.nlspini) '_nlspgap' num2str(Par.nlspgap) '_delta' num2str(Par.delta) '_l1' num2str(Par.lambda1) '_l2' num2str(Par.lambda2) '.mat']);
    save(name,'nSig','PSNR','SSIM','mPSNR','mSSIM','mT512','sT512','mT256','sT256');
end

author: csjunxu


text:
classdef chebguiExporter
%CHEBGUIEXPORTER   Class to export a problem from CHEBGUI.
%   This is a an abstract class that defines interfaces for export problems from
%   CHEBGUI to .m-files, to the workspace, or to a .chebgui file. It is not
%   intended to be called directly by the end user.
%
% See also CHEBGUI, CHEBGUIBVPEXPORTER, CHEBGUIEIGEXPORTER, CHEBGUIPDEEXPOERTER.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEVELOPER NOTE:
%   The CHEBGUIEXPORTER class defines a number of abstract methods, used to
%   export problems from CHEBGUI. In v4, this functionality used to live in the
%   @chebgui folder, but to increase modularity, it has been spun off to its own
%   class and subclasses.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% ABSTRACT METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties ( Abstract = true )

        % Default file name to be saved to
        defaultFileName

        % Description of exported .m files
        description

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% ABSTRACT STATIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Abstract = true, Static = true )

        % Extract information from the CHEBGUI object to a struct:
        expInfo = exportInfo(guifile)

        % Print problem description:
        printDescription(fid, expInfo)

        % Print options for solving the problem:
        printOptions(fid, expInfo)

        % Print lines for setting up the problem:
        printSetup(fid, expInfo, guifile)

        % Print the actual lines for calling the solver method:
        printSolver(fid, expInfo)

        % Print steps taken after the solver finishes:
        printPostSolver(fid, expInfo)

        % Export to a .mat-file
        toMat(handles)

        % Export to workspace
        toWorkspace(handles)

        % Export to workspace, solution only
        toWorkspaceSolutionOnly(handles)

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CONCRETE METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function toFile(exporter, guifile, fileName, pathName)
        %TOFILE   Export a CHEBGUI to an m-file
        %   TOFILE(EXPORTER, GUIFILE, FILENAME, PATHNAME) exports the CHEBGUI
        %   object GUIFILE to the file 'PATHNAME/FILENAME' using the
        %   CHEBGUIEXPORTER object EXPORTER (of BVP, EIG, or PDE kind).

        % Concatenate the pathName and the fileName to get the full path:
        fullFileName = fullfile(pathName, fileName);

        try
            % Open a stream to write to a file:
            fid = fopen(fullFileName, 'wt');

            % Extract the necessary info for export to an .m file from the
            % GUIFILE object:
            expInfo = exporter.exportInfo(guifile);

            % Write the header information:
            writeHeader(exporter, fid, fileName)

            % Print description of the problem:
            exporter.printDescription(fid, expInfo)

            % Print the problem set-up:
            exporter.printSetup(fid, expInfo, guifile)

            % Print the options set-up:
            exporter.printOptions(fid, expInfo)

            % Print the solution step:
            exporter.printSolver(fid, expInfo)

            % Print the post-solution process:
            exporter.printPostSolver(fid, expInfo)

            % Close the file writing stream:
            fclose(fid);
        catch ME
            % Make sure to tidy up first
            fclose(fid);

            % Rethrow the error:
            rethrow(ME)
        end

        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CONCRETE STATIC METHODS
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = true )
        function obj = constructor(type)
        %CONSTRUCTOR   Constructor for the CHEBGUIEXPORTER class.
        %   OBJ = CONSTRUCTOR(TYPE), where TYPE = 'bvp', 'eig' or 'pde' returns
        %   an object that is a concrete implementation of CHEBGUIEXPORTER
        %   object OBJ. The type of OBJ is determined by the TYPE input
        %   argument.
        %
        %   This method allows one to construct concrete subclasses of the
        %   CHEBGUIEXPORTER class.

            % Return the desired type of CHEBGUIEXPORTER
            switch type
                case 'bvp'
                    obj = chebguiExporterBVP;

                case 'eig'
                    obj = chebguiExporterEIG;

                case 'ivp'
                    obj = chebguiExporterIVP;

                case 'pde'
                    obj = chebguiExporterPDE;

                otherwise
                    error('CHEBFUN:CHEBGUIEXPORTER:chebguiExporter:constructor', ...
                        'Unknown type for CHEBGUIEXPORTER constructor.')
            end

        end

        function toChebgui(guifile)
        %TOCHEBGUI   Export CHEBGUI to a CHEBGUI object in the workspace.

            % Setup for input dialog
            prompt = 'Enter the name of the chebgui variable:';
            name = 'Export GUI';
            numlines = 1;
            defaultAnswer ='chebg';
            options.Resize ='on';
            options.WindowStyle ='modal';

            % Get user input
            answer = inputdlg(prompt, name, numlines, {defaultAnswer}, options);

            % If user did not press cancel, assign to a variable in the base
            % workspace.
            if ( ~isempty(answer) )
                assignin('base', answer{1}, guifile);
            end
        end

        function str = prettyPrintFevalString(str, varnames)
        %PRETTYPRINTFEVALSTRING   Tidy up strings printed to .m files.

            % Do various string manipulations to make the output look nicer.
            for k = 1:numel(varnames)
                oldstr = ['feval(' varnames{k} ','];
                newstr = [varnames{k} '('];
                str = strrep(str, oldstr, newstr);
                oldstr = [varnames{k} '(''end'''];
                newstr = [varnames{k} '(end'];
                str = strrep(str, oldstr, newstr);
                oldstr = [varnames{k} '(''right'''];
                newstr = [varnames{k} '(end'];
                str = strrep(str, oldstr, newstr);
                oldstr = [varnames{k} '(''start'''];
                newstr = [varnames{k} '(' varnames{k} '.ends(1)'];
                str = strrep(str, oldstr, newstr);
                oldstr = [varnames{k} '(''left'''];
                newstr = [varnames{k} '(' varnames{k} '.ends(1)'];
                str = strrep(str, oldstr, newstr);
            end

        end

        function disc = discOption(isPeriodic, dom, opt)
            % Do we want to use TRIGCOLLOC for the discretization? This will be
            % the case if we're solving a periodic problem, and have
            % 'collocation' specified as the discretization option. However, we
            % don't want to use TRIGCOLLOC if there are breakpoints in the
            % domain.
            if ( isPeriodic && strcmp(opt, 'collocation') &&  ...
                    ( length(str2num(dom)) == 2 ) )
                disc = 'periodic';
            else
                disc = opt;
            end
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CONCRETE PROTECTED METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = protected )

        function writeHeader(e, fid, filename)
            %WRITEHEADER   The first few lines of exported .m-files.
            %
            % Inputs:
            %   E:          a CHEBGUIEXPORT object
            %   FID:        a valid file ID
            %   FILENAME:   name of the file being written to

            % Access usernames in different ways, depending on whether are
            % running on PCs or not:
            if ( ispc )
                userName = getenv('UserName');
            else
                userName = getenv('USER');
            end

            % Print first few lines of the .m-file:
            fprintf(fid, ['%%%% %s -- an executable m-file for solving ', ...
                '%s\n'], filename, e.description);
            fprintf(fid, ['%% Automatically created in CHEBGUI ', ...
                'by user %s.\n'], userName);
            fprintf(fid, '%% Created on %s', ...
                datestr(now, 'mmmm dd, yyyy at HH:MM.\n\n'));

            fprintf(fid, '%%%% Problem description.\n');

        end

    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function icol = r8col_find ( m, n, a, x )

%*****************************************************************************80
%
%% R8COL_FIND seeks a column value in an R8COL.
%
%  Example:
%
%    Input:
%
%      M = 3,
%      N = 4,
%
%      A = (
%        1.  2.  3.  4.
%        5.  6.  7.  8.
%        9. 10. 11. 12. )
%
%      x = ( 3.,
%            7.,
%           11. )
%
%    Output:
%
%      ICOL = 3
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), a table of numbers, regarded as
%    N columns of vectors of length M.
%
%    Input, real X(M), a vector to be matched with a column of A.
%
%    Output, integer ICOL, the index of the first column of A
%    which exactly matches every entry of X, or -1 if no match
%    could be found.
%
  icol = -1;

  for j = 1 : n

    if ( x(1:m) == a(1:m,j)' )
      icol = j;
      return
    end

  end

  return
end


text:
% GP_COV_PP - Piecewise polynomial covariance function with compact support
% for Gaussian processes.
%
% [K, DK_LOGTHETA, DK_X2] = GP_COV_CS(X1, X2, LOGTHETA)

% Last modified 2010-10-06
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function covfunc = gp_cov_pp_init(x1,x2)

% Dimensions
d = rows(x2); % dimensionality of inputs
m = cols(x1); % number of other inputs
n = cols(x2); % number of inputs

q = 2; % PP parameter, fix to 2 for now..

% Distance matrix
D = sqrt(sq_dist(X1,X2));

covfunc = @set_hyperparameters;

  function [K,L] = set_hyperparameters(logtheta)

  % Threshold
  thres = exp(logtheta(1)); % length scale, threshold
  R = D./thres;
  I = sparse(R<1);
  NNZ = nnz(I); % number of nonzero elements

  j = floor(rows(x1)) + q + 1;
  switch q
    %%%%%%%%%%
   case 0
    error('not yet implemented for q=0')
    K(I) = (1-R(I)) .^ j;
    %%%%%%%%%%
   case 1
    error('not yet implemented for q=1')
    %%%%%%%%%%

   case 2
    a = j^2 + 4*j + 3;
    b = 3*j + 6;
    c = 3;
    z = 1/3;

    % Helper
    R2 = spalloc(m,n,NNZ);
    R2(I) = D2(I) ./ (thres^2);

    % Covariance matrix
    K = spalloc(m,n,NNZ);
    K(I) = z * (1-R(I)).^(j+2) .* (a*R2(I) + b*R(I) + c);

    dK_dR = spalloc(m,n,NNZ); % copy sparseness
    dK_dR(I) = -z * (j+2) * (1-R(I)).^(j+1) .* (a*R2(I) + b*R(I) + c) ...
        + z * (1-R(I)).^(j+2) .* (2*a*R(I) + b);

    % Gradient for hyperparameters
    if nargout >= 2
      % TODO 3D SPARSE MATRICES DO NOT WORK. USE CELL ARRAYS?
      dK_dlogtheta = spalloc(m,n,NNZ);
      dK_dlogtheta(I) = dK_dR(I) .* (-D(I) ./ thres);
      dK_dlogtheta = full(dK_dlogtheta); % blaaah.. :(
    end

    % Gradients for inputs x2
    if nargout >= 3
      if isempty(x2)
        error('Can''t calculate gradient: x2 not given');
      end
      % TODO 3D SPARSE MATRICES DO NOT WORK. USE CELL ARRAYS?
      % TODO: you have dK_dR and dD_dx but NOT dR_dD!!
      dK_dx2 = bsxfun(@times, reshape(full(dK_dR./thres),[1,m,n]), dD_dx2);
      % blaah the need for fullness.. :(
    end

    %%%%%%%%%%

   case 3
    error('not yet implemented for q=3')
    %%%%%%%%%%
   otherwise
    error('q not valid')
  end

  if ~issparse(K)
    error('K not sparse, it should, wtf?!');
  end

  if nargout >= 2
    LD = ldlchol(K);
    L = ldlsplit(LD);
  end
  end

end

x = [];  % data inputs
D = [];  % distance matrix for data
K = [];  % data covariance
LD = []; % LDL for data covariance
hyperparameters = [];

covstruct.set_data_inputs = @set_data_inputs;

  function set_data_inputs(inputs)
  % Evaluate Euclidean distances
  if ~isequal(x,inputs)
    x = inputs;
    D = sqrt(sq_dist(x));
    K = [];
    LD = [];
  end
  end

  function set_hyperparameters(params)
  if ~isequal(hyperparameters,params)
    hyperparameters = params;
    K = [];
    LD = [];
  end
  end

if nargin < 3
  x2 = [];
end
if nargin < 4
  distfun = @sqdistEuclidean;
end
if nargin < 5
  q = 2;
end

% 2-norm distances
if ~isempty(x2)
  if nargout <= 1
    D = feval(distfun, x1, x2, false);
  else
    [D, dD_dx2] = feval(distfun, x1, x2, false);
  end
  D2 = D.^2; % get this from the distfun?
else
  % Distances for variances
  D = spalloc(cols(x1),1,0); %zeros(cols(x1),1);
  D2 = D; %zeros(cols(x1),1);
end

% Dimensions
d = rows(x2); % dimensionality of inputs
m = cols(x1); % number of other inputs
n = cols(x2); % number of inputs

% Threshold
thres = exp(logtheta(1)); % length scale, threshold
R = D./thres;
I = sparse(R<1);
NNZ = nnz(I); % number of nonzero elements
%K = spones(I); % copy sparseness
% $$$ Dsp = spalloc(m,n,NNZ);
% $$$ Dsp(I) = D(I);
% $$$ D2sp = spalloc(m,n,NNZ);
% $$$ D2sp(I) = D2(I);

j = floor(rows(x1)) + q + 1;
switch q
  %%%%%%%%%%
 case 0
  error('not yet implemented for q=0')
  K(I) = (1-R(I)) .^ j;
  %%%%%%%%%%
 case 1
  error('not yet implemented for q=1')
  %%%%%%%%%%

 case 2
  a = j^2 + 4*j + 3;
  b = 3*j + 6;
  c = 3;
  z = 1/3;

  % Helper
  R2 = spalloc(m,n,NNZ);
  R2(I) = D2(I) ./ (thres^2);

  % Covariance matrix
  K = spalloc(m,n,NNZ);
  K(I) = z * (1-R(I)).^(j+2) .* (a*R2(I) + b*R(I) + c);

  dK_dR = spalloc(m,n,NNZ); % copy sparseness
  dK_dR(I) = -z * (j+2) * (1-R(I)).^(j+1) .* (a*R2(I) + b*R(I) + c) ...
      + z * (1-R(I)).^(j+2) .* (2*a*R(I) + b);

  % Gradient for hyperparameters
  if nargout >= 2
    % TODO 3D SPARSE MATRICES DO NOT WORK. USE CELL ARRAYS?
    dK_dlogtheta = spalloc(m,n,NNZ);
    dK_dlogtheta(I) = dK_dR(I) .* (-D(I) ./ thres);
    dK_dlogtheta = full(dK_dlogtheta); % blaaah.. :(
  end

  % Gradients for inputs x2
  if nargout >= 3
    if isempty(x2)
      error('Can''t calculate gradient: x2 not given');
    end
    % TODO 3D SPARSE MATRICES DO NOT WORK. USE CELL ARRAYS?
    % TODO: you have dK_dR and dD_dx but NOT dR_dD!!
    dK_dx2 = bsxfun(@times, reshape(full(dK_dR./thres),[1,m,n]), dD_dx2);
    % blaah the need for fullness.. :(
  end

  %%%%%%%%%%

 case 3
  error('not yet implemented for q=3')
  %%%%%%%%%%
 otherwise
  error('q not valid')
end

if ~issparse(K)
  error('K not sparse, it should, wtf?!');
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% This script defines a project shortcut.
%
% To get a handle to the current project use the following function:
%
% project = simulinkproject();
%
% You can use the fields of project to get information about the currently
% loaded project.
%
% See: help simulinkproject

Simulink.fileGenControl('reset')

close all;
clear all;

author: TUMFTM


text:
% load textbook.mat
S=[-1  0  1 1 0;
    1 -1  0 0 0;
    0  1 -1 0 -1];

model.S=S;
model.lb=[0 0 -10 5 5]';
model.ub=[5 5 0 0 0]';
model.c=zeros(5,1);
model.c(5,1)=1;

v = checkThermodynamicConsistency(model);

disp(norm(v))

author: opencobra

repo: cobratoolbox


text:
function ind = inpolygon(ebsd,xy)
% checks which ebsd data are within given polygon
%
% Syntax
%   ind = inpolygon(ebsd,[xmin,ymin,dx,dy]) % select indices by rectangle
%   ind = inpolygon(ebsd,[x1 y1; x2 y2; x3 y3; x4 y4]) % select indices by poylgon
%   ebsd = ebsd(ind) % select EBSD data by indices
%
% Input
%  ebsd    - @EBSD
%  xmin, xmax - lower left corner of a rectangle
%  dx, dy - extend of a rectangle
%  x, y  - vertices of a polygon
%
% Output
%  ind - logical
%
% See also
% inpolygon

% shortcut for simple rectangles
if numel(xy)==4

  xy(3:4) = xy(1:2) + xy(3:4);
  corners = [1 2; 3 2; 3 4; 1 4; 1 2];
  xy = xy(corners);

end

% check for inside
if ~getMTEXpref('insidepoly')
  ind = inpolygon(ebsd.prop.x,ebsd.prop.y,xy(:,1),xy(:,2));
else
  ind = insidepoly(ebsd.prop.x,ebsd.prop.y,xy(:,1),xy(:,2));
end

end


text:
function d = dot(g1,g2,varargin)
% inner product of quaternions g1 and g2
%
% Input
%  q1, q2 - @quaternion
%
% Output
%  d - double
%

d = g1.a .* g2.a + g1.b .* g2.b + g1.c .* g2.c + g1.d .* g2.d;