Invited talk

 We were invited to a session of the 181th meeting of ASA.

https://www.researchgate.net/publication/356367868_High-performance_computing_for_long-range_underwater_acoustics

Guest Investigator @Woods Hole Oceanographic Institution

 From this December, I am also a Guest Investigator at Woods Hole Oceanographic Institution.

I will work on HPC in underwater acoustics.

Hope I can provide some good achievements and make any feed-back to CTBT.

heFFTe

 Installation of heFFTe with MKL on SPACK;

% spack install heffte%intel@2021.4.0+mkl ^intel-oneapi-mpi@2021.4.0 ^intel-mkl threads=openmp

To use installed Intel OneAPI and CMAKE on SPACK

In order to avoid multiple installations of Intel-OneAPI and CMAKE, we can use installed ones.
Please install Install OneAPI and CMAKE via APT and put the packages.yaml under ~/.spack/ whose contents are as follows;

 $ cat ~/.spack/packages.yaml

packages:

intel-oneapi-mkl:

buildable: false

externals:

- spec: intel-oneapi-mkl@2021.4.0

prefix: /opt/intel/oneapi/

intel-oneapi-mpi:

buildable: false

externals:

- spec: intel-oneapi-mpi@2021.4.0

prefix: /opt/intel/oneapi/

cmake:

buildable: false

externals:

- spec: cmake@3.16.3

prefix: /usr

buildable: false prevents us from installing packages via SPACK.

MPI/OpenMP hybrid with Intel MPI

  mpirun -np 8 -genv OMP_NUM_THREADS=16 -hosts node2-ib,node3-ib,node4-ib,node5-ib -ppn 2 ./sol

To run an OpenMP/MPI hybrid code with IntelMPI, the above command is needed.

The meanings are;

8 MPI processes (-np 8),

16 OpenMP threads (-genv OMP_NUM_THREADS=16) ,

2 MPI processes per node (-ppn 2).

Reverse vertical (y) axis of Google spreadsheet graph

Although you don't expect it, plotting the below graph is a bit tricky on Google Spreadsheet, since the vertical axis goes from the largest to smallest as the bottom to top.

To do this, we need to change the "appearance" of the values in a column.

This can be achieved by applying, format -> number -> more formats -> custom formats on Google Spreadsheet, after you select a column of your spreadsheet.

 For the above plot, I used "-0;0".

Trip to Gumunden and Attersee

 ちょっとした夏休みの旅行ということで、2泊3日でグムンデンとアッター湖に行ってきました。
グムンデンは妻が食器を見たいということで、Gmundner Keramic を目標に、アッター湖はクリムトに興味があるとのことで Gustav Klimt Zentrumに。

食事は、初日の夕方はイワナを食べに、トラウン湖の西湖岸のGegrillte Steckerlfische Trawöger-Dorfnerにいき（食事は魚しか選択肢がない）、

二日目は私の要望で、Brot und Weinで（ハウスワインが超おすすめ。食事は日替わり3種類）。
GmundenからGustav Klimt Zentumは電車だと遠いのでタクシーを使いました。20分弱で40ユーロ。

(www.taxigmunden.at)

