review all files, write more comments

Author: deng-cy

fluid_coarse/bba.py

@@ -1,22 +1,24 @@
-# https://github.com/Sanbongawa/binary_swarm_intelligence
+# reference: https://github.com/Sanbongawa/binary_swarm_intelligence
 
 import math
 
 import torch
 import tqdm
 
 
+@torch.no_grad()
 def bba(func, x0=None, n=200, maxiter=2000, topk=100, is_min=True, fmin=0, fmax=2, loud_A=0.25, r0=0.4, alpha=0.999, gamma=10., progress=False):
     """
+    Binary Bat Algorithm
     input:{ func: Evaluate_Function, type is nn.Module
             x0: initial guess of the solution
-            n: Number of population, default=20
-            maxiter: Number of max iteration, default=300
+            n: Number of population
+            maxiter: Number of max iteration
             topk: the number of states to output (see output best_s_list)
             is_min: minimizing (True) or maximizing (False) func
             fmin: a hyperparamter, frequency minimum to step
             fmax: a hyperparamter, frequency maximum to step
-            loud_A: a hyperparamter, value of Loudness, default=0.25
+            loud_A: a hyperparamter, value of Loudness
             r0: a hyperparamter, pulse rate, Probability to relocate near the best position
             alpha: a hyperparamter, decay rate of loud_A
             gamma: a hyperparamter, decay rate of gamma
@@ -29,28 +31,28 @@ def bba(func, x0=None, n=200, maxiter=2000, topk=100, is_min=True, fmin=0, fmax=
             }
     """
 
-    device = next(func.parameters()).device
+    device = next(func.parameters()).device  # get device(CPU/GPU) from `func`
     dim = x0.shape[1]
 
+    # initialize
     v = torch.zeros(n, dim, device=device)
-
     x = torch.rand(n - 1, dim)
-    x = (x > 0.9).type_as(x0)
+    x = (x > 0.9).type_as(x0)  # convert to binary
     x = torch.cat([x, x0], dim=0)
 
-    torch.set_grad_enabled(False)
     fit = func(x) if is_min else -func(x)
     min_fit = fit.min()
     best_v = min_fit
     best_s = x[fit.argmin()].reshape((1, -1))
-    best_list = []
+    best_v_list = []
 
     if progress:
         m = tqdm.tqdm(range(maxiter))
     else:
         m = range(maxiter)
     for t in m:
 
+        # update parameters
         r = r0 * (1 - torch.exp(torch.tensor(-gamma * t, device=device)))
         loud_A *= alpha
 
@@ -61,10 +63,10 @@ def bba(func, x0=None, n=200, maxiter=2000, topk=100, is_min=True, fmin=0, fmax=
         x_new = x.clone()
 
         # flip
-        accept = torch.rand(n, dim, device=device) < vstf + 1 / 160
+        accept = torch.rand(n, dim, device=device) < vstf + 1 / dim
         x_new[accept] = 1. - x_new[accept]
 
-        # local opt
+        # accept local optimum
         accept = torch.rand(n, dim, device=device) > r
         x_new[accept] = best_s.expand(n, dim)[accept]
 
@@ -77,14 +79,15 @@ def bba(func, x0=None, n=200, maxiter=2000, topk=100, is_min=True, fmin=0, fmax=
         x[accept, :] = x_new[accept, :]
         fit[accept] = fit_new[accept]
 
+        # update best solution
         min_fit, argmin_f = fit.min(0)
-
         if min_fit < best_v:
             best_s = x[argmin_f].reshape((1, -1)).clone()
             best_v = min_fit.clone()
 
-        best_list.append(best_v if is_min else - best_v)
+        best_v_list.append(best_v if is_min else - best_v)
 
-    best_s_list = torch.unique(torch.cat((fit, x), dim=1), dim=0)[0:topk, :]
+    # postprcess, get ready to output
+    best_s_list = torch.unique(torch.cat((fit, x), dim=1), dim=0)[0:topk, :]  # delete repeated solutions
 
-    return (best_v if is_min else - best_v), best_s, best_s_list, best_list
+    return (best_v if is_min else - best_v), best_s, best_s_list, best_v_list

fluid_coarse/func_inputs_gen.m

@@ -1,16 +1,22 @@
 function inputs = func_inputs_gen(n,seed_gen,base,nx,ny,p,inputs_ori,n_coeff)
-% generate inputs by three ways based on args
+% generate new inputs (samples) by three ways based on args
 % 1) if `base` is empty, meaning the inital step, then 161 arrays will be
 % generated
 % 2) if `base` is not empty, `p` is empty, employ SOLO-G, return base
