Emergency Research: Alleviating Spiciness

Let me introduce you to sambal, extremely spicy sauce with a rich seafood taste. Also for the quick research on dealing with spiciness skip to the bottom paragraph.

“Sambal is a hot sauce or paste typically made from a mixture of a variety of chili peppers with secondary ingredients such as shrimp paste, fish sauce, garlic, ginger, shallot, scallion, palm sugar, lime juice, and rice vinegar or other vinegar.”

Sambal in a glass jar. (Jaclynn Seah, 2015, The Occasional Traveler)

 

Because of its high fat content, in mild climate like Bend, Oregon it solidified. This solidified chunk looked like two large meat piece with sambal on top.

Well I was hungry so I thought “lets microwave these two tasty meat piece, and add rice on it.” Which I did and for some reason not realizing the sauce has melted completely (because there’s no meat, just pure sambal). I added rice, mixed it and in my blind hunger ate the rice with a chunk of the sauce.

The first bite made me realize that it was extremely spicy and my stomach immediately got hot. Then what had transpired just clicked in my mind and a tinge of regret for using the valuable sauce in one go. I finally finished it, the entire experience was a combination of pleasure and pain. Well mostly me burning my mouth and stomach. So I researched on journals on how to deal with the pain swiftly.

The spiciness of sambal or other chili sauce variants is capsaicin which is a vanilloid compound. This causes heat hyperalgesia a painful feeling akin to burning and heat sensitivity although mechanical hyperalgesia was not linked (Baad-Hansen et al, 2003). Because vanilloid compound which binds to transient receptor potential vanilloid type 1  (TRPV 1) that responds to heat and acidic pH stimuli (Pingle et al, 2007).  So I thought it was acidic and looked to see if antacid will help. Unfortunately since capsaicin does not easily dissociate in water, this is difficult to determine.  So I needed more research and looked up capsaicin interaction with antacids. I did not find much info on oral intake, so I will generalize findings from dermal exposure of capsaicin. Turns out there is some conflicting information regarding the effectiveness of antacids in dealing with pain, one study by Barry et al (2008) compares Mg-Al antacid, lidocaine gel, baby shampoo, milk, and water which results in no difference in pain relief for capsaicin exposure and noted that time is the best cure. While another study by Lee and Ryan (2003) noted that Mg-Al antacid helps a little bit. I was wondering then how fast does capsaicin go away if left alone? Well a pharmacokinetic analysis of capsaicin showed rapid plasma level reduction, with a mean population elimination half-life of 1.64 hours (Babbar et al, 2009). So pretty quick. According to Capsaicin MSDS Sheet it has oil/water partition coefficient of 3 which meant capsaicin is far more soluble in oil than water, which explains the popularity of milk (with casein explanation being suspect and without supporting evidence) as it contains ~3% fat. Sucrose water are also known to alleviate the pain (Nasrawi & Pangborn, 1990).

TLDR; swish your mouth with warm milk + butter mix  or oil + water mix and spit. Then swish with crushed antacid and spit again. Drink sugary stuff. After 1 hour eat antacid to alleviate some stomach pain if needed.

References

Pingle, S. C., Matta, J. A., & Ahern, G. P. (2007). Capsaicin receptor: TRPV1 a promiscuous TRP channel. In Transient Receptor Potential (TRP) Channels (pp. 155-171). Springer, Berlin, Heidelberg. https://link.springer.com/chapter/10.1007/978-3-540-34891-7_9

Barry, J. D., Hennessy, R., & McManus Jr, J. G. (2008). A randomized controlled trial comparing treatment regimens for acute pain for topical oleoresin capsaicin (pepper spray) exposure in adult volunteers. Prehospital Emergency Care, 12(4), 432-437. https://www.ncbi.nlm.nih.gov/m/pubmed/18924005/?i=5&from=/20579556/related

Lee, D. C., & Ryan, J. R. (2003). Magnesium–Aluminum Hydroxide suspension for the treatment of dermal Capsaicin exposures. Academic emergency medicine, 10(6), 688-690. https://www.ncbi.nlm.nih.gov/m/pubmed/12782534/?i=12&from=/20579556/related

Babbar, S., Marier, J. F., Mouksassi, M. S., Beliveau, M., Vanhove, G. F., Chanda, S., & Bley, K. (2009). Pharmacokinetic analysis of capsaicin after topical administration of a high-concentration capsaicin patch to patients with peripheral neuropathic pain. Therapeutic drug monitoring, 31(4), 502-510. https://www.ncbi.nlm.nih.gov/m/pubmed/19494795/?i=14&from=/20579556/related

Baad-Hansen, L., Jensen, T. S., & Svensson, P. (2003). A human model of intraoral pain and heat hyperalgesia. Journal of orofacial pain, 17(4). https://www.ncbi.nlm.nih.gov/m/pubmed/14737878/?i=28&from=/20579556/related

MSDS Sheet of Capsaicin http://www.sciencelab.com/msds.php?msdsId=9923296

Nasrawi, C. W., & Pangborn, R. M. (1990). Temporal effectiveness of mouth-rinsing on capsaicin mouth-burn. Physiology & behavior, 47(4), 617-623. https://www.sciencedirect.com/science/article/abs/pii/003193849090067E

Advertisements

Seneca’s Letter and Modern Ideas 1

I thought of this idea as I rewrote and summarized Letter 2,3,5,6,7,8,9,11. Some, if not all, of the ideas are very modern.

Richard Feynman did mention this in The Meaning of It All when he recognized that while the norms, rules, and knowledge of the time changes. Some philosophical ideas remain the same and relatively unchanged.

It was very interesting for me that I am reading letters older than dirt that has modern ideas about people I encountered in psychology, I will recollect what I thought of the Letters I rewrote. Some people will call these collective “common sense”. As an aside “common sense” is a bull-crap and useless idea, actually it is only useful if you want other people to agree with you on status quo without expediting any energy.

Seneca Letter II mentions of achieving constancy, sticking to a subject, book, or a person. Just like how it was in Dao De Jing’s philosophy ingrained in my Chinese culture that a master of one is superior to dabbler in many.

Seneca Letter V, the usual universal idea of trying to be better but not self-advertise, the quintessential collectivist community principle. Also to outwardly blend in with the people around you for practicality. Unless you’re Diogenes, then you can just never take a shower, say GTFO to Alexander The Great, and still be alive.

Seneca Letter VII, where he discourages us from crowds, as peer pressure will get to you. Which is true back then and true today. Especially in schools.

Seneca Letter IX, the Stoic philosopher Hecato mentioned “If you wish to be loved, love” the good ol’ Rule of Reciprocity mentioned by Cialdini and probably many others. Also highly effective.

Again On Seneca Letter IX. It is a mine of gold. If you remember painting an art and just being happily self-absorbed, he noted of it too. Which is nowadays called state of Flow by Csikszentmihalyi.

Another contemporary idea Seneca mentioned was of how what people say didn’t matter, but what they feel matters. A powerful variation of this was quoted to modern poet Maya Angelou “They won’t remember what you say, but they will remember how you make them feel”. For me with my terrible memory, this intuitively makes sense. I never forget how other people make me feel.

Letter XI, where he basically said to get a coach, a moral guide-person of sorts.
Which is really how world-class people gets to improve. In business they keep talking of the importance of mentors, in academics its a given that PhD candidate is mentored by the Professors. Athletes, top musician, and artists alike get coached continually. Like Itzhak Perlman coached by his wife, a quite demanding proposition she had to quit her main career path in the beginning. Even the well-known neurosurgeon Atul Gawande is trying to make this a thing, a continued coaching by requesting his retired Professor to coach him. Then proceeds to see an improvement in his surgical skills where he had stalled before.

No references for today, because this is me procrastinating from what I have to actually do.

Analytical Solution of Colebrook Equation

Colebrook Puzzle

In my introduction to fluid dynamics class (ME331) at OSU. I encountered a pesky implicit formula called the Colebrook equation [1] which is important to solve friction factor in turbulent flows, defined as

.

I put it in the back of my mind and forget about it. This is my first principle of problem-solving, if it’s too troublesome and not urgent, ignore it for a while. In other words, let the subconscious do it for you.

Then after a few weeks, I googled LambertW Wk(x) [2] function as I was thinking I have seen it somewhere. By this point, I have all the tools necessary to solve the Colebrook equation, but I still didn’t realize it. Until in my procrastination I meet a cousin of LambertW function, the Wright Omega ω(x) [3] function defined as the solution to

.

But this is very close to the form of which Colebrook equation is presented, and I was ready to solve the equation. The first thing that I think about was how to make the Colebrook equation in a similar form to the above. Getting Y is the hard part because I have to sync both the expression inside the log and an additive term. The X is easy, it’s just whatever multiplicative term is leftover after the process of “syncing” the Y has succeeded.

The first step is to rearrange the equation so that all f variable is on one side and ln(10) factor is pulled out so that

becomes

.

Now the solution looks closer, and I have to sync the variable inside the natural logarithm with the additive term. Here we exploit the fact that multiplying an expression in a logarithm yields an addition outside of it, and thus both expression is almost exactly of the same form minus an additive term that we can dump into X.

Thus the form that we wanted is finally in sight, the Wright Omega solution in the form of

.

Which by definition yields to

.

After further simplifications, the closed form of the friction factor in terms of Wright Omega is

.

As a bonus, using the definition of Wright Omega function in terms of LambertW function expressed as

.