宿はGumundenにある、Ferienwohnun Gailer。AirBnBみたいな感じで、一部屋をそのまま借りました。8月5日はGumndenのお祭り？(https://www.mondscheinbummel.at/)で、ちょっと盛り上がっておりました。

トラウン湖もアッター湖も非常に水が綺麗でいい意味で裏切られました。もう一度行くなら泳いでもいいかもしれない。

トラウン湖

イワナ

 Brot und Wein

Klimt Zentrumの庭

宿近くの橋からトラウン川とKasbergを望む

Intall Octave with amdblis and amdlibflame with Spack

 spack install octave+hdf5+fftw ^hdf5~mpi ^fftw~mpi ^amdblis threads=openmp ^amdlibflame threads=openmp

To enable threading you need to set OMP_NUN_THREADS=4 (e.g. with 4 cores).

Using pyplot@Plots on Julia@spack

 Using pyplot@Plots on Julia installed via spack is a bit complicated.

We first need to setup Python environment
(actually, you may want to use conda provided by Julia, so, you may do
Julia> ENV["PYTHON"]=""
)

as described in the following link;
https://github.com/JuliaPy/PyCall.jl

and then,

Julia> using Plots
Julia> pyplot()

Phase estimation

 Since we have the inverse QFFT, we can move to the phase estimation.

In the example, we are going to estimate the eigenvalues of the following matrix:

.

This matrix (or operator) can be written as a tensor product of 

and

.

.

This matrix is obviously unitary and can be used as the U operator in the diagram.
This matrix must be used with a control qubit, and I created a new "operation" matrixOperation in the following code.
The results agree with the example, and it seems that no need to repeat the measurement many times. The following table shows the number of measured bit strings out of 1,000 trials.

BTW, MS-Q# has QFT function which is "Adjoint"able (invertible), and the results are completely identical to mine. And also, I realize that the R1 gate should also take Adjoint, and I should use it instead of multiplying -1.0 to theta.

namespace phaseEstimation {

    open Microsoft.Quantum.Canon;

    open Microsoft.Quantum.Intrinsic;

    open Microsoft.Quantum.Convert;

    open Microsoft.Quantum.Arrays;

    open Microsoft.Quantum.Measurement;

    open Microsoft.Quantum.Math;

    open Microsoft.Quantum.Arithmetic as Arithmetic;

    

    @EntryPoint()

    operation SayHello() : Unit {

        let nQubit = 3; // Number of control qubits

        let maxNum = PowI(2,nQubit);

        mutable counts = new Int[maxNum];

        mutable theta  = new Double[maxNum];

        mutable revSeq = new Int[nQubit];

        for iQubit in 0..(nQubit-1){

            set theta w/= iQubit <- calcTheta(IntAsDouble(iQubit+1));

        }

        for iQubit in 0..(nQubit-1){

            set revSeq w/= iQubit <- nQubit-1-iQubit;

        }

        for trial in 1..1000 {

        use  (resultBits,vec) = (Qubit[nQubit],Qubit[2]) {

            ApplyToEach(Reset,resultBits); // Zero-clear the control qubits

            ApplyToEach(H,resultBits);

            //X(vec[0]);

            //X(vec[1]);

            for iQubit in 0..(nQubit-1){

                (Controlled matrixOperation)([resultBits[iQubit]],(iQubit,vec));

            }

            // Start Inverse QFFT

            let numLoopSwap = Floor(IntAsDouble(nQubit)*0.5);

            for iQubit in revSeq{

                for jQubit in 0..(nQubit-1)-iQubit-1{

                    let kQubit = revSeq[(nQubit-1)-jQubit];

                    let controlQubit = nQubit-kQubit-1;

                    let rotationIndex= nQubit-(kQubit+iQubit)-1;

                    (Controlled R1)([resultBits[controlQubit]],(-1.0*theta[rotationIndex],resultBits[iQubit]));

                }

                H(resultBits[iQubit]);

            }

            //for iQubit in 0..(numLoopSwap-1) {

            //    SWAP(resultBits[iQubit],resultBits[nQubit-1-iQubit]);

            //}              

            //let register = Arithmetic.BigEndian(resultBits);

            //(Adjoint QFT)(register);         

            let res = MultiM(resultBits);

            let bit = ResultArrayAsInt(res);

            set counts w/= bit <- counts[bit]+ 1;

            ResetAll(resultBits);

            ResetAll(vec);

        }

        }

        for j in 0..(maxNum-1)  {

           Message(IntAsString(counts[j]));

        }

    }

    function calcTheta(l : Double) : Double {

        return 2.0*PI()/(PowD(2.0,l));

    }

    operation matrixOperation(iQubit : Int, q : Qubit[]): Unit is Ctl{

       //Message(IntAsString(iQubit)+" "+IntAsString(PowI(2,iQubit)));

        for i in 0..(PowI(2,iQubit)-1){

            S(q[0]);

            T(q[1]);

        }

        return ();

    }

}

コロナウイルス

流石にちょっとぐらい理解しておいたほうが良いかなと思ってネットを巡っていたら良い資料が見つかったので。 

新型コロナウイルスのウイルス学的特徴

https://www.eiken.co.jp/uploads/modern_media/literature/P11-18.pdf

Compile spack packages for generic X86_64

 It seems that Spack automatically generates very optimized binaries (probably with something like gcc -march=native).

This is usually a good thing, but if a cluster consists of several generations of CPUs, this would cause a problem. In my case, the front node is newer than compute nodes, and the spack modules die on compute nodes.

To avoid this, you may want to add the following in ./etc/spack/defaults/packages.yaml.

Then the compiled binaries are for generic x86_64.

packages:
  all:
    target: [x86_64]

Appended on 30/July/2021

It should be better to have a user level configuration rather than changing the global file 
(the global config file will be overwritten with git pull or you will be complained).

$ cat ~/.spack/packages.yaml
packages:
  all:
    target: [x86_64]

https://spack.readthedocs.io/en/latest/configuration.html

Profiling Oct-file

Since Oct file of Octave is used to get a speedup, I think profiling of oct-files is useful, but I have not been able to do so. Oct-file is an executable binary called by Octave, and so, the profilers might be confused.

Recently, I have just found that "Intel Advisor" can get the information of my oct file. In addition, it really gives me the advice to get better performance.

We just need to launch Octave script to get the performance, namely you need to specify the following on Intel Advisor,
Binary: octave or octave-cli,
Option: an Octave script that calls your Oct-file.

That's it!

Compile Octave with intel compiler

I finally found a way to build Octave with intel compiler and mkl.

$ cat conf_ifort.sh

CPPFLAGS=' -I${MKLROOT}/include -I${MKLROOT}/include/fftw' \

LDFLAGS=' -L${MKLROOT}/lib/intel64' \

CFLAGS="-O3 -fPIC -DM_PI=3.1415926535897932384  -I${MKLROOT}/include -I${MKLROOT}/include/fftw" \

CC=icx CXX=icpx F77="ifort" \

../octave-6.2.0/configure --prefix=/path/to/OCTAVE/6.2.0 \

--with-blas="-lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm -ldl" \

--with-lapack="-lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm -ldl" \

--with-fftw3="-lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm -ldl" \

--with-fftw3f="-lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm -ldl" \

--disable-java --enable-static --enable-fortran-calling-convention=gfortran 

bad memory

 I sometimes recall the memory that my former supervisor rejected my proposal (GPU computation), but my stolen idea became vital with one of my colleagues after half a year.

OK, if I did not have any proven record while he had, I should have been convinced, but the truth is the other way round.

Or, some disgusting my colleagues contributed to disturbing my scientific activities.

Even UN is not a fair place.

Illegal instruction by an empty locale????

 I have just got access to a hand-made PC cluster, and I encounter an interesting problem.

A compute node complains of "illegal instruction", while the architecture is almost the same as that of the head node (Sandybridge, and Haswell). Furthermore, I experienced the problem even with binaries compiled on the compute node. That is totally clueless.

When I login (yeah, it is allowed, luckily) to a compute node, I realize that the terminal says

"-bash: warning: setlocale: LC_CTYPE: cannot change locale (ja_JP.utf8): No such file or directory".

I did the following half in doubt, and the "illegal instruction" disappears!;
"export LC_ALL=en_US.utf8"

I do not know what is really going on, but this is my story.

offsetarray as a member of structure in Julia

 I do feel that accessing offsetarray, which is a member of a structure in Julia, is extremely slow.

In my case, I use a staggered grid in a FDTD code, and defined a structure as;

struct grid

    p :: OffsetArray{Float64,2}

    vx:: OffsetArray{Float64,2}

    vy:: OffsetArray{Float64,2}

end

then accessed like (output, input, resOfCalcF are all "grid".)

for j in -m+1:gridParam.ny+m

    for i in -m+1:gridParam.nx+m

         output.p[i,j] = input.p[i,j] + coef*resOfCalcF.p[i,j]

    end

end

It seems that Julia cannot determine the type of ".p" and this causes the type-instability.

To avoid this, I needed to cast type as 

op = parent(output.p)::Matrix{Float64}

ip = parent(input.p)::Matrix{Float64}

rp = parent(resOfCalcF.p)::Matrix{Float64}

and the array starts from 1 (the loop starts from 2, because the outer most layer is just a dummy)

for j in 1+1:gridParam.ny+m+m+1

    for i in 1+1:gridParam.nx+m+m+1

         op[i,j] = ip[i,j] + coef*rp[i,j] :: Float64

    end

end

I have not measured the speed-up, but probably x10 - x100 speedup is obtained.

ジム再開/Re-opening of gyms in Vienna

 5月19日から3G(陰性証明書/ワクチン接種証明書/罹患証明書の略らしい)を持っていればジムに入れるようになりました。おかげで首の調子もだいぶ戻ったし、血圧も一気に下がった。これから徐々に制限が解除されていくようだけれど、3G所持じゃないとダメな行動が広がるということでもあり、早めにワクチン接種したい。

From the 19th of May, the gyms in Vienna were allowed to re-open, but for those who have 3G certificates (certificates of negative/vaccinated/recovery). Because of the lack of exercise, my neck started complaining to me, and I got hypertension. after a couple of weeks, they are getting better. importance of exercise :)

