text:
function y = norms( x, p, dim )

%NORMS   Internal cvx version.

%
% Size check
%

narginchk(1,3);
try
    if nargin < 3, dim = []; end
    [ x, sx, sy, zx, zy, nx, nv, perm ] = cvx_reduce_size( x, dim ); %#ok
catch exc
    error( exc.message );
end

%
% Check second argument
%

if nargin < 2 || isempty( p ),
    p = 2;
elseif ~isnumeric( p ) || numel( p ) ~= 1 || ~isreal( p ),
    error( 'Second argument must be a real number.' );
elseif p < 1 || isnan( p ),
    error( 'Second argument must be between 1 and +Inf, inclusive.' );
end

%
% Quick exit for empty matrices
%

if isempty( x ),
    y = zeros( zy );
    return
end

%
% Type check
%

persistent remap1 remap2
if isempty( remap2 ),
    remap1 = cvx_remap( 'constant', 'log-convex' );
    remap2 = cvx_remap( 'affine', 'log-convex' );
end
xc = reshape( cvx_classify( x ), sx );
if ~all( remap2( xc( : ) ) ),
    error( 'Disciplined convex programming error:\n   Invalid computation: norms( {%s}, ... )', cvx_class( x, true, true ) );
end

%
% Compute norms
%

if nx == 1,
	p = 0;
end
switch p,
	case 0,
		y = abs( x );
    case 1,
        y = sum( abs(x) );
    case Inf,
        y = max( abs(x) );
    otherwise,
        tt = all( remap1( xc ) );
        if all( tt( : ) ),
            y = sum( abs(x) .^ p ) .^ (1/p);
        elseif any( tt( : ) ),
            y  = cvx( [ 1, nv ], [] );
            xt = cvx_subsref( x, ':', tt );
            y  = cvx_subsasgn( y, tt, norms( xt, p ) );
            tt = ~tt;
            xt = cvx_subsref( x, ':', tt );
            y  = cvx_subsasgn( y, tt, norms( xt, p ) );
        elseif p == 2,
        	y = [];
            cvx_begin
                epigraph variable y( 1, nv )
                { cvx_accept_convex(x), y } == lorentz( [ nx, nv ], 1, ~isreal( x ) ); %#ok
            cvx_end
		else
			z = []; y = [];
            cvx_begin
                variable z( nx, nv )
                epigraph variable y( 1, nv )
                if isreal(x), cmode = 'abs'; else cmode = 'cabs'; end
                { cat( 3, z, y( ones(nx,1), : ) ), cvx_accept_convex(x) } ...
                    == geo_mean_cone( sw, 3, [1/p,1-1/p], cmode ); %#ok
                sum( z ) == y; %#ok
            cvx_end
        end
end

%
% Reverse the reshaping and permutation steps
%

y = reshape( y, sy );
if ~isempty( perm ),
    y = ipermute( y, perm );
end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [embedding_layer_state, hidden_layer_state, output_layer_state] = ...
  fprop(input_batch, word_embedding_weights, embed_to_hid_weights,...
  hid_to_output_weights, hid_bias, output_bias)
% This method forward propagates through a neural network.
% Inputs:
%   input_batch: The input data as a matrix of size numwords X batchsize where,
%     numwords is the number of words, batchsize is the number of data points.
%     So, if input_batch(i, j) = k then the ith word in data point j is word
%     index k of the vocabulary.
%
%   word_embedding_weights: Word embedding as a matrix of size
%     vocab_size X numhid1, where vocab_size is the size of the vocabulary
%     numhid1 is the dimensionality of the embedding space.
%
%   embed_to_hid_weights: Weights between the word embedding layer and hidden
%     layer as a matrix of size numhid1*numwords X numhid2, numhid2 is the
%     number of hidden units.
%
%   hid_to_output_weights: Weights between the hidden layer and output softmax
%               unit as a matrix of size numhid2 X vocab_size
%
%   hid_bias: Bias of the hidden layer as a matrix of size numhid2 X 1.
%
%   output_bias: Bias of the output layer as a matrix of size vocab_size X 1.
%
% Outputs:
%   embedding_layer_state: State of units in the embedding layer as a matrix of
%     size numhid1*numwords X batchsize
%
%   hidden_layer_state: State of units in the hidden layer as a matrix of size
%     numhid2 X batchsize
%
%   output_layer_state: State of units in the output layer as a matrix of size
%     vocab_size X batchsize
%

[numwords, batchsize] = size(input_batch);
[vocab_size, numhid1] = size(word_embedding_weights);
numhid2 = size(embed_to_hid_weights, 2);

%% COMPUTE STATE OF WORD EMBEDDING LAYER.
% Look up the inputs word indices in the word_embedding_weights matrix.
embedding_layer_state = reshape(...
  word_embedding_weights(reshape(input_batch, 1, []),:)',...
  numhid1 * numwords, []);

%% COMPUTE STATE OF HIDDEN LAYER.
% Compute inputs to hidden units.
inputs_to_hidden_units = embed_to_hid_weights' * embedding_layer_state + ...
  repmat(hid_bias, 1, batchsize);

% Apply logistic activation function.
% FILL IN CODE. Replace the line below by one of the options.
hidden_layer_state = zeros(numhid2, batchsize);
% Options
% (a) hidden_layer_state = 1 ./ (1 + exp(inputs_to_hidden_units));
% (b) hidden_layer_state = 1 ./ (1 - exp(-inputs_to_hidden_units));
hidden_layer_state = 1 ./ (1 + exp(-inputs_to_hidden_units));
% (d) hidden_layer_state = -1 ./ (1 + exp(-inputs_to_hidden_units));

%% COMPUTE STATE OF OUTPUT LAYER.
% Compute inputs to softmax.
% FILL IN CODE. Replace the line below by one of the options.
inputs_to_softmax = zeros(vocab_size, batchsize);
% Options
inputs_to_softmax = hid_to_output_weights' * hidden_layer_state +  repmat(output_bias, 1, batchsize);
% (b) inputs_to_softmax = hid_to_output_weights' * hidden_layer_state +  repmat(output_bias, batchsize, 1);
% (c) inputs_to_softmax = hidden_layer_state * hid_to_output_weights' +  repmat(output_bias, 1, batchsize);
% (d) inputs_to_softmax = hid_to_output_weights * hidden_layer_state +  repmat(output_bias, batchsize, 1);

% Subtract maximum.
% Remember that adding or subtracting the same constant from each input to a
% softmax unit does not affect the outputs. Here we are subtracting maximum to
% make all inputs <= 0. This prevents overflows when computing their
% exponents.
inputs_to_softmax = inputs_to_softmax...
  - repmat(max(inputs_to_softmax), vocab_size, 1);

% Compute exp.
output_layer_state = exp(inputs_to_softmax);

% Normalize to get probability distribution.
output_layer_state = output_layer_state ./ repmat(...
  sum(output_layer_state, 1), vocab_size, 1);

author: BradNeuberg

repo: hinton-coursera


text:
function unicycle_check ( n, p )