We can write the solution using the principal branch of LambertW as

.

I never really looked if anyone had done this before, and did not want to find spoilers because it is such a perfect puzzle. Not too easy, not too hard.

Now I have finished the puzzle, I think it is decent to acknowledge prior works that do the similar/same approach to derivation. So today I delved into Google scholar to find who did it using a similar approach.

The earliest analytical solution of Colebrook equation that uses Wright Omega is a paper by Clamond in 2008 [4], this paper actually discussed the iterative solution in depth and compared different speed using different derivations. In order of rapidness, from the most to the least is a shifted Wright Omega, Wright Omega, and principal branch LambertW function.

Edit: Just to be sure, I did some test on the function using some random numbers on MATLAB iterative solver and the difference magnitude is much lower (~1e-11) than FSOLVE tolerance (1e-6). Which is good.

a = 3.7; b = 2.51;
iterative_colebrook = @(Re,rf) fsolve(@(f) (-2.*log10(rf./a+b./(Re.*sqrt(f)))-1./sqrt(f)),0.02,optimset('Display','off'));

wright_colebrook = @(Re,rf) (2/log(10)*wrightOmega(log(10)*Re*rf/(2*a*b)-log(2*b/(log(10)*Re)))-Re*rf/(a*b))^-2;
lambert_colebrook = @(Re,rf) (2/log(10)*lambertw(exp(log(10)*Re*rf/(2*a*b)-log(2*b/(log(10)*Re))))-Re*rf/(a*b))^-2;

wright_log_error = @(Re,rf) log10(abs(wright_colebrook(Re,rf)-iterative_colebrook(Re,rf)));
lambert_log_error = @(Re,rf) log10(abs(lambert_colebrook(Re,rf)-iterative_colebrook(Re,rf)));

I wrote this for a sort of fun writing and math practice. But if anyone happens to read this and have any corrections/questions/clarifications please email me at agustinuslaw@gmail.com.

References

[1] C. F. Colebrook, “Turbulent Flow in Pipes, with Particular Reference to the Transition between the Smooth and Rough Pipe Laws,” J. Inst. Civ. Eng. Lond., vol. 11, 1938–1939, pp. 133–156.

[2] Corless, R. M., Gonnet, G. H., Hare, D. E. G., Jeffrey, D. J., and Knuth, D. E.,1996, “On the LambertW Function,” Adv Comput Math, 5(1), pp. 329–359. https://www.researchgate.net/publication/2458364_On_the_Lambert_W_Function

[3] Corless, R., and Jeffrey, D., 2002, “The Wright Omega Function.,” pp. 76–89. https://www.researchgate.net/publication/221352071_The_Wright_omega_Function

[4] Clamond, D., 2008, “Efficient Resolution of the Colebrook Equation,” arXiv:0810.5564 [physics].

Eclipse from Mt Strawberry and Some Calculations

Mt. Strawberry, OR (alt. 2700 m | 9000 ft).

Just moments from the total solar eclipse I saw the shadow of the moon coming from the horizon. It crept slowly but then at one point it sprinted so fast I didn’t realize when the shadow hit me. It was a surreal moment. The darkness took hold and frigid wind just sends chills to my bones. Then I saw total solar eclipse for the first time in my life, a ring of light with an organic life-like corona dancing slowly about the dark sphere that used to be the sun. All life stood still at that moment. Time stops. And of course, everyone boomed in excitement. I was too mesmerized by the scenery to care. It was worth going up the mountains. I loved every moment of it.

I could swear there was a temperature drop of 10 C/20 F since I had to actually wear my double layer and everyone else started wearing jackets, or maybe it was the wind?

Back Of the Envelope Calculations

Since I am back home, tired, and did not wish to deal with real life yet. I made some rough calculations regarding the event I have seen.

At the top of the mountain (2700 m | 9000 ft) I was able to see about 195 km | 120 mi away (Bohren and Fraser, 1986). Since the speed of the shadow provided by Rao averaged around 3600 km/h | 2240 mi/h which is about 1 km/s | 0.62 mi/s this would mean it took roughly 3 minute for the shadow to travel from the horizon to my place (Rao, 2017). Of course the speed of the moon shadow itself actually varies because it was projected on an oblate spheroid surface of Earth. Becoming progressively slower with distance. But we can safely assume the speed to be constant because of the short distances involved relative to the total distance the moon shadow travels. Therefore I should see a constant speed moon shadow no?

Well no, since human perception, and in extension human vision does not operate linearly but is logarithmic, this is called Weber-Fechner’s Law which relates the relationship between actual change and perceived change (Weber, 1996). Similar to how a person looking at oncoming train as slow from afar, but it is actually progressing quicker than what it seems to be. This will only be apparent to the poor lad who tried to pass the train way too close and it’s already too late. So far it all makes sense, but I’m not sure it will still make sense when I wake up tomorrow.

I wrote this for a sort of journaling and fun writing practice. But if anyone happen to read this and have any corrections/questions/clarifications please email me at agustinuslaw@gmail.com.

References

Bohren, C. F., & Fraser, A. B. (1986). At what altitude does the horizon cease to be visible?. American Journal of Physics, 54(3), 222-227.
d ~ 3.57*sqrt(h). Where d is distance of sight in km, h is the altitude in m.
This approximates how far we can see from a given altitude disregarding the effect of refraction. Though that was also discussed in the paper. 

Rao, J. 2017. “How Long Will the 2017 Solar Eclipse Last? Depends Where You Are”. https://www.space.com/36388-total-solar-eclipse-2017-duration.html. This gave the average speed of moon’s shadow. In Oregon that will be about 3600 km/h ~1 km/s and 2240 mi/h ~ 0.62 mi/s. 

Weber, E. H. (1996). EH Weber on the tactile senses. Psychology Press.
This elaborates on the Weber-Fechner’s Law on logarithmic human perception which I used qualitatively here.

Research: Practice on Matlab Day 5 (MIT 6.094)

Motivation is a hard thing to keep, I did finish the homework by day 5 but then I put off publishing this because my mind was preoccupied by something else. I bought a car, and now I have to fix it and practice on it so I can pass driving exam.

For now I’ll just upload without much description. I also skipped the optional homework.. Someday I’ll do it..

Resources:

Lecture 4 and Homework 4

LECTURE 4                                                                                  

Statistics

Simple histogram

close all;
%populate the scores
scores = 100*rand(1,100);

%histogram with bin centered at 5,15,...,95
hist(scores, 5:10:95);

%returns the number of occurences between specified bin edges
   %[0,10), [10,20), [30,40), ..., [90,100],
N=histc(scores, 0:10:100);

%plot barchart
bar(0:10:100, N, 'r');

Brownian motion

Simulation of Brownian random walk, then visualize it as histogram.

a=zeros(0,10000);
x=0;
for n=1:10000
    if rand()<0.5
        %move left
        x=x-1;
    else
        %move right
        x=x+1;
    end

    a(n)=x;
end
figure(1);
%plot a 50 bin histogram
hist(a,50);

Sparse matrices

Making matrices sparse so they don’t take up a lot of memory, and more efficient to process. whos details the matrices size, type and attributes.

a=zeros(100);
a(1,3)=10;
a(21,5)=pi;
whos a;
b=sparse(a);
whos b;

Data structures

Practice using cell and struct as a way to organize data

clc;
clear all;

%initializing empty cell
a=cell(1,2);
%intializing a cell
c={'hello',[1 2 3 4], rand(3,2)};
%accessing a cell element
c{1,2};

%initialize empty struct
s=struct([]);
%add fields to the struct
m(1).name = 'hell joe';
m(1).scores = [100 98 50];
m(1).year = 'G3';

people=struct('name',{'John','Mary','Leo'},...
    'age',{32,27,18},'childAge',{[2,4],1,[]});
people(3).name;

Sentence generation 

Illustrates the convenience of cell for dealing with string.

clc;
%create cell
data=cell(3,2);
data(:,1)={'andy';'warhol';'zumba'};
data(:,2)={'stinky';'smelly';'shiny'};

%display sentence
disp([data{ceil(rand*3),1} ' is ' data{ceil(rand*3),2}]);

Figure handles

Changing the figure setting using handles. Anything in the figure and plot can be customized via handles.

clc; close all;
%handle for the plotted line
L=plot(1:10,rand(1,10));
%handle for the current axis
A=gca;
%handle for current figure
F=gcf;
%current property values
get(L);
%change property values
set(A, 'FontName', 'Cambria', 'XScale', 'log');
set(L,'LineWidth',1.5,'Marker','*');

Read write images

Matlab is quite good for processing images.

%write RGB image
imwrite(rand(500,500,3),'test1.jpg');
%write indices and colormap
imwrite(ceil(rand(200)*256),jet(256),'test2.jpg');
%import image to matlab
im=imread('test1.jpg');

Making animations

The script uses movie2avi which has been deprecated, use videowriter instead.

%making animation by plotting frame by frame
    ...with pause in between
close all;
for t=1:100
    imagesc(rand(500));
    colormap(gray);
    pause(.01);
end

close all;
for n=1:100
    imagesc(rand(500));
    colormap(jet);
    M(n)=getframe;
end

%play the movie 2 times at 100fps
movie(M,2,100);

%save as avi files
movie2avi(M,'testMovie.avi','FrameRate',100,'VideoCompresssionMethod','Motion JPEG AVI');

Animated gifs

%using imwrite to save animated gif
temp=ceil(rand(300,300,1,10)*256);