Getting the negative certificate is not difficult whilst a bit annoying. I may wish to get vaccinated soon.

Octave 6.2.0 + MKL (intel-oneapi-mkl 2021)

The following is the configure option to include MKL in Octave.

MKL is installed via Spack, and the version is intel-oneapi-mkl@2021.2.0

CPPFLAGS=" -I${MKLROOT}/include -I${MKLROOT}/include/fftw" \

LDFLAGS=" -L${MKLROOT}/lib/intel64" \

CFLAGS="-O3 -fPIC -DM_PI=3.1415926535897932384  -I${MKLROOT}/include -I${MKLROOT}/include/fftw" \

F77="gfortran" \

../octave-6.2.0/configure --prefix=/path/to/install \

--with-blas="-lmkl_gf_lp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread" \

--with-lapack="-lmkl_gf_lp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread" \

--with-fftw3="-lmkl_gf_lp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread -lm" \

--with-fftw3f="-lmkl_gf_lp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread -lm" \

--disable-java --enable-static

It seems that using intel compiler with mkl is problematic.

GeoTracker and Googe Earth Virtual Tour

 I am creating this kind of movie, and finally, I have succeeded to find the set of tools useful on Linux.

(1) GeoTracker on Android
This allows me to get the positions at each time point located with GPS.
https://play.google.com/store/apps/details?id=com.ilyabogdanovich.geotracker&hl=ja&gl=US

