text:
function all_sal = computeSaliency(options)

    if( ~exist( options.outfolder, 'dir' ) ), mkdir( options.outfolder ), end;
    if( ~exist( fullfile( options.outfolder, 'intra-frame_saliency'), 'dir' ) )
        mkdir(fullfile( options.outfolder, 'intra-frame_saliency'));
    end
    if( ~exist( fullfile( options.outfolder, 'inter-frame_saliency'), 'dir' ) )
        mkdir(fullfile( options.outfolder, 'inter-frame_saliency'));
    end
    if( ~exist( fullfile( options.outfolder, 'final_saliency'), 'dir' ) )
        mkdir(fullfile( options.outfolder, 'final_saliency'));
    end
    % Cache all frames in memory
    [data.frames,data.names,height,width,nframe ]= readAllFrames( options );
     % Load optical flow (or compute if file is not found)
    data.flow = loadFlow( options );
    if( isempty( data.flow ) )
        data.flow = computeOpticalFlow( options, data.frames );
    end
    % Load superpixels (or compute if not found)
    data.superpixels = loadSuperpixels( options );
    if( isempty( data.superpixels ) )
        data.superpixels = computeSuperpixels(  options, data.frames );
    end
    % Load Boundary (or compute if not found)
    % computeBoundary
    %

    [ superpixels, ~, bounds, labels ] = makeSuperpixelIndexUnique( data.superpixels );
    [ colours, centres, ~ ] = getSuperpixelStats( data.frames(1:nframe-1), superpixels, double(labels) );%
    valLAB = [];
    for index = 1:nframe-1
        valLAB = [valLAB;data.superpixels{index}.Sup1, data.superpixels{index}.Sup2, data.superpixels{index}.Sup3];
    end

    k =1:6:nframe-1;
    k(end) = [];
    global_saliency_track = zeros(height,width);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%computing global and local location cues %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i = 1:size(k,2)
        all_flowmagnitude = zeros(height,width);
        if i > 1
            for index = k(i):k(i)+5
                flowgradient = getFlowGradient( data.flow{index} );
                flowmagnitude = getMagnitude( flowgradient );
                flowmagnitude = (0.5+saliency_trac).*reshape(flowmagnitude(:),height,width);
                all_flowmagnitude = max(all_flowmagnitude,flowmagnitude);
            end
        else
            for index = k(i):k(i)+5
                flowgradient = getFlowGradient( data.flow{index} );
                flowmagnitude = getMagnitude( flowgradient );
                flowmagnitude = reshape(flowmagnitude(:),height,width);
                all_flowmagnitude = max(all_flowmagnitude,flowmagnitude);
            end
        end

        all_flowmagnitude = imresize(all_flowmagnitude,0.1,'bilinear');
        [h,w] = size(all_flowmagnitude);
        all_flowmagnitudex = all_flowmagnitude(:);
        all_flowmagnitudex(h*w+1)=0;
        all_Label = zeros(h,w);
        all_Label(:) = 1:size(all_flowmagnitude,1)*size(all_flowmagnitude,2);
        all_Labelx = zeros(h+2,w+2);
        all_Labelx(2:end-1,2:end-1) = all_Label;
        all_Labelx(1,:)=h*w+1;
        all_Labelx(end,:)=h*w+1;
        all_Labelx(:,1)=h*w+1;
        all_Labelx(:,end)=h*w+1;
        [ConSPix,~]= find_connect_superpixel(all_Labelx, h*w+1, h+2 ,w+2 );
        ConSPix = ConSPix + eye(size(ConSPix,1));
        all_fDistM = squareform(pdist(all_flowmagnitudex(:)));
        bdIds = GetBndPatchIds(all_Label,1);
        clipVal = EstimateDynamicParas(ConSPix,all_fDistM);
        geoDist = GeodesicSaliency(ConSPix, double(bdIds), all_fDistM, clipVal,true,[]);
        geoDist = geoDist(1:end-1);
        saliency_trac = reshape(geoDist,h,w);
        saliency_trac = imresize(saliency_trac,[height,width],'bilinear');
        global_saliency_track = max(global_saliency_track,saliency_trac);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%computing saliency via intra-frame graph %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if( options.vocal )
          disp( 'compute saliency via intra-frame graph:');
    end
    all_sal = cell(1,nframe-1);
    for index = 1:nframe-1
        if( options.vocal )
                fprintf( 'Processing frame: %i of %i... \n', ...
                    index, nframe-1);
        end
        FGR = [];
        frame = data.frames{index};
        frameName = data.names{index};
        nLabel = double(max(data.superpixels{index}.Label(:)));
        Label = data.superpixels{index}.Label;
        G = edge_detect(imfilter(uint8(frame),fspecial('average',3),'same','replicate'));
        flowgradient = getFlowGradient( data.flow{index} );
        flowmagnitude = getMagnitude( flowgradient );
        gradBoundary = 1 - exp( -flowmagnitude);
        flowmagnitude = G.*( gradBoundary +0.1);
        flowmagnitude = reshape(flowmagnitude(:),height,width);
        for  i = 1:nLabel
             flowmagnitude_R = flowmagnitude(Label==i);
             [flowmagnitude_R,~] = sort(flowmagnitude_R, 'descend');
             FGR(i)=mean(flowmagnitude_R(1:10));
             flowmagnitude(Label==i) = FGR(i);
        end

        [ConSPix,~]= find_connect_superpixel( Label, nLabel, height ,width );
        ConSPix = ConSPix + eye(nLabel);
        FGRM = squareform(pdist(FGR'));
        bdIds = GetBndPatchIds(Label,1);
        clipVal = EstimateDynamicParas(ConSPix,FGRM);
        geoDist = GeodesicSaliency(ConSPix, double(bdIds), FGRM, clipVal,true,[]);
        geo_sal = geoDist(data.superpixels{index}.Label);
        all_sal{index}=geo_sal.*(global_saliency_track+0.3);
        all_sal{index}=all_sal{index}/max(all_sal{index}(:));

        imwrite(all_sal{index}, [options.outfolder '\intra-frame_saliency\' frameName  '.bmp']);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%computing saliency via inter-frame graph %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if( options.vocal )
          disp( 'compute saliency via inter-frame graph:');
    end
    all_geoDist = cell(1,nframe-1);
    for index = 1:nframe-2
        if( options.vocal )
                fprintf( 'Processing frame: %i of %i... \n', ...
                    index, nframe-1);
        end
        nLabel = double(max(data.superpixels{index}.Label(:)));

        colDistM = squareform(pdist(valLAB( bounds(index):bounds(index+2)-1,:)));
        Conedge = [];
        Label_1 = data.superpixels{index}.Label;
        [~,conedge]= find_connect_superpixel( Label_1, max(Label_1(:)), height ,width );
        Conedge = [Conedge;conedge];

        [x,y] = meshgrid(1:bounds(index+1)-bounds(index),1:bounds(index+2)-bounds(index+1));
        conedge = [x(:),y(:)+bounds(index+1)-bounds(index)];
        connect = sum((centres(conedge(:,1)+bounds(index)-1,:) - centres(conedge(:,2)+bounds(index)-1,:)).^2,2 );
        cross_po_dis = conedge(connect<800,:);
        Conedge = [Conedge;cross_po_dis];

        Label_2 = data.superpixels{index+1}.Label;
        [~,conedge]= find_connect_superpixel( Label_2, max(Label_2(:)), height ,width );
        Conedge = [Conedge;conedge + bounds(index+1)-bounds(index)];
        ConSPix=sparse([Conedge(:,1);Conedge(:,2)],[Conedge(:,2);Conedge(:,1)], ...
         [ones(size(Conedge(:,1)));ones(size(Conedge(:,1)))],bounds(index+2)-bounds(index),bounds(index+2)-bounds(index));
        ConSPix = full(ConSPix);
        ConSPix = ConSPix +eye(size(ConSPix));

        firstmap = all_sal{index}>mean(all_sal{index}(:));
        str = strel('disk',1);  firstmap = imdilate(firstmap,str);
        fd = int32(firstmap).*data.superpixels{index}.Label;
        fd = unique(fd(:));
        fd(fd==0) = [];
        bd = setdiff(unique(data.superpixels{index}.Label(:)),fd);

        bdIds = bd;

        secondmap = all_sal{index+1}>mean(all_sal{index+1}(:));
        secondmap = ~((~secondmap).*(~firstmap));
        str = strel('disk',1);  secondmap = imdilate(secondmap,str);
        fd = int32(secondmap).*data.superpixels{index+1}.Label;
        fd = unique(fd(:));
        fd(fd==0) = [];
        bd = setdiff(unique(data.superpixels{index+1}.Label(:)),fd);
        bdIds = [bdIds;bd+bounds(index+1)-bounds(index)];

        clipVal = EstimateDynamicParas(ConSPix,colDistM);
        geoDist = GeodesicSaliency(ConSPix, double(bdIds), colDistM, clipVal, false,[]);
        geoDist_1 = geoDist(1:bounds(index+1)-bounds(index));
        geoDist_2 = geoDist(bounds(index+1)-bounds(index)+1:end);

        tmp = sort(geoDist_1, 'descend');
        pos = round(options.topRate * length(tmp));
        maxVal = tmp(pos);
        geoDist_1 = geoDist_1 / maxVal;
        geoDist_1(geoDist_1 > 1) = 1;

        tmp = sort(geoDist_2, 'descend');
        pos = round(options.topRate * length(tmp));
        maxVal = tmp(pos);
        geoDist_2 = geoDist_2 / maxVal;
        geoDist_2(geoDist_2 > 1) = 1;

        all_geoDist{index} = geoDist_1;
        all_geoDist{index+1} = geoDist_2;
        geo_sal = geoDist_1(data.superpixels{index}.Label);

        all_sal{index} = geo_sal*0.5+all_sal{index}*0.5;
        all_sal{index} = all_sal{index}./max( all_sal{index}(:));

        L{1} = uint32(data.superpixels{index}.Label);
        S{1} = repmat(all_sal{index},[1 3]);
        [ R, ~, ~ ] = getSuperpixelStats(S(1:1),L, nLabel );
        R = double(R(:,1));
        sR = sort(R);
        t = sum(sR(end-9:end))/10;
        R = (R-min(R))/(t-min(R));
        R(R>1)=1;
        all_geoDist{index} = R';

        geo_sal = geoDist_2(data.superpixels{index+1}.Label);

        all_sal{index+1} = geo_sal*0.5+all_sal{index+1}*0.5;
        all_sal{index+1} = all_sal{index+1}./max( all_sal{index+1}(:));

        nLabel = double(max(data.superpixels{index+1}.Label(:)));
        L{1} = uint32(data.superpixels{index+1}.Label);
        S{1} = repmat(all_sal{index+1},[1 3]);
        [ R, ~, ~ ] = getSuperpixelStats(S(1:1),L, nLabel );
        R = double(R(:,1));
        sR = sort(R);
        t = sum(sR(end-9:end))/10;
        R = (R-min(R))/(t-min(R));
        R(R>1)=1;
        all_geoDist{index+1} = R';

    end

    RegionSal = [];
    for index = 1:nframe-1
        nLabel = double(max(data.superpixels{index}.Label(:)));
        all_sal{index} = all_sal{index}.*(global_saliency_track+0.1);
        all_sal{index} = all_sal{index}/max(all_sal{index}(:));

        sal = reshape( all_sal{index},height*width,1);
        L{1} = uint32(data.superpixels{index}.Label);
        S{1} = repmat(sal,[1 3]);
        [ R, ~, ~ ] = getSuperpixelStats(S(1:1),L, nLabel );
        R = double(R(:,1));
        sR = sort(R);
        t = sum(sR(end-9:end))/10;
        R = (R-min(R))/(t-min(R));
        R(R>1)=1;
        all_geoDist{index} = R';
        RegionSal = [RegionSal;all_geoDist{index}(:)];
        all_sal{index} = double(R(data.superpixels{index}.Label));
        imwrite(all_sal{index}, [options.outfolder '\inter-frame_saliency\' data.names{index} '.bmp']);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Spatio-temporal consistance optimization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%you can uncomment this part for computation efficiency %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if( options.vocal )
          disp( 'spatio-temporal consistance optimization');
    end
    Conedge = [];
    for index = 1:nframe-1
        Label = data.superpixels{index}.Label;
        [~,conedge]= find_connect_superpixel( Label, max(Label(:)), height ,width );
        Conedge = [Conedge;conedge + bounds(index)-1];
    end
    intralength = size(Conedge,1);
    for index = 1:nframe-2
        [x,y] = meshgrid(1:bounds(index+1)-bounds(index),1:bounds(index+2)-bounds(index+1));
        conedge = [x(:)+bounds(index)-1,y(:)+bounds(index+1)-1];
        connect = sum((centres(conedge(:,1),:) - centres(conedge(:,2),:)).^2,2 );
        Conedge = [Conedge;conedge(connect<800,:)];
    end
    valDistances=sqrt(sum((valLAB(Conedge(:,1),:)-valLAB(Conedge(:,2),:)).^2,2));
    valDistances(intralength+1:end)=valDistances(intralength+1:end)/5;
    valDistances=normalize(valDistances);
    weights=exp(-options.valScale*valDistances)+ 1e-5;
    weights=sparse([Conedge(:,1);Conedge(:,2)],[Conedge(:,2);Conedge(:,1)], ...
    [weights;weights],labels,labels);
    E = sparse(1:labels,1:labels,ones(labels,1)); iD = sparse(1:labels,1:labels,1./sum(weights));
    P = iD*weights;
    RegionSal = (E-P+10*options.alpha*E)\RegionSal;
    for index = 1:nframe-1
        R = RegionSal(bounds(index):bounds(index+1)-1);
        sR = sort(R);
        t = sum(sR(end-9:end))/10;
        R = (R-min(R))/(t-min(R));
        R(R>1)=1;
        all_geoDist{index} = R';
        all_sal{index} = 0.6*double(R(data.superpixels{index}.Label))+0.4*all_sal{index};
        imwrite(all_sal{index}, [options.outfolder '\final_saliency\' data.names{index}  '.bmp']);
    end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ rotImg, R ] = rotatePanorama( img, vp, R )
%ROTATEPANORAMA Rotate panorama
%   if R is given, vp (vanishing point) will be overlooked
%   otherwise R is computed from vp

[sphereH, sphereW, C] = size(img);
% rotImg = zeros( sphereH, sphereW, C);

%% new uv coordinates
[TX TY] = meshgrid(1:sphereW, 1:sphereH);
TX = TX(:);
TY = TY(:);
ANGx = (TX- sphereW/2 -0.5)/sphereW * pi *2 ;
ANGy = -(TY- sphereH/2 -0.5)/sphereH * pi;
uvNew = [ANGx ANGy];
xyzNew = uv2xyzN(uvNew,1);

%% rotation matrix
if nargin<3
    R = diag([1 1 1])/(vp');
end

xyzOld = (R\xyzNew')';
uvOld = xyz2uvN(xyzOld, 1);

% Px = uvOld(:,1)/2/pi*sphereW + 0.5 + sphereW/2;
% Py = -uvOld(:,2)/pi*sphereH + 0.5 + sphereH/2;
Px = (uvOld(:,1)+pi) / (2*pi) * sphereW + 0.5;
Py = (-uvOld(:,2) + pi/2) / pi * sphereH + 0.5;

Px = reshape(Px, [sphereH sphereW]);
Py = reshape(Py, [sphereH sphereW]);

% boundary
imgNew = double(zeros(sphereH+2, sphereW+2, C));
imgNew(2:end-1, 2:end-1, :) = img;
imgNew(2:end-1,1,:) = img(:,end,:);
imgNew(2:end-1,end,:) = img(:,1,:);
imgNew(1,2:sphereW/2+1,:) = img(1,sphereW:-1:sphereW/2+1,:);
imgNew(1,sphereW/2+2:end-1,:) = img(1,sphereW/2:-1:1,:);
imgNew(end,2:sphereW/2+1,:) = img(end,sphereW:-1:sphereW/2+1,:);
imgNew(end,sphereW/2+2:end-1,:) = img(1,sphereW/2:-1:1,:);
imgNew(1,1,:) = img(1,1,:);
imgNew(end,end,:) = img(end,end,:);
imgNew(1,end,:) = img(1,end,:);
imgNew(end,1,:) = img(end,1,:);

rotImg = warpImageFast(imgNew, Px+1, Py+1);
% rotImg = warpImageFast(img, Px, Py);

end

author: zouchuhang

repo: LayoutNet


text:
function [vec, outnames] = contrast_code(vec)
% :Usage:
% ::
%
%     [vec,outnames] = contrast_code(vec)
%
% Changes values to 1, -1, or 0 for contrast coding

outnames = [];

if iscell(vec) && ~isempty(vec)
    [indic, names]  = string2indicator(vec);

    for i = 2:length(names)
        vec = indic(:, i - 1) - indic(:, i);
        outnames{i - 1} = [names{i - 1} ' - ' names{i}];
    end

elseif ~isempty(vec)
    wh = find(vec > 0);

    vec(wh) = 1;

    wh = find(vec < 0);

    vec(wh) = -1;

end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function d = sqdist(x, y)
% L2 dist from x to y
% x: [1 x d]
% y: [n x d]
% d: [n x 1]

d = sqrt(sum((repmat(x, size(y,1), 1) - y).^2, 2));


text:
function [P, G] = getPoints(q,p)
% [P, G] = getPoints(q,p)
%
% This function computes the joint positions (P) and center of mass
% positions (G) for the five-link biped, in configuration given by q with
% physical parameters in p.
%
% INPUTS:
%   q = [5, n] = link configuration
%   p = parameters struct
%
% OUTPUTS:
%   P = [10, n] = joint positions [x;y;x;y;...]
%   G = [10, n] = CoM positions {x;y;x;y;...]
%

q1 = q(1,:);  %stance leg tibia angle
q2 = q(2,:);  %stance leg femur angle
q3 = q(3,:);  %torso angle
q4 = q(4,:);  %swing leg femur angle
q5 = q(5,:);  %swing leg tibia angle

[P,G] = autoGen_getPoints(...
    q1,q2,q3,q4,q5,...
    p.l1 ,p.l2 ,p.l3 ,p.l4 ,p.l5 ,p.c1 ,p.c2 ,p.c3 ,p.c4 ,p.c5);

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function D = EuDist2(fea_a,fea_b,bSqrt)
%EUDIST2 Efficiently Compute the Euclidean Distance Matrix by Exploring the
%Matlab matrix operations.
%
%   D = EuDist(fea_a,fea_b)
%   fea_a:    nSample_a * nFeature
%   fea_b:    nSample_b * nFeature
%   D:      nSample_a * nSample_a
%       or  nSample_a * nSample_b
%
%    Examples:
%
%       a = rand(500,10);
%       b = rand(1000,10);
%
%       A = EuDist2(a); % A: 500*500
%       D = EuDist2(a,b); % D: 500*1000
%
%   version 2.1 --November/2011
%   version 2.0 --May/2009
%   version 1.0 --November/2005
%
%   Written by Deng Cai (dengcai AT gmail.com)

if ~exist('bSqrt','var')
    bSqrt = 1;
end

if (~exist('fea_b','var')) || isempty(fea_b)
    aa = sum(fea_a.*fea_a,2);
    ab = fea_a*fea_a';

    if issparse(aa)
        aa = full(aa);
    end

    D = bsxfun(@plus,aa,aa') - 2*ab;
    D(D<0) = 0;
    if bSqrt
        D = sqrt(D);
    end
    D = max(D,D');
else
    aa = sum(fea_a.*fea_a,2);
    bb = sum(fea_b.*fea_b,2);
    ab = fea_a*fea_b';

    if issparse(aa)
        aa = full(aa);
        bb = full(bb);
    end

    D = bsxfun(@plus,aa,bb') - 2*ab;
    D(D<0) = 0;
    if bSqrt
        D = sqrt(D);
    end
end

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function [cmps2] = ftps22cmps2(ftps2)
% Convert acceleration from feet per square-second to nanometers per second
% squared.
% Chad A. Greene 2012
cmps2 = ftps2*30.4800000;
end

author: Sable

repo: mcbench-benchmarks


text:
function [kl1,kl2] = mvnkl(Mu1,Sigma1,Mu2,Sigma2)
%MVNKL Kullback-Leibler divergence between two multivariate normal pdfs.

D = numel(Mu1);

Mu1 = Mu1(:);
Mu2 = Mu2(:);

dmu = Mu2 - Mu1;
detq1 = det(Sigma1);
detq2 = det(Sigma2);
lndet = log(detq2 / detq1);

kl1 = 0.5*(trace(Sigma2\Sigma1) + dmu'*(Sigma2\dmu) - D + lndet);
if nargout > 1
    kl2 = 0.5*(trace(Sigma1\Sigma2) + dmu'*(Sigma1\dmu) - D - lndet);
end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function [ area, radin, side ] = polygon_outrad_data_2d ( n, radout )

%*****************************************************************************80
%
%% POLYGON_OUTRAD_DATA_2D determines polygonal data from its outer radius in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of sides of the polygon.
%    N must be at least 3.
%
%    Input, real RADOUT, the outer radius of the polygon, that is,
%    the radius of the smallest circle that can be described
%    around the polygon.
%
%    Output, real AREA, the area of the regular polygon.
%
%    Output, real RADIN, the inner radius of the polygon, that is,
%    the radius of the largest circle that can be inscribed
%    within the polygon.
%
%    Output, real SIDE, the length of one side of the polygon.
%
  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'POLYGON_OUTRAD_DATA_2D - Fatal error!\n' );
    fprintf ( 1, '  Input value of N must be at least 3.\n' );
    fprintf ( 1, '  but your input value was N = %d\n', n );
    error ( 'POLYGON_OUTRAD_DATA_2D - Fatal error!' );
  end

  angle = pi / n;
  area = 0.5 * n * radout * radout * sin ( 2.0 * angle );
  side = 2.0 * radout * sin ( angle );
  radin = 0.5 * side / tan ( angle );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [childId, packetId] = calcVariantId(parentOri,childOri,p2c,varargin)
% compute variantIds and packetId from parent / child orientation pairs
%
% Syntax
%
%   variantId = calcVariantId(parentOri,childOri,p2c)
%
%   % compute packetIds
%   [variantId,packetId] = calcVariantId(parentOri,childOri,p2c,...
%     {hklParent,hklChild})
%
%   % packet determination
%   hklParent = Miller({1,1,1},{1,-1,1},{-1,1,1},{1,1,-1},p2c.CS);
%   hklChild  = Miller(1,0,1,p2c.SS);
%
%   [variantId,packetId] = calcVariantId(parentOri,childOri,p2c,...
%     {hklParent,hklChild})
%
% Input
%  parentOri - parent @orientation
%  childOri  - child @orientation
%  p2c       - parent to child mis@orientation
%  hklParent, hklChild - correspondent planes between parent and child
%
% Output
%  variantId - variant id
%  packetId  - packet id
%

% all child variants
childVariants  = variants(p2c, parentOri);

if size(childVariants,1) == 1
  childVariants = repmat(childVariants,length(childOri),1);
end

% compute distance to all possible variants
d = dot(childVariants,repmat(childOri(:),1,size(childVariants,2)));

% take the best fit
[~,childId] = max(d,[],2);

% compute packetId if required
if nargout == 2
  % Get packet definition
  tmp = getClass(varargin,'cell');
  isMiller = [];
  for ii = 1:length(tmp); isMiller(ii) = ~isempty(getClass(tmp(ii),'Miller')); end

  if sum(isMiller) == 2 % definition given
    ind = find(isMiller);
    h1 = tmp{ind(1)};
    h2 = tmp{ind(2)};
  else % definition assumed
    warning('Packet ID calculation assuming {111}_p||{110}_c');
    h1 = Miller({1,1,1},{1,-1,1},{-1,1,1},{1,1,-1},p2c.CS);
    h2 = Miller(1,0,1,p2c.SS);
  end

  omega = dot(variants(p2c,h1),h2);

  [~,packetId] = max(omega,[],1);

  packetId = packetId(childId);

end

end

author: mtex-toolbox

repo: mtex

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


text:
function res = lt(a,b)
%LT           Implements  a < b  elementwise for intervals a and b
%
%  if true,  a  is definitely less than  b
%

% written  10/16/98     S.M. Rump
% modified 11/30/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if ~isa(a,'intval')
    a = intval(a);
  end
  if ~isa(b,'intval')
    b = intval(b);
  end

  if a.complex | b.complex
    res = real(sup(a)) < real(inf(b)) & imag(sup(a)) < imag(inf(b)) ;
  else
    res = sup(a) < inf(b) ;
  end

repo: OpenMAS


text:
function c8_mag_test ( )

%*****************************************************************************80
%
%% C8_MAG_TEST tests C8_MAG.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2015
%
%  Author:
%
%    John Burkardt
%
  seed = 123456678;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8_MAG_TEST\n' );
  fprintf ( 1, '  C8_MAG computes the magnitude of a C8.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '       C1=C8_UNIFORM_01          R2=C8_MAG(C1)             R3=ABS(C1))\n' );
  fprintf ( 1, '     ---------------------     ---------------------     ---------------------\n' );
  fprintf ( 1, '\n' );

  for test = 1 : 10

    [ c1, seed ] = c8_uniform_01 ( seed );
    r2 = c8_mag ( c1 );
    r3 = abs ( c1 );

    fprintf ( 1, '  (%12f  %12f)  %12f  %12f\n', ...
      real ( c1 ), imag ( c1 ), r2, r3 );

  end

  return
end


text:
function [parent] = multiNonUnifMutation(parent,bounds,Ops)
% Multi-Non uniform mutation changes all of the parameters of the parent
% based on a non-uniform probability distribution.  This Gaussian
% distribution starts wide, and narrows to a point distribution as the
% current generation approaches the maximum generation.
%
% function [newSol] = multiNonUnifMutate(parent,bounds,Ops)
% parent  - the first parent ( [solution string function value] )
% bounds  - the bounds matrix for the solution space
% Ops     - Options for multiNonUnifMutation
%          [gen #MultiNonUnifMutations maxGen b]

% Binary and Real-Valued Simulation Evolution for Matlab
% Copyright (C) 1996 C.R. Houck, J.A. Joines, M.G. Kay
%
% C.R. Houck, J.Joines, and M.Kay. A genetic algorithm for function
% optimization: A Matlab implementation. ACM Transactions on Mathmatical
% Software, Submitted 1996.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 1, or (at your option)
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. A copy of the GNU
% General Public License can be obtained from the
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

cg=Ops(1); 				% Current Generation
mg=Ops(3); 				% Maximum Number of Generations
b=Ops(4);                               % Shape parameter
df = bounds(:,2) - bounds(:,1); 	% Range of the variables
numVar = size(parent,2)-1; 		% Get the number of variables
% Now mutate that point
md = round(rand(1,numVar));
for i = 1:numVar
  if md(i)
    parent(i)=parent(i)+delta(cg,mg,bounds(i,2)-parent(i),b);
  else
    parent(i)=parent(i)-delta(cg,mg,parent(i)-bounds(i,1),b);
  end
end

author: HuangCongQing


text:
function x = product_grid ( m, a, b, n, c, d )

%*****************************************************************************80
%
%% PRODUCT_GRID produces a product grid of points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real M, the number of grid points in X.
%
%    Input, real A, B, the X limits.
%
%    Input, real N, the number of grid points in Y.
%
%    Input, real C, D, the Y limits.
%
%    Output, real X(M*N,2), the points of the product grid.
%
  x = zeros ( m * n, 2 );

  k = 0;
  for i = 1 : m
    for j = 1 : n
      k = k + 1;
      x(k,1) = ( ( m - i ) * a + ( i - 1 ) * b ) / ( m - 1 );
      x(k,2) = ( ( n - j ) * c + ( j - 1 ) * d ) / ( n - 1 );
    end
  end

  return
end


text:
function [text] = writeTXT(para,fname_out)
% Write a txt file for online PD map to highlight specific reaction nodes.
%
% USAGE:
%
%    [text] = writeTXT(para, fname_out)
%
% INPUTS:
%    para:           A variable that stores two columns: the first column
%                    contains a list of reaction names,whereas the second
%                    column contains a corresponding list of the colours (Hex
%                    triplet, e.g.,
%                    https://closedxml.codeplex.com/wikipage?title=Excel%20Indexed%20Colors)
%    fanme_out:      The name of the output text file name.
%
% OPTIONAL OUTPUT:
%    text:           The lines of the text file
%
% .. Author: - Longfei Mao Oct/2014

if nargin<2 || isempty(fname_out)

    [fname_out, fpath]=uiputfile('*.txt','for Online PD map');
    if(fname_out==0)
        return;
    end
    f_out=fopen([fpath,fname_out],'w');
else
    f_out=fopen(fname_out,'w');
end

%%%%%%%%%%%% read

% blank={' '}
[r,c]=size(para)

if c==1;
    para(:,2)={'#FF0000FF'}; % By default, the colour hex code is set to be red.
end

if ~isempty(strfind(para{1,2},'#'))||~isempty(strfind(para{1,2},'FF'))
    r=r+1;
    s=2;
    text{1,1}=['name', ' ', 'color']; % add title row;
else
    s=1;r=r+1;
end

for rr=s:r; % content rows + one title row
    name=para{rr-1,1};

    if ~isempty(para{rr-1,2});
        color=para{rr-1,2};

        if length(color)>6;
            st=strfind(color,'FF');  % remove 'FF' necessary for CellDesigner files
            if ~isempty(st)
                color=color(st+2:end);
            end
            if isempty(strfind(color,'#'));  % add # in front of the color code.
                color=['#',color];
                %errordlg('The format of the colour code is not correct');
            end

        end

    end

    %     if rr==s;
    %         text{1,1}=['name',' ', 'color'];
    %
    %     else
    text{rr,1}=[name,' ',color];

end

%%%%%%% Write to the file %%%%

for d=1:length(text);

    fprintf(f_out,'%s\n',text{d});

end

fclose(f_out);

end

author: opencobra

repo: cobratoolbox


text:
% a tutorial adapted from the fipy diffusion 1D example
% see: http://www.ctcms.nist.gov/fipy/examples/diffusion/index.html
% How to convert a PDE into an ODE and solve it using a matlab ode solver
% needs more documentation
% If you want to have better time-stepping with Matlab ODE solver, use this example
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc
clear

%% define the domain
L = 50;  % domain length
Nx = 20; % number of cells
meshstruct = createMesh1D(Nx, L);
BC = createBC(meshstruct); % all Neumann boundary condition structure
BC.left.a = 0; BC.left.b=1; BC.left.c=1; % left boundary
BC.right.a = 0; BC.right.b=1; BC.right.c=0; % right boundary
x = meshstruct.cellcenters.x;
%% define the transfer coeffs
D_val = 1;
D = createFaceVariable(meshstruct, D_val);
alfa = createCellVariable(meshstruct, 1);
%% define initial values
c_old = createCellVariable(meshstruct, 0); % initial values
c_value = c_old;
%% loop
dt = 0.1; % time step
final_t = 100;
Mdiff = diffusionTerm(D);
[M, RHS] = combineBC(BC, Mdiff, zeros(Nx+2,1));
% eq: dcdt = D d2c/dx2
dcdt = @(t,c)(M*c-RHS);
[t_temp, c_temp] = ode45(dcdt, [0 final_t], internalCells(c_old));
c_analytical = 1-erf(x./(2*sqrt(D_val*final_t)));
plot(x, c_temp(end,:), x, c_analytical, 'o')

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function r = get_consensus_set_rank(CS, E, mode, T_noise_squared)

% r = get_consensus_set_rank(CS, E, mode, T_noise_squared)
%
% DESC:
% get the rank of a consensus set
%
% AUTHOR
% Marco Zuliani - marco.zuliani@gmail.com
%
% VERSION:
% 1.0.2
%
% INPUT:
% CS                = logical array identifyng the CS
% E                 = error associated to each data element
% mode              = specify the type of ranking
%                     0 -> RANSAC (cardinality of the CS)
%                     1 -> MSAC
% T_noise_squared   = noise threshold to discriminate inliers vs outliers
%
% OUTPUT:
% r                 = consensus set rank

% HISTORY:
%
% 1.0.0             - 01/13/08 - Initial version
% 1.0.1             - 01/17/08 - Added different M-estimators
% 1.0.2             - 02/22/08 - Removed the M-estimators
%                                Added different ranking modes
% 1.0.3             - 05/26/14 - Removed MLESAC

N = length(CS);

switch mode

    case 'RANSAC'

        r = sum(CS);

    case 'MSAC'

        rho = E;
        rho(rho >= T_noise_squared) = T_noise_squared;

        % set to negative sothat also this rank is to be
        % maximized
        r = -sum(rho)/N;

    otherwise

        error('RANSACToolbox:optionError', 'Unknown ranking mode');

end

return

author: RANSAC

repo: RANSAC-Toolbox


text:
% Calculate the value of lambda so that if lambda >= lambdamax, the TVD
% functional solved by l1pwc is minimized by the trivial constant
% solution x = mean(y). This can then be used to determine a useful range
% of values of lambda, for example.
%
% Usage:
% lambdamax = l1pwclmax(y)
%
% Input arguments:
% - y          Original signal to denoise, size N x 1.
%
% Output arguments:
% - lambdamax  Value of at which x = mean(y) is the output of the l1pwc
%              function.
%
% (c) Max Little, 2010. If you use this code for your research, please
% cite:
% M.A. Little, Nick S. Jones (2010)
% "Sparse Bayesian Step-Filtering for High-Throughput Analysis of Molecular
% Machine Dynamics", in 2010 IEEE International Conference on Acoustics,
% Speech and Signal Processing, 2010, ICASSP 2010 Proceedings.
%

function lambdamax = ML_l1pwclmax(y)

narginchk(1,1);
y = y(:);
N = length(y);
M = N - 1;

% Construct sparse operator matrices
I1 = speye(M,M);
O1 = spalloc(M,1,M);
D = [I1 O1]-[O1 I1];

DDT = D*D';
Dy  = D*y;

lambdamax = max(abs(DDT\Dy));

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [fusion,w0] = qfuser_v4(w,scores,wfuse)
% qindx: index set for rows of scores.scores which are per-trial quality
%        measures.
%
% sindx: index set for rows of scores.scores which are normal discriminative
%        scores.

if nargin==0
    test_this();
    return;
end

sindx = scores.sindx;
qindx = scores.qindx;
m =length(sindx);

% Create building blocks
[Cal,w1] = parallel_cal([],scores.scores(sindx,:),wfuse);
[whead,wtail] = splitvec_fh(length(w1));
Cal = Cal(whead);

[LLH1,w2] = QQtoLLH([],scores.modelQ,scores.segQ,m);
[whead,wtail] = splitvec_fh(length(w2),wtail);
LLH1 = LLH1(whead);
W2 = reshape(w2,[],m);
W2(:) = 0;
W2(end,:) = 0.5/(m+1);
w2 = W2(:);

[LLH2,w3] = QtoLLH([],scores.scores(qindx,:),m);
LLH2 = LLH2(wtail);

LLH =  sum_of_functions([],[1,1],LLH1,LLH2);
%LLH = LLH1;

P = LLH;
%P = exp_mv2df(logsoftmax_trunc_mv2df(LLH,m));

w0 = [w1;w2;w3];

% Assemble building blocks

%  modulate linear fusion with quality
fusion = sumcolumns_fh(m,dottimes_of_functions(w,P,Cal));

end

function test_this()

m = 5;
k = 2;
n1 = 4;
n2 = 5;

scores.sindx = [1,2,3];
scores.qindx = [4,5];

scores.scores = randn(m,n1*n2);
scores.modelQ = randn(k,n1);
scores.segQ = randn(k,n2);

wfuse = [1,2,3,4]';

[fusion,w0] = qfuser_v4([],scores,wfuse);

%test_MV2DF(fusion,w0);

[fusion(w0),linear_fuser(wfuse,scores.scores(scores.sindx,:))]

%fusion(w0)

end

author: nesl


text:
function ft_postamble(cmd, varargin)

% FT_POSTAMBLE is a helper function that is included in many of the FieldTrip
% functions and which takes care of some general settings and operations at the end
% of the function.
%
% This ft_postamble m-file is a function, but internally it executes a number of
% private scripts in the callers workspace. This allows the private script to access
% the variables in the callers workspace and behave as if the script were included as
% a header file in C-code.
%
% See also FT_PREAMBLE

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

% ideally this would be a script, because the local variables would then be
% shared with the calling function. Instead, this is a function which then
% passes the variables explicitely to another script which is eval'ed.

% the following section ensures that these scripts are included as
% dependencies when using the MATLAB compiler
%
%#function ft_postamble_debug
%#function ft_postamble_provenance
%#function ft_postamble_previous
%#function ft_postamble_history
%#function ft_postamble_savevar
%#function ft_postamble_savefig
%#function ft_postamble_randomseed
%#function ft_postamble_hastoolbox

% this is a trick to pass the input arguments into the ft_postamble_xxx script
assignin('caller', 'postamble_argin', varargin);

full_cmd = ['ft_postamble_' cmd];
cmd_exists = false;

if exist(full_cmd, 'file')
  % Matlab can find commands in a private subdirectory; Octave cannot.
  % If pwd is already the private directory, or if using Matlab then
  % the command can be evaluated directly
  cmd_exists = true;

elseif ~ft_platform_supports('exists-in-private-directory')
  % Octave does not find files by name in a private directory, so the full filename must be specified.
  private_dir = fullfile(fileparts(which(mfilename)),'private');
  full_path   = fullfile(private_dir,[full_cmd '.m']);
  cmd_exists  = exist(full_path, 'file');
  full_cmd_parts = {'ft_tmp_orig_pwd = pwd(); ',...
                    'ft_tmp_orig_pwd_cleaner = onCleanup(@()cd(ft_tmp_orig_pwd)); ',...
                    sprintf('cd(''%s''); ',private_dir),...
                    [full_cmd '; '],...
                    'clear ft_tmp_orig_pwd_cleaner;'};
  full_cmd = sprintf('%s', full_cmd_parts{:});
end

if ~cmd_exists
  % XXX earlier versions would not do anything if ~cmd_exists,
  % but fail silently (without raising an error or warning).
  % Should that behavior be kept?
  ft_error('Could not run %s - does not seem to exist', full_cmd);
end

evalin('caller', full_cmd);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%
% This script will take a given set of fiber groups loaded in dtiFiberUI,
% which cross the callosum, and find 11 points along each fiber that are
% cented around the midsagital point of each path (5 steps to each side
% of the center point).
%
% For each of the fiber groups loaded it will then return the means for FA, MD,
% RD, AD, and selected background image value, in a Nx5 array stored in the
% variable mnVal. Currently this script is set up for use with fmap background
% images.
%
% The standard deviation will also be computed for each of the 5 vals and
% stored in a Nx5 array in the variable sdVal.
%
% The names and sequence of fiber grous can be found in fg.name
%
%
% HISTORY:
% 2008.08.08 RFD Wrote it.

%dataDir = '/biac3/wandell4/data/reading_longitude/dti_adults/ah080521_sense/dtiss06';
dataDir = '/biac3/wandell4/data/reading_longitude/dti_adults/rfd050504_SENSE/dti23';
roiDir = fullfile(fileparts(dataDir),'ROIs');

handles = guidata(gcf);
nSteps = 5;

fg = handles.fiberGroups;
nFg = numel(fg);

% Create a legend for the fiber groups
for(ii=1:nFg)
    fgCol(ii,:) = fg(ii).colorRgb;
    fgStr{ii} = fg(ii).name;
end
legImg = ones(18*nFg-1,16,3);
for(ii=1:nFg)
    for(jj=1:3)
    	yPos = (ii-1)*18+1;
    	legImg(yPos:yPos+16,:,jj) = fgCol(ii,jj)./255;
    end
end
figure(88); image(legImg); axis equal off tight;
set(gca,'units','pixels','position',[8 10 size(legImg,2) size(legImg,1)]);
for(ii=1:nFg)
  text(18,(ii-1)*18+9, strrep(fgStr{ii},'_',' '),'FontSize',10);
end

% Get the callosal ROI
ccRoi = dtiReadRoi(fullfile(roiDir,'CC'));
minDist = 0.87;

h = mrvWaitbar(0,'Processing fibers...');
msCoords = [];
for(ii=1:nFg)
    fiberCoords{ii} = [];
    for(jj=1:numel(fg(ii).fibers))
        fc = fg(ii).fibers{jj};
        % first find those points that are within the CC ROI
        [indices, bestSqDist] = nearpoints(fc, ccRoi.coords');
        keepAll = bestSqDist<=minDist^2;
        if(any(keepAll))
            midSagPos = min(abs(fc(1,keepAll)));
            midSagInd = find(abs(fc(1,:))==midSagPos & keepAll);
            if(numel(midSagInd)~=1)
                disp('ignoring fiber');
            else
                msCoords = horzcat(msCoords, fc(:,midSagInd));
                midSagInds = [midSagInd-nSteps:midSagInd+nSteps];
                if(midSagInds(1)>1 && midSagInds(end)<size(fc,2))
                    fiberCoords{ii} = horzcat(fiberCoords{ii}, fc(:,midSagInds));
                end
            end
        end
    end
    mrvWaitbar(ii/nFg,h);
end
close(h);
figure; subplot(2,1,1);
c = horzcat(fiberCoords{:});
plot(c(2,:),c(3,:),'b.');
axis equal;
subplot(2,1,2);
plot(msCoords(2,:),msCoords(3,:),'r.'); hold on;
plot(ccRoi.coords(:,2),ccRoi.coords(:,3),'ko');
hold off; axis equal tight;

% Extract the values and summarize the callosal segments
bg = handles.bg(dtiGet(handles,'curbgnum'));
bg.img = bg.img.*(bg.maxVal-bg.minVal)+bg.minVal;
for(ii=1:nFg)
    [ev1,ev2,ev3] = dtiGetValFromTensors(handles.dt6, fiberCoords{ii}, ...
        inv(handles.xformToAcpc), 'eigvals', 'nearest');
    [fa,md,rd,ad] = dtiComputeFA(horzcat(ev1,ev2,ev3));
    gv = fa<1;
    bgCoords = round(mrAnatXformCoords(inv(bg.mat),fiberCoords{ii}));
    bgInds = sub2ind(size(bg.img),bgCoords(:,1),bgCoords(:,2),bgCoords(:,3));
    bgVals = bg.img(bgInds);
    % get rid of bad values (values that are exactly zero)
    bgVals = bgVals(bgVals>0 & gv);
    mnVal(ii,:) = [mean(fa(gv)),mean(md(gv)),mean(rd(gv)),mean(ad(gv)),mean(bgVals)];
    sdVal(ii,:) = [std(fa(gv)),std(md(gv)),std(rd(gv)),std(ad(gv)),std(bgVals)];
end

error('stop here');

handles = guidata(gcf);
ccfg = handles.fiberGroups(end);
% concatenate all fibers into one big list of points
ccCoords = horzcat(ccfg.fibers{:});
[ev1,ev2,ev3] = dtiGetValFromTensors(handles.dt6, ccCoords, inv(handles.xformToAcpc), 'eigvals', 'nearest');
[fa,md,rd,ad] = dtiComputeFA(horzcat(ev1,ev2,ev3));
bgCoords = round(mrAnatXformCoords(inv(bg.mat),ccCoords));
bgInds = sub2ind(size(bg.img),bgCoords(:,1),bgCoords(:,2),bgCoords(:,3));
bgVals = bg.img(bgInds);
% get rid of bad values (values that are exactly zero)
gv = fa<1 & md>0.5 & md<1.2 & bgVals>0 & bgVals<0.27;
figure;plot(md(gv),bgVals(gv),'.')

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Doyen Sahoo
% Contributors: Bin LI, Steven C.H. Hoi
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ ] = experimenterMenu( job )

    % The job must contain:
    % 1. Array of algorithm IDs
    % 2. Default parameters for each algorithm
    % 3. The dataset choice

    % Read configuration
    load ../GUI/config/config.mat;

    menuId = 3; % Experimenter Menu
    displayMenu(menuId); %

    prompt = 'Please enter your choice (1-5):';
    choice = input(prompt);

    switch(choice),
        case 1,
            % Select Algorithms
            job     = selectAlgorithms (job);
            experimenterMenu( job );

        case 2,
            % Set Parameters
            disp('The current job is:');
            displayExperimenterJob(job);
            prompt = 'Which Algorithm Parameters do you want to change: ';
            choice = input(prompt);
            [ parameters ] = selectParameters(job.algorithmId(choice));

            p = length(parameters);
            for i = 1:1:p
                job.parameters(choice,i) = num2cell(parameters(i));
            end
            experimenterMenu(job);

        case 3
            % Select dataset
            [ dataId ] = dataSelectMenu(dataList );
            if dataId ~= -1
                job.datasetId   = dataId;
                DL              = char(dataList);
                job.dataset     = DL(job.datasetId,:);
            end
            experimenterMenu( job );

        case 4,
            % Look at current job set by the user
            displayExperimenterJob(job);
            experimenterMenu( job);

        case 5,
            % START EXECUTION
            fprintf('Executing Algorithms...\n');
            nAlgos  = length(job.algorithmId);

            % Load the configuration
            load ../GUI/config/config.mat;
            % Add the strategy library
            addpath('../GUI/lib');
            addpath('../Strategy');
            % Confirm the dataset
            load(strcat('../Data/', char(dataName(job.datasetId))));

            % Setting options
            opts.quiet_mode = 1; opts.display_interval = 500;
            opts.log_mode = 1; opts.mat_mode = 1;
            opts.analyze_mode = 1;
            opts.his = 0;
            opts.progress = 1;

            for i = 1:1:nAlgos
                % Construct input for the i-th algo, and execute

                msg = 'Running the ';
                msg = strcat(msg, num2str(i));
                msg = strcat(msg, '-th Algorithm ...');
                disp(msg);

                % Confirm the algorithm
                strategy_name = char(algorithmName(job.algorithmId(i)));
                % Confirm the parameters
                inputs = (job.parameters(i,:));

                strategy_fun = [strategy_name '(1, data, inputs, opts)'];
                [stats(:, i), cumprod_ret(:,i), daily_ret(:,i), daily_portfolio{i}] = eval(strategy_fun);
                returns = daily_ret - 1;
                portfolio = daily_portfolio;
                chosenStrategy = algorithmList{job.algorithmId};
            end

            % Construct the results data structure
            r.daily_ret             = daily_ret-1;
            r.cum_ret               = cumprod(daily_ret);
            r.daily_portfolio       = daily_portfolio;
            r.selectedAlgorithms    = job.algorithmId;
            r.dataFrequency         = dataFrequency(job.datasetId);
            displayTable2(r);
            resultManager2(r,job);

        case 6,
            disp('Exiting Experimenter --> to Home');
            homeMenu;

        otherwise,
            disp('ERROR: Please enter a valid input');
            experimenterMenu;
    end
end

%% Menu for Selecting Dataset
function [ dataId ] = dataSelectMenu(dataList )
% Display algorithm lsit, and ask user to select the method

    menuId = 5; % DataList
    displayMenu(menuId);
    numberOfDatasets = length(dataList);

    prompt  = strcat('Please enter your choice of Dataset (1-', num2str(numberOfDatasets));
    prompt  = strcat(prompt, '):');
    choice  = input(prompt);

    if choice >= 1 && choice <= numberOfDatasets
        dataId  = choice;
        msg     = strcat('You have selected: ', dataList(dataId));
        fprintf('\n');
        disp(char(msg));
        disp('Returning back to Algorithm Analyser...');
    else
        disp('ERROR: Please enter a valid input');
        dataId = dataSelectMenu( dataList );
    end
end

function [ job ] = selectAlgorithms( job )

    load ../GUI/config/config.mat;
    menuId = 9; % AlgorithmList for Experimenter
    displayMenu(menuId);
    numberOfAlgorithms  = length(algorithmList);
    numberSelected      = length(job.algorithmId);
    disp('Currently Selected Algorithms:');
    disp(job.algorithmId);

    msg         = 'Please Enter the Strategy IDs that you want to compare (1:';
    msg         = strcat(msg, num2str(numberOfAlgorithms));
    msg         = strcat(msg, ')');
    msg         = strcat(msg,'\nEg. To compare Strategies 1, 3, 4, 5, type: [1, 3, 4, 5]\n');
    selected    = input(msg);

    job.algorithmId = sort(selected);
    job.algorithmId = unique(selected);
    job.parameters  = {};
    % Add parameters
    numberSelected = length(job.algorithmId);
    for i = 1:1:numberSelected
        job.parameters(i,:) = defaultParameters(job.algorithmId(i),:);
    end
end

function [ ] = displayExperimenterJob( job )

    load ../GUI/config/config.mat;
    fprintf('\n');
    fprintf('________________________________________');
    fprintf('\n\t\t\t CURRENT JOB \n');

    nAlgos = length(job.algorithmId);
    fprintf('Algorithms: \n');
    fprintf('______________\n');
    for i = 1:1:nAlgos
        id = job.algorithmId(i);
        msg = strcat(num2str(i), '.');
        msg = strcat(msg,algorithmList{id});
        fprintf(msg);
        fprintf('\n');
        displayParameters(job, i);
        fprintf('______________\n');
    end

    fprintf('\n');
    fprintf(strcat('Dataset: \t', char(dataList(job.datasetId))));
    fprintf('\n');

    fprintf('________________________________________');
    fprintf('\n');
end

function [ ] = displayParameters( job, index )
    load ../GUI/config/config.mat;

    names           = algorithmParameters(job.algorithmId(index),:);
    numParameters   = length(cell2mat((job.parameters(index,:))));

    for i = 1:1:numParameters
        line = char(names(i));
        line = strcat(line, '\t=\t');
        line = strcat(line, num2str(cell2mat(job.parameters(index,i))));
        line = strcat(line, '\n');
        fprintf(line);
    end
end

function [ parameters ] = selectParameters(algoId)
% The interface for the user to input parameters for the algorithm
% currently selected.

    load ../GUI/config/config.mat;
    disp('________________________________________________________');
    disp('Please Enter the new values of parameters for Algorithm:');
    algo    = char(algorithmList(algoId));
    fprintf(strcat(strcat('Algorithm:\t',algo), '\n'));

    defaultParameters   = cell2mat(defaultParameters(algoId,:));
    numberOfParmeters   = length(defaultParameters);

    for i = 1:1:numberOfParmeters

        parameterName   = algorithmParameters(algoId, i);
        parameterName   = parameterName{1};

        msg = strcat(num2str(i), ')');
        msg = strcat(msg, parameterName);
        msg = strcat(msg, '(current value = ');
        msg = strcat(msg, num2str(defaultParameters(i)));
        msg = strcat(msg, '): ');

        parameters(i) = input(msg);

    end
    disp('________________________________________________________');
end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
clear all
close all
clc

% csidata;
% csi2batch;
getbatchdata;
gettestbatchdata;
fprintf(1,'finished');
%%
load batchdatalong;
[numcases numdims numbatches]=size(batchdata);
maxepoch=4; %In the Science paper we use maxepoch=50, but it works just fine.
numhid=100; numpen=50; numpen2=30; numopen=30;

fprintf(1,'Pretraining Layer 1 with RBM: %d-%d \n',numdims,numhid);
restart=1;
rbm;
hidrecbiases=hidbiases;
save mnistvh vishid hidrecbiases visbiases;
%%

fprintf(1,'\nPretraining Layer 2 with RBM: %d-%d \n',numhid,numpen);
batchdata=batchposhidprobs;
numhid=numpen;
restart=1;
rbm;
hidpen=vishid; penrecbiases=hidbiases; hidgenbiases=visbiases;
save mnisthp hidpen penrecbiases hidgenbiases;
%%
fprintf(1,'\nPretraining Layer 3 with RBM: %d-%d \n',numpen,numpen2);
batchdata=batchposhidprobs;
numhid=numpen2;
restart=1;
rbm;
hidpen2=vishid; penrecbiases2=hidbiases; hidgenbiases2=visbiases;
save mnisthp2 hidpen2 penrecbiases2 hidgenbiases2;
%%
fprintf(1,'\nPretraining Layer 4 with RBM: %d-%d \n',numpen2,numopen);
batchdata=batchposhidprobs;
numhid=numopen;
restart=1;
rbmhidlinear;
hidtop=vishid; toprecbiases=hidbiases; topgenbiases=visbiases;
save mnistpo hidtop toprecbiases topgenbiases;
% clear all
fprintf(1,'finished');
%%
load batchdata55;
load testbatchdata55;
% load testbatchdata44;
backprop;
% testerror('testbatchdata13');

repo: DeepFi


text:
% prune isolated, noisy neurons

f_r2d = @(r) sqrt(2 * (1-r));
f_d2r = @(d) 1- d.^2 / 2;
f_2r = @(r) f_d2r(f_r2d(r) * 2);

%%
FishClusterList = {11, 'all_thr50'};

stimCase = 2;
flag_truncate_after_transition = 1;

rx_f = 0.85;
q_f =  1/3; % final quantile

iFish = FishClusterList{1};
clusterName = FishClusterList{2};
if getappdata(hfig, 'i_fish') ~= iFish
    f.LoadFullFish(hfig,iFish);
end

[cIXt,gIXt,numKt]=load_cluster(hfig, clusterName);
cgk0 = {cIXt,gIXt,numKt};
cgk_current = cgk0;
push_cgk(hfig, cgk_current);

setStimRange(hfig, f, s2stimStr(stimCase));
if flag_truncate_after_transition == 1
    updateTRange_remove_transition(hfig,f);  % truncate after transition
end

%% step 1, norm prune
p_var = 0.95;
norm_th = 0.85;

tic
display('step 1: prune by norm (PCA)');
[cgk1, cgk1_rest] = cgk_norm_prune(hfig, norm_th);
toc

push_cgk(hfig, cgk1);
% push_cgk(hfig, cgk1_rest);

cgk_current = cgk1;

%% step 2, prune isolated neurons
p_var = 0.95;
% rx2 = f_d2r(f_r2d(rx_f) * sqrt(2)) % triangle inequality (conservative)
rx2 = rx_f .^ 2 % high dim asymptotic
nc_th = 8;
nc_bin = 1e3;
niter = 20;

tic
display('step2, prune isolated neurons');

[cgk2, cgk2_rest] = cgk_isolated_prune_dist(hfig, rx2, nc_th, niter);

push_cgk_k10(hfig, cgk2_rest);

push_cgk_k10(hfig, cgk2);
cgk_current = cgk2;

%% step 3-pre, parameter search
ncutoff = 4;
qx = q_f;
nc_th = 8;
niter = 20;

flag_linkage_method = 'complete';

% rx2 = f_d2r(f_r2d(rx_f) * 2) % center rcos threshold
% rx2 = rx_f .^ 2 % high dim asymptotic
d1 = 1 - f_d2r(f_r2d(rx_f))
d2 = min(2, 1 - f_d2r(f_r2d(rx_f) * 4))
rx = rx_f;
cutoff_dist_ls = linspace(d1, d2, ncutoff);
nC_prune_ls = zeros(1, ncutoff);
ng_ls = zeros(1, ncutoff);
ng0_ls = zeros(1, ncutoff); % raw result from linkage
cgk3_ls = cell(1, ncutoff);
cgk3_rest_ls = cell(1, ncutoff);

M = getappdata(hfig, 'M');
[nC, nT] = size(M);
M = zscore(M, 0, 2);
[pca_dim_M, k_var_M, U_M, S_M, V_M] = PCA_MP_fit(M, p_var, 0);

Mp = U_M(:, 1:k_var_M) * S_M(1:k_var_M, 1:k_var_M);
nTp = size(Mp, 2);

tic
display('linkage parameter search')
for i_c = 1 : ncutoff
    cutoff_dist = cutoff_dist_ls(i_c);
    switch flag_linkage_method
        case 'center'
            Z_link = linkage_cent_cos_v2(Mp);
        case 'average'
            Z_link = linkage(Mp, 'average', 'cosine');
        case  'single'
            Z_link = linkage(Mp, 'single', 'cosine');
        case  'complete'
            Z_link = linkage(Mp, 'complete', 'cosine');
        case 'centroid'
            norm_Mp = normM(Mp);
            Mp_norm = Mp ./ repmat(norm_Mp, [1, nTp]);
            %         d_scale = 2 * median(norm_Mp);
            %         cutoff_dist = d_scale * cutoff_dist;
%             cutoff_dist = sqrt(2 * cutoff_dist);
            Z_link = linkage(Mp_norm, 'centroid', 'euclidean');
    end
c_link_temp = cluster_u(Z_link, 'cutoff', cutoff_dist, 'criterion', 'distance');
cgk_current_temp = cgk_current;
cgk_current_temp{2} = c_link_temp;
cgk_current_temp{3} = max(c_link_temp); % not yet pushed to GUI

[cgk3_ls{i_c}, cgk3_rest_ls{i_c}] = cgkMp_rx_qx_prune(cgk_current_temp, Mp, rx, qx, nc_th, niter);

nC_prune_ls(i_c) = length(cgk3_rest_ls{i_c}{1});
ng0_ls(i_c) = length(unique(cgk_current_temp{2}));
ng_ls(i_c) = length(unique(cgk3_ls{i_c}{2}));
end
toc

%
x_ls = cutoff_dist_ls;
figure;
subplot(211);
plot(x_ls, ng0_ls, 'bo-');
hold on;
plot(x_ls, ng_ls, 'ro-');
hold off;
legend({'linkage', 'pruned'});
title('number of clusteres')
subplot(212);
plot(x_ls, nC_prune_ls, 'bo-');
title('pruned cells')

%%
save(fullfile(working_data_dir, 'step_prune_data_temp_01.mat'), 'cgk2', 'cgk3_ls', 'cgk3_rest_ls');

% %% step 3, linkage clustering
% p_var = 0.95;
% rx2 = f_d2r(f_r2d(rx_0) * 2) % center rcos threshold
% rx = rx_f;
% qx = q_f;
% nc_th = 8;
% niter = 20;
%
% cutoff_dist = 1 - rx2;
% flag_linkage_method = 'complete';
%
% tic
% display('step3, complete linkage');
%
% M = getappdata(hfig, 'M');
% [nC, nT] = size(M);
% M = zscore(M, 0, 2);
% [pca_dim_M, k_var_M, U_M, S_M, V_M] = PCA_MP_fit(M, p_var, 0);
%
% Mp = U_M(:, 1:k_var_M) * S_M(1:k_var_M, 1:k_var_M);
% nTp = size(Mp, 2);
%
%
% switch flag_linkage_method
%         case 'center'
%             Z_link = linkage_cent_cos_v2(Mp);
%         case 'average'
%             Z_link = linkage(Mp, 'average', 'cosine');
%         case  'single'
%             Z_link = linkage(Mp, 'single', 'cosine');
%         case  'complete'
%             Z_link = linkage(Mp, 'complete', 'cosine');
%         case 'centroid'
%             norm_Mp = normM(Mp);
%             Mp_norm = Mp ./ repmat(norm_Mp, [1, nTp]);
%             %         d_scale = 2 * median(norm_Mp);
%             %         cutoff_dist = d_scale * cutoff_dist;
% %             cutoff_dist = sqrt(2 * cutoff_dist);
%             Z_link = linkage(Mp_norm, 'centroid', 'euclidean');
% end
% c_link = cluster_u(Z_link, 'cutoff', cutoff_dist, 'criterion', 'distance');
% cgk_current{2} = c_link;
% cgk_current{3} = max(c_link); % not yet pushed to GUI
%
% [cgk3, cgk3_rest] = cgkMp_rx_qx_prune(cgk_current, Mp, rx, qx, nc_th, niter);
%
%
% push_cgk_k10(hfig, cgk3_rest);
%
% push_cgk(hfig, cgk3);
% cgk_current = cgk3;

%%

%
% %%
% push_cgk(hfig, cgk_link);
% rx = 0.9;
% rq = 0.30;
% nsize_th = 5;
%
% M = getappdata(hfig, 'M');
% M = zscore(M, 0, 2);
% [nC, nT] = size(M);
% p_var = 0.95;
% flag_plot = 0;
% [pca_dim_M, k_var_M, U_M, S_M, V_M] = PCA_MP_fit(M, p_var, flag_plot);
% Mp = U_M(:, 1:k_var_M) * S_M(1:k_var_M, 1:k_var_M);
% nTp = size(Mp, 2);
%
% gls = sort(unique(cgk_link{2}), 'ascend');
% ngc = length(gls);
% tf_core = zeros(nC, 1);
% for i = 1 : ngc
%     tf_gc_temp = cgk_link{2} == gls(i);
%     id_gc_temp = find(tf_gc_temp);
%     mM_temp = mean(Mp(tf_gc_temp, :), 1);
%     rc_temp = rcosM(Mp(tf_gc_temp, :), mM_temp);
%     rc_temp_sort = sort(rc_temp, 'descend');
%     i_th_temp = round(rq * length(rc_temp) + 1);
%     i_th_temp = max(i_th_temp, 1);
%     i_th_temp = min(i_th_temp, length(rc_temp));
%     rq_th = rc_temp_sort(i_th_temp);
%     if i_th_temp >= nsize_th
%         r_th_temp = max(rx, rq_th);
%     else
%         r_th_temp = rx;
%     end
%     if sum(rc_temp > r_th_temp) >= nsize_th
%         tf_core(id_gc_temp(rc_temp > rx)) = 1;
%     end
% end
%
% tf_core = logical(tf_core);
% cgk_core = cgk_link;
% cgk_core{1} = cgk_core{1}(tf_core);
% cgk_core{2} = cgk_core{2}(tf_core);
%
% push_cgk(hfig, cgk_core);

% %%
%
%
% ns=length(s);
% % plot eigenvalues
% figure;
% set(gcf,'position',[1117         536         1120         420]);
% subplot(1,3,1);
% scatter(1:ns,s,'ro');
% title('eigenvalues')
%
% subplot(1,3,2);
% scatter(1:ns,log(s),'ro');
% title('log eigenvalues')
%
% sum_res=flipud(cumsum(flipud(s.^2)));
% total_var=sum_res(1);
% sum_res=[sum_res(2:end);0]/total_var;
% subplot(1,3,3);
% scatter(1:ns,sum_res*100,'ro');
% title('residue variance (percentage)')
%
% nX=4;
% X=V(:,1:nX);
% Y=X;
% % for i=1:nC
% %     Y(i,:)=Y(i,:)/norm(X(i,:),2);
% % end
%
% markerSize=16;
%
% u1=Y(:,1);
% u2=Y(:,2);
% u3=Y(:,3);
%
% if ismember(2,dims)
%     figure;
%     set(gcf,'position',[1131   43   560   420]);
%     scatter(u1,u2,markerSize*ones(1,nC), cmap(gIX,:),'fill');
%     axis equal;
% end;
%
% if ismember(3,dims)
%     figure;
%     set(gcf,'position',[87   386   733   548]);
%     scatter3(u1,u2,u3,markerSize*ones(1,nC), cmap(gIX,:),'fill');
%     axis equal;
% end;
%
%
% %%
% nkmean = 5;
% idx = kmeans(Y, nkmean, 'replicate', 20);
%
% cmap_kmean = cmap_cluster(hfig, nkmean);
%
% if ismember(2,dims)
%     figure;
%     set(gcf,'position',[1131   43   560   420]);
%     scatter(u1,u2,markerSize*ones(1,nC), cmap_kmean(idx,:),'fill');
%     axis equal;
% end;
%
% if ismember(3,dims)
%     figure;
%     set(gcf,'position',[87   386   733   548]);
%     scatter3(u1,u2,u3,markerSize*ones(1,nC), cmap_kmean(idx,:),'fill');
%     axis equal;
% end;
%
% % cgk0_temp = get_cgk(hfig);
% cgk1_temp = cgk0_temp;
% cgk1_temp{2} = idx;
% cgk1_temp{3} = max(idx);
% push_cgk(hfig, cgk1_temp);

%%
% cgk_s12_diff = cgk_s1;
% [cgk_s1{1}, id_sort_temp] = sort(cgk_s1{1});
% cgk_s1{2} = cgk_s1{2}(id_sort_temp);
% [cgk_s2{1}, id_sort_temp] = sort(cgk_s2{1});
% cgk_s2{2} = cgk_s2{2}(id_sort_temp);
%
% tf_diff = cgk_s1{2} ~= cgk_s2{2};
% cgk_s12_diff{1} = cgk_s12_diff{1}(tf_diff);
% cgk_s12_diff{2} = cgk_s12_diff{2}(tf_diff);
%
% push_cgk(hfig, cgk_s12_diff);

% %% parameter  search of inconsistency
% ncutoff = 60;
% q_f =  1/3; % final quantile
% cutoff_dist_ls = linspace(0.2, 1, ncutoff);
% nC_rest1_ls = zeros(1, ncutoff);
% ngc1_ls = zeros(1, ncutoff);
% nC_rest2_ls = zeros(1, ncutoff);
% ngc2_ls = zeros(1, ncutoff);
% cgk_link1_ls = cell(1, ncutoff);
% cgk_link2_ls = cell(1, ncutoff);
% cgk_link1_rest_ls = cell(1, ncutoff);
% cgk_link2_rest_ls = cell(1, ncutoff);
%
% tic
% display('linkage clustering')
% for i_p = 1 : ncutoff
%     cutoff_dist = cutoff_dist_ls(i_p);
%     nc_size_th = nc_th_f;
%     p_var = 0.95;
% %     flag_linkage_method = 'average';
%     flag_linkage_method = 'centroid';
%     % flag_linkage_method = 'center';
%     % flag_linkage_method = 'complete';
%
%     push_cgk(hfig, cgk3_temp);
%     M = getappdata(hfig, 'M');
%     M = zscore(M, 0, 2);
%     [nC, nT] = size(M);
%
%     [pca_dim_M, k_var_M, U_M, S_M, V_M] = PCA_MP_fit(M, p_var, 0);
%     Mp = U_M(:, 1:k_var_M) * S_M(1:k_var_M, 1:k_var_M);
%     nTp = size(Mp, 2);
%
%     switch flag_linkage_method
%         case 'center'
%             Z_link = linkage_cent_cos_v2(Mp);
%         case 'average'
%             Z_link = linkage(Mp, 'average', 'cosine');
%         case  'single'
%             Z_link = linkage(Mp, 'single', 'cosine');
%         case  'complete'
%             Z_link = linkage(Mp, 'complete', 'cosine');
%         case 'centroid'
%             norm_Mp = normM(Mp);
%             Mp_norm = Mp ./ repmat(norm_Mp, [1, nTp]);
%             %         d_scale = 2 * median(norm_Mp);
%             %         cutoff_dist = d_scale * cutoff_dist;
%             cutoff_dist = sqrt(2 * cutoff_dist);
%             Z_link = linkage(Mp_norm, 'centroid', 'euclidean');
%     end
%     % c_link = cluster(Z_link, 'cutoff', cutoff_dist, 'criterion', 'distance');
%     c_link = cluster_u(Z_link, 'cutoff', cutoff_dist, 'criterion', 'inconsistent');
%
%     nc_link = max(c_link);
%     nsize_c_link = zeros(nc_link, 1);
%     for i = 1 : nc_link
%         nsize_c_link(i) = sum(c_link == i);
%     end
%     id_c_size_th = find(nsize_c_link >= nc_size_th);
%     nid_c = length(id_c_size_th);
%
%     % update loop parameter
%     ngc1_ls(i_p) = nid_c;
%     nC_rest1_ls(i_p) = length(cgk3_temp{1}) - sum(nsize_c_link(id_c_size_th));
%
%
%     Mcent_link = zeros(nid_c, nTp);
%     nsize_c_link_new = zeros(nid_c, 1);
%     for i = 1 : nid_c
%         Mcent_link(i, :) = mean(Mp(c_link == id_c_size_th(i), :), 1);
%         nsize_c_link_new(i) = sum(c_link == id_c_size_th(i));
%     end
%
%     if size(Mcent_link, 1) > 1
%         switch flag_linkage_method
%             case 'center'
%                 D_temp = pdist(Mcent_link, 'cosine');
%                 tree_temp = linkage_cent_cos_v2(Mcent_link, nsize_c_link_new);
%             case 'average'
%                 D_temp = pdist(Mcent_link, 'cosine');
%                 tree_temp = linkage(Mcent_link, 'average', 'cosine');
%             case  'single'
%                 D_temp = pdist(Mcent_link, 'cosine');
%                 tree_temp = linkage(Mcent_link, 'single', 'cosine');
%             case  'complete'
%                 D_temp = pdist(Mcent_link, 'cosine');
%                 tree_temp = linkage(Mcent_link, 'complete', 'cosine');
%             case 'centroid'
%                 D_temp = pdist(Mcent_link, 'euclidean');
%                 tree_temp = linkage(Mcent_link, 'centroid', 'euclidean');
%         end
%         leafOrder_temp = optimalleaforder(tree_temp, D_temp);
%     else
%         leafOrder_temp = 1;
%     end
%
%     cgk_link_rest = cgk3_temp;
%     cgk_link = cgk3_temp;
%     cgk_link{2}(:) = 0;
%     for i = 1 : nid_c
%         cgk_link{2}(c_link == id_c_size_th(leafOrder_temp(i))) = i;
%     end
%     cgk_link_rest{1} = cgk_link_rest{1}(cgk_link{2} == 0);
%     cgk_link_rest{2} = cgk_link_rest{2}(cgk_link{2} == 0);
%     cgk_link{1} = cgk_link{1}(cgk_link{2} > 0);
%     cgk_link{2} = cgk_link{2}(cgk_link{2} > 0);
%     cgk_link{3} = nid_c;
%
%     % update loop parameter
%     cgk_link1_ls{i_p} = cgk_link;
%     cgk_link1_rest_ls{i_p} = cgk_link_rest;
%
%     if length(cgk_link{1}) > 1
%         % core selection
%         push_cgk(hfig, cgk_link);
%         M = getappdata(hfig, 'M');
%         M = zscore(M, 0, 2);
%         [nC, nT] = size(M);
%         [pca_dim_M, k_var_M, U_M, S_M, V_M] = PCA_MP_fit(M, p_var, 0);
%         Mp = U_M(:, 1:k_var_M) * S_M(1:k_var_M, 1:k_var_M);
%         nTp = size(Mp, 2);
%
%         nc_link = cgk_link{3};
%         nsize_c_link = zeros(nc_link, 1);
%         tf_core = false(length(cgk_link{1}), 1);
%         for i = 1 : nc_link
%             mM_temp = mean(Mp(cgk_link{2} == i, :), 1);
%             id_cluster_temp = find(cgk_link{2} == i);
%             rcos_v_temp = rcosM(Mp(cgk_link{2} == i, :), mM_temp);
%             r_sort_temp = sort(rcos_v_temp, 'descend');
%             i_th_temp = round(q_f * length(r_sort_temp) + 1);
%             i_th_temp = max(i_th_temp, 1);
%             i_th_temp = min(i_th_temp, length(r_sort_temp));
%             rq_th = r_sort_temp(i_th_temp);
%             if i_th_temp >= nc_th_f
%                 r_th_temp = max(rq_th, rx_f);
%             else
%                 r_th_temp = rx_f;
%             end
%             nsize_c_link(i) = sum(rcos_v_temp > r_th_temp);
%             if nsize_c_link(i) >= nc_th_f
%                 tf_core(id_cluster_temp(rcos_v_temp > r_th_temp)) = true;
%             end
%         end
%         id_c_size_th = find(nsize_c_link >= nc_th_f);
%         nid_c = length(id_c_size_th);
%         Mcent_link = zeros(nid_c, nTp);
%         for i = 1 : nid_c
%             tf_valid_temp = cgk_link{2} == id_c_size_th(i) & tf_core;
%             Mcent_link(i, :) = mean(Mp(tf_valid_temp, :), 1);
%         end
%         if size(Mcent_link, 1) > 1
%             flag_linkage_method_step2 = 'average';
%             switch flag_linkage_method_step2
%                 case 'center'
%                     D_temp = pdist(Mcent_link, 'cosine');
%                     tree_temp = linkage_cent_cos_v2(Mcent_link, nsize_c_link_new);
%                 case 'average'
%                     D_temp = pdist(Mcent_link, 'cosine');
%                     tree_temp = linkage(Mcent_link, 'average', 'cosine');
%                 case  'single'
%                     D_temp = pdist(Mcent_link, 'cosine');
%                     tree_temp = linkage(Mcent_link, 'single', 'cosine');
%                 case  'complete'
%                     D_temp = pdist(Mcent_link, 'cosine');
%                     tree_temp = linkage(Mcent_link, 'complete', 'cosine');
%                 case 'centroid'
%                     D_temp = pdist(Mcent_link, 'euclidean');
%                     tree_temp = linkage(Mcent_link, 'centroid', 'euclidean');
%             end
%             leafOrder_temp = optimalleaforder(tree_temp, D_temp);
%         else
%             leafOrder_temp = 1;
%         end
%
%         cgk_link2_rest = cgk_link;
%         cgk_link2 = cgk_link;
%         cgk_link2{2}(:) = 0;
%         for i = 1 : nid_c
%             cgk_link2{2}(cgk_link{2} == id_c_size_th(leafOrder_temp(i))) = i;
%         end
%         cgk_link2_rest{1} = cgk_link2_rest{1}(cgk_link2{2} == 0);
%         cgk_link2_rest{2} = cgk_link2_rest{2}(cgk_link2{2} == 0);
%         cgk_link2{1} = cgk_link2{1}(cgk_link2{2} > 0);
%         cgk_link2{2} = cgk_link2{2}(cgk_link2{2} > 0);
%         cgk_link2{3} = nid_c;
%
%         % update loop parameter
%         cgk_link2_ls{i_p} = cgk_link2;
%         cgk_link2_rest_ls{i_p} = cgk_link2_rest;
%
%         ngc2_ls(i_p) = nid_c;
%         nC_rest2_ls(i_p) = length(cgk_link{1}) - sum(nsize_c_link(id_c_size_th));
%     else
%         % update loop parameter
%         cgk_link2_ls{i_p} = nan;
%         cgk_link2_rest_ls{i_p} = nan;
%
%         ngc2_ls(i_p) = 0;
%         nC_rest2_ls(i_p) = 0;
%     end
% end
%
% toc
%
% %
% figure;
% subplot(211);
% plot(cutoff_dist_ls, ngc1_ls, 'bo-');
% hold on;
% plot(cutoff_dist_ls, ngc2_ls, 'ro-');
% hold off;
% title('number of clusteres')
% subplot(212);
% plot(cutoff_dist_ls, nC_rest1_ls, 'bo-');
% hold on;
% plot(cutoff_dist_ls, nC_rest2_ls, 'ro-');
% hold off;
% title('pruned cells')
%
%
% beep
% beep
%
% %%
% push_cgk(hfig, cgk_link2_ls{35});

author: xiuyechen

repo: FishExplorer


text:
%QUESTION NO:3

% d)Design and Train a feedforward network for the following problem:
%   Addition: Consider a 4-input and 3-output problem, where the output
%   should be the result of the sum of two 2-bit input numbers.

clear
inp=[0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1;0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1;...
    0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1;0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1];
out=[0 0 0 0 0 0 0 1 0 0 1 1 0 1 1 1; 0 0 1 1 0 1 1 0 1 1 0 0 1 0 0 1; ...
     0 1 0 1 1 0 1 0 0 1 0 1 1 0 1 0];
network=newff([0 1;0 1; 0 1; 0 1],[6 3],{'logsig','logsig'});
network=init(network);
y=sim(network,inp);
network.trainParam.epochs = 500;
network=train(network,inp,out);
y=sim(network,inp);
Layer1_Weights=network.iw{1};
Layer1_Bias=network.b{1};
Layer2_Weights=network.lw{2};
Layer2_Bias=network.b{2};
Layer1_Weights
Layer1_Bias
Layer2_Weights
Layer2_Bias
Actual_Desired=[y' out'];
Actual_Desired


text:
function [ out ] = func_projection_filterbank( dat, w )

CSP_W=w;

for ii=1:length(CSP_W)
    SMT{ii}=func_projection(dat{ii}, CSP_W{ii});
end

out=SMT;
end

author: PatternRecognition

repo: OpenBMI


text:
classdef Momentum < handle

    properties (Access = private)
        iterator
        designVariable
        momentumParameter
        tOld
        tOldOld
    end

    methods (Access = public)

        function obj = Momentum(cParams)
            obj.designVariable = cParams.designVariable;
            obj.iterator       = cParams.iterator;
        end

        function apply(obj)
            obj.computeMomentumParameter();
            obj.addIntertialTerm();
        end

    end

    methods (Access = private)

        function computeMomentumParameter(obj)
            i = obj.iterator.value;
            if i < 3
                obj.tOld = 1;
                obj.tOldOld = 1;
                beta = 1;
            else
                t = (1+sqrt(1+4*obj.tOld^2))/2;
                beta = (obj.tOldOld-1)/t;  % beta = i/(i+3);
                obj.tOldOld = obj.tOld;
                obj.tOld = t;
            end
            obj.momentumParameter = beta;
        end

        function addIntertialTerm(obj)
            beta    = obj.momentumParameter;
            xOldOld = obj.designVariable.valueOldOld;
            xOld    = obj.designVariable.valueOld;
            x = xOld + beta*(xOld - xOldOld);
            obj.designVariable.value = x;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function fx = p03_fun ( n, x )

%*****************************************************************************80
%
%% P03_FUN evaluates the integrand for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gwynne Evans,
%    Practical Numerical Integration,
%    Wiley, 1993,
%    ISBN: 047193898X,
%    LC: QA299.3E93.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(2,N), the evaluation points.
%
%    Output, real FX(N,1), the integrand values.
%
  fx(1:n,1) = 1.0 ./ sqrt ( 2.0 - x(1,1:n) - x(2,1:n) );

  return
end


text:
function x = quasigeometric_cdf_inv ( cdf, a, b )

%*****************************************************************************80
%
%% QUASIGEOMETRIC_CDF_INV inverts the Quasigeometric CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%    0.0 <= CDF <= 1.0
%
%    Input, real A, the probability of 0 successes.
%    0.0 <= A <= 1.0.
%
%    Input, real B, the depreciation constant.
%    0.0 <= B < 1.0.
%
%    Output, integer X, the corresponding value of X.
%
  if ( cdf < 0.0 | 1.0 < cdf )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'QUASIGEOMETRIC_CDF_INV - Fatal error!\n' );
    fprintf ( 1, '  CDF < 0 or 1 < CDF.\n' );
    error ( 'QUASIGEOMETRIC_CDF_INV - Fatal error!' );
  end

  if ( cdf < a )
    x = 0;
  elseif ( b == 0.0 )
    x = 1;
  else
    x = 1 + floor ( ( log ( 1.0 - cdf ) - log ( 1.0 - a ) ) / log ( b ) );
  end

  return
end


text:
function g = sigmoid(z)
%SIGMOID Compute sigmoid functoon
%   J = SIGMOID(z) computes the sigmoid of z.

% You need to return the following variables correctly
g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the sigmoid of each value of z (z can be a matrix,
%               vector or scalar).

g = 1 ./ (1 + (1 ./ exp(z)));

% =============================================================

end

author: sfvsfv

repo: Mathematical-modeling


text:
function value = r4_si ( x )

%*****************************************************************************80
%
%% R4_SI evaluates the sine integral Si of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VASLUE, the sine integral Si evaluated at X.
%
  persistent nsi
  persistent sics
  persistent xsml

  if ( isempty ( nsi ) )

    sics = [ ...
      -0.1315646598184841929, ...
      -0.2776578526973601892, ...
       0.0354414054866659180, ...
      -0.0025631631447933978, ...
       0.0001162365390497009, ...
      -0.0000035904327241606, ...
       0.0000000802342123706, ...
      -0.0000000013562997693, ...
       0.0000000000179440722, ...
      -0.0000000000001908387, ...
       0.0000000000000016670, ...
      -0.0000000000000000122 ]';

    nsi = r4_inits ( sics, 12, 0.1 * r4_mach ( 3 ) );
    xsml = sqrt ( eps );

  end

  absx = abs ( x );

  if ( absx < xsml )

    value = x;

  elseif ( absx <= 4.0 )

    value = x * ( 0.75 + r4_csevl ( ( x * x - 8.0 ) * 0.125, sics, nsi ) );

  else

    [ f, g ] = r4_sifg ( absx );
    cosx = cos ( absx );

    if ( x < 0.0 )
      value = - 0.5 * pi + f * cosx + g * sin ( x );
    else
      value = 0.5 * pi - f * cosx - g * sin ( x );
    end

  end

  return
end


text:
function Y = multiherm(X)
% Returns the Hermitian parts of the matrices in a 3D array
%
% function Y = multiherm(X)
%
% Y is a 3D array the same size as X. Each slice Y(:, :, i) is the
% Hermitian part of the slice X(:, :, i), that is,
%
%   Y(:, :, i) = .5*(X(:, :, i) + X(:, :, i)')
%
% See also: multisym multiskew multiskewh multihconj multiprod multitransp multiscale

% This file is part of Manopt: www.manopt.org.
% Original author: Hiroyuki Sato, April 27, 2015.
% Contributors:
% Change log:

    Y = .5*(X + multihconj(X));

end

author: NicolasBoumal

repo: manopt


text:
function [wedge_data wedge_density] = wedge();

% wedge coordinates
Y = [-7.0 -4.35    0  3.0  0 -4.35 -7.0 -7.0];
Z = [-5.8 -5.8  -2.46    0  0   0    0  -5.8];
figure; line(Y,Z, 'LineStyle', '-', 'Marker', 'o');
disp('This wedge can be modeled as a rectangle and a triangle, as shown below')
Y = [-7.0 -4.35   -4.35 -7.0 -7.0];
Z = [-5.8 -5.8  0    0  -5.8];
figure; line(Y,Z, 'Color', 'b', 'LineStyle', '-', 'Marker', 'o');
Y = [-4.35  3.0  -4.35 -4.35];
Z = [-5.8    0     0   -5.8];
hold on; line(Y,Z, 'Color', 'r', 'LineStyle', '-', 'Marker', 'o'); axis equal;
xlabel('Y/cm'); ylabel('Z/cm'); title('Elekta Universal Wedge');

% From the NIST website, get the attenuation coefficients.
% Components of the wedge:
Pb = 0.96     %Fraction by weigth for Lead
Sb = 0.04     %Fraction by weight for Antimony

% Ref:
%http://physics.nist.gov/PhysRefData/XrayMassCoef/ElemTab/z82.html
%Energy       mu over rho        mu e n over rho
%      (MeV)      (cm2/g)     (cm2/g)
%____________________________________
%
Pb_data =   [1.00000E-03  5.210E+03  5.197E+03
   1.50000E-03  2.356E+03  2.344E+03
   2.00000E-03  1.285E+03  1.274E+03
   2.48400E-03  8.006E+02  7.895E+02
   2.48400E-03  1.397E+03  1.366E+03
   2.53429E-03  1.726E+03  1.682E+03
   2.58560E-03  1.944E+03  1.895E+03
   2.58560E-03  2.458E+03  2.390E+03
   3.00000E-03  1.965E+03  1.913E+03
   3.06640E-03  1.857E+03  1.808E+03
   3.06640E-03  2.146E+03  2.090E+03
   3.30130E-03  1.796E+03  1.748E+03
   3.55420E-03  1.496E+03  1.459E+03
   3.55420E-03  1.585E+03  1.546E+03
   3.69948E-03  1.442E+03  1.405E+03
   3.85070E-03  1.311E+03  1.279E+03
   3.85070E-03  1.368E+03  1.335E+03
   4.00000E-03  1.251E+03  1.221E+03
   5.00000E-03  7.304E+02  7.124E+02
   6.00000E-03  4.672E+02  4.546E+02
   8.00000E-03  2.287E+02  2.207E+02
   1.00000E-02  1.306E+02  1.247E+02
   1.30352E-02  6.701E+01  6.270E+01
   1.30352E-02  1.621E+02  1.291E+02
   1.50000E-02  1.116E+02  9.100E+01
   1.52000E-02  1.078E+02  8.807E+01
   1.52000E-02  1.485E+02  1.131E+02
   1.55269E-02  1.416E+02  1.083E+02
   1.58608E-02  1.344E+02  1.032E+02
   1.58608E-02  1.548E+02  1.180E+02
   2.00000E-02  8.636E+01  6.899E+01
   3.00000E-02  3.032E+01  2.536E+01
   4.00000E-02  1.436E+01  1.211E+01
   5.00000E-02  8.041E+00  6.740E+00
   6.00000E-02  5.021E+00  4.149E+00
   8.00000E-02  2.419E+00  1.916E+00
   8.80045E-02  1.910E+00  1.482E+00
   8.80045E-02  7.683E+00  2.160E+00
   1.00000E-01  5.549E+00  1.976E+00
   1.50000E-01  2.014E+00  1.056E+00
   2.00000E-01  9.985E-01  5.870E-01
   3.00000E-01  4.031E-01  2.455E-01
   4.00000E-01  2.323E-01  1.370E-01
   5.00000E-01  1.614E-01  9.128E-02
   6.00000E-01  1.248E-01  6.819E-02
   8.00000E-01  8.870E-02  4.644E-02
   1.00000E+00  7.102E-02  3.654E-02
   1.25000E+00  5.876E-02  2.988E-02
   1.50000E+00  5.222E-02  2.640E-02
   2.00000E+00  4.606E-02  2.360E-02
   3.00000E+00  4.234E-02  2.322E-02
   4.00000E+00  4.197E-02  2.449E-02
   5.00000E+00  4.272E-02  2.600E-02
   6.00000E+00  4.391E-02  2.744E-02
   8.00000E+00  4.675E-02  2.989E-02
   1.00000E+01  4.972E-02  3.181E-02
   1.50000E+01  5.658E-02  3.478E-02
   2.00000E+01  6.206E-02  3.595E-02];

%Antimony
%Z = 51
%
%ASCII format
%
%____________________________________
%
%     Energy       mu over rho        mu e n over rho
%      (MeV)      (cm2/g)     (cm2/g)
%____________________________________

Sb_data = [1.00000E-03  8.582E+03  8.568E+03
   1.50000E-03  3.491E+03  3.481E+03
   2.00000E-03  1.767E+03  1.759E+03
   3.00000E-03  6.536E+02  6.469E+02
   4.00000E-03  3.169E+02  3.113E+02
   4.13220E-03  2.918E+02  2.863E+02
   4.13220E-03  8.691E+02  8.252E+02
   4.25449E-03  8.308E+02  7.865E+02
   4.38040E-03  7.776E+02  7.392E+02
   4.38040E-03  1.050E+03  9.906E+02
   4.53657E-03  9.743E+02  9.178E+02
   4.69830E-03  8.939E+02  8.457E+02
   4.69830E-03  1.029E+03  9.725E+02
   5.00000E-03  8.846E+02  8.377E+02
   6.00000E-03  5.569E+02  5.305E+02
   8.00000E-03  2.631E+02  2.518E+02
   1.00000E-02  1.459E+02  1.396E+02
   1.50000E-02  4.923E+01  4.657E+01
   2.00000E-02  2.268E+01  2.105E+01
   3.00000E-02  7.631E+00  6.755E+00
   3.04912E-02  7.307E+00  6.452E+00
   3.04912E-02  4.073E+01  1.391E+01
   4.00000E-02  2.027E+01  9.789E+00
   5.00000E-02  1.120E+01  6.400E+00
   6.00000E-02  6.879E+00  4.311E+00
   8.00000E-02  3.176E+00  2.173E+00
   1.00000E-01  1.758E+00  1.237E+00
   1.50000E-01  6.361E-01  4.312E-01
   2.00000E-01  3.381E-01  2.084E-01
   3.00000E-01  1.677E-01  8.504E-02
   4.00000E-01  1.172E-01  5.288E-02
   5.00000E-01  9.453E-02  4.061E-02
   6.00000E-01  8.153E-02  3.465E-02
   8.00000E-01  6.670E-02  2.896E-02
   1.00000E+00  5.797E-02  2.608E-02
   1.25000E+00  5.086E-02  2.378E-02
   1.50000E+00  4.628E-02  2.230E-02
   2.00000E+00  4.105E-02  2.081E-02
   3.00000E+00  3.686E-02  2.043E-02
   4.00000E+00  3.567E-02  2.118E-02
   5.00000E+00  3.559E-02  2.219E-02
   6.00000E+00  3.598E-02  2.321E-02
   8.00000E+00  3.745E-02  2.509E-02
   1.00000E+01  3.921E-02  2.664E-02
   1.50000E+01  4.351E-02  2.920E-02
   2.00000E+01  4.704E-02  3.038E-02];

% put data on the same, unique intervals.
[A, s] = unique(Pb_data(:,1));
[B,I] = unique(Sb_data(:,1));
Pb_data = [B interp1(Pb_data(s,1), Pb_data(s,2), B) interp1(Pb_data(s,1), Pb_data(s,3), B)];
Sb_data = [B interp1(Sb_data(I,1), Sb_data(I,2), B) interp1(Sb_data(I,1), Sb_data(I,3), B)];

wedge_data = Pb*Pb_data + Sb*Sb_data;
figure; loglog(Pb_data(:,1), Pb_data(:,2), Sb_data(:,1), Sb_data(:,2), wedge_data(:,1), wedge_data(:,2));
wedge_density = 11.16;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ y, m, d, f] = frac_borrow_julian ( y, m, d, f )

%*****************************************************************************80
%
%% FRAC_BORROW_JULIAN borrows fractions from days in a Julian YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, integer M, integer D, real F,
%    a YMDF date.
%
%    Output, integer Y, integer M, integer D, real F,
%    a YMDF date.
%
  while ( f < 0.0 )

    f = f + 1.0;

    d = d - 1;

  end

  [ y, m, d ] = day_borrow_julian ( y, m, d );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function proposals = calcselective_searchForIm( input,ssconfig)

colorTypes=ssconfig.params.colorTypes;
simFunctionHandles=ssconfig.params.simFunctionHandles;
ks=ssconfig.params.ks;
im_width=ssconfig.params.imWidth;
sigma=ssconfig.params.sigma;
minBoxWidth =ssconfig.params.minBoxWidth;

if(isstr(input))
        im = im2uint8(imread(input));
else
        im = im2uint8(input); % Just to make input consistent
end
if(size(im, 3) == 1)
        im=repmat(im,[1,1,3]);
end
if ~isfield(ssconfig.params, 'imWidth')
	im_width = [];
        scale = 1;
else
        scale = size(im, 2) / im_width;
end

if scale ~= 1
        im = imresize(im, [NaN im_width]);
end
idx = 1;
for j = 1:length(ks)
        k = ks(j); % Segmentation threshold k
        minSize = k; % We set minSize = k
        for n = 1:length(colorTypes)
        	colorType = colorTypes{n};
                [boxesT{idx} blobIndIm blobBoxes hierarchy priorityT{idx}] = ...
                Image2HierarchicalGrouping(im, sigma, k, minSize, colorType, simFunctionHandles);
                idx = idx + 1;
        end
end
boxes = cat(1, boxesT{:}); % Concatenate boxes from all hierarchies
priority = cat(1, priorityT{:}); % Concatenate priorities

% Do pseudo random sorting as in paper
priority = priority .* rand(size(priority));
[priority sortIds] = sort(priority, 'ascend');
boxes = boxes(sortIds,:);

boxes = FilterBoxesWidth(boxes, minBoxWidth);
boxes = BoxRemoveDuplicates(boxes);

if scale ~= 1
	boxes = (boxes - 1) * scale + 1;
end

if(isfield(ssconfig.opts,'numProposals'))
        numProposals=ssconfig.opts.numProposals;
        if(size(boxes,1)>=numProposals)
        	boxes=boxes(1:numProposals,:);
        else
                fprintf('Only %d proposals were generated for input image \n',size(boxes,1));
        end
end

% reset boxes to xmin ymin xmax ymnx
boxes=[boxes(:,2) boxes(:,1) boxes(:,4) boxes(:,3)];
proposals.boxes=boxes;

author: Cloud-CV

repo: object-proposals


text:
function pde = HodgeLaplacianEdata1
%% HODGELAPLACIANEDATA1
%
% u = [cos(x), -sin(y)];
% sigma = -div u = sin(x) + cos(y);
% f = -grad div u + curl rot u = - Lap u;
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

pde.f = @f;
pde.u = @exactu;
pde.sigma = @exactsigma;
pde.gu = @exactu;
pde.gsigma = @exactsigma;
pde.gun = @gun;
pde.grotu = 0;

    function s = f(p)
    s = exactu(p);
    end

    function u = exactu(p)
    x = p(:,1); y = p(:,2);
    u = [cos(x), -sin(y)];
    end

    function s = exactsigma(p)
    x = p(:,1); y = p(:,2);
    s = sin(x) + cos(y);
    end

    function s = gun(p) % for unit square [0,1]^2
    s = zeros(size(p,1),1);
    x = p(:,1); y = p(:,2);
    u = exactu(p);
    leftbd = (abs(x)<eps);  % n = (-1,0);
    s(leftbd) = - u(leftbd,1);
    rightbd = (abs(x-1)<eps); % n = (1,0);
    s(rightbd) = u(rightbd,1);
    topbd = (abs(y-1)<eps);   % n = (0,1)
    s(topbd) = u(topbd,2);
    bottombd = (abs(y)<eps);% n = (0,-1)
    s(bottombd) = - u(bottombd,2);
    end

end

repo: ifem


text:
function hypercube_integrals_test ( )

%*****************************************************************************80
%
%% HYPERCUBE_INTEGRALS_TEST tests the HYPERCUBE_INTEGRALS library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HYPERCUBE_INTEGRALS_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the HYPERCUBE_INTEGRALS library.\n' );

  hypercube_integrals_test01 ( );
  hypercube_integrals_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HYPERCUBE_INTEGRALS_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function [eparams, errors] = msCalibrateEdgeClassifier(efeatures, adjlist, imsegs, eclassifier, labels, ncv)
% [eparams, errors] = calibrateEdgeClassifier(efeatures, adjlist, imsegs,
% eclassifier, ncv)

nimages = numel(imsegs);
for k = 1:ncv
    if ncv > 1
        testind = [(k-1)*nimages/ncv+1:k*nimages/ncv];
        trainind = setdiff([1:nimages], testind);
    else
        trainind = (1:nimages);
    end
    [edata{k}, elab{k}] = formatData(efeatures(trainind), adjlist(trainind), labels(trainind));
    econf{k} = test_boosted_dt_mc(eclassifier, edata{k});
%    econf{k} = 1 ./ (1+exp(-econf{k}));
end

for k = 1:ncv
    disp(['iter: ' num2str(k)])
    if ncv>1
        traink = setdiff([1:ncv], k);
    else
        traink = k;
    end
    eparams{k} = fminunc(@(x) objective(x, cat(1, econf{traink}), cat(1, elab{traink})), [-1 0], optimset('TolFun', 0.001));
end

for k = 1:ncv
    econf{k} = 1 ./ (1+exp(eparams{k}(1)*econf{k}+eparams{k}(2)));
end

elab = cat(1, elab{:});
econf = cat(1, econf{:});

eerror = mean((econf>0.5)~=elab);

econf2 = 1-abs(elab-econf);

ind1 = find(elab==0);
ind2 = find(elab==1);
px = [0.025:0.05:0.975];
f1 = ksdensity(econf(ind1), px, 'support', [0 1]);
f2 = ksdensity(econf(ind2), px, 'support', [0 1]);
fc = ksdensity(econf2, px, 'support', [0 1]);
%fc = fc;

errors.err = eerror;
errors.pneg = f1;
errors.ppos = f2;
errors.conf = fc;
errors.px = px;

medFS = 18;
bigFS = 20;

figure(1), hold on, plot(px, fc, 'y', 'LineWidth', 2);
%axis([0 1 0 1])
xlabel('Confidence in True Label', 'FontSize', medFS)
ylabel('Frequency', 'FontSize', medFS)
title('Same Label Confidence', 'FontSize', bigFS)
set(gca, 'FontSize', medFS)

figure(2), hold on, plot(px, f2 ./ (f1+f2), 'y', 'LineWidth', 2)
hold on, plot(px, px, '--k')
axis([0 1 0 1])
xlabel('Estimated Probability', 'FontSize', medFS)
ylabel('Empirical Probability', 'FontSize', medFS)
%title('Same Label Confidence', 'FontSize', bigFS)
set(gca, 'FontSize', medFS)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function err = objective(param, econf, elab)
econf = 1./(1+exp(param(1)*econf+param(2)));
px = [0.025:0.05:0.975];
f1 = ksdensity(econf(elab==0), px, 'support', [0 1])+1E-10;
f2 = ksdensity(econf(elab==1), px, 'support', [0 1])+1E-10;
f1 = f1 / sum(f1+f2);
f2 = f2 / sum(f1+f2);
err = sum((f1+f2).*(px - f2./(f1+f2)).^2);
disp(num2str([sum((f1+f2).*abs(px - f2./(f1+f2))) param]))

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

function [data, lab] = formatData(features, adjlist, labels)
% concatenate data and select ndata random datapoints

nimages = numel(features);

[tmp, nvars] = size(features{1});

% count edges
ne = 0;
for f = 1:nimages
    ne = ne + size(features{f}, 1);
end

data = zeros(ne, nvars);
lab = zeros(ne, 1);

% concatenate data
c = 0;
for f = 1:nimages
    cf = size(features{f}, 1);
    data(c+1:c+cf, :) = features{f};
    s1 = adjlist{f}(:, 1);
    s2 = adjlist{f}(:, 2);
    lab(c+1:c+cf) = (labels{f}(s1)==labels{f}(s2));
    c = c + cf;
end

author: zouchuhang

repo: LayoutNet


text:
% test for patchwise denoising
%
%   Copyright (c) 2007 Gabriel Peyre

path(path,'toolbox/');

%% load the image
rep = 'images/';    % directory where you can find the image
name = 'barb';
n = 80;
M = load_image([rep name]);
M = rescale( crop(M,n) );

sigma = 0.03; % variance of additional noise
if sigma>0
    % avoid saturation
    M = clamp( rescale(M,sigma,1-sigma) + sigma * randn(size(M)) );
end

%% options of NL means
options.k = 3;          % half size for the windows
options.T = 0.03;       % width of the gaussian, relative to max(M(:))  (=1 here)
options.max_dist = 15;  % search width, the smaller the faster the algorithm will be
options.ndims = 30;     % number of dimension used for distance computation (PCA dim.reduc. to speed up)
options.do_patchwise = 1;

%% do denoising
tic;
[M1,Wx,Wy] = perform_nl_means(M, options);
toc;

%% display results
ax = [];
clf;
ax(1) = subplot(2,2,1);
imagesc(M); axis image; axis off;
title('Original image');
ax(2) = subplot(2,2,2);
imagesc(clamp(M1)); axis image; axis off;
title('Denoised');
ax(3) = subplot(2,2,3);
imagesc(rescale(M-M1));
title('Removed noise');
axis image; axis off;
if size(M,3)>1
    colormap gray(256);
end
linkaxes(ax,'xy');

author: gpeyre

repo: matlab-toolboxes


text:
function img_list = batch_imread(batch)
% -------------------------------------------------------------------------
%   Description:
%       read a batch of images
%
%   Input:
%       - batch : array of ID to fetch
%
%   Output:
%       - img_list: batch of images
%
%   Citation:
%       Deep Laplacian Pyramid Networks for Fast and Accurate Super-Resolution
%       Wei-Sheng Lai, Jia-Bin Huang, Narendra Ahuja, and Ming-Hsuan Yang
%       IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017
%
%   Contact:
%       Wei-Sheng Lai
%       wlai24@ucmerced.edu
%       University of California, Merced
% -------------------------------------------------------------------------

    img_list = cell(length(batch), 1);

    for i = 1:length(batch)
         img = imread(batch{i});

         if( size(img, 3) > 1 )
             img = rgb2ycbcr(img);
             img = img(:, :, 1);
         end

         img_list{i} = im2single(img);
    end

end

repo: LapSRN


text:
function [ y, m, d, f ] = jed_to_ymdf_julian3 ( jed )

%*****************************************************************************80
%
%% JED_TO_YMDF_JULIAN3 converts a JED to a Julian YMDF date.
%
%  Discussion:
%
%    The theory behind this routine is very clean.  The Julian
%    calendar has cycles of 1 and 4 years, and we can analyze a date
%    by determining where it lies within these cycles.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Reingold, Nachum Dershowitz, Stewart Clamen,
%    Calendrical Calculations, II: Three Historical Calendars,
%    Software - Practice and Experience,
%    Volume 23, Number 4, pages 383-404, April 1993.
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, integer Y, M, D, real F,
%    the YMDF date.
%
  cycle_length = 1461;
  year_length = 365;

  jed_epoch = epoch_to_jed_julian ( );

  j = floor ( jed - jed_epoch );
  f1 = ( jed - jed_epoch ) - j;

  if ( f1 < 0.0 )
    f1 = f1 + 1.0;
    j = j - 1;
  end

  d0 = j;
  n4 = 0;

  while ( d0 <= 0 )
    d0 = d0 + cycle_length;
    n4 = n4 - 1;
  end

  n4 = n4 + floor ( d0 / cycle_length );
  d1 = i4_modp ( d0, cycle_length );
  n1 = floor ( d1 / year_length );
  d2 = i4_modp ( d1, year_length );

  if ( n1 == 4 )
    j1 = 366;
    y1 = 4 * n4 + n1;
  else
    j1 = d2 + 1;
    y1 = 4 * n4 + n1 + 1;
  end
%
%  Any year before 1 AD must be moved one year further back, since
%  this calendar does not include a year 0.
%
  y1 = y_astronomical_to_common ( y1 );

  [ y, m, d, f ] = yjf_to_ymdf_julian ( y1, j1, f1 );

  return
end


text:
function mbasis = basis_matrix_beta_uni ( beta1, beta2 )

%*****************************************************************************80
%
%% BASIS_MATRIX_BETA_UNI sets up the uniform beta spline basis matrix.
%
%  Discussion:
%
%    If BETA1 = 1 and BETA2 = 0, then the beta spline reduces to
%    the B spline.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Foley, van Dam, Feiner, Hughes,
%    Computer Graphics: Principles and Practice,
%    page 505.
%
%  Parameters:
%
%    Input, real BETA1, the skew or bias parameter.
%    BETA1 = 1 for no skew or bias.
%
%    Input, real BETA2, the tension parameter.
%    BETA2 = 0 for no tension.
%
%    Output, real MBASIS(4,4), the basis matrix.
%
  mbasis = zeros(4,4);

  mbasis(1,1) = - 2.0 * beta1 * beta1 * beta1;
  mbasis(1,2) =   2.0 * beta2 ...
    + 2.0 * beta1 * ( beta1 * beta1 + beta1 + 1.0 );
  mbasis(1,3) = - 2.0 * ( beta2 + beta1 * beta1 + beta1 + 1.0 );
  mbasis(1,4) =   2.0;

  mbasis(2,1) =   6.0 * beta1 * beta1 * beta1;
  mbasis(2,2) = - 3.0 * beta2 ...
    - 6.0 * beta1 * beta1 * ( beta1 + 1.0 );
  mbasis(2,3) =   3.0 * beta2 + 6.0 * beta1 * beta1;
  mbasis(2,4) =   0.0;

  mbasis(3,1) = - 6.0 * beta1 * beta1 * beta1;
  mbasis(3,2) =   6.0 * beta1 * ( beta1 - 1.0 ) * ( beta1 + 1.0 );
  mbasis(3,3) =   6.0 * beta1;
  mbasis(3,4) =   0.0;

  mbasis(4,1) =   2.0 * beta1 * beta1 * beta1;
  mbasis(4,2) =   4.0 * beta1 * ( beta1 + 1.0 ) + beta2;
  mbasis(4,3) =   2.0;
  mbasis(4,4) =   0.0;

  delta = ( ( 2.0   ...
    * beta1 + 4.0 ) ...
    * beta1 + 4.0 ) ...
    * beta1 + 2.0 + beta2;

  mbasis(1:4,1:4) = mbasis(1:4,1:4) / delta;

  return
end


text:
classdef mme_reserve_gen < mp.mm_element

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

%     properties
%     end

    methods
        function name = name(obj)
            name = 'reserve_gen';
        end

        function obj = add_vars(obj, mm, nm, dm, mpopt)
            dme = obj.data_model_element(dm);
            mm.add_var('R', dme.n, 0, 0, dme.r_ub);
        end

        function obj = add_costs(obj, mm, nm, dm, mpopt)
            dme = obj.data_model_element(dm);
            c = dme.tab.cost(dme.on) * dm.base_mva;    %% p.u. cost coeffs
            mm.add_quad_cost('Rcost', [], c, 0, {'R'});
        end

        function obj = add_constraints(obj, mm, nm, dm, mpopt)
            dme = obj.data_model_element(dm);
            gen_dme = dm.elements.gen;
            ng = gen_dme.n; %% number of online generators
            ngr = dme.n;    %% number of online reserve gens
            igr = dme.gen;  %% index of online gen for each online reserve gen

            Ar = [sparse(1:ngr, igr, 1, ngr, ng) speye(ngr)];
            ur = gen_dme.pg_ub(igr);

            mm.add_lin_constraint('Pg_plus_R', Ar, [], ur, {'Pg', 'R'});
        end

        function obj = data_model_update(obj, mm, nm, dm, mpopt)
            dme = obj.data_model_element(dm);
            rz_dme = dm.elements.reserve_zone;

            %% get reserve allocations, bounds, shadow prices, cost
            [vv, ll] = mm.get_idx('var', 'lin');
            r = mm.soln.x(vv.i1.R:vv.iN.R) * dm.base_mva;
            r_lb = 0;
            r_ub = dme.r_ub * dm.base_mva;
            zprc = zeros(rz_dme.nr, 1);
            zprc(rz_dme.on) = mm.soln.lambda.mu_l(ll.i1.Rreq:ll.iN.Rreq) / dm.base_mva;
            prc = rz_dme.tab.zones(:, dme.tab.gen)' * zprc;
            mu_lb = mm.soln.lambda.lower(vv.i1.R:vv.iN.R) / dm.base_mva;
            mu_ub = mm.soln.lambda.upper(vv.i1.R:vv.iN.R) / dm.base_mva;
            mu_pg_ub = mm.soln.lambda.mu_u(ll.i1.Pg_plus_R:ll.iN.Pg_plus_R) / dm.base_mva;
            total_cost = mm.eval_quad_cost(mm.soln.x, 'Rcost');

            %% update in the data model
            dme.tab.r(dme.on) = r;
            dme.tab.r_lb(dme.on) = r_lb;
            dme.tab.r_ub(dme.on) = r_ub;
            dme.tab.prc = prc;
            dme.tab.mu_lb(dme.on) = mu_lb;
            dme.tab.mu_ub(dme.on) = mu_ub;
            dme.tab.mu_pg_ub(dme.on) = mu_pg_ub;
            dme.tab.total_cost(dme.on) = total_cost;
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function r = not(p)
%UMINUS: this function applies a unary minus to the x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

r=p;
r.xvalue = ~p.xvalue;
r.yvalue = ~p.yvalue;
r.zvalue = ~p.zvalue;

author: simulkade

repo: FVTool


text:
function x=pcmu2lin(p,s)
%PCMU2LIN Convert Mu-law PCM to linear X=(P,S)
%	lin = pcmu2lin(pcmu) where pcmu contains a vector
%	of mu-law values in the range 0 to 255.
%	No checking is performed to see that numbers are in this range.
%
%	Output values are divided by the scale factor s:
%
%		   s		Output Range
%
%		   1		+-8031	(integer values)
%		4004.2	+-2.005649 (default)
%		8031		+-1
%		8159		+-0.9843118 (+-1 nominal full scale)
%
%	The default scaling factor 4004.189931 is equal to
%	sqrt((2207^2 + 5215^2)/2) this follows ITU standard G.711.
%	The sine wave with PCM-Mu values [158 139 139 158 30 11 11 30]
%	has a mean square value of unity corresponding to 0 dBm0.

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: pcmu2lin.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<2
  t=9.98953613E-4;
else
  t=4/s;
end

m=15-rem(p,16);
q=floor(p/128);
e=(127-p-m+128*q)/16;
x=(q-0.5).*(pow2(m+16.5,e)-16.5)*t;

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD