text:
%mmf writer
function [keyframes,nframes]=mmfwritertrans(inputfile,outputfile,obj,firstframe,lastframe)

%obj=VideoReader(inputfile);
foutput=fopen(outputfile,'wb');
%load text
ftexto=fopen('/groups/visitors/home/albam/Desktop/mmfreader/text.text');

texto=[];
for i=1:15
linea=fgets(ftexto);
texto=[texto linea];
end

%write header text
fwrite(foutput,texto,'char')
% end of the header
fwrite(foutput,0,'char')
%id_code
d=hex2dec('a3d2d45d');
fwrite(foutput,d,'ulong')
%header size
fwrite(foutput,10240,'int')
%calculate keyframes
keyframes=fix((lastframe-firstframe+1)/180)+1;
%keyframe
fwrite(foutput,keyframes,'int')
%other recording parameters TODO
% zero padded to reach header size
for i=ftell(foutput)+1:10240
fwrite(foutput,0,'uint8');
end
%calculate nframes
%
nframes=zeros(1,keyframes);
nframes(1:end-1)=180;
nframes(end)=(lastframe-firstframe+1)-(keyframes-1)*180;

isfirst=1;
% for the stacks
for stack=1:keyframes

%Image Stacks
begining_header=ftell(foutput);
%idcode
d=hex2dec('bb67ca20');
fwrite(foutput,d,'ulong')
%header size
fwrite(foutput,512,'int')
%total size of stack in disk TODO
fwrite(foutput,0,'int')
%nframes in this stack default=180 but the last one can be smaller
fwrite(foutput,nframes(stack),'int');
% zero padded to reach header size
for i=ftell(foutput)-begining_header+1:512
fwrite(foutput,0,'uint8');
end

