text:
% FT_POSTAMBLE_SAVEFIG is a helper script that optionally saves the output
% MATLAB figure to a *.fig and to a *.png file on disk. This is useful for
% batching and for distributed processing. This makes use of the
% cfg.outputfile variable.
%
% Use as
%   ft_postamble savefig
%
% See also FT_PREAMBLE, FT_POSTAMBLE, FT_POSTAMBLE_SAVEVAR

% Copyright (C) 2019, Robert Oostenveld, DCCN
%
% 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$

% the output data should be saved to a MATLAB file
if (isfield(cfg, 'outputfile') && ~isempty(cfg.outputfile)) || exist('Fief7bee_reproducescript', 'var')

  if isfield(cfg, 'outputlock') && ~isempty(cfg.outputlock)
    mutexlock(cfg.outputlock);
  end

  if exist('Fief7bee_reproducescript', 'var')
    % write the output figure to a MATLAB file
    iW1aenge_now = datestr(now, 30);
    cfg.outputfile = fullfile(Fief7bee_reproducescript, sprintf('%s_%s_output',...
      iW1aenge_now, FjmoT6aA_highest_ft)); % the file extension is added later

    % write the large configuration fields to a MATLAB file
    % this applies to layout, event, sourcemodel, headmodel, grad, etc.
    % note that this is here, rather than in the (seemingly more logical)
    % ft_preamble_loadvar, because this code depends on cfg.callinfo (which
    % is only present at postamble stage)
    cfg = save_large_cfg_fields(cfg, FjmoT6aA_highest_ft, Fief7bee_reproducescript, iW1aenge_now);

    % write a snippet of MATLAB code with the user-specified configuration and function call
    reproducescript(Fief7bee_reproducescript, fullfile(Fief7bee_reproducescript, 'script.m'),...
      FjmoT6aA_highest_ft, cfg, false);

  elseif (isfield(cfg, 'outputfile') && ~isempty(cfg.outputfile))
    % keep the output file as it is

  else
    % don't write to an output file
    cfg.outputfile = {};
  end

  % save the current figure to a MATLAB .fig and to a .png file
  if ~isempty(cfg.outputfile)
    [Fief7bee_p, Fief7bee_f, Fief7bee_x] = fileparts(cfg.outputfile);
    Fief7bee_outputfile = fullfile(Fief7bee_p, [Fief7bee_f, '.fig']);
    ft_info('writing output figure to file ''%s''\n', Fief7bee_outputfile);
    savefig(gcf, cfg.outputfile, 'compact')
    Fief7bee_outputfile = fullfile(Fief7bee_p, [Fief7bee_f, '.png']);
    ft_info('writing screenshot to file ''%s''\n', Fief7bee_outputfile);
    set(gcf, 'PaperOrientation', 'portrait');
    print(Fief7bee_outputfile, '-dpng');
  end

  if isfield(cfg, 'outputlock') && ~isempty(cfg.outputlock)
    mutexunlock(cfg.outputlock);
  end

end


text:
function areaIntersection = computeIntersectionArea(bb1,bb2)
%compute intersection anrea of bb1 and bb2
%bb1 and bb2 - bounding boxes
%bbi = [xmin ymin xmax ymax] for i=1,2

xmin = max(bb1(1),bb2(1));
xmax = min(bb1(3),bb2(3));
ymin = max(bb1(2),bb2(2));
ymax = min(bb1(4),bb2(4));

areaIntersection = computeArea([xmin ymin xmax ymax]);

end

author: Cloud-CV

repo: object-proposals


text:
function p00_poly_write ( test, file_name )

%*****************************************************************************80
%
%% P00_POLY_WRITE collects data and writes it to a POLY file.
%
%  Discussion:
%
%    This routine collects information about the boundary for a given
%    problem, and writes that data to a POLY file that can be read by
%    TRIANGLE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TEST, the index of the test problem
%
%    Input, string FILE_NAME, the name of the file to create.
%
  dim_num = 2;
%
%  Get a length scale.
%
  [ lo, hi ] = p00_box ( test, dim_num );

  scale = max ( abs ( hi(1:dim_num) - lo(1:dim_num) ) );
%
%  How many boundary segments are there?
%
  boundary_segment_num = p00_boundary_segment_num ( test );
%
%  Choose H so that we would expect to get about 100 boundary points if our
%  region were a square of any size.
%
  h = 0.04 * scale;

  node_num = 0;
  for segment_index = 1 : boundary_segment_num
    segment_length = p00_boundary_segment_length ( test, segment_index, h );
    node_num = node_num + segment_length;
  end
%
%  Now collect all the boundary nodes into one array.
%
  next = 1;

  for segment_index = 1 : boundary_segment_num

    segment_length = p00_boundary_segment_length ( test, segment_index, h );

    segment(1:dim_num,next:next+segment_length-1) = p00_boundary_segment ( ...
      test, segment_index, dim_num, segment_length );

    next = next + segment_length;

  end
%
%  How many edges are there?
%
  edge_num = 0;

  for segment_index = 1 : boundary_segment_num

    segment_length = p00_boundary_segment_length ( test, segment_index, h );

    edge_num = edge_num + segment_length - 1;

  end
%
%  Now collect the edges.
%
  e = 0;
  next = 1;

  for segment_index = 1 : boundary_segment_num

    segment_length = p00_boundary_segment_length ( test, segment_index, h );

    for j = 1 : segment_length - 1
      edge_nodes(1,e+j) = next + j - 1;
      edge_nodes(2,e+j) = next + j;
    end

    next = next + segment_length;
    e = e + segment_length - 1;

  end
%
%  Handle the holes.
%
  hole_num = p00_hole_num ( test );
  hole = zeros ( dim_num, hole_num );

  for hole_index = 1 : hole_num
    hole(1:dim_num,hole_index) = ( p00_hole_point ( test, hole_index, dim_num ) )';
  end
%
%  Write the POLY file.
%
  poly_write ( file_name, node_num, segment, edge_num, ...
    edge_nodes, hole_num, hole );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [newdata] = mne_compensate_to(data,to)
%
% [newdata] = mne_compensate_to(data,to)
%
% Apply compensation to the data as desired
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%   Revision 1.4  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.3  2006/04/18 20:44:46  msh
%   Added reading of forward solution.
%   Use length instead of size when appropriate
%
%   Revision 1.2  2006/04/14 15:49:49  msh
%   Improved the channel selection code and added ch_names to measurement info.
%
%   Revision 1.1  2006/04/12 10:51:19  msh
%   Added projection writing and compensation routines
%
%

me='MNE:mne_compensate_to';

newdata = data;
now     = mne_get_current_comp(newdata.info);
%
%   Are we there already?
%
if now == to
    fprintf(1,'Data are already compensated as desired\n');
end
%
%   Make the compensator and apply it to all data sets
%
comp = mne_make_compensator(newdata.info,now,to);
for k = 1:length(newdata.evoked)
   newdata.evoked(k).epochs = comp*newdata.evoked(k).epochs;
end
%
%  Update the compensation info in the channel descriptors
%
newdata.info.chs = mne_set_current_comp(newdata.info.chs,to);
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [fig_name_tca]=ps_mean_v(ifg_list,n_boot,subtract_switches,use_small_baselines,aps_flag)
%PS_MEAN_V Calculate mean velocities and their standard deviations
%   PS_MEAN_V(IFG_LIST,N_BOOT,SUBTRACT_SWITCHES) where IFG_LIST gives indices
%   of interferograms to be included in the calculation,N_BOOT specifies number
%   of bootstrap iterations and RAMP_FLAG is 1 to subtract orbital ramps
%
%   Andy Hooper, Jan 2008
%
%   ======================================================================
%   06/2009 AH: Subtract orbital ramp if present
%   06/2009 AH: Process smaller chunks to reduce memory needs
%   02/2010 AH: Revert to v3.1 version
%   02/2010 AH: Replace unwrap_ifg_index with drop_ifg_index
%   03/2010 AH: 3rd input changed to 'subtract_switches'
%   03/2010 AH: Use small baselines option added
%   03/2010 AH: var/cov added to inversion
%   05/2010 AH: Include master if there are not ifgs before and after
%   06/2010 AH: small change to ps_mean_v.m
%   03/2014 AH: -a etc added
%   04/2014 DB: fix fig_name_tca output in case no aps correction
%   05/2014 DB: Fix v-doa to v-dao
%   05/2014 DB: For APS related corrections include deramping on the fly
%   05/2014 DB: big fix in case of a nan
%   07/2014 EH: Fix for nanmean in case of single element compared to vector
%   03/2015 DB: Do all deramping on the fly, make consistent for TRAIN release
%   05/2017 DB: use stamps save to save larger variables than 2GB.
%   06/2017 DB: Include the option to choose between chol decompostion for
%               bootstrapping (not to invert covariance often), or conventional
%               bootstrapping requiring to invert covariance matrix for each iteration.
%   ======================================================================

fprintf('Calculating standard deviation of mean velocity...\n')

chol_flag = 'n';

if nargin<1
    ifg_list=[];
end
if nargin<2
    n_boot=100;
end

if nargin<3 | isempty(subtract_switches)
   subtract_switches='';
end

if nargin<4
    use_small_baselines=0;
end

if nargin<5
    aps_flag=1;
end

load psver
psname=['./ps',num2str(psver)];
phuwname=['./phuw',num2str(psver)];
phuwsbresname=['./phuw_sb_res',num2str(psver)];
ifgstdname=['./ifgstd',num2str(psver)];
sclaname=['./scla',num2str(psver)];
apsname=['./tca',num2str(psver)];
tidename=['./tide',num2str(psver)];

mvname=['mv',num2str(psver)];

ps=load(psname);

drop_ifg_index=getparm('drop_ifg_index');
fig_name_tca = '';
if strcmpi(getparm('small_baseline_flag'),'y')
    if use_small_baselines==0
        phuw=load(phuwname,'unwrap_ifg_index_sm');
        if isfield(phuw,'unwrap_ifg_index_sm');
            unwrap_ifg_index=phuw.unwrap_ifg_index_sm;
        else
            unwrap_ifg_index=[1:ps.n_image];
        end
        phuwres=load(phuwsbresname,'sm_cov');
        if isfield(phuwres,'sm_cov');
            ifg_cov=phuwres.sm_cov;
        else
            ifg_cov=eye(ps.n_image);
        end
    else
        unwrap_ifg_index=setdiff([1:ps.n_ifg],drop_ifg_index);
        phuwname=['./phuw_sb',num2str(psver)];
        sclaname=['./scla_sb',num2str(psver)];
        apsname=['./tca_sb',num2str(psver)];
        tidename=['./tide_sb',num2str(psver)];

        phuwres=load(phuwsbresname,'sb_cov');
        if isfield(phuwres,'sb_cov');
            ifg_cov=phuwres.sb_cov;
        else
            ifg_cov=eye(ps.n_ifg);
        end
    end
else
    use_small_baselines=0;
    unwrap_ifg_index=setdiff([1:ps.n_ifg],drop_ifg_index);
    if ~exist([ifgstdname,'.mat'],'file')
        ifg_cov=eye(ps.n_ifg);
    else ifgstd=load(ifgstdname);
      if isfield(ifgstd,'ifg_std');
        ifgvar=(ifgstd.ifg_std*pi/181).^2;
        ifg_cov=diag(ifgvar);
      else
        ifg_cov=eye(ps.n_ifg);
      end
    end
end

uw=load(phuwname);
ph_uw=uw.ph_uw;
clear uw

switch(subtract_switches)
case('')
case('d')
    scla=load(sclaname);
    ph_uw=ph_uw - scla.ph_scla;
    clear scla
case('o')
    [ph_uw] = ps_deramp(ps,ph_uw);
case('a')
    aps=load(apsname);
    [aps_corr,fig_name_tca] = ps_plot_tca(aps,aps_flag);
    ph_uw=ph_uw - aps_corr;
    clear aps aps_corr
case('ao')
    aps=load(apsname);
    [aps_corr,fig_name_tca] = ps_plot_tca(aps,aps_flag);
    ph_uw=ph_uw - aps_corr;
    [ph_uw] = ps_deramp(ps,ph_uw);
    clear aps aps_corr
case('do')
    scla=load(sclaname);
    ph_uw=ph_uw - scla.ph_scla;
    [ph_uw] = ps_deramp(ps,ph_uw);
case('da')
    scla=load(sclaname);
    ph_uw=ph_uw - scla.ph_scla;
    clear scla
    aps=load(apsname);
    [aps_corr,fig_name_tca] = ps_plot_tca(aps,aps_flag);
    ph_uw=ph_uw - aps_corr;
    clear aps aps_corr
case('dat')
    tide=load(tidename);
    scla=load(sclaname);
    ph_uw=ph_uw - scla.ph_scla;
    clear scla
    aps=load(apsname);
    [aps_corr,fig_name_tca] = ps_plot_tca(aps,aps_flag);
    ph_uw=ph_uw - aps_corr-tide.ph_tide;
    clear aps aps_corr
case('dao')
    scla=load(sclaname);
    ph_uw=ph_uw - scla.ph_scla ;
    [ph_uw] = ps_deramp(ps,ph_uw);
    aps=load(apsname);
    [aps_corr,fig_name_tca] = ps_plot_tca(aps,aps_flag);
    ph_uw=ph_uw - aps_corr;
    clear aps aps_corr
