ans =
Latest Contributions
This stems purely from some play on my part. Suppose I asked you to work with the sequence formed as 2*n*F_n + 1, where F_n is the n'th Fibonacci number? Part of me would not be surprised to find there is nothing simple we could do. But, then it costs nothing to try, to see where MATLAB can take me in an explorative sense.
n = sym(0:100).';
Fn = fibonacci(n);
Sn = 2*n.*Fn + 1;
Sn(1:10) % A few elements
For kicks, I tried asking ChatGPT. Giving it nothing more than the first 20 members of thse sequence as integers, it decided this is a Perrin sequence, and gave me a recurrence relation, but one that is in fact incorrect. Good effort from the Ai, but a fail in the end.
Is there anything I can do? Try null! (Look carefully at the array generated by Toeplitz. It is at least a pretty way to generate the matrix I needed.)
X = toeplitz(Sn,[1,zeros(1,4)]);
rank(X(5:end,:))
Hmm. So there is no linear combination of those columns that yields all zeros, since the resulting matrix was full rank.
X = toeplitz(Sn,[1,zeros(1,5)]);
rank(X(6:end,:))
But if I take it one step further, we see the above matrix is now rank deficient. What does that tell me? It says there is some simple linear combination of the columns of X(6:end,:) that always yields zero. The previous test tells me there is no shorter constant coefficient recurrence releation, using fewer terms.
null(X(6:end,:))
Let me explain what those coefficients tell me. In fact, they yield a very nice recurrence relation for the sequence S_n, not unlike the original Fibonacci sequence it was based upon.
S(n+1) = 3*S(n) - S_(n-1) - 3*S(n-2) + S(n-3) + S(n-4)
where the first 5 members of that sequence are given as [1 3 5 13 25]. So a 6 term linear constant coefficient recurrence relation. If it reminds you of the generating relation for the Fibonacci sequence, that is good, because it should. (Remember I started the sequence at n==0, IF you decide to test it out.) We can test it out, like this:
SfunM = memoize(@(N) Sfun(N));
SfunM(25)
2*25*fibonacci(sym(25)) + 1
And indeed, it works as expected.
function Sn = Sfun(n)
switch n
case 0
Sn = 1;
case 1
Sn = 3;
case 2
Sn = 5;
case 3
Sn = 13;
case 4
Sn = 25;
otherwise
Sn = Sfun(n-5) + Sfun(n-4) - 3*Sfun(n-3) - Sfun(n-2) +3*Sfun(n-1);
end
end
A beauty of this, is I started from nothing but a sequence of integers, derived from an expression where I had no rational expectation of finding a formula, and out drops something pretty. I might call this explorational mathematics.
The next step of course is to go in the other direction. That is, given the derived recurrence relation, if I substitute the formula for S_n in terms of the Fibonacci numbers, can I prove it is valid in general? (Yes.) After all, without some proof, it may fail for n larger than 100. (I'm not sure how much I can cram into a single discussion, so I'll stop at this point for now. If I see interest in the ideas here, I can proceed further. For example, what was I doing with that sequence in the first place? And of course, can I prove the relation is valid? Can I do so using MATLAB?)
(I'll be honest, starting from scratch, I'm not sure it would have been obvious to find that relation, so null was hugely useful here.)
Function Syntax Design Conundrum
As a MATLAB enthusiast, I particularly enjoy Steve Eddins' blog and the cool things he explores. MATLAB's new argument blocks are great, but there's one frustrating limitation that Steve outlined beautifully in his blog post "Function Syntax Design Conundrum": cases where an argument should accept both enumerated values AND other data types.
Steve points out this could be done using the input parser, but I prefer having tab completions and I'm not a fan of maintaining function signature JSON files for all my functions.
Personal Context on Enumerations
To be clear: I honestly don't like enumerations in any way, shape, or form. One reason is how awkward they are. I've long suspected they're simply predefined constructor calls with a set argument, and I think that's all but confirmed here. This explains why I've had to fight the enumeration system when trying to take arguments of many types and normalize them to enumerated members, or have numeric values displayed as enumerated members without being recast to the superclass every operation.
The Discovery
While playing around extensively with metadata for another project, I realized (and I'm entirely unsure why it took so long) that the properties of a metaclass object are just, in many cases, the attributes of the classdef. In this realization, I found a solution to Steve's and my problem.
To be clear: I'm not in love with this solution. I would much prefer a better approach for allowing variable sets of membership validation for arguments. But as it stands, we don't have that, so here's an interesting, if incredibly hacky, solution.
If you call struct() on a metaclass object to view its hidden properties, you'll notice that in addition to the public "Enumeration" property, there's a hidden "Enumerable" property. They're both logicals, which implies they're likely functionally distinct. I was curious about that distinction and hoped to find some functionality by intentionally manipulating these values - and I did, solving the exact problem Steve mentions.
The Problem Statement
We have a function with an argument that should allow "dual" input types: enumerated values (Steve's example uses days of the week, mine uses the "all" option available in various dimension-operating functions) AND integers. We want tab completion for the enumerated values while still accepting the numeric inputs.
A Solution for Tab-Completion Supported Arguments
Rather than spoil Steve's blog post, let me use my own example: implementing a none() function. The definition is simple enough tf = ~any(A, dim); but when we wrap this in another function, we lose the tab-completion that any() provides for the dim argument (which gives you "all"). There's no great way to implement this as a function author currently - at least, that's well documented.
So here's my solution:
%% Example Function Implementation
% This is a simple implementation of the DimensionArgument class for implementing dual type inputs that allow enumerated tab-completion.
function tf = none(A, dim)
arguments(Input)
A logical;
dim DimensionArgument = DimensionArgument(A, true);
end
% Simple example (notice the use of uplus to unwrap the hidden property)
tf = ~any(A, +dim);
end
I like this approach because the additional work required to implement it, once the enumeration class is initialized, is minimal. Here are examples of function calls, note that the behavior parallels that of the MATLAB native-style tab-completion:
%% Test Data
% Simple logical array for testing
A = randi([0, 1], [3, 5], "logical");
%% Example function calls
tf = none(A, "all"); % This is the tab-completion it's 1:1 with MATLABs behavior
tf = none(A, [1, 2]); % We can still use valid arguments (validated in the constructor)
tf = none(A); % Showcase of the constructors use as a default argument generator
How It Works
What makes this work is the previously mentioned Enumeration attribute. By setting Enumeration = false while still declaring an enumeration block in the classdef file, we get the suggested members as auto-complete suggestions. As I hinted at, the value of enumerations (if you don't subclass a builtin and define values with the someMember (1) syntax) are simply arguments to constructor calls.
We also get full control over the storage and handling of the class, which means we lose the implicit storage that enumerations normally provide and are responsible for doing so ourselves - but I much prefer this. We can implement internal validation logic to ensure values that aren't in the enumerated set still comply with our constraints, and store the input (whether the enumerated member or alternative type) in an internal property.
As seen in the example class below, this maintains a convenient interface for both the function caller and author the only particuarly verbose portion is the conversion methods... Which if your willing to double down on the uplus unwrapping context can be avoided. What I have personally done is overload the uplus function to return the input (or perform the identity property) this allowss for the uplus to be used universally to unwrap inputs and for those that cant, and dont have a uplus definition, the value itself is just returned:
classdef(Enumeration = false) DimensionArgument % < matlab.mixin.internal.MatrixDisplay
%DimensionArgument Enumeration class to provide auto-complete on functions needing the dimension type seen in all()
% Enumerations are just macros to make constructor calls with a known set of arguments. Declaring the 'all'
% enumeration member means this class can be set as the type for an input and the auto-completion for the given
% argument will show the enumeration members, allowing tab-completion. Declaring the Enumeration attribute of
% the class as false gives us control over the constructor and internal implementation. As such we can use it
% to validate the numeric inputs, in the event the 'all' option was not used, and return an object that will
% then work in place of valid dimension argument options.
%% Enumeration members
% These are the auto-complete options you'd like to make available for the function signature for a given
% argument.
enumeration(Description="Enumerated value for the dimension argument.")
all
end
%% Properties
% The internal property allows the constructor's input to be stored; this ensures that the value is store and
% that the output of the constructor has the class type so that the validation passes.
% (Constructors must return the an object of the class they're a constructor for)
properties(Hidden, Description="Storage of the constructor input for later use.")
Data = [];
end
%% Constructor method
% By the magic of declaring (Enumeration = false) in our class def arguments we get full control over the
% constructor implementation.
%
% The second argument in this specific instance is to enable the argument's default value to be set in the
% arguments block itself as opposed to doing so in the function body... I like this better but if you didn't
% you could just as easily keep the constructor simple.
methods
function obj = DimensionArgument(A, Adim)
%DimensionArgument Initialize the dimension argument.
arguments
% This will be the enumeration member name from auto-completed entries, or the raw user input if not
% used.
A = [];
% A flag that indicates to create the value using different logic, in this case the first non-singleton
% dimension, because this matches the behavior of functions like, all(), sum() prod(), etc.
Adim (1, 1) logical = false;
end
if(Adim)
% Allows default initialization from an input to match the aforemention function's behavior
obj.Data = firstNonscalarDim(A);
else
% As a convenience for this style of implementation we can validate the input to ensure that since we're
% suppose to be an enumeration, the input is valid
DimensionArgument.mustBeValidMember(A);
% Store the input in a hidden property since declaring ~Enumeration means we are responsible for storing
% it.
obj.Data = A;
end
end
end
%% Conversion methods
% Applies conversion to the data property so that implicit casting of functions works. Unfortunately most of
% the MathWorks defined functions use a different system than that employed by the arguments block, which
% defers to the class defined converter methods... Which is why uplus (+obj) has been defined to unwrap the
% data for ease of use.
methods
function obj = uplus(obj)
obj = obj.Data;
end
function str = char(obj)
str = char(obj.Data);
end
function str = cellstr(obj)
str = cellstr(obj.Data);
end
function str = string(obj)
str = string(obj.Data);
end
function A = double(obj)
A = double(obj.Data);
end
function A = int8(obj)
A = int8(obj.Data);
end
function A = int16(obj)
A = int16(obj.Data);
end
function A = int32(obj)
A = int32(obj.Data);
end
function A = int64(obj)
A = int64(obj.Data);
end
end
%% Validation methods
% These utility methods are for input validation
methods(Static, Access = private)
function tf = isValidMember(obj)
%isValidMember Checks that the input is a valid dimension argument.
tf = (istext(obj) && all(obj == "all", "all")) || (isnumeric(obj) && all(isint(obj) & obj > 0, "all"));
end
function mustBeValidMember(obj)
%mustBeValidMember Validates that the input is a valid dimension argument for the dim/dimVec arguments.
if(~DimensionArgument.isValidMember(obj))
exception("JB:DimensionArgument:InvalidInput", "Input must be an integer value or the term 'all'.")
end
end
end
%% Convenient data display passthrough
methods
function disp(obj, name)
arguments
obj DimensionArgument
name string {mustBeScalarOrEmpty} = [];
end
% Dispatch internal data's display implementation
display(obj.Data, char(name));
end
end
end
In the event you'd actually play with theres here are the function definitions for some of the utility functions I used in them, including my exception would be a pain so i wont, these cases wont use it any...
% Far from my definition isint() but is consistent with mustBeInteger() for real numbers but will suffice for the example
function tf = isint(A)
arguments
A {mustBeNumeric(A)};
end
tf = floor(A) == A
end
% Sort of the same but its fine
function dim = firstNonscalarDim(A)
arguments
A
end
dim = [find(size(A) > 1, 1), 0];
dim(1) = dim(1);
end
Hello MATLAB Central, this is my first article.
My name is Yann. And I love MATLAB.
I also love HTTP (i know, weird fetish)
So i started a conversation with ChatGPT about it:
gitclone('https://github.com/yanndebray/HTTP-with-MATLAB');
cd('HTTP-with-MATLAB')
http_with_MATLAB
I'm not sure that this platform is intended to clone repos from github, but i figured I'd paste this shortcut in case you want to try out my live script http_with_MATLAB.m
A lot of what i program lately relies on external web services (either for fetching data, or calling LLMs).
So I wrote a small tutorial of the 7 or so things I feel like I need to remember when making HTTP requests in MATLAB.
Let me know what you think
Summary:
Dynamically accessing variable names can negatively impact the readability of your code and can cause it to run slower by preventing MATLAB from optimizing it as well as it could if you used alternate techniques. The most common alternative is to use simple and efficient indexing.
Explanation:
Sometimes beginners (and some self-taught professors) think it would be a good idea to dynamically create or access variable names, the variables are often named something like these:
- matrix1, matrix2, matrix3, matrix4, ...
- test_20kmh, test_50kmh, test_80kmh, ...
- nameA, nameB, nameC, nameD,...
Good reasons why dynamic variable names should be avoided:
- Slow
- Buggy
- Security Risk
- Difficult to Work With
- Obfuscated Code Intent
- Confuses Data with Code
- Code Helper Tools do not Work
- Magically Making Variables Appear in a Workspace is Risky
There are much better alternatives to accessing dynamic variable names:
- Indexing into Cell Array or ND-Array
- Non-scalar Structures (with Indexing)
- Dynamic Field-Names in a Structure
- load into a Structure, not into the Workspace
- save the Fields of a Scalar Structure
- Use a table or timetable Array
- Use more Efficient Ways to Pass Variables Between Workspaces
Note that avoiding eval (and assignin, etc.) is not some esoteric MATLAB restriction, it also applies to many other programming languages as well:
MATLAB Documentation:
If you are not interested in reading the answers below then at least read MATLAB's own documentation on this topic Alternatives to the eval Function, which states "A frequent use of the eval function is to create sets of variables such as A1, A2, ..., An, but this approach does not use the array processing power of MATLAB and is not recommended. The preferred method is to store related data in a single array." Data in a single array can be accessed very efficiently using indexing.
Note that all of these problems and disadvantages also apply to functions load (without an output variable), assignin, evalin, and evalc, and the MATLAB documentation explicitly recommends to "Avoid functions such as eval, evalc, evalin, and feval(fname)".
The official MATLAB blogs explain why eval should be avoided, the better alternatives to eval, and clearly recommend against magically creating variables. Using eval comes out at position number one on this list of Top 10 MATLAB Code Practices That Make Me Cry. Experienced MATLAB users recommend avoiding using eval for trivial code, and have written extensively on this topic.
Did you know that function double with string vector input significantly outperforms str2double with the same input:
x = rand(1,50000);
t = string(x);
tic; str2double(t); toc
tic; I1 = str2double(t); toc
tic; I2 = double(t); toc
isequal(I1,I2)
Recently I needed to parse numbers from text. I automatically tried to use str2double. However, profiling revealed that str2double was the main bottleneck in my code. Than I realized that there is a new note (since R2024a) in the documentation of str2double:
"Calling string and then double is recommended over str2double because it provides greater flexibility and allows vectorization. For additional information, see Alternative Functionality."
t = turtle(); % Start a turtle
t.forward(100); % Move forward by 100
t.backward(100); % Move backward by 100
t.left(90); % Turn left by 90 degrees
t.right(90); % Tur right by 90 degrees
t.goto(100, 100); % Move to (100, 100)
t.turnto(90); % Turn to 90 degrees, i.e. north
t.speed(1000); % Set turtle speed as 1000 (default: 500)
t.pen_up(); % Pen up. Turtle leaves no trace.
t.pen_down(); % Pen down. Turtle leaves a trace again.
t.color('b'); % Change line color to 'b'
t.begin_fill(FaceColor, EdgeColor, FaceAlpha); % Start filling
t.end_fill(); % End filling
t.change_icon('person.png'); % Change the icon to 'person.png'
t.clear(); % Clear the Axes
classdef turtle < handle
properties (GetAccess = public, SetAccess = private)
x = 0
y = 0
q = 0
end
properties (SetAccess = public)
speed (1, 1) double = 500
end
properties (GetAccess = private)
speed_reg = 100
n_steps = 20
ax
l
ht
im
is_pen_up = false
is_filling = false
fill_color
fill_alpha
end
methods
function obj = turtle()
figure(Name='MATurtle', NumberTitle='off')
obj.ax = axes(box="on");
hold on,
obj.ht = hgtransform();
icon = flipud(imread('turtle.png'));
obj.im = imagesc(obj.ht, icon, ...
XData=[-30, 30], YData=[-30, 30], ...
AlphaData=(255 - double(rgb2gray(icon)))/255);
obj.l = plot(obj.x, obj.y, 'k');
obj.ax.XLim = [-500, 500];
obj.ax.YLim = [-500, 500];
obj.ax.DataAspectRatio = [1, 1, 1];
obj.ax.Toolbar.Visible = 'off';
disableDefaultInteractivity(obj.ax);
end
function home(obj)
obj.x = 0;
obj.y = 0;
obj.ht.Matrix = eye(4);
end
function forward(obj, dist)
obj.step(dist);
end
function backward(obj, dist)
obj.step(-dist)
end
function step(obj, delta)
if numel(delta) == 1
delta = delta*[cosd(obj.q), sind(obj.q)];
end
if obj.is_filling
obj.fill(delta);
else
obj.move(delta);
end
end
function goto(obj, x, y)
dx = x - obj.x;
dy = y - obj.y;
obj.turnto(rad2deg(atan2(dy, dx)));
obj.step([dx, dy]);
end
function left(obj, q)
obj.turn(q);
end
function right(obj, q)
obj.turn(-q);
end
function turnto(obj, q)
obj.turn(obj.wrap_angle(q - obj.q, -180));
end
function pen_up(obj)
if obj.is_filling
warning('not available while filling')
return
end
obj.is_pen_up = true;
end
function pen_down(obj, go)
if obj.is_pen_up
if nargin == 1
obj.l(end+1) = plot(obj.x, obj.y, Color=obj.l(end).Color);
else
obj.l(end+1) = go;
end
uistack(obj.ht, 'top')
end
obj.is_pen_up = false;
end
function color(obj, line_color)
if obj.is_filling
warning('not available while filling')
return
end
obj.pen_up();
obj.pen_down(plot(obj.x, obj.y, Color=line_color));
end
function begin_fill(obj, FaceColor, EdgeColor, FaceAlpha)
arguments
obj
FaceColor = [.6, .9, .6];
EdgeColor = [0 0.4470 0.7410];
FaceAlpha = 1;
end
if obj.is_filling
warning('already filling')
return
end
obj.fill_color = FaceColor;
obj.fill_alpha = FaceAlpha;
obj.pen_up();
obj.pen_down(patch(obj.x, obj.y, [1, 1, 1], ...
EdgeColor=EdgeColor, FaceAlpha=0));
obj.is_filling = true;
end
function end_fill(obj)
if ~obj.is_filling
warning('not filling now')
return
end
obj.l(end).FaceColor = obj.fill_color;
obj.l(end).FaceAlpha = obj.fill_alpha;
obj.is_filling = false;
end
function change_icon(obj, filename)
icon = flipud(imread(filename));
obj.im.CData = icon;
obj.im.AlphaData = (255 - double(rgb2gray(icon)))/255;
end
function clear(obj)
obj.x = 0;
obj.y = 0;
delete(obj.ax.Children(2:end));
obj.l = plot(0, 0, 'k');
obj.ht.Matrix = eye(4);
end
end
methods (Access = private)
function animated_step(obj, delta, q, initFcn, updateFcn)
arguments
obj
delta
q
initFcn = @() []
updateFcn = @(~, ~) []
end
dx = delta(1)/obj.n_steps;
dy = delta(2)/obj.n_steps;
dq = q/obj.n_steps;
pause_duration = norm(delta)/obj.speed/obj.speed_reg;
initFcn();
for i = 1:obj.n_steps
updateFcn(dx, dy);
obj.ht.Matrix = makehgtform(...
translate=[obj.x + dx*i, obj.y + dy*i, 0], ...
zrotate=deg2rad(obj.q + dq*i));
pause(pause_duration)
drawnow limitrate
end
obj.x = obj.x + delta(1);
obj.y = obj.y + delta(2);
end
function obj = turn(obj, q)
obj.animated_step([0, 0], q);
obj.q = obj.wrap_angle(obj.q + q, 0);
end
function move(obj, delta)
initFcn = @() [];
updateFcn = @(dx, dy) [];
if ~obj.is_pen_up
initFcn = @() initializeLine();
updateFcn = @(dx, dy) obj.update_end_point(obj.l(end), dx, dy);
end
function initializeLine()
obj.l(end).XData(end+1) = obj.l(end).XData(end);
obj.l(end).YData(end+1) = obj.l(end).YData(end);
end
obj.animated_step(delta, 0, initFcn, updateFcn);
end
function obj = fill(obj, delta)
initFcn = @() initializePatch();
updateFcn = @(dx, dy) obj.update_end_point(obj.l(end), dx, dy);
function initializePatch()
obj.l(end).Vertices(end+1, :) = obj.l(end).Vertices(end, :);
obj.l(end).Faces = 1:size(obj.l(end).Vertices, 1);
end
obj.animated_step(delta, 0, initFcn, updateFcn);
end
end
methods (Static, Access = private)
function update_end_point(l, dx, dy)
l.XData(end) = l.XData(end) + dx;
l.YData(end) = l.YData(end) + dy;
end
function q = wrap_angle(q, min_angle)
q = mod(q - min_angle, 360) + min_angle;
end
end
end
I am thrilled python interoperability now seems to work for me with my APPLE M1 MacBookPro and MATLAB V2025a. The available instructions are still, shall we say, cryptic. Here is a summary of my interaction with GPT 4o to get this to work.
===========================================================
MATLAB R2025a + Python (Astropy) Integration on Apple Silicon (M1/M2/M3 Macs)
===========================================================
Author: D. Carlsmith, documented with ChatGPT
Last updated: July 2025
This guide provides full instructions, gotchas, and workarounds to run Python 3.10 with MATLAB R2025a (Apple Silicon/macOS) using native ARM64 Python and calling modules like Astropy, Numpy, etc. from within MATLAB.
===========================================================
Overview
===========================================================
- MATLAB R2025a on Apple Silicon (M1/M2/M3) runs as "maca64" (native ARM64).
- To call Python from MATLAB, the Python interpreter must match that architecture (ARM64).
- Using Intel Python (x86_64) with native MATLAB WILL NOT WORK.
- The cleanest solution: use Miniforge3 (Conda-forge's lightweight ARM64 distribution).
===========================================================
1. Install Miniforge3 (ARM64-native Conda)
===========================================================
In Terminal, run:
curl -LO https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-MacOSX-arm64.sh
bash Miniforge3-MacOSX-arm64.sh
Follow prompts:
- Press ENTER to scroll through license.
- Type "yes" when asked to accept the license.
- Press ENTER to accept the default install location: ~/miniforge3
- When asked:
Do you wish to update your shell profile to automatically initialize conda? [yes|no]
Type: yes
===========================================================
2. Restart Terminal and Create a Python Environment for MATLAB
===========================================================
Run the following:
conda create -n matlab python=3.10 astropy numpy -y
conda activate matlab
Verify the Python path:
which python
Expected output:
/Users/YOURNAME/miniforge3/envs/matlab/bin/python
===========================================================
3. Verify Python + Astropy From Terminal
===========================================================
Run:
python -c "import astropy; print(astropy.__version__)"
Expected output:
6.x.x (or similar)
===========================================================
4. Configure MATLAB to Use This Python
===========================================================
In MATLAB R2025a (Apple Silicon):
clear classes
pyenv('Version', '/Users/YOURNAME/miniforge3/envs/matlab/bin/python')
py.sys.version
You should see the Python version printed (e.g. 3.10.18). No error means it's working.
===========================================================
5. Gotchas and Their Solutions
===========================================================
❌ Error: Python API functions are not available
→ Cause: Wrong architecture or broken .dylib
→ Fix: Use Miniforge ARM64 Python. DO NOT use Intel Anaconda.
❌ Error: Invalid text character (↑ points at __version__)
→ Cause: MATLAB can’t parse double underscores typed or pasted
→ Fix: Use: py.getattr(module, '__version__')
❌ Error: Unrecognized method 'separation' or 'sec'
→ Cause: MATLAB can't reflect dynamic Python methods
→ Fix: Use: py.getattr(obj, 'method')(args)
===========================================================
6. Run Full Verification in MATLAB
===========================================================
Paste this into MATLAB:
% Set environment
clear classes
pyenv('Version', '/Users/YOURNAME/miniforge3/envs/matlab/bin/python');
% Import modules
coords = py.importlib.import_module('astropy.coordinates');
time_mod = py.importlib.import_module('astropy.time');
table_mod = py.importlib.import_module('astropy.table');
% Astropy version
ver = char(py.getattr(py.importlib.import_module('astropy'), '__version__'));
disp(['Astropy version: ', ver]);
% SkyCoord angular separation
c1 = coords.SkyCoord('10h21m00s', '+41d12m00s', pyargs('frame', 'icrs'));
c2 = coords.SkyCoord('10h22m00s', '+41d15m00s', pyargs('frame', 'icrs'));
sep_fn = py.getattr(c1, 'separation');
sep = sep_fn(c2);
arcsec = double(sep.to('arcsec').value);
fprintf('Angular separation = %.3f arcsec\n', arcsec);
% Time difference in seconds
Time = time_mod.Time;
t1 = Time('2025-01-01T00:00:00', pyargs('format','isot','scale','utc'));
t2 = Time('2025-01-02T00:00:00', pyargs('format','isot','scale','utc'));
dt = py.getattr(t2, '__sub__')(t1);
seconds = double(py.getattr(dt, 'sec'));
fprintf('Time difference = %.0f seconds\n', seconds);
% Astropy table display
tbl = table_mod.Table(pyargs('names', {'a','b'}, 'dtype', {'int','float'}));
tbl.add_row({1, 2.5});
tbl.add_row({2, 3.7});
disp(tbl);
===========================================================
7. Optional: Automatically Configure Python in startup.m
===========================================================
To avoid calling pyenv() every time, edit your MATLAB startup:
edit startup.m
Add:
try
pyenv('Version', '/Users/YOURNAME/miniforge3/envs/matlab/bin/python');
catch
warning("Python already loaded.");
end
===========================================================
8. Final Notes
===========================================================
- This setup avoids all architecture mismatches.
- It uses a clean, minimal ARM64 Python that integrates seamlessly with MATLAB.
- Do not mix Anaconda (Intel) with Apple Silicon MATLAB.
- Use py.getattr for any Python attribute containing underscores or that MATLAB can't resolve.
You can now run NumPy, Astropy, Pandas, Astroquery, Matplotlib, and more directly from MATLAB.
===========================================================
I would like to zoom directly on the selected region when using 
on my image created with image or imagesc. First of all, I would recommend using image or imagesc and not imshow for this case, see comparison here: Differences between imshow() and image()? However when zooming Stretch-to-Fill behavior happens and I don't want that. Try range zoom to image generated by this code:
on my image created with image or imagesc. First of all, I would recommend using image or imagesc and not imshow for this case, see comparison here: Differences between imshow() and image()? However when zooming Stretch-to-Fill behavior happens and I don't want that. Try range zoom to image generated by this code:fig = uifigure;
ax = uiaxes(fig);
im = imread("peppers.png");
h = imagesc(im,"Parent",ax);
axis(ax,'tight', 'off')
I can fix that with manualy setting data aspect ratio:
daspect(ax,[1 1 1])
However, I need this code to run automatically after zooming. So I create zoom object and ActionPostCallback which is called everytime after I zoom, see zoom - ActionPostCallback.
z = zoom(ax);
z.ActionPostCallback = @(fig,ax) daspect(ax.Axes,[1 1 1]);
If you need, you can also create ActionPreCallback which is called everytime before I zoom, see zoom - ActionPreCallback.
z.ActionPreCallback = @(fig,ax) daspect(ax.Axes,'auto');
Code written and run in R2025a.
The beautiful and elegant chord diagrams were all created using MATLAB?
Indeed, they were all generated using the chord diagram plotting toolkit that I developed myself:
- - Chord chart: [chord chart](https://www.mathworks.com/matlabcentral/fileexchange/116550-chord-chart)
- - Directed graph chord chart: [digraph chord chart]:(https://www.mathworks.com/matlabcentral/fileexchange/121043-digraph-chord-chart)
You can download these toolkits from the provided links.
The reason for writing this article is that many people have started using the chord diagram plotting toolkit that I developed. However, some users are unsure about customizing certain styles. As the developer, I have a good understanding of the implementation principles of the toolkit and can apply it flexibly. This has sparked the idea of challenging myself to create various styles of chord diagrams. Currently, the existing code is quite lengthy. In the future, I may integrate some of this code into the toolkit, enabling users to achieve the effects of many lines of code with just a few lines.
Without further ado, let's see the extent to which this MATLAB toolkit can currently perform.
demo 1
rng(2)
dataMat = randi([0,5], [11,5]);
dataMat(1:6,1) = 0;
dataMat([11,7],1) = [45,25];
dataMat([1,4,5,7],2) = [20,20,30,30];
dataMat(:,3) = 0;
dataMat(6,3) = 45;
dataMat(1:5,4) = 0;
dataMat([6,7],4) = [25,25];
dataMat([5,6,9],5) = [25,25,25];
colName = {'Fly', 'Beetle', 'Leaf', 'Soil', 'Waxberry'};
rowName = {'Bartomella', 'Bradyrhizobium', 'Dysgomonas', 'Enterococcus',...
'Lactococcus', 'norank', 'others', 'Pseudomonas', 'uncultured',...
'Vibrionimonas', 'Wolbachia'};
figure('Units','normalized', 'Position',[.02,.05,.6,.85])
CC = chordChart(dataMat, 'rowName',rowName, 'colName',colName, 'Sep',1/80);
CC = CC.draw();
% 修改上方方块颜色(Modify the color of the blocks above)
CListT = [0.7765 0.8118 0.5216; 0.4431 0.4706 0.3843; 0.5804 0.2275 0.4549;
0.4471 0.4039 0.6745; 0.0157 0 0 ];
for i = 1:size(dataMat, 2)
CC.setSquareT_N(i, 'FaceColor',CListT(i,:))
end
% 修改下方方块颜色(Modify the color of the blocks below)
CListF = [0.5843 0.6863 0.7843; 0.1098 0.1647 0.3255; 0.0902 0.1608 0.5373;
0.6314 0.7961 0.2118; 0.0392 0.2078 0.1059; 0.0157 0 0 ;
0.8549 0.9294 0.8745; 0.3882 0.3255 0.4078; 0.5020 0.7216 0.3843;
0.0902 0.1843 0.1804; 0.8196 0.2314 0.0706];
for i = 1:size(dataMat, 1)
CC.setSquareF_N(i, 'FaceColor',CListF(i,:))
end
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
CC.setChordMN(i,j, 'FaceColor',CListT(j,:), 'FaceAlpha',.5)
end
end
CC.tickState('on')
CC.labelRotate('on')
CC.setFont('FontSize',17, 'FontName','Cambria')
% CC.labelRotate('off')
% textHdl = findobj(gca,'Tag','ChordLabel');
% for i = 1:length(textHdl)
% if textHdl(i).Position(2) < 0
% if abs(textHdl(i).Position(1)) > .7
% textHdl(i).Rotation = textHdl(i).Rotation + 45;
% textHdl(i).HorizontalAlignment = 'right';
% if textHdl(i).Rotation > 90
% textHdl(i).Rotation = textHdl(i).Rotation + 180;
% textHdl(i).HorizontalAlignment = 'left';
% end
% else
% textHdl(i).Rotation = textHdl(i).Rotation + 10;
% textHdl(i).HorizontalAlignment = 'right';
% end
% end
% end
demo 2
rng(3)
dataMat = randi([1,15], [7,22]);
dataMat(dataMat < 11) = 0;
dataMat(1, sum(dataMat, 1) == 0) = 15;
colName = {'A2M', 'FGA', 'FGB', 'FGG', 'F11', 'KLKB1', 'SERPINE1', 'VWF',...
'THBD', 'TFPI', 'PLAT', 'SERPINA5', 'SERPIND1', 'F2', 'PLG', 'F12',...
'SERPINC1', 'SERPINA1', 'PROS1', 'SERPINF2', 'F13A1', 'PROC'};
rowName = {'Lung', 'Spleen', 'Liver', 'Heart',...
'Renal cortex', 'Renal medulla', 'Thyroid'};
figure('Units','normalized', 'Position',[.02,.05,.6,.85])
CC = chordChart(dataMat, 'rowName',rowName, 'colName',colName, 'Sep',1/80, 'LRadius',1.21);
CC = CC.draw();
CC.labelRotate('on')
% 单独设置每一个弦末端方块(Set individual end blocks for each chord)
% Use obj.setEachSquareF_Prop
% or obj.setEachSquareT_Prop
% F means from (blocks below)
% T means to (blocks above)
CListT = [173,70,65; 79,135,136]./255;
% Upregulated:1 | Downregulated:2
Regulated = rand([7, 22]);
Regulated = (Regulated < .8) + 1;
for i = 1:size(Regulated, 1)
for j = 1:size(Regulated, 2)
CC.setEachSquareT_Prop(i, j, 'FaceColor', CListT(Regulated(i,j),:))
end
end
% 绘制图例(Draw legend)
H1 = fill([0,1,0] + 100, [1,0,1] + 100, CListT(1,:), 'EdgeColor','none');
H2 = fill([0,1,0] + 100, [1,0,1] + 100, CListT(2,:), 'EdgeColor','none');
lgdHdl = legend([H1,H2], {'Upregulated','Downregulated'}, 'AutoUpdate','off', 'Location','best');
lgdHdl.ItemTokenSize = [12,12];
lgdHdl.Box = 'off';
lgdHdl.FontSize = 13;
% 修改下方方块颜色(Modify the color of the blocks below)
CListF = [128,108,171; 222,208,161; 180,196,229; 209,150,146; 175,201,166;
134,156,118; 175,175,173]./255;
for i = 1:size(dataMat, 1)
CC.setSquareF_N(i, 'FaceColor',CListF(i,:))
end
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
CC.setChordMN(i,j, 'FaceColor',CListF(i,:), 'FaceAlpha',.45)
end
end
demo 3
dataMat = rand([15,15]);
dataMat(dataMat > .15) = 0;
CList = [ 75,146,241; 252,180, 65; 224, 64, 10; 5,100,146; 191,191,191;
26, 59,105; 255,227,130; 18,156,221; 202,107, 75; 0, 92,219;
243,210,136; 80, 99,129; 241,185,168; 224,131, 10; 120,147,190]./255;
figure('Units','normalized', 'Position',[.02,.05,.6,.85])
BCC = biChordChart(dataMat, 'Arrow','on', 'CData',CList);
BCC = BCC.draw();
% 添加刻度
BCC.tickState('on')
% 修改字体，字号及颜色
BCC.setFont('FontName','Cambria', 'FontSize',17, 'Color',[0,0,.8])
demo 4
rng(5)
dataMat = randi([1,20], [5,5]);
dataMat(1,1) = 110;
dataMat(2,2) = 40;
dataMat(3,3) = 50;
dataMat(5,5) = 50;
CList1 = [164,190,158; 216,213,153; 177,192,208; 238,238,227; 249,217,153]./255;
CList2 = [247,204,138; 128,187,185; 245,135,124; 140,199,197; 252,223,164]./255;
CList = CList2;
NameList={'CHORD','CHART','MADE','BY','SLANDARER'};
figure('Units','normalized', 'Position',[.02,.05,.6,.85])
BCC = biChordChart(dataMat, 'Arrow','on', 'CData',CList, 'Sep',1/30, 'Label',NameList, 'LRadius',1.33);
BCC = BCC.draw();
% 添加刻度
BCC.tickState('on')
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
if dataMat(i,j) > 0
BCC.setChordMN(i,j, 'FaceAlpha',.7, 'EdgeColor',CList(i,:)./1.1)
end
end
end
% 修改方块颜色(Modify the color of the blocks)
for i = 1:size(dataMat, 1)
BCC.setSquareN(i, 'EdgeColor',CList(i,:)./1.7)
end
% 修改字体，字号及颜色
BCC.setFont('FontName','Cambria', 'FontSize',17)
BCC.tickLabelState('on')
BCC.setTickFont('FontName','Cambria', 'FontSize',9)
demo 5
dataMat=randi([1,20], [14,3]);
dataMat(11:14,1) = 0;
dataMat(6:10,2) = 0;
dataMat(1:5,3) = 0;
colName = compose('C%d', 1:3);
rowName = [compose('A%d', 1:7), compose('B%d', 7:-1:1)];
figure('Units','normalized', 'Position',[.02,.05,.6,.85])
CC = chordChart(dataMat, 'rowName',rowName, 'colName',colName, 'Sep',1/80);
CC = CC.draw();
% 修改上方方块颜色(Modify the color of the blocks above)
for i = 1:size(dataMat, 2)
CC.setSquareT_N(i, 'FaceColor',[190,190,190]./255)
end
% 修改下方方块颜色(Modify the color of the blocks below)
CListF=[255,244,138; 253,220,117; 254,179, 78; 253,190, 61;
252, 78, 41; 228, 26, 26; 178, 0, 36; 4, 84,119;
1,113,137; 21,150,155; 67,176,173; 68,173,158;
123,204,163; 184,229,162]./255;
for i = 1:size(dataMat, 1)
CC.setSquareF_N(i, 'FaceColor',CListF(i,:))
end
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
CC.setChordMN(i,j, 'FaceColor',CListF(i,:), 'FaceAlpha',.5)
end
end
CC.tickState('on')
CC.tickLabelState('on')
demo 6
rng(2)
dataMat = randi([0,40], [20,4]);
dataMat(rand([20,4]) < .2) = 0;
dataMat(1,3) = 500;
dataMat(20,1:4) = [140; 150; 80; 90];
colName = compose('T%d', 1:4);
rowName = compose('SL%d', 1:20);
figure('Units','normalized', 'Position',[.02,.05,.6,.85])
CC = chordChart(dataMat, 'rowName',rowName, 'colName',colName, 'Sep',1/80, 'LRadius',1.23);
CC = CC.draw();
% 修改上方方块颜色(Modify the color of the blocks above)
CListT = [0.62,0.49,0.27; 0.28,0.57,0.76
0.25,0.53,0.30; 0.86,0.48,0.34];
for i = 1:size(dataMat, 2)
CC.setSquareT_N(i, 'FaceColor',CListT(i,:))
end
% 修改下方方块颜色(Modify the color of the blocks below)
CListF = [0.94,0.84,0.60; 0.16,0.50,0.67; 0.92,0.62,0.49;
0.48,0.44,0.60; 0.48,0.44,0.60; 0.71,0.79,0.73;
0.96,0.98,0.98; 0.51,0.82,0.95; 0.98,0.70,0.82;
0.97,0.85,0.84; 0.55,0.64,0.62; 0.94,0.93,0.60;
0.98,0.90,0.85; 0.72,0.84,0.81; 0.85,0.45,0.49;
0.76,0.76,0.84; 0.59,0.64,0.62; 0.62,0.14,0.15;
0.75,0.75,0.75; 1.00,1.00,1.00];
for i = 1:size(dataMat, 1)
CC.setSquareF_N(i, 'FaceColor',CListF(i,:))
end
CC.setSquareF_N(size(dataMat, 1), 'EdgeColor','k', 'LineWidth',1)
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
CC.setChordMN(i,j, 'FaceColor',CListT(j,:), 'FaceAlpha',.46)
end
end
CC.tickState('on')
CC.labelRotate('on')
CC.setFont('FontSize',17, 'FontName','Cambria')
demo 7
dataMat = randi([10,10000], [10,10]);
dataMat(6:10,:) = 0;
dataMat(:,1:5) = 0;
NameList = {'BOC', 'ICBC', 'ABC', 'BOCM', 'CCB', ...
'yama', 'nikoto', 'saki', 'koto', 'kawa'};
CList = [0.63,0.75,0.88
0.67,0.84,0.75
0.85,0.78,0.88
1.00,0.92,0.93
0.92,0.63,0.64
0.57,0.67,0.75
1.00,0.65,0.44
0.72,0.73,0.40
0.65,0.57,0.58
0.92,0.94,0.96];
figure('Units','normalized', 'Position',[.02,.05,.6,.85])
BCC = biChordChart(dataMat, 'Arrow','on', 'CData',CList, 'Label',NameList);
BCC = BCC.draw();
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
if dataMat(i,j) > 0
BCC.setChordMN(i,j, 'FaceAlpha',.85, 'EdgeColor',CList(i,:)./1.5, 'LineWidth',.8)
end
end
end
for i = 1:size(dataMat, 1)
BCC.setSquareN(i, 'EdgeColor',CList(i,:)./1.5, 'LineWidth',1)
end
% 添加刻度、修改字体
BCC.tickState('on')
BCC.setFont('FontName','Cambria', 'FontSize',17)
demo 8
dataMat = rand([11,4]);
dataMat = round(10.*dataMat.*((11:-1:1).'+1))./10;
colName = {'A','B','C','D'};
rowName = {'Acidobacteriota', 'Actinobacteriota', 'Proteobacteria', ...
'Chloroflexi', 'Bacteroidota', 'Firmicutes', 'Gemmatimonadota', ...
'Verrucomicrobiota', 'Patescibacteria', 'Planctomyetota', 'Others'};
figure('Units','normalized', 'Position',[.02,.05,.8,.85])
CC = chordChart(dataMat, 'colName',colName, 'Sep',1/80, 'SSqRatio',30/100);% -30/100
CC = CC.draw();
% 修改上方方块颜色(Modify the color of the blocks above)
CListT = [0.93,0.60,0.62
0.55,0.80,0.99
0.95,0.82,0.18
1.00,0.81,0.91];
for i = 1:size(dataMat, 2)
CC.setSquareT_N(i, 'FaceColor',CListT(i,:))
end
% 修改下方方块颜色(Modify the color of the blocks below)
CListF = [0.75,0.73,0.86
0.56,0.83,0.78
0.00,0.60,0.20
1.00,0.49,0.02
0.78,0.77,0.95
0.59,0.24,0.36
0.98,0.51,0.45
0.96,0.55,0.75
0.47,0.71,0.84
0.65,0.35,0.16
0.40,0.00,0.64];
for i = 1:size(dataMat, 1)
CC.setSquareF_N(i, 'FaceColor',CListF(i,:))
end
% 修改弦颜色(Modify chord color)
CListC = [0.55,0.83,0.76
0.75,0.73,0.86
0.00,0.60,0.19
1.00,0.51,0.04];
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
CC.setChordMN(i,j, 'FaceColor',CListC(j,:), 'FaceAlpha',.4)
end
end
% 单独设置每一个弦末端方块(Set individual end blocks for each chord)
% Use obj.setEachSquareF_Prop
% or obj.setEachSquareT_Prop
% F means from (blocks below)
% T means to (blocks above)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
CC.setEachSquareT_Prop(i,j, 'FaceColor', CListF(i,:))
end
end
% 添加刻度
CC.tickState('on')
% 修改字体，字号及颜色
CC.setFont('FontName','Cambria', 'FontSize',17)
% 隐藏下方标签
textHdl = findobj(gca, 'Tag','ChordLabel');
for i = 1:length(textHdl)
if textHdl(i).Position(2) < 0
set(textHdl(i), 'Visible','off')
end
end
% 绘制图例(Draw legend)
scatterHdl = scatter(10.*ones(size(dataMat,1)),10.*ones(size(dataMat,1)), ...
55, 'filled');
for i = 1:length(scatterHdl)
scatterHdl(i).CData = CListF(i,:);
end
lgdHdl = legend(scatterHdl, rowName, 'Location','best', 'FontSize',16, 'FontName','Cambria', 'Box','off');
set(lgdHdl, 'Position',[.7482,.3577,.1658,.3254])
demo 9
dataMat = randi([0,10], [5,5]);
CList1 = [0.70,0.59,0.67
0.62,0.70,0.62
0.81,0.75,0.62
0.80,0.62,0.56
0.62,0.65,0.65];
CList2 = [0.02,0.02,0.02
0.59,0.26,0.33
0.38,0.49,0.38
0.03,0.05,0.03
0.29,0.28,0.32];
CList = CList2;
NameList={'CHORD','CHART','MADE','BY','SLANDARER'};
figure('Units','normalized', 'Position',[.02,.05,.6,.85])
BCC = biChordChart(dataMat, 'Arrow','on', 'CData',CList, 'Sep',1/30, 'Label',NameList, 'LRadius',1.33);
BCC = BCC.draw();
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
BCC.setChordMN(i,j, 'FaceAlpha',.5)
end
end
% 修改方块颜色(Modify the color of the blocks)
for i = 1:size(dataMat, 1)
BCC.setSquareN(i, 'EdgeColor',[0,0,0], 'LineWidth',5)
end
% 添加刻度
BCC.tickState('on')
% 修改字体，字号及颜色
BCC.setFont('FontSize',17, 'FontWeight','bold')
BCC.tickLabelState('on')
BCC.setTickFont('FontSize',9)
demo 10
rng(2)
dataMat = rand([14,5]) > .3;
colName = {'phosphorylation', 'vasculature development', 'blood vessel development', ...
'cell adhesion', 'plasma membrane'};
rowName = {'THY1', 'FGF2', 'MAP2K1', 'CDH2', 'HBEGF', 'CXCR4', 'ECSCR',...
'ACVRL1', 'RECK', 'PNPLA6', 'CDH5', 'AMOT', 'EFNB2', 'CAV1'};
figure('Units','normalized', 'Position',[.02,.05,.9,.85])
CC = chordChart(dataMat, 'colName',colName, 'rowName',rowName, 'Sep',1/80, 'LRadius',1.2);
CC = CC.draw();
% 修改上方方块颜色(Modify the color of the blocks above)
CListT1 = [0.5686 0.1961 0.2275
0.2275 0.2863 0.3765
0.8431 0.7882 0.4118
0.4275 0.4510 0.2706
0.3333 0.2706 0.2510];
CListT2 = [0.4941 0.5490 0.4118
0.9059 0.6510 0.3333
0.8980 0.6157 0.4980
0.8902 0.5137 0.4667
0.4275 0.2824 0.2784];
CListT3 = [0.4745 0.5843 0.7569
0.4824 0.5490 0.5843
0.6549 0.7216 0.6510
0.9412 0.9216 0.9059
0.9804 0.7608 0.6863];
CListT = CListT3;
for i = 1:size(dataMat, 2)
CC.setSquareT_N(i, 'FaceColor',CListT(i,:), 'EdgeColor',[0,0,0])
end
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
CC.setChordMN(i,j, 'FaceColor',CListT(j,:), 'FaceAlpha',.9, 'EdgeColor',[0,0,0])
end
end
% 修改下方方块颜色(Modify the color of the blocks below)
logFC = sort(rand(1,14))*6 - 3;
for i = 1:size(dataMat, 1)
CC.setSquareF_N(i, 'CData',logFC(i), 'FaceColor','flat', 'EdgeColor',[0,0,0])
end
CMap = [ 0 0 1.0000; 0.0645 0.0645 1.0000; 0.1290 0.1290 1.0000; 0.1935 0.1935 1.0000
0.2581 0.2581 1.0000; 0.3226 0.3226 1.0000; 0.3871 0.3871 1.0000; 0.4516 0.4516 1.0000
0.5161 0.5161 1.0000; 0.5806 0.5806 1.0000; 0.6452 0.6452 1.0000; 0.7097 0.7097 1.0000
0.7742 0.7742 1.0000; 0.8387 0.8387 1.0000; 0.9032 0.9032 1.0000; 0.9677 0.9677 1.0000
1.0000 0.9677 0.9677; 1.0000 0.9032 0.9032; 1.0000 0.8387 0.8387; 1.0000 0.7742 0.7742
1.0000 0.7097 0.7097; 1.0000 0.6452 0.6452; 1.0000 0.5806 0.5806; 1.0000 0.5161 0.5161
1.0000 0.4516 0.4516; 1.0000 0.3871 0.3871; 1.0000 0.3226 0.3226; 1.0000 0.2581 0.2581
1.0000 0.1935 0.1935; 1.0000 0.1290 0.1290; 1.0000 0.0645 0.0645; 1.0000 0 0];
colormap(CMap);
try clim([-3,3]),catch,end
try caxis([-3,3]),catch,end
CBHdl = colorbar();
CBHdl.Position = [0.74,0.25,0.02,0.2];
% =========================================================================
% 交换XY轴(Swap XY axis)
patchHdl = findobj(gca, 'Type','patch');
for i = 1:length(patchHdl)
tX = patchHdl(i).XData;
tY = patchHdl(i).YData;
patchHdl(i).XData = tY;
patchHdl(i).YData = - tX;
end
txtHdl = findobj(gca, 'Type','text');
for i = 1:length(txtHdl)
txtHdl(i).Position([1,2]) = [1,-1].*txtHdl(i).Position([2,1]);
if txtHdl(i).Position(1) < 0
txtHdl(i).HorizontalAlignment = 'right';
else
txtHdl(i).HorizontalAlignment = 'left';
end
end
lineHdl = findobj(gca, 'Type','line');
for i = 1:length(lineHdl)
tX = lineHdl(i).XData;
tY = lineHdl(i).YData;
lineHdl(i).XData = tY;
lineHdl(i).YData = - tX;
end
% =========================================================================
txtHdl = findobj(gca, 'Type','text');
for i = 1:length(txtHdl)
if txtHdl(i).Position(1) > 0
txtHdl(i).Visible = 'off';
end
end
text(1.25,-.15, 'LogFC', 'FontSize',16)
text(1.25,1, 'Terms', 'FontSize',16)
patchHdl = [];
for i = 1:size(dataMat, 2)
patchHdl(i) = fill([10,11,12],[10,13,13], CListT(i,:), 'EdgeColor',[0,0,0]);
end
lgdHdl = legend(patchHdl, colName, 'Location','best', 'FontSize',14, 'FontName','Cambria', 'Box','off');
lgdHdl.Position = [.735,.53,.167,.27];
lgdHdl.ItemTokenSize = [18,8];
demo 11
rng(2)
dataMat = rand([12,12]);
dataMat(dataMat < .85) = 0;
dataMat(7,:) = 1.*(rand(1,12)+.1);
dataMat(11,:) = .6.*(rand(1,12)+.1);
dataMat(12,:) = [2.*(rand(1,10)+.1), 0, 0];
CList = [repmat([49,49,49],[10,1]); 235,28,34; 19,146,241]./255;
figure('Units','normalized', 'Position',[.02,.05,.6,.85])
BCC = biChordChart(dataMat, 'Arrow','off', 'CData',CList);
BCC = BCC.draw();
% 添加刻度
BCC.tickState('on')
% 修改字体，字号及颜色
BCC.setFont('FontName','Cambria', 'FontSize',17)
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
if dataMat(i,j) > 0
BCC.setChordMN(i,j, 'FaceAlpha',.78, 'EdgeColor',[0,0,0])
end
end
end
% 修改方块颜色(Modify the color of the blocks)
for i = 1:size(dataMat, 1)
BCC.setSquareN(i, 'EdgeColor',[0,0,0], 'LineWidth',2)
end
demo 12
dataMat = rand([9,9]);
dataMat(dataMat > .7) = 0;
dataMat(eye(9) == 1) = (rand([1,9])+.2).*3;
CList = [0.85,0.23,0.24
0.96,0.39,0.18
0.98,0.63,0.22
0.99,0.80,0.26
0.70,0.76,0.21
0.24,0.74,0.71
0.27,0.65,0.84
0.09,0.37,0.80
0.64,0.40,0.84];
figure('Units','normalized', 'Position',[.02,.05,.6,.85])
BCC = biChordChart(dataMat, 'Arrow','on', 'CData',CList);
BCC = BCC.draw();
% 添加刻度、刻度标签
BCC.tickState('on')
% 修改字体，字号及颜色
BCC.setFont('FontName','Cambria', 'FontSize',17)
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
if dataMat(i,j) > 0
BCC.setChordMN(i,j, 'FaceAlpha',.7)
end
end
end
demo 13
rng(2)
dataMat = randi([1,40], [7,4]);
dataMat(rand([7,4]) < .1) = 0;
colName = compose('MATLAB%d', 1:4);
rowName = compose('SL%d', 1:7);
figure('Units','normalized', 'Position',[.02,.05,.7,.85])
CC = chordChart(dataMat, 'rowName',rowName, 'colName',colName, 'Sep',1/80, 'LRadius',1.32);
CC = CC.draw();
% 修改上方方块颜色(Modify the color of the blocks above)
CListT = [0.49,0.64,0.53
0.75,0.39,0.35
0.80,0.74,0.42
0.40,0.55,0.66];
for i = 1:size(dataMat, 2)
CC.setSquareT_N(i, 'FaceColor',CListT(i,:))
end
% 修改下方方块颜色(Modify the color of the blocks below)
CListF = [0.91,0.91,0.97
0.62,0.95,0.66
0.91,0.61,0.20
0.54,0.45,0.82
0.99,0.76,0.81
0.91,0.85,0.83
0.53,0.42,0.43];
for i = 1:size(dataMat, 1)
CC.setSquareF_N(i, 'FaceColor',CListF(i,:))
end
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
CC.setChordMN(i,j, 'FaceColor',CListT(j,:), 'FaceAlpha',.46)
end
end
CC.tickState('on')
CC.tickLabelState('on')
CC.setFont('FontSize',17, 'FontName','Cambria')
CC.setTickFont('FontSize',8, 'FontName','Cambria')
% 绘制图例(Draw legend)
scatterHdl = scatter(10.*ones(size(dataMat,1)),10.*ones(size(dataMat,1)), ...
55, 'filled');
for i = 1:length(scatterHdl)
scatterHdl(i).CData = CListF(i,:);
end
lgdHdl = legend(scatterHdl, rowName, 'Location','best', 'FontSize',16, 'FontName','Cambria', 'Box','off');
set(lgdHdl, 'Position',[.77,.38,.1658,.27])
demo 14
rng(6)
dataMat = randi([1,20], [8,8]);
dataMat(dataMat > 5) = 0;
dataMat(1,:) = randi([1,15], [1,8]);
dataMat(1,8) = 40;
dataMat(8,8) = 60;
dataMat = dataMat./sum(sum(dataMat));
CList = [0.33,0.53,0.86
0.94,0.50,0.42
0.92,0.58,0.30
0.59,0.47,0.45
0.37,0.76,0.82
0.82,0.68,0.29
0.75,0.62,0.87
0.43,0.69,0.57];
NameList={'CHORD', 'CHART', 'AND', 'BICHORD',...
'CHART', 'MADE', 'BY', 'SLANDARER'};
figure('Units','normalized', 'Position',[.02,.05,.6,.85])
BCC = biChordChart(dataMat, 'Arrow','on', 'CData',CList, 'Sep',1/12, 'Label',NameList, 'LRadius',1.33);
BCC = BCC.draw();
% 添加刻度
BCC.tickState('on')
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
if dataMat(i,j) > 0
BCC.setChordMN(i,j, 'FaceAlpha',.7, 'EdgeColor',CList(i,:)./1.1)
end
end
end
% 修改方块颜色(Modify the color of the blocks)
for i = 1:size(dataMat, 1)
BCC.setSquareN(i, 'EdgeColor',CList(i,:)./1.7)
end
% 修改字体，字号及颜色
BCC.setFont('FontName','Cambria', 'FontSize',17)
BCC.tickLabelState('on')
BCC.setTickFont('FontName','Cambria', 'FontSize',9)
% 调整数值字符串格式
% Adjust numeric string format
BCC.setTickLabelFormat(@(x)[num2str(round(x*100)),'%'])
demo 15
CList = [0.81,0.72,0.83
0.69,0.82,0.89
0.17,0.44,0.64
0.70,0.85,0.55
0.03,0.57,0.13
0.97,0.67,0.64
0.84,0.09,0.12
1.00,0.80,0.46
0.98,0.52,0.01
];
figure('Units','normalized', 'Position',[.02,.05,.53,.85], 'Color',[1,1,1])
% =========================================================================
ax1 = axes('Parent',gcf, 'Position',[0,1/2,1/2,1/2]);
dataMat = rand([9,9]);
dataMat(dataMat > .4) = 0;
BCC = biChordChart(dataMat, 'Arrow','on', 'CData',CList);
BCC = BCC.draw();
BCC.tickState('on')
BCC.setFont('Visible','off')
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
if dataMat(i,j) > 0
BCC.setChordMN(i,j, 'FaceAlpha',.6)
end
end
end
text(-1.2,1.2, 'a', 'FontName','Times New Roman', 'FontSize',35)
% =========================================================================
ax2 = axes('Parent',gcf, 'Position',[1/2,1/2,1/2,1/2]);
dataMat = rand([9,9]);
dataMat(dataMat > .4) = 0;
dataMat = dataMat.*(1:9);
BCC = biChordChart(dataMat, 'Arrow','on', 'CData',CList);
BCC = BCC.draw();
BCC.tickState('on')
BCC.setFont('Visible','off')
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
if dataMat(i,j) > 0
BCC.setChordMN(i,j, 'FaceAlpha',.6)
end
end
end
text(-1.2,1.2, 'b', 'FontName','Times New Roman', 'FontSize',35)
% =========================================================================
ax3 = axes('Parent',gcf, 'Position',[0,0,1/2,1/2]);
dataMat = rand([9,9]);
dataMat(dataMat > .4) = 0;
dataMat = dataMat.*(1:9).';
BCC = biChordChart(dataMat, 'Arrow','on', 'CData',CList);
BCC = BCC.draw();
BCC.tickState('on')
BCC.setFont('Visible','off')
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
if dataMat(i,j) > 0
BCC.setChordMN(i,j, 'FaceAlpha',.6)
end
end
end
text(-1.2,1.2, 'c', 'FontName','Times New Roman', 'FontSize',35)
% =========================================================================
ax4 = axes('Parent',gcf, 'Position',[1/2,0,1/2,1/2]);
ax4.XColor = 'none'; ax4.YColor = 'none';
ax4.XLim = [-1,1]; ax4.YLim = [-1,1];
hold on
NameList = {'Food supply', 'Biodiversity', 'Water quality regulation', ...
'Air quality regulation', 'Erosion control', 'Carbon storage', ...
'Water retention', 'Recreation', 'Soil quality regulation'};
patchHdl = [];
for i = 1:size(dataMat, 2)
patchHdl(i) = fill([10,11,12],[10,13,13], CList(i,:), 'EdgeColor',[0,0,0]);
end
lgdHdl = legend(patchHdl, NameList, 'Location','best', 'FontSize',14, 'FontName','Cambria', 'Box','off');
lgdHdl.Position = [.625,.11,.255,.27];
lgdHdl.ItemTokenSize = [18,8];
demo 16
dataMat = rand([15,15]);
dataMat(dataMat > .2) = 0;
CList = [ 75,146,241; 252,180, 65; 224, 64, 10; 5,100,146; 191,191,191;
26, 59,105; 255,227,130; 18,156,221; 202,107, 75; 0, 92,219;
243,210,136; 80, 99,129; 241,185,168; 224,131, 10; 120,147,190]./255;
CListC = [54,69,92]./255;
CList = CList.*.6 + CListC.*.4;
figure('Units','normalized', 'Position',[.02,.05,.6,.85])
BCC = biChordChart(dataMat, 'Arrow','on', 'CData',CList);
BCC = BCC.draw();
% 添加刻度
BCC.tickState('on')
% 修改字体，字号及颜色
BCC.setFont('FontName','Cambria', 'FontSize',17, 'Color',[0,0,0])
% 修改弦颜色(Modify chord color)
for i = 1:size(dataMat, 1)
for j = 1:size(dataMat, 2)
if dataMat(i,j) > 0
BCC.setChordMN(i,j, 'FaceColor',CListC ,'FaceAlpha',.07)
end
end
end
[~, N] = max(sum(dataMat > 0, 2));
for j = 1:size(dataMat, 2)
BCC.setChordMN(N,j, 'FaceColor',CList(N,:) ,'FaceAlpha',.6)
end
You need to download following tools:
- - Chord chart: [chord chart](https://www.mathworks.com/matlabcentral/fileexchange/116550-chord-chart)
- - Directed graph chord chart: [digraph chord chart]:(https://www.mathworks.com/matlabcentral/fileexchange/121043-digraph-chord-chart)
I wanted to turn a Markdown nested list of text labels:
- A
- B
- C
- D
- G
- H
- E
- F
- Q
into a directed graph, like this:
Here is my blog post with some related tips for doing this, including text I/O, text processing with patterns, and directed graph operations and visualization.
The topic recently came up in a MATLAB Central Answers forum thread, where community members discussed how to programmatically control when the end user can close a custom app. Imagine you need to prevent app closure during a critical process but want to allow the end user to close the app afterwards. This article will guide you through the steps to add this behavior to your app.
A demo is attached containing an app with a state button that, when enabled, disables the ability to close the app.
Steps
1. Add a property that stores the state of the closure as a scalar logical value. In this example, I named the property closeEnabled. The default value in this example is true, meaning that closing is enabled. -- How to add a property to an app in app designer
properties (Access = private)
closeEnabled = true % Flag that controls ability to close app
end
2. Add a CloseRequest function to the app figure. This function is called any time there is an attempt to close the app. Within the CloseRequest function, add a condition that deletes the app when closure is enabled. -- How to add a CloseRequest function to an app figure in app designer
function UIFigureCloseRequest(app, event)
if app.closeEnabled
delete(app)
end
3. Toggle the value of the closeEnabled property as needed in your code. Imagine you have a "Process" button that initiates a process where it is crucial for the app to remain open. Set the closeEnabled flag to false (closure is disabled) at the beginning of the button's callback function and then set it to true at the end (closure is enabled).
function ProcessButtonPress(app, event)
app.closeEnabled = false;
% MY PROCESS CODE
app.closeEnabled = true;
end
Handling Errors: There is one trap to keep in mind in the example above. What if something in the callback function breaks before the app.closeEnabled is returned to true? That leaves the app in a bad state where closure is blocked. A pro move would be to use a cleanupObj to manage returning the property to true. In the example below, the task to return the closeEnabled property to true is managed by the cleanup object, which will execute that command when execution is terminated in the ProcessButtonPress function—whether execution was terminated by error or by gracefully exiting the function.
function ProcessButtonPress(app, event)
app.closeEnabled = false;
cleanupClosure = onCleanup(@()set(app,'closeEnabled',true));
% MY CODE
end
Force Closure: If the CloseRequest function is preventing an app from closing, here are a couple of ways to force a closure.
- If you have the app's handle, use delete(app) or close(app,'force'). This will also work on the app's figure handle.
- If you do not have the app's handle, you can use close('all','force') to close all figures or use findall(groot,'type','figure') to find the app's figure handle.
Me: If you have parallel code and you apply this trick that only requires changing one line then it might go faster.
Reddit user: I did and it made my code 3x faster
Not bad for just one line of code!
Which makes me wonder. Could it make your MATLAB program go faster too? If you have some MATLAB code that makes use of parallel constructs like parfor or parfeval then start up your parallel pool like this
parpool("Threads")
before running your program.
The worst that will happen is you get an error message and you'll send us a bug report....or maybe it doesn't speed up much at all....
....or maybe you'll be like the Reddit user and get 3x speed-up for 10 seconds work. It must be worth a try...after all, you're using parallel computing to make your code faster right? May as well go all the way.
In an artificial benchmark I tried, I got 10x speedup! More details in my recent blog post: Parallel computing in MATLAB: Have you tried ThreadPools yet? » The MATLAB Blog - MATLAB & Simulink
Give it a try and let me know how you get on.
What is a rough number? What can they be used for? Today I'll take you down a journey into the land of prime numbers (in MATLAB). But remember that a journey is not always about your destination, but about what you learn along the way. And so, while this will be all about primes, and specifically large primes, before we get there we need some background. That will start with rough numbers.
Rough numbers are what I would describe as wannabe primes. Almost primes, and even sometimes prime, but often not prime. They could've been prime, but may not quite make it to the top. (If you are thinking of Marlon Brando here, telling us he "could've been a contender", you are on the right track.)
Mathematically, we could call a number k-rough if it is evenly divisible by no prime smaller than k. (Some authors will use the term k-rough to denote a number where the smallest prime factor is GREATER than k. The difference here is a minor one, and inconsequential for my purposes.) And there are also smooth numbers, numerical antagonists to the rough ones, those numbers with only small prime factors. They are not relevant to the topic today, even though smooth numbers are terribly valuable tools in mathematics. Please forward my apologies to the smooth numbers.
Have you seen rough numbers in use before? Probably so, at least if you ever learned about the sieve of Eratosthenes for prime numbers, though probably the concept of roughness was never explicitly discussed at the time. The sieve is simple. Suppose you wanted a list of all primes less than 100? (Without using the primes function itself.)
% simple sieve of Eratosthenes
Nmax = 100;
N = true(1,Nmax); % A boolean vector which when done, will indicate primes
N(1) = false; % 1 is not a prime by definition
nextP = find(N,1,'first'); % the first prime is 2
while nextP <= sqrt(Nmax)
% flag multiples of nextP as not prime
N(nextP*nextP:nextP:end) = false;
% find the first element after nextP that remains true
nextP = nextP + find(N(nextP+1:end),1,'first');
end
primeList = find(N)
Indeed, that is the set of all 25 primes not exceeding 100. If you think about how the sieve worked, it first found 2 is prime. Then it discarded all integer multiples of 2. The first element after 2 that remains as true is 3. 3 is of course the second prime. At each pass through the loop, the true elements that remain correspond to numbers which are becoming more and more rough. By the time we have eliminated all multiples of 2, 3, 5, and finally 7, everything else that remains below 100 must be prime! The next prime on the list we would find is 11, but we have already removed all multiples of 11 that do not exceed 100, since 11^2=121. For example, 77 is 11*7, but we already removed it, because 77 is a multiple of 7.
Such a simple sieve to find primes is great for small primes. However is not remotely useful in terms of finding primes with many thousands or even millions of decimal digits. And that is where I want to go, eventually. So how might we use roughness in a useful way? You can think of roughness as a way to increase the relative density of primes. That is, all primes are rough numbers. In fact, they are maximally rough. But not all rough numbers are primes. We might think of roughness as a necessary, but not sufficient condition to be prime.
How many primes lie in the interval [1e6,2e6]?
numel(primes(2e6)) - numel(primes(1e6))
There are 70435 primes greater than 1e6, but less than 2e6. Given there are 1 million natural numbers in that set, roughly 7% of those numbers were prime. Next, how many 100-rough numbers lie in that same interval?
N = (1e6:2e6)';
roughInd = all(mod(N,primes(100)) > 0,2);
sum(roughInd)
That is, there are 120571 100-rough numbers in that interval, but all those 70435 primes form a subset of the 100-rough numbers. What does this tell us? Of the 1 million numbers in that interval, approximately 12% of them were 100-rough, but 58% of the rough set were prime.
The point being, if we can efficiently identify a number as being rough, then we can substantially increase the chance it is also prime. Roughness in this sense is a prime densifier. (Is that even a word? It is now.) If we can reduce the number of times we need to perform an explicit isprime test, that will gain greatly because a direct test for primality is often quite costly in CPU time, at least on really large numbers.
In my next post, I'll show some ways we can employ rough numbers to look for some large primes.
The GCD approach to identify rough numbers is a terribly useful one, well worth remembering. But at some point, I expect someone to notice that all work done with these massively large symbolic numbers uses only one of the cores on your computer. And, having spent so much money on those extra cores in your CPU, surely we can find a way to use them all? The problem is, computations done on symbolic integers never use more than 1 core. (Sad, unhappy face.)
In order to use all of the power available to your computer using MATLAB, you need to work in double precision, or perhaps int64 or uint64. To do that, I'll next search for primes among the family 3^n+4. In fact, they seem pretty common, at least if we look at the first few such examples.
F = @(n) sym(3).^n + 4;
F(0:16)
ans =
[5, 7, 13, 31, 85, 247, 733, 2191, 6565, 19687, 59053, 177151, 531445, 1594327, 4782973, 14348911, 43046725]
isprime(F(0:16))
ans =
1×17 logical array
1 1 1 1 0 0 1 0 0 1 1 0 0 0 0 0 0
Of the first 11 members of that sequence, 7 of them were prime. Naturally, primes will become less frequent in this sequence as we look further out. The members of this family grow rapidly in size. F(10000) has 4771 decimal digits, and F(100000) has 47712 decimal digits. We certainly don't want to directly test every member of that sequence for primality. However, what I will call a partial or incomplete sieve can greatly decrease the work needed.
Consider there are roughly 5.7 million primes less than 1e8.
numel(primes(1e8))
ans =
5761455
F(17) is the first member of our sequence that exceeds 1e8. So we can start there, since we already know the small-ish primes in this sequence.
roughlim = 1e8;
primes1e8 = primes(roughlim);
primes1e8([1 2]) = []; % F(n) is never divisible by 2 or 3
F_17 = double(F(17));
Fremainders = mod(F_17,primes1e8);
nmax = 100000;
FnIsRough = false(1,nmax);
for n = 17:nmax
if all(Fremainders)
FnIsRough(n) = true;
end
% update the remainders for the next term in the sequence
% This uses the recursion: F(n+1) = 3*F(n) - 8
Fremainders = mod(Fremainders*3 - 8,primes1e8);
end
sum(FnIsRough)
ans =
6876
These will be effectively trial divides, even though we use mod for the purpose. The result is 6876 1e8-rough numbers, far less than that total set of 99984 values for n. One thing of great importance is to recognize this sequence of tests will use an approximately constant time per test regardless of the size of the numbers because each test works off the remainders from the previous one. And that works as long as we can update those remainders in some simple, direct, and efficient fashion. All that matters is the size of the set of primes to test against. Remember, the beauty of this scheme is that while I did what are implicitly trial divides against 5.76 million primes at each step, ALL of the work was done in double precision. That means I used all 8 of the cores on my computer, pushing them as hard as I could. I never had to go into the realm of big integer arithmetic to identify the rough members in that sequence, and by staying in the realm of doubles, MATLAB will automatically use all the cores you have available.
The first 10 values of n (where n is at least 17), such that F(n) is 1e8-rough were
FnIsRough = find(FnIsRough);
FnIsRough(1:10)
ans =
22 30 42 57 87 94 166 174 195 198
How well does the roughness test do to eliminate composite members of this sequence?
isprime(F(FnIsRough(1:10)))
ans =
1×10 logical array
1 1 1 1 1 0 0 1 1 1
As you can see, 8 of those first few 1e8-rough members were actually prime, so only 2 of those eventual isprime tests were effectively wasted. That means the roughness test was quite useful indeed as an efficient but relatively weak pre-test for possible primality. More importantly it is a way to quickly eliminate those values which can be known to be composite.
You can apply a similar set of tests on many families of numbers. For example, repunit primes are a great case. A rep-digit number is any number composed of a sequence of only a single digit, like 11, 777, and 9999999999999.
However, you should understand that only rep-digit numbers composed of entirely ones can ever be prime. Naturally, any number composed entirely of the digit D, will always be divisible by the single digit number D, and so only rep-unit numbers can be prime. Repunit numbers are a subset of the rep-digit family, so numbers composed only of a string of ones. 11 is the first such repunit prime. We can write them in MATLAB as a simple expression:
RU = @(N) (sym(10).^N - 1)/9;
RU(N) is a number composed only of the digit 1, with N decimal digits. This family also follows a recurrence relation, and so we could use a similar scheme as was used to find rough members of the set 3^N-4.
RU(N+1) == 10*RU(N) + 1
However, repunit numbers are rarely prime. Looking out as far as 500 digit repunit numbers, we would see primes are pretty scarce in this specific family.
find(isprime(RU(1:500)))
ans =
2 19 23 317
There are of course good reasons why repunit numbers are rarely prime. One of them is they can only ever be prime when the number of digits is also prime. This is easy to show, as you can always factor any repunit number with a composite number of digits in a simple way:
1111 (4 digits) = 11*101
111111111 (9 digits) = 111*1001001
Finally, I'll mention that Mersenne primes are indeed another example of repunit primes, when expressed in base 2. A fun fact: a Mersenne number of the form 2^n-1, when n is prime, can only have prime factors of the form 1+2*k*n. Even the Mersenne number itself will be of the same general form. And remember that a Mersenne number M(n) can only ever be prime when n is itself prime. Try it! For example, 11 is prime.
Mn = @(n) sym(2).^n - 1;
Mn(11)
ans =
2047
Note that 2047 = 1 + 186*11. But M(11) is not itself prime.
factor(Mn(11))
ans =
[23, 89]
Looking carefully at both of those factors, we see that 23 == 1+2*11, and 89 = 1+8*11.
How does this help us? Perhaps you may see where this is going. The largest known Mersenne prime at this date is Mn(136279841). This is one seriously massive prime, containing 41,024,320 decimal digits. I have no plans to directly test numbers of that size for primality here, at least not with my current computing capacity. Regardless, even at that realm of immensity, we can still do something.
If the largest known Mersenne prime comes from n=136279841, then the next such prime must have a larger prime exponent. What are the next few primes that exceed 136279841?
np = NaN(1,11); np(1) = 136279841;
for i = 1:10
np(i+1) = nextprime(np(i)+1);
end
np(1) = [];
np
np =
Columns 1 through 8
136279879 136279901 136279919 136279933 136279967 136279981 136279987 136280003
Columns 9 through 10
136280009 136280051
The next 10 candidates for Mersenne primality lie in the set Mn(np), though it is unlikely that any of those Mersenne numbers will be prime. But ... is it possible that any of them may form the next Mersenne prime? At the very least, we can exclude a few of them.
for i = 1:10
2*find(powermod(sym(2),np(i),1+2*(1:50000)*np(i))==1)
end
ans =
18 40 64
ans =
1×0 empty double row vector
ans =
2
ans =
1×0 empty double row vector
ans =
1×0 empty double row vector
ans =
1×0 empty double row vector
ans =
1×0 empty double row vector
ans =
1×0 empty double row vector
ans =
1×0 empty double row vector
ans =
2
Even with this quick test which took only a few seconds to run on my computer, we see that 3 of those Mersenne numbers are clearly not prime. In fact, we already know three of the factors of M(136279879), as 1+[18,40,64]*136279879.
You might ask, when is the MOD style test, using a large scale test for roughness against many thousands or millions of small primes, when is it better than the use of GCD? The answer here is clear. Use the large scale mod test when you can easily move from one member of the family to the next, typically using a linear recurrence. Simple such examples of this are:
1. Repunit numbers
General form: R(n) = (10^n-1)/9
Recurrence: R(n+1) = 10*R(n) + 1, R(0) = 1, R(1) = 11
2. Fibonacci numbers.
Recurrence: F(n+1) = F(n) + F(n-1), F(0) = 0, F(1) = 1
3. Mersenne numbers.
General form: M(n) = 2^n - 1
Recurrence: M(n+1) = 2*M(n) + 1
4. Cullen numbers, https://en.wikipedia.org/wiki/Cullen_number
General form: C(n) = n*2^n + 1
Recurrence: C(n+1) = 4*C(n) + 4*C(n-1) + 1
5. Hampshire numbers: (My own choice of name)
General form: H(n,b) = (n+1)*b^n - 1
Recurrence: H(n+1,b) = 2*b*H(n-1,b) - b^2*H(n-2,b) - (b-1)^2, H(0,b) = 0, H(1,b) = 2*b-1
6. Tin numbers, so named because Sn is the atomic symbol for tin.
General form: S(n) = 2*n*F(n) + 1, where F(n) is the nth Fibonacci number.
Recurrence: S(n) = S(n-5) + S(n-4) - 3*S(n-3) - S(n-2) +3*S(n-1);
To wrap thing up, I hope you have enjoyed this beginning of a journey into large primes and non-primes. I've shown a few ways we can use roughness, first in a constructive way to identify numbers which may harbor primes in a greater density than would otherwise be expected. Next, using GCD in a very pretty way, and finally by use of MOD and the full power of MATLAB to test elements of a sequence of numbers for potential primality.
My next post will delve into the world of Fermat and his little theorem, showing how it can be used as a stronger test for primality (though not perfect.)
Yes, some readers might now argue that I used roughness in a crazy way in my last post, in my approach to finding a large twin prime pair. That is, I deliberately constructed a family of integers that were known to be a-priori rough. But, suppose I gave you some large, rather arbitrarily constructed number, and asked you to tell me if it is prime? For example, to pull a number out of my hat, consider
P = sym(2)^122397 + 65;
floor(vpa(log10(P) + 1))
36846 decimal digits is pretty large. And in fact, large enough that sym/isprime in R2024b will literally choke on it. But is it prime? Can we efficiently learn if it is at least not prime?
A nice way to learn the roughness of even a very large number like this is to use GCD.
gcd(P,prod(sym(primes(10000))))
If the greatest common divisor between P and prod(sym(primes(10000))) is 1, then P is NOT divisible by any small prime from that set, since they have no common divisors. And so we can learn that P is indeed fairly rough, 10000-rough in fact. That means P is more likely to be prime than most other large integers in that domain.
gcd(P,prod(sym(primes(100000))))
However, this rather efficiently tells us that in fact, P is not prime, as it has a common factor with some integer greater than 1, and less then 1e5.
I suppose you might think this is nothing different from doing trial divides, or using the mod function. But GCD is a much faster way to solve the problem. As a test, I timed the two.
timeit(@() gcd(P,prod(sym(primes(100000)))))
timeit(@() any(mod(P,primes(100000)) == 0))
Even worse, in the first test, much if not most of that time is spent in merely computing the product of those primes.
pprod = prod(sym(primes(100000)));
timeit(@() gcd(P,pprod))
So even though pprod is itself a huge number, with over 43000 decimal digits, we can use it quite efficiently, especially if you precompute that product if you will do this often.
How might I use roughness, if my goal was to find the next larger prime beyond 2^122397? I'll look fairly deeply, looking only for 1e7-rough numbers, because these numbers are pretty seriously large. Any direct test for primality will take some serious time to perform.
pprod = prod(sym(primes(10000000)));
find(1 == gcd(sym(2)^122397 + (1:2:199),pprod))*2 - 1
2^122397 plus any one of those numbers is known to be 1e7-rough, and therefore very possibly prime. A direct test at this point would surely take hours and I don't want to wait that long. So I'll back off just a little to identify the next prime that follows 2^10000. Even that will take some CPU time.
What is the next prime that follows 2^10000? In this case, the number has a little over 3000 decimal digits. But, even with pprod set at the product of primes less than 1e7, only a few seconds were needed to identify many numbers that are 1e7-rough.
P10000 = sym(2)^10000;
k = find(1 == gcd(P10000 + (1:2:1999),pprod))*2 - 1
k =
Columns 1 through 8
15 51 63 85 165 171 177 183
Columns 9 through 16
253 267 273 295 315 421 427 451
Columns 17 through 24
511 531 567 601 603 675 687 717
Columns 25 through 32
723 735 763 771 783 793 795 823
Columns 33 through 40
837 853 865 885 925 955 997 1005
Columns 41 through 48
1017 1023 1045 1051 1071 1075 1095 1107
Columns 49 through 56
1261 1285 1287 1305 1371 1387 1417 1497
Columns 57 through 64
1507 1581 1591 1593 1681 1683 1705 1771
Columns 65 through 69
1773 1831 1837 1911 1917
Among the 1000 odd numbers immediately following 2^10000, there are exactly 69 that are 1e7-rough. Every other odd number in that sequence is now known to be composite, and even though we don't know the full factorization of those 931 composite numbers, we don't care in the context as they are not prime. I would next apply a stronger test for primality to only those few candidates which are known to be rough. Eventually after an extensive search, we would learn the next prime succeeding 2^10000 is 2^10000+13425.
In my next post, I show how to use MOD, and all the cores in your CPU to test for roughness.
How can we use roughness in an effective context to identify large primes? I can quickly think of quite a few examples where we might do so. Again, remember I will be looking for primes with not just hundreds of decimal digits, or even only a few thousand digits. The eventual target is higher than that. Forget about targets for now though, as this is a journey, and what matters in this journey is what we may learn along the way.
I think the most obvious way to employ roughness is in a search for twin primes. Though not yet proven, the twin prime conjecture:
If it is true, it tells us there are infinitely many twin prime pairs. A twin prime pair is two integers with a separation of 2, such that both of them are prime. We can find quite a few of them at first, as we have {3,5}, {5,7}, {11,13}, etc. But there is only ONE pair of integers with a spacing of 1, such that both of them are prime. That is the pair {2,3}. And since primes are less and less common as we go further out, possibly there are only a finite number of twins with a spacing of exactly 2? Anyway, while I'm fairly sure the twin prime conjecture will one day be shown to be true, it can still be interesting to search for larger and larger twin prime pairs. The largest such known pair at the moment is
2996863034895*2^1290000 +/- 1
This is a pair with 388342 decimal digits. And while seriously large, it is still in range of large integers we can work with in MATLAB, though certainly not in double precision. In my own personal work on my own computer, I've done prime testing on integers (in MATLAB) with considerably more than 100,000 decimal digits.
But, again you may ask, just how does roughness help us here? In fact, this application of roughness is not new with me. You might want to read about tools like NewPGen {https://t5k.org/programs/NewPGen/} which sieves out numbers known to be composite, before any direct tests for primality are performed.
Before we even try to talk about numbers with thousands or hundreds of thousands of decimal digits, look at 6=2*3. You might observe
isprime([-1,1] + 6)
shows that both 5 and 7 are prime. This should not be a surprise, but think about what happens, about why it generated a twin prime pair. 6 is divisible by both 2 and 3, so neither 5 or 7 can possibly be divisible by either small prime as they are one more or one less than a multiple of both 2 and 3. We can try this again, pushing the limits just a bit.
isprime([-1,1] + 2*3*5)
That is again interesting. 30=2*3*5 is evenly divisible by 2, 3, and 5. The result is both 29 and 31 are prime, because adding 1 or subtracting 1 from a multiple of 2, 3, or 5 will always result in a number that is not divisible by any of those small primes. The next larger prime after 5 is 7, but it cannot be a factor of 29 or 31, since it is greater than both sqrt(29) and sqrt(31).
We have quite efficiently found another twin prime pair. Can we take this a step further? 210=2*3*5*7 is the smallest such highly composite number that is divisible by all primes up to 7. Can we use the same trick once more?
isprime([-1,1] + 2*3*5*7)
And here the trick fails, because 209=11*19 is not in fact prime. However, can we use the large twin prime trick we saw before? Consider numbers of the form [-1,1]+a*210, where a is itself some small integer?
a = 2;
isprime([-1,1] + a*2*3*5*7)
I did not need to look far, only out to a=2, because both 419 and 421 are prime. You might argue we have formed a twin prime "factory", of sorts. Next, I'll go out as far as the product of all primes not exceeding 60. This is a number with 22 decimal digits, already too large to represent as a double, or even as uint64.
prod(sym(primes(60)))
a = find(all(isprime([-1;1] + prod(sym(primes(60)))*(1:100)),1))
That easily identifies 3 such twin prime pairs, each of which has roughly 23 decimal digits, each of which have the form a*1922760350154212639070+/-1. The twin prime factory is still working well. Going further out to integers with 37 decimal digits, we can easily find two more such pairs that employ the product of all primes not exceeding 100.
prod(sym(primes(100)))
a = find(all(isprime([-1;1] + prod(sym(primes(100)))*(1:100)),1))
This is in fact an efficient way of identifying large twin prime pairs, because it chooses a massively composite number as the product of many distinct small primes. Adding or subtracting 1 from such a number will result always in a rough number, not divisible by any of the primes employed. With a little more CPU time expended, now working with numbers with over 1000 decimal digits, I will claim this next pair forms a twin prime pair, and is the smallest such pair we can generate in this way from the product of the primes not exceeding 2500.
isprime(7826*prod(sym(primes(2500))) + [-1 1])
ans =
logical
1
Unfortunately, 1000 decimal digits is at or near the limit of what the sym/isprime tool can do for us. It does beg the question, asking if there are alternatives to the sym/isprime tool, as an isProbablePrime test, usually based on Miller-Rabin is often employed. But this is gist for yet another set of posts.
Anyway, I've done a search for primes of the form
a*prod(sym(primes(10000))) +/- 1
having gone out as far as a = 600000, with no success as of yet. (My estimate is I will find a pair by the time I get near 5e6 for a.) Anyway, if others can find a better way to search for large twin primes in MATLAB, or if you know of a larger twin prime pair of this extended form, feel free to chime in.
My next post shows how to use GCD in a very nice way to identify roughness, on a large scale.
It's frustrating when a long function or script runs and prints unexpected outputs to the command window. The line producing those outputs can be difficult to find.
Run this line of code before running the script or function. Execution will pause when the line is hit and the file will open to that line. Outputs that are intentionaly displayed by functions such as disp() or fprintf() will be ignored.
dbstop if unsuppressed output
To turn this off,
dbclear if unsuppressed output
Christmas is coming, here are two dynamic Christmas tree drawing codes:
Crystal XMas Tree
function XmasTree2024_1
fig = figure('Units','normalized', 'Position',[.1,.1,.5,.8],...
'Color',[0,9,33]/255, 'UserData',40 + [60,65,75,72,0,59,64,57,74,0,63,59,57,0,1,6,45,75,61,74,28,57,76,57,1,1]);
axes('Parent',fig, 'Position',[0,-1/6,1,1+1/3], 'UserData',97 + [18,11,0,13,3,0,17,4,17],...
'XLim',[-1.5,1.5], 'YLim',[-1.5,1.5], 'ZLim',[-.2,3.8], 'DataAspectRatio', [1,1,1], 'NextPlot','add',...
'Projection','perspective', 'Color',[0,9,33]/255, 'XColor','none', 'YColor','none', 'ZColor','none')
%% Draw Christmas tree
F = [1,3,4;1,4,5;1,5,6;1,6,3;...
2,3,4;2,4,5;2,5,6;2,6,3];
dP = @(V) patch('Faces',F, 'Vertices',V, 'FaceColor',[0 71 177]./255,...
'FaceAlpha',rand(1).*0.2+0.1, 'EdgeColor',[0 71 177]./255.*0.8,...
'EdgeAlpha',0.6, 'LineWidth',0.5, 'EdgeLighting','gouraud', 'SpecularStrength',0.3);
r = .1; h = .8;
V0 = [0,0,0; 0,0,1; 0,r,h; r,0,h; 0,-r,h; -r,0,h];
% Rotation matrix
Rx = @(V, theta) V*[1 0 0; 0 cos(theta) sin(theta); 0 -sin(theta) cos(theta)];
Rz = @(V, theta) V*[cos(theta) sin(theta) 0;-sin(theta) cos(theta) 0; 0 0 1];
N = 180; Vn = zeros(N, 3); eval(char(fig.UserData))
for i = 1:N
tV = Rz(Rx(V0.*(1.2 - .8.*i./N + rand(1).*.1./i^(1/5)), pi/3.*(1 - .6.*i./N)), i.*pi/8.1 + .001.*i.^2) + [0,0,.016.*i];
dP(tV); Vn(i,:) = tV(2,:);
end
scatter3(Vn(:,1).*1.02,Vn(:,2).*1.02,Vn(:,3).*1.01, 30, 'w', 'Marker','*', 'MarkerEdgeAlpha',.5)
%% Draw Star of Bethlehem
w = .3; R = .62; r = .4; T = (1/8:1/8:(2 - 1/8)).'.*pi;
V8 = [ 0, 0, w; 0, 0,-w;
1, 0, 0; 0, 1, 0; -1, 0, 0; 0,-1,0;
R, R, 0; -R, R, 0; -R,-R, 0; R,-R,0;
cos(T).*r, sin(T).*r, T.*0];
F8 = [1,3,25; 1,3,11; 2,3,25; 2,3,11; 1,7,11; 1,7,13; 2,7,11; 2,7,13;
1,4,13; 1,4,15; 2,4,13; 2,4,15; 1,8,15; 1,8,17; 2,8,15; 2,8,17;
1,5,17; 1,5,19; 2,5,17; 2,5,19; 1,9,19; 1,9,21; 2,9,19; 2,9,21;
1,6,21; 1,6,23; 2,6,21; 2,6,23; 1,10,23; 1,10,25; 2,10,23; 2,10,25];
V8 = Rx(V8.*.3, pi/2) + [0,0,3.5];
patch('Faces',F8, 'Vertices',V8, 'FaceColor',[255,223,153]./255,...
'EdgeColor',[255,223,153]./255, 'FaceAlpha', .2)
%% Draw snow
sXYZ = rand(200,3).*[4,4,5] - [2,2,0];
sHdl1 = plot3(sXYZ(1:90,1),sXYZ(1:90,2),sXYZ(1:90,3), '*', 'Color',[.8,.8,.8]);
sHdl2 = plot3(sXYZ(91:200,1),sXYZ(91:200,2),sXYZ(91:200,3), '.', 'Color',[.6,.6,.6]);
annotation(fig,'textbox',[0,.05,1,.09], 'Color',[1 1 1], 'String','Merry Christmas Matlaber',...
'HorizontalAlignment','center', 'FontWeight','bold', 'FontSize',48,...
'FontName','Times New Roman', 'FontAngle','italic', 'FitBoxToText','off','EdgeColor','none');
% Rotate the Christmas tree and let the snow fall
for i=1:1e8
sXYZ(:,3) = sXYZ(:,3) - [.05.*ones(90,1); .06.*ones(110,1)];
sXYZ(sXYZ(:,3)<0, 3) = sXYZ(sXYZ(:,3) < 0, 3) + 5;
sHdl1.ZData = sXYZ(1:90,3); sHdl2.ZData = sXYZ(91:200,3);
view([i,30]); drawnow; pause(.05)
end
end
Curved XMas Tree
function XmasTree2024_2
fig = figure('Units','normalized', 'Position',[.1,.1,.5,.8],...
'Color',[0,9,33]/255, 'UserData',40 + [60,65,75,72,0,59,64,57,74,0,63,59,57,0,1,6,45,75,61,74,28,57,76,57,1,1]);
axes('Parent',fig, 'Position',[0,-1/6,1,1+1/3], 'UserData',97 + [18,11,0,13,3,0,17,4,17],...
'XLim',[-6,6], 'YLim',[-6,6], 'ZLim',[-16, 1], 'DataAspectRatio', [1,1,1], 'NextPlot','add',...
'Projection','perspective', 'Color',[0,9,33]/255, 'XColor','none', 'YColor','none', 'ZColor','none')
%% Draw Christmas tree
[X,T] = meshgrid(.4:.1:1, 0:pi/50:2*pi);
XM = 1 + sin(8.*T).*.05;
X = X.*XM; R = X.^(3).*(.5 + sin(8.*T).*.02);
dF = @(R, T, X) surf(R.*cos(T), R.*sin(T), -X, 'EdgeColor',[20,107,58]./255,...
'FaceColor', [20,107,58]./255, 'FaceAlpha',.2, 'LineWidth',1);
CList = [254,103,110; 255,191,115; 57,120,164]./255;
for i = 1:5
tR = R.*(2 + i); tT = T+i; tX = X.*(2 + i) + i;
SFHdl = dF(tR, tT, tX);
[~, ind] = sort(SFHdl.ZData(:)); ind = ind(1:8);
C = CList(randi([1,size(CList,1)], [8,1]), :);
scatter3(tR(ind).*cos(tT(ind)), tR(ind).*sin(tT(ind)), -tX(ind), 120, 'filled',...
'CData', C, 'MarkerEdgeColor','none', 'MarkerFaceAlpha',.3)
scatter3(tR(ind).*cos(tT(ind)), tR(ind).*sin(tT(ind)), -tX(ind), 60, 'filled', 'CData', C)
end
%% Draw Star of Bethlehem
Rx = @(V, theta) V*[1 0 0; 0 cos(theta) sin(theta); 0 -sin(theta) cos(theta)];
% Rz = @(V, theta) V*[cos(theta) sin(theta) 0;-sin(theta) cos(theta) 0; 0 0 1];
w = .3; R = .62; r = .4; T = (1/8:1/8:(2 - 1/8)).'.*pi;
V8 = [ 0, 0, w; 0, 0,-w;
1, 0, 0; 0, 1, 0; -1, 0, 0; 0,-1,0;
R, R, 0; -R, R, 0; -R,-R, 0; R,-R,0;
cos(T).*r, sin(T).*r, T.*0];
F8 = [1,3,25; 1,3,11; 2,3,25; 2,3,11; 1,7,11; 1,7,13; 2,7,11; 2,7,13;
1,4,13; 1,4,15; 2,4,13; 2,4,15; 1,8,15; 1,8,17; 2,8,15; 2,8,17;
1,5,17; 1,5,19; 2,5,17; 2,5,19; 1,9,19; 1,9,21; 2,9,19; 2,9,21;
1,6,21; 1,6,23; 2,6,21; 2,6,23; 1,10,23; 1,10,25; 2,10,23; 2,10,25];
V8 = Rx(V8.*.8, pi/2) + [0,0,-1.3];
patch('Faces',F8, 'Vertices',V8, 'FaceColor',[255,223,153]./255,...
'EdgeColor',[255,223,153]./255, 'FaceAlpha', .2)
annotation(fig,'textbox',[0,.05,1,.09], 'Color',[1 1 1], 'String','Merry Christmas Matlaber',...
'HorizontalAlignment','center', 'FontWeight','bold', 'FontSize',48,...
'FontName','Times New Roman', 'FontAngle','italic', 'FitBoxToText','off','EdgeColor','none');
%% Draw snow
sXYZ = rand(200,3).*[12,12,17] - [6,6,16];
sHdl1 = plot3(sXYZ(1:90,1),sXYZ(1:90,2),sXYZ(1:90,3), '*', 'Color',[.8,.8,.8]);
sHdl2 = plot3(sXYZ(91:200,1),sXYZ(91:200,2),sXYZ(91:200,3), '.', 'Color',[.6,.6,.6]);
for i=1:1e8
sXYZ(:,3) = sXYZ(:,3) - [.1.*ones(90,1); .12.*ones(110,1)];
sXYZ(sXYZ(:,3)<-16, 3) = sXYZ(sXYZ(:,3) < -16, 3) + 17.5;
sHdl1.ZData = sXYZ(1:90,3); sHdl2.ZData = sXYZ(91:200,3);
view([i,30]); drawnow; pause(.05)
end
end
I wish all MATLABers a Merry Christmas in advance!
About Tips & Tricks
Looking to improve your skills and efficiency? Reach your coding goals quickly and more easily with how-tos, tutorials, and shortcuts posted by our staff and community power users.Do you have your own tips and tricks to share? Help other community users with your insights and experience.