imwrite(temp,jet(256),'l4-4-testGif.gif',...
    'delaytime',1/60,'loopcount',100);

Revisiting random surface as animation

This isn’t in the lecture, but I was thinking of using what I’ve learned so far by making a smooth surface using random numbers and  interp2, then have a smooth time-transition between surfaces using interp, and then customize the plot using handles such as gca, gcf, and set.

clc;close all;

%PARAMETER
%Peaks and Valleys
N=5;
%Length
T=3;
%Graph interpolation resolution
rg=0.1;
%Time interpolation resolution
rt=0.01;
%GIF settings
gifDelay=0;
gifLoop=Inf;

%INITIALIZATION
%Time dimension for Z and sZ
t0=1:T;
ts=1:rt:T;
%meshgrid for Z0
[X0,Y0]=meshgrid(1:N,1:N);
%meshgrid for interpolated Zi
[X,Y]=meshgrid(1:rg:N,1:rg:N);
%initialize for speed
Z=zeros(length(X),length(Y),T);

%CREATE GRAPH INTERPOLATION
for i=1:T
    %Peaks and valleys
    Z0=rand(N);
    %interpolate values
    Z(:,:,i)=interp2(X0,Y0,Z0,X,Y,'cubic');
end

%SMOOTHEN TRANSITION BETWEEN GRAPH
%Permute to get interpolated dimension first:
permutedZ = permute(Z,[3 1 2]);
%create sZ for smooth transition
iPermutedZ = interp1(t0,permutedZ,ts);
%permute to return to original dimension arrangement
sZ = permute(iPermutedZ,[2 3 1]);

%ANIMATE
h=figure;
for i=1:length(ts)
    %create surface
    surf(X,Y,sZ(:,:,i))
    hold on
    contour(X,Y,sZ(:,:,i))
    hold off

    %figure settings
    colormap(jet);
    set(gca,'color','k');
    set(gcf,'color','k');
    shading interp;
    zlim([0,1]);
    %colorbar;

%             %PLOT->GIF IMPLEMENTATION 1
%             % Capture the plot as an image
%             frame = getframe(h);
%             im = frame2im(frame);
%             [imind,cm] = rgb2ind(im,256);
%             % Write to the GIF File
%             if i==1
%                 imwrite(imind,cm,'revisitSurfaceImwrite.gif','gif','DelayTime',gifDelay,'Loopcount',gifLoop);
%             else
%                 imwrite(imind,cm,'revisitSurfaceImwrite.gif','gif','DelayTime',gifDelay,'WriteMode','append');
%             end

    %PLOT->GIF IMPLEMENTATION 2, (better but still choppy)
    %add to gif
    if i==1
        gif('l4-revisitSurfaceGif.gif','DelayTime',gifDelay, ...
            'LoopCount',gifLoop,'frame',gcf);
    else
        gif
    end

    %progress
    disp([num2str(round(i/length(ts)*100),4) ' %']);
end

%view the animation
web('revisitSurfaceGif.gif');
web('revisitSurfaceImwrite.gif');

Plot ten lines

the legend function didn’t work for me for the longest time until I realized I made a variable with the same name, Matlab prioritizes variables so functions don’t work anymore.

function l4plotTenlines
    %Exercise in debugging
    clc;
    % Click on – next to line numbers in MATLAB files
    % ¾ Each red dot that appears is a breakpoint
    % ¾ Run the program
    % ¾ The program pauses when it reaches a breakpoint
    % ¾ Use the command window to probe variables
    % ¾ Use the debugging buttons to control debugger
    clear all;
    close all;

    %make x vector
    x=0:.01:1;
    %legendNames=char(ones(10,8))

    %make a figure and plot 10 random lines, labeling each one
    figure;
    %legendNames = {};   %initialize
    for n=1:10
        plot(x,polyval(rand(3,1),x),'color',rand(1,3));
        hold on;
        legendN=['Line ' num2str(n)];
        legendNames(n,1:length(legendN))=legendN;
    end
    %label the graph
    xlabel('X');
    ylabel('Y');
    title('Ten Lines Plot');

    legendNames
    legend(legendNames);

end

HOMEWORK 4

Problem 1 – random numbers

% Homework 4 - 1
% Random variables. Make a vector of 500 random numbers from a
    % Normal distribution with mean 2 and standard deviation 5 (randn).

% After you generate the vector, verify that the sample mean and standard
    % deviation of the vector are close to 2
    % and 5 respectively (mean, std).

%I played to just see the most extreme values after 100000 runs
clc;
sigma=5;
mu=2;
v=randn(500,1).*sigma+mu;
for i=1:100000
    avg=mean(v);
    if abs(avg-mu)>abs(avgD-mu)
        avgD=avg
    end
    sd=std(v);
    if abs(sd-sigma)>abs(sdD-sigma)
        sdD=sd
    end
end

Problem 2 – flip coins

Although simple, I thought this was a good illustration to show how probability works considering the course was for lower division students.

% Homework 4 - problem 2
% Flipping a coin.
% Write a script called coinTest.m to simulate sequentially flipping a coin 5000 times.
    % I will define head as 1 and tail as 0
close all;
rep=10000;
flips=round(rand(1,rep));
% Keep track of every time you get ‘heads’ and plot the running estimate of the probability of getting ‘heads’ with this coin.
headEvent=cumsum(flips);
totalEvent=1:rep;
headProbability=headEvent./totalEvent;
% Plot this running estimate along with a horizontal line at the expected value of 0.5, as below.
plot(totalEvent,headProbability,totalEvent,0.5*ones(1,rep));
% This is most easily done without a loop (useful functions: rand, round, cumsum).

Problem 3 – Poisson distribution

Plotting Poisson random numbers (poissrnd) vs the probability density function (poisspdf). Also hist is deprecated in favor of histogram.

% Homework 4 Problem 3
% Histogram. Generate 1000 Poisson distributed random numbers with parameter ? = 5(poissrnd).
close all;
clc;
lambda=5;
count=1000;
x=poissrnd(lambda,1,count);
xMin=min(x);
xMax=max(x);
% Get the histogram of the data and normalize the counts so that the histogram sums to 1
    ...(hist – the version that returns 2 outputs N and X, sum).
H=histogram(x,'BinEdges',[xMin:xMax],'Normalization','probability');
% Plot the normalized histogram (which is now a probability mass function) as a bar graph (bar).


% Hold on and also plot the actualPoisson probability mass function with ? = 5 as a line (poisspdf).
hold on;
index=linspace(xMin,xMax,count);
plot(index,poisspdf(lambda,index),'r');
% You can try doing this with more than 1000 samples from the Poisson distribution to get better agreement between the two.

Problem 4 – cells

Cells are interesting that they can fit anything inside one of their element. Also accessing cell element as a cell requires the use of round brackets (). While accessing cell element as their respective data types requires the use of curly brackets {}.

% Practice with cells.
% Usually, cells are most useful for storing strings, because the length of each string can be unique.
% a. Make a 3x3 cell where the first column contains the names: ‘Joe’, ’Sarah’, and ’Pat’, the
% second column contains their last names: ‘Smith’, ‘Brown’, ‘Jackson’, and the third
% column contains their salaries: $30,000, $150,000, and $120,000. Display the cell using
% disp.
employee={'Joe','Smith',30000;'Sarah','Brown',150000;'Pat','Jackson',120000};
disp(employee)

clear employee2;
employee2=cell(3,3);
employee2(:,1)={'Joe';'Sarah';'Pat'};
employee2(:,2)={'Smith';'Brown';'Jackson'};
employee2(:,3)={30000,150000,120000};
disp(employee2)

% b. Sarah gets married and changes her last name to ‘Meyers’. Make this change in the cell
% you made in a. Display the cell using disp.
employee(2,2)={'Meyers'}
% c. Pat gets promoted and gets a raise of $50,000. Change his salary by adding this amount
% to his current salary. Display the cell using disp.
employee{3,3}=employee{3,3}+50000
% The output to parts a-c should look like this:

Problem 5 – structs

Illustrates the utility of structs by using dir to list the file properties within the current work directory.

function hw4p5structs
% Using Structs. Structs are useful in many situations when dealing with diverse data. For
% example, get the contents of your current directory by typing a=dir;
A=dir;
% a. a is a struct array. What is its size? What are the names of the fields in a?
disp(['size of A: ' num2str(size(A))]);

% b. write a loop to go through all the elements of a, and if the element is not a directory,
% display the following sentence ‘File filename contains X bytes’, where filename is the
% name of the file and X is the number of bytes.

for i=1:length(A)
   disp(['File ' A(i).name ' contains ' num2str(A(i).bytes) ' bytes']);
end

% c. Write a function called displayDir.m, which will display the sizes of the files in the
% current directory when run, as below:
end

Problem 6 – handles 

Customizing plots using handles ~ gca, gcf, set, and grid

