text:
classdef UniformGraph < AbstractGraph
%UNIFORMGRAPH is a class which generates graphs with seeds and unaries
% which are not dependent on image information (for instance unaries are
% only based on superpixel size). This class is derived from AbstractGraph.
% To undestand how unaries are generated, its worth seeing how
% AbstractGraph.prepare_graphs() calls create_unaries() in this class to
% compute unary values.
%
% @authors:     Ahmad Humayun
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

    properties
    end

    methods
        function obj = UniformGraph(seg_obj, segm_index)
            obj.init(seg_obj, segm_index);
        end

        function create_unaries(obj, graph_sub_method, seed_sets)
            switch graph_sub_method
                case 'internal'
                    obj.generate_uniform_unaries('internal', seed_sets);
                case 'external'
                    obj.generate_uniform_unaries('external', seed_sets);
                case 'subframe'
                    error('UniformGraph:compute_graph', ...
                          '''subframe'' not implemented yet');
                otherwise
                    error('UniformGraph:compute_graph', ...
                          '''%s'' is an invalid graph method', ...
                          graph_sub_method);
            end
        end

        function generate_uniform_unaries(obj, internal_external, ...
                                          seed_sets)
            inside_frame = ~obj.seg_obj.sp_frame_set;

            funcs = create_unary_aux_funcs(obj, internal_external);

            % find the offset according to the number of unary graphs
            % already added to graph_unaries_all
            ofst = size(obj.graph_unaries_all.nonlambda_s,2);

            % initialize the matrices for holding all the unaries (see the
            % comment for this structure in AbstractGraph.prepare_graphs())
            obj.graph_unaries_all.nonlambda_s = ...
                [obj.graph_unaries_all.nonlambda_s, zeros(size(seed_sets))];
            obj.graph_unaries_all.nonlambda_t = ...
                [obj.graph_unaries_all.nonlambda_t, zeros(size(seed_sets))];
            obj.graph_unaries_all.lambda_s = ...
                [obj.graph_unaries_all.lambda_s, zeros(size(seed_sets))];
            obj.graph_unaries_all.lambda_t = ...
                [obj.graph_unaries_all.lambda_t, zeros(size(seed_sets))];

            % iterate over all seed locations
            for i = 1:size(seed_sets, 2)
                if strcmp(internal_external, 'internal')
                    [nonlambda_s, lambda_t] = internal_uniform(obj, ...
                                                     seed_sets(:,i), ...
                                                     funcs.parametric_func);

                    obj.graph_unaries_all.nonlambda_s(:,ofst+i) = ...
                                                               nonlambda_s;
                    obj.graph_unaries_all.lambda_t(:,ofst+i) = lambda_t;
                else
                    [nonlambda_s, nonlambda_t, lambda_s, lambda_t] = ...
                        external_uniform(obj, seed_sets(:,i), ...
                                         inside_frame, ...
                                         funcs.parametric_func, ...
                                         funcs.frame_cost_func);

                    obj.graph_unaries_all.nonlambda_s(:,ofst+i) = ...
                                                               nonlambda_s;
                    obj.graph_unaries_all.nonlambda_t(:,ofst+i) = ...
                                                               nonlambda_t;
                    obj.graph_unaries_all.lambda_s(:,ofst+i) = lambda_s;
                    obj.graph_unaries_all.lambda_t(:,ofst+i) = lambda_t;
                end
            end
        end

        function [nonlambda_s, lambda_t] = internal_uniform(obj, ...
                                         the_rectangle, parametric_func)
            % see the comment for these output matrices in
            % AbstractGraph.prepare_graphs()
            nonlambda_s = zeros(length(the_rectangle),1);
            lambda_t = zeros(length(the_rectangle),1);

            bground_pixels = ~the_rectangle;
            nonlambda_s(the_rectangle) = inf;
            lambda_t(bground_pixels) = parametric_func(bground_pixels);
        end

        function [nonlambda_s, nonlambda_t, lambda_s, lambda_t] = ...
                external_uniform(obj, the_rectangle, inside_frame, ...
                                 parametric_func, frame_cost_func)
            % see the comment for these output matrices in
            % AbstractGraph.prepare_graphs()
            nonlambda_s = zeros(length(the_rectangle),1);
            nonlambda_t = zeros(length(the_rectangle),1);
            lambda_s = zeros(length(the_rectangle),1);
            lambda_t = zeros(length(the_rectangle),1);

            inside_frame_but_rect = inside_frame & ~the_rectangle;
            lambda_weights = parametric_func(inside_frame_but_rect);
            bground_pixels = ~(inside_frame | obj.seg_obj.sp_frame_set);

            nonlambda_s(the_rectangle) = inf;
            nonlambda_t(bground_pixels) = inf;
            lambda_s(inside_frame_but_rect) = lambda_weights;

            frame_cost = frame_cost_func(the_rectangle);

            nonlambda_t(obj.seg_obj.sp_frame_set) = ...
                frame_cost * obj.graph_seed_frame_weight;
            %hyp_conns{1,5} = {inside_frame_but_rect, obj.Background, 0, 1};
            lambda_t(inside_frame_but_rect) = lambda_weights;
        end
    end
end

author: Cloud-CV

repo: object-proposals


text:
function eeg_save_graphics(F,type,p,confirm)

% eeg_save_graphics - Save figure as graphics file
%
% Usage: eeg_save_graphics(F,type,p,confirm)
%
% 'F' is a figure handle, when empty uses 'gcf'
%
% 'type' is the format of the graphic file, eg:
%
%       'png'  - portable network graphics (default)
%       'jpg'  - JPEG image, quality level of 90
%       'tiff' - TIFF with packbits (lossless run-length encoding)
%       'eps'  - encapsulated postscript
%
% Note that tiff and eps files can be very large.
% The default resolution is 300 DPI.
%
% 'p' is the data struct of eeg_toolbox
%
% 'confirm' is a boolean to control direct or GUI
% based saving.  0 = direct, 1 = GUI confirmation.
%
% The output file is saved to p.volt.path with
% the filename of the voltage file and the timing
% of the graphic appended (assumes that the graphic
% is a topographic map).
%

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

% Licence:  GNU GPL, no implied or express warranties
% History:  04/2002, Darren.Weber_at_radiology.ucsf.edu
%           - extracted out of eeg_contours_engine
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if exist('F','var'),
    if isempty(F), F = gcf; end
else
    F = gcf;
end
if exist('type','var'),
    if isempty(type), type = 'png'; end
else
    type = 'png';
end
if exist('confirm','var'),
    if isempty(confirm), confirm = 0; end
else
    confirm = 0;
end

[path,file,ext] = fileparts(strcat(p.volt.path,filesep,p.volt.file));
if p.topoView,
    file = strcat(file, sprintf('_%s_%08.2f',p.topoView,p.volt.sampleTime),'.',type);
else
    file = strcat(file, sprintf('_%08.2f',p.volt.sampleTime),'.',type);
end
file = fullfile(path,file);

if confirm,
    [filename, filepath] = uiputfile(file, 'EEG Save Graphic');
    if ~isequal(filename,0),
        file = fullfile(filepath,filename);
        fprintf('EEG_SAVE_GRAPHIC: Saving to:...\n...%s\n',file);
        saveas(F,file);
    end

else
    type = lower(type);
    if strmatch(type,'png'),
        driver = '-dpng';
        option1 = '-noui';
        option2 = '-r300';
        %option3 = '-zbuffer';
        fprintf('EEG_SAVE_GRAPHIC: Saving to:...\n...%s\n',file);
        print(F,driver,option1,option2,file);
    elseif strmatch(type,'jpeg'),
        driver  = '-djpeg90'; % JPEG image, quality level of nn (90 here)
        option1 = '-r300';
        fprintf('EEG_SAVE_GRAPHIC: Saving to:...\n...%s\n',file);
        print(F,driver,option1,file);
    elseif strmatch(type,'eps'),
        driver = '-depsc2'; % Encapsulated Level 2 Color PostScript
        option1 = '-r300';
        option2 = '-tiff';
        fprintf('EEG_SAVE_GRAPHIC: Saving to:...\n...%s\n',file);
        print(F,driver,option1,option2,file);
    elseif strmatch(type,'tiff'),
        driver = '-dtiff'; % TIFF with packbits (lossless run-length encoding) compression
        option1 = '-r300';
        fprintf('EEG_SAVE_GRAPHIC: Saving to:...\n...%s\n',file);
        print(F,driver,option1,file);
    else
        fprintf('EEG_SAVE_GRAPHIC: Saving to:...\n...%s\n',file);
        saveas(F,file);
    end
end

return

author: PatternRecognition

repo: OpenBMI


text:
function w = rcoswin(N)
%RCOSWIN Raised cosine window.
%
%   IOSR.DSP.RCOSWIN(N) returns an N-point raised cosine window.
%
%   See also GAUSSWIN, CHEBWIN, KAISER, TUKEYWIN, WINDOW.

%   Copyright 2016 University of Surrey.

    assert(isscalar(N) && round(N)==N, 'iosr:rcoswin:invalidN', 'N must be a scalar and whole number.')
    assert(N > 0, 'iosr:rcoswin:invalidN', 'N must be greater than 0.')

    switch N
        case 1
            w = 1;
        case 2
            w = [.5; .5];
        otherwise
            r = linspace(-pi,pi,N)';
            w = (cos(r)+1)/2;
    end

end

author: IoSR-Surrey

repo: MatlabToolbox


text:
classdef Var < Param
%VAR Defines a network variable explicitly
%   Var(VALUE) or Var('value', VALUE) defines an explicit variable in the
%   network, with a given value.
%
%   Most of the time this is not needed; variables are created implicitly
%   to hold the values of any Inputs, Params, and output values computed by
%   Layers (typically instantiated with overloaded function calls).
%   Constants such as strings or Matlab arrays are stored as arguments to
%   any function call that uses them (in a Layer).
%
%   An explicit variable (Var) acts like a constant argument to a Layer,
%   but may still be useful for the following reasons:
%   1. Its value can be moved to and from the GPU using NET.USEGPU.
%   2. Its value can be set after network construction, with NET.SETVALUE.
%   3. Its derivative can be retrieved with NET.GETDER.
%
%   Internally, a Var is just a Param with trainMethod set to 'none'.

% Copyright (C) 2016 Joao F. Henriques.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

  methods
    function obj = Var(varargin)
      if isscalar(varargin)
        % special syntax, just passed in the value
        varargin = {'value', varargin{1}} ;
      end

      obj@Param('trainMethod','none', 'weightDecay',0, ...
        'learningRate',0, varargin{:})
    end
  end
end

repo: MFT


function price = CrankNicolsonFD_BlackScholes_func( S_0, K, r, T, sigma, call, dS, dt, Smax, Smin)
% Description: Crank-Nicolson PDE Finite Difference method to price European Option in Black Scholes Model
% Author: Justin Kirkby
M = round((Smax - Smin)/dS); % grid points
N = round(T/dt); % time steps

dS = (Smax - Smin) / M; % readjust
dt = T / N;  % readjust

vals = zeros(M+1, N+1);
vS = linspace(Smin, Smax, M+1)';
vI = vS / dS;
vJ = 0:N;

% Boundary Conditions
if call == 1  % call option
    vals(:, N+1) = max(vS - K, 0);
    vals(1, :) = 0;
    vals(M+1, :) = Smax - K*exp(-r*dt*(N - vJ));
else % put option
    vals(:, N+1) = max(K - vS, 0);
    vals(1, :) = K*exp(-r*dt*(N - vJ));
    vals(M+1, :) = 0;
end

% Tridiagonal Coefficients
a = 0.25 * dt * (sigma^2 *(vI.^2) - r*vI);
b = -dt*0.5*(sigma^2*(vI.^2) + r);
c = 0.25*dt*(sigma^2*(vI.^2) + r*vI);

M1 = -diag(a(3:M), -1) + diag(1 - b(2:M)) - diag(c(2:M-1),1);
[L,U] = lu(M1);
M2 = diag(a(3:M), -1) + diag(1 + b(2:M)) + diag(c(2:M-1),1);

% Solve systems (backward in time)
for j = N:-1:1
    vals(2:M,j) = U \ (L \ (M2*vals(2:M,j+1)));
