HONR269*************************************************************************************
Quantum Computing with Superconducting Circuits
2/15
Began learning about resonance
Learned how to properly use SMA ports and channels
Determined quality factor of resonator chamber
Q = ResonanceFrequency/deltaFrequency
Q = 1671.83 for this resonator
*************************************************************************************
2/21
Steps for calibrating network analyzer:
For signal ports 1 and 2, run
Open circuit
Short circuit
Loaded circuit
Combined transmission
2-port load isolation
When calibration is on, interpolation is acceptable to a certain degree, but extrapolation is very unreliable for this analyzer
How to calculate resonance of cubical chambers in a vacuum:
Fr = c/(2pi) * sqrt((mpi/a)^2 + (npi/b)^2 + (Lpi/d)^2)
Where Fr is resonant frequency
c is speed of light
m, n, l = "Nodes" of resonance
a, b, d = dimensions of cavity
Detectable frequencies of resonator using this equation:
7.0966GHz
9.7634 GHz
12.509 GHz
*************************************************************************************
2/22
Began learning how to use the evaporator
Cleaned a bunch of gold out from the bottom
Started obtaining magnitude and phase return plots for different resonant frequencies of chamber
Calculated difference between expected resonances of the cavity versus what the network analyzer picked up
Ideal: 7.0966
Calculated: 7.34625
%Diff = 3.4%
Ideal: 9.7634
Calculated: 10.0375
%Diff = 2.8%
*************************************************************************************
2/25
Continued collecting data from plots
Will put plots in as soon as they finish
Note that the power output signal had to be lowered from 10 mdB to 7 mdB for reasons unknown to prevent the analyzer from breaking, need to fix
*************************************************************************************
3/1
This function and subsequent script automatically run the network analyzer for all four different measurements
function [sData] = runMultipleMeasurements(power, centerFrequency, span, date)
%
clear sData;
NetAnzAddr = 17; %Network analyzer has GPIB address 17.
SNames = {'S11', 'S12'; 'S21', 'S22'};
for n = 1:2
for m = 1:2
NetAnzObj = gpib('ni', 0, 17); %Define object for GPIB communication.
NetAnzObj.InputBufferSize = 1602*40; %Make input buffer longer than the longest possible data trace (1601 points, 2 values per point, 20 ascii bytes per value)
fopen(NetAnzObj);
fprintf(NetAnzObj, [':CALC1:PAR1:DEF S' num2str(n) num2str(m)]);
fprintf(NetAnzObj, [':CALC1:PAR2:DEF S' num2str(n) num2str(m)]);
fclose(NetAnzObj);
delete(NetAnzObj);
clear NetAnzObj;
[data, ~] = E5071C_multiplepowertraces_singlespan_header(power, power, 1, centerFrequency, span, 1601, 1000, 10, [date '-S' num2str(n) num2str(m) '-' num2str(centerFrequency)]);
sData.(SNames{n, m}).freq = data(:, 1);
sData.(SNames{n, m}).logMag = data(:, 2);
sData.(SNames{n, m}).phase = data(:, 3);
end
end
end
f = 7.345e9;
span = 30e6;
power = 7;
date = '2019-3-1';
sData = runMultipleMeasurements(power, f, span, date);
save(['sDataCal30Span' num2str(f)], '-struct', 'sData')
*************************************************************************************
3/2
This function graphs subplots of all the different measurements the analyzer takes
function graphSubplots(sDataCal, sDataNoCal)
figure
subplot(2, 2, 1);
plot(sDataCal.S11.freq, sDataCal.S11.logMag);
title('Reflection w/Calibration');
xlabel('Frequency (Hz)');
ylabel('LogMag (dB)');
hold on;
plot(sDataCal.S22.freq, sDataCal.S22.logMag);
subplot(2, 2, 2);
plot(sDataCal.S11.freq, sDataCal.S11.phase);
title('Reflection w/Calibration');
xlabel('Frequency (Hz)');
ylabel('Phase (radians)');
hold on;
plot(sDataCal.S22.freq, sDataCal.S22.phase);
subplot(2, 2, 3);
plot(sDataCal.S21.freq, sDataCal.S21.logMag);
title('Transmission w/Calibration');
xlabel('Frequency (Hz)');
ylabel('LogMag(dB)');
hold on;
plot(sDataCal.S12.freq, sDataCal.S12.logMag);
subplot(2, 2, 4);
plot(sDataCal.S21.freq, sDataCal.S21.phase);
title('Transmission w/Calibration');
xlabel('Frequency (Hz)');
ylabel('Phase (radians)');
hold on;
plot(sDataCal.S12.freq, sDataCal.S12.phase);
if nargin == 2
figure
subplot(2, 2, 1);
plot(sDataNoCal.S11.freq, sDataNoCal.S11.logMag);
title('Reflection No Calibration');
xlabel('Frequency (Hz)');
ylabel('LogMag (dB)');
hold on;
plot(sDataNoCal.S22.freq, sDataNoCal.S22.logMag);
subplot(2, 2, 2);
plot(sDataNoCal.S11.freq, sDataNoCal.S11.phase);
title('Reflection No Calibration');
xlabel('Frequency (Hz)');
ylabel('Phase (radians)');
hold on;
plot(sDataNoCal.S22.freq, sDataNoCal.S22.phase);
subplot(2, 2, 3);
plot(sDataNoCal.S21.freq, sDataNoCal.S21.logMag);
title('Transmission No Calibration');
xlabel('Frequency (Hz)');
ylabel('LogMag (dB)');
hold on;
plot(sDataNoCal.S12.freq, sDataNoCal.S12.logMag);
subplot(2, 2, 4);
plot(sDataNoCal.S21.freq, sDataNoCal.S21.phase);
title('Transmission No Calibration');
xlabel('Frequency (Hz)');
ylabel('Phase (radians)');
hold on;
plot(sDataNoCal.S12.freq, sDataNoCal.S12.phase);
hold off;
end
end
Some example plots:
The resonance frequency is the peak of the LogMag transmission spikes (Example above is 7.348 GHz)
This function finds a Q-value of the chamber based on data it collects
function [Q] = getQValue(sData)
[~, xRes] = max(sData.S21.logMag);
resFreq = sData.S21.freq(xRes);
resLogMag = sData.S21.logMag(xRes);
x1 = find(sData.S21.logMag >= (resLogMag - 3), 1, 'first');
x1Freq = sData.S21.freq(x1);
x2 = find(sData.S21.logMag >= (resLogMag - 3), 1, 'last');
x2Freq = sData.S21.freq(x2);
xDif = x2Freq - x1Freq;
Q = resFreq / xDif;
end
This function finds the resonance frequency of a box chamber
function [f] = getFrequency(m, n, l, a, b, c)
%GETFREQUENCY
%m, n, l are modes
%a, b, c are dimensions of cavity
%Assumes MuR * EpsilonR = 1 (vacuum conditions)
f = 1.18e10/2/pi*sqrt((a*pi/m)^2+(b*pi/n)^2+(c*pi/l)^2);
end
*************************************************************************************
3/4
Took more data at higher resonance frequencies
Notice the insane amount of error on one of the traces
This is because something is wrong with the network analyzer
I ran a diagnostics test and it seems to be the port 2 receiver
Looking into buying a new one and ensuring that this is the problem
*************************************************************************************
3/7
Collected more data, assembled into nice powerpoint (will update with slides)
*************************************************************************************
3/8
Started fitting data to Lorentzian function
Wrote multiple functions to do this
fitDataToLorentz
clear fit;
data = sData;
ft = 'lorentzFit(x, w, f, m)';
power = 10 .^(data.S12.logMag(:) ./ 10) ./ 1000;
normPower = power ./ max(power);
x1 = find(normPower >= .5, 1, 'first');
x1Freq = data.S12.freq(x1);
x2 = find(normPower >= .5, 1, 'last');
x2Freq = data.S12.freq(x2);
xDifHalfPower = x2Freq - x1Freq;
%Note that this has to be in ALPHABETICAL ORDER for some reasofn
startPoint = [data.S12.freq(round(length(data.S12.freq)/2)), 1, xDifHalfPower / 2];
[curve, goodness, output] = fit(data.S12.freq(:), normPower(:), ft, 'StartPoint', startPoint);
plot(curve, data.S12.freq, normPower);
lorentzFit
function [y] = lorentzFit(x, w, f, m)
%Where n is number of points, w is half width, and f is resonance frequency
y = zeros(length(x), 1);
for l = 1:length(x)
y(l) = m * w ^ 2 / ((x(l) - f)^2 + w^2);
end
end
*************************************************************************************
3/11
Started fitting to thesis function
Need to learn complex and real analysis with MATLAB
Signed up for lab training for the Spinner
*************************************************************************************
3/14
edited function to fit resonance function - adapted from postdoc Sudeep
function fitstr = fitresonanceEdit(Spower, varargin)
% fitstr = fitresonance(Spower, varargin)
% This fits the transmission coefficient for a resonator. 'Spower' should
% be a three-column matrix [frequency (Hz), linmag, phase (degrees)]
% that corresponds to, for example, S21. The second column is assumed to
% have linear units, so any conversion from dB needs to be done outside of
% this function. 'Spower' can also have two columns, if only the magnitude
% is being fit.
%
% After the required input, optional parameter-value pairs can be
% passed in:
%
% 'fitmode': This is an integer that specifies the fit:
% 1: Fit the magnitude of the power to a simple Lorentzian (default),
% with parameters f0, Q, A (amplitude), using Florentzian.m. This is
% just for a generic resonance. Right now, neither the phase nor
% complex components can be fit.
% 2: Fit the voltage transmission coefficient to Eq. (3.7) in Rangga's
% thesis, with parameters Q, Qep, f0, phi, Cr, Ci using FSresonator.m.
% Most of the parameters are described there, with the exception of
% 'Qep', which is sqrt(Qin Qout)/2, assuming Qin = Qout.
%
% 'fitrange': This specifies the fit range in terms of the initial guess
% for the FWHM. The default is 3.
%
% 'comp': This is a vector of integers to indicate which components to fit
% to. It can contain 1 (magnitude), 2 (phase), 3 (real), 4 (imaginary).
% The default for 'fitmode' 1 is [1] and for 'fitmode' 2 is [3 4]. Right
% now the phase is not unwrapped correctly, so a fit to it probably won't
% work.
%
% 'graph': This is a string that specified the type of plot to make to show
% the fit result. It can be 'plot', 'semilogx', 'semilogy', or 'loglog'.
% If the first letter is a capital, then a new figure is opened first;
% otherwise, the plot is made on the current figure. (This option is for
% calling functions that may generate many plots.) The default is
% 'Plot'. Pass in an empty to skip the graph.
%
% 'residual': If this is true, then the fit residuals are also plotted
% (if graphing is enabled). The default is false. NOT WORKING
%
% 'quiet': If this is true (default), then the updates and warnings from
% LMfit.m are suppressed.
%
% calls: LMfit, Florentzian, Dlorentzian, FSresonator
% created 2/11/16 modified 4/8/16vs
% To do:
% - In fitmode 2, the phase guess is way off; get better HW point, the
% current assumes even df spacing
% - Get Qi and errors as in SACF2
% - Figure out how to get a more reproducible fit range
% - Add residual to fit
% - Figure out how to unwrap phase
% - For the plots, indicate which components have been fit
% - Update callers
% - Add option to estimate noise and refit to get uncertainties? Or can
% same answer come the generic uniform fit?
% Check the inputs with Matlab's parser. At the end, just put the
% optionals in the workspace.
p = inputParser;
p.addRequired('Spower', @(x) isnumeric(x) && ismatrix(x) && ....
ismember(size(x, 2), [2 3 5]));
p.addParameter('fitmode', 1, @(x) ismember(x, [1 2]));
p.addParameter('fitrange', 3, @(x) isnumeric(x) && isscalar(x) && (x > 0));
p.addParameter('comp', [], @(x) isempty(x) || ...
(isnumeric(x) && isvector(x) && all(ismember(x, 1:4))));
p.addParameter('graph', 'Plot', @(x) isempty(x) || ...
ismember(x, {'plot', 'semilogx', 'semilogy', 'loglog', 'Plot', ...
'Semilogx', 'Semilogy', 'Loglog'}));
p.addParameter('residual', false, @(x) islogical(x) || ismember(x, [0 1]));
p.addParameter('quiet', true, @(x) islogial(x) || ismember(x, [0 1]));
p.parse(Spower, varargin{:});
fitmode = p.Results.fitmode;
fitrange = p.Results.fitrange;
comp = p.Results.comp;
graph = p.Results.graph;
residual = p.Results.residual;
quiet = p.Results.quiet;
% Take care of some defaults
if isempty(comp)
switch fitmode
case 1
comp = 1;
case 2
comp = [3 4];
end
end
if (fitmode == 1) && (~isequal(comp, 1))
error('Can only fit to magnitude in ''fitmode'' 1');
end
if ~isempty(graph)
if int8(graph(1)) <= 90
figure;
graph(1) = char(int8(graph(1)) + 32);
end
graph = eval(['@' graph]);
end
% Suppress the warnings for LMfit.m if requested. Does this take time?
if quiet
s1 = warning('off', 'MATLAB:nearlySingularMatrix');
s2 = warning('off', 'MATLAB:singularMatrix');
LMparams = {'Display', 'off'};
else
LMparams = {};
end
% Sort by frequency and get rid of any NaN's.
keep = all(~isnan(Spower), 2);
Spower = Spower(keep, :);
[f, order] = sort(Spower(:,1));
Spower = Spower(order, :);
Plin = Spower(:, 2);
% For the initial guess for the offset, average the lowest 10% on the
% assumption that part of the background is included. If this is not true,
% it shouldn't matter much. For the max, average the top 3 points, in
% case there is some noise. Should make this more robust.
[Psort, order] = sort(Plin);
Poff = mean(Psort( 1 : round(0.1*length(f)) ));
Pmax = mean(Psort(end-2:end));
sub0 = round(mean(order(end-2:end)));
%[Pmax, sub0] = max(Plin);
h0 = Pmax - Poff;
f00 = f(sub0);
df = mean(diff(f));
fwhm0 = sum((Plin - Poff) > h0/2) * df; % rough guess at FWHM
% Get the fit range
fitmax = f00 + fitrange * fwhm0 / 2;
fitmin = f00 - fitrange * fwhm0 / 2;
keep = (f > fitmin) & (f < fitmax);
switch fitmode
case 1
params0 = [Poff, h0, f00, fwhm0/2/sqrt(2)];
params = LMfit(f(keep), Plin(keep), [], params0, [1 1 1 1], ...
@Florentzian, 'Dfit', @Dlorentzian, LMparams{:});
fitstr.f0 = params(3);
fitstr.Q = params(3) / (2 * sqrt(2) * params(4));
fitstr.A = params(2);
fitstr.Ffit = @Florentzian;
fitstr.params = params;
if ~isempty(graph)
hold off
ffull = min(f) : df/3 : max(f);
xl(1) = min(f(keep)) - 0.15 * (max(f(keep)) - min(f(keep)));
xl(1) = max(xl(1), min(f));
xl(2) = max(f(keep)) + 0.15 * (max(f(keep)) - min(f(keep)));
xl(2) = min(xl(2), max(f));
graph(f / 1e9, Plin, '.', 'Color', [0.3 0.3 0.3]); hold on
graph(f(keep) / 1e9, Plin(keep), 'r.');
graph(ffull / 1e9, Florentzian(ffull, params), 'b');
xlim(xl / 1e9);
xlabel('f (GHz)');
ylabel('S');
end
case 2
% Pphase = mod(Spower(:, 3), 360) * pi/180; % radians
Pphase = (mod(Spower(:, 3) + 180, 360) - 180) * pi/180;
VRe = sqrt(Plin) .* cos(Pphase);
VIm = sqrt(Plin) .* sin(Pphase);
Q0 = f00 / fwhm0;
Qin0 = 2 * Q0 / sqrt(Pmax); % Set Qout0 = Qin0 with a NaN
hwsub = sub0 - round(fwhm0/df/2); % phase at half max?
hwsub = min(max(1, hwsub), length(f));
phi0 = Pphase(hwsub);
Cr0 = sqrt(min(abs(Plin)));
Ci0 = Cr0/10;
params0 = [Q0, Qin0, NaN, f00, phi0, Cr0, Ci0];
% Cycle through the components to fit. The first column specifies the
% component. This is probably very inefficient.
x = []; y = [];
xc(:, 2) = f(keep);
for j = 1 : length(comp)
xc(:, 1) = comp(j); x = [x; xc];
switch comp(j)
case 1
y = [y; sqrt(Plin(keep))];
case 2
y = [y; Pphase(keep)];
case 3
y = [y; VRe(keep)];
case 4
y = [y; VIm(keep)];
end
end
% If only the real or imaginary components are fit, don't vary the other
% offset. Don't vary the NaN in Qout.
if isequal(comp, 3)
vary = [1 1 0 1 1 1 0];
elseif isequal(comp, 4)
vary = [1 1 0 1 1 0 1];
else
vary = [1 1 0 1 1 1 1];
end
if isempty(x)
params = NaN(1, 6);
else
params = LMfit(x, y, [], params0, vary, @FSresonator, LMparams{:});
end
fitstr.Q = params(1);
fitstr.Qep = params(2) / 2;
fitstr.f0 = params(4);
fitstr.phi = params(5);
fitstr.Cr = params(6);
fitstr.Ci = params(7);
fitstr.Ffit = @FSresonator;
fitstr.params = params;
% Just plot all of the components for now
if ~isempty(graph)
comps = {sqrt(Plin), Pphase, VRe, VIm};
labels = {'|V|', 'Phase', 'V_R', 'V_I'};
clear xc
xc(:, 2) = min(f) : df/3 : max(f);
xl(1) = min(f(keep)) - 0.15 * (max(f(keep)) - min(f(keep)));
xl(1) = max(xl(1), min(f));
xl(2) = max(f(keep)) + 0.15 * (max(f(keep)) - min(f(keep)));
xl(2) = min(xl(2), max(f));
for j = 1 : 4
subplot(2, 2, j)
hold off
graph(f / 1e9, comps{j}, '.', 'Color', [0.3 0.3 0.3]); hold on
graph(f(keep) / 1e9, comps{j}(keep), 'r.');
xc(:, 1) = j;
graph(xc(:,2) / 1e9, FSresonator(xc, params), 'b');
xlim(xl / 1e9);
xlabel('f (GHz)'); ylabel(labels{j});
end
%Added section to plot imaginary vs. real components
figure;
plot(comps{3}, comps{4}, 'r.');
hold on;
xc2 = xc;
xc(:, 1) = 3;
xc2(:, 1) = 4;
plot(FSresonator(xc, params), FSresonator(xc2, params), 'b');
xlabel(labels{3}); ylabel(labels{4});
axis equal;
%Get Q value
params(1)
end
end
% Restore the warnings to the value at the start. This won't work if this
% function crashes before getting to this point. However, don't want
% to explicitly turn warnings on though, in case they were disabled
% externally...
if quiet
warning(s1); warning(s2);
end
*************************************************************************************
3/25
Collected all current data into a presentation
*************************************************************************************
3/28, 3/29, 4/1
Moving into simulation of the cavity
My measurements will confirm whether the simulation is reliable for our purposes
If so, the cavity can be modeled accurately and a newer one can be engineered
Took a couple of days to learn how the simulator worked, and now I've created a cavity model in the simulation with boundary conditions applied
Tests will be applied on Thursday
Completed cavity simulation
*************************************************************************************
4/4, 4/5, 4/8
After finishing the cavity simulation, I ran tests at each of the different modes
The corresponding Q-values returned by the simulation are similar enough to experimental values that we believe the simulation works
Made lots of graphs
Discovered two new modes, one was (3, 0, 1) which we originally didn't anticipate would be observable, and another is an asymmetric mode that we have no clue existed
We're looking into the asymmetric mode as we might learn something new from it
(1, 0, 1)
(2, 0, 1)
Higher modes - includes (1, 0, 2), (3, 0, 1) and mystery asymmetric mode (note the dip in reflection graphs)
*************************************************************************************
4/11 - 4/12
For these two days, I focused on re-dimensioning the cavity to be exactly aligned with the test cavity we modeled off of. After making these corrections, we adjusted the error on the resonance frequency down to consistently less than 50 MegaHertz, which is incredibly close for our purposes. We also remeasured the Q values, and found that we are now consistently within 2 times the measured amount, which is also close enough for our purposes. This data has given us confirmation that the simulation will work precisely for our intended purpose, which is developing a new cavity that we can more accurately predict the Q-values of. This week I also figured out exactly what the asymmetric mode was by modeling the transformation between the (1, 0,2) mode and the (3,0,1) mode, which showed that the dip in transmission was actually just created by a combination of interference between the two modes creating a strong electric field right in front of the pins. This was really neat to discover, and I have included a picture below to show that strong electric field in front of the pin.
*************************************************************************************
4/15, 4/16, 4/18, 4/19
My primary task this week was to make the cavity superconducting and adjust pin dimensions to study the relationship between these geometries and the returned Q-values. I was able to do a lot of this, and I made some neat progress. On Monday I started by remeasuring the cavity after making it superconducting, which resulted in getting insanely high Q values (as expected) - which were on the order of 30,000 - 60,000 for each mode. The advantage to creating a superconducting cavity is that it ensures that all of the Q will be as a result of the pin dimensions, where as when the cavity is modeled at room temperature the Q is significantly more dependent on the material properties (room temp aluminum). On Tuesday, I began switching around pin dimensions. I started by making the pin diameter larger and smaller - and learned that a larger diameter pin has a much lower Q value than a smaller diameter pin. We also confirmed that the simulation works by using a simple equation - that the inverse of the total Q value at a superconducting temperature is equal to the sum of the inverses of each pin's intrinsic Q value. I confirmed this by measuring Q value with two small pins, two thick pins, and a mix of one thick and one small pin. On Thursday and Friday I went back to my old MATLAB functions and modeled the simulation data to confirm that it was accurate with superconducting measurements - and it was perfect. The graphs were exactly what we wanted to see, and fit the resonance function perfectly. This gives us confidence that I can now model our actual cavities that we use in the refrigerators and ensure that we will get the exactly Q value that we want after I readjust cavity and pin dimensions in the simulation.
*************************************************************************************
4/22, 4/25, 4/26
This week, I've been focusing on modeling our new resonance chamber. Thankfully, this was partially done for me, and I was able to import most of the design straight from Solidworks into ANSYS. The chamber looks like this:
The main difference between this cavity and the previous one is that this one contains a shelf on the outside that attempts to mitigate outside radiation. This is what we are simulating for, as we are trying to see if we can actually get a more controlled Q with this chamber like we predicted.
This chamber simulation will also include a sapphire chip that holds the qubit and several aluminum lines that can be used to communicate with it. For now, these aluminum lines are not very important, as we are looking for the total Q of the chamber and the effects of these lines is minimal. As a result, only one line is on the chamber, and it comes directly through the chamber and out through the shelf.
In these pictures, the sapphire chip is the small blue rectangle, and directly on top of it is the input signal port. The output signal port is coming out of the back of the chamber. In the picture above, you can also see the top-down view of the resonance cavity. This cavity has slightly smaller dimensions, so it returns a significantly smaller first mode resonance frequency, as can be seen in the following plot.
The dip is at approximately 6.3 GHz. This matches exactly what we see in experiments, which is good.
Currently, I am having difficulty getting the input and output ports to work correctly. If I can get this working correctly next week, I should be set to run a lot of tests to look for new Q values, and then I should be able to determine whether this chamber will be more effective than our previous ones.
4/29
I spent a couple of hours trying to get the input pin to work correctly. I now know that the output pin is fine, it's just a matter of ensuring that the input pin has the correct impedance and is not interfering with the aluminum of the cavity.
5/2, 5/3
I have consistently been having issues getting accurate S-parameters for all of the pins. The problem is that they continue to change, which I believe means that the simulation is not converging. As a result, I spent a lot of time trying to get the simulation to do a more refined mesh around the model. Turns out when the model is a lot more complicated the amount of time it takes to mesh increases significantly (definitely not linearly) so I'm looking into ways to localize a very high defined mesh and leave it less defined in unimportant areas. For reference, a mesh is composed of two and three dimensional polygons that the computer solves inside. The smaller the polygons, the more accurate of a simulation will run, but of course it will take more time. I'm still working on getting this mesh to be exactly how I want it.
5/6
Turns out that in order to get the meshing that I wanted all I had to do was increase the amount of runs that the simulation completed. Apparently, the simulation was exceeding the maximum amount of times that it could re-mesh, so the simulation was not converging to a solution, resulting in slightly different measurements every time. I lowered the threshold for convergence (made it more accurate) and increased the cap on maximum simulation reiterations. Even with this, I was getting data that looked quite odd. The cavity that we were modeling based off of had a Q factor of about 300, and we were finding about 2000-3000. As a result, I'm going to go back and ensure that I have modeled everything correctly. There's a chance that this could be correct as well, and perhaps the reason that the cavity was getting such a low Q in real life was due to the qubit chip instead of the resonator itself. It's unclear as of now.
5/9, 5/10
I spent these two days doing two different tasks: continuing to model flux loops and beginning to model qubits. The flux loop is interesting, because it is not behaving as we believed it ought to. For example, shorting the flux loop by attaching it to the superconducting chamber should be giving low Q factors, but instead we’re seeing much higher ones than when the flux loop isn’t shorted. We are unsure why this is, and looking more into theory as to why. Also, we found a huge difference in the quality of our pins and flux loops by shortening them just a couple thousandths of an inch. This is useful to know and we’re containing to look for the trend between pin length and Q factor by taking more data in the simulation. Meanwhile, I’ve also been modeling quints. This is a difficult process, as the qubits are tiny and hard to create on HFSS. It’s also impossible to model quantum effects, so instead we’re introducing them inductively afterwards in Mathematica. I’m still working on getting the qubits to mesh properly in the simulation. Since everything is so small sometimes the simulation doesn’t recognize that parts of the qubit are touching. Modeling these qubits will give us a good way to simulate the relaxation time of our qubit - I.e. how long they can hold their higher energy states. This is the first time my group has tried this and should prove useful if I can get it to work before I leave in two weeks.
*************************************************************************************
Index
AutoCAD
"Hitchhiker's Guide to AutoCAD"
In the process of learning how to design chips
Will update once I begin building them
Fabrication
Finished training for the Nanotechnology Fabrication Laboratory in the Kim Building
Learned how to use the Raith for etching chips
Network Analyzer
Varies signal outputs to range across a wide spectrum of frequencies.
Measures returned phase and amplitude of signals at those frequencies
Has advanced capabilities for measuring and adjusting these output signals, with a total of four ports.
Measures resonance frequencies and responses of devices, such as superconducting qubits in my case
Josephson Junctions
Made out of a sandwich of superconductors with a thin insulator in between.
At cold enough temperatures, electron pairs can move through the insulator with zero resistance and create a super current
This is due to a lack of lattice resistance in the metals and a slight attraction between
HONR268*************************************************************************************
Beginning stuff
Setting up CERN VM
Download virtual box
Download CERN VM
Launch and configure VM
Set up CMS environment - building local release structure, creating a new directory
Create a new analyzer
Ensure that root is configured and that it is possible to grab data from CERN
Some useful commands -
Echo $0 list bash
Echo $SHELL list your shell
Standard format for commands
command -option1 -option2 -option3... target1 target2
Pwd - print working directory
Ls - list files and folders in current folder, options can specify other locations
Cd takes you to new directories
.. goes back a directory
. current
Rmdir - remove directory and everything in it
File-specific commands
ls –CFx /usr/bin
ls –l /bin
find /usr –name “*g++*” –print
which g++
Try typing “ls /home/cms” and then hitting the “tab” key (before hitting return).
Whoami lists user
~ references home directory
There are different file extensions, such as .txt, .x, .exe, and .sh
Mkdir makes directories
Touch will update the time that a certain file or directory was last accessed by “touching” it
HW2
The “touch” command updates the timestamp of last access and modification times to a file to current time. It essentially “touches” the file but doesn’t do anything. If a file name doesn’t exist it’ll make a new file
The ‘.’ is used to explicitly specify that you are searching through your current directory for a file
Putting the * after a string will allow you to find any file name containing that string in the name in the current directory (example above uses the string “1”)
The rm command removes both files and directories, and can be extended to recursively remove all subsequent directories with the -rf command after it. The rmdir command will only work on an empty directory
Since the * wildcard will grab all files with the string “” in the name (which happens to be all of them), this command will delete literally everything in the current directory, which would likely be problematic
The > command will replace everything in a file with the output of the previous command, while the >> command will add the output to the file. They’re similar because you can access files with both of them.
The > command is likely slightly dangerous because it deletes everything in a file when you run it. It’s even more dangerous because it’s easy to accidentally type “>” when you mean “>>”
This command ^ takes the list of everything in the current directory and will replace everything in the log.txt file with it
uname -a is useful for quickly finding the machine’s name and other information about it
cal gives a calendar
cp lets you copy a file to a destination
rev lets you reverse a string (not very useful just kind of funny)
tree lets you print all the file names in a very intuitive way
expr lets you do calculations
history shows history of commands
clear clears the terminal
aspell -c “file name” lets you spell check a file
sed -i ‘s/original/new/g’ file.txt lets you replace words in a text file
This is for text editing
Paragraph on nobel prize:
HW3
This script searches through the usr directory and finds all files with the string argument in the name. Adding the $ also ensures that you are using the value of an argument instead of just the argument name.
It appears as though the chmod command grants permission to all users to use the script.
Running this script
Returns
The output of this file is the testfile itself, because testfile is the only file in the home directory that contains the string “testfile”.
In linux, the piping command | allows you to take outputs from other commands and put them into different commands as inputs. This lets you “layer” different commands into one.
For example:
For the replace script:
Output:
Additional homework problems:
HW4
//hello world script
#include <iostream>
using namespace std;
int main() {
cout <<"Hello World!" << endl; //Print hello world to screen followed by end line (endl)
return 0; //Exit the program
}
g++ compiles the code, and ./a.out executes it
ls /usr/include/c++/4.4.4/
ls /usr/include/math.h
ls /usr/include/bits/*math*.h
These commands will list some of the functions available for use on this VM
Variables intro
#include <iostream>
using namespace std;
int main() {
cout << "hello world" << end;
int i=2;
cout << "i = " <<i<<endl;
double a=3.3;
cout << "a = " <<a<<endl;
int j = a*i;
cout << "a*i = "<<j<<endl;
return 0;
}
Numbers
#include <iostream>
using namespace std;
int main() {
cout << "hello world" << end;
int i=2;
cout << "i = " <<i<<endl;
double a=3.3;
cout << "a = " <<a<<endl;
int j = a*i;
cout << "a*i = "<<j<<endl;
return 0;
}
Numbers 2
#include <iostream>
using namespace std;
int main() {
int n=10;
cout << "n is "<<n<<endl;
n--;
cout<<"n is now "<<n<<endl;
n++;
cout<<n is now "<<n<<endl;
return 0;
}
non-numeric
#include <iostream>
using namespace std;
int main() {
bool prop;
prop = (5>1);
cout<<"prop is "<<prop<<endl;
prop = (1>5);
cout<<"prop is "<<prop<<endl;
prop = (1 != 5);
cout << "prop is " <<prop<<endl;
return 0;
}
Loops 1
#include <iostream>
using namespace std;
int main() {
int n=10;
while(n>0) {
cout<<"n is "<<n<<endl;
n--;
}
return 0;
}
Loops 2
This:
#include <iostream>
using namespace std;
int main() {
int n=0, m=0;
while(n<10) {
// this is the slow (or outer) loop
cout << "n is " << n << ": ";
m=0;
while(m<=n) {
// this is the fast (or inner) loop
// in this loop, the slow loop variable (n) is a constant
// this loop must run to completion before the slow loop
// can progress (during every iteration of the slow loop!)
cout << m;
m++;
}
// now the fast loop has finished and the slow loop can
// continue with the current iteration
cout << endl;
n++;
}
return 0 ;
}
Rewritten with for loops
HW5
Logic statements
commented code
#include <iostream> //required for printing taking in inputs and printing outputs to the command line
using namespace std; //specifies you will be using the C++ standard library, which includes the definitions for some features (such as string)
int main() //sets up the main function to execute code
{
int n = 10; //assign a variable n to have a value 10
while (n>=10) { //only run when n>=10
if(n>5) { //only run if n is greater than 5
cout<<"n is "<<n<<endl; //print out "n is " (value of n)
}
else { //if n is not 5
cout<<"n = "<<n<<endl; //output "n =" (value of n)
}
n--; //subtract 1 from value of n
}
return 0 ; //end program
}
outputs
Pointers
#include <iostream>
using namespace std;
int main() {
int i = 10;
cout << "The memory address of i is " << &i << "\n"; //output memory address of i
cout << "The data stored at memory address " << &i << " is " << i << "\n"; //output actual value of i and its corresponding memory address
int* p = &i; //store a variable p to be the memory address of i
cout << "The value of p is " << p << "\n"; //output memory address of i
cout << "We say that p 'points at' the memory location referenced by address " << p << "\n"; //also output memory address of i, but this time label it
cout << "The data stored at memory address " << p << " is " << *p << "\n"; //output the value of p (which is memory address of i) and then the actual value of i, which is 10
return 0;
}
output
nice
Program 1
/*PROGRAM 1*/
#include <iostream>
using namespace std;
int main(){
int i = 10;
int j = i;
cout << "i= " << i << " and j= " << j << "\n"; //output the values of j and I, which are each 10
i=5; //redefine i as 5
cout << "i= " << i << " and j= " << j << "\n"; //once again output the values of i and j
j=1; //redefine j as 1
cout << "i= " << i << " and j= " << j << "\n"; //output new values of i and j
return 0; //end program
}
program 2
/*PROGRAM 2*/
#include <iostream>
using namespace std;
int main(){
int i = 10; //define i = 10
int *p = &i; //define the memory address of p to be the memory address of i
cout << "i= " << i << " and *p= " << *p << "\n"; //output the value of i and the actual value of p, which is i
i=5; //redefine i
cout << "i= " << i << " and *p= " << *p << "\n"; //output i and then the. value of p, which happens to be the old value of i
*p=1; //redefine the value that p points to
cout << "i= " << i << " and *p= " << *p << "\n"; //output new value of I and the new pointed value to from p
return 0; //end the program
}
Program 2 redefines
Example with "new"
#include <iostream>
using namespace std;
int main(){
int* p = new int(5); //define p as an address that points to the value of 5
cout << "p points at address " << p << "\n"; //output p
cout << "The data stored in address " << p << " is " << *p << "\n"; //output data in p
*p = 10; //redefine the actual value that is pointed to by p
cout << "Now the data stored in address " << p << " is " << *p << "\n"; //output the new value pointed at by p
return 0;
}
output
HOMEWORK
Physics
1) The projected range for copper is 4.739 * 10^2 g/cm^2 with a 1 GeV kinetic energy. The density of copper is 8.96 g/c^3, so multiplying projected range by 1/density of copper will return the length of copper necessary in cm to stop a proton with 1 GeV of kinetic energy. This numbers is 4.739*10^2 / 8.96 = 52.89 cm of copper.
2)
Computing
1)
#include <iostream>
using namespace std;
int main()
{
int i = 0;
cout << "This is the value of i " << i << endl;
int j = 20;
cout << "This is the value of j " << j << endl;
while (i < j)
{
int dif = j - i;
cout << "i is now " << dif << " away from j" << endl;
if (i % 2 == 0)
{
cout << "i is even" << endl;
}
else
{
cout << "i is odd" <<endl;
}
i++;
}
return 0;
}
Output
2)
#!/bin/tcsh
set month = `date | awk '{print $2}'`
echo $month
ls -lat -R ../| grep "$month"
Outputs
HW6
Setting up arrays
#include <iostream>
using namespace std;
int main() {
int ii[3] = {1,2,3}; //set up array
int j=0; //define j to be 0
while (j<3) { //create while loop
cout <<" ii of "<<j<<" is "<<ii[j]<<endl; //output all values of the array ii
j++; //do this by incrementing j
}
int LL[2][3] = {1,2,3,4,5,6}; //create a new matrix of six values with 2 rows and 3 columns
j=0; //define j
int k; //define integer variable k
while (j<2) { //create a while loop to increment j, restart k value every loop, incrementing the rows
k=0;
while (k<3) { //create another while loop to increment k (columns) and end up printing out all values of the matrix LL
cout<<" LL of "<<j<<" "<<k<<" is "<<LL[j][k]<<endl;
k++;
}
j++;
}
return 0; //end
}
#include <iostream>
using namespace std;
/************************************************\
* *
* Arrays *
* This program demonstrates arrays *
* *
\************************************************/
int main() {
// a loop to demonstrate 1D arrays
int ii[3] = {1,2,3};
int j=0;
while (j<3) {
cout <<" ii of "<<j<<" is "<<ii[j]<<endl;
j++;
}
// a loop ot demonstrate 2D arrays
int LL[2][3] = {1,2,3,4,5,6};
j=0;
int k;
while (j<2) {
k=0; // do not forget to initialize k here
while (k<3) {
cout<<" LL of "<<j<<" "<<k<<" is "<<LL[j][k]<<endl;
k++;
}
j++;
}
return 0;
}
Same script, just with helpful comments along the way
Examples with files
#include <iostream>
#include <fstream> //this allows you to execute commands involving files
using namespace std;
int main() {
ofstream myfile; //declare an object to be used for writing files
myfile.open("example.txt"); //declare a file name and open said file
myfile<<"write some junk."; //put words in the file
myfile.close(); //close the file
return 0; //end program
}
Making two files
Main
/* MAIN.CPP */
#include <iostream>
#include <fstream>
// include the program dotprod.cpp so that we can find the dot_prod function
#include "dotprod.cpp"
using namespace std;
int main () {
// declare the vectors
double vector1[3]; //create an array of doubles
double vector2[3]; //another one
// open the input file
ifstream infile; //open a file to take input
infile.open("vectors.txt"); //name it "vectors.txt"
// store the values from the vector file and print the vectors for the user
infile>>vector1[0]>>vector1[1]>>vector1[2]; //put intputs into vector file
cout<<" Vector 1 is ("<<vector1[0]<<","<<vector1[1]<<","<<vector1[2]<<")"<<endl; //list values of each respective inputted vector
infile>>vector2[0]>>vector2[1]>>vector2[2];
cout<<" Vector 2 is ("<<vector2[0]<<","<<vector2[1]<<","<<vector2[2]<<")"<<endl;
// close the input file
infile.close();
// call the dot_prod function from dotprod.cpp
dot_prod(vector1,vector2);
return 0;
}
Dotprod
/* DOTPROD.CPP */
#include <iostream>
using namespace std;
double dot_prod(double v1[3],double v2[3]) { //define a function that takes in v1 and v2 vectors as inputs
double dotdot; //define value for dot product answer
dotdot = v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]; //compute dot product
cout<<" The dot product is "<<dotdot<<endl; //output dot product
return 0; //end function
}
Modified code to do scalar operations -
Main
/* MAIN.CPP */
#include <iostream>
#include <fstream>
// include the program dotprod.cpp so that we can find the dot_prod function
#include "dotprod.cpp"
#include "scalar_mult.cpp"
using namespace std;
int main () {
// declare the vectors
double vector1[3]; //create an array of doubles
double vector2[3]; //another one
double scalar; //create scalar variable
// open the input file
ifstream infile; //open a file to take input
infile.open("vectors.txt"); //take from name "vectors.txt"
// store the input in the vectors and print the vectors for the user
infile>>vector1[0]>>vector1[1]>>vector1[2]; //put inputs into vector file
cout<<" Vector 1 is ("<<vector1[0]<<","<<vector1[1]<<","<<vector1[2]<<")"<<endl; //list names of each respective inputted vector
infile>>vector2[0]>>vector2[1]>>vector2[2];
cout<<" Vector 2 is ("<<vector2[0]<<","<<vector2[1]<<","<<vector2[2]<<")"<<endl;
infile>>scalar;
cout <<"The scalar is " <<scalar<<endl;
// close the input file
infile.close();
// call the dot_prod function from dotprod.cpp
dot_prod(vector1,vector2);
scalar_mult(vector1, scalar);
scalar_mult(vector2, scalar);
return 0;
}
Scalar_mult
/*SCALAR_MULT.CPP*/
#include <iostream>
using namespace std;
double scalar_mult(double v1[3], double scalar) {
double scaled_vector[3];
scaled_vector[0]=v1[0]*scalar;
scaled_vector[1]=v1[1]*scalar;
scaled_vector[2]=v1[2]*scalar;
cout<<" The scaled vector is ("<<scaled_vector[0]<<","<<scaled_vector[1]<<","<<scaled_vector[2]<<")"<<endl;
return 0;
}
Output
Random numbers
#include <iostream>
#include <math.h>
using namespace std;
// this function is the actual random number generator
// this code is stolen from the book numerical recipes for fortran
// it relies on random generated from overflow of memory locations
// and is a pseudo random number generator
const int a = 7141;
const int c = 54773;
const int mmod=256200;
double getFlatRandom(int& inew) { //define a new function to take in the pointer of a integer variable
cout << inew <<endl;
double mranflat = 0.;
inew = inew%mmod;
double aa = double(inew)/double(mmod);
mranflat=aa;
inew = a*inew+c;
return mranflat;
}
// in this code, we will call the pseudo-random number generator and learn some things
// about its properties by filling and then displaying a histogram
int main() {
int num;
cout << "Enter the number of loop iterations: ";
cin >> num;
int inew = 2345; // This is the "seed" for the random number generator
// we will put the results from the call into a histogram that we can look at, to learn some of its
// properties. This histogram has 10 bins.
int histo[10] = {0,0,0,0,0,0,0,0,0,0};
double atmp; //
// call the random number generator 1000 times and fill a histogram
for (int i = 0; i<num; i++){
atmp=getFlatRandom(inew); // call the random number generator
histo[int(atmp*10)]++; // increment the histogram bin the number falls within
}
// print the histogram to the screen
for (int i = 0; i<10; i++){
cout<<i<<": ";
for (int j=0; j<int((double(100*histo[i])/double(num))+0.5); j++) {
cout << "=";
}
cout << endl;
}
return 0;
}
Note that this is seeded, so the generated numbers are the same as long as the seed remains the same. Also note that this generator does not work when the memory access allower (&) is removed
Demonstrating cout and how distribution evens out over time:
N=10
N=100
N=1000
Homework part
Code
#include <iostream>
#include <math.h>
using namespace std;
// this function is the actual random number generator
// this code is stolen from the book numerical recipes for fortran
// it relies on random generated from overflow of memory locations
// and is a pseudo random number generator
const int a = 7141;
const int c = 54773;
const int mmod=256200;
double getFlatRandom(int& inew) { //define a new function to take in the pointer of a integer variable
cout << inew <<endl;
double mranflat = 0.;
inew = inew%mmod;
double aa = double(inew)/double(mmod);
mranflat=aa;
inew = a*inew+c;
return mranflat;
}
// in this code, we will call the pseudo-random number generator and learn some things
// about its properties by filling and then displaying a histogram
int main() {
int num;
cout << "Enter the number of loop iterations: ";
cin >> num;
int inew = 2345; // This is the "seed" for the random number generator
// we will put the results from the call into a histogram that we can look at, to learn some of its
// properties. This histogram has 10 bins.
int histo[10] = {0,0,0,0,0,0,0,0,0,0};
double atmp; //
// call the random number generator 1000 times and fill a histogram
for (int i = 0; i<num; i++){
atmp=getFlatRandom(inew); // call the random number generator
histo[int(atmp*10)]++; // increment the histogram bin the number falls within
}
// print the histogram to the screen
for (int i = 0; i<10; i++){
cout<<i<<": ";
for (int j=0; j<int((double(100*histo[i])/double(num))+0.5); j++) {
cout << "=";
}
cout << endl;
}
return 0;
}
Removing the & results in code that returns only 0’s. The & allows you to redefine the variable inside the function every time the loop runs
Each byte contains 8 bits. Since each bit has two combinations, each byte has 2^8 combinations. Then when combining multiple bytes you get 2^8 * 2^n -1 = 2^(8*n) -1. (The minus 1 is because the value of binary numbers starts at 0, so it’s shifted by one).
The numbers are the same each loop because the starting seed is remaining the same and every number is being “generated” by the same formula with respect to that seed. Changing the seed results in different numbers being generated
As you increase the number of loops being run, the histogram levels out more and more
Calorimetry
With N=1, I would guess that the mass of the mother particle is approximately 111 eV. It’s hard to tell how far away the true mass would be, because with only one data point there is not a lot of consensus on what the average value actually is.
Changing N to 10 yields this plot, which at least starts to give an idea of what the middle of the data is
Changing N to 100 yields this plot, which begins to show a Gaussian distribution
Changing N to 1000 gives a clear picture of the actual mass, which appears to be around 90.
After changing the true mass to 1 and adjusting the amount of points taken, it’s clear that the distribution is no longer Gaussian, rather it is skewed to the right. Oddly enough, it’s not even centered at 1, and instead sits around 10. This is because mass cannot be negative, so the distribution cannot remain Gaussian about 1.
Picture above with N=10000.
Energy
HW7
CODING PART
The above code runs madgraph and generates some collisions
Below specifies the output destination
Launching this event begins the simulation
You can then view results
Lots of data is available after unzipping the data file and converting it to a root file using a python script
HW PART
Beam pipe - Photons travel through and collide in the beam pipe.
Tracker - The tracker is used to determine the paths of particles after they initially collide. About 1 meter radius
Calorimeter - The calorimeter measures energy of particles as they pass through it. About a 2.5 meter radius
Solenoid - The solenoid generates a uniform 4 Tesla magnetic field. It has a 3.5 meter radius
Muon tracker - Muons don't get stopped by the calorimeter, so the extra muon chambers are used to measure the energies of muons. Has a radius of 7 meters.
Display - Press the button that focuses the z display coming out of the screen to look along the beam pipe
Electrons - Bright green lines
Muons - Red lines
Missing energy - Dashed blue lines
Jets - Jets are represented as green cones
Direction - Since the magnetic field is going into the page, the particle will be negatively charged if it has a clockwise curvature.
4Lepton
Run 1 - Run / Event / LS: 178424 / 666626491 / 585
Electrons - 0
Muons - 12
Photons - 0
Jets - 54
Missing Energy - 2
Vertices - 6
Run 2 - Run / Event / LS: 193575 / 400912970 / 523
Electrons - 4
Muons - 3
Photons - 4
Jets - 106
Missing Energy - 2
Vertices - 11
Diphotons
Run 3 - Run / Event / LS: 199319 / 641436592 / 464
Electrons - 0
Muons - 2
Photons - 2
Jets - 31
Missing Energy - 1
Vertices - 11
Run 4 - Run / Event / LS: 199699 / 336259924 / 272
Electrons - 0
Muons - 0
Photons - 2
Jets - 29
Missing Energy - 1
Vertices - 6
Electrons
Run 5 - Events/Run_146644/Event_718940510 [1 of 25
Electrons - 1
Muons - 0
Photons - 2
Jets - 10
Missing Energy - 1
Vertices - 1
Run 6 - Events/Run_146644/Event_718948926 [2 of 25]
Electrons - 1
Muons - 0
Photons - 1
Jets - 18
Missing Energy - 1
Vertices - 4
Muon
Run 7 - Events/Run_146436/Event_90625265 [1 of 25]
Electrons - 0
Muons - 2
Photons - 0
Jets - 0
Missing Energy - 1
Vertices - 1
Run 8 - Events/Run_146436/Event_90626440 [2 of 25]
Electrons - 3
Muons - 2
Photons - 0
Jets - 17
Missing Energy - 0
Vertices - 0
MinimumBias
Run 9 - Events/Run_148031/Event_441572355 [1 of 25]
Electrons - 0
Muons - 0
Photons - 0
Jets - 16
Missing Energy - 1
Vertices - 1
Run 10 - Events/Run_148031/Event_441584555 [2 of 25]
Electrons - 0
Muons - 0
Photons - 0
Jets - 0
Missing Energy - 1
Vertices - 1
HW8
Adding 3rd and 4th vectors and setting parameters
Add them to get zCandidate2
Add zCandidate and zCandidate2 to get Higgs
Create and write to the Higgs histogram
Output of histograms
HW9
Here’s a link to colab if that would be easier to see https://colab.research.google.com/drive/1kuZw4_0xpEusxc0AE9aMtB21kX5BS-c_
HW10
Deep Dream code
# -*- coding: utf-8 -*- """deep_dream.ipynb Automatically generated by Colaboratory. Original file is located at https://colab.research.google.com/drive/12P9vBkS58sj6JKVJ3Xf3reJGMQuIvxoC """ from __future__ import print_function from keras.preprocessing.image import load_img, save_img, img_to_array import numpy as np import scipy import argparse from keras.applications import inception_v3 from keras import backend as K parser = argparse.ArgumentParser(description='Deep Dreams with Keras.') parser.add_argument('base_image_path', metavar='base', type=str, help='Path to the image to transform.') parser.add_argument('result_prefix', metavar='res_prefix', type=str, help='Prefix for the saved results.') args = parser.parse_args() base_image_path = args.base_image_path result_prefix = args.result_prefix # These are the names of the layers # for which we try to maximize activation, # as well as their weight in the final loss # we try to maximize. # You can tweak these setting to obtain new visual effects. settings = { 'features': { 'mixed2': 0.2, 'mixed3': 0.5, 'mixed4': 2., 'mixed5': 1.5, }, } def preprocess_image(image_path): # Util function to open, resize and format pictures # into appropriate tensors. img = load_img(image_path) img = img_to_array(img) img = np.expand_dims(img, axis=0) img = inception_v3.preprocess_input(img) return img def deprocess_image(x): # Util function to convert a tensor into a valid image. if K.image_data_format() == 'channels_first': x = x.reshape((3, x.shape[2], x.shape[3])) x = x.transpose((1, 2, 0)) else: x = x.reshape((x.shape[1], x.shape[2], 3)) x /= 2. x += 0.5 x *= 255. x = np.clip(x, 0, 255).astype('uint8') return x K.set_learning_phase(0) # Build the InceptionV3 network with our placeholder. # The model will be loaded with pre-trained ImageNet weights. model = inception_v3.InceptionV3(weights='imagenet', include_top=False) dream = model.input print('Model loaded.') # Get the symbolic outputs of each "key" layer (we gave them unique names). layer_dict = dict([(layer.name, layer) for layer in model.layers]) # Define the loss. loss = K.variable(0.) for layer_name in settings['features']: # Add the L2 norm of the features of a layer to the loss. if layer_name not in layer_dict: raise ValueError('Layer ' + layer_name + ' not found in model.') coeff = settings['features'][layer_name] x = layer_dict[layer_name].output # We avoid border artifacts by only involving non-border pixels in the loss. scaling = K.prod(K.cast(K.shape(x), 'float32')) if K.image_data_format() == 'channels_first': loss += coeff * K.sum(K.square(x[:, :, 2: -2, 2: -2])) / scaling else: loss += coeff * K.sum(K.square(x[:, 2: -2, 2: -2, :])) / scaling # Compute the gradients of the dream wrt the loss. grads = K.gradients(loss, dream)[0] # Normalize gradients. grads /= K.maximum(K.mean(K.abs(grads)), K.epsilon()) # Set up function to retrieve the value # of the loss and gradients given an input image. outputs = [loss, grads] fetch_loss_and_grads = K.function([dream], outputs) def eval_loss_and_grads(x): outs = fetch_loss_and_grads([x]) loss_value = outs[0] grad_values = outs[1] return loss_value, grad_values def resize_img(img, size): img = np.copy(img) if K.image_data_format() == 'channels_first': factors = (1, 1, float(size[0]) / img.shape[2], float(size[1]) / img.shape[3]) else: factors = (1, float(size[0]) / img.shape[1], float(size[1]) / img.shape[2], 1) return scipy.ndimage.zoom(img, factors, order=1) def gradient_ascent(x, iterations, step, max_loss=None): for i in range(iterations): loss_value, grad_values = eval_loss_and_grads(x) if max_loss is not None and loss_value > max_loss: break print('..Loss value at', i, ':', loss_value) x += step * grad_values return x """Process: - Load the original image. - Define a number of processing scales (i.e. image shapes), from smallest to largest. - Resize the original image to the smallest scale. - For every scale, starting with the smallest (i.e. current one): - Run gradient ascent - Upscale image to the next scale - Reinject the detail that was lost at upscaling time - Stop when we are back to the original size. To obtain the detail lost during upscaling, we simply take the original image, shrink it down, upscale it, and compare the result to the (resized) original image. """ # Playing with these hyperparameters will also allow you to achieve new effects step = 0.01 # Gradient ascent step size num_octave = 3 # Number of scales at which to run gradient ascent octave_scale = 1.4 # Size ratio between scales iterations = 20 # Number of ascent steps per scale max_loss = 10. img = preprocess_image(base_image_path) if K.image_data_format() == 'channels_first': original_shape = img.shape[2:] else: original_shape = img.shape[1:3] successive_shapes = [original_shape] for i in range(1, num_octave): shape = tuple([int(dim / (octave_scale ** i)) for dim in original_shape]) successive_shapes.append(shape) successive_shapes = successive_shapes[::-1] original_img = np.copy(img) shrunk_original_img = resize_img(img, successive_shapes[0]) for shape in successive_shapes: print('Processing image shape', shape) img = resize_img(img, shape) img = gradient_ascent(img, iterations=iterations, step=step, max_loss=max_loss) upscaled_shrunk_original_img = resize_img(shrunk_original_img, shape) same_size_original = resize_img(original_img, shape) lost_detail = same_size_original - upscaled_shrunk_original_img img += lost_detail shrunk_original_img = resize_img(original_img, shape) save_img(result_prefix + '.png', deprocess_image(np.copy(img)))
Example filtered image