% Handles. We’ll use handles to set various properties of a figure in order to make it look like this:
% One sine wave from 0 to 2? X values (in terms of ?) Sin(X)
% a. Do all the following in a script named handlesPractice.m
% no
% b. First, make a variable x that goes from 0 to 2? , and then make y=sin(x).
res=200;
x=linspace(0,2*pi,res);
y=sin(x);
% c. Make a new figure and do plot(x,y,’r’)
figure;
plot(x,y,'r');
% d. Set the x limit to go from 0 to 2? (xlim)
xlim([0,2*pi]);
% e. Set the xtick property of the axis to be just the values [0 pi 2*pi], and set
% xticklabel to be {‘0’,’1’,’2’}. Use set and gca
set(gca,'xtick',[0,pi,2*pi],'xticklabel',{'0','1','2'});
% f. Set the ytick property of the axis to be just the values -1:.5:1. Use set and gca
set(gca,'ytick',-1:0.5:1);
% g. Turn on the grid by doing grid on.
grid on;
% h. Set the ycolor property of the axis to green, the xcolor property to cyan, and the
% color property to black (use set and gca)
set(gca,'xcolor','cyan','ycolor','black','color','black');
% i. Set the color property of the figure to a dark gray (I used [.3 .3 .3]). Use set and gcf
set(gcf,'color',[.3 .3 .3]);
% j. Add a title that says ‘One sine wave from 0 to 2?’ with fontsize 14, fontweight
% bold, and color white. Hint: to get the ? to display properly, use \pi in your string.
% Matlab uses a Tex or Latex interpreter in xlabel, ylabel, and title. You can do all this just
% by using title, no need for handles.
title('One sine wave from 0 to 2\pi', 'fontsize', 14,'fontweight', 'bold','color','white');

% k. Add the appropriate x and y labels (make sure the ? shows up that way in the x label)
% using a fontsize of 12 and color cyan for x and green for y. Use xlabel and ylabel
xlabel('multiples of \pi','fontsize',12,'color','cyan');
ylabel('sin(x\pi)','fontsize',12,'color','green');

% l. Before you copy the figure to paste it into word, look at copy options (in the figure’s Edit
% menu) and under ‘figure background color’ select ‘use figure color’.

Problem 7 – image

This question was problematic for me because at first I didn’t really understand how meshgrid really works and I flipped rows and columns (as x and y variables). and had it fine and dandy the whole time because I was dealing with square matrices/images. Until I had to deal with non-square images in this exercise. Also this exercise taught me how some image processors resizes images and results in smoothed but blurry picture. It was 2D interpolation all along.

% Image processing
% Write a function to display a color image, as well as its red, green, and blue layers separately
function im=hw4p7image(filename)
    % The function declaration should be im=displayRGB(filename)
    % filename should be the name of the image (make the function work for jpg images only)
    clc;

    %if function is called with no argument
    if nargin == 0 %
        %filename = 'test1.jpg'      %SUCCESS: This is a square image nxn
        filename = 'https://upload.wikimedia.org/wikipedia/commons/d/da/The_Parthenon_in_Athens.jpg'  %FAIL: This is a rectangular image mxn
    end

    %this function can only accept JPG
    [pathstr,name,ext] = fileparts(filename);
    %somehow using short-circuit
    if length(ext) ~= 4 | ext ~= '.jpg'
        error('The file received is not jpg');
    end

    % im should be the final image returned as a matrix
    % To test the function, you should put a jpg file into the same directory as the function and run it with the filename
    % (include the extension, for example im=displayRGB(‘testImage jpg’))
    % You can use any picture you like, from your files or off the internet
    % Useful functions: imread, meshgrid, interp2, uint8, image, axis equal, axis tight
    imageOriginal=imread(filename);
    % a To make the program work efficiently with all image sizes, first interpolate each color layer of the original image
    % so that the larger dimension ends up with 800 pixels
    % The smaller dimension should be appropriately scaled so that the length:width ratio stays the same
    % Use interp2 with cubic interpolation to resample the image
    imageSize=size(imageOriginal);

    %convert to double
    imageOriginal=double(imageOriginal);
    %the parameter for new image
    newLength=800;

    if max(imageSize) ~= newLength
        %get all information from original image
        imageSize=size(imageOriginal);
        x=1:imageSize(2);
        y=1:imageSize(1);
        [X,Y]=meshgrid(x,y);

        %new dimension
        newSize=floor(imageSize/max(imageSize)*newLength);
        xn=linspace(1,imageSize(2),newSize(2));
        yn=linspace(1,imageSize(1),newSize(1));
        [Xn,Yn]=meshgrid(xn,yn);
        %new image
        for i=1:3
            imageResampled(:,:,i)=interp2(X,Y,imageOriginal(:,:,i),Xn,Yn,'cubic');
        end

    else
        imageResampled=imageOriginal;
    end

    %reconvert both images as uint8
    imageOriginal=uint8(imageOriginal);
    imageResampled=uint8(imageResampled);

    close all;
    %for comparison
    figure(1);image(imageOriginal);title('before resampling');
%         figure(2);image(imresize(imageOriginal,newSize));title('built in resize function');
    % b Create a composite image that is 2 times as tall as the original, and 2 times as wide
    % Place the original image in the top left, the red layer in the top right, the green layer in the bottom left, and the blue layer in the bottom right parts of this composite image
    % The function should return the composite image matrix in case you want to save it as a jpg again (before displaying or returning, convert the values to unsigned 8-bit integers using uint8)
    % Hint: To get just a single color layer, all you have to do is set the other two layers to zero
    % For example if X is an MxNx3 image, then X(:,:,2)=0; X(:,:,3)=0; will retain just the red layer
    % Include your code and the final image in your homework writeup
    imageRed=imageResampled;
    imageRed(:,:,2)=0;imageRed(:,:,3)=0;
    imageGreen=imageResampled;
    imageGreen(:,:,1)=0;imageGreen(:,:,3)=0;
    imageBlue=imageResampled;
    imageBlue(:,:,1)=0;imageBlue(:,:,2)=0;
    imageDecomposition=[imageResampled,imageRed;imageGreen,imageBlue];

    figure(3);
    image(imageDecomposition);title('image decomposition');
end

 

Problem 7 troubleshooting – interpExperiment

To troubleshoot 7. I had to realize first I might be doing something sketchy with meshgrid or interp2 so I made a smaller function to deal with those two. Voila, found my problem.

clc; clear all; close all;
%original dimensions
setOriginalSize=[3,5,3];       %param1

%random numbers
imageOriginal=rand(setOriginalSize);

%get all information from originalMat
imageSize=size(imageOriginal);
x=1:imageSize(2);
y=1:imageSize(1);
[X,Y]=meshgrid(x,y);

%new dimension
newLength=100;
newSize=floor(imageSize/max(imageSize)*newLength);
xn=linspace(1,imageSize(2),newSize(2));
yn=linspace(1,imageSize(1),newSize(1));
[Xn,Yn]=meshgrid(xn,yn);
%new image
for i=1:3
    imageResampled(:,:,i)=interp2(X,Y,imageOriginal(:,:,i),Xn,Yn,'cubic');
end
%--------------------------------------------------------------------------
figure(1);
image(imageOriginal);
figure(2);
image(imageResampled);

Problem 8 – Brownian diffusion

I thought this was pretty neat. Illustrating Brownian diffusion only using normal distribution (randn).

% Animation: Brownian motion
% Write a function with the following declaration:
function hw4p8browniandiffusion(N,iter,delay,sigma)
    clc;
% The function takes in a single input N, which is an integer specifying the number of points in the simulation
    if nargin == 0
        N=100;
        iter=1000;
        delay=0.01;
        sigma=0.005;
    elseif nargin == 1
        iter=1000;
        delay=0.01;
        sigma=0.005;
    elseif nargin == 2
        delay=0.01;
        sigma=0.005;
    elseif nargin == 3
        sigma=0.005;
    end
% All the points should initially start at the origin (0,0)
    x=zeros(1,N);
    y=zeros(1,N);
% Plot all the points on a figure using ‘.’ markers, and set the axis to be
    ...square and have limits from -1 to 1 in both the x and y direction
    close all;
    P=plot(x,y,'.');
    xlim([-1,1]);ylim([-1,1]);
% To simulate Brownian motion of the points, write a 1000-iteration loop which will calculate a new x and y
    ...position for each point and will display these new positions as an animation

    for i=1:iter
% The new position of each point is obtained by adding a normally distributed random variable with a standard deviation of 0
    ... 005 to each x and y value (use randn; if you have 100 points, you need to add 100 distinct random values to
    ...the x values and 100 distinct random values to the y values)
        x=x+randn(1,N)*sigma;
        y=y+randn(1,N)*sigma;
% Each time that the new positions of all the points are calculated, plot them on the figure and pause for 0
    ... 01 seconds (it’s best to use set and the line object handle in order to
    ...update the xdata and ydata properties of the points, as mentioned in lecture)
        set(P,'XData',x,'YData',y,'Marker','.');
        pause(delay);
    end
% What you will see is a simulation of diffusion, wherein the particles randomly move away from the center of the figure
% For example, 100 points look like the figures below at the start, middle, and end of the simulation:
end

Optional Problem – Julia Set animation

 

Research: Practice on Matlab Day 3 (MIT 6.094)

Wednesday was my third day of practice, I finished all the code but couldn’t get a solution on the problem regarding Hodgkin-Huxley differential equations using Matlab ode45 solver. It was so slow. Well I did get a solution for timescale 0 to 0.5 millisecond. But as I increased it to 20 ms, I left it overnight and when I woke up it still ran!  In fact it went up to about 70.5 million matrix element, that couldn’t be right. So I researched about Hodgkin-Huxley and found that a paper (Chance & Clore, 1981)  on it describes the HH equation partial differential equation as a stiff differential equations. Which means there are terms that can lead to rapid variation and operates on large and small timescales. With the usual numerical methods the solution set will become unstable. Or solved maybe very slowly.

I thought this was the problem, so maybe I should use ode23s or  ode15s solver, but because I don’t want to just brute force it and leave it overnight again I went to campus to ask Dr. Kyle Webb on my problem. So he did agree that if it was stiff I should try and use the ode15s stiff solver but since the question explicitly says I am recommended to use ode45 it is far more likely that I made a mistake in my implementation.

So he did what I should’ve done in the beginning. He checked all my equations and compared with the problem. Lo and behold, I forgot a negative sign on one of the equations.

The lecture (link) exercises and homework (link) associated with Wednesday’s practice is posted below

Lectures

polyFittingLecture.m – this illustrates polynomial fitting, representation, and evaluation

%Matlab Lecture 3
clc;
close all;

%y=x^2
p=[1,0,0];
x=-4:0.1:4;
%evaluate polynomial
y=polyval(p,x);

%add noise
yNoisy=y+rand(size(x));
%plot noisy data
plot(x,yNoisy,'.');
%curve fit the noisy data
pFit=polyfit(x,yNoisy,2);
%evaluate fitted polynomial
yFit=polyval(pFit,x);
hold on;
plot(x,yFit,'r-');

chemode.m – on setting up ordinary differential equations (ODE)

%Matlab Lecture 3 - ODE
function dydt=chemode(t,y)
    dydt=zeros(2,1);
    dydt(1)=-10*y(1)+50*y(2);
    dydt(2)=10*y(1)-50*y(2);
end

solveChemode.m – on solving ODEs

%Matlab Lecture 3 - ODE Solver
clc;
close all;
%solve the chemode equations from t=0 to 0.5 and
%initial value y(1)=0 and y(2)=1
[t,y]=ode45('chemode',[0,0.5],[0,1]);
%plotting the amount of y(1) or substance A over time
plot (t,y(:,1),'k','LineWidth',1.5); hold on;
%plotting the amount of y(2) or substance B over time
plot (t,y(:,2),'r','LineWidth',1.5);
%formatting
legend('A','B');
xlabel('Time (s)');
ylabel('Amount of chemical (g)');
title('Chem reaction');

pendulum.m – tricks to deal with higher order ODEs

basically rephrase the higher order term as multiple variables and equations with only single order terms

%Matlab Lecture 3 - higher order ODE

%d(dtheta) + g/L sin(theta) = 0
%has to be broken down to a system of first order equations

function dxdt = pendulum(t,x)
    %0. Constant
    L=1;        %length of pendulum
    g=9.81;     %gravity

    %1. Input setup
    theta=x(1); %the angle
    gamma=x(2); %derivative of the angle, defined below

    %2. ODE
    dtheta=gamma;
    dgamma=-(g/L)*sin(theta);

    %3. The derivatives setup
    dxdt=zeros(2,1);    %derivatives initialization
    dxdt(1)=dtheta;
    dxdt(2)=dgamma;

end

Homework

Problem 1 – Linear System of Equations

% Homework 3 Problem 1
% Linear system of equations. Solve the following system of equations using \. Compute and
% display the error vector
% 3a + 6b + 4c = 1
% a + 5b = 2
% 7b + 7c = 3

%rephrase knowns
A=[3,6,4;1,5,0;0,7,7];
b=[1;2;3];

%check the rank
rankProblem=rank(A);
disp('rankProblem');
disp(rankProblem);

%Ax=b -> x=A\b
x=A\b;
disp('x');disp(x);

%Error vector abs(b-Ax)
Error = abs(b-A*x);
disp('Error');disp(Error);

Problem 2 – Numerical Integration Using Trapezoids and Adaptive Quadrature

% Homework 3 Problem 2
clc;close all;
% what is the value of int_0^5 xe^(-x/3)dx
% using adaptive quadrature
quadInt=integral(@(x) x.*exp(-x./3),0,5);
% using trapezoidal rule
x=linspace(0,5,100);
y=x.*exp(-x./3);
trapzInt=trapz(x,y);
figure;plot(x,y,'.-');
% compute difference with -24e^(-5/3)+9
quadError=quadInt-(-24*exp(-5/3)+9);
trapzError=trapzInt-(-24*exp(-5/3)+9);
%display results
disp(['quadInt=' num2str(quadInt) ', trapzInt=' num2str(trapzInt)]);
disp(['quadError' num2str(quadError) ', trapzError=' num2str(trapzError)]);

Problem 3 – Inverse

% Homework 3 Problem 3
% Computing the inverse. Calculate the inverse of A and verify that when you multiply the
% original matrix by the inverse, you get the identity matrix (inv). Display the inverse matrix as well
% as the result of the multiplication of the original matrix by its inverse

A=[1,2;3,4];
invA=inv(A)
I=A*invA

Problem 4 – Polynomial Fit

% Homework 3 Problem 4
% Fitting polynomials. Write a script to load the data file randomData.mat (which contains
    % variables x and y) and
clc; clear all;
load randomData.mat;
% fit first, second, third, fourth, and fifth degree polynomials to it.
[p1,s1,mu1]=polyfit(x,y,1);
[p2,s2,mu2]=polyfit(x,y,2);
[p3,s3,mu3]=polyfit(x,y,3);
[p4,s4,mu4]=polyfit(x,y,4);
[p5,s5,mu5]=polyfit(x,y,5);
%evaluate fitted polynomial
[y1,delta1] = polyval(p1,x,s1,mu1);
[y2,delta2] = polyval(p2,x,s2,mu2);
[y3,delta3] = polyval(p3,x,s3,mu3);
[y4,delta4] = polyval(p4,x,s4,mu4);
[y5,delta5] = polyval(p5,x,s5,mu5);

% Plot the data as blue dots on a figure, and plot all five polynomial fits using lines of different colors on
    % the same axes.
close all;
figure;
plot(x,y,'k.',x,y1,'c',x,y2,'b',x,y3,'m',x,y4,'r',x,y5,'y');
% Label the figure appropriately.
xlabel('x');
ylabel('y');
title('polynomial fit');
legend('data','p1','p2','p3','p4','p5');
% To get good fits, you’ll have to use the centering and scaling version of polyfit
    % (the one that returns three arguments, see help) and its counterpart in polyval (the one that accepts
    % the centering and scaling parameters).

Problem 5 – ODEs – Hodgkin Huxley Model of Neurons*

This is is my favorite problem, because I got to learn some neat stuff about how neurons behaves, about differential equations, numerical methods associated with them, and numerical problems.

This is the ODE setup

% Homework 3 Problem 5 - Hodgkin Huxley ODE
function dydt=HHode(t,x)
    % Constants
    C = 1 ;
    GK = 36 ;
    GNa = 120 ;
    GL = 0.3 ;
    EK = -72 ;
    ENa = 55 ;
    EL = -49.4 ;

    %Input Setup
    n=x(1);
    m=x(2);
    h=x(3);
    V=x(4);

    % ODE
    dndt=(1-n)*alphan(V)-n*betan(V);
    dmdt=(1-m)*alpham(V)-m*betam(V);
    dhdt=(1-h)*alphah(V)-h*betah(V);
    dVdt=-1/C*(GK*n^4*(V-EK)+GNa*m^3*h*(V-ENa)+GL*(V-EL)); %K.Webb: fix negative sign mistake

    % Derivative Setup
    dydt=[dndt;dmdt;dhdt;dVdt];
end

This is the solver

% Homework 3 Problem 5 - Main
% Hodgkin-Huxley model of the neuron
clc; close all;
% You will write an ODE file to describe the spiking of a neuron, based on the equations developed by
    %Hodgkin and Huxley in 1952 (they received a Nobel Prize for this work)
% The main idea behind the model is that ion channels in the neuron’s membrane have voltage-sensitive gates
    %that open or close as the transmembrane voltage changes
% Once the gates are open, charged ions can flow through them, affecting the transmembrane voltage
% The equations are nonlinear and coupled, so they must be solved numerically

% a. Download the HH.zip file from the class website and unzip its contents into your homework folder
% This zip folder contains 6 m-files: alphah.m, alpham.m, alphan.m, betah.m, betam.m, betan.m
% These functions return the voltagedependent opening ? (V ) ) and closing ( ? (V ) ) rate constants
    % for the h, m, and n ?( V ) gates: ? ( V )

% b. Write an ODE file to return the following derivatives
% 	and the following constants ( C is membrane capacitance, G are the conductances and E are the reversal potentials of the
% 	potassium ( K ), sodium ( Na ), and leak ( L ) channels): C = 1 GK = 36 GNa = 120 GL = 0.3 EK = ?72 ENa = 55 EL = ?49.4

% c. Write a script called HH.m which will solve this system of ODEs and plot the transmembrane voltage
% First, we’ll run the system to steady state
% Run the simulation for 20ms (the timescale of the equations is ms) with initial values:
    %n = 0.5; m = 0.5;h = 0.5;V = ?60 (ode45)
[t,y]=ode45('HHode',[0,20],[0.5,0.5,0.5,-60]);

% Store the steady-state value of all 4 parameters in a vector called ySS
n=y(:,1);
m=y(:,2);
h=y(:,3);
V=y(:,4);
ySS=y(end,:);