-% 3) if `base` is not empty, `p` is not empty, employ SOLO-R, add
-% disturbance
+% 3) if `base` is not empty, `p` is not empty, employ SOLO-R, add disturbance
+
+% Inputs:
+%   iter: number of total samples that needs to output
+%   base: base design, i.e., the optimized solution from DNN in SOLO
+%   p: empty or [a1,a2,a3,...,an] where `ai` (i<n) denotes the probability to mutate a ixi matrix, and `an` denotes convolution possbility
+%   inputs_ori: existing samples, used to remove repeated ones of new samples  
+%   n_coeff: a coefficient to generate more samples than needed. Since repeated samples will be deleted later, more samples are needed at the generation stage
 rng(seed_gen);
 base_empty=isempty(base);
 alpha=4;
 
 if ~exist('n_coeff','var')
-    n_coeff=100; % used when p exist
+    n_coeff=100; % set a default value, used when p exist
 end
 
 

fluid_coarse/main_greedy.m

@@ -12,15 +12,14 @@
 n=10; % number of additional samples per iteration
 n_0=161; % number of samples in the initial batch
 results(:,1)=(n:n:n*iter)';
-seed_train=3; % DNN training random seed
-seed_opt=123; % BBA optimization random seed
+seed_train=3; % random seed for DNN training
+seed_opt=123; % random seed for BBA optimization
 nx=a/delta; % number of design variables in x direction
 ny=b/delta; % number of design variables in y direction
 %% loops
 
 for i=1:iter
     % get input/output file names
-    testfile=['./',folder,'/',num2str(n),'ntrain_test_step',num2str(i),'.mat'];
     optfile=['./',folder,'/',num2str(n),'ntrain_opt_step',num2str(i),'.mat'];
     trainfile=['./',folder,'/',num2str(n),'ntrain_train_step',num2str(i),'.mat'];
 
@@ -31,7 +30,8 @@
     % eliminate invalid data
     inputs(outputs==Inf,:,:)=[];
     outputs(outputs==Inf)=[];
-    if i>1
+    
+    if i>1 % if not initial batch, combine data
         mydata=load(['./',folder,'/',num2str(n),'ntrain_train_step',num2str(i-1),'.mat']);
         inputs=[mydata.inputs;inputs];
         outputs=[mydata.outputs;outputs];

fluid_coarse/main_regular.m

@@ -17,10 +17,9 @@
 nx=a/delta; % number of design variables in x direction
 ny=b/delta; % number of design variables in y direction
 %% loops
-p=[0.1,0.1,0.2,0.2,0.2];
+p=[0.1,0.1,0.2,0.2,0.2]; % [a1,a2,a3,...,an] where `ai` (i<n) denotes the probability to mutate a ixi matrix, and `an` denotes crossover possbility
 for i=1:iter
     % get input/output file names
-    testfile=['./',folder,'/',num2str(n),'ntrain_test_step',num2str(i),'.mat'];
     optfile=['./',folder,'/',num2str(n),'ntrain_opt_step',num2str(i),'.mat'];
     trainfile=['./',folder,'/',num2str(n),'ntrain_train_step',num2str(i),'.mat'];
 

fluid_coarse/plot_b_d_e.m

@@ -8,6 +8,7 @@
 %% set parameters and preprocessing
 g=round(p0./g,6); % there are equivalent solutions with very slight different outputs, use round to eliminate numerical error
 r=p0./r;
+font='Times';
 %% generate some beautiful colors
 colors=func_linspecer(8);
 %% energy (objective)
@@ -18,16 +19,16 @@
 plot(func_best_min(r),'-','Color',colors(7,:),'LineWidth',3);hold on;
 
 xl=xlabel('$n_{train}$','Interpreter','latex');
-yl=ylabel('$\tilde{P}$','Interpreter','latex','Rotation',0);
+yl=ylabel('$\widetilde{P}$','Interpreter','latex','Rotation',0);
 [y,x]=min(g);
 scatter(x,y,300,colors(6,:),'x','LineWidth',2); hold on;
 [y,x]=min(r);
 scatter(x,y,300,colors(7,:),'x','LineWidth',2); hold on;
 
-lgd=legend('Gradient-based','SOLO-G','SOLO-R','Interpreter','latex');
+lgd=legend('Gradient-based','SOLO-G','SOLO-R');
 lgd.Position=([0.6529 0.7008 0.3305 0.2692]);
 
-set(gca,'FontSize',fontsize,'Fontname', 'Times New Roman', 'Position',[0.1,0.13,0.7750,0.7950])
+set(gca,'FontSize',fontsize,'Fontname', font, 'Position',[0.1,0.13,0.7750,0.7950])
 xl.Position=[2.7169e+03 0.9605 -1];
 yl.Position=[-222.5806 0.9971 -1];
 ylim([0.955,1])
@@ -69,7 +70,7 @@
 [val, idx]=min(g);
 inputs=gdata.inputs;
 input=inputs(idx,:,:);
-h=heatmap(squeeze(1-input)','CellLabelColor','none','FontName','Times New Roman','FontSize',20);
+h=heatmap(squeeze(1-input)','CellLabelColor','none','FontName',font,'FontSize',20);
 colormap(gray)
 colorbar off
 h.Position=[0.1,0.2,0.8,0.8];
@@ -82,7 +83,7 @@
 [val, idx]=min(r);
 inputs=rdata.inputs;
 input=inputs(idx,:,:);
-h=heatmap(squeeze(1-input)','CellLabelColor','none','FontName','Times New Roman','FontSize',20);
+h=heatmap(squeeze(1-input)','CellLabelColor','none','FontName',font,'FontSize',20);
 colormap(gray)
 colorbar off
 h.Position=[0.1,0.2,0.8,0.8];

fluid_coarse/plot_c.m

@@ -1,11 +1,13 @@
 %% comsol server required, you can jump to `plot` if data is downloaded
-model = mphload('fluid_optimized.mph');
-data=mpheval(model,{'dtopo1.theta'},'selection',2,'edim','domain','dataset','dset2');
 a=2;
 b=0.8;
-delta=0.1;
+delta=0.05;
 nx=a/delta;
 ny=b/delta;
+
+% extract nodal values
+model = mphload('fluid_optimized.mph');
+data=mpheval(model,{'dtopo1.theta'},'selection',2,'edim','domain','dataset','dset2');
 indices=round(data.p./delta+1);
 input_node=zeros(nx+1,ny+1);
 for i=1:size(indices,2)
@@ -33,11 +35,13 @@
 outputs = func_outputs_nowrite(inputs_unique,nx,ny);
 [val, idx]=max(outputs);
 save('data_gb.mat','idx','outputs','inputs_unique');
-% 1.557373167442459
+% val=1.557373167442459
 %% plot
+font='Times';
+load('data_gb.mat');
 figure('Position',[0,0,500,200])
 input=inputs_unique(idx,:,:);
-h=heatmap(squeeze(1-input)','CellLabelColor','none','FontName','Times New Roman','FontSize',20);
+h=heatmap(squeeze(1-input)','CellLabelColor','none','FontName',font,'FontSize',20);
 colormap(gray)
 colorbar off
 h.Position=[0.1,0.2,0.8,0.8];

fluid_coarse/readme.md

@@ -53,11 +53,13 @@ from [Google Drive](https://drive.google.com/drive/folders/1f6Xrd9e-RAUsh9vqIqUX
 ## How to run
 
 1. Connect MATLAB with COMSOL server
-2. Run __fluid.m__ by COMSOL to generate __fluid.mph__
+2. Run __fluid.m__ by COMSOL to generate __fluid.mph__, or download it
+   from [Google Drive](https://drive.google.com/drive/folders/1f6Xrd9e-RAUsh9vqIqUXbEw8F1_2Qg_5?usp=sharing)
 3. Set up MATLAB such that it can call python scripts ([instructions](https://www.mathworks.com/help/matlab/call-python-libraries.html))
 4. Run __main_greedy.m__ and __main_regular.m__
 5. [optional] Run __plot_b_d_e.m__ and __plot_s\[4-6\].m__
-6. Run __fluid_optimized.m__ by _COMSOL_ to generate __fluid_optimized.mph__
+6. Run __fluid_optimized.m__ by _COMSOL_ to generate __fluid_optimized.mph__, or download it
+   from [Google Drive](https://drive.google.com/drive/folders/1f6Xrd9e-RAUsh9vqIqUXbEw8F1_2Qg_5?usp=sharing)
 7. [optional] Run __plot_c.m__
 
 As long as a __*.m__ file says `Comsol server required` in the comment at the top, you need to connect the MATLAB with a COMSOL server to run the

fluid_fine/bba.py

@@ -1,22 +1,24 @@
-# https://github.com/Sanbongawa/binary_swarm_intelligence
+# reference: https://github.com/Sanbongawa/binary_swarm_intelligence
 
 import math
 
 import torch
 import tqdm
 
 
+@torch.no_grad()
 def bba(func, x0=None, n=200, maxiter=2000, topk=100, is_min=True, fmin=0, fmax=2, loud_A=0.25, r0=0.4, alpha=0.999, gamma=10., progress=False):
     """
+    Binary Bat Algorithm
     input:{ func: Evaluate_Function, type is nn.Module
             x0: initial guess of the solution
-            n: Number of population, default=20
-            maxiter: Number of max iteration, default=300
+            n: Number of population
+            maxiter: Number of max iteration
             topk: the number of states to output (see output best_s_list)
             is_min: minimizing (True) or maximizing (False) func
             fmin: a hyperparamter, frequency minimum to step
             fmax: a hyperparamter, frequency maximum to step
-            loud_A: a hyperparamter, value of Loudness, default=0.25
+            loud_A: a hyperparamter, value of Loudness
             r0: a hyperparamter, pulse rate, Probability to relocate near the best position
             alpha: a hyperparamter, decay rate of loud_A
             gamma: a hyperparamter, decay rate of gamma
@@ -29,28 +31,28 @@ def bba(func, x0=None, n=200, maxiter=2000, topk=100, is_min=True, fmin=0, fmax=
             }
     """
 
-    device = next(func.parameters()).device
+    device = next(func.parameters()).device  # get device(CPU/GPU) from `func`
     dim = x0.shape[1]
 
+    # initialize
     v = torch.zeros(n, dim, device=device)
-
     x = torch.rand(n - 1, dim)
-    x = (x > 0.9).type_as(x0)
+    x = (x > 0.9).type_as(x0)  # convert to binary
     x = torch.cat([x, x0], dim=0)
 
-    torch.set_grad_enabled(False)
     fit = func(x) if is_min else -func(x)
     min_fit = fit.min()
     best_v = min_fit
     best_s = x[fit.argmin()].reshape((1, -1))
-    best_list = []
+    best_v_list = []
 
     if progress:
         m = tqdm.tqdm(range(maxiter))
     else:
         m = range(maxiter)
     for t in m:
 
+        # update parameters
         r = r0 * (1 - torch.exp(torch.tensor(-gamma * t, device=device)))
         loud_A *= alpha
 
@@ -62,9 +64,11 @@ def bba(func, x0=None, n=200, maxiter=2000, topk=100, is_min=True, fmin=0, fmax=
 
         # flip
         accept = torch.rand(n, dim, device=device) < vstf + 1 / 160
+        # TODO: I made a mistake here, `160` above should be `dim` to match the paper (although I do not think it matters).
+        #  I kept this mistake to help repeat my result, but you can change it for new problems
         x_new[accept] = 1. - x_new[accept]
 
-        # local opt
+        # accept local optimum
         accept = torch.rand(n, dim, device=device) > r
         x_new[accept] = best_s.expand(n, dim)[accept]
 
@@ -77,14 +81,15 @@ def bba(func, x0=None, n=200, maxiter=2000, topk=100, is_min=True, fmin=0, fmax=
         x[accept, :] = x_new[accept, :]
         fit[accept] = fit_new[accept]
 
+        # update best solution
         min_fit, argmin_f = fit.min(0)
-
         if min_fit < best_v:
             best_s = x[argmin_f].reshape((1, -1)).clone()
             best_v = min_fit.clone()
 
-        best_list.append(best_v if is_min else - best_v)
+        best_v_list.append(best_v if is_min else - best_v)
 
-    best_s_list = torch.unique(torch.cat((fit, x), dim=1), dim=0)[0:topk, :]
+    # postprcess, get ready to output
+    best_s_list = torch.unique(torch.cat((fit, x), dim=1), dim=0)[0:topk, :]  # delete repeated solutions
 
-    return (best_v if is_min else - best_v), best_s, best_s_list, best_list
+    return (best_v if is_min else - best_v), best_s, best_s_list, best_v_list

fluid_fine/main_greedy.m

@@ -7,20 +7,19 @@
 delta=0.05; % grid size of design variables
 
 base=[]; % initial guess of solution
-iter=200; % total number of iterations of the algorithm
+iter=200; % number of total iterations of the algorithm
 results=zeros(iter,3); % initialize result table
 n=10; % number of additional samples per iteration
 n_0=641; % number of samples in the initial batch
 results(:,1)=(n:n:n*iter)';
-seed_train=3; % DNN training random seed
-seed_opt=123; % BBA optimization random seed
+seed_train=3; % random seed for DNN training
+seed_opt=123; % random seed for BBA optimization
 nx=a/delta; % number of design variables in x direction
 ny=b/delta; % number of design variables in y direction
 
 %% loops
 for i=1:iter
     % get input/output file names
-    testfile=['./',folder,'/',num2str(n),'ntrain_test_step',num2str(i),'.mat'];
     optfile=['./',folder,'/',num2str(n),'ntrain_opt_step',num2str(i),'.mat'];
     trainfile=['./',folder,'/',num2str(n),'ntrain_train_step',num2str(i),'.mat'];
 
@@ -31,8 +30,8 @@
     % eliminate invalid data
     inputs(outputs==Inf,:,:)=[];
     outputs(outputs==Inf)=[];
-    %
-    if i>1
+
+    if i>1 % if not initial batch, combine data
         mydata=load(['./',folder,'/',num2str(n),'ntrain_train_step',num2str(i-1),'.mat']);
         inputs=[mydata.inputs;inputs];
         outputs=[mydata.outputs;outputs];