text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = INVERSEKINEMATIC_KR_1000_1300_TITAN(robot, T)
%   Solves the inverse kinematic problem for the KUKA KR 1000 1300 TITAN robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=INVERSEKINEMATIC_KR_1000_1300_TITAN returns 8 possible solutions, thus,
%   Q is a 6x8 matrix where each column stores 6 feasible joint values.
%
%
%   Example code:
%
%   robot=load_robot('kuka', 'KR_1000_1300_TITAN');
%   q = [0 0 0 0 0 0];
%   T = directkinematic(robot, q);
%   %Call the inversekinematic for this robot
%   qinv = inversekinematic(robot, T);
%   check that all of them are feasible solutions!
%   and every Ti equals T
%   for i=1:8,
%        Ti = directkinematic(robot, qinv(:,i))
%   end
%	See also DIRECTKINEMATIC.
%   Authors:   Javier Martinez Gonzalez
%              Jose Francisco Munoz Sempere
%              Silvia Carretero Monasor
%              Marcos Gomez Parres
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.

function q = inversekinematic_KR_1000_1300_TITAN(robot, T)

%initialize q,
%eight possible solutions are generally feasible
q=zeros(6,8);

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

L6=d(6); %Distancia de la mu�eca al efector final.

A1 = a(1);

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Compute the position of the wrist, being W the Z component of the end effector's system
%Tomamos el vector "a"
W = T(1:3,3); %ax, ay, az

% Pm: wrist position (Posici�n de la mu�eca)
Pm = [Px Py Pz]' - L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution
%Obtenemos "theta 1" mediante m�todos geom�tricos.
q1=atan2(-Pm(2), Pm(1));

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);

q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);

q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%Arrange solutions, there are 8 possible solutions so far.
% if q1 is a solution, q1* = q1 + pi is also a solution.
% For each (q1, q1*) there are two possible solutions
% for q2 and q3 (namely, elbow up and elbow up solutions)
% So far, we have 4 possible solutions. Howefer, for each triplet (theta1, theta2, theta3),
% there exist two more possible solutions for the last three joints, generally
% called wrist up and wrist down solutions. For this reason,
%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];

%At this point, we want to asure, that, at least, the function returns
%4 real solutions q(1:4). If any of the solutions q(1:4) is complex, only
%the real part will be returned. If any of the solutions q(5:8) is complex,
%it will not be considered and removed
q = arrange_solutions(q);

%leave only the real part of the solutions
q=real(q);

%Note that in this robot, the joint q3 has a non-simmetrical range. In this
%case, the joint ranges from 60 deg to -219 deg, thus, the typical normalizing
%step is avoided in this angle (the next line is commented). When solving
%for the orientation, the solutions are normalized to the [-pi, pi] range
%only for the theta4, theta5 and theta6 joints.

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));

% % solve for the last three joints
% % for any of the possible combinations (theta1, theta2, theta3)
% %Please note the special orientation of the wrist of this robot. In this
% %case, we employ an 'ad hoc' function to solve for the final orientaion
for i=1:2:size(q,2),
    qtemp = solve_for_last_three_joints(robot, q(:,i), T, 1); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    qtemp = solve_for_last_three_joints(robot, q(:,i), T, -1); %wrist down
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
d = eval(robot.DH.d);
a = eval(robot.DH.a);

%See geometry
%Par�metros con los que calculamos "theta 2" y "theta 3", tambi�n mediante
%m�todos geom�tricos.
%Tenemos en cuenta el desfase de 65mm entre los centros de los sistemas
%de referencia de los eslabones 2 y 3.
L2=abs(a(2));
L3=abs(d(4));
A2 = abs(a(3)); %ESTE ELEMENTO ES EL DESFASE CITADO
%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

L4 = sqrt(A2^2 + L3^2); %LONGITUD DEL ESLAB�N TENIENDO EN CUENTA EL DESFASE
%(DISTANCIA REAL ENTRE MU�ECA Y ESLAB�N 2)
%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2); %r ES LA DISTANCIA DEL SIST. 1 A LA MU�ECA

beta = atan2(-p1(2), p1(1)); %BETA=ARCTG(-Y/X)
gamma = (acos((L2^2+r^2-L4^2)/(2*r*L2))); %TEOREMA DEL COSENO

if ~isreal(gamma)
    disp('WARNING:inversekinematic_KUKA_1000_1300_TITAN: the point is not reachable for this configuration, imaginary solutions');
    %gamma = real(gamma);
end

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta
q2(1) = pi/2-gamma-beta; %elbow up (CODO ARRIBA)
q2(2) = pi/2+gamma-beta; %elbow down (CODO ABAJO)
%DEBEMOS A�ADIR UN �NGULO DE PI/2 DEBIDO A DESFASES EN NUESTROS SISTEMAS DE
%REFERENCIA

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
d = eval(robot.DH.d);
a = eval(robot.DH.a);

%See geometry
L2=abs(a(2));
L3=abs(d(4));
A2 = abs(a(3));

L4 = sqrt(A2^2 + L3^2); %VOLVEMOS A TENER EN CUENTA EL DESFASE

phi=acos((A2^2+L4^2-L3^2)/(2*A2*L4)); %TEOMERA DEL COSENO
%PHI, ES EL �NGULO QUE FORMAN LAS L�NEA QUE DEFINEN EL DESFASE ENTRE 2 Y 3,
%Y EL PROPIO ESLAB�N 4
%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2); %DE NUEVO, DISTANCIA ENTRE SIST.1 Y MU�ECA

eta = (acos((L2^2 + L4^2 - r^2)/(2*L2*L4))); %TEOREMA DEL COSENO
%ETA, ES EL �NGULO REAL ENTRE LOS ESLABONES 3 Y 4

if ~isreal(eta)
   disp('WARNING:inversekinematic_KUKA_1000_1300_TITAN: the point is not reachable for this configuration, imaginary solutions');
   %eta = real(eta);
end

%return two possible solutions
%elbow up and elbow down solutions
%the order here is important
q3(1) = pi - phi - eta; %CODO ARRIBA
q3(2) = pi - phi + eta; %CODO ABAJO
%SE A�ADEN LOS �NGULOS "PI" EN AMBAS SOLUCIONES DEBIDO A DESFASES ENTRE
%NUESTROS SISTEMAS DE REFERENCIA.

% Solve for the last three joints asuming an spherical wrist
function q = solve_for_last_three_joints(robot, q, T, wrist)

% T is the noa matrix defining the position/orientation of the end
% effector's reference system
vx6=T(1:3,1);
vz6=T(1:3,3);
vz5=T(1:3,3); % The vector a T(1:3,3) is coincident with z5

% Obtain the position and orientation of the system 3
% using the already computed joints q1, q2 and q3
T01=dh(robot, q, 1);
T12=dh(robot, q, 2);
T23=dh(robot, q, 3);
T03=T01*T12*T23;

vx3=T03(1:3,1);
vy3=T03(1:3,2);
vz3=T03(1:3,3);

% find z4 normal to the plane formed by z3 and a
vz4=cross(vz3, vz6);	% end effector's vector a: T(1:3,3)

% in case of degenerate solution,
% when vz3 and vz6 are parallel--> then z4=0 0 0, choose q(4)=0 as solution
if norm(vz4) <= 0.00001
    if wrist == 1 %wrist up
        q(4)=0;
    else
        q(4)=-pi; %wrist down
    end
else
    %this is the normal and most frequent solution
    cosq4=wrist*dot(-vy3,vz4);
    sinq4=wrist*dot(vx3,vz4);
    q(4)=atan2(sinq4, cosq4);
end
%propagate the value of q(4) to compute the system 4
T34=dh(robot, q, 4);
T04=T03*T34;
vx4=T04(1:3,1);
vy4=T04(1:3,2);

% solve for q5
cosq5=dot(vx4,vz5);
sinq5=dot(vy4,vz5);
q(5)=atan2(sinq5, cosq5);

%propagate now q(5) to compute T05
T45=dh(robot, q, 5);
T05=T04*T45;
vx5=T05(1:3,1);
vy5=T05(1:3,2);

% solve for q6
cosq6=dot(vx6,vx5);
sinq6=dot(vx6,vy5);
q(6)=atan2(sinq6, cosq6);

%remove complex solutions for q for the q1+pi solutions
function  qreal = arrange_solutions(q)
qreal=q(:,1:4);

%sum along rows if any angle is complex, for any possible solutions, then v(i) is complex
v = sum(q, 1);

for i=5:8,
    if isreal(v(i))
        qreal=[qreal q(:,i)]; %store the real solutions
    end
end

qreal = real(qreal);

repo: ARTE


text:
function [S] = lukfPropagation(dt,chi,chiAnt,omega_b,a_b,S,omega,acc,Qc,g)
%Left-UKF on Lie Groups
q = length(S); % state size
S_aug = blkdiag(S,Qc);
N_aug = length(S_aug);

%scaled unsented transform
W0 = 1-N_aug/3;
Wj = (1-W0)/(2*N_aug);
gamma = sqrt(N_aug/(1-W0));

% Prediction
ichi = invSE3(chi);
omega = omega-omega_b; % unbiased inputs
acc = acc-a_b;
X = gamma*[zeros(N_aug,1) S_aug' -S_aug'];% sigma-points
X(10:15,:) = X(10:15,:)+X(q+7:q+12,:)*dt; % add bias noise
for j = 2:2*N_aug+1
    xi_j = X([1:9,16:q],j); %do not take bias
    w_j = X(q+1:N_aug,j);
    omega_bj = X(10:12,j);
    a_bj = X(13:15,j);
    chi_j = chiAnt*exp_multiSE3(xi_j);
    Rot = chi_j(1:3,1:3)*expSO3((omega+w_j(1:3)-omega_bj)*dt);
    v = chi_j(1:3,4)+(Rot*(acc+w_j(4:6)-a_bj)+g)*dt;
    x = chi_j(1:3,5)+v*dt;
    chi = state2chi(Rot,v,x,chi_j(1:3,6:end));
    Xi_j = ichi*chi; % can be more time efficient
    X([1:9,16:q],j) = log_multiSE3(Xi_j); %propagated sigma points