%*****************************************************************************80
%
%% UNICYCLE_CHECK checks that a vector represents a unicycle.
%
%  Discussion:
%
%    A unicycle is a permutation with a single cycle.  This might be called
%    a cyclic permutation, except that that name is used with at least two
%    other meanings.  So, to be clear, a unicycle is a permutation of N
%    objects in which each object is returned to itself precisely after
%    N applications of the permutation.
%
%    This routine verifies that each of the integers from 1
%    to N occurs among the N entries of the permutation.
%
%    Any permutation of the integers 1 to N describes a unicycle.
%    The permutation ( 3, 4, 2, 1 ) indicates that the unicycle
%    sends 3 to 4, 4 to 2, 2 to 1 and 1 to 3.  This is the sequential
%    description of a unicycle.
%
%    The standard sequence "rotates" the permutation so that it begins
%    with 1.  The above sequence becomes a standard sequence when
%    written as ( 1, 3, 4, 2 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries.
%
%    Input, integer P(N), the unicycle sequence vector.
%
  ierror = 0;

  for iseek = 1 : n

    ierror = iseek;

    for ifind = 1 : n
      if ( p(ifind) == iseek )
        ierror = 0
        break
      end
    end

    if ( ierror ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'UNICYCLE_CHECK - Fatal error!\n' );
      fprintf ( 1, '  The input array does not represent\n' );
      fprintf ( 1, '  a unicycle.  In particular, the\n' );
      fprintf ( 1, '  array is missing the value %d\n', ierror );
      error ( 'UNICYCLE_CHECK - Fatal error!' )
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [disp_row, disp_col, sind] = resp_newton(response, responsef, iterations, ky, kx, use_sz)

[max_resp_row, max_row] = max(response, [], 1);
[init_max_response, max_col] = max(max_resp_row, [], 2);
max_row_perm = permute(max_row, [2 3 1]);
col = max_col(:)';
row = max_row_perm(sub2ind(size(max_row_perm), col, 1:size(response,3)));

trans_row = mod(row - 1 + floor((use_sz(1)-1)/2), use_sz(1)) - floor((use_sz(1)-1)/2);
trans_col = mod(col - 1 + floor((use_sz(2)-1)/2), use_sz(2)) - floor((use_sz(2)-1)/2);
init_pos_y = permute(2*pi * trans_row / use_sz(1), [1 3 2]);
init_pos_x = permute(2*pi * trans_col / use_sz(2), [1 3 2]);
max_pos_y = init_pos_y;
max_pos_x = init_pos_x;

% pre-compute complex exponential
exp_iky = exp(bsxfun(@times, 1i * ky, max_pos_y));
exp_ikx = exp(bsxfun(@times, 1i * kx, max_pos_x));

% gradient_step_size = gradient_step_size / prod(use_sz);

ky2 = ky.*ky;
kx2 = kx.*kx;

iter = 1;
while iter <= iterations
    % Compute gradient
    ky_exp_ky = bsxfun(@times, ky, exp_iky);
    kx_exp_kx = bsxfun(@times, kx, exp_ikx);
    y_resp = mtimesx(exp_iky, responsef, 'speed');
    resp_x = mtimesx(responsef, exp_ikx, 'speed');
    grad_y = -imag(mtimesx(ky_exp_ky, resp_x, 'speed'));
    grad_x = -imag(mtimesx(y_resp, kx_exp_kx, 'speed'));
    ival = 1i * mtimesx(exp_iky, resp_x, 'speed');
    H_yy = real(-mtimesx(bsxfun(@times, ky2, exp_iky), resp_x, 'speed') + ival);
    H_xx = real(-mtimesx(y_resp, bsxfun(@times, kx2, exp_ikx), 'speed') + ival);
    H_xy = real(-mtimesx(ky_exp_ky, mtimesx(responsef, kx_exp_kx, 'speed'), 'speed'));
    det_H = H_yy .* H_xx - H_xy .* H_xy;

    % Compute new position using newtons method
    max_pos_y = max_pos_y - (H_xx .* grad_y - H_xy .* grad_x) ./ det_H;
    max_pos_x = max_pos_x - (H_yy .* grad_x - H_xy .* grad_y) ./ det_H;

    % Evaluate maximum
    exp_iky = exp(bsxfun(@times, 1i * ky, max_pos_y));
    exp_ikx = exp(bsxfun(@times, 1i * kx, max_pos_x));

    iter = iter + 1;
end
max_response = 1 / prod(use_sz) * real(mtimesx(mtimesx(exp_iky, responsef, 'speed'), exp_ikx, 'speed'));

% check for scales that have not increased in score
ind = max_response < init_max_response;
max_response(ind) = init_max_response(ind);
max_pos_y(ind) = init_pos_y(ind);
max_pos_x(ind) = init_pos_x(ind);

[max_scale_response, sind] = max(max_response(:));
disp_row = (mod(max_pos_y(1,1,sind) + pi, 2*pi) - pi) / (2*pi) * use_sz(1);
disp_col = (mod(max_pos_x(1,1,sind) + pi, 2*pi) - pi) / (2*pi) * use_sz(2);
end

repo: STRCF


text:
function success=testMissingMovieTimeStamp()

% This tests to make sure that if an experiment has no movie, the
% JLabelData object represents that appropriately.

jabFileName='/groups/branson/bransonlab/projects/JAABA/test_data/larva_mwt_rolling.jab';
gtMode=false;
data=JLabelData('setstatusfn',@(str)(fprintf('%s\n',str)), ...
                'clearstatusfn',@()(nop()));
data.openJabFile(jabFileName,gtMode);
fileTypes=data.filetypes;
iFileType=whichstr('movie',fileTypes);
nExps=data.nexps;
for iExp=1:nExps
  if data.fileexists(iExp,iFileType) ,
    error('testMissingMovieTimeStamp:fileexistsSetWrong', ...
          'filexists is true where it should be false');
  end
  if ~isinf(data.filetimestamps(iExp,iFileType)) ,
    error('testMissingMovieTimeStamp:filetimestampSetWrong', ...
          'filetimestamps is finite where it should be -inf');
  end
end
data.closeJabFile();
data=[];  %#ok
success=true;

end

author: kristinbranson

repo: JAABA


text:
function out = DatesOption(fo_year,lo_year,frequency,fo_period,lo_period)
% =======================================================================
% For given initial observation, final observation and frequency this
% function computes the number of observations between them and creates a
% cell aray with the dates (of the type 1979M2).
% =======================================================================
% out = DatesOption(fo_year,lo_year,frequency,fo_period,lo_period)
% -----------------------------------------------------------------------
% INPUT
%	- fo_year: initial year of the timeline
%	- lo_year: final year of the timeline
%	- frequency: quarterly 'q' [def], monthly 'm', yearly 'y'
% -----------------------------------------------------------------------
% OPTIONAL INPUT
%   - fo_period: initial quarter/month of the timeline (not for
%       yearly frequency)
%   - lo_period: last quarter/month of the timeline (not for yearly
%       frequency)
% -----------------------------------------------------------------------
% OUTPUT
%	- out: structuree with all info
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com

% Check inputs
if ~exist('frequency','var')
    frequency = 'q';
end

if strcmp(frequency,'m')
    if ~exist('fo_period','var'), error('You need to provide first period'), end
    if ~exist('lo_period','var'), error('You need to provide last period'), end
    fo = fo_year + fo_period*(1/12);
    lo = lo_year + lo_period*(1/12);
    aux = fo:(1/12):lo;
    nobs = length(aux);
elseif strcmp(frequency,'q')
    if ~exist('fo_period','var'), error('You need to provide first period'), end
    if ~exist('lo_period','var'), error('You need to provide last period'), end
    fo = fo_year + fo_period*(1/4);
    lo = lo_year + lo_period*(1/4);
    aux = fo:(1/4):lo;
    nobs = length(aux);
elseif strcmp(frequency,'y')
    fo = fo_year;
    lo = lo_year;
    aux = fo:lo;
    nobs = length(aux);
end

% Get first
if strcmp(frequency,'y')
    [out.dates, out.dates_short] = DatesCreate(fo_year,nobs,frequency);
else
    [out.dates, out.dates_short] = DatesCreate(fo_year,nobs,frequency,fo_period);
end

out.nobs = nobs;
out.fo_year = fo_year;
out.lo_year = lo_year;
out.frequency = frequency;
out.fo_period = fo_period;
out.lo_period = lo_period;

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function new_regions = split_regions(regions, bndinfo, sp_areas)

spLR = bndinfo.edges.spLR;

n_sps = bndinfo.nseg;

if(iscell(regions))
   % get overlapping regions
   overlapping = [];

   sp_map = zeros(numel(regions), n_sps);

   for prop = 1:length(regions)
      sp_map(prop,regions{prop}) = 1;
   end
else
   sp_map = regions;
end

inds = double([spLR(:,[1 2]); spLR(:, [2 1]); repmat([1:n_sps]', 1, 2)]);

%sp_map_sparse = sparse(sp_map);
%sparse_adjacency = sparse(inds(:,1), inds(:,2), ones(2*size(spLR,1), 1));
adjacency = (accumarray(inds, ones(size(inds,1), 1), [n_sps n_sps]));
%adjacency_s = sparse(accumarray(inds, ones(size(inds,1), 1), [n_sps n_sps]));

[adjacency_x adjacency_y adjacency_val] = find(adjacency);

new_regions = cell(size(sp_map,1), 1);
new_r_ind = 1;
tic;
for i = 1:size(sp_map,1)
   sp_inds = find(sp_map(i, :));
%   sub_graph = sparse(adjacency(sp_inds, sp_inds));

   ok = ismember_sorted(adjacency_x, sp_inds) & ismember_sorted(adjacency_y, sp_inds);
   ind = zeros(max(adjacency_x),1);
   ind(sp_inds) = 1:length(sp_inds);
   sub_graph = sparse(ind(adjacency_x(ok)), ind(adjacency_y(ok)), adjacency_val(ok), length(sp_inds), length(sp_inds));

   [conn1 comps1] = graphconncomp(sub_graph, 'Directed', true);

   if(conn1 > 1)
      if(toc>1)
         fprintf('%d (%d/%d)\n',new_r_ind, i, size(sp_map,1) );
         tic;
      end

      for j = 1:conn1
         if(new_r_ind>=numel(new_regions))
            fprintf('growing\n')
            new_regions{2*end}= [];
         end
         if(sum(sp_areas(sp_inds(comps1==j)))>(24*24))
            new_regions{new_r_ind} = sp_inds(comps1==j);
            new_r_ind = new_r_ind+1;
         end
      end
%      pause
   else
      if(new_r_ind>=numel(new_regions))
         fprintf('growing\n')
         new_regions{2*end}= [];
      end

      new_regions{new_r_ind} = sp_inds;
      new_r_ind = new_r_ind+1;
   %      pause(0.3)
   end

end

new_regions(new_r_ind+1:end) = [];

author: Cloud-CV

repo: object-proposals


text:
function p01_story ( )

%*****************************************************************************80
%
%% P01_STORY prints the "story" for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This example is due to deBoor.\n' );
  fprintf ( 1, '  For this example, X is measured in days,\n' );
  fprintf ( 1, '  and Y records the observed position of Mars in a heliocentric\n' );
  fprintf ( 1, '  coordinate system.\n' );

  return
end


text:
function M=Euler1Ang2RotMat(theta,series,handed)
%%EULER1ANG2ROTMAT  Obtain the Euler-angle rotation matrix for a rotation
%            of angle theta about single axis (x, y, or z). The rotation is
%            either right-handed (counterclockwise, the default) or
%            left-handed depending on handed. The components of vectors to
%            be rotated are assumed ordered [x;y;z].
%
%INPUTS: theta The angle in radians about which the rotation matrix should
%              rotate a vector around the axis given by series. The
%              handedness of the angle is given by handed.
%       series A text string specifying the angle about which the rotation
%              is performed. This can be 'x', 'y', or 'z'.
%       handed The handedness of the rotation angle. If omitted, it is
%              assumed that the rotation is right-handed (the standard).
%              Possible values are:
%              'right' The default if omitted. The rotation is right-
%                      handed.
%              'left'  The rotation is left-handed. The rotation angle is
%                      clockwise when one is looking into the rotation
%                      axis.
%
%OUTPUTS: M The rotation matrix such that M*v rotates a vector v by theta
%           about the specified axis.
%
%Euler angles are discussed in [1].
%
%REFERENCES:
%[1] M. D. Shuster, "A survey of attitude representations," The Journal of
%    Astronautical Sciences, vol. 41, no. 4, pp. 439-517, Oct. - Dec. 1993.
%
%September 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

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

switch(handed)
    case 'right'
    case 'left'
        theta=-theta;
    otherwise
        error('Invalid handedness provided.')
end

cosT=cos(theta);
sinT=sin(theta);

switch(series)
    case 'x'
        M=[1, 0,       0;
           0, cosT, -sinT;
           0, sinT,  cosT];
    case 'y'
        M=[cosT, 0,  sinT;
              0, 1,     0;
          -sinT, 0, cosT];
    case 'z'
        M=[cosT, -sinT, 0;
           sinT,  cosT, 0;
              0,    0, 1];
    otherwise
        error('Invalid rotation axis specified')
end

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function example2 ( )

%*****************************************************************************80
%
%% EXAMPLE2 demonstrates the use of PDEPE on a scalar PDE.
%
%  Discussion:
%
%    Solve the convection-diffusion equation.
%
%    PDE:
%      ut + (a(x)*u)x = uxx
%    BC:
%      u(t,-oo) = 0, u(t,+oo) = 0
%    IC:
%      u(0,x) = 1 / (1+(x-5)^2)
%
%    Here, the quantity a(x) is a bit complicated.  It is defined by
%      a(x) = 3 (ubar(x))^2 + 2 * ubar(x)
%    where ubar is defined implicitly by the nonlinear equation
%      1/ubar + log ( abs ( (1-ubar)/ubar ) ) = x
%
%    Ubar(x) is the equilibrium solution to the conservation law
%      ut + (u^3-u^2)x = uxx.
%
%    Although the mathematical problem has boundary conditions at infinity,
%    we approximate this by the interval [-50,+50].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2013
%
%  Author:
%
%    Original formulation by P Howard, 2005.
%    This version by John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXAMPLE2:\n' );
  fprintf ( 1, '  The convection-diffusion equation.\n' );
  fprintf ( 1, '  ut + (a(x)*u)x = uxx\n' );
  fprintf ( 1, '  u(t,-oo) = 0, u(t,+oo) = 0\n' );
  fprintf ( 1, '  u(0,x) = 1 / (1+(x-5)^2)\n' );
%
%  M defines the coordinate system:
%  0: cartesian
%  1: cylindrical
%  2: spherical
%
  m = 0;
%
%  Define the spatial mesh.
%
  nx = 201;
  xmesh = linspace ( -50.0, +50.0, nx );
%
%  Define the time mesh.
%
  nt = 101;
  tspan = linspace ( 0.0, 10.0, nt );
%
%  Call PDEPE() for the solution.
%
  sol = pdepe ( m, @pdefun, @icfun, @bcfun, xmesh, tspan );
%
%  Even though our solution is "really" a 2D array, PDEPE stores it
%  in a 3D array SOL(:,:,:).  The surf() command needs a 2D array to plot,
%  so let's copy U out of SOL.
%
  u = sol(:,:,1);

  figure ( 1 )
  surf ( xmesh, tspan, u, 'EdgeColor', 'None' );
  title ( 'Example 2: Solution Over Time', 'Fontsize', 16 );
  xlabel ( '<--- X --->' )
  ylabel ( '<--- T --->' );
  zlabel ( '<---U(X,T)--->' );
  filename = 'example2.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved solution plot in file "%s"\n', filename );
%
%  Plot the initial condition, U at time 0.
%
  figure ( 2 )
  plot ( xmesh, u(1,:), 'LineWidth', 3 );
  grid on
  title ( 'Example 2: Initial Condition', 'Fontsize', 16 );
  xlabel ( '<--- X --->' )
  ylabel ( '<--- U(X,T0) --->' );
  filename = 'example2_ic.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Saved initial condition plot in file "%s"\n', filename );
%
%  Plot the solution U at a fixed point, with time varying.
%
  figure ( 3 )
  mid = 1 + floor ( 55 * ( nx - 1 ) / 100 );
  plot ( tspan, u(:,mid), 'LineWidth', 3 );
  grid on
  title ( 'Example 2: Time evolution of solution at X=5.0', 'Fontsize', 16 );
  xlabel ( '<--- T --->' )
  ylabel ( '<--- U(5.0,T) --->' );
  filename = 'example2_profile.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Saved time evolution plot in file "%s"\n', filename );
%
%  Animate the profile.
%  I wish I could also display the running value of time, but it
%  does not seem possible.
%
  figure ( 4 )
  fig = plot ( xmesh, u(1,:), 'erase', 'xor' );
  title ( 'Profile animation', 'Fontsize', 16 );
  grid on
  for k = 2 : length ( tspan )
    set ( fig, 'xdata', xmesh, 'ydata', u(k,:) );
    pause ( 0.4 );
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXAMPLE2:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function value = degwave ( x )

%*****************************************************************************80
%
%% DEGWAVE determines the value of UBAR in the equation.
%
%  Discussion:
%
%    MATLAB's zero finder "fzero()" must be called in order to determine
%    the value of UBAR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the spatial location.
%
%    Output, real VALUE, the value of A(X).
%
  if ( x < -35.0 )
    value = 1.0;
  elseif ( 2.0 < x )
    guess = 1.0 / x;
    value = fzero ( @f, guess, [], x );
  elseif ( -2.5 < x )
    guess = 0.6;
    value = fzero ( @f, guess, [], x );
  else
    guess = 1.0 - exp ( - 2.0 ) * exp ( x );
    value = fzero ( @f, guess, [], x );
  end

  return
end
function value = f ( u, x )

%*****************************************************************************80
%
%% F evaluates a function, which should be zero if U = UBAR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real U, the estimated solution value at X.
%
%    Input, real X, the spatial location.
%
%    Output, real VALUE, the function value.
%
  value = ( 1.0 / u ) + log ( ( 1.0 - u ) / u ) - x;

  return
end
function [ c, f, s ] = pdefun ( x, t, u, dudx )

%*****************************************************************************80
%
%% PDEFUN defines the components of the PDE.
%
%  Discussion:
%
%    The PDE has the form:
%
%      c * du/dt = x^(-m) d/dx ( x^m f ) + s
%
%    where m is 0, 1 or 2,
%    c, f and s are functions of x, t, u, and dudx,
%    and most typically, f = dudx.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the spatial location.
%
%    Input, real T, the current time.
%
%    Input, real U(:,1), the estimated solution at T and X.
%
%    Input, real DUDX(:,1), the estimated spatial derivative of U at T and X.
%
%    Output, real C(:,1), the coefficients of du/dt.
%
%    Output, real F(:,1), the flux terms.
%
%    Output, real S(:,1), the source terms.
%
  c = 1.0;
  ubar = degwave ( x );
  f = dudx - ( 3.0 * ubar.^2 - 2.0 * ubar ) * u;
  s = 0.0;

  return
end
function u0 = icfun ( x )

%*****************************************************************************80
%
%% ICFUN defines the initial conditions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the spatial location.
%
%    Output, real U0(:,1), the value of the solution at the initial time,
%    and location X.
%
  u0 = 1.0 / ( 1.0 + ( x - 5.0 ).^2 );

  return
end
function [ pl, ql, pr, qr ] = bcfun ( xl, ul, xr, ur, t )

%*****************************************************************************80
%
%% BCFUN defines the boundary conditions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real XL, the spatial coordinate of the left boundary.
%
%    Input, real UL(:,1), the solution estimate at the left boundary.
%
%    Input, real XR, the spatial coordinate of the right boundary.
%
%    Input, real UR(:,1), the solution estimate at the right boundary.
%
%    Output, real PL(:,1), the Dirichlet portion of the left boundary condition.
%
%    Output, real QL(:,1), the coefficient of the flux portion of the left
%    boundary condition.
%
%    Output, real PR(:,1), the Dirichlet portion of the right boundary condition.
%
%    Output, real QR(:,1), the coefficient of the flux portion of the right
%    boundary condition.
%
  pl = ul;
  ql = 0.0;
  pr = ur;
  qr = 0.0;

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% 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


text:
function fileList = GetAllFiles(dirName)

dirData = dir(dirName);
dirIndex = [dirData.isdir];
fileList = {dirData(~dirIndex).name}';
if ~isempty(fileList)
    fileList = cellfun(@(x) fullfile(dirName,x),...
        fileList,'UniformOutput',false);
end
subDirs = {dirData(dirIndex).name};
validIndex = ~ismember(subDirs,{'.','..'});
for iDir = find(validIndex)
    nextDir = fullfile(dirName,subDirs{iDir});
    fileList = [fileList; GetAllFiles(nextDir)];
end

end

author: decouples


text:
function [mov ok aviPath] = aviSave(M, aviPath, fps, varargin);
% Export a 3D matrix to an AVI file.
%
% USAGE:
%	[mov ok aviPath] = aviSave(M, [aviPath], [fps], [options]);
%
% INPUTS:
% M: image matrix. Can either be a 3D matrix, in which every slice is a
% frame of the movie (converted to uint8), or a 4D M x N x 3 x frames matrix,
% in which each 3D subvolume is a truecolor ([R G B] slice) image.
%
% aviPath: save path for the .avi file. Will append .avi if it's not there.
% [pops up a dialog if omitted.]
%
% fps: frames per second for the movie. [Default 6]
%
% Options include:
%	'compression': select a compression for the movie. [Defaults based on
%	computer type: 'Indeo5' for windows OS, 'none' otherwise.]
%
%	'cmap': colormap for converting 3D matrices to color.
%
%	'description': descriptive name for the AVI file (will be saved as
%	metadata.)
%
%	other options found in HELP MOVIE2AVI will be supported here.
%
% OUTPUTS:
%	mov: matlab movie structure of the movie.
%
%	ok: flag indicating whether the save was successful.
%
%	aviPath: save path (useful in case a dialog is brought up.)
%
% ras, 07/18/2008.
if notDefined('M'),		error('Need image matrix.');		end
if notDefined('fps'),	fps = 6;							end
if notDefined('aviPath')
	aviPath = mrvSelectFile('w', 'avi', 'Save .AVI file as...');
end

%% default params
quality = 100;
description = aviPath;
keyframe = fps;
if ispc
	compression = 'Indeo3';
	cmap = gray(236);  % max allowed by compressor
else
	compression = 'None';
	cmap = gray(256);
end

%% parse options
for ii = 1:2:length(varargin)
	switch lower(varargin{ii})
		case 'compression', compression = varargin{ii+1};
		case 'description', description = varargin{ii+1};
		case {'cmap' 'colormap'}, cmap = varargin{ii+1};
		case 'quality', quality = varargin{ii+1};
		case 'keyframe', keyframe = varargin{ii+1};
	end
end

%% make the movie structure
if ndims(M)<=3
	nFrames = size(M, 3);

	% enforce uint8 data type
	nCmap = size(cmap, 1);
	if min(M(:)) < 1 | max(M(:)) > nCmap | any(mod(M(:), 1) ~= 0)
		warning('Matrix not uint8: rescaling to be uint8');
		M = uint8( rescale2(M, [], [0 nCmap-1]) );
	end

elseif ndims(M)==4
	if size(M, 3) ~= 3
		error('For 4-D matrices, third dimension must be [R G B] slices.');
	end
	nFrames = size(M, 4);

else
	error('Inappropriate movie size.');

end

for f = 1:nFrames
	if ndims(M) <= 3
		mov(f) = im2frame(M(:,:,f), cmap);
	else
		mov(f) = im2frame(M(:,:,:,f));
	end
end

%% export to AVI
% enfore .avi extension
[p f ext] = fileparts(aviPath);
if ~isequal( lower(ext), '.avi' )
	aviPath = fullfile(p, [f '.avi']);
end

% export
try
	movie2avi(mov, aviPath, 'FPS', fps, 'Compression', compression, ...
		  'Quality', quality, 'Videoname', description, 'Keyframe', keyframe);
	ok = 1;
	fprintf('[%s]: Saved movie as %s.\n', mfilename, aviPath);
catch
	disp(lasterr);
	ok = 0;
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function handle = plot_sequence_strip(sequence, trajectories, varargin)
% plot_sequence_strip Generate an preview of a sequence as a strip of frames
%
% The function generates dedicated tracker legend plot. Tracker labels and
% their symbols are ordered in a grid.
%
% Input:
% - trackers (cell): A cell array of tracker structures.
% - varargin[Handle] (handle): Plot on existing figure handle.
% - varargin[Visible] (boolean): Is the figure visible on the display.
% - varargin[Window] (double): Figure width hint.
% - varargin[Samples] (double): Figure height hint.
% - varargin[Scale] (integer): Number of rows in a grid.
% - varargin[Columns] (integer): Number of columns in a grid.
%
% Output:
% - handle (handle): A figure handle.
%

handle = -1;
visible = false;
trajectories_markers = {};
window = 120;
samples = 12;
scale = 1;
frame_numbers = false;

groundtruth_color = [1, 1, 1];
trajectories_colors = mat2cell(repmat([1, 0, 0], length(trajectories), 1), ...
    ones(length(trajectories), 1));

for i = 1:2:length(varargin)
    switch lower(varargin{i})
        case 'handle'
            handle = varargin{i+1};
        case 'visible'
            visible = varargin{i+1};
        case 'window'
            window = varargin{i+1};
        case 'samples'
            samples = varargin{i+1};
        case 'scale'
            scale = varargin{i+1};
        case 'framenumbers'
            frame_numbers = varargin{i+1};
        case 'groundtruthcolor'
            groundtruth_color = varargin{i+1};
        case 'trajectorycolor'
            trajectories_colors = varargin{i+1};
        case 'trajectorymarkers'
            trajectories_markers = varargin{i+1};
        otherwise
            error(['Unknown switch ', varargin{i},'!']) ;
    end
end

if numel(samples) > 1
    indices = samples(samples > 0 & samples <= sequence.length);
else
    indices = round(linspace(1, sequence.length, samples));
end;

if ishandle(handle)
    set(0, 'CurrentFigure', handle);
else
    if visible
        handle = figure();
    else
        handle = figure('Visible', 'off');
    end;
end;
handles = tight_subplots(1, length(indices), 0, 0, 0);

for i = 1:length(indices)

	if isinf(window)

		offset = [0, 0];

		patch = imread(sequence_get_image(sequence, indices(i)));

	else

	    patch = zeros(window, window, 3);

		region = region_convert(sequence_get_region(sequence, indices(i)), 'rectangle');

		offset = region(1:2) + (region(3:4) - window) / 2;

		source = imread(sequence_get_image(sequence, indices(i)));

		patch = patch_operation(patch, source, -offset([2, 1]), '=');

	end;

    axes(handles(i)); %#ok<LAXES>

    if ~visible
        set(handle, 'Visible', 'off');
    end;

    axis tight;

    imshow(uint8(patch));

    hold on;

    for t = 1:length(trajectories)

        region = trajectories{t}{indices(i)};

        if numel(region) < 2;
            continue;
        end;

        region = region_offset(region, -offset);

        region_draw(region, trajectories_colors{t}, 2);

        if ~isempty(trajectories_markers)
            bounds = region_convert(region, 'rectangle');
            center = bounds(1:2) + bounds(3:4) / 2;

            plot(center(1), center(2), trajectories_markers{t}, 'MarkerSize', 7, 'LineWidth', 1.5, 'Color', trajectories_colors{t});

        end;

    end;

    if ~isempty(groundtruth_color)

        region = sequence_get_region(sequence, indices(i));

        region = region_offset(region, -offset);

        region_draw(region, groundtruth_color, 2);

    end;

    if frame_numbers
        text(10, 10, sprintf('%d', indices(i)), 'Color', 'w', 'BackgroundColor', [0, 0, 0]);

    end;

    hold off;
end

width = numel(indices);
height = 1;

set(handle, 'PaperUnits', 'inches', 'PaperSize', [width, height] * scale, 'PaperPosition', [0, 0, width, height] * scale);

author: votchallenge

repo: toolkit-legacy


text:
function [fnames, defaults] = mysubs_fields

% function [fnames, defaults] = mysubs_fields
% Additional fields for class cfg_menu. See help of
% @cfg_item/subs_fields for general help about this function.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: mysubs_fields.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

fnames={'labels','values'};
defaults={{},{}};

author: spm


text:
function [ kSpace, measPara ] = convertKSpace_main( data )
% convert from ISMRMD to processable k-space format
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

iFreq = unique(double(data.head.number_of_samples));
iY = unique(double(data.head.idx.kspace_encode_step_1));
iZ = unique(double(data.head.idx.kspace_encode_step_2));
iPha = unique(double(data.head.idx.phase));
iSli = unique(double(data.head.idx.slice));
iRep = unique(double(data.head.idx.repetition));
iAvg = unique(double(data.head.idx.average));
iEcho = unique(double(data.head.idx.contrast));
iCha = unique(double(data.head.active_channels));
measPara.dim = [length(iY), iFreq, length(iZ), length(iPha), iCha];
measPara.LCall = [length(iSli), length(iAvg), length(iEcho), length(iRep)];

if(measPara.dim(3) == 1)
    measPara.dimension = '2D';
else
    if(measPara.dim(4) == 1)
        measPara.dimension = '3D';
    else
        measPara.dimension = '4D';
    end
end

% begin conversion
kSpace = cell(measPara.LCall(1), measPara.dim(5), measPara.LCall(4), measPara.LCall(2));
if(strcmp(measPara.dimension,'2D'))
    [kSpace{:}] = deal(complex(zeros(measPara.dim(1),iFreq,measPara.dim(4),'single'),zeros(measPara.dim(1),iFreq,measPara.dim(4),'single')));
elseif(strcmp(measPara.dimension,'3D'))
    [kSpace{:}] = deal(complex(zeros(measPara.dim(1),iFreq,measPara.dim(3),'single'),zeros(measPara.dim(1),iFreq,measPara.dim(3),'single')));
elseif(strcmp(measPara.dimension,'4D'))
    [kSpace{:}] = deal(complex(zeros(measPara.dim(1),iFreq,measPara.dim(3),measPara.dim(4),'single'),zeros(measPara.dim(1),iFreq,measPara.dim(3),measPara.dim(4),'single')));
end

for idxRep = 1:measPara.LCall(4)
    for idxAvg = 1:measPara.LCall(2)

        if(strcmp(measPara.dimension,'2D'))
            for idxSli = 1:measPara.LCall(1)
                for idxPha = 1:measPara.dim(4)
                    for idxY = 1:measPara.dim(1)
                        lIdx = double(data.head.idx.repetition) == iRep(idxRep) & ...
                            double(data.head.idx.average) == iAvg(idxAvg) & ...
                            double(data.head.idx.slice) == iSli(idxSli) & ...
                            double(data.head.idx.phase) == iPha(idxPha) & ...
                            double(data.head.idx.kspace_encode_step_1) == iY(idxY);
                        if(nnz(lIdx) > 1)
                            lIdx(find(lIdx,nnz(lIdx)-1,'last')) = false;
                        end
                        dataLine = data.data{lIdx};

                        for idxCha = 1:measPara.dim(5)
                            kSpace{idxSli,idxCha,idxRep,idxAvg}(idxY,:,idxPha) = dataLine(2*iFreq*(idxCha-1)+1:2:2*iFreq*idxCha).'+ 1i * dataLine(2*iFreq*(idxCha-1)+2:2:2*iFreq*idxCha).';
                        end
                    end
                end
            end

        elseif(strcmp(measPara.dimension,'3D') || strcmp(measPara.dimension,'4D'))
            for idxZ = 1:measPara.dim(3)
                for idxPha = 1:measPara.dim(4)
                    for idxY = 1:measPara.dim(1)
                        lIdx = double(data.head.idx.repetition) == iRep(idxRep) & ...
                            double(data.head.idx.average) == iAvg(idxAvg) & ...
                            double(data.head.idx.kspace_encode_step_2) == iZ(idxZ) & ...
                            double(data.head.idx.phase) == iPha(idxPha) & ...
                            double(data.head.idx.kspace_encode_step_1) == iY(idxY);
                        if(nnz(lIdx) > 1)
                            lIdx(find(lIdx,nnz(lIdx)-1,'last')) = false;
                        end
                        dataLine = data.data{lIdx};

                        for idxCha = 1:measPara.dim(5)
                            kSpace{1,idxCha,idxRep,idxAvg}(idxY,:,idxZ,idxPha) = dataLine(2*iFreq*(idxCha-1)+1:2:2*iFreq*idxCha).' + 1i * dataLine(2*iFreq*(idxCha-1)+2:2:2*iFreq*idxCha).';
                        end
                    end
                end
            end
        end
    end
end
end

author: thomaskuestner


text:
function rgb = white

% This returns a predefined color as [red green blue] values
%   red               = [255   0   0]/255;
%   green             = [  0 192   0]/255;
%   blue              = [  0   0 255]/255;
%   magenta           = [255 255   0]/255;
%   cyan              = [  0 255 255]/255;
%   yellow            = [255 255   0]/255;
%   white             = [255 255 255]/255;
%   black             = [  0   0   0]/255;
%
%   skull             = [140  85  85]/255
%   cortex            = [255 213 119]/255;
%   cortex_light      = [199 194 169]/255;
%   cortex_dark       = [100  97  85]/255;
%   skin              = [249 223 192]/255;
%   skin_light        = [249 223 192]/255;
%   skin_medium_light = [225 194 158]/255;
%   skin_medium       = [188 142 106]/255;
%   skin_medium_dark  = [155 102	65]/255;
%   skin_dark         = [ 91  71  61]/255;

rgb = [255 255 255]/255;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function best_acc = COPAR_top(dataset, N_train, k, k0, lambda, eta)
% * function `COPAR_top(dataset, N_train, k, k0, lambda, eta)`
% * The top function of COPAR
% * INPUT:
%   + `dataset`: name of the dataset stored in `.mat` file in `data` folder.
%     Note that `dataset` is the file name of the `.mat`, excluding `.mat`.
%   + `N_train`: number of training samples in each class
%   + `k`: number of bases in EACH PARTICULAR dictionary
%   + `k0`: number of bases in the COMMON dictionary
%   + `lambda, eta`: regularization parameters.
% * To run an small example, type `COPAR_top` without input in
%     MATLAB command window.
% -----------------------------------------------
% Author: Tiep Vu, thv102@psu.edu, 5/11/2016
%         (http://www.personal.psu.edu/thv102/)
% -----------------------------------------------
    addpath(genpath('utils'));
    addpath(genpath('DLSI'));
    addpath(genpath('COPAR'));
    addpath('ODL')
    %% test mode
    if nargin == 0
        dataset = 'myYaleB';
        N_train = 10;
        k = 8;
        k0 = 5;
        lambda = 0.001;
        eta = 0.01;
    end
    %% get data
    t = getTimeStr();
    [dataset, Y_train, Y_test, label_train, label_test] = ...
        train_test_split(dataset, N_train);
    %% main
    [acc, rt] = COPAR_wrapper(Y_train, label_train, Y_test , label_test, ...
                            k, k0, lambda, eta);
    disp(rt);
    %% output filename
    if ~exist('results', 'dir')
        mkdir('results');
    end
    if ~exist(fullfile('results', 'COPAR'), 'dir')
        mkdir('results', 'COPAR');
    end
    fn = fullfile('results', 'COPAR', strcat(dataset, '_N_', ...
        num2str(N_train), '_k_', num2str(k), '_k0_', num2str(k0), ...
        '_l_', num2str(lambda), '_e_', num2str(eta), '_', t, '.mat'));
    disp(fn);
    save(fn, 'acc', 'rt');
    best_acc = max(acc);
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function y = sinvchi2_pdf(x,nu,s2)
%SINVCHI2_PDF  Scaled inverse-chi probability density function.
%   Y = SINVCHI2_PDF(X,NU,S2) returns the scaled inverse-chi2 probability
%   density function with parameters NU and S2, at the values in X.
%
%   Note: Parameterization as in (Gelman et al, 1996).
%      NU is degrees of freedom
%      S2 is scale

% Copyright (c) 2003 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin < 2,
   error('Requires at least two input arguments.');
end

y = log(nu/2).*(nu/2) -gammaln(nu/2) + log(s2)/2*nu - log(x).*(nu/2+1) -nu.*s2/2./x;
y=exp(y);


text:
%%*****************************************************************
%% sqlp: main solver
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

  function [obj,X,y,Z,info,runhist] = sqlpmain(blk,At,C,b,par,parbarrier,X0,y0,Z0);

   global spdensity smallblkdim  printlevel msg
   global solve_ok  use_LU  exist_analytic_term numpertdiagschur
   global schurfun  schurfun_par
%%
   matlabversion = par.matlabversion;
   vers          = par.vers;
   predcorr      = par.predcorr;
   gam           = par.gam;
   expon         = par.expon;
   gaptol        = par.gaptol;
   inftol        = par.inftol;
   steptol       = par.steptol;
   maxit         = par.maxit;
   printlevel    = par.printlevel;
   stoplevel     = par.stoplevel;
   scale_data    = par.scale_data;
   spdensity     = par.spdensity;
   rmdepconstr   = par.rmdepconstr;
   smallblkdim   = par.smallblkdim;
   schurfun      = par.schurfun;
   schurfun_par  = par.schurfun_par;
   ublksize      = par.ublksize;
%%
   tstart = clock;
   X = X0; y = y0; Z = Z0;
   for p = 1:size(blk,1)
      if strcmp(blk{p,1},'u'); Z{p} = zeros(blk{p,2},1); end
   end
%%
%%-----------------------------------------
%% convert unrestricted blk to linear blk.
%%-----------------------------------------
%%
   convertlen = 0;
   [blk,At,C,X,Z,u2lblk,ublkidx] = sqlpu2lblk(blk,At,C,X,Z,par,convertlen);
   for p = 1:size(blk,1)
      pblk = blk(p,:);
      if (u2lblk(p) == 1)
         n = 2*blk{p,2};
         blk{p,1} = 'l';  blk{p,2} = n;
         parbarrier{p} = zeros(1,n);
         At{p} = [At{p}; -At{p}];
         tau = max(1,norm(C{p}));
         C{p} = [C{p}; -C{p}];
         msg = 'convert ublk to lblk';
         if (printlevel); fprintf(' *** %s',msg); end
         b2 = 1 + abs(b');
         normCtmp = 1+norm(C{p});
         normAtmp = 1+sqrt(sum(At{p}.*At{p}));
         if (n > 1000)
            const = sqrt(n);
         else
	    const = n;
         end
         if (par.startpoint == 1)
            X{p} = const* max([1,b2./normAtmp]) *ones(n,1);
            Z{p} = const* max([1,normAtmp/sqrt(n),normCtmp/sqrt(n)]) *ones(n,1);
            X{p} = X{p}.*(1+1e-10*randmat(n,1,0,'u'));
            Z{p} = Z{p}.*(1+1e-10*randmat(n,1,0,'u'));
	 else
            const = max(abs(X{p})) + 100;
            X{p} = [X{p}+const; const*ones(n/2,1)];
            %%old: const = 100; Z{p} = [const*ones(n/2,1); const*ones(n/2,1)];
            Z{p} = [abs(Z0{p}); abs(Z0{p})] + 1e-4;
         end
      end
   end
%%-----------------------------------------
%% check whether {A1,...,Am} is
%% linearly independent.
%%-----------------------------------------
%%
   m0 = length(b);
   [At,b,y,indeprows,par.depconstr,feasible,par.AAt] = ...
    checkdepconstr(blk,At,b,y,rmdepconstr);
   if (~feasible)
      obj = []; X = cell(size(blk,1),1); y = []; Z = cell(size(blk,1),1);
      runhist = [];
      msg = 'SQLP is not feasible';
      if (printlevel); fprintf('\n %s \n',msg); end
      return;
   end
   par.normAAt = norm(par.AAt,'fro');
%%
%%-----------------------------------------
%% scale SQLP data. Note: must be done only
%% after checkdepconstr
%%-----------------------------------------
%%
   normA2 = 1+ops(At,'norm');
   normb2 = 1+norm(b);
   normC2 = 1+ops(C,'norm');
   normX0 = 1+ops(X0,'norm');
   normZ0 = 1+ops(Z0,'norm');
   if (scale_data)
      [At,C,b,normA,normC,normb,X,y,Z] = scaling(blk,At,C,b,X,y,Z);
   else
      normA = 1; normC = 1; normb = 1;
   end
%%
%%-----------------------------------------
%% find the combined list of non-zero
%% elements of Aj, j = 1:k, for each k.
%% IMPORTANT NOTE: Ak, C are permuted.
%%-----------------------------------------
%%
   par.numcolAt = length(b);
   [At,C,X,Z,par.permA,par.permZ] = sortA(blk,At,C,b,X,Z);
   [par.isspA,par.nzlistA,par.nzlistAsum,par.isspAy,par.nzlistAy] = nzlist(blk,At,par);
%%
%%-----------------------------------------
%% create an artifical non-negative block
%% for a purely log-barrier problem
%%-----------------------------------------
%%
   numblkold = size(blk,1);
   nn = 0;
   for p = 1:size(blk,1);
      pblk = blk(p,:);
      idx = find(parbarrier{p}==0);
      if ~isempty(idx);
         if strcmp(pblk{1},'l')
            nn = nn + length(idx);
         elseif strcmp(pblk{1},'q')
            nn = nn + sum(pblk{2}(idx));
         elseif strcmp(pblk{1},'s')
            nn = nn + sum(pblk{2}(idx));
         end
      end
   end
   if (nn==0)
      analytic_prob = 1;
      numblk = size(blk,1)+1;
      blk{numblk,1} = 'l'; blk{numblk,2} = 1;
      At{numblk,1} = sparse(1,length(b));
      C{numblk,1} = 1;
      X{numblk,1} = 1e3;
      Z{numblk,1} = 1e3;
      parbarrier{numblk,1} = 0;
      u2lblk(numblk,1) = 0;
      nn = nn + 1;
   else
      analytic_prob = 0;
   end
%%
   exist_analytic_term = 0;
   for p = 1:size(blk,1);
      idx = find(parbarrier{p} > 0);
      if ~isempty(idx);
         exist_analytic_term = 1;
      end
   end
%%-----------------------------------------
%% initialization
%%-----------------------------------------
%%
   EE = ops(blk,'identity');
   normE2 = ops(EE,'norm'); Zpertold = 1;
   for p = 1:size(blk,1)
      normCC(p) = 1+ops(C(p),'norm');
      normEE(p) = 1+ops(EE(p),'norm');
   end
   [Xchol,indef(1)] = blkcholfun(blk,X);
   [Zchol,indef(2)] = blkcholfun(blk,Z);
   if any(indef)
      msg = 'stop: X or Z not positive definite';
      if (printlevel); fprintf('\n  %s\n',msg); end
      info.termcode = -3;
      info.msg1 = msg;
      obj = []; X = cell(size(blk,1),1); y = []; Z = cell(size(blk,1),1);
      runhist = [];
      return;
   end
   AX = AXfun(blk,At,par.permA,X);
   rp = b-AX;
   ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,y));
   ZpATynorm = ops(ZpATy,'norm');
   Rd = ops(C,'-',ZpATy);
   objadd0 = 0;
   if (scale_data)
      for p = 1:size(blk,1)
         pblk = blk(p,:);
         objadd0 = objadd0 + sum(parbarrier{p}.*pblk{2})*log(normA{p});
      end
   end
   objadd = blkbarrier(blk,X,Z,Xchol,Zchol,parbarrier) + objadd0;
   obj = (normb*normC)*[blktrace(blk,C,X), b'*y] + objadd;
   gap = (normb*normC)*blktrace(blk,X,Z) - diff(objadd);
   relgap = gap/(1+sum(abs(obj)));
   prim_infeas = norm(rp)/normb2;
   dual_infeas = ops(Rd,'norm')/normC2;
   infeas = max(prim_infeas,dual_infeas);
   if (scale_data)
      infeas_org(1) = prim_infeas*normb;
      infeas_org(2) = dual_infeas*normC;
   else
      infeas_org = [0,0];
   end
   trXZ = blktrace(blk,X,Z,parbarrier);
   if (nn > 0); mu  = trXZ/nn; else; mu = gap/ops(X,'getM'); end
   normX = ops(X,'norm');
%%
   termcode = 0; restart = 0;
   pstep = 1; dstep = 1; pred_convg_rate = 1; corr_convg_rate = 1;
   prim_infeas_min  = prim_infeas;
   dual_infeas_min  = dual_infeas;
   prim_infeas_best = prim_infeas;
   dual_infeas_best = dual_infeas;
   infeas_best = infeas;
   relgap_best = relgap;
   homRd = inf; homrp = inf; dy = zeros(length(b),1);
   msg = []; msg2 = []; msg3 = [];
   runhist.pobj    = obj(1);
   runhist.dobj    = obj(2);
   runhist.gap     = gap;
   runhist.relgap  = relgap;
   runhist.pinfeas = prim_infeas;
   runhist.dinfeas = dual_infeas;
   runhist.infeas  = infeas;
   runhist.pstep   = 0;
   runhist.dstep   = 0;
   runhist.step    = 0;
   runhist.normX   = normX;
   runhist.cputime = etime(clock,tstart);
   ttime.preproc   = runhist.cputime;
   ttime.pred = 0; ttime.pred_pstep = 0; ttime.pred_dstep = 0;
   ttime.corr = 0; ttime.corr_pstep = 0; ttime.corr_dstep = 0;
   ttime.pchol = 0; ttime.dchol = 0; ttime.misc = 0;
%%
%%-----------------------------------------
%% display parameters and initial info
%%-----------------------------------------
%%
   if (printlevel >= 2)
      fprintf('\n********************************************');
      fprintf('***********************\n');
      fprintf('   SDPT3: Infeasible path-following algorithms');
      fprintf('\n********************************************');
      fprintf('***********************\n');
      [hh,mm,ss] = mytime(ttime.preproc);
      if (printlevel>=3)
         fprintf(' version  predcorr  gam  expon  scale_data\n');
         if (vers == 1); fprintf('   HKM '); elseif (vers == 2); fprintf('    NT '); end
         fprintf('     %1.0f      %4.3f',predcorr,gam);
         fprintf('   %1.0f        %1.0f    %1.0f\n',expon,scale_data);
         fprintf('\nit pstep dstep pinfeas dinfeas  gap')
         fprintf('      prim-obj      dual-obj    cputime\n');
         fprintf('------------------------------------------------');
         fprintf('-------------------\n');
         fprintf('%2.0f|%4.3f|%4.3f|%2.1e|%2.1e|',0,0,0,prim_infeas,dual_infeas);
         fprintf('%2.1e|%- 7.6e %- 7.6e| %s:%s:%s|',gap,obj(1),obj(2),hh,mm,ss);
      end
   end
%%
%%---------------------------------------------------------------
%% start main loop
%%---------------------------------------------------------------
%%
   param.termcode    = termcode;
   param.iter        = 0;
   param.obj         = obj;
   param.relgap      = relgap;
   param.prim_infeas = prim_infeas;   param.dual_infeas = dual_infeas;
   param.homRd       = homRd;         param.homrp       = homrp;
   param.AX          = AX;            param.ZpATynorm   = ZpATynorm;
   param.normA       = normA;
   param.normb       = normb;         param.normC       = normC;
   param.normX0      = normX0;        param.normZ0      = normZ0;
   param.m0          = m0;            param.indeprows   = indeprows;
   param.prim_infeas_bad = 0;
   param.dual_infeas_bad = 0;
   param.prim_infeas_min = prim_infeas;
   param.dual_infeas_min = dual_infeas;
   param.gaptol      = gaptol;
   param.inftol      = inftol;
   param.maxit       = maxit;
   param.scale_data  = scale_data;
   param.printlevel  = printlevel;
   param.ublksize    = ublksize;
   Xbest = X; ybest = y; Zbest = Z;
%%
   for iter = 1:maxit;
      tstart  = clock;
      timeold = tstart;
      update_iter = 0; breakyes = 0;
      pred_slow = 0; corr_slow = 0; step_short = 0;
      par.parbarrier = parbarrier;
      par.iter    = iter;
      par.obj     = obj;
      par.relgap  = relgap;
      par.pinfeas = prim_infeas;
      par.dinfeas = dual_infeas;
      par.rp      = rp;
      par.y       = y;
      par.dy      = dy;
      par.normX   = normX;
      par.ZpATynorm = ZpATynorm;
      %%if (printlevel > 2); fprintf(' %2.1e',par.normX); end
      if (iter == 1 | restart); Cpert = min(1,normC2/ops(EE,'norm')); end
      if (runhist.dinfeas(1) > 1e-3) & (~exist_analytic_term) ...
         & (relgap > 1e-4)
         if (par.normX > 5e3 & iter < 20)
            Cpert = Cpert*0.5;
         elseif (par.normX > 5e2 & iter < 20);
            Cpert = Cpert*0.3;
         else;
            Cpert = Cpert*0.1;
         end
         Rd = ops(Rd,'+',EE,Cpert);
         %%if (printlevel > 2); fprintf('|%2.1e',Cpert); end
      end
%%---------------------------------------------------------------
%% predictor step.
%%---------------------------------------------------------------
%%
      if (predcorr)
         sigma = 0;
      else
         sigma = 1-0.9*min(pstep,dstep);
         if (iter == 1); sigma = 0.5; end;
      end
      sigmu = cell(size(blk,1),1);
      for p = 1:size(blk,1)
         sigmu{p} = max(sigma*mu, parbarrier{p}');
      end
      invXchol = cell(size(blk,1),1);
      invZchol = ops(Zchol,'inv');
      if (vers == 1);
         [par,dX,dy,dZ,coeff,L,hRd] = ...
          HKMpred(blk,At,par,rp,Rd,sigmu,X,Z,invZchol);
      elseif (vers == 2);
         [par,dX,dy,dZ,coeff,L,hRd] = ...
          NTpred(blk,At,par,rp,Rd,sigmu,X,Z,Zchol,invZchol);
      end
      if (solve_ok <= 0)
         msg = 'stop: difficulty in computing predictor directions';
         if (printlevel); fprintf('\n  %s',msg); end
         runhist.pinfeas(iter+1) = runhist.pinfeas(iter);
         runhist.dinfeas(iter+1) = runhist.dinfeas(iter);
         runhist.relgap(iter+1)  = runhist.relgap(iter);
         runhist.cputime(iter+1) = etime(clock,tstart);
         termcode = -4;
         break; %% do not ues breakyes = 1
      end
      timenew = clock;
      ttime.pred = ttime.pred + etime(timenew,timeold); timeold = timenew;
%%
%%-----------------------------------------
%% step-lengths for predictor step
%%-----------------------------------------
%%
      if (gam == 0)
         gamused = 0.9 + 0.09*min(pstep,dstep);
      else
         gamused = gam;
      end
      [Xstep,invXchol] = steplength(blk,X,dX,Xchol,invXchol);
      pstep = min(1,gamused*full(Xstep));
      timenew = clock;
      ttime.pred_pstep = ttime.pred_pstep + etime(timenew,timeold); timeold = timenew;
      Zstep = steplength(blk,Z,dZ,Zchol,invZchol);
      dstep = min(1,gamused*full(Zstep));
      trXZnew = trXZ + pstep*blktrace(blk,dX,Z,parbarrier) ...
                 + dstep*blktrace(blk,X,dZ,parbarrier) ...
                 + pstep*dstep*blktrace(blk,dX,dZ,parbarrier);
      if (nn > 0); mupred  = trXZnew/nn; else; mupred = 1e-16; end
      mupredhist(iter) = mupred;
      timenew = clock;
      ttime.pred_dstep = ttime.pred_dstep + etime(timenew,timeold); timeold = timenew;
%%
%%-----------------------------------------
%%  stopping criteria for predictor step.
%%-----------------------------------------
%%
      if (min(pstep,dstep) < steptol) & (stoplevel) & (iter > 10)
         msg = 'stop: steps in predictor too short';
         if (printlevel)
            fprintf('\n  %s',msg);
            fprintf(': pstep = %3.2e,  dstep = %3.2e\n',pstep,dstep);
         end
         runhist.cputime(iter+1) = etime(clock,tstart);
         termcode = -2;
         breakyes = 1;
      end
      if (~predcorr)
         if (iter >= 2)
            idx = [max(2,iter-2) : iter];
            pred_slow = all(mupredhist(idx)./mupredhist(idx-1) > 0.4);
            idx = [max(2,iter-5) : iter];
            pred_convg_rate = mean(mupredhist(idx)./mupredhist(idx-1));
            pred_slow = pred_slow + (mupred/mu > 5*pred_convg_rate);
         end
         if (max(mu,infeas) < 1e-6) & (pred_slow) & (stoplevel)
            msg = 'stop: lack of progress in predictor';
            if (printlevel)
               fprintf('\n  %s',msg);
               fprintf(': mupred/mu = %3.2f, pred_convg_rate = %3.2f.',...
               mupred/mu,pred_convg_rate);
            end
            runhist.cputime(iter+1) = etime(clock,tstart);
            termcode = -2;
            breakyes = 1;
         else
            update_iter = 1;
         end
      end
%%---------------------------------------------------------------
%% corrector step.
%%---------------------------------------------------------------
%%
      if (predcorr) & (~breakyes)
         step_pred = min(pstep,dstep);
         if (mu > 1e-6)
            if (step_pred < 1/sqrt(3));
               expon_used = 1;
            else
               expon_used = max(expon,3*step_pred^2);
            end
         else
            expon_used = max(1,min(expon,3*step_pred^2));
         end
         if (nn==0)
             sigma = 0.2;
         elseif (mupred < 0)
             sigma = 0.8;
         else
            sigma = min(1, (mupred/mu)^expon_used);
         end
         sigmu = cell(size(blk,1),1);
         for p = 1:size(blk,1)
            sigmu{p} = max(sigma*mu, parbarrier{p}');
         end
         if (vers == 1)
            [dX,dy,dZ] = HKMcorr(blk,At,par,rp,Rd,sigmu,hRd,...
             dX,dZ,coeff,L,X,Z);
         elseif (vers == 2)
            [dX,dy,dZ] = NTcorr(blk,At,par,rp,Rd,sigmu,hRd,...
             dX,dZ,coeff,L,X,Z);
         end
         if (solve_ok <= 0)
            msg = 'stop: difficulty in computing corrector directions';
            if (printlevel); fprintf('\n  %s',msg); end
            runhist.pinfeas(iter+1) = runhist.pinfeas(iter);
            runhist.dinfeas(iter+1) = runhist.dinfeas(iter);
            runhist.relgap(iter+1)  = runhist.relgap(iter);
            runhist.cputime(iter+1) = etime(clock,tstart);
            termcode = -4;
            break; %% do not ues breakyes = 1
         end
         timenew = clock;
         ttime.corr = ttime.corr + etime(timenew,timeold); timeold = timenew;
%%
%%-----------------------------------
%% step-lengths for corrector step
%%-----------------------------------
%%
         if (gam == 0)
            gamused = 0.9 + 0.09*min(pstep,dstep);
         else
            gamused = gam;
         end
         Xstep = steplength(blk,X,dX,Xchol,invXchol);
         pstep = min(1,gamused*full(Xstep));
         timenew = clock;
         ttime.corr_pstep = ttime.corr_pstep+etime(timenew,timeold); timeold = timenew;
         Zstep = steplength(blk,Z,dZ,Zchol,invZchol);
         dstep = min(1,gamused*full(Zstep));
         trXZnew = trXZ + pstep*blktrace(blk,dX,Z,parbarrier) ...
                    + dstep*blktrace(blk,X,dZ,parbarrier)...
                    + pstep*dstep*blktrace(blk,dX,dZ,parbarrier);
         if (nn > 0); mucorr  = trXZnew/nn; else; mucorr = 1e-16; end
         timenew = clock;
         ttime.corr_dstep = ttime.corr_dstep+etime(timenew,timeold); timeold = timenew;
%%
%%-----------------------------------------
%%  stopping criteria for corrector step
%%-----------------------------------------
         if (iter >= 2)
            idx = [max(2,iter-2) : iter];
            corr_slow = all(runhist.gap(idx)./runhist.gap(idx-1) > 0.8);
            idx = [max(2,iter-5) : iter];
            corr_convg_rate = mean(runhist.gap(idx)./runhist.gap(idx-1));
            corr_slow = corr_slow + (mucorr/mu > max(min(1,5*corr_convg_rate),0.8));
         end
	 if (max(relgap,infeas) < 1e-6) & (iter > 20) ...
            & (corr_slow > 1) & (stoplevel)
            msg = 'stop: lack of progress in corrector';
   	    if (printlevel)
               fprintf('\n  %s',msg);
               fprintf(': mucorr/mu = %3.2f, corr_convg_rate = %3.2f',...
               mucorr/mu,corr_convg_rate);
            end
            runhist.cputime(iter+1) = etime(clock,tstart);
            termcode = -2;
            breakyes = 1;
         else
            update_iter = 1;
         end
      end
%%---------------------------------------------------------------
%% udpate iterate
%%---------------------------------------------------------------
      indef = [1,1];
      if (update_iter)
         for t = 1:5
            [Xchol,indef(1)] = blkcholfun(blk,ops(X,'+',dX,pstep));
            timenew = clock;
            ttime.pchol = ttime.pchol + etime(timenew,timeold); timeold = timenew;
            if (indef(1)); pstep = 0.8*pstep; else; break; end
         end
	 if (t > 1); pstep = gamused*pstep; end
	 for t = 1:5
            [Zchol,indef(2)] = blkcholfun(blk,ops(Z,'+',dZ,dstep));
            timenew = clock;
            ttime.dchol = ttime.dchol + etime(timenew,timeold); timeold = timenew;
            if (indef(2)); dstep = 0.8*dstep; else; break; end
         end
	 if (t > 1); dstep = gamused*dstep; end
         %%-------------------------------------------
         AXtmp = AX + pstep*AXfun(blk,At,par.permA,dX);
         prim_infeasnew = norm(b-AXtmp)/normb2;
         if (relgap < 5*infeas); alpha = 1e2; else; alpha = 1e3; end
         if any(indef)
            if indef(1); msg = 'stop: X not positive definite'; end
            if indef(2); msg = 'stop: Z not positive definite'; end
            if (printlevel); fprintf('\n  %s',msg); end
            termcode = -3;
            breakyes = 1;
         elseif (prim_infeasnew > max([1e-8,relgap,20*prim_infeas]) & iter > 10) ...
            | (prim_infeasnew > max([1e-7,1e3*prim_infeas,0.1*relgap]) & relgap < 1e-2) ...
            | (prim_infeasnew > alpha*max([1e-9,param.prim_infeas_min]) ...
               & (prim_infeasnew > max([3*prim_infeas,0.1*relgap])) ...
               & (iter > 25) & (dual_infeas < 1e-6) & (relgap < 0.1)) ...
            | ((prim_infeasnew > 1e3*prim_infeas & prim_infeasnew > 1e-12) ...
               & (max(relgap,dual_infeas) < 1e-8))
            if (stoplevel)
               msg = 'stop: primal infeas has deteriorated too much';
               if (printlevel); fprintf('\n  %s, %2.1e',msg,prim_infeasnew); end
               termcode = -7;
               breakyes = 1;
            end
         elseif (trXZnew > 1.05*runhist.gap(iter)) & (~exist_analytic_term) ...
	    & ((infeas < 1e-5) & (relgap < 1e-4) & (iter > 20) ...
	       | (max(infeas,relgap) < 1e-7) & (iter > 10))
            if (stoplevel)
               msg = 'stop: progress in duality gap has deteriorated';
               if (printlevel); fprintf('\n  %s, %2.1e',msg,trXZnew); end
               termcode = -8;
               breakyes = 1;
            end
         else
            X = ops(X,'+',dX,pstep);
            y = y + dstep*dy;
            Z = ops(Z,'+',dZ,dstep);
         end
      end
%%---------------------------------------------------------------
%% adjust linear blk arising from unrestricted blk
%%---------------------------------------------------------------
      if (~breakyes)
         for p = 1:size(blk,1)
            if (u2lblk(p) == 1)
               len = blk{p,2}/2;
               xtmp = min(X{p}([1:len]),X{p}(len+[1:len]));
               alpha = 0.8;
               X{p}([1:len])     = X{p}([1:len]) - alpha*xtmp;
               X{p}(len+[1:len]) = X{p}(len+[1:len]) - alpha*xtmp;
               if (mu < 1e-4) %% old: (mu < 1e-7)
                  Z{p} = 0.5*mu./max(1,X{p}); %% good to keep this step
               else
                  ztmp = min(1,max(Z{p}([1:len]),Z{p}(len+[1:len])));
                  if (dual_infeas > 1e-4 & dstep < 0.2)
                     beta = 0.3;
                  else
                     beta = 0.0;
                  end
                  %% important to set beta = 0 at later stage.
                  Z{p}([1:len])     = Z{p}([1:len]) + beta*ztmp;
                  Z{p}(len+[1:len]) = Z{p}(len+[1:len]) + beta*ztmp;
               end
            end
         end
      end
%%--------------------------------------------------
%% perturb Z: do this step before checking for break
%%--------------------------------------------------
      if (~breakyes) & (~exist_analytic_term)
         trXZtmp = blktrace(blk,X,Z);
         trXE  = blktrace(blk,X,EE);
         Zpert = max(1e-12,0.2*min(relgap,prim_infeas)).*normC2./normE2;
         Zpert = min(Zpert,0.1*trXZtmp./trXE);
         Zpert = min([1,Zpert,1.5*Zpertold]);
         if (infeas < 0.1)
            Z = ops(Z,'+',EE,Zpert);
            [Zchol,indef(2)] = blkcholfun(blk,Z);
            if any(indef(2))
               msg = 'stop: Z not positive definite';
               if (printlevel); fprintf('\n  %s',msg); end
               termcode = -3;
               breakyes = 1;
            end
            %%if (printlevel > 2); fprintf(' %2.1e',Zpert); end
         end
         Zpertold = Zpert;
      end
%%---------------------------------------------------------------
%% compute rp, Rd, infeasibities, etc
%%---------------------------------------------------------------
%%
      AX  = AXfun(blk,At,par.permA,X);
      rp  = b-AX;
      ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,y));
      ZpATynorm = ops(ZpATy,'norm');
      Rd  = ops(C,'-',ZpATy);
      objadd = blkbarrier(blk,X,Z,Xchol,Zchol,parbarrier) + objadd0;
      obj = (normb*normC)*[blktrace(blk,C,X), b'*y] + objadd;
      gap = (normb*normC)*blktrace(blk,X,Z) - diff(objadd);
      relgap = gap/(1+sum(abs(obj)));
      prim_infeas = norm(rp)/normb2;
      dual_infeas = ops(Rd,'norm')/normC2;
      infeas = max(prim_infeas,dual_infeas);
      if (scale_data)
         infeas_org(1) = prim_infeas*normb;
         infeas_org(2) = dual_infeas*normC;
      end
      homRd = inf; homrp = inf;
      if (ops(parbarrier,'norm') == 0)
         if (obj(2) > 0); homRd = ZpATynorm/(obj(2)); end
         if (obj(1) < 0); homrp = norm(AX)/(-obj(1))/(normC); end
      end
      trXZ = blktrace(blk,X,Z,parbarrier);
      if (nn > 0); mu = trXZ/nn; else; mu = gap/ops(X,'getM'); end
      normX = ops(X,'norm');
%%
      runhist.pobj(iter+1)  = obj(1);
      runhist.dobj(iter+1)  = obj(2);
      runhist.gap(iter+1)   = gap;
      runhist.relgap(iter+1)  = relgap;
      runhist.pinfeas(iter+1) = prim_infeas;
      runhist.dinfeas(iter+1) = dual_infeas;
      runhist.infeas(iter+1)  = infeas;
      runhist.pstep(iter+1)   = pstep;
      runhist.dstep(iter+1)   = dstep;
      runhist.step(iter+1)    = min(pstep,dstep);
      runhist.normX(iter+1)   = normX;
      runhist.cputime(iter+1) = etime(clock,tstart);
      timenew = clock;
      ttime.misc = ttime.misc + etime(timenew,timeold); timeold = timenew;
      [hh,mm,ss] = mytime(sum(runhist.cputime));
      if (printlevel>=3)
         fprintf('\n%2.0f|%4.3f|%4.3f',iter,pstep,dstep);
         fprintf('|%2.1e|%2.1e|%2.1e|',prim_infeas,dual_infeas,gap);
         fprintf('%- 7.6e %- 7.6e| %s:%s:%s|',obj(1),obj(2),hh,mm,ss);
      end
%%--------------------------------------------------
%% check convergence
%%--------------------------------------------------
      param.use_LU      = use_LU;
      param.stoplevel   = stoplevel;
      param.termcode    = termcode;
      param.iter        = iter;
      param.obj         = obj;
      param.gap         = gap;
      param.relgap      = relgap;
      param.prim_infeas = prim_infeas;
      param.dual_infeas = dual_infeas;
      param.mu        = mu;
      param.homRd     = homRd;
      param.homrp     = homrp;
      param.AX        = AX;
      param.ZpATynorm = ZpATynorm;
      param.normX     = ops(X,'norm');
      param.normZ     = ops(Z,'norm');
      param.numpertdiagschur = numpertdiagschur;
      if (~breakyes)
         [param,breakyes,restart,msg2] = sqlpcheckconvg(param,runhist);
      end
      if (restart)
         [X,y,Z] = infeaspt(blk,At,C,b,2,1e5);
         rp  = b-AXfun(blk,At,par.permA,X);
         ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,y));
         Rd  = ops(C,'-',ZpATy);
         trXZ = blktrace(blk,X,Z,parbarrier);
         mu   = trXZ/nn;
         gap  =  (normb*normC)*blktrace(blk,X,Z) - diff(objadd);
         prim_infeas = norm(rp)/normb2;
         dual_infeas = ops(Rd,'norm')/normC2;
         infeas = max(prim_infeas,dual_infeas);
         [Xchol,indef(1)] = blkcholfun(blk,X);
         [Zchol,indef(2)] = blkcholfun(blk,Z);
         stoplevel = 3;
      end
%%--------------------------------------------------
%% check for break
%%--------------------------------------------------
      if ((prim_infeas < 1.5*prim_infeas_best) ...
         | (max(relgap,infeas) < 0.8*max(relgap_best,infeas_best))) ...
         & (max(relgap,dual_infeas) < 0.8*max(relgap_best,dual_infeas_best))
         Xbest = X; ybest = y; Zbest = Z;
         prim_infeas_best = prim_infeas;
         dual_infeas_best = dual_infeas;
         relgap_best = relgap; infeas_best = infeas;
         update_best(iter+1) = 1;
         %%fprintf('#')
      else
         update_best(iter+1) = 0;
      end
      if (max(relgap_best,infeas_best) < 1e-4 ...
          & norm(update_best(max(1,iter-1):iter+1)) == 0)
         msg = 'lack of progress in infeas';
         if (printlevel); fprintf('\n  %s',msg); end
         termcode = -9;
         breakyes = 1;
      end
      if (breakyes); break; end
   end
%%---------------------------------------------------------------
%% end of main loop
%%---------------------------------------------------------------
%%
   use_bestiter = 1;
   if (use_bestiter) & (param.termcode <= 0)
      X = Xbest; y = ybest; Z = Zbest;
      Xchol = blkcholfun(blk,X);
      Zchol = blkcholfun(blk,Z);
      AX = AXfun(blk,At,par.permA,X);
      rp = b-AX;
      ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,y));
      Rd = ops(C,'-',ZpATy);
      objadd = blkbarrier(blk,X,Z,Xchol,Zchol,parbarrier) + objadd0;
      obj = (normb*normC)*[blktrace(blk,C,X), b'*y] + objadd;
      gap = (normb*normC)*blktrace(blk,X,Z) - diff(objadd);
      relgap = gap/(1+sum(abs(obj)));
      prim_infeas = norm(rp)/normb2;
      dual_infeas = ops(Rd,'norm')/normC2;
      infeas = max(prim_infeas,dual_infeas);
      runhist.pobj(iter+1)  = obj(1);
      runhist.dobj(iter+1)  = obj(2);
      runhist.gap(iter+1)   = gap;
      runhist.relgap(iter+1)  = relgap;
      runhist.pinfeas(iter+1) = prim_infeas;
      runhist.dinfeas(iter+1) = dual_infeas;
      runhist.infeas(iter+1)  = infeas;
   end
%%---------------------------------------------------------------
%% unscale and produce infeasibility certificates if appropriate
%%---------------------------------------------------------------
   if (iter >= 1)
      [X,y,Z,termcode,resid,reldist,msg3] = ...
      sqlpmisc(blk,At,C,b,X,y,Z,par.permZ,param);
   end
%%---------------------------------------------------------------
%% recover unrestricted blk from linear blk
%%---------------------------------------------------------------
%%
   for p = 1:size(blk,1)
      if (u2lblk(p) == 1)
         n = blk{p,2}/2;
         X{p} = X{p}(1:n)-X{p}(n+[1:n]);
         Z{p} = Z{p}(1:n);
      end
   end
   for p = 1:size(ublkidx,1)
      if ~isempty(ublkidx{p,2})
         n0 = ublkidx{p,1}; idxB = setdiff([1:n0]',ublkidx{p,2});
         tmp = zeros(n0,1); tmp(idxB) = X{p}; X{p} = tmp;
         tmp = zeros(n0,1); tmp(idxB) = Z{p}; Z{p} = tmp;
      end
   end
   if (analytic_prob)
      X = X(1:numblkold); Z = Z(1:numblkold);
   end
%%---------------------------------------------------------------
%% print summary
%%---------------------------------------------------------------
%%
   maxC = 1+ops(ops(C,'abs'),'max');
   maxb = 1+max(abs(b));
   if (scale_data)
      dimacs = [infeas_org(1)*normb2/maxb; 0; infeas_org(2)*normC2/maxC; 0];
   else
      dimacs = [prim_infeas*normb2/maxb; 0; dual_infeas*normC2/maxC; 0];
   end
   dimacs = [dimacs; [-diff(obj); gap]/(1+sum(abs(obj)))];
   info.dimacs   = dimacs;
   info.termcode = termcode;
   info.iter     = iter;
   info.obj      = obj;
   info.gap      = gap;
   info.relgap   = relgap;
   info.pinfeas  = prim_infeas;
   info.dinfeas  = dual_infeas;
   info.cputime  = sum(runhist.cputime);
   info.time     = ttime;
   info.resid    = resid;
   info.reldist  = reldist;
   info.normX    = ops(X,'norm');
   info.normy    = norm(y);
   info.normZ    = ops(Z,'norm');
   info.normb    = normb2; info.maxb = maxb;
   info.normC    = normC2; info.maxC = maxC;
   info.normA    = normA2;
   info.msg1     = msg;
   info.msg2     = msg2;
   info.msg3     = msg3;
   sqlpsummary(info,ttime,infeas_org,printlevel);
%%*****************************************************************************

author: yu-jiang

repo: radpbook

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


text:
function heated_plate ( epsilon, output_filename )

%*****************************************************************************80
%
%  Purpose:
%
%    MAIN is the main program for HEATED_PLATE.
%
%  Discussion:
%
%    This code solves the steady state heat equation on a rectangular region.
%
%    The sequential version of this program needs approximately
%    18/eps iterations to complete.
%
%
%    The physical region, and the boundary conditions, are suggested
%    by this diagram;
%
%                   W = 0
%             +------------------+
%             |                  |
%    W = 100  |                  | W = 100
%             |                  |
%             +------------------+
%                   W = 100
%
%    The region is covered with a grid of M by N nodes, and an N by N
%    array W is used to record the temperature.  The correspondence between
%    array indices and locations in the region is suggested by giving the
%    indices of the four corners:
%
%                  I = 0
%          [0][0]-------------[0][N-1]
%             |                  |
%      J = 0  |                  |  J = N-1
%             |                  |
%        [M-1][0]-----------[M-1][N-1]
%                  I = M-1
%
%    The steady state solution to the discrete heat equation satisfies the
%    following condition at an interior grid point:
%
%      W[Central] = (1/4) * ( W[North] + W[South] + W[East] + W[West] )
%
%    where "Central" is the index of the grid point, "North" is the index
%    of its immediate neighbor to the "north", and so on.
%
%    Given an approximate solution of the steady state heat equation, a
%    "better" solution is given by replacing each interior point by the
%    average of its 4 neighbors - in other words, by using the condition
%    as an ASSIGNMENT statement:
%
%      W[Central]  <=  (1/4) * ( W[North] + W[South] + W[East] + W[West] )
%
%    If this process is repeated often enough, the difference between successive
%    estimates of the solution will go to zero.
%
%    This program carries out such an iteration, using a tolerance specified by
%    the user, and writes the final estimate of the solution to a file that can
%    be used for graphic processing.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 July 2010
%
%  Author:
%
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Michael Quinn,
%    Parallel Programming in C with MPI and OpenMP,
%    McGraw-Hill, 2004,
%    ISBN13: 978-0071232654,
%    LC: QA76.73.C15.Q55.
%
%  Parameters:
%
%    Commandline argument 1, real EPSILON, the error tolerance.
%
%    Commandline argument 2, string OUTPUT_FILENAME, the name of the file into which
%    the steady state solution is written when the program has completed.
%
%  Local parameters:
%
%    Local, real DIFF, the norm of the change in the solution from
%    one iteration to the next.
%
%    Local, real MEAN, the average of the boundary values, used
%    to initialize the values of the solution in the interior.
%
%    Local, real U(M,N), the solution at the previous iteration.
%
%    Local, real W(M,N), the solution computed at the latest
%    iteration.
%
  m = 500;
  n = 500;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HEATED_PLATE\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  A program to solve for the steady state temperature distribution\n' );
  fprintf ( 1, '  over a rectangular plate.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial grid of %d by %d points.\n', m, n );
%
%  Read EPSILON from the command line or the user.
%
  if ( nargin < 1 )
    epsilon = input ( '  Enter EPSILON, the error tolerance.' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The iteration will repeat until the change is <= %e\n', epsilon );
%
%  Read OUTPUT_FILENAME from the command line or the user.
%
  if ( nargin < 2 )
    output_filename = input ( ...
      '  Enter OUTPUT_FILENAME, the name of the output file.' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The steady state solution will be written to "%s".\n', ...
    output_filename );
%
%  Set the boundary values, which don't change.
%
  w(2:m-1,1) = 100.0;
  w(2:m-1,n) = 100.0;
  w(m,1:n) = 100.0;
  w(1,1:n) =   0.0;
%
%  Average the boundary values, to come up with a reasonable
%  initial value for the interior.
%
  mean = ( ...
      sum ( w(2:m-1,1) ) ...
    + sum ( w(2:m-1,n) ) ...
    + sum ( w(m,1:n)   ) ...
    + sum ( w(1,1:n)   ) ) ...
    / ( 2 * m + 2 * n - 4 );
%
%  Initialize the interior solution to the mean value.
%
  w(2:m-1,2:n-1) = mean;
%
%  iterate until the  new solution W differs from the old solution U
%  by no more than EPSILON.
%
  iterations = 0;
  iterations_print = 1;
  diff = epsilon;

  fprintf ( 1, '\n' );
  fprintf ( 1, ' Iteration  Change\n' );
  fprintf ( 1, '\n' );

  tic;

  while ( epsilon <= diff )

    u(1:m,1:n) = w(1:m,1:n);

    w(2:m-1,2:n-1) = 0.25 * ( ...
        u(1:m-2,2:n-1) ...
      + u(3:m,2:n-1) ...
      + u(2:m-1,1:n-2) ...
      + u(2:m-1,3:n) );

    diff = max ( max ( abs ( u(1:m,1:n) - w(1:m,1:n) ) ) );

    iterations = iterations + 1;

    if ( iterations == iterations_print )
      fprintf ( 1, '  %8d  %14f\n', iterations, diff );
      iterations_print = 2 * iterations_print;
    end

  end

  wtime = toc;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  %8d  %14f\n', iterations, diff );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Error tolerance achieved.\n' );
  fprintf ( 1, '  Wallclock time = %f\n', wtime );
%
%  Write the solution to the output file.
%
  output_unit = fopen ( output_filename, 'wt' );

  fprintf ( output_unit, '%d\n', m );
  fprintf ( output_unit, '%d\n', n );
  for i = 1 : m
    for j = 1 : n
      fprintf (output_unit, '  %f', w(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end

  fclose ( output_unit );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solution written to the output file "%s".\n', output_filename );
%
%  Terminate.
%

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HEATED_PLATE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end


text:
function F = asech(F, varargin)
%ASECH   Inverse hyperbolic secant of a CHEBFUN.
%   ASECH(F) computes the inverse hyperbolic secant of the CHEBFUN F.
%
%   ASECH(F, PREF) does the same but uses the CHEBFUNPREF object PREF when
%   computing the composition.
%
% See also SECH.

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

% Call the compose method:
F = compose(F, @asech, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function out = angles(n)
%ANGLES   Return the angles of the Chebyshev points of 1st kind in [-1, 1].
%   CHEBTECH1.ANGLES(N) returns ACOS(X), where X are the N Chebyshev points of
%   the 1st kind in [-1, 1].
%
% See also POINTS, CHEBPTS, LENGTH.
%
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

out = (n-.5:-1:.5).'*pi/n;

end


text:
function linpack_z_test35 ( )

%*****************************************************************************80
%
%% TEST35 tests ZTRCO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST35\n' );
  fprintf ( 1, '  For a double precision complex (C)\n' );
  fprintf ( 1, '  triangular matrix (TR),\n' );
  fprintf ( 1, '  ZTRCO estimates the condition.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix order is N = %d\n', n );
%
%  Set the matrix.
%
  seed = 123456789;

  for i = 1 : n
    for j = 1 : i
      [ a(i,j), seed ] = c8_uniform_01 ( seed );
    end
    a(i,i+1:n) = 0.0;
  end
%
%  Get the condition of the lower triangular matrix.
%
  job = 0;
  rcond = ztrco ( a, lda, n, job );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimated reciprocal condition RCOND = %f\n', rcond );

  return
end


text:
function mCol=grColEdge(E)
% function mCol=grColEdge(E) solve the color graph problem
% for edges of the graph.
% Input parameter:
%   E(m,2) - the edges of graph;
%     1st and 2nd elements of each row is numbers of vertexes;
%     m - number of edges.
% Output parameter:
%   mCol(m,1) - the list of the colors of edges.
% Uses the sequential deleting of the maximal matching sets.
% Required the Optimization Toolbox v.3.0.1 or over.
% Author: Sergiy Iglin
% e-mail: siglin@yandex.ru
% personal page: http://iglin.exponenta.ru

% ============= Input data validation ==================
if nargin<1,
  error('There are no input data!')
end
[m,n,E] = grValidation(E); % E data validation
E=[E(:,1:2),[1:m]']; % numbers of vertexes and numbers of edges
mCol=zeros(m,1); % initial value
% ============= Main cycle with MaxMatch deleting ====
while any(mCol==0),
  ne=find(mCol==0); % uncolored edges
  E1=E(ne,:); % it's edges
  nMM=grMaxMatch(E1(:,1:2)); % the maximal matching
  mCol(E1(nMM,3))=max(mCol)+1; % the next colorend
end
return

author: HuangCongQing


% params
totalNum = 1502861;
bValue = 3000;
lamda = 0.1;
qFile = '/home/yinhuan/Data/mapModel/park/weightVector.txt';
visFilesDir = '/home/yinhuan/Data/mapModel/park/visMatrix/';
maxQ = 110;
minQ = 1;

splitLength_0 = 100;
saveResultsDir_0 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/0/';
saveReIndexFile_0 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/index/0.txt';
saveNewQFile_0 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/weight/0.txt';

splitLength_1 = 200;
saveResultsDir_1 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/1/';
saveReIndexFile_1 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/index/1.txt';
saveNewQFile_1 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/weight/1.txt';

splitLength_2 = 400;
saveResultsDir_2 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/2/';
saveReIndexFile_2 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/index/2.txt';
saveNewQFile_2 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/weight/2.txt';

splitLength_3 = 800;
saveResultsDir_3 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/3/';
saveReIndexFile_3 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/index/3.txt';
saveNewQFile_3 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/weight/3.txt';

splitLength_4 = 800;
saveResultsDir_4 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/4/';
saveReIndexFile_4 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/index/4.txt';
saveNewQFile_4 = '/home/yinhuan/Data/mapModel/park/iter_b_3000/weight/4.txt';

%% iteration of optimization

visCells_origin = fromVisDirToCells(visFilesDir);

[ epsilon_soft_0, time_sum_0  ] = Uniform_loopCompress_Cells( lamda, qFile, visCells_origin, splitLength_0, totalNum, bValue, saveResultsDir_0 );
min_cost_0 = get_min_cost(minQ, maxQ, saveResultsDir_0, qFile, lamda, epsilon_soft_0);
[totalNum_0, visCells_0] = deleteZeroPoints_Cells( qFile, visCells_origin, saveResultsDir_0, saveReIndexFile_0, saveNewQFile_0 );

[ epsilon_soft_1, time_sum_1  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_0, visCells_0, splitLength_1, totalNum_0, bValue, saveResultsDir_1 );
min_cost_1 = get_min_cost(minQ, maxQ, saveResultsDir_1, saveNewQFile_0, lamda, epsilon_soft_1);
[totalNum_1, visCells_1] = deleteZeroPoints_Cells( saveNewQFile_0, visCells_0, saveResultsDir_1, saveReIndexFile_1, saveNewQFile_1 );

[ epsilon_soft_2, time_sum_2  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_1, visCells_1, splitLength_2, totalNum_1, bValue, saveResultsDir_2 );
min_cost_2 = get_min_cost(minQ, maxQ, saveResultsDir_2, saveNewQFile_1, lamda, epsilon_soft_2);
[totalNum_2, visCells_2] = deleteZeroPoints_Cells( saveNewQFile_1, visCells_1, saveResultsDir_2, saveReIndexFile_2, saveNewQFile_2 );

[ epsilon_soft_3, time_sum_3  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_2, visCells_2, splitLength_3, totalNum_2, bValue, saveResultsDir_3 );
min_cost_3 = get_min_cost(minQ, maxQ, saveResultsDir_3, saveNewQFile_2, lamda, epsilon_soft_3);
[totalNum_3, visCells_3] = deleteZeroPoints_Cells( saveNewQFile_2, visCells_2, saveResultsDir_3, saveReIndexFile_3, saveNewQFile_3 );

[ epsilon_soft_4, time_sum_4  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_3, visCells_3, splitLength_4, totalNum_3, bValue, saveResultsDir_4 );
min_cost_4 = get_min_cost(minQ, maxQ, saveResultsDir_4, saveNewQFile_3, lamda, epsilon_soft_4);
[totalNum_4, visCells_4] = deleteZeroPoints_Cells( saveNewQFile_3, visCells_3, saveResultsDir_4, saveReIndexFile_4, saveNewQFile_4 );

%% save the results

compressIndex_3 = anal_reindex_last(saveResultsDir_4, saveReIndexFile_3);
compressIndex_2 = anal_reindex_middle(compressIndex_3, saveReIndexFile_2);
compressIndex_1 = anal_reindex_middle(compressIndex_2, saveReIndexFile_1);
compressIndex_0 = anal_reindex_middle(compressIndex_1, saveReIndexFile_0);

dlmwrite('/home/yinhuan/Data/mapModel/park/iter_b_3000/compressIndex_final.txt', compressIndex_0, 'precision', '%d');


text:
function [Rtilt,R] = rectify(XYZ)
    %% XYZ is HxWx3 matrix
    % X = XYZ(:,:,1);Y = XYZ(:,:,2);Z = XYZ(:,:,3);
    % XYZnew = Rtilt*[X(:),Y(:),Z(:)]'
    [Rtilt,R,world_center] = dominantAxes([eye(3) zeros(3,1)],XYZ);

function [Rtilt,R,world_center] = dominantAxes(cameraRt, pts)

XYZ = pts;
S = 1;
points = [reshape(XYZ(:,:,1),1,[]);reshape(XYZ(:,:,2),1,[]);reshape(XYZ(:,:,3),1,[])];
pointsOK = points(:,sum(isnan(points),1)==0);
pointsOK = pointsOK(:,1:S:end);
%tic;normals = points2normals_radius(pointsOK);toc;
normals = points2normals(pointsOK);

%{
figure,
s =1;
quiver3(pointsOK(1,1:S*s:end),pointsOK(2,1:S*s:end),pointsOK(3,1:S*s:end),normals(1,1:S*s:end),normals(2,1:S*s:end),normals(3,1:S*s:end));

quiver3(pointsOK(1,1:S*s:end),pointsOK(2,1:S*s:end),pointsOK(3,1:S*s:end),normals(1,1:s:end),normals(2,1:s:end),normals(3,1:s:end));
figure,
indxxx = B == b;
pointsOKxxx = pointsOK(:,1:S:end);
quiver3(pointsOKxxx(1,indxxx),pointsOKxxx(2,indxxx),pointsOKxxx(3,indxxx),normals(1,indxxx),normals(2,indxxx),normals(3,indxxx));
hold on
quiver3(pointsOK(1,1:S*s:end),pointsOK(2,1:S*s:end),pointsOK(3,1:S*s:end),nrm(1,1:s:end),nrm(2,1:s:end),nrm(3,1:s:end),'-.r');

figure,
plot3(sphere(1,:),sphere(2,:),sphere(3,:),'.')

%}

% approximately 1313 bins
sphere = icosahedron2sphere(12)';
bins = sphere(:, sphere(1, :) >= 0);

%NSAMPLE = 1e5;
%sampleind = randsample(1 : size(normals, 2), min(size(normals, 2), NSAMPLE));
%normals = normals(:,sampleind );
[D, B] = max(abs(bins' * normals), [], 1);

H = accumarray(cat(2, B', repmat(1, [length(B) 1])), repmat(1, [length(B) 1]));

A = eye(3);
[~, I] = sort(-H);
for j = 1 : 3
    if ~isempty(I)
          b = I(1);
          % choose mean normal that falls into the biggest bin
          in_bin = normals(:, B == b);
          % flip mirrored normals
          dots = sum(in_bin .* repmat(bins(:, b), [1 size(in_bin, 2)]), 1);
          in_bin(:, (dots < 0)) = -in_bin(:, (dots < 0));
          v = mean(in_bin, 2);
          v = v / norm(v);

          A(:, j) = v;

          fprintf('Bin: %d, Normal: %f %f %f. Contains %d points. Mean vector: %f %f %f\n', b, bins(:, b), H(b), v);
          % remove bins that are not ~90 degrees away
          dots = sum(bins(:, I) .* repmat(v, [1 length(I)]), 1);
          I = I((dots >= cos(deg2rad(110))) & (dots <= cos(deg2rad(70))));
    end
end

axisI = A(:,1);
axisII = A(:,2);
axisII = axisII - (axisI'*axisII)*axisI;
axisII =axisII/norm(axisII);
axisIII = cross(axisI,axisII);

AA =[axisI,axisII,axisIII -1*[axisI,axisII,axisIII]];
[~, zi] = max(squeeze(cameraRt(1:3, 3, :))'*AA);
ZZ = AA(:, zi);
[~, xi] = max(squeeze(cameraRt(1:3, 1, :))'*AA);
XX = AA(:, xi);
[~, yi] = max(squeeze(cameraRt(1:3, 2, :))'*AA);
YY = AA(:, yi);
%{
for i =1:3,
    hold on;
    quiver3(1,1,1,AA(1,i),AA(2,i),AA(3,i));
quiver3(0,0,0,A(1,i),A(2,i),A(3,i));
pause;
end
axis tight;
%}
R = [XX YY ZZ]';
q = quaternion.rotateutov(ZZ, [0;0;1]);
Rtilt = RotationMatrix(q);
world_center = nanmean(reshape(pts,3,[]),2);

function rad = deg2rad(deg)

rad = deg*pi/180;

return;

function [coor,tri] = icosahedron2sphere(level)

% copyright by Jianxiong Xiao http://mit.edu/jxiao
% this function use a icosahedron to sample uniformly on a sphere
%{
Please cite this paper if you use this code in your publication:
J. Xiao, T. Fang, P. Zhao, M. Lhuillier, and L. Quan
Image-based Street-side City Modeling
ACM Transaction on Graphics (TOG), Volume 28, Number 5
Proceedings of ACM SIGGRAPH Asia 2009
%}

a= 2/(1+sqrt(5));
M=[
    0 a -1 a 1 0 -a 1 0
    0 a 1 -a 1 0 a 1 0
    0 a 1 0 -a 1 -1 0 a
    0 a 1 1 0 a 0 -a 1
    0 a -1 0 -a -1 1 0 -a
    0 a -1 -1 0 -a 0 -a -1
    0 -a 1 a -1 0 -a -1 0
    0 -a -1 -a -1 0 a -1 0
    -a 1 0 -1 0 a -1 0 -a
    -a -1 0 -1 0 -a -1 0 a
    a 1 0 1 0 -a 1 0 a
    a -1 0 1 0 a 1 0 -a
    0 a 1 -1 0 a -a 1 0
    0 a 1 a 1 0 1 0 a
    0 a -1 -a 1 0 -1 0 -a
    0 a -1 1 0 -a a 1 0
    0 -a -1 -1 0 -a -a -1 0
    0 -a -1 a -1 0 1 0 -a
    0 -a 1 -a -1 0 -1 0 a
    0 -a 1 1 0 a a -1 0
    ];

coor = reshape(M',3,60)';
%[M(:,[1 2 3]); M(:,[4 5 6]); M(:,[7 8 9])];

[coor, ~, idx] = unique(coor,'rows');

tri = reshape(idx,3,20)';

%{
for i=1:size(tri,1)
    x(1)=coor(tri(i,1),1);
    x(2)=coor(tri(i,2),1);
    x(3)=coor(tri(i,3),1);
    y(1)=coor(tri(i,1),2);
    y(2)=coor(tri(i,2),2);
    y(3)=coor(tri(i,3),2);
    z(1)=coor(tri(i,1),3);
    z(2)=coor(tri(i,2),3);
    z(3)=coor(tri(i,3),3);
    patch(x,y,z,'r');
end

axis equal
axis tight
%}

% extrude
coor = coor ./ repmat(sqrt(sum(coor .* coor,2)),1, 3);

for i=1:level
    m = 0;
    for t=1:size(tri,1)
        n = size(coor,1);
        coor(n+1,:) = ( coor(tri(t,1),:) + coor(tri(t,2),:) ) / 2;
        coor(n+2,:) = ( coor(tri(t,2),:) + coor(tri(t,3),:) ) / 2;
        coor(n+3,:) = ( coor(tri(t,3),:) + coor(tri(t,1),:) ) / 2;

        triN(m+1,:) = [n+1     tri(t,1)    n+3];
        triN(m+2,:) = [n+1     tri(t,2)    n+2];
        triN(m+3,:) = [n+2     tri(t,3)    n+3];
        triN(m+4,:) = [n+1     n+2         n+3];

        n = n+3;

        m = m+4;

    end
    tri = triN;

    % uniquefy
    [coor, ~, idx] = unique(coor,'rows');
    tri = idx(tri);

    % extrude
    coor = coor ./ repmat(sqrt(sum(coor .* coor,2)),1, 3);
end

% vertex number: 12  42  162  642

author: jianxiongxiao

repo: ProfXkit


text:
%% rastrigrin function
[x,y]=meshgrid(-5:0.1:5,-5:0.1:5);
     z=x.^2+y.^2-10*cos(2*pi*x)-10*cos(2*pi*y)+20;
     mesh(x,y,z)


text:
function g = gaussianPriorGradient(prior, x)

% GAUSSIANPRIORGRADIENT Gradient wrt x of the log Gaussian prior.

% PRIOR

% Compute gradient of prior
g = -prior.precision*x;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% The COBRAToolbox: testReactionEq.m
%
% Purpose:
%     - testReactionEq tests the functionality of ReactionEq in the
%     rBioNet extension
%
% Authors:
%     - Stefania Magnusdottir April 2017
%

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testReactionEq'));
cd(fileDir);

% load the reaction database
load([fileDir, filesep 'rxn.mat'])

% check match for made up reaction
matched = ReactionEq('etoh[e] + 2pg[c] -> atp[c] + h[e]', rxn);

% test that no reaction was matched
assert(isempty(matched))

% check match for existing reaction
matched = ReactionEq(rxn{1, 3}, rxn);

% test that reaction was matched
assert(~isempty(matched))

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function gui_check_scalar(edithandle,resetstr)
% check whether the the given edit field contains a scalar; warn and reset if not
if ~exist('resetstr','var')
    resetstr = ''; end

str = get(edithandle,'String'); num = str2num(str);
if ~isempty(str) && numel(num) ~= 1
    errordlg2('This field must be a number.');
    set(edithandle,'String',resetstr);
end

repo: REST


text:
function inter = learn_struct_dbn_reveal(seqs, ns, max_fan_in, penalty)
% LEARN_STRUCT_DBN_REVEAL Learn inter-slice adjacency matrix given fully observable discrete time series
% inter = learn_struct_dbn_reveal(seqs, node_sizes, max_fan_in, penalty)
%
% seqs{l}{i,t} = value of node i in slice t of time-series l.
%   If you have a single time series in an N*T array D, use
%      seqs = { num2cell(D) }.
%   If you have L time series, each of length T, in an N*T*L array D, use
%      seqs= cell(1,L); for l=1:L, seqs{l} = num2cell(D(:,:,l)); end
%   or, in vectorized form,
%      seqs = squeeze(num2cell(num2cell(D),[1 2]));
% Currently the data is assumed to be discrete (1,2,...)
%
% node_sizes(i) is the number of possible values for node i
% max_fan_in is the largest number of parents we allow per node (default: N)
% penalty is weight given to the complexity penalty (default: 0.5)
%  A penalty of 0.5 gives the BIC score.
%  A penalty of 0 gives the ML score.
%  Maximizing likelihood is equivalent to maximizing mutual information between parents and child.
%
% inter(i,j) = 1 iff node in slice t connects to node j in slice t+1
%
% The parent set for each node in slice 2 is computed by evaluating all subsets of nodes in slice 1,
% and picking the largest scoring one. This takes O(n^k) time per node, where n is the num. nodes
% per slice, and k <= n is the max fan in.
% Since all the nodes are observed, we do not need to use an inference engine.
% And since we are only learning the inter-slice matrix, we do not need to check for cycles.
%
% This algorithm is described in
% - "REVEAL: A general reverse engineering algorithm for inference of genetic network
%      architectures", Liang et al. PSB 1998
% - "Extended dependency analysis of large systems",
%       Roger Conant, Intl. J. General Systems, 1988, vol 14, pp 97-141
% - "Learning the structure of DBNs", Friedman, Murphy and Russell, UAI 1998.

n = length(ns);

if nargin < 3, max_fan_in = n; end
if nargin < 4, penalty = 0.5; end

inter = zeros(n,n);

if ~iscell(seqs)
  data{1} = seqs;
end

nseq = length(seqs);
nslices = 0;
data = cell(1, nseq);
for l=1:nseq
  nslices = nslices + size(seqs{l}, 2);
  data{l} = cell2num(seqs{l})'; % each row is a case
end
ndata = nslices - nseq; % subtract off the initial slice of each sequence

% We concatenate the sequences as in the following example.
% Let there be 2 sequences of lengths 4 and 5, with n nodes per slice,
% and let i be the target node.
% Then we construct following matrix D
%
% s{1}{1,1} ... s{1}{1,3}     s{2}{1,1} ... s{2}{1,4}
% ....
% s{1}{n,1} ... s{1}{n,3}     s{2}{n,1} ... s{2}{n,4}
% s{1}{i,2} ... s{1}{i,4}     s{2}{i,2} ... s{2}{i,5}
%
% D(1:n, i) is the i'th input and D(n+1, i) is the i'th output.
%
% We concatenate each sequence separately to avoid treating the transition
% from the end of one sequence to the beginning of another as a "normal" transition.

for i=1:n
  D = [];
  for l=1:nseq
    T = size(seqs{l}, 2);
    A = cell2num(seqs{l}(:, 1:T-1));
    B = cell2num(seqs{l}(i, 2:T));
    C = [A;B];
    D = [D C];
  end
  SS = subsets(1:n, max_fan_in, 1); % skip the empty set
  nSS = length(SS);
  bic_score = zeros(1, nSS);
  ll_score = zeros(1, nSS);
  target = n+1;
  ns2 = [ns ns(i)];
  for h=1:nSS
    ps = SS{h};
    dom = [ps target];
    counts = compute_counts(D(dom, :), ns2(dom));
    CPT = mk_stochastic(counts);
    [bic_score(h), ll_score(h)] = bic_score_family(counts, CPT, ndata);
  end
  if penalty == 0
    h = argmax(ll_score);
  else
    h = argmax(bic_score);
  end
  ps = SS{h};
  inter(ps, i) = 1;
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt