text:
######### MachO.jl - An implementation of the MachO File format ################
#
## General design philosophy
#
# All methods should operate on seekable IO object. In particular this means
# that the basic interface buffers as little as possible (e.g. the MachOHandle)
# object does not contain the mach_header.
#
################################################################################

module MachO

# For printing
import Base: show, print, bytestring

# This package implements the ObjFileBase interface
import ObjFileBase
import ObjFileBase: sectionsize, sectionoffset, readheader, ObjectHandle, readmeta

# Reexports from ObjFileBase
export sectionsize, sectionoffset, readheader, readmeta

export readmeta, readheader, LoadCmds, Sections, Symbols, symname, segname,
    debugsections

############################ Data Structures ###################################

#
# Contains all constants defined in the MachO standard with values taken
# from the appropriate headers on MacOS 10.9
#
include("constants.jl")

#
# StrPack is used mostly for endianness handling
#
using StrPack

# MachOHandle

#
# Represents the actual MachO file
#
immutable MachOHandle{T<:IO} <: ObjectHandle
    # The IO object. This field is speciallized on to avoid dispatch performance
    # hits, especially when operating on an IOBuffer, which is an important
    # usecase for in-memory files
    io::T
    # position(io) of the start of the file in the io stream.
    start::Int
    # Whether or not the data is bswap'ed in memory (i.e. has different
    # endianness)
    bswapped::Bool
    # Whether or not the file is 64bit
    is64::Bool
end

endianness(oh::MachOHandle) = oh.bswapped ? :SwappedEndian : :NativeEndian

MachOHandle{T<:IO}(io::T,start::Int,bswapped::Bool,is64::Bool) =
    MachOHandle{T}(io,start,bswapped,is64)

function show(io::IO,h::MachOHandle)
    print(io,"MachO handle (")
    print(io,h.is64?"64-bit":"32-bit")
    h.bswapped && print(io,",swapped")
    print(io,")")
end

############################ Data Structures ###################################
#
# This section contains data structures as defined by the MachO specification.
# They are used below to acutally read in the data and may be ocassionally
# referenced from interface structs where this is convenient and no other
# interface exists
#
################################################################################

abstract MachOLC

@struct immutable mach_header
    magic::Uint32
    cputype::Uint32
    cpusubtype::Uint32
    filetype::Uint32
    ncmds::Uint32
    sizeofcmds::Uint32
    flags::Uint32
end

@struct immutable mach_header_64
    magic::Uint32
    cputype::Uint32
    cpusubtype::Uint32
    filetype::Uint32
    ncmds::Uint32
    sizeofcmds::Uint32
    flags::Uint32
    reserved::Uint32
end

@struct immutable load_command
    cmd::Uint32
    cmdsize::Uint32
end

# A 16 byte string, represented as a Uint128, but shown as a string
@struct immutable small_fixed_string
    string::Uint128
end

@struct immutable uuid_command <: MachOLC
    uuid::Uint128
end

@struct immutable segment_command <: MachOLC
    segname::small_fixed_string
    vmaddr::Uint32
    vmsize::Uint32
    fileoff::Uint32
    filesize::Uint32
    maxprot::Uint32
    initprot::Uint32
    nsects::Uint32
    flags::Uint32
end

@struct immutable segment_command_64 <: MachOLC
    segname::small_fixed_string
    vmaddr::Uint64
    vmsize::Uint64
    fileoff::Uint64
    filesize::Uint64
    maxprot::Cint
    initprot::Cint
    nsects::Uint32
    flags::Uint32
end

@struct immutable section <: ObjFileBase.Section{MachOHandle}
    sectname::small_fixed_string
    segname::small_fixed_string
    addr::Uint32
    size::Uint32
    offset::Uint32
    align::Uint32
    reloff::Uint32
    nreloc::Uint32
    flags::Uint32
    reserved1::Uint32
    reserved2::Uint32
end

@struct immutable section_64 <: ObjFileBase.Section{MachOHandle}
    sectname::small_fixed_string
    segname::small_fixed_string
    addr::Uint64
    size::Uint64
    offset::Uint32
    align::Uint32
    reloff::Uint32
    nreloc::Uint32
    flags::Uint32
    reserved1::Uint32
    reserved2::Uint32
end

@struct immutable symtab_command <: MachOLC
    symoff::Uint32
    nsyms::Uint32
    stroff::Uint32
    strsize::Uint32
end

@struct immutable dysymtab_command <: MachOLC
    ilocalsym::Uint32
    nlocalsym::Uint32
    iextdefsym::Uint32
    nextdefsym::Uint32
    iundefsym::Uint32
    nundefsym::Uint32
    tocoff::Uint32
    ntoc::Uint32
    modtaboff::Uint32
    nmodtab::Uint32
    extrefsymoff::Uint32
    nextrefsyms::Uint32
    indirectsymoff::Uint32
    nindirectsyms::Uint32
    extreloff::Uint32
    nextrel::Uint32
    locreloff::Uint32
    nlocrel::Uint32
end

@struct immutable version_min_macosx_command <: MachOLC
    version::Uint32
    sdk::Uint32
end

@struct immutable lc_str
    offset::Uint32
end

@struct immutable dylib
    name::lc_str
    timestamp::Uint32
    current_version::Uint32
    compatibilty::Uint32
end

@struct immutable dylib_command <: MachOLC
    dylib::dylib
end

@struct immutable dyld_info_command <: MachOLC
    rebase_off::Uint32
    rebase_size::Uint32
    bind_off::Uint32
    bind_size::Uint32
    weak_bind_off::Uint32
    weak_bind_size::Uint32
    lazy_bind_off::Uint32
    lazy_bind_size::Uint32
    export_off::Uint32
    export_size::Uint32
end

@struct immutable source_version_command <: MachOLC
    version::Uint64
end

@struct immutable linkedit_data_command <: MachOLC
    offset::Uint32
    size::Uint32
end

@struct immutable sub_framework_command <: MachOLC
    umbrella::lc_str
end

@struct immutable nlist <: ObjFileBase.SymtabEntry{MachOHandle}
    n_strx::Uint32
    n_type::Uint8
    n_sect::Uint8
    n_desc::Uint16
    n_value::Uint32
end

@struct immutable nlist_64 <: ObjFileBase.SymtabEntry{MachOHandle}
    n_strx::Uint32
    n_type::Uint8
    n_sect::Uint8
    n_desc::Uint16
    n_value::Uint64
end

@struct immutable fat_header
    nfat_arch::Uint32
end

@struct immutable fat_arch
    cputype::Uint32
    cpusubtype::Uint32
    offset::Uint32
    size::Uint32
    align::Uint32
end

########################### Printing Data Structures ###########################
#
# This prints the basic Mach-O data structures above. Where there is no good
# reason not to, the output matches that of otool.
#
################################################################################

function bytestring(x::small_fixed_string)
    a8 = reinterpret(Uint8,[x.string])
    z = findfirst(a8,0)
    ASCIIString(a8[1:(z == 0 ? length(a8) : z-1)])
end
show(io::IO,x::small_fixed_string) = show(io,bytestring(x))
print(io::IO,x::small_fixed_string) = print(io,bytestring(x))

==(x::small_fixed_string,y::String) = bytestring(x) == y
==(x::String,y::small_fixed_string) = y==x

*(a::ASCIIString,b::small_fixed_string) = a*bytestring(b)

function printfield(io::IO,string,fieldlength)
    print(io," "^max(fieldlength-length(string),0))
    print(io,string)
end

# TODO: Implement
printflags(io,flags) = nothing

function show(io::IO,h::Union(mach_header,mach_header_64))
    print(io,"Mach header (",isa(h,mach_header_64)?"64-bit":"32-bit",")\n")
    print(io,"")
    print(io,"""
          magic  cputype cpusubtype     filetype ncmds sizeofcmds    flags
    """)
    printfield(io,MAGICS[h.magic],11)
    printfield(io,CPUTYPES[h.cputype],9)
    print(io," "^11)
    # TODO: cpusubtype printing
    #printfield(io,CPUSUBTYPES[h.cpusubtype],8)
    # otool has caps == h.cpusubtype & CPU_SUBTYPE_MASK here, I don't consider
    # it all that useful. If somebody has a strong case why it's useful to have
    # that extra field, prlease submit a PR!
    printfield(io,FILETYPES[h.filetype],13)
    printfield(io,dec(h.ncmds),6)
    printfield(io,dec(h.sizeofcmds),11)
    printflags(io,h.flags)
end

const VM_PROT_READ      =   0x01 # read permission
const VM_PROT_WRITE     =   0x02 # write permission
const VM_PROT_EXECUTE   =   0x04 # execute permission

protstr(prot) = string(
    (prot & VM_PROT_READ)    > 0 ? "r" : "",
    (prot & VM_PROT_WRITE)   > 0 ? "w" : "",
    (prot & VM_PROT_EXECUTE) > 0 ? "x" : "")

function segflags(flags)
    if flags == 0
        return "(none)"
    else
        return "(unimplemented)"
    end
end

#     printentry(io,"offset",dec(l.offset))
#    printentry(io,"align","2^",l.align," (",2^l.align,")")

printentry(io::IO,header,values...) = (printfield(io,header,15);println(io," ",values...))
# Showing Load commands
function show(io::IO,l::Union(segment_command_64,segment_command))
    println(io,"Load Command (",isa(l,segment_command_64)?"SEGMENT_64":
        "SEGMENT","):")

    printentry(io,"name",l.segname)
    printentry(io,"addr",l.vmaddr)
    printentry(io,"size",l.vmsize)
    printentry(io,"fileoff",dec(l.fileoff))
    printentry(io,"filesize",dec(l.filesize))
    printentry(io,"maxprot",protstr(l.maxprot))
    printentry(io,"initprot",protstr(l.initprot))
    printentry(io,"nsects",dec(l.nsects))
    printentry(io,"flags",segflags(l.flags))
end

function sattrs(attributes)
    strings = ASCIIString[]
    for (k,v) in SECATTRS
        if attributes & k > 0
            push!(strings, v)
        end
    end
    join(strings,",")
end

function show(io::IO, s::Union(section,section_64))
    #function body
    println(io,"  Section:")

    printentry(io,"sectname",s.sectname)
    printentry(io,"segname",s.segname)
    printentry(io,"addr","0x",hex(s.addr,2*sizeof(s.addr)))
    printentry(io,"size","0x",hex(s.size,2*sizeof(s.addr)))
    printentry(io,"offset",dec(s.offset))
    printentry(io,"align","2^",s.align," (",2^s.align,")")
    printentry(io,"reloff",s.reloff)
    printentry(io,"nreloc",s.nreloc)
    printentry(io,"type",SECTYPES[s.flags&SECTION_TYPE])
    printentry(io,"attributes",sattrs(s.flags&SECTION_ATTRIBUTES))
end

function show(io::IO,l::symtab_command)
    println(io,"Load Command (SYMTAB):")
    printentry(io,"symoff",dec(l.symoff))
    printentry(io,"nsyms",dec(l.nsyms))
    printentry(io,"stroff",dec(l.stroff))
    printentry(io,"strsize",dec(l.strsize))
end

function show(io::IO,l::dysymtab_command)
    println(io,"Load Command (DYSYMTAB):")
    printentry(io,"ilocalsym",dec(l.ilocalsym))
    printentry(io,"nlocalsym",dec(l.nlocalsym))
    printentry(io,"iextdefsym",dec(l.iextdefsym))
    printentry(io,"nextdefsym",dec(l.nextdefsym))
    printentry(io,"iundefsym",dec(l.iundefsym))
    printentry(io,"nundefsym",dec(l.nundefsym))
    printentry(io,"tocoff",dec(l.tocoff))
    printentry(io,"ntoc",dec(l.ntoc))
    printentry(io,"modtaboff",dec(l.modtaboff))
    printentry(io,"nmodtab",dec(l.nmodtab))
    printentry(io,"extrefsymoff",dec(l.extrefsymoff))
    printentry(io,"nextrefsyms",dec(l.nextrefsyms))
    printentry(io,"indirectsymoff",dec(l.indirectsymoff))
    printentry(io,"nindirectsyms",dec(l.nindirectsyms))
    printentry(io,"extreloff",dec(l.extreloff))
    printentry(io,"nextrel",dec(l.nextrel))
    printentry(io,"locreloff",dec(l.locreloff))
    printentry(io,"nlocrel",dec(l.nlocrel))
end

function decodeversion(io::IO, version::Uint32)
    print(io,version>>16)
    print(io,'.')
    print(io,(version>>8)&0xFF)
    print(io,'.')
    print(io,version&0xFF)
end

function show(io::IO,l::version_min_macosx_command)
    println(io,"Load Command (VERSION_MIN_MACOSX):")
    printentry(io,"version",sprint(decodeversion,l.version))
    printentry(io,"sdk",sprint(decodeversion,l.sdk))
end

################################ Interface #####################################

import Base: show,
    # IO methods
    read, write, seek, seekstart, position, readuntil, readbytes,
    # Iteration
    start, next, done,
    # Indexing
    length, getindex

import StrPack: unpack, pack

#
# Note that this function is different from ObjFileBase.readmeta
# Constructs and initializes the MachOHandle object
#
function readmeta(io::IO,::Type{MachOHandle})
    start = position(io)
    magic = read(io,Uint32)
    if magic == MH_MAGIC
        return MachOHandle(io,start,false,false)
    elseif magic == MH_CIGAM
        return MachOHandle(io,start,true,false)
    elseif magic == MH_MAGIC_64
        return MachOHandle(io,start,false,true)
    elseif magic == MH_CIGAM_64
        return MachOHandle(io,start,true,true)
    elseif magic == FAT_CIGAM || magic == FAT_MAGIC
        return FatMachOHandle(io,start)
    else
        error("Invalid Magic ($(hex(magic)))!")
    end
end

function readloadcmd(h::MachOHandle)
    cmd = unpack(h,load_command)
    ccmd = cmd.cmd & ~LC_REQ_DYLD
    if ccmd == LC_UUID
        return (cmd,unpack(h, uuid_command))
    elseif ccmd == LC_SEGMENT
        return (cmd,unpack(h, segment_command))
    elseif ccmd == LC_SEGMENT_64
        return (cmd,unpack(h, segment_command_64))
    elseif ccmd == LC_SYMTAB
        return (cmd,unpack(h, symtab_command))
    elseif ccmd == LC_DYSYMTAB
        return (cmd,unpack(h, dysymtab_command))
    elseif ccmd == LC_VERSION_MIN_MACOSX
        return (cmd,unpack(h, version_min_macosx_command))
    elseif ccmd == LC_ID_DYLIB || ccmd == LC_LOAD_DYLIB ||
        cmd.cmd == LC_REEXPORT_DYLIB || ccmd == LC_LOAD_UPWARD_DYLIB
        return (cmd,unpack(h, dylib_command))
    elseif ccmd == LC_DYLD_INFO
        return (cmd,unpack(h, dyld_info_command))
    elseif ccmd == LC_SOURCE_VERSION
        return (cmd,unpack(h, source_version_command))
    elseif ccmd == LC_CODE_SIGNATURE || ccmd == LC_SEGMENT_SPLIT_INFO ||
            ccmd == LC_FUNCTION_STARTS || ccmd == LC_DATA_IN_CODE ||
            ccmd == LC_DYLIB_CODE_SIGN_DRS
        return (cmd,unpack(h, linkedit_data_command))
    elseif ccmd == LC_SUB_FRAMEWORK
        return (cmd,unpack(h,sub_framework_command))
    else
        error("Unimplemented load command $(LCTYPES[cmd.cmd]) (0x$(hex(cmd.cmd)))")
    end
end

## Iteration

import Base: eltype

# Iterate over the load commands of a file
# If the header has
# already been read, it may be passed to the constructor. Optionally
# the start address of the load commands may also be passed in
# passed in. If a header is passed, it is assumed that the start address
# is the current position of the stream. Otherwise, it will be determined when
# reading the header.
immutable LoadCmds
    h::MachOHandle
    start::Uint64
    ncmds::Uint32
    sizeofcmds::Uint32
end

immutable LoadCmd{T<:MachOLC}
    h::MachOHandle
    off::Uint64
    cmd::T
end

show{T}(io::IO, x::LoadCmd{T}) = (print(io,"0x",hex(x.off,8),":\n "); show(io,x.cmd); print(io,'\n'))

LoadCmd{T<:MachOLC}(h::MachOHandle, off::Uint64, cmd::T) = LoadCmd{T}(h,off,cmd)
eltype{T<:MachOLC}(::LoadCmd{T}) = T

function LoadCmds(h::MachOHandle, header = nothing, start = -1)
    if header === nothing
        header = readheader(h)
    end
    if start == -1
        start = position(h)
    end
    LoadCmds(h,start,header.ncmds,header.sizeofcmds)
end

# A tuple of the position before the current load command,
# the number of the load command and the size of the current load
# command. I.e. the next load command will be found at state[1]+state[3]
start(l::LoadCmds) = (l.start,0,0)
seek(l::LoadCmds,state) = seek(l.h,state[1]+state[3])
function next(l::LoadCmds,state)
    seek(l,state)
    cmdh,cmd = readloadcmd(l.h)
    (LoadCmd(l.h,state[1]+state[3],cmd),(state[1]+state[3],state[2]+1,cmdh.cmdsize))
end
done(l::LoadCmds,state) = state[2] >= l.ncmds

# Access to sections

typealias segment_commands Union(segment_command_64,segment_command)

immutable Sections
    h::MachOHandle
    command::segment_commands
    start::Int
    function Sections(h::MachOHandle, segment::segment_commands, start=-1)
        if start==-1
            start = position(h)
        end
        new(h,segment,start)
    end
    function Sections(segment::LoadCmd)
        if !(eltype(segment) <: segment_commands)
            error("Load command is not a segment")
        end
        start = segment.off + sizeof(eltype(segment)) + sizeof(load_command)
        new(segment.h,segment.cmd,start)
    end
end

length(s::Sections) = s.command.nsects
function getindex(s::Sections,n)
    if n < 1 || n > length(s)
        throw(BoundsError())
    end
    sT = isa(s.command,segment_command_64) ? section_64 : section
    seek(s.h,s.start + (n-1)*sizeof(sT))
    unpack(s.h, sT)
end

start(s::Sections) = 1
done(s::Sections,n) = n > length(s)
next(s::Sections,n) = (s[n],n+1)

for f in (:read,:readuntil,:write)
    @eval $(f){T<:IO}(io::MachOHandle{T},args...) = $(f)(io.io,args...)
end
readbytes{T<:IO}(io::MachOHandle{T},num::Integer) = readbytes(io.io,num)

seek{T<:IO}(io::MachOHandle{T},pos::Integer) = seek(io.io,io.start+pos)
seekstart(io::MachOHandle) = seek(io.io,io.start)
position{T<:IO}(io::MachOHandle{T}) = position(io.io)-io.start

unpack{T,ioT<:IO}(h::MachOHandle{ioT},::Type{T}) =
    unpack(h.io,T,h.bswapped ? :SwappedEndian : :NativeEndian)

pack{T,ioT<:IO}(h::MachOHandle{ioT},::Type{T}) =
    pack(h.io,T,h.bswapped ? :SwappedEndian : :NativeEndian)

function readheader(h::MachOHandle)
    seekstart(h)
    unpack(h,h.is64 ? mach_header_64 : mach_header)
end

sectionsize(sect::Union(section,section_64)) = sect.size
sectionoffset(sect::Union(section,section_64)) = sect.offset

# Access to Symbols
immutable Symbols
    h::MachOHandle
    command::symtab_command
    function Symbols(h::MachOHandle, segment::symtab_command)
        new(h,segment,start)
    end
    function Symbols(lc::LoadCmd{symtab_command})
        new(lc.h,lc.cmd)
    end
end

start(s::Symbols) = 1
done(s::Symbols,n) = n > length(s)
next(s::Symbols,n) = (s[n],n+1)

length(s::Symbols) = s.command.nsyms
function getindex(s::Symbols,n)
    if n < 1 || n > length(s)
        throw(BoundsError())
    end
    sT = s.h.is64 ? nlist_64 : nlist
    seek(s.h,s.command.symoff + (n-1)*sizeof(sT))
    unpack(s.h, sT)
end

function strtable_lookup(io::MachOHandle,command::symtab_command,offset)
    seek(io,command.stroff+offset)
    strip(readuntil(io,'\0'),'\0')
end

symname(io::MachOHandle,command::symtab_command,sym) = strtable_lookup(io, command, sym.n_strx)
symname(x::LoadCmd{symtab_command}, sym) = symname(x.h, x.cmd, sym)
segname(x::Union(segment_command_64,section_64)) = x.segname
segname(x::LoadCmd{segment_command_64}) = segname(x.cmd)
sectname(x::section_64) = x.sectname

### Fat Handle
immutable FatMachOHandle
    io::IO
    start::Uint64
    header::fat_header
    archs::Vector{fat_arch}
end

function FatMachOHandle(io,start)
    header = unpack(io,fat_header,:BigEndian)
    archs = Array(fat_arch,header.nfat_arch)
    for i in 1:header.nfat_arch
        archs[i] = unpack(io,fat_arch,:BigEndian)
    end
    FatMachOHandle(io,start,header,archs)
end

function show(io::IO,h::fat_arch)
    printentry(io,"cputype",CPUTYPES[h.cputype])
    # TODO: subtype printing
    #printentry(io,"cputype",CPUSUBTYPES[h.cputype])
    printentry(io,"offset","0x",hex(h.offset,2sizeof(Uint32)))
    printentry(io,"size","0x",hex(h.size,2sizeof(Uint32)))
    printentry(io,"align",dec(h.align))
end

function show(io::IO,h::FatMachOHandle)
    print(io, "Fat Mach Handle (",length(h.archs), " architectures)")
    for (i,arch) in enumerate(h.archs)
        println(io)
        println(io, "architecture ", i)
        print(io, arch)
    end
end

function getindex(h::FatMachOHandle,i)
    seek(h.io,h.start + h.archs[i].offset)
    readmeta(h.io,MachOHandle)
end

### Compact Unwind Support

include("compact_unwind.jl")

### DWARF support

using DWARF

function debugsections{T<:segment_commands}(seg::LoadCmd{T})
    sects = collect(Sections(seg))
    snames = map(sectname,sects)
    sections = Dict{ASCIIString,section_64}()
    for i in 1:length(snames)
        # remove leading "__"
        ind = findfirst(DWARF.DEBUG_SECTIONS,bytestring(snames[i])[3:end])
        if ind != 0
            sections[DWARF.DEBUG_SECTIONS[ind]] = sects[i]
        end
    end
    sections
end

end # module


text:
# Metida
# Copyright © 2019-2020 Vladimir Arnautov aka PharmCat <mail@pharmcat.net>

module MetidaNCA

using RecipesBase

import Base: ht_keyindex, length, length, push!
import MetidaBase
import MetidaBase: Tables, StatsBase,
PrettyTables,
AbstractIdData,
AbstractSubject,
DataSet,
AbstractSubjectResult,
AbstractResultData,
isnanormissing,
getid,
getdata,
metida_table, metida_table_, MetidaTable,
uniqueidlist,
indsdict!,
subset

export pkimport, nca!, nca, DoseTime, ElimRange, LimitRule,
setdosetime!, setkelauto!, setkelrange!, applylimitrule!,
pkplot,
getkeldata, getkelauto, getkelrange, getdosetime, subset

const LOG2 = log(2)

    include("types.jl")
    include("setkelauto.jl")
    include("setkelrange.jl")
    include("setdosetime.jl")
    include("getkeldata.jl")
    include("applylimitrule.jl")
    include("show.jl")
    include("import.jl")
    include("nca.jl")
    include("plots.jl")
    include("metidatable.jl")

end # module

ext: jl


text:
module TestSum

using Test
using Folds

using ..Utils: check_rrule

function test_rrule()
    @testset for ex in [SequentialEx(), ThreadedEx()], f in [identity, abs]
        check_rrule(Folds.sum, f, 1:10, ex)
    end
end

end  # module


text:
using Test, Flux

@testset "DeepONet" begin
    @testset "dimensions" begin
        # Test the proper construction
        # Branch net
        @test size(DeepONet((32,64,72), (24,48,72), σ, tanh).branch_net.layers[end].weight) == (72,64)
        @test size(DeepONet((32,64,72), (24,48,72), σ, tanh).branch_net.layers[end].bias) == (72,)
        # Trunk net
        @test size(DeepONet((32,64,72), (24,48,72), σ, tanh).trunk_net.layers[end].weight) == (72,48)
        @test size(DeepONet((32,64,72), (24,48,72), σ, tanh).trunk_net.layers[end].bias) == (72,)
    end

    # Accept only Int as architecture parameters
    @test_throws MethodError DeepONet((32.5,64,72), (24,48,72), σ, tanh)
    @test_throws MethodError DeepONet((32,64,72), (24.1,48,72))
end

#Just the first 16 datapoints from the Burgers' equation dataset
a = [0.83541104, 0.83479851, 0.83404712, 0.83315711, 0.83212979, 0.83096755, 0.82967374, 0.82825263, 0.82670928, 0.82504949, 0.82327962, 0.82140651, 0.81943734, 0.81737952, 0.8152405, 0.81302771]
sensors = collect(range(0, 1, length=16))'

model = DeepONet((16, 22, 30), (1, 16, 24, 30), σ, tanh; init_branch=Flux.glorot_normal, bias_trunk=false)

model(a,sensors)

#forward pass
@test size(model(a, sensors)) == (1, 16)

mgrad = Flux.Zygote.gradient((x,p)->sum(model(x,p)),a,sensors)

#gradients
@test !iszero(Flux.Zygote.gradient((x,p)->sum(model(x,p)),a,sensors)[1])
@test !iszero(Flux.Zygote.gradient((x,p)->sum(model(x,p)),a,sensors)[2])

#Output size of branch and trunk subnets should be same
branch = Chain(Dense(16, 22), Dense(22, 30))
trunk = Chain(Dense(1, 16), Dense(16, 24), Dense(24, 32))
m = DeepONet(branch, trunk)
@test_throws AssertionError DeepONet((32,64,70), (24,48,72), σ, tanh)
@test_throws DimensionMismatch m(a, sensors)


text:
using FiniteVolumeMethod
using BenchmarkTools

operators = ((gradient, grid -> CellVar(x -> sum(x), grid)), (divergence, grid -> FaceVar(x -> x, grid)))

suite = BenchmarkGroup()
refinement_steps = 5

for dim = 1
    suite[dim] = BenchmarkGroup()

    # Create a uniform grid
    uniform = TreeGrid(zeros(dim))
    for step=1:refinement_steps
        FiniteVolumeMethod.refine!(uniform, [uniform.tree], recurse=true)
    end

    # Create a locally refined grid
    refined = TreeGrid(zeros(dim))
    child = refined.tree
    for step=1:refinement_steps
        FiniteVolumeMethod.refine!(refined, [child])
        child = child.children[1+mod(step-1, 1<<dim)]
    end

    cartesian = CartesianGrid(1 << refinement_steps, dim)

    grids = [refined, uniform, cartesian]
    gridNames = ["refined", "uniform", "cartesian"]

    for (oper, constructor) ∈ operators
        suite[dim][string(oper)] = BenchmarkGroup()

        for (i, grid) ∈ enumerate(grids)
            var = constructor(grid)
            suite[dim][string(oper)][gridNames[i]] = @benchmarkable $(oper)($var, $grid)
        end
    end
end


text:
using Revise
using DifferenceEquations
using Distributions
using Optim
using Random

function ar2_transition(u, p, t) # f
    A = [p[1] p[2]; 1 0]
    return A * u
end

function ar2_noise(u, p, t, noise::Nothing) # g
    z = [p[3] * randn(); 0]
    println(z)
    return [p[3] * randn(); 0]
end

function ar2_noise(u, p, t, noise) # g
    return [noise; 0]
end

function ar2_observation(u, p, t, noise) # h
    # @info "Observation" u
    return [1 0] * u
end

Random.seed!(1)

p = [0.5, -0.25]
tspan = (1, 15)

phi1 = 0.5
phi2 = -0.25
sigma = 5.0

Y = zeros(tspan[2])
n = sigma .* randn(length(Y))
Y += n

for t in 3:length(Y)
    Y[t] += phi1 * Y[t-1] + phi2 * Y[t-2]
end

u0 = [Y[2]; Y[1]]

# prob = StateSpaceProblem(
#     ar2_transition,
#     ar2_noise,
#     ar2_observation,
#     MvNormal(1,1),
#     [0.0, 0.0],
#     tspan,
#     nothing,
#     Y
# )

prob = StateSpaceProblem(
    ar2_transition,
    ar2_noise,
    ar2_observation,
    MvNormal(1,1),
    [0.0, 0.0],
    tspan,
    nothing, #n,
    Y
)

sol = DifferenceEquations._solve(prob, [phi1, phi2, sigma])
# target(x) = -DifferenceEquations._solve(prob, x).likelihood
# res = optimize(target, [0.0, 0.0, 1.0], LBFGS())


text:
using Pkg

Pkg.activate("FibonacciOptimization")
Pkg.test()


text:
"""
    fathom!(y::BnBModel)

Fathoms the stack by value dominance.
"""
function fathom!(y::BnBModel)
  loc = find(z->(z>y.UBDg),y.LBD)
  deleteat!(y.LBD,loc)
  deleteat!(y.UBD,loc)
  deleteat!(y.box,loc)
  deleteat!(y.id,loc)
end


text:
struct History{T<:Real}
    # model setings
    model::Symbol
    params::NamedTuple

    # kernel settings
    kernel::Symbol
    γ::T
    scale::Bool
    mlkernels::Bool
    precomputed::Bool

    # solver settings
    maxiter::Int
    at::Int
    seed::Int
    pupdate::T
    ε::T

    # history
    iter::Vector{Int}
    k::Vector{Int}
    l::Vector{Int}
    Δ::Vector{T}
    Lprimal::Vector{T}
    Ldual::Vector{T}
    gap::Vector{T}
    train::Dict{Int, NamedTuple}
    solution::Dict{Symbol, NamedTuple}
end

function Base.show(io::IO, h::History)
    print(io, "History: $(h.model)", values(h.params), " with ", h.kernel, " kernel")
end

function Base.Dict(h::History)
    d = Dict(key => getproperty(h, key) for key in fieldnames(History))
    for key in keys(d[:params])
        d[key] = getproperty(d[:params], key)
    end
    return d
end

function History(
    model::M,
    ker::KernelType{T, D},
    seed;
    maxiter = 20000,
    at = max(1, round(Int, maxiter/100)),
    pupdate = 0.9,
    ε::Real = 1e-4,
    kwargs...
) where {T<:Real, M<:Model, D<:Kernel}

    bar = Progress(maxiter, 1, "⋅ Training: ")

    return bar, History{T}(
        M.name.name,
        parameters(model),
        D.name.name,
        ker.γ,
        ker.scale,
        ker.mlkernels,
        ker.precomputed,
        maxiter,
        at,
        seed,
        pupdate,
        ε,
        Int[],
        Int[],
        Int[],
        Int[],
        T[],
        T[],
        T[],
        Dict{Int, NamedTuple}(),
        Dict{Symbol, NamedTuple}(),
    )
end

function log!(h::History, model, K, iter, k, l, Δ)
    append!(h.iter, iter)
    append!(h.k, k)
    append!(h.l, l)
    append!(h.Δ, Δ)
    if Δ != 0 || iter == 0
        Lprimal, Ldual, gap = objective(model, K)
    else
        Lprimal, Ldual, gap = h.Lprimal[end], h.Ldual[end], h.gap[end]
    end
    append!(h.Lprimal, Lprimal)
    append!(h.Ldual, Ldual)
    append!(h.gap, gap)
    if mod(iter, h.at) == 0
        h.train[iter] = (; s = extract_scores(model, K), extract_state(model, K)...)
    end
    return [
        (:L_primal, Lprimal),
        (:L_dual, Ldual),
        (:gap, gap),
        (:iter, iter),
        (:stop, gap/h.gap[1]),
    ]
end

function terminate(h::History)
    g0, g = h.gap[1], h.gap[end]
    return isnan(g) || isinf(g) || g/g0 <= h.ε
end

function solve!(
    model,
    X::Matrix{T},
    y,
    ker::KernelType{T, D};
    seed = 1234,
    kwargs...
) where {T<:Real, D <: Kernel}

    # initialization
    printstyled("$(model): \n"; color=:magenta, bold=true)
    printstyled("⋅ Initialization: "; color=:green)
    tm = @timed begin
        Random.seed!(seed)
        kernel = init(ker, size(X, 2))
        K = KernelMatrix(model, X, y, kernel; precomputed = ker.precomputed)
        initialization!(model, K)
    end
    printstyled(ProgressMeter.durationstring(tm.time), " \n"; color=:green)

    # progress bar and history
    bar, history = History(model, ker, seed; kwargs...)
    k, l, Δ = 0, 0, zero(T)
    vals = log!(history, model, K, 0, k, l, Δ)

    # train
    for iter in 1:history.maxiter
        k = if Δ == 0 || rand() > history.pupdate
            rand(1:K.n)
        else
            l
        end
        Δ, k, l = update!(model, K; k)

        # update progress bar
        vals = log!(history, model, K, iter, k, l, Δ)
        next!(bar; showvalues = vals)

        # stop condition
        terminate(history) && break
    end
    finish!(bar; showvalues = vals)

    # evaluation
    printstyled("⋅ Evaluation: "; color=:green)
    tm = @timed begin
        history.solution[:train] = (;y, s = extract_scores(model, K), extract_state(model, K)...)
    end
    printstyled(ProgressMeter.durationstring(tm.time), " \n"; color=:green)
    return history
end

function predict(
    model,
    X::Matrix{T},
    y,
    ker::KernelType,
    Xtest;
    chunksize = 100,
) where {T<:Real}

    y = BitVector(y)
    αβ = copy(model.state.αβ)
    nα, ~, perm = permutation(model, y)
    αβ[(nα + 1):end] .*= -1
    kernel = init(ker, size(X, 2))

    s = zeros(T, size(Xtest, 1))
    for rows in partition(1:size(Xtest, 1), chunksize)
        K = kernelmatrix(kernel, X[perm, :], Xtest[rows, :])
        s[rows] .= vec(αβ' * K)
    end
    return s
end

# SVM
function solve!(
    model::SVM,
    X::Matrix{T},
    y,
    ker::KernelType{T, D};
    seed = 1234,
    kwargs...
) where {T<:Real, D <: Kernel}

    # initialization
    printstyled("$(model): \n"; color=:magenta, bold=true)
    printstyled("⋅ Initialization: "; color=:green)
    tm = @timed begin
        Random.seed!(seed)
        if ker.precomputed
            kernel = init(ker, size(X, 2))
            K = kernelmatrix(kernel, X)
            kernel_type = LIBSVM.Kernel.Precomputed
        else
            K = X'
            kernel_type = get_kernel(D)
        end
    end
    printstyled(ProgressMeter.durationstring(tm.time), " \n"; color=:green)

    # train
    printstyled("⋅ Training: "; color=:green)
    tm = @timed begin
        ~, history = History(model, ker, seed; kwargs...)
        svm_model = LIBSVM.svmtrain(
            K,
            Int.(y);
            probability = true,
            cost = Float64(model.C),
            gamma = Float64(compute_gamma(ker, size(X, 2))),
            kernel = kernel_type,
            epsilon = Float64(history.ε),
        )
        model.state = svm_model
    end
    printstyled(ProgressMeter.durationstring(tm.time), " \n"; color=:green)

    # evaluation
    printstyled("⋅ Evaluation: "; color=:green)
    tm = @timed begin
        history.solution[:train] = (;y, s = extract_scores(model, K))
    end
    printstyled(ProgressMeter.durationstring(tm.time), " \n"; color=:green)
    return history
end

function predict(
    model::SVM,
    X::Matrix{T},
    y,
    ker::KernelType,
    Xtest;
    chunksize = 100,
) where {T<:Real}

    if ker.precomputed
        kernel = init(ker, size(X, 2))

        s = zeros(T, size(Xtest, 1))
        for rows in partition(1:size(Xtest, 1), chunksize)
            K = kernelmatrix(kernel, X, Xtest[rows, :])
            s[rows] .= extract_scores(model, K)
        end
    else
        s = extract_scores(model, Xtest')
    end
    return s
end


text:
n = rand(1:10)
@test matrixdepot("lehmer", Float64, n) == matrixdepot("lehmer", n)

A = ones(n, 1) * [1:n;]'
A = A./A'
A = tril(A) + tril(A)' - eye(n)

@test_approx_eq matrixdepot("lehmer", n) A
println("'lehmer' passed test...")


text:
"""
    get_momenta(
        rx  :: NTuple{2, Float64},
        ry  :: NTuple{2, Float64},
        rz  :: NTuple{2, Float64},
        num :: NTuple{3, Int64}
        )   :: Matrix{Float64}

Generate a uniform momentum space discretization within a cuboid. rx, ry and rz are the respective cartesian boundaries.
num = (num_x, num_y, num_z) contains the desired number of steps (spacing h = (r[2] - r[1]) / num) along the respective axis.
Returns a matrix k, with k[:, n] being the n-th momentum vector.
"""
function get_momenta(
    rx  :: NTuple{2, Float64},
    ry  :: NTuple{2, Float64},
    rz  :: NTuple{2, Float64},
    num :: NTuple{3, Int64}
    )   :: Matrix{Float64}

    # allocate mesh
    momenta = zeros(Float64, 3, (num[1] + 1) * (num[2] + 1) * (num[3] + 1))

    # fill mesh
    for nx in 0 : num[1]
        for ny in 0 : num[2]
            for nz in 0 : num[3]
                # compute linear index
                idx = nz + 1 + (num[3] + 1) * ny + (num[3] + 1) * (num[2] + 1) * nx

                # compute kx
                if num[1] > 0
                    momenta[1, idx] = rx[1] + nx * (rx[2] - rx[1]) / num[1]
                end

                # compute ky
                if num[2] > 0
                    momenta[2, idx] = ry[1] + ny * (ry[2] - ry[1]) / num[2]
                end

                # compute kz
                if num[3] > 0
                    momenta[3, idx] = rz[1] + nz * (rz[2] - rz[1]) / num[3]
                end
            end
        end
    end

    return momenta
end

"""
    get_path(
        nodes :: Vector{Vector{Float64}},
        nums  :: Vector{Int64}
        )     :: Tuple{Vector{Float64}, Matrix{Float64}}

Generate a discrete path in momentum space linearly connecting the given nodes (passed via their cartesian coordinates [kx, ky, kz]).
nums[i] is the desired number of points between node[i] and node[i + 1], including node[i] and excluding node[i + 1].
Returns a tuple (l, k) where k[:, n] is the n-th momentum vector and l[n] is the distance to node[1] along the generated path.
"""
function get_path(
    nodes :: Vector{Vector{Float64}},
    nums  :: Vector{Int64}
    )     :: Tuple{Vector{Float64}, Matrix{Float64}}

    # sanity check: need some number of points for each path increment between two nodes
    @assert length(nums) == length(nodes) - 1 "For N nodes 'nums' must be of length N - 1."

    # allocate output arrays
    num  = sum(nums) + 1
    dist = zeros(num)
    path = zeros(3, num)

    # iterate over path increments between nodes
    idx = 0

    for i in 1 : length(nums)
        dif   = nodes[i + 1] .- nodes[i]
        step  = dif ./ nums[i]
        width = norm(step)

        # fill path from node i to node i + 1 (excluding node i + 1)
        for j in 1 : nums[i]
            dist[idx + j + 1]  = dist[idx + 1] + width * j
            path[:, idx + j]  .= nodes[i] .+ step .* (j - 1)
        end

        idx += nums[i]
    end

    # set last point in the path to be the last node
    path[:, end] = nodes[end]

    return dist, path
end

"""
    compute_structure_factor(
        χ :: Vector{Float64},
        k :: Matrix{Float64},
        l :: Lattice,
        r :: Reduced_lattice
        ) :: Vector{Float64}

Compute the structure factor for given real space correlations χ on irreducible lattice sites.
k[:, n] is the n-th discrete momentum vector.
Return structure factor s, where s[n] is the value for the n-th momentum.
"""
function compute_structure_factor(
    χ :: Vector{Float64},
    k :: Matrix{Float64},
    l :: Lattice,
    r :: Reduced_lattice
    ) :: Vector{Float64}

    # allocate structure factor
    s = zeros(Float64, size(k, 2))

    # compute all contributions from reference sites in the unitcell
    for b in eachindex(l.uc.basis)
        vec = l.uc.basis[b]
        int = get_site(vec, l)

        # compute structure factor for all momenta
        Threads.@threads for n in eachindex(s)
            @inbounds @fastmath for i in eachindex(l.sites)
                # consider only those sites in range of reference site
                index = r.project[int, i]

                if index > 0
                    val   = k[1, n] * (vec[1] - l.sites[i].vec[1])
                    val  += k[2, n] * (vec[2] - l.sites[i].vec[2])
                    val  += k[3, n] * (vec[3] - l.sites[i].vec[3])
                    s[n] += cos(val) * χ[index]
                end
            end
        end
    end

    s ./= length(l.uc.basis)

    return s
end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    ExecutableProduct(prefix, "hybrid-Lambda", :hybridlambda),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/cecileane/HybridLambdaBuilder/releases/download/v0.6.2-beta"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc) => ("$bin_prefix/HybridLambda.v0.6.2-beta.aarch64-linux-gnu.tar.gz", "cb9723304dfc7a4d67ba305d6a7c3041c665c8db53b170fde320d9533f7cb367"),
    Linux(:aarch64, libc=:musl) => ("$bin_prefix/HybridLambda.v0.6.2-beta.aarch64-linux-musl.tar.gz", "6291d48b2e6b569dbc2bc011f5d768eb44e5611aa1943faae720cb1ce281b5d3"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf) => ("$bin_prefix/HybridLambda.v0.6.2-beta.arm-linux-gnueabihf.tar.gz", "807aafa70fb1d022b0b796abc9e2a89a030ef79357cc3880e1a8dcae3823a667"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf) => ("$bin_prefix/HybridLambda.v0.6.2-beta.arm-linux-musleabihf.tar.gz", "381f40c0d6892bd7ceafdf2cce0d35b06617388ff9d090dcc6a7fa93bc2bc006"),
    Linux(:i686, libc=:glibc) => ("$bin_prefix/HybridLambda.v0.6.2-beta.i686-linux-gnu.tar.gz", "c49e49736aa7cc5b82f9d40fdffd7effba2f2aaaf77a8090dee5c3c064c0895d"),
    Linux(:powerpc64le, libc=:glibc) => ("$bin_prefix/HybridLambda.v0.6.2-beta.powerpc64le-linux-gnu.tar.gz", "8643c0708d1cce0080984e08168a68f1361a2ae1876cf95da8ab39a0cd1a2423"),
    MacOS(:x86_64) => ("$bin_prefix/HybridLambda.v0.6.2-beta.x86_64-apple-darwin14.tar.gz", "4ed16d6775ddcef45e63c1a769d9bdf731fc9eeeb5d44572a486548f53c27f71"),
    Linux(:x86_64, libc=:glibc) => ("$bin_prefix/HybridLambda.v0.6.2-beta.x86_64-linux-gnu.tar.gz", "cf27dbf2ced3b6d167112a8ca6c7c2aa818cc622494bfc5bca8ede03c5da6ea8"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    @warn "Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!"
end

if dl_info !== nothing
    # If we have a download, and we are unsatisfied (or the version we're
    # trying to install is not itself installed) then load it up!
    if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
        # Download and install binaries
        install(dl_info...; prefix=prefix, force=true, verbose=verbose)
    end
end

if dl_info !== nothing || !unsatisfied
    # Write out deps.jl that will contain mappings for our products:
    #   import Libdl
    #   const hybridlambda = joinpath(dirname(@__FILE__), "usr/bin/hybrid-Lambda")
    #   definition of function check_deps()
    write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)
end


text:
# using GSL.sf_legendre_sphPlm

# function sphharm(l::Integer, m::Integer, theta::Real, phi::Real)
# 	@assert l >= 0 && abs(m) <= l
# 	if m >= 0
# 		return exp(m*phi*im)*sf_legendre_sphPlm(l,m,cos(theta))
# 	else
#     	return (-1.)^m*exp(m*phi*im)*sf_legendre_sphPlm(l,-m,cos(theta))
#     end
# end

doc"""
`sphharm(l,m,θ,ϕ)`

computes the spherical harmonic function (standard definition used in quantum mechanics)

$ Y_l^m(\theta,\phi) = (-1)^m \sqrt{{2l+1 \over 4\pi} {(l-m)!\over (l+m)!}} P_l^m(cos \theta) e^{im\phi} $

which are normalized such that

$ \int Y_l^m(\theta,\phi)^* Y_{l'}^{m'}(\theta,\phi) d\Omega = \delta_{ll'} \delta_{mm'} $
"""
# sphharm(l::Int, m::Int, θ::Real, ϕ::Real) = sphharmtab[l,m](θ,ϕ)::Complex{Float64}
sphharm(l::Int, m::Int, θ::Real, ϕ::Real) = sphharmtab[i_lm(l,m)](θ,ϕ)::Complex{Float64}

doc"""
`realsphharm(l,m,θ,ϕ)`

computes the real spherical harmonic function

$ Y_{lm} = {i \over \sqrt{2}} \left(Y_l^m - (-1)^m Y_l^{-m}\right), \quad m<0 $

$ Y_{lm} = Y_l^m, \quad m=0 $

$ Y_{lm} = {1 \over \sqrt{2}} \left(Y_l^{-m} + (-1)^m Y_l^m\right), \quad m>0 $

"""
# realsphharm(l::Int, m::Int, θ::Real, ϕ::Real) = realsphharmtab[l,m](θ,ϕ)::Float64
realsphharm(l::Int, m::Int, θ::Real, ϕ::Real) = realsphharmtab[i_lm(l,m)](θ,ϕ)::Float64

sphharm00(θ::Float64, ϕ::Float64) = complex(1/(2*sqrt(pi)))
sphharm1m1(θ::Float64, ϕ::Float64) = ((1/2)*sqrt(3/(2*pi))*sin(θ))/exp(ϕ*im)
sphharm10(θ::Float64, ϕ::Float64) = complex((1/2)*sqrt(3/pi)*cos(θ))
sphharm11(θ::Float64, ϕ::Float64) = (-(1/2))*exp(ϕ*im)*sqrt(3/(2*pi))*sin(θ)
sphharm2m2(θ::Float64, ϕ::Float64) = ((1/4)*sqrt(15/(2*pi))*sin(θ)^2)/exp(2*ϕ*im)
sphharm2m1(θ::Float64, ϕ::Float64) = ((1/2)*sqrt(15/(2*pi))*cos(θ)*sin(θ))/exp(ϕ*im)
sphharm20(θ::Float64, ϕ::Float64) = complex((1/8)*sqrt(5/pi)*(1 + 3*cos(2*θ)))
sphharm21(θ::Float64, ϕ::Float64) = (-(1/2))*exp(ϕ*im)*sqrt(15/(2*pi))*cos(θ)*sin(θ)
sphharm22(θ::Float64, ϕ::Float64) = (1/4)*exp(2*ϕ*im)*sqrt(15/(2*pi))*sin(θ)^2
sphharm3m3(θ::Float64, ϕ::Float64) = ((1/8)*sqrt(35/pi)*sin(θ)^3)/exp(3*ϕ*im)
sphharm3m2(θ::Float64, ϕ::Float64) = ((1/4)*sqrt(105/(2*pi))*cos(θ)*sin(θ)^2)/exp(2*ϕ*im)
sphharm3m1(θ::Float64, ϕ::Float64) = ((1/16)*sqrt(21/pi)*(3 + 5*cos(2*θ))*sin(θ))/exp(ϕ*im)
sphharm30(θ::Float64, ϕ::Float64) = complex((1/16)*sqrt(7/pi)*(3*cos(θ) + 5*cos(3*θ)))
sphharm31(θ::Float64, ϕ::Float64) = (-(1/16))*exp(ϕ*im)*sqrt(21/pi)*(3 + 5*cos(2*θ))*sin(θ)
sphharm32(θ::Float64, ϕ::Float64) = (1/4)*exp(2*ϕ*im)*sqrt(105/(2*pi))*cos(θ)*sin(θ)^2
sphharm33(θ::Float64, ϕ::Float64) = (-(1/8))*exp(3*ϕ*im)*sqrt(35/pi)*sin(θ)^3
sphharm4m4(θ::Float64, ϕ::Float64) = ((3/16)*sqrt(35/(2*pi))*sin(θ)^4)/exp(4*ϕ*im)
sphharm4m3(θ::Float64, ϕ::Float64) = ((3/8)*sqrt(35/pi)*cos(θ)*sin(θ)^3)/exp(3*ϕ*im)
sphharm4m2(θ::Float64, ϕ::Float64) = ((3/16)*sqrt(5/(2*pi))*(5 + 7*cos(2*θ))*sin(θ)^2)/exp(2*ϕ*im)
sphharm4m1(θ::Float64, ϕ::Float64) = ((3/32)*sqrt(5/pi)*(1 + 7*cos(2*θ))*sin(2*θ))/exp(ϕ*im)
sphharm40(θ::Float64, ϕ::Float64) = complex((3*(9 + 20*cos(2*θ) + 35*cos(4*θ)))/(128*sqrt(pi)))
sphharm41(θ::Float64, ϕ::Float64) = (-(3/32))*exp(ϕ*im)*sqrt(5/pi)*(1 + 7*cos(2*θ))*sin(2*θ)
sphharm42(θ::Float64, ϕ::Float64) = (3/16)*exp(2*ϕ*im)*sqrt(5/(2*pi))*(5 + 7*cos(2*θ))*sin(θ)^2
sphharm43(θ::Float64, ϕ::Float64) = (-(3/8))*exp(3*ϕ*im)*sqrt(35/pi)*cos(θ)*sin(θ)^3
sphharm44(θ::Float64, ϕ::Float64) = (3/16)*exp(4*ϕ*im)*sqrt(35/(2*pi))*sin(θ)^4
sphharm5m5(θ::Float64, ϕ::Float64) = ((3/32)*sqrt(77/pi)*sin(θ)^5)/exp(5*ϕ*im)
sphharm5m4(θ::Float64, ϕ::Float64) = ((3/16)*sqrt(385/(2*pi))*cos(θ)*sin(θ)^4)/exp(4*ϕ*im)
sphharm5m3(θ::Float64, ϕ::Float64) = ((1/64)*sqrt(385/pi)*(7 + 9*cos(2*θ))*sin(θ)^3)/exp(3*ϕ*im)
sphharm5m2(θ::Float64, ϕ::Float64) = ((1/16)*sqrt(1155/(2*pi))*cos(θ)*(1 + 3*cos(2*θ))*sin(θ)^2)/exp(2*ϕ*im)
sphharm5m1(θ::Float64, ϕ::Float64) = ((1/16)*sqrt(165/(2*pi))*(1 - 14*cos(θ)^2 +21*cos(θ)^4)*sin(θ))/exp(ϕ*im)
sphharm50(θ::Float64, ϕ::Float64) = complex((1/256)*sqrt(11/pi)*(30*cos(θ) + 35*cos(3*θ) +63*cos(5*θ)))
sphharm51(θ::Float64, ϕ::Float64) = (-(1/16))*exp(ϕ*im)*sqrt(165/(2*pi))*(1 - 14*cos(θ)^2 + 21*cos(θ)^4)*sin(θ)
sphharm52(θ::Float64, ϕ::Float64) = (1/16)*exp(2*ϕ*im)*sqrt(1155/(2*pi))*cos(θ)*(1 + 3*cos(2*θ))*sin(θ)^2
sphharm53(θ::Float64, ϕ::Float64) = (-(1/64))*exp(3*ϕ*im)*sqrt(385/pi)*(7 + 9*cos(2*θ))*sin(θ)^3
sphharm54(θ::Float64, ϕ::Float64) = (3/16)*exp(4*ϕ*im)*sqrt(385/(2*pi))*cos(θ)*sin(θ)^4
sphharm55(θ::Float64, ϕ::Float64) = (-(3/32))*exp(5*ϕ*im)*sqrt(77/pi)*sin(θ)^5
sphharm6m6(θ::Float64, ϕ::Float64) = ((1/64)*sqrt(3003/pi)*sin(θ)^6)/exp(6*ϕ*im)
sphharm6m5(θ::Float64, ϕ::Float64) = ((3/32)*sqrt(1001/pi)*cos(θ)*sin(θ)^5)/exp(5*ϕ*im)
sphharm6m4(θ::Float64, ϕ::Float64) = ((3/64)*sqrt(91/(2*pi))*(9 + 11*cos(2*θ))*sin(θ)^4)/exp(4*ϕ*im)
sphharm6m3(θ::Float64, ϕ::Float64) = ((1/64)*sqrt(1365/pi)*cos(θ)*(5 + 11*cos(2*θ))*sin(θ)^3)/exp(3*ϕ*im)
sphharm6m2(θ::Float64, ϕ::Float64) = ((1/64)*sqrt(1365/pi)*(1 - 18*cos(θ)^2 +33*cos(θ)^4)*sin(θ)^2)/exp(2*ϕ*im)
sphharm6m1(θ::Float64, ϕ::Float64) = ((1/16)*sqrt(273/(2*pi))*cos(θ)*(5 - 30*cos(θ)^2 +33*cos(θ)^4)*sin(θ))/exp(ϕ*im)
sphharm60(θ::Float64, ϕ::Float64) = complex((1/32)*sqrt(13/pi)*(-5 + 105*cos(θ)^2 -315*cos(θ)^4 + 231*cos(θ)^6))
sphharm61(θ::Float64, ϕ::Float64) = (-(1/16))*exp(ϕ*im)*sqrt(273/(2*pi))*cos(θ)*(5 - 30*cos(θ)^2 + 33*cos(θ)^4)*sin(θ)
sphharm62(θ::Float64, ϕ::Float64) = (1/64)*exp(2*ϕ*im)*sqrt(1365/pi)*(1 - 18*cos(θ)^2 +33*cos(θ)^4)*sin(θ)^2
sphharm63(θ::Float64, ϕ::Float64) = (-(1/64))*exp(3*ϕ*im)*sqrt(1365/pi)*cos(θ)*(5 + 11*cos(2*θ))*sin(θ)^3
sphharm64(θ::Float64, ϕ::Float64) = (3/64)*exp(4*ϕ*im)*sqrt(91/(2*pi))*(9 + 11*cos(2*θ))*sin(θ)^4
sphharm65(θ::Float64, ϕ::Float64) = (-(3/32))*exp(5*ϕ*im)*sqrt(1001/pi)*cos(θ)*sin(θ)^5
sphharm66(θ::Float64, ϕ::Float64) = (1/64)*exp(6*ϕ*im)*sqrt(3003/pi)*sin(θ)^6
sphharm7m7(θ::Float64, ϕ::Float64) = ((3/64)*sqrt(715/(2*pi))*sin(θ)^7)/exp(7*ϕ*im)
sphharm7m6(θ::Float64, ϕ::Float64) = ((3/64)*sqrt(5005/pi)*cos(θ)*sin(θ)^6)/exp(6*ϕ*im)
sphharm7m5(θ::Float64, ϕ::Float64) = ((3/64)*sqrt(385/(2*pi))*(-1 + 13*cos(θ)^2)*sin(θ)^5)/exp(5*ϕ*im)
sphharm7m4(θ::Float64, ϕ::Float64) = ((3/64)*sqrt(385/(2*pi))*cos(θ)*(7 + 13*cos(2*θ))*sin(θ)^4)/exp(4*ϕ*im)
sphharm7m3(θ::Float64, ϕ::Float64) = ((3/64)*sqrt(35/(2*pi))*(3 - 66*cos(θ)^2 +143*cos(θ)^4)*sin(θ)^3)/exp(3*ϕ*im)
sphharm7m2(θ::Float64, ϕ::Float64) = ((3/64)*sqrt(35/pi)*cos(θ)*(15 - 110*cos(θ)^2 +143*cos(θ)^4)*sin(θ)^2)/exp(2*ϕ*im)
sphharm7m1(θ::Float64, ϕ::Float64) = ((1/64)*sqrt(105/(2*pi))*(-5 + 135*cos(θ)^2 -495*cos(θ)^4 + 429*cos(θ)^6)*sin(θ))/exp(ϕ*im)
sphharm70(θ::Float64, ϕ::Float64) = complex((1/32)*sqrt(15/pi)*cos(θ)*(-35 + 315*cos(θ)^2 -693*cos(θ)^4 + 429*cos(θ)^6))
sphharm71(θ::Float64, ϕ::Float64) = (-(1/64))*exp(ϕ*im)*sqrt(105/(2*pi))*(-5 + 135*cos(θ)^2 - 495*cos(θ)^4 + 429*cos(θ)^6)*sin(θ)
sphharm72(θ::Float64, ϕ::Float64) = (3/64)*exp(2*ϕ*im)*sqrt(35/pi)*cos(θ)*(15 - 110*cos(θ)^2 + 143*cos(θ)^4)*sin(θ)^2
sphharm73(θ::Float64, ϕ::Float64) = (-(3/64))*exp(3*ϕ*im)*sqrt(35/(2*pi))*(3 - 66*cos(θ)^2 + 143*cos(θ)^4)*sin(θ)^3
sphharm74(θ::Float64, ϕ::Float64) = (3/64)*exp(4*ϕ*im)*sqrt(385/(2*pi))*cos(θ)*(7 + 13*cos(2*θ))*sin(θ)^4
sphharm75(θ::Float64, ϕ::Float64) = (-(3/64))*exp(5*ϕ*im)*sqrt(385/(2*pi))*(-1 + 13*cos(θ)^2)*sin(θ)^5
sphharm76(θ::Float64, ϕ::Float64) = (3/64)*exp(6*ϕ*im)*sqrt(5005/pi)*cos(θ)*sin(θ)^6
sphharm77(θ::Float64, ϕ::Float64) = (-(3/64))*exp(7*ϕ*im)*sqrt(715/(2*pi))*sin(θ)^7
sphharm8m8(θ::Float64, ϕ::Float64) = ((3/256)*sqrt(12155/(2*pi))*sin(θ)^8)/exp(8*ϕ*im)
sphharm8m7(θ::Float64, ϕ::Float64) = ((3/64)*sqrt(12155/(2*pi))*cos(θ)*sin(θ)^7)/exp(7*ϕ*im)
sphharm8m6(θ::Float64, ϕ::Float64) = ((1/256)*sqrt(7293/pi)*(13 + 15*cos(2*θ))*sin(θ)^6)/exp(6*ϕ*im)
sphharm8m5(θ::Float64, ϕ::Float64) = ((3/128)*sqrt(17017/(2*pi))*cos(θ)*(3 + 5*cos(2*θ))*sin(θ)^5)/exp(5*ϕ*im)
sphharm8m4(θ::Float64, ϕ::Float64) = ((3/128)*sqrt(1309/(2*pi))*(1 - 26*cos(θ)^2 +65*cos(θ)^4)*sin(θ)^4)/exp(4*ϕ*im)
sphharm8m3(θ::Float64, ϕ::Float64) = ((1/64)*sqrt(19635/(2*pi))*cos(θ)*(3 - 26*cos(θ)^2 + 39*cos(θ)^4)*sin(θ)^3)/exp(3*ϕ*im)
sphharm8m2(θ::Float64, ϕ::Float64) = ((3/128)*sqrt(595/pi)*(-1 + 33*cos(θ)^2 - 143*cos(θ)^4 + 143*cos(θ)^6)*sin(θ)^2)/exp(2*ϕ*im)
sphharm8m1(θ::Float64, ϕ::Float64) = (1/4096)*((3*sqrt(17/(2*pi))*(178 + 869*cos(2*θ) +286*cos(4*θ) + 715*cos(6*θ))*sin(2*θ))/exp(ϕ*im))
sphharm80(θ::Float64, ϕ::Float64) = complex((1/256)*sqrt(17/pi)*(35 - 1260*cos(θ)^2 + 6930*cos(θ)^4 -12012*cos(θ)^6 + 6435*cos(θ)^8))
sphharm81(θ::Float64, ϕ::Float64) = -((1/4096)*(3*exp(ϕ*im)*sqrt(17/(2*pi))*(178 + 869*cos(2*θ) + 286*cos(4*θ) +715*cos(6*θ))*sin(2*θ)))
sphharm82(θ::Float64, ϕ::Float64) = (3/128)*exp(2*ϕ*im)*sqrt(595/pi)*(-1 + 33*cos(θ)^2 -143*cos(θ)^4 + 143*cos(θ)^6)*sin(θ)^2
sphharm83(θ::Float64, ϕ::Float64) = (-(1/64))*exp(3*ϕ*im)*sqrt(19635/(2*pi))*cos(θ)*(3 - 26*cos(θ)^2 + 39*cos(θ)^4)*sin(θ)^3
sphharm84(θ::Float64, ϕ::Float64) = (3/128)*exp(4*ϕ*im)*sqrt(1309/(2*pi))*(1 - 26*cos(θ)^2 + 65*cos(θ)^4)*sin(θ)^4
sphharm85(θ::Float64, ϕ::Float64) = (-(3/128))*exp(5*ϕ*im)*sqrt(17017/(2*pi))*cos(θ)*(3 + 5*cos(2*θ))*sin(θ)^5
sphharm86(θ::Float64, ϕ::Float64) = (1/256)*exp(6*ϕ*im)*sqrt(7293/pi)*(13 + 15*cos(2*θ))*sin(θ)^6
sphharm87(θ::Float64, ϕ::Float64) = (-(3/64))*exp(7*ϕ*im)*sqrt(12155/(2*pi))*cos(θ)*sin(θ)^7
sphharm88(θ::Float64, ϕ::Float64) = (3/256)*exp(8*ϕ*im)*sqrt(12155/(2*pi))*sin(θ)^8

realsphharm00(θ::Float64, ϕ::Float64) = 1/(2*sqrt(pi))
realsphharm1m1(θ::Float64, ϕ::Float64) = (1/2)*sqrt(3/pi)*sin(ϕ)*sin(θ)
realsphharm10(θ::Float64, ϕ::Float64) = (1/2)*sqrt(3/pi)*cos(θ)
realsphharm11(θ::Float64, ϕ::Float64) = (1/2)*sqrt(3/pi)*cos(ϕ)*sin(θ)
realsphharm2m2(θ::Float64, ϕ::Float64) = (-(1/4))*sqrt(15/pi)*sin(2*ϕ)*sin(θ)^2
realsphharm2m1(θ::Float64, ϕ::Float64) = (1/2)*sqrt(15/pi)*cos(θ)*sin(ϕ)*sin(θ)
realsphharm20(θ::Float64, ϕ::Float64) = (1/8)*sqrt(5/pi)*(1 + 3*cos(2*θ))
realsphharm21(θ::Float64, ϕ::Float64) = (1/2)*sqrt(15/pi)*cos(ϕ)*cos(θ)*sin(θ)
realsphharm22(θ::Float64, ϕ::Float64) = (1/4)*sqrt(15/pi)*cos(2*ϕ)*sin(θ)^2
realsphharm3m3(θ::Float64, ϕ::Float64) = (1/4)*sqrt(35/(2*pi))*sin(3*ϕ)*sin(θ)^3
realsphharm3m2(θ::Float64, ϕ::Float64) = (-(1/4))*sqrt(105/pi)*cos(θ)*sin(2*ϕ)*sin(θ)^2
realsphharm3m1(θ::Float64, ϕ::Float64) = (1/16)*sqrt(21/(2*pi))*sin(ϕ)*(sin(θ) + 5*sin(3*θ))
realsphharm30(θ::Float64, ϕ::Float64) = (1/8)*sqrt(7/pi)*cos(θ)*(-1 + 5*cos(2*θ))
realsphharm31(θ::Float64, ϕ::Float64) = (1/16)*sqrt(21/(2*pi))*cos(ϕ)*(sin(θ) + 5*sin(3*θ))
realsphharm32(θ::Float64, ϕ::Float64) = (1/4)*sqrt(105/pi)*cos(2*ϕ)*cos(θ)*sin(θ)^2
realsphharm33(θ::Float64, ϕ::Float64) = (1/4)*sqrt(35/(2*pi))*cos(3*ϕ)*sin(θ)^3
realsphharm4m4(θ::Float64, ϕ::Float64) = (-(3/16))*sqrt(35/pi)*sin(4*ϕ)*sin(θ)^4
realsphharm4m3(θ::Float64, ϕ::Float64) = (3/4)*sqrt(35/(2*pi))*cos(θ)*sin(3*ϕ)*sin(θ)^3
realsphharm4m2(θ::Float64, ϕ::Float64) = (3/8)*sqrt(5/pi)*(1 - 7*cos(θ)^2)*sin(2*ϕ)*sin(θ)^2
realsphharm4m1(θ::Float64, ϕ::Float64) = (3/32)*sqrt(5/(2*pi))*sin(ϕ)*(2*sin(2*θ) + 7*sin(4*θ))
realsphharm40(θ::Float64, ϕ::Float64) = (3*(9 + 20*cos(2*θ) + 35*cos(4*θ)))/(128*sqrt(pi))
realsphharm41(θ::Float64, ϕ::Float64) = (3/32)*sqrt(5/(2*pi))*cos(ϕ)*(2*sin(2*θ) + 7*sin(4*θ))
realsphharm42(θ::Float64, ϕ::Float64) = (3/16)*sqrt(5/pi)*cos(2*ϕ)*(5 + 7*cos(2*θ))*sin(θ)^2
realsphharm43(θ::Float64, ϕ::Float64) = (3/4)*sqrt(35/(2*pi))*cos(3*ϕ)*cos(θ)*sin(θ)^3
realsphharm44(θ::Float64, ϕ::Float64) = (3/16)*sqrt(35/pi)*cos(4*ϕ)*sin(θ)^4
realsphharm5m5(θ::Float64, ϕ::Float64) = (3/16)*sqrt(77/(2*pi))*sin(5*ϕ)*sin(θ)^5
realsphharm5m4(θ::Float64, ϕ::Float64) = (-(3/16))*sqrt(385/pi)*cos(θ)*sin(4*ϕ)*sin(θ)^4
realsphharm5m3(θ::Float64, ϕ::Float64) = (1/32)*sqrt(385/(2*pi))*(7 + 9*cos(2*θ))*sin(3*ϕ)*sin(θ)^3
realsphharm5m2(θ::Float64, ϕ::Float64) = (-(1/32))*sqrt(1155/pi)*(5*cos(θ) + 3*cos(3*θ))*sin(2*ϕ)*sin(θ)^2
realsphharm5m1(θ::Float64, ϕ::Float64) = (1/256)*sqrt(165/pi)*sin(ϕ)*(2*sin(θ) + 7*(sin(3*θ) + 3*sin(5*θ)))
realsphharm50(θ::Float64, ϕ::Float64) = (1/256)*sqrt(11/pi)*(30*cos(θ) + 35*cos(3*θ) + 63*cos(5*θ))
realsphharm51(θ::Float64, ϕ::Float64) = (1/256)*sqrt(165/pi)*cos(ϕ)*(2*sin(θ) + 7*(sin(3*θ) + 3*sin(5*θ)))
realsphharm52(θ::Float64, ϕ::Float64) = (1/16)*sqrt(1155/pi)*cos(2*ϕ)*cos(θ)*(1 + 3*cos(2*θ))*sin(θ)^2
realsphharm53(θ::Float64, ϕ::Float64) = (1/32)*sqrt(385/(2*pi))*cos(3*ϕ)*(7 + 9*cos(2*θ))*sin(θ)^3
realsphharm54(θ::Float64, ϕ::Float64) = (3/16)*sqrt(385/pi)*cos(4*ϕ)*cos(θ)*sin(θ)^4
realsphharm55(θ::Float64, ϕ::Float64) = (3/16)*sqrt(77/(2*pi))*cos(5*ϕ)*sin(θ)^5
realsphharm6m6(θ::Float64, ϕ::Float64) = (-(1/32))*sqrt(3003/(2*pi))*sin(6*ϕ)*sin(θ)^6
realsphharm6m5(θ::Float64, ϕ::Float64) = (3/16)*sqrt(1001/(2*pi))*cos(θ)*sin(5*ϕ)*sin(θ)^5
realsphharm6m4(θ::Float64, ϕ::Float64) = (-(3/64))*sqrt(91/pi)*(9 + 11*cos(2*θ))*sin(4*ϕ)*sin(θ)^4
realsphharm6m3(θ::Float64, ϕ::Float64) = (1/64)*sqrt(1365/(2*pi))*(21*cos(θ) + 11*cos(3*θ))*sin(3*ϕ)*sin(θ)^3
realsphharm6m2(θ::Float64, ϕ::Float64) = (-(1/256))*sqrt(1365/(2*pi))*(35 + 60*cos(2*θ) + 33*cos(4*θ))*sin(2*ϕ)*sin(θ)^2
realsphharm6m1(θ::Float64, ϕ::Float64) = (1/512)*sqrt(273/pi)*sin(ϕ)*(5*sin(2*θ) + 12*sin(4*θ) + 33*sin(6*θ))
realsphharm60(θ::Float64, ϕ::Float64) = (sqrt(13/pi)*(50 + 105*cos(2*θ) + 126*cos(4*θ) + 231*cos(6*θ)))/1024
realsphharm61(θ::Float64, ϕ::Float64) = (1/512)*sqrt(273/pi)*cos(ϕ)*(5*sin(2*θ) + 12*sin(4*θ) + 33*sin(6*θ))
realsphharm62(θ::Float64, ϕ::Float64) = (1/256)*sqrt(1365/(2*pi))*cos(2*ϕ)*(35 + 60*cos(2*θ) + 33*cos(4*θ))*sin(θ)^2
realsphharm63(θ::Float64, ϕ::Float64) = (1/64)*sqrt(1365/(2*pi))*cos(3*ϕ)*(21*cos(θ) + 11*cos(3*θ))*sin(θ)^3
realsphharm64(θ::Float64, ϕ::Float64) = (3/64)*sqrt(91/pi)*cos(4*ϕ)*(9 + 11*cos(2*θ))*sin(θ)^4
realsphharm65(θ::Float64, ϕ::Float64) = (3/16)*sqrt(1001/(2*pi))*cos(5*ϕ)*cos(θ)*sin(θ)^5
realsphharm66(θ::Float64, ϕ::Float64) = (1/32)*sqrt(3003/(2*pi))*cos(6*ϕ)*sin(θ)^6
realsphharm7m7(θ::Float64, ϕ::Float64) = (3/64)*sqrt(715/pi)*sin(7*ϕ)*sin(θ)^7
realsphharm7m6(θ::Float64, ϕ::Float64) = (-(3/32))*sqrt(5005/(2*pi))*cos(θ)*sin(6*ϕ)*sin(θ)^6
realsphharm7m5(θ::Float64, ϕ::Float64) = (3/128)*sqrt(385/pi)*(11 + 13*cos(2*θ))*sin(5*ϕ)*sin(θ)^5
realsphharm7m4(θ::Float64, ϕ::Float64) = (-(3/64))*sqrt(385/pi)*cos(θ)*(7 + 13*cos(2*θ))*sin(4*ϕ)*sin(θ)^4
realsphharm7m3(θ::Float64, ϕ::Float64) = (3/512)*sqrt(35/pi)*(189 + 308*cos(2*θ) + 143*cos(4*θ))*sin(3*ϕ)*sin(θ)^3
realsphharm7m2(θ::Float64, ϕ::Float64) = (-(3/512))*sqrt(35/(2*pi))*(350*cos(θ) + 275*cos(3*θ) + 143*cos(5*θ))*sin(2*ϕ)*sin(θ)^2
realsphharm7m1(θ::Float64, ϕ::Float64) = (sqrt(105/pi)*sin(ϕ)*(25*sin(θ) + 81*sin(3*θ) + 165*sin(5*θ) + 429*sin(7*θ)))/4096
realsphharm70(θ::Float64, ϕ::Float64) = (sqrt(15/pi)*(175*cos(θ) + 189*cos(3*θ) + 231*cos(5*θ) + 429*cos(7*θ)))/2048
realsphharm71(θ::Float64, ϕ::Float64) = (sqrt(105/pi)*cos(ϕ)*(25*sin(θ) + 81*sin(3*θ) + 165*sin(5*θ) + 429*sin(7*θ)))/4096
realsphharm72(θ::Float64, ϕ::Float64) = (3/512)*sqrt(35/(2*pi))*cos(2*ϕ)*(350*cos(θ) + 275*cos(3*θ) + 143*cos(5*θ))*sin(θ)^2
realsphharm73(θ::Float64, ϕ::Float64) = (3/512)*sqrt(35/pi)*cos(3*ϕ)*(189 + 308*cos(2*θ) + 143*cos(4*θ))*sin(θ)^3
realsphharm74(θ::Float64, ϕ::Float64) = (3/128)*sqrt(385/pi)*cos(4*ϕ)*(27*cos(θ) + 13*cos(3*θ))*sin(θ)^4
realsphharm75(θ::Float64, ϕ::Float64) = (3/128)*sqrt(385/pi)*cos(5*ϕ)*(11 + 13*cos(2*θ))*sin(θ)^5
realsphharm76(θ::Float64, ϕ::Float64) = (3/32)*sqrt(5005/(2*pi))*cos(6*ϕ)*cos(θ)*sin(θ)^6
realsphharm77(θ::Float64, ϕ::Float64) = (3/64)*sqrt(715/pi)*cos(7*ϕ)*sin(θ)^7
realsphharm8m8(θ::Float64, ϕ::Float64) = (-(3/256))*sqrt(12155/pi)*sin(8*ϕ)*sin(θ)^8
realsphharm8m7(θ::Float64, ϕ::Float64) = (3/64)*sqrt(12155/pi)*cos(θ)*sin(7*ϕ)*sin(θ)^7
realsphharm8m6(θ::Float64, ϕ::Float64) = (-(1/128))*sqrt(7293/(2*pi))*(13 + 15*cos(2*θ))*sin(6*ϕ)*sin(θ)^6
realsphharm8m5(θ::Float64, ϕ::Float64) = (3/256)*sqrt(17017/pi)*(11*cos(θ) + 5*cos(3*θ))*sin(5*ϕ)*sin(θ)^5
realsphharm8m4(θ::Float64, ϕ::Float64) = -((3*sqrt(1309/pi)*(99 + 156*cos(2*θ) + 65*cos(4*θ))*sin(4*ϕ)*sin(θ)^4)/1024)
realsphharm8m3(θ::Float64, ϕ::Float64) = (sqrt(19635/pi)*(126*cos(θ) + 91*cos(3*θ) + 39*cos(5*θ))*sin(3*ϕ)*sin(θ)^3)/1024
realsphharm8m2(θ::Float64, ϕ::Float64) = -((3*sqrt(595/(2*pi))*(210 + 385*cos(2*θ) + 286*cos(4*θ) + 143*cos(6*θ))*sin(2*ϕ)*sin(θ)^2)/2048)
realsphharm8m1(θ::Float64, ϕ::Float64) = (3*sqrt(17/pi)*sin(ϕ)*(70*sin(2*θ) + 154*sin(4*θ) + 286*sin(6*θ) + 715*sin(8*θ)))/8192
realsphharm80(θ::Float64, ϕ::Float64) = (sqrt(17/pi)*(1225 + 2520*cos(2*θ) + 2772*cos(4*θ) + 3432*cos(6*θ) + 6435*cos(8*θ)))/32768
realsphharm81(θ::Float64, ϕ::Float64) = (3*sqrt(17/pi)*cos(ϕ)*(70*sin(2*θ) + 154*sin(4*θ) + 286*sin(6*θ) + 715*sin(8*θ)))/8192
realsphharm82(θ::Float64, ϕ::Float64) = (3*sqrt(595/(2*pi))*cos(2*ϕ)*(210 + 385*cos(2*θ) + 286*cos(4*θ) + 143*cos(6*θ))*sin(θ)^2)/2048
realsphharm83(θ::Float64, ϕ::Float64) = (sqrt(19635/pi)*cos(3*ϕ)*(126*cos(θ) + 91*cos(3*θ) + 39*cos(5*θ))*sin(θ)^3)/1024
realsphharm84(θ::Float64, ϕ::Float64) = (3*sqrt(1309/pi)*cos(4*ϕ)*(99 + 156*cos(2*θ) + 65*cos(4*θ))*sin(θ)^4)/1024
realsphharm85(θ::Float64, ϕ::Float64) = (3/256)*sqrt(17017/pi)*cos(5*ϕ)*(11*cos(θ) + 5*cos(3*θ))*sin(θ)^5
realsphharm86(θ::Float64, ϕ::Float64) = (1/128)*sqrt(7293/(2*pi))*cos(6*ϕ)*(13 + 15*cos(2*θ))*sin(θ)^6
realsphharm87(θ::Float64, ϕ::Float64) = (3/64)*sqrt(12155/pi)*cos(7*ϕ)*cos(θ)*sin(θ)^7
realsphharm88(θ::Float64, ϕ::Float64) = (3/256)*sqrt(12155/pi)*cos(8*ϕ)*sin(θ)^8

i_lm(l::Int, m::Int) = l*(l+1)+1+m

sphharmtab = [
    sphharm00,
    sphharm1m1,
    sphharm10,
    sphharm11,
    sphharm2m2,
    sphharm2m1,
    sphharm20,
    sphharm21,
    sphharm22,
    sphharm3m3,
    sphharm3m2,
    sphharm3m1,
    sphharm30,
    sphharm31,
    sphharm32,
    sphharm33,
    sphharm4m4,
    sphharm4m3,
    sphharm4m2,
    sphharm4m1,
    sphharm40,
    sphharm41,
    sphharm42,
    sphharm43,
    sphharm44,
    sphharm5m5,
    sphharm5m4,
    sphharm5m3,
    sphharm5m2,
    sphharm5m1,
    sphharm50,
    sphharm51,
    sphharm52,
    sphharm53,
    sphharm54,
    sphharm55,
    sphharm6m6,
    sphharm6m5,
    sphharm6m4,
    sphharm6m3,
    sphharm6m2,
    sphharm6m1,
    sphharm60,
    sphharm61,
    sphharm62,
    sphharm63,
    sphharm64,
    sphharm65,
    sphharm66,
    sphharm7m7,
    sphharm7m6,
    sphharm7m5,
    sphharm7m4,
    sphharm7m3,
    sphharm7m2,
    sphharm7m1,
    sphharm70,
    sphharm71,
    sphharm72,
    sphharm73,
    sphharm74,
    sphharm75,
    sphharm76,
    sphharm77,
    sphharm8m8,
    sphharm8m7,
    sphharm8m6,
    sphharm8m5,
    sphharm8m4,
    sphharm8m3,
    sphharm8m2,
    sphharm8m1,
    sphharm80,
    sphharm81,
    sphharm82,
    sphharm83,
    sphharm84,
    sphharm85,
    sphharm86,
    sphharm87,
    sphharm88
]

realsphharmtab = [
    realsphharm00,
    realsphharm1m1,
    realsphharm10,
    realsphharm11,
    realsphharm2m2,
    realsphharm2m1,
    realsphharm20,
    realsphharm21,
    realsphharm22,
    realsphharm3m3,
    realsphharm3m2,
    realsphharm3m1,
    realsphharm30,
    realsphharm31,
    realsphharm32,
    realsphharm33,
    realsphharm4m4,
    realsphharm4m3,
    realsphharm4m2,
    realsphharm4m1,
    realsphharm40,
    realsphharm41,
    realsphharm42,
    realsphharm43,
    realsphharm44,
    realsphharm5m5,
    realsphharm5m4,
    realsphharm5m3,
    realsphharm5m2,
    realsphharm5m1,
    realsphharm50,
    realsphharm51,
    realsphharm52,
    realsphharm53,
    realsphharm54,
    realsphharm55,
    realsphharm6m6,
    realsphharm6m5,
    realsphharm6m4,
    realsphharm6m3,
    realsphharm6m2,
    realsphharm6m1,
    realsphharm60,
    realsphharm61,
    realsphharm62,
    realsphharm63,
    realsphharm64,
    realsphharm65,
    realsphharm66,
    realsphharm7m7,
    realsphharm7m6,
    realsphharm7m5,
    realsphharm7m4,
    realsphharm7m3,
    realsphharm7m2,
    realsphharm7m1,
    realsphharm70,
    realsphharm71,
    realsphharm72,
    realsphharm73,
    realsphharm74,
    realsphharm75,
    realsphharm76,
    realsphharm77,
    realsphharm8m8,
    realsphharm8m7,
    realsphharm8m6,
    realsphharm8m5,
    realsphharm8m4,
    realsphharm8m3,
    realsphharm8m2,
    realsphharm8m1,
    realsphharm80,
    realsphharm81,
    realsphharm82,
    realsphharm83,
    realsphharm84,
    realsphharm85,
    realsphharm86,
    realsphharm87,
    realsphharm88
]

max_stars_repo_name: johncsnyder/Element


text:
@testset "trace" begin
    @testset "basic" begin
        @test isnothing(current_span())
        @test isnothing(span_context())

        @testset "TraceState" begin
            ts = TraceState(
                "1a-2f@foo" => "bar1",
                "1a-_*/2b@foo" => "bar2",
                "foo" => "bar3",
                "foo-_*/bar" => "bar4",
                "^foo-_*/bar" => "bar5",  # invalid key
                "hello" => "世界！",  # invalid value
            )

            @test string(ts) == "1a-2f@foo=bar1,1a-_*/2b@foo=bar2,foo=bar3,foo-_*/bar=bar4"
            @test length(ts) == 4
            @test haskey(ts, :foo)
            @test !haskey(ts, "^foo-_*/bar")
            @test ts[:foo] == "bar3"
        end
    end

    @testset "TracerProvider" begin
        global_tracer_provider!(global_tracer_provider())  # smoke test

        tracer = Tracer()
        with_span("test", tracer) do
            s = current_span()
            @test resource(s) == resource(global_tracer_provider())
            @test span_context() === INVALID_SPAN_CONTEXT

            s["foo"] = "bar"

            @test is_recording() == false
            @test haskey(s, "foo") == false  # not recording
            @test_throws Exception s["foo"]
            @test length(span_events()) == 0
            @test length(span_links()) == 0
            @test span_status().code == SPAN_STATUS_UNSET

            span_name!("XYZ")
            @test span_name() == "test"
            @info "Other APIs" span_kind() parent_span_context()
        end

        with_span("foo", tracer) do
            with_span("bar", tracer) do
                with_span("baz", tracer) do
                    println("hello")
                end
                with_span("baz", tracer) do
                    println("world!")
                end
            end
        end
    end
end


text:
using DelimitedFiles
using Plots

#open("./examples/ES_Fourier-DG/array1.txt","r") do io
convarr = readdlm("./examples/ES_Fourier-DG/conv/trans_vortex_y/array_conv.txt",'\t',Float64,'\n')
convarr = reshape(convarr,4,1,5)

N_P_arr = [1;2;3;4;5]
Np_F_arr = [8]
K1D_arr = [2;4;10;20]

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
function devsetup()
    Base.depwarn("WebIO.devsetup() is no longer required. You can simply run WebIO.bundlejs().", :webio_devesetup)
end

function bundlejs()
    include(joinpath(@__DIR__, "..", "deps", "_bundlejs.jl"))
end


text:
# Avoid variables, constraints and objective being added when loading
@enum State Loading Modeling

mutable struct Data
    variables::Set{SetVariableRef}
    constraints::Dict{ConstraintIndex, SetConstraint}
    transformed_constraints::Dict{ConstraintIndex, JuMP.ConstraintRef}
    names::Dict{ConstraintIndex, String}
    last_index::Int
    objective_sense::MOI.OptimizationSense
    objective::Union{Nothing, AbstractScalarFunction}
    objective_variable::Union{Nothing, JuMP.VariableRef}
    perspective_polyvar::Union{Nothing, SpaceVariable}
    space::Space
    state::State
    spaces::Union{Nothing, Spaces}
end

function JuMP.copy_extension_data(::SetProg.Data, ::JuMP.Model, ::JuMP.Model)
end

function lift_space_variables(data::Data,
                              x::Vector{SpaceVariable})
    @assert data.perspective_polyvar > x[1]
    SpaceVariable[data.perspective_polyvar; x]
end

function set_space(cur::Space, space::Space)
    if space == Undecided || cur == Undecided || cur == space
        return space
    else
        error("Incompatible constraints/objective, some require to do the modeling in the primal space and some in the dual space.")
    end
end
function set_space(d::Data, model::JuMP.Model)
    space = Undecided
    for (index, constraint) in d.constraints
        space = set_space(space, constraint)
    end
    if d.objective !== nothing
        space = set_space(space, d.objective, model)
    end
    if space == Undecided
        space = PrimalSpace
    end
    d.space = space
end

# Remove spaces creating in earlier `optimize!`
function clear_spaces(d::Data)
    for vref in d.variables
        clear_spaces(vref)
    end
    for (index, constraint) in d.constraints
        clear_spaces(constraint)
    end
end
# No space index stored in constants
function clear_spaces(::Union{Polyhedra.Rep, Sets.AbstractSet}) end

# No perspective variable
function Sets.perspective_variable(::Polyhedra.Rep) end

synchronize_perspective(::Nothing, ::Nothing) = nothing
synchronize_perspective(::Nothing, z::SpaceVariable) = z
synchronize_perspective(z::SpaceVariable, ::Nothing) = z
function synchronize_perspective(z1::SpaceVariable, z2::SpaceVariable)
    if z1 !== z2
        throw(ArgumentError("Perspective variables do not match"))
    end
    return z1
end
function synchronize_perspective(d::Data)
    z = nothing
    for (index, constraint) in d.constraints
        z = synchronize_perspective(z, Sets.perspective_variable(constraint))
    end
    if z === nothing
        DynamicPolynomials.@polyvar z
    end
    d.perspective_polyvar = z
end

function create_spaces(set::Polyhedra.Rep, spaces::Spaces)
    return new_space(spaces, Polyhedra.fulldim(set))
end
function create_spaces(set::Union{Sets.AbstractSet, Sets.Projection}, spaces::Spaces)
    polyvars = Sets.space_variables(set)
    if polyvars === nothing
        return new_space(spaces, Sets.dimension(set))
    else
        return new_space(spaces, polyvars)
    end
end

# Create spaces and identify objects lying the the same space
function create_spaces(d::Data)
    spaces = Spaces()
    for vref in d.variables
        create_spaces(vref, spaces)
    end
    for (index, constraint) in d.constraints
        create_spaces(constraint, spaces)
    end
    d.spaces = spaces
end

function data(model::JuMP.Model)
    if !haskey(model.ext, :SetProg)
        model.ext[:SetProg] = Data(Set{SetVariableRef}(),
                                   Dict{ConstraintIndex, SetConstraint}(),
                                   Dict{ConstraintIndex, JuMP.ConstraintRef}(),
                                   Dict{ConstraintIndex, String}(), 0,
                                   MOI.FEASIBILITY_SENSE, nothing, nothing,
                                   nothing, Undecided, Modeling, nothing)
        model.optimize_hook = optimize_hook
    end
    return model.ext[:SetProg]
end


text:
using LinearAlgebra

"""Optimize some function f
# Args:
    ∇f: Gradient of f
    x_0: Initial guess
    η: learning rate
    ε: cancellation tolerance
    max_iters: maximum number of iterations
"""
function gradient_descent(∇f::Function,
  x_0,
  η::Number,
  ε::Number,
  max_iters::Integer)
    x_k = x_0
    count = 0
    for _ = 1:max_iters
        if any(isnan.(x_k))
            error("Encoutered NaN")
        end
        if any(isinf.(x_k))
            error("Encoutered Inf")
        end
        count += 1
        p_k = -∇f(x_k)
        if norm(∇f(p_k)) < ε
            break
        end
        x_k += η * p_k
    end
    x_k
end

#numeric_differentiation(f, h) = x -> (f(x + h) - f(x)) / h
numeric_differentiation(f, h) = x -> (f(x + h) - f(x - h)) / (2 * h)

"""BFGS Optimization Algorithm
Broyden–Fletcher–Goldfarb–Shanno algorithm
for details see:
    https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm

# Args:
    ∇f: Gradient of function to optimize
    x_0: Initial guess for optimal value
    iters: Maximum number of iterations before cancellation
    line_search: Line search function that's used
    ε: Optimization stops once the norm of the gradient is below this value
    ls_x_0: Initial value for the line search
    ls_η: Learning rate for the line search
    ε_ls: same as ε but for the line search
"""
function BFGS(f::Function,
    ∇f::Function,
    x_0::Vector{<:Number},
    iters::Integer,
    line_search::Function,
    ε = 10e-12::Real,
    )
    n = size(x_0)[1]
    x_k = x_0
    B_k = Matrix{typeof(x_0[1])}(I, n, n)
    did_iters = 0
    for i = 1:iters
        did_iters += 1
        if norm(∇f(x_k)) < ε
            break
        end
        # Step 1: obtain direction p_k by solving B_k ∙ p_k = - (gradient of f at x_k)
        p_k = B_k \ -∇f(x_k)

        # Step 2.: Find stepsize α_k such that α_k = arg min ∂f(x_k + α_k * p_k)/∂α
        α_k = line_search(numeric_differentiation(α->(f(x_k + α * p_k)), 10e-10))
        if isnan(α_k)
            break
        end

        # Step 3.
        s_k = α_k * p_k
        x_k_prime = x_k + s_k

        # Step 4.
        y_k = ∇f(x_k_prime) - ∇f(x_k)

        # Step 5.
        B_k += (y_k * y_k') / (y_k' * s_k) - (B_k * s_k * s_k' * B_k) / (s_k' * B_k * s_k)
        x_k = x_k_prime
    end
    println(did_iters)
    x_k
end

# f(x) = -5 * x[1] + 10 * x[1]^2
# ∇f(x) = [-5 + 20 * x[1]]
f(x) = x[1]^5 / 5000 + 21 * x[1]^4 / 4000 + 17 * x[1]^3 / 375 + 293 * x[1]^2 / 1000 + 521 * x[1] / 1000
∇f(x) = [x[1]^4 / 1000 + 21 * x[1]^3 / 1000 + 17 * x[1]^2 / 125 + 293 * x[1] / 500 + 521 / 1000]
a = 1
b = 100
r(X) = begin
    x = X[1]
    y = X[2]
    (a - x)^2 + b*(y - x^2)^2
end
∇r(X) = begin
    x = X[1]
    y = X[2]
    [-2*a+4*b*x^3-4*b*x*y+2*x, 2*b*(y-x^2)]
end
#f(X) = begin
#    x = X[1]
#    y = X[2]
#    x^2 + y^2
#end
#∇f(X) = begin
#    x = X[1]
#    y = X[2]
#    [2*x, 2*y]
#end
println("optimal x = ", BFGS(
    f,
    ∇f,
    [0],
    500,
    ∇f -> gradient_descent(∇f, 10e-10, 1, 1e-10, 1000),
    10e-5))

# f(w, x) = w[1] * x[1] + w[2] * x[1]^2

# ∇f(w, x) = [w[1] + 2 * x[1] * w[2]] # [x; x^2]

# # initialize parameters
# w_0 = [-5, 10]
# x_0 = [-5]
# B_0 = Matrix{Float64}(I, 1, 1) # initialize approximate hessian with identity matrix
# ε = 1e-10 # cancellation parameter

# w_k = w_0
# x_k = x_0
# B_k = B_0

# while begin
#     x = norm(∇f(w_k, x_k))
#     println("∇f = ", x)
#     x > ε
# end
#     global B_k
#     global x_k
#     # Step 1.: obtain direction p_k by solving B_k*p_k = - (gradient of f at x_k)

#     p_k = B_k \ -∇f(w_k, x_k)

#     # Step 2.: Find stepsize α_k such that α_k minimizes f(x_k + α_k * p_k)

#     g_k(α) = w_k[1] + 2 * (x_k[1] + α * p_k[1]) * w_k[2]
#     α_k = line_search_gradient_descent(g_k, 0, 0.000001, 1e-10)
#     println("α_k = ", α_k)

#     # Step 3.
#     s_k = α_k * p_k
#     x_k_prime = x_k + s_k

#     # Step 4.
#     y_k = ∇f(w_k, x_k_prime) - ∇f(w_k, x_k)

#     # Step 5.
#     B_k += (y_k * y_k') / (y_k' * s_k) - (B_k * s_k * s_k' * B_k) / (s_k' * B_k * s_k)
#     x_k = x_k_prime
# end

# println(x_k)

"""
begin
    inv_B = inv(B_k)
    (s_k' * y_k + y_k' * inv_B * y_k) * (s_k * s_k') / (s_k' * y_k)^2 - (inv_B * y_k * s_k' + s_k * y_k' * inv_B) / (s_k' * y_k)
end
"""


text:
using DataDeps
using Test
import DataDeps.get_filename

@testset "post processing" begin
    phf(x) = DataDeps.process_header_filename(match(r"(.*)",x))

    @test phf("foo.txt") == "foo.txt"
    @test phf("\"foo.txt\"") == "foo.txt"
    @test phf("\"bar \\\"zap zing\\\" foo.txt\"") == "bar \"zap zing\" foo.txt"
end

@testset "using headers" begin
    # From http://test.greenbytes.de/tech/tc2231/

    @test "foo.html" == get_filename("http://test.greenbytes.de/tech/tc2231/attwithasciifilename.asis")

    @test "foo.html" == get_filename("http://test.greenbytes.de/tech/tc2231/attwithasciifilenamenq.asis") # Technically invalid, but github API does this

    @test "attonly.asis" == get_filename("http://test.greenbytes.de/tech/tc2231/attonly.asis")

    @test "\"quoting\" tested.html" == get_filename("http://test.greenbytes.de/tech/tc2231/attwithasciifnescapedquote.asis")

    @test_broken "ä-%41.html" == get_filename("http://test.greenbytes.de/tech/tc2231/attwithfilenamepctandiso.asis")

end

@testset "not requiring headers" begin
    @test get_filename("https://www.angio.net/pi/digits/10000.txt") == "10000.txt"
end


text:
module DashBio
using Dash

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "0.7.0"

include("dashbio_alignmentchart.jl")
include("dashbio_fornacontainer.jl")
include("dashbio_ideogram.jl")
include("dashbio_igv.jl")
include("dashbio_molecule2dviewer.jl")
include("dashbio_molecule3dviewer.jl")
include("dashbio_needleplot.jl")
include("dashbio_nglmoleculeviewer.jl")
include("dashbio_oncoprint.jl")
include("dashbio_sequenceviewer.jl")
include("dashbio_speck.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "dash_bio",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "async-alignment.js",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-alignment.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-circos.js",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-circos.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-ideogram.js",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-ideogram.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-igv.js",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-igv.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-pileup.js",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-pileup.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-moleculeviewer2.js",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-moleculeviewer2.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-moleculeviewer3.js",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-moleculeviewer3.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-needle.js",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-needle.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-nglmoleculeviewer.js",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-nglmoleculeviewer.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-onco.js",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-onco.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-sequence.js",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-sequence.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-speck.js",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-speck.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-alignment.js.map",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-alignment.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-circos.js.map",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-circos.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-ideogram.js.map",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-ideogram.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-igv.js.map",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-igv.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-pileup.js.map",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-pileup.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-moleculeviewer2.js.map",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-moleculeviewer2.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-moleculeviewer3.js.map",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-moleculeviewer3.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-needle.js.map",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-needle.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-nglmoleculeviewer.js.map",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-nglmoleculeviewer.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-onco.js.map",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-onco.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-sequence.js.map",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-sequence.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-speck.js.map",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/async-speck.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "bundle.js",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/bundle.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "bundle.js.map",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/bundle.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_bio-shared.js",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/dash_bio-shared.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_bio-shared.js.map",
    external_url = "https://unpkg.com/dash-bio@0.7.0/dash_bio/dash_bio-shared.js.map",
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end


text:
using DashBootstrapComponents

progress = dbc_progress(label = "25%", value = 25);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
""" Cálculo de movimento retilíneo
# Prepare um algoritmo para calcular
# o espaço percorrido por um móvel em movimento retilíneo uniforme dada a seguinte fórmula:
# S=So+V*T.
# Inicialmente o algoritmo deverá perguntar ao usuário
# os valores do espaço inicial(So),
# da velocidade(V)
# e do tempo(T).
# Após a entrada dos devidos valores estes devem ser substituídos na fórmula,
# resultando assim no espaço percorrido(S) que deverá ser informado na tela do usuário.
"""

println("informe o valor da velocidade")
mru = parse(Float32, readline())

println("informe o valor do tempo")
mru *= parse(Float32, readline())

println("informe o valor do espaço inicial")
mru += parse(Float32, readline())

println("espaço percorrido = ", mru)

max_stars_repo_name: codinginbrazil/functional-programming

max_issues_repo_name: codinginbrazil/functional-programming

max_forks_repo_name: codinginbrazil/functional-programming


text:
using BaseExtensions
using Test
using Documenter
using Aqua

@testset "project quality" begin
    Aqua.test_all(BaseExtensions)
end

DocMeta.setdocmeta!(BaseExtensions, :DocTestSetup, :(using BaseExtensions); recursive=true)

@testset "docstrings" begin
    doctest(BaseExtensions)
end

@testset "anynonzero + allzero" begin
    a = [1,2,3];
    @test anynonzero(a)
    @test !allzero(a)

    a = [1,2,0];
    @test anynonzero(a)
    @test !allzero(a)

    a = zeros(3);
	@test !anynonzero(a)
	@test allzero(a)
end

@testset "zero!" begin
    a = [1, 2, 3]
    zero!(a)
    @test allzero(a)
end

@testset "realtype" begin
    c = ComplexF64(3,2)
    @test realtype(c) == Float64
    @test realtype(ComplexF64) == Float64

    c = ComplexF32(3,2)
    @test realtype(c) == Float32
    @test realtype(ComplexF32) == Float32
end

@testset "reinterpretreal" begin
	for T in [Int, Float64]
	    a = Complex{T}[1 + im, 2 + 2im]
	    @test reinterpretreal(a) == T[1, 1, 2, 2]
	end
end

@testset "reinterpretcomplex" begin
	for T in [Int, Float64]
	    a = T[1, 1, 2, 2]
	    @test reinterpretcomplex(a) == Complex{T}[1 + im, 2 + 2im]
	end
end

@testset "squeeze" begin
    A = rand(2,1,2)
    B = squeeze(A)
    @test B == reshape(A, 2, 2)
end

@testset "flipdims" begin
    A = rand(1,2,3)
    B = flipdims(A)
    @test B == permutedims(A, [3,2,1])
    @test axes(B) == reverse(axes(A))
end

@testset "filldeepcopy" begin
    A = Vector{Vector{Vector{Int}}}(undef, 3);
    filldeepcopy!(A, [ones(Int, 2) for i = 1:2]);

    A[1][1] .= 3:4;

    for i in 1:3, j = 1:2
        if i == j == 1
            @test all(A[i][j] .== 3:4)
        else
            @test all(A[i][j] .== 1)
        end
    end
end

@testset "fillcopy" begin
    A = Vector{Vector{Vector{Int}}}(undef, 3);
    fillcopy!(A, [ones(Int, 2) for i = 1:2]);

    A[1][1] .= 3:4

    for i in 1:3, j = 1:2
        if j == 1
            @test all(A[i][j] .== 3:4)
        else
            @test all(A[i][j] .== 1)
        end
    end
end


text:
porkchop_layout(
    ) = Layout(# autosize=true, width=800, height=600,
            autosize=true,
            paper_bgcolor = "rgb(30, 30, 30)",
            plot_bgcolor="rgb(30, 30, 30)",
            font = attr(family = "Courier New, monospace",
                        size = 10,
                        color = "rgb(200, 200, 200)"
            ),
            xaxis=attr(title_text="Transfer start (day #)",
                       showspikes=true,
                       spikemode="across",
                       spikecolor="rgb(200, 200, 200)",
                       spikedash="solid",
                       spikethickness=-1,
                       showgrid=false,
            ),
            yaxis=attr(title_text="Transfer duration (days)",
                       showspikes=true,
                       spikemode="across",
                       spikecolor="rgb(200, 200, 200)",
                       spikedash="solid",
                       spikethickness=-1,
                       showgrid=false,
            ),
)

"""
    trace = draw_porkchop(pc; kwargs...)

Creates a trace of a porkchop plot.
"""
function draw_porkchop(pc::Porkchop; levelscale=1.1, nlevels=16, timedef=KerbalTime, plottype=:Δv)
    minΔv = minimum(getfield(pc, plottype))
    lvls = minΔv * [levelscale .^ i for i=0:nlevels]
    logΔv = transpose(broadcast(x-> log(x)/log(levelscale), getfield(pc, plottype)))
    loglvls = [log(lvl)/log(levelscale) for lvl in lvls]
    colorlabels = ["$(Int(floor(lvl)))" for lvl in lvls]

    day_to_sec = timedef["day"] * timedef["hour"] * timedef["minute"]

    trace = contour(
        z = logΔv,
        x = pc.starttimes ./ day_to_sec,
        y = pc.flighttimes ./ day_to_sec,
        contours = Dict(
            :start => loglvls[1]+0.05,
            :end => loglvls[end],
            :size => (loglvls[end]-loglvls[1])/nlevels,
            # :showlines => false,
        ),
        colorscale = "Viridis",
        reversescale = true,
        colorbar = attr(
            tickvals = loglvls,
            ticktext = colorlabels,
            title = attr(
                text = "Δv (m/s)",
                font = attr(
                    family = "Courier New, monospace",
                ),
            ),
            tickfont = attr(
                family = "Courier New, monospace",
            ),
        ),
        customdata = getfield(pc, plottype),
        hovertemplate = "Δv = %{customdata:.2f} m/s<extra></extra>",
    )
    return trace
end


text:
function serial_connect(syss::Vector)
    n_sys = length(syss)
    @assert n_sys > 1 "Connect at least two components"

    connections = [
        connect(syss[i].out, syss[i+1].in)
        for i in 1:(n_sys-1)
    ]

    return connections
end

serial_connect(syss...) = serial_connect([syss...])


text:
module Calibration

using Roots
using NLsolve
using Distributions
using Macros
using CustomTypes

export calibrate_y, calibrate_m, calibrate_ssξ, calibrate_u, calibrate_c

@def cal_common begin
    μ_z = 0.
    Ez(α) = exp(μ_z + α^2/2)
    G(x,α) = cdf(LogNormal(μ_z,α),x)
    I(x,α) = exp(μ_z+α^2/2)*cdf(Normal(0,1),((μ_z+α^2-log(max(x,0)))/α))
    I0(x,α) = I(x,α) - x*(1-G(x,α))

    zp(α) = fzero( zp -> (ξ-s)/λ - G(zp,α) , 0, Ez(α)+3*α )
    rU(α) = zp(α)+(r+s)*F+λ*I0(zp(α),α)/(r+s+λ)
    zf(α) = rU(α) - λ*(Ez(α)-rU(α))/(r+δ)
    zs(α) = (zp(α)/(r+s+λ) + F - zf(α)/(r+δ+λ))/(1/(r+s+λ) - 1/(r+δ+λ))

    share_tc_jc(α) = (G(zs(α),α)-G(zf(α),α))/(1-G(zf(α),α))

    α = fzero(α->share_tc_jc(α) - μᶠsμ, 0.01, 1.)

    Ez_cal = Ez(α)
    G_cal(x) = G(x,α)
    I_cal(x) = I(x,α)
    I0_cal(x) = I0(x,α)

    zp_cal = zp(α)
    rU_cal = rU(α)
    zf_cal = zf(α)
    zs_cal = zs(α)

    p = μᵖ/(1-G_cal(zs_cal))
    γ = q*(1-η)*( I0_cal(zs_cal)/(r+s+λ) + (I0_cal(zf_cal)-I0_cal(zs_cal))/(r+δ+λ) )
    θ = p/q
    b = rU_cal - η*γ*θ/(1-η)
    m = q*θ^σ

    wf_bar = η*( λ*Ez_cal + δ*(I_cal(zf_cal) - I_cal(zs_cal))/(G_cal(zs_cal)-G_cal(zf_cal)) )/(δ+λ) + (1-η)*rU_cal

    E_wp = η*(I_cal(zp_cal)/(1-G_cal(zp_cal)) + (r+s)*F) + (1-η)*rU_cal
    E_wp_0 = η*(I_cal(zs_cal)/(1-G_cal(zs_cal)) - λ*F) + (1-η)*rU_cal
    wp_bar = (λ*(1-G_cal(zp_cal))*E_wp + ξ*E_wp_0)/(s+λ)

end

"""
    calibrate_y()

Returns a dictionary linking parameters and their values when the average time interval between two productivity shocks is one year
"""
function calibrate_y()
    η = 0.6 #Worker's bargaining power
    σ = 0.6 #Elasticity of the matching function w.r.t vacancies
    r = 0.05/12 #Monthly interest rate
    δ = 1/1.5 #Probability of temporary jobs' destruction
    rho_F = 4/3 #Firing-cost-to-average-permanent-contract-wage ratio
    λ = 1/12 #Probability of occurence of productivity shocks

    u = 0.26 #Unemployment
    nfsn = 0.12 #Temporary-to-total-employment ratio
    q = 1-(1-0.7)^(1/3) #Firms' worker-meeting probability
    ssξ = 0.7 #Exogenous to total permanent splits
    μᶠsμ = 0.83 #Share of temporary jobs in job creation

    n = 1-u
    nf = nfsn*n
    np = n-nf

    μᶠ = δ*nf/u
    μ = μᶠ/μᶠsμ
    μᵖ = μ-μᶠ

    ξ = μᵖ*u/np
    s = ssξ*ξ

    n = 1-u
    nf = nfsn*n
    np = n-nf

    μᶠ = δ*nf/u
    μ = μᶠ/μᶠsμ
    μᵖ = μ-μᶠ

    ξ = μᵖ*u/np
    s = ssξ*ξ

    function Eq(F)
        @cal_common
        return (F/wp_bar)-rho_F
    end

    F = fzero(Eq,1.2,2.)

    @cal_common

    cal = Parameters()
    Macros.@save_struct cal η r s λ α m b γ σ F δ
    return cal

end

"""

    calibrate_m()

Returns a dictionary linking parameters and their values when the average time interval between two productivity shocks is one month.
"""
function calibrate_m()

    η = 0.6 #Worker's bargaining power
    σ = 0.6 #Elasticity of the matching function w.r.t vacancies
    r = 0.05/12 #Monthly interest rate
    δ = 1/1.5 #Probability of temporary jobs' destruction
    rho_F = 4/3 #Firing-cost-to-average-permanent-contract-wage ratio
    λ = 1 #Probability of occurence of productivity shocks

    u = 0.26 #Unemployment
    nfsn = 0.12 #Temporary-to-total-employment ratio
    q = 1-(1-0.7)^(1/3) #Firms' worker-meeting probability
    ssξ = 0.7 #Exogenous to total permanent splits
    μᶠsμ = 0.83 #Share of temporary jobs in job creation

    n = 1-u
    nf = nfsn*n
    np = n-nf

    μᶠ = δ*nf/u
    μ = μᶠ/μᶠsμ
    μᵖ = μ-μᶠ

    ξ = μᵖ*u/np
    s = ssξ*ξ

    n = 1-u
    nf = nfsn*n
    np = n-nf

    μᶠ = δ*nf/u
    μ = μᶠ/μᶠsμ
    μᵖ = μ-μᶠ

    ξ = μᵖ*u/np
    s = ssξ*ξ

    function Eq(F)
        @cal_common
        return (F/wp_bar)-rho_F
    end

    F = fzero(Eq,1.2,2.)

    @cal_common

    cal = Parameters()
    Macros.@save_struct cal η r s λ α m b γ σ F δ
    return cal

end

"""
    calibrate_c()

Returns a dictionary linking parameters and their values when the average time interval between two productivity shocks is three years.
"""
function calibrate_c()

    η = 0.6 #Worker's bargaining power
    σ = 0.6 #Elasticity of the matching function w.r.t vacancies
    r = 0.05/12 #Monthly interest rate
    δ = 1/1.5 #Probability of temporary jobs' destruction
    rho_F = 4/3 #Firing-cost-to-average-permanent-contract-wage ratio
    λ = 1/36 #Probability of occurence of productivity shocks

    u = 0.26 #Unemployment
    nfsn = 0.12 #Temporary-to-total-employment ratio
    q = 1-(1-0.7)^(1/3) #Firms' worker-meeting probability
    ssξ = 0.7 #Exogenous to total permanent splits
    μᶠsμ = 0.83 #Share of temporary jobs in job creation

    n = 1-u
    nf = nfsn*n
    np = n-nf

    μᶠ = δ*nf/u
    μ = μᶠ/μᶠsμ
    μᵖ = μ-μᶠ

    ξ = μᵖ*u/np
    s = ssξ*ξ

    n = 1-u
    nf = nfsn*n
    np = n-nf

    μᶠ = δ*nf/u
    μ = μᶠ/μᶠsμ
    μᵖ = μ-μᶠ

    ξ = μᵖ*u/np
    s = ssξ*ξ

    function Eq(F)
        @cal_common
        return (F/wp_bar)-rho_F
    end

    F = fzero(Eq,1.2,2.)

    @cal_common

    cal = Parameters()
    Macros.@save_struct cal η r s λ α m b γ σ F δ
    return cal

end

"""
    calibrate_ssξ()

Returns a dictionary linking parameters and their values when the average time interval between two productivity shocks is one year and s = 0.
"""
function calibrate_ssξ()

    η = 0.6 #Worker's bargaining power
    σ = 0.6 #Elasticity of the matching function w.r.t vacancies
    r = 0.05/12 #Monthly interest rate
    δ = 1/1.5 #Probability of temporary jobs' destruction
    rho_F = 4/3 #Firing-cost-to-average-permanent-contract-wage ratio
    λ = 1/12 #Probability of occurence of productivity shocks

    u = 0.26 #Unemployment
    nfsn = 0.12 #Temporary-to-total-employment ratio
    q = 1-(1-0.7)^(1/3) #Firms' worker-meeting probability
    ssξ = 0. #Exogenous to total permanent splits
    μᶠsμ = 0.83 #Share of temporary jobs in job creation

    n = 1-u
    nf = nfsn*n
    np = n-nf

    μᶠ = δ*nf/u
    μ = μᶠ/μᶠsμ
    μᵖ = μ-μᶠ

    ξ = μᵖ*u/np
    s = ssξ*ξ

    n = 1-u
    nf = nfsn*n
    np = n-nf

    μᶠ = δ*nf/u
    μ = μᶠ/μᶠsμ
    μᵖ = μ-μᶠ

    ξ = μᵖ*u/np
    s = ssξ*ξ

    function Eq(F)
        @cal_common
        return (F/wp_bar)-rho_F
    end

    F = fzero(Eq,1.2,2.)

    @cal_common

    cal = Parameters()
    Macros.@save_struct cal η r s λ α m b γ σ F δ
    return cal

end

"""
    calibrate_u()

Returns a dictionary linking parameters and their values when the average time interval between two productivity shocks is one year and the job seekers' pool does not include inactive workers.
"""
function calibrate_u()
    η = 0.6 #Worker's bargaining power
    σ = 0.6 #Elasticity of the matching function w.r.t vacancies
    r = 0.05/12 #Monthly interest rate
    δ = 1/1.5 #Probability of temporary jobs' destruction
    rho_F = 4/3 #Firing-cost-to-average-permanent-contract-wage ratio
    λ = 1/12 #Probability of occurence of productivity shocks

    u = 0.1 #Unemployment
    nfsn = 0.12 #Temporary-to-total-employment ratio
    q = 1-(1-0.7)^(1/3) #Firms' worker-meeting probability
    ssξ = 0.7 #Exogenous to total permanent splits
    μᶠsμ = 0.83 #Share of temporary jobs in job creation

    n = 1-u
    nf = nfsn*n
    np = n-nf

    μᶠ = δ*nf/u
    μ = μᶠ/μᶠsμ
    μᵖ = μ-μᶠ

    ξ = μᵖ*u/np
    s = ssξ*ξ

    n = 1-u
    nf = nfsn*n
    np = n-nf

    μᶠ = δ*nf/u
    μ = μᶠ/μᶠsμ
    μᵖ = μ-μᶠ

    ξ = μᵖ*u/np
    s = ssξ*ξ

    function Eq(F)
        @cal_common
        return (F/wp_bar)-rho_F
    end

    F = fzero(Eq,1.2,2.)

    @cal_common

    cal = Parameters()
    Macros.@save_struct cal η r s λ α m b γ σ F δ
    return cal

end
end


text:
# linear stochastic differential equation
using LinearAlgebra
using ZigZagBoomerang
const ZZB = ZigZagBoomerang

n = 10
X  = -1.0*I(n) + zeros(n,n)
for i in 1:n
    j = rand(collect(1:n)[1:end .!= i]) # pick a random off diagonal element
    X[i,j] = 1.0
end
X = 2.0*X
X
real.(eigvals(X))

# Simulate data y
function simulate_data(X)
    n = size(X)[1]
    σ = 0.1
    b(y) = X*y
    sigma(x) = fill(σ, n)
    δt = 0.001
    T = 10.0
    tt = range(0.0, step = δt, stop = T)
    Y = zeros(n, length(tt))
    y0 = randn(n)
    Y[:, 1] = y0
# euler forward simulation
    for i in eachindex(tt)[2:end]
        dt = tt[i] - tt[i-1]
        y = Y[:, i-1]
        Y[:, i] = y + dt*b(y) + sigma(y).*randn(n)
    end
    Y, tt
end
Y, tt = simulate_data(X)
using Plots
plot(tt, Y[1, :])

## helper for ∫y^(i)_t dy^(j)_t for each i,j
# with ito approximation: see https://en.wikipedia.org/wiki/It%C3%B4_calculus
function tilde_Psi(XX, i, j, tt)
    res = 0.0
    for k in eachindex(tt)[2:end]
        yi = XX[i, k-1]
        dyj = (XX[j, k] - XX[j, k-1])
        res += yi*dyj
    end
    res
end

## helper for ∫y^(i)_t y^(j)_t dt for each i,j
function Psi(XX, i, j, tt)
    res = 0.0
    for k in eachindex(tt)[2:end]
        yi = (XX[i, k-1] + XX[i, k])/2
        yj = (XX[j, k-1] + XX[j, k])/2
        dt = tt[k] - tt[k-1]
        res += yi*yj*dt
    end
    res
end

#pre-computing matrices
tilde_Ψ = [tilde_Psi(Y, i[1], i[2], tt) for i in CartesianIndices(X)]
tilde_Ψt = tilde_Ψ'

Ψ = [Psi(Y, i[1], i[2], tt) for i in CartesianIndices(X)]
Ψ - Ψ' ≈ zeros(n,n) ? println("symmetric matrix") : error("not symmetric matrix")

# x sparse, y whatever
function sdot(x, y)
    r = 0.0
    @inbounds for i in 1:length(x)
        if x[i] == 0
            continue
        else
            r += x[i]*y[i]
        end
    end
    r
end

# matrix to vector

# vector to matrix

# stupid: getting column of matrix from vector indexation of matrix
col(k,n) = div(k-1, n) + 1
row(k,n) = rem(k-1, n)+1

firstcol(k, n) = div(k-1,n)*n + 1
firstrow(k, n) = rem(k-1, n)*n + 1

# partial ϕi (Gaussian)
function ∇ϕ(x, k, tilde_Ψ, tilde_Ψt, Ψ, γ0, n)
    i, j = row(k,n), col(k,n)
    x1 = reshape(x, n, n)
    return -tilde_Ψt[k] + sdot(x1[i,:], Ψ[:,j]) + γ0*x[k]
end

function ∇ϕ(x, k, tilde_Ψ, tilde_Ψt, Ψ, γ0, n)
    i, j = row(k,n), col(k,n)
    x1 = reshape(x, n, n)
    return -tilde_Ψt[k] + sdot(x1[i,:], Ψ[:,j]) + γ0*x[k]
    # return -tilde_Ψt[k] +  γ0*x[k]
end

## Override ab
function ZZB.ab(G, k, x, θ, c, F::ZigZag,  tilde_Ψ, tilde_Ψt, Ψ, γ0, n)
    i, j = row(k,n), col(k,n)
    x1 = reshape(x, n,n)
    θ1 = reshape(θ, n,n)
    a = c + (θ[k]*(-tilde_Ψt[k] + sdot(x1[i,:], Ψ[:,j]) + γ0*x[k]))
    # a = c + (θ[k]*(-tilde_Ψt[k] + γ0*x[k]))
    b =  θ[k]*(sdot(θ1[i,:], Ψ[:,j]) + γ0*θ[k])
    # return c + θ[k]*γ0*(x[k] - 1.0), θ[k]^2*γ0
    return a, b
end

# Prior σ0^(-1)
γ0 = 0.2
#kappa given the prior (here Gaussian) and w
nz = sum(X[:] .!= 0) #non zero elements
w = nz/(n*n) # proportion of non zero elemenets
κ = (γ0/sqrt(2π))/(1/w -1)
p = n*n
μ = zeros(p)

# TODO
# Γ = ones(p,p)
Γ = ones(p,p)

Z = ZigZag(Γ, μ)
c = 0.00001
t0, x0, T = 0.0, randn(p), 1000.0
θ0 = rand([-1.0,1.0], p)
su = false #strong upperbounds
adapt = false
traj, (t, x, θ), (acc, num), c = @time ZZB.sspdmp(∇ϕ, t0, x0, θ0, T, c, Z, κ,
                                                 tilde_Ψ, tilde_Ψt, Ψ, γ0, n;
                                                strong_upperbounds = su ,
                                                adapt = adapt)


text:
using LinearAlgebra
using ContinuousDPs
using QuantEcon: PFI, VFI, solve, stationary_values
using BasisMatrices: Basis, ChebParams, SplineParams, LinParams, nodes
using Test

include("test_cdp.jl")
include("test_lq_approx.jl")


text:
using DataFrames
using CSV
using GLM
using Lathe

#df_gas = DataFrame(CSV.File("../resources/data1.csv"))
df_gas = DataFrame(CSV.File("../resources/data1_v1.csv"))

colnames = Symbol[]
for i in string.(names(df_gas))
    i = replace(strip(i), " " => "_")
    i = replace(i, "-" => "")
    i = replace(i, "/" => "")
    i = replace(i, "(" => "")
    i = replace(i, ")" => "")
    i = replace(i, "%" => "")
    i = replace(i, "," => "")
    i = replace(i, "*" => "")
    println(i)
    push!(colnames, Symbol(i))
    #push!(colnames, Symbol(replace(replace(replace(strip(i), " " => "_"), "-" => "_"), "/" => "_")))
end

rename!(df_gas, colnames)
println(names(df_gas))
formulasGasT = FormulaTerm[]

fm0 = @formula(Electrical_Power_MW ~ Percent_of_Gas_Turbine_Electrical_Load)
push!(formulasGasT, fm0)
#fm1 = @formula(Natural_Gas_Flow_lbhr ~ Percent_of_Gas_Turbine_Electrical_Load)
fm1 = @formula(HHV_Heat_Input_MMBTUhr ~ Percent_of_Gas_Turbine_Electrical_Load)
push!(formulasGasT, fm1)
fm2 = @formula(CO2_Emissions_lbhr ~ Percent_of_Gas_Turbine_Electrical_Load)
push!(formulasGasT, fm2)
fm3 = @formula(HP_Turbine_Power_kW ~ Percent_of_Gas_Turbine_Electrical_Load)
push!(formulasGasT, fm3)
fm4 = @formula(IP_Turbine_Power_kW ~ Percent_of_Gas_Turbine_Electrical_Load)
push!(formulasGasT, fm4)
fm5 = @formula(Total_kW ~ Percent_of_Gas_Turbine_Electrical_Load) # auxiliary
push!(formulasGasT, fm5)

coefReg = Vector{Float64}[]
coef_gas_df = DataFrame(y = String[], b = Float64[], a = Float64[])
for i in formulasGasT
    linReg = lm(i, df_gas)
    y = string(terms(i)[1])
    b = coef(linReg)[1]
    a = coef(linReg)[2]
    push!(coef_gas_df, (y, b, a))
    println(linReg)
    println(deviance(linReg))
    println(coef(linReg))
end

CSV.write("gas_coeffs.csv", coef_gas_df)

df_steam_t = DataFrame(CSV.File("../resources/data2.csv"))
coln_steam = Symbol[]
for i in string.(names(df_steam_t))
    i = replace(strip(i), " " => "_")
    i = replace(i, "-" => "")
    i = replace(i, "/" => "")
    i = replace(i, "(" => "")
    i = replace(i, ")" => "")
    i = replace(i, "%" => "")
    i = replace(i, "," => "")
    i = replace(i, "*" => "")
    i = replace(i, "+" => "")
    println(i)
    push!(coln_steam, Symbol(i))
end

rename!(df_steam_t, coln_steam)
#LP_Turbine_Power_kW
#Condenser_Duty_MMBTUhr
#DAC_Steam_Duty_MMBTUhr
#Total_ST_and_BOP_Auxiliary_Loads_kW
formulasSteamT = FormulaTerm[]
#
fm0st = @formula(LP_Turbine_Power_kW ~ Percent_of_Gas_Turbine_Electrical_Load)
push!(formulasSteamT, fm0st)
fm1st = @formula(Duty_of_CC_Reboiler_MMBTUhr ~ Percent_of_Gas_Turbine_Electrical_Load)
push!(formulasSteamT, fm1st)
fm2st = @formula(DAC_Steam_Duty_MMBTUhr ~ Percent_of_Gas_Turbine_Electrical_Load)
push!(formulasSteamT, fm2st)
fm3st = @formula(Total_ST_and_BOP_Auxiliary_Loads_kW ~ Percent_of_Gas_Turbine_Electrical_Load)
push!(formulasSteamT, fm3st)
coef_steam_df = DataFrame(y = String[], b = Float64[], a = Float64[])
for i in formulasSteamT
    linReg = lm(i, df_steam_t)
    y = string(terms(i)[1])
    b = coef(linReg)[1]
    a = coef(linReg)[2]
    push!(coef_steam_df, (y, b, a))
    println(linReg)
    println(deviance(linReg))
    println(coef(linReg))
end

CSV.write("steam_coeffs.csv", coef_steam_df)

max_stars_repo_name: dthierry/princetonDacLti

max_issues_repo_name: dthierry/princetonDacLti

max_forks_repo_name: dthierry/princetonDacLti