text:
##include("Crystal.jl")
###include("DFToutMod.jl")
#using XMLDict

####################### QE specific
"""
    module AtomicProj

Analyze projections from QE projwfc.x
"""
module AtomicProj
"""
Scripts to analyze atomic projections
"""

#export load_xml
using LinearAlgebra
using DelimitedFiles

#include("Atomdata.jl")
using ..Atomdata:atoms
#using JLD
#using ExXML
using XMLDict
using GZip

using FFTW
using ..DFToutMod
#using ..CrystalMod:crystal
#using ..Utility:arr2str
#using ..Utility:str_w_spaces
using ..Utility:parse_str_ARR_float
using ..Utility:parse_str_ARR_complex
using ..Utility:write_to_file
using ..DFToutMod:bandstructure
using ..TB:make_tb
using ..TB:make_tb_crys
using ..BandTools:band_energy
#using ..TB:orbital_index
using ..CrystalMod:orbital_index

using ..TB:write_tb_crys
#using ..TB:orbital_index
using ..TB:trim
using ..QE:loadXML_bs
using ..QE:loadXML
using ..TB:calc_energy
using ..TB:types_energy
using ..CrystalMod:crystal
using ..DFT:runSCF
using ..CrystalMod:get_grid
using ..TB:myfft
using ..TB:tb_indexes
using ..TB:symm_by_orbitals

using ..TB:tb_crys_kspace
using ..TB:make_tb_crys_kspace
using ..TB:make_tb_k
using ..TB:write_tb_crys_kspace

include("Commands.jl")
using ..TightlyBound:TEMPLATEDIR

"""
    mutable struct proj_dat

Holds data from projwfc.x

- `bs::bandstructure`
- `natwfc::Int64` number of atomic wavefunctions, from dft projection
- `proj::Array{Complex{Float64}, 3}`  atomc projections `nk × natwfc × nbnd` , where `nbnd` is the number of bands in DFT
- `overlaps::Array{Complex{Float64}, 3} overlap matrix from dft `nk × natwfc × natwfc`

This is created from `loadXML_proj`, which calls `make_proj`
"""
mutable struct proj_dat

    bs::bandstructure
    natwfc::Int64
    proj::Array{Complex{Float64}, 3}
    overlaps::Array{Complex{Float64}, 3}
end

#print proj dat
Base.show(io::IO, d::proj_dat) = begin
    println(io,"projection data: nbnd = ", d.bs.nbnd, "; nkpts = ", d.bs.nks, "; natwfc = ", d.natwfc)
end

"""
    function makedict_proj(savedir)

Return xml proj data from QE savedir
"""
function makedict_proj(savedir)

    if isfile(savedir*"/atomic_proj.xml")
        filename=savedir*"/atomic_proj.xml"
    elseif isfile(savedir*"/atomic_proj.xml.gz")
        filename=savedir*"/atomic_proj.xml.gz"
    else
        println("error warning no atomic_proj.xml.gz or atomic_proj.xml")
        filename = missing
    end

    f = gzopen(filename, "r")
    fs = read(f, String)
    close(f)

    d = xml_dict(fs)

    return d
end

"""
    function make_projwfcx(prefix, tmpdir)

Create the inputfile for running projwfc.x from QE
"""
function make_projwfcx(prefix, tmpdir)
"""
makes the proj input file
"""
    c_dict = make_commands(1)

    template_file=open("$TEMPLATEDIR/template.proj")
    temp = read(template_file, String)
    close(template_file)

    temp = replace(temp, "TMPDIR" => tmpdir)
    temp = replace(temp, "PREFIX" => prefix)

    return temp

end

"""
    function run_projwfcx(projfile="proj.in"; directory="./", nprocs=1)

Run projwfc.x from QE code
"""
function run_projwfcx(projfile="proj.in"; directory="./", nprocs=1)
"""
run projwfc.x QE command
"""

    c_dict = make_commands(nprocs)
    proj = c_dict["proj"]
    command = `$proj $directory/$projfile `
    println("projwfc.x command")
    println(command)
    println()

    try
        println("Running projwfc.x")
        s = read(command, String)
        f = open(directory*"/"*projfile*".out", "w")
        write(f, s)
        close(f)

        println("Ran projwfc.x")
        println()
    catch
        println("Failed to run projwfc.x ")
        #println(s)
        return -1
    end

end

"""
    function makeOG(prefix, tmpdir )

open_grid.x inputfile. My workflow doesn't currently use this, as open_grid.x was unreliable, and
I've rewritten the code to avoid it and work directly with independent k-points only.
"""
function makeOG(prefix, tmpdir )
"""
makes the OG file
"""
    template_file=open("$TEMPLATEDIR/template_og.in")
    temp = read(template_file, String)
    close(template_file)

    temp = replace(temp, "TMPDIR" => tmpdir)
    temp = replace(temp, "PREFIX" => prefix)

    return temp

end

"""
    function run_og(filename="og.in";  directory="./", nprocs=1)

runs open_grid.x
"""
function run_og(filename="og.in";  directory="./", nprocs=1)
"""
run open_grid.x command
"""
    c_dict = make_commands(nprocs)
    og = c_dict["og"]
    command = `$og $directory/$filename`
    println("opengrid command")
    println(command)
    try
        println("Running open_grid.x")

        s = read(command, String)
        f = open(directory*"/"*filename*".out", "w")
        write(f, s)
        close(f)

        println("Ran open_grid.x")
        return 0
    catch
        println("Failed to run open_grid.x")
        return -1
    end

end

"""
    function projwfx_workf(dft::dftout)

This is the main workflow for the creation of TB matrix elements from QE DFT calculations.

Starting from a converged QE scf calculation...

1) Run NSCF calculation with extra empty bands. If you want a real space TB object, these need to be a full k-point grid not using symmtery. k-space only can use irreducible k-points

2) Run projwfc.x

3) Construct the TB hamiltonian in k-space

4) (Optional) FT to real space

5) (Optional) cleanup wavefunctions.

# Arguments

- `dft::dftout` The starting scf calculation
- `directory="./"`
- `nprocs=1` number of processors
- `freeze=true` Keep occupied eigenvalues fixed to exact DFT values
- `writefile="projham.xml"` output file for real-space TB
- `writefilek="projham_K.xml"` output file for k-space TB
- `skip_og=true`  Not used anymore
- `skip_proj=true` If projections are already run, don't run them again, load from file.
- `shift_energy=true` Shift energy of eigenvalues s.t. total energy equal band energy.
- `cleanup=true` Remove the large wavefunction files from disk, keeping nscf/projection files.
- `skip_nscf=true` If nscf calculation is already done, load results from file.
- `localized_factor = 0.15` Adjust extent of overlap matrix. 0.0 uses full atomic wavefunctions, which can be overly delocalized. 1.0 is fully localized.
- `only_kspace=false` Do not create real-space tb. Usually true in current code, as I can fit directly from k-space tb only.
- `screening = 1.0` If use a screening factor to reduce value of U in Ewald calculation. Usually leave at 1.0.
"""
function projwfx_workf(dft::dftout; directory="./", nprocs=1, freeze=true, writefile="projham.xml",writefilek="projham_K.xml", skip_og=true, skip_proj=true, shift_energy=true, cleanup=true, skip_nscf=true, localized_factor = 0.15, only_kspace=false, screening = 1.0)
"""

Steps:

2) run projwfc.x
2a) load output
3) make projected ham in k-space
4) do fourier transform
5) optionally write results to file
"""
    prefix=dft.prefix
    outdir=dft.outdir
    println()
    println("projwfc_workflow---------------------------------------------------------------------------")
    println()
    println("Step #0 nscf---------------------------------------------------------------------------------")
    s=""
    prefix_orig = deepcopy(prefix)

#    try

    olddir = "$directory/$prefix.save"
    nscfdir = "$directory/$prefix.nscf.save"

    function run_nscf()

        println("doing nscf")
        if isdir(nscfdir)
            println("removing original nscf directory $nscfdir")
            rm(nscfdir,recursive=true)
        end
        mkpath(nscfdir)
#        println("after mkdpath")

        files = readdir(olddir)
        for f in files
            if length(f) > 3 && f[end-2] == 'U' && f[end-1] == 'P' && f[end] == 'F'
#                println("found $f , copy")
                cp("$olddir/$f", "$nscfdir/$f")
            end
        end

        try
            if isfile("$olddir/charge-density.dat")
                cp("$olddir/charge-density.dat",  "$nscfdir/charge-density.dat")
            end
            if isfile("$olddir/data-file-schema.xml")
                cp("$olddir/data-file-schema.xml", "$nscfdir/data-file-schema.xml" )
            end
            if isfile("$olddir/data-file-schema.xml.gz")
                cp("$olddir/data-file-schema.xml.gz", "$nscfdir/data-file-schema.xml.gz" )

                tounzip = "$nscfdir/data-file-schema.xml.gz"
                command = `gunzip $tounzip`
                s = read(command, String)
                println("gunzipped $tounzip")
            end
            if isfile("$olddir/charge-density.dat.gz")
                cp("$olddir/charge-density.dat.gz", "$nscfdir/charge-density.dat.gz" )

                tounzip = "$nscfdir/charge-density.dat.gz"
                command = `gunzip $tounzip`
                s = read(command, String)
                println("gunzipped $tounzip")

            end

        catch
            println("missing charge density or xml file, cannot run nscf")
        end

        crys = dft.crys
        tot_charge = dft.tot_charge
        grid=dft.bandstruct.kgrid

        if only_kspace
            calc = "nscf-sym"
        else
            calc = "nscf"
        end

        try
            dft_nscf = runSCF(crys, prefix="$prefix.nscf", directory=directory,tmpdir=directory, wannier=2, nprocs=nprocs, skip=false, calculation=calc, tot_charge=tot_charge, grid=grid)
        catch
            println()
            println("first nscf failed, trying backup nscf with fewer extra bands")
            println()
            try
                dft_nscf = runSCF(crys, prefix="$prefix.nscf", directory=directory,tmpdir=directory, wannier=1, nprocs=nprocs, skip=false, calculation=calc, tot_charge=tot_charge, use_backup=true, grid=grid)
            catch
                println("try 2")
                dft_nscf = runSCF(crys, prefix="$prefix.nscf", directory=directory,tmpdir=directory, wannier=-1, nprocs=nprocs, skip=false, calculation=calc, tot_charge=tot_charge, use_backup=true, grid=grid)
            end
        end

        dft_nscf.energy = dft.energy
        dft_nscf.energy_smear = dft.energy_smear
        dft_nscf.forces = dft.forces
        dft_nscf.stress = dft.stress

        prefix = prefix*".nscf"

        return dft_nscf

    end

    if !(skip_nscf) || !(isdir(nscfdir)) ||  ( !isfile(nscfdir*"/atomic_proj.xml") && !isfile(nscfdir*"/atomic_proj.xml.gz"))

        dft_nscf = run_nscf()

    else
        if (isdir(nscfdir))
            crys = dft.crys
            #            dft_nscf = runSCF(crys, prefix="$prefix.nscf", directory=directory,tmpdir=directory, wannier=true, nprocs=nprocs, skip=false, calculation="nscf")
            dft_nscf = loadXML("$directory/$prefix.nscf.save")
            dft_nscf.energy = dft.energy
            dft_nscf.forces = dft.forces
            dft_nscf.stress = dft.stress
            prefix = prefix*".nscf"
            println("loaded nscf aaaaaaaaaaaaaaaaaaaaaa")
        else
            crys = dft.crys
            dft_nscf = dft
        end

    end

#    catch
#        println("failed nscf")
#        println(s)
#    end

#    return

    newprefix = prefix

    useog = false
    if useog
        println("Step #1 og---------------------------------------------------------------------------------")
        println()

        newprefix=prefix*"_open"

        function og()
            og= makeOG(prefix, outdir)
            write_to_file(og, "og.in", directory)
            run_og("og.in", directory=directory, nprocs=nprocs)
        end

        if skip_og
            try
                dftopen_bs = loadXML_bs("$directory/$newprefix.save") #get unfolded k-points and bs
                println("using precomputed opengrid dir")
            catch
                println("precomputinga opengrid dir")
                og()
                dftopen_bs = loadXML_bs("$directory/$newprefix.save") #get unfolded k-points and bs

            end
        else
            og()
            dftopen_bs = loadXML_bs("$directory/$newprefix.save") #get unfolded k-points and bs

        end

