text:
function [tinyImage, colorHist, param] = LMcolor(D, HOMEIMAGES, param)
%
% [tinyImage, colorHist, param] = LMcolor(D, HOMEIMAGES, param);
%
% param.tinySize
% param.colorHist.imagesize
% param.colotHist.
%
% Use transformations from "color indexing", Swain & Ballard, IJCV 91
%
%   wb = (r+g+b)/3
%   rg = r - g
%   by = b-r/2-g/2

if nargin<3
    % Default parameters
    param.tinySize = [16 16];
    param.colorHist.imagesize = [256 256];
    param.colorHist.nbins = [8 16 16];
    param.colorHist.margins.wb = [0 255];
    param.colorHist.margins.rg = [-255 255];
    param.colorHist.margins.by = [-255 255];
end

% Precompute filter transfert functions (only need to do this once, unless
% image size is changes):
Nfeatures = prod(param.colorHist.nbins);

if isstruct(D)
    % [gist, param] = LMcolor(D, HOMEIMAGES, param);
    Nscenes = length(D);
    typeD = 1;
end
if iscell(D)
    % [gist, param] = LMcolor(filename, HOMEIMAGES, param);
    Nscenes = length(D);
    typeD = 2;
end
if isnumeric(D)
    % [gist, param] = LMcolor(img, HOMEIMAGES, param);
    Nscenes = size(D,4);
    typeD = 3;
end

% Loop: Compute gist features for all scenes
tinyImage = zeros([param.tinySize 3 Nscenes], 'single');
colorHist = zeros([Nscenes Nfeatures], 'single');
for n = 1:Nscenes
    disp([n Nscenes])

    % load image
    try
        switch typeD
            case 1
                img = LMimread(D, n, HOMEIMAGES);
            case 2
                img = imread(fullfile(HOMEIMAGES, D{n}));
            case 3
                img = D(:,:,:,n);
        end
    catch
        disp(D(n).annotation.folder)
        disp(D(n).annotation.filename)
        rethrow(lasterror)
    end

    % resize and crop image to make it square
    tiny = imresizecrop(img, param.tinySize+2, 'bilinear');
    tinyImage(:,:,:,n) = tiny(2:end-1, 2:end-1,:);

    subplot(211)
    imagesc(img); axis('off'); axis('equal')
    subplot(212)
    imagesc(tiny); axis('on'); axis('equal')

    % resize and crop image to make it square
    img = single(imresizecrop(img, param.imagesize, 'bilinear'));

    %   wb = (r+g+b)/3
    %   rg = r - g
    %   by = b-r/2-g/2
    wb = mean(img,3);
    rg = img(:,:,1)-img(:,:,2);
    by = img(:,:,3)-img(:,:,1)/2-img(:,:,2)/2;

    wb = (wb(:)-param.colorHist.margins.wb(1))/(param.colorHist.margins.wb(2)-param.colorHist.margins.wb(1));
    rg = (rg(:)-param.colorHist.margins.rg(1))/(param.colorHist.margins.rg(2)-param.colorHist.margins.rg(1));
    by = (by(:)-param.colorHist.margins.by(1))/(param.colorHist.margins.by(2)-param.colorHist.margins.by(1));

    wb = fix(wb*param.colorHist.nbins(1));
    rg = fix(rg*param.colorHist.nbins(2));
    by = fix(by*param.colorHist.nbins(3));

    wb = min(max(0,wb), param.colorHist.nbins(1)-1);
    rg = min(max(0,rg), param.colorHist.nbins(2)-1);
    by = min(max(0,by), param.colorHist.nbins(3)-1);

    h = by + rg*param.colorHist.nbins(3) + wb*param.colorHist.nbins(3)*param.colorHist.nbins(2);
    H = hist(h, [0:Nfeatures-1]);

    % store
    colorHist(n,:) = H;
    drawnow
end

author: CSAILVision

repo: LabelMeToolbox


text:
function L=intesselation(X,TES,XI)

t = tsearchn(X,TES,XI);
L=~isnan(t);

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 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, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%% stochasticMicrostructure
% Below is a demonstration of the features of the |stochasticMicrostructure| function

%%
clear; close all; clc;

%% Syntax
% |[F,V,C]=stochasticMicrostructure(inputStruct);|

%% Description
% This function generates Stochastic Bicontinuous Microstructures
%
% Input structure and default values:
%
%   inputStruct.L=1; % characteristic length
%   inputStruct.Ns=80; % number of sampling points
%   inputStruct.Nw=120; % number of waves
%   inputStruct.q0=55; % wave number
%   inputStruct.relD=0.5; % relative density
%   inputStruct.anisotropyFactors=[1 1 1]; %Anisotropy factors
%   inputStruct.isocap=1; %Option to cap the isosurface
%
% Based on: Soyarslan et al. "3D stochastic bicontinuous
% microstructures: Generation, topology and elasticity"
% https://doi.org/10.1016/j.actamat.2018.01.005
%
% Original author: Sebastien Callens, September 2020

%% Examples

%%
% Plot settings
cMap=parula(250);
faceAlpha1=1;
faceAlpha2=0.5;
edgeColor1='none';
edgeColor2='none';
fontSize=15;

%% Example 1
inputStruct.L=1; % characteristic length
inputStruct.Ns=75; % number of sampling points
inputStruct.Nw=500; % number of waves
inputStruct.q0=55; % wave number
inputStruct.relD=0.55; % relative density
inputStruct.anisotropyFactors=[1 1 1]; %Anisotropy factors
inputStruct.isocap=1; %Option to cap the isosurface

%%
% Create stochastic structure

[F,V,C]=stochasticMicrostructure(inputStruct);

%%
% Using grouping to keep only largest group
groupOptStruct.outputType='label';
[G,~,groupSize]=tesgroup(F,groupOptStruct); %Group connected faces
[~,indKeep]=max(groupSize); %Index of largest group

%Keep only largest group
F=F(G==indKeep,:); %Trim faces
C=C(G==indKeep,:); %Trim color data
[F,V]=patchCleanUnused(F,V); %Remove unused nodes

%%
% Visualize surface

cFigure;
gpatch(F,V,C,'none');
axisGeom; camlight headlight;
colormap gjet; icolorbar;
gdrawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 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, see <http://www.gnu.org/licenses/>.


text:
function out = traj_cost_lgr(opt_vars,traj_par,baseQR)
% -------------------------------------------------------------------
% This function computes cost in terms of condition number for
% trajectory optimization needed for dynamic parameter identification
% The computation of regressor matrix is obtained using screw
% theory methods.
% -------------------------------------------------------------------

% Trajectory parameters
N = traj_par.N;
wf = traj_par.wf;
T = traj_par.T;
t = traj_par.t;

% As paramters of the trajectory are in a signle vector we reshape them as
% to feed the function that computes the trajectory
ab = reshape(opt_vars,[12,N]);
a = ab(1:6,:); % sin coeffs
b = ab(7:12,:); % cos coeffs

% To guarantee that positions, velocities and accelerations are zero in the
% beginning and at time T, we add fifth order polynomial to fourier
% series. The parameters of the polynomial depends on the parameters of
% fourier series. Here we compute them.
c_pol = getPolCoeffs(T, a, b, wf, N, traj_par.q0);

% Compute trajectory (Fouruer series + fifth order polynomail)
[q,qd,q2d] = mixed_traj(t, c_pol, a, b, wf, N);

% Obtain observation matrix by computing regressor for each sampling time
E1 = baseQR.permutationMatrix(:,1:baseQR.numberOfBaseParameters);
W = [];
for i = 1:length(t)
%     Y = base_regressor_UR10E(q(:,i),qd(:,i),q2d(:,i));
    if baseQR.motorDynamicsIncluded
        Y = [regressorWithMotorDynamics(q(:,i),qd(:,i),q2d(:,i))*E1, ...
             frictionRegressor(qd(:,i))];
    else
        Y = [full_regressor_UR10E(q(:,i),qd(:,i),q2d(:,i))*E1, ...
             frictionRegressor(qd(:,i))];
    end
    W = vertcat(W,Y);
end

out = cond(W);

author: shamilmamedov


text:
function spline_test07 ( )

%*****************************************************************************80
%
%% TEST07 tests BASIS_MATRIX_OVERHAUSER_UNI, BASIS_MATRIX_TMP.
%
%  Discussion:
%
%   YDATA(1:NDATA) = ( TDATA(1:NDATA) + 2 )**2 + 3
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author
%
%    John Burkardt
%
  n = 4;
  ndata = 4;
  nsample = 4;
  tdata = [ -1.0E+00, 0.0E+00, 1.0E+00, 2.0E+00 ];
  ydata = [ 4.0E+00, 7.0E+00, 12.0E+00, 19.0E+00 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST07\n' );
  fprintf ( 1, '  BASIS_MATRIX_OVERHAUSER_UNI sets up the basis\n' );
  fprintf ( 1, '    matrix for the uniform Overhauser spline.\n' );

  mbasis = basis_matrix_overhauser_uni ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '    TDATA, YDATA\n' );
  fprintf ( 1, '\n' );
  for i = 1 : ndata
    fprintf ( 1, '%14f  %14f\n', tdata(i), ydata(i) );
  end

  left = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, '    T, Spline(T)\n' );
  fprintf ( 1, '\n' );

  for i = 0 : ndata

    if ( i == 0 )
      tlo = tdata(1) - 0.5E+00 * ( tdata(2) - tdata(1) );
      thi = tdata(1);
    elseif ( i < ndata )
      tlo = tdata(i);
      thi = tdata(i+1);
    elseif ( ndata <= i )
      tlo = tdata(ndata);
      thi = tdata(ndata) + 0.5E+00 * ( tdata(ndata) - tdata(ndata-1) );
    end

    if ( i < ndata )
      jhi = nsample - 1;
    else
      jhi = nsample;
    end

    for j = 0 : jhi

      tval = ( ( nsample - j ) * tlo   ...
             + (           j ) * thi ) ...
             /   nsample;

      yval = basis_matrix_tmp ( left, n, mbasis, ndata, tdata, ydata, tval );

      if ( 0 < i & j == 0 )
        mark = '*';
      else
        mark = ' ';
      end

      fprintf ( 1, '%c  %14f  %14f\n', mark, tval, yval );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [C, H, W, M] = FindLargestRectangles(I, crit, minSize)
% finds largest rectangle regions within all points set to 1.
% input: I       - B/W boolean matrix or output of FindLargestSquares
%        minSize - [height width] - minimum width and height on regions of
%                  interest (used to restrict final choise)
%        crit    - Optimazation Criteria parameters to optimize:
%                   crit(1)*height + crit(2)*width + crit(3)*height*width
% output:
%         C    - value of the optimization criteria "crit" calculated for
%                each pixel
%         W, H - for each pixel I(r,c) return height and width of the largest
%                all-white rectangle with its upper-left corner at I(r,c)
%         M    - Mask the largest all-white rectangle of the image

if (nargin<2)
  crit = [1 1 0];
end
if (nargin<3)
  minSize = [1 1];
end
p = crit;
[nR nC] = size(I);
if (minSize(1)<1), minSize(1)= floor(minSize(1)*nR); end
if (minSize(2)<1), minSize(2)= floor(minSize(2)*nC); end
if (max(I(:)) - min(I(:))==1),
  S = FindLargestSquares(I);
else
  S = I;
end
n = max(S(:));
W = S; % make a carbon copy of the matrix data
H = S;
C = ((p(1)+p(2)) + p(3)*S) .* S; % p(1)*width + p(2)*height + p(3)*height*width for height=width=S;
d = round((3*n)/4);
minH = max(min(minSize(1), d),1);
minW = max(min(minSize(2), d),1);

%% look for rectangles with width>height
hight2width = zeros(n+1,1);  % Store array with largest widths aviable for a given height
for r = 1 : nR               % each row is processed independently
  hight2width(:) = 0;        % reset the List
  for c = nC: -1 : 1         % go through all pixels in a row right to left
    s = S(r,c);              % s is a size of a square with its corner at (r,c)
    if (s>0)                 % if pixel I(r,c) is true
      MaxCrit = C(r,c);      % initialize the Max Criteria using square
      for hight = s:-1:1     % go through all possible width&hight combinations. Start with more likely to be the best
        width = hight2width(hight); % look up width for a given hight
        width = max(width+1,s);
        hight2width(hight) = width;
        Crit = p(1)*hight + p(2)*width + p(3)*width*hight;
        if (Crit>MaxCrit),   % check if it produces larger Criteria
          MaxCrit = Crit;    % if it does than save the results
          W(r,c)  = width;
          H(r,c)  = hight;
        end % if Crit
      end % for hight
      C(r,c)  = MaxCrit;
    end % if s
    hight2width((s+1):end) = 0;    % hights>s will not be aviable for the next pixel
  end % for c
