text:
const element = ["H",  "He", "Li", "Be", "B",  "C",  "N",  "O",  "F",  "Ne",
                 "Na", "Mg", "Al", "Si", "P",  "S",  "Cl", "Ar", "K",  "Ca",
                 "Sc", "Ti", "V",  "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",
                 "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y",  "Zr",
                 "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn",
                 "Sb", "Te", "I",  "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd",
                 "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",
                 "Lu", "Hf", "Ta", "W",  "Re", "Os", "Ir", "Pt", "Au", "Hg",
                 "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th",
                 "Pa", "U",  "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm",
                 "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds",
                 "Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og"]

const reg = r"[0-9]+"

struct bit2b
    a::Int64
    b::Int64
    c::Int64
    d::Int64
end

struct ifph
    i::Int64
    f::Int64
    phase::Bool
end

struct Jpninfo
    fac::Float64
    pjump::Array{ifph,1}
    njump::Array{ifph,1}
end

struct T1info
    f::Int64
    coef::Float64
end

struct MiMf
    Mi::Int64
    Mf::Int64
    fac::Float64
end

"""
    main(sntf,target_nuc,num_ev,target_J;
         save_wav=false,q=1,is_block=false,is_show=false,num_history=3,lm=100,ls=20,tol=1.e-6,
         in_wf="",mdimmode=false,calc_moment = false,gfactors = [1.0,0.0,5.586,-3.826],effcharge=[1.5,0.5])

Digonalize the model-space Hamiltonian

# Arguments:

- `sntf`:       path to input interaction file (.snt fmt)
- `target_nuc`: target nucleus
- `num_ev`:     number of eigenstates to be evaluated
- `target_J`:   target total J (specified by e.g. [0]). Set to [] if you want lowest states with any J.
   Note that J should be doubled (J=0=>[0], J=1/2=>[1], J=1=>[2],...)

# Optional arguments:
- `q=1`              block size for Block-Lanczos methods
- `is_block=false`   whether or not to use Block algorithm
- `save_wav=false`   whether or not to save wavefunction file
- `is_show = true`   to show elapsed time & allocations
- `lm = 100`         number of Lanczos vectors to store
- `ls = 15`          number of vectors to be used for Thick-Restart
- `tol= 1.e-6`       tolerance for convergence check in the Lanczos method
- `in_wf=""`      path to initial w.f. (for preprocessing)
- `mdimmode=false`   `true` => calculate only the M-scheme dimension
- `calc_moment=false`  `true` => calculate mu&Q moments
- `calc_entropy=false`  `true` => calculate entanglement entropy
- `gfactors=[1.0,0.0,5.586,-3.826]` angular momentum and spin g-factors
- `effcgarge=[1.5,0.5]` effective charges
"""
function main(sntf,target_nuc,num_ev,target_J;save_wav=false,
              q=1,is_block=false,is_show=false,
              num_history=3,lm=100,ls=20,tol=1.e-6,
              in_wf="",mdimmode=false,
              calc_moment = false,
              calc_entropy = false,
              gfactors = [1.0,0.0,5.586,-3.826],
              effcharge=[1.5,0.5])
    to = TimerOutput()
    Anum = parse(Int64, match(reg,target_nuc).match)
    lp,ln,cp,cn,massop,Aref,pow,p_sps,n_sps,SPEs,olabels,oTBMEs,labels,TBMEs = readsnt(sntf,Anum)
    massformula = 1
    if 16<= Anum <= 40; massformula = 2;end
    ## massformula=2: J. Blomqvist and A. Molinari, Nucl. Phys. A106, 545 (1968).
    ## we use this for the sd-shell nuclei
    hw, bpar = init_ho_by_mass(Anum,massformula)

    if length(target_J) > 1;
        println("warn! Multiple J is not supported now.");exit()
    end
    Mtot = 0;if Anum % 2 != 0; Mtot = 1;end
    tJ = -1; eval_jj = -1.0
    if length(target_J) > 0; Mtot = minimum(target_J);tJ=target_J[1];
        eval_jj = 0.5*tJ*(tJ/2+1) ; end
    if Anum % 2 != Mtot % 2; println("invalid targetJ $tJ");exit();end
    target_el = replace.(target_nuc, string(Anum)=>"")
    Z,N,vp,vn = getZNA(target_el,Anum,cp,cn)
    mstates_p, mstates_n,mz_p,mz_n = def_mstates(p_sps,n_sps)
    pbits,nbits,jocc_p,jocc_n,Mps,Mns,tdims = occ(p_sps,mstates_p,mz_p,vp,
                                                  n_sps,mstates_n,mz_n,vn,Mtot)
    lblock=length(pbits)
    mdim = tdims[end]; if mdim==0;exit();end

    mdim_print(target_nuc,Z,N,cp,cn,vp,vn,mdim,tJ)
    if mdimmode; return nothing;end

    @timeit to "prep. 1bjumps" begin
        ## bit representation of Hamiltonian operators
        bV1,V1 = HbitT1(p_sps,n_sps,mstates_p,mstates_n,labels,TBMEs)
        bVpn,Vpn,delMs = Hbitpn(p_sps,n_sps,mstates_p,mstates_n,labels[3],TBMEs[3])

        ## storing two-body jumps for pp/nn 2b interaction
        ppinfo = prep_pp(mstates_p,pbits,bV1[1],V1[1])
        nninfo = prep_nn(mstates_n,nbits,bV1[2],V1[2])
        bV1 = nothing

        ## storing one-body jumps for pn 2b interaction
        l_pbit = length(mstates_p);l_nbit = length(mstates_n)
        bis,bfs,p_NiNfs,n_NiNfs,num_task = prep_pn(lblock,tdims,l_pbit,l_nbit,
                                                   pbits,nbits,Mps,delMs,bVpn,Vpn)
        bVpn=nothing
    end
    ## distribute task
    block_tasks = make_distribute(num_task)

    @timeit to "Jcalc." begin
        oPP,oNN,oPNu,oPNd = prep_J(tdims,p_sps,n_sps,mstates_p,mstates_n,
                                   pbits,nbits)
        Js = [ 0.5*Mtot*(0.5*Mtot+1) for i = 1:num_ev]
        Jtasks = zeros(Int64,lblock)
        for i = 1:lblock
            Jtasks[i] = length(pbits[i])*length(nbits[i])
        end
        Jidxs = make_distribute_J(Jtasks)
    end

    @timeit to "Lanczos" begin
        en =[ [1.e+4 for j=1:num_ev] for i = 1:num_history]
        Rvecs = [ zeros(Float64,mdim) for i=1:num_ev]
        Tmat = zeros(Float64,lm,lm)
        vks =[]; uks=[];itmin = 1; elit=1
        doubleLanczos = false
        if tJ !=-1; doubleLanczos = true;end

        if is_block #Thick-Restart double Block Lanczos: TRdBL
            ls_sub = div(ls,q)
            vks = [ zeros(Float64,q,mdim) for i=1:div(lm,q)]
            uks = [ zeros(Float64,q,mdim) for i=1:ls_sub*q]
            Beta_H = zeros(Float64,q,q)
            V = vks[1]
            if in_wf !=""
                try
                    read_appwav(in_wf,mdim,V,q,true)
                    bl_QR!(V',Beta_H,mdim,q)
                catch
                    println("error @preprocessing: failed to read appwav")
                    initialize_bl_wav(mdim,q,vks[1])
                    bl_QR!(V',Beta_H,mdim,q)
                end
            else
                initialize_bl_wav(mdim,q,vks[1])
                bl_QR!(V',Beta_H,mdim,q)
            end
            elit = TRBL(q,vks,uks,Tmat,Beta_H,pbits,nbits,jocc_p,jocc_n,SPEs,
                         ppinfo,nninfo,bis,bfs,block_tasks,
                         p_NiNfs,n_NiNfs,Mps,delMs,Vpn,tdims,
                         eval_jj,oPP,oNN,oPNu,oPNd,Jidxs,
                         num_ev,num_history,lm,ls_sub,en,tol,to,doubleLanczos)
        else #Thick Restart (double) Lanczos: TRL
            vks = [ zeros(Float64,mdim) for i=1:lm]
            uks = [ zeros(Float64,mdim) for i=1:ls]
            initialize_wf(vks[1],"rand",tJ,mdim)
            elit = TRL(vks,uks,Tmat,itmin,
                       pbits,nbits,jocc_p,jocc_n,SPEs,
                       ppinfo,nninfo,bis,bfs,block_tasks,
                       p_NiNfs,n_NiNfs,Mps,delMs,Vpn,
                       eval_jj,oPP,oNN,oPNu,oPNd,Jidxs,
                       tdims,num_ev,num_history,lm,ls,en,tol,to,doubleLanczos)
        end
    end

    @timeit to "Rvecs" begin
        vals,vecs = eigen(@views Tmat[1:elit*q,1:elit*q])
        @inbounds for (nth,Rvec) in enumerate(Rvecs)
            if is_block == false
                @inbounds for k=1:length(vals)
                    Rvec .+= vecs[k,nth] .* vks[k]
                end
            else
                @inbounds for k=1:length(vals)
                    it = div(k-1,q)
                    b = k - q*it
                    Rvec .+= @views vecs[k,nth] .* vks[it+1][b,:]
                end
            end
            Rvec .*= 1.0/sqrt(dot(Rvec,Rvec))
        end
    end
    vt = zeros(Float64,mdim)
    for (nth,Rv) in enumerate(Rvecs)
        vt .= 0.0
        operate_J!(Rv,vt,pbits,nbits,tdims,
                   Jidxs,oPP,oNN,oPNu,oPNd)
        Js[nth] += dot(Rv,vt)
    end
    totJs = J_from_JJ1.(Js)
    ### entanglement entropy
    if calc_entropy
        entropy(Rvecs,pbits,nbits,tdims,to)
    end
    ### mu&Q-moment
    tx_mom =""
    if calc_moment
        tx_mom = eval_moment(Mtot,Rvecs,totJs,p_sps,n_sps,
                             mstates_p,mstates_n,tdims,
                             jocc_p,jocc_n,pbits,nbits,bpar,
                             gfactors,effcharge)
    end
    if save_wav
        @timeit to "I/O" begin
            csnt = split(split(sntf,"/")[end],".")[1]
            oupf="./"*target_nuc*"_"*csnt*".wav"
            if tJ != -1;oupf="./"*target_nuc*"_"*csnt*"_j"*string(tJ)*".wav";end
            writeRitzvecs(mdim,Mtot,en[1],totJs,Rvecs,oupf)
        end
    end
    if is_show
        show(to, allocations = true,compact = true);println("")
    end
    println("sntf: $sntf")
    println("J $totJs")
    print("En. ");map(x -> @printf("%9.3f ",x), en[1])
    print("\nEx. ");map(x -> @printf("%9.3f ",x),[en[1][i]-en[1][1] for i=1:num_ev])
    print("\n")
    if tx_mom != ""
        println(tx_mom)
    end
    return nothing
end

function rm_comment(lines)
    nlines = []
    for line in lines
        line = strip(line)
        if length(line)>1
            if startswith(line,"!")
                continue
            end
        end
        push!(nlines,line)
    end
    return nlines
end

function rm_nan(array)
    na = []
    for tmp in array
        if tmp != "";push!(na,tmp); end
    end
    return na
end

"""
    readsnt(sntf,Anum)

To read interaction file in ".snt" format.
- `sntf`: path to the interaction file
- `Anum`: mass number (used for "scaling" of TBMEs)

!!! note
    The current version supports only "properly ordered" interaction file in .snt format.

    A .snt file can be ordered to be a<=b,c<=d,a<=c for V(abcd;J) by the Python script "ShellModel.jl/src/make_ordered_snt.py"(, which will be replaced by Julia implementation...).

"""
function readsnt(sntf,Anum)
    f = open(sntf,"r");tlines = readlines(f);close(f)
    lines = rm_comment(tlines)
    line = lines[1]
    lp,ln,cp,cn = map(x->parse(Int,x),rm_nan(split(line," ")))
    p_sps = [[0,0]];deleteat!(p_sps,1)
    n_sps = [[0,0]];deleteat!(n_sps,1)
    @inbounds for i = 1:lp
        ith,n,l,j,tz = map(x->parse(Int,x),rm_nan(split(lines[1+i]," "))[1:5])
        push!(p_sps,[n,l,j,tz])
    end
    @inbounds for i = 1:ln
        ith, n,l,j,tz = map(x->parse(Int,x),rm_nan(split(lines[1+i+ln]," "))[1:5])
        push!(n_sps,[n,l,j,tz])
    end
    nsp,zero = map(x->parse(Int,x),rm_nan(split(lines[1+ln+lp+1]," "))[1:2])
    SPEs = [ [0.0 for i=1:lp],[0.0 for i=1:ln]]
    @inbounds for i = 1:nsp
        idx=0; j=i
        if i<=lp;idx=1;else;idx=2;j-=lp;end
        SPEs[idx][j] =parse(Float64,rm_nan(split(lines[1+ln+lp+1+i]," "))[3])
    end
    ntbme,massop,Aref,p = rm_nan(split(lines[1+ln+lp+1+nsp+1]," "))
    ntbme = parse(Int,ntbme); massop=parse(Int,massop)
    Aref = parse(Int,Aref); p = parse(Float64,p)
    labels = [ [ [0,0] ] for i=1:3]
    olabels = [ [0,0] ]
    for i=1:3; deleteat!(labels[i],1);end
    deleteat!(olabels,1)
    TBMEs=[ Float64[] for i =1:3] #pp/nn/pn
    oTBMEs= Float64[]
    @inbounds for ith = 1:ntbme
        tmp = rm_nan(split(lines[1+ln+lp+1+nsp+1+ith], " "))
        i = tmp[1]; j = tmp[2]; k = tmp[3]; l =tmp[4]; totJ = tmp[5];  TBME= tmp[6]
        i = parse(Int,i);j = parse(Int,j);k = parse(Int,k);l = parse(Int,l);
        nth = 0
        if i<=lp && j<=lp
            nth = 1
        elseif i>lp && j > lp
            nth = 2
        elseif i<=lp && j>lp
            nth = 3
        else
            println("err");exit()
        end
        ## snt file must be "ordered"; a<=b & c=d & a<=c
        TBME = parse(Float64,TBME)
        if massop == 1; TBME*= (Anum/Aref)^(p);end
        if unique([i,j]) != unique([k,l])
            push!(labels[nth],[k,l,i,j,parse(Int,totJ),ith])
            push!(TBMEs[nth],TBME)
        end
        push!(labels[nth],[i,j,k,l,parse(Int,totJ),ith])
        push!(TBMEs[nth],TBME)

        push!(olabels,[i,j,k,l,parse(Int,totJ),ith])
        push!(oTBMEs,TBME)
    end
    return lp,ln,cp,cn,massop,Aref,p,p_sps,n_sps,SPEs,olabels,oTBMEs,labels,TBMEs
end

"""
    def_mstate(p_sps, n_sps)

To define the single particle states specified by n,l,j,mz
"""
function def_mstates(p_sps,n_sps)
    mstates_p = [[1]]; mstates_n = [[1]]; mz_p = Int64[]; mz_n = Int64[]
    deleteat!(mstates_p,1); deleteat!(mstates_n,1)
    for (pidx,tsps) in enumerate(p_sps)
        n,l,j,tz = tsps
        for mz = -j:2:j;push!(mstates_p,[n,l,j,tz,mz,pidx]);push!(mz_p,mz);end
    end
    for (nidx,tsps) in enumerate(n_sps)
        n,l,j,tz = tsps
        for mz = -j:2:j;push!(mstates_n,[n,l,j,tz,mz,nidx]);push!(mz_n,mz);end
    end
    return mstates_p, mstates_n,mz_p,mz_n
end

"""
    all_perm!(ln::Int64,num_valence::Int64,occs::Array{Array{Bool,1}})

make all possible permutation of 'bits'

Example:
If 2 protons and 1 neutron are in the 0p-shell space,
valence orbits(0p1/2,0p3/2) => -1/2, 1/2, -3/2, -1/2, 1/2, 3/2

configurations are represented like:

 proton: 000011, 000101, ..., 110000

neutron: 000001, 000010, ..., 100000
"""
function all_perm!(ln::Int64,num_valence::Int64,
                   occs::Array{Array{Bool,1}})
    for (i,tcomb) in enumerate(collect(combinations(collect(1:ln),num_valence)))
        @inbounds for nth in tcomb
            occs[i][nth] = 1
        end
    end
    return nothing
end

function Mcount!(ln::Int64,mzs::Array{Int64,1},
                 occ::Array{Bool,1},
                 Mret::Array{Int64,1})
    Mret[1] = 0
    @inbounds for i = 1:ln
        if occ[i]
            Mret[1] += mzs[i]
        end
    end
    return nothing
end

function possible_mz(nljtz,mstates)
    n,l,j,tz = nljtz
    mzs = Int64[]; midxs=Int64[]
    for mz = -j:2:j
        push!(mzs,mz)
        for k = 1:length(mstates)
            if @views mstates[k][1:5] == [n,l,j,tz,mz]
                push!(midxs,k)
                break
            end
        end
    end
    return j,mzs, midxs
end

function initialize_tvec!(tvec::Array{Bool,1})
    tvec .= false
    return nothing
end

"""
    function HbitT1(p_sps::Array{Array{Int64,1}},n_sps::Array{Array{Int64,1}},
                mstates_p::Array{Array{Int64,1},1},mstates_n::Array{Array{Int64,1},1},
                labels::Array{Array{Array{Int64,1},1},1},TBMEs::Array{Array{Float64,1}})

make bit representation of T=1 (proton-proton&neutron-neutron) interactions for each {m_z}
"""
function HbitT1(p_sps::Array{Array{Int64,1}},
                n_sps::Array{Array{Int64,1}},
                mstates_p::Array{Array{Int64,1},1},
                mstates_n::Array{Array{Int64,1},1},
                labels::Array{Array{Array{Int64,1},1},1},
                TBMEs::Array{Array{Float64,1}})
    lp = length(mstates_p)
    ln = length(mstates_n)
    bV1 = [ bit2b[] for i=1:2 ]
    V1 = [ Float64[] for i=1:2]
    mstates = [mstates_p,mstates_n]
    sps = [p_sps,n_sps]
    loffs = [ 0, length(p_sps)]

    for vrank =1:2 #pp:1, nn:2
        loff = loffs[vrank]
        vecs= [ [ [ false for i = 1:lp] for j=1:2],
                [ [ false for i = 1:ln] for j=1:2]]
        blist = bit2b[]
        Vs=Float64[]
        @inbounds for (i,ME) in enumerate(TBMEs[vrank])
            a,b,c,d,totJ,dummy = labels[vrank][i]
            J2  = 2*totJ
            ja,ma_s,ma_idxs = possible_mz(sps[vrank][a-loff],mstates[vrank])
            jb,mb_s,mb_idxs = possible_mz(sps[vrank][b-loff],mstates[vrank])
            jc,mc_s,mc_idxs = possible_mz(sps[vrank][c-loff],mstates[vrank])
            jd,md_s,md_idxs = possible_mz(sps[vrank][d-loff],mstates[vrank])
            @inbounds for (ic,mc) in enumerate(mc_s)
                @inbounds for (id,md) in enumerate(md_s)
                    if c == d && mc >= md; continue;end
                    if abs(mc + md) > J2; continue;end
                    M_ani = mc + md
                    initialize_tvec!(vecs[vrank][1]);vecs[vrank][1][mc_idxs[ic]] = true
                    bit_c = bitarr_to_int(vecs[vrank][1])
                    initialize_tvec!(vecs[vrank][1]); vecs[vrank][1][md_idxs[id]] = true
                    bit_d = bitarr_to_int(vecs[vrank][1])
                    @inbounds for (ia,ma) in enumerate(ma_s)
                        @inbounds for (ib,mb) in enumerate(mb_s)
                            if a == b && ma>=mb; continue;end
                            if ma + mb != M_ani;continue;end
                            initialize_tvec!(vecs[vrank][2]);vecs[vrank][2][ma_idxs[ia]] = true
                            bit_a = bitarr_to_int(vecs[vrank][2])
                            initialize_tvec!(vecs[vrank][2]);vecs[vrank][2][mb_idxs[ib]] = true
                            bit_b = bitarr_to_int(vecs[vrank][2])
                            CG1 = clebschgordan(Float64,ja//2,ma//2,jb//2,mb//2, J2//2, M_ani//2)
                            CG2 = clebschgordan(Float64,jc//2,mc//2,jd//2,md//2, J2//2, M_ani//2)
                            tl = bit2b(bit_a,bit_b,bit_c,bit_d)
                            if ( tl in blist) == false
                                push!(blist,tl)
                                push!(Vs, ME * sqrt( (1.0+deltaf(a,b)) *(1.0+deltaf(c,d)) ) * CG1 * CG2)
                                continue
                            end
                            @inbounds for kk = 1:length(blist)
                                if blist[kk] == tl
                                    Vs[kk] += ME * sqrt( (1.0+deltaf(a,b)) *(1.0+deltaf(c,d)) ) * CG1 * CG2
                                    break
                                end
                            end
                        end
                    end
                end
            end
        end
        bV1[vrank] = blist
        V1[vrank] = Vs
    end
    return bV1,V1 #bVpp,Vpp,bVnn,Vnn
end

"""
    Hbitpn(p_sps::Array{Array{Int64,1}},n_sps::Array{Array{Int64,1}},
           mstates_p::Array{Array{Int64,1},1},mstates_n::Array{Array{Int64,1},1},
           labels::Array{Array{Int64,1}},TBMEs::Array{Float64,1},zeroME=false)

make bit representation of T=0 (proton-neutron) interactions for each {m_z}
"""
function Hbitpn(p_sps::Array{Array{Int64,1}},
                n_sps::Array{Array{Int64,1}},
                mstates_p::Array{Array{Int64,1},1},
                mstates_n::Array{Array{Int64,1},1},
                labels::Array{Array{Int64,1}},
                TBMEs::Array{Float64,1},
                zeroME=false)
    lp = length(mstates_p); ln = length(mstates_n)
    loff = length(p_sps)
    Mzs = Int64[]
    for j = 1:lp
        for i = 1:lp
            push!(Mzs,mstates_p[i][5]-mstates_p[j][5])
        end
    end
    unique!(Mzs);sort!(Mzs,rev=true)
    lenMz = length(Mzs)

    bVpn= [ [[1,1]] for i=1:lenMz ]
    for i = 1:lenMz;deleteat!(bVpn[i],1);end
    Vpn=[ Float64[ ] for i=1:lenMz]

    ret = [-1,-1,-1]
    vec_ani_p = [false for i = 1:lp];vec_cre_p = [false for i = 1:lp]
    vec_ani_n = [false for i = 1:ln];vec_cre_n = [false for i = 1:ln]

    @inbounds for (i,ME) in enumerate(TBMEs)
        a,b,c,d,totJ,dummy = labels[i]
        J2  = 2*totJ
        ja,ma_s,ma_idxs = possible_mz(p_sps[a],mstates_p)
        jc,mc_s,mc_idxs = possible_mz(p_sps[c],mstates_p)
        jb,mb_s,mb_idxs = possible_mz(n_sps[b-loff],mstates_n)
        jd,md_s,md_idxs = possible_mz(n_sps[d-loff],mstates_n)
        ja = ja//2; jb=jb//2;jc = jc//2; jd=jd//2
        for (ic,mc) in enumerate(mc_s)
            initialize_tvec!(vec_ani_p); vec_ani_p[mc_idxs[ic]] = true
            bit_c = bitarr_to_int(vec_ani_p)
            for (ia,ma) in enumerate(ma_s)
                initialize_tvec!(vec_cre_p); vec_cre_p[ma_idxs[ia]] = true
                bit_a = bitarr_to_int(vec_cre_p)
                Mp = ma - mc
                bisearch!(Mzs,Mp,ret);idx = ret[1]
                tV = Vpn[idx];  bV = bVpn[idx]
                for (id,md) in enumerate(md_s)
                    if abs(mc + md) > J2; continue;end
                    initialize_tvec!(vec_ani_n); vec_ani_n[md_idxs[id]] = true
                    bit_d = bitarr_to_int(vec_ani_n)
                    CG1 = clebschgordan(Float64,jc,mc//2,jd,md//2,J2//2,(mc+md)//2)
                    tfac = ME .* CG1
                    for (ib,mb) in enumerate(mb_s)
                        if mb - md + Mp != 0; continue;end
                        initialize_tvec!(vec_cre_n); vec_cre_n[mb_idxs[ib]] = true
                        bit_b = bitarr_to_int(vec_cre_n)
                        fac = tfac .* clebschgordan(Float64,ja,ma//2,jb,mb//2,J2//2,(ma+mb)//2)
                        tl = [bit_a,bit_b,bit_c,bit_d]
                        if (tl in bV) == false
                            push!(bV,copy(tl))
                            push!(tV,fac)
                            continue
                        end
                        @inbounds for kk = 1:length(bV)
                            if bV[kk] == tl
                                tV[kk] += fac
                                break
                            end
                        end
                    end
                end
            end
        end
    end
    return bVpn,Vpn,Mzs
end

function bitarr_to_int(arr::Array{Bool,1}, val = 0)
    v = 2^(length(arr)-1)
    for i in eachindex(arr)
        val += v*arr[i]
        v >>= 1
    end
    return val
end
function bitarr_to_int!(arr::Array{Bool,1},ret)
    ret[1] = 0
    ret[2] = 2^(length(arr)-1)
    for i in eachindex(arr)
        ret[1] += ret[2]*arr[i]
        ret[2] >>= 1
    end
    return nothing
end

function deltaf(i::Int64,j::Int64)
    ifelse(i==j,1.0,0.0)
end

"""
    occ(p_sps::Array{Array{Int64,1}},mstates_p::Array{Array{Int64,1}},mzp::Array{Int64,1},num_vp::Int64,
        n_sps::Array{Array{Int64,1}},mstates_n::Array{Array{Int64,1}},mzn::Array{Int64,1},num_vn::Int64,Mtot::Int64)

prepare bit representations of proton/neutron Slater determinants => pbits/nbits

jocc_p, jocc_n: corresponding occupation numbers for a "j" state,
which is used for one-body operation and OBTDs.

Mps/Mns: total M for proton/neutron "blocks"

For 6Li in the p shell and M=0, Mps = [-3,-1,1,3] & Mns = [3,1,-1,-3]
blocks => [ (Mp,Mn)=(-3,3),(Mp,Mn)=(-1,1),...]

tdims: array of cumulative number of M-scheme dimensions for "blocks"

tdims =[ # of possible configurations of (-3,3),
         # of possible configurations of (-1,1),...]
"""
function occ(p_sps::Array{Array{Int64,1}},
             mstates_p::Array{Array{Int64,1}},
             mzp::Array{Int64,1},num_vp::Int64,
             n_sps::Array{Array{Int64,1}},
             mstates_n::Array{Array{Int64,1}},
             mzn::Array{Int64,1},num_vn::Int64,
             Mtot::Int64)
    lp = length(mzp); ln = length(mzn)
    pDim = binomial(lp,num_vp)
    occs_p =[ [false for j=1:lp] for i = 1:pDim]
    all_perm!(lp,num_vp,occs_p)

    nDim = binomial(ln,num_vn)
    occs_n =[ [false for j=1:ln] for i = 1:nDim]
    all_perm!(ln,num_vn,occs_n)

    mdim = 0;tdims=Int64[];push!(tdims,0)
    Mret = Int64[0]; Mps = Int64[]; Mns = Int64[]
    for i=1:pDim
        Mcount!(lp,mzp,occs_p[i],Mret); push!(Mps, Mret[1])
    end
    for i=1:nDim
        Mcount!(ln,mzn,occs_n[i],Mret); push!(Mns, Mret[1])
    end
    Mps = unique(Mps); Mns = unique(Mns)
    sort!(Mps,rev=false); sort!(Mns,rev=true)
    possidxs = []
    for i = 1:length(Mps)
        for j = 1:length(Mns)
            if Mps[i] + Mns[j] == Mtot
                push!(possidxs,[i,j])
            end
        end
    end
    pbits = [ Int64[] for i=1:length(possidxs) ]
    nbits = [ Int64[] for i=1:length(possidxs) ]
    occ_p_j = [ [[0.0,0.0]] for i=1:length(possidxs) ]
    occ_n_j = [ [[0.0,0.0]] for i=1:length(possidxs) ]
    for i =1:length(possidxs)
        deleteat!(occ_p_j[i],1);deleteat!(occ_n_j[i],1)
    end
    tocc_p_j = zeros(Float64,length(p_sps))
    tocc_n_j = zeros(Float64,length(n_sps))
    for ith = 1:length(possidxs)
        ni,nj = possidxs[ith]
        Mp = Mps[ni]; Mn = Mns[nj]
        if Mp + Mn != Mtot;println("warn");end
        for i = 1:pDim
            Mcount!(lp,mzp,occs_p[i],Mret)
            if Mp != Mret[1]; continue;end
            push!(pbits[ith], bitarr_to_int(occs_p[i]))
            count_jocc!(p_sps,mstates_p,occs_p[i],tocc_p_j)
            push!(occ_p_j[ith],copy(tocc_p_j))
            for j=1:nDim
                Mcount!(ln,mzn,occs_n[j],Mret)
                if Mret[1] + Mp == Mtot
                    mdim += 1
                end
            end
        end
        push!(tdims,mdim)
    end
    for ith = 1:length(possidxs)
        ni,nj = possidxs[ith]
        Mp = Mps[ni]; Mn = Mns[nj]
        for j=1:nDim
            Mcount!(ln,mzn,occs_n[j],Mret)
            if Mret[1] == Mn
                push!(nbits[ith], bitarr_to_int(occs_n[j]))
                count_jocc!(n_sps,mstates_n,occs_n[j],tocc_n_j)
                push!(occ_n_j[ith],copy(tocc_n_j))
            end
        end
    end
    return pbits,nbits,occ_p_j,occ_n_j,Mps,Mns,tdims
end

function TF_connectable(Phi::Int64,bVs::bit2b,TF::Array{Bool,1})
    # bVs = [ bit(Ca^+), bit(Cb^+), bit(Cc), bit(Cd) ]
    if bVs.c != (bVs.c & Phi);TF[1]=false;return nothing;end
    if bVs.d != (bVs.d & Phi);TF[1]=false;return nothing;end
    APhi =  (~bVs.d) & ( (~bVs.c) & Phi)
    if bVs.a != (bVs.a & (~APhi));TF[1]=false;return nothing;end
    if bVs.b != (bVs.b & (~APhi));TF[1]=false;return nothing;end
    TF[1] = true
    return nothing
end

function TF_connectable_1(Phi::Int64,a::Int64,c::Int64,TF::Array{Bool,1})
    # a: creation c: anihilation
    if c != c & Phi;TF[1]=false;return nothing;end
    if a != (a & (~((~c) & Phi)));TF[1]=false;return nothing;end
    TF[1] = true
    return nothing
end

function calc_phase!(Phi::Int64,bVs::bit2b,
                     ln::Int64,ret::Array{Int64,1})
    ret[1] = 0; ret[2] = -1 # ret=[phase,nPhi]
    ## anihilation c
    ret[1] += count_ones(Phi& (~(bVs.c-1)))
    ret[2] = (~bVs.c) & Phi #|phi>:= Cc|v>
    ## anihilation d
    ret[1] += count_ones(ret[2]& (~(bVs.d-1)))
    ret[2] = (~bVs.d) & ret[2] #|phi>:=CdCc|v>
    ## creation b
    ret[1] += count_ones(ret[2]& (~(bVs.b-1)))
    ret[2] = bVs.b | ret[2] #|phi>:=C^+bCdCc|v>
    ## creation a
    ret[1] += count_ones(ret[2]& (~(bVs.a-1)))
    ret[2] = bVs.a | ret[2] #|phi>:=C^+aC^+bCdCc|v>
    ret[1] = (-1)^ret[1]
    return nothing
end

function calc_phase_1!(Phi::Int64,cre::Int64,ani::Int64,ret::Array{Int64,1})
    ret[1] = 1; ret[2] = -1
    ## anihilation #|phi>:=Cc|v>
    ret[1] += count_ones(Phi & (~(ani-1)))
    ret[2] = (~ani) & Phi
    ## creation  #|phi>:=C^+aCc|v>
    ret[1] += count_ones(ret[2] & (~(cre-1)))
    ret[2] = cre | ret[2]
    ret[1] = (-1)^ret[1]
    return nothing
end

function count_jocc!(sps,mstates,occ_bit,tocc)
    tocc .= 0.0
    @simd for i = 1:length(occ_bit)
        if occ_bit[i]
            tocc[mstates[i][6]] += 1.0
        end
    end
    return nothing
end
function getZNA(target_el,Anum,cp,cn)
    Z = 0
    for i = 1:length(element)
        if element[i] == target_el;Z = i; break; end
    end
    N = Anum - Z
    vp = Z - cp; vn = N - cn
    return Z,N,vp,vn
end
function ThickRestart(vks,uks,Tmat,lm,ls)
    vals,vecs = eigen(@views Tmat[1:lm-1,1:lm-1])
    beta = Tmat[lm-1,lm]
    Tmat .= 0.0
    @inbounds for k = 1:ls
        Tmat[k,k] = vals[k]
    end
    @inbounds for (k,uk) in enumerate(uks)
        tmp = beta .* vecs[lm-1,k]
        Tmat[ls+1,k] = tmp; Tmat[k,ls+1] = tmp
        uk .= 0.0
        @inbounds for j=1:lm-1
            axpy!(vecs[j,k],vks[j],uk)
        end
    end
    @inbounds for (k,uk) in enumerate(uks)
        vks[k] .= uk .* (1.0 /sqrt(dot(uk,uk)))
    end
    vks[ls+1] .= vks[lm]
    @inbounds for k = ls+2:lm
        vks[k] .= 0.0
    end
    return nothing
end

function ThickRestart_J(vks,uks,Tmat,lm,ls,eval_jj,Jtol)
    vals,vecs = eigen(@views Tmat[1:lm-1,1:lm-1])
    k = argmin(abs.(vals.-eval_jj))
    beta = Tmat[lm-1,lm]
    Tmat .= 0.0
    Tmat[1,1] = vals[k]
    uk=uks[1]
    tmp = beta .* vecs[lm-1,k]
    Tmat[ls+1,k] = tmp; Tmat[k,ls+1] = tmp
    uk .= 0.0
    @inbounds for j=1:lm-1
        axpy!(vecs[j,k],vks[j],uk)
    end
    vks[1] .= uk .* (1.0/sqrt(dot(uk,uk)))
    vks[2] .= vks[lm]
    @inbounds for k = ls+2:length(vks)
        vks[k] .= 0.0
    end
    return nothing
end

function initialize_wf(v1,method,tJ,mdim)
#    Random.seed!(123)
    if method=="rand" || tJ == -1
        v1 .= randn(mdim,)
        v1 ./= sqrt(dot(v1,v1))
    elseif method == "one"
        v1 .= 1.0 / sqrt(mdim)
    end
    return nothing
end

function initialize_bl_wav(mdim,q,vks)
    Random.seed!(123)
    for i=1:q
        v = @views vks[i,:]
        v .= randn(mdim,)
        if i>1
            for k=i-1:-1:1
                tv = @views vks[k,:]
                v .-= dot(v,tv) .* tv
            end
        end
        tmp = 1.0/sqrt(dot(v,v))
        v.*= tmp
    end
    return nothing
end

function read_wav(inpf,mdim,n;all=false,verbose=false)
    fid =open(inpf)
    neig = read(fid,Int32)
    mtot = read(fid,Int32)
    Es = [read(fid,Float64) for i = 1:neig]
    jj = [read(fid,Int32) for i=1:neig]
    V = [ [read(fid,Float64) for i=1:mdim] for nth=1:neig]
    close(fid)
    if all
        return V,jj
    else
        return V[n],jj[n]
    end
end

function read_appwav(inpf,mdim,V,q,verbose=false)
    fid =open(inpf)
    neig = read(fid,Int32)
    if neig < q; println("#(appwav) must be >= q");exit();end
    mtot = read(fid,Int32)
    Es = [read(fid,Float64) for i = 1:neig]
    jj = [read(fid,Int32) for i=1:neig]
    if verbose
        println("appwav: $inpf")
        print_vec("Es(n=$neig)",Es)
    end
    for j=1:q
        V[j,:] .= [read(fid,Float64) for i=1:mdim]
    end
    close(fid)
    return nothing
end

function ReORTH(it,vtarget,vks)
    @inbounds for l = it:-1:1
        v = vks[l]
        axpy!(-dot(v,vtarget),v,vtarget)
    end
    return nothing
end

function myQR!(Q,R::FA2,
               d1::Int64,d2::Int64) where{FA2<:Array{Float64,2}}
    R .= 0.0
    @inbounds for j = 1:d2
        @inbounds @simd for i = 1:d1
            R[j,j] += Q[i,j] .* Q[i,j]
        end
        R[j,j] = sqrt(R[j,j])
        @inbounds @simd for i = 1:d1
            Q[i,j] /= R[j,j]
        end
        @inbounds for k = j+1:d2
            @simd for i = 1:d1
                R[j,k] += Q[i,j] .* Q[i,k]
            end
            @simd for i = 1:d1
                Q[i,k] -= Q[i,j] .* R[j,k]
            end
        end
    end
    return nothing
end

function bl_QR!(Q,R::FA2,
               d1::Int64,d2::Int64) where{FA2<:Array{Float64,2}}
    R .= 0.0
    @inbounds for j = 1:d2
        q = @views Q[:,j]
        t = sqrt(dot(q,q))
        R[j,j] = t
        q .*= 1.0/t
        @inbounds for k = j+1:d2
            @simd for i = 1:d1
                R[j,k] += q[i] .* Q[i,k]
            end
            @simd for i = 1:d1
                Q[i,k] -= q[i] .* R[j,k]
            end
        end
    end
    return nothing
end

function bisearch!(v::Array{Int64,1}, target::Int64,
                  ret::Array{Int64,1})
    hi = length(v); lo = 1
    ret[2] = 1; ret[3]=length(v)
    @inbounds while ret[2] <= ret[3]
        ret[1] = div(ret[2]+ret[3],2)
        if v[ret[1]] > target
            ret[2] = ret[1]+1
        elseif v[ret[1]] < target
            ret[3] = ret[1]-1
        else
            return nothing
        end
    end
    ret[1] = 0
    return nothing
end

function bisearch_ord!(v::Array{Int64,1}, target::Int64,
                       ret::Array{Int64,1})
    ret[3] = length(v); ret[2] = 1
    @inbounds while ret[2] <= ret[3]
        ret[1] = div(ret[2]+ret[3],2)
        if v[ret[1]] > target
            ret[3] = ret[1]-1
        elseif v[ret[1]] < target
            ret[2] = ret[1]+1
        else
            return nothing
        end
    end
    ret[1] = 0
    return nothing
end

function func_j(j1::I,j2::I,m1::I,m2::I) where{I<:Int64}
    sqrt( (0.5*j1*(0.5*j1+1.0)-0.5*m1*(0.5*m1-1.0))*(0.5*j2*(0.5*j2+1.0)-0.5*m2*(0.5*m2+1.0)) )
end

function J_from_JJ1(JJ,tol=1.e-4)
    for J = 0:100 ## ad hoc J<=50
        hJ = 0.5*J
        if abs(hJ*(hJ+1.0)-JJ) <tol
            return hJ
        end
    end
    return -1.0
end

function prep_pp(mstates_p::Array{Array{Int64,1},1},
                 pbits::Array{Array{Int64,1}},
                 bVpp::Array{bit2b,1},
                 Vpp::Array{Float64,1}) where {IA<:Array{Int64,1}}
    lMp=length(pbits)
    lmstates_p = length(mstates_p)
    ppinfo = [ [ T1info[] for j=1:length(pbits[i]) ] for i = 1:lMp]

    @inbounds @threads for pblock_i = 1:lMp # number of proton subblock
        TF=[true]; ret=[1,-1]; ridx=[-1,-1,-1]
        @inbounds for (Npi,Phi) in enumerate(pbits[pblock_i])
            @inbounds for nth = 1:length(Vpp)
                if Vpp[nth] == 0.0; continue;end
                TF_connectable(Phi,bVpp[nth],TF)
                if TF[1]==false; continue;end
                calc_phase!(Phi,bVpp[nth],lmstates_p,ret)#ret=>phase,fpbit
                bisearch!(pbits[pblock_i],ret[2],ridx)#ridx=Npf
                if Npi <= ridx[1]
                    coef = Vpp[nth]*ret[1]*ifelse(Npi==ridx[1],0.5,1.0)
                    push!(ppinfo[pblock_i][Npi], T1info(ridx[1],coef))
                end
            end
        end
    end
    return ppinfo
end

function prep_nn(mstates_n::Array{Array{Int64,1},1},
                 nbits::Array{Array{Int64,1}},
                 bVnn::Array{bit2b,1},
                 Vnn::Array{Float64,1}) where {IA<:Array{Int64,1}}
    lMn=length(nbits)
    lmstates_n = length(mstates_n)
    nninfo = [ [ T1info[ ] for j=1:length(nbits[i]) ] for i = 1:lMn]
    @inbounds @threads for nblock_i = 1:lMn
        TF=[true]; ret=[0,-1];ridx=[-1,-1,-1]
        @inbounds for (Nni,Phi) in enumerate(nbits[nblock_i])
            @inbounds for nth = 1:length(Vnn)
                if Vnn[nth] == 0.0; continue;end
                TF_connectable(Phi,bVnn[nth],TF)
                if TF[1]==false; continue;end
                calc_phase!(Phi,bVnn[nth],lmstates_n,ret)
                bisearch!(nbits[nblock_i],ret[2],ridx)
                if Nni <= ridx[1]
                    coef = Vnn[nth]*ret[1]*ifelse(Nni==ridx[1],0.5,1.0)
                    push!(nninfo[nblock_i][Nni],T1info(ridx[1],coef))
                end
            end
        end
    end
    return nninfo
end
function prep_pn(lblock::Int64,tdims::IA,
                 l_pbit::Int64,l_nbit::Int64,
                 pbits::Array{Array{Int64,1}},
                 nbits::Array{Array{Int64,1}},
                 Mps::IA,delMs::IA,
                 bVpn::Array{Array{Array{Int64,1},1},1},
                 Vpn::Array{Array{Float64,1},1}) where {IA<:Array{Int64,1}}
    #to = TimerOutput()
    maxDeltaM = maximum(delMs)
    bfs = [ Int64[ ] for i = 1:lblock]
    bis = [ Int64[ ] for i = 1:lblock]
    for i = 1:lblock
        for j = i:lblock
            if Mps[j] - Mps[i] in delMs
                push!(bfs[i],j); push!(bis[j],i)
            end
        end
    end

    p_NiNfs = [ [ [ ifph[] ] for j=1:length(bfs[i]) ] for i=1:lblock]
    n_NiNfs = [ [ [ ifph[] ] for j=1:length(bfs[i]) ] for i=1:lblock]
    hits = [ [0,0] for i=1:lblock]
    @threads for bi = 1:lblock
        ret = [0,0,0]
        for j = 1:length(bfs[bi])
            bf = bfs[bi][j]
            deltaM = Mps[bf] - Mps[bi]
            if (deltaM in delMs) == false;continue;end
            bisearch!(delMs,deltaM,ret) # Vidx=ret[1]
            hits[bi] += calc_1b_jumps!(bi,bf,j,lblock,tdims,
                                   l_pbit,l_nbit,pbits,nbits,
                                   bVpn[ret[1]],Vpn[ret[1]],
                                   p_NiNfs,n_NiNfs)
        end
    end
    return bis,bfs,p_NiNfs,n_NiNfs, hits
end

function print_vec(s,v;ine=false)
    s *= " "
    for i = 1:length(v)
        if ine
            s *= @sprintf "%9.1e" v[i]
        else
            s *= @sprintf "%10.4f" v[i]
            #s *= @sprintf "%25.15f" v[i]
        end
    end
    println(s)
end

function calc_1b_jumps!(bi::Int64,bf::Int64,j::Int64,
                        lblock::Int64,tdims::IA,
                        l_pbit::Int64,l_nbit::Int64,
                        pbits::AAI,nbits::AAI,
                        bVpn::AAI,Vpn::Array{Float64,1},
                        p_NiNfs::Aifph,
                        n_NiNfs::Aifph) where {IA<:Array{Int64,1},
                                               AAI<:Array{Array{Int64,1},1},
                                               Aifph<:Array{Array{Array{Array{ifph,1},1},1},1}}
    TF=[true]; ret=[1,-1] # [phase,possible]
    ridx=[-1,-1,-1]
    plist =  [ ifph[ ] for i =1:length(Vpn) ]
    nlist =  [ ifph[ ] for i =1:length(Vpn) ]
    hits = [0,0]
    @inbounds for nth = 1:length(Vpn)
        a,b,c,d = bVpn[nth]; if Vpn[nth] == 0.0; continue;end
        @inbounds for (Npi,pPhi) in enumerate(pbits[bi])# Npi:idx pPhi:bit pSD
            TF_connectable_1(pPhi,a,c,TF)
            if TF[1]==false; continue;end
            calc_phase_1!(pPhi,a,c,ret)# ret<=[phase,bit]
            bisearch!(pbits[bf],ret[2],ridx)# ridx=[Npf]
            if ridx[1] ==0 ;continue;end
            push!(plist[nth],ifph(Npi,ridx[1],ret[1]==-1))
            hits[1] += 1
        end
        @inbounds for (Nni,nPhi) in enumerate(nbits[bi]) # Nni:idx
            TF_connectable_1(nPhi,b,d,TF);if TF[1]==false; continue;end
            calc_phase_1!(nPhi,b,d,ret)
            bisearch!(nbits[bf],ret[2],ridx) # ridx=[Nnf]
            if ridx[1] == 0;continue;end
            push!(nlist[nth],ifph(Nni,ridx[1],ret[1]==-1))
            hits[2] += 1
        end
    end
    p_NiNfs[bi][j] = plist
    n_NiNfs[bi][j] = nlist
    return hits
end

function add_bl_T!(q::Int64,k::Int64,
                   Tmat::FA2,R::FA2) where{FA2<:Array{Float64,2}}
    Tmat[q*k+1:q*k+q,q*k-q+1:q*k] .= R
    Tmat[q*k-q+1:q*k,q*k+1:q*k+q] .= R'
    return nothing
end

function vec_to_block(vecs,V)
    for (b,v) in enumerate(vecs)
        V[:,b] .= v
    end
    return nothing
end

function writeRitzvecs(mdim,mtot,vals,totJs,Rvecs,oupf)
    num_ev = length(vals)
    io = open(oupf,"w")
    write(io,Int32(num_ev))
    write(io,Int32(mtot))
    for i = 1:num_ev
        write(io,vals[i])
    end
    js= [ Int32(2*totJs[i]) for i=1:num_ev]
    for i = 1:num_ev
        write(io,Int32(js[i]))
    end
    for nth = 1:num_ev; write(io,Rvecs[nth]); end
    close(io)
end
"""
    entropy(Rvecs,pbits,nbits,tdims,to)

To calculate entanglement entropy (not optimized yet).

```math
S  = -\\mathrm{Tr} \\rho \\log (\\rho)
```
"""
function entropy(Rvecs,pbits,nbits,tdims,to)
    pdims = [0];for bi=1:length(pbits);push!(pdims,pdims[end]+length(pbits[bi]));end
    ndims = [0];for bi=1:length(nbits);push!(ndims,ndims[end]+length(nbits[bi]));end
    pdim = pdims[end]; ndim = ndims[end]
    rho_p = zeros(Float64,pdim,pdim)
    rho_n = zeros(Float64,ndim,ndim)
    for (n,Rvec) in enumerate(Rvecs)
        rho_p .= 0.0#; rho_n .= 0.0
        for bi = 1:length(pbits)
            idim = tdims[bi]
            nbit = nbits[bi]
            ln = length(nbit)
            pofst = pdims[bi]
            for (pid_i,pbit_i) in enumerate(pbits[bi])
                tMi = idim + (pid_i-1)*ln
                ri = pofst + pid_i
                for (pid_j,pbit_j) in enumerate(pbits[bi])
                    if pid_j <pid_i;continue;end
                    tMj = idim + (pid_j-1)*ln
                    rj = pofst + pid_j
                    for nidx =1:ln
                        Mi = tMi+nidx; Mj = tMj+nidx
                        tmp = Rvec[Mi] .* Rvec[Mj]
                        rho_p[ri,rj] += tmp
                        if ri != rj;rho_p[rj,ri] += tmp;end
                    end
                end
            end
        end
        s = -tr( rho_p * log(rho_p) )
        println("n= $n entropy $s")
    end
    return nothing
end

function make_Hermite!(H,Dim)
    @inbounds @simd for i = 1:Dim
        for j=i+1:Dim
            tmp = 0.5 * (H[i,j] + H[j,i])
            H[i,j] = tmp; H[j,i]=tmp
        end
    end
    return nothing
end

function make_distribute(num_task)
    lblock = length(num_task)
    n = nthreads()
    r = div(lblock,n)
    tasks = [ num_task[bi][1]*num_task[bi][2] for bi=1:lblock]
    idxs = sortperm(tasks,rev=true)
    tmp = [ Int64[ ] for i=1:n]
    imin = 1;order = true

    while length([(tmp...)...]) < lblock
        imax = minimum([imin+n-1,lblock])
        if order == true
            for i = imin:imax
                push!(tmp[i-imin+1],idxs[i])
            end
            order = false
        else
            #for i = imax:-1:imin
            for i = n:-1:1
                j = imax -(n-i)
                if j < imin
                    push!(tmp[i],0)
                else
                    push!(tmp[i],idxs[j])
                end
            end
            order = true
        end
        imin += n
    end
    task_idxs =[(tmp...)...]
    # println("idxs $idxs")
    # println("tmp $tmp")
    # println("task_idxs $task_idxs")
    # println("uniq ", unique(task_idxs), " lblock: $lblock")
    # am = zeros(Float64,n)
    # for i = 1:n
    #     for j = 1:length(tmp[i])
    #         if tmp[i][j] != 0
    #             am[i] += tasks[tmp[i][j]]
    #         end
    #     end
    # end
    # println("task amount ", log10.(am))
    return task_idxs
end

function make_distribute_J(Jtasks)
    lblock = length(Jtasks)
    n = nthreads()
    r = div(lblock,n)
    idxs = sortperm(Jtasks,rev=true)
    tmp = [ Int64[ ] for i=1:n]
    imin = 1;order = true

    while length([(tmp...)...]) < lblock
        imax = minimum([imin+n-1,lblock])
        if order == true
            for i = imin:imax
                push!(tmp[i-imin+1],idxs[i])
            end
            order = false
        else
            for i = n:-1:1
                j = imax -(n-i)
                if j < imin
                    push!(tmp[i],0)
                else
                    push!(tmp[i],idxs[j])
                end
            end
            order = true
        end
        imin += n
    end
    Jidxs = [(tmp...)...]
    return Jidxs
end

function idx_from_ij(i::I,j::I,Dim::I) where{I<:Int64}
    return Int64((Dim+1)*(i-1) -(i*(i-1))/2+j-i+1)
end

function sparse_mat_plot(Mat::Array{Float64,2},ba="")
    n1,n2 = size(Mat)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    for j = 1:n2
        for i=1:n1
            if Mat[i,j] != 0.0
                ax.scatter(i,j,marker=",",color="green", alpha=0.6)
            end
        end
    end
    plt.savefig("./mat_plot_"*ba*".pdf")
    plt.close()
end

function mdim_print(target_nuc,Z,N,cp,cn,vp,vn,mdim,tJ=-1)
    s0 = @sprintf "%6s %14s %10s %10s" target_nuc "Z,N=($Z,$N)" "c($cp,$cn)" "v($vp,$vn)"
    if tJ != -1; s0 * " 2J=$tJ ";end
    s1 = @sprintf "%s %12i %s %5.2f %s " "  mdim:" mdim "(" log10(mdim) ")"
    println(s0,s1)
end

function JT1(bi,oPP,oNN,
             p_sps,n_sps,
             mstates_p,mstates_n,
             pbits,nbits,
             tdims)
    pMat = zeros(Float64,length(pbits[bi]),length(pbits[bi]))
    nMat = zeros(Float64,length(nbits[bi]),length(nbits[bi]))
    lp = length(mstates_p);ln = length(mstates_n)
    TF=[true];ridx=[-1,-1,-1]
    ret_a = [-1,-1];ret_b = [-1,-1]
    ret_c = [-1,-1];ret_d = [-1,-1]
    idim = tdims[bi]
    l_Nn = length(nbits[bi])
    #### pp
    vec_ani = [false for i = 1:lp]; vec_cre=[false for i = 1:lp]
    for i2 = 1:length(p_sps)
        j2 = p_sps[i2][3]
        for m2 = -j2:2:j2-2
            vec_ani .= false; vec_cre .= false
            @inbounds for k = 1:lp-1
                if @views mstates_p[k][1:4] != p_sps[i2];continue;end
                if mstates_p[k][5] == m2
                    vec_ani[k] = true; vec_cre[k+1] = true
                    break
                end
            end
            bitarr_to_int!(vec_ani,ret_d); bitarr_to_int!(vec_cre,ret_c)
            for (SDidx,Phi) in enumerate(pbits[bi])
                TF_connectable_1(Phi,ret_c[1],ret_d[1],TF)
                if TF[1]==false;continue;end
                APhi = ret_c[1] | ((~ret_d[1]) & Phi)
                for i1 = 1:length(p_sps)
                    j1 = p_sps[i1][3]
                    @inbounds for m1 = -j1+2:2:j1
                        vec_ani .= false; vec_cre .= false
                        @inbounds for k = 2:lp
                            if @views mstates_p[k][1:4] != p_sps[i1];continue;end
                            if mstates_p[k][5] == m1
                                vec_ani[k] = true; vec_cre[k-1] = true
                                break
                            end
                        end
                        bitarr_to_int!(vec_ani,ret_b)
                        bitarr_to_int!(vec_cre,ret_a)
                        TF_connectable_1(APhi,ret_a[1],ret_b[1],TF)
                        if TF[1]==false;continue;end
                        NPhi = ret_a[1] | ((~ret_b[1]) & APhi)
                        bisearch!(pbits[bi],NPhi,ridx) #ridx <=Npf
                        Npi = SDidx
                        Npf = ridx[1]
                        if Npi > Npf; continue;end
                        fac = ifelse(Npi==Npf,0.5,1.0) .* func_j(j1,j2,m1,m2)
                        pMat[Npi,Npf] += fac
                    end
                end
            end
        end
    end
    #### nn
    vec_ani = [false for i = 1:ln];vec_cre = [false for i = 1:ln]
    vec_ani_2 = [false for i = 1:ln];vec_cre_2 = [false for i = 1:ln]
    for i2 = 1:length(n_sps)
        j2 = n_sps[i2][3]
        for m2 = -j2:2:j2-2
            vec_ani_2 .= false ; vec_cre_2 .= false
            @inbounds for k = 1:ln-1
                if @views mstates_n[k][1:4] != n_sps[i2];continue;end
                if mstates_n[k][5] == m2
                    vec_ani_2[k] = true; vec_cre_2[k+1] = true;break
                end
            end
            bitarr_to_int!(vec_ani_2,ret_d); bitarr_to_int!(vec_cre_2,ret_c)
            for (SDidx,Phi) in enumerate(nbits[bi])
                TF_connectable_1(Phi,ret_c[1],ret_d[1],TF)
                if TF[1]==false;continue;end
                APhi = ret_c[1] | ((~ret_d[1]) & Phi)
                @inbounds for i1 = 1:length(n_sps)
                    j1 = n_sps[i1][3]
                    @inbounds for m1 = -j1+2:2:j1
                        vec_ani .= false; vec_cre .= false
                        @inbounds for k = 2:ln
                            if @views mstates_n[k][1:4] != n_sps[i1];continue;end
                            if mstates_n[k][5] == m1
                                vec_ani[k] = true; vec_cre[k-1]=true;break
                            end
                        end
                        bitarr_to_int!(vec_ani,ret_b)
                        bitarr_to_int!(vec_cre,ret_a)
                        TF_connectable_1(APhi,ret_a[1],ret_b[1],TF)
                        if TF[1]==false;continue;end
                        NPhi = ret_a[1] | ((~ret_b[1]) & APhi)
                        bisearch!(nbits[bi],NPhi,ridx) #ridx <=Npf
                        Nni = SDidx
                        Nnf = ridx[1]
                        if Nni > Nnf; continue;end
                        fac = ifelse(Nni==Nnf,0.5,1.0) .* func_j(j1,j2,m1,m2)
                        nMat[Nni,Nnf] += fac
                    end
                end
            end
        end
    end
    op_p=MiMf[]; op_n=MiMf[]
    m = size(pMat)[1]
    for j=1:m
        for i=1:j
            if pMat[i,j] != 0.0; push!(op_p,MiMf(i,j,pMat[i,j]));end
        end
    end
    m = size(nMat)[1]
    for j=1:m
        for i=1:j
            if nMat[i,j] != 0.0; push!(op_n,MiMf(i,j,nMat[i,j]));end
        end
    end
    oPP[bi] = op_p; oNN[bi] = op_n
    return nothing
end

function prep_J(tdims,p_sps,n_sps,
                mstates_p::Array{Array{Int64,1},1},
                mstates_n::Array{Array{Int64,1},1},
                pbits::Array{Array{Int64,1},1},
                nbits::Array{Array{Int64,1},1})
    #to = TimerOutput()
    lp=length(mstates_p); ln=length(mstates_n)
    lblock = length(nbits)
    oPP = [ MiMf[] for i =1:lblock]
    oNN = [ MiMf[] for i =1:lblock]
    #@timeit to "nn//pp" begin
    @inbounds @threads for bi=1:lblock
        JT1(bi,oPP,oNN,p_sps,n_sps,
            mstates_p,mstates_n,pbits,nbits,tdims)
    end
    oPNu = [ Jpninfo[]  for i =1:lblock]
    oPNd = [ Jpninfo[]  for i =1:lblock]
    #end
    #@timeit to "pn" begin
    vec_p_ani = [false for i = 1:lp]; vec_p_cre = [false for i = 1:lp]
    vec_n_ani = [false for i = 1:ln]; vec_n_cre = [false for i = 1:ln]
    ret_a = [-1,-1];ret_b = [-1,-1]
    ret_c = [-1,-1];ret_d = [-1,-1]
    birange = [1:lblock-1,2:lblock]
    for pidx = 1:length(p_sps)
        jp = p_sps[pidx][3]
        mp_range = [-jp:2:jp-2,-jp+2:2:jp]
        for nidx = 1:length(n_sps)
            jn = n_sps[nidx][3]
            mn_range = [-jn+2:2:jn,-jn:2:jn-2]
            for ud = 1:2 # up:1 down:2
                @inbounds for mp in mp_range[ud]
                    @inbounds for mn in mn_range[ud]
                        fac = 1.0
                        if ud == 1
                            fac = func_j(jn,jp,mn,mp)
                        else
                            fac = func_j(jp,jn,mp,mn)
                        end
                        if fac==0.0;continue;end
                        vec_p_ani .= false; vec_p_cre .= false
                        vec_n_ani .= false; vec_n_cre .= false

                        if ud == 2  ### for pn down
                            @inbounds for k = 2:lp
                                if mstates_p[k][1] != p_sps[pidx][1];continue;end
                                if mstates_p[k][2] != p_sps[pidx][2];continue;end
                                if mstates_p[k][3] != p_sps[pidx][3];continue;end
                                if mstates_p[k][4] != p_sps[pidx][4];continue;end
                                if mstates_p[k][5] == mp
                                    vec_p_ani[k] = true; vec_p_cre[k-1] = true
                                    break
                                end
                            end
                            @inbounds for k = 1:ln-1
                                if mstates_n[k][1] != n_sps[nidx][1];continue;end
                                if mstates_n[k][2] != n_sps[nidx][2];continue;end
                                if mstates_n[k][3] != n_sps[nidx][3];continue;end
                                if mstates_n[k][4] != n_sps[nidx][4];continue;end
                                if mstates_n[k][5] == mn
                                    vec_n_ani[k] = true; vec_n_cre[k+1] = true
                                    break
                                end
                            end
                        else  ### for p up n down
                            @inbounds for k = 1:lp-1
                                if mstates_p[k][1] != p_sps[pidx][1];continue;end
                                if mstates_p[k][2] != p_sps[pidx][2];continue;end
                                if mstates_p[k][3] != p_sps[pidx][3];continue;end
                                if mstates_p[k][4] != p_sps[pidx][4];continue;end
                                if mstates_p[k][5] == mp
                                    vec_p_ani[k] = true;vec_p_cre[k+1] = true;break
                                end
                            end
                            @inbounds for k = 2:ln
                                if mstates_n[k][1] != n_sps[nidx][1];continue;end
                                if mstates_n[k][2] != n_sps[nidx][2];continue;end
                                if mstates_n[k][3] != n_sps[nidx][3];continue;end
                                if mstates_n[k][4] != n_sps[nidx][4];continue;end
                                if mstates_n[k][5] == mn
                                    vec_n_ani[k] = true;vec_n_cre[k-1] = true;break
                                end
                            end
                        end
                        bitarr_to_int!(vec_p_ani,ret_c)
                        bitarr_to_int!(vec_p_cre,ret_a)
                        bitarr_to_int!(vec_n_ani,ret_d)
                        bitarr_to_int!(vec_n_cre,ret_b)

                        @inbounds @threads for bi in birange[ud]
                            TF=[true]; ret_p=[0,0];ret_n=[0,0]
                            ridx_p=[-1,-1,-1];ridx_n=[-1,-1,-1]
                            bf = bi + Int64(2*(1.5-ud))
                            pjump = ifph[]; njump=ifph[]
                            for (Npi,pPhi) in enumerate(pbits[bi])
                                TF_connectable_1(pPhi,ret_a[1],ret_c[1],TF)
                                if TF[1]==false;continue;end
                                calc_phase_1!(pPhi,ret_a[1],ret_c[1],ret_p)
                                bisearch!(pbits[bf],ret_p[2],ridx_p)
                                phase_p = ret_p[1]
                                Npf= ridx_p[1]
                                push!(pjump,ifph(Npi,Npf,phase_p==-1))
                            end
                            for (Nni,nPhi) in enumerate(nbits[bi])
                                TF_connectable_1(nPhi,ret_b[1],ret_d[1],TF)
                                if TF[1]==false;continue;end
                                calc_phase_1!(nPhi,ret_b[1],ret_d[1],ret_n)
                                bisearch!(nbits[bf],ret_n[2],ridx_n)
                                phase_n = ret_n[1]
                                Nnf= ridx_n[1]
                                push!(njump,ifph(Nni,Nnf,phase_n==-1))
                            end
                            if length(pjump)*length(njump) == 0; continue;end
                            if ud == 1
                                push!(oPNu[bi],Jpninfo(fac,pjump,njump))
                            else
                                push!(oPNd[bi],Jpninfo(fac,pjump,njump))
                            end
                        end
                    end
                end
            end
        end
    end
    #end
    #show(to, allocations = true,compact = false)
    #println("")
    return oPP,oNN,oPNu,oPNd
end

function operate_J!(Rvec,Jv,pbits,nbits,tdims,Jidxs,
                    oPP,oNN,oPNu,oPNd,beta_J=1.0)
    #to = TimerOutput()
    lblock=length(pbits)
    #@timeit to "pp/nn " begin
    @inbounds @threads for bi in Jidxs
        if bi==0;continue;end
        idim = tdims[bi]
        lNn = length(nbits[bi])
        lNp = length(pbits[bi])
        opPP = oPP[bi]
        opNN = oNN[bi]
        offset = idim-lNn
        @inbounds for tmp in opPP
            Npi =tmp.Mi; Npf=tmp.Mf; fac=tmp.fac .* beta_J
            tMi = offset + Npi*lNn
            tMf = offset + Npf*lNn
            @inbounds @simd for nidx = 1:lNn
                Mi = tMi+nidx; Mf = tMf+nidx
                Jv[Mf] += fac .* Rvec[Mi]
                Jv[Mi] += fac .* Rvec[Mf]
            end
        end
        @inbounds for tmp in opNN #nn
            Nni =tmp.Mi; Nnf=tmp.Mf; fac=tmp.fac .* beta_J
            tMi = offset + Nni
            tMf = offset + Nnf
            @inbounds @simd for pidx = 1:lNp
                Mi = tMi+pidx*lNn; Mf = tMf+pidx*lNn
                Jv[Mf] += fac .* Rvec[Mi]
                Jv[Mi] += fac .* Rvec[Mf]
            end
        end
    end
    #end
    #@timeit to "pn down" begin
    @inbounds @threads for bi = 2:lblock
        operator = oPNd[bi]
        bf = bi-1
        l_Nn_i = length(nbits[bi])
        l_Nn_f = length(nbits[bf])
        off_i = tdims[bi] - l_Nn_i
        off_f = tdims[bf] - l_Nn_f
        @inbounds for top in operator
            pj = top.pjump
            nj = top.njump
            fac =top.fac .* beta_J
            @inbounds for tmp_p in pj
                phase_p=tmp_p.phase
                tMi = off_i + tmp_p.i * l_Nn_i
                tMf = off_f + tmp_p.f * l_Nn_f
                @inbounds for tmp_n in nj
                    phase_n=tmp_n.phase
                    Mi = tMi + tmp_n.i
                    Mf = tMf + tmp_n.f
                    Jv[Mf] += ifelse(phase_p!=phase_n,-fac,fac) .* Rvec[Mi]
                end
            end
        end
    end
    @inbounds @threads for bi = 1:lblock-1
        operator = oPNu[bi]
        bf = bi+1
        l_Nn_i = length(nbits[bi])
        l_Nn_f = length(nbits[bf])
        off_i = tdims[bi] - l_Nn_i
        off_f = tdims[bf] - l_Nn_f
        @inbounds for top in operator
            pj  = top.pjump
            nj  = top.njump
            fac = top.fac .* beta_J
            @inbounds for tmp_p in pj
                phase_p=tmp_p.phase
                tMi = off_i + tmp_p.i * l_Nn_i
                tMf = off_f + tmp_p.f * l_Nn_f
                @inbounds for tmp_n in nj
                    Nni = tmp_n.i; Nnf = tmp_n.f; phase_n=tmp_n.phase
                    Mi = tMi + tmp_n.i
                    Mf = tMf + tmp_n.f
                    Jv[Mf] += ifelse(phase_p!=phase_n,-fac,fac) .* Rvec[Mi]
                end
            end
        end
    end
    #end
    #show(to, allocations = true,compact = false)
    #println("")
    return nothing
end

"""
    samplerun()

Sample scripts to calculate
- (a) 10 lowest states of 28Si in sd shell
- (b) 10 lowest states of 28Si with J=0
- (c) EC estimates of 10 lowest J=0 states of 28Si
- (d) (b) with the preprocessing

!!! note
    To run ```samplerun()```, you need a copy of ShellModel.jl in your environment.
    Try e.g.,
    ```sh
    git clone https://github.com/SotaYoshida/ShellModel.jl
    ```
    and then, execute `samplerun()` in the repository:
    ```julia
    >julia using ShellModel
    >julia samplerun()
    ```
    or
    ```sh
    julia -t 12 sample_run.jl
    ```
"""
function samplerun()
    ### for compilation
    println("For JIT compilation...")
    @time main("./snts/x_mass.snt","Be8",5,[0])
    @time main("./snts/x_mass.snt","Be8",5,[0];is_block=true,q=3)
    println("Done: JIT compilation")
    ### input
    tJ = 0 # target J
    tJs =[ tJ ]
    num_ev = 4; n_block = 4
    num_ev_target = 10
    sntf = "./snts/usdb.snt"; target_nuc = "Si28"
    ###

    println("\n### sample (a). 10-lowest states for $target_nuc")
    @time main(sntf,target_nuc,10,[])

    println("\n### sample (b). 10-lowest states with J=$tJs")
    @time main(sntf,target_nuc,10,tJs)

    println("\n### sample (c). calc. EC estimates for $num_ev_target J=$tJs states ")
    Hs = ["./snts/usdb.snt"]
    tJNs = [ [tJ,num_ev_target] ]
    solveEC(Hs,target_nuc,tJNs;
            verbose=true, exact_logf = "runlog_sigma1_valid.dat")

    println("\n### sample (d).  (b). with the preprocessing")
    println("Block Lanczos method with n_block=$n_block")
    @time main(sntf,target_nuc,num_ev,tJs;is_block=true,q=n_block)
    println("\n+ preprocessing")
    @time main(sntf,target_nuc,num_ev,tJs;is_block=true,q=n_block,
               in_wf="./appwavs/"*target_nuc*"_usdb_j"*string(tJ)*".appwav")
end


text:
using ThreadPools

function heavy_task(n_iter::Float64)
    x = 1.0
    for _ in 1:n_iter
        x += rand()
    end
end

function mytask()
    println(Threads.threadid())
    heavy_task(1e10)
end

function mytask(i::Int64, j::Int64)
    n_max = 1e10
    n_iter = floor(rand() * n_max)
    println("$(Threads.threadid()) on $i $j for $n_iter")
    heavy_task(n_iter)
end

n = Threads.nthreads()

for k in 1:2 # Say we do two iterations
    println("Iteration $k")
    n_enco = 2
    n_cont = 2
    indexes = [(i, j) for i in 1:n_enco for j in 1:n_cont]
    Threads.@threads for l in 1:(n_enco + n_cont)
        i, j = indexes[l]
        mytask(i, j)
    end
end

ext: jl


text:
using DSGE, Distributions, HDF5, MAT, ModelConstructors, NLsolve, Nullables, OrderedCollections, Random
import DSGE: eqcond, measurement, init_observable_mappings!, init_pseudo_observable_mappings!,
    steadystate!, init_parameters!, init_subspec!, augment_states, model_settings!, description

# src directory
include("util.jl")

# julia_models/THANK
include("julia_models/COVIDHANK/covidhank.jl")
include("julia_models/COVIDHANK/eqcond.jl")
include("julia_models/COVIDHANK/measurement.jl")
include("julia_models/COVIDHANK/observables.jl")
include("julia_models/COVIDHANK/subspecs.jl")
include("julia_models/COVIDHANK/augment_states.jl")


text:
using Random
using ProximalOperators
using Test

@testset "NormL1plusL2 standard case" begin
    f = NormL1(1.0)
    g = NormL2(2.0)
    fplusg = NormL1plusL2(1.0, 2.0)

    x = randn(50)

    y1, f1 = prox(f, x)
    y2, f2 = prox(g, y1)

    y3, f3 = prox(fplusg, x)

    @test f3 ≈ f(y3)+g(y3)
    @test y3 ≈ y2
end

@testset "NormL1plusL2 norm constructor" begin
    f = NormL1(1.0)
    g = NormL2(2.0)
    fplusg = NormL1plusL2(f, g)

    x = randn(50)

    y1, f1 = prox(f, x)
    y2, f2 = prox(g, y1)

    y3, f3 = prox(fplusg, x)

    @test f3 ≈ f(y3)+g(y3)
    @test y3 ≈ y2
end

@testset "NormL1plusL2 vector case" begin
    λ1 = abs.(randn(50))

    f = NormL1(λ1)
    g = NormL2(2.0)

    fplusg = NormL1plusL2(λ1, 2.0)

    x = randn(50)

    y1, f1 = prox(f, x)
    y2, f2 = prox(g, y1)

    y3, f3 = prox(fplusg, x)

    @test f3 ≈ f(y3)+g(y3)
    @test y3 ≈ y2
end


text:
using FractalFlow
using Base.Test

include("test_curve.jl")


text:
# function tcm_estimateMMNL_kcv(N, K_list, nfolds, transaction_folds, transactions, numInitializations, time_limit, perTol)

# 	unique_fold_ids = unique(transaction_folds)

# 	KL_avg_holdout_fold = zeros(length(K_list), nfolds)

# 	start_time = time()

# 	for fold_ind in 1:nfolds
# 		fold = unique_fold_ids[fold_ind]

# 		@show fold

# 		@show unique_fold_ids

# 		println("Pre-subsetting ")

# 		#@show transaction_folds
# 		#temp = find(transaction_folds .!= fold)
# 		#@show temp
# 		@show size(transaction_folds)
# 		@show size(transactions)
# 		@show fold

# 		temp = find(transaction_folds .!= fold)
# 		@show temp[1:5]
# 		trainfolds_transactions = transactions[ find(transaction_folds .!= fold), :]
# 		println("Post-subsetting")
# 		testfold_transactions = transactions[ find(transaction_folds .== fold), :]

# 		trainfolds_assortments, trainfolds_transaction_counts = tcm_transactionsToCounts(N, trainfolds_transactions)
# 		testfold_assortments, testfold_transaction_counts = tcm_transactionsToCounts(N, testfold_transactions)

# 		for Kind in 1:length(K_list)
# 			K = K_list[Kind]

# 			u, class_probs, loglik, grand_elapsed_time = tcm_estimateMMNL(N, K, trainfolds_assortments, trainfolds_transaction_counts, numInitializations, time_limit, perTol)

# 			testfold_MMNL_predict, testfold_MMNL_predict_by_class = tcm_predictMMNL(N, K, class_probs, u, testfold_assortments)

# 			testfold_KL, testfold_noinf_count = tcm_evaluateKL(N, testfold_MMNL_predict, testfold_transaction_counts)

# 			KL_avg_holdout_fold[Kind, fold_ind] = testfold_KL / testfold_noinf_count # sum(testfold_transaction_counts)
# 		end
# 	end

# 	end_time = time()
# 	elapsed_time = end_time - start_time

# 	KL_avg_by_K = mean(KL_avg_holdout_fold, 2)

# 	best_Kind = indmin(KL_avg_by_K)
# 	best_K = K_list[best_Kind]
# 	best_KL_avg_cv = KL_avg_by_K[best_Kind]

# 	return best_K, best_KL_avg_cv, elapsed_time
# end

function tcm_estimateMMNL_kcv_weighted(N, K_list, nfolds, transaction_folds, transactions, numInitializations, time_limit, perTol)

	unique_fold_ids = unique(transaction_folds)

	KL_total_holdout_fold = zeros(length(K_list), nfolds)
	noinf_count_holdout_fold = zeros(length(K_list), nfolds)

	start_time = time()

	for fold_ind in 1:nfolds
		fold = unique_fold_ids[fold_ind]

		@show fold

		@show unique_fold_ids

		println("Pre-subsetting ")

		#@show transaction_folds
		#temp = find(transaction_folds .!= fold)
		#@show temp
		@show size(transaction_folds)
		@show size(transactions)
		@show fold

		temp = find(transaction_folds .!= fold)
		@show temp[1:5]
		trainfolds_transactions = transactions[ find(transaction_folds .!= fold), :]
		println("Post-subsetting")
		testfold_transactions = transactions[ find(transaction_folds .== fold), :]

		trainfolds_assortments, trainfolds_transaction_counts = tcm_transactionsToCounts(N, trainfolds_transactions)
		testfold_assortments, testfold_transaction_counts = tcm_transactionsToCounts(N, testfold_transactions)

		for Kind in 1:length(K_list)
			K = K_list[Kind]

			u, class_probs, loglik, grand_elapsed_time = tcm_estimateMMNL(N, K, trainfolds_assortments, trainfolds_transaction_counts, numInitializations, time_limit, perTol)

			testfold_MMNL_predict, testfold_MMNL_predict_by_class = tcm_predictMMNL(N, K, class_probs, u, testfold_assortments)

			testfold_KL, testfold_noinf_count, testfold_KL_avg = tcm_evaluateKL(N, testfold_MMNL_predict, testfold_transaction_counts)

			@show testfold_KL
			@show testfold_noinf_count
			@show testfold_KL_avg

			KL_total_holdout_fold[Kind, fold_ind] = testfold_KL #testfold_KL / sum(testfold_transaction_counts)
			noinf_count_holdout_fold[Kind, fold_ind] = testfold_noinf_count

			@show KL_total_holdout_fold
			@show noinf_count_holdout_fold
		end
	end

	end_time = time()
	elapsed_time = end_time - start_time

	KL_avg_by_K = sum(KL_total_holdout_fold, 2) ./ sum(noinf_count_holdout_fold, 2)  * 100

	best_Kind = indmin(KL_avg_by_K)
	best_K = K_list[best_Kind]
	best_KL_avg_cv = KL_avg_by_K[best_Kind]

	@show KL_avg_by_K
	@show KL_total_holdout_fold
	@show noinf_count_holdout_fold

	return best_K, best_KL_avg_cv, elapsed_time
end

function tcm_estimateMMNL_kcv_weighted_ass(N, K_list, nfolds, assortment_folds, assortments, transaction_counts, numInitializations, time_limit, perTol)

	unique_fold_ids = unique(assortment_folds)

	KL_total_holdout_fold = zeros(length(K_list), nfolds)
	noinf_count_holdout_fold = zeros(length(K_list), nfolds)

	start_time = time()

	for fold_ind in 1:nfolds
		fold = unique_fold_ids[fold_ind]

		# trainfolds_assortments, trainfolds_transaction_counts = tcm_transactionsToCounts(N, trainfolds_transactions)
		# testfold_assortments, testfold_transaction_counts = tcm_transactionsToCounts(N, testfold_transactions)

		trainfolds_assortments = assortments[ find(assortment_folds .!= fold), :]
		testfold_assortments = assortments[ find(assortment_folds .== fold), :]

		trainfolds_transaction_counts = transaction_counts[ find(assortment_folds .!= fold), :]
		testfold_transaction_counts = transaction_counts[ find(assortment_folds .== fold), :]

		for Kind in 1:length(K_list)
			K = K_list[Kind]

			u, class_probs, loglik, grand_elapsed_time = tcm_estimateMMNL(N, K, trainfolds_assortments, trainfolds_transaction_counts, numInitializations, time_limit, perTol)

			testfold_MMNL_predict, testfold_MMNL_predict_by_class = tcm_predictMMNL(N, K, class_probs, u, testfold_assortments)

			testfold_KL, testfold_noinf_count, testfold_KL_avg = tcm_evaluateKL(N, testfold_MMNL_predict, testfold_transaction_counts)

			KL_total_holdout_fold[Kind, fold_ind] = testfold_KL #testfold_KL / sum(testfold_transaction_counts)
			noinf_count_holdout_fold[Kind, fold_ind] = testfold_noinf_count
		end
	end

	end_time = time()
	elapsed_time = end_time - start_time

	KL_avg_by_K = sum(KL_total_holdout_fold, 2) ./ sum(noinf_count_holdout_fold, 2)  * 100

	best_Kind = indmin(KL_avg_by_K)
	best_K = K_list[best_Kind]
	best_KL_avg_cv = KL_avg_by_K[best_Kind]

	return best_K, best_KL_avg_cv, elapsed_time
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
using Underscores
using Utils

function parse_file(f)
    @_ eachline(f) |> map(let t = split(_); t[1], parse(Int, t[2]) end, __)
end

function load()
    open(joinpath(@__DIR__, "08_bootcode.txt"), "r") do f
        parse_file(f)
    end
end

assembly = load()

struct Infinite
    acc::Int
end

struct Finished
    acc::Int
end

function ex1(assembly)::Union{Infinite,Finished}
    valid = eachindex(assembly)
    ran = zeros(Bool, valid.stop)
    acc = 0
    pc = 1
    while pc in valid
        if ran[pc]
            return Infinite(acc)
        end
        ins, arg = assembly[pc]
        ran[pc] = true
        if ins == "acc"
            acc += arg
        elseif ins == "jmp"
            pc += arg
            continue
        end
        pc += 1
    end
    Finished(acc)
end

function ex2(assembly)
    for i in eachindex(assembly)
        before = assembly[i]
        ins, arg = before
        if ins == "jmp"
            assembly[i] = "nop", arg
            acc = ex1(assembly)
            assembly[i] = before
            if isa(acc, Finished)
                return acc
            end
        elseif ins == "nop"
            assembly[i] = "jmp", arg
            acc = ex1(assembly)
            assembly[i] = before
            if isa(acc, Finished)
                return acc
            end
        end
    end
end

println("Accumulator before infinite loop: ", ex1(assembly))
println("Accumulator after  finished  run: ", ex2(assembly))

using Test

@testset "Handy Haversacks" begin
    input = """
        nop +0
        acc +1
        jmp +4
        acc +3
        jmp -3
        acc -99
        acc +1
        jmp -4
        acc +6
        """

    example = parse_file(IOBuffer(input))

    @testset "example 1" begin
        @test ex1(example) == Infinite(5)
        @test ex2(example) == Finished(8)
    end

    @testset "results" begin
        @test ex1(assembly) == Infinite(1928)
        @test ex2(assembly) == Finished(1319)
    end
end

max_stars_repo_name: CmdQ/AoC


text:
# kuramoto
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, OrdinaryDiffEq, Random, Statistics

function plot_circle!(ax, phi_vec,omega_vec, add_colorbar=true; s = 1)
    cmap = "viridis"
    # plot circle
    alpha_vec = range(0, 2. * pi; length = 200)
    x_vec = cos.(alpha_vec)
    y_vec = sin.(alpha_vec)
    ax.plot(x_vec,y_vec, linewidth=1.0,color="k") # unit circle  # linestyle="--",

    # plot markers
    u_vec = cos.(phi_vec)
    v_vec = sin.(phi_vec)
    implt = ax.scatter(u_vec, v_vec, marker = "o" , s=70, c=omega_vec,
    cmap = cmap, zorder = 3, edgecolors = "k", linewidths = 1, vmin = -2, vmax = +2)

    # plot arrows
    D = length(u_vec)
    ax.quiver(zeros(D), zeros(D), u_vec, v_vec, omega_vec;
        # color = COLORS[5],
        cmap = cmap, # color = omega_vec,
        scale = 1.5, scale_units = "xy"
    )

    # Plot mean field
    mf = sum(exp(im*φ) for φ in phi_vec)/length(phi_vec)
    ax.quiver([0.0], [0.0], [s*real(mf)], [s*imag(mf)];
        color = COLORS[4], width = 0.02,
        scale = 1.5, scale_units = "xy"
    )

    ax.set_xticks([])
    ax.set_yticks([])
    ax.set_aspect("equal")
    if add_colorbar
        cbar = plt.colorbar(implt; ax, orientation = "horizontal", fraction=0.046, pad=0.04)
        cbar.set_ticks([-2,2])
        cbar.set_label(L"\omega_n", labelpad = -30)
    end

end

nosc1 = 50     # number of oscillators for (a), (b), (c)
nosc2 = 2000    # number of oscillators for second curve in (c)
sigma = 1.   # 0.1
T = 100.0
Ttr = 100.0
Δt = 0.1
t = 0:Δt:T
nts = length(t)

# open plot
fig = plt.figure(figsize=(figx,figy))

# plot identifiers first, otherwise they appear with colorbars
ax1 = subplot(1,3,1)
ax2 = subplot(1,3,2)
ax3 = subplot(1,3,3)

Kmin = 0.0
Kmax = 8.0   # 0.8
nK = 81
K_vec = LinRange(Kmin,Kmax,nK)
Rmean_vec = zeros(nK)
Rstd_vec = zeros(nK)

for (nj, nosc) in enumerate((nosc1, nosc2))

    rng = MersenneTwister(1234)
    omega_vec = sigma .* randn(rng,nosc)
    u0 = zeros(nosc)
    ds = Systems.kuramoto(nosc, u0; K = 0.0, ω = omega_vec)

    for iK = 1:nK
        @show iK
        K = K_vec[iK]
        set_parameter!(ds, :K, K)
        tr = trajectory(ds, T; Ttr, Δt)
        phi_vec = tr[end]

        Rts = zeros(nts)
        for j = 1:nts  # R time series
            Rts[j] = abs( sum( exp.(im .* tr[j]) ) ) ./ nosc
        end
        Rmean_vec[iK] = mean(Rts)
        Rstd_vec[iK] = std(Rts)

        if iK == 11 && nj == 1
            plot_circle!(ax1, phi_vec, omega_vec)
        end
        if iK == 81 && nj == 1  # sync
            plot_circle!(ax2, phi_vec,omega_vec)
        end

    end
    ax3.plot(K_vec,Rmean_vec, linewidth=2.0,color=COLORS[nj])
    ax3.fill_between(K_vec, Rmean_vec - Rstd_vec, Rmean_vec + Rstd_vec;
        color = COLORS[nj], alpha = 0.15
    )
end

ax3.set_ylabel("\$R\$")
ax3.set_xticks([0, 2.5, 5, 7.5])
ax3.set_xticklabels(["0", "2.5", "5", "7.5"])
ax3.set_xlabel("\$K\$", labelpad = -20)

add_identifiers!(fig, (ax1,ax2,ax3))
fig.tight_layout(pad=0.38)
wsave(plotsdir("9", "kuramoto"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
import Base: expand
IMPLEMENT_ONE_ARG_FUNC(:expand, :expand)


text:
## Aggregate batch results
using EstimateSDE
using Serialization
using Pipe: @pipe
using DataFrames

# name = "experiment01"

## Load batch settings and results
results_dir = settings_resultsdir(name)
S = deserialize(joinpath(results_dir,"settings.jld"))
df_eval = deserialize(joinpath(results_dir,"model-evaluation.jld")).df
df_data = deserialize(joinpath(results_dir,"data.jld")).df

df_eval = @pipe(
  innerjoin(df_eval,df_data,on=:index) |>
  rename(_,:label=>:kld_label,:T=>:kld_T))

## Aggregate
df_kld_mean = @pipe(
  groupby(df_eval,[:proclabel,:engine,:init_actual,:method,:pmethod,:kld_label,:kld_T]) |>
  combine(_) do df
    x = df.kld
    c = length(x)
    μ = mean(x); σ = std(x)/√c
    (kld=μ,kld_std=σ,count=c)
  end
)

## Save to file
f = joinpath(results_dir,"aggregate.jld")
R = (settings=nothing,df=df_kld_mean,status=systemstatus())
serialize(f,R)
printstyled("\nSaved to file:\n$f",color=:blue,bold=true)


text:
using Test
using Suppressor
using LogicCircuits

include("../helper/little_circuits.jl")

@testset "Queries test" begin

    r1 = fully_factorized_circuit(PlainLogicCircuit,10)

    @test isdecomposable(r1)
    @test isdecomposable(compile(PlainLogicCircuit, Lit(1)))

    @test isdeterministic(r1)
    @test isdeterministic(compile(PlainLogicCircuit, Lit(1)))

    @test !isstruct_decomposable(r1)
    @test isstruct_decomposable(compile(PlainLogicCircuit, Lit(1)))

    @test variables(r1) == BitSet(1:10)

    @test num_variables(r1) == 10
    @test issmooth(r1)

    @test isone(sat_prob(r1))
    @test model_count(r1) == BigInt(2)^10

    @test r1(true, false, true, false, true, false, true, false, true, false)
    @test r1([true, false, true, false, true, false, true, false, true, false])

    #####################
    ors = map(1:10) do v
        pos = compile(PlainLogicCircuit, var2lit(Var(v)))
        neg = compile(PlainLogicCircuit, -var2lit(Var(v)))
        pos | neg
    end
    and1 = conjoin(ors[1], ors[2], ors[1], ors[3])
    and2 = conjoin(ors[5:10])
    or1 = and1 | and2
    @test !isdecomposable(or1)
    @test !isdeterministic(or1)
    # @test !isstruct_decomposable(or1)

    #######################
    ors = map(1:10) do v
        pos = compile(PlainLogicCircuit, var2lit(Var(v)))
        neg = compile(PlainLogicCircuit, -var2lit(Var(v)))
        pos | neg
    end
    and1 = conjoin(ors[1:3])
    and2 = conjoin(ors[3:6])
    and3 = conjoin(ors[6:8])
    and4 = conjoin(ors[8:10])
    or1 = and1 | and2
    or2 = and3 | and4
    and5 = or1 & or2
    @test !isdecomposable(and5)
    @test !isdeterministic(and5)
    # @test !isstruct_decomposable(and5)

    #######################
    leaf1 = compile(PlainLogicCircuit, Lit(1))
    leaf2 = compile(PlainLogicCircuit, Lit(-1))
    and = leaf1 & leaf2
    @test !isdecomposable(and)
    @test isdeterministic(and)
    @test !isstruct_decomposable(and)

    #######################
    n0c = little_3var_constants()
    @test isdecomposable(n0c)
    @test isdeterministic(n0c)
    @test isstruct_decomposable(n0c)
    @test !isdecomposable(n0c & little_2var())
    @test !isstruct_decomposable(n0c & little_2var())
    @test !issmooth(n0c | little_2var())
    @test !isdeterministic(n0c | little_2var())
    @test issmooth(n0c)
    @test @suppress_out !iscanonical(n0c & little_3var_constants(), 5, verbose=true)

    #######################
    lits = map(1:3) do var
        compile(PlainLogicCircuit, Lit(var))
    end

    or1 = lits[1] & lits[2] | lits[1] & lits[2]
    or2 = lits[1] & lits[3] | lits[1] & lits[3]
    and1 = or1 & lits[3]
    and2 = or2 & lits[2]
    circuit = and1 | and2
    @test isstruct_decomposable(or1)
    @test isstruct_decomposable(and2)
    @test isdecomposable(circuit)
    @test !isstruct_decomposable(circuit)
    @test isstruct_decomposable(or1 & lits[3] | lits[3] & or1)
    @test isstruct_decomposable(or1 & lits[3] | or1 & lits[3])

end

@testset "Infer Vtrees Test" begin
    random = zoo_sdd("random.sdd")
    random_vtree = zoo_vtree("random.vtree")
    @test nothing == infer_vtree(random)
    @test respects_vtree(random, random_vtree)
end

@testset "Implied Literals Test" begin

    # Load and compile as a bdd
    c17 = zoo_cnf("easy/C17_mince.cnf")
    vtr = PlainVtree(num_variables(c17), :rightlinear)
    mgr = SddMgr(vtr)
    circuit = compile(mgr, c17)
    plc = PlainLogicCircuit(circuit)

    # store implied literals in each data field
    data = Dict()
    implied_literals(plc, data)

    for ornode in or_nodes(plc)
        @test num_children(ornode) == 2
        # If there's a nothing just continue, it'll always work
        if data[ornode.children[1]] === nothing ||
            data[ornode.children[2]] === nothing
            continue
        end
        implied1 = data[ornode.children[1]]
        implied2 = data[ornode.children[2]]
        neg_implied2 = BitSet(map(x -> -x, collect(implied2)))
        @test !isempty(intersect(implied1, neg_implied2))
    end
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
struct DecPolicy{P <: Policy, M <: Union{MDP, POMDP}, A} <: Policy
    policy::P # the single agent policy
    problem::M # the pomdp definition
    action_map::Vector{A}
    op # the reduction operator for utiliy fusion (e.g. sum or min)
 end

 function POMDPs.action(p::DecPolicy, b::Dict)
    vals = action_values(p, b)
    ai = findmax(vals)
    return p.action_map[ai[2]]
 end

function action_values(policy::DecPolicy, dec_belief::Dict)
    return reduce(policy.op, action_values(policy.policy, b) for (_,b) in dec_belief)
 end

 function action_values(p::AlphaVectorPolicy, b::SparseCat)

    num_vectors = length(p.alphas)
    max_values = -Inf*ones(n_actions(p.pomdp))
    for i = 1:num_vectors
        temp_value = sparse_cat_dot(p.pomdp, p.alphas[i], b)
        ai = actionindex(p.pomdp, p.action_map[i])
        if ( temp_value > max_values[ai])
            max_values[ai] = temp_value
        end

    end
    return max_values
 end

 # perform dot product between an alpha vector and a sparse cat object
 function sparse_cat_dot(problem::POMDP, alpha::Vector{Float64}, b::SparseCat)
    val = 0.
    for (s, p) in weighted_iterator(b)
        si = stateindex(problem, s)
        val += alpha[si]*p
    end
    return val
 end

 function AutomotivePOMDPs.action(policy::AlphaVectorPolicy, b::SingleOCFBelief)
    alphas = policy.alphas
    util = zeros(n_actions(policy.pomdp))
    for i=1:n_actions(policy.pomdp)
        res = 0.0
        for (j,s) in enumerate(b.vals)
            si = stateindex(policy.pomdp, s)
            res += alphas[i][si]*b.probs[j]
        end
        util[i] = res
    end
    ihi = findmax(util)[2]
    return policy.action_map[ihi]
end