#        kgrid = copy(dft.bandstruct.kgrid)
#        dft_nscf.bandstruct = dftopen_bs
#        dft_nscf.bandstruct.kgrid = kgrid

    else
        println("DON'T USE OG, since it works poorly")
        dftopen_bs = loadXML_bs("$directory/$newprefix.save")
    end

    println()
    println("Step #2 run projwfc.x----------------------------------------------------------------------")
    println()

    println("A dft")
    println(dft.crys.A)
    println("A dft nscf")
    println(dft_nscf.crys.A)

    A = dft.crys.A             #needed to convert kpoints to crystal units
    a1=sum(A[1,:].^2)^0.5
    B = inv(A./a1)'

    function proj()
        projstr = make_projwfcx(newprefix, directory)
        write_to_file(projstr, "proj.in", directory)
        run_projwfcx("proj.in", directory=directory, nprocs=nprocs)
    end

    if skip_proj
        try
            p = loadXML_proj("$directory/$newprefix.save", B)
            println("using precomputed proj file")
        catch
            println("computing proj file")
            proj()
            p = loadXML_proj("$directory/$newprefix.save", B)
        end
    else
        proj()
        p = loadXML_proj("$directory/$newprefix.save", B)

    end

    if sum(abs.(dft_nscf.crys.A - dft.crys.A)) > 1e-4
        println("Warning, dft and dft_nscf may not match")
        println(dft.crys)
        println(dft_nscf.crys)
    end

    p.bs.kpts[:,:] = dft_nscf.bandstruct.kpts #fix kpoint rounding issue

    println("P OVERLAPS ", size(p.overlaps))

    #check if p matches dft_nscf
    if sum(abs.(dft_nscf.bandstruct.eigs[1,:] - p.bs.eigs[1,:])) > 1e-5
        println("warning, projected eigs and dft_nscf eigs do not match ", sum(abs.(dft_nscf.bandstruct.eigs[1,:] - p.bs.eigs[1,:])))
        println("rerun nscf")
        prefix = deepcopy(prefix_orig)
        dft_nscf = run_nscf()
        println("rerun proj")
        proj()
        p = loadXML_proj("$directory/$newprefix.save", B)

    end

    println()
    println("Step #3 ham_k------------------------------------------------------------------------------")
    println()
    projection_warning = false
    if freeze
        band_froz = Int64(round(dft_nscf.bandstruct.nelec/2))+1
        en_froz = minimum(dft_nscf.bandstruct.eigs[:,band_froz])
        en_froz = max(en_froz, dft.bandstruct.efermi + 0.05)
        println("en_froz: ", en_froz, " band froz $band_froz")
        println("efermi energy ", dft.bandstruct.efermi)
        println("nk ", size(dft_nscf.bandstruct.kpts), "     ", size(dft_nscf.bandstruct.kweights))
        ham_k, EIG, Pmat, Nmat, VAL, projection_warning = AtomicProj.create_tb(p, dft_nscf, energy_froz=en_froz+.05, shift_energy=shift_energy);
        #A,B,C = AtomicProj.create_tb(p, dft, energy_froz=en_froz+.01);
        #return A,B,C
    else
        ham_k, EIG, Pmat, Nmat, VAL, projection_warning = AtomicProj.create_tb(p, dft_nscf);
    end
    println()
    println("Step #4 ham_r------------------------------------------------------------------------------")
    println()

    println("after create P OVERLAPS ", size(p.overlaps))

    #get tight binding
    tbck = prepare_ham_k(p, dft_nscf, dft_nscf.bandstruct.kgrid ,ham_k, nonorth=true, localized_factor = localized_factor, screening=screening)

    println("after prepare P OVERLAPS ", size(p.overlaps))

    if !ismissing(writefilek)
        println("Step #5-0 write_tb_k-------")
        write_tb_crys_kspace("$directory/$writefilek", tbck);
    end

    if !only_kspace
        tbc = AtomicProj.get_ham_r(tbck)
        if !ismissing(writefile)
            println()
            println("Step #5 write_tb---------------------------------------------------------------------------")
            println()
            write_tb_crys("$directory/$writefile", tbc);
        end
    else
        tbc = missing
    end

#    tbc = AtomicProj.get_ham_r(p, dft_nscf, dft_nscf.bandstruct.kgrid ,ham_k, nonorth=true, localized_factor = localized_factor);

    #tbc = AtomicProj.get_ham_r( p, dft_nscf, dft_nscf.bandstruct.kgrid ,ham_k, nonorth=false);

    #tbc = AtomicProj.get_ham_r_slow(p, dft_nscf, dft.bandstruct.kgrid ,ham_k, nonorth=true);

#    trim(tbc.tb, 0.00005)

    println()
    println("Done proj----------------------------------------------------------------------------------")
    println()

    if cleanup
        println("clean: warning, removing wavefunctions to save space")

        newdir = "$directory/$newprefix.save"
        olddir = "$directory/$prefix.save"
        parentdir = "$directory"

        if (prefix_orig != prefix) && (newprefix != prefix_orig)
            origdir = "$directory/$prefix_orig.save"
            toclean = [newdir,olddir,origdir]
        else
            toclean = [newdir,olddir]
        end

        if (isdir(nscfdir))
            toclean = [toclean; nscfdir]
        end
        if (isdir(parentdir))
            toclean = [toclean; parentdir]
        end

        tot = 0
        for d in toclean
            files=read(`ls $d/`, String)
            for f in split(files, "\n")
                if occursin("wfc", f)
                    if tot <= 2
                        println("removing $d/$f")
                    elseif tot==3
                        println("removing $d/$f")
                        println("...")
                    end
                    cmd=`\rm $d/$f`
                    s = read(cmd, String)
                    tot += 1
                end
            end
        end
        println("clean-up done: $tot wfc files deleted")
    else
        println("don't clean wfcs")
    end

    return tbc, tbck, projection_warning

end

"""
    function loadXML_proj(savedir, B=missing)

Load proj from QE output file. Need the QE save dir like "qe.save". `B` are reciprocal lattice vectors.
"""
function loadXML_proj(savedir, B=missing)

    d= makedict_proj(savedir)

    da = d["ATOMIC_PROJECTIONS"] #has all the real data in the xml

    nbnd = parse(Int, da["HEADER"]["NUMBER_OF_BANDS"][""])
    nk = parse(Int, da["HEADER"]["NUMBER_OF_K-POINTS"][""])
    nspin = parse(Int, da["HEADER"]["NUMBER_OF_SPIN_COMPONENTS"][""])
    natwfc = parse(Int, da["HEADER"]["NUMBER_OF_ATOMIC_WFC"][""])
    nelec = parse(Float64, da["HEADER"]["NUMBER_OF_ELECTRONS"][""])
    efermi = parse(Float64, da["HEADER"]["FERMI_ENERGY"][""])

    units_energy = da["HEADER"]["UNITS_FOR_ENERGY"][:UNITS]
    units_kpt = da["HEADER"]["UNITS_FOR_K-POINTS"][:UNITS]

    println(units_energy, " " , units_kpt)
    println("nbands: ", nbnd," nk: ", nk)

    kpts = parse_str_ARR_float(da["K-POINTS"][""])
    kpts = reshape(kpts, 3,nk)'

    if !(ismissing(B))
        kpts = kpts * inv(B)
    end

    weights = parse_str_ARR_float(da["WEIGHT_OF_K-POINTS"][""])

    dae = da["EIGENVALUES"]

    eigs = zeros(nk,nbnd)
    for k in 1:nk
        eigs[k, :] = parse_str_ARR_float(dae["K-POINT."*string(k)]["EIG"][""])

    end

    dap = da["PROJECTIONS"]

    proj = zeros(Complex{Float64}, nk, natwfc, nbnd)
    for k in 1:nk
        t = dap["K-POINT."*string(k)]
        for a in 1:natwfc
            proj[k, a, :] =  parse_str_ARR_complex(t["ATMWFC."*string(a)][""])
        end
    end

    dao = da["OVERLAPS"]
    overlaps = zeros(Complex{Float64}, nk, natwfc, natwfc)

    for k in 1:nk
        t = dao["K-POINT."*string(k)]
        overlaps[k, :,:] = reshape(parse_str_ARR_complex(t["OVERLAP.1"][""]), natwfc, natwfc)'
    end

#    println(overlaps)

    bs = DFToutMod.makebs(nelec, efermi, kpts, weights, [0,0,0], eigs)

    return make_proj(bs, proj, overlaps)

end

"""
    function make_proj(bs, proj, overlaps)

Constructor for proj_dat.
"""
function make_proj(bs, proj, overlaps)

    if bs.nks != size(proj)[1] || bs.nks != size(overlaps)[1]
        error("make_proj something wrong nks ", bs.nks," ",size(proj)[1]," ",size(overlaps)[1])
    end
    if bs.nbnd != size(proj)[3]
        error("make_proj something wrong nbnd ", bs.nbnd," ",size(proj)[3])
    end
    if size(proj)[2] != size(overlaps)[2] ||  size(proj)[2] !=	size(overlaps)[3]
        error("make_proj something wrong in natwfc ",  size(proj)[2]," ",size(overlaps)[2]," " ,size(overlaps)[3])
    end

    return proj_dat(bs, size(proj)[2], proj, overlaps)

end

#function create_temp(p::proj_dat, d::dftout; energy_froz=missing, nfroz=0, shift_energy=true)
##
#
#    return
#end

"""
    function create_tb(p::proj_dat, d::dftout; energy_froz=missing, nfroz=0, shift_energy=true)

Does the main creation of TB hamiltonian from DFT projection data in k-space.

# Arguments
- `p::proj_dat` Projection data
- `d::dftout` DFT scf data
- `energy_froz=missing` Energy to start freezing eigenvalues
- `nfroz=0` number of frozen bands.
- `shift_energy=true` if `true` shift eigenvalues so band energy `==` total energy
"""
function create_tb(p::proj_dat, d::dftout; energy_froz=missing, nfroz=0, shift_energy=true)

    wan, semicore, nwan, nsemi, wan_atom, atom_wan = tb_indexes(d)

    efermi_dft = d.bandstruct.efermi

    #decide semicore
    p2 = zeros(p.bs.nbnd)
    INDSEMI = zeros(Int64, p.bs.nks, nsemi)
    for k = 1:p.bs.nks
        p2[:] = real.(sum(p.proj[k,semicore,:] .* conj.(p.proj[k,semicore,:]) , dims=1))
        INDSEMI[k, :] = sortperm(p2, rev=true)[1:nsemi]
        if sum(p2) < nsemi - 0.5
            println("warning, identify semicore")
        end
        if k < 10
            println("p2 ", p2)
            println(k, " indsemi " , INDSEMI[k, :])
        end

    end
    println("INDSEMI ", INDSEMI)

    NBND = p.bs.nbnd - nsemi
    EIGS = zeros(p.bs.nks, NBND)
    PROJ = zeros(Complex{Float64}, p.bs.nks, nwan, NBND)

    # setup data
    for k = 1:p.bs.nks
        counter = 0
        for n = 1:p.bs.nbnd
            if !(n in INDSEMI[k, :])
                counter += 1
                EIGS[k,counter] = p.bs.eigs[k,n]
                PROJ[k,:,counter] = p.proj[k,wan,n]

            end
        end
#        if k < 20
#            println("$k EIGS ", EIGS[k,1:5])
#        end

#        if k < 2
#            println("ALLEIGS k $k ", p.bs.eigs[k,1:4])
#            println("EIGS    k $k ", EIGS[k,1:4])#
#
#        end
    end

    if shift_energy

        println("shifting eigenvalues to match dft atomization energy")
        ind2orb, orb2ind, etotal_atoms, nval =  orbital_index(d.crys)
#        band_en = band_energy(d.bandstruct.eigs[:,nsemi+1:end], d.bandstruct.kweights, nval)
        band_en = band_energy(EIGS, d.bandstruct.kweights, nval)

#        println("d.crys")
#        println(d.crys)

        etypes = types_energy(d.crys)

        etot_dft = d.energy
        e_smear = d.energy_smear

        atomization_energy = etot_dft - etotal_atoms - etypes  - e_smear

        println("atomization_energy $atomization_energy")

        band_en = band_en
        shift = (atomization_energy - band_en  )/nval

        EIGS = EIGS .+ shift

        if !(ismissing(energy_froz))
            println("energy_froz before $energy_froz")
#            energy_froz = energy_froz + shift
            energy_froz = efermi_dft + shift

        end
        println("energy_froz after $energy_froz")
        println("sum shifted EIGS ", sum(EIGS))
    else
        println("no shift: match dft eigenvals")
    end

    P = zeros(Complex{Float64}, NBND, NBND)  #p.bs.nbnd-nsemi,p.bs.nbnd-nsemi)

    Pmat = zeros(Float64, p.bs.nks, NBND)
    Nmat = zeros(Float64, p.bs.nks, nwan)

    Btilde = zeros(Complex{Float64}, nwan, NBND)
    B = zeros(Complex{Float64}, nwan, NBND)

    ham_dft = zeros(Float64, NBND)

    ham_k = zeros(Complex{Float64},  p.natwfc-nsemi, p.natwfc-nsemi,p.bs.nks)

    htemp = zeros(Complex{Float64}, nwan, nwan)

    max_c = minimum(EIGS[:,end])
    min_c = max_c - 0.1

    min_c = max(maximum(EIGS[:,1:nwan])+.001, min_c) #ensure we don't cut off needed states

    #simple cutoff function
    function cutoff(num, min_c, max_c)
        if num < min_c
            return 1.0
        elseif num > max_c
            return 0.0
        else
            t = (num - min_c)/(max_c-min_c)
            return 1.0 - 10.0 * t^3 + 15.0 *  t^4  - 6.0 * t^5
        end
    end

    Palt_min = 1.0

    badk=0
    for k = 1:p.bs.nks
#    for k = 1:nks

        #this avoids breaking symmetry by cutting off a symmetrically equivalent pair/triplet, etc because we include a finite number of bands

        energy_cutoff = EIGS[k,end]
        max_ind = 0

        for i=max(1,NBND-5):NBND
#            if p.bs.eigs[k,i] < energy_cutoff - 1e-3
            if EIGS[k,i] < energy_cutoff - 1e-3
                max_ind = i
            end
        end

#        B[:,1:max_ind-nsemi] = p.proj[k,wan, 1+nsemi:max_ind]
        B[:,1:max_ind] = PROJ[k,:, 1:max_ind]

        #smoothing
        for i in 1:max_ind