% Make a new figure and plot the timecourse of V (t) to verify that it reaches steady state by the end of the simulation
figure(1);
plot(t,V);
xlabel('t [ms]');ylabel('Transmembrane Voltage [mV]');title('Approaching Steady State');
% d. Next, we’ll explore the trademark feature of the system: the all-or-none action potential
% Neurons are known to ‘fire’ only when their membrane surpasses a certain voltage threshold
% To find the threshold of the system, solve the system 10 times, each time using ySS as ...
    %the initial condition but increasing the initial value of V by 1, 2, … 10 mV from its steady state value
figure(2);
clear V;
Vthres=10;
for i=1:10
   y0=ySS+[0,0,0,i];
   [t,y]=ode45('HHode',[0,20],y0);

   % After each simulation, check whether the peak voltage surpassed 0mV, and if it did, plot the voltage with a red line,
    %but if it didn’t, plot it with a black line
   if max(y(:,4))>0
       plot(t,y(:,4),'r');
       hold on;
       Vthres=min(Vthres,i);    %record the voltage threshold
   else
       plot(t,y(:,4),'k');
       hold on;
   end
end

%Figure formatting
xlabel('t [ms]');ylabel('Transmembrane Voltage [mV]');title('Threshold Behavior');

% Plot all the membrane voltage trajectories on the same figure, like below
% We see that if the voltage threshold is surpassed, then the neuron ‘fires’ an action potential;
    % otherwise it just returns to the steady state value
% You can zoom in on your figure to see the threshold value (the voltage that separates the red lines from the black lines).

Problem 6 Optional – Recursion – Julia Set*

My second favorite problem. Just because it was made using a very simple process over and over again to produce intricate patterns. This particular problem uses a complex quadratic polynomial with complex parameter c, the simplest nontrivial case is:

The Julia set of f(z) for c =-0.297491 + i0.641051 is shown below

The escape velocity function

% a. It has been shown that if the modulus of zn becomes larger than 2 for some n then it is guaranteed
    % that the orbit will tend to infinity.
% The value of n for which this becomes true is called the ‘escape velocity’ of a particular z0 .
% Write a function that returns the escape velocity of a given z0 and c .
% The function declaration should be: n=escapeVelocity(z0,c,N) where N is the maximum allowed
% escape velocity (basically, if the modulus of zn does not exceed 2 for n<N return N as the escape velocity.
% This will prevent infinite loops). Use abs to calculate the modulus of a complex number
function n=escapeVelocity(z0,c,N)
    clear i;        %remove potential problems with variables vs imaginary number
    z=z0;
    n=0;            %invariant n : subcript of z
    while (abs(z)<2) && (n ~= N)
        zLast = z;      %store previous number, for speed
        z = z^2+c;      %the dynamic recursion
        n = n + 1;
        if z==zLast     %premature end if it becomes constant, for speed
            n = N;
        end
    end
end

The Julia Set

% Homework 3 Optional Problem 6 - Julia Sets
% Optional, but highly recommended: Julia Sets
% In this problem you will generate quadratic Julia Sets. The following description is adapted from Wikipedia at http://en.wikipedia.org/wiki/Julia
% For more information about Julia Sets please read the entire article there.
% Given two complex numbers, c and z0 , we define the following recursion:
% z_n = (z_(n?1))^2 + c
% This is a dynamical system known as a quadratic map
% Given a specific choice of c and z0 , the above recursion leads to a sequence of complex numbers z1, z2 , z3 … called the orbit of z0
% Depending on the exact choice of c and z0 , a large range of orbit patterns are possible
% For a given fixed c , most choices of z0 yield orbits that tend towards infinity
% (That is, the modulus z grows without limit as n increases.)
% For some values of c certain choices of z0 yield orbits n that eventually go into a periodic loop
% Finally, some starting values yield orbits that appear to dance around the complex plane, apparently at random
% (This is an example of chaos.) These starting values, z0 , make up the Julia set of the map, denoted Jc
% In this problem, you will write a MATLAB script that visualizes a slightly different set,
% called the filled-in Julia set (or Prisoner Set), denoted Kc , which is the set of all z0 with orbits which do not tend towards infinity
% The "normal" Julia set Jc is the edge of the filled-in Julia set
% The figure below illustrates a Julia Set for one particular value of c
% You will write MATLAB code that can generate such fractals in this problem.

% a. It has been shown that if the modulus of zn becomes larger than 2 for some n then it is guaranteed that the orbit will tend to infinity.
% The value of n for which this becomes true is called the ‘escape velocity’ of a particular z0 . Write a function that returns the escape velocity of a given z0 and c .
% The function declaration should be: n=escapeVelocity(z0,c,N) where N is the maximum allowed escape velocity (basically, if the modulus of zn does not exceed 2 for n<N return N as the escape velocity.
% This will prevent infinite loops). Use abs to calculate the modulus of a complex number

% b. To generate the filled in Julia Set, write the following function M=julia(zMax,c,N).
% zMax will be the maximum of the imaginary and complex parts of the various z0 ’s for which we will compute escape velocities.
% c and N are the same as defined above, and M is the matrix that contains the escape velocity of various z0 ‘s.
function M=hw3p6o(zMax,c,iteration,resolution)
    clear i;clc;
    % i. In this function, you first want to make a 500x500 matrix that contains complex numbers with real part between –zMax and zMax, and imaginary part between –zMax and zMax.
        % Call this matrix Z. Make the imaginary part vary along the y axis of this matrix. You can most easily do this by using linspace and meshgrid, but you can also do it with a loop.
    if nargin
    re=linspace(-zMax,zMax,resolution);
    im=1i*re;
    [Re,Im]=meshgrid(re,im);
    Z=Re+Im;
    % ii. For each element of Z, compute the escape velocity (by calling your escapeVelocity) and store it in the same location in a matrix M.
        % When done, the matrix M should be the same size as Z and contain escape velocities with values between 1 and N.
    M=zeros(size(Z)); %allocate M
    for x=1:size(Z,1)
        for y=1:size(Z,2)
            M(x,y)=escapeVelocity(Z(x,y),c,iteration);
        end
    end

    % iii. Run your julia function with various zMax, c, and N values to generate various fractals. To display the fractal nicely, use imagesc to visualize atan(0.1*M),
        % (taking the arctangent of M makes the image look nicer; you may also want to use axis xy so the y values aren’t flipped). WARNING: this function may take a while to run.
    close all;
    figure;
    imagesc(atan(0.1*M));
    xlabel('Re(z)');
    ylabel('Im(z)');
    axis xy;
    colormap(jet);
end

References

Chance, E. M., Clore, G. M., Curtis, A. R., & Shephard, E. P. (1981). Numerical solution of the Hodgkin-Huxley equations in a moving coordinate system. Simulation of nerve impulse transmission over long distances. Journal of Computational Physics, 40(2), 318-326.

Research: Practice on Matlab Day 2 (MIT 6.094)

It’s Tuesday! although technically I’m uploading this on Wednesday.

So Lecture 1 + Homework 1 was refreshingly full of information and good practices, it was worth a good day’s work. Continuing from yesterday. Today I finished Lecture 2 (MIT OCW) and the mandatory part of Homework 2 (HW). Which is perfect because this is the second day of the week.

Today a major chunk of my time was spent researching about cars, reasonable prices, and my lost international driving permit. Big time waster. Well I thought the deal was on Wednesday, but then we agreed to move the car inspection to Friday afternoon.

For my own documentation purposes, and I’ll upload my answers of Homework 2

I might add the optional homework later on, all the question is embedded in the comments. I’ll also add more remarks on how my thought process goes on

IAP January 2010 — MIT 6.094 — Assignment 2

Problem 1

Simple warm-up questions designed to remind me of homework 1

% Homework 2 - 1 Semilog Plot
% Semilog plot. Over the past 5 years, the number of students in 6.094 has been 15, 25, 55, 115,
    % 144. Class size seems like it’s growing exponentially. To verify this, plot these values on a plot
    % with a log y scale and label it (semilogy, xlabel, ylabel, title). Use magenta square symbols of
    % marker size 10 and line width 4, and no line connecting them. You may have to change the x
    % limits to see all 5 symbols (xlim). If the relationship really is exponential, it will look linear on a
    % log plot.

student=[15,25,55,115,144];
semilogy(student,'ms','MarkerSize',10,'LineWidth',4);
title('Student enrollment in 6.094 growth');
xlabel('Years');
ylabel('Students');
%it really does look linear, so it is exponential!

Problem 3

Very easy – An exercise in looking over special graph documentation

% Homework 2 - 3 Bar Graph, problem 2 is skipped due to copyright issues
% Bar graph. Make a vector of 5 random values and plot them on a bar graph using red bars,
doc specgraph 
bar(rand(1,5),'r');

Problem 4

Using the interpolation functions of Matlab with advanced data visualization

% Homework 2 - 4 Interpolation and Surface
% Interpolation and surface plots. Write a script called randomSurface.m to do the following
    clear all;
    clc;
    close all;
    % a. To make a random surface, make Z0 a 5x5 matrix of random values on the range [0,1]
        % (rand).
    Z0=rand(5);
    % b. Make an X0 and Y0 using meshgrid and the vector 1:5 (use the same vector for both
        % inputs into meshgrid). Now, X0, Y0, and Z0 define 25 points on a surface.
    [X0,Y0]=meshgrid(1:5,1:5);
    % c. We are going to interpolate intermediate values to make the surface seem smooth.
        % Make X1 and Y1 using meshgrid and the vector 1:.1:5 (again use the same vector for
        % both inputs into meshgrid).
    [X1,Y1]=meshgrid(1:0.1:5,1:0.1:5);
    % d. Make Z1 by interpolating X0, Y0, and Z0 at the positions in X1 and Y1 using cubic
    % interpolation (interp2, specify cubic as the interpolation method).
    Z1=interp2(X0,Y0,Z0,X1,Y1,'cubic');
    % e. Plot a surface plot of Z1. Set the colormap to hsv and the shading property to interp
    % (surf, colormap, shading).
    surf(X1,Y1,Z1);
    colormap(hsv);
    shading interp;
    % f. Hold on to the axes and plot the 15-line contour on the same axes (contour).
    hold on;
    contour(X1,Y1,Z1);
    % g. Add a colorbar (colorbar).
    colorbar;
    % h. Set the color axis to be from 0 to 1 (caxis).
    caxis([0,1]);

