text:
function data = loadubjson(fname,varargin)
%
% data=loadubjson(fname,opt)
%    or
% data=loadubjson(fname,'param1',value1,'param2',value2,...)
%
% parse a JSON (JavaScript Object Notation) file or string
%
% authors:Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% created on 2013/08/01
%
% $Id: loadubjson.m 460 2015-01-03 00:30:45Z fangq $
%
% input:
%      fname: input file name, if fname contains "{}" or "[]", fname
%             will be interpreted as a UBJSON string
%      opt: a struct to store parsing options, opt can be replaced by
%           a list of ('param',value) pairs - the param string is equivallent
%           to a field in opt. opt can have the following
%           fields (first in [.|.] is the default)
%
%           opt.SimplifyCell [0|1]: if set to 1, loadubjson will call cell2mat
%                         for each element of the JSON data, and group
%                         arrays based on the cell2mat rules.
%           opt.IntEndian [B|L]: specify the endianness of the integer fields
%                         in the UBJSON input data. B - Big-Endian format for
%                         integers (as required in the UBJSON specification);
%                         L - input integer fields are in Little-Endian order.
%
% output:
%      dat: a cell array, where {...} blocks are converted into cell arrays,
%           and [...] are converted to arrays
%
% examples:
%      obj=struct('string','value','array',[1 2 3]);
%      ubjdata=saveubjson('obj',obj);
%      dat=loadubjson(ubjdata)
%      dat=loadubjson(['examples' filesep 'example1.ubj'])
%      dat=loadubjson(['examples' filesep 'example1.ubj'],'SimplifyCell',1)
%
% license:
%     BSD, see LICENSE_BSD.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

global pos inStr len  esc index_esc len_esc isoct arraytoken fileendian systemendian

if(regexp(fname,'[\{\}\]\[]','once'))
   string=fname;
elseif(exist(fname,'file'))
   fid = fopen(fname,'rb');
   string = fread(fid,inf,'uint8=>char')';
   fclose(fid);
else
   error('input file does not exist');
end

pos = 1; len = length(string); inStr = string;
isoct=exist('OCTAVE_VERSION','builtin');
arraytoken=find(inStr=='[' | inStr==']' | inStr=='"');
jstr=regexprep(inStr,'\\\\','  ');
escquote=regexp(jstr,'\\"');
arraytoken=sort([arraytoken escquote]);

