text:
function [eg] = dag_to_essential_graph(dag)
cpdag = dag_to_cpdag(dag);
eg = dag + dag .* (cpdag + cpdag');

return;

% Coverts a DAG into Essential Graph where edges are coded by 2 and 3, 2 is
% directed edge and 3 is bidirected edge and is at one (the same as the original DAG) of the two
% symetrical places.

% Is implemented by the algorithm of Max Chickering in D.M.Chickering (1995).
% A transformational characterization of equivalent Bayesian network structures.
% In Proceedings of Eleventh Conference on Uncertainty in Artificial Intelligence, Montreal, QU,
% pages 87-98. Morgan Kaufmann
% http://research.microsoft.com/~dmax/publications/uai95.pdf

% Implemented by Tomas Kocka, AAU.

function [eg] = dag_to_essential_graph(dagx)

%print_dag(dagx); % Just checking input

order = topological_sort(dagx); % get the topological order of nodes and their number

% fprintf('the topological order is: %d',order);
% fprintf('\n');

[nx,ny] = size(dagx); % gets the number of nodes, note that nx == ny
[I,J] = find(dagx); % finds all nonzero elements in the adjacency matrix, i.e. arcs in the DAG - however we will overwrite it in a special order
% we will sort the arcs from lowest possible y and highest possible x, arcs are x->y
e = 1;
for y = 1:ny
    for x = nx:-1:1
        %fprintf('x %d ',order(x)); fprintf('y %d ',order(y));
        if dagx(order(x),order(y)) == 1
            I(e) = order(x);
            J(e) = order(y);
            e = e + 1;
            %fprintf('x order %d',x);
            %fprintf('y order %d',y);
            %fprintf('\n');
        end
    end
end

% fprintf('the arcs are: %d',I);
% fprintf('\n');
% fprintf('the arcs are: %d',J);
% fprintf('\n');

% Now we have to decide which arcs are part of the essential graph and
% which are undirected edges in the essential graph.
% Undecided arc in the DAG are 1, directed in EG are 2 and undirected in EG
% are 3.

for e = 1:length(I)
    if dagx(I(e),J(e)) == 1
        cont = true;
        for w = 1:nx
            if dagx(w,I(e)) == 2
                if dagx(w,J(e)) ~= 0
                    dagx(w,J(e)) = 2;
                else
                    for ww = 1:nx
                        if dagx(ww,J(e)) ~= 0
                           dagx(ww,J(e)) = 2;
                        end
                    end % and now skip the rest and start with another arc from the list
                    w = nx;
                    cont = false;
                end
            end
        end
        if cont
           exists = false;
           for z = 1:nx
               %fprintf('test %d',dagx(z,J(e)));
               if dagx(z,J(e)) ~= 0 & z ~= I(e) & dagx(z,I(e)) == 0
                  exists = true;
                  for ww = 1:nx
                        if dagx(ww,J(e)) == 1
                           dagx(ww,J(e)) = 2;
                        end
                  end
               end
           end
           if ~ exists
               for ww = 1:nx
                   if dagx(ww,J(e)) == 1
                      dagx(ww,J(e)) = 3;
                   end
               end
           end
        end
    end
end

%print_dag(dagx); % Just checking output

author: bayesnet

repo: bnt


text:
function geometry_test007 ( )

%*****************************************************************************80
%
%% TEST007 tests BALL_UNIT_SAMPLE_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  n_sample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST007\n' );
  fprintf ( 1, '  For the unit ball in 3 dimensions:\n' );
  fprintf ( 1, '  BALL_UNIT_SAMPLE_3D samples;\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  A few sample values:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 5
    [ x, seed ] = ball_unit_sample_3d ( seed );
    fprintf ( 1, '  %10f  %10f  %10f\n', x(1:dim_num) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of sample points = %d\n', n_sample );

  average(1:dim_num) = 0.0;

  for i = 1 : n_sample
    [ x, seed ] = ball_unit_sample_3d ( seed );
    average(1:dim_num) = average(1:dim_num) + x(1:dim_num);
  end

  average(1:dim_num) = average(1:dim_num) / n_sample;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now average the points, which should get a value\n' );
  fprintf ( 1, '  close to zero, and closer as N_SAMPLE increases.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Average:        %10f  %10f  %10f\n', average(1:dim_num) );

  average_r = 0.0;

  for i = 1 : n_sample
    [ x, seed ] = ball_unit_sample_3d ( seed );
    r = sqrt ( sum ( x(1:dim_num).^2 ) );
    average_r = average_r + r;
  end

  average_r = average_r / n_sample;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now average the distance of the points from\n' );
  fprintf ( 1, '  the center, which should be the \n' );
  fprintf ( 1, '  1/2**(1/dim_num) = %f\n', 0.5^( 1.0 / dim_num ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Average:        %f\n', average_r );

  average_theta = 0.0;

  for i = 1 : n_sample
    [ x, seed ] = ball_unit_sample_3d ( seed );
    theta = r8_atan ( x(2), x(1) );
    average_theta = average_theta + theta;
  end

  average_theta = average_theta / n_sample;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now average the angle THETA,\n' );
  fprintf ( 1, '  which should be PI.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Average:        %f\n', average_theta );

  average_phi = 0.0;

  for i = 1 : n_sample
    [ x, seed ] = ball_unit_sample_3d ( seed );
    r = sqrt ( sum ( x(1:dim_num).^2 ) );
    phi = acos ( x(3) / r );
    average_phi = average_phi + phi;
  end

  average_phi = average_phi / n_sample;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now average the angle PHI,\n' );
  fprintf ( 1, '  which should be PI/2.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Average:        %f\n', average_phi );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef IQR < Algorithm

    methods (Access = public)

        function obj = IQR()
            obj.name = 'IQR';
            obj.inputPort = DataType.kSignal;
            obj.outputPort = DataType.kFeature;
        end

        function result = compute(~,signal)
            result = iqr(signal);
        end

        function metrics = computeMetrics(~,input)
            n = size(input,1);
            flops = 57 * n;
            memory = n;
            outputSize = Constants.kFeatureBytes;
            metrics = Metric(flops,memory,outputSize);
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function [ ap, rcond, z, info ] = dppco ( ap, n )

%*****************************************************************************80
%
%% DPPCO factors a real symmetric positive definite matrix in packed form.
%
%  Discussion:
%
%    DPPCO also estimates the condition of the matrix.
%
%    If RCOND is not needed, DPPFA is slightly faster.
%
%    To solve A*X = B, follow DPPCO by DPPSL.
%
%    To compute inverse(A)*C, follow DPPCO by DPPSL.
%
%    To compute determinant(A), follow DPPCO by DPPDI.
%
%    To compute inverse(A), follow DPPCO by DPPDI.
%
%  Packed storage:
%
%    The following program segment will pack the upper triangle of
%    a symmetric matrix.
%
%      k = 0
%      do j = 1, n
%        do i = 1, j
%          k = k + 1
%          ap(k) = a(i,j)
%        end
%      end
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real AP(N*(N+1)/2), the packed form of a symmetric matrix A.
%    The columns of the upper triangle are stored sequentially in a
%    one-dimensional array.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real AP(N*(N+1)/2), an upper riangular matrix R, stored
%    in packed form, so that A = R'*R.  If INFO /= 0, the factorization
%    is not complete.
%
%    Output, real RCOND, an estimate of the reciprocal condition
%    of A.  For the system A*X = B, relative perturbations in A and B of size
%    EPSILON may cause relative perturbations in X of size EPSILON/RCOND.
%    If RCOND is so small that the logical expression
%      1.0 + RCOND == 1.0D+00
%    is true, then A may be singular to working precision.  In particular,
%    RCOND is zero if exact singularity is detected or the estimate underflows.
%
%    Output, real Z(N), a work vector whose contents are usually
%    unimportant.  If A is singular to working precision, then Z is an
%    approximate null vector in the sense that
%      norm(A*Z) = RCOND * norm(A) * norm(Z).
%    If INFO /= 0, Z is unchanged.
%
%    Output, integer INFO, error flag.
%    0, for normal return.
%    K, signals an error condition.  The leading minor of order K is
%    not positive definite.
%

%
%  Find the norm of A.
%
  j1 = 1;
  for j = 1 : n
    z(j) = dasum ( j, ap(j1:j1+j-1), 1 );
    ij = j1;
    j1 = j1 + j;
    for i = 1 : j-1
      z(i) = z(i) + abs ( ap(ij) );
      ij = ij + 1;
    end
  end

  anorm = max ( z(1:n) );
%
%  Factor.
%
  [ ap, info ] = dppfa ( ap, n );

  if ( info ~= 0 )
    return
  end
%
%  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
%
%  Estimate = norm(Z)/norm(Y) where A * Z = Y and A * Y = E.
%
%  The components of E are chosen to cause maximum local
%  growth in the elements of W where R'*W = E.
%
%  The vectors are frequently rescaled to avoid overflow.
%
%  Solve R' * W = E.
%
  ek = 1.0;
  z(1:n) = 0.0;

  kk = 0;

  for k = 1 : n

    kk = kk + k;

    if ( z(k) ~= 0.0 )
      ek = -abs ( ek ) * r8_sign ( z(k) );
    end

    if ( ap(kk) < abs ( ek - z(k) ) )
      s = ap(kk) / abs ( ek - z(k) );
      z(1:n) = s * z(1:n);
      ek = s * ek;
    end

    wk = ek - z(k);
    wkm = -ek - z(k);
    s = abs ( wk );
    sm = abs ( wkm );
    wk = wk / ap(kk);
    wkm = wkm / ap(kk);
    kj = kk + k;

    if ( k + 1 <= n )

      for j = k + 1 : n
        sm = sm + abs ( z(j) + wkm * ap(kj) );
        z(j) = z(j) + wk * ap(kj);
        s = s + abs ( z(j) );
        kj = kj + j;
      end

      if ( s < sm )

        t = wkm - wk;
        wk = wkm;
        kj = kk + k;

        for j = k+1 : n
          z(j) = z(j) + t * ap(kj);
          kj = kj + j;
        end

      end

    end

    z(k) = wk;

  end

  z(1:n) = z(1:n) / dasum ( n, z(1:n), 1 );
%
%  Solve R * Y = W.
%
  for k = n : -1 : 1

    if ( ap(kk) < abs ( z(k) ) );
      s = ap(kk) / abs ( z(k) );
      z(1:n) = s * z(1:n);
    end

    z(k) = z(k) / ap(kk);
    kk = kk - k;
    t = -z(k);
    z(1:k-1) = daxpy ( k-1, t, ap(kk+1:kk+k-1), 1, z(1:k-1), 1 );

  end

  z(1:n) = z(1:n) / dasum ( n, z(1:n), 1 );

  ynorm = 1.0;
%
%  Solve R' * V = Y.
%
  for k = 1 : n

    z(k) = z(k) - ddot ( k-1, ap(kk+1:kk+k-1), 1, z(1:k-1), 1 );
    kk = kk + k;

    if ( ap(kk) < abs ( z(k) ) )
      s = ap(kk) / abs ( z(k) );
      z(1:n) = s * z(1:n);
      ynorm = s * ynorm;
    end

    z(k) = z(k) / ap(kk);

  end

  s = 1.0 / dasum ( n, z(1:n), 1 );
  z(1:n) = s * z(1:n);
  ynorm = s * ynorm;
%
%  Solve R * Z = V.
%
  for k = n : -1 : 1

    if ( ap(kk) < abs ( z(k) ) )
      s = ap(kk) / abs ( z(k) );
      z(1:n) = s * z(1:n);
      ynorm = s * ynorm;
    end

    z(k) = z(k) / ap(kk);
    kk = kk - k;
    t = -z(k);
    z(1:k-1) = daxpy ( k-1, t, ap(kk+1:kk+k-1), 1, z(1:k-1), 1 );

  end
%
%  Make ZNORM = 1.0.
%
  s = 1.0 / dasum ( n, z(1:n), 1 );
  z(1:n) = s * z(1:n);
  ynorm = s * ynorm;

  if ( anorm ~= 0.0 )
    rcond = ynorm / anorm;
  else
    rcond = 0.0;
  end

  return
end


text:
% This script generates randomly several testing
% Quadratic Programming tasks.
%
% The generated tasks are stored to ../data/ folder.
%

randn('state',0);
rand('state',0);

%=======================================================
% Quadratic Programming task with simplex constraints
% solved by LIBQP_SPLX
%=======================================================

% TASK 1: splxcon_qp_random100x10
% number of variables 100
% number of equlity and inequality constraints 10 (not considering constraints x >= 0)
I = [];
for i=1:10;
  I = [I i*ones(1,10)];
end
Z = randn(length(I),length(I)*2); H=Z*Z';
f=randn(length(I),1);
b = rand(10,1)*10;

% 5 randomly selected constraints will be equlities and inequlities
S = zeros(10,1);
idx = randperm(10);
S(idx(1:5)) = 1;

save ../data/splx_qp_random100x10 H f b I S;

% TASK 2: splxcon_qp_random300x30
% number of variables 300
% number of equlity and inequality constraints 30 (not considering constraints x >= 0)
I = [];
for i=1:30;
  I = [I i*ones(1,10)];
end
Z = randn(length(I),length(I)*2); H=Z*Z';
f=randn(length(I),1);
b = rand(30,1)*10;

% 5 randomly selected constraints will be equlity and inequlity
S = zeros(30,1);
idx = randperm(30);
S(idx(1:15)) = 1;

save ../data/splx_qp_random300x30 H f b I S;

%=======================================================
% Quadratic Programming task solved by GSMO algorithm
% solved by LIBQP_GSMO
%=======================================================

% TASK 1: gsmo_qp_random100
% number of variables 100

f=randn(100,1)*20;
Z = randn(length(f),length(f)*2); H=Z*Z';

a = randn(1,length(f)); a(find(a==0)) = 1;
b = randn(1);
UB = rand(length(f),1);
LB = zeros(length(f),1);

save ../data/gsmo_qp_random100 H f a b UB LB;

% TASK 2: gsmo_qp_random300
% number of variables 300

f=randn(300,1)*20;
Z = randn(length(f),length(f)*2); H=Z*Z';

a = randn(1,length(f)); a(find(a==0)) = 1;
b = randn(1);
UB = zeros(length(f),1);
LB = -rand(length(f),1);

save ../data/gsmo_qp_random300 H f a b UB LB;

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function val = lval(fun)
%LVAL  Obtain the value of a FUN at its left endpoint.
%   LVAL(F) returns the value of F at F.domain(1).

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

% Call LVAL() of the ONEFUN:
val = lval(fun.onefun);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
sim('pend2');

a = logsout.angles;
f = logsout.force;

figure(5)
plot(a.time, a.data*180/pi)
xlabel('time')
ylabel('angle [degree]')
grid on;

figure(6)
plot(f.time, f.data)
xlabel('time')
ylabel('force [N]')
grid on;

author: Sable

repo: mcbench-benchmarks


text:
function L = eq(p1,p2)
%PREAL/EQ Overloaded EQ (==) operator for class preal.

global useUnitsFlag

if ~(useUnitsFlag) % If physunits is disabled...
    L=eq(double(p1),double(p2)); % ... treat as double.
    return
end

p1=preal(p1);
p2=preal(p2);
tol=0.001;
if isscalar(p1)
    p1=preal(p1.value*ones(size(p2)),p1.units);
elseif isscalar(p2)
    p2=preal(p2.value*ones(size(p1)),p2.units);
elseif ndims(p1)==ndims(p2)&size(p1)==size(p2)
    % do nothing
else
    error('Array dimensions must agree')
end

L=false(size(p1));
for k=1:numel(L)
    if any(abs(p1(k).units-p2(k).units)>tol)
        error('Attempt to compare dimensionally inconsistent preals')
    end
    L(k)=eq(p1(k).value,p2(k).value);
end


text:
function result=plotbox(Target,Template,M);
%
%By Alaa Eleyan May,2008
%*********************************************************
[r1,c1]=size(Target);
[r2,c2]=size(Template);

[r,c]=max(M);
[r3,c3]=max(max(M));

i=c(c3);
j=c3;
result=Target;
for x=i:i+r2-1
   for y=j
       result(x,y)=255;
   end
end
for x=i:i+r2-1
   for y=j+c2-1
       result(x,y)=255;
   end
end
for x=i
   for y=j:j+c2-1
       result(x,y)=255;
   end
end
for x=i+r2-1
   for y=j:j+c2-1
       result(x,y)=255;
   end
end


text:
clear all; close all; clc

t=0:0.01:100;
x0=[5 5 5];

[t,Y]=ode45('lrz_rhs',t,x0);

x=Y(:,1); y=Y(:,2); z=Y(:,3);
plot3(x,y,z,'Linewidth',[2]), grid on

figure(2)
plot(t,x,t,y,t,z,'Linewidth',[2])

%%
figure(2)
plot(t,x,'Linewidth',[2])

%%
H=[x(1:5000).'
   x(2:5001).'
   x(3:5002).'
   x(4:5003).'
   x(5:5004).'
   x(6:5005).'
   x(7:5006).'
   x(8:5007).'];

[u,s,v]=svd(H,'econ');
figure(3), plot(diag(s)/(sum(diag(s))),'ro','Linewidth',[3])

%%
figure(4), plot(u(:,1:5),'Linewidth',[2])
figure(5), plot3(v(:,1),v(:,2),v(:,3),'Linewidth',[2])

figure(6), plot(v(:,4),'Linewidth',[2])

author: dynamicslab


text:
function value = r8_tanh ( x )

%*****************************************************************************80
%
%% R8_TANH evaluates the hyperbolic tangent of an R8 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the hyperbolic tangent of X.
%
  persistent nterms
  persistent sqeps
  persistent tanhcs
  persistent xmax

  if ( isempty ( nterms ) )
    tanhcs = [ ...
      -0.25828756643634710438338151450605, ...
      -0.11836106330053496535383671940204, ...
      +0.98694426480063988762827307999681E-02, ...
      -0.83579866234458257836163690398638E-03, ...
      +0.70904321198943582626778034363413E-04, ...
      -0.60164243181207040390743479001010E-05, ...
      +0.51052419080064402965136297723411E-06, ...
      -0.43320729077584087216545467387192E-07, ...
      +0.36759990553445306144930076233714E-08, ...
      -0.31192849612492011117215651480953E-09, ...
      +0.26468828199718962579377758445381E-10, ...
      -0.22460239307504140621870997006196E-11, ...
      +0.19058733768288196054319468396139E-12, ...
      -0.16172371446432292391330769279701E-13, ...
      +0.13723136142294289632897761289386E-14, ...
      -0.11644826870554194634439647293781E-15, ...
      +0.98812684971669738285540514338133E-17, ...
      -0.83847933677744865122269229055999E-18, ...
      +0.71149528869124351310723506176000E-19, ...
      -0.60374242229442045413288837119999E-20, ...
      +0.51230825877768084883404663466666E-21, ...
      -0.43472140157782110106047829333333E-22, ...
      +0.36888473639031328479423146666666E-23, ...
      -0.31301874774939399883325439999999E-24, ...
      +0.26561342006551994468488533333333E-25, ...
      -0.22538742304145029883494399999999E-26, ...
      +0.19125347827973995102208000000000E-27, ...
      -0.16228897096543663117653333333333E-28, ...
      +0.13771101229854738786986666666666E-29, ...
      -0.11685527840188950118399999999999E-30, ...
      +0.99158055384640389120000000000000E-32 ]';
    nterms = r8_inits ( tanhcs, 31, 0.1 * r8_mach ( 3 ) );
    sqeps = sqrt ( 3.0 * r8_mach ( 3 ) );
    xmax = - 0.5 * log ( r8_mach ( 3 ) );
  end

  y = abs ( x );

  if ( y <= sqeps )

    value = x;

  elseif ( y <= 1.0 )

    value = x * ( 1.0 + r8_csevl ( 2.0 * x * x - 1.0, tanhcs, nterms ) );

  elseif ( y <= xmax )

    y = exp ( y );
    yrec = 1.0 / y;
    value = ( y - yrec ) / ( y + yrec );

    if ( x < 0.0 )
      value = - value;
    end

  else

    if ( x < 0.0 )
      value = - 1.0;
    else
      value = + 1.0;
    end

  end

  return
end


text:
function R = subsref(logger, ss)
%SUBSREF Get the properties by subscription
%
% $ Syntax $
%   - R = logger.<propname>
%   - R = logger.('<propname>')
%
% $ History $
%   - Created by Dahua Lin on Aug 12nd, 2006
%

if ss(1).type == '.' && ischar(ss(1).subs)
    R = get(logger, ss(1).subs);
    if (length(ss) >= 2)
        R = subsref(R, ss(2:end));
    end
else
    error('sltoolbox:invalidarg', 'Invalid subscription for dataset');
end

author: lmthang


text:
%% Example:
% Compute the tridiagonal matrix used to solve implicitly a laplace (heat)
% equation.

m = 5; % domain size.
h = 0.1; % delta x value.

I = eye(m);
e = ones(1,m);
x = spdiags([e' -4*e' e'],[-1 0 1],m,m);
%x = zeros(m);
%x = x+spdiags([e' -4*e' e'],[-1 0 1],m,m);
s = spdiags([e',e'],[-1 1],m,m);
s = zeros(m);
s = s+spdiags([e',e'],[-1 1],m,m);
%kron(I,x);
%kron(s,I);
A = (kron(I,x)+kron(s,I))/h^2;
%size(A);

repo: Aero-matlab


text:
function child  = crossoverParents(p1,p2)
% Oren Rosen
% The MathWorks
% 8/29/2007
%
% This custom crossover function is written to work on a population of
% vectors of zeros and ones with the same amount of ones in each vector.
% The children that are produced from 2 parents will have the same genes
% for every element they agree on, and random choices of zerps and ones for
% the elements they don't agree on. The number of each is set so that all
% children have the same number of ones as their parents. All children
% automatically satisfy this constraint so there is no need to impose these
% constraints.

% *** Calculate dimensions of p1 (assume p2 matches) ***
numBits = length(p1);
num1s = sum(p1);

% *** Initialize child to all zeros ***
child = zeros(1,numBits);

% *** Need this later ***
indexVec = 1:numBits;

% *** Find Matching 1's and 0's ***
% Ex: If           p1 == [ 1 0 1 0 0 1 1 0 0 0 ]
%                  p2 == [ 1 0 0 1 0 0 1 0 1 0 ]
%     Then matching1s == [ 1 0 0 0 0 0 1 0 0 0 ]
%     Then matching0s == [ 0 1 0 0 1 0 0 1 0 1 ]
matching1s = ~xor(p1,p2) & (p1 == 1);
matching0s = ~xor(p1,p2) & (p1 == 0);

% *** Find Matching Indices ***
%     If       matching1s == [ 1 0 0 0 0 0 1 0 0 0 ]
%              matching0s == [ 0 1 0 0 1 0 0 1 0 1 ]
%     Then matching1sIndx == [ 1 7 ]
%          matching0sIndx == [ 2 5 8 10 ]
%         nonmatchingIndx == [ 3 4 6 9 ]
matching1sIndx = indexVec(matching1s);
matching0sIndx = indexVec(matching0s);
nonmatchingIndx = setdiff(indexVec,[matching0sIndx,matching1sIndx]);

% *** Create Child ***
% Ex: If   num1s == 4
%          matching1sIndx == [ 1 7 ]
%          nonmatchingIndx == [ 3 4 6 9 ]
%     Then numMatching1s == 2
%     num1sToFill == 2
%     Indx1sToFill == 2 random choices from [ 3 4 6 9 ]
numMatching1s = numel(matching1sIndx);
num1sToFill = num1s - numMatching1s;
Indx1sToFill = randsample(nonmatchingIndx,num1sToFill);

% *** Fill in 1s ***
% Ex: If      p1 == [ 1 0 1 0 0 1 1 0 0 0 ]
%             p2 == [ 1 0 0 1 0 0 1 0 1 0 ]
%     Then child == [ 1 0 ? ? 0 ? 1 0 ? 0 ]
%     With exactly 2 of the '?' equal to 1, the rest 0.
child(matching1sIndx) = 1;
child(Indx1sToFill) = 1;

% *** Display results ***
disp(' ');
disp(['Given: parent1 = [ ',num2str(p1),' ]']);
disp(['       parent2 = [ ',num2str(p2),' ]']);
disp(' ');
disp(['         child = [ ',num2str(child),' ]']);


text:
function areRPrime=areRelativelyPrime(a)
%%ARERELATIVELYPRIME Determine whether a set of numbers is relatively
%               prime (are coprime). A collections of numbers are
%               relatively prime if the greatest common divisor between all
%               pairs of the numbers is 1. Relatively prime numbers are
%               also known as coprime numbers.
%
%INPUTS: a An nXnumSets set of numSets sets of n-integers where one wishes
%          to determine whether the n numbers in each set are relatively
%          prime.
%
%OUTPUTS: areRPrime A numSetsX1 boolean array indicating whether the
%                   numbers in each set are coprime.
%
%Coprime numbers arise in a number of areas, such as in the Chinese
%remainder theorem. This function jst goes through all pairs of numbers and
%uses the gcd function to see whether the greatest common divisor between
%all pairs is 1.
%
%EXAMPLE:
% a=zeros(3,2);
% %This set of numbers is coprime even though 4 is not a prime number.
% a(:,1)=[7;9;4];
% %This set of numbers is not coprime
% a(:,2)=[3;2;9];
% areRPrime=areRelativelyPrime(a)
%The functions returns [1;0]
%
%September 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

N=size(a,1);
numSets=size(a,2);

areRPrime=true(numSets,1);

for curSet=1:numSets
    for n1=1:(N-1)
        for n2=(n1+1):N
            if(gcd(a(n1,curSet),a(n2,curSet))~=1)
                areRPrime(curSet)=false;
                break;
            end
        end

        if(areRPrime(curSet)~=true)
            break;
        end
    end
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 colors2 = cluster_orthviews_classes(cl,classes,overlay,myview,domontage, varargin)
% :Usage:
% ::
%
%    colors2 = cluster_orthviews_classes(cl,classes,overlay,myview,domontage, [orthview axis], [colors cell])
%
% Makes montage and cluster orthviews of classes of clusters in different
% colors (hard coded colors right now).
%
% :Inputs:
%
%   **cl:**
%        clusters structure with regions
%
%   **classes:**
%        integers with classes (i.e., networks) to be color-coded on plot
%
%   **overlay:**
%        anatomical underlay image
%
%   **myview:**
%        if entered, shows centers on plot
%
%   **domontage:**
%        if non-zero, make montages of networks
%
%   **[orthview axis]:**
%        Optional; integer for which orthviews axis to use, 1:kk
%
% :Output:
%
%   Cell of colors, in order, for use in other functions
%
% :Examples:
% ::
%
%    classes = c.ClusterSolution.classes;
%    overlay = EXPT.overlay;
%    cluster_orthviews_classes(cl,c.ClusterSolution.classes, EXPT.overlay, 'saggital', 0);
%    cluster_orthviews_classes(cl,c.ClusterSolution.classes, EXPT.overlay, 'axial', 0);
%    cluster_orthviews_classes(cl,c.ClusterSolution.classes, EXPT.overlay, 'coronal', 0);
%    cluster_orthviews_classes(cl,c.ClusterSolution.classes, EXPT.overlay, [], 0);
%
% ..
%    tor wager
%    Colors update, minor improvements, Jan 2010
% ..

whichorth = 1;
if ~isempty(varargin), whichorth = varargin{1}; end

if length(varargin) > 1
    colors2 = varargin{2};
    disp('Using input colors: '); %colors2{:}
else
    % text colors no longer needed -- can handle string colors
    %colors = {'yo' 'bo' 'go' 'ro' 'co' 'mo' 'ko' 'r^' 'g^' 'b^' 'y^' 'c^' 'm^' 'k^'};
    disp('Using default colors. ');
    colors2 = {[1 1 0] [0 0 1] [0 1 0] [1 0 0] [1 .5 0] [0 1 1] [1 0 1] [.5 1 0] [0 .5 1] [0 1 .5] [1 0 .5]};
end

while length(colors2) < max(classes)
    colors2tmp = colors2;
    for jj = 1:length(colors2tmp)
        colors2tmp{jj}(colors2tmp{jj} > .8) = colors2tmp{jj}(colors2tmp{jj} > .8) - .2;
        colors2tmp{jj}(colors2tmp{jj} < .2) = colors2tmp{jj}(colors2tmp{jj} < .2) + .2;
    end
    colors2 = [colors2 colors2tmp];
end
colors2 = colors2(1:max(classes));

% orthviews
cluster_orthviews(cl(classes==1),colors2(1),'overlay',overlay,'solid');

for i = 2:max(classes)
    cluster_orthviews(cl(classes==i),colors2(i),'add','solid', 'handle', whichorth);
end

if ~isempty(myview) && ischar(myview)
    cluster_orthviews_showcenters(cl,myview,overlay,0);
end

if domontage

% montage: build function call
str = ['montage_clusters(overlay,'];
for i = 1:max(classes)
    clc{i} = cl(classes==i);
    str = [str 'clc{' num2str(i) '},'];
end
str = [str '{'];
for i = 1:max(classes)
    %str = [str '''' colors2{i} ''''];  % no longer need solid colors
    str = [str '[' num2str(colors2{i}) '] '];
    %if i ~= max(classes), str = [str ' ']; end
end
str = [str '}, ''nooverlap'');'];
eval(str)
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [ a, seed ] = r83s_random ( n, seed )

%*****************************************************************************80
%
%% R83S_RANDOM randomizes an R83S matrix.
%
%  Discussion:
%
%    The R83S storage format is used for a tridiagonal scalar matrix.
%    The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
%    values that occur on every row.
%
%  Example:
%
%    Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
%    be interpreted:
%
%      A2  A3   0   0   0
%      A1  A2  A3   0   0
%       0  A1  A2  A3   0
%       0   0  A1  A2  A3
%       0   0   0  A1  A2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the linear system.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real A(3), the R83 matrix.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ a(1:3), seed ] = r8vec_uniform_01 ( 3, seed );

  return
end


text:
% Demonstration of Rank Aware Order Recursive Matching Pursuit.

close all;
clear all;
clc;
rng('default');
png_export = true;
pdf_export = false;
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');

mf = spx.graphics.Figures();

% Signal space
N = 256;
% Number of measurements
M = 64;
% Sparsity level
K = 8;
% Number of signals
S = 4;
% Construct the signal generator.
gen  = spx.data.synthetic.SparseSignalGenerator(N, K, S);
% Generate bi-uniform signals
X = gen.biUniform(1, 2);
% Sensing matrix
Phi = spx.dict.simple.gaussian_dict(M, N);
% Measurement vectors
Y = Phi.apply(X);
% Rank Aware ORMP MMV solver instance
solver = spx.pursuit.joint.RankAwareORMP(Phi, K);
% Solve the sparse recovery problem
result = solver.solve(Y);
% Solution vector
Z = result.Z;

for s=1:S
    mf.new_figure(sprintf('Rank Aware ORMP signal: %d', s));
    subplot(411);
    stem(X(:, s), '.');
    title('Sparse vector');
    subplot(412);
    stem(Z(:, s), '.');
    title('Recovered sparse vector');
    subplot(413);
    stem(abs(X(:, s) - Z(:, s)), '.');
    title('Recovery error');
    subplot(414);
    stem(Y(:, s), '.');
    title('Measurement vector');
end

author: indigits

repo: sparse-plex


text:
function h = plot_gaussian_ellipsoid(m, C, sdwidth, npts, axh,alphaChannel,color)
% PLOT_GAUSSIAN_ELLIPSOIDS plots 2-d and 3-d Gaussian distributions
%
% H = PLOT_GAUSSIAN_ELLIPSOIDS(M, C) plots the distribution specified by
%  mean M and covariance C. The distribution is plotted as an ellipse (in
%  2-d) or an ellipsoid (in 3-d).  By default, the distributions are
%  plotted in the current axes. H is the graphics handle to the plotted
%  ellipse or ellipsoid.
%
% PLOT_GAUSSIAN_ELLIPSOIDS(M, C, SD) uses SD as the standard deviation
%  along the major and minor axes (larger SD => larger ellipse). By
%  default, SD = 1. Note:
%  * For 2-d distributions, SD=1.0 and SD=2.0 cover ~ 39% and 86%
%     of the total probability mass, respectively.
%  * For 3-d distributions, SD=1.0 and SD=2.0 cover ~ 19% and 73%
%     of the total probability mass, respectively.
%
% PLOT_GAUSSIAN_ELLIPSOIDS(M, C, SD, NPTS) plots the ellipse or
%  ellipsoid with a resolution of NPTS (ellipsoids are generated
%  on an NPTS x NPTS mesh; see SPHERE for more details). By
%  default, NPTS = 50 for ellipses, and 20 for ellipsoids.
%
% PLOT_GAUSSIAN_ELLIPSOIDS(M, C, SD, NPTS, AX) adds the plot to the
%  axes specified by the axis handle AX.
%
% Examples:
% -------------------------------------------
%  % Plot three 2-d Gaussians
%  figure;
%  h1 = plot_gaussian_ellipsoid([1 1], [1 0.5; 0.5 1]);
%  h2 = plot_gaussian_ellipsoid([2 1.5], [1 -0.7; -0.7 1]);
%  h3 = plot_gaussian_ellipsoid([0 0], [1 0; 0 1]);
%  set(h2,'color','r');
%  set(h3,'color','g');
%
%  % "Contour map" of a 2-d Gaussian
%  figure;
%  for sd = [0.3:0.4:4],
%    h = plot_gaussian_ellipsoid([0 0], [1 0.8; 0.8 1], sd);
%  end
%
%  % Plot three 3-d Gaussians
%  figure;
%  h1 = plot_gaussian_ellipsoid([1 1  0], [1 0.5 0.2; 0.5 1 0.4; 0.2 0.4 1]);
%  h2 = plot_gaussian_ellipsoid([1.5 1 .5], [1 -0.7 0.6; -0.7 1 0; 0.6 0 1]);
%  h3 = plot_gaussian_ellipsoid([1 2 2], [0.5 0 0; 0 0.5 0; 0 0 0.5]);
%  set(h2,'facealpha',0.6);
%  view(129,36); set(gca,'proj','perspective'); grid on;
%  grid on; axis equal; axis tight;
% -------------------------------------------
%
%  Gautam Vallabha, Sep-23-2007, Gautam.Vallabha@mathworks.com

%  Revision 1.0, Sep-23-2007
%    - File created
%  Revision 1.1, 26-Sep-2007
%    - NARGOUT==0 check added.
%    - Help added on NPTS for ellipsoids

if ~exist('sdwidth', 'var'), sdwidth = 1; end
if ~exist('npts', 'var'), npts = []; end
if ~exist('axh', 'var'), axh = gca; end
if ~exist('alphaChannel','var'),alphaChannel=1;end
if ~exist('color','var'),color=[0 0 1];end

if numel(m) ~= length(m),
    error('M must be a vector');
end
if ~( all(numel(m) == size(C)) )
    error('Dimensionality of M and C must match');
end
if ~(isscalar(axh) && ishandle(axh) && strcmp(get(axh,'type'), 'axes'))
    error('Invalid axes handle');
end

set(axh, 'nextplot', 'add');

switch numel(m)
%   case 2, if alphaChannel,  h=show2d_noalpha(m(:),C,sdwidth,npts,axh,color); else h=show2d(m(:),C,sdwidth,npts,axh,alphaChannel,color); end
    case 2, h=show2d(m(:),C,sdwidth,npts,axh,alphaChannel,color);
    case 3, h=show3d(m(:),C,sdwidth,npts,axh);
   otherwise
      error('Unsupported dimensionality');
end

if nargout==0,
    clear h;
end

%-----------------------------
function h = show2d(means, C, sdwidth, npts, axh,alphaChannel,color)
if isempty(npts), npts=50; end
% plot the gaussian fits
tt=linspace(0,2*pi,npts)';
x = cos(tt); y=sin(tt);
ap = [x(:) y(:)]';
[v,d]=eig(C);
d = sdwidth * sqrt(d); % convert variance to sdwidth*sd
bp = (v*d*ap) + repmat(means, 1, size(ap,2));
%if alphaChannel == 1
%h = plot(bp(1,:), bp(2,:), '-', 'parent', axh);
%set(h,'color',color);
%else
h = patchline(bp(1,:),bp(2,:),'linewidth',2,'linestyle','-','edgealpha',alphaChannel,'edgecolor',color,'LineSmoothing','on');
%end

function h = show2d_noalpha(means, C, sdwidth, npts, axh,color)
if isempty(npts), npts=50; end
% plot the gaussian fits
tt=linspace(0,2*pi,npts)';
x = cos(tt); y=sin(tt);
ap = [x(:) y(:)]';
[v,d]=eig(C);
d = sdwidth * sqrt(d); % convert variance to sdwidth*sd
bp = (v*d*ap) + repmat(means, 1, size(ap,2));
h = plot(bp(1,:), bp(2,:), '-', 'parent', axh);
set(h,'color',color);

%-----------------------------
function h = show3d(means, C, sdwidth, npts, axh)
if isempty(npts), npts=20; end
[x,y,z] = sphere(npts);
ap = [x(:) y(:) z(:)]';
[v,d]=eig(C);
if any(d(:) < 0)
   fprintf('warning: negative eigenvalues\n');
   d = max(d,0);
end
d = sdwidth * sqrt(d); % convert variance to sdwidth*sd
bp = (v*d*ap) + repmat(means, 1, size(ap,2));
xp = reshape(bp(1,:), size(x));
yp = reshape(bp(2,:), size(y));
zp = reshape(bp(3,:), size(z));
h = surf(axh, xp,yp,zp);

author: epfl-lasa


text:
function G = gsp_stochastic_block_graph(N, k, params)
%GSP_STOCHASTIC_BLOCK_GRAPH  Create a stochastic block graph
%   Usage:  G = gsp_stochastic_block_graph( N );
%           G = gsp_stochastic_block_graph(N , k);
%
%   Input parameters:
%         N     : Number of nodes (default 1024)
%         k     : Number of clusters (default 5)
%         params: Structure of optional parameters
%   Output parameters:
%         G     : Graph structure.
%
%   *param* is an optional structure with the following fields
%
%   * *params.p*        : Intra-cluster edge probability (default 0.7)
%   * *params.q*        : Inter-cluster edge probability (default 0.3/k)
%   * *params.z*        : Assignment vector of nodes (default uniform random)
%   * *params.M*        : Link probability matrix between clusters (default uses p and q)
%   * *params.directed* : Flag the graph as directed or not (default false)
%
%   Use the stochastic block model to create a graph.

% Author: Pierre Vandergheynst, Nathanael Perraudin
% Date  : 2 novemeber 2015 (revision: 6 october 2016 -- Lionel Martin)

%% Stochastic Block Model generator
if nargin<1
    N = 1024; % number of nodes
end
if nargin<2
    k = 5; % number of clusters
end
if nargin < 3
    params = struct;
end

if ~isfield(params, 'p')
    params.p = 0.7;
end

if ~isfield(params, 'q')
    params.q = (1-params.p) / k;
end

if ~isfield(params, 'force_full')
    params.force_full = 0;
end

if ~isfield(params, 'auto_gen_M')
    params.auto_gen_M = 0;
end

if (~isfield(params, 'z') || length(params.z) ~= N || max(params.z) > k || min(params.z) < 1)
    params.z = randi(k, 1, N);
end

if (~isfield(params, 'M') || size(params.M) ~= [k, k])
    if params.force_full
        params.M = params.q * ones(k);
        params.M(1:k+1:end) = params.p;
    end
    params.auto_gen_M = 1;
end

if ~isfield(params, 'directed')
    params.directed = false;
end

%% partition with clusters of homogenous sizes
%
% THIS PART IS NOT NECESSARY BECAUSE SLOW AND REPLACED
% BY THE SIMPLE "UNIFORM" CLUSTER ASSIGNMENT ABOVE
%
% L = ceil(N/k);
% for i=1:k-1
%     z((i-1)*L + 1:i*L) = i;
% end
% z((k-1)*L + 1:end) = k;

%% Generate adjacency
z = params.z;

% for i=1:N
%     for j=i+1:N
%         W(i, j) = ( rand <= M(z(i), z(j)) );
%         if params.directed
%             W(j, i) = ( rand <= M(z(j), z(i)) );
%         else
%             W(j, i) = W(i, j);
%         end
%     end
% end
if params.auto_gen_M && ~params.force_full
    [val_z, idx_z] = sort(z);
    counts = diff(find(diff([0, val_z, N])));

    if length(counts) ~= k
        error('There is at least one empty class. Check your z.');
    end

    W = sprandsym(counts(k), params.p);
    nb_cols_rect = 0;

    for i=k:-1:2
        nb_cols_rect = nb_cols_rect + counts(i);
        rect = sprand(counts(i-1), nb_cols_rect, params.q);
        top_left = sprandsym(counts(i-1), params.p);
        W = vertcat(horzcat(top_left, rect), horzcat(rect', W));
    end

    W(1:N+1:end) = 0;
    G.W(idx_z, idx_z) = abs(W) > 0;

else
    M = params.M;

    if params.directed
        W = rand(N) <= M(z, z);
    else
        A = rand(N);
        A(logical(triu(ones(N)))) = 1;
        W = A <= M(z, z);
        W = W + W';
    end

    G.W = sparse(W);
end

G = gsp_graph_default_parameters(G);
G.info.node_com = z;

G.coords = ones(N, 2);
com_coords = sqrt(N) * [-cos(2*pi*(1:k)/k)', sin(2*pi*(1:k)/k)'];

% create uniformly random points in the unit disc
for ii = 1:N
    % use rejection sampling to sample from a unit disc (probability = pi/4)
    while norm(G.coords(ii, :)) >= 1/2
        % sample from the square and reject anything outside the circle
        G.coords(ii, :) = [rand-.5, rand-.5];
    end
end

% add the offset for each node depending on which community it belongs to
for ii = 1:k
    idx_ii = find(z==ii);
    rad_com = sqrt(numel(idx_ii));
    G.coords(idx_ii, :) = bsxfun(@plus, rad_com * G.coords(idx_ii, :), com_coords(ii, :));
end

end

repo: gspbox


text:
%% Large NLP
clc
clear
%Number of decision vars AND number of linear constraints AND number of
%nonlinear constraints
n = 2000;

% Objective & Gradient
fun = @(x) -sum(x);
grad = @(x) -ones(n,1);

% Linear constraints
A = randn(n,n);
rl = zeros(n,1);
ru = zeros(n,1);
%%
% Nonlinear Constraint, Jacobian & Structure
nlcon = @(x) x(:).^4;
nljac = @(x) sparse(diag(4*x(:).^3));
jacstr = @() speye(n);
hessian = @(x,sigma,lambda) sparse(diag(12*x(:).^2));
hessianstructure = @() speye(n);
cl = -Inf(n,1);
cu = ones(n,1);

% Starting Guess
x0 = 1*ones(n,1);
% Build Function Structure
funcs.objective = fun;
funcs.gradient = grad;
funcs.constraints = nlcon;
funcs.jacobian = nljac;
funcs.jacobianstructure = jacstr;
% funcs.hessian = hessian;
% funcs.hessianstructure = hessianstructure;

% Build Options Structure
lb = -10*ones(n,1);
ub = 40*ones(n,1);
opts = [];
opts.lb = lb;
opts.ub = ub;
opts.rl = rl;
opts.ru = ru;
opts.cl = cl;
opts.cu = cu;
opts.A = sparse(A);
opts.ipopt.hessian_approximation = 'limited-memory';
opts.ipopt.linear_solver = 'pardiso';
% Call IPOPT
[~,output] = ipopt(x0,funcs,opts);

%% Modify to call BONMIN
clc
opts.display = 1;
opts.var_type = zeros(n,1); opts.var_type(1:10) = 1;
[~,output] = bonmin(x0,funcs,opts);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%kkcp 'Copy a toolbox, software object, or file'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kcp.pane file
%
% Parameters:
% Toggle: force 'Force Output', default: 0: 'force output?'
% String: otb 'Old Toolbox', default: ' ': 'origination toolbox'
% String: ntb 'New Toolbox', default: ' ': 'destination toolbox'
% String: oname 'Object Name', default: ' ': 'software object to be copied'
% InputFile: i 'Old File', optional: 'original file'
% OutputFile: o 'New File', optional: 'copied file'
%
% Example: o = kkcp(i, {'force',0;'otb',' ';'ntb',' ';'oname',' ';'i','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kcp - Copy a toolbox, software object, or file
%
%  DESCRIPTION
% When a program or library object is copied from one toolbox to another,
% the directory structure of the program or library object must be duplicated
% in the second toolbox.  Then, all the files associated with the object must be
% copied from the first toolbox to the second.  Note that the relative paths
% from the toolbox top level source are preserved during the copy.  After all
% the files are copied, the new object is registered in the toolbox.
%
% If the object already exists in the destination toolbox,
% "\s-1kcpobj\s+1" will
% prompt whether to cancel the operation, remove the existing object in the
% destination toolbox before re-creating it, or to overwrite the existing object
% in the destination toolbox.
%
% When the [-force] option is used, if the object already exists in the
% destination toolbox, the object will be removed completely and then re-created.
%
% This program allows you to copy your object out of your distribution toolbox
% to a work toolbox so you can work on it.  This method of work is recommended,
% if a program needs to be worked on for an extended period of time.
%
% It is not recommended that the overwrite mode of kcpobj is used after
% major enhancements are made to the object, especially if one or more filenames
% have changed or files have been deleted (because over-write mode will not
% delete files, only add them).
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% SOLVE_KIRCHHOFF_LOVE_SHELL: Solve the Kirchhoff-Love shell model in a NURBS domain.
%
% USAGE:
%
%  [geometry, msh, space, u] = solve_kirchhoff_love_shell (problem_data, method_data)
%
% INPUT:
%
%  problem_data: a structure with data of the problem. It contains the fields:
%    - geo_name:     name of the file containing the geometry
%    - drchlt_sides: sides with Dirichlet boundary condition
%    - drchlt_components: cell-array, the components that are set to zero for each drchlt_side
%    - E_coeff:      function handle for Young's modulus
%    - nu_coeff:     function handle for Poisson's ratio
%    - thickness:    scalar value, thickness of the shell
%    - f:            source term, distributed load
%
%  method_data : a structure with discretization data. Its fields are:
%    - degree:     degree (>=2) of the spline functions.
%    - regularity: continuity (>=1) of the spline functions.
%    - nsub:       number of subelements with respect to the geometry mesh
%                   (nsub=1 leaves the mesh unchanged)
%    - nquad:      number of points for Gaussian quadrature rule
%
% OUTPUT:
%
%  geometry: geometry structure (see geo_load)
%  msh:      mesh object that defines the quadrature rule (see msh_cartesian)
%  space:    space object that defines the discrete basis functions (see sp_vector)
%  u:        the computed degrees of freedom
%
% NOTE: only homogeneous Dirichlet conditions implemented so far
%
% See also EX_KL_SHELL_SCORDELIS_LO_ROOF for an example.
%
% Copyright (C) 2017-2019 Pablo Antolin, Luca Coradello, Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function [geometry, msh, space, u] = solve_kirchhoff_love_shell (problem_data, method_data)

% Extract the fields from the data structures into local variables
data_names = fieldnames (problem_data);
for iopt  = 1:numel (data_names)
  eval ([data_names{iopt} '= problem_data.(data_names{iopt});']);
end
data_names = fieldnames (method_data);
for iopt  = 1:numel (data_names)
  eval ([data_names{iopt} '= method_data.(data_names{iopt});']);
end

if (any(degree <= 1) || any(regularity == 0))
  error ('The degree must be at least two, and the regularity at least C^1')
end

geometry = geo_load (geo_name);
degelev  = max (degree - (geometry.nurbs.order-1), 0);
nurbs    = nrbdegelev (geometry.nurbs, degelev);
[rknots, ~, nknots] = kntrefine (nurbs.knots, nsub-1, nurbs.order-1, regularity);

nurbs    = nrbkntins (nurbs, nknots);
geometry = geo_load (nurbs);

% Construct msh structure
rule     = msh_gauss_nodes (nquad);
[qn, qw] = msh_set_quad_nodes (geometry.nurbs.knots, rule);
msh      = msh_cartesian (geometry.nurbs.knots, qn, qw, geometry);

% Construct space structure
sp_scalar = sp_nurbs (geometry.nurbs, msh);
scalar_spaces = repmat ({sp_scalar}, 1, msh.rdim);
space = sp_vector (scalar_spaces, msh);

% Assemble the stiffness matrix and right-hand side
K = op_KL_shells_tp (space, space, msh, E_coeff, nu_coeff, thickness);
rhs = op_f_v_tp (space, msh, f);

u = zeros (space.ndof, 1);

% Apply boundary conditions
drchlt_dofs = [];
for iside = 1:numel(drchlt_sides)
  side = drchlt_sides(iside);
  if (~exist('drchlt_components','var'))
    components = 1:3;
  else
    components = drchlt_components{iside};
  end
  for icomp = components
    drchlt_dofs = union (drchlt_dofs, space.boundary(side).dofs(space.boundary(side).comp_dofs{icomp}));
  end
end

int_dofs = setdiff (1:space.ndof, drchlt_dofs);
%rhs(int_dofs) = rhs(int_dofs) - K(int_dofs, drchlt_dofs)*u(drchlt_dofs);

% Solve the linear system
u(int_dofs) = K(int_dofs, int_dofs) \ rhs(int_dofs);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
filename='Sphere_Tetrahedra_Linear_Unstructured';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'sphere';
fracRadius = 1-1e-6;
cost = {'compliance'};
weights = [1];
constraint = {'volumeConstraint'};
optimizer = 'HAMILTON-JACOBI';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';

nsteps = 1;
Vfrac_final = 1;
Perimeter_target=3.5;
optimality_final =1e-5;
constr_final =1e-5;

BCscale_factor = 0.3;
HJiter0 = 1;
e2 = 5;

Vfrac_initial = 1;
optimality_initial = 5e-2;
constr_initial = 5e-2;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;

% maxiter = 10;
% maxiter = 1;
TOL.nu_plus = 0.3;
TOL.nu_minus = 0.3;

plotting = 0;
printing = 0;
monitoring = 0;
monitoring_interval = 0;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [prediction,  weight,  RF] = rfMakePrediction(params, id);
% rfMakePrediction - makes time-series based on rf-profile and image
% sequence.
%
% [prediction,  weight,  rf] = rfMakePrediction(params, id);
%
% input:
%  params  : parameter struct
%  id      : id for grid search [n] or rf-parameters
%            [sigmaMajor,  sigmaMinor,  sigmaTheta,  x0, 0]
%
%

% 2006/01 SOD: wrote it.
% 2006/03 SOD: vectorized.
% 2006/06 SOD: moved Hrf convolution out of this function to rmMakeStimulus.
% 2007/04 RAS: revert to version 1.15 in repository; this restores the
% making of the RF, based on the id params. It seems that removing that
% part made the function obsolete (it was just a multiplication), and also
% broke the RF visualization code, so for now I'm restoring it.

% Programming note: This function is at the heart of the retinotopy
% model program and will be called lots. It is now vectorized so we
% can make all RFs in advance.

% Make (several) receptive field profiles.
% Ugly input check allowing different kinds of inputs.
if numel(id)==1,
    % if id == 0 then we make all the RFs
    if ~id,
        RF =  rfGaussian2d(params.analysis.X, params.analysis.Y, ...
                            params.analysis.sigmaMajor, ...
                            params.analysis.sigmaMajor, ...
                            0, ...
                            params.analysis.x0,  ...
                            params.analysis.y0);

    else,  % we make just a particular one
        RF =  rfGaussian2d(params.analysis.X, params.analysis.Y, ...
                            params.analysis.sigmaMajor(id), ...
                            params.analysis.sigmaMajor(id), ...
                            0, ...
                            params.analysis.x0(id),  ...
                            params.analysis.y0(id));

    end;

elseif numel(id)==5,
    RF =  rfGaussian2d(params.analysis.X, params.analysis.Y, ...
        id(1), id(2), id(3), id(4), id(5));

else,  % make lots
    RF =  rfGaussian2d(params.analysis.X, params.analysis.Y, ...
        id(:, 1), id(:, 2), id(:, 3), id(:, 4), id(:, 5));
end;

% Now we have to loop over each stimulus to convolve it with the RF. The
% stimulus is already convolved with the Hrf so we don't have to do that
% here,  also removing of initial time frames and time averaging is
% moved to the stimulus creation.
prediction = params.analysis.allstimimages * RF;

% compute the amount within the stimulus window to weight the
% fit
if nargout > 1,
    weight = sum(RF);
end;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% WAVELET_3D_PYRAMID Compute the roto-translation wavelet transfrom
%
% Usage
%   [y_Phi, y_Psi, meta_Phi, meta_Psi] = wavelet_3d_pyramid(y, filters, filters_rot, options)
%
% Input
%   y (numeric): a 3d matrix whose first two dimension corresponds to spatial
%       postion and third dimension corresponds to orientation.
%   filters (struct): a 2d pyramid filter bank (applied along spatial variable)
%   filters_rot (struct): a 1d filter bank (applied along orientation)
%   options (struct): containing the following optional fields
%       x_resolution (int): the log of spatial resolution
%       psi_mask (bool array): a mask for determining which filter to apply
%       oversampling (int): the log of spatial oversampling
%       oversampling_rot (int): the log of orientation oversampling
%
% Output
%   y_Phi (numeric): the roto-translation convolution y * Phi
%   y_Psi (cell): containing all the roto-translation convolution y * Psi
%   meta_phi (struct): meta associated to y_Phi
%   meta_psi (struct): meta associated to y_Psi
%
% Description
%   This is a pyramid implementation of the roto-translation wavelet
%   transform. Another (slower, FFT-based) implementation is available as
%   WAVELET_3D.
%	This function computes the roto-translation convolution of a three dimensional
%	signal y, with roto-translation wavelets defined as the separable product
%	low pass :
%		PHI(U,V) * PHI(THETA)
%	high pass :
%		PHI(U,V) * PHI(THETA)
%		PSI(U,V) * PHI(THETA)
%		PSI(U,V) * PSI(THETA)
%
% Reference
%	Rotation, Scaling and Deformation Invariant Scattering for Texture
%	Discrimination, Laurent Sifre, Stephane Mallat
%	Proc of CVPR 2013
%	http://www.cmapx.polytechnique.fr/~sifre/research/cvpr_13_sifre_mallat_final.pdf
%
% See also
%   WAVELET_LAYER_3D, WAVELET_FACTORY_3D

function [y_Phi, y_Psi, meta_Phi, meta_Psi] = wavelet_3d_pyramid(y,...
        filters, filters_rot, options)

    % retrieve info on input
    Q = filters.meta.Q;
    L = filters_rot.meta.size_filter / 2;
    % do not compute any convolution with high pass if y_Psi is not
    % outputed
    calculate_psi = (nargout>=2);
    [N, M, T] = size(y);

    % retrieve options
    white_list = {'J', 'angular_range', 'oversampling_rot', 'j_min', 'q_mask'};
    check_options_white_list(options, white_list);
    options = fill_struct(options, 'J', 4);
    options = fill_struct(options, 'angular_range', 'non_specified');
    options = fill_struct(options, 'oversampling_rot', -1);
    options = fill_struct(options, 'j_min', 0);
    options = fill_struct(options, 'q_mask', ones(1, Q));

    % angular resolution
    switch (options.angular_range)
        case 'zero_pi'
            angular_res = floor(log2(2*L/(2*T)));
        case 'zero_2pi'
            angular_res = floor(log2(2*L/T));
        otherwise
            error('unspecified angular range');
    end

    if (isa(filters.h.filter.coefft, 'single'))
        cast = @single;
    else
        cast = @double;
    end

    %%%%%%%%%%%%%%%%%%
    %%% low passes %%%
    %%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%
    %%% phi - phi %%%%
    %%%%%%%%%%%%%%%%%%

    %% low pass filter spatial : cascade with h for every slice
    % initialize structure
    hy.signal{1} = y;
    hy.meta.j(1) = 0;

    % low pass spatial using a pyramid
    h = filters.h.filter;
    for j = 1:options.J
        for theta = 1:T % filter each slice
            if (theta == 1) % allocate when size is known
                tmp_slice = pad_conv_sub_unpad(hy.signal{j}(:,:,theta), h);
                tmp = cast(zeros([size(tmp_slice), T]));
                tmp(:,:,theta) = tmp_slice;
            else
                clear tmp_slice;
                tmp(:,:,theta) = ...
                    pad_conv_sub_unpad(hy.signal{j}(:,:,theta), h);
            end
        end
        hy.signal{j+1} = tmp;
        hy.meta.j(j+1) = j;
    end
    y_phi.signal{1} = hy.signal{options.J+1};
    y_phi.meta.j(1) = hy.meta.j(options.J+1);

    %% low pass filter anguler : fft along anguler
    % initialize structure (recopy if range is zero_pi)
    ds = max(filters_rot.meta.J/filters_rot.meta.Q - options.oversampling_rot, 0);
    if (2^ds == 2*L)
        % in this case it is faster to compute the sum along the angle
        y_Phi = sum(y_phi.signal{1},3) / 2^(ds/2);
    else
        if (strcmp(options.angular_range, 'zero_pi'))
            % divide by sqrt(2) for energy preservation
            y_phi.signal{1} = repmat(y_phi.signal{1}, [1 1 2]) / sqrt(2);
        end
        y_phi_f = fft(y_phi.signal{1}, [], 3);
        phi_angle = filters_rot.phi.filter;
        y_Phi = real(...
            sub_conv_1d_along_third_dim_simple(y_phi_f, phi_angle, ds));
    end
    meta_Phi.j2 = options.J;
    meta_Phi.k2 = filters_rot.meta.J;

    %%

    %%%%%%%%%%%%%%%%%%%
    %%% high passes %%%
    %%%%%%%%%%%%%%%%%%%

    if (calculate_psi)
        p = 1;
        %%%%%%%%%%%%%%%%%%
        %%% phi - psi %%%%
        %%%%%%%%%%%%%%%%%%

        %% low pass spatial - already computed
        %% high pass angular
        % NOTE : y_phi_f
        % (the (angular) fourier transform of y * (spatial) phi)
        % might have already been computed
        if ~exist('y_phi_f', 'var')
            % fourier angle
            if (strcmp(options.angular_range, 'zero_pi'))
                % divide by sqrt(2) for energy preservation
                y_phi.signal{1} = repmat(y_phi.signal{1}, [1 1 2]) / sqrt(2);
            end
            y_phi_f = fft(y_phi.signal{1}, [], 3);
        end
        for k2 = 0:numel(filters_rot.psi.filter)-1
            psi_angle = filters_rot.psi.filter{k2+1};
            ds = max(k2/filters_rot.meta.Q - options.oversampling_rot, 0);
            y_Psi{p} = ...
                sub_conv_1d_along_third_dim_simple(y_phi_f, psi_angle, ds);
            meta_Psi.j2(p) = options.J;
            meta_Psi.q2(p) = 0;
            meta_Psi.theta2(p) = 0;
            meta_Psi.k2(p) = k2;
            p = p + 1;
        end

        %%%%%%%%%%%%%%%%%%
        %%% psi - phi %%%%
        %%%%%%%%%%%%%%%%%%
        %      AND       %
        %%%%%%%%%%%%%%%%%%
        %%% psi - phi %%%%
        %%%%%%%%%%%%%%%%%%

        %% high pass spatial - with pyramid
        if (angular_res == 0)
            g = filters.g.filter;

            for j2 = options.j_min:options.J-1
                for q = find(options.q_mask==1)-1
                    for theta2 = 1:L
                        if (strcmp(options.angular_range, 'zero_pi'))
                            for theta = 1:L
                                % convolution with psi_{j1, theta + theta2}
                                theta_sum_mod2L =  1 + mod(theta + theta2 - 2, 2*L);
                                theta_sum_modL =  1 + mod(theta + theta2 - 2, L);
                                tmp_slice = pad_conv_unpad(hy.signal{j2+1}(:,:,theta), g{theta_sum_modL + L*q});
                                if (theta == 1) % allocate
                                    tmp = cast(zeros([size(tmp_slice), 2*L]));
                                end
                                if (theta_sum_mod2L <= L)
                                    tmp(:,:,theta) = tmp_slice;
                                    tmp(:,:,theta+L) = conj(tmp_slice);
                                else
                                    tmp(:,:,theta) = conj(tmp_slice);
                                    tmp(:,:,theta+L) = tmp_slice;
                                end
                            end
                        else % options.angular_range is 'zero_2pi'
                            for theta = 1:2*L
                                % convolution with psi_{j1, theta + theta2}
                                theta_sum_mod2L =  1 + mod(theta + theta2 - 2, 2*L);
                                theta_sum_modL =  1 + mod(theta + theta2 - 2, L);
                                if (theta_sum_mod2L <= L)
                                    curr_g =  g{theta_sum_modL + L*q};
                                else
                                    curr_g =  conjugate_filter(g{theta_sum_modL + L*q});
                                end
                                tmp_slice = pad_conv_unpad(hy.signal{j2+1}(:,:,theta), curr_g);
                                if (theta == 1) % allocate
                                    tmp = cast(zeros([size(tmp_slice), 2*L]));
                                end
                                tmp(:,:,theta) = tmp_slice;

                            end
                        end

                        % tmp can now be filtered along the orientation
                        tmp_f = fft(tmp, [], 3);
                        %% low pass angle
                        ds = min(2*L, max(filters_rot.meta.J/filters_rot.meta.Q - options.oversampling_rot, 0));
                        if (2^ds == 2*L)
                            % faster to compute the sum along the angle
                            y_Psi{p} = sum(tmp, 3) / 2^(ds/2);
                        else
                            y_Psi{p} = ...
                                sub_conv_1d_along_third_dim_simple(tmp_f, phi_angle, ds);
                        end
                        meta_Psi.j2(p) = j2;
                        meta_Psi.theta2(p) = theta2;
                        meta_Psi.k2(p) = filters_rot.meta.J;
                        meta_Psi.q2(p) = q;
                        p = p + 1;

                        %% high pass angle
                        for k2 = 0:numel(filters_rot.psi.filter)-1
                            psi_angle = filters_rot.psi.filter{k2+1};
                            ds = max(k2/filters_rot.meta.Q - options.oversampling_rot, 0);
                            y_Psi{p} = ...
                                sub_conv_1d_along_third_dim_simple(tmp_f, psi_angle, ds);
                            meta_Psi.j2(p) = j2;
                            meta_Psi.theta2(p) = theta2;
                            meta_Psi.k2(p) = k2;
                            meta_Psi.q2(p) = q;
                            p = p + 1;
                        end

                    end
                end
            end
        else
            error('not yet supported');
        end

    end

    function out = pad_conv_sub_unpad(in, filter)
        Npad = size(in) + filters.meta.size_filter - 1;
        out = pad_signal(in, Npad, 'symm', 1);
        out = conv_sub_2d(out, filter, 1);
        out = unpad_signal(out, 1, size(in), filters.meta.offset);
    end

    function out = pad_conv_unpad(in, filter)
        Npad = size(in) + filters.meta.size_filter - 1;
        out = pad_signal(in, Npad, 'symm', 1);
        out = conv_sub_2d(out, filter, 0);
        out = unpad_signal(out, 0, size(in), filters.meta.offset);
    end

    function conj_filter = conjugate_filter(filter)
        conj_filter = filter;
        conj_filter.coefft = conj(filter.coefft);
    end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [v,usediters] = projfunc( s, k1, k2, nn )

% Solves the following problem:
% Given a vector s, find the vector v having sum(abs(v))=k1
% and sum(v.^2)=k2 which is closest to s in the euclidian sense.
% If the binary flag nn is set, the vector v is additionally
% restricted to being non-negative (v>=0).
%
% Written 2.7.2004 by Patrik O. Hoyer
%

% Problem dimension
N = length(s);

% If non-negativity flag not set, record signs and take abs
if ~nn,
    isneg = s<0;
    s = abs(s);
end

% Start by projecting the point to the sum constraint hyperplane
v = s + (k1-sum(s))/N;

% Initialize zerocoeff (initially, no elements are assumed zero)
zerocoeff = [];

j = 0;
while 1,

    % This does the proposed projection operator
    midpoint = ones(N,1)*k1/(N-length(zerocoeff));
    midpoint(zerocoeff) = 0;
    w = v-midpoint;
    a = sum(w.^2);
    b = 2*w'*v;
    c = sum(v.^2)-k2;
    alphap = (-b+real(sqrt(b^2-4*a*c)))/(2*a);
    v = alphap*w + v;

    if all(v>=0),
	% We've found our solution
	usediters = j+1;
	break;
    end

    j = j+1;

    % Set negs to zero, subtract appropriate amount from rest
    zerocoeff = find(v<=0);
    v(zerocoeff) = 0;
    tempsum = sum(v);
    v = v + (k1-tempsum)/(N-length(zerocoeff));
    v(zerocoeff) = 0;

end

% If non-negativity flag not set, return signs to solution
if ~nn,
    v = (-2*isneg + 1).*v;
end

% Check for problems
if max(max(abs(imag(v))))>1e-10,
    error('Somehow got imaginary values!');
end


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

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

% Last Modified by GUIDE v2.5 03-Nov-2016 22:15:06

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @camera_pipeline_OpeningFcn, ...
                   'gui_OutputFcn',  @camera_pipeline_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

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

% Choose default command line output for camera_pipeline
handles.output = hObject;
%TODO
% maximize the GUI
set(0, 'Unit', 'pixel');
winsize = get(0, 'ScreenSize');

% Set the Gui position
winsize(1) = winsize(1) + 100;
winsize(2) = winsize(2) + 100;
winsize(3) = winsize(3) - 150;
winsize(4) = winsize(4) - 150;
set(gcf, 'Position', winsize);

step = winsize(4) / 16;
% Set the Panel position
set(handles.uipanel1, 'Position', [3/4* winsize(3) + 100, 10, 1/4*winsize(3) - 20, winsize(4) + 50]);
% Set the Axes position disappear
set(handles.axes1, 'Position', [-10, -10, 5, 5]);
% Update handles structure
handles.myData.winsize = winsize;

%current stage in pipeline
handles.myData.current_stage = zeros(1,1);

pipeline_image = im2double(imread(fullfile('.','Adobe_DNG_SDK_Pipeline_v4.bmp')));
% set the slider range and step size
numSteps = 11;
set(handles.slider1, 'Min', 1);
set(handles.slider1, 'Max', numSteps);
set(handles.slider1, 'Value', 1);
set(handles.slider1, 'SliderStep', [1/(numSteps-1) , 1/(numSteps-1) ]);
% save the current/last slider value
handles.lastSliderVal = get(handles.slider1,'Value');

 % Update handles structure
handles.pushbutton1Pressed = 0;

guidata(hObject, handles);
% UIWAIT makes camera_pipeline wait for user response (see UIRESUME)
% uiwait(handles.figure1);
showImage(handles, pipeline_image);
% prompt = {'Enter matrix size:','Enter colormap name:'};
% dlg_title = 'Input';
% num_lines = 1;
% defaultans = {'20','hsv'};
% answer = inputdlg(prompt,dlg_title,num_lines,defaultans);

% Get all the handles to everything we want to set in a single array.
handleArray = [handles.pushbutton1, handles.pushbutton12, handles.radiobutton3, handles.radiobutton8, handles.radiobutton15, handles.radiobutton16];
% Set them all disabled.
set(handleArray, 'Enable', 'off');

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

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

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
%%
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[filename, pathname] = uigetfile( ...
{'*.*',  'All Files (*.*)';
'*.cr2;*.nef;*.pef;','RAW Image Files (*.cr2,*.nef, *.pef)';
'*.dat','Spectral data files (*.dat)'}, ...
   'Pick a file');
if isequal(filename,0) || isequal(pathname,0)
else
    handles.myData.fileName = [pathname filename];

    a = [get(handles.radiobutton15,'Value') get(handles.radiobutton16,'Value') get(handles.radiobutton3,'Value') get(handles.radiobutton8,'Value')];
    fileID = fopen('wbAndGainSettings.txt','w');
    fprintf(fileID,'%d\n',a);
    fclose(fileID);

    rwSettings = [0];
    if (get(handles.radiobutton17,'Value') == 1)
        rwSettings(1) = 1;
    end
    if ((get(handles.radiobutton17,'Value') == 0) && (get(handles.radiobutton18,'Value') == 0))
        rwSettings(1) = 2;
    end
    fileID = fopen('rwSettings.txt','w');
    fprintf(fileID,'%d\n',rwSettings);
    fclose(fileID);

    cam_settings = [0];
    if (get(handles.radiobutton20,'Value') == 1)
        cam_settings(1) = 1;
    end
    if ((get(handles.radiobutton20,'Value') == 0) && (get(handles.radiobutton19,'Value') == 0))
        cam_settings(1) = 2;
    end
    fileID = fopen('cam_settings.txt','w');
    fprintf(fileID,'%d\n',cam_settings);
    fclose(fileID);

    my_waitbar = waitbar(0, 'Rendering Image...');

    inputFileName = [pathname filename];
    if (handles.lastSliderVal == 1)
        stage1OutputFileName = [inputFileName(1:end-4) '_stage1.tif'];
        system_command = [fullfile('.','dngOneExeSDK','dng_validate.exe -16 -1 ') stage1OutputFileName ' ' inputFileName];
        system(system_command);
        handles.myData.image = imread(stage1OutputFileName);
        %imtool(uint8(handles.myData.image/256))
        delete(stage1OutputFileName);
        delete(my_waitbar)       % DELETE the waitbar; don't try to CLOSE it.
    end

    if (handles.lastSliderVal == 2)

        a = [0 1 0 1]; % without gain [0 1 ...
        fileID = fopen('wbAndGainSettings.txt','w');
        fprintf(fileID,'%d\n',a);
        fclose(fileID);

        stageSettings_indicator = [2];
        fileID = fopen('stageSettings.txt','w');
        fprintf(fileID,'%d\n',stageSettings_indicator);
        fclose(fileID);

        if (get(handles.radiobutton17,'Value') == 0)
            fileID = fopen('lastStage.txt','w');
            fprintf(fileID,'%d\n',stageSettings_indicator);
            fclose(fileID);
        end

        stage2OutputFileName = [inputFileName(1:end-4) '_stage2.tif'];
        system_command = [fullfile('.','dngOneExeSDK','dng_validate.exe -16 -2 ')  stage2OutputFileName ' ' inputFileName];
        system(system_command);
        handles.myData.image = imread(stage2OutputFileName);
        delete(stage2OutputFileName);
        delete(my_waitbar)       % DELETE the waitbar; don't try to CLOSE it.
    end

    if (handles.lastSliderVal == 3)

        a = [1 0 0 1]; % with gain [1 0 ...
        fileID = fopen('wbAndGainSettings.txt','w');
        fprintf(fileID,'%d\n',a);
        fclose(fileID);

        stageSettings_indicator = [3];
        fileID = fopen('stageSettings.txt','w');
        fprintf(fileID,'%d\n',stageSettings_indicator);
        fclose(fileID);

        if (get(handles.radiobutton17,'Value') == 0)
            fileID = fopen('lastStage.txt','w');
            fprintf(fileID,'%d\n',stageSettings_indicator);
            fclose(fileID);
        end

        stage3OutputFileName = [inputFileName(1:end-4) '_stage3.tif'];
        system_command = [fullfile('.','dngOneExeSDK','dng_validate.exe -16 -2 ')  stage3OutputFileName ' ' inputFileName];
        system(system_command);
        handles.myData.image = imread(stage3OutputFileName);
        delete(stage3OutputFileName);
        delete(my_waitbar)       % DELETE the waitbar; don't try to CLOSE it.
    end

    if (handles.lastSliderVal == 4)
        a = [4];
        fileID = fopen('stageSettings.txt','w');
        fprintf(fileID,'%d\n',a);
        fclose(fileID);

        stage4OutputFileName = [inputFileName(1:end-4) '_stage4.tif'];

        if (get(handles.radiobutton17,'Value') == 0 && get(handles.radiobutton18,'Value') == 1 )
            fileID = fopen('lastStage.txt','w');
            fprintf(fileID,'%d\n',a);
            fclose(fileID);
            system_command = [fullfile('.','dngOneExeSDK','dng_validate.exe -16 -cs1 -tif ')  stage4OutputFileName ' ' inputFileName];
            system(system_command);
            delete(stage4OutputFileName);
        end

        system_command = [fullfile('.','dngOneExeSDK','dng_validate.exe -16 -3 ')  stage4OutputFileName ' ' inputFileName];
        system(system_command);
        handles.myData.image = imread(stage4OutputFileName);
        delete(stage4OutputFileName);
        delete(my_waitbar)       % DELETE the waitbar; don't try to CLOSE it.
    end

    if (handles.lastSliderVal == 5)

        a = [5];
        fileID = fopen('stageSettings.txt','w');
        fprintf(fileID,'%d\n',a);
        fclose(fileID);

        if (get(handles.radiobutton17,'Value') == 0)
            fileID = fopen('lastStage.txt','w');
            fprintf(fileID,'%d\n',a);
            fclose(fileID);
        end
        %
        stage5OutputFileName = [inputFileName(1:end-4) '_stage5.tif'];
        system_command = [fullfile('.','dngOneExeSDK','dng_validate.exe -16 -cs1 -tif ')  stage5OutputFileName ' ' inputFileName];
        system(system_command);
        handles.myData.image = imread(stage5OutputFileName);
        delete(stage5OutputFileName);
        delete(my_waitbar)       % DELETE the waitbar; don't try to CLOSE it.
    end

    if (handles.lastSliderVal == 6)

        a = [6];
        fileID = fopen('stageSettings.txt','w');
        fprintf(fileID,'%d\n',a);
        fclose(fileID);

        if (get(handles.radiobutton17,'Value') == 0)
            fileID = fopen('lastStage.txt','w');
            fprintf(fileID,'%d\n',a);
            fclose(fileID);
        end

        stage6OutputFileName = [inputFileName(1:end-4) '_stage6.tif'];
        system_command = [fullfile('.','dngOneExeSDK','dng_validate.exe -16 -cs1 -tif ')  stage6OutputFileName ' ' inputFileName];
        system(system_command);
        handles.myData.image = imread(stage6OutputFileName);
        delete(stage6OutputFileName);
        delete(my_waitbar)       % DELETE the waitbar; don't try to CLOSE it.
    end

    if (handles.lastSliderVal == 7)

        a = [7];
        fileID = fopen('stageSettings.txt','w');
        fprintf(fileID,'%d\n',a);
        fclose(fileID);

        if (get(handles.radiobutton17,'Value') == 0)
            fileID = fopen('lastStage.txt','w');
            fprintf(fileID,'%d\n',a);
            fclose(fileID);
        end

        stage7OutputFileName = [inputFileName(1:end-4) '_stage7.tif'];
        system_command = [fullfile('.','dngOneExeSDK','dng_validate.exe -16 -cs1 -tif ') stage7OutputFileName ' ' inputFileName];
        system(system_command);
        handles.myData.image = imread(stage7OutputFileName);
        delete(stage7OutputFileName);
        delete(my_waitbar)       % DELETE the waitbar; don't try to CLOSE it.
    end

    if (handles.lastSliderVal == 8)

        a = [8];
        fileID = fopen('stageSettings.txt','w');
        fprintf(fileID,'%d\n',a);
        fclose(fileID);

        if (get(handles.radiobutton17,'Value') == 0)
            fileID = fopen('lastStage.txt','w');
            fprintf(fileID,'%d\n',a);
            fclose(fileID);
        end

        stage8OutputFileName = [inputFileName(1:end-4) '_stage8.tif'];
        system_command = [fullfile('.','dngOneExeSDK','dng_validate.exe -16 -cs1 -tif ') stage8OutputFileName ' ' inputFileName];
        system(system_command);
        handles.myData.image = imread(stage8OutputFileName);
        delete(stage8OutputFileName);
        delete(my_waitbar)       % DELETE the waitbar; don't try to CLOSE it.
    end

    if (handles.lastSliderVal == 9)

        a = [9];
        fileID = fopen('stageSettings.txt','w');
        fprintf(fileID,'%d\n',a);
        fclose(fileID);

        if (get(handles.radiobutton17,'Value') == 0)
            fileID = fopen('lastStage.txt','w');
            fprintf(fileID,'%d\n',a);
            fclose(fileID);
        end

        stage9OutputFileName = [inputFileName(1:end-4) '_stage9.tif'];
        system_command = [fullfile('.','dngOneExeSDK','dng_validate.exe -16 -cs1 -tif ') stage9OutputFileName ' ' inputFileName];
        system(system_command);
        handles.myData.image = imread(stage9OutputFileName);
        delete(stage9OutputFileName);
        delete(my_waitbar)       % DELETE the waitbar; don't try to CLOSE it.
    end

    if (handles.lastSliderVal == 10)

        a = [10];
        fileID = fopen('stageSettings.txt','w');
        fprintf(fileID,'%d\n',a);
        fclose(fileID);

        if (get(handles.radiobutton17,'Value') == 0)
            fileID = fopen('lastStage.txt','w');
            fprintf(fileID,'%d\n',a);
            fclose(fileID);
        end

        stage10OutputFileName = [inputFileName(1:end-4) '_stage10.tif'];%comment here when dng save tests.
        system_command = [fullfile('.','dngOneExeSDK','dng_validate.exe -16 -cs1 -tif ') stage10OutputFileName ' ' inputFileName];%comment here when dng save tests.
%         stage10OutputFileName = [inputFileName(1:end-4) '_stage10_V3.dng']; %comment out here when dng save tests.
%         system_command = [fullfile('.','dngOneExeSDK','dng_validate.exe -dng ') stage10OutputFileName ' ' inputFileName]; %comment out here when dng save tests.

        system(system_command);
        handles.myData.image = imread(stage10OutputFileName);%comment here when dng save tests.
        delete(stage10OutputFileName);        %comment here when dng save tests.
        delete(my_waitbar)       % DELETE the waitbar; don't try to CLOSE it.
    end

    if (handles.lastSliderVal == 11)

        a = [11];
        fileID = fopen('stageSettings.txt','w');
        fprintf(fileID,'%d\n',a);
        fclose(fileID);

        if (get(handles.radiobutton17,'Value') == 0)
            fileID = fopen('lastStage.txt','w');
            fprintf(fileID,'%d\n',a);
            fclose(fileID);
        end

        stage11OutputFileName = [inputFileName(1:end-4) '_stage11.tif'];%comment here when dng save tests.
        system_command = [fullfile('.','dngOneExeSDK','dng_validate.exe -16 -cs1 -tif ') stage11OutputFileName ' ' inputFileName];%comment here when dng save tests.
%         stage11OutputFileName = [inputFileName(1:end-4) '_stage11_V3.dng']; %comment out here when dng save tests.
%         system_command = [fullfile('.','dngOneExeSDK','dng_validate.exe -dng ') stage11OutputFileName ' ' inputFileName]; %comment out here when dng save tests.
        system(system_command);
        handles.myData.image = imread(stage11OutputFileName);
        delete(stage11OutputFileName);   %comment here when dng save tests.
        delete(my_waitbar)       % DELETE the waitbar; don't try to CLOSE it.
    end

%     if (handles.lastSliderVal <= 9)
%         set(handles.pushbutton3, 'Enable', 'on');
%     end
    set(handles.pushbutton1, 'Enable', 'off');
    set(handles.pushbutton12, 'Enable', 'on');
    set(handles.slider1, 'Enable', 'off');
    % Update handles structure
    guidata(hObject, handles);
    if (handles.lastSliderVal >= 4)
        showImage(handles, uint16(handles.myData.image));
    else
        showImage(handles, handles.myData.image);
    end
end

function showImage(handles, I)
% I = handles.myData.image;
winsize = handles.myData.winsize;
[wi, hi, ~] = size(I);
ws = winsize(4);
hs = 3/4*winsize(3);
if (ws/hs > wi/hi)
    ws = hs * (wi/hi);
else
    hs = ws * (hi/wi);
end

set(handles.axes1, 'Position', [10, 100, hs - 20, ws - 20]);
set(handles.pushbutton1, 'Parent', handles.uipanel1);
axes(handles.axes1);
% I(:,:,1) = I(:,:,1) fullfile('.', max(max(I(:,:,1)));
% I(:,:,2) = I(:,:,2) fullfile('.', max(max(I(:,:,2)));
% I(:,:,3) = I(:,:,3) fullfile('.', max(max(I(:,:,3)));
imshow(I);
% guidata(hObject, handles);

% --- Executes on button press in radiobutton3.
function radiobutton3_Callback(hObject, eventdata, handles)
%%
% hObject    handle to radiobutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in pushbutton12.
function pushbutton12_Callback(hObject, eventdata, handles)
%%
% hObject    handle to pushbutton12 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

saved_image =handles.myData.image;

% outputFileName = 'current_result.dng';
outputFileName = 'current_result.tif';
t = Tiff(outputFileName,'w');
% 8 bit save works well
output_unit16 = saved_image;

sliderValue =  handles.lastSliderVal;

tagstruct.ImageLength = size(output_unit16,1);
tagstruct.ImageWidth = size(output_unit16,2);
tagstruct.BitsPerSample = 16;
if sliderValue == 1 || sliderValue == 2 || sliderValue == 3
    tagstruct.SamplesPerPixel = 1;
    tagstruct.Photometric = Tiff.Photometric.MinIsBlack;
else
    tagstruct.SamplesPerPixel = 3;
    tagstruct.Photometric = Tiff.Photometric.RGB;
end
% tagstruct.RowsPerStrip = 16;
tagstruct.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky;
tagstruct.Software = 'MATLAB';
% tagstruct.DNGVersion = [1; 4; 0; 0];
t.setTag(tagstruct);

t.write(output_unit16);
t.close();

% 16 bit save does not work well
% output_unit16 = uint16(saved_image*255);
%
% tagstruct.ImageLength = size(output_unit16,1);
% tagstruct.ImageWidth = size(output_unit16,2);
% tagstruct.Photometric = Tiff.Photometric.RGB;
% tagstruct.BitsPerSample = 16;
% tagstruct.SamplesPerPixel = 3;
% tagstruct.RowsPerStrip = 8; %16 and 32 are also tried. cannot save 16
% bit.
% tagstruct.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky;
% tagstruct.Software = 'MATLAB';
% t.setTag(tagstruct);
%
% t.write(output_unit16);
% t.close();

% --- Executes on button press in radiobutton8.
function radiobutton8_Callback(hObject, eventdata, handles)
%%
% hObject    handle to radiobutton8 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton8

% --- Executes on slider movement.
function slider1_Callback(hObject, eventdata, handles)
% hObject    handle to slider1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

% set(handles.edit32,'String',num2str( get( get(handles.slider1,'Value'))));
%

  handles = guidata(hObject);
 % get the slider value and convert it to the nearest integer that is less
 % than this value
 newVal = floor(get(hObject,'Value'));
 % set the slider value to this integer which will be in the set {1,2,3,...,12,13}
 set(hObject,'Value',newVal);
 % now only do something in response to the slider movement if the
 % new value is different from the last slider value
 if newVal ~= handles.lastSliderVal
     % it is different, so we have moved up or down from the previous integer
     % save the new value
     handles.lastSliderVal = newVal;
     guidata(hObject,handles);
    % display the current value of the slider
    disp(['at slider value ' num2str(get(hObject,'Value'))]);
 end

 sliderValue =  handles.lastSliderVal;
 set(handles.text3,'String',num2str(sliderValue));

if sliderValue >= 4
    % Get all the handles to everything we want to set in a single array.
    handleArray = [handles.radiobutton15, handles.radiobutton16];
    % Set them all disabled.
    set(handleArray, 'Enable', 'on');
    if sliderValue >= 5
        % Get all the handles to everything we want to set in a single array.
        handleArray = [handles.radiobutton3, handles.radiobutton8, handles.radiobutton15, handles.radiobutton16];
        % Set them all disabled.
        set(handleArray, 'Enable', 'on');
    end
end
if sliderValue < 4
    % Get all the handles to everything we want to set in a single array.
    handleArray = [handles.radiobutton3, handles.radiobutton8, handles.radiobutton15, handles.radiobutton16];
    % Set them all disabled.
    set(handleArray, 'Enable', 'off');
end
if sliderValue < 5
    % Get all the handles to everything we want to set in a single array.
    handleArray = [handles.radiobutton3, handles.radiobutton8];
    % Set them all disabled.
    set(handleArray, 'Enable', 'off');
end

% if sliderValue >= 4
%     % Set them all disabled.
%     set(handles.pushbutton3, 'Enable', 'off');
% end
if sliderValue >= 1
    % Set them all disabled.
    set(handles.pushbutton1, 'Enable', 'on');
end
% --- Executes during object creation, after setting all properties.
function slider1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slider1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on button press in radiobutton16.
function radiobutton16_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton16 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton16

% --- Executes on button press in radiobutton15.
function radiobutton15_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton15 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton15

% --- Executes on button press in radiobutton18.
function radiobutton18_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton18 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton18

% --- Executes on button press in radiobutton17.
function radiobutton17_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton17 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton17

% --- Executes on button press in radiobutton20.
function radiobutton20_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton20 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton20

% --- Executes on button press in radiobutton19.
function radiobutton19_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton19 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton19

author: karaimer

repo: camera-pipeline-UI


text:
function rankVal=rankNestedParenth(pString)
%%UNRANKNESTEDPARENTH Obtain the lexicographic order rank of a set of
%           nested parentheses, where in each pair, '(' always comes before
%           ')' for it to be valid, but ')' is considered lexicographically
%           smaller than ')'.
%
%INPUTS: pString A string or character array holding a valid set of nested
%                parenthese, such as pString='((()()()))'. Invalid string
%                will either cause an error or an incorrect result to be
%                produced.
%
%OUTPUTS: rankVal The lexicographic rank of the string, string from 1.
%
%This function implements the modification of Algorithm U of Section
%7.2.1.6 of [1] that is given in Problem 50 in Section 7.2.1.6.
%
%REFERENCES
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4A: Combinatorial
%    Algorithms, Part I, Boston: Addison-Wesley, 2011.
%
%October 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(isempty(pString))
    rankVal=[];
    return
end

pString=char(pString);
n=length(pString)/2;

q=n;
m=1;
p=1;
c=1;
rankVal=1;

while(p<n)
    p=p+1;
    c=((4*p-2)*c)/(p+1);
end

while(q~=0)
    while(1)
        if(m>2*n)
            error('Invalid parenthesis character string provided.')
        end

        cPrime=((q+1)*(q-p)*c)/((q+p)*(q-p+1));
        if(pString(m)==')')
            q=q-1;
            c=cPrime;
            m=m+1;
            break;
        end

        p=p-1;
        c=c-cPrime;
        rankVal=rankVal+cPrime;
        m=m+1;
    end
end
if(m~=2*n+1)
    error('Invalid parenthesis character string provided.')
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.


text:
function r=rot(psi,theta,phi);
% function r=rot(psi,theta,phi)
%
% Description
%  Computes the rotarion matrix corresponding to a oritentation vector
% The input data are:
% - psi   .- the psi angle (radians)
% - theta .- the theta angle (radians)
% - phi   .- the phi angle (radians)
%
% The return value is
% the rotation matrix (dim 3x3)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%               Grupo de Robotica                                          %%
%%               Departamento de Informatica e Ingenieria de Sistemas       %%
%%               Universidad de Zaragoza                                    %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Name         : .m                                              %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Programmer   : Jose Maria Martinez                                       %%
%% Languaje     : Matlab 4.2 c                                              %%
%% Date         : February 1996                                             %%
%% Status       : Prueba                                                    %%
%% History	: 16-8-95 creation                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cosphi=cos(phi);
sinphi=sin(phi);
costheta=cos(theta);
sintheta=sin(theta);
cospsi=cos(psi);
sinpsi=sin(psi);

r=zeros(3,3);
r(1,1)=cosphi*costheta;
r(1,2)=cosphi*sintheta*sinpsi-sinphi*cospsi;
r(1,3)=cosphi*sintheta*cospsi+sinphi*sinpsi;
r(2,1)=sinphi*costheta;
r(2,2)=sinphi*sintheta*sinpsi+cosphi*cospsi;
r(2,3)=sinphi*sintheta*cospsi-cosphi*sinpsi;
r(3,1)=-sintheta;
r(3,2)=costheta*sinpsi;
r(3,3)=costheta*cospsi;
return;

repo: ekfmonoslam


text:
function calpak_test013 ( )

%*****************************************************************************80
%
%% CALPAK_TEST013 tests JED_TO_YMDF_ARMENIAN and YMDF_TO_JED_ARMENIAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST013\n' );
  fprintf ( 1, '  For the Armenian calendar:\n' );
  fprintf ( 1, '  JED_TO_YMDF_ARMENIAN: JED -> YMDF.\n' );
  fprintf ( 1, '  YMDF_TO_JED_ARMENIAN: YMDF -> JED.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  JED (in)    YMDF               JED (out)\n' );
  fprintf ( 1, '\n' );

  jed_epoch = epoch_to_jed_armenian ( );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0.0 )
      break
    end

    if ( jed_epoch <= jed1 )

      [ y2, m2, d2, f2 ] = jed_to_ymdf_armenian ( jed1 );

      s2 = ymdf_to_s_numeric ( y2, m2, d2, f2 );

      jed3 = ymdf_to_jed_armenian ( y2, m2, d2, f2 );

      fprintf ( 1, '  %11.2f  %20s  %11.2f\n', jed1, s2, jed3 );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m