#            B[:,i] *= cutoff(p.bs.eigs[k,i+nsemi], min_c, max_c)
            B[:,i] *= cutoff(EIGS[k,i], min_c, max_c)
        end

        P[:,:] .= 0.0

        P[1:max_ind,1:max_ind] = B[:,1:max_ind]'*B[:,1:max_ind]

        Palt = B[:,1:max_ind]*B[:,1:max_ind]'

        for i in 1:nwan

            if real(Palt[i,i]) < Palt_min
                Palt_min = real(Palt[i,i])
            end

            #warn user if projection is bad, i.e. projectors aren't projecting onto anything.
            if real(Palt[i,i]) < 0.80
                println("Warning, difficult to project atomic wfc ", i," ",  real(Palt[i,i]), " k ", k, " tr(P) = " , tr(real(P)))
            end
            if real(Palt[i,i]) < 0.45
                badk += 1
                println("bad k ", badk)
            end
            if real(Palt[i,i]) < 0.35
                badk += 1
                println("very bad k ", badk)
            end

        end

        #eigenvalues of projection matrix are key to this method
        val, vect = eigen( 0.5 * (P[1:max_ind,1:max_ind] + P[1:max_ind,1:max_ind]')     ) #

        good_proj = (max_ind - nwan + 1 ) : max_ind

        Btilde[:,1:max_ind] = vect[:,good_proj]'

        #approximate hamiltonian using highest eigenvalues of projection matrix.

        htemp[:,:] = Btilde[:,1:max_ind] * Diagonal(EIGS[k,1:max_ind])  * Btilde[:,1:max_ind]'

        neweigs, newvect = eigen((htemp+htemp')/2.0)

        Pmat[k, :] = real(diag(P))
        Nmat[k, :] = neweigs

        ham_k[:,:,k] = B[:,1:max_ind] * Btilde[:,1:max_ind]' * htemp * Btilde[:,1:max_ind] *  B[:,1:max_ind]'

        ham_k[:,:,k] = (ham_k[:,:,k]  + ham_k[:,:,k]')/2.0

    end

    #decide if need to send serious warning to user.
    if badk / p.bs.nks > 0.15
        println("warning lots of bad kpoints ", badk, " of ", p.bs.nks)
        warn_badk=true
    else
        warn_badk=false
    end

    println("Min Palt: $Palt_min")

    #here we shift the eigenvalues around to match DFT eigenvalues below a cutoff.
    #this requires identifying which bands are supposed to match which eigenvlues, which
    #can be tricky. I use a projection heuriestic, but it can fail at band crossings that mix bands.
    if !(ismissing(energy_froz))
        energy_froz2 = energy_froz+0.25
        println("energy_froz: $energy_froz , $energy_froz2")
        for k = 1:p.bs.nks
            val_tbt, vect = eigen(ham_k[:,:,k] )
            val_tb = real(val_tbt)
            val_tb_new = deepcopy(val_tb)
#            val_pw = p.bs.eigs[k,nsemi+1:nsemi+nwan]
            val_pw = EIGS[k,1:nwan]

            order = Dict()

            score_mat = zeros(nwan, nwan)

            for n2 = 1:nwan
                cd_vect = real(vect[:,n2].*conj(vect[:,n2]))

                dist_min = 1000.0
                nmin = 0

                for n1 = 1:nwan

#                    t = p.proj[k,wan, nsemi + n1]
                    t = PROJ[k,:,n1]

                    cd_dft = real(t .* conj(t))

                    dist_en = (val_tb[n2] - val_pw[n1]).^2 * 10.0
                    dist_cd = sum((cd_dft - cd_vect).^2)
                    dist = dist_en + dist_cd + 0.1*(n1 - n2)^2
                    score_mat[n1,n2] = dist
#                    if dist < dist_min
#                        dist_min = dist
#                        nmin = n1
#                    end
                end
            end
            for n = 1:nwan
                s = sortperm(score_mat[n,:])
                for ss in s
                    if !(ss in keys(order))
                        order[ss] = n
                        break
                    end
                end
            end

#acutally do the change
            for n in 1:nwan
                if val_pw[order[n]] < energy_froz
                    val_tb_new[n] = val_pw[order[n]]
                elseif val_pw[order[n]] < energy_froz2
                    x=cutoff(val_tb[n], energy_froz, energy_froz2)
                    val_tb_new[n] = val_pw[order[n]] * x + val_tb[n]*(1.0-x)
                end

            end

#resym
            for n = 1:nwan-1
                c= [n]
                for n2 = n+1:nwan
                    if abs(val_tb_new[n] - val_tb_new[n2]) < 1e-4
                        push!(c, n2)
                    end
                end
                if length(c) > 1
                    t = sum(val_tb_new[c]) / length(c)
                    val_tb_new[c] .= t
                end
            end

            ham_k[:,:,k] = vect*Diagonal(val_tb_new)*vect'
            ham_k[:,:,k] = (ham_k[:,:,k]  + ham_k[:,:,k]')/2.0
            val_tb_new, vect = eigen(ham_k[:,:,k] )

        end
    end

#alternate eigenvalue fix method.
    if nfroz >= 1 && (ismissing(energy_froz))
        println("nfroz: ", nfroz)
        for k = 1:p.bs.nks
            val, vect = eigen(ham_k[:,:,k] )
            val[1:nfroz] = EIGS[k,1:nfroz]
            if (abs(EIGS[k, nfroz+1] - EIGS[k, nfroz])< 1e-5) && nwan >= nfroz+1
                val[nfroz+1] = EIGS[k,nfroz+1]
                if (abs(EIGS[k, nfroz+2] - EIGS[k, nfroz])< 1e-5) && nwan >= nfroz+2
                    val[nfroz+2] = EIGS[k,nfroz+2]
                end
            end

            ham_k[:,:,k] = vect*Diagonal(val)*vect'
            ham_k[:,:,k] = (ham_k[:,:,k]  + ham_k[:,:,k]')/2.0

        end
    end

    #reshift to match dft total energy

    VAL = zeros(Float64, p.bs.nks, nwan)
    VECT = zeros(Complex{Float64}, p.bs.nks, nwan, nwan)
    for k = 1:p.bs.nks
        val, vect = eigen(ham_k[:,:,k])
        VAL[k,:] = val
        VECT[k,:,:] = vect
    end

    if shift_energy
        println("reshift")
#        println("size VAL ", size(VAL), " size kweights ", size(d.bandstruct.kweights))
        band_en = band_energy(VAL, d.bandstruct.kweights, nval)

        println("band_en_old " , band_en)

        shift = (atomization_energy - band_en)/nval
        VAL = VAL .+ shift
        for k = 1:p.bs.nks

            val = VAL[k,:]
            vect = VECT[k,:,:]
            ham_k[:,:,k] = vect*Diagonal(val)*vect'

        end
        println("done reshift")

        for k = 1:p.bs.nks
            val, vect = eigen((ham_k[:,:,k] + ham_k[:,:,k]')/2.0)
            VAL[k,:] = val
            VECT[k,:,:] = vect
        end

        band_en = band_energy(VAL, d.bandstruct.kweights, nval)

        println("band_en_new " , band_en, " " , band_en+etypes)

    end

    return ham_k, EIGS, Pmat, Nmat, VAL, warn_badk

end

"""
    function get_ham_r_slow(p::proj_dat, d::dftout, grid, ham_k::Array{Complex{Float64}, 3}; nonorth=true, K=missing)

Slow method for doing fourier transform. Uses standard ft formula, not fft.

Use fast version instead.
"""
function get_ham_r_slow(p::proj_dat, d::dftout, grid, ham_k::Array{Complex{Float64}, 3}; nonorth=true, K=missing)
"""
DOESN'T USE FFTW. MOSTLY FOR DEBUGGING / missing fft libraries

Do Fourier transform k->R. grid is the size of the k-space grid.
Optionally you can include the k space grid explictly, but normally
K is already inside p. K must match ham_k
"""

    wan, semicore, nwan, nsemi, wan_atom, atom_wan = tb_indexes(d)

#    grid = d.bandstruct.kgrid

    println("grid ", grid)
    println("ham_k size " , size(ham_k))
    println("nonorth: ", nonorth)
    nwan = size(ham_k)[2]

    if ismissing(K)
        K = p.bs.kpts
        nks = size(K)[1]
    else
        nks = size(K)[1]
    end

    if nonorth
        #nonorthogonal!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ham_kS = zeros(Complex{Float64},  nwan, nwan,nks)
        S = zeros(Complex{Float64}, nwan, nwan)
        Sk = zeros(Complex{Float64}, nwan, nwan,nks)
        Sk2 = zeros(Complex{Float64},  nwan, nwan,nks)
        for k in 1:nks
            S[:,:] = p.overlaps[k, wan,wan]
            S = (S+S')/2.0
            Sk[:,:,k] =  S
        end
        for i=1:nwan #normalize overlaps
            sumi = sum(Sk[i,i,:])
            Sk[i,i,:] = Sk[i,i,:] / real(sumi) * Float64(nks)
        end
        for k in 1:nks
            S[:,:] = Sk[:,:,k]
            S = (S + S')/2.0
            Sk2[:,:,k] = sqrt(S)

            ham_kS[:,:,k] = Sk2[:,:,k]*ham_k[:,:,k]*Sk2[:,:,k]
            ham_kS[:,:,k] = (ham_kS[:,:,k]  + ham_kS[:,:,k]')/2.0

        end
        #end nonorthogonal!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    end

    r=zeros(Float64,3)
    twopi_i = 1.0im*2.0*pi

    grid2 = [0,0,0]
    for i = 1:3
        if grid[i]%2 == 0
            grid2[i] = Int(grid[i]/2)
        else
            grid2[i] = Int((grid[i]-1)/2)
        end
    end

    r_dict = Dict()

    ngrid2 = (grid2[1]*2+1)*(grid2[2]*2+1)*(grid2[3]*2+1)
    println("new real space grid: " , -grid2, " to " , grid2, " " , ngrid2)

    ind_arr = zeros(ngrid2,3)
    ham_r = zeros(Complex{Float64},  nwan, nwan,ngrid2)

    if nonorth
        S_r = zeros(Complex{Float64},  nwan, nwan,ngrid2)
    end

#    return 0

    sym = zero(Float64)
    exp_val = zero(Complex{Float64})
    c=0

    sym_mat = zeros(Float64, ngrid2)

    for r1 = -grid2[1]:(grid2[1])

        if (r1 + grid[1] == grid2[1]) ||  (r1 - grid[1] == -grid2[1])
            sym1 = 2.0
        else
            sym1 = 1.0
        end

        for r2 = -grid2[2]:(grid2[2])
            if (r2 + grid[2] == grid2[2]) ||  (r2 - grid[2] == -grid2[2])
                sym2 = 2.0
            else
                sym2 = 1.0
            end

            for r3 = -grid2[3]:(grid2[3])

                if (r3 + grid[3] == grid2[3]) ||  (r3 - grid[3] == -grid2[3])
                    sym3 = 2.0
                else
                    sym3 = 1.0
                end

#                println("rs $r1 $r2 $r3 $sym1 $sym2 $sym3 $sym")

                c+=1
                sym_mat[c] = sym1 * sym2 * sym3
                r[:] = [r1 r2 r3]
                ind_arr[c,:] = r[:]
                r_dict[ind_arr[c,:]] = c
            end
        end
    end

    if nonorth
        for k in 1:nks
            for c = 1:ngrid2
                exp_val = exp(twopi_i * (ind_arr[c,:]'*K[k,:])) / sym_mat[c]
                ham_r[:,:, c ] .+= ham_kS[:,:,k] .* exp_val
                S_r[:,:, c] .+=  Sk[:,:,k] .* exp_val

            end
        end

    else
        for k in 1:nks
            for c = 1:ngrid2
                exp_val = exp(twopi_i * (ind_arr[c,:]'*K[k,:])) / sym_mat[c]
                ham_r[ :,:, c] .+= ham_k[:,:,k] .* exp_val
            end
        end
    end

    #turn ham_r into tight-binding object
    if nonorth
        tb = make_tb(ham_r/ Float64(nks), ind_arr, r_dict, S_r / Float64(nks))
    else
        tb = make_tb(ham_r/ Float64(nks), ind_arr, r_dict)
    end

#    renormalize_tb(d, tb)

    nelec=d.bandstruct.nelec - nsemi * 2.0

    ind2orb, orb2ind, etotal_atoms, nval =  orbital_index(d.crys)
    dftenergy=d.energy - etotal_atoms

    model = make_tb_crys(tb, d.crys, nelec, dftenergy)

    return model

end

"""
    function prepare_ham_k(p::proj_dat, d::dftout, grid, ham_k::Array{Complex{Float64}, 3}; nonorth=true, K=missing, localized_factor = 0.05, screening=1.0)

Constructs the actual k-space hamiltonain, which involves lots of putting matricies in the correct form and de-orthogonalizing if desired.

#arguments
- `p::proj_dat` projection data
- `d::dftout` dft data
- `grid` kpoint grid spacing
- `ham_k::Array{Complex{Float64}, 3}` Hamiltonian
- `nonorth=true` make a non-orthogonal TB
- `K=missing` k-point array, usually get from `p`
- `localized_factor = 0.15` increase localization of overlaps.
- `screening=1.0` mulitply U by this factor. usually not used.
"""
function prepare_ham_k(p::proj_dat, d::dftout, grid, ham_k::Array{Complex{Float64}, 3}; nonorth=true, K=missing, localized_factor = 0.15, screening=1.0)

    wan, semicore, nwan, nsemi, wan_atom, atom_wan = tb_indexes(d)

    println("grid ", grid)
    println("ham_k size " , size(ham_k))
    println("nonorth: ", nonorth)
    nwan = size(ham_k)[2]

    if ismissing(K)
        K = p.bs.kpts
        nks = size(K)[1]
    else
        nks = size(K)[1]
    end

    #fix mysterious sign error K-space. QE overlaps have weird signs in x and y direction
    # I don't know why, you have to ask them.

    ind2orb, orb2ind, etotal_atoms, nval =  orbital_index(d.crys)

    OVERLAPS = deepcopy(p.overlaps)
    println("wan ", wan)
    println("so ", size(OVERLAPS))
    OVERLAPS = OVERLAPS[:,wan, wan]
    println("so2 ", size(OVERLAPS))
    for n = 1:nwan
        a1,t1,orb = ind2orb[n]
        if orb == :px || orb == :py
            println("fix mysterious sign issue K-space $n, $a1, $t1, $orb")
            ham_k[n,:,:] = -1.0 * ham_k[n,:,:]
            ham_k[:,n,:] = -1.0 * ham_k[:,n,:]

            if nonorth
                OVERLAPS[:,n,:] = -1.0*OVERLAPS[:,n,:]
                OVERLAPS[:,:,n] = -1.0*OVERLAPS[:,:,n]

            end
        elseif orb == :dxz || orb == :dyz
            println("fix mysterious sign issue K-space D ORBITALS $n, $a1, $t1, $orb")
            ham_k[n,:,:] = -1.0 * ham_k[n,:,:]
            ham_k[:,n,:] = -1.0 * ham_k[:,n,:]

            if nonorth
                OVERLAPS[:,n,:] = -1.0*OVERLAPS[:,n,:]
                OVERLAPS[:,:,n] = -1.0*OVERLAPS[:,:,n]

            end

        end
    end

    if nonorth
        #nonorthogonal!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ham_kS = zeros(Complex{Float64},  nwan, nwan,nks)
        S = zeros(Complex{Float64}, nwan, nwan)
        Sk = zeros(Complex{Float64}, nwan, nwan,nks)
        Sk2 = zeros(Complex{Float64},  nwan, nwan,nks)

        for k in 1:nks
#            S[:,:] = p.overlaps[k, wan,wan]
            S[:,:] = OVERLAPS[k, :,:]

            S = (S+S')/2.0           #+ I(nwan) * localized_factor
            Sk[:,:,k] =  S
        end

#        println("fourteen")
#        println(Sk[:,:,14])

        wtot = sum(p.bs.kweights)
        SUMI = zeros(nwan)
        for i=1:nwan #normalize overlaps
            SUMI[i] = sum(Sk[i,i,:] .* p.bs.kweights) / wtot
        end
        println("old SUMI ", SUMI)
        SUMI = symm_by_orbitals(d.crys, SUMI) #this handles the fact the k-grid can break orbital symmetry
        println("new SUMI ", SUMI)
        for i=1:nwan #normalize overlaps
            Sk[i,i,:] = Sk[i,i,:] / SUMI[i]
        end

        for k in 1:nks
            S = Sk[:,:,k]
            S = (S + S')/2.0
            Sk[:,:,k] = S
        end

        if localized_factor > 1e-5
            for k = 1:nks
                Sk[:,:,k] = Sk[:,:,k] * (1.0 - localized_factor) + I(nwan) * localized_factor
            end
        end

        #does the de-orthogonalization
        for k in 1:nks
            S[:,:] = Sk[:,:,k]
            S = (S+S')/2.0
            Sk2[:,:,k] = sqrt(S)

            ham_kS[:,:,k] = Sk2[:,:,k]*ham_k[:,:,k]*Sk2[:,:,k]
            ham_kS[:,:,k] = (ham_kS[:,:,k]  + ham_kS[:,:,k]')/2.0

        end

        tbk = make_tb_k(ham_kS, K, d.bandstruct.kweights, Sk, grid=grid, nonorth=true)
        tbck = make_tb_crys_kspace(tbk, d.crys, nval, d.energy - etotal_atoms, scf=false, screening=screening)

        return tbck

        #end nonorthogonal!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    else

        Sk = zeros(Complex{Float64}, nwan, nwan,nks)
        for i = 1:nks
            Sk[i,:,:] = I(nwan)
        end

        tbk = make_tb_k(ham_k, K, d.bandstruct.kweights, Sk, grid=grid, nonorth=false)
        tbck = make_tb_crys_kspace(tbk, d.crys, nval, d.energy - etotal_atoms, scf=false)

        return tbck

    end

    ###########

end

"""
    function get_ham_r(tbck::tb_crys_kspace)

Do fft to get the real-space ham. Requires tbck to be on a regular grid centered at Gamma with no symmetry.
"""
function get_ham_r(tbck::tb_crys_kspace)
"""
Do Fourier transform k->R.
"""

    nonorth = tbck.tb.nonorth
    ham_k = tbck.tb.Hk
    Sk = tbck.tb.Sk
    K = tbck.tb.K
    grid = tbck.tb.grid
    nwan = tbck.tb.nwan

    if nonorth
        ham_r, S_r, r_dict, ind_arr = myfft(tbck.crys, nonorth, grid, K,ham_k, Sk)
    else
        ham_r, r_dict, ind_arr =      myfft(tbck.crys, nonorth, grid, K,ham_k, missing)
    end

    ind2orb, orb2ind, etotal_atoms, nval =  orbital_index(tbck.crys)

    #turn ham_r into tight-binding object
    if nonorth
        tb = make_tb(ham_r, ind_arr, r_dict, S_r )
    else
        S_r = zeros(size(ham_r))

        println(size(S_r))

        r = r_dict[[0,0,0]]
        for i =  1:nwan
            S_r[i,i,r] = 1.0  #trivial overlaps
        end
        tb = make_tb(ham_r, ind_arr, r_dict, S_r)
    end

#    renormalize_tb(d, tb)
#    nelec=.bandstruct.nelec ##- nsemi * 2.0
    nelec = nval

    dftenergy=tbck.dftenergy
    ind2orb, orb2ind, etotal_atoms, nval =  orbital_index(tbck.crys)

    model = make_tb_crys(tb, tbck.crys, nelec, dftenergy )

    return model

end

##function run_proj(inputstr, outputstr, nprocs=1, directory="./")
##"""
##run command
##"""
##    directory=rstrip(directory, '/')
##
##    #get commandline
##    c_dict = make_commands(nprocs)
##    qe = c_dict["qe"]
##
##
##    command = `$qe $directory/$inputstr`
##    println("actual command")
##    println(command)
##    s=""
##    try
##        println("Running DFT")
##
##        s = read(command, String)
##
##        println("Writing output")
##        f = open(directory*"/"*outputstr, "w")
##        write(f, s)
##        close(f)
##
##        println("Ran DFT")
##        return 0
##    catch
##        println("Failed to run DFT")
##        println(s)
##        return -1
##    end
##
##end
##

end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdpopconfirm

"""
    'feffery'_antdpopconfirm(;kwargs...)
    'feffery'_antdpopconfirm(children::Any;kwargs...)
    'feffery'_antdpopconfirm(children_maker::Function;kwargs...)

An AntdPopconfirm component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): The content of the tab - will only be displayed if this tab is selected
- `id` (String; optional)
- `cancelButtonProps` (Dict; optional)
- `cancelCounts` (Real; optional)
- `cancelText` (String; optional)
- `className` (String; optional)
- `confirmCounts` (Real; optional)
- `disabled` (Bool; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `locale` (a value equal to: 'zh-cn', 'en-us'; optional)
- `mouseEnterDelay` (Real; optional)
- `mouseLeaveDelay` (Real; optional)
- `okButtonProps` (Dict; optional)
- `okText` (String; optional)
- `overlayInnerStyle` (Dict; optional)
- `overlayStyle` (Dict; optional)
- `placement` (a value equal to: 'top', 'left', 'right', 'bottom', 'topLeft', 'topRight', 'bottomLeft', 'bottomRight', 'leftTop', 'leftBottom', 'rightTop', 'rightBottom'; optional)
- `style` (Dict; optional)
- `title` (String; optional)
- `trigger` (a value equal to: 'hover', 'focus', 'click' | Array of a value equal to: 'hover', 'focus', 'click's; optional)
"""
function 'feffery'_antdpopconfirm(; kwargs...)
        available_props = Symbol[:children, :id, :cancelButtonProps, :cancelCounts, :cancelText, :className, :confirmCounts, :disabled, :loading_state, :locale, :mouseEnterDelay, :mouseLeaveDelay, :okButtonProps, :okText, :overlayInnerStyle, :overlayStyle, :placement, :style, :title, :trigger]
        wild_props = Symbol[]
        return Component("'feffery'_antdpopconfirm", "AntdPopconfirm", "feffery_antd_components", available_props, wild_props; kwargs...)
end

'feffery'_antdpopconfirm(children::Any; kwargs...) = 'feffery'_antdpopconfirm(;kwargs..., children = children)
'feffery'_antdpopconfirm(children_maker::Function; kwargs...) = 'feffery'_antdpopconfirm(children_maker(); kwargs...)

ext: jl

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
module TreeTensors

importall Base
using Tensors
include("Trees.jl")
include("ModeTrees.jl")

# Tree tensor networks

export TreeTensor
typealias TensorDict{T} Dict{Tree,Tensor{T}}
immutable TreeTensor{T}
    mtree::AbstractModeTree
    tensors::TensorDict{T}
end
(::Type{TreeTensor{T}}){T}(mtree) = TreeTensor(mtree, TensorDict{T}())

# Basic functions

export modetree, tree
Tensors.scalartype{T}(::Type{TreeTensor{T}}) = T
modetree(x::TreeTensor) = x.mtree
tree(x::TreeTensor) = tree(modetree(x))
copy(x::TreeTensor) = TreeTensor(x.mtree, TensorDict{scalartype(x)}(x.tensors))
convert{T}(::Type{TreeTensor{T}}, x::TreeTensor) = TreeTensor(x.mtree, convert(TensorDict{T},x.tensors))
copyconvert{T}(::Type{TreeTensor{T}}, x::TreeTensor{T}) = copy(x)
copyconvert{T}(::Type{TreeTensor{T}}, x::TreeTensor) = convert(TreeTensor{T}, x)

# Indexing and iteration

for f in (
    :getindex, :setindex!,
    :start, :next, :done,
    :eltype, :length,
    :keys, :values, :keytype, :valtype
)
    @eval $f(x::TreeTensor, args...) = $f(x.tensors, args...)
end
function index(x::TreeTensor, s::Symbol)
    if s == :root return tree(x)
    else throw(KeyError(s)) end
end
getindex(x::TreeTensor, s::Symbol) = x[index(x,s)]
setindex!(x::TreeTensor, xv, s::Symbol) = x[index(x,s)] = xv

# Initialization

ones{T}(::Type{T}, mtree::AbstractModeTree) = TreeTensor(mtree, Dict{Tree,Tensor{T}}(
    v => ones(T, [[Mode(e,1) for e in neighbor_edges(v)]; mtree[v]])
    for v in vertices(mtree, root_to_leaves)
))
eye{T}(::Type{T}, mtree::AbstractModeTree) = TreeTensor(square(mtree), Dict{Tree,Tensor{T}}(
    v => begin
        t = eye(T, mtree[v]);
        t.modes = [t.modes; [Mode(e,1) for e in neighbor_edges(v)]];
        t
    end
    for v in vertices(mtree, root_to_leaves)
))
function rand{T}(::Type{T}, mtree::AbstractModeTree, r)
    evaluate(r::Int,e) = r
    evaluate(r::Dict,e) = r[e]
    return TreeTensor(mtree, Dict{Tree,Tensor{T}}(
        v => rand(T, [[Mode(e,evaluate(r,e)) for e in neighbor_edges(v)]; mtree[v]])
        for v in vertices(mtree, root_to_leaves)
    ))
end
for f in (:ones, :rand, :eye)
    @eval $f(mtree::AbstractModeTree, args...) = $f(Float64, mtree, args...)
end
for f in (:ones, :rand, :eye)
    @eval $f(x::TreeTensor, args...) = $f(scalartype(x), modetree(x), args...)
end

# Conversion to and from full

export decompose
function decompose(x::Tensor, mtree, rank)
    y = TreeTensor{scalartype(x)}(mtree)
    for (v,p) in edges(y, leaves_to_root)
        b,s,y[v] = svd(x, [child_edges(v,p);mlabel(mtree[v])], PairSet(v,p), rank)
        x = scale!(b,s)
    end
    y[:root] = x
    return y
end

export contract!, contract
function contract!(x::TreeTensor)
    for (v,p) in edges(x, leaves_to_root)
        x[p] *= x[v]
    end
    return x[:root]
end
contract(x::TreeTensor) = contract!(copy(x))

# Arithmetic

function +(x::TreeTensor, y::TreeTensor)
    @assert modetree(x) == modetree(y) "x and y must have the same mode tree"
    mtree = modetree(x)
    return TreeTensor(
        mtree,
        Dict{Tree,Tensor{promote_type(scalartype(x), scalartype(y))}}(
            v => padcat(x[v],y[v], neighbor_edges(v))
            for v in vertices(mtree, root_to_leaves)
        )
    )
end

function *(x::TreeTensor, y::TreeTensor)
    @assert unsquare(modetree(x)) == unsquare(modetree(y)) "x and y must have the same mode tree"
    mtree = modetree(x)
    return TreeTensor(
        mtree,
        Dict{Tree,Tensor{promote_type(scalartype(x),scalartype(y))}}(
            v => mergem!(tag(x[v], 1,neighbor_edges(v))*tag(y[v], 2,neighbor_edges(v)), Dict([tag(1,e),tag(2,e)] => e for e in neighbor_edges(v)))
            for v in vertices(mtree, root_to_leaves)
        )
    )
end

*(a::Number, x::TreeTensor) = (y = copyconvert(TreeTensor{promote_type(typeof(a),scalartype(x))}, x); y[:root] *= a; return y)
*(x::TreeTensor, a::Number) = a*x
scale!(a::Number, x::TreeTensor) = (x[:root] *= a; return x)
scale!(x::TreeTensor, a::Number) = scale!(a,x)
-(x::TreeTensor, y::TreeTensor) = x + (-1)*y

# Transposition and conjugation

for f in (:conj,:transpose,:ctranspose)
    f! = Symbol(string(f)*"!")
    @eval begin
        function Base.$f(t::TreeTensor)
            return TreeTensor(
                modetree(t),
                Dict{Tree,Tensor{scalartype(t)}}(
                    v => $f(tv)
                    for (v,tv) in t
                )
            )
        end
        function Base.$f!(t::TreeTensor)
            for v in keys(t)
                t[v] = $f!(t[v])
            end
            return t
        end
    end
end

# Orthogonalisation and truncation

export orthogonalize!, orthogonalize
function orthogonalize!(x::TreeTensor)
    for (v,p) in edges(x, leaves_to_root)
        x[v],c = qr(x[v], PairSet(v,p))
        x[p] = c*x[p]
    end
    return x
end
orthogonalize(x::TreeTensor) = orthogonalize!(copy(x))

export truncate!, truncate
function truncate!(x::TreeTensor, rank)
    orthogonalize!(x)
    s = Dict{PairSet{Tree}, Tensor{real(scalartype(x))}}()
    for (v,p) in edges_with_root(x, root_to_leaves)
        for u in children(v,p)
            e = PairSet(u,v)
            b,s[e],d = svd(x[v], e, maxrank())
            x[u] = scale!(s[e],d)*x[u]
            x[v] = scale!(b,s[e])
        end
        x[v] = resize(x[v], Dict(e => rank(s[e].data) for e in neighbor_edges(v)))
        scale!(x[v], [resize(1./s[e], Dict(e => msize(x[v],e))) for e in child_edges(v,p)]...)
    end
    return x
end
truncate(x::TreeTensor, rank) = truncate!(copy(x), rank)

# Contracted subtrees

function contracted_subtree(v::Tree,p::Tree, c::TreeTensor, args::TreeTensor...)
    cv = retag!(tag!(conj(args[1][v]), 1, neighbor_edges(v)), :C => :_, :R => :C)
    for u in children(v,p) cv *= c[u] end
    for i = 2:length(args)-1 cv *= tag(args[i][v], i, neighbor_edges(v)) end
    cv *= retag!(tag(args[end][v], length(args), neighbor_edges(v)), :C => :_)
    return cv
end

function contracted_subtrees(args::TreeTensor...)
    c = TreeTensor(
        modetree(args[1]),
        TensorDict{scalartype(args[1])}()
    )
    for (v,p) in edges(c, leaves_to_root)
        c[v] = contracted_subtree(v,p,c,args...)
    end
    return c
end

# Norm and dot

export norm!
norm!(x::TreeTensor) = norm(orthogonalize!(x)[:root])
norm( x::TreeTensor) = norm(orthogonalize( x)[:root])

function dot(x::TreeTensor, y::TreeTensor)
    return scalar(
        contracted_subtree(
           tree(x),tree(x),
           contracted_subtrees(x,y),
           x,y
        )
    )
end

# Local solvers

function localmatrix(A, v, xAx)
    lA = tag(A[v], 2, neighbor_edges(v))
    for u in neighbors(v)
        lA *= xAx[u]
    end
    retag!(lA, 1 => :R, 3 => :C)
    return lA
end

function localrhs(b, v, xb)
    lb = tag(b[v], 2, neighbor_edges(v))
    for u in neighbors(v)
        lb *= xb[u]
    end
    untag!(lb, 1)
    return lb
end

function localmatvecfunc(A, v, xAx)
    Av = tag(A[v],2,neighbor_edges(v))
    c = [retag(xAx[u], 3 => :C, 1 => :R) for u in neighbors(v)]
    if length(c) == 3
        return xv -> begin
            (c[3]*((c[1]*Av)*(c[2]*xv)))
        end
    elseif length(c) == 2
        return xv -> begin
            (c[2]*(Av*(c[1]*xv)))
        end
    elseif length(c) == 1
        return xv -> begin
            (c[1]*Av*xv)
        end
    else
        error("Local matvec structure not implemented!")
    end
end

function localproblem(x,A,b, v, xAx,xb)
    modes = x[v].modes
    matvec = localmatvecfunc(A, v, xAx)
    return x[v][mlabel(modes)],
           xv -> (xv = Tensor(modes,xv); matvec(xv)[mlabel(modes)]),
           localrhs(b,v,xb)[mlabel(modes)]
end

# Local solvers

using IterativeSolvers

abstract LocalSolver

export DenseSolver
immutable DenseSolver <: LocalSolver end
function localsolve!(x,A,b, solver::DenseSolver, v, xAx, xb)
    x[v] = localmatrix(A,v,xAx)\localrhs(b,v,xb)
end

export GMRES
immutable GMRES <: LocalSolver
    tol::Float64
    maxiter::Int
    restart::Int
    warn::Bool
end
GMRES(; tol = sqrt(eps()), maxiter = 1, restart = 20, warn = false) = GMRES(tol, maxiter, restart, warn)
function localsolve!(x,A,b, solver::GMRES, v, xAx,xb)
    _,hist = gmres!(
        localproblem(x,A,b, v, xAx,xb)...;
        tol = solver.tol,
        maxiter = solver.maxiter,
        restart = solver.restart
    )
    if !hist.isconverged && solver.warn
        warn(
            "Local GMRES iteration did not convergence.\n"*
            "  LSE size:  "*string(length(x[v]))*"\n"*
            "  Threshold: "*string(hist.threshold)*"\n"*
            "  Residuals: "*string(hist.residuals)
        )
    end
end

export CG
immutable CG <: LocalSolver
    tol::Float64
    maxiter::Int
    warn::Bool
end
CG(; tol = sqrt(eps()), maxiter = 20, warn = false) = CG(tol, maxiter, warn)
function localsolve!(x,A,b, solver::CG, v, xAx,xb)
    _,hist = cg!(
        localmatvec(x,A,b, v, xAx,xb)...;
        tol = solver.tol,
        maxiter = solver.maxiter
    )
    if !hist.isconverged && solver.warn
        warn(
            "Local CG iteration did not convergence.\n"*
            "  LSE size:  "*string(length(x[v]))*"\n"*
            "  Threshold: "*string(hist.threshold)*"\n"*
            "  Residuals: "*string(hist.residuals)
        )
    end
end

# ALS linear solver

export als_solve!
function als_solve!(
    x,A,b, solver = GMRES();
    maxiter::Int = 20, tol = sqrt(eps(real(scalartype(x))))
)
    updatenorms = zeros(real(scalartype(x)), maxiter)

    orthogonalize!(x)
    xAx = contracted_subtrees(x,A,x)
    xb = contracted_subtrees(x,b)
    for i = 1:maxiter
        for (u,v) in edges(x, both_ways)
            # Solve
            xu = copy(x[u])
            localsolve!(x,A,b, solver, u, xAx,xb)
            updatenorms[i] = max(updatenorms[i], norm(xu - x[u]))

            # Orthogonalize
            x[u],r = qr(x[u],PairSet(u,v))
            x[v] = r*x[v]

            # Compute contracted subtrees
            xAx[u] = contracted_subtree(u,v, xAx, x,A,x)
            xb[u] = contracted_subtree(u,v, xb, x,b)
        end

        # Convergence check
        if updatenorms[i] < tol*norm(x[:root])
            return x, ConvergenceHistory(true, tol*norm(x[:root]), i, updatenorms[1:i])
        end
    end
    return x, ConvergenceHistory(false, tol*norm(x[:root]), maxiter, updatenorms)
end

# ALS sum

export als_sum!
function als_sum!{T}(
    y,x::AbstractVector{TreeTensor{T}};
    maxiter::Int = 20, tol = sqrt(eps(real(T)))
)
    updatenorms = zeros(real(T), maxiter)

    orthogonalize!(y)
    yx = [contracted_subtrees(y,x) for x in x]
    for i = 1:maxiter
        for (u,v) in edges(y, both_ways)
            # Project
            yunew = sum([localrhs(x, u, yx) for (x,yx) in zip(x,yx)])
            updatenorms[i] = max(updatenorms[i], norm(yunew - y[u]))
            y[u] = yunew

            # Orthogonalize
            y[u],r = qr(y[u],PairSet(u,v))
            y[v] = r*y[v]

            # Compute contracted subtrees
            for (x,yx) in zip(x,yx)
                local x,yx
                yx[u] = contracted_subtree(u,v, yx, y,x)
            end
        end

        # Convergence check
        if updatenorms[i] < tol*norm(y[:root])
            return y, ConvergenceHistory(true, tol*norm(y[:root]), i, updatenorms[1:i])
        end
    end
    return y, ConvergenceHistory(false, tol*norm(y[:root]), maxiter, updatenorms)
end

# ALS AXPY

export als_axpy!
function als_axpy!(
    z,A,x,y;
    maxiter::Int = 20, tol = sqrt(eps(real(scalartype(x))))
)
    updatenorms = zeros(real(scalartype(x)), maxiter)

    orthogonalize!(z)
    zAx = contracted_subtrees(z,A,x)
    zy = contracted_subtrees(z,y)
    for i = 1:maxiter
        for (u,v) in edges(z, both_ways)
            # Project
            zunew = localmatvecfunc(A, u, zAx)(x[u]) + localrhs(y, u, zy)
            updatenorms[i] = max(updatenorms[i], norm(zunew - z[u]))
            z[u] = zunew

            # Orthogonalize
            z[u],r = qr(z[u],PairSet(u,v))
            z[v] = r*z[v]

            # Compute contracted subtrees
            zAx[u] = contracted_subtree(u,v, zAx, z,A,x)
            zy[u] = contracted_subtree(u,v, zy, z,y)
        end

        # Convergence check
        if updatenorms[i] < tol*norm(z[:root])
            return z, ConvergenceHistory(true, tol*norm(z[:root]), i, updatenorms[1:i])
        end
    end
    return z, ConvergenceHistory(false, tol*norm(z[:root]), maxiter, updatenorms)
end

## ALSSD linear solver
#
#function alssd_solve!(
#    x,A,b, solver = GMRES();
#    residualrank = 4, maxiter::Int = 20, tol = sqrt(eps(real(scalartype(x))))
#)
#    residualnorms = zeros(real(scalartype(x)), maxiter)
#    tol *= norm(b)
#
#    z = rand(x, residualrank)
#    for i = 1:maxiter
#        als_axpy!(z, -A,x,b; maxiter = 1)
#
#        residualnorms[i] = norm!(z)
#        if residualnorms[i] < tol
#            return x,ConvergenceHistory(true, tol, i, residualnorms[1:i])
#        end
#
#        x += z
#        als_solve!(x,A,b, solver; maxiter = 1)
#        truncate!(x, adaptive(tol))
#    end
#    return x
#end

## ALS operator inversion
#
#function als_inv!(
#    X,A, solver = GMRES();
#    maxiter::Int = 20, tol = sqrt(eps(real(scalartype(X))))
#)
#    updatenorms = zeros(real(scalartype(X)), maxiter)
#
#    orthogonalize!(X)
#    XAX = contracted_subtrees(X,A,X)
#    XAXt = contracted_subtrees(X,A,X)
#    xb = contracted_subtrees(x,b)
#    for i = 1:maxiter
#        for (u,v) in edges(x, both_ways)
#            # Solve
#            xu = copy(x[u])
#            localsolve!(x,A,b, solver, u, xAx,xb)
#            updatenorms[i] = max(updatenorms[i], norm(xu - x[u]))
#
#            # Orthogonalize
#            x[u],r = qr(x[u],PairSet(u,v))
#            x[v] = r*x[v]
#
#            # Compute contracted subtrees
#            xAx[u] = contracted_subtree(u,v, xAx, x,A,x)
#            xb[u] = contracted_subtree(u,v, xb, x,b)
#        end
#
#        # Convergence check
#        if updatenorms[i] < tol*norm(x[:root])
#            return x, ConvergenceHistory(true, tol*norm(x[:root]), i, updatenorms[1:i])
#        end
#    end
#    return x, ConvergenceHistory(false, tol*norm(x[:root]), maxiter, updatenorms)
#end

end # module


text:
__precompile__()

using Formatting
using Plots
# pyplot(reuse=true, size=(900,300))
pyplot()
theme(:ggplot2)

using ArgParse
import ArgParse: parse_item

import TimeSeries
import Distributions
import Dates
import FracFin

function parse_commandline()
    settings = ArgParseSettings("Plot results of estimation.")

    @add_arg_table settings begin
        "infile"
        help = "input .csv file containing results of estimation"
        required = true
        "outdir"
        help = "output directory"
        required = true
        "--tfmt"
        help = "time format for parsing csv file"
        arg_type = String
        default = "yyyy-mm-ddTHH:MM:SS"
        "--ci"
        help = "size of the confident interval"
        default = 0.75
        "--intraday", "-i"
        help = "plot intraday results (if available)"
        action = :store_true
        "--verbose", "-v"
        help = "print message"
        action = :store_true
    end

    return parse_args(settings)
end

function main()
    parsed_args = parse_commandline()
    # println("Parsed args:")
    # for (arg,val) in parsed_args
    #     println("  $arg  =>  $val")
    # end

    # recover global parameters and options
    infile = parsed_args["infile"]
    outdir = parsed_args["outdir"]

    tfmt = parsed_args["tfmt"]  # time format
    intraday = parsed_args["intraday"]  # intraday estimation
    ci = parsed_args["ci"]  # size of the confident interval
    verbose = parsed_args["verbose"]  # print messages

    # load data
    if verbose
        printfmtln("Loading file {}...", infile)
    end
    toto = TimeSeries.readtimearray(infile, format=tfmt, delim=',')
    X0 = toto[:Log_Price]
    Ht0 = toto[:Hurst]
    Hs0 = toto[:Hurst_std]
    σt0 = toto[:σ]
    σs0 = toto[:σ_std]

    # make the output folder
    try
        mkpath(outdir)
    catch
    end

    # plot overall estimate
    q = Distributions.quantile(Distributions.Normal(), (1+ci)/2)  # quantile

    fig1 = plot(X0, ylabel="Log Price", label="")
    fig2 = plot(Ht0, ribbon=(2*Hs0, 2*Hs0), ylim=[0,1], shape=:circle, ylabel="Hurst", label="")
    hline!(fig2, 0.5, w=3, color=:red, label="")
    fig3 = plot(σt0, ribbon=(q*σs0, q*σs0), shape=:circle, ylabel="σ", label="")
    plot(fig1, fig2, fig3, layout=(3,1), size=(1200,900), legend=true)
    savefig(format("{}/estimate.pdf",outdir))

    if intraday
        sdata = FracFin.window_split_timearray(toto, Dates.Hour(24); fillmode=:o, endpoint=true)

        for Y0 in sdata
            T0 = TimeSeries.timestamp(Y0)
            fig1 = plot(Y0[:Log_Price], ylabel="Log Price", label="")
            fig2 = plot(Y0[:Hurst], shape=:circle, ylabel="Hurst", label="")
            fig3 = plot(Y0[:σ], shape=:circle, ylabel="σ", label="")
            fig4 = plot(Y0[:Hurst], ylim=[0,1], shape=:circle, ylabel="Hurst", label="")
            hline!(fig4, 0.5, w=3, color=:red, label="")
            plot(fig1, fig2, fig3, fig4, layout=(4,1), size=(1200,900), xticks=T0[1]:Dates.Hour(1):T0[end], legend=true)
            savefig(format("{}/{}.pdf",outdir, Dates.Date(T0[1])))
            # TimeSeries.writetimearray(Yt, format("{}/csv/{}.csv",outdir, day))
        end
    end

    if verbose
        printfmtln("Outputs saved in {}", outdir)
    end
end

main()

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
function CompiledBasicBlock(b::BasicBlock)
    return CompiledBasicBlock(b.instructions.entries, b.inputs, b.mutable_locals, b.immutable_locals, b.outputs)
end

@inline input_variables(bb::CompiledBasicBlock) = bb.inputs
@inline output_variables(bb::CompiledBasicBlock) = bb.outputs
@inline is_input(bb::CompiledBasicBlock, var::ImmutableVariable) = var ∈ bb.inputs
@inline is_input(bb::CompiledBasicBlock, var::MutableVariable) = false
@inline is_output(bb::CompiledBasicBlock, var::Variable) = var ∈ bb.outputs
@inline inputindex(bb::CompiledBasicBlock, op::Variable) = indexof(bb.inputs, op)

@inline required_memory_size(program::CompiledBasicBlock) = length(program.instructions)
@inline required_input_values_size(program::CompiledBasicBlock) = length(input_variables(program))
@inline required_memory_size(program::BasicBlock) = length(program.instructions)
@inline required_input_values_size(program::BasicBlock) = length(input_variables(program))

"Used to compile to a function that is interpreted when executed."
struct BasicBlockInterpreter{T} <: AbstractMachine
    program::CompiledBasicBlock
    memory::Vector{T}
    input_values::Vector{T}
    runtime_bindings::Dict{MutableVariable, T}

    function BasicBlockInterpreter(
                program::CompiledBasicBlock,
                memory_buffer::Vector{T},
                input_values_buffer::Vector{T};
                auto_resize_buffers::Bool=true
            ) where {T}

        #@assert !hasbranches(b) "BasicBlockInterpreter does not support branches"

        if auto_resize_buffers
            if length(memory_buffer) < required_memory_size(program)
                resize!(memory_buffer, required_memory_size(program))
            end

            if length(input_values_buffer) < required_input_values_size(program)
                resize!(input_values_buffer, required_input_values_size(program))
            end
        else
            @assert length(memory_buffer) >= required_memory_size(program)
            @assert length(input_values_buffer) >= required_input_values_size(program)
        end

        return new{T}(program, memory_buffer, input_values_buffer, Dict{MutableVariable, T}())
    end
end

function BasicBlockInterpreter(
            b::BasicBlock,
            memory_buffer::Vector{T},
            input_values_buffer::Vector{T};
            auto_resize_buffers::Bool=true
        ) where {T}

    return BasicBlockInterpreter(CompiledBasicBlock(b), memory_buffer, input_values_buffer, auto_resize_buffers=auto_resize_buffers)
end

function BasicBlockInterpreter(b::BasicBlock; word_type::Type{T}=Any) where {T}
    return BasicBlockInterpreter(b, Vector{word_type}(undef, required_memory_size(b)), Vector{word_type}(undef, required_input_values_size(b)), auto_resize_buffers=false)
end

function BasicBlockInterpreter(b::CompiledBasicBlock; word_type::Type{T}=Any) where {T}
    return BasicBlockInterpreter(b, Vector{word_type}(undef, required_memory_size(b)), Vector{word_type}(undef, required_input_values_size(b)), auto_resize_buffers=false)
end

compile(::Type{BasicBlockInterpreter}, program::Program) = BasicBlockInterpreter(program)
compile(::Type{BasicBlockInterpreter{T}}, program::Program) where {T} = BasicBlockInterpreter(program, word_type=T)

function set_input!(machine::BasicBlockInterpreter, input)

    input_variables_count = length(input_variables(machine.program))

    @assert input_variables_count <= length(machine.input_values)
    @assert input_variables_count == length(input) "Expected $input_variables_count arguments. Got $(length(input))."

    for i in 1:input_variables_count
        @inbounds machine.input_values[i] = input[i]
    end

    nothing
end

# `machine(args...)` will run the program
function (machine::BasicBlockInterpreter)(input...)
    set_input!(machine, input)

    for (i, instruction) in enumerate(machine.program.instructions)
        @inbounds machine.memory[i] = execute_op(machine, instruction)
    end

    return return_values(machine)
end

function execute_op(machine::BasicBlockInterpreter, op::Const{T}) :: T where {T}
    return op.val
end

deref(::BasicBlockInterpreter, arg::Const) = arg.val
deref(machine::BasicBlockInterpreter, arg::SSAValue) = machine.memory[arg.address]
deref(machine::BasicBlockInterpreter, arg::MutableVariable) = machine.runtime_bindings[arg]

function deref(machine::BasicBlockInterpreter, arg::ImmutableVariable)
    if is_input(machine.program, arg)
        @inbounds machine.input_values[ inputindex(machine.program, arg) ]
    else
        deref(machine, machine.program.immutable_locals[arg])
    end
end

deref(machine::BasicBlockInterpreter, args::AbstractValue...) = map(ssa -> deref(machine, ssa), args)

execute_op(machine::BasicBlockInterpreter, op::PureInstruction) = execute_op(machine, op.call)
execute_op(machine::BasicBlockInterpreter, op::ImpureInstruction) = execute_op(machine, op.call)
execute_op(machine::BasicBlockInterpreter, op::CallUnary) = op.op(deref(machine, op.arg))
execute_op(machine::BasicBlockInterpreter, op::CallBinary) = op.op(deref(machine, op.arg1), deref(machine, op.arg2))
execute_op(machine::BasicBlockInterpreter, op::Call3Args) = op.op(deref(machine, op.arg1), deref(machine, op.arg2), deref(machine, op.arg3))
execute_op(machine::BasicBlockInterpreter, op::Call4Args) = op.op(deref(machine, op.arg1), deref(machine, op.arg2), deref(machine, op.arg3), deref(machine, op.arg4))
execute_op(machine::BasicBlockInterpreter, op::Call5Args) = op.op(deref(machine, op.arg1), deref(machine, op.arg2), deref(machine, op.arg3), deref(machine, op.arg4), deref(machine, op.arg5))
execute_op(machine::BasicBlockInterpreter, op::Call6Args) = op.op(deref(machine, op.arg1), deref(machine, op.arg2), deref(machine, op.arg3), deref(machine, op.arg4), deref(machine, op.arg5), deref(machine, op.arg6))

function execute_op(::BasicBlockInterpreter, ::Assignment{V}) where {V<:ImmutableVariable}
    # no-op
    nothing
end

function execute_op(machine::BasicBlockInterpreter, assignment::Assignment{V}) where {V<:MutableVariable}
    runtime_bind!(machine, assignment.lhs, deref(machine, assignment.rhs))
end

function runtime_bind!(machine::BasicBlockInterpreter, var::MutableVariable, val::Any)
    machine.runtime_bindings[var] = val
end

eachvariable(machine::BasicBlockInterpreter) = eachvariable(machine.program)

function return_values(machine::BasicBlockInterpreter)
    out_vars = output_variables(machine.program)

    if isempty(out_vars)
        return nothing
    elseif length(out_vars) == 1
        return deref(machine, out_vars[1])
    else
        return Tuple( deref(machine, variable) for variable in output_variables(machine.program) )
    end
end


text:
struct Cauchy <: Distribution{Float64} end

"""
    cauchy(x0::Real, gamma::Real)

Sample a `Float64` value from a Cauchy distribution with location x0 and scale gamma.
"""
const cauchy = Cauchy()

function logpdf(::Cauchy, x::Real, x0::Real, gamma::Real)
    return Distributions.logpdf(Distributions.Cauchy(x0, gamma), x)
end

function logpdf_grad(::Cauchy, x::Real, x0::Real, gamma::Real)
    x_x0 = x - x0
    x_x0_sq = x_x0^2
    gamma_sq = gamma^2
    deriv_x0 =  2 * x_x0 / (gamma_sq + x_x0_sq)
    deriv_x = - deriv_x0
    deriv_gamma = (x_x0_sq - gamma_sq) / (gamma * (gamma_sq + x_x0_sq))
    (deriv_x, deriv_x0, deriv_gamma)
end

is_discrete(::Cauchy) = false

random(::Cauchy, x0::Real, gamma::Real) = rand(Distributions.Cauchy(x0, gamma))

(::Cauchy)(x0::Real, gamma::Real) = random(Cauchy(), x0, gamma)

has_output_grad(::Cauchy) = true
has_argument_grads(::Cauchy) = (true, true)

export cauchy

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 60941eaa-1aea-11eb-1277-97b991548781
using Revise,PlutoUI,ExtendableGrids

# ╔═╡ df8aefc0-5273-11eb-1a67-3552732deae0
using PyPlot

# ╔═╡ 08fa72d6-5274-11eb-03bd-ef885fdd216d
P=PyPlot

# ╔═╡ e5cec5c8-5273-11eb-0e0b-03960befcc71
X=collect(0:0.5:10)

# ╔═╡ f088d080-5273-11eb-2085-e5d4734512ff
g1=simplexgrid(X)

# ╔═╡ ff091d22-5273-11eb-2645-89630b592fa0
gridplot(g1,Plotter=P,resolution=(500,500))

# ╔═╡ 2dd685fe-5274-11eb-173c-6f2ccaf41aae
gridplot(g1,map(x->sin(x),g1),Plotter=P)

# ╔═╡ 8d1857ae-5274-11eb-265d-09604d9dba94
g2=simplexgrid(X,X)

# ╔═╡ 9526886c-5274-11eb-33a6-d95a16547472
gridplot(g2,Plotter=P)

# ╔═╡ ae81ae16-5274-11eb-36e6-a1dd77e51abb
gridplot(g2,map((x,y)->(sin(x)*exp(-(y-5)^2)),g2),Plotter=P,colormap=:hot)

# ╔═╡ eb380f5a-5274-11eb-2a7a-87c75706a2ae
p=GridPlotContext(Plotter=P,layout=(1,2),fignum=2,resolution=(800,400))

# ╔═╡ 0114cb24-5275-11eb-3013-0f3b170f25de
begin
	gridplot!(p[1,1],g1,map(x->sin(x),g1),title="1D")
	gridplot!(p[1,2],g2,map((x,y)->(sin(x)*exp(-(y-5)^2)),g2),colormap=:hot,title="2D")
	reveal(p)
end

# ╔═╡ 677d67f4-5275-11eb-17d6-1954c7cd0131
g3=simplexgrid(X,X,X)

# ╔═╡ aaab7994-5275-11eb-159b-fdd7f6bee8c6
@bind zplane Slider(0:0.51:10,default=4)

# ╔═╡ 6dbcc196-5275-11eb-0bd5-2941757247b3
gridplot(g3,Plotter=P,zplane=zplane)

# ╔═╡ d6f7eb72-5275-11eb-2fa5-8b7b03a4743b
@bind flevel Slider(0:0.2:20,default=0.5)

# ╔═╡ 77e28b4c-5275-11eb-199f-253cae56c928
gridplot(g3,map((x,y,z)->(sin(x)*exp(0.1*y)*z),g3),Plotter=P,flevel=flevel,zplane=zplane,colormap=:brg)

# ╔═╡ Cell order:
# ╠═60941eaa-1aea-11eb-1277-97b991548781
# ╠═df8aefc0-5273-11eb-1a67-3552732deae0
# ╠═08fa72d6-5274-11eb-03bd-ef885fdd216d
# ╠═e5cec5c8-5273-11eb-0e0b-03960befcc71
# ╠═f088d080-5273-11eb-2085-e5d4734512ff
# ╠═ff091d22-5273-11eb-2645-89630b592fa0
# ╠═2dd685fe-5274-11eb-173c-6f2ccaf41aae
# ╠═8d1857ae-5274-11eb-265d-09604d9dba94
# ╠═9526886c-5274-11eb-33a6-d95a16547472
# ╠═ae81ae16-5274-11eb-36e6-a1dd77e51abb
# ╠═eb380f5a-5274-11eb-2a7a-87c75706a2ae
# ╠═0114cb24-5275-11eb-3013-0f3b170f25de
# ╠═677d67f4-5275-11eb-17d6-1954c7cd0131
# ╠═6dbcc196-5275-11eb-0bd5-2941757247b3
# ╠═aaab7994-5275-11eb-159b-fdd7f6bee8c6
# ╠═77e28b4c-5275-11eb-199f-253cae56c928
# ╠═d6f7eb72-5275-11eb-2fa5-8b7b03a4743b


text:
@testset "Plane Waves" begin
    @test abs(gpw(0,π/8)(0.5, π/4)) ≈ 1 atol = 1e-14
end

max_stars_repo_name: EmileParolin/STrAW

max_issues_repo_name: EmileParolin/STrAW

max_forks_repo_name: EmileParolin/STrAW


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SLEEF_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SLEEF")
JLLWrappers.@generate_main_file("SLEEF", UUID("63e82ce6-3d80-5af4-a84c-484b71ab98bb"))
end  # module SLEEF_jll


text:
function insertionsort!(v)
    for (i, key) in enumerate(v)

        while i > 1 && v[i-1] > key
            v[i] = v[i-1]
            i -= 1
        end
        v[i] = key
    end
    return v
end

v = [144 89 4 9 95 12 86 25]
println("unsorted array = ",v)

insertionsort!(v)

println("insertion sorted array = ", v)


text:
import DomainSets: convert_space, restrict_space, IntegerSpace, RationalSpace, RealSpace, ComplexSpace,
                VectorSpace, AnySpace, EmptySpace, FullSpace, GeometricSpace, promote_space_type, promote_space

# test_spaces.jl

######################################
# Helper functions and generic tests
######################################

nonzero_element(::Type{T}) where {T <: Number} = one(T)
nonzero_element(::Type{SVector{N,T}}) where {N,T<:Number} = ones(SVector{N,T})
nonzero_element(::Type{Tuple{T,S}}) where {T,S} = (nonzero_element(T),nonzero_element(S))
nonzero_element(::Type{Tuple{T,S,U}}) where {T,S,U} = (nonzero_element(T),nonzero_element(S),nonzero_element(U))
nonzero_element(::Type{Tuple{T,S,U,V}}) where {T,S,U,V} = (nonzero_element(T),nonzero_element(S),nonzero_element(U),nonzero_element(V))

# Generic tests for isomorphisms
function test_isomorphism(A, B)
    @test isomorphic(A(), B())
    @test isomorphic(B, A)
    test_embedding(A, B)
    test_embedding(B, A)
    T = eltype(A)
    S = eltype(B)
    x = nonzero_element(T)
    y = nonzero_element(S)
    # Is conversion invertible?
    @test convert_space(A, convert_space(B, x)) == x
    @test convert_space(B, convert_space(A, y)) == y
end

# Generic tests for embeddings
function test_embedding(A, B)
    @test embedded(A(), B())
    T = eltype(A)
    S = eltype(B)
    x = nonzero_element(T)
    y = nonzero_element(S)
    # Does conversion work? We verify if it yields the right type.
    @test typeof(convert_space(B, x)) == S
    # And we test the left inverse
    @test restrict_space(A, convert_space(B, x)) == x
end

##################
# Specific tests
##################

@testset "Spaces" begin
    @testset "basic spaces" begin
        Z = IntegerSpace{Int}
        Q = RationalSpace{Int}
        R = RealSpace{Float64}
        C = ComplexSpace{Float64}
        R1 = VectorSpace{1,Float64}
        R2 = VectorSpace{2,Float64}
        R3 = VectorSpace{3,Float64}
        @test ℤ == Z
        @test ℚ == Q
        @test ℝ == R
        @test ℝ1 == R1
        @test ℝ2 == R2
        @test ℝ3 == R3

        # The basic spaces are not isomorphic to each other
        @test !isomorphic(Z, Q)
        @test !isomorphic(Q, R)
        @test !isomorphic(R1, R2)
        @test !isomorphic(R2, R3)
        test_embedding(Z, Q)
        test_embedding(Z, R)
        test_embedding(Z, C)
        test_embedding(Q, R)
        test_embedding(Q, C)
        test_embedding(R, C)
        test_embedding(R, R1)
        test_embedding(R1, R2)
        test_embedding(R1, R3)

        @test embedded(AnySpace, AnySpace)
        @test embedded(Z, AnySpace)
        @test !embedded(AnySpace, Z)

        # Isomorphism between T and SVector{1,T}
        test_isomorphism(Z, VectorSpace{1,Int})
        test_isomorphism(Q, VectorSpace{1,Rational{Int}})
        test_isomorphism(R, R1)

        # Isomorphism between C and R2
        test_isomorphism(C, R2)

        # Some tests with BigFloat and BigInt
        @test widen(Z) == IntegerSpace{widen(Int)}
        @test widen(Q) == RationalSpace{widen(Int)}
        @test widen(C) == ComplexSpace{widen(Float64)}
        @test widen(R) == RealSpace{widen(Float64)}
        @test widen(R3) == VectorSpace{3,widen(Float64)}
        test_embedding(Z, widen(Z))
        test_embedding(Q, widen(Q))
        test_embedding(C, widen(C))
        test_embedding(R, widen(R))
        test_embedding(R1, widen(R3))

        test_isomorphism(widen(R), widen(R1))
        test_isomorphism(widen(C), widen(R2))

        @testset "Combination of isomorphism and widening" begin
            test_embedding(R, widen(R1))
            test_embedding(R1, widen(R))
            test_embedding(C, widen(R2))
            test_embedding(R2, widen(C))
            @test !embedded(widen(R), R1)
            @test !embedded(widen(C), R2)
            @test !embedded(widen(R1), R2)
            @test !embedded(widen(R1), R3)
        end

        @testset "Some promotions" begin
            @test promote_space() == ()
            @test typeof(promote_space(1.0)) == Tuple{typeof(1.0)}
            @test typeof(promote_space(1, 1.0)) == Tuple{typeof(1.0),typeof(1.0)}
        end

        @test promote_space_type(RealSpace{Float64}, RealSpace{BigFloat}) == RealSpace{BigFloat}
        @test promote_space_type(VectorSpace{1,BigFloat}, VectorSpace{2,Float64}) == VectorSpace{2,BigFloat}
        @test promote_space_type(VectorSpace{1,Float64}, VectorSpace{2,BigFloat}) == VectorSpace{2,BigFloat}

        # Preference for scalars over 1D vectors
        @test promote_space_type(RealSpace{Float64}, VectorSpace{1,Float64}) == RealSpace{Float64}
        @test promote_space_type(VectorSpace{1,Float64}, RealSpace{Float64}) == RealSpace{Float64}
        # - even when one space is bigger
        @test promote_space_type(RealSpace{Float64}, VectorSpace{1,BigFloat}) == RealSpace{BigFloat}
        @test promote_space_type(RealSpace{BigFloat}, VectorSpace{1,Float64}) == RealSpace{BigFloat}

        # Preference for 2D vectors over complex numbers
        @test promote_space_type(ComplexSpace{Float64}, VectorSpace{2,Float64}) == VectorSpace{2,Float64}
        @test promote_space_type(VectorSpace{2,Float64}, ComplexSpace{Float64}) == VectorSpace{2,Float64}
        # - except when one space is bigger
        @test promote_space_type(ComplexSpace{BigFloat}, VectorSpace{2,Float64}) == VectorSpace{2,BigFloat}
        @test promote_space_type(ComplexSpace{Float64}, VectorSpace{2,BigFloat}) == VectorSpace{2,BigFloat}

        @test promote_space_type(AnySpace, VectorSpace{2,Float64}) == AnySpace
        @test promote_space_type(ComplexSpace{Float64}, AnySpace) == AnySpace
        @test promote_space_type(AnySpace, AnySpace) == AnySpace

        # some basic tests

        @test_throws InexactError convert_space(RealSpace{Float64}, RealSpace)
        @test_throws InexactError restrict_space(RealSpace{Float64}, RealSpace)

        @test DomainSets._promote_via_embedding_reduction(AnySpace, AnySpace, AnySpace, ℤ, 1) == AnySpace
        @test DomainSets._promote_via_embedding_reduction(AnySpace, AnySpace, 1, AnySpace, ℤ) == AnySpace

        @test !(zero(Z) ∈ R())
        @test !(zero(Z) ∈ R())
        @test zero(Z) ∈ Z

        @test origin(Z()) == zero(Z)
        @test typeof(origin(Z())) == typeof(zero(Z))

        @test spaceof(1) == ℤ
        @test superspaceof(1) == AnySpace

        @test issubspace(Z, Z)
        @test issubspace(Z, AnySpace)
        @test !issubspace(AnySpace, Z)

        @test DomainSets.result(DomainSets.isomorphism_reduction_result(GeometricSpace{Float64}, GeometricSpace{Float64}, GeometricSpace{Float64}, GeometricSpace{Float64}, GeometricSpace{Float64}, GeometricSpace{Float64}))
        @test DomainSets.result(DomainSets.isomorphism_reduction_result(GeometricSpace{Float64}, GeometricSpace{Float64}, GeometricSpace{SVector{2,Float64}}, GeometricSpace{ComplexF64}, GeometricSpace{ComplexF64}, GeometricSpace{SVector{2,Float64}}))
        @test DomainSets.result(DomainSets.isomorphism_reduction_result(GeometricSpace{Float64}, GeometricSpace{Float64}, GeometricSpace{ComplexF64}, GeometricSpace{SVector{2,Float64}}, GeometricSpace{ComplexF64}, GeometricSpace{SVector{2,Float64}}))
        @test_throws ErrorException  DomainSets.isomorphism_reduction_result(GeometricSpace{Float64}, GeometricSpace{Float64}, GeometricSpace{Float64}, GeometricSpace{Float64}, GeometricSpace{SVector{1,Float64}}, GeometricSpace{Float64})
    end

    @testset "product spaces" begin
        T = Float64
        z = T(0)
        Z = IntegerSpace{Int}
        R = RealSpace{T}
        P = Z × R
        @test eltype(P) == Tuple{Int,T}

        R2 = R × R
        @test eltype(R2) == Tuple{T,T}
        @test zero(R2) == (z,z)
        test_isomorphism(R2, VectorSpace{2,T})

        R3 = cartesianproduct(R, R, R)
        @test eltype(R3) == Tuple{T,T,T}
        @test zero(R3) == (z,z,z)
        test_isomorphism(R3, VectorSpace{3,T})

        R4 = cartesianproduct(R, R, R, R)
        @test eltype(R4) == Tuple{T,T,T,T}
        @test zero(R4) == (z,z,z,z)
        test_isomorphism(R4, VectorSpace{4,T})

        RR = R2 × R
        @test eltype(RR) == Tuple{Tuple{T,T},T}
        test_isomorphism(RR, VectorSpace{3,T})

        @test R == cartesianproduct(R)

        @test cartesianproduct() == nothing
        @test eltype(cross(R())) == T
        @test eltype(cross(R(),R())) == NTuple{2,T}
        @test eltype(cross(R(),R(),R())) == NTuple{3,T}
        @test eltype(cross(R(),R(),R(),R())) == NTuple{4,T}
    end
end


text:
module Blue

include("api.jl")

end # module


text:
ndgrid(v::AbstractVector) = (copy(v),)

function ndgrid(v1::AbstractVector{T}, v2::AbstractVector{T}) where {T}
    m, n = length(v1), length(v2)
    v1 = reshape(v1, m, 1)
    v2 = reshape(v2, 1, n)
    (repeat(v1, 1, n), repeat(v2, m, 1))
end

function ndgrid_fill(a, v, s, snext)
    for j = 1:length(a)
        a[j] = v[div(rem(j - 1, snext), s) + 1]
    end
end

function ndgrid(vs::AbstractVector{T}...) where {T}
    n = length(vs)
    sz = map(length, vs)
    out = ntuple(i -> Array{T}(undef, sz), n)
    s = 1
    for i = 1:n
        a = out[i]::Array
        v = vs[i]
        snext = s * size(a, i)
        ndgrid_fill(a, v, s, snext)
        s = snext
    end
    out
end

function xchenodes(N)
    # Extended Chebyshev nodes on the interval [-1, 1]

    x = 0.0;
    if N>0
        n = N+1;
        k = Array(1:n);
        x = -cos.((2*k .- 1)*pi/(2*n))/cos(pi/(2*n));
        #x = 0.5 .+ 0.5*x;
    end
    return x
end

function tensornodes(N::Vector{Int64})

    dim = length(N)
    if dim==1
        x = xchenodes(N[1])
        return x
    end

    n = prod(N .+ 1)
    x1 = xchenodes(N[1])
    x2 = xchenodes(N[2])

    if dim==2
        xt = ndgrid(x1, x2)
    elseif dim==3
        x3 = xchenodes(N[3])
        xt = ndgrid(x1, x2, x3)
    elseif dim==4
        x3 = xchenodes(N[3])
        x4 = xchenodes(N[4])
        xt = ndgrid(x1, x2, x3, x4)
    elseif dim==5
        x3 = xchenodes(N[3])
        x4 = xchenodes(N[4])
        x5 = xchenodes(N[5])
        xt = ndgrid(x1, x2, x3, x4, x5)
    elseif dim==6
        x3 = xchenodes(N[3])
        x4 = xchenodes(N[4])
        x5 = xchenodes(N[5])
        x6 = xchenodes(N[6])
        xt = ndgrid(x1, x2, x3, x4, x5, x6)
    elseif dim==7
        xi = Array{Any}(undef, dim)
        for i = 1:dim
            xi[i] = xchenodes(N[i])
        end
        xt = ndgrid(xi[1], xi[2], xi[3], xi[4], xi[5], xi[6], xi[7])
    elseif dim==8
        xi = Array{Any}(undef, dim)
        for i = 1:dim
            xi[i] = xchenodes(N[i])
        end
        xt = ndgrid(xi[1], xi[2], xi[3], xi[4], xi[5], xi[6], xi[7], xi[8])
    else
        error("Dimension must not exceed 8")
    end

    x = zeros(n, dim)
    for i = 1:dim
        x[:,i] = xt[i][:]
    end

    return x
end

function legendrepoly(x::Vector{Float64}, N::Int64)

    m = length(x)
    k = N + 1
    p = ones(m,k)

    if N > 0
        p[:,2] = x
    end

    if N > 1
        for i = 3:k
            p[:,i] = ((2*(i-2)+1)/(i-1))*(x.*p[:,i-1]) - ((i-2)/(i-1))*p[:,i-2]
        end
    end

    return p
end

function legendrepolyder(x::Vector{Float64}, N::Int64)

    m = length(x)
    k = N + 1
    p = ones(m,k)
    dp = zeros(m, k)

    if N > 0
        p[:,2] = x
        dp[:,2] .= 1.0
    end

    if N > 1
        for i = 3:k
            p[:,i] = ((2*(i-2)+1)/(i-1))*(x.*p[:,i-1]) - ((i-2)/(i-1))*p[:,i-2]
            dp[:,i] = ((2*(i-2)+1)/(i-1))*(p[:,i-1] + x.*dp[:,i-1]) - ((i-2)/(i-1))*dp[:,i-2]
        end
    end

    return p, dp
end

function tensorpoly(x::Matrix{Float64}, N::Vector{Int64})

    n,dim = size(x)
    if dim==1
        p = legendrepoly(x[:,1], N[1])
        return p
     end

    p = zeros(n, prod(N .+ 1))
    p1 = legendrepoly(x[:,1], N[1])
    p2 = legendrepoly(x[:,2], N[2])
    if dim==2
        for i = 1:n
            p[i,:] = kron(p1[i,:], p2[i,:])
        end
    elseif dim==3
        p3 = legendrepoly(x[:,3], N[3])
        for i = 1:n
            p[i,:] = kron(kron(p1[i,:], p2[i,:]), p3[i,:])
        end
    elseif dim==4
        p3 = legendrepoly(x[:,3], N[3])
        p4 = legendrepoly(x[:,4], N[4])
        for i = 1:n
            p[i,:] = kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), p4[i,:])
        end
    elseif dim==5
        p3 = legendrepoly(x[:,3], N[3])
        p4 = legendrepoly(x[:,4], N[4])
        p5 = legendrepoly(x[:,5], N[5])
        for i = 1:n
            p[i,:] = kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), p5[i,:])
        end
    elseif dim==6
        p3 = legendrepoly(x[:,3], N[3])
        p4 = legendrepoly(x[:,4], N[4])
        p5 = legendrepoly(x[:,5], N[5])
        p6 = legendrepoly(x[:,6], N[6])
        for i = 1:n
            p[i,:] = kron(kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), p5[i,:]), p6[:,i])
        end
    elseif dim==7
        p3 = legendrepoly(x[:,3], N[3])
        p4 = legendrepoly(x[:,4], N[4])
        p5 = legendrepoly(x[:,5], N[5])
        p6 = legendrepoly(x[:,6], N[6])
        p7 = legendrepoly(x[:,7], N[7])
        for i = 1:n
            p[i,:] = kron(kron(kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), p5[i,:]), p6[:,i]), p7[:,i])
        end
    else
        error("Dimension must not exceed 7")
    end

    return p
end

function tensorpolyder(x::Matrix{Float64}, N::Vector{Int64})

    n,dim = size(x)
    if dim==1
        p, dp = legendrepolyder(x[:,1], N[1])
        return p, dp
     end

    p = zeros(n, prod(N .+ 1))
    dp = zeros(n, dim, prod(N .+ 1))

    p1, dpx1 = legendrepolyder(x[:,1], N[1])
    p2, dpx2 = legendrepolyder(x[:,2], N[2])
    if dim==2
        for i = 1:n
            p[i,:] = kron(p1[i,:], p2[i,:])
            dp[i,1,:] = kron(dpx1[i,:], p2[i,:])
            dp[i,2,:] = kron(p1[i,:], dpx2[i,:])
        end
    elseif dim==3
        p3, dpx3 = legendrepolyder(x[:,3], N[3])
        for i = 1:n
            p[i,:] = kron(kron(p1[i,:], p2[i,:]), p3[i,:])
            dp[i,1,:] = kron(kron(dpx1[i,:], p2[i,:]), p3[i,:])
            dp[i,2,:] = kron(kron(p1[i,:], dpx2[i,:]), p3[i,:])
            dp[i,3,:] = kron(kron(p1[i,:], p2[i,:]), dpx3[i,:])
        end
    elseif dim==4
        p3, dpx3 = legendrepolyder(x[:,3], N[3])
        p4, dpx4 = legendrepolyder(x[:,4], N[4])
        for i = 1:n
            p[i,:] = kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), p4[i,:])
            dp[i,1,:] = kron(kron(kron(dpx1[i,:], p2[i,:]), p3[i,:]), p4[i,:])
            dp[i,2,:] = kron(kron(kron(p1[i,:], dpx2[i,:]), p3[i,:]), p4[i,:])
            dp[i,3,:] = kron(kron(kron(p1[i,:], p2[i,:]), dpx3[i,:]), p4[i,:])
            dp[i,4,:] = kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), dpx4[i,:])
        end
    elseif dim==5
        p3, dpx3 = legendrepolyder(x[:,3], N[3])
        p4, dpx4 = legendrepolyder(x[:,4], N[4])
        p5, dpx5 = legendrepolyder(x[:,5], N[5])
        for i = 1:n
            p[i,:] = kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), p5[i,:])
            dp[i,1,:] = kron(kron(kron(kron(dpx1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), p5[i,:])
            dp[i,2,:] = kron(kron(kron(kron(p1[i,:], dpx2[i,:]), p3[i,:]), p4[i,:]), p5[i,:])
            dp[i,3,:] = kron(kron(kron(kron(p1[i,:], p2[i,:]), dpx3[i,:]), p4[i,:]), p5[i,:])
            dp[i,4,:] = kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), dpx4[i,:]), p5[i,:])
            dp[i,5,:] = kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), dpx5[i,:])
        end
    elseif dim==6
        p3, dpx3 = legendrepolyder(x[:,3], N[3])
        p4, dpx4 = legendrepolyder(x[:,4], N[4])
        p5, dpx5 = legendrepolyder(x[:,5], N[5])
        p6, dpx6 = legendrepolyder(x[:,6], N[6])
        for i = 1:n
            p[i,:] = kron(kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), p5[i,:]), p6[i,:])
            dp[i,1,:] = kron(kron(kron(kron(kron(dpx1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), p5[i,:]), p6[i,:])
            dp[i,2,:] = kron(kron(kron(kron(kron(p1[i,:], dpx2[i,:]), p3[i,:]), p4[i,:]), p5[i,:]), p6[i,:])
            dp[i,3,:] = kron(kron(kron(kron(kron(p1[i,:], p2[i,:]), dpx3[i,:]), p4[i,:]), p5[i,:]), p6[i,:])
            dp[i,4,:] = kron(kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), dpx4[i,:]), p5[i,:]), p6[i,:])
            dp[i,5,:] = kron(kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), dpx5[i,:]), p6[i,:])
            dp[i,6,:] = kron(kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), p5[i,:]), dpx6[i,:])
        end
    elseif dim==7
        p3, dpx3 = legendrepolyder(x[:,3], N[3])
        p4, dpx4 = legendrepolyder(x[:,4], N[4])
        p5, dpx5 = legendrepolyder(x[:,5], N[5])
        p6, dpx6 = legendrepolyder(x[:,6], N[6])
        p6, dpx6 = legendrepolyder(x[:,7], N[7])
        for i = 1:n
            p[i,:] = kron(kron(kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), p5[i,:]), p6[i,:]), p7[i,:])
            dp[i,1,:] = kron(kron(kron(kron(kron(kron(dpx1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), p5[i,:]), p6[i,:]), p7[i,:])
            dp[i,2,:] = kron(kron(kron(kron(kron(kron(p1[i,:], dpx2[i,:]), p3[i,:]), p4[i,:]), p5[i,:]), p6[i,:]), p7[i,:])
            dp[i,3,:] = kron(kron(kron(kron(kron(kron(p1[i,:], p2[i,:]), dpx3[i,:]), p4[i,:]), p5[i,:]), p6[i,:]), p7[i,:])
            dp[i,4,:] = kron(kron(kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), dpx4[i,:]), p5[i,:]), p6[i,:]), p7[i,:])
            dp[i,5,:] = kron(kron(kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), dpx5[i,:]), p6[i,:]), p7[i,:])
            dp[i,6,:] = kron(kron(kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), p5[i,:]), dpx6[i,:]), p7[i,:])
            dp[i,7,:] = kron(kron(kron(kron(kron(kron(p1[i,:], p2[i,:]), p3[i,:]), p4[i,:]), p5[i,:]), p6[i,:]), dpx7[i,:])
        end
    else
        error("Dimension must not exceed 7")
    end

    return p, dp