% String delimiters and escape chars identified to improve speed:
esc = find(inStr=='"' | inStr=='\' ); % comparable to: regexp(inStr, '["\\]');
index_esc = 1; len_esc = length(esc);

opt=varargin2struct(varargin{:});
fileendian=upper(jsonopt('IntEndian','B',opt));
[os,maxelem,systemendian]=computer;

jsoncount=1;
while pos <= len
    switch(next_char)
        case '{'
            data{jsoncount} = parse_object(opt);
        case '['
            data{jsoncount} = parse_array(opt);
        otherwise
            error_pos('Outer level structure must be an object or an array');
    end
    jsoncount=jsoncount+1;
end % while

jsoncount=length(data);
if(jsoncount==1 && iscell(data))
    data=data{1};
end

if(~isempty(data))
      if(isstruct(data)) % data can be a struct array
          data=jstruct2array(data);
      elseif(iscell(data))
          data=jcell2array(data);
      end
end

%%
function newdata=parse_collection(id,data,obj)

if(jsoncount>0 && exist('data','var'))
    if(~iscell(data))
       newdata=cell(1);
       newdata{1}=data;
       data=newdata;
    end
end

%%
function newdata=jcell2array(data)
len=length(data);
newdata=data;
for i=1:len
      if(isstruct(data{i}))
          newdata{i}=jstruct2array(data{i});
      elseif(iscell(data{i}))
          newdata{i}=jcell2array(data{i});
      end
end

%%-------------------------------------------------------------------------
function newdata=jstruct2array(data)
fn=fieldnames(data);
newdata=data;
len=length(data);
for i=1:length(fn) % depth-first
    for j=1:len
        if(isstruct(getfield(data(j),fn{i})))
            newdata(j)=setfield(newdata(j),fn{i},jstruct2array(getfield(data(j),fn{i})));
        end
    end
end
if(~isempty(strmatch('x0x5F_ArrayType_',fn)) && ~isempty(strmatch('x0x5F_ArrayData_',fn)))
  newdata=cell(len,1);
  for j=1:len
    ndata=cast(data(j).x0x5F_ArrayData_,data(j).x0x5F_ArrayType_);
    iscpx=0;
    if(~isempty(strmatch('x0x5F_ArrayIsComplex_',fn)))
        if(data(j).x0x5F_ArrayIsComplex_)
           iscpx=1;
        end
    end
    if(~isempty(strmatch('x0x5F_ArrayIsSparse_',fn)))
        if(data(j).x0x5F_ArrayIsSparse_)
            if(~isempty(strmatch('x0x5F_ArraySize_',fn)))
                dim=double(data(j).x0x5F_ArraySize_);
                if(iscpx && size(ndata,2)==4-any(dim==1))
                    ndata(:,end-1)=complex(ndata(:,end-1),ndata(:,end));
                end
                if isempty(ndata)
                    % All-zeros sparse
                    ndata=sparse(dim(1),prod(dim(2:end)));
                elseif dim(1)==1
                    % Sparse row vector
                    ndata=sparse(1,ndata(:,1),ndata(:,2),dim(1),prod(dim(2:end)));
                elseif dim(2)==1
                    % Sparse column vector
                    ndata=sparse(ndata(:,1),1,ndata(:,2),dim(1),prod(dim(2:end)));
                else
                    % Generic sparse array.
                    ndata=sparse(ndata(:,1),ndata(:,2),ndata(:,3),dim(1),prod(dim(2:end)));
                end
            else
                if(iscpx && size(ndata,2)==4)
                    ndata(:,3)=complex(ndata(:,3),ndata(:,4));
                end
                ndata=sparse(ndata(:,1),ndata(:,2),ndata(:,3));
            end
        end
    elseif(~isempty(strmatch('x0x5F_ArraySize_',fn)))
        if(iscpx && size(ndata,2)==2)
             ndata=complex(ndata(:,1),ndata(:,2));
        end
        ndata=reshape(ndata(:),data(j).x0x5F_ArraySize_);
    end
    newdata{j}=ndata;
  end
  if(len==1)
      newdata=newdata{1};
  end
end

%%-------------------------------------------------------------------------
function object = parse_object(varargin)
    parse_char('{');
    object = [];
    type='';
    count=-1;
    if(next_char == '$')
        type=inStr(pos+1); % TODO
        pos=pos+2;
    end
    if(next_char == '#')
        pos=pos+1;
        count=double(parse_number());
    end
    if next_char ~= '}'
        num=0;
        while 1
            str = parseStr(varargin{:});
            if isempty(str)
                error_pos('Name of value at position %d cannot be empty');
            end
            %parse_char(':');
            val = parse_value(varargin{:});
            num=num+1;
            eval( sprintf( 'object.%s  = val;', valid_field(str) ) );
            if next_char == '}' || (count>=0 && num>=count)
                break;
            end
            %parse_char(',');
        end
    end
    if(count==-1)
        parse_char('}');
    end

%%-------------------------------------------------------------------------
function [cid,len]=elem_info(type)
id=strfind('iUIlLdD',type);
dataclass={'int8','uint8','int16','int32','int64','single','double'};
bytelen=[1,1,2,4,8,4,8];
if(id>0)
    cid=dataclass{id};
    len=bytelen(id);
else
    error_pos('unsupported type at position %d');
end
%%-------------------------------------------------------------------------

function [data adv]=parse_block(type,count,varargin)
global pos inStr isoct fileendian systemendian
[cid,len]=elem_info(type);
datastr=inStr(pos:pos+len*count-1);
if(isoct)
    newdata=int8(datastr);
else
    newdata=uint8(datastr);
end
id=strfind('iUIlLdD',type);
if(id<=5 && fileendian~=systemendian)
    newdata=swapbytes(typecast(newdata,cid));
end
data=typecast(newdata,cid);
adv=double(len*count);

%%-------------------------------------------------------------------------

function object = parse_array(varargin) % JSON array is written in row-major order
global pos inStr isoct
    parse_char('[');
    object = cell(0, 1);
    dim=[];
    type='';
    count=-1;
    if(next_char == '$')
        type=inStr(pos+1);
        pos=pos+2;
    end
    if(next_char == '#')
        pos=pos+1;
        if(next_char=='[')
            dim=parse_array(varargin{:});
            count=prod(double(dim));
        else
            count=double(parse_number());
        end
    end
    if(~isempty(type))
        if(count>=0)
            [object adv]=parse_block(type,count,varargin{:});
            if(~isempty(dim))
                object=reshape(object,dim);
            end
            pos=pos+adv;
            return;
        else
            endpos=matching_bracket(inStr,pos);
            [cid,len]=elem_info(type);
            count=(endpos-pos)/len;
            [object adv]=parse_block(type,count,varargin{:});
            pos=pos+adv;
            parse_char(']');
            return;
        end
    end
    if next_char ~= ']'
         while 1
            val = parse_value(varargin{:});
            object{end+1} = val;
            if next_char == ']'
                break;
            end
            %parse_char(',');
         end
    end
    if(jsonopt('SimplifyCell',0,varargin{:})==1)
      try
        oldobj=object;
        object=cell2mat(object')';
        if(iscell(oldobj) && isstruct(object) && numel(object)>1 && jsonopt('SimplifyCellArray',1,varargin{:})==0)
            object=oldobj;
        elseif(size(object,1)>1 && ndims(object)==2)
            object=object';
        end
      catch
      end
    end
    if(count==-1)
        parse_char(']');
    end

%%-------------------------------------------------------------------------

function parse_char(c)
    global pos inStr len
    skip_whitespace;
    if pos > len || inStr(pos) ~= c
        error_pos(sprintf('Expected %c at position %%d', c));
    else
        pos = pos + 1;
        skip_whitespace;
    end

%%-------------------------------------------------------------------------

function c = next_char
    global pos inStr len
    skip_whitespace;
    if pos > len
        c = [];
    else
        c = inStr(pos);
    end

%%-------------------------------------------------------------------------

function skip_whitespace
    global pos inStr len
    while pos <= len && isspace(inStr(pos))
        pos = pos + 1;
    end

%%-------------------------------------------------------------------------
function str = parseStr(varargin)
    global pos inStr esc index_esc len_esc
 % len, ns = length(inStr), keyboard
    type=inStr(pos);
    if type ~= 'S' && type ~= 'C' && type ~= 'H'
        error_pos('String starting with S expected at position %d');
    else
        pos = pos + 1;
    end
    if(type == 'C')
        str=inStr(pos);
        pos=pos+1;
        return;
    end
    bytelen=double(parse_number());
    if(length(inStr)>=pos+bytelen-1)
        str=inStr(pos:pos+bytelen-1);
        pos=pos+bytelen;
    else
        error_pos('End of file while expecting end of inStr');
    end

%%-------------------------------------------------------------------------

function num = parse_number(varargin)
    global pos inStr len isoct fileendian systemendian
    id=strfind('iUIlLdD',inStr(pos));
    if(isempty(id))
        error_pos('expecting a number at position %d');
    end
    type={'int8','uint8','int16','int32','int64','single','double'};
    bytelen=[1,1,2,4,8,4,8];
    datastr=inStr(pos+1:pos+bytelen(id));
    if(isoct)
        newdata=int8(datastr);
    else
        newdata=uint8(datastr);
    end
    if(id<=5 && fileendian~=systemendian)
        newdata=swapbytes(typecast(newdata,type{id}));
    end
    num=typecast(newdata,type{id});
    pos = pos + bytelen(id)+1;

%%-------------------------------------------------------------------------

function val = parse_value(varargin)
    global pos inStr len
    true = 1; false = 0;

    switch(inStr(pos))
        case {'S','C','H'}
            val = parseStr(varargin{:});
            return;
        case '['
            val = parse_array(varargin{:});
            return;
        case '{'
            val = parse_object(varargin{:});
            if isstruct(val)
                if(~isempty(strmatch('x0x5F_ArrayType_',fieldnames(val), 'exact')))
                    val=jstruct2array(val);
                end
            elseif isempty(val)
                val = struct;
            end
            return;
        case {'i','U','I','l','L','d','D'}
            val = parse_number(varargin{:});
            return;
        case 'T'
            val = true;
            pos = pos + 1;
            return;
        case 'F'
            val = false;
            pos = pos + 1;
            return;
        case {'Z','N'}
            val = [];
            pos = pos + 1;
            return;
    end
    error_pos('Value expected at position %d');
%%-------------------------------------------------------------------------

function error_pos(msg)
    global pos inStr len
    poShow = max(min([pos-15 pos-1 pos pos+20],len),1);
    if poShow(3) == poShow(2)
        poShow(3:4) = poShow(2)+[0 -1];  % display nothing after
    end
    msg = [sprintf(msg, pos) ': ' ...
    inStr(poShow(1):poShow(2)) '<error>' inStr(poShow(3):poShow(4)) ];
    error( ['JSONparser:invalidFormat: ' msg] );

%%-------------------------------------------------------------------------

function str = valid_field(str)
global isoct
% From MATLAB doc: field names must begin with a letter, which may be
% followed by any combination of letters, digits, and underscores.
% Invalid characters will be converted to underscores, and the prefix
% "x0x[Hex code]_" will be added if the first character is not a letter.
    pos=regexp(str,'^[^A-Za-z]','once');
    if(~isempty(pos))
        if(~isoct)
            str=regexprep(str,'^([^A-Za-z])','x0x${sprintf(''%X'',unicode2native($1))}_','once');
        else
            str=sprintf('x0x%X_%s',char(str(1)),str(2:end));
        end
    end
    if(isempty(regexp(str,'[^0-9A-Za-z_]', 'once' ))) return;  end
    if(~isoct)
        str=regexprep(str,'([^0-9A-Za-z_])','_0x${sprintf(''%X'',unicode2native($1))}_');
    else
        pos=regexp(str,'[^0-9A-Za-z_]');
        if(isempty(pos)) return; end
        str0=str;
        pos0=[0 pos(:)' length(str)];
        str='';
        for i=1:length(pos)
            str=[str str0(pos0(i)+1:pos(i)-1) sprintf('_0x%X_',str0(pos(i)))];
        end
        if(pos(end)~=length(str))
            str=[str str0(pos0(end-1)+1:pos0(end))];
        end
    end
    %str(~isletter(str) & ~('0' <= str & str <= '9')) = '_';

%%-------------------------------------------------------------------------
function endpos = matching_quote(str,pos)
len=length(str);
while(pos<len)
    if(str(pos)=='"')
        if(~(pos>1 && str(pos-1)=='\'))
            endpos=pos;
            return;
        end
    end
    pos=pos+1;
end
error('unmatched quotation mark');
%%-------------------------------------------------------------------------
function [endpos e1l e1r maxlevel] = matching_bracket(str,pos)
global arraytoken
level=1;
maxlevel=level;
endpos=0;
bpos=arraytoken(arraytoken>=pos);
tokens=str(bpos);
len=length(tokens);
pos=1;
e1l=[];
e1r=[];
while(pos<=len)
    c=tokens(pos);
    if(c==']')
        level=level-1;
        if(isempty(e1r)) e1r=bpos(pos); end
        if(level==0)
            endpos=bpos(pos);
            return
        end
    end
    if(c=='[')
        if(isempty(e1l)) e1l=bpos(pos); end
        level=level+1;
        maxlevel=max(maxlevel,level);
    end
    if(c=='"')
        pos=matching_quote(tokens,pos+1);
    end
    pos=pos+1;
end
if(endpos==0)
    error('unmatched "]"');
end

author: Ayatans

repo: Machine-Learning-homework


text:
function timestamp ( )

%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [rmse, primal_objs, dual_objs] = run_mbsdca(Xtrain, Ytrain, Xtest, Ytest, lambda, opts)
% Mini-batch SDCA
% Inputs
%   Xtrain: input training data
%   Ytrain: output training data
%   Xtest: input test data
%   Ytest: output test data
%   lambda: regularization parameter
%   opts: optional arguments
% Output
%   Average RMSE across tasks, primal and dual objectives

%% intialize variables
fprintf('Running Mb-SDCA\n');
m = length(Xtrain); % # of tasks
d = size(Xtrain{1}, 2); % # of features
W = zeros(d, m); alpha = cell(m,0);
Sigma = eye(m) * (1/m);
Omega = inv(Sigma);
totaln = 0; n = zeros(m, 1);
for t = 1:m
    n(t) = length(Ytrain{t});
    totaln = totaln + n(t);
    alpha{t} = zeros(n(t), 1);
end

%% intialize counters
rmse = zeros(opts.mbsdca_inner_iters, 1);
dual_objs = zeros(opts.mbsdca_inner_iters, 1);
primal_objs = zeros(opts.mbsdca_inner_iters, 1);

for h = 1:opts.mbsdca_outer_iters
    % update W
    for hh = 1:opts.mbsdca_inner_iters
        rng(hh * 1000);
        if(opts.sys_het)
            sys_iters = (opts.top - opts.bottom) .* rand(m,1) + opts.bottom;
        end

        % compute RMSE
        rmse(hh) = compute_rmse(Xtest, Ytest, W, opts);
        primal_objs(hh) = compute_primal(Xtrain, Ytrain, W, Omega, lambda);
        dual_objs(hh) = compute_dual(alpha, Ytrain, W, Omega, lambda);

        % loop over tasks (in parallel)
        deltaB = zeros(d, m);
        alpha_prev = alpha;
        local_iters = cell(m, 0);
        for t = 1:m
            tperm = randperm(n(t));
            alpha_t = alpha{t};
            curr_sig = Sigma(t,t);
            if(opts.sys_het)
                local_iters{t} = n(t) * sys_iters(t);
            else
                local_iters{t} = n(t) * opts.mocha_sdca_frac;
            end

            % run SDCA locally
            for s=1:local_iters{t}
                % select random coordinate
                idx = tperm(mod(s, n(t)) + 1);
                alpha_old = alpha_t(idx);
                curr_y = Ytrain{t}(idx);
                curr_x = Xtrain{t}(idx, :);

                % compute update
                update = (curr_y * curr_x * W(:,t));
                grad = lambda * n(t) * (1.0 - update) / (curr_sig * (curr_x * curr_x')) + (alpha_old * curr_y);
                alpha_t(idx) = curr_y * max(0.0, min(1.0, grad));
                deltaB(:, t) = deltaB(:, t) + (alpha_t(idx) - alpha_old) * curr_x' / n(t);
                alpha{t} = alpha_t;
            end
        end

        scaling = opts.mbsdca_scaling / sum(cell2mat(local_iters(:)));
        % combine updates globally
        for t = 1:m
            alpha{t} = alpha_prev{t} + scaling * (alpha{t} - alpha_prev{t});
            for tt = 1:m
                W(:, t) = W(:, t) + scaling * deltaB(:, tt) * Sigma(t, tt) * (1.0 / lambda);
            end
        end

    end

    %% make sure eigenvalues are positive
    A = W'*W;
    if(any(eig(A) < 0))
        [V,Dmat] = eig(A);
        dm= diag(Dmat);
        dm(dm <= 1e-7) = 1e-7;
        D_c = diag(dm);
        A = V*D_c*V';
    end

    %% update Omega, Sigma
    sqm = sqrtm(A);
    Sigma = sqm / trace(sqm);
    Omega = inv(Sigma);

end

end

author: gingsmith

repo: fmtl

save_path: github-repos/MATLAB/gingsmith-fmtl


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

% used market data
load 'MODEL_Hist';
% load 'MODEL_Hist_new';

% Discretization parameters
NTime = 12;                    % number of time steps
NSim = 100;                  % number of simulations
NBatches = 1;                  % number of batches

NHist = 251;                     % number of historic scenarios

NModels = 10;                  % number of models

% Options
%K = 6800;                              % Strike
%C = 1;                                 % call (1) or put (0)
LF = 0; LC = 0.04; GF = 0; GC = 0.15;   % cliquet params
U = Spot * 0.9;
L = Spot * 1.1;
T= 1;                                   % maturity of option

% Comment in the option you wish to price
option1 = 'Arithmetic Asian Call';
option2 = 'Arithmetic Asian Put';
exotic_price_c = @(x,y,r) Price_ArithmeticAsian(x,y,1,r,T);
exotic_price_p = @(x,y,r) Price_ArithmeticAsian(x,y,0,r,T);

%option1 = 'FixedStrike Lookback Call';
%option2 = 'FixedStrike Lookback Put';
%exotic_price_c = @(x,y,r) Price_LookBackFixedStrike(x,y,1,r,T);
%exotic_price_p = @(x,y,r) Price_LookBackFixedStrike(x,y,0,r,T);

%option1 = 'FloatingStrike Lookback Call';
%option2 = 'FloatingStrike Lookback Put';
%exotic_price_c = @(x,y,r) Price_LookBackFloatingStrike(x,y,1,r,T);
%exotic_price_p = @(x,y,r) Price_LookBackFloatingStrike(x,y,0,r,T);

% option1 = 'Arithmetic Cliquet';
% exotic_price_c = @(x,r) Price_ArithmeticCliquet(x,LF,LC,GF,GC,r,T);

% option1 = 'Knock-Out Call';
% option2 = 'Knock-Out Put';
% exotic_price_c = @(x,y,r) Price_KnockOut(x,x(1,1),y,1,r,T);
% exotic_price_p = @(x,y,r) Price_KnockOut(x,x(1,1),y,0,r,T);

vanilla_price_c = @(x,y,r) Price_CallPut(x,y,1,r,T);
vanilla_price_p = @(x,y,r) Price_CallPut(x,y,0,r,T);

% output
ep1c = zeros(NHist,NModels); %ep2c = ep1c; ep3c = ep1c;
ep1p = ep1c; %ep2p = ep1c; ep3p = ep1c;

vp_c = ep1c; vp_p = ep1c;
qp_c = ep1c; qp_p = ep1c;

vc = zeros(NHist,1); vp = vc;

% Black Scholes
for l = 1: NHist
    pathS = MC_B(Spot(l),r1y(l),0,T,BS_Par_Hist(l),NTime,NSim,NBatches);
    ep1c(l,1) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,1) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,1) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,1) = vanilla_price_p(pathS,Spot(l),r1y(l));
    [vc(l), vp(l)] = blsprice(Spot(l), Spot(l), r1y(l), T, BS_Par_Hist(l),0);
    qp_c(l,1) = ep1c(l,1)/vp_c(l,1);
    qp_p(l,1) = ep1p(l,1)/vp_p(l,1);
end

% Merton
for l = 1: NHist
    pathS = MC_M(Spot(l),r1y(l),0,T,...
        MJ_Par_Hist(l,1),MJ_Par_Hist(l,2),MJ_Par_Hist(l,3),MJ_Par_Hist(l,4),...
        NTime,NSim,NBatches);
    ep1c(l,2) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,2) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,2) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,2) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,2) = ep1c(l,2)/vp_c(l,2);
    qp_p(l,2) = ep1p(l,2)/vp_p(l,2);
end

% Heston
for l = 1: NHist
    [pathS, pathV] = MC_QE(Spot(l),r1y(l),0,T,...
        HESTON_Par_Hist(l,1),HESTON_Par_Hist(l,2),HESTON_Par_Hist(l,4),...
        HESTON_Par_Hist(l,5),HESTON_Par_Hist(l,3),...
        NTime,NSim,NBatches);
    ep1c(l,3) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,3) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,3) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,3) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,3) = ep1c(l,3)/vp_c(l,3);
    qp_p(l,3) = ep1p(l,3)/vp_p(l,3);
end

% Bates
for l = 1: NHist
    [pathS, pathV] = MC_QE_j(Spot(l),r1y(l),0,T,...
        BATES_Par_Hist(l,1),BATES_Par_Hist(l,2),BATES_Par_Hist(l,4),BATES_Par_Hist(l,5),...
        BATES_Par_Hist(l,3),BATES_Par_Hist(l,7),BATES_Par_Hist(l,7),BATES_Par_Hist(l,8),...
        NTime,NSim,NBatches);
    ep1c(l,4) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,4) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,4) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,4) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,4) = ep1c(l,4)/vp_c(l,4);
    qp_p(l,4) = ep1p(l,4)/vp_p(l,4);
end

% Variance Gamma CGM
for l = 1:NHist
    pathS  = MC_VG_CGM(Spot(l),r1y(l),0,T,...
        VG_Par_Hist(l,1),VG_Par_Hist(l,2),VG_Par_Hist(l,3),...
        NTime,NSim,NBatches);
    ep1c(l,5) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,5) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,5) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,5) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,5) = ep1c(l,5)/vp_c(l,1);
    qp_p(l,5) = ep1p(l,5)/vp_p(l,1);
end

% Normal Inverse Gaussian
for l = 1:NHist
    pathS  = MC_NIG(Spot(l),r1y(l),0,T,...
        NIG_Par_Hist(l,1),NIG_Par_Hist(l,2),NIG_Par_Hist(l,3),...
        NTime,NSim,NBatches);
    ep1c(l,6) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,6) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,6) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,6) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,6) = ep1c(l,1)/vp_c(l,1);
    qp_p(l,6) = ep1p(l,6)/vp_p(l,1);
end

% Variance Gamma - GOU
for l = 1:NHist
    pathS  = MC_VGGOU(Spot(l),r1y(l),0,T,...
        VGGOU_Par_Hist(l,1), VGGOU_Par_Hist(l,2), VGGOU_Par_Hist(l,3),...
        VGGOU_Par_Hist(l,6), VGGOU_Par_Hist(l,4), VGGOU_Par_Hist(l,5),...
        NTime,NSim,NBatches);
    ep1c(l,7) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,7) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,1) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,7) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,7) = ep1c(l,7)/vp_c(l,1);
    qp_p(l,7) = ep1p(l,7)/vp_p(l,1);
end

% % Variance Gamma - CIR
for l = 1: NHist
    pathS =  MC_VGCIR(Spot(l),r1y(l),0,T,...
        VGCIR_Par_Hist(l,1),VGCIR_Par_Hist(l,2),VGCIR_Par_Hist(l,3), ...
        VGCIR_Par_Hist(l,4),VGCIR_Par_Hist(l,5),VGCIR_Par_Hist(l,6),...
        NTime,NSim,NBatches);
    ep1c(l,8) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,8) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,8) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,8) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,8) = ep1c(l,8)/vp_c(l,1);
    qp_p(l,8) = ep1p(l,8)/vp_p(l,1);
end

% NIG - GOU
for l = 1: NHist
    pathS =  MC_NIGGOU(Spot(l),r1y(l),0,T,...
        NIGGOU_Par_Hist(l,1),NIGGOU_Par_Hist(l,2),NIGGOU_Par_Hist(l,3), ...
        NIGGOU_Par_Hist(l,4),NIGGOU_Par_Hist(l,5),NIGGOU_Par_Hist(l,6),NTime,NSim,NBatches);
    ep1c(l,9) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,9) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,9) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,9) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,9) = ep1c(l,9)/vp_c(l,1);
    qp_p(l,9) = ep1p(l,9)/vp_p(l,1);
end

% NIG - CIR
for l = 1: NHist
    pathS =  MC_NIGCIR(Spot(l),r1y(l),0,T,...
        NIGCIR_Par_Hist(l,1),NIGCIR_Par_Hist(l,2),NIGCIR_Par_Hist(l,3), ...
        NIGCIR_Par_Hist(l,5),NIGCIR_Par_Hist(l,6),NIGCIR_Par_Hist(l,4),NTime,NSim,NBatches);
    ep1c(l,10) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,10) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,10) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,10) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,10) = ep1c(l,10)/vp_c(l,1);
    qp_p(l,10) = ep1p(l,10)/vp_p(l,1);
end
m1 = 'BS'; m2 = 'MJ'; m3 = 'Heston'; m4 = 'Bates';
m5 = 'VG'; m6 = 'NIG'; m7 = 'VGGOU'; m8 = 'VGCIR';
m9 = 'NIGGOU'; m10 =  'NIGCIR';

% plot the results
plotit(ep1c,option1,m1, m2, m3, m4, m5, m6 ,m7, m8, m9, m10);
plotit(ep1p,option2,m1, m2, m3, m4, m5, m6 ,m7, m8, m9, m10);
plotit(qp_c,strcat('Normalized ',option1),m1, m2, m3, m4, m5, m6 ,m7, m8, m9, m10);
plotit(qp_p,strcat('Normalized ',option2),m1, m2, m3, m4, m5, m6 ,m7, m8, m9, m10);

modelriskc = repmat(max(ep1c,[],2),1,NModels)-ep1c;
modelriskp = -repmat(max(ep1p,[],2),1,NModels)+ep1p;

author: Sable

repo: mcbench-benchmarks


text:
function s = mean(f, dim)
%MEAN   Average or mean value of a DISKFUN.
%   MEAN(F) takes the mean in the angular direction (default), i.e.,
%          MEAN(F) = sum(F).
%
%   MEAN(F, DIM) takes the mean along the direction DIM. If DIM = 1 it is the
%   radial direction, and if DIM = 2 then it is the angular direction,
%   i.e., MEAN(F,2) = 1/(2*pi)*sum(F,2).
%
% See also DISKFUN/MEAN2, DISKFUN/STD2.

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

% Empty check:
if ( isempty(f) )
    s = chebfun;
    return
end

if ( nargin == 1)
    % Default to the angular direction.
    dim = 1;
end

s = sum(f, dim);
if ( dim == 1 )
    return; % Mean in the angular direction (default)
elseif ( dim == 2 )
    s = s / (2*pi); % Mean in the radial direction
else
    error('CHEBFUN:DISKFUN:mean:dim', ...
        'dim must be 1 or 2.')
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function showCoinSeg(img, coinseg, edgeColor, lineStyle)
% function showCoinSeg(CoinSeg, coinseg, edgeColor, lineStyle)
%
% purpose:
%   Utility to display segmentation parameters found by functions like
%   SegScaleAccHT on an image. A circle is draw at the segmentation.
% parameters:
%  IN:
%    img       image
%    coinseg   segmentation triplet [radius, centerY, centerX]
%    edgeColor color of the coircle to draw (optional)
%              default 'r' (red)
%    edgeColor line style of the coircle to draw (optional)
%              default '--'
%
% Example: (read coin image, segment coin and show segmentation)
%              img = imread('img1.png');
%              segdata = showCoinSeg(img, segScaleAccHT(img));
%

%
% (C) Christian Kotz 2007
%
% $Id: showCoinSeg.m,v 1.1 2008/04/25 17:39:40 ck Exp $
% author:  $Author: ck $
% date:    $Date: 2008/04/25 17:39:40 $
% version: $Revision: 1.1 $
% history:
% $Log: showCoinSeg.m,v $
% Revision 1.1  2008/04/25 17:39:40  ck
% added example
%
% Revision 1.2  2008/04/22 12:30:50  ck
% _
%
% Revision 1.1  2008/04/22 12:29:25  ck
% _
%
%

if nargin < 3
    edgeColor = 'r';
end

if nargin < 4
    lineStyle = '--';
end

r = coinseg(1);
y = coinseg(2);
x = coinseg(3);

imshow(img);
hold on;

rectangle('Position', [x-abs(r),y-abs(r),abs(2*r),abs(2*r)], 'Curvature', [1,1], ...
    'EdgeColor', edgeColor, 'LineStyle', lineStyle);
hold off


text:
function [v,lambda] = trideigs ( a, b, TOL, Nmax, vects )

%QRST           determine all of the eigenvalues (and optionally all of the
%               eigenvectors) of a symmetric tridiagonal matrix using the
%               QR algorithm with Wilkinson shift
%
%     calling sequences:
%             [lambda, v] = qrst ( a, b, TOL, Nmax, vects )
%             [lambda, v] = qrst ( a, b, TOL, Nmax )
%             lambda = qrst ( a, b, TOL, Nmax )
%             qrst ( a, b, TOL, Nmax, vects )
%             qrst ( a, b, TOL, Nmax )
%
%     inputs:
%             a       vector containing elements along the main diagonal
%                     of the symmetric tridiagonal matrix whose eigenvalues
%                     are to be determined
%             b       vector containing elements along the off diagonal
%                     of the symmetric tridiagonal matrix whose eigenvalues
%                     are to be determined
%             TOL     convergence tolerance
%             Nmax    maximum number of iterations
%             vects   optional input argument
%                     matrix containing eigenvector information produced
%                     during the reduction of the original symmetric
%                     matrix to symmetric tridiagonal form
%                     - this input is needed only if computation of the
%                       eigenvectors is requested (by including the second
%                       output argument) and the original matrix was not in
%                       symmetric tridiagonal form
%
%     output:
%             lambda  vector containing the eigenvalues of the symmetric
%                     tridiagonal matrix determined by the vectors a and b
%             v       optional output argument
%                     matrix containing the eigenvectors of the symmetric
%                     tridiagonal matrix determined by the vectors a and b
%                     - the i-th column of this matrix is an eigenvector
%                       which corresponds to the i-th eigenvalue in the
%                       vector lambda
%                     - eigenvectors will be not computed if this second
%                       output argument is omitted
%
%     NOTE:
%             if the maximum number of iterations is exceeded, a message
%             to this effect will be displayed, along with the number of
%             eigenvalues which had been determined - these eigenvalues
%             will be returned in the last entries of the output vector
%             lambda
%

n = length(a);
if ( length(b) == n-1 ) b(2:n) = b(1:n-1); end;

c = zeros ( 1, n );
s = zeros ( 1, n );
shift = 0;
togo = n;

if ( nargout >= 2 )
   if ( nargin >= 5 )
      v = vects;
   else
      v = eye(n);
   end;
end;