(2) RouteConverterLinux
This allows me to edit the points measured above without losing the timing information.
I tried with Google Maps, but it automatically deletes the timing information.

On Ubuntu 20.04, I needed to install OpenJDK 14, and the command is;
/usr/lib/jvm/java-14-openjdk-amd64/bin/java -jar RouteConverterLinux.jar 

https://www.routeconverter.de/downloads/de

(3) Google Earth
You should be able to google "Virtual Tour", and scrutinize ;)

The movie I am creating:

Poco F3

 I have not written this, but I succeeded to get Xiaomi's Poco F3 on the launch date in Europe via Amazon (8G RAM/ 256 storage ver).

So far so good, and love 48M pixel camera!

PS:
I did not know, but recent android can transfer the environment very easily. The thing I love most is that Google authenticator entries can be transferred using QR code with just one step.

Pixel Experience (Android 11) for Xiaomi A1

 I have recently bought the Poco F3, and my old Xiaomi A1 can be bricked!
So, I became to want to test some custom ROMs.

Finally, I realized the pixel experience only offers Android 11 as an official one, and installed it.

I just followed the instruction and it worked!

I may not use this intensively, but it seems the ROM is quite stable.

So far, the only defect I have found is that there is no way to use the tele-lens.

Inverse QFT

Based on the last QFT , I tried the inverse QFT.

I got the result of  [640, 654, 573, 666, 615, 601, 635, 654, 616, 603, 646, 602, 629, 605, 593, 668].
This means that the initial state of QFFT is recovered by iQFFT.

The codes of QFT and iQFT are as follows (as far as I understand).