end
X = sqrt(Wj)*X;
[~,Rs] = qr(X(1:q,2:2*N_aug+1)');
S = Rs(1:q,1:q);
end

author: mbrossar


text:
% Defines a protein complex
%
% Author: Jonathan Karr
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 5/7/2009
classdef ProteinComplex < edu.stanford.covert.cell.kb.PhysicalObject & edu.stanford.covert.cell.kb.Protein
    properties
        prostheticGroups            = edu.stanford.covert.cell.kb.Metabolite.empty(0, 0);
        prostheticGroupCompartments = edu.stanford.covert.cell.kb.Compartment.empty(0, 0);
        prostheticGroupCoefficients = [];

        chaperoneSubstrates         = edu.stanford.covert.cell.kb.ProteinMonomer.empty(0, 0);
        chaperoneCompartments       = edu.stanford.covert.cell.kb.Compartment.empty(0, 0);
        chaperoneCoefficients       = [];

        metabolites                 = edu.stanford.covert.cell.kb.Metabolite.empty(0, 0);
        metaboliteCompartments      = edu.stanford.covert.cell.kb.Compartment.empty(0, 0);
        metaboliteCoefficients      = [];

        rnas                        = edu.stanford.covert.cell.kb.Gene.empty(0, 0);
        rnaCompartments             = edu.stanford.covert.cell.kb.Compartment.empty(0, 0);
        rnaCoefficients             = [];

        proteinMonomers             = edu.stanford.covert.cell.kb.ProteinMonomer.empty(0, 0);
        proteinMonomerCompartments  = edu.stanford.covert.cell.kb.Compartment.empty(0, 0);
        proteinMonomerCoefficients  = [];

        proteinComplexs             = edu.stanford.covert.cell.kb.ProteinComplex.empty(0, 0);
        proteinComplexCompartments  = edu.stanford.covert.cell.kb.Compartment.empty(0, 0);
        proteinComplexCoefficients  = [];

        stableModificationReactions = edu.stanford.covert.cell.kb.Reaction.empty(0, 0);

        regulatedTranscriptionUnits = edu.stanford.covert.cell.kb.TranscriptionUnit.empty(0, 0);

        compartment                 = edu.stanford.covert.cell.kb.Compartment.empty(0, 0);

        stimuliRegulators           = edu.stanford.covert.cell.kb.Stimuli.empty(0, 0);
        metaboliteRegulators        = edu.stanford.covert.cell.kb.Metabolite.empty(0, 0);
        proteinMonomerRegulators    = edu.stanford.covert.cell.kb.ProteinMonomer.empty(0, 0);
        proteinComplexRegulators    = edu.stanford.covert.cell.kb.ProteinComplex.empty(0, 0);

        regulatedProteinMonomers    = edu.stanford.covert.cell.kb.ProteinMonomer.empty(0, 0);
        regulatedProteinComplexs    = edu.stanford.covert.cell.kb.ProteinComplex.empty(0, 0);

        parameters                  = edu.stanford.covert.cell.kb.Parameter.empty(0, 0);

        complexFormationProcess      = edu.stanford.covert.cell.kb.Process.empty(0, 0); %process where complex is formed
    end

    properties %(SetAccess = protected)
        empiricalFormula
        smiles
        charge
        pKa
        halfLife
        activationRule
        disulfideBondMonomers = cell(0, 1);
        disulfideBondCysteines = zeros(0, 2);
    end

    %computed properties
    properties %(SetAccess = protected)
        rnaComposition
        proteinMonomerComposition

        molecularWeight
        density
        volume
        pI
        extinctionCoefficient
        absorbanceFactor

        negAA
        posAA
        instabilityIndex
        stable
        aliphaticIndex
        gravy

        baseCount
        cumulativeBaseCount
        decayReaction
        tRNACount

        matureBaseCount
        matureMolecularWeight
        matureHalfLife
        matureDecayReaction

        numDistinctSubunits
        numSubunits

        numNucleicAcids
        numAminoAcids
        numRNASubunits
        numProteinMonomerSubunits
    end

    methods
        function this = ProteinComplex(knowledgeBase, wid, wholeCellModelID, name,...
                dnaFootprint, dnaFootprintBindingStrandedness, dnaFootprintRegionStrandedness, molecularInteraction, chemicalRegulation, ...
                subsystem, generalClassification, proteaseClassification, ...
                transporterClassification, ...
                activationRule,...
                disulfideBonds,...
                comments, crossReferences)

            if nargin == 0; return; end;

            this = edu.stanford.covert.cell.kb.ProteinComplex.empty(size(wid, 1), 0);
            this(size(wid, 1), 1) = edu.stanford.covert.cell.kb.ProteinComplex;
            for i = 1:size(wid, 1)
                this(i, 1).idx = i;
                this(i, 1).knowledgeBase = knowledgeBase;
                this(i, 1).wid = wid(i);
                this(i, 1).wholeCellModelID = wholeCellModelID{i};
                this(i, 1).name = name{i};
                if exist('comments','var') && ~isempty(comments); this(i,1).comments = comments{i}; end;
                if exist('crossReferences','var')
                    if size(crossReferences, 1) > 1
                        this(i, 1).crossReferences = crossReferences(i);
                    else
                        this(i, 1).crossReferences = struct;
                        fields = fieldnames(crossReferences);
                        for j = 1:size(fields, 1)
                            values = crossReferences.(fields{j});
                            this(i,1).crossReferences.(fields{j}) = values(i);
                        end
                    end
                end

                this(i, 1).dnaFootprint = dnaFootprint(i);
                this(i, 1).dnaFootprintBindingStrandedness = dnaFootprintBindingStrandedness{i};
                this(i, 1).dnaFootprintRegionStrandedness = dnaFootprintRegionStrandedness{i};
                this(i, 1).molecularInteraction = molecularInteraction{i};
                this(i, 1).chemicalRegulation = chemicalRegulation{i};
                this(i, 1).subsystem = subsystem{i};
                this(i, 1).generalClassification = generalClassification{i};
                this(i, 1).proteaseClassification = proteaseClassification{i};
                this(i, 1).transporterClassification = transporterClassification{i};
                this(i, 1).activationRule = activationRule{i};

                if disulfideBonds{i}
                    disulfideBondsArr = strsplit(';', disulfideBonds{i});

                    this(i, 1).disulfideBondMonomers = cell(length(disulfideBondsArr), 1);
                    this(i, 1).disulfideBondCysteines = zeros(length(disulfideBondsArr), 2);
                    for j = 1:length(disulfideBondsArr)
                        tmp = strsplit(': ', disulfideBondsArr{j});
                        this(i, 1).disulfideBondMonomers{j, 1} = tmp{1};
                        tmp = strsplit('-', tmp{2});
                        this(i, 1).disulfideBondCysteines(j, :) = [str2double(tmp{1}(2:end)) str2double(tmp{2}(2:end))];
                    end
                end
            end
        end

        function serializeLinks(this)
            for i = 1:numel(this)
                this(i).prostheticGroups            = this.serializeLinksHelper(this(i).prostheticGroups);
                this(i).prostheticGroupCompartments = this.serializeLinksHelper(this(i).prostheticGroupCompartments);
                this(i).prostheticGroupCoefficients = this.serializeLinksHelper(this(i).prostheticGroupCoefficients);

                this(i).chaperoneSubstrates         = this.serializeLinksHelper(this(i).chaperoneSubstrates);
                this(i).chaperoneCompartments       = this.serializeLinksHelper(this(i).chaperoneCompartments);
                this(i).chaperoneCoefficients       = this.serializeLinksHelper(this(i).chaperoneCoefficients);

                this(i).metabolites                 = this.serializeLinksHelper(this(i).metabolites);
                this(i).metaboliteCompartments      = this.serializeLinksHelper(this(i).metaboliteCompartments);
                this(i).metaboliteCoefficients      = this.serializeLinksHelper(this(i).metaboliteCoefficients);

                this(i).rnas                        = this.serializeLinksHelper(this(i).rnas);
                this(i).rnaCompartments             = this.serializeLinksHelper(this(i).rnaCompartments);
                this(i).rnaCoefficients             = this.serializeLinksHelper(this(i).rnaCoefficients);

                this(i).proteinMonomers             = this.serializeLinksHelper(this(i).proteinMonomers);
                this(i).proteinMonomerCompartments  = this.serializeLinksHelper(this(i).proteinMonomerCompartments);
                this(i).proteinMonomerCoefficients  = this.serializeLinksHelper(this(i).proteinMonomerCoefficients);

                this(i).proteinComplexs             = this.serializeLinksHelper(this(i).proteinComplexs);
                this(i).proteinComplexCompartments  = this.serializeLinksHelper(this(i).proteinComplexCompartments);
                this(i).proteinComplexCoefficients  = this.serializeLinksHelper(this(i).proteinComplexCoefficients);

                this(i).stableModificationReactions = this.serializeLinksHelper(this(i).stableModificationReactions);

                this(i).regulatedTranscriptionUnits = this.serializeLinksHelper(this(i).regulatedTranscriptionUnits);

                this(i).compartment                 = this.serializeLinksHelper(this(i).compartment);

                this(i).stimuliRegulators           = this.serializeLinksHelper(this(i).stimuliRegulators);
                this(i).metaboliteRegulators        = this.serializeLinksHelper(this(i).metaboliteRegulators);
                this(i).proteinMonomerRegulators    = this.serializeLinksHelper(this(i).proteinMonomerRegulators);
                this(i).proteinComplexRegulators    = this.serializeLinksHelper(this(i).proteinComplexRegulators);

                this(i).regulatedProteinMonomers    = this.serializeLinksHelper(this(i).regulatedProteinMonomers);
                this(i).regulatedProteinComplexs    = this.serializeLinksHelper(this(i).regulatedProteinComplexs);

                this(i).parameters                  = this.serializeLinksHelper(this(i).parameters);

                this(i).complexFormationProcess     = this.serializeLinksHelper(this(i).complexFormationProcess);

                serializeLinks@edu.stanford.covert.cell.kb.PhysicalObject(this(i));
                serializeLinks@edu.stanford.covert.cell.kb.Protein(this(i));
            end
        end

        function deserializeLinks(this, kb)
            for i = 1:numel(this)
                this(i).prostheticGroups            = this.deserializeLinksHelper(this(i).prostheticGroups, kb.metabolites);
                this(i).prostheticGroupCompartments = this.deserializeLinksHelper(this(i).prostheticGroupCompartments, kb.compartments);
                this(i).prostheticGroupCoefficients = this.deserializeLinksHelper(this(i).prostheticGroupCoefficients);

                this(i).chaperoneSubstrates         = this.deserializeLinksHelper(this(i).chaperoneSubstrates, kb.proteinMonomers);
                this(i).chaperoneCompartments       = this.deserializeLinksHelper(this(i).chaperoneCompartments, kb.compartments);
                this(i).chaperoneCoefficients       = this.deserializeLinksHelper(this(i).chaperoneCoefficients);

                this(i).metabolites                 = this.deserializeLinksHelper(this(i).metabolites, kb.metabolites);
                this(i).metaboliteCompartments      = this.deserializeLinksHelper(this(i).metaboliteCompartments, kb.compartments);
                this(i).metaboliteCoefficients      = this.deserializeLinksHelper(this(i).metaboliteCoefficients);

                this(i).rnas                        = this.deserializeLinksHelper(this(i).rnas, kb.genes);
                this(i).rnaCompartments             = this.deserializeLinksHelper(this(i).rnaCompartments, kb.compartments);
                this(i).rnaCoefficients             = this.deserializeLinksHelper(this(i).rnaCoefficients);

                this(i).proteinMonomers             = this.deserializeLinksHelper(this(i).proteinMonomers, kb.proteinMonomers);
                this(i).proteinMonomerCompartments  = this.deserializeLinksHelper(this(i).proteinMonomerCompartments, kb.compartments);
                this(i).proteinMonomerCoefficients  = this.deserializeLinksHelper(this(i).proteinMonomerCoefficients);

                this(i).proteinComplexs             = this.deserializeLinksHelper(this(i).proteinComplexs, kb.proteinComplexs);
                this(i).proteinComplexCompartments  = this.deserializeLinksHelper(this(i).proteinComplexCompartments, kb.compartments);
                this(i).proteinComplexCoefficients  = this.deserializeLinksHelper(this(i).proteinComplexCoefficients);

                this(i).stableModificationReactions = this.deserializeLinksHelper(this(i).stableModificationReactions, kb.reactions);

                this(i).regulatedTranscriptionUnits = this.deserializeLinksHelper(this(i).regulatedTranscriptionUnits, kb.transcriptionUnits);

                this(i).compartment                 = this.deserializeLinksHelper(this(i).compartment, kb.compartments);

                this(i).stimuliRegulators           = this.deserializeLinksHelper(this(i).stimuliRegulators, kb.stimulis);
                this(i).metaboliteRegulators        = this.deserializeLinksHelper(this(i).metaboliteRegulators, kb.metabolites);
                this(i).proteinMonomerRegulators    = this.deserializeLinksHelper(this(i).proteinMonomerRegulators, kb.proteinMonomers);
                this(i).proteinComplexRegulators    = this.deserializeLinksHelper(this(i).proteinComplexRegulators, kb.proteinComplexs);

                this(i).regulatedProteinMonomers    = this.deserializeLinksHelper(this(i).regulatedProteinMonomers, kb.proteinMonomers);
                this(i).regulatedProteinComplexs    = this.deserializeLinksHelper(this(i).regulatedProteinComplexs, kb.proteinComplexs);

                this(i).parameters                  = this.deserializeLinksHelper(this(i).parameters, kb.parameters);

                this(i).complexFormationProcess     = this.deserializeLinksHelper(this(i).complexFormationProcess, kb.processes);

                deserializeLinks@edu.stanford.covert.cell.kb.PhysicalObject(this(i), kb);
                deserializeLinks@edu.stanford.covert.cell.kb.Protein(this(i), kb);
            end
        end

        function deleteLinks(this)
            for i = 1:numel(this)
                this(i).prostheticGroups            = [];
                this(i).prostheticGroupCompartments = [];
                this(i).prostheticGroupCoefficients = [];

                this(i).chaperoneSubstrates         = [];
                this(i).chaperoneCompartments       = [];
                this(i).chaperoneCoefficients       = [];

                this(i).metabolites                 = [];
                this(i).metaboliteCompartments      = [];
                this(i).metaboliteCoefficients      = [];

                this(i).rnas                        = [];
                this(i).rnaCompartments             = [];
                this(i).rnaCoefficients             = [];

                this(i).proteinMonomers             = [];
                this(i).proteinMonomerCompartments  = [];
                this(i).proteinMonomerCoefficients  = [];

                this(i).proteinComplexs             = [];
                this(i).proteinComplexCompartments  = [];
                this(i).proteinComplexCoefficients  = [];

                this(i).stableModificationReactions = [];

                this(i).regulatedTranscriptionUnits = [];

                this(i).compartment                 = [];

                this(i).stimuliRegulators           = [];
                this(i).metaboliteRegulators        = [];
                this(i).proteinMonomerRegulators    = [];
                this(i).proteinComplexRegulators    = [];

                this(i).regulatedProteinMonomers    = [];
                this(i).regulatedProteinComplexs    = [];

                this(i).parameters                  = [];

                deleteLinks@edu.stanford.covert.cell.kb.PhysicalObject(this(i));
                deleteLinks@edu.stanford.covert.cell.kb.Protein(this(i));
            end
        end

        function value = get.empiricalFormula(this)
            throw(MException('ProteinComplex:error', 'property is not defined'));
        end

        function value = get.smiles(this)
            throw(MException('ProteinComplex:error', 'property is not defined'));
        end

        function value = get.charge(this)
            throw(MException('ProteinComplex:error', 'property is not defined'));
        end

        function value = get.pKa(this)
            throw(MException('ProteinComplex:error', 'property is not defined'));
        end

        function value = get.halfLife(this)
            %retrieve
            if ~isempty(this.halfLife)
                value = this.halfLife;
                return;
            end

            %calculate
            %rnaHalfLives = zeros(1,length(this.rnas));
            %for i = 1:length(this.rnas)
            %    rnaHalfLives(i) = this.rnas(i).halfLife;
            %end

            monomerHalfLives = zeros(1,length(this.proteinMonomers));
            for i = 1:length(this.proteinMonomers)
                monomerHalfLives(i) = this.proteinMonomers(i).matureSequenceHalfLife;
            end

            complexHalfLives = zeros(1,length(this.proteinComplexs));
            for i = 1:length(this.proteinComplexs)
                complexHalfLives(i) = this.proteinComplexs(i).matureHalfLife;
            end

            halfLives = [monomerHalfLives complexHalfLives];
            coefficients = [this.proteinMonomerCoefficients; this.proteinComplexCoefficients];

            if ~isempty(coefficients)
                value = halfLives * coefficients / sum(coefficients);
            else
                value = Inf;
            end

            %store
            this.halfLife = value;
        end

        function value = get.matureHalfLife(this)
            %retrieve
            if ~isempty(this.matureHalfLife)
                value = this.matureHalfLife;
                return;
            end

            %calculate
            value = this.halfLife;

            %store
            this.matureHalfLife = value;
        end

        %leave on composition tree
        function value = get.rnaComposition(this)
            %retrieve
            if ~isempty(this.rnaComposition)
                value = this.rnaComposition;
                return;
            end

            %calculate
            value = this.rnaComposition_Helper();

            %store
            this.rnaComposition = value;
        end

        %helper method because MATLAB won't let you call another objects
        %get.rnaComposition method from within get.rnaComposition
        function value = rnaComposition_Helper(this)
            value = zeros(1, this.knowledgeBase.numGenes, this.knowledgeBase.numCompartments);

            if ~isempty(this.rnas)
                value(sub2ind(size(value),...
                    ones(1, numel(this.rnas)),...
                    double([this.rnas.idx]),...
                    double([this.rnaCompartments.idx]))) = ...
                    this.rnaCoefficients;
            end

            if ~isempty(this.proteinComplexs)
                value = value + ...
                    sum(permute(reshape([this.proteinComplexs.rnaComposition],this.knowledgeBase.numGenes, [], this.knowledgeBase.numCompartments),[2 1 3]) .* ...
                    repmat(this.proteinComplexCoefficients, [1 this.knowledgeBase.numGenes this.knowledgeBase.numCompartments]),1);
            end
        end

        %leave on composition tree
        function value = get.proteinMonomerComposition(this)
            %retrieve
            if ~isempty(this.proteinMonomerComposition)
                value = this.proteinMonomerComposition;
                return;
            end

            %calculate
            value = this.proteinMonomerComposition_Helper();

            %store
            this.proteinMonomerComposition = value;
        end

        %helper method because MATLAB won't let you call another objects
        %get.proteinMonomerComposition method from within
        %get.proteinMonomerComposition
        function value = proteinMonomerComposition_Helper(this)
            value = zeros(1, this.knowledgeBase.numProteinMonomers, this.knowledgeBase.numCompartments);

            if ~isempty(this.proteinMonomers)
                value(sub2ind(size(value),...
                    ones(1, numel(this.proteinMonomers)),...
                    double([this.proteinMonomers.idx]),...
                    double([this.proteinMonomerCompartments.idx]))) = ...
                    this.proteinMonomerCoefficients;
            end

            if ~isempty(this.proteinComplexs)
                value = value + ...
                    sum(permute(reshape([this.proteinComplexs.proteinMonomerComposition],this.knowledgeBase.numProteinMonomers, [], this.knowledgeBase.numCompartments),[2 1 3]) .* ...
                    repmat(this.proteinComplexCoefficients, [1 this.knowledgeBase.numProteinMonomers this.knowledgeBase.numCompartments]),1);
            end
        end

        function value = get.molecularWeight(this)
            %retrieve
            if ~isempty(this.molecularWeight)
                value = this.molecularWeight;
                return;
            end

            %calculate
            metaboliteMolecularWeights = [this.metabolites.molecularWeight];

            rnaMolecularWeights = zeros(1, length(this.rnas));
            for i = 1:length(this.rnas)
                rnaMolecularWeights(i) = this.rnas(i).matureMolecularWeight;
            end

            monomerMolecularWeights = zeros(1, length(this.proteinMonomers));
            for i = 1:length(this.proteinMonomers)
                monomerMolecularWeights(i) = this.proteinMonomers(i).matureSequenceMolecularWeight;
            end

            complexMolecularWeights = zeros(1, length(this.proteinComplexs));
            for i = 1:length(this.proteinComplexs)
                complexMolecularWeights(i) = this.proteinComplexs(i).matureMolecularWeight;
            end

            molecularWeights = [metaboliteMolecularWeights rnaMolecularWeights monomerMolecularWeights complexMolecularWeights];
            coefficients = [this.metaboliteCoefficients; this.rnaCoefficients; this.proteinMonomerCoefficients; this.proteinComplexCoefficients];

            if ~isempty(coefficients)
                value = molecularWeights * coefficients;
            else
                value = 0;
            end

            %store
            this.molecularWeight = value;
        end

        function value = get.matureMolecularWeight(this)
            %retrieve
            if ~isempty(this.matureMolecularWeight)
                value = this.matureMolecularWeight;
                return;
            end

            %calculate
            value = this.molecularWeight;
            if ~isempty(this.prostheticGroups)
                value = value + [this.prostheticGroups.molecularWeight] * ...
                    max(1, this.prostheticGroupCoefficients');
            end

            %store
            this.matureMolecularWeight = value;
        end

        function value = get.density(this)
            %retrieve
            if ~isempty(this.density)
                value = this.density;
                return;
            end

            %calculate
            value = 1.35;

            %store
            this.density = value;
        end

        function value = get.volume(this)
            %retrieve
            if ~isempty(this.volume)
                value = this.volume;
                return;
            end

            %calculate
            value = this.molecularWeight / this.density;

            %store
            this.volume = value;
        end

        function value = get.pI(this)
            %retrieve
            if ~isempty(this.pI)
                value = this.pI;
                return;
            end

            %calculate
            value = this.meanSubunitPropertyValue('pI', true);

            %store
            this.pI = value;
        end

        function value = get.extinctionCoefficient(this)
            %retrieve
            if ~isempty(this.extinctionCoefficient)
                value = this.extinctionCoefficient;
                return;
            end

            %calculate
            value = this.meanSubunitPropertyValue('extinctionCoefficient', true);

            %store
            this.extinctionCoefficient = value;
        end

        function value = get.absorbanceFactor(this)
            %retrieve
            if ~isempty(this.absorbanceFactor)
                value = this.absorbanceFactor;
                return;
            end

            %calculate
            value = this.extinctionCoefficient/this.molecularWeight;

            %store
            this.absorbanceFactor = value;
        end

        function value = get.negAA(this)
            %retrieve
            if ~isempty(this.negAA)
                value = this.negAA;
                return;
            end

            %calculate
            value = this.meanSubunitPropertyValue('negAA', false);

            %store
            this.negAA = value;
        end

        function value = get.posAA(this)
            %retrieve
            if ~isempty(this.posAA)
                value = this.posAA;
                return;
            end

            %calculate
            value = this.meanSubunitPropertyValue('posAA', false);

            %store
            this.posAA = value;
        end

        %Instability index (II)
        function value = get.instabilityIndex(this)
            %retrieve
            if ~isempty(this.instabilityIndex)
                value = this.instabilityIndex;
                return;
            end

            %calculate
            value = this.meanSubunitPropertyValue('instabilityIndex', true);

            %store
            this.instabilityIndex = value;
        end

        function value = get.stable(this)
            %retrieve
            if ~isempty(this.stable)
                value = this.stable;
                return;
            end

            %calculate
            value = (this.instabilityIndex<40);

            %store
            this.stable = value;
        end

        %Aliphatic index
        function value = get.aliphaticIndex(this)
            %retrieve
            if ~isempty(this.aliphaticIndex)
                value = this.aliphaticIndex;
                return;
            end

            %calculate
            value = this.meanSubunitPropertyValue('aliphaticIndex', true);

            %store
            this.aliphaticIndex = value;
        end

        %GRAVY (Grand Average of Hydropathy)
        function value = get.gravy(this)
            %retrieve
            if ~isempty(this.gravy)
                value = this.gravy;
                return;
            end

            %calculate
            value = this.meanSubunitPropertyValue('gravy', true);

            %store
            this.gravy = value;
        end

        %base count (numbers of amino acids)
        function value = get.baseCount(this)
            %retrieve
            if ~isempty(this.baseCount)
                value = this.baseCount;
                return;
            end

            %calculate
            rnaBaseCounts = zeros(length(this.rnas), this.knowledgeBase.numMetabolites);
            for i = 1:length(this.rnas)
                rnaBaseCounts(i, :) = this.rnas(i).matureBaseCount;
            end

            monomerBaseCounts = zeros(length(this.proteinMonomers), this.knowledgeBase.numMetabolites);
            for i = 1:length(this.proteinMonomers)
                monomerBaseCounts(i, :) = this.proteinMonomers(i).matureSequenceBaseCount;
            end

            complexBaseCounts = zeros(length(this.proteinComplexs), this.knowledgeBase.numMetabolites);
            for i = 1:length(this.proteinComplexs)
                complexBaseCounts(i, :) = this.proteinComplexs(i).matureBaseCount;
            end

            baseCounts = [rnaBaseCounts; monomerBaseCounts; complexBaseCounts];
            coefficients = [this.rnaCoefficients; this.proteinMonomerCoefficients; this.proteinComplexCoefficients];

            if ~isempty(coefficients)
                value = coefficients' * baseCounts;
            else
                value = zeros(1, this.knowledgeBase.numMetabolites);
            end

            if ~isempty(this.metabolites)
                value([this.metabolites.idx]) = ...
                    value([this.metabolites.idx]) + ...
                    this.metaboliteCoefficients';
            end

            %store
            this.baseCount = value;
        end

        function value = get.matureBaseCount(this)
            %retrieve
            if ~isempty(this.matureBaseCount)
                value = this.matureBaseCount;
                return;
            end

            %calculate
            value = this.baseCount;
            if ~isempty(this.prostheticGroups)
                value([this.prostheticGroups.idx]) = ...
                    value([this.prostheticGroups.idx]) + ...
                    max(1, this.prostheticGroupCoefficients);
            end

            %store
            this.matureBaseCount = value;
        end

        %base count (numbers of amino acids)
        function value = get.cumulativeBaseCount(this)
            throw(MException('ProteinComplex:error','property is not defined'));
        end

        %codon count
        function value = get.tRNACount(this)
            %retrieve
            if ~isempty(this.tRNACount)
                value = this.tRNACount;
                return;
            end

            %calculate
            monomerTRNACounts = zeros(length(this.proteinMonomers),this.knowledgeBase.numMetabolites);
            for i = 1:length(this.proteinMonomers)
                switch this.proteinMonomers(i).gene.type
                    case 'mRNA', monomerTRNACounts(i,:) = this.proteinMonomers(i).matureSequenceBaseCount;
                end
            end

            complexTRNACounts = zeros(length(this.proteinComplexs),this.knowledgeBase.numMetabolites);
            for i = 1:length(this.proteinComplexs)
                complexTRNACounts(i,:) = this.proteinComplexs.matureBaseCount;
            end

            TRNACounts = [monomerTRNACounts;complexTRNACounts];
            coefficients = [this.proteinMonomerCoefficients;this.proteinComplexCoefficients];

            if ~isempty(coefficients)
                value = coefficients'*TRNACounts;
            else
                value = zeros(1,this.knowledgeBase.numMetabolites);
            end

            %store
            this.tRNACount = value;
        end

        function value = get.decayReaction(this)
            %retrieve
            if ~isempty(this.decayReaction)
                value = this.decayReaction;
                return;
            end

            %calculate
            value = this.decayReaction_Helper();

            %store
            this.decayReaction = value;
        end

        function value = decayReaction_Helper(this)
            value = zeros(1, this.knowledgeBase.numMetabolites);

            for i = 1:length(this.proteinComplexs)
                value = value + this.proteinComplexCoefficients(i) * this.proteinComplexs(i).decayReaction;
            end

            if ~isempty(this.metabolites)
                value([this.metabolites.idx]) =   ...
                    + value([this.metabolites.idx]) ...
                    + this.metaboliteCoefficients';
            end
        end

        function value = get.matureDecayReaction(this)
            %retrieve
            if ~isempty(this.matureDecayReaction)
                value = this.matureDecayReaction;
                return;
            end

            %calculate
            value = this.decayReaction;
            if ~isempty(this.prostheticGroups)
                value([this.prostheticGroups.idx]) = ...
                    + value([this.prostheticGroups.idx]) ...
                    + max(1, this.prostheticGroupCoefficients);
            end

            %store
            this.matureDecayReaction = value;
        end

        %number of distinct subunits
        function value = get.numDistinctSubunits(this)
            %retrieve
            if ~isempty(this.numDistinctSubunits)
                value = this.numDistinctSubunits;
                return;
            end

            %calculate
            value = length(this.rnas)+length(this.proteinMonomers)+length(this.proteinComplexs);

            %store
            this.numDistinctSubunits = value;
        end

        %number of subunits
        function value = get.numSubunits(this)
            %retrieve
            if ~isempty(this.numSubunits)
                value = this.numSubunits;
                return;
            end

            %calculate
            value = sum([this.proteinMonomerCoefficients;this.proteinComplexCoefficients]);
            if ~isempty(value)
                value = 0;
            end

            %store
            this.numSubunits = value;
        end

        function value = get.numNucleicAcids(this)
            %retrieve
            if ~isempty(this.numNucleicAcids)
                value = this.numNucleicAcids;
                return;
            end

            %calculate
           value = this.computeNumNucleicAcids();

           %store
            this.numNucleicAcids = value;
        end

        function value = computeNumNucleicAcids(this)
            value = 0;

            for i = 1:length(this.rnas)
                value = value + this.rnas(i).sequenceLength * this.rnaCoefficients(i);
            end

            for i = 1:length(this.proteinComplexs)
                value = value + this.proteinComplexs(i).computeNumNucleicAcids() * this.proteinComplexCoefficients(i);
            end
        end

        function value = get.numAminoAcids(this)
            %retrieve
            if ~isempty(this.numAminoAcids)
                value = this.numAminoAcids;
                return;
            end

            %calculate
            value = this.computeNumAminoAcids();

            %store
            this.numAminoAcids = value;
        end

        function value = computeNumAminoAcids(this)
            value = 0;

            for i = 1:length(this.proteinMonomers)
                value = value + this.proteinMonomers(i).matureSequenceLength * this.proteinMonomerCoefficients(i);
            end

            for i = 1:length(this.proteinComplexs)
                value = value + this.proteinComplexs(i).computeNumAminoAcids() * this.proteinComplexCoefficients(i);
            end
        end

        function value = get.numRNASubunits(this)
            %retrieve
            if ~isempty(this.numRNASubunits)
                value = this.numRNASubunits;
                return;
            end

            %calculate
            value = this.computeNumRNASubunits();

            %store
            this.numRNASubunits = value;
        end

        function value = computeNumRNASubunits(this)
            value = sum(this.rnaCoefficients);

            for i = 1:length(this.proteinComplexs)
                value = value + this.proteinComplexs(i).computeNumRNASubunits() * this.proteinComplexCoefficients(i);
            end
        end

        function value = get.numProteinMonomerSubunits(this)
            %retrieve
            if ~isempty(this.numProteinMonomerSubunits)
                value = this.numProteinMonomerSubunits;
                return;
            end

            %calculate
            value = this.computeNumProteinMonomerSubunits();

            %store
            this.numProteinMonomerSubunits = value;
        end

        function value = computeNumProteinMonomerSubunits(this)
            value = sum(this.proteinMonomerCoefficients);

            for i = 1:length(this.proteinComplexs)
                value = value + this.proteinComplexs(i).computeNumProteinMonomerSubunits() * this.proteinComplexCoefficients(i);
            end
        end

        %average/sum property value over subunits
        function value = meanSubunitPropertyValue(this, property, average)
            values = [];

            if ~isempty(this.rnas)
                values = [values this.rnas.(property)];
            end

            if ~isempty(this.proteinMonomers)
                values = [values this.proteinMonomers.(property)];
            end

            if ~isempty(this.proteinComplexs)
                values = [values this.proteinComplexs.(property)];
            end

            coefficients = [this.rnaCoefficients;this.proteinMonomerCoefficients;this.proteinComplexCoefficients];
            if ~isempty(coefficients)
                if average
                    value = values*coefficients/sum(coefficients);
                else
                    value = values*coefficients;
                end
            else
                if average
                    value = NaN;
                else
                    value = 0;
                end
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [mapping, hs, axes_h, fig_h] = scatter_plot(Y,names,options)
%ENDMEMBER_SCATTER_PLOT Summary of this function goes here
%   Detailed explanation goes here
M = length(names);
Xs = [];
colors = distinguishable_colors(M);
dimension_num = 2;

if nargin > 2 && isstruct(options)
    arg_set = fieldnames(options);
    for i = 1:length(arg_set)
        eval([arg_set{i},'=options.',arg_set{i},';']);
    end
end

[mappedX, mapping] = pca(Y,dimension_num);

fig_h = figure;
hs = cell(0);
if dimension_num == 2
    scatter_h = scatter(mappedX(:,1), mappedX(:,2), 10, '.', ...
        'MarkerEdgeColor', [0.6 0.6 0.6]);
elseif dimension_num == 3
    scatter_h = scatter3(mappedX(:,1), mappedX(:,2), mappedX(:,3), 10, '.', ...
        'MarkerEdgeColor', [0.6 0.6 0.6]);
end

hs(end+1) = {scatter_h};
hold on; axis equal;
axes_h = gca;

for i = 1:length(Xs)
    marker_color = colors(i,:);
    mappedX = (Xs{i} - repmat(mapping.mean, size(Xs{i},1), 1)) * mapping.M;
    if dimension_num == 2
        scatter_h = scatter(axes_h, mappedX(:,1), mappedX(:,2), 5, '+', ...
            'MarkerEdgeColor', marker_color);
    elseif dimension_num == 3
        scatter_h = scatter3(axes_h, mappedX(:,1), mappedX(:,2), ...
            mappedX(:,3), 5, '+', 'MarkerEdgeColor', marker_color);
    end
    hs(end+1) = {scatter_h};
end

% show legend
names = {'Pixels',names{:}};

hs_legend = zeros(1,length(hs));
for j = 1:length(hs)
    hs_legend(j) = hs{j}(1);
end
legend(hs_legend, names,'Location','best');

author: zhouyuanzxcv

repo: Hyperspectral


text:
% This program compares the execution time of mgs.m to mgsfact.m
% on matrices with 200 rows and various numbers of columns ranging
% from 50 to 200 in steps of 10.
%
% The results are saved in a plot named timing.eps.
%
% Dianne O'Leary 09/2005

m = 200;
nsizes = 50:10:m;

disp('Comparing the new algorithm to the old.')
disp('Ratios of old execution times to new:')

for j=1:length(nsizes)
  n = nsizes(j)
  C = rand(m,n);
  tic
  [q,r] = mgs(C);
  t1(j) = toc;
  tic
  [q,r] = mgsfact(C);
  t2(j) = toc;
  disp(sprintf('For %d columns, ratio = %f',n,t1(j)/t2(j)))
end

semilogy(nsizes,t1,nsizes,t2)
xlabel('number of columns')
ylabel('time (sec)')
title('Times for matrices with 200 rows')
legend('Original algorithm','Modified algorithm')

print -depsc timing.eps

author: johannesgerer

repo: jburkardt-m


text:
function test_bug1925

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY surface_nesting ft_headmodel_bemcp

[ftver, ftpath] = ft_version;
cd(fullfile(ftpath, 'forward/private')); % this is where the surface_nesting function is located

[pos, tri] = mesh_sphere(162);

bnd10.id  = 10;
bnd10.pos = pos*10;
bnd10.tri = tri;

bnd20.id  = 20;
bnd20.pos = pos*20;
bnd20.tri = tri;

bnd30.id  = 30;
bnd30.pos = pos*30;
bnd30.tri = tri;

bnd40.id  = 40;
bnd40.pos = pos*40;
bnd40.tri = tri;

bnd50.id  = 50;
bnd50.pos = pos*50;
bnd50.tri = tri;

bnd    = bnd10;
bnd(2) = bnd20;
bnd(3) = bnd30;
bnd(4) = bnd40;
bnd(5) = bnd50;
assert(equalorder(surface_nesting(bnd, 'insidefirst'), 1:5));
assert(equalorder(surface_nesting(bnd, 'outsidefirst'), fliplr(1:5)));

bnd    = bnd10;
bnd(3) = bnd20;
bnd(2) = bnd30;
bnd(5) = bnd40;
bnd(4) = bnd50;
assert(equalorder(surface_nesting(bnd, 'insidefirst'), [1 3 2 5 4]));
assert(equalorder(surface_nesting(bnd, 'outsidefirst'), fliplr([1 3 2 5 4])));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to deal with row and column comparisons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function bool = equalorder(a, b)
bool = isequal(a(:), b(:));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [P] = spm_fp_fmin(M)
% optimises the parameters with respect to an equilibrium density
% FORMAT [P] = spm_fp_fmin(M)
%
% M   - model structure with desired density specified by M(1).fq and
%       support specified by M(1).X = spm_ndgrid(x)
%
% P   - optimised parameters
%
%--------------------------------------------------------------------------
% This routine uses EM (spm_nlsi_NG) and the Fokker Planck formulation to
% minimise the difference between the flow and dispersion terms induced by
% the free parameters of the flow (M(1),f).
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_fp_fmin.m 4136 2010-12-09 22:22:28Z guillaume $

% specify function returning the flow-dependent part of dp/dt 'spm_fp_fun'
%--------------------------------------------------------------------------
M      = M(1);
try, M = rmfield(M,'hE'); end
try, M = rmfield(M,'hC'); end
U      = [];
M.IS   = 'spm_fp_fun';

% Dispersion
%--------------------------------------------------------------------------
N     = size(M.X,1);
D     = inv(M.W)/2;
for i = 1:N
    Y.y(i,1) = trace(D*spm_cat(spm_diff(M.fq,M.X(i,:),[1 1])')');
end

% Optimise
%--------------------------------------------------------------------------
P     = spm_nlsi_GN(M,U,Y);

author: spm


text:
function works = test_leap()
%TEST_LEAP Checks whether LEAP is properly installed.
% Usage:
%   test_leap
%   works = test_leap % returns true/false
%
% See also: install_leap

works = true;

% Check if we can import the LEAP package from anywhere
cdCmd = ['cd "' matlabroot '"'];
if ispc(); cdCmd = ['cd /D "' matlabroot '"']; end
[status,msg] = system([cdCmd ' && python -c "import leap"']);

if status ~= 0
    works = false;
end

if nargout == 0
    if works
        printf('Test LEAP successful!')
    else
        disp('Unable to import LEAP python package. Make sure LEAP and its dependencies are installed.')
        disp('Go to the base LEAP directory containing setup.py and run from MATLAB:')
        disp('    !pip install -e .')
        disp('Or try the MATLAB installer:')
        disp('    install_leap')
    end
    clear works
end

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function zbar = eval2dPoly( x, y, coeffs )
% Given the coefficients of a polynomial as returned by fit2dPolySVD,
% calculates the values z for input values (x,y).
% x, y are column vectors specifying the points to be calculated.
% The vectors must be the same length.
% Coeffs is the coefficients array returned by fit2dPolySVD.

[sizexR, sizexC] = size(x);
[sizeyR, sizeyC] = size(y);

if (sizexC ~= 1) || (sizeyC ~= 1)
    fprintf( 'Inputs of eval2dPoly must be column vectors' );
    return;
end

if (sizeyR ~= sizexR)
    fprintf( 'Inputs vectors of eval2dPoly must be the same length' );
    return;
end

numVals = sizexR;

order = 0.5 * (sqrt(8*length(coeffs)+1) - 3);

zbar = zeros(numVals,1);
column = 1;
for xpower = 0:order
    for ypower = 0:(order-xpower)
        zbar = zbar + (coeffs(column) .* x.^xpower .* y.^ypower);
        column = column + 1;
    end
end

author: Sable

repo: mcbench-benchmarks


text:
options = optimset('GradObj','on');
[x,y]=fminunc('fun3',rand(1,2),options)

author: Eurus-Holmes


text:
    L = iread('rocks2-l.png', 'reduce', 2);
R = iread('rocks2-r.png', 'reduce', 2);
stdisp(L, R)
d = istereo(L, R, [40, 90], 3, 'interp');
idisp(d, 'bar');
Z = 3740*0.160 ./ d;

clf;
surf(Z)
shading interp; view(-150, 75)
set(gca, 'ZDir', 'reverse'); set(gca, 'XDir', 'reverse');
colormap(flipud(hot))

anaglyph(L,R)

author: petercorke

repo: machinevision-toolbox-matlab


text:
%  Generate an n-column convolution matirx of
%  An n-column convolution matrix of p is the matrix for the
%  linear transformation L(f) = p*f  for polynomial f of degree n-1.
%  If f is the coefficient vector of f, and C is the convolution matrix,
%  then C*f is the coefficient vector of p*f.
%
%    Syntax:  C = ConvolutionMatrix(p,n)
%
%    Input:   p --- (string or vector) the univariate polynomial
%             n --- (integer) the column dimension of the convolution matrix
%
%   Output:   C --- the convolution matrix
%
%   Example:  >> p = '1+3*x^4-6*x^8';
%             >> C = ConvolutionMatrix(p,3)
%             C =
%             -6     0     0
%              0    -6     0
%              0     0    -6
%              0     0     0
%              3     0     0
%              0     3     0
%              0     0     3
%              0     0     0
%              1     0     0
%              0     1     0
%              0     0     1

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% plot_lognorm_scat2_1d: Plots the log normalized 2nd order scattering
% coefficients for 1d signals.
% Usage
%    plot_lognorm_scat2_1d(S)
% Input
%    S: The log normalized scattering coefficients.
% Output
%    N/A
% Description
%   Plots the second order log normalized scattering coefficients as a
%   function of j2-j1 (where 2^(-j1) is the first frequency, 2^(-j2) is the
%   second frequency) for various values of j1.

function plot_lognorm_scat2_1d(S)

% Settings

% The lower limit of j2-j1
lowlimit = -2;

% Stop plotting j2-j1 at the max-earlystop
earlystop = 2;

% The values of j1 go up to max-earlyend
earlyend = 1;

%% Initialize

% Adjust settings as needed
earlyend = max(earlyend,earlystop);

% Second order coefficients
S2 = S{3};

% Number of frequencies minus one
maxj1 = max(S2.meta.j(1,:));

% Colormap
C = jet(maxj1-earlyend+1);

% Legend
leg = cell(maxj1-earlyend+1,1);

%% Loop through values of j1 and plot in terms of j2-j1

figure,hold on;
for j1=0:(maxj1-earlyend)
    indj1 = find(S2.meta.j(1,:)==j1);
    j2 = S2.meta.j(2,indj1);

    S2j1j2 = S2.signal(indj1);
    S2j1j2 = cat(2,S2j1j2{:});
    S2j1j2 = S2j1j2(1,:);

    [j2,indj2] = sort(j2,'ascend');
    S2j1j2 = S2j1j2(indj2);

    j2minj1 = j2-j1;
    indj2j1 = j2minj1 >= lowlimit;
    indj2j1((end-earlystop+1):end) = false;

    plot(j2minj1(indj2j1),S2j1j2(indj2j1),'Color',C(j1+1,:));
    leg{j1+1} = sprintf('%d',j1);
end
legend(leg);
hold off;

end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function z0 = bicinter(x, y, z, x0, y0)
N = 16;
M = 4;
[xn, yn, zn] = findnearest(x, y, z, x0, y0, N);
F = zeros(N, N);
for p = 1:N
    Fr = zeros(1, N);
    for i = 1:M
    for j = 1:M
        Fr( map_ijk(i, j, M) ) = (xn(p)^(i-1))*(yn(p)^(j-1));
    end
    end
    F(p, :) = Fr;
end
a = F\zn;
z0 = 0;
for i = 1:M
for j = 1:M
    z0 = z0 + a( map_ijk(i, j, M) )*(x0^(i-1))*(y0^(j-1));
end
end

function k = map_ijk (i, j, M)
    k = (i-1)*M + j;


text:
function [x,y,l]=simulatej(a,b,c,d,g,h,lambda,x0,T)
%Simulation of dynamical system
x(1)=x0;
for t=2:T
if rand>lambda
l(t)=0;
x(t)=random('normal',mx(a,x(t-1)),b);
else
l(t)=1;
x(t)=random('normal',mx([a,g],x(t-1)),sqrt((1+h))*b);
end
y(t)=random('normal',my(c,x(t)),d);
end


text:
function checkFieldNames(structure, field_names)
%CHECKFIELDNAMES   Checks structure field names for existance.
%
% DESCRIPTION:
%       checkFieldNames checks the field names of a MATLAB structure
%       against a cell array of allowable field names. If a field name is
%       found within the structure that is not contained within the cell
%       array, an error is thrown with the field name given in the error
%       message. The function is useful for checking for spelling errors in
%       field names for structure inputs.
%
% USAGE:
%       checkFieldNames(structure, field_names)
%
% INPUTS:
%       structure   - MATLAB structure
%       field_names - cell array of allowable field names
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 13th October 2009
%       last update - 13th October 2009
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also fieldnames

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

names = fieldnames(structure);
for names_index = 1:length(names)
    name_found = false;
    for field_index = 1:length(field_names)
        if strcmp(names{names_index}, field_names{field_index})
            name_found = true;
        end
    end
    if ~name_found
    	error([names{names_index} ' is not a valid field for the structure ' inputname(1)]);
    end
end

repo: Aero-matlab


text:
function trx = genTrack(nframes,vsz)

trx = struct('x',repmat(vsz(2)/2,1,nframes),...
  'y',repmat(vsz(1)/2,1,nframes),...
  'theta',repmat(0,1,nframes),...
  'a',repmat(5,1,nframes),...
  'b',repmat(2,1,nframes),...
  'x_mm',repmat(vsz(2)/2,1,nframes),...
  'y_mm',repmat(vsz(1)/2,1,nframes),...
  'theta_mm',repmat(0,1,nframes),...
  'a_mm',repmat(5,1,nframes),...
  'b_mm',repmat(2,1,nframes),...
  'id',0,...
  'moviename','',....
  'annname','',...
  'off',0,...
  'firstframe',1,...
  'nframes',nframes,...
  'endframe',nframes,...
  'timestamps',(1:nframes)/30,...
  'matname',-1,...
  'dt',repmat(1/30,1,nframes),...
  'fps',30,...
  'pxpermm',1);

author: kristinbranson

repo: JAABA


text:
function tagS = general_equipment_module_tags
%"general_equipment_module_tags"
%   Return the tags used to represent a general equipment module as
%   specified by C.7.5.1 in PS3.3 of 2006 DICOM specification.
%
%   Tags are returned in a struct array with 3 fields:
%
%   Tag: String containing hex DICOM tag of a field.
%  Type: String describing type of field, with 5 options:
%         '1' Field must exist, data must exist and be valid.
%         '2' Field must exist, data can exist and be valid, or be NULL.
%         '3' Field is optional, if the field exists data can exist and be
%             valid or be NULL.
%         '1C' Field must exist under certain conditions and contain valid
%             data.
%         '2C' Field must exist under certain conditions and can contain
%             valid data or be NULL.
%Children: For sequences, a tagS with the same format as this struct array.
%
%JRA 06/06/06
%
%Usage:
%   tagS = general_equipment_module_tags
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Initialize the tagS structure.
tagS = struct('tag', {}, 'tagdec', {}, 'type', {}, 'children', {});

%Create an empty tagS template for sequence creation.
%template = tagS;

%Add tags based on PS3.3 attribute lists.

%Manufacturer
tagS(end+1) = struct('tag', '00080070', 'tagdec', 524400, 'type', '2', 'children', []);

%Institution Name
tagS(end+1) = struct('tag', '00080080', 'tagdec', 524416, 'type', '2', 'children', []);

%Institution Address
tagS(end+1) = struct('tag', '00080081', 'tagdec', 524417, 'type', '3', 'children', []);

%Station Name
tagS(end+1) = struct('tag', '00081010', 'tagdec', 528400, 'type', '3', 'children', []);

%Institutional Department Name
tagS(end+1) = struct('tag', '00081040', 'tagdec', 528448, 'type', '3', 'children', []);

%Manufacturer's Model Name
tagS(end+1) = struct('tag', '00081090', 'tagdec', 528528, 'type', '3', 'children', []);

%Device Serial Number
tagS(end+1) = struct('tag', '00181000', 'tagdec', 1576960, 'type', '3', 'children', []);

%Software Versions
tagS(end+1) = struct('tag', '00181020', 'tagdec', 1576992, 'type', '3', 'children', []);

%Spatial Resolution
tagS(end+1) = struct('tag', '00181050', 'tagdec', 1577040, 'type', '3', 'children', []);

%Date of Last Calibration
tagS(end+1) = struct('tag', '00181200', 'tagdec', 1577472, 'type', '3', 'children', []);

%Time of Last Calibration
tagS(end+1) = struct('tag', '00181201', 'tagdec', 1577473, 'type', '3', 'children', []);

%Pixel Padding Value
tagS(end+1) = struct('tag', '00280120', 'tagdec', 2621728, 'type', '3', 'children', []);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function m_utmgrid(varargin)
% M_UTMGRID Draws a UTM grid on a map.
%        M_UTMGRID('parameter','value',...) with any number (or no)
%        optional parameters is used to draw a UTM grid (i.e. with
%        grid lines at 1 km intervals) for a map using the UTM
%        coordinate system. This can be used instead of, or in addition
%        to, the latitude/longitude grid drawn by M_GRID.
%
%        The optional parameters allow the user
%        to control the look of the grid. These parameters are listed
%        by M_UTMGRID('get'), with default parameters in M_UTMGRID('set');
%
%        This is probably useful only for maps less than a few km across.
%
%        see also M_PROJ, M_GRID

% Rich Pawlowicz (rich@eoas.ubc.ca) 25/April/2018

% These structures are initialized by m_proj()

global MAP_PROJECTION MAP_VAR_LIST

% Have to have initialized a map first

if isempty(MAP_PROJECTION)
  disp('No Map Projection initialized - call M_PROJ(''UTM'',... first!');
  return;
elseif ~strcmp(MAP_PROJECTION.name,'UTM')
  disp('Not a UTM projection');
  return;
elseif strcmp(MAP_VAR_LIST.ellipsoid,'normal')
  disp('Cannot show UTM coordinates using the ``normal'' (earth radius=1) ellipsoid in the m_proj call');
  return;
end

gcolor=get(gca,'gridcolor');
gxcolor=get(gca,'xcolor');
gycolor=get(gca,'ycolor');
glinestyle=get(gca,'gridlinestyle');
glinewidth=get(gca,'linewidth');
gfontsize=get(gca,'fontsize');
gfontname=get(gca,'fontname');
gxaxisloc=get(gca,'xaxislocation');
gyaxisloc=get(gca,'yaxislocation');
gtickdir=get(gca,'tickdir');
gticklen=get(gca,'ticklength'); gticklen=gticklen(1);

k=1;
while k<=length(varargin)
  switch lower(varargin{k}(1:3))
    case 'gri'
      gcolor=varargin{k+1};
    case 'lin'
      switch lower(varargin{k}(1:5))
         case 'linew'
           glinewidth=varargin{k+1};
         case 'lines'
           glinestyle=varargin{k+1};
      end
    case 'fon'
       switch lower(varargin{k}(1:5))
         case 'fonts'
           gfontsize=varargin{k+1};
         case 'fontn'
           gfontname=varargin{k+1};
       end
    case 'xco'
      gxcolor=varargin{k+1};
    case 'yco'
      gycolor=varargin{k+1};
    case 'xax'
      gxaxisloc=varargin{k+1};
    case 'yax'
      gyaxisloc=varargin{k+1};
    case 'tic'
      switch lower(varargin{k}(1:5))
        case 'tickl'
           gticklen=varargin{k+1};
        case 'tickd'
           gtickdir=varargin{k+1};
      end
    case {'get','usa'}

      disp('      ''ticklength'',value');
      disp('      ''tickdir'',( ''in'' | ''out'' )');
      disp('      ''gridcolor'',colorspec');
      disp('      ''linewidth'', value');
      disp('      ''linestyle'', ( linespec | ''none'' )');
      disp('      ''fontsize'',value');
      disp('      ''fontname'',name');
      disp('      ''Xcolor'',colorspec');
      disp('      ''Ycolor'',colorspec');
      disp('      ''XaxisLocation'',( ''bottom'' | ''top'' ) ');
      disp('      ''YaxisLocation'',( ''left'' | ''right'' ) ');
      return;
    case 'set'
      disp(['      ticklength = ' num2str(gticklen)]);
      disp(['      tickdir = ' gtickdir]);
      disp(['      gridcolor = ' gcolor]);
      disp(['      linewidth = ' num2str(glinewidth)]);
      disp(['      linestyle = ' glinestyle]);
      disp(['      fontsize = ' num2str(gfontsize)]);
      disp(['      fontname = ' gfontname]);
      disp(['      Xcolor = ' gxcolor]);
      disp(['      Ycolor = ' gycolor]);
      disp(['      XaxisLocation = ' gxaxisloc]);
      disp(['      YaxisLocation = ' gyaxisloc]);
      return;
  end
  k=k+2;
end

% Bring the grid back!

set(gca,'visible','on','layer','top','tickdir',gtickdir,...
   'gridcolor',gcolor,'gridlinestyle',glinestyle,'linewidth',glinewidth,...
   'xaxislocation',gxaxisloc,'yaxislocation',gyaxisloc,'xcolor',gxcolor',...
   'ycolor',gycolor,'fontname',gfontname,'ticklength',[gticklen .025],...
   'fontsize',gfontsize,'xgrid','on','ygrid','on');

fnt=gfontsize;
fnt2=fnt*0.8;

% Do the x axis

tk=get(gca,'xtick');
% Nothing smaller than 1 km squares
tk(rem(tk,1e3)~=0)=[];

tk1=fix(tk/1e5);
tk2=fix((tk-1e5*tk1)/1e3);
tk3=rem(tk,1e3);

for k=1:length(tk)-1,
    lab{k}=sprintf('\\fontsize{%d}%02d',fnt,tk2(k));
end
k=length(tk);
lab{k}=sprintf('{\\fontsize{%d}%2d}\\fontsize{%d}%02d{\\fontsize{%d}%03d} E',fnt2,tk1(k),fnt,tk2(k),fnt2,tk3(k));
set(gca,'xtick',tk,'xticklabel',lab);

% Do the y axis

tk=get(gca,'ytick');
% Nothing smaller than 1 km squares
tk(rem(tk,1e3)~=0)=[];

tk1=fix(tk/1e5);
tk2=fix((tk-1e5*tk1)/1e3);
tk3=rem(tk,1e3);

for k=1:length(tk)-1,
    lab{k}=sprintf('\\fontsize{%d}%02d',fnt,tk2(k));
end
k=length(tk);
lab{k}=sprintf('{\\fontsize{%d}%2d}\\fontsize{%d}%02d{\\fontsize{%d}%03d} N',fnt2,tk1(k),fnt,tk2(k),fnt2,tk3(k));
set(gca,'ytick',tk,'yticklabel',lab','yticklabelrotation',90);

set(gca,'tag','m_utmgrid');

author: goGPS-Project


text:
function subpak_test11 ( )

%*****************************************************************************80
%
%% TEST11 tests GRID1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 5;
  nstep = 11;

  x1 = [ 1.0,  0.0, 20.0, -5.0, 1.0 ];
  x2 = [ 1.0, 10.0,  0.0,  5.0, 2.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST11\n' );
  fprintf ( 1, '  GRID1 computes a 1D grid between\n' );
  fprintf ( 1, '  two DIM_NUM dimensional points X1 and X2.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here, we will use %d steps\n', nstep );
  fprintf ( 1, '  going from\n' );
  for i = 1 : dim_num
    fprintf ( 1, '  %12f', x1(i) );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  to\n' );
  for i = 1 : dim_num
    fprintf ( 1, '  %12f', x2(i) );
  end
  fprintf ( 1, '\n' );

  x = grid1 ( dim_num, nstep, x1, x2 );

  r8mat_transpose_print ( dim_num, nstep, x, '  The grid matrix:' );

  return
end


text:
classdef SettingsTranslator < handle

    properties (GetAccess = public, SetAccess = private)
        fileName
        filePath
    end

    properties (Access = private)
        oldSettings
        propList
        nProps
    end

    methods (Access = public)

        function translate(obj,oldSettings)
            obj.init(oldSettings);
            for i = 1:obj.nProps
                prop = obj.propList{i};
                value = oldSettings.(prop);
                if strcmp(prop,'filename')
                    s.problemData.femData.fileName = value;
                elseif strcmp(prop,'homegenizedVariablesComputer')
                    s.homogenizedVarComputerSettings.type = value;
                    s = obj.translateHomogenizedVariablesComputer(s);
                elseif strcmp(prop,'designVariable')
                    s.designVarSettings.type = value;
                elseif strcmp(prop,'initial_case')
                    s.designVarSettings.initialCase = value;
                elseif strcmp(prop,'m1')
                    if ~isempty(value)
                        s.designVarSettings.creatorSettings.m1 = value;
                    end
                elseif strcmp(prop,'m2')
                    if ~isempty(value)
                        s.designVarSettings.creatorSettings.m2 = value;
                    end
                elseif strcmp(prop,'alpha0')
                    if ~isempty(value)
                        s.designVarSettings.creatorSettings.alpha0 = value;
                    end
                elseif strcmp(prop,'rho0')
                    if ~isempty(value)
                        s.designVarSettings.creatorSettings.rho = value;
                    end
                elseif strcmp(prop,'levelSetDataBase')
                    s = obj.translateLevelSetCreator(s);
                elseif strcmp(prop,'nsteps')
                    s.incrementalSchemeSettings.nSteps = value;
                elseif strcmp(prop,'printIncrementalIter')
                    s.incrementalSchemeSettings.shallPrintIncremental = value;
                elseif strcmp(prop,'Vfrac_initial')
                    s.incrementalSchemeSettings.targetParamsSettings.VfracInitial = value;
                elseif strcmp(prop,'Vfrac_final')
                    s.incrementalSchemeSettings.targetParamsSettings.VfracFinal = value;
                elseif strcmp(prop,'optimality_initial')
                    s.incrementalSchemeSettings.targetParamsSettings.optimalityInitial = value;
                elseif strcmp(prop,'optimality_final')
                    s.incrementalSchemeSettings.targetParamsSettings.optimalityFinal = value;
                elseif strcmp(prop,'constr_initial')
                    s.incrementalSchemeSettings.targetParamsSettings.constrInitial = value;
                elseif strcmp(prop,'constr_final')
                    s.incrementalSchemeSettings.targetParamsSettings.constrFinal = value;
                elseif strcmp(prop,'stressNormExponent_initial')
                    s.incrementalSchemeSettings.targetParamsSettings.stressNormExponentInitial = value;
                elseif strcmp(prop,'stressNormExponent_final')
                    s.incrementalSchemeSettings.targetParamsSettings.stressNormExponentFinal = value;
                elseif strcmp(prop,'cost')
                    for k = 1:length(value)
                        s.costSettings.shapeFuncSettings{k}.type = value{k};
                        if strcmp(value{k},'chomog_alphabeta') || strcmp(value{k},'chomog_fraction')
                            if isprop(oldSettings,'micro')
                                s.costSettings.shapeFuncSettings{k}.alpha = oldSettings.micro.alpha;
                                s.costSettings.shapeFuncSettings{k}.beta = oldSettings.micro.beta;
                            end
                        elseif strcmp(value{k},'enforceCh_CCstar_L2')
                            s.costSettings.shapeFuncSettings{k}.ChTarget.type = oldSettings.selectiveC_Cstar;

                        elseif strcmp(value{k},'perimeterConstraint')
                            if isprop(oldSettings,'Perimeter_target')
                                s.costSettings.shapeFuncSettings{k}.perimeterTarget = oldSettings.Perimeter_target;
                            end
                        end
                    end
                elseif strcmp(prop,'weights')
                    s.costSettings.weights = value;
                elseif strcmp(prop,'filter')
                    s.costSettings.filterType = value;
                    s.constraintSettings.filterType = value;
                elseif strcmp(prop,'constraint')
                    for k = 1:length(value)
                        s.constraintSettings.shapeFuncSettings{k}.type = value{k};
                        if strcmp(value{k},'perimeterConstraint')
                            if isprop(oldSettings,'Perimeter_target')
                                s.constraintSettings.shapeFuncSettings{k}.perimeterTarget = oldSettings.Perimeter_target;
                            end
                        end
                    end
                elseif strcmp(prop,'constraintDomainNotOptimizable')
                    for k = 1:numel(value)
                        s.constraintSettings.shapeFuncSettings{k}.domainNotOptimizable = value{k};
                    end
                elseif strcmp(prop,'costDomainNotOptimizable')
                    for k = 1:numel(value)
                        s.costSettings.shapeFuncSettings{k}.domainNotOptimizable = value{k};
                    end
                elseif strcmp(prop,'optimizer')
                    s.optimizerSettings.type = value;
                elseif strcmp(prop,'constraint_case')
                    s.optimizerSettings.constraintCase = value;
                elseif strcmp(prop,'optimizerUnconstrained')
                    s.optimizerSettings.uncOptimizerSettings.type = value;
                elseif strcmp(prop,'e2')
                    s.optimizerSettings.uncOptimizerSettings.e2 = value;
                elseif strcmp(prop,'ub')
                    s.optimizerSettings.uncOptimizerSettings.ub = value;
                elseif strcmp(prop,'lb')
                    s.optimizerSettings.uncOptimizerSettings.lb = value;
                elseif strcmp(prop,'line_search_initiator')
                    s.optimizerSettings.uncOptimizerSettings.lineSearchSettings.lineSearchInitiatorSettings.type = value;
                elseif strcmp(prop,'HJiter0')
                    s.optimizerSettings.uncOptimizerSettings.lineSearchSettings.HJiter0 = value;
                elseif strcmp(prop,'incrementFactor')
                    s.optimizerSettings.uncOptimizerSettings.lineSearchSettings.lineSearchInitiatorSettings.incrementFactor = value;
                elseif strcmp(prop,'rate')
                    s.optimizerSettings.uncOptimizerSettings.lineSearchSettings.rate = value;
                elseif strcmp(prop,'maxiter')
                    s.optimizerSettings.maxIter = value;
                elseif strcmp(prop,'printing')
                    s.optimizerSettings.shallPrint = value;
                elseif strcmp(prop,'printMode')
                    s.optimizerSettings.printMode = value;
                elseif strcmp(prop,'monitoring')
                    s.optimizerSettings.monitoringDockerSettings.showOptParams = value;
                elseif strcmp(prop,'monitoring_interval')
                    s.optimizerSettings.monitoringDockerSettings.refreshInterval = value;
                elseif strcmp(prop,'ptype')
                    s.optimizerSettings.monitoringDockerSettings.scale = value;
                elseif strcmp(prop,'plotting')
                    s.optimizerSettings.monitoringDockerSettings.shallDisplayDesignVar = value;
                elseif strcmp(prop,'showBC')
                    s.optimizerSettings.monitoringDockerSettings.shallShowBoundaryConditions = value;
                elseif strcmp(prop,'isDesignVariableFixed')
                    s.designVarSettings.isFixed = value;
                end
            end
            if strcmp(s.designVarSettings.type,'MicroParams')
                s.designVarSettings.creatorSettings.homogSettings = s.homogenizedVarComputerSettings;
            end

            obj.exportFile(s);
        end

    end

    methods (Access = private)

        function init(obj,oldSettings)
            obj.oldSettings = oldSettings;
            obj.propList = fieldnames(oldSettings);
            obj.nProps = length(obj.propList);
        end

        function exportFile(obj,s)
            obj.getFilePath();
            str = jsonencode(s);
            fid = fopen([obj.filePath,'.json'],'w+');
            fprintf(fid,str);
            fclose(fid);
        end

        function getFileName(obj)
            old = obj.oldSettings.case_file;
            new = old;
            %new = [old '.json'];
            obj.fileName = new;
        end

        function getFilePath(obj)
            obj.getFileName();
            obj.filePath = fullfile('.','tests','Applications',obj.fileName);
        end

        function s = translateLevelSetCreator(obj,s)
            prop = obj.oldSettings.levelSetDataBase;
            if ~isempty(prop)
                fields = fieldnames(prop);
                for i = 1:length(fields)
                    field = fields{i};
                    value = prop.(field);
                    s.designVarSettings.creatorSettings.(field) = value;
                end
            end
        end

        function s = translateHomogenizedVariablesComputer(obj,s)
            old = obj.oldSettings;
            switch s.homogenizedVarComputerSettings.type
                case 'ByInterpolation'
                    if isprop(old,'materialInterpolation')
                        s.homogenizedVarComputerSettings.interpolation = old.materialInterpolation;
                    end
                    if isprop(old,'material')
                        s.homogenizedVarComputerSettings.typeOfMaterial = old.material;
                    end
                    if isprop(old,'TOL')
                        s.homogenizedVarComputerSettings.constitutiveProperties = old.TOL;
                    end
                case 'ByVademecum'
                    if isprop(old,'vademecumFileName')
                        s.homogenizedVarComputerSettings.fileName = old.vademecumFileName;
                    end
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function print_error_red(msg, tab_before)
% show error in red

    if nargin < 1
        disp('usage: print_error_red(msg)');
        return;
    end

    if nargin < 2
        tab_before = 1;
    end

    if ~iscell(msg)
        disp('error message must be a cell array of strings.');
        return;
    end

    % make sure we have a column vector
    if size(msg, 1) > size(msg, 2)
        msg = msg';
    end

    %% create a fake error output
    if tab_before == 1
        for i = 1 : numel(msg)
            fprintf(2, '\t %s \n', msg{1, i});
        end
    else
        for i = 1 : numel(msg)
            fprintf(2, ' %s \n', msg{1, i});
        end
    end

end
%% EOF

author: thomaskuestner


text:
function varargout = vals2coeffs( U, varargin )
%VALS2COEFFS    Convert matrix of values to Chebyshev-Fourier coefficients.
%
% V = VALS2COEFFS( U ) converts a matrix U of values representing
% samples of a function from a tensor Chebyshev-Fourier grid
% to a matrix V of Chebyshev-Fourier coefficients for the corresponding
% interpolant.
%
% [U, S, V] = VALS2COEFFS( U, S, V ) the same as above but keeps
% everything in low rank form.
%
% See also COEFFS2VALS

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

if ( nargin == 1 )
    U = chebtech2.vals2coeffs( U );
    U = trigtech.vals2coeffs( U.' ).';
    varargout = {U};
elseif ( nargin == 3 )
    S = varargin{1};
    V = varargin{2};
    U = chebtech2.vals2coeffs( U );
    V = trigtech.vals2coeffs( V );
    varargout = { U S V };
else
    error('CHEBFUN:DISKFUN:vals2coeffs:inputs', ...
        'The number of input arguments should be one or two.');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function handles = plotDataMatrix(hObject, handles)

plotDimensions = get(handles.lstPlotDimensions, 'Value');

FigureTitle = 'Data';
[handles, handles.figData] = openPlotFigure(hObject, handles, ...
    'Data Plot (Matrix)', FigureTitle);

% OLD --
%plotmatrix((handles.Data(plotDimensions, :))', ...
%    (handles.Data(plotDimensions, :))');
% -- OLD

plotmatrix((handles.Data(plotDimensions, :))');

% Avoid bug
% (plotmatrix takes away the title created by openPlotFigure)
title(FigureTitle)

hold off;

guidata(hObject, handles);

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
%FIG_DESCI_CASSI_BIRD Demonstrate decompress snapshot compressive imaging
%(DeSCI) results of simulated coded aperture snapshot spectral imaging
%(CASSI) `bird` dataset.
% Note: Run this demonstration script after completing the reconstruction
%       process TEST_DESCI_CASSI_BIRD and obtaining the saved .mat file.
%   See also TEST_DESCI_CASSI_BIRD, GAPDENOISE_CACTI, GAPDENOISE.
clear; clc;
% close all
addpath('../'); % repository root
addpath(genpath('../packages/')); % packages
% [1] run the corresponding test file for recovery or load the saved
% results
% test_desci_cassi_bird; % run the corresponding test file for recovery
load('../results/savedmat/desci_cassi_bird.mat'); % load the saved results

methods   = {'gaptv','desci'}; % all methods for comparison
methnames = {'GAP-TV','DeSCI'}; % corresponding names of the methods
nim = nframe*nmask; % number of images

%% [2.1] demonstrate the reconstructed spectum
colors = {'Yellow bird', 'Green bird', 'Blue bird', 'Purple bird'}; % colors of the four birds
crops  =  [ 79 413 60 60;  % yellow [x y width height] as insertShape
           356 500 60 60;  % green
           780 460 60 60;  % blue
           950 600 60 60]; % purple
locos = {'southeast','southeast','southeast','northwest'};

markers   = {'bx--','r*-'}; % markers

linewidth = 2;
markersize = 4;
textfontsize = 14;
legendfontsize = 12;

% Change default axes fonts.
set(0,'DefaultAxesFontName','Arial');
set(0,'DefaultAxesFontSize',textfontsize);
set(0,'DefaultAxesFontWeight','normal');
% Change default text fonts.
set(0,'DefaultAxesFontName','Arial');
set(0,'DefaultTextFontSize',textfontsize);
set(0,'DefaultTextFontWeight','normal');

% legends
lgall = 'lgtruth';
lgtruth = 'Ground truth';
for imeth = 1:length(methods)
    meth = methods{imeth}; % method
    eval(sprintf('lg%s=''%s'';',meth,methnames{imeth}));
    lgall = [lgall ',lg' meth];
end

f = figure('position',[25 50 800 900]); % PSNR plot
h = tight_subplot(3,2,[0.07 0.08],[.07 .035],[.1 .02]);
axes(h(1));
imshow(insertShape(orig_rgb,'Rectangle',crops,'Color','cyan','Linewidth',8));
title('Original');
axes(h(2));
imshow(meas,[]);
title('Coded frame');
for ic = 1:length(colors)
    color = colors{ic};
    crop = crops(ic,:);
    sp_truth = squeeze(sum(sum(orig(crop(2):crop(2)+crop(4)-1,crop(1):crop(1)+crop(3)-1,:)/MAXB,2),1));
    axes(h(ic+2));
    plot(wavelength,sp_truth/max(sp_truth),'k-','linewidth',linewidth,'markersize',markersize); hold on;

    for imeth = 1:length(methods)
        meth = methods{imeth};
        eval(sprintf('sp_%s = squeeze(sum(sum(v%s(crop(2):crop(2)+crop(4)-1,crop(1):crop(1)+crop(3)-1,:),2),1));',meth,meth));
        eval(sprintf('plot(wavelength,rescale(sp_%s,min(sp_truth)/max(sp_truth),1),markers{imeth},''linewidth'',linewidth,''markersize'',markersize);',meth));
        eval(sprintf('corrall(ic,imeth)=corr(sp_%s,sp_truth);',meth));
        % legend with correlation
        eval(sprintf('lg%s=''%s, corr: %.5f'';',meth,methnames{imeth},corrall(ic,imeth)));
    end
    xlim([390 700]); xticks(400:100:700);
    title(color);
    if ic==3 || ic==4
        xlabel('Wavelength (nm)');
    end
    if ic==1 || ic==3
        ylabel('Intensity (a.u.)');
    end
    eval(sprintf('hlg=legend(%s);',lgall));
    set(hlg,'FontSize',legendfontsize);
    legend('location',locos{ic});
    legend('boxoff');
end

%% [2.2] show exemplar frames (wavelength)
bands = [9 11 16 21];

nb = length(bands);
nmeth = length(methods);
f = figure('Position',[50 50 1010 950]);
h = tight_subplot(nb,nmeth+1,[.002 .002],[.03 .03],[.02 .02]);
set(f, 'Color', 'white');
for ib = 1:nb
    band = bands(ib);
    lambda = wavelength(band);
    wlmap = (gray*kron(ones(3,1),spectrumRGB(lambda))); % colormap with the corresponding RGB wavelength
    wlmap = wlmap/max(wlmap(:));

    axes(h(1+(ib-1)*(nmeth+1)));
    imshow(orig(:,:,band)/MAXB,'colormap',wlmap);
    if ib==1
        title('Ground truth');
    end
    for imeth = 1:nmeth
        meth = methods{imeth};
        axes(h(imeth+1+(ib-1)*(nmeth+1)));
        eval(sprintf('imshow(v%s(:,:,band),''colormap'',wlmap);',meth));
        if ib==1
            title(sprintf('%s',methnames{imeth}));
        end
    end
end

%% [2.3] show all the frames wavelength-by-wavelength
%         via vshowSpectralData
addpath(genpath('../packages/'));

opts = [];
  opts.npcol      = 6; % number of columns for subplots
  % opts.magsize    = 0.7; % magnification of the image size
  opts.width      = 1080; % width of the plot
  opts.textpos    = [600 600]; % position of the text on the sub-images
  opts.fontname   = 'Myriad Pro'; % font name of the text on sub-images
  opts.fontsize   = 15; % fone size of the text on sub-images
  opts.fontweight = 'bold'; % font weight of the text on sub-images
% show ground truth
vshowSpectralData(orig/MAXB,wavelength,opts);
for imeth = 1:length(methods)
    meth = methods{imeth};
    eval(sprintf('vshowSpectralData(v%s,wavelength,opts);',meth));
end

repo: DeSCI


text:
function len = strlen(str)
%  len = strlen(str)
% compute the # of characters in str (ignoring 0s at the end)

%
% strlen.m
%
% Original Author: Bruce Fischl
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:13 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

len = length(str) ;
for i=length(str):-1:1
	if (str(i) ~= 0)
		break ;
	end

	len = len-1;
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MRiLab auto generated file: DO NOT EDIT!     %
% Generated by MRiLab "DoWriteXML2m" Generator %
% MRiLab Version 1.3                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dB0=Mag_GaussianHead
%====================================
p.DeltaX=0.08;
p.DeltaY=0.08;
p.DeltaZ=0.08;
p.PosX=0;
p.PosY=0;
p.PosZ=0;
p.Scale=5e-6;
dB0t=MagGaussian(p);
dB0=dB0t;
p=[];
%--------------------
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
%%

cs = crystalSymmetry('m-3m');
ss = specimenSymmetry('1');

odf = unimodalODF(quaternion.id,cs,ss,'halfwidth',1.5*degree);

%%

[~,q] = max(odf)

%%
v = [];
r = linspace(1*degree,10*degree,10);
for i = 1:length(r)
  fprintf('.');
  v(i,1) = volume(odf,q,r(i));
  v(i,2) = volume(odf,quaternion.id,r(i));
  %v(i,1) = volume(uniformODF(cs,ss),quaternion.id,r(i));
end
fprintf('\n');
%v(:,2) = length(cs)*(r - sin(r))./pi;

plot(r/degree,v);

%%
clear v;
cs = crystalSymmetry('m-3m');
ss = specimenSymmetry('mmm');
r = linspace(0*degree,60*degree,20);
omega = linspace(0,45*degree,4);
for i = 1:length(r)
  for j = 1:length(omega)
    %disp(r(i));
    %tic
    v(i,j,1) = volume(uniformODF(cs,ss),axis2quat(vector3d(1,1,1),omega(j)),r(i));
    %toc
    %tic
    %v(i,j,2) = volume(uniformODF(cs,ss),axis2quat(xvector,omega(j)),r(i),'local');
    %toc
    fprintf('.');
  end
end
fprintf('\n');
plot(r/degree,reshape(v,size(v,1),[]));

%%

%r = linspace(0,100*degree,20);
%for i = 1:length(r)
%  v(i,1) = volume(odf,calcModes(odf),r(i));
%  v(i,2) = volume(uniformODF(cs,ss),quaternion.id,r(i));
%end

%plot(r/degree,v);

%%
cs = crystalSymmetry('m-3m');
ss = specimenSymmetry('mmm');
S3G = SO3Grid(1*degree,cs,ss);
r = plotS2Grid('resolution',10*degree,'hemisphere','upper','maxrho',90*degree,'RESTRICT2MINMAX');
rv = vector3d(r);
d1 = zeros(size(rv));
d2 = zeros(size(rv));
for i = 1:length(rv)
  q = axis2quat(rv(i),20*degree);
  d1(i) = length(find(dot_outer(S3G,q,'epsilon',20*degree)));
  d2(i) = volume(uniformODF(cs,ss),q,15*degree);
  fprintf('.');
end
fprintf('\n');
figure(1)
plot(r,d1./d1(1),'smooth')
figure(2)
plot(r,d2./d2(1),'smooth')
colorbar
%%

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function pass = test_norm( )
% Test spherefun norm() command

tol = 1000*chebfunpref().cheb2Prefs.chebfun2eps;

f = spherefun( @(x,y,z) 1+0*x );
pass(1) = abs( sqrt(sum2( f )) - norm( f ) ) < tol;

f = spherefun( @(x,y,z) cos(x.*y.*z) );
s = svd( f );
pass(2) = abs( sum(s.^2) - norm(f).^2 ) < tol;

f = spherefun( @(x,y,z) x + y + z );
pass(3) = abs( norm(f,inf) - sqrt(3) ) < tol;

end


text:
function tfocs_ver2 = tfocs_version

%TFOCS_VERSION   Version information.
%    Prints information about the version of TFOCS being used, and the
%    system on which it is running. When submitting a bug report, please
%    run this function and include its output in your report.

tfocs_ver = '1.4';
if nargout == 0,
    fprintf( 'TFOCS v%s\n', tfocs_ver );
    verd = ver('MATLAB');
    fprintf( 'MATLAB version %s %s on %s\n', verd.Version, verd.Release, computer );
else
    tfocs_ver2 = tfocs_ver;
end

% TFOCS v1.4 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2015 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function y = scale_cols(x, s)
% SCALE_COLS       Scale each column of a matrix.
% SCALE_COLS(x,s) returns matrix y, same size as x, such that
% y(:,i) = x(:,i)*s(i)
% It is more efficient than x*diag(s), but consumes a similar amount of memory.
% Warning: It consumes a lot of memory when x is sparse.

y = x.*repmat(s(:).', rows(x), 1);
%y = x.*(ones(rows(x),1)*s(:)');

author: andrewssobral

repo: lrslibrary


text:
function params = tSeriesParamsDefault;
% Default signal-processing parameters for time series fMRI data.
%
% params = tSeriesParamsDefault;
%
% This is a rough adaptation of the mrVista 1.0 params. But,
% I expect this to be updated significantly. It should, add least,
% include more explicit specification of desired high- and
% low-pass frequency filtering cutoffs. (Currently this is only
% achieved by setting detrend = 1 and setting the detrend frames
% to 1/[desired high-pass cutoff] for an approximate high-pass detrend
% using multiple boxcar smoothing.)
%
% Conversely, if the preprocessing done by mrDetrend is satisfactory,
% then we may want to ditch the whole idea of having separate tSeries
% methods, and assume all mr objects w/ time series data have already
% been pre-processed.
%
% ras, 08/05

% Detrend options: integer w/ following values
% 0: Don't Detrend
% 1: Detrend w/ multiple boxcar smoothing
% 2: Linear Trend Removal
% 3: Quartic Trend Removal
params.detrend = 1;

% frames to use as period for boxcar detrending
params.detrendFrames = 20;

% Options for how to compensate for distance from the coil, depending
% on the value of inhomoCorrection
%   0 do nothing
%   1 divide by the mean, independently at each voxel
%   2 divide by null condition
%   3 divide by anything you like, e.g., robust estimate of intensity
%   inhomogeneity
params.inhomoCorrect = 1;

% temporal normalization:
% if 1, normalize each frame to have same mean intensity
params.temporalNormalization = 0;

% subtract the mean:
% if 1, subtract the mean, setting the overall mean to 0,
% if 0, don't do this
params.subtractMean = 1;

return


text:
%% This function will generate the PDE function based on calculation result of the symbolic expression.
% This file has been modified to perform the parallel computing.
% Last Update: 2019/06/20
% Coded By: K

function Generate_PDE_RHS(Sindy_PDEs,Vars,name)
f_x=Sindy_PDEs;
matlabFunction(f_x,'File',name,'Optimize',true,'Vars',{Vars});

author: dynamicslab

repo: SINDy-PI