end
clear hight2width

%% look for rectangles with width<height
width2hight = zeros(n+1,1);  % Store array with largest widths aviable for a given height
for c = 1 : nC               % each column is processed independently
  width2hight(:) = 0;        % reset the List
  for r = nR: -1 : 1         % go through all pixels in a column bottom to top
    s = S(r,c);              % s is a size of a square with its corner at (r,c)
    if (s>0)                 % if pixel I(r,c) is true
      MaxCrit = C(r,c);      % initialize the Max Criteria using square
      for width = s:-1:1     % go through all possible width&hight combinations. Start with more likely to be the best
        hight = width2hight(width); % look up hight for a given width
        hight = max(hight+1,s);
        width2hight(width) = hight;
        Crit = p(1)*hight + p(2)*width + p(3)*width*hight;
        if (Crit>MaxCrit),   % check if it produces larger Criteria
          MaxCrit = Crit;    % if it does than save the results
          W(r,c)  = width;
          H(r,c)  = hight;
        end % if Crit
      end % for width
      C(r,c)  = MaxCrit;
    end % if s
    width2hight((s+1):end) = 0;    % hights>s will not be aviable for the next pixel
  end % for r
end

%% Create container mask
CC = C;
CC( H<minH | W<minW ) = 0; % first try to obey size restrictions
[~, pos] = max(CC(:));
if (isempty(pos)), [~, pos] = max(C(:)); end % but when it fails than drop them
[r c] = ind2sub(size(C), pos);
M = false(size(C));
M( r:(r+H(r,c)-1), c:(c+W(r,c)-1) ) = 1;

author: Sable

repo: mcbench-benchmarks


text:
function a = sqrt(a)
%SQRT         Gradient square root sqrt(a)
%

% written  10/16/98     S.M. Rump
% modified 10/14/00     S.M. Rump  use Tony's trick
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/08/08     S.M. Rump  improved sparse multiplication: not using intval data type
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  N = getappdata(0,'INTLAB_GRADIENT_NUMVAR');

  wng = warning;
  warning off

  a.x = sqrt(a.x);
  % use full(a.x(:)): cures Matlab V6.0 bug
  % a=7; i=[1 1]; x=a(i), b=sparse(a); y=b(i)  yields row vector x but column vector y
  % ax is full anyway
  ax = 1./(2*full(a.x(:)));
  if issparse(a.dx)
    sizeax = size(a.dx,1);
    [ia,ja,sa] = find(a.dx);
    if isa(a.x,'intval')
      adx = times(ax(ia),sa,0);
      if adx.complex
        a.dx = intval( sparse(ia,ja,adx.mid,sizeax,N) , sparse(ia,ja,adx.rad,sizeax,N) , 'midrad' );
      else
        a.dx = intval( sparse(ia,ja,adx.inf,sizeax,N) , sparse(ia,ja,adx.sup,sizeax,N) , 'infsup' );
      end
    else
      a.dx = sparse(ia,ja,ax(ia).*sa,sizeax,N);
    end
  else
    a.dx = a.dx .* ax(:,ones(1,N));
  end

  warning(wng)

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function [ DCM ] = tapas_ceode_synData_erp()
% [ DCM ] = tapas_ceode_synData_erp()
%
% Sets parameters and function handles to create a synthetic DCM
% structure. The resulting structure can then be used to generate
% artificial data with tapas_ceode_gen_erp.m
%
% INPUT
%
% OUTPUT
%   DCM         struct          DCM structure with all necessary fields to
%                               generate synthetic data.
%
% -------------------------------------------------------------------------
%
% Author: Dario Schöbi
% Created: 2020-08-10
% Copyright (C) 2020 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS ceode Toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

% ---------------------- MODEL SPECIFICATION ------------------------------
M = struct();

% Handles to specify integration scheme, dynamical equation, forward model
% and input
M.IS = 'tapas_ceode_gen_erp';
M.f = 'spm_fx_cmc';
M.int = 'spm_int_L';
M.G = 'spm_lx_erp';
M.FS = 'spm_fy_erp';
M.fu = 'spm_erp_u';

% Design/Parameter values of the driving input (see spm_erp_u.m)
M.dur = 16;
M.ons = 64;
M.ns = 500;
M.x = zeros(2, 9);

% ---------------------- DESIGN SPECIFICATION -----------------------------
xU = struct();

% No condition specific effec
xU.X = [0]';
xU.name = {'std'};

% Sampling rate of the artificial signal
xU.dt = 1e-03;

% ---------------------- PARAMETER VALUES ---------------------------------
% Neuronal parameters (also see spm_fx_erp.m / tapas_ceode_fx_erp.m)
Ep = struct();

Ep.T = [0, 0; 0 0];
Ep.G = [0, 0; 0 0];
Ep.S = [0, 0];
Ep.A{1} = [-4, -4; 1, -4];
Ep.A{2} = [-4, 0; -4, -4];
Ep.A{3} = -4 * ones(2, 2);
Ep.B{1} = zeros(2, 2);
Ep.C = [0; -4];
Ep.H = [0 0 0 0];
Ep.D = -10 * ones(2, 2);
Ep.R = [0 0];

% Leadfield parameters (also see spm_lx_erp.m)
Eg.L = [1, 2];
Eg.J = [1, zeros(1, 5), 1, 0, 1];

% ---------------------- CREATE OUTPUT ------------------------------------
% Combine all structures into single DCM structure
DCM.Ep = Ep;
DCM.Eg = Eg;
DCM.M = M;
DCM.xU = xU;

end

author: translationalneuromodeling

repo: tapas


text:
function [feature1, feature2, feature3] = SFAfeature(im_lists, layer_name, options)
% SFA features

% written by Dingquan Li
% dingquanli@pku.edu.cn
% IDM, SMS, PKU
% Last update: Aug. 9, 2017

if ~exist('options','var')
    options = 'none';
end

% multi-patch representation
patchSize = [224 224]; %
stride = floor(patchSize/2); %

% Initialize the DCNN
model = [fileparts(which(mfilename)) '/models/ResNet-50-deploy.prototxt']; %
weights = [fileparts(which(mfilename)) '/models/ResNet-50-model.caffemodel']; %
caffe.set_mode_gpu(); %
net = caffe.Net(model, weights, 'test');

l = net.blobs(layer_name).shape;
l = l(3);
feature1 = zeros(length(im_lists),2*l);
feature2 = zeros(length(im_lists),5*l);
feature3 = zeros(length(im_lists),4*l);
for k = 1:length(im_lists)
    fprintf(['Extracting the feature of the ' num2str(k) 'th image ...\n']);

    im = imread(im_lists{k});

    patches = im2patches(im, patchSize, stride);

    % feature extraction
    patch_features = DCNN(patches, net, layer_name);

    % feature aggregation
    feature_mean = mean(patch_features,2);
    feature_std = std(patch_features,0,2);

    feature_q4 = prctile(patch_features,100,2);
    feature_q3 = prctile(patch_features,75,2);
    feature_q2 = prctile(patch_features,50,2);
    feature_q1 = prctile(patch_features,25,2);
    feature_q0 = prctile(patch_features,0,2);

    feature_Moment2 = moment(patch_features,2,2);
    feature_Moment3 = moment(patch_features,3,2);
    feature_Moment4 = moment(patch_features,4,2);

    feature1(k,:) = [feature_mean;feature_std]'; % mean&std aggregation
    feature2(k,:) = [feature_q4;feature_q3;feature_q2;feature_q1;feature_q0]'; % quantile aggregation
    feature3(k,:) = [feature_mean; nthroot(feature_Moment2,2);nthroot(feature_Moment3,3);...
        nthroot(feature_Moment4,4)]'; % moment aggregation

    switch options
        case 'flipH'
            patches = im2patches(im(:,end:-1:1,:), patchSize, stride);

            % feature extraction
            patch_features = DCNN(patches, net, layer_name);

            % feature aggregation
            feature_mean = mean(patch_features,2);
            feature_std = std(patch_features,0,2);

            feature_q4 = prctile(patch_features,100,2);
            feature_q3 = prctile(patch_features,75,2);
            feature_q2 = prctile(patch_features,50,2);
            feature_q1 = prctile(patch_features,25,2);
            feature_q0 = prctile(patch_features,0,2);

            feature_Moment2 = moment(patch_features,2,2);
            feature_Moment3 = moment(patch_features,3,2);
            feature_Moment4 = moment(patch_features,4,2);

            feature1(k+length(im_lists),:) = [feature_mean;feature_std]'; % mean&std aggregation
            feature2(k+length(im_lists),:) = [feature_q4;feature_q3;feature_q2;feature_q1;feature_q0]'; % quantile aggregation
            feature3(k+length(im_lists),:) = [feature_mean; nthroot(feature_Moment2,2);nthroot(feature_Moment3,3);...
                nthroot(feature_Moment4,4)]'; % moment aggregation
        case 'clipUL'
            T = 6;
            for i = 0:T-1
                for j = 0:T-1
                    if i^2+j^2==0
                        continue
                    end
                    patches = im2patches(im(i+1:end,j+1:end,:), patchSize, stride);

                    % feature extraction
                    patch_features = DCNN(patches, net, layer_name);

                    % feature aggregation'res5c'
                    feature_mean = mean(patch_features,2);
                    feature_std = std(patch_features,0,2);

                    feature_q4 = prctile(patch_features,100,2);
                    feature_q3 = prctile(patch_features,75,2);
                    feature_q2 = prctile(patch_features,50,2);
                    feature_q1 = prctile(patch_features,25,2);
                    feature_q0 = prctile(patch_features,0,2);

                    feature_Moment2 = moment(patch_features,2,2);
                    feature_Moment3 = moment(patch_features,3,2);
                    feature_Moment4 = moment(patch_features,4,2);

                    feature1(k+(i*T+j)*length(im_lists),:) = [feature_mean;feature_std]'; % mean&std aggregation
                    feature2(k+(i*T+j)*length(im_lists),:) = [feature_q4;feature_q3;feature_q2;feature_q1;feature_q0]'; % quantile aggregation
                    feature3(k+(i*T+j)*length(im_lists),:) = [feature_mean; nthroot(feature_Moment2,2);nthroot(feature_Moment3,3);...
                        nthroot(feature_Moment4,4)]'; % moment aggregation
                end
            end
    end
end
caffe.reset_all();

function patch_features = DCNN(patches, net, layer_name)
% Extract high-level semantic features from pool5 of ResNet-50

% Preprocess the patches
im_data = patches(:,:,[3 2 1],:); % RGB2BGR
im_data = permute(im_data,[2 1 3 4]); % HWC2WHC
im_data = single(im_data); % single

cropped_dim = size(patches,1);
mean_data = caffe.io.read_mean([fileparts(which(mfilename)) '/models/imagenet_mean.binaryproto']); %
topleft = floor((size(mean_data,1)-cropped_dim)/2)+1;
mean_data = mean_data(topleft:topleft+cropped_dim-1,topleft:topleft+cropped_dim-1,:);

crops_data = arrayfun(@(i)im_data(:,:,:,i)-mean_data,1:size(patches,4),'UniformOutput',false);

% Extract features
l = net.blobs(layer_name).shape; % 2048
l = l(3);
patch_features = zeros(l,length(crops_data));
for n = 1:length(crops_data)
    net.forward(crops_data(n));
%     patch_features(:,n) = std(reshape(net.blobs('res5c').get_data(),[],l,1),0,1); % pool5
    patch_features(:,n) = mean(reshape(net.blobs(layer_name).get_data(),[],l,1)); % pool5
%     patch_features(:,n) = max(reshape(net.blobs('res5c').get_data(),[],l,1)); % pool5
end
patch_features = double(patch_features);

function [patches,imcol,I] = im2patches(I, patchSize, stride)
% IM2PATCHES Extract rectangular patches from input image.
%
%   patches = IM2PATCHES(I,patchSize) I can be a MxN matrix (e.g. gray-
%   scale image) or a MxNxK array (e.g. RGB image). patchSize can be:
%   a) a scalar, in which case: patchHeight = patchWidth = patchSize.
%   b) a vector [patchHeight,patchWidth]
%   c) a Mx4 matrix with bounding box coordinates in the form
%   [xmin,ymin,xmax,ymax].
%   patches is a patchHeight x patchWidth x K x nPatches array (K >= 1).
%
%   patches = IM2PATCHES(I,patchSize,stride) stride between consecutive
%   patches. Stride can be used to extract overlapping patchesand can be
%   either a scalar, or a 2x1 vector to define different strides across
%   axes x and y.
%
%   [patches,imcol,I] = IM2PATCHES(...) also returns imcol, which is a
%   nPixelsPerPatch x nPatches matrix whose columns are the elements of
%   each patch, and I, which is the input image I after padding with zeros.
%
%   USAGE EXAMPLES:
%   patches = im2patches(I,[h,w]);   % equal to: im2col(I,[h,w],'distinct')
%   patches = im2patches(I,[h,w],1); % equal to: im2col(I,[h,w],'sliding')
%
%   NOTE: IM2PATCHES extracts patches padding with zeros when necessary.
%   The only exception is when stride == 1, or when the stride has a value
%   that the last patches horizontally and vertically fit precicely in the
%   image, without crossing the borders.
%
% See also: patches2im, im2col, col2im
%
% Stavros Tsogkas, <stavros.tsogkas@ecp.fr>
% Last update: August 2015