namespace QFT_naive {

open Microsoft.Quantum.Canon;

open Microsoft.Quantum.Intrinsic;

open Microsoft.Quantum.Convert;

open Microsoft.Quantum.Arrays;

open Microsoft.Quantum.Measurement;

open Microsoft.Quantum.Math;

@EntryPoint()

operation SayHello() : Unit {

let nQubit = 4;

let maxNum = PowI(2,nQubit);

mutable counts = new Int[maxNum];

mutable theta = new Double[maxNum];

mutable revSeq = new Int[nQubit];

for iQubit in 0..(nQubit-1){

set theta w/= iQubit <- calcTheta(IntAsDouble(iQubit+1));

}

for iQubit in 0..(nQubit-1){

set revSeq w/= iQubit <- nQubit-1-iQubit;

}

for trial in 1..10000 {

use input = Qubit[nQubit] {

// 1/sqrt(8) (|000> + |001> + ,,, + |111>)

for iQubit in 0..(nQubit-1){

H(input[iQubit]);

}

// Star QFT

for iQubit in 0..(nQubit-1){

H(input[iQubit]);

for jQubit in 1..(nQubit-1-iQubit){

(Controlled R1)([input[jQubit+iQubit]],(theta[jQubit],input[iQubit]));

}

}

let numLoopSwap = Floor(IntAsDouble(nQubit)*0.5);

for iQubit in 0..(numLoopSwap-1) {

SWAP(input[iQubit],input[nQubit-1-iQubit]);

}

// Start Inverse QFFT

for iQubit in 0..(numLoopSwap-1) {

SWAP(input[iQubit],input[nQubit-1-iQubit]);

}

for iQubit in revSeq{

for jQubit in 0..(nQubit-1)-iQubit-1{

let kQubit = revSeq[(nQubit-1)-jQubit];

let controlQubit = nQubit-kQubit-1;

let rotationIndex= nQubit-(kQubit+iQubit)-1;

(Controlled R1)([input[controlQubit]],(-1.0*theta[rotationIndex],input[iQubit]));

}

H(input[iQubit]);

}

let res = MultiM(input);

let bit = ResultArrayAsInt(res);

set counts w/= bit <- counts[bit]+ 1;

ResetAll(input);

}

}

for j in 0..(maxNum-1) {

Message(IntAsString(counts[j]));

}

}

function calcTheta(l : Double) : Double {

return 2.0*PI()/(PowD(2.0,l));

}

}

Normal distribution of dice rolls: Python animation

 Since my wife started learning statistics, and I wanted to give her some intuitive understanding of the normal distribution, standard deviation, 95% confidence region etc, I made a small python code which generates an animation of the histogram of dice rolls.

This could help others time ;)

(Distribute under MIT license)

import random

import matplotlib.pyplot as plt

import numpy as np

from scipy.stats import norm

def main():

numTrial = 1000000

numDice = 20

sumResults = []

bin=np.arange(numDice-0.5, numDice*6+1.5, 1)

print(bin)

fig = plt.figure()

iFig = 1

for iTrial in range(numTrial):

sumDice = 0

for iDice in range(numDice):

sumDice += random.randint(1,6)

sumResults.append(sumDice)

if iTrial % 100 == 0:

ax = fig.add_subplot()

ax.set_title(str(numDice)+" dices case")

ax.set_xlabel("Sum of rolls")

ax.set_ylabel("Frequency")

n, bins, patches = plt.hist(sumResults,bins=bin)

loc,scale = norm.fit(sumResults)

plt.plot(bin,sum(n)*norm.pdf(bin,loc,scale),color="y")

txt = "STD:"+str(round(scale,4))

plt.text(5.0,0.01,txt,size=30)

p25 =norm.ppf(0.025 ,loc,scale)

p975=norm.ppf(0.975,loc,scale)

plt.axvline(x=p25, ymin=0.1, ymax=0.9,color="r")

plt.axvline(x=p975, ymin=0.1, ymax=0.9,color="r")

plt.savefig(str('{:05d}'.format(iFig))+".png")

iFig += 1

plt.pause(.001)

fig.clear()

if __name__ == "__main__":

main()

キルシェンハイン Kirschenhain, 桜の森

 ウィーンのドナウインゼルの北東に桜の木がたくさん植えてあるキルシェンハインという場所があり、そろそろ満開になりそうだったので行ってきました。まだ木が若いので少し迫力にはかけますが、春の気配を感じるには十分でした。

https://goo.gl/maps/4fJwEzJJn8LTLhb89

"lto1: error: ‘-fcf-protection=full’ is not supported for this target" on GCC-10

Probably someone has updated the DGX station system, and g++-10 compiler complains as 

lto1: error: ‘-fcf-protection=full’ is not supported for this target.

