/matlabcentral/discussions/channelsConversaciones en Channels2024-07-21T02:51:11Ztag:es.mathworks.com,2005:Topic/8667112024-06-18T21:19:39Z2024-07-21T02:51:11ZAsk Me Anything about image analysis or the Mathworks community<p>Hello, everyone! I’m Mark Hayworth, but you might know me better in the community as Image Analyst. I've been using MATLAB since 2006 (18 years). My background spans a rich career as a former senior scientist and inventor at The Procter & Gamble Company (HQ in Cincinnati). I hold both master’s & Ph.D. degrees in optical sciences from the College of Optical Sciences at the University of Arizona, specializing in imaging, image processing, and image analysis. I have 40+ years of military, academic, and industrial experience with image analysis programming and algorithm development. I have experience designing custom light booths and other imaging systems. I also work with color and monochrome imaging, video analysis, thermal, ultraviolet, hyperspectral, CT, MRI, radiography, profilometry, microscopy, NIR, and Raman spectroscopy, etc. on a huge variety of subjects.
I'm thrilled to participate in MATLAB Central's Ask Me Anything (AMA) session, a fantastic platform for knowledge sharing and community engagement. Following Adam Danz’s insightful AMA on staff contributors in the Answers forum, I’d like to discuss topics in the area of image analysis and processing. I invite you to ask me anything related to this field, whether you're seeking recommendations on tools, looking for tips and tricks, my background, or career development advice. Additionally, I'm more than willing to share insights from my experiences in the MATLAB Answers community, File Exchange, and my role as a member of the Community Advisory Board. If you have questions related to your specific images or your custom MATLAB code though, I'll invite you to ask those in the Answers forum. It's a more appropriate forum for those kinds of questions, plus you can get the benefit of other experts offering their solutions in addition to me.
For the coming weeks, I'll be here to engage with your questions and help shed light on any topics you're curious about.</p>Image Analysthttps://es.mathworks.com/matlabcentral/profile/authors/1343420tag:es.mathworks.com,2005:Topic/8441012024-02-02T17:41:51Z2024-07-15T15:17:31ZRead this before posting<p>Hello and a warm welcome to all! We're thrilled to have you visit our community. MATLAB Central is a place for learning, sharing, and connecting with others who share your passion for MATLAB and Simulink. To ensure you have the best experience, here are some tips to get you started:
Read the Community Guidelines: Understanding our community standards is crucial. Please take a moment to familiarize yourself with them. Keep in mind that posts not adhering to these guidelines may be flagged by moderators or other community members.
Ask Technical Questions at MATLAB Answers: If you have questions related to MathWorks products, head over to MATLAB Answers (new question form - Ask the community). It's the go-to spot for technical inquiries, with responses often provided within an hour, depending on the complexity of the question and volunteer availability. To increase your chances of a speedy reply, check out our tips on how to craft a good question (link to post on asking good questions).
Choosing the Right Channel: We offer a variety of discussion channels tailored to different contexts. Select the one that best fits your post. If you're unsure, the General channel is always a safe bet. If you feel there's a need for a new channel, we encourage you to suggest it in the Ideas channel.
Reporting Issues: If you encounter posts that violate our guidelines, please use the 🚩Flag/Report feature (found in the 3-dot menu) to bring them to our attention.
Quality Control: We strive to maintain a high standard of discussion. Accounts that post spam or too much nonsense may be subject to moderation, which can include temporary suspensions or permanent bans.
Share Your Ideas: Your feedback is invaluable. If you have suggestions on how we can improve the community or MathWorks products, the Ideas channel is the perfect place to voice your thoughts.
Enjoy yourself and have fun! We're committed to fostering a supportive and educational environment. Dive into discussions, share your expertise, and grow your knowledge. We're excited to see what you'll contribute to the community!</p>Davidhttps://es.mathworks.com/matlabcentral/profile/authors/4480925tag:es.mathworks.com,2005:Topic/8648412024-06-03T21:22:40Z2024-06-03T21:22:40ZWelcome to the Cody Discussion Channel! Please Read Before Posting<p>Hello and a warm welcome to everyone! We're excited to have you in the Cody Discussion Channel. To ensure the best possible experience for everyone, it's important to understand the types of content that are most suitable for this channel.</p><p>Content that belongs in the Cody Discussion Channel:
Tips & tricks: Discuss strategies for solving Cody problems that you've found effective.
Ideas or suggestions for improvement: Have thoughts on how to make Cody better? We'd love to hear them.
Issues: Encountering difficulties or bugs with Cody? Let us know so we can address them.
Requests for guidance: Stuck on a Cody problem? Ask for advice or hints, but make sure to show your efforts in attempting to solve the problem first.
General discussions: Anything else related to Cody that doesn't fit into the above categories.</p><p>Content that does not belong in the Cody Discussion Channel:
Comments on specific Cody problems: Examples include unclear problem descriptions or incorrect testing suites.
Comments on specific Cody solutions: For example, you find a solution creative or helpful.
Please direct such comments to the Comments section on the problem or solution page itself.</p><p>We hope the Cody discussion channel becomes a vibrant space for sharing expertise, learning new skills, and connecting with others.</p>Chen Linhttps://es.mathworks.com/matlabcentral/profile/authors/6682740tag:es.mathworks.com,2005:Topic/8314622023-11-07T18:43:52Z2023-11-09T13:55:07Z🚀 The MATLAB AI Chat Playground Has Launched<p>The MATLAB AI Chat Playground is open to everyone!
Check it out here on the community: https://www.mathworks.com/matlabcentral/playground</p><p>I just published a blog post announcing the release.</p>Hans Scharlerhttps://es.mathworks.com/matlabcentral/profile/authors/5863695tag:es.mathworks.com,2005:Topic/8705862024-07-20T17:37:30Z2024-07-20T17:37:30ZCode help<p>Hello, I am a beginner in MATLAB.I am trying to solve this problem, can someone help.
A space rocket has 3 stages:
stage 1, s1;
stage 2, s2;
stage 3, s3.
If s1 burns 3 x as long as s2 which burns 2 x as long as s3 then how long did s3 burn if the total burn time was tt minutes? How long did s2 burn? s1?
This is the code I wrote:
function s = rocketburntime(totaltime,r1,r2)
t1=totaltime/1.5
t2=totaltime/4.5
t3=totaltime/9 s=[t1 t2 t3]
end
But on MATLAB cody for input (32,4,3) , its giving error</p>Zuha Altafhttps://es.mathworks.com/matlabcentral/profile/authors/15222140tag:es.mathworks.com,2005:Topic/8591662024-05-07T00:14:14Z2024-07-20T16:47:10ZHow many code files are typically in your MATLAB Projects (scripts, functions, classes, tests)?goc3https://es.mathworks.com/matlabcentral/profile/authors/5349647tag:es.mathworks.com,2005:Topic/8660312024-06-13T19:18:16Z2024-07-19T15:22:05Z+1 for Backwards Loops!<p>Base case:
Suppose you need to do a computation many times. We are going to assume that this computation cannot be vectorized. The simplest case is to use a for loop:
number_of_elements = 1e6;
test_fcn = @(x) sqrt(x) / x;
tic
for i = 1:number_of_elements
x(i) = test_fcn(i);
end
t_forward = toc;
disp(t_forward + " seconds")
Preallocation:
This can easily be sped up by preallocating the variable that houses results:
tic
x = zeros(number_of_elements, 1);
for i = 1:number_of_elements
x(i) = test_fcn(i);
end
t_forward_prealloc = toc;
disp(t_forward_prealloc + " seconds")
In this example, preallocation speeds up the loop by a factor of about three to four (running in R2024a). Comment below if you get dramatically different results.
disp(sprintf("%.1f", t_forward / t_forward_prealloc))
Run it in reverse:
Is there a way to skip the explicit preallocation and still be fast? Indeed, there is.
clear x
tic
for i = number_of_elements:-1:1
x(i) = test_fcn(i);
end
t_backward = toc;
disp(t_backward + " seconds")
By running the loop backwards, the preallocation is implicitly performed during the first iteration and the loop runs in about the same time (within statistical noise):
disp(sprintf("%.2f", t_forward_prealloc / t_backward))
Do you get similar results when running this code? Let us know your thoughts in the comments below.
Beneficial side effect:
Have you ever had to use a for loop to delete elements from a vector? If so, keeping track of index offsets can be tricky, as deleting any element shifts all those that come after. By running the for loop in reverse, you don't need to worry about index offsets while deleting elements.</p>goc3https://es.mathworks.com/matlabcentral/profile/authors/5349647tag:es.mathworks.com,2005:Topic/8473212024-02-28T20:25:39Z2024-07-19T10:25:29ZSet colormap limits when creating m colors<p>Several of the colormaps are great for a 256 color surface plot, but aren't well optimized for extracting m colors for plotting several independent lines. The issue is that many colormaps have start/end colors that are too similar or are suboptimal colors for lines. There are certainly many workarounds for this, but it would be a great quality of life to adjust that directly when calling this.
Example:
x = linspace(0,2*pi,101)';
y = [1:6].*cos(x);
figure; plot(x,y,'LineWidth',2); grid on; axis tight;
And now if I wanted to color these lines, I could use something like turbo(6) or gray(6) and then apply it using colororder.
colororder(turbo(6))
But my issue is that the ends of the colormap are too similar. For other colormaps, you may get lines that are too light to be visible against the white background. There are plenty of workarounds, with my preference being to create extra colors and truncate that before using colororder.
cmap = turbo(8); cmap = cmap(2:end-1,:); % Truncate the end colors
figure; plot(x,y,'LineWidth',2); grid on; axis tight;
colororder(cmap)
I think it would be really awesome to add some name-argument input pair to these colormaps that can specify the range you want so this could even be done inside the colororder calling if desired. An example of my proposed solution would look something like this:
cmap = turbo(6,'Range',[0.1 0.8]); % Proposed idea to add functionality
Where in this scenario, the resulting colormap would be 6 equally spaced colors that range from 10% to 80% of the total color range. This would be especially nice because you could more quickly modify the range of colors, or you could set the limits regardless of whether you need to plot 3, 6, or 20 lines.</p>Emma Farnanhttps://es.mathworks.com/matlabcentral/profile/authors/4779029tag:es.mathworks.com,2005:Topic/8693332024-07-09T22:11:17Z2024-07-18T17:03:27ZMatlab Custom Font<p>I've noticed is that the highly rated fonts for coding (e.g. Fira Code, Inconsolata, etc.) seem to overlook one issue that is key for coding in Matlab. While these fonts make 0 and O, as well as the 1 and l easily distinguishable, the brackets are not. Quite often the curly bracket looks similar to the curved bracket, which can lead to mistakes when coding or reviewing code.
So I was thinking: Could Mathworks put together a team to review good programming fonts, and come up with their own custom font designed specifically and optimized for Matlab syntax?</p>Honzikhttps://es.mathworks.com/matlabcentral/profile/authors/219831tag:es.mathworks.com,2005:Topic/8483362024-03-06T17:18:42Z2024-07-16T13:33:19ZHelp Document of Chord Chart<p>Chord diagrams are very common in Python and R, but there are no related functions in MATLAB before. It is not easy to draw chord diagrams of the same quality as R language, But I created a MATLAB tool that could almost do it.
↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓
https://www.mathworks.com/matlabcentral/fileexchange/116550-chord-chart
↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑
Here is the help document:</p><p>1 Data Format
The data requirement is a numerical matrix with all values greater than or equal to 0, or a table array, or a numerical matrix and cell array for names. First, give an example of a numerical matrix:
1.1 Numerical Matrix
dataMat=randi([0,5],[5,4]);</p><p>% 绘图(draw)
CC=chordChart(dataMat);
CC=CC.draw();</p><p>Since each object is not named, it will be automatically named Rn and Cn
1.2 Numerical Matrix and Cell Array for Names
dataMat=[2 0 1 2 5 1 2;
3 5 1 4 2 0 1;
4 0 5 5 2 4 3];
colName={'G1','G2','G3','G4','G5','G6','G7'};
rowName={'S1','S2','S3'};</p><p>CC=chordChart(dataMat,'rowName',rowName,'colName',colName);
CC=CC.draw();</p><p>RowName should be the same size as the rows of the matrix
ColName should be the same size as the columns of the matrix
For this example, if the value in the second row and third column is 1, it indicates that there is an energy flow from S2 to G3, and a chord with a width of 1 is needed between these two.
1.3 Table Array
A table array in the following format is required:</p><p>2 Decorate Chord
2.1 Batch modification of chords
Batch modification of chords can be done using the setChordProp function, and all properties of the Patch object can be modified. For example, modifying the color of the string, edge color, edge line sstyle, etc.:
CC.setChordProp('EdgeColor',[.3,.3,.3],'LineStyle','--',...
'LineWidth',.1,'FaceColor',[.3,.3,.3])</p><p>2.2 Individual Modification of Chord
The individual modification of chord can be done using the setChordMN function, where the values of m and n correspond exactly to the rows and columns of the original numerical matrix. For example, changing the color of the strings flowing from S2 to G4 to red:
CC.setChordMN(2,4,'FaceColor',[1,0,0])</p><p>2.3 Color Mapping of Chords
Just use function colormap to do so:</p><p>% version 1.7.0更新
% 可使用colormap函数直接修改颜色
% Colors can be adjusted directly using the function colormap(demo4)</p><p>colormap(flipud(pink))</p><p>3 Arc Shaped Block Decoration
3.1 Batch Decoration of Arc-Shaped Blocks</p><p>use:
setSquareT_Prop
setSquareF_Prop
to modify the upper and lower blocks separately, and all attributes of the Patch object can be modified. For example, batch modify the upper blocks (change to black):
CC.setSquareT_Prop('FaceColor',[0,0,0])</p><p>3.2 Arc-Shaped Blocks Individually Decoration
use:
setSquareT_N
setSquareF_N
to modify the upper and lower blocks separately. For example, modify the second block above separately (changed to red):
CC.setSquareT_N(2,'FaceColor',[.8,0,0])</p><p>4 Font Adjustment
Use the setFont function to adjust the font, and all properties of the text object can be modified. For example, changing the font size, font, and color of the text:
CC.setFont('FontSize',25,'FontName','Cambria','Color',[0,0,.8])</p><p>5 Show and Hide Ticks
Usage:
CC.tickState('on')
% CC.tickState('off')</p><p>6 Attribute 'Sep' with Adjustable Square Spacing
If the matrix size is large, the drawing will be out of scale:
dataMat=randi([0,1],[20,10]);
CC=chordChart(dataMat);
CC=CC.draw();</p><p>% CC.tickState('on')</p><p>We can modify its Sep attribute:
dataMat=randi([0,1],[20,10]);
% use Sep to decrease space (separation)
% 使用 sep 减小空隙
CC=chordChart(dataMat,'Sep',1/120);
CC=CC.draw();</p><p>7 Modify Text Direction
dataMat=randi([0,1],[20,10]);
% use Sep to decrease space (separation)
% 使用 sep 减小空隙
CC=chordChart(dataMat,'Sep',1/120);
CC=CC.draw();</p><p>CC.tickState('on')</p><p>% version 1.7.0更新
% 函数labelRatato用来旋转标签
% The function labelRatato is used to rotate the label
CC.labelRotate('on')</p><p>8 Add Tick Labels
dataMat=[2 0 1 2 5 1 2;
3 5 1 4 2 0 1;
4 0 5 5 2 4 3];
colName={'G1','G2','G3','G4','G5','G6','G7'};
rowName={'S1','S2','S3'};</p><p>CC=chordChart(dataMat,'rowName',rowName,'colName',colName);
CC=CC.draw();
CC.setFont('FontSize',17,'FontName','Cambria')</p><p>% 显示刻度和数值
% Displays scales and numeric values
CC.tickState('on')
CC.tickLabelState('on')</p><p>% 调节标签半径
% Adjustable Label radius
CC.setLabelRadius(1.3);</p><p>% figure()
% dataMat=[2 0 1 2 5 1 2;
% 3 5 1 4 2 0 1;
% 4 0 5 5 2 4 3];
% dataMat=dataMat+rand(3,7);
% dataMat(dataMat<1)=0;
%
% CC=chordChart(dataMat,'rowName',rowName,'colName',colName);
% CC=CC.draw();
% CC.setFont('FontSize',17,'FontName','Cambria')
%
% % 显示刻度和数值
% % Displays scales and numeric values
% CC.tickState('on')
% CC.tickLabelState('on')
%
% % 调节标签半径
% % Adjustable Label radius
% CC.setLabelRadius(1.4);</p><p>9 Custom Tick Label Format
A function handle is required to input numeric output strings. The format can be set through the setTickLabelFormat function, such as Scientific notation:
dataMat=[2 0 1 2 5 1 2;
3 5 1 4 2 0 1;
4 0 5 5 2 4 3];
dataMat=dataMat+rand(3,7);
dataMat(dataMat<1)=0;
dataMat=dataMat.*1000;</p><p>CC=chordChart(dataMat);
CC=CC.draw();
CC.setFont('FontSize',17,'FontName','Cambria')</p><p>% 显示刻度和数值
% Displays scales and numeric values
CC.tickState('on')
CC.tickLabelState('on')</p><p>% 调节标签半径
% Adjustable Label radius
CC.setLabelRadius(1.4);</p><p>% 调整数值字符串格式
% Adjust numeric string format
CC.setTickLabelFormat(@(x)sprintf('%0.1e',x))</p><p>10 A Demo
rng(2)</p><p>dataMat=randi([1,7],[11,5]);
colName={'Fly','Beetle','Leaf','Soil','Waxberry'};
rowName={'Bartomella','Bradyrhizobium','Dysgomonas','Enterococcus',...
'Lactococcus','norank','others','Pseudomonas','uncultured',...
'Vibrionimonas','Wolbachia'};</p><p>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:5
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:11
CC.setSquareF_N(i,'FaceColor',CListF(i,:))
end
% 修改弦颜色(Modify chord color)
for i=1:5
for j=1:11
CC.setChordMN(j,i,'FaceColor',CListT(i,:),'FaceAlpha',.5)
end
end</p><p>CC.tickState('on')
CC.labelRotate('on')
CC.setFont('FontSize',17,'FontName','Cambria')</p><p>Hope to have your Reviews and Stars!!!
↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓
https://www.mathworks.com/matlabcentral/fileexchange/116550-chord-chart
↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑</p>Zhaoxu Liu / slandarerhttps://es.mathworks.com/matlabcentral/profile/authors/18192500tag:es.mathworks.com,2005:Topic/8688312024-07-05T17:56:04Z2024-07-14T04:24:52ZNeed help to get started with matlab<p>Hello everyone, i hope you all are in good health. i need to ask you about the help about where i should start to get indulge in matlab. I am an electrical engineer but having experience of construction field. I am new here. Please do help me. I shall be waiting forward to hear from you. I shall be grateful to you. Need recommendations and suggestions from experienced members. Thank you.</p>Muhammadhttps://es.mathworks.com/matlabcentral/profile/authors/34384838tag:es.mathworks.com,2005:Topic/8698932024-07-13T06:25:53Z2024-07-13T14:31:05ZMore easily enable CAD-like view mode<p>Something that had bothered me ever since I became an FEA analyst (2012) was the apparent inability of the "camera" in Matlab's 3D plot to function like the "cameras" in CAD/CAE packages.
For instance, load the ForearmLink.stl model that ships with the PDE Toolbox in Matlab and ParaView and try rotating the model.
clear
close all
gm = importGeometry( "ForearmLink.stl" );
pdegplot(gm)</p><p>To provide talking points, here's a YouTube video I recorded.
Things to observe:
Note that I cant seem to rotate continuously around the x-axis. It appears to only support rotations from [0, 360] as opposed to [-inf, inf]. So, for example, if I'm looking in the Y+ direction and rotate around X so that I'm looking at the Z- direction, and then want to look in the Y- direction, I can't simply keep rotating around the X axis... instead have to rotate 180 degrees around the Z axis and then around the X axis. I'm not aware of any data visualization applications (e.g., ParaView, VisIt, EnSight) or CAD/CAE tools with such an interaction.
Note that at the 50 second mark, I set a view in ParaView: looking in the [X-, Y-, Z-] direction with Y+ up. Try as I might in Matlab, I'm unable to achieve that same view perspective.
Today I discovered that if one turns on the Camera Toolbar from the View menubar, then clicks the Orbit Camera icon, then the No Principal Axis icon:</p><p>That then it acts in the manner I've long desired. Oh, and also, for the interested, it is programmatically available: https://www.mathworks.com/help/matlab/ref/cameratoolbar.html
I might humbly propose this mode either be made more discoverable, similar to the little interaction widgets that pop up in figures:</p><p>Or maybe use the middle-mouse button to temporarily use this mode (a mouse setting in, e.g., Abaqus/CAE).</p>Gregory Vernonhttps://es.mathworks.com/matlabcentral/profile/authors/14170832tag:es.mathworks.com,2005:Topic/8606162024-05-13T12:01:17Z2024-07-10T05:24:10ZWhich of the following is not a built-in string-checking function?goc3https://es.mathworks.com/matlabcentral/profile/authors/5349647tag:es.mathworks.com,2005:Topic/8690112024-07-07T20:34:22Z2024-07-09T18:48:38ZRequest for Assistance with Calculating Controller Gain Values<p>could you explain me how to calculate the gain values for different types of controllers (Conventional Sliding Mode Control, Third Order Sliding Mode Control, Variable Gain Super Twisting Algorithm.
Could you, assist me in providing a mathematical method, for example, to calculate the gains of the above-mentioned controllers?
Thank you
M. Itouchene</p>itouchenehttps://es.mathworks.com/matlabcentral/profile/authors/28540679tag:es.mathworks.com,2005:Topic/7966772023-08-20T02:57:58Z2024-07-05T19:20:23ZMy preferred vacation is to visitImage Analysthttps://es.mathworks.com/matlabcentral/profile/authors/1343420tag:es.mathworks.com,2005:Topic/8481262011-02-23T00:18:46Z2024-07-05T06:16:16ZSome Foreign Matlab forums<p>Which Matlab related forums and newsgroups do you use beside MATLAB Answers? Which languages do they use? Which advantages and unique features do they have?</p><p>Do you think that these forums complement or compete against MathWorks and its communication platform?</p><p>Actually <b>all</b> answers are accepted.</p>Janhttps://es.mathworks.com/matlabcentral/profile/authors/869888tag:es.mathworks.com,2005:Topic/8664112024-06-16T23:28:03Z2024-07-04T03:12:44ZStochastic simulation of a novel pathogen<p>We are modeling the introduction of a novel pathogen into a completely susceptible population. In the cells below, I have provided you with the Matlab code for a simple stochastic SIR model, implemented using the "GillespieSSA" function</p><p>Simulating the stochastic model 100 times for
Since is 0.4 per day, per day
% Define the parameters
beta = 0.36;
gamma = 0.4;
n_sims = 100;
tf = 100; % Time frame changed to 100</p><p>% Calculate R0
R0 = beta / gamma</p><p>% Initial state values
initial_state_values = [1000000; 1; 0; 0]; % S, I, R, cum_inc</p><p>% Define the propensities and state change matrix
a = @(state) [beta * state(1) * state(2) / 1000000, gamma * state(2)];
nu = [-1, 0; 1, -1; 0, 1; 0, 0];</p><p>% Define the Gillespie algorithm function
function [t_values, state_values] = gillespie_ssa(initial_state, a, nu, tf)
t = 0;
state = initial_state(:); % Ensure state is a column vector
t_values = t;
state_values = state';</p><pre> while t < tf
rates = a(state);
rate_sum = sum(rates);
if rate_sum == 0
break;
end</pre><pre> tau = -log(rand) / rate_sum;
t = t + tau;</pre><pre> r = rand * rate_sum;
cum_sum_rates = cumsum(rates);
reaction_index = find(cum_sum_rates >= r, 1);</pre><pre> state = state + nu(:, reaction_index);</pre><pre> % Update cumulative incidence if infection occurred
if reaction_index == 1
state(4) = state(4) + 1; % Increment cumulative incidence
end</pre><pre> t_values = [t_values; t];
state_values = [state_values; state'];
end
end</pre><p>% Function to simulate the stochastic model multiple times and plot results
function simulate_stoch_model(beta, gamma, n_sims, tf, initial_state_values, R0, plot_type)
% Define the propensities and state change matrix
a = @(state) [beta * state(1) * state(2) / 1000000, gamma * state(2)];
nu = [-1, 0; 1, -1; 0, 1; 0, 0];</p><pre> % Set random seed for reproducibility
rng(11);</pre><pre> % Initialize plot
figure;
hold on;</pre><pre> for i = 1:n_sims
[t, output] = gillespie_ssa(initial_state_values, a, nu, tf);</pre><pre> % Check if the simulation had only one step and re-run if necessary
while length(t) == 1
[t, output] = gillespie_ssa(initial_state_values, a, nu, tf);
end</pre><pre> if strcmp(plot_type, 'cumulative_incidence')
plot(t, output(:, 4), 'LineWidth', 2, 'Color', rand(1, 3));
elseif strcmp(plot_type, 'prevalence')
plot(t, output(:, 2), 'LineWidth', 2, 'Color', rand(1, 3));
end
end</pre><pre> xlabel('Time (days)');</pre><pre> if strcmp(plot_type, 'cumulative_incidence')
ylabel('Cumulative Incidence');
ylim([0 inf]);
elseif strcmp(plot_type, 'prevalence')
ylabel('Prevalence of Infection');
ylim([0 50]);
end</pre><pre> title(['Stochastic model output for R0 = ', num2str(R0)]);
subtitle([num2str(n_sims), ' simulations']);
xlim([0 tf]);
grid on;
hold off;
end</pre><p>% Simulate the model 100 times and plot cumulative incidence
simulate_stoch_model(beta, gamma, n_sims, tf, initial_state_values, R0, 'cumulative_incidence');</p><p>% Simulate the model 100 times and plot prevalence
simulate_stoch_model(beta, gamma, n_sims, tf, initial_state_values, R0, 'prevalence');</p>Athanasios Paraskevopouloshttps://es.mathworks.com/matlabcentral/profile/authors/30623616tag:es.mathworks.com,2005:Topic/8607012024-05-13T19:06:20Z2024-07-03T16:41:09ZNorthern Lights<p>Northern lights captured from this weekend at MathWorks campus ✨</p><p>Did you get a chance to see lights and take some photos?</p>Chen Linhttps://es.mathworks.com/matlabcentral/profile/authors/6682740tag:es.mathworks.com,2005:Topic/8683862024-07-02T13:01:53Z2024-07-02T13:01:53Zdocument on solving ODEs and PDEs<p>I recently wrote up a document which addresses the solution of ordinary and partial differential equations in Matlab (with some Python examples thrown in for those who are interested). For ODEs, both initial and boundary value problems are addressed. For PDEs, it addresses parabolic and elliptic equations. The emphasis is on finite difference approaches and built-in functions are discussed when available. Theory is kept to a minimum. I also provide a discussion of strategies for checking the results, because I think many students are too quick to trust their solutions. For anyone interested, the document can be found at https://blanchard.neep.wisc.edu/SolvingDifferentialEquationsWithMatlab.pdf</p>James Blanchardhttps://es.mathworks.com/matlabcentral/profile/authors/6002112tag:es.mathworks.com,2005:Topic/8576412024-04-27T21:03:01Z2024-07-01T20:54:05ZHow many spaces per tab do you prefer?goc3https://es.mathworks.com/matlabcentral/profile/authors/5349647