Dylan Poppert

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

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.

Output:

Additional homework problems:

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

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

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

#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

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.

• • Picture above with N=10000.

Energy

HW7

CODING PART

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

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