text:
function p = predictOneVsAll(all_theta, X)
%PREDICT Predict the label for a trained one-vs-all classifier. The labels
%are in the range 1..K, where K = size(all_theta, 1).
%  p = PREDICTONEVSALL(all_theta, X) will return a vector of predictions
%  for each example in the matrix X. Note that X contains the examples in
%  rows. all_theta is a matrix where the i-th row is a trained logistic
%  regression theta vector for the i-th class. You should set p to a vector
%  of values from 1..K (e.g., p = [1; 3; 1; 2] predicts classes 1, 3, 1, 2
%  for 4 examples)

m = size(X, 1);
num_labels = size(all_theta, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters (one-vs-all).
%               You should set p to a vector of predictions (from 1 to
%               num_labels).
%
% Hint: This code can be done all vectorized using the max function.
%       In particular, the max function can also return the index of the
%       max element, for more information see 'help max'. If your examples
%       are in rows, then, you can use max(A, [], 2) to obtain the max
%       for each row.
%

[max_value, p] = max(sigmoid(X * all_theta'), [], 2);

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

end

author: gopaczewski

repo: coursera-ml


text:
% Extract the coefficient of a monomial t in the polynomial p
%
% Syntax:   >> c = PolynomialCoefficient(p,t)
% Syntax:   >> c = pcoef(p,t)
%
%   Input:    p --- (string) polynomial
%             t --- (string) monomial
%
%  Output:    c --- (numeric or string) the coefficient of t in p
%
% Example:   >> p = '-2*x + 9*x*y^2 - 9*x^2*y^2 + 8*x^3*y^2';
%            >> c = Coeff(p,'y^2')
%
%            c =
%
%            9*x - 9*x^2 + 8*x^3
%
%            >> d = Coeff(p,'x*y^2')
%
%            d =
%
%            9

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function x = PQgetData (WAV, i, N)
% Get data from internal buffer or file
% i - file position
% N - number of samples
% x - output data (scaled to the range -32768 to +32767)

% Only two files can be "active" at a time.
% N = 0 resets the buffer

% P. Kabal $Revision: 1.1 $  $Date: 2003/12/07 13:34:10 $

persistent Buff

iB = WAV.iB + 1;
if (N == 0)
    Buff(iB).N = 20 * 1024;     % Fixed size
    Buff(iB).x = PQ_ReadWAV (WAV, i, Buff(iB).N);
    Buff(iB).i = i;
end

if (N > Buff(iB).N)
    error ('>>> PQgetData: Request exceeds buffer size');
end

% Check if requested data is not already in the buffer
is = i - Buff(iB).i;
if (is < 0 | is + N - 1 > Buff(iB).N - 1)
    Buff(iB).x = PQ_ReadWAV (WAV, i, Buff(iB).N);
    Buff(iB).i = i;
end

% Copy the data
Nchan = WAV.Nchan;
is = i - Buff(iB).i;
x = Buff(iB).x(1:Nchan,is+1:is+N-1+1);

%------
function x = PQ_ReadWAV (WAV, i, N)
% This function considers the data to extended with zeros before and
% after the data in the file. If the starting offset i is negative,
% zeros are filled in before the data starts at offset 0. If the request
% extends beyond the end of data in the file, zeros are appended.

Amax = 32768;
Nchan = WAV.Nchan;

x = zeros (Nchan, N);

Nz = 0;
if (i < 0)
    Nz = min (-i, N);
    i = i + Nz;
end

Ns = min (N - Nz, WAV.Nframe - i);
if (i >= 0 & Ns > 0)
    x(1:Nchan,Nz+1:Nz+Ns-1+1) = Amax * (wavread (WAV.Fname, [i+1 i+Ns-1+1]))';
end

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
function [ pn, dist ] = quad_point_near_2d ( q, p )

%*****************************************************************************80
%
%% QUAD_POINT_NEAR_2D computes the nearest point on a quadrilateral in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real Q(2,4), the quadrilateral vertices.
%
%    Input, real P(2), the point whose nearest quadrilateral point
%    is to be determined.
%
%    Output, real PN(2), the nearest point to P.
%
%    Output, real DIST, the distance from the point to the
%    quadrilateral.
%
  dim_num = 2;
  nside = 4;
%
%  Find the distance to each of the line segments that make up the edges
%  of the quadrilateral.
%
  dist = Inf;
  pn(1:dim_num) = 0.0;

  for j = 1 : nside

    jp1 = i4_wrap ( j+1, 1, nside );

    [ pn2, dist2, tval ] = segment_point_near_2d ( q(1:dim_num,j)', ...
      q(1:dim_num,jp1)', p );

    if ( dist2 < dist )
      dist = dist2;
      pn(1:dim_num) = pn2(1:dim_num);
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [index,numBeforeDeg]=monomial2DegNegLexOrderIdx(monomial,numBeforeDeg)
%%MONOMIAL2DEGNEGLEXORDERIDX Given the nonegative integer exponents of a
%          monomial term (e.g. x^5*y^3*z^6 is [5;3;6]), determine the index
%          of that monomial term in a list of all monomial terms in degree
%          negative lexicographic ordering. The degree of a monomial is the
%          sum of the exponents of the terms.
%
%INPUTS: monomial A numDimXnumMonomials set of numMonomials monomial terms
%                 consisting of nonnegative integer vectors of exponents of
%                 the elements Alternatively, if one just wants the
%                 numBeforeDegVec vector, then this should be an empty
%                 matrix and numBeforeDegVec should contain the
%                 dimensionality and the maximum degree desired.
%    numBeforeDeg An optional parameter. Providing this parameter will
%                 speed up repeated calls to this function by avoiding the
%                 recalculation of certain values. This vector is such that
%                 numBeforeDeg(i+1) is the number of monomial terms
%                 before degree i starting with degree 0, where
%                 numBeforeDeg(0+1)=0. This parameter is returned by the
%                 function. Note that numBeforeDeg depends on numDim, so
%                 vectors cannot be reused between cals where numDim
%                 changes. If one just wants this parameter (for example
%                 in advance to speed up many repeated calls), then set
%                 monomial=[] and numBeforeDeg=[numDim;maxDeg], where
%                 maxDeg the maximum degree desired, >=0.
%
%OUTPUTS: index A numMonomialsX1 vector of the position of the monomials in
%              the ordering of monomials in degree negative lexicographic
%              order, starting from 1. If monomial=[] and numBeforeDeg
%              is set to the dimensionality and maximum degree, then this
%              is an empty matrix.
% numBeforeDeg The same vector as the optional input. If the input is
%              not provided, then numBeforeDeg will be filled with
%              values up to the degree of the monomial provided having the
%              highest degree. If monomials=[] and numBeforeDeg on the
%              input is set to the maximum degree desired, then
%              numBeforeDeg will be filled all the way up to that
%              degree.
%
%Various monomial orderings are used in algorithms that manipulate
%polynomials, such as methods for solving or minimizing multivariate
%polynomials as in [1] and [2].
%
%The degree negative lexicographic ordering, for a fixed degree, coincides
%with the ordering of the terms provided by unrankTComposition with
%firstElMostSig=true. However, that rank is only an offset for a given
%degree. The number of terms before that degree is the sum of the number of
%compositions of lower degrees. Thus, one sees the need for the
%numBeforeDeg vector. The number of compositions of a given degree into
%n parts, (each part starting from 0, not 1) is
%binomial(degree+numDim-1,numDim-1).
%
%EXAMPLE:
%Consider monomials like x^1*y^0*z^0, x^12*y^4*z^8;, x*y*z, x^0*y^0*z^0(=1)
% monomial=[1,12,1,0;
%           0, 4,1,0;
%           0, 8,1,0];
% [index,numBeforeDeg]=monomial2DegNegLexOrderIdx(monomial)
% %One will get the indices index=[2;2687;15;1].
% %The value numBeforeDeg can be reused with monomials having the same
% %number of variables, such as
% monomial=[0;1;1];
% index=monomial2DegNegLexOrderIdx(monomial,numBeforeDeg)
% %On the other hand, if one only wants the numBeforeDeg vector for up to
% %degree 24 for three variables, then one
% [~,numBeforeDeg]=monomial2DegNegLexOrderIdx([],[3,24])
%
%REFERENCES:
%[1] P. Dreesen, "Back to the roots: Polynomial system solving using linear
%    algebra," Ph.D. dissertation, Katholieke Universiteit Leuven, Leuven,
%    Flanders, Belgium, Sep. 2013.
%[2] A. Papachristodoulou and S. Prajna, "A tutorial on sum of squares
%    techniques for systems analysis," in Proceedings of the American
%    Control Conference, Portland, OR, 8-10 Jun. 2005, pp. 2686-2700.
%
%February 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(isempty(monomial))
    numDim=numBeforeDeg(1);
    maxDeg=numBeforeDeg(2);
else
    numDim=size(monomial,1);
    degs=sum(monomial,1);
    maxDeg=max(degs);
end

%Fill in the numBeforeDegVec vector.
if(isempty(monomial)||nargin<2||isempty(numBeforeDeg))
    numBeforeDeg=zeros(maxDeg+1,1);

    if(maxDeg>0)
        numBeforeDeg(1+1)=1;
    end

    for curDeg=2:maxDeg
        numBeforeDeg(curDeg+1)=numBeforeDeg(curDeg-1+1)+binomial(curDeg-1+numDim-1,numDim-1);
    end

    %If just the list numBeforeDegVec is desired.
    if(isempty(monomial))
        index=[];
        return;
    end
end

numMonomials=size(monomial,2);
%Allocate space.
index=zeros(1,numMonomials);

for curMonomial=1:numMonomials
    offset=numBeforeDeg(degs(curMonomial)+1);
    index(curMonomial)=rankTComposition(monomial(:,curMonomial)+1,true)+offset+1;
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 [L,W] = get_l(n,d)
%GET_L Compute discrete derivative operators.
%
% [L,W] = get_l(n,d)
%
% Computes the discrete approximation L to the derivative operator
% of order d on a regular grid with n points, i.e. L is (n-d)-by-n.
%
% L is stored as a sparse matrix.
%
% Also computes W, an orthonormal basis for the null space of L.

% Per Christian Hansen, IMM, 02/05/98.

% Initialization.
if (d<0), error ('Order d must be nonnegative'), end

% Zero'th derivative.
if (d==0), L = speye(n); W = zeros(n,0); return, end

% Compute L.
c = [-1,1,zeros(1,d-1)];
nd = n-d;
for i=2:d, c = [0,c(1:d)] - [c(1:d),0]; end
L = sparse(nd,n);
for i=1:d+1
  L = L + sparse(1:nd,(1:nd)+i-1,c(i)*ones(1,nd),nd,n);
end

% If required, compute the null vectors W via modified Gram-Schmidt.
if (nargout==2)
  W = zeros(n,d);
  W(:,1) = ones(n,1);
  for i=2:d, W(:,i) = W(:,i-1).*(1:n)'; end
  for k=1:d
     W(:,k) = W(:,k)/norm(W(:,k));
     W(:,k+1:d) = W(:,k+1:d) - W(:,k)*(W(:,k)'*W(:,k+1:d));
  end
end

repo: REST


text:
function ierror = perm0_check ( n, p )

%*****************************************************************************80
%
%% PERM0_CHECK checks a permutation of ( 0, ..., N-1).
%
%  Discussion:
%
%    The routine verifies that each of the integers from 0 to
%    to N-1 occurs among the N entries of the permutation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 May 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries.
%
%    Input, integer P(N), the array to check.
%
%    Output, integer IERROR:
%    0, P is a legal permutation of (0,...,N-1).
%    1, P is not a legal permuation of ( 0,...,N-1).
%
  ierror = 0;

  for value = 0 : n - 1

    ierror = 1;

    for location = 1 : n
      if ( p(location) == value )
        ierror = 0;
        break;
      end
    end

    if ( ierror ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PERM0_CHECK - Warning!\n' );
      fprintf ( 1, '  Permutation is missing the value %d.\n', value );
      break;
    end

  end

  return
end


text:
%UT_WEIGHTS - Generate unscented transformation weights
%
% Syntax:
%   [WM,WC,c] = ut_weights(n,alpha,beta,kappa)
%
% In:
%   n     - Dimensionality of random variable
%   alpha - Transformation parameter  (optional, default 0.5)
%   beta  - Transformation parameter  (optional, default 2)
%   kappa - Transformation parameter  (optional, default 3-n)
%
% Out:
%   WM - Weights for mean calculation
%   WC - Weights for covariance calculation
%    c - Scaling constant
%
% Description:
%   Computes unscented transformation weights.
%
% See also UT_MWEIGHTS UT_TRANSFORM UT_SIGMAS
%

% Copyright (C) 2006 Simo S�rkk�
%
% $Id$
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.

function [WM,WC,c] = ut_weights(n,alpha,beta,kappa)

%
% Check which arguments are there
%
if nargin < 1
  error('At least dimensionality n required.');
end
if nargin < 2
  alpha = [];
end
if nargin < 3
  beta = [];
end
if nargin < 4
  kappa = [];
end

%
% Apply default values
%
if isempty(alpha)
  alpha = 1;
end
if isempty(beta)
  beta = 0;
end
if isempty(kappa)
  kappa = 3 - n;
end

%
% Compute the normal weights
%
lambda = alpha^2 * (n + kappa) - n;

WM = zeros(2*n+1,1);
WC = zeros(2*n+1,1);
for j=1:2*n+1
  if j==1
    wm = lambda / (n + lambda);
    wc = lambda / (n + lambda) + (1 - alpha^2 + beta);
  else
    wm = 1 / (2 * (n + lambda));
    wc = wm;
  end
  WM(j) = wm;
  WC(j) = wc;
end

c = n + lambda;

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
filename = 'CantileverSquareSym';%'CantileverSquareNewFine';%'LshapeTriFine';%'CantileverSquareNew';'ArchTriFine';%'CantileverSquare';%'ArchTriFine';%'CantileverSquare';%'Lshape';%'LshapeTriFine';%'CantileverSquareSmall';%'';%'Lshape';'CantileverSquare';'Lshape';%'CantileverSquare';%'LshapeFine';%'Bridge_quad_coarse';%'Arch_quad_coarse';%'BridgeCool_Quadrilateral_Bilinear_Structured_Coarse';%'Bridge';%'CantileverSquareSmall';
ptype = 'MACRO';
initial_case = 'given';
% m1 = 0.0101;
% m2 = 0.0101;
m1 = 0.8;
m2 = 0.8;
%cost = {'compliance'};
cost = {'stressNorm'};
%cost = {'stressNorm','compliance'};
%weights = [0.55,0.45];
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'P1';
constraint_case = 'EQUALITY';

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

Vfrac_final = 0.3;
optimality_final = 1e-5;
constr_final = 1e-5;

stressNormExponent_initial = 2;
stressNormExponent_final = 32;

optimizer = 'DualNestedInPrimal';
optimizerUnconstrained = 'PROJECTED GRADIENT';

designVariable = 'MicroParams';
ub = 0.989;
lb = 0.011;
homegenizedVariablesComputer = 'ByVademecum';
%
vademecumFileName = 'SuperEllipseQMax';
%vademecumFileName = 'SuperEllipseQ2';
%vademecumFileName = 'SuperEllipseQOptAnalytic';
%
% designVariable = 'Density';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';

kfrac = 2;
nsteps = 50;

plotting = false;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 500;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function p_xt= Motion_Model_Velocity_Test(u, x_p, a, gridsize)
%% Test to figure out which 'a' parameters suit you best
    % u: control command
    % x_p: initial position with orientation
    % a: [a_1 ... a_6 ] of tested parameters

    %the test will calculate the probability of the robot moving to every
    %possible cell in the grid generating a banana-shaped density function
    p_xt=zeros(gridsize(1), gridsize(2));
    for i=1:gridsize(1)
        for j=1:gridsize(2)-1
            p_xt(i,j) = Motion_Model_Velocity([j i x_p(3)], u, x_p, a) ;
        end
    end

    % Visualize results
    figure(66)
    clf;
    hold on;
    plot(x_p(1), x_p(2), 'dr');
    surf(p_xt);
    view([0 -90]);
    hold off;

end

author: Allopart

repo: rbpf-gmapping


text:
% VL_TEST_AIBHIST
function vl_test_aibhist

D   = 4 ;
K   = 20 ;

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

X1 = randn(2,300) ; X1(1,:) = X1(1,:) + 2 ;
X2 = randn(2,300) ; X2(1,:) = X2(1,:) - 2 ;
X3 = randn(2,300) ; X3(2,:) = X3(2,:) + 2 ;

C   = 1:K*K ;
Pcx = zeros(3,K*K) ;

f1 = quantize(X1,D,K) ;
f2 = quantize(X2,D,K) ;
f3 = quantize(X3,D,K) ;

Pcx(1,:) = vl_binsum(Pcx(1,:), ones(size(f1)), f1) ;
Pcx(2,:) = vl_binsum(Pcx(2,:), ones(size(f2)), f2) ;
Pcx(3,:) = vl_binsum(Pcx(3,:), ones(size(f3)), f3) ;

Pcx = Pcx / sum(Pcx(:)) ;

[parents_, cost_] = vl_aib(Pcx) ;
[parents,  cost ] = vl_aib(Pcx,'clusternull') ;

% find a null node for testing purposes
anull = min(find(parents_==0)) ;
f1 = [f1 repmat(anull,1,10)] ;

figure(100);  clf ;
subplot(1,2,1) ; hold on ;
plot(parents_, 'r.-') ;
plot(parents,  'g') ;
legend('signal null', 'cluster null') ;
subplot(1,2,2) ; hold on ;
plot(cost_, 'r.-') ;
plot(cost,  'g') ;
legend('signal null', 'cluster null') ;

range = [1 10 K*K-10 K*K] ;
for c=1:length(range)

  cut_size = range(c) ;

  % compare two methods of getting the same cut histogram
  [cut_,map_]  = vl_aibcut(parents_, cut_size) ;
  hist_        = vl_aibcuthist(map_, f1, 'nulls', 'append') ;
  histtree_    = vl_aibhist(parents_, f1) ;
  thist_       = histtree_(cut_) ;

  [cut,map]    = vl_aibcut(parents, cut_size) ;
  hist         = vl_aibcuthist(map, f1, 'nulls', 'append') ;
  histtree     = vl_aibhist(parents, f1) ;
  thist        = histtree(cut) ;

  figure(100 + c) ; clf ;
  subplot(2,2,1) ; hold on ; plot(hist_,'g.-') ; plot(thist_,'r') ;
  legend('cut+cuthist', 'hist+cut') ;
  title('vl_aibcuthist vs aibhist');
  subplot(2,2,2) ; hold on ; plot(histtree_) ; title('aibtree') ;

  subplot(2,2,3) ; hold on ; plot(hist,'g.-') ; plot(thist,'r') ;
  legend('cut+cuthist', 'hist+cut') ;
  title('vl_aibcuthist vs vl_aibhist (clust null)');
  subplot(2,2,4) ; hold on ; plot(histtree) ; title('aibtree (clust null)') ;
end

% --------------------------------------------------------------------
function f = quantize(X,D,K)
% --------------------------------------------------------------------
d = 2*D / K ;
j = round((X(1,:) + D) / d) ;
i = round((X(2,:) + D) / d) ;
j = max(min(j,K),1) ;
i = max(min(i,K),1) ;
f = sub2ind([K K],i,j) ;

author: jianxiongxiao

repo: ProfXkit


text:
function Write_ContinuousDynamics(FileWritingSetup)
% Write_ContinuousDynamics(FileWritingSetup)
%
% FUNCTION:
%
%   This function takes the symbolically derived equations of motion for
%   the retractable double pendulum biped model and writes them to function
%   files.
%
%
% INPUTS: FileWritingSetup.
%
%   States = [2N x 2] cell array where the first column is the name of
%       the states, and the second column is comments about that states.
%       The index of a state in this array will match the index in the
%       dynamics function to be written. This function assumes a second
%       order system that is being written in first order form, so if
%       States{i,1} = z, then States{i+N,1} = dz; Furthermore, the output
%       (dStates) will have the form: dStates(i,:) = dz,
%       and dStates(i+N,:) = ddz.
%
%   Actuators = [P x 2] cell array giving the names of each of the
%       actuators and their comments.
%
%   Contacts = [M x 2] cell array giving all of the contact forces to
%       output. The first column is names, second column is comments
%
%   Parameters = [K x 2] cell array with all of the input parameters. The
%       first column is for names, the second for comments
%
%   FileData = [L x 3+?] cell array with the names of each field in Dyn in
%       the first column, and the file name to write those dynamics
%       equations to. The final column is used for comments that are
%       specific to those dynamics. Seperate multiple lines by commas:
%       'comment line one', 'comment line two', ...
%
%   Dyn = A struct, with a field for each dynamics equation to write. Every
%       field must contain a symbolic expression for each of the state
%       accelerations and contact forces.
%
%   Directory = the name of the directory to write the files in.
%
%
% OUTPUTS:
%   This function will write a file for the dynamics equations in each
%   field of Dyn.
%
%
% NOTE:
%
%   1)  It is assumed that the symbols in all of the inputs are self
%       consistent; that is - if Contacts{1,1} = H1, then there must be a field
%       in Dyn.(FileNames{i,1}) called H1.
%
%
% Written by Matthew Kelly
% December 6, 2013
% Cornell University
%
% See also DERIVE_EOM

States = FileWritingSetup.States;
Dyn = FileWritingSetup.Dyn;
Contacts = FileWritingSetup.Contacts;
Actuators = FileWritingSetup.Actuators;
Parameters = FileWritingSetup.Parameters;
FileData = FileWritingSetup.FileData;
Directory = FileWritingSetup.Directory;

CurrDir = cd;
if ~exist(Directory,'dir')
    mkdir(Directory);
end
cd(Directory);

for idxPhase = 1:size(FileData,1)
    fid = fopen([FileData{idxPhase,2} '.m'],'w');

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
    %                         Function Header                             %
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

    fprintf(fid,['function [dStates, contactForces] = ' ...
        FileData{idxPhase,2} '(States, Actuators, Parameters)\n']);
    fprintf(fid,['%% function [dStates, Actuators] = ' ...
        FileData{idxPhase,2} '(States, Actuators, Parameters)\n']);
    fprintf(fid,'%%\n%% Computer Generated File -- DO NOT EDIT \n');
    fprintf(fid,['%%\n%% This function was created by the function '...
        'Write_ContinuousDynamics()\n']);
    fprintf(fid,['%% ' datestr(now) '\n%%\n']);
    for idxComment = 3:length(FileData(idxPhase,:))
        fprintf(fid,['%% ' FileData{idxPhase,idxComment} '\n']);
    end %idxComment
    fprintf(fid,'%%\n');
    fprintf(fid,'%% Matthew Kelly \n%% Cornell University \n%% \n\n');
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
    %                           Parse Inputs                              %
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

    for idxState=1:size(States,1)
        fprintf(fid,[States{idxState,1} ' = States(' num2str(idxState) ',:); %% ' ...
            States{idxState,2} '\n']);
    end %idxState
    fprintf(fid,'\n');

    for idxAct=1:size(Actuators,1)
        rhs = [' = Actuators(' num2str(idxAct) ',:); %% ' Actuators{idxAct,2} '\n'];
        %Need to make sure that we disconnect the ankle motors when their
        %respective feed are not in contact with the ground!
        if strcmp(FileData{idxPhase,1},'Flight') || strcmp(FileData{idxPhase,1},'SingleTwo')
            if strcmp(Actuators{idxAct,1},'T1')
                rhs = ' = 0; %%Foot One not in contact with ground!\n';
            end
        end
        if strcmp(FileData{idxPhase,1},'Flight') || strcmp(FileData{idxPhase,1},'SingleOne')
            if strcmp(Actuators{idxAct,1},'T2')
                rhs = ' = 0; %%Foot Two not in contact with ground!\n';
            end
        end
        fprintf(fid,[Actuators{idxAct,1} rhs]);
    end %idxAct
    fprintf(fid,'\n');

    for idxParam=1:size(Parameters,1)
        fprintf(fid,[Parameters{idxParam,1} ' = Parameters.' ...
            Parameters{idxParam,1} '; %% ' ...
            Parameters{idxParam,2} '\n']);
    end %idxState
    fprintf(fid,'\n');

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
    %                   Explicitly write constraints                      %
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

    Ndof = size(States,1)/2;
    fprintf(fid,'%% Constraints for this phase: \n');
    for i=1:Ndof
        writeIndex = i + Ndof;
        writeName = ['d' States{writeIndex,1}];
        writeData = Dyn.(FileData{idxPhase,1}).(writeName);
        if (writeData==0)
            fprintf(fid,[writeName ' = 0; %%' States{i,2} '\n']);
        end
    end
    for i=1:size(Contacts,1)
        writeName = Contacts{i,1};
        writeData = Dyn.(FileData{idxPhase,1}).(writeName);
        if (writeData==0)
            fprintf(fid,[writeName ' = 0; %%' Contacts{i,2} '\n']);
        end
    end
    fprintf(fid,'\n');

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
    %                       Write Accelerations                           %
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

     N = num2str(Ndof);
    fprintf(fid,'dStates = zeros(size(States));\n');
    fprintf(fid,['dStates(1:' N ',:) = '...
        'States((1+' N '):(' N '+' N '),:);\n']);
    for idxAcc=1:Ndof
        writeIndex = idxAcc + Ndof;
        writeName = ['d' States{writeIndex,1}];
        writeData = Dyn.(FileData{idxPhase,1}).(writeName);
        vectorizedExpression = vectorize(writeData);
        fprintf(fid,['dStates(' num2str(writeIndex) ',:) = ' ...
            vectorizedExpression ';\n']);
    end %idxAcc
    fprintf(fid,'\n');

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
    %                       Write Contact Forces                          %
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

    for idxContact=1:size(Contacts,1)
        fprintf(fid,['%% contactForces(' num2str(idxContact) ',:) == '...
            Contacts{idxContact,1} ' == ' Contacts{idxContact,2} '\n']);
    end %idxContact
    fprintf(fid,['contactForces = zeros(' num2str(size(Contacts,1))...
        ',size(States,2));\n']);
    for idxContact=1:size(Contacts,1)
        writeIndex = num2str(idxContact);
        writeName = Contacts{idxContact,1};
        writeData = Dyn.(FileData{idxPhase,1}).(writeName);
        vectorizedExpression = vectorize(writeData);
        fprintf(fid,['contactForces(' writeIndex ',:) = ' ...
            vectorizedExpression ';\n']);
    end %idxContact
    fprintf(fid,'\n');

    fprintf(fid,'end\n');
    fclose(fid);
end %idxPhase
cd(CurrDir);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  ~~~~~~~~~~~~~~~~~~~ Sub Functions  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function characterExpression = vectorize(symbolicExpression)

characterExpression = char(symbolicExpression);
characterExpression = strrep(characterExpression,'*','.*');
characterExpression = strrep(characterExpression,'/','./');
characterExpression = strrep(characterExpression,'^','.^');
characterExpression = strrep(characterExpression,...
    '0','zeros(1,size(States,2))');

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function varargout = opensignals_txt(filename, hdr, begsample, endsample, chanindx)

% OPENSIGNALS_TXT reads time series data from a Bitalino OpenSignals txt file
%
% Use as
%   hdr = opensignals_txt(filename);
%   dat = opensignals_txt(filename, hdr, begsample, endsample, chanindx);
%   evt = opensignals_txt(filename, hdr);
%
% See also FT_FILETYPE, FT_READ_HEADER, FT_READ_DATA, FT_READ_EVENT, QUALISYS_TSV, MOTION_C3D

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

needhdr = (nargin==1);
needevt = (nargin==2);
needdat = (nargin==5);

% use the full filename including path to distinguish between similarly named files in different directories
[p, f, x] = fileparts(filename);
if isempty(p)
  % no path was specified
  fullname = which(filename);
elseif startsWith(p, ['.' filesep])
  % a relative path was specified
  fullname = fullfile(pwd, p(3:end), [f, x]);
else
  fullname = filename;
end

fid = fopen_or_error(fullname, 'r');
line1 = fgetl(fid);
line2 = fgetl(fid);
fclose(fid);
% convert the json string into a structure
orig = jsondecode(line2(2:end));
fn = fieldnames(orig);
if numel(fn)==1
  % it contains a single field with the date and time
  orig = orig.(fn{1});
end

dat = readtable(fullname, 'FileType', 'text');
keep = false(1,size(dat,2));
for i=1:size(dat,2)
  keep(i) = isnumeric(dat{:,i});
end
dat = table2array(dat(:,keep))';

if needhdr
  %% parse the header
  hdr.Fs = orig.samplingRate;
  hdr.nChans = length(orig.column);
  hdr.nTrials = 1;
  hdr.nSamplesPre = 0;
  hdr.nSamples = size(dat,2);
  hdr.label = orig.column;

  % remember the details
  hdr.orig = orig;

  % return the header details
  varargout = {hdr};

elseif needdat
  %% parse the data
  dat = dat(chanindx, begsample:endsample);

  % return the data
  varargout = {dat};

elseif needevt
  %% parse the events
  ft_warning('cannot read events from %s', filename);

  varargout = {[]};
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%Generates state evolution and compares to observed MSEs for D-AMP and D-IT algorithms

addpath(genpath('..'));

%Parameters
denoiser='BM3D';
SamplingRate=.2;
noise_sig=0;
imsize=128;
filename='barbara.png';
iters=30;
N=5;%Number of tests to run.  Must be at least 2

Im=double(imread(filename));
x_0=imresize(Im,imsize/size(Im,1));
[height, width]=size(x_0);
n=length(x_0(:));
m=round(n*SamplingRate);

%Generate State Evolution and compute intermediate MSEs of D-AMP and D-IT
True_DAMP_MSE_array=zeros(N,iters);
True_DIT_MSE_array=zeros(N,iters);
Predicted_MSE_array=zeros(N,iters);
True_DAMP_MSE_array(:,1)=mean(x_0(:).^2);
True_DIT_MSE_array(:,1)=mean(x_0(:).^2);
Predicted_MSE_array(:,1)=mean(x_0(:).^2);
for i=1:N
    i
    M=randn(m,n);
    for j = 1:n
        M(:,j) = M(:,j) ./ sqrt(sum(abs(M(:,j)).^2));
    end
    y=M*x_0(:)+noise_sig*randn(m,1);
    x_t=zeros(n,1);x_t2=zeros(n,1);
    z_t=y;
    z_t2=y;
    for iter=2:iters
            [x_tplus1,z_tplus1,NA] = DAMP_oneIter(y,z_t,x_t,width,height,denoiser,M);
            z_t=z_tplus1; x_t=x_tplus1;
            [x_tplus12,z_tplus12,NA] = DIT_oneIter(y,z_t2,x_t2,width,height,denoiser,M);
            z_t2=z_tplus12; x_t2=x_tplus12;
            True_DAMP_MSE_array(i,iter)=mean((x_0(:)-x_tplus1(:)).^2);
            True_DIT_MSE_array(i,iter)=mean((x_0(:)-x_tplus12(:)).^2);
            Predicted_MSE_array(i,iter)= DAMP_SE_Prediction(x_0(:), Predicted_MSE_array(i,iter-1), m,n,noise_sig,denoiser,width,height);
    end
end
True_DAMP_MSE=mean(True_DAMP_MSE_array)';
True_DAMP_MSE_std=std(True_DAMP_MSE_array)';
True_DIT_MSE=mean(True_DIT_MSE_array)';
True_DIT_MSE_std=std(True_DIT_MSE_array)';
Predicted_MSE=mean(Predicted_MSE_array)';
Predicted_MSE_std=std(Predicted_MSE_array)';

%Plot Results
h=figure;
hold
errorbar(0:29,Predicted_MSE,Predicted_MSE_std,'-.b');
errorbar(0:29,True_DAMP_MSE,True_DAMP_MSE_std,'--g');
errorbar(0:29,True_DIT_MSE,True_DIT_MSE_std,'-r');
title([denoiser,'-AMP and ', denoiser '-IT State Evolution']);
legend(['Predicted ',num2str(100*SamplingRate),'%'], [denoiser,'-AMP ', num2str(100*SamplingRate),'%'], [denoiser,'-IT ', num2str(100*SamplingRate),'%']);
xlabel('Iteration');
ylabel('MSE');

author: ricedsp


text:
classdef StressVoigt3DRotator < StressVoigtRotator

    methods (Access = public)

        function obj = StressVoigt3DRotator(angle,dir)
            obj.compute(angle,dir)
        end
    end

    methods (Access = protected)

        function generateRotator(obj)
            theta = obj.angle;
            d  = obj.dir.getValue();
            u1 = d(1);
            u2 = d(2);
            u3 = d(3);

            A(1,1) = ((1 - cos(theta))*u1^2 + cos(theta))^2;
            A(2,2) = ((1 - cos(theta))*u2^2 + cos(theta))^2;
            A(3,3) = ((1 - cos(theta))*u3^2 + cos(theta))^2;

            A(4,4) = ((1 - cos(theta))*u2^2 + cos(theta))*((1 - cos(theta))*u3^2 + cos(theta)) - (u1*sin(theta) + u2*u3*(cos(theta) - 1))*(u1*sin(theta) - u2*u3*(cos(theta) - 1));
            A(5,5) = ((1 - cos(theta))*u1^2 + cos(theta))*((1 - cos(theta))*u3^2 + cos(theta)) - (u2*sin(theta) + u1*u3*(cos(theta) - 1))*(u2*sin(theta) - u1*u3*(cos(theta) - 1));
            A(6,6) = ((1 - cos(theta))*u1^2 + cos(theta))*((1 - cos(theta))*u2^2 + cos(theta)) - (u3*sin(theta) + u1*u2*(cos(theta) - 1))*(u3*sin(theta) - u1*u2*(cos(theta) - 1));

            A(1,2) = (u3*sin(theta) - u1*u2*(cos(theta) - 1))^2;
            A(2,1) = (u3*sin(theta) + u1*u2*(cos(theta) - 1))^2;

            A(1,3) = (u2*sin(theta) + u1*u3*(cos(theta) - 1))^2;
            A(3,1) = (u2*sin(theta) - u1*u3*(cos(theta) - 1))^2;

            A(1,4) = -2*(u2*sin(theta) + u1*u3*(cos(theta) - 1))*(u3*sin(theta) - u1*u2*(cos(theta) - 1));
            A(4,1) = -(u3*sin(theta) + u1*u2*(cos(theta) - 1))*(u2*sin(theta) - u1*u3*(cos(theta) - 1));

            A(1,5) = -2*(u2*sin(theta) + u1*u3*(cos(theta) - 1))*((1 - cos(theta))*u1^2 + cos(theta));
            A(5,1) = (u2*sin(theta) - u1*u3*(cos(theta) - 1))*((1 - cos(theta))*u1^2 + cos(theta));

            A(1,6) = 2*(u3*sin(theta) - u1*u2*(cos(theta) - 1))*((1 - cos(theta))*u1^2 + cos(theta));
            A(6,1) = -(u3*sin(theta) + u1*u2*(cos(theta) - 1))*((1 - cos(theta))*u1^2 + cos(theta));

            A(2,3) = (u1*sin(theta) - u2*u3*(cos(theta) - 1))^2;
            A(3,2) = (u1*sin(theta) + u2*u3*(cos(theta) - 1))^2;

            A(2,4) = 2*(u1*sin(theta) - u2*u3*(cos(theta) - 1))*((1 - cos(theta))*u2^2 + cos(theta));
            A(4,2) = -(u1*sin(theta) + u2*u3*(cos(theta) - 1))*((1 - cos(theta))*u2^2 + cos(theta));

            A(2,5) = -2*(u3*sin(theta) + u1*u2*(cos(theta) - 1))*(u1*sin(theta) - u2*u3*(cos(theta) - 1));
            A(5,2) = -(u1*sin(theta) + u2*u3*(cos(theta) - 1))*(u3*sin(theta) - u1*u2*(cos(theta) - 1));

            A(2,6) = -2*(u3*sin(theta) + u1*u2*(cos(theta) - 1))*((1 - cos(theta))*u2^2 + cos(theta));...
                A(6,2) = (u3*sin(theta) - u1*u2*(cos(theta) - 1))*((1 - cos(theta))*u2^2 + cos(theta));

            A(3,4) = -2*(u1*sin(theta) + u2*u3*(cos(theta) - 1))*((1 - cos(theta))*u3^2 + cos(theta));
            A(4,3) = (u1*sin(theta) - u2*u3*(cos(theta) - 1))*((1 - cos(theta))*u3^2 + cos(theta));

            A(3,5) = 2*(u2*sin(theta) - u1*u3*(cos(theta) - 1))*((1 - cos(theta))*u3^2 + cos(theta));
            A(5,3) = -(u2*sin(theta) + u1*u3*(cos(theta) - 1))*((1 - cos(theta))*u3^2 + cos(theta));

            A(3,6) = -2*(u1*sin(theta) + u2*u3*(cos(theta) - 1))*(u2*sin(theta) - u1*u3*(cos(theta) - 1));
            A(6,3) = -(u2*sin(theta) + u1*u3*(cos(theta) - 1))*(u1*sin(theta) - u2*u3*(cos(theta) - 1));

            A(4,5) = (u1*sin(theta) - u2*u3*(cos(theta) - 1))*(u2*sin(theta) - u1*u3*(cos(theta) - 1)) - (u3*sin(theta) + u1*u2*(cos(theta) - 1))*((1 - cos(theta))*u3^2 + cos(theta));
            A(5,4) = (u1*sin(theta) + u2*u3*(cos(theta) - 1))*(u2*sin(theta) + u1*u3*(cos(theta) - 1)) + (u3*sin(theta) - u1*u2*(cos(theta) - 1))*((1 - cos(theta))*u3^2 + cos(theta));

            A(4,6) = (u1*sin(theta) + u2*u3*(cos(theta) - 1))*(u3*sin(theta) + u1*u2*(cos(theta) - 1)) + (u2*sin(theta) - u1*u3*(cos(theta) - 1))*((1 - cos(theta))*u2^2 + cos(theta));
            A(6,4) = (u1*sin(theta) - u2*u3*(cos(theta) - 1))*(u3*sin(theta) - u1*u2*(cos(theta) - 1)) - (u2*sin(theta) + u1*u3*(cos(theta) - 1))*((1 - cos(theta))*u2^2 + cos(theta));

            A(5,6) = (u2*sin(theta) - u1*u3*(cos(theta) - 1))*(u3*sin(theta) - u1*u2*(cos(theta) - 1)) - (u1*sin(theta) + u2*u3*(cos(theta) - 1))*((1 - cos(theta))*u1^2 + cos(theta));
            A(6,5) = (u2*sin(theta) + u1*u3*(cos(theta) - 1))*(u3*sin(theta) + u1*u2*(cos(theta) - 1)) + (u1*sin(theta) - u2*u3*(cos(theta) - 1))*((1 - cos(theta))*u1^2 + cos(theta));

            obj.rotationMatrix = A;

        end
    end

end


text:
function [ a, isplit ] = r8vec_split ( n, a, split )

%*****************************************************************************80
%
%% R8VEC_SPLIT "splits" an unsorted R8VEC based on a splitting value.
%
%  Discussion:
%
%    If the vector is already sorted, it is simpler to do a binary search
%    on the data than to call this routine.
%
%    The vector is not assumed to be sorted before input, and is not
%    sorted during processing.  If sorting is not needed, then it is
%    more efficient to use this routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of elements of A.
%
%    Input, real A(N), the array to split.
%
%    Input, real SPLIT, the value used to split the vector.
%    It is not necessary that any value of A actually equal SPLIT.
%
%    Output, real A(N), all the entries of A that are less than or
%    equal to SPLIT are in A(1:ISPLIT).
%
%    Output, integer ISPLIT, indicates the position of the last
%    entry of the split vector that is less than or equal to SPLIT.
%

%
%  Partition the vector into A1, A2, A3, where
%    A1 = A(I1:J1) holds values <= SPLIT,
%    A2 = A(I2:J2) holds untested values,
%    A3 = A(I3:J3) holds values > SPLIT.
%
  i1 = 1;
  j1 = 0;

  i2 = 1;
  j2 = n;

  i3 = n + 1;
  j3 = n;
%
%  Pick the next item from A2, and move it into A1 or A3.
%  Adjust indices appropriately.
%
  for i = 1 : n

    if ( a(i2) <= split )

      i2 = i2 + 1;
      j1 = j1 + 1;

    else

      temp = a(i2);
      a(i2) = a(i3-1);
      a(i3-1) = temp;

      i3 = i3 - 1;
      j2 = j2 - 1;

    end

  end

  isplit = j1;

  return
end


text:
function sparse_grid_mixed_weight_tests ( )

%*****************************************************************************80
%
%% SPARSE_GRID_MIXED_WEIGHT_TESTS calls SPARSE_GRID_MIXED_WEIGHT_TEST with various arguments.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2011
%
%  Author:
%
%    John Burkardt
%
%  Local Parameters:
%
%    Local, real TOL, a tolerance for point equality.
%    A value of sqrt ( eps ) is reasonable, and will allow the code to
%    consolidate points which are equal, or very nearly so.  A value of
%    -1.0, on the other hand, will force the code to use every point, regardless
%    of duplication.
%
  addpath ( '../sandia_rules' );

  tol = sqrt ( eps );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_MIXED_WEIGHT_TESTS\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Call SPARSE_GRID_MIXED_WEIGHT with various arguments.\n' );
  fprintf ( 1, '  All tests will use a point equality tolerance of %e\n', tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 1 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 3 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 4 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 7 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 8 ]';
  alpha = [ 0.0, 1.5 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 2, 9 ]';
  alpha = [ 0.0, 0.5 ]';
  beta = [ 0.0, 1.5 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 6, 4 ]';
  alpha = [ 2.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 3;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 2, 5 ]';
  alpha = [ 0.0, 0.0, 0.0 ]';
  beta = [ 0.0, 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );
%
%  Dimension 2, Level 4, Rule 3.
%
  dim_num = 2;
  level_max_min = 0;
  level_max_max = 4;
  rule = [ 3, 3 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );
%
%  Dimension 2, Level 4, Rule 13.
%
  dim_num = 2;
  level_max_min = 0;
  level_max_max = 4;
  rule = [ 13, 13 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );
%
%  Dimension 2, Level 4, Rule 16.
%
  dim_num = 2;
  level_max_min = 0;
  level_max_max = 4;
  rule = [ 16, 16 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );
%
%  Dimension 2, Level 4, Rule 17.
%
  dim_num = 2;
  level_max_min = 0;
  level_max_max = 4;
  rule = [ 17, 17 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_MIXED_WEIGHT_TESTS\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../sandia_rules' );

  return
end


text:
function [Hx,Hxtest] = computeHx()
    Hx = eye(3,6);
    Hxtest = eye(6);
end

author: ccalas

repo: mpc

save_path: github-repos/MATLAB/ccalas-mpc


text:
function [net, data] = gtminit2(net, options, data, init, samp_type, varargin)
%GTMINIT Initialise the weights and latent sample in a GTM given reference vectors.
%
%	Description
%	NET = GTMINIT(NET, OPTIONS, DATA, INIT, SAMPTYPE) takes a GTM NET and
%	generates a sample of latent data points and sets the centres (and
%	widths if appropriate) of NET.RBFNET.
%
% INIT defines initial locations of reference vectors (Gaussian means).
%
%	If the SAMPTYPE is 'REGULAR', then regular grids of latent data
%	points and RBF centres are created.  The dimension of the latent data
%	space must be 1 or 2.  For one-dimensional latent space, the
%	LSAMPSIZE parameter gives the number of latent points and the
%	RBFSAMPSIZE parameter gives the number of RBF centres.  For a two-
%	dimensional latent space, these parameters must be vectors of length
%	2 with the number of points in each of the x and y directions to
%	create a rectangular grid.  The widths of the RBF basis functions are
%	set by a call to RBFSETFW passing OPTIONS(7) as the scaling
%	parameter.
%
%	If the SAMPTYPE is 'UNIFORM' or 'GAUSSIAN' then the latent data is
%	found by sampling from a uniform or Gaussian distribution
%	correspondingly.  The RBF basis function parameters are set by a call
%	to RBFSETBF with the DATA parameter as dataset and the OPTIONS
%	vector.
%
%	Finally, the output layer weights of the RBF are initialised by
%	mapping the mean of the latent variable to the mean of the target
%	variable, and the L-dimensional latent variale variance to the
%	variance of the targets along the first L principal components.
%
%	See also
%	GTMINIT, GTM, GTMEM, PCA, RBFSETBF, RBFSETFW
%

%	Copyright (c) Ian T Nabney (1996-2001)
% Modified by Tommi Vatanen (2012)

% Check for consistency
errstring = consist(net, 'gtm', data);
if ~isempty(errstring)
  error(errstring);
end

% Check type of sample
stypes = {'regular', 'uniform', 'gaussian', 'special'};
if (strcmp(samp_type, stypes)) == 0
  error('Undefined sample type.')
end

if net.dim_latent > size(data, 2)
  error('Latent space dimension must not be greater than data dimension')
end
nlatent = net.gmmnet.ncentres;
nhidden = net.rbfnet.nhidden;

% Create latent data sample and set RBF centres

switch samp_type
case 'regular'
   if nargin ~= 7
      error('Regular type must specify latent and RBF shapes');
   end
   l_samp_size = varargin{1};
   rbf_samp_size = varargin{2};
   if round(l_samp_size) ~= l_samp_size
      error('Latent sample specification must contain integers')
   end
   % Check existence and size of rbf specification
   if any(size(rbf_samp_size) ~= [1 net.dim_latent]) || ...
         prod(rbf_samp_size) ~= nhidden
      error('Incorrect specification of RBF centres')
   end
   % Check dimension and type of latent data specification
   if any(size(l_samp_size) ~= [1 net.dim_latent]) || ...
         prod(l_samp_size) ~= nlatent
      error('Incorrect dimension of latent sample spec.')
   end
   if net.dim_latent == 1
      net.X = [-1:2/(l_samp_size-1):1]';
      net.rbfnet.c = [-1:2/(rbf_samp_size-1):1]';
      net.rbfnet = rbfsetfw(net.rbfnet, options(7));
   elseif net.dim_latent == 2
      net.X = gtm_rctg(l_samp_size);
      net.rbfnet.c = gtm_rctg(rbf_samp_size);
      net.rbfnet = rbfsetfw(net.rbfnet, options(7));
   else
      error('For regular sample, input dimension must be 1 or 2.')
   end

case {'uniform', 'gaussian'}
   if strcmp(samp_type, 'uniform')
      net.X = 2 * (rand(nlatent, net.dim_latent) - 0.5);
   else
      % Sample from N(0, 0.25) distribution to ensure most latent
      % data is inside square
      net.X = randn(nlatent, net.dim_latent)/2;
   end
   net.rbfnet = rbfsetbf(net.rbfnet, options, net.X);

%   case {'special'}
%     net.X = [-0.5 sqrt(3)/2;0.5 sqrt(3)/2;0 0];
%     net.rbfnet.c = [-0.5 sqrt(3)/2;0.5 sqrt(3)/2;0 0];
%     net.rbfnet = rbfsetfw(net.rbfnet, options(7));
otherwise
   % Shouldn't get here
   error('Invalid sample type');
end

% use given unit vectors as the initial reference vectors
[~, Phi] = rbffwd(net.rbfnet, net.X);
% initparam are the reference vectors
net.rbfnet.w2 = Phi \ init;
net.rbfnet.b2 = zeros(1,size(data,2));

net.gmmnet.centres = rbffwd(net.rbfnet, net.X);
d = dist2(net.gmmnet.centres, net.gmmnet.centres) + ...
  diag(ones(net.gmmnet.ncentres, 1)*realmax);
sigma = mean(min(d))/2;

net.gmmnet.covars = sigma*ones(1, net.gmmnet.ncentres);

% Sub-function to create the sample data in 2d
function sample = gtm_rctg(samp_size)

xDim = samp_size(1);
yDim = samp_size(2);
% Produce a grid with the right number of rows and columns
[X, Y] = meshgrid([0:1:(xDim-1)], [(yDim-1):-1:0]);

% Change grid representation
sample = [X(:), Y(:)];

% Shift grid to correct position and scale it
maxXY= max(sample);
sample(:,1) = 2*(sample(:,1) - maxXY(1)/2)./maxXY(1);
sample(:,2) = 2*(sample(:,2) - maxXY(2)/2)./maxXY(2);

return;

author: ilarinieminen

repo: SOM-Toolbox


text:
classdef EnsembleClassifier < Algorithm

    properties (Access = public)
        shouldTrain = false;
        nLearners = 30;
    end

    properties (Access = private)
        classifier;
    end

    methods (Access = public)

        function obj = EnsembleClassifier()
            obj.name = 'Ensemble';
            obj.inputPort = DataType.kTable;
            obj.outputPort = DataType.kTable;
        end

        function dataOut = compute(obj,data)
            if obj.shouldTrain
                obj.train(data);
                dataOut = [];
            else
                dataOut = obj.test(data);
            end
        end

        function train(obj,table)
            template = templateTree(...
                'MaxNumSplits', 1852);

            if ~isempty(template)
                obj.classifier = fitcensemble(...
                    table.features, ...
                    table.label, ...
                    'Method', 'Bag', ...
                    'NumLearningCycles', obj.nLearners, ...
                    'Learners', template);
            end
        end

        function metrics = computeMetrics(obj,table)
            flops = timeit(@()obj.test(table)) / Constants.kReferenceComputingTime;
            memory = Helper.ComputeObjectSize(obj.classifier);
            outputSize = table.height * Constants.kClassificationResultBytes;
            metrics = Metric(flops,memory,outputSize);
        end

        function labels = test(obj,table)
            labels = predict(obj.classifier,table.features);
        end

        function str = toString(obj)
            str = sprintf('%s_%d_%d',obj.name,obj.nLearners);
        end

        function editableProperties = getEditableProperties(obj)
            editableProperties = Property('nLearners',obj.nLearners,1,100,PropertyType.kNumber);
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function q=dcm2q(R)
% DCM2Q(R) converts direction cosine matrices into quaternions.
%
%     The resultant quaternion(s) will perform the equivalent vector
%     transformation as the input DCM(s), i.e.:
%
%       qconj(Q)*V*Q = R*V
%
%     where R is the DCM, V is a vector, and Q is the quaternion.  Note that
%     for purposes of quaternion-vector multiplication, a vector is treated
%     as a quaterion with a scalar element of zero.
%
%     If the input is a 3x3xN array, the output will be a vector of
%     quaternions where input direction cosine matrix R(:,:,k) corresponds
%     to the output quaternion Q(k,:).
%
% See also Q2DCM.

% Release: $Name: quaternions-1_2_2 $
% $Revision: 1.9 $
% $Date: 2002/01/21 06:46:20 $

% Copyright (C) 2000-2002, Jay A. St. Pierre.  All rights reserved.

if nargin~=1
  error('dcm2q() requires one input argument');
else
  size_R=size(R);
  if ( size_R(1)~=3 | size_R(2)~=3 | length(size_R)>3 )
    error(['Invalid input: must be a 3x3xN array'])
  end
end

q(1,:)=0.5*sqrt(1 + R(1,1,:) - R(2,2,:) - R(3,3,:)).*sgn(R(2,3,:)-R(3,2,:));
q(2,:)=0.5*sqrt(1 - R(1,1,:) + R(2,2,:) - R(3,3,:)).*sgn(R(3,1,:)-R(1,3,:));
q(3,:)=0.5*sqrt(1 - R(1,1,:) - R(2,2,:) + R(3,3,:)).*sgn(R(1,2,:)-R(2,1,:));
q(4,:)=0.5*sqrt(1 + R(1,1,:) + R(2,2,:) + R(3,3,:));

% Make quaternion vector a column of quaternions
q=q.';

q=real(q);

% Signum function (for our purposes)
function s=sgn(x)
  s=sign(x);
  s=s+(s==0);

author: christianwengert


text:
function [A_hat, E_hat, Y, mu, iter] = inexact_alm_rpca(D, lambda, tol, maxIter)

% Oct 2009
% This matlab code implements the inexact augmented Lagrange multiplier
% method for Robust PCA.
%
% D - m x n matrix of observations/data (required input)
%
% lambda - weight on sparse error term in the cost function
%
% tol - tolerance for stopping criterion.
%     - DEFAULT 1e-7 if omitted or -1.
%
% maxIter - maximum number of iterations
%         - DEFAULT 1000, if omitted or -1.
%
% Initialize A,E,Y,u
% while ~converged
%   minimize (inexactly, update A and E only once)
%     L(A,E,Y,u) = |A|_* + lambda * |E|_1 + <Y,D-A-E> + mu/2 * |D-A-E|_F^2;
%   Y = Y + \mu * (D - A - E);
%   \mu = \rho * \mu;
% end
%
% Minming Chen, October 2009. Questions? v-minmch@microsoft.com ;
% Arvind Ganesh (abalasu2@illinois.edu)
%
% Copyright: Perception and Decision Laboratory, University of Illinois, Urbana-Champaign
%            Microsoft Research Asia, Beijing

[m,n] = size(D);

if(nargin < 2) lambda = 1 / sqrt(m); end
if(nargin < 3) tol = 1e-7; elseif(tol == -1) tol = 1e-7; end
if(nargin < 4) maxIter = 1000; elseif(maxIter == -1) maxIter = 1000; end

% initialize
Y = D;
norm_two = norm(Y, 2);
norm_inf = norm( Y(:), inf) / lambda;
dual_norm = max(norm_two, norm_inf);
Y = Y / dual_norm;

A_hat = zeros( m, n);
E_hat = zeros( m, n);
mu = 1.25/norm_two; % this one can be tuned
mu_bar = mu * 1e7;
rho = 1.5;          % this one can be tuned
d_norm = norm(D, 'fro');

iter = 0;
total_svd = 0;
converged = false;
stopCriterion = 1;
sv = 10;
while ~converged
    iter = iter + 1;

    temp_T = D - A_hat + (1/mu)*Y;
    E_hat = max(temp_T - lambda/mu, 0);
    E_hat = E_hat+min(temp_T + lambda/mu, 0);

    %[U,S,V] = svd(D - E_hat + (1/mu)*Y, 'econ');
    [U,S,V] = svdecon(D - E_hat + (1/mu)*Y); % fastest

    diagS = diag(S);
    svp = length(find(diagS > 1/mu));
    if svp < sv
        sv = min(svp + 1, n);
    else
        sv = min(svp + round(0.05*n), n);
    end

    A_hat = U(:, 1:svp) * diag(diagS(1:svp) - 1/mu) * V(:, 1:svp)';

    total_svd = total_svd + 1;

    Z = D - A_hat - E_hat;

    Y = Y + mu*Z;
    mu = min(mu*rho, mu_bar);

    %% stop Criterion
    stopCriterion = norm(Z, 'fro') / d_norm;
    if stopCriterion < tol
        converged = true;
    end

    if mod( total_svd, 10) == 0
        disp(['#svd ' num2str(total_svd) ' r(A) ' num2str(rank(A_hat))...
            ' |E|_0 ' num2str(length(find(abs(E_hat)>0)))...
            ' stopCriterion ' num2str(stopCriterion)]);
    end

    if ~converged && iter >= maxIter
        disp('Maximum iterations reached') ;
        converged = 1 ;
    end
end

author: andrewssobral

repo: lrslibrary


text:
function state = pbsGetJobState(scheduler, job, state)
%pbsGetJobState Gets the state of a job on a cluster.
%
% Set your schedulers's GetJobStateFcn to this function using the following
% command (see README):
%     set(sched, 'GetJobStateFcn', @pbsGetJobState);

%  Copyright 2006-2008 The MathWorks, Inc.

mlock;
persistent jobsToMonitorNames;

if isempty(jobsToMonitorNames)
    jobsToMonitorNames = {};
end

if isempty(scheduler.UserData)
    if ~iscell(scheduler.SubmitFcn) || length(scheduler.SubmitFcn) < 3
        error('distcomp:genericscheduler:SubmitFcnError',...
              'SubmitFcn must include clusterHost and remoteDataLocation as extra arguments.');
    end
    scheduler.UserData={ scheduler.SubmitFcn{2} ; scheduler.SubmitFcn{3} };
end

jobName = job.pGetEntityLocation;
if strcmp(state, 'finished')
% emc 6 March 2010 replace setxor with setdiff - behavior with empty set
      jobsToMonitorNames = setdiff(jobsToMonitorNames, {jobName});
else
%   timerName as suggested by Ray Norris TMW 5 March 2010
%   timerName = [job.Name '_' num2str(job.ID)];
%   timerName emc 6 March 2010
    timerName = [job.Name '_' strrep(job.SubmitTime, ' ','_')];
    if strcmp(state, 'queued') || strcmp(state, 'running')
        if ~ismember(jobName, jobsToMonitorNames) || ...
                                    isempty(timerfind('Name',timerName))
            jobsToMonitorNames = union(jobsToMonitorNames, {jobName});
            jobStateTimer = ...
                      timer('Period', 10.0, 'ExecutionMode', 'fixedRate');
            set(jobStateTimer, 'TimerFcn', ...
                         { @copyJobFilesIfFinished, scheduler, job }, ...
                                                        'Name',timerName);
            start(jobStateTimer);
        end
    end
end


text:
% Load data "bias60.mat".
clc;
clear;
load bias60;

% Create a figure.
figure('DefaultAxesFontSize',16);
ax = gca;

% Create a heatmap.
mat = bias{2}.*ones(61,144)+bias{3}'.*ones(61,144);
imagesc(mat);
c = colorbar;
colormap hot;
ylabel(c,'Speed (km/h)');
ax.XTick = [1,13,25,37,49,61,73,85,97,109,121,133,144];
ax.XTickLabel = ['00:00';'02:00';'04:00';'06:00';'08:00';'10:00';'12:00';...
    '14:00';'16:00';'18:00';'20:00';'22:00';'23:50'];
xtickangle(45);
ax.YTick = [1,6,11,16,21,26,31,36,41,46,51,56,61];
ax.YTickLabel = ['08-01';'08-06';'08-11';'08-16';'08-21';'08-26';'08-31';...
    '09-05';'09-10';'09-15';'09-20';'09-25';'09-30'];

% Save the figure as "heat_map60.pdf".
set(gcf, 'PaperSize', [6 4.5]);
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperPosition', [0 0 6 4.5]);
saveas(gcf,'heat_map60','pdf');

author: xinychen

repo: academic-drawing


text:
init_fig;

% c
axes('position', [.05, .57, .95, .37]);
hold on;
plot(y, 'color', col{8}/255);
alpha(.7);
plot(true_c, 'color', col{3}/255, 'linewidth', 1.5);
plot(c_oasis, '-.', 'color', col{1}/255);
axis tight;
xlim([0, 2000]);
set(gca, 'xtick', [0, 25, 50, 75]*30);
set(gca, 'xticklabel', []);
set(gca, 'ytick', 0:2);
ylabel('Fluor.');
box off;
legend('Data', 'Truth', 'OASIS', 'location', 'northeast', 'orientation', 'horizental');
% s
axes('position', [.05, .18, .95, .37]);
hold on;
plot(true_s, 'color', col{3}/255, 'linewidth', 1.5);
plot(s_oasis, '-.', 'color', col{1}/255);
axis tight;
xlim([0, 2000]);
set(gca, 'xtick', [0, 25, 50, 75]*30);
set(gca, 'xticklabel', get(gca, 'xtick')/30);
set(gca, 'ytick', [0,1]);
xlabel('Time [s]');
ylabel('Activity.');

author: zhoupc


text:
% This file load eq catalogs in hypoellipes format. Since the formated reading
% is fast, this is the best way to load data for
% large catalogs.
% Earthquakes within a rectengular box (latmin latmax lonmin lonmax)
% a time window and a magnitude threshold can be selected.
% Note: lat/Lon are assumed to be in minutes in the original catalog
%
% last modified: March 96   stefan wiemer

report_this_filefun(mfilename('fullpath'));

%  close the file in case it is still open
safe_fclose(fid);

% selection criteria
tmin = 60.0001;
tmax = 97.00900;
lonmin = -180.0;
lonmax =  180.0;
latmin = -90.0;
latmax = 90.0 ;
Mmin = 3.1
deperr = 1000.5;

a = []; b = [];
n2 = 0;

% open the file and read 2000 lines at a time
[file1,path1] = uigetfile([ '*.mat'],' Earthquake Datafile');
fid = fopen([path1 file1],'r') ;;

while  ferror(fid) == ''
    l =fscanf(fid,'%2d%2d%2d%2d%2d%4d%2d%*1c%4d%3d%*1c%4d%5f%2d%*31c%2d%*15c%4d',[14 2000]) ;
    n2 = n2 +  length(l(1,:));

    b = [ -l(9,:)-l(10,:)/6000 ; l(7,:)+l(8,:)/6000 ;...
        l(1,:);l(2,:);l(3,:);l(13,:)/10;l(11,:)/100;l(4,:);...
        l(5,:); l(14,:)/100;l(12,:)/10];
    b = b';
    l =  b(:,6) >= Mmin & b(:,1) >= lonmin & b(:,1) <= lonmax & ...
        b(:,2) >= latmin & b(:,2) <= latmax & b(:,3) <= tmax  & ...
        b(:,3) >= tmin & b(:,10) <= deperr;
    a = [a ; b(l,:)];
    disp([ num2str(n2) ' earthquakes scanned, ' num2str(length(a)) ' EQ found'])
    if max(b(:,3)) >  tmax ; break; end
end
ferror(fid)
fclose(fid);

dep1 = 0.3*(max(a.Depth)-min(a.Depth))+min(a.Depth);
dep2 = 0.6*(max(a.Depth)-min(a.Depth))+min(a.Depth);
dep3 = max(a.Depth);

stri1 = [file1];
tim1 = minti;
tim2 = maxti;
minma2 = minma;
maxma2 = maxma;
minde = min(a.Depth);
maxde = max(a.Depth);
rad = 50.;
ic = 0;
ya0 = 0.;
xa0 = 0.;
iwl3 = 1.;
step = 3;
t1p(1) = 80.;
t2p(1) = 85.;
t3p(1) = 90.;
t4p(1) = 93.;
tresh = 10;

% save the catalog
save newcat.mat a

author: CelsoReyes


text:
% VL_DEMO_QUICKSHIFT  Demo: Quick shift: basic functionality

pfx = fullfile(vl_root,'figures','demo') ;
randn('state',0) ;
rand('state',0) ;
figure(1) ; clf ;

% --------------------------------------------------------------------
%                                                        Load a figure
% --------------------------------------------------------------------
I = imread(fullfile(vl_root,'data','roofs1.jpg')) ;

image(I) ;
axis equal off tight ;
vl_demo_print('quickshift_image') ;

% --------------------------------------------------------------------
%                                     Create a quickshift segmentation
% --------------------------------------------------------------------

ratio = 0.5;
kernelsize = 2;
maxdist = 10;

Iseg = vl_quickseg(I, ratio, kernelsize, maxdist);
image(Iseg);
axis equal off tight;
vl_demo_print('quickshift_qseg_1') ;

kernelsize = 2;
maxdist = 20;
Iseg = vl_quickseg(I, ratio, kernelsize, maxdist);
image(Iseg);
axis equal off tight;
vl_demo_print('quickshift_qseg_2') ;

maxdist = 50;
ndists = 10;
Iedge = vl_quickvis(I, ratio, kernelsize, maxdist, ndists);
imagesc(Iedge);
axis equal off tight;
colormap gray;
vl_demo_print('quickshift_qvis') ;

author: Cloud-CV

repo: object-proposals


text:
Require Import HoareDef IntroHeader IntroF0 IntroF1 SimModSem.
Require Import Coqlib.
Require Import ImpPrelude.
Require Import Skeleton.
Require Import PCM.
Require Import ModSem Behavior.
Require Import Relation_Definitions.

(*** TODO: export these in Coqlib or Universe ***)
Require Import Relation_Operators.
Require Import RelationPairs.
From ITree Require Import
     Events.MapDefault.
From ExtLib Require Import
     Core.RelDec
     Structures.Maps
     Data.Map.FMapAList.

Require Import HTactics.

Set Implicit Arguments.

Local Open Scope nat_scope.

Section SIMMODSEM.

  Context `{Σ: GRA.t}.

  Let W: Type := Any.t * Any.t.

  Let wf: unit -> W -> Prop :=
    fun _ '(mrps_src0, mrps_tgt0) =>
      (<<SRC: mrps_src0 = tt↑>>) /\
      (<<TGT: mrps_tgt0 = tt↑>>)
  .

  Theorem correct:
    refines2 [IntroF0.F] [IntroF1.F].
  Proof.
    eapply adequacy_local2. econs; ss. i.
    econstructor 1 with (wf:=wf) (le:=top2); et; ss.
    econs; ss. init. unfold cfunU.
    unfold fF, IntroF0.fF.
    steps. des. clarify. ss. steps.
    force_r. { eapply max_intrange; ss. } steps.
    destruct (Z.of_nat x <? 0)%Z eqn:T; try lia.
    unfold Ncall.
    destruct (Z.of_nat x =? 0)%Z eqn:U.
    { steps. assert(x = 0) by lia. steps. force_l. exists false. steps. force_l. eexists (Vint 0). steps.
      force_l; ss. steps. force_l. eexists (Vint _). steps. force_l; ss. steps. r. esplits; et.
    }
    steps. destruct x0.
    - force_l. exists true. unfold ccallU. steps.
      rewrite Z.eqb_neq in *. force_l; ss. { splits; try lia. eapply OrdArith.lt_from_nat. lia. } steps.
      force_r. { f_equal. lia. } steps.
      force_l. eexists. steps. force_l; ss. steps. r. esplits; et. do 2 f_equal. lia.
    - steps. force_l. exists false. steps. force_l. eexists. steps. force_l; ss. steps.
      force_l. esplits. steps. force_l; ss. steps. r. esplits; et. do 2 f_equal. lia.
  Unshelve. all: try exact 0. all: ss.
  Qed.

End SIMMODSEM.

author: alxest

repo: CCR

save_path: github-repos/coq/alxest-CCR


text:
Require Import sflib.

Require Import Axioms.
Require Import Basic.
Require Import DataStructure.
Require Import Loc.
Require Import Language.

Require Import Time.
Require Import Event.
Require Import View.
Require Import Cell.
Require Import Memory.
Require Import TView.
Require Import Local.
Require Import Thread.

Set Implicit Arguments.

(** * The Program Step in the Promising Semantics 2.1 *)

(** This file defines the program steps (machine steps) in promising semantics 2.1. *)
(** The program steps in PS2.1 defined in this file corresponds to
    the definitions in Figure 9 in our paper. *)

Module Threads.

  (** ** Thread Pool  *)
  (** It is a set of threads and corresponds to (ThrdPool) in paper. *)
  Definition t := IdentMap.t ({lang:language & (Language.state lang)} * Local.t).

  (** Add new threads in the thread pool (for initialization). *)
  Definition add_th {lang: language}
             (ths:t) (tid:IdentMap.key) (f: Language.fid) (code: Language.syntax lang) : option t :=
    match ((Language.init lang) code f) with
    | Some s => Some (IdentMap.add tid (existT _ _ s, Local.init) ths)
    | _ => None
    end.

  Definition tid_incr := POrderedType.Positive_as_DT.succ.
  Fixpoint add_ths {lang: language}
           (ths: t) (tid: IdentMap.key) (fs: list Language.fid) (code: Language.syntax lang) : option t :=
    match fs with
    | f::fs' => match add_ths ths (tid_incr tid) fs' code with
               | Some ths' => add_th ths' tid f code
               | _ => None
               end
    | _ => Some ths
    end.

  Definition empty_ths := IdentMap.empty ({lang:language & (Language.state lang)} * Local.t).

  Definition init {lang: language} (fs: list Language.fid) (code: Language.syntax lang) :=
    add_ths empty_ths BinNums.xH fs code.

  Definition is_terminal (ths:t): Prop :=
    forall tid lang st lc (FIND: IdentMap.find tid ths = Some (existT _ lang st, lc)),
      <<STATE: (Language.is_terminal lang) st>> /\
      <<THREAD: Local.is_terminal lc>>.

  Inductive wf (ths:t) (mem:Memory.t): Prop :=
  | wf_intro
      (DISJOINT:
         forall tid1 lang1 st1 lc1
           tid2 lang2 st2 lc2
           (TID: tid1 <> tid2)
           (TH1: IdentMap.find tid1 ths = Some (existT _ lang1 st1, lc1))
           (TH2: IdentMap.find tid2 ths = Some (existT _ lang2 st2, lc2)),
           Local.disjoint lc1 lc2)
      (THREADS: forall tid lang st lc
                  (TH: IdentMap.find tid ths = Some (existT _ lang st, lc)),
          Local.wf lc mem)
  .

  Definition consistent (ths:t) (sc:TimeMap.t) (mem:Memory.t) (lo: Ordering.LocOrdMap): Prop :=
    forall tid lang st lc
      (TH: IdentMap.find tid ths = Some (existT _ lang st, lc)),
      Thread.consistent (Thread.mk lang st lc sc mem) lo.

  Definition consistent_nprm (ths:t) (sc:TimeMap.t) (mem:Memory.t) (lo: Ordering.LocOrdMap): Prop :=
    forall tid lang st lc
      (TH: IdentMap.find tid ths = Some (existT _ lang st, lc)),
      Thread.consistent_nprm (Thread.mk lang st lc sc mem) lo.

  Inductive is_promised tid (loc:Loc.t) (to:Time.t) (threads:t): Prop :=
  | is_promised_intro
      lang st lc from msg
      (TID: IdentMap.find tid threads = Some (existT _ lang st, lc))
      (PROMISES: Memory.get loc to (Local.promises lc) = Some (from, msg))
  .

  (* tids *)

  Definition tids (ths: t): IdentSet.t :=
    List.fold_right (fun p s => IdentSet.add (fst p) s) IdentSet.empty (IdentMap.elements ths).
End Threads.

Module Configuration.
  (** ** Program Configuration *)
  (** The program configuration in promising semantics 2.1. It includes:
      - [threads]: thread pool;
      - [tid]: the identifier of the current thread;
      - [sc]: the timemap for SC fence;
      - [memory]: memory shared by all threads;

      It corresponds to (World) in paper. *)
  Structure t := mk {
    threads: Threads.t;
    tid: IdentMap.key;
    sc: TimeMap.t;
    memory: Memory.t;
  }.

  (** ** Initialization of the Program in the Promising Semantics 2.1 *)
  Definition init {lang: language} (fs: list Language.fid) (code: Language.syntax lang) (ctid: IdentMap.key) :=
    match (Threads.init fs code) with
    | Some ths => Some (mk ths ctid TimeMap.bot Memory.init)
    | _ => None
    end.

  Definition is_terminal (conf:t): Prop := Threads.is_terminal (threads conf).

  Inductive wf (conf:t): Prop :=
  | wf_intro
      (WF: Threads.wf (threads conf) (memory conf))
      (SC: Memory.closed_timemap (sc conf) (memory conf))
      (MEM: Memory.closed (memory conf))
  .

  Definition sw_tid (c1: t) (new_tid: IdentMap.key): t := mk (threads c1) new_tid (sc c1) (memory c1).
  Lemma tid_eq_dec: forall (tid1 tid2: IdentMap.key),
      {tid1 = tid2} + {tid1 <> tid2}.
  Proof. eapply IdentSet.Raw.L.MO.eq_dec; eauto. Qed.

  Lemma sw_to_self
        c:
    c = sw_tid c (tid c).
  Proof.
    destruct c; ss.
  Qed.

  Lemma sw_tid_twice
        c tid1 tid2:
    Configuration.sw_tid (Configuration.sw_tid c tid1) tid2 =
    Configuration.sw_tid c tid2.
  Proof.
    destruct c; ss.
  Qed.

  Definition consistent (lo: Ordering.LocOrdMap) (c: t): Prop :=
    Threads.consistent (threads c) (sc c) (memory c) lo.

  Definition no_promise (c: t): Prop :=
    forall tid lang st lc (FIND: IdentMap.find tid (threads c) = Some (existT _ lang st, lc)),
      (Local.promises lc) = Memory.bot.

  (** ** Program Step *)
  (** Program step (or machine step) in promising semantics 2.1.
      It takes four parameters:
      - [e]: machine event;
      - [lo]: a set of atomic locations;
      - [c1] and [c2]: program configurations before and after transition. *)
  Inductive step: forall (e:MachineEvent.t) (lo: Ordering.LocOrdMap) (c1 c2:t), Prop :=
  | step_tau
      lang lo c1 c2 st1 st2 lc1 lc2 e1 tid1
      (CTID: tid c1 = tid1)
      (TID1: IdentMap.find tid1 (threads c1) = Some (existT _ lang st1, lc1))
      (STEPS: rtc (@Thread.tau_step _ lo) (Thread.mk _ st1 lc1 (sc c1) (memory c1)) e1)
      (STEP: Thread.tau_step lo e1 (Thread.mk _ st2 lc2 (sc c2) (memory c2)))
      (** The execution of current thread must reach a promise certified state. *)
      (CONSISTENT: Thread.consistent (Thread.mk _ st2 lc2 (sc c2) (memory c2)) lo)
      (C2: c2 = mk (IdentMap.add tid1 (existT _ _ st2, lc2) (threads c1)) tid1 (sc c2) (memory c2)):
      step MachineEvent.silent lo c1 c2
  | step_out
      lang lo c1 c2 st1 st2 lc1 lc2 e1 e tid1
      (CTID: tid c1 = tid1)
      (TID1: IdentMap.find tid1 (threads c1) = Some (existT _ lang st1, lc1))
      (STEPS: rtc (@Thread.tau_step _ lo) (Thread.mk _ st1 lc1 (sc c1) (memory c1)) e1)
      (STEP_OUT: Thread.out_step lo e e1 (Thread.mk _ st2 lc2 (sc c2) (memory c2)))
      (CONSISTENT: Thread.consistent (Thread.mk _ st2 lc2 (sc c2) (memory c2)) lo)
      (C2: c2 = mk (IdentMap.add tid1 (existT _ _ st2, lc2) (threads c1)) tid1 (sc c2) (memory c2)):
      step (MachineEvent.syscall e) lo c1 c2
  | step_sw
      tid' lang c1 c2 st2 lc2 lo
      (TID2_OK: IdentMap.find tid' (threads c1) = Some (existT _ lang st2, lc2))
      (C2: c2 = Configuration.mk (threads c1) tid' (sc c1) (memory c1)):
      step MachineEvent.switch lo c1 c2
  | step_thread_term
      lang c1 c2 st1 lc1 st2 lc2 lo tid1 tid'
      (CTID: tid c1 = tid1)
      (OLD_TID: IdentMap.find tid1 (threads c1) = Some (existT _ lang st1, lc1))
      (THRD_DONE: Thread.is_done (Thread.mk _ st1 lc1 (sc c1) (memory c1)))
      (THRDS_REMOVE: ((IdentMap.remove tid1 (threads c1)) = (threads c2)))
      (NEW_TID_OK: IdentMap.find tid' (threads c2) = Some (existT _ lang st2, lc2))
      (C2: c2 = mk (threads c2) tid' (sc c1) (memory c1)):
      step MachineEvent.switch lo c1 c2
  .
  Hint Constructors step.

  (** ** Done Configuration *)
  (** All of threads terminate. *)
  Definition is_done (c: t): Prop :=
    exists lang st lc,
      (IdentMap.find (tid c) (threads c) = Some (existT _ lang st, lc)) /\
      Thread.is_done (Thread.mk _ st lc (sc c) (memory c)) /\
      IdentMap.is_empty (IdentMap.remove (tid c) (threads c)) = true.

  Definition is_abort (c: t) (lo: Ordering.LocOrdMap): Prop :=
      exists lang st1 lc1  e',
        IdentMap.find (tid c) (threads c) = Some (existT _ lang st1, lc1) /\
        rtc (@Thread.tau_step _ lo) (Thread.mk _ st1 lc1 (sc c) (memory c)) e' /\
        Thread.is_abort e' lo.

  Inductive tau_step: Ordering.LocOrdMap -> t -> t -> Prop :=
  | tau_step_intro
      e lo c1 c2
      (WP_STEP: step e lo c1 c2)
      (TAU: ~(exists e0, e = MachineEvent.syscall e0)):
      tau_step lo c1 c2.

  Inductive all_step: Ordering.LocOrdMap -> t -> t -> Prop :=
  | all_step_intro
      e lo c1 c2
      (WP_STEP: step e lo c1 c2):
      all_step lo c1 c2.

  (** ** Abort Configuration *)
  Inductive abort_config: forall (lo: Ordering.LocOrdMap) (c: t), Prop :=
  | Abort_config_intro
      lo c c'
      (STEPS: rtc (all_step lo) c c')
      (ABORT: is_abort c' lo):
      abort_config lo c.

  (** ** Safe Program *)
  (** A program is safe, if its executions will never reach an abort configuration. *)
  Inductive safe {lang: language} :
    forall (lo: Ordering.LocOrdMap) (fs: list Language.fid) (code: Language.syntax lang) (ctid: IdentMap.key), Prop :=
  | Safe_intro
      lo fs code ctid
      (SAFR_EXEC: forall c, init fs code ctid = Some c -> ~(abort_config lo c)):
      safe lo fs code ctid.

  (** Aux promising semantics (machine step).
    1. same init/term/sw/done/abort rule, done/abort definition
    2. we have na/prc/atm step now (substitute for single tau step)
    *)
  Inductive aux_step: forall (e:AuxEvent.t) (lo: Ordering.LocOrdMap) (c1 c2:t), Prop :=
  | aux_step_na
    lang lo c1 c2 st1 st2 lc1 lc2 e1 tid1
    (CTID: tid c1 = tid1)
    (TID1: IdentMap.find tid1 (threads c1) = Some (existT _ lang st1, lc1))
    (STEPS: rtc (@Thread.na_step _ lo) (Thread.mk _ st1 lc1 (sc c1) (memory c1)) e1)
    (STEP: Thread.na_step lo e1 (Thread.mk _ st2 lc2 (sc c2) (memory c2)))
    (CONSISTENT: Thread.consistent (Thread.mk _ st2 lc2 (sc c2) (memory c2)) lo)
    (C2: c2 = mk (IdentMap.add tid1 (existT _ _ st2, lc2) (threads c1)) tid1 (sc c2) (memory c2)):
    aux_step AuxEvent.na lo c1 c2
  | aux_step_prc
    lang lo c1 c2 st1 st2 lc1 lc2 e1 tid1
    (CTID: tid c1 = tid1)
    (TID1: IdentMap.find tid1 (threads c1) = Some (existT _ lang st1, lc1))
    (STEPS: rtc (@Thread.prc_step _ lo) (Thread.mk _ st1 lc1 (sc c1) (memory c1)) e1)
    (STEP: Thread.prc_step lo e1 (Thread.mk _ st2 lc2 (sc c2) (memory c2)))
    (CONSISTENT: Thread.consistent (Thread.mk _ st2 lc2 (sc c2) (memory c2)) lo)
    (C2: c2 = mk (IdentMap.add tid1 (existT _ _ st2, lc2) (threads c1)) tid1 (sc c2) (memory c2)):
    aux_step AuxEvent.prc lo c1 c2
  | aux_step_at
    lang lo c1 c2 st1 st2 lc1 lc2 e1 e1' tid1
    (CTID: tid c1 = tid1)
    (TID1: IdentMap.find tid1 (threads c1) = Some (existT _ lang st1, lc1))
    (STEPS_PRC: rtc (@Thread.prc_step _ lo) (Thread.mk _ st1 lc1 (sc c1) (memory c1)) e1)
    (STEPS_ATM: rtc (@Thread.atmblk_step _ lo) e1 e1')
    (STEP_ATM: Thread.atmblk_step lo e1' (Thread.mk _ st2 lc2 (sc c2) (memory c2)))
    (CONSISTENT: Thread.consistent (Thread.mk _ st2 lc2 (sc c2) (memory c2)) lo)
    (C2: c2 = mk (IdentMap.add tid1 (existT _ _ st2, lc2) (threads c1)) tid1 (sc c2) (memory c2)):
    aux_step AuxEvent.atm lo c1 c2
  | aux_step_out
    lang lo c1 c2 st1 st2 lc1 lc2 e1 e tid1
    (CTID: tid c1 = tid1)
    (TID1: IdentMap.find tid1 (threads c1) = Some (existT _ lang st1, lc1))
    (STEPS: rtc (@Thread.na_step _ lo) (Thread.mk _ st1 lc1 (sc c1) (memory c1)) e1)
    (STEP_OUT: Thread.out_step lo e e1 (Thread.mk _ st2 lc2 (sc c2) (memory c2)))
    (CONSISTENT: Thread.consistent (Thread.mk _ st2 lc2 (sc c2) (memory c2)) lo)
    (C2: c2 = mk (IdentMap.add tid1 (existT _ _ st2, lc2) (threads c1)) tid1 (sc c2) (memory c2)):
    aux_step (AuxEvent.out e) lo c1 c2
  | aux_step_sw
      tid' lang c1 c2 st2 lc2 lo
      (TID2_OK: IdentMap.find tid' (threads c1) = Some (existT _ lang st2, lc2))
      (C2: c2 = Configuration.mk (threads c1) tid' (sc c1) (memory c1)):
      aux_step AuxEvent.sw lo c1 c2
  | aux_step_thread_term
      lang c1 c2 st1 lc1 st2 lc2 lo tid1 tid'
      (CTID: tid c1 = tid1)
      (OLD_TID: IdentMap.find tid1 (threads c1) = Some (existT _ lang st1, lc1))
      (THRD_DONE: Thread.is_done (Thread.mk _ st1 lc1 (sc c1) (memory c1)))
      (THRDS_REMOVE: ((IdentMap.remove tid1 (threads c1)) = (threads c2)))
      (NEW_TID_OK: IdentMap.find tid' (threads c2) = Some (existT _ lang st2, lc2))
      (C2: c2 = mk (threads c2) tid' (sc c1) (memory c1)):
      aux_step AuxEvent.tterm lo c1 c2
  .
  Hint Constructors aux_step.

  Inductive aux_tau_step: Ordering.LocOrdMap -> t -> t -> Prop :=
  | aux_tau_step_intro
      e lo c1 c2
      (WP_STEP: aux_step e lo c1 c2)
      (TAU: ~(exists e0, e = AuxEvent.out e0)):
    aux_tau_step lo c1 c2.

  Inductive aux_all_step: Ordering.LocOrdMap -> t -> t -> Prop :=
  | aux_all_step_intro
      e lo c1 c2
      (WP_STEP: aux_step e lo c1 c2):
    aux_all_step lo c1 c2.

End Configuration.

author: Hughshine

repo: promising-comp


text:
Require Import Auxiliaries.
Require Import Nat List Bool Arith Omega.

From QuickChick Require Import QuickChick.
Import ListNotations.
Import QcDefaultNotation. Open Scope qc_scope.
Import GenLow GenHigh.

Require Import Coq.Structures.Equalities.

(**
  Just to be sure, the correctness of some of our definitions and lemmas are
  "proved" with QuickChick
*)

Module NatType <: Typ with Definition t := nat.
Definition t := nat.
End NatType.

Module InserterRemoverOverNat := (NthInserterRemover NatType).
Import InserterRemoverOverNat.

Definition genInput : G (nat * (list nat)) :=
  genPair (choose (0, 100)) (listOf (choose (0, 100))).

Definition nth_remove_length_check n (l : list nat) :=
  ((n <? length l) && ((length (nth_remove n l)) =? (pred (length l))))
  ||
  ((length l <=? n) && ((length (nth_remove n l)) =? (length l))).

Definition nth_insert_length_check n (l: list nat) x :=
  ((n <=? length l) && (length (nth_insert n l x) =? (succ (length l))))
  ||
  ((length l <? n) && ((length (nth_insert n l x)) =? (length l))).

Definition nth_insert_prop n (l : list nat) x :=
  ((n <=? length l) &&
   (0 <? length (filter (fun x' => (x' =? x)) (nth_insert n l x))))
  || (length l <? n).

Definition qc_test_remove_length (nl: prod nat (list nat)) :=
  nth_remove_length_check (fst nl) (snd nl).
Definition qc_test_insert_length (nl: prod nat (list nat)) :=
  nth_insert_length_check (fst nl) (snd nl).
Definition qc_test_insert_prop (nl: prod nat (list nat)) :=
  nth_insert_prop (fst nl) (snd nl).

QuickChick (forAll genInput qc_test_remove_length).
QuickChick (forAll genInput qc_test_insert_length).
QuickChick (forAll genInput qc_test_insert_prop).

author: holmuk

repo: coq-zipper

save_path: github-repos/coq/holmuk-coq-zipper


text:
From Coq Require Import Extraction.
From Coq Require Import List.
Declare ML Module "coq_spec".
From Hammer Require Import Hammer.
From Hammer Require Import Reconstr.
Add LoadPath "~/Documents/cpdt/src" as Cpdt.
Require Import Cpdt.CpdtTactics.
Open Scope list_scope.

MaxTermSize 5.
MaxTestSize 20.

From Coq Require Import ZArith.BinIntDef.
Module ZZ := BinIntDef.Z.
Print ZZ.

Extract Inductive bool => "Prelude.Bool" [ "Prelude.True" "Prelude.False" ].
Extract Inductive option => "Prelude.Maybe" [ "Prelude.Just" "Prelude.Nothing" ].
Extract Inductive unit => "()" [ "()" ].
Extract Inductive list => "([])" [ "([])" "(:)" ].
Extract Inductive prod => "(,)" [ "(,)" ].

Extract Inductive sumbool => "Prelude.Bool" [ "Prelude.True" "Prelude.False" ].
Extract Inductive sumor => "Prelude.Maybe" [ "Prelude.Just" "Prelude.Nothing" ].
Extract Inductive sum => "Prelude.Either" [ "Prelude.Left" "Prelude.Right" ].

Extract Inlined Constant andb => "(Prelude.&&)".
Extract Inlined Constant orb => "(Prelude.||)".
Extract Inlined Constant negb => "Prelude.not".

(* not needed to prove the laws of groups, but can be useful for clarity *)
Definition isUnit (n : Z) := ZZ.eqb n ZZ.zero.

Definition zadd := ZZ.add.

(* DiscoverLemmas "ZGroupTest" zadd isUnit. *)

(* result with zadd:
== Laws ==
  1. zadd x y = zadd y x
  2. zadd x Z0 = x
  3. zadd (zadd x y) z = zadd x (zadd y z)
  4. zadd (Zneg x) (Zneg x) = Zneg (XO x)
  5. zadd (Zneg x) (Zpos x) = Z0
  6. zadd (Zpos x) (Zpos x) = Zpos (XO x)

  3. is associativity law for group
  2. and 1. implies identity element law for group
  5. and 1. implies inverse element law for group
  Note: the Closure law doesn't make sense here because it is internalized in the underlying type theory.
 *)

(* results with isUnit included:
  1. isUnit x & isUnit y => x = y
  2. isUnit x => Z0 = x
  3. zadd x y = zadd y x
  4. isUnit y => zadd x y = x
  5. isUnit (Zneg x) = False
  6. isUnit (Zpos x) = False
  7. isUnit (zadd x x) = isUnit x
  8. zadd (zadd x y) z = zadd x (zadd y z)
  9. zadd (Zneg x) (Zneg x) = Zneg (XO x)
 10. isUnit y => zadd (Zneg x) (Zpos x) = y
 11. zadd (Zpos x) (Zpos x) = Zpos (XO x)

  1. is the law of uniqueness of unit
  4. is another way of formulating the identity law.
  10. is another way of formulating the inverse element law.
*)

(*  Lemma lemma_1 :
  forall (x : Z) (y : Z),
  (zadd y x) = (zadd x y).
  Proof. crush. Qed.
Lemma lemma_2 :
  forall (x : Z),
  (zadd x Z0) = x.
  Proof. crush. Qed.
Lemma lemma_3 :
  forall (x : Z) (y : Z) (z : Z),
  (zadd (zadd x y) z) = (zadd x (zadd y z)).
  Proof. crush. Qed.
Lemma lemma_4 :
  forall (x : positive),
  (zadd (Zneg x) (Zneg x)) = (Zneg (xO x)).
  Proof. hammer. Qed.
Lemma lemma_5 :
  forall (x : positive),
  (zadd (Zneg x) (Zpos x)) = Z0.
  Proof. crush. hammer. Qed.
Lemma lemma_6 :
  forall (x : positive),
  (zadd (Zpos x) (Zpos x)) = (Zpos (xO x)).
  Proof. crush. hammer. Qed. *)

author: mikkelmilo

repo: rooster-spec


text:
Require Import UniMath.Foundations.All.
Require Import UniMath.MoreFoundations.All.
Require Import UniMath.CategoryTheory.Core.Categories.
Require Import UniMath.CategoryTheory.Core.Functors.
Require Import UniMath.CategoryTheory.Core.Isos.
Require Import UniMath.CategoryTheory.Monoidal.AlternativeDefinitions.MonoidalCategoriesCurried.

Local Open Scope cat.

Section local_helper_lemmas.
  (* I would assume that the following lemmas should already exists in the "iso file", but I can't find it (I need is_iso_stable_undertransportation)  *)
  Lemma iso_stable_under_equality {C : category} {x y : C} {f g : C⟦x,y⟧} : (g = f) → (is_z_isomorphism f) → (is_z_isomorphism g).
  Proof.
    intros pe pi.
    induction pe.
    exact pi.
  Qed.

  Lemma iso_stable_under_tranportation {C : category} {x y z : C} {f : C⟦x,y⟧} {pf : y=z} : (is_z_isomorphism f) → (is_z_isomorphism (transportf _ pf f)).
  Proof.
    intro pfi.
    induction pf.
    use pfi.
  Qed.

  Lemma iso_stable_under_equalitytransportation {C : category} {x y z : C} {f : C⟦x,y⟧} {g : C⟦x,z⟧} {pf : y=z} :
    (g = transportf _ pf f) -> (is_z_isomorphism f) -> (is_z_isomorphism g).
  Proof.
    intros p isof.
    use (iso_stable_under_equality p).
    use (iso_stable_under_tranportation).
    exact isof.
  Qed.
End local_helper_lemmas.

Section Curried_Monoidal_Functors.

  Context (C D : category) (M : monoidalcategory_data C) (N : monoidalcategory_data D) (F : functor C D).

  Notation "x ⊗_{ M } y" := ((tensoronobjects_from_tensordata M) x y) (at level 31).
  Notation "f ⊗^{ M } g" := ((tensoronmorphisms_from_tensordata M) _ _ _ _ f g) (at level 31).
  Notation "I_{ M }" := (unit_from_monoidalcatdata M).
  Notation "lu^{ M }_{ x }" := (leftunitordata_from_monoidalcatdata M x ).
  Notation "ru^{ M }_{ x }" := (rightunitordata_from_monoidalcatdata M x ).
  Notation "α^{ M }_{ x , y , z }" := ((associatordata_from_monoidalcatdata M) x y z).

  (** (Weak) Monoidal functors **)
  (* Monoidal functor data *)
  Definition tensor_preserving_data : UU
    := ∏ (x y : C), D⟦(F x) ⊗_{N} (F y), F (x ⊗_{M} y)⟧.

  Definition unit_preserving_data : UU := D ⟦ I_{N} , F I_{M} ⟧.

  Definition monoidalfunctor_data :=
    tensor_preserving_data × unit_preserving_data.

  Definition tensorpreservingdata_from_monoidalfunctordata (mfd : monoidalfunctor_data) : tensor_preserving_data := pr1 mfd.
  Coercion tensorpreservingdata_from_monoidalfunctordata : monoidalfunctor_data >-> tensor_preserving_data.

  Definition unitpreservingdata_from_monoidalfunctordata (mfd : monoidalfunctor_data) : unit_preserving_data := pr2 mfd.
  Coercion unitpreservingdata_from_monoidalfunctordata : monoidalfunctor_data >-> unit_preserving_data.

  (* Weak monoidal functor properties *)
  Definition tensor_preserving_data_is_natural (tpd : tensor_preserving_data) :=
    ∏ (x y x' y' : C) (f : C⟦x,x'⟧) (g : C⟦y,y'⟧),
      (#F (f ⊗^{M} g))∘(tpd x y) = (tpd x' y')∘((#F f) ⊗^{N} (#F g)).

  Definition preserves_associativity (tpd : tensor_preserving_data) :=
    ∏ (x y z : C), (#F (α^{M}_{x,y,z})) ∘ (tpd (x ⊗_{M} y) z) ∘ ((tpd x y) ⊗^{N} (identity (F z)))
                   = (tpd x (y ⊗_{M} z)) ∘ ((identity (F x)) ⊗^{N} (tpd y z)) ∘ α^{N}_{F x, F y, F z}.

  Definition preserves_leftunitality (tpd : tensor_preserving_data) (upd : unit_preserving_data) : UU :=
    ∏ (x : C), ((#F (lu^{M}_{x})) ∘ (tpd I_{M} x) ∘ (upd ⊗^{N} (identity (F x)))) = lu^{N}_{F x}.

  Definition preserves_rightunitality (tpd : tensor_preserving_data) (upd : unit_preserving_data) : UU :=
    ∏ (x : C), (#F (ru^{M}_{x})) ∘ (tpd x I_{M}) ∘ ((identity (F x)) ⊗^{N} upd) = ru^{N}_{F x}.

  Definition monoidalfunctor_laws (mfd : monoidalfunctor_data) : UU :=
    (tensor_preserving_data_is_natural mfd) × (preserves_associativity mfd) ×
    (preserves_leftunitality mfd mfd) × (preserves_rightunitality mfd mfd).

  Definition monoidalfunctor : UU :=
    ∑ (mfd : monoidalfunctor_data), monoidalfunctor_laws mfd.

  (** Strong and strict monoidal properties *)

  Definition is_stronglytensorpreserving (tpd : tensor_preserving_data) : UU :=
    ∏ (x y : C), is_z_isomorphism (tpd x y).

  Definition is_strictlytensorpreserving (tpd : tensor_preserving_data) : UU :=
      ∏ (x y : C), ∑ (pf : (F x) ⊗_{N} (F y) = F (x ⊗_{M} y)),
      (tpd x y ) = transportf _ pf (identity ((F x) ⊗_{N} (F y))).

  Lemma strictlytensorpreserving_is_strong {tpd : tensor_preserving_data} (pfstrict : is_strictlytensorpreserving tpd) : is_stronglytensorpreserving tpd.
  Proof.
    intros x y.
    use (iso_stable_under_equalitytransportation (pr2 (pfstrict x y)) (is_z_isomorphism_identity ((F x) ⊗_{N} (F y)))).
  Defined.
  Coercion strictlytensorpreserving_is_strong : is_strictlytensorpreserving >-> is_stronglytensorpreserving.

  Definition is_stronglyunitpreserving (upd : unit_preserving_data) : UU := is_z_isomorphism upd.

  Definition is_strictlyunitpreserving (upd : unit_preserving_data) : UU :=
    ∑ (pf : I_{N} = (F I_{M})), upd = transportf _ pf (identity I_{N}).

  Definition strictlyunitpreserving_is_strong {upd : unit_preserving_data} (pfstrict : is_strictlyunitpreserving upd) : is_stronglyunitpreserving upd.
  Proof.
    use (iso_stable_under_equalitytransportation (pr2 pfstrict) (is_z_isomorphism_identity I_{N})).
  Defined.
  Coercion strictlyunitpreserving_is_strong : is_strictlyunitpreserving >-> is_stronglyunitpreserving.

End Curried_Monoidal_Functors.

author: UniMath

repo: UniMath

save_path: github-repos/coq/UniMath-UniMath


text:
Definition Seq{A:Type} := nat -> A.

Definition sin{A : Type}(a : A)(q: Seq(A:=A)) := exists i, a = q i.

Lemma sin_nth{A:Type}(q: Seq(A:=A)): forall n,  sin (q n) q.
Proof.
intro n.
now exists n.
Qed.

Definition increasing {A: Type} (R: A -> A -> Prop)(q : Seq) : Prop  :=
  forall n,  R (q n) (q (S n)).

Definition strictly_increasing{A : Type}(R: A -> A -> Prop)(q : Seq(A := A)) :=
  increasing R q /\
forall n,  ~R (q (S n)) (q n).

Lemma increasing_snth{A : Type}
      (R : A -> A -> Prop)
      (RRefl : forall a,  R a a)
      (Rtrans : forall a b c,  R a b -> R b c -> R a c)
  : forall (q: Seq (A := A)),
    increasing R q  -> forall (i j : nat),  i <= j -> R (q i) (q j).
Proof.
intros q Hr i j Hle.
induction Hle.
*
  apply RRefl.
*
  specialize (Hr m).
  eapply Rtrans  ; eauto.
Qed.

Definition sup{A: Type} (R : A -> A-> Prop) (sup_val: A) (q : Seq(A:=A)) :=
  (forall a, sin a q -> R a sup_val) /\
  (forall sup_val',  (forall a, sin a q -> R a sup_val')  -> R sup_val sup_val').

author: hidden-author


text:
Load LFindLoad.
From lfind Require Import LFind.
From QuickChick Require Import QuickChick.
From adtind Require Import goal33.

Derive Show for natural.

              Derive Arbitrary for natural.

              Instance Dec_Eq_natural : Dec_Eq natural.

              Proof. dec_eq. Qed.

Lemma conj1synthconj6 : forall  (lv0 : natural) (lv1 : natural), (@eq natural (lv0) lv1).
Admitted.
QuickChick conj1synthconj6.

author: yalhessi

repo: lemmaranker

save_path: github-repos/coq/yalhessi-lemmaranker


(** * Linearize: Place all and only operations in let binders *)
Require Import Crypto.Compilers.Syntax.
Require Import Crypto.Compilers.Wf.
Require Import Crypto.Compilers.ExprInversion.
Require Import Crypto.Compilers.WfInversion.
Require Import Crypto.Compilers.WfProofs.
Require Import Crypto.Compilers.Linearize.
Require Import Crypto.Util.Sigma.
Require Import Crypto.Util.Tactics.BreakMatch.
Require Import Crypto.Util.Tactics.DestructHead.

Local Open Scope ctype_scope.
Section language.
  Context {base_type_code : Type}
          {op : flat_type base_type_code -> flat_type base_type_code -> Type}.

  Local Notation flat_type := (flat_type base_type_code).
  Local Notation type := (type base_type_code).
  Local Notation Tbase := (@Tbase base_type_code).
  Local Notation exprf := (@exprf base_type_code op).
  Local Notation expr := (@expr base_type_code op).
  Local Notation Expr := (@Expr base_type_code op).
  Local Notation wff := (@wff base_type_code op).
  Local Notation wf := (@wf base_type_code op).

  Section with_var.
    Context {var1 var2 : base_type_code -> Type}.

    Local Ltac t_fin_step tac :=
      match goal with
      | _ => assumption
      | _ => progress simpl in *
      | _ => progress subst
      | _ => progress inversion_sigma
      | _ => progress destruct_head'_sig
      | _ => progress destruct_head'_and
      | _ => setoid_rewrite List.in_app_iff
      | [ H : context[List.In _ (_ ++ _)] |- _ ] => setoid_rewrite List.in_app_iff in H
      | _ => progress intros
      | _ => solve [ eauto ]
      | _ => solve [ intuition (subst; eauto) ]
      | [ H : forall (x : prod _ _) (y : prod _ _), _ |- _ ] => specialize (fun x x' y y' => H (x, x') (y, y'))
      | _ => rewrite !List.app_assoc
      | [ H : _ \/ _ |- _ ] => destruct H
      | [ H : _ |- _ ] => apply H
      | [ H : forall G : list _, _ |- _ ] => apply (H nil)
      | _ => eapply wff_in_impl_Proper; [ solve [ eauto using wff_SmartVarf ] | solve [ repeat t_fin_step tac ] ]
      | _ => progress tac
      | [ |- wff _ _ _ ] => constructor
      | [ |- wf _ _ _ ] => constructor
      | _ => break_innermost_match_step
      | _ => progress inversion_expr
      | _ => congruence
      | _ => progress destruct_head' sum
      | _ => progress unfold invert_Op in *
      | _ => break_innermost_match_hyps_step
      end.
    Local Ltac t_fin tac := repeat t_fin_step tac.

    Local Hint Constructors Wf.wff.
    Local Hint Resolve List.in_app_or List.in_or_app.

    Section gen1.
      Context (let_bind_op_args : bool).

      Fixpoint wff_under_letsf' G {t} e1 e2 {tC} eC1 eC2
               let_bind_pairs
               (wf : @wff var1 var2 G t e1 e2)
               (H : forall (x1 : interp_flat_type var1 t) (x2 : interp_flat_type var2 t),
                   wff (flatten_binding_list x1 x2 ++ G) (eC1 (inl x1)) (eC2 (inl x2)))
               (H' : forall G' (x y : exprf t),
                   wff (G' ++ G) x y
                   -> wff (G' ++ G) (eC1 (inr x)) (eC2 (inr y)))
               {struct e1}
        : @wff var1 var2 G tC (under_letsf' let_bind_op_args let_bind_pairs e1 eC1) (under_letsf' let_bind_op_args let_bind_pairs e2 eC2).
      Proof using Type.
        revert H.
        set (e1v := e1) in *.
        destruct e1 as [ | | ? ? ? args | tx ex tC0 eC0 | ? ex ? ey ];
          [ clear wff_under_letsf'
          | clear wff_under_letsf'
          | clear wff_under_letsf'
          | generalize (fun G => match e1v return match e1v with LetIn _ _ _ _ => _ | _ => _ end with
                                 | LetIn _ ex _ eC => wff_under_letsf' G _ ex
                                 | _ => I
                                 end);
            generalize (fun G => match e1v return match e1v with LetIn tx0 _ tC1 e0 => _ | _ => _ end with
                                 | LetIn _ ex tC' eC => fun x => wff_under_letsf' G tC' (eC x)
                                 | _ => I
                                 end);
            clear wff_under_letsf'
          | generalize (fun G => match e1v return match e1v with Pair _ _ _ _ => _ | _ => _ end with
                                 | Pair _ ex _ ey => wff_under_letsf' G _ ex
                                 | _ => I
                                 end);
            generalize (fun G => match e1v return match e1v with Pair _ _ _ _ => _ | _ => _ end with
                                 | Pair _ ex _ ey => wff_under_letsf' G _ ey
                                 | _ => I
                                 end);
            clear wff_under_letsf' ];
          subst e1v;
          cbv beta iota;
          (invert_one_expr e2 || destruct e2); intros; try break_innermost_match_step; try exact I; intros;
            inversion_wf;
            t_fin idtac.
      Qed.

      Lemma  wff_under_letsf G {t} e1 e2 {tC} eC1 eC2
             (wf : @wff var1 var2 G t e1 e2)
             (H' : forall G' (x y : exprf t),
                 wff (G' ++ G) x y
                 -> wff (G' ++ G) (eC1 x) (eC2 y))
        : @wff var1 var2 G tC (under_letsf let_bind_op_args e1 eC1) (under_letsf let_bind_op_args e2 eC2).
      Proof using Type.
        apply wff_under_letsf'; t_fin idtac.
      Qed.
    End gen1.

    Local Hint Resolve wff_under_letsf.
    Local Hint Constructors or.
    Local Hint Extern 1 => progress unfold List.In in *.
    Local Hint Resolve wff_in_impl_Proper.
    Local Hint Resolve wff_SmartVarf.

    Section gen2.
      Context (let_bind_op_args : bool).

      Lemma wff_linearizef_gen G {t} e1 e2
        : @wff var1 var2 G t e1 e2
          -> @wff var1 var2 G t (linearizef_gen let_bind_op_args e1) (linearizef_gen let_bind_op_args e2).
      Proof using Type.
        induction 1; t_fin ltac:(apply wff_under_letsf).
      Qed.

      Local Hint Resolve wff_linearizef_gen.

      Lemma wf_linearize_gen {t} e1 e2
        : @wf var1 var2 t e1 e2
          -> @wf var1 var2 t (linearize_gen let_bind_op_args e1) (linearize_gen let_bind_op_args e2).
      Proof using Type.
        destruct 1; constructor; auto.
      Qed.
    End gen2.
  End with_var.

  Section gen.
    Context (let_bind_op_args : bool).

    Lemma Wf_Linearize_gen {t} (e : Expr t) : Wf e -> Wf (Linearize_gen let_bind_op_args e).
    Proof using Type.
      intros wf var1 var2; apply wf_linearize_gen, wf.
    Qed.
  End gen.

  Definition wff_linearizef {var1 var2} G {t} e1 e2
    : @wff var1 var2 G t e1 e2
      -> @wff var1 var2 G t (linearizef e1) (linearizef e2)
    := wff_linearizef_gen _ G e1 e2.
  Definition wff_a_normalf {var1 var2} G {t} e1 e2
    : @wff var1 var2 G t e1 e2
      -> @wff var1 var2 G t (a_normalf e1) (a_normalf e2)
    := wff_linearizef_gen _ G e1 e2.

  Definition wf_linearize {var1 var2 t} e1 e2
    : @wf var1 var2 t e1 e2
      -> @wf var1 var2 t (linearize e1) (linearize e2)
    := wf_linearize_gen _ e1 e2.
  Definition wf_a_normal {var1 var2 t} e1 e2
    : @wf var1 var2 t e1 e2
      -> @wf var1 var2 t (a_normal e1) (a_normal e2)
    := wf_linearize_gen _ e1 e2.

  Definition Wf_Linearize {t} (e : Expr t) : Wf e -> Wf (Linearize e)
    := Wf_Linearize_gen _ e.
  Definition Wf_ANormal {t} (e : Expr t) : Wf e -> Wf (ANormal e)
    := Wf_Linearize_gen _ e.
End language.

Hint Resolve Wf_Linearize_gen Wf_Linearize Wf_ANormal : wf.


text:
Require Export euclidean__axioms.
Require Export lemma__3__7b.
Require Export lemma__betweennotequal.
Require Export lemma__extensionunique.
Require Export lemma__inequalitysymmetric.
Require Export logic.
Definition lemma__partnotequalwhole : forall A B C, (euclidean__axioms.BetS A B C) -> (~(euclidean__axioms.Cong A B A C)).
Proof.
intro A.
intro B.
intro C.
intro H.
assert (* Cut *) (euclidean__axioms.neq A B) as H0.
- assert (* Cut *) ((euclidean__axioms.neq B C) /\ ((euclidean__axioms.neq A B) /\ (euclidean__axioms.neq A C))) as H0.
-- apply (@lemma__betweennotequal.lemma__betweennotequal A B C H).
-- destruct H0 as [H1 H2].
destruct H2 as [H3 H4].
exact H3.
- assert (* Cut *) (euclidean__axioms.neq B A) as H1.
-- apply (@lemma__inequalitysymmetric.lemma__inequalitysymmetric A B H0).
-- assert (* Cut *) (exists D, euclidean__axioms.BetS B A D) as H2.
--- apply (@euclidean__axioms.postulate__Euclid2 B A H1).
--- destruct H2 as [D H3].
assert (* Cut *) (euclidean__axioms.BetS D A B) as H4.
---- apply (@euclidean__axioms.axiom__betweennesssymmetry B A D H3).
---- assert (* Cut *) (euclidean__axioms.BetS D A C) as H5.
----- apply (@lemma__3__7b.lemma__3__7b D A B C H4 H).
----- assert (* Cut *) (euclidean__axioms.neq B C) as H6.
------ assert (* Cut *) ((euclidean__axioms.neq B C) /\ ((euclidean__axioms.neq A B) /\ (euclidean__axioms.neq A C))) as H6.
------- apply (@lemma__betweennotequal.lemma__betweennotequal A B C H).
------- destruct H6 as [H7 H8].
destruct H8 as [H9 H10].
exact H7.
------ assert (* Cut *) (~(euclidean__axioms.Cong A B A C)) as H7.
------- intro H7.
assert (* Cut *) (B = C) as H8.
-------- apply (@lemma__extensionunique.lemma__extensionunique D A B C H4 H5 H7).
-------- apply (@H6 H8).
------- exact H7.
Qed.

author: Karnaj

repo: dktactgeo

save_path: github-repos/coq/Karnaj-dktactgeo