To prevent this I needed to give the linker options as;

LD= g++-10 -fopenacc -foffload="-fcf-protection=none -fno-stack-protector"

Colour scheme for scientific computing visualization

 Looks nice, and the author provides Python codes (and some other languages).

I also have a kind of colour blindness (red and green). It would be nice to have this.

https://ai.googleblog.com/2019/08/turbo-improved-rainbow-colormap-for.html

Delete Linux from Dual Boot Environment

I have not known that GPT boot loaders can be deleted in BIOS setting.
So, if you want to delete dual boot Linux, you just need to delete the bootloader and remove the Linux partitions (and expand the Windows partition). 

https://itsfoss.com/uninstall-ubuntu-linux-windows-dual-boot/

Data transfer from Oracle To PostgreSQL via Pandas

 I needed to transfer data from Oracle to PostgreSQL, which are on different machines.
This time, I used Pandas and Pickle on Python 2.7.

Reading from Oracle:

import pickle

import cx_Oracle

import pandas as pd

def main():

conn = cx_Oracle.connect("ID/PW@DB")

tables = ["a","b"]

data = {}

for iTab in tables:

query = "select * from " + iTab

curDF = pd.read_sql(query, con=conn)

data[iTab] = curDF

conn.close()

with open("data.dat","wb") as f:

pickle.dump(data,f)

Writing to PostgreSQL:
(Note1, table name, schema and column names should be in lower case. Otherwise, you will find or need strange double quotation marks.
Note2, you should give chunksize. Otherwise, the python instance will eat up all the memory.
Note3, you need to set up schemas, but tables will be automatically created)

import pickle

import pandas as pd

from sqlalchemy import create_engine

def main():

with open("data.dat","rb") as f:

arrData = pickle.load(f)

engine = create_engine('postgresql://ID:PW@localhost:5432/server')

for iTab in arrData:

print(iTab)

arrData[iTab].to_sql(iTab.lower(), engine,schema="schema".lower(),chunksize=100)

客員共同研究員/Visiting Researcher @東大/The University of Tokyo

As of today, I work for/with Prof. Hiroshi Okuda at the University of Tokyo as a visiting researcher.

Hope that I can bridge the Japanese academia and CTBTO.

Well, it is obviously at no cost. So, if you feel you want to support my research, please feel free to donate with my wishlist on Amazon :)

今日付けで東大　新領域創成科学　奥田研究室の客員共同研究員になりました。

大したことはできませんが、少しでも世の中の役に立てると良いと思っています。

QFT part 5: adjustable number of qubits

As continued from QFT part 4 , the number of qubits is now adjustable.

NQubit = 3 and 4 are working correctly.

I think I made a mistake and updated on 1st of May, 2021.

namespace QFT_naive {

open Microsoft.Quantum.Canon;

open Microsoft.Quantum.Intrinsic;

open Microsoft.Quantum.Convert;

open Microsoft.Quantum.Arrays;

open Microsoft.Quantum.Measurement;

open Microsoft.Quantum.Math;

@EntryPoint()

operation SayHello() : Unit {

let nQubit = 4;

let maxNum = PowI(2,nQubit);

mutable counts = new Int[maxNum];

mutable theta = new Double[maxNum];

for iQubit in 0..(nQubit-1){

set theta w/= iQubit <- calcTheta(IntAsDouble(iQubit+1));

}

for trian in 1..1000 {

use input = Qubit[nQubit] {

// 1/sqrt(8) (|000> + |001> + ,,, + |111>)

for iQubit in 0..(nQubit-1){

H(input[iQubit]);

}

// Star QFT

for iQubit in 0..(nQubit-1){

H(input[iQubit]);

for jQubit in 1..(nQubit-1-iQubit){

(Controlled R1)([input[jQubit+iQubit]],(theta[jQubit],input[iQubit]));

}

}

//H(input[0]);

//(Controlled R1)([input[1]],(theta[1],input[0]));

//(Controlled R1)([input[2]],(theta[2],input[0]));

//H(input[1]);

//(Controlled R1)([input[2]],(theta[1],input[1]));

//H(input[2]);

let numLoopSwap = Floor(IntAsDouble(nQubit)*0.5);

for iQubit in 0..(numLoopSwap-1) {

SWAP(input[iQubit],input[nQubit-1-iQubit]);

}

//SWAP(input[2],input[0]);

let res = MultiM(input);

let bit = ResultArrayAsInt(res);

set counts w/= bit <- counts[bit]+ 1;

ResetAll(input);

}

}

for j in 0..(maxNum-1) {

Message(IntAsString(counts[j]));

}

}

function calcTheta(l : Double) : Double {

return 2.0*PI()/(PowD(2.0,l));

}

}