end

function ref2dom(x::Vector{Float64}, ymin::Float64, ymax::Float64)

    y = ymin .+ (0.5*(ymax-ymin))*(x .+ 1)
    return y
end

function ref2dom(x::Matrix{Float64}, ymin::Vector{Float64}, ymax::Vector{Float64})

    y = 0.0*x

    n = size(x,2)
    for i = 1:n
        y[:,i] = ymin[i] .+ (0.5*(ymax[i]-ymin[i]))*(x[:,i] .+ 1)
    end

    return y
end

function tensornodes(ymin::Vector{Float64}, ymax::Vector{Float64}, N::Vector{Int64})

    x = tensornodes(N)
    if length(N) == 1
        y = ref2dom(x, ymin[1], ymax[1])
        return reshape(y,(N[1]+1,1))
    end
    y = ref2dom(x, ymin, ymax)

    return y
end

function dom2ref(y::Vector{Float64}, ymin::Float64, ymax::Float64)

    x = (2/(ymax-ymin))*(y .- ymin) .- 1

    return x
end

function dom2ref(y::Matrix{Float64}, ymin::Vector{Float64}, ymax::Vector{Float64})

    x = 0.0*y

    n = size(y,2)
    for i = 1:n
        x[:,i] = (2/(ymax[i] - ymin[i]))*(y[:,i] .- ymin[i]) .- 1
    end

    return x
