forked from BarcelonaMedia-Audio/idhoa
-
Notifications
You must be signed in to change notification settings - Fork 0
/
main.py
314 lines (261 loc) · 13.7 KB
/
main.py
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
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
#!/usr/bin/python
'''
* This file is part of IDHOA software.
*
* Copyright (C) 2013 Barcelona Media - www.barcelonamedia.org
* (Written by Davide Scaini <[email protected]> for Barcelona Media)
*
* The clone https://github.com/davrandom/idhoa and its modifications are under
* Copyright (C) 2014 Davide Scaini <[email protected]>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>,
* or write to the Free Software Foundation, Inc.,
* 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
'''
import nlopt
# on Mac:
# export PYTHONPATH=$PYTHONPATH:/usr/local/lib/python2.7/site-packages
import time
from argparse import ArgumentParser
import logging
import numpy as np
import pickle
import os
import ini_parser as prs
import functions as fct
import plotting as ptt
# parse cli arguments
parser = ArgumentParser(description="Generate Higher Order Ambisonics decoding for Irregular Layouts")
parser.add_argument('configfile', metavar='configfile', type=str,
help="Give a configuration file to calculate your coefficients. (Look at the example.ini)")
args = parser.parse_args()
# logging
LOGGING_LEVEL = logging.DEBUG
logging.basicConfig(
format='%(filename)s: %(lineno)d: %(levelname)s: %(message)s',
level=LOGGING_LEVEL)
# initialize configuration
configfile = args.configfile
cfg = prs.ConfigConstants(configfile)
ambicl = fct.Ambisonic(cfg)
support = fct.Support(cfg, ambicl)
logging.info("You choose %s at order: %d " % (cfg.DEC, cfg.DEG))
logging.info("The number of speakers is %d " % cfg.NSPK)
logging.info("Number of sampling points %d " % cfg.NPOINTS)
# INITIAL PARAMETERS
# calculating ambisonics coefficients in test directions
coffs_in_test_directions = ambicl.get_ambisonic_coeffs(phi=cfg.phiTest, theta=cfg.thetaTest, DEC='basic', DEG=cfg.DEG, inversion=0)
support.coffs_in_test_directions = coffs_in_test_directions
initial_guess_projection = ambicl.get_ambisonic_coeffs(inversion=0) # initial guess
initial_guess_pseudoinv = ambicl.get_ambisonic_coeffs(inversion=1) # initial guess
obj_function_at_proj_guess = support.objective_function(initial_guess_projection, coffs_in_test_directions[:])
obj_function_at_pinv_guess = support.objective_function(initial_guess_pseudoinv, coffs_in_test_directions[:])
initial_best_guess = initial_guess_projection
if obj_function_at_proj_guess > obj_function_at_pinv_guess and np.asarray([abs(initial_guess_pseudoinv) < 1.]).all():
initial_best_guess = initial_guess_pseudoinv # This way we choose the starting point closest to the minimum
sij = support.Sij(initial_best_guess, coffs_in_test_directions[:], cfg.NPOINTS)
pressure0, V0, energyD0, J0, Vradial0, Jradial0, Vtang0, Jtang0 = support.directional_components(sij, cfg.phiTest, cfg.thetaTest)
######################
# Initial PLOTTING #
######################
fct.controls(cfg.NSPK, cfg.DEG) # checks before running the main program
# Plot speakers
ptt.threed_polar_plot(cfg.PHI, cfg.THETA, 1, True)
# Plot evaluations points grid
if cfg.exclude_with_theta_binary_mask:
ptt.threed_polar_plot(cfg.phiTest * cfg.WbinVec, cfg.thetaTest * cfg.WbinVec, 1)
if cfg.autoexclude_regions_with_no_spkrs_smooth or cfg.autoexclude_regions_with_no_spkrs_binary:
ptt.threed_polar_plot(cfg.phiTest * cfg.WremVec, cfg.thetaTest * cfg.WremVec, 1)
# print "\nIf you are happy with the distribution of speakers press a key \nOtherwise abort pressing ctrl+c"
# wait = raw_input()
phi_plot, theta_plot = ptt.points_over_sphere_plotting(cfg.SEED, cfg)
Npt = len(phi_plot)
phi = np.asarray(phi_plot)
theta = np.asarray(theta_plot)
coeffs_plot = ambicl.get_ambisonic_coeffs(phi_plot, theta_plot, 'basic', cfg.DEG, 0)
sij = support.Sij(initial_guess_projection, coeffs_plot, Npt)
pressureZ, VZ, energyDZ, JZ, VradialZ, JradialZ, VtangZ, JtangZ = support.directional_components(sij, phi_plot, theta_plot)
if cfg.DEC == 'basic':
ptt.polar_plot(cfg, "Naive-Horizontal", phi[theta == 0.], pressureZ[theta == 0.], VradialZ[theta == 0.], VtangZ[theta == 0.],
('pressure', 'velocity L', 'velocity T'))
else:
ptt.polar_plot(cfg, "Naive-Horizontal", phi[theta == 0.], energyDZ[theta == 0.], JradialZ[theta == 0.], JtangZ[theta == 0.],
('energy', 'intensity L', 'intensity T'))
if cfg.nD == '3D':
if cfg.DEC == 'basic':
ptt.polar_plot(cfg, "Naive-Vertical", theta[phi == 0], pressureZ[phi == 0], VradialZ[phi == 0], VtangZ[phi == 0],
('pressure', 'velocity L', 'velocity T'))
else:
ptt.polar_plot(cfg, "Naive-Vertical", theta[phi == 0], energyDZ[phi == 0], JradialZ[phi == 0], JtangZ[phi == 0],
('energy', 'intensity L', 'intensity T'))
##################
# MINIMIZATION #
##################
start = time.time()
# here we have to shrink the initGuess
initial_best_guess = support.downscalematrix(initial_best_guess, cfg.MATCHED)
initvect = support.mat2vec(initial_best_guess)
# Global Optimization
# GN_DIRECT, GN_DIRECT_L, GN_ORIG_DIRECT, GN_ORIG_DIRECT_L, GN_DIRECT_NOSCAL, GN_DIRECT_L_NOSCAL, GN_DIRECT_L_RAND_NOSCAL
# GN_CRS2_LM
# G_MLSL_LDS, G_MLSL
# GD_STOGO, GD_STOGO_RAND
# GN_ISRES
prog = []
tprogress = []
nonzeroprog = []
run = 0
minstart = time.time()
while True:
# Local Optimization
# LN_COBYLA, LN_BOBYQA, LN_NEWUOA, LN_NEWUOA_BOUND, LN_PRAXIS, LN_NELDERMEAD, LN_SBPLX
logging.info("ITERATION n. %d" % run)
if run > 0:
ResCoeff = support.downscalematrix(ResCoeff, cfg.MATCHED)
initvect = support.mat2vec(ResCoeff)
opt = nlopt.opt(nlopt.LN_SBPLX, len(initvect))
opt.set_min_objective(support.objective_function)
tol = np.asarray([0.1] * len(initvect))
# opt.add_equality_mconstraint(equality_constr, tol)
# opt.add_inequality_mconstraint(inequality_constr, tol)
opt.set_initial_step([0.2] * len(initvect))
if run > 0: opt.set_initial_step([0.9] * len(initvect))
# opt.set_initial_step(np.random.rand(len(initvect)));
upbound = [1.] * len(initvect)
lowbound = [-1.] * len(initvect)
nodes = ambicl.get_ambisonic_coeffs(cfg.PHI, cfg.THETA, 'basic', cfg.DEG, 0)
nodes = support.downscalematrix(nodes, cfg.MATCHED)
# calculates the coefficients in the basic or naive scheme.
# These values are used to figure out which speakers lay in the nodes of some spherical harmonic
nodes, upbound, lowbound = support.zeroing_bounds(nodes, upbound, lowbound, run)
if run > 0: nodes, upbound, lowbound = support.zeroing_bounds(res, upbound, lowbound, run)
# number of non-zero elements
logging.info( "%d non-zero elements" % len(np.nonzero(initvect)[0]) )
nonzeroprog.append(len(np.nonzero(initvect)[0]))
opt.set_upper_bounds(upbound)
opt.set_lower_bounds(lowbound)
opt.set_xtol_abs(10e-6)
opt.set_ftol_abs(10e-6)
if run > 0:
opt.set_xtol_abs(10e-8)
opt.set_ftol_abs(10e-8)
# getting the NEW minimization matrix
res = opt.optimize(initvect)
result = opt.last_optimize_result() # alternative way to get the result
ResCoeff = support.vec2mat(res)
ResCoeff = support.upscalematrix(ResCoeff, cfg.MATCHED)
obj_function_at_NL_minimum = support.objective_function(ResCoeff, coffs_in_test_directions[:])
logging.info("Objective function value for Naive: %s" % obj_function_at_proj_guess)
logging.info("Objective function value for Pinv: %s" % obj_function_at_pinv_guess)
logging.info("Objective function value for non-linear minimization: %s" % obj_function_at_NL_minimum)
logging.info("Elapsed time: %s" % (time.time() - minstart))
prog.append(support.objective_function(ResCoeff, coffs_in_test_directions[:]))
tprogress.append(time.time() - minstart)
minstart = time.time()
#####################
# exit condition
# if (run>0 and (str(prog[run-1])[0:6]==str(prog[run])[0:6] or prog[run]>min(prog)+1)): break # for PRAXIS use this
if (run > 0 and ("{:7.3f}".format(prog[run - 1]) == "{:7.3f}".format(prog[run]) or
prog[run] > min(prog) + 1) or not cfg.mute_small_coeffs):
break # for SBPLX use this
run += 1
#############
# RESULTS #
#############
logging.info( "Minimization Results:" )
if result == nlopt.SUCCESS: logging.info("Successful minimization")
if result == nlopt.STOPVAL_REACHED: logging.info("Optimization stopped because stopval was reached")
if result == nlopt.FTOL_REACHED: logging.info("Optimization stopped because ftol_rel or ftol_abs was reached")
if result == nlopt.XTOL_REACHED: logging.info("Optimization stopped because xtol_rel or xtol_abs was reached")
if result == nlopt.MAXEVAL_REACHED: logging.info("Optimization stopped because maxeval was reached")
if result == nlopt.MAXTIME_REACHED: logging.info("Optimization stopped because maxtime was reached")
if result == nlopt.FAILURE: logging.error("Minimization FAILED")
if result == nlopt.INVALID_ARGS: logging.error("Invalid arguments (e.g. lower bounds are bigger than upper bounds, "
"an unknown algorithm was specified, etcetera)")
if result == nlopt.OUT_OF_MEMORY: logging.error("Ran out of memory.")
if result == nlopt.ROUNDOFF_LIMITED:logging.error("Halted because roundoff errors limited progress. "
"(In this case, the optimization still typically returns a useful result.)")
if result == nlopt.FORCED_STOP: logging.error("Halted because of a forced termination: "
"the user called nlopt_force_stop(opt) on the optimization's "
"nlopt_opt object opt from the user's objective function or constraints.")
ResCoeff = support.vec2mat(res)
# np.reshape(res, ((DEG+1)**2,NSPKmatch))
ResCoeff = support.upscalematrix(ResCoeff, cfg.MATCHED)
logging.info("LF matrix \n%s" % initial_guess_pseudoinv.T )
logging.info("Coefficients \n%s" % ResCoeff.T)
if ResCoeff.T[abs(ResCoeff.T > 1.)].any(): logging.warning("You reached a bad minimum.")
##############
# PLOTTING #
##############
# results plots
coeffs_plot = ambicl.get_ambisonic_coeffs(phi_plot, theta_plot, 'basic', cfg.DEG, 0)
sij = support.Sij(ResCoeff, coeffs_plot, Npt)
pressure, V, energyD, J, Vradial, Jradial, Vtang, Jtang = support.directional_components(sij, phi_plot, theta_plot)
if cfg.DEC == 'basic':
ptt.polar_plot(cfg, "Horizontal", phi[theta == 0.], pressure[theta == 0.], Vradial[theta == 0.], Vtang[theta == 0.],
('pressure', 'velocity L', 'velocity T'))
else:
ptt.polar_plot(cfg, "Horizontal", phi[theta == 0.], energyD[theta == 0.], Jradial[theta == 0.], Jtang[theta == 0.],
('energy', 'intensity L', 'intensity T'))
if cfg.nD == '3D':
if cfg.DEC == 'basic':
ptt.polar_plot(cfg, "Vertical", theta[phi == 0], pressure[phi == 0], Vradial[phi == 0], Vtang[phi == 0],
('pressure', 'velocity L', 'velocity T'))
else:
ptt.polar_plot(cfg, "Vertical", theta[phi == 0], energyD[phi == 0], Jradial[phi == 0], Jtang[phi == 0],
('energy', 'intensity L', 'intensity T'))
###threed_polar_plot(phi, theta, Jradial)
logging.info("Elapsed time %s" % (time.time() - start))
###############################
# output some files
filename = "%s-%s-rem%s" % (cfg.DEG, cfg.DEC, cfg.autoexclude_regions_with_no_spkrs_binary)
np.savetxt(filename + ".txt", ResCoeff.T, fmt="%f", delimiter=" ")
np.savetxt(filename + "-Gpinv.txt", initial_guess_pseudoinv.T, fmt="%f", delimiter=" ")
np.savetxt(filename + "-G0.txt", initial_guess_projection.T, fmt="%f", delimiter=" ")
# Save state for later reuse with utilities
cfg.projection_guess_matrix = initial_guess_projection.T
cfg.pseudoinv_guess_matrix = initial_guess_pseudoinv.T
cfg.obj_minimization_matrix = ResCoeff.T
outfile = os.path.split(configfile)[-1]
outfile = outfile.split('.')[0]
outfile = "results_for_%s_%s_%s.idhoa" % (outfile, cfg.DEG, cfg.DEC)
pickle.dump(cfg, open(outfile, "wb"))
# save MATLAB output
# patched by Aaron Heller for compatibility with adt
# bitbucket.org/ambidecodertoolbox/adt.git
if cfg.is_mat_out_file:
import scipy.io
import os.path
scipy.io.savemat(os.path.expanduser(cfg.matlab_filename),
{'ResCoeff': ResCoeff,
'Naive': initial_guess_projection,
'Pinv': initial_guess_pseudoinv,
'PHI': cfg.PHI,
'THETA': cfg.THETA,
'DEC': cfg.DEC,
'DEG': cfg.DEG,
'NSPK': cfg.NSPK,
'thetaThreshold': cfg.thetaThreshold,
'phiTest': cfg.phiTest,
'thetaTest': cfg.thetaTest,
'NPOINTS': cfg.NPOINTS,
'WfrontVec': cfg.WfrontVec,
'WplaneVec': cfg.WplaneVec,
'WbinVec': cfg.WbinVec,
'WremVec': cfg.WremVec})
obj_function_at_NL_minimum = support.objective_function(ResCoeff, coffs_in_test_directions[:])
logging.info("Objective function value for Naive: %s" % obj_function_at_proj_guess)
logging.info("Objective function value for Pinv: %s" % obj_function_at_pinv_guess)
logging.info("Objective function value for non-linear minimization: %s" % obj_function_at_NL_minimum)
# wait = raw_input()