QFT part 4: Using R1 gate, instead of S and T.

Since I do want to extend the number of qubits from three to more, I have re-written the code with R1 gate, in stead of S and T.

Lessons learned: Functions and operations can be defined. Here a function which returns Double is created. Operations might be used for qubits. Not sure now.

namespace QFT_naive {

open Microsoft.Quantum.Canon;

open Microsoft.Quantum.Intrinsic;

open Microsoft.Quantum.Convert;

open Microsoft.Quantum.Arrays;

open Microsoft.Quantum.Measurement;

open Microsoft.Quantum.Math;

@EntryPoint()

operation SayHello() : Unit {

mutable counts = new Int[8];

let theta_1 = calcTheta(1.0);

let theta_2 = calcTheta(2.0);

let theta_3 = calcTheta(3.0);

for trian in 1..1000 {

use input = Qubit[3] {

// 1/sqrt(8) (|000> + |001> + ,,, + |111>)

H(input[0]);

H(input[1]);

H(input[2]);

// Star QFT

// R1 = Z, R2=S, R3=T

H(input[0]);

//(Controlled S)([input[1]],input[0]);

//(Controlled T)([input[2]],input[0]);

(Controlled R1)([input[1]],(theta_2,input[0]));

(Controlled R1)([input[2]],(theta_3,input[0]));

H(input[1]);

//(Controlled S)([input[2]],input[1]);

(Controlled R1)([input[2]],(theta_2,input[1]));

H(input[2]);

SWAP(input[2],input[0]);

let res = MultiM(input);

let bit = ResultArrayAsInt(res);

set counts w/= bit <- counts[bit]+ 1;

ResetAll(input);

}

}

for j in 0..7 {

Message(IntAsString(counts[j]));

}

}

function calcTheta(l : Double) : Double {

return 2.0*PI()/(PowD(2.0,l));

}

}

 

QFT part 3. With ResultArrayAsInt and MultiM.

Still playing around QFFT. Improved the last one with ResultArrayAsInt and MultiM.
Those allow me to shorten the code a bit, and may allow me to go the 4 qubit case.

One thing I am not quite sure is, which end is the LSB?
I suspect 0 is the LSB, but no clue.

namespace QFT_naive {

open Microsoft.Quantum.Canon;

open Microsoft.Quantum.Intrinsic;

open Microsoft.Quantum.Convert;

open Microsoft.Quantum.Arrays;

open Microsoft.Quantum.Measurement;

@EntryPoint()

operation SayHello() : Unit {

Message("qFFT");

mutable counts = new Int[8];

for trian in 1..1000 {

use input = Qubit[3] {

// 1/sqrt(8) (|000> + |001> + ,,, + |111>)

H(input[0]);

H(input[1]);

H(input[2]);

// Star QFT

// R1 = Z, R2=S, R3=T

H(input[0]);

(Controlled S)([input[1]],input[0]);

(Controlled T)([input[2]],input[0]);

H(input[1]);

(Controlled S)([input[2]],input[1]);

H(input[2]);

SWAP(input[2],input[0]);

let res = MultiM(input);

let bit = ResultArrayAsInt(res);

set counts w/= bit <- counts[bit]+ 1;

ResetAll(input);

}

}

for j in 0..7 {

Message(IntAsString(counts[j]));

}

}

}

I think 0 is LSB because when I modify the code as follows; 5th element (namely 100=4) of "count" has "1000"

mutable res = MultiM(input);

set res w/= 2 <- One;

let bit = ResultArrayAsInt(res);

DOI for my newly published paper

 I have just notified that my newly published paper got a DOI.
https://doi.org/10.1142/S0219876221500274

Actually, I got the notification from ORCID, not from the journal. Things are getting well organized, but feel some skew.

Jug, Parallel framework in Python

 Minor tips for problems I was into in Jug, a parallel computing framework in Python.