otherwise
        error('unknown subtract flags')
end

if use_small_baselines==0
    %AH2CHECK %%%%%%% & unwrap_ifg_index(1)~=ps.master_ix & unwrap_ifg_index(end)~=ps.master_ix
    unwrap_ifg_index=setdiff(unwrap_ifg_index,ps.master_ix);
end

ph_all=zeros(ps.n_ps,1);
ref_ps=ps_setref;
if ~isempty(ifg_list)
    unwrap_ifg_index=intersect(unwrap_ifg_index,ifg_list);
    ifg_list=[];
end
ph_uw=ph_uw(:,unwrap_ifg_index);
ifg_cov=ifg_cov(unwrap_ifg_index,unwrap_ifg_index);
if rank(ifg_cov)<size(ifg_cov,1)
    error('Interferograms included for which there is no unwrapped solution')
end

if use_small_baselines==0
    day=ps.day(unwrap_ifg_index)-ps.master_day;
else
    day=ps.day(ps.ifgday_ix(unwrap_ifg_index,2))-ps.day(ps.ifgday_ix(unwrap_ifg_index,1));
end
N=length(day);
ph_uw=ph_uw-repmat(nanmean(ph_uw(ref_ps,:),1),ps.n_ps,1);
lambda=getparm('lambda');
ph_uw=double(ph_uw/4/pi*lambda*1000)';
G=[ones(N,1),double(day)/365.25] ;

% computing mean velocity
mean_v=lscov(G,ph_uw,ifg_cov)';
mean_v=mean_v(:,2);

% perform the bootstrapping
n=10000; % process n PS at a time
i=1;
mean_v_std=NaN(ps.n_ps,1,'single');
rand_init=sum(100*clock);
n_ifg = size(ph_uw,1);

if strcmpi(chol_flag,'y')
    % cholesky decomposition of the set of linear equations including
    % covariance. This allows to use the decomposition values in a uniform
    % LSQ inversion. i.e. the covariance matrix does not need to be
    % inverted each time.
    ph_uw=chol(inv(ifg_cov))*ph_uw;
    G=chol(inv(ifg_cov))*G;

    % looping over the PS segments
    while i<ps.n_ps
        if i+n< ps.n_ps
           i_end=i+n-1;
        else
           i_end=ps.n_ps;
        end
        ph_bit=ph_uw(:,i:i_end);

        % in-case of nans and only a subset being processed
        ix_non_nan = sum(isnan(ph_bit),1)==0;

        % fix the random generator, such its the same for each segment
        rand('twister',rand_init);

        % skip the segment in case of nan's
        if sum(ix_non_nan)>0
            if sum(ix_non_nan)== length(ix_non_nan)
                [mean_v_dist,boot_ix] = bootstrp(n_boot, @(x) single(lscov(G(x,:),ph_bit(x,:))), [1:N]);
                 temp = std(mean_v_dist(:,2:2:end))';
            else
                [mean_v_dist_temp,boot_ix] = bootstrp(n_boot, @(x) single(lscov(G(x,:),ph_bit(x,ix_non_nan))), [1:N]);
                temp = NaN([size(ph_bit,2) 1]);
                temp(ix_non_nan) = std(mean_v_dist_temp(:,2:2:end))';
                clear mean_v_dist_temp
            end
            mean_v_std(i:i_end) =temp;
        end
        i=i+n;
        fprintf('%d PS processed\n',i-1)
    end

else
    % conventional bootstrapping where the coveriance matrix is inverted on
    % each set PS segments.
    while i<ps.n_ps
        if i+n< ps.n_ps
           i_end=i+n-1;
        else
           i_end=ps.n_ps;
        end
        ph_bit=ph_uw(:,i:i_end);

        % in-case of nans and only a subset being processed
        ix_non_nan = sum(isnan(ph_bit),1)==0;

        % fix the random generator, such its the same for each segment
        rand('twister',rand_init)

        % skip the segment in case of nan's
        if sum(ix_non_nan)>0
            if sum(ix_non_nan)== length(ix_non_nan)
                [mean_v_dist,boot_ix] = bootstrp(n_boot, @(x) single(lscov(G(x,:),ph_bit(x,:),ifg_cov)), [1:N]);
                temp = std(mean_v_dist(:,2:2:end))';
            else
               fprintf('test\n')
               [mean_v_dist_temp,boot_ix] = bootstrp(n_boot, @(x) single(lscov(G(x,:),ph_bit(x,ix_non_nan),ifg_cov(ix_non_nan,ix_non_nan))), [1:N]);
               temp = NaN([size(ph_bit,2) 1]);
               temp(ix_non_nan) = std(mean_v_dist_temp(:,2:2:end))';
               clear mean_v_dist_temp
            end
            mean_v_std(i:i_end) =temp;
        end
        i=i+n;
        fprintf('%d PS processed\n',i-1)
    end
end

stamps_save(mvname,n_boot,subtract_switches,mean_v,mean_v_std)

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function sin_power_int_test ( )