for its = 1 : Nmax

    if ( togo == 1 )
	   lambda(1) = a(1) + shift;
	   disp ( its );
	   return;
	end;

    trace = a(togo-1) + a(togo);
	det   = a(togo-1)*a(togo) - b(togo)*b(togo);
	disc  = sqrt ( trace*trace - 4*det );
	mu1 = (1/2) * ( trace + disc );
	mu2 = (1/2) * ( trace - disc );
	if ( abs ( mu1 - a(togo) ) < abs ( mu2 - a(togo) ) )
	   s = mu1;
	else
	   s = mu2;
	end;

    shift = shift + s;
	for i = 1:togo
	    a(i) = a(i) - s;
	end;

	oldb = b(2);
    for i = 2:togo
        j = i-1;
	    r = sqrt ( a(j)^2 + oldb^2 );
	    c(i) = a(j) / r;
	    s(i) = oldb / r;
	    a(j) = r;
	    temp1 = c(i)*b(i) + s(i)*a(i);
	    temp2 = -s(i)*b(i) + c(i)*a(i);
	    b(i) = temp1;
	    a(i) = temp2;
	    if ( i ~= togo ) oldb = b(i+1); b(i+1) = c(i)*b(i+1); end;
    end;

    a(1) = c(2)*a(1) + s(2)*b(2);
    b(2) = s(2)*a(2);
    for i = 2:togo-1
        a(i) = s(i+1)*b(i+1) + c(i)*c(i+1)*a(i);
	    b(i+1) = s(i+1)*a(i+1);
    end;
    a(togo) = c(togo)*a(togo);

	if ( nargout >= 2 )
	   for i = 2 : togo
	       col1 = v(:,i-1) * c(i) + v(:,i) * s(i);
		   v(:,i) = -s(i) * v(:,i-1) + c(i) * v(:,i);
		   v(:,i-1) = col1;
	   end;
	end;

	if ( abs(b(togo)) < TOL )
	   lambda(togo) = a(togo) + shift;
	   disp([lambda(togo) its]);
	   togo = togo - 1;
	end;

end;

disp ( 'qrst error: Maximum number of iterations exceeded' );
disp ( sprintf ( '%d eigenvalues determined \n', n-togo ) );

repo: Aero-matlab


text:
settings.ptype = 'MACRO';
% settings.filename = 'CantileverBeam_Triangle_Linear_Fine';
% settings.filename = 'topopt_quad';
% settings.filename = 'GrippingNew';
% settings.ptype = 'MICRO';
% settings.filename = 'RVE_Square_Triangle';
% settings.filename = 'RVE_Square_Triangle_Fine';

settings.method = 'SIMPALL';
%settings.method = 'SIMP_P3';
%settings.method = 'SIMP_Adaptative';

settings.material = 'ISOTROPIC';
settings.initial_case = 'full';
% settings.initial_case = 'circle';
% settings.initial_case = 'horizontal';
% settings.initial_case = 'square';
% settings.initial_case = 'feasible';
% settings.initial_case = 'rand';

settings.cost = {'compliance';'perimeter'}; %'chomog_fraction';'compliance';'perimeter';'chomog_alphabeta';'nonadjoint_compliance';
% settings.weights = []; %all 1
settings.weights = [1 0.1]; %compl+lambda*perimeter
settings.constraint = {'volumeConstraint'};

settings.optimizer = 'SLERP';
%settings.optimizer = 'PROJECTED GRADIENT';settings.incrementFactor = 1;
%settings.optimizer = 'MMA';
%settings.optimizer = 'IPOPT';

settings.filterType = 'P1';
%settings.filter = 'PDE';

settings.TOL.rho_plus = 1;
settings.TOL.rho_minus = 0;
settings.TOL.E_plus = 1;
settings.TOL.E_minus = 1e-3;
settings.TOL.nu_plus = 1/3;
settings.TOL.nu_minus = 1/3;

settings.nsteps = 1;
settings.Vfrac_final = settings.target_parameters.Vfrac;
settings.optimality_final = settings.target_parameters.optimality_tol;
settings.constr_final = settings.target_parameters.constr_tol;
settings.Vfrac_initial = 1;

settings.optimality_initial = 1e-3;
settings.constr_initial = 1e-3;

settings.micro.alpha = [1 1 0]';
settings.micro.beta = [1 1 0]';

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function this = coregister(this)
% Affine coregistraton of given images to a stationary image (geometry only)
%
%   Y = MrSeries()
%   Y.coregister(inputs)
%
% This is a method of class MrSeries.
%
% IN
%       parameters.coregister ->
%           .nameStationaryImage            (e.g. 'mean')
%           .nameTransformedImage           (e.g. 'anatomy')
%           .nameEquallyTransformedImages   (e.g.'tissueProbabilityMap*')
% OUT
%
% EXAMPLE
%   coregister
%
%   See also MrSeries

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

nameStationaryImage = this.parameters.coregister.nameStationaryImage;
nameTransformedImage = this.parameters.coregister.nameTransformedImage;
nameEquallyTransformedImages = ...
    this.parameters.coregister.nameEquallyTransformedImages;

stationaryImage = this.find('MrImage', 'name', ['^' nameStationaryImage '*']);
stationaryImage = stationaryImage{1};
transformedImage = this.find('MrImage', 'name', ['^' nameTransformedImage '*']);
transformedImage = transformedImage{1};
equallyTransformedImages = this.find('MrImage', 'name',...
    ['^' nameEquallyTransformedImages '*']);

% TODO: if transformed image has more than one volume, make sure to apply
% to all volumes!

% set paths of all images correctly
this.init_processing_step('coregister', transformedImage, ...
    equallyTransformedImages);

[~, affineCoregistrationGeometry] = transformedImage.coregister_to(...
    stationaryImage, 'applyTransformation', 'geometry');
this.parameters.coregister.affineCoregistrationGeometry = ...
    affineCoregistrationGeometry;

% now apply geometry Update to all other listed images that shall be
% transformed
% inverse transformation used, since coregister_to gives mapping from
% stationary to transformed image
nImages = numel(equallyTransformedImages);
for iImage = 1:nImages
    equallyTransformedImages{iImage}.affineTransformation.apply_inverse_transformation(...
        affineCoregistrationGeometry);
    equallyTransformedImages{iImage}.save;
end

this.finish_processing_step('coregister', transformedImage, ...
    equallyTransformedImages);

author: translationalneuromodeling

repo: tapas


text:
function [data, opt] = RawPixelExtractor(im, tmpl, opt)
sz = opt.FeatureExtractor.tmplsize;

if isfield(opt, 'pixel_im')
    im = opt.pixel_im;
    % disp('lab space image saved.');
else
    if (ndims(im) == 3)
%         if (norm(im(:,:,1) - im(:,:,2)) > 1e-6)
            im = RGB2Lab(im);
%         else
%             im = rgb2gray(im);
%         end
    end
end
data.tmpl = tmpl;

features = zeros(prod(sz)*size(im, 3), size(tmpl, 1));

% pad
minW = min(round(tmpl(:, 1) - tmpl(:, 3) / 2)) - 1;
maxW = max(round(tmpl(:, 1) + tmpl(:, 3) / 2)) + 1;
minH = min(round(tmpl(:, 2) - tmpl(:, 4) / 2)) - 1;
maxH = max(round(tmpl(:, 2) + tmpl(:, 4) / 2)) + 1;
[h, w, c] = size(im);
if (minW < 1)
    im_new = zeros(h, w + abs(minW) + 1, c);
    im_new(:, abs(minW) + 2:end, :) = im;
    im = im_new;
    tmpl(:, 1) = tmpl(:, 1) + abs(minW) + 1;
end
if (maxW > w)
    im_new = zeros(h, size(im, 2) + maxW - w, c);
    im_new(:, 1:size(im, 2), :) = im;
    im = im_new;
end
if (minH < 1)
    im_new = zeros(h + abs(minH) + 1, size(im, 2), c);
    im_new(abs(minH) + 2:end, :, :) = im;
    im = im_new;
    tmpl(:, 2) = tmpl(:, 2) + abs(minH) + 1;
end
if (maxH > h)
    im_new = zeros(size(im, 1) + maxH - h, size(im, 2), c);
    im_new(1:size(im, 1), :, :) = im;
    im = im_new;
end

for i = 1:size(tmpl, 1)
    midW    = tmpl(i, 1);
    midH    = tmpl(i, 2);
    w       = tmpl(i, 3);
    h       = tmpl(i, 4);

    tempIm = im(round(midH-h/2) : round(midH+h/2),...
                round(midW-w/2) : round(midW+w/2), :);
%     tempIm = imresize(tempIm, sz);
    tempIm = mexResize(tempIm, sz, 'auto');
    features(:, i) = tempIm(:);
    if (norm(features(:, i)) > 1e-6)
        features(:, i) = features(:, i) / norm(features(:, i));
    end
end
% features = features - 0.5;

data.feat = features;

author: flyers

repo: drone-tracking


text:
function varargout = segymat(varargin)
% segymat : Garphical User Interface for SegyMAT
%
%   (C) 2001-2004 Thomas Mejer Hansen, tmh@gfy.ku.dk/thomas@cultpenguin.com
%

% SEGYMAT M-file for segymat.fig
%      SEGYMAT, by itself, creates a new SEGYMAT or raises the existing
%      singleton*.
%
%      H = SEGYMAT returns the handle to a new SEGYMAT or the handle to
%      the existing singleton*.
%
%      SEGYMAT('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in SEGYMAT.M with the given input arguments.
%
%      SEGYMAT('Property','Value',...) creates a new SEGYMAT or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before segymat_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to segymat_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help segymat

% Last Modified by GUIDE v2.5 05-Jan-2009 17:29:20

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @segymat_OpeningFcn, ...
                   'gui_OutputFcn',  @segymat_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);

% Next line is edited
% ORIGINAL : if nargin & isstr(varargin{1})
if nargin>1 & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end

if nargin==1 & isstr(varargin{1})
    [path,file,suffix]=fileparts(varargin{1});
    segyfile.pathname=path;
    if length(suffix)>1,
        segyfile.filename=[file,'.',suffix];
    else
        segyfile.filename=file;
    end
    % FIGURE OUT HOW TO GET THE hObject ??
    % fReadSegy_Callback(hObject, eventdata, handles,segyfile);
end

% End initialization code - DO NOT EDIT

% --- Executes just before segymat is made visible.
function segymat_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to segymat (see VARARGIN)

  SegymatVerbose(['GUI : OPENING GUI'],20)
% Choose default command line output for segymat
handles.output = hObject;

% Initialize PlotPref
handles.PlotPref.Show=0;

% Update handles structure
guidata(hObject, handles);

% Resize Figure
fMain_ResizeFcn(hObject, eventdata, handles)

%  read in segy file
segyfile.pathname='c:\MaTLAB6p5\work\';
segyfile.filename='test.sgy';
%segyfile.filename='peaksh.segy';
%fReadSegy_Callback(hObject, eventdata, handles,segyfile);

% disable edit menu
UpdateMenus(hObject, eventdata, handles)

% UIWAIT makes segymat wait for user response (see UIRESUME)
% uiwait(handles.fMain);