end

function mpolyinterp(y::Matrix{Float64}, f, ymin::Vector{Float64}, ymax::Vector{Float64}, N::Vector{Int64})

    x = dom2ref(y, ymin, ymax)
    A = tensorpoly(x, N)
    c = A\f

    return c
end

function mpolyinterp(y::Matrix{Float64}, f, yspace)

    ymin = yspace[:,1]
    ymax = yspace[:,2]
    N = Int64.(yspace[:,3])
    x = dom2ref(y, ymin, ymax)
    A = tensorpoly(x, N)
    c = A\f

    return c
end

function mpolyeval(y::Matrix{Float64}, c::Vector{Float64}, ymin::Vector{Float64}, ymax::Vector{Float64}, N::Vector{Int64})

    x = dom2ref(y, ymin, ymax)
    A = tensorpoly(x, N)
    f = A*c

    return f
end

function mpolyeval(yspace, c)

    ymin = yspace[:,1]
    ymax = yspace[:,2]
    Nin = Int64.(yspace[:,3])
    Nout = Int64.(yspace[:,4])
    y = tensornodes(ymin, ymax, Nout)
    x = dom2ref(y, ymin, ymax)
    A = tensorpoly(x, Nin)
    f = A*c

    return y, f
end

function mpolyder(y::Matrix{Float64}, c::Vector{Float64}, ymin::Vector{Float64}, ymax::Vector{Float64}, N::Vector{Int64})

    x = dom2ref(y, ymin, ymax)
    A, D = tensorpolyder(x, N)
    f = A*c

    n, dim = size(x)
    m = length(c)

    df = reshape(D,(n*dim, m))*c
    df = reshape(df, (n, dim))
    for i = 1:dim
        df[:,i] = df[:,i]*(2.0/(ymax[i]-ymin[i]))
    end

    return f, df
end

function tensorpolyinterp(f,  N)

    N = Int64.(N)
    d = length(N)

    A = Array{Any}(undef, d)
    B = Array{Any}(undef, d)
    for i = 1:d
        Ni = N[i]
        xi = tensornodes([Ni])
        A[i] = legendrepoly(xi, Ni)
        B[i] = inv(A[i])
    end

    # naive method
    D = B[1]
    for i = 2:d
        D = kron(D, B[i])
    end
    c = D*f

    # if d==1
    #     c = B[1]*f[:]
    #     return c
    # end
    # sum factorization method
    # n = N .+ 1
    # f = reshape(f, tuple(n...))
    # c = B[1]*reshape(f, (n[1], prod(n[2:end])))
    # ind = Array(1:d)
    # for i = 2:d
    #     perm = vcat([i, ind[(i+1):end], ind[1:i-1]]...)
    #     c = permutedims(reshape(c, tuple(n...)), tuple(perm...))
    #     c = B[i]*reshape(c, (n[i], prod(n[perm[2:end]])))
    # end
    # perm = [Array(2:d); 1]
    # c = permutedims(reshape(c, tuple(n...)), tuple(perm...))
    # c = c[:]

    return c
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP