直接二进制搜索（Direct Binary Search, DBS）方法是一种暴力式搜索方法。由于未引入随机化过程，对于同一初始化结构，经过相同的迭代次数，其优化结果一致。因此，使用该方法想要收敛至全局最优解与初始化结构息息相关。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
import numpy as np
import matplotlib as mpl
mpl.rcParams['figure.dpi']=100
mpl.use('agg')
import matplotlib.pylab as plt
import ceviche
from skimage.draw import disk
from ceviche.modes import insert_mode
from ceviche import fdfd_ez, jacobian
import autograd.numpy as npa
from scipy import constants
# from ceviche import fdtd as fdfd_ez
import multiprocessing
import collections
import json
# Create a container for our slice coords to be used for sources and probes
Slice = collections.namedtuple('Slice', 'x y')

# class opt_structure():

def viz_sim_init(epsr, source1, source2):
    """Solve and visualize a simulation with permittivity 'epsr'
    """   
    simulation1 = fdfd_ez(omega1, dl, epsr, [Npml, Npml])
    _, _, Ez1 = simulation1.solve(source1)
    simulation2 = fdfd_ez(omega2, dl, epsr, [Npml, Npml])
    _, _, Ez2 = simulation2.solve(source2)
    E01 = mode_overlap(Ez1, probe1)
    E02 = mode_overlap(Ez2, probe2)
    return (E01,E02)

def viz_sim(epsr, source1, source2, poles_list):
    """Solve and visualize a simulation with permittivity 'epsr'
    """   
    simulation1 = fdfd_ez(omega1, dl, epsr, [Npml, Npml])
    _, _, Ez1 = simulation1.solve(source1)
    simulation2 = fdfd_ez(omega2, dl, epsr, [Npml, Npml])
    _, _, Ez2 = simulation2.solve(source2)
    try:
        pid = multiprocessing.current_process().pid
    except:
        pid = "none_pid"

    # fig, ax = plt.subplots(1, 3, constrained_layout=True, figsize=(9,3))
    current_time = time.localtime()
    str_name = "WDM_{}_{}_{}_{}_{}_{}_{}_{}_{}".format(
                                                 int(np.around(lambda1*10e8)),
                                                 int(np.around(lambda2*10e8)),
                                                 current_time[0],
                                                 current_time[1],
                                                 current_time[2],
                                                 current_time[3],
                                                 current_time[4],
                                                 current_time[5],
                                                 pid,
                                                 )
    ceviche.viz.abs(Ez1, outline=epsr, cbar=False)
    plt.axis("off")
    plt.savefig(r"./png/{}_1_0.png".format(str_name))
    plt.close()
    ceviche.viz.abs(Ez1, cbar=False)
    plt.axis("off")
    plt.savefig(r"./png/{}_1.png".format(str_name))
    plt.close()
    ceviche.viz.abs(Ez2, outline=epsr, cbar=False)
    plt.axis("off")
    plt.savefig(r"./png/{}_2_0.png".format(str_name))
    plt.close()
    ceviche.viz.abs(Ez2, cbar=False)
    plt.axis("off")
    plt.savefig(r"./png/{}_2.png".format(str_name))
    plt.close()
    # ceviche.viz.abs(epsr, ax=ax[2], cmap='Greys')
    ceviche.viz.abs(epsr, cmap='Greys')
    plt.axis("off")
    plt.savefig(r"./png/{}_3.png".format(str_name))
    plt.close()
    # for sl in slices:
    #     ax[0].plot(sl.x*np.ones(len(sl.y)), sl.y, 'w-', alpha=0.5)
    #     ax[1].plot(sl.x*np.ones(len(sl.y)), sl.y, 'w-', alpha=0.5)
    
    # ax[0].set_title(r'$\lambda_1$ = %.2f $\mu$m' % (constants.c/(omega1/2/np.pi)/1e-6))
    # ax[1].set_title(r'$\lambda_2$ = %.2f $\mu$m' % (constants.c/(omega2/2/np.pi)/1e-6))
    fom = 1/(mode_overlap(Ez1, probe1) / E01 * mode_overlap(Ez2,probe2) / E02)
    recover_json  = open(r"./png/{}.json".format(str_name), "w", encoding="utf-8")
    recover_Chrom = {}
    recover_Chrom[str_name] = [list(poles_list), float(fom)]
    json.dump(recover_Chrom, recover_json, indent = 4)
    recover_json.close()
    return fom

def init_structure(Nx, Ny, Nwg, Nwd, Nox, Noy, Npml, edge_N, index_font, index_background):
    epsr = np.ones((Nx, Ny)) * (index_font ** 2)
    epsr[0:Nwg, (Ny-Nwd)//2:(Ny+Nwd)//2] = index_background ** 2
    epsr[Nwg:Nwg+Nox, (edge_N) * Npml:(edge_N) * Npml+Noy] = index_background ** 2
    epsr[Nwg+Nox:2*Nwg+Nox,(edge_N) * Npml: (edge_N) * Npml+Nwd] = index_background ** 2
    epsr[Nwg+Nox:2*Nwg+Nox,(edge_N) * Npml+Noy-Nwd:(edge_N) * Npml+Noy] = index_background ** 2
    input_slice = Slice(x=np.array(Npml+Nwg//2), 
        y=np.arange((Ny-Nwd)//2- Nwd, (Ny+Nwd)//2+Nwd))
    output_slice1 = Slice(x=np.array(Nwg+Nox+Nwg//2), 
        y=np.arange(edge_N*Npml+Noy-2*Nwd, edge_N*Npml+Noy+Nwd))
    output_slice2 = Slice(x=np.array(Nwg+Nox+Nwg//2), 
        y=np.arange( edge_N*Npml-Nwd, edge_N*Npml+2*Nwd))
    return epsr, input_slice, output_slice1, output_slice2

def init_opt_structure(Nwg, circule_outer_radius, Npml, circule_inter_radius, epsr, index_font, index_background):
    m = Nox // (circule_outer_radius * 2)
    n = Noy // (circule_outer_radius * 2)
    poles = []
    for i in range(m):
        for j in range(n):
            rr, cc = disk((Nwg+2*circule_outer_radius*(j+0.5), (edge_N)*Npml+2*circule_outer_radius*(i+0.5)), circule_inter_radius)
            poles.append([rr,cc])
            epsr[rr,cc] = (index_font ** 2)
    return epsr, poles

def poles2list(poles, epsr, index_font, index_background):
    poles_list = []
    for i in range(len(poles)):
        if epsr[poles[i][0], poles[i][1]][0]   == index_font **2:
            poles_list.append(1)
        elif epsr[poles[i][0], poles[i][1]][0] == index_background **2:
            poles_list.append(0)
    return poles_list

def list2poles(poles_list, pole_index, poles, epsr, index_font, index_background):
    if pole_index == "full":
        for i in range(len(poles_list)):
            if poles_list[i] == 1:
                epsr[poles[i][0], poles[i][1]] = index_font ** 2
            elif poles_list[i] == 0:
                epsr[poles[i][0], poles[i][1]] = index_background ** 2
    else: 
        if poles_list[pole_index] == 1:
            epsr[poles[pole_index][0], poles[pole_index][1]] = index_font ** 2
        elif poles_list[pole_index] == 0:
            epsr[poles[pole_index][0], poles[pole_index][1]] = index_background ** 2
    return epsr

# DBS 算法翻转
def reverse_opt_structure(poles_list, poles, pole_index, epsr, index_font, index_background):
    # if epsr[poles[pole_index][0], poles[pole_index][1]][0]   == index_font **2:
    #     epsr[poles[pole_index][0], poles[pole_index][1]]     = index_background ** 2
    # elif epsr[poles[pole_index][0], poles[pole_index][1]][0] == index_background **2:
    #     epsr[poles[pole_index][0], poles[pole_index][1]]     = index_font ** 2
    if poles_list[pole_index] == 1:
        epsr[poles[pole_index][0], poles[pole_index][1]]     = index_background ** 2
        poles_list[pole_index] = 0
    elif poles_list[pole_index] == 0:
        epsr[poles[pole_index][0], poles[pole_index][1]]     = index_font ** 2
        poles_list[pole_index] = 1
    return epsr

def mode_overlap(E1, E2):
    return npa.abs(npa.sum(npa.conj(E1)*E2))*1e6

def objective(epsr, i="test", j="test"):
    Ez1, Ez2, str_name = viz_sim(epsr, source1, source2)
    # ceviche.viz.abs(epsr)
    plt.savefig(r"./log/test_{}_{}.png".format(i, j))
    plt.close()
    return mode_overlap(Ez1, probe1) / E01 * mode_overlap(Ez2,probe2) / E02

import time
def obj_cfunc(p):
    global epsr
    poles_list = p
    epsr = list2poles(poles_list, "full", poles, epsr, index_font, index_background)
    fom = viz_sim(epsr, source1, source2, poles_list)
    return fom

# if __name__ == "__main__":
# print("test")
# user_define
lambda1 = 1550e-9
lambda2 = 980e-9
omega1  = 2 * np.pi * (constants.c / lambda1)
omega2  = 2 * np.pi * (constants.c / lambda2)
dl = 40e-9
opt_size_x = 4800e-9
opt_size_y = 4800e-9
wg_len   = 2000e-9
wg_width = 500e-9
index_background = 3.47
index_font       = 1.22
outer_radius = 120e-9
inter_radius = 100e-9
pml_width    = 400e-9
edge_width   = 80e-9

Npml                 = int(np.around(pml_width/dl))
edge_N               = int(np.around(edge_width/dl))
circule_outer_radius = int(np.around(outer_radius/dl))
circule_inter_radius = int(np.around(inter_radius/dl))
Nx = int((opt_size_x + 2 * wg_len)*10e9 / (dl*10e9))
Ny = int((opt_size_y)*10e9 / (dl*10e9) + 2 * edge_N * Npml)
Nox = int((opt_size_x*10e9) / (dl*10e9))
Noy = int((opt_size_y*10e9) / (dl*10e9))
Nwg = int(wg_len*10e9 / (dl*10e9))
Nwd = int(wg_width*10e9 / (dl*10e9))

epsr, input_slice, output_slice1, output_slice2 = init_structure(Nx, Ny, Nwg, Nwd, Nox, Noy, Npml, edge_N, index_font, index_background)
epsr, poles = init_opt_structure(Nwg, circule_outer_radius, Npml, circule_inter_radius, epsr, index_font, index_background)
poles_list = poles2list(poles, epsr, index_font, index_background)

# ceviche.viz.abs(epsr, cbar=True)

source1 = insert_mode(omega1, dl, input_slice.x, input_slice.y, epsr, m=1)
source2 = insert_mode(omega2, dl, input_slice.x, input_slice.y, epsr, m=1)

probe1 = insert_mode(omega1, dl, output_slice1.x, output_slice1.y, epsr, m=1)
probe2 = insert_mode(omega2, dl, output_slice2.x, output_slice2.y, epsr, m=1)

E01, E02 = viz_sim_init(epsr, source1, source2)

fom0 = viz_sim(epsr, source1, source2, poles_list)

# DBS 算法采集数据
iter_num = 100
for j in range(iter_num):
    for i in range(len(poles)):
        epsr = reverse_opt_structure(poles_list, poles, i, epsr, index_font, index_background)
        fom1 = objective(epsr, i , j)
        if fom1 <= fom0:
            epsr = reverse_opt_structure(poles_list, poles, i, epsr, index_font, index_background)
        else:
            fom0 = fom1
        print(fom0)

基于随机化的直接二进制搜索方法可能会降低对初始化结构的需求，并且提升了算法收敛至全局最优解的可能性。

⚓ Carl Zhao
🏢 逍遥科技有限公司
💭 曾经也是追光少年，然而少年归来已不再是少年，但依然在追光的路上。
📧 邮箱：1005513510@qq.com