% --- Outputs from this function are returned to the command line.
function varargout = segymat_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes when fMain window is resized.
function fMain_ResizeFcn(hObject, eventdata, handles)
% hObject    handle to fMain (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if (isfield(handles,'fMain')==0)
    return
end

SegymatVerbose(['GUI : Resizing'],20)

data=guidata(hObject);

% Get Size of Main Window
[mainPos]=get(handles.fMain,'Position');

d=4; % Distance to edges
d_ax=.2; % distance between axes
if isfield(data,'PlotPref')==0,
    data.PlotPref.Show=1;
    guidata(hObject,data);
end
if isfield(data.PlotPref,'Show')==0,
    data.PlotPref.Show=1;
    guidata(hObject,data);
end

% fHeader Location
if data.PlotPref.Show==0,
    fHeaderW=d;
    set(handles.fHeader,'Visible','Off');
    set(handles.text7,'Visible','Off');
    set(handles.txtGain,'Visible','Off');
    set(handles.pbGainUp,'Visible','Off');
    set(handles.pbGainDown,'Visible','Off');
    set(handles.eGainMin,'Visible','Off');
    set(handles.eGainMax,'Visible','Off');
    set(handles.popTop,'Visible','Off');
    set(handles.popBot,'Visible','Off');
    set(handles.txtTop,'Visible','Off');
    set(handles.txtBot,'Visible','Off');
    set(handles.txtColormap,'Visible','Off');
    set(handles.popColormap,'Visible','Off');
    set(handles.txtStyle,'Visible','Off');
    set(handles.popStyle,'Visible','Off');
else
    fHeaderW=30; % WIDTH OF INFO PAGE
    x0=d;y0=d;
    w=fHeaderW-d;
    h=mainPos(4)-2*d;
    x1=x0+w;
    y1=y0+h;

    Wuic=8; % WIDTH OF UICONTROL
    Huic=1.4; %    HEIGHT OF UICONTROL

    set(handles.fHeader,'Position',[x0 y0 w h])
    d_frame=1;
    x0=x0+d_frame;
    wframe=w-2*d_frame;
    % TXT HEADER
    ytop=mainPos(4)-d;
%    set(handles.text7,'Position',[d+d ytop fHeaderW-3*d 2]);
    set(handles.text7,'Position',[x0 ytop wframe 2]);

    % GAIN
    ytop=ytop-d_ax-Huic;
    set(handles.txtGain,'Position',[x0,ytop,wframe,Huic])
    ytop=ytop-d_ax-Huic;
    w_sign=4;
    set(handles.eGainMin,'Position',[x0,ytop,wframe-w_sign,Huic])
    set(handles.pbGainDown,'Position',[x0+wframe-w_sign,ytop,w_sign,Huic])
    ytop=ytop-d_ax-Huic;
    set(handles.eGainMax,'Position',[x0,ytop,Wuic*1.5,Huic])
    set(handles.pbGainUp,'Position',[x0+Wuic*1.5,ytop,Wuic/2,Huic])
    set(handles.eGainMax,'Position',[x0,ytop,wframe-w_sign,Huic])
    set(handles.pbGainUp,'Position',[x0+wframe-w_sign,ytop,w_sign,Huic])

    % TopSub
    ytop=ytop-2*d_ax-Huic;
    set(handles.txtTop,'Position',[x0,ytop,wframe,Huic])
    ytop=ytop-d_ax-Huic;
    set(handles.popTop,'Position',[x0,ytop,wframe,Huic])

    % BaseSub
    ytop=ytop-2*d_ax-Huic;
    set(handles.txtBot,'Position',[x0,ytop,wframe,Huic])
    ytop=ytop-d_ax-Huic;
    set(handles.popBot,'Position',[x0,ytop,wframe,Huic])

    % Style
    ytop=ytop-2*d_ax-Huic;
    set(handles.txtStyle,'Position',[x0,ytop,wframe,Huic])
    ytop=ytop-d_ax-Huic;
    set(handles.popStyle ,'Position',[x0,ytop,wframe,Huic])

    % Colormap
    ytop=ytop-2*d_ax-Huic;
    set(handles.txtColormap,'Position',[x0,ytop,wframe,Huic])
    ytop=ytop-d_ax-Huic;
    set(handles.popColormap ,'Position',[x0,ytop,wframe,Huic])

    %

    % MAKE ALL VISIBLE
    set(handles.fHeader,'Visible','On');
    set(handles.text7,'Visible','On');
    set(handles.txtGain,'Visible','On');
    set(handles.pbGainUp,'Visible','On');
    set(handles.pbGainDown,'Visible','On');
    set(handles.eGainMin,'Visible','On');
    set(handles.eGainMax,'Visible','On');
    set(handles.popTop,'Visible','On');
    set(handles.popBot,'Visible','On');
    set(handles.txtTop,'Visible','On');
    set(handles.txtBot,'Visible','On');
    set(handles.txtColormap,'Visible','On');
    set(handles.popColormap,'Visible','On');
    set(handles.txtStyle,'Visible','On');
    set(handles.popStyle,'Visible','On');
end

% AXIS
axLeft=fHeaderW+2*d;
Hsmallax=4;

% axBot
axBotX=axLeft;
axBotY=d;
axBotW=mainPos(3)-axLeft-2*d;
axBotH=Hsmallax;
axB=get(handles.axBot,'Position');
set(handles.axBot,'Position',[axBotX axBotY axBotW axBotH])

% axTop
axTopX=axLeft;
axTopY=mainPos(4)-Hsmallax-d;
axTopW=axBotW;
axTopH=Hsmallax;
axT=get(handles.axTop,'Position');
set(handles.axTop,'Position',[axTopX axTopY axTopW axTopH])

% axMain
axMainX=axLeft;
axMainY=axBotY+Hsmallax+d_ax;
axMainW=axBotW;
axMainH=mainPos(4)-axMainY-Hsmallax-d-d_ax;
set(handles.axMain,'Position',[axMainX axMainY axMainW axMainH])

% --------------------------------------------------------------------
function UpdatePrefs(hObject, eventdata, handles)
  SegymatVerbose(['GUI : UpdatePrefs'],20)
  data=guidata(hObject);

if isfield(data,'SegyData')==0,
    set(handles.eGainMin,'Enable','off');
    set(handles.eGainMin,'String','');
    set(handles.eGainMax,'Enable','off');
    set(handles.eGainMax,'String','');
   return
else
    set(handles.eGainMin,'Enable','on');
    set(handles.eGainMax,'Enable','on');
end

if isfield(data.PlotPref,'caxis')==0
    data.PlotPref.caxis=[min(data.SegyData(:)) max(data.SegyData(:))];
end
set(handles.eGainMin,'String',data.PlotPref.caxis(1));
set(handles.eGainMax,'String',data.PlotPref.caxis(2));

if isfield(data.PlotPref,'TraceHeaders')==0,
    data.PlotPref.TraceHeaders=fieldnames(data.SegyTraceHeaders);
end
if isfield(data.PlotPref,'TopPlot')==0, data.PlotPref.TopPlot=7; end
if isfield(data.PlotPref,'BotPlot')==0, data.PlotPref.BotPlot=13; end
set(handles.popTop,'String',data.PlotPref.TraceHeaders);
set(handles.popTop,'value',data.PlotPref.TopPlot);
set(handles.popTop,'Enable','On');
set(handles.popBot,'String',data.PlotPref.TraceHeaders);
set(handles.popBot,'value',data.PlotPref.BotPlot);
set(handles.popBot,'Enable','On');
guidata(hObject,data);

% --------------------------------------------------------------------
function UpdatePlots(hObject, eventdata, handles)
  SegymatVerbose(['GUI : Update Plots'],20)
  %data=guidata(hObject);

  UpdateMainPlot(hObject, eventdata, handles)

  UpdateTopPlot(hObject, eventdata, handles)
  UpdateBotPlot(hObject, eventdata, handles)

axes(handles.axMain);%%%zoom on;

% --------------------------------------------------------------------
function UpdateMainPlot(hObject, eventdata, handles)
  SegymatVerbose(['GUI : Update Main Plot'],20)
  data=guidata(hObject);
  Style=get(handles.popStyle,'Value'); % GET PLOTTING STYLE
  axes(handles.axMain);

  if ((Style==1)|(Style==4)|(Style==5))
    imagesc(data.SegyTrace,data.SegyTime,data.SegyData);
    caxis(data.PlotPref.caxis)
    hold on
  end

  gain=max(abs(data.PlotPref.caxis));
  if ((Style==2)|(Style==4))
    wiggle(data.SegyTrace,data.SegyTime,data.SegyData,'wiggle',gain);
  end
  if ((Style==3)|(Style==5))
    wiggle(data.SegyTrace,data.SegyTime,data.SegyData,'VA',gain);
  end
  hold off
  set(gca,'YaxisLocation','Right');
  set(gca,'XtickLabel',[]);
  ylabel('TWT [ms]')

  %% RESET ZOOM SUCH THAT THIS STATE IS THE 'ZOOM OUT' STATE.
  zoom reset;

  %%%zoom on;

% --------------------------------------------------------------------
function UpdateTopPlot(hObject, eventdata, handles)
  SegymatVerbose(['GUI : Update TOP Plots'],60)
  data=guidata(hObject);
  axes(handles.axMain);axMain=axis;

  pop=get(handles.popTop,'String');
  ipop=get(handles.popTop,'value');
  hname=pop{ipop};

  axes(handles.axTop)
  if isfield(data,'dTopPlot')
		if ipop~=data.dTopPlot_ipop
			data.dTopPlot=bar([data.SegyTraceHeaders.(hname)]);
			data.dTopPlot_ipop=ipop;
			guidata(hObject,data)
		end
	else
	  data.dTopPlot=bar([data.SegyTraceHeaders.(hname)]);
		data.dTopPlot_ipop=ipop;
	  guidata(hObject,data)
	end
  ax=axis;axis([axMain(1) axMain(2) ax(3) ax(4)])
  set(gca,'XaxisLocation','Top');
  set(gca,'YaxisLocation','Right');
  xlabel(hname)
  axes(handles.axMain);%%%zoom on;

function UpdateBotPlot(hObject, eventdata, handles)
  SegymatVerbose(['GUI : Update BOT Plots'],60)
  data=guidata(hObject);
  axes(handles.axMain);axMain=axis;

  pop=get(handles.popBot,'String');
  ipop=get(handles.popBot,'value');
  hname=pop{ipop};

  axes(handles.axBot)
	if isfield(data,'dBotPlot')
		if ipop~=data.dBotPlot_ipop
			data.dBotPlot=plot([data.SegyTraceHeaders.(hname)]);
			data.dBotPlot_ipop=ipop;
			guidata(hObject,data)
		end
	else
		data.dBotPlot=plot([data.SegyTraceHeaders.(hname)]);
	  data.dBotPlot_ipop=ipop;
		guidata(hObject,data)
	end

  % bar([data.SegyTraceHeaders.(hname)]);
  ax=axis;axis([axMain(1) axMain(2) ax(3) ax(4)])
  set(gca,'XaxisLocation','Bot');
  set(gca,'YaxisLocation','Right');
  xlabel(hname)
  axes(handles.axMain);%%%zoom on;

% --------------------------------------------------------------------
function UpdateGain(hObject, eventdata, handles)
  SegymatVerbose(['GUI : UpdateGain'],20)
  data=guidata(hObject);
  set(handles.eGainMin,'String',data.PlotPref.caxis(1));
  set(handles.eGainMax,'String',data.PlotPref.caxis(2));
  axes(handles.axMain);
  caxis(data.PlotPref.caxis);
  drawnow;

% --------------------------------------------------------------------
function KeyPressFcn_Callback(hObject, eventdata, handles)
data=guidata(hObject);
Key=get(gcf,'CurrentCharacter');
%SegymatVerbose(['GUI : KeyPressFcn :  ',num2str(double(Key)),' ',char(Key)],20)
%disp(char(Key));
%disp(double(Key));

%%%%%%%%%%%
% GAIN
%if (Key=='+'),
if ((double(Key)==30)|(double(Key)==43)|(double(Key)==115))
    data.PlotPref.caxis=data.PlotPref.caxis./(1.2);
    guidata(hObject,data);
    UpdateGain(hObject, eventdata, handles) ;
end
%if (Key=='-');
if ((double(Key)==31)|(double(Key)==45)|(double(Key)==120))
    data.PlotPref.caxis=data.PlotPref.caxis.*(1.2);
    guidata(hObject,data);
    UpdateGain(hObject, eventdata, handles) ;
end

%% SHOW PREFS
if (lower(Key)=='h')
    if isfield(data.PlotPref,'Show')==0, data.PlotPref.Show=1;end
    data.PlotPref.Show=1-data.PlotPref.Show;
    guidata(hObject,data);
    fMain_ResizeFcn(hObject,[],handles);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ZOOM ON MAIN AX
if isfield(data,'zoomfac')==0
    data.zoomfac=1.5;
    guidata(hObject,data);
end
zoom_in_keys=[29,97]; % ARROW RIGHT OR 'a'
zoom_out_keys=[28,122]; % ARROW LEFT OR 'z'
if isempty(Key)==1
  Key=0;
end

try
if find(double(Key)==zoom_in_keys);
	zoom(data.zoomfac)
	UpdateTopPlot(hObject, eventdata, handles)
	UpdateBotPlot(hObject, eventdata, handles)
end
catch
SegymatVerbose(sprintf('%s : failed to zoom',mfilename))
%keyboard;
end
if find(double(Key)==zoom_out_keys);
    zoom(-data.zoomfac)
	UpdateTopPlot(hObject, eventdata, handles)
	UpdateBotPlot(hObject, eventdata, handles)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MOVE
if isfield(data,'movefac')==0
    % Movefac=1; move in stepos of size of window
    % Movefac=0.5; move in steps of half-size of window
    data.movefac=.02;
    guidata(hObject,data);
end
if (double(Key)>=49)&(double(Key)<=117), % NUM PAD

  go_left_keys=[52,117]; % 2,j
	go_right_keys=[54,111]; % 2,j
	go_down_keys=[50,107]; % 2,j
	go_up_keys=[56]; % 2,j

  xlim=get(handles.axMain,'Xlim');
  ylim=get(handles.axMain,'Ylim');
  wx=(xlim(2)-xlim(1)).*data.movefac;
  wy=(ylim(2)-ylim(1)).*data.movefac;

  if (find(go_up_keys==double(Key))), set(handles.axMain,'Ylim',ylim-wy); end % 8
  if (find(go_down_keys==double(Key))), set(handles.axMain,'Ylim',ylim+wy); end % 2
  if (find(go_left_keys==double(Key))), set(handles.axMain,'Xlim',xlim-wx); end % 4
  if (find(go_right_keys==double(Key))), set(handles.axMain,'Xlim',xlim+wx); end % 6

  if (double(Key)==49), % 1
    set(handles.axMain,'Ylim',ylim+wy);
    set(handles.axMain,'Xlim',xlim-wx);
  end % 1
  if (double(Key)==51),  % 3
    set(handles.axMain,'Ylim',ylim+wy);
    set(handles.axMain,'Xlim',xlim+wx);
  end % 3
  if (double(Key)==55),  % 7
    set(handles.axMain,'Ylim',ylim-wy);
    set(handles.axMain,'Xlim',xlim-wx);
  end % 7
  if (double(Key)==57), % 9
    set(handles.axMain,'Ylim',ylim-wy);
    set(handles.axMain,'Xlim',xlim+wx);
  end % 9
  if (double(Key)==53),  % 5
    axes(handles.axMain)
    zoom out;

  end % 5
  xlim_after=get(handles.axMain,'Xlim');
  if (sum(xlim==xlim_after)<2)
    UpdateTopPlot(hObject, eventdata, handles)
    UpdateBotPlot(hObject, eventdata, handles)
  end
end

%
% IO FUNCTIONS
%

% --------------------------------------------------------------------
function mFileOpen_Callback(hObject, eventdata, handles)
% hObject    handle to mFileOpen (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
try
    data=guidata(hObject);
    if isfield(data,'SegyFile');
        segyfile=DIAOpenSEGY(data.SegyFile);
    else
        segyfile=DIAOpenSEGY;
    end
	fReadSegy_Callback(hObject, eventdata, handles,segyfile);
catch
end
% --------------------------------------------------------------------
function mOpenFileFast_Callback(hObject, eventdata, handles)
% hObject    handle to mOpenFileFast (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
try
	[segyfile.filename,segyfile.pathname]=uigetfile( ...
    {'*.segy;*.SEGY;*.SEGY;*.sgy','All Segy files'; ...
     '*.su;*.SU;*.sU;*.Su','All SU files'; ...
     '*','All Files'},...
    'Pick A SEGY file');
	fReadSegy_Callback(hObject, eventdata, handles,segyfile);
catch
end

% --------------------------------------------------------------------
function fReadSegy_Callback(hObject, eventdata, handles,segyfile);
data=guidata(hObject);
try
	[dpath,dfile,dsuffix]=fileparts(segyfile.filename);
	if (strcmp(lower(dsuffix),'.su'))
		[data.SegyData,data.SegyTraceHeaders,data.SegyHeader]=ReadSu(fullfile(segyfile.pathname,segyfile.filename));
    else
        if isfield(segyfile,'JUMP');
            k=0;
            try
            if segyfile.JUMP.enable==1
                k=k+1;v{k}='jump';
                k=k+1;v{k}=segyfile.JUMP.jump;
            end
            if segyfile.TIMERANGE.enable==1
                k=k+1;v{k}='trange';
                k=k+1;v{k}=segyfile.TIMERANGE.min;
                k=k+1;v{k}=segyfile.TIMERANGE.max;
            end
            if segyfile.TRACEHEADER.enable==1
                k=k+1;v{k}='minmax';
                k=k+1;v{k}=segyfile.TRACEHEADER.name;
                k=k+1;v{k}=segyfile.TRACEHEADER.min;
                k=k+1;v{k}=segyfile.TRACEHEADER.max;
            end
            catch
                SegymatVerbose(sprintf('%s : failed to make use of headers for reading',mfilename))
                %keyboard
            end
            if k>0
                [data.SegyData,data.SegyTraceHeaders,data.SegyHeader]=ReadSegy(fullfile(segyfile.pathname,segyfile.filename),v);
            else
                [data.SegyData,data.SegyTraceHeaders,data.SegyHeader]=ReadSegy(fullfile(segyfile.pathname,segyfile.filename));
            end
        else
            [data.SegyData,data.SegyTraceHeaders,data.SegyHeader]=ReadSegy(fullfile(segyfile.pathname,segyfile.filename));
        end
    end
    if isempty(data.SegyData)
          f = warndlg('No data was read !', 'ReadSegy warning', 'modal');
    else
        set(handles.fMain,'name',['SegyMAT : ',segyfile.filename])
        data.SegyFile=segyfile;
        data.SegyTime=[1:1:data.SegyHeader.ns].*data.SegyHeader.dt/1e+6;
        data.SegyTrace=[1:length(data.SegyTraceHeaders)];
        data.PlotPref.Show=1;
        data.PlotPref.caxis=[min(data.SegyData(:)) max(data.SegyData(:))];
        guidata(hObject,data);
        fMain_ResizeFcn(hObject, eventdata, handles);
        UpdatePrefs(hObject, eventdata, handles);
        UpdatePlots(hObject, eventdata, handles);
    end
catch
    errordlg('An error occured while reading the SEGY file','Error reading SGY file','modal')
    %keyboard
end
UpdateMenus(hObject, eventdata, handles);

% --------------------------------------------------------------------
function mFileSave_Callback(hObject, eventdata, handles)
% hObject    handle to mFileSave (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
try
	data=guidata(hObject);
	file=fullfile(data.SegyFile.pathname,data.SegyFile.filename);
	ButtonName=questdlg(['Are you sure you want to override :',file], ...
                       'Warning !!!', ...
                       'Yes','No','No');
  if strcmp(ButtonName,'Yes')
		[dpath,dfile,dsuffix]=fileparts(file);
		if (strcmp(lower(dsuffix),'.su'))
			WriteSuStructure(file,data.SegyHeader,data.SegyTraceHeaders,data.SegyData);
		else
			WriteSegyStructure(file,data.SegyHeader,data.SegyTraceHeaders,data.SegyData);
		end
	end
catch
	errordlg('An error occured while writing SEGY file','Error writing SGY file','modal')
end

% --------------------------------------------------------------------
function mFileSaveAs_Callback(hObject, eventdata, handles)
% hObject    handle to mFileSaveAs (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
try
	data=guidata(hObject);
	[segyfile.filename,segyfile.pathname]=uiputfile( ...
    {'*.segy;*.SEGY;*.SEGY;*.sgy','All Segy files'; ...
     '*.su;*.SU;*.sU;*.Su','All SU files'; ...
     '*','All Files'},...
    'Save as ');
	file=fullfile(segyfile.pathname,segyfile.filename);
	[dpath,dfile,dsuffix]=fileparts(file);
	if (strcmp(lower(dsuffix),'.su'))
		WriteSuStructure(file,data.SegyHeader,data.SegyTraceHeaders,data.SegyData);
	else
		WriteSegyStructure(file,data.SegyHeader,data.SegyTraceHeaders,data.SegyData);
	end
catch
	errordlg('An error occured while writing SEGY file','Error writing SGY file','modal')
end

% --------------------------------------------------------------------
function mHelp_Callback(hObject, eventdata, handles)
% hObject    handle to mHelp (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function mHelpHelp_Callback(hObject, eventdata, handles)
% hObject    handle to mHelpHelp (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
SegymatHelp('index')

% --------------------------------------------------------------------
function mHelpAbout_Callback(hObject, eventdata, handles)
% hObject    handle to mHelpAbout (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
DIAAboutSegymat;

% --- Executes on button press in pbGainUp.
function pbGainUp_Callback(hObject, eventdata, handles)
% hObject    handle to pbGainUp (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
  data=guidata(hObject);
  data.PlotPref.caxis=data.PlotPref.caxis./(1.2);
  guidata(hObject,data);
  UpdateGain(hObject, eventdata, handles) ;

	Style=get(handles.popStyle,'Value'); % GET PLOTTING STYLE
  if ((Style==2)|(Style==4|(Style==3)|(Style==5)))
		UpdateMainPlot(hObject, eventdata, handles)
	end

% --- Executes on button press in pbGainDown.
function pbGainDown_Callback(hObject, eventdata, handles)
% hObject    handle to pbGainDown (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
  data=guidata(hObject);
  data.PlotPref.caxis=data.PlotPref.caxis.*(1.2);
  guidata(hObject,data);
  UpdateGain(hObject, eventdata, handles) ;

	Style=get(handles.popStyle,'Value'); % GET PLOTTING STYLE
  if ((Style==2)|(Style==4|(Style==3)|(Style==5)))
		UpdateMainPlot(hObject, eventdata, handles)
	end

% --- Executes during object creation, after setting all properties.
function eGainMin_CreateFcn(hObject, eventdata, handles)
% hObject    handle to eGainMin (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function eGainMin_Callback(hObject, eventdata, handles)
% hObject    handle to eGainMin (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of eGainMin as text
%        str2double(get(hObject,'String')) returns contents of eGainMin as a double
data=guidata(hObject);
data.PlotPref.caxis(1)=str2double(get(hObject,'String'));
guidata(hObject,data);
UpdateGain(hObject, eventdata, handles);

% --- Executes during object creation, after setting all properties.
function eGainMax_CreateFcn(hObject, eventdata, handles)
% hObject    handle to eGainMax (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function eGainMax_Callback(hObject, eventdata, handles)
% hObject    handle to eGainMax (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of eGainMax as text
%        str2double(get(hObject,'String')) returns contents of eGainMax as a double
data=guidata(hObject);
data.PlotPref.caxis(2)=str2double(get(hObject,'String'));
guidata(hObject,data);
UpdateGain(hObject, eventdata, handles);

% --- Executes during object creation, after setting all properties.
function popTop_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popTop (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in popTop.
function popTop_Callback(hObject, eventdata, handles)
% hObject    handle to popTop (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns popTop contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popTop
  SegymatVerbose(['GUI : Update popTOP'],20)
  UpdateTopPlot(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function popBot_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popBot (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in popBot.
function popBot_Callback(hObject, eventdata, handles)
% hObject    handle to popBot (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns popBot contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popBot
  SegymatVerbose(['GUI : Update popBOT'],20)
  UpdateBotPlot(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function popColormap_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popColormap (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in popColormap.
function popColormap_Callback(hObject, eventdata, handles)
% hObject    handle to popColormap (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns popColormap contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popColormap
  SegymatVerbose(['GUI : Update popColormap'],20)
  data=guidata(hObject);
  cmap=get(hObject,'String');
  icmap=get(hObject,'value');
  axes(handles.axMain);
  colormap(cmap{icmap});
  %%% zoom on;

% --- Executes during object creation, after setting all properties.
function popStyle_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popStyle (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in popStyle.
function popStyle_Callback(hObject, eventdata, handles)
% hObject    handle to popStyle (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns popStyle contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popStyle
  SegymatVerbose(['GUI : popStyle'],20)
  UpdateMainPlot(hObject, eventdata, handles);
  UpdateTopPlot(hObject, eventdata, handles)
  UpdateBotPlot(hObject, eventdata, handles)

%
% MENUS
%
% --------------------------------------------------------------------
function UpdateMenus(hObject, eventdata, handles)
  SegymatVerbose(['GUI : Update Menus'],20)
	data=guidata(hObject);

	if isfield(data,'SegyData')==0
		set(handles.mEdit,'Visible','Off')
	else
		set(handles.mEdit,'Visible','On')
	end

% --------------------------------------------------------------------
function mFile_Callback(hObject, eventdata, handles)
% hObject    handle to mFile (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function mEdit_Callback(hObject, eventdata, handles)
% hObject    handle to mEdit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function mEditSH_Callback(hObject, eventdata, handles)
% hObject    handle to mEditSH (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
data=guidata(hObject);
try
	data.SegyHeader=GUIEditSegyHeader(data.SegyHeader);
	guidata(hObject,data)
catch
end
% --------------------------------------------------------------------
function mEditSTH_Callback(hObject, eventdata, handles)
% hObject    handle to mEditSTH (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
data=guidata(hObject);
try
	data.SegyTraceHeaders=GUIEditSegyTraceHeader(data.SegyTraceHeaders);
	guidata(hObject,data)
	UpdateTopPlot(hObject, eventdata, handles)
	UpdateBotPlot(hObject, eventdata, handles)
catch
  SegymatVerbose('Something Went wrong calling GUISegyMat')
end

% --------------------------------------------------------------------
function mPlotXY_Callback(hObject, eventdata, handles)
% hObject    handle to mPlotXY (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
data=guidata(hObject);
try
    GUIPlotXY(data.SegyTraceHeaders);
catch
  SegymatVerbose('Something Went wrong calling GUISegyMat')
end

% --------------------------------------------------------------------
function mplot_Callback(hObject, eventdata, handles)
% hObject    handle to mplot (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function mEditTextualHeader_Callback(hObject, eventdata, handles)
% hObject    handle to mEditTextualHeader (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
data=guidata(hObject);
try
	data.SegyHeader=GUIEditTextualFileHeader(data.SegyHeader);
	guidata(hObject,data)
catch
end

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
% make sure that running JAABADetect one classifier at a time gives the
% same results as running them in batch

function maxerr = FlyBowlJAABADetect_SanityCheck(expdir,classifierparamsfiles,varargin)

[classifierparams,isrelativepath] = ...
  myparse(varargin,'classifierparams',[],'isrelativepath',true);

if isempty(classifierparams),
  classifierparams = ReadClassifierParamsFile(classifierparamsfiles,'isrelativepath',isrelativepath);
end
nbehaviors = numel(classifierparams);

allScores = cell(1,nbehaviors);
maxerr = zeros(1,nbehaviors);
for i = 1:nbehaviors,

  behavior = classifierparams(i).behaviors.names{1};

  scorefile = fullfile(expdir,classifierparams(i).file.scorefilename);
  if ~exist(scorefile,'file'),
    warning('Score file %s does not exist',scorefile);
    continue;
  end

  allScores(i) = JAABADetect(expdir,'classifierfiles',{classifierparams(i).classifierfile},...
    'configfiles',{classifierparams(i).configfile},'debug',true,'forcecompute',true);

  scoresload = load(scorefile);
  nflies = numel(allScores{i});
  nflies2 = numel(scoresload.allScores.scores);
  if nflies ~= nflies2,
    warning('Number of flies does not match for %s',behavior);
    maxerr(i) = inf;
  end
  nflies = min(nflies,nflies2);
  for fly = 1:nflies,
    if numel(scoresload.allScores.scores{fly}) ~= numel(allScores{i}{fly}),
      warning('Number of frames does not match for %s fly %d\n',behavior,fly);
      maxerr(i) = inf;
    end
    realidx = ~isnan(scoresload.allScores.scores{fly});
    if ~all(~realidx == isnan(allScores{i}{fly})),
      fprintf('%s fly %d: nans do not match\n',behavior,fly);
      maxerr(i) = inf;
      continue;
    end
    maxdiff = max(abs(scoresload.allScores.scores{fly}(realidx) - allScores{i}{fly}(realidx)));
    if maxdiff > 0,
      fprintf('%s fly %d: scores do not match, maxdiff = %f\n',behavior,fly,maxdiff);
      maxerr(i) = max(maxerr(i),maxdiff);
    end
  end

end

author: kristinbranson

repo: JAABA


text:
close all
clear

robot_num = 5;
% reference
angle_central = 2*pi/robot_num;
if ~mod(robot_num,2)
    pos_ref = [cos([0 kron(1:(robot_num-1)/2,[1,-1]) robot_num/2]*angle_central)',...
        sin([0 kron(1:(robot_num-1)/2,[1,-1]) robot_num/2]*angle_central)']/2;
else
    pos_ref = [cos([0 kron(1:(robot_num-1)/2,[1,-1])]*angle_central)',...
        sin([0 kron(1:(robot_num-1)/2,[1,-1])]*angle_central)']/2;
end
% topology
graph = bigraph(robot_num,pos_ref);
% circle centers and initials
pos_circ = pos_ref;
% pos_rob = pos_circ.*rand(robot_num,2)*2;
pos_rob = pos_ref*rot2(2*pi/3);
color_list = ['b','g','m','r','c'];
for i=1:robot_num
    plot(pos_rob(i,1),pos_rob(i,2),[color_list(i) 'o']); hold on
    circle_(pos_circ(i,:),1,'color',color_list(i));
    plot(pos_ref(i,1),pos_ref(i,2),[color_list(i) 'd']);
end
hold off
% initial control
uqp_rob = zeros(size(pos_rob));
% reference distance
dis_ref = distances(graph,pos_ref);
ld = dis_ref-0.2; ud = dis_ref+0.1;

% circular constraint sets
zbf_circ = @(i,pos_rob) 1-(pos_rob(i,:)-pos_circ(i,:))*(pos_rob(i,:)-pos_circ(i,:))';
% nominal control
knom = 2;
unom_rob = @(i,pos_rob) -knom*(pos_rob(i,:)-pos_ref(i,:));
% extended K-class function
akcf = 1;
ekcf = @(h) akcf*h;
% control barrier function (less or equal)
cbfa_circ = @(i,pos_rob) pos_rob(i,:)-pos_circ(i,:);
cbfb_circ = @(i,pos_rob) ekcf(zbf_circ(i,pos_rob));

pos_data = pos_rob;
dt = 0.01;
T = 10;
loop = 0;
for t=0:dt:T
    loop = loop+1;
    % control barrier function - quadratic program
    for i=1:robot_num
        fqp = -unom_rob(i,pos_rob)';
        aqp = cbfa_circ(i,pos_rob);
        bqp = cbfb_circ(i,pos_rob);
        [aqp,bqp] = bdzcbf(graph,i,pos_rob,ld,ud,aqp,bqp,uqp_rob);
        [uqp,~,is_solved] = quadprog(eye(2),fqp,aqp,bqp);
        if is_solved>0
            uqp_rob(i,:) = uqp';
        else
            uqp_rob(i,:) = [0,0];
        end
    end
    % update
    pos_rob = pos_rob+dt*uqp_rob;
    pos_data(:,:,loop) = pos_rob;
    % plot
    for i=1:robot_num
        plot(pos_rob(i,1),pos_rob(i,2),[color_list(i) 'o']); hold on
        circle_(pos_circ(i,:),1,'color',color_list(i));
        plot(pos_ref(i,1),pos_ref(i,2),[color_list(i) 'd']);
        quiver(pos_rob(i,1),pos_rob(i,2),uqp_rob(i,1)/knom,uqp_rob(i,2)/knom,'color',color_list(i));
    end
    hold off
    drawnow
end

repo: paper-simulation


text:
% fastregress - perform fast regression and return p-value
%
% Usage:
% [ypred, alpha, rsq, B] = myregress(x, y, plotflag);
%
% Inputs
%  y - y values
%  x - x values
%  plotflag - [0|1] plot regression
%
% Outputs
%  ypred - y prediction
%  alpha - significance level
%  R^2   - r square
%  slope - slope of the fit
%
% Arnaud Delorme, 25 Feb 2003

function [ypred, alpha, rsq, B, BINT] = fastregress(x, y, plotflag)

    if nargin < 1
        help fastregress; return;
    end;

    % this part is useless but still works
    %B=polyfit(x, y, 1);         % B is the slope
    %ypred = polyval(B,x);       % predictions
    %dev = y - mean(y);          % deviations - measure of spread
    %SST = sum(dev.^2);          % total variation to be accounted for
    %resid = y - ypred;          % residuals - measure of mismatch
    %SSE = sum(resid.^2);        % variation NOT accounted for
    %rsq = 1 - SSE/SST;          % percent of error explained
    % see the page  http://www.facstaff.bucknell.edu/maneval/help211/fitting.html

    [B,BINT,R,RINT,STATS] = regress(y(:), [ones(length(x),1) x(:)]);
    alpha = STATS(3);
    rsq   = STATS(1);

    %note that we also have
    %ypred =  [ones(size(x,2),1) x(:)]*B;
    ypred =  x*B(2) + B(1);

    % B(1) contain the offset, B(2) the slope
    B = B(2);

    if nargin > 2
        hold on;
        [ynew tmp] = sort(ypred);
        xnew       = x(tmp);
        plot(xnew, ynew, 'r');
        legend(sprintf('R^2=%f', rsq), sprintf('p  =%f', alpha));
    end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function varargout = diff(varargin)
%DIFF   Derivative of a CHEBFUN2 object.
%   DIFF(F) is the derivative of F along the y direction.
%
%   DIFF(F, N) is the Nth derivative of F in the y direction.
%
%   DIFF(F, N, DIM) is the Nth derivative of F along the dimension DIM.
%     DIM = 1 (default) is the derivative in the y-direction.
%     DIM = 2 is the derivative in the x-direction.
%
%   DIFF(F, [NX NY]) is the partial derivative of NX of F in the first variable,
%   and NY of F in the second derivative. For example, DIFF(F,[1 2]) is
%   d^3F/dxd^2y.
%
% See also GRADIENT, SUM, PROD.

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

[varargout{1:nargout}] = diff@separableApprox(varargin{:});

end


text:
%kktextu 'Compute "texture" extracted from the surface surroundings '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros ktextu.pane file
%
% Parameters:
% InputFile: i1 '3D Data', required: 'Input scene (gray or binary scene)'
% InputFile: i2 'Coordinates', required: 'Coordinates at z-buffer distance (surface)'
% OutputFile: o '"Texture"', required: '"Texture" information which can be incorporated on shading'
% Double: step 'increment step', default: 1: 'Increment step during "texture" computation'
% Toggle: ave_type 'average value', default: 0: 'Meaning of "texture" is the average of voxel values'
% Toggle: wei_type 'weighted average value', default: 0: 'Meaning of "texture" is the weighted average of voxel values'
% Toggle: max_type 'maximum value', default: 0: 'Meaning of "texture" is the maximum voxel value'
% Toggle: min_type 'minimum value', default: 0: 'Meaning of "texture" is the minimum voxel value'
% Toggle: thi_type 'thickness (specify threshold)', default: 0: 'Meaning of "texture" is the length traced while below threshold'
% Toggle: norm_dir 'normal direction', default: 0: 'Follows the normal direction to extract "texture" information'
% InputFile: i3 'Normal Vectors', optional: 'Surface normal vectors'
% Toggle: view_dir 'view direction', default: 0: 'follows the viewer (projection) direction to extract "texture" information'
% Double: alpha 'alpha', default: 0: 'View plane rotation angle around Z axis'
% Double: beta 'beta ', default: 0: 'View plane rotation angle around X axis'
% Double: out_range 'outside surface', default: 0: 'Length of the region outside the surface where to compute "texture"'
% Double: in_range 'inside surface ', default: 10: 'Length of the region inside the surface where to compute "texture"'
% Double: thres_range 'threshold >=', default: 0: 'Threshold value which determines the region inside the surface where to compute "texture"'
%
% Example: o = kktextu({i1, i2, i3}, {'i1','';'i2','';'o','';'step',1;'ave_type',0;'wei_type',0;'max_type',0;'min_type',0;'thi_type',0;'norm_dir',0;'i3','';'view_dir',0;'alpha',0;'beta',0;'out_range',0;'in_range',10;'thres_range',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% ktextu - Compute texture information for each point lying on the visible surface
%
%  DESCRIPTION
%
% creates an image file based on volumetric information of a 3D object. This
% information can be of several types and is called here TEXTURE. It can can be
% obtained specifying the walk direction, the range of walking around the object
% surface and the type of texture. The output image is of type KFLOAT.
%
% This operator requires as inputs the volume (3D object) file, and the
% coordinate file, which is an object w x h x 1 x 1 x 3
% that can be generated optionally by the operator kzbuff, containing the
% coordinates X, Y and Z of each point of the object surface.
%
% The "ktextu" simulates that, for each point of the object's surface,
% there is a
% ray going inside the object in the specified direction. Then it walks (visits
% voxels) in that ray in the given range and gets the information according
% to the texture type.
%
% There are two options for the direction of "texture" extraction: Normal
% direction, which
% corresponds to the directions of the normal vectors to the object's surface,
% and View direction, which corresponds to the view direction of the z-buffer.
% The first option requires the input file of normal vectors, which is an
% object w x h x 1 x 1 x 3 that can be generated optionally by the operators
% kvsnorm or kisnorm,
% containing the X, Y and Z components of the vectors, stored with inverted
% signals (vectors pointing to inside the object).  The second option requires
% the angles alpha and beta, used to obtain the z-buffer.
%
% There are two options for the range of "texture" extraction: along a region
% around the
% surface, specified as how many steps to walk backwards (outside the object)
% and forwards (inside the object); and
% along the thickness, which is determined by the specified threshold. The
% later option requires the direction of "texture" extraction to be the normal
% direction.
%
% The increment step can be specified and represents the sampling step while
% extracting the "texture". An increment step equal to 1 visits the voxels one
% by one.
%
% There are five types of texture: 1) the voxels' average value, 2) the voxels'
% weighted average value, 3) the voxels' maximal value, 4) the voxels' minimal
% value and 5) the thickness.
%
% The first option calculates the average value of all the visited voxels. The
% second option is still not implemented. The third option calculates the
% maximal voxel value of all the visited voxels. The forth option calculates
% the minimal voxel value of all the visited voxels. Finally, the last option
% estimate the thickness of the surface based on the length of the ray within
% the object. This last option requires the direction of extraction to be the
% normal direction and the range of extraction to be given by the threshold.
%
% The images generated by "ktextu" are useful for analysis by themselves,
% or can be used as input to the operator
% kshad, to compose a shaded image based not only on the object surface form,
% but also on its volumetric structure information.
%
%
%
%  EXAMPLES
% ktextu -i1 volume.viff -i2 coord.viff -o texture.viff -step 1.0 -norm_dir
%        -i3 normal.viff -out_range 10.0 -in_range 10.0 -max_type
%
% This example will create a texture file texture.viff from the volume file
% volume.viff and the coordinate file coord.viff by walking 10 steps
% backwards (outside the object) and 10 steps forwards (inside the object)
% around the surface, in the normal direction
% (which is defined by the normal vector file normal.viff), with a sampling
% step
% equal to 1. The specified texture type is the voxels' maximal value.
%
% ktextu -i1 volume.viff -i2 coord.viff -o texture.viff -step 2.0 -view_dir
%        -alpha 90.0 -beta 15.0 -out_range 10.0 -in_range 10.0 -ave_type
%
% This example will create a texture file texture.viff from the volume file
% volume.viff and the coodinate file coord.viff by walking 10 steps
% backwards and 10 steps forwards around the surface, in the view direction
% (which is a plane rotated 90 degrees in alpha and 15 degrees in beta), with a
% sampling step equal to 2. The specified texture type is the voxels' average
% value.
%
% ktextu -i1 volume.viff -i2 coord.viff -o texture.viff -step 1.0 -norm_dir
%        -i3 normal.viff -thres 70 -thi_type
%
% This example will create a texture file texture.viff from the volume file
% volume.viff and the coodinate file coord.viff by walking
% in the normal direction (which is defined by the normal vector
% file normal.viff), with a sampling step equal to 1. The specified texture type is the thickness, while the threshold is less than 70.
%
%  "SEE ALSO"
% kzbuff, kshad, kisnorm, kvsnorm, kvoxext
%
%  RESTRICTIONS
% The input objects must have only the value segment.
%
% The
% input object scene can not have dimention e > 1. The input objects coord and
% norm must have dimention e=3.
%
% The input object scene can not be of data types KCOMPLEX and KDCOMPLEX. The
% input object coord can not be of data types KBIT, KFLOAT, KDOUBLE, KCOMPLEX
% and KDCOMPLEX. The input object scene can not be of data types KBIT, KCOMPLEX
% and KDCOMPLEX.
%
% The input objects coord and norm must match in their dimentions w and h.
%
% In case of t > 1 in the input objects, the operator will be applied to the time
% t=0 only.
%
% None of the input and output objects are referenced, therefore some attributes
% may change, as the VALUE_POSITION, for example.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993, 1994, 1995 UNICAMP, R A Lotufo,  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function proposals = gen_proposals(img, map, resizeRatio, minmaxRF)
    global param
    global globalVar

    %% bw prob map
    f = fspecial('gaussian',[5 5],7);
    map = imfilter(map,f,'same');

    bwmap = map > param.minRegionProb;

    SE1=strel('disk',3);
    bwmap = imerode(bwmap,SE1);
    bwmap = imdilate(bwmap,SE1);
    bwmap = bwmap>0;
    bwmap = imfill(bwmap,'holes');

    %% process each region
    [L, L_num] = bwlabel(bwmap, 8);
    proposals = cell(L_num, 1);
    estimated_orientations = zeros(L_num, 1);
    regions = cell(L_num, 1);
    regionPerims = cell(L_num, 1);
    totalTimeOfGetCompOfRegion = 0;
    totalTimeOfEstimateOrientation = 0;
    totalTimeOfGenProposalsByOrientation = 0;
    for i = 1 : L_num
        regionMap = L == i;

        %% construct regions
        [y, x] = find(regionMap);
        regions{i} = [y, x];

        %% Substruct sub img and sub region map by regions
        regionMinX= min(x);
        regionMaxX = max(x);
        regionMinY = min(y);
        regionMaxY = max(y);
        regionHeight = regionMaxY - regionMinY + 1;
        regionWidth = regionMaxX - regionMinX + 1;

        extendRegionMinX = max(1, regionMinX - regionWidth * 0.005);
        extendRegionMaxX = min(size(img, 2), regionMaxX + regionWidth * 0.005);
        extendRegionMinY = max(1, regionMinY - regionHeight * 0.005);
        extendRegionMaxY = min(size(img, 1), regionMaxY + regionHeight * 0.005);

        subImg = img(extendRegionMinY : extendRegionMaxY, extendRegionMinX : extendRegionMaxX, :);
        subRegionMap = regionMap(extendRegionMinY : extendRegionMaxY, extendRegionMinX : extendRegionMaxX);
        [subY, subX] = find(subRegionMap);
        subRegions = [subY, subX];

        %% Get comp info from subImg
        %comp_infos = normal_mser3(subImg, param.mser_info, resizeRatio, minmaxRF);
        comp_infos = normal_mser(subImg, resizeRatio, minmaxRF, [globalVar.imgName, '_', num2str(i)]);

        %% ---- for debug ----
        if(param.debug)
            boxes = zeros(length(comp_infos), 4);
            for n = 1 : length(comp_infos)
                boxes(n,:) = comp_infos{n}.box;
            end
            show_bbox(subImg, boxes);
        end

        %% construct regionPerims
        subRegionMap = imfill(subRegionMap, 'holes');
        perimMap = bwperim(subRegionMap);
        [y, x] = find(perimMap);
        subRegionPerims = [y, x];

        tic;
        [region_comp_infos] = getCompOfRegion(subImg, subRegions, comp_infos, param.minRegionCompCoveredArea, param.maxRegionCompArea);
        tmpTime = toc;
        totalTimeOfGetCompOfRegion = totalTimeOfGetCompOfRegion + tmpTime;

        %% generate secondary_region_comp_infos
%         [~, secondary_region_comp_idx] = getCompOfRegion(img, regions{i}, comp_infos, param.secondaryMinRegionCompArea);
%         secondary_region_comp_flags = false(length(comp_infos), 1);
%         secondary_region_comp_flags(secondary_region_comp_idx) = true;
%         secondary_region_comp_flags(region_comp_idx) = false;
%         secondary_region_comp_infos = num2cell(comp_infos(secondary_region_comp_flags));

        if(param.debug)
            debug_bbox = zeros(length(region_comp_infos), 4);
            for d_i = 1 : length(region_comp_infos)
                debug_bbox(d_i,:) = region_comp_infos{d_i}.box;
            end
            show_bbox(subImg, debug_bbox);
        end

        if(false && param.debug)
            debug_bbox = zeros(length(secondary_region_comp_infos), 4);
            for d_i = 1 : length(secondary_region_comp_infos)
                debug_bbox(d_i,:) = secondary_region_comp_infos{d_i}.box;
            end
            show_bbox(subImg, debug_bbox);
        end

        tic;
        estimated_orientations(i) = estimateOrientation(subImg, subRegions, region_comp_infos, param.orient_param);
        tmpTime = toc;
        totalTimeOfEstimateOrientation = totalTimeOfEstimateOrientation + tmpTime;

        tic;
        proposals{i} = genProposalsByOrientation(subImg, subRegions, subRegionPerims, estimated_orientations(i), region_comp_infos);
        tmpTime = toc;
        totalTimeOfGenProposalsByOrientation = totalTimeOfGenProposalsByOrientation + tmpTime;

        proposals{i}(:, 1 : 2 : 8) = proposals{i}(:, 1 : 2 : 8) + extendRegionMinX;
        proposals{i}(:, 2 : 2 : 8) = proposals{i}(:, 2 : 2 : 8) + extendRegionMinY;
    end

    if(false)
        fprintf('TotalTime: %.2f, %.2f, %.2f\n', ...
            totalTimeOfGetCompOfRegion, ...
            totalTimeOfEstimateOrientation, ...
            totalTimeOfGenProposalsByOrientation);
    end

    proposals = cell2mat(proposals);
end

author: stupidZZ


text:
function model = circulant_rand(d)
% random circulant embedding
% d dim of feature
% model.r  circulant vector from gaussian distribution
% model.bernoulli random bernoulli vector for sign flipping

    model.r = randn(d,1);
    rr = randn(1,d);
    rr(rr > 0) = 1;
    rr(rr <= 0 ) = -1;
    model.bernoulli = rr;
end

author: willard-yuan

repo: hashing-baseline-for-image-retrieval