Problem 5

Finding a value/index closest to the desired value requires only 1 line of code and 2 functions. Getting all equally closest values/indices would require a loop, usually, but not on Matlab. The find() function saves the day.

% Homework 2 - 5 Fun with Find
% Fun with find. Write a function to return the index of the value that is nearest to a desired
    % value. The function declaration should be: ind=findNearest(x, desiredVal). x is a
    % vector or matrix of values, and desiredVal is the scalar value you want to find. Do not
    % assume that desiredVal exists in x, rather find the value that is closest to desiredVal. If
    % multiple values are the same distance from desiredVal, return all of their indices. Test your
    % function to make sure it works on a few vectors and matrices. Useful functions are abs, min, and
    % find. Hint: You may have some trouble using min when x is a matrix. To convert a matrix Q into
    % a vector you can do something like y=Q(:). Then, doing m=min(y) will give you the
    % minimum value in Q. To find where this minimum occurs in Q, do inds=find(Q==m);.

function ind=findNearest(x, desiredVal)
    %first attempt at converting to vector
    %x=reshape(x,1,prod(size(x)));
    %second attempt after reading the question
        %make the matrix into a vector and calculate the distance
    x=abs(x(:)-desiredVal);
    %find the minimum distance
    m=min(x);
    %find all index that has the minimum distance calculated
    ind=find(x==m)';
end

Problem 6

Just a simple practice on else if structure.

% Homework 2 - 6 Loops and flow control
% Loops and flow control. Make function called loopTest(N) that loops through the values 1
    % through N and for each number n it should display ‘n is divisible by 2’, ‘n is divisible by 3’, ‘n is
    % divisible by 2 AND 3’ or ‘n is NOT divisible by 2 or 3’. Use a for loop, the function mod or rem to
    % figure out if a number is divisible by 2 or 3, and num2str to convert each number to a string for
    % displaying. You can use any combination of if, else, and elseif.
function loopTest(N)
    for n=1:N
        if mod(n,6)==0
            disp([num2str(n) ' is divisible by 2 AND 3']);
        elseif mod(n,2)==0
            disp([num2str(n) ' is divisible by 2']);
        elseif mod(n,3)==0
            disp([num2str(n) ' is divisible by 3']);
        else
            disp([num2str(n) ' is NOT divisible by 2 or 3']);
    end
end

Problem 7 – fun

This one was fun because I had to actually think, on how to smooth the kinks at the end, and how convolution actually works. After thinking about the question ‘… correct edge effects so that the smoothed function doesn’t deviate from the original at the start or the end…‘ I finally get it. They wanted the averaged/convolved data points to stay the same at the start or end stays the same which can only mean there is a buffer as large as half the width on either side. I’ll put both before and after plot.

Left is the result of blindly convolving, right is after a ‘buffer’ is added.

% Homework 2 - 7 Smoothing filter
% Smoothing filter. Although it is a really useful function, Matlab does not contain an easy to use
% smoothing filter. Write a function with the declaration: smoothed=rectFilt(x,width).
% The filter should take a vector of noisy data (x) and smooth it by doing a symmetric moving
% average with a window of the specified width. For example if width=5, then smoothed(n)
% should equal mean(x(n-2:n+2)). Note that you may have problems around the edges:
% when n<3 and n>length(x)-2.
function smoothed=rectFilt(x,width)
    close all;
    % a. The lengths of x and smoothed should be equal.
    % b. For symmetry to work, make sure that width is odd. If it isn’t, increase it by 1 to make
        % it odd and display a warning, but still do the smoothing.
    if mod(width,2)==0
        disp('warning: width is increased by one for symmetry');
        width=width+1;
    end
    % c. Make sure you correct edge effects so that the smoothed function doesn’t deviate from
        % the original at the start or the end. Also make sure you don’t have any horizontal offset
        % between the smoothed function and the original (since we’re using a symmetric moving
        % average, the smoothed values should lie on top of the original data).
    half=(width-1)/2;
    xBuffer=[x(1)*ones(1,half) x x(end)*ones(1,half)];
    % d. You can do this using a loop and mean (which should be easy but may be slow), or more
        % efficiently by using conv (if you are familiar with convolution).
    k=ones(1,width)/width;
    xConvolved=conv(xBuffer,k,'same');
    smoothed=xConvolved(half+1:end-half);

    %smoothed=conv(x,k,'same'); %this is the 'before' version with kinks at
        %the edge

    % e. Load the mat file called noisyData.mat. It contains a variable x which is a noisy line. Plot
        % the noisy data as well as your smoothed version, like below (a width of 11 was used):
    figure;
    plot(x,'.');
    hold on;
    plot(smoothed,'r-')

end

Research: Practice on Matlab Day 1 (MIT 6.094)

Since I had to research water desalination and make a model using Matlab-Simulink. Today I decided to learn Matlab right and take on all the assignment from the online IAP course MIT 6.094.

This course is perfect for me as it covers all the essential topics including Simulink.

The first lecture (MIT OCW) is concise and easily understandable, unfortunately there is no solution available for the first assignment. Although the first homework (HW1) was simple enough that it is easy to check the correctness.

For my own documentation purposes, and I’ll upload my answers of Homework 1:

IAP January 2010 — MIT 6.094 — Assignment 1
calculateGrades.m
load classGrades;
totalStudent = size(namesAndGrades,1);

%extract only the grades, col 2 to col end
grades = namesAndGrades(:,2:end);
disp(‘original grades’);
disp(grades);
%take a assignment-wise mean
meanAssignment = nanmean(grades,1);
disp(‘mean grade before curve’);
disp(meanAssignment);
%normalize mean grade to 3.5
meanMatrix = ones(totalStudent,1)*meanAssignment;
%disp(‘meanMatrix’);
%disp(meanMatrix);
curvedGrades=(3.5*grades./meanMatrix);
%mean of curved grade
meanCurved = nanmean(curvedGrades,1);
disp(‘mean grade after curve’);
disp(meanCurved);

%fix values above 5, the A cutoff
curvedGrades(find(curvedGrades>5))=5;
%curved result
disp(‘curvedGrades’);
disp(curvedGrades);
%letter grading
letters=(‘FDCBA’);
studentGrade = ceil(nanmean(curvedGrades,2));
disp(‘studentGrade’);
disp(studentGrade);
studentGradeLetter=letters(studentGrade);
disp(‘studentGradeLetter’);
disp(studentGradeLetter);

seriesConvergence.m
clear all;
clc;

p=0.99;
k=0:1:1000;
geomSeries=p.^k;
%disp(geomSeries(1:10:100));
G=(1-p)^(-1);
disp([‘G=’ num2str(G)]);
figure;
plot(k,G+0*k,’r’,k,cumsum(geomSeries)); %red line at the max k
hold on;
%put limits on the plot
xlim([0,1000]);
ylim([0,150]);
%label the k and function
xlabel(‘index’);
ylabel(‘sum’);
title([‘convergence of infinire sum with p=’ num2str(p)]);
legend(‘Infinite sum’,’finite sum’);

%——————————————————-

p=2;
n=1:1:500;
pSeries=n.^(-p);
disp(pSeries(1:10:100));
P=pi^2/6;
%disp([‘P=’ num2str(P)]);
figure;
plot(n,P+0*n,’r’,n,cumsum(pSeries)); %red line at the max k
hold on;
%put limits on the plot
xlim([0,500]);
ylim([1.2,1.8]);
%label the k and function
xlabel(‘index’);
ylabel(‘sum’);
title([‘convergence of infinite sum with p=’ num2str(p)]);
legend(‘Infinite sum’,’finite sum’);

throwBall.m
%THROWBALL
clc;
%b. At the top of the file, define some constants (you can pick your own variable names)
h=1.5; %[m] initial height
g=9.8; %[m/s^2] gravitational acceleration
v=4; %[m/s] velocity of ball release
dir=45; %[degree] angle of release

%c. Next, make a time vector that has 1000 linearly spaced values between 0 and 1, inclusive.
t=linspace(0,1,1000);

%d. If x is distance and y is height, the equations below describe their dependence on
%time and all the other parameters (initial height h , gravitational acceleration g , initial
%ball velocity v , angle of velocity vector in degrees θ ). Solve for x and y
rad=pi/180*dir; %convert deg to rad
x=v*cos(rad)*t; %x position overtime
y=h+v*sin(rad)*t-0.5*g*t.^2; %y position overtime

%e. Approximate when the ball hits the ground.
%i. Find the index when the height first becomes negative (use find).
[error,groundInd]=min(abs(y)); %deviation from instruction: I minimized the function to be close to 0 instead
%ii. The distance at which the ball hits the ground is value of x at that index
xGround=x(groundInd);
%iii. Display the words: The ball hits the ground at a distance of X meters. (where X is
%the distance you found in part ii above)
disp([‘The ball hits the ground at a distance of ‘ num2str(xGround) ‘ meters.’]);