%Backgorund image
begining_header=ftell(foutput);
%header size
fwrite(foutput,112,'int')
%imageid
fwrite(foutput,0,'int')
%numchannels
fwrite(foutput,1,'int')
%alphachannel
fwrite(foutput,0,'int')
%depth
fwrite(foutput,8,'int')
%colormodel
fwrite(foutput,[0 0 0 0],'char')
%channel seq
fwrite(foutput,[0 0 0 0],'char')
%data order
fwrite(foutput,0,'int')
% origin
fwrite(foutput,0,'int')
%align should 4 or 8 TODO
fwrite(foutput,0,'int')
%width
fwrite(foutput,obj.width,'int')
%height
fwrite(foutput,obj.height,'int')
%roi 0=> full image
fwrite(foutput,0,'int')
%mask roi
fwrite(foutput,0,'int')
%im id
fwrite(foutput,0,'int')
%tile info
fwrite(foutput,0,'int')
%images size (heightxwidth
fwrite(foutput,obj.width*obj.height,'int')
%imagedata
fwrite(foutput,0,'int')
%widthstep
fwrite(foutput,obj.width,'int')
%bordermode
fwrite(foutput,0,'int')
%borderconst
fwrite(foutput,0,'int')
%imagedataorigin
fwrite(foutput,0,'int')
% zero padded to reach header size
for i=ftell(foutput)-begining_header+1:112
fwrite(foutput,0,'uint8');
end
s=0;
im=zeros(obj.width,obj.height,nframes(stack));
if isfirst==1
    for i=firstframe:firstframe+nframes(stack);
        s=s+1;
        imrgbtemp=read(obj,i);
        imrgbtemp=double(imrgbtemp);
        imtemp=sum(imrgbtemp,3)/3;
        im(:,:,s)=imtemp';
    end
    isfirst=0;
else
    for i=firstframe+sum(nframes(1:stack-1)):firstframe+sum(nframes(1:stack))
        s=s+1;
        imrgbtemp=read(obj,i);
        imrgbtemp=double(imrgbtemp);
        imtemp=sum(imrgbtemp,3)/3;
        im(:,:,s)=imtemp';
    end
end

%%creting mean backgfix
bkgim=median(im,3);
%writing background image
fwrite(foutput,reshape(bkgim',1,obj.width*obj.height))

%Extracting larvae out of the images

 frames=zeros(size(im));
 SE=strel('disk',8);
 bounding=cell(1,nframes(stack));

for i=1:nframes(stack)
   begining_header=ftell(foutput);
   %id code
   d=hex2dec('f80921af');
   fwrite(foutput,d,'ulong')
   %header size
   fwrite(foutput,1024,'int')
   %depth
   fwrite(foutput,8,'int')
   %nchannels
   fwrite(foutput,1,'int')

   %calculating numims
    frames(:,:,i)=abs(im(:,:,i)-bkgim);
    %threshold automatically set
    %level(i)=graythresh(frames(:,:,i));
     bw=zeros(size(frames(:,:,1)));
     bw(frames(:,:,i)*10>max(max(frames(:,:,i))))=1;
     bw2=bwareaopen(bw,7);
     dilated=imdilate(bw2,SE);
     cc=bwconncomp(dilated);

     %numims
     fwrite(foutput,cc.NumObjects,'int')
     % zero padded to reach header size
        for k=ftell(foutput)-begining_header+1:1024
        fwrite(foutput,0,'uint8');
        end
%     %Creating bounding box
     STATS=regionprops(cc,'BoundingBox');
      bounding{i}=zeros(cc.NumObjects,4);
     for j=1:cc.NumObjects
         bounding{i}(j,:)=STATS(j,1).BoundingBox;
         %temp=zeros(fix(bounding{i}(j,4)),fix(bounding{i}(j,3)));
         if fix(bounding{i}(j,2))==0 && fix(bounding{i}(j,1))==0
            temp=im(1:fix(bounding{i}(j,4)),1:fix(bounding{i}(j,3)),i);
        elseif fix(bounding{i}(j,2))==0
            temp=im(1:fix(bounding{i}(j,4)),fix(bounding{i}(j,1)):fix(bounding{i}(j,1))+fix(bounding{i}(j,3))-1,i);
        elseif fix(bounding{i}(j,1))==0
            temp=im(fix(bounding{i}(j,2)):fix(bounding{i}(j,2))+fix(bounding{i}(j,4))-1,1:fix(bounding{i}(j,3)),i);
        else
            temp=im(fix(bounding{i}(j,2)):fix(bounding{i}(j,2))+fix(bounding{i}(j,4))-1,fix(bounding{i}(j,1)):fix(bounding{i}(j,1))+fix(bounding{i}(j,3))-1,i);
         end
        fwrite(foutput,fix(bounding{i}(j,[2 1 4 3])),'int')
        fwrite(foutput,reshape(temp',1,fix(bounding{i}(j,3))*fix(bounding{i}(j,4))))
      end

end

end
fclose(foutput)

author: kristinbranson

repo: JAABA


text:
function plotStats(obj, epochs, stats, plotAccuracy) %#ok<INUSL>
% plotStats(obj, epochs, stats, plotAccuracy)
%
% Plot the results of different metrics on train and validation set.
%
% Copyright by Holger Caesar, 2016

% Define color map
cmap = parula(3);
cmap(3, :) = [255, 140, 0] / 255;

if true
    figure(1);
    clf;
    if plotAccuracy
        hsp1 = subplot(2, 1, 1);
    end
    hold on;
    leg = {};
    datasetModes = {'train', 'val'};
    for datasetModeIdx = 1 : numel(datasetModes)
        datasetMode = datasetModes{datasetModeIdx};
        field = 'objective';
        fieldValues = [stats.(datasetMode).(field)];
        if datasetModeIdx == 1
            marker = '-';
        else
            marker = '--';
        end

        % For i.e. objective with 1 value
        leg{end + 1} = sprintf('%s (%s)', field, datasetMode); %#ok<AGROW>
        values = fieldValues(1, :)';
        plot(epochs, values, 'Color', cmap(1, :), 'LineStyle', marker);
    end
    hleg1 = legend(leg, 'Location', 'NorthEastOutSide');

    xlabel('epoch');
    ylabel('objective');
    grid on;
end

if plotAccuracy
    hsp2 = subplot(2, 1, 2);
    hold on;

    leg = {};
    datasetModes = {'train', 'val'};
    for datasetModeIdx = 1 : numel(datasetModes)
        datasetMode = datasetModes{datasetModeIdx};
        field = 'accuracy';
        fieldValues = [stats.(datasetMode).(field)];
        if datasetModeIdx == 1
            marker = '-';
        else
            marker = '--';
        end

        leg{end + 1} = sprintf('Pix. Acc. (%s)', datasetMode); %#ok<AGROW>
        values = fieldValues(1, :)';
        plot(epochs, values, 'Color', cmap(1, :), 'LineStyle', marker);

        leg{end + 1} = sprintf('Class. Acc. (%s)', datasetMode); %#ok<AGROW>
        values = fieldValues(2, :)';
        plot(epochs, values, 'Color', cmap(2, :), 'LineStyle', marker);

        leg{end + 1} = sprintf('Mean IU (%s)', datasetMode); %#ok<AGROW>
        values = fieldValues(3, :)';
        plot(epochs, values, 'Color', cmap(3, :), 'LineStyle', marker);
    end

    % Define legend
    hleg2 = legend(leg, 'Location', 'NorthEastOutSide');

    % Align both plots
    if false
        plotPos1 = get(hsp1, 'Position');
        plotPos2 = get(hsp2, 'Position');
        set(hsp2, 'Position', [plotPos1(1), plotPos2(2), plotPos1(3:4)]);
        legPos1 = get(hleg1, 'Position');
        legPos2 = get(hleg2, 'Position');
        set(hleg2, 'Position', [legPos1(1), legPos2(2), legPos1(3:4)]);
    end

    xlabel('epoch');
    ylabel('accuracy');
    grid on;
end

author: nightrome

repo: matconvnet-calvin


text:
function lebesgue_test ( )

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

  lebesgue_test01 ( );
  lebesgue_test02 ( );
  lebesgue_test03 ( );
  lebesgue_test04 ( );
  lebesgue_test05 ( );
  lebesgue_test06 ( );
  lebesgue_test07 ( );
  lebesgue_test08 ( );
  lebesgue_test09 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LEBESGUE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [scales,weights,covar]=realized_quantile_variance_scale(samplesperbin,quantiles,simulations,symmetric)
% Computes the scales needed for estimating the integrated variance using Realized Quantile
% Variance.  Also computes the weights of the optimal combination and non-scaled covariance from
% which the weights are derived.
%
% USAGE:
%   [SCALES,WEIGHTS,COVAR]=realized_quantile_variance_scale(SAMPLESPERBIN,QUANTILES,SIMULATIONS)
%
% INPUTS:
%   SAMPLESPERBIN    - Number of returns to use in each bin when computing the quantiles. NOTE: The
%                        number of returns produced by filtering according to SAMPLINGTYPE and
%                        SAMPLINGINTERVAL must be an integer multiple of SAMPLESPERBIN.
%   QUANTILES        - k by 1 vector of quantile values to use when computing RQ, must satisfy 0.5<QUANTILES<=1.
%                        Quantiles must produce values such that QUANTILES*SAMPLESPERBIN is an
%                        integer.  The simplest method to accomplish this is to specify QUANTILES as
%                        the ratio of the index number of the return to SAMPLESPERBIN (e.g.
%                        QUANTILES = [13 15 19]/20)
%   SIMULATIONS       - [OPTIONAL] Scalar integer.  Should usually be at least 1,000,000.  Default
%                         value is 10,000,000.
%   SYMMETRIC         - [OPTIONAL] Logical value indicating whether the symmetric estimator is being
%                         used. Default is false.
%
% OUTPUTS:
%   SCALES            - k by 1 vector of scales for standardized quantile based estimates of variance
%   WEIGHTS           - k by 1 vector of weights that produce the optimal combination (lowest variance)
%   COVAR             - k by k non-scaled covariance matrix of the quantile realized variance - one
%                         for each quantile used.  Forms the actual (scaled) covariance when
%                         multiplied by the integrated quarticity.
%
% COMMENTS:
%   Uses Monte Carlo integration with 10,000,000 simulations.  Trade offs between accuracy and time
%   can be made by altering SIMULATIONS.
%
%  See also REALIZED_QUANTILE_VARIANCE
%

%
% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 5/1/2008

if nargin<4
    symmetric = false;
end
if nargin ==2
    simulations = 10000000;
end

% Save the state and set it to a specific value so that the simulated results will be reproducable
state0=randn('state');
randn('state',datenum('MAR-26-1974'))

k = length(quantiles);

if symmetric
    indices = round(samplesperbin*quantiles);
else
    indicesHigh = round( samplesperbin*quantiles);
    indicesLow = round(samplesperbin*(1-quantiles)+1 );
end

% Need to block eventually incase someone tries to compute this with a
% large number of samplesperbin
rows = min(floor(2^13/samplesperbin),simulations);
iter = ceil(simulations/rows);
remaining = simulations;

scales = zeros(1,k);
opMatrix  = zeros(k);
% The MC integration
for i=1:iter
    numThisIter = min(remaining,rows);
    remaining = remaining - rows;

    x=randn(numThisIter,samplesperbin);
    if symmetric
        x=sort(x.^2,2);
        x2=x(:,indices);
    else
        x=sort(x,2);
        x2=x(:,indicesHigh).^2 + x(:,indicesLow).^2;
    end
    scales = scales + sum(x2)/simulations;
    opMatrix = opMatrix + x2'*x2/simulations;
end

% Compute the covariance
covar = opMatrix - scales'*scales;
covar = samplesperbin*diag(1./scales)*covar*diag(1./scales);
% And the weights according to the GMVP
weights=covar^(-1)*ones(k,1)/(ones(1,k)*covar^(-1)*ones(k,1));

% Reset the state
randn('state',state0);

author: bashtage

repo: mfe-toolbox


text:
function [wout, hout] = nomalize_wh(win,hin)
% [wout, hout] = nomalize_wh(win,hin,k)

sumw = sum(win,1);
wout = win./repmat(sumw,size(win,1),1);                %normalization of each component
hout = hin.*repmat(sumw',1,size(hin,2));    %to keep the multiplication equal

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef Mesh_Composite < AbstractMesh

    properties (GetAccess = public, SetAccess = public)
        meshes
        nMeshes
        nActiveMeshes
        activeMeshesList

        globalConnectivities
    end

    methods (Access = public)

        function obj = Mesh_Composite()
            obj.createMeshes();
        end

        function aMeshes = getActiveMeshes(obj)
            aMeshes = cell(1,obj.nActiveMeshes);
            for iActive = 1:obj.nActiveMeshes
                iMesh = obj.activeMeshesList(iActive);
                aMeshes{iActive} = obj.meshes{iMesh};
            end
        end
    end

    methods (Access = protected)

        function append(obj,mesh)
           obj.meshes{end+1} = mesh;
        end

    end

    methods (Access = private)

        function createMeshes(obj)

        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% function qs = equantile(X, ps, dim)
%
% Calculates empirical (non-smoothed) quantiles of multidimensional
% arrays.  Ignores NaN entries.  Faster than, but similar to, the
% quantile() function in the MATLAB Statistics toolbox.
%
% Input:
%  X: multidimensional vector/matrix over which to calculate
%     empirical quantiles.
%  ps: vector of quantile values (between 0 and 1)
%  dim: (optional) dimension along which to calculate quantiles
%
% Output:
%  qs: vector/matrix of empirical quantiles
%
%---------------------------------------------------------------------------------
%    Author: Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/)
%---------------------------------------------------------------------------------
function qs = equantile(X,ps,dim)

ndX = ndims(X);
sX = size(X);
psn = length(ps);

assert(all(ps >= 0 & ps <= 1));

needsort = 0;
if nargin<3
  % Sort properly depending on data
  if isvector(X) && ~issorted(X)
    needsort = 1;
  elseif ndX <= 2 && ~issorted(X,'rows')
    needsort = 1;
  elseif ndX > 2
    needsort = 1;
  end

  % Which dimension are we calculating quantiles along?
  if isvector(X)
    [tmpvar,mdim] = max(sX);
    dim = mdim;
  else
    dim = 1;
  end
else
  needsort = 1;
end

% Reshape data to make ensuing steps easier (quantile-dim becomes first dim)
dperm = 1:ndX;
dperm([1 dim]) = dperm([dim 1]);
if (dim ~= 1)
  X = permute(X, dperm);
end
sXperm = sX(dperm);

% Reshape to 2D; we calculate quantiles along first dimension
X = reshape(X, sX(dim), prod(sX([1:dim-1, dim+1:ndX])));

if (needsort),
  X = sort(X,1);
end

% If we see NaNs, need to calculate quantiles row-by-row.
% Otherwise can do it in one fell swoop.
Xnan = isnan(X);
if any(any(Xnan,2))
  hasnan = 1;
else
  hasnan = 0;
end

if (~hasnan)
  ks = round(ps*(sX(dim)-1) + 1);
  assert(all(ks >= 1 & ks <= sX(dim)));
  qs = X(ks,:);
else % if (hasnan)
  ns = zeros(size(X,2), 1);
  qs = nan(psn, size(X,2));

  % Limit ourselves to non-NaN values
  for ci=1:size(X,2)
    ni = find(Xnan(:, ci), 1, 'first')-1;
    % If don't find any NaNs, use full dim
    if isempty(ni), ni = sX(dim); end

    ns(ci) = ni;
  end

  % Calculate quantiles in each case
  Xcol = 1:size(X,2);
  for pi=1:psn
    p = ps(pi);
    ks = round(p*(ns-1)+1);
    % Pull out the proper elements of sorted X
    qs(pi,ks>0) = X(sub2ind(size(X), ks(ks>0)', Xcol(ks>0)));
  end
end

% Convert back to ndim array
qs = reshape(qs, [psn, sX(1:dim-1), sX(dim+1:ndX)]);

% Move coordinates back to their proper places
if (dim ~= 1)
  qs = ipermute(qs, dperm);
end

author: HeLiangHIT


text:
function [Population,Fitness] = Second_Stage_EnvironmentalSelection(Population,N,isOrigin)
% The environmental selection of SPEA2

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

% This function is written by Kangjia Qiao

    %% Calculate the fitness of each solution
    if isOrigin==1
    %     con = Population.cons;
    %     con(find(con<=0)) = 0;
        Fitness = CalFitness(Population.objs,Population.cons);
    else
        Fitness = CalFitness(Population.objs);
    end

    %% Environmental selection
    Next = Fitness < 1;
    if sum(Next) < N
        [~,Rank] = sort(Fitness);
        Next(Rank(1:N)) = true;
    elseif sum(Next) > N
        Del  = Truncation(Population(Next).objs,sum(Next)-N);
        Temp = find(Next);
        Next(Temp(Del)) = false;
    end
    % Population for next generation
    Population = Population(Next);
    Fitness    = Fitness(Next);
    % Sort the population
    [Fitness,rank] = sort(Fitness);
    Population = Population(rank);
end

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

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

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [f,fu,fx] = f_fu_fx_Wrapper(u,x,p,parIdx)

% dx = f(u,x,p)
    coder.cinclude('iiwa14.h');
    [xDim,~] = size(x);
    [uDim,~] = size(u);
    f = f_Wrapper(u,x,p,parIdx);
    fu = zeros(xDim,uDim);
    fx = zeros(xDim,xDim);

    q = x(1:7,1);
    qd = x(8:end,1);
    tau = u(1:7,1);

    dq  = zeros(7,7);
    dqd = zeros(7,7);
    dtau= zeros(7,7);

    coder.ceval('derivatives_cal', ...
                coder.ref(q),...
                coder.ref(qd),...
                coder.ref(tau),...
                coder.ref(dq),...
                coder.ref(dqd),...
                coder.ref(dtau),...
                parIdx);
    fu(8:end,1:7) = dtau;
    fx = [zeros(7,7),eye(7);dq,dqd];
end

author: deng-haoyang

repo: ParNMPC


text:
function[R]=radearth
%RADEARTH  The radius of the earth in kilometers.
%
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2006 J.M. Lilly --- type 'help jlab_license' for details

R=6371;

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [] = test_trace_norm_mc()

    clc;
    clear;
    close all;

    rng('default');

    %% Set algorithms
    if 0
        algorithms = gd_solver_list('ALL');
    else
        %algorithms = {'PG-BKT', 'PG-TFOCS-BKT', 'APG-BKT', 'APG-TFOCS-BKT', 'L-BFGS-BKT'};
        algorithms = {'PG-TFOCS-BKT', 'APG-TFOCS-BKT', 'L-BFGS-BKT'};
    end

    %% prepare dataset
    if 1
        % generate synthtic data
        n = 100;
        m = 50;
        r = 10;
        density = 0.2;
        lambda = 5;
        M = randn(m,r)*randn(r,n);
        mask = (rand(m,n)<density);
    else
    end

    %% define problem definitions
    problem = trace_norm_matrix_completion(M, mask, lambda);

    %% initialize
    w_init = randn(m*n, 1);
    w_list = cell(length(algorithms),1);
    info_list = cell(length(algorithms),1);

    %% perform algorithms
    for alg_idx=1:length(algorithms)
        fprintf('\n\n### [%02d] %s ###\n\n', alg_idx, algorithms{alg_idx});

        clear options;
        % general options for optimization algorithms
        options.w_init = w_init;
        options.tol_gnorm = 1e-10;
        options.max_iter = 100;
        options.verbose = true;
        options.f_opt = 0;
        options.store_w = false;

        switch algorithms{alg_idx}
            case {'PG-BKT'}

                options.step_alg = 'backtracking';
                options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'PG-TFOCS-BKT'}

                options.step_alg = 'tfocs_backtracking';
                options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'APG-BKT'}

                options.step_alg = 'backtracking';
                options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = ag(problem, options);

            case {'APG-TFOCS-BKT'}

                options.step_alg = 'tfocs_backtracking';
                options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = ag(problem, options);

            case {'L-BFGS-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'L-BFGS-WOLFE'}

                options.step_alg = 'strong_wolfe';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            otherwise
                warn_str = [algorithms{alg_idx}, ' is not supported.'];
                warning(warn_str);
                w_list{alg_idx} = '';
                info_list{alg_idx} = '';
        end

    end

    fprintf('\n\n');

    %% plot all
    close all;

    % display iter vs. cost
    display_graph('iter','cost', algorithms, w_list, info_list);
    % display iter vs. trace (nuclear) norm
    display_graph('iter','reg', algorithms, w_list, info_list);
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function output = callipopt(model)

options = [];
try
    options.ipopt = optiRemoveDefaults(model.options.ipopt,ipoptset());
catch
    options.ipopt = model.options.ipopt;
end
options.ipopt.print_level = 2+model.options.verbose;

model = yalmip2nonlinearsolver(model);

if ~model.derivative_available
    disp('Derivate-free call to ipopt not yet implemented')
    error('Derivate-free call to ipopt not yet implemented')
end
if model.options.savedebug
    save ipoptdebug model
end
showprogress('Calling IPOPT',model.options.showprogress);

% Figure out which variables are artificially introduced to normalize
% arguments in callback operators to simplify chain rules etc. We can do
% our function evaluations and gradient computations in our lifted world,
% but only expose the model in the original variables to the nonlinear
% solver.
% model = compressLifted(model);

Fupp = [ repmat(0,length(model.bnonlinineq)+length(model.K.q)*(model.K.q(1)>0),1);
    repmat(0,length(model.bnonlineq),1);
    repmat(0,length(model.b),1);
    repmat(0,length(model.beq),1)];

Flow = [ repmat(-inf,length(model.bnonlinineq)+length(model.K.q)*(model.K.q(1)>0),1);
    repmat(0,length(model.bnonlineq),1);
    repmat(-inf,length(model.b),1);
    repmat(0,length(model.beq),1)];

if isempty(Flow)
    Flow = [];
    Fupp = [];
end

% Since ipopt react strangely on lb>ub, we should bail if that is detected
% (ipopt creates an exception)
if ~isempty(model.lb)
    if any(model.lb>model.ub)
        problem = 1;
        solverinput = [];
        solveroutput = [];
        output = createoutput(model.c*0,[],[],problem,'IPOPT',solverinput,solveroutput,0);
        return
    end
end

% These are needed to avoid recomputation due to ipopts double call to get
% f and df, and g and dg
global latest_x_f
global latest_x_g
global latest_df
global latest_f
global latest_G
global latest_g
global latest_xevaled
global latest_x_xevaled
latest_G = [];
latest_g = [];
latest_x_f = [];
latest_x_g = [];
latest_xevaled = [];
latest_x_xevaled = [];

funcs.objective = @(x)ipopt_callback_f(x,model);
funcs.gradient = @(x)ipopt_callback_df(x,model);
if ~isempty(Fupp)
    funcs.constraints = @(x)ipopt_callback_g(x,model);
    funcs.jacobian  = @(x)ipopt_callback_dg(x,model);
end

options.lb = model.lb(:)';
options.ub = model.ub(:)';
if ~isempty(Fupp)
    options.cl = Flow;
    options.cu = Fupp;
end

if ~isempty(Fupp)
    Z = jacobiansparsityfromnonlinear(model);
    funcs.jacobianstructure = @() Z;
end

if ~model.options.usex0
    model.x0 = (options.lb+options.ub)/2;
    model.x0(isinf(options.ub)) = options.lb(isinf(options.ub))+1;
    model.x0(isinf(options.lb)) = options.ub(isinf(options.lb))-1;
    model.x0(isinf(model.x0)) = 0;
    if any(model.variabletype == 4)
        problematic = find(any(model.monomtable(:,model.linearindicies) < 0 ,1));
        if ~isempty(problematic)
            problematic = problematic(find(model.x0(problematic)==0));
            Oneisfeas = problematic(find(model.ub(problematic) > 1));
            model.x0(Oneisfeas) = 1;
        end
    end
    model.x0(find(model.lb==model.ub)) = model.lb(find(model.lb==model.ub));
end

% If quadratic objective and no nonlinear constraints, we can supply an
% Hessian of the Lagrangian
usedinObjective = find(model.c | any(model.Q,2));
if ~any(model.variabletype(usedinObjective)) & any(model.Q)
    if  length(model.bnonlinineq)==0 & length(model.bnonlineq)==0
        H = model.Q(:,model.linearindicies);
        H = H(model.linearindicies,:);
        funcs.hessian = @(x,s,l) tril(2*H);
        funcs.hessianstructure = @()tril(sparse(double(H | H)));
        options.ipopt.hessian_approximation = 'exact';
    end
end

solvertime = tic;
[xout,info] = ipopt(model.x0,funcs,options);
solvertime = toc(solvertime);

% Duals currently not supported
lambda = [];

if ~isempty(xout) && ~isempty(model.lift);
    x = zeros(length(model.linearindicies),1);
    x(model.lift.linearIndex) = xout(:);
    x(model.lift.liftedIndex) = model.lift.T*xout(:) + model.lift.d;
    x = RecoverNonlinearSolverSolution(model,x);
else
    x = RecoverNonlinearSolverSolution(model,xout);
end

switch info.status
    case {0,1}
        problem = 0;
    case {2}
        problem = 1;
    case {-1}
        problem = 3;
    case {3,4,-2,-3}
        problem = 4;
    case {-11,-12,-13}
        problem = 7;
    case {-10,-100,-101,-102,-199}
        problem = 11;
    otherwise
        problem = -1;
end

% Internal format for duals
D_struc = [];

% Save all data sent to solver?
if model.options.savesolverinput
    solverinput.model = model;
else
    solverinput = [];
end

% Save all data from the solver?
if model.options.savesolveroutput
    solveroutput.x = xout;
    solveroutput.info = info;
else
    solveroutput = [];
end

% Standard interface
output = createoutput(x,D_struc,[],problem,'IPOPT',solverinput,solveroutput,solvertime);

author: shamilmamedov


text:
function view=xformFlatLevelMenu(view)
%
% view = xformFlatLevelMenu(view)
%
% Set up the callbacks for the xformView menu in the Flat Level view
%
% djh, 1/9/98
% rmk, 1/10/99 added xformParMap
% rmk, 1/15/99 added xformAllROIs
% ras, 10/04 created xformFlatLevelMenu

mrGlobals;

xformMenu = uimenu('Label','Xform','separator','on');

% Xform ROI callback:
%   gray=checkSelectedGray;
%   view=vol2flatCurROILevels(gray,view);
%   view=refreshScreen(view,0);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatCurROILevels(gray,',view.name,'); ',...
        view.name,'=refreshScreen(',view.name,',0);'];
uimenu(xformMenu,'Label','ROI (selected)','Separator','off',...
    'CallBack',callBackstr);

% Xform AllROIs callback:
%   gray=checkSelectedGray;
%   view=vol2flatAllROIsLevels(gray,view);
%   view=refreshScreen(view,0);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatAllROIsLevels(gray,',view.name,'); ',...
        view.name,'=refreshScreen(',view.name,',0);'];
uimenu(xformMenu,'Label','ROIs (all)','Separator','on',...
    'CallBack',callBackstr);

% Xform CorAnal all scans callback:
%   gray=checkSelectedGray;
%   view=vol2flatCorAnalLevels(gray,view,0);
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatCorAnalLevels(gray,',view.name,',0); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','CorAnal (all scans)','Separator','on',...
    'CallBack',callBackstr);

% Xform CorAnal current scan callback:
%   gray=checkSelectedGray;
%   view=vol2flatCorAnalLevels(gray,view,getCurScan(view));
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatCorAnalLevels(gray,',view.name,',getCurScan(',view.name,')); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','CorAnal (current scan)','Separator','on',...
    'CallBack',callBackstr);

% Xform CorAnal select scans callback:
%   gray=checkSelectedGray;
%   view=vol2flatCorAnalLevels(gray,view);
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatCorAnalLevels(gray,',view.name,'); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','CorAnal (select scans)','Separator','on',...
    'CallBack',callBackstr);

% Xform ParMap all scans callback:
%   gray=checkSelectedGray;
%   view=vol2flatParMapLevels(gray,view,0);
%   view=setDisplayMode(view,'map');
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatParMapLevels(gray,',view.name,',0); ',...
        view.name,'=setDisplayMode(',view.name,',''map''); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','Parameter Map (all scans)','Separator','on',...
    'CallBack',callBackstr);