%*****************************************************************************80
%
%% SIN_POWER_INT_TEST demonstrates the use of SIN_POWER_INT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SIN_POWER_INT_TEST\n' );
  fprintf ( 1, '  SIN_POWER_INT returns values of\n' );
  fprintf ( 1, '  the integral of SIN(X)^N from A to B.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      A         B          N      Exact           Computed\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, a, b, n, fx ] = sin_power_int_values ( n_data );

    if ( n_data == 0 )
      break;
    end

    fx2 = sin_power_int ( a, b, n );

    fprintf ( 1, '  %8f  %8f  %6d  %14e  %14e\n', a, b, n, fx, fx2 );

  end

  return
end


text:
function test_suite = test_isLeftOriented
% One-line description here, please.
%   output = test_isLeftOriented(input)
%
%   Example
%   test_isLeftOriented
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2009-04-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function test_HorizLeft(testCase) %#ok<*DEFNU>

line = [10 20 3 0];
p1 = [15 25];

testCase.assertTrue(isLeftOriented(p1, line));

function test_HorizRight(testCase) %#ok<*DEFNU>

line = [10 20 3 0];
p1 = [15 5];

testCase.assertFalse(isLeftOriented(p1, line));

function test_PointArray(testCase) %#ok<*DEFNU>

% one line, three points
line = [10 20 3 0];
pts = [15 25; 15 15; 15 05];

% expect 3-by-1 array
exp = [true; false; false];
testCase.assertEqual(exp, isLeftOriented(pts, line));

function test_LineArray(testCase)

% four lines, one point
lines = [....
    10 10 3 0; ...  % left
    10 10 0 3; ...  % right
    10 10 1.5 1.5; ...  % left
    10 10 1.5 -1.5; ... % left
];
pt = [21.1 32.2];

% expect 1-by-4 array
exp = [true false true true];

testCase.assertEqual(exp, isLeftOriented(pt, lines));

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
classdef ShapeOptimizationSolver < handle

    properties (Access = public)
        tV
        JV
        betaV
        incXvalues
        designVariable
    end

    properties (Access = private)
        topOpt
        cost
        plotter
        TOL
        maxIter
        momentumParameter
        momentumParams
    end

    methods (Access = public)

        function obj = ShapeOptimizationSolver(cParams)
            obj.init(cParams);
        end

        function solve(obj)
            xNew = obj.computeInitialValue();
            xOld = xNew;
            incX = obj.computeIncX(xOld,xNew);
            iter = 1;
            while ~obj.hasConverged(iter,incX)
                [J,dJ] = obj.cost.computeValueAndGradient(xNew);
                beta = obj.computeBeta(iter);
                x = obj.addMomentumTerm(xOld,xNew,beta);
                t = obj.computeLineSearch(x,dJ);
                x = obj.computeGradientStep(x,dJ,t);
                x = obj.computeProjection(x);
                [xOld,xNew] = obj.updateXnewXold(xNew,x);
                incX = obj.computeIncX(xOld,xNew);
                obj.plotCostAndLineSearch(iter,J,t,beta,incX);
                iter = iter + 1;
            end
        end

    end

    methods (Access = private, Static)

        function xNewNew = addMomentumTerm(xOld,xNew,beta)
            xNewNew = xNew + beta*(xNew - xOld);
        end

        function x = computeGradientStep(x,dJ,t)
            x = x - t*dJ;
        end

        function x = computeProjection(x)
            x = max(min(x,1),0);
        end

        function [xOld,xNew] = updateXnewXold(xNew,xNewNew)
            xOld = xNew;
            xNew = xNewNew;
        end

    end

    methods (Access = private)

        function beta = computeBeta(obj,iter)
            s.iter = iter;
            beta = obj.momentumParameter.computeValue(s);
        end

        function init(obj,cParams)
            obj.momentumParams = cParams.momentumParams;
            obj.TOL = cParams.TOL;
            obj.maxIter = cParams.maxIter;
            obj.createSettings();
            obj.createDesignVariable();
            obj.createCost();
            obj.createPlotter();
            obj.createMomentumParameter();
        end

        function createMomentumParameter(obj)
            s = obj.momentumParams;
            obj.momentumParameter = MomentumParameter.create(s);
        end

        function createSettings(obj)
            settings = Settings('Example1');
            translator = SettingsTranslator();
            translator.translate(settings);
            fileName = translator.fileName;
            settingsTopOpt = SettingsTopOptProblem(fileName);
            obj.topOpt = TopOpt_Problem(settingsTopOpt);
        end

        function createDesignVariable(obj)
            obj.designVariable = obj.topOpt.designVariable;
        end

        function createCost(obj)
            s.topOpt = obj.topOpt;
            s.designVariable = obj.designVariable;
            obj.cost = CostComplianceVolume(s);
        end

        function createPlotter(obj)
            s.designVariable = obj.designVariable;
            obj.plotter = PlotterDensity(s);
        end

        function x0 = computeInitialValue(obj)
            x0 = obj.designVariable.value;
        end

        function incX = computeIncX(obj,xOld,xNew)
            incX  = obj.computeNorm(xOld - xNew);
            xNorm = obj.computeNorm(xOld);
            incX = incX/xNorm;
        end

        function itHas = hasConverged(obj,iter,incX)
            if iter == 1
                itHas = false;
            else
                itHas = iter >= obj.maxIter || incX < obj.TOL;
            end
        end

        function t = computeLineSearch(obj,x,dJ)
            tC = 100;
            incT = 1;
            tA = obj.computeAdimensionalLineSearch(x,dJ);
            t = max(tA,incT*tC);
        end

        function t0 = computeAdimensionalLineSearch(obj,x,dJ)
            nX  = obj.computeNorm(x);
            ndJ = obj.computeNorm(dJ);
            t0 = ndJ/nX;
        end

        function plotCostAndLineSearch(obj,iter,J,t,beta,incX)
            obj.JV(iter) = J;
            obj.tV(iter) = t;
            obj.betaV(iter) = beta;
            obj.incXvalues(iter) = incX;
            obj.plotter.plot(obj.JV,obj.tV,obj.betaV,obj.incXvalues);
        end

        function n = computeNorm(obj,x)
            sc = obj.designVariable.scalarProduct;
            s  = sc.computeSP_M(x,x);
            n  = sqrt(s);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function ui = mrViewMeshMenu(ui);
%
% ui = mrViewMeshMenu(ui);
%
% Attaches a mesh menu to a mrViewer UI. Initializes
% it to be hidden, until a segmentation is loaded.
%
% ras 07/05.
if ~exist('ui', 'var') | isempty(ui),  ui = get(gcf, 'UserData'); end

ui.menus.mesh = uimenu(ui.fig, 'Label', 'Mesh');

%%%%%mesh build, display, close options
hFile = uimenu(ui.menus.mesh, 'Label', 'Build \ Open \ Close', ...
                              'Separator', 'off');

cb = 'mrViewLoad(gcf, [], ''Mesh''); mrViewDisplayMesh; ';
ui.menus.meshLoadDisp = uimenu(hFile, 'Label', 'Load and Display', ...
                        'Callback', cb);

ui.menus.meshBuild = uimenu(hFile, 'Label', 'Build Mesh', ...
                        'Separator', 'on', 'Callback', 'mrViewBuildMesh; ');

ui.menus.meshSave = uimenu(hFile, 'Label', 'Save Selected Mesh', ...
                        'Callback', 'mrViewSave(gcf, [], ''Mesh''); ');

ui.menus.meshDisp = uimenu(hFile, 'Label', 'Display Selected Mesh', ...
                        'Callback', 'mrViewDisplayMesh;');

cb = 'mrmSet( mrViewGet(gcf, ''CurMesh''), ''Close'' ); ';
ui.menus.meshCloseDisp = uimenu(hFile, 'Label', 'Close Mesh Display', ...
                        'Callback', cb);

ui.menus.meshClose = uimenu(hFile, 'Label', 'Close Mesh', ...
                        'Callback', 'mrViewCloseMesh(gcf);');

%%%%%project data onto mesh options
ui.menus.meshUpdate = uimenu(ui.menus.mesh, 'Label', 'Project Data Onto Mesh', ...
                            'Separator', 'on', 'Callback', 'mrViewMesh(gcf); ');
ui.menus.meshUpdateAll = uimenu(ui.menus.mesh, 'Label', 'Update All Meshes', ...
                            'Separator', 'off', 'Callback', 'mrViewUpdateAllMeshes; ', ...
                            'Accelerator', '0');

%%%%%mesh smooth, revert options
hInflate = uimenu(ui.menus.mesh, 'Label', 'Inflation', 'Separator', 'on');

cb = ['tmp = mrViewGet(gcf, ''CurMesh''); ' ...
      'tmp = meshSmooth(tmp); ' ...
      'mrViewSet(gcf, ''CurMesh'', tmp); ' ...
      'clear tmp '];
ui.menus.meshSmooth = uimenu(hInflate, 'Label', 'Inflate Mesh', ...
                        'Separator', 'off', 'Callback', cb);

cb = ['tmp = mrViewGet(gcf, ''CurMesh''); ' ...
      'tmp = meshSmooth(tmp, 1); ' ...
      'mrViewSet(gcf, ''CurMesh'', tmp); ' ...
      'clear tmp '];
ui.menus.meshSmooth2 = uimenu(hInflate, 'Label', 'Inflate (Set Params)', ...
                                        'Callback', cb);

cb = ['tmp = mrViewGet(gcf, ''CurMesh''); ' ...
      'tmp = meshSet(tmp, ''Vertices'', meshGet(tmp, ''InitialVertices'')); ' ...
      'tmp = meshSet(tmp, ''Smooth_Relaxation'', 0); ' ...
      'mrmSet(tmp, ''Vertices''); ' ...
      'mrViewSet(gcf, ''CurMesh'', tmp); ' ...
      'clear tmp '];
ui.menus.meshRevert = uimenu(hInflate, 'Label', 'Revert to InitVertices', ...
                                       'Callback', cb);

%%%%%display tools (screenshot, multi-angle, movie)
hDisp = uimenu(ui.menus.mesh, 'Label', 'Display Tools', 'Separator', 'off');

ui.menus.meshPrefs = uimenu(hDisp, 'Label', 'Set Mesh Preferences', ...
                            'Accelerator', '.', 'Callback', 'mrmPreferences');

ui.menus.meshToggleCursor = uimenu(hDisp, 'Label', 'Toggle Cursor', ...
						'Accelerator', '/', ...
						'Callback', 'meshToggleCursor(mrViewGet(gcf, ''mesh''));');

% recolor sulcul/gyral curvature
cb = ['msh = mrViewGet(gcf, ''CurMesh''); ' ...
	 'msh = meshColor(msh, 1); ' ...
	 'mrViewSet(gcf, ''CurMesh'', msh); ' ...
	 'mrViewMesh(gcf); '...
	 'clear msh '];
uimenu(hDisp, 'Label', 'Set Sulcus/Gyrus Modulation', 'Callback', cb);

% Set lighting parameters
uimenu(hDisp, 'Label', 'Set Lighting Levels', 'Callback', 'mrViewMeshLighting;');

% Take mesh snapshot
cb = ['img = mrmGet( mrViewGet(gcf, ''Mesh''), ''screenshot'' ) ./ 255; ' ...
	  'hTmp = figure(''Color'', ''w''); ' ...
	  'imagesc(img); axis image; axis off; ' ...
	  'try, figToCB(hTmp); end; clear hTmp; '...
	  'disp(''Image copied to clipboard as well as variable "img".'')'];
uimenu(hDisp, 'Label', 'Copy Mesh Image to Clipboard', 'Callback', cb);

ui.menus.meshBgColor = uimenu(hDisp, 'Label', 'Set background color...');

% individual callbacks for different colors
colors = {[0 0 0], [.3 .3 .3], [1 1 1], [0 .1 .8]};
names = {'black' 'gray' 'white' 'blue'};
for i = 1:length(colors)
    ui.menus.meshBgColor(i+1) = uimenu(ui.menus.meshBgColor(1), ...
        'Label', names{i}, 'UserData', colors{i}, 'Separator', 'off', ...
        'Callback', 'mrViewSet(gcf, ''MeshBackground'', gcbo); ');
end

cb = 'meshMultiAngle2(mrViewGet(gcf, ''CurMesh''), ''dialog''); ';
ui.menus.meshMultiAngle = uimenu(hDisp, 'Label', 'Multi-Angle Snapshot', ...
                                    'Separator', 'off', 'Callback', cb);

cb = ['VOLUME{1}.mesh{1} = mrViewGet(gcf, ''mesh''); ' ...
	  'meshMovie(VOLUME{1}); ' ...
	  'VOLUME{1}.mesh = {}; '];
ui.menus.meshMovie = uimenu(hDisp, 'Label', 'Mesh Movie', ...
						'Separator', 'off', 'Callback', cb);

%%%%%ROI tools
hRoi = uimenu(ui.menus.mesh, 'Label', 'Mesh ROIs', 'Separator', 'off');

ui.menus.meshCursor = uimenu(hRoi, 'Label', 'Get Cursor Position From Mesh', ...
						'Callback', 'mrViewSetCursorFromMesh(gcf); ');

ui.menus.meshROI1 = uimenu(hRoi, 'Label', 'Get From Mesh (Layer 1)', ...
						'Callback', 'mrViewROI(''mesh'', gcf, ''layer1'');');

ui.menus.meshROI2 = uimenu(hRoi, 'Label', 'Get From Mesh (All Layers)', ...
						'Callback', 'mrViewROI(''mesh'', gcf, ''all'');');

ui.menus.diskROI = uimenu(hRoi, 'Label', 'Create Disk ROI', ...
							'Callback', 'mrViewROI(''disk'', gcf);');

%%%%%select segmentation submenu
% (will be added to in mrViewLoad([], [], 'segmentation'))
ui.menus.meshSeg = uimenu(ui.menus.mesh, 'Label', 'Selected Segmentation', ...
                          'Separator', 'on');

%%%%%select mesh submenu (will be added to in mrViewLoad([], [], 'mesh'))
ui.menus.meshSelect = uimenu(ui.menus.mesh, 'Label', 'Selected Mesh', ...
                             'Separator', 'off');

% initialize to invisible, until a segmentation is loaded
set(ui.menus.mesh, 'Visible', 'off');

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function Xform = inplane2VolXform(rot,trans,scaleFac)
%
% Xform = inplane2VolXform(rot,trans,scaleFac)
%
% Returns 4x4 homogeneous tranform that tranforms from inplane to
% volume.
%
% djh/gmb, '97
%
% Modification:
% - Flip first 2 rows and cols so that it deals with
% (y,x,z) coords instead of (x,y,z).  DJH, 7/98.

A=diag(scaleFac(2,:))*rot*diag(1./scaleFac(1,:));
b = (scaleFac(2,:).*trans)';

Xform = zeros(4,4);
Xform(1:3,1:3)=A;
Xform(1:3,4)=b;
Xform(4,4)=1;

Xform([1 2],:) = Xform([2 1],:);
Xform(:,[1 2]) = Xform(:,[2 1]);


text:
function [model,rxnNames] = addDemandReactionOri(model,metaboliteNameList, addRxnGeneMat)
% addDemandReaction adds demand reactions for a set of metabolites
% The reaction names for the demand reactions will be DM_[metaboliteName]
%
% model = addDemandReaction(model,metaboliteNameList)
%
% INPUTS
% model                 COBRA model structure
% metaboliteNameList    List of metabolite names (cell array)
% addRxnGeneMat         Adds rxnGeneMat to model structure (default = true)
%
% OUTPUTS
% model                 COBRA model structure with added demand reactions
% rxnNames              List of added reactions
%
% Markus Herrgard 5/8/07
% Ines Thiele 03/09 - Corrected reaction coefficient for demand reaction
% Ines Thiele 08/03/2015, made rxnGeneMat optional
if ~exist('addRxnGeneMat','var')
    addRxnGeneMat = 1;
end

if (~iscell(metaboliteNameList))
    tmp = metaboliteNameList;
    clear metaboliteNameList;
    metaboliteNameList{1} = tmp;
end

for i = 1:length(metaboliteNameList)
    rxnName = ['DM_' metaboliteNameList{i}];
    rxnNames{i}=rxnName;
    metaboliteList = {metaboliteNameList{i}};
%     [model,rxnIDexists] = addReaction(model,rxnName,metaboliteList,stoichCoeffList,revFlag,lowerBound,upperBound,objCoeff,subSystem,grRule,geneNameList,systNameList,checkDuplicate, addRxnGeneMat)
    model = addReactionOri(model,rxnName,metaboliteList,-1,false,0,1000,0,'Demand','',[],[],0, addRxnGeneMat);
end

author: opencobra

repo: cobratoolbox


text:
function [knot] = kmph2knot(kmph)
% Convert speed from miles per hour to knots
% Chad A. Greene 2012
knot = kmph*0.5399568034557;

author: Sable

repo: mcbench-benchmarks


text:
function [] = Test_D_Gen_S(block,signal)

% Callback processing function for mask parameter
% - number of <Requirement name> used in the System
set_param(block, 'LinkStatus', 'none');
load_system('Simulink'); % Load Simulink library (background)
load_system('MIL_Test');
if signal < 0
No = str2double(get_param(block,'So')); % Get # of VF-s
else
    No = signal;
end

T1 = strcat(block, '/<Requirement name>/<Preconditions name>');
T2 = strcat(block, '/<Requirement name>');
OldNo = length(find_system(T1, 'regexp', 'on','LookUnderMasks','all','Name','Feature generation'));
OldNoReq = length(find_system(block, 'regexp', 'on','LookUnderMasks','all','Name','<Requirement name>'));
for i = 1:max(No,OldNo)
    if OldNo > No && i > No && i <= OldNo % Remove
        % delete lines

        % block

        SigAtTLines2 = get_param([block '/switch <signal' num2str(i) '>'],'LineHandles');
        delete_line(SigAtTLines2.Outport(1));
        for j = 1:(OldNoReq+1)
         delete_line(SigAtTLines2.Inport(j));
        end
        delete_block([block '/<signal' num2str(i) '>']);
        delete_block([block '/switch <signal' num2str(i) '>']);
        %T2
        for k = 1:OldNoReq
            if k == 1
           T1 = strcat(block, '/<Requirement name>/<Preconditions name>');
           T2 = strcat(block, '/<Requirement name>');
        SigAtTLines1 = get_param([T2 '/switch <signal' num2str(i) '>'],'LineHandles');
        delete_line(SigAtTLines1.Outport(1));
        delete_line(SigAtTLines1.Inport(1));
        delete_line(SigAtTLines1.Inport(2));
        delete_line(SigAtTLines1.Inport(3));

        % T1
        SigAtTLines = get_param([T1 '/Feature generation' num2str(i)],'LineHandles');
        %delete_line(SigAtTLines.Inport(1));
        delete_line(SigAtTLines.Outport(1));
        % delete block
        delete_block([T2 '/<signal' num2str(i) '>']);
        delete_block([T2 '/switch <signal' num2str(i) '>']);
        delete_block([T1 '/Feature generation' num2str(i)]);
        delete_block([T1 '/<signal' num2str(i) '>']);
            else
                 T1 = strcat(block, '/<Requirement name>', num2str(k), '/<Preconditions name>');
                 T2 = strcat(block, '/<Requirement name>', num2str(k));
                 SigAtTLines1 = get_param([T2 '/switch <signal' num2str(i) '>'],'LineHandles');
        delete_line(SigAtTLines1.Outport(1));
        delete_line(SigAtTLines1.Inport(1));
        delete_line(SigAtTLines1.Inport(2));
        delete_line(SigAtTLines1.Inport(3));

        % T1
        SigAtTLines = get_param([T1 '/Feature generation' num2str(i)],'LineHandles');
        %delete_line(SigAtTLines.Inport(1));
        delete_line(SigAtTLines.Outport(1));
        % delete block
        delete_block([T2 '/<signal' num2str(i) '>']);
        delete_block([T2 '/switch <signal' num2str(i) '>']);
        delete_block([T1 '/Feature generation' num2str(i)]);
        delete_block([T1 '/<signal' num2str(i) '>']);
            end
        end

    elseif No > OldNo && i > OldNo && i <= No % Add
         % Block
                add_block('Simulink/Sinks/Out1',[block '/<signal' num2str(i) '>']);
                fix_pos([block '/<signal' num2str(i-1) '>'],[block '/<signal' num2str(i) '>'],0,145);
                add_block([block '/switch <signal1>'],[block '/switch <signal' num2str(i) '>']);
                fix_pos([block '/switch <signal' num2str(i-1) '>'],[block '/switch <signal' num2str(i) '>'],0,60);
                add_line(block,['switch <signal' num2str(i) '>' '/1'],['<signal' num2str(i) '>/' num2str(1)],'autorouting','on');
                add_line(block,['Bus Selector/1'],['switch <signal' num2str(i) '>/' num2str(1)],'autorouting','on');
        % Add blocks
        % <Test data generator>/<Requirement name>/<Preconditions name>
        for j = 1:OldNoReq
           if j == 1
               T1 = strcat(block, '/<Requirement name>/<Preconditions name>');
               T2 = strcat(block, '/<Requirement name>');
               add_block('Simulink/Sinks/Out1',[block '/<Requirement name>/<Preconditions name>/<signal' num2str(i) '>']);
               fix_pos([block '/<Requirement name>/<Preconditions name>/<signal' num2str(i-1) '>'],[block '/<Requirement name>/<Preconditions name>/<signal' num2str(i) '>'],0,-70);
               add_block([block '/<Requirement name>/<Preconditions name>/Feature generation'],[block '/<Requirement name>/<Preconditions name>/Feature generation' num2str(i)]);
               fix_pos([block '/<Requirement name>/<Preconditions name>/Feature generation'],[block '/<Requirement name>/<Preconditions name>/Feature generation' num2str(i)],0,-i*50);
               add_line(T1,['Feature generation' num2str(i) '/1'],['<signal' num2str(i) '>/' num2str(1)],'autorouting','on');
               % T2
               add_block('Simulink/Sinks/Out1',[T2 '/<signal' num2str(i) '>']);
               fix_pos([T2 '/<signal' num2str(i-1) '>'],[T2 '/<signal' num2str(i) '>'],0,145);
               add_block([T2 '/switch <signal1>'],[T2 '/switch <signal' num2str(i) '>']);
               fix_pos([T2 '/switch <signal' num2str(i-1) '>'],[T2 '/switch <signal' num2str(i) '>'],0,60);
               add_line(T2,['switch <signal' num2str(i) '>' '/1'],['<signal' num2str(i) '>/' num2str(1)],'autorouting','on');
               add_line(T2,['Sequencing of test data in time due to Preconditions/1'],['switch <signal' num2str(i) '>/' num2str(1)],'autorouting','on');
               add_line(T2,['Initialisation & Stabilisation/1'],['switch <signal' num2str(i) '>/' num2str(3)],'autorouting','on');
               add_line(T2,['<Preconditions name>/' num2str(i+1)],['switch <signal' num2str(i) '>/' num2str(2)],'autorouting','on');
               % Block
               add_line(block,['<Requirement name>/' num2str(i)],['switch <signal' num2str(i) '>/' num2str(2)],'autorouting','on');
           else

                T1 = strcat(block, '/<Requirement name>', num2str(j), '/<Preconditions name>');
                T2 = strcat(block, '/<Requirement name>', num2str(j));

                add_block('Simulink/Sinks/Out1',[T1 '/<signal' num2str(i) '>']);
                fix_pos([T1 '/<signal' num2str(i-1) '>'],[T1 '/<signal' num2str(i) '>'],0,-70);
                add_block([T1 '/Feature generation'],[T1 '/Feature generation' num2str(i)]);
                fix_pos([T1 '/Feature generation'],[T1 '/Feature generation' num2str(i)],0,-i*50);
                add_line(T1,['Feature generation' num2str(i) '/1'],['<signal' num2str(i) '>/' num2str(1)],'autorouting','on');
                % T2
                add_block('Simulink/Sinks/Out1',[T2 '/<signal' num2str(i) '>']);
                fix_pos([T2 '/<signal' num2str(i-1) '>'],[T2 '/<signal' num2str(i) '>'],0,145);
                add_block([T2 '/switch <signal1>'],[T2 '/switch <signal' num2str(i) '>']);
                fix_pos([T2 '/switch <signal' num2str(i-1) '>'],[T2 '/switch <signal' num2str(i) '>'],0,60);
                add_line(T2,['switch <signal' num2str(i) '>' '/1'],['<signal' num2str(i) '>/' num2str(1)],'autorouting','on');
                add_line(T2,['Sequencing of test data in time due to Preconditions/1'],['switch <signal' num2str(i) '>/' num2str(1)],'autorouting','on');
                add_line(T2,['Initialisation & Stabilisation/1'],['switch <signal' num2str(i) '>/' num2str(3)],'autorouting','on');
                add_line(T2,['<Preconditions name>/' num2str(i+1)],['switch <signal' num2str(i) '>/' num2str(2)],'autorouting','on');
                add_line(block,['<Requirement name>' num2str(j) '/' num2str(i)],['switch <signal' num2str(i) '>/' num2str(j+1)],'autorouting','on');
           end
        end

    end
end


text:
function [newgt, newdet] = dropObjectsInIgr(gt, det, imgHeight, imgWidth)
%% drop objects in ignored region
newgt = [];
newdet = [];
numFr = max(gt(:,1));
for fr = 1:numFr
    % parse objects
    idxGtFr = gt(:,1) == fr & gt(:, 8) ~= 0;
    curgt = gt(idxGtFr,:);
    idxDetFr = det(:,1) == fr;
    curdet = det(idxDetFr,:);
    % parse ignored regions
    idxIgr = gt(:,1) == fr & gt(:, 8) == 0;
    igrRegion = gt(idxIgr, 3:6);
    if(~isempty(igrRegion))
        igrMap = zeros(imgHeight, imgWidth);
        numIgr = size(igrRegion,1);
        for j = 1:numIgr
            igrMap(igrRegion(j,2):min(imgHeight,igrRegion(j,2)+igrRegion(j,4)),igrRegion(j,1):min(imgWidth,igrRegion(j,1)+igrRegion(j,3))) = 1;
        end
        intIgrMap = createIntImg(double(igrMap));
        idxLeft = [];
        for i = 1:size(curgt, 1)
            pos = max(1,round(curgt(i,3:6)));
            x = max(1, min(imgWidth, pos(1)));
            y = max(1, min(imgHeight, pos(2)));
            w = pos(3);
            h = pos(4);
            tl = intIgrMap(y, x);
            tr = intIgrMap(y, min(imgWidth,x+w));
            bl = intIgrMap(max(1,min(imgHeight,y+h)), x);
            br = intIgrMap(max(1,min(imgHeight,y+h)), min(imgWidth,x+w));
            igrVal = tl + br - tr - bl;
            if(igrVal/(h*w)<0.5)
                idxLeft = cat(1, idxLeft, i);
            end
        end

        curdet = curdet(idxLeft, :);
        idxLeft = [];
        for i = 1:size(curdet, 1)
            pos = max(1,round(curdet(i,3:6)));
            x = max(1, min(imgWidth, pos(1)));
            y = max(1, min(imgHeight, pos(2)));
            w = pos(3);
            h = pos(4);
            tl = intIgrMap(y, x);
            tr = intIgrMap(y, min(imgWidth,x+w));
            bl = intIgrMap(max(1,min(imgHeight,y+h)), x);
            br = intIgrMap(max(1,min(imgHeight,y+h)), min(imgWidth,x+w));
            igrVal = tl + br - tr - bl;
            if(igrVal/(h*w)<0.5)
                idxLeft = cat(1, idxLeft, i);
            end
        end
        curdet = curdet(idxLeft, :);
    end
    newgt = cat(1, newgt, curgt);
    newdet = cat(1, newdet, curdet);

end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function net = resnet52_new_hope()

%----------------------img cnn----------------------
netStruct = load('./data/imagenet-vgg-verydeep-19.mat') ;
net = dagnn.DagNN.loadobj(netStruct) ;
net.removeLayer('fc8');
net.removeLayer('relu7');
net.removeLayer('prob');

for i = 1:numel(net.params)
    if(mod(i,2)==0)
        net.params(i).learningRate= 0; %0.02;
    else net.params(i).learningRate= 0; %0.001;
    end
    net.params(i).weightDecay=0; %1;
end

net.params(1).learningRate = 0; %1e-5;

fc1Block = dagnn.Conv('size',[1 1 4096 4096],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc1',fc1Block,{'x40'},{'fc1'},{'fc1f'});
net.addLayer('fc1bn',dagnn.BatchNorm(),{'fc1'},{'fc1bn'},...
    {'fc1bn_w','fc1bn_b','fc1bn_m'});
net.addLayer('fc1x',dagnn.ReLU(),{'fc1bn'},{'fc1bnx'});
fc1Block = dagnn.Conv('size',[1 1 4096 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc1_1',fc1Block,{'fc1bnx'},{'fc1_1'},{'fc1_1f'});
net.addLayer('fc1_1bn',dagnn.BatchNorm(),{'fc1_1'},{'fc1_1bn'},...
    {'fc1_1bn_w','fc1_1bn_b','fc1_1bn_m'});
net.addLayer('fc1_1x',dagnn.ReLU(),{'fc1_1bn'},{'fc1_1bnx'});
net.addLayer('dropout',dagnn.DropOut('rate',0.5),{'fc1_1bnx'},{'fc1_1bnxd'});

%----------------------char cnn----------------------
% input is 1*32*20074*16
fc2Block = dagnn.Conv('size',[1 1 29972 300],'hasBias',true,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2',fc2Block,{'data2'},{'fc2'},{'fc2f','fc2b'});
net.addLayer('fc2bn',dagnn.BatchNorm(),{'fc2'},{'fc2bn'},...
    {'fc2bn_w','fc2bn_b','fc2bn_m'});
%net.addLayer('fc2x',dagnn.ReLU(),{'fc2bn'},{'fc2bnx'});
%net.addLayer('dropout_diction',dagnn.DropOut('rate',0.5),{'fc2bn'},{'fc2bnd'});
% 32*256
convBlock = dagnn.Conv('size',[1 1 300 128],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_1_1',convBlock,{'fc2bn'},{'fc2_1_1'},{'fc2_1_1f'});
net.addLayer('fc2_1_1bn',dagnn.BatchNorm(),{'fc2_1_1'},{'fc2_1_1bn'},...
    {'fc2_1_1bn_w','fc2_1_1bn_b','fc2_1_1bn_m'});
net.addLayer('fc2_1_1x',dagnn.ReLU(),{'fc2_1_1bn'},{'fc2_1_1bnx'});

convBlock = dagnn.Conv('size',[1 2 128 128],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
net.addLayer('fc2_1_2',convBlock,{'fc2_1_1bnx'},{'fc2_1_2'},{'fc2_1_2f'});
net.addLayer('fc2_1_2bn',dagnn.BatchNorm(),{'fc2_1_2'},{'fc2_1_2bn'},...
    {'fc2_1_2bn_w','fc2_1_2bn_b','fc2_1_2bn_m'});
net.addLayer('fc2_1_2x',dagnn.ReLU(),{'fc2_1_2bn'},{'fc2_1_2bnx'});

convBlock = dagnn.Conv('size',[1 1 128 256],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_1_3',convBlock,{'fc2_1_2bnx'},{'fc2_1_3'},{'fc2_1_3f'});
net.addLayer('fc2_1_3bn',dagnn.BatchNorm(),{'fc2_1_3'},{'fc2_1_3bn'},...
    {'fc2_1_3bn_w','fc2_1_3bn_b','fc2_1_3bn_m'});

convBlock = dagnn.Conv('size',[1 1 300 256],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_1b',convBlock,{'fc2bn'},{'fc2_1b'},{'fc2_1bf'});
net.addLayer('fc2_1bbn',dagnn.BatchNorm(),{'fc2_1b'},{'fc2_1bbn'},...
    {'fc2_1bbn_w','fc2_1bbn_b','fc2_1bbn_m'});

net.addLayer('fc2_1sum',dagnn.Sum(),{'fc2_1_3bn','fc2_1bbn'},...
    {'fc2_1sum'});
net.addLayer('fc2_1x',dagnn.ReLU(),{'fc2_1sum'},{'fc2_1sumx'});

% 32*256
for i  = 2:3
    convBlock = dagnn.Conv('size',[1 1 256 64],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc2_%d_1',i),convBlock,{sprintf('fc2_%dsumx',i-1)},{sprintf('fc2_%d_1',i)}, ...
        {sprintf('fc2_%d_1f',i)});
    net.addLayer(sprintf('fc2_%d_1bn_1',i),dagnn.BatchNorm(),{sprintf('fc2_%d_1',i)},{sprintf('fc2_%d_1bn',i)},...
        {sprintf('fc2_%d_1bn_w',i),sprintf('fc2_%d_1bn_b',i),sprintf('fc2_%d_1bn_m',i)});
    net.addLayer(sprintf('fc2_%d_1x',i),dagnn.ReLU(),{sprintf('fc2_%d_1bn',i)},{sprintf('fc2_%d_1bnx',i)});

    convBlock = dagnn.Conv('size',[1 2 64 64],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
    net.addLayer( sprintf('fc2_%d_2',i),convBlock,{sprintf('fc2_%d_1bnx',i)},{sprintf('fc2_%d_2',i)}, ...
        {sprintf('fc2_%d_2f',i)});
    net.addLayer(sprintf('fc2_%d_2bn',i),dagnn.BatchNorm(),{sprintf('fc2_%d_2',i)},{sprintf('fc2_%d_2bn',i)},...
        {sprintf('fc2_%d_2bn_w',i),sprintf('fc2_%d_2bn_b',i),sprintf('fc2_%d_2bn_m',i)});
    net.addLayer(sprintf('fc2_%d_2x',i),dagnn.ReLU(),{sprintf('fc2_%d_2bn',i)},{sprintf('fc2_%d_2bnx',i)});

    convBlock = dagnn.Conv('size',[1 1 64 256],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc2_%d_3',i),convBlock,{sprintf('fc2_%d_2bnx',i)},{sprintf('fc2_%d_3',i)}, ...
        {sprintf('fc2_%d_3f',i)});
    net.addLayer(sprintf('fc2_%d_3bn',i),dagnn.BatchNorm(),{sprintf('fc2_%d_3',i)},{sprintf('fc2_%d_3bn',i)},...
        {sprintf('fc2_%d_3bn_w',i),sprintf('fc2_%d_3bn_b',i),sprintf('fc2_%d_3bn_m',i)});

    net.addLayer(sprintf('fc2_%dsum',i),dagnn.Sum(),{sprintf('fc2_%dsumx',i-1),sprintf('fc2_%d_3bn',i)},...
        {sprintf('fc2_%dsum',i)});
    net.addLayer(sprintf('fc2_%dx',i),dagnn.ReLU(),{sprintf('fc2_%dsum',i)},{sprintf('fc2_%dsumx',i)});
end

%32*256
convBlock = dagnn.Conv('size',[1 1 256 512],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_4a_1',convBlock,{'fc2_3sumx'},{'fc2_4a_1'},{'fc2_4a_1f'});
net.addLayer('fc2_4a_1bn',dagnn.BatchNorm(),{'fc2_4a_1'},{'fc2_4a_1bn'},...
    {'fc2_4a_1bn_w','fc2_4a_1bn_b','fc2_4a_1bn_m'});
net.addLayer('fc2_4a_1x',dagnn.ReLU(),{'fc2_4a_1bn'},{'fc2_4a_1bnx'});

convBlock = dagnn.Conv('size',[1 2 512 512],'hasBias',false,'stride',[2,2],'pad',[0,0,1,0]);
net.addLayer('fc2_4a_2',convBlock,{'fc2_4a_1bnx'},{'fc2_4a_2'},{'fc2_4a_2f'});
net.addLayer('fc2_4a_2bn',dagnn.BatchNorm(),{'fc2_4a_2'},{'fc2_4a_2bn'},...
    {'fc2_4a_2bn_w','fc2_4a_2bn_b','fc2_4a_2bn_m'});
net.addLayer('fc2_4a_2x',dagnn.ReLU(),{'fc2_4a_2bn'},{'fc2_4a_2bnx'});

convBlock = dagnn.Conv('size',[1 1 512 512],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_4a_3',convBlock,{'fc2_4a_2bnx'},{'fc2_4a_3'},{'fc2_4a_3f'});
net.addLayer('fc2_4a_3bn',dagnn.BatchNorm(),{'fc2_4a_3'},{'fc2_4a_3bn'},...
    {'fc2_4a_3bn_w','fc2_4a_3bn_b','fc2_4a_3bn_m'});

convBlock = dagnn.Conv('size',[1 1 256 512],'hasBias',false,'stride',[2,2],'pad',[0,0,0,0]);
net.addLayer('fc2_4b',convBlock,{'fc2_3sumx'},{'fc2_4b'},{'fc2_4bf'});
net.addLayer('fc2_4bbn',dagnn.BatchNorm(),{'fc2_4b'},{'fc2_4bbn'},...
    {'fc2_4bbn_w','fc2_4bbn_b','fc2_4bbn_m'});

%16*512
net.addLayer('fc2_4sum',dagnn.Sum(),{'fc2_4a_3bn','fc2_4bbn'},...
    {'fc2_4sum'});
net.addLayer('fc2_4x',dagnn.ReLU(),{'fc2_4sum'},{'fc3_1sumx'});

for i  = 2:4
    convBlock = dagnn.Conv('size',[1 1 512 128],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc3_%d_1',i),convBlock,{sprintf('fc3_%dsumx',i-1)},{sprintf('fc3_%d_1',i)}, ...
        {sprintf('fc3_%d_1f',i)});
    net.addLayer(sprintf('fc3_%d_1bn_1',i),dagnn.BatchNorm(),{sprintf('fc3_%d_1',i)},{sprintf('fc3_%d_1bn',i)},...
        {sprintf('fc3_%d_1bn_w',i),sprintf('fc3_%d_1bn_b',i),sprintf('fc3_%d_1bn_m',i)});
    net.addLayer(sprintf('fc3_%d_1x',i),dagnn.ReLU(),{sprintf('fc3_%d_1bn',i)},{sprintf('fc3_%d_1bnx',i)});

    convBlock = dagnn.Conv('size',[1 2 128 128],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
    net.addLayer( sprintf('fc3_%d_2',i),convBlock,{sprintf('fc3_%d_1bnx',i)},{sprintf('fc3_%d_2',i)}, ...
        {sprintf('fc3_%d_2f',i)});
    net.addLayer(sprintf('fc3_%d_2bn',i),dagnn.BatchNorm(),{sprintf('fc3_%d_2',i)},{sprintf('fc3_%d_2bn',i)},...
        {sprintf('fc3_%d_2bn_w',i),sprintf('fc3_%d_2bn_b',i),sprintf('fc3_%d_2bn_m',i)});
    net.addLayer(sprintf('fc3_%d_2x',i),dagnn.ReLU(),{sprintf('fc3_%d_2bn',i)},{sprintf('fc3_%d_2bnx',i)});

    convBlock = dagnn.Conv('size',[1 1 128 512],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc3_%d_3',i),convBlock,{sprintf('fc3_%d_2bnx',i)},{sprintf('fc3_%d_3',i)}, ...
        {sprintf('fc3_%d_3f',i)});
    net.addLayer(sprintf('fc3_%d_3bn',i),dagnn.BatchNorm(),{sprintf('fc3_%d_3',i)},{sprintf('fc3_%d_3bn',i)},...
        {sprintf('fc3_%d_3bn_w',i),sprintf('fc3_%d_3bn_b',i),sprintf('fc3_%d_3bn_m',i)});

    net.addLayer(sprintf('fc3_%dsum',i),dagnn.Sum(),{sprintf('fc3_%dsumx',i-1),sprintf('fc3_%d_3bn',i)},...
        {sprintf('fc3_%dsum',i)});
    net.addLayer(sprintf('fc3_%dx',i),dagnn.ReLU(),{sprintf('fc3_%dsum',i)},{sprintf('fc3_%dsumx',i)});
end

convBlock = dagnn.Conv('size',[1 1 512 1024],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc3_5a_1',convBlock,{'fc3_4sumx'},{'fc3_5a_1'},{'fc3_5a_1f'});
net.addLayer('fc3_5a_1bn',dagnn.BatchNorm(),{'fc3_5a_1'},{'fc3_5a_1bn'},...
    {'fc3_5a_1bn_w','fc3_5a_1bn_b','fc3_5a_1bn_m'});
net.addLayer('fc3_5a_1x',dagnn.ReLU(),{'fc3_5a_1bn'},{'fc3_5a_1bnx'});

convBlock = dagnn.Conv('size',[1 2 1024 1024],'hasBias',false,'stride',[2,2],'pad',[0,0,1,0]);
net.addLayer('fc3_5a_2',convBlock,{'fc3_5a_1bnx'},{'fc3_5a_2'},{'fc3_5a_2f'});
net.addLayer('fc3_5a_2bn',dagnn.BatchNorm(),{'fc3_5a_2'},{'fc3_5a_2bn'},...
    {'fc3_5a_2bn_w','fc3_5a_2bn_b','fc3_5a_2bn_m'});
net.addLayer('fc3_5a_2x',dagnn.ReLU(),{'fc3_5a_2bn'},{'fc3_5a_2bnx'});

convBlock = dagnn.Conv('size',[1 1 1024 1024],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc3_5a_3',convBlock,{'fc3_5a_2bnx'},{'fc3_5a_3'},{'fc3_5a_3f'});
net.addLayer('fc3_5a_3bn',dagnn.BatchNorm(),{'fc3_5a_3'},{'fc3_5a_3bn'},...
    {'fc3_5a_3bn_w','fc3_5a_3bn_b','fc3_5a_3bn_m'});

convBlock = dagnn.Conv('size',[1 1 512 1024],'hasBias',false,'stride',[2,2],'pad',[0,0,0,0]);
net.addLayer('fc3_5b',convBlock,{'fc3_4sumx'},{'fc3_5b'},{'fc3_5bf'});
net.addLayer('fc3_5bbn',dagnn.BatchNorm(),{'fc3_5b'},{'fc3_5bbn'},...
    {'fc3_5bbn_w','fc3_5bbn_b','fc3_5bbn_m'});

%8*1024
net.addLayer('fc3_5sum',dagnn.Sum(),{'fc3_5a_3bn','fc3_5bbn'},...
    {'fc3_5sum'});
net.addLayer('fc3_5x',dagnn.ReLU(),{'fc3_5sum'},{'fc4_1sumx'});

for i  = 2:6
    convBlock = dagnn.Conv('size',[1 1 1024 256],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc4_%d_1',i),convBlock,{sprintf('fc4_%dsumx',i-1)},{sprintf('fc4_%d_1',i)}, ...
        {sprintf('fc4_%d_1f',i)});
    net.addLayer(sprintf('fc4_%d_1bn_1',i),dagnn.BatchNorm(),{sprintf('fc4_%d_1',i)},{sprintf('fc4_%d_1bn',i)},...
        {sprintf('fc4_%d_1bn_w',i),sprintf('fc4_%d_1bn_b',i),sprintf('fc4_%d_1bn_m',i)});
    net.addLayer(sprintf('fc4_%d_1x',i),dagnn.ReLU(),{sprintf('fc4_%d_1bn',i)},{sprintf('fc4_%d_1bnx',i)});

    convBlock = dagnn.Conv('size',[1 2 256 256],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
    net.addLayer( sprintf('fc4_%d_2',i),convBlock,{sprintf('fc4_%d_1bnx',i)},{sprintf('fc4_%d_2',i)}, ...
        {sprintf('fc4_%d_2f',i)});
    net.addLayer(sprintf('fc4_%d_2bn',i),dagnn.BatchNorm(),{sprintf('fc4_%d_2',i)},{sprintf('fc4_%d_2bn',i)},...
        {sprintf('fc4_%d_2bn_w',i),sprintf('fc4_%d_2bn_b',i),sprintf('fc4_%d_2bn_m',i)});
    net.addLayer(sprintf('fc4_%d_2x',i),dagnn.ReLU(),{sprintf('fc4_%d_2bn',i)},{sprintf('fc4_%d_2bnx',i)});

    convBlock = dagnn.Conv('size',[1 1 256 1024],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc4_%d_3',i),convBlock,{sprintf('fc4_%d_2bnx',i)},{sprintf('fc4_%d_3',i)}, ...
        {sprintf('fc4_%d_3f',i)});
    net.addLayer(sprintf('fc4_%d_3bn',i),dagnn.BatchNorm(),{sprintf('fc4_%d_3',i)},{sprintf('fc4_%d_3bn',i)},...
        {sprintf('fc4_%d_3bn_w',i),sprintf('fc4_%d_3bn_b',i),sprintf('fc4_%d_3bn_m',i)});

    net.addLayer(sprintf('fc4_%dsum',i),dagnn.Sum(),{sprintf('fc4_%dsumx',i-1),sprintf('fc4_%d_3bn',i)},...
        {sprintf('fc4_%dsum',i)});
    net.addLayer(sprintf('fc4_%dx',i),dagnn.ReLU(),{sprintf('fc4_%dsum',i)},{sprintf('fc4_%dsumx',i)});
end

convBlock = dagnn.Conv('size',[1 1 1024 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc4_7a_1',convBlock,{'fc4_6sumx'},{'fc4_7a_1'},{'fc4_7a_1f'});
net.addLayer('fc4_7a_1bn',dagnn.BatchNorm(),{'fc4_7a_1'},{'fc4_7a_1bn'},...
    {'fc4_7a_1bn_w','fc4_7a_1bn_b','fc4_7a_1bn_m'});
net.addLayer('fc4_7a_1x',dagnn.ReLU(),{'fc4_7a_1bn'},{'fc4_7a_1bnx'});

convBlock = dagnn.Conv('size',[1 2 2048 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
net.addLayer('fc4_7a_2',convBlock,{'fc4_7a_1bnx'},{'fc4_7a_2'},{'fc4_7a_2f'});
net.addLayer('fc4_7a_2bn',dagnn.BatchNorm(),{'fc4_7a_2'},{'fc4_7a_2bn'},...
    {'fc4_7a_2bn_w','fc4_7a_2bn_b','fc4_7a_2bn_m'});
net.addLayer('fc4_7a_2x',dagnn.ReLU(),{'fc4_7a_2bn'},{'fc4_7a_2bnx'});

convBlock = dagnn.Conv('size',[1 1 2048 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc4_7a_3',convBlock,{'fc4_7a_2bnx'},{'fc4_7a_3'},{'fc4_7a_3f'});
net.addLayer('fc4_7a_3bn',dagnn.BatchNorm(),{'fc4_7a_3'},{'fc4_7a_3bn'},...
    {'fc4_7a_3bn_w','fc4_7a_3bn_b','fc4_7a_3bn_m'});

convBlock = dagnn.Conv('size',[1 1 1024 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc4_7b',convBlock,{'fc4_6sumx'},{'fc4_7b'},{'fc4_7bf'});
net.addLayer('fc4_7bbn',dagnn.BatchNorm(),{'fc4_7b'},{'fc4_7bbn'},...
    {'fc4_7bbn_w','fc4_7bbn_b','fc4_7bbn_m'});

%8*2048
net.addLayer('fc4_7sum',dagnn.Sum(),{'fc4_7a_3bn','fc4_7bbn'},...
    {'fc4_7sum'});
net.addLayer('fc4_7x',dagnn.ReLU(),{'fc4_7sum'},{'fc5_1sumx'});

poolBlock = dagnn.Pooling('poolSize',[1 8]);
net.addLayer('fc5_1',poolBlock,{'fc5_1sumx'},{'fc5_1bnx'});
%net.addLayer('fc5_1bn',dagnn.BatchNorm(),{'fc5_1'},{'fc5_1bn'},...
 %   {'fc5_1bn_w','fc5_1bn_b','fc5_1bn_m'});
%net.addLayer('fc5_1x',dagnn.ReLU(),{'fc5_1bn'},{'fc5_1bnx'});

fc5_2Block = dagnn.Conv('size',[1 1 2048 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc5_2',fc5_2Block,{'fc5_1bnx'},{'fc5_2'},{'fc5_2f'});
net.addLayer('fc5_2bn',dagnn.BatchNorm(),{'fc5_2'},{'fc5_2bn'},...
    {'fc5_2bn_w','fc5_2bn_b','fc5_2bn_m'});
net.addLayer('fc5_2x',dagnn.ReLU(),{'fc5_2bn'},{'fc5_2bnx'});
net.addLayer('dropout2',dagnn.DropOut('rate',0.5),{'fc5_2bnx'},{'fc5_2bnxd'});

%----------------------add share layer----------------------
%1

fc_imgBlock = dagnn.Conv('size',[1 1 2048 113287],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc_img',fc_imgBlock,{'fc1_1bnxd'},{'prediction_img'},{'fcsharef'});
net.addLayer('softmaxloss_img',dagnn.Loss('loss','softmaxlog'),{'prediction_img','label_img'},'objective_img');
net.addLayer('top1err_img', dagnn.Loss('loss', 'classerror'), ...
    {'prediction_img','label_img'}, 'top1err_img') ;
net.addLayer('top5err_img', dagnn.Loss('loss', 'topkerror', ...
    'opts', {'topK',5}), ...
    {'prediction_img','label_img'}, 'top5err_img') ;
%2
fc_txtBlock = dagnn.Conv('size',[1 1 2048 113287],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc_txt',fc_txtBlock,{'fc5_2bnxd'},{'prediction_txt'},{'fcsharef'});
net.addLayer('softmaxloss_txt',dagnn.Loss('loss','softmaxlog'),{'prediction_txt','label_txt'},'objective_txt');
net.addLayer('top1err_txt', dagnn.Loss('loss', 'classerror'), ...
    {'prediction_txt','label_txt'}, 'top1err_txt') ;
net.addLayer('top5err_txt', dagnn.Loss('loss', 'topkerror', ...
    'opts', {'topK',5}), ...
    {'prediction_txt','label_txt'}, 'top5err_txt') ;
%}

%}
net.initParams();
%----------NOTICE--------------
first = net.getParamIndex('fc2f');

net.params(first).learningRate = 1e-3;  %w
net.params(first+1).learningRate = 1e-3;  %b

load('./dataset/MSCOCO-prepare/COCO_dictionary.mat');
net.params(first).value = reshape(single(subset.features'),1,1,29972,300);

end

author: layumi

repo: Image-Text-Embedding


text:
function [nodes, edges] = Connections2Graph(uniqueSparse,vertices);
%
%  [nodes, edges] = Connections2Graph(uniqueSparse,vertices);
%
% AUTHOR:  MK/BW
% DATE:    09.06.99
% PURPOSE:
%
%   Compute nodes and edges for gray graph data.  The graph is
% computed from the uniqueSparse matrix created in
% Mesh2GrayGraph, along with the mesh vertices (floats).
% In this routine, the vertices are rounded to be integer
% valued "nodes."
%
% TODO:
% 1.  Figure out which routines break if nodes are not integers.
%
% Vertices must be integers in mrGray and other places.  But,
% the mesh falls on float values.  When we round, some of the
% mesh vertices that are close, fall on the same integer
% node.  This may not happen if we use cube faces.
%
% Possibly, we should scale up the precision by 10 or so.
% Or, possibly, we can keep the precision. I am not sure what
% consequences this will have for mapping data from gLocs3d
% ...  Needs to be thought through.  Who breaks if nodes are
% not integers?
%
% BW 09.05.99
%

% Convert the mrGray format vertices to mrLoadRet format
% Perhaps the rounding has already happened, but that shouldn't
% matter.  If it is arounded by now, then this won't hurt.
%
vertices = mrGray2mrLoadRet(round(vertices));

% The rounding now collapses some of the nodes to the same
% position.  What should we do to remove the case of multiple
% nodes at the same position?
%
% One possibility is to keep them all, say by scaling the
% positions by 10 and holding things as integers.  Another
% possibility is to collapse the vertices once more and combine
% their connections
%
disp('Connections2Graph isn't right yet.');

% NEW numNodes
%
numNodes = size(uniqueSparse,1);

nodes = zeros(8,numNodes);
nodes(1,:) = vertices(:,1)';
nodes(2,:) = vertices(:,2)';
nodes(3,:) = vertices(:,3)';
nodes(4,:) = full(sum(uniqueSparse,2))';
nodes(5,:) = [0 cumsum(nodes(4,1:(numNodes-1)))] +1;
nodes(6,:) = ones(1,numNodes);
%% 7,8 are junk.  No need to fill them.

% The edges are defined with Matlab uses numbering from 1 to N.
%
% This code performs the conversion for all nodes all at once.
% We use the mod operation to find the column number.
% We first subtracting 1 and then mod by numNodes.  This produces
% a column number that runs between ([0,numNodes-1]).  We add 1
% to get back to Matlab numbering.
%
edges = mod((find(uniqueSparse)-1),numNodes)+1;

% keepNodes, and other mrGray functions want the edges to be a
% row, not a column.
%
edges = edges(:)';

return;


text:
function msm_to_mm_array_real_skewsymmetric ( output_filename, a )

%*****************************************************************************80
%
%% MSM_TO_MM_ARRAY_REAL_SKEWSYMMETRIC writes a "matrix array real skew-symmetric" Matrix Market file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the name of the file to which the information
%    should be written.
%
%    Input, sparse matrix A, the NROW by NCOL matrix, stored in MATLAB sparse
%    matrix format, which is to be written to the file.
%
  [ nrow, ncol ] = size ( a );

  fid = fopen ( output_filename, 'wt+' );

  if ( fid < 0 );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MSM_TO_MM_ARRAY_REAL_SKEWSYMMETRIC - Fatal error!\n' );
    fprintf ( 1, '  Cannot open the output file.\n' );
    error ( 'MSM_TO_MM_ARRAY_REAL_SKEWSYMMETRIC - Fatal error!' );
  end

  fprintf ( fid, '%%%%MatrixMarket matrix array real skew-symmetric\n');
  fprintf ( fid, '%%%%  Created by MSM_TO_MM_ARRAY_REAL_SKEWSYMMETRIC.M\n' );
  fprintf ( fid, '  %d  %d\n', nrow, ncol );

  for j = 1 : ncol
    for i = j + 1 : nrow
      fprintf ( fid, '  %f\n', a(i,j) );
    end
  end

  fclose ( fid );

  return
end


text:
function yi = nearest_interp_1d ( nd, xd, yd, ni, xi )

%*****************************************************************************80
%
%% NEAREST_INTERP_1D evaluates the nearest neighbor interpolant.
%
%  Discussion:
%
%    The nearest neighbor interpolant L(ND,XD,YD)(X) is the piecewise
%    constant function which interpolates the data (XD(I),YD(I)) for I = 1
%    to ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ND, the number of data points.
%    ND must be at least 1.
%
%    Input, real XD(ND,1), the data points.
%
%    Input, real YD(ND,1), the data values.
%
%    Input, integer NI, the number of interpolation points.
%
%    Input, real XI(NI,1), the interpolation points.
%
%    Output, real YI(NI,1), the interpolated values.
%

%
%  KNNSEARCH is a built-in Matlab function which, with this call, returns
%  the indices in XD of the nearest elements to the elements in XI.
%
  idx = knnsearch ( xd, xi );

  yi = yd(idx);

  return
end


text:
function [corr,t,p, fdrsig, fdrthresh] = correlation(meth,x,varargin)
% Multiple types of correlations, including Spearman's rho
% (nonparametric) and phi (dichotomous)
%
% :Usage:
% ::
%
%     [corr,t,p,fdrp, fdrthresh] = correlation(method,x,[y],['matrix'])
%
%     IN PROGRESS : Warning : Use at your own risk.
%     Some methods are not adequately tested yet.
%     Spearman's rho does not correct for ties
%
% :Inputs:
%
%   **Methods:**
%        String indicating the method for computing the correlation
%        coefficient
%        - Pearson's r.          Enter: {'r','pearson',[]}
%        - IRLS                  Enter: {'irls','robust'}
%        - Phi                   Enter: {'phi'}
%        - Spearman's rho        Enter: {'rho','spearman'}
%        - Kendall's Tau (a)     Enter: {'taua','kendalla'}
%        - Tau (b)               Enter: {'tau','kendall','taub','kendallb'}
%        - Gamma                 Enter: {'gamma','kruskal'}
%
%   **x:**
%        Matrix of observations (n instances by p varianbles)
%
% :Optional Inputs:
%
%   **varargin:**
%        To be documented
%
% :Outputs:
%
%   **OUT:**
%        Output stats structure
%
% :Examples:
% ::
%
%    % Corelation between two variables, Pearson's
%    x = rand(10,1); y = rand(10,1);
%    [corr,t,p] = correlation('r',x,y);
%
%    % Correlation matrix of 10 variables, phi correlation:
%    studybyroi = magic(10);
%    [corr,t,p] = correlation('phi',studybyroi);
%
% :See Also: correlation_fast_series.m
%
% ..
%    tor wager, november 2006, Jan 2007
% ..

corr = []; t = []; p = [];

% get y data, if entered
% --------------------------------
for i = 1:length(varargin)
    if ~ischar(varargin{i}) && ~any(size(varargin{i}) - size(x))
        y = varargin{i};
    end
end

% get flag for whether to compute correls among all possible pairs
% --------------------------------
matrixFormFlag = 0;
if any(strcmp(varargin,'matrix')) || ~exist('y','var')
    matrixFormFlag = 1;
end

[n,nvars] = size(x);

% Get indices
% --------------------------------
if matrixFormFlag
    % input is matrix
    [rows,cols,npairs] = corrcoef_indices(nvars);

else
    % input is columns of consecutive pairs to correlate
    npairs = size(x,2); % number of pairs of correlations to compute
end

doverbose = 0;
if npairs > 1000, doverbose = 1; fprintf('%03d%%', 0); end

% Compute correlations
% --------------------------------
for i = 1:npairs

    if doverbose, fprintf('\b\b\b\b%03d%%', round(i*100 ./ npairs)), end

    % get data
    if matrixFormFlag
        x1 = x(:,rows(i)); x2 = x(:,cols(i));
    else
        x1 = x(:,i); x2 = y(:,i);
    end

    if nargout > 1
        [c,tt,pp] = compute_single_correl(meth,x1,x2,n);
    else
        c = compute_single_correl(meth,x1,x2,n);
    end

    corr(i,1) = c;

    if nargout > 1
        if ~isempty(tt), t(i,1) = tt; end
        if ~isempty(pp), p(i,1) = pp; end
    end
end

% Reconstruct into matrix, if needed
% --------------------------------
if matrixFormFlag
    corr = reconstruct(corr,nvars,npairs,rows,cols);

    corr = corr + eye(nvars);

    if ~isempty(t), t = reconstruct(t,nvars,npairs,rows,cols); end

    if ~isempty(p), p = reconstruct(p,nvars,npairs,rows,cols); end

    % FDR correction, if requested
    % only works for matrices, otherwise error
    if nargout > 3
        [fdrthresh,fdrsig] = fdr_correct_pvals(p,corr);
    end

elseif nargout > 3
    error('FDR output is only allowed for matrix form output.  Try requesting fewer outputs')
end

return

% sub-functions

function [est,t,p] = compute_single_correl(meth,x,y,n)

est = [];
t = [];
p = [];

[wasnan, x, y] = nanremove(x, y);

switch lower(meth)

    case {'irls','robust'}
        [b,stats]=robustfit(x,y,'bisquare'); %,[],'off');
        est = weighted_corrcoef([x y], stats.w);
        est = est(1,2);

        t = stats.t(2);
        p = stats.p(2);

    case {'r','pearson',[]}

        if islogical(x)
            warning('Logical vectors not appropriate for Pearson''s r.');
            x = double(x);
            y = double(y);
        end

        [est,p] = corrcoef([x y]);
        est = est(1,2);
        p = p(1,2);

    case 'phi'
        % Phi is for dichotomous (2-level) variables
        % formula from Robert Yaffee, NYU, online.

        tab = make_crosstabs(x,y);
        num = det(tab);                 % same as: num = tab(1,1)*tab(2,2) - tab(1,2)*tab(2,1);
        den = ( prod(sum(tab,2)) .* prod(sum(tab,1)) ).^.5;
        est = num ./ den;

        % F = varexp/dfexp / varunexp/dferror; t = sqrt(F)
        if nargout > 1 % for speed
            t = est .* sqrt((n - 2) ./ (1 - est.^2));
            p = 2 .* (1 - tcdf(abs(t),n-2));        % two-tailed p-value
        end
        %wts = [];
        %nonpar = 0;
        % this works,but is unsigned
        %[chi2,df,p,sig] = chi2test([x y],'obs'); %,wts,nonpar);
        %est = (chi2 ./ n) .^ .5;

    case {'rho','spearman'}

        % This method uses midrank method to handle ties; needs checking
        % vs. SPSS
        D = rankdata(x) - rankdata(y);
        est = 1 - (6*(D'*D)) ./ (n * (n^2-1));

        % same as : corrcoef(rankdata(x),rankdata(y));
        % but faster
        % tested against SPSS 11.04
        % problem with large n?

        % % n = 1000;
        % % tic, for i = 1:100, x = rand(n,1); y = x + rand(n,1);
        % % corrcoef(rankdata(x),rankdata(y));
        % % end, toc
        % % tic, for i = 1:100, x = rand(n,1); y = x + rand(n,1);
        % %     n = size(x,1);
        % % D = rankdata(x) - rankdata(y);
        % % 1 - (6*(D'*D)) ./ (n * ((n^2)-1));
        % % end, toc

    case {'taua','kendalla'}

        % checked 1/13/07 against http://www.wessa.net/rwasp_kendall.wasp
        % should probably be checked against SPSS

        [r,i] = sort(x);
        y = y(i);               % get ordered y; don't really need to rank

        % count number of ranks above each successive value of y
        for i = 1:(n - 1)
            P(i) = sum( y(i+1:end) > y(i) );
        end
        est = 4 * sum(P) ./ (n * (n - 1)) - 1;

        if nargout > 1 % for speed
            % actually a z-score
            t = 3 * est * (n * (n - 1)) .^.5 ./ (2 * (2 * n + 5)) .^.5;
            p = 1 - normcdf(abs(t));

            if p == 0, p = 10*eps; end
        end

    case {'tau','kendall','taub','kendallb'}

        % checked 1/13/07 against http://www.wessa.net/rwasp_kendall.wasp
        % should probably be checked against SPSS

        [r,i] = sort(x);
        y = y(i);               % get ordered y

        % count number of ranks above each successive value of y
        for i = 1:(n - 1)
            wh = r > r(i);  % find pts with x value greater than r(i); handle x ties
            P(i) = sum( y(wh) > y(i) ); % concordances
            N(i) = sum( y(wh) < y(i) ); % discordances
        end

        %% *** maybe ties have to include all values that are tied"??
        % http://www.statsdirect.com/help/nonparametric_methods/kend.htm
        xties = get_ties(x);
        yties = get_ties(y);

        num = sum(P) - sum(N);

        n1 = n * (n - 1) ./ 2;

        % From http://www.unesco.org/webworld/idams/advguide/Chapt4_2.htm
        % and http://www.statsdirect.com/help/nonparametric_methods/kend.htm
        den = sqrt(  (n1 - sum(xties)) * (n1 - sum(yties))  );

        % Tau b : from http://www.unesco.org/webworld/idams/advguide/Chapt4_2.htm
        %den = sqrt((sum(P) + sum(N) + (n - sum(xties))) * (sum(P) + sum(N) + (n - sum(yties))));

        est = num ./ den;

        if nargout > 1 % for speed
            % actually a z-score
            t = 3 * est * (n * (n - 1)) .^.5 ./ (2 * (2 * n + 5)) .^.5;
            p = 1 - normcdf(abs(t));
        end

    case {'gamma','kruskal'}

        [r,i] = sort(x);
        y = y(i);               % get ranks of y

        % count number of ranks above each successive value of y
        for i = 1:(n - 1)
            wh = r > r(i);  % find pts with x value greater than r(i); handle x ties
            P(i) = sum( y(wh) > y(i) ); % concordances
            N(i) = sum( y(wh) < y(i) ); % discordances
        end

        num = sum(P) - sum(N);
        den = sum(P) + sum(N);
        est = num ./ den;

    otherwise
        error('Unknown correlation option');
end

return

function tab = make_crosstabs(x,y)
% make table (crosstabs)
u1 = unique(x);
u2 = unique(y);
len1 = length(u1);
len2 = length(u2);
tab = zeros(len1,len2);
w = 1;                  % weights

for i = 1:length(u1)                   % rows are 0 then 1 on the first var, "Nos" then "Yesses"
    for j = 1:length(u2)               % for each column
        tab(i,j) = sum( (x == u1(i) & y == u2(j)) .* w );
    end
end

return

function ties = get_ties(x)
r = rankdata(x, 'nomidrank');
for i = 1:length(unique(r))
    ties(i) = sum(r == i);
end
ties = sum(ties .* (ties - 1) ./ 2);
return

function [rows,cols,npairs] = corrcoef_indices(nvars)

    if nvars > 1000
        fprintf('Setting up indices of matrix.');
        % added to deal with large datasets
        tmp = triu(true(nvars));

        for i = 1:nvars
            tmp(i, i) = 0;
        end

%         % create logical identity matrix of large size
%         t2 = eye(1000);
%         t2 = logical(t2);
%         eyemtx = t2;
%         while size(eyemtx, 2) < nvars
%             eyemtx = blkdiag(eyemtx, t2);
%         end
%         eyemtx = eyemtx(1:nvars, 1:nvars);
%
%         tmp = tmp - eyemtx;

        fprintf(1,'Done.\n');

    else
        % upper triangle only
        tmp = triu(ones(nvars));
        tmp = tmp - eye(nvars);
    end

    [rows,cols] = find(tmp);
    npairs = length(rows);

    return

function valmat = reconstruct(vals,nvars,npairs,rows,cols)

valmat = zeros(nvars);
for i = 1:npairs
    valmat(rows(i),cols(i)) = vals(i);
end
valmat = valmat + valmat';

return

function [pthr,sig] = fdr_correct_pvals(p,r)

    psq = p; psq(find(eye(size(p,1)))) = 0;
    psq = squareform(psq);
    pthr = FDR(p,.05);
    if isempty(pthr), pthr = 0; end

    sig = sign(r) .* (p < pthr);

    return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%--- help for generic/evaluate_general_restrictions ---
%
%  Evaluate general restrictions
%
%  ::
%
%    g = evaluate_general_restrictions(obj)
%
%  Args:
%
%     obj (rise | dsge | rfvar | svar): scalar of vector or RISE model objects
%
%  Returns:
%     :
%
%     - **g** [cell array] : value of restrictions for each object.
%
%  Note:
%
%     - The restrictions will be processed as g(x)<=0. But all the user has to
%       do is to put zero where the restrictions are not violated!!!
%
%

author: jmaih


text: %% SO3Fun

author: mtex-toolbox

repo: mtex

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


text:
%%
% \documentclass[12pt]{article}
%
% \title{ODEbox: A Toolbox for Ordinary Differential Equations\\
% Example 1}
%
% \author{Matthew Harker and Paul O'Leary\\
% Institute for Automation\\
% University of Leoben\\
% A-8700 Leoben,
% Austria\\
% URL: automation.unileoben.ac.at\\
% \\
% Original: January 9, 2013\\
% $\copyright$ 2013\\
% \\
% Last Modified: \today}
%%
%\section{Initial Value Problem}
%
% This file demonstrated the solution of initial value problems using the
% DOPBox and ODEBox toolbaxes. The example being considered here is:
%
% \begin{equation}
%    \ddot{y} + 6\,\dot{y} + 9 y = 0
%    \hspace{5mm}
%    \text{with}
%    \hspace{5mm}
%    y(0) = 10
%   \hspace{5mm}
%   \text{and}
%   \hspace{5mm}
%   \dot{y}(0) = -75
% \end{equation}
%
% The analytical solution to this problem is:
%
% \begin{equation}
%   y(x) = 10 \, e^{-3\,x} - 45 \,x {e^{-3\,x}}.
% \end{equation}
%
% This example has been taken from~\cite{adams}.
%
% This M-file demonstrated the new approach to solving this problem and
% compared the result with a classical runga-kutta solution as provided by
% MATLAB.
%
%%
% \section{Prepare the workspace}
%
close all;
clear all;
setUpGraphics(12) ;
%%
% Define the equation and solution as a strings for documentation
%
Eq = '$$\ddot{y} - 6\,\dot{y} - 9 y = 0$$, with, $$y(0) = 10$$ and $$\dot{y}(1) = -75$$' ;
Sol = '$$y(x) = 10 \, e^{-3\,x} - 45 \,x {e^{-3\,x}}$$';
%
%------------------------------------------------------------------
%%
% \section{Define the Matrix Linear Differential Operator}
%
% Define the number of points used in the solution
%
noPts = 85 ;
%
% Define the interval in which the problem is to be solved
%
xMin = 0 ;
xMax = 3 ;
%
% Use a nonlinear node placement. This ensures a higher density of nodes
% where the solution has a highed derivative.
%
x = linspace(0,1,noPts)';
x = 5 * x.^2;
%
% Setup the differentiating matrix, with support length 13.
%
ls = 13;
D = dopDiffLocal( x, ls, ls );
%
% Define the matrix C of constraints, i.e. the initial conditions
%
C = zeros(noPts,2) ;
%
C(1,1) = 1 ;
C(:,2) = D(1,:)' ;
%
d = [ 10 ; -75]  ;
%
% Compute the linear differential operator.
%
L = D*D + 6*D + 9*eye(noPts) ;
%%
%\section{Compute the Solutions}
%
% Solve the algebraic system of equations as a least squares problem. This
% is the actual computation of the solution of the differential equation.
%
y = odeLSE( L, [] , C, d );
%%
% For comparison solve using the MATLAB ode45 method.
% Runge-Kutta Solution:
%
[xm,Y] = ode45(@eulerODE01,[xMin xMax],d');
%
ym = Y(:,1);
%%
% Compute the Analytic Solution as an inline to compare the results
%
f = inline('10*exp(-3*t) - 45*t.*exp(-3*t)') ;
%
ya = f(x);
%%
% \section{Plot the Results}
%
% for plot purposes solve beyond the ends
%
xp = linspace(xMin-1/4,xMax+1/4,1000) ;
yp = f(xp);
%
% Present the solutions
%
setUpGraphics(12)
FigureSize=[1 1 10 6];
set(0,'DefaultFigureUnits','centimeters');
set(0,'DefaultFigurePosition',FigureSize);
set(0,'DefaultFigurePaperUnits','centimeters');
set(0,'DefaultFigurePaperPosition',FigureSize);
MyAxesPosition=[0.13 0.15 0.86 0.82];
set(0,'DefaultaxesPosition',MyAxesPosition);
%
fig1 = figure;
plot( xp, yp, 'k' ) ;
hold on
plot( x, y, 'ko', 'MarkerFaceColor', 'w' ) ;
%
plot( xm, ym,'kv');
legend('Analytical','New','Runge-Kutta','Location','NorthEast');
xlabel('$$x$$');
ylabel('$$y$$');
grid on;
range = axis;
plot( range(1:2), [0,0],'k');
axis([-0.5,3.5,-5,15]);
plot( x, -4.5*ones(size(x)), 'k.');
%
%\caption{Comparison of the analytical solution, the new numberical
% method and the Runga-Kutta solution.}
%%
%
setUpGraphics(12)
FigureSize=[1 1 10 6];
set(0,'DefaultFigureUnits','centimeters');
set(0,'DefaultFigurePosition',FigureSize);
set(0,'DefaultFigurePaperUnits','centimeters');
set(0,'DefaultFigurePaperPosition',FigureSize);
MyAxesPosition=[0.13 0.13 0.86 0.82];
set(0,'DefaultaxesPosition',MyAxesPosition);
%
fig2 = figure;
subplot(2,1,1)
plot( x, y - ya, 'k', 'MarkerFaceColor', 'w' ) ;
hold on
plot( x, y - ya, 'k.', 'MarkerFaceColor', 'w' ) ;
range = axis;
plot( range(1:2), [0,0],'k');
ylabel('$$e_{L}$$');
grid on;
%title(Sol);
%
subplot(2,1,2)
plot( xm, ym - f(xm), 'k', 'MarkerFaceColor', 'w' ) ;
hold on
plot( xm, ym - f(xm), 'k.', 'MarkerFaceColor', 'w' ) ;
hold on;
range = axis;
plot( range(1:2), [0,0],'k');
xlabel('$$x$$');
ylabel('$$e_{RK}$$');
grid on;
%
% \caption{Comparison the residual error for the new method (top)  (Bottom)and
% for the Runga-Kutta method. Note the residual is orders of magnitude smaller for the new method.}
%
%%
%\section{Save the figures to disk.}
%
fileType = 'eps';
printFigure( fig1, 'ivpExp3Sol', fileType);
printFigure( fig2, 'ivpExp3Err', fileType);
%
%% Define the Bibliography
%
% \bibliographystyle{IEEETran}
% \bibliography{odebib}


text:
function exact = p32_exact ( )

%*****************************************************************************80
%
%% P32_EXACT returns the exact integral for problem 32.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  exact = - 0.5 * ( exp ( pi ) + 1.0 );

  return
end


text:
function [indices] = sequence_query_tag(sequence, tag, subset)
% sequence_query_tag Find tag occurences in sequence
%
% The function returns indices of all frames that contain the specified tag.
%
% Input:
% - sequence (structure): A valid sequence structure.
% - tag (string): Name of the tag.
% - subset (vector): A set of frame numbers or a binary mask to be used
% when looking for a tag. Frame indices will be returned relative to that
% subset of frames.
%
% Output:
% - indices (integer): A vector of indices of frames that contain the given tag.

if nargin < 3
   subset = true(sequence.length, 1);
end

if strcmp(tag, 'all')
    if nargin < 3
        indices = 1:sequence.length;
    else
        indices = 1:numel(subset);
    end;
    return;
end;

if strcmp(tag, 'empty')
    indices = find(all(~sequence.tags.data, 2));
    return;
end;

tag_index = find(strcmp(sequence.tags.names, tag), 1);

if isempty(tag_index)
    indices = [];
else
    indices = find(sequence.tags.data(subset, tag_index));
end;

author: votchallenge

repo: toolkit-legacy


text:
function padua_test05 ( )

%*****************************************************************************80
%
%% PADUA_TEST05 tests PADUA_WEIGHTS and PADUA_WEIGHTS_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 June 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PADUA_TEST05\n' );
  fprintf ( 1, '  PADUA_WEIGHTS computes the weights of a Padua rule.\n' );
  fprintf ( 1, '  PADUA_WEIGHTS_SET looks them up in a table.\n' );

  for l = 3 : 4
    n = padua_order ( l );
    w1 = padua_weights ( l );
    w2 = padua_weights_set ( l );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Level %d  Padua weights\n', l );
    fprintf ( 1, '\n' );
    diff = 0.0;
    for j = 1 : n
      diff = max ( diff, abs ( w1(j) - w2(j) ) );
      fprintf ( 1, '  %4d  %14.6g  %14.6g\n', j, w1(j), w2(j) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Maximum difference = %g\n', diff );
  end

  return
end


text:
function x = AngNormalize(x)
%ANGNORMALIZE  Reduce angle to range [-180, 180)
%
%   X = ANGNORMALIZE(X) reduces angles in [-540, 540) to the range
%   [-180, 180).  X can be any shape.

  x(x >= 180) = x(x >= 180) - 360;
  x(x < -180) = x(x < -180) + 360;
end


text:
function [image,index]=find_color(image,map,tr,colors)
%[IMAGE,IND]=FIND_COLOR(A,MAP,THRESHOLD,COLOR VECTOR)
%Finds a color of an Image A and returns IMAGE, the
%original image with the pixels with the specified color set
%to zero, and IND, an index to the pixels that are were
%equal to the selected color.
%
%MAP is the input if the image is indexed, if the image is not indexed
%enter empty brackerts ( [] ), TR is the threshold tolerance for the
%color range, COLOR VECTOR is an additional argument that takes
%in the hard-coded values for the color, use empty brackects ( [] )
%if you wish to select the color from the origimal image using the mouse.
%
%Composite image, EXAMPLE:
%x=imread('football.jpg');
%subplot(2,1,1);imshow(x);
%title('Before');
%[image,ind]=find_color(x,[],90,[27 80 111]);
%subplot(2,1,2);imshow(image);
%title('After')
%figure
% y=imread('saturn.tif');
% image(328,438,:)=0;
% i= ~(im2bw(image,0.01));
% c=immultiply(y,i);
% c(:,:,2)=c;
% c(:,:,3)=c(:,:,2);
% y(:,:,2)=y;
% y(:,:,3)=y(:,:,1);
% C=imadd(image,c);
% imshow(C)

%by is

I=make_rgb(image,map);

if isempty(colors)
disp('Right click on the color you want off the image with the mouse:') ;

colors=impixel(I);
close
end

[i j]=size(colors);

disp('Wait....');

if i>1,
    colors=colors(1,:);
end
[image,index]= getcolor(I,colors,tr);

disp('Done')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [image,index]=getcolor(I,color,tr)

[i j k]=size(I);
R=color(1,1);
G=color(1,2);
B=color(1,3);
I=double(I);

mask=( abs(I(:,:,1)-R) <tr ) & ( abs(I(:,:,2)-G) <tr ) &( abs(I(:,:,3)-B) <tr );

I(:,:,1)=I(:,:,1).*(~mask);
I(:,:,2)=I(:,:,2).*(~mask);
I(:,:,3)=I(:,:,3).*(~mask);

image=uint8(I);
index=find(mask==1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function I=make_rgb(image,map)

[i j k]=size(image);

if (~isempty(map)),
    I=ind2rgb(image,map);
    I=im2uint8(I);

elseif (isempty(k)),
    I(:,:,2)=I;
    I(:,:,3)=I(:,:,1);

else
I=image;
end


text:
function vol=surfvolume(node,face,option)
%
% vol=surfvolume(node,face,option)
%
% calculate the enclosed volume for a closed surface
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%    node:  node coordinates
%    face:  surface triangle list
%
% output:
%    vol:   total volume of the enclosed space
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

face=face(:,1:3);

ed=surfedge(face);
if(~isempty(ed))
   error('open surface is detected, you have to close it first, consider meshcheckrepair() with meshfix option');
end

[no,el]=fillsurf(node,face);

vol=elemvolume(no,el);
vol=sum(vol);


text:
function h = bc_test02 ( x_num, x, t, h )

%*****************************************************************************80
%
%% BC_TEST02 evaluates the boundary conditions for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X_NUM, the number of nodes.
%
%    Input, real X(X_NUM,1), the node coordinates.
%
%    Input, real T, the current time.
%
%    Input, real H(X_NUM,1), the current heat values.
%
%    Output, real H(X_NUM,1), the current heat values, after boundary
%    conditions have been imposed.
%
  h(1,1)     = exact_test02 ( 1, x(1),     t );
  h(x_num,1) = exact_test02 ( 1, x(x_num), t );

  return
end


text:
function [x, infos] = div_admm_nmf(V, rank, in_options)
% Divergence-based ADMM algorithm for non-negative matrix factorization (KL-FPA-NMF).
%
% Inputs:
%       matrix      V
%       rank        rank
%       options     options
%           d_beta  : parameter of beta divergence
%                       (only beta=0 (IS) and beta=1 (KL) are supported)
%           rho     : smothing parameter
%           fixed   : vector containing the indices of the basis vectors in
%                       W to hold fixed (e.g., when W is known a priori)
% Output:
%       w           solution of w
%       infos       information
%
% References:
%       D.L. Sun and C. Fvotte,
%       "Alternating direction method of multipliers for non-negative matrix
%       factorization with the beta divergence,"
%       ICASSP 2014.
%
%
% This file is part of NMFLibrary
%
% This file has been ported from
%   nmf_kl_fpa.m at https://github.com/felipeyanez/nmf
%   by Felipe Yanez
%
%   Copyright (c) 2014-2016 Felipe Yanez
%
%   Permission is hereby granted, free of charge, to any person obtaining a
%   copy of this software and associated documentation files (the "Software"),
%   to deal in the Software without restriction, including without limitation
%   the rights to use, copy, modify, merge, publish, distribute, sublicense,
%   and/or sell copies of the Software, and to permit persons to whom the
%   Software is furnished to do so, subject to the following conditions:
%
%   The above copyright notice and this permission notice shall be included
%   in all copies or substantial portions of the Software.
%
%   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
%   OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
%   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
%   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
%   OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
%   ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
%   OTHER DEALINGS IN THE SOFTWARE.
%
%
% Ported by M.Horie and H.Kasai on June 30, 2022
%
% Change log:
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.rho = 1;
    local_options.metric_type = 'beta-div';
    local_options.d_beta = 1; % parameter of beta divergence (only beta=0 (IS) and beta=1 (KL) are supported)

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    W = init_factors.W;
    H = init_factors.H;

    % initialize
    method_name = sprintf('Div-ADMM (%s=%.1f)', options.metric_type, options.d_beta);
    epoch = 0;
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % initialize for this algorithm
    fixed=[];
    % get the vector of indices to update
    free = setdiff(1:rank, fixed);

    X = W*H;
    Wplus = W;
    Hplus = H;
    alphaX = zeros(size(X));
    alphaW = zeros(size(W));
    alphaH = zeros(size(H));

    if options.d_beta == 0 && ~isfield(in_options, 'rho') % when IS divergence (d_beta == 0), rho should be much higher.
        local_options.rho = 1000;
    end

    % store initial info
    clear infos;

    %[options.metric_type, options.metric.param] = check_divergence(options);

    [infos, f_val, optgap] = store_nmf_info(V, W, H, [], options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('%s: Epoch = 0000, cost = %.16e, optgap = %.4e\n', method_name, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        % update H
        H = (W'*W + eye(rank)) \ (W'*X + Hplus + 1/options.rho*(W'*alphaX - alphaH));

        % update W
        P = H*H' + eye(rank);
        Q = H*X' + Wplus' + 1/options.rho*(H*alphaX' - alphaW');
        W(:,free) = ( P(:,free) \ (Q - P(:,fixed)*W(:,fixed)') )';

        % update X (this is the only step that depends on beta)
        X_ap = W*H;
        if options.d_beta == 1

            b = options.rho*X_ap - alphaX - 1;
            X = (b + sqrt(b.^2 + 4*options.rho*V))/(2*options.rho);

        elseif options.d_beta == 0

            A = alphaX/options.rho - X_ap;
            B = 1/(3*options.rho) - A.^2/9;
            C = - A.^3/27 + A/(6*options.rho) + V/(2*options.rho);
            D = B.^3 + C.^2;

            X(D>=0) = nthroot(C(D>=0)+sqrt(D(D>=0)),3) + ...
                nthroot(C(D>=0)-sqrt(D(D>=0)),3) - ...
                A(D>=0)/3;

            phi = acos(C(D<0) ./ ((-B(D<0)).^1.5));
            X(D<0) = 2*sqrt(-B(D<0)).*cos(phi/3) - A(D<0)/3;

        else
            error('beta is not currently supported.')
        end

        % update for H_+ and W_+
        Hplus = max(H + 1/options.rho * alphaH, 0);
        Wplus = max(W + 1/options.rho * alphaW, 0);

        % update for dual variables
        alphaX = alphaX + options.rho * (X - X_ap);
        alphaH = alphaH + options.rho * (H - Hplus);
        alphaW = alphaW + options.rho * (W - Wplus);

        % measure gradient calc count
        grad_calc_count = grad_calc_count + m*n;

        % measure elapsed time
        elapsed_time = toc(start_time);

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, W, H, [], options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = W;
    x.W(:,free) = Wplus(:,free);
    x.H = Hplus;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [summa, index] = max_sum(v, n)
    [~, len] = size(v);
    summa = -inf;
    if n > len
        summa = 0;
    end
    index = -1;
    for i = 1:(len - n + 1)
        s = sum(v(i : i + n - 1));
        if s > summa
            summa = s;
            index = i;
        end
    end
end

author: anishLearnsToCode

repo: introduction-to-programming-with-matlab


text:
function [linecolors, linestyles, linewidths] = SetDataPlotLineStyles(handles, iCh)
global maingui

datatype    = GetDatatype(handles);
ml          = GetMeasurementList(handles);

linecolorsAvailable  = maingui.axesData.SDPairColors;
linestylesAvailable  = {'-',':','-.','--'};
linewidthsAvailable   = [1.5, 2.0, 1.5, 1.5, 1.5];

linecolors  = zeros(length(iCh), 3);
linestyles  = cell(1, length(iCh));
linewidths  = zeros(1, length(iCh));

for ii = 1:length(iCh)

    if datatype == maingui.buttonVals.RAW   ||   datatype == maingui.buttonVals.OD
        iSD = GetSelectedSDPairIndex(ml(iCh(ii),1), ml(iCh(ii),2));
        linecolors(ii, :)  = linecolorsAvailable(iSD, :);
        linestyles{ii}     = linestylesAvailable{ml(iCh(ii),4)};
        linewidths(ii)     = linewidthsAvailable(ml(iCh(ii),4));

    elseif datatype == maingui.buttonVals.CONC

        iHb = ml(iCh(ii),4);
        iSD = GetSelectedSDPairIndex(ml(iCh(ii),1), ml(iCh(ii),2));
        linecolors(ii, :)  = linecolorsAvailable(iSD, :);
        linestyles{ii}     = linestylesAvailable{iHb};
        linewidths(ii)     = linewidthsAvailable(iHb);

    elseif datatype == maingui.buttonVals.OD_HRF

        iSD = GetSelectedSDPairIndex(ml(iCh(ii),1), ml(iCh(ii),2));
        linecolors(ii, :)  = linecolorsAvailable(iSD, :);
        linestyles{ii}     = linestylesAvailable{ml(iCh(ii),4)};
        linewidths(ii)     = linewidthsAvailable(ml(iCh(ii),4));

    elseif datatype == maingui.buttonVals.CONC_HRF

        iHb = ml(iCh(ii),4);
        iSD = GetSelectedSDPairIndex(ml(iCh(ii),1), ml(iCh(ii),2));
        linecolors(ii, :)  = linecolorsAvailable(iSD, :);
        linestyles{ii}     = linestylesAvailable{iHb};
        linewidths(ii)     = linewidthsAvailable(iHb);

    end

end

% --------------------------------------------------------------------------
function iHb = HbLabel2Index(label)
hbTypes     = {'hbo','hbr','hbt'};
iHb = [];
for jj = 1:length(hbTypes)
    if ~includes(lower(label), lower(hbTypes{jj}))
        continue;
    end
    iHb = jj;
    break;
end

author: BUNPC