I may be wrong, and if you know I misunderstand, please leave comments:

1. Tasks can be created only outside functions.
   Probably, this is common to create the main function with, if __name__ = __main__: main()
   I did this, and called functions decorated with @TaskGenerator. However, jug could not launch tasks that are created in the main function. Instead, when I remove the main function, and write lines on the ground level (no indent level), the same code worked.
2. Multi-node parallelism can be achieved by just running "jug execute XXX.py" if the directory which stores the python code is shared by NFS. I could not find any clear comment on this, and I wondered how I can achieve multi-node parallelism.
3. An atomic operation for File I/O can be achieved by creating a function decorated with @TaskGenerator. According to the summary of Jug, such a function is executed once. So, we can use it for an atomic operation.

Otherwise, I think Jug is excellent, and easy to use!

Still a question:

I do not know what jugfile can do.
On page17 of the documentation(https://jug.readthedocs.io/_/downloads/en/latest/pdf/), a post-process script import jugfile. I think the tasks in the main script make outputs to the jugfile, but do not know how to do.

Upper and lower limit in matplotlib contourf

When the data range in matplotlib contourf does not correspond to the range of colorbar, matplotlib automatically adjust the range, even when set_clim is set manually.

(In other words, I needed a fixed colorbar regardless of the data)

The solution is set "levels" as follows;

interval = np.linspace(2.5,4.5, num=32, endpoint=True)

img=plt.contourf(LON, LAT, DATA,

vmin=2.5,vmax=4.5,

cmap = "jet",

levels=interval,

transform=ccrs.PlateCarree())

Q# without VS code

 Not sure the reason, but I cannot install Q# template with newly install Ubuntu 20.04 and VS code.
I have just learned that CLI also works.

% dotnet new console -lang Q# -o [ProjectName] % cd [ProjectName]

% dotnet build % dotnet run

I needed to install .NET-SDK 3.1, instead of 5.0.

福島の事故がもたらしたもの

福島の事故を軽く見せたい人が一定数おりますが、住めない土地を作ったことについてはどう考えるのか聞きたいものです。

 【福島から伝えたい】”ハゲと白髪”集う「じじい部隊」に問う「原発があって幸せだった?」記者が見た10年

http://web.archive.org/web/20210130064836/https://news.yahoo.co.jp/articles/e05797830e256aece6056c6a4cd8b343d29037a8

Give select grants on all tables in my schema

 The following code gives the select grant on all the tables in your schema to "test2" user.

BEGIN

FOR t IN (SELECT * FROM user_tables)

LOOP

EXECUTE IMMEDIATE 'GRANT SELECT ON ' || t.table_name || ' TO test2';

END LOOP;

END;

/

(last slash "/" is important)

Just a copy from;
https://sites.google.com/site/nazmulhudadba/grant-select-on-all-tables-in-a-specific-schema-to-a-user

Pandas on Python 2.6

We are still using Python2.6 on RH6, and packages are getting out-dated.

Pandas is one of those, and need the following command to be installed.

% easy_install 'pandas==0.17.1'

This seems the latest version for Python 2.6

And also numexpr is necessary for the "query" but the latest version numexpr cannot be installed.

% easy_install 'pandas==2.4'

at least works.

Convert datetime to epoch in a old Python (2.6)

 Since I always forget how to convert datetime to epoch, I just want to write a memo.

In Python
>>> import datetime
>>> d = datetime.datetime(2019,1,1,0,0,0)
>>> print(d)

2019-01-01 00:00:00

>>> epoch = d.strftime("%s")
>>> print(epoch)

1546300800

On Bash

% date -d @1546300800

Tue Jan  1 00:00:00 GMT 2019

OpenMP Offloading on Intel Core i7

 I am thinking if OpenMP Offloading works with XcalableMP, and desire to play around it.

Luckily, my laptop has Core i7-8665U, and Intel's OneAPI compiler provides us with openmp-offloading.

Software installation is quite straightforward;
https://software.intel.com/content/www/us/en/develop/articles/installing-intel-oneapi-toolkits-via-apt.html#aptpkg

And the first step guide is;
https://software.intel.com/content/www/us/en/develop/documentation/get-started-with-cpp-fortran-compiler-openmp/top.html

Surely working!