assert(ismatrix(I) || ndims(I) == 3, 'Input must be a 2D or 3D array');
if ismatrix(patchSize) && size(patchSize,2) == 4
    warning('This part of the function has not been tested')
    % patchSize is a Mx4 matrix containing the [xmin,ymin,xmax,ymax]
    % coordinates of bounding boxes that are fully contained in the image
    bb      = patchSize;
    bb(:,1) = max(1,bb(:,1));
    bb(:,2) = max(1,bb(:,2));
    bb(:,3) = min(size(I,2),bb(:,3));
    bb(:,4) = min(size(I,1),bb(:,4));
    patches = cell(size(bb,1),1);
    for i=1:size(bb,1)
        patches{i} = I(bb(i,2):bb(i,4),bb(i,1):bb(i,3),:);
    end
else
    assert(all(patchSize > 0), 'Patch size cannot be negative or zero')
    if isscalar(patchSize)          % Square patches
        patchHeight = patchSize;
        patchWidth  = patchSize;
    elseif numel(patchSize) == 2    % Rectangular patches
        patchHeight = patchSize(1);
        patchWidth  = patchSize(2);
    else
        error('patchSize can be either a scalar or a [patchHeight, patchWidth] vector.')
    end
    if nargin < 3                   % Identical strides for X and Y axis.
        strideX = patchWidth;       % Default is equivalent to 'distinct'
        strideY = patchHeight;      % option for Matlab's im2col.
    elseif isscalar(stride)
        strideY = stride;
        strideX = stride;
    elseif numel(stride) == 2       % Different strides for X and Y axis.
        strideY = stride(1);
        strideX = stride(2);
    else
        error('Stride can be either a scalar or a [strideX, strideY] vector.')
    end
    [hin,win,din] = size(I);
    hout = hin - mod(hin-patchHeight,strideY) + strideY*(strideY > 1);
    wout = win - mod(win-patchWidth, strideX) + strideX*(strideX > 1);
    nPixelsPerPatch = patchWidth*patchHeight*din;
    I(end+1:hout,end+1:wout,:) = 0; % pad with zeros if necessary
    % x-y indices for a single (possibly N-dimensional) patch
    [x,y,z] = meshgrid(1:patchWidth,1:patchHeight,1:din);
    % pixel indices for all patches
    [xstart,ystart] = meshgrid(0:strideX:(wout-patchWidth),0:strideY:(hout-patchHeight));
    inds    = bsxfun(@plus, reshape(y,nPixelsPerPatch,[]), ystart(:)');
    inds    = inds + (bsxfun(@plus, hout*reshape(x-1,nPixelsPerPatch,[]), hout*xstart(:)'));
    inds    = bsxfun(@plus, inds, (z(:)-1)*(hout*wout));
    imcol   = I(inds);
    patches = reshape(imcol,patchHeight,patchWidth,din,[]);
end
% =========================================================================
% Alternative versions for Matlab's im2col 'distinct' and 'sliding' modes.
% These versions are a little faster than the code used in im2patches for
% creating the imcol image, but they cannot handle an arbitrary stride
% between neighboring patches.
% =========================================================================
% % -------------------------------------------------------------------------
% function out = im2col_distinct(A,blocksize)
% % -----------------------------------------------------------------------
% nrows = blocksize(1);
% ncols = blocksize(2);
% nele = nrows*ncols;
%
% row_ext = mod(size(A,1),nrows);
% col_ext = mod(size(A,2),ncols);
%
% padrowlen = (row_ext~=0)*(nrows - row_ext);
% padcollen = (col_ext~=0)*(ncols - col_ext);
%
% A1 = zeros(size(A,1)+padrowlen,size(A,2)+padcollen);
% A1(1:size(A,1),1:size(A,2)) = A;
%
% t1 = reshape(A1,nrows,size(A1,1)/nrows,[]);
% t2 = reshape(permute(t1,[1 3 2]),size(t1,1)*size(t1,3),[]);
% t3 =  permute(reshape(t2,nele,size(t2,1)/nele,[]),[1 3 2]);
% out = reshape(t3,nele,[]);
%
% return;
%
% % -------------------------------------------------------------------------
% function out = im2col_sliding(A,blocksize)
% % -------------------------------------------------------------------------
%
% nrows = blocksize(1);
% ncols = blocksize(2);
%
% %// Get sizes for later usages
% [m,n] = size(A);
%
% %// Start indices for each block
% start_ind = reshape(bsxfun(@plus,[1:m-nrows+1]',[0:n-ncols]*m),[],1); %//'
%
% %// Row indices
% lin_row = permute(bsxfun(@plus,start_ind,[0:nrows-1])',[1 3 2]);  %//'
%
% %// Get linear indices based on row and col indices and get desired output
% out = A(reshape(bsxfun(@plus,lin_row,[0:ncols-1]*m),nrows*ncols,[]));
%
% return;

repo: SFA


text:
%% chkarg
% Throw an exception if the argument is |false|.

%%% Syntax
%  chkarg(truth, errmsg, varargin)

%%% Description
% |chkarg(truth, errmsg, varargin)| throws an exception of |truth| is |false|.
% |errmsg| is an error message embedded in the exception.  The error message
% can be generated by formatting additional arguments |varargin|; in this case,
% |errmsg| should be a format string.

%%% Example
%   chkarg(isreal(3+4i), '3+4i is not real.');

function chkarg(truth, errmsg, varargin)

errmsg = sprintf(errmsg, varargin{:});

if ~truth
	exception = MException('Maxwell:argChk', 'Wrong argument: %s', errmsg);
	throwAsCaller(exception);
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function test_ft_networkanalysis

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_networkanalysis

data = [];
for k = 1:5
  x = randn(7,15);
  cx = x*x';
  cx = cx./sqrt(diag(x)*diag(x)');
  data.cohspctrm(:,:,k) = cx;
end
data.freq  = 1:5;
data.label = {'chan1';'chan2';'chan3';'chan4';'chan5';'chan6';'chan7'};
data.dimord = 'chan_chan_freq';
data.cfg    = 'this is the cfg';

tmp = data;
tmp.cohspctrm = data.cohspctrm>0.3;

% at present just checks for undirected binary and weighted graphs
cfg           = [];
cfg.parameter = 'cohspctrm';

cfg.method    = 'assortativity';
stat1 = ft_networkanalysis(cfg, tmp);
stat2 = ft_networkanalysis(cfg, data);

cfg.method    = 'betweenness';
stat3 = ft_networkanalysis(cfg, tmp);
stat4 = ft_networkanalysis(cfg, data);

cfg.method    = 'clustering_coef';
stat5 = ft_networkanalysis(cfg, tmp);
stat6 = ft_networkanalysis(cfg, data);

cfg.method    = 'degrees';
stat7 = ft_networkanalysis(cfg, tmp);
stat8 = ft_networkanalysis(cfg, data);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = dual( varargin )

%DUAL Declares one or more dual variables.
%
%   DUAL VARIABLE x
%   where x is a valid MATLAB variable name, declares a dual variable
%   for the current cvx problem. A dual variable with that name is added
%   to the problem and a cvxdual object with that name is created in
%   the current worksapce. An error is generated if a cvx problem is not
%   in the current workspace.
%
%   Note that a dual variable initially has no size; i.e., SIZE(x) returns
%   [0,0]. This is because the size of a dual variable is determined by the
%   constraint to which it is attached. To attach a dual variable to a
%   constraint, use the colon notation as follows:
%
%      variable x(n)
%      dual variable z
%      A * x == b : z
%
%   In this example, the dual variable z is attached to the equality
%   constraint A * x == b. Assuming SIZE(A) = [m,n] and SIZE(B) = [m,1],
%   the size of z will be set at [m,1] by this operation.
%
%   DUAL VARIABLES x1 x2 ... xk
%   where x1, x2, ..., xk are valid MATLAB variable names, declares
%   multiple dual variables.
%
%   DUAL VARIABLE x{s1,s2,...,sk}
%   DUAL VARIABLES x1{s1,s2,...,sk} x2{...
%   creates a cell array of dual variables, of dimension [s1,s2,..,sk].
%   This is useful if you want to create a sequence of related constraints,
%   using a FOR loop for example, and assign a separate dual variable to
%   each one; for example:
%
%       variable x(n)
%       dual variables z{10}
%       for k = 1 : 10,
%           A(:,:,k) * x <= b(:,k) : z{k}
%       end
%
%   Assuming SIZE(A) = [m,n,10] and SIZE(B) = [m,10], this loop assigns
%   each of the dual variables z{k} to a constraint and sets their sizes
%   to [m,1]. Note that the elements of the dual variable cell array need
%   not end up the same size.
%
%   See also VARIABLE, VARIABLES.

if nargin < 2,
    error( 'Incorrect syntax for DUAL VARIABLE(S). Type HELP DUAL for details.' );
elseif nargout && nargout ~= length( varargin ) - 1,
    error( 'Incorrect number of output arguments.' );
elseif ~iscellstr( varargin ),
    error( 'All arguments must be strings.' );
end

global cvx___
prob = evalin( 'caller', 'cvx_problem', '[]' );
if ~isa( prob, 'cvxprob' ),
    error( 'No CVX model exists in this scope.' );
elseif isempty( cvx___.problems ) || cvx___.problems( end ).self ~= prob,
    error( 'Internal CVX data corruption. Please CLEAR ALL and rebuild your model.' );
end

if strcmp( varargin{1}, 'variable' ),
    if nargin > 2,
        error( 'Too many input arguments.\nTrying to declare multiple dual variables? Use the DUAL VARIABLES command instead.', 1 ); %#ok
    end
elseif ~strcmp( varargin{1}, 'variables' ),
    error( 'Incorrect syntax for DUAL VARIABLE(S). Type HELP DUAL for details.' );
end

for k = 2 : nargin,
    arg = varargin{k};
    toks = regexp( arg, '^\s*([a-zA-Z]\w*)\s*({.*})?\s*$', 'tokens' );
    if isempty( toks ),
        error( 'Invalid dual variable specification: %s', arg );
    end
    tok = toks{1};
    nam = tok{1};
    if length(tok) < 2,
        siz = {};
    else
        siz = tok{2};
    end
    if nam(end) == '_',
        error( 'Invalid dual variable specification: %s\n   Variables ending in underscores are reserved for internal use.', arg );
    end
    if ~isempty( siz ),
        try
            siz = evalin( 'caller', [ '[', siz(2:end-1), ']' ] );
        catch exc
            error( exc.identifier, 'Error attempting to determine size of: %s\n   %s', arg, exc.message );
        end
        [ temp, siz ] = cvx_check_dimlist( siz, true );
        if ~temp,
            error( 'Invalid dual variable specification: %s\n   Dimension list must be a vector of finite nonnegative integers.', arg );
        end
    end
    temp = newdual( prob, nam, siz );
    if nargout,
        varargout{k-1} = temp; %#ok
    end
    assignin( 'caller', nam, temp );
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function demoNapoleon(varargin)
%DEMONAPOLEON  Small demo of a theorem on triangles
%
%   output = demoNapoleon(input)
%
%   Example
%   demoNapoleon
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2013-08-19,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2013 INRA - Cepia Software Platform.

%% base triangle

% choose 3 points on the plane
A = [10 8];
B = [3 2];
C = [12 3];

% create a polygon for the base triangle
ABC = [A ; B ; C];

% Draw the base triangle
figure;
drawPolygon(ABC, 'color', 'k', 'lineWidth', 2);
axis equal;
axis([-5 20 -8 12]);
hold on;

%% Equliateral triangles

% create an equilateral triangle from each side, pointing outwards
tAB = createTriangle(B, A);
tBC = createTriangle(C, B);
tAC = createTriangle(A, C);

% draw each equilateral triangle
drawPolygon(tAB, 'color', 'b');
drawPolygon(tBC, 'color', 'b');
drawPolygon(tAC, 'color', 'b');

%% Compute centroid of each triangle

% compute centroids
c1 = polygonCentroid(tAB);
c2 = polygonCentroid(tBC);
c3 = polygonCentroid(tAC);

% draw the triangle formed by the centroids
triC = [c1; c2; c3];
drawPolygon(triC, 'lineWidth', 2, 'color', [0 .8 0]);
drawPoint(triC, 'marker', 'o', 'markerSize', 8, 'markerfacecolor', [0 .8 0]);

function tri = createTriangle(p1, p2)

d12 = distancePoints(p1, p2);
c1 = [p1 d12];
c2 = [p2 d12];
inters = intersectCircles(c1, c2);

indPos = isCounterClockwise(p1, p2, inters) > 0;
p3 = inters(indPos, :);

tri = [p1; p2; p3];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% OP_GRADGRADU_GRADGRADV: assemble the matrix A = [a(i,j)], a(i,j) = (epsilon grad grad u_j, grad grad v_i).
%
%   mat = op_gradgradu_gradgradv (spu, spv, msh, epsilon);
%   [rows, cols, values] = op_gradgradu_gradgradv (spu, spv, msh, epsilon);
%
% INPUT:
%
%   spu:   structure representing the space of trial functions (see sp_scalar/sp_evaluate_col)
%   spv:   structure representing the space of test functions (see sp_scalar/sp_evaluate_col)
%   msh:   structure containing the domain partition and the quadrature rule (see msh_cartesian/msh_evaluate_col)
%   epsilon: diffusion coefficient
%
% OUTPUT:
%
%   mat:    assembled matrix
%   rows:   row indices of the nonzero entries
%   cols:   column indices of the nonzero entries
%   values: values of the nonzero entries
%
% Copyright (C) 2009, 2010 Carlo de Falco
% Copyright (C) 2011, 2017 Rafael Vazquez
% Copyright (C) 2013, Marco Pingaro
%
%    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 3 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, see <http://www.gnu.org/licenses/>.

function varargout = op_gradgradu_gradgradv (spu, spv, msh, coeff)

  der2u = reshape (spu.shape_function_hessians, spu.ncomp, [], msh.nqn, spu.nsh_max, msh.nel);
  der2v = reshape (spv.shape_function_hessians, spv.ncomp, [], msh.nqn, spv.nsh_max, msh.nel);

  ndir = size (der2u, 2);

  rows = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  cols = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  values = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);

  jacdet_weights = msh.jacdet .* msh.quad_weights .* coeff;

  ncounter = 0;
  for iel = 1:msh.nel
    if (all (msh.jacdet(:, iel)))
      der2u_iel = reshape (der2u(:,:,:,:,iel), spu.ncomp*ndir, msh.nqn, 1, spu.nsh_max);
      der2v_iel = reshape (der2v(:,:,:,:,iel), spv.ncomp*ndir, msh.nqn, spv.nsh_max, 1);

      jacdet_iel = reshape (jacdet_weights(:,iel), [1,msh.nqn,1,1]);

      jacdet_der2u = bsxfun (@times, jacdet_iel, der2u_iel);
      tmp1 = sum (bsxfun (@times, jacdet_der2u, der2v_iel), 1);
      elementary_values = reshape (sum (tmp1, 2), spv.nsh_max, spu.nsh_max);

      [rows_loc, cols_loc] = ndgrid (spv.connectivity(:,iel), spu.connectivity(:,iel));
      indices = rows_loc & cols_loc;
      rows(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = rows_loc(indices);
      cols(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = cols_loc(indices);
      values(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = elementary_values(indices);
      ncounter = ncounter + spu.nsh(iel)*spv.nsh(iel);

    else
      warning ('geopdes:jacdet_zero_at_quad_node', 'op_gradgradu_gradgradv: singular map in element number %d', iel)
    end
  end

  if (nargout == 1 || nargout == 0)
    varargout{1} = sparse (rows(1:ncounter), cols(1:ncounter), ...
                           values(1:ncounter), spv.ndof, spu.ndof);
  elseif (nargout == 3)
    varargout{1} = rows(1:ncounter);
    varargout{2} = cols(1:ncounter);
    varargout{3} = values(1:ncounter);
  else
    error ('op_gradgradu_gradgradv: wrong number of output arguments')
  end

end

%% COPY OF THE FIRST VERSION OF THE FUNCTION (MORE UNDERSTANDABLE)
% function mat = op_gradgradu_gradgradv (spu, spv, msh, coeff)
%   mat = spalloc (spv.ndof, spu.ndof, 1);
%
%   der2u = reshape (spu.shape_function_hessians, spu.ncomp, [], msh.nqn, spu.nsh_max, msh.nel);
%   der2v = reshape (spv.shape_function_hessians, spv.ncomp, [], msh.nqn, spv.nsh_max, msh.nel);
%
%   ndir = size (der2u, 2);
%
%   for iel = 1:msh.nel
%     if (all (msh.jacdet(:,iel)))
%       mat_loc = zeros (spv.nsh(iel), spu.nsh(iel));
%       for idof = 1:spv.nsh(iel)
%           ishh = reshape(der2v(:,:,:,idof,iel), spv.ncomp * ndir, []);
%           for jdof = 1:spu.nsh(iel)
%               jshh = reshape(der2u(:,:,:,jdof,iel), spu.ncomp * ndir, []);
%           % The cycle on the quadrature points is vectorized
%           %for inode = 1:msh.nqn
%               mat_loc(idof, jdof) = mat_loc(idof, jdof) + ...
%                   sum (msh.jacdet(:,iel) .* msh.quad_weights(:, iel) .* ...
%                   sum (ishh .* jshh, 1).' .* coeff(:,iel));
%           %end
%           end
%       end
%       mat(spv.connectivity(:, iel), spu.connectivity(:, iel)) = ...
%           mat(spv.connectivity(:, iel), spu.connectivity(:, iel)) + mat_loc;
%     else
%       warning ('geopdes:jacdet_zero_at_quad_node',...
%           'op_gradgradu_gradgradv: singular map in element number %d', iel)
%     end
%   end
% end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
%% RATE OF CONVERGENCE OF PRIMAL-DUAL WEAK-GALERKIN for DIV-CURL SYSTEM
%
% This example is to show the rate of convergence on an L-shaped domain with a singular solution.
%
% See also DivCurl3PDWG
%
% Reference: A New Numerical Method for Div-Curl Systems with Low Regularity Assumptions
% S. Cao, C. Wang, J. Wang, https://arxiv.org/abs/2101.03466
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

clear; close all;
%%
maxIt = 4;
N = zeros(maxIt,1);
h = zeros(maxIt,1);
errL2UTotal = zeros(maxIt,1);
errqlambdaTotal = zeros(maxIt,1);
errSTotal = zeros(maxIt,1);
errL2QuTotal = zeros(maxIt,1); % \|Q_h u - u \|

%% Generate an initial mesh
[node,elem] = cubemesh([-1,1,-1,1,-1,1],1/2);
% [node,elem] = delmesh(node,elem,'(x < 0 & y > 0) | (z<0.5)');

% [node,elem] = delmesh(node,elem,'(x < 0 & y < 0) | (z<0) | (z>0.5)'); %
% for gradient

[node,elem] = delmesh(node,elem,'(x > 0 & y < 0) | (z<0) | (z>0.5)');
% for curl
bdFlag = setboundary3(node,elem,'Dirichlet');

%% data
pde = DataDivCurlSingular4;
option.inverseh = true;

%% Finite Element Method
for k = 1:maxIt
    % solve the equation
    [soln,eqn,info] = DivCurl3PDWG(node,elem,bdFlag,pde,option);
    N(k) = size(elem,1);
%     h(k) = (36./N(k)).^(1/3);
    h(k) = 1/2^k;

    errL2UTotal(k) = eqn.errorU;
    errqlambdaTotal(k) = eqn.errorqlambda;
    errSTotal(k) = eqn.errorS;

    [errL2QuTotal(k), ~] = getL2error3vec(node,elem,soln.u,pde);
%     [errL2QuTotal(k), ~] = getL2error3(node,elem,pde.exactu, soln.u);

    fprintf('PDWG method with h:     %3g \n', h(k));
    fprintf('Number of Dof:          %d \n', sum(eqn.freeDof));
    fprintf('L2 error in u:          %6g \n', errL2UTotal(k));
    fprintf('L2 error in Qu-u_h:     %6g \n', errL2QuTotal(k));
    fprintf('L2 error in (q,lambda): %6g \n', errqlambdaTotal(k));
    fprintf('L2 error in s:          %6g \n\n\n', errSTotal(k));

    % refine mesh
   if k < maxIt; [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag); end
end

%%
visualizationPDwgDivCurl;

repo: ifem


text:
function [data] = ft_megrealign(cfg, data)

% FT_MEGREALIGN interpolates MEG data towards standard gradiometer locations by
% projecting the individual timelocked data towards a coarse source reconstructed
% representation and computing the magnetic field on the standard gradiometer
% locations.
%
% Use as
%   [interp] = ft_megrealign(cfg, data)
%
% Required configuration options:
%   cfg.template
%   cfg.inwardshift
%
% The new gradiometer definition is obtained from a template dataset,
% or can be constructed by averaging the gradiometer positions over
% multiple datasets.
%   cfg.template       = single dataset that serves as template
%   cfg.template(1..N) = datasets that are averaged into the standard
%
% The realignment is done by computing a minumum norm estimate using a
% large number of dipoles that are placed in the upper layer of the brain
% surface, followed by a forward computation towards the template
% gradiometer array. This requires the specification of a volume conduction
% model of the head and of a source model.
%
% A volume conduction model of the head should be specified with
%   cfg.headmodel   = structure, see FT_PREPARE_HEADMODEL
%
% A source model (i.e. a superficial layer with distributed sources) can be
% constructed from a headshape file, or from inner surface of the volume conduction
% model using FT_PREPARE_SOIURCEMODEL using the following options
%   cfg.spheremesh  = number of dipoles in the source layer (default = 642)
%   cfg.inwardshift = depth of the source layer relative to the headshape
%                     surface or volume conduction model (no default
%                     supplied, see below)
%   cfg.headshape   = a filename containing headshape, a structure containing a
%                     single triangulated boundary, or a Nx3 matrix with surface
%                     points
%
% If you specify a headshape and it describes the skin surface, you should specify an
% inward shift of 2.5 cm.
%
% For a single-sphere or a local-spheres volume conduction model based on the skin
% surface, an inward shift of 2.5 cm is reasonable.
%
% For a single-sphere or a local-spheres volume conduction model based on the brain
% surface, you should probably use an inward shift of about 1 cm.
%
% For a realistic single-shell volume conduction model based on the brain surface, you
% should probably use an inward shift of about 1 cm.
%
% Other options are
% cfg.pruneratio  = for singular values, default is 1e-3
% cfg.verify      = 'yes' or 'no', show the percentage difference (default = 'yes')
% cfg.feedback    = 'yes' or 'no' (default = 'no')
% cfg.channel     =  Nx1 cell-array with selection of channels (default = 'MEG'),
%                      see FT_CHANNELSELECTION for details
% cfg.trials      = 'all' or a selection given as a 1xN vector (default = 'all')
%
% This implements the method described by T.R. Knosche, Transformation
% of whole-head MEG recordings between different sensor positions.
% Biomed Tech (Berl). 2002 Mar;47(3):59-62. For more information and
% related methods, see Stolk et al., Online and offline tools for head
% movement compensation in MEG. NeuroImage, 2012.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_PREPARE_LOCALSPHERES, FT_PREPARE_SINGLESHELL

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data
ft_preamble trackconfig

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed',    {'plot3d',      'feedback'});
cfg = ft_checkconfig(cfg, 'renamedval', {'headshape',   'headmodel', []});
cfg = ft_checkconfig(cfg, 'required',   {'inwardshift', 'template'});
cfg = ft_checkconfig(cfg, 'renamed',    {'hdmfile',     'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'vol',         'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'grid',        'sourcemodel'});

% set the default configuration
cfg.headshape  = ft_getopt(cfg, 'headshape', []);
cfg.pruneratio = ft_getopt(cfg, 'pruneratio', 1e-3);
cfg.spheremesh = ft_getopt(cfg, 'spheremesh', 642);
cfg.verify     = ft_getopt(cfg, 'verify',     'yes');
cfg.feedback   = ft_getopt(cfg, 'feedback',   'yes');
cfg.trials     = ft_getopt(cfg, 'trials',     'all', 1);
cfg.channel    = ft_getopt(cfg, 'channel',    'MEG');
cfg.topoparam  = ft_getopt(cfg, 'topoparam',  'rms');

% store original datatype
dtype = ft_datatype(data);

% check if the input data is valid for this function
data = ft_checkdata(data, 'datatype', 'raw', 'feedback', 'yes', 'hassampleinfo', 'yes', 'ismeg', 'yes');

% do realignment per trial
pertrial = all(ismember({'nasX';'nasY';'nasZ';'lpaX';'lpaY';'lpaZ';'rpaX';'rpaY';'rpaZ'}, data.label));

% put the low-level options pertaining to the dipole grid in their own field
cfg = ft_checkconfig(cfg, 'renamed', {'tightgrid', 'tight'}); % this is moved to cfg.sourcemodel.tight by the subsequent createsubcfg
cfg = ft_checkconfig(cfg, 'renamed', {'sourceunits', 'unit'}); % this is moved to cfg.sourcemodel.unit by the subsequent createsubcfg

% put the low-level options pertaining to the sourcemodel in their own field
cfg = ft_checkconfig(cfg, 'createsubcfg', {'sourcemodel'});
% move some fields from cfg.sourcemodel back to the top-level configuration
cfg = ft_checkconfig(cfg, 'createtopcfg', {'sourcemodel'});

if isstruct(cfg.template)
  % this should be a cell-array
  cfg.template = {cfg.template};
end

% retain only the MEG channels in the data and temporarily store
% the rest, these will be added back to the transformed data later.

% select trials and channels of interest
tmpcfg = [];
tmpcfg.trials  = cfg.trials;
tmpcfg.channel = setdiff(data.label, ft_channelselection(cfg.channel, data.label));
rest = ft_selectdata(tmpcfg, data);

tmpcfg.channel = ft_channelselection(cfg.channel, data.label);
data = ft_selectdata(tmpcfg, data);

% restore the provenance information
[cfg, data] = rollback_provenance(cfg, data);

Ntrials = length(data.trial);

% cfg.channel = ft_channelselection(cfg.channel, data.label);
% dataindx = match_str(data.label, cfg.channel);
% restindx = setdiff(1:length(data.label),dataindx);
% if ~isempty(restindx)
%   fprintf('removing %d non-MEG channels from the data\n', length(restindx));
%   rest.label = data.label(restindx);    % first remember the rest
%   data.label = data.label(dataindx);    % then reduce the data
%   for i=1:Ntrials
%     rest.trial{i} = data.trial{i}(restindx,:);  % first remember the rest
%     data.trial{i} = data.trial{i}(dataindx,:);  % then reduce the data
%   end
% else
%   rest.label = {};
%   rest.trial = {};
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% construct the average template gradiometer array
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
template = struct([]); % initialize as 0x0 empty struct array with no fields
for i=1:length(cfg.template)
  if ischar(cfg.template{i})
    fprintf('reading template sensor position from %s\n', cfg.template{i});
    tmp = ft_read_sens(cfg.template{i}, 'senstype', 'meg');
  elseif isstruct(cfg.template{i}) && isfield(cfg.template{i}, 'coilpos') && isfield(cfg.template{i}, 'coilori') && isfield(cfg.template{i}, 'tra')
    tmp = cfg.template{i};
  elseif isstruct(cfg.template{i}) && isfield(cfg.template{i}, 'pnt') && isfield(cfg.template{i}, 'ori') && isfield(cfg.template{i}, 'tra')
    % it seems to be a pre-2011v1 type gradiometer structure, update it
    tmp = ft_datatype_sens(cfg.template{i});
  else
    ft_error('unrecognized template input');
  end
  % prevent "Subscripted assignment between dissimilar structures" error
  template = appendstruct(template, tmp); clear tmp
end

grad = ft_average_sens(template);

% construct the final template gradiometer definition
template = [];
template.grad = grad;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FT_PREPARE_VOL_SENS will match the data labels, the gradiometer labels and the
% volume model labels (in case of a localspheres model) and result in a gradiometer
% definition that only contains the gradiometers that are present in the data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

volcfg = [];
volcfg.headmodel = cfg.headmodel;
volcfg.grad      = data.grad;
volcfg.channel   = data.label; % this might be a subset of the MEG channels
gradorig         = data.grad;  % this is needed later on for plotting. As of
% yet the next step is not entirely correct, because it does not keep track
% of the balancing of the gradiometer array. FIXME this may require some
% thought because the leadfields are computed with low level functions and
% do not easily accommodate for matching the correct channels with each
% other (in order to compute the projection matrix).
[volold, data.grad] = prepare_headmodel(volcfg);

% note that it is necessary to keep the two volume conduction models
% separate, since the single-shell Nolte model contains gradiometer specific
% precomputed parameters. Note that this is not guaranteed to result in a
% good projection for local sphere models.
volcfg.grad    = template.grad;
volcfg.channel = 'MEG'; % include all MEG channels
[volnew, template.grad] = prepare_headmodel(volcfg);

if strcmp(ft_senstype(data.grad), ft_senstype(template.grad))
  [id, it] = match_str(data.grad.label, template.grad.label);
  fprintf('mean distance towards template gradiometers is %.2f %s\n', mean(sum((data.grad.chanpos(id,:)-template.grad.chanpos(it,:)).^2, 2).^0.5), template.grad.unit);
else
  % the projection is from one MEG system to another MEG system, which makes a comparison of the data difficult
  cfg.feedback = 'no';
  cfg.verify = 'no';
end

% copy all options that are potentially used in ft_prepare_sourcemodel
tmpcfg           = keepfields(cfg, {'sourcemodel', 'mri', 'headshape', 'symmetry', 'smooth', 'threshold', 'spheremesh', 'inwardshift', 'xgrid' 'ygrid', 'zgrid', 'resolution', 'tight', 'warpmni', 'template', 'showcallinfo'});
tmpcfg.headmodel = volold;
tmpcfg.grad      = data.grad;
% create the source positions on which the data will be projected
sourcemodel = ft_prepare_sourcemodel(tmpcfg);
pos = sourcemodel.pos;

% sometimes some of the dipole positions are nan, due to problems with the headsurface triangulation
% remove them to prevent problems with the forward computation
sel = find(any(isnan(pos(:,1)),2));
pos(sel,:) = [];

% compute the forward model for the new gradiometer positions
fprintf('computing forward model for %d dipoles\n', size(pos,1));
lfnew = ft_compute_leadfield(pos, template.grad, volnew);
if ~pertrial
  %this needs to be done only once
  lfold = ft_compute_leadfield(pos, data.grad, volold);
  [realign, noalign, bkalign] = computeprojection(lfold, lfnew, cfg.pruneratio, cfg.verify);
else
  %the forward model and realignment matrices have to be computed for each trial
  %this also goes for the singleshell volume conductor model
  %x = which('rigidbodyJM'); %this function is needed
  %if isempty(x)
  %  ft_error('you are trying out experimental code for which you need some extra functionality which is currently not in the release version of FieldTrip. if you are interested in trying it out, contact Jan-Mathijs');
  %end
end

% interpolate the data towards the template gradiometers
for i=1:Ntrials
  fprintf('realigning trial %d\n', i);
  if pertrial
    %warp the gradiometer array according to the motiontracking data
    sel   = match_str(rest.label, {'nasX';'nasY';'nasZ';'lpaX';'lpaY';'lpaZ';'rpaX';'rpaY';'rpaZ'});
    hmdat = rest.trial{i}(sel,:);
    if ~all(hmdat==repmat(hmdat(:,1),[1 size(hmdat,2)]))
      ft_error('only one position per trial is at present allowed');
    else
      %M    = rigidbodyJM(hmdat(:,1))
      M    = ft_headcoordinates(hmdat(1:3,1),hmdat(4:6,1),hmdat(7:9,1));
      grad = ft_transform_geometry(M, data.grad);
    end

    volcfg.grad = grad;
    %compute volume conductor
    [volold, grad] = prepare_headmodel(volcfg);
    %compute forward model
    lfold = ft_compute_leadfield(pos, grad, volold);
    %compute projection matrix
    [realign, noalign, bkalign] = computeprojection(lfold, lfnew, cfg.pruneratio, cfg.verify);
  end
  data.realign{i} = realign * data.trial{i};
  if strcmp(cfg.verify, 'yes')
    % also compute the residual variance when interpolating
    [id,it]   = match_str(data.grad.label, template.grad.label);
    rvrealign = rv(data.trial{i}(id,:), data.realign{i}(it,:));
    fprintf('original -> template             RV %.2f %%\n', 100 * mean(rvrealign));
    datnoalign = noalign * data.trial{i};
    datbkalign = bkalign * data.trial{i};
    rvnoalign = rv(data.trial{i}, datnoalign);
    rvbkalign = rv(data.trial{i}, datbkalign);
    fprintf('original             -> original RV %.2f %%\n', 100 * mean(rvnoalign));
    fprintf('original -> template -> original RV %.2f %%\n', 100 * mean(rvbkalign));
  end
end

% plot the topography before and after the realignment
if strcmp(cfg.feedback, 'yes')

  ft_warning('showing MEG topography (RMS value over time) in the first trial only');
  Nchan = length(data.grad.label);
  [id,it]   = match_str(data.grad.label, template.grad.label);
  pos1 = data.grad.chanpos(id,:);
  pos2 = template.grad.chanpos(it,:);
  prj1 = elproj(pos1); tri1 = delaunay(prj1(:,1), prj1(:,2));
  prj2 = elproj(pos2); tri2 = delaunay(prj2(:,1), prj2(:,2));

  switch cfg.topoparam
    case 'rms'
      p1 = sqrt(mean(data.trial{1}(id,:).^2, 2));
      p2 = sqrt(mean(data.realign{1}(it,:).^2, 2));
    case 'svd'
      [u, s, v] = svd(data.trial{1}(id,:)); p1 = u(:,1);
      [u, s, v] = svd(data.realign{1}(it,:)); p2 = u(:,1);
    otherwise
      ft_error('unsupported cfg.topoparam');
  end

  X = [pos1(:,1) pos2(:,1)]';
  Y = [pos1(:,2) pos2(:,2)]';
  Z = [pos1(:,3) pos2(:,3)]';

  % show figure with old an new helmets, volume model and source positions
  figure
  hold on
  ft_plot_headmodel(volold);
  plot3(sourcemodel.pos(:,1),sourcemodel.pos(:,2),sourcemodel.pos(:,3),'b.');
  plot3(pos1(:,1), pos1(:,2), pos1(:,3), 'r.') % original positions
  plot3(pos2(:,1), pos2(:,2), pos2(:,3), 'g.') % template positions
  line(X,Y,Z, 'color', 'black');
  view(-90, 90);

  % show figure with data on old helmet location
  figure
  hold on
  plot3(pos1(:,1), pos1(:,2), pos1(:,3), 'r.') % original positions
  plot3(pos2(:,1), pos2(:,2), pos2(:,3), 'g.') % template positions
  line(X,Y,Z, 'color', 'black');
  axis equal; axis vis3d
  bnd1 = [];
  bnd1.pos = pos1;
  bnd1.tri = tri1;
  ft_plot_mesh(bnd1,'vertexcolor',p1,'edgecolor','none')
  title('RMS, before realignment')
  view(-90, 90)

  % show figure with data on new helmet location
  figure
  hold on
  plot3(pos1(:,1), pos1(:,2), pos1(:,3), 'r.') % original positions
  plot3(pos2(:,1), pos2(:,2), pos2(:,3), 'g.') % template positions
  line(X,Y,Z, 'color', 'black');
  axis equal; axis vis3d
  bnd2 = [];
  bnd2.pos = pos2;
  bnd2.tri = tri2;
  ft_plot_mesh(bnd2,'vertexcolor',p2,'edgecolor','none')
  title('RMS, after realignment')
  view(-90, 90)
end

% store the realigned data in a new structure
interp.label   = template.grad.label;
interp.grad    = template.grad;   % replace with the template gradiometer array
interp.trial   = data.realign;    % remember the processed data
interp.fsample = data.fsample;
interp.time    = data.time;

% add the rest channels back to the data, these were not interpolated
if ~isempty(rest.label)
  fprintf('adding %d non-MEG channels back to the data (', length(rest.label));
  fprintf('%s, ', rest.label{1:end-1});
  fprintf('%s)\n', rest.label{end});
  for trial=1:length(rest.trial)
    interp.trial{trial} = [interp.trial{trial}; rest.trial{trial}];
  end
  interp.label = [interp.label; rest.label];
end

% copy the trial specific information into the output
if isfield(data, 'trialinfo')
  interp.trialinfo = data.trialinfo;
end

% copy the sampleinfo field as well
if isfield(data, 'sampleinfo')
  interp.sampleinfo = data.sampleinfo;
end

% convert back to input type if necessary
switch dtype
  case 'timelock'
    interp = ft_checkdata(interp, 'datatype', 'timelock');
  otherwise
    % keep the output as it is
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous data

% rename the output variable to accomodate the savevar postamble
data = interp;

ft_postamble provenance data
ft_postamble history    data
ft_postamble savevar    data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction that computes the projection matrix(ces)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [realign, noalign, bkalign] = computeprojection(lfold, lfnew, pruneratio, verify)

% compute this inverse only once, although it is used twice
tmp = prunedinv(lfold, pruneratio);
% compute the three interpolation matrices
fprintf('computing interpolation matrix #1\n');
realign = lfnew * tmp;
if strcmp(verify, 'yes')
  fprintf('computing interpolation matrix #2\n');
  noalign = lfold * tmp;
  fprintf('computing interpolation matrix #3\n');
  bkalign = lfold * prunedinv(lfnew, pruneratio) * realign;
else
  noalign = [];
  bkalign = [];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction that computes the inverse using a pruned SVD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [lfi] = prunedinv(lf, r)
[u, s, v] = svd(lf);
if r<1
  % treat r as a ratio
  p = find(s<(s(1,1)*r) & s~=0);
else
  % treat r as the number of spatial components to keep
  diagels = 1:(min(size(s))+1):(min(size(s)).^2);
  p       = diagels((r+1):end);
end
fprintf('pruning %d from %d, i.e. removing the %d smallest spatial components\n', length(p), min(size(s)), length(p));
s(p) = 0;
s(find(s~=0)) = 1./s(find(s~=0));
lfi = v * s' * u';

author: spm


text:
function s = SaturationPouli(C, I)
%
%       s = SaturationPouli(C, I)
%
%       This computes the saturation using channel C and I from LCh color
%       space
%
%       input:
%         - C: chroma channel from LCh
%         - I: intensity channel from LCh
%
%       output:
%         - S: saturation channel
%
%     Copyright (C) 2013  Francesco Banterle
%
%     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 3 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, see <http://www.gnu.org/licenses/>.
%
%     The paper describing this technique is:
%     "Color Correction for Tone Reproduction"
% 	  by Tania Pouli1, Alessandro Artusi, Francesco Banterle,
%     Ahmet Oguz Akyuz, Hans-Peter Seidel and Erik Reinhard
%     in the Twenty-first Color and Imaging Conference (CIC21), Albuquerque, Nov. 2013
%
%

D = sqrt(C.^2 + I.^2);
s = C ./ D;

end

author: banterle


text:
function Q1 = mpc1(T_meas,T_sp)

s = 'http://byu.apmonitor.com';
c = 'mpc';

% input measurement
apm_meas(s,c,'TC',T_meas);

% input setpoint with deadband +/- DT
DT = 0.1;
apm_option(s,c,'TC.sphi',T_sp+DT);
apm_option(s,c,'TC.splo',T_sp-DT);

% solve MPC
output = apm(s,c,'solve');

% test for successful solution
if (apm_tag(s,c,'apm.appstatus')==1)
    % retrieve the first Q value
    Q1 = apm_tag(s,c,'Q1.Newval');
else
    % display output for debugging
    print(output)
    % not successful, set voltage to zero
    Q1 = 0;
end

end

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function value = sin_power_int ( a, b, n )

%*****************************************************************************80
%
%% SIN_POWER_INT evaluates the sine power integral.
%
%  Discussion:
%
%    The function is defined by
%
%      SIN_POWER_INT(A,B,N) = Integral ( A <= T <= B ) ( sin ( t ))^n dt
%
%    The algorithm uses the following fact:
%
%      Integral sin^n ( t ) = (1/n) * (
%        sin^(n-1)(t) * cos(t) + ( n-1 ) * Integral sin^(n-2) ( t ) dt )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters
%
%    Input, real A, B, the limits of integration.
%
%    Input, integer N, the power of the sine function.
%
%    Output, real SIN_POWER_INT, the value of the integral.
%
  if ( n < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SIN_POWER_INT - Fatal error!\n' );
    fprintf ( 1, '  Power N < 0.\n' );
    value = 0.0;
    error ( 'SIN_POWER_INT - Fatal error!' );
  end

  sa = sin ( a );
  sb = sin ( b );
  ca = cos ( a );
  cb = cos ( b );

  if ( mod ( n, 2 ) == 0 )
    value = b - a;
    mlo = 2;
  else
    value = ca - cb;
    mlo = 3;
  end

  for m = mlo : 2 : n
    value = ( ( m - 1 ) * value + sa^(m-1) * ca - sb^(m-1) * cb ) / m;
  end

  return
end


text:
classdef uicontext < handle
%uicontext - Context menu with push buttons.

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

    properties(SetAccess = private)
        fig;        % The figure
        panel;      % The main panel
        items;      % The menu items
        gaps;       % Labels for covering the borders of items
        listener;	% Listener
    end
    methods
        %% Constructor
        function obj = uicontext(fig,width)
            obj.fig   = fig;
            obj.panel = uipanel(obj.fig,'Position',[0 0 width 0],'Visible',false);
        end
        %% Add a new item
        function add(obj,str,icon,cb)
        	obj.items = [obj.items,uibutton(obj.panel,'Position',[-2 0 obj.panel.Position(3)+5 25],'Text',str,'HorizontalAlignment','left','BackgroundColor',[.94 .94 .94],'Icon',icon,'ButtonPushedFcn',cb)];
            obj.gaps  = [obj.gaps,uipanel(obj.panel,'Position',[-2 0 obj.panel.Position(3)+5 3],'BorderType','none')];
            obj.panel.Children = obj.panel.Children([1,3:end,2]);
        end
        %% Flush the menu
        function flush(obj)
            if ~isempty(obj.items)
                obj.panel.Position(4) = 23*length(obj.items) - 2;
                loc = obj.panel.Position(4) + 1;
                for i = 1 : length(obj.items)
                    loc = loc - 23;
                    obj.items(i).Position(2) = loc;
                    obj.gaps(i).Position(2)  = loc - 1;
                end
            end
        end
        %% Show the menu
        function show(obj)
            if ~isempty(obj.items)
                if obj.fig.CurrentPoint(1)-10+obj.panel.Position(3) < obj.fig.Position(3)
                    obj.panel.Position(1) = obj.fig.CurrentPoint(1) - 10;
                else
                    obj.panel.Position(1) = obj.fig.CurrentPoint(1) + 10 - obj.panel.Position(3);
                end
                if obj.fig.CurrentPoint(2)+10-obj.panel.Position(4) > 0
                    obj.panel.Position(2) = obj.fig.CurrentPoint(2) + 10 - obj.panel.Position(4);
                else
                    obj.panel.Position(2) = obj.fig.CurrentPoint(2) - 10;
                end
                obj.panel.Visible = true;
                obj.listener      = obj.fig.addlistener('CurrentPoint','PostSet',@obj.cb_motion);
            end
        end
    end
    methods(Access = private)
        %% Hide the menu when moving out of it
        function cb_motion(obj,~,~)
            if obj.fig.CurrentPoint(1)<obj.panel.Position(1) || obj.fig.CurrentPoint(1)>sum(obj.panel.Position([1,3])) || obj.fig.CurrentPoint(2)<obj.panel.Position(2) || obj.fig.CurrentPoint(2)>sum(obj.panel.Position([2,4]))
                obj.panel.Visible = false;
                delete(obj.listener);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ result ] = infer_events_stove(result, events, times, threshold)

    appliance_name = 'Stove';
    phases = length(events);
    % check events that could be created by the stove of phase 1
    could_be_stove = abs(events{1}(:,1)) > threshold;
    times_could_be_stove = times{1}(could_be_stove);
    events_could_be_stove = events{1}(could_be_stove,1);
    for phase = 2:phases
        % check events that could be created by the stove of phase
        % "phase"
        could_be_stove = abs(events{phase}(:,1)) > threshold;
        times_phase = times{phase}(could_be_stove);
        % for each possible event on phase check if time lies within event on
        % old phase
        tmp = arrayfun(@(t) any(times_phase > t-5 & times_phase < t+5), times_could_be_stove);
        times_could_be_stove = times_could_be_stove(tmp);
        events_could_be_stove = events_could_be_stove(tmp);
    end

    result.appliance_names = {appliance_name};

    % events inferred from algorithm
    % col1: event time
    % col2: appliance id
    % col3: delta
    num_events = length(times_could_be_stove);
    result.events = [times_could_be_stove, ones(num_events,1), events_could_be_stove];
end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
% l1dantzig_pd.m
%
% Solves
% min_x  ||x||_1  subject to  ||A'(Ax-b)||_\infty <= epsilon
%
% Recast as linear program
% min_{x,u}  sum(u)  s.t.  x - u <= 0
%                         -x - u <= 0
%             A'(Ax-b) - epsilon <= 0
%            -A'(Ax-b) - epsilon <= 0
% and use primal-dual interior point method.
%
% Usage: xp = l1dantzig_pd(x0, A, At, b, epsilon, pdtol, pdmaxiter, cgtol, cgmaxiter)
%
% x0 - Nx1 vector, initial point.
%
% A - Either a handle to a function that takes a N vector and returns a K
%     vector , or a KxN matrix.  If A is a function handle, the algorithm
%     operates in "largescale" mode, solving the Newton systems via the
%     Conjugate Gradients algorithm.
%
% At - Handle to a function that takes a K vector and returns an N vector.
%      If A is a KxN matrix, At is ignored.
%
% b - Kx1 vector of observations.
%
% epsilon - scalar or Nx1 vector of correlation constraints
%
% pdtol - Tolerance for primal-dual algorithm (algorithm terminates if
%     the duality gap is less than pdtol).
%     Default = 1e-3.
%
% pdmaxiter - Maximum number of primal-dual iterations.
%     Default = 50.
%
% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
%     Default = 1e-8.
%
% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
%     if A is a matrix.
%     Default = 200.
%
% Written by: Justin Romberg, Caltech
% Email: jrom@acm.caltech.edu
% Created: October 2005
%

function xp = l1dantzig_pd(x0, A, At, b, epsilon, pdtol, pdmaxiter, cgtol, cgmaxiter)

largescale = isa(A,'function_handle');

if (nargin < 6), pdtol = 1e-3;  end
if (nargin < 7), pdmaxiter = 50;  end
if (nargin < 8), cgtol = 1e-8;  end
if (nargin < 9), cgmaxiter = 200;  end

N = length(x0);

alpha = 0.01;
beta = 0.5;
mu = 10;

gradf0 = [zeros(N,1); ones(N,1)];

% starting point --- make sure that it is feasible
if (largescale)
  if (max( abs(At(A(x0) - b)) - epsilon ) > 0)
    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
    AAt = @(z) A(At(z));
    [w, cgres] = cgsolve(AAt, b, cgtol, cgmaxiter, 0);
    if (cgres > 1/2)
      disp('A*At is ill-conditioned: cannot find starting point');
      xp = x0;
      return;
    end
    x0 = At(w);
  end
else
  if (max(abs(A'*(A*x0 - b)) - epsilon ) > 0)
    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
    opts.POSDEF = true; opts.SYM = true;
    [w, hcond] = linsolve(A*A', b, opts);
    if (hcond < 1e-14)
      disp('A*At is ill-conditioned: cannot find starting point');
      xp = x0;
      return;
    end
    x0 = A'*w;
  end
end
x = x0;
u = (0.95)*abs(x0) + (0.10)*max(abs(x0));

% set up for the first iteration
if (largescale)
  Atr = At(A(x) - b);
else
  Atr = A'*(A*x - b);
end
fu1 = x - u;
fu2 = -x - u;
fe1 = Atr - epsilon;
fe2 = -Atr - epsilon;
lamu1 = -(1./fu1);
lamu2 = -(1./fu2);
lame1 = -(1./fe1);
lame2 = -(1./fe2);
if (largescale)
  AtAv = At(A(lame1-lame2));
else
  AtAv = A'*(A*(lame1-lame2));
end

% sdg = surrogate duality gap
sdg = -[fu1; fu2; fe1; fe2]'*[lamu1; lamu2; lame1; lame2];
tau = mu*(4*N)/sdg;

% residuals
rdual = gradf0 + [lamu1-lamu2 + AtAv; -lamu1-lamu2];
rcent = -[lamu1.*fu1; lamu2.*fu2; lame1.*fe1; lame2.*fe2] - (1/tau);
resnorm = norm([rdual; rcent]);

% iterations
pditer = 0;
done = (sdg < pdtol) | (pditer >= pdmaxiter);
dispProgress('Newton', 0, pdmaxiter);
while (~done)

  % solve for step direction
  w2 = - 1 - (1/tau)*(1./fu1 + 1./fu2);

  sig11 = -lamu1./fu1 - lamu2./fu2;
  sig12 = lamu1./fu1 - lamu2./fu2;
  siga = -(lame1./fe1 + lame2./fe2);
  sigx = sig11 - sig12.^2./sig11;

  if (largescale)
    w1 = -(1/tau)*( At(A(1./fe2-1./fe1)) + 1./fu2 - 1./fu1 );
    w1p = w1 - (sig12./sig11).*w2;
    hpfun = @(z) At(A(siga.*At(A(z)))) + sigx.*z;
    [dx, cgres, cgiter] = cgsolve(hpfun, w1p, cgtol, cgmaxiter, 0);
    if (cgres > 1/2)
      disp('Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)');
      xp = x;
      return
    end
    AtAdx = At(A(dx));
  else
    w1 = -(1/tau)*( A'*(A*(1./fe2-1./fe1)) + 1./fu2 - 1./fu1 );
    w1p = w1 - (sig12./sig11).*w2;
    Hp = A'*(A*sparse(diag(siga))*A')*A + diag(sigx);
    opts.POSDEF = true; opts.SYM = true;
    [dx, hcond] = linsolve(Hp, w1p,opts);
    if (hcond < 1e-14)
      disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');
      xp = x;
      return
    end
    AtAdx = A'*(A*dx);
  end
  du = w2./sig11 - (sig12./sig11).*dx;

  dlamu1 = -(lamu1./fu1).*(dx-du) - lamu1 - (1/tau)*1./fu1;
  dlamu2 = -(lamu2./fu2).*(-dx-du) - lamu2 - (1/tau)*1./fu2;
  dlame1 = -(lame1./fe1).*(AtAdx) - lame1 - (1/tau)*1./fe1;
  dlame2 = -(lame2./fe2).*(-AtAdx) - lame2 - (1/tau)*1./fe2;
  if (largescale)
    AtAdv = At(A(dlame1-dlame2));
  else
    AtAdv = A'*(A*(dlame1-dlame2));
  end

  % find minimal step size that keeps ineq functions < 0, dual vars > 0
  iu1 = find(dlamu1 < 0); iu2 = find(dlamu2 < 0);
  ie1 = find(dlame1 < 0); ie2 = find(dlame2 < 0);
  ifu1 = find((dx-du) > 0); ifu2 = find((-dx-du) > 0);
  ife1 = find(AtAdx > 0); ife2 = find(-AtAdx > 0);
  smax = min(1,min([...
    -lamu1(iu1)./dlamu1(iu1); -lamu2(iu2)./dlamu2(iu2); ...
    -lame1(ie1)./dlame1(ie1); -lame2(ie2)./dlame2(ie2); ...
    -fu1(ifu1)./(dx(ifu1)-du(ifu1)); -fu2(ifu2)./(-dx(ifu2)-du(ifu2)); ...
    -fe1(ife1)./AtAdx(ife1); -fe2(ife2)./(-AtAdx(ife2)) ]));
  s = 0.99*smax;

  % backtracking line search
  suffdec = 0;
  backiter = 0;
  while (~suffdec)
    xp = x + s*dx;  up = u + s*du;
    Atrp = Atr + s*AtAdx;  AtAvp = AtAv + s*AtAdv;
    fu1p = fu1 + s*(dx-du);  fu2p = fu2 + s*(-dx-du);
    fe1p = fe1 + s*AtAdx;  fe2p = fe2 + s*(-AtAdx);
    lamu1p = lamu1 + s*dlamu1;  lamu2p = lamu2 + s*dlamu2;
    lame1p = lame1 + s*dlame1; lame2p = lame2 + s*dlame2;
    rdp = gradf0 + [lamu1p-lamu2p + AtAvp; -lamu1p-lamu2p];
    rcp = -[lamu1p.*fu1p; lamu2p.*fu2p; lame1p.*fe1p; lame2p.*fe2p] - (1/tau);
    suffdec = (norm([rdp; rcp]) <= (1-alpha*s)*resnorm);
    s = beta*s;
    backiter = backiter+1;
    if (backiter > 32)
      disp('Stuck backtracking, returning last iterate.  (See Section 4 of notes for more information.)')
      xp = x;
      return
    end
  end

  % setup for next iteration
  x = xp;  u = up;
  Atr = Atrp;  AtAv = AtAvp;
  fu1 = fu1p; fu2 = fu2p;
  fe1 = fe1p; fe2 = fe2p;
  lamu1 = lamu1p; lamu2 = lamu2p;
  lame1 = lame1p; lame2 = lame2p;

  sdg = -[fu1; fu2; fe1; fe2]'*[lamu1; lamu2; lame1; lame2];
  tau = mu*(4*N)/sdg;

  rdual = rdp;
  rcent = -[lamu1.*fu1; lamu2.*fu2; lame1.*fe1; lame2.*fe2] - (1/tau);
  resnorm = norm([rdual; rcent]);

  pditer = pditer+1;
  done = (sdg < pdtol) | (pditer >= pdmaxiter);
  dispProgress('Newton', pditer/pdmaxiter);

%   disp(sprintf('Iteration = %d, tau = %8.3e, Primal = %8.3e, PDGap = %8.3e, Dual res = %8.3e',...
%     pditer, tau, sum(u), sdg, norm(rdual)));
%   if (largescale)
%     disp(sprintf('                CG Res = %8.3e, CG Iter = %d', cgres, cgiter));
%   else
%     disp(sprintf('                  H11p condition number = %8.3e', hcond));
%   end

end
dispProgress('Newton', 'Close');
end

author: thomaskuestner


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
    = corn_run(fid, data, w, c, tc, opts)
% This program simulates the BK strategy on data
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
%    = corn_run(fid, data, w, c, tc, opts)
%
% cum_ret: a number representing the final cumulative wealth.
% cumprod_ret: cumulative return until each trading period
% daily_ret: individual returns for each trading period
% daily_portfolio: individual portfolio for each trading period
% exp_ret: experts' return
%
% data: market sequence vectors
% fid: handle for write log file
% w: window size
% c: correlation coefficient threshold
% tc: transaction cost rate parameter
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
%            = corn_run(fid, data, w, c, tc, opts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%[T, N]=size(data);
[n, m] = size(data);

% Variables for return, start with uniform weight
% cumprod_ret = 1;
% daily_ret = 1;
% weight = ones(nStocks, 1)/nStocks;

cum_ret = 1;
cumprod_ret = ones(n, 1);
daily_ret = ones(n, 1);
day_weight = ones(m, 1)/m;  %#ok<*NASGU>
day_weight_o = zeros(m, 1);
daily_portfolio = zeros(n, m);

% print file head
fprintf(fid, '-------------------------------------\n');
if (~opts.quiet_mode)
    fprintf(fid, 'Parameters [w=%d, c=%f, tc=%f\n]', w, c, tc);
    fprintf(fid, 'day\t Daily Return\t Total return\n');
end

for t = 1:1:n,
    % Calculate t's portfolio
    if (t >=2)
        [day_weight] = corn_kernel(data(1:t-1, :), w, c);
    end

    % Normalize the constraint
    day_weight = day_weight./sum(day_weight);
    daily_portfolio(t, :) = day_weight';

    % Cal t's return and total return
    daily_ret(t, 1) = (data(t, :)*day_weight)*(1-tc/2*sum(abs(day_weight-day_weight_o)));
    cum_ret = cum_ret * daily_ret(t, 1);
    cumprod_ret(t, 1) = cum_ret;

    day_weight_o = day_weight.*data(t, :)'/daily_ret(t, 1);

    % Debug information
    fprintf(fid, '%d\t%f\t%f\n', t, daily_ret(t, 1), cum_ret);
    if (~opts.quiet_mode)
        if (~mod(t, opts.display_interval)),
            fprintf(1, '%d\t%f\t%f\n', t, daily_ret(t, 1), cum_ret);
        end
    end
end

% Debug Information
fprintf(fid, 'CORN(w:%d, c:%.2f, tc:%.4f), Final return: %.2f\n', ...
    w, c, tc, cum_ret);
fprintf(fid, '-------------------------------------\n');

fprintf(1, 'CORN(w:%d, c:%.2f, tc:%.4f), Final return: %.2f\n', ...
    w, c, tc, cum_ret);
fprintf(fid, '-------------------------------------\n');
end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function score = score_eigenvalues(~,~,D)
%SCORE_EIGENVALUES - Eigenvalues as score

score = diag(D);
score = score(:);

author: bbci


text:
function val=jsonopt(key,default,varargin)
%
% val=jsonopt(key,default,optstruct)
%
% setting options based on a struct. The struct can be produced
% by varargin2struct from a list of 'param','value' pairs
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% $Id: loadjson.m 371 2012-06-20 12:43:06Z fangq $
%
% input:
%      key: a string with which one look up a value from a struct
%      default: if the key does not exist, return default
%      optstruct: a struct where each sub-field is a key
%
% output:
%      val: if key exists, val=optstruct.key; otherwise val=default
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONlab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

val=default;
if(nargin<=2) return; end
opt=varargin{1};
if(isstruct(opt))
    if(isfield(opt,key))
       val=getfield(opt,key);
    elseif(isfield(opt,lower(key)))
       val=getfield(opt,lower(key));
    end
end


text:
function z = aatx(a,x)
%AATX Implicitly compute A * A' * x for sptenmat.
%
%   Z = AATX(A,X) takes a sptenmat object A and computes A * A' *
%   X. This is done without converting A to a standard MATLAB sparse
%   matrix.
%
%   This function is likely most useful as an argument to a routine
%   such as EIGS.
%
%   See also SPTENMAT, SPTENSOR/EIGS.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

subs = a.subs;
s1 = subs(:,1);
s2 = subs(:,2);
m = size(a,1);
n = size(a,2);
vals = a.vals;

v1 = x(s1);
v1 = vals .* v1;
y = accumarray(s2, v1, [n 1]);

v2 = y(s2);
v2 = vals .* v2;
z = accumarray(s1, v2, [m 1]);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function c = mod(a,b)
	checkinputs(a,b);
	if numel(b)>1
		ind = b~=0;
		c = uint64mod(int64(a(ind)),int64(b(ind)));
		c(~ind) = a(~ind);
	else
		if b == 0
			c = a;
		else
			c = uint64mod(a,b);
		end
	end
end


text:
function [V,F] = png2mesh( ...
  filename, laplacian_smoothness_iterations, max_points_on_boundary)
  % PNG2MESH
  %
  % [V,F] = png2mesh(filename, laplacian_smoothness_iterations,
  % max_points_on_boundary)
  %
  % Mesh the non transparent part of an image. Best to have the entire shape
  % surrounded by at least one transparent pixel. Also the outer boundary and
  % boundary of holes should be roughly the same size as resampling to meet the
  % max_points_on_boundary parameter is done by edge collapsing on the boundary
  % (small boundaries will end up with too many samples and therefor the
  % triangulation will be dense there).
  %
  % Inputs:
  %   filename  path to .png file
  %   laplacian_smoothness_iterations  Number of iterations of laplacian
  %     smoothing of the positions of boundary curve points before handing
  %     curve to Triangle
  %   max_points_on_boundary  maximum number of points in the input curve
  %     before handing to Triangle
  % Outputs:
  %   F  #faces by 3, list of face indices
  %   V  #V by 2 list of polygon vertices, note that y-coordinates will be
  %     "flipped" with respect to the image if you think of the image rows as
  %     counting from 0 at the top left corner to size(im,1) at the bottom
  %     right corner
  %
  % Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch)
  %
  % See also: bwmesh, png2poly, triangle
  %

  warning('obsolete: use bwmesh');

  basename = regexprep(filename,'\.png$','');
  % use alpha channel of png image to define poly
  [V,E,H] = png2poly(...
    filename,...
    laplacian_smoothness_iterations, ...
    max_points_on_boundary);
  if ~isempty(H)
    warning('Holes non-empty, but I know holes sometimes come out broken');
  end

  unr = setdiff(1:size(V,1),E(:));
  if(~isempty(unr))
    warning('Unreferenced vertices in outline...\n');
  end
  % get average squared edge length as a guess at the maximum area constraint
  % for the triangulation
  avg_sqr_edge_length = mean(sum((V(E(:,1),:)-V(E(:,2),:)).^2,2));
  % triangulate the polygon
  [V,F] = triangle(V,E,H,'MaxArea',avg_sqr_edge_length/2.0,'Quality',30);
  %[V,F] = triangle(V,E,[]);
  unr = setdiff(1:size(V,1),F(:));
  if(~isempty(unr))
    warning('Removing unreferenced vertices in mesh...');
    [V,~,F] = faces_first(V,[],F);
    V=V(1:max(F(:)),:);
  end

end

author: alecjacobson

repo: gptoolbox


text:
%% Random Search Algorithm (Pure Random Search Algorithm)

% This code finds the minimum of f(x) = x(1)^2 + x(2)^2
% in which -5 < x(i) < 5
% This function is a convex and the minimum is at (0,0)
% The RSA is the simplest algorithm to solve optimization problem
% it is not efficient and it sometimes cannot solve the problem

clc
close all
clear all
dim=2;
popsize=100;
ftarget=0.01;
numIter=100;
ObjFun=@(x) sum(x.^2);
for i=1:numIter
    candidate=10*rand(dim,popsize)-5;
    best=min(feval(ObjFun,candidate));
    if best <= ftarget
        break;
    end
end
disp(best);

author: TheAlgorithms

repo: MATLAB-Octave


text:
clear;
netStruct = load('../data/res52_batch32_Rankloss_2:1:0.1_margin0.5_img0.75_shift_hard_256/net-epoch-60.mat');
net = dagnn.DagNN.loadobj(netStruct.net);
clear netStruct;
net.mode = 'test' ;
net.move('gpu') ;
net.removeLayer('RankLoss');  % If you test Stage I model, please comment it.
net.conserveMemory = true;
im_mean = reshape(net.meta.normalization.averageImage,1,1,3);
load('../dataset/Flickr30k-prepare/url_data.mat');

% val 2 test 3
p = imdb.images.data(imdb.images.set==3);
ff = [];
%%------------------------------
for i = 1:1000
    disp(i);
    str = [ p{i}];
    %im = imread(['./Flickr30k/flickr30k-images-300/',str(34:end)]);
    im = imread(['.',str]);
    oim = im;   % or oim = im;
    f = getFeature2(net,oim,im_mean,'data','fc1_1bn');
    f = sum(sum(f,1),2);
    f2 = getFeature2(net,fliplr(oim),im_mean,'data','fc1_1bn');
    f2 = sum(sum(f2,1),2);
    f = f+f2;
    size4 = size(f,4);
    f = reshape(f,[],size4)';
    f = norm_zzd(f);
    ff = cat(1,ff,f);
end

save('../test/resnet_flikr30k_pool5_img.mat','ff','-v7.3');
%}

ff = [];
load('../dataset/Flickr30k-prepare/dense_feature_word2.1.mat');
% val 2 test 3
test_set = find(imdb.images.set==3);

for i = 1:1000
    disp(i);
    for j=1:5
        index = 5*test_set(i)-5+j;
        %get vector
        content = wordcnn(:,index);
        len = sum(wordcnn(:,index)>0);
        txtinput = zeros(len,20074,'single');
        for k=1:len %32
            txtinput(k,content(k))=1;
        end
        %transfer it to different location
        win = 33-len;
        input = zeros(32,20074,win,'single');
        for kk = 1:win
            input(kk:kk+len-1,:,kk) = txtinput;
        end

        input = reshape(input,1,32,20074,[]);
        f = getFeature2(net,input,[],'data2','fc5_2bn');
        f = sum(f,4);
        size4 = size(f,4);
        f = reshape(f,[],size4)';
        f = norm_zzd(f);
        ff = cat(1,ff,f);
    end
end
save('../test/resnet_flikr30k_pool5_txt.mat','ff','-v7.3');
%}
evaluate;

author: layumi

repo: Image-Text-Embedding


text:
function varargout= bbci_acquire_lsl(varargin)
%BBCI_ACQUIRE_lsl - Online data acquisition from labstreaming layer (LSL)
%   This function acquires small blocks of signals including meta
%   information from the LSL. The LSL can hold streams from any common
%   devices that will be accessed by type: 'eeg' or 'marker'. For details
%   see https://code.google.com/p/labstreaminglayer/.
%
%Synopsis:
%  STATE= bbci_acquire_XYZ('init', <PARAM>)
%  [CNTX, MRKTIME, MRKDESC, STATE]= bbci_acquire_XYZ(STATE)
%  bbci_acquire_XYZ('close')
%  bbci_acquire_XYZ('close', STATE)
%
%Arguments:
%  PARAM - Optional arguments, specific to XYZ.
%
%Output:
%  STATE - Structure characterizing the incoming signals; fields:
%     'fs', 'clab', and intern stuff
%  CNTX - 'acquired' signals [Time x Channels]
%  The following variables hold the markers that have been 'acquired' within
%  the current block (if any).
%  MRKTIME - DOUBLE: [1 nMarkers] position [msec] within data block.
%      A marker occurrence within the first sample would give
%      MARTIME= 1/STATE.fs.
%  MRKDESC - CELL {1 nMarkers} descriptors like 'S 52'

% 11-2015 Jan Boelts
% --- --- --- ---

global BTB

% initialization of state structure and LSL streams
if isequal(varargin{1}, 'init'),

    % check whether lsl toolbox is on the path
    if ~isdir('liblsl-Matlab')
        error('LSL Toolbox is not on the path. add it via addpath(genpath(''path_to_LSL/liblsl-Matlab''))')
    end

    % use default electrode setting
    state= opt_proplistToStruct(varargin{2:end});
    default_clab= ...
        {'AF5' 'AF3' 'AF1' 'AFz' 'AF2' 'AF4' 'AF6' ...
        'F5' 'F3', 'F1' 'Fz' 'F2' 'F4' 'F6' ...
        'FC7' 'FC5' 'FC3' 'FC1' 'FCz' 'FC2' 'FC4' 'FC6' 'FC8' ...
        'T7' 'C5' 'C3' 'C1' 'Cz' 'C2' 'C4' 'C6' 'T8' ...
        'CP7' 'CP5' 'CP3' 'CP1' 'CPz' 'CP2' 'CP4' 'CP6' 'CP8' ...
        'P7' 'P5' 'P3' 'P1' 'Pz' 'P2' 'P4' 'P6' 'P8' ...
        'PO5' 'PO3' 'PO1' 'POz' 'PO2' 'PO4' 'PO6' ...
        'O5' 'O3' 'O1' 'Oz' 'O2' 'O4' 'O6' ...
        };
    % set default parameters that may adapted from the stream later on
    props= {'fs'            100            '!DOUBLE[1]'
        'clab'          default_clab   'CELL{CHAR}'
        'blocksize'     40             '!DOUBLE[1]'
        'port'          'COM11'        '!CHAR'
        'timeout'       3              '!DOUBLE[1]'
        'filtHd'        []             'STRUCT'
        'verbose'       true           '!BOOL'
        };
    [state, isdefault]= opt_setDefaults(state, props, 1);

    % set default filter coeffs
    if isdefault.filtHd,
        % Fs/4
        filt1.b= [0.85 0 0.85];
        filt1.a= [1 0 0.7];
        %Fs/2
        filt2.b= [0.8 0.8];
        filt2.a= [1 0.6];
        state.filtHd= procutil_catFilters(filt1, filt2);
        state.filtHd.PersistentMemory= true;
    end
    % set number of channels and corresponding values
    state.nChans= length(state.clab);
    state.nBytesPerPacket= 2+3*state.nChans+4;
    nPacketsPerPoll= ceil(state.blocksize/1000*state.fs);
    state.nBytesPerPoll= nPacketsPerPoll*state.nBytesPerPacket;

    %%%%%  resolve eeg stream from LSL %%%%%
    eeg = {};
    % load a lsl library
    state.lib = lsl_loadlib();
    % look for the stream several times
    for i=1:3
        % look for stream on the network
        eeg = lsl_resolve_byprop(state.lib,'type','EEG');
        if ~isempty(eeg)
            break
        end
        pause(0.1)
    end
    if isempty(eeg)
        error('No LSL EEG stream on the network')
    else
        % create a new inlet
        % save lsl structures to state structure
        state.inlet.x = lsl_inlet(eeg{1});
        [~, state.starttime] = state.inlet.x.pull_sample();
        state.running = 1;
        % get the stream info object for eeg stream
        eeg_info = state.inlet.x.info();

        % set sampling rate of current stream
        if not(state.fs==eeg_info.nominal_srate)
            state.fs = eeg_info.nominal_srate;
            warning('EEG sampling rate is different from default: was adapted from 100 to %i', state.fs)
        end
        % check number of channels
        if not(state.nChans==eeg_info.channel_count())
            state.nChans = eeg_info.channel_count();
            warning('EEG nChans is different from default')
        end
        state.packetNo=[];
        state.buffer= [];
        state.lastx= zeros(1, state.nChans);
        state.scale= 1000000/2^24;
    end

    % resolve marker stream, try several times
    mrks = {};
    for i=1:3
        mrks = lsl_resolve_byprop(state.lib, 'name', 'MyMarkerStream', 1, 1);
        if ~isempty(mrks)
            break
        end
        pause(0.1)
    end
    if isempty(mrks)
        error('No LSL marker stream with name ''MyMarkerStream'' on the network, did set up the marker stream?')
    else
        state.inlet.mrk = lsl_inlet(mrks{1});
        state.lastMrkDesc= 256;
    end

    if isempty(state.filtHd),
        reset(state.filtHd);
    end
    output= {state};
% close condition needs the 'state' structure.
elseif isequal(varargin{1}, 'close') && length(varargin)==1
    error('Please use ''close'' option with ''state'' variable as second argument: bbci_lsl_acquire(close, state)');
elseif isequal(varargin{1}, 'close') && isstruct(varargin{2}),
    % close inlets and libraries
    state = varargin{2};
    state.inlet.x.delete();
    state.inlet.mrk.delete();
    %state.lib.delete();
    output= {state};

elseif length(varargin)~=1,
    error('Except for INIT/CLOSE case, only one input argument expected');
else % this is the running condition that receives and returns the samples
    if ~isstruct(varargin{1}),
        error('First input argument must be ''init'', ''close'', or a struct');
    end
    state= varargin{1};

    % get data sample from the inlet
    % set timeout to zero to prevent blocking if there is no marker in the
    % current sample.
    timeout = 0; % in seconds
    [cntx, ~] = state.inlet.x.pull_sample();

    % get marker and time
    [mrkDesc, mrkTime] = state.inlet.mrk.pull_sample(timeout);
    % lsl gives mrkTime in relation to the total time of the stream, not
    % the recording. therefore we reset the time to the beginning of the
    % recording. Maybe unneccessary and might be removed?
    mrkTime = mrkTime-state.starttime;

    state.lastMrkDesc= mrkDesc;

    % check whether streams are still on
    % if the timeout was exceeded mrkTime will be empty
    state.running = not(isempty(cntx));

    % save most recent sample
    state.lastx= cntx;

    output = {cntx, mrkTime, mrkDesc, state};
end
varargout= output(1:nargout);


text:
function imgOut = CIECAM02_ChromaticAdaptation(imgXYZ, img_white_XYZ)
%
%
%       imgOut = CIECAM02_ChromaticAdaptation(imgXYZ, img_white_XYZ)
%
%       Input:
%           -imgXYZ: an image in the XYZ color space.
%           -img_white_XYZ: is imgXYZ white point image in the XYZ color space.
%
%       Output:
%           -imgOut: is imgXYZ in the XYZ color space with chromatic
%           adaption.
%
%     Copyright (C) 2015  Francesco Banterle
%
%     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 3 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, see <http://www.gnu.org/licenses/>.
%
%     The paper describing this technique is:
%     "The CIECAM02 color appearance model"
% 	  by Nathan Moroney , Mark D. Fairchild , Robert W. G. Hunt ,
%     Changjun Li , M. Ronnier Luo , Todd Newman
%     in IS&T/SID 10 th Color Imaging Conference
%

M_CAT02 = [ 0.7328 0.4296 -0.1624;...
           -0.7036 1.6975  0.0061;...
            0.003  0.0136  0.9834];

if(exist('img_white_XYZ', 'var'))
    img_RGB_w = ConvertLinearSpace(img_white_XYZ, M_CAT02);
    L_A = 0.2 * img_white_XYZ(:,:,2); %adaptation luminance
else
    img_RGB_w = ones(1,1,3);
    L_A = 1.0;
end

imgRGB = ConvertLinearSpace(imgXYZ, M_CAT02);

D = CIECAM02_DegreeAdaptation(L_A);

wp_D65_XYZ = [96.047, 100, 108.883]; %D65 white point in XYZ
wp_D65_RGB = ConvertLinearSpace(reshape(wp_D65_XYZ, 1, 1, 3), M_CAT02);

imgRGB_c = zeros(size(imgRGB));

for i=1:size(imgXYZ, 3)
    imgRGB_c(:,:,i) = imgRGB(:,:,i) .* (wp_D65_RGB(i) .* D ./ img_RGB_w(:,:,i) + (1.0 - D));
end

imgOut = ConvertLinearSpace(imgRGB_c, inv(M_CAT02));

end