%f. Plot the ball’s trajectory
%i. Open a new figure (use figure)
figure;
%ii. Plot the ball’s height on the y axis and the distance on the x axis (plot)
plot(x,y);
%iii. Label the axes meaningfully and give the figure a title (use xlabel, ylabel, and
%title)
xlabel(‘distance’);
ylabel(‘height’);
title(‘parametric curve of ball trajectory’);
%iv. Hold on to the figure (use hold on)
hold on;
%v. Plot the ground as a dashed black line. This should be a horizontal line going
%from 0 to the maximum value of x (use max). The height of this line should be
%0. (see help plot for line colors and styles)
plot(x,0+0*x,’k–‘);
ylim([-0.5,ceil(max(y))]);

encrypt.m
%ENCRYPTION
clc;
%b. At the top of the script, define the original string to be: This is my top secret message!
msg=’This is my top secret message!’;
%c. Next, let’s shuffle the indices of the letters. To do this, we need to make a string of encoding indices
%i. Make a vector that has the indices from 1 to the length of the original string in a
…randomly permuted order. Use randperm and length
randInd=randperm(length(msg));
%ii. Encode the original string by using your encoding vector as indices into
…original . Name the encoded message encoded .
encoded=msg(randInd);
%d. Now, we need to figure out the decoding key to match the encoding key we just made.
%i. Assemble a temporary matrix where the first column is the encoding vector you
…made in the previous part and the second column are the integers from 1 to the
…length of the original string in order. Use length, and you may need to transpose
…some vectors to make them columns using ‘.
tempMat=[transpose(randInd),transpose(1:1:length(msg))];
%ii. Next, we want to sort the rows of this temporary matrix according to the values
…in the first column. Use sortrows.
tempMat=sortrows(tempMat,1);
%iii. After it’s been sorted, extract the second column of the temporary matrix. This
…is your decoding vector.
decodeInd=tempMat(:,2);
%iv. To make the decoded message, use the decoding vector as indices into
…encoded .
decoded=encoded(decodeInd);
%e. Display the original, encoded, and decoded messages
%i. Display the following three strings, where : original , encoded , and decoded
…are the strings you made above. Use disp
%Original: original
disp(‘original:’);disp(msg);
%Encoded: encoded
disp(‘encoded:’);disp(encoded);
%Decoded: decoded
disp(‘decoded:’);disp(decoded);

%f. Compare the original and decoded strings to make sure they’re identical and display the
…result
%i. Use strcmp to compare the original and decoded strings. Name the output
…of this operation correct . correct will have the value 1 if the strings match
…and the value 0 if they don’t
correct=strcmp(msg,decoded);

%ii. Display the following string: Decoded correctly (1 true, 0 false): correct
…use disp and num2str
disp([‘Decoded correctly (1 true, 0 false):’ num2str(correct)]);
%debugger

Starting to Code

My task working as undergraduate research assistant in Energy Systems Lab is to help program an Arduino board so it can communicate effectively with an engine control unit and battery control.

At first I was very uncomfortable and felt at a loss at the task in hand. As I’ve never considered programming to be what I can do well. But I can’t refuse such an interesting and hard challenge. So I mainly stared at the code given to me wondering what to do with it. This is quite unlike anything I’ve done before and it takes a jump in ability. Although it does help that I had programming course 2 years ago.

So I decided to revisit the fundamentals to help for me to brush up and practice.

Well I know there’s a community that loves coding and Arduino in my campus, the Tech Club so that’s a big resource. I met and chatted with the current president, Adam, and he gave me tips on just do Python for PC tasks and C++ instead of learning C for Arduino. I also joined the club to learn more from people and ask questions. That weekend I spent the whole time learning Python syntax and usage from Codecademy up to Functions. It worked wonderfully because Python is so user friendly, I created my first program to write and read to a serial communication by using Python. It was a huge boost in motivation.

But next tasks wasn’t quite so easy as I had to create stuff in Arduino language (which is a subset of C++).  So after googling I found Accelerated C++ by Andrew Koenig always came up and had a very good review as being easy to understand and gets me up to speed quickly. Every weekend I sit with it and did all the exercises. It does wonders on making me feel more at ease, and also I have better understanding about what is going on in a code.

Once again I am reminded that ordered practice (i.e. doing vs reading) is much more useful for understanding than trying to read various random snippets online and just skimming a book. Though there is a time-trade off.

Accelerated C++ – Andrew Koenig – Chapter 4

I practiced CPP today with Chapter 4 – Working with Batches of data:

This is my solution, they worked in Code::Block following  C++14 convention.

 

//Exercise 4-1
    string::size_type maxlen;
    StudentInfo s;
    max(s.name.size(), maxlen);


  //Exercise 4-2
    for(int i = 0; i!=101; i++){
        std::cout << std::setw(3) << i << std::setw(6) << i*i << std::endl;
    }


  //Exercise 4-3
    int n = 1000;
    std::streamsize width = (unsigned int)(log10(n)) + 1;

    for(int i = 0; i!= n+1; i++){
        std::cout << std::setw(width) << i << std::setw(width*2) << i*i << std::endl;
    }

    std::setw(0);


  //Exercise 4-4
    int n = 1000;
    std::streamsize digits = (unsigned int)(log10(n)) + 1;
    std::streamsize p = digits*2;
    std::streamsize originalPrec = std::cout.precision();
    std::streamsize originalWidth = std::cout.width();

    for(int i = 0; i!= n+1; i++){
        double d = i;
        std::cout << std::setprecision(p) << std::setw(digits) << d
            << std::setw(digits*2) << d*d <<std::setw(originalWidth)
            << std::setprecision(originalPrec) << std::endl;
    }


  //Exercise 4-5
    struct Words{
    std::string name;
    int tally ;
    };

    using vecword_t = std::vector<Words>;
    using vecsz_t   = std::vector<Words>::size_type;
    using strsz_t   = std::string::size_type;

    std::istream& readTally(std::istream& in, vecword_t word){
        if(in){
            //empties word vector
            word.clear();

            std::string wordInput;
            strsz_t maxlen = 0;

            while(in >> wordInput){
                //initialize duplicate flag
                bool duplicate = 0;
                //get the longest word
                maxlen = std::max(maxlen, wordInput.size());
                //check duplicate
                for(vecsz_t i = 0; i != word.size(); i++){
                    if(wordInput == word[i].name){
                        word[i].tally++;

                        duplicate = 1;
                        break;
                    }
                }

                //if no duplicates found add vector
                if(duplicate == 0){
                    Words wordPush = {.name = wordInput, .tally = 1};
                    word.push_back(wordPush);
                }
            }
            //clear error state in istream
            in.clear();
            //print the result
            std::cout << "Word count" << std::endl;
            for(vecsz_t i = 0; i != word.size(); i++){
                std::cout << std::setw(maxlen + 1) << word[i].name
                    << ' ' << std::setw(4) << word[i].tally
                    << std::setw(0) << std::endl;
            }
        }
        return in;
    }
    std::istream& readTally(std::istream& in){
        vecword_t word = {};
        return readTally(in,word);
    }


  //Exercise 4-6
	struct Student_info{
		std::string name;
		double grade;
	};

	bool compare(const Student_info& x, const Student_info& y){
		return x.name < y.name;
	}

	std::istream& readHw(std::istream& , std::vector<double>& );
	double medianVec(std::vector<double>);
	double gradeCalc(double midGrade, double finGrade, const std::vector<double>& hwGradeList);

	std::istream& read(std::istream& is, Student_info& s){
		double mid, fin;
		std::vector<double> homeWork;
		is >> s.name >> mid >> fin;
		readHw(is, homeWork);

		s.grade = gradeCalc(mid,fin,homeWork);

		return is;

	}

	void printPrecision(double inputNum, std::streamsize sigFig){
		std::streamsize prec = std::cout.precision();
		std::cout << std::setprecision(sigFig) << inputNum
			<< std::setprecision(prec) << std::endl;
	}

	int main()
	{
		Student_info buffer;

		std::vector<Student_info> students;
		std::string::size_type maxlen = 0;

		//get input and final grade
		while(read(std::cin, buffer)){
				maxlen = std::max(maxlen, buffer.name.size());
				students.push_back(buffer);
		}

		std::cin.clear();

		//alphabetize student list
		std::sort(students.begin(),students.end(),compare);

		std::cout << "The result is ordered in alphabetical order.." << std::endl;
		// calculate grade
		for(unsigned int i=0; i != students.size(); ++i){
			std::cout << students[i].name << std::string(maxlen + 1 - students[i].name.size(), ' ');
			printPrecision(students[i].grade,3);
		}
		return 0;
	}



  //Exercise 4-7
    typedef std::vector<double>::size_type vecsz_t;
    std::vector<double> vecNum = {1.0, 1.0, 1.0, 2.0};

    double sum = 0;

    for(vecsz_t i = 0; i != vecNum.size(); i++ ){
        sum += vecNum[i];
        std::cout << vecNum[i] <<std::endl;
    }


    double avg = sum/(double)vecNum.size();
    std::cout<<"sum: " <<sum <<" avg: " << avg;



  //Exercise 4-8
    //if statement below is legal then f has return type of std::vector<double>
    double d = f()[n];