end

price = interp1(vS, vals(:,1), S_0);

end


text:
h=[4,-4;-4,8];
f=[-6;-3];
a=[1,1;4,1];
b=[3;9];
[x,value]=quadprog(h,f,a,b,[],[],zeros(2,1))

author: Eurus-Holmes


text:
function [xstar,fstar]=Bisection(a0,b0)

global counterf
counterf=0;

% initialization
epsilon=0.00001;
L=1e-4;

ak=a0;
bk=b0;

% iteration
while 1
    % stop criteria
    len=bk-ak;
    if len<L
        xstar=(ak+bk)/2;
        fstar=obj1(xstar);
        break;
    end

    mean=(ak+bk)/2;
    lambdak=mean-epsilon;
    muk=mean+epsilon;
    %compute two points every iteration
    flambdak=obj1(lambdak);
    fmuk=obj1(muk);

    if flambdak>fmuk
        ak=lambdak;
    else
        bk=muk;
    end
end

counterf

repo: Optimization-Theory-and-Algorithm


text:
%kcpRGBNTSCtoRC 'converts an image from RGB NTSC to Retinal Cone '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros cpRGBNTSCtoRC.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kcpRGBNTSCtoRC(i, {'i','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% cpRGBNTSCtoRC - converts an image from RGB NTSC to Retinal Cone
%
%  DESCRIPTION
% call ckfromRGBNTSC with the parameter -retinalcone to convert an image from RGB NTSC to Retinal Cone
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
% All references for the COLOR toolbox are on the file $COLOR/Readme.html
%
%  COPYRIGHT
% Copyright (C) 1997 Rafael Santos. Khoros (C) Khoral Research, Inc.
%

function varargout = kcpRGBNTSCtoRC(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kcpRGBNTSCtoRC(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output'};
maxval={0,0};
minval={0,0};
istoggle=[0,0];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'ckfromRGBNTSC"  -retinalcone'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [e, varargout] = neterr(w, net, x, t)
%NETERR	Evaluate network error function for generic optimizers
%
%	Description
%
%	E = NETERR(W, NET, X, T) takes a weight vector W and a network data
%	structure NET, together with the matrix X of input vectors and the
%	matrix T of target vectors, and returns the value of the error
%	function evaluated at W.
%
%	[E, VARARGOUT] = NETERR(W, NET, X, T) also returns any additional
%	return values from the error function.
%
%	See also
%	NETGRAD, NETHESS, NETOPT
%

%	Copyright (c) Ian T Nabney (1996-2001)

errstr = [net.type, 'err'];
net = netunpak(net, w);

[s{1:nargout}] = feval(errstr, net, x, t);
e = s{1};
if nargout > 1
  for i = 2:nargout
    varargout{i-1} = s{i};
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [mmHg] = Pa2mmHg(Pa)
% Convert pressure from pascals to millimeters of mercury.
% Chad Greene 2012
mmHg = Pa*0.00750062;

author: Sable

repo: mcbench-benchmarks


text:
function [kPa] = Torr2kPa(Torr)
% Convert pressure from torr to kilopascals
% Chad Greene 2012
kPa = Torr*.133322;


text:
function path=VMR2mrGray(VMRfile,mmPerPix);
% function path=VMR2mrGray(VMRfile,mmPerPix);
% Function to convert between Brainvoyager VMR file format and mrGray
% Last edited : $Date: 2007/07/05 19:51:58 $
% ARW 100300
if ~exist('mmPerPix','var')
   mmPerPix = [240/256 240/256 1.2];
   disp(['mmPerPix defaulting to [ ' num2str(mmPerPix,'%.4f ') ...
         	'].  I hope this is correct!']);
end

fid=fopen(VMRfile,'r');
header=fread(fid,3,'int16');
mainImg=fread(fid,'uchar');
fclose(fid);
mainImg=reshape(mainImg,[header(1),header(2),header(3)]);
disp('Flipping images...');

for thisIm=1:header(3)
   mainImg(:,:,thisIm)=rot90(fliplr(mainImg(:,:,thisIm)));
end

path = writeVolAnat(mainImg, mmPerPix);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function MatingPool = MatingSelection(Obj,RefPoint,Range)
% The mating selection of AR-MOEA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Calculate the fitness of each feasible solution based on IGD-NS
        % Calculate the distance between each solution and point
        N = size(Obj,1);
        Distance    = CalDistance(Obj-repmat(Range(1,:),N,1),RefPoint);%N*NR
        Convergence = min(Distance,[],2);% The minimum Angle between each solution and the reference point
        [dis,rank]  = sort(Distance,1);% In columns, for reference points
        % Calculate the fitness of noncontributing solutions
        Noncontributing = true(1,N);
        Noncontributing(rank(1,:)) = false;
        METRIC   = sum(dis(1,:)) + sum(Convergence(Noncontributing));
        % The sum of the minimum angles of all reference points plus the sum of the minimum angles of all non-contributing solutions
        fitness  = inf(1,N);
        fitness(Noncontributing) = METRIC - Convergence(Noncontributing);
        % Calculate the fitness of contributing solutions
        for p = find(~Noncontributing)
            temp = rank(1,:) == p;
            noncontributing = false(1,N);
            noncontributing(rank(2,temp)) = true;% After removing this contribution point, the new contribution point Index
            noncontributing = noncontributing & Noncontributing;% Judge whether the new contribution point is the original non-contribution point or not
            fitness(p) = METRIC - sum(dis(1,temp)) + sum(dis(2,temp)) - sum(Convergence(noncontributing));
            % When this contribution point is removed, the METRIC adjusts
        end

    %% Combine the fitness of feasible solutions with the fitness of infeasible solutions
    Fitness = fitness;

    %% Binary tournament selection
    CV=zeros(N,1);
    MatingPool = TournamentSelection(2,ceil(N/2)*2,CV,-Fitness);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Max manipulability index ALONG A LINE.
% Start by maximum global manipulability at a pose
% Start by maximizing manipulability locally at each time step
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [qq, manips]=path_planning_local_max%(direction)
close all;
global robot

experiment_name = 'min_manip_local.mat';

height1 = 2; %m
height2 = 2.5; %m
x1 = -1;
x2 = 2;

%the orientation needed
phi = pi/2;
p0 = [x1 height1 0]';
pf = [x2 height2 0]';

T0 = build_T(p0, phi);
Tf = build_T(pf, phi);

%this defines the line to follow in direction
deltaV = pf - p0;
deltaV = deltaV/norm(deltaV);
%and total length of the movement
total_length = norm(pf - p0); %m
%the movement step in m
ds=0.1;
mov = 0.0:ds:total_length;

line0 = [p0 pf];
plot(line0(1,:),line0(2,:),'k')

%initial pose and manipulability
qq = [];

q0 = [pi/2 -pi -pi pi/4]';
q = q0;
Ti = T0;
for i=1:length(mov)
    fprintf('Move %d out of %d\n', i, length(mov))
    %update to next point in trajectory
    Ti(1:3,4) = p0 + mov(i)*deltaV;
    %(robot,  Tf, q0, restriction)
    q = inversekinematic_4dofplanar(robot, Ti, q);
    drawrobot3d(robot, q)
    q = max_manipulability_local(robot, q, -1);
    %drawrobot3d(robot, q)
    draw_axes(Ti, 'Xpiece', 'Ypiece', 'Zpiece', 1.2);
    plot3(line0(1,:),line0(2,:),line0(3,:),'k')
    qq = [qq q];
end

animate_local(robot, qq, line0)

manips = compute_manip(robot, qq);
figure,
plot(manips)
title('manipulability index at each movement')
figure,
plot(qq')
title('joint positions')
legend('q_1', 'q_2','q_3','q_4','q_5','q_6','q_7')

save(experiment_name)

% Simple iteration along the null space to increase manipulability
function [q] = max_manipulability_local(robot, q, sign_max)

delta_time = 0.01;
while 1
%    [qd_null] = compute_null_space(robot, q);
[qd_null] = null_space_4dof(robot, q);
    delta_manip = compute_delta_manip(robot, q, qd_null, 0.01);
    % add a sign for maximization/minimization
    % the constant is based on the rate of delta_manip also
    x = sign_max*sign(delta_manip)*qd_null;
    q = q + x*delta_time;
    if abs(delta_manip) < 1e-2
        break
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computes a gradient so as to let the manipulability be improved
% q is the current joint position
% qd is the instantaneous speed that lies in the null space
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function delta_manip = compute_delta_manip(robot, q, qd, delta)
%compute first manipulability index
m0 = compute_manip(robot, q);
%move differentially along the null space.
q = q + delta*qd;
%compute second manipulability index
m1 = compute_manip(robot, q);
%delta_manip = trace(inv(J*J')*(Jd*J'+J*Jd'));
%return difference
delta_manip=(m1-m0)/delta;

%
% Computes manipulability of poses given in qs
% Caution: only taking into account vx, vy and wz
%
function manips = compute_manip(robot, qs)
manips = [];
for i=1:size(qs,2)
    %compute current manip
    J = manipulator_jacobian(robot, qs(:,i));
    J = [J(1:2,:); J(6,:)];
    manip = sqrt(det(J*J'));
    manips = [manips manip];
end

function T = build_T(p, phi)
T = [cos(phi) -sin(phi) 0 p(1);
     sin(phi) cos(phi) 0 p(2);
     0            0     1  p(3);
     0             0    0   1];

function animate_local(robot, qq, line_work)
for i=1:1:size(qq,2)
     drawrobot3d(robot, qq(:,i))
     plot(line_work(1,:),line_work(2,:))
end

repo: ARTE


text:
function y = sinfun1(M)
%Sample function used in Chapter 2.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt
x = 0:M - 1;
for k = 1:numel(x)
   y(k) = sin(x(k)/(100*pi));
end

author: dipum

repo: dipum-toolbox

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


text:
function offset = find_offset(sig, sig_fs, ref_signal, ref_signal_fs)
    % Find the lag between two signals using cross-correlation
    sig_us = resample(sig, ref_signal_fs, sig_fs);
    sig_us = normalization(sig_us);
    [c, lags] = xcorr(normalization(ref_signal), sig_us);
    [~, ind] = max(c);
    offset = lags(ind) + 1; % add 1 for Matlab indexing
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function cvx_setnpos( x )

global cvx___
x  = x.basis_;
xc = x(1,:);
x  = x(2:end,:);
tt = ( xc <= 0 ) & ( sum( x ~= 0 ) == 1 );
[ tx, dummy ] = find( x(:,tt) ); %#ok
tx = tx + 1;
persistent mneg;
if isempty( mneg ),
    mneg = int8([1,2,2,22,5,5,6,8,8,9,11,11,12,22,22,22,22,22,22,22,22,22]);
end
cvx___.classes(tx) = mneg(cvx___.classes(tx));

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

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function g = heavisideNoiseGradientParam(noise, mu, varsigma, y)

% HEAVISIDENOISEGRADIENTPARAM Gradient of the heaviside noise model's parameters.

% IVM

c = y./sqrt(varsigma);
denom = zeros(size(c));
u = zeros(size(c));
for i = 1:size(mu, 2)
  u(:, i) = c(:, i).*(mu(:, i) + noise.bias(i));
  denom(:, i) = ((1-2*noise.eta)*cumGaussian(u(:, i))+noise.eta);
end

gnoise.bias = (1-2*noise.eta).*sum(c.*ngaussian(u)./denom, 1);
gnoise.eta = sum(sum((-2*cumGaussian(u)+1)./denom, 1));
gnoise.eta = gnoise.eta.*(noise.eta.*(1-2*noise.eta));
g = [gnoise.eta gnoise.bias];

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [cst,cstJac] = autoGen_cst_footVel(q1p,q2p,q4p,q5p,q1m,q2m,q4m,q5m,dq1p,dq2p,dq4p,dq5p,dq1m,dq2m,dq4m,dq5m,l1,l2,l4,l5)
%AUTOGEN_CST_FOOTVEL
%    [CST,CSTJAC] = AUTOGEN_CST_FOOTVEL(Q1P,Q2P,Q4P,Q5P,Q1M,Q2M,Q4M,Q5M,DQ1P,DQ2P,DQ4P,DQ5P,DQ1M,DQ2M,DQ4M,DQ5M,L1,L2,L4,L5)

%    This function was generated by the Symbolic Math Toolbox version 6.2.
%    22-Oct-2015 19:14:34

t2 = sin(q1p);
t3 = sin(q2p);
t4 = sin(q4p);
t5 = sin(q5p);
t6 = sin(q1m);
t7 = sin(q2m);
t8 = sin(q4m);
t9 = sin(q5m);
cst = [dq1p.*l1.*t2+dq2p.*l2.*t3-dq4p.*l4.*t4-dq5p.*l5.*t5;-dq1m.*l1.*t6-dq2m.*l2.*t7+dq4m.*l4.*t8+dq5m.*l5.*t9];
if nargout > 1
    cstJac = reshape([0.0,0.0,dq1p.*l1.*cos(q1p),0.0,dq2p.*l2.*cos(q2p),0.0,0.0,0.0,-dq4p.*l4.*cos(q4p),0.0,-dq5p.*l5.*cos(q5p),0.0,l1.*t2,0.0,l2.*t3,0.0,0.0,0.0,-l4.*t4,0.0,-l5.*t5,0.0,0.0,0.0,0.0,-dq1m.*l1.*cos(q1m),0.0,-dq2m.*l2.*cos(q2m),0.0,0.0,0.0,dq4m.*l4.*cos(q4m),0.0,dq5m.*l5.*cos(q5m),0.0,-l1.*t6,0.0,-l2.*t7,0.0,0.0,0.0,l4.*t8,0.0,l5.*t9],[2, 22]);
end

author: MatthewPeterKelly

repo: OptimTraj


text:
function x = is_search(x)
% Determine whether some object is a search() expression.
x = all(isfield(x,{'head','parts'})) && numel(x)==1 && strcmp(char(x.head),'search');

repo: REST


text:
function [TCER, TCER_rand] = gsp_learn_tcer(G, X, params)
%GSP_LEARN_TCER Learn the Total Cumulative Energy Residual
%   Usage: TCER = gsp_learn_tcer(G, X);
%          TCER = gsp_learn_tcer(G, X, params);
%          [TCER, TCER_rand] = gsp_learn_tcer(...);
%
%   Input parameters:
%       G           : the graph
%       X           : a data matrix
%   Output parameters:
%       TCER        : the computed total commulative energy residual for G
%       TCER_rand   : the computed total commulative energy residual for a random basis
%
%   Total Cumulative Energy Residual (TCER): This is a number from [0 to 1]
%   describing how well a graph fits a given data matrix X, or
%   distribution.
%
%   Optional parameters
%   -------------------
%
%   * *params.s*      : s = svd(X);
%   * *params.verbose*: 0 = nothing, 1 = plot cum energy, 2 = plot against
%     random basis as a baseline (default 0)
%   * *params.sort*   : sort basis columns to get the minimum residual
%     possible (default 0)
%
%   See also: gsp_good_graph_index, gsp_stationarity_ratio
%

% Author  : Andreas Loukas, Vassilis Kalofolias
% Date    : 15 Nov 2016

% Handle input
if nargin < 3, params = struct(); end
if not(isfield(params, 'verbose')); params.verbose = 0; end;
if not(isfield(params, 'sort'));    params.sort = 0;    end;
if not(isfield(params, 's'));       params.s = svd(full(X));  end;

s = params.s;

if numel(s)< G.N
    s = [s; zeros(G.N-numel(s),1)];
end

if not(isfield(G, 'U'))
    G = gsp_compute_fourier_basis(G);
end

% frobenius norm of X
norm_X_fro = norm(s);

% sorting already done by svd.
% var_cumsum_svd = cumsum(sort(s.^2, 'descend')) ./ norm_X_fro^2;
var_cumsum_svd = cumsum(s.^2) ./ norm_X_fro^2;

% frobenius norm of X
norm_X_fro = norm(X, 'fro');

% for the graph
if params.sort
    var_cumsum_s_G = cumsum(sort(sum(abs(G.U' * X).^2, 2) ./ norm_X_fro^2, 'descend'));
else
    var_cumsum_s_G = cumsum(sum(abs(G.U' * X).^2, 2) ./ norm_X_fro^2);
end

% area under the curve
AUC_SVD   = sum(var_cumsum_svd);
AUC_G     = sum(var_cumsum_s_G);
TCER      = 1 - AUC_G/AUC_SVD;

% compute also the result for a random basis
AUC_RB    = sum(linspace(1/G.N, 1, G.N));
TCER_rand = (AUC_SVD - AUC_RB)/(AUC_SVD);

if params.verbose >= 1

    figure; hold on;

    plot(var_cumsum_svd)
    plot(var_cumsum_s_G)
    [Urand, ~] = qr(randn(G.N));        % random basis
    var_cumsum_s_rand = cumsum(sum(abs(Urand' * X).^2, 2) ./ norm_X_fro^2);
    plot(var_cumsum_s_rand)

    xlabel('k (number of basis vectors used)');
    ylabel('cumulative variance');
    legend('svd', 'G', 'random basis')
end

% plot good quality residual figure
if params.verbose >= 2
    figure;
    plot(var_cumsum_svd);
    hold on;
    %     plot(cum_energy_graph);
    plot(var_cumsum_s_G, 'color', [0.8500    0.3250    0.0980]);
    n_nodes = length(var_cumsum_svd);
    xf = [1:n_nodes, fliplr(1:n_nodes)];
    yf = [zeros(n_nodes,1); flipud(var_cumsum_svd)];
    fill(xf, yf, [     0    0.4470    0.7410],'facealpha', .2);

    xf = [1:n_nodes, fliplr(1:n_nodes)];
    yf = [zeros(n_nodes,1); flipud(var_cumsum_s_G)];
    fill(xf, yf, [0.8500    0.3250    0.0980],'facealpha', .2);

    plot(var_cumsum_s_G, 'color', [0.8500    0.3250    0.0980]);
    plot(var_cumsum_svd, 'color', [     0    0.4470    0.7410]);
    legend('Ground truth covariance', 'Graph', 'Ground truth total C.E.', 'Graph total C.E.','location', 'se')
%     title(['graph from ', num2str(k), ' samples'])
    ylabel('Expected cumulative energy');
    text(n_nodes/2, .5, ['Residual: ', num2str(100*TCER, 3), '%']);
    xlabel('k');
    ylim([0,1]);
end

repo: gspbox


text:
function [h,filtertype] = gsp_jtv_filter_array(G,g,filtertype)
%GSP_JTV_FILTER_ARRAY Convert ts/js filters to -array filters
%   Usage: [h,filtertype] = gsp_jtv_filter_array(G,g, filtertype)
%
%   Input parameters:
%         G          : Graph
%         g          : Cell array of time-vertex filters
%         filtertype : Filter domain (ts,js)
%   Output parameters:
%         h          : Cell array of graph filterbank
%         filtertype : Filter domain (ts-array,js-array)
%
%   Convert ts/js filters to -array filters
%

% Author :  Francesco Grassi, Nathanael Perraudin
% Date : September 2016

if ~iscell(g)
   g = {g};
end

if ~gsp_check_filtertype(filtertype,{'ts','js'})
    error('Invalid filtertype.');
end

T  = G.jtv.T;
Nf = numel(g);

switch filtertype
    case 'ts'
        v = gsp_jtv_ta(G);
    case 'js'
        v = gsp_jtv_fa(G);
end

h  = cell(Nf,T);
for n=1:Nf
    for ii = 1:T
        h{n,ii} = @(x) g{n}(x,v(ii));
    end
end

filtertype = [filtertype '-array'];

end


text:
function s = msize(a);

    asdfdsa;

repo: OpenBMI


text:
function [BFtest] = BFtest(X,alpha)
%Brown-Forsythe's Test for Homogeneity of Variances.
%[In the Brown-Forsythe's test the data are transforming to yij = abs[xij - median(xj)]
%and uses the F distribution performing an one-way ANOVA using y as the
%dependent variable. The Brown-Frosythe statistic is corrected for artificial zeros
%occurring in odd sized samples.
%
%   Syntax: function [BFtest] = BFtest(X,alpha)
%
%     Inputs:
%          X - data matrix (Size of matrix must be n-by-2; data=column 1, sample=column 2).
%       alpha - significance level (default = 0.05).
%     Outputs:
%          - Sample variances vector.
%          - Whether or not the homoscedasticity was met.
%
%    Example: From the example 10.1 of Zar (1999, p.180), to test the Brown-Forsythe's
%             homoscedasticity of data with a significance level = 0.05.
%
%                                 Diet
%                   ---------------------------------
%                       1       2       3       4
%                   ---------------------------------
%                     60.8    68.7   102.6    87.9
%                     57.0    67.7   102.1    84.2
%                     65.0    74.0   100.2    83.1
%                     58.6    66.3    96.5    85.7
%                     61.7    69.8            90.3
%                   ---------------------------------
%
%           Data matrix must be:
%            X=[60.8 1;57.0 1;65.0 1;58.6 1;61.7 1;68.7 2;67.7 2;74.0 2;66.3 2;69.8 2;
%            102.6 3;102.1 3;100.2 3;96.5 3;87.9 4;84.2 4;83.1 4;85.7 4;90.3 4];
%
%     Calling on Matlab the function:
%             BFtest(X)
%
%       Answer is:
%
% The number of samples are: 4
%
% ----------------------------
% Sample    Size      Variance
% ----------------------------
%   1        5         9.3920
%   2        5         8.5650
%   3        4         7.6567
%   4        5         8.3880
% ----------------------------
%
% Brown-Forsythe's Test for Equality of Variances F=0.0831, df1= 3, df2=15
% Probability associated to the F statistic = 0.9682
% The associated probability for the F test is larger than 0.05
% So, the assumption of homoscedasticity was met.
%

%  Created by A. Trujillo-Ortiz and R. Hernandez-Walls
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             atrujo@uabc.mx
%
%  April 19, 2003.
%
%  To cite this file, this would be an appropriate format:
%  Trujillo-Ortiz, A. and R. Hernandez-Walls. (2003). BFtest: Brown-Forsythe's test for homogeneity of
%    variances. A MATLAB file. [WWW document]. URL http://www.mathworks.com/matlabcentral/fileexchange/
%    loadFile.do?objectId=3412&objectType=FILE
%
%  References:
%
%  Brown, M. B. and Forsythe, A. B. (1974), Robust Tests for
%           the Equality of Variances. Journal of the American
%           Statistical Association, 69:364-367.
%  Zar, J. H. (1999), Biostatistical Analysis (2nd ed.).
%           NJ: Prentice-Hall, Englewood Cliffs. p. 180.
%

if nargin < 2,
   alpha = 0.05;
end

Y=X;
k=max(Y(:,2));
fprintf('The number of samples are:%2i\n\n', k);

%Brown-Forsythe's Procedure.
n=[];s2=[];Z=[];
indice=Y(:,2);
for i=1:k
   Ye=find(indice==i);
   eval(['Y' num2str(i) '=Y(Ye,1);']);
   eval(['mdY' num2str(i) '=median(Y(Ye,1));']);
   eval(['n' num2str(i) '=length(Y' num2str(i) ') ;']);
   eval(['s2' num2str(i) '=(std(Y' num2str(i) ').^2) ;']);
   eval(['Z' num2str(i) '= abs((Y' num2str(i) ') - mdY' num2str(i) ');']);
   eval(['xn= n' num2str(i) ';']);
   eval(['xs2= s2' num2str(i) ';']);
   eval(['x= Z' num2str(i) ';']);
   n=[n;xn];s2=[s2;xs2];Z=[Z;x];
end

fprintf('-----------------------------\n');
disp(' Sample    Size      Variance')
fprintf('-----------------------------\n');
for i=1:k
   fprintf('   %d       %2i         %.4f\n',i,n(i),s2(i))
end
fprintf('-----------------------------\n');
disp(' ')

Y=[Z Y(:,2)];

%Correction for artificial zeros occurring in odd sized samples.
for i=1:k
   Ye=find(Y(:,2)==i);
   ncero=find((Y(:,1)==0)&(Y(:,2)==i));
   Y(ncero,1)=999;
   Y(ncero,1)=min(Y(Ye,1));
end

%Analysis of variance procedure.
C=(sum(Y(:,1)))^2/length(Y(:,1)); %correction term.
SST=sum(Y(:,1).^2)-C; %total sum of squares.
dfT=length(Y(:,1))-1; %total degrees of freedom.

indice=Y(:,2);
for i=1:k
   Ye=find(indice==i);
   eval(['A' num2str(i) '=Y(Ye,1);']);
end

A=[];
for i=1:k
   eval(['x =((sum(A' num2str(i) ').^2)/length(A' num2str(i) '));']);
   A=[A,x];
end

SSA=sum(A)-C; %sample sum of squares.
dfA=k-1; %sample degrees of freedom.
SSE=SST-SSA; %error sum of squares.
dfE=dfT-dfA; %error degrees of freedom.
MSA=SSA/dfA; %sample mean squares.
MSE=SSE/dfE; %error mean squares.
F=MSA/MSE; %Brown-Forsythe's F-statistic.
v1=dfA;df1=v1;
v2=dfE;df2=v2;

P = 1 - fcdf(F,v1,v2);  %probability associated to the F-statistic.

fprintf('Brown-Forsythe''s Test for Equality of Variances F=%3.4f, df1=%2i, df2=%2i\n', F,df1,df2);
fprintf('Probability associated to the F statistic = %3.4f\n', P);

if P >= alpha;
  fprintf('The associated probability for the F test is equal or larger than% 3.2f\n', alpha);
  fprintf('So, the assumption of homoscedasticity was met.\n');
else
  fprintf('The associated probability for the F test is smaller than% 3.2f\n', alpha);
  fprintf('So, the assumption of homoscedasticity was not met.\n');
end


text:
function [pot, MEU] = marginal_nodes(engine, d)

C = length(cliques);
%clpot = init_clpot(limid, cliques, d, clq_ass_to_node);
clpot = init_clpot(limid, cliques, [], clq_ass_to_node);

% collect to root
if 1
  % HUGIN
  seppot = cell(C, C);    % separators are implicitely initialized to 1s
  for n=postorder{di}(1:end-1)
    for p=parents(rooted_jtree{di}, n)
      %clpot{p} = divide_by_pot(clpot{n}, seppot{p,n}); % dividing by 1 is redundant
      seppot{p,n} = marginalize_pot(clpot{n}, separator{p,n});
      clpot{p} = multiply_by_pot(clpot{p}, seppot{p,n});
    end
  end
else
  % Shafer-Shenoy
  msg = cell(C,C);
  for n=postorder{di}(1:end-1)
    for c=children(rooted_jtree{di}, n)
      clpot{n} = multiply_by_pot(clpot{n}, msg{c,n});
    end
    p = parents(rooted_jtree{di}, n);
    %msg{n,p} = marginalize_pot(clpot{n}, cliques{p});
    msg{n,p} = marginalize_pot(clpot{n}, separator{n,p});
  end
  root = clq_ass_to_node(d);
  n=postorder{di}(end);
  assert(n == root);
  for c=children(rooted_jtree{di}, n)
    clpot{n} = multiply_by_pot(clpot{n}, msg{c,n});
  end
end

fam = family(limid.dag, d);
pot = marginalize_pot(clpot{root}, fam);

%%%%%%%
jpot = compute_joint_pot_limid(limid);
pot2 = marginalize_pot(jpot, fam);
assert(approxeq_pot(pot, pot2))
%%%%%%

[policy, score] = extract_policy(pot);

e = limid.equiv_class(d);
limid.CPD{e} = set_params(limid.CPD{e}, 'policy', policy);


text:
clear all; close all; clc
%%
M_dir = GetFishDirectories();
data_masterdir = GetNestedDataDir();
save_masterdir = GetCurrentDataDir();
%%
tiffname = 'C:\Janelia2015\Elavl3-H2BRFP_6dpf_MeanImageOf10Fish.tif';
info = imfinfo(tiffname,'tiff');
nPlanes = length(info);
s1 = info(1).Height;
s2 = info(1).Width;
anat_stack_norm = zeros(s1,s2,nPlanes);
for i=1:nPlanes,
    anat_stack_norm(:,:,i) = imread(tiffname,i);
end
%%
% im = max(anat_stack_norm,[],3);
% out=imNormalize99(im);
% anat_yx_norm = repmat(out,[1 1 3]);
%
% % y-z view
% im = squeeze(max(anat_stack_norm,[],2));
% out=imNormalize99(im);
% anat_yz_norm = repmat(out,[1 1 3]);
%
% % x-z view
% im = squeeze(max(anat_stack_norm,[],1));
% out=imNormalize99(im);
% out = flipud(out'); %%%% empirically necessary...
% anat_zx_norm = repmat(out,[1 1 3]);
%%
range_fish = 2:11;

for i_fish = range_fish,
    tic
    disp(['i_fish = ', num2str(i_fish)]);

    data_dir = fullfile(data_masterdir,['subject_' num2str(i_fish)]);
%     save_dir = fullfile(save_masterdir,['subject_' num2str(i_fish)]);

    load(fullfile(data_dir,'OptionalInfo.mat'));
%     load(fullfile(save_dir,'data.mat'));
%     absIX = data.absIX;
%     absIX = (1:data.numcell_full)';

%     filename = fullfile(save_dir,'TimeSeries.h5');
%     h5create(filename,'/absIX',size(absIX));
%     h5write(filename,'/absIX',absIX);

%     CellXYZ = data.CellXYZ;
%     CellXYZ_ref = data.CellXYZ_ref;
%
%     load(fullfile(save_dir,'data.mat'));
%
%     data.CellXYZ = CellXYZ;
%     data.CellXYZ_ref = CellXYZ_ref;

%%
filename = fullfile('C:\Janelia2015\norm cell coord',['F' num2str(i_fish) '_XYZ_norm.txt']);
delimiter = ' ';
formatSpec = '%f%f%f%s%[^\n\r]';
fileID = fopen(filename,'r');
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'EmptyValue' ,NaN, 'ReturnOnError', false);
fclose(fileID);

VarName1 = dataArray{:, 1};
VarName2 = dataArray{:, 2};
VarName3 = dataArray{:, 3};
FAILED = dataArray{:, 4};

temp = zeros(length(FAILED),1);
for i = 1:length(FAILED),
    temp(i) = length(FAILED{i});
end
IX_inval_norm = find(temp~=0);

Y = round(VarName1/0.798);
X = round(VarName2/0.798);
Z = round(VarName3/2);
CellXYZ_norm = horzcat(X,Y,Z);

    %%

    save(fullfile(data_dir,'OptionalInfo.mat'),'Behavior_raw','CellXYZ_norm','IX_inval','numcell_full','IX_inval_norm');%,'anat_stack_norm'

    toc
end

%%
range_fish = 2:11;
for i_fish = range_fish,
    %% load data
    disp(['load fish ' num2str(i_fish) '...']);
    save_dir = fullfile(save_masterdir,['subject_' num2str(i_fish)]);
    load(fullfile(save_dir,'data_full.mat'),'data'); % struct with many fields
    %%
    data_full = data;
    names = fieldnames(data_full); % cell of strings
    for i = 1:length(names),
        eval([names{i} ' = data_full.' names{i} ';']);
    end

    data_dir = fullfile(data_masterdir,['subject_' num2str(i_fish)]);
    load(fullfile(data_dir,'OptionalInfo.mat'));

    %%
    data = [];
    names = {'periods','timelists_names','stimuluskey_raw','CellXYZ','anat_stack','fpsec',...
        'Behavior_raw','numcell_full','CellXYZ_norm','IX_inval_norm','IX_inval_anat',...
        'anat_yx','anat_yz','anat_zx',...
        'timelists','stim_full','stimAvr','Behavior_full','BehaviorAvr'};
%         'absIX'};
    if length(timelists_names)>1, % M_stimset(i_fish) > 1,
        names = [names,{'stimset'}];
    end

    for i = 1:length(names), % use loop to save variables into fields of 'data'
        eval(['data.',names{i},'=', names{i},';']);
    end

    save(fullfile(save_dir,'data_full.mat'),'data');
end
%%
filename = fullfile(save_dir,'TimeSeries.h5');
fileattrib(filename,'+w');
h5writeatt(filename,'/','absIX',absIX);
h5disp(filename);

%%
% srcFile = fullfile(matlabroot,'toolbox','matlab','demos','example.h5');
copyfile(filename,'myfile.h5','f');
fileattrib('myfile.h5','+w');
h5writeatt('myfile.h5','/','creation_date',datestr(now));

h5writeatt('myfile.h5','/','absIX',absIX(1:10));

%%
filename = fullfile(save_dir,'TimeSeries_full.h5');
%%
tic
CellResp = h5read(filename,'/CellResp');
CellRespZ = h5read(filename,'/CellRespZ');
CellRespAvr = h5read(filename,'/CellRespAvr');
CellRespAvrZ = h5read(filename,'/CellRespAvrZ');
toc
%%
tic
CellResp_ = CellResp(absIX,:);
CellRespZ_ = CellRespZ(absIX,:);
CellRespAvr_ = CellRespAvr(absIX,:);
CellRespAvrZ_ = CellRespAvrZ(absIX,:);
toc

%%
absIX = (1:data.numcell_full)';

filename = fullfile(save_dir,'TimeSeries_full.h5');
h5create(filename,'/CellResp',size(CellResp),'Datatype','single','ChunkSize',[1000 100]);
h5write(filename,'/CellResp',CellResp);

h5create(filename,'/CellRespZ',size(CellRespZ),'Datatype','single','ChunkSize',[1000 100]);
h5write(filename,'/CellRespZ',CellRespZ);

h5create(filename,'/CellRespAvr',size(CellRespAvr),'Datatype','single','ChunkSize',[1000 100]);
h5write(filename,'/CellRespAvr',CellRespAvr);

h5create(filename,'/CellRespAvrZ',size(CellRespAvrZ),'Datatype','single','ChunkSize',[1000 100]);
h5write(filename,'/CellRespAvrZ',CellRespAvrZ);

h5create(filename,'/absIX',size(absIX));
h5write(filename,'/absIX',absIX);

%% need to delete 'absIX' from data_full.mat

author: xiuyechen

repo: FishExplorer


text:
function fh = plot4D(this, varargin)
%plots an MR image per slice
%
%   Y  = MrImage
%   fh = Y.plot4D('ParameterName', ParameterValue)
%
% IN
%   varargin    'ParameterName', 'ParameterValue'-pairs for the following
%               properties:
%               'plotType'          Type of plot that is created
%                                       'montage'   images are plotted as
%                                                   montages of slices or
%                                                   volumes
%                                       'labeledMontage'
%                                                   as montage, but with
%                                                   slice/volume labels
%                                                   (default)
%                                       'overlay'   overlays of images are
%                                                   plotted with different
%                                                   colormaps (e.g. for
%                                                   activation maps, mask
%                                                   visualization)
%                                       'spm'       uses display functions
%                                                   from SPM (spm_display/
%                                                   spm_check_registration)
%                                                   to visualize 3D volumes
%                                                   with header information
%                                                   applied ("world space")
%                                                   Note: if multiple
%                                                   selected volumes are
%                                                   specified,
%                                                   spm_check_registration
%                                                   is used
%                                       'spmInteractive' /'spmi'
%                                                   same as SPM, but keeps
%                                                   temporary nifti files to
%                                                   allow clicking into spm
%                                                   figure
%                                       '3D'/'3d'/'ortho'
%                                                   See also tapas_uniqc_view3d plot3
%                                                   Plots 3 orthogonal
%                                                   sections
%                                                   (with CrossHair) of
%                                                   3D image interactively
%
%               'displayRange'      [1,2] vector for pixel value = black and
%                                                    pixel value = white
%               'signalPart'        for complex data, defines which signal
%                                   part shall be extracted for plotting
%                                       'all'       - take signal as is
%                                                     (default for
%                                                     real-valued data)
%                                       'abs'       - absolute value
%                                                     (default for complex
%                                                     data)
%                                       'phase'     - phase of signal
%                                       'real'      - real part of signal
%                                       'imag'      - imaginary part of
%                                                     signal
%               'plotMode'          transformation of data before plotting
%                                   'linear' (default), 'log'
%               'selectedVolumes'   [1,nVols] vector of selected volumes to
%                                             be displayed
%               'selectedSlices'    [1,nSlices] vector of selected slices to
%                                               be displayed
%                                   choose Inf to display all volumes
%               'sliceDimension'    (default: 3) determines which dimension
%                                   shall be plotted as a slice
%               'rotate90'          default: 0; 0,1,2,3; rotates image
%                                   by multiple of 90 degrees AFTER
%                                   flipping slice dimensions
%               'useSlider'         true or false
%                                   provides interactive slider for
%                                   slices/volumes;
%                                   assumes default:    selectedSlices = Inf
%                                                       selectedVolumes = Inf
%               'fixedWithinFigure' determines what dimension is plotted in
%                                   (subplots of) 1 figure
%                                   'slice(s)'    all slices in 1 figure;
%                                   new figure for each volume
%                                   'volume(s)'   all volumes in 1 figurel
%                                   new figure for each slice
%               'colorMap'          string, any matlab colormap name
%                                   e.g. 'jet', 'gray'
%               'colorBar',         'on' or 'off' (default)
%                                   where applicable, determines whether
%                                   colorbar with displayRange shall be plotted
%                                   in figure;
%               'overlayImages'     (cell of) MrImages that will be
%                                   overlayed
%               'overlayMode'       'edge', 'mask', 'map'
%                                   'edge'  only edges of overlay are
%                                           displayed
%                                   'mask'  every non-zero voxel is
%                                           displayed (different colors for
%                                           different integer values, i.e.
%                                           clusters'
%                                   'map'   thresholded map in one colormap
%                                           is displayed (e.g. spmF/T-maps)
%                                           thresholds from
%                                           overlayThreshold
%               'overlayThreshold'  [minimumThreshold, maximumThreshold]
%                                   tresholds for overlayMode 'map'
%                                   default: [-Inf, Inf] = [minValue, maxValue]
%                                   everything below minValue will not be
%                                   displayed;
%                                   everything above maxValue
%                                   will have brightest color
%               'overlayAlpha'      transparency value of overlays
%                                   (0 = transparent; 1 = opaque; default: 0.1)
%
%
% OUT
%   fh          [nFigures,1] vector of figure handles
%
% EXAMPLE
%
%   Y.plot4D('selectedVolumes', [6:10])
%   Y.plot4D('displayRange', [0 1000])
%   Y.plot4D('useSlider', true, 'selectedVolumes', Inf);
%
%   See also

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-05-21
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% for complex data, plot absolute value per default
if isreal(this)
    defaults.signalPart         = 'all';
else
    defaults.signalPart         = 'abs';
end

defaults.plotType               = 'labeledMontage';
defaults.plotMode               = 'linear';
defaults.selectedVolumes        = 1;
defaults.selectedX              = Inf;
defaults.selectedY              = Inf;
defaults.selectedSlices         = Inf;
defaults.sliceDimension         = 3;
defaults.rotate90               = 0;
defaults.displayRange           = [];
defaults.useSlider              = false;
defaults.fixedWithinFigure      = 'volume';
defaults.colorMap               = 'gray';
defaults.colorBar               = 'off';
defaults.overlayImages          = {};
defaults.overlayMode            = 'mask';
defaults.overlayThreshold       = [];
defaults.overlayAlpha           = 0.1;

args = tapas_uniqc_propval(varargin, defaults);
tapas_uniqc_strip_fields(args);

doPlotColorBar = strcmpi(colorBar, 'on');
doPlotOverlays = any(strcmpi(plotType, {'overlay', 'overlays'})) || ...
    ~isempty(overlayImages);

% slider enables output of all Slices and Volumes per default, strip data
% again under this assumption, if slider is used for display
if useSlider
    defaults.selectedVolumes = Inf;
    defaults.selectedSlices = Inf;
    args = tapas_uniqc_propval(varargin, defaults);
    tapas_uniqc_strip_fields(args);
end

% Assemble parameters for data extraction into one structure
argsExtract = struct('sliceDimension', sliceDimension, ...
    'selectedX', selectedX, 'selectedY', selectedY, ...
    'selectedSlices', selectedSlices, 'selectedVolumes', selectedVolumes, ...
    'plotMode', plotMode, 'rotate90', rotate90, 'signalPart', signalPart);

if isempty(this.data)
    error('tapas:uniqc:MrImage:plot4D:EmptyDataMatrix', ...
        'Data matrix empty for MrImage-object %s', this.name);
end

% retrieve plot data without actually plotting...
if doPlotOverlays
    argsOverlays = struct('sliceDimension', sliceDimension, ...
        'selectedSlices', selectedSlices, 'selectedVolumes', selectedVolumes, ...
        'plotMode', plotMode, 'rotate90', rotate90, 'signalPart', signalPart, ...
        'overlayMode', overlayMode, 'overlayThreshold',  overlayThreshold, ...
        'doPlot', true);

    [fh, dataPlot] = this.plot_overlays(overlayImages, argsOverlays);
    return
else
    if isempty(displayRange)
        [dataPlot, displayRange, resolution_mm] = this.extract_plot4D_data(argsExtract);
    else
        [dataPlot, ~, resolution_mm] = this.extract_plot4D_data(argsExtract);
    end
end

nVolumes = size(dataPlot,4);
nSlices = size(dataPlot,3);

if isinf(selectedSlices)
    selectedSlices = 1:nSlices;
end

if isinf(selectedVolumes)
    selectedVolumes = 1:nVolumes;
end

% slider view
if useSlider
    % useSlider is not a plotType, since it shall be combined with all
    % plot-types (overlays, montages) in a later version of this code

    % tapas_uniqc_slider4d(dataPlot, @(varargin) ...
    %      tapas_uniqc_plot_abs_image(varargin{:}, colorMap), ...
    %     nSlices);

    tapas_uniqc_slider4d(dataPlot, @(Y,iDynSli, fh, yMin, yMax) ...
        tapas_uniqc_plot_abs_image(Y,iDynSli, fh, yMin, yMax, colorMap, colorBar), ...
        nSlices, displayRange(1), displayRange(2), this.name);

    % to also plot phase:
    %    tapas_uniqc_slider4d(dataPlot, @tapas_uniqc_plot_image_diagnostics, ...
    %        nSlices);

else
    switch lower(plotType)
        case {'3d', 'ortho'}
            this.plot3d(argsExtract);
        case {'spm', 'spminteractive', 'spmi'}
            % calls spm_image-function (for single volume) or
            % spm_check_registration (multiple volumes)

            % get current filename, make sure it is nifti-format
            fileName = this.parameters.save.fileName;
            fileNameNifti = fullfile(this.parameters.save.path, ...
                regexprep(fileName, '\..*$', '\.nii'));

            % create nifti file, if not existing
            % TODO: how about saved objects with other file names
            if ~exist(fileNameNifti, 'file')
                this.save(fileNameNifti);
            end

            % select Volumes
            fileNameVolArray = tapas_uniqc_get_vol_filenames(fileNameNifti);

            % display image
            if numel(selectedVolumes) > 1
                spm_check_registration( ...
                    fileNameVolArray{selectedVolumes});
            else
                spm_image('Display', fileNameVolArray{selectedVolumes});
            end

            % delete temporary files for display
            if strcmpi(this.parameters.save.keepCreatedFiles, 'none')

                switch lower(plotType)
                    case {'spminteractive', 'spmi'}
                        input('Press Enter to leave interactive mode');
                end

                delete(fileNameNifti);
                [stat, mess, id] = rmdir(this.parameters.save.path);
            end

            % restore original file name
            this.parameters.save.fileName = fileName;

        case {'montage', 'labeledmontage'}

            if strcmpi(plotType, 'labeledMontage')

                stringLabelSlices = cellfun(@(x) num2str(x), ...
                    num2cell(selectedSlices), 'UniformOutput', false);
                stringLabelVolumes = cellfun(@(x) num2str(x), ...
                    num2cell(selectedVolumes), 'UniformOutput', false);
            else
                stringLabelSlices = [];
                stringLabelVolumes = [];
            end

            switch lower(fixedWithinFigure);
                case {'volume', 'volumes'}
                    for iVol = 1:nVolumes

                        nFrames = size(dataPlot, 3);

                        stringTitle = sprintf('%s - volume %d', this.name, ...
                            selectedVolumes(iVol));
                        fh(iVol,1) = figure('Name', stringTitle, 'WindowStyle', 'docked');

                        tapas_uniqc_labeled_montage(permute(dataPlot(:,:,:,iVol), [1, 2, 4, 3]), ...
                            'DisplayRange', displayRange, ...
                            'LabelsIndices', stringLabelSlices);

                        set(gca, 'DataAspectRatio', abs([resolution_mm(1) resolution_mm(2), 1]));

                        title(tapas_uniqc_str2label(stringTitle));
                        if doPlotColorBar
                            colorbar;
                        end
                        colormap(colorMap);

                    end

                case {'slice', 'slices'}

                    for iSlice = 1:nSlices
                        stringTitle = sprintf('%s - slice %d', this.name, ...
                            selectedSlices(iSlice));
                        fh(iSlice,1) = figure('Name', stringTitle, 'WindowStyle', 'docked');
                        tapas_uniqc_labeled_montage(dataPlot(:,:,iSlice,:), ...
                            'DisplayRange', displayRange, ...
                            'LabelsIndices', stringLabelVolumes);

                        set(gca, 'DataAspectRatio', [resolution_mm(1) resolution_mm(2), 1]);
                        title(tapas_uniqc_str2label(stringTitle));
                        if doPlotColorBar
                            colorbar;
                        end
                        colormap(colorMap);
                    end

            end % fixedWithinFigure
    end % plotType
end % use Slider
end

author: translationalneuromodeling

repo: tapas


text:
function top = startofPOWCone(K)

top = 1 + K.f + K.l + sum(K.q) + 3*K.e;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function model = ivmRun(XTrain, yTrain, options, kernelTieStructure, noiseTieStructure)

% IVMRUN Run the IVM on a given data set.
% FORMAT
% DESC takes an input data matrix, a target data matrix and an options
% structure (as defined by ivmOptions). It then creates and
% optimises an IVM using ivmOptimise.
% ARG X : the input training data.
% ARG y : the target training data.
% ARG options : options structure (as defined by ivmOptions).
%
% FORMAT
% DESC is as above but includes the facility to 'tie' some of the
% kernel an noise parameters together.
% ARG X : the input training data.
% ARG y : the target training data.
% ARG options : options structure (as defined by ivmOptions).
% ARG kernelTie : structure which speficies which kernel parameters
% are forced to be equal to each other (see cmpndTieParameters).
% ARG noiseTie : structure which speficies which noise parameters
% are forced to be equal to each other (see cmpndTieParameters).
%
% SEEALSO : ivmCreate, ivmOptimise, ivmOptions, cmpndTieParameters
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006, 2007

% IVM

% Intitalise IVM
model = ivmCreate(size(XTrain, 1), size(yTrain, 2), XTrain, yTrain, options);

if nargin > 8 & ~isempty(noiseTieStructure);
  % Some of the noise parameters are constrained equal to each other
  model.noise = cmpndTieParameters(model.noise, noiseTieStructure);
end
if nargin > 7 & ~isempty(kernelTieStructure);
  % Some of the kernel parameters are constrained to equal each other.
  model.kern = cmpndTieParameters(model.kern, kernelTieStructure);
end

% Optimise the parameters
model = ivmOptimise(model, options);

% Select data-points in an IVM with the given parameters.
model = ivmOptimiseIvm(model, options.display);


text:
%% OPT_POOL uses the MATLABPOOL command to run the QUAD code.
%
%  Discussion:
%
%    Output printed by the function appears directly on the screen.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2010
%
%  Author:
%
%    John Burkardt
%
  n = 16;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'OPT_POOL\n' );
  fprintf ( 1, '  Solve optimization problem with discretization parameter N = %d\n', n );

  matlabpool open local 4

  [ z_star, PAR, exit ] = opt_fun ( @zermelo, n );

  matlabpool close
%
%  On successful conclusion, display the plot.
%
  if ( 0 <= exit && isfield ( PAR, 'plot' ) )
    feval ( PAR.plot, z_star, PAR )
  end

author: johannesgerer

repo: jburkardt-m


text:
function y = GUSFT_simple(x,shift,boxlen,center,w);
% GUSFT_simple -- Gram matrix of the USFT
% Usage:
%   y = GUSFT_simple(x,shift,boxlen,center,w)
% Inputs:
%   x	    vector of length n
%   shift   vector of shifts
%   boxlen  half-length of the window associated with shift
%   center  boolean variable
%   w       window of size 2*boxlen
% Outputs:
%   y      vector le length 2*n
% Description:
%  Evaluates the A'A where A is the nonuniform FT at the points omega_k =
%  shift + k, -boxlen <= k <boxlen
%  See Also
%    USFT_simple, USFFT, MakeFourierDiagonal, GUSFT_Toeplitz
%
% By Emmanuel candes, 2003-2004

if nargin < 5,
  w = ones(1,2*boxlen);
end

if nargin < 4,
  center = 0;
end

       tx  =  USFT_simple(x,shift,boxlen,center,w);
       y   =  Adj_USFT_simple(tx,shift,boxlen,center,w);

author: thomaskuestner


text:
%% Unit Cube
[node,elem] = cubemesh([-1,1,-1,1,-1,1],1);
subplot(1,2,1); showmesh3(node,elem);

%% Lshape domain
% mesh
[node,elem] = cubemesh([-1,1,-1,1,-1,1],1);
[node,elem] = delmesh(node,elem,'x<0 & y<0');
bdFlag = setboundary3(node,elem,'Dirichlet');
% for adaptive mesh refinement, special ordering of elem and HB is needed
[elem,bdFlag,HB] = label3(node,elem,'all',bdFlag);
showmesh3(node,elem);

repo: ifem


text:
%*******************************************************************
%Summary: This file runs the specified model and logs the data from the listed signals
%and plots the results.
%
%In this file you specify the name of the signal
%that you want to log.  You can also specify a decimation for the logged
%signals.  This models uses a variation of the Simulink Demo model named
%'sldemo_househeat'
% JAD - 1/6/06
%*********************************************************************

%**********Prepare workspace************
    clear
    clc
    close all
    name_mod='sldemo_househeat1'; %Name of model to simulate
%***************end***************************

    % Simulation Time Calculationss *******

    Sub_Sample_Time=3; % Log every third value

%Structure file that will contain the name of the signals to log
%Comment all names to disable signal logging

    log_struct{1}='Set_Temp';
    log_struct{2}='blower_cmd';
    log_struct{3}='Tout';
    log_struct{4}='Temp_in';
    log_struct{5}='Heat_Cost';
    log_struct{6}='Tout2';
    log_struct{7}='time';  %The time value can also be obtained from the logged Structure

    log_struct{8}='House/Sub_Temp';
%     log_struct{9}='';
%     log_struct{10}='';
%     log_struct{11}='';
%     log_struct{12}='';
%     log_struct{13}='';
%
%     log_struct{14}='';
%     log_struct{15}='';
%     log_struct{16}='';
%     log_struct{17}='';
%     log_struct{18}='';
%     log_struct{19}='';

    %Turn on logging for specified signals in the model
    if exist('log_struct')
                   %log_signals(name of model,name_of_log_var, enable_decimation, decimation_time)
        signal_name=log_signals(name_mod,log_struct,'off',Sub_Sample_Time);
    end

    % ********* Run Sim  ***********************************************
    tic
     sim(name_mod)
    toc
%***************end***************************

%*************Save logged values in the workspace************************
if exist('log_struct')

    for i=1:length(log_struct)
        index=findstr(log_struct{i},'/');

        if isempty(index)
           temp=[signal_name{i},'=logsout.',signal_name{i},'.Data;'];
           eval(temp);

        elseif (length(index)>1) % If this signal is in a subsystem
            sys_name=log_struct{i}(1:index(1)-1); %Grab the name of the system
            sub_sys_name=log_struct{i}(index(1)+1:index(2)-1); %Grab the name of the subsystem
            temp=[signal_name{i},'=logsout.(''',sys_name,''').',sub_sys_name,'.',signal_name{i},'.Data;'];
            eval(temp);
        else  %This signal is not in a subsystem
           temp1=log_struct{i}(1:index(1)-1); %Grab the name of the system
           temp=[signal_name{i},'=logsout.(''',temp1,''').',signal_name{i},'.Data;'];
           eval(temp);
        end
    end
end
%***************end***************************

%Close model and do not save model changes
%       close_system(name_mod,0);
 %Plot results

    figure(1)
    subplot(2,1,1)
    plot(time,Tout,time,Tout2)
    subplot(2,1,2)
    plot(Heat_Cost,time)

repo: mcbench-benchmarks


text:
filename='CantileverSquare';
ptype = 'MACRO';
initial_case = 'full';
%cost = {'compliance'};
cost = {'stressNorm'};
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'P1';

optimizer = 'DualNestedInPrimal';
optimizerUnconstrained = 'PROJECTED GRADIENT';
designVariable = 'MicroParams';
homegenizedVariablesComputer = 'ByVademecum';
vademecumFileName = 'Rectangle';%'SmoothRectangle';
%vademecumFileName = 'SmoothRectangle';

nsteps = 1;
maxiter = 500;

Vfrac_initial = 0.3;
Vfrac_final = 0.3;

optimality_initial = 1e-5;
optimality_final = 1e-5;
constr_initial = 1e-3;
constr_final = 1e-3;

printing = false;
monitoring_interval = 1;

ub = 0.989;
lb = 0.011;
rate = 0.5;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% define global paths
global CBTDIR
global GUROBI_PATH
global ILOG_CPLEX_PATH
global TOMLAB_PATH

c = clock;
fprintf('\n\n      _____   _____   _____   _____     _____     |\n');
fprintf('     /  ___| /  _  \\ |  _  \\ |  _  \\   / ___ \\    |   COnstraint-Based Reconstruction and Analysis\n');
fprintf(['     | |     | | | | | |_| | | |_| |  | |___| |   |   The COBRA Toolbox - ' num2str(c(1)) '\n']);
fprintf('     | |     | | | | |  _  { |  _  /  |  ___  |   |\n');
fprintf('     | |___  | |_| | | |_| | | | \\ \\  | |   | |   |   Test Suite\n');
fprintf(['     \\_____| \\_____/ |_____/ |_|  \\_\\ |_|   |_|   |\n']);
fprintf('                                                  | \n\n');

% request explicitly from the user to launch test suite locally
if contains(getenv('HOME'), 'vmhadmin') || contains(getenv('HOME'), 'jenkins')
    % Running in CI environment
    fprintf('Running test in Jenkins/CI environment\n');

    % on the CI, always reset the path to make absolutely sure, that we test
    % the current version
    restoredefaultpath;
    launchTestSuite = true;
else
    reply = '';
    while isempty(reply)
        reply = input([' -> Do you want to launch the test suite locally? Time estimate: more than 60 minutes Y/N: '], 's');
    end

    if strcmpi(reply, 'y') || strcmpi(reply, 'yes')
        launchTestSuite = true;
    else
        launchTestSuite = false;
    end
end

% save the current folder
origDir = pwd;

% if the location of initCobraToolbox is not yet known
if isempty(which('initCobraToolbox.m'))
    % move back to the root of the repository
    cd([fileparts(which('testAll.m')) filesep '..'])

    % assign the path
    CBTDIR = pwd;
else
    CBTDIR = fileparts(which('initCobraToolbox.m'));
    cd(CBTDIR);
end

% change to the root folder of The COBRA TOolbox
cd(CBTDIR);

% run the official initialisation script
if launchTestSuite

    if ~isempty(getenv('MOCOV_PATH')) && ~isempty(getenv('JSONLAB_PATH'))
        addpath(genpath(getenv('MOCOV_PATH')))
        addpath(genpath(getenv('JSONLAB_PATH')))
        COVERAGE = true;
        fprintf('MoCov and JsonLab are on path, coverage will be computed.\n')
    else
        COVERAGE = false;
    end

    initCobraToolbox;

    % initialize the cleanup
    currentDir = cd('test');
    testDirContent = getFilesInDir('type', 'all');  % Get all currently present files in the folder.
    testDirPath = pwd;
    cd(currentDir);
%{
    if ~isempty(strfind(getenv('HOME'), 'jenkins')) || ~isempty(strfind(getenv('USERPROFILE'), 'jenkins'))
        WAITBAR_TYPE = 0;

        % add CellNetAnalyzer for testing purposes
        addpath(genpath(getenv('CNA_PATH')));

        % check the CNA installation
        checkCNAinstallation(0);
    else
        WAITBAR_TYPE = 1;
    end
%}
    if verLessThan('matlab', '8.2')
        error('The testsuite of The COBRA Toolbox can only be run with MATLAB R2014b+.')
    end

    % define a success exit code
    exit_code = 0;

    % enable profiler
    profile on;

    if COVERAGE
        % Get the ignored Files from gitIgnore
        % only retain the lines that end with .txt and .m and
        % are not comments and point to files in the /src folder
        ignoredPatterns = {'^.{0,3}$', ...  % Is smaller than four.
                        ['^[^s][^r][^c][^' regexptranslate('escape', filesep) ']']};  % does not start with src/
        filterPatterns = {'\.txt$', '\.m$'};  % Is either a .m file or a .txt file.
        ignoreFiles = getIgnoredFiles(ignoredPatterns, filterPatterns);

        % check the code quality
        listFiles = getFilesInDir('type', 'tracked', 'restrictToPattern', '^.*\.m$', 'checkSubFolders', true);

        % count the number of failed code quality checks per file
        nMsgs = 0;
        nCodeLines = 0;
        nEmptyLines = 0;
        nCommentLines = 0;

        for i = 1:length(listFiles)
            nMsgs = nMsgs + length(checkcode(listFiles{i}));
            fid = fopen(listFiles{i});

            while ~feof(fid)
                lineOfFile = strtrim(char(fgetl(fid)));
                if length(lineOfFile) > 0 && length(strfind(lineOfFile(1), '%')) ~= 1  ...
                    && length(strfind(lineOfFile, 'end')) ~= 1 && length(strfind(lineOfFile, 'otherwise')) ~= 1 ...
                    && length(strfind(lineOfFile, 'switch')) ~= 1 && length(strfind(lineOfFile, 'else')) ~= 1  ...
                    && length(strfind(lineOfFile, 'case')) ~= 1 && length(strfind(lineOfFile, 'function')) ~= 1
                        nCodeLines = nCodeLines + 1;

                elseif length(lineOfFile) == 0
                    nEmptyLines = nEmptyLines + 1;

                elseif length(strfind(lineOfFile(1), '%')) == 1
                    nCommentLines = nCommentLines + 1;
                end
            end
            fclose(fid);
        end

        % average number of messages per codeLines
        avMsgsPerc = floor(nMsgs / nCodeLines * 100);

        grades = {'A', 'B', 'C', 'D', 'E', 'F'};
        intervals = [0, 3;
                    3, 6;
                    6, 9;
                    9, 12;
                    12, 15;
                    15, 100];

        grade = 'F';
        for i = 1:length(intervals)
            if avMsgsPerc >= intervals(i, 1) && avMsgsPerc < intervals(i, 2)
                grade = grades{i};
            end
        end

        fprintf('\n\n -> The code grade is %s (%1.2f%%).\n\n', grade, avMsgsPerc);

        % set the new badge
        if ~isempty(strfind(getenv('HOME'), 'jenkins'))
            coverageBadgePath = [getenv('ARTENOLIS_DATA_PATH') filesep 'cobratoolbox' filesep 'codegrade' filesep];
            system(['cp ' coverageBadgePath 'codegrade-', grade, '.svg '  coverageBadgePath 'codegrade.svg']);
        end
    end
end

fprintf ('force-setting launchTestSuite = true\n');
launchTestSuite = true;
fprintf ('launchTestSuite=%s\n', string(launchTestSuite));

try
    if launchTestSuite
        % save the userpath
        originalUserPath = path;

        % run the tests in the subfolder verifiedTests/ recursively
        [result, resultTable] = runTestSuite();

        sumSkipped = sum(resultTable.Skipped);
        sumFailed = sum(resultTable.Failed);

        fprintf(['\n > ', num2str(sumFailed), ' tests failed. ', num2str(sumSkipped), ' tests were skipped due to missing requirements.\n\n']);

        % count the number of covered lines of code
        if COVERAGE
            % write coverage based on profile('info')
            fprintf('Running MoCov ... \n')
            mocov('-cover', 'src', ...
                '-profile_info', ...
                '-cover_json_file', 'coverage.json', ...
                '-cover_html_dir', 'coverage_html', ...
                '-cover_method', 'profile', ...
                '-verbose');

            % load the coverage file
            data = loadjson('coverage.json', 'SimplifyCell', 1);

            sf = data.source_files;
            clFiles = zeros(length(sf), 1);
            tlFiles = zeros(length(sf), 1);

            for i = 1:length(sf)
                clFiles(i) = nnz(sf(i).coverage);
                tlFiles(i) = length(sf(i).coverage);
            end

            % average the values for each file
            cl = sum(clFiles);
            tl = sum(tlFiles);

            % print out the coverage
            fprintf('Covered Lines: %i, Total Lines: %i, Coverage: %f%%.\n', cl, tl, cl / tl * 100);
        end

        % print out a summary table
        resultTable

        % Print some information on failed and skipped tests.
        skippedTests = find(resultTable.Skipped);
        if sum(skippedTests > 0)
            fprintf('The following tests were skipped:\n%s\n\n', strjoin(resultTable.TestName(skippedTests), '\n'));
            fprintf('The reasons were as follows:\n')
            for i = 1:numel(skippedTests)
                fprintf('------------------------------------------------\n')
                fprintf('%s:\n', resultTable.TestName{skippedTests(i)});
                fprintf('%s\n', resultTable.Details{skippedTests(i)});
                fprintf('------------------------------------------------\n')
            end
            fprintf('\n\n')
        end

        failedTests = find(resultTable.Failed & ~resultTable.Skipped);
        if sum(failedTests > 0)
            fprintf('The following tests failed:\n%s\n\n', strjoin(resultTable.TestName(failedTests), '\n'));
            fprintf('The reasons were as follows:\n')
            for i = 1:numel(failedTests)
                fprintf('------------------------------------------------\n')
                fprintf('%s:\n', resultTable.TestName{failedTests(i)});
                trace = result(failedTests(i)).Error.getReport();
                tracePerLine = strsplit(trace, '\n');
                testSuitePosition = find(cellfun(@(x) ~isempty(strfind(x, 'runTestSuite')), tracePerLine));
                trace = sprintf(strjoin(tracePerLine(1:(testSuitePosition - 7)), '\n'));  % Remove the testSuiteTrace.
                fprintf('%s\n', trace);
                fprintf('------------------------------------------------\n')
            end
            fprintf('\n\n')
        end

        % restore the original path
        restoredefaultpath;
        addpath(originalUserPath);

        if sumFailed > 0
            exit_code = 1;
        end

        fprintf(['\n > The exit code is ', num2str(exit_code), '.\n\n']);

        % clean up temporary files.
        removeTempFiles(testDirPath, testDirContent);

        % ensure that we ALWAYS call exit
        if ~isempty(strfind(getenv('HOME'), 'jenkins')) || ~isempty(strfind(getenv('USERPROFILE'), 'jenkins'))
            %exit(exit_code);
        end
    end
catch ME

    fprintf('exception %s while running test: %s\n', ME.identifier, ME.message);

    % Also clean up temporary files in case of an error.
    removeTempFiles(testDirPath, testDirContent);
    if ~isempty(strfind(getenv('HOME'), 'jenkins')) || ~isempty(strfind(getenv('USERPROFILE'), 'jenkins'))
        % only exit on jenkins.
        exit(1);
    else
        % switch back to the folder we were in and rethrow the error
        cd(origDir);
        rethrow(ME);
    end
end

% switch back to the original directory
cd(origDir)
if contains(getenv('HOME'), 'vmhadmin') || contains(getenv('HOME'), 'jenkins')
    % Running in CI environment
    fprintf('Running test in Jenkins/CI environment\n');
    % explicit 'exit' required for R2018b in non-interactive mode to avoid SEGV near end of test
    exit
end

author: opencobra

repo: cobratoolbox


text:
classdef STFTNode < GraphNode
    properties
        fft_len=512;
        frame_len=400;
        frame_shift=160;
        win_type = 'hamming';
        removeDC = 1;
        doDithering = 1;
    end

    methods
        function obj = STFTNode(dimOut)
            obj = obj@GraphNode('STFT',dimOut);
        end

        function obj = forward(obj,prev_layers)
            input = prev_layers{1}.a;
            [nCh,D2,T,N] = size(input);
            assert(D2==1, sprintf('%s:forward, Error: dimension 2 (%d) is not equal to 1', obj.name, D2));
            input = permute(input, [1 3 4 2]);

            precision = class(gather(input(1)));
            useGPU = IsInGPU(input);

            if N==1
                fft_x = sfft_multi(input',obj.frame_len,obj.frame_shift,obj.fft_len, obj.win_type,obj.removeDC, useGPU, obj.doDithering);
                fft_x = fft_x(1:obj.fft_len/2+1,:,:);
                [d1,d2,d3] = size(fft_x);
                fft_x = reshape(fft_x, d1*d2, d3);
                obj.mask = [];
            else
                [mask, variableLength] = CheckTrajectoryLength(input);
                input2 = PadShortTrajectory(input, mask, 0);
                input2 = reshape(permute(input2, [2 1 3]), T, nCh*N);

                fft_x = sfft_multi(input2,obj.frame_len,obj.frame_shift,obj.fft_len, obj.win_type,obj.removeDC, useGPU, obj.doDithering);
                fft_x = fft_x(1:obj.fft_len/2+1,:,:);
                fft_x = permute(fft_x, [1 3 2]);
                [nBin,nFr,d3] = size(fft_x);
                fft_x = reshape(fft_x, nBin, nFr, nCh, N);
                fft_x = permute(fft_x, [1 3 2 4]);
                fft_x = reshape(fft_x, nBin*nCh, nFr, N);

                % now build a mask for spectrogram
                nSampleChannel = sum(mask==0);
                nSampleChannel = gather(nSampleChannel);
                if useGPU
                    obj.mask = gpuArray.zeros(nFr, N, precision);
                else
                    obj.mask = zeros(nFr, N, precision);
                end
                for i=1:N
                    nFrChannel = enframe_decide_frame_number(nSampleChannel(i), obj.frame_len, obj.frame_shift, 0);
                    obj.mask(nFrChannel+1:end,i) = 1;
                end
            end
            obj.a = permute(fft_x, [1 4 2 3]);
            obj = forward@GraphNode(obj, prev_layers);
        end

        function obj = backward(obj,future_layers, prev_layers)
            % to be implemented
        end

    end

end

author: singaxiong

repo: SignalGraph


text:
function figOffscreenCheck(fig);
% figOffscreenCheck([fig]): check that a figure
% doesn't get moved off the edges of the screen.
%
% This function checks the figure fig [default gcf]
% and, if it's running off the edge of the screen, moves / reshapes
% it to fit within the screen.
%
% ras, 07/05/05.
if ~exist('fig','var') | isempty(fig), fig = gcf;     end

% record the existing units of the figure, then set to normalized
% for now:
exUnits = get(fig,'Units');
set(fig,'Units','normalized');
pos = get(fig,'Position');

% first, check if either size argument is larger
% than the screen size (normalized value > 1):
if pos(3)>1, pos(3)=1; end
if pos(4)>1, pos(4)=1; end

% check if right edge out of bounds
if pos(1)+pos(3)>1,  pos(1) = 1-pos(3);    end

% check if left edge out of bounds
if pos(1)<0,    pos(1) = 0;    end

% check if upper edge out of bounds
% (allow extra 0.05 for toolbars which matlab sometimes adds)
if pos(2)+pos(4)>0.95,  pos(2) = 0.95-pos(4);    end

% check if lower edge out of bounds
if pos(2)<0,    pos(2) = 0;    end

% set new position
set(fig,'Position',pos);

% restore old unit convention
set(fig,'Units',exUnits);

return

repo: vistasoft


text:
function s = gsp_jtv_frame_synthesis(F,c)
%GSP_JTV_FRAME_SYNTHESIS Synthesis operator in time-vertex domain
%   Usage:  s = gsp_jtv_frame_synthesis(F,c,param);
%
%   Input parameters:
%         F          : Frame matrix (vertex_loc x time_loc x vertex x time)
%         c          : Coefficients matrix
%   Output parameters:
%         s          : Time-Vertex signal
%
%   This function compute the synthesis operator using frame matrix

% Author :  Francesco Grassi
% Date   : July 2016

[N,lag,~,T]=size(F);
s = reshape(reshape(F,N*lag,[])'*c(:),N,[]);


text:
%IM_THRESHOLD Fixed mapping thresholding images (DIP_Image)
%
%	B = IM_THRESHOLD(A,TYPE,PAR,INV)
%	B = A*IM_THRESHOLD([],TYPE,PAR,INV)
%	B = A*IM_THRESHOLD(TYPE,PAR,INV)
%
% INPUT
%   A        Dataset with object images (possibly multi-band)
%   TYPE     Type of procedure, see below
%   PAR      Related parameter
%   INV      If INV = 1, result inverted, default INV = 0.
%
% OUTPUT
%   B        Dataset with thresholded images
%
% DESCRIPTION
% The following procedures are supported (TYPE)
% 'isodata':      Thresholding using the Isodata algorithm
% 'otsu':         See IM2BW
% 'triangle':    *Thresholding using chord method
%                 (a.k.a. skewed bi-modality, maximum distance to triangle)
%                 by Zack, Rogers and Latt (1973).
% 'background':  *Thresholding using unimodal background-symmetry method.
% 'fixed':        Thresholding at a fixed value.
% 'double':      *Thresholding between two fixed values.
% 'volume':      *Thresholding to obtain a given volume fraction.
% 'hysteresis':  *From the binary image (in>low) only those regions are
%                 selected for which at least one pixel is (in>high)
%
% The following parameters are related to these procedures (PAR):
% ('background'): Distance to the peak where we cut-off, in
%                 terms of the half-width at half the maximum.
%                 Inf selects the default value, which is 2.
% ('fixed'):      Threshold value. Inf means halfway between
%                 minimum and maximum value.
% ('double'):     Two threshold values. Inf means min+[1/3,2/3]*(max-min).
% ('volume'):     Parameter = the volume fraction (Inf means 0.5)
% ('hysteresis'): Two values: low, high (see above)
%
% * routine needs DIP_IMAGE
%
% This routine is for smart thresholding. Simple thresholding can also be
% performed by B = A > PAR, which is the same as using the type 'fixed'.
%
% SEE ALSO (PRTools Guide)
% DATASETS, DATAFILES, DIP_IMAGE, THRESHOLD

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function b = im_threshold(varargin)

	argin = shiftargin(varargin,'char');
  argin = setdefaults(argin,[],'isodata',inf,0);
  if mapping_task(argin,'definition')
    b = define_mapping(argin,'fixed');
    b = setname(b,'Image threshold');
  else
    [a,type,par,inv] = deal(argin{:});
    if isa(a,'prdataset') % allows datafiles too
      isobjim(a);
      b = filtim(a,mfilename,{type,par,inv});
      b = setfeatsize(b,getfeatsize(a));
    elseif isa(a,'double') | isa(a,'dip_image') % here we have a single image
      switch lower(type)
        case('isodata')
          mina = min(a(:));  maxa = max(a(:));
          b = 255*(a-mina)/(maxa-mina);
          t0 = 128;
          for j=1:25
            U = b >  t0;
            L = b <= t0;
            n2 = mean(b(U)); n1 = mean(b(L));
            t1 = (n1+n2)/2;
            if t1 == t0;
              break;
            else
              t0 = t1;
            end
          end
          b = b > t0;
        case('fixed')
          b = a > par;
        case('otsu')
          b = im2bw(a);
        case {'triangle','background','double','volume','hysteresis'}
          checktoolbox('dipimage');
          a = 1.0*dip_image(a);
          b = threshold(a,type,par);
        otherwise
          error 'Illegal thresholding type'
      end
      if inv
        b = 1-b;
      end
    end
  end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [N, D] = ratpolyfit(x,y,kn,kd)
% Ratpolyfit Rational Polynomial Fitting
%
%  This programs finds 2 polynomials N(x) and D(x),
%  of user given order kn and kd respectively,
%  such that N(xi)/D(xi) ~= y(xi) in a least squares sense.
%
%usage: [N, D] = ratpolyfit(x, y, kn, kd)
%
%note: kn and kd must be large enough to get a good fit.
%      usually, kn=kd gives good results
%
%note: If you "overfit" the data, then you will usually have pole-zero
%      cancellations and/or poles and zeros with a very large magnitude.
%      If that happens, then reduce the values of kn and/or kd
%
%note: Often, if you have a good fit, you will find that your polynomials
%      have roots where the real function has zeros and poles.
%
%note: Polynomial curve fitting becomes ill conditioned
%      as the range of x increases and as kn and kd increase
%
%note: If you think that your function goes to infinity at some x value
%      then make sure y(xi) is set equal to Inf at that point.
%      The program will compensate for all +/- Inf values
%
%For example, here's a rational polynomial approximation to the Gamma function:
%   x=[-5:1/32:5]'; y=gamma(x);
%   [N, D]=ratpolyfit(x,y,10,10);
%   figure(1); plot(roots(N),'ob'); hold on; plot(roots(D),'xr'); grid on
%   yy=polyval(N,x)./polyval(D,x);
%   figure(2);plot(x,y,'b', x,yy,'dr'); grid on; axis([min(x) max(x) -25 25]);
%
%other demos are in the text listing at the end of this program
%
%tested under version 6.5 (R13)
%
%see also: polyfit, padefit, polyval, vander
%

%Paul Godfrey
%pgodfrey@conexant.com
%May 31, 2006

%default length if none given
if exist('kn','var')
   kn=round(real(kn));
   if kn<0, kn=0; end
else
   kn=5;
end
if exist('kd','var')
   kd=round(real(kd));
   if kd<0, kd=0;end
else
   kd=5;
end

x=x(:);
y=y(:);

% we must remove +/- Inf values from y first
% and insert those as separate poles at the end
p=find(~isfinite(y)); % find Nan 0/0 and Inf a/0 values
% or find abs(y) values > than some really big number

dinf=[];
while length(p)>0
   y=y.*(x-x(p(1))); %adjust remaining y values
   y(p(1))=[]; % remove bad y value, now a NaN
   dinf=[dinf; x(p(1))]; % remember where pole was
   x(p(1))=[]; % now remove that x value too
   if kd>0, kd=kd-1; end; %reduce expected order of den
   p=find(~isfinite(y)); % have all Inf values been removed yet?
end

yy=length(y);
an=ones(yy,kn+1);
for k=kn:-1:1
    an(:,k)=x.*an(:,k+1); % form vandermonde matrix
end

ad=ones(yy,kd+1);
for k=kd:-1:1
    ad(:,k)=x.*ad(:,k+1);
end
for k=1:yy
    ad(k,:)=y(k)*ad(k,:);
end

% A is basically N-y*D
A=[an -ad]; % LS solution is in the null space of A

[u,s,v]=svd(A); % null space is in the cols of V
ND=v(:,end); % use the "most null" vector

N=ND(1:kn+1).';
D=ND(kn+2:end).';

D1=D(1);
if D1==0; D1=1; end
% we have to make the D polynomial monic since thats
% what poly makes, so we have to first adjust N
N=N/D1;
D=D/D1;
% and then add the removed +/- Inf poles back in
D=poly([dinf; roots(D)]);

Dmax=max(abs(D));
if Dmax==0; Dmax=1; end
% normalize max Den value to be +/- 1
N=N/Dmax;
D=D/Dmax;

return

%a demo of this program is
clc
clear all
close all
x=[-5:1/16:5]';

y=gamma(x);
[N, D]=ratpolyfit(x,y,10,10);
yy=polyval(N,x)./polyval(D,x);
figure(1);plot(x,y,'b', x,yy,'dr');
grid on; axis([min(x) max(x) -15 15]);

y=erf(x);
[N, D]=ratpolyfit(x,y,10,10);
yy=polyval(N,x)./polyval(D,x);
figure(2);plot(x,y,'b', x,yy,'dr');
grid on; axis([min(x) max(x) -1.5 1.5]);

y=1./(1+x.*x);
[N, D]=ratpolyfit(x,y,10,10);
yy=polyval(N,x)./polyval(D,x);
figure(3);plot(x,y,'b', x,yy,'dr');
grid on; axis([min(x) max(x) 0 1]);

y=tan(x);
[N, D]=ratpolyfit(x,y,10,10);
yy=polyval(N,x)./polyval(D,x);
figure(4);plot(x,y,'b', x,yy,'dr');
grid on; axis([min(x) max(x) -5 5]);

y=exp(-x.*x);
[N, D]=ratpolyfit(x,y,10,10);
yy=polyval(N,x)./polyval(D,x);
figure(5);plot(x,y,'b', x,yy,'dr');
grid on; axis([min(x) max(x) 0 1]);

y=zeta(x);
[N, D]=ratpolyfit(x,y,10,10);
yy=polyval(N,x)./polyval(D,x);
figure(6);plot(x,y,'b', x,yy,'dr');
grid on; axis([min(x) max(x) -2 2]);

y=psi(x);
[N, D]=ratpolyfit(x,y,10,10);
yy=polyval(N,x)./polyval(D,x);
figure(7);plot(x,y,'b', x,yy,'dr');
grid on; axis([min(x) max(x) -5 5]);

x=x-min(x);

y=log(x);
[N, D]=ratpolyfit(x,y,10,10);
yy=polyval(N,x)./polyval(D,x);
figure(8);plot(x,y,'b', x,yy,'dr');
grid on; axis([min(x) max(x) -5 5]);

x=(x-min(x))*pi/2;

y=besselj(0,x);
[N, D]=ratpolyfit(x,y,10,10);
yy=polyval(N,x)./polyval(D,x);
figure(9);plot(x,y,'b', x,yy,'dr');
grid on; axis([min(x) max(x) -0.5 1]);

y=sin(x)./x;
[N, D]=ratpolyfit(x,y,10,10);
yy=polyval(N,x)./polyval(D,x);
figure(10);plot(x,y,'b', x,yy,'dr');
grid on; axis([min(x) max(x) -0.5 1]);

return

author: Sable


text:
function log(logger,level,message,varargin)
% Writes a message to all log files (and possibly standard out and
% standard error).
% Inputs:
%   level: The logging level.  See the constants defined in the
%     class header.  Several constants can be combined by bit oring
%     them.

assert(logger.validate())

for ii=1:length(logger.fidlist)
    if bitand(logger.fidlevels(ii),level)
	fprintf(logger.fidlist(ii),message,varargin{:});
    end
end
for ii=1:length(logger.filenamefids)
    if bitand(logger.filenamelevels(ii),level)
	fprintf(logger.filenamefids(ii),message,varargin{:});
    end
end
end

author: nesl


text:
function varargout = selectdata(varargin)

% This function is a backward compatibility wrapper. It allows existing
% MATLAB scripts that do not use the new FieldTrip ft_xxx function naming
% scheme to work with recent versions of the FieldTrip toolbox.
%
% Please look in ft_xxx for the help of the function that you are looking
% for, where xxx is the name of the function that you were looking for.

% Copyright (C) 2009, Robert Oostenveld

prefix    = 'ft_';
funname   = mfilename;
warning([upper(mfilename) ' is only a compatibility wrapper, which will soon be removed. Please instead call ' upper(prefix) upper(funname) '.']);
funhandle = str2func([prefix funname]);
[varargout{1:nargout}] = funhandle(varargin{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function X = mrdivide(A, B)
%/   Right matrix divide for an UNBNDFUN.
%   A/B divides the UNBNDFUN A by a scalar B. More generally, it gives the
%   least-squares solution (with respect to the continuous L^2 norm) to X*B = A
%   when either A or B is an UNBNDFUN. Note that in the latter case, formally it
%   is X.' that is returned, as UNBNDFUN objects are always columns.
%
% See also QR, RDIVIDE, MLDIVIDE.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Case B is a matrix: X*B = A  ==> X = A/B = (Q*R)/B = Q*(R/B)
%
% Case A is a matrix: X*B = X*(Q*R) = A ==> X = (A/R)*Q' ==> X' = Q*(A/R)'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ( (size(B, 2) ~= size(A, 2)) && ~(isa(B, 'double') && isscalar(B)) )
    error('CHEBFUN:UNBNDFUN:mrdivide:size', 'Matrix dimensions must agree.');

elseif ( isa(B, 'double') )  % UNBNDFUN / double

    if ( isscalar(B) )
        % Scalar case is easy:
        X = A;                              % Copy A to X
        X.onefun = X.onefun/B;              % mrdivide of the onefun
    else
        % Call MRDIVIDE at the ONEFUN level:
        X = A;
        X.onefun = (A.onefun/B);

%         % Alternatively, we could call QR() at the UNBNDFUN level
%         % For matrix case, we do least squares via QR:
%         [Q, R] = qr(A, 0);
%         X = Q*(R/B);
    end

elseif ( isa(A, 'double') )  % double / UNBNDFUN
    error('CHEBFUN:UNBNDFUN:mrdivide:doubleDivUnbndfun', ...
        ['/ does not support the division of a numerical matrix divided by ' ...
        'by an UNBNDFUN.']);

elseif ( isa(B, 'unbndfun') && isa(A, 'unbndfun') )
    error('CHEBFUN:UNBNDFUN:mrdivide:unbndfunDivUnbndfun', ...
        'Use ./ to divide UNBNDFUN by an UNBNDFUN.');

else
    error('CHEBFUN:UNBNDFUN:mrdivide:badArg', '%s/%s is not well-defined.', ...
        class(A), class(B));

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function y = abs( x )

%Disciplined convex/geometric programming information for ABS:
%   ABS(X) is convex and nonmonotonic in X. Therefore, when used in
%   DCPs, X must be affine. ABS(X) is not useful in DGPs, since all
%   log-convex and log-concave expressions are already positive.

%
% Determine the expression types
%

% 0 : convex, concave, invalid
% 1 : constant
% 2 : real affine
% 3 : complex affine
persistent remap
if isempty( remap ),
    remap_1 = cvx_remap( 'constant' );
    remap_2 = cvx_remap( 'real-affine' );
    remap_3 = cvx_remap( 'complex-affine' );
    remap_4 = cvx_remap( 'log-valid' );
    remap = remap_1 + ( 2 * remap_2 + 3 * remap_3 + 4 * remap_4 ) .* ~remap_1;
end
v = remap( cvx_classify( x ) );

%
% Process each type of expression one piece at a time
%

vu = sort( v(:) );
vu = vu([true;diff(vu)~=0]);
nv = length( vu );
sx = x.size_;
if nv ~= 1,
    y = cvx( sx, [] );
end
for k = 1 : nv,

    %
    % Select the category of expression to compute
    %

    vk = vu( k );
    if nv == 1,
        xt = x;
    else
        t = v == vk;
        xt = cvx_subsref( x, t );
    end

    %
    % Perform the computations
    %

    switch vk,
        case 0,
            % Invalid
            error( 'Disciplined convex programming error:\n    Illegal operation: abs( {%s} ).', cvx_class( xt ) );
        case 1,
            % Constant
            cvx_optval = cvx( builtin( 'abs', cvx_constant( xt ) ) );
        case 2,
            % Real affine
            w = [];
            st = size( xt );
            cvx_begin
                epigraph variable w( st )
                { xt, w } == lorentz( st, 0 ); %#ok
            cvx_end
        case 3,
            % Complex affine
            w = [];
            st = size( xt );
            cvx_begin
                epigraph variable w( st )
                { xt, w } == complex_lorentz( st, 0 ); %#ok
            cvx_end
        case 4,
            % log-affine, log-convex
            cvx_optval = xt;
        otherwise,
            error( 'Shouldn''t be here.' );
    end

    %
    % Store the results
    %

    if nv == 1,
        y = cvx_optval;
    else
        y = cvx_subsasgn( y, t, cvx_optval );
    end

end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX