|
Thread Rules 1. This is not a "do my homework for me" thread. If you have specific questions, ask, but don't post an assignment or homework problem and expect an exact solution. 2. No recruiting for your cockamamie projects (you won't replace facebook with 3 dudes you found on the internet and $20) 3. If you can't articulate why a language is bad, don't start slinging shit about it. Just remember that nothing is worse than making CSS IE6 compatible. 4. Use [code] tags to format code blocks. |
I started with JS: The Good Parts. I'm not an expert though, so I don't know what's a good place to start.
Definitely don't start with AngularJS though IMO, I feel like it's pretty confusing in the beginning, especially if you haven't learned JavaScript yet.
|
|
On September 30 2016 14:46 Aerisky wrote: Yeah, normally I've applied to a ton of places for internships (100% agree that it's a numbers game). Now that I'm a soon-to-be new grad, though, I have a return offer from the company at which I interned this summer; I really like the offer so I've only been applying to really cool places where I'd be interested in working/prestigious places.
Unless you are 100% sure that the company at which you interned will hire you, don't change your application strategy. I did this mistake. I was an intern in one company for a year, even did the research for my Master's thesis there, and I got accepted for the next wave of their graduate program after I would graduate (I still had one year of university after finishing the internship) or alternatively I could start a PhD at one of the universities they had an agreement with.
Then when I was back at the university - a month or so before graduation last Summer - I was told by my previous supervisor that the new management decided to discontinue both programs... Never put all your eggs in one basket. If you get accepted by more companies, you can always decline.
|
I'm running into more and more issues with Ubuntu Mate on Raspberry Pi 3 and it's a real pain to troubleshoot a new issue every day.
Does anyone know if Raspian is an adequate OS to code on? I'm doing a web development programme called The Odin Project and have to learn imperative programming in Pascal (...) for a course at uni.
|
On September 30 2016 20:19 DickMcFanny wrote: I'm running into more and more issues with Ubuntu Mate on Raspberry Pi 3 and it's a real pain to troubleshoot a new issue every day.
Does anyone know if Raspian is an adequate OS to code on? I'm doing a web development programme called The Odin Project and have to learn imperative programming in Pascal (...) for a course at uni. Raspbian is an adequate OS to code on. I built a LAMP stack on Raspbian and Raspberry Pi 2 with zero experience on any part of the LAMP stack in an hour. It is 95% compatible with Ubuntu too.
|
On September 30 2016 20:19 DickMcFanny wrote: I'm running into more and more issues with Ubuntu Mate on Raspberry Pi 3 and it's a real pain to troubleshoot a new issue every day.
Does anyone know if Raspian is an adequate OS to code on? I'm doing a web development programme called The Odin Project and have to learn imperative programming in Pascal (...) for a course at uni. I've run a Pi 2 as a web server with Raspian and it worked fine. Getting scientific python libraries to work was a real bitch, but other than that, it was all pretty easy. Apache and nginx both work. I didn't have problems with node or php either.
Whether it's a good platform to actually code on? I wouldn't really, but as long as your favourite IDE works, it should be okay. Vim works fine, as does emacs.
Also, what the hell uni still teaches programming in pascal?! Please change to a uni that isn't stuck in the 90s.
|
.NET file moving and permissions question.
This may or may not be a silly question, but I just want to make sure there isn't something I haven't considered. I have a network share that can only be accessed by a windows service account. I have a piece of software that at times will need to copy a file to the network share. The regular user that is using my software won't have access to the share, but I can use impersonation to impersonate the service account to do it (http://stackoverflow.com/questions/125341/how-do-you-do-impersonation-in-net).
The problem is that the service account doesn't have access to the user's documents directory. That means I can't do a file copy as the user because the user doesn't have access to the share, nor can I do a file copy as the windows service account because the account doesn't have access to the source file.
Is there any way to copy a file to a network share the current logged in user may not have access to using a service account? It seems like there should be since in Windows you can map a network share as a different user.
Edit: I did a bit more googling and was able to find this: http://stackoverflow.com/questions/8999373/passing-data-between-threads-in-c-sharp. I think I should be able to use that to do what I need to do.
|
On September 27 2016 22:57 maybenexttime wrote:Let me first introduce myself briefly. I am a mechanical engineering graduate who'd like to learn some Python and MATLAB programming, because I consider programming an important skill in every engineer's repertoire. For the past two months I have studied both languages and worked on some private projects of mine since that seems to be the best way to learn useful things. I'm currently working on a MATLAB program that is supposed to quantify the degree of banding in the microstructures I simulated with an earlier script (the output is basically a set of (x,y) coordinates), and I've encountered some performance issues when compared with the program whose analysis I am trying to recreate (PASSaGE, link below). My understanding of how the analysis is supposed to work is based on the paper written by the software's author (link) and the program's manual, but some sections were quite misleading in my opinion, so I had to do several modifications based on some experiments with analysis. http://www.passagesoftware.net/http://rosenberglab.net/Pubs/JVegSci2004v15p277.pdfInitially, my MATLAB script took roughly 15 minutes (more specifically, the function that runs the analysis, not the overarching script) to analyze a set of 1000 points, whereas PASSaGE takes around 1-2 minutes. It relied mostly on loops. After a few days of optimizing, I managed to decrease the calculation time to 160-170 seconds, but I seem to have hit a dead end. Here's the overarching script (just in case): + Show Spoiler +clc clear
%%% PARAMETERS %%%
answer = input('Should the default value be used for the edge region parameter? Y/N [Y]: ','s'); if isempty(answer) % If the user simply presses 'Enter', the answer is assumed to be 'Y'. answer = 'Y'; end
disp(' ') % Empty line for clarity.
if strcmp(answer, 'Y') == 1 edge_parameter = 0.25; disp('The default value of the edge region parameter (0.25) will be used.') else edge_parameter = input('Please, specify the value of the edge region parameter: '); end
%%% IMPORT OF DATA %%%
filename = 'banding_coordinates.txt'; delimiterIn = ' '; headerlinesIn = 1; D = importdata(filename,delimiterIn,headerlinesIn);
%%% EXTRACTION OF COORDINATES %%%
x_extracted = D.data(:,1); y_extracted = D.data(:,2);
% Transposition into vectors:
x = transpose(x_extracted); y = transpose(y_extracted);
% Verification of vector lengths:
if length(x) ~= length(y) disp('The imported data is corrupted: the numbers of the x and y coordinates do not match') end
%%% BOUNDARY CONDITIONS %%%
% Limits of the analyzed space:
x_min = min(x); x_max = max(x); y_min = min(y); y_max = max(y);
% Edge region parameters for each coordinate:
edge_parameter_x = edge_parameter*(x_max - x_min); edge_parameter_y = edge_parameter*(y_max - y_min);
% Limits of the focal point space:
x_foc_min = x_min + edge_parameter_x; x_foc_max = x_max - edge_parameter_x; y_foc_min = y_min + edge_parameter_y; y_foc_max = y_max - edge_parameter_y;
%%% EXTRACTION OF FOCAL POINTS %%%
% Generation of focal points:
[in,on] = inpolygon(x, y, [x_foc_min, x_foc_max], [y_foc_min, y_foc_max] ;
x_foc = x(in); y_foc = y(in);
% Focal points plot:
answer = input('Would you like to display the generated focal points? Y/N [N]: ','s'); if isempty(answer) % If the user simply presses 'Enter', the answer is assumed to be 'Y'. answer = 'N'; end
disp(' ') % Empty line for clarity.
if strcmp(answer, 'Y') == 1 figure
scatter(x,y,'.') % all points axis equal
hold on scatter(x_foc,y_foc,'r+') % focal points hold off end
%%% ANALYSIS %%% tic overall_variance = calculate_variance_test6_3(x, y, x_foc, y_foc, ... x_min, x_max, y_min, y_max); toc %%% OUTPUT %%%
% Original variance peak:
disp('Original variance peak:') max(overall_variance)
% Plot:
figure
plot(1:1:180, overall_variance) And here's the function that does the actual analysis (in its most optimal form): + Show Spoiler +function variance_avg_foc = calculate_variance_test6_3(x, y, x_foc, y_foc, ... x_min, x_max, y_min, y_max)
% Calculation of the number of focal points:
m = length(x_foc); % Number of focal points = number of iterations in the first loop.
% Calculation of the number of scales:
scale = 1:1:45; % All scales. l = length(scale);
% Preallocation of memory for all focal points:
variance_avg_scale = zeros(m, 180); % Calculation of the number of analyzed points: n = length(x) - 1; % Calculation of transect properties:
p = 1:1:180; transect_angles = p - 1; slopes = tand((transect_angles)); % Selection of the focal point
for j = 1:m % Removal of the focal point from the pool of analyzed points: index = find( (x(1, == x_foc(j)) & (y(1, == y_foc(j)) ); if length(index) ~= 1 % If there is more than once point with the x_foc(j), y_foc(j) index = index(1); % coordinates, only one point will be subtracted from the future end % local coordinates pool. x_without_foc = x(1:end ~= index); y_without_foc = y(1:end ~= index); % Creation of local coordinates /vectorized:
i = 1:1:n; % One iteration for each (x, y) point. x_local = x_without_foc(i) - x_foc(j); y_local = y_without_foc(i) - y_foc(j);
% Limits of the studied area in local coordinates: x_min_local = x_min - x_foc(j); x_max_local = x_max - x_foc(j); y_min_local = y_min - y_foc(j); y_max_local = y_max - y_foc(j); % Translation to local polar coordinates:
[angle,radius] = cart_to_pol(x_local,y_local); angle = rad2deg(angle); % Consolidation of coordinates data:
data = [angle; radius];
% Preallocation of memory for all scales:
transect_area = zeros(1, 180); % Calculation of the transect area: for p = 1:180 transect_angle = transect_angles(p); slope = slopes(p); transect_area(p) = calculate_transect_area(x_min_local, ... x_max_local, y_min_local, y_max_local, transect_angle, slope); end % Preallocation of memory for all scales:
variance_scale = zeros(l, 180);
% Selection of scale:
for k = 1:l
% Transect width: transect_width = scale; % Preallocation of memory for all positions:
wavelet_transform = zeros(1, 180); variance_normalized = zeros(1, 180);
% Selection of angular position:
for p = 1:180 % Finding the points inside the specified transect:
points_observed = observe_points(data, p, transect_width(k));
% Extraction of the polar coordinates of the observed points:
angle_observed = points_observed(1, ; radius_observed = points_observed(2, ;
% Calculation of the number of observed points:
observations = length(angle_observed); % Calculation of scaled wavelets for each point within the transect /vectorized:
if observations == 0 scaled_wavelet = 0; else o = 1:1:observations; t = (angle_observed(o) - transect_angles(p))./scale(k); scaled_wavelet = radius_observed.*wavelet_function(t); end wavelet_transform(p) = sum(scaled_wavelet)/scale(k); variance_normalized(p) = wavelet_transform(p)^2/transect_area(p); end variance_scale(k, = variance_normalized; end variance_avg_scale(j, = sum(variance_scale)./l; % All specified scales. end variance_avg_foc = sum(variance_avg_scale)./m; end
%-------------------------------------------------------------------------- function [angle, radius] = cart_to_pol(x, y)
angle = arctan(y, x); radius = hypot(x,y);
end
%-------------------------------------------------------------------------- function angle = arctan(y, x) n = length(x); angle = zeros(1, n); parfor i = 1:n % /vectorized angle(i) = atan2(y(i), x(i)); if angle(i) < 0 angle(i) = angle(i) + pi; end end end
%-------------------------------------------------------------------------- function points_observed = observe_points(data, p, transect_width)
% Finding the points inside the specified transect:
step_size = 1;
log_ind = ( (p - 1 - (transect_width/2))*step_size <= data(1, & ... data(1, <= (p - 1 + (transect_width/2))*step_size ) | ... ( ((p - 1 - (transect_width/2))*step_size + 180) <= data(1, & ... data(1, <= ((p - 1 + (transect_width/2))*step_size + 180) ); points_observed = data(:,log_ind);
end
%-------------------------------------------------------------------------- function transect_area = calculate_transect_area(x_min_local, ... x_max_local, y_min_local, y_max_local, transect_angle, slope) % Slope variants: if (transect_angle == 0) || (transect_angle == 180) % Horizontal line. x_lim_1 = x_min_local; x_lim_2 = x_max_local; y_lim_1 = 0; y_lim_2 = 0; elseif transect_angle == 90 % Vertical line. x_lim_1 = 0; x_lim_2 = 0; y_lim_1 = y_min_local; y_lim_2 = y_max_local; elseif transect_angle > 0 && transect_angle < 90 % Boundary coordinates:
if y_min_local/slope >= x_min_local x_lim_1 = y_min_local/slope; y_lim_1 = y_min_local; else x_lim_1 = x_min_local; y_lim_1 = slope*x_min_local; end
if y_max_local/slope <= x_max_local x_lim_2 = y_max_local/slope; y_lim_2 = y_max_local; else x_lim_2 = x_max_local; y_lim_2 = slope*x_max_local; end else
% Boundary coordinates:
if y_min_local/slope >= x_max_local x_lim_1 = x_max_local; y_lim_1 = slope*x_max_local; else x_lim_1 = y_min_local/slope; y_lim_1 = y_min_local; end
if y_max_local/slope <= x_min_local x_lim_2 = x_min_local; y_lim_2 = slope*x_min_local; else x_lim_2 = y_max_local/slope; y_lim_2 = y_max_local; end end transect_area = (x_lim_1)^2 + (y_lim_1)^2 + (x_lim_2)^2 + (y_lim_2)^2; % This is not real area of the transect. % It's area divided by omega/2, but that holds % for all transects and omega is constant. end
%-------------------------------------------------------------------------- function mexican_hat = wavelet_function(t) % /vectorized
mexican_hat = 2./(sqrt(3)).*pi.^(-1/4).*(1 - 4.*t.^2).*exp(-2.*t.^2); % Wavelet function normalized to have unit energy. % Vectorizing with "vectorize()" does not work. % It cannot tell if I want to use ^ for matrices (mpower) % or scalars (power). end According to the MATLAB profiler, there are two lines/functions whose vectorization could yield significant improvements in calculation time, but I failed at vectorizing them properly... Any help would be appreciated. This line currently requires roughly 95 seconds (it's being called 4 million times): + Show Spoiler +points_observed = observe_points(data, p, transect_width(k)); And these two lines require roughly 30 seconds (they're being called 700 thousand times): + Show Spoiler +t = (angle_observed(o) - transect_angles(p))./scale(k); scaled_wavelet = radius_observed.*wavelet_function(t); Additionally, I am looking for a way to make scaled_wavelet call wavelet_function only for non-zero values of radius_observed. I tried the solution below, but it only made the calculations longer: + Show Spoiler +% Calculation of scaled wavelets for each point within the transect /vectorized:
if observations == 0 scaled_wavelet = 0; else log_ind = (radius_observed ~= 0); %o = 1:1:observations; scaled_wavelet = radius_observed(log_ind).*wavelet_function((angle_observed(log_ind) - transect_angles(p))./scale(k)); end I also found that spfun() does something similar, but I don't know how to make it work with a custom function like wavelet_function().
So could anyone lend me a helping hand with this?
|
Sorry, haven't touched MATLAB since EE like 8+ yrs ago.
On September 30 2016 14:46 Aerisky wrote: Yeah, normally I've applied to a ton of places for internships (100% agree that it's a numbers game). Now that I'm a soon-to-be new grad, though, I have a return offer from the company at which I interned this summer; I really like the offer so I've only been applying to really cool places where I'd be interested in working/prestigious places. Fair enough, but there's enough places that are prestigious or what have you that you can still shotgun approach it. Depends on how good your offer is too. If it looks like it'd be comparable to e.g. making >200k/yr after a few years at BigCo, then w/e.
|
On October 01 2016 00:40 maybenexttime wrote:Show nested quote +On September 27 2016 22:57 maybenexttime wrote:Let me first introduce myself briefly. I am a mechanical engineering graduate who'd like to learn some Python and MATLAB programming, because I consider programming an important skill in every engineer's repertoire. For the past two months I have studied both languages and worked on some private projects of mine since that seems to be the best way to learn useful things. I'm currently working on a MATLAB program that is supposed to quantify the degree of banding in the microstructures I simulated with an earlier script (the output is basically a set of (x,y) coordinates), and I've encountered some performance issues when compared with the program whose analysis I am trying to recreate (PASSaGE, link below). My understanding of how the analysis is supposed to work is based on the paper written by the software's author (link) and the program's manual, but some sections were quite misleading in my opinion, so I had to do several modifications based on some experiments with analysis. http://www.passagesoftware.net/http://rosenberglab.net/Pubs/JVegSci2004v15p277.pdfInitially, my MATLAB script took roughly 15 minutes (more specifically, the function that runs the analysis, not the overarching script) to analyze a set of 1000 points, whereas PASSaGE takes around 1-2 minutes. It relied mostly on loops. After a few days of optimizing, I managed to decrease the calculation time to 160-170 seconds, but I seem to have hit a dead end. Here's the overarching script (just in case): + Show Spoiler +clc clear
%%% PARAMETERS %%%
answer = input('Should the default value be used for the edge region parameter? Y/N [Y]: ','s'); if isempty(answer) % If the user simply presses 'Enter', the answer is assumed to be 'Y'. answer = 'Y'; end
disp(' ') % Empty line for clarity.
if strcmp(answer, 'Y') == 1 edge_parameter = 0.25; disp('The default value of the edge region parameter (0.25) will be used.') else edge_parameter = input('Please, specify the value of the edge region parameter: '); end
%%% IMPORT OF DATA %%%
filename = 'banding_coordinates.txt'; delimiterIn = ' '; headerlinesIn = 1; D = importdata(filename,delimiterIn,headerlinesIn);
%%% EXTRACTION OF COORDINATES %%%
x_extracted = D.data(:,1); y_extracted = D.data(:,2);
% Transposition into vectors:
x = transpose(x_extracted); y = transpose(y_extracted);
% Verification of vector lengths:
if length(x) ~= length(y) disp('The imported data is corrupted: the numbers of the x and y coordinates do not match') end
%%% BOUNDARY CONDITIONS %%%
% Limits of the analyzed space:
x_min = min(x); x_max = max(x); y_min = min(y); y_max = max(y);
% Edge region parameters for each coordinate:
edge_parameter_x = edge_parameter*(x_max - x_min); edge_parameter_y = edge_parameter*(y_max - y_min);
% Limits of the focal point space:
x_foc_min = x_min + edge_parameter_x; x_foc_max = x_max - edge_parameter_x; y_foc_min = y_min + edge_parameter_y; y_foc_max = y_max - edge_parameter_y;
%%% EXTRACTION OF FOCAL POINTS %%%
% Generation of focal points:
[in,on] = inpolygon(x, y, [x_foc_min, x_foc_max], [y_foc_min, y_foc_max] ;
x_foc = x(in); y_foc = y(in);
% Focal points plot:
answer = input('Would you like to display the generated focal points? Y/N [N]: ','s'); if isempty(answer) % If the user simply presses 'Enter', the answer is assumed to be 'Y'. answer = 'N'; end
disp(' ') % Empty line for clarity.
if strcmp(answer, 'Y') == 1 figure
scatter(x,y,'.') % all points axis equal
hold on scatter(x_foc,y_foc,'r+') % focal points hold off end
%%% ANALYSIS %%% tic overall_variance = calculate_variance_test6_3(x, y, x_foc, y_foc, ... x_min, x_max, y_min, y_max); toc %%% OUTPUT %%%
% Original variance peak:
disp('Original variance peak:') max(overall_variance)
% Plot:
figure
plot(1:1:180, overall_variance) And here's the function that does the actual analysis (in its most optimal form): + Show Spoiler +function variance_avg_foc = calculate_variance_test6_3(x, y, x_foc, y_foc, ... x_min, x_max, y_min, y_max)
% Calculation of the number of focal points:
m = length(x_foc); % Number of focal points = number of iterations in the first loop.
% Calculation of the number of scales:
scale = 1:1:45; % All scales. l = length(scale);
% Preallocation of memory for all focal points:
variance_avg_scale = zeros(m, 180); % Calculation of the number of analyzed points: n = length(x) - 1; % Calculation of transect properties:
p = 1:1:180; transect_angles = p - 1; slopes = tand((transect_angles)); % Selection of the focal point
for j = 1:m % Removal of the focal point from the pool of analyzed points: index = find( (x(1, == x_foc(j)) & (y(1, == y_foc(j)) ); if length(index) ~= 1 % If there is more than once point with the x_foc(j), y_foc(j) index = index(1); % coordinates, only one point will be subtracted from the future end % local coordinates pool. x_without_foc = x(1:end ~= index); y_without_foc = y(1:end ~= index); % Creation of local coordinates /vectorized:
i = 1:1:n; % One iteration for each (x, y) point. x_local = x_without_foc(i) - x_foc(j); y_local = y_without_foc(i) - y_foc(j);
% Limits of the studied area in local coordinates: x_min_local = x_min - x_foc(j); x_max_local = x_max - x_foc(j); y_min_local = y_min - y_foc(j); y_max_local = y_max - y_foc(j); % Translation to local polar coordinates:
[angle,radius] = cart_to_pol(x_local,y_local); angle = rad2deg(angle); % Consolidation of coordinates data:
data = [angle; radius];
% Preallocation of memory for all scales:
transect_area = zeros(1, 180); % Calculation of the transect area: for p = 1:180 transect_angle = transect_angles(p); slope = slopes(p); transect_area(p) = calculate_transect_area(x_min_local, ... x_max_local, y_min_local, y_max_local, transect_angle, slope); end % Preallocation of memory for all scales:
variance_scale = zeros(l, 180);
% Selection of scale:
for k = 1:l
% Transect width: transect_width = scale; % Preallocation of memory for all positions:
wavelet_transform = zeros(1, 180); variance_normalized = zeros(1, 180);
% Selection of angular position:
for p = 1:180 % Finding the points inside the specified transect:
points_observed = observe_points(data, p, transect_width(k));
% Extraction of the polar coordinates of the observed points:
angle_observed = points_observed(1, ; radius_observed = points_observed(2, ;
% Calculation of the number of observed points:
observations = length(angle_observed); % Calculation of scaled wavelets for each point within the transect /vectorized:
if observations == 0 scaled_wavelet = 0; else o = 1:1:observations; t = (angle_observed(o) - transect_angles(p))./scale(k); scaled_wavelet = radius_observed.*wavelet_function(t); end wavelet_transform(p) = sum(scaled_wavelet)/scale(k); variance_normalized(p) = wavelet_transform(p)^2/transect_area(p); end variance_scale(k, = variance_normalized; end variance_avg_scale(j, = sum(variance_scale)./l; % All specified scales. end variance_avg_foc = sum(variance_avg_scale)./m; end
%-------------------------------------------------------------------------- function [angle, radius] = cart_to_pol(x, y)
angle = arctan(y, x); radius = hypot(x,y);
end
%-------------------------------------------------------------------------- function angle = arctan(y, x) n = length(x); angle = zeros(1, n); parfor i = 1:n % /vectorized angle(i) = atan2(y(i), x(i)); if angle(i) < 0 angle(i) = angle(i) + pi; end end end
%-------------------------------------------------------------------------- function points_observed = observe_points(data, p, transect_width)
% Finding the points inside the specified transect:
step_size = 1;
log_ind = ( (p - 1 - (transect_width/2))*step_size <= data(1, & ... data(1, <= (p - 1 + (transect_width/2))*step_size ) | ... ( ((p - 1 - (transect_width/2))*step_size + 180) <= data(1, & ... data(1, <= ((p - 1 + (transect_width/2))*step_size + 180) ); points_observed = data(:,log_ind);
end
%-------------------------------------------------------------------------- function transect_area = calculate_transect_area(x_min_local, ... x_max_local, y_min_local, y_max_local, transect_angle, slope) % Slope variants: if (transect_angle == 0) || (transect_angle == 180) % Horizontal line. x_lim_1 = x_min_local; x_lim_2 = x_max_local; y_lim_1 = 0; y_lim_2 = 0; elseif transect_angle == 90 % Vertical line. x_lim_1 = 0; x_lim_2 = 0; y_lim_1 = y_min_local; y_lim_2 = y_max_local; elseif transect_angle > 0 && transect_angle < 90 % Boundary coordinates:
if y_min_local/slope >= x_min_local x_lim_1 = y_min_local/slope; y_lim_1 = y_min_local; else x_lim_1 = x_min_local; y_lim_1 = slope*x_min_local; end
if y_max_local/slope <= x_max_local x_lim_2 = y_max_local/slope; y_lim_2 = y_max_local; else x_lim_2 = x_max_local; y_lim_2 = slope*x_max_local; end else
% Boundary coordinates:
if y_min_local/slope >= x_max_local x_lim_1 = x_max_local; y_lim_1 = slope*x_max_local; else x_lim_1 = y_min_local/slope; y_lim_1 = y_min_local; end
if y_max_local/slope <= x_min_local x_lim_2 = x_min_local; y_lim_2 = slope*x_min_local; else x_lim_2 = y_max_local/slope; y_lim_2 = y_max_local; end end transect_area = (x_lim_1)^2 + (y_lim_1)^2 + (x_lim_2)^2 + (y_lim_2)^2; % This is not real area of the transect. % It's area divided by omega/2, but that holds % for all transects and omega is constant. end
%-------------------------------------------------------------------------- function mexican_hat = wavelet_function(t) % /vectorized
mexican_hat = 2./(sqrt(3)).*pi.^(-1/4).*(1 - 4.*t.^2).*exp(-2.*t.^2); % Wavelet function normalized to have unit energy. % Vectorizing with "vectorize()" does not work. % It cannot tell if I want to use ^ for matrices (mpower) % or scalars (power). end According to the MATLAB profiler, there are two lines/functions whose vectorization could yield significant improvements in calculation time, but I failed at vectorizing them properly... Any help would be appreciated. This line currently requires roughly 95 seconds (it's being called 4 million times): + Show Spoiler +points_observed = observe_points(data, p, transect_width(k)); And these two lines require roughly 30 seconds (they're being called 700 thousand times): + Show Spoiler +t = (angle_observed(o) - transect_angles(p))./scale(k); scaled_wavelet = radius_observed.*wavelet_function(t); Additionally, I am looking for a way to make scaled_wavelet call wavelet_function only for non-zero values of radius_observed. I tried the solution below, but it only made the calculations longer: + Show Spoiler +% Calculation of scaled wavelets for each point within the transect /vectorized:
if observations == 0 scaled_wavelet = 0; else log_ind = (radius_observed ~= 0); %o = 1:1:observations; scaled_wavelet = radius_observed(log_ind).*wavelet_function((angle_observed(log_ind) - transect_angles(p))./scale(k)); end I also found that spfun() does something similar, but I don't know how to make it work with a custom function like wavelet_function(). So could anyone lend me a helping hand with this? Well, the biggest issue is the number of times the line is hit. Any way to reduce it? Unroll the loops? Maybe store the needed info in a better structure to later calculate points? Make more functions (and name them aptly) to improve readability.
Next, and I don't know if it helps, but you can use more local variables in the observe_points function. I think you can optimize the calculation in observe_points? Perhaps something like
function points_observed = observe_points(data, p, transect_width) % Finding the points inside the specified transect: t = transect_width/2; plow = p - 1 - t; phigh = p - 1 + t log_ind = (plow <= data(1,:) & data(1,:) <= phigh) | ((plow + 180) <= data(1,:) & data(1,:) <= (phigh + 180)); points_observed = data(:,log_ind); end
I'm not familiar with mathlab at all though...
|
On October 01 2016 00:40 maybenexttime wrote:Show nested quote +On September 27 2016 22:57 maybenexttime wrote:Let me first introduce myself briefly. I am a mechanical engineering graduate who'd like to learn some Python and MATLAB programming, because I consider programming an important skill in every engineer's repertoire. For the past two months I have studied both languages and worked on some private projects of mine since that seems to be the best way to learn useful things. I'm currently working on a MATLAB program that is supposed to quantify the degree of banding in the microstructures I simulated with an earlier script (the output is basically a set of (x,y) coordinates), and I've encountered some performance issues when compared with the program whose analysis I am trying to recreate (PASSaGE, link below). My understanding of how the analysis is supposed to work is based on the paper written by the software's author (link) and the program's manual, but some sections were quite misleading in my opinion, so I had to do several modifications based on some experiments with analysis. http://www.passagesoftware.net/http://rosenberglab.net/Pubs/JVegSci2004v15p277.pdfInitially, my MATLAB script took roughly 15 minutes (more specifically, the function that runs the analysis, not the overarching script) to analyze a set of 1000 points, whereas PASSaGE takes around 1-2 minutes. It relied mostly on loops. After a few days of optimizing, I managed to decrease the calculation time to 160-170 seconds, but I seem to have hit a dead end. Here's the overarching script (just in case): + Show Spoiler +clc clear
%%% PARAMETERS %%%
answer = input('Should the default value be used for the edge region parameter? Y/N [Y]: ','s'); if isempty(answer) % If the user simply presses 'Enter', the answer is assumed to be 'Y'. answer = 'Y'; end
disp(' ') % Empty line for clarity.
if strcmp(answer, 'Y') == 1 edge_parameter = 0.25; disp('The default value of the edge region parameter (0.25) will be used.') else edge_parameter = input('Please, specify the value of the edge region parameter: '); end
%%% IMPORT OF DATA %%%
filename = 'banding_coordinates.txt'; delimiterIn = ' '; headerlinesIn = 1; D = importdata(filename,delimiterIn,headerlinesIn);
%%% EXTRACTION OF COORDINATES %%%
x_extracted = D.data(:,1); y_extracted = D.data(:,2);
% Transposition into vectors:
x = transpose(x_extracted); y = transpose(y_extracted);
% Verification of vector lengths:
if length(x) ~= length(y) disp('The imported data is corrupted: the numbers of the x and y coordinates do not match') end
%%% BOUNDARY CONDITIONS %%%
% Limits of the analyzed space:
x_min = min(x); x_max = max(x); y_min = min(y); y_max = max(y);
% Edge region parameters for each coordinate:
edge_parameter_x = edge_parameter*(x_max - x_min); edge_parameter_y = edge_parameter*(y_max - y_min);
% Limits of the focal point space:
x_foc_min = x_min + edge_parameter_x; x_foc_max = x_max - edge_parameter_x; y_foc_min = y_min + edge_parameter_y; y_foc_max = y_max - edge_parameter_y;
%%% EXTRACTION OF FOCAL POINTS %%%
% Generation of focal points:
[in,on] = inpolygon(x, y, [x_foc_min, x_foc_max], [y_foc_min, y_foc_max] ;
x_foc = x(in); y_foc = y(in);
% Focal points plot:
answer = input('Would you like to display the generated focal points? Y/N [N]: ','s'); if isempty(answer) % If the user simply presses 'Enter', the answer is assumed to be 'Y'. answer = 'N'; end
disp(' ') % Empty line for clarity.
if strcmp(answer, 'Y') == 1 figure
scatter(x,y,'.') % all points axis equal
hold on scatter(x_foc,y_foc,'r+') % focal points hold off end
%%% ANALYSIS %%% tic overall_variance = calculate_variance_test6_3(x, y, x_foc, y_foc, ... x_min, x_max, y_min, y_max); toc %%% OUTPUT %%%
% Original variance peak:
disp('Original variance peak:') max(overall_variance)
% Plot:
figure
plot(1:1:180, overall_variance) And here's the function that does the actual analysis (in its most optimal form): + Show Spoiler +function variance_avg_foc = calculate_variance_test6_3(x, y, x_foc, y_foc, ... x_min, x_max, y_min, y_max)
% Calculation of the number of focal points:
m = length(x_foc); % Number of focal points = number of iterations in the first loop.
% Calculation of the number of scales:
scale = 1:1:45; % All scales. l = length(scale);
% Preallocation of memory for all focal points:
variance_avg_scale = zeros(m, 180); % Calculation of the number of analyzed points: n = length(x) - 1; % Calculation of transect properties:
p = 1:1:180; transect_angles = p - 1; slopes = tand((transect_angles)); % Selection of the focal point
for j = 1:m % Removal of the focal point from the pool of analyzed points: index = find( (x(1, == x_foc(j)) & (y(1, == y_foc(j)) ); if length(index) ~= 1 % If there is more than once point with the x_foc(j), y_foc(j) index = index(1); % coordinates, only one point will be subtracted from the future end % local coordinates pool. x_without_foc = x(1:end ~= index); y_without_foc = y(1:end ~= index); % Creation of local coordinates /vectorized:
i = 1:1:n; % One iteration for each (x, y) point. x_local = x_without_foc(i) - x_foc(j); y_local = y_without_foc(i) - y_foc(j);
% Limits of the studied area in local coordinates: x_min_local = x_min - x_foc(j); x_max_local = x_max - x_foc(j); y_min_local = y_min - y_foc(j); y_max_local = y_max - y_foc(j); % Translation to local polar coordinates:
[angle,radius] = cart_to_pol(x_local,y_local); angle = rad2deg(angle); % Consolidation of coordinates data:
data = [angle; radius];
% Preallocation of memory for all scales:
transect_area = zeros(1, 180); % Calculation of the transect area: for p = 1:180 transect_angle = transect_angles(p); slope = slopes(p); transect_area(p) = calculate_transect_area(x_min_local, ... x_max_local, y_min_local, y_max_local, transect_angle, slope); end % Preallocation of memory for all scales:
variance_scale = zeros(l, 180);
% Selection of scale:
for k = 1:l
% Transect width: transect_width = scale; % Preallocation of memory for all positions:
wavelet_transform = zeros(1, 180); variance_normalized = zeros(1, 180);
% Selection of angular position:
for p = 1:180 % Finding the points inside the specified transect:
points_observed = observe_points(data, p, transect_width(k));
% Extraction of the polar coordinates of the observed points:
angle_observed = points_observed(1, ; radius_observed = points_observed(2, ;
% Calculation of the number of observed points:
observations = length(angle_observed); % Calculation of scaled wavelets for each point within the transect /vectorized:
if observations == 0 scaled_wavelet = 0; else o = 1:1:observations; t = (angle_observed(o) - transect_angles(p))./scale(k); scaled_wavelet = radius_observed.*wavelet_function(t); end wavelet_transform(p) = sum(scaled_wavelet)/scale(k); variance_normalized(p) = wavelet_transform(p)^2/transect_area(p); end variance_scale(k, = variance_normalized; end variance_avg_scale(j, = sum(variance_scale)./l; % All specified scales. end variance_avg_foc = sum(variance_avg_scale)./m; end
%-------------------------------------------------------------------------- function [angle, radius] = cart_to_pol(x, y)
angle = arctan(y, x); radius = hypot(x,y);
end
%-------------------------------------------------------------------------- function angle = arctan(y, x) n = length(x); angle = zeros(1, n); parfor i = 1:n % /vectorized angle(i) = atan2(y(i), x(i)); if angle(i) < 0 angle(i) = angle(i) + pi; end end end
%-------------------------------------------------------------------------- function points_observed = observe_points(data, p, transect_width)
% Finding the points inside the specified transect:
step_size = 1;
log_ind = ( (p - 1 - (transect_width/2))*step_size <= data(1, & ... data(1, <= (p - 1 + (transect_width/2))*step_size ) | ... ( ((p - 1 - (transect_width/2))*step_size + 180) <= data(1, & ... data(1, <= ((p - 1 + (transect_width/2))*step_size + 180) ); points_observed = data(:,log_ind);
end
%-------------------------------------------------------------------------- function transect_area = calculate_transect_area(x_min_local, ... x_max_local, y_min_local, y_max_local, transect_angle, slope) % Slope variants: if (transect_angle == 0) || (transect_angle == 180) % Horizontal line. x_lim_1 = x_min_local; x_lim_2 = x_max_local; y_lim_1 = 0; y_lim_2 = 0; elseif transect_angle == 90 % Vertical line. x_lim_1 = 0; x_lim_2 = 0; y_lim_1 = y_min_local; y_lim_2 = y_max_local; elseif transect_angle > 0 && transect_angle < 90 % Boundary coordinates:
if y_min_local/slope >= x_min_local x_lim_1 = y_min_local/slope; y_lim_1 = y_min_local; else x_lim_1 = x_min_local; y_lim_1 = slope*x_min_local; end
if y_max_local/slope <= x_max_local x_lim_2 = y_max_local/slope; y_lim_2 = y_max_local; else x_lim_2 = x_max_local; y_lim_2 = slope*x_max_local; end else
% Boundary coordinates:
if y_min_local/slope >= x_max_local x_lim_1 = x_max_local; y_lim_1 = slope*x_max_local; else x_lim_1 = y_min_local/slope; y_lim_1 = y_min_local; end
if y_max_local/slope <= x_min_local x_lim_2 = x_min_local; y_lim_2 = slope*x_min_local; else x_lim_2 = y_max_local/slope; y_lim_2 = y_max_local; end end transect_area = (x_lim_1)^2 + (y_lim_1)^2 + (x_lim_2)^2 + (y_lim_2)^2; % This is not real area of the transect. % It's area divided by omega/2, but that holds % for all transects and omega is constant. end
%-------------------------------------------------------------------------- function mexican_hat = wavelet_function(t) % /vectorized
mexican_hat = 2./(sqrt(3)).*pi.^(-1/4).*(1 - 4.*t.^2).*exp(-2.*t.^2); % Wavelet function normalized to have unit energy. % Vectorizing with "vectorize()" does not work. % It cannot tell if I want to use ^ for matrices (mpower) % or scalars (power). end According to the MATLAB profiler, there are two lines/functions whose vectorization could yield significant improvements in calculation time, but I failed at vectorizing them properly... Any help would be appreciated. This line currently requires roughly 95 seconds (it's being called 4 million times): + Show Spoiler +points_observed = observe_points(data, p, transect_width(k)); And these two lines require roughly 30 seconds (they're being called 700 thousand times): + Show Spoiler +t = (angle_observed(o) - transect_angles(p))./scale(k); scaled_wavelet = radius_observed.*wavelet_function(t); Additionally, I am looking for a way to make scaled_wavelet call wavelet_function only for non-zero values of radius_observed. I tried the solution below, but it only made the calculations longer: + Show Spoiler +% Calculation of scaled wavelets for each point within the transect /vectorized:
if observations == 0 scaled_wavelet = 0; else log_ind = (radius_observed ~= 0); %o = 1:1:observations; scaled_wavelet = radius_observed(log_ind).*wavelet_function((angle_observed(log_ind) - transect_angles(p))./scale(k)); end I also found that spfun() does something similar, but I don't know how to make it work with a custom function like wavelet_function(). So could anyone lend me a helping hand with this?
Do you know linear algebra? Looking at your implementation, it seems that you are doing everything using loops instead of vectorized forms. Vectorized forms in matlab are far more efficient than using loops. I saw a comment that said "vectorized", but I am not sure you know what this means.
Writing some test code to help you understand the difference: If you start making optimizations like this in your programs, your program to should start running a lot faster.
Vecotorized versions: (ex) A = ones(5,5) B = eye(1,5) C = A + B
Loop version: (ex)
A = ones(5,5) B = eye(1,5) for r = 1 : 5 for c = 1 : 5 if (r == c) c(a,b) = a(r,c) + b(r,c)
Additionally, use markers to pinpoint the code that is taking the longest to run. A measure of the count of a function call alone is not a good indicator of bad performance. Try to figure out using print statements. Log the time taken in your program at each step and keep breaking it down and figuring how what is really taking time.
|
On October 01 2016 00:40 maybenexttime wrote:Show nested quote +On September 27 2016 22:57 maybenexttime wrote:Let me first introduce myself briefly. I am a mechanical engineering graduate who'd like to learn some Python and MATLAB programming, because I consider programming an important skill in every engineer's repertoire. For the past two months I have studied both languages and worked on some private projects of mine since that seems to be the best way to learn useful things. I'm currently working on a MATLAB program that is supposed to quantify the degree of banding in the microstructures I simulated with an earlier script (the output is basically a set of (x,y) coordinates), and I've encountered some performance issues when compared with the program whose analysis I am trying to recreate (PASSaGE, link below). My understanding of how the analysis is supposed to work is based on the paper written by the software's author (link) and the program's manual, but some sections were quite misleading in my opinion, so I had to do several modifications based on some experiments with analysis. http://www.passagesoftware.net/http://rosenberglab.net/Pubs/JVegSci2004v15p277.pdfInitially, my MATLAB script took roughly 15 minutes (more specifically, the function that runs the analysis, not the overarching script) to analyze a set of 1000 points, whereas PASSaGE takes around 1-2 minutes. It relied mostly on loops. After a few days of optimizing, I managed to decrease the calculation time to 160-170 seconds, but I seem to have hit a dead end. Here's the overarching script (just in case): + Show Spoiler +clc clear
%%% PARAMETERS %%%
answer = input('Should the default value be used for the edge region parameter? Y/N [Y]: ','s'); if isempty(answer) % If the user simply presses 'Enter', the answer is assumed to be 'Y'. answer = 'Y'; end
disp(' ') % Empty line for clarity.
if strcmp(answer, 'Y') == 1 edge_parameter = 0.25; disp('The default value of the edge region parameter (0.25) will be used.') else edge_parameter = input('Please, specify the value of the edge region parameter: '); end
%%% IMPORT OF DATA %%%
filename = 'banding_coordinates.txt'; delimiterIn = ' '; headerlinesIn = 1; D = importdata(filename,delimiterIn,headerlinesIn);
%%% EXTRACTION OF COORDINATES %%%
x_extracted = D.data(:,1); y_extracted = D.data(:,2);
% Transposition into vectors:
x = transpose(x_extracted); y = transpose(y_extracted);
% Verification of vector lengths:
if length(x) ~= length(y) disp('The imported data is corrupted: the numbers of the x and y coordinates do not match') end
%%% BOUNDARY CONDITIONS %%%
% Limits of the analyzed space:
x_min = min(x); x_max = max(x); y_min = min(y); y_max = max(y);
% Edge region parameters for each coordinate:
edge_parameter_x = edge_parameter*(x_max - x_min); edge_parameter_y = edge_parameter*(y_max - y_min);
% Limits of the focal point space:
x_foc_min = x_min + edge_parameter_x; x_foc_max = x_max - edge_parameter_x; y_foc_min = y_min + edge_parameter_y; y_foc_max = y_max - edge_parameter_y;
%%% EXTRACTION OF FOCAL POINTS %%%
% Generation of focal points:
[in,on] = inpolygon(x, y, [x_foc_min, x_foc_max], [y_foc_min, y_foc_max] ;
x_foc = x(in); y_foc = y(in);
% Focal points plot:
answer = input('Would you like to display the generated focal points? Y/N [N]: ','s'); if isempty(answer) % If the user simply presses 'Enter', the answer is assumed to be 'Y'. answer = 'N'; end
disp(' ') % Empty line for clarity.
if strcmp(answer, 'Y') == 1 figure
scatter(x,y,'.') % all points axis equal
hold on scatter(x_foc,y_foc,'r+') % focal points hold off end
%%% ANALYSIS %%% tic overall_variance = calculate_variance_test6_3(x, y, x_foc, y_foc, ... x_min, x_max, y_min, y_max); toc %%% OUTPUT %%%
% Original variance peak:
disp('Original variance peak:') max(overall_variance)
% Plot:
figure
plot(1:1:180, overall_variance) And here's the function that does the actual analysis (in its most optimal form): + Show Spoiler +function variance_avg_foc = calculate_variance_test6_3(x, y, x_foc, y_foc, ... x_min, x_max, y_min, y_max)
% Calculation of the number of focal points:
m = length(x_foc); % Number of focal points = number of iterations in the first loop.
% Calculation of the number of scales:
scale = 1:1:45; % All scales. l = length(scale);
% Preallocation of memory for all focal points:
variance_avg_scale = zeros(m, 180); % Calculation of the number of analyzed points: n = length(x) - 1; % Calculation of transect properties:
p = 1:1:180; transect_angles = p - 1; slopes = tand((transect_angles)); % Selection of the focal point
for j = 1:m % Removal of the focal point from the pool of analyzed points: index = find( (x(1, == x_foc(j)) & (y(1, == y_foc(j)) ); if length(index) ~= 1 % If there is more than once point with the x_foc(j), y_foc(j) index = index(1); % coordinates, only one point will be subtracted from the future end % local coordinates pool. x_without_foc = x(1:end ~= index); y_without_foc = y(1:end ~= index); % Creation of local coordinates /vectorized:
i = 1:1:n; % One iteration for each (x, y) point. x_local = x_without_foc(i) - x_foc(j); y_local = y_without_foc(i) - y_foc(j);
% Limits of the studied area in local coordinates: x_min_local = x_min - x_foc(j); x_max_local = x_max - x_foc(j); y_min_local = y_min - y_foc(j); y_max_local = y_max - y_foc(j); % Translation to local polar coordinates:
[angle,radius] = cart_to_pol(x_local,y_local); angle = rad2deg(angle); % Consolidation of coordinates data:
data = [angle; radius];
% Preallocation of memory for all scales:
transect_area = zeros(1, 180); % Calculation of the transect area: for p = 1:180 transect_angle = transect_angles(p); slope = slopes(p); transect_area(p) = calculate_transect_area(x_min_local, ... x_max_local, y_min_local, y_max_local, transect_angle, slope); end % Preallocation of memory for all scales:
variance_scale = zeros(l, 180);
% Selection of scale:
for k = 1:l
% Transect width: transect_width = scale; % Preallocation of memory for all positions:
wavelet_transform = zeros(1, 180); variance_normalized = zeros(1, 180);
% Selection of angular position:
for p = 1:180 % Finding the points inside the specified transect:
points_observed = observe_points(data, p, transect_width(k));
% Extraction of the polar coordinates of the observed points:
angle_observed = points_observed(1, ; radius_observed = points_observed(2, ;
% Calculation of the number of observed points:
observations = length(angle_observed); % Calculation of scaled wavelets for each point within the transect /vectorized:
if observations == 0 scaled_wavelet = 0; else o = 1:1:observations; t = (angle_observed(o) - transect_angles(p))./scale(k); scaled_wavelet = radius_observed.*wavelet_function(t); end wavelet_transform(p) = sum(scaled_wavelet)/scale(k); variance_normalized(p) = wavelet_transform(p)^2/transect_area(p); end variance_scale(k, = variance_normalized; end variance_avg_scale(j, = sum(variance_scale)./l; % All specified scales. end variance_avg_foc = sum(variance_avg_scale)./m; end
%-------------------------------------------------------------------------- function [angle, radius] = cart_to_pol(x, y)
angle = arctan(y, x); radius = hypot(x,y);
end
%-------------------------------------------------------------------------- function angle = arctan(y, x) n = length(x); angle = zeros(1, n); parfor i = 1:n % /vectorized angle(i) = atan2(y(i), x(i)); if angle(i) < 0 angle(i) = angle(i) + pi; end end end
%-------------------------------------------------------------------------- function points_observed = observe_points(data, p, transect_width)
% Finding the points inside the specified transect:
step_size = 1;
log_ind = ( (p - 1 - (transect_width/2))*step_size <= data(1, & ... data(1, <= (p - 1 + (transect_width/2))*step_size ) | ... ( ((p - 1 - (transect_width/2))*step_size + 180) <= data(1, & ... data(1, <= ((p - 1 + (transect_width/2))*step_size + 180) ); points_observed = data(:,log_ind);
end
%-------------------------------------------------------------------------- function transect_area = calculate_transect_area(x_min_local, ... x_max_local, y_min_local, y_max_local, transect_angle, slope) % Slope variants: if (transect_angle == 0) || (transect_angle == 180) % Horizontal line. x_lim_1 = x_min_local; x_lim_2 = x_max_local; y_lim_1 = 0; y_lim_2 = 0; elseif transect_angle == 90 % Vertical line. x_lim_1 = 0; x_lim_2 = 0; y_lim_1 = y_min_local; y_lim_2 = y_max_local; elseif transect_angle > 0 && transect_angle < 90 % Boundary coordinates:
if y_min_local/slope >= x_min_local x_lim_1 = y_min_local/slope; y_lim_1 = y_min_local; else x_lim_1 = x_min_local; y_lim_1 = slope*x_min_local; end
if y_max_local/slope <= x_max_local x_lim_2 = y_max_local/slope; y_lim_2 = y_max_local; else x_lim_2 = x_max_local; y_lim_2 = slope*x_max_local; end else
% Boundary coordinates:
if y_min_local/slope >= x_max_local x_lim_1 = x_max_local; y_lim_1 = slope*x_max_local; else x_lim_1 = y_min_local/slope; y_lim_1 = y_min_local; end
if y_max_local/slope <= x_min_local x_lim_2 = x_min_local; y_lim_2 = slope*x_min_local; else x_lim_2 = y_max_local/slope; y_lim_2 = y_max_local; end end transect_area = (x_lim_1)^2 + (y_lim_1)^2 + (x_lim_2)^2 + (y_lim_2)^2; % This is not real area of the transect. % It's area divided by omega/2, but that holds % for all transects and omega is constant. end
%-------------------------------------------------------------------------- function mexican_hat = wavelet_function(t) % /vectorized
mexican_hat = 2./(sqrt(3)).*pi.^(-1/4).*(1 - 4.*t.^2).*exp(-2.*t.^2); % Wavelet function normalized to have unit energy. % Vectorizing with "vectorize()" does not work. % It cannot tell if I want to use ^ for matrices (mpower) % or scalars (power). end According to the MATLAB profiler, there are two lines/functions whose vectorization could yield significant improvements in calculation time, but I failed at vectorizing them properly... Any help would be appreciated. This line currently requires roughly 95 seconds (it's being called 4 million times): + Show Spoiler +points_observed = observe_points(data, p, transect_width(k)); And these two lines require roughly 30 seconds (they're being called 700 thousand times): + Show Spoiler +t = (angle_observed(o) - transect_angles(p))./scale(k); scaled_wavelet = radius_observed.*wavelet_function(t); Additionally, I am looking for a way to make scaled_wavelet call wavelet_function only for non-zero values of radius_observed. I tried the solution below, but it only made the calculations longer: + Show Spoiler +% Calculation of scaled wavelets for each point within the transect /vectorized:
if observations == 0 scaled_wavelet = 0; else log_ind = (radius_observed ~= 0); %o = 1:1:observations; scaled_wavelet = radius_observed(log_ind).*wavelet_function((angle_observed(log_ind) - transect_angles(p))./scale(k)); end I also found that spfun() does something similar, but I don't know how to make it work with a custom function like wavelet_function(). So could anyone lend me a helping hand with this?
I can take a look later tonight, did a bunch of Matlab for a couple months, Matlab optimization is hella whack, involves lots of parallelism instead.
|
On October 01 2016 01:44 supereddie wrote:Well, the biggest issue is the number of times the line is hit. Any way to reduce it? Unroll the loops? Maybe store the needed info in a better structure to later calculate points? Make more functions (and name them aptly) to improve readability.
Next, and I don't know if it helps, but you can use more local variables in the observe_points function. I think you can optimize the calculation in observe_points? Perhaps something like
I'm not familiar with mathlab at all though...
That's the essence of my question: how to vectorize (and thus remove the loop altogether) the following line in order to reduce the number of times it is called 180 times.
points_observed = observe_points(data, p, transect_width(k));
I tried your suggestion about optimizing calculations in observe_points, but it made no difference. As for making more functions, I might do it in the final version of the code. Vectorization might require some major overhauls, so more functions would make implementing changes harder.
On October 01 2016 02:56 Neshapotamus wrote:Do you know linear algebra? Looking at your implementation, it seems that you are doing everything using loops instead of vectorized forms. Vectorized forms in matlab are far more efficient than using loops. I saw a comment that said "vectorized", but I am not sure you know what this means.
Writing some test code to help you understand the difference: If you start making optimizations like this in your programs, your program to should start running a lot faster.
Vecotorized versions: (ex) A = ones(5,5) B = eye(1,5) C = A + B
Loop version: (ex)
A = ones(5,5) B = eye(1,5) for r = 1 : 5 for c = 1 : 5 if (r == c) c(a,b) = a(r,c) + b(r,c)
Additionally, use markers to pinpoint the code that is taking the longest to run. A measure of the count of a function call alone is not a good indicator of bad performance. Try to figure out using print statements. Log the time taken in your program at each step and keep breaking it down and figuring how what is really taking time.
I know what vectorization means. I didn't when I first started working on this code (although I used it unknowingly a couple of times), since I only had one semester of MATLAB fundamentals and we mostly relied on loops. In the course of working on this problem, I learned what vectorization is and implemented it wherever I could.
The problem is that it's rather trivial when you're dealing with one loop, but it gets very tricky when you're dealing with multiple loops and have to aggregate data and call different functions at various stages. Then you're dealing with multidimensional matrices and have to make sure you meet various criteria MATLAB functions have, such as matching matrix dimensions and so on.
E.g. here's an example of how I tried to vectorize the p forloop:
+ Show Spoiler + % Selection of angular position:
p = 1:1:180;
% Transect angle: transect_angle = p - 1; % Finding the points inside the specified transect: [angles_observed, radii_observed] = observe_points(angle, radius, p, transect_width(k)); for p = 1:180 % Calculation of the transect area:
transect_area = zeros(1, 180); transect_area(p) = calculate_transect_area(x_min_local, ... x_max_local, y_min_local, y_max_local, transect_angle(p));
% Preallocation of memory:
scaled_wavelet = zeros(1, n); t = zeros(p, n);
% Calculation of scaled wavelets for each point within the transect /vectorized: o = 1:1:n; t = (angles_observed(p, o) - transect_angle(p))./scale(k); scaled_wavelet = radii_observed(p,:).*wavelet_function(t); wavelet_transform(p) = sum(scaled_wavelet)/scale(k); % No density measure. variance_normalized(p) = wavelet_transform(p)^2/transect_area(p); end variance_scale(k,:) = variance_normalized;
(...)
function [angles_observed, radii_observed] = observe_points(angle, radius, p, transect_width)
% Finding the points inside the specified transect:
i = 1:1:180; step_size = 1; angle_array = repmat((angle(1,:)), length(p), 1); radius_array = repmat((radius(1,:)), length(p), 1); p = transpose(p); transect_lim_min = repmat(((p(i,1) - 1 - (transect_width/2)).*step_size),1,length(angle)); transect_lim_max = repmat(((p(i,1) - 1 + (transect_width/2)).*step_size),1,length(angle));
log_ind = ( transect_lim_min <= angle_array & ... angle_array <= transect_lim_max ) | ... ( (transect_lim_min + 180) <= angle_array & ... angle_array <= (transect_lim_max + 180) );
angles_observed = angle_array.*log_ind; radii_observed = radius_array.*log_ind; end
While the loop was successfully vectorized to a large extent, the calculation time of observe_points increased by an order of magnitude because of those extra lines that were required to format the data. Perhaps my method of vectorization is not optimal, but I couldn't find a way to use logical indexing for multidimensional data without resorting to comparing corresponding elements of the matrices. That is why I asked for help.
I did not vectorize transect_area because I did not know how to properly vectorize all the if conditions in calculate_transect_area.
I did not vectorize "scaled_wavelet = radius_observed.*wavelet_function(t);" because I have absolutely no idea how to apply a custom function like wavelet_function to all non-zero elements in a radius_observed vector/matrix.
Is there any algorithm/guide to loop flattening/vectorization?
|
On September 30 2016 16:09 Aerisky wrote: I started with JS: The Good Parts. I'm not an expert though, so I don't know what's a good place to start.
Definitely don't start with AngularJS though IMO, I feel like it's pretty confusing in the beginning, especially if you haven't learned JavaScript yet.
thanks for the tips!
|
28078 Posts
On October 01 2016 09:49 necrosexy wrote:Show nested quote +On September 30 2016 16:09 Aerisky wrote: I started with JS: The Good Parts. I'm not an expert though, so I don't know what's a good place to start.
Definitely don't start with AngularJS though IMO, I feel like it's pretty confusing in the beginning, especially if you haven't learned JavaScript yet. thanks for the tips! I started with Eloquent Javascript (there's a free PDF on his website which is nice) and I think it was one of the best introduction books I've read. The Definitive Guide is great as well, although the purpose of that book is quite a bit different.
You need to figure out what you want out of the resource before picking one I guess. Most people recommend reading something like Eloquent Javascript first since it's more of a true intro to coding as well as an intro to JS. The Definitive Guide is more of an official JS documentation/reference book. It's like 1000+ pages and has literally everything about JS in it, although it doesn't necessarily take you on a step-by-step process like other books might. If you already know multiple languages and just want to learn the syntax and how JS works overall, then The Definitive Guide is probably better.
After you are done reading one of those books (or even while you are) I'd recommend taking advantage of something like codeacademy to further ingrain the basics
|
On September 30 2016 18:33 maybenexttime wrote:Show nested quote +On September 30 2016 14:46 Aerisky wrote: Yeah, normally I've applied to a ton of places for internships (100% agree that it's a numbers game). Now that I'm a soon-to-be new grad, though, I have a return offer from the company at which I interned this summer; I really like the offer so I've only been applying to really cool places where I'd be interested in working/prestigious places. Unless you are 100% sure that the company at which you interned will hire you, don't change your application strategy. I did this mistake. I was an intern in one company for a year, even did the research for my Master's thesis there, and I got accepted for the next wave of their graduate program after I would graduate (I still had one year of university after finishing the internship) or alternatively I could start a PhD at one of the universities they had an agreement with. Then when I was back at the university - a month or so before graduation last Summer - I was told by my previous supervisor that the new management decided to discontinue both programs... Never put all your eggs in one basket. If you get accepted by more companies, you can always decline. Wow, that sucks. That's actually really fucked up on that company's part, damn.
And yeah, I'll be applying to more places as well I suppose
|
On October 01 2016 07:23 maybenexttime wrote:Show nested quote +On October 01 2016 01:44 supereddie wrote:Well, the biggest issue is the number of times the line is hit. Any way to reduce it? Unroll the loops? Maybe store the needed info in a better structure to later calculate points? Make more functions (and name them aptly) to improve readability.
Next, and I don't know if it helps, but you can use more local variables in the observe_points function. I think you can optimize the calculation in observe_points? Perhaps something like
I'm not familiar with mathlab at all though... That's the essence of my question: how to vectorize (and thus remove the loop altogether) the following line in order to reduce the number of times it is called 180 times. points_observed = observe_points(data, p, transect_width(k)); I tried your suggestion about optimizing calculations in observe_points, but it made no difference. As for making more functions, I might do it in the final version of the code. Vectorization might require some major overhauls, so more functions would make implementing changes harder. Show nested quote +On October 01 2016 02:56 Neshapotamus wrote:Do you know linear algebra? Looking at your implementation, it seems that you are doing everything using loops instead of vectorized forms. Vectorized forms in matlab are far more efficient than using loops. I saw a comment that said "vectorized", but I am not sure you know what this means.
Writing some test code to help you understand the difference: If you start making optimizations like this in your programs, your program to should start running a lot faster.
Vecotorized versions: (ex) A = ones(5,5) B = eye(1,5) C = A + B
Loop version: (ex)
A = ones(5,5) B = eye(1,5) for r = 1 : 5 for c = 1 : 5 if (r == c) c(a,b) = a(r,c) + b(r,c)
Additionally, use markers to pinpoint the code that is taking the longest to run. A measure of the count of a function call alone is not a good indicator of bad performance. Try to figure out using print statements. Log the time taken in your program at each step and keep breaking it down and figuring how what is really taking time.
I know what vectorization means. I didn't when I first started working on this code (although I used it unknowingly a couple of times), since I only had one semester of MATLAB fundamentals and we mostly relied on loops. In the course of working on this problem, I learned what vectorization is and implemented it wherever I could. The problem is that it's rather trivial when you're dealing with one loop, but it gets very tricky when you're dealing with multiple loops and have to aggregate data and call different functions at various stages. Then you're dealing with multidimensional matrices and have to make sure you meet various criteria MATLAB functions have, such as matching matrix dimensions and so on. E.g. here's an example of how I tried to vectorize the p forloop: + Show Spoiler + % Selection of angular position:
p = 1:1:180;
% Transect angle: transect_angle = p - 1; % Finding the points inside the specified transect: [angles_observed, radii_observed] = observe_points(angle, radius, p, transect_width(k)); for p = 1:180 % Calculation of the transect area:
transect_area = zeros(1, 180); transect_area(p) = calculate_transect_area(x_min_local, ... x_max_local, y_min_local, y_max_local, transect_angle(p));
% Preallocation of memory:
scaled_wavelet = zeros(1, n); t = zeros(p, n);
% Calculation of scaled wavelets for each point within the transect /vectorized: o = 1:1:n; t = (angles_observed(p, o) - transect_angle(p))./scale(k); scaled_wavelet = radii_observed(p, .*wavelet_function(t); wavelet_transform(p) = sum(scaled_wavelet)/scale(k); % No density measure. variance_normalized(p) = wavelet_transform(p)^2/transect_area(p); end variance_scale(k, = variance_normalized;
(...)
function [angles_observed, radii_observed] = observe_points(angle, radius, p, transect_width)
% Finding the points inside the specified transect:
i = 1:1:180; step_size = 1; angle_array = repmat((angle(1, ), length(p), 1); radius_array = repmat((radius(1, ), length(p), 1); p = transpose(p); transect_lim_min = repmat(((p(i,1) - 1 - (transect_width/2)).*step_size),1,length(angle)); transect_lim_max = repmat(((p(i,1) - 1 + (transect_width/2)).*step_size),1,length(angle));
log_ind = ( transect_lim_min <= angle_array & ... angle_array <= transect_lim_max ) | ... ( (transect_lim_min + 180) <= angle_array & ... angle_array <= (transect_lim_max + 180) );
angles_observed = angle_array.*log_ind; radii_observed = radius_array.*log_ind; end While the loop was successfully vectorized to a large extent, the calculation time of observe_points increased by an order of magnitude because of those extra lines that were required to format the data. Perhaps my method of vectorization is not optimal, but I couldn't find a way to use logical indexing for multidimensional data without resorting to comparing corresponding elements of the matrices. That is why I asked for help. I did not vectorize transect_area because I did not know how to properly vectorize all the if conditions in calculate_transect_area. I did not vectorize "scaled_wavelet = radius_observed.*wavelet_function(t);" because I have absolutely no idea how to apply a custom function like wavelet_function to all non-zero elements in a radius_observed vector/matrix. Is there any algorithm/guide to loop flattening/vectorization?
There isn't any guide that I have seen.
Just to give you some insight into the matrix multiplication algorithm.
A naive implementation would be something like this: int i_init =2 int j_init =5 int k_init = 10 A = rand(i_init,k_init) B = rand(k_init,j_init) C = zero(k_init,k_init) for i = 1 : i_init for j = 1 : j_init sum = 0 for k = 1 to k_init sum = sum + A(i,k) * B(k,j) C(i,j) = sum
If we were to ask the running time of this algorithm, it would be i * j * k We can see that this is cubic in nature if we set i=j=k. so n^3
Now, there exists a much faster implementation called Strassen_algorithm: https://en.wikipedia.org/wiki/Strassen_algorithm This runs in N^2.8074
The reason why I am telling you this is because matlab does these optimizations for you when you use the multiplication operator (*)
Also to add, there are other algorithms to determine if you have a sparse matrix( meaning, lots of zeros) to speed up your algorithms as well. All this is taken care for you if you can use matrix multiplication. ex: Section 3.5 (sparse vectors) (http://www.cs.princeton.edu/courses/archive/fall16/cos226/lectures/34HashTables+35SearchingApplications.pdf)
As for a guide to reducing loops, here is some advice. (It seems like you are doing it already in some sense) 1. Match the matrix dimensions. Use the matrix multiplication. This is just linear algebra. 2. Use the transpose operator(') to get the right dimensions 3. Use zero, ones, diag, and identity matrix for augmenting your data. 4. Matrix with multiple dim A * Matrix with multiple B. You have to iterate and do it at a vector level. ex psudeo code) foreach v: Vector in A:Matrix { v * B }
You can do the following, send me a private chat and we can discuss more in detail about what you are doing. Check in your code to github, so I can see your current implementation.
|
On October 01 2016 13:16 Aerisky wrote:Show nested quote +On September 30 2016 18:33 maybenexttime wrote:On September 30 2016 14:46 Aerisky wrote: Yeah, normally I've applied to a ton of places for internships (100% agree that it's a numbers game). Now that I'm a soon-to-be new grad, though, I have a return offer from the company at which I interned this summer; I really like the offer so I've only been applying to really cool places where I'd be interested in working/prestigious places. Unless you are 100% sure that the company at which you interned will hire you, don't change your application strategy. I did this mistake. I was an intern in one company for a year, even did the research for my Master's thesis there, and I got accepted for the next wave of their graduate program after I would graduate (I still had one year of university after finishing the internship) or alternatively I could start a PhD at one of the universities they had an agreement with. Then when I was back at the university - a month or so before graduation last Summer - I was told by my previous supervisor that the new management decided to discontinue both programs... Never put all your eggs in one basket. If you get accepted by more companies, you can always decline. Wow, that sucks. That's actually really fucked up on that company's part, damn. And yeah, I'll be applying to more places as well I suppose 
Well, the decision was made several levels of management above my former supervisor and in another country to boot - by the new CTO. The latter wouldn't even be aware that I existed. All in all, I am mostly to blame here. I was naive in thinking the new management wouldn't make such drastic changes (out of four waves of graduates, there was only one person whose hiring they considered a mistake, and everyone else was very successful within the company). Like you said, it's a numbers game.
|
I have a major brainfart at the moment and can not think about this issue for the life of me: Is there any way I can implement this function without a loop?
public int f(int x) { int cost = 5; int count = 0; while (x >= cost) { x -= cost; cost += 5; count++; } return count; } It gives diminishing returns for growing X.
|
On October 01 2016 20:31 RoomOfMush wrote:I have a major brainfart at the moment and can not think about this issue for the life of me: Is there any way I can implement this function without a loop? public int f(int x) { int cost = 5; int count = 0; while (x >= cost) { x -= cost; cost += 5; count++; } return count; } It gives diminishing returns for growing X.
You can use recursion.
int f(int x, int cost, int count) { if (x < cost) { return count; }
return f(x - cost, cost + 5, count + 1); }
|
|
|
|