% Xform ParMap current scan callback:
%   gray=checkSelectedGray;
%   view=vol2flatParMapLevels(gray,view,getCurScan(view));
%   view=setDisplayMode(view,'map');
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatParMapLevels(gray,',view.name,',getCurScan(',view.name,')); ',...
        view.name,'=setDisplayMode(',view.name,',''map''); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','Parameter Map (current scan)','Separator','on',...
    'CallBack',callBackstr);

% Xform ParMap select scans callback:
%   gray=checkSelectedGray;
%   view=vol2flatParMapLevels(gray,view);
%   view=setDisplayMode(view,'map');
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatParMapLevels(gray,',view.name,'); ',...
        view.name,'=setDisplayMode(',view.name,',''map''); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','Parameter Map (select scans)','Separator','on',...
    'CallBack',callBackstr);

% Xform tSeries all scans callback:
%   gray=checkSelectedGray;
%   view=vol2flatTSeriesLevels(gray,view,0);
%   view=setDisplayMode(view,'map');
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatTSeriesLevels(gray,',view.name,',0); ',...
        view.name,'=setDisplayMode(',view.name,',''map''); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','tSeries (all scans)','Separator','on',...
    'CallBack',callBackstr);

% Xform tSeries current scan callback:
%   gray=checkSelectedGray;
%   view=vol2flatTSeriesLevels(gray,view,getCurScan(view));
%   view=setDisplayMode(view,'map');
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatTSeriesLevels(gray,',view.name,',getCurScan(',view.name,')); ',...
        view.name,'=setDisplayMode(',view.name,',''map''); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','tSeries (current scan)','Separator','on',...
    'CallBack',callBackstr);

% Xform tSeries select scans callback:
%   gray=checkSelectedGray;
%   view=vol2flatTSeriesLevels(gray,view);
%   view=setDisplayMode(view,'map');
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatTSeriesLevels(gray,',view.name,'); ',...
        view.name,'=setDisplayMode(',view.name,',''map''); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','tSeries (select scans)','Separator','on',...
    'CallBack',callBackstr);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function out = gaussian_kernel(x,d1,d2,bandWidth)

% taken from a reference on the web to generate a Gaussian kernal
ns = 1000; % resolution
xs1 = linspace(0,bandWidth(1,1),ns+1); % spatial
xs2 = linspace(0,bandWidth(1,2),ns+1); %range
kfun1 = exp(-(xs1.^2)/(2*bandWidth(1,1)^2));
kfun2 = exp(-(xs2.^2)/(2*bandWidth(1,2)^2));
w1 = kfun1(1,1:size(d1)).*(round(d1/bandWidth(1,1)*ns)+1);
w2=kfun2(1,1:size(d2)).*(round(d2/bandWidth(1,2)*ns)+1);
w=w1+w2;
w = w/sum(w); % normalise
out = sum( bsxfun(@times, x, w ), 2 );
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%DEMO_POLYGONSKELETON  One-line description here, please.
%
%   output = demo_polygonSkeleton(input)
%
%   Example
%   demo_polygonSkeleton
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
% Created: 2022-02-18,    using Matlab 9.10.0.1739362 (R2021a) Update 5
% Copyright 2022 INRAE.

% Note: this example require the "matGraphs" toolbox
% (https://github.com/mattools/matGraphs)
% start from a binary shape
img = imread('circles.png');
img = imfill(img, 'holes');
% figure; imshow(img); hold on;

% compute a smooth contour
cntList = bwboundaries(img);
poly = cntList{1}(2:end,:);
poly = smoothPolygon(poly, 5);

% draw the polygon
figure; hold on; axis equal; axis([0 250 0 250]);
drawPolygon(poly, 'color', 'k', 'linewidth', 1);

% compute skeleton
[vertices, adjList] = polygonSkeleton(poly);
% convert adjacency list to an edge array
edges = adjacencyListToEdges(adjList);
% draw the skeleton graph
drawGraphEdges(vertices, edges, 'color', 'b');

print(gcf, 'demo_polygonSkeleton.png', '-dpng');

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [xopt,fopt,niter,gnorm,dx] = grad_descent(varargin)
% grad_descent.m demonstrates how the gradient descent method can be used
% to solve a simple unconstrained optimization problem. Taking large step
% sizes can lead to algorithm instability. The variable alpha below
% specifies the fixed step size. Increasing alpha above 0.32 results in
% instability of the algorithm. An alternative approach would involve a
% variable step size determined through line search.
%
% This example was used originally for an optimization demonstration in ME
% 149, Engineering System Design Optimization, a graduate course taught at
% Tufts University in the Mechanical Engineering Department. A
% corresponding video is available at:
%
% http://www.youtube.com/watch?v=cY1YGQQbrpQ
%
% Author: James T. Allison, Assistant Professor, University of Illinois at
% Urbana-Champaign
% Date: 3/4/12

if nargin==0
    % define starting point
    x0 = [3 3]';
elseif nargin==1
    % if a single input argument is provided, it is a user-defined starting
    % point.
    x0 = varargin{1};
else
    error('Incorrect number of input arguments.')
end

% termination tolerance
tol = 1e-6;

% maximum number of allowed iterations
maxiter = 1000;

% minimum allowed perturbation
dxmin = 1e-6;

% step size ( 0.33 causes instability, 0.2 quite accurate)
alpha = 0.1;

% initialize gradient norm, optimization vector, iteration counter, perturbation
gnorm = inf; x = x0; niter = 0; dx = inf;

% define the objective function:
f = @(x1,x2) x1.^2 + x1.*x2 + 3*x2.^2;

% plot objective function contours for visualization:
figure(1); clf; ezcontour(f,[-5 5 -5 5]); axis equal; hold on

% redefine objective function syntax for use with optimization:
f2 = @(x) f(x(1),x(2));

% gradient descent algorithm:
while and(gnorm>=tol, and(niter <= maxiter, dx >= dxmin))
    % calculate gradient:
    g = grad(x);
    gnorm = norm(g);
    % take step:
    xnew = x - alpha*g;
    % check step
    if ~isfinite(xnew)
        display(['Number of iterations: ' num2str(niter)])
        error('x is inf or NaN')
    end
    % plot current point
    plot([x(1) xnew(1)],[x(2) xnew(2)],'ko-')
    refresh
    % update termination metrics
    niter = niter + 1;
    dx = norm(xnew-x);
    x = xnew;

end
xopt = x;
fopt = f2(xopt);
niter = niter - 1;

% define the gradient of the objective
function g = grad(x)
g = [2*x(1) + x(2)
    x(1) + 6*x(2)];

author: Sable

repo: mcbench-benchmarks


text:
function scatterbar3(X,Y,Z,width)
%SCATTERBAR3   3-D scatter bar graph.
%   SCATTERBAR3(X,Y,Z,WIDTH) draws 3-D bars of height Z at locations X and Y with width WIDTH.
%
%   X, Y and Z must be of equal size.  If they are vectors, than bars are placed
%   in the same fashion as the SCATTER3 or PLOT3 functions.
%
%   If they are matrices, then bars are placed in the same fashion as the SURF
%   and MESH functions.
%
%   The colors of each bar read from the figure's colormap according to the bar's height.
%
%   NOTE:  For best results, you should use the 'zbuffer' renderer.  To set the current
%   figure renderer to 'zbuffer' use the following command:
%
%       set(gcf,'renderer','zbuffer')
%
%    % EXAMPLE 1:
%    y=[1 2 3 1 2 3 1 2 3];
%    x=[1 1 1 2 2 2 3 3 3];
%    z=[1 2 3 6 5 4 7 8 9];
%    scatterbar3(x,y,z,1)
%    colorbar
%
%    % EXAMPLE 2:
%    [X,Y]=meshgrid(-1:0.25:1);
%    Z=2-(X.^2+Y.^2);
%    scatterbar3(X,Y,Z,0.2)
%    colormap(hsv)
%
%    % EXAMPLE 3:
%    t=0:0.1:(2*pi);
%    x=cos(t);
%    y=sin(t);
%    z=sin(t);
%    scatterbar3(x,y,z,0.07)

% By Mickey Stahl - 2/25/02
% Engineering Development Group
% Aspiring Developer

[r,c]=size(Z);
for j=1:r,
    for k=1:c,
        if ~isnan(Z(j,k))
            drawbar(X(j,k),Y(j,k),Z(j,k),width/2)
        end
    end
end

zlim=[min(Z(:)) max(Z(:))];
if zlim(1)>0,zlim(1)=0;end
if zlim(2)<0,zlim(2)=0;end
axis([min(X(:))-width max(X(:))+width min(Y(:))-width max(Y(:))+width zlim])
caxis([min(Z(:)) max(Z(:))])

function drawbar(x,y,z,width)

h(1)=patch([-width -width width width]+x,[-width width width -width]+y,[0 0 0 0],'b');
h(2)=patch(width.*[-1 -1 1 1]+x,width.*[-1 -1 -1 -1]+y,z.*[0 1 1 0],'b');
h(3)=patch(width.*[-1 -1 -1 -1]+x,width.*[-1 -1 1 1]+y,z.*[0 1 1 0],'b');
h(4)=patch([-width -width width width]+x,[-width width width -width]+y,[z z z z],'b');
h(5)=patch(width.*[-1 -1 1 1]+x,width.*[1 1 1 1]+y,z.*[0 1 1 0],'b');
h(6)=patch(width.*[1 1 1 1]+x,width.*[-1 -1 1 1]+y,z.*[0 1 1 0],'b');
set(h,'facecolor','flat','FaceVertexCData',z)


text:
% lssat;
% check saturation condition
%

cont=saturated;

if cont,
  % check boundary minimizer
  [fmin,i]=min(flist);
  if i==1 | i==s, cont=0; end;
end;

if cont,
  % select three points for parabolic interpolation
  aa=alist([i-1:i+1]);ff=flist([i-1:i+1]);

  % get divided differences
  f12=(ff(2)-ff(1))/(aa(2)-aa(1));
  f23=(ff(3)-ff(2))/(aa(3)-aa(2));
  f123=(f23-f12)/(aa(3)-aa(1));

  if f123>0,
    % parabolic minimizer
    alp=0.5*(aa(2)+aa(3)-f23/f123);
    alp=max(amin,min(alp,amax));
    alptol=small*(aa(3)-aa(1));
    saturated= ( abs(alist(i)-alp)<=alptol );
  else
    saturated=0;
  end;

  if prt>1 & ~saturated,
    disp('saturation check negative')
  end;
end;

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function [nst nst_TI nst_TV]=st_size(errdefs)

nsen=length(errdefs);

nst=0;
nst_TV=0;
nst_TI=0;
for in=1:nsen
    nst_TI=nst_TI+size(errdefs(in).A,1);
end

for in=1:nsen
    if (~isempty(errdefs(in).tparam))
        nst_TV=nst_TV+2;
    end
end

nst=nst_TV+nst_TI;

repo: ekfmonoslam


text:
function [surface_handles, ax] = canlab_canonical_brain_surface_cutaways(method_keyword, varargin)
% Create one of a number of pre-set 3D brain views, and return surface handles for rendering
%
% :Usage:
% ::
%
% [surface_handles, axis_handles] = canlab_canonical_brain_surface_cutaways(method_keyword, ['noverbose'])
%
% - Uses fmri_data.isosurface method, which is a flexible way of creating surfaces
% - You can then render activation blobs on these surfaces using the
%   surface() object methods for CANlab objects or cluster_surf
%
% - Uses a pre-set group-average anatomical image tuned for quality visual display
%   ('keuken_2014_enhanced_for_underlay.img') from Keuken et al. 7T atlas
%
%   Keuken et al. (2014). Quantifying inter-individual anatomical variability in the subcortex using 7T structural MRI
%   Forstmann et al. (2014). Forstmann, Birte U., Max C. Keuken, Andreas Schafer, Pierre-Louis Bazin, Anneke Alkemade, and Robert Turner. 2014. ?Multi-Modal Ultra-High Resolution Structural 7-Tesla MRI Data Repository.? Scientific Data 1 (December): 140050.
%
% - This function is called from addbrain, which has many options for
%   surface rendering.
%
% - fmri_data.isosurface() can be used to create many more custom surfaces
%   See this code for examples.
%
% :Inputs:
%
%   **method_keyword:** One of the following:
%     'left_cutaway'
%     'right_cutaway'
%     'left_insula_slab'
%     'right_insula_slab'
%     'accumbens_slab'
%
% :Outputs:
%   **A figure display with rendering**
%
%   **surface_handles:**
%        A vector of handles to surface objects (isosurfaces and isocaps)
%
%  % To recreate these isosurfaces, do this:
% anat = fmri_data(which('keuken_2014_enhanced_for_underlay.img'), 'noverbose');
% [isosurf, isocap] = deal({});
% [p, ~, isosurf{1}, isocap{1}] = isosurface(anat, 'thresh', 140, 'nosmooth', 'ylim', [-Inf -30]);
% [p2, ~, isosurf{2}, isocap{2}] = isosurface(anat, 'thresh', 140, 'nosmooth', 'xlim', [-Inf 0], 'YLim', [-30 Inf]);
%
% :See also:
%   fmri_data.isosurface, addbrain, cluster_cutaways, fmri_data.surface,
%   region.surface
%

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

% load this to create new isosurfaces, but not needed to load existing saved ones:
% anat = fmri_data(which('keuken_2014_enhanced_for_underlay.img'), 'noverbose');

ax = gca;

doverbose = true;
if any(strcmp(varargin, 'noverbose')), doverbose = false; end

anat.source_notes = {};
anat.source_notes{1, 1} = 'Anatomical image: Adapted from the 7T high-resolution atlas of Keuken, Forstmann et al.';
anat.source_notes{2, 1} = 'Keuken et al. (2014). Quantifying inter-individual anatomical variability in the subcortex using 7T structural MRI. ';
anat.source_notes{3, 1} = 'Forstmann, Birte U., Max C. Keuken, Andreas Schafer, Pierre-Louis Bazin, Anneke Alkemade, and Robert Turner. 2014. ?Multi-Modal Ultra-High Resolution Structural 7-Tesla MRI Data Repository.? Scientific Data 1 (December): 140050.';

switch method_keyword

    case {'cutaway' 'cutaways' 'left_cutaway'}

        % To recreate these isosurfaces, do this:
        %         [isosurf, isocap] = deal({});
        %         [p, ~, isosurf{1}, isocap{1}] = isosurface(anat, 'thresh', 140, 'nosmooth', 'ylim', [-Inf -30]);
        %         [p2, ~, isosurf{2}, isocap{2}] = isosurface(anat, 'thresh', 140, 'nosmooth', 'xlim', [-Inf 0], 'YLim', [-30 Inf]);
        %
        %         save('canlab_canonical_brain_surface_left_cutaway', 'isosurf', 'isocap');

        load('canlab_canonical_brain_surface_left_cutaway', 'isosurf', 'isocap');
        p = draw_isosurface_patches(isosurf, isocap, [.5 .5 .5]);

        p3 = addbrain('limbic hires');
        alpha 1
        delete(p3(3)); p3(3) = [];
        set(p3, 'FaceAlpha', .6, 'FaceColor', [.5 .5 .5]);

        set_slab_view_properties(p, [135 15])

        surface_handles = [p p3];

    case 'right_cutaway'

        %         [isosurf, isocap] = deal({});
        %         [p, ~, isosurf{1}, isocap{1}] = isosurface(anat, 'thresh', 140, 'nosmooth', 'ylim', [-Inf -30]);
        %         [p2, ~, isosurf{2}, isocap{2}] = isosurface(anat, 'thresh', 140, 'nosmooth', 'xlim', [0 Inf], 'yLim', [-30 Inf]);
        %         save('canlab_canonical_brain_surface_right_cutaway', 'isosurf', 'isocap');

        load('canlab_canonical_brain_surface_right_cutaway', 'isosurf', 'isocap');
        p = draw_isosurface_patches(isosurf, isocap, [.5 .5 .5]);

        p3 = addbrain('limbic hires');
        alpha 1
        delete(p3(3)); p3(3) = [];
        set(p3, 'FaceAlpha', .6, 'FaceColor', [.5 .5 .5]);

        set_slab_view_properties(p, [-137 14])

        surface_handles = [p p3];

    case 'left_insula_slab'

        %         [isosurf, isocap] = deal({});
        %         [p, ~, isosurf{1}, isocap{1}] = isosurface(anat, 'thresh', 140, 'nosmooth', 'xlim', [15 44], 'yLim', [-45 Inf]);
        %          save('canlab_canonical_brain_surface_left_insula', 'isosurf', 'isocap');
        %
        load('canlab_canonical_brain_surface_left_insula', 'isosurf', 'isocap');
        p = draw_isosurface_patches(isosurf, isocap, [.5 .5 .5]);

        set_slab_view_properties(p, [271 0])
        surface_handles = p;

    case 'right_insula_slab'

        %         [isosurf, isocap] = deal({});
        %         [p, ~, isosurf{1}, isocap{1}] = isosurface(anat, 'thresh', 140, 'nosmooth', 'xlim', [-15 -44], 'YLim', [-45 Inf]);
        %         save('canlab_canonical_brain_surface_right_insula', 'isosurf', 'isocap');
        %
        load('canlab_canonical_brain_surface_right_insula', 'isosurf', 'isocap');
        p = draw_isosurface_patches(isosurf, isocap, [.5 .5 .5]);

        set_slab_view_properties(p, [91 0])
        surface_handles = p;

    case 'accumbens_slab'

        %         [isosurf, isocap] = deal({});
        %         [p, ~, isosurf{1}, isocap{1}] = isosurface(anat, 'thresh', 140, 'nosmooth', 'yLim', [-18 8]);
        %         save('canlab_canonical_brain_surface_accumbens_slab', 'isosurf', 'isocap');
        % %
        load('canlab_canonical_brain_surface_accumbens_slab', 'isosurf', 'isocap');
        p = draw_isosurface_patches(isosurf, isocap, [.5 .5 .5]);

        set_slab_view_properties(p, [156 0])
        surface_handles = p;

    case 'coronal_slabs'

        % surface_handles = render_coronal_slabs(anat, 6);

        [surface_handles, ax] = draw_coronal_slab_patches('canlab_canonical_brain_surface_coronal_slabs6');

    case 'coronal_slabs_4'

        % surface_handles = render_coronal_slabs(anat, 4);

        [surface_handles, ax] = draw_coronal_slab_patches('canlab_canonical_brain_surface_coronal_slabs4');

    case 'coronal_slabs_5'

        % surface_handles = render_coronal_slabs(anat, 5);

        [surface_handles, ax] = draw_coronal_slab_patches('canlab_canonical_brain_surface_coronal_slabs5');

    otherwise
        error('illegal method keyword.');

end % switch

if doverbose
    disp(char(anat.source_notes));
end

end % main function

function set_slab_view_properties(p, az_el)

set(p, 'FaceAlpha', 1);
colormap(gca, gray);
axis vis3d image
view(az_el(1), az_el(2));
lightRestoreSingle;

end

function [p, isosurf, isocap] = render_coronal_slabs(anat, nslabs)
% Use this to create custom isosurface slabs

start_mm = -80;
end_mm = 50;
slab_thickness_mm = 20;

st = linspace(start_mm, end_mm - slab_thickness_mm, nslabs); % define starting points

axwid = 2*0.8 / (nslabs + 1);
xstart = .1;
ystart = .9 - axwid;
xpos = linspace(xstart, .65, nslabs);
ypos = linspace(ystart, .7 - axwid, nslabs);

p = {};

for i = 1:nslabs
    ax = axes('Position', [xpos(i) ypos(i) axwid axwid]);
    [p{i}, ~, isosurf{i}, isocap{i}] = isosurface(anat, 'thresh', 140, 'nosmooth', 'yLim', [st(i) st(i)+slab_thickness_mm]);
    set(ax, 'Color', 'none');
    set_slab_view_properties(p{i}, [135 12])
    axis off
    drawnow
end

p = cat(2, p{:});

% *pause and save file here if creating stock slabs*
end

function [p, ax] = draw_coronal_slab_patches(isosurfacefilename)
% Use this to draw patches from saved isosurface slabs

load(isosurfacefilename, 'isosurf', 'isocap')
nslabs = length(isosurf);

start_mm = -80;
end_mm = 50;
slab_thickness_mm = 20;

st = linspace(start_mm, end_mm - slab_thickness_mm, nslabs); % define starting points

axwid = 2*0.8 / (nslabs + 1);
xstart = .1;
ystart = .9 - axwid;
xpos = linspace(xstart, .65, nslabs);
ypos = linspace(ystart, .7 - axwid, nslabs);

p = {};

for i = 1:nslabs
    ax = axes('Position', [xpos(i) ypos(i) axwid axwid]);
    %[p{i}, ~, isosurf{i}, isocap{i}] = isosurface(anat, 'thresh', 140, 'nosmooth', 'yLim', [st(i) st(i)+slab_thickness_mm]);
    p{i} = draw_isosurface_patches(isosurf(i),  isocap(i), [.5 .5 .5]);
    set(ax, 'Color', 'none');
    set_slab_view_properties(p{i}, [135 12])
    axis off
    drawnow
end

p = cat(2, p{:});
end

function p = draw_isosurface_patches(my_isosurface,  my_isocap, mycolor)

p = [];

for i = 1:length(my_isosurface) % cell
    p(end + 1) = patch('Faces',my_isosurface{i}.faces,'Vertices',my_isosurface{i}.vertices,'FaceColor', mycolor, ...
        'EdgeColor','none','SpecularStrength', .2,'FaceAlpha', .3,'SpecularExponent', 200);
end

for i = 1:length(my_isocap) % cell
    p(end + 1) = patch(my_isocap{i}, 'FaceColor', 'interp','EdgeColor', 'none', 'FaceAlpha',1);
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function output = exportAsOpenQuakeCatalogueCsv(catalog, filename)
    % export a Catalog to a csv file readable as a Catalogue
    if isa(catalog, 'ZmapCatalog')
        catalog = catalog.table();
    end
    if ~istable(catalog)
        error('expected a ZmapCatalog')
    end
    ymdHMS = datevec(catalog.Date);
    output = table;
    output.year = ymdHMS(:,1);
    output.month = ymdHMS(:,2);
    output.day = ymdHMS(:,3);
    output.hour = ymdHMS(:,4);
    output.minute = ymdHMS(:,5);
    output.second = ymdHMS(:,6);
    output.latitude = catalog.Latitude;
    output.longitude = catalog.Longitude;
    output.depth = catalog.Depth;
    output.magnitude = catalog.Magnitude;
    % output.moment = 10.^(1.5*magnitude + 9.05);
    output.eventID = "id"+string(1:height(output))';
    output.magnitudeType = catalog.MagnitudeType;
    if exist('filename','var') && (~exist(filename,'file') ||...
            questdlg(['File ',filename,' exists, overwrite?'],'export catalog')=="Yes")
        writetable(output,filename);
    end
end

author: CelsoReyes


text:
3     % problem
4     % grid parameter
1     % uniform/stretched grid
2     % domain size
1     % approximation

%% Data file for test problem P3

author: johannesgerer


text:
function value = r4_sin_deg ( x )

%*****************************************************************************80
%
%% R4_SIN_DEG evaluates the sine of an R4 argument in degrees.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 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 in degrees.
%
%    Output, real VALUE, the sine of X.
%
  raddeg = 0.017453292519943296;

  value = sin ( raddeg * x );

  if ( mod ( x, 90.0 ) == 0.0 )

    n = floor ( abs ( x ) / 90.0 + 0.5 );
    n = mod ( n, 2 );

    if ( n == 0 )
      value = 0.0;
    elseif ( value < 0.0 )
      value = - 1.0;
    else
      value = + 1.0;
    end

  end

  return
end


text:
% emacs, this is -*-Matlab-*- mode

% The main script. Performes the self calibration
% The point coordinates are expected to be known
% see the directory FindingPoints
%
% $Author: svoboda $
% $Revision: 2.7 $
% $Id: gocal.m,v 2.7 2005/05/24 09:15:11 svoboda Exp $
% $State: Exp $

clear variables globals

Octave = exist('OCTAVE_VERSION', 'builtin') ~= 0;

% add necessary paths
addpath ('../CommonCfgAndIO')
addpath ('../RadialDistortions')
addpath ('./CoreFunctions')
addpath ('./OutputFunctions')
addpath ('./BlueCLocal')
addpath ('./LocalAlignments')
addpath ('../CalTechCal')
addpath ('../RansacM'); % ./Ransac for mex functions (it is significantly faster for noisy data)
% get the configuration
config = read_configuration();
disp('Multi-Camera Self-Calibration, Tomas Svoboda et al., 07/2003')
disp('************************************************************')
disp(sprintf('Experiment name: %s',expname))

%%%
% how many cameras to be filled
% if 0 then only points visible in all cameras will be used for selfcal
% the higher, the more points -> better against Gaussian noise
% however the higher probability of wrong filling
% Then the iterative search for outliers takes accordingly longer
% However, typically no more than 5-7 iterations are needed
% this number should correspond to the total number of the cameras
NUM_CAMS_FILL  = config.cal.NUM_CAMS_FILL;
%%%
% tolerance for inliers. The higher uncorrected radial distortion
% the higher value. For BlueC cameras set to 2 for the ViRoom
% plastic cams, set to 4 (see FINDINL)
INL_TOL = config.cal.INL_TOL;
%%%
% Use Bundle Adjustment to refine the final (after removing outliers) results
% It is often not needed at all
DO_BA = config.cal.DO_BA;

UNDO_RADIAL = config.cal.UNDO_RADIAL;		% undo radial distortion, parameters are expected to be available
BA_RADIAL = config.cal.BA_RADIAL;		% bundle adjustment
                                                % also finds
                                                % non-linear parameters
SAVE_STEPHI = 1;		% save calibration parameters in Stephi's Carve/BlueC compatible form
SAVE_PGUHA	= 1;		% save calib pars in Prithwijit's compatible form
USED_MULTIPROC = 0;		% was the multipropcessing used?
						% if yes then multiple IdMat.dat and points.dat have to be loaded
						% setting to 1 it forces to read the multiprocessor data against the
						% monoprocessor see the IM2POINTS, IM2PMULTIPROC.PL

%%%
% Data structures
% lin.* corrected values which obey linear model
% in.* inliers, detected by a chain application of Ransac

if strfind(expname,'oscar')
  % add a projector idx to the cameras
  % they are handled the same
  config.files.cams2use = [config.files.idxcams,config.files.idxproj];
end

selfcal.par2estimate = config.cal.nonlinpar;
selfcal.iterate = 1;
selfcal.count	= 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load initial distortion parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nl_params_all_cams = [];
if BA_RADIAL
  for i=1:CAMS,
    if UNDO_RADIAL
      [K,kc] = ...
	  readradfile(sprintf(config.files.rad,config.cal.cams2use(i)));
    else
      % no radial distortion
      K = [ 1 0 config.cal.Res(i,1)/2; ...
	    0 1 config.cal.Res(i,2)/2;
	0 0 1];
      kc = [0,0,0,0];
    end
    cam_pvec = rad2pvec(K,kc); % convert all NL params to row vector
    nl_params_all_cams(i) = cam_pvec; % append row vector to matrix
  end
end

if ~isempty(nl_params_all_cams)
  % This error is here to prevent a long run just to be confronted
  % with it in bundle_PX_proj. Remove this test when support is
  % implemented.
  error('Not implemented: radial distortion in bundle adjustment');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main global cycle begins
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while selfcal.iterate && selfcal.count < config.cal.GLOBAL_ITER_MAX,
  % read the input data
  loaded = loaddata(config);
  linear = loaded;		% initalize the linear structure

  CAMS = size(config.cal.cams2use,2);
  FRAMES = size(loaded.IdMat,2);

  if CAMS < 3 || FRAMES < 20
	error('gocal: Not enough cameras or images -> Problem in loading data?')
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%% correct the required amount of cameras to be filled
  if CAMS-NUM_CAMS_FILL < 3
	NUM_CAMS_FILL = CAMS-3;
  end

  config.cal.Res= loaded.Res;
  config.cal.pp = reshape([loaded.Res./2,zeros(size(loaded.Res(:,1)))]',CAMS*3,1);
  config.cal.pp = [loaded.Res./2,zeros(size(loaded.Res(:,1)))];

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % See the README how to compute data
  % for undoing of the radial distortion
  if UNDO_RADIAL
	for i=1:CAMS,
	  [K,kc] = readradfile(sprintf(config.files.rad,config.cal.cams2use(i)));
	  xn	 = undoradial(loaded.Ws(i*3-2:i*3,:),K,[kc,0]);
	  linear.Ws(i*3-2:i*3,:) = xn;
	end
	linear.Ws = linear.Ws - repmat(reshape(config.cal.pp',CAMS*3,1), 1, FRAMES);
  elseif config.cal.UNDO_HEIKK,
	for i=1:CAMS,
	  heikkpar = load(sprintf(config.files.heikkrad,config.cal.cams2use(i)),'-ASCII');
	  xn = undoheikk(heikkpar(1:4),heikkpar(5:end),loaded.Ws(i*3-2:i*3-1,:)');
	  linear.Ws(i*3-2:i*3-1,:) = xn';
	end
	linear.Ws = linear.Ws - repmat(reshape(config.cal.pp',CAMS*3,1), 1, FRAMES);
  else
	linear.Ws = loaded.Ws - repmat(reshape(config.cal.pp',CAMS*3,1), 1, FRAMES);
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Detection of outliers
  % RANSAC is pairwise applied
  disp(' ')
  disp(sprintf('********** After %d iteration *******************************************',selfcal.count))
  %  disp('****************************************************************************')
  disp(sprintf('RANSAC validation step running with tolerance threshold: %2.2f ...',INL_TOL));

  inliers.IdMat = findinl(linear.Ws,linear.IdMat,INL_TOL);

  addpath ('./MartinecPajdla');
  setpaths;		% set paths for M&P algorithms

  % remove zero-columns or just 1 point columns
  % create packed represenatation
  % it is still a bit tricky, the number of the minimum number of cameras
  % are specified here, may be some automatic method would be useful
  packed.idx	  = find(sum(inliers.IdMat)>=size(inliers.IdMat,1)-NUM_CAMS_FILL);
  packed.IdMat  = inliers.IdMat(:,packed.idx);
  packed.Ws	  = linear.Ws(:,packed.idx);

  if size(packed.Ws,2)<20
	error(sprintf('Only %d points survived RANSAC validation and packing: probably not enough points for reliable selfcalibration',size(packed.Ws,2)));
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%% fill cam(i) structures
  for i=1:CAMS,
	cam(i).camId  = config.cal.cams2use(i);
	cam(i).idlin  = find(linear.IdMat(i,:)); % loaded structure
	cam(i).idin	  = find(inliers.IdMat(i,:));	% survived initial pairwise validation
	cam(i).xdist  = loaded.Ws(3*i-2:3*i,cam(i).idlin);  % original distorted coordinates
	cam(i).xgt	  = linear.Ws(3*i-2:3*i,cam(i).idlin);
	cam(i).xgtin  = linear.Ws(3*i-2:3*i,cam(i).idin);
	% convert the ground truth coordinates by using the known principal point
	cam(i).xgt    = cam(i).xgt + repmat(config.cal.pp(i,:)', 1, size(cam(i).xgt,2));
	cam(i).xgtin  = cam(i).xgtin + repmat(config.cal.pp(i,:)', 1, size(cam(i).xgtin,2));
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%% options for the Martinec-Pajdla filling procedure
  options.verbose = 0;
%  options.strategy = 1; % force a particular central camera
  options.no_BA = 1;
  options.iter = 5;
  options.detection_accuracy = 2;
  options.consistent_number  = 9;
  options.consistent_number_min = 6;
  options.samples = 1000; %1000; %10000;
  options.sequence = 0;
  options.create_nullspace.trial_coef = 10; %20;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%% start of the *compute and remove outliers cycle*
  outliers = 1;
  inliers.idx = packed.idx;
  inliers.reprerr = [9e9];
  while outliers
	disp(sprintf('%d points/frames have survived validations so far',size(inliers.idx,2)))
	disp('Filling of missing points is running ...')
	[P,X, u1,u2, info] = fill_mm_bundle(linear.Ws(:,inliers.idx),config.cal.pp(:,1:2)',nl_params_all_cams,options);
	%
	Rmat = P*X;
	Lambda = Rmat(3:3:end,:);
	%
	[Pe,Xe,Ce,Re] = euclidize(Rmat,Lambda,P,X,config);
	disp('************************************************************')
	%
	% compute reprojection errors
	cam = reprerror(cam,Pe,Xe,FRAMES,inliers);
	%
	% detect outliers in cameras
	[outliers,inliers] = findoutl(cam,inliers,INL_TOL,NUM_CAMS_FILL);
	%
	disp(sprintf('Number of detected outliers: %3d',outliers))
	disp('About cameras (Id, 2D reprojection error, #inliers):')
	dispcamstats(cam,inliers);
	% [[cam(:).camId]',[cam(:).std2Derr]',[cam(:).mean2Derr]', sum(inliers.IdMat')']
	disp('***************************************************************')
	%
	% do BA after removing very bad outliers or if the process starts to diverge
	% and only if required config.cal.START_BA
	inliers.reprerr = [inliers.reprerr, mean([cam(:).mean2Derr])];
	if inliers.reprerr(end)<5*INL_TOL || inliers.reprerr(end-1)<inliers.reprerr(end),
	   try, options.no_BA = ~config.cal.START_BA; catch, options.no_BA = 1; end % 1 0
	end
  end
  %%% end of the cycle
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%% Do the final refinement through the BA if required and if not
  %%% performed during the iteration steps
  if (DO_BA && ~config.cal.START_BA) || (DO_BA && config.cal.START_BA && size(inliers.reprerr,2)<3)
	disp('**************************************************************')
	disp('Refinement by using Bundle Adjustment')
	options.no_BA = 0;
	[P,X, u1,u2, info] = fill_mm_bundle(linear.Ws(:,inliers.idx),config.cal.pp(:,1:2)',nl_params_all_cams,options);
	Rmat = P*X;
	Lambda = Rmat(3:3:end,:);
	[in.Pe,in.Xe,in.Ce,in.Re] = euclidize(Rmat,Lambda,P,X,config);
	cam = reprerror(cam,in.Pe,in.Xe,FRAMES,inliers);
	% [outliers,inliers] = findoutl(cam,inliers,INL_TOL,NUM_CAMS_FILL);
  else
	in.Pe = Pe; in.Xe = Xe; in.Ce = Ce; in.Re = Re;
  end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if ~Octave,
	% plot reconstructed cameras and points
	drawscene(Xe,Ce,Re,3,'cloud','reconstructed points/camera setup');
	drawscene(in.Xe,in.Ce,in.Re,4,'cloud','reconstructed points/camera setup only inliers are used',config.cal.cams2use);
	% plot measured and reprojected 2D points
	for i=1:CAMS
	  in.xe		= in.Pe(((3*i)-2):(3*i),:)*in.Xe;
	  cam(i).inxe	= in.xe./repmat(in.xe(3,:),3,1);
	  figure(i+10)
	  clf
	  plot(cam(i).xdist(1,:),cam(i).xdist(2,:),'go');
	  hold on, grid on
	  plot(cam(i).xgt(1,:),cam(i).xgt(2,:),'ro');
	  plot(cam(i).xgtin(1,:),cam(i).xgtin(2,:),'bo');
	  % plot(cam(i).xe(1,:),cam(i).xe(2,:),'r+')
	  plot(cam(i).inxe(1,:),cam(i).inxe(2,:),'k+','MarkerSize',7)
	  %plot(xe(1,:),xe(2,:),'r+','linewidth',3,'MarkerSize',10)
	  title(sprintf('o: measured (green=distorted, blue=inlier, red=discarded). +: reprojected. (camera: %d)',config.cal.cams2use(i)));
	  for j=1:size(cam(i).visandrec,2); % plot the reprojection errors
		line([cam(i).xgt(1,cam(i).visandrec(j)), ...
		      cam(i).inxe(1,cam(i).recandvis(j))], ...
		    [cam(i).xgt(2,cam(i).visandrec(j)), ...
		      cam(i).inxe(2,cam(i).recandvis(j))],'Color','g');
	  end
	  % draw the image border
	  line([0 0 0 2*config.cal.pp(i,1) 2*config.cal.pp(i,1) 2*config.cal.pp(i,1) 2*config.cal.pp(i,1) 0],[0 2*config.cal.pp(i,2) 2*config.cal.pp(i,2) 2*config.cal.pp(i,2) 2*config.cal.pp(i,2) 0 0 0],'Color','k','LineWidth',2,'LineStyle','--')
	  axis('equal')
	  drawnow
	  eval(['print -depsc ', config.paths.data, sprintf('%s%d.reprojection.eps',config.files.basename,cam(i).camId)])
	end
  end

  %%%
  % SAVE camera matrices
  P = in.Pe;
  X = in.Xe;
  R = in.Re;
  C = in.Ce;
  if Octave
    % all Octave data in ASCII format
    save(config.files.Pmats,'P');
    save(config.files.Xe,'X');
    save(config.files.Re,'R');
    save(config.files.Ce,'C');
  else
    save(config.files.Pmats,'P','-ASCII');
    save(config.files.Xe,'X','-ASCII');
    save(config.files.Re,'R','-ASCII');
    save(config.files.Ce,'C','-ASCII');
  end

  % save normal data
  if SAVE_STEPHI || SAVE_PGUHA
	 [in.Cst,in.Rot] = savecalpar(in.Pe,config);
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % local routines for the BlueC installations
  % The main functionality of these functions that
  % they align the coordinate frame from the selfcalibration
  % with the pre-defined world frame
  % it is assumed the necessary informations are avialable
  if strfind(expname,'BlueCRZ')
	  [align] = bluecrz(in,config);
  end
  if strfind(expname,'Hoengg')
	  [align] = bluechoengg(in,config);
  end
  if strfind(expname,'Erlangen')
	  [align] = erlangen(in,config);
	  % [align] = planarmove(in,config);
  end
  if strfind(expname,'G9')
	  [align] = g9(in,config);
  end
  if strfind(expname,'flydra')
	  [align] = flydra(in,config);
  end
  if strfind(expname,'mamarama')
	  [align] = mamarama(in,config);
  end
  if strfind(expname,'humdra')
	  [align] = humdra(in,config);
  end

  if config.cal.ALIGN_EXISTING
      align_existing_camera_centers(in,config);
  end

  % planar alignement if knowledge available
  % [align,cam] = planarmove(in,cam,config);
  % try, [align,cam] = planarcams(in,cam,config,config.cal.planarcams); end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Evaluate reprojection error
  %%%
  cam = evalreprerror(cam,config);

  %%%%
  % Save the 2D-3D correpondences for further processing
  for i=1:CAMS,
	xe = cam(i).xdist(1:3,cam(i).visandrec);	% save the original distorted coordinates
 	% save the reconstructed points (aligned if available)
	try, Xe = align.X(1:4,cam(i).recandvis); catch, Xe = in.Xe(1:4,cam(i).recandvis); end;
	% Xe = in.Xe(1:4,cam(i).recandvis);
	corresp = [Xe',xe'];
	if Octave
	  % all Octave data in ASCII format
          save(sprintf(config.files.points4cal,config.cal.cams2use(i)),'corresp');
	else
          save(sprintf(config.files.points4cal,config.cal.cams2use(i)),'corresp','-ASCII');
	end
  end

%%%
% TO-DO:
% - find a suitable end condition for the global iteration.
%   This threshold may very depend on local conditions
%
% - how to check meaningful number of iterations
%   typically only few iterations are needed
%
% - The precision of the non-linear estimation should be somehow taken into account

  selfcal.count = selfcal.count+1;

  if max([cam.mean2Derr])>config.cal.GLOBAL_ITER_THR && config.cal.DO_GLOBAL_ITER && selfcal.count < config.cal.GLOBAL_ITER_MAX
    % if the maximal reprojection error is still bigger
    % than acceptable estimate radial distortion and
    % iterate further
    cd ../CalTechCal
    selfcalib = goradf(config,selfcal.par2estimate,INL_TOL);
    cd ../MultiCamSelfCal
    selfcal.iterate = 1;
    UNDO_RADIAL = 1;
    if ~selfcalib.goradproblem;
      % if all non-linear parameters estimated reliable
      % we can reduce the tolerance threshold
      INL_TOL = max([(2/3)*INL_TOL,config.cal.GLOBAL_ITER_THR]);
      % add the second radial distortion parameter
      if config.cal.NL_UPDATE(4), selfcal.par2estimate(4) = 1; end
      % estimate also the principal point
      if selfcal.count > 1 && config.cal.NL_UPDATE(2), selfcal.par2estimate(2) = 1; end
      % estimate also the tangential distortion
      if selfcal.count > 3 && all(config.cal.NL_UPDATE(5:6)), selfcal.par2estimate(5:6) = 1; end
    else
      INL_TOL = min([3/2*INL_TOL,config.cal.INL_TOL]);
    end
  else
    % ends the iteration
    % the last computed parameters will be taken as valid
    selfcal.iterate = 0;
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of the main global cycle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: strawlab

repo: MultiCamSelfCal


text:
function declus_SingleCatalog
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Example: mDeclusIndex=declus_SingleCatalog
%
% This function to calculate declustering probabilities
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Author:
% van Stiphout, Thomas, vanstiphout@sed.ethz.ch
%
% Created on 02.07.2008
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

initialize
disp('~/zmap/src/thomas/decluster/declus_SingleCatalog');

% load catalog to work with (can also be selected later
load anss-1981-19923-M15-Wiemer1994.mat
% load '~/zmap/eq_data/lsh_1.12/lsh1.12WiemerWyss1994.mat';
%  output filename
sFilename='mDeclus08072800.mat';

if ~exist('mCatalog')
  mCatalog=a;
  clear a;
end
bHypo=true;
bMag=true;

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

nMCS=2;            % number of simulation for hypocenter and/or magnitude and/or declustering

nDeclusMode=2;   % Applied declustering algorithm
% 1: Reasenberg Matlab
% 2: Gardner Knopoff
% 3: Stochastic
% 4: Reasenberg Cluster2000 (fortran)
% 5: Marsan
% 6: cluster GK written by Annemarie
% 7: cluster Uhrhammer written by Annemarie
% 8: cluster Utsu written by Annemarie

% parameter range for reasenberg declustering (nDeclusMode=4)
% Taumin, Taumax, P1, Xk, Xmeff, Rfact, Err, Derr
    mReasenParam=[[1 1];[10 10];[0.95 0.95];[0.5 0.5];...
        [1.6 1.6];[10 10];[2 2];[4 4]];
% mReasenParam=[[0.5 2.5];[03 15];[0.90 0.99];[0.0 1.0];...
%    [1.6 1.8];[05 20];[2 2];[4 4]];

fMc=3.5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% preallocate matrices
mCatalog=mCatalog(mCatalog(:,3)<1992.3,:);
mNumDeclus=zeros(size(mCatalog,1),nMCS);
mDeclusLon=nan(size(mCatalog,1),nMCS);
mDeclusLat=nan(size(mCatalog,1),nMCS);
mDeclusDepth=nan(size(mCatalog,1),nMCS);
mDeclusMag=nan(size(mCatalog,1),nMCS);

nMCS
for i=1:nMCS
    mCat=mCatalog;
    %% include hypocenter uncertainties
    if bHypo
%         mHypo=[2 2 5];mDelta=repmat(mHypo,size(mCatalog,1),1);
        mDelta=[mCat(:,11) mCat(:,11) mCat(:,12)];
        [mCat]=calc_hyposhift(mCat,mDelta,true);
    end

    %% include magnitude uncertainties
    if bMag
%         mMag=[0.1];mDeltaMag=repmat(mMag,size(mCatalog,1),1);
        mDeltaMag=mCat(:,13);
        [mCat]=calc_magerr(mCat,mDeltaMag);
    end

    %% cut catalog
    vSel=(mCat(:,6)>fMc);
%     mCat=mCat{i}(vSel{i},:)
    [declusCat, mNumDeclus(vSel,i)] = MonteDeclus(mCat(vSel,:),...
        1,nDeclusMode,mReasenParam);
    % save results of i-th simulation in matrices
    if bHypo
%         mNumDeclus(isnan(mNumDeclus(~vSel,i)),i)=0;
        mDeclusLon(logical(mNumDeclus(:,i)),i)=declusCat(:,1);
        mDeclusLat(logical(mNumDeclus(:,i)),i)=declusCat(:,2);
        mDeclusDepth(logical(mNumDeclus(:,i)),i)=declusCat(:,7);
    end
    if bMag
%         mNumDeclus(isnan(mNumDeclus(vSel,i)),i)=0;
        mDeclusMag(logical(mNumDeclus(:,i)),i)=declusCat(:,6);
    end
    disp(i);
end
mDeclusIndex=sum(mNumDeclus,2)./nMCS;
clear params vSel declusCat mCat  mDelta mDeltaMag ans i sString

sString=sprintf('save %s * -mat',sFilename);
eval(sString);
disp(sString);


text:
function setproxy(proxyname, port)
% SETPROXY - Set the proxy to use when making all remote calls with
% NCTOOLBOX (https://github.com/nctoolbox/nctoolbox)
%
% Use as: setproxy(proxyname, port)
%
% Inputs:
%  proxyname = the proxy name that all connections will be routed to
%  port      = The port on the proxy to connect through
%

% Brian Schlining
% 2010-02-16
%
% httpClient = org.apache.commons.httpclient.HttpClient();
% httpClient.getHostConfiguration().setProxy(proxyname, port);
% ucar.nc2.dataset.NetcdfDataset.setHttpClient(httpClient);

%ucar.nc2.util.net.HTTPSession.setGlobalSimpleProxy(proxyname, port)
ucar.nc2.util.net.HTTPSession.setGlobalProxy(proxyname, port)

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
report_this_filefun(mfilename('fullpath'));

% Now lets plot the color-map of the z-value
%
figure

set(gca,'visible','off','FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold',...
    'FontWeight','bold','LineWidth',1.5,...
    'Box','on','SortMethod','childorder')

rect = [0.18,  0.10, 0.7, 0.75];
rect1 = rect;

% set values greater tresh = nan
%
re4 = avm;l = r > tresh; re4(l) = zeros(1,length(find(l)))*nan;

% plot image
%
orient portrait
set(gcf,'PaperPosition', [2. 1 7.0 5.0])

axes('position',rect);hold on
pco1 = pcolor(gx,gy,(re4));
axis([ min(gx) max(gx) min(gy) max(gy)]); axis image; hold on
shading interp;
caxis([2.6 4.00])
if exist('pro', 'var')
    l = pro > 0;
    pro2 = pro;
    pro2(l) = pro2(l)*nan;
    %cs =contour(gx,gy,pro,[ 99 100],'w--');

    % cs =contour(gx,gy,pro,[95 99],'k');
    %[cs, hc] =contour(gx,gy,pro,[ 99 100],'w-');
end % if exist pro
h = hot(64);h = [  h(64:-1:1,:)]; colormap(h)
%end
if fre == 1
    caxis([fix1 fix2])
end

%brighten(0.5)
%xlabel('Distance in [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
%ylabel('depth in [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
set(gca,'XTickLabels',[])
set(gca,'YTickLabels',[  10 8 6 4 2])
set(gca,'Color',[0.9 0.9 0.9])

if exist('maex', 'var')
    pl = plot(maex,-maey,'*k');
    set(pl,'MarkerSize',6,'LineWidth',2)
end
overlay

set(gca,'visible','on','FontSize',10,'FontWeight','normal',...
    'FontWeight','normal','LineWidth',1.0,...
    'Box','on','TickDir','out')
h1 = gca;hzma = gca;

% Create a colobar
%
%h5 = colorbar('vertical');
%set(h5,'Pos',[0.25 0.25 0.03 0.25],...
%'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m,'TickDir','out')

% Make the